Fix slope calculation; add alternative methods

.gitignore

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

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 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<double> one_axis{6, 10, 15};
@@ -39,37 +32,119 @@ 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<double, double> 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 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<std::vector<double>> vector_grid = {{6, 10, 15},
                                                  {4, 6, 9}};
-std::vector<std::vector<double>> values = {...};
+std::vector<std::vector<double>> 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<double> target{12.5, 5.1};
 my_interpolator.set_target(target);
 std::vector<double> result = my_interpolator.get_values_at_target();
 
-// or set the target when requesting values:
-std::vector<double> new_target{11.7, 6.1};
-result = my_interpolator.get_values_at_target(new_target);
+// 2. set the target when requesting values:
+std::vector<double> target{12.5, 5.1};
+std::vector<double> result = my_interpolator.get_values_at_target(target);
+
+// 3. set the target by applying the () operator:
+std::vector<double> 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}
 ```
+using $\mu_{0}\rightarrow\mu$ and $w_{0}\rightarrow w$. 
 
-We have overloaded the () operator on the RegularGridInterpolator to perform the interpolation, enabling the following:
+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.
 
-```c++
-std::vector<double> target{12.5, 5.1};
-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);
-```
+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*}
+```

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()