source: src/vector.cpp@ b2531f

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

#include "Helpers/MemDebug.hpp"

#include "vector.hpp"
#include "verbose.hpp"
#include "World.hpp"
#include "Helpers/Assert.hpp"
#include "Helpers/fast_functions.hpp"

#include <iostream>

using namespace std;


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

/** Constructor of class vector.
 */
Vector::Vector()
{
  content = gsl_vector_calloc (NDIM);
};

/**
 * Copy constructor
 */

Vector::Vector(const Vector& src)
{
  content = gsl_vector_alloc(NDIM);
  gsl_vector_set(content,0,src[0]);
  gsl_vector_set(content,1,src[1]);
  gsl_vector_set(content,2,src[2]);
}

/** Constructor of class vector.
 */
Vector::Vector(const double x1, const double x2, const double x3)
{
  content = gsl_vector_alloc(NDIM);
  gsl_vector_set(content,0,x1);
  gsl_vector_set(content,1,x2);
  gsl_vector_set(content,2,x3);
};

/**
 * Assignment operator
 */
Vector& Vector::operator=(const Vector& src){
  // check for self assignment
  if(&src!=this){
    gsl_vector_set(content,0,src[0]);
    gsl_vector_set(content,1,src[1]);
    gsl_vector_set(content,2,src[2]);
  }
  return *this;
}

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

/** 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
{
  double res = 0.;
  for (int i=NDIM;i--;)
    res += (at(i)-y[i])*(at(i)-y[i]);
  return (res);
};

/** 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)));
};

Vector Vector::getClosestPoint(const Vector &point) const{
  // the closest point to a single point space is always the single point itself
  return *this;
}

/** 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
{
  double res = distance(y), tmp, matrix[NDIM*NDIM];
    Vector Shiftedy, TranslationVector;
    int N[NDIM];
    matrix[0] = cell_size[0];
    matrix[1] = cell_size[1];
    matrix[2] = cell_size[3];
    matrix[3] = cell_size[1];
    matrix[4] = cell_size[2];
    matrix[5] = cell_size[4];
    matrix[6] = cell_size[3];
    matrix[7] = cell_size[4];
    matrix[8] = cell_size[5];
    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    for (N[0]=-1;N[0]<=1;N[0]++)
      for (N[1]=-1;N[1]<=1;N[1]++)
        for (N[2]=-1;N[2]<=1;N[2]++) {
          // create the translation vector
          TranslationVector.Zero();
          for (int i=NDIM;i--;)
            TranslationVector[i] = (double)N[i];
          TranslationVector.MatrixMultiplication(matrix);
          // add onto the original vector to compare with
          Shiftedy = y + TranslationVector;
          // get distance and compare with minimum so far
          tmp = distance(Shiftedy);
          if (tmp < res) res = tmp;
        }
    return (res);
};

/** 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
{
  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
    Vector Shiftedy, TranslationVector;
    int N[NDIM];
    matrix[0] = cell_size[0];
    matrix[1] = cell_size[1];
    matrix[2] = cell_size[3];
    matrix[3] = cell_size[1];
    matrix[4] = cell_size[2];
    matrix[5] = cell_size[4];
    matrix[6] = cell_size[3];
    matrix[7] = cell_size[4];
    matrix[8] = cell_size[5];
    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    for (N[0]=-1;N[0]<=1;N[0]++)
      for (N[1]=-1;N[1]<=1;N[1]++)
        for (N[2]=-1;N[2]<=1;N[2]++) {
          // create the translation vector
          TranslationVector.Zero();
          for (int i=NDIM;i--;)
            TranslationVector[i] = (double)N[i];
          TranslationVector.MatrixMultiplication(matrix);
          // add onto the original vector to compare with
          Shiftedy = y + TranslationVector;
          // get distance and compare with minimum so far
          tmp = DistanceSquared(Shiftedy);
          if (tmp < res) res = tmp;
        }
    return (res);
};

/** 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)
{
  //  int N[NDIM];
  //  bool flag = false;
    //vector Shifted, TranslationVector;
  //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
  //  Log() << Verbose(2) << "Vector is: ";
  //  Output(out);
  //  Log() << Verbose(0) << endl;
    InverseMatrixMultiplication(matrix);
    for(int i=NDIM;i--;) { // correct periodically
      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
        at(i) += ceil(at(i));
      } else {
        at(i) -= floor(at(i));
      }
    }
    MatrixMultiplication(matrix);
  //  Log() << Verbose(2) << "New corrected vector is: ";
  //  Output(out);
  //  Log() << Verbose(0) << endl;
  //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
};

/** 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
{
  double res = 0.;
  for (int i=NDIM;i--;)
    res += at(i)*y[i];
  return (res);
};


/** 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)
{
  Vector tmp;
  for(int i=NDIM;i--;)
    tmp[i] = at((i+1)%NDIM)*y[(i+2)%NDIM] - at((i+2)%NDIM)*y[(i+1)%NDIM];
  (*this) = tmp;
};


/** 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)
{
  Vector tmp;
  tmp = y;
  tmp.Normalize();
  tmp.Scale(ScalarProduct(tmp));
  *this -= tmp;
};

/** Calculates the minimum distance of this vector to the plane.
 * \sa Vector::GetDistanceVectorToPlane()
 * \param *out output stream for debugging
 * \param *PlaneNormal normal of plane
 * \param *PlaneOffset offset of plane
 * \return distance to plane
 */
