Merge pull request #904 from peastman/select

Implemented select() in custom expressions

docs-source/usersguide/theory.rst

@@ -1040,9 +1040,9 @@ The following operators are supported: + (add), - (subtract), * (multiply), /
 
 The following standard functions are supported: sqrt, exp, log, sin, cos, sec,
 csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs,
-floor, ceil, step, delta. step(x) = 0 if x < 0, 1 otherwise.  delta(x) = 1 if x is 0, 0
-otherwise.  Some custom forces allow additional functions to be defined from
-tabulated values.
+floor, ceil, step, delta, select. step(x) = 0 if x < 0, 1 otherwise.
+delta(x) = 1 if x is 0, 0 otherwise.  select(x,y,z) = z if x = 0, y otherwise.
+Some custom forces allow additional functions to be defined from tabulated values.
 
 Numbers may be given in either decimal or exponential form.  All of the
 following are valid numbers: 5, -3.1, 1e6, and 3.12e-2.

libraries/lepton/include/lepton/Operation.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) 2009-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -64,7 +64,7 @@ public:
      */
     enum Id {CONSTANT, VARIABLE, CUSTOM, ADD, SUBTRACT, MULTIPLY, DIVIDE, POWER, NEGATE, SQRT, EXP, LOG,
              SIN, COS, SEC, CSC, TAN, COT, ASIN, ACOS, ATAN, SINH, COSH, TANH, ERF, ERFC, STEP, DELTA, SQUARE, CUBE, RECIPROCAL,
-             ADD_CONSTANT, MULTIPLY_CONSTANT, POWER_CONSTANT, MIN, MAX, ABS, FLOOR, CEIL};
+             ADD_CONSTANT, MULTIPLY_CONSTANT, POWER_CONSTANT, MIN, MAX, ABS, FLOOR, CEIL, SELECT};
     /**
      * Get the name of this Operation.
      */
@@ -155,6 +155,7 @@ public:
     class Abs;
     class Floor;
     class Ceil;
+    class Select;
 };
 
 class LEPTON_EXPORT Operation::Constant : public Operation {
@@ -1137,6 +1138,28 @@ public:
     ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const;
 };
 
+class LEPTON_EXPORT Operation::Select : public Operation {
+public:
+    Select() {
+    }
+    std::string getName() const {
+        return "select";
+    }
+    Id getId() const {
+        return SELECT;
+    }
+    int getNumArguments() const {
+        return 3;
+    }
+    Operation* clone() const {
+        return new Select();
+    }
+    double evaluate(double* args, const std::map<std::string, double>& variables) const {
+        return (args[0] != 0.0 ? args[1] : args[2]);
+    }
+    ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const;
+};
+
 } // namespace Lepton
 
 #endif /*LEPTON_OPERATION_H_*/

libraries/lepton/src/Operation.cpp

@@ -7,7 +7,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2009 Stanford University and the Authors.           *
+ * Portions copyright (c) 2009-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -325,3 +325,11 @@ ExpressionTreeNode Operation::Floor::differentiate(const std::vector<ExpressionT
 ExpressionTreeNode Operation::Ceil::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
     return ExpressionTreeNode(new Operation::Constant(0.0));
 }
+
+ExpressionTreeNode Operation::Select::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
+    vector<ExpressionTreeNode> derivChildren;
+    derivChildren.push_back(children[0]);
+    derivChildren.push_back(childDerivs[1]);
+    derivChildren.push_back(childDerivs[2]);
+    return ExpressionTreeNode(new Operation::Select(), derivChildren);
+}

libraries/lepton/src/Parser.cpp

@@ -328,6 +328,7 @@ Operation* Parser::getFunctionOperation(const std::string& name, const map<strin
         opMap["abs"] = Operation::ABS;
         opMap["floor"] = Operation::FLOOR;
         opMap["ceil"] = Operation::CEIL;
+        opMap["select"] = Operation::SELECT;
     }
     string trimmed = name.substr(0, name.size()-1);
 
@@ -397,6 +398,8 @@ Operation* Parser::getFunctionOperation(const std::string& name, const map<strin
             return new Operation::Floor();
         case Operation::CEIL:
             return new Operation::Ceil();
+        case Operation::SELECT:
+            return new Operation::Select();
         default:
             throw Exception("unknown function");
     }

