source: trunk/yat/utility/matrix.cc @ 1027

// $Id: matrix.cc 1027 2008-02-02 21:29:29Z peter $

/*
  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
  Copyright (C) 2005, 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
  Copyright (C) 2007 Jari Häkkinen, Peter Johansson

  This file is part of the yat library, http://trac.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 2 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 this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "matrix.h"
#include "vector.h"
#include "VectorBase.h"
#include "VectorConstView.h"
#include "VectorView.h"
#include "utility.h"

#include <cassert>
#include <cmath>
#include <sstream>
#include <vector>

#include <gsl/gsl_blas.h>

namespace theplu {
namespace yat {
namespace utility {


  matrix::matrix(void)
    : blas_result_(NULL), m_(NULL), view_(NULL), view_const_(NULL),
      proxy_m_(NULL)
  {
  }


  matrix::matrix(const size_t& r, const size_t& c, double init_value)
    : blas_result_(NULL), m_(gsl_matrix_alloc(r,c)), view_(NULL),
      view_const_(NULL), proxy_m_(m_)
  {
    if (!m_)
      throw utility::GSL_error("matrix::matrix failed to allocate memory");
    all(init_value);
  }


  matrix::matrix(const matrix& o)
    : blas_result_(NULL), m_(o.create_gsl_matrix_copy()), view_(NULL),
      view_const_(NULL), proxy_m_(m_)
  {
  }


  matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
                 size_t n_row, size_t n_column)
    : blas_result_(NULL), view_const_(NULL)
  {
    view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
                                                     offset_row,offset_column,
                                                     n_row,n_column));
    if (!view_)
      throw utility::GSL_error("matrix::matrix failed to setup view");
    proxy_m_ = m_ = &(view_->matrix);
  }


  // Constructor that gets data from istream
  matrix::matrix(std::istream& is, char sep) 
    throw (utility::IO_error,std::exception)
    : blas_result_(NULL), view_(NULL), view_const_(NULL)
  {
    // read the data file and store in stl vectors (dynamically
    // expandable)
    std::vector<std::vector<double> > data_matrix;
    u_int nof_columns=0;
    u_int nof_rows = 0;
    std::string line;
    while(getline(is, line, '\n')){
      // Ignoring empty lines
      if (!line.size()) {
        continue;
      }
      nof_rows++;
      std::vector<double> v;
      std::string element;
      std::stringstream ss(line);
      
      bool ok=true;
      while(ok) {
        if(sep=='\0')
          ok=(ss>>element);
        else
          ok=getline(ss, element, sep);
        if(!ok)
          break;
        
        if(utility::is_double(element)) {
          v.push_back(atof(element.c_str()));
        }
        else if (!element.size() || utility::is_nan(element)) {
          v.push_back(std::numeric_limits<double>::quiet_NaN());
        }
        else {
          std::stringstream ss("Warning: '");
          ss << element << "' is not accepted as a matrix element.";
          throw IO_error(ss.str());
        }
      }           
      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
          v.push_back(std::numeric_limits<double>::quiet_NaN());
      if (!nof_columns)
        nof_columns=v.size();
      else if (v.size()!=nof_columns) {
        std::ostringstream s;
        s << "matrix::matrix(std::istream&, char) data file error: "
          << "line " << nof_rows << " has " << v.size()
          << " columns; expected " << nof_columns << " columns.";
        throw utility::IO_error(s.str());
      }
      data_matrix.push_back(v);
    }

    // manipulate the state of the stream to be good
    is.clear(std::ios::goodbit);
    // convert the data to a gsl matrix
    proxy_m_ = m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
    if (!m_)
      throw utility::GSL_error("matrix::matrix failed to allocate memory");

    // if gsl error handler disabled, out of bounds index will not
    // abort the program.
    for(u_int i=0;i<nof_rows;i++)
      for(u_int j=0;j<nof_columns;j++)
        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
  }


  matrix::~matrix(void)
  {
    delete_allocated_memory();
    if (blas_result_)
      gsl_matrix_free(blas_result_);
  }


  const matrix& matrix::clone(const matrix& other)
  {
    if (this!=&other) {

      delete_allocated_memory();

      if (other.view_) {
        view_ = new gsl_matrix_view(*other.view_);
        proxy_m_ = m_ = &(view_->matrix);
      }
      else if (other.view_const_) {
        view_const_ = new gsl_matrix_const_view(*other.view_const_);
        proxy_m_ = &(view_const_->matrix);
      }
      else if (other.m_)
        proxy_m_ = m_ = other.create_gsl_matrix_copy();

      // no need to delete blas_result_ if the number of rows fit, it
      // may be useful later.
      if (blas_result_ && (blas_result_->size1!=rows())) {
        gsl_matrix_free(blas_result_);
        blas_result_=NULL;
      }
    }
    return *this;
  } 


  size_t matrix::columns(void) const
  {
    if (!proxy_m_)
      return 0;
    return proxy_m_->size2;
  }


  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
  {
    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
    if (!m)
      throw utility::GSL_error("matrix::create_gsl_matrix_copy failed to allocate memory");
    if (gsl_matrix_memcpy(m,proxy_m_))
      throw utility::GSL_error("matrix::create_gsl_matrix_copy dimension mis-match");
    return m;
  }


  void matrix::delete_allocated_memory(void)
  {
    if (view_)
      delete view_;
    else if (view_const_)
      delete view_const_;
    else if (m_)
      gsl_matrix_free(m_);
    blas_result_=NULL;
    proxy_m_=m_=NULL;
  }


  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(std::string("matrix::div_elements",status));
  }


  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++)
      for (size_t j=0; j<columns(); j++)
        // The two last condition checks are needed for NaN detection
        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
          return false;
    return true;
  }


  const gsl_matrix* matrix::gsl_matrix_p(void) const
  {
    return proxy_m_;
  }


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


  bool matrix::isview(void) const
  {
    return view_ || view_const_;
  }


  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(std::string("matrix::mul_elements",status));
  }


  void matrix::resize(size_t r, size_t c, double init_value)
  {
    delete_allocated_memory();

    proxy_m_ = m_ = gsl_matrix_alloc(r,c);
    if (!m_)
      throw utility::GSL_error("matrix::matrix failed to allocate memory");
    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())) {
      gsl_matrix_free(blas_result_);
      blas_result_=NULL;
    }
  }


  size_t matrix::rows(void) const
  {
    if (!proxy_m_)
      return 0;
    return proxy_m_->size1;
  }


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


  VectorView matrix::column_vec(size_t col)
  {
    VectorView res(*this, col, false);
    return res;
  }


  const VectorConstView matrix::column_const_vec(size_t col) const
  {
    return VectorConstView(*this, col, false);
  }


  const VectorConstView matrix::row_const_vec(size_t col) const
  {
    return VectorConstView(*this, col, true);
  }


  VectorView matrix::row_vec(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(std::string("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(std::string("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(std::string("matrix::swap_rows",status));
  }


  void matrix::transpose(void)
  {
    assert(m_);
    if (columns()==rows())
      gsl_matrix_transpose(m_); // this never fails
    else {
      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
      if (!transposed)
        throw utility::GSL_error("matrix::transpose failed to allocate memory");
      // next line never fails if allocation above succeeded.
      gsl_matrix_transpose_memcpy(transposed,m_);
      gsl_matrix_free(m_);
      proxy_m_ = m_ = transposed;
      if (blas_result_) {
        gsl_matrix_free(blas_result_);
        blas_result_=NULL;
      }
    }
  }


  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(proxy_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 )
  {
    assert(m_);
    if (this!=&other)
      if (gsl_matrix_memcpy(m_, other.gsl_matrix_p()))
        throw utility::GSL_error("matrix::create_gsl_matrix_copy dimension mis-match");
    return *this;
  }


  const matrix& matrix::operator+=(const matrix& other)
  {
    assert(m_);
    int status=gsl_matrix_add(m_, other.proxy_m_);
    if (status)
      throw utility::GSL_error(std::string("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.proxy_m_);
    if (status)
      throw utility::GSL_error(std::string("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_);
    if ( blas_result_ && ((blas_result_->size1!=rows()) ||
                          (blas_result_->size2!=other.columns())) ) {
      gsl_matrix_free(blas_result_);
      blas_result_=NULL;
    }
    if (!blas_result_) {
      blas_result_ = gsl_matrix_alloc(rows(),other.columns());
      if (!blas_result_)
        throw utility::GSL_error("matrix::operator*= failed to allocate memory");
    }
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.proxy_m_, 0.0,
                   blas_result_);
    gsl_matrix* tmp=m_;
    proxy_m_ = m_ = blas_result_;
    blas_result_=tmp;
    return *this;
  }


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


  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)
  {
    size_t rows=templat.rows();
    size_t columns=templat.columns();
    if (rows!=flag.rows() && columns!=flag.columns())
      flag.clone(matrix(rows,columns,1.0));
    else
      flag.all(1.0);
    bool nan=false;
    for (size_t i=0; i<rows; i++)
      for (size_t j=0; j<columns; j++) 
        if (std::isnan(templat(i,j))) {
          flag(i,j)=0;
          nan=true;
        }
    return nan;
  }


  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(std::string("swap(matrix&,matrix&)",status));
  }


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


  vector operator*(const matrix& m, const VectorBase& v)
  {
    utility::vector res(m.rows());
    for (size_t i=0; i<res.size(); ++i)
      res(i) = VectorConstView(m,i) * v;
    return res;
  }


  vector operator*(const VectorBase& v, const matrix& m)
  {
    utility::vector res(m.columns());
    for (size_t i=0; i<res.size(); ++i)
      res(i) = v * VectorConstView(m,i,false);
    return res;
  }

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