// Copyright (C) 2006 Davis E. King (davisking@users.sourceforge.net) // License: Boost Software License See LICENSE.txt for the full license. #ifndef DLIB_MATRIx_ #define DLIB_MATRIx_ #include "matrix_abstract.h" #include "../algs.h" #include "../serialize.h" #include "../enable_if.h" #include <sstream> #include <algorithm> #include "../memory_manager.h" #include "../is_kind.h" #ifdef _MSC_VER // Disable the following warnings for Visual Studio // This warning is: // "warning C4355: 'this' : used in base member initializer list" // Which we get from this code but it is not an error so I'm turning this // warning off and then turning it back on at the end of the file. #pragma warning(disable : 4355) #endif namespace dlib { // ---------------------------------------------------------------------------------------- template < typename T, long num_rows = 0, long num_cols = 0, typename mem_manager = memory_manager<char>::kernel_1a > class matrix; // ---------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------- template < typename T, long num_rows, long num_cols, typename mem_manager > class matrix_ref { public: typedef T type; typedef matrix_ref ref_type; typedef mem_manager mem_manager_type; const static long NR = num_rows; const static long NC = num_cols; matrix_ref ( const matrix<T,num_rows,num_cols,mem_manager>& m_ ) : m(m_) {} matrix_ref ( const matrix_ref& i_ ) : m(i_.m) {} const T& operator() ( long r, long c ) const { return m(r,c); } long nr ( ) const { return m.nr(); } long nc ( ) const { return m.nc(); } long size ( ) const { return m.size(); } template <typename U, long iNR, long iNC, typename mm > bool aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return false; } template <typename U, long iNR, long iNC, typename mm> bool destructively_aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return false; } bool aliases ( const matrix<T,num_rows,num_cols,mem_manager>& item ) const { return (&m == &item); } const matrix_ref ref( ) const { return *this; } private: // no assignment operator matrix_ref& operator=(const matrix_ref&); const matrix<T,num_rows,num_cols,mem_manager>& m; // This is the item contained by this expression. }; // ---------------------------------------------------------------------------------------- // this is a hack to avoid a compile time error in visual studio 8. I would just // use sizeof(T) and be done with it but that won't compile. The idea here // is to avoid using the stack allocation of the matrix_data object if it // is going to contain another matrix and also avoid asking for the sizeof() // the contained matrix. template <typename T> struct get_sizeof_helper { const static std::size_t val = sizeof(T); }; template <typename T, long NR, long NC, typename mm> struct get_sizeof_helper<matrix<T,NR,NC,mm> > { const static std::size_t val = 1000000; }; template < typename T, long num_rows, long num_cols, typename mem_manager, int val = static_switch < // when the sizes are all non zero and small (num_rows*num_cols*get_sizeof_helper<T>::val <= 64) && (num_rows != 0 && num_cols != 0), // when the sizes are all non zero and big (num_rows*num_cols*get_sizeof_helper<T>::val >= 65) && (num_rows != 0 && num_cols != 0), num_rows == 0 && num_cols != 0, num_rows != 0 && num_cols == 0, num_rows == 0 && num_cols == 0 >::value > class matrix_data ; /*! WHAT THIS OBJECT REPRESENTS This object represents the actual allocation of space for a matrix. Small matrices allocate all their data on the stack and bigger ones use a memory_manager to get their memory. !*/ // ---------------------------------------------------------------------------------------- template < typename T, long num_rows, long num_cols, typename mem_manager > class matrix_data<T,num_rows,num_cols,mem_manager,1> : noncopyable // when the sizes are all non zero and small { public: const static long NR = num_rows; const static long NC = num_cols; matrix_data() {} T& operator() ( long r, long c ) { return data[r][c]; } const T& operator() ( long r, long c ) const { return data[r][c]; } T& operator() ( long i ) { return *(*data + i); } const T& operator() ( long i ) const { return *(*data + i); } void swap( matrix_data& item ) { for (long r = 0; r < num_rows; ++r) { for (long c = 0; c < num_cols; ++c) { exchange((*this)(r,c),item(r,c)); } } } long nr ( ) const { return num_rows; } long nc ( ) const { return num_cols; } void set_size ( long nr, long nc ) { } private: T data[num_rows][num_cols]; }; // ---------------------------------------------------------------------------------------- template < typename T, long num_rows, long num_cols, typename mem_manager > class matrix_data<T,num_rows,num_cols,mem_manager,2> : noncopyable // when the sizes are all non zero and big { public: const static long NR = num_rows; const static long NC = num_cols; matrix_data ( ) { data = pool.allocate_array(num_rows*num_cols); } ~matrix_data () { pool.deallocate_array(data); } T& operator() ( long r, long c ) { return data[r*num_cols + c]; } const T& operator() ( long r, long c ) const { return data[r*num_cols + c]; } T& operator() ( long i ) { return data[i]; } const T& operator() ( long i ) const { return data[i]; } void swap( matrix_data& item ) { std::swap(item.data,data); pool.swap(item.pool); } long nr ( ) const { return num_rows; } long nc ( ) const { return num_cols; } void set_size ( long nr, long nc ) { } private: T* data; typename mem_manager::template rebind<T>::other pool; }; // ---------------------------------------------------------------------------------------- template < typename T, long num_rows, long num_cols, typename mem_manager > class matrix_data<T,num_rows,num_cols,mem_manager,3> : noncopyable // when num_rows == 0 && num_cols != 0, { public: const static long NR = num_rows; const static long NC = num_cols; matrix_data ( ):data(0), nr_(0) { } ~matrix_data () { if (data) pool.deallocate_array(data); } T& operator() ( long r, long c ) { return data[r*num_cols + c]; } const T& operator() ( long r, long c ) const { return data[r*num_cols + c]; } T& operator() ( long i ) { return data[i]; } const T& operator() ( long i ) const { return data[i]; } void swap( matrix_data& item ) { std::swap(item.data,data); std::swap(item.nr_,nr_); pool.swap(item.pool); } long nr ( ) const { return nr_; } long nc ( ) const { return num_cols; } void set_size ( long nr, long nc ) { if (data) { pool.deallocate_array(data); } data = pool.allocate_array(nr*nc); nr_ = nr; } private: T* data; long nr_; typename mem_manager::template rebind<T>::other pool; }; // ---------------------------------------------------------------------------------------- template < typename T, long num_rows, long num_cols, typename mem_manager > class matrix_data<T,num_rows,num_cols,mem_manager,4> : noncopyable // when num_rows != 0 && num_cols == 0 { public: const static long NR = num_rows; const static long NC = num_cols; matrix_data ( ):data(0), nc_(0) { } ~matrix_data () { if (data) { pool.deallocate_array(data); } } T& operator() ( long r, long c ) { return data[r*nc_ + c]; } const T& operator() ( long r, long c ) const { return data[r*nc_ + c]; } T& operator() ( long i ) { return data[i]; } const T& operator() ( long i ) const { return data[i]; } void swap( matrix_data& item ) { std::swap(item.data,data); std::swap(item.nc_,nc_); pool.swap(item.pool); } long nr ( ) const { return num_rows; } long nc ( ) const { return nc_; } void set_size ( long nr, long nc ) { if (data) { pool.deallocate_array(data); } data = pool.allocate_array(nr*nc); nc_ = nc; } private: T* data; long nc_; typename mem_manager::template rebind<T>::other pool; }; // ---------------------------------------------------------------------------------------- template < typename T, long num_rows, long num_cols, typename mem_manager > class matrix_data<T,num_rows,num_cols,mem_manager,5> : noncopyable // when num_rows == 0 && num_cols == 0 { public: const static long NR = num_rows; const static long NC = num_cols; matrix_data ( ):data(0), nr_(0), nc_(0) { } ~matrix_data () { if (data) { pool.deallocate_array(data); } } T& operator() ( long r, long c ) { return data[r*nc_ + c]; } const T& operator() ( long r, long c ) const { return data[r*nc_ + c]; } T& operator() ( long i ) { return data[i]; } const T& operator() ( long i ) const { return data[i]; } void swap( matrix_data& item ) { std::swap(item.data,data); std::swap(item.nc_,nc_); std::swap(item.nr_,nr_); pool.swap(item.pool); } long nr ( ) const { return nr_; } long nc ( ) const { return nc_; } void set_size ( long nr, long nc ) { if (data) { pool.deallocate_array(data); } data = pool.allocate_array(nr*nc); nr_ = nr; nc_ = nc; } private: T* data; long nr_; long nc_; typename mem_manager::template rebind<T>::other pool; }; // ---------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------- // We want to return the compile time constant if our NR and NC dimensions // aren't zero but if they are then we want to call ref_.nx() and return // the correct values. template < typename ref_type, long NR > struct get_nr_helper { static inline long get(const ref_type&) { return NR; } }; template < typename ref_type > struct get_nr_helper<ref_type,0> { static inline long get(const ref_type& m) { return m.nr(); } }; template < typename ref_type, long NC > struct get_nc_helper { static inline long get(const ref_type&) { return NC; } }; template < typename ref_type > struct get_nc_helper<ref_type,0> { static inline long get(const ref_type& m) { return m.nc(); } }; // the matrix_exp for statically sized matrices template < typename EXP > class matrix_exp { public: typedef typename EXP::type type; typedef typename EXP::ref_type ref_type; typedef typename EXP::mem_manager_type mem_manager_type; const static long NR = EXP::NR; const static long NC = EXP::NC; typedef matrix<type,NR,NC,mem_manager_type> matrix_type; matrix_exp ( const EXP& exp ) : ref_(exp.ref()) {} inline const type operator() ( long r, long c ) const { DLIB_ASSERT(r < nr() && c < nc() && r >= 0 && c >= c, "\tconst type matrix_exp::operator(r,c)" << "\n\tYou must give a valid row and column" << "\n\tr: " << r << "\n\tc: " << c << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tthis: " << this ); return ref_(r,c); } const type operator() ( long i ) const { COMPILE_TIME_ASSERT(NC == 1 || NC == 0 || NR == 1 || NR == 0); DLIB_ASSERT(nc() == 1 || nr() == 1, "\tconst type matrix_exp::operator(i)" << "\n\tYou can only use this operator on column or row vectors" << "\n\ti: " << i << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tthis: " << this ); DLIB_ASSERT( ((nc() == 1 && i < nr()) || (nr() == 1 && i < nc())) && i >= 0, "\tconst type matrix_exp::operator(i)" << "\n\tYou must give a valid row/column number" << "\n\ti: " << i << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tthis: " << this ); if (nc() == 1) return ref_(i,0); else return ref_(0,i); } long size ( ) const { return nr()*nc(); } long nr ( ) const { return get_nr_helper<ref_type,NR>::get(ref_); } long nc ( ) const { return get_nc_helper<ref_type,NC>::get(ref_); } template <typename U, long iNR, long iNC, typename mm > bool aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return ref_.aliases(item); } template <typename U, long iNR, long iNC , typename mm> bool destructively_aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return ref_.destructively_aliases(item); } const ref_type& ref ( ) const { return ref_; } inline operator const type ( ) const { COMPILE_TIME_ASSERT(NC == 1 || NC == 0); COMPILE_TIME_ASSERT(NR == 1 || NR == 0); DLIB_ASSERT(nr() == 1 && nc() == 1, "\tmatrix_exp::operator const type&() const" << "\n\tYou can only use this operator on a 1x1 matrix" << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tthis: " << this ); return ref_(0,0); } private: const ref_type ref_; }; // ---------------------------------------------------------------------------------------- template <typename T> struct is_matrix<matrix_exp<T> > { static const bool value = true; }; template <typename T, long NR, long NC, typename mm> struct is_matrix<matrix_ref<T,NR,NC,mm> > { static const bool value = true; }; template <typename T, long NR, long NC, typename mm> struct is_matrix<matrix<T,NR,NC,mm> > { static const bool value = true; }; template <typename T> struct is_matrix<T&> { static const bool value = is_matrix<T>::value; }; template <typename T> struct is_matrix<const T&> { static const bool value = is_matrix<T>::value; }; template <typename T> struct is_matrix<const T> { static const bool value = is_matrix<T>::value; }; /* is_matrix<T>::value == 1 if T is a matrix type else 0 */ // ---------------------------------------------------------------------------------------- // This template will perform the needed loop for element multiplication using whichever // dimension is provided as a compile time constant (if one is at all). template < typename LHS, typename RHS, long lhs_nc = LHS::NC, long rhs_nr = RHS::NR > struct matrix_multiply_helper { typedef typename LHS::type type; inline const static type eval ( const RHS& rhs, const LHS& lhs, long r, long c ) { type temp = type(); for (long i = 0; i < rhs.nr(); ++i) { temp += lhs(r,i)*rhs(i,c); } return temp; } }; template < typename LHS, typename RHS, long lhs_nc > struct matrix_multiply_helper <LHS,RHS,lhs_nc,0> { typedef typename LHS::type type; inline const static type eval ( const RHS& rhs, const LHS& lhs, long r, long c ) { type temp = type(); for (long i = 0; i < lhs.nc(); ++i) { temp += lhs(r,i)*rhs(i,c); } return temp; } }; template < typename LHS, typename RHS, unsigned long count = 0 > class matrix_multiply_exp { /*! REQUIREMENTS ON LHS AND RHS - they must be matrix_exp or matrix_ref objects (or objects with a compatible interface). !*/ public: typedef typename LHS::type type; typedef matrix_multiply_exp ref_type; typedef typename LHS::mem_manager_type mem_manager_type; const static long NR = LHS::NR; const static long NC = RHS::NC; matrix_multiply_exp ( const matrix_multiply_exp& item ) : lhs(item.lhs), rhs(item.rhs) {} inline matrix_multiply_exp ( const LHS& lhs_, const RHS& rhs_ ) : lhs(lhs_), rhs(rhs_) { // You are trying to multiply two incompatible matrices together. The number of columns // in the matrix on the left must match the number of rows in the matrix on the right. COMPILE_TIME_ASSERT(LHS::NC == RHS::NR || LHS::NC*RHS::NR == 0); DLIB_ASSERT(lhs.nc() == rhs.nr(), "\tconst matrix_exp operator*(const matrix_exp& lhs, const matrix_exp& rhs)" << "\n\tYou are trying to multiply two incompatible matrices together" << "\n\tlhs.nr(): " << lhs.nr() << "\n\tlhs.nc(): " << lhs.nc() << "\n\trhs.nr(): " << rhs.nr() << "\n\trhs.nc(): " << rhs.nc() << "\n\t&lhs: " << &lhs << "\n\t&rhs: " << &rhs ); // You can't multiply matrices together if they don't both contain the same type of elements. COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true)); } inline const type operator() ( long r, long c ) const { return matrix_multiply_helper<LHS,RHS>::eval(rhs,lhs,r,c); } long nr ( ) const { return lhs.nr(); } long nc ( ) const { return rhs.nc(); } template <typename U, long iNR, long iNC, typename mm > bool aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return lhs.aliases(item) || rhs.aliases(item); } template <typename U, long iNR, long iNC , typename mm> bool destructively_aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return aliases(item); } const ref_type& ref( ) const { return *this; } const LHS lhs; const RHS rhs; }; template < typename T, long NR, long NC, typename EXP1, typename EXP2, typename MM > inline const matrix_exp<matrix_multiply_exp<EXP1, matrix_multiply_exp<EXP2,typename matrix<T,NR,NC,MM>::ref_type >,0 > > operator* ( const matrix_exp<matrix_multiply_exp<EXP1,EXP2,1> >& m1, const matrix<T,NR,NC,MM>& m2 ) { // We are going to reorder the order of evaluation of the terms here. This way the // multiplication will go faster. typedef matrix_multiply_exp<EXP2,typename matrix<T,NR,NC,MM>::ref_type > exp_inner; typedef matrix_multiply_exp<EXP1, exp_inner,0 > exp_outer; return matrix_exp<exp_outer>(exp_outer(m1.ref().lhs,exp_inner(m1.ref().rhs,m2))); } template < typename EXP1, typename EXP2 > inline const matrix_exp<matrix_multiply_exp<EXP1, EXP2 > > operator* ( const matrix_exp<EXP1>& m1, const matrix_exp<EXP2>& m2 ) { typedef matrix_multiply_exp<EXP1, EXP2> exp; return matrix_exp<exp>(exp(m1.ref(),m2.ref())); } template < typename T, long NR, long NC, typename EXP, typename MM > inline const matrix_exp<matrix_multiply_exp<typename matrix<T,NR,NC,MM>::ref_type, matrix_exp<EXP> > > operator* ( const matrix<T,NR,NC,MM>& m1, const matrix_exp<EXP>& m2 ) { typedef matrix_multiply_exp<typename matrix<T,NR,NC,MM>::ref_type, matrix_exp<EXP> > exp; return matrix_exp<exp>(exp(m1,m2)); } template < typename T, long NR, long NC, typename EXP, typename MM > inline const matrix_exp<matrix_multiply_exp< matrix_exp<EXP>, typename matrix<T,NR,NC,MM>::ref_type, 1> > operator* ( const matrix_exp<EXP>& m1, const matrix<T,NR,NC,MM>& m2 ) { typedef matrix_multiply_exp< matrix_exp<EXP>, typename matrix<T,NR,NC,MM>::ref_type, 1 > exp; return matrix_exp<exp>(exp(m1,m2)); } template < typename T, long NR1, long NC1, long NR2, long NC2, typename MM1, typename MM2 > inline const matrix_exp<matrix_multiply_exp<typename matrix<T,NR1,NC1,MM1>::ref_type,typename matrix<T,NR2,NC2,MM2>::ref_type > > operator* ( const matrix<T,NR1,NC1,MM1>& m1, const matrix<T,NR2,NC2,MM2>& m2 ) { typedef matrix_multiply_exp<typename matrix<T,NR1,NC1,MM1>::ref_type, typename matrix<T,NR2,NC2,MM2>::ref_type > exp; return matrix_exp<exp>(exp(m1,m2)); } // ---------------------------------------------------------------------------------------- template < typename LHS, typename RHS > class matrix_add_expression { /*! REQUIREMENTS ON LHS AND RHS - they must be matrix_exp or matrix_ref objects (or objects with a compatible interface). !*/ public: typedef typename LHS::type type; typedef typename LHS::mem_manager_type mem_manager_type; typedef matrix_add_expression ref_type; const static long NR = (RHS::NR > LHS::NR) ? RHS::NR : LHS::NR; const static long NC = (RHS::NC > LHS::NC) ? RHS::NC : LHS::NC; matrix_add_expression ( const matrix_add_expression& item ) : lhs(item.lhs), rhs(item.rhs) {} matrix_add_expression ( const LHS& lhs_, const RHS& rhs_ ) : lhs(lhs_), rhs(rhs_) { // You can only add matrices together if they both have the same number of rows and columns. COMPILE_TIME_ASSERT(LHS::NR == RHS::NR || LHS::NR == 0 || RHS::NR == 0); COMPILE_TIME_ASSERT(LHS::NC == RHS::NC || LHS::NC == 0 || RHS::NC == 0); DLIB_ASSERT(lhs.nc() == rhs.nc() && lhs.nr() == rhs.nr(), "\tconst matrix_exp operator+(const matrix_exp& lhs, const matrix_exp& rhs)" << "\n\tYou are trying to add two incompatible matrices together" << "\n\tlhs.nr(): " << lhs.nr() << "\n\tlhs.nc(): " << lhs.nc() << "\n\trhs.nr(): " << rhs.nr() << "\n\trhs.nc(): " << rhs.nc() << "\n\t&lhs: " << &lhs << "\n\t&rhs: " << &rhs ); // You can only add matrices together if they both contain the same types of elements. COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true)); } const type operator() ( long r, long c ) const { return lhs(r,c) + rhs(r,c); } template <typename U, long iNR, long iNC , typename mm> bool aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return lhs.aliases(item) || rhs.aliases(item); } template <typename U, long iNR, long iNC, typename mm > bool destructively_aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return lhs.destructively_aliases(item) || rhs.destructively_aliases(item); } const ref_type& ref( ) const { return *this; } long nr ( ) const { return lhs.nr(); } long nc ( ) const { return lhs.nc(); } const LHS lhs; const RHS rhs; }; template < typename EXP1, typename EXP2 > inline const matrix_exp<matrix_add_expression<EXP1, EXP2 > > operator+ ( const matrix_exp<EXP1>& m1, const matrix_exp<EXP2>& m2 ) { typedef matrix_add_expression<EXP1, EXP2 > exp; return matrix_exp<exp>(exp(m1.ref(),m2.ref())); } // ---------------------------------------------------------------------------------------- template < typename LHS, typename RHS > class matrix_subtract_exp { /*! REQUIREMENTS ON LHS AND RHS - they must be matrix_exp or matrix_ref objects (or objects with a compatible interface). !*/ public: typedef typename LHS::type type; typedef typename LHS::mem_manager_type mem_manager_type; typedef matrix_subtract_exp ref_type; const static long NR = (RHS::NR > LHS::NR) ? RHS::NR : LHS::NR; const static long NC = (RHS::NC > LHS::NC) ? RHS::NC : LHS::NC; matrix_subtract_exp ( const LHS& lhs_, const RHS& rhs_ ) : lhs(lhs_), rhs(rhs_) { // You can only subtract one matrix from another if they both have the same number of rows and columns. COMPILE_TIME_ASSERT(LHS::NR == RHS::NR || LHS::NR == 0 || RHS::NR == 0); COMPILE_TIME_ASSERT(LHS::NC == RHS::NC || LHS::NC == 0 || RHS::NC == 0); DLIB_ASSERT(lhs.nc() == rhs.nc() && lhs.nr() == rhs.nr(), "\tconst matrix_exp operator-(const matrix_exp& lhs, const matrix_exp& rhs)" << "\n\tYou are trying to add two incompatible matrices together" << "\n\tlhs.nr(): " << lhs.nr() << "\n\tlhs.nc(): " << lhs.nc() << "\n\trhs.nr(): " << rhs.nr() << "\n\trhs.nc(): " << rhs.nc() << "\n\t&lhs: " << &lhs << "\n\t&rhs: " << &rhs ); // You can only subtract one matrix from another if they both contain elements of the same type. COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true)); } const type operator() ( long r, long c ) const { return lhs(r,c) - rhs(r,c); } template <typename U, long iNR, long iNC, typename mm > bool aliases ( const matrix<U,iNR,iNC, mm>& item ) const { return lhs.aliases(item) || rhs.aliases(item); } template <typename U, long iNR, long iNC , typename mm> bool destructively_aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return lhs.destructively_aliases(item) || rhs.destructively_aliases(item); } const ref_type& ref( ) const { return *this; } long nr ( ) const { return lhs.nr(); } long nc ( ) const { return lhs.nc(); } const LHS lhs; const RHS rhs; }; template < typename EXP1, typename EXP2 > inline const matrix_exp<matrix_subtract_exp<EXP1, EXP2 > > operator- ( const matrix_exp<EXP1>& m1, const matrix_exp<EXP2>& m2 ) { typedef matrix_subtract_exp<EXP1, EXP2 > exp; return matrix_exp<exp>(exp(m1.ref(),m2.ref())); } // ---------------------------------------------------------------------------------------- template < typename M, typename S > class matrix_divscal_exp { /*! REQUIREMENTS ON M - must be a matrix_exp or matrix_ref object (or an object with a compatible interface). REQUIREMENTS ON S - must be some kind of scalar type !*/ public: typedef typename M::type type; typedef typename M::mem_manager_type mem_manager_type; typedef matrix_divscal_exp ref_type; const static long NR = M::NR; const static long NC = M::NC; matrix_divscal_exp ( const M& m_, const S& s_ ) : m(m_), s(s_) {} const type operator() ( long r, long c ) const { return m(r,c)/s; } template <typename U, long iNR, long iNC, typename mm > bool aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return m.aliases(item); } template <typename U, long iNR, long iNC, typename mm > bool destructively_aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return m.destructively_aliases(item); } const ref_type& ref( ) const { return *this; } long nr ( ) const { return m.nr(); } long nc ( ) const { return m.nc(); } const M m; const S s; }; template < typename EXP, typename S > inline const matrix_exp<matrix_divscal_exp<matrix_exp<EXP>, S> > operator/ ( const matrix_exp<EXP>& m, const S& s ) { typedef matrix_divscal_exp<matrix_exp<EXP>,S > exp; return matrix_exp<exp>(exp(m,s)); } // ---------------------------------------------------------------------------------------- template < typename M, typename S > class matrix_mulscal_exp { /*! REQUIREMENTS ON M - must be a matrix_exp or matrix_ref object (or an object with a compatible interface). REQUIREMENTS ON S - must be some kind of scalar type !*/ public: typedef typename M::type type; typedef typename M::mem_manager_type mem_manager_type; typedef matrix_mulscal_exp ref_type; const static long NR = M::NR; const static long NC = M::NC; matrix_mulscal_exp ( const M& m_, const S& s_ ) : m(m_), s(s_) {} const type operator() ( long r, long c ) const { return m(r,c)*s; } template <typename U, long iNR, long iNC , typename mm> bool aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return m.aliases(item); } template <typename U, long iNR, long iNC, typename mm > bool destructively_aliases ( const matrix<U,iNR,iNC,mm>& item ) const { return m.destructively_aliases(item); } const ref_type& ref( ) const { return *this; } long nr ( ) const { return m.nr(); } long nc ( ) const { return m.nc(); } const M m; const S s; }; template < typename EXP, typename S > inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, S> > operator* ( const matrix_exp<EXP>& m, const S& s ) { typedef matrix_mulscal_exp<matrix_exp<EXP>,S > exp; return matrix_exp<exp>(exp(m,s)); } template < typename EXP, typename S > inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, S> > operator* ( const S& s, const matrix_exp<EXP>& m ) { typedef matrix_mulscal_exp<matrix_exp<EXP>,S > exp; return matrix_exp<exp>(exp(m,s)); } template < typename EXP > inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, float> > operator/ ( const matrix_exp<EXP>& m, const float& s ) { typedef matrix_mulscal_exp<matrix_exp<EXP>,float > exp; return matrix_exp<exp>(exp(m,1.0f/s)); } template < typename EXP > inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, double> > operator/ ( const matrix_exp<EXP>& m, const double& s ) { typedef matrix_mulscal_exp<matrix_exp<EXP>,double > exp; return matrix_exp<exp>(exp(m,1.0/s)); } template < typename EXP > inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, long double> > operator/ ( const matrix_exp<EXP>& m, const long double& s ) { typedef matrix_mulscal_exp<matrix_exp<EXP>,long double > exp; return matrix_exp<exp>(exp(m,1.0/s)); } template < typename EXP > inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, int> > operator- ( const matrix_exp<EXP>& m ) { typedef matrix_mulscal_exp<matrix_exp<EXP>,int > exp; return matrix_exp<exp>(exp(m,-1)); } // ---------------------------------------------------------------------------------------- template < typename EXP1, typename EXP2 > bool operator== ( const matrix_exp<EXP1>& m1, const matrix_exp<EXP2>& m2 ) { if (m1.nr() == m2.nr() && m1.nc() == m2.nc()) { for (long r = 0; r < m1.nr(); ++r) { for (long c = 0; c < m1.nc(); ++c) { if (m1(r,c) != m2(r,c)) return false; } } return true; } return false; } template < typename EXP1, typename EXP2 > inline bool operator!= ( const matrix_exp<EXP1>& m1, const matrix_exp<EXP2>& m2 ) { return !(m1 == m2); } // ---------------------------------------------------------------------------------------- template < typename matrix_dest_type, typename src_exp > void matrix_assign ( matrix_dest_type& dest, const matrix_exp<src_exp>& src, const long row_offset = 0, const long col_offset = 0 ) /*! requires - src.destructively_aliases(dest) == false - dest.nr() == src.nr()-row_offset - dest.nc() == src.nc()-col_offset ensures - #subm(dest, row_offset, col_offset, src.nr(), src.nc()) == src - the part of dest outside the above sub matrix remains unchanged !*/ { for (long r = 0; r < src.nr(); ++r) { for (long c = 0; c < src.nc(); ++c) { dest(r+row_offset,c+col_offset) = src(r,c); } } } // ---------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------- template < typename T, long num_rows, long num_cols, typename mem_manager > class matrix : public matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager> > { COMPILE_TIME_ASSERT(num_rows >= 0 && num_cols >= 0); public: typedef T type; typedef matrix_ref<T,num_rows,num_cols,mem_manager> ref_type; typedef mem_manager mem_manager_type; const static long NR = num_rows; const static long NC = num_cols; matrix () : matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager> >(ref_type(*this)) { } explicit matrix ( long length ) : matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager> >(ref_type(*this)) { // This object you are trying to call matrix(length) on is not a column or // row vector. COMPILE_TIME_ASSERT(NR == 1 || NC == 1); DLIB_ASSERT( length >= 0, "\tmatrix::matrix(length)" << "\n\tlength must be at least 0" << "\n\tlength: " << length << "\n\tNR: " << NR << "\n\tNC: " << NC << "\n\tthis: " << this ); if (NR == 1) { DLIB_ASSERT(NC == 0 || NC == length, "\tmatrix::matrix(length)" << "\n\tSince this is a statically sized matrix length must equal NC" << "\n\tlength: " << length << "\n\tNR: " << NR << "\n\tNC: " << NC << "\n\tthis: " << this ); data.set_size(1,length); } else { DLIB_ASSERT(NR == 0 || NR == length, "\tvoid matrix::set_size(length)" << "\n\tSince this is a statically sized matrix length must equal NR" << "\n\tlength: " << length << "\n\tNR: " << NR << "\n\tNC: " << NC << "\n\tthis: " << this ); data.set_size(length,1); } } matrix ( long rows, long cols ) : matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager> >(ref_type(*this)) { DLIB_ASSERT( (NR == 0 || NR == rows) && ( NC == 0 || NC == cols) && rows >= 0 && cols >= 0, "\tvoid matrix::matrix(rows, cols)" << "\n\tYou have supplied conflicting matrix dimensions" << "\n\trows: " << rows << "\n\tcols: " << cols << "\n\tNR: " << NR << "\n\tNC: " << NC ); data.set_size(rows,cols); } template <typename EXP> matrix ( const matrix_exp<EXP>& m ): matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager> >(ref_type(*this)) { // You get an error on this line if the matrix m contains a type that isn't // the same as the type contained in the target matrix. COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true) || (is_matrix<typename EXP::type>::value == true)); // The matrix you are trying to assign m to is a statically sized matrix and // m's dimensions don't match that of *this. COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0); COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0); DLIB_ASSERT((NR == 0 || NR == m.nr()) && (NC == 0 || NC == m.nc()), "\tmatrix& matrix::matrix(const matrix_exp& m)" << "\n\tYou are trying to assign a dynamically sized matrix to a statically sized matrix with the wrong size" << "\n\tNR: " << NR << "\n\tNC: " << NC << "\n\tm.nr(): " << m.nr() << "\n\tm.nc(): " << m.nc() << "\n\tthis: " << this ); data.set_size(m.nr(),m.nc()); matrix_assign(*this, m); } matrix ( const matrix& m ): matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager> >(ref_type(*this)) { data.set_size(m.nr(),m.nc()); matrix_assign(*this, m); } template <typename U, size_t len> matrix ( U (&array)[len] ): matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager> >(ref_type(*this)) { COMPILE_TIME_ASSERT(NR*NC == len && len > 0); size_t idx = 0; for (long r = 0; r < NR; ++r) { for (long c = 0; c < NC; ++c) { data(r,c) = static_cast<T>(array[idx]); ++idx; } } } T& operator() ( long r, long c ) { DLIB_ASSERT(r < nr() && c < nc() && r >= 0 && c >= 0, "\tT& matrix::operator(r,c)" << "\n\tYou must give a valid row and column" << "\n\tr: " << r << "\n\tc: " << c << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tthis: " << this ); return data(r,c); } const T& operator() ( long r, long c ) const { DLIB_ASSERT(r < nr() && c < nc() && r >= 0 && c >= 0, "\tconst T& matrix::operator(r,c)" << "\n\tYou must give a valid row and column" << "\n\tr: " << r << "\n\tc: " << c << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tthis: " << this ); return data(r,c); } T& operator() ( long i ) { // You can only use this operator on column vectors. COMPILE_TIME_ASSERT(NC == 1 || NC == 0 || NR == 1 || NR == 0); DLIB_ASSERT(nc() == 1 || nr() == 1, "\tconst type matrix::operator(i)" << "\n\tYou can only use this operator on column or row vectors" << "\n\ti: " << i << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tthis: " << this ); DLIB_ASSERT( ((nc() == 1 && i < nr()) || (nr() == 1 && i < nc())) && i >= 0, "\tconst type matrix::operator(i)" << "\n\tYou must give a valid row/column number" << "\n\ti: " << i << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tthis: " << this ); return data(i); } const T& operator() ( long i ) const { // You can only use this operator on column vectors. COMPILE_TIME_ASSERT(NC == 1 || NC == 0 || NR == 1 || NR == 0); DLIB_ASSERT(nc() == 1 || nr() == 1, "\tconst type matrix::operator(i)" << "\n\tYou can only use this operator on column or row vectors" << "\n\ti: " << i << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tthis: " << this ); DLIB_ASSERT( ((nc() == 1 && i < nr()) || (nr() == 1 && i < nc())) && i >= 0, "\tconst type matrix::operator(i)" << "\n\tYou must give a valid row/column number" << "\n\ti: " << i << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tthis: " << this ); return data(i); } inline operator const type ( ) const { COMPILE_TIME_ASSERT(NC == 1 || NC == 0); COMPILE_TIME_ASSERT(NR == 1 || NR == 0); DLIB_ASSERT( nr() == 1 && nc() == 1 , "\tmatrix::operator const type" << "\n\tYou can only attempt to implicit convert a matrix to a scalar if" << "\n\tthe matrix is a 1x1 matrix" << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tthis: " << this ); return data(0); } void set_size ( long rows, long cols ) { DLIB_ASSERT( (NR == 0 || NR == rows) && ( NC == 0 || NC == cols) && rows >= 0 && cols >= 0, "\tvoid matrix::set_size(rows, cols)" << "\n\tYou have supplied conflicting matrix dimensions" << "\n\trows: " << rows << "\n\tcols: " << cols << "\n\tNR: " << NR << "\n\tNC: " << NC << "\n\tthis: " << this ); if (nr() != rows || nc() != cols) data.set_size(rows,cols); } void set_size ( long length ) { // This object you are trying to call set_size(length) on is not a column or // row vector. COMPILE_TIME_ASSERT(NR == 1 || NC == 1); DLIB_ASSERT( length >= 0, "\tvoid matrix::set_size(length)" << "\n\tlength must be at least 0" << "\n\tlength: " << length << "\n\tNR: " << NR << "\n\tNC: " << NC << "\n\tthis: " << this ); if (NR == 1) { DLIB_ASSERT(NC == 0 || NC == length, "\tvoid matrix::set_size(length)" << "\n\tSince this is a statically sized matrix length must equal NC" << "\n\tlength: " << length << "\n\tNR: " << NR << "\n\tNC: " << NC << "\n\tthis: " << this ); if (nc() != length) data.set_size(1,length); } else { DLIB_ASSERT(NR == 0 || NR == length, "\tvoid matrix::set_size(length)" << "\n\tSince this is a statically sized matrix length must equal NR" << "\n\tlength: " << length << "\n\tNR: " << NR << "\n\tNC: " << NC << "\n\tthis: " << this ); if (nr() != length) data.set_size(length,1); } } long nr ( ) const { return data.nr(); } long nc ( ) const { return data.nc(); } long size ( ) const { return data.nr()*data.nc(); } template <typename U, size_t len> matrix& operator= ( U (&array)[len] ) { COMPILE_TIME_ASSERT(NR*NC == len && len > 0); size_t idx = 0; for (long r = 0; r < NR; ++r) { for (long c = 0; c < NC; ++c) { data(r,c) = static_cast<T>(array[idx]); ++idx; } } return *this; } template <typename EXP> matrix& operator= ( const matrix_exp<EXP>& m ) { // You get an error on this line if the matrix you are trying to // assign m to is a statically sized matrix and m's dimensions don't // match that of *this. COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0); COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0); DLIB_ASSERT((NR == 0 || nr() == m.nr()) && (NC == 0 || nc() == m.nc()), "\tmatrix& matrix::operator=(const matrix_exp& m)" << "\n\tYou are trying to assign a dynamically sized matrix to a statically sized matrix with the wrong size" << "\n\tnr(): " << nr() << "\n\tnc(): " << nc() << "\n\tm.nr(): " << m.nr() << "\n\tm.nc(): " << m.nc() << "\n\tthis: " << this ); // You get an error on this line if the matrix m contains a type that isn't // the same as the type contained in the target matrix. COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true) || (is_matrix<typename EXP::type>::value == true)); if (m.destructively_aliases(*this) == false) { set_size(m.nr(),m.nc()); matrix_assign(*this, m); } else { // we have to use a temporary matrix object here because // *this is aliased inside the matrix_exp m somewhere. matrix temp; temp.set_size(m.nr(),m.nc()); matrix_assign(temp, m); temp.swap(*this); } return *this; } template <typename EXP> matrix& operator += ( const matrix_exp<EXP>& m ) { // The matrix you are trying to assign m to is a statically sized matrix and // m's dimensions don't match that of *this. COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0); COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0); DLIB_ASSERT(this->nr() == m.nr() && this->nc() == m.nc(), "\tmatrix& matrix::operator+=(const matrix_exp& m)" << "\n\tYou are trying to add a dynamically sized matrix to a statically sized matrix with the wrong size" << "\n\tthis->nr(): " << nr() << "\n\tthis->nc(): " << nc() << "\n\tm.nr(): " << m.nr() << "\n\tm.nc(): " << m.nc() << "\n\tthis: " << this ); COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true)); if (m.destructively_aliases(*this) == false) { matrix_assign(*this, m + *this); } else { // we have to use a temporary matrix_data object here because // this->data is aliased inside the matrix_exp m somewhere. matrix temp; temp.set_size(m.nr(),m.nc()); matrix_assign(temp, m + *this); temp.swap(*this); } return *this; } template <typename EXP> matrix& operator -= ( const matrix_exp<EXP>& m ) { // The matrix you are trying to assign m to is a statically sized matrix and // m's dimensions don't match that of *this. COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0); COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0); DLIB_ASSERT(this->nr() == m.nr() && this->nc() == m.nc(), "\tmatrix& matrix::operator-=(const matrix_exp& m)" << "\n\tYou are trying to subtract a dynamically sized matrix from a statically sized matrix with the wrong size" << "\n\tthis->nr(): " << nr() << "\n\tthis->nc(): " << nc() << "\n\tm.nr(): " << m.nr() << "\n\tm.nc(): " << m.nc() << "\n\tthis: " << this ); COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true)); if (m.destructively_aliases(*this) == false) { matrix_assign(*this, *this - m); } else { // we have to use a temporary matrix_data object here because // this->data is aliased inside the matrix_exp m somewhere. matrix temp; temp.set_size(m.nr(),m.nc()); matrix_assign(temp, *this - m); temp.swap(*this); } return *this; } matrix& operator += ( const matrix& m ) { const long size = m.nr()*m.nc(); for (long i = 0; i < size; ++i) data(i) += m.data(i); return *this; } matrix& operator -= ( const matrix& m ) { const long size = m.nr()*m.nc(); for (long i = 0; i < size; ++i) data(i) -= m.data(i); return *this; } matrix& operator *= ( const T& a ) { const long size = data.nr()*data.nc(); for (long i = 0; i < size; ++i) data(i) *= a; return *this; } matrix& operator /= ( const T& a ) { const long size = data.nr()*data.nc(); for (long i = 0; i < size; ++i) data(i) /= a; return *this; } matrix& operator= ( const matrix& m ) { if (this != &m) { set_size(m.nr(),m.nc()); const long size = m.nr()*m.nc(); for (long i = 0; i < size; ++i) data(i) = m.data(i); } return *this; } void swap ( matrix& item ) { data.swap(item.data); } private: matrix_data<T,NR,NC, mem_manager> data; }; // ---------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------- template < typename T, long NR, long NC, typename mm > void swap( matrix<T,NR,NC,mm>& a, matrix<T,NR,NC,mm>& b ) { a.swap(b); } template < typename T, long NR, long NC, typename mm > void serialize ( const matrix<T,NR,NC,mm>& item, std::ostream& out ) { try { serialize(item.nr(),out); serialize(item.nc(),out); for (long r = 0; r < item.nr(); ++r) { for (long c = 0; c < item.nc(); ++c) { serialize(item(r,c),out); } } } catch (serialization_error& e) { throw serialization_error(e.info + "\n while serializing dlib::matrix"); } } template < typename T, long NR, long NC, typename mm > void deserialize ( matrix<T,NR,NC,mm>& item, std::istream& in ) { try { long nr, nc; deserialize(nr,in); deserialize(nc,in); if (NR != 0 && nr != NR) throw serialization_error("Error while deserializing a dlib::matrix. Invalid rows"); if (NC != 0 && nc != NC) throw serialization_error("Error while deserializing a dlib::matrix. Invalid columns"); item.set_size(nr,nc); for (long r = 0; r < nr; ++r) { for (long c = 0; c < nc; ++c) { deserialize(item(r,c),in); } } } catch (serialization_error& e) { throw serialization_error(e.info + "\n while deserializing a dlib::matrix"); } } template < typename EXP > std::ostream& operator<< ( std::ostream& out, const matrix_exp<EXP>& m ) { using namespace std; const streamsize old = out.width(); // first figure out how wide we should make each field string::size_type w = 0; ostringstream sout; for (long r = 0; r < m.nr(); ++r) { for (long c = 0; c < m.nc(); ++c) { sout << m(r,c); w = std::max(sout.str().size(),w); sout.str(""); } } // now actually print it for (long r = 0; r < m.nr(); ++r) { for (long c = 0; c < m.nc(); ++c) { out.width(static_cast<streamsize>(w)); out << m(r,c) << " "; } out << "\n"; } out.width(old); return out; } // ---------------------------------------------------------------------------------------- } #ifdef _MSC_VER // put that warning back to its default setting #pragma warning(default : 4355) #endif #endif // DLIB_MATRIx_