Add expression templates (#4)

  * using underlying Eigen expression templates
  * scalar operations with expression templates
  * save DataContainer meta info in expressions
  * add unary operators
  * adding enum for DataHandlerMapCPU type
  * added expression templates readme
  * removed operators between DataContainers, scalars and DataHandlers
  * in-place operations using expressions
  * test cases
  * benchmark script

CMakeLists.txt

@@ -19,6 +19,7 @@ endif ()
 option(ELSA_TESTING "Enable building the unit tests" ${ELSA_MASTER_PROJECT})
 option(ELSA_CREATE_JUNIT_REPORTS "Enable creating JUnit style reports when running tests" ON)
 option(ELSA_COVERAGE "Enable test coverage computation and reporting" OFF)
+option(ELSA_BENCHMARKS "Enable elsa benchmark test cases" OFF)
 
 option(GIT_SUBMODULE "Enable updating the submodules during build" ${ELSA_MASTER_PROJECT})
 option(ELSA_INSTALL "Enable generating the install targets for make install" ${ELSA_MASTER_PROJECT})
@@ -87,7 +88,7 @@ if(ELSA_TESTING)
     enable_testing()
     add_subdirectory(thirdparty/Catch2)
 
-	add_custom_target(tests
+    add_custom_target(tests
         COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure
         WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
         COMMENT "Build and run all the tests.")

CONTRIBUTING.md

@@ -45,6 +45,18 @@ We use a testing style described as [Behaviour-Driven
 Development](https://github.com/catchorg/Catch2/blob/master/docs/tutorial.md#bdd-style) (BDD). Please
 follow this style when adding new tests.
 
+## Benchmarking
+You can use the catch testing framework to do [benchmarking
+](https://github.com/catchorg/Catch2/blob/master/docs/benchmarks.md). If so, add your benchmarking
+case following this template
+```cmake
+if(ELSA_BENCHMARKS)
+    ELSA_TEST(BenchmarkName)
+endif()
+```
+which ensures that the test case is only registered and build if the cmake option was
+enabled.
+
 ## Style Guide
 We use the tool `clang-format` to autoformat our code with the [given style
 file](.clang-format). Please make sure that your code is formatted accordingly, otherwise the CI

docs/modules/core.rst

@@ -9,7 +9,6 @@ Core Type Declarations
 
 .. doxygenfile:: elsa.h
 
-
 DataDescriptor
 ==============
 
@@ -23,17 +22,14 @@ BlockDescriptor
 DataContainer
 =============
 
-
 .. doxygenclass:: elsa::DataContainer
 
 LinearOperator
 ==============
 
-
 .. doxygenclass:: elsa::LinearOperator
 
 
-
 Implementation Details
 ======================
 
@@ -42,7 +38,6 @@ Cloneable
 
 .. doxygenclass:: elsa::Cloneable
 
-
 DataHandler
 -----------
 
@@ -52,3 +47,9 @@ DataHandlerCPU
 --------------
 
 .. doxygenclass:: elsa::DataHandlerCPU
+
+Expression
+----------
+.. mdinclude:: expression_templates.md
+
+.. doxygenclass:: elsa::Expression

docs/modules/expression_templates.md

+In elsa, we are using expression templates for fast and efficient mathematical operations while at the same time having intuitive syntax. This technique is also known as lazy-evaluation which means that computations are delayed until the results are actually needed.
+
+Take the example
+```cpp
+auto expression = dataContainer1 + dataContainer2;
+```
+where the type of `expression` is `elsa::Expression<operator+, DataContainer<real_t>, DataContainer<real_t>>`.
+
+The `expression` contains the work *to be done* rather than the actual result. Only when assigning to a `DataContainer` or constructing a new DataContainer, the expression gets evaluated automatically
+```cpp
+DataContainer result = dataContainer1 + dataContainer2;
+```
+into `result`.
+
+Nesting is easily possible like
+```cpp
+auto expression = dataContainer1 + dataContainer2 * expression;
+```
+which will store an expression tree as the type information.
+
+Please note also that the operators/member functions available on an `Expression` type are different from a `DataContainer`. 
+
+An expression has a member function `eval()`. However, calling `eval()` will *not* return a `DataContainer` but depending on the current `DataHandler` the underlying raw data computation result.
+
+If single element-wise access is necessary, it is possible to call `expression.eval()[index]`. Note that this is computational very expensive as the whole expression gets evaluated at every index individually.

elsa/core/CMakeLists.txt

@@ -17,7 +17,9 @@ set(MODULE_HEADERS
         DataHandler.h
         DataHandlerCPU.h
         DataHandlerMapCPU.h
-        LinearOperator.h)
+        LinearOperator.h
+        Expression.h
+        ExpressionPredicates.h)
 
 # list all the code files of the module
 set(MODULE_SOURCES
@@ -65,4 +67,4 @@ endif(ELSA_TESTING)
 registerComponent(${ELSA_MODULE_NAME})
 
 # install the module
-InstallElsaModule(${ELSA_MODULE_NAME} ${ELSA_MODULE_TARGET_NAME} ${ELSA_MODULE_EXPORT_TARGET})
\ No newline at end of file
+InstallElsaModule(${ELSA_MODULE_NAME} ${ELSA_MODULE_TARGET_NAME} ${ELSA_MODULE_EXPORT_TARGET})

elsa/core/DataContainer.cpp

@@ -13,7 +13,8 @@ namespace elsa
     DataContainer<data_t>::DataContainer(const DataDescriptor& dataDescriptor,
                                          DataHandlerType handlerType)
         : _dataDescriptor{dataDescriptor.clone()},
-          _dataHandler{createDataHandler(handlerType, _dataDescriptor->getNumberOfCoefficients())}
+          _dataHandler{createDataHandler(handlerType, _dataDescriptor->getNumberOfCoefficients())},
+          _dataHandlerType{handlerType}
     {
     }
 
@@ -22,7 +23,8 @@ namespace elsa
                                          const Eigen::Matrix<data_t, Eigen::Dynamic, 1>& data,
                                          DataHandlerType handlerType)
         : _dataDescriptor{dataDescriptor.clone()},
-          _dataHandler{createDataHandler(handlerType, _dataDescriptor->getNumberOfCoefficients())}
+          _dataHandler{createDataHandler(handlerType, _dataDescriptor->getNumberOfCoefficients())},
+          _dataHandlerType{handlerType}
     {
         if (_dataHandler->getSize() != data.size())
             throw std::invalid_argument("DataContainer: initialization vector has invalid size");
@@ -33,7 +35,9 @@ namespace elsa
 
     template <typename data_t>
     DataContainer<data_t>::DataContainer(const DataContainer<data_t>& other)
-        : _dataDescriptor{other._dataDescriptor->clone()}, _dataHandler{other._dataHandler->clone()}
+        : _dataDescriptor{other._dataDescriptor->clone()},
+          _dataHandler{other._dataHandler->clone()},
+          _dataHandlerType{other._dataHandlerType}
     {
     }
 
@@ -48,6 +52,8 @@ namespace elsa
             } else {
                 _dataHandler = other._dataHandler->clone();
             }
+
+            _dataHandlerType = other._dataHandlerType;
         }
 
         return *this;
@@ -56,7 +62,8 @@ namespace elsa
     template <typename data_t>
     DataContainer<data_t>::DataContainer(DataContainer<data_t>&& other) noexcept
         : _dataDescriptor{std::move(other._dataDescriptor)},
-          _dataHandler{std::move(other._dataHandler)}
+          _dataHandler{std::move(other._dataHandler)},
+          _dataHandlerType{std::move(other._dataHandlerType)}
     {
         // leave other in a valid state
         other._dataDescriptor = nullptr;
@@ -67,12 +74,15 @@ namespace elsa
     DataContainer<data_t>& DataContainer<data_t>::operator=(DataContainer<data_t>&& other)
     {
         _dataDescriptor = std::move(other._dataDescriptor);
+
         if (_dataHandler) {
             *_dataHandler = std::move(*other._dataHandler);
         } else {
             _dataHandler = std::move(other._dataHandler);
         }
 
+        _dataHandlerType = std::move(other._dataHandlerType);
+
         // leave other in a valid state
         other._dataDescriptor = nullptr;
         other._dataHandler = nullptr;
@@ -147,30 +157,6 @@ namespace elsa
         return _dataHandler->sum();
     }
 
-    template <typename data_t>
-    DataContainer<data_t> DataContainer<data_t>::square() const
-    {
-        return DataContainer<data_t>(*_dataDescriptor, _dataHandler->square());
-    }
-
-    template <typename data_t>
-    DataContainer<data_t> DataContainer<data_t>::sqrt() const
-    {
-        return DataContainer<data_t>(*_dataDescriptor, _dataHandler->sqrt());
-    }
-
-    template <typename data_t>
-    DataContainer<data_t> DataContainer<data_t>::exp() const
-    {
-        return DataContainer<data_t>(*_dataDescriptor, _dataHandler->exp());
-    }
-
-    template <typename data_t>
-    DataContainer<data_t> DataContainer<data_t>::log() const
-    {
-        return DataContainer<data_t>(*_dataDescriptor, _dataHandler->log());
-    }
-
     template <typename data_t>
     DataContainer<data_t>& DataContainer<data_t>::operator+=(const DataContainer<data_t>& dc)
     {
@@ -242,6 +228,8 @@ namespace elsa
         switch (handlerType) {
             case DataHandlerType::CPU:
                 return std::make_unique<DataHandlerCPU<data_t>>(std::forward<Args>(args)...);
+            case DataHandlerType::MAP_CPU:
+                return std::make_unique<DataHandlerCPU<data_t>>(std::forward<Args>(args)...);
             default:
                 throw std::invalid_argument("DataContainer: unknown handler type");
         }
@@ -249,8 +237,11 @@ namespace elsa
 
     template <typename data_t>
     DataContainer<data_t>::DataContainer(const DataDescriptor& dataDescriptor,
-                                         std::unique_ptr<DataHandler<data_t>> dataHandler)
-        : _dataDescriptor{dataDescriptor.clone()}, _dataHandler{std::move(dataHandler)}
+                                         std::unique_ptr<DataHandler<data_t>> dataHandler,
+                                         DataHandlerType dataType)
+        : _dataDescriptor{dataDescriptor.clone()},
+          _dataHandler{std::move(dataHandler)},
+          _dataHandlerType{dataType}
     {
     }
 
@@ -286,7 +277,8 @@ namespace elsa
         const auto& ithDesc = blockDesc->getDescriptorOfBlock(i);
         index_t blockSize = ithDesc.getNumberOfCoefficients();
 
-        return DataContainer<data_t>{ithDesc, _dataHandler->getBlock(startIndex, blockSize)};
+        return DataContainer<data_t>{ithDesc, _dataHandler->getBlock(startIndex, blockSize),
+                                     DataHandlerType::MAP_CPU};
     }
 
     template <typename data_t>
@@ -305,7 +297,8 @@ namespace elsa
 
         // getBlock() returns a pointer to non-const DH, but that's fine as it gets wrapped in a
         // constant container
-        return DataContainer<data_t>{ithDesc, _dataHandler->getBlock(startIndex, blockSize)};
+        return DataContainer<data_t>{ithDesc, _dataHandler->getBlock(startIndex, blockSize),
+                                     DataHandlerType::MAP_CPU};
     }
 
     template <typename data_t>
@@ -314,7 +307,8 @@ namespace elsa
         if (dataDescriptor.getNumberOfCoefficients() != getSize())
             throw std::invalid_argument("DataContainer: view must have same size as container");
 
-        return DataContainer<data_t>{dataDescriptor, _dataHandler->getBlock(0, getSize())};
+        return DataContainer<data_t>{dataDescriptor, _dataHandler->getBlock(0, getSize()),
+                                     DataHandlerType::MAP_CPU};
     }
 
     template <typename data_t>
