dlib/examples/matrix_expressions_ex.cpp

// The contents of this file are in the public domain. See LICENSE_FOR_EXAMPLE_PROGRAMS.txt

/*
    This example contains a detailed discussion of the template expression
    technique used to implement the matrix tools in the dlib C++ library.

    It also gives examples showing how a user can create their own custom
    matrix expressions.

    Note that you should be familiar with the dlib::matrix before reading
    this example.  So if you haven't done so already you should read the
    matrix_ex.cpp example program.
*/


#include <iostream>
#include <dlib/matrix.h>

using namespace dlib;
using namespace std;

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

void custom_matrix_expressions_example();

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

int main()
{

    // Declare some variables used below
    matrix<double,3,1> y;
    matrix<double,3,3> M;
    matrix<double> x;

    // set all elements to 1
    y = 1;
    M = 1;
    

    // -------------------------  Template Expressions -----------------------------
    // Now I will discuss the "template expressions" technique and how it is 
    // used in the matrix object.  First consider the following expression:
    x = y + y;

    /*
        Normally this expression results in machine code that looks, at a high 
        level, like the following:
            temp = y + y;
            x = temp

        Temp is a temporary matrix returned by the overloaded + operator.  
        temp then contains the result of adding y to itself.  The assignment 
        operator copies the value of temp into x and temp is then destroyed while 
        the blissful C++ user never sees any of this.

        This is, however, totally inefficient.  In the process described above 
        you have to pay for the cost of constructing a temporary matrix object
        and allocating its memory.  Then you pay the additional cost of copying
        it over to x.  It also gets worse when you have more complex expressions
        such as x = round(y + y + y + M*y) which would involve the creation and copying 
        of 5 temporary matrices.  
        
        All these inefficiencies are removed by using the template expressions 
        technique.  The basic idea is as follows, instead of having operators and 
        functions return temporary matrix objects you return a special object that 
        represents the expression you wish to perform.

        So consider the expression x = y + y again.  With dlib::matrix what happens
        is the expression y+y returns a matrix_exp object instead of a temporary matrix.
        The construction of a matrix_exp does not allocate any memory or perform any 
        computations.  The matrix_exp however has an interface that looks just like a 
        dlib::matrix object and when you ask it for the value of one of its elements 
        it computes that value on the spot.  Only in the assignment operator does
        someone ask the matrix_exp for these values so this avoids the use of any
        temporary matrices.  Thus the statement x = y + y is equivalent to the following 
        code:
            // loop over all elements in y matrix
            for (long r = 0; r < y.nr(); ++r)
                for (long c = 0; c < y.nc(); ++c)
                    x(r,c) = y(r,c) + y(r,c);  
                
       
        This technique works for expressions of arbitrary complexity.  So if you typed
        x = round(y + y + y + M*y) it would involve no temporary matrices being created
        at all.  Each operator takes and returns only matrix_exp objects.  Thus, no
        computations are performed until the assignment operator requests the values
        from the matrix_exp it receives as input.  This also means that statements such as:
            auto x = round(y + y + y + M*y)
        will not work properly because x would be a matrix expression that references
        parts of the expression round(y + y + y + M*y) but those expression parts will
        immediately go out of scope so x will contain references to non-existing sub
        matrix expressions.  This is very bad, so you should never use auto to store
        the result of a matrix expression.  Always store the output in a matrix object
        like so:
            matrix<double> x = round(y + y + y + M*y)




        In terms of implementation, there is a slight complication in all of this.  It
        is for statements that involve the multiplication of a complex matrix_exp such
        as the following:
    */
        x = M*(M+M+M+M+M+M+M);
    /*
        According to the discussion above, this statement would compute the value of 
        M*(M+M+M+M+M+M+M) totally without any temporary matrix objects.  This sounds 
        good but we should take a closer look.  Consider that the + operator is 
        invoked 6 times.  This means we have something like this:

        x = M * (matrix_exp representing M+M+M+M+M+M+M);

        M is being multiplied with a quite complex matrix_exp.  Now recall that when 
        you ask a matrix_exp what the value of any of its elements are it computes 
        their values *right then*.  
        
        If you think on what is involved in performing a matrix multiply you will 
        realize that each element of a matrix is accessed M.nr() times.  In the 
        case of our above expression the cost of accessing an element of the 
        matrix_exp on the right hand side is the cost of doing 6 addition operations. 

        Thus, it would be faster to assign M+M+M+M+M+M+M to a temporary matrix and then
        multiply that by M.  This is exactly what the dlib::matrix does under the covers.  
        This is because it is able to spot expressions where the introduction of a 
        temporary is needed to speed up the computation and it will automatically do this 
        for you.  



        
        Another complication that is dealt with automatically is aliasing.  All matrix 
        expressions are said to "alias" their contents.  For example, consider
        the following expressions:
           M + M
           M * M 

        We say that the expressions (M + M) and (M * M) alias M.  Additionally, we say that
        the expression (M * M) destructively aliases M.    

        To understand why we say (M * M) destructively aliases M consider what would happen
        if we attempted to evaluate it without first assigning (M * M) to a temporary matrix.
        That is, if we attempted to perform the following:
            for (long r = 0; r < M.nr(); ++r)
                for (long c = 0; c < M.nc(); ++c)
                    M(r,c) = rowm(M,r)*colm(M,c);  

        It is clear that the result would be corrupted and M wouldn't end up with the right
        values in it.  So in this case we must perform the following:
            temp = M*M;
            M = temp;

        This sort of interaction is what defines destructive aliasing.  Whenever we are
        assigning a matrix expression to a destination that is destructively aliased by
        the expression we need to introduce a temporary.  The dlib::matrix is capable of 
        recognizing the two forms of aliasing and introduces temporary matrices only when 
        necessary.
    */



    // Next we discuss how to create custom matrix expressions.   In what follows we 
    // will define three different matrix expressions and show their use.  
    custom_matrix_expressions_example();
}

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

