Update interpolation input and output handling

mrmd/assert/verbose.hpp

@@ -92,12 +92,6 @@ std::string generateAssertionMessage(const T& lhs,
         assumption::log(ss.str()); \
     }
 
-template <typename FLOAT>
-bool isFloatEqual(const FLOAT& lhs, const FLOAT& rhs)
-{
-    return std::abs(lhs - rhs) < FLOAT(1e-10);
-}
-
 // macro overloading (-> https://stackoverflow.com/a/24028231)
 
 #define GLUE(x, y) x y

mrmd/datatypes.hpp

@@ -125,6 +125,15 @@ KOKKOS_INLINE_FUNCTION constexpr float float_c(const T& val)
     return static_cast<float>(val);
 }
 
+template <typename FLOAT>
+KOKKOS_INLINE_FUNCTION bool isFloatEqual(
+    const FLOAT& lhs,
+    const FLOAT& rhs,
+    const FLOAT& epsilon = std::numeric_limits<FLOAT>::epsilon())
+{
+    return std::abs(lhs - rhs) < epsilon;
+}
+
 enum class AXIS : idx_t
 {
     X = 0,  ///< index of the x coordinate in vector views

mrmd/util/interpolation.cpp

@@ -18,49 +18,33 @@ namespace mrmd
 {
 namespace util
 {
-data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid)
+void updateInterpolate(const data::MultiHistogram& target, const data::MultiHistogram& input)
 {
-    real_t gridSpacing = grid(1) - grid(0);
-    real_t gridMin = grid(0) - 0.5_r * gridSpacing;
-    real_t gridMax = grid(grid.extent(0) - 1) + 0.5_r * gridSpacing;
+    MRMD_HOST_ASSERT_EQUAL(target.numHistograms, input.numHistograms);
 
-    data::MultiHistogram output(
-        "interpolated-profile", gridMin, gridMax, idx_c(grid.extent(0)), input.numHistograms);
-
-    MRMD_HOST_ASSERT_EQUAL(output.numBins, idx_c(grid.extent(0)), "Output grid size mismatch!");
-    for (idx_t idx = 0; idx < idx_c(output.numBins); ++idx)
-    {
-        MRMD_HOST_ASSERT_EQUAL(output.getBinPosition(idx), grid(idx), "Output grid mismatch!");
-    }
-
-    auto policy = Kokkos::MDRangePolicy<Kokkos::Rank<2>>(
-        {idx_t(0), idx_t(0)}, {idx_c(grid.extent(0)), input.numHistograms});
+    auto policy = Kokkos::MDRangePolicy<Kokkos::Rank<2>>({idx_t(0), idx_t(0)},
+                                                         {target.numBins, input.numHistograms});
     auto kernel = KOKKOS_LAMBDA(const idx_t binIdx, const idx_t histogramIdx)
     {
         // find the two enclosing bins in the input histogram
-        real_t outputBinPosition = grid(binIdx);
+        real_t outputBinPosition = target.getBinPosition(binIdx);
         idx_t leftBinIdx = input.getBin(outputBinPosition - 0.5_r * input.binSize);
         idx_t rightBinIdx = leftBinIdx + 1;
 
-        // handle boundaries
-        if (leftBinIdx < 0 || rightBinIdx >= input.numBins)
+        // only update if within bounds of input histogram
+        if (leftBinIdx >= 0 && rightBinIdx < input.numBins)
         {
-            output.data(binIdx, histogramIdx) = 0_r;  // out of bounds, set to zero
-            return;
-        }
+            auto inputLeft = input.data(leftBinIdx, histogramIdx);
+            auto inputRight = input.data(rightBinIdx, histogramIdx);
 
-        auto inputDataLeft = input.data(leftBinIdx, histogramIdx);
-        auto inputDataRight = input.data(rightBinIdx, histogramIdx);
-
-        output.data(binIdx, histogramIdx) =
-            lerp(inputDataLeft,
-                 inputDataRight,
-                 (outputBinPosition - input.getBinPosition(leftBinIdx)) * input.inverseBinSize);
+            target.data(binIdx, histogramIdx) +=
+                lerp(inputLeft,
+                     inputRight,
+                     (outputBinPosition - input.getBinPosition(leftBinIdx)) * input.inverseBinSize);
+        }
     };
-    Kokkos::parallel_for("MultiHistogram::interpolate", policy, kernel);
+    Kokkos::parallel_for("MultiHistogram::updateInterpolate", policy, kernel);
     Kokkos::fence();
-
-    return output;
 }
 
 }  // namespace util

