Created reference implementation of Continuous2DFunction

openmmapi/include/openmm/TabulatedFunction.h

@@ -70,7 +70,7 @@ public:
      * Create a Continuous1DFunction f(x) based on a set of tabulated values.
      * 
      * @param values         the tabulated values of the function f(x) at uniformly spaced values of x between min
-     *                       and max.  A natural cubic spline is used to interpolated between the tabulated values.
+     *                       and max.  A natural cubic spline is used to interpolate between the tabulated values.
      *                       The function is assumed to be zero for x < min or x > max.
      * @param min            the value of x corresponding to the first element of values
      * @param max            the value of x corresponding to the last element of values
@@ -80,7 +80,7 @@ public:
      * Get the parameters for the tabulated function.
      *
      * @param values         the tabulated values of the function f(x) at uniformly spaced values of x between min
-     *                       and max.  A natural cubic spline is used to interpolated between the tabulated values.
+     *                       and max.  A natural cubic spline is used to interpolate between the tabulated values.
      *                       The function is assumed to be zero for x < min or x > max.
      * @param min            the value of x corresponding to the first element of values
      * @param max            the value of x corresponding to the last element of values
@@ -90,7 +90,7 @@ public:
      * Set the parameters for the tabulated function.
      *
      * @param values         the tabulated values of the function f(x) at uniformly spaced values of x between min
-     *                       and max.  A natural cubic spline is used to interpolated between the tabulated values.
+     *                       and max.  A natural cubic spline is used to interpolate between the tabulated values.
      *                       The function is assumed to be zero for x < min or x > max.
      * @param min            the value of x corresponding to the first element of values
      * @param max            the value of x corresponding to the last element of values
@@ -101,6 +101,62 @@ private:
     double min, max;
 };
 
