gsl_matrix.hh
Go to the documentation of this file.
1 /* -*- mia-c++ -*-
2  *
3  * This file is part of MIA - a toolbox for medical image analysis
4  * Copyright (c) Leipzig, Madrid 1999-2016 Gert Wollny
5  *
6  * MIA is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with MIA; if not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 #ifndef GSLPP_MATRIX_HH
22 #define GSLPP_MATRIX_HH
23 
24 #include <cassert>
25 #include <iterator>
26 #include <gsl/gsl_matrix.h>
27 #include <mia/core/gsl_vector.hh>
28 
29 namespace gsl {
30 
32  friend class const_matrix_iterator;
33 public:
34  matrix_iterator(gsl_matrix *m, bool begin):
35  m_matrix(m),
36  m_current(begin ? m->data : m->data + m->size1 * m->tda),
37  m_current_column(0),
38  m_row_jump(m->tda - m->size2)
39  {
40  }
41 
43  m_matrix(nullptr),
44  m_current(nullptr),
45  m_current_column(0)
46  {
47  }
48 
49 
50  matrix_iterator(const matrix_iterator& other) = default;
51 
52  double& operator *(){
53  assert(m_current);
54  return *m_current;
55  }
56 
57  matrix_iterator& operator ++() {
58  assert(m_current);
59  ++m_current;
60  ++m_current_column;
61  if(m_current_column == m_matrix->size2) {
62  m_current += m_row_jump;
63  m_current_column = 0;
64  }
65  return *this;
66  }
67 
68 
69  matrix_iterator operator ++(int) {
70  matrix_iterator result(*this);
71  ++(*this);
72  return result;
73  }
74 
75  friend bool operator == (const matrix_iterator& lhs, const matrix_iterator& rhs);
76  friend bool operator != (const matrix_iterator& lhs, const matrix_iterator& rhs);
77 private:
78  gsl_matrix *m_matrix;
79  double *m_current;
80  size_t m_current_column;
81  size_t m_row_jump;
82 };
83 
84 bool EXPORT_GSL operator == (const matrix_iterator& lhs, const matrix_iterator& rhs);
85 bool EXPORT_GSL operator != (const matrix_iterator& lhs, const matrix_iterator& rhs);
86 
87 
89 public:
90  const_matrix_iterator(const gsl_matrix *m, bool begin):
91  m_matrix(m),
92  m_current(begin ? m->data : m->data + m->size1 * m->tda),
93  m_current_column(0),
94  m_row_jump(m->tda - m->size2)
95  {
96  }
97 
99  m_matrix(other.m_matrix),
100  m_current(other.m_current),
101  m_current_column(other.m_current_column),
102  m_row_jump(other.m_row_jump)
103  {
104  }
105 
107  m_matrix(nullptr),
108  m_current(nullptr),
109  m_current_column(0)
110  {
111  }
112 
113 
114  const_matrix_iterator(const const_matrix_iterator& other) = default;
115 
116  double operator *() const{
117  assert(m_current);
118  return *m_current;
119  }
120 
121  const_matrix_iterator& operator ++() {
122  assert(m_current);
123  ++m_current;
124  ++m_current_column;
125  if(m_current_column == m_matrix->size2) {
126  m_current += m_row_jump;
127  m_current_column = 0;
128  }
129  return *this;
130  }
131 
132  const_matrix_iterator operator ++(int) {
133  const_matrix_iterator result(*this);
134  ++(*this);
135  return result;
136  }
137 
138  friend bool operator == (const const_matrix_iterator& lhs, const const_matrix_iterator& rhs);
139  friend bool operator != (const const_matrix_iterator& lhs, const const_matrix_iterator& rhs);
140 private:
141  const gsl_matrix *m_matrix;
142  const double *m_current;
143  size_t m_current_column;
144  size_t m_row_jump;
145 };
146 
149 
150 
151 
152 
153 
159 public:
162 
163  Matrix();
164 
171  Matrix(size_t rows, size_t columns, bool clean);
172 
179  Matrix(size_t rows, size_t columns, double init);
180 
188  Matrix(size_t rows, size_t columns, const double *init);
189 
193  Matrix(const Matrix& other);
194 
199  Matrix(gsl_matrix* m);
200 
205  Matrix(const gsl_matrix* m);
206 
211  Matrix transposed() const;
212 
213 
217  Matrix& operator =(const Matrix& other);
218 
225  void reset(size_t rows, size_t columns, bool clean);
226 
233  void reset(size_t rows, size_t columns, double init);
234 
235  ~Matrix();
236 
237 
241  size_t rows()const;
242 
246  size_t cols()const;
247 
254  void set(size_t i, size_t j, double x);
255 
263  double operator ()(size_t i, size_t j) const;
264 
268  operator gsl_matrix *();
269 
273  operator const gsl_matrix *() const;
274 
279  Matrix column_covariance() const;
280 
285  Matrix row_covariance() const;
286 
292  matrix_iterator begin();
293 
299  matrix_iterator end();
300 
306  const_matrix_iterator begin() const;
307 
313  const_matrix_iterator end() const;
314 
321  void set_row(int r, const Vector& row);
322 
328  VectorView get_row(int r);
329 
335  ConstVectorView get_row(int r) const;
336 
343  void set_column(int c, const Vector& col);
344 
350  VectorView get_column(int c);
351 
357  ConstVectorView get_column(int c) const;
358 
366  double dot_row(int r, const Vector& v) const;
367 
376  double dot_column(int c, const Vector& col) const;
377 
382  void print(std::ostream& os) const;
383 
389  Matrix& operator -=(const Matrix& rhs) {
390  gsl_matrix_sub(*this, rhs);
391  return *this;
392  }
393 
399  Matrix& operator +=(const Matrix& rhs) {
400  gsl_matrix_add(*this, rhs);
401  return *this;
402  }
403 
409  Matrix& operator *=(double rhs) {
410  gsl_matrix_scale(*this, rhs);
411  return *this;
412  }
413 
414 private:
415  gsl_matrix *m_matrix;
416  const gsl_matrix *m_const_matrix;
417  bool m_owner;
418 };
419 
420 
421 Matrix EXPORT_GSL operator * (const Matrix& lhs, const Matrix& rhs);
422 Matrix EXPORT_GSL operator + (const Matrix& lhs, const Matrix& rhs);
423 Matrix EXPORT_GSL operator - (const Matrix& lhs, const Matrix& rhs);
424 
425 inline std::ostream& operator << (std::ostream& os, const Matrix& m)
426 {
427  m.print(os);
428  return os;
429 }
430 
436 
439 };
440 
446 
447 
448 } // end namespace
449 
450 namespace std {
451 
452 template <>
453 class iterator_traits< gsl::const_matrix_iterator > {
454 public:
455  typedef size_t difference_type;
456  typedef double value_type;
457  typedef const double* pointer;
458  typedef const double& reference;
459  typedef forward_iterator_tag iterator_category;
460 };
461 
462 template <>
463 class iterator_traits< gsl::matrix_iterator > {
464 public:
465  typedef size_t difference_type;
466  typedef double value_type;
467  typedef double* pointer;
468  typedef double& reference;
469  typedef forward_iterator_tag iterator_category;
470 };
471 
472 }
473 
474 
475 #endif
void print(std::ostream &os) const
void EXPORT_GSL matrix_inv_sqrt(Matrix &m)
const_matrix_iterator const_iterator
Definition: gsl_matrix.hh:161
matrix_iterator iterator
Definition: gsl_matrix.hh:160
const_matrix_iterator(const gsl_matrix *m, bool begin)
Definition: gsl_matrix.hh:90
EXPORT_2D C2DFVectorfield & operator+=(C2DFVectorfield &a, const C2DFVectorfield &b)
bool EXPORT_GSL operator!=(const matrix_iterator &lhs, const matrix_iterator &rhs)
vector_iterator operator+(const vector_iterator &it, int dist)
vector_iterator operator-(const vector_iterator &it, int dist)
#define EXPORT_GSL
Definition: gsl_defines.hh:38
matrix_iterator(gsl_matrix *m, bool begin)
Definition: gsl_matrix.hh:34
std::ostream & operator<<(std::ostream &os, const Matrix &m)
Definition: gsl_matrix.hh:425
const_matrix_iterator(const matrix_iterator &other)
Definition: gsl_matrix.hh:98
Matrix EXPORT_GSL operator*(const Matrix &lhs, const Matrix &rhs)
bool EXPORT_GSL operator==(const matrix_iterator &lhs, const matrix_iterator &rhs)