openmmapi/include/openmm/CustomAngleForce.h

@@ -65,8 +65,9 @@ namespace OpenMM {
  * </pre></tt>
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
- * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta.  All trigonometric functions
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
  * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
+ * select(x,y,z) = z if x = 0, y otherwise.
  */
 
 class OPENMM_EXPORT CustomAngleForce : public Force {

openmmapi/include/openmm/CustomBondForce.h

@@ -65,8 +65,9 @@ namespace OpenMM {
  * </pre></tt>
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
- * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta.  All trigonometric functions
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
  * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
+ * select(x,y,z) = z if x = 0, y otherwise.
  */
 
 class OPENMM_EXPORT CustomBondForce : public Force {

openmmapi/include/openmm/CustomCompoundBondForce.h

@@ -89,8 +89,9 @@ namespace OpenMM {
  * </pre></tt>
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
- * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta.  All trigonometric functions
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
  * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
+ * select(x,y,z) = z if x = 0, y otherwise.
  *
  * In addition, you can call addTabulatedFunction() to define a new function based on tabulated values.  You specify the function by
  * creating a TabulatedFunction object.  That function can then appear in the expression.

openmmapi/include/openmm/CustomExternalForce.h

@@ -68,8 +68,9 @@ namespace OpenMM {
  * </pre></tt>
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
- * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta.  All trigonometric functions
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
  * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
+ * select(x,y,z) = z if x = 0, y otherwise.
  */
 
 class OPENMM_EXPORT CustomExternalForce : public Force {

openmmapi/include/openmm/CustomGBForce.h

@@ -129,9 +129,9 @@ namespace OpenMM {
  * particular piece of the computation.
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
- * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta.  All trigonometric functions
- * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.  In expressions for
- * particle pair calculations, the names of per-particle parameters and computed values
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
+ * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
+ * select(x,y,z) = z if x = 0, y otherwise.  In expressions for particle pair calculations, the names of per-particle parameters and computed values
  * have the suffix "1" or "2" appended to them to indicate the values for the two interacting particles.  As seen in the above example,
  * an expression may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
  *

openmmapi/include/openmm/CustomHbondForce.h

@@ -89,8 +89,9 @@ namespace OpenMM {
  * </pre></tt>
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
- * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta.  All trigonometric functions
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
  * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
+ * select(x,y,z) = z if x = 0, y otherwise.
  *
  * In addition, you can call addTabulatedFunction() to define a new function based on tabulated values.  You specify the function by
  * creating a TabulatedFunction object.  That function can then appear in the expression.

openmmapi/include/openmm/CustomIntegrator.h

@@ -204,9 +204,9 @@ namespace OpenMM {
  * the force from a single force group, or a random number.
  * 
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
- * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta.  All trigonometric functions
- * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.  An expression
- * may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
+ * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
+ * select(x,y,z) = z if x = 0, y otherwise.  An expression may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
  */
 
 class OPENMM_EXPORT CustomIntegrator : public Integrator {

openmmapi/include/openmm/CustomManyParticleForce.h

@@ -149,10 +149,10 @@ namespace OpenMM {
  * still only be evaluated once for each triplet, so it must still be symmetric with respect to p2 and p3.
  * 
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
- * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta.  All trigonometric functions
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
  * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
- * The names of per-particle parameters have the suffix "1", "2", etc. appended to them to indicate the values for the multiple interacting particles.
- * For example, if you define a per-particle parameter called "charge", then the variable "charge2" is the charge of particle p2.
+ * select(x,y,z) = z if x = 0, y otherwise.  The names of per-particle parameters have the suffix "1", "2", etc. appended to them to indicate the values for
+ * the multiple interacting particles. For example, if you define a per-particle parameter called "charge", then the variable "charge2" is the charge of particle p2.
  * As seen above, the expression may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
  *
  * In addition, you can call addTabulatedFunction() to define a new function based on tabulated values.  You specify the function by

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -120,10 +120,11 @@ namespace OpenMM {
  * frequently, the long range correction can be very slow.  For this reason, it is disabled by default.
  * 
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
- * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta.  All trigonometric functions
- * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.  The names of per-particle parameters
- * have the suffix "1" or "2" appended to them to indicate the values for the two interacting particles.  As seen in the above example,
- * the expression may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
+ * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
+ * select(x,y,z) = z if x = 0, y otherwise.  The names of per-particle parameters have the suffix "1" or "2" appended to them to indicate the values for the
+ * two interacting particles.  As seen in the above example, the expression may also involve intermediate quantities that are defined following the main expression,
+ * using ";" as a separator.
  *
  * In addition, you can call addTabulatedFunction() to define a new function based on tabulated values.  You specify the function by
  * creating a TabulatedFunction object.  That function can then appear in the expression.

openmmapi/include/openmm/CustomTorsionForce.h

@@ -65,8 +65,9 @@ namespace OpenMM {
  * </pre></tt>
  *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
- * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta.  All trigonometric functions
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
  * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
+ * select(x,y,z) = z if x = 0, y otherwise.
  */
 
 class OPENMM_EXPORT CustomTorsionForce : public Force {

platforms/cuda/src/CudaExpressionUtilities.cpp

@@ -463,6 +463,9 @@ void CudaExpressionUtilities::processExpression(stringstream& out, const Express
         case Operation::CEIL:
             out << "ceil(" << getTempName(node.getChildren()[0], temps) << ")";
             break;
+        case Operation::SELECT:
+            out << "(" << getTempName(node.getChildren()[0], temps) << " != 0 ? " << getTempName(node.getChildren()[1], temps) << " : " << getTempName(node.getChildren()[2], temps) << ")";
+            break;
         default:
             throw OpenMMException("Internal error: Unknown operation in user-defined expression: "+node.getOperation().getName());
     }

platforms/cuda/tests/TestCudaCustomIntegrator.cpp

@@ -334,7 +334,7 @@ void testMonteCarlo() {
     integrator.addComputePerDof("oldx", "x");
     integrator.addComputePerDof("x", "x+dt*gaussian");
     integrator.addComputeGlobal("accept", "step(exp((oldE-energy)/kT)-uniform)");
-    integrator.addComputePerDof("x", "accept*x + (1-accept)*oldx");
+    integrator.addComputePerDof("x", "select(accept, x, oldx)");
     HarmonicBondForce* forceField = new HarmonicBondForce();
     forceField->addBond(0, 1, 2.0, 10.0);
     system.addForce(forceField);

platforms/opencl/src/OpenCLExpressionUtilities.cpp

@@ -455,6 +455,9 @@ void OpenCLExpressionUtilities::processExpression(stringstream& out, const Expre
         case Operation::CEIL:
             out << "ceil(" << getTempName(node.getChildren()[0], temps) << ")";
             break;
+        case Operation::SELECT:
+            out << "(" << getTempName(node.getChildren()[0], temps) << " != 0 ? " << getTempName(node.getChildren()[1], temps) << " : " << getTempName(node.getChildren()[2], temps) << ")";
+            break;
         default:
             throw OpenMMException("Internal error: Unknown operation in user-defined expression: "+node.getOperation().getName());
     }

platforms/opencl/tests/TestOpenCLCustomIntegrator.cpp

@@ -334,7 +334,7 @@ void testMonteCarlo() {
     integrator.addComputePerDof("oldx", "x");
     integrator.addComputePerDof("x", "x+dt*gaussian");
     integrator.addComputeGlobal("accept", "step(exp((oldE-energy)/kT)-uniform)");
-    integrator.addComputePerDof("x", "accept*x + (1-accept)*oldx");
+    integrator.addComputePerDof("x", "select(accept, x, oldx)");
     HarmonicBondForce* forceField = new HarmonicBondForce();
     forceField->addBond(0, 1, 2.0, 10.0);
     system.addForce(forceField);

platforms/reference/tests/TestReferenceCustomIntegrator.cpp

@@ -334,7 +334,7 @@ void testMonteCarlo() {
     integrator.addComputePerDof("oldx", "x");
     integrator.addComputePerDof("x", "x+dt*gaussian");
     integrator.addComputeGlobal("accept", "step(exp((oldE-energy)/kT)-uniform)");
-    integrator.addComputePerDof("x", "accept*x + (1-accept)*oldx");
+    integrator.addComputePerDof("x", "select(accept, x, oldx)");
     HarmonicBondForce* forceField = new HarmonicBondForce();
     forceField->addBond(0, 1, 2.0, 10.0);
     system.addForce(forceField);