matrix.h

// Copyright (C) 2006  Davis E. King (davisking@users.sourceforge.net)
// License: Boost Software License   See LICENSE.txt for the full license.
#ifndef DLIB_MATRIx_
#define DLIB_MATRIx_

#include "matrix_fwd.h"
#include "matrix_abstract.h"
#include "../algs.h"
#include "../serialize.h"
#include "../enable_if.h"
#include <sstream>
#include <algorithm>
#include "../memory_manager.h"
#include "../is_kind.h"
#include "matrix_data_layout.h"

#ifdef _MSC_VER
// Disable the following warnings for Visual Studio

// This warning is:
//    "warning C4355: 'this' : used in base member initializer list"
// Which we get from this code but it is not an error so I'm turning this
// warning off and then turning it back on at the end of the file.
#pragma warning(disable : 4355)

#endif

namespace dlib
{

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

    template <
        typename T,
        long num_rows,
        long num_cols,
        typename mem_manager,
        typename layout
        >
    class matrix_ref
    {
    public:
        typedef T type;
        typedef matrix_ref ref_type;
        typedef mem_manager mem_manager_type;
        typedef layout layout_type;
        const static long NR = num_rows;
        const static long NC = num_cols;

        matrix_ref (
            const matrix<T,num_rows,num_cols,mem_manager,layout>& m_
        ) : m(m_) {}

        matrix_ref (
            const matrix_ref& i_
        ) : m(i_.m) {}

        const T& operator() (
            long r,
            long c
        ) const { return m(r,c); }

        long nr (
        ) const { return m.nr(); }

        long nc (
        ) const { return m.nc(); }

        long size (
        ) const { return m.size(); }

        template <typename U, long iNR, long iNC, typename mm, typename l >
        bool aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const  { return false; }

        template <typename U, long iNR, long iNC, typename mm, typename l>
        bool destructively_aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return false; }

        bool aliases (
            const matrix<T,num_rows,num_cols,mem_manager,layout>& item
        ) const { return (&m == &item); }

        const matrix_ref ref(
        ) const { return *this; }

    private:
        // no assignment operator
        matrix_ref& operator=(const matrix_ref&);

