Ignore:
Timestamp:
Jul 17, 2008, 12:36:23 AM (13 years ago)
Author:
Peter
Message:

working on #363

File:
1 copied

Legend:

Unmodified
Added
Removed
  • trunk/yat/utility/MatrixWeighted.cc

    r1360 r1380  
    2525*/
    2626
    27 #include "Matrix.h"
    28 #include "Vector.h"
    29 #include "VectorBase.h"
    30 #include "VectorConstView.h"
    31 #include "VectorView.h"
    32 #include "utility.h"
     27#include "MatrixWeighted.h"
    3328
    3429#include <cassert>
    35 #include <cmath>
    36 #include <sstream>
    3730#include <vector>
    38 
    39 #include <gsl/gsl_blas.h>
    4031
    4132namespace theplu {
     
    4435
    4536
    46   Matrix::Matrix(void)
    47     : blas_result_(NULL), m_(NULL)
     37  MatrixWeighted::MatrixWeighted(void)
     38    : columns_(0)
    4839  {
    4940  }
    5041
    5142
    52   Matrix::Matrix(const size_t& r, const size_t& c, double init_value)
    53     : blas_result_(NULL), m_(gsl_matrix_alloc(r,c))
    54   {
    55     if (!m_)
    56       throw utility::GSL_error("Matrix::Matrix failed to allocate memory");
    57     all(init_value);
    58   }
    59 
    60 
    61   Matrix::Matrix(const Matrix& o)
    62     : blas_result_(NULL), m_(o.create_gsl_matrix_copy())
    63   {
    64   }
    65 
    66 
    67   // Constructor that gets data from istream
    68   Matrix::Matrix(std::istream& is, char sep)
    69     throw (utility::IO_error,std::exception)
    70     : blas_result_(NULL)
    71   {
    72     if (!is.good())
    73       throw utility::IO_error("Matrix: istream is not good");
    74 
    75     // read the data file and store in stl vectors (dynamically
    76     // expandable)
    77     std::vector<std::vector<double> > data_matrix;
    78     unsigned int nof_columns=0;
    79     unsigned int nof_rows = 0;
    80     std::string line;
    81     while(getline(is, line, '\n')){
    82       // Ignoring empty lines
    83       if (!line.size()) {
    84         continue;
    85       }
    86       nof_rows++;
    87       std::vector<double> v;
    88       std::string element;
    89       std::stringstream ss(line);
    90      
    91       bool ok=true;
    92       while(ok) {
    93         if(sep=='\0')
    94           ok=(ss>>element);
    95         else
    96           ok=getline(ss, element, sep);
    97         if(!ok)
    98           break;
    99        
    100         if (!element.size())
    101           v.push_back(std::numeric_limits<double>::quiet_NaN());
    102         else {
    103           try {
    104             v.push_back(convert<double>(element));
    105           }
    106           catch (std::runtime_error& e) {
    107             std::stringstream ss(e.what());
    108             ss << "\nMatrix.cc: " << element
    109                << " is not accepted as a matrix element\n";
    110             throw IO_error(ss.str());
    111           }
    112         }
    113       }           
    114       if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
    115           v.push_back(std::numeric_limits<double>::quiet_NaN());
    116       if (!nof_columns)
    117         nof_columns=v.size();
    118       else if (v.size()!=nof_columns) {
    119         std::ostringstream s;
    120         s << "Matrix::Matrix(std::istream&, char) data file error: "
    121           << "line " << nof_rows << " has " << v.size()
    122           << " columns; expected " << nof_columns << " columns.";
    123         throw utility::IO_error(s.str());
    124       }
    125       data_matrix.push_back(v);
    126     }
    127 
    128     // manipulate the state of the stream to be good
    129     is.clear(std::ios::goodbit);
    130    
    131     // if stream was empty, create nothing
    132     if (!nof_columns || !nof_rows)
    133       return;
    134 
    135     // convert the data to a gsl matrix
    136     m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
    137     if (!m_)
    138       throw utility::GSL_error("Matrix::Matrix failed to allocate memory");
    139 
    140     // if gsl error handler disabled, out of bounds index will not
    141     // abort the program.
    142     for(size_t i=0;i<nof_rows;i++)
    143       for(size_t j=0;j<nof_columns;j++)
    144         gsl_matrix_set( m_, i, j, data_matrix[i][j] );
    145   }
    146 
    147 
    148   Matrix::~Matrix(void)
    149   {
    150     delete_allocated_memory();
    151     if (blas_result_)
    152       gsl_matrix_free(blas_result_);
    153     blas_result_=NULL;
    154   }
    155 
    156 
    157   void Matrix::all(const double value)
    158   {
    159     assert(m_);
    160     gsl_matrix_set_all(m_, value);
    161   }
    162 
    163 
    164   Matrix::iterator Matrix::begin(void)
    165   {
    166     return iterator(&(*this)(0,0), 1);
    167   }
    168 
    169 
    170   Matrix::const_iterator Matrix::begin(void) const
    171   {
    172     return const_iterator(&(*this)(0,0), 1);
    173   }
    174 
    175 
    176   Matrix::column_iterator Matrix::begin_column(size_t i)
    177   {
    178     return iterator(&(*this)(0,i), this->columns());
    179   }
    180 
    181 
    182   Matrix::const_column_iterator Matrix::begin_column(size_t i) const
    183   {
    184     return const_iterator(&(*this)(0,i), this->columns());
    185   }
    186 
    187 
    188   Matrix::row_iterator Matrix::begin_row(size_t i)
    189   {
    190     return iterator(&(*this)(i,0), 1);
    191   }
    192 
    193 
    194   Matrix::const_row_iterator Matrix::begin_row(size_t i) const
    195   {
    196     return const_iterator(&(*this)(i,0), 1);
    197   }
    198 
    199 
    200   VectorView Matrix::column_view(size_t col)
    201   {
    202     VectorView res(*this, col, false);
    203     return res;
    204   }
    205 
    206 
    207   const VectorConstView Matrix::column_const_view(size_t col) const
    208   {
    209     return VectorConstView(*this, col, false);
    210   }
    211 
    212 
    213   size_t Matrix::columns(void) const
    214   {
    215     return (m_ ? m_->size2 : 0);
    216   }
    217 
    218 
    219   gsl_matrix* Matrix::create_gsl_matrix_copy(void) const
    220   {
    221     gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
    222     if (!m)
    223       throw utility::GSL_error("Matrix::create_gsl_matrix_copy failed to allocate memory");
    224     if (gsl_matrix_memcpy(m,m_))
    225       throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mis-match");
    226     return m;
    227   }
    228 
    229 
    230   void Matrix::delete_allocated_memory(void)
    231   {
    232     if (m_)
    233       gsl_matrix_free(m_);
    234     m_=NULL;
    235   }
    236 
    237 
    238   void Matrix::div(const Matrix& other)
    239   {
    240     assert(m_);
    241     int status=gsl_matrix_div_elements(m_, other.gsl_matrix_p());
    242     if (status)
    243       throw utility::GSL_error(std::string("matrix::div_elements",status));
    244   }
    245 
    246 
    247   Matrix::iterator Matrix::end(void)
    248   {
    249     return iterator(&(*this)(0,0)+rows()*columns(), 1);
    250   }
    251 
    252 
    253   Matrix::const_iterator Matrix::end(void) const
    254   {
    255     return const_iterator(&(*this)(0,0)+rows()*columns(), 1);
    256   }
    257 
    258 
    259   Matrix::column_iterator Matrix::end_column(size_t i)
    260   {
    261     return column_iterator(&(*this)(0,i)+rows()*columns(), this->columns());
    262   }
    263 
    264 
    265   Matrix::const_column_iterator Matrix::end_column(size_t i) const
    266   {
    267     return const_column_iterator(&(*this)(0,i)+rows()*columns(),this->columns());
    268   }
    269 
    270 
    271   Matrix::row_iterator Matrix::end_row(size_t i)
    272   {
    273     return row_iterator(&(*this)(i,0)+columns(), 1);
    274   }
    275 
    276 
    277   Matrix::const_row_iterator Matrix::end_row(size_t i) const
    278   {
    279     return const_row_iterator(&(*this)(i,0)+columns(), 1);
    280   }
    281 
    282 
    283   bool Matrix::equal(const Matrix& other, const double d) const
    284   {
    285     if (this==&other)
    286       return true;
    287     if (columns()!=other.columns() || rows()!=other.rows())
    288       return false;
    289     for (size_t i=0; i<rows(); i++)
    290       for (size_t j=0; j<columns(); j++)
    291         // The two last condition checks are needed for NaN detection
    292         if (fabs( (*this)(i,j)-other(i,j) ) > d ||
    293             (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
    294           return false;
    295     return true;
    296   }
    297 
    298 
    299   const gsl_matrix* Matrix::gsl_matrix_p(void) const
    300   {
    301     return m_;
    302   }
    303 
    304 
    305   gsl_matrix* Matrix::gsl_matrix_p(void)
    306   {
    307     return m_;
    308   }
    309 
    310 
    311   void Matrix::mul(const Matrix& other)
    312   {
    313     assert(m_);
    314     int status=gsl_matrix_mul_elements(m_, other.gsl_matrix_p());
    315     if (status)
    316       throw utility::GSL_error(std::string("Matrix::mul_elements",status));
    317   }
    318 
    319 
    320   void Matrix::resize(size_t r, size_t c, double init_value)
    321   {
    322     delete_allocated_memory();
    323 
    324     assert( (r&&c) || (!r&&!c) );
    325     if (r && c){
    326       m_ = gsl_matrix_alloc(r,c);
    327       if (!m_)
    328         throw utility::GSL_error("Matrix::Matrix failed to allocate memory");
    329       all(init_value);
    330     }
    331 
    332     // no need to delete blas_result_ if the number of rows fit, it
    333     // may be useful later.
    334     if (blas_result_ && (blas_result_->size1!=rows())) {
    335       gsl_matrix_free(blas_result_);
    336       blas_result_=NULL;
    337     }
    338   }
    339 
    340 
    341   size_t Matrix::rows(void) const
    342   {
    343     return (m_ ? m_->size1 : 0);
    344   }
    345 
    346 
    347   const VectorConstView Matrix::row_const_view(size_t col) const
    348   {
    349     return VectorConstView(*this, col, true);
    350   }
    351 
    352 
    353   VectorView Matrix::row_view(size_t row)
    354   {
    355     VectorView res(*this, row, true);
    356     return res;
    357   }
    358 
    359 
    360   void Matrix::swap_columns(const size_t i, const size_t j)
    361   {
    362     assert(m_);
    363     int status=gsl_matrix_swap_columns(m_, i, j);
    364     if (status)
    365       throw utility::GSL_error(std::string("Matrix::swap_columns",status));
    366   }
    367 
    368 
    369   void Matrix::swap_rowcol(const size_t i, const size_t j)
    370   {
    371     assert(m_);
    372     int status=gsl_matrix_swap_rowcol(m_, i, j);
    373     if (status)
    374       throw utility::GSL_error(std::string("Matrix::swap_rowcol",status));
    375   }
    376 
    377 
    378   void Matrix::swap_rows(const size_t i, const size_t j)
    379   {
    380     assert(m_);
    381     int status=gsl_matrix_swap_rows(m_, i, j);
    382     if (status)
    383       throw utility::GSL_error(std::string("Matrix::swap_rows",status));
    384   }
    385 
    386 
    387   void Matrix::transpose(void)
    388   {
    389     assert(m_);
    390     if (columns()==rows())
    391       gsl_matrix_transpose(m_); // this never fails
    392     else {
    393       gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
    394       if (!transposed)
    395         throw utility::GSL_error("Matrix::transpose failed to allocate memory");
    396       // next line never fails if allocation above succeeded.
    397       gsl_matrix_transpose_memcpy(transposed,m_);
    398       gsl_matrix_free(m_);
    399       m_ = transposed;
    400       if (blas_result_) {
    401         gsl_matrix_free(blas_result_);
    402         blas_result_=NULL;
    403       }
    404     }
    405   }
    406 
    407 
    408   double& Matrix::operator()(size_t row, size_t column)
    409   {
    410     assert(m_);
    411     assert(row<rows());
    412     assert(column<columns());
    413     double* d=gsl_matrix_ptr(m_, row, column);
    414     if (!d)
    415       throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
    416     return *d;
    417   }
    418 
    419 
    420   const double& Matrix::operator()(size_t row, size_t column) const
    421   {
    422     assert(row<rows());
    423     assert(column<columns());
    424     const double* d=gsl_matrix_const_ptr(m_, row, column);
    425     if (!d)
    426       throw utility::GSL_error("Matrix::operator()",GSL_EINVAL);
    427     return *d;
    428   }
    429 
    430 
    431   bool Matrix::operator==(const Matrix& other) const
    432   {
    433     return equal(other);
    434   }
    435 
    436 
    437   bool Matrix::operator!=(const Matrix& other) const
    438   {
    439     return !equal(other);
    440   }
    441 
    442 
    443   const Matrix& Matrix::operator=( const Matrix& other )
    444   {
    445     assert(other.m_);
    446     if (this!=&other) {
    447       if ( !m_ || (other.m_->size1!=m_->size1) || (other.m_->size2!=m_->size2) )
    448         resize(other.m_->size1,other.m_->size2);
    449       if (gsl_matrix_memcpy(m_, other.gsl_matrix_p()))
    450         throw utility::GSL_error("Matrix::create_gsl_matrix_copy dimension mis-match");
    451     }
    452     return *this;
    453   }
    454 
    455 
    456   const Matrix& Matrix::operator+=(const Matrix& other)
    457   {
    458     assert(m_);
    459     int status=gsl_matrix_add(m_, other.m_);
    460     if (status)
    461       throw utility::GSL_error(std::string("Matrix::operator+=", status));
    462     return *this;
    463   }
    464 
    465 
    466   const Matrix& Matrix::operator+=(const double d)
    467   {
    468     assert(m_);
    469     gsl_matrix_add_constant(m_, d);
    470     return *this;
    471   }
    472 
    473 
    474   const Matrix& Matrix::operator-=(const Matrix& other)
    475   {
    476     assert(m_);
    477     int status=gsl_matrix_sub(m_, other.m_);
    478     if (status)
    479       throw utility::GSL_error(std::string("Matrix::operator-=", status));
    480     return *this;
    481   }
    482 
    483 
    484   const Matrix& Matrix::operator-=(const double d)
    485   {
    486     assert(m_);
    487     gsl_matrix_add_constant(m_, -d);
    488     return *this;
    489   }
    490 
    491 
    492   const Matrix& Matrix::operator*=(const Matrix& other)
    493   {
    494     assert(m_);
    495     if ( blas_result_ && ((blas_result_->size1!=rows()) ||
    496                           (blas_result_->size2!=other.columns())) ) {
    497       gsl_matrix_free(blas_result_);
    498       blas_result_=NULL;
    499     }
    500     if (!blas_result_) {
    501       blas_result_ = gsl_matrix_alloc(rows(),other.columns());
    502       if (!blas_result_)
    503         throw utility::GSL_error("Matrix::operator*= failed to allocate memory");
    504     }
    505     gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, blas_result_);
    506     gsl_matrix* tmp=m_;
    507     m_ = blas_result_;
    508     blas_result_=tmp;
    509     return *this;
    510   }
    511 
    512 
    513   const Matrix& Matrix::operator*=(const double d)
    514   {
    515     assert(m_);
    516     gsl_matrix_scale(m_, d);
    517     return *this;
    518   }
    519 
    520 
    521   bool isnull(const Matrix& other)
    522   {
    523     return gsl_matrix_isnull(other.gsl_matrix_p());
    524   }
    525 
    526 
    527   double max(const Matrix& other)
    528   {
    529     return gsl_matrix_max(other.gsl_matrix_p());
    530   }
    531 
    532 
    533   double min(const Matrix& other)
    534   {
    535     return gsl_matrix_min(other.gsl_matrix_p());
    536   }
    537 
    538 
    539   void minmax_index(const Matrix& other,
    540                     std::pair<size_t,size_t>& min, std::pair<size_t,size_t>& max)
    541   {
    542     gsl_matrix_minmax_index(other.gsl_matrix_p(), &min.first, &min.second,
    543                             &max.first, &max.second);
    544   }
    545 
    546 
    547   bool nan(const Matrix& templat, Matrix& flag)
    548   {
    549     size_t rows=templat.rows();
    550     size_t columns=templat.columns();
    551     if (rows!=flag.rows() && columns!=flag.columns())
    552       flag.resize(rows,columns,1.0);
    553     else
    554       flag.all(1.0);
    555     bool nan=false;
    556     for (size_t i=0; i<rows; i++)
    557       for (size_t j=0; j<columns; j++)
    558         if (std::isnan(templat(i,j))) {
    559           flag(i,j)=0;
    560           nan=true;
    561         }
    562     return nan;
    563   }
    564 
    565 
    566   void swap(Matrix& a, Matrix& b)
    567   {
    568     assert(a.gsl_matrix_p()); assert(b.gsl_matrix_p());
    569     int status=gsl_matrix_swap(a.gsl_matrix_p(), b.gsl_matrix_p());
    570     if (status)
    571       throw utility::GSL_error(std::string("swap(Matrix&,Matrix&)",status));
    572   }
    573 
    574 
    575   std::ostream& operator<<(std::ostream& s, const Matrix& m)
    576   {
    577     s.setf(std::ios::dec);
    578     s.precision(12);
    579     for(size_t i=0, j=0; i<m.rows(); i++)
    580       for (j=0; j<m.columns(); j++) {
    581         s << m(i,j);
    582         if (j<m.columns()-1)
    583           s << s.fill();
    584         else if (i<m.rows()-1)
    585           s << "\n";
    586       }
    587     return s;
    588   }
    589 
    590 
    591   Vector operator*(const Matrix& m, const VectorBase& v)
    592   {
    593     utility::Vector res(m.rows());
    594     for (size_t i=0; i<res.size(); ++i)
    595       res(i) = VectorConstView(m,i) * v;
    596     return res;
    597   }
    598 
    599 
    600   Vector operator*(const VectorBase& v, const Matrix& m)
    601   {
    602     utility::Vector res(m.columns());
    603     for (size_t i=0; i<res.size(); ++i)
    604       res(i) = v * VectorConstView(m,i,false);
    605     return res;
    606   }
    607 
    60843}}} // of namespace utility, yat and thep
Note: See TracChangeset for help on using the changeset viewer.