//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the 
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
// 
#ifndef VECTOR_HPP
#define VECTOR_HPP
/**
 * \file Vector.hpp
 *
 * Defines the Vector class template
 * This file is shared by cuda, do not include headers that nvcc can't handle (i.e. Qt)
 */

#include <Config.hpp>

#include <cuda/CudaGlobal.hpp>
#include <StaticAssert.hpp>

#include <cassert>
#include <cmath>
#include <cstdarg>
#include <iostream>
#include <math.h>

#ifndef __CUDACC__
  #include <QTextStream>
#endif

namespace mgx 
{
  /**
   * \brief Namespace containing all the utility classes.
   */
  /**
   * \class Vector Vector.hpp <vector.h>
   *
   * Vector class supporting all classic classic vector operations.
   *
   * \note The addition and subtraction with a scalar is possible and is
   * equivalent to adding/subtracting a vector filled with the scalar.
   */
  template <size_t dim, class T = float> class Vector {
  protected:
    T elems[dim];
  
  public:
    typedef T value_type;
    typedef T& reference_type;
    typedef const T& const_reference_type;
    typedef T* pointer_type;
    typedef const T* const_pointer_type;
    typedef T* iterator;
    typedef const T* const_iterator;

    static const size_t numElems = dim; 
// This defeats tr1::has_trivial_copy macro, and it's trivial anyway
//    /**
//     * Copy another vector
//     */
//    CU_HOST_DEVICE
//    Vector(const Vector& vec)
//    {
//      for(size_t i = 0; i < dim; i++)
//        elems[i] = vec[i];
//    }
//  
//    /**
//     * Copy another vector with different number of elements
//     */
//    template <class T1> CU_HOST_DEVICE Vector(const Vector<dim, T1>& vec)
//    {
//      for(size_t i = 0; i < dim; ++i)
//        elems[i] = vec[i];
//    }
  
    /**
     * Copy another vector with different number of elements
     */
    template <size_t d1, class T1> CU_HOST_DEVICE explicit Vector(const Vector<d1, T1>& vec)
    {
      if(d1 < dim) {
        for(size_t i = 0; i < d1; ++i)
          elems[i] = vec[i];
        for(size_t i = d1; i < dim; ++i)
          elems[i] = 0;
      } else {
        for(size_t i = 0; i < dim; ++i)
          elems[i] = vec[i];
      }
    }
  
  #ifndef COMPILE_CUDA
    /**
     * Initialize a vector from any object behaving like an array.
     *
     * The only constraints are:
     *  - the type has to be convertible to \c T
     *  - the vector needs a [] operator
     *  - the size of the vector has to be at least \c dim
     */
    template <class Vec> CU_HOST_DEVICE 
    explicit Vector(const Vec& el)
    {
      for(size_t i = 0; i < dim; i++)
        elems[i] = T(el[i]);
    }
  #endif
  
    /**
     * Initialize a vector with all values to \c x.
     */
    CU_HOST_DEVICE
    explicit Vector(const T& x = T())
    {
      for(size_t i = 0; i < dim; ++i)
        elems[i] = x;
    }
  
    /**
     * Initialize a 2D vector
     */
    CU_HOST_DEVICE
    explicit Vector(const T& x, const T& y)
    {
      STATIC_ASSERT(dim == 2);
      elems[0] = x;
      elems[1] = y;
    }
  
    /**
     * Initialize a 3D vector
     */
    CU_HOST_DEVICE
    explicit Vector(const T& x, const T& y, const T& z)
    {
      STATIC_ASSERT(dim == 3);
      elems[0] = x;
      elems[1] = y;
      elems[2] = z;
    }
  
    /**
     * Initialize a 4D vector
     */
    CU_HOST_DEVICE
    explicit Vector(const T& x, const T& y, const T& z, const T& t)
    {
      STATIC_ASSERT(dim == 4);
      elems[0] = x;
      elems[1] = y;
      elems[2] = z;
      elems[3] = t;
    }
  
    /**
     * Initialize a 5D vector
     */
    CU_HOST_DEVICE
    Vector(const T& x, const T& y, const T& z, const T& a, const T& b)
    {
      STATIC_ASSERT(dim == 5);
      elems[0] = x;
      elems[1] = y;
      elems[2] = z;
      elems[3] = a;
      elems[4] = b;
    }
  
    /**
     * Initialize a 6D vector
     */
    CU_HOST_DEVICE
    Vector(const T& x, const T& y, const T& z, const T& a, const T& b, const T& c)
    {
      STATIC_ASSERT(dim == 6);
      elems[0] = x;
      elems[1] = y;
      elems[2] = z;
      elems[3] = a;
      elems[4] = b;
      elems[5] = c;
    }
  
    /**
     * Initialize a 7D vector
     */
    CU_HOST_DEVICE
    Vector(const T& x, const T& y, const T& z, const T& a, const T& b, const T& c, const T& d)
    {
      STATIC_ASSERT(dim == 7);
      elems[0] = x;
      elems[1] = y;
      elems[2] = z;
      elems[3] = a;
      elems[4] = b;
      elems[5] = c;
      elems[6] = d;
    }

    /**
     * Initialize a 8D vector
     */
    CU_HOST_DEVICE
    Vector(const T& a, const T& b, const T& c, const T& d, const T& e, const T& f, const T& g, const T& h)
    {
      STATIC_ASSERT(dim == 8);
      elems[0] = a;
      elems[1] = b;
      elems[2] = c;
      elems[3] = d;
      elems[4] = e;
      elems[5] = f;
      elems[6] = g;
      elems[7] = h;
    } 

    /**
     * Initialize a 9D vector
     */
    CU_HOST_DEVICE
    Vector(const T& a, const T& b, const T& c, const T& d, const T& e, const T& f, 
                                                             const T& g, const T& h, const T& i)
    {
      STATIC_ASSERT(dim == 9);
      elems[0] = a;
      elems[1] = b;
      elems[2] = c;
      elems[3] = d;
      elems[4] = e;
      elems[5] = f;
      elems[6] = g;
      elems[7] = h;
      elems[8] = i;
    }

    /**
     * Initialize a 12D vector
     */
    CU_HOST_DEVICE
    Vector(const T& a, const T& b, const T& c, const T& d, const T& e, const T& f, const T& g, 
										             const T& h, const T& i, const T& j, const T& k, const T& l)
    {
      STATIC_ASSERT(dim == 12);
      elems[0] = a;
      elems[1] = b;
      elems[2] = c;
      elems[3] = d;
      elems[4] = e;
      elems[5] = f;
      elems[6] = g;
      elems[7] = h;
      elems[8] = i;
      elems[9] = j;
      elems[10] = k;
      elems[11] = l;
    }

    /**
     * Initialize a 18D vector
     */
    CU_HOST_DEVICE
    Vector(const T& a, const T& b, const T& c, const T& d, const T& e, const T& f, 
           const T& g, const T& h, const T& i, const T& j, const T& k, const T& l, 
           const T& m, const T& n, const T& o, const T& p, const T& q, const T& r)
    {
      STATIC_ASSERT(dim == 18);
      elems[0] = a;
      elems[1] = b;
      elems[2] = c;
      elems[3] = d;
      elems[4] = e;
      elems[5] = f;
      elems[6] = g;
      elems[7] = h;
      elems[8] = i;
      elems[9] = j;
      elems[10] = k;
      elems[11] = l;
      elems[12] = m;
      elems[13] = n;
      elems[14] = o;
      elems[15] = p;
      elems[16] = q;
      elems[17] = r;
    } 

    /**
     * Returns the size of the vector (i.e. the number of elements)
     */
    CU_HOST_DEVICE
    static size_t size() { return dim; }
  
    /**
     * Returns a raw pointer on the data
     */
    CU_HOST_DEVICE
    T* data() { return elems; }
  
    /**
     * STL-iteration begin
     */
    CU_HOST_DEVICE
    iterator begin() { return elems; }
    /**
     * Stl-iteration constant begin
     */
    CU_HOST_DEVICE
    const_iterator begin() const { return elems; }
  
    /**
     * STL-iteration end
     */
    CU_HOST_DEVICE
    iterator end() { return elems + dim; }
    /**
     * Stl-iteration constant end
     */
    CU_HOST_DEVICE
    const_iterator end() const { return elems + dim; }
  
    /**
     * Returns a constant raw pointer on the data
     */
    CU_HOST_DEVICE
    const T* c_data() const { return elems; }
  
    /**
     * Vector negation
     */
    CU_HOST_DEVICE
    Vector operator-(void) const
    {
      Vector ans;
      for(size_t i = 0; i < dim; i++)
        ans[i] = -elems[i];
  
      return ans;
    }
  
    /**
     * Vector addition
     */
    CU_HOST_DEVICE
    Vector operator+(const Vector& vec) const
    {
      Vector ans(*this);
      ans += vec;
      return ans;
    }
  
    /**
     * Vector subtraction
     */
    CU_HOST_DEVICE
    Vector operator-(const Vector& vec) const
    {
      Vector ans(*this);
      ans -= vec;
      return ans;
    }
  
    /**
     * Element-wise multiplcation
     */
    CU_HOST_DEVICE
    Vector mult(const Vector& vec) const
    {
      Vector ans;
      for(size_t i = 0; i < dim; i++)
        ans[i] = elems[i] * vec.elems[i];
      return ans;
    }
  
    /**
     * Multiplication by a scalar
     */
    CU_HOST_DEVICE
    Vector operator*(const T& scalar) const
    {
      Vector ans(*this);
      ans *= scalar;
      return ans;
    }
  
    /**
     * Division by a scalar
     */
    CU_HOST_DEVICE
    Vector operator/(const T& scalar) const
    {
      Vector ans(*this);
      ans /= scalar;
      return ans;
    }
  
    /**
     * Element-wise division
     */
    CU_HOST_DEVICE
    Vector operator/(const Vector& vec) const
    {
      Vector ans = *this;
      ans /= vec;
      return ans;
    }
  
    /**
     * In-place element-wise division by a scalar
     */
    CU_HOST_DEVICE
    Vector& operator/=(const Vector& vec)
    {
      for(size_t i = 0; i < dim; ++i)
        elems[i] /= vec.elems[i];
      return *this;
    }
  
    /**
     * Multiplication by a scalar
     */
    CU_HOST_DEVICE
    friend Vector operator*(const T& scalar, const Vector& vec)
    {
      Vector ans;
      for(size_t i = 0; i < dim; i++)
        ans[i] = scalar * vec.elems[i];
  
      return ans;
    }
  
    /**
     * Dot product
     */
    CU_HOST_DEVICE
    T operator*(const Vector& vec) const
    {
      T ans = 0;
      for(size_t i = 0; i < dim; i++)
        ans += elems[i] * vec.elems[i];
  
      return ans;
    }
  
    /**
     * Vector copy
     */
//    CU_HOST_DEVICE
//    Vector& operator=(const Vector& vec)
//    {
//      for(size_t i = 0; i < dim; i++)
//        elems[i] = vec.elems[i];
//  
//      return (*this);
//    }
  
    /**
     * In-place vector addition
     */
    CU_HOST_DEVICE
    Vector& operator+=(const Vector& vec)
    {
      for(size_t i = 0; i < dim; i++)
        elems[i] += vec.elems[i];
      return *this;
    }
  
    /**
     * In-place constant addition
     */
    CU_HOST_DEVICE
    Vector& operator+=(const T& val)
    {
      for(size_t i = 0; i < dim; i++)
        elems[i] += val;
      return *this;
    }
  
    /**
     * In-place vector subtraction
     */
    CU_HOST_DEVICE
    Vector& operator-=(const Vector& vec)
    {
      for(size_t i = 0; i < dim; i++)
        elems[i] -= vec.elems[i];
      return *this;
    }
  
    /**
     * In-place value subtraction
     */
    CU_HOST_DEVICE
    Vector& operator-=(const T& val)
    {
      for(size_t i = 0; i < dim; i++)
        elems[i] -= val;
      return *this;
    }
  
    /**
     * In-place multiplication by a scalar
     */
    CU_HOST_DEVICE
    Vector& operator*=(const T& scalar)
    {
      for(size_t i = 0; i < dim; i++)
        elems[i] *= scalar;
      return *this;
    }
  
    /**
     * In-place multiplication by a scalar
     */
    template <typename T1> CU_HOST_DEVICE 
		Vector& operator*=(const T1& scalar)
    {
      for(size_t i = 0; i < dim; i++)
        elems[i] = (T)(elems[i] * scalar);
      return *this;
    }
  
    /**
     * In-place division by a scalar
     */
    CU_HOST_DEVICE
    Vector& operator/=(const T& scalar)
    {
      for(size_t i = 0; i < dim; ++i)
        elems[i] /= scalar;
      return *this;
    }
  
    /**
     * In-place division by a scalar
     */
    template <typename T1> CU_HOST_DEVICE 
	  Vector& operator/=(const T1& scalar)
    {
      for(size_t i = 0; i < dim; ++i)
        elems[i] = (T)(elems[i] / scalar);
      return *this;
    }
  
    /**
     * Element-wise equality
     */
    CU_HOST_DEVICE
    bool operator==(const Vector& vec) const
    {
      for(size_t i = 0; i < dim; i++)
        if(elems[i] != vec.elems[i])
          return false;
  
      return true;
    }
  
    /**
     * Element-wise inequality
     */
    CU_HOST_DEVICE
    bool operator!=(const Vector& vec) const {
      return (!((*this) == vec));
    }
  
    /**
     * Access to the element \c idx
     */
    CU_HOST_DEVICE
    T& operator[](size_t idx) {
      return elems[idx];
    }
  
    /**
     * Access to the element \c idx
     */
    CU_HOST_DEVICE
    const T& operator[](size_t idx) const {
      return elems[idx];
    }
  
    /**
     * Euclidean norm of the vector
     */
    CU_HOST_DEVICE
    T norm() const
    {
      #ifdef _MSC_VER
			  // Incorrect, but stupid compiler can't compile correct code
        return std::sqrt((float)normsq()); 
      #else
        return std::sqrt(normsq());
      #endif
    }
  
    /**
     * Square of the Euclidean norm of the vector
     */
    CU_HOST_DEVICE
    T normsq() const
    {
      T ans = 0;
      for(size_t i = 0; i < dim; i++)
        ans += elems[i] * elems[i];
  
      return ans;
    }
  
    /**
     * Normalize the vector
     */
    CU_HOST_DEVICE
    Vector& normalize(void)
    {
      T sz = norm();
      return ((*this) /= sz);
    }
  
    /**
     * Returns a normalized version of the vector
     */
    CU_HOST_DEVICE
    Vector normalized(void) const
    {
      Vector ans(*this);
      return ans.normalize();
    }
  
    bool iszero(void)
    {
      for(size_t i = 0; i < dim; i++)
        if(elems[i] != 0)
          return false;
      return true;
    }
  
    Vector& zero(void)
    {
      for(size_t i = 0; i < dim; i++)
        elems[i] = 0;
      return (*this);
    }
  
    /**
     * Set the values of a 1-D vector
     */
    CU_HOST_DEVICE
    void set(const T& x)
    {
      STATIC_ASSERT(dim == 1);
      elems[0] = x;
    }
  
    /**
     * Set the values of a 2-D vector
     */
    CU_HOST_DEVICE
    void set(const T& x, const T& y)
    {
      STATIC_ASSERT(dim == 2);
      elems[0] = x;
      elems[1] = y;
    }
  
    /**
     * Set the values of a 3-D vector
     */
    CU_HOST_DEVICE
    void set(const T& x, const T& y, const T& z)
    {
      STATIC_ASSERT(dim == 3);
      elems[0] = x;
      elems[1] = y;
      elems[2] = z;
    }
  
    /**
     * Set the values of a 4-D vector
     */
    CU_HOST_DEVICE
    void set(const T& x, const T& y, const T& z, const T& t)
    {
      STATIC_ASSERT(dim == 4);
      elems[0] = x;
      elems[1] = y;
      elems[2] = z;
      elems[3] = t;
    }
  
    /**
     * Set all the elements to \c value
     */
    CU_HOST_DEVICE
    Vector& operator=(const T& value)
    {
      for(size_t i = 0; i < dim; ++i) {
        elems[i] = value;
      }
      return *this;
    }
  
    /**
     * Compute the cross product as \c this x \c other
     */
    CU_HOST_DEVICE
    Vector cross(const Vector& other) const
    {
      STATIC_ASSERT(dim == 3);
      return (*this) ^ other;
    }
  
    /**
     * Short access to the first element
     */
    CU_HOST_DEVICE
    void x(const T& v)
    {
      STATIC_ASSERT(dim > 0);
      elems[0] = v;
    }
    /**
     * Short access to the second element
     */
    CU_HOST_DEVICE
    void y(const T& v)
    {
      STATIC_ASSERT(dim > 1);
      elems[1] = v;
    }
    /**
     * Short access to the third element
     */
    CU_HOST_DEVICE
    void z(const T& v)
    {
      STATIC_ASSERT(dim > 2);
      elems[2] = v;
    }
    /**
     * Short access to the fourth element
     */
    CU_HOST_DEVICE
    void t(const T& v)
    {
      STATIC_ASSERT(dim > 3);
      elems[3] = v;
    }
    /**
     * Short access to the first element
     */
    CU_HOST_DEVICE
    T& x()
    {
      STATIC_ASSERT(dim > 0);
      return elems[0];
    }
    /**
     * Short access to the second element
     */
    CU_HOST_DEVICE
    T& y()
    {
      STATIC_ASSERT(dim > 1);
      return elems[1];
    }
    /**
     * Short access to the third element
     */
    CU_HOST_DEVICE
    T& z()
    {
      STATIC_ASSERT(dim > 2);
      return elems[2];
    }
    /**
     * Short access to the fourth element
     */
    CU_HOST_DEVICE
    T& t()
    {
      STATIC_ASSERT(dim > 3);
      return elems[3];
    }
    /**
     * Short access to the first element
     */
    CU_HOST_DEVICE
    const T& x() const
    {
      STATIC_ASSERT(dim > 0);
      return elems[0];
    }
    /**
     * Short access to the second element
     */
    CU_HOST_DEVICE
    const T& y() const
    {
      STATIC_ASSERT(dim > 1);
      return elems[1];
    }
    /**
     * Short access to the third element
     */
    CU_HOST_DEVICE
    const T& z() const
    {
      STATIC_ASSERT(dim > 2);
      return elems[2];
    }
    /**
     * Short access to the fourth element
     */
    CU_HOST_DEVICE
    const T& t() const
    {
      STATIC_ASSERT(dim > 3);
      return elems[3];
    }
  
    /**
     * Short access to the first element
     */
    CU_HOST_DEVICE
    void i(const T& v)
    {
      STATIC_ASSERT(dim > 0);
      elems[0] = v;
    }
    /**
     * Short access to the second element
     */
    CU_HOST_DEVICE
    void j(const T& v)
    {
      STATIC_ASSERT(dim > 1);
      elems[1] = v;
    }
    /**
     * Short access to the third element
     */
    CU_HOST_DEVICE
    void k(const T& v)
    {
      STATIC_ASSERT(dim > 2);
      elems[2] = v;
    }
    /**
     * Short access to the fourth element
     */
    CU_HOST_DEVICE
    void l(const T& v)
    {
      STATIC_ASSERT(dim > 3);
      elems[3] = v;
    }
    /**
     * Short access to the first element
     */
    CU_HOST_DEVICE
    T& i()
    {
      STATIC_ASSERT(dim > 0);
      return elems[0];
    }
    /**
     * Short access to the second element
     */
    CU_HOST_DEVICE
    T& j()
    {
      STATIC_ASSERT(dim > 1);
      return elems[1];
    }
    /**
     * Short access to the third element
     */
    CU_HOST_DEVICE
    T& k()
    {
      STATIC_ASSERT(dim > 2);
      return elems[2];
    }
    /**
     * Short access to the fourth element
     */
    CU_HOST_DEVICE
    T& l()
    {
      STATIC_ASSERT(dim > 3);
      return elems[3];
    }
    /**
     * Short access to the first element
     */
    CU_HOST_DEVICE
    const T& i() const
    {
      STATIC_ASSERT(dim > 0);
      return elems[0];
    }
    /**
     * Short access to the second element
     */
    CU_HOST_DEVICE
    const T& j() const
    {
      STATIC_ASSERT(dim > 1);
      return elems[1];
    }
    /**
     * Short access to the third element
     */
    CU_HOST_DEVICE
    const T& k() const
    {
      STATIC_ASSERT(dim > 2);
      return elems[2];
    }
    /**
     * Short access to the fourth element
     */
    CU_HOST_DEVICE
    const T& l() const
    {
      STATIC_ASSERT(dim > 3);
      return elems[3];
    }
  
    /**
     * Extract the two first elements of the vector
     */
    CU_HOST_DEVICE
    Vector<2, T> projectXY(void)
    {
      STATIC_ASSERT(dim > 1);
      return Vector<2, T>(elems[0], elems[1]);
    }
  
    /**
     * Comparison operator.
     *
     * Compare the axis in order.
     */
    CU_HOST_DEVICE
    bool operator<(const Vector& other) const
    {
      for(size_t i = 0; i < dim; ++i) {
        if(elems[i] < other.elems[i])
          return true;
        if(elems[i] > other.elems[i])
          return false;
      }
      return false;
    }
  
    /**
     * Comparison operator.
     *
     * Compare the axis in order.
     */
    CU_HOST_DEVICE
    bool operator<=(const Vector& other) const
    {
      for(size_t i = 0; i < dim; ++i) {
        if(elems[i] < other.elems[i])
          return true;
        if(elems[i] > other.elems[i])
          return false;
      }
      return true;
    }
  
    /**
     * Comparison operator.
     *
     * Compare the axis in order.
     */
    CU_HOST_DEVICE
    bool operator>(const Vector& other) const
    {
      for(size_t i = 0; i < dim; ++i) {
        if(elems[i] > other.elems[i])
          return true;
        if(elems[i] < other.elems[i])
          return false;
      }
      return false;
    }
  
    /**
     * Comparison operator.
     *
     * Compare the axis in order.
     */
    CU_HOST_DEVICE
    bool operator>=(const Vector& other) const
    {
      for(size_t i = 0; i < dim; ++i) {
        if(elems[i] > other.elems[i])
          return true;
        if(elems[i] < other.elems[i])
          return false;
      }
      return true;
    }

    friend std::ostream& operator<<(std::ostream& out, const Vector& vec)
    {
      for(size_t i = 0; i < dim; i++) {
        out << vec.elems[i];
        if(i != (dim - 1))
          out << " ";
      }
      return out;
    }
  
    friend std::istream& operator>>(std::istream& in, Vector& vec)
    {
      in >> vec[0];
      for(size_t i = 1; i < dim && in; i++)
        in >> std::ws >> vec[i];
      return in;
    }
  
  #ifndef __CUDACC__
    friend QTextStream& operator<<(QTextStream& out, const Vector& vec)
    {
      for(size_t i = 0; i < dim; i++) {
        out << vec.elems[i];
        if(i != (dim - 1))
          out << " ";
      }
      return out;
    }
  
    friend QTextStream& operator>>(QTextStream& in, Vector& vec)
    {
      in >> vec[0];
      for(size_t i = 1; i < dim && !in.atEnd(); i++)
        in >> ws >> vec[i];
      return in;
    }
  #endif
  };
  
