mixmod  3.2.0
Mixture models for clustering and classification
 All Classes Namespaces Files Functions Variables Enumerations Friends
NEWMAT.h
1 /***************************************************************************
2  SRC/mixmod/Utilities/maths/NEWMAT.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_NEWMAT_H
26 #define XEM_MATH_NEWMAT_H
27 
28 #include "NEWMAT/newmatap.h"
29 #include "NEWMAT/newmatio.h"
30 
31 namespace XEM {
32 namespace MATH {
33 
34 // TODO: copy constructors ?
35 
36 class DiagonalMatrix {
37 
38 public:
39 
40  DiagonalMatrix(int dim) {
41  _value = new NEWMAT::DiagonalMatrix(dim);
42  }
43 
44  ~DiagonalMatrix() {
45  if (_value)
46  delete _value;
47  }
48 
49  double* Store() {
50  return _value->Store();
51  }
52 
53  int Nrow() {
54  return _value->Nrows();
55  }
56 
57 private:
58 
59  NEWMAT::DiagonalMatrix* _value;
60 };
61 
62 class Matrix {
63 
64 public:
65 
66  Matrix(int nrow, int ncol) {
67  _value = new NEWMAT::Matrix(nrow, ncol);
68  }
69 
70  ~Matrix() {
71  if (_value)
72  delete _value;
73  }
74 
75  double* Store() {
76  return _value->Store();
77  }
78 
79  double* GetRow(int index) {
80  return _value->Store() + index * _value->Ncols();
81  }
82 
83  int Nrow() {
84  return _value->Nrows();
85  }
86 
87  int Ncol() {
88  return _value->Ncols();
89  }
90 
91 private:
92 
93  NEWMAT::Matrix* _value;
94 };
95 
96 class SymmetricMatrix {
97 
98 public:
99 
100  // nrow == ncol
101  SymmetricMatrix(int nrow) {
102  _value = new NEWMAT::SymmetricMatrix(nrow);
103  }
104  // better 4 Inverse(), it avoids some memory leaks...
105  SymmetricMatrix(NEWMAT::SymmetricMatrix* mx) {
106  _value = mx;
107  }
108 
109  ~SymmetricMatrix() {
110  if (_value)
111  delete _value;
112  }
113 
114  double* Store() {
115  return _value->Store();
116  }
117 
118  // return log(abs(det(M)))
119  // NOTE: NEWMAT also return the sign of det(M)
120  double LogDeterminant() {
121  return NEWMAT::LogDeterminant(*_value).Value();
122  }
123 
124  // get inverse
125  SymmetricMatrix* Inverse() {
126  int _s_pbDimension = _value->Nrows();
127  NEWMAT::SymmetricMatrix* value_Inv = new NEWMAT::SymmetricMatrix(_s_pbDimension);
128  *value_Inv << _value->i();
129  //SymmetricMatrix* invMat = new SymmetricMatrix(_s_pbDimension);
130  //invMat->_value = value_Inv; //buggy, it generates memory leaks :-(
131  SymmetricMatrix* invMat = new SymmetricMatrix(value_Inv);
132  return invMat;
133  }
134 
135  // compute SVD (only matrices U and D, not V)
136  void computeSVD(DiagonalMatrix* D, Matrix* U) {
137  NEWMAT::DiagonalMatrix * D_NEWMAT = new NEWMAT::DiagonalMatrix(D->Nrow());
138  for (int64_t i=0; i<D->Nrow(); i++)
139  D_NEWMAT->Store()[i] = D->Store()[i];
140  NEWMAT::Matrix * U_NEWMAT = new NEWMAT::Matrix(U->Nrow(),U->Ncol());
141  for (int64_t i=0; i<U->Nrow()*U->Ncol(); i++)
142  U_NEWMAT->Store()[i] = U->Store()[i];
143  NEWMAT::SVD(*_value, *D_NEWMAT, *U_NEWMAT);
144  for (int64_t i=0; i<D->Nrow(); i++)
145  D->Store()[i] = D_NEWMAT->Store()[i];
146  for (int64_t i=0; i<U->Nrow()*U->Ncol(); i++)
147  U->Store()[i] = U_NEWMAT->Store()[i];
148  delete D_NEWMAT;
149  delete U_NEWMAT;
150  }
151 
152 private:
153 
154  NEWMAT::SymmetricMatrix* _value;
155 };
156 
157 }
158 }
159 
160 #endif