diff options
Diffstat (limited to 'ml/dlib/dlib/matrix/matrix_utilities.h')
| -rw-r--r-- | ml/dlib/dlib/matrix/matrix_utilities.h | 4544 |
1 files changed, 4544 insertions, 0 deletions
diff --git a/ml/dlib/dlib/matrix/matrix_utilities.h b/ml/dlib/dlib/matrix/matrix_utilities.h new file mode 100644 index 000000000..0c5091a4b --- /dev/null +++ b/ml/dlib/dlib/matrix/matrix_utilities.h @@ -0,0 +1,4544 @@ +// Copyright (C) 2006 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_MATRIx_UTILITIES_ +#define DLIB_MATRIx_UTILITIES_ + +#include "matrix_utilities_abstract.h" +#include "matrix.h" +#include <cmath> +#include <complex> +#include <limits> +#include "../pixel.h" +#include "../stl_checked.h" +#include <vector> +#include <algorithm> +#include "../std_allocator.h" +#include "matrix_expressions.h" +#include "matrix_math_functions.h" +#include "matrix_op.h" +#include "../general_hash/random_hashing.h" +#include "matrix_mat.h" + + +namespace dlib +{ + +// ---------------------------------------------------------------------------------------- + + /*!A is_complex + This is a template that can be used to determine if a type is a specialization + of the std::complex template class. + + For example: + is_complex<float>::value == false + is_complex<std::complex<float> >::value == true + !*/ + + template <typename T> + struct is_complex { static const bool value = false; }; + + template <typename T> + struct is_complex<std::complex<T> > { static const bool value = true; }; + template <typename T> + struct is_complex<std::complex<T>& > { static const bool value = true; }; + template <typename T> + struct is_complex<const std::complex<T>& > { static const bool value = true; }; + template <typename T> + struct is_complex<const std::complex<T> > { static const bool value = true; }; + +// ---------------------------------------------------------------------------------------- + + template <typename EXP> + inline bool is_row_vector ( + const matrix_exp<EXP>& m + ) { return m.nr() == 1; } + + template <typename EXP> + inline bool is_col_vector ( + const matrix_exp<EXP>& m + ) { return m.nc() == 1; } + + template <typename EXP> + inline bool is_vector ( + const matrix_exp<EXP>& m + ) { return is_row_vector(m) || is_col_vector(m); } + +// ---------------------------------------------------------------------------------------- + + template <typename EXP> + inline bool is_finite ( + const matrix_exp<EXP>& m + ) + { + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + if (!is_finite(m(r,c))) + return false; + } + } + return true; + } + +// ---------------------------------------------------------------------------------------- + + namespace impl + { + template <typename T> + const T& magnitude (const T& item) { return item; } + template <typename T> + T magnitude (const std::complex<T>& item) { return std::norm(item); } + } + + template < + typename EXP + > + void find_min_and_max ( + const matrix_exp<EXP>& m, + typename EXP::type& min_val, + typename EXP::type& max_val + ) + { + DLIB_ASSERT(m.size() > 0, + "\ttype find_min_and_max(const matrix_exp& m, min_val, max_val)" + << "\n\tYou can't ask for the min and max of an empty matrix" + << "\n\tm.size(): " << m.size() + ); + typedef typename matrix_exp<EXP>::type type; + + min_val = m(0,0); + max_val = min_val; + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + type temp = m(r,c); + if (dlib::impl::magnitude(temp) > dlib::impl::magnitude(max_val)) + max_val = temp; + if (dlib::impl::magnitude(temp) < dlib::impl::magnitude(min_val)) + min_val = temp; + } + } + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + point max_point ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(m.size() > 0, + "\tpoint max_point(const matrix_exp& m)" + << "\n\tm can't be empty" + << "\n\tm.size(): " << m.size() + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + ); + typedef typename matrix_exp<EXP>::type type; + + point best_point(0,0); + type val = m(0,0); + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + type temp = m(r,c); + if (dlib::impl::magnitude(temp) > dlib::impl::magnitude(val)) + { + val = temp; + best_point = point(c,r); + } + } + } + return best_point; + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + point min_point ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(m.size() > 0, + "\tpoint min_point(const matrix_exp& m)" + << "\n\tm can't be empty" + << "\n\tm.size(): " << m.size() + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + ); + typedef typename matrix_exp<EXP>::type type; + + point best_point(0,0); + type val = m(0,0); + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + type temp = m(r,c); + if (dlib::impl::magnitude(temp) < dlib::impl::magnitude(val)) + { + val = temp; + best_point = point(c,r); + } + } + } + return best_point; + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + long index_of_max ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(m.size() > 0 && is_vector(m) == true, + "\tlong index_of_max(const matrix_exp& m)" + << "\n\tm must be a row or column matrix" + << "\n\tm.size(): " << m.size() + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + ); + typedef typename matrix_exp<EXP>::type type; + + type val = m(0); + long best_idx = 0; + for (long i = 1; i < m.size(); ++i) + { + type temp = m(i); + if (dlib::impl::magnitude(temp) > dlib::impl::magnitude(val)) + { + val = temp; + best_idx = i; + } + } + return best_idx; + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + long index_of_min ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(m.size() > 0 && is_vector(m), + "\tlong index_of_min(const matrix_exp& m)" + << "\n\tm must be a row or column matrix" + << "\n\tm.size(): " << m.size() + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + ); + typedef typename matrix_exp<EXP>::type type; + + type val = m(0); + long best_idx = 0; + for (long i = 1; i < m.size(); ++i) + { + type temp = m(i); + if (dlib::impl::magnitude(temp) < dlib::impl::magnitude(val)) + { + val = temp; + best_idx = i; + } + } + return best_idx; + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + const typename matrix_exp<EXP>::type max ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(m.size() > 0, + "\ttype max(const matrix_exp& m)" + << "\n\tYou can't ask for the max() of an empty matrix" + << "\n\tm.size(): " << m.size() + ); + typedef typename matrix_exp<EXP>::type type; + + type val = m(0,0); + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + type temp = m(r,c); + if (dlib::impl::magnitude(temp) > dlib::impl::magnitude(val)) + val = temp; + } + } + return val; + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + const typename matrix_exp<EXP>::type min ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(m.size() > 0, + "\ttype min(const matrix_exp& m)" + << "\n\tYou can't ask for the min() of an empty matrix" + << "\n\tm.size(): " << m.size() + ); + typedef typename matrix_exp<EXP>::type type; + + type val = m(0,0); + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + type temp = m(r,c); + if (dlib::impl::magnitude(temp) < dlib::impl::magnitude(val)) + val = temp; + } + } + return val; + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2> + struct op_binary_min : basic_op_mm<M1,M2> + { + op_binary_min( const M1& m1_, const M2& m2_) : basic_op_mm<M1,M2>(m1_,m2_){} + + typedef typename M1::type type; + typedef const type const_ret_type; + const static long cost = M1::cost + M2::cost + 1; + + const_ret_type apply ( long r, long c) const + { return std::min(this->m1(r,c),this->m2(r,c)); } + }; + + template < + typename EXP1, + typename EXP2 + > + inline const matrix_op<op_binary_min<EXP1,EXP2> > min_pointwise ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0); + COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NC == 0 || EXP2::NC == 0); + DLIB_ASSERT(a.nr() == b.nr() && + a.nc() == b.nc(), + "\t const matrix_exp min_pointwise(const matrix_exp& a, const matrix_exp& b)" + << "\n\ta.nr(): " << a.nr() + << "\n\ta.nc(): " << a.nc() + << "\n\tb.nr(): " << b.nr() + << "\n\tb.nc(): " << b.nc() + ); + typedef op_binary_min<EXP1,EXP2> op; + return matrix_op<op>(op(a.ref(),b.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2, typename M3> + struct op_min_pointwise3 : basic_op_mmm<M1,M2,M3> + { + op_min_pointwise3( const M1& m1_, const M2& m2_, const M3& m3_) : + basic_op_mmm<M1,M2,M3>(m1_,m2_,m3_){} + + typedef typename M1::type type; + typedef const typename M1::type const_ret_type; + const static long cost = M1::cost + M2::cost + M3::cost + 2; + + const_ret_type apply (long r, long c) const + { return std::min(this->m1(r,c),std::min(this->m2(r,c),this->m3(r,c))); } + }; + + template < + typename EXP1, + typename EXP2, + typename EXP3 + > + inline const matrix_op<op_min_pointwise3<EXP1,EXP2,EXP3> > + min_pointwise ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b, + const matrix_exp<EXP3>& c + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + COMPILE_TIME_ASSERT((is_same_type<typename EXP2::type,typename EXP3::type>::value == true)); + COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0); + COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NR == 0 || EXP2::NC == 0); + COMPILE_TIME_ASSERT(EXP2::NR == EXP3::NR || EXP2::NR == 0 || EXP3::NR == 0); + COMPILE_TIME_ASSERT(EXP2::NC == EXP3::NC || EXP2::NC == 0 || EXP3::NC == 0); + DLIB_ASSERT(a.nr() == b.nr() && + a.nc() == b.nc() && + b.nr() == c.nr() && + b.nc() == c.nc(), + "\tconst matrix_exp min_pointwise(a,b,c)" + << "\n\tYou can only make a do a pointwise min between equally sized matrices" + << "\n\ta.nr(): " << a.nr() + << "\n\ta.nc(): " << a.nc() + << "\n\tb.nr(): " << b.nr() + << "\n\tb.nc(): " << b.nc() + << "\n\tc.nr(): " << c.nr() + << "\n\tc.nc(): " << c.nc() + ); + + typedef op_min_pointwise3<EXP1,EXP2,EXP3> op; + return matrix_op<op>(op(a.ref(),b.ref(),c.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2> + struct op_binary_max : basic_op_mm<M1,M2> + { + op_binary_max( const M1& m1_, const M2& m2_) : basic_op_mm<M1,M2>(m1_,m2_){} + + typedef typename M1::type type; + typedef const type const_ret_type; + const static long cost = M1::cost + M2::cost + 1; + + const_ret_type apply ( long r, long c) const + { return std::max(this->m1(r,c),this->m2(r,c)); } + }; + + template < + typename EXP1, + typename EXP2 + > + inline const matrix_op<op_binary_max<EXP1,EXP2> > max_pointwise ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0); + COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NC == 0 || EXP2::NC == 0); + DLIB_ASSERT(a.nr() == b.nr() && + a.nc() == b.nc(), + "\t const matrix_exp max_pointwise(const matrix_exp& a, const matrix_exp& b)" + << "\n\ta.nr(): " << a.nr() + << "\n\ta.nc(): " << a.nc() + << "\n\tb.nr(): " << b.nr() + << "\n\tb.nc(): " << b.nc() + ); + typedef op_binary_max<EXP1,EXP2> op; + return matrix_op<op>(op(a.ref(),b.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2, typename M3> + struct op_max_pointwise3 : basic_op_mmm<M1,M2,M3> + { + op_max_pointwise3( const M1& m1_, const M2& m2_, const M3& m3_) : + basic_op_mmm<M1,M2,M3>(m1_,m2_,m3_){} + + typedef typename M1::type type; + typedef const typename M1::type const_ret_type; + const static long cost = M1::cost + M2::cost + M3::cost + 2; + + const_ret_type apply (long r, long c) const + { return std::max(this->m1(r,c),std::max(this->m2(r,c),this->m3(r,c))); } + }; + + template < + typename EXP1, + typename EXP2, + typename EXP3 + > + inline const matrix_op<op_max_pointwise3<EXP1,EXP2,EXP3> > + max_pointwise ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b, + const matrix_exp<EXP3>& c + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + COMPILE_TIME_ASSERT((is_same_type<typename EXP2::type,typename EXP3::type>::value == true)); + COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0); + COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NR == 0 || EXP2::NC == 0); + COMPILE_TIME_ASSERT(EXP2::NR == EXP3::NR || EXP2::NR == 0 || EXP3::NR == 0); + COMPILE_TIME_ASSERT(EXP2::NC == EXP3::NC || EXP2::NC == 0 || EXP3::NC == 0); + DLIB_ASSERT(a.nr() == b.nr() && + a.nc() == b.nc() && + b.nr() == c.nr() && + b.nc() == c.nc(), + "\tconst matrix_exp max_pointwise(a,b,c)" + << "\n\tYou can only make a do a pointwise max between equally sized matrices" + << "\n\ta.nr(): " << a.nr() + << "\n\ta.nc(): " << a.nc() + << "\n\tb.nr(): " << b.nr() + << "\n\tb.nc(): " << b.nc() + << "\n\tc.nr(): " << c.nr() + << "\n\tc.nc(): " << c.nc() + ); + + typedef op_max_pointwise3<EXP1,EXP2,EXP3> op; + return matrix_op<op>(op(a.ref(),b.ref(),c.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + typename enable_if_c<std::numeric_limits<typename EXP::type>::is_integer, double>::type length ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(is_vector(m) == true, + "\ttype length(const matrix_exp& m)" + << "\n\tm must be a row or column vector" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + ); + + return std::sqrt(static_cast<double>(sum(squared(m)))); + } + + template < + typename EXP + > + typename disable_if_c<std::numeric_limits<typename EXP::type>::is_integer, const typename EXP::type>::type length ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(is_vector(m) == true, + "\ttype length(const matrix_exp& m)" + << "\n\tm must be a row or column vector" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + ); + return std::sqrt(sum(squared(m))); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + const typename matrix_exp<EXP>::type length_squared ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(is_vector(m) == true, + "\ttype length_squared(const matrix_exp& m)" + << "\n\tm must be a row or column vector" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + ); + return sum(squared(m)); + } + +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_trans + { + op_trans( const M& m_) : m(m_){} + + const M& m; + + const static long cost = M::cost; + const static long NR = M::NC; + const static long NC = M::NR; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + + const_ret_type apply (long r, long c) const { return m(c,r); } + + long nr () const { return m.nc(); } + long nc () const { return m.nr(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + + }; + + template < + typename M + > + const matrix_op<op_trans<M> > trans ( + const matrix_exp<M>& m + ) + { + typedef op_trans<M> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + +// don't to anything at all for diagonal matrices + template < + typename M + > + const matrix_diag_exp<M>& trans ( + const matrix_diag_exp<M>& m + ) + { + return m; + } + +// ---------------------------------------------------------------------------------------- + +// I introduced this struct because it avoids an inane compiler warning from gcc + template <typename EXP> + struct is_not_ct_vector{ static const bool value = (EXP::NR != 1 && EXP::NC != 1); }; + + template < + typename EXP1, + typename EXP2 + > + typename enable_if_c<(is_not_ct_vector<EXP1>::value) || (is_not_ct_vector<EXP2>::value), + typename EXP1::type>::type + dot ( + const matrix_exp<EXP1>& m1, + const matrix_exp<EXP2>& m2 + ) + { + // You are getting an error on this line because you are trying to + // compute the dot product between two matrices that aren't both vectors (i.e. + // they aren't column or row matrices). + COMPILE_TIME_ASSERT(EXP1::NR*EXP1::NC == 0 || + EXP2::NR*EXP2::NC == 0); + + DLIB_ASSERT(is_vector(m1) && is_vector(m2) && m1.size() == m2.size() && + m1.size() > 0, + "\t type dot(const matrix_exp& m1, const matrix_exp& m2)" + << "\n\t You can only compute the dot product between non-empty vectors of equal length." + << "\n\t is_vector(m1): " << is_vector(m1) + << "\n\t is_vector(m2): " << is_vector(m2) + << "\n\t m1.size(): " << m1.size() + << "\n\t m2.size(): " << m2.size() + ); + + if (is_col_vector(m1) && is_col_vector(m2)) return (trans(m1)*m2)(0); + if (is_col_vector(m1) && is_row_vector(m2)) return (m2*m1)(0); + if (is_row_vector(m1) && is_col_vector(m2)) return (m1*m2)(0); + + //if (is_row_vector(m1) && is_row_vector(m2)) + return (m1*trans(m2))(0); + } + + template < typename EXP1, typename EXP2 > + typename enable_if_c<EXP1::NR == 1 && EXP2::NR == 1 && EXP1::NC != 1 && EXP2::NC != 1, typename EXP1::type>::type + dot ( const matrix_exp<EXP1>& m1, const matrix_exp<EXP2>& m2) + { + DLIB_ASSERT(m1.size() == m2.size(), + "\t type dot(const matrix_exp& m1, const matrix_exp& m2)" + << "\n\t You can only compute the dot product between vectors of equal length" + << "\n\t m1.size(): " << m1.size() + << "\n\t m2.size(): " << m2.size() + ); + + return m1*trans(m2); + } + + template < typename EXP1, typename EXP2 > + typename enable_if_c<EXP1::NR == 1 && EXP2::NC == 1 && EXP1::NC != 1 && EXP2::NR != 1, typename EXP1::type>::type + dot ( const matrix_exp<EXP1>& m1, const matrix_exp<EXP2>& m2) + { + DLIB_ASSERT(m1.size() == m2.size(), + "\t type dot(const matrix_exp& m1, const matrix_exp& m2)" + << "\n\t You can only compute the dot product between vectors of equal length" + << "\n\t m1.size(): " << m1.size() + << "\n\t m2.size(): " << m2.size() + ); + + return m1*m2; + } + + template < typename EXP1, typename EXP2 > + typename enable_if_c<EXP1::NC == 1 && EXP2::NR == 1 && EXP1::NR != 1 && EXP2::NC != 1, typename EXP1::type>::type + dot ( const matrix_exp<EXP1>& m1, const matrix_exp<EXP2>& m2) + { + DLIB_ASSERT(m1.size() == m2.size(), + "\t type dot(const matrix_exp& m1, const matrix_exp& m2)" + << "\n\t You can only compute the dot product between vectors of equal length" + << "\n\t m1.size(): " << m1.size() + << "\n\t m2.size(): " << m2.size() + ); + + return m2*m1; + } + + template < typename EXP1, typename EXP2 > + typename enable_if_c<EXP1::NC == 1 && EXP2::NC == 1 && EXP1::NR != 1 && EXP2::NR != 1, typename EXP1::type>::type + dot ( const matrix_exp<EXP1>& m1, const matrix_exp<EXP2>& m2) + { + DLIB_ASSERT(m1.size() == m2.size(), + "\t type dot(const matrix_exp& m1, const matrix_exp& m2)" + << "\n\t You can only compute the dot product between vectors of equal length" + << "\n\t m1.size(): " << m1.size() + << "\n\t m2.size(): " << m2.size() + ); + + return trans(m1)*m2; + } + + template < typename EXP1, typename EXP2 > + typename enable_if_c<(EXP1::NC*EXP1::NR == 1) || (EXP2::NC*EXP2::NR == 1), typename EXP1::type>::type + dot ( const matrix_exp<EXP1>& m1, const matrix_exp<EXP2>& m2) + { + DLIB_ASSERT(m1.size() == m2.size(), + "\t type dot(const matrix_exp& m1, const matrix_exp& m2)" + << "\n\t You can only compute the dot product between vectors of equal length" + << "\n\t m1.size(): " << m1.size() + << "\n\t m2.size(): " << m2.size() + ); + + return m1(0)*m2(0); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M, long R, long C> + struct op_removerc + { + op_removerc( const M& m_) : m(m_){} + + const M& m; + + const static long cost = M::cost+2; + const static long NR = (M::NR==0) ? 0 : (M::NR - 1); + const static long NC = (M::NC==0) ? 0 : (M::NC - 1); + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply (long r, long c) const + { + if (r < R) + { + if (c < C) + return m(r,c); + else + return m(r,c+1); + } + else + { + if (c < C) + return m(r+1,c); + else + return m(r+1,c+1); + } + } + + long nr () const { return m.nr() - 1; } + long nc () const { return m.nc() - 1; } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template <typename M> + struct op_removerc2 + { + op_removerc2( const M& m_, const long R_, const long C_) : m(m_), R(R_), C(C_){} + const M& m; + const long R; + const long C; + + const static long cost = M::cost+2; + const static long NR = (M::NR==0) ? 0 : (M::NR - 1); + const static long NC = (M::NC==0) ? 0 : (M::NC - 1); + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply (long r, long c) const + { + if (r < R) + { + if (c < C) + return m(r,c); + else + return m(r,c+1); + } + else + { + if (c < C) + return m(r+1,c); + else + return m(r+1,c+1); + } + } + + long nr () const { return m.nr() - 1; } + long nc () const { return m.nc() - 1; } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template < + long R, + long C, + typename EXP + > + const matrix_op<op_removerc<EXP,R,C> > removerc ( + const matrix_exp<EXP>& m + ) + { + // you can't remove a row from a matrix with only one row + COMPILE_TIME_ASSERT((EXP::NR > R && R >= 0) || EXP::NR == 0); + // you can't remove a column from a matrix with only one column + COMPILE_TIME_ASSERT((EXP::NC > C && C >= 0) || EXP::NR == 0); + DLIB_ASSERT(m.nr() > R && R >= 0 && m.nc() > C && C >= 0, + "\tconst matrix_exp removerc<R,C>(const matrix_exp& m)" + << "\n\tYou can't remove a row/column from a matrix if it doesn't have that row/column" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\tR: " << R + << "\n\tC: " << C + ); + typedef op_removerc<EXP,R,C> op; + return matrix_op<op>(op(m.ref())); + } + + template < + typename EXP + > + const matrix_op<op_removerc2<EXP> > removerc ( + const matrix_exp<EXP>& m, + long R, + long C + ) + { + DLIB_ASSERT(m.nr() > R && R >= 0 && m.nc() > C && C >= 0, + "\tconst matrix_exp removerc(const matrix_exp& m,R,C)" + << "\n\tYou can't remove a row/column from a matrix if it doesn't have that row/column" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\tR: " << R + << "\n\tC: " << C + ); + typedef op_removerc2<EXP> op; + return matrix_op<op>(op(m.ref(),R,C)); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M, long C> + struct op_remove_col + { + op_remove_col( const M& m_) : m(m_){} + const M& m; + + const static long cost = M::cost+2; + const static long NR = M::NR; + const static long NC = (M::NC==0) ? 0 : (M::NC - 1); + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply ( long r, long c) const + { + if (c < C) + { + return m(r,c); + } + else + { + return m(r,c+1); + } + } + + long nr () const { return m.nr(); } + long nc () const { return m.nc() - 1; } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template <typename M> + struct op_remove_col2 + { + op_remove_col2( const M& m_, const long C_) : m(m_), C(C_){} + const M& m; + const long C; + + const static long cost = M::cost+2; + const static long NR = M::NR; + const static long NC = (M::NC==0) ? 0 : (M::NC - 1); + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply ( long r, long c) const + { + if (c < C) + { + return m(r,c); + } + else + { + return m(r,c+1); + } + } + + long nr () const { return m.nr(); } + long nc () const { return m.nc() - 1; } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template < + long C, + typename EXP + > + const matrix_op<op_remove_col<EXP, C> > remove_col ( + const matrix_exp<EXP>& m + ) + { + // You can't remove the given column from the matrix because the matrix doesn't + // have a column with that index. + COMPILE_TIME_ASSERT((EXP::NC > C && C >= 0) || EXP::NC == 0); + DLIB_ASSERT(m.nc() > C && C >= 0 , + "\tconst matrix_exp remove_col<C>(const matrix_exp& m)" + << "\n\tYou can't remove a col from a matrix if it doesn't have it" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\tC: " << C + ); + typedef op_remove_col<EXP,C> op; + return matrix_op<op>(op(m.ref())); + } + + template < + typename EXP + > + const matrix_op<op_remove_col2<EXP> > remove_col ( + const matrix_exp<EXP>& m, + long C + ) + { + DLIB_ASSERT(m.nc() > C && C >= 0 , + "\tconst matrix_exp remove_col(const matrix_exp& m,C)" + << "\n\tYou can't remove a col from a matrix if it doesn't have it" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\tC: " << C + ); + typedef op_remove_col2<EXP> op; + return matrix_op<op>(op(m.ref(),C)); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M, long R> + struct op_remove_row + { + op_remove_row( const M& m_) : m(m_){} + const M& m; + + const static long cost = M::cost+2; + const static long NR = (M::NR==0) ? 0 : (M::NR - 1); + const static long NC = M::NC; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply ( long r, long c) const + { + if (r < R) + { + return m(r,c); + } + else + { + return m(r+1,c); + } + } + + long nr () const { return m.nr() - 1; } + long nc () const { return m.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template <typename M> + struct op_remove_row2 + { + op_remove_row2( const M& m_, const long R_) : m(m_), R(R_){} + const M& m; + const long R; + + const static long cost = M::cost+2; + const static long NR = (M::NR==0) ? 0 : (M::NR - 1); + const static long NC = M::NC; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply ( long r, long c) const + { + if (r < R) + { + return m(r,c); + } + else + { + return m(r+1,c); + } + } + + long nr () const { return m.nr() - 1; } + long nc () const { return m.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template < + long R, + typename EXP + > + const matrix_op<op_remove_row<EXP,R> > remove_row ( + const matrix_exp<EXP>& m + ) + { + // You can't remove the given row from the matrix because the matrix doesn't + // have a row with that index. + COMPILE_TIME_ASSERT((EXP::NR > R && R >= 0) || EXP::NR == 0); + DLIB_ASSERT(m.nr() > R && R >= 0, + "\tconst matrix_exp remove_row<R>(const matrix_exp& m)" + << "\n\tYou can't remove a row from a matrix if it doesn't have it" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\tR: " << R + ); + typedef op_remove_row<EXP,R> op; + return matrix_op<op>(op(m.ref())); + } + + template < + typename EXP + > + const matrix_op<op_remove_row2<EXP> > remove_row ( + const matrix_exp<EXP>& m, + long R + ) + { + DLIB_ASSERT(m.nr() > R && R >= 0, + "\tconst matrix_exp remove_row(const matrix_exp& m, long R)" + << "\n\tYou can't remove a row from a matrix if it doesn't have it" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\tR: " << R + ); + typedef op_remove_row2<EXP> op; + return matrix_op<op>(op(m.ref(),R)); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_diagm + { + op_diagm( const M& m_) : m(m_){} + const M& m; + + const static long cost = M::cost+2; + const static long N = M::NC*M::NR; + const static long NR = N; + const static long NC = N; + typedef typename M::type type; + typedef const typename M::type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply ( long r, long c) const + { + if (r==c) + return m(r); + else + return 0; + } + + long nr () const { return (m.nr()>m.nc())? m.nr():m.nc(); } + long nc () const { return (m.nr()>m.nc())? m.nr():m.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template < + typename EXP + > + const matrix_diag_op<op_diagm<EXP> > diagm ( + const matrix_exp<EXP>& m + ) + { + // You can only make a diagonal matrix out of a row or column vector + COMPILE_TIME_ASSERT(EXP::NR == 0 || EXP::NR == 1 || EXP::NC == 1 || EXP::NC == 0); + DLIB_ASSERT(is_vector(m), + "\tconst matrix_exp diagm(const matrix_exp& m)" + << "\n\tYou can only apply diagm() to a row or column matrix" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + ); + typedef op_diagm<EXP> op; + return matrix_diag_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2> + struct op_diagm_mult : basic_op_mm<M1,M2> + { + op_diagm_mult( const M1& m1_, const M2& m2_) : basic_op_mm<M1,M2>(m1_,m2_){} + + typedef typename M1::type type; + typedef const type const_ret_type; + const static long cost = M1::cost + M2::cost + 1; + + const_ret_type apply ( long r, long c) const + { + if (r == c) + return this->m1(r,c)*this->m2(r,c); + else + return 0; + } + }; + + template < + typename EXP1, + typename EXP2 + > + inline const matrix_diag_op<op_diagm_mult<EXP1,EXP2> > operator* ( + const matrix_diag_exp<EXP1>& a, + const matrix_diag_exp<EXP2>& b + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type, typename EXP2::type>::value)); + COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0); + COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NC == 0 || EXP2::NC == 0); + DLIB_ASSERT(a.nr() == b.nr() && + a.nc() == b.nc(), + "\tconst matrix_exp operator(const matrix_diag_exp& a, const matrix_diag_exp& b)" + << "\n\tYou can only multiply diagonal matrices together if they are the same size" + << "\n\ta.nr(): " << a.nr() + << "\n\ta.nc(): " << a.nc() + << "\n\tb.nr(): " << b.nr() + << "\n\tb.nc(): " << b.nc() + ); + typedef op_diagm_mult<EXP1,EXP2> op; + return matrix_diag_op<op>(op(a.ref(),b.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_diag + { + op_diag( const M& m_) : m(m_){} + const M& m; + + const static long cost = M::cost; + const static long NR = tmin<M::NR,M::NC>::value; + const static long NC = 1; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply ( long r, long ) const { return m(r,r); } + + long nr () const { return std::min(m.nc(),m.nr()); } + long nc () const { return 1; } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template < + typename EXP + > + const matrix_op<op_diag<EXP> > diag ( + const matrix_exp<EXP>& m + ) + { + typedef op_diag<EXP> op; + return matrix_op<op>(op(m.ref())); + } + + template <typename EXP> + struct diag_exp + { + typedef matrix_op<op_diag<EXP> > type; + }; + +// ---------------------------------------------------------------------------------------- + + template <typename M, typename target_type> + struct op_cast + { + op_cast( const M& m_) : m(m_){} + const M& m; + + const static long cost = M::cost+2; + const static long NR = M::NR; + const static long NC = M::NC; + typedef target_type type; + typedef const target_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply ( long r, long c) const { return static_cast<target_type>(m(r,c)); } + + long nr () const { return m.nr(); } + long nc () const { return m.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.destructively_aliases(item); } + }; + + template < + typename target_type, + typename EXP + > + const matrix_op<op_cast<EXP, target_type> > matrix_cast ( + const matrix_exp<EXP>& m + ) + { + typedef op_cast<EXP, target_type> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- + + namespace impl + { + template <typename type, typename S> + inline type lessthan(const type& val, const S& s) + { + if (val < s) + return 1; + else + return 0; + } + + } + DLIB_DEFINE_OP_MS(op_lessthan, impl::lessthan, 1); + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_lessthan<EXP,S> > >::type operator< ( + const matrix_exp<EXP>& m, + const S& s + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_lessthan<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_lessthan<EXP,S> > >::type operator> ( + const S& s, + const matrix_exp<EXP>& m + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_lessthan<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + +// ---------------------------------------------------------------------------------------- + + namespace impl + { + template <typename type, typename S> + inline type lessthan_eq(const type& val, const S& s) + { + if (val <= s) + return 1; + else + return 0; + } + + } + DLIB_DEFINE_OP_MS(op_lessthan_eq, impl::lessthan_eq, 1); + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_lessthan_eq<EXP,S> > >::type operator<= ( + const matrix_exp<EXP>& m, + const S& s + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_lessthan_eq<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_lessthan_eq<EXP,S> > >::type operator>= ( + const S& s, + const matrix_exp<EXP>& m + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_lessthan_eq<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + +// ---------------------------------------------------------------------------------------- + + namespace impl + { + template <typename type, typename S> + inline type greaterthan(const type& val, const S& s) + { + if (val > s) + return 1; + else + return 0; + } + + } + DLIB_DEFINE_OP_MS(op_greaterthan, impl::greaterthan, 1); + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_greaterthan<EXP,S> > >::type operator> ( + const matrix_exp<EXP>& m, + const S& s + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_greaterthan<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_greaterthan<EXP,S> > >::type operator< ( + const S& s, + const matrix_exp<EXP>& m + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_greaterthan<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + +// ---------------------------------------------------------------------------------------- + + namespace impl + { + template <typename type, typename S> + inline type greaterthan_eq(const type& val, const S& s) + { + if (val >= s) + return 1; + else + return 0; + } + + } + DLIB_DEFINE_OP_MS(op_greaterthan_eq, impl::greaterthan_eq, 1); + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_greaterthan_eq<EXP,S> > >::type operator>= ( + const matrix_exp<EXP>& m, + const S& s + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_greaterthan_eq<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_greaterthan_eq<EXP,S> > >::type operator<= ( + const S& s, + const matrix_exp<EXP>& m + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_greaterthan_eq<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + +// ---------------------------------------------------------------------------------------- + + namespace impl + { + template <typename type, typename S> + inline type equal_to(const type& val, const S& s) + { + if (val == s) + return 1; + else + return 0; + } + + } + DLIB_DEFINE_OP_MS(op_equal_to, impl::equal_to, 1); + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_equal_to<EXP,S> > >::type operator== ( + const matrix_exp<EXP>& m, + const S& s + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT( is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_equal_to<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_equal_to<EXP,S> > >::type operator== ( + const S& s, + const matrix_exp<EXP>& m + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT( is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_equal_to<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + +// ---------------------------------------------------------------------------------------- + + namespace impl + { + template <typename type, typename S> + inline type not_equal_to(const type& val, const S& s) + { + if (val != s) + return 1; + else + return 0; + } + + } + DLIB_DEFINE_OP_MS(op_not_equal_to, impl::not_equal_to, 1); + + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_not_equal_to<EXP,S> > >::type operator!= ( + const matrix_exp<EXP>& m, + const S& s + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_not_equal_to<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + + template < + typename EXP, + typename S + > + const typename enable_if<is_built_in_scalar_type<S>, matrix_op<op_not_equal_to<EXP,S> > >::type operator!= ( + const S& s, + const matrix_exp<EXP>& m + ) + { + // you can only use this relational operator with the built in scalar types like + // long, float, etc. + COMPILE_TIME_ASSERT(is_built_in_scalar_type<typename EXP::type>::value); + + typedef op_not_equal_to<EXP,S> op; + return matrix_op<op>(op(m.ref(),s)); + } + +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- + + template < + typename T, + long NR, + long NC, + typename MM, + typename U, + typename L + > + typename disable_if<is_matrix<U>,void>::type set_all_elements ( + matrix<T,NR,NC,MM,L>& m, + const U& value + ) + { + // The value you are trying to assign to each element of the m matrix + // doesn't have the appropriate type. + COMPILE_TIME_ASSERT(is_matrix<T>::value == is_matrix<U>::value); + + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + m(r,c) = static_cast<T>(value); + } + } + } + +// ---------------------------------------------------------------------------------------- + + template < + typename T, + long NR, + long NC, + typename MM, + typename U, + typename L + > + typename enable_if<is_matrix<U>,void>::type set_all_elements ( + matrix<T,NR,NC,MM,L>& m, + const U& value + ) + { + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + m(r,c) = value; + } + } + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + inline const typename matrix_exp<EXP>::matrix_type tmp ( + const matrix_exp<EXP>& m + ) + { + return typename matrix_exp<EXP>::matrix_type (m); + } + +// ---------------------------------------------------------------------------------------- + + template <typename EXP> + constexpr bool is_row_major ( + const matrix_exp<EXP>& + ) + { + return is_same_type<typename EXP::layout_type,row_major_layout>::value; + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + const typename lazy_disable_if<is_matrix<typename EXP::type>, EXP>::type sum ( + const matrix_exp<EXP>& m + ) + { + typedef typename matrix_exp<EXP>::type type; + + type val = 0; + if (is_row_major(m)) + { + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + val += m(r,c); + } + } + } + else + { + for (long c = 0; c < m.nc(); ++c) + { + for (long r = 0; r < m.nr(); ++r) + { + val += m(r,c); + } + } + } + return val; + } + + template < + typename EXP + > + const typename lazy_enable_if<is_matrix<typename EXP::type>, EXP>::type sum ( + const matrix_exp<EXP>& m + ) + { + typedef typename matrix_exp<EXP>::type type; + + type val; + if (m.size() > 0) + val.set_size(m(0,0).nr(),m(0,0).nc()); + set_all_elements(val,0); + + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + val += m(r,c); + } + } + return val; + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_sumr + { + op_sumr(const M& m_) : m(m_) {} + const M& m; + + const static long cost = M::cost+10; + const static long NR = 1; + const static long NC = M::NC; + typedef typename M::type type; + typedef const typename M::type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply ( long , long c) const + { + type temp = m(0,c); + for (long r = 1; r < m.nr(); ++r) + temp += m(r,c); + return temp; + } + + long nr () const { return 1; } + long nc () const { return m.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template < + typename EXP + > + const matrix_op<op_sumr<EXP> > sum_rows ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(m.size() > 0 , + "\tconst matrix_exp sum_rows(m)" + << "\n\t The matrix can't be empty" + << "\n\t m.size(): " << m.size() + ); + typedef op_sumr<EXP> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_sumc + { + op_sumc(const M& m_) : m(m_) {} + const M& m; + + const static long cost = M::cost + 10; + const static long NR = M::NR; + const static long NC = 1; + typedef typename M::type type; + typedef const typename M::type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply ( long r, long ) const + { + type temp = m(r,0); + for (long c = 1; c < m.nc(); ++c) + temp += m(r,c); + return temp; + } + + long nr () const { return m.nr(); } + long nc () const { return 1; } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template < + typename EXP + > + const matrix_op<op_sumc<EXP> > sum_cols ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(m.size() > 0 , + "\tconst matrix_exp sum_cols(m)" + << "\n\t The matrix can't be empty" + << "\n\t m.size(): " << m.size() + ); + typedef op_sumc<EXP> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + inline const typename disable_if<is_complex<typename EXP::type>, typename matrix_exp<EXP>::type>::type mean ( + const matrix_exp<EXP>& m + ) + { + return sum(m)/(m.nr()*m.nc()); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + inline const typename enable_if<is_complex<typename EXP::type>, typename matrix_exp<EXP>::type>::type mean ( + const matrix_exp<EXP>& m + ) + { + typedef typename EXP::type::value_type type; + return sum(m)/(type)(m.nr()*m.nc()); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + const typename matrix_exp<EXP>::type variance ( + const matrix_exp<EXP>& m + ) + { + using std::pow; + using dlib::pow; + const typename matrix_exp<EXP>::type avg = mean(m); + + typedef typename matrix_exp<EXP>::type type; + + type val; + val = 0; + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + val += pow(m(r,c) - avg,2); + } + } + + if (m.nr() * m.nc() <= 1) + { + return val; + } + else + { + // Note, for some reason, in gcc 4.1 performing this division using a + // double instead of a long value avoids a segmentation fault. That is, + // using 1.0 instead of 1 does the trick. + return val/(m.nr()*m.nc() - 1.0); + } + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + const typename matrix_exp<EXP>::type stddev ( + const matrix_exp<EXP>& m + ) + { + using std::sqrt; + using dlib::sqrt; + return sqrt(variance(m)); + } + +// ---------------------------------------------------------------------------------------- + +// this is a workaround for a bug in visual studio 7.1 + template <typename EXP> + struct visual_studio_sucks_cov_helper + { + typedef typename EXP::type inner_type; + typedef matrix<typename inner_type::type, inner_type::NR, inner_type::NR, typename EXP::mem_manager_type> type; + }; + + template < + typename EXP + > + const typename visual_studio_sucks_cov_helper<EXP>::type covariance ( + const matrix_exp<EXP>& m + ) + { + // perform static checks to make sure m is a column vector + COMPILE_TIME_ASSERT(EXP::NR == 0 || EXP::NR > 1); + COMPILE_TIME_ASSERT(EXP::NC == 1 || EXP::NC == 0); + + // perform static checks to make sure the matrices contained in m are column vectors + COMPILE_TIME_ASSERT(EXP::type::NC == 1 || EXP::type::NC == 0 ); + + DLIB_ASSERT(m.size() > 1 && is_col_vector(m), + "\tconst matrix covariance(const matrix_exp& m)" + << "\n\tYou can only apply covariance() to a column matrix" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + ); +#ifdef ENABLE_ASSERTS + for (long i = 0; i < m.nr(); ++i) + { + DLIB_ASSERT(m(0).size() == m(i).size() && m(i).size() > 0 && is_col_vector(m(i)), + "\tconst matrix covariance(const matrix_exp& m)" + << "\n\tYou can only apply covariance() to a column matrix of column matrices" + << "\n\tm(0).size(): " << m(0).size() + << "\n\tm(i).size(): " << m(i).size() + << "\n\tis_col_vector(m(i)): " << (is_col_vector(m(i)) ? "true" : "false") + << "\n\ti: " << i + ); + } +#endif + + // now perform the actual calculation of the covariance matrix. + typename visual_studio_sucks_cov_helper<EXP>::type cov(m(0).nr(),m(0).nr()); + set_all_elements(cov,0); + + const typename EXP::type avg = mean(m); + + for (long r = 0; r < m.nr(); ++r) + { + cov += (m(r) - avg)*trans(m(r) - avg); + } + + cov *= 1.0 / (m.nr() - 1.0); + return cov; + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + const typename matrix_exp<EXP>::type prod ( + const matrix_exp<EXP>& m + ) + { + typedef typename matrix_exp<EXP>::type type; + + type val = 1; + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + val *= m(r,c); + } + } + return val; + } + +// ---------------------------------------------------------------------------------------- + + template < + typename T + > + struct op_uniform_matrix_3 : does_not_alias + { + op_uniform_matrix_3(const long& rows_, const long& cols_, const T& val_ ) : + rows(rows_), cols(cols_), val(val_) {} + + const long rows; + const long cols; + const T val; + + const static long cost = 1; + const static long NR = 0; + const static long NC = 0; + typedef default_memory_manager mem_manager_type; + typedef row_major_layout layout_type; + typedef T type; + typedef const T& const_ret_type; + const_ret_type apply (long, long ) const { return val; } + + long nr() const { return rows; } + long nc() const { return cols; } + }; + + template < + typename T + > + const matrix_op<op_uniform_matrix_3<T> > uniform_matrix ( + long nr, + long nc, + const T& val + ) + { + DLIB_ASSERT(nr >= 0 && nc >= 0, + "\tconst matrix_exp uniform_matrix<T>(nr, nc, val)" + << "\n\tnr and nc have to be bigger than 0" + << "\n\tnr: " << nr + << "\n\tnc: " << nc + ); + typedef op_uniform_matrix_3<T> op; + return matrix_op<op>(op(nr, nc, val)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename T + > + const matrix_op<op_uniform_matrix_3<T> > zeros_matrix ( + long nr, + long nc + ) + { + DLIB_ASSERT(nr >= 0 && nc >= 0, + "\tconst matrix_exp zeros_matrix<T>(nr, nc)" + << "\n\tnr and nc have to be >= 0" + << "\n\tnr: " << nr + << "\n\tnc: " << nc + ); + typedef op_uniform_matrix_3<T> op; + return matrix_op<op>(op(nr, nc, 0)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + const matrix_op<op_uniform_matrix_3<typename EXP::type> > zeros_matrix ( + const matrix_exp<EXP>& mat + ) + { + DLIB_ASSERT(mat.nr() >= 0 && mat.nc() >= 0, + "\tconst matrix_exp zeros_matrix(mat)" + << "\n\t nr and nc have to be >= 0" + << "\n\t mat.nr(): " << mat.nr() + << "\n\t mat.nc(): " << mat.nc() + ); + typedef typename EXP::type T; + typedef op_uniform_matrix_3<T> op; + return matrix_op<op>(op(mat.nr(), mat.nc(), 0)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename T + > + const matrix_op<op_uniform_matrix_3<T> > ones_matrix ( + long nr, + long nc + ) + { + DLIB_ASSERT(nr >= 0 && nc >= 0, + "\tconst matrix_exp ones_matrix<T>(nr, nc)" + << "\n\tnr and nc have to be >= 0" + << "\n\tnr: " << nr + << "\n\tnc: " << nc + ); + typedef op_uniform_matrix_3<T> op; + return matrix_op<op>(op(nr, nc, 1)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP + > + const matrix_op<op_uniform_matrix_3<typename EXP::type> > ones_matrix ( + const matrix_exp<EXP>& mat + ) + { + DLIB_ASSERT(mat.nr() >= 0 && mat.nc() >= 0, + "\tconst matrix_exp ones_matrix(mat)" + << "\n\t nr and nc have to be >= 0" + << "\n\t mat.nr(): " << mat.nr() + << "\n\t mat.nc(): " << mat.nc() + ); + typedef typename EXP::type T; + typedef op_uniform_matrix_3<T> op; + return matrix_op<op>(op(mat.nr(), mat.nc(), 1)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename T, + long NR_, + long NC_ + > + struct op_uniform_matrix_2 : does_not_alias + { + op_uniform_matrix_2( const T& val_ ) : val(val_) {} + const T val; + + const static long cost = 1; + const static long NR = NR_; + const static long NC = NC_; + typedef default_memory_manager mem_manager_type; + typedef row_major_layout layout_type; + typedef T type; + typedef const T& const_ret_type; + + const_ret_type apply (long , long ) const { return val; } + + long nr() const { return NR; } + long nc() const { return NC; } + }; + + template < + typename T, + long NR, + long NC + > + const matrix_op<op_uniform_matrix_2<T,NR,NC> > uniform_matrix ( + const T& val + ) + { + COMPILE_TIME_ASSERT(NR > 0 && NC > 0); + + typedef op_uniform_matrix_2<T,NR,NC> op; + return matrix_op<op>(op(val)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename T, + long NR_, + long NC_, + T val + > + struct op_uniform_matrix : does_not_alias + { + const static long cost = 1; + const static long NR = NR_; + const static long NC = NC_; + typedef default_memory_manager mem_manager_type; + typedef row_major_layout layout_type; + typedef T type; + typedef const T const_ret_type; + const_ret_type apply ( long , long ) const { return val; } + + long nr() const { return NR; } + long nc() const { return NC; } + }; + + template < + typename T, + long NR, + long NC, + T val + > + const matrix_op<op_uniform_matrix<T,NR,NC,val> > uniform_matrix ( + ) + { + COMPILE_TIME_ASSERT(NR > 0 && NC > 0); + typedef op_uniform_matrix<T,NR,NC,val> op; + return matrix_op<op>(op()); + } + +// ---------------------------------------------------------------------------------------- + + struct op_gaussian_randm : does_not_alias + { + op_gaussian_randm ( + long nr_, + long nc_, + unsigned long seed_ + ) :_nr(nr_), _nc(nc_), seed(seed_){} + + const long _nr; + const long _nc; + const unsigned long seed; + + const static long cost = 100; + const static long NR = 0; + const static long NC = 0; + typedef default_memory_manager mem_manager_type; + typedef row_major_layout layout_type; + typedef double type; + typedef double const_ret_type; + const_ret_type apply ( long r, long c) const { return gaussian_random_hash(r,c,seed); } + + long nr() const { return _nr; } + long nc() const { return _nc; } + }; + + inline const matrix_op<op_gaussian_randm> gaussian_randm ( + long nr, + long nc, + unsigned long seed = 0 + ) + { + DLIB_ASSERT(nr >= 0 && nc >= 0, + "\tmatrix_exp gaussian_randm(nr, nc, seed)" + << "\n\tInvalid inputs to this function" + << "\n\tnr: " << nr + << "\n\tnc: " << nc + ); + + typedef op_gaussian_randm op; + return matrix_op<op>(op(nr,nc,seed)); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_add_diag + { + op_add_diag( const M& m_, const typename M::type& value_) : m(m_), value(value_){} + const M& m; + const typename M::type value; + + const static long cost = M::cost+1; + const static long NR = M::NR; + const static long NC = M::NC; + typedef typename M::type type; + typedef const typename M::type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply ( long r, long c) const + { + if (r==c) + return m(r,c)+value; + else + return m(r,c); + } + + long nr () const { return m.nr(); } + long nc () const { return m.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.destructively_aliases(item); } + }; + +// ---------------------------------------------------------------------------------------- + + template < + typename T + > + struct op_identity_matrix_2 : does_not_alias + { + op_identity_matrix_2(const long& size_) : size(size_) {} + + const long size; + + const static long cost = 1; + const static long NR = 0; + const static long NC = 0; + typedef default_memory_manager mem_manager_type; + typedef row_major_layout layout_type; + typedef T type; + typedef const T const_ret_type; + const_ret_type apply (long r, long c) const { return static_cast<type>(r == c); } + + long nr() const { return size; } + long nc() const { return size; } + }; + + template < + typename T, + typename U + > + const matrix_diag_op<op_identity_matrix_2<T> > identity_matrix ( + const U& size + ) + { + // the size argument must be some scalar value, not a matrix! + COMPILE_TIME_ASSERT(is_matrix<U>::value == false); + + DLIB_ASSERT(size > 0, + "\tconst matrix_exp identity_matrix<T>(size)" + << "\n\tsize must be bigger than 0" + << "\n\tsize: " << size + ); + typedef op_identity_matrix_2<T> op; + return matrix_diag_op<op>(op(size)); + } + + template < + typename EXP + > + const matrix_diag_op<op_identity_matrix_2<typename EXP::type> > identity_matrix ( + const matrix_exp<EXP>& mat + ) + { + DLIB_ASSERT(mat.nr() == mat.nc(), + "\tconst matrix_exp identity_matrix(mat)" + << "\n\t mat must be a square matrix." + << "\n\t mat.nr(): " << mat.nr() + << "\n\t mat.nc(): " << mat.nc() + ); + typedef typename EXP::type T; + typedef op_identity_matrix_2<T> op; + return matrix_diag_op<op>(op(mat.nr())); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP, + typename T + > + const matrix_op<op_add_diag<EXP> > operator+ ( + const matrix_exp<EXP>& lhs, + const matrix_exp<matrix_diag_op<op_identity_matrix_2<T> > >& DLIB_IF_ASSERT(rhs) + ) + { + // both matrices must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<T,typename EXP::type>::value == true)); + + // You can only add matrices together if they both have the same number of rows and columns. + 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 + ); + + + typedef op_add_diag<EXP> op; + return matrix_op<op>(op(lhs.ref(),1)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP, + typename T + > + const matrix_op<op_add_diag<EXP> > operator+ ( + const matrix_exp<matrix_diag_op<op_identity_matrix_2<T> > >& DLIB_IF_ASSERT(lhs), + const matrix_exp<EXP>& rhs + ) + { + // both matrices must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<T,typename EXP::type>::value == true)); + + // You can only add matrices together if they both have the same number of rows and columns. + 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 + ); + + + typedef op_add_diag<EXP> op; + return matrix_op<op>(op(rhs.ref(),1)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename T, + long N + > + struct op_const_diag_matrix : does_not_alias + { + op_const_diag_matrix(const long& size_, const T& value_) : size(size_),value(value_) {} + + const long size; + const T value; + + const static long cost = 1; + const static long NR = N; + const static long NC = N; + typedef default_memory_manager mem_manager_type; + typedef row_major_layout layout_type; + typedef T type; + typedef const T const_ret_type; + const_ret_type apply (long r, long c) const + { + if (r == c) + return value; + else + return 0; + } + + long nr() const { return size; } + long nc() const { return size; } + }; + + template < + typename T, + typename U + > + const typename disable_if<is_matrix<U>, matrix_diag_op<op_const_diag_matrix<T,0> > >::type operator* ( + const matrix_exp<matrix_diag_op<op_identity_matrix_2<T> > >& m, + const U& value + ) + { + typedef op_const_diag_matrix<T,0> op; + return matrix_diag_op<op>(op(m.nr(), value)); + } + + template < + typename T, + typename U + > + const typename disable_if<is_matrix<U>, matrix_diag_op<op_const_diag_matrix<T,0> > >::type operator* ( + const U& value, + const matrix_exp<matrix_diag_op<op_identity_matrix_2<T> > >& m + ) + { + typedef op_const_diag_matrix<T,0> op; + return matrix_diag_op<op>(op(m.nr(), value)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP, + typename T, + long N + > + const matrix_op<op_add_diag<EXP> > operator+ ( + const matrix_exp<EXP>& lhs, + const matrix_exp<matrix_diag_op<op_const_diag_matrix<T,N> > >& rhs + ) + { + // both matrices must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<T,typename EXP::type>::value == true)); + + // You can only add matrices together if they both have the same number of rows and columns. + 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 + ); + + + typedef op_add_diag<EXP> op; + return matrix_op<op>(op(lhs.ref(),rhs.ref().op.value)); + } + + template < + typename EXP, + typename T, + long N + > + const matrix_op<op_add_diag<EXP> > operator+ ( + const matrix_exp<matrix_diag_op<op_const_diag_matrix<T,N> > >& lhs, + const matrix_exp<EXP>& rhs + ) + { + // both matrices must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<T,typename EXP::type>::value == true)); + + // You can only add matrices together if they both have the same number of rows and columns. + 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 + ); + + + typedef op_add_diag<EXP> op; + return matrix_op<op>(op(rhs.ref(),lhs.ref().op.value)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename T, + long N + > + struct op_identity_matrix : does_not_alias + { + const static long cost = 1; + const static long NR = N; + const static long NC = N; + typedef default_memory_manager mem_manager_type; + typedef row_major_layout layout_type; + typedef T type; + typedef const T const_ret_type; + const_ret_type apply ( long r, long c) const { return static_cast<type>(r == c); } + + long nr () const { return NR; } + long nc () const { return NC; } + }; + + template < + typename T, + long N + > + const matrix_diag_op<op_identity_matrix<T,N> > identity_matrix ( + ) + { + COMPILE_TIME_ASSERT(N > 0); + + typedef op_identity_matrix<T,N> op; + return matrix_diag_op<op>(op()); + } + + template < + typename T, + typename U, + long N + > + const typename disable_if<is_matrix<U>, matrix_diag_op<op_const_diag_matrix<T,N> > >::type operator* ( + const matrix_exp<matrix_diag_op<op_identity_matrix<T,N> > >& m, + const U& value + ) + { + typedef op_const_diag_matrix<T,N> op; + return matrix_diag_op<op>(op(m.nr(), value)); + } + + template < + typename T, + typename U, + long N + > + const typename disable_if<is_matrix<U>, matrix_diag_op<op_const_diag_matrix<T,N> > >::type operator* ( + const U& value, + const matrix_exp<matrix_diag_op<op_identity_matrix<T,N> > >& m + ) + { + typedef op_const_diag_matrix<T,N> op; + return matrix_diag_op<op>(op(m.nr(), value)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP, + typename T, + long N + > + const matrix_op<op_add_diag<EXP> > operator+ ( + const matrix_exp<matrix_diag_op<op_identity_matrix<T,N> > >& DLIB_IF_ASSERT(lhs), + const matrix_exp<EXP>& rhs + ) + { + // both matrices must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<T,typename EXP::type>::value == true)); + + // You can only add matrices together if they both have the same number of rows and columns. + 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 + ); + + + typedef op_add_diag<EXP> op; + return matrix_op<op>(op(rhs.ref(),1)); + } + + template < + typename EXP, + typename T, + long N + > + const matrix_op<op_add_diag<EXP> > operator+ ( + const matrix_exp<EXP>& lhs, + const matrix_exp<matrix_diag_op<op_identity_matrix<T,N> > >& DLIB_IF_ASSERT(rhs) + ) + { + // both matrices must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<T,typename EXP::type>::value == true)); + + // You can only add matrices together if they both have the same number of rows and columns. + 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 + ); + + + typedef op_add_diag<EXP> op; + return matrix_op<op>(op(lhs.ref(),1)); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M, long R, long C> + struct op_rotate + { + op_rotate(const M& m_) : m(m_) {} + const M& m; + + const static long cost = M::cost + 2; + const static long NR = M::NR; + const static long NC = M::NC; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + const_ret_type apply ( long r, long c) const { return m((r+R)%m.nr(),(c+C)%m.nc()); } + + long nr () const { return m.nr(); } + long nc () const { return m.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template < + long R, + long C, + typename EXP + > + const matrix_op<op_rotate<EXP,R,C> > rotate ( + const matrix_exp<EXP>& m + ) + { + typedef op_rotate<EXP,R,C> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + + namespace impl + { + // A template to tell me if two types can be multiplied together in a sensible way. Here + // I'm saying it is ok if they are both the same type or one is the complex version of the other. + template <typename T, typename U> struct compatible { static const bool value = false; typedef T type; }; + template <typename T> struct compatible<T,T> { static const bool value = true; typedef T type; }; + template <typename T> struct compatible<std::complex<T>,T> { static const bool value = true; typedef std::complex<T> type; }; + template <typename T> struct compatible<T,std::complex<T> > { static const bool value = true; typedef std::complex<T> type; }; + } + + + template <typename M1, typename M2> + struct op_pointwise_multiply : basic_op_mm<M1,M2> + { + op_pointwise_multiply( const M1& m1_, const M2& m2_) : basic_op_mm<M1,M2>(m1_,m2_){} + + typedef typename impl::compatible<typename M1::type, typename M2::type>::type type; + typedef const type const_ret_type; + const static long cost = M1::cost + M2::cost + 1; + + const_ret_type apply ( long r, long c) const + { return this->m1(r,c)*this->m2(r,c); } + }; + + template < + typename EXP1, + typename EXP2 + > + inline const matrix_op<op_pointwise_multiply<EXP1,EXP2> > pointwise_multiply ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b + ) + { + COMPILE_TIME_ASSERT((impl::compatible<typename EXP1::type,typename EXP2::type>::value == true)); + COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0); + COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NC == 0 || EXP2::NC == 0); + DLIB_ASSERT(a.nr() == b.nr() && + a.nc() == b.nc(), + "\tconst matrix_exp pointwise_multiply(const matrix_exp& a, const matrix_exp& b)" + << "\n\tYou can only make a do a pointwise multiply with two equally sized matrices" + << "\n\ta.nr(): " << a.nr() + << "\n\ta.nc(): " << a.nc() + << "\n\tb.nr(): " << b.nr() + << "\n\tb.nc(): " << b.nc() + ); + typedef op_pointwise_multiply<EXP1,EXP2> op; + return matrix_op<op>(op(a.ref(),b.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2, typename M3> + struct op_pointwise_multiply3 : basic_op_mmm<M1,M2,M3> + { + op_pointwise_multiply3( const M1& m1_, const M2& m2_, const M3& m3_) : + basic_op_mmm<M1,M2,M3>(m1_,m2_,m3_){} + + typedef typename M1::type type; + typedef const typename M1::type const_ret_type; + const static long cost = M1::cost + M2::cost + M3::cost + 2; + + const_ret_type apply (long r, long c) const + { return this->m1(r,c)*this->m2(r,c)*this->m3(r,c); } + }; + + template < + typename EXP1, + typename EXP2, + typename EXP3 + > + inline const matrix_op<op_pointwise_multiply3<EXP1,EXP2,EXP3> > + pointwise_multiply ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b, + const matrix_exp<EXP3>& c + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + COMPILE_TIME_ASSERT((is_same_type<typename EXP2::type,typename EXP3::type>::value == true)); + COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0); + COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NR == 0 || EXP2::NC == 0); + COMPILE_TIME_ASSERT(EXP2::NR == EXP3::NR || EXP2::NR == 0 || EXP3::NR == 0); + COMPILE_TIME_ASSERT(EXP2::NC == EXP3::NC || EXP2::NC == 0 || EXP3::NC == 0); + DLIB_ASSERT(a.nr() == b.nr() && + a.nc() == b.nc() && + b.nr() == c.nr() && + b.nc() == c.nc(), + "\tconst matrix_exp pointwise_multiply(a,b,c)" + << "\n\tYou can only make a do a pointwise multiply between equally sized matrices" + << "\n\ta.nr(): " << a.nr() + << "\n\ta.nc(): " << a.nc() + << "\n\tb.nr(): " << b.nr() + << "\n\tb.nc(): " << b.nc() + << "\n\tc.nr(): " << c.nr() + << "\n\tc.nc(): " << c.nc() + ); + + typedef op_pointwise_multiply3<EXP1,EXP2,EXP3> op; + return matrix_op<op>(op(a.ref(),b.ref(),c.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2, typename M3, typename M4> + struct op_pointwise_multiply4 : basic_op_mmmm<M1,M2,M3,M4> + { + op_pointwise_multiply4( const M1& m1_, const M2& m2_, const M3& m3_, const M4& m4_) : + basic_op_mmmm<M1,M2,M3,M4>(m1_,m2_,m3_,m4_){} + + typedef typename M1::type type; + typedef const typename M1::type const_ret_type; + const static long cost = M1::cost + M2::cost + M3::cost + M4::cost + 3; + + const_ret_type apply (long r, long c) const + { return this->m1(r,c)*this->m2(r,c)*this->m3(r,c)*this->m4(r,c); } + }; + + + template < + typename EXP1, + typename EXP2, + typename EXP3, + typename EXP4 + > + inline const matrix_op<op_pointwise_multiply4<EXP1,EXP2,EXP3,EXP4> > pointwise_multiply ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b, + const matrix_exp<EXP3>& c, + const matrix_exp<EXP4>& d + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + COMPILE_TIME_ASSERT((is_same_type<typename EXP2::type,typename EXP3::type>::value == true)); + COMPILE_TIME_ASSERT((is_same_type<typename EXP3::type,typename EXP4::type>::value == true)); + COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0); + COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NC == 0 || EXP2::NC == 0 ); + COMPILE_TIME_ASSERT(EXP2::NR == EXP3::NR || EXP2::NR == 0 || EXP3::NR == 0); + COMPILE_TIME_ASSERT(EXP2::NC == EXP3::NC || EXP2::NC == 0 || EXP3::NC == 0); + COMPILE_TIME_ASSERT(EXP3::NR == EXP4::NR || EXP3::NR == 0 || EXP4::NR == 0); + COMPILE_TIME_ASSERT(EXP3::NC == EXP4::NC || EXP3::NC == 0 || EXP4::NC == 0); + DLIB_ASSERT(a.nr() == b.nr() && + a.nc() == b.nc() && + b.nr() == c.nr() && + b.nc() == c.nc() && + c.nr() == d.nr() && + c.nc() == d.nc(), + "\tconst matrix_exp pointwise_multiply(a,b,c,d)" + << "\n\tYou can only make a do a pointwise multiply between equally sized matrices" + << "\n\ta.nr(): " << a.nr() + << "\n\ta.nc(): " << a.nc() + << "\n\tb.nr(): " << b.nr() + << "\n\tb.nc(): " << b.nc() + << "\n\tc.nr(): " << c.nr() + << "\n\tc.nc(): " << c.nc() + << "\n\td.nr(): " << d.nr() + << "\n\td.nc(): " << d.nc() + ); + + typedef op_pointwise_multiply4<EXP1,EXP2,EXP3,EXP4> op; + return matrix_op<op>(op(a.ref(),b.ref(),c.ref(),d.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename P, + int type = static_switch< + pixel_traits<P>::grayscale, + pixel_traits<P>::rgb, + pixel_traits<P>::hsi, + pixel_traits<P>::rgb_alpha, + pixel_traits<P>::lab + >::value + > + struct pixel_to_vector_helper; + + template <typename P> + struct pixel_to_vector_helper<P,1> + { + template <typename M> + static void assign ( + M& m, + const P& pixel + ) + { + m(0) = static_cast<typename M::type>(pixel); + } + }; + + template <typename P> + struct pixel_to_vector_helper<P,2> + { + template <typename M> + static void assign ( + M& m, + const P& pixel + ) + { + m(0) = static_cast<typename M::type>(pixel.red); + m(1) = static_cast<typename M::type>(pixel.green); + m(2) = static_cast<typename M::type>(pixel.blue); + } + }; + + template <typename P> + struct pixel_to_vector_helper<P,3> + { + template <typename M> + static void assign ( + M& m, + const P& pixel + ) + { + m(0) = static_cast<typename M::type>(pixel.h); + m(1) = static_cast<typename M::type>(pixel.s); + m(2) = static_cast<typename M::type>(pixel.i); + } + }; + + template <typename P> + struct pixel_to_vector_helper<P,4> + { + template <typename M> + static void assign ( + M& m, + const P& pixel + ) + { + m(0) = static_cast<typename M::type>(pixel.red); + m(1) = static_cast<typename M::type>(pixel.green); + m(2) = static_cast<typename M::type>(pixel.blue); + m(3) = static_cast<typename M::type>(pixel.alpha); + } + }; + + template <typename P> + struct pixel_to_vector_helper<P,5> + { + template <typename M> + static void assign ( + M& m, + const P& pixel + ) + { + m(0) = static_cast<typename M::type>(pixel.l); + m(1) = static_cast<typename M::type>(pixel.a); + m(2) = static_cast<typename M::type>(pixel.b); + } + }; + + + template < + typename T, + typename P + > + inline const matrix<T,pixel_traits<P>::num,1> pixel_to_vector ( + const P& pixel + ) + { + COMPILE_TIME_ASSERT(pixel_traits<P>::num > 0); + matrix<T,pixel_traits<P>::num,1> m; + pixel_to_vector_helper<P>::assign(m,pixel); + return m; + } + +// ---------------------------------------------------------------------------------------- + + template < + typename P, + int type = static_switch< + pixel_traits<P>::grayscale, + pixel_traits<P>::rgb, + pixel_traits<P>::hsi, + pixel_traits<P>::rgb_alpha, + pixel_traits<P>::lab + >::value + > + struct vector_to_pixel_helper; + + template <typename P> + struct vector_to_pixel_helper<P,1> + { + template <typename M> + static void assign ( + P& pixel, + const M& m + ) + { + pixel = static_cast<unsigned char>(m(0)); + } + }; + + template <typename P> + struct vector_to_pixel_helper<P,2> + { + template <typename M> + static void assign ( + P& pixel, + const M& m + ) + { + pixel.red = static_cast<unsigned char>(m(0)); + pixel.green = static_cast<unsigned char>(m(1)); + pixel.blue = static_cast<unsigned char>(m(2)); + } + }; + + template <typename P> + struct vector_to_pixel_helper<P,3> + { + template <typename M> + static void assign ( + P& pixel, + const M& m + ) + { + pixel.h = static_cast<unsigned char>(m(0)); + pixel.s = static_cast<unsigned char>(m(1)); + pixel.i = static_cast<unsigned char>(m(2)); + } + }; + + template <typename P> + struct vector_to_pixel_helper<P,4> + { + template <typename M> + static void assign ( + P& pixel, + const M& m + ) + { + pixel.red = static_cast<unsigned char>(m(0)); + pixel.green = static_cast<unsigned char>(m(1)); + pixel.blue = static_cast<unsigned char>(m(2)); + pixel.alpha = static_cast<unsigned char>(m(3)); + } + }; + + template <typename P> + struct vector_to_pixel_helper<P,5> + { + template <typename M> + static void assign ( + P& pixel, + const M& m + ) + { + pixel.l = static_cast<unsigned char>(m(0)); + pixel.a = static_cast<unsigned char>(m(1)); + pixel.b = static_cast<unsigned char>(m(2)); + } + }; + + template < + typename P, + typename EXP + > + inline void vector_to_pixel ( + P& pixel, + const matrix_exp<EXP>& vector + ) + { + COMPILE_TIME_ASSERT(pixel_traits<P>::num == matrix_exp<EXP>::NR); + COMPILE_TIME_ASSERT(matrix_exp<EXP>::NC == 1); + vector_to_pixel_helper<P>::assign(pixel,vector); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M, long lower, long upper> + struct op_clamp : basic_op_m<M> + { + op_clamp( const M& m_) : basic_op_m<M>(m_){} + + typedef typename M::type type; + typedef const typename M::type const_ret_type; + const static long cost = M::cost + 2; + + const_ret_type apply ( long r, long c) const + { + const type temp = this->m(r,c); + if (temp > static_cast<type>(upper)) + return static_cast<type>(upper); + else if (temp < static_cast<type>(lower)) + return static_cast<type>(lower); + else + return temp; + } + }; + + template < + long l, + long u, + typename EXP + > + const matrix_op<op_clamp<EXP,l,u> > clamp ( + const matrix_exp<EXP>& m + ) + { + typedef op_clamp<EXP,l,u> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_clamp2 : basic_op_m<M> + { + typedef typename M::type type; + + op_clamp2( const M& m_, const type& l, const type& u) : + basic_op_m<M>(m_), lower(l), upper(u){} + + const type& lower; + const type& upper; + + typedef const typename M::type const_ret_type; + const static long cost = M::cost + 2; + + const_ret_type apply ( long r, long c) const + { + const type temp = this->m(r,c); + if (temp > upper) + return upper; + else if (temp < lower) + return lower; + else + return temp; + } + }; + + template < + typename EXP + > + const matrix_op<op_clamp2<EXP> > clamp ( + const matrix_exp<EXP>& m, + const typename EXP::type& lower, + const typename EXP::type& upper + ) + { + typedef op_clamp2<EXP> op; + return matrix_op<op>(op(m.ref(),lower, upper)); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2, typename M3> + struct op_clamp_m : basic_op_mmm<M1,M2,M3> + { + op_clamp_m( const M1& m1_, const M2& m2_, const M3& m3_) : + basic_op_mmm<M1,M2,M3>(m1_,m2_,m3_){} + + typedef typename M1::type type; + typedef const typename M1::type const_ret_type; + const static long cost = M1::cost + M2::cost + M3::cost + 2; + + const_ret_type apply (long r, long c) const + { + const type val = this->m1(r,c); + const type lower = this->m2(r,c); + const type upper = this->m3(r,c); + if (val <= upper) + { + if (lower <= val) + return val; + else + return lower; + } + else + { + return upper; + } + } + }; + + template < + typename EXP1, + typename EXP2, + typename EXP3 + > + const matrix_op<op_clamp_m<EXP1,EXP2,EXP3> > + clamp ( + const matrix_exp<EXP1>& m, + const matrix_exp<EXP2>& lower, + const matrix_exp<EXP3>& upper + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + COMPILE_TIME_ASSERT((is_same_type<typename EXP2::type,typename EXP3::type>::value == true)); + COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || EXP1::NR == 0 || EXP2::NR == 0); + COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || EXP1::NR == 0 || EXP2::NC == 0); + COMPILE_TIME_ASSERT(EXP2::NR == EXP3::NR || EXP2::NR == 0 || EXP3::NR == 0); + COMPILE_TIME_ASSERT(EXP2::NC == EXP3::NC || EXP2::NC == 0 || EXP3::NC == 0); + DLIB_ASSERT(m.nr() == lower.nr() && + m.nc() == lower.nc() && + m.nr() == upper.nr() && + m.nc() == upper.nc(), + "\tconst matrix_exp clamp(m,lower,upper)" + << "\n\t Invalid inputs were given to this function." + << "\n\t m.nr(): " << m.nr() + << "\n\t m.nc(): " << m.nc() + << "\n\t lower.nr(): " << lower.nr() + << "\n\t lower.nc(): " << lower.nc() + << "\n\t upper.nr(): " << upper.nr() + << "\n\t upper.nc(): " << upper.nc() + ); + + typedef op_clamp_m<EXP1,EXP2,EXP3> op; + return matrix_op<op>(op(m.ref(),lower.ref(),upper.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_lowerbound : basic_op_m<M> + { + typedef typename M::type type; + + op_lowerbound( const M& m_, const type& thresh_) : + basic_op_m<M>(m_), thresh(thresh_){} + + const type& thresh; + + typedef const typename M::type const_ret_type; + const static long cost = M::cost + 2; + + const_ret_type apply ( long r, long c) const + { + const type temp = this->m(r,c); + if (temp >= thresh) + return temp; + else + return thresh; + } + }; + + template < + typename EXP + > + const matrix_op<op_lowerbound<EXP> > lowerbound ( + const matrix_exp<EXP>& m, + const typename EXP::type& thresh + ) + { + typedef op_lowerbound<EXP> op; + return matrix_op<op>(op(m.ref(), thresh)); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_upperbound : basic_op_m<M> + { + typedef typename M::type type; + + op_upperbound( const M& m_, const type& thresh_) : + basic_op_m<M>(m_), thresh(thresh_){} + + const type& thresh; + + typedef const typename M::type const_ret_type; + const static long cost = M::cost + 2; + + const_ret_type apply ( long r, long c) const + { + const type temp = this->m(r,c); + if (temp <= thresh) + return temp; + else + return thresh; + } + }; + + template < + typename EXP + > + const matrix_op<op_upperbound<EXP> > upperbound ( + const matrix_exp<EXP>& m, + const typename EXP::type& thresh + ) + { + typedef op_upperbound<EXP> op; + return matrix_op<op>(op(m.ref(), thresh)); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_reshape + { + op_reshape(const M& m_, const long& rows_, const long& cols_) : m(m_),rows(rows_),cols(cols_) {} + const M& m; + const long rows; + const long cols; + + const static long cost = M::cost+2; + const static long NR = 0; + const static long NC = 0; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + + const_ret_type apply ( long r, long c) const + { + const long idx = r*cols + c; + return m(idx/m.nc(), idx%m.nc()); + } + + long nr () const { return rows; } + long nc () const { return cols; } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template < + typename EXP + > + const matrix_op<op_reshape<EXP> > reshape ( + const matrix_exp<EXP>& m, + const long& rows, + const long& cols + ) + { + DLIB_ASSERT(m.size() == rows*cols && rows > 0 && cols > 0, + "\tconst matrix_exp reshape(m, rows, cols)" + << "\n\t The size of m must match the dimensions you want to reshape it into." + << "\n\t m.size(): " << m.size() + << "\n\t rows*cols: " << rows*cols + << "\n\t rows: " << rows + << "\n\t cols: " << cols + ); + + typedef op_reshape<EXP> op; + return matrix_op<op>(op(m.ref(), rows, cols)); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename EXP1, + typename EXP2 + > + typename disable_if<is_complex<typename EXP1::type>,bool>::type equal ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b, + const typename EXP1::type eps = 100*std::numeric_limits<typename EXP1::type>::epsilon() + ) + { + // check if the dimensions don't match + if (a.nr() != b.nr() || a.nc() != b.nc()) + return false; + + for (long r = 0; r < a.nr(); ++r) + { + for (long c = 0; c < a.nc(); ++c) + { + if (std::abs(a(r,c)-b(r,c)) > eps) + return false; + } + } + + // no non-equal points found so we return true + return true; + } + + template < + typename EXP1, + typename EXP2 + > + typename enable_if<is_complex<typename EXP1::type>,bool>::type equal ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b, + const typename EXP1::type::value_type eps = 100*std::numeric_limits<typename EXP1::type::value_type>::epsilon() + ) + { + // check if the dimensions don't match + if (a.nr() != b.nr() || a.nc() != b.nc()) + return false; + + for (long r = 0; r < a.nr(); ++r) + { + for (long c = 0; c < a.nc(); ++c) + { + if (std::abs(real(a(r,c)-b(r,c))) > eps || + std::abs(imag(a(r,c)-b(r,c))) > eps) + return false; + } + } + + // no non-equal points found so we return true + return true; + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2> + struct op_scale_columns + { + op_scale_columns(const M1& m1_, const M2& m2_) : m1(m1_), m2(m2_) {} + const M1& m1; + const M2& m2; + + const static long cost = M1::cost + M2::cost + 1; + typedef typename M1::type type; + typedef const typename M1::type const_ret_type; + typedef typename M1::mem_manager_type mem_manager_type; + typedef typename M1::layout_type layout_type; + const static long NR = M1::NR; + const static long NC = M1::NC; + + const_ret_type apply ( long r, long c) const { return m1(r,c)*m2(c); } + + long nr () const { return m1.nr(); } + long nc () const { return m1.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item) ; } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const + { return m1.destructively_aliases(item) || m2.aliases(item); } + }; + + template < + typename EXP1, + typename EXP2 + > + const matrix_op<op_scale_columns<EXP1,EXP2> > scale_columns ( + const matrix_exp<EXP1>& m, + const matrix_exp<EXP2>& v + ) + { + // Both arguments to this function must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + // The v argument must be a row or column vector. + COMPILE_TIME_ASSERT((EXP2::NC == 1 || EXP2::NC == 0) || (EXP2::NR == 1 || EXP2::NR == 0)); + + // figure out the compile time known length of v + const long v_len = ((EXP2::NR)*(EXP2::NC) == 0)? 0 : (tmax<EXP2::NR,EXP2::NC>::value); + + // the length of v must match the number of columns in m + COMPILE_TIME_ASSERT(EXP1::NC == v_len || EXP1::NC == 0 || v_len == 0); + + DLIB_ASSERT(is_vector(v) == true && v.size() == m.nc(), + "\tconst matrix_exp scale_columns(m, v)" + << "\n\tv must be a row or column vector and its length must match the number of columns in m" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\tv.nr(): " << v.nr() + << "\n\tv.nc(): " << v.nc() + ); + typedef op_scale_columns<EXP1,EXP2> op; + return matrix_op<op>(op(m.ref(),v.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2> + struct op_scale_columns_diag + { + op_scale_columns_diag(const M1& m1_, const M2& m2_) : m1(m1_), m2(m2_) {} + const M1& m1; + const M2& m2; + + const static long cost = M1::cost + M2::cost + 1; + typedef typename M1::type type; + typedef const typename M1::type const_ret_type; + typedef typename M1::mem_manager_type mem_manager_type; + typedef typename M1::layout_type layout_type; + const static long NR = M1::NR; + const static long NC = M1::NC; + + const_ret_type apply ( long r, long c) const { return m1(r,c)*m2(c,c); } + + long nr () const { return m1.nr(); } + long nc () const { return m1.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item) ; } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const + { return m1.destructively_aliases(item) || m2.aliases(item); } + }; + +// turn expressions of the form mat*diagonal_matrix into scale_columns(mat, d) + template < + typename EXP1, + typename EXP2 + > + const matrix_op<op_scale_columns_diag<EXP1,EXP2> > operator* ( + const matrix_exp<EXP1>& m, + const matrix_diag_exp<EXP2>& d + ) + { + // Both arguments to this function must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + + // figure out the compile time known length of d + const long v_len = ((EXP2::NR)*(EXP2::NC) == 0)? 0 : (tmax<EXP2::NR,EXP2::NC>::value); + + // the length of d must match the number of columns in m + COMPILE_TIME_ASSERT(EXP1::NC == v_len || EXP1::NC == 0 || v_len == 0); + + DLIB_ASSERT(m.nc() == d.nr(), + "\tconst matrix_exp operator*(m, d)" + << "\n\tmatrix dimensions don't match" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\td.nr(): " << d.nr() + << "\n\td.nc(): " << d.nc() + ); + typedef op_scale_columns_diag<EXP1,EXP2> op; + return matrix_op<op>(op(m.ref(),d.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2> + struct op_scale_rows + { + op_scale_rows(const M1& m1_, const M2& m2_) : m1(m1_), m2(m2_) {} + const M1& m1; + const M2& m2; + + const static long cost = M1::cost + M2::cost + 1; + typedef typename M1::type type; + typedef const typename M1::type const_ret_type; + typedef typename M1::mem_manager_type mem_manager_type; + typedef typename M1::layout_type layout_type; + const static long NR = M1::NR; + const static long NC = M1::NC; + + const_ret_type apply ( long r, long c) const { return m1(r,c)*m2(r); } + + long nr () const { return m1.nr(); } + long nc () const { return m1.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item) ; } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const + { return m1.destructively_aliases(item) || m2.aliases(item); } + }; + + template < + typename EXP1, + typename EXP2 + > + const matrix_op<op_scale_rows<EXP1,EXP2> > scale_rows ( + const matrix_exp<EXP1>& m, + const matrix_exp<EXP2>& v + ) + { + // Both arguments to this function must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + // The v argument must be a row or column vector. + COMPILE_TIME_ASSERT((EXP2::NC == 1 || EXP2::NC == 0) || (EXP2::NR == 1 || EXP2::NR == 0)); + + // figure out the compile time known length of v + const long v_len = ((EXP2::NR)*(EXP2::NC) == 0)? 0 : (tmax<EXP2::NR,EXP2::NC>::value); + + // the length of v must match the number of rows in m + COMPILE_TIME_ASSERT(EXP1::NR == v_len || EXP1::NR == 0 || v_len == 0); + + DLIB_ASSERT(is_vector(v) == true && v.size() == m.nr(), + "\tconst matrix_exp scale_rows(m, v)" + << "\n\tv must be a row or column vector and its length must match the number of rows in m" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\tv.nr(): " << v.nr() + << "\n\tv.nc(): " << v.nc() + ); + typedef op_scale_rows<EXP1,EXP2> op; + return matrix_op<op>(op(m.ref(),v.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2> + struct op_scale_rows_diag + { + op_scale_rows_diag(const M1& m1_, const M2& m2_) : m1(m1_), m2(m2_) {} + const M1& m1; + const M2& m2; + + const static long cost = M1::cost + M2::cost + 1; + typedef typename M1::type type; + typedef const typename M1::type const_ret_type; + typedef typename M1::mem_manager_type mem_manager_type; + typedef typename M1::layout_type layout_type; + const static long NR = M1::NR; + const static long NC = M1::NC; + + const_ret_type apply ( long r, long c) const { return m1(r,c)*m2(r,r); } + + long nr () const { return m1.nr(); } + long nc () const { return m1.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item) ; } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const + { return m1.destructively_aliases(item) || m2.aliases(item); } + }; + +// turn expressions of the form diagonal_matrix*mat into scale_rows(mat, d) + template < + typename EXP1, + typename EXP2 + > + const matrix_op<op_scale_rows_diag<EXP1,EXP2> > operator* ( + const matrix_diag_exp<EXP2>& d, + const matrix_exp<EXP1>& m + ) + { + // Both arguments to this function must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + + // figure out the compile time known length of d + const long v_len = ((EXP2::NR)*(EXP2::NC) == 0)? 0 : (tmax<EXP2::NR,EXP2::NC>::value); + + // the length of d must match the number of rows in m + COMPILE_TIME_ASSERT(EXP1::NR == v_len || EXP1::NR == 0 || v_len == 0); + + DLIB_ASSERT(d.nc() == m.nr(), + "\tconst matrix_exp operator*(d, m)" + << "\n\tThe dimensions of the d and m matrices don't match." + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\td.nr(): " << d.nr() + << "\n\td.nc(): " << d.nc() + ); + typedef op_scale_rows_diag<EXP1,EXP2> op; + return matrix_op<op>(op(m.ref(),d.ref())); + } + +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- + + /* + The idea here is to catch expressions of the form d*M*d where d is diagonal and M + is some square matrix and turn them into something equivalent to + pointwise_multiply(diag(d)*trans(diag(d)), M). + + The reason for this is that doing it this way is more numerically stable. In particular, + doing 2 matrix multiplies as suggested by d*M*d could result in an asymmetric matrix even + if M is symmetric to begin with. + */ + + template <typename M1, typename M2, typename M3> + struct op_diag_m_diag + { + // This operator represents M1*M2*M3 where M1 and M3 are diagonal + + op_diag_m_diag(const M1& m1_, const M2& m2_, const M3& m3_) : m1(m1_), m2(m2_), m3(m3_) {} + const M1& m1; + const M2& m2; + const M3& m3; + + const static long cost = M1::cost + M2::cost + M3::cost + 1; + typedef typename M2::type type; + typedef const typename M2::type const_ret_type; + typedef typename M2::mem_manager_type mem_manager_type; + typedef typename M2::layout_type layout_type; + const static long NR = M2::NR; + const static long NC = M2::NC; + + const_ret_type apply ( long r, long c) const { return (m1(r,r)*m3(c,c))*m2(r,c); } + + long nr () const { return m2.nr(); } + long nc () const { return m2.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item) || m3.aliases(item) ; } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const + { return m2.destructively_aliases(item) || m1.aliases(item) || m3.aliases(item) ; } + }; + + // catch d*(M*d) = EXP1*EXP2*EXP3 + template < + typename EXP1, + typename EXP2, + typename EXP3 + > + const matrix_op<op_diag_m_diag<EXP1,EXP2,EXP3> > operator* ( + const matrix_diag_exp<EXP1>& d, + const matrix_exp<matrix_op<op_scale_columns_diag<EXP2,EXP3> > >& m + ) + { + // Both arguments to this function must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + + // figure out the compile time known length of d + const long v_len = ((EXP1::NR)*(EXP1::NC) == 0)? 0 : (tmax<EXP1::NR,EXP1::NC>::value); + + // the length of d must match the number of rows in m + COMPILE_TIME_ASSERT(EXP2::NR == v_len || EXP2::NR == 0 || v_len == 0); + + DLIB_ASSERT(d.nc() == m.nr(), + "\tconst matrix_exp operator*(d, m)" + << "\n\tmatrix dimensions don't match" + << "\n\td.nr(): " << d.nr() + << "\n\td.nc(): " << d.nc() + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + ); + typedef op_diag_m_diag<EXP1,EXP2,EXP3> op; + return matrix_op<op>(op(d.ref(), m.ref().op.m1, m.ref().op.m2)); + } + + // catch (d*M)*d = EXP1*EXP2*EXP3 + template < + typename EXP1, + typename EXP2, + typename EXP3 + > + const matrix_op<op_diag_m_diag<EXP1,EXP2,EXP3> > operator* ( + const matrix_exp<matrix_op<op_scale_rows_diag<EXP2,EXP1> > >& m, + const matrix_diag_exp<EXP3>& d + ) + { + // Both arguments to this function must contain the same type of element + COMPILE_TIME_ASSERT((is_same_type<typename EXP3::type,typename EXP2::type>::value == true)); + + // figure out the compile time known length of d + const long v_len = ((EXP3::NR)*(EXP3::NC) == 0)? 0 : (tmax<EXP3::NR,EXP3::NC>::value); + + // the length of d must match the number of columns in m + COMPILE_TIME_ASSERT(EXP2::NC == v_len || EXP2::NC == 0 || v_len == 0); + + DLIB_ASSERT(m.nc() == d.nr(), + "\tconst matrix_exp operator*(m, d)" + << "\n\tmatrix dimensions don't match" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\td.nr(): " << d.nr() + << "\n\td.nc(): " << d.nc() + ); + typedef op_diag_m_diag<EXP1,EXP2,EXP3> op; + return matrix_op<op>(op(m.ref().op.m2, m.ref().op.m1, d.ref())); + } + +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- + + struct sort_columns_sort_helper + { + template <typename T> + bool operator() ( + const T& item1, + const T& item2 + ) const + { + return item1.first < item2.first; + } + }; + + template < + typename T, long NR, long NC, typename mm, typename l1, + long NR2, long NC2, typename mm2, typename l2 + > + void sort_columns ( + matrix<T,NR,NC,mm,l1>& m, + matrix<T,NR2,NC2,mm2,l2>& v + ) + { + COMPILE_TIME_ASSERT(NC2 == 1 || NC2 == 0); + COMPILE_TIME_ASSERT(NC == NR2 || NC == 0 || NR2 == 0); + + DLIB_ASSERT(is_col_vector(v) == true && v.size() == m.nc(), + "\tconst matrix_exp sort_columns(m, v)" + << "\n\tv must be a column vector and its length must match the number of columns in m" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\tv.nr(): " << v.nr() + << "\n\tv.nc(): " << v.nc() + ); + + + + // Now we have to sort the given vectors in the m matrix according + // to how big their corresponding v(column index) values are. + typedef std::pair<T, matrix<T,0,1,mm> > col_pair; + typedef std_allocator<col_pair, mm> alloc; + std::vector<col_pair,alloc> colvalues; + col_pair p; + for (long r = 0; r < v.nr(); ++r) + { + p.first = v(r); + p.second = colm(m,r); + colvalues.push_back(p); + } + std::sort(colvalues.begin(), colvalues.end(), sort_columns_sort_helper()); + + for (long i = 0; i < v.nr(); ++i) + { + v(i) = colvalues[i].first; + set_colm(m,i) = colvalues[i].second; + } + + } + +// ---------------------------------------------------------------------------------------- + + template < + typename T, long NR, long NC, typename mm, typename l1, + long NR2, long NC2, typename mm2, typename l2 + > + void rsort_columns ( + matrix<T,NR,NC,mm,l1>& m, + matrix<T,NR2,NC2,mm2,l2>& v + ) + { + COMPILE_TIME_ASSERT(NC2 == 1 || NC2 == 0); + COMPILE_TIME_ASSERT(NC == NR2 || NC == 0 || NR2 == 0); + + DLIB_ASSERT(is_col_vector(v) == true && v.size() == m.nc(), + "\tconst matrix_exp rsort_columns(m, v)" + << "\n\tv must be a column vector and its length must match the number of columns in m" + << "\n\tm.nr(): " << m.nr() + << "\n\tm.nc(): " << m.nc() + << "\n\tv.nr(): " << v.nr() + << "\n\tv.nc(): " << v.nc() + ); + + + + // Now we have to sort the given vectors in the m matrix according + // to how big their corresponding v(column index) values are. + typedef std::pair<T, matrix<T,0,1,mm> > col_pair; + typedef std_allocator<col_pair, mm> alloc; + std::vector<col_pair,alloc> colvalues; + col_pair p; + for (long r = 0; r < v.nr(); ++r) + { + p.first = v(r); + p.second = colm(m,r); + colvalues.push_back(p); + } + std::sort(colvalues.rbegin(), colvalues.rend(), sort_columns_sort_helper()); + + for (long i = 0; i < v.nr(); ++i) + { + v(i) = colvalues[i].first; + set_colm(m,i) = colvalues[i].second; + } + + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2> + struct op_tensor_product + { + op_tensor_product(const M1& m1_, const M2& m2_) : m1(m1_),m2(m2_) {} + const M1& m1; + const M2& m2; + + const static long cost = M1::cost + M2::cost + 1; + const static long NR = M1::NR*M2::NR; + const static long NC = M1::NC*M2::NC; + typedef typename M1::type type; + typedef const typename M1::type const_ret_type; + typedef typename M1::mem_manager_type mem_manager_type; + typedef typename M1::layout_type layout_type; + + const_ret_type apply ( long r, long c) const + { + return m1(r/m2.nr(),c/m2.nc())*m2(r%m2.nr(),c%m2.nc()); + } + + long nr () const { return m1.nr()*m2.nr(); } + long nc () const { return m1.nc()*m2.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item); } + }; + + template < + typename EXP1, + typename EXP2 + > + inline const matrix_op<op_tensor_product<EXP1,EXP2> > tensor_product ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + typedef op_tensor_product<EXP1,EXP2> op; + return matrix_op<op>(op(a.ref(),b.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_make_symmetric : basic_op_m<M> + { + op_make_symmetric ( const M& m_) : basic_op_m<M>(m_){} + + const static long cost = M::cost+1; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + const_ret_type apply ( long r, long c) const + { + if (r >= c) + return this->m(r,c); + else + return this->m(c,r); + } + }; + + template < + typename EXP + > + const matrix_op<op_make_symmetric<EXP> > make_symmetric ( + const matrix_exp<EXP>& m + ) + { + DLIB_ASSERT(m.nr() == m.nc(), + "\tconst matrix make_symmetric(m)" + << "\n\t m must be a square matrix" + << "\n\t m.nr(): " << m.nr() + << "\n\t m.nc(): " << m.nc() + ); + + typedef op_make_symmetric<EXP> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_lowerm : basic_op_m<M> + { + op_lowerm( const M& m_) : basic_op_m<M>(m_){} + + const static long cost = M::cost+2; + typedef typename M::type type; + typedef const typename M::type const_ret_type; + const_ret_type apply ( long r, long c) const + { + if (r >= c) + return this->m(r,c); + else + return 0; + } + }; + + template <typename M> + struct op_lowerm_s : basic_op_m<M> + { + typedef typename M::type type; + op_lowerm_s( const M& m_, const type& s_) : basic_op_m<M>(m_), s(s_){} + + const type s; + + const static long cost = M::cost+2; + typedef const typename M::type const_ret_type; + const_ret_type apply ( long r, long c) const + { + if (r > c) + return this->m(r,c); + else if (r==c) + return s; + else + return 0; + } + }; + + template < + typename EXP + > + const matrix_op<op_lowerm<EXP> > lowerm ( + const matrix_exp<EXP>& m + ) + { + typedef op_lowerm<EXP> op; + return matrix_op<op>(op(m.ref())); + } + + template < + typename EXP + > + const matrix_op<op_lowerm_s<EXP> > lowerm ( + const matrix_exp<EXP>& m, + typename EXP::type s + ) + { + typedef op_lowerm_s<EXP> op; + return matrix_op<op>(op(m.ref(),s)); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_upperm : basic_op_m<M> + { + op_upperm( const M& m_) : basic_op_m<M>(m_){} + + const static long cost = M::cost+2; + typedef typename M::type type; + typedef const typename M::type const_ret_type; + const_ret_type apply ( long r, long c) const + { + if (r <= c) + return this->m(r,c); + else + return 0; + } + }; + + template <typename M> + struct op_upperm_s : basic_op_m<M> + { + typedef typename M::type type; + op_upperm_s( const M& m_, const type& s_) : basic_op_m<M>(m_), s(s_){} + + const type s; + + const static long cost = M::cost+2; + typedef const typename M::type const_ret_type; + const_ret_type apply ( long r, long c) const + { + if (r < c) + return this->m(r,c); + else if (r==c) + return s; + else + return 0; + } + }; + + template < + typename EXP + > + const matrix_op<op_upperm<EXP> > upperm ( + const matrix_exp<EXP>& m + ) + { + typedef op_upperm<EXP> op; + return matrix_op<op>(op(m.ref())); + } + + template < + typename EXP + > + const matrix_op<op_upperm_s<EXP> > upperm ( + const matrix_exp<EXP>& m, + typename EXP::type s + ) + { + typedef op_upperm_s<EXP> op; + return matrix_op<op>(op(m.ref(),s)); + } + +// ---------------------------------------------------------------------------------------- + + template <typename rand_gen> + inline const matrix<double> randm( + long nr, + long nc, + rand_gen& rnd + ) + { + DLIB_ASSERT(nr >= 0 && nc >= 0, + "\tconst matrix randm(nr, nc, rnd)" + << "\n\tInvalid inputs to this function" + << "\n\tnr: " << nr + << "\n\tnc: " << nc + ); + + matrix<double> m(nr,nc); + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + m(r,c) = rnd.get_random_double(); + } + } + + return m; + } + +// ---------------------------------------------------------------------------------------- + + inline const matrix<double> randm( + long nr, + long nc + ) + { + DLIB_ASSERT(nr >= 0 && nc >= 0, + "\tconst matrix randm(nr, nc)" + << "\n\tInvalid inputs to this function" + << "\n\tnr: " << nr + << "\n\tnc: " << nc + ); + + matrix<double> m(nr,nc); + // make a double that contains RAND_MAX + the smallest number that still + // makes the resulting double slightly bigger than static_cast<double>(RAND_MAX) + double max_val = RAND_MAX; + max_val += std::numeric_limits<double>::epsilon()*RAND_MAX; + + for (long r = 0; r < m.nr(); ++r) + { + for (long c = 0; c < m.nc(); ++c) + { + m(r,c) = std::rand()/max_val; + } + } + + return m; + } + +// ---------------------------------------------------------------------------------------- + + inline const matrix_range_exp<double> linspace ( + double start, + double end, + long num + ) + { + DLIB_ASSERT(num >= 0, + "\tconst matrix_exp linspace(start, end, num)" + << "\n\tInvalid inputs to this function" + << "\n\tstart: " << start + << "\n\tend: " << end + << "\n\tnum: " << num + ); + + return matrix_range_exp<double>(start,end,num,false); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_linpiece + { + op_linpiece(const double val_, const M& joints_) : joints(joints_), val(val_){} + + const M& joints; + const double val; + + const static long cost = 10; + + const static long NR = (M::NR*M::NC==0) ? (0) : (M::NR*M::NC-1); + const static long NC = 1; + typedef typename M::type type; + typedef default_memory_manager mem_manager_type; + typedef row_major_layout layout_type; + + typedef type const_ret_type; + const_ret_type apply (long i, long ) const + { + if (joints(i) < val) + return std::min<type>(val,joints(i+1)) - joints(i); + else + return 0; + } + + long nr () const { return joints.size()-1; } + long nc () const { return 1; } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return joints.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return joints.aliases(item); } + }; + + template < typename EXP > + const matrix_op<op_linpiece<EXP> > linpiece ( + const double val, + const matrix_exp<EXP>& joints + ) + { + // make sure requires clause is not broken + DLIB_ASSERT(is_vector(joints) && joints.size() >= 2, + "\t matrix_exp linpiece()" + << "\n\t Invalid inputs were given to this function " + << "\n\t is_vector(joints): " << is_vector(joints) + << "\n\t joints.size(): " << joints.size() + ); +#ifdef ENABLE_ASSERTS + for (long i = 1; i < joints.size(); ++i) + { + DLIB_ASSERT(joints(i-1) < joints(i), + "\t matrix_exp linpiece()" + << "\n\t Invalid inputs were given to this function " + << "\n\t joints("<<i-1<<"): " << joints(i-1) + << "\n\t joints("<<i<<"): " << joints(i) + ); + } +#endif + + typedef op_linpiece<EXP> op; + return matrix_op<op>(op(val,joints.ref())); + } + +// ---------------------------------------------------------------------------------------- + + inline const matrix_log_range_exp<double> logspace ( + double start, + double end, + long num + ) + { + DLIB_ASSERT(num >= 0, + "\tconst matrix_exp logspace(start, end, num)" + << "\n\tInvalid inputs to this function" + << "\n\tstart: " << start + << "\n\tend: " << end + << "\n\tnum: " << num + ); + + return matrix_log_range_exp<double>(start,end,num); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2> + struct op_cart_prod + { + op_cart_prod(const M1& m1_, const M2& m2_) : m1(m1_),m2(m2_) {} + const M1& m1; + const M2& m2; + + const static long cost = M1::cost+M2::cost+1; + typedef typename M1::type type; + typedef const typename M1::const_ret_type const_ret_type; + + typedef typename M1::mem_manager_type mem_manager_type; + typedef typename M1::layout_type layout_type; + const static long NR = M1::NR+M2::NR; + const static long NC = M1::NC*M2::NC; + + const_ret_type apply ( long r, long c) const + { + if (r < m1.nr()) + return m1(r, c/m2.nc()); + else + return m2(r-m1.nr(), c%m2.nc()); + } + + long nr () const { return m1.nr() + m2.nr(); } + long nc () const { return m1.nc() * m2.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item); } + }; + + template < + typename EXP1, + typename EXP2 + > + const matrix_op<op_cart_prod<EXP1,EXP2> > cartesian_product ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + + typedef op_cart_prod<EXP1,EXP2> op; + return matrix_op<op>(op(a.ref(),b.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_mat_to_vect + { + op_mat_to_vect(const M& m_) : m(m_) {} + const M& m; + + const static long cost = M::cost+2; + const static long NR = M::NC*M::NR; + const static long NC = 1; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + + const_ret_type apply ( long r, long ) const { return m(r/m.nc(), r%m.nc()); } + + long nr () const { return m.size(); } + long nc () const { return 1; } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template < + typename EXP + > + const matrix_op<op_mat_to_vect<EXP> > reshape_to_column_vector ( + const matrix_exp<EXP>& m + ) + { + typedef op_mat_to_vect<EXP> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template < + typename T, + long NR_, + long NC_, + typename MM + > + struct op_mat_to_vect2 + { + typedef matrix<T,NR_,NC_,MM,row_major_layout> M; + op_mat_to_vect2(const M& m_) : m(m_) {} + const M& m; + + const static long cost = M::cost+2; + const static long NR = M::NC*M::NR; + const static long NC = 1; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + + const_ret_type apply ( long r, long ) const { return (&m(0,0))[r]; } + + long nr () const { return m.size(); } + long nc () const { return 1; } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + }; + + template < + typename T, + long NR, + long NC, + typename MM + > + const matrix_op<op_mat_to_vect2<T,NR,NC,MM> > reshape_to_column_vector ( + const matrix<T,NR,NC,MM,row_major_layout>& m + ) + { + typedef op_mat_to_vect2<T,NR,NC,MM> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2> + struct op_join_rows + { + op_join_rows(const M1& m1_, const M2& m2_) : m1(m1_),m2(m2_),_nr(std::max(m1.nr(),m2.nr())) {} + const M1& m1; + const M2& m2; + const long _nr; + + template <typename T, typename U, bool selection> + struct type_selector; + template <typename T, typename U> + struct type_selector<T,U,true> { typedef T type; }; + template <typename T, typename U> + struct type_selector<T,U,false> { typedef U type; }; + + // If both const_ret_types are references then we should use them as the const_ret_type type + // but otherwise we should use the normal type. + typedef typename M1::const_ret_type T1; + typedef typename M1::type T2; + typedef typename M2::const_ret_type T3; + typedef typename type_selector<T1, T2, is_reference_type<T1>::value && is_reference_type<T3>::value>::type const_ret_type; + + const static long cost = M1::cost + M2::cost + 1; + const static long NR = tmax<M1::NR, M2::NR>::value; + const static long NC = (M1::NC*M2::NC != 0)? (M1::NC+M2::NC) : (0); + typedef typename M1::type type; + typedef typename M1::mem_manager_type mem_manager_type; + typedef typename M1::layout_type layout_type; + + const_ret_type apply (long r, long c) const + { + if (c < m1.nc()) + return m1(r,c); + else + return m2(r,c-m1.nc()); + } + + long nr () const { return _nr; } + long nc () const { return m1.nc()+m2.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item); } + }; + + template < + typename EXP1, + typename EXP2 + > + inline const matrix_op<op_join_rows<EXP1,EXP2> > join_rows ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + // You are getting an error on this line because you are trying to join two matrices that + // don't have the same number of rows + COMPILE_TIME_ASSERT(EXP1::NR == EXP2::NR || (EXP1::NR*EXP2::NR == 0)); + + DLIB_ASSERT(a.nr() == b.nr() || a.size() == 0 || b.size() == 0, + "\tconst matrix_exp join_rows(const matrix_exp& a, const matrix_exp& b)" + << "\n\tYou can only use join_rows() if both matrices have the same number of rows" + << "\n\ta.nr(): " << a.nr() + << "\n\tb.nr(): " << b.nr() + << "\n\ta.nc(): " << a.nc() + << "\n\tb.nc(): " << b.nc() + ); + + typedef op_join_rows<EXP1,EXP2> op; + return matrix_op<op>(op(a.ref(),b.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M1, typename M2> + struct op_join_cols + { + op_join_cols(const M1& m1_, const M2& m2_) : m1(m1_),m2(m2_),_nc(std::max(m1.nc(),m2.nc())) {} + const M1& m1; + const M2& m2; + const long _nc; + + template <typename T, typename U, bool selection> + struct type_selector; + template <typename T, typename U> + struct type_selector<T,U,true> { typedef T type; }; + template <typename T, typename U> + struct type_selector<T,U,false> { typedef U type; }; + + // If both const_ret_types are references then we should use them as the const_ret_type type + // but otherwise we should use the normal type. + typedef typename M1::const_ret_type T1; + typedef typename M1::type T2; + typedef typename M2::const_ret_type T3; + typedef typename type_selector<T1, T2, is_reference_type<T1>::value && is_reference_type<T3>::value>::type const_ret_type; + + + + const static long cost = M1::cost + M2::cost + 1; + const static long NC = tmax<M1::NC, M2::NC>::value; + const static long NR = (M1::NR*M2::NR != 0)? (M1::NR+M2::NR) : (0); + typedef typename M1::type type; + typedef typename M1::mem_manager_type mem_manager_type; + typedef typename M1::layout_type layout_type; + + const_ret_type apply ( long r, long c) const + { + if (r < m1.nr()) + return m1(r,c); + else + return m2(r-m1.nr(),c); + } + + long nr () const { return m1.nr()+m2.nr(); } + long nc () const { return _nc; } + + + template <typename U> bool aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const + { return m1.aliases(item) || m2.aliases(item); } + }; + + template < + typename EXP1, + typename EXP2 + > + inline const matrix_op<op_join_cols<EXP1,EXP2> > join_cols ( + const matrix_exp<EXP1>& a, + const matrix_exp<EXP2>& b + ) + { + COMPILE_TIME_ASSERT((is_same_type<typename EXP1::type,typename EXP2::type>::value == true)); + // You are getting an error on this line because you are trying to join two matrices that + // don't have the same number of columns + COMPILE_TIME_ASSERT(EXP1::NC == EXP2::NC || (EXP1::NC*EXP2::NC == 0)); + + DLIB_ASSERT(a.nc() == b.nc() || a.size() == 0 || b.size() == 0, + "\tconst matrix_exp join_cols(const matrix_exp& a, const matrix_exp& b)" + << "\n\tYou can only use join_cols() if both matrices have the same number of columns" + << "\n\ta.nr(): " << a.nr() + << "\n\tb.nr(): " << b.nr() + << "\n\ta.nc(): " << a.nc() + << "\n\tb.nc(): " << b.nc() + ); + + typedef op_join_cols<EXP1,EXP2> op; + return matrix_op<op>(op(a.ref(),b.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_fliplr + { + op_fliplr( const M& m_) : m(m_){} + + const M& m; + + const static long cost = M::cost; + const static long NR = M::NR; + const static long NC = M::NC; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + + const_ret_type apply (long r, long c) const { return m(r,m.nc()-c-1); } + + long nr () const { return m.nr(); } + long nc () const { return m.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + + }; + + template < + typename M + > + const matrix_op<op_fliplr<M> > fliplr ( + const matrix_exp<M>& m + ) + { + typedef op_fliplr<M> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_flipud + { + op_flipud( const M& m_) : m(m_){} + + const M& m; + + const static long cost = M::cost; + const static long NR = M::NR; + const static long NC = M::NC; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + + const_ret_type apply (long r, long c) const { return m(m.nr()-r-1,c); } + + long nr () const { return m.nr(); } + long nc () const { return m.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + + }; + + template < + typename M + > + const matrix_op<op_flipud<M> > flipud ( + const matrix_exp<M>& m + ) + { + typedef op_flipud<M> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename M> + struct op_flip + { + op_flip( const M& m_) : m(m_){} + + const M& m; + + const static long cost = M::cost; + const static long NR = M::NR; + const static long NC = M::NC; + typedef typename M::type type; + typedef typename M::const_ret_type const_ret_type; + typedef typename M::mem_manager_type mem_manager_type; + typedef typename M::layout_type layout_type; + + const_ret_type apply (long r, long c) const { return m(m.nr()-r-1, m.nc()-c-1); } + + long nr () const { return m.nr(); } + long nc () const { return m.nc(); } + + template <typename U> bool aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); } + + }; + + template < + typename M + > + const matrix_op<op_flip<M> > flip ( + const matrix_exp<M>& m + ) + { + typedef op_flip<M> op; + return matrix_op<op>(op(m.ref())); + } + +// ---------------------------------------------------------------------------------------- + + template <typename T, long NR, long NC, typename MM, typename L> + uint32 hash ( + const matrix<T,NR,NC,MM,L>& item, + uint32 seed = 0 + ) + { + DLIB_ASSERT_HAS_STANDARD_LAYOUT(T); + + if (item.size() == 0) + return 0; + else + return murmur_hash3(&item(0,0), sizeof(T)*item.size(), seed); + } + +// ---------------------------------------------------------------------------------------- + +} + +#endif // DLIB_MATRIx_UTILITIES_ + |