template <typename M>
struct example_op_trans 
{
    /*!
        This object defines a matrix expression that represents a transposed matrix.
        As discussed above, constructing this object doesn't compute anything.  It just
        holds a reference to a matrix and presents an interface which defines
        matrix transposition.   
    !*/

    // Here we simply hold a reference to the matrix we are supposed to be the transpose of.
    example_op_trans( const M& m_) : m(m_){}
    const M& m;

    // The cost field is used by matrix multiplication code to decide if a temporary needs to 
    // be introduced.  It represents the computational cost of evaluating an element of the
    // matrix expression.  In this case we say that the cost of obtaining an element of the
    // transposed matrix is the same as obtaining an element of the original matrix (since
    // transpose doesn't really compute anything).
    const static long cost = M::cost;

    // Here we define the matrix expression's compile-time known dimensions.  Since this
    // is a transpose we define them to be the reverse of M's dimensions.
    const static long NR = M::NC;
    const static long NC = M::NR;

    // Define the type of element in this matrix expression.  Also define the 
    // memory manager type used and the layout.  Usually we use the same types as the
    // input matrix.
    typedef typename M::type type;
    typedef typename M::mem_manager_type mem_manager_type;
    typedef typename M::layout_type layout_type;

    // This is where the action is.  This function is what defines the value of an element of
    // this matrix expression.  Here we are saying that m(c,r) == trans(m)(r,c) which is just
    // the definition of transposition.  Note also that we must define the return type from this
    // function as a typedef.  This typedef lets us either return our argument by value or by
    // reference.  In this case we use the same type as the underlying m matrix.  Later in this
    // example program you will see two other options.
    typedef typename M::const_ret_type const_ret_type;
    const_ret_type apply (long r, long c) const { return m(c,r); }

    // Define the run-time defined dimensions of this matrix.  
    long nr () const { return m.nc(); }
    long nc () const { return m.nr(); }

    // Recall the discussion of aliasing.  Each matrix expression needs to define what
    // kind of aliasing it introduces so that we know when to introduce temporaries.  The 
    // aliases() function indicates that the matrix transpose expression aliases item if
    // and only if the m matrix aliases item.  
    template <typename U> bool aliases               ( const matrix_exp<U>& item) const { return m.aliases(item); }
    // This next function indicates that the matrix transpose expression also destructively 
    // aliases anything m aliases.  I.e. transpose has destructive aliasing.
    template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.aliases(item); }

}; 


// Here we define a simple function that creates and returns transpose expressions.  Note that the
// matrix_op<> template is a matrix_exp object and exists solely to reduce the amount of boilerplate
// you have to write to create a matrix expression.
template < typename M >
const matrix_op<example_op_trans<M> > example_trans (
    const matrix_exp<M>& m
)
{
    typedef example_op_trans<M> op;
    // m.ref() returns a reference to the object of type M contained in the matrix expression m.
    return matrix_op<op>(op(m.ref()));
}

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

template <typename T>
struct example_op_vector_to_matrix  
{
    /*!
        This object defines a matrix expression that holds a reference to a std::vector<T>
        and makes it look like a column vector.  Thus it enables you to use a std::vector
        as if it was a dlib::matrix.

    !*/
    example_op_vector_to_matrix( const std::vector<T>& vect_) : vect(vect_){}