        const matrix<T,num_rows,num_cols,mem_manager,layout>& m; // This is the item contained by this expression.
    };

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

    // We want to return the compile time constant if our NR and NC dimensions
    // aren't zero but if they are then we want to call ref_.nx() and return
    // the correct values. 
    template < typename ref_type, long NR >
    struct get_nr_helper
    {
        static inline long get(const ref_type&) { return NR; }
    };

    template < typename ref_type >
    struct get_nr_helper<ref_type,0>
    {
        static inline long get(const ref_type& m) { return m.nr(); }
    };

    template < typename ref_type, long NC >
    struct get_nc_helper
    {
        static inline long get(const ref_type&) { return NC; }
    };

    template < typename ref_type >
    struct get_nc_helper<ref_type,0>
    {
        static inline long get(const ref_type& m) { return m.nc(); }
    };


    // the matrix_exp for statically sized matrices 
    template <
        typename EXP
        >
    class matrix_exp
    {
    public:
        typedef typename EXP::type type;
        typedef typename EXP::ref_type ref_type;
        typedef typename EXP::mem_manager_type mem_manager_type;
        const static long NR = EXP::NR;
        const static long NC = EXP::NC;
        typedef matrix<type,NR,NC,mem_manager_type> matrix_type;

        matrix_exp (
            const EXP& exp
        ) : ref_(exp.ref()) {}

        inline const type operator() (
            long r,
            long c
        ) const 
        { 
            DLIB_ASSERT(r < nr() && c < nc() && r >= 0 && c >= c, 
                "\tconst type matrix_exp::operator(r,c)"
                << "\n\tYou must give a valid row and column"
                << "\n\tr:    " << r 
                << "\n\tc:    " << c
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc() 
                << "\n\tthis: " << this
                );
            return ref_(r,c); 
        }

        const type operator() (
            long i
        ) const 
        {
            COMPILE_TIME_ASSERT(NC == 1 || NC == 0 || NR == 1 || NR == 0);
            DLIB_ASSERT(nc() == 1 || nr() == 1, 
                "\tconst type matrix_exp::operator(i)"
                << "\n\tYou can only use this operator on column or row vectors"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            DLIB_ASSERT( ((nc() == 1 && i < nr()) || (nr() == 1 && i < nc())) && i >= 0, 
                "\tconst type matrix_exp::operator(i)"
                << "\n\tYou must give a valid row/column number"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            if (nc() == 1)
                return ref_(i,0);
            else
                return ref_(0,i);
        }

        long size (
        ) const { return nr()*nc(); }

        long nr (
        ) const { return get_nr_helper<ref_type,NR>::get(ref_); }

        long nc (
        ) const { return get_nc_helper<ref_type,NC>::get(ref_); }

        template <typename U, long iNR, long iNC, typename mm, typename l >
        bool aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return ref_.aliases(item); }

        template <typename U, long iNR, long iNC , typename mm, typename l>
        bool destructively_aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return ref_.destructively_aliases(item); }

        const ref_type& ref (
        ) const { return ref_; }

        inline operator const type (
        ) const 
        {
            COMPILE_TIME_ASSERT(NC == 1 || NC == 0);
            COMPILE_TIME_ASSERT(NR == 1 || NR == 0);
            DLIB_ASSERT(nr() == 1 && nc() == 1, 
                "\tmatrix_exp::operator const type&() const"
                << "\n\tYou can only use this operator on a 1x1 matrix"
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            return ref_(0,0);
        }


    private:


        const ref_type ref_;
    };

// ----------------------------------------------------------------------------------------

    template <typename T>
    struct is_matrix<matrix_exp<T> > { static const bool value = true; }; 
    template <typename T, long NR, long NC, typename mm, typename l>
    struct is_matrix<matrix_ref<T,NR,NC,mm,l> > { static const bool value = true; }; 
    template <typename T, long NR, long NC, typename mm, typename l>
    struct is_matrix<matrix<T,NR,NC,mm,l> > { static const bool value = true; }; 
    template <typename T>
    struct is_matrix<T&> { static const bool value = is_matrix<T>::value; }; 
    template <typename T>
    struct is_matrix<const T&> { static const bool value = is_matrix<T>::value; }; 
    template <typename T>
    struct is_matrix<const T> { static const bool value = is_matrix<T>::value; }; 
    /*
        is_matrix<T>::value == 1 if T is a matrix type else 0
    */

 // ----------------------------------------------------------------------------------------

    // This template will perform the needed loop for element multiplication using whichever
    // dimension is provided as a compile time constant (if one is at all).
    template <
        typename LHS,
        typename RHS,
        long lhs_nc = LHS::NC,
        long rhs_nr = RHS::NR
        >
    struct matrix_multiply_helper 
    {
        typedef typename LHS::type type;
        inline const static type  eval (
            const RHS& rhs,
            const LHS& lhs,
            long r, 
            long c
        )  
        { 
            type temp = type();
            for (long i = 0; i < rhs.nr(); ++i)
            {
                temp += lhs(r,i)*rhs(i,c);
            }
            return temp;
        }
    };

    template <
        typename LHS,
        typename RHS,
        long lhs_nc 
        >
    struct matrix_multiply_helper <LHS,RHS,lhs_nc,0>
    {
        typedef typename LHS::type type;
        inline const static type  eval (
            const RHS& rhs,
            const LHS& lhs,
            long r, 
            long c
        )  
        { 
            type temp = type();
            for (long i = 0; i < lhs.nc(); ++i)
            {
                temp += lhs(r,i)*rhs(i,c);
            }
            return temp;
        }
    };

    template <
        typename LHS,
        typename RHS,
        unsigned long count = 0
        >
    class matrix_multiply_exp 
    {
        /*!
            REQUIREMENTS ON LHS AND RHS
                - they must be matrix_exp or matrix_ref objects (or
                  objects with a compatible interface).
        !*/
    public:
        typedef typename LHS::type type;
        typedef matrix_multiply_exp ref_type;
        typedef typename LHS::mem_manager_type mem_manager_type;
        const static long NR = LHS::NR;
        const static long NC = RHS::NC;

        matrix_multiply_exp (
            const matrix_multiply_exp& item
        ) : lhs(item.lhs), rhs(item.rhs) {}

        inline matrix_multiply_exp (
            const LHS& lhs_,
            const RHS& rhs_
        ) :
            lhs(lhs_),
            rhs(rhs_)
        {
            // You are trying to multiply two incompatible matrices together.  The number of columns 
            // in the matrix on the left must match the number of rows in the matrix on the right.
            COMPILE_TIME_ASSERT(LHS::NC == RHS::NR || LHS::NC*RHS::NR == 0);
            DLIB_ASSERT(lhs.nc() == rhs.nr(), 
                "\tconst matrix_exp operator*(const matrix_exp& lhs, const matrix_exp& rhs)"
                << "\n\tYou are trying to multiply two incompatible matrices together"
                << "\n\tlhs.nr(): " << lhs.nr()
                << "\n\tlhs.nc(): " << lhs.nc()
                << "\n\trhs.nr(): " << rhs.nr()
                << "\n\trhs.nc(): " << rhs.nc()
                << "\n\t&lhs: " << &lhs 
                << "\n\t&rhs: " << &rhs 
                );

            // You can't multiply matrices together if they don't both contain the same type of elements.
            COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true));
        }

        inline const type operator() (
            long r, 
            long c
        ) const 
        { 
            return matrix_multiply_helper<LHS,RHS>::eval(rhs,lhs,r,c);
        }

        long nr (
        ) const { return lhs.nr(); }

        long nc (
        ) const { return rhs.nc(); }

        template <typename U, long iNR, long iNC, typename mm, typename l >
        bool aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return lhs.aliases(item) || rhs.aliases(item); }

        template <typename U, long iNR, long iNC , typename mm, typename l>
        bool destructively_aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return aliases(item); }

        const ref_type& ref(
        ) const { return *this; }

        const LHS lhs;
        const RHS rhs;
    };

    template <
        typename T,
        long NR,
        long NC,
        typename EXP1,
        typename EXP2,
        typename MM,
        typename L
        >
    inline const matrix_exp<matrix_multiply_exp<EXP1, matrix_multiply_exp<EXP2,typename matrix<T,NR,NC,MM,L>::ref_type >,0 > > operator* (
        const matrix_exp<matrix_multiply_exp<EXP1,EXP2,1> >& m1,
        const matrix<T,NR,NC,MM,L>& m2
    )
    {
        // We are going to reorder the order of evaluation of the terms here.  This way the
        // multiplication will go faster.
        typedef matrix_multiply_exp<EXP2,typename matrix<T,NR,NC,MM,L>::ref_type > exp_inner;
        typedef matrix_multiply_exp<EXP1, exp_inner,0 >  exp_outer;
        return matrix_exp<exp_outer>(exp_outer(m1.ref().lhs,exp_inner(m1.ref().rhs,m2)));
    }

    template <
        typename EXP1,
        typename EXP2
        >
    inline const matrix_exp<matrix_multiply_exp<EXP1, EXP2 > >  operator* (
        const matrix_exp<EXP1>& m1,
        const matrix_exp<EXP2>& m2
    )
    {
        typedef matrix_multiply_exp<EXP1, EXP2>  exp;
        return matrix_exp<exp>(exp(m1.ref(),m2.ref()));
    }

    template <
        typename T,
        long NR,
        long NC,
        typename EXP,
        typename MM,
        typename L
        >
    inline const matrix_exp<matrix_multiply_exp<typename matrix<T,NR,NC,MM,L>::ref_type, matrix_exp<EXP> > >  operator* (
        const matrix<T,NR,NC,MM,L>& m1,
        const matrix_exp<EXP>& m2
    )
    {
        typedef matrix_multiply_exp<typename matrix<T,NR,NC,MM,L>::ref_type, matrix_exp<EXP> >  exp;
        return matrix_exp<exp>(exp(m1,m2));
    }

    template <
        typename T,
        long NR,
        long NC,
        typename EXP,
        typename MM,
        typename L
        >
    inline const matrix_exp<matrix_multiply_exp< matrix_exp<EXP>, typename matrix<T,NR,NC,MM,L>::ref_type, 1> >  operator* (
        const matrix_exp<EXP>& m1,
        const matrix<T,NR,NC,MM,L>& m2
    )
    {
        typedef matrix_multiply_exp< matrix_exp<EXP>, typename matrix<T,NR,NC,MM,L>::ref_type, 1 >  exp;
        return matrix_exp<exp>(exp(m1,m2));
    }

    template <
        typename T,
        long NR1,
        long NC1,
        long NR2,
        long NC2,
        typename MM1,
        typename MM2,
        typename L1,
        typename L2
        >
    inline const matrix_exp<matrix_multiply_exp<typename matrix<T,NR1,NC1,MM1,L1>::ref_type,typename matrix<T,NR2,NC2,MM2,L2>::ref_type > >  operator* (
        const matrix<T,NR1,NC1,MM1,L1>& m1,
        const matrix<T,NR2,NC2,MM2,L2>& m2
    )
    {
        typedef matrix_multiply_exp<typename matrix<T,NR1,NC1,MM1,L1>::ref_type, typename matrix<T,NR2,NC2,MM2,L2>::ref_type >  exp;
        return matrix_exp<exp>(exp(m1,m2));
    }

