/*
 * MatrixContent.cpp
 *
 *  Created on: Nov 14, 2010
 *      Author: heber
 */


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

#include "Helpers/MemDebug.hpp"

#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "Helpers/Assert.hpp"
#include "Exceptions/NotInvertibleException.hpp"
#include "Helpers/fast_functions.hpp"
#include "Helpers/Assert.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/VectorContent.hpp"
#include "LinearAlgebra/MatrixContent.hpp"

#include <gsl/gsl_blas.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_vector.h>
#include <cmath>
#include <iostream>

using namespace std;


/** Constructor for class MatrixContent.
 * \param rows number of rows
 * \param columns number of columns
 */
MatrixContent::MatrixContent(size_t _rows, size_t _columns) :
    rows(_rows),
    columns(_columns)
{
  content = gsl_matrix_calloc(rows, columns);
}

/** Constructor for class MatrixContent.
 * \param rows number of rows
 * \param columns number of columns
 * \param *src array with components to initialize matrix with
 */
MatrixContent::MatrixContent(size_t _rows, size_t _columns, const double *src) :
        rows(_rows),
        columns(_columns)
{
  content = gsl_matrix_calloc(rows, columns);
  set(0,0, src[0]);
  set(1,0, src[1]);
  set(2,0, src[2]);

  set(0,1, src[3]);
  set(1,1, src[4]);
  set(2,1, src[5]);

  set(0,2, src[6]);
  set(1,2, src[7]);
  set(2,2, src[8]);
}

/** Constructor for class MatrixContent.
 * We embed the given gls_matrix pointer within this class and set it to NULL
 * afterwards.
 * \param *src source gsl_matrix vector to embed within this class
 */
MatrixContent::MatrixContent(gsl_matrix *&src) :
        rows(src->size1),
        columns(src->size2)
{
  content = gsl_matrix_alloc(src->size1, src->size2);
  gsl_matrix_memcpy(content,src);
//  content = src;
//  src = NULL;
}

/** Copy constructor for class MatrixContent.
 * \param &src reference to source MatrixContent
 */
MatrixContent::MatrixContent(const MatrixContent &src) :
        rows(src.rows),
        columns(src.columns)
{
  content = gsl_matrix_alloc(src.rows, src.columns);
  gsl_matrix_memcpy(content,src.content);
}

/** Copy constructor for class MatrixContent.
 * \param *src pointer to source MatrixContent
 */
MatrixContent::MatrixContent(const MatrixContent *src) :
        rows(src->rows),
        columns(src->columns)
{
  ASSERT(src != NULL, "MatrixContent::MatrixContent - pointer to source matrix is NULL!");
  content = gsl_matrix_alloc(src->rows, src->columns);
  gsl_matrix_memcpy(content,src->content);
}

/** Destructor for class MatrixContent.
 */
MatrixContent::~MatrixContent()
{
  gsl_matrix_free(content);
}

/** Set matrix to identity.
 */
void MatrixContent::setIdentity()
{
  for(int i=rows;i--;){
    for(int j=columns;j--;){
      set(i,j,i==j);
    }
  }
}

/** Set all matrix components to zero.
 */
void MatrixContent::setZero()
{
  for(int i=rows;i--;){
    for(int j=columns;j--;){
      set(i,j,0.);
    }
  }
}

/** Copy operator for MatrixContent with self-assignment check.
 * \param &src matrix to compare to
 * \return reference to this
 */
MatrixContent &MatrixContent::operator=(const MatrixContent &src)
{
  if(&src!=this){
    gsl_matrix_memcpy(content,src.content);
  }
  return *this;
}

/** Addition operator.
 * \param &rhs matrix to add
 * \return reference to this
 */
const MatrixContent &MatrixContent::operator+=(const MatrixContent &rhs)
{
  gsl_matrix_add(content, rhs.content);
  return *this;
}

/** Subtraction operator.
 * \param &rhs matrix to subtract
 * \return reference to this
 */
const MatrixContent &MatrixContent::operator-=(const MatrixContent &rhs)
    {
  gsl_matrix_sub(content, rhs.content);
  return *this;
}

/** Multiplication operator.
 * Note that here matrix have to have same dimensions.
 * \param &rhs matrix to multiply with
 * \return reference to this
 */
const MatrixContent &MatrixContent::operator*=(const MatrixContent &rhs)
{
  ASSERT(rows == rhs.rows,
      "MatrixContent::operator*=() - row dimension differ: "+toString(rows)+" != "+toString(rhs.rows)+".");
  ASSERT(columns == rhs.columns,
      "MatrixContent::operator*=() - columns dimension differ: "+toString(columns)+" != "+toString(rhs.columns)+".");
  (*this) = (*this)*rhs;
  return *this;
}

/** Multiplication with copy operator.
 * \param &rhs matrix to multiply with
 * \return reference to newly allocated MatrixContent
 */