mrmd/util/interpolation.hpp

@@ -14,6 +14,8 @@
 
 #pragma once
 
+#include <limits>
+
 #include "data/MultiHistogram.hpp"
 #include "util/IsInSymmetricSlab.hpp"
 
@@ -35,13 +37,14 @@ real_t lerp(const real_t& left, const real_t& right, const real_t& factor)
 }
 
 /**
- * Linear interpolation of data contained in input MultiHistogram onto given grid.
- * Data for grid points outside of the grid range of the input MultiHistogram are set to zero.
- * @param input input MultiHistogram containing data to interpolate
- * @param grid grid to interpolate data onto
- * @return MultiHistogram containing interpolated data on given grid
+ * Linear interpolation of data contained in input MultiHistogram onto grid of target histogram,
+ * adding the interpolated data to the data already contained in target. Data for target grid points
+ * outside of the grid range of the input MultiHistogram is not updated.
+ * @param input MultiHistogram containing data to interpolate on coarse grid.
+ * @param target MultiHistogram defining the grid to interpolate onto and containing the data
+ * to be updated by interpolating the data from input.
  */
-data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid);
+void updateInterpolate(const data::MultiHistogram& target, const data::MultiHistogram& input);
 
 }  // namespace util
 }  // namespace mrmd

mrmd/util/interpolation.test.cpp