// ----------------------------------------------------------------------------------------

    template <
        typename LHS,
        typename RHS
        >
    class matrix_add_expression 
    {
        /*!
            REQUIREMENTS ON LHS AND RHS
                - they must be matrix_exp or matrix_ref objects (or
                  objects with a compatible interface).
        !*/
    public:
        typedef typename LHS::type type;
        typedef typename LHS::mem_manager_type mem_manager_type;
        typedef matrix_add_expression ref_type;
        const static long NR = (RHS::NR > LHS::NR) ? RHS::NR : LHS::NR;
        const static long NC = (RHS::NC > LHS::NC) ? RHS::NC : LHS::NC;

        matrix_add_expression (
            const matrix_add_expression& item
        ) : lhs(item.lhs), rhs(item.rhs) {}

        matrix_add_expression (
            const LHS& lhs_,
            const RHS& rhs_
        ) :
            lhs(lhs_),
            rhs(rhs_)
        {
            // You can only add matrices together if they both have the same number of rows and columns.
            COMPILE_TIME_ASSERT(LHS::NR == RHS::NR || LHS::NR == 0 || RHS::NR == 0);
            COMPILE_TIME_ASSERT(LHS::NC == RHS::NC || LHS::NC == 0 || RHS::NC == 0);
            DLIB_ASSERT(lhs.nc() == rhs.nc() &&
                   lhs.nr() == rhs.nr(), 
                "\tconst matrix_exp operator+(const matrix_exp& lhs, const matrix_exp& rhs)"
                << "\n\tYou are trying to add two incompatible matrices together"
                << "\n\tlhs.nr(): " << lhs.nr()
                << "\n\tlhs.nc(): " << lhs.nc()
                << "\n\trhs.nr(): " << rhs.nr()
                << "\n\trhs.nc(): " << rhs.nc()
                << "\n\t&lhs: " << &lhs 
                << "\n\t&rhs: " << &rhs 
                );

            // You can only add matrices together if they both contain the same types of elements.
            COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true));
        }

        const type operator() (
            long r, 
            long c
        ) const { return lhs(r,c) + rhs(r,c); }

        template <typename U, long iNR, long iNC , typename mm, typename l>
        bool aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return lhs.aliases(item) || rhs.aliases(item); }

        template <typename U, long iNR, long iNC, typename mm, typename l >
        bool destructively_aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return lhs.destructively_aliases(item) || rhs.destructively_aliases(item); }

        const ref_type& ref(
        ) const { return *this; }

        long nr (
        ) const { return lhs.nr(); }

        long nc (
        ) const { return lhs.nc(); }

        const LHS lhs;
        const RHS rhs;
    };

    template <
        typename EXP1,
        typename EXP2
        >
    inline const matrix_exp<matrix_add_expression<EXP1, EXP2 > > operator+ (
        const matrix_exp<EXP1>& m1,
        const matrix_exp<EXP2>& m2
    )
    {
        typedef matrix_add_expression<EXP1, EXP2 >  exp;
        return matrix_exp<exp>(exp(m1.ref(),m2.ref()));
    }

