interpolant positivity and smooth extrapolation (#37)

* switch to linear interpolation if result is negative

* use positive constraints for spline coefficients

* Gaussian tail for extrapolation

inst/include/interpolation.hpp

@@ -11,9 +11,7 @@
 class InterpolationGrid1d
 {
 public:
-    InterpolationGrid1d()
-    {
-    }
+    InterpolationGrid1d() {}
 
     InterpolationGrid1d(const Eigen::VectorXd& grid_points,
                         const Eigen::VectorXd& values,
@@ -187,44 +185,44 @@ inline double InterpolationGrid1d::cubic_integral(const double& lower,
 //! @param vals length 4 vector of function values.
 //! @param grid length 4 vector of grid points.
 inline Eigen::VectorXd
-    InterpolationGrid1d::find_coefs(const Eigen::VectorXd& vals,
-                                    const Eigen::VectorXd& grid)
-    {
-        Eigen::VectorXd a(4);
-
-        double dt0 = grid(1) - grid(0);
-        double dt1 = grid(2) - grid(1);
-        double dt2 = grid(3) - grid(2);
-
-        /* check for repeated points (important for boundaries) */
-        if (dt1 <= 0)
-            dt1 = 1.0;
-        if (dt0 <= 0)
-            dt0 = dt1;
-        if (dt2 <= 0)
-            dt2 = dt1;
-
-        // compute tangents when parameterized in (t1,t2)
-        double dx1 = (vals(1) - vals(0)) / dt0;
+InterpolationGrid1d::find_coefs(const Eigen::VectorXd& vals,
+                                const Eigen::VectorXd& grid)
+{
+    double dt0 = grid(1) - grid(0);
+    double dt1 = grid(2) - grid(1);
+    double dt2 = grid(3) - grid(2);
+
+    // compute tangents when parameterized in (t1,t2)
+    // for smooth extrapolation, derivative is set to zero at boundary
+    double dx1 = 0.0, dx2 = 0.0;
+    if (dt0 > 0) {
+        dx1 = (vals(1) - vals(0)) / dt0;
         dx1 -= (vals(2) - vals(0)) / (dt0 + dt1);
         dx1 += (vals(2) - vals(1)) / dt1;
-        double dx2 = (vals(2) - vals(1)) / dt1;
+    }
+    if (dt2 > 0) {
+        dx2 = (vals(2) - vals(1)) / dt1;
         dx2 -= (vals(3) - vals(1)) / (dt1 + dt2);
         dx2 += (vals(3) - vals(2)) / dt2;
+    }
 
-        // rescale tangents for parametrization in (0,1)
-        dx1 *= dt1;
-        dx2 *= dt1;
+    // rescale tangents for parametrization in (0,1)
+    dx1 *= dt1;
+    dx2 *= dt1;
 
-        // compute coefficents
-        a(0) = vals(1);
-        a(1) = dx1;
-        a(2) = -3 * vals(1) + 3 * vals(2) - 2 * dx1 - dx2;
-        a(3) = 2 * vals(1) - 2 * vals(2) + dx1 + dx2;
+    // ensure positivity (Schmidt and Hess, DOI:10.1007/bf01934097)
+    dx1 = std::max(dx1, -3 * vals(1));
+    dx2 = std::min(dx2, 3 * vals(2));
 
-        return a;
-    }
+    // compute coefficents
+    Eigen::VectorXd a(4);
+    a(0) = vals(1);
+    a(1) = dx1;
+    a(2) = -3 * (vals(1) - vals(2)) - 2 * dx1 - dx2;
+    a(3) = 2 * (vals(1) - vals(2)) + dx1 + dx2;
 
+    return a;
+}
 //! Interpolate on 4 points
 //!
 //! @param x evaluation point.
@@ -234,8 +232,15 @@ inline double InterpolationGrid1d::interp_on_grid(const double& x,
                                                   const Eigen::VectorXd& vals,
                                                   const Eigen::VectorXd& grid)
 {
+    double xev = (x - grid(1)) / (grid(2) - grid(1));
+    // use Gaussian tail for extrapolation
+    if (xev <= 0) {
+        return vals(1) * std::exp(-0.5 * xev * xev);
+    } else if (xev >= 1) {
+        return vals(2) * std::exp(-0.5 * xev * xev);
+    }
+
     Eigen::VectorXd a = find_coefs(vals, grid);
-    double xev = std::fmax((x - grid(1)), 0) / (grid(2) - grid(1));
     return cubic_poly(xev, a);
 }
 

inst/include/lpdens.hpp

@@ -319,46 +319,33 @@ inline Eigen::VectorXd LPDens1d::construct_grid_points(
         const Eigen::VectorXd& x,
         const Eigen::VectorXd& weights)
 {
-    // set up grid
-    size_t grid_size = 52;
-    Eigen::VectorXd grid_points(grid_size);
-
-    // determine "inner" grid by sample quantiles
-    // (need to leave room for boundary extensions)
-    if (std::isnan(xmin_))
-        grid_size -= 3;
-    if (std::isnan(xmax_))
-        grid_size -= 3;
-
     // hybrid grid: quantiles and two equally spaced points between them
     auto qgrid = stats::quantile(
-        x, Eigen::VectorXd::LinSpaced((grid_size - 1)/3 + 1, 0, 1), weights);
+        x, Eigen::VectorXd::LinSpaced(14, 0, 1), weights);
+    size_t grid_size = 53;
     Eigen::VectorXd inner_grid(grid_size);
     for (unsigned i = 0; i < qgrid.size() - 1; i++) {
-        inner_grid.segment(i * 3, 4) =
-            Eigen::VectorXd::LinSpaced(4, qgrid(i), qgrid(i + 1));
+        inner_grid.segment(i * 4, 5) =
+            Eigen::VectorXd::LinSpaced(5, qgrid(i), qgrid(i + 1));
     }
 
     // extend grid where there's no boundary
     double range = inner_grid.maxCoeff() - inner_grid.minCoeff();
     Eigen::VectorXd lowr_ext, uppr_ext;
     if (std::isnan(xmin_)) {
         // no left boundary -> add a few points to the left
-        lowr_ext = Eigen::VectorXd(3);
-        double step = inner_grid[1] - inner_grid[0];
-        lowr_ext[2] = inner_grid[0] - step;
-        lowr_ext[1] = inner_grid[0] - 2 * step;
-        lowr_ext[0] = lowr_ext[1] - std::max(0.25 * range, step);
+        lowr_ext = Eigen::VectorXd(2);
+        lowr_ext[1] = inner_grid[0] - 1 * bw_;
+        lowr_ext[0] = inner_grid[0] - 2 * bw_;
     }
     if (std::isnan(xmax_)) {
         // no right boundary -> add a few points to the right
-        uppr_ext = Eigen::VectorXd(3);
-        double step = inner_grid[grid_size - 1] - inner_grid[grid_size - 2];
-        uppr_ext[0] = inner_grid[grid_size - 1] + step;
-        uppr_ext[1] = inner_grid[grid_size - 1] + 2 * step;
-        uppr_ext[2] = uppr_ext[1] + std::max(0.25 * range, step);
+        uppr_ext = Eigen::VectorXd(2);
+        uppr_ext[0] = inner_grid[grid_size - 1] + 1 * bw_;
+        uppr_ext[1] = inner_grid[grid_size - 1] + 2 * bw_;
     }
 
+    Eigen::VectorXd grid_points(grid_size + uppr_ext.size() + lowr_ext.size());
     grid_points << lowr_ext, inner_grid, uppr_ext;
     return grid_points;
 }