/** \file vector.cpp
 *
 * Function implementations for the class vector.
 *
 */


#include "SingleVector.hpp"
#include "Helpers/Assert.hpp"

#include <iostream>

using namespace std;


/************************************ Functions for class vector ************************************/

/** Constructor of class vector.
 */
Vector::Vector() :
  rep(new SingleVector())
{};

Vector::Vector(Baseconstructor)  // used by derived objects to construct their bases
{}

Vector::Vector(Baseconstructor,const Vector* v) :
  rep(v->clone())
{}

Vector Vector::VecFromRep(const Vector* v){
  return Vector(Baseconstructor(),v);
}

/** Constructor of class vector.
 */
Vector::Vector(const double x1, const double x2, const double x3) :
  rep(new SingleVector(x1,x2,x3))
{};

/**
 * Copy constructor
 */
Vector::Vector(const Vector& src) :
  rep(src.rep->clone())
{}

/**
 * Assignment operator
 */
Vector& Vector::operator=(const Vector& src){
  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
  // check for self assignment
  if(&src!=this){
    rep.reset(src.rep->clone());
  }
  return *this;
}

/** Desctructor of class vector.
 */
Vector::~Vector() {};

/** Calculates square of distance between this and another vector.
 * \param *y array to second vector
 * \return \f$| x - y |^2\f$
 */
double Vector::DistanceSquared(const Vector &y) const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->DistanceSquared(y);
};

/** Calculates distance between this and another vector.
 * \param *y array to second vector
 * \return \f$| x - y |\f$
 */
double Vector::Distance(const Vector &y) const
{
  return (sqrt(DistanceSquared(y)));
};

/** Calculates distance between this and another vector in a periodic cell.
 * \param *y array to second vector
 * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
 * \return \f$| x - y |\f$
 */
double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->PeriodicDistance(y,cell_size);
};

/** Calculates distance between this and another vector in a periodic cell.
 * \param *y array to second vector
 * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
 * \return \f$| x - y |^2\f$
 */
double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->PeriodicDistanceSquared(y,cell_size);
};

/** Keeps the vector in a periodic cell, defined by the symmetric \a *matrix.
 * \param *out ofstream for debugging messages
 * Tries to translate a vector into each adjacent neighbouring cell.
 */
void Vector::KeepPeriodic(const double * const matrix)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  rep->KeepPeriodic(matrix);
};

/** Calculates scalar product between this and another vector.
 * \param *y array to second vector
 * \return \f$\langle x, y \rangle\f$
 */
double Vector::ScalarProduct(const Vector &y) const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->ScalarProduct(y);
};


/** Calculates VectorProduct between this and another vector.
 *  -# returns the Product in place of vector from which it was initiated
 *  -# ATTENTION: Only three dim.
 *  \param *y array to vector with which to calculate crossproduct
 *  \return \f$ x \times y \f&
 */
void Vector::VectorProduct(const Vector &y)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  rep->VectorProduct(y);
};


/** projects this vector onto plane defined by \a *y.
 * \param *y normal vector of plane
 * \return \f$\langle x, y \rangle\f$
 */
void Vector::ProjectOntoPlane(const Vector &y)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  rep->ProjectOntoPlane(y);
};

/** Calculates the minimum distance of this vector to the plane.
 * \param *out output stream for debugging
 * \param *PlaneNormal normal of plane
 * \param *PlaneOffset offset of plane
 * \return distance to plane
 */
double Vector::DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->DistanceToPlane(PlaneNormal,PlaneOffset);
};

/** Calculates the projection of a vector onto another \a *y.
 * \param *y array to second vector
 */
void Vector::ProjectIt(const Vector &y)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  rep->ProjectIt(y);
};

/** Calculates the projection of a vector onto another \a *y.
 * \param *y array to second vector
 * \return Vector
 */
Vector Vector::Projection(const Vector &y) const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->Projection(y);
};

/** Calculates norm of this vector.
 * \return \f$|x|\f$
 */
double Vector::Norm() const
{
  return (sqrt(NormSquared()));
};

/** Calculates squared norm of this vector.
 * \return \f$|x|^2\f$
 */
double Vector::NormSquared() const
{
  return (ScalarProduct(*this));
};

/** Normalizes this vector.
 */
void Vector::Normalize()
{
  double factor = Norm();
  (*this) *= 1/factor;
};

/** Zeros all components of this vector.
 */
void Vector::Zero()
{
  rep.reset(new SingleVector());
};

/** Zeros all components of this vector.
 */
void Vector::One(const double one)
{
  rep.reset(new SingleVector(one,one,one));
};

/** Checks whether vector has all components zero.
 * @return true - vector is zero, false - vector is not
 */
bool Vector::IsZero() const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->IsZero();
};

/** Checks whether vector has length of 1.
 * @return true - vector is normalized, false - vector is not
 */
bool Vector::IsOne() const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->IsOne();
};

/** Checks whether vector is normal to \a *normal.
 * @return true - vector is normalized, false - vector is not
 */
bool Vector::IsNormalTo(const Vector &normal) const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->IsNormalTo(normal);
};

/** Checks whether vector is normal to \a *normal.
 * @return true - vector is normalized, false - vector is not
 */
bool Vector::IsEqualTo(const Vector &a) const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->IsEqualTo(a);
};

/** Calculates the angle between this and another vector.
 * \param *y array to second vector
 * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
 */
double Vector::Angle(const Vector &y) const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->Angle(y);
};


double& Vector::operator[](size_t i){
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return (*rep)[i];
}

