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.cpp b/mrmd/util/interpolation.cpp new file mode 100644 index 0000000..59ce613 --- /dev/null +++ b/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>( + {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 \ No newline at end of file diff --git a/mrmd/util/interpolation.hpp b/mrmd/util/interpolation.hpp new file mode 100644 index 0000000..41987b2 --- /dev/null +++ b/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 \ 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..e136467 --- /dev/null +++ b/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 + +#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 \ No newline at end of file