Changeset 1380


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

working on #363

Location:
trunk/yat/utility
Files:
3 edited
2 copied

Legend:

Unmodified
Added
Removed

  • trunk/yat/utility/DataWeight.cc

    r1379 r1380  
    2323
    2424#include "DataWeight.h"
    25 
    26 #include <cassert>
    2725
    2826namespace theplu {

  • trunk/yat/utility/DataWeight.h

    r1379 r1380  
    2424  02111-1307, USA.
    2525*/
     26
     27#include "iterator_traits.h"
     28
     29#include <vector>
    2630
    2731namespace theplu {
    … …  
    113117
    114118
     119  /**
     120     Specialization for std::vector<DataWeight>::iterator
     121   */
     122  template<>
     123  struct weighted_iterator_traits<std::vector<DataWeight>::iterator> {
     124    /**
     125       iterator to vector<DataWeight> is weighted
     126    */
     127    typedef weighted_iterator_tag type;
     128  };
     129
     130
     131  /**
     132     Specialization for std::vector<DataWeight>::const_iterator
     133   */
     134  template<>
     135  struct weighted_iterator_traits<std::vector<DataWeight>::const_iterator> {
     136    /**
     137       const_iterator to vector<DataWeight> is weighted
     138    */
     139    typedef weighted_iterator_tag type;
     140  };
     141
    115142}}} // of namespace utility, yat, and theplu
    116143

  • trunk/yat/utility/Makefile.am

    r1379 r1380  
    2929  Alignment.cc ColumnStream.cc CommandLine.cc DataWeight.cc \
    3030  FileUtil.cc Index.cc kNNI.cc \
    31   Matrix.cc NNI.cc Option.cc OptionFile.cc OptionInFile.cc OptionOutFile.cc \
     31  Matrix.cc MatrixWeighted.cc NNI.cc Option.cc \
     32  OptionFile.cc OptionInFile.cc OptionOutFile.cc \
    3233  OptionHelp.cc OptionSwitch.cc \
    3334  PCA.cc stl_utility.cc SVD.cc TypeInfo.cc utility.cc Vector.cc \
    … …  
    4243  Exception.h FileUtil.h Index.h \
    4344  IteratorPolicy.h iterator_traits.h \
    44   kNNI.h Matrix.h NNI.h \
     45  kNNI.h Matrix.h MatrixWeighted.h NNI.h \
    4546  Option.h OptionArg.h OptionFile.h OptionInFile.h OptionOutFile.h \
    4647  OptionHelp.h OptionSwitch.h \

  • 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

  • trunk/yat/utility/MatrixWeighted.h

    r1360 r1380  
    1 #ifndef _theplu_yat_utility_matrix_
    2 #define _theplu_yat_utility_matrix_
     1#ifndef _theplu_yat_utility_matrix_weighted_
     2#define _theplu_yat_utility_matrix_weighted_
    33
    44// $Id$
    … …  
    2828*/
    2929
    30 #include "Exception.h"
     30#include "DataWeight.h"
    3131#include "StrideIterator.h"
    32 #include "Vector.h"
    33 #include "VectorConstView.h"
    34 #include "VectorView.h"
    35 
    36 #include <gsl/gsl_matrix.h>
    37 #include <iostream>
    38 #include <utility>
     32
     33#include <vector>
    3934
    4035namespace theplu {
    … …  
    4237namespace utility {
    4338
    44   class VectorBase;
    45 
    4639  ///
    47   /// @brief Interface to GSL matrix.
     40  /// @brief Weighted Matrix
    4841  ///
    49   /// For the time being 'double' is the only type supported.
    50   ///
    51   /// \par[File streams] Reading and writing vectors to file streams
    52   /// are of course supported. These are implemented without using GSL
    53   /// functionality, and thus binary read and write to streams are not
    54   /// supported.
    55   ///
    56   /// @note All GSL matrix related functions are not implement but
    57   /// most functionality defined for GSL matrices can be achieved with
    58   /// this interface class. Most notable GSL functionality not
    59   /// supported are no binary file support and views on arrays,
    60   /// utility::vectors, gsl_vectors, diagonals, subdiagonals, and
    61   /// superdiagonals.
    62   ///
    63   class Matrix
     42  class MatrixWeighted
    6443  {
    6544  public:
    … …  
    6746       Mutable iterator that iterates over all elements
    6847     */
    69     typedef StrideIterator<double*> iterator;
     48    typedef StrideIterator<std::vector<DataWeight>::iterator> iterator;
    7049
    7150    /**
    7251       Read-only iterator that iterates over all elements
    7352     */
    74     typedef StrideIterator<const double*> const_iterator;
     53    typedef StrideIterator<std::vector<DataWeight>::const_iterator>
     54    const_iterator;
    7555
    7656    /**
    7757       Mutable iterator that iterates over one column
    7858     */
    79     typedef StrideIterator<double*> column_iterator;
     59    typedef StrideIterator<std::vector<DataWeight>::iterator> column_iterator;
    8060
    8161    /**
    8262       Read-only iterator that iterates over one column
    8363     */
    84     typedef StrideIterator<const double*> const_column_iterator;
     64    typedef StrideIterator<std::vector<DataWeight>::const_iterator>
     65    const_column_iterator;
    8566
    8667    /**
    8768       Mutable iterator that iterates over one row
    8869     */
    89     typedef StrideIterator<double*> row_iterator;
     70    typedef StrideIterator<std::vector<DataWeight>::iterator> row_iterator;
    9071
    9172    /**
    9273       Read-only iterator that iterates over one row
    9374     */
    94     typedef StrideIterator<const double*> const_row_iterator;
     75    typedef StrideIterator<std::vector<DataWeight>::const_iterator>
     76    const_row_iterator;
    9577
    9678    /**
    … …  
    10082       structures.
    10183    */
    102     Matrix(void);
     84    MatrixWeighted(void);
    10385
    10486    /**
    10587       \brief Constructor allocating memory space for \a r times \a c
    106        elements, and sets all elements to \a init_value.
    107 
    108        \throw GSL_error if memory allocation fails.
    109     */
    110     Matrix(const size_t& r, const size_t& c, double init_value=0);
     88       elements, and sets all elements to
     89       DataWeight(\a init_value, \a init_weight)
     90    */
     91    MatrixWeighted(const size_t& r, const size_t& c, double init_value=0,
     92                 double init_weight=1.0);
    11193
    11294    /**
    11395       \brief The copy constructor.
    114 
    115        \throw A GSL_error is indirectly thrown if memory allocation
    116        fails.
    117     */
    118     Matrix(const Matrix&);
     96    */
     97    MatrixWeighted(const MatrixWeighted&);
    11998
    12099    /**
    121100       \brief The istream constructor.
    122101
    123        Matrix rows are sepearated with the new line character, and
    124        column elements in a row must be separated either with white
    125        space characters except the new line (\\n) character (default),
    126        or by the delimiter \a sep.  Rows, and columns, are read
    127        consecutively and only rectangular matrices are
    128        supported. Empty lines are ignored. End of input is at end of
    129        file marker.
    130 
    131        \throw GSL_error if memory allocation fails, IO_error if
    132        unexpected input is found in the input stream.
    133     */
    134     explicit Matrix(std::istream &, char sep='\0')
    135       throw(utility::IO_error, std::exception);
    136 
    137     /**
    138        \brief The destructor.
    139     */
    140     ~Matrix(void);
    141 
    142     ///
    143     /// Set all elements to \a value.
    144     ///
    145     void all(const double value);
     102       Constructor is reading data values using corresponding
     103       constructor in Matrix class. The weights are calculate dusing
     104       the nan function.
     105    */
     106    explicit MatrixWeighted(std::istream &, char sep='\0');
    146107
    147108    /**
    … …  
    190151
    191152    /**
    192        \return Vector view of column \a i
    193      */
    194     VectorView column_view(size_t i);
    195 
    196     /**
    197        \return const vector view of column \a i
    198      */
    199     const VectorConstView column_const_view(size_t) const;
    200 
    201     ///
    202     /// @return The number of columns in the matrix.
    203     ///
     153       \return The number of columns in the matrix.
     154    */
    204155    size_t columns(void) const;
    205156
    206157    /**
    207        Elementwise division of the elemnts of the calling matrix by
    208        the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
    209        \forall i,j \f$. The result is stored into the calling matrix.
    210 
    211        \throw GSL_error if dimensions mis-match.
    212     */
    213     void div(const Matrix& b);
    214 
    215     /**
    216158       \return iterator pointing to end of matrix
    217159     */
    … …  
    244186
    245187    /**
    246        \brief Check whether matrices are equal within a user defined
    247        precision, set by \a precision.
    248 
    249        \return True if each element deviates less or equal than \a
    250        d. If any matrix contain a NaN, false is always returned.
    251 
    252        \see operator== and operator!=
    253     */
    254     bool equal(const Matrix&, const double precision=0) const;
    255 
    256     ///
    257     /// @return A const pointer to the internal GSL matrix.
    258     ///
    259     const gsl_matrix* gsl_matrix_p(void) const;
    260 
    261     ///
    262     /// @return A pointer to the internal GSL matrix.
    263     ///
    264     gsl_matrix* gsl_matrix_p(void);
    265 
    266     /**
    267        Multiply the elements of Matrix \a b with the elements of the
    268        calling Matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall i,j
    269        \f$. The result is stored into the calling Matrix.
    270 
    271        \throw GSL_error if dimensions mis-match.
    272     */
    273     void mul(const Matrix& b);
    274 
    275     /**
    276188       \brief Resize Matrix
    277        
    278        All elements are set to @a init_value.
    279 
    280        \note underlying GSL matrix is destroyed and views into this
    281        Matrix becomes invalid.
    282     */
    283     void resize(size_t, size_t, double init_value=0);
    284 
    285     ///
    286     /// @return The number of rows in the matrix.
    287     ///
     189
     190       \note this function may invalidate iterators.
     191    */
     192    void resize(size_t, size_t);
     193
     194    /**
     195       \return The number of rows in the matrix.
     196    */
    288197    size_t rows(void) const;
    289198
    290199    /**
    291        \return Vector view of \a i
    292      */
    293     VectorView row_view(size_t);
    294 
    295     /**
    296        \return const vector view of \a i
    297      */
    298     const VectorConstView row_const_view(size_t) const;
    299 
    300     /**
    301        \brief Swap columns \a i and \a j.
    302 
    303        \throw GSL_error if either index is out of bounds.
    304     */
    305     void swap_columns(const size_t i,const size_t j);
    306 
    307     /**
    308        \brief Swap row \a i and column \a j.
    309 
    310        The Matrix must be square.
    311 
    312        \throw GSL_error if either index is out of bounds, or if matrix
    313        is not square.
    314     */
    315     void swap_rowcol(const size_t i,const size_t j);
    316 
    317     /**
    318        \brief Swap rows \a i and \a j.
    319 
    320        \throw GSL_error if either index is out of bounds.
    321     */
    322     void swap_rows(const size_t i, const size_t j);
    323 
    324     /**
    325200       \brief Transpose the matrix.
    326201
    327        \throw GSL_error if memory allocation fails for the new
    328        transposed matrix.
     202       \note This function may invalidate iterators.
    329203    */
    330204    void transpose(void);
    … …  
    334208
    335209       \return Reference to the element position (\a row, \a column).
    336 
    337        \throw If GSL range checks are enabled in the underlying GSL
    338        library a GSL_error exception is thrown if either index is out
    339        of range.
    340     */
    341     double& operator()(size_t row,size_t column);
     210    */
     211    DataWeight& operator()(size_t row,size_t column);
    342212
    343213    /**
    … …  
    346216       \return Const reference to the element position (\a row, \a
    347217       column).
    348 
    349        \throw If GSL range checks are enabled in the underlying GSL
    350        library a GSL_error exception is thrown if either index is out
    351        of range.
    352     */
    353     const double& operator()(size_t row,size_t column) const;
    354 
    355     /**
    356        \brief Comparison operator. Takes squared time.
    357 
    358        Checks are performed with exact matching, i.e., rounding off
    359        effects may destroy comparison. Use the equal function for
    360        comparing elements within a user defined precision.
    361 
    362        \return True if all elements are equal otherwise false.
    363 
    364        \see equal
    365     */
    366     bool operator==(const Matrix& other) const;
    367 
    368     /**
    369        \brief Comparison operator. Takes squared time.
    370 
    371        Checks are performed with exact matching, i.e., rounding off
    372        effects may destroy comparison. Use the equal function for
    373        comparing elements within a user defined precision.
    374 
    375        \return False if all elements are equal otherwise true.
    376 
    377        \see equal
    378     */
    379     bool operator!=(const Matrix& other) const;
     218    */
     219    const DataWeight& operator()(size_t row,size_t column) const;
     220
    380221
    381222    /**
    … …  
    384225       \return A const reference to the resulting Matrix.
    385226    */
    386     const Matrix& operator=(const Matrix& other);
    387 
    388     /**
    389        \brief Add and assign operator.
    390 
    391        Elementwise addition of the elements of Matrix \a b to the left
    392        hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall i,j
    393        \f$.
    394 
    395        \return A const reference to the resulting Matrix.
    396 
    397        \throw GSL_error if dimensions mis-match.
    398     */
    399     const Matrix& operator+=(const Matrix& b);
    400 
    401     /**
    402        \brief Add and assign operator
    403 
    404        Add the scalar value \a d to the left hand side Matrix, \f$
    405        a_{ij} = a_{ij} + d \; \forall i,j \f$.
    406     */
    407     const Matrix& operator+=(const double d);
    408 
    409     /**
    410        \brief Subtract and assign operator.
    411 
    412        Elementwise subtraction of the elements of Matrix \a b to the
    413        left hand side Matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \; \forall
    414        i,j \f$.
    415 
    416        \return A const reference to the resulting Matrix.
    417 
    418        \throw GSL_error if dimensions mis-match.
    419     */
    420     const Matrix& operator-=(const Matrix&);
    421 
    422     /**
    423        \brief Subtract and assign operator
    424 
    425        Subtract the scalar value \a d to the left hand side Matrix,
    426        \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$.
    427     */
    428     const Matrix& operator-=(const double d);
    429 
    430     /**
    431        \brief Multiply and assigment operator.
    432 
    433        \return Const reference to the resulting Matrix.
    434 
    435        \throw GSL_error if memory allocation fails.
    436     */
    437     const Matrix& operator*=(const Matrix&);
    438 
    439     /**
    440        \brief Multiply and assignment operator
    441 
    442        Multiply the elements of the left hand side Matrix with a
    443        scalar \a d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$.
    444 
    445        \throw GSL_error if memory allocation fails.
    446     */
    447     const Matrix& operator*=(double d);
     227    const MatrixWeighted& operator=(const MatrixWeighted& other);
    448228
    449229  private:
    450 
    451     /**
    452        \brief Create a new copy of the internal GSL matrix.
    453 
    454        Necessary memory for the new GSL matrix is allocated and the
    455        caller is responsible for freeing the allocated memory.
    456 
    457        \return A pointer to a copy of the internal GSL matrix.
    458 
    459        \throw GSL_error if memory cannot be allocated for the new
    460        copy, or if dimensions mis-match.
    461     */
    462     gsl_matrix* create_gsl_matrix_copy(void) const;
    463 
    464     /**
    465        \brief Clear all dynamically allocated memory.
    466 
    467        Internal utility function.
    468     */
    469     void delete_allocated_memory(void);
    470 
    471     // blas_result_ is used to temporarily store result in BLAS calls
    472     // and when not NULL it should always have the same size as
    473     // m_. Memory is not allocated for blas_result_ until it is
    474     // needed.
    475     gsl_matrix* blas_result_;
    476     gsl_matrix* m_;
     230    std::vector<DataWeight> vec_;
     231    size_t columns_;
     232
    477233  };
    478234
    479235  /**
    480      \brief Check if all elements of the Matrix are zero.
    481 
    482      \return True if all elements in the Matrix is zero, false
    483      othwerwise.
     236     \brief Exchange all elements.
     237
     238     Takes constant time.
     239
     240     \throw std::runtime_error if the two matrices disagree in size.
    484241  */
    485   bool isnull(const Matrix&);
    486 
    487   /**
    488      \brief Get the maximum value of the Matrix.
    489 
    490      \return The maximum value of the Matrix.
    491   */
    492   double max(const Matrix&);
    493 
    494   /**
    495      \brief Get the minimum value of the Matrix.
    496 
    497      \return The minimum value of the Matrix.
    498   */
    499   double min(const Matrix&);
    500 
    501   /**
    502      \brief Locate the maximum and minumum element in the Matrix.
    503 
    504      \return The indecies to the element with the minimum and maximum
    505      values of the Matrix, respectively.
    506 
    507      \note Lower index has precedence (searching in row-major order).
    508   */
    509   void minmax_index(const Matrix&,
    510                     std::pair<size_t,size_t>& min,
    511                     std::pair<size_t,size_t>& max);
    512 
    513   /**
    514      \brief Create a Matrix \a flag indicating NaN's in another Matrix
    515      \a templat.
    516 
    517      The \a flag Matrix is changed to contain 1's and 0's only. A 1
    518      means that the corresponding element in the \a templat Matrix is
    519      valid and a zero means that the corresponding element is a NaN.
    520 
    521      \note Space for Matrix \a flag is reallocated to fit the size of
    522      Matrix \a templat if sizes mismatch.
    523 
    524      \return True if the \a templat Matrix contains at least one NaN.
    525   */
    526   bool nan(const Matrix& templat, Matrix& flag);
    527 
    528   /**
    529      \brief Exchange all elements between the matrices by copying.
    530 
    531      The two matrices must have the same size.
    532 
    533      \throw GSL_error if either index is out of bounds.
    534   */
    535   void swap(Matrix&, Matrix&);
    536 
    537   /**
    538      \brief The output operator for the Matrix class.
    539   */
    540   std::ostream& operator<< (std::ostream& s, const Matrix&);
    541 
    542   /**
    543      \brief Vector Matrix multiplication
    544    */
    545   Vector operator*(const Matrix&, const VectorBase&);
    546 
    547   /**
    548      \brief Matrix Vector multiplication
    549    */
    550   Vector operator*(const VectorBase&, const Matrix&);
     242  void swap(MatrixWeighted&, MatrixWeighted&);
    551243
    552244}}} // of namespace utility, yat, and theplu
Note: See TracChangeset for help on using the changeset viewer.