mixmod  3.2.0
Mixture models for clustering and classification
 All Classes Namespaces Files Functions Variables Enumerations Friends
Eigen.h
1 /***************************************************************************
2  SRC/mixmod/Utilities/maths/Eigen.h description
3  copyright : (C) MIXMOD Team - 2001-2016
4  email : contact@mixmod.org
5  ***************************************************************************/
6 
7 /***************************************************************************
8  This file is part of MIXMOD
9 
10  MIXMOD is free software: you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  MIXMOD is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with MIXMOD. If not, see <http://www.gnu.org/licenses/>.
22 
23  All informations available on : http://www.mixmod.org
24 ***************************************************************************/
25 #ifndef XEM_MATH_EIGEN_H
26 #define XEM_MATH_EIGEN_H
27 
28 #include <Eigen/Dense>
29 
30 namespace XEM {
31 namespace MATH {
32 
33 // TODO: copy constructors ?
34 // TODO: fix this, especially inverse() and SVD() [segmentation fault on NR tests]
35 
37 
38 public:
39 
40  DiagonalMatrix(int dim) {
41  _dim = dim;
42  _value = new double[dim];
43  }
44 
45  ~DiagonalMatrix() {
46  if (_value)
47  delete [] _value;
48  }
49 
50  double* Store() {
51  return _value;
52  }
53 
54  int Nrow() {
55  return _dim;
56  }
57 
58 private:
59 
60  int64_t _dim;
61  double* _value;
62 };
63 
64 class Matrix {
65 
66 public:
67 
68  Matrix(int nrow, int ncol) {
69  _value = new Eigen::MatrixXd(nrow, ncol);
70  }
71 
72  ~Matrix() {
73  if (_value)
74  delete _value;
75  }
76 
77  double* Store() {
78  return _value->data();
79  }
80 
81  double* GetRow(int index) {
82  return _value->data() + index * _value->cols();
83  }
84 
85  int Nrow() {
86  return _value->rows();
87  }
88 
89  int Ncol() {
90  return _value->cols();
91  }
92 
93 private:
94 
95  Eigen::MatrixXd* _value;
96 };
97 
99 
100 public:
101 
102  // nrow == ncol
103  SymmetricMatrix(int nrow) {
104  _value = new Eigen::MatrixXd(nrow, nrow);
105  }
106 
107  ~SymmetricMatrix() {
108  if (_value)
109  delete _value;
110  }
111 
112  double* Store() {
113  return _value->data();
114  }
115 
116  // return log(abs(det(M)))
117  double LogDeterminant() {
118  return log(_value->determinant());
119  }
120 
121  // get inverse
122  SymmetricMatrix* Inverse() {
123  int _s_pbDimension = _value->rows();
124  SymmetricMatrix* invMat = new SymmetricMatrix(_s_pbDimension);
125  auto eigenInverse = (Eigen::MatrixXd)(_value->inverse());
126  for (int64_t i=0; i<eigenInverse.rows()*eigenInverse.cols(); i++)
127  invMat->Store()[i] = eigenInverse.data()[i];
128  return invMat;
129  }
130 
131  // compute SVD (only matrices U and D, not V)
132  void computeSVD(DiagonalMatrix* D, Matrix* U) {
133  // TODO: Eigen::ComputeThinU ?
134  auto svd = new Eigen::JacobiSVD<Eigen::MatrixXd>(*_value, Eigen::ComputeFullU);
135  auto eigenD = svd->singularValues();
136  auto eigenU = svd->matrixU();
137  for (int64_t i=0; i<eigenD.rows(); i++)
138  D->Store()[i] = eigenD.data()[i];
139  for (int64_t i=0; i<eigenU.rows()*eigenU.cols(); i++)
140  U->Store()[i] = eigenU.data()[i];
141  delete svd;
142  }
143 
144 private:
145 
146  Eigen::MatrixXd* _value;
147 };
148 
149 }
150 }
151 
152 #endif
Definition: Eigen.h:98
Definition: Eigen.h:64
Definition: Eigen.h:36