@@ -325,7 +319,8 @@ namespace elsa
 
         // getBlock() returns a pointer to non-const DH, but that's fine as it gets wrapped in a
         // constant container
-        return DataContainer<data_t>{dataDescriptor, _dataHandler->getBlock(0, getSize())};
+        return DataContainer<data_t>{dataDescriptor, _dataHandler->getBlock(0, getSize()),
+                                     DataHandlerType::MAP_CPU};
     }
 
     template <typename data_t>
@@ -400,6 +395,28 @@ namespace elsa
         return const_reverse_iterator(cbegin());
     }
 
+    template <typename data_t>
+    typename DataContainer<data_t>::HandlerTypes_t DataContainer<data_t>::getHandlerPtr() const
+    {
+        DataContainer<data_t>::HandlerTypes_t handler;
+
+        if (_dataHandlerType == DataHandlerType::CPU) {
+            handler = static_cast<DataHandlerCPU<data_t>*>(_dataHandler.get());
+        }
+
+        if (_dataHandlerType == DataHandlerType::MAP_CPU) {
+            handler = static_cast<DataHandlerMapCPU<data_t>*>(_dataHandler.get());
+        }
+
+        return handler;
+    }
+
+    template <typename data_t>
+    DataHandlerType DataContainer<data_t>::getDataHandlerType() const
+    {
+        return _dataHandlerType;
+    }
+
     // ------------------------------------------
     // explicit template instantiation
     template class DataContainer<float>;

