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