Periodic 2D and 3D continuous functions in Reference platform

platforms/reference/include/ReferenceTabulatedFunction.h

@@ -78,6 +78,7 @@ private:
     const Continuous2DFunction& function;
     int xsize, ysize;
     double xmin, xmax, ymin, ymax;
+    bool periodic;
     std::vector<double> x, y, values;
     std::vector<std::vector<double> > c;
 };
@@ -97,6 +98,7 @@ private:
     const Continuous3DFunction& function;
     int xsize, ysize, zsize;
     double xmin, xmax, ymin, ymax, zmin, zmax;
+    bool periodic;
     std::vector<double> x, y, z, values;
     std::vector<std::vector<double> > c;
 };

platforms/reference/src/ReferenceTabulatedFunction.cpp

@@ -51,6 +51,12 @@ static int round(double x) {
 #include <cmath>
 #endif
 
+static double wrap(double t, double min, double max) {
+    double L = max - min;
+    double s = (t - min)/L;
+    return min + L*(s - floor(s));
+}
+
 using namespace OpenMM;
 using namespace std;
 using Lepton::CustomFunction;
@@ -75,21 +81,18 @@ extern "C" OPENMM_EXPORT CustomFunction* createReferenceTabulatedFunction(const
 }
 
 ReferenceContinuous1DFunction::ReferenceContinuous1DFunction(const Continuous1DFunction& function) : function(function) {
-    function.getFunctionParameters(values, min, max);
     periodic = function.getPeriodic();
+    function.getFunctionParameters(values, min, max);
     int numValues = values.size();
     x.resize(numValues);
     for (int i = 0; i < numValues; i++)
         x[i] = min+i*(max-min)/(numValues-1);
-    if (periodic)
-        SplineFitter::createPeriodicSpline(x, values, derivs);
-    else
-        SplineFitter::createNaturalSpline(x, values, derivs);
+    SplineFitter::createSpline(x, values, periodic, derivs);
 }
 
 ReferenceContinuous1DFunction::ReferenceContinuous1DFunction(const ReferenceContinuous1DFunction& other) : function(other.function) {
-    function.getFunctionParameters(values, min, max);
     periodic = function.getPeriodic();
+    function.getFunctionParameters(values, min, max);
     x = other.x;
     values = other.values;
     derivs = other.derivs;
@@ -100,30 +103,16 @@ int ReferenceContinuous1DFunction::getNumArguments() const {
 }
 
 double ReferenceContinuous1DFunction::evaluate(const double* arguments) const {
-    double t = arguments[0];
-    if (t < min || t > max) {
-        if (periodic) {
-            double L = max-min;
-            double s = (t-min)/L;
-            t = min + L*(s-floor(s));
-        }
-        else
+    double t = periodic ? wrap(arguments[0], min, max) : arguments[0];
+    if (t < min || t > max)
         return 0.0;
-    }
     return SplineFitter::evaluateSpline(x, values, derivs, t);
 }
 
 double ReferenceContinuous1DFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
-    double t = arguments[0];
-    if (t < min || t > max) {
-        if (periodic) {
-            double L = max-min;
-            double s = (t-min)/L;
-            t = min + L*(s-floor(s));
-        }
-        else
+    double t = periodic ? wrap(arguments[0], min, max) : arguments[0];
+    if (t < min || t > max)
         return 0.0;
-    }
     return SplineFitter::evaluateSplineDerivative(x, values, derivs, t);
 }
 
@@ -132,6 +121,7 @@ CustomFunction* ReferenceContinuous1DFunction::clone() const {
 }
 
 ReferenceContinuous2DFunction::ReferenceContinuous2DFunction(const Continuous2DFunction& function) : function(function) {
+    periodic = function.getPeriodic();
     function.getFunctionParameters(xsize, ysize, values, xmin, xmax, ymin, ymax);
     x.resize(xsize);
     y.resize(ysize);
@@ -139,10 +129,11 @@ ReferenceContinuous2DFunction::ReferenceContinuous2DFunction(const Continuous2DF
         x[i] = xmin+i*(xmax-xmin)/(xsize-1);
     for (int i = 0; i < ysize; i++)
         y[i] = ymin+i*(ymax-ymin)/(ysize-1);
-    SplineFitter::create2DNaturalSpline(x, y, values, c);
+    SplineFitter::create2DSpline(x, y, values, periodic, c);
 }
 
 ReferenceContinuous2DFunction::ReferenceContinuous2DFunction(const ReferenceContinuous2DFunction& other) : function(other.function) {
+    periodic = function.getPeriodic();
     function.getFunctionParameters(xsize, ysize, values, xmin, xmax, ymin, ymax);
     x = other.x;
     y = other.y;
@@ -155,20 +146,20 @@ int ReferenceContinuous2DFunction::getNumArguments() const {
 }
 
 double ReferenceContinuous2DFunction::evaluate(const double* arguments) const {
-    double u = arguments[0];
+    double u = periodic ? wrap(arguments[0], xmin, xmax) : arguments[0];
     if (u < xmin || u > xmax)
         return 0.0;
-    double v = arguments[1];
+    double v = periodic ? wrap(arguments[1], ymin, ymax) : arguments[1];
     if (v < ymin || v > ymax)
         return 0.0;
     return SplineFitter::evaluate2DSpline(x, y, values, c, u, v);
 }
 
 double ReferenceContinuous2DFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
-    double u = arguments[0];
+    double u = periodic ? wrap(arguments[0], xmin, xmax) : arguments[0];
     if (u < xmin || u > xmax)
         return 0.0;
-    double v = arguments[1];
+    double v = periodic ? wrap(arguments[1], ymin, ymax) : arguments[1];
     if (v < ymin || v > ymax)
         return 0.0;
     double dx, dy;
@@ -185,6 +176,7 @@ CustomFunction* ReferenceContinuous2DFunction::clone() const {
 }
 
 ReferenceContinuous3DFunction::ReferenceContinuous3DFunction(const Continuous3DFunction& function) : function(function) {
+    periodic = function.getPeriodic();
     function.getFunctionParameters(xsize, ysize, zsize, values, xmin, xmax, ymin, ymax, zmin, zmax);
     x.resize(xsize);
     y.resize(ysize);
@@ -195,10 +187,11 @@ ReferenceContinuous3DFunction::ReferenceContinuous3DFunction(const Continuous3DF
         y[i] = ymin+i*(ymax-ymin)/(ysize-1);
     for (int i = 0; i < zsize; i++)
         z[i] = zmin+i*(zmax-zmin)/(zsize-1);
-    SplineFitter::create3DNaturalSpline(x, y, z, values, c);
+    SplineFitter::create3DSpline(x, y, z, values, periodic, c);
 }
 
 ReferenceContinuous3DFunction::ReferenceContinuous3DFunction(const ReferenceContinuous3DFunction& other) : function(other.function) {
+    periodic = function.getPeriodic();
     function.getFunctionParameters(xsize, ysize, zsize, values, xmin, xmax, ymin, ymax, zmin, zmax);
     x = other.x;
     y = other.y;
@@ -212,26 +205,26 @@ int ReferenceContinuous3DFunction::getNumArguments() const {
 }
 
 double ReferenceContinuous3DFunction::evaluate(const double* arguments) const {
-    double u = arguments[0];
+    double u = periodic ? wrap(arguments[0], xmin, xmax) : arguments[0];
     if (u < xmin || u > xmax)
         return 0.0;
-    double v = arguments[1];
+    double v = periodic ? wrap(arguments[1], ymin, ymax) : arguments[1];
     if (v < ymin || v > ymax)
         return 0.0;
-    double w = arguments[2];
+    double w = periodic ? wrap(arguments[2], zmin, zmax) : arguments[2];
     if (w < zmin || w > zmax)
         return 0.0;
     return SplineFitter::evaluate3DSpline(x, y, z, values, c, u, v, w);
 }
 
 double ReferenceContinuous3DFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
-    double u = arguments[0];
+    double u = periodic ? wrap(arguments[0], xmin, xmax) : arguments[0];
     if (u < xmin || u > xmax)
         return 0.0;
-    double v = arguments[1];
+    double v = periodic ? wrap(arguments[1], ymin, ymax) : arguments[1];
     if (v < ymin || v > ymax)
         return 0.0;
-    double w = arguments[2];
+    double w = periodic ? wrap(arguments[2], zmin, zmax) : arguments[2];
     if (w < zmin || w > zmax)
         return 0.0;
     double dx, dy, dz;