source: src/LinearAlgebra/VectorContent.cpp@ 9b410d

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * VectorContent.cpp
 *
 *  Created on: Nov 15, 2010
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "CodePatterns/MemDebug.hpp"

#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <cmath>
#include <iostream>
#include <limits>

#include "CodePatterns/Assert.hpp"
#include "Helpers/defs.hpp"
#include "LinearAlgebra/defs.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/VectorContent.hpp"

/** Constructor of class VectorContent.
 * Allocates GSL structures
 * \param _dim number of components
 */
VectorContent::VectorContent(size_t _dim) :
    dimension(_dim)
{
  content = gsl_vector_calloc(dimension);
}

/** Constructor of class VectorContent.
 * We need this VectorBaseCase for the VectorContentView class.
 * There no content should be allocated, as it is just a view with an internal
 * gsl_vector_view. Hence, VectorBaseCase is just dummy class to give the
 * constructor a unique signature.
 * \param VectorBaseCase
 */
VectorContent::VectorContent(VectorBaseCase)
{}

/** Copy constructor of class VectorContent.
 * Allocates GSL structures and copies components from \a *src.
 * \param *src source vector
 */
VectorContent::VectorContent(const VectorContent * const src) :
    dimension(src->dimension)
{
  content = gsl_vector_alloc(dimension);
  gsl_vector_memcpy (content, src->content);
};

/** Copy constructor of class VectorContent.
 * Allocates GSL structures and copies components from \a *src.
 * \param *src source vector
 */
VectorContent::VectorContent(const VectorContent & src) :
    dimension(src.dimension)
{
  content = gsl_vector_alloc(dimension);
  gsl_vector_memcpy (content, src.content);
};

/** Copy constructor of class VectorContent.
 * No allocation, we just take over gsl_vector.
 * \param *src source gsl_vector
 */
VectorContent::VectorContent(gsl_vector * _src) :
  dimension(_src->size)
{
  content = _src;
}

/** Destructor of class VectorContent.
 * Frees GSL structures
 */
VectorContent::~VectorContent()
{
  if(content != NULL){
    gsl_vector_free(content);
    content = NULL;
  }
}

/* ============================ Accessing =============================== */
/** Accessor for manipulating component (i).
 * \param i component number
 * \return reference to component (i)
 */
double &VectorContent::at(size_t i)
{
  ASSERT((i>=0) && (i<dimension),
      "VectorContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(dimension)+"]");
  return *gsl_vector_ptr (content, i);
}

/** Constant accessor for (value of) component (i).
 * \param i component number
 * \return const component (i)
 */
const double VectorContent::at(size_t i) const
{
  ASSERT((i>=0) && (i<dimension),
      "VectorContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(dimension)+"]");
  return gsl_vector_get(content, i);
}

/** These functions return a pointer to the \a m-th element of a vector.
 *  If \a m lies outside the allowed range of 0 to VectorContent::dimension-1 then the error handler is invoked and a null pointer is returned.
 * \param m m-th element
 * \return pointer to \a m-th element
 */
double *VectorContent::Pointer(size_t m) const
{
  return gsl_vector_ptr (content, m);
};

/** These functions return a constant pointer to the \a m-th element of a vector.
 *  If \a m lies outside the allowed range of 0 to VectorContent::dimension-1 then the error handler is invoked and a null pointer is returned.
 * \param m m-th element
 * \return const pointer to \a m-th element
 */
const double *VectorContent::const_Pointer(size_t m) const
{
  return gsl_vector_const_ptr (content, m);
};

/** Assignment operator.
 * \param &src source vector to assign \a *this to
 * \return reference to \a *this
 */
VectorContent& VectorContent::operator=(const VectorContent& src)
{
  ASSERT(dimension == src.dimension,
      "Dimensions have to be the same to assign VectorContent onto another:"
      +toString(dimension)+" != "+toString(src.dimension)+"!");
  // check for self assignment
  if(&src!=this){
    gsl_vector_memcpy(content, src.content);
  }
  return *this;
}

/* ========================== Initializing =============================== */
/** This function sets all the elements of the vector to the value \a x.
 * \param x value to set to
 */
void VectorContent::setValue(double x)
{
  gsl_vector_set_all (content, x);
};

/** This function sets the vector from a double array.
 * Creates a vector view of the array and performs a memcopy.
 * \param *x array of values (no dimension check is performed)
 */
