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

// $Id: matrix.cc 732 2007-01-06 16:11:22Z 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

  This file is part of the yat library, http://lev.thep.lu.se/trac/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 "utility.h"

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

#include <gsl/gsl_blas.h>

namespace theplu {
namespace yat {
namespace utility {


  matrix::matrix(void)
    : m_(NULL), view_(NULL)
  {
  }


  matrix::matrix(const size_t& r, const size_t& c, double init_value)
    : view_(NULL)
  {
    m_ = gsl_matrix_alloc(r,c);
    set_all(init_value);
  }


  matrix::matrix(const matrix& o)
    : view_(NULL)
  {
    m_ = o.create_gsl_matrix_copy();
  }


  matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
                 size_t n_row, size_t n_column)
  {
    // Jari, exception handling needed here. Failure in setting up a
    // proper gsl_matrix_view is communicated by NULL pointer in the
    // view structure (cf. GSL manual). How about GSL error state?
    view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
                                                     offset_row,offset_column,
                                                     n_row,n_column));
    m_ = &(view_->matrix);
  }


  // Constructor that gets data from istream
  matrix::matrix(std::istream& is, char sep) 
    throw (utility::IO_error,std::exception)
    : view_(NULL)
  {
    // Markus to Jari, somewhere we should check that quiet_NaNs are supported
    // std::numeric_limits<double>::has_quiet_NaN has to be true.
    // Also in vector

    // 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 {
          // Jari, this should be communicated with as an exception.
          std::cerr << "Warning: '" << element 
                    << "' is not accepted as a matrix element." << std::endl;
        }
      }           
      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
    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
    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)
  {
    if (view_)
      delete view_;
    else if (m_)
      gsl_matrix_free(m_);
    m_=NULL;
  }


  int matrix::add(const matrix& b)
  {
    return gsl_matrix_add(m_, b.m_);
  }


  int matrix::add_constant(const double d)
  {
    return gsl_matrix_add_constant(m_, d);
  }


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


  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
  {
    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
    gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
    return m;
  }


  int matrix::div_elements(const matrix& b)
  {
    return gsl_matrix_div_elements(m_, b.m_);
  }


  bool matrix::equal(const matrix& other, const double d) const
  {
    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 m_;
  }


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


  bool matrix::isnull(void) const
  {
    return gsl_matrix_isnull(m_);
  }


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


  double matrix::max(void) const
  {
    return gsl_matrix_max(m_);
  }


  double matrix::min(void) const
  {
    return gsl_matrix_min(m_);
  }


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


  bool matrix::nan(matrix &m) const
  {
    m=matrix(rows(),columns(),1.0);
    bool nan=false;
    for (size_t i=0; i<rows(); i++)
      for (size_t j=0; j<columns(); j++) 
        if (std::isnan(operator()(i,j))) {
          m(i,j)=0;
          nan=true;
        }
    return nan;
  }


  int matrix::mul_elements(const matrix& b)
  {
    return gsl_matrix_mul_elements(m_, b.m_);
  }


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


  int matrix::scale(const double d)
  {
    return gsl_matrix_scale(m_, d);
  }


  int matrix::set(const matrix& mat)
  {
    return gsl_matrix_memcpy(m_, mat.m_);
  }


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


  int matrix::set_column(const size_t column, const vector& vec)
  {
    return gsl_matrix_set_col(m_, column, vec.gsl_vector_p());
  }


  int matrix::set_row(const size_t row, const vector& vec)
  {
    return gsl_matrix_set_row(m_, row, vec.gsl_vector_p());
  }


  int matrix::sub(const matrix& b)
  {
    return gsl_matrix_sub(m_, b.m_);
  }


  int matrix::swap(matrix& other)
  {
    return gsl_matrix_swap(m_, other.m_);
  }


  int matrix::swap_columns(const size_t i, const size_t j)
  {
    return gsl_matrix_swap_columns(m_, i, j);
  }


  int matrix::swap_rowcol(const size_t i, const size_t j)
  {
    return gsl_matrix_swap_rowcol(m_, i, j);
  }


  int matrix::swap_rows(const size_t i, const size_t j)
  {
    return gsl_matrix_swap_rows(m_, i, j);
  }


  // Jari, checkout GSL transpose support in GSL manual 8.4.9
  void matrix::transpose(void)
  {
    if (columns()==rows())
      gsl_matrix_transpose(m_);
    else {
      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
      gsl_matrix_transpose_memcpy(transposed,m_);
      gsl_matrix_free(m_);
      m_=transposed;
    }
  }


  double& matrix::operator()(size_t row, size_t column)
  {
    return (*gsl_matrix_ptr(m_, row, column));
  }


  const double& matrix::operator()(size_t row, size_t column) const
  {
    return (*gsl_matrix_const_ptr(m_, row, column));
  }


  const vector matrix::operator*(const vector&) const
  {
    std::cerr << "Not implemented:" << std::endl
              << "   const vector matrix::operator*(const vector&) const"
              << std::endl;
    return vector(0);
  }


  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 ) {
      if (view_)
        delete view_;
      else if (m_)
        gsl_matrix_free(m_);
      m_ = other.create_gsl_matrix_copy();
    }
    return *this;
  }


  const matrix& matrix::operator+=(const matrix& m)
  {
    add(m);
    return *this;
  }


  const matrix& matrix::operator+=(const double d)
  {
    add_constant(d);
    return *this;
  }


  const matrix& matrix::operator-=(const matrix& m)
  {
    sub(m);
    return *this;
  }


  const matrix& matrix::operator*=(const matrix& other)
  {
    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
    gsl_matrix_free(m_);
    m_=result;
    return *this;
  }


  const matrix& matrix::operator*=(const double d)
  {
    scale(d);
    return *this;
  }


  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;
  }

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