Added 2D spline functions to SplineFitter

openmmapi/include/openmm/internal/SplineFitter.h

@@ -67,7 +67,7 @@ public:
      */
     static void createPeriodicSpline(const std::vector<double>& x, const std::vector<double>& y, std::vector<double>& deriv);
     /**
-     * Evaluate a spline generated by one of the other methods in this class.
+     * Evaluate a 1D spline generated by one of the other methods in this class.
      *
      * @param x     the values of the independent variable at the data points to interpolate
      * @param y     the values of the dependent variable at the data points to interpolate
@@ -77,7 +77,7 @@ public:
      */
     static double evaluateSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& deriv, double t);
     /**
-     * Evaluate the derivative of a spline generated by one of the other methods in this class.
+     * Evaluate the derivative of a 1D spline generated by one of the other methods in this class.
      *
      * @param x     the values of the independent variable at the data points to interpolate
      * @param y     the values of the dependent variable at the data points to interpolate
@@ -86,6 +86,44 @@ 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 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.
+     *
+     * @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 c      on exit, this contains the spline coefficients at each of the data points
+     */
+    static void create2DNaturalSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& values, std::vector<std::vector<double> >& c);
+    /**
+     * Evaluate a 2D 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 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
+     * @return the value of the spline at the specified point
+     */
+    static double evaluate2DSpline(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);
+    /**
+     * Evaluate the derivatives of a 2D 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 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 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
+     */
+    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);
 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);
 };

openmmapi/src/SplineFitter.cpp

@@ -190,3 +190,183 @@ void SplineFitter::solveTridiagonalMatrix(const vector<double>& a, const vector<
     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) {
+    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 (values.size() != xsize*ysize)
+        throw OpenMMException("create2DNaturalSpline: incorrect number of values");
+    vector<double> d1(xsize*ysize), d2(xsize*ysize), d12(xsize*ysize);
+    vector<double> t(xsize), deriv(xsize);
+
+    // Compute derivatives with respect to x.
+
+    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);
+        for (int j = 0; j < xsize; j++)
+            d1[j+xsize*i] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[j]);
+    }
+
+    // Compute derivatives with respect to y.
+
+    t.resize(ysize);
+    deriv.resize(ysize);
+    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);
+        for (int j = 0; j < ysize; j++)
+            d2[i+xsize*j] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[j]);
+    }
+
+    // Compute cross derivatives.
+
+    t.resize(xsize);
+    deriv.resize(xsize);
+    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);
+        for (int j = 0; j < xsize; j++)
+            d12[j+xsize*i] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[j]);
+    }
+
+    // Now compute the coefficients.
+
+    const int wt[] = {
+        1, 0, -3, 2, 0, 0, 0, 0, -3, 0, 9, -6, 2, 0, -6, 4,
+        0, 0, 0, 0, 0, 0, 0, 0, 3, 0, -9, 6, -2, 0, 6, -4,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -6, 0, 0, -6, 4,
+        0, 0, 3, -2, 0, 0, 0, 0, 0, 0, -9, 6, 0, 0, 6, -4,
+        0, 0, 0, 0, 1, 0, -3, 2, -2, 0, 6, -4, 1, 0, -3, 2,
+        0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 3, -2, 1, 0, -3, 2,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2, 0, 0, 3, -2,
+        0, 0, 0, 0, 0, 0, 3, -2, 0, 0, -6, 4, 0, 0, 3, -2,
+        0, 1, -2, 1, 0, 0, 0, 0, 0, -3, 6, -3, 0, 2, -4, 2,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -6, 3, 0, -2, 4, -2,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 2, -2,
+        0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 3, -3, 0, 0, -2, 2,
+        0, 0, 0, 0, 0, 1, -2, 1, 0, -2, 4, -2, 0, 1, -2, 1,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 1, -2, 1,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, -1, 1,
+        0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 2, -2, 0, 0, -1, 1
+    };
+    vector<double> rhs(16);
+    c.resize(xsize*ysize);
+    for (int i = 0; i < xsize; i++) {
+        for (int j = 0; j < ysize; j++) {
+            // Compute the 16 coefficients for patch (i, j).
+
+            int nexti = i+1;
+            int nextj = j+1;
+            double deltax = x[nexti]-x[i];
+            double deltay = y[nextj]-y[j];
+            double e[] = {values[i+j*xsize], values[nexti+j*xsize], values[nexti+nextj*xsize], values[i+nextj*xsize]};
+            double e1[] = {d1[i+j*xsize], d1[nexti+j*xsize], d1[nexti+nextj*xsize], d1[i+nextj*xsize]};
+            double e2[] = {d2[i+j*xsize], d2[nexti+j*xsize], d2[nexti+nextj*xsize], d2[i+nextj*xsize]};
+            double e12[] = {d12[i+j*xsize], d12[nexti+j*xsize], d12[nexti+nextj*xsize], d12[i+nextj*xsize]};
+
+            for (int k = 0; k < 4; k++) {
+                rhs[k] = e[k];
+                rhs[k+4] = e1[k]*deltax;
+                rhs[k+8] = e2[k]*deltay;
+                rhs[k+12] = e12[k]*deltax*deltay;
+            }
+            vector<double>& coeff = c[i+j*xsize];
+            coeff.resize(16);
+            for (int k = 0; k < 16; k++) {
+                double sum = 0.0;
+                for (int m = 0; m < 16; m++)
+                    sum += wt[k+16*m]*rhs[m];
+                coeff[k] = sum;
+            }
+        }
+    }
+}
+
+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();
+    if (u < x[0] || u > x[xsize-1] || v < y[0] || v > y[ysize-1])
+        throw OpenMMException("evaluate2DSpline: specified point is outside the range defined by the spline");
+
+    // Perform a binary search to identify the interval containing the point to evaluate.
+
+    int lowerx = 0;
+    int upperx = xsize-1;
+    while (upperx-lowerx > 1) {
+        int middle = (upperx+lowerx)/2;
+        if (x[middle] > u)
+            upperx = middle;
+        else
+            lowerx = middle;
+    }
+    int lowery = 0;
+    int uppery = ysize-1;
+    while (uppery-lowery > 1) {
+        int middle = (uppery+lowery)/2;
+        if (y[middle] > v)
+            uppery = middle;
+        else
+            lowery = middle;
+    }
+    double deltax = x[upperx]-x[lowerx];
+    double deltay = y[uppery]-y[lowery];
+    double da = (u-x[lowerx])/deltax;
+    double db = (v-y[lowery])/deltay;
+    const vector<double>& coeff = c[lowerx+xsize*lowery];
+
+    // Evaluate the spline to determine the value and gradients.
+
+    double value = 0;
+    for (int i = 3; i >= 0; i--)
+        value = da*value + ((coeff[i*4+3]*db + coeff[i*4+2])*db + coeff[i*4+1])*db + coeff[i*4+0];
+    return value;
+}
+
+void SplineFitter::evaluate2DSplineDerivatives(const vector<double>& x, const vector<double>& y, const vector<double>& values, const vector<vector<double> >& c, double u, double v, double& dx, double &dy) {
+    int xsize = x.size();
+    int ysize = y.size();
+    if (u < x[0] || u > x[xsize-1] || v < y[0] || v > y[ysize-1])
+        throw OpenMMException("evaluate2DSplineDerivatives: specified point is outside the range defined by the spline");
+
+    // Perform a binary search to identify the interval containing the point to evaluate.
+
+    int lowerx = 0;
+    int upperx = xsize-1;
+    while (upperx-lowerx > 1) {
+        int middle = (upperx+lowerx)/2;
+        if (x[middle] > u)
+            upperx = middle;
+        else
+            lowerx = middle;
+    }
+    int lowery = 0;
+    int uppery = ysize-1;
+    while (uppery-lowery > 1) {
+        int middle = (uppery+lowery)/2;
+        if (y[middle] > v)
+            uppery = middle;
+        else
+            lowery = middle;
+    }
+    double deltax = x[upperx]-x[lowerx];
+    double deltay = y[uppery]-y[lowery];
+    double da = (u-x[lowerx])/deltax;
+    double db = (v-y[lowery])/deltay;
+    const vector<double>& coeff = c[lowerx+xsize*lowery];
+
+    // Evaluate the spline to determine the value and gradients.
+
+    dx = 0;
+    dy = 0;
+    for (int i = 3; i >= 0; i--) {
+        dx = db*dx + (3.0*coeff[i+3*4]*da + 2.0*coeff[i+2*4])*da + coeff[i+1*4];
+        dy = da*dy + (3.0*coeff[i*4+3]*db + 2.0*coeff[i*4+2])*db + coeff[i*4+1];
+    }
+    dx /= deltax;
+    dy /= deltay;
+}