const MatrixContent MatrixContent::operator*(const MatrixContent &rhs) const
{
  gsl_matrix *res = gsl_matrix_alloc(rows, columns);
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, content, rhs.content, 0.0, res);
  // gsl_matrix is taken over by constructor, hence no free
  MatrixContent tmp(res);
  gsl_matrix_free(res);
  return tmp;
}

/** Accessor for manipulating component (i,j).
 * \param i row number
 * \param j column number
 * \return reference to component (i,j)
 */
double &MatrixContent::at(size_t i, size_t j)
{
  ASSERT((i>=0) && (i<rows),
      "MatrixContent::at() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
  ASSERT((j>=0) && (j<columns),
      "MatrixContent::at() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
  return *gsl_matrix_ptr (content, i, j);
}

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

/** Setter for component (i,j).
 * \param i row numbr
 * \param j column numnber
 * \param value value to set componnt (i,j) to
 */
void MatrixContent::set(size_t i, size_t j, const double value)
{
  ASSERT((i>=0) && (i<rows),
      "MatrixContent::set() - Index i="+toString(i)+" for Matrix access out of range [0,"+toString(rows)+"]");
  ASSERT((j>=0) && (j<columns),
      "MatrixContent::set() - Index j="+toString(j)+" for Matrix access out of range [0,"+toString(columns)+"]");
  gsl_matrix_set(content,i,j,value);
}

/** Return transposed matrix.
 * \return new matrix that is transposed of this.
 */
MatrixContent MatrixContent::transpose() const
{
  std::cout << "const MatrixContent::transpose()." << std::endl;
  gsl_matrix *res = gsl_matrix_alloc(rows, columns);
  gsl_matrix_transpose_memcpy(res, content);
  MatrixContent newContent(res);
  gsl_matrix_free(res);
  return newContent;
}

/** Turn this matrix into its transposed.
 * Note that this is only possible if rows == columns.
 */
MatrixContent &MatrixContent::transpose()
{
  std::cout << "MatrixContent::transpose()." << std::endl;
  ASSERT( rows == columns,
      "MatrixContent::transpose() - cannot transpose onto itself as matrix not square: "+toString(rows)+"!="+toString(columns)+"!");
  double tmp;
  for (size_t i=0;i<rows;i++)
    for (size_t j=i+1;j<rows;j++) {
      tmp = at(j,i);
      at(j,i) = at(i,j);
      at(i,j) = tmp;
    }
  return *this;
}

/** Transform the matrix to its eigenbasis ans return resulting eigenvalues.
 * \warn return vector has to be freed'd
 * \return gsl_vector pointer to vector of eigenvalues
 */
gsl_vector* MatrixContent::transformToEigenbasis()
{
  ASSERT (rows == columns,
      "MatrixContent::transformToEigenbasis() - only implemented for square matrices: "+toString(rows)+"!="+toString(columns)+"!");
  gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(rows);
  gsl_vector *eval = gsl_vector_alloc(rows);
  gsl_matrix *evec = gsl_matrix_alloc(rows, rows);
  gsl_eigen_symmv(content, eval, evec, T);
  gsl_eigen_symmv_free(T);
  gsl_matrix_memcpy(content, evec);
  gsl_matrix_free(evec);
  return eval;
}

/** Scalar multiplication operator.
 * \param factor factor to scale with
 */
const MatrixContent &MatrixContent::operator*=(const double factor)
{
  gsl_matrix_scale(content, factor);
  return *this;
}

/** Scalar multiplication and copy operator.
 * \param factor factor to scale with
 * \param &mat MatrixContent to scale
 * \return copied and scaled MatrixContent
 */
const MatrixContent operator*(const double factor,const MatrixContent& mat)
{
  MatrixContent tmp = mat;
  tmp*=factor;
  return tmp;
}

/** Scalar multiplication and copy operator (with operands exchanged).
 * \param &mat MatrixContent to scale
 * \param factor factor to scale with
 * \return copied and scaled MatrixContent
 */
const MatrixContent operator*(const MatrixContent &mat,const double factor)
{
  return factor*mat;
}

/** Equality operator.
 * Note that we use numerical sensible checking, i.e. with threshold MYEPSILON.
 * \param &rhs MatrixContent to checks against
 */
bool MatrixContent::operator==(const MatrixContent &rhs) const
    {
  if ((rows == rhs.rows) && (columns == rhs.columns)) {
    for(int i=rows;i--;){
      for(int j=columns;j--;){
        if(fabs(at(i,j)-rhs.at(i,j))>MYEPSILON){
          return false;
        }
      }
    }
    return true;
  }
  return false;
}

Vector operator*(const MatrixContent &mat,const Vector &vec)
{
  Vector res;
  gsl_blas_dgemv( CblasNoTrans, 1.0, mat.content, vec.content->content, 0.0, res.content->content);
  return res;
}