elsa/core/DataContainer.h

@@ -3,7 +3,10 @@
 #include "elsaDefines.h"
 #include "DataDescriptor.h"
 #include "DataHandler.h"
+#include "DataHandlerCPU.h"
+#include "DataHandlerMapCPU.h"
 #include "DataContainerIterator.h"
+#include "Expression.h"
 
 #include <memory>
 #include <type_traits>
@@ -18,6 +21,7 @@ namespace elsa
      * \author Tobias Lasser - rewrite, modularization, modernization
      * \author David Frank - added DataHandler concept, iterators
      * \author Nikola Dinev - add block support
+     * \author Jens Petit - expression templates
      *
      * \tparam data_t - data type that is stored in the DataContainer, defaulting to real_t.
      *
@@ -25,17 +29,12 @@ namespace elsa
      * be n-dimensional, and will be stored in memory in a linearized fashion. The information
      * on how this linearization is performed is provided by an associated DataDescriptor.
      */
-    template <typename data_t = real_t>
+    template <typename data_t>
     class DataContainer
     {
     public:
-        /**
-         * \brief type of the DataHandler used to store the actual data
-         *
-         * The following handler types are currently supported:
-         *    - CPU: data is stored as an Eigen::Matrix in CPU main memory
-         */
-        enum class DataHandlerType { CPU };
+        /// union of all possible handler raw pointers
+        using HandlerTypes_t = std::variant<DataHandlerCPU<data_t>*, DataHandlerMapCPU<data_t>*>;
 
         /// delete default constructor (without metadata there can be no valid container)
         DataContainer() = delete;
@@ -96,6 +95,37 @@ namespace elsa
          */
         DataContainer<data_t>& operator=(DataContainer<data_t>&& other);
 
+        /**
+         * \brief Expression evaluation assignment for DataContainer
+         *
+         * \param[in] source expression to evaluate
+         *
+         * This evaluates an expression template term into the underlying data member of
+         * the DataHandler in use.
+         */
+        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
+        DataContainer<data_t>& operator=(Source const& source)
+        {
+            _dataHandler->accessData() = source.eval();
+
+            return *this;
+        }
+
+        /**
+         * \brief Expression constructor
+         *
+         * \param[in] source expression to evaluate
+         *
+         * It creates a new DataContainer out of an expression. For this the meta information which
+         * is saved in the expression is used.
+         */
+        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
+        DataContainer<data_t>(Source const& source)
+            : DataContainer<data_t>(source.getDataMetaInfo().first, source.eval(),
+                                    source.getDataMetaInfo().second)
+        {
+        }
+
         /// return the current DataDescriptor
         const DataDescriptor& getDataDescriptor() const;
 
@@ -140,6 +170,13 @@ namespace elsa
         /// return the dot product of this signal with the one from container other
         data_t dot(const DataContainer<data_t>& other) const;
 
+        /// return the dot product of this signal with the one from an expression
+        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
+        data_t dot(const Source& source) const
+        {
+            return (*this * source).eval().sum();
+        }
+
         /// return the squared l2 norm of this signal (dot product with itself)
         data_t squaredL2Norm() const;
 
@@ -152,30 +189,50 @@ namespace elsa
         /// return the sum of all elements of this signal
         data_t sum() const;
 
-        /// return a new DataContainer with element-wise squared values of this one
-        DataContainer<data_t> square() const;
-
-        /// return a new DataContainer with element-wise square roots of this one
-        DataContainer<data_t> sqrt() const;
-
-        /// return a new DataContainer with element-wise exponentials of this one
-        DataContainer<data_t> exp() const;
-
-        /// return a new DataContainer with element-wise logarithms of this one
-        DataContainer<data_t> log() const;
-
         /// compute in-place element-wise addition of another container
         DataContainer<data_t>& operator+=(const DataContainer<data_t>& dc);
 
+        /// compute in-place element-wise addition with another expression
+        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
+        DataContainer<data_t>& operator+=(Source const& source)
+        {
+            *this = *this + source;
+            return *this;
+        }
+
         /// compute in-place element-wise subtraction of another container
         DataContainer<data_t>& operator-=(const DataContainer<data_t>& dc);
 
+        /// compute in-place element-wise subtraction with another expression
+        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
+        DataContainer<data_t>& operator-=(Source const& source)
+        {
+            *this = *this - source;
+            return *this;
+        }
+
         /// compute in-place element-wise multiplication with another container
         DataContainer<data_t>& operator*=(const DataContainer<data_t>& dc);
 
+        /// compute in-place element-wise multiplication with another expression
+        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
+        DataContainer<data_t>& operator*=(Source const& source)
+        {
+            *this = *this * source;
+            return *this;
+        }
+
         /// compute in-place element-wise division by another container
         DataContainer<data_t>& operator/=(const DataContainer<data_t>& dc);
 
+        /// compute in-place element-wise division with another expression
+        template <typename Source, typename = std::enable_if_t<isExpression<Source>>>
+        DataContainer<data_t>& operator/=(Source const& source)
+        {
+            *this = *this / source;
+            return *this;
+        }
+
         /// compute in-place addition of a scalar
         DataContainer<data_t>& operator+=(data_t scalar);
 
@@ -275,12 +332,26 @@ namespace elsa
         /// difference type for iterators
         using difference_type = std::ptrdiff_t;
 
+        /// returns the type of the DataHandler in use
+        DataHandlerType getDataHandlerType() const;
+
+        /// friend constexpr function to implement expression templates
+        template <class Operand, std::enable_if_t<isDataContainer<Operand>, int>>
+        friend constexpr auto evaluateOrReturn(Operand const& operand);
+
     private:
+        /// returns the underlying derived handler as a raw pointer in a std::variant
+        HandlerTypes_t getHandlerPtr() const;
+
         /// the current DataDescriptor
         std::unique_ptr<DataDescriptor> _dataDescriptor;
+
         /// the current DataHandler
         std::unique_ptr<DataHandler<data_t>> _dataHandler;
 
+        /// the current DataHandlerType
+        DataHandlerType _dataHandlerType;
+
         /// factory method to create DataHandlers based on handlerType with perfect forwarding of
         /// constructor arguments
         template <typename... Args>
@@ -289,121 +360,136 @@ namespace elsa
 
         /// private constructor accepting a DataDescriptor and a DataHandler
         explicit DataContainer(const DataDescriptor& dataDescriptor,
-                               std::unique_ptr<DataHandler<data_t>> dataHandler);
+                               std::unique_ptr<DataHandler<data_t>> dataHandler,
+                               DataHandlerType dataType = DataHandlerType::CPU);
     };
 
