/** \brief 2D/3D point class. -*- c++ -*- */
/***************************************************************************
    begin                : 2006-02-02
    copyright            : (C) 2006 by Gunter Winkler
    email                : guwi17@gmx.de
***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef _FTL_POINT_H_
#define _FTL_POINT_H_

#include <boost/numeric/ublas/storage.hpp>
//#include <boost/numeric/ublas/storage_adaptors.hpp>
#include "storage_adaptors.hpp"
#include <iosfwd>

namespace ftl 
{

  template <int dim, class T = double> 
  class Point 
  {
  public:
    //types
    typedef Point<dim, T>   self_type;

    typedef Point& closure_type;
    typedef const Point& const_closure_type;
    typedef T *pointer;
    typedef const T *const_pointer;
    typedef T &reference;
    typedef const T &const_reference;
    typedef T value_type;

    /* iterators are simple pointers */
    typedef const_pointer const_iterator;
    typedef pointer       iterator;

    typedef size_t   size_type;
    typedef int      difference_type;

    enum { dimension = dim };

    // reduced point type (without first/last component)
    typedef Point<dim-1, T> reduced_point_type;

    // data
    T coords[dim];

    // constructors
    Point() {}

    // general initializing constructor
    explicit Point(const T& x) 
    {
      for (size_t i=0; i<dim; ++i) coords[i] = x;
    }
    // abbreviation for dim == 2 && dim == 3
    Point(const T& x, const T&y) 
    {
      assert( dim == 2 );
      coords[0] = x;
      coords[1] = y;
    }
    Point(const T& x, const T&y, const T& z) 
    {
      assert( dim == 3 );
      coords[0] = x;
      coords[1] = y;
      coords[2] = z;
    }
    // creation from data pointer - warning: no range checks!
    Point(const T * px) 
    {
      for (size_t i=0; i<dim; ++i) coords[i] = px[i];
    }
    
    
    // assignement
    Point & operator= (const Point& x)
    {
      for (size_t i=0; i<dim; ++i) coords[i] = x[i];
      return (*this);
    }
    Point & operator= (const T& x)
    {
      for (size_t i=0; i<dim; ++i) coords[i] = x;
      return (*this);
    }

    // access
    const T & operator[](const size_t index) const
    {
      assert( (0 <= index) && (index < dim) );
      return coords[index];
    }
    T & operator[](const size_t index) 
    {
      assert( (0 <= index) && (index < dim) );
      return coords[index];
    }
    const T & operator()(const size_t index) const
    {
      assert( (0 <= index) && (index < dim) );
      return coords[index];
    }
    T & operator()(const size_t index) 
    {
      assert( (0 <= index) && (index < dim) );
      return coords[index];
    }
    // size
    std::size_t size() const { return std::size_t( dim ); }
    // scalar operation
    Point & operator += (const T& x)
    {
      for (size_t i=0; i<dim; ++i) coords[i] += x;
      return (*this);
    }
    Point & operator -= (const T& x)
    {
      for (size_t i=0; i<dim; ++i) coords[i] -= x;
      return (*this);
    }
    Point & operator *= (const T& x)
    {
      for (size_t i=0; i<dim; ++i) coords[i] *= x;
      return (*this);
    }
    Point & operator /= (const T& x)
    {
      for (size_t i=0; i<dim; ++i) coords[i] /= x;
      return (*this);
    }
    // vector operations
    Point & operator += (const Point& x)
    {
      for (size_t i=0; i<dim; ++i) coords[i] += x[i];
      return (*this);
    }
    Point & operator -= (const Point& x)
    {
      for (size_t i=0; i<dim; ++i) coords[i] -= x[i];
      return (*this);
    }
    // unary minus
    void negate () 
    {
      for (size_t i=0; i<dim; ++i) coords[i] = -coords[i];
      return;
    }
    // inner product
    T inner_prod(const Point& p) const
    {
      T temp(0);
      for (size_t i=0; i<dim; ++i) temp += coords[i] * p[i];
      return temp;      
    }
    // iterators
    iterator begin()
    {
      return &(coords[0]);
    }
    iterator end()
    {
      return &(coords[dim]);
    }
    // const iterators
    const_iterator begin() const
    {
      return &(coords[0]);
    }
    const_iterator end() const
    {
      return &(coords[dim]);
    }
    // raw access
    pointer data() 
    {
      return &(coords[0]);
    }
    const_pointer data() const
    {
      return &(coords[0]);
    }
    // reduce size
    reduced_point_type without_last() const
    {
      return reduced_point_type( data() );
    }
    reduced_point_type without_first() const
    {
      return reduced_point_type( ++ data() );
    }
    template<int k>
    Point<k, value_type> first_k() const
    {
      return Point<k, value_type>( data() );
    }
    template<int k>
    Point<k, value_type> last_k() const
    {
      return Point<k, value_type>( data() + (dim-k) );
    }
    // i/o
    std::istream & read(std::istream & is)
    {
      for (size_t i=0; i<dim; ++i) {
        is >> coords[i];
      }
      return is;
    }
    std::ostream & write(std::ostream & os) const
    {
      for (size_t i=0; i<dim; ++i) {
        os << coords[i] << " ";
      }
      return os;
    }
  };

  // scalar arithmetic
  template <int dim, class T>
  Point<dim, T> operator + (const T& x, const Point<dim, T>& p)
  {
    Point<dim, T> temp = p;
    return (temp += x);
  }
  template <int dim, class T>
  Point<dim, T> operator - (const T& x, const Point<dim, T>& p)
  {
    Point<dim, T> temp = p;
    temp.negate();
    return (temp += x);
  }
  template <int dim, class T>
  Point<dim, T> operator + (const Point<dim, T>& p, const T& x)
  {
    Point<dim, T> temp = p;
    return (temp += x);
  }
  template <int dim, class T>
  Point<dim, T> operator - (const Point<dim, T>& p, const T& x)
  {
    Point<dim, T> temp = p;
    return (temp -= x);
  }
  template <int dim, class T>
  Point<dim, T> operator * (const T& x, const Point<dim, T>& p)
  {
    Point<dim, T> temp = p;
    return (temp *= x);
  }
  template <int dim, class T>
  Point<dim, T> operator * (const Point<dim, T>& p, const T& x)
  {
    Point<dim, T> temp = p;
    return (temp *= x);
  }
  // vector arithmetic
  template <int dim, class T>
  Point<dim, T> operator + (const Point<dim, T>& p, const Point<dim, T>& q)
  {
    Point<dim, T> temp = p;
    return (temp += q);
  }
  template <int dim, class T>
  Point<dim, T> operator - (const Point<dim, T>& p, const Point<dim, T>& q)
  {
    Point<dim, T> temp = p;
    return (temp -= q);
  }
  // inner product
  template <int dim, class T>
  T operator * (const Point<dim, T>& x, const Point<dim, T>& p)
  {
    return p.inner_prod(x);
  }
  // i/o
  template <int dim, class T>
  std::ostream & operator<< (std::ostream & os, const Point<dim, T>& p)
  {
    return p.write(os);
  }
  template <int dim, class T>
  std::istream & operator>> (std::istream & os, Point<dim, T>& p)
  {
    return p.read(os);
  }
  // norms
  template <int dim, class T>
  T norm_2 (const Point<dim, T>& p)
  {
    return sqrt( p.inner_prod(p) );
  }
  template <int dim, class T>
  T norm_inf (const Point<dim, T>& p)
  {
    T temp = std::abs(p[0]);
    for (size_t i=1; i<dim; ++i) {
      if (temp < std::abs(p[i])) temp = std::abs(p[i]);
    }
    return temp;
  }

  // some ublas compatibility routines
  template <int dim, class T>
  Point<dim, T>& noalias (Point<dim, T>& p)
  {
    return p;
  }
  // only quadratic matrices of size dim-by-dim supported
  template <int dim, class T, class E>
  Point<dim, T> prod( const boost::numeric::ublas::matrix_expression<E> & e, const Point<dim, T>& p) 
  {
    Point<dim, T> temp;
    assert( e().size2() == dim && e().size1() == dim);
    for (size_t i=0; i<dim; ++i) {
      T t = 0.0;
      for (size_t j=0; j<dim; ++j) {
        t += e()(i,j) * p[j];
      }
      temp[i] = t;
    }
    return temp;
  }
  template <int dim, class T, class E>
  Point<dim, T> prod( const Point<dim, T>& p, const boost::numeric::ublas::matrix_expression<E> & e ) 
  {
    Point<dim, T> temp;
    assert( e().size2() == dim && e().size1() == dim);
    for (size_t i=0; i<dim; ++i) {
      T t = 0.0;
      for (size_t j=0; j<dim; ++j) {
        t += e()(j,i) * p[j];
      }
      temp[i] = t;
    }
    return temp;
  }
  template <int dim, class T>
  T inner_prod(const Point<dim, T>& p1, const Point<dim, T>& p2)
  {
    return p1.inner_prod(p2);
  }

  template <int dim, class T>
  boost::numeric::ublas::vector<const T, boost::numeric::ublas::readonly_array_adaptor<T> >
  make_vector(const Point<dim, T>& p)
  {
    return boost::numeric::ublas::make_vector_from_pointer(dim, p.data());
  }
}


#endif