// ----------------------------------------------------------------------------------------

    template <
        typename LHS,
        typename RHS
        >
    class matrix_subtract_exp 
    {
        /*!
            REQUIREMENTS ON LHS AND RHS
                - they must be matrix_exp or matrix_ref objects (or
                  objects with a compatible interface).
        !*/
    public:
        typedef typename LHS::type type;
        typedef typename LHS::mem_manager_type mem_manager_type;
        typedef matrix_subtract_exp ref_type;
        const static long NR = (RHS::NR > LHS::NR) ? RHS::NR : LHS::NR;
        const static long NC = (RHS::NC > LHS::NC) ? RHS::NC : LHS::NC;

        matrix_subtract_exp (
            const LHS& lhs_,
            const RHS& rhs_
        ) :
            lhs(lhs_),
            rhs(rhs_)
        {
            // You can only subtract one matrix from another if they both have the same number of rows and columns.
            COMPILE_TIME_ASSERT(LHS::NR == RHS::NR || LHS::NR == 0 || RHS::NR == 0);
            COMPILE_TIME_ASSERT(LHS::NC == RHS::NC || LHS::NC == 0 || RHS::NC == 0);
            DLIB_ASSERT(lhs.nc() == rhs.nc() &&
                   lhs.nr() == rhs.nr(), 
                "\tconst matrix_exp operator-(const matrix_exp& lhs, const matrix_exp& rhs)"
                << "\n\tYou are trying to add two incompatible matrices together"
                << "\n\tlhs.nr(): " << lhs.nr()
                << "\n\tlhs.nc(): " << lhs.nc()
                << "\n\trhs.nr(): " << rhs.nr()
                << "\n\trhs.nc(): " << rhs.nc()
                << "\n\t&lhs: " << &lhs 
                << "\n\t&rhs: " << &rhs 
                );

            // You can only subtract one matrix from another if they both contain elements of the same type.
            COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true));
        }

        const type operator() (
            long r, 
            long c
        ) const { return lhs(r,c) - rhs(r,c); }

        template <typename U, long iNR, long iNC, typename mm, typename l >
        bool aliases (
            const matrix<U,iNR,iNC, mm,l>& item
        ) const { return lhs.aliases(item) || rhs.aliases(item); }

        template <typename U, long iNR, long iNC , typename mm, typename l>
        bool destructively_aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return lhs.destructively_aliases(item) || rhs.destructively_aliases(item); }

        const ref_type& ref(
        ) const { return *this; }

        long nr (
        ) const { return lhs.nr(); }

        long nc (
        ) const { return lhs.nc(); }

        const LHS lhs;
        const RHS rhs;
    };

    template <
        typename EXP1,
        typename EXP2
        >
    inline const matrix_exp<matrix_subtract_exp<EXP1, EXP2 > > operator- (
        const matrix_exp<EXP1>& m1,
        const matrix_exp<EXP2>& m2
    )
    {
        typedef matrix_subtract_exp<EXP1, EXP2 >  exp;
        return matrix_exp<exp>(exp(m1.ref(),m2.ref()));
    }

// ----------------------------------------------------------------------------------------

    template <
        typename M,
        typename S
        >
    class matrix_divscal_exp  
    {
        /*!
            REQUIREMENTS ON M 
                - must be a matrix_exp or matrix_ref object (or
                  an object with a compatible interface).

            REQUIREMENTS ON S
                - must be some kind of scalar type
        !*/
    public:
        typedef typename M::type type;
        typedef typename M::mem_manager_type mem_manager_type;
        typedef matrix_divscal_exp ref_type;
        const static long NR = M::NR;
        const static long NC = M::NC;

        matrix_divscal_exp (
            const M& m_,
            const S& s_
        ) :
            m(m_),
            s(s_)
        {}

        const type operator() (
            long r, 
            long c
        ) const { return m(r,c)/s; }

        template <typename U, long iNR, long iNC, typename mm , typename l>
        bool aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return m.aliases(item); }

        template <typename U, long iNR, long iNC, typename mm, typename l >
        bool destructively_aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return m.destructively_aliases(item); }

        const ref_type& ref(
        ) const { return *this; }

        long nr (
        ) const { return m.nr(); }

        long nc (
        ) const { return m.nc(); }

        const M m;
        const S s;
    };

    template <
        typename EXP,
        typename S 
        >
    inline const matrix_exp<matrix_divscal_exp<matrix_exp<EXP>, S> > operator/ (
        const matrix_exp<EXP>& m,
        const S& s
    )
    {
        typedef matrix_divscal_exp<matrix_exp<EXP>,S >  exp;
        return matrix_exp<exp>(exp(m,s));
    }