void VectorContent::setFromDoubleArray(double * x)
{
  gsl_vector_view m = gsl_vector_view_array (x, dimension);
  gsl_vector_memcpy (content, &m.vector);
};

/**
 * This function sets the GSLvector from an ordinary vector.
 *
 * Takes access to the internal gsl_vector and copies it
 */
void VectorContent::setFromVector(Vector &v)
{
  gsl_vector_memcpy (content, v.get()->content);
}

/** This function sets all the elements of the vector to zero.
 */
void VectorContent::setZero()
{
  gsl_vector_set_zero (content);
};

/** This function makes a basis vector by setting all the elements of the vector to zero except for the i-th element which is set to one.
 * \param i i-th component to set to unity (all other to zero)
 * \return vector set
 */
int VectorContent::setBasis(size_t i)
{
  return gsl_vector_set_basis (content, i);
};

/* ====================== Exchanging elements ============================ */
/** This function exchanges the \a i-th and \a j-th elements of the vector in-place.
 * \param i i-th element to swap with ...
 * \param j ... j-th element to swap against
 */
int VectorContent::SwapElements(size_t i, size_t j)
{
  return gsl_vector_swap_elements (content, i, j);
};

/** This function reverses the order of the elements of the vector.
 */
int VectorContent::Reverse()
{
  return gsl_vector_reverse (content);
};


/* ========================== Operators =============================== */
/** Accessor for manipulating component (i).
 * \param i component number
 * \return reference to component (i)
 */
double &VectorContent::operator[](size_t i)
{
  ASSERT((i>=0) && (i<dimension),
      "VectorContent::operator[]() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(dimension)+"]");
  return *gsl_vector_ptr (content, i);
}

/** Constant accessor for (value of) component (i).
 * \param i component number
 * \return const component (i)
 */
const double VectorContent::operator[](size_t i) const
{
  ASSERT((i>=0) && (i<dimension),
      "VectorContent::operator[]() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(dimension)+"]");
  return gsl_vector_get(content, i);
}


/** Compares VectorContent \a to VectorContent \a b component-wise.
 * \param a base VectorContent
 * \param b VectorContent components to add
 * \return a == b
 */
bool VectorContent::operator==(const VectorContent& b) const
{
  bool status = true;
  ASSERT(dimension == b.dimension, "Dimenions of VectorContents to compare differ");
  for (size_t i=0;i<dimension;i++)
    status = status && (fabs(at(i) - b.at(i)) <= LINALG_MYEPSILON);
  return status;
};

/** Sums VectorContent \a to this lhs component-wise.
 * \param a base VectorContent
 * \param b VectorContent components to add
 * \return lhs + a
 */
const VectorContent& VectorContent::operator+=(const VectorContent& b)
{
  ASSERT(dimension == b.dimension, "Dimenions of VectorContents to compare differ");
  for (size_t i=0;i<dimension;i++)
    at(i) = at(i)+b.at(i);
  return *this;
};

/** Subtracts VectorContent \a from this lhs component-wise.
 * \param a base VectorContent
 * \param b VectorContent components to add
 * \return lhs - a
 */
const VectorContent& VectorContent::operator-=(const VectorContent& b)
{
  ASSERT(dimension == b.dimension, "Dimenions of VectorContents to compare differ");
  for (size_t i=0;i<dimension;i++)
    at(i) = at(i)-b.at(i);
  return *this;
};

/** factor each component of \a a times a double \a m.
 * \param a base VectorContent
 * \param m factor
 * \return lhs.Get(i) * m
 */
const VectorContent& VectorContent::operator*=(const double m)
{
  for (size_t i=0;i<dimension;i++)
    at(i) = at(i)*m;
  return *this;
};

/** Sums two VectorContents \a  and \b component-wise.
 * \param a first VectorContent
 * \param b second VectorContent
 * \return a + b
 */
VectorContent const operator+(const VectorContent& a, const VectorContent& b)
{
  ASSERT(a.dimension == b.dimension, "VectorContent::operator+() - dimensions have to match: "+toString(a.dimension)+" != "+toString(b.dimension)+"!");
  VectorContent x(a);
  for (size_t i=0;i<a.dimension;i++)
    x.at(i) = a.at(i)+b.at(i);
  return x;
};

/** Subtracts VectorContent \a from \b component-wise.
 * \param a first VectorContent
 * \param b second VectorContent
 * \return a - b
 */
