Added 3D spline functions to SplineFitter

openmmapi/include/openmm/internal/SplineFitter.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010 Stanford University and the Authors.           *
+ * Portions copyright (c) 2010-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -124,6 +124,52 @@ 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 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.
+     *
+     * @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 c      on exit, this contains the spline coefficients at each of the data points
+     */
+    static void create3DNaturalSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, const std::vector<double>& values, std::vector<std::vector<double> >& c);
+    /**
+     * Evaluate a 3D spline generated by one of the other methods in this class.
+     *
+     * @param x      the values of the first independent variable at the data points to interpolate
+     * @param y      the values of the second independent variable at the data points to interpolate
+     * @param z      the values of the third independent variable at the data points to interpolate
+     * @param values the values of the dependent variable at the data points to interpolate
+     * @param c      the vector of spline coefficients that was calculated by one of the other methods
+     * @param u      the value of the first independent variable at which to evaluate the spline
+     * @param v      the value of the second independent variable at which to evaluate the spline
+     * @param w      the value of the third independent variable at which to evaluate the spline
+     * @return the value of the spline at the specified point
+     */
+    static double evaluate3DSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, const std::vector<double>& values, const std::vector<std::vector<double> >& c, double u, double v, double w);
+    /**
+     * Evaluate the derivatives of a 3D spline generated by one of the other methods in this class.
+     *
+     * @param x      the values of the first independent variable at the data points to interpolate
+     * @param y      the values of the second independent variable at the data points to interpolate
+     * @param z      the values of the third independent variable at the data points to interpolate
+     * @param values the values of the dependent variable at the data points to interpolate
+     * @param c      the vector of spline coefficients that was calculated by one of the other methods
+     * @param u      the value of the first independent variable at which to evaluate the spline
+     * @param v      the value of the second independent variable at which to evaluate the spline
+     * @param w      the value of the third independent variable at which to evaluate the spline
+     * @param dx     on exit, the x derivative of the spline at the specified point
+     * @param dy     on exit, the y derivative of the spline at the specified point
+     * @param dz     on exit, the z derivative of the spline at the specified point
+     */
+    static void evaluate3DSplineDerivatives(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, const std::vector<double>& values, const std::vector<std::vector<double> >& c, double u, double v, double w, double& dx, double& dy, double &dz);
 private:
     static void solveTridiagonalMatrix(const std::vector<double>& a, const std::vector<double>& b, const std::vector<double>& c, const std::vector<double>& rhs, std::vector<double>& sol);
 };

tests/TestSplineFitter.cpp

@@ -106,8 +106,8 @@ void test2DSpline() {
         }
     for (int i = 0; i < 10; i++) {
         for (int j = 0; j < 10; j++) {
-            double s = x[0]+(i+1)*(x[xsize-1]-x[0])/11.0;
-            double t = y[0]+(j+1)*(y[ysize-1]-y[0])/11.0;
+            double s = x[0]+(i+1)*(x[xsize-1]-x[0])/12.0;
+            double t = y[0]+(j+1)*(y[ysize-1]-y[0])/12.0;
             double value = SplineFitter::evaluate2DSpline(x, y, f, c, s, t);
             ASSERT_EQUAL_TOL(sin(s)*cos(0.4*t), value, 0.02);
             double dx, dy;
@@ -118,11 +118,57 @@ void test2DSpline() {
     }
 }
 
+void test3DSpline() {
+    const int xsize = 8;
+    const int ysize = 9;
+    const int zsize = 10;
+    vector<double> x(xsize);
+    vector<double> y(ysize);
+    vector<double> z(zsize);
+    vector<double> f(xsize*ysize*zsize);
+    for (int i = 0; i < xsize; i++)
+        x[i] = 0.2*i+0.02*sin(0.4*double(i));
+    for (int i = 0; i < ysize; i++)
+        y[i] = 0.2*i+0.02*sin(0.45*double(i));
+    for (int i = 0; i < zsize; i++)
+        z[i] = 0.2*i+0.02*sin(0.5*double(i));
+    for (int i = 0; i < xsize; i++)
+        for (int j = 0; j < ysize; j++)
+            for (int k = 0; k < zsize; k++)
+                f[i+j*xsize+k*xsize*ysize] = sin(x[i])*cos(0.4*y[j])*(1+z[k]);
+    vector<vector<double> > c;
+    SplineFitter::create3DNaturalSpline(x, y, z, f, c);
+    for (int i = 0; i < xsize; i++)
+        for (int j = 0; j < ysize; j++) {
+            for (int k = 0; k < zsize; k++) {
+                double value = SplineFitter::evaluate3DSpline(x, y, z, f, c, x[i], y[j], z[k]);
+                ASSERT_EQUAL_TOL(f[i+j*xsize+k*xsize*ysize], value, 1e-6);
+            }
+        }
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            for (int k = 0; k < 10; k++) {
+                double s = x[0]+(i+1)*(x[xsize-1]-x[0])/12.0;
+                double t = y[0]+(j+1)*(y[ysize-1]-y[0])/12.0;
+                double u = z[0]+(k+1)*(z[zsize-1]-z[0])/12.0;
+                double value = SplineFitter::evaluate3DSpline(x, y, z, f, c, s, t, u);
+                ASSERT_EQUAL_TOL(sin(s)*cos(0.4*t)*(1+u), value, 0.02);
+                double dx, dy, dz;
+                SplineFitter::evaluate3DSplineDerivatives(x, y, z, f, c, s, t, u, dx, dy, dz);
+                ASSERT_EQUAL_TOL(cos(s)*cos(0.4*t)*(1+u), dx, 0.1);
+                ASSERT_EQUAL_TOL(-0.4*sin(s)*sin(0.4*t)*(1+u), dy, 0.1);
+                ASSERT_EQUAL_TOL(sin(s)*cos(0.4*t), dz, 0.1);
+            }
+        }
+    }
+}
+
 int main() {
     try {
         testNaturalSpline();
         testPeriodicSpline();
         test2DSpline();
+        test3DSpline();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;