Remove default zero init of DataContainer (#50)

For performance reasons this was removed, includes:
- removing functionality from DataHandlerCPU
- added simple benchmark
- adapting clang-tidy to not show warnings from catch2 testing

.clang-tidy

@@ -5,12 +5,14 @@ Checks:          '-*,
                   llvm-namespace-comment,
                   modernize-*,
                   -modernize-use-nodiscard,
+                  -modernize-avoid-c-arrays,
                   readability-*,
                   -readability-else-after-return,
                   -readability-braces-around-statements,
                   -readability-magic-numbers,
                   -readability-uppercase-literal-suffix,
-                  -readability-isolate-declaration'
+                  -readability-isolate-declaration,
+                  -readability-function-size'
 HeaderFilterRegex: 'elsa/elsa'
 AnalyzeTemporaryDtors: false
 CheckOptions:

elsa/core/DataHandlerCPU.cpp

@@ -5,11 +5,9 @@ namespace elsa
 {
 
     template <typename data_t>
-    DataHandlerCPU<data_t>::DataHandlerCPU(index_t size, bool initialize)
+    DataHandlerCPU<data_t>::DataHandlerCPU(index_t size)
         : _data(std::make_shared<DataVector_t>(size))
     {
-        if (initialize)
-            _data->setZero();
     }
 
     template <typename data_t>

elsa/core/DataHandlerCPU.h

@@ -63,15 +63,15 @@ namespace elsa
         ~DataHandlerCPU() override;
 
         /**
-         * \brief Constructor initializing an appropriately sized vector with zeros
+         * \brief Constructor initializing an appropriately sized vector
+         *
+         * Note that the values will not be initialized and therefore contain undefined values.
          *
          * \param[in] size of the vector
-         * \param[in] initialize - set to false if you do not need initialization with zeros
-         * (default: true)
          *
          * \throw std::invalid_argument if the size is non-positive
          */
-        explicit DataHandlerCPU(index_t size, bool initialize = true);
+        explicit DataHandlerCPU(index_t size);
 
         /**
          * \brief Constructor initializing a data vector with a given vector

elsa/core/tests/CMakeLists.txt

@@ -18,4 +18,5 @@ ELSA_TEST(ExpressionTemplates)
 
 if(ELSA_BENCHMARKS)
     ELSA_TEST(BenchmarkExpressionTemplates)
+    ELSA_TEST(BenchmarkDataContainerConstruction)
 endif()

elsa/core/tests/test_BenchmarkDataContainerConstruction.cpp

+/**
+ * \file test_BenchmarkDataContainerConstruction.cpp
+ *
+ * \brief Benchmarks for constructing an empty DataContainer
+ *
+ * \author Jens Petit
+ */
+
+#define CATCH_CONFIG_ENABLE_BENCHMARKING
+
+#include <catch2/catch.hpp>
+#include "DataContainer.h"
+#include <string>
+#include <cstdlib>
+
+using namespace elsa;
+
+TEST_CASE("DataContainer construction benchmark")
+{
+    const index_t dim = 1024;
+    IndexVector_t numCoeff(3);
+    numCoeff << dim, dim, dim;
+    DataDescriptor desc(numCoeff);
+
+    BENCHMARK("only allocating")
+    {
+        DataContainer dc(desc);
+        return dc;
+    };
+
+    BENCHMARK("allocating and assigning")
+    {
+        DataContainer dc(desc);
+        return dc = 0;
+    };
+
+    BENCHMARK("Eigen directly only allocating ")
+    {
+        Eigen::Matrix<float, Eigen::Dynamic, 1> vec(dim * dim * dim);
+        return vec;
+    };
+
+    BENCHMARK("Eigen directly allocating and assigning")
+    {
+        Eigen::Matrix<float, Eigen::Dynamic, 1> vec(dim * dim * dim);
+        return vec.setZero();
+    };
+}

elsa/core/tests/test_DataContainer.cpp

@@ -147,10 +147,6 @@ SCENARIO("Element-wise access of DataContainers")
 
             THEN("it works as expected when using indices/coordinates")
             {
-                REQUIRE(dc[index] == 0.0_a);
-                REQUIRE(dc(coord) == 0.0_a);
-                REQUIRE(dc(17, 4) == 0.0_a);
-
                 dc[index] = 2.2f;
                 REQUIRE(dc[index] == 2.2_a);
                 REQUIRE(dc(coord) == 2.2_a);

elsa/core/tests/test_DataHandlers.cpp

@@ -30,17 +30,11 @@ TEMPLATE_PRODUCT_TEST_CASE("Scenario: Constructing DataHandler", "", (DataHandle
     {
         index_t size = 314;
 
-        WHEN("constructing with zeros")
+        WHEN("constructing")
         {
             const TestType dh{size};
 
             THEN("it has the correct size") { REQUIRE(size == dh.getSize()); }
-
-            THEN("it has a zero vector")
-            {
-                for (index_t i = 0; i < size; ++i)
-                    REQUIRE(dh[i] == static_cast<data_t>(0));
-            }
         }
 
         WHEN("constructing with a given vector")
@@ -55,7 +49,9 @@ TEMPLATE_PRODUCT_TEST_CASE("Scenario: Constructing DataHandler", "", (DataHandle
 
         WHEN("copy constructing")
         {
-            const TestType dh{size};
+            Eigen::Matrix<data_t, Eigen::Dynamic, 1> randVec{size};
+            randVec.setRandom();
+            const TestType dh{randVec};
             const auto dhView = dh.getBlock(0, size);
 
             TestType dh2 = dh;
@@ -70,7 +66,7 @@ TEMPLATE_PRODUCT_TEST_CASE("Scenario: Constructing DataHandler", "", (DataHandle
                 {
                     dh2[0] = data_t(1);
                     REQUIRE(dh2[0] == data_t(1));
-                    REQUIRE((*dhView)[0] == data_t(0));
+                    REQUIRE((*dhView)[0] == randVec[0]);
                     REQUIRE((*dh2View)[0] == data_t(1));
                 }
             }

elsa/operators/FiniteDifferences.cpp

@@ -61,6 +61,8 @@ namespace elsa
     {
         Timer<> timeguard("FiniteDifferences", "apply");
 
+        Ax = 0;
+
         switch (_type) {
             case DiffType::FORWARD:
                 applyHelper(x, Ax, DiffType::FORWARD);
@@ -82,6 +84,8 @@ namespace elsa
     {
         Timer<> timeguard("FiniteDifferences", "applyAdjoint");
 
+        Aty = 0;
+
         switch (_type) {
             case DiffType::FORWARD:
                 applyAdjointHelper(y, Aty, DiffType::FORWARD);

elsa/problems/Problem.cpp

@@ -36,6 +36,8 @@ namespace elsa
                 throw std::invalid_argument(
                     "Problem: one of the reg terms' domain does not match the data term's");
         }
+
+        _currentSolution = 0;
     }
 
     template <typename data_t>
@@ -67,6 +69,8 @@ namespace elsa
             != regTerm.getFunctional().getDomainDescriptor().getNumberOfCoefficients())
             throw std::invalid_argument(
                 "Problem: one of the reg terms' domain does not match the data term's");
+
+        _currentSolution = 0;
     }
 
     template <typename data_t>
@@ -83,6 +87,7 @@ namespace elsa
     Problem<data_t>::Problem(const Functional<data_t>& dataTerm)
         : _dataTerm{dataTerm.clone()}, _currentSolution{dataTerm.getDomainDescriptor()}
     {
+        _currentSolution = 0;
     }
 
     template <typename data_t>

elsa/problems/WLSProblem.cpp

@@ -93,9 +93,10 @@ namespace elsa
             rangeDescList.push_back(regTermResidual->getRangeDescriptor().clone());
         }
 
-        // problem is convertible, allocate memory and start building block op
+        // problem is convertible, allocate memory, set to zero and start building block op
         RandomBlocksDescriptor dataVecDesc{rangeDescList};
         DataContainer<data_t> dataVec{dataVecDesc};
+        dataVec = 0;
         std::vector<std::unique_ptr<LinearOperator<data_t>>> opList(0);
 
         // add block corresponding to data term

elsa/problems/tests/test_Problem.cpp

@@ -51,6 +51,7 @@ SCENARIO("Testing Problem without regularization")
             THEN("the problem behaves as expected")
             {
                 DataContainer dcZero(dd);
+                dcZero = 0;
                 REQUIRE(prob.getCurrentSolution() == dcZero);
 
                 REQUIRE(prob.evaluate() == Approx(0.5 * dataVec.squaredNorm()));
@@ -140,6 +141,7 @@ SCENARIO("Testing Problem with one regularization term")
             THEN("the problem behaves as expected")
             {
                 DataContainer dcZero(dd);
+                dcZero = 0;
                 REQUIRE(prob.getCurrentSolution() == dcZero);
 
                 REQUIRE(prob.evaluate() == Approx(0.5 * dataVec.squaredNorm()));
@@ -236,6 +238,7 @@ SCENARIO("Testing Problem with several regularization terms")
             THEN("the problem behaves as expected")
             {
                 DataContainer dcZero(dd);
+                dcZero = 0;
                 REQUIRE(prob.getCurrentSolution() == dcZero);
 
                 REQUIRE(prob.evaluate() == Approx(0.5 * dataVec.squaredNorm()));

elsa/problems/tests/test_QuadricProblem.cpp

@@ -59,6 +59,7 @@ TEMPLATE_TEST_CASE("Scenario: Testing QuadricProblem", "", QuadricProblem<float>
             THEN("the problem behaves as expected")
             {
                 DataContainer<data_t> dcZero(dd);
+                dcZero = 0;
                 REQUIRE(prob.getCurrentSolution() == dcZero);
 
                 REQUIRE(prob.evaluate() == 0);
@@ -128,6 +129,7 @@ TEMPLATE_TEST_CASE("Scenario: Testing QuadricProblem", "", QuadricProblem<float>
             THEN("the problem behaves as expected")
             {
                 DataContainer<data_t> dcZero(dd);
+                dcZero = 0;
                 REQUIRE(prob.getCurrentSolution() == dcZero);
 
                 REQUIRE(prob.evaluate() == 0);
@@ -197,6 +199,7 @@ TEMPLATE_TEST_CASE("Scenario: Testing QuadricProblem", "", QuadricProblem<float>
             THEN("the problem behaves as expected")
             {
                 DataContainer<data_t> dcZero(dd);
+                dcZero = 0;
                 REQUIRE(prob.getCurrentSolution() == dcZero);
 
                 REQUIRE(prob.evaluate() == 0);
@@ -433,13 +436,14 @@ TEMPLATE_TEST_CASE("Scenario: Testing QuadricProblem", "", QuadricProblem<float>
 
                     AND_THEN("the problem behaves as expected")
                     {
-                        REQUIRE(prob.getCurrentSolution() == DataContainer<data_t>{dd});
+                        DataContainer<data_t> zero{dd};
+                        zero = 0;
+                        REQUIRE(prob.getCurrentSolution() == zero);
 
                         REQUIRE(prob.evaluate() == 0);
 
-                        DataContainer<data_t> grad = res.hasDataVector()
-                                                         ? static_cast<data_t>(-1.0) * dc
-                                                         : DataContainer<data_t>{dd};
+                        DataContainer<data_t> grad =
+                            res.hasDataVector() ? static_cast<data_t>(-1.0) * dc : zero;
                         if (res.hasDataVector() && res.hasOperator())
                             grad *= scaleFactor;
                         if (res.hasDataVector() && isWeighted)
@@ -594,7 +598,9 @@ TEMPLATE_TEST_CASE("Scenario: Testing QuadricProblem", "", QuadricProblem<float>
 
                         AND_THEN("the problem behaves as expected")
                         {
-                            REQUIRE(prob.getCurrentSolution() == DataContainer<data_t>{dd});
+                            DataContainer<data_t> zero{dd};
+                            zero = 0;
+                            REQUIRE(prob.getCurrentSolution() == zero);
 
                             REQUIRE(prob.evaluate() == 0);
 

elsa/problems/tests/test_TikhonovProblem.cpp

@@ -79,6 +79,7 @@ TEMPLATE_TEST_CASE("Scenario: Testing TikhonovProblem with one regularization te
             THEN("the TikhonovProblem behaves as expected")
             {
                 DataContainer<TestType> dcZero(dd);
+                dcZero = 0;
                 REQUIRE(prob.getCurrentSolution() == dcZero);
 
                 REQUIRE(prob.evaluate() == Approx(0.5 * dataVec.squaredNorm()));
@@ -205,6 +206,7 @@ TEMPLATE_TEST_CASE("Scenario: Testing TikhonovProblem with several regularizatio
             THEN("the TikhonovProblem behaves as expected")
             {
                 DataContainer<TestType> dcZero(dd);
+                dcZero = 0;
                 REQUIRE(prob.getCurrentSolution() == dcZero);
 
                 REQUIRE(prob.evaluate() == Approx(0.5 * dataVec.squaredNorm()));

elsa/problems/tests/test_WLSProblem.cpp

@@ -52,6 +52,7 @@ TEMPLATE_TEST_CASE("Scenario: Testing WLSProblem", "", float, double)
             THEN("the problem behaves as expected")
             {
                 DataContainer<TestType> dcZero(dd);
+                dcZero = 0;
                 REQUIRE(prob.getCurrentSolution() == dcZero);
 
                 REQUIRE(prob.evaluate() == Approx(0.5 * bVec.squaredNorm()));
@@ -123,6 +124,7 @@ TEMPLATE_TEST_CASE("Scenario: Testing WLSProblem", "", float, double)
             THEN("the problem behaves as expected")
             {
                 DataContainer<TestType> dcZero(dd);
+                dcZero = 0;
                 REQUIRE(prob.getCurrentSolution() == dcZero);
 
                 REQUIRE(prob.evaluate()
@@ -153,6 +155,7 @@ TEMPLATE_TEST_CASE("Scenario: Testing WLSProblem", "", float, double)
             THEN("the problem behaves as expected")
             {
                 DataContainer<TestType> dcZero(dd);
+                dcZero = 0;
                 REQUIRE(prob.getCurrentSolution() == dcX0);
 
                 REQUIRE(

elsa/projectors/tests/test_BinaryMethod.cpp

@@ -500,6 +500,7 @@ SCENARIO("Output DataContainer is not zero initialized")
                 op.apply(volume, sino);
 
                 DataContainer zero(sinoDescriptor);
+                zero = 0;
                 REQUIRE(sino == zero);
             }
         }
@@ -514,6 +515,7 @@ SCENARIO("Output DataContainer is not zero initialized")
                 op.applyAdjoint(sino, volume);
 
                 DataContainer zero(volumeDescriptor);
+                zero = 0;
                 REQUIRE(volume == zero);
             }
         }
@@ -548,6 +550,7 @@ SCENARIO("Output DataContainer is not zero initialized")
                 op.apply(volume, sino);
 
                 DataContainer zero(sinoDescriptor);
+                zero = 0;
                 REQUIRE(sino == zero);
             }
         }
@@ -562,6 +565,7 @@ SCENARIO("Output DataContainer is not zero initialized")
                 op.applyAdjoint(sino, volume);
 
                 DataContainer zero(volumeDescriptor);
+                zero = 0;
                 REQUIRE(volume == zero);
             }
         }
@@ -600,6 +604,7 @@ SCENARIO("Rays not intersecting the bounding box are present")
                 op.apply(volume, sino);
 
                 DataContainer zeroSino(sinoDescriptor);
+                zeroSino = 0;
                 REQUIRE(sino == zeroSino);
 
                 AND_THEN("Result of backprojection is zero")
@@ -607,6 +612,7 @@ SCENARIO("Rays not intersecting the bounding box are present")
                     op.applyAdjoint(sino, volume);
 
                     DataContainer zeroVolume(volumeDescriptor);
+                    zeroVolume = 0;
                     REQUIRE(volume == zeroVolume);
                 }
             }
@@ -625,6 +631,7 @@ SCENARIO("Rays not intersecting the bounding box are present")
                 op.apply(volume, sino);
 
                 DataContainer zeroSino(sinoDescriptor);
+                zeroSino = 0;
                 REQUIRE(sino == zeroSino);
 
                 AND_THEN("Result of backprojection is zero")
@@ -632,6 +639,7 @@ SCENARIO("Rays not intersecting the bounding box are present")
                     op.applyAdjoint(sino, volume);
 
                     DataContainer zeroVolume(volumeDescriptor);
+                    zeroVolume = 0;
                     REQUIRE(volume == zeroVolume);
                 }
             }
@@ -649,6 +657,7 @@ SCENARIO("Rays not intersecting the bounding box are present")
                 op.apply(volume, sino);
 
                 DataContainer zeroSino(sinoDescriptor);
+                zeroSino = 0;
                 REQUIRE(sino == zeroSino);
 
                 AND_THEN("Result of backprojection is zero")
@@ -656,6 +665,7 @@ SCENARIO("Rays not intersecting the bounding box are present")
                     op.applyAdjoint(sino, volume);
 
                     DataContainer zeroVolume(volumeDescriptor);
+                    zeroVolume = 0;
                     REQUIRE(volume == zeroVolume);
                 }
             }
@@ -674,6 +684,7 @@ SCENARIO("Rays not intersecting the bounding box are present")
                 op.apply(volume, sino);
 
                 DataContainer zeroSino(sinoDescriptor);
+                zeroSino = 0;
                 REQUIRE(sino == zeroSino);
 
                 AND_THEN("Result of backprojection is zero")
@@ -681,6 +692,7 @@ SCENARIO("Rays not intersecting the bounding box are present")
                     op.applyAdjoint(sino, volume);
 
                     DataContainer zeroVolume(volumeDescriptor);
+                    zeroVolume = 0;
                     REQUIRE(volume == zeroVolume);
                 }
             }
@@ -731,6 +743,7 @@ SCENARIO("Rays not intersecting the bounding box are present")
                     op.apply(volume, sino);
 
                     DataContainer zeroSino(sinoDescriptor);
+                    zeroSino = 0;
                     REQUIRE(sino == zeroSino);
 
                     AND_THEN("Result of backprojection is zero")
@@ -738,6 +751,7 @@ SCENARIO("Rays not intersecting the bounding box are present")
                         op.applyAdjoint(sino, volume);
 
                         DataContainer zeroVolume(volumeDescriptor);
+                        zeroVolume = 0;
                         REQUIRE(volume == zeroVolume);
                     }
                 }
@@ -860,6 +874,7 @@ SCENARIO("Axis-aligned rays are present")
                     op.applyAdjoint(sino, volume);
 
                     DataContainer zero(volumeDescriptor);
+                    zero = 0;
                     REQUIRE(volume == zero);
                 }
             }
@@ -1078,8 +1093,9 @@ SCENARIO("Axis-aligned rays are present")
                     {
                         sino[0] = 1;
                         op.applyAdjoint(sino, volume);
-
-                        REQUIRE(volume == DataContainer(volumeDescriptor));
+                        DataContainer zero(volumeDescriptor);
+                        zero = 0;
+                        REQUIRE(volume == zero);
                     }
                 }
             }
@@ -1260,6 +1276,7 @@ SCENARIO("Projection under an angle")
                 op.apply(volume, sino);
 
                 DataContainer sZero(sinoDescriptor);
+                sZero = 0;
                 REQUIRE(sino == sZero);
 
                 AND_THEN("The correct weighting is applied")
@@ -1322,6 +1339,7 @@ SCENARIO("Projection under an angle")
 
                 op.apply(volume, sino);
                 DataContainer sZero(sinoDescriptor);
+                sZero = 0;
                 REQUIRE(sino == sZero);
 
                 AND_THEN("The correct weighting is applied")
@@ -1373,6 +1391,7 @@ SCENARIO("Projection under an angle")
 
                 op.apply(volume, sino);
                 DataContainer sZero(sinoDescriptor);
+                sZero = 0;
                 REQUIRE(sino == sZero);
 
                 AND_THEN("The correct weighting is applied")
@@ -1416,6 +1435,7 @@ SCENARIO("Projection under an angle")
 
                 op.apply(volume, sino);
                 DataContainer sZero(sinoDescriptor);
+                sZero = 0;
                 REQUIRE(sino == sZero);
 
                 AND_THEN("The correct weighting is applied")
@@ -1469,6 +1489,7 @@ SCENARIO("Projection under an angle")
 
                 op.apply(volume, sino);
                 DataContainer sZero(sinoDescriptor);
+                sZero = 0;
                 REQUIRE(sino == sZero);
 
                 AND_THEN("The correct weighting is applied")
@@ -1531,6 +1552,7 @@ SCENARIO("Projection under an angle")
 
                 op.apply(volume, sino);
                 DataContainer sZero(sinoDescriptor);
+                sZero = 0;
                 REQUIRE(sino == sZero);
 
                 AND_THEN("The correct weighting is applied")
@@ -1583,6 +1605,7 @@ SCENARIO("Projection under an angle")
 
                 op.apply(volume, sino);
                 DataContainer sZero(sinoDescriptor);
+                sZero = 0;
                 REQUIRE(sino == sZero);
 
                 AND_THEN("The correct weighting is applied")
@@ -1628,6 +1651,7 @@ SCENARIO("Projection under an angle")
 
                 op.apply(volume, sino);
                 DataContainer sZero(sinoDescriptor);
+                sZero = 0;
                 REQUIRE(sino == sZero);
 
                 AND_THEN("The correct weighting is applied")

elsa/projectors/tests/test_JosephsMethod.cpp

@@ -229,7 +229,9 @@ SCENARIO("Calls to functions of super class")
         DataDescriptor volumeDescriptor(volumeDims);
         DataDescriptor sinoDescriptor(sinoDims);
         DataContainer volume(volumeDescriptor);
+        volume = 0;
         DataContainer sino(sinoDescriptor);
+        sino = 0;
         std::vector<Geometry> geom;
         for (std::size_t i = 0; i < numImgs; i++) {
             real_t angle = static_cast<real_t>(i * 2) * pi_t / 50;
@@ -291,6 +293,7 @@ SCENARIO("Output DataContainer is not zero initialized")
                 op.apply(volume, sino);
 
                 DataContainer zero(sinoDescriptor);
+                zero = 0;
                 REQUIRE(sino == zero);
             }
         }
@@ -305,6 +308,7 @@ SCENARIO("Output DataContainer is not zero initialized")
                 op.applyAdjoint(sino, volume);
 
                 DataContainer zero(volumeDescriptor);
+                zero = 0;
                 REQUIRE(volume == zero);
             }
         }
@@ -337,6 +341,7 @@ SCENARIO("Output DataContainer is not zero initialized")
             {
                 op.apply(volume, sino);
                 DataContainer zero(sinoDescriptor);
+                zero = 0;
                 REQUIRE(sino == zero);
             }
         }
@@ -350,6 +355,7 @@ SCENARIO("Output DataContainer is not zero initialized")
             {
                 op.applyAdjoint(sino, volume);
                 DataContainer zero(volumeDescriptor);
+                zero = 0;
                 REQUIRE(volume == zero);
             }
         }
@@ -386,12 +392,14 @@ SCENARIO("Rays not intersecting the bounding box are present")
             {
                 op.apply(volume, sino);
                 DataContainer zero(sinoDescriptor);
+                zero = 0;
                 REQUIRE(sino == zero);
 
                 AND_THEN("Result of backprojection is zero")
                 {
                     op.applyAdjoint(sino, volume);
                     DataContainer zero(volumeDescriptor);
+                    zero = 0;
                     REQUIRE(volume == zero);
                 }
             }
@@ -410,12 +418,14 @@ SCENARIO("Rays not intersecting the bounding box are present")
             {
                 op.apply(volume, sino);
                 DataContainer zero(sinoDescriptor);
+                zero = 0;
                 REQUIRE(sino == zero);
 
                 AND_THEN("Result of backprojection is zero")
                 {
                     op.applyAdjoint(sino, volume);
                     DataContainer zero(volumeDescriptor);
+                    zero = 0;
                     REQUIRE(volume == zero);
                 }
             }
@@ -433,12 +443,14 @@ SCENARIO("Rays not intersecting the bounding box are present")
             {
                 op.apply(volume, sino);
                 DataContainer zero(sinoDescriptor);
+                zero = 0;
                 REQUIRE(sino == zero);
 
                 AND_THEN("Result of backprojection is zero")
                 {
                     op.applyAdjoint(sino, volume);
                     DataContainer zero(volumeDescriptor);
+                    zero = 0;
                     REQUIRE(volume == zero);
                 }
             }
@@ -457,12 +469,14 @@ SCENARIO("Rays not intersecting the bounding box are present")
             {
                 op.apply(volume, sino);
                 DataContainer zero(sinoDescriptor);
+                zero = 0;
                 REQUIRE(sino == zero);
 
                 AND_THEN("Result of backprojection is zero")
                 {
                     op.applyAdjoint(sino, volume);
                     DataContainer zero(volumeDescriptor);
+                    zero = 0;
                     REQUIRE(volume == zero);