@@ -22,38 +22,101 @@ namespace mrmd
 {
 namespace util
 {
-TEST(interpolate, testInterpolate)
+TEST(interpolate, preZeroInner)
 {
-    data::MultiHistogram histogramCoarse("histogram", 0_r, 10_r, 10, 3);
-    data::MultiHistogram histogramFine("histogram", 0.5_r, 9.5_r, 30, 3);
+    data::MultiHistogram histogramInput("histogramInput", 0_r, 10_r, 10, 3);
+    data::MultiHistogram histogramTarget("histogramTarget", 0.5_r, 9.5_r, 30, 3);
+    data::MultiHistogram histogramRef("histogramRef", histogramTarget);
 
-    auto h_dataCoarse = Kokkos::create_mirror_view(histogramCoarse.data);
+    auto h_dataInput = Kokkos::create_mirror_view(histogramInput.data);
     for (auto idx = 0; idx < 10; ++idx)
     {
-        h_dataCoarse(idx, 0) = 0_r;
-        h_dataCoarse(idx, 1) = 1_r;
-        h_dataCoarse(idx, 2) = histogramCoarse.getBinPosition(idx);
+        h_dataInput(idx, 0) = 0_r;
+        h_dataInput(idx, 1) = 1_r;
+        h_dataInput(idx, 2) = histogramInput.getBinPosition(idx);
     }
-    Kokkos::deep_copy(histogramCoarse.data, h_dataCoarse);
+    Kokkos::deep_copy(histogramInput.data, h_dataInput);
 
-    auto h_dataFine = Kokkos::create_mirror_view(histogramFine.data);
+    auto h_dataTarget = Kokkos::create_mirror_view(histogramTarget.data);
     for (auto idx = 0; idx < 30; ++idx)
     {
-        h_dataFine(idx, 0) = 0_r;
-        h_dataFine(idx, 1) = 1_r;
-        h_dataFine(idx, 2) = histogramFine.getBinPosition(idx);
+        h_dataTarget(idx, 0) = 0_r;
+        h_dataTarget(idx, 1) = 0_r;
+        h_dataTarget(idx, 2) = 0_r;
     }
-    Kokkos::deep_copy(histogramFine.data, h_dataFine);
+    Kokkos::deep_copy(histogramTarget.data, h_dataTarget);
 
-    auto histogramInterp = util::interpolate(histogramCoarse, createGrid(histogramFine));
-    auto h_dataInterp =
-        Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramInterp.data);
+    auto h_dataRef = Kokkos::create_mirror_view(histogramRef.data);
+    for (auto idx = 0; idx < 30; ++idx)
+    {
+        h_dataRef(idx, 0) = 0_r;
+        h_dataRef(idx, 1) = 1_r;
+        h_dataRef(idx, 2) = histogramRef.getBinPosition(idx);
+    }
+    Kokkos::deep_copy(histogramRef.data, h_dataRef);
+
+    util::updateInterpolate(histogramTarget, histogramInput);
+    h_dataTarget = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramTarget.data);
+
+    for (auto idx = 0; idx < 30; ++idx)
+    {
+        EXPECT_FLOAT_EQ(h_dataTarget(idx, 0), h_dataRef(idx, 0));
+        EXPECT_FLOAT_EQ(h_dataTarget(idx, 1), h_dataRef(idx, 1));
+        EXPECT_FLOAT_EQ(h_dataTarget(idx, 2), h_dataRef(idx, 2));
+    }
+}
+
+TEST(interpolate, nonZeroWithBoundary)
+{
+    data::MultiHistogram histogramInput("histogramInput", 0_r, 10_r, 10, 3);
+    data::MultiHistogram histogramTarget("histogramTarget", 0_r, 10_r, 30, 3);
+    data::MultiHistogram histogramRef("histogramRef", histogramTarget);
+
+    auto h_dataInput = Kokkos::create_mirror_view(histogramInput.data);
+    for (auto idx = 0; idx < 10; ++idx)
+    {
+        h_dataInput(idx, 0) = 0_r;
+        h_dataInput(idx, 1) = 1_r;
+        h_dataInput(idx, 2) = histogramInput.getBinPosition(idx);
+    }
+    Kokkos::deep_copy(histogramInput.data, h_dataInput);
+
+    auto h_dataTarget = Kokkos::create_mirror_view(histogramTarget.data);
+    for (auto idx = 0; idx < 30; ++idx)
+    {
+        h_dataTarget(idx, 0) = 1_r;
+        h_dataTarget(idx, 1) = 1_r;
+        h_dataTarget(idx, 2) = 1_r;
+    }
+    Kokkos::deep_copy(histogramTarget.data, h_dataTarget);
+
+    auto h_dataRef = Kokkos::create_mirror_view(histogramRef.data);
+    for (auto idx = 0; idx < 30; ++idx)
+    {
+        if (histogramRef.getBinPosition(idx) >= histogramInput.getBinPosition(0) &&
+            histogramRef.getBinPosition(idx) < histogramInput.getBinPosition(9))
+        {
+            h_dataRef(idx, 0) = 1_r + 0_r;
+            h_dataRef(idx, 1) = 1_r + 1_r;
+            h_dataRef(idx, 2) = 1_r + histogramRef.getBinPosition(idx);
+        }
+        else
+        {
+            h_dataRef(idx, 0) = 1_r + 0_r;
+            h_dataRef(idx, 1) = 1_r + 0_r;
+            h_dataRef(idx, 2) = 1_r + 0_r;
+        }
+    }
+    Kokkos::deep_copy(histogramRef.data, h_dataRef);
+
+    util::updateInterpolate(histogramTarget, histogramInput);
+    h_dataTarget = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramTarget.data);
 
     for (auto idx = 0; idx < 30; ++idx)
     {
-        EXPECT_FLOAT_EQ(h_dataFine(idx, 0), h_dataInterp(idx, 0));
-        EXPECT_FLOAT_EQ(h_dataFine(idx, 1), h_dataInterp(idx, 1));
-        EXPECT_FLOAT_EQ(h_dataFine(idx, 2), h_dataInterp(idx, 2));
+        EXPECT_FLOAT_EQ(h_dataTarget(idx, 0), h_dataRef(idx, 0));
+        EXPECT_FLOAT_EQ(h_dataTarget(idx, 1), h_dataRef(idx, 1));
+        EXPECT_FLOAT_EQ(h_dataTarget(idx, 2), h_dataRef(idx, 2));
     }
 }