-    /// element-wise addition of two DataContainers
-    template <typename data_t>
-    DataContainer<data_t> operator+(const DataContainer<data_t>& left,
-                                    const DataContainer<data_t>& right)
-    {
-        DataContainer<data_t> result(left);
-        result += right;
-        return result;
-    }
-
-    /// element-wise subtraction of two DataContainers
-    template <typename data_t>
-    DataContainer<data_t> operator-(const DataContainer<data_t>& left,
-                                    const DataContainer<data_t>& right)
+    /// Multiplying two operands (including scalars)
+    template <typename LHS, typename RHS, typename = std::enable_if_t<isBinaryOpOk<LHS, RHS>>>
+    auto operator*(LHS const& lhs, RHS const& rhs)
     {
-        DataContainer<data_t> result(left);
-        result -= right;
-        return result;
-    }
-
-    /// element-wise multiplication of two DataContainers
-    template <typename data_t>
-    DataContainer<data_t> operator*(const DataContainer<data_t>& left,
-                                    const DataContainer<data_t>& right)
-    {
-        DataContainer<data_t> result(left);
-        result *= right;
-        return result;
-    }
-
-    /// element-wise division of two DataContainers
-    template <typename data_t>
-    DataContainer<data_t> operator/(const DataContainer<data_t>& left,
-                                    const DataContainer<data_t>& right)
-    {
-        DataContainer<data_t> result(left);
-        result /= right;
-        return result;
+        if constexpr (isDcOrExpr<LHS> && isDcOrExpr<RHS>) {
+            auto multiplication = [](auto const& left, auto const& right) {
+                return (left.array() * right.array()).matrix();
+            };
+            return Expression{multiplication, lhs, rhs};
+        } else if constexpr (isArithmetic<LHS>) {
+            auto multiplication = [](auto const& left, auto const& right) {
+                return (left * right.array()).matrix();
+            };
+            return Expression{multiplication, lhs, rhs};
+        } else if constexpr (isArithmetic<RHS>) {
+            auto multiplication = [](auto const& left, auto const& right) {
+                return (left.array() * right).matrix();
+            };
+            return Expression{multiplication, lhs, rhs};
+        } else {
+            auto multiplication = [](auto const& left, auto const& right) { return left * right; };
+            return Expression{multiplication, lhs, rhs};
+        }
     }
 
