source: trunk/yat/utility/Matrix.cc @ 3743

// $Id: Matrix.cc 3743 2018-07-12 00:43:25Z peter $

/*
  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
  Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009, 2010, 2011, 2012, 2017, 2018 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "Matrix.h"

#include "DiagonalMatrix.h"
#include "Exception.h"
#include "SVD.h"
#include "Vector.h"
#include "VectorBase.h"
#include "VectorConstView.h"
#include "VectorView.h"
#include "utility.h"

#include <gsl/gsl_blas.h>

#include <algorithm>
#include <cassert>
#include <climits>
#include <cmath>
#include <iostream>
#include <limits>
#include <sstream>
#include <vector>

namespace theplu {
namespace yat {
namespace utility {


  Matrix::Matrix(void)
    : blas_result_(NULL), m_(NULL)
  {
  }


  Matrix::Matrix(const size_t& r, const size_t& c, double init_value)
    : blas_result_(NULL), m_(NULL)
  {
    resize(r,c,init_value);
  }


  Matrix::Matrix(const Matrix& o)
    : blas_result_(NULL), m_(o.create_gsl_matrix_copy())
  {
  }


#ifdef YAT_HAVE_RVALUE
  Matrix::Matrix(Matrix&& other) noexcept
    : blas_result_(NULL), m_(NULL)
  {
    std::swap(m_, other.m_);
  }
#endif


  // Constructor that gets data from istream
  Matrix::Matrix(std::istream& is, char sep)
    throw (utility::IO_error,std::exception)
    : blas_result_(NULL)
  {
    if (!is.good())
      throw utility::IO_error("Matrix: istream is not good");

    // read the data file and store in stl vectors (dynamically
    // expandable)
    std::vector<std::vector<double> > data_matrix;
    try {
      load(is, data_matrix, sep, '\n', true);
    }
    catch (utility::IO_error& e) {
      std::stringstream ss(e.what());
      ss << "\nMatrix(std::istream&): invalid dimensions\n";
      throw IO_error(ss.str());
    }
    catch (runtime_error& e) {
      std::stringstream ss(e.what());
      ss << "\nMatrix(std::istream&): invalid matrix element\n";
      throw IO_error(ss.str());
    }

    unsigned int nof_rows = data_matrix.size();
    // if stream was empty, create nothing
    if (!nof_rows)
      return;

    unsigned int nof_columns=data_matrix[0].size();

    // convert the data to a gsl matrix
    m_ = detail::allocate_matrix(nof_rows, nof_columns);

    for (size_t i=0; i<data_matrix.size(); ++i) {
      assert(data_matrix[i].size()==columns());
      assert(i<rows());
      std::copy(data_matrix[i].begin(), data_matrix[i].end(), begin_row(i));
    }
  }


  Matrix::~Matrix(void)
  {
    delete_allocated_memory();
    detail::deallocate(blas_result_);
  }


  void Matrix::all(const double value)
  {
    assert(m_);
    gsl_matrix_set_all(m_, value);
  }


  Matrix::iterator Matrix::begin(void)
  {
    return iterator(&(*this)(0,0), 1);
  }


  Matrix::const_iterator Matrix::begin(void) const
  {
    return const_iterator(&(*this)(0,0), 1);
  }


  Matrix::column_iterator Matrix::begin_column(size_t i)
  {
    return column_iterator(&(*this)(0,i), this->columns());
  }


  Matrix::const_column_iterator Matrix::begin_column(size_t i) const
  {
    return const_column_iterator(&(*this)(0,i), this->columns());
  }


  Matrix::row_iterator Matrix::begin_row(size_t i)
  {
    return row_iterator(&(*this)(i,0), 1);
  }


  Matrix::const_row_iterator Matrix::begin_row(size_t i) const
  {
    return const_row_iterator(&(*this)(i,0), 1);
  }


  VectorView Matrix::column_view(size_t col)
  {
    VectorView res(*this, col, false);
    return res;
  }


  const VectorConstView Matrix::column_const_view(size_t col) const
  {
    return VectorConstView(*this, col, false);
  }


  size_t Matrix::columns(void) const
  {
    return detail::columns(m_);
  }


  gsl_matrix* Matrix::create_gsl_matrix_copy(void) const
  {
    return detail::create_copy(m_);
  }


  void Matrix::delete_allocated_memory(void)
  {
    detail::deallocate(m_);
  }


  void Matrix::div(const Matrix& other)
  {
    assert(m_);
    int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p());
    if (status)
      throw utility::GSL_error("matrix::div_elements",status);
  }


  Matrix::iterator Matrix::end(void)
  {
    return iterator(&(*this)(0,0)+rows()*columns(), 1);
  }


  Matrix::const_iterator Matrix::end(void) const
  {
    return const_iterator(&(*this)(0,0)+rows()*columns(), 1);
  }


  Matrix::column_iterator Matrix::end_column(size_t i)
  {
    return column_iterator(&(*this)(0,i)+rows()*columns(), this->columns());
  }


  Matrix::const_column_iterator Matrix::end_column(size_t i) const
  {
    return const_column_iterator(&(*this)(0,i)+rows()*columns(),this->columns());
  }


  Matrix::row_iterator Matrix::end_row(size_t i)
  {
    return row_iterator(&(*this)(i,0)+columns(), 1);
  }


  Matrix::const_row_iterator Matrix::end_row(size_t i) const
  {
    return const_row_iterator(&(*this)(i,0)+columns(), 1);
  }


  bool Matrix::equal(const Matrix& other, const double d) const
  {
    if (this==&other)
      return true;
    if (columns()!=other.columns() || rows()!=other.rows())
      return false;
    for (size_t i=0; i<rows(); i++)
      if (!row_const_view(i).equal(other.row_const_view(i), d))
        return false;
    return true;
  }


  const gsl_matrix* Matrix::gsl_matrix_p(void) const
  {
    return m_;
  }


  gsl_matrix* Matrix::gsl_matrix_p(void)
  {
    return m_;
  }


  void Matrix::mul(const Matrix& other)
  {
    assert(m_);
    int status=gsl_matrix_mul_elements(m_, other.gsl_matrix_p());
    if (status)
      throw utility::GSL_error("Matrix::mul_elements",status);
  }


  void Matrix::resize(size_t r, size_t c, double init_value)
  {
    detail::reallocate(m_, r, c);

    if (m_)
      all(init_value);

    // no need to delete blas_result_ if the number of rows fit, it
    // may be useful later.
    if (blas_result_ && (blas_result_->size1!=rows())) {
      detail::deallocate(blas_result_);
    }
  }


  size_t Matrix::rows(void) const
  {
    return detail::rows(m_);
  }


  const VectorConstView Matrix::row_const_view(size_t col) const
  {
    return VectorConstView(*this, col, true);
  }


  VectorView Matrix::row_view(size_t row)
  {
    VectorView res(*this, row, true);
    return res;
  }


  void Matrix::swap_columns(const size_t i, const size_t j)
  {
    assert(m_);
    int status=gsl_matrix_swap_columns(m_, i, j);
    if (status)
      throw utility::GSL_error("Matrix::swap_columns",status);
  }


  void Matrix::swap_rowcol(const size_t i, const size_t j)
  {
    assert(m_);
    int status=gsl_matrix_swap_rowcol(m_, i, j);
    if (status)
      throw utility::GSL_error("Matrix::swap_rowcol",status);
  }


  void Matrix::swap_rows(const size_t i, const size_t j)
  {
    assert(m_);
    int status=gsl_matrix_swap_rows(m_, i, j);
    if (status)
      throw utility::GSL_error("Matrix::swap_rows",status);
  }


  void Matrix::transpose(void)
  {
    assert(m_);
    if (columns()==rows())
      gsl_matrix_transpose(m_); // this never fails
    else {
      gsl_matrix* transposed = detail::allocate_matrix(columns(), rows());
      // next line never fails if allocation above succeeded.
      gsl_matrix_transpose_memcpy(transposed,m_);
      gsl_matrix_free(m_);
      m_ = transposed;
      detail::deallocate(blas_result_);
    }
  }


  double& Matrix::operator()(size_t row, size_t column)
  {
    assert(m_);
    assert(row<rows());
    assert(column<columns());
    double* d=gsl_matrix_ptr(m_, row, column);
    if (!d)
      throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
    return *d;
  }


  const double& Matrix::operator()(size_t row, size_t column) const
  {
    assert(row<rows());
    assert(column<columns());
    const double* d=gsl_matrix_const_ptr(m_, row, column);
    if (!d)
      throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
    return *d;
  }


  bool Matrix::operator==(const Matrix& other) const
  {
    return equal(other);
  }


  bool Matrix::operator!=(const Matrix& other) const
  {
    return !equal(other);
  }


  const Matrix& Matrix::operator=( const Matrix& other )
  {
    if (this!=&other) {
      detail::copy(m_, other.gsl_matrix_p());
    }
    return *this;
  }


#ifdef YAT_HAVE_RVALUE
  Matrix& Matrix::operator=(Matrix&& other)
  {
    std::swap(m_, other.m_);
    return *this;
  }
#endif


  const Matrix& Matrix::operator+=(const Matrix& other)
  {
    assert(m_);
    int status=gsl_matrix_add(m_, other.m_);
    if (status)
      throw utility::GSL_error("Matrix::operator+=", status);
    return *this;
  }


  const Matrix& Matrix::operator+=(const double d)
  {
    assert(m_);
    gsl_matrix_add_constant(m_, d);
    return *this;
  }


  const Matrix& Matrix::operator-=(const Matrix& other)
  {
    assert(m_);
    int status=gsl_matrix_sub(m_, other.m_);
    if (status)
      throw utility::GSL_error("Matrix::operator-=", status);
    return *this;
  }


  const Matrix& Matrix::operator-=(const double d)
  {
    assert(m_);
    gsl_matrix_add_constant(m_, -d);
    return *this;
  }


  const Matrix& Matrix::operator*=(const Matrix& other)
  {
    assert(m_);
    assert(other.columns());
    detail::reallocate(blas_result_, rows(), other.columns());
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0,
                   blas_result_);
    std::swap(m_, blas_result_);
    return *this;
  }


  const Matrix& Matrix::operator*=(const double d)
  {
    assert(m_);
    gsl_matrix_scale(m_, d);
    return *this;
  }


  void inverse_svd(const Matrix& A, Matrix& result)
  {
    assert(A.rows() == A.columns());
    SVD svd(A);
    // A = U S V'
    svd.decompose();

    // result = A^-1 = V * S^-1 * U'

    DiagonalMatrix D(A.rows(), A.columns());
    // If eigenvalue is zero, A is not invertable, perhaps that should
    // trigger an exception
    for (size_t i=0; i<D.rows(); ++i)
      D(i) = 1.0 / svd.s()(i);

    result = svd.V() * D * transpose(svd.U());

    assert(result.columns() == result.rows());
    assert(result.rows() == A.rows());
  }


  bool isnull(const Matrix& other)
  {
    return gsl_matrix_isnull(other.gsl_matrix_p());
  }


  double max(const Matrix& other)
  {
    return gsl_matrix_max(other.gsl_matrix_p());
  }


  double min(const Matrix& other)
  {
    return gsl_matrix_min(other.gsl_matrix_p());
  }


  void minmax_index(const Matrix& other,
                    std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max)
  {
    gsl_matrix_minmax_index(other.gsl_matrix_p(), &min.first, &min.second,
                            &max.first, &max.second);
  }


  bool nan(const Matrix& templat, Matrix& flag)
  {
    if (templat.rows()!=flag.rows() || templat.columns()!=flag.columns())
      flag.resize(templat.rows(),templat.columns(),1.0);
    return binary_weight(templat.begin(), templat.end(), flag.begin());
  }


  void swap(Matrix& a, Matrix& b)
  {
    assert(a.gsl_matrix_p());
    assert(b.gsl_matrix_p());
    int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());
    if (status)
      throw utility::GSL_error("swap(Matrix&,Matrix&)",status);
  }


  std::ostream& operator<<(std::ostream& s, const Matrix& m)
  {
    s.setf(std::ios::dec);
    std::streamsize precision = s.precision();
    s.precision(std::numeric_limits<double>().digits10);
    for(size_t i=0, j=0; i<m.rows(); i++) {
      s << m(i,0);
      for (j=1; j<m.columns(); ++j)
        s << s.fill() << m(i,j);
      if (i<m.rows()-1)
        s << "\n";
    }
    s.precision(precision);
    return s;
  }


  double trace(const Matrix& m)
  {
    size_t n = std::min(m.rows(), m.columns());
    double sum = 0.0;
    for (size_t i=0; i<n; ++i)
      sum += m(i,i);
    return sum;
  }


}}} // of namespace utility, yat and thep