double Vector::DistanceToSpace(const Space &space) const
{
  return space.distance(*this);
};

/** Calculates the projection of a vector onto another \a *y.
 * \param *y array to second vector
 */
void Vector::ProjectIt(const Vector &y)
{
  (*this) += (-ScalarProduct(y))*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
{
  Vector helper = y;
  helper.Scale((ScalarProduct(y)/y.NormSquared()));

  return helper;
};

/** 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()
{
  at(0)=at(1)=at(2)=0;
};

/** Zeros all components of this vector.
 */
void Vector::One(const double one)
{
  at(0)=at(1)=at(2)=one;
};

/** Checks whether vector has all components zero.
 * @return true - vector is zero, false - vector is not
 */
bool Vector::IsZero() const
{
  return (fabs(at(0))+fabs(at(1))+fabs(at(2)) < MYEPSILON);
};

/** Checks whether vector has length of 1.
 * @return true - vector is normalized, false - vector is not
 */
bool Vector::IsOne() const
{
  return (fabs(Norm() - 1.) < MYEPSILON);
};

/** Checks whether vector is normal to \a *normal.
 * @return true - vector is normalized, false - vector is not
 */
bool Vector::IsNormalTo(const Vector &normal) const
{
  if (ScalarProduct(normal) < MYEPSILON)
    return true;
  else
    return false;
};

/** Checks whether vector is normal to \a *normal.
 * @return true - vector is normalized, false - vector is not
 */
bool Vector::IsEqualTo(const Vector &a) const
{
  bool status = true;
  for (int i=0;i<NDIM;i++) {
    if (fabs(at(i) - a[i]) > MYEPSILON)
      status = false;
  }
  return status;
};

/** 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
{
  double norm1 = Norm(), norm2 = y.Norm();
  double angle = -1;
  if ((fabs(norm1) > MYEPSILON) && (fabs(norm2) > MYEPSILON))
    angle = this->ScalarProduct(y)/norm1/norm2;
  // -1-MYEPSILON occured due to numerical imprecision, catch ...
  //Log() << Verbose(2) << "INFO: acos(-1) = " << acos(-1) << ", acos(-1+MYEPSILON) = " << acos(-1+MYEPSILON) << ", acos(-1-MYEPSILON) = " << acos(-1-MYEPSILON) << "." << endl;
  if (angle < -1)
    angle = -1;
  if (angle > 1)
    angle = 1;
  return acos(angle);
};


double& Vector::operator[](size_t i){
  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
  return *gsl_vector_ptr (content, i);
}

const double& Vector::operator[](size_t i) const{
  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
  return *gsl_vector_ptr (content, i);
}

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

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

gsl_vector* Vector::get(){
  return content;
}

/** 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
{
  return IsEqualTo(b);
};

bool Vector::operator!=(const Vector& b) const
{
  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
{
  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
{
  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)
{
  for (int i=NDIM;i--;)
    at(i) *= factor[i];
};



void Vector::Scale(const double factor)
{
  for (int i=NDIM;i--;)
    at(i) *= 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)
{
  MatrixMultiplication(Minv);
  // truncate to [0,1] for each axis
  for (int i=0;i<NDIM;i++) {
    //at(i) += 0.5;  // set to center of box
    while (at(i) >= 1.)
      at(i) -= 1.;
    while (at(i) < 0.)
      at(i) += 1.;
  }
  MatrixMultiplication(M);
};

std::pair<Vector,Vector> Vector::partition(const Vector &rhs) const{
  double factor = ScalarProduct(rhs)/rhs.NormSquared();
  Vector res= factor * rhs;
  return make_pair(res,(*this)-res);
}

std::pair<pointset,Vector> Vector::partition(const pointset &points) const{
  Vector helper = *this;
  pointset res;
  for(pointset::const_iterator iter=points.begin();iter!=points.end();++iter){
    pair<Vector,Vector> currPart = helper.partition(*iter);
    res.push_back(currPart.first);
    helper = currPart.second;
  }
  return make_pair(res,helper);
}

/** Do a matrix multiplication.
 * \param *matrix NDIM_NDIM array
 */
void Vector::MatrixMultiplication(const double * const M)
{
  Vector tmp;
  // do the matrix multiplication
  for(int i=NDIM;i--;)
    tmp[i] = M[i]*at(0)+M[i+3]*at(1)+M[i+6]*at(2);

  (*this) = tmp;
};

/** Do a matrix multiplication with the \a *A' inverse.
 * \param *matrix NDIM_NDIM array
 */
bool Vector::InverseMatrixMultiplication(const double * const A)
{
  double B[NDIM*NDIM];
  double detA = RDET3(A);
  double detAReci;

  // calculate the inverse B
  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
    detAReci = 1./detA;
    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33

    MatrixMultiplication(B);

    return true;
  } else {
    return false;
  }
};