const double& Vector::operator[](size_t i) const{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return (*rep)[i];
}

double& Vector::at(size_t i){
  return (*this)[i];
}

const double& Vector::at(size_t i) const{
  return (*this)[i];
}

double* Vector::get(){
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->get();
}

/** Compares vector \a to vector \a b component-wise.
 * \param a base vector
 * \param b vector components to add
 * \return a == b
 */
bool Vector::operator==(const Vector& b) const
{
  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
  return IsEqualTo(b);
};

/** Sums vector \a to this lhs component-wise.
 * \param a base vector
 * \param b vector components to add
 * \return lhs + a
 */
const Vector& Vector::operator+=(const Vector& b)
{
  this->AddVector(b);
  return *this;
};

/** Subtracts vector \a from this lhs component-wise.
 * \param a base vector
 * \param b vector components to add
 * \return lhs - a
 */
const Vector& Vector::operator-=(const Vector& b)
{
  this->SubtractVector(b);
  return *this;
};

/** factor each component of \a a times a double \a m.
 * \param a base vector
 * \param m factor
 * \return lhs.x[i] * m
 */
const Vector& operator*=(Vector& a, const double m)
{
  a.Scale(m);
  return a;
};

/** Sums two vectors \a  and \b component-wise.
 * \param a first vector
 * \param b second vector
 * \return a + b
 */
Vector const Vector::operator+(const Vector& b) const
{
  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
  Vector x = *this;
  x.AddVector(b);
  return x;
};

/** Subtracts vector \a from \b component-wise.
 * \param a first vector
 * \param b second vector
 * \return a - b
 */
Vector const Vector::operator-(const Vector& b) const
{
  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
  Vector x = *this;
  x.SubtractVector(b);
  return x;
};

/** Factors given vector \a a times \a m.
 * \param a vector
 * \param m factor
 * \return m * a
 */
Vector const operator*(const Vector& a, const double m)
{
  Vector x(a);
  x.Scale(m);
  return x;
};

/** Factors given vector \a a times \a m.
 * \param m factor
 * \param a vector
 * \return m * a
 */
Vector const operator*(const double m, const Vector& a )
{
  Vector x(a);
  x.Scale(m);
  return x;
};

ostream& operator<<(ostream& ost, const Vector& m)
{
  ost << "(";
  for (int i=0;i<NDIM;i++) {
    ost << m[i];
    if (i != 2)
      ost << ",";
  }
  ost << ")";
  return ost;
};


void Vector::ScaleAll(const double *factor)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  rep->ScaleAll(factor);
};



void Vector::Scale(const double factor)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  rep->Scale(factor);
};

/** Given a box by its matrix \a *M and its inverse *Minv the vector is made to point within that box.
 * \param *M matrix of box
 * \param *Minv inverse matrix
 */
void Vector::WrapPeriodically(const double * const M, const double * const Minv)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  rep->WrapPeriodically(M,Minv);
};

/** Do a matrix multiplication.
 * \param *matrix NDIM_NDIM array
 */
void Vector::MatrixMultiplication(const double * const M)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  rep->MatrixMultiplication(M);
};

/** Do a matrix multiplication with the \a *A' inverse.
 * \param *matrix NDIM_NDIM array
 */
bool Vector::InverseMatrixMultiplication(const double * const A)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->InverseMatrixMultiplication(A);
};


/** Creates this vector as the b y *factors' components scaled linear combination of the given three.
 * this vector = x1*factors[0] + x2* factors[1] + x3*factors[2]
 * \param *x1 first vector
 * \param *x2 second vector
 * \param *x3 third vector
 * \param *factors three-component vector with the factor for each given vector
 */
void Vector::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
{
  (*this) = (factors[0]*x1) +
            (factors[1]*x2) +
            (factors[2]*x3);
};

/** Mirrors atom against a given plane.
 * \param n[] normal vector of mirror plane.
 */
void Vector::Mirror(const Vector &n)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  rep->Mirror(n);
};


/** Calculates orthonormal vector to one given vector.
 * Just subtracts the projection onto the given vector from this vector.
 * The removed part of the vector is Vector::Projection()
 * \param *x1 vector
 * \return true - success, false - vector is zero
 */
bool Vector::MakeNormalTo(const Vector &y1)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->MakeNormalTo(y1);
};

/** Creates this vector as one of the possible orthonormal ones to the given one.
 * Just scan how many components of given *vector are unequal to zero and
 * try to get the skp of both to be zero accordingly.
 * \param *vector given vector
 * \return true - success, false - failure (null vector given)
 */
bool Vector::GetOneNormalVector(const Vector &GivenVector)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->GetOneNormalVector(GivenVector);
};

/** Adds vector \a *y componentwise.
 * \param *y vector
 */
void Vector::AddVector(const Vector &y)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  rep->AddVector(y);
}

/** Adds vector \a *y componentwise.
 * \param *y vector
 */
void Vector::SubtractVector(const Vector &y)
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  rep->SubtractVector(y);
}

/**
 * Checks whether this vector is within the parallelepiped defined by the given three vectors and
 * their offset.
 *
 * @param offest for the origin of the parallelepiped
 * @param three vectors forming the matrix that defines the shape of the parallelpiped
 */
bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
{
  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
  return rep->IsInParallelepiped(offset, parallelepiped);
}

bool Vector::isBaseClass() const{
  return true;
}

Vector* Vector::clone() const{
  ASSERT(false, "Cannot clone a base Vector object");
  return 0;
}