Fix slope calculation; add alternative methods

.gitignore

@@ -2,4 +2,4 @@ builds/
 build/
 .idea/
 .vs/
-.vscode/
+.vscode/

README.md

@@ -72,4 +72,4 @@ my_interpolator.set_target(target);
 std::vector<double> result = my_interpolator();
 std::vector<double> new_target{11.7, 6.1};
 result = my_interpolator(new_target);
-```
+```

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, quadratic };
 
 class GridAxis {
     // A single input dimension of the grid
@@ -31,6 +32,8 @@ class GridAxis {
         Method interpolation_method = Method::linear,
         Method extrapolation_method = Method::constant,
         std::pair<double, double> extrapolation_limits = {-DBL_MAX, DBL_MAX},
+        SlopeMethod slope_method_in = SlopeMethod::quadratic,
+        double slope_reduction_in = 0.0,
         const std::shared_ptr<Courierr::Courierr>& logger = std::make_shared<BtwxtLogger>());
 
     // Setters
@@ -41,6 +44,8 @@ class GridAxis {
         extrapolation_limits = limits;
         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<Courierr::Courierr> logger_in)
     {
@@ -58,20 +63,22 @@ class GridAxis {
         return extrapolation_limits;
     }
 
-    [[nodiscard]] const std::vector<double>&
-    get_cubic_spacing_ratios(std::size_t floor_or_ceiling) const;
+    [[nodiscard]] const std::vector<std::pair<double,double>>&
+        get_cubic_spacing_ratios(std::size_t elem_index) const;
+
     std::string name;
 
   private:
     std::vector<double> values;
     Method extrapolation_method {Method::constant};
     Method interpolation_method {Method::linear};
+    SlopeMethod slope_method {SlopeMethod::quadratic};
+    double slope_reduction = 0.0;
     std::pair<double, double> extrapolation_limits {-DBL_MAX, DBL_MAX};
-    std::vector<std::vector<double>>
-        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<std::vector<std::pair<double,double>>>
+        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<Courierr::Courierr> logger;
     void calculate_cubic_spacing_ratios();
     void check_grid_sorted();
@@ -107,4 +114,4 @@ std::vector<double> linspace(double start, double stop, std::size_t number_of_po
 
 } // namespace Btwxt
 
-#endif // define GRID_AXIS_H_
+#endif // define GRID_AXIS_H_

src/grid-axis.cpp

@@ -11,23 +11,29 @@ GridAxis::GridAxis(std::vector<double> values_in,
                    Method interpolation_method,
                    Method extrapolation_method,
                    std::pair<double, double> extrapolation_limits,
+                   SlopeMethod slope_method_in,
+                   double slope_reduction_in,
                    const std::shared_ptr<Courierr::Courierr>& logger_in)
     : name(name)
     , values(std::move(values_in))
     , extrapolation_method(extrapolation_method)
     , interpolation_method(interpolation_method)
+    , slope_method(slope_method_in)
+    , slope_reduction(slope_reduction_in)
     , extrapolation_limits(std::move(extrapolation_limits))
     , cubic_spacing_ratios(
-          2, std::vector<double>(std::max(static_cast<int>(values.size()) - 1, 0), 1.0))
-    , logger(logger_in)
-{
+        std::max(static_cast<int>(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();
     }
@@ -74,6 +80,16 @@ 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::set_slope_reduction(double slope_reduction_in)
+{
+    slope_reduction = slope_reduction_in;
+}
+
 void GridAxis::calculate_cubic_spacing_ratios()
 {
     if (get_length() == 1) {
@@ -86,23 +102,80 @@ 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]);
+    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 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:{
+                    c_0 = 1 / (1 + t_0);
+                    break;
+                }
+                case SlopeMethod::quadratic:
+                default:{
+                    c_0 = t_0 / (1 + t_0);
+                    break;
+                }
+            }
+
+            //general
+            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);
+            ratio[2].first = s_1_m;
+            ratio[3].first = 0.0;
         }
-        if (i + 2 != values.size()) {
-            cubic_spacing_ratios[ceiling][i] = center_spacing / (values[i + 2] - values[i]);
+        if (i_elem + 2 < values.size())
+        {
+            double w_1 = values[i_elem + 2] - values[i_elem + 1];
+            double t_1 = w_1 / w_0;
+
+            double c_1(0.0);
+            switch (slope_method){
+                case SlopeMethod::finite_diff:
+                {
+                    c_1 = 0.5;
+                    break;
+                }
+                case SlopeMethod::cardinal:{
+                    c_1 = 1 / (1 + t_1);
+                    break;
+                }
+                case SlopeMethod::quadratic:
+                default:{
+                    c_1 = t_1 / (1 + t_1);
+                    break;
+                }
+            }
+
+            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;
+            ratio[2].second = -(s_0_p + s_2_p);
+            ratio[3].second = s_2_p;
         }
     }
 }
 
-const std::vector<double>&
-GridAxis::get_cubic_spacing_ratios(const std::size_t floor_or_ceiling) const
+const std::vector<std::pair<double,double>>&
+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()

src/regular-grid-interpolator-implementation.h

@@ -137,11 +137,11 @@ class RegularGridInterpolatorImplementation {
         return hypercube;
     };
 
-    [[nodiscard]] const std::vector<double>&
-    get_axis_cubic_spacing_ratios(std::size_t axis_index, std::size_t floor_or_ceiling) const
+    [[nodiscard]] const std::vector<std::pair<double,double>>&
+    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<double>& get_grid_point_data(const std::vector<std::size_t>& coordinates);

src/regular-grid-interpolator.cpp

@@ -25,6 +25,8 @@ std::vector<GridAxis> construct_grid_axes(const std::vector<std::vector<double>>
                                Method::linear,
                                Method::constant,
                                std::pair<double, double> {-DBL_MAX, DBL_MAX},
+                               SlopeMethod::quadratic,
+                               0.0,
                                logger_in);
     }
     return grid_axes;
@@ -97,7 +99,7 @@ RegularGridInterpolator::RegularGridInterpolator(const std::vector<GridAxis>& gr
 {
 }
 
-RegularGridInterpolator::~RegularGridInterpolator() = default; 
+RegularGridInterpolator::~RegularGridInterpolator() = default;
 
 RegularGridInterpolatorImplementation::RegularGridInterpolatorImplementation(
     const std::vector<GridAxis>& grid_axes,
@@ -119,7 +121,7 @@ RegularGridInterpolatorImplementation::RegularGridInterpolatorImplementation(
     , weighting_factors(number_of_grid_axes, std::vector<double>(4, 0.))
     , results(number_of_grid_point_data_sets, 0.)
     , interpolation_coefficients(number_of_grid_axes, std::vector<double>(2, 0.))
-    , cubic_slope_coefficients(number_of_grid_axes, std::vector<double>(2, 0.))
+    , cubic_slope_coefficients(number_of_grid_axes, std::vector<double>(4, 0.))
     , logger(logger)
 {
     set_axis_sizes();
@@ -151,24 +153,25 @@ 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;
     }
 }
 
 // Public manipulation methods
-
 std::size_t
 RegularGridInterpolator::add_grid_point_data_set(const std::vector<double>& grid_point_data_vector,
                                                  const std::string& name)
@@ -678,39 +681,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<double> linear_interpolation_factors = {0.0, 1.0 - mu, mu, 0.0};
+
+            double coef = 1- 2 * mu; //(1 - mu) - mu;
+            std::vector<double> 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<double> 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<double> 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)
     }
 }