  /**
   * Cross product \c v1 x \c v2
   *
   * \relates Vector
   */
  template <class T> CU_HOST_DEVICE 
	T operator%(const Vector<2, T>& v1, const Vector<2, T>& v2) 
	{
    return v1 ^ v2;
  }
  
  /**
   * Cross product \c v1 x \c v2 (French notation)
   *
   * \relates Vector
   */
  template <class T> CU_HOST_DEVICE 
	T operator^(const Vector<2, T>& v1, const Vector<2, T>& v2) 
	{
    return ((v1[0] * v2[1]) - (v1[1] * v2[0]));
  }
  
  /**
   * Cross product \c v1 x \c v2 (French notation)
   *
   * \relates Vector
   */
  template <class T> CU_HOST_DEVICE 
	T operator^(const Vector<1, T>&, const Vector<1, T>& ) 
	{ 
		return 0; 
	}
  
  /**
   * Cross product \c v1 x \c v2
   *
   * \relates Vector
   */
  template <class T>
  CU_HOST_DEVICE Vector<3, T> operator%(const Vector<3, T>& v1, const Vector<3, T>& v2)
  {
    return v1 ^ v2;
  }
  
  /**
   * Cross product \c v1 x \c v2 (French notation)
   *
   * \relates Vector
   */
  template <class T> CU_HOST_DEVICE 
	Vector<3, T> operator^(const Vector<3, T>& v1, const Vector<3, T>& v2) {
    Vector<3, T> ans;
    ans[0] = v1[1] * v2[2] - v1[2] * v2[1];
    ans[1] = v1[2] * v2[0] - v1[0] * v2[2];
    ans[2] = v1[0] * v2[1] - v1[1] * v2[0];
    return ans;
  }
  