    const std::vector<T>& vect;

    // This expression wraps direct memory accesses so we use the lowest possible cost. 
    const static long cost = 1; 

    const static long NR = 0; // We don't know the length of the vector until runtime.  So we put 0 here.
    const static long NC = 1; // We do know that it only has one column (since it's a vector)
    typedef T type;
    // Since the std::vector doesn't use a dlib memory manager we list the default one here.
    typedef default_memory_manager mem_manager_type;
    // The layout type also doesn't really matter in this case.  So we list row_major_layout
    // since it is a good default.
    typedef row_major_layout layout_type;

    // Note that we define const_ret_type to be a reference type.  This way we can
    // return the contents of the std::vector by reference.
    typedef const T& const_ret_type;
    const_ret_type apply (long r, long ) const { return vect[r]; }

    long nr () const { return vect.size(); }
    long nc () const { return 1; }

    // This expression never aliases anything since it doesn't contain any matrix expression (it 
    // contains only a std::vector which doesn't count since you can't assign a matrix expression
    // to a std::vector object).
    template <typename U> bool aliases               ( const matrix_exp<U>& ) const { return false; }
    template <typename U> bool destructively_aliases ( const matrix_exp<U>& ) const { return false; }
}; 

template < typename T >
const matrix_op<example_op_vector_to_matrix<T> > example_vector_to_matrix (
    const std::vector<T>& vector
)
{
    typedef example_op_vector_to_matrix<T> op;
    return matrix_op<op>(op(vector));
}

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

template <typename M, typename T>
struct example_op_add_scalar
{
    /*!
        This object defines a matrix expression that represents a matrix with a single
        scalar value added to all its elements.  
    !*/

    example_op_add_scalar( const M& m_, const T& val_) : m(m_), val(val_){}

    // A reference to the matrix 
    const M& m;
    // A copy of the scalar value that should be added to each element of m
    const T val;

    // This time we add 1 to the cost since evaluating an element of this 
    // expression means performing 1 addition operation.
    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 typename M::mem_manager_type mem_manager_type;
    typedef typename M::layout_type layout_type;

    // Note that we declare const_ret_type to be a non-reference type.  This is important
    // since apply() computes a new temporary value and thus we can't return a reference
    // to it.
    typedef type const_ret_type; 
    const_ret_type apply (long r, long c) const { return m(r,c) + val; }

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

    // This expression aliases anything m aliases.
    template <typename U> bool aliases               ( const matrix_exp<U>& item) const { return m.aliases(item); }
    // Unlike the transpose expression.  This expression only destructively aliases something if m does. 
    // So this expression has the regular non-destructive kind of aliasing.
    template <typename U> bool destructively_aliases ( const matrix_exp<U>& item) const { return m.destructively_aliases(item); }

}; 

template < typename M, typename T >
const matrix_op<example_op_add_scalar<M,T> > add_scalar (
    const matrix_exp<M>& m,
    T val
)
{
    typedef example_op_add_scalar<M,T> op;
    return matrix_op<op>(op(m.ref(), val));
}

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

void custom_matrix_expressions_example(
)
{
    matrix<double> x(2,3);
    x = 1, 1, 1,
        2, 2, 2;

    cout << x << endl;

    // Finally, lets use the matrix expressions we defined above.

    // prints the transpose of x
    cout << example_trans(x) << endl;

    // prints this:
    //   11 11 11
    //   12 12 12
    cout << add_scalar(x, 10) << endl;


    // matrix expressions can be nested, even the user defined ones.
    // the following statement prints this:
    //   11 12 
    //   11 12 
    //   11 12 
    cout << example_trans(add_scalar(x, 10)) << endl;

    // Since we setup the alias detection correctly we can even do this:
    x = example_trans(add_scalar(x, 10));
    cout << "new x:\n" << x << endl;

    cout << "Do some testing with the example_vector_to_matrix() function: " << endl;
    std::vector<float> vect;
    vect.push_back(1);
    vect.push_back(3);
    vect.push_back(5);

    // Now lets treat our std::vector like a matrix and print some things.
    cout << example_vector_to_matrix(vect) << endl;
    cout << add_scalar(example_vector_to_matrix(vect), 10) << endl;



    /*
        As an aside, note that dlib contains functions equivalent to the ones we 
        defined above.  They are:
            - dlib::trans()
            - dlib::mat() (converts things into matrices)
            - operator+ (e.g. you can say my_mat + 1)


        Also, if you are going to be creating your own matrix expression you should also
        look through the matrix code in the dlib/matrix folder.  There you will find 
        many other examples of matrix expressions. 
    */
}

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