1 | // This file is part of the PACPUS framework distributed under the
|
---|
2 | // CECILL-C License, Version 1.0.
|
---|
3 | //
|
---|
4 | /// @file
|
---|
5 | /// @author Firstname Surname <firstname.surname@utc.fr>
|
---|
6 | /// @date Month, Year
|
---|
7 | /// @version $Id: ublas.hpp 73 2013-01-10 16:56:42Z kurdejma $
|
---|
8 | /// @copyright Copyright (c) UTC/CNRS Heudiasyc 2006 - 2013. All rights reserved.
|
---|
9 | /// @brief Brief description.
|
---|
10 | ///
|
---|
11 | /// Detailed description.
|
---|
12 |
|
---|
13 | #ifndef __BOOST_UBLAS__
|
---|
14 | #define __BOOST_UBLAS__
|
---|
15 |
|
---|
16 | #include <boost/numeric/ublas/vector.hpp>
|
---|
17 | #include <boost/numeric/ublas/matrix.hpp>
|
---|
18 |
|
---|
19 | #include <boost/numeric/ublas/lu.hpp>
|
---|
20 | #include <boost/numeric/ublas/vector_proxy.hpp>
|
---|
21 | #include <boost/numeric/ublas/triangular.hpp>
|
---|
22 | #include <boost/numeric/ublas/lu.hpp>
|
---|
23 | #include <boost/numeric/ublas/io.hpp>
|
---|
24 |
|
---|
25 | #include <cmath>
|
---|
26 | #include "math_exception.hpp"
|
---|
27 |
|
---|
28 | namespace math {
|
---|
29 | namespace ublas {
|
---|
30 |
|
---|
31 | /// Definition of a matrix using double precision
|
---|
32 | typedef boost::numeric::ublas::matrix<double> Matrix;
|
---|
33 | /// Definition of a vector using double precision
|
---|
34 | typedef boost::numeric::ublas::vector<double> Vector;
|
---|
35 | /// Definition of empty vector using double precision
|
---|
36 | typedef boost::numeric::ublas::zero_vector<double> ZeroVector;
|
---|
37 | /// Definition of empty matrix using double precision
|
---|
38 | typedef boost::numeric::ublas::zero_matrix<double> ZeroMatrix;
|
---|
39 |
|
---|
40 | /// Multiplication of two matrices.
|
---|
41 | ///
|
---|
42 | /// @tparam T @todo Documentation
|
---|
43 | /// @param m1 ublas matrix
|
---|
44 | /// @param m2 ublas matrix
|
---|
45 | /// @returns ublas matrix
|
---|
46 | template<class T>
|
---|
47 | inline boost::numeric::ublas::matrix<T> operator *(const boost::numeric::ublas::matrix<T> & m1, const boost::numeric::ublas::matrix<T> & m2)
|
---|
48 | {
|
---|
49 | return prod(m1,m2);
|
---|
50 | }
|
---|
51 |
|
---|
52 | /// product of a vector by a matrix
|
---|
53 | ///
|
---|
54 | /// @tparam T @todo Documentation
|
---|
55 | /// @param m ublas matrix
|
---|
56 | /// @param v ublas vector
|
---|
57 | /// @returns ublas vector
|
---|
58 | template<class T>
|
---|
59 | inline boost::numeric::ublas::vector<T> operator *(const boost::numeric::ublas::matrix<T> & m, const boost::numeric::ublas::vector<T> & v)
|
---|
60 | {
|
---|
61 | return prod(m,v);
|
---|
62 | }
|
---|
63 |
|
---|
64 | /// addition of two vectors
|
---|
65 | ///
|
---|
66 | /// @tparam T @todo Documentation
|
---|
67 | /// @param v1 ublas vector
|
---|
68 | /// @param v2 ublas vector
|
---|
69 | /// @returns ublas vector
|
---|
70 | template<class T>
|
---|
71 | inline boost::numeric::ublas::vector<T> operator +(const boost::numeric::ublas::vector<T> & v1, const boost::numeric::ublas::vector<T> & v2)
|
---|
72 | {
|
---|
73 | boost::numeric::ublas::vector<T> tmp = v1;
|
---|
74 | return tmp+=v2;
|
---|
75 | }
|
---|
76 |
|
---|
77 | /// subtraction of two vectors
|
---|
78 | ///
|
---|
79 | /// @tparam T @todo Documentation
|
---|
80 | /// @param v1 ublas vector
|
---|
81 | /// @param v2 ublas vector
|
---|
82 | /// @returns ublas vector
|
---|
83 | template<class T>
|
---|
84 | inline boost::numeric::ublas::vector<T> operator -(const boost::numeric::ublas::vector<T> & v1, const boost::numeric::ublas::vector<T> & v2)
|
---|
85 | {
|
---|
86 | boost::numeric::ublas::vector<T> tmp = v1;
|
---|
87 | return tmp-=v2;
|
---|
88 | }
|
---|
89 |
|
---|
90 | /// Transposes a matrix
|
---|
91 | ///
|
---|
92 | /// @tparam T @todo Documentation
|
---|
93 | /// @param m ublas matrix
|
---|
94 | /// @returns ublas matrix
|
---|
95 | template<class T>
|
---|
96 | inline boost::numeric::ublas::matrix<T> Trans(boost::numeric::ublas::matrix<T> &m)
|
---|
97 | {
|
---|
98 | return trans(m);
|
---|
99 | }
|
---|
100 |
|
---|
101 | /// Converts a vector to a matrix.
|
---|
102 | ///
|
---|
103 | /// @param v ublas vector
|
---|
104 | /// @returns ublas matrix
|
---|
105 | template <class T>
|
---|
106 | inline boost::numeric::ublas::matrix<T> vector2matrix(const boost::numeric::ublas::vector<T> &v)
|
---|
107 | {
|
---|
108 | boost::numeric::ublas::matrix<T> tmp(v.size(),1);
|
---|
109 | for(size_t i=0;i<v.size();i++) {
|
---|
110 | tmp(i,0)=v[i];
|
---|
111 | }
|
---|
112 | return tmp;
|
---|
113 | }
|
---|
114 |
|
---|
115 | /// compute the norm of a vector
|
---|
116 | ///
|
---|
117 | /// @param v ublas vector
|
---|
118 | /// @returns norm value
|
---|
119 | template <class T>
|
---|
120 | inline double Norm(const boost::numeric::ublas::vector<T> & v){
|
---|
121 | double norm =0;
|
---|
122 | for(typename boost::numeric::ublas::vector<T>::const_iterator I=v.begin(); I != v.end(); ++I) {
|
---|
123 | norm+=(*I)*(*I);
|
---|
124 | }
|
---|
125 | return std::sqrt(norm);
|
---|
126 | }
|
---|
127 |
|
---|
128 | /// term by term multiplication of two vectors
|
---|
129 | ///
|
---|
130 | /// @param v1 : ublas vector
|
---|
131 | /// @param v2 : ublas vector
|
---|
132 | /// @returns ublas vector
|
---|
133 | template <class T>
|
---|
134 | inline boost::numeric::ublas::vector<T> Mult(const boost::numeric::ublas::vector<T> & v1, const boost::numeric::ublas::vector<T> & v2)
|
---|
135 | {
|
---|
136 | if(v1.size()!=v2.size()) {
|
---|
137 | throw math_error("Dot(v1,v2) : vectors must have the same size");
|
---|
138 | }
|
---|
139 | boost::numeric::ublas::vector<T> v(v1.size());
|
---|
140 | for(size_t i=0;i<v1.size();i++) {
|
---|
141 | v[i]=v1[i]*v2[i];
|
---|
142 | }
|
---|
143 | return v;
|
---|
144 | }
|
---|
145 |
|
---|
146 | /// dot product
|
---|
147 | ///
|
---|
148 | /// @param v1 ublas vector
|
---|
149 | /// @param v2 ublas vector
|
---|
150 | /// @returns dot product value
|
---|
151 | template <class T>
|
---|
152 | inline double Dot(const boost::numeric::ublas::vector<T> & v1, const boost::numeric::ublas::vector<T> & v2)
|
---|
153 | {
|
---|
154 | if(v1.size()!=v2.size()) {
|
---|
155 | throw math_error("Dot(v1,v2) : vectors must have the same size");
|
---|
156 | }
|
---|
157 | double dot=0;
|
---|
158 | for(size_t i=0;i<v1.size();i++) {
|
---|
159 | dot+=v1[i]*v2[i];
|
---|
160 | }
|
---|
161 | return dot;
|
---|
162 | }
|
---|
163 |
|
---|
164 | /// matrix inversion using LU decomposition
|
---|
165 | /// @param m ublas matrix
|
---|
166 | /// @returns ublas matrix
|
---|
167 | inline Matrix InvLU(const Matrix &m)
|
---|
168 | throw(math_error)
|
---|
169 | {
|
---|
170 | using namespace boost::numeric::ublas;
|
---|
171 | typedef permutation_matrix<std::size_t> pmatrix;
|
---|
172 |
|
---|
173 | if(m.size1() != m.size2()) throw math_error("Inv(m): matrix must be square");
|
---|
174 |
|
---|
175 | // create a working copy of the input
|
---|
176 | Matrix A(m);
|
---|
177 | // create a permutation matrix for the LU-factorization
|
---|
178 | pmatrix pm(A.size1());
|
---|
179 | // perform LU-factorization
|
---|
180 | int res = lu_factorize(A,pm);
|
---|
181 | if( res != 0 ) throw math_error("Inv(m) : singular matrix");
|
---|
182 | // create identity matrix of "inverse"
|
---|
183 | Matrix inverse = identity_matrix<double>(A.size1());
|
---|
184 | // backsubstitute to get the inverse
|
---|
185 | lu_substitute(A, pm, inverse);
|
---|
186 |
|
---|
187 | return inverse;
|
---|
188 | }
|
---|
189 |
|
---|
190 | ////////////////////////////////////////////////////////////////////////////////
|
---|
191 | ///////////////////////// QR DECOMPOSITION ////////////////////////////////////
|
---|
192 | ///////////////////////////////////////////////////////////////////////////////
|
---|
193 |
|
---|
194 | /// @todo Documentation
|
---|
195 | /// @tparam T @todo Documentation
|
---|
196 | template<class T>
|
---|
197 | bool InvertMatrix (const boost::numeric::ublas::matrix<T>& input, boost::numeric::ublas::matrix<T>& inverse)
|
---|
198 | {
|
---|
199 | using namespace boost::numeric::ublas;
|
---|
200 | typedef permutation_matrix<std::size_t> pmatrix;
|
---|
201 | // create a working copy of the input
|
---|
202 | matrix<T> A(input);
|
---|
203 | // create a permutation matrix for the LU-factorization
|
---|
204 | pmatrix pm(A.size1());
|
---|
205 |
|
---|
206 | // perform LU-factorization
|
---|
207 | int res = lu_factorize(A,pm);
|
---|
208 | if( res != 0 ) return false;
|
---|
209 |
|
---|
210 | // create identity matrix of "inverse"
|
---|
211 | inverse.assign(boost::numeric::ublas::identity_matrix<T>(A.size1()));
|
---|
212 |
|
---|
213 | // backsubstitute to get the inverse
|
---|
214 | lu_substitute(A, pm, inverse);
|
---|
215 |
|
---|
216 | return true;
|
---|
217 | }
|
---|
218 |
|
---|
219 | /// @todo Documentation
|
---|
220 | /// @tparam T @todo Documentation
|
---|
221 | template<class T>
|
---|
222 | void TransposeMultiply(const boost::numeric::ublas::vector<T>& vector,
|
---|
223 | boost::numeric::ublas::matrix<T>& result,
|
---|
224 | size_t size)
|
---|
225 | {
|
---|
226 | result.resize (size,size);
|
---|
227 | result.clear ();
|
---|
228 | for(unsigned int row=0; row< vector.size(); ++row)
|
---|
229 | {
|
---|
230 | for(unsigned int col=0; col < vector.size(); ++col)
|
---|
231 | result(row,col) = vector(col) * vector(row);
|
---|
232 |
|
---|
233 | }
|
---|
234 | }
|
---|
235 |
|
---|
236 | /// @todo Documentation
|
---|
237 | /// @tparam T @todo Documentation
|
---|
238 | template<class T>
|
---|
239 | void HouseholderCornerSubstraction (boost::numeric::ublas::matrix<T>& LeftLarge,
|
---|
240 | const boost::numeric::ublas::matrix<T>& RightSmall)
|
---|
241 | {
|
---|
242 | using namespace boost::numeric::ublas;
|
---|
243 | using namespace std;
|
---|
244 | if(
|
---|
245 | !(
|
---|
246 | (LeftLarge.size1() >= RightSmall.size1())
|
---|
247 | && (LeftLarge.size2() >= RightSmall.size2())
|
---|
248 | )
|
---|
249 | )
|
---|
250 | {
|
---|
251 | cerr << "invalid matrix dimensions" << endl;
|
---|
252 | return;
|
---|
253 | }
|
---|
254 |
|
---|
255 | size_t row_offset = LeftLarge.size2() - RightSmall.size2();
|
---|
256 | size_t col_offset = LeftLarge.size1() - RightSmall.size1();
|
---|
257 |
|
---|
258 | for(unsigned int row = 0; row < RightSmall.size2(); ++row )
|
---|
259 | for(unsigned int col = 0; col < RightSmall.size1(); ++col )
|
---|
260 | LeftLarge(col_offset+col,row_offset+row) -= RightSmall(col,row);
|
---|
261 | }
|
---|
262 |
|
---|
263 | /// @todo Documentation
|
---|
264 | /// @tparam T @todo Documentation
|
---|
265 | template<class T>
|
---|
266 | void QR (const boost::numeric::ublas::matrix<T>& M,
|
---|
267 | boost::numeric::ublas::matrix<T>& Q,
|
---|
268 | boost::numeric::ublas::matrix<T>& R)
|
---|
269 | {
|
---|
270 | using namespace boost::numeric::ublas;
|
---|
271 | using namespace std;
|
---|
272 |
|
---|
273 | if(
|
---|
274 | !(
|
---|
275 | (M.size1() == M.size2())
|
---|
276 | )
|
---|
277 | )
|
---|
278 | {
|
---|
279 | cerr << "invalid matrix dimensions" << endl;
|
---|
280 | return;
|
---|
281 | }
|
---|
282 | size_t size = M.size1();
|
---|
283 |
|
---|
284 | // init Matrices
|
---|
285 | matrix<T> H, HTemp;
|
---|
286 | HTemp = identity_matrix<T>(size);
|
---|
287 | Q = identity_matrix<T>(size);
|
---|
288 | R = M;
|
---|
289 |
|
---|
290 | // find Householder reflection matrices
|
---|
291 | for(unsigned int col = 0; col < size-1; ++col)
|
---|
292 | {
|
---|
293 | // create X vector
|
---|
294 | boost::numeric::ublas::vector<T> RRowView = boost::numeric::ublas::column(R,col);
|
---|
295 | vector_range< boost::numeric::ublas::vector<T> > X2 (RRowView, range (col, size));
|
---|
296 | boost::numeric::ublas::vector<T> X = X2;
|
---|
297 |
|
---|
298 | // X -> U~
|
---|
299 | if(X(0) >= 0)
|
---|
300 | X(0) += norm_2(X);
|
---|
301 | else
|
---|
302 | X(0) += -1*norm_2(X);
|
---|
303 |
|
---|
304 | HTemp.resize(X.size(),X.size(),true);
|
---|
305 |
|
---|
306 | TransposeMultiply(X, HTemp, X.size());
|
---|
307 |
|
---|
308 | // HTemp = the 2UUt part of H
|
---|
309 | HTemp *= ( 2 / inner_prod(X,X) );
|
---|
310 |
|
---|
311 | // H = I - 2UUt
|
---|
312 | H = identity_matrix<T>(size);
|
---|
313 | HouseholderCornerSubstraction(H,HTemp);
|
---|
314 |
|
---|
315 | // add H to Q and R
|
---|
316 | Q = prod(Q,H);
|
---|
317 | R = prod(H,R);
|
---|
318 | }
|
---|
319 | }
|
---|
320 | //////////////////////////////////////////////////////////////////////////////////////////::
|
---|
321 |
|
---|
322 | /// Matrix inversion using QR decomposition
|
---|
323 | /// @param m ublas matrix
|
---|
324 | /// @returns ublas matrix
|
---|
325 | inline Matrix InvQR(const Matrix &m)
|
---|
326 | throw(math_error)
|
---|
327 | {
|
---|
328 | using namespace boost::numeric::ublas;
|
---|
329 |
|
---|
330 | if(m.size1() != m.size2()) throw math_error("Inv(m): matrix must be square");
|
---|
331 | Matrix Q(m), R(m), Rinv(m);
|
---|
332 | QR (m,Q,R);
|
---|
333 | for( int i = 0 ; i < R.size1() ; i++ )
|
---|
334 | for( int j = 0 ; j < R.size2() ; j++ )
|
---|
335 | if( R(i,j) < 1e-10 )
|
---|
336 | R(i,j) = 0;
|
---|
337 | InvertMatrix(R,Rinv);
|
---|
338 | return Rinv*Trans(Q);
|
---|
339 | }
|
---|
340 |
|
---|
341 |
|
---|
342 | /// compute matrix determinant using LU decomposition
|
---|
343 | /// @param m ublas matrix
|
---|
344 | /// @returns ublas matrix
|
---|
345 | inline double DetLU(const Matrix & m)
|
---|
346 | throw(math_error)
|
---|
347 | {
|
---|
348 | using namespace boost::numeric::ublas;
|
---|
349 | typedef permutation_matrix<std::size_t> pmatrix;
|
---|
350 |
|
---|
351 |
|
---|
352 | if(m.size1() != m.size2()) throw math_error("Determinant(m): matrix must be square");
|
---|
353 |
|
---|
354 | // create a working copy of the input
|
---|
355 | Matrix A(m);
|
---|
356 | // create a permutation matrix for the LU-factorization
|
---|
357 | pmatrix pm(m.size1());
|
---|
358 | // perform LU-factorization
|
---|
359 | int res = lu_factorize(A, pm);
|
---|
360 | if( res != 0 ) throw math_error("Determinant(m) : singular matrix");
|
---|
361 | //compute determinant
|
---|
362 | double det = 1.0;
|
---|
363 | for (std::size_t i=0; i < pm.size(); ++i) {
|
---|
364 | if (pm(i) != i)
|
---|
365 | det *= -1.0;
|
---|
366 | det *= A(i,i);
|
---|
367 | }
|
---|
368 | return det;
|
---|
369 | }
|
---|
370 |
|
---|
371 |
|
---|
372 | /// output stream function
|
---|
373 | inline std::ostream& operator << (std::ostream& ostrm, const Matrix & m)
|
---|
374 | {
|
---|
375 | for (size_t i=0; i < m.size1(); i++)
|
---|
376 | {
|
---|
377 | ostrm << '['<<'\t';
|
---|
378 | for (size_t j=0; j < m.size2(); j++)
|
---|
379 | {
|
---|
380 | double x = m(i,j);
|
---|
381 | ostrm << x << '\t';
|
---|
382 | }
|
---|
383 | ostrm << ']'<< std::endl;
|
---|
384 | }
|
---|
385 | return ostrm;
|
---|
386 | }
|
---|
387 |
|
---|
388 | /// output stream function
|
---|
389 | inline std::ostream& operator << (std::ostream& ostrm, const Vector & v)
|
---|
390 | {
|
---|
391 | for (size_t i=0; i < v.size(); i++)
|
---|
392 | {
|
---|
393 | ostrm << '['<<'\t';
|
---|
394 | double x = v(i);
|
---|
395 | ostrm << x << '\t';
|
---|
396 | ostrm << ']'<< std::endl;
|
---|
397 | }
|
---|
398 | return ostrm;
|
---|
399 | }
|
---|
400 |
|
---|
401 | } // namespace ublas
|
---|
402 | } // namespace math
|
---|
403 |
|
---|
404 | #endif // __BOOST_UBLAS__
|
---|