Periodic 2D and 3D spline filters

openmmapi/include/openmm/internal/SplineFitter.h

@@ -43,6 +43,18 @@ namespace OpenMM {
 
 class OPENMM_EXPORT SplineFitter {
 public:
+    /**
+     * Fit a cubic spline to a set of data points.  The resulting spline interpolates all the
+     * data points and has a continuous second derivative everywhere. The second derivatives are
+     * identical at the end points if periodic=true or 0 at the end points if periodic=false.
+     *
+     * @param x        the values of the independent variable at the data points to interpolate.  They must
+     *                 be strictly increasing: x[i] > x[i-1].
+     * @param y        the values of the dependent variable at the data points to interpolate
+     * @param periodic whether the interpolated function is periodic
+     * @param deriv    on exit, this contains the second derivative of the spline at each of the data points
+     */
+    static void createSpline(const std::vector<double>& x, const std::vector<double>& y, bool periodic, std::vector<double>& deriv);
     /**
      * Fit a natural cubic spline to a set of data points.  The resulting spline interpolates all the
      * data points, has a continuous second derivative everywhere, and has a second derivative of 0 at
@@ -86,6 +98,21 @@ public:
      * @return the value of the spline's derivative  at the specified point
      */
     static double evaluateSplineDerivative(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& deriv, double t);
+    /**
+     * Fit a cubic spline surface f(x,y) to a 2D set of data points.  The resulting spline interpolates all the
+     * data points and has a continuous second derivative everywhere. The second derivatives are identical at
+     * the boundary if periodic=true or 0 at the boundary if periodic=false.
+     *
+     * @param x        the values of the first independent variable at the data points to interpolate.  They must
+     *                 be strictly increasing: x[i] > x[i-1].
+     * @param y        the values of the second independent variable at the data points to interpolate.  They must
+     *                 be strictly increasing: y[i] > y[i-1].
+     * @param values   the values of the dependent variable at the data points to interpolate.  They must be ordered
+     *                 so that values[i+xsize*j] = f(x[i],y[j]), where xsize is the length of x.
+     * @param periodic whether the interpolated function is periodic
+     * @param c        on exit, this contains the spline coefficients at each of the data points
+     */
+    static void create2DSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& values, bool periodic, std::vector<std::vector<double> >& c);
     /**
      * Fit a natural cubic spline surface f(x,y) to a 2D set of data points.  The resulting spline interpolates all the
      * data points, has a continuous second derivative everywhere, and has a second derivative of 0 at the boundary.
@@ -124,6 +151,24 @@ public:
      * @param dy     on exit, the y derivative of the spline at the specified point
      */
     static void evaluate2DSplineDerivatives(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& values, const std::vector<std::vector<double> >& c, double u, double v, double& dx, double& dy);
+    /**
+     * Fit a cubic spline surface f(x,y,z) to a 3D set of data points.  The resulting spline interpolates all the
+     * data points and has a continuous second derivative everywhere. The second derivatives are identical at
+     * the boundary if periodic=true or 0 at the boundary if periodic=false.
+     *
+     * @param x        the values of the first independent variable at the data points to interpolate.  They must
+     *                 be strictly increasing: x[i] > x[i-1].
+     * @param y        the values of the second independent variable at the data points to interpolate.  They must
+     *                 be strictly increasing: y[i] > y[i-1].
+     * @param z        the values of the third independent variable at the data points to interpolate.  They must
+     *                 be strictly increasing: z[i] > z[i-1].
+     * @param values   the values of the dependent variable at the data points to interpolate.  They must be ordered
+     *                 so that values[i+xsize*j+xsize*ysize*k] = f(x[i],y[j],z[k]), where xsize is the length of x
+     *                 and ysize is the length of y.
+     * @param periodic whether the interpolated function is periodic
+     * @param c        on exit, this contains the spline coefficients at each of the data points
+     */
+    static void create3DSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, const std::vector<double>& values, bool periodic, std::vector<std::vector<double> >& c);
     /**
      * Fit a natural cubic spline surface f(x,y,z) to a 3D set of data points.  The resulting spline interpolates all the
      * data points, has a continuous second derivative everywhere, and has a second derivative of 0 at the boundary.

openmmapi/src/SplineFitter.cpp

@@ -37,6 +37,13 @@
 using namespace OpenMM;
 using namespace std;
 
+void SplineFitter::createSpline(const vector<double>& x, const vector<double>& y, bool periodic, vector<double>& deriv) {
+    if (periodic)
+        SplineFitter::createPeriodicSpline(x, y, deriv);
+    else
+        SplineFitter::createNaturalSpline(x, y, deriv);
+}
+
 void SplineFitter::createNaturalSpline(const vector<double>& x, const vector<double>& y, vector<double>& deriv) {
     int n = x.size();
     if (y.size() != n)
@@ -183,16 +190,16 @@ void SplineFitter::solveTridiagonalMatrix(const vector<double>& a, const vector<
         sol[i] = (rhs[i]-a[i]*sol[i-1])/beta;
     }
 
-    // Perform backsubstitation.
+    // Perform backsubstitution.
 
     for (int i = n-2; i >= 0; i--)
         sol[i] -= gamma[i+1]*sol[i+1];
 }
 
-void SplineFitter::create2DNaturalSpline(const vector<double>& x, const vector<double>& y, const vector<double>& values, vector<vector<double> >& c) {
+void SplineFitter::create2DSpline(const vector<double>& x, const vector<double>& y, const vector<double>& values, bool periodic, vector<vector<double> >& c) {
     int xsize = x.size(), ysize = y.size();
-    if (xsize < 2 || ysize < 2)
-        throw OpenMMException("create2DNaturalSpline: must have at least two points along each axis");
+    // if (xsize < 2 || ysize < 2)
+    //     throw OpenMMException("create2DNaturalSpline: must have at least two points along each axis");
     if (values.size() != xsize*ysize)
         throw OpenMMException("create2DNaturalSpline: incorrect number of values");
     vector<double> d1(xsize*ysize), d2(xsize*ysize), d12(xsize*ysize);
@@ -203,7 +210,7 @@ void SplineFitter::create2DNaturalSpline(const vector<double>& x, const vector<d
     for (int i = 0; i < ysize; i++) {
         for (int j = 0; j < xsize; j++)
             t[j] = values[j+xsize*i];
-        SplineFitter::createNaturalSpline(x, t, deriv);
+        SplineFitter::createSpline(x, t, periodic, deriv);
         for (int j = 0; j < xsize; j++)
             d1[j+xsize*i] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[j]);
     }
@@ -215,7 +222,7 @@ void SplineFitter::create2DNaturalSpline(const vector<double>& x, const vector<d
     for (int i = 0; i < xsize; i++) {
         for (int j = 0; j < ysize; j++)
             t[j] = values[i+xsize*j];
-        SplineFitter::createNaturalSpline(y, t, deriv);
+        SplineFitter::createSpline(y, t, periodic, deriv);
         for (int j = 0; j < ysize; j++)
             d2[i+xsize*j] = SplineFitter::evaluateSplineDerivative(y, t, deriv, y[j]);
     }
@@ -227,7 +234,7 @@ void SplineFitter::create2DNaturalSpline(const vector<double>& x, const vector<d
     for (int i = 0; i < ysize; i++) {
         for (int j = 0; j < xsize; j++)
             t[j] = d2[j+xsize*i];
-        SplineFitter::createNaturalSpline(x, t, deriv);
+        SplineFitter::createSpline(x, t, periodic, deriv);
         for (int j = 0; j < xsize; j++)
             d12[j+xsize*i] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[j]);
     }
@@ -285,6 +292,10 @@ void SplineFitter::create2DNaturalSpline(const vector<double>& x, const vector<d
     }
 }
 
+void SplineFitter::create2DNaturalSpline(const vector<double>& x, const vector<double>& y, const vector<double>& values, vector<vector<double> >& c) {
+    SplineFitter::create2DSpline(x, y, values, false, c);
+}
+
 double SplineFitter::evaluate2DSpline(const vector<double>& x, const vector<double>& y, const vector<double>& values, const vector<vector<double> >& c, double u, double v) {
     int xsize = x.size();
     int ysize = y.size();
@@ -369,11 +380,11 @@ void SplineFitter::evaluate2DSplineDerivatives(const vector<double>& x, const ve
     dy /= deltay;
 }
 
-void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<double>& y, const vector<double>& z, const vector<double>& values, vector<vector<double> >& c) {
+void SplineFitter::create3DSpline(const vector<double>& x, const vector<double>& y, const vector<double>& z, const vector<double>& values, bool periodic, vector<vector<double> >& c) {
     int xsize = x.size(), ysize = y.size(), zsize = z.size();
     int xysize = xsize*ysize;
-    if (xsize < 2 || ysize < 2 || zsize < 2)
-        throw OpenMMException("create2DNaturalSpline: must have at least two points along each axis");
+    // if (xsize < 2 || ysize < 2 || zsize < 2)
+    //     throw OpenMMException("create2DNaturalSpline: must have at least two points along each axis");
     if (values.size() != xsize*ysize*zsize)
         throw OpenMMException("create2DNaturalSpline: incorrect number of values");
     vector<double> d1(xsize*ysize*zsize), d2(xsize*ysize*zsize), d3(xsize*ysize*zsize);
@@ -386,7 +397,7 @@ void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<d
         for (int j = 0; j < zsize; j++) {
             for (int k = 0; k < xsize; k++)
                 t[k] = values[k+xsize*i+xysize*j];
-            SplineFitter::createNaturalSpline(x, t, deriv);
+            SplineFitter::createSpline(x, t, periodic, deriv);
             for (int k = 0; k < xsize; k++)
                 d1[k+xsize*i+xysize*j] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[k]);
         }
@@ -400,7 +411,7 @@ void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<d
         for (int j = 0; j < zsize; j++) {
             for (int k = 0; k < ysize; k++)
                 t[k] = values[i+xsize*k+xysize*j];
-            SplineFitter::createNaturalSpline(y, t, deriv);
+            SplineFitter::createSpline(y, t, periodic, deriv);
             for (int k = 0; k < ysize; k++)
                 d2[i+xsize*k+xysize*j] = SplineFitter::evaluateSplineDerivative(y, t, deriv, y[k]);
         }
@@ -414,7 +425,7 @@ void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<d
         for (int j = 0; j < ysize; j++) {
             for (int k = 0; k < zsize; k++)
                 t[k] = values[i+xsize*j+xysize*k];
-            SplineFitter::createNaturalSpline(z, t, deriv);
+            SplineFitter::createSpline(z, t, periodic, deriv);
             for (int k = 0; k < zsize; k++)
                 d3[i+xsize*j+xysize*k] = SplineFitter::evaluateSplineDerivative(z, t, deriv, z[k]);
         }
@@ -428,7 +439,7 @@ void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<d
         for (int j = 0; j < zsize; j++) {
             for (int k = 0; k < xsize; k++)
                 t[k] = d2[k+xsize*i+xysize*j];
-            SplineFitter::createNaturalSpline(x, t, deriv);
+            SplineFitter::createSpline(x, t, periodic, deriv);
             for (int k = 0; k < xsize; k++)
                 d12[k+xsize*i+xysize*j] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[k]);
         }
@@ -442,7 +453,7 @@ void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<d
         for (int j = 0; j < xsize; j++) {
             for (int k = 0; k < ysize; k++)
                 t[k] = d3[j+xsize*k+xysize*i];
-            SplineFitter::createNaturalSpline(y, t, deriv);
+            SplineFitter::createSpline(y, t, periodic, deriv);
             for (int k = 0; k < ysize; k++)
                 d23[j+xsize*k+xysize*i] = SplineFitter::evaluateSplineDerivative(y, t, deriv, y[k]);
         }
@@ -456,7 +467,7 @@ void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<d
         for (int j = 0; j < ysize; j++) {
             for (int k = 0; k < zsize; k++)
                 t[k] = d1[i+xsize*j+xysize*k];
-            SplineFitter::createNaturalSpline(z, t, deriv);
+            SplineFitter::createSpline(z, t, periodic, deriv);
             for (int k = 0; k < zsize; k++)
                 d13[i+xsize*j+xysize*k] = SplineFitter::evaluateSplineDerivative(z, t, deriv, z[k]);
         }
@@ -470,7 +481,7 @@ void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<d
         for (int j = 0; j < zsize; j++) {
             for (int k = 0; k < xsize; k++)
                 t[k] = d23[k+xsize*i+xysize*j];
-            SplineFitter::createNaturalSpline(x, t, deriv);
+            SplineFitter::createSpline(x, t, periodic, deriv);
             for (int k = 0; k < xsize; k++)
                 d123[k+xsize*i+xysize*j] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[k]);
         }
@@ -599,6 +610,10 @@ void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<d
     }
 }
 
+void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<double>& y, const vector<double>& z, const vector<double>& values, vector<vector<double> >& c) {
+    SplineFitter::create3DSpline(x, y, z, values, false, c);
+}
+
 double SplineFitter::evaluate3DSpline(const vector<double>& x, const vector<double>& y, const vector<double>& z, const vector<double>& values, const vector<vector<double> >& c, double u, double v, double w) {
     int xsize = x.size();
     int ysize = y.size();