source: trunk/lib/gslapi/matrix.cc @ 516

Last change on this file since 516 was 516, checked in by Peter, 16 years ago

added function generating weight matrix with respect to not-a-number vs number

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.0 KB
Line 
1// $Id: matrix.cc 516 2006-02-20 18:37:42Z peter $
2
3#include <c++_tools/gslapi/matrix.h>
4
5#include <c++_tools/gslapi/vector.h>
6#include <c++_tools/utility/stl_utility.h>
7#include <c++_tools/utility/utility.h>
8
9#include <cmath>
10#include <sstream>
11#include <vector>
12
13#include <gsl/gsl_blas.h>
14
15
16namespace theplu {
17namespace gslapi {
18
19
20
21  matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
22                 size_t n_row, size_t n_column)
23  {
24    // Jari, exception handling needed here. Failure in setting up a
25    // proper gsl_matrix_view is communicated by NULL pointer in the
26    // view structure (cf. GSL manual). How about GSL error state?
27    view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
28                                                     offset_row,offset_column,
29                                                     n_row,n_column));
30    m_ = &(view_->matrix);
31  }
32
33
34  // Constructor that gets data from istream
35  matrix::matrix(std::istream& is, char sep) 
36    throw (utility::IO_error,std::exception)
37    : view_(NULL)
38  {
39    // Markus to Jari, somewhere we should check that quiet_NaNs are supported
40    // std::numeric_limits<double>::has_quiet_NaN has to be true.
41    // Also in vector
42
43    // read the data file and store in stl vectors (dynamically
44    // expandable)
45    std::vector<std::vector<double> > data_matrix;
46    u_int nof_columns=0;
47    u_int nof_rows = 0;
48    std::string line;
49    while(getline(is, line, '\n')){
50      // Ignoring empty lines
51      if (!line.size()) {
52        continue;
53      }
54      nof_rows++;
55      std::vector<double> v;
56      std::string element;
57      std::stringstream ss(line);
58     
59      bool ok=true;
60      while(ok) {
61        if(sep=='\0')
62          ok=(ss>>element);
63        else
64          ok=getline(ss, element, sep);
65        if(!ok)
66          break;
67
68        if(utility::is_double(element)) {
69          v.push_back(atof(element.c_str()));
70        }
71        else if (!element.size() || utility::is_nan(element)) {
72          v.push_back(std::numeric_limits<double>::quiet_NaN());
73        }
74        else {
75          // Jari, this should be communicated with as an exception.
76          //          std::cerr << "Warning: '" << vec_str[i]
77          //                    << "' is not an integer." << std::endl;
78        }
79      }           
80      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
81          v.push_back(std::numeric_limits<double>::quiet_NaN());
82      if (!nof_columns)
83        nof_columns=v.size();
84      else if (v.size()!=nof_columns) {
85        std::ostringstream s;
86        s << "matrix::matrix(std::istream&, char) data file error: "
87          << "line " << nof_rows << " has " << v.size()
88          << " columns; expected " << nof_columns << " columns.";
89        throw utility::IO_error(s.str());
90      }
91      data_matrix.push_back(v);
92    }
93
94    // manipulate the state of the stream to be good
95    is.clear(std::ios::goodbit);
96    // convert the data to a gsl matrix
97    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
98    for(u_int i=0;i<nof_rows;i++)
99      for(u_int j=0;j<nof_columns;j++)
100        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
101  }
102
103
104  matrix::~matrix(void)
105  {
106    if (view_)
107      delete view_;
108    else if (m_)
109      gsl_matrix_free(m_);
110    m_=NULL;
111  }
112
113
114
115  bool matrix::equal(const matrix& other, const double d) const
116  {
117    if (columns()!=other.columns() || rows()!=other.rows())
118      return false;
119    for (size_t i=0; i<rows(); i++)
120      for (size_t j=0; j<columns(); j++)
121        // The two last condition checks are needed for NaN detection
122        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
123            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
124          return false;
125    return true;
126  }
127
128
129
130  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
131  {
132    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
133    gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
134    return m;
135  }
136
137
138
139  matrix matrix::nan(void) const
140  {
141    matrix m(rows(),columns(),1.0);
142    for (size_t i=0; i<rows(); i++)
143      for (size_t j=0; j<columns(); j++) 
144        if (std::isnan(operator()(i,j)))
145          m(i,j)=0;
146    return m;
147  }
148
149
150
151  // Jari, checkout GSL transpose support in GSL manual 8.4.9
152  void matrix::transpose(void)
153  {
154    if (columns()==rows())
155      gsl_matrix_transpose(m_);
156    else {
157      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
158      gsl_matrix_transpose_memcpy(transposed,m_);
159      gsl_matrix_free(m_);
160      m_=transposed;
161    }
162  }
163
164
165
166  const vector matrix::operator*(const vector&) const
167  {
168    std::cerr << "Not implemented:" << std::endl
169              << "   const vector matrix::operator*(const vector&) const"
170              << std::endl;
171    return vector(0);
172  }
173
174
175
176  const matrix& matrix::operator=( const matrix& other )
177  {
178    if ( this != &other ) {
179      if (view_)
180        delete view_;
181      else if (m_)
182        gsl_matrix_free(m_);
183      m_ = other.create_gsl_matrix_copy();
184    }
185    return *this;
186  }
187
188
189
190  const matrix& matrix::operator*=(const matrix& other)
191  {
192    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
193    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
194    gsl_matrix_free(m_);
195    m_=result;
196    return *this;
197  }
198
199
200
201  std::ostream& operator<<(std::ostream& s, const matrix& m)
202  {
203    s.setf(std::ios::dec);
204    s.precision(12);
205    for(size_t i=0, j=0; i<m.rows(); i++)
206      for (j=0; j<m.columns(); j++) {
207        s << m(i,j);
208        if (j<m.columns()-1)
209          s << "\t";
210        else if (i<m.rows()-1)
211          s << "\n";
212      }
213    return s;
214  }
215
216
217}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.