VectorContent const operator-(const VectorContent& a, const VectorContent& b)
{
  ASSERT(a.dimension == b.dimension, "VectorContent::operator-() - dimensions have to match: "+toString(a.dimension)+" != "+toString(b.dimension)+"!");
  VectorContent x(a);
  for (size_t i=0;i<a.dimension;i++)
    x.at(i) = a.at(i)-b.at(i);
  return x;
};

/** Factors given VectorContent \a a times \a m.
 * \param a VectorContent
 * \param m factor
 * \return m * a
 */
VectorContent const operator*(const VectorContent& a, const double m)
{
  VectorContent x(a);
  for (size_t i=0;i<a.dimension;i++)
    x.at(i) = a.at(i)*m;
  return x;
};

/** Factors given VectorContent \a a times \a m.
 * \param m factor
 * \param a VectorContent
 * \return m * a
 */
VectorContent const operator*(const double m, const VectorContent& a )
{
  VectorContent x(a);
  for (size_t i=0;i<a.dimension;i++)
    x.at(i) = a.at(i)*m;
  return x;
};

ostream& operator<<(ostream& ost, const VectorContent& m)
{
  ost << "[";
  for (size_t i=0;i<m.dimension;i++) {
    ost << m.at(i);
    if (i != m.dimension-1)
      ost << " ";
  }
  ost << "]";
  return ost;
};

/* ====================== Checking State ============================ */
/** Checks whether vector has all components zero.
 * TODO: This might see some numerical instabilities for huge dimension and small number.
 * For stability one should sort the order of summing!
 * @return true - vector is zero, false - vector is not
 */
bool VectorContent::IsZero() const
{
  double result = 0.;
  for (size_t i = dimension; i--; )
    result += fabs(at(i));
  return (result <= LINALG_MYEPSILON);
};

/** Checks whether vector has length of 1.
 * @return true - vector is normalized, false - vector is not
 */
bool VectorContent::IsOne() const
{
  double NormValue = 0.;
  for (size_t i=dimension;--i;)
    NormValue += at(i)*at(i);
  return (fabs(NormValue - 1.) <= LINALG_MYEPSILON);
};

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

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

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

VectorContent VectorContent::getNormalized() const{
  VectorContent res= *this;
  res.Normalize();
  return res;
}

/* ======================== Properties ============================= */
/** Calculates the squared distance to some other vector.
 * @param y other vector
 * @return \f$|(\text{*this})-y|^2\f$
 */
double VectorContent::DistanceSquared(const VectorContent &y) const
{
  double res = 0.;
  for (int i=dimension;i--;)
    res += (at(i)-y[i])*(at(i)-y[i]);
  return (res);
}

/** Calculates scalar product between \a *this and \a b.
 * @param b other vector
 * @return \f$| (*this) \cdot (b)|^2\f$
 */
double VectorContent::ScalarProduct(const VectorContent &y) const
{
  return ((*this)*y);
}

/** Calculates the angle between \a *this and \a y.
 *
 * @param y other vector
 * @return \f$\acos\bigl(frac{\langle \text{*this}, y \rangle}{|\text{*this}||y|}\bigr)\f$
 */
double VectorContent::Angle(const VectorContent &y) const
{
  double norm1 = Norm(), norm2 = y.Norm();
  double angle = -1;
  if ((fabs(norm1) > LINALG_MYEPSILON) && (fabs(norm2) > LINALG_MYEPSILON))
    angle = this->ScalarProduct(y)/norm1/norm2;
  // -1-LINALG_MYEPSILON occured due to numerical imprecision, catch ...
  //Log() << Verbose(2) << "INFO: acos(-1) = " << acos(-1) << ", acos(-1+LINALG_MYEPSILON) = " << acos(-1+LINALG_MYEPSILON) << ", acos(-1-LINALG_MYEPSILON) = " << acos(-1-LINALG_MYEPSILON) << "." << endl;
  if (angle < -1)
    angle = -1;
  if (angle > 1)
    angle = 1;
  return acos(angle);
}

/* ======================== Properties ============================= */
/** Calculates scalar product between \a *this and \a b.
 * @param b other vector
 * @return \f$| (*this) \cdot (b)|^2\f$
 */
const double VectorContent::operator*(const VectorContent& b) const
{
  double res = 0.;
  gsl_blas_ddot(content, b.content, &res);
  return res;
};

/** Getter for VectorContent::dimension.
 * @return VectorContent::dimension
 */
const size_t VectorContent::getDimension() const
{
  return dimension;
}