// ----------------------------------------------------------------------------------------

    template <
        typename M,
        typename S
        >
    class matrix_mulscal_exp  
    {
        /*!
            REQUIREMENTS ON M 
                - must be a matrix_exp or matrix_ref object (or
                  an object with a compatible interface).

            REQUIREMENTS ON S
                - must be some kind of scalar type
        !*/
    public:
        typedef typename M::type type;
        typedef typename M::mem_manager_type mem_manager_type;
        typedef matrix_mulscal_exp ref_type;
        const static long NR = M::NR;
        const static long NC = M::NC;

        matrix_mulscal_exp (
            const M& m_,
            const S& s_
        ) :
            m(m_),
            s(s_)
        {}

        const type operator() (
            long r, 
            long c
        ) const { return m(r,c)*s; }

        template <typename U, long iNR, long iNC , typename mm, typename l>
        bool aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return m.aliases(item); }

        template <typename U, long iNR, long iNC, typename mm, typename l >
        bool destructively_aliases (
            const matrix<U,iNR,iNC,mm,l>& item
        ) const { return m.destructively_aliases(item); }

        const ref_type& ref(
        ) const { return *this; }

        long nr (
        ) const { return m.nr(); }

        long nc (
        ) const { return m.nc(); }

        const M m;
        const S s;
    };

    template <
        typename EXP,
        typename S 
        >
    inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, S> > operator* (
        const matrix_exp<EXP>& m,
        const S& s
    )
    {
        typedef matrix_mulscal_exp<matrix_exp<EXP>,S >  exp;
        return matrix_exp<exp>(exp(m,s));
    }

    template <
        typename EXP,
        typename S 
        >
    inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, S> > operator* (
        const S& s,
        const matrix_exp<EXP>& m
    )
    {
        typedef matrix_mulscal_exp<matrix_exp<EXP>,S >  exp;
        return matrix_exp<exp>(exp(m,s));
    }

    template <
        typename EXP 
        >
    inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, float> > operator/ (
        const matrix_exp<EXP>& m,
        const float& s
    )
    {
        typedef matrix_mulscal_exp<matrix_exp<EXP>,float >  exp;
        return matrix_exp<exp>(exp(m,1.0f/s));
    }

    template <
        typename EXP
        >
    inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, double> > operator/ (
        const matrix_exp<EXP>& m,
        const double& s
    )
    {
        typedef matrix_mulscal_exp<matrix_exp<EXP>,double >  exp;
        return matrix_exp<exp>(exp(m,1.0/s));
    }

    template <
        typename EXP
        >
    inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, long double> > operator/ (
        const matrix_exp<EXP>& m,
        const long double& s
    )
    {
        typedef matrix_mulscal_exp<matrix_exp<EXP>,long double >  exp;
        return matrix_exp<exp>(exp(m,1.0/s));
    }

    template <
        typename EXP
        >
    inline const matrix_exp<matrix_mulscal_exp<matrix_exp<EXP>, int> > operator- (
        const matrix_exp<EXP>& m
    )
    {
        typedef matrix_mulscal_exp<matrix_exp<EXP>,int >  exp;
        return matrix_exp<exp>(exp(m,-1));
    }

// ----------------------------------------------------------------------------------------

    template <
        typename EXP1,
        typename EXP2
        >
    bool operator== (
        const matrix_exp<EXP1>& m1,
        const matrix_exp<EXP2>& m2
    )
    {
        if (m1.nr() == m2.nr() && m1.nc() == m2.nc())
        {
            for (long r = 0; r < m1.nr(); ++r)
            {
                for (long c = 0; c < m1.nc(); ++c)
                {
                    if (m1(r,c) != m2(r,c))
                        return false;
                }
            }
            return true;
        }
        return false;
    }

    template <
        typename EXP1,
        typename EXP2
        >
    inline bool operator!= (
        const matrix_exp<EXP1>& m1,
        const matrix_exp<EXP2>& m2
    ) { return !(m1 == m2); }

