// 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; };
// ----------------------------------------------------------------------------------------
/*!A remove_complex
This is a template that can be used to remove std::complex from the underlying type.
For example:
remove_complex<float>::type == float
remove_complex<std::complex<float> >::type == float
!*/
template <typename T>
struct remove_complex {typedef T type;};
template <typename T>
struct remove_complex<std::complex<T> > {typedef T type;};
template<typename T>
using remove_complex_t = typename remove_complex<T>::type;
// ----------------------------------------------------------------------------------------
/*!A add_complex
This is a template that can be used to add std::complex to the underlying type if it isn't already complex.
For example:
add_complex<float>::type == std::complex<float>
add_complex<std::complex<float> >::type == std::complex<float>
!*/
template <typename T>
struct add_complex {typedef std::complex<T> type;};
template <typename T>
struct add_complex<std::complex<T> > {typedef std::complex<T> type;};
template<typename T>
using add_complex_t = typename add_complex<T>::type;
// ----------------------------------------------------------------------------------------
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 M1, typename M2>
struct op_pointwise_divide : basic_op_mm<M1,M2>
{
op_pointwise_divide( 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_divide<EXP1,EXP2> > pointwise_divide (
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_divide(const matrix_exp& a, const matrix_exp& b)"
<< "\n\tYou can only make a do a pointwise divide 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_divide<EXP1,EXP2> op;
return matrix_op<op>(op(a.ref(),b.ref()));
}
// ----------------------------------------------------------------------------------------
template <typename M1, typename M2, typename M3>
struct op_pointwise_divide3 : basic_op_mmm<M1,M2,M3>
{
op_pointwise_divide3( 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_divide3<EXP1,EXP2,EXP3> >
pointwise_divide (
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_divide(a,b,c)"
<< "\n\tYou can only make a do a pointwise divide 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_divide3<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_divide4 : basic_op_mmmm<M1,M2,M3,M4>
{
op_pointwise_divide4( 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_divide4<EXP1,EXP2,EXP3,EXP4> > pointwise_divide (
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_divide(a,b,c,d)"
<< "\n\tYou can only make a do a pointwise divide 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_divide4<EXP1,EXP2,EXP3,EXP4> op;
return matrix_op<op>(op(a.ref(),b.ref(),c.ref(),d.ref()));
}
// ----------------------------------------------------------------------------------------
template <typename M1, typename M2>
struct op_pointwise_pow : basic_op_mm<M1, M2>
{
op_pointwise_pow(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 + 7;
const_ret_type apply(long r, long c) const
{ return std::pow(this->m1(r, c), this->m2(r, c)); }
};
template <
typename EXP1,
typename EXP2
>
inline const matrix_op<op_pointwise_pow<EXP1, EXP2>> pointwise_pow (
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_pow(const matrix_exp& a, const matrix_exp& b)"
<< "\n\tYou can only make a do a pointwise power 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_pow<EXP1, EXP2> op;
return matrix_op<op>(op(a.ref(), b.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<P>(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_