Merge pull request #2629 from craabreu/periodic_spline

Bug fix in periodic spline filter (solves #2627)

openmmapi/src/SplineFitter.cpp

@@ -86,8 +86,8 @@ void SplineFitter::createPeriodicSpline(const vector<double>& x, const vector<do
 
     // Create the system of equations to solve.
 
-    vector<double> a(n), b(n), c(n), rhs(n);
-    a[0] = 0.0;
+    vector<double> a(n-1), b(n-1), c(n-1), rhs(n-1);
+    a[0] = x[n-1]-x[n-2];
     b[0] = 2.0*(x[1]-x[0]+x[n-1]-x[n-2]);
     c[0] = x[1]-x[0];
     rhs[0] = 6.0*((y[1]-y[0])/(x[1]-x[0]) - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
@@ -97,17 +97,14 @@ void SplineFitter::createPeriodicSpline(const vector<double>& x, const vector<do
         c[i] = x[i+1]-x[i];
         rhs[i] = 6.0*((y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]));
     }
-    a[n-1] = 0.0;
-    b[n-1] = 1.0;
-    c[n-1] = 0.0;
-    rhs[n-1] = 0.0;
-    double beta = x[n-1]-x[n-2];
-    double alpha = -1.0;
+    double beta = a[0];
+    double alpha = c[n-2];
     double gamma = -b[0];
 
     // This is a cyclic tridiagonal matrix.  We solve it using the Sherman-Morrison method,
     // which involves solving two tridiagonal systems.
 
+    n--;
     b[0] -= gamma;
     b[n-1] -= alpha*beta/gamma;
     solveTridiagonalMatrix(a, b, c, rhs, deriv);
@@ -118,6 +115,7 @@ void SplineFitter::createPeriodicSpline(const vector<double>& x, const vector<do
     double scale = (deriv[0]+beta*deriv[n-1]/gamma)/(1.0+z[0]+beta*z[n-1]/gamma);
     for (int i = 0; i < n; i++)
         deriv[i] -= scale*z[i];
+    deriv[n] = deriv[0];
 }
 
 double SplineFitter::evaluateSpline(const vector<double>& x, const vector<double>& y, const vector<double>& deriv, double t) {

tests/TestSplineFitter.cpp

@@ -82,6 +82,22 @@ void testPeriodicSpline() {
         ASSERT_EQUAL_TOL(sin((double)i), SplineFitter::evaluateSpline(x, y, deriv, i), 0.05);
         ASSERT_EQUAL_TOL(cos((double)i), SplineFitter::evaluateSplineDerivative(x, y, deriv, i), 0.05);
     }
+    for (unsigned int i = 0; i < x.size(); i++)
+        x[i] = i/(x.size()-1.0);
+    double ya[] = {15.579, 16.235, 17.325, 18.741, 20.454, 22.517, 24.944, 27.554, 29.942, 31.657,
+                   32.486, 32.612, 32.494, 32.532, 32.785, 32.917, 32.402, 30.842, 28.229, 24.989,
+                   21.762, 19.074, 17.147, 15.970, 15.467, 15.579};
+    // scipy.interpolate.CubicSpline solution:
+    double sol[] = { 345.520,  271.991,  194.015,  174.449, 221.940, 250.291, 141.895, -131.620,
+                    -447.916, -600.465, -472.723, -144.892, 137.290, 180.733, -53.971, -418.600,
+                    -697.879, -708.635, -416.330,   22.704, 374.262, 501.498, 473.496,  417.019,
+                     385.928,  345.520};
+    y.assign(begin(ya), end(ya));
+    SplineFitter::createPeriodicSpline(x, y, deriv);
+    ASSERT_EQUAL_TOL(SplineFitter::evaluateSplineDerivative(x, y, deriv, x[0]),
+                     SplineFitter::evaluateSplineDerivative(x, y, deriv, x[x.size()-1]), 1e-6);
+    for (int i = 0; i < x.size(); i++)
+        ASSERT_EQUAL_TOL(deriv[i], sol[i], 1e-3);
 }
 
 void test2DSpline() {
@@ -177,4 +193,3 @@ int main() {
     cout << "Done" << endl;
     return 0;
 }
-