// ----------------------------------------------------------------------------------------

    template <
        typename matrix_dest_type,
        typename src_exp 
        >
    void matrix_assign (
        matrix_dest_type& dest,
        const matrix_exp<src_exp>& src
    )
    /*!
        requires
            - src.destructively_aliases(dest) == false
        ensures
            - #dest == src
            - the part of dest outside the above sub matrix remains unchanged
    !*/
    {
        for (long r = 0; r < src.nr(); ++r)
        {
            for (long c = 0; c < src.nc(); ++c)
            {
                dest(r,c) = src(r,c);
            }
        }
    }

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

    template <
        typename T,
        long num_rows,
        long num_cols,
        typename mem_manager,
        typename layout
        >
    class matrix : public matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager,layout> > 
    {

        COMPILE_TIME_ASSERT(num_rows >= 0 && num_cols >= 0); 

    public:
        typedef T type;
        typedef matrix_ref<T,num_rows,num_cols,mem_manager,layout> ref_type;
        typedef mem_manager mem_manager_type;
        typedef layout layout_type;
        const static long NR = num_rows;
        const static long NC = num_cols;

        matrix () : matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager, layout> >(ref_type(*this)) 
        {
        }

        explicit matrix (
            long length 
        ) : matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager, layout> >(ref_type(*this)) 
        {
            // This object you are trying to call matrix(length) on is not a column or 
            // row vector.
            COMPILE_TIME_ASSERT(NR == 1 || NC == 1);
            DLIB_ASSERT( length >= 0, 
                "\tmatrix::matrix(length)"
                << "\n\tlength must be at least 0"
                << "\n\tlength: " << length 
                << "\n\tNR:     " << NR 
                << "\n\tNC:     " << NC 
                << "\n\tthis:   " << this
                );

            if (NR == 1)
            {
                DLIB_ASSERT(NC == 0 || NC == length,
                    "\tmatrix::matrix(length)"
                    << "\n\tSince this is a statically sized matrix length must equal NC"
                    << "\n\tlength: " << length 
                    << "\n\tNR:     " << NR 
                    << "\n\tNC:     " << NC 
                    << "\n\tthis:   " << this
                    );

                data.set_size(1,length);
            }
            else
            {
                DLIB_ASSERT(NR == 0 || NR == length,
                    "\tvoid matrix::set_size(length)"
                    << "\n\tSince this is a statically sized matrix length must equal NR"
                    << "\n\tlength: " << length 
                    << "\n\tNR:     " << NR 
                    << "\n\tNC:     " << NC 
                    << "\n\tthis:   " << this
                    );

                data.set_size(length,1);
            }
        }

        matrix (
            long rows,
            long cols 
        ) : matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager, layout> >(ref_type(*this)) 
        {
            DLIB_ASSERT( (NR == 0 || NR == rows) && ( NC == 0 || NC == cols) && 
                    rows >= 0 && cols >= 0, 
                "\tvoid matrix::matrix(rows, cols)"
                << "\n\tYou have supplied conflicting matrix dimensions"
                << "\n\trows: " << rows
                << "\n\tcols: " << cols
                << "\n\tNR:   " << NR 
                << "\n\tNC:   " << NC 
                );
            data.set_size(rows,cols);
        }

        template <typename EXP>
        matrix (
            const matrix_exp<EXP>& m
        ): matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager, layout> >(ref_type(*this)) 
        {
            // You get an error on this line if the matrix m contains a type that isn't
            // the same as the type contained in the target matrix.
            COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true) ||
                                (is_matrix<typename EXP::type>::value == true));

            // The matrix you are trying to assign m to is a statically sized matrix and 
            // m's dimensions don't match that of *this. 
            COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0);
            COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0);
            DLIB_ASSERT((NR == 0 || NR == m.nr()) && (NC == 0 || NC == m.nc()), 
                "\tmatrix& matrix::matrix(const matrix_exp& m)"
                << "\n\tYou are trying to assign a dynamically sized matrix to a statically sized matrix with the wrong size"
                << "\n\tNR:     " << NR
                << "\n\tNC:     " << NC
                << "\n\tm.nr(): " << m.nr()
                << "\n\tm.nc(): " << m.nc()
                << "\n\tthis:   " << this
                );

            data.set_size(m.nr(),m.nc());

            matrix_assign(*this, m);
        }

        matrix (
            const matrix& m
        ): matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager, layout> >(ref_type(*this)) 
        {
            data.set_size(m.nr(),m.nc());
            matrix_assign(*this, m);
        }

        template <typename U, size_t len>
        matrix (
            U (&array)[len]
        ): matrix_exp<matrix_ref<T,num_rows,num_cols, mem_manager, layout> >(ref_type(*this)) 
        {
            COMPILE_TIME_ASSERT(NR*NC == len && len > 0);
            size_t idx = 0;
            for (long r = 0; r < NR; ++r)
            {
                for (long c = 0; c < NC; ++c)
                {
                    data(r,c) = static_cast<T>(array[idx]);
                    ++idx;
                }
            }
        }

        T& operator() (
            long r, 
            long c
        ) 
        { 
            DLIB_ASSERT(r < nr() && c < nc() &&
                   r >= 0 && c >= 0, 
                "\tT& matrix::operator(r,c)"
                << "\n\tYou must give a valid row and column"
                << "\n\tr:    " << r 
                << "\n\tc:    " << c
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc() 
                << "\n\tthis: " << this
                );
            return data(r,c); 
        }

        const T& operator() (
            long r, 
            long c
        ) const 
        { 
            DLIB_ASSERT(r < nr() && c < nc() &&
                   r >= 0 && c >= 0, 
                "\tconst T& matrix::operator(r,c)"
                << "\n\tYou must give a valid row and column"
                << "\n\tr:    " << r 
                << "\n\tc:    " << c
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc() 
                << "\n\tthis: " << this
                );
            return data(r,c);
        }

        T& operator() (
            long i
        ) 
        {
            // You can only use this operator on column vectors.
            COMPILE_TIME_ASSERT(NC == 1 || NC == 0 || NR == 1 || NR == 0);
            DLIB_ASSERT(nc() == 1 || nr() == 1, 
                "\tconst type matrix::operator(i)"
                << "\n\tYou can only use this operator on column or row vectors"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            DLIB_ASSERT( ((nc() == 1 && i < nr()) || (nr() == 1 && i < nc())) && i >= 0, 
                "\tconst type matrix::operator(i)"
                << "\n\tYou must give a valid row/column number"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            return data(i);
        }

        const T& operator() (
            long i
        ) const
        {
            // You can only use this operator on column vectors.
            COMPILE_TIME_ASSERT(NC == 1 || NC == 0 || NR == 1 || NR == 0);
            DLIB_ASSERT(nc() == 1 || nr() == 1, 
                "\tconst type matrix::operator(i)"
                << "\n\tYou can only use this operator on column or row vectors"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            DLIB_ASSERT( ((nc() == 1 && i < nr()) || (nr() == 1 && i < nc())) && i >= 0, 
                "\tconst type matrix::operator(i)"
                << "\n\tYou must give a valid row/column number"
                << "\n\ti:    " << i
                << "\n\tnr(): " << nr()
                << "\n\tnc(): " << nc()
                << "\n\tthis: " << this
                );
            return data(i);
        }

        inline operator const type (
        ) const 
        {
            COMPILE_TIME_ASSERT(NC == 1 || NC == 0);
            COMPILE_TIME_ASSERT(NR == 1 || NR == 0);
            DLIB_ASSERT( nr() == 1 && nc() == 1 , 
                "\tmatrix::operator const type"
                << "\n\tYou can only attempt to implicit convert a matrix to a scalar if"
                << "\n\tthe matrix is a 1x1 matrix"
                << "\n\tnr(): " << nr() 
                << "\n\tnc(): " << nc() 
                << "\n\tthis: " << this
                );
            return data(0);
        }

        void set_size (
            long rows,
            long cols
        )
        {
            DLIB_ASSERT( (NR == 0 || NR == rows) && ( NC == 0 || NC == cols) &&
                    rows >= 0 && cols >= 0, 
                "\tvoid matrix::set_size(rows, cols)"
                << "\n\tYou have supplied conflicting matrix dimensions"
                << "\n\trows: " << rows
                << "\n\tcols: " << cols
                << "\n\tNR:   " << NR 
                << "\n\tNC:   " << NC 
                << "\n\tthis: " << this
                );
            if (nr() != rows || nc() != cols)
                data.set_size(rows,cols);
        }

        void set_size (
            long length
        )
        {
            // This object you are trying to call set_size(length) on is not a column or 
            // row vector.
            COMPILE_TIME_ASSERT(NR == 1 || NC == 1);
            DLIB_ASSERT( length >= 0, 
                "\tvoid matrix::set_size(length)"
                << "\n\tlength must be at least 0"
                << "\n\tlength: " << length 
                << "\n\tNR:     " << NR 
                << "\n\tNC:     " << NC 
                << "\n\tthis:   " << this
                );

            if (NR == 1)
            {
                DLIB_ASSERT(NC == 0 || NC == length,
                    "\tvoid matrix::set_size(length)"
                    << "\n\tSince this is a statically sized matrix length must equal NC"
                    << "\n\tlength: " << length 
                    << "\n\tNR:     " << NR 
                    << "\n\tNC:     " << NC 
                    << "\n\tthis:   " << this
                    );

                if (nc() != length)
                    data.set_size(1,length);
            }
            else
            {
                DLIB_ASSERT(NR == 0 || NR == length,
                    "\tvoid matrix::set_size(length)"
                    << "\n\tSince this is a statically sized matrix length must equal NR"
                    << "\n\tlength: " << length 
                    << "\n\tNR:     " << NR 
                    << "\n\tNC:     " << NC 
                    << "\n\tthis:   " << this
                    );

                if (nr() != length)
                    data.set_size(length,1);
            }
        }

        long nr (
        ) const { return data.nr(); }

        long nc (
        ) const { return data.nc(); }

        long size (
        ) const { return data.nr()*data.nc(); }

        template <typename U, size_t len>
        matrix& operator= (
            U (&array)[len]
        )
        {
            COMPILE_TIME_ASSERT(NR*NC == len && len > 0);
            size_t idx = 0;
            for (long r = 0; r < NR; ++r)
            {
                for (long c = 0; c < NC; ++c)
                {
                    data(r,c) = static_cast<T>(array[idx]);
                    ++idx;
                }
            }
            return *this;
        }

        template <typename EXP>
        matrix& operator= (
            const matrix_exp<EXP>& m
        )
        {
            // You get an error on this line if the matrix you are trying to 
            // assign m to is a statically sized matrix and  m's dimensions don't 
            // match that of *this. 
            COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0);
            COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0);
            DLIB_ASSERT((NR == 0 || nr() == m.nr()) && 
                   (NC == 0 || nc() == m.nc()), 
                "\tmatrix& matrix::operator=(const matrix_exp& m)"
                << "\n\tYou are trying to assign a dynamically sized matrix to a statically sized matrix with the wrong size"
                << "\n\tnr():   " << nr()
                << "\n\tnc():   " << nc()
                << "\n\tm.nr(): " << m.nr()
                << "\n\tm.nc(): " << m.nc()
                << "\n\tthis:   " << this
                );

            // You get an error on this line if the matrix m contains a type that isn't
            // the same as the type contained in the target matrix.
            COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true) ||
                                (is_matrix<typename EXP::type>::value == true));
            if (m.destructively_aliases(*this) == false)
            {
                set_size(m.nr(),m.nc());
                matrix_assign(*this, m);
            }
            else
            {
                // we have to use a temporary matrix object here because
                // *this is aliased inside the matrix_exp m somewhere.
                matrix temp;
                temp.set_size(m.nr(),m.nc());
                matrix_assign(temp, m);
                temp.swap(*this);
            }
            return *this;
        }

        template <typename EXP>
        matrix& operator += (
            const matrix_exp<EXP>& m
        )
        {
            // The matrix you are trying to assign m to is a statically sized matrix and 
            // m's dimensions don't match that of *this. 
            COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0);
            COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0);
            DLIB_ASSERT(this->nr() == m.nr() && this->nc() == m.nc(), 
                "\tmatrix& matrix::operator+=(const matrix_exp& m)"
                << "\n\tYou are trying to add a dynamically sized matrix to a statically sized matrix with the wrong size"
                << "\n\tthis->nr(): " << nr()
                << "\n\tthis->nc(): " << nc()
                << "\n\tm.nr():     " << m.nr()
                << "\n\tm.nc():     " << m.nc()
                << "\n\tthis:       " << this
                );
            COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true));
            if (m.destructively_aliases(*this) == false)
            {
                matrix_assign(*this, m + *this);
            }
            else
            {
                // we have to use a temporary matrix object here because
                // this->data is aliased inside the matrix_exp m somewhere.
                matrix temp;
                temp.set_size(m.nr(),m.nc());
                matrix_assign(temp, m + *this);
                temp.swap(*this);
            }
            return *this;
        }


        template <typename EXP>
        matrix& operator -= (
            const matrix_exp<EXP>& m
        )
        {
            // The matrix you are trying to assign m to is a statically sized matrix and 
            // m's dimensions don't match that of *this. 
            COMPILE_TIME_ASSERT(EXP::NR == NR || NR == 0 || EXP::NR == 0);
            COMPILE_TIME_ASSERT(EXP::NC == NC || NC == 0 || EXP::NC == 0);
            DLIB_ASSERT(this->nr() == m.nr() && this->nc() == m.nc(), 
                "\tmatrix& matrix::operator-=(const matrix_exp& m)"
                << "\n\tYou are trying to subtract a dynamically sized matrix from a statically sized matrix with the wrong size"
                << "\n\tthis->nr(): " << nr()
                << "\n\tthis->nc(): " << nc()
                << "\n\tm.nr():     " << m.nr()
                << "\n\tm.nc():     " << m.nc()
                << "\n\tthis:       " << this
                );
            COMPILE_TIME_ASSERT((is_same_type<typename EXP::type,type>::value == true));
            if (m.destructively_aliases(*this) == false)
            {
                matrix_assign(*this, *this - m);
            }
            else
            {
                // we have to use a temporary matrix object here because
                // this->data is aliased inside the matrix_exp m somewhere.
                matrix temp;
                temp.set_size(m.nr(),m.nc());
                matrix_assign(temp, *this - m);
                temp.swap(*this);
            }
            return *this;
        }

        matrix& operator += (
            const matrix& m
        )
        {
            const long size = m.nr()*m.nc();
            for (long i = 0; i < size; ++i)
                data(i) += m.data(i);
            return *this;
        }

        matrix& operator -= (
            const matrix& m
        )
        {
            const long size = m.nr()*m.nc();
            for (long i = 0; i < size; ++i)
                data(i) -= m.data(i);
            return *this;
        }

        matrix& operator *= (
            const T& a
        )
        {
            const long size = data.nr()*data.nc();
            for (long i = 0; i < size; ++i)
                data(i) *= a;
            return *this;
        }

        matrix& operator /= (
            const T& a
        )
        {
            const long size = data.nr()*data.nc();
            for (long i = 0; i < size; ++i)
                data(i) /= a;
            return *this;
        }

        matrix& operator= (
            const matrix& m
        )
        {
            if (this != &m)
            {
                set_size(m.nr(),m.nc());
                const long size = m.nr()*m.nc();
                for (long i = 0; i < size; ++i)
                    data(i) = m.data(i);
            }
            return *this;
        }

        void swap (
            matrix& item
        )
        {
            data.swap(item.data);
        }

    private:
        struct literal_assign_helper
        {
            /*
                This struct is a helper struct returned by the operator<<() function below.  It is
                used primarily to enable us to put DLIB_CASSERT statements on the usage of the
                operator<< form of matrix assignment.
            */

            literal_assign_helper(const literal_assign_helper& item) : m(item.m), r(item.r), c(item.c), is_copy(true) {}
            literal_assign_helper(matrix* m_): m(m_), r(0), c(0),is_copy(false) {}
            ~literal_assign_helper()
            {
                DLIB_CASSERT(!is_copy || r == m->nr(),
                             "You have used the matrix operator<< assignment incorrectly by failing to\n"
                             "supply a full set of values for every element of a matrix object.\n");
            }

            template <typename U>
            const literal_assign_helper& operator, (
                const U& val
            ) const
            {
                DLIB_CASSERT(r < m->nr() && c < m->nc(),
                             "You have used the matrix operator<< assignment incorrectly by attempting to\n" <<
                             "supply more values than there are elements in the matrix object being assigned to.\n\n" <<
                             "Did you forget to call set_size()?\n");
                (*m)(r,c) = static_cast<T>(val);
                ++c;
                if (c == m->nc())
                {
                    c = 0;
                    ++r;
                }
                return *this;
            }


            matrix* m;
            mutable long r;
            mutable long c;
            bool is_copy;
        };

    public:

        template <typename U>
        const literal_assign_helper operator << (
            const U& val
        ) { return literal_assign_helper(this) , val; }

    private:


        typename layout::template layout<T,NR,NC,mem_manager> data;
    };

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

    template <
        typename T,
        long NR,
        long NC,
        typename mm,
        typename l
        >
    void swap(
        matrix<T,NR,NC,mm,l>& a,
        matrix<T,NR,NC,mm,l>& b
    ) { a.swap(b); }

    template <
        typename T,
        long NR,
        long NC,
        typename mm,
        typename l
        >
    void serialize (
        const matrix<T,NR,NC,mm,l>& item, 
        std::ostream& out
    )
    {
        try
        {
            serialize(item.nr(),out);
            serialize(item.nc(),out);
            for (long r = 0; r < item.nr(); ++r)
            {
                for (long c = 0; c < item.nc(); ++c)
                {
                    serialize(item(r,c),out);
                }
            }
        }
        catch (serialization_error& e)
        {
            throw serialization_error(e.info + "\n   while serializing dlib::matrix");
        }
    }

    template <
        typename T,
        long NR,
        long NC,
        typename mm,
        typename l
        >
    void deserialize (
        matrix<T,NR,NC,mm,l>& item, 
        std::istream& in
    )
    {
        try
        {
            long nr, nc;
            deserialize(nr,in); 
            deserialize(nc,in); 

            if (NR != 0 && nr != NR)
                throw serialization_error("Error while deserializing a dlib::matrix.  Invalid rows");
            if (NC != 0 && nc != NC)
                throw serialization_error("Error while deserializing a dlib::matrix.  Invalid columns");

            item.set_size(nr,nc);
            for (long r = 0; r < nr; ++r)
            {
                for (long c = 0; c < nc; ++c)
                {
                    deserialize(item(r,c),in);
                }
            }
        }
        catch (serialization_error& e)
        {
            throw serialization_error(e.info + "\n   while deserializing a dlib::matrix");
        }
    }

    template <
        typename EXP
        >
    std::ostream& operator<< (
        std::ostream& out,
        const matrix_exp<EXP>& m
    )
    {
        using namespace std;
        const streamsize old = out.width();

        // first figure out how wide we should make each field
        string::size_type w = 0;
        ostringstream sout;
        for (long r = 0; r < m.nr(); ++r)
        {
            for (long c = 0; c < m.nc(); ++c)
            {
                sout << m(r,c); 
                w = std::max(sout.str().size(),w);
                sout.str("");
            }
        }

        // now actually print it
        for (long r = 0; r < m.nr(); ++r)
        {
            for (long c = 0; c < m.nc(); ++c)
            {
                out.width(static_cast<streamsize>(w));
                out << m(r,c) << " ";
            }
            out << "\n";
        }
        out.width(old);
        return out;
    }

// ----------------------------------------------------------------------------------------

}

#ifdef _MSC_VER
// put that warning back to its default setting
#pragma warning(default : 4355)
#endif

#endif // DLIB_MATRIx_