From 92caeac68776f454fe274ef31dd741328b5ff50e Mon Sep 17 00:00:00 2001 From: julianhille Date: Mon, 19 Jan 2026 15:54:36 +0100 Subject: [PATCH 01/15] added interpolation routines, Ref. #047 --- mrmd/util/interpolation.cpp | 81 +++++++++++++++++++++++++++++++++++++ mrmd/util/interpolation.hpp | 67 ++++++++++++++++++++++++++++++ 2 files changed, 148 insertions(+) create mode 100644 mrmd/util/interpolation.cpp create mode 100644 mrmd/util/interpolation.hpp diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp new file mode 100644 index 0000000..083c46e --- /dev/null +++ b/mrmd/util/interpolation.cpp @@ -0,0 +1,81 @@ +// 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) +{ + data::MultiHistogram output( + "interpolated-profile", input.min, input.max, idx_c(grid.extent(0)), input.numHistograms); + const ScalarView& inputGrid = input.createGrid(); + + auto policy = Kokkos::MDRangePolicy>( + {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 + auto rightBin = findRightBin(inputGrid, grid(binIdx)); + auto leftBin = rightBin - 1; + + // handle boundaries + if (leftBin < 0 || rightBin >= input.numBins) + { + output.data(binIdx, histogramIdx) = 0.0_r; // out of bounds, set to zero + return; + } + + auto inputDataLeft = input.data(leftBin, histogramIdx); + auto inputDataRight = input.data(rightBin, histogramIdx); + + output.data(binIdx, histogramIdx) = + lerp(inputDataLeft, + inputDataRight, + (grid(binIdx) - inputGrid(leftBin)) * input.inverseBinSize); + }; + Kokkos::parallel_for("MultiHistogram::interpolate", policy, kernel); + Kokkos::fence(); + + return output; +} + +data::MultiHistogram constrainToSymmetricSlab(const data::MultiHistogram& input, + const util::IsInSymmetricSlab& applicationRegion) +{ + data::MultiHistogram constrainedProfile( + "constrained-profile", input.min, input.max, input.numBins, input.numHistograms); + + auto policy = Kokkos::MDRangePolicy>({idx_t(0), idx_t(0)}, + {input.numBins, input.numHistograms}); + auto kernel = KOKKOS_LAMBDA(const idx_t binIdx, const idx_t histogramIdx) + { + if (!applicationRegion(input.getBinPosition(binIdx), 0_r, 0_r)) + { + constrainedProfile.data(binIdx, histogramIdx) = 0_r; + } + else + { + constrainedProfile.data(binIdx, histogramIdx) = input.data(binIdx, histogramIdx); + } + }; + Kokkos::parallel_for("MultiHistogram::constrainToApplicationRegion", policy, kernel); + Kokkos::fence(); + + return constrainedProfile; +} +} // namespace util +} // namespace mrmd \ No newline at end of file diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp new file mode 100644 index 0000000..80ebab1 --- /dev/null +++ b/mrmd/util/interpolation.hpp @@ -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. + +#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; +} + +/** + * Find the index of the first bin in the grid that is greater than or equal to the value. + * @param grid grid to search in + * @param value value to find + * @return index of the right bin + */ +KOKKOS_INLINE_FUNCTION +idx_t findRightBin(const ScalarView& grid, const real_t& value) +{ + idx_t idx = 0; + for (; idx < idx_c(grid.extent(0)) && grid(idx) < value; ++idx) + { + // empty loop body + } + + return idx; +} + +/** + * Interpolate data values of MultiHistogram to a new grid. + * @param input input MultiHistogram + * @param grid grid to interpolate to + * @return MultiHistogram with interpolated values + */ +data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid); + +data::MultiHistogram constrainToSymmetricSlab(const data::MultiHistogram& input, + const util::IsInSymmetricSlab& applicationRegion); + +} // namespace util +} // namespace mrmd \ No newline at end of file From 214a0b4d08687ac991e95d97d7c5f71a19dcdac1 Mon Sep 17 00:00:00 2001 From: julianhille Date: Tue, 20 Jan 2026 14:02:58 +0100 Subject: [PATCH 02/15] made interpolation consistent with changes to createGrid, Ref. #047 --- mrmd/util/interpolation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index 083c46e..9694563 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -22,7 +22,7 @@ data::MultiHistogram interpolate(const data::MultiHistogram& input, const Scalar { data::MultiHistogram output( "interpolated-profile", input.min, input.max, idx_c(grid.extent(0)), input.numHistograms); - const ScalarView& inputGrid = input.createGrid(); + const ScalarView& inputGrid = createGrid(input); auto policy = Kokkos::MDRangePolicy>( {idx_t(0), idx_t(0)}, {idx_c(grid.extent(0)), input.numHistograms}); From f0731a8640ce476626c95eb9ddc30e31b004c767 Mon Sep 17 00:00:00 2001 From: julianhille Date: Tue, 20 Jan 2026 14:03:21 +0100 Subject: [PATCH 03/15] added test for interpolation, Ref. #047 --- mrmd/util/CMakeLists.txt | 3 ++ mrmd/util/interpolation.test.cpp | 57 ++++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+) create mode 100644 mrmd/util/interpolation.test.cpp diff --git a/mrmd/util/CMakeLists.txt b/mrmd/util/CMakeLists.txt index b17b54c..91cd91e 100644 --- a/mrmd/util/CMakeLists.txt +++ b/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) \ No newline at end of file diff --git a/mrmd/util/interpolation.test.cpp b/mrmd/util/interpolation.test.cpp new file mode 100644 index 0000000..a1ead45 --- /dev/null +++ b/mrmd/util/interpolation.test.cpp @@ -0,0 +1,57 @@ +// 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 + +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 < 10; ++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 \ No newline at end of file From 5ed24008c0a31a0cd3ed6acaabacad5475f7f015 Mon Sep 17 00:00:00 2001 From: julianhille Date: Tue, 20 Jan 2026 14:06:12 +0100 Subject: [PATCH 04/15] renamed kokkos parallel_for label, Ref. #047 --- mrmd/util/interpolation.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index 9694563..0cf830f 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -54,7 +54,7 @@ data::MultiHistogram interpolate(const data::MultiHistogram& input, const Scalar } data::MultiHistogram constrainToSymmetricSlab(const data::MultiHistogram& input, - const util::IsInSymmetricSlab& applicationRegion) + const util::IsInSymmetricSlab& applicationRegion) { data::MultiHistogram constrainedProfile( "constrained-profile", input.min, input.max, input.numBins, input.numHistograms); @@ -72,7 +72,7 @@ data::MultiHistogram constrainToSymmetricSlab(const data::MultiHistogram& input, constrainedProfile.data(binIdx, histogramIdx) = input.data(binIdx, histogramIdx); } }; - Kokkos::parallel_for("MultiHistogram::constrainToApplicationRegion", policy, kernel); + Kokkos::parallel_for("MultiHistogram::constrainToSymmetricSlab", policy, kernel); Kokkos::fence(); return constrainedProfile; From 3749eabc9190962dd3e76ffa9d73589d6f05a727 Mon Sep 17 00:00:00 2001 From: julianhille Date: Tue, 20 Jan 2026 14:07:08 +0100 Subject: [PATCH 05/15] clang format fixes, Ref. #047 --- mrmd/util/interpolation.hpp | 2 +- mrmd/util/interpolation.test.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp index 80ebab1..328e1ef 100644 --- a/mrmd/util/interpolation.hpp +++ b/mrmd/util/interpolation.hpp @@ -61,7 +61,7 @@ idx_t findRightBin(const ScalarView& grid, const real_t& value) data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid); data::MultiHistogram constrainToSymmetricSlab(const data::MultiHistogram& input, - const util::IsInSymmetricSlab& applicationRegion); + const util::IsInSymmetricSlab& applicationRegion); } // 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 a1ead45..fbc17f3 100644 --- a/mrmd/util/interpolation.test.cpp +++ b/mrmd/util/interpolation.test.cpp @@ -44,7 +44,8 @@ TEST(interpolate, testInterpolate) 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); + auto h_dataInterp = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramInterp.data); for (auto idx = 0; idx < 10; ++idx) { From 07f199dd3d7b58924398a9659b4cf68f65600e0b Mon Sep 17 00:00:00 2001 From: julianhille Date: Tue, 20 Jan 2026 14:27:59 +0100 Subject: [PATCH 06/15] added test for constrainToSymmetricSlab function, Ref. #047 --- mrmd/util/interpolation.test.cpp | 46 +++++++++++++++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/mrmd/util/interpolation.test.cpp b/mrmd/util/interpolation.test.cpp index fbc17f3..16ddbdf 100644 --- a/mrmd/util/interpolation.test.cpp +++ b/mrmd/util/interpolation.test.cpp @@ -16,6 +16,8 @@ #include +#include "util/IsInSymmetricSlab.hpp" + namespace mrmd { namespace util @@ -47,12 +49,54 @@ TEST(interpolate, testInterpolate) auto h_dataInterp = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramInterp.data); - for (auto idx = 0; idx < 10; ++idx) + 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)); } } + +TEST(interpolate, testConstrainToSymmetricSlab) +{ + data::MultiHistogram histogramBare("histogram", 0_r, 10_r, 30, 3); + + auto gridBare = createGrid(histogramBare); + + auto h_dataBare = Kokkos::create_mirror_view(histogramBare.data); + for (auto idx = 0; idx < 30; ++idx) + { + h_dataBare(idx, 0) = 0_r; + h_dataBare(idx, 1) = 1_r; + h_dataBare(idx, 2) = histogramBare.getBinPosition(idx); + } + Kokkos::deep_copy(histogramBare.data, h_dataBare); + + // constrain data to slab automatically in decive space + auto applicationRegion = util::IsInSymmetricSlab(Point3D{5_r, 0_r, 0_r}, 2_r, 4_r); + auto histogramConstr = constrainToSymmetricSlab(histogramBare, applicationRegion); + auto h_dataConstr = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramConstr.data); + + // constrain data to slab manually in host space + for (auto idx = 0; idx < 30; ++idx) + { + real_t binPosition = gridBare(idx); + if (binPosition <= 1_r || (binPosition >= 3_r && binPosition <= 7_r) || binPosition >= 9_r) + { + h_dataBare(idx, 0) = 0_r; + h_dataBare(idx, 1) = 0_r; + h_dataBare(idx, 2) = 0_r; + } + } + + for (auto idx = 0; idx < 30; ++idx) + { + EXPECT_FLOAT_EQ(h_dataBare(idx, 0), h_dataConstr(idx, 0)); + EXPECT_FLOAT_EQ(h_dataBare(idx, 1), h_dataConstr(idx, 1)); + EXPECT_FLOAT_EQ(h_dataBare(idx, 2), h_dataConstr(idx, 2)); + } +} + } // namespace util } // namespace mrmd \ No newline at end of file From 0704aeadd8ba1d9020f63dc9d2755145b72a5b55 Mon Sep 17 00:00:00 2001 From: julianhille Date: Wed, 21 Jan 2026 11:47:14 +0100 Subject: [PATCH 07/15] added assertion that target grid matches output grid, Ref. #047 --- mrmd/util/interpolation.cpp | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index 0cf830f..b37701a 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -20,8 +20,21 @@ 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", input.min, input.max, idx_c(grid.extent(0)), input.numHistograms); + "interpolated-profile", gridMin, gridMax, idx_c(grid.extent(0)), input.numHistograms); + + const ScalarView& outputGrid = createGrid(output); + + MRMD_HOST_ASSERT_EQUAL(outputGrid.extent(0), grid.extent(0), "Output grid size mismatch!"); + for (idx_t idx = 0; idx < idx_c(outputGrid.extent(0)); ++idx) + { + MRMD_HOST_ASSERT_EQUAL(outputGrid(idx), grid(idx), "Output grid mismatch!"); + } + const ScalarView& inputGrid = createGrid(input); auto policy = Kokkos::MDRangePolicy>( From 536d70c1e31832d3e776a27db5149572c564867b Mon Sep 17 00:00:00 2001 From: julianhille Date: Wed, 21 Jan 2026 12:03:08 +0100 Subject: [PATCH 08/15] clarified documentation, 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 328e1ef..1debff5 100644 --- a/mrmd/util/interpolation.hpp +++ b/mrmd/util/interpolation.hpp @@ -35,10 +35,11 @@ real_t lerp(const real_t& left, const real_t& right, const real_t& factor) } /** - * Find the index of the first bin in the grid that is greater than or equal to the value. + * Find the index of the first bin in the grid that has a position greater than or equal to the + * value. * @param grid grid to search in * @param value value to find - * @return index of the right bin + * @return index of the right-hand bin */ KOKKOS_INLINE_FUNCTION idx_t findRightBin(const ScalarView& grid, const real_t& value) @@ -48,7 +49,6 @@ idx_t findRightBin(const ScalarView& grid, const real_t& value) { // empty loop body } - return idx; } From 277d61f5b001edde995e793803217a65b0da1f4c Mon Sep 17 00:00:00 2001 From: julianhille Date: Thu, 22 Jan 2026 14:08:53 +0100 Subject: [PATCH 09/15] rationalized away getRightBin function, Ref. #047 --- mrmd/util/interpolation.cpp | 17 ++++++++++------- mrmd/util/interpolation.hpp | 18 ------------------ 2 files changed, 10 insertions(+), 25 deletions(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index b37701a..312157b 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -42,23 +42,26 @@ data::MultiHistogram interpolate(const data::MultiHistogram& input, const Scalar auto kernel = KOKKOS_LAMBDA(const idx_t binIdx, const idx_t histogramIdx) { // find the two enclosing bins in the input histogram - auto rightBin = findRightBin(inputGrid, grid(binIdx)); - auto leftBin = rightBin - 1; + real_t outputBinPosition = output.getBinPosition(binIdx); + idx_t inputBinIdx = input.getBin(outputBinPosition); + idx_t leftBinIdx = + inputBinIdx + idx_c(std::floor(outputBinPosition - input.getBinPosition(inputBinIdx))); + idx_t rightBinIdx = leftBinIdx + 1; // handle boundaries - if (leftBin < 0 || rightBin >= input.numBins) + if (leftBinIdx < 0 || rightBinIdx >= input.numBins) { - output.data(binIdx, histogramIdx) = 0.0_r; // out of bounds, set to zero + output.data(binIdx, histogramIdx) = 0_r; // out of bounds, set to zero return; } - auto inputDataLeft = input.data(leftBin, histogramIdx); - auto inputDataRight = input.data(rightBin, histogramIdx); + auto inputDataLeft = input.data(leftBinIdx, histogramIdx); + auto inputDataRight = input.data(rightBinIdx, histogramIdx); output.data(binIdx, histogramIdx) = lerp(inputDataLeft, inputDataRight, - (grid(binIdx) - inputGrid(leftBin)) * input.inverseBinSize); + (outputBinPosition - inputGrid(leftBinIdx)) * input.inverseBinSize); }; Kokkos::parallel_for("MultiHistogram::interpolate", policy, kernel); Kokkos::fence(); diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp index 1debff5..fffcfbb 100644 --- a/mrmd/util/interpolation.hpp +++ b/mrmd/util/interpolation.hpp @@ -34,24 +34,6 @@ real_t lerp(const real_t& left, const real_t& right, const real_t& factor) return left + (right - left) * factor; } -/** - * Find the index of the first bin in the grid that has a position greater than or equal to the - * value. - * @param grid grid to search in - * @param value value to find - * @return index of the right-hand bin - */ -KOKKOS_INLINE_FUNCTION -idx_t findRightBin(const ScalarView& grid, const real_t& value) -{ - idx_t idx = 0; - for (; idx < idx_c(grid.extent(0)) && grid(idx) < value; ++idx) - { - // empty loop body - } - return idx; -} - /** * Interpolate data values of MultiHistogram to a new grid. * @param input input MultiHistogram From fbbe05abfc99a2f1eb4b4afef1f16f5b8dce582a Mon Sep 17 00:00:00 2001 From: julianhille Date: Thu, 22 Jan 2026 14:32:20 +0100 Subject: [PATCH 10/15] simplified search for enclosing bins, Ref. #047 --- mrmd/util/interpolation.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index 312157b..a78f580 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -43,9 +43,7 @@ data::MultiHistogram interpolate(const data::MultiHistogram& input, const Scalar { // find the two enclosing bins in the input histogram real_t outputBinPosition = output.getBinPosition(binIdx); - idx_t inputBinIdx = input.getBin(outputBinPosition); - idx_t leftBinIdx = - inputBinIdx + idx_c(std::floor(outputBinPosition - input.getBinPosition(inputBinIdx))); + idx_t leftBinIdx = input.getBin(outputBinPosition - 0.5_r * input.binSize); idx_t rightBinIdx = leftBinIdx + 1; // handle boundaries From df840f72cbf2c4b585f972d6a662f874bc081f7e Mon Sep 17 00:00:00 2001 From: julianhille Date: Thu, 22 Jan 2026 14:41:06 +0100 Subject: [PATCH 11/15] used grid rather than output histogram for outputBinPosition, Ref. #047 --- mrmd/util/interpolation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index a78f580..d364d31 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -42,7 +42,7 @@ data::MultiHistogram interpolate(const data::MultiHistogram& input, const Scalar 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); + real_t outputBinPosition = grid(binIdx); idx_t leftBinIdx = input.getBin(outputBinPosition - 0.5_r * input.binSize); idx_t rightBinIdx = leftBinIdx + 1; From 44727bb61990c70ec46778b74327791bd97ce0ba Mon Sep 17 00:00:00 2001 From: julianhille Date: Thu, 22 Jan 2026 14:49:16 +0100 Subject: [PATCH 12/15] added harsh boundary conditions to function docstring, Ref. #047 --- mrmd/util/interpolation.hpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp index fffcfbb..22b2005 100644 --- a/mrmd/util/interpolation.hpp +++ b/mrmd/util/interpolation.hpp @@ -35,10 +35,11 @@ real_t lerp(const real_t& left, const real_t& right, const real_t& factor) } /** - * Interpolate data values of MultiHistogram to a new grid. - * @param input input MultiHistogram - * @param grid grid to interpolate to - * @return MultiHistogram with interpolated values + * 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); From 8325d1fa4986fcba1a15fe1ecd97617f50e9dbc6 Mon Sep 17 00:00:00 2001 From: julianhille Date: Thu, 22 Jan 2026 14:55:10 +0100 Subject: [PATCH 13/15] got rid of additional grids, Ref. #047 --- mrmd/util/interpolation.cpp | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index d364d31..495c7ff 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -27,16 +27,12 @@ data::MultiHistogram interpolate(const data::MultiHistogram& input, const Scalar data::MultiHistogram output( "interpolated-profile", gridMin, gridMax, idx_c(grid.extent(0)), input.numHistograms); - const ScalarView& outputGrid = createGrid(output); - - MRMD_HOST_ASSERT_EQUAL(outputGrid.extent(0), grid.extent(0), "Output grid size mismatch!"); - for (idx_t idx = 0; idx < idx_c(outputGrid.extent(0)); ++idx) + MRMD_HOST_ASSERT_EQUAL(output.numBins, grid.extent(0), "Output grid size mismatch!"); + for (idx_t idx = 0; idx < idx_c(output.numBins); ++idx) { - MRMD_HOST_ASSERT_EQUAL(outputGrid(idx), grid(idx), "Output grid mismatch!"); + MRMD_HOST_ASSERT_EQUAL(output.getBinPosition(idx), grid(idx), "Output grid mismatch!"); } - const ScalarView& inputGrid = createGrid(input); - auto policy = Kokkos::MDRangePolicy>( {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) @@ -59,7 +55,7 @@ data::MultiHistogram interpolate(const data::MultiHistogram& input, const Scalar output.data(binIdx, histogramIdx) = lerp(inputDataLeft, inputDataRight, - (outputBinPosition - inputGrid(leftBinIdx)) * input.inverseBinSize); + (outputBinPosition - input.getBinPosition(leftBinIdx)) * input.inverseBinSize); }; Kokkos::parallel_for("MultiHistogram::interpolate", policy, kernel); Kokkos::fence(); From 9a97003888dca97217811681feb0926bd44fccfc Mon Sep 17 00:00:00 2001 From: julianhille Date: Thu, 22 Jan 2026 15:03:38 +0100 Subject: [PATCH 14/15] debug signedness of integers, Ref. #047 --- mrmd/util/interpolation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index 495c7ff..8f0c258 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -27,7 +27,7 @@ data::MultiHistogram interpolate(const data::MultiHistogram& input, const Scalar data::MultiHistogram output( "interpolated-profile", gridMin, gridMax, idx_c(grid.extent(0)), input.numHistograms); - MRMD_HOST_ASSERT_EQUAL(output.numBins, grid.extent(0), "Output grid size mismatch!"); + 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!"); From bc3fb564d0d26a4e1895ec41efef55cfb3e4b3bf Mon Sep 17 00:00:00 2001 From: julianhille Date: Thu, 22 Jan 2026 15:12:41 +0100 Subject: [PATCH 15/15] got rid of redundant constrainToSymmetricSlab function, Ref. #047 --- mrmd/util/interpolation.cpp | 24 ------------------- mrmd/util/interpolation.hpp | 3 --- mrmd/util/interpolation.test.cpp | 41 -------------------------------- 3 files changed, 68 deletions(-) diff --git a/mrmd/util/interpolation.cpp b/mrmd/util/interpolation.cpp index 8f0c258..59ce613 100644 --- a/mrmd/util/interpolation.cpp +++ b/mrmd/util/interpolation.cpp @@ -63,29 +63,5 @@ data::MultiHistogram interpolate(const data::MultiHistogram& input, const Scalar return output; } -data::MultiHistogram constrainToSymmetricSlab(const data::MultiHistogram& input, - const util::IsInSymmetricSlab& applicationRegion) -{ - data::MultiHistogram constrainedProfile( - "constrained-profile", input.min, input.max, input.numBins, input.numHistograms); - - auto policy = Kokkos::MDRangePolicy>({idx_t(0), idx_t(0)}, - {input.numBins, input.numHistograms}); - auto kernel = KOKKOS_LAMBDA(const idx_t binIdx, const idx_t histogramIdx) - { - if (!applicationRegion(input.getBinPosition(binIdx), 0_r, 0_r)) - { - constrainedProfile.data(binIdx, histogramIdx) = 0_r; - } - else - { - constrainedProfile.data(binIdx, histogramIdx) = input.data(binIdx, histogramIdx); - } - }; - Kokkos::parallel_for("MultiHistogram::constrainToSymmetricSlab", policy, kernel); - Kokkos::fence(); - - return constrainedProfile; -} } // namespace util } // namespace mrmd \ No newline at end of file diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp index 22b2005..41987b2 100644 --- a/mrmd/util/interpolation.hpp +++ b/mrmd/util/interpolation.hpp @@ -43,8 +43,5 @@ real_t lerp(const real_t& left, const real_t& right, const real_t& factor) */ data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid); -data::MultiHistogram constrainToSymmetricSlab(const data::MultiHistogram& input, - const util::IsInSymmetricSlab& applicationRegion); - } // 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 16ddbdf..e136467 100644 --- a/mrmd/util/interpolation.test.cpp +++ b/mrmd/util/interpolation.test.cpp @@ -57,46 +57,5 @@ TEST(interpolate, testInterpolate) } } -TEST(interpolate, testConstrainToSymmetricSlab) -{ - data::MultiHistogram histogramBare("histogram", 0_r, 10_r, 30, 3); - - auto gridBare = createGrid(histogramBare); - - auto h_dataBare = Kokkos::create_mirror_view(histogramBare.data); - for (auto idx = 0; idx < 30; ++idx) - { - h_dataBare(idx, 0) = 0_r; - h_dataBare(idx, 1) = 1_r; - h_dataBare(idx, 2) = histogramBare.getBinPosition(idx); - } - Kokkos::deep_copy(histogramBare.data, h_dataBare); - - // constrain data to slab automatically in decive space - auto applicationRegion = util::IsInSymmetricSlab(Point3D{5_r, 0_r, 0_r}, 2_r, 4_r); - auto histogramConstr = constrainToSymmetricSlab(histogramBare, applicationRegion); - auto h_dataConstr = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), histogramConstr.data); - - // constrain data to slab manually in host space - for (auto idx = 0; idx < 30; ++idx) - { - real_t binPosition = gridBare(idx); - if (binPosition <= 1_r || (binPosition >= 3_r && binPosition <= 7_r) || binPosition >= 9_r) - { - h_dataBare(idx, 0) = 0_r; - h_dataBare(idx, 1) = 0_r; - h_dataBare(idx, 2) = 0_r; - } - } - - for (auto idx = 0; idx < 30; ++idx) - { - EXPECT_FLOAT_EQ(h_dataBare(idx, 0), h_dataConstr(idx, 0)); - EXPECT_FLOAT_EQ(h_dataBare(idx, 1), h_dataConstr(idx, 1)); - EXPECT_FLOAT_EQ(h_dataBare(idx, 2), h_dataConstr(idx, 2)); - } -} - } // namespace util } // namespace mrmd \ No newline at end of file