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.LGPL ** 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 <ICLUtils/Macros.h> 00034 #include <ICLUtils/Exception.h> 00035 00036 #include <iterator> 00037 #include <algorithm> 00038 #include <numeric> 00039 #include <functional> 00040 #include <vector> 00041 #include <cmath> 00042 00043 #ifdef ICL_HAVE_IPP 00044 #include <ipp.h> 00045 #endif 00046 00047 // Intel Math Kernel Library 00048 #ifdef ICL_HAVE_MKL 00049 #include "mkl_cblas.h" 00050 #endif 00051 00052 namespace icl{ 00053 namespace math{ 00054 00056 struct InvalidMatrixDimensionException :public utils::ICLException{ 00057 InvalidMatrixDimensionException(const std::string &msg):utils::ICLException(msg){} 00058 }; 00059 00061 struct IncompatibleMatrixDimensionException :public utils::ICLException{ 00062 IncompatibleMatrixDimensionException(const std::string &msg):utils::ICLException(msg){} 00063 }; 00064 00066 struct InvalidIndexException : public utils::ICLException{ 00067 InvalidIndexException(const std::string &msg):utils::ICLException(msg){} 00068 }; 00069 00071 struct SingularMatrixException : public utils::ICLException{ 00072 SingularMatrixException(const std::string &msg):utils::ICLException(msg){} 00073 }; 00074 00076 00080 template<class T> 00081 struct DynMatrix{ 00082 00084 class DynMatrixColumn; 00087 00088 DynMatrix(const DynMatrixColumn &column); 00089 00091 inline DynMatrix():m_rows(0),m_cols(0),m_data(0),m_ownData(true){} 00092 00094 inline DynMatrix(unsigned int cols,unsigned int rows,const T &initValue=0) throw (InvalidMatrixDimensionException) : 00095 m_rows(rows),m_cols(cols),m_ownData(true){ 00096 if(!dim()) throw InvalidMatrixDimensionException("matrix dimensions must be > 0"); 00097 m_data = new T[cols*rows]; 00098 std::fill(begin(),end(),initValue); 00099 } 00100 00102 00104 inline DynMatrix(unsigned int cols,unsigned int rows, T *data, bool deepCopy=true) throw (InvalidMatrixDimensionException) : 00105 m_rows(rows),m_cols(cols),m_ownData(deepCopy){ 00106 if(!dim()) throw InvalidMatrixDimensionException("matrix dimensions must be > 0"); 00107 if(deepCopy){ 00108 m_data = new T[dim()]; 00109 std::copy(data,data+dim(),begin()); 00110 }else{ 00111 m_data = data; 00112 } 00113 } 00114 00116 inline DynMatrix(unsigned int cols,unsigned int rows,const T *data) throw (InvalidMatrixDimensionException) : 00117 m_rows(rows),m_cols(cols),m_ownData(true){ 00118 if(!dim()) throw InvalidMatrixDimensionException("matrix dimensions must be > 0"); 00119 m_data = new T[dim()]; 00120 std::copy(data,data+dim(),begin()); 00121 } 00122 00124 inline DynMatrix(const DynMatrix &other): 00125 m_rows(other.m_rows),m_cols(other.m_cols),m_data(dim() ? new T[dim()] : 0),m_ownData(true){ 00126 std::copy(other.begin(),other.end(),begin()); 00127 } 00128 00130 00131 inline DynMatrix(const std::string &filename):m_rows(0),m_cols(0),m_data(0),m_ownData(true){ 00132 *this = loadCSV(filename); 00133 } 00134 00136 00140 static DynMatrix<T> loadCSV(const std::string &filename) throw (utils::ICLException); 00141 00143 00144 void saveCSV(const std::string &filename) throw (utils::ICLException); 00145 00147 inline bool isNull() const { return !m_data; } 00148 00150 inline ~DynMatrix(){ 00151 if(m_data && m_ownData) delete [] m_data; 00152 } 00153 00155 00159 inline DynMatrix &operator=(const DynMatrix &other){ 00160 if(!m_data && !other.m_ownData){ 00161 m_data = other.m_data; 00162 m_ownData = false; 00163 m_rows = other.m_rows; 00164 m_cols = other.m_cols; 00165 }else{ 00166 if(dim() != other.dim()){ 00167 delete[] m_data; 00168 m_data = other.dim() ? new T[other.dim()] : 0; 00169 } 00170 m_cols = other.m_cols; 00171 m_rows = other.m_rows; 00172 00173 std::copy(other.begin(),other.end(),begin()); 00174 } 00175 return *this; 00176 } 00177 00179 inline void setBounds(unsigned int cols, unsigned int rows, bool holdContent=false, const T &initializer=0) throw (InvalidMatrixDimensionException){ 00180 if((int)cols == m_cols && (int)rows==m_rows) return; 00181 if(!(cols*rows)) throw InvalidMatrixDimensionException("matrix dimensions must be > 0"); 00182 DynMatrix M(cols,rows,initializer); 00183 if(holdContent){ 00184 unsigned int min_cols = iclMin(cols,(unsigned int)m_cols); 00185 unsigned int min_rows = iclMin(rows,(unsigned int)m_rows); 00186 for(unsigned int i=0;i<min_cols;++i){ 00187 for(unsigned int j=0;j<min_rows;++j){ 00188 M(i,j) = (*this)(i,j); 00189 } 00190 } 00191 } 00192 m_cols = cols; 00193 m_rows = rows; 00194 if(m_data && m_ownData) delete [] m_data; 00195 m_data = M.begin(); 00196 m_ownData = true; 00197 M.set_data(0); 00198 } 00199 00201 inline bool isSimilar(const DynMatrix &other, T tollerance=0.0001) const{ 00202 if(other.cols() != cols() || other.rows() != rows()) return false; 00203 for(unsigned int i=0;i<dim();++i){ 00204 T diff = m_data[i] - other.m_data[i]; 00205 if((diff>0?diff:-diff) > tollerance) return false; 00206 } 00207 return true; 00208 } 00209 00211 inline bool operator==(const DynMatrix &other) const{ 00212 if(other.cols() != cols() || other.rows() != rows()) return false; 00213 for(unsigned int i=0;i<dim();++i){ 00214 if(m_data[i] != other.m_data[i]) return false; 00215 } 00216 return true; 00217 } 00218 00220 inline bool operator!=(const DynMatrix &other) const{ 00221 if(other.cols() != cols() || other.rows() != rows()) return false; 00222 for(unsigned int i=0;i<dim();++i){ 00223 if(m_data[i] != other.m_data[i]) return true; 00224 } 00225 return false; 00226 } 00227 00228 00230 inline DynMatrix operator*(T f) const{ 00231 DynMatrix dst(cols(),rows()); 00232 return mult(f,dst); 00233 } 00234 00236 inline DynMatrix &mult(T f, DynMatrix &dst) const{ 00237 dst.setBounds(cols(),rows()); 00238 std::transform(begin(),end(),dst.begin(),std::bind2nd(std::multiplies<T>(),f)); 00239 return dst; 00240 } 00241 00243 inline DynMatrix &operator*=(T f){ 00244 std::transform(begin(),end(),begin(),std::bind2nd(std::multiplies<T>(),f)); 00245 return *this; 00246 } 00247 00249 inline DynMatrix operator/(T f) const{ 00250 return this->operator*(1/f); 00251 } 00252 00254 inline DynMatrix &operator/=(T f){ 00255 return this->operator*=(1/f); 00256 } 00257 00259 inline DynMatrix &mult(const DynMatrix &m, DynMatrix &dst) const throw (IncompatibleMatrixDimensionException){ 00260 if( cols() != m.rows() ) throw IncompatibleMatrixDimensionException("A*B : cols(A) must be rows(B)"); 00261 dst.setBounds(m.cols(),rows()); 00262 for(unsigned int c=0;c<dst.cols();++c){ 00263 for(unsigned int r=0;r<dst.rows();++r){ 00264 dst(c,r) = std::inner_product(row_begin(r),row_end(r),m.col_begin(c),T(0)); 00265 } 00266 } 00267 return dst; 00268 } 00269 00271 inline DynMatrix &elementwise_mult(const DynMatrix &m, DynMatrix &dst) const throw (IncompatibleMatrixDimensionException){ 00272 if((m.cols() != cols()) || (m.rows() != rows())) throw IncompatibleMatrixDimensionException("A.*B dimension mismatch"); 00273 dst.setBounds(cols(),rows()); 00274 for(unsigned int i=0;i<dim();++i){ 00275 dst[i] = m_data[i] * m[i]; 00276 } 00277 return dst; 00278 } 00279 00281 inline DynMatrix elementwise_mult(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){ 00282 DynMatrix dst(cols(),rows()); 00283 return elementwise_mult(m,dst); 00284 } 00285 00287 inline DynMatrix &elementwise_div(const DynMatrix &m, DynMatrix &dst) const throw (IncompatibleMatrixDimensionException){ 00288 if((m.cols() != cols()) || (m.rows() != rows())) throw IncompatibleMatrixDimensionException("A./B dimension mismatch"); 00289 dst.setBounds(cols(),rows()); 00290 for(int i=0;i<dim();++i){ 00291 dst[i] = m_data[i] / m[i]; 00292 } 00293 return dst; 00294 } 00295 00297 inline DynMatrix elementwise_div(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){ 00298 DynMatrix dst(cols(),rows()); 00299 return elementwise_div(m,dst); 00300 } 00301 00302 00303 00304 00306 inline DynMatrix operator*(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){ 00307 DynMatrix d(m.cols(),rows()); 00308 return mult(m,d); 00309 } 00310 00312 inline DynMatrix &operator*=(const DynMatrix &m) throw (IncompatibleMatrixDimensionException){ 00313 return *this=((*this)*m); 00314 } 00315 00317 inline DynMatrix operator/(const DynMatrix &m) const 00318 throw (IncompatibleMatrixDimensionException, 00319 InvalidMatrixDimensionException, 00320 SingularMatrixException){ 00321 return this->operator*(m.inv()); 00322 } 00323 00325 inline DynMatrix &operator/=(const DynMatrix &m) const 00326 throw (IncompatibleMatrixDimensionException, 00327 InvalidMatrixDimensionException, 00328 SingularMatrixException){ 00329 return *this = this->operator*(m.inv()); 00330 } 00331 00333 inline DynMatrix operator+(const T &t) const{ 00334 DynMatrix d(cols(),rows()); 00335 std::transform(begin(),end(),d.begin(),std::bind2nd(std::plus<T>(),t)); 00336 return d; 00337 } 00338 00340 inline DynMatrix operator-(const T &t) const{ 00341 DynMatrix d(cols(),rows()); 00342 std::transform(begin(),end(),d.begin(),std::bind2nd(std::minus<T>(),t)); 00343 return d; 00344 } 00345 00347 inline DynMatrix &operator+=(const T &t){ 00348 std::transform(begin(),end(),begin(),std::bind2nd(std::plus<T>(),t)); 00349 return *this; 00350 } 00351 00353 inline DynMatrix &operator-=(const T &t){ 00354 std::transform(begin(),end(),begin(),std::bind2nd(std::minus<T>(),t)); 00355 return *this; 00356 } 00357 00359 inline DynMatrix operator+(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){ 00360 if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)"); 00361 DynMatrix d(cols(),rows()); 00362 std::transform(begin(),end(),m.begin(),d.begin(),std::plus<T>()); 00363 return d; 00364 } 00365 00367 inline DynMatrix operator-(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){ 00368 if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)"); 00369 DynMatrix d(cols(),rows()); 00370 std::transform(begin(),end(),m.begin(),d.begin(),std::minus<T>()); 00371 return d; 00372 } 00373 00375 inline DynMatrix &operator+=(const DynMatrix &m) throw (IncompatibleMatrixDimensionException){ 00376 if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)"); 00377 std::transform(begin(),end(),m.begin(),begin(),std::plus<T>()); 00378 return *this; 00379 } 00380 00382 inline DynMatrix &operator-=(const DynMatrix &m) throw (IncompatibleMatrixDimensionException){ 00383 if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)"); 00384 std::transform(begin(),end(),m.begin(),begin(),std::minus<T>()); 00385 return *this; 00386 } 00387 00389 inline T &operator()(unsigned int col,unsigned int row){ 00390 #ifdef DYN_MATRIX_INDEX_CHECK 00391 if((int)col >= m_cols || (int)row >= m_rows) ERROR_LOG("access to "<<m_cols<<'x'<<m_rows<<"-matrix index (" << col << "," << row << ")"); 00392 #endif 00393 return m_data[col+cols()*row]; 00394 } 00395 00397 inline const T &operator() (unsigned int col,unsigned int row) const{ 00398 #ifdef DYN_MATRIX_INDEX_CHECK 00399 if((int)col >= m_cols || (int)row >= m_rows) ERROR_LOG("access to "<<m_cols<<'x'<<m_rows<<"-matrix index (" << col << "," << row << ")"); 00400 #endif 00401 return m_data[col+cols()*row]; 00402 } 00403 00405 inline T &at(unsigned int col,unsigned int row) throw (InvalidIndexException){ 00406 if(col>=cols() || row >= rows()) throw InvalidIndexException("row or col index too large"); 00407 return m_data[col+cols()*row]; 00408 } 00409 00411 inline const T &at(unsigned int col,unsigned int row) const throw (InvalidIndexException){ 00412 return const_cast<DynMatrix*>(this)->at(col,row); 00413 } 00414 00415 00416 00418 inline T &operator[](unsigned int idx) { 00419 idx_check(idx); 00420 if(idx >= dim()) ERROR_LOG("access to "<<m_cols<<'x'<<m_rows<<"-matrix index [" << idx<< "]"); 00421 00422 return m_data[idx]; 00423 00424 } 00425 00426 00428 inline const T &operator[](unsigned int idx) const { 00429 idx_check(idx); 00430 return m_data[idx]; 00431 } 00432 00434 inline T norm(double l=2) const{ 00435 double accu = 0; 00436 for(unsigned int i=0;i<dim();++i){ 00437 accu += ::pow(double(m_data[i]),l); 00438 } 00439 return ::pow(double(accu),1.0/l); 00440 } 00441 00443 private: 00444 static T diff_power_two(const T&a, const T&b){ 00445 T d = a-b; 00446 return d*d; 00447 } 00448 public: 00451 00452 inline T sqrDistanceTo(const DynMatrix &other) const throw (InvalidMatrixDimensionException){ 00453 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::sqrDistanceTo: dimension missmatch")); 00454 return std::inner_product(begin(),end(),other.begin(),T(0), std::plus<T>(), diff_power_two); 00455 } 00456 00458 inline T distanceTo(const DynMatrix &other) const throw (InvalidMatrixDimensionException){ 00459 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::distanceTo: dimension missmatch")); 00460 return ::sqrt( distanceTo(other) ); 00461 } 00462 00463 00465 typedef T* iterator; 00466 00468 typedef const T* const_iterator; 00469 00471 typedef T* row_iterator; 00472 00474 typedef const T* const_row_iterator; 00475 00477 unsigned int rows() const { return m_rows; } 00478 00480 unsigned int cols() const { return m_cols; } 00481 00483 T *data() { return m_data; } 00484 00486 const T *data() const { return m_data; } 00487 00489 unsigned int dim() const { return m_rows*m_cols; } 00490 00492 int stride0() const { return sizeof(T) * dim(); } 00493 00495 int stride1() const { return sizeof(T) * cols(); } 00496 00498 int stride2() const { return sizeof(T); } 00499 00501 struct col_iterator : public std::iterator<std::random_access_iterator_tag,T>{ 00502 typedef unsigned int difference_type; 00503 mutable T *p; 00504 unsigned int stride; 00505 inline col_iterator(T *col_begin,unsigned int stride):p(col_begin),stride(stride){} 00506 00507 00509 inline col_iterator &operator++(){ 00510 p+=stride; 00511 return *this; 00512 } 00514 inline const col_iterator &operator++() const{ 00515 p+=stride; 00516 return *this; 00517 } 00519 inline col_iterator operator++(int){ 00520 col_iterator tmp = *this; 00521 ++(*this); 00522 return tmp; 00523 } 00525 inline const col_iterator operator++(int) const{ 00526 col_iterator tmp = *this; 00527 ++(*this); 00528 return tmp; 00529 } 00530 00532 inline col_iterator &operator--(){ 00533 p-=stride; 00534 return *this; 00535 } 00536 00538 inline const col_iterator &operator--() const{ 00539 p-=stride; 00540 return *this; 00541 } 00542 00544 inline col_iterator operator--(int){ 00545 col_iterator tmp = *this; 00546 --(*this); 00547 return tmp; 00548 } 00549 00551 inline const col_iterator operator--(int) const{ 00552 col_iterator tmp = *this; 00553 --(*this); 00554 return tmp; 00555 } 00556 00558 inline col_iterator &operator+=(difference_type n){ 00559 p += n * stride; 00560 return *this; 00561 } 00562 00564 inline const col_iterator &operator+=(difference_type n) const{ 00565 p += n * stride; 00566 return *this; 00567 } 00568 00569 00571 inline col_iterator &operator-=(difference_type n){ 00572 p -= n * stride; 00573 return *this; 00574 } 00575 00577 inline const col_iterator &operator-=(difference_type n) const{ 00578 p -= n * stride; 00579 return *this; 00580 } 00581 00582 00584 inline col_iterator operator+(difference_type n) { 00585 col_iterator tmp = *this; 00586 tmp+=n; 00587 return tmp; 00588 } 00589 00591 inline const col_iterator operator+(difference_type n) const{ 00592 col_iterator tmp = *this; 00593 tmp+=n; 00594 return tmp; 00595 } 00596 00598 inline col_iterator operator-(difference_type n) { 00599 col_iterator tmp = *this; 00600 tmp-=n; 00601 return tmp; 00602 } 00603 00605 inline const col_iterator operator-(difference_type n) const { 00606 col_iterator tmp = *this; 00607 tmp-=n; 00608 return tmp; 00609 } 00610 00611 inline difference_type operator-(const col_iterator &other) const{ 00612 return (p-other.p)/stride; 00613 } 00614 00615 00617 inline T &operator*(){ 00618 return *p; 00619 } 00620 00622 inline T operator*() const{ 00623 return *p; 00624 } 00625 00627 inline bool operator==(const col_iterator &i) const{ return p == i.p; } 00628 00630 inline bool operator!=(const col_iterator &i) const{ return p != i.p; } 00631 00633 inline bool operator<(const col_iterator &i) const{ return p < i.p; } 00634 00636 inline bool operator<=(const col_iterator &i) const{ return p <= i.p; } 00637 00639 inline bool operator>=(const col_iterator &i) const{ return p >= i.p; } 00640 00642 inline bool operator>(const col_iterator &i) const{ return p > i.p; } 00643 }; 00644 00646 typedef const col_iterator const_col_iterator; 00647 00649 class ICLMath_API DynMatrixColumn{ 00650 public: 00651 #ifdef DYN_MATRIX_INDEX_CHECK 00652 #define DYN_MATRIX_COLUMN_CHECK(C,E) if(C) ERROR_LOG(E) 00653 #else 00654 #define DYN_MATRIX_COLUMN_CHECK(C,E) 00655 #endif 00656 00657 DynMatrix<T> *matrix; 00658 00660 unsigned int column; 00661 00663 inline DynMatrixColumn(const DynMatrix<T> *matrix, unsigned int column): 00664 matrix(const_cast<DynMatrix<T>*>(matrix)),column(column){ 00665 DYN_MATRIX_COLUMN_CHECK(column >= matrix->cols(),"invalid column index"); 00666 } 00667 00669 inline DynMatrixColumn(const DynMatrix<T> &matrix): 00670 matrix(const_cast<DynMatrix<T>*>(&matrix)),column(0){ 00671 DYN_MATRIX_COLUMN_CHECK(matrix->cols() != 1,"source matrix must have exactly ONE column"); 00672 } 00674 inline DynMatrixColumn(const DynMatrixColumn &c): 00675 matrix(c.matrix),column(c.column){} 00676 00678 inline col_iterator begin() { return matrix->col_begin(column); } 00679 00681 inline col_iterator end() { return matrix->col_end(column); } 00682 00684 inline const col_iterator begin() const { return matrix->col_begin(column); } 00685 00687 inline const col_iterator end() const { return matrix->col_end(column); } 00688 00690 inline unsigned int dim() const { return matrix->rows(); } 00691 00693 inline DynMatrixColumn &operator=(const DynMatrixColumn &c){ 00694 DYN_MATRIX_COLUMN_CHECK(dim() != c.dim(),"dimension missmatch"); 00695 std::copy(c.begin(),c.end(),begin()); 00696 return *this; 00697 } 00698 00700 inline DynMatrixColumn &operator=(const DynMatrix &src){ 00701 DYN_MATRIX_COLUMN_CHECK(dim() != src.dim(),"dimension missmatch"); 00702 std::copy(src.begin(),src.end(),begin()); 00703 return *this; 00704 } 00705 00707 inline DynMatrixColumn &operator+=(const DynMatrixColumn &c){ 00708 DYN_MATRIX_COLUMN_CHECK(dim() != c.dim(),"dimension missmatch"); 00709 std::transform(c.begin(),c.end(),begin(),begin(),std::plus<T>()); 00710 return *this; 00711 } 00712 00714 inline DynMatrixColumn &operator-=(const DynMatrixColumn &c){ 00715 DYN_MATRIX_COLUMN_CHECK(dim() != c.dim(),"dimension missmatch"); 00716 std::transform(c.begin(),c.end(),begin(),begin(),std::minus<T>()); 00717 return *this; 00718 } 00719 00721 inline DynMatrixColumn &operator+=(const DynMatrix &m){ 00722 DYN_MATRIX_COLUMN_CHECK(dim() != m.dim(),"dimension missmatch"); 00723 std::transform(m.begin(),m.end(),begin(),begin(),std::plus<T>()); 00724 return *this; 00725 } 00727 inline DynMatrixColumn &operator-=(const DynMatrix &m){ 00728 DYN_MATRIX_COLUMN_CHECK(dim() != m.dim(),"dimension missmatch"); 00729 std::transform(m.begin(),m.end(),begin(),begin(),std::minus<T>()); 00730 return *this; 00731 } 00732 00734 inline DynMatrixColumn &operator*=(const T&val){ 00735 std::for_each(begin(),end(),std::bind2nd(std::multiplies<T>(),val)); 00736 return *this; 00737 } 00739 inline DynMatrixColumn &operator/=(const T&val){ 00740 std::for_each(begin(),end(),std::bind2nd(std::divides<T>(),val)); 00741 return *this; 00742 } 00743 00744 }; 00745 00746 inline DynMatrix &operator=(const DynMatrixColumn &col){ 00747 DYN_MATRIX_COLUMN_CHECK(dim() != col.dim(),"dimension missmatch"); 00748 std::copy(col.begin(),col.end(),begin()); 00749 return *this; 00750 } 00751 00752 #undef DYN_MATRIX_COLUMN_CHECK 00753 00754 00755 00757 inline iterator begin() { return m_data; } 00758 00760 inline iterator end() { return m_data+dim(); } 00761 00763 inline const_iterator begin() const { return m_data; } 00764 00766 inline const_iterator end() const { return m_data+dim(); } 00767 00769 inline col_iterator col_begin(unsigned int col) { 00770 col_check(col); 00771 return col_iterator(m_data+col,cols()); 00772 } 00773 00775 inline col_iterator col_end(unsigned int col) { 00776 col_check(col); 00777 return col_iterator(m_data+col+dim(),cols()); 00778 } 00779 00781 inline const_col_iterator col_begin(unsigned int col) const { 00782 col_check(col); 00783 return col_iterator(m_data+col,cols()); 00784 } 00785 00787 inline const_col_iterator col_end(unsigned int col) const { 00788 col_check(col); 00789 return col_iterator(m_data+col+dim(),cols()); 00790 } 00791 00793 inline row_iterator row_begin(unsigned int row) { 00794 row_check(row); 00795 return m_data+row*cols(); 00796 } 00797 00799 inline row_iterator row_end(unsigned int row) { 00800 row_check(row); 00801 return m_data+(row+1)*cols(); 00802 } 00803 00805 inline const_row_iterator row_begin(unsigned int row) const { 00806 row_check(row); 00807 return m_data+row*cols(); 00808 } 00809 00811 inline const_row_iterator row_end(unsigned int row) const { 00812 row_check(row); 00813 return m_data+(row+1)*cols(); 00814 } 00815 00817 inline DynMatrix row(int row){ 00818 row_check(row); 00819 return DynMatrix(m_cols,1,row_begin(row),false); 00820 } 00821 00823 inline const DynMatrix row(int row) const{ 00824 row_check(row); 00825 return DynMatrix(m_cols,1,const_cast<T*>(row_begin(row)),false); 00826 } 00827 00829 inline DynMatrixColumn col(int col){ 00830 return DynMatrixColumn(this,col); 00831 } 00832 00833 inline const DynMatrixColumn col(int col) const{ 00834 return DynMatrixColumn(this,col); 00835 } 00836 00838 void decompose_QR(DynMatrix &Q, DynMatrix &R) const; 00839 00841 void decompose_RQ(DynMatrix &R, DynMatrix &Q) const; 00842 00844 00846 void decompose_LU(DynMatrix &L, DynMatrix &U, T zeroThreshold = 1E-16) const; 00847 00849 DynMatrix solve_upper_triangular(const DynMatrix &b) const throw(InvalidMatrixDimensionException); 00850 00852 DynMatrix solve_lower_triangular(const DynMatrix &b) const throw(InvalidMatrixDimensionException); 00853 00855 00902 DynMatrix solve(const DynMatrix &b, const std::string &method = "lu", T zeroThreshold = 1E-16) 00903 throw(InvalidMatrixDimensionException, utils::ICLException, SingularMatrixException); 00904 00905 00907 DynMatrix inv() const throw (InvalidMatrixDimensionException, SingularMatrixException); 00908 00910 00923 void eigen(DynMatrix &eigenvectors, DynMatrix &eigenvalues) const throw(InvalidMatrixDimensionException, utils::ICLException); 00924 00926 00933 void svd(DynMatrix &U, DynMatrix &S, DynMatrix &V) const throw (utils::ICLException); 00934 00936 00960 DynMatrix pinv(bool useSVD = false, T zeroThreshold = 1E-16) const 00961 throw (InvalidMatrixDimensionException,SingularMatrixException,utils::ICLException); 00962 00964 00970 DynMatrix big_matrix_pinv(T zeroThreshold = 1E-16) const 00971 throw (InvalidMatrixDimensionException,SingularMatrixException,utils::ICLException); 00972 00973 #ifdef ICL_HAVE_MKL 00974 typedef void(*GESDD)(const char*, const int*, const int*, T*, const int*, T*, T*, const int*, T*, const int*, T*, const int*, int*, int*); 00975 typedef void(*CBLAS_GEMM)(CBLAS_ORDER,CBLAS_TRANSPOSE,CBLAS_TRANSPOSE,int,int,int,T,const T*,int,const T*,int,T,T*,int); 00976 DynMatrix big_matrix_pinv(T zeroThreshold, GESDD gesdd, CBLAS_GEMM cblas_gemm) const 00977 throw (InvalidMatrixDimensionException,SingularMatrixException,utils::ICLException); 00978 #endif 00979 00981 T det() const throw (InvalidMatrixDimensionException); 00982 00984 inline DynMatrix transp() const{ 00985 DynMatrix d(rows(),cols()); 00986 for(unsigned int x=0;x<cols();++x){ 00987 for(unsigned int y=0;y<rows();++y){ 00988 d(y,x) = (*this)(x,y); 00989 } 00990 } 00991 return d; 00992 } 00993 00995 00996 inline const DynMatrix<T> shallowTransposed() const{ 00997 return DynMatrix<T>(m_rows,m_cols,const_cast<T*>(m_data),false); 00998 } 00999 01001 const DynMatrix<T> shallowTransposed() { 01002 return DynMatrix<T>(m_rows,m_cols,const_cast<T*>(m_data),false); 01003 } 01004 01006 01015 inline void reshape(int newCols, int newRows) throw (InvalidMatrixDimensionException){ 01016 if((cols() * rows()) != (newCols * newRows)){ 01017 throw InvalidMatrixDimensionException("DynMatrix<T>::reshape: source dimension and destination dimension differs!"); 01018 } 01019 m_cols = newCols; 01020 m_rows = newRows; 01021 } 01022 01024 01025 T element_wise_inner_product(const DynMatrix<T> &other) const { 01026 return std::inner_product(begin(),end(),other.begin(),T(0)); 01027 } 01028 01029 01031 01034 DynMatrix<T> dot(const DynMatrix<T> &M) const throw(InvalidMatrixDimensionException){ 01035 return this->transp() * M; 01036 } 01037 01038 01040 DynMatrix<T> diag() const{ 01041 ICLASSERT_RETURN_VAL(cols()==rows(),DynMatrix<T>()); 01042 DynMatrix<T> d(1,rows()); 01043 for(int i=0;i<rows();++i){ 01044 d[i] = (*this)(i,i); 01045 } 01046 return d; 01047 } 01048 01050 T trace() const{ 01051 ICLASSERT_RETURN_VAL(cols()==rows(),0); 01052 double accu = 0; 01053 for(unsigned int i=0;i<dim();i+=cols()+1){ 01054 accu += m_data[i]; 01055 } 01056 return accu; 01057 } 01058 01060 static DynMatrix<T> cross(const DynMatrix<T> &x, const DynMatrix<T> &y) throw(InvalidMatrixDimensionException){ 01061 if(x.cols()==1 && y.cols()==1 && x.rows()==3 && y.rows()==3){ 01062 DynMatrix<T> r(1,x.rows()); 01063 r(0,0) = x(0,1)*y(0,2)-x(0,2)*y(0,1); 01064 r(0,1) = x(0,2)*y(0,0)-x(0,0)*y(0,2); 01065 r(0,2) = x(0,0)*y(0,1)-x(0,1)*y(0,0); 01066 return r; 01067 }else{ 01068 ICLASSERT_RETURN_VAL(x.rows() == 3 && y.rows() == 3,DynMatrix<T>()); 01069 return DynMatrix<T>(); 01070 } 01071 } 01072 01074 T cond(const double p=2) const { 01075 if(cols() == 3 && rows() == 3){ 01076 DynMatrix<T> M_inv = (*this).inv(); 01077 return (*this).norm(p) * M_inv.norm(p); 01078 } else { 01079 DynMatrix<T> U,S,V; 01080 (*this).svd(U,S,V); 01081 if(S[S.rows()-1]){ 01082 return S[0]/S[S.rows()-1]; 01083 } else { 01084 return S[0]; 01085 } 01086 } 01087 } 01088 01090 inline T *set_data(T *newData){ 01091 T *old_data = m_data; 01092 m_data = newData; 01093 return old_data; 01094 } 01095 01097 static inline DynMatrix id(unsigned int dim) { 01098 DynMatrix M(dim,dim); 01099 for(unsigned int i=0;i<dim;++i) M(i,i) = 1; 01100 return M; 01101 } 01102 01103 private: 01104 inline void row_check(unsigned int row) const{ 01105 #ifdef DYN_MATRIX_INDEX_CHECK 01106 if((int)row >= m_rows) ERROR_LOG("access to row index " << row << " on a "<<m_cols<<'x'<<m_rows<<"-matrix"); 01107 #else 01108 (void)row; 01109 #endif 01110 } 01111 inline void col_check(unsigned int col) const{ 01112 #ifdef DYN_MATRIX_INDEX_CHECK 01113 if((int)col >= m_cols) ERROR_LOG("access to column index " << col << " on a "<<m_cols<<'x'<<m_rows<<"-matrix"); 01114 #else 01115 (void)col; 01116 #endif 01117 } 01118 inline void idx_check(unsigned int col, unsigned int row) const{ 01119 col_check(col); 01120 row_check(row); 01121 } 01122 01123 inline void idx_check(unsigned int idx) const{ 01124 #ifdef DYN_MATRIX_INDEX_CHECK 01125 if(idx >= dim()) ERROR_LOG("access to linear index " << idx << " on a "<<m_cols<<'x'<<m_rows<<"-matrix"); 01126 #else 01127 (void)idx; 01128 #endif 01129 } 01130 01131 int m_rows; 01132 int m_cols; 01133 T *m_data; 01134 bool m_ownData; 01135 }; 01136 01138 01139 template<class T> 01140 DynMatrix<T>::DynMatrix(const typename DynMatrix<T>::DynMatrixColumn &column) : 01141 m_rows(column.dim()),m_cols(1),m_data(new T[column.dim()]),m_ownData(true){ 01142 std::copy(column.begin(),column.end(),begin()); 01143 } 01146 01147 template<class T> ICLMath_IMP 01148 std::ostream &operator<<(std::ostream &s, const DynMatrix<T> &m); 01149 01151 template<class T> ICLMath_IMP 01152 std::istream &operator>>(std::istream &s, DynMatrix<T> &m); 01153 01154 01155 #ifdef ICL_HAVE_IPP 01156 01157 template<> 01158 inline float DynMatrix<float>::sqrDistanceTo(const DynMatrix<float> &other) const throw (InvalidMatrixDimensionException){ 01159 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::sqrDistanceTo: dimension missmatch")); 01160 float norm = 0 ; 01161 ippsNormDiff_L2_32f(begin(), other.begin(), dim(), &norm); 01162 return norm*norm; 01163 } 01164 01165 template<> 01166 inline double DynMatrix<double>::sqrDistanceTo(const DynMatrix<double> &other) const throw (InvalidMatrixDimensionException){ 01167 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::sqrDistanceTo: dimension missmatch")); 01168 double norm = 0 ; 01169 ippsNormDiff_L2_64f(begin(), other.begin(), dim(), &norm); 01170 return norm*norm; 01171 } 01172 01173 template<> 01174 inline float DynMatrix<float>::distanceTo(const DynMatrix<float> &other) const throw (InvalidMatrixDimensionException){ 01175 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::distanceTo: dimension missmatch")); 01176 float norm = 0 ; 01177 ippsNormDiff_L2_32f(begin(), other.begin(), dim(), &norm); 01178 return norm; 01179 } 01180 01181 template<> 01182 inline double DynMatrix<double>::distanceTo(const DynMatrix<double> &other) const throw (InvalidMatrixDimensionException){ 01183 ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::distanceTo: dimension missmatch")); 01184 double norm = 0 ; 01185 ippsNormDiff_L2_64f(begin(), other.begin(), dim(), &norm); 01186 return norm; 01187 } 01188 01189 #define DYN_MATRIX_MULT_SPECIALIZE(IPPT) \ 01190 template<> \ 01191 inline DynMatrix<Ipp##IPPT> &DynMatrix<Ipp##IPPT>::mult( \ 01192 const DynMatrix<Ipp##IPPT> &m, DynMatrix<Ipp##IPPT> &dst) const \ 01193 throw (IncompatibleMatrixDimensionException){ \ 01194 if(cols() != m.rows() ) throw IncompatibleMatrixDimensionException("A*B : cols(A) must be row(B)"); \ 01195 dst.setBounds(m.cols(),rows()); \ 01196 ippmMul_mm_##IPPT(data(),sizeof(Ipp##IPPT)*cols(),sizeof(Ipp##IPPT),cols(),rows(), \ 01197 m.data(),sizeof(Ipp##IPPT)*m.cols(),sizeof(Ipp##IPPT),m.cols(),m.rows(), \ 01198 dst.data(),m.cols()*sizeof(Ipp##IPPT),sizeof(Ipp##IPPT)); \ 01199 return dst; \ 01200 } 01201 01202 DYN_MATRIX_MULT_SPECIALIZE(32f) 01203 DYN_MATRIX_MULT_SPECIALIZE(64f) 01204 #undef DYN_MATRIX_MULT_SPECIALIZE 01205 01206 01207 #define DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE(IPPT) \ 01208 template<> \ 01209 inline DynMatrix<Ipp##IPPT> &DynMatrix<Ipp##IPPT>::elementwise_div( \ 01210 const DynMatrix<Ipp##IPPT> &m, DynMatrix<Ipp##IPPT> &dst) const \ 01211 throw (IncompatibleMatrixDimensionException){ \ 01212 if((m.cols() != cols()) || (m.rows() != rows())){ \ 01213 throw IncompatibleMatrixDimensionException("A./B dimension mismatch"); \ 01214 } \ 01215 dst.setBounds(cols(),rows()); \ 01216 ippsDiv_##IPPT(data(),m.data(),dst.data(),dim()); \ 01217 return dst; \ 01218 } 01219 01220 DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE(32f) 01221 DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE(64f) 01222 01223 #undef DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE 01224 01225 01226 01227 01228 01229 01230 #define DYN_MATRIX_MULT_BY_CONSTANT(IPPT) \ 01231 template<> \ 01232 inline DynMatrix<Ipp##IPPT> &DynMatrix<Ipp##IPPT>::mult( \ 01233 Ipp##IPPT f, DynMatrix<Ipp##IPPT> &dst) const{ \ 01234 dst.setBounds(cols(),rows()); \ 01235 ippsMulC_##IPPT(data(),f, dst.data(),dim()); \ 01236 return dst; \ 01237 } 01238 01239 DYN_MATRIX_MULT_BY_CONSTANT(32f) 01240 DYN_MATRIX_MULT_BY_CONSTANT(64f) 01241 01242 #undef DYN_MATRIX_MULT_BY_CONSTANT 01243 01244 #define DYN_MATRIX_NORM_SPECIALZE(T,IPPT) \ 01245 template<> \ 01246 inline T DynMatrix<T> ::norm(double l) const{ \ 01247 if(l==1){ \ 01248 T val; \ 01249 ippsNorm_L1_##IPPT(m_data,dim(),&val); \ 01250 return val; \ 01251 }else if(l==2){ \ 01252 T val; \ 01253 ippsNorm_L2_##IPPT(m_data,dim(),&val); \ 01254 return val; \ 01255 } \ 01256 double accu = 0; \ 01257 for(unsigned int i=0;i<dim();++i){ \ 01258 accu += ::pow(double(m_data[i]),l); \ 01259 } \ 01260 return ::pow(accu,1.0/l); \ 01261 } 01262 01263 DYN_MATRIX_NORM_SPECIALZE(float,32f) 01264 // DYN_MATRIX_NORM_SPECIALZE(double,64f) 01265 01266 #undef DYN_MATRIX_NORM_SPECIALZE 01267 01270 #endif // ICL_HAVE_IPP 01271 01273 01274 template<class T> 01275 inline DynMatrix<T> operator,(const DynMatrix<T> &left, const DynMatrix<T> &right){ 01276 int w = left.cols() + right.cols(); 01277 int h = iclMax(left.rows(),right.rows()); 01278 DynMatrix<T> result(w,h,float(0)); 01279 for(unsigned int y=0;y<left.rows();++y){ 01280 std::copy(left.row_begin(y), left.row_end(y), result.row_begin(y)); 01281 } 01282 for(unsigned int y=0;y<right.rows();++y){ 01283 std::copy(right.row_begin(y), right.row_end(y), result.row_begin(y) + left.cols()); 01284 } 01285 return result; 01286 } 01287 01289 01290 template<class T> 01291 inline DynMatrix<T> operator%(const DynMatrix<T> &top, const DynMatrix<T> &bottom){ 01292 int w = iclMax(top.cols(),bottom.cols()); 01293 int h = top.rows() + bottom.rows(); 01294 DynMatrix<T> result(w,h,float(0)); 01295 for(unsigned int y=0;y<top.rows();++y){ 01296 std::copy(top.row_begin(y), top.row_end(y), result.row_begin(y)); 01297 } 01298 for(unsigned int y=0;y<bottom.rows();++y){ 01299 std::copy(bottom.row_begin(y), bottom.row_end(y), result.row_begin(y+top.rows())); 01300 } 01301 return result; 01302 } 01303 } // namespace math 01304 }