tests/TestSplineFitter.cpp

@@ -84,10 +84,47 @@ void testPeriodicSpline() {
     }
 }
 
+void test2DSpline() {
+    const int xsize = 15;
+    const int ysize = 17;
+    vector<double> x(xsize);
+    vector<double> y(ysize);
+    vector<double> f(xsize*ysize);
+    for (int i = 0; i < xsize; i++)
+        x[i] = 0.5*i+0.1*sin(double(i));
+    for (int i = 0; i < ysize; i++)
+        y[i] = 0.6*i+0.1*sin(double(i));
+    for (int i = 0; i < xsize; i++)
+        for (int j = 0; j < ysize; j++)
+            f[i+j*xsize] = sin(x[i])*cos(0.4*y[j]);
+    vector<vector<double> > c;
+    SplineFitter::create2DNaturalSpline(x, y, f, c);
+    for (int i = 0; i < xsize; i++)
+        for (int j = 0; j < ysize; j++) {
+            double value = SplineFitter::evaluate2DSpline(x, y, f, c, x[i], y[j]);
+            ASSERT_EQUAL_TOL(f[i+j*xsize], value, 1e-6);
+            double dx, dy;
+            SplineFitter::evaluate2DSplineDerivatives(x, y, f, c, x[i], y[j], dx, dy);
+        }
+    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 value = SplineFitter::evaluate2DSpline(x, y, f, c, s, t);
+            ASSERT_EQUAL_TOL(sin(s)*cos(0.4*t), value, 0.05);
+            double dx, dy;
+            SplineFitter::evaluate2DSplineDerivatives(x, y, f, c, s, t, dx, dy);
+            ASSERT_EQUAL_TOL(cos(s)*cos(0.4*t), dx, 0.1);
+            ASSERT_EQUAL_TOL(-0.4*sin(s)*sin(0.4*t), dy, 0.1);
+        }
+    }
+}
+
 int main() {
     try {
         testNaturalSpline();
         testPeriodicSpline();
+        test2DSpline();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;