  /**
   * Angle of the vector with (0,1)
   *
   * \relates Vector
   */
  template <class T> CU_HOST_DEVICE 
	float angle(const Vector<2, T>& v)
  {
    return atan2(v.y(), v.x());
  }
  
  /**
   * Non-oriented angle between \c v1 and \c v2
   *
   * \relates Vector
   */
  template <class T> CU_HOST_DEVICE 
	float angle(const Vector<3, T>& v1, const Vector<3, T>& v2)
  {
    float x = v1 * v2;
    float y = norm(v1 ^ v2);
    return atan2(y, x);
  }
  
  /**
   * Oriented angle between \c v1 and \c v2
   *
   * \relates Vector
   */
  template <class T> CU_HOST_DEVICE 
	float angle(const Vector<2, T>& v1, const Vector<2, T>& v2)
  {
    float x = v1 * v2;
    float y = v1 ^ v2;
    return atan2(y, x);
  }
  
  /**
   * Oriented angle between \c v1 and \c v2
   *
   * \relates Vector
   */
  template <class T> CU_HOST_DEVICE 
	float angle(const Vector<1, T>& v1, const Vector<1, T>& v2)
  {
    return (v1 * v2 < 0) ? -1 : 1;
  }
  
  /**
   * Oriented angle between \c v1 and \c v2 with \c ref to orient the space.
   *
   * \relates Vector
   */
  template <class T> CU_HOST_DEVICE 
	float angle(const Vector<3, T>& v1, const Vector<3, T>& v2, const Vector<3, T>& ref)
  {
    float x = v1 * v2;
    Vector<3, T> n = v1 ^ v2;
    float y = norm(n);
    if(n * ref < 0)
      return atan2(-y, x);
    else
      return atan2(y, x);
  }
  