+/**
+ * This is a TabulatedFunction that computes a continuous two dimensional function.
+ */
+class OPENMM_EXPORT Continuous2DFunction : public TabulatedFunction {
+public:
+    /**
+     * Create a Continuous2DFunction f(x,y) based on a set of tabulated values.
+     * 
+     * @param values     the tabulated values of the function f(x,y) at xsize uniformly spaced values of x between xmin
+     *                   and xmax, and ysize values of y between ymin and ymax.  A natural cubic spline is used to interpolate between the tabulated values.
+     *                   The function is assumed to be zero when x or y is outside its specified range.  The values should be ordered so that
+     *                   values[i+xsize*j] = f(x_i,y_j), where x_i is the i'th uniformly spaced value of x.  This must be of length xsize*ysize.
+     * @param xsize      the number of table elements along the x direction
+     * @param ysize      the number of table elements along the y direction
+     * @param xmin       the value of x corresponding to the first element of values
+     * @param xmax       the value of x corresponding to the last element of values
+     * @param ymin       the value of y corresponding to the first element of values
+     * @param ymax       the value of y corresponding to the last element of values
+     */
+    Continuous2DFunction(int xsize, int ysize, const std::vector<double>& values, double xmin, double xmax, double ymin, double ymax);
+    /**
+     * Get the parameters for the tabulated function.
+     *
+     * @param values     the tabulated values of the function f(x,y) at xsize uniformly spaced values of x between xmin
+     *                   and xmax, and ysize values of y between ymin and ymax.  A natural cubic spline is used to interpolate between the tabulated values.
+     *                   The function is assumed to be zero when x or y is outside its specified range.  The values should be ordered so that
+     *                   values[i+xsize*j] = f(x_i,y_j), where x_i is the i'th uniformly spaced value of x.  This must be of length xsize*ysize.
+     * @param xsize      the number of table elements along the x direction
+     * @param ysize      the number of table elements along the y direction
+     * @param xmin       the value of x corresponding to the first element of values
+     * @param xmax       the value of x corresponding to the last element of values
+     * @param ymin       the value of y corresponding to the first element of values
+     * @param ymax       the value of y corresponding to the last element of values
+     */
+    void getFunctionParameters(int& xsize, int& ysize, std::vector<double>& values, double& xmin, double& xmax, double& ymin, double& ymax) const;
+    /**
+     * Set the parameters for the tabulated function.
+     *
+     * @param values     the tabulated values of the function f(x,y) at xsize uniformly spaced values of x between xmin
+     *                   and xmax, and ysize values of y between ymin and ymax.  A natural cubic spline is used to interpolate between the tabulated values.
+     *                   The function is assumed to be zero when x or y is outside its specified range.  The values should be ordered so that
+     *                   values[i+xsize*j] = f(x_i,y_j), where x_i is the i'th uniformly spaced value of x.  This must be of length xsize*ysize.
+     * @param xsize      the number of table elements along the x direction
+     * @param ysize      the number of table elements along the y direction
+     * @param xmin       the value of x corresponding to the first element of values
+     * @param xmax       the value of x corresponding to the last element of values
+     * @param ymin       the value of y corresponding to the first element of values
+     * @param ymax       the value of y corresponding to the last element of values
+     */
+    void setFunctionParameters(int xsize, int ysize, const std::vector<double>& values, double xmin, double xmax, double ymin, double ymax);
+private:
+    std::vector<double> values;
+    int xsize, ysize;
+    double xmin, xmax, ymin, ymax;
+};
+
 /**
  * This is a TabulatedFunction that computes a discrete one dimensional function f(x).
  * To evaluate it, x is rounded to the nearest integer and the table element with that

openmmapi/src/SplineFitter.cpp

@@ -219,7 +219,7 @@ void SplineFitter::create2DNaturalSpline(const vector<double>& x, const vector<d
             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]);
+            d2[i+xsize*j] = SplineFitter::evaluateSplineDerivative(y, t, deriv, y[j]);
     }
 
     // Compute cross derivatives.

openmmapi/src/TabulatedFunction.cpp

@@ -35,7 +35,7 @@
 using namespace OpenMM;
 using namespace std;
 
-Continuous1DFunction::Continuous1DFunction(const std::vector<double>& values, double min, double max) {
+Continuous1DFunction::Continuous1DFunction(const vector<double>& values, double min, double max) {
     if (max <= min)
         throw OpenMMException("Continuous1DFunction: max <= min for a tabulated function.");
     if (values.size() < 2)
@@ -45,13 +45,13 @@ Continuous1DFunction::Continuous1DFunction(const std::vector<double>& values, do
     this->max = max;
 }
 
-void Continuous1DFunction::getFunctionParameters(std::vector<double>& values, double& min, double& max) const {
+void Continuous1DFunction::getFunctionParameters(vector<double>& values, double& min, double& max) const {
     values = this->values;
     min = this->min;
     max = this->max;
 }
 
-void Continuous1DFunction::setFunctionParameters(const std::vector<double>& values, double min, double max) {
+void Continuous1DFunction::setFunctionParameters(const vector<double>& values, double min, double max) {
     if (max <= min)
         throw OpenMMException("Continuous1DFunction: max <= min for a tabulated function.");
     if (values.size() < 2)
@@ -61,19 +61,65 @@ void Continuous1DFunction::setFunctionParameters(const std::vector<double>& valu
     this->max = max;
 }
 
-Discrete1DFunction::Discrete1DFunction(const std::vector<double>& values) {
+Continuous2DFunction::Continuous2DFunction(int xsize, int ysize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax) {
+    if (xsize < 2 || ysize < 2)
+        throw OpenMMException("Continuous2DFunction: must have at least two points along each axis");
+    if (values.size() != xsize*ysize)
+        throw OpenMMException("Continuous2DFunction: incorrect number of values");
+    if (xmax <= xmin)
+        throw OpenMMException("Continuous2DFunction: xmax <= xmin for a tabulated function.");
+    if (ymax <= ymin)
+        throw OpenMMException("Continuous2DFunction: ymax <= ymin for a tabulated function.");
+    this->values = values;
+    this->xsize = xsize;
+    this->ysize = ysize;
+    this->xmin = xmin;
+    this->xmax = xmax;
+    this->ymin = ymin;
+    this->ymax = ymax;
+}
+
+void Continuous2DFunction::getFunctionParameters(int& xsize, int& ysize, vector<double>& values, double& xmin, double& xmax, double& ymin, double& ymax) const {
+    values = this->values;
+    xsize = this->xsize;
+    ysize = this->ysize;
+    xmin = this->xmin;
+    xmax = this->xmax;
+    ymin = this->ymin;
+    ymax = this->ymax;
+}
+
+void Continuous2DFunction::setFunctionParameters(int xsize, int ysize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax) {
+    if (xsize < 2 || ysize < 2)
+        throw OpenMMException("Continuous2DFunction: must have at least two points along each axis");
+    if (values.size() != xsize*ysize)
+        throw OpenMMException("Continuous2DFunction: incorrect number of values");
+    if (xmax <= xmin)
+        throw OpenMMException("Continuous2DFunction: xmax <= xmin for a tabulated function.");
+    if (ymax <= ymin)
+        throw OpenMMException("Continuous2DFunction: ymax <= ymin for a tabulated function.");
+    this->values = values;
+    this->xsize = xsize;
+    this->ysize = ysize;
+    this->xmin = xmin;
+    this->xmax = xmax;
+    this->ymin = ymin;
+    this->ymax = ymax;
+}
+
+Discrete1DFunction::Discrete1DFunction(const vector<double>& values) {
     this->values = values;
 }
 
-void Discrete1DFunction::getFunctionParameters(std::vector<double>& values) const {
+void Discrete1DFunction::getFunctionParameters(vector<double>& values) const {
     values = this->values;
 }
 
-void Discrete1DFunction::setFunctionParameters(const std::vector<double>& values) {
+void Discrete1DFunction::setFunctionParameters(const vector<double>& values) {
     this->values = values;
 }
 
-Discrete2DFunction::Discrete2DFunction(int xsize, int ysize, const std::vector<double>& values) {
+Discrete2DFunction::Discrete2DFunction(int xsize, int ysize, const vector<double>& values) {
     if (values.size() != xsize*ysize)
         throw OpenMMException("Discrete2DFunction: incorrect number of values");
     this->xsize = xsize;
@@ -81,13 +127,13 @@ Discrete2DFunction::Discrete2DFunction(int xsize, int ysize, const std::vector<d
     this->values = values;
 }
 
-void Discrete2DFunction::getFunctionParameters(int& xsize, int& ysize, std::vector<double>& values) const {
+void Discrete2DFunction::getFunctionParameters(int& xsize, int& ysize, vector<double>& values) const {
     xsize = this->xsize;
     ysize = this->ysize;
     values = this->values;
 }
 
-void Discrete2DFunction::setFunctionParameters(int xsize, int ysize, const std::vector<double>& values) {
+void Discrete2DFunction::setFunctionParameters(int xsize, int ysize, const vector<double>& values) {
     if (values.size() != xsize*ysize)
         throw OpenMMException("Discrete2DFunction: incorrect number of values");
     this->xsize = xsize;
@@ -95,7 +141,7 @@ void Discrete2DFunction::setFunctionParameters(int xsize, int ysize, const std::
     this->values = values;
 }
 
-Discrete3DFunction::Discrete3DFunction(int xsize, int ysize, int zsize, const std::vector<double>& values) {
+Discrete3DFunction::Discrete3DFunction(int xsize, int ysize, int zsize, const vector<double>& values) {
     if (values.size() != xsize*ysize*zsize)
         throw OpenMMException("Discrete3DFunction: incorrect number of values");
     this->xsize = xsize;
@@ -104,14 +150,14 @@ Discrete3DFunction::Discrete3DFunction(int xsize, int ysize, int zsize, const st
     this->values = values;
 }
 
-void Discrete3DFunction::getFunctionParameters(int& xsize, int& ysize, int& zsize, std::vector<double>& values) const {
+void Discrete3DFunction::getFunctionParameters(int& xsize, int& ysize, int& zsize, vector<double>& values) const {
     xsize = this->xsize;
     ysize = this->ysize;
     zsize = this->zsize;
     values = this->values;
 }
 
-void Discrete3DFunction::setFunctionParameters(int xsize, int ysize, int zsize, const std::vector<double>& values) {
+void Discrete3DFunction::setFunctionParameters(int xsize, int ysize, int zsize, const vector<double>& values) {
     if (values.size() != xsize*ysize*zsize)
         throw OpenMMException("Discrete3DFunction: incorrect number of values");
     this->xsize = xsize;

platforms/reference/include/ReferenceTabulatedFunction.h

@@ -60,6 +60,24 @@ private:
     std::vector<double> x, values, derivs;
 };
 
+/**
+ * This class adapts a Continuous2DFunction into a Lepton::CustomFunction.
+ */
+class OPENMM_EXPORT ReferenceContinuous2DFunction : public Lepton::CustomFunction {
+public:
+    ReferenceContinuous2DFunction(const Continuous2DFunction& function);
+    int getNumArguments() const;
+    double evaluate(const double* arguments) const;
+    double evaluateDerivative(const double* arguments, const int* derivOrder) const;
+    CustomFunction* clone() const;
+private:
+    const Continuous2DFunction& function;
+    int xsize, ysize;
+    double xmin, xmax, ymin, ymax;
+    std::vector<double> x, y, values;
+    std::vector<std::vector<double> > c;
+};
+
 /**
  * This class adapts a Discrete1DFunction into a Lepton::CustomFunction.
  */

