Moved the discussion of matrix expressions into its own example file. Also

expanded it with examples of how to create new matrix expressions.

--HG--
extra : convert_revision : svn%3Afdd8eb12-d10e-0410-9acb-85c331704f74/trunk%403753

examples/CMakeLists.txt

@@ -47,6 +47,7 @@ add_example(linear_manifold_regularizer_ex)
 add_example(logger_ex)
 add_example(logger_ex_2)
 add_example(matrix_ex)
+add_example(matrix_expressions_ex)
 add_example(member_function_pointer_ex)
 add_example(mlp_ex)
 add_example(model_selection_ex)

examples/matrix_ex.cpp

@@ -3,10 +3,6 @@
 /*
     This is an example illustrating the use of the matrix object 
     from the dlib C++ Library.
-
-    This file also contains a discussion of the template expression
-    technique and how it is used by this library.
-
 */
 
 
@@ -248,113 +244,6 @@ int main()
 
     // MATLAB: var = min(min(A))
     var = min(A);
-
-
-
-    // -------------------------  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 exact details of how the technique is performed are well
-        outside the scope of this example but 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. 
-            
-
-
-
-
-        There is, however, 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.  Consider
-        the following expressions:
-           (1)  M = M + M
-           (2)  B = M * M. 
-           (3)  M = M * M. 
-
-        Expressions (1) and (3) are an example of aliasing and expression (3) is also
-        an example of destructive aliasing.  
-
-        Expression (1) can and does operate without introducing any temporaries even though
-        there is aliasing present in the expression.  The result is loaded straight into M 
-        using the template expression techniques described above.  Expression (2) also 
-        operates without any temporaries being introduced since there isn't any aliasing at all.
-        Expression (3) however contains destructive aliasing.  This is because we can't 
-        change any of the values in the M matrix without corrupting the ultimate result of
-        the matrix multiply.  So we need to introduce a temporary.  These situations are
-        dealt with by dlib::matrix automatically.  Moreover, it can tell the different between
-        simple aliasing and destructive aliasing and will only introduce temporaries when
-        they are necessary.
-    */
 }
 
 // ----------------------------------------------------------------------------------------

examples/matrix_expressions_ex.cpp

@@ -0,0 +1,398 @@
+// 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;
+    x = 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. 
+            
+
+
+
+
+        There is, however, 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 memory_manager<char>::kernel_1a 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 example_trans() and
+        example_vector_to_matrix() already.  They are:
+            - dlib::trans()
+            - dlib::vector_to_matrix()
+
+
+        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. 
+    */
+}
+
+// ----------------------------------------------------------------------------------------
+