  /**
   * Normalized real
   *
   * Just return 1
   *
   * \relates Vector
   */
  CU_HOST_DEVICE
  float normalized(float) 
	{
    return 1;
  }
  
  /**
   * Euclidean square norm of a real
   *
   * Just the square of the real
   *
   * \relates Vector
   */
  CU_HOST_DEVICE
  float normsq(float s) 
	{
    return s * s;
  }
  
  /**
   * Euclidean norm of a real
   *
   * Just the absolute value of that real
   *
   * \relates Vector
   */
  CU_HOST_DEVICE
  float norm(float s) 
	{
    return (s < 0) ? -s : s;
  }
  
  /**
   * Function-version of the norm
   *
   * \relates Vector
   * \see Vector::norm()
   */
  template <size_t dim, typename T> CU_HOST_DEVICE T 
	norm(const Vector<dim, T>& v) {
    return v.norm();
  }
  
  /**
   * Function-version of the square norm
   *
   * \relates Vector
   * \see Vector::normsq()
   */
  template <size_t dim, typename T> CU_HOST_DEVICE 
	T normsq(const Vector<dim, T>& v) {
    return v.normsq();
  }
  
  /**
   * Function-version of the square norm
   *
   * \relates Vector
   * \see Vector::normsq()
   */
  template <size_t dim, typename T> CU_HOST_DEVICE 
	Vector<dim, T> normalized(const Vector<dim, T>& v)
  {
    return v.normalized();
  }
  
