/*
 * gslvector.cpp
 *
 *  Created on: Jan 8, 2010
 *      Author: heber
 */

#include <cassert>
#include <cmath>

#include "gslvector.hpp"
#include "defs.hpp"

/** Constructor of class GSLVector.
 * Allocates GSL structures
 * \param m dimension of vector
 */
GSLVector::GSLVector(size_t m) : dimension(m)
{
  vector = gsl_vector_calloc(dimension);
};

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

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

/** Destructor of class GSLVector.
 * Frees GSL structures
 */
GSLVector::~GSLVector()
{
  gsl_vector_free(vector);
};

/* ============================ Accessing =============================== */
/** 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 GSLVector::SetFromDoubleArray(double * x)
{
  gsl_vector_view m = gsl_vector_view_array (x, dimension);
  gsl_vector_memcpy (vector, &m.vector);
};

/** This function returns the i-th element of a vector.
 * If \a m lies outside the allowed range of 0 to GSLVector::dimension-1 then the error handler is invoked and 0 is returned.
 * \param m m-th element
 * \return  m-th element of vector
 */
double GSLVector::Get(size_t m) const
{
  return gsl_vector_get (vector, m);
};

/** This function sets the value of the \a m -th element of a vector to \a x.
 *  If \a m lies outside the allowed range of 0 to GSLVector::dimension-1 then the error handler is invoked.
 * \param m-th element to set
 * \param x value to set to
 */
void GSLVector::Set(size_t m, double x)
{
  gsl_vector_set (vector, m, x);
};

/** 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 GSLVector::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 *GSLVector::Pointer(size_t m) const
{
  return gsl_vector_ptr (vector, 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 GSLVector::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 *GSLVector::const_Pointer(size_t m) const
{
  return gsl_vector_const_ptr (vector, m);
};

/** Returns the dimension of the vector.
 * \return dimension of vector
 */
#ifdef HAVE_INLINE
inline
#endif
size_t GSLVector::GetDimension() const
{
  return dimension;
};

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

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

/** 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 GSLVector::SetBasis(size_t i)
{
  return gsl_vector_set_basis (vector, 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 GSLVector::SwapElements(size_t i, size_t j)
{
  return gsl_vector_swap_elements (vector, i, j);
};

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


/* ========================== Operators =============================== */
/** Compares GSLVector \a to GSLVector \a b component-wise.
 * \param a base GSLVector
 * \param b GSLVector components to add
 * \return a == b
 */
bool operator==(const GSLVector& a, const GSLVector& b)
{
  bool status = true;
  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
  for (size_t i=0;i<a.GetDimension();i++)
    status = status && (fabs(a.Get(i) - b.Get(i)) < MYEPSILON);
  return status;
};

/** Sums GSLVector \a to this lhs component-wise.
 * \param a base GSLVector
 * \param b GSLVector components to add
 * \return lhs + a
 */
const GSLVector& operator+=(GSLVector& a, const GSLVector& b)
{
  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
  for (size_t i=0;i<a.GetDimension();i++)
    a.Set(i,a.Get(i)+b.Get(i));
  return a;
};

/** Subtracts GSLVector \a from this lhs component-wise.
 * \param a base GSLVector
 * \param b GSLVector components to add
 * \return lhs - a
 */
const GSLVector& operator-=(GSLVector& a, const GSLVector& b)
{
  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
  for (size_t i=0;i<a.GetDimension();i++)
    a.Set(i,a.Get(i)-b.Get(i));
  return a;
};

/** factor each component of \a a times a double \a m.
 * \param a base GSLVector
 * \param m factor
 * \return lhs.Get(i) * m
 */
const GSLVector& operator*=(GSLVector& a, const double m)
{
  for (size_t i=0;i<a.GetDimension();i++)
    a.Set(i,a.Get(i)*m);
  return a;
};

/** Sums two GSLVectors \a  and \b component-wise.
 * \param a first GSLVector
 * \param b second GSLVector
 * \return a + b
 */
GSLVector const operator+(const GSLVector& a, const GSLVector& b)
{
  GSLVector x(a);
  for (size_t i=0;i<a.GetDimension();i++)
    x.Set(i,a.Get(i)+b.Get(i));
  return x;
};

/** Subtracts GSLVector \a from \b component-wise.
 * \param a first GSLVector
 * \param b second GSLVector
 * \return a - b
 */
GSLVector const operator-(const GSLVector& a, const GSLVector& b)
{
  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
  GSLVector x(a);
  for (size_t i=0;i<a.GetDimension();i++)
    x.Set(i,a.Get(i)-b.Get(i));
  return x;
};

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

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

ostream& operator<<(ostream& ost, const GSLVector& m)
{
  ost << "(";
  for (size_t i=0;i<m.GetDimension();i++) {
    ost << m.Get(i);
    if (i != 2)
      ost << ",";
  }
  ost << ")";
  return ost;
};

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

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