-    /// addition of DataContainer and scalar
-    template <typename data_t>
-    DataContainer<data_t> operator+(const DataContainer<data_t>& left, data_t right)
+    /// Adding two operands (including scalars)
+    template <typename LHS, typename RHS, typename = std::enable_if_t<isBinaryOpOk<LHS, RHS>>>
+    auto operator+(LHS const& lhs, RHS const& rhs)
     {
-        DataContainer<data_t> result(left);
-        result += right;
-        return result;
+        if constexpr (isDcOrExpr<LHS> && isDcOrExpr<RHS>) {
+            auto addition = [](auto const& left, auto const& right) { return left + right; };
+            return Expression{addition, lhs, rhs};
+        } else if constexpr (isArithmetic<LHS>) {
+            auto addition = [](auto const& left, auto const& right) {
+                return (left + right.array()).matrix();
+            };
+            return Expression{addition, lhs, rhs};
+        } else if constexpr (isArithmetic<RHS>) {
+            auto addition = [](auto const& left, auto const& right) {
+                return (left.array() + right).matrix();
+            };
+            return Expression{addition, lhs, rhs};
+        } else {
+            auto addition = [](auto const& left, auto const& right) { return left + right; };
+            return Expression{addition, lhs, rhs};
+        }
     }
 
-    /// addition of scalar and DataContainer
-    template <typename data_t>
-    DataContainer<data_t> operator+(data_t left, const DataContainer<data_t>& right)
+    /// Subtracting two operands (including scalars)
+    template <typename LHS, typename RHS, typename = std::enable_if_t<isBinaryOpOk<LHS, RHS>>>
+    auto operator-(LHS const& lhs, RHS const& rhs)
     {
-        DataContainer<data_t> result(right);
-        result += left;
-        return result;
-    }
-
-    /// subtraction of scalar from DataContainer
-    template <typename data_t>
-    DataContainer<data_t> operator-(const DataContainer<data_t>& left, data_t right)
-    {
-        DataContainer<data_t> result(left);
-        result -= right;
-        return result;
+        if constexpr (isDcOrExpr<LHS> && isDcOrExpr<RHS>) {
+            auto subtraction = [](auto const& left, auto const& right) { return left - right; };
+            return Expression{subtraction, lhs, rhs};
+        } else if constexpr (isArithmetic<LHS>) {