  using ::fabs;
  
  /**
   * Return the vector whose component is the absolute value of the input vector.
   *
   * \relates Vector
   */
  template <size_t dim, typename T> CU_HOST_DEVICE 
	Vector<dim, T> fabs(const Vector<dim, T>& v)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = fabs(v[i]);
    }
    return result;
  }
  
  /**
   * Return the vector whose component is the max of the two input vectors
   * components.
   *
   * \relates Vector
   */
  template <size_t dim, typename T> CU_HOST_DEVICE 
	Vector<dim, T> max(const Vector<dim, T>& v1, const Vector<dim, T>& v2)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = (v1[i] >= v2[i] ? v1[i] : v2[i]);
    }
    return result;
  }
  
  /**
   * Return the vector whose component is the min of the two input vectors
   * components.
   *
   * \relates Vector
   */
  template <size_t dim, typename T> CU_HOST_DEVICE 
	Vector<dim, T> min(const Vector<dim, T>& v1, const Vector<dim, T>& v2)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = (v1[i] <= v2[i] ? v1[i] : v2[i]);
    }
    return result;
  }
  
  /**
   * Return the vector with components clipped to min and max vectors
   *
   * \relates Vector
   */
  template <size_t dim, typename T>
  CU_HOST_DEVICE Vector<dim, T> trim(const Vector<dim, T>& v, const Vector<dim, T>& minv, const Vector<dim, T>& maxv)
  {
    Vector<dim, T> result;
    result = max(minv, min(maxv, v));
    return result;
  }
  
  /**
   * Return the vector whose component is the product of the two input vectors
   * components.
   *
   * \relates Vector
   */
  template <size_t dim, typename T>
  CU_HOST_DEVICE Vector<dim, T> multiply(const Vector<dim, T>& v1, const Vector<dim, T>& v2)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = v1[i] * v2[i];
    }
    return result;
  }
  
  //  /**
  //   * Return the vector whose components are clipped to min/max
  //   * components.
  //   *
  //   * \relates Vector
  //   */
  //  template <size_t dim, typename T>
  //  CU_HOST_DEVICE
  //  Vector<dim,T> multiply( const Vector<dim,T>& v, const Vector<dim,T>& minv, const Vector<dim,T>& maxv)
  //  {
  //    Vector<dim,T> result;
  //    result = min(minv, max(maxv, v));
  //    return result;
  //  }
  
  /**
   * Return the vector whose component is the ratio of the two input vectors
   * components.
   *
   * \relates Vector
   */
  template <size_t dim, typename T>
  CU_HOST_DEVICE Vector<dim, T> divide(const Vector<dim, T>& v1, const Vector<dim, T>& v2)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = v1[i] / v2[i];
    }
    return result;
  }
  
  /**
   * Find a vector orthogonal to v
   *
   * \relates Vector
   */
  template <typename T> CU_HOST_DEVICE 
	Vector<3, T> orthogonal(const Vector<3, T>& v)
  {
    const float ratio = 1 - 1e-8;
    if((std::abs(v.y()) >= ratio * std::abs(v.x())) && (std::abs(v.z()) >= ratio * std::abs(v.x())))
      return Vector<3, T>(0, -v.z(), v.y());
    else if((std::abs(v.x()) >= ratio * std::abs(v.y())) && (std::abs(v.z()) >= ratio * std::abs(v.y())))
      return Vector<3, T>(-v.z(), 0, v.x());
    else
      return Vector<3, T>(-v.y(), v.x(), 0);
  }
  
  /*
   * Transform all components of a vector with a function
   *
   * \relates Vector
   */
  template <size_t dim, typename T, typename T1>
  CU_HOST_DEVICE Vector<dim, T> map(const T& (*fct)(const T1 &), const Vector<dim, T1>& v)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = (*fct)(v[i]);
    }
    return result;
  }
  
  /*
   * Transform all components of a vector with a function
   *
   * \relates Vector
   */
  template <size_t dim, typename T, typename T1>
  CU_HOST_DEVICE Vector<dim, T> map(T (*fct)(const T1&), const Vector<dim, T1>& v)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = (*fct)(v[i]);
    }
    return result;
  }
  
  /**
   * Map a function to all elements of a vector
   *
   * \relates Vector
   */
  template <size_t dim, typename T, typename T1> CU_HOST_DEVICE 
	Vector<dim, T> map(T (*fct)(T1), const Vector<dim, T1>& v)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = (*fct)(v[i]);
    }
    return result;
  }
  
  /**
   * Map a function to all elements of a vector
   *
   * \relates Vector
   */
  template <size_t dim, typename T>
  CU_HOST_DEVICE Vector<dim, T> map(const T& (*fct)(const T &, const T &), const Vector<dim, T>& v1,
                                    const Vector<dim, T>& v2)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = (*fct)(v1[i], v2[i]);
    }
    return result;
  }
  
  /**
   * Map a function to all elements of a vector
   *
   * \relates Vector
   */
  template <size_t dim, typename T>
  CU_HOST_DEVICE Vector<dim, T> map(T (*fct)(const T&, const T&), const Vector<dim, T>& v1, const Vector<dim, T>& v2)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = (*fct)(v1[i], v2[i]);
    }
    return result;
  }
  
  /**
   * Map a function to all elements of a vector
   *
   * \relates Vector
   */
  template <size_t dim, typename T>
  CU_HOST_DEVICE Vector<dim, T> map(T (*fct)(T, T), const Vector<dim, T>& v1, const Vector<dim, T>& v2)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = (*fct)(v1[i], v2[i]);
    }
    return result;
  }
  
  /**
   * Map a function to all elements of a vector
   *
   * \relates Vector
   */
  template <size_t dim, typename T, typename T1, typename T2>
  CU_HOST_DEVICE Vector<dim, T> map(const T& (*fct)(const T1 &, const T2 &), const Vector<dim, T1>& v1,
                                    const Vector<dim, T2>& v2)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = (*fct)(v1[i], v2[i]);
    }
    return result;
  }
  
  /**
   * Map a function to all elements of a vector
   *
   * \relates Vector
   */
  template <size_t dim, typename T, typename T1, typename T2> CU_HOST_DEVICE 
	Vector<dim, T> map(T (*fct)(const T1&, const T2&), const Vector<dim, T1>& v1, 
									                                                 const Vector<dim, T2>& v2)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = (*fct)(v1[i], v2[i]);
    }
    return result;
  }
  
  /**
   * Map a function to all elements of a vector
   *
   * \relates Vector
   */
  template <size_t dim, typename T, typename T1, typename T2> CU_HOST_DEVICE 
  Vector<dim, T> map(T (*fct)(T1, T2), const Vector<dim, T1>& v1, const Vector<dim, T2>& v2)
  {
    Vector<dim, T> result;
    for(size_t i = 0; i < dim; ++i) {
      result[i] = (*fct)(v1[i], v2[i]);
    }
    return result;
  }
  
  /**
   * Add a value to all elements of a vector
   *
   * \relates Vector
   */
  template <size_t dim, typename T> CU_HOST_DEVICE 
  Vector<dim, T> operator+(const Vector<dim, T>& v, const T& value)
  {
    Vector<dim, T> res;
    for(size_t i = 0; i < dim; ++i)
      res[i] = v[i] + value;
    return res;
  }
  
  /**
   * Add a value to all elements of a vector
   *
   * \relates Vector
   */
  template <size_t dim, typename T> CU_HOST_DEVICE 
	Vector<dim, T> operator+(const T& value, const Vector<dim, T>& v)
  {
    Vector<dim, T> res;
    for(size_t i = 0; i < dim; ++i)
      res[i] = v[i] + value;
    return res;
  }
  
  /**
   * Substact a value to all elements of a vector
   *
   * \relates Vector
   */
  template <size_t dim, typename T> CU_HOST_DEVICE 
	Vector<dim, T> operator-(const Vector<dim, T>& v, const T& value)
  {
    Vector<dim, T> res;
    for(size_t i = 0; i < dim; ++i)
      res[i] = v[i] - value;
    return res;
  }
  
  /**
   * Equivalent to substracting a vector with all component the same to another
   * one.
   *
   * \relates Vector
   */
  template <size_t dim, typename T> CU_HOST_DEVICE 
	Vector<dim, T> operator-(const T& value, const Vector<dim, T>& v)
  {
    Vector<dim, T> res;
    for(size_t i = 0; i < dim; ++i)
      res[i] = value - v[i];
    return res;
  }
  
  /**
   * Create a homogeneous coordinate vector from a cartesian one.
   *
   * This adds one dimension, and set it to 1.
   *
   * \relates Vector
   */
  template <size_t dim, typename T> CU_HOST_DEVICE 
	Vector<dim + 1, T> homogeneous(const Vector<dim, T>& v)
  {
    Vector<dim + 1, T> res(v);
    res[dim] = T(1);
    return res;
  }
  
  /**
   * Extract the cartesion coordinates from a homogeneous vector.
   *
   * \warning The homogenous coordinate shouldn't be zero or there will be a division by 0!
   *
   * \relates Vector
   */
  template <size_t dim, typename T> CU_HOST_DEVICE 
	Vector<dim - 1, T> cartesian(const Vector<dim, T>& v)
  {
    Vector<dim - 1, T> res(v);
    res /= v[dim - 1];
    return res;
  }
}
  
#ifndef __CUDACC__
// Add hashes for STL
#include <Common.hpp>
#include <boost/functional/hash.hpp>
  
namespace std 
{
  HASH_NS_ENTER
  template <size_t dim, typename T> struct hash<mgx::Vector<dim, T> > {
    size_t operator()(const mgx::Vector<dim, T>& v) const
    {
      std::size_t seed = 0;
      for(size_t i = 1; i < dim; ++i)
        boost::hash_combine(seed, v[i]);
      return seed;
//
//      boost::hash_combine(seed, 2);
//      static const hash<T> hash_value = hash<T>();
//      size_t acc = hash_value(v[0]);
//      for(size_t i = 1; i < dim; ++i)
//        acc ^= hash_value(v[i]);
//      return acc;
    }
  };
  HASH_NS_EXIT
}

#endif

#endif
