From 9bfc1ba0b618b8e6c3ab2c32be62738a72f87444 Mon Sep 17 00:00:00 2001 From: julianhille Date: Mon, 26 Jan 2026 19:17:57 +0100 Subject: [PATCH 1/8] changed input for interpolation, Ref. #047 --- mrmd/util/interpolation.cpp | 29 ++++++++++------------------- mrmd/util/interpolation.hpp | 6 +++--- mrmd/util/interpolation.test.cpp | 2 +- 3 files changed, 14 insertions(+), 23 deletions(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index 59ce613..7de80f2 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -18,44 +18,35 @@ namespace mrmd { namespace util { -data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid) +data::MultiHistogram interpolate(const data::MultiHistogram& inputCoarse, const data::MultiHistogram& inputFine) { - 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(inputFine.numHistograms, inputCoarse.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!"); - } + data::MultiHistogram output("interpolated-profile", inputFine); auto policy = Kokkos::MDRangePolicy>( - {idx_t(0), idx_t(0)}, {idx_c(grid.extent(0)), input.numHistograms}); + {idx_t(0), idx_t(0)}, {output.numBins, inputCoarse.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); + real_t outputBinPosition = output.getBinPosition(binIdx); + idx_t leftBinIdx = inputCoarse.getBin(outputBinPosition - 0.5_r * inputCoarse.binSize); idx_t rightBinIdx = leftBinIdx + 1; // handle boundaries - if (leftBinIdx < 0 || rightBinIdx >= input.numBins) + if (leftBinIdx < 0 || rightBinIdx >= inputCoarse.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); + auto inputDataLeft = inputCoarse.data(leftBinIdx, histogramIdx); + auto inputDataRight = inputCoarse.data(rightBinIdx, histogramIdx); output.data(binIdx, histogramIdx) = lerp(inputDataLeft, inputDataRight, - (outputBinPosition - input.getBinPosition(leftBinIdx)) * input.inverseBinSize); + (outputBinPosition - inputCoarse.getBinPosition(leftBinIdx)) * inputCoarse.inverseBinSize); }; Kokkos::parallel_for("MultiHistogram::interpolate", policy, kernel); Kokkos::fence(); diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp index 41987b2..3b77fd5 100644 --- a/mrmd/util/interpolation.hpp +++ b/mrmd/util/interpolation.hpp @@ -37,11 +37,11 @@ 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 + * @param inputCoarse input MultiHistogram containing data to interpolate on coarse grid. + * @param inputFine MultiHistogram defining the grid to interpolate onto. * @return MultiHistogram containing interpolated data on given grid */ -data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid); +data::MultiHistogram interpolate(const data::MultiHistogram& inputCoarse, const data::MultiHistogram& inputFine); } // namespace util } // namespace mrmd \ No newline at end of file diff --git a/mrmd/util/interpolation.test.cpp b/mrmd/util/interpolation.test.cpp index e136467..7c5cdd3 100644 --- a/mrmd/util/interpolation.test.cpp +++ b/mrmd/util/interpolation.test.cpp @@ -45,7 +45,7 @@ TEST(interpolate, testInterpolate) } Kokkos::deep_copy(histogramFine.data, h_dataFine); - auto histogramInterp = util::interpolate(histogramCoarse, createGrid(histogramFine)); + auto histogramInterp = util::interpolate(histogramCoarse, histogramFine); auto h_dataInterp = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramInterp.data); From fbec185004ecc67c4e798697c1de44ef2cfbfa17 Mon Sep 17 00:00:00 2001 From: julianhille Date: Mon, 26 Jan 2026 19:19:06 +0100 Subject: [PATCH 2/8] clang-format fixes, Ref. #047 --- mrmd/util/interpolation.cpp | 6 ++++-- mrmd/util/interpolation.hpp | 3 ++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index 7de80f2..3f31d6a 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -18,7 +18,8 @@ namespace mrmd { namespace util { -data::MultiHistogram interpolate(const data::MultiHistogram& inputCoarse, const data::MultiHistogram& inputFine) +data::MultiHistogram interpolate(const data::MultiHistogram& inputCoarse, + const data::MultiHistogram& inputFine) { MRMD_HOST_ASSERT_EQUAL(inputFine.numHistograms, inputCoarse.numHistograms); @@ -46,7 +47,8 @@ data::MultiHistogram interpolate(const data::MultiHistogram& inputCoarse, const output.data(binIdx, histogramIdx) = lerp(inputDataLeft, inputDataRight, - (outputBinPosition - inputCoarse.getBinPosition(leftBinIdx)) * inputCoarse.inverseBinSize); + (outputBinPosition - inputCoarse.getBinPosition(leftBinIdx)) * + inputCoarse.inverseBinSize); }; Kokkos::parallel_for("MultiHistogram::interpolate", policy, kernel); Kokkos::fence(); diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp index 3b77fd5..ac7d9e5 100644 --- a/mrmd/util/interpolation.hpp +++ b/mrmd/util/interpolation.hpp @@ -41,7 +41,8 @@ real_t lerp(const real_t& left, const real_t& right, const real_t& factor) * @param inputFine MultiHistogram defining the grid to interpolate onto. * @return MultiHistogram containing interpolated data on given grid */ -data::MultiHistogram interpolate(const data::MultiHistogram& inputCoarse, const data::MultiHistogram& inputFine); +data::MultiHistogram interpolate(const data::MultiHistogram& inputCoarse, + const data::MultiHistogram& inputFine); } // namespace util } // namespace mrmd \ No newline at end of file From b8c110ac72860c26fa29dc9588fbc264038a336a Mon Sep 17 00:00:00 2001 From: julianhille Date: Tue, 27 Jan 2026 15:09:56 +0100 Subject: [PATCH 3/8] updated interpolation input and output handling, Ref. #047 --- mrmd/util/interpolation.cpp | 34 ++++++++++++++------------------ mrmd/util/interpolation.hpp | 14 ++++++------- mrmd/util/interpolation.test.cpp | 20 +++++++++---------- 3 files changed, 32 insertions(+), 36 deletions(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index 3f31d6a..206773e 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -18,42 +18,38 @@ namespace mrmd { namespace util { -data::MultiHistogram interpolate(const data::MultiHistogram& inputCoarse, - const data::MultiHistogram& inputFine) +void updateInterpolate(const data::MultiHistogram& inputData, data::MultiHistogram& inputTarget) { - MRMD_HOST_ASSERT_EQUAL(inputFine.numHistograms, inputCoarse.numHistograms); - - data::MultiHistogram output("interpolated-profile", inputFine); + MRMD_HOST_ASSERT_EQUAL(inputTarget.numHistograms, inputData.numHistograms); auto policy = Kokkos::MDRangePolicy>( - {idx_t(0), idx_t(0)}, {output.numBins, inputCoarse.numHistograms}); + {idx_t(0), idx_t(0)}, {inputTarget.numBins, inputData.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 = output.getBinPosition(binIdx); - idx_t leftBinIdx = inputCoarse.getBin(outputBinPosition - 0.5_r * inputCoarse.binSize); + real_t outputBinPosition = inputTarget.getBinPosition(binIdx); + idx_t leftBinIdx = inputData.getBin(outputBinPosition - 0.5_r * inputData.binSize); idx_t rightBinIdx = leftBinIdx + 1; // handle boundaries - if (leftBinIdx < 0 || rightBinIdx >= inputCoarse.numBins) + if (leftBinIdx < 0 || rightBinIdx >= inputData.numBins) { - output.data(binIdx, histogramIdx) = 0_r; // out of bounds, set to zero + inputTarget.data(binIdx, histogramIdx) += 0_r; // out of bounds, set to zero return; } - auto inputDataLeft = inputCoarse.data(leftBinIdx, histogramIdx); - auto inputDataRight = inputCoarse.data(rightBinIdx, histogramIdx); + auto inputDataLeft = inputData.data(leftBinIdx, histogramIdx); + auto inputDataRight = inputData.data(rightBinIdx, histogramIdx); - output.data(binIdx, histogramIdx) = - lerp(inputDataLeft, - inputDataRight, - (outputBinPosition - inputCoarse.getBinPosition(leftBinIdx)) * - inputCoarse.inverseBinSize); + inputTarget.data(binIdx, histogramIdx) += lerp( + inputDataLeft, + inputDataRight, + (outputBinPosition - inputData.getBinPosition(leftBinIdx)) * inputData.inverseBinSize); }; - Kokkos::parallel_for("MultiHistogram::interpolate", policy, kernel); + Kokkos::parallel_for("MultiHistogram::updateInterpolate", policy, kernel); Kokkos::fence(); - return output; + return; } } // namespace util diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp index ac7d9e5..5ec9152 100644 --- a/mrmd/util/interpolation.hpp +++ b/mrmd/util/interpolation.hpp @@ -35,14 +35,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 inputCoarse input MultiHistogram containing data to interpolate on coarse grid. - * @param inputFine MultiHistogram defining the grid to interpolate onto. - * @return MultiHistogram containing interpolated data on given grid + * Linear interpolation of data contained in input MultiHistogram onto grid of target histogram, + * updating the data of target. 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 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& inputCoarse, - const data::MultiHistogram& inputFine); +void updateInterpolate(const data::MultiHistogram& input, data::MultiHistogram& target); } // namespace util } // namespace mrmd \ No newline at end of file diff --git a/mrmd/util/interpolation.test.cpp b/mrmd/util/interpolation.test.cpp index 7c5cdd3..9689ab0 100644 --- a/mrmd/util/interpolation.test.cpp +++ b/mrmd/util/interpolation.test.cpp @@ -24,30 +24,30 @@ 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); + data::MultiHistogram histogramData("histogramData", 0_r, 10_r, 10, 3); + data::MultiHistogram histogramTarget("histogramTarget", 0.5_r, 9.5_r, 30, 3); - auto h_dataCoarse = Kokkos::create_mirror_view(histogramCoarse.data); + auto h_dataCoarse = Kokkos::create_mirror_view(histogramData.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_dataCoarse(idx, 2) = histogramData.getBinPosition(idx); } - Kokkos::deep_copy(histogramCoarse.data, h_dataCoarse); + Kokkos::deep_copy(histogramData.data, h_dataCoarse); - auto h_dataFine = Kokkos::create_mirror_view(histogramFine.data); + auto h_dataFine = 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_dataFine(idx, 2) = histogramTarget.getBinPosition(idx); } - Kokkos::deep_copy(histogramFine.data, h_dataFine); + Kokkos::deep_copy(histogramTarget.data, h_dataFine); - auto histogramInterp = util::interpolate(histogramCoarse, histogramFine); + util::updateInterpolate(histogramData, histogramTarget); auto h_dataInterp = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramInterp.data); + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramTarget.data); for (auto idx = 0; idx < 30; ++idx) { From de44bf1bc6e89fd1acd3e380a3a6379104b6805a Mon Sep 17 00:00:00 2001 From: julianhille Date: Tue, 27 Jan 2026 15:18:52 +0100 Subject: [PATCH 4/8] updated updateInterpolate test, Ref. #047 --- mrmd/util/interpolation.test.cpp | 43 +++++++++++++++++++------------- 1 file changed, 26 insertions(+), 17 deletions(-) diff --git a/mrmd/util/interpolation.test.cpp b/mrmd/util/interpolation.test.cpp index 9689ab0..acc14c6 100644 --- a/mrmd/util/interpolation.test.cpp +++ b/mrmd/util/interpolation.test.cpp @@ -24,36 +24,45 @@ namespace util { TEST(interpolate, testInterpolate) { - data::MultiHistogram histogramData("histogramData", 0_r, 10_r, 10, 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(histogramData.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) = histogramData.getBinPosition(idx); + h_dataInput(idx, 0) = 0_r; + h_dataInput(idx, 1) = 1_r; + h_dataInput(idx, 2) = histogramInput.getBinPosition(idx); } - Kokkos::deep_copy(histogramData.data, h_dataCoarse); + Kokkos::deep_copy(histogramInput.data, h_dataInput); - auto h_dataFine = Kokkos::create_mirror_view(histogramTarget.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) = histogramTarget.getBinPosition(idx); + h_dataTarget(idx, 0) = 0_r; + h_dataTarget(idx, 1) = 0_r; + h_dataTarget(idx, 2) = 0_r; } - Kokkos::deep_copy(histogramTarget.data, h_dataFine); + Kokkos::deep_copy(histogramTarget.data, h_dataTarget); - util::updateInterpolate(histogramData, histogramTarget); - auto h_dataInterp = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramTarget.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) = histogramTarget.getBinPosition(idx); + } + Kokkos::deep_copy(histogramRef.data, h_dataRef); + + util::updateInterpolate(histogramInput, histogramTarget); + 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)); } } From 60f8089fe9580178c03318f48c4a5f0ed4956eaa Mon Sep 17 00:00:00 2001 From: julianhille Date: Tue, 27 Jan 2026 17:48:24 +0100 Subject: [PATCH 5/8] deleted redundant return, Ref. #047 --- mrmd/util/interpolation.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index 206773e..72a4fad 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -48,8 +48,6 @@ void updateInterpolate(const data::MultiHistogram& inputData, data::MultiHistogr }; Kokkos::parallel_for("MultiHistogram::updateInterpolate", policy, kernel); Kokkos::fence(); - - return; } } // namespace util From 1cd2b38ac2b3a461946fd0fd5c642c69c6161856 Mon Sep 17 00:00:00 2001 From: julianhille Date: Wed, 4 Feb 2026 09:20:58 +0100 Subject: [PATCH 6/8] changed order of parameters in updateInterpolate, Ref. #047 --- mrmd/util/interpolation.cpp | 28 ++++++++++++++-------------- mrmd/util/interpolation.hpp | 2 +- mrmd/util/interpolation.test.cpp | 2 +- 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index 72a4fad..c81dacf 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -18,33 +18,33 @@ namespace mrmd { namespace util { -void updateInterpolate(const data::MultiHistogram& inputData, data::MultiHistogram& inputTarget) +void updateInterpolate(const data::MultiHistogram& target, const data::MultiHistogram& input) { - MRMD_HOST_ASSERT_EQUAL(inputTarget.numHistograms, inputData.numHistograms); + MRMD_HOST_ASSERT_EQUAL(target.numHistograms, input.numHistograms); - auto policy = Kokkos::MDRangePolicy>( - {idx_t(0), idx_t(0)}, {inputTarget.numBins, inputData.numHistograms}); + auto policy = Kokkos::MDRangePolicy>({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 = inputTarget.getBinPosition(binIdx); - idx_t leftBinIdx = inputData.getBin(outputBinPosition - 0.5_r * inputData.binSize); + 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 >= inputData.numBins) + if (leftBinIdx < 0 || rightBinIdx >= input.numBins) { - inputTarget.data(binIdx, histogramIdx) += 0_r; // out of bounds, set to zero + target.data(binIdx, histogramIdx) += 0_r; // out of bounds, set to zero return; } - auto inputDataLeft = inputData.data(leftBinIdx, histogramIdx); - auto inputDataRight = inputData.data(rightBinIdx, histogramIdx); + auto inputLeft = input.data(leftBinIdx, histogramIdx); + auto inputRight = input.data(rightBinIdx, histogramIdx); - inputTarget.data(binIdx, histogramIdx) += lerp( - inputDataLeft, - inputDataRight, - (outputBinPosition - inputData.getBinPosition(leftBinIdx)) * inputData.inverseBinSize); + target.data(binIdx, histogramIdx) += + lerp(inputLeft, + inputRight, + (outputBinPosition - input.getBinPosition(leftBinIdx)) * input.inverseBinSize); }; Kokkos::parallel_for("MultiHistogram::updateInterpolate", policy, kernel); Kokkos::fence(); diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp index 5ec9152..bded04a 100644 --- a/mrmd/util/interpolation.hpp +++ b/mrmd/util/interpolation.hpp @@ -42,7 +42,7 @@ real_t lerp(const real_t& left, const real_t& right, const real_t& factor) * @param target MultiHistogram defining the grid to interpolate onto and containing the data * to be updated by interpolating the data from input. */ -void updateInterpolate(const data::MultiHistogram& input, data::MultiHistogram& target); +void updateInterpolate(const data::MultiHistogram& target, const data::MultiHistogram& input); } // namespace util } // namespace mrmd \ No newline at end of file diff --git a/mrmd/util/interpolation.test.cpp b/mrmd/util/interpolation.test.cpp index acc14c6..1bfc2d9 100644 --- a/mrmd/util/interpolation.test.cpp +++ b/mrmd/util/interpolation.test.cpp @@ -55,7 +55,7 @@ TEST(interpolate, testInterpolate) } Kokkos::deep_copy(histogramRef.data, h_dataRef); - util::updateInterpolate(histogramInput, histogramTarget); + util::updateInterpolate(histogramTarget, histogramInput); h_dataTarget = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramTarget.data); for (auto idx = 0; idx < 30; ++idx) From e886e783da9dfba49e639cfe605925000ad141b7 Mon Sep 17 00:00:00 2001 From: julianhille Date: Mon, 9 Feb 2026 17:31:06 +0100 Subject: [PATCH 7/8] debugged updateInterpolate and added test with out-of-bounds grid positions, Ref. #047 --- mrmd/util/interpolation.cpp | 29 +++++++++------- mrmd/util/interpolation.hpp | 10 ++++++ mrmd/util/interpolation.test.cpp | 58 ++++++++++++++++++++++++++++++-- 3 files changed, 83 insertions(+), 14 deletions(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index c81dacf..3ddc338 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -31,20 +31,25 @@ void updateInterpolate(const data::MultiHistogram& target, const data::MultiHist 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 or exactly at boundary + if (leftBinIdx >= 0 && rightBinIdx < input.numBins) { - target.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 inputLeft = input.data(leftBinIdx, histogramIdx); + auto inputRight = input.data(rightBinIdx, histogramIdx); - target.data(binIdx, histogramIdx) += - lerp(inputLeft, - inputRight, - (outputBinPosition - input.getBinPosition(leftBinIdx)) * input.inverseBinSize); + target.data(binIdx, histogramIdx) += + lerp(inputLeft, + inputRight, + (outputBinPosition - input.getBinPosition(leftBinIdx)) * input.inverseBinSize); + } + else if (isFloatEQ(outputBinPosition, input.getBinPosition(0))) + { + target.data(binIdx, histogramIdx) += input.data(0, histogramIdx); + } + else if (isFloatEQ(outputBinPosition, input.getBinPosition(input.numBins - 1))) + { + target.data(binIdx, histogramIdx) += input.data(input.numBins - 1, histogramIdx); + } }; Kokkos::parallel_for("MultiHistogram::updateInterpolate", policy, kernel); Kokkos::fence(); diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp index bded04a..0a2fdee 100644 --- a/mrmd/util/interpolation.hpp +++ b/mrmd/util/interpolation.hpp @@ -14,6 +14,8 @@ #pragma once +#include + #include "data/MultiHistogram.hpp" #include "util/IsInSymmetricSlab.hpp" @@ -34,6 +36,14 @@ real_t lerp(const real_t& left, const real_t& right, const real_t& factor) return left + (right - left) * factor; } +KOKKOS_INLINE_FUNCTION +bool isFloatEQ(const real_t& a, + const real_t& b, + const real_t& epsilon = std::numeric_limits::epsilon()) +{ + return std::abs(a - b) <= epsilon; +} + /** * Linear interpolation of data contained in input MultiHistogram onto grid of target histogram, * updating the data of target. Data for grid points outside of the grid range of the input diff --git a/mrmd/util/interpolation.test.cpp b/mrmd/util/interpolation.test.cpp index 1bfc2d9..fdba058 100644 --- a/mrmd/util/interpolation.test.cpp +++ b/mrmd/util/interpolation.test.cpp @@ -22,7 +22,7 @@ namespace mrmd { namespace util { -TEST(interpolate, testInterpolate) +TEST(interpolate, preZeroInner) { data::MultiHistogram histogramInput("histogramInput", 0_r, 10_r, 10, 3); data::MultiHistogram histogramTarget("histogramTarget", 0.5_r, 9.5_r, 30, 3); @@ -51,7 +51,61 @@ TEST(interpolate, testInterpolate) { h_dataRef(idx, 0) = 0_r; h_dataRef(idx, 1) = 1_r; - h_dataRef(idx, 2) = histogramTarget.getBinPosition(idx); + 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); From 563863a31f395416210c7d34134a8f461be0f72e Mon Sep 17 00:00:00 2001 From: julianhille Date: Mon, 9 Feb 2026 17:40:42 +0100 Subject: [PATCH 8/8] clarified documentation for updateInterpolate, Ref. #047 --- mrmd/util/interpolation.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp index 0a2fdee..bb20dc0 100644 --- a/mrmd/util/interpolation.hpp +++ b/mrmd/util/interpolation.hpp @@ -46,9 +46,9 @@ bool isFloatEQ(const real_t& a, /** * Linear interpolation of data contained in input MultiHistogram onto grid of target histogram, - * updating the data of target. 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 on coarse grid. + * 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. */