From cbe2a37c3c468a19ea9ad9ca3239ab4bd3a513fd Mon Sep 17 00:00:00 2001 From: Phil Ahrenkiel Date: Tue, 20 Jun 2023 16:34:59 -0600 Subject: [PATCH 1/8] Update .gitignore --- .gitignore | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index c87f188..d8aa5cb 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,5 @@ builds/ build/ .idea/ .vs/ -.vscode/ \ No newline at end of file +.vscode/ +vendor/ \ No newline at end of file From a6ccf0ccbb9985a16d81140c17cd474dda7ebd69 Mon Sep 17 00:00:00 2001 From: Phil Ahrenkiel Date: Wed, 21 Jun 2023 09:07:10 -0600 Subject: [PATCH 2/8] Remove extaneous files. --- include/btwxt/grid-axis.h | 13 +- src/grid-axis.cpp | 42 ++-- ...regular-grid-interpolator-implementation.h | 6 +- src/regular-grid-interpolator.cpp | 72 +++--- test/btwxt-tests.cpp | 71 ++++-- test/grid-axis-tests.cpp | 43 +++- test/implementation-tests.cpp | 213 ++++++++++++------ 7 files changed, 315 insertions(+), 145 deletions(-) diff --git a/include/btwxt/grid-axis.h b/include/btwxt/grid-axis.h index 79d443b..08e6160 100644 --- a/include/btwxt/grid-axis.h +++ b/include/btwxt/grid-axis.h @@ -58,8 +58,9 @@ class GridAxis { return extrapolation_limits; } - [[nodiscard]] const std::vector& - get_cubic_spacing_ratios(std::size_t floor_or_ceiling) const; + [[nodiscard]] const std::vector>& + get_cubic_spacing_ratios(std::size_t elem_index) const; + std::string name; private: @@ -67,11 +68,7 @@ class GridAxis { Method extrapolation_method {Method::constant}; Method interpolation_method {Method::linear}; std::pair extrapolation_limits {-DBL_MAX, DBL_MAX}; - std::vector> - cubic_spacing_ratios; // Used for cubic interpolation. Outer vector is size 2: 0: spacing - // for the floor, 1: spacing for the ceiling. Inner vector is length - // of axis values, but the floor vector doesn't use the first entry - // and the ceiling doesn't use the last entry. + std::vector>> cubic_spacing_ratios; std::shared_ptr logger; void calculate_cubic_spacing_ratios(); void check_grid_sorted(); @@ -107,4 +104,4 @@ std::vector linspace(double start, double stop, std::size_t number_of_po } // namespace Btwxt -#endif // define GRID_AXIS_H_ \ No newline at end of file +#endif // define GRID_AXIS_H_ diff --git a/src/grid-axis.cpp b/src/grid-axis.cpp index 2a5519b..844be95 100644 --- a/src/grid-axis.cpp +++ b/src/grid-axis.cpp @@ -18,7 +18,8 @@ GridAxis::GridAxis(std::vector values_in, , interpolation_method(interpolation_method) , extrapolation_limits(std::move(extrapolation_limits)) , cubic_spacing_ratios( - 2, std::vector(std::max(static_cast(values.size()) - 1, 0), 1.0)) + std::max(static_cast(values.size()) - 1, 0), + {{-1.0, 0.0}, {0.0, 1.0}, {1.0, 0.0}, {0.0, -1.0}}) , logger(logger_in) { if (values.empty()) { @@ -86,23 +87,38 @@ void GridAxis::calculate_cubic_spacing_ratios() if (interpolation_method == Method::linear) { return; } - static constexpr std::size_t floor = 0; - static constexpr std::size_t ceiling = 1; - for (std::size_t i = 0; i < values.size() - 1; i++) { - double center_spacing = values[i + 1] - values[i]; - if (i != 0) { - cubic_spacing_ratios[floor][i] = center_spacing / (values[i + 1] - values[i - 1]); - } - if (i + 2 != values.size()) { - cubic_spacing_ratios[ceiling][i] = center_spacing / (values[i + 2] - values[i]); + for (std::size_t i_elem = 0; i_elem < values.size() - 1; i_elem++) + { + double w_0 = values[i_elem + 1] - values[i_elem]; + if (i_elem > 0) + { + double w_m1 = values[i_elem] - values[i_elem - 1]; + double s0_m = (w_0 - w_m1) / w_m1; + double s1_m = w_m1 / (w_0 + w_m1); + + cubic_spacing_ratios[i_elem][0].first = -(s0_m + s1_m); + cubic_spacing_ratios[i_elem][1].first = s0_m; + cubic_spacing_ratios[i_elem][2].first = s1_m; + cubic_spacing_ratios[i_elem][3].first = 0.0; + } + if (i_elem + 2 < values.size()) + { + double w_1 = values[i_elem + 2] - values[i_elem + 1]; + double s0_p = w_1 / (w_0 + w_1); + double s1_p = (w_0 - w_1) / w_1; + + cubic_spacing_ratios[i_elem][0].second = 0.0; + cubic_spacing_ratios[i_elem][1].second = s0_p; + cubic_spacing_ratios[i_elem][2].second = s1_p; + cubic_spacing_ratios[i_elem][3].second = -(s0_p + s1_p); } } } -const std::vector& -GridAxis::get_cubic_spacing_ratios(const std::size_t floor_or_ceiling) const +const std::vector>& +GridAxis::get_cubic_spacing_ratios(const std::size_t elem_index) const { - return cubic_spacing_ratios[floor_or_ceiling]; + return cubic_spacing_ratios[elem_index]; } void GridAxis::check_grid_sorted() diff --git a/src/regular-grid-interpolator-implementation.h b/src/regular-grid-interpolator-implementation.h index e8eb189..327a58b 100644 --- a/src/regular-grid-interpolator-implementation.h +++ b/src/regular-grid-interpolator-implementation.h @@ -137,11 +137,11 @@ class RegularGridInterpolatorImplementation { return hypercube; }; - [[nodiscard]] const std::vector& - get_axis_cubic_spacing_ratios(std::size_t axis_index, std::size_t floor_or_ceiling) const + [[nodiscard]] const std::vector>& + get_axis_cubic_spacing_ratios(std::size_t axis_index, std::size_t elem_index) const { check_axis_index(axis_index, "get axis cubic spacing ratios"); - return grid_axes[axis_index].get_cubic_spacing_ratios(floor_or_ceiling); + return grid_axes[axis_index].get_cubic_spacing_ratios(elem_index); }; const std::vector& get_grid_point_data(const std::vector& coordinates); diff --git a/src/regular-grid-interpolator.cpp b/src/regular-grid-interpolator.cpp index 96a10a9..b79f59d 100644 --- a/src/regular-grid-interpolator.cpp +++ b/src/regular-grid-interpolator.cpp @@ -97,7 +97,7 @@ RegularGridInterpolator::RegularGridInterpolator(const std::vector& gr { } -RegularGridInterpolator::~RegularGridInterpolator() = default; +RegularGridInterpolator::~RegularGridInterpolator() = default; RegularGridInterpolatorImplementation::RegularGridInterpolatorImplementation( const std::vector& grid_axes, @@ -119,7 +119,7 @@ RegularGridInterpolatorImplementation::RegularGridInterpolatorImplementation( , weighting_factors(number_of_grid_axes, std::vector(4, 0.)) , results(number_of_grid_point_data_sets, 0.) , interpolation_coefficients(number_of_grid_axes, std::vector(2, 0.)) - , cubic_slope_coefficients(number_of_grid_axes, std::vector(2, 0.)) + , cubic_slope_coefficients(number_of_grid_axes, std::vector(4, 0.)) , logger(logger) { set_axis_sizes(); @@ -678,39 +678,55 @@ void RegularGridInterpolatorImplementation::calculate_interpolation_coefficients { static constexpr std::size_t floor = 0; static constexpr std::size_t ceiling = 1; - for (std::size_t axis_index = 0; axis_index < number_of_grid_axes; axis_index++) { + + for (std::size_t axis_index = 0; axis_index < number_of_grid_axes; axis_index++) + { double mu = floor_to_ceiling_fractions[axis_index]; - if (methods[axis_index] == Method::cubic) { - interpolation_coefficients[axis_index][floor] = 2 * mu * mu * mu - 3 * mu * mu + 1; - interpolation_coefficients[axis_index][ceiling] = -2 * mu * mu * mu + 3 * mu * mu; - cubic_slope_coefficients[axis_index][floor] = - (mu * mu * mu - 2 * mu * mu + mu) * - get_axis_cubic_spacing_ratios(axis_index, - floor)[floor_grid_point_coordinates[axis_index]]; - cubic_slope_coefficients[axis_index][ceiling] = - (mu * mu * mu - mu * mu) * - get_axis_cubic_spacing_ratios(axis_index, - ceiling)[floor_grid_point_coordinates[axis_index]]; + + if (methods[axis_index] == Method::cubic) + { + std::vector linear_interpolation_factors = {0.0, 1.0 - mu, mu, 0.0}; + + double coef = 1- 2 * mu; //(1 - mu) - mu; + std::vector cubic_interpolation_factors = {0.0, coef, -coef, 0.0}; + + interpolation_coefficients[axis_index][floor] = linear_interpolation_factors[1] + + (1 - mu) * mu * cubic_interpolation_factors[1]; + interpolation_coefficients[axis_index][ceiling] = linear_interpolation_factors[2] + + (1 - mu) * mu * cubic_interpolation_factors[2]; + + std::vector cubic_slope_factors(4, 0.0); + auto& ratios = get_axis_cubic_spacing_ratios(axis_index, floor_grid_point_coordinates[axis_index]); + for (std::size_t i = 0; i < 4; ++i) + { + cubic_slope_factors[i] = (1 - mu) * ratios[i].first + mu * ratios[i].second; + cubic_slope_coefficients[axis_index][i] = (1 - mu) * mu * cubic_slope_factors[i]; + } + for (std::size_t i = 0; i < 4; ++i) + { + weighting_factors[axis_index][i] = linear_interpolation_factors[i] + + (1 - mu) * mu * (cubic_interpolation_factors[i] + cubic_slope_factors[i]); + } + continue; } - else { - if (methods[axis_index] == Method::constant) { + + { + if (methods[axis_index] == Method::constant) + { mu = mu < 0 ? 0 : 1; } + + std::vector linear_interpolation_factors = {0.0, 1.0 - mu, mu, 0.0}; + interpolation_coefficients[axis_index][floor] = 1 - mu; interpolation_coefficients[axis_index][ceiling] = mu; - cubic_slope_coefficients[axis_index][floor] = 0.0; - cubic_slope_coefficients[axis_index][ceiling] = 0.0; + + for (std::size_t i = 0; i < 4; ++i) + { + cubic_slope_coefficients[axis_index][i] = 0.0; + weighting_factors[axis_index][i] = linear_interpolation_factors[i]; + } } - weighting_factors[axis_index][0] = - -cubic_slope_coefficients[axis_index][floor]; // point below floor (-1) - weighting_factors[axis_index][1] = - interpolation_coefficients[axis_index][floor] - - cubic_slope_coefficients[axis_index][ceiling]; // floor (0) - weighting_factors[axis_index][2] = - interpolation_coefficients[axis_index][ceiling] + - cubic_slope_coefficients[axis_index][floor]; // ceiling (1) - weighting_factors[axis_index][3] = - cubic_slope_coefficients[axis_index][ceiling]; // point above ceiling (2) } } diff --git a/test/btwxt-tests.cpp b/test/btwxt-tests.cpp index 398e906..d70aece 100644 --- a/test/btwxt-tests.cpp +++ b/test/btwxt-tests.cpp @@ -16,6 +16,40 @@ namespace Btwxt { +double fCubic(double x) +{ + return 1.+x*x*x; +} + +TEST(BasicSplineTest, cubic_interp) +{ + std::vector x_values={1., 2., 4., 5.}; + std::vector y_values; + + for(auto &x_value:x_values) + { + y_values.push_back(fCubic(x_value)); + } + std::vector> values; + values.push_back(y_values); + + Method extrapolation_method = Method::linear; + Method interpolation_method = Method::cubic; + + std::pair extrapolation_limits{-DBL_MAX, DBL_MAX}; + + GridAxis axis(x_values, "x", interpolation_method, extrapolation_method, extrapolation_limits); + RegularGridInterpolator my_interpolator({axis}, values); + + const double epsilon = 0.0001; + double x(3); + double expected_value=fCubic(x); + std::vector target{x}; + my_interpolator.set_target(target); + double result = my_interpolator.get_values_at_target()[0]; + EXPECT_NEAR(result, expected_value, epsilon); +} + TEST_F(FunctionFixture, scipy_3d_grid) { // Based on @@ -100,15 +134,23 @@ TEST(Constructors, test_all_constructors) TEST_F(GridFixture, four_point_1d_cubic_interpolate) { + grid = {{1., 2., 4., 5.}}; + + data_sets.clear(); + std::vector data_values; + for(auto &grid_value: grid[0]) + { + data_values.push_back(fCubic(grid_value)); + } + std::vector> values; + data_sets.push_back(data_values); - grid = {{0, 2, 5, 10}}; - data_sets = {{6, 5, 4, 3}}; - target = {2.5}; + target = {3.}; setup(); interpolator.set_axis_interpolation_method(0, Method::cubic); - const double expected_value = 4.804398; + const double expected_value = fCubic(target[0]); const double epsilon = 0.0001; double result = interpolator.get_values_at_target(target)[0]; EXPECT_NEAR(result, expected_value, epsilon); @@ -310,7 +352,7 @@ TEST_F(Grid2DFixture, cubic_interpolate) // All values, current target std::vector result = interpolator.get_values_at_target(); - EXPECT_THAT(result, testing::ElementsAre(testing::DoubleEq(4.416), testing::DoubleEq(8.832))); + EXPECT_THAT(result, testing::ElementsAre(testing::DoubleEq(4.308), testing::DoubleEq(8.616))); } TEST_F(Grid2DFixture, normalize) @@ -389,42 +431,43 @@ TEST_F(Function4DFixture, calculate) TEST_F(Function4DFixture, verify_linear) { - // no matter what we do, result[1] should always be 11! + const double epsilon = 0.0001; + // no matter what we do, result[1] should always be 11! std::vector result; interpolator.set_target(target); result = interpolator.get_values_at_target(); - EXPECT_DOUBLE_EQ(result[1], 11); + EXPECT_NEAR(result[1], 11, epsilon); interpolator.set_axis_interpolation_method(0, Method::cubic); interpolator.set_target(target); result = interpolator.get_values_at_target(); - EXPECT_DOUBLE_EQ(result[1], 11); + EXPECT_NEAR(result[1], 11, epsilon); interpolator.set_axis_interpolation_method(3, Method::cubic); interpolator.set_target(target); result = interpolator.get_values_at_target(); - EXPECT_DOUBLE_EQ(result[1], 11); + EXPECT_NEAR(result[1], 11, epsilon); interpolator.set_axis_interpolation_method(0, Method::linear); interpolator.set_target(target); result = interpolator.get_values_at_target(); - EXPECT_DOUBLE_EQ(result[1], 11); + EXPECT_NEAR(result[1], 11, epsilon); interpolator.set_axis_interpolation_method(2, Method::cubic); interpolator.set_target(target); result = interpolator.get_values_at_target(); - EXPECT_DOUBLE_EQ(result[1], 11); + EXPECT_NEAR(result[1], 11, epsilon); interpolator.set_axis_interpolation_method(0, Method::cubic); interpolator.set_target(target); result = interpolator.get_values_at_target(); - EXPECT_DOUBLE_EQ(result[1], 11); + EXPECT_NEAR(result[1], 11, epsilon); interpolator.set_axis_interpolation_method(1, Method::cubic); interpolator.set_target(target); result = interpolator.get_values_at_target(); - EXPECT_DOUBLE_EQ(result[1], 11); + EXPECT_NEAR(result[1], 11, epsilon); } TEST_F(Function4DFixture, timer) @@ -501,4 +544,4 @@ TEST(CartesianProduct, cartesian_product) EXPECT_THAT(result[3 * 2 * 4 - 1], testing::ElementsAre(3, 5, 9)); } -} // namespace Btwxt \ No newline at end of file +} // namespace Btwxt diff --git a/test/grid-axis-tests.cpp b/test/grid-axis-tests.cpp index 8703809..ee4b259 100644 --- a/test/grid-axis-tests.cpp +++ b/test/grid-axis-tests.cpp @@ -41,19 +41,44 @@ TEST(GridAxis, sorting) TEST(GridAxis, calculate_cubic_spacing_ratios) { - static constexpr std::size_t floor = 0; - static constexpr std::size_t ceiling = 1; - - GridAxis grid_axis({6., 10., 15., 20., 22.}, + static constexpr std::size_t i_interval = 2; + + std::vector grid_values={6., 10., 15., 20., 22.}; + GridAxis grid_axis(grid_values, "", Method::cubic, Method::constant, {-DBL_MAX, DBL_MAX}, std::make_shared()); - EXPECT_THAT(grid_axis.get_cubic_spacing_ratios(floor), - testing::ElementsAre(1, 5.0 / 9, 0.5, 2.0 / 7)); - EXPECT_THAT(grid_axis.get_cubic_spacing_ratios(ceiling), - testing::ElementsAre(4.0 / 9, 0.5, 5.0 / 7, 1)); + + std::vector> result(4, {0., 0.}); + + if ((i_interval > 0) && (i_interval + 1 < grid_values.size())) + { + double w_m1 = grid_values[i_interval] - grid_values[i_interval - 1]; + double w_0 = grid_values[i_interval + 1] - grid_values[i_interval]; + result[0].first = -w_0 * w_0 / w_m1 / (w_0 + w_m1); + result[1].first = (w_0 - w_m1) / w_m1; + result[2].first = w_m1 / (w_0 + w_m1); + } + + if (i_interval + 2 < grid_values.size()) + { + double w_0 = grid_values[i_interval + 1] - grid_values[i_interval]; + double w_1 = grid_values[i_interval + 2] - grid_values[i_interval + 1]; + result[1].second = w_1 / (w_0 + w_1); + result[2].second = -(w_1 - w_0) / w_1; + result[3].second = -w_0 * w_0 / w_1 / (w_0 + w_1); + } + + auto& expected=grid_axis.get_cubic_spacing_ratios(i_interval); + + for (std::size_t i = 0; i<4; ++i) + { + EXPECT_NEAR(expected[i].first, result[i].first, 0.0001); + EXPECT_NEAR(expected[i].second, result[i].second, 0.0001); + } + //EXPECT_THAT(expected[0].first, testing::ElementsAre(result[0].first)); } TEST(GridAxis, bad_limits) @@ -74,4 +99,4 @@ TEST(GridAxis, bad_limits) EXPECT_STDOUT(my_grid_axis.set_extrapolation_limits(extrapolation_limits);, expected_out) EXPECT_EQ(my_grid_axis.get_extrapolation_limits().second, 15); } -} // namespace Btwxt \ No newline at end of file +} // namespace Btwxt diff --git a/test/implementation-tests.cpp b/test/implementation-tests.cpp index 5905d0e..a70ef2a 100644 --- a/test/implementation-tests.cpp +++ b/test/implementation-tests.cpp @@ -16,23 +16,25 @@ namespace Btwxt { TEST_F(CubicImplementationFixture, spacing_multiplier) { - static constexpr std::size_t floor = 0; - static constexpr std::size_t ceiling = 1; + static constexpr std::size_t axis_index = 0; + static constexpr std::size_t elem_index = 1; + const double epsilon = 0.0001; + double result; - result = interpolator.get_axis_cubic_spacing_ratios(0, floor)[0]; - EXPECT_DOUBLE_EQ(result, 1.0); + result = interpolator.get_axis_cubic_spacing_ratios(axis_index, elem_index)[0].first; + EXPECT_NEAR(result, -0.694444, epsilon); - result = interpolator.get_axis_cubic_spacing_ratios(0, ceiling)[0]; - EXPECT_DOUBLE_EQ(result, (10. - 6.) / (15.0 - 6.0)); + result = interpolator.get_axis_cubic_spacing_ratios(axis_index, elem_index)[1].first; + EXPECT_NEAR(result, 0.25, epsilon); - result = interpolator.get_axis_cubic_spacing_ratios(0, floor)[1]; - EXPECT_DOUBLE_EQ(result, (15. - 10.) / (15.0 - 6.0)); + result = interpolator.get_axis_cubic_spacing_ratios(axis_index, elem_index)[1].second; + EXPECT_NEAR(result, 0.5, epsilon); - result = interpolator.get_axis_cubic_spacing_ratios(0, ceiling)[2]; - EXPECT_DOUBLE_EQ(result, 1.0); + result = interpolator.get_axis_cubic_spacing_ratios(axis_index, elem_index)[2].first; + EXPECT_NEAR(result, 0.444444, epsilon); - result = interpolator.get_axis_cubic_spacing_ratios(1, floor)[0]; - EXPECT_DOUBLE_EQ(result, 1.0); + result = interpolator.get_axis_cubic_spacing_ratios(axis_index, elem_index)[3].second; + EXPECT_NEAR(result, -0.5, epsilon); } TEST_F(CubicImplementationFixture, switch_interp_method) @@ -59,29 +61,30 @@ TEST_F(CubicImplementationFixture, interpolate) BtwxtLogger message_display; message_display.info( fmt::format("Time to do cubic interpolation: {} microseconds", duration.count())); - EXPECT_THAT(result, testing::ElementsAre(testing::DoubleEq(4.158), testing::DoubleEq(11.836))); + EXPECT_THAT(result, testing::ElementsAre(testing::DoubleEq(4.1475), testing::DoubleEq(11.82))); } TEST_F(CubicImplementationFixture, grid_point_interp_coeffs) { - static constexpr std::size_t floor = 0; - static constexpr std::size_t ceiling = 1; + static constexpr std::size_t axis_index = 0; + static constexpr std::size_t rel_ctrl_index = 0; + const double epsilon = 0.0001; interpolator.set_target(target); double mu = interpolator.get_floor_to_ceiling_fractions()[0]; std::size_t floor_grid_point_index = interpolator.get_floor_grid_point_coordinates()[0]; - EXPECT_EQ(interpolator.get_interpolation_coefficients()[0][0], - 2 * mu * mu * mu - 3 * mu * mu + 1); - EXPECT_EQ(interpolator.get_interpolation_coefficients()[0][1], -2 * mu * mu * mu + 3 * mu * mu); + EXPECT_NEAR(interpolator.get_interpolation_coefficients()[axis_index][0], + 2 * mu * mu * mu - 3 * mu * mu + 1, epsilon); + EXPECT_NEAR(interpolator.get_interpolation_coefficients()[axis_index][1], + -2 * mu * mu * mu + 3 * mu * mu, epsilon); - EXPECT_EQ(interpolator.get_cubic_slope_coefficients()[0][0], - (mu * mu * mu - 2 * mu * mu + mu) * - interpolator.get_axis_cubic_spacing_ratios(0, floor)[floor_grid_point_index]); - EXPECT_EQ(interpolator.get_cubic_slope_coefficients()[0][1], - (mu * mu * mu - mu * mu) * - interpolator.get_axis_cubic_spacing_ratios(0, ceiling)[floor_grid_point_index]); + EXPECT_NEAR(interpolator.get_cubic_slope_coefficients()[axis_index][rel_ctrl_index], + (1 - mu) * mu * + ((1 - mu) * interpolator.get_axis_cubic_spacing_ratios(axis_index, floor_grid_point_index)[rel_ctrl_index].first + + mu * interpolator.get_axis_cubic_spacing_ratios(axis_index, floor_grid_point_index)[rel_ctrl_index].second), + epsilon); } TEST_F(CubicImplementationFixture, hypercube_weigh_one_vertex) @@ -89,46 +92,85 @@ TEST_F(CubicImplementationFixture, hypercube_weigh_one_vertex) interpolator.set_axis_interpolation_method(1, Method::cubic); interpolator.set_target(target); std::vector methods = interpolator.get_current_methods(); - std::vector mus = interpolator.get_floor_to_ceiling_fractions(); - double mx = mus[0]; - double my = mus[1]; - double c0x = 2 * mx * mx * mx - 3 * mx * mx + 1; - double c0y = 2 * my * my * my - 3 * my * my + 1; - // double c1x = -2*mx*mx*mx + 3*mx*mx; - double c1y = -2 * my * my * my + 3 * my * my; - double d0x = mx * mx * mx - 2 * mx * mx + mx; - double d0y = my * my * my - 2 * my * my + my; - double d1x = mx * mx * mx - mx * mx; - double d1y = my * my * my - my * my; - double s1x = 5.0 / 10; - double s1y = 2.0 / 4; - double s0x = 5.0 / 9; - double s0y = 2.0 / 4; - - std::vector this_vertex = {0, 0}; - double weight = interpolator.get_grid_point_weighting_factor(this_vertex); - double expected_result = c0x * c0y; - expected_result += -1 * c0x * d1y * s1y; - expected_result += -1 * d1x * s1x * c0y; - expected_result += d1x * s1x * d1y * s1y; + + // get interpolation factors + std::vector> linear_interpolation_factors; + std::vector> cubic_interpolation_factors; + for (std::size_t i=0; i < 2; ++i) + { + linear_interpolation_factors.push_back({0.0, 1.0 - mus[i], mus[i], 0.0}); + + double coef = 1 - 2 * mus[i]; + cubic_interpolation_factors.push_back({0.0, coef, -coef, 0.0}); + } + + // get cubic spacing ratios + std::vector>> ratios(2, {{-1.0, 0.0}, {0.0, 1.0}, {1.0, 0.0}, {0.0, -1.0}}); + for (std::size_t i=0; i < 2; ++i) + { + const std::vector& values = interpolator.get_grid_axis(i).get_values(); + std::size_t floor_index = interpolator.get_floor_grid_point_coordinates()[i]; + + double w_m1 = values[floor_index] - values[floor_index - 1]; + double w_0 = values[floor_index + 1] - values[floor_index]; + double w_1 = values[floor_index + 2] - values[floor_index + 1]; + { + double s_0 = (w_0 - w_m1) / w_m1; + double s_1 = w_m1 / (w_0 + w_m1); + double s_m1 = -(s_0 + s_1); + + ratios[i][0].first = s_m1; + ratios[i][1].first = s_0; + ratios[i][2].first = s_1; + } + { + double s_0 = w_1 / (w_0 + w_1); + double s_1 = -(w_1 - w_0) / w_1; + double s_2 = -(s_0 + s_1); + + ratios[i][1].second = s_0; + ratios[i][2].second = s_1; + ratios[i][3].second = s_2; + } + } + + // get slope factors + std::vector> cubic_slope_factors(2, {0.0, 0.0, 0.0, 0.0}); + for (std::size_t i=0; i < 2; ++i) + for (std::size_t j = 0; j < 4; ++j) + { + cubic_slope_factors[i][j] = (1 - mus[i]) * ratios[i][j].first + mus[i] * ratios[i][j].second; + } + + std::vector> weighting_factors(2, {0.0, 0.0, 0.0, 0.0}); + for (std::size_t i=0; i < 2; ++i) + { + for (std::size_t j = 0; j < 4; ++j) + { + weighting_factors[i][j] = linear_interpolation_factors[i][j] + + (1 - mus[i]) * mus[i] * (cubic_interpolation_factors[i][j] + cubic_slope_factors[i][j]); + } + } + + std::vector hypercube_indices = {0, 0}; + double weight = interpolator.get_grid_point_weighting_factor(hypercube_indices); + double expected_result = weighting_factors[0][hypercube_indices[0] + 1] * weighting_factors[1][hypercube_indices[1] + 1]; EXPECT_DOUBLE_EQ(weight, expected_result); - this_vertex = {-1, 1}; - weight = interpolator.get_grid_point_weighting_factor(this_vertex); - expected_result = -1 * d0x * s0x * c1y; - expected_result += -1 * d0x * s0x * d0y * s0y; + hypercube_indices = {-1, 1}; + weight = interpolator.get_grid_point_weighting_factor(hypercube_indices); + expected_result = weighting_factors[0][hypercube_indices[0] + 1] * weighting_factors[1][hypercube_indices[1] + 1]; EXPECT_DOUBLE_EQ(weight, expected_result); - this_vertex = {2, 0}; - weight = interpolator.get_grid_point_weighting_factor(this_vertex); - expected_result = d1x * s1x * c0y; - expected_result += -1 * d1x * s1x * d1y * s1y; + hypercube_indices = {2, 0}; + weight = interpolator.get_grid_point_weighting_factor(hypercube_indices); + expected_result = weighting_factors[0][hypercube_indices[0] + 1] * weighting_factors[1][hypercube_indices[1] + 1]; EXPECT_DOUBLE_EQ(weight, expected_result); - this_vertex = {2, 2}; - weight = interpolator.get_grid_point_weighting_factor(this_vertex); - expected_result = d1x * s1x * d1y * s1y; + hypercube_indices = {2, 2}; + weight = interpolator.get_grid_point_weighting_factor(hypercube_indices); + expected_result = weighting_factors[0][hypercube_indices[0] + 1] * weighting_factors[1][hypercube_indices[1] + 1]; EXPECT_DOUBLE_EQ(weight, expected_result); } @@ -138,26 +180,57 @@ TEST_F(CubicImplementationFixture, hypercube_calculations) interpolator.set_target(target); auto result = interpolator.get_results(); - EXPECT_NEAR(result[0], 4.1953, 0.0001); - EXPECT_NEAR(result[1], 11.9271, 0.0001); + EXPECT_NEAR(result[0], 4.1859, 0.0001); + EXPECT_NEAR(result[1], 11.9070, 0.0001); } TEST_F(CubicImplementationFixture, get_cubic_spacing_ratios) { - // for cubic dimension 0: {6, 10, 15, 20}, multipliers should be: - // floor: {4/4, 5/9, 5/10} - // ceiling: {4/9, 5/10, 5/5} - std::vector> expected_results {{4.0 / 4, 5.0 / 9, 5.0 / 10}, - {4.0 / 9, 5.0 / 10, 5.0 / 5}}; + static constexpr std::size_t axis_index = 0; + const std::vector& axis_values = interpolator.get_grid_axis(axis_index).get_values(); + + double w_m1 = axis_values[1] - axis_values[0]; + double w_0 = axis_values[2] - axis_values[1]; + double w_1 = axis_values[3] - axis_values[2]; + + double sm_0 = (w_0 - w_m1) / w_m1; + double sm_1 = w_m1 / (w_0 + w_m1); + double sm_m1 = -(sm_0 + sm_1); + + double sp_0 = w_1 / (w_0 + w_1); + double sp_1 = -(w_1 - w_0) / w_1; + double sp_2 = -(sp_0 + sp_1); + + static constexpr std::size_t i_ratio = 1; + + std::vector> expected_results {{sm_m1, 0}, {sm_0, sp_0}, {sm_1, sp_1}, {0, sp_2}}; + double result; - for (std::size_t floor_or_ceiling = 0; floor_or_ceiling <= 1; floor_or_ceiling++) { - for (std::size_t index = 0; index < 3; index++) { - result = interpolator.get_axis_cubic_spacing_ratios(0, floor_or_ceiling)[index]; - EXPECT_DOUBLE_EQ(result, expected_results[floor_or_ceiling][index]); - } + for (std::size_t index = 0; index < 4; index++) { + result = interpolator.get_axis_cubic_spacing_ratios(0, i_ratio)[index].first; + EXPECT_DOUBLE_EQ(result, expected_results[index].first); + + result = interpolator.get_axis_cubic_spacing_ratios(0, i_ratio)[index].second; + EXPECT_DOUBLE_EQ(result, expected_results[index].second); } } +TEST_F(CubicImplementationFixture, null_checking_calculations) +{ + std::vector table_with_null = {std::numeric_limits::quiet_NaN(), 3, 1.5, 1, + 5, 4, 2, 1, + 8, 6, 3, 2, + 10, 8, 4, 2}; + + GridPointDataSet dataset_with_null(table_with_null); + interpolator.add_grid_point_data_set(dataset_with_null); + + target = {7, 3}; + interpolator.set_target(target); + auto result = interpolator.get_results(); + EXPECT_TRUE(isnan(result[2])); +} + TEST_F(EmptyGridImplementationFixture, locate_coordinates) { grid = {{1, 2, 3, 4, 5}, {1, 2, 3, 4, 5, 6, 7}}; @@ -357,4 +430,4 @@ TEST_F(Grid3DImplementationFixture, make_linear_hypercube) EXPECT_THAT(hypercube[5], testing::ElementsAre(1, 0, 1)); } -} // namespace Btwxt \ No newline at end of file +} // namespace Btwxt From ad032709075c8cd8a1691f7e805a82a4af341292 Mon Sep 17 00:00:00 2001 From: Phil Ahrenkiel Date: Fri, 23 Jun 2023 14:55:10 -0600 Subject: [PATCH 3/8] Enable selection of slope method --- .gitignore | 1 - README.md | 62 +++++++------- include/btwxt/grid-axis.h | 9 ++- src/grid-axis.cpp | 98 ++++++++++++++++++----- src/regular-grid-interpolator.cpp | 14 ++-- test/btwxt-tests.cpp | 3 +- test/grid-axis-tests.cpp | 129 ++++++++++++++++++++++-------- test/implementation-tests.cpp | 63 +++++++++------ 8 files changed, 264 insertions(+), 115 deletions(-) diff --git a/.gitignore b/.gitignore index d8aa5cb..ce2694d 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,3 @@ build/ .idea/ .vs/ .vscode/ -vendor/ \ No newline at end of file diff --git a/README.md b/README.md index e60365d..5b5928c 100644 --- a/README.md +++ b/README.md @@ -1,35 +1,28 @@ [![Build and Test](https://github.com/bigladder/btwxt/actions/workflows/build-and-test.yml/badge.svg)](https://github.com/bigladder/btwxt/actions/workflows/build-and-test.yml) [![codecov](https://codecov.io/gh/bigladder/btwxt/branch/master/graph/badge.svg)](https://codecov.io/gh/bigladder/btwxt) -# Btwxt +# **btwxt** ## General-purpose, N-dimensional interpolation library... -Btwxt is a free and open source C++ library to perform numerical interpolation on a regular grid data set. The primary -class is a RegularGridInterpolator constructed from: +**btwxt** is a free and open source C++ library to perform multivariate interpolation of one or more functions represented by known values on provided data sets associated with coordinates on a regular (rectilinear) grid. The primary class is a RegularGridInterpolator constructed from: -1. a set of N grid axes representing the values of the independent variables, and -2. a collection of one or more grid point data sets aligned to the grid axes. +1. a set of one or more *grid axes*, representing the unique variates for interpolation, each specfied through a set of *grid-axis coordinates*, and +2. a collection of one or more *data sets*, repesenting the function values at each of the *control points* specified by ordered pairs of the grid-axis coordinates. -A RegularGridInterpolator object can then be queried repeatedly for interpolated values that -fall inside its grid points, or extrapolated values beyond the grid points (up to a limit defined by the user). +A RegularGridInterpolator object can then be queried repeatedly for function values within the space of the grid axes. Interpolation is performed for when a grid coordinate is inside the range of a particular grid axis; extrapolation is performed (up to a user-defined limit) when the grid coordinate is outside that range. A general query of a multivariate function may entail any combination of interpolation and extrapolation on the collection of grid axes. -Btwxt supports linear and cubic spline (Catmull-Rom) interpolation methods. Different methods can be specified for each -axis. The API also allows specification of preferred extrapolation methods (constant or linear)--again independently for -each axis--and extrapolation limits. +**btwxt** supports various methods for both interpolation and extrpolation; different methods can be specifed for each axis. Interpolation can be performed using either linear or cubic (Catmull-Rom) spline. Whereas both linear and cubic splines maintain continuity of a function along the associated axis, cubic splines also maintain continuty of the function's first derivative on that axis. Extrapolation can be either constant or linear, with the option to specify finite extrapolation limits for a given axis. -The grid is required to be regular (rectilinear), but the axes are not required to have uniform spacing. Each axis need -only be strictly increasing. +The grid is required to be rectilinear, but the axes are not required to have regular (uniform) spacing. However, the grid coordinates provided for each axis must be strictly increasing. -The grid point data sets are imported as a single long vector (i.e., a flattened array). It is the user's responsibility -to ensure that the grid point data is ordered according to the order of grid axes. Btwxt ensures that the number of -grid points aligns with the total combined points of the grid axes. +Data sets are imported as a single long vector (i.e., a flattened array). The user must ensure that values in a data set are ordered, with respect to the orders of both grid coordinates and grid axes, such that each value is correctly associated with the appropriate control point. **btwxt** verifies that the number of values with each data set corresponds with the total number of control points. ## How to Use ### API -An individual axis can be instantiated with: +The snippet below illustrates instantiation of an individual axis. ```c++ std::vector one_axis{6, 10, 15}; @@ -43,33 +36,42 @@ std::pair extrapolation_limits{0, 20}; GridAxis(one_axis, "x", interpolation_method, extrapolation_method, extrapolation_limits); ``` -A RegularGridInterpolator holds a collection of axes and the corresponding grid point data at all permutations of the grid axes' values. +A RegularGridInterpolator holds a collection of axes, which specify the control points, and a collection of data sets containing function values for each of the control points. The snippet below illustrates instantiation of a two-dimensional grid, an associated data set, and their assocation with a RegularGridInterpolator. ```c++ std::vector> vector_grid = {{6, 10, 15}, {4, 6, 9}}; -std::vector> values = {...}; +std::vector> values = {2, 4, 7, 1, 3, 6}; RegularGridInterpolator my_interpolator(vector_grid, values); ``` -Once you have an RegularGridInterpolator object, you can set a point to interpolate to: +Once a RegularGridInterpolator has been synthesized, the interpolated can be queried for function values at target grid-axis coordinates. Typical syntax to perform such a query on a 2-D grid is shown in the snippet below. ```c++ +// 1. set the target before requesting values: std::vector target{12.5, 5.1}; my_interpolator.set_target(target); std::vector result = my_interpolator.get_values_at_target(); -// or set the target when requesting values: -std::vector new_target{11.7, 6.1}; -result = my_interpolator.get_values_at_target(new_target); +// 2. set the target when requesting values: +std::vector target{12.5, 5.1}; +std::vector result = my_interpolator.get_values_at_target(target); + +// 3. set the target by applying the () operator: +std::vector result = my_interpolator({12.5, 5.1}); ``` +### Single-Axis Interpolation +#### Linear case +The *linear* interpolation of a function $f(x)$ on a single axis within the interval between control points $P_{0}$ and $P_{1}$, with axis coordinates $x_{0}$ and $x_{1}$, respectively, can be expressed as +$$h(\mu) = (1-\mu)\cdot f_{0}+\mu \cdot f_{1}$$ -We have overloaded the () operator on the RegularGridInterpolator to perform the interpolation, enabling the following: +where +$$\mu = {x-x_{0}\over x_{1}-x_{0}}$$ -```c++ -std::vector target{12.5, 5.1}; -my_interpolator.set_target(target); -std::vector result = my_interpolator(); -std::vector new_target{11.7, 6.1}; -result = my_interpolator(new_target); -``` \ No newline at end of file +and $f_{0}=f(x_{0})$ and $f_{1}=f(x_{1})$. + +#### Cubic case +We can also ensure continuity of the function slope through the use of *cubic* interpolation, by incorporating the neighboring control points, $P_{-1}$, and $P_{2}$, exterior to the interval in either direction. Various methods have been considered for specifying the slopes at $P_{0}$ and $P_{1}$. +The three points $P_{-1}$, $P_{0}$, and $P_{1}$ specify a quadractic curve $g_{0}(\mu)$; $P_{0}$, $P_{1}$, and $P_{2 }$ specify another quadractic curve $g_{1}(\mu)$. Note that these two functions both contain the points $P_{0}$ and $P_{1}$ that define our interval of interest, so the linear interpolation +$$h(\mu) = (1-\mu)\cdot g_{0}(\mu)+\mu \cdot g_{1}(\mu)$$ + will also contain these points. Further analysis reveals that the first derivatives at these points are equal to those of $g_{0}$ and $g_{1}$, respectively. Continuity of both the value and first derivative are therefore maintained. \ No newline at end of file diff --git a/include/btwxt/grid-axis.h b/include/btwxt/grid-axis.h index 08e6160..2baff76 100644 --- a/include/btwxt/grid-axis.h +++ b/include/btwxt/grid-axis.h @@ -18,6 +18,7 @@ namespace Btwxt { enum class Method { undefined, constant, linear, cubic }; +enum class SlopeMethod { undefined, finite_diff, cardinal_0, cardinal_1, quadratic }; class GridAxis { // A single input dimension of the grid @@ -31,6 +32,7 @@ class GridAxis { Method interpolation_method = Method::linear, Method extrapolation_method = Method::constant, std::pair extrapolation_limits = {-DBL_MAX, DBL_MAX}, + SlopeMethod slope_method = SlopeMethod::quadratic, const std::shared_ptr& logger = std::make_shared()); // Setters @@ -41,6 +43,7 @@ class GridAxis { extrapolation_limits = limits; check_extrapolation_limits(); } + void set_slope_method(SlopeMethod slope_method_in); void set_logger(std::shared_ptr logger_in) { @@ -67,8 +70,12 @@ class GridAxis { std::vector values; Method extrapolation_method {Method::constant}; Method interpolation_method {Method::linear}; + SlopeMethod slope_method {SlopeMethod::quadratic}; std::pair extrapolation_limits {-DBL_MAX, DBL_MAX}; - std::vector>> cubic_spacing_ratios; + std::vector>> + cubic_spacing_ratios; // Used for cubic interpolation. Outer vector is length + // of axis values; inner vector is size 4 (points in relative positions + // -1, 0, 1, 2). Pair elements are applied to (-) and (+) coefficients. std::shared_ptr logger; void calculate_cubic_spacing_ratios(); void check_grid_sorted(); diff --git a/src/grid-axis.cpp b/src/grid-axis.cpp index 844be95..0a1c796 100644 --- a/src/grid-axis.cpp +++ b/src/grid-axis.cpp @@ -11,24 +11,27 @@ GridAxis::GridAxis(std::vector values_in, Method interpolation_method, Method extrapolation_method, std::pair extrapolation_limits, + SlopeMethod slope_method, const std::shared_ptr& logger_in) : name(name) , values(std::move(values_in)) , extrapolation_method(extrapolation_method) , interpolation_method(interpolation_method) + , slope_method(slope_method) , extrapolation_limits(std::move(extrapolation_limits)) , cubic_spacing_ratios( - std::max(static_cast(values.size()) - 1, 0), - {{-1.0, 0.0}, {0.0, 1.0}, {1.0, 0.0}, {0.0, -1.0}}) - , logger(logger_in) -{ + std::max(static_cast(values.size()) - 1, 0), + { { -1.0, 0.0 }, { 0.0, -1.0 }, { 1.0, 0.0 }, { 0.0, 1.0 } }) + , logger(logger_in) { if (values.empty()) { throw BtwxtException( - fmt::format("Cannot create grid axis (name=\"{}\") from a zero-length vector.", name), - *logger); + fmt::format("Cannot create grid axis (name=\"{}\") from a zero-length vector.", name), + *logger); } + check_grid_sorted(); check_extrapolation_limits(); + if (interpolation_method == Method::cubic) { calculate_cubic_spacing_ratios(); } @@ -75,6 +78,11 @@ void GridAxis::set_extrapolation_method(Method extrapolation_method_in) extrapolation_method = extrapolation_method_in; } +void GridAxis::set_slope_method(SlopeMethod slope_method_in) +{ + slope_method = slope_method_in; +} + void GridAxis::calculate_cubic_spacing_ratios() { if (get_length() == 1) { @@ -89,28 +97,78 @@ void GridAxis::calculate_cubic_spacing_ratios() } for (std::size_t i_elem = 0; i_elem < values.size() - 1; i_elem++) { + auto& ratio = cubic_spacing_ratios[i_elem]; double w_0 = values[i_elem + 1] - values[i_elem]; if (i_elem > 0) { double w_m1 = values[i_elem] - values[i_elem - 1]; - double s0_m = (w_0 - w_m1) / w_m1; - double s1_m = w_m1 / (w_0 + w_m1); - - cubic_spacing_ratios[i_elem][0].first = -(s0_m + s1_m); - cubic_spacing_ratios[i_elem][1].first = s0_m; - cubic_spacing_ratios[i_elem][2].first = s1_m; - cubic_spacing_ratios[i_elem][3].first = 0.0; - } + double t_0 = w_0 / w_m1; + + double c_0(0.0); + switch (slope_method) + { + case SlopeMethod::finite_diff:{ + c_0 = 0.5; + break; + } + case SlopeMethod::cardinal_0:{ + c_0 = 1 / (1 + t_0); + break; + } + case SlopeMethod::cardinal_1:{ + c_0 = 0.0; + break; + } + case SlopeMethod::quadratic: + default:{ + c_0 = t_0 / (1 + t_0); + break; + } + } + + //general + double s_m1_m = -t_0 * c_0;; + double s_1_m = 1 - c_0; + + ratio[0].first = s_m1_m; + ratio[1].first = -(s_m1_m + s_1_m); + ratio[2].first = s_1_m; + ratio[3].first = 0.0; + } if (i_elem + 2 < values.size()) { double w_1 = values[i_elem + 2] - values[i_elem + 1]; - double s0_p = w_1 / (w_0 + w_1); - double s1_p = (w_0 - w_1) / w_1; + double t_1 = w_0 / w_1; + + double c_1(0.0); + switch (slope_method){ + case SlopeMethod::finite_diff: + { + c_1 = 0.5; + break; + } + case SlopeMethod::cardinal_0:{ + c_1 = 1 / (1 + t_1); + break; + } + case SlopeMethod::cardinal_1:{ + c_1 = 0.0; + break; + } + case SlopeMethod::quadratic: + default:{ + c_1 = t_1 / (1 + t_1); + break; + } + } + + double s_0_p = -(1 - c_1); + double s_2_p = t_1 * c_1; - cubic_spacing_ratios[i_elem][0].second = 0.0; - cubic_spacing_ratios[i_elem][1].second = s0_p; - cubic_spacing_ratios[i_elem][2].second = s1_p; - cubic_spacing_ratios[i_elem][3].second = -(s0_p + s1_p); + ratio[0].second = 0.0; + ratio[1].second = s_0_p; + ratio[2].second = -(s_0_p + s_2_p); + ratio[3].second = s_2_p; } } } diff --git a/src/regular-grid-interpolator.cpp b/src/regular-grid-interpolator.cpp index b79f59d..1f0d4b8 100644 --- a/src/regular-grid-interpolator.cpp +++ b/src/regular-grid-interpolator.cpp @@ -25,6 +25,7 @@ std::vector construct_grid_axes(const std::vector> Method::linear, Method::constant, std::pair {-DBL_MAX, DBL_MAX}, + SlopeMethod::quadratic, logger_in); } return grid_axes; @@ -151,16 +152,18 @@ RegularGridInterpolator& RegularGridInterpolator::operator=(const RegularGridInt return *this; } -void RegularGridInterpolatorImplementation::set_axis_sizes() -{ +void RegularGridInterpolatorImplementation::set_axis_sizes() { number_of_grid_points = 1; + for (std::size_t axis_index = number_of_grid_axes; axis_index-- > 0;) { std::size_t length = grid_axes[axis_index].get_length(); + if (length == 0) { throw BtwxtException( - fmt::format("Grid axis (name=\"{}\") has zero length.", grid_axes[axis_index].name), - *logger); + fmt::format("Grid axis (name=\"{}\") has zero length.", grid_axes[axis_index].name), + *logger); } + grid_axis_lengths[axis_index] = length; grid_axis_step_size[axis_index] = number_of_grid_points; number_of_grid_points *= length; @@ -168,7 +171,6 @@ void RegularGridInterpolatorImplementation::set_axis_sizes() } // Public manipulation methods - std::size_t RegularGridInterpolator::add_grid_point_data_set(const std::vector& grid_point_data_vector, const std::string& name) @@ -699,7 +701,7 @@ void RegularGridInterpolatorImplementation::calculate_interpolation_coefficients auto& ratios = get_axis_cubic_spacing_ratios(axis_index, floor_grid_point_coordinates[axis_index]); for (std::size_t i = 0; i < 4; ++i) { - cubic_slope_factors[i] = (1 - mu) * ratios[i].first + mu * ratios[i].second; + cubic_slope_factors[i] = (1 - mu) * ratios[i].first - mu * ratios[i].second; cubic_slope_coefficients[axis_index][i] = (1 - mu) * mu * cubic_slope_factors[i]; } for (std::size_t i = 0; i < 4; ++i) diff --git a/test/btwxt-tests.cpp b/test/btwxt-tests.cpp index d70aece..da21883 100644 --- a/test/btwxt-tests.cpp +++ b/test/btwxt-tests.cpp @@ -35,10 +35,11 @@ TEST(BasicSplineTest, cubic_interp) Method extrapolation_method = Method::linear; Method interpolation_method = Method::cubic; + SlopeMethod slope_method = SlopeMethod::quadratic; std::pair extrapolation_limits{-DBL_MAX, DBL_MAX}; - GridAxis axis(x_values, "x", interpolation_method, extrapolation_method, extrapolation_limits); + GridAxis axis(x_values, "x", interpolation_method, extrapolation_method, extrapolation_limits, slope_method); RegularGridInterpolator my_interpolator({axis}, values); const double epsilon = 0.0001; diff --git a/test/grid-axis-tests.cpp b/test/grid-axis-tests.cpp index ee4b259..29d1d81 100644 --- a/test/grid-axis-tests.cpp +++ b/test/grid-axis-tests.cpp @@ -44,41 +44,106 @@ TEST(GridAxis, calculate_cubic_spacing_ratios) static constexpr std::size_t i_interval = 2; std::vector grid_values={6., 10., 15., 20., 22.}; - GridAxis grid_axis(grid_values, - "", - Method::cubic, - Method::constant, - {-DBL_MAX, DBL_MAX}, - std::make_shared()); - - std::vector> result(4, {0., 0.}); - - if ((i_interval > 0) && (i_interval + 1 < grid_values.size())) + std::vector slope_method_vector + ({ SlopeMethod::quadratic, + SlopeMethod::finite_diff, + SlopeMethod::cardinal_0, + SlopeMethod::cardinal_1 + + }); + for (auto slope_method: slope_method_vector) { - double w_m1 = grid_values[i_interval] - grid_values[i_interval - 1]; + GridAxis grid_axis(grid_values, + "", + Method::cubic, + Method::constant, + {-DBL_MAX, DBL_MAX}, + slope_method, + std::make_shared()); + + std::vector> result(4, {0., 0.}); + double w_0 = grid_values[i_interval + 1] - grid_values[i_interval]; - result[0].first = -w_0 * w_0 / w_m1 / (w_0 + w_m1); - result[1].first = (w_0 - w_m1) / w_m1; - result[2].first = w_m1 / (w_0 + w_m1); - } - - if (i_interval + 2 < grid_values.size()) - { - double w_0 = grid_values[i_interval + 1] - grid_values[i_interval]; - double w_1 = grid_values[i_interval + 2] - grid_values[i_interval + 1]; - result[1].second = w_1 / (w_0 + w_1); - result[2].second = -(w_1 - w_0) / w_1; - result[3].second = -w_0 * w_0 / w_1 / (w_0 + w_1); + if ((i_interval > 0) && (i_interval + 1 < grid_values.size())){ + double w_m1 = grid_values[i_interval] - grid_values[i_interval - 1]; + double t_0 = w_0 / w_m1; + + double c_0(0.0); + switch (slope_method){ + case SlopeMethod::finite_diff: + { + c_0 = 0.5; + break; + } + + case SlopeMethod::cardinal_0:{ + c_0 = 1 / (1 + t_0); + break; + } + + case SlopeMethod::cardinal_1:{ + c_0 = 0.0; + break; + } + + case SlopeMethod::quadratic: + default:{ + c_0 = t_0 / (1 + t_0); + break; + } + } + + //general + double s_m1_m = -t_0 * c_0;; + double s_1_m = 1 - c_0; + result[0].first = s_m1_m; + result[1].first = -(s_m1_m + s_1_m); + result[2].first = s_1_m; + } + + if (i_interval + 2 < grid_values.size()){ + double w_1 = grid_values[i_interval + 2] - grid_values[i_interval + 1]; + double t_1 = w_0 / w_1; + + double c_1(0.0); + switch (slope_method){ + case SlopeMethod::finite_diff:{ + c_1 = 0.5; + break; + } + + case SlopeMethod::cardinal_0:{ + c_1 = 1 / (1 + t_1);; + break; + } + + case SlopeMethod::cardinal_1:{ + c_1 = 0.0; + break; + } + + case SlopeMethod::quadratic: + default:{ + c_1 = t_1 / (1 + t_1); + break; + } + } + + double s_0_p = -(1 - c_1); + double s_2_p = t_1 * c_1; + + result[1].second = s_0_p; + result[2].second = -(s_0_p + s_2_p); + result[3].second = s_2_p; + } + + auto& expected=grid_axis.get_cubic_spacing_ratios(i_interval); + + for (std::size_t i = 0; i < 4; ++i){ + EXPECT_NEAR(expected[i].first, result[i].first, 0.0001); + EXPECT_NEAR(expected[i].second, result[i].second, 0.0001); + } } - - auto& expected=grid_axis.get_cubic_spacing_ratios(i_interval); - - for (std::size_t i = 0; i<4; ++i) - { - EXPECT_NEAR(expected[i].first, result[i].first, 0.0001); - EXPECT_NEAR(expected[i].second, result[i].second, 0.0001); - } - //EXPECT_THAT(expected[0].first, testing::ElementsAre(result[0].first)); } TEST(GridAxis, bad_limits) diff --git a/test/implementation-tests.cpp b/test/implementation-tests.cpp index a70ef2a..aa37b60 100644 --- a/test/implementation-tests.cpp +++ b/test/implementation-tests.cpp @@ -28,13 +28,13 @@ TEST_F(CubicImplementationFixture, spacing_multiplier) EXPECT_NEAR(result, 0.25, epsilon); result = interpolator.get_axis_cubic_spacing_ratios(axis_index, elem_index)[1].second; - EXPECT_NEAR(result, 0.5, epsilon); + EXPECT_NEAR(result, -0.5, epsilon); result = interpolator.get_axis_cubic_spacing_ratios(axis_index, elem_index)[2].first; EXPECT_NEAR(result, 0.444444, epsilon); result = interpolator.get_axis_cubic_spacing_ratios(axis_index, elem_index)[3].second; - EXPECT_NEAR(result, -0.5, epsilon); + EXPECT_NEAR(result, 0.5, epsilon); } TEST_F(CubicImplementationFixture, switch_interp_method) @@ -106,31 +106,40 @@ TEST_F(CubicImplementationFixture, hypercube_weigh_one_vertex) } // get cubic spacing ratios - std::vector>> ratios(2, {{-1.0, 0.0}, {0.0, 1.0}, {1.0, 0.0}, {0.0, -1.0}}); - for (std::size_t i=0; i < 2; ++i) + std::vector>> ratios(2, {{-1.0, 0.0}, {0.0, -1.0}, {1.0, 0.0}, {0.0, 1.0}}); + for (std::size_t i = 0; i < 2; ++i) { const std::vector& values = interpolator.get_grid_axis(i).get_values(); std::size_t floor_index = interpolator.get_floor_grid_point_coordinates()[i]; - double w_m1 = values[floor_index] - values[floor_index - 1]; double w_0 = values[floor_index + 1] - values[floor_index]; - double w_1 = values[floor_index + 2] - values[floor_index + 1]; { - double s_0 = (w_0 - w_m1) / w_m1; - double s_1 = w_m1 / (w_0 + w_m1); - double s_m1 = -(s_0 + s_1); + double w_m1 = values[floor_index] - values[floor_index - 1]; + double t_0 = w_0 / w_m1; + + double c_0(0.5); + c_0 = t_0 / (1 + t_0); // quadratic + //general + double s_m1 = -t_0 * c_0; + double s_1 = 1 - c_0; + ratios[i][0].first = s_m1; - ratios[i][1].first = s_0; + ratios[i][1].first = -(s_m1 + s_1); ratios[i][2].first = s_1; } { - double s_0 = w_1 / (w_0 + w_1); - double s_1 = -(w_1 - w_0) / w_1; - double s_2 = -(s_0 + s_1); + double w_1 = values[floor_index + 2] - values[floor_index + 1]; + double t_1 = w_0 / w_1; + + double c_1(0.5); + c_1 = t_1 / (1 + t_1); // quadratic + + double s_0 = -(1 - c_1); + double s_2 = t_1 * c_1; ratios[i][1].second = s_0; - ratios[i][2].second = s_1; + ratios[i][2].second = -(s_0 + s_2); ratios[i][3].second = s_2; } } @@ -140,7 +149,7 @@ TEST_F(CubicImplementationFixture, hypercube_weigh_one_vertex) for (std::size_t i=0; i < 2; ++i) for (std::size_t j = 0; j < 4; ++j) { - cubic_slope_factors[i][j] = (1 - mus[i]) * ratios[i][j].first + mus[i] * ratios[i][j].second; + cubic_slope_factors[i][j] = (1 - mus[i]) * ratios[i][j].first - mus[i] * ratios[i][j].second; } std::vector> weighting_factors(2, {0.0, 0.0, 0.0, 0.0}); @@ -193,17 +202,23 @@ TEST_F(CubicImplementationFixture, get_cubic_spacing_ratios) double w_0 = axis_values[2] - axis_values[1]; double w_1 = axis_values[3] - axis_values[2]; - double sm_0 = (w_0 - w_m1) / w_m1; - double sm_1 = w_m1 / (w_0 + w_m1); - double sm_m1 = -(sm_0 + sm_1); - - double sp_0 = w_1 / (w_0 + w_1); - double sp_1 = -(w_1 - w_0) / w_1; - double sp_2 = -(sp_0 + sp_1); + double t_0 = w_0 / w_m1; + double t_1 = w_0 / w_1; + + double c_0(0.5), c_1(0.5); + c_0 = t_0 / (1 + t_0); // quadratic + c_1 = t_1 / (1 + t_1); // quadratic + + //general + double s_m1_m = -t_0 * c_0; + double s_1_m = 1 - c_0; + + double s_0_p = -(1 - c_1); + double s_2_p = t_1 * c_1; static constexpr std::size_t i_ratio = 1; - std::vector> expected_results {{sm_m1, 0}, {sm_0, sp_0}, {sm_1, sp_1}, {0, sp_2}}; + std::vector> expected_results {{s_m1_m, 0}, {-(s_m1_m + s_1_m), s_0_p}, {s_1_m, -(s_0_p + s_2_p)}, {0, s_2_p}}; double result; for (std::size_t index = 0; index < 4; index++) { @@ -343,7 +358,7 @@ TEST_F(Grid2DImplementationFixture, construct_from_axes) GridAxis ax0 = GridAxis({0, 10, 15}, "ax0"); auto logger = ax0.get_logger(); GridAxis ax1 = - GridAxis({4, 6}, "ax1", Method::linear, Method::constant, {-DBL_MAX, DBL_MAX}, logger); + GridAxis({4, 6}, "ax1", Method::linear, Method::constant, {-DBL_MAX, DBL_MAX}, SlopeMethod::quadratic, logger); std::vector test_axes = {ax0, ax1}; interpolator = RegularGridInterpolatorImplementation(test_axes, logger); EXPECT_EQ(interpolator.get_number_of_grid_axes(), 2u); From bb1f1e2a893aef0cc6d3415f41a4bda1fddbd431 Mon Sep 17 00:00:00 2001 From: Phil Ahrenkiel Date: Fri, 23 Jun 2023 18:08:44 -0600 Subject: [PATCH 4/8] Update README.md --- README.md | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 5b5928c..aa05449 100644 --- a/README.md +++ b/README.md @@ -62,16 +62,20 @@ std::vector result = my_interpolator({12.5, 5.1}); ``` ### Single-Axis Interpolation #### Linear case -The *linear* interpolation of a function $f(x)$ on a single axis within the interval between control points $P_{0}$ and $P_{1}$, with axis coordinates $x_{0}$ and $x_{1}$, respectively, can be expressed as -$$h(\mu) = (1-\mu)\cdot f_{0}+\mu \cdot f_{1}$$ +A *linear* interpolation over the interval [$x_{k}$, $x_{k+1}$] can be expressed in terms of the data values $f_{k}$ and $f_{k+1}$ as +$$h_{k}(\mu_{k}) = (1-\mu)\cdot f_{k}+\mu \cdot f_{k+1}$$ where -$$\mu = {x-x_{0}\over x_{1}-x_{0}}$$ +$$\mu_{k} = {x-x_{k}\over w_{k}}$$ -and $f_{0}=f(x_{0})$ and $f_{1}=f(x_{1})$. +using the interval width $w_{k}=x_{k+1}-x_{k}$. The resulting interpolation will consist of straight line segments spanning each interval between control points. #### Cubic case -We can also ensure continuity of the function slope through the use of *cubic* interpolation, by incorporating the neighboring control points, $P_{-1}$, and $P_{2}$, exterior to the interval in either direction. Various methods have been considered for specifying the slopes at $P_{0}$ and $P_{1}$. -The three points $P_{-1}$, $P_{0}$, and $P_{1}$ specify a quadractic curve $g_{0}(\mu)$; $P_{0}$, $P_{1}$, and $P_{2 }$ specify another quadractic curve $g_{1}(\mu)$. Note that these two functions both contain the points $P_{0}$ and $P_{1}$ that define our interval of interest, so the linear interpolation -$$h(\mu) = (1-\mu)\cdot g_{0}(\mu)+\mu \cdot g_{1}(\mu)$$ - will also contain these points. Further analysis reveals that the first derivatives at these points are equal to those of $g_{0}$ and $g_{1}$, respectively. Continuity of both the value and first derivative are therefore maintained. \ No newline at end of file +A *cubic* interpolation generates a smoother funciton that can be constructed to have continuous first derivative (slope), which is useful for many types of data analysis. We apply the formulation used above to a pair of functions $g_{k}(\mu_{k})$ and $g_{k+1}(\mu_{k})$, up to quadratic in order, associated with each control point, rather than the data values, themselves. The interpolation becomes +$$h(\mu_{k}) = (1-\mu_{k})\cdot g_{k}(\mu_{k})+\mu_{k} \cdot g_{k+1}(\mu_{k})$$ + +(Note that the representations of these functions may depend on the argument indicated.) To enforce continuity, these two functions must intersect and have the appropriate data value at both $x_{k}$ and $x_{k+1}$. Although the slope is not uniquely defined for a discrete data set, it can be estimated by inspection of neighboring control points. It is useful to express the slope at $x_{k}$ as +$$f'_{k}=(1-c_{k})\cdot r_{k}+c_{k}\cdot r_{k-1}$$ +where +$$r_{k}=\frac{f_{k+1}-f_{k}}{w_{k}}$$ +In **btwxt**, the user can select from among various methods to specify $c_{k}$, which may be applicable in different contexts. \ No newline at end of file From 872b1da2307716a51ec9ab3bfd52352df22e4a59 Mon Sep 17 00:00:00 2001 From: Phil Ahrenkiel Date: Mon, 26 Jun 2023 09:58:26 -0600 Subject: [PATCH 5/8] added slope_reduction and updated markdown --- README.md | 73 ++++++++++++++++++++++++------- include/btwxt/grid-axis.h | 7 ++- src/grid-axis.cpp | 33 +++++++------- src/regular-grid-interpolator.cpp | 1 + test/btwxt-tests.cpp | 3 +- test/grid-axis-tests.cpp | 50 ++++++++------------- test/implementation-tests.cpp | 29 ++++++------ 7 files changed, 114 insertions(+), 82 deletions(-) diff --git a/README.md b/README.md index aa05449..71a805a 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ A RegularGridInterpolator object can then be queried repeatedly for function values within the space of the grid axes. Interpolation is performed for when a grid coordinate is inside the range of a particular grid axis; extrapolation is performed (up to a user-defined limit) when the grid coordinate is outside that range. A general query of a multivariate function may entail any combination of interpolation and extrapolation on the collection of grid axes. -**btwxt** supports various methods for both interpolation and extrpolation; different methods can be specifed for each axis. Interpolation can be performed using either linear or cubic (Catmull-Rom) spline. Whereas both linear and cubic splines maintain continuity of a function along the associated axis, cubic splines also maintain continuty of the function's first derivative on that axis. Extrapolation can be either constant or linear, with the option to specify finite extrapolation limits for a given axis. +**btwxt** supports various methods for both interpolation and extrpolation; different methods can be specifed for each axis. Interpolation can be performed using either linear or cubic spline. Whereas both linear and cubic splines maintain continuity of a function along the associated axis, cubic splines also maintain continuty of the function's first derivative on that axis. Extrapolation can be either constant or linear, with the option to specify finite extrapolation limits for a given axis. The grid is required to be rectilinear, but the axes are not required to have regular (uniform) spacing. However, the grid coordinates provided for each axis must be strictly increasing. @@ -32,8 +32,11 @@ GridAxis(one_axis, "x"); // Construct a new axis named "x" // extrapolation limits default to {-DBL_MAX, DBL_MAX} Method extrapolation_method = Method::linear; Method interpolation_method = Method::cubic; +SlopeMethod slope_method = SlopeMethod::quadratic; +double slope_reduction = 0.0; std::pair extrapolation_limits{0, 20}; -GridAxis(one_axis, "x", interpolation_method, extrapolation_method, extrapolation_limits); +GridAxis(one_axis, "x", interpolation_method, extrapolation_method, extrapolation_limits, +slope_method, slope_reduction); ``` A RegularGridInterpolator holds a collection of axes, which specify the control points, and a collection of data sets containing function values for each of the control points. The snippet below illustrates instantiation of a two-dimensional grid, an associated data set, and their assocation with a RegularGridInterpolator. @@ -60,22 +63,60 @@ std::vector result = my_interpolator.get_values_at_target(target); // 3. set the target by applying the () operator: std::vector result = my_interpolator({12.5, 5.1}); ``` -### Single-Axis Interpolation +### Single-Axis Spline Interpolation +An interpolation along a single axis generates a function that spans the intervals between control points, typically corresponding to a linear combination of the data values. A spline is essentially a chain of such functions, such that a different set of coefficients is identified for each interval. For example, in the interval [$x_{k}$, $x_{k+1}$] we write +$$h_k(\mu_k) = \sum^{n+1}_{i=-n}\sigma_{i}^{(k)} (\mu_k)\cdot f_{k+i}$$ +where $f_{k}$ is the data value at $x_{k}$ and +$$\mu_{k}={x-x_{k}\over w_{k}}$$ +The interval width is given by $w_{k}=x_{k+1}-x_{k}$. It is useful to limit the range of the summation for practical computational. Without loss of generality, we can focus on the case $k=0$, simplifying the notation: +$$h(\mu) = \sum^{n+1}_{i=-n}\sigma_{i} (\mu)\cdot f_{i}$$ +using $\mu_{0}\rightarrow\mu$ and $w_{0}\rightarrow w$. + +Various orders of differentiability may berequired of the spline. Continuity of the $N^{th}$-order derivative of $h(\mu)$ can be satisfied with polynomials of order $2N+1$, motivating the prominence of linear and cubic splines. #### Linear case -A *linear* interpolation over the interval [$x_{k}$, $x_{k+1}$] can be expressed in terms of the data values $f_{k}$ and $f_{k+1}$ as -$$h_{k}(\mu_{k}) = (1-\mu)\cdot f_{k}+\mu \cdot f_{k+1}$$ - -where -$$\mu_{k} = {x-x_{k}\over w_{k}}$$ - -using the interval width $w_{k}=x_{k+1}-x_{k}$. The resulting interpolation will consist of straight line segments spanning each interval between control points. +A *linear* spline can be generated with $n=0$, such that data values $f_0$ and $f_1$ are reproduced precisely at $x_0$ and $x_1$, respectively: +$$h(\mu)=(1-\mu)\cdot f_{0}+\mu\cdot f_{1}$$ +The coefficients are clearly +$$\sigma_{0}(\mu)=1-\mu, \qquad \sigma_{1}(\mu)=\mu$$ +The resulting interpolation preserve will consist of straight line segments spanning each interval between control points. #### Cubic case -A *cubic* interpolation generates a smoother funciton that can be constructed to have continuous first derivative (slope), which is useful for many types of data analysis. We apply the formulation used above to a pair of functions $g_{k}(\mu_{k})$ and $g_{k+1}(\mu_{k})$, up to quadratic in order, associated with each control point, rather than the data values, themselves. The interpolation becomes -$$h(\mu_{k}) = (1-\mu_{k})\cdot g_{k}(\mu_{k})+\mu_{k} \cdot g_{k+1}(\mu_{k})$$ - -(Note that the representations of these functions may depend on the argument indicated.) To enforce continuity, these two functions must intersect and have the appropriate data value at both $x_{k}$ and $x_{k+1}$. Although the slope is not uniquely defined for a discrete data set, it can be estimated by inspection of neighboring control points. It is useful to express the slope at $x_{k}$ as -$$f'_{k}=(1-c_{k})\cdot r_{k}+c_{k}\cdot r_{k-1}$$ +A *cubic* interpolation generates a smoother function that can be constructed to have continuous first derivative (slope), which is useful for many types of data analysis. The general form can be written +$$h(\mu)=a_0(\mu)\cdot f_0+a_1(\mu)\cdot f_1+b_0(\mu)\cdot w_0\cdot f'_0++b_1(\mu)\cdot w_0\cdot f'_1$$ +with +$$\begin{align*}a_0(\mu)&=(1-\mu)\cdot [1+\mu\cdot (1-2\mu)]\\ +a_1(\mu)&=\mu\cdot[1-(1-\mu)\cdot(1-2\mu)]\\ +b_0(\mu)&=(1-\mu)^2\cdot \mu\\ +b_1(\mu)&=-(1-\mu)\cdot \mu^2\\ +\end{align*}$$ +Nonetheless, the slope is not uniquely defined for a discrete data set, so a reasonable estimation must be generated by inspection of neighboring control points. +In fact, we need consider only the additional data values at $x_{k-1}$ and $x_{k+1}$ to develop an expression for $f'_{k}$: +$$f'_{k}=(1-\alpha)\cdot [(1-\beta_{k})\cdot r_{k}+\beta_{k}\cdot r_{k-1}]$$ where $$r_{k}=\frac{f_{k+1}-f_{k}}{w_{k}}$$ -In **btwxt**, the user can select from among various methods to specify $c_{k}$, which may be applicable in different contexts. \ No newline at end of file +The factor $1-\alpha$ uniformly scales the slope at control points. A default value of $\alpha =0$ is assigned. Various criteria for the specification of $\beta _{k}$ may be applicable in different contexts; **btwxt** allows provides options for slope evaluation. The coefficient may depend on the relative widths of neighboring intervals. We define +$$t_{k}=\frac{w_{k}}{w_{k-1}}$$ +The default setting computes the slope by a local fit to a *quadratic* function about each control point, which corresponds to +$$\beta_{k}=\frac{t_{k}}{1 + t_{k}}$$ +The *cardinal* method is realized with +$$\beta_{k}=\frac{1}{1 + t_{k}}$$ +The slope in *finite-difference* methods is computed with $\beta _{k}=0.5$\, imposing no width dependence. + +We can now represent the cubic spline over [$x_{0}$, $x_{1}$] as the sum +$$h(\mu)=\sigma_{-1}(\mu)\cdot f_{-1}+\sigma_{0}(\mu)\cdot f_{0}+\sigma_{1}(\mu)\cdot f_{1}+\sigma_{2}(\mu)\cdot f_{2}$$ +where +$$\begin{align*} +\sigma_{-1}(\mu)=&b_0(\mu)\cdot S_{-1}^{(-)}\\ +\sigma_{0}(\mu)=&a_0(\mu)+b_0(\mu)\cdot S_{0}^{(-)}+b_1(\mu)\cdot S_{0}^{(+)}\\ +\sigma_{1}(\mu)=&a_1(\mu)+b_0(\mu)\cdot S_{1}^{(-)}+b_1(\mu)\cdot S_{1}^{(+)}\\ +\sigma_{2}(\mu)=&b_1(\mu)\cdot S_{2}^{(+)}\\ +\end{align*}$$ +The six parameters $S$, referred to as *cubic-slope coefficients*, are found to be +$$\begin{align*} +S_{-1}^{(-)}=&-t_0\cdot \beta_0\\ +S_0^{(-)}=&(1+t_0)\cdot \beta_0-1,\qquad +S_{0}^{(+)}=-\beta_1\\ +S_1^{(-)}=&1-\beta_0,\qquad +S_{1}^{(+)}=-1/t_1+(1+1/t_1)\cdot\beta_1\\ +S_{2}^{(+)}=&(1-\beta_1)/t_1\\ +\end{align*}$$ diff --git a/include/btwxt/grid-axis.h b/include/btwxt/grid-axis.h index 2baff76..03325ab 100644 --- a/include/btwxt/grid-axis.h +++ b/include/btwxt/grid-axis.h @@ -18,7 +18,7 @@ namespace Btwxt { enum class Method { undefined, constant, linear, cubic }; -enum class SlopeMethod { undefined, finite_diff, cardinal_0, cardinal_1, quadratic }; +enum class SlopeMethod { undefined, finite_diff, cardinal, quadratic }; class GridAxis { // A single input dimension of the grid @@ -32,7 +32,8 @@ class GridAxis { Method interpolation_method = Method::linear, Method extrapolation_method = Method::constant, std::pair extrapolation_limits = {-DBL_MAX, DBL_MAX}, - SlopeMethod slope_method = SlopeMethod::quadratic, + SlopeMethod slope_method_in = SlopeMethod::quadratic, + double slope_reduction_in = 0.0, const std::shared_ptr& logger = std::make_shared()); // Setters @@ -44,6 +45,7 @@ class GridAxis { check_extrapolation_limits(); } void set_slope_method(SlopeMethod slope_method_in); + void set_slope_reduction(double slope_reduction_in); void set_logger(std::shared_ptr logger_in) { @@ -71,6 +73,7 @@ class GridAxis { Method extrapolation_method {Method::constant}; Method interpolation_method {Method::linear}; SlopeMethod slope_method {SlopeMethod::quadratic}; + double slope_reduction = 0.0; std::pair extrapolation_limits {-DBL_MAX, DBL_MAX}; std::vector>> cubic_spacing_ratios; // Used for cubic interpolation. Outer vector is length diff --git a/src/grid-axis.cpp b/src/grid-axis.cpp index 0a1c796..c76d99d 100644 --- a/src/grid-axis.cpp +++ b/src/grid-axis.cpp @@ -11,13 +11,15 @@ GridAxis::GridAxis(std::vector values_in, Method interpolation_method, Method extrapolation_method, std::pair extrapolation_limits, - SlopeMethod slope_method, + SlopeMethod slope_method_in, + double slope_reduction_in, const std::shared_ptr& logger_in) : name(name) , values(std::move(values_in)) , extrapolation_method(extrapolation_method) , interpolation_method(interpolation_method) - , slope_method(slope_method) + , slope_method(slope_method_in) + , slope_reduction(slope_reduction_in) , extrapolation_limits(std::move(extrapolation_limits)) , cubic_spacing_ratios( std::max(static_cast(values.size()) - 1, 0), @@ -83,6 +85,11 @@ void GridAxis::set_slope_method(SlopeMethod slope_method_in) slope_method = slope_method_in; } +void GridAxis::set_slope_reduction(double slope_reduction_in) +{ + slope_reduction = slope_reduction_in; +} + void GridAxis::calculate_cubic_spacing_ratios() { if (get_length() == 1) { @@ -111,14 +118,10 @@ void GridAxis::calculate_cubic_spacing_ratios() c_0 = 0.5; break; } - case SlopeMethod::cardinal_0:{ + case SlopeMethod::cardinal:{ c_0 = 1 / (1 + t_0); break; } - case SlopeMethod::cardinal_1:{ - c_0 = 0.0; - break; - } case SlopeMethod::quadratic: default:{ c_0 = t_0 / (1 + t_0); @@ -127,8 +130,8 @@ void GridAxis::calculate_cubic_spacing_ratios() } //general - double s_m1_m = -t_0 * c_0;; - double s_1_m = 1 - c_0; + double s_m1_m = (1 - slope_reduction) * (-t_0 * c_0); + double s_1_m = (1 - slope_reduction) * (1 - c_0); ratio[0].first = s_m1_m; ratio[1].first = -(s_m1_m + s_1_m); @@ -138,7 +141,7 @@ void GridAxis::calculate_cubic_spacing_ratios() if (i_elem + 2 < values.size()) { double w_1 = values[i_elem + 2] - values[i_elem + 1]; - double t_1 = w_0 / w_1; + double t_1 = w_1 / w_0; double c_1(0.0); switch (slope_method){ @@ -147,14 +150,10 @@ void GridAxis::calculate_cubic_spacing_ratios() c_1 = 0.5; break; } - case SlopeMethod::cardinal_0:{ + case SlopeMethod::cardinal:{ c_1 = 1 / (1 + t_1); break; } - case SlopeMethod::cardinal_1:{ - c_1 = 0.0; - break; - } case SlopeMethod::quadratic: default:{ c_1 = t_1 / (1 + t_1); @@ -162,8 +161,8 @@ void GridAxis::calculate_cubic_spacing_ratios() } } - double s_0_p = -(1 - c_1); - double s_2_p = t_1 * c_1; + double s_0_p = (1 - slope_reduction) * (-c_1); + double s_2_p = (1 - slope_reduction) * ((1 - c_1) / t_1); ratio[0].second = 0.0; ratio[1].second = s_0_p; diff --git a/src/regular-grid-interpolator.cpp b/src/regular-grid-interpolator.cpp index 1f0d4b8..9d3a6dd 100644 --- a/src/regular-grid-interpolator.cpp +++ b/src/regular-grid-interpolator.cpp @@ -26,6 +26,7 @@ std::vector construct_grid_axes(const std::vector> Method::constant, std::pair {-DBL_MAX, DBL_MAX}, SlopeMethod::quadratic, + 0.0, logger_in); } return grid_axes; diff --git a/test/btwxt-tests.cpp b/test/btwxt-tests.cpp index da21883..cdc7b9c 100644 --- a/test/btwxt-tests.cpp +++ b/test/btwxt-tests.cpp @@ -36,10 +36,11 @@ TEST(BasicSplineTest, cubic_interp) Method extrapolation_method = Method::linear; Method interpolation_method = Method::cubic; SlopeMethod slope_method = SlopeMethod::quadratic; + double slope_reduction = 0.0; std::pair extrapolation_limits{-DBL_MAX, DBL_MAX}; - GridAxis axis(x_values, "x", interpolation_method, extrapolation_method, extrapolation_limits, slope_method); + GridAxis axis(x_values, "x", interpolation_method, extrapolation_method, extrapolation_limits, slope_method, slope_reduction); RegularGridInterpolator my_interpolator({axis}, values); const double epsilon = 0.0001; diff --git a/test/grid-axis-tests.cpp b/test/grid-axis-tests.cpp index 29d1d81..50036b5 100644 --- a/test/grid-axis-tests.cpp +++ b/test/grid-axis-tests.cpp @@ -47,10 +47,9 @@ TEST(GridAxis, calculate_cubic_spacing_ratios) std::vector slope_method_vector ({ SlopeMethod::quadratic, SlopeMethod::finite_diff, - SlopeMethod::cardinal_0, - SlopeMethod::cardinal_1 - + SlopeMethod::cardinal, }); + double slope_reduction = 0.0; for (auto slope_method: slope_method_vector) { GridAxis grid_axis(grid_values, @@ -59,6 +58,7 @@ TEST(GridAxis, calculate_cubic_spacing_ratios) Method::constant, {-DBL_MAX, DBL_MAX}, slope_method, + slope_reduction, std::make_shared()); std::vector> result(4, {0., 0.}); @@ -68,34 +68,27 @@ TEST(GridAxis, calculate_cubic_spacing_ratios) double w_m1 = grid_values[i_interval] - grid_values[i_interval - 1]; double t_0 = w_0 / w_m1; - double c_0(0.0); + double beta_0(0.0); switch (slope_method){ case SlopeMethod::finite_diff: { - c_0 = 0.5; - break; - } - - case SlopeMethod::cardinal_0:{ - c_0 = 1 / (1 + t_0); + beta_0 = 0.5; break; } - - case SlopeMethod::cardinal_1:{ - c_0 = 0.0; + case SlopeMethod::cardinal:{ + beta_0 = 1 / (1 + t_0); break; } - case SlopeMethod::quadratic: default:{ - c_0 = t_0 / (1 + t_0); + beta_0 = t_0 / (1 + t_0); break; } } //general - double s_m1_m = -t_0 * c_0;; - double s_1_m = 1 - c_0; + double s_m1_m = (1 - slope_reduction) * (-t_0 * beta_0); + double s_1_m = (1 - slope_reduction) * (1 - beta_0); result[0].first = s_m1_m; result[1].first = -(s_m1_m + s_1_m); result[2].first = s_1_m; @@ -103,34 +96,27 @@ TEST(GridAxis, calculate_cubic_spacing_ratios) if (i_interval + 2 < grid_values.size()){ double w_1 = grid_values[i_interval + 2] - grid_values[i_interval + 1]; - double t_1 = w_0 / w_1; + double t_1 = w_1 / w_0; - double c_1(0.0); + double beta_1(0.0); switch (slope_method){ case SlopeMethod::finite_diff:{ - c_1 = 0.5; + beta_1 = 0.5; break; } - - case SlopeMethod::cardinal_0:{ - c_1 = 1 / (1 + t_1);; + case SlopeMethod::cardinal:{ + beta_1 = 1 / (1 + t_1); break; } - - case SlopeMethod::cardinal_1:{ - c_1 = 0.0; - break; - } - case SlopeMethod::quadratic: default:{ - c_1 = t_1 / (1 + t_1); + beta_1 = t_1 / (1 + t_1); break; } } - double s_0_p = -(1 - c_1); - double s_2_p = t_1 * c_1; + double s_0_p = (1 - slope_reduction) * (-beta_1); + double s_2_p = (1 - slope_reduction) * (1 - beta_1) / t_1; result[1].second = s_0_p; result[2].second = -(s_0_p + s_2_p); diff --git a/test/implementation-tests.cpp b/test/implementation-tests.cpp index aa37b60..9b5a013 100644 --- a/test/implementation-tests.cpp +++ b/test/implementation-tests.cpp @@ -107,6 +107,7 @@ TEST_F(CubicImplementationFixture, hypercube_weigh_one_vertex) // get cubic spacing ratios std::vector>> ratios(2, {{-1.0, 0.0}, {0.0, -1.0}, {1.0, 0.0}, {0.0, 1.0}}); + double slope_reduction = 0.0; for (std::size_t i = 0; i < 2; ++i) { const std::vector& values = interpolator.get_grid_axis(i).get_values(); @@ -121,8 +122,8 @@ TEST_F(CubicImplementationFixture, hypercube_weigh_one_vertex) c_0 = t_0 / (1 + t_0); // quadratic //general - double s_m1 = -t_0 * c_0; - double s_1 = 1 - c_0; + double s_m1 = (1 - slope_reduction) * (-t_0 * c_0); + double s_1 = (1 - slope_reduction) * (1 - c_0); ratios[i][0].first = s_m1; ratios[i][1].first = -(s_m1 + s_1); @@ -130,13 +131,13 @@ TEST_F(CubicImplementationFixture, hypercube_weigh_one_vertex) } { double w_1 = values[floor_index + 2] - values[floor_index + 1]; - double t_1 = w_0 / w_1; + double t_1 = w_1 / w_0; double c_1(0.5); c_1 = t_1 / (1 + t_1); // quadratic - double s_0 = -(1 - c_1); - double s_2 = t_1 * c_1; + double s_0 = (1 - slope_reduction) * (-c_1); + double s_2 = (1 - slope_reduction) * (t_1 * (1 - c_1)); ratios[i][1].second = s_0; ratios[i][2].second = -(s_0 + s_2); @@ -203,18 +204,18 @@ TEST_F(CubicImplementationFixture, get_cubic_spacing_ratios) double w_1 = axis_values[3] - axis_values[2]; double t_0 = w_0 / w_m1; - double t_1 = w_0 / w_1; + double t_1 = w_1 / w_0; - double c_0(0.5), c_1(0.5); - c_0 = t_0 / (1 + t_0); // quadratic - c_1 = t_1 / (1 + t_1); // quadratic + double beta_0(0.0), beta_1(0.0); + beta_0 = t_0 / (1 + t_0); // quadratic + beta_1 = t_1 / (1 + t_1); // quadratic //general - double s_m1_m = -t_0 * c_0; - double s_1_m = 1 - c_0; + double s_m1_m = -t_0 * beta_0; + double s_1_m = 1 - beta_0; - double s_0_p = -(1 - c_1); - double s_2_p = t_1 * c_1; + double s_0_p = -beta_1; + double s_2_p = (1 - beta_1) / t_1; static constexpr std::size_t i_ratio = 1; @@ -358,7 +359,7 @@ TEST_F(Grid2DImplementationFixture, construct_from_axes) GridAxis ax0 = GridAxis({0, 10, 15}, "ax0"); auto logger = ax0.get_logger(); GridAxis ax1 = - GridAxis({4, 6}, "ax1", Method::linear, Method::constant, {-DBL_MAX, DBL_MAX}, SlopeMethod::quadratic, logger); + GridAxis({4, 6}, "ax1", Method::linear, Method::constant, {-DBL_MAX, DBL_MAX}, SlopeMethod::quadratic, 1.0, logger); std::vector test_axes = {ax0, ax1}; interpolator = RegularGridInterpolatorImplementation(test_axes, logger); EXPECT_EQ(interpolator.get_number_of_grid_axes(), 2u); From d16ea5ea049812ea835895dc1b01a368116163b6 Mon Sep 17 00:00:00 2001 From: Phil Ahrenkiel Date: Mon, 26 Jun 2023 10:10:24 -0600 Subject: [PATCH 6/8] Update README.md --- README.md | 66 +++++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 47 insertions(+), 19 deletions(-) diff --git a/README.md b/README.md index 71a805a..59ad06c 100644 --- a/README.md +++ b/README.md @@ -64,59 +64,87 @@ std::vector result = my_interpolator.get_values_at_target(target); std::vector result = my_interpolator({12.5, 5.1}); ``` ### Single-Axis Spline Interpolation -An interpolation along a single axis generates a function that spans the intervals between control points, typically corresponding to a linear combination of the data values. A spline is essentially a chain of such functions, such that a different set of coefficients is identified for each interval. For example, in the interval [$x_{k}$, $x_{k+1}$] we write -$$h_k(\mu_k) = \sum^{n+1}_{i=-n}\sigma_{i}^{(k)} (\mu_k)\cdot f_{k+i}$$ +An interpolation along a single axis generates a function that spans the intervals between control points, typically corresponding to a linear combination of the data values. A spline is essentially a chain of such functions, such that a different set of coefficients is identified for each interval. For example, in the interval [ $x_{k}$, $x_{k+1}$] we write +```math +h_{k}(\mu_{k}) = \sum^{n+1}_{i=-n}\sigma_{i}^{(k)} (\mu_{k})\cdot f_{k+i} +``` where $f_{k}$ is the data value at $x_{k}$ and -$$\mu_{k}={x-x_{k}\over w_{k}}$$ +```math +\mu_{k}={x-x_{k}\over w_{k}} +``` The interval width is given by $w_{k}=x_{k+1}-x_{k}$. It is useful to limit the range of the summation for practical computational. Without loss of generality, we can focus on the case $k=0$, simplifying the notation: -$$h(\mu) = \sum^{n+1}_{i=-n}\sigma_{i} (\mu)\cdot f_{i}$$ +```math +h(\mu) = \sum^{n+1}_{i=-n}\sigma_{i} (\mu)\cdot f_{i} +``` using $\mu_{0}\rightarrow\mu$ and $w_{0}\rightarrow w$. Various orders of differentiability may berequired of the spline. Continuity of the $N^{th}$-order derivative of $h(\mu)$ can be satisfied with polynomials of order $2N+1$, motivating the prominence of linear and cubic splines. #### Linear case A *linear* spline can be generated with $n=0$, such that data values $f_0$ and $f_1$ are reproduced precisely at $x_0$ and $x_1$, respectively: -$$h(\mu)=(1-\mu)\cdot f_{0}+\mu\cdot f_{1}$$ +```math +h(\mu)=(1-\mu)\cdot f_{0}+\mu\cdot f_{1} +``` The coefficients are clearly -$$\sigma_{0}(\mu)=1-\mu, \qquad \sigma_{1}(\mu)=\mu$$ +```math +\sigma_{0}(\mu)=1-\mu, \qquad \sigma_{1}(\mu)=\mu +``` The resulting interpolation preserve will consist of straight line segments spanning each interval between control points. #### Cubic case A *cubic* interpolation generates a smoother function that can be constructed to have continuous first derivative (slope), which is useful for many types of data analysis. The general form can be written $$h(\mu)=a_0(\mu)\cdot f_0+a_1(\mu)\cdot f_1+b_0(\mu)\cdot w_0\cdot f'_0++b_1(\mu)\cdot w_0\cdot f'_1$$ with -$$\begin{align*}a_0(\mu)&=(1-\mu)\cdot [1+\mu\cdot (1-2\mu)]\\ +```math +\begin{align*}a_0(\mu)&=(1-\mu)\cdot [1+\mu\cdot (1-2\mu)]\\ a_1(\mu)&=\mu\cdot[1-(1-\mu)\cdot(1-2\mu)]\\ b_0(\mu)&=(1-\mu)^2\cdot \mu\\ b_1(\mu)&=-(1-\mu)\cdot \mu^2\\ -\end{align*}$$ +\end{align*} +``` Nonetheless, the slope is not uniquely defined for a discrete data set, so a reasonable estimation must be generated by inspection of neighboring control points. In fact, we need consider only the additional data values at $x_{k-1}$ and $x_{k+1}$ to develop an expression for $f'_{k}$: -$$f'_{k}=(1-\alpha)\cdot [(1-\beta_{k})\cdot r_{k}+\beta_{k}\cdot r_{k-1}]$$ +```math +f'_{k}=(1-\alpha)\cdot [(1-\beta_{k})\cdot r_{k}+\beta_{k}\cdot r_{k-1}] +``` where -$$r_{k}=\frac{f_{k+1}-f_{k}}{w_{k}}$$ +```math +r_{k}=\frac{f_{k+1}-f_{k}}{w_{k}} +``` The factor $1-\alpha$ uniformly scales the slope at control points. A default value of $\alpha =0$ is assigned. Various criteria for the specification of $\beta _{k}$ may be applicable in different contexts; **btwxt** allows provides options for slope evaluation. The coefficient may depend on the relative widths of neighboring intervals. We define -$$t_{k}=\frac{w_{k}}{w_{k-1}}$$ +```math +t_{k}=\frac{w_{k}}{w_{k-1}} +``` The default setting computes the slope by a local fit to a *quadratic* function about each control point, which corresponds to -$$\beta_{k}=\frac{t_{k}}{1 + t_{k}}$$ +```math +\beta_{k}=\frac{t_{k}}{1 + t_{k}} +``` The *cardinal* method is realized with -$$\beta_{k}=\frac{1}{1 + t_{k}}$$ +```math +\beta_{k}=\frac{1}{1 + t_{k}} +``` The slope in *finite-difference* methods is computed with $\beta _{k}=0.5$\, imposing no width dependence. -We can now represent the cubic spline over [$x_{0}$, $x_{1}$] as the sum -$$h(\mu)=\sigma_{-1}(\mu)\cdot f_{-1}+\sigma_{0}(\mu)\cdot f_{0}+\sigma_{1}(\mu)\cdot f_{1}+\sigma_{2}(\mu)\cdot f_{2}$$ +We can now represent the cubic spline over [ $x_{0}$, $x_{1}$] as the sum +```math +h(\mu)=\sigma_{-1}(\mu)\cdot f_{-1}+\sigma_{0}(\mu)\cdot f_{0}+\sigma_{1}(\mu)\cdot f_{1}+\sigma_{2}(\mu)\cdot f_{2} +``` where -$$\begin{align*} +```math +\begin{align*} \sigma_{-1}(\mu)=&b_0(\mu)\cdot S_{-1}^{(-)}\\ \sigma_{0}(\mu)=&a_0(\mu)+b_0(\mu)\cdot S_{0}^{(-)}+b_1(\mu)\cdot S_{0}^{(+)}\\ \sigma_{1}(\mu)=&a_1(\mu)+b_0(\mu)\cdot S_{1}^{(-)}+b_1(\mu)\cdot S_{1}^{(+)}\\ \sigma_{2}(\mu)=&b_1(\mu)\cdot S_{2}^{(+)}\\ -\end{align*}$$ +\end{align*} +``` The six parameters $S$, referred to as *cubic-slope coefficients*, are found to be -$$\begin{align*} +```math +\begin{align*} S_{-1}^{(-)}=&-t_0\cdot \beta_0\\ S_0^{(-)}=&(1+t_0)\cdot \beta_0-1,\qquad S_{0}^{(+)}=-\beta_1\\ S_1^{(-)}=&1-\beta_0,\qquad S_{1}^{(+)}=-1/t_1+(1+1/t_1)\cdot\beta_1\\ S_{2}^{(+)}=&(1-\beta_1)/t_1\\ -\end{align*}$$ +\end{align*} +``` From b0feceaf3aca1a1a9110fd594b6632e6dea01e61 Mon Sep 17 00:00:00 2001 From: Phil Ahrenkiel Date: Mon, 26 Jun 2023 12:13:34 -0600 Subject: [PATCH 7/8] Added scope for std::isnan --- test/implementation-tests.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/implementation-tests.cpp b/test/implementation-tests.cpp index 9b5a013..c47a622 100644 --- a/test/implementation-tests.cpp +++ b/test/implementation-tests.cpp @@ -244,7 +244,7 @@ TEST_F(CubicImplementationFixture, null_checking_calculations) target = {7, 3}; interpolator.set_target(target); auto result = interpolator.get_results(); - EXPECT_TRUE(isnan(result[2])); + EXPECT_TRUE(std::isnan(result[2])); } TEST_F(EmptyGridImplementationFixture, locate_coordinates) From 5045a2e8aa45ade1d4d6c02f13afc5c96541d836 Mon Sep 17 00:00:00 2001 From: Phil Ahrenkiel Date: Tue, 27 Jun 2023 12:49:49 -0600 Subject: [PATCH 8/8] Update README.md --- README.md | 135 ++++++++++++------------------------------------------ 1 file changed, 30 insertions(+), 105 deletions(-) diff --git a/README.md b/README.md index 59ad06c..5ab196a 100644 --- a/README.md +++ b/README.md @@ -1,28 +1,35 @@ [![Build and Test](https://github.com/bigladder/btwxt/actions/workflows/build-and-test.yml/badge.svg)](https://github.com/bigladder/btwxt/actions/workflows/build-and-test.yml) [![codecov](https://codecov.io/gh/bigladder/btwxt/branch/master/graph/badge.svg)](https://codecov.io/gh/bigladder/btwxt) -# **btwxt** +# Btwxt ## General-purpose, N-dimensional interpolation library... -**btwxt** is a free and open source C++ library to perform multivariate interpolation of one or more functions represented by known values on provided data sets associated with coordinates on a regular (rectilinear) grid. The primary class is a RegularGridInterpolator constructed from: +Btwxt is a free and open source C++ library to perform numerical interpolation on a regular grid data set. The primary +class is a RegularGridInterpolator constructed from: -1. a set of one or more *grid axes*, representing the unique variates for interpolation, each specfied through a set of *grid-axis coordinates*, and -2. a collection of one or more *data sets*, repesenting the function values at each of the *control points* specified by ordered pairs of the grid-axis coordinates. +1. a set of N grid axes representing the values of the independent variables, and +2. a collection of one or more grid point data sets aligned to the grid axes. -A RegularGridInterpolator object can then be queried repeatedly for function values within the space of the grid axes. Interpolation is performed for when a grid coordinate is inside the range of a particular grid axis; extrapolation is performed (up to a user-defined limit) when the grid coordinate is outside that range. A general query of a multivariate function may entail any combination of interpolation and extrapolation on the collection of grid axes. +A RegularGridInterpolator object can then be queried repeatedly for interpolated values that +fall inside its grid points, or extrapolated values beyond the grid points (up to a limit defined by the user). -**btwxt** supports various methods for both interpolation and extrpolation; different methods can be specifed for each axis. Interpolation can be performed using either linear or cubic spline. Whereas both linear and cubic splines maintain continuity of a function along the associated axis, cubic splines also maintain continuty of the function's first derivative on that axis. Extrapolation can be either constant or linear, with the option to specify finite extrapolation limits for a given axis. +Btwxt supports linear and cubic spline (Catmull-Rom) interpolation methods. Different methods can be specified for each +axis. The API also allows specification of preferred extrapolation methods (constant or linear)--again independently for +each axis--and extrapolation limits. -The grid is required to be rectilinear, but the axes are not required to have regular (uniform) spacing. However, the grid coordinates provided for each axis must be strictly increasing. +The grid is required to be regular (rectilinear), but the axes are not required to have uniform spacing. Each axis need +only be strictly increasing. -Data sets are imported as a single long vector (i.e., a flattened array). The user must ensure that values in a data set are ordered, with respect to the orders of both grid coordinates and grid axes, such that each value is correctly associated with the appropriate control point. **btwxt** verifies that the number of values with each data set corresponds with the total number of control points. +The grid point data sets are imported as a single long vector (i.e., a flattened array). It is the user's responsibility +to ensure that the grid point data is ordered according to the order of grid axes. Btwxt ensures that the number of +grid points aligns with the total combined points of the grid axes. ## How to Use ### API -The snippet below illustrates instantiation of an individual axis. +An individual axis can be instantiated with: ```c++ std::vector one_axis{6, 10, 15}; @@ -32,119 +39,37 @@ GridAxis(one_axis, "x"); // Construct a new axis named "x" // extrapolation limits default to {-DBL_MAX, DBL_MAX} Method extrapolation_method = Method::linear; Method interpolation_method = Method::cubic; -SlopeMethod slope_method = SlopeMethod::quadratic; -double slope_reduction = 0.0; std::pair extrapolation_limits{0, 20}; -GridAxis(one_axis, "x", interpolation_method, extrapolation_method, extrapolation_limits, -slope_method, slope_reduction); +GridAxis(one_axis, "x", interpolation_method, extrapolation_method, extrapolation_limits); ``` -A RegularGridInterpolator holds a collection of axes, which specify the control points, and a collection of data sets containing function values for each of the control points. The snippet below illustrates instantiation of a two-dimensional grid, an associated data set, and their assocation with a RegularGridInterpolator. +A RegularGridInterpolator holds a collection of axes and the corresponding grid point data at all permutations of the grid axes' values. ```c++ std::vector> vector_grid = {{6, 10, 15}, {4, 6, 9}}; -std::vector> values = {2, 4, 7, 1, 3, 6}; +std::vector> values = {...}; RegularGridInterpolator my_interpolator(vector_grid, values); ``` -Once a RegularGridInterpolator has been synthesized, the interpolated can be queried for function values at target grid-axis coordinates. Typical syntax to perform such a query on a 2-D grid is shown in the snippet below. +Once you have an RegularGridInterpolator object, you can set a point to interpolate to: ```c++ -// 1. set the target before requesting values: std::vector target{12.5, 5.1}; my_interpolator.set_target(target); std::vector result = my_interpolator.get_values_at_target(); -// 2. set the target when requesting values: -std::vector target{12.5, 5.1}; -std::vector result = my_interpolator.get_values_at_target(target); - -// 3. set the target by applying the () operator: -std::vector result = my_interpolator({12.5, 5.1}); -``` -### Single-Axis Spline Interpolation -An interpolation along a single axis generates a function that spans the intervals between control points, typically corresponding to a linear combination of the data values. A spline is essentially a chain of such functions, such that a different set of coefficients is identified for each interval. For example, in the interval [ $x_{k}$, $x_{k+1}$] we write -```math -h_{k}(\mu_{k}) = \sum^{n+1}_{i=-n}\sigma_{i}^{(k)} (\mu_{k})\cdot f_{k+i} -``` -where $f_{k}$ is the data value at $x_{k}$ and -```math -\mu_{k}={x-x_{k}\over w_{k}} -``` -The interval width is given by $w_{k}=x_{k+1}-x_{k}$. It is useful to limit the range of the summation for practical computational. Without loss of generality, we can focus on the case $k=0$, simplifying the notation: -```math -h(\mu) = \sum^{n+1}_{i=-n}\sigma_{i} (\mu)\cdot f_{i} +// or set the target when requesting values: +std::vector new_target{11.7, 6.1}; +result = my_interpolator.get_values_at_target(new_target); ``` -using $\mu_{0}\rightarrow\mu$ and $w_{0}\rightarrow w$. -Various orders of differentiability may berequired of the spline. Continuity of the $N^{th}$-order derivative of $h(\mu)$ can be satisfied with polynomials of order $2N+1$, motivating the prominence of linear and cubic splines. -#### Linear case -A *linear* spline can be generated with $n=0$, such that data values $f_0$ and $f_1$ are reproduced precisely at $x_0$ and $x_1$, respectively: -```math -h(\mu)=(1-\mu)\cdot f_{0}+\mu\cdot f_{1} -``` -The coefficients are clearly -```math -\sigma_{0}(\mu)=1-\mu, \qquad \sigma_{1}(\mu)=\mu -``` -The resulting interpolation preserve will consist of straight line segments spanning each interval between control points. - -#### Cubic case -A *cubic* interpolation generates a smoother function that can be constructed to have continuous first derivative (slope), which is useful for many types of data analysis. The general form can be written -$$h(\mu)=a_0(\mu)\cdot f_0+a_1(\mu)\cdot f_1+b_0(\mu)\cdot w_0\cdot f'_0++b_1(\mu)\cdot w_0\cdot f'_1$$ -with -```math -\begin{align*}a_0(\mu)&=(1-\mu)\cdot [1+\mu\cdot (1-2\mu)]\\ -a_1(\mu)&=\mu\cdot[1-(1-\mu)\cdot(1-2\mu)]\\ -b_0(\mu)&=(1-\mu)^2\cdot \mu\\ -b_1(\mu)&=-(1-\mu)\cdot \mu^2\\ -\end{align*} -``` -Nonetheless, the slope is not uniquely defined for a discrete data set, so a reasonable estimation must be generated by inspection of neighboring control points. -In fact, we need consider only the additional data values at $x_{k-1}$ and $x_{k+1}$ to develop an expression for $f'_{k}$: -```math -f'_{k}=(1-\alpha)\cdot [(1-\beta_{k})\cdot r_{k}+\beta_{k}\cdot r_{k-1}] -``` -where -```math -r_{k}=\frac{f_{k+1}-f_{k}}{w_{k}} -``` -The factor $1-\alpha$ uniformly scales the slope at control points. A default value of $\alpha =0$ is assigned. Various criteria for the specification of $\beta _{k}$ may be applicable in different contexts; **btwxt** allows provides options for slope evaluation. The coefficient may depend on the relative widths of neighboring intervals. We define -```math -t_{k}=\frac{w_{k}}{w_{k-1}} -``` -The default setting computes the slope by a local fit to a *quadratic* function about each control point, which corresponds to -```math -\beta_{k}=\frac{t_{k}}{1 + t_{k}} -``` -The *cardinal* method is realized with -```math -\beta_{k}=\frac{1}{1 + t_{k}} -``` -The slope in *finite-difference* methods is computed with $\beta _{k}=0.5$\, imposing no width dependence. +We have overloaded the () operator on the RegularGridInterpolator to perform the interpolation, enabling the following: -We can now represent the cubic spline over [ $x_{0}$, $x_{1}$] as the sum -```math -h(\mu)=\sigma_{-1}(\mu)\cdot f_{-1}+\sigma_{0}(\mu)\cdot f_{0}+\sigma_{1}(\mu)\cdot f_{1}+\sigma_{2}(\mu)\cdot f_{2} -``` -where -```math -\begin{align*} -\sigma_{-1}(\mu)=&b_0(\mu)\cdot S_{-1}^{(-)}\\ -\sigma_{0}(\mu)=&a_0(\mu)+b_0(\mu)\cdot S_{0}^{(-)}+b_1(\mu)\cdot S_{0}^{(+)}\\ -\sigma_{1}(\mu)=&a_1(\mu)+b_0(\mu)\cdot S_{1}^{(-)}+b_1(\mu)\cdot S_{1}^{(+)}\\ -\sigma_{2}(\mu)=&b_1(\mu)\cdot S_{2}^{(+)}\\ -\end{align*} -``` -The six parameters $S$, referred to as *cubic-slope coefficients*, are found to be -```math -\begin{align*} -S_{-1}^{(-)}=&-t_0\cdot \beta_0\\ -S_0^{(-)}=&(1+t_0)\cdot \beta_0-1,\qquad -S_{0}^{(+)}=-\beta_1\\ -S_1^{(-)}=&1-\beta_0,\qquad -S_{1}^{(+)}=-1/t_1+(1+1/t_1)\cdot\beta_1\\ -S_{2}^{(+)}=&(1-\beta_1)/t_1\\ -\end{align*} +```c++ +std::vector target{12.5, 5.1}; +my_interpolator.set_target(target); +std::vector result = my_interpolator(); +std::vector new_target{11.7, 6.1}; +result = my_interpolator(new_target); ```