plugins/base1/se.lu.thep.wenni/trunk/lib/c++_tools/gslapi/matrix.cc

Code
Comments
Other
Rev Date Author Line
69 11 Feb 06 jari 1 // $Id$
69 11 Feb 06 jari 2
95 05 Apr 06 jari 3 /*
95 05 Apr 06 jari 4   Copyright (C) 2003 Daniel Dalevi, Peter Johansson
95 05 Apr 06 jari 5   Copyright (C) 2004 Jari Häkkinen, Peter Johansson
95 05 Apr 06 jari 6   Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
95 05 Apr 06 jari 7
95 05 Apr 06 jari 8   This file is part of the thep c++ tools library,
95 05 Apr 06 jari 9                                 http://lev.thep.lu.se/trac/c++_tools
95 05 Apr 06 jari 10
95 05 Apr 06 jari 11   The c++ tools library is free software; you can redistribute it
95 05 Apr 06 jari 12   and/or modify it under the terms of the GNU General Public License
824 26 Nov 08 jari 13   as published by the Free Software Foundation; either version 3 of
95 05 Apr 06 jari 14   the License, or (at your option) any later version.
95 05 Apr 06 jari 15
95 05 Apr 06 jari 16   The c++ tools library is distributed in the hope that it will be
95 05 Apr 06 jari 17   useful, but WITHOUT ANY WARRANTY; without even the implied warranty
95 05 Apr 06 jari 18   of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
95 05 Apr 06 jari 19   General Public License for more details.
95 05 Apr 06 jari 20
95 05 Apr 06 jari 21   You should have received a copy of the GNU General Public License
824 26 Nov 08 jari 22   along with WeNNI. If not, see <http://www.gnu.org/licenses/>.
95 05 Apr 06 jari 23 */
95 05 Apr 06 jari 24
69 11 Feb 06 jari 25 #include <c++_tools/gslapi/matrix.h>
69 11 Feb 06 jari 26
69 11 Feb 06 jari 27 #include <c++_tools/gslapi/vector.h>
69 11 Feb 06 jari 28 #include <c++_tools/utility/stl_utility.h>
110 13 Jun 06 jari 29 #include <c++_tools/utility/utility.h>
69 11 Feb 06 jari 30
69 11 Feb 06 jari 31 #include <cmath>
69 11 Feb 06 jari 32 #include <sstream>
69 11 Feb 06 jari 33 #include <vector>
69 11 Feb 06 jari 34
69 11 Feb 06 jari 35 #include <gsl/gsl_blas.h>
69 11 Feb 06 jari 36
69 11 Feb 06 jari 37
69 11 Feb 06 jari 38 namespace theplu {
69 11 Feb 06 jari 39 namespace gslapi {
69 11 Feb 06 jari 40
69 11 Feb 06 jari 41
69 11 Feb 06 jari 42
110 13 Jun 06 jari 43   matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
110 13 Jun 06 jari 44                  size_t n_row, size_t n_column)
69 11 Feb 06 jari 45   {
110 13 Jun 06 jari 46     // Jari, exception handling needed here. Failure in setting up a
110 13 Jun 06 jari 47     // proper gsl_matrix_view is communicated by NULL pointer in the
110 13 Jun 06 jari 48     // view structure (cf. GSL manual). How about GSL error state?
110 13 Jun 06 jari 49     view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
110 13 Jun 06 jari 50                                                      offset_row,offset_column,
110 13 Jun 06 jari 51                                                      n_row,n_column));
110 13 Jun 06 jari 52     m_ = &(view_->matrix);
69 11 Feb 06 jari 53   }
69 11 Feb 06 jari 54
69 11 Feb 06 jari 55
110 13 Jun 06 jari 56   // Constructor that gets data from istream
110 13 Jun 06 jari 57   matrix::matrix(std::istream& is, char sep) 
110 13 Jun 06 jari 58     throw (utility::IO_error,std::exception)
110 13 Jun 06 jari 59     : view_(NULL)
69 11 Feb 06 jari 60   {
110 13 Jun 06 jari 61     // Markus to Jari, somewhere we should check that quiet_NaNs are supported
110 13 Jun 06 jari 62     // std::numeric_limits<double>::has_quiet_NaN has to be true.
110 13 Jun 06 jari 63     // Also in vector
69 11 Feb 06 jari 64
69 11 Feb 06 jari 65     // read the data file and store in stl vectors (dynamically
69 11 Feb 06 jari 66     // expandable)
69 11 Feb 06 jari 67     std::vector<std::vector<double> > data_matrix;
819 24 Nov 08 jari 68     unsigned int nof_columns=0;
819 24 Nov 08 jari 69     unsigned int nof_rows = 0;
110 13 Jun 06 jari 70     std::string line;
110 13 Jun 06 jari 71     while(getline(is, line, '\n')){
69 11 Feb 06 jari 72       // Ignoring empty lines
110 13 Jun 06 jari 73       if (!line.size()) {
69 11 Feb 06 jari 74         continue;
69 11 Feb 06 jari 75       }
110 13 Jun 06 jari 76       nof_rows++;
110 13 Jun 06 jari 77       std::vector<double> v;
110 13 Jun 06 jari 78       std::string element;
110 13 Jun 06 jari 79       std::stringstream ss(line);
69 11 Feb 06 jari 80       
110 13 Jun 06 jari 81       bool ok=true;
110 13 Jun 06 jari 82       while(ok) {
110 13 Jun 06 jari 83         if(sep=='\0')
110 13 Jun 06 jari 84           ok=(ss>>element);
110 13 Jun 06 jari 85         else
110 13 Jun 06 jari 86           ok=getline(ss, element, sep);
110 13 Jun 06 jari 87         if(!ok)
110 13 Jun 06 jari 88           break;
110 13 Jun 06 jari 89         
110 13 Jun 06 jari 90         if(utility::is_double(element)) {
110 13 Jun 06 jari 91           v.push_back(atof(element.c_str()));
110 13 Jun 06 jari 92         }
110 13 Jun 06 jari 93         else if (!element.size() || utility::is_nan(element)) {
110 13 Jun 06 jari 94           v.push_back(std::numeric_limits<double>::quiet_NaN());
110 13 Jun 06 jari 95         }
110 13 Jun 06 jari 96         else {
110 13 Jun 06 jari 97           // Jari, this should be communicated with as an exception.
110 13 Jun 06 jari 98           std::cerr << "Warning: '" << element 
110 13 Jun 06 jari 99                     << "' is not accepted as a matrix element." << std::endl;
110 13 Jun 06 jari 100         }
110 13 Jun 06 jari 101       }             
110 13 Jun 06 jari 102       if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
110 13 Jun 06 jari 103           v.push_back(std::numeric_limits<double>::quiet_NaN());
110 13 Jun 06 jari 104       if (!nof_columns)
69 11 Feb 06 jari 105         nof_columns=v.size();
110 13 Jun 06 jari 106       else if (v.size()!=nof_columns) {
110 13 Jun 06 jari 107         std::ostringstream s;
110 13 Jun 06 jari 108         s << "matrix::matrix(std::istream&, char) data file error: "
110 13 Jun 06 jari 109           << "line " << nof_rows << " has " << v.size()
110 13 Jun 06 jari 110           << " columns; expected " << nof_columns  << " columns.";
110 13 Jun 06 jari 111         throw utility::IO_error(s.str());
69 11 Feb 06 jari 112       }
69 11 Feb 06 jari 113       data_matrix.push_back(v);
69 11 Feb 06 jari 114     }
69 11 Feb 06 jari 115
69 11 Feb 06 jari 116     // manipulate the state of the stream to be good
69 11 Feb 06 jari 117     is.clear(std::ios::goodbit);
69 11 Feb 06 jari 118     // convert the data to a gsl matrix
69 11 Feb 06 jari 119     m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
819 24 Nov 08 jari 120     for(unsigned int i=0;i<nof_rows;i++)
819 24 Nov 08 jari 121       for(unsigned int j=0;j<nof_columns;j++)
69 11 Feb 06 jari 122         gsl_matrix_set( m_, i, j, data_matrix[i][j] );
69 11 Feb 06 jari 123   }
69 11 Feb 06 jari 124
69 11 Feb 06 jari 125
69 11 Feb 06 jari 126   matrix::~matrix(void)
69 11 Feb 06 jari 127   {
110 13 Jun 06 jari 128     if (view_)
110 13 Jun 06 jari 129       delete view_;
110 13 Jun 06 jari 130     else if (m_)
69 11 Feb 06 jari 131       gsl_matrix_free(m_);
110 13 Jun 06 jari 132     m_=NULL;
69 11 Feb 06 jari 133   }
69 11 Feb 06 jari 134
69 11 Feb 06 jari 135
69 11 Feb 06 jari 136
69 11 Feb 06 jari 137   bool matrix::equal(const matrix& other, const double d) const
69 11 Feb 06 jari 138   {
69 11 Feb 06 jari 139     if (columns()!=other.columns() || rows()!=other.rows())
69 11 Feb 06 jari 140       return false;
110 13 Jun 06 jari 141     for (size_t i=0; i<rows(); i++)
110 13 Jun 06 jari 142       for (size_t j=0; j<columns(); j++)
110 13 Jun 06 jari 143         // The two last condition checks are needed for NaN detection
110 13 Jun 06 jari 144         if (fabs( (*this)(i,j)-other(i,j) ) > d ||
69 11 Feb 06 jari 145             (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
69 11 Feb 06 jari 146           return false;
69 11 Feb 06 jari 147     return true;
69 11 Feb 06 jari 148   }
69 11 Feb 06 jari 149
69 11 Feb 06 jari 150
69 11 Feb 06 jari 151
110 13 Jun 06 jari 152   gsl_matrix* matrix::create_gsl_matrix_copy(void) const
69 11 Feb 06 jari 153   {
110 13 Jun 06 jari 154     gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
110 13 Jun 06 jari 155     gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
110 13 Jun 06 jari 156     return m;
69 11 Feb 06 jari 157   }
69 11 Feb 06 jari 158
69 11 Feb 06 jari 159
69 11 Feb 06 jari 160
110 13 Jun 06 jari 161   matrix matrix::nan(void) const
69 11 Feb 06 jari 162   {
110 13 Jun 06 jari 163     matrix m(rows(),columns(),1.0);
110 13 Jun 06 jari 164     for (size_t i=0; i<rows(); i++)
110 13 Jun 06 jari 165       for (size_t j=0; j<columns(); j++) 
110 13 Jun 06 jari 166         if (std::isnan(operator()(i,j)))
110 13 Jun 06 jari 167           m(i,j)=0;
110 13 Jun 06 jari 168     return m;
69 11 Feb 06 jari 169   }
69 11 Feb 06 jari 170
69 11 Feb 06 jari 171
69 11 Feb 06 jari 172
110 13 Jun 06 jari 173   // Jari, checkout GSL transpose support in GSL manual 8.4.9
110 13 Jun 06 jari 174   void matrix::transpose(void)
69 11 Feb 06 jari 175   {
110 13 Jun 06 jari 176     if (columns()==rows())
110 13 Jun 06 jari 177       gsl_matrix_transpose(m_);
110 13 Jun 06 jari 178     else {
110 13 Jun 06 jari 179       gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
110 13 Jun 06 jari 180       gsl_matrix_transpose_memcpy(transposed,m_);
110 13 Jun 06 jari 181       gsl_matrix_free(m_);
110 13 Jun 06 jari 182       m_=transposed;
110 13 Jun 06 jari 183     }
69 11 Feb 06 jari 184   }
69 11 Feb 06 jari 185
69 11 Feb 06 jari 186
69 11 Feb 06 jari 187
110 13 Jun 06 jari 188   const vector matrix::operator*(const vector&) const
69 11 Feb 06 jari 189   {
110 13 Jun 06 jari 190     std::cerr << "Not implemented:" << std::endl
110 13 Jun 06 jari 191               << "   const vector matrix::operator*(const vector&) const"
110 13 Jun 06 jari 192               << std::endl;
110 13 Jun 06 jari 193     return vector(0);
69 11 Feb 06 jari 194   }
69 11 Feb 06 jari 195
69 11 Feb 06 jari 196
69 11 Feb 06 jari 197
110 13 Jun 06 jari 198   const matrix& matrix::operator=( const matrix& other )
69 11 Feb 06 jari 199   {
69 11 Feb 06 jari 200     if ( this != &other ) {
110 13 Jun 06 jari 201       if (view_)
110 13 Jun 06 jari 202         delete view_;
110 13 Jun 06 jari 203       else if (m_)
110 13 Jun 06 jari 204         gsl_matrix_free(m_);
110 13 Jun 06 jari 205       m_ = other.create_gsl_matrix_copy();
110 13 Jun 06 jari 206     }
69 11 Feb 06 jari 207     return *this;
69 11 Feb 06 jari 208   }
69 11 Feb 06 jari 209
69 11 Feb 06 jari 210
69 11 Feb 06 jari 211
110 13 Jun 06 jari 212   const matrix& matrix::operator*=(const matrix& other)
69 11 Feb 06 jari 213   {
110 13 Jun 06 jari 214     gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
110 13 Jun 06 jari 215     gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
110 13 Jun 06 jari 216     gsl_matrix_free(m_);
110 13 Jun 06 jari 217     m_=result;
110 13 Jun 06 jari 218     return *this;
69 11 Feb 06 jari 219   }
69 11 Feb 06 jari 220
69 11 Feb 06 jari 221
69 11 Feb 06 jari 222
69 11 Feb 06 jari 223   std::ostream& operator<<(std::ostream& s, const matrix& m)
69 11 Feb 06 jari 224   {
110 13 Jun 06 jari 225     s.setf(std::ios::dec);
69 11 Feb 06 jari 226     s.precision(12);
110 13 Jun 06 jari 227     for(size_t i=0, j=0; i<m.rows(); i++)
69 11 Feb 06 jari 228       for (j=0; j<m.columns(); j++) {
69 11 Feb 06 jari 229         s << m(i,j);
69 11 Feb 06 jari 230         if (j<m.columns()-1)
69 11 Feb 06 jari 231           s << "\t";
69 11 Feb 06 jari 232         else if (i<m.rows()-1)
69 11 Feb 06 jari 233           s << "\n";
69 11 Feb 06 jari 234       }
110 13 Jun 06 jari 235     return s;
69 11 Feb 06 jari 236   }
69 11 Feb 06 jari 237
69 11 Feb 06 jari 238
69 11 Feb 06 jari 239 }} // of namespace gslapi and namespace thep

Generated on Tue Nov 21 03:42:20 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)