Add interpolation

mrmd/util/CMakeLists.txt

@@ -24,5 +24,8 @@ mrmd_add_test(NAME test.util.IsInSymmetricSlab
 mrmd_add_test(NAME test.util.math
         FILES math.test.cpp)
 
+mrmd_add_test(NAME test.util.interpolation
+        FILES interpolation.test.cpp)
+
 mrmd_add_test(NAME test.util.Random
         FILES Random.test.cpp)

mrmd/util/interpolation.cpp

@@ -0,0 +1,67 @@
+// Copyright 2024 Sebastian Eibl
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "util/interpolation.hpp"
+
+namespace mrmd
+{
+namespace util
+{
+data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid)
+{
+    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;
+
+    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 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);
+        idx_t leftBinIdx = input.getBin(outputBinPosition - 0.5_r * input.binSize);
+        idx_t rightBinIdx = leftBinIdx + 1;
+
+        // handle boundaries
+        if (leftBinIdx < 0 || rightBinIdx >= input.numBins)
+        {
+            output.data(binIdx, histogramIdx) = 0_r;  // out of bounds, set to zero
+            return;
+        }
+
+        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);
+    };
+    Kokkos::parallel_for("MultiHistogram::interpolate", policy, kernel);
+    Kokkos::fence();
+
+    return output;
+}
+
+}  // namespace util
+}  // namespace mrmd

mrmd/util/interpolation.hpp

@@ -0,0 +1,47 @@
+// Copyright 2024 Sebastian Eibl
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#pragma once
+
+#include "data/MultiHistogram.hpp"
+#include "util/IsInSymmetricSlab.hpp"
+
+namespace mrmd
+{
+namespace util
+{
+/**
+ * Linear interpolation between two values.
+ * @param left left value
+ * @param right right value
+ * @param factor interpolation factor in [0, 1]
+ * @return interpolated value
+ */
+KOKKOS_INLINE_FUNCTION
+real_t lerp(const real_t& left, const real_t& right, const real_t& factor)
+{
+    return left + (right - left) * 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
+ */
+data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid);
+
+}  // namespace util
+}  // namespace mrmd

mrmd/util/interpolation.test.cpp

@@ -0,0 +1,61 @@
+// Copyright 2024 Sebastian Eibl
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "util/interpolation.hpp"
+
+#include <gtest/gtest.h>
+
+#include "util/IsInSymmetricSlab.hpp"
+
+namespace mrmd
+{
+namespace util
+{
+TEST(interpolate, testInterpolate)
+{
+    data::MultiHistogram histogramCoarse("histogram", 0_r, 10_r, 10, 3);
+    data::MultiHistogram histogramFine("histogram", 0.5_r, 9.5_r, 30, 3);
+
+    auto h_dataCoarse = Kokkos::create_mirror_view(histogramCoarse.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);
+    }
+    Kokkos::deep_copy(histogramCoarse.data, h_dataCoarse);
+
+    auto h_dataFine = Kokkos::create_mirror_view(histogramFine.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);
+    }
+    Kokkos::deep_copy(histogramFine.data, h_dataFine);
+
+    auto histogramInterp = util::interpolate(histogramCoarse, createGrid(histogramFine));
+    auto h_dataInterp =
+        Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramInterp.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));
+    }
+}
+
+}  // namespace util
+}  // namespace mrmd