Image Component Library (ICL)
|
00001 /******************************************************************** 00002 ** Image Component Library (ICL) ** 00003 ** ** 00004 ** Copyright (C) 2006-2013 CITEC, University of Bielefeld ** 00005 ** Neuroinformatics Group ** 00006 ** Website: www.iclcv.org and ** 00007 ** http://opensource.cit-ec.de/projects/icl ** 00008 ** ** 00009 ** File : ICLMath/src/ICLMath/DynMatrix.h ** 00010 ** Module : ICLMath ** 00011 ** Authors: Christof Elbrechter ** 00012 ** ** 00013 ** ** 00014 ** GNU LESSER GENERAL PUBLIC LICENSE ** 00015 ** This file may be used under the terms of the GNU Lesser General ** 00016 ** Public License version 3.0 as published by the ** 00017 ** ** 00018 ** Free Software Foundation and appearing in the file LICENSE.GPL ** 00019 ** included in the packaging of this file. Please review the ** 00020 ** following information to ensure the license requirements will ** 00021 ** be met: http://www.gnu.org/licenses/lgpl-3.0.txt ** 00022 ** ** 00023 ** The development of this software was supported by the ** 00024 ** Excellence Cluster EXC 277 Cognitive Interaction Technology. ** 00025 ** The Excellence Cluster EXC 277 is a grant of the Deutsche ** 00026 ** Forschungsgemeinschaft (DFG) in the context of the German ** 00027 ** Excellence Initiative. ** 00028 ** ** 00029 ********************************************************************/ 00030 00031 #pragma once 00032 00033 #include <iterator> 00034 #include <algorithm> 00035 #include <numeric> 00036 #include <functional> 00037 #include <iostream> 00038 #include <vector> 00039 #include <cmath> 00040 00041 #include <ICLUtils/Exception.h> 00042 #include <ICLUtils/Macros.h> 00043 00044 #ifdef HAVE_IPP 00045 #include <ipp.h> 00046 #endif 00047 00048 // Intel Math Kernel Library 00049 #ifdef HAVE_MKL 00050 #include "mkl_cblas.h" 00051 #endif 00052 00053 namespace icl{ 00054 namespace math{ 00055 00057 struct InvalidMatrixDimensionException :public utils::ICLException{ 00058 InvalidMatrixDimensionException(const std::string &msg):utils::ICLException(msg){} 00059 }; 00060 00062 struct IncompatibleMatrixDimensionException :public utils::ICLException{ 00063 IncompatibleMatrixDimensionException(const std::string &msg):utils::ICLException(msg){} 00064 }; 00065 00067 struct InvalidIndexException : public utils::ICLException{ 00068 InvalidIndexException(const std::string &msg):utils::ICLException(msg){} 00069 }; 00070 00072 struct SingularMatrixException : public utils::ICLException{ 00073 SingularMatrixException(const std::string &msg):utils::ICLException(msg){} 00074 }; 00075 00077 00081 template<class T> 00082 struct DynMatrix{ 00083 00085 class DynMatrixColumn; 00088 00089 DynMatrix(const DynMatrixColumn &column); 00090 00092 inline DynMatrix():m_rows(0),m_cols(0),m_data(0),m_ownData(true){} 00093 00095 inline DynMatrix(unsigned int cols,unsigned int rows,const T &initValue=0) throw (InvalidMatrixDimensionException) : 00096 m_rows(rows),m_cols(cols),m_ownData(true){ 00097 if(!dim()) throw InvalidMatrixDimensionException("matrix dimensions must be > 0"); 00098 m_data = new T[cols*rows]; 00099 std::fill(begin(),end(),initValue); 00100 } 00101 00103 00105 inline DynMatrix(unsigned int cols,unsigned int rows, T *data, bool deepCopy=true) throw (InvalidMatrixDimensionException) : 00106 m_rows(rows),m_cols(cols),m_ownData(deepCopy){ 00107 if(!dim()) throw InvalidMatrixDimensionException("matrix dimensions must be > 0"); 00108 if(deepCopy){ 00109 m_data = new T[dim()]; 00110 std::copy(data,data+dim(),begin()); 00111 }else{ 00112 m_data = data; 00113 } 00114 } 00115 00117 inline DynMatrix(unsigned int cols,unsigned int rows,const T *data) throw (InvalidMatrixDimensionException) : 00118 m_rows(rows),m_cols(cols),m_ownData(true){ 00119 if(!dim()) throw InvalidMatrixDimensionException("matrix dimensions must be > 0"); 00120 m_data = new T[dim()]; 00121 std::copy(data,data+dim(),begin()); 00122 } 00123 00125 inline DynMatrix(const DynMatrix &other): 00126 m_rows(other.m_rows),m_cols(other.m_cols),m_data(dim() ? new T[dim()] : 0),m_ownData(true){ 00127 std::copy(other.begin(),other.end(),begin()); 00128 } 00129 00131 00132 inline DynMatrix(const std::string &filename):m_rows(0),m_cols(0),m_data(0),m_ownData(true){ 00133 *this = loadCSV(filename); 00134 } 00135 00137 00141 static DynMatrix<T> loadCSV(const std::string &filename) throw (utils::ICLException); 00142 00144 00145 void saveCSV(const std::string &filename) throw (utils::ICLException); 00146 00148 inline bool isNull() const { return !m_data; } 00149 00151 inline ~DynMatrix(){ 00152 if(m_data && m_ownData) delete [] m_data; 00153 } 00154 00156 00160 inline DynMatrix &operator=(const DynMatrix &other){ 00161 if(!m_data && !other.m_ownData){ 00162 m_data = other.m_data; 00163 m_ownData = false; 00164 m_rows = other.m_rows; 00165 m_cols = other.m_cols; 00166 }else{ 00167 if(dim() != other.dim()){ 00168 delete[] m_data; 00169 m_data = other.dim() ? new T[other.dim()] : 0; 00170 } 00171 m_cols = other.m_cols; 00172 m_rows = other.m_rows; 00173 00174 std::copy(other.begin(),other.end(),begin()); 00175 } 00176 return *this; 00177 } 00178 00180 inline void setBounds(unsigned int cols, unsigned int rows, bool holdContent=false, const T &initializer=0) throw (InvalidMatrixDimensionException){ 00181 if((int)cols == m_cols && (int)rows==m_rows) return; 00182 if(!(cols*rows)) throw InvalidMatrixDimensionException("matrix dimensions must be > 0"); 00183 DynMatrix M(cols,rows,initializer); 00184 if(holdContent){ 00185 unsigned int min_cols = iclMin(cols,(unsigned int)m_cols); 00186 unsigned int min_rows = iclMin(rows,(unsigned int)m_rows); 00187 for(unsigned int i=0;i<min_cols;++i){ 00188 for(unsigned int j=0;j<min_rows;++j){ 00189 M(i,j) = (*this)(i,j); 00190 } 00191 } 00192 } 00193 m_cols = cols; 00194 m_rows = rows; 00195 if(m_data && m_ownData) delete [] m_data; 00196 m_data = M.begin(); 00197 m_ownData = true; 00198 M.set_data(0); 00199 } 00200 00202 inline bool isSimilar(const DynMatrix &other, T tollerance=0.0001) const{ 00203 if(other.cols() != cols() || other.rows() != rows()) return false; 00204 for(unsigned int i=0;i<dim();++i){ 00205 T diff = m_data[i] - other.m_data[i]; 00206 if((diff>0?diff:-diff) > tollerance) return false; 00207 } 00208 return true; 00209 } 00210 00212 inline bool operator==(const DynMatrix &other) const{ 00213 if(other.cols() != cols() || other.rows() != rows()) return false; 00214 for(unsigned int i=0;i<dim();++i){ 00215 if(m_data[i] != other.m_data[i]) return false; 00216 } 00217 return true; 00218 } 00219 00221 inline bool operator!=(const DynMatrix &other) const{ 00222 if(other.cols() != cols() || other.rows() != rows()) return false; 00223 for(unsigned int i=0;i<dim();++i){ 00224 if(m_data[i] != other.m_data[i]) return true; 00225 } 00226 return false; 00227 } 00228 00229 00231 inline DynMatrix operator*(T f) const{ 00232 DynMatrix dst(cols(),rows()); 00233 return mult(f,dst); 00234 } 00235 00237 inline DynMatrix &mult(T f, DynMatrix &dst) const{ 00238 dst.setBounds(cols(),rows()); 00239 std::transform(begin(),end(),dst.begin(),std::bind2nd(std::multiplies<T>(),f)); 00240 return dst; 00241 } 00242 00244 inline DynMatrix &operator*=(T f){ 00245 std::transform(begin(),end(),begin(),std::bind2nd(std::multiplies<T>(),f)); 00246 return *this; 00247 } 00248 00250 inline DynMatrix operator/(T f) const{ 00251 return this->operator*(1/f); 00252 } 00253 00255 inline DynMatrix &operator/=(T f){ 00256 return this->operator*=(1/f); 00257 } 00258 00260 inline DynMatrix &mult(const DynMatrix &m, DynMatrix &dst) const throw (IncompatibleMatrixDimensionException){ 00261 if( cols() != m.rows() ) throw IncompatibleMatrixDimensionException("A*B : cols(A) must be rows(B)"); 00262 dst.setBounds(m.cols(),rows()); 00263 for(unsigned int c=0;c<dst.cols();++c){ 00264 for(unsigned int r=0;r<dst.rows();++r){ 00265 dst(c,r) = std::inner_product(row_begin(r),row_end(r),m.col_begin(c),T(0)); 00266 } 00267 } 00268 return dst; 00269 } 00270 00272 inline DynMatrix &elementwise_mult(const DynMatrix &m, DynMatrix &dst) const throw (IncompatibleMatrixDimensionException){ 00273 if((m.cols() != cols()) || (m.rows() != rows())) throw IncompatibleMatrixDimensionException("A.*B dimension mismatch"); 00274 dst.setBounds(cols(),rows()); 00275 for(unsigned int i=0;i<dim();++i){ 00276 dst[i] = m_data[i] * m[i]; 00277 } 00278 return dst; 00279 } 00280 00282 inline DynMatrix elementwise_mult(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){ 00283 DynMatrix dst(cols(),rows()); 00284 return elementwise_mult(m,dst); 00285 } 00286 00288 inline DynMatrix &elementwise_div(const DynMatrix &m, DynMatrix &dst) const throw (IncompatibleMatrixDimensionException){ 00289 if((m.cols() != cols()) || (m.rows() != rows())) throw IncompatibleMatrixDimensionException("A./B dimension mismatch"); 00290 dst.setBounds(cols(),rows()); 00291 for(int i=0;i<dim();++i){ 00292 dst[i] = m_data[i] / m[i]; 00293 } 00294 return dst; 00295 } 00296 00298 inline DynMatrix elementwise_div(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){ 00299 DynMatrix dst(cols(),rows()); 00300 return elementwise_div(m,dst); 00301 } 00302 00303 00304 00305 00307 inline DynMatrix operator*(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){ 00308 DynMatrix d(m.cols(),rows()); 00309 return mult(m,d); 00310 } 00311 00313 inline DynMatrix &operator*=(const DynMatrix &m) throw (IncompatibleMatrixDimensionException){ 00314 return *this=((*this)*m); 00315 } 00316 00318 inline DynMatrix operator/(const DynMatrix &m) const 00319 throw (IncompatibleMatrixDimensionException, 00320 InvalidMatrixDimensionException, 00321 SingularMatrixException){ 00322 return this->operator*(m.inv()); 00323 } 00324 00326 inline DynMatrix &operator/=(const DynMatrix &m) const 00327 throw (IncompatibleMatrixDimensionException, 00328 InvalidMatrixDimensionException, 00329 SingularMatrixException){ 00330 return *this = this->operator*(m.inv()); 00331 } 00332 00334 inline DynMatrix operator+(const T &t) const{ 00335 DynMatrix d(cols(),rows()); 00336 std::transform(begin(),end(),d.begin(),std::bind2nd(std::plus<T>(),t)); 00337 return d; 00338 } 00339 00341 inline DynMatrix operator-(const T &t) const{ 00342 DynMatrix d(cols(),rows()); 00343 std::transform(begin(),end(),d.begin(),std::bind2nd(std::minus<T>(),t)); 00344 return d; 00345 } 00346 00348 inline DynMatrix &operator+=(const T &t){ 00349 std::transform(begin(),end(),begin(),std::bind2nd(std::plus<T>(),t)); 00350 return *this; 00351 } 00352 00354 inline DynMatrix &operator-=(const T &t){ 00355 std::transform(begin(),end(),begin(),std::bind2nd(std::minus<T>(),t)); 00356 return *this; 00357 } 00358 00360 inline DynMatrix operator+(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){ 00361 if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)"); 00362 DynMatrix d(cols(),rows()); 00363 std::transform(begin(),end(),m.begin(),d.begin(),std::plus<T>()); 00364 return d; 00365 } 00366 00368 inline DynMatrix operator-(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){ 00369 if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)"); 00370 DynMatrix d(cols(),rows()); 00371 std::transform(begin(),end(),m.begin(),d.begin(),std::minus<T>()); 00372 return d; 00373 } 00374 00376 inline DynMatrix &operator+=(const DynMatrix &m) throw (IncompatibleMatrixDimensionException){ 00377 if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)"); 00378 std::transform(begin(),end(),m.begin(),begin(),std::plus<T>()); 00379 return *this; 00380 } 00381 00383 inline DynMatrix &operator-=(const DynMatrix &m) throw (IncompatibleMatrixDimensionException){ 00384 if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)"); 00385 std::transform(begin(),end(),m.begin(),begin(),std::minus<T>()); 00386 return *this; 00387 } 00388 00390 inline T &operator()(unsigned int col,unsigned int row){ 00391 #ifdef DYN_MATRIX_INDEX_CHECK 00392 if((int)col >= m_cols || (int)row >= m_rows) ERROR_LOG("access to "<<m_cols<<'x'<<m_rows<<"-matrix index (" << col << "," << row << ")"); 00393 #endif 00394 return m_data[col+cols()*row]; 00395 } 00396 00398 inline const T &operator() (unsigned int col,unsigned int row) const{ 00399 #ifdef DYN_MATRIX_INDEX_CHECK 00400 if((int)col >= m_cols || (int)row >= m_rows) ERROR_LOG("access to "<<m_cols<<'x'<<m_rows<<"-matrix index (" << col << "," << row << ")"); 00401 #endif 00402 return m_data[col+cols()*row]; 00403 } 00404 00406 inline T &at(unsigned int col,unsigned int row) throw (InvalidIndexException){ 00407 if(col>=cols() || row >= rows()) throw InvalidIndexException("row or col index too large"); 00408 return m_data[col+cols()*row]; 00409 } 00410 00412 inline const T &at(unsigned int col,unsigned int row) const throw (InvalidIndexException){ 00413 return const_cast<DynMatrix*>(this)->at(col,row); 00414 } 00415 00416 00417 00419 inline T &operator[](unsigned int idx) { 00420 idx_check(idx); 00421 if(idx >= dim()) ERROR_LOG("access to "<<m_cols<<'x'<<m_rows<<"-matrix index [" << idx<< "]"); 00422 00423 return m_data[idx]; 00424 00425 } 00426 00427 00429 inline const T &operator[](unsigned int idx) const { 00430 idx_check(idx); 00431 return m_data[idx]; 00432 } 00433 00435 inline T norm(double l=2) const{ 00436 double accu = 0; 00437 for(unsigned int i=0;i<dim();++i){ 00438 accu += ::pow(double(m_data[i]),l); 00439 } 00440 return ::pow(double(accu),1.0/l); 00441 } 00442 00444 private: 00445 static T diff_power_two(const T&a, const T&b){ 00446 T d = a-b; 00447 return d*d; 00448 } 00449 public: 00452 00453 inline T sqrDistanceTo(const DynMatrix &other) const throw (InvalidMatrixDimensionException){ 00454 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::sqrDistanceTo: dimension missmatch")); 00455 return std::inner_product(begin(),end(),other.begin(),T(0), std::plus<T>(), diff_power_two); 00456 } 00457 00459 inline T distanceTo(const DynMatrix &other) const throw (InvalidMatrixDimensionException){ 00460 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::distanceTo: dimension missmatch")); 00461 return ::sqrt( distanceTo(other) ); 00462 } 00463 00464 00466 typedef T* iterator; 00467 00469 typedef const T* const_iterator; 00470 00472 typedef T* row_iterator; 00473 00475 typedef const T* const_row_iterator; 00476 00478 unsigned int rows() const { return m_rows; } 00479 00481 unsigned int cols() const { return m_cols; } 00482 00484 T *data() { return m_data; } 00485 00487 const T *data() const { return m_data; } 00488 00490 unsigned int dim() const { return m_rows*m_cols; } 00491 00493 int stride0() const { return sizeof(T) * dim(); } 00494 00496 int stride1() const { return sizeof(T) * cols(); } 00497 00499 int stride2() const { return sizeof(T); } 00500 00502 struct col_iterator : public std::iterator<std::random_access_iterator_tag,T>{ 00503 typedef unsigned int difference_type; 00504 mutable T *p; 00505 unsigned int stride; 00506 inline col_iterator(T *col_begin,unsigned int stride):p(col_begin),stride(stride){} 00507 00508 00510 inline col_iterator &operator++(){ 00511 p+=stride; 00512 return *this; 00513 } 00515 inline const col_iterator &operator++() const{ 00516 p+=stride; 00517 return *this; 00518 } 00520 inline col_iterator operator++(int){ 00521 col_iterator tmp = *this; 00522 ++(*this); 00523 return tmp; 00524 } 00526 inline const col_iterator operator++(int) const{ 00527 col_iterator tmp = *this; 00528 ++(*this); 00529 return tmp; 00530 } 00531 00533 inline col_iterator &operator--(){ 00534 p-=stride; 00535 return *this; 00536 } 00537 00539 inline const col_iterator &operator--() const{ 00540 p-=stride; 00541 return *this; 00542 } 00543 00545 inline col_iterator operator--(int){ 00546 col_iterator tmp = *this; 00547 --(*this); 00548 return tmp; 00549 } 00550 00552 inline const col_iterator operator--(int) const{ 00553 col_iterator tmp = *this; 00554 --(*this); 00555 return tmp; 00556 } 00557 00559 inline col_iterator &operator+=(difference_type n){ 00560 p += n * stride; 00561 return *this; 00562 } 00563 00565 inline const col_iterator &operator+=(difference_type n) const{ 00566 p += n * stride; 00567 return *this; 00568 } 00569 00570 00572 inline col_iterator &operator-=(difference_type n){ 00573 p -= n * stride; 00574 return *this; 00575 } 00576 00578 inline const col_iterator &operator-=(difference_type n) const{ 00579 p -= n * stride; 00580 return *this; 00581 } 00582 00583 00585 inline col_iterator operator+(difference_type n) { 00586 col_iterator tmp = *this; 00587 tmp+=n; 00588 return tmp; 00589 } 00590 00592 inline const col_iterator operator+(difference_type n) const{ 00593 col_iterator tmp = *this; 00594 tmp+=n; 00595 return tmp; 00596 } 00597 00599 inline col_iterator operator-(difference_type n) { 00600 col_iterator tmp = *this; 00601 tmp-=n; 00602 return tmp; 00603 } 00604 00606 inline const col_iterator operator-(difference_type n) const { 00607 col_iterator tmp = *this; 00608 tmp-=n; 00609 return tmp; 00610 } 00611 00612 inline difference_type operator-(const col_iterator &other) const{ 00613 return (p-other.p)/stride; 00614 } 00615 00616 00618 inline T &operator*(){ 00619 return *p; 00620 } 00621 00623 inline T operator*() const{ 00624 return *p; 00625 } 00626 00628 inline bool operator==(const col_iterator &i) const{ return p == i.p; } 00629 00631 inline bool operator!=(const col_iterator &i) const{ return p != i.p; } 00632 00634 inline bool operator<(const col_iterator &i) const{ return p < i.p; } 00635 00637 inline bool operator<=(const col_iterator &i) const{ return p <= i.p; } 00638 00640 inline bool operator>=(const col_iterator &i) const{ return p >= i.p; } 00641 00643 inline bool operator>(const col_iterator &i) const{ return p > i.p; } 00644 }; 00645 00647 typedef const col_iterator const_col_iterator; 00648 00650 class DynMatrixColumn{ 00651 public: 00652 #ifdef DYN_MATRIX_INDEX_CHECK 00653 #define DYN_MATRIX_COLUMN_CHECK(C,E) if(C) ERROR_LOG(E) 00654 #else 00655 #define DYN_MATRIX_COLUMN_CHECK(C,E) 00656 #endif 00657 00658 DynMatrix<T> *matrix; 00659 00661 unsigned int column; 00662 00664 inline DynMatrixColumn(const DynMatrix<T> *matrix, unsigned int column): 00665 matrix(const_cast<DynMatrix<T>*>(matrix)),column(column){ 00666 DYN_MATRIX_COLUMN_CHECK(column >= matrix->cols(),"invalid column index"); 00667 } 00668 00670 inline DynMatrixColumn(const DynMatrix<T> &matrix): 00671 matrix(const_cast<DynMatrix<T>*>(&matrix)),column(0){ 00672 DYN_MATRIX_COLUMN_CHECK(matrix->cols() != 1,"source matrix must have exactly ONE column"); 00673 } 00675 inline DynMatrixColumn(const DynMatrixColumn &c): 00676 matrix(c.matrix),column(c.column){} 00677 00679 inline col_iterator begin() { return matrix->col_begin(column); } 00680 00682 inline col_iterator end() { return matrix->col_end(column); } 00683 00685 inline const col_iterator begin() const { return matrix->col_begin(column); } 00686 00688 inline const col_iterator end() const { return matrix->col_end(column); } 00689 00691 inline unsigned int dim() const { return matrix->rows(); } 00692 00694 inline DynMatrixColumn &operator=(const DynMatrixColumn &c){ 00695 DYN_MATRIX_COLUMN_CHECK(dim() != c.dim(),"dimension missmatch"); 00696 std::copy(c.begin(),c.end(),begin()); 00697 return *this; 00698 } 00699 00701 inline DynMatrixColumn &operator=(const DynMatrix &src){ 00702 DYN_MATRIX_COLUMN_CHECK(dim() != src.dim(),"dimension missmatch"); 00703 std::copy(src.begin(),src.end(),begin()); 00704 return *this; 00705 } 00706 00708 inline DynMatrixColumn &operator+=(const DynMatrixColumn &c){ 00709 DYN_MATRIX_COLUMN_CHECK(dim() != c.dim(),"dimension missmatch"); 00710 std::transform(c.begin(),c.end(),begin(),begin(),std::plus<T>()); 00711 return *this; 00712 } 00713 00715 inline DynMatrixColumn &operator-=(const DynMatrixColumn &c){ 00716 DYN_MATRIX_COLUMN_CHECK(dim() != c.dim(),"dimension missmatch"); 00717 std::transform(c.begin(),c.end(),begin(),begin(),std::minus<T>()); 00718 return *this; 00719 } 00720 00722 inline DynMatrixColumn &operator+=(const DynMatrix &m){ 00723 DYN_MATRIX_COLUMN_CHECK(dim() != m.dim(),"dimension missmatch"); 00724 std::transform(m.begin(),m.end(),begin(),begin(),std::plus<T>()); 00725 return *this; 00726 } 00728 inline DynMatrixColumn &operator-=(const DynMatrix &m){ 00729 DYN_MATRIX_COLUMN_CHECK(dim() != m.dim(),"dimension missmatch"); 00730 std::transform(m.begin(),m.end(),begin(),begin(),std::minus<T>()); 00731 return *this; 00732 } 00733 00735 inline DynMatrixColumn &operator*=(const T&val){ 00736 std::for_each(begin(),end(),std::bind2nd(std::multiplies<T>(),val)); 00737 return *this; 00738 } 00740 inline DynMatrixColumn &operator/=(const T&val){ 00741 std::for_each(begin(),end(),std::bind2nd(std::divides<T>(),val)); 00742 return *this; 00743 } 00744 00745 }; 00746 00747 inline DynMatrix &operator=(const DynMatrixColumn &col){ 00748 DYN_MATRIX_COLUMN_CHECK(dim() != col.dim(),"dimension missmatch"); 00749 std::copy(col.begin(),col.end(),begin()); 00750 return *this; 00751 } 00752 00753 #undef DYN_MATRIX_COLUMN_CHECK 00754 00755 00756 00758 inline iterator begin() { return m_data; } 00759 00761 inline iterator end() { return m_data+dim(); } 00762 00764 inline const_iterator begin() const { return m_data; } 00765 00767 inline const_iterator end() const { return m_data+dim(); } 00768 00770 inline col_iterator col_begin(unsigned int col) { 00771 col_check(col); 00772 return col_iterator(m_data+col,cols()); 00773 } 00774 00776 inline col_iterator col_end(unsigned int col) { 00777 col_check(col); 00778 return col_iterator(m_data+col+dim(),cols()); 00779 } 00780 00782 inline const_col_iterator col_begin(unsigned int col) const { 00783 col_check(col); 00784 return col_iterator(m_data+col,cols()); 00785 } 00786 00788 inline const_col_iterator col_end(unsigned int col) const { 00789 col_check(col); 00790 return col_iterator(m_data+col+dim(),cols()); 00791 } 00792 00794 inline row_iterator row_begin(unsigned int row) { 00795 row_check(row); 00796 return m_data+row*cols(); 00797 } 00798 00800 inline row_iterator row_end(unsigned int row) { 00801 row_check(row); 00802 return m_data+(row+1)*cols(); 00803 } 00804 00806 inline const_row_iterator row_begin(unsigned int row) const { 00807 row_check(row); 00808 return m_data+row*cols(); 00809 } 00810 00812 inline const_row_iterator row_end(unsigned int row) const { 00813 row_check(row); 00814 return m_data+(row+1)*cols(); 00815 } 00816 00818 inline DynMatrix row(int row){ 00819 row_check(row); 00820 return DynMatrix(m_cols,1,row_begin(row),false); 00821 } 00822 00824 inline const DynMatrix row(int row) const{ 00825 row_check(row); 00826 return DynMatrix(m_cols,1,const_cast<T*>(row_begin(row)),false); 00827 } 00828 00830 inline DynMatrixColumn col(int col){ 00831 return DynMatrixColumn(this,col); 00832 } 00833 00834 inline const DynMatrixColumn col(int col) const{ 00835 return DynMatrixColumn(this,col); 00836 } 00837 00839 void decompose_QR(DynMatrix &Q, DynMatrix &R) const; 00840 00842 void decompose_RQ(DynMatrix &R, DynMatrix &Q) const; 00843 00845 00847 void decompose_LU(DynMatrix &L, DynMatrix &U, T zeroThreshold=1E-16) const; 00848 00850 DynMatrix solve_upper_triangular(const DynMatrix &b) const throw(InvalidMatrixDimensionException); 00851 00853 DynMatrix solve_lower_triangular(const DynMatrix &b) const throw(InvalidMatrixDimensionException); 00854 00856 00903 DynMatrix solve(const DynMatrix &b, const std::string &method="lu",T zeroThreshold=1E-16) 00904 throw(InvalidMatrixDimensionException, utils::ICLException, SingularMatrixException); 00905 00906 00908 DynMatrix inv() const throw (InvalidMatrixDimensionException,SingularMatrixException); 00909 00911 00924 void eigen(DynMatrix &eigenvectors, DynMatrix &eigenvalues) const throw(InvalidMatrixDimensionException, utils::ICLException); 00925 00927 00934 void svd(DynMatrix &U, DynMatrix &S, DynMatrix &V) const throw (utils::ICLException); 00935 00937 00961 DynMatrix pinv(bool useSVD=false, T zeroThreshold=1E-16) const 00962 throw (InvalidMatrixDimensionException,SingularMatrixException,utils::ICLException); 00963 00965 00971 DynMatrix big_matrix_pinv(T zeroThreshold=1E-16) const 00972 throw (InvalidMatrixDimensionException,SingularMatrixException,utils::ICLException); 00973 00974 #ifdef HAVE_MKL 00975 typedef void(*GESDD)(const char*, const int*, const int*, T*, const int*, T*, T*, const int*, T*, const int*, T*, const int*, int*, int*); 00976 typedef void(*CBLAS_GEMM)(CBLAS_ORDER,CBLAS_TRANSPOSE,CBLAS_TRANSPOSE,int,int,int,T,const T*,int,const T*,int,T,T*,int); 00977 DynMatrix big_matrix_pinv(T zeroThreshold, GESDD gesdd, CBLAS_GEMM cblas_gemm) const 00978 throw (InvalidMatrixDimensionException,SingularMatrixException,utils::ICLException); 00979 #endif 00980 00982 T det() const throw (InvalidMatrixDimensionException); 00983 00985 inline DynMatrix transp() const{ 00986 DynMatrix d(rows(),cols()); 00987 for(unsigned int x=0;x<cols();++x){ 00988 for(unsigned int y=0;y<rows();++y){ 00989 d(y,x) = (*this)(x,y); 00990 } 00991 } 00992 return d; 00993 } 00994 00996 00997 inline const DynMatrix<T> shallowTransposed() const{ 00998 return DynMatrix<T>(m_rows,m_cols,const_cast<T*>(m_data),false); 00999 } 01000 01002 const DynMatrix<T> shallowTransposed() { 01003 return DynMatrix<T>(m_rows,m_cols,const_cast<T*>(m_data),false); 01004 } 01005 01007 01016 inline void reshape(int newCols, int newRows) throw (InvalidMatrixDimensionException){ 01017 if((cols() * rows()) != (newCols * newRows)){ 01018 throw InvalidMatrixDimensionException("DynMatrix<T>::reshape: source dimension and destination dimension differs!"); 01019 } 01020 m_cols = newCols; 01021 m_rows = newRows; 01022 } 01023 01025 01026 T element_wise_inner_product(const DynMatrix<T> &other) const { 01027 return std::inner_product(begin(),end(),other.begin(),T(0)); 01028 } 01029 01030 01032 01035 DynMatrix<T> dot(const DynMatrix<T> &M) const throw(InvalidMatrixDimensionException){ 01036 return this->transp() * M; 01037 } 01038 01039 01041 DynMatrix<T> diag() const{ 01042 ICLASSERT_RETURN_VAL(cols()==rows(),DynMatrix<T>()); 01043 DynMatrix<T> d(1,rows()); 01044 for(int i=0;i<rows();++i){ 01045 d[i] = (*this)(i,i); 01046 } 01047 return d; 01048 } 01049 01051 T trace() const{ 01052 ICLASSERT_RETURN_VAL(cols()==rows(),0); 01053 double accu = 0; 01054 for(unsigned int i=0;i<dim();i+=cols()+1){ 01055 accu += m_data[i]; 01056 } 01057 return accu; 01058 } 01059 01061 static DynMatrix<T> cross(const DynMatrix<T> &x, const DynMatrix<T> &y) throw(InvalidMatrixDimensionException){ 01062 if(x.cols()==1 && y.cols()==1 && x.rows()==3 && y.rows()==3){ 01063 DynMatrix<T> r(1,x.rows()); 01064 r(0,0) = x(0,1)*y(0,2)-x(0,2)*y(0,1); 01065 r(0,1) = x(0,2)*y(0,0)-x(0,0)*y(0,2); 01066 r(0,2) = x(0,0)*y(0,1)-x(0,1)*y(0,0); 01067 return r; 01068 }else{ 01069 ICLASSERT_RETURN_VAL(x.rows() == 3 && y.rows() == 3,DynMatrix<T>()); 01070 return DynMatrix<T>(); 01071 } 01072 } 01073 01075 T cond(const double p=2) const { 01076 if(cols() == 3 && rows() == 3){ 01077 DynMatrix<T> M_inv = (*this).inv(); 01078 return (*this).norm(p) * M_inv.norm(p); 01079 } else { 01080 DynMatrix<T> U,S,V; 01081 (*this).svd(U,S,V); 01082 if(S[S.rows()-1]){ 01083 return S[0]/S[S.rows()-1]; 01084 } else { 01085 return S[0]; 01086 } 01087 } 01088 } 01089 01091 inline T *set_data(T *newData){ 01092 T *old_data = m_data; 01093 m_data = newData; 01094 return old_data; 01095 } 01096 01098 static inline DynMatrix id(unsigned int dim) { 01099 DynMatrix M(dim,dim); 01100 for(unsigned int i=0;i<dim;++i) M(i,i) = 1; 01101 return M; 01102 } 01103 01104 private: 01105 inline void row_check(unsigned int row) const{ 01106 #ifdef DYN_MATRIX_INDEX_CHECK 01107 if((int)row >= m_rows) ERROR_LOG("access to row index " << row << " on a "<<m_cols<<'x'<<m_rows<<"-matrix"); 01108 #else 01109 (void)row; 01110 #endif 01111 } 01112 inline void col_check(unsigned int col) const{ 01113 #ifdef DYN_MATRIX_INDEX_CHECK 01114 if((int)col >= m_cols) ERROR_LOG("access to column index " << col << " on a "<<m_cols<<'x'<<m_rows<<"-matrix"); 01115 #else 01116 (void)col; 01117 #endif 01118 } 01119 inline void idx_check(unsigned int col, unsigned int row) const{ 01120 col_check(col); 01121 row_check(row); 01122 } 01123 01124 inline void idx_check(unsigned int idx) const{ 01125 #ifdef DYN_MATRIX_INDEX_CHECK 01126 if(idx >= dim()) ERROR_LOG("access to linear index " << idx << " on a "<<m_cols<<'x'<<m_rows<<"-matrix"); 01127 #else 01128 (void)idx; 01129 #endif 01130 } 01131 01132 int m_rows; 01133 int m_cols; 01134 T *m_data; 01135 bool m_ownData; 01136 }; 01137 01139 01140 template<class T> 01141 DynMatrix<T>::DynMatrix(const DynMatrix<T>::DynMatrixColumn &column): 01142 m_rows(column.dim()),m_cols(1),m_data(new T[column.dim()]),m_ownData(true){ 01143 std::copy(column.begin(),column.end(),begin()); 01144 } 01147 01148 template<class T> 01149 std::ostream &operator<<(std::ostream &s,const DynMatrix<T> &m); 01150 01152 template<class T> 01153 std::istream &operator>>(std::istream &s,DynMatrix<T> &m); 01154 01155 01156 #ifdef HAVE_IPP 01157 01158 template<> 01159 inline float DynMatrix<float>::sqrDistanceTo(const DynMatrix<float> &other) const throw (InvalidMatrixDimensionException){ 01160 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::sqrDistanceTo: dimension missmatch")); 01161 float norm = 0 ; 01162 ippsNormDiff_L2_32f(begin(), other.begin(), dim(), &norm); 01163 return norm*norm; 01164 } 01165 01166 template<> 01167 inline double DynMatrix<double>::sqrDistanceTo(const DynMatrix<double> &other) const throw (InvalidMatrixDimensionException){ 01168 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::sqrDistanceTo: dimension missmatch")); 01169 double norm = 0 ; 01170 ippsNormDiff_L2_64f(begin(), other.begin(), dim(), &norm); 01171 return norm*norm; 01172 } 01173 01174 template<> 01175 inline float DynMatrix<float>::distanceTo(const DynMatrix<float> &other) const throw (InvalidMatrixDimensionException){ 01176 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::distanceTo: dimension missmatch")); 01177 float norm = 0 ; 01178 ippsNormDiff_L2_32f(begin(), other.begin(), dim(), &norm); 01179 return norm; 01180 } 01181 01182 template<> 01183 inline double DynMatrix<double>::distanceTo(const DynMatrix<double> &other) const throw (InvalidMatrixDimensionException){ 01184 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::distanceTo: dimension missmatch")); 01185 double norm = 0 ; 01186 ippsNormDiff_L2_64f(begin(), other.begin(), dim(), &norm); 01187 return norm; 01188 } 01189 01190 #define DYN_MATRIX_MULT_SPECIALIZE(IPPT) \ 01191 template<> \ 01192 inline DynMatrix<Ipp##IPPT> &DynMatrix<Ipp##IPPT>::mult( \ 01193 const DynMatrix<Ipp##IPPT> &m, DynMatrix<Ipp##IPPT> &dst) const \ 01194 throw (IncompatibleMatrixDimensionException){ \ 01195 if(cols() != m.rows() ) throw IncompatibleMatrixDimensionException("A*B : cols(A) must be row(B)"); \ 01196 dst.setBounds(m.cols(),rows()); \ 01197 ippmMul_mm_##IPPT(data(),sizeof(Ipp##IPPT)*cols(),sizeof(Ipp##IPPT),cols(),rows(), \ 01198 m.data(),sizeof(Ipp##IPPT)*m.cols(),sizeof(Ipp##IPPT),m.cols(),m.rows(), \ 01199 dst.data(),m.cols()*sizeof(Ipp##IPPT),sizeof(Ipp##IPPT)); \ 01200 return dst; \ 01201 } 01202 01203 DYN_MATRIX_MULT_SPECIALIZE(32f) 01204 DYN_MATRIX_MULT_SPECIALIZE(64f) 01205 #undef DYN_MATRIX_MULT_SPECIALIZE 01206 01207 01208 #define DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE(IPPT) \ 01209 template<> \ 01210 inline DynMatrix<Ipp##IPPT> &DynMatrix<Ipp##IPPT>::elementwise_div( \ 01211 const DynMatrix<Ipp##IPPT> &m, DynMatrix<Ipp##IPPT> &dst) const \ 01212 throw (IncompatibleMatrixDimensionException){ \ 01213 if((m.cols() != cols()) || (m.rows() != rows())){ \ 01214 throw IncompatibleMatrixDimensionException("A./B dimension mismatch"); \ 01215 } \ 01216 dst.setBounds(cols(),rows()); \ 01217 ippsDiv_##IPPT(data(),m.data(),dst.data(),dim()); \ 01218 return dst; \ 01219 } 01220 01221 DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE(32f) 01222 DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE(64f) 01223 01224 #undef DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE 01225 01226 01227 01228 01229 01230 01231 #define DYN_MATRIX_MULT_BY_CONSTANT(IPPT) \ 01232 template<> \ 01233 inline DynMatrix<Ipp##IPPT> &DynMatrix<Ipp##IPPT>::mult( \ 01234 Ipp##IPPT f, DynMatrix<Ipp##IPPT> &dst) const{ \ 01235 dst.setBounds(cols(),rows()); \ 01236 ippsMulC_##IPPT(data(),f, dst.data(),dim()); \ 01237 return dst; \ 01238 } 01239 01240 DYN_MATRIX_MULT_BY_CONSTANT(32f) 01241 DYN_MATRIX_MULT_BY_CONSTANT(64f) 01242 01243 #undef DYN_MATRIX_MULT_BY_CONSTANT 01244 01245 #define DYN_MATRIX_NORM_SPECIALZE(T,IPPT) \ 01246 template<> \ 01247 inline T DynMatrix<T> ::norm(double l) const{ \ 01248 if(l==1){ \ 01249 T val; \ 01250 ippsNorm_L1_##IPPT(m_data,dim(),&val); \ 01251 return val; \ 01252 }else if(l==2){ \ 01253 T val; \ 01254 ippsNorm_L2_##IPPT(m_data,dim(),&val); \ 01255 return val; \ 01256 } \ 01257 double accu = 0; \ 01258 for(unsigned int i=0;i<dim();++i){ \ 01259 accu += ::pow(double(m_data[i]),l); \ 01260 } \ 01261 return ::pow(accu,1.0/l); \ 01262 } 01263 01264 DYN_MATRIX_NORM_SPECIALZE(float,32f) 01265 // DYN_MATRIX_NORM_SPECIALZE(double,64f) 01266 01267 #undef DYN_MATRIX_NORM_SPECIALZE 01268 01271 #endif // HAVE_IPP 01272 01274 01275 template<class T> 01276 inline DynMatrix<T> operator,(const DynMatrix<T> &left, const DynMatrix<T> &right){ 01277 int w = left.cols() + right.cols(); 01278 int h = iclMax(left.rows(),right.rows()); 01279 DynMatrix<T> result(w,h,float(0)); 01280 for(unsigned int y=0;y<left.rows();++y){ 01281 std::copy(left.row_begin(y), left.row_end(y), result.row_begin(y)); 01282 } 01283 for(unsigned int y=0;y<right.rows();++y){ 01284 std::copy(right.row_begin(y), right.row_end(y), result.row_begin(y) + left.cols()); 01285 } 01286 return result; 01287 } 01288 01290 01291 template<class T> 01292 inline DynMatrix<T> operator%(const DynMatrix<T> &top, const DynMatrix<T> &bottom){ 01293 int w = iclMax(top.cols(),bottom.cols()); 01294 int h = top.rows() + bottom.rows(); 01295 DynMatrix<T> result(w,h,float(0)); 01296 for(unsigned int y=0;y<top.rows();++y){ 01297 std::copy(top.row_begin(y), top.row_end(y), result.row_begin(y)); 01298 } 01299 for(unsigned int y=0;y<bottom.rows();++y){ 01300 std::copy(bottom.row_begin(y), bottom.row_end(y), result.row_begin(y+top.rows())); 01301 } 01302 return result; 01303 } 01304 } // namespace math 01305 }