/** 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);
};

/** Calculates orthonormal vector to one given vectors.
 * 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)
{
  bool result = false;
  double factor = y1.ScalarProduct(*this)/y1.NormSquared();
  Vector x1 = factor * y1;
  SubtractVector(x1);
  for (int i=NDIM;i--;)
    result = result || (fabs(at(i)) > MYEPSILON);

  return result;
};

/** 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)
{
  int Components[NDIM]; // contains indices of non-zero components
  int Last = 0;   // count the number of non-zero entries in vector
  int j;  // loop variables
  double norm;

  for (j=NDIM;j--;)
    Components[j] = -1;

  // in two component-systems we need to find the one position that is zero
  int zeroPos = -1;
  // find two components != 0
  for (j=0;j<NDIM;j++){
    if (fabs(GivenVector[j]) > MYEPSILON)
      Components[Last++] = j;
    else
      // this our zero Position
      zeroPos = j;
  }

  switch(Last) {
    case 3:  // threecomponent system
      // the position of the zero is arbitrary in three component systems
      zeroPos = Components[2];
    case 2:  // two component system
      norm = sqrt(1./(GivenVector[Components[1]]*GivenVector[Components[1]]) + 1./(GivenVector[Components[0]]*GivenVector[Components[0]]));
      at(zeroPos) = 0.;
      // in skp both remaining parts shall become zero but with opposite sign and third is zero
      at(Components[1]) = -1./GivenVector[Components[1]] / norm;
      at(Components[0]) = 1./GivenVector[Components[0]] / norm;
      return true;
      break;
    case 1: // one component system
      // set sole non-zero component to 0, and one of the other zero component pendants to 1
      at((Components[0]+2)%NDIM) = 0.;
      at((Components[0]+1)%NDIM) = 1.;
      at(Components[0]) = 0.;
      return true;
      break;
    default:
      return false;
  }
};

/** Adds vector \a *y componentwise.
 * \param *y vector
 */
void Vector::AddVector(const Vector &y)
{
  for(int i=NDIM;i--;)
    at(i) += y[i];
}

/** Adds vector \a *y componentwise.
 * \param *y vector
 */
void Vector::SubtractVector(const Vector &y)
{
  for(int i=NDIM;i--;)
    at(i) -= y[i];
}

/**
 * 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
{
  Vector a = (*this)-offset;
  a.InverseMatrixMultiplication(parallelepiped);
  bool isInside = true;

  for (int i=NDIM;i--;)
    isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));

  return isInside;
}


// some comonly used vectors
const Vector zeroVec(0,0,0);
const Vector e1(1,0,0);
const Vector e2(0,1,0);
const Vector e3(0,0,1);