platforms/reference/src/ReferenceTabulatedFunction.cpp

@@ -41,6 +41,8 @@ using Lepton::CustomFunction;
 extern "C" CustomFunction* createReferenceTabulatedFunction(const TabulatedFunction& function) {
     if (dynamic_cast<const Continuous1DFunction*>(&function) != NULL)
         return new ReferenceContinuous1DFunction(dynamic_cast<const Continuous1DFunction&>(function));
+    if (dynamic_cast<const Continuous2DFunction*>(&function) != NULL)
+        return new ReferenceContinuous2DFunction(dynamic_cast<const Continuous2DFunction&>(function));
     if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL)
         return new ReferenceDiscrete1DFunction(dynamic_cast<const Discrete1DFunction&>(function));
     if (dynamic_cast<const Discrete2DFunction*>(&function) != NULL)
@@ -81,6 +83,51 @@ CustomFunction* ReferenceContinuous1DFunction::clone() const {
     return new ReferenceContinuous1DFunction(function);
 }
 
+ReferenceContinuous2DFunction::ReferenceContinuous2DFunction(const Continuous2DFunction& function) : function(function) {
+    function.getFunctionParameters(xsize, ysize, values, xmin, xmax, ymin, ymax);
+    x.resize(xsize);
+    y.resize(ysize);
+    for (int i = 0; i < xsize; i++)
+        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);
+}
+
+int ReferenceContinuous2DFunction::getNumArguments() const {
+    return 2;
+}
+
+double ReferenceContinuous2DFunction::evaluate(const double* arguments) const {
+    double u = arguments[0];
+    if (u < xmin || u > xmax)
+        return 0.0;
+    double v = 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];
+    if (u < xmin || u > xmax)
+        return 0.0;
+    double v = arguments[1];
+    if (v < ymin || v > ymax)
+        return 0.0;
+    double dx, dy;
+    SplineFitter::evaluate2DSplineDerivatives(x, y, values, c, u, v, dx, dy);
+    if (derivOrder[0] == 1 && derivOrder[1] == 0)
+        return dx;
+    if (derivOrder[0] == 0 && derivOrder[1] == 1)
+        return dy;
+    throw OpenMMException("ReferenceContinuous2DFunction: Unsupported derivative order");
+}
+
+CustomFunction* ReferenceContinuous2DFunction::clone() const {
+    return new ReferenceContinuous2DFunction(function);
+}
+
 ReferenceDiscrete1DFunction::ReferenceDiscrete1DFunction(const Discrete1DFunction& function) : function(function) {
     function.getFunctionParameters(values);
 }

platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp

@@ -244,7 +244,6 @@ void testContinuous1DFunction() {
     Context context(system, integrator, platform);
     vector<Vec3> positions(2);
     positions[0] = Vec3(0, 0, 0);
-    double tol = 0.01;
     for (int i = 1; i < 30; i++) {
         double x = (7.0/30.0)*i;
         positions[1] = Vec3(x, 0, 0);
@@ -267,6 +266,55 @@ void testContinuous1DFunction() {
     }
 }
 
+void testContinuous2DFunction() {
+    const int xsize = 20;
+    const int ysize = 21;
+    const double xmin = 0.4;
+    const double xmax = 1.5;
+    const double ymin = 0.0;
+    const double ymax = 2.1;
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* forceField = new CustomNonbondedForce("fn(r,a)+1");
+    forceField->addGlobalParameter("a", 0.0);
+    forceField->addParticle(vector<double>());
+    forceField->addParticle(vector<double>());
+    vector<double> table(xsize*ysize);
+    for (int i = 0; i < xsize; i++) {
+        for (int j = 0; j < ysize; j++) {
+            double x = xmin + i*(xmax-xmin)/xsize;
+            double y = ymin + j*(ymax-ymin)/ysize;
+            table[i+xsize*j] = sin(0.25*x)*cos(0.33*y);
+        }
+    }
+    forceField->addFunction("fn", new Continuous2DFunction(xsize, ysize, table, xmin, xmax, ymin, ymax));
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    for (double x = xmin-0.15; x < xmax+0.2; x += 0.1) {
+        for (double y = ymin-0.15; y < ymax+0.2; y += 0.1) {
+            positions[1] = Vec3(x, 0, 0);
+            context.setParameter("a", y);
+            context.setPositions(positions);
+            State state = context.getState(State::Forces | State::Energy);
+            const vector<Vec3>& forces = state.getForces();
+            double energy = 1;
+            double force = 0;
+            if (x >= xmin && x <= xmax && y >= ymin && y <= ymax) {
+                energy = sin(0.25*x)*cos(0.33*y)+1.0;
+                force = -0.25*cos(0.25*x)*cos(0.33*y);
+            }
+            ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], 0.1);
+            ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], 0.1);
+            ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.02);
+        }
+    }
+}
+
 void testDiscrete1DFunction() {
     ReferencePlatform platform;
     System system;
@@ -757,6 +805,7 @@ int main() {
         testCutoff();
         testPeriodic();
         testContinuous1DFunction();
+        testContinuous2DFunction();
         testDiscrete1DFunction();
         testDiscrete2DFunction();
         testDiscrete3DFunction();

tests/TestSplineFitter.cpp

@@ -103,19 +103,17 @@ void test2DSpline() {
         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);
+            ASSERT_EQUAL_TOL(sin(s)*cos(0.4*t), value, 0.02);
             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);
+            ASSERT_EQUAL_TOL(cos(s)*cos(0.4*t), dx, 0.05);
+            ASSERT_EQUAL_TOL(-0.4*sin(s)*sin(0.4*t), dy, 0.05);
         }
     }
 }