Update Lepton to version 2019-06-02 02/5202/2
authorJérôme Hénin <heninj@gmail.com>
Thu, 20 Jun 2019 11:10:01 +0000 (13:10 +0200)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Tue, 25 Jun 2019 16:21:09 +0000 (11:21 -0500)
Fixes and optimizations, added atan2() function.

Change-Id: Iec8d73f974fb57fd48e58206dcac9a6d1afa0094

lepton/include/lepton/CompiledExpression.h
lepton/include/lepton/CustomFunction.h
lepton/include/lepton/ExpressionProgram.h
lepton/include/lepton/Operation.h
lepton/src/CompiledExpression.cpp
lepton/src/ExpressionProgram.cpp
lepton/src/MSVC_erfc.h
lepton/src/Operation.cpp
lepton/src/ParsedExpression.cpp
lepton/src/Parser.cpp

index 67442e0..c7e393e 100644 (file)
@@ -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) 2013-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -52,9 +52,9 @@ class ParsedExpression;
  * A CompiledExpression is a highly optimized representation of an expression for cases when you want to evaluate
  * it many times as quickly as possible.  You should treat it as an opaque object; none of the internal representation
  * is visible.
- * 
+ *
  * A CompiledExpression is created by calling createCompiledExpression() on a ParsedExpression.
- * 
+ *
  * WARNING: CompiledExpression is NOT thread safe.  You should never access a CompiledExpression from two threads at
  * the same time.
  */
@@ -99,10 +99,11 @@ private:
     mutable std::vector<double> workspace;
     mutable std::vector<double> argValues;
     std::map<std::string, double> dummyVariables;
-    void* jitCode;
+    double (*jitCode)();
 #ifdef LEPTON_USE_JIT
     void generateJitCode();
-    void generateSingleArgCall(asmjit::X86Compiler& c, asmjit::X86XmmVar& dest, asmjit::X86XmmVar& arg, double (*function)(double));
+    void generateSingleArgCall(asmjit::X86Compiler& c, asmjit::X86Xmm& dest, asmjit::X86Xmm& arg, double (*function)(double));
+    void generateTwoArgCall(asmjit::X86Compiler& c, asmjit::X86Xmm& dest, asmjit::X86Xmm& arg1, asmjit::X86Xmm& arg2, double (*function)(double, double));
     std::vector<double> constants;
     asmjit::JitRuntime runtime;
 #endif
index 5c55861..fbb0ddd 100644 (file)
@@ -72,6 +72,38 @@ public:
     virtual CustomFunction* clone() const = 0;
 };
 
+/**
+ * This class is an implementation of CustomFunction that does no computation.  It just returns
+ * 0 for the value and derivatives.  This is useful when using the parser to analyze expressions
+ * rather than to evaluate them.  You can just create PlaceholderFunctions to represent any custom
+ * functions that may appear in expressions.
+ */
+
+class LEPTON_EXPORT PlaceholderFunction : public CustomFunction {
+public:
+    /**
+     * Create a Placeholder function.
+     *
+     * @param numArgs    the number of arguments the function expects
+     */
+    PlaceholderFunction(int numArgs) : numArgs(numArgs) {
+    }
+    int getNumArguments() const {
+        return numArgs;
+    }
+    double evaluate(const double* arguments) const {
+        return 0.0;
+    }
+    double evaluateDerivative(const double* arguments, const int* derivOrder) const {
+        return 0.0;
+    }
+    CustomFunction* clone() const {
+        return new PlaceholderFunction(numArgs);
+    };
+private:
+    int numArgs;
+};
+
 } // namespace Lepton
 
 #endif /*LEPTON_CUSTOM_FUNCTION_H_*/
index 94d37f4..a49a909 100644 (file)
@@ -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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2009-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -65,6 +65,14 @@ public:
      * Get an Operation in this program.
      */
     const Operation& getOperation(int index) const;
+    /**
+     * Change an Operation in this program.
+     *
+     * The Operation must have been allocated on the heap with the "new" operator.
+     * The ExpressionProgram assumes ownership of it and will delete it when it
+     * is no longer needed.
+     */
+    void setOperation(int index, Operation* operation);
     /**
      * Get the size of the stack needed to execute this program.  This is the largest number of elements present
      * on the stack at any point during evaluation.
index f7a8b78..1ddde0b 100644 (file)
@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -63,7 +63,7 @@ public:
      * can be used when processing or analyzing parsed expressions.
      */
     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,
+             SIN, COS, SEC, CSC, TAN, COT, ASIN, ACOS, ATAN, ATAN2, SINH, COSH, TANH, ERF, ERFC, STEP, DELTA, SQUARE, CUBE, RECIPROCAL,
              ADD_CONSTANT, MULTIPLY_CONSTANT, POWER_CONSTANT, MIN, MAX, ABS, FLOOR, CEIL, SELECT};
     /**
      * Get the name of this Operation.
@@ -137,6 +137,7 @@ public:
     class Asin;
     class Acos;
     class Atan;
+    class Atan2;
     class Sinh;
     class Cosh;
     class Tanh;
@@ -226,6 +227,11 @@ class LEPTON_EXPORT Operation::Custom : public Operation {
 public:
     Custom(const std::string& name, CustomFunction* function) : name(name), function(function), isDerivative(false), derivOrder(function->getNumArguments(), 0) {
     }
+    Custom(const std::string& name, CustomFunction* function, const std::vector<int>& derivOrder) : name(name), function(function), isDerivative(false), derivOrder(derivOrder) {
+        for (int order : derivOrder)
+            if (order != 0)
+                isDerivative = true;
+    }
     Custom(const Custom& base, int derivIndex) : name(base.name), function(base.function->clone()), isDerivative(true), derivOrder(base.derivOrder) {
         derivOrder[derivIndex]++;
     }
@@ -684,6 +690,28 @@ public:
     ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const;
 };
 
+class LEPTON_EXPORT Operation::Atan2 : public Operation {
+public:
+    Atan2() {
+    }
+    std::string getName() const {
+        return "atan2";
+    }
+    Id getId() const {
+        return ATAN2;
+    }
+    int getNumArguments() const {
+        return 2;
+    }
+    Operation* clone() const {
+        return new Atan2();
+    }
+    double evaluate(double* args, const std::map<std::string, double>& variables) const {
+        return std::atan2(args[0], args[1]);
+    }
+    ExpressionTreeNode differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const;
+};
+
 class LEPTON_EXPORT Operation::Sinh : public Operation {
 public:
     Sinh() {
@@ -989,7 +1017,7 @@ public:
     double evaluate(double* args, const std::map<std::string, double>& variables) const {
         if (isIntPower) {
             // Integer powers can be computed much more quickly by repeated multiplication.
-            
+
             int exponent = intValue;
             double base = args[0];
             if (exponent < 0) {
index 84f827d..1ad348b 100644 (file)
-/* -------------------------------------------------------------------------- *\r
- *                                   Lepton                                   *\r
- * -------------------------------------------------------------------------- *\r
- * This is part of the Lepton expression parser originating from              *\r
- * Simbios, the NIH National Center for Physics-Based Simulation of           *\r
- * Biological Structures at Stanford, funded under the NIH Roadmap for        *\r
- * Medical Research, grant U54 GM072970. See https://simtk.org.               *\r
- *                                                                            *\r
- * Portions copyright (c) 2013-2016 Stanford University and the Authors.      *\r
- * Authors: Peter Eastman                                                     *\r
- * Contributors:                                                              *\r
- *                                                                            *\r
- * Permission is hereby granted, free of charge, to any person obtaining a    *\r
- * copy of this software and associated documentation files (the "Software"), *\r
- * to deal in the Software without restriction, including without limitation  *\r
- * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *\r
- * and/or sell copies of the Software, and to permit persons to whom the      *\r
- * Software is furnished to do so, subject to the following conditions:       *\r
- *                                                                            *\r
- * The above copyright notice and this permission notice shall be included in *\r
- * all copies or substantial portions of the Software.                        *\r
- *                                                                            *\r
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *\r
- * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *\r
- * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *\r
- * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *\r
- * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *\r
- * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *\r
- * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *\r
- * -------------------------------------------------------------------------- */\r
-\r
-#include "lepton/CompiledExpression.h"\r
-#include "lepton/Operation.h"\r
-#include "lepton/ParsedExpression.h"\r
-#include <utility>\r
-\r
-using namespace Lepton;\r
-using namespace std;\r
-#ifdef LEPTON_USE_JIT\r
-    using namespace asmjit;\r
-#endif\r
-\r
-CompiledExpression::CompiledExpression() : jitCode(NULL) {\r
-}\r
-\r
-CompiledExpression::CompiledExpression(const ParsedExpression& expression) : jitCode(NULL) {\r
-    ParsedExpression expr = expression.optimize(); // Just in case it wasn't already optimized.\r
-    vector<pair<ExpressionTreeNode, int> > temps;\r
-    compileExpression(expr.getRootNode(), temps);\r
-    int maxArguments = 1;\r
-    for (int i = 0; i < (int) operation.size(); i++)\r
-        if (operation[i]->getNumArguments() > maxArguments)\r
-            maxArguments = operation[i]->getNumArguments();\r
-    argValues.resize(maxArguments);\r
-#ifdef LEPTON_USE_JIT\r
-    generateJitCode();\r
-#endif\r
-}\r
-\r
-CompiledExpression::~CompiledExpression() {\r
-    for (int i = 0; i < (int) operation.size(); i++)\r
-        if (operation[i] != NULL)\r
-            delete operation[i];\r
-}\r
-\r
-CompiledExpression::CompiledExpression(const CompiledExpression& expression) : jitCode(NULL) {\r
-    *this = expression;\r
-}\r
-\r
-CompiledExpression& CompiledExpression::operator=(const CompiledExpression& expression) {\r
-    arguments = expression.arguments;\r
-    target = expression.target;\r
-    variableIndices = expression.variableIndices;\r
-    variableNames = expression.variableNames;\r
-    workspace.resize(expression.workspace.size());\r
-    argValues.resize(expression.argValues.size());\r
-    operation.resize(expression.operation.size());\r
-    for (int i = 0; i < (int) operation.size(); i++)\r
-        operation[i] = expression.operation[i]->clone();\r
-    setVariableLocations(variablePointers);\r
-    return *this;\r
-}\r
-\r
-void CompiledExpression::compileExpression(const ExpressionTreeNode& node, vector<pair<ExpressionTreeNode, int> >& temps) {\r
-    if (findTempIndex(node, temps) != -1)\r
-        return; // We have already processed a node identical to this one.\r
-    \r
-    // Process the child nodes.\r
-    \r
-    vector<int> args;\r
-    for (int i = 0; i < node.getChildren().size(); i++) {\r
-        compileExpression(node.getChildren()[i], temps);\r
-        args.push_back(findTempIndex(node.getChildren()[i], temps));\r
-    }\r
-    \r
-    // Process this node.\r
-    \r
-    if (node.getOperation().getId() == Operation::VARIABLE) {\r
-        variableIndices[node.getOperation().getName()] = (int) workspace.size();\r
-        variableNames.insert(node.getOperation().getName());\r
-    }\r
-    else {\r
-        int stepIndex = (int) arguments.size();\r
-        arguments.push_back(vector<int>());\r
-        target.push_back((int) workspace.size());\r
-        operation.push_back(node.getOperation().clone());\r
-        if (args.size() == 0)\r
-            arguments[stepIndex].push_back(0); // The value won't actually be used.  We just need something there.\r
-        else {\r
-            // If the arguments are sequential, we can just pass a pointer to the first one.\r
-            \r
-            bool sequential = true;\r
-            for (int i = 1; i < args.size(); i++)\r
-                if (args[i] != args[i-1]+1)\r
-                    sequential = false;\r
-            if (sequential)\r
-                arguments[stepIndex].push_back(args[0]);\r
-            else\r
-                arguments[stepIndex] = args;\r
-        }\r
-    }\r
-    temps.push_back(make_pair(node, (int) workspace.size()));\r
-    workspace.push_back(0.0);\r
-}\r
-\r
-int CompiledExpression::findTempIndex(const ExpressionTreeNode& node, vector<pair<ExpressionTreeNode, int> >& temps) {\r
-    for (int i = 0; i < (int) temps.size(); i++)\r
-        if (temps[i].first == node)\r
-            return i;\r
-    return -1;\r
-}\r
-\r
-const set<string>& CompiledExpression::getVariables() const {\r
-    return variableNames;\r
-}\r
-\r
-double& CompiledExpression::getVariableReference(const string& name) {\r
-    map<string, double*>::iterator pointer = variablePointers.find(name);\r
-    if (pointer != variablePointers.end())\r
-        return *pointer->second;\r
-    map<string, int>::iterator index = variableIndices.find(name);\r
-    if (index == variableIndices.end())\r
-        throw Exception("getVariableReference: Unknown variable '"+name+"'");\r
-    return workspace[index->second];\r
-}\r
-\r
-void CompiledExpression::setVariableLocations(map<string, double*>& variableLocations) {\r
-    variablePointers = variableLocations;\r
-#ifdef LEPTON_USE_JIT\r
-    // Rebuild the JIT code.\r
-    \r
-    if (workspace.size() > 0)\r
-        generateJitCode();\r
-#else\r
-    // Make a list of all variables we will need to copy before evaluating the expression.\r
-    \r
-    variablesToCopy.clear();\r
-    for (map<string, int>::const_iterator iter = variableIndices.begin(); iter != variableIndices.end(); ++iter) {\r
-        map<string, double*>::iterator pointer = variablePointers.find(iter->first);\r
-        if (pointer != variablePointers.end())\r
-            variablesToCopy.push_back(make_pair(&workspace[iter->second], pointer->second));\r
-    }\r
-#endif\r
-}\r
-\r
-double CompiledExpression::evaluate() const {\r
-#ifdef LEPTON_USE_JIT\r
-    return ((double (*)()) jitCode)();\r
-#else\r
-    for (int i = 0; i < variablesToCopy.size(); i++)\r
-        *variablesToCopy[i].first = *variablesToCopy[i].second;\r
-\r
-    // Loop over the operations and evaluate each one.\r
-    \r
-    for (int step = 0; step < operation.size(); step++) {\r
-        const vector<int>& args = arguments[step];\r
-        if (args.size() == 1)\r
-            workspace[target[step]] = operation[step]->evaluate(&workspace[args[0]], dummyVariables);\r
-        else {\r
-            for (int i = 0; i < args.size(); i++)\r
-                argValues[i] = workspace[args[i]];\r
-            workspace[target[step]] = operation[step]->evaluate(&argValues[0], dummyVariables);\r
-        }\r
-    }\r
-    return workspace[workspace.size()-1];\r
-#endif\r
-}\r
-\r
-#ifdef LEPTON_USE_JIT\r
-static double evaluateOperation(Operation* op, double* args) {\r
-    map<string, double>* dummyVariables = NULL;\r
-    return op->evaluate(args, *dummyVariables);\r
-}\r
-\r
-void CompiledExpression::generateJitCode() {\r
-    X86Compiler c(&runtime);\r
-    c.addFunc(kFuncConvHost, FuncBuilder0<double>());\r
-    vector<X86XmmVar> workspaceVar(workspace.size());\r
-    for (int i = 0; i < (int) workspaceVar.size(); i++)\r
-        workspaceVar[i] = c.newXmmVar(kX86VarTypeXmmSd);\r
-    X86GpVar argsPointer(c);\r
-    c.mov(argsPointer, imm_ptr(&argValues[0]));\r
-    \r
-    // Load the arguments into variables.\r
-    \r
-    for (set<string>::const_iterator iter = variableNames.begin(); iter != variableNames.end(); ++iter) {\r
-        map<string, int>::iterator index = variableIndices.find(*iter);\r
-        X86GpVar variablePointer(c);\r
-        c.mov(variablePointer, imm_ptr(&getVariableReference(index->first)));\r
-        c.movsd(workspaceVar[index->second], x86::ptr(variablePointer, 0, 0));\r
-    }\r
-\r
-    // Make a list of all constants that will be needed for evaluation.\r
-    \r
-    vector<int> operationConstantIndex(operation.size(), -1);\r
-    for (int step = 0; step < (int) operation.size(); step++) {\r
-        // Find the constant value (if any) used by this operation.\r
-        \r
-        Operation& op = *operation[step];\r
-        double value;\r
-        if (op.getId() == Operation::CONSTANT)\r
-            value = dynamic_cast<Operation::Constant&>(op).getValue();\r
-        else if (op.getId() == Operation::ADD_CONSTANT)\r
-            value = dynamic_cast<Operation::AddConstant&>(op).getValue();\r
-        else if (op.getId() == Operation::MULTIPLY_CONSTANT)\r
-            value = dynamic_cast<Operation::MultiplyConstant&>(op).getValue();\r
-        else if (op.getId() == Operation::RECIPROCAL)\r
-            value = 1.0;\r
-        else if (op.getId() == Operation::STEP)\r
-            value = 1.0;\r
-        else if (op.getId() == Operation::DELTA)\r
-            value = 1.0;\r
-        else\r
-            continue;\r
-        \r
-        // See if we already have a variable for this constant.\r
-        \r
-        for (int i = 0; i < (int) constants.size(); i++)\r
-            if (value == constants[i]) {\r
-                operationConstantIndex[step] = i;\r
-                break;\r
-            }\r
-        if (operationConstantIndex[step] == -1) {\r
-            operationConstantIndex[step] = constants.size();\r
-            constants.push_back(value);\r
-        }\r
-    }\r
-    \r
-    // Load constants into variables.\r
-    \r
-    vector<X86XmmVar> constantVar(constants.size());\r
-    if (constants.size() > 0) {\r
-        X86GpVar constantsPointer(c);\r
-        c.mov(constantsPointer, imm_ptr(&constants[0]));\r
-        for (int i = 0; i < (int) constants.size(); i++) {\r
-            constantVar[i] = c.newXmmVar(kX86VarTypeXmmSd);\r
-            c.movsd(constantVar[i], x86::ptr(constantsPointer, 8*i, 0));\r
-        }\r
-    }\r
-    \r
-    // Evaluate the operations.\r
-    \r
-    for (int step = 0; step < (int) operation.size(); step++) {\r
-        Operation& op = *operation[step];\r
-        vector<int> args = arguments[step];\r
-        if (args.size() == 1) {\r
-            // One or more sequential arguments.  Fill out the list.\r
-            \r
-            for (int i = 1; i < op.getNumArguments(); i++)\r
-                args.push_back(args[0]+i);\r
-        }\r
-        \r
-        // Generate instructions to execute this operation.\r
-        \r
-        switch (op.getId()) {\r
-            case Operation::CONSTANT:\r
-                c.movsd(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);\r
-                break;\r
-            case Operation::ADD:\r
-                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                c.addsd(workspaceVar[target[step]], workspaceVar[args[1]]);\r
-                break;\r
-            case Operation::SUBTRACT:\r
-                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                c.subsd(workspaceVar[target[step]], workspaceVar[args[1]]);\r
-                break;\r
-            case Operation::MULTIPLY:\r
-                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                c.mulsd(workspaceVar[target[step]], workspaceVar[args[1]]);\r
-                break;\r
-            case Operation::DIVIDE:\r
-                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                c.divsd(workspaceVar[target[step]], workspaceVar[args[1]]);\r
-                break;\r
-            case Operation::NEGATE:\r
-                c.xorps(workspaceVar[target[step]], workspaceVar[target[step]]);\r
-                c.subsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                break;\r
-            case Operation::SQRT:\r
-                c.sqrtsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                break;\r
-            case Operation::EXP:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], exp);\r
-                break;\r
-            case Operation::LOG:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], log);\r
-                break;\r
-            case Operation::SIN:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], sin);\r
-                break;\r
-            case Operation::COS:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], cos);\r
-                break;\r
-            case Operation::TAN:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], tan);\r
-                break;\r
-            case Operation::ASIN:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], asin);\r
-                break;\r
-            case Operation::ACOS:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], acos);\r
-                break;\r
-            case Operation::ATAN:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], atan);\r
-                break;\r
-            case Operation::SINH:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], sinh);\r
-                break;\r
-            case Operation::COSH:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], cosh);\r
-                break;\r
-            case Operation::TANH:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], tanh);\r
-                break;\r
-            case Operation::STEP:\r
-                c.xorps(workspaceVar[target[step]], workspaceVar[target[step]]);\r
-                c.cmpsd(workspaceVar[target[step]], workspaceVar[args[0]], imm(18)); // Comparison mode is _CMP_LE_OQ = 18\r
-                c.andps(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);\r
-                break;\r
-            case Operation::DELTA:\r
-                c.xorps(workspaceVar[target[step]], workspaceVar[target[step]]);\r
-                c.cmpsd(workspaceVar[target[step]], workspaceVar[args[0]], imm(16)); // Comparison mode is _CMP_EQ_OS = 16\r
-                c.andps(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);\r
-                break;\r
-            case Operation::SQUARE:\r
-                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                c.mulsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                break;\r
-            case Operation::CUBE:\r
-                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                c.mulsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                c.mulsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                break;\r
-            case Operation::RECIPROCAL:\r
-                c.movsd(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);\r
-                c.divsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                break;\r
-            case Operation::ADD_CONSTANT:\r
-                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                c.addsd(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);\r
-                break;\r
-            case Operation::MULTIPLY_CONSTANT:\r
-                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);\r
-                c.mulsd(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);\r
-                break;\r
-            case Operation::ABS:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], fabs);\r
-                break;\r
-            case Operation::FLOOR:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], floor);\r
-                break;\r
-            case Operation::CEIL:\r
-                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], ceil);\r
-                break;\r
-            default:\r
-                // Just invoke evaluateOperation().\r
-                \r
-                for (int i = 0; i < (int) args.size(); i++)\r
-                    c.movsd(x86::ptr(argsPointer, 8*i, 0), workspaceVar[args[i]]);\r
-                X86GpVar fn(c, kVarTypeIntPtr);\r
-                c.mov(fn, imm_ptr((void*) evaluateOperation));\r
-                X86CallNode* call = c.call(fn, kFuncConvHost, FuncBuilder2<double, Operation*, double*>());\r
-                call->setArg(0, imm_ptr(&op));\r
-                call->setArg(1, imm_ptr(&argValues[0]));\r
-                call->setRet(0, workspaceVar[target[step]]);\r
-        }\r
-    }\r
-    c.ret(workspaceVar[workspace.size()-1]);\r
-    c.endFunc();\r
-    jitCode = c.make();\r
-}\r
-\r
-void CompiledExpression::generateSingleArgCall(X86Compiler& c, X86XmmVar& dest, X86XmmVar& arg, double (*function)(double)) {\r
-    X86GpVar fn(c, kVarTypeIntPtr);\r
-    c.mov(fn, imm_ptr((void*) function));\r
-    X86CallNode* call = c.call(fn, kFuncConvHost, FuncBuilder1<double, double>());\r
-    call->setArg(0, arg);\r
-    call->setRet(0, dest);\r
-}\r
-#endif\r
+/* -------------------------------------------------------------------------- *
+ *                                   Lepton                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the Lepton expression parser originating from              *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2013-2019 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "lepton/CompiledExpression.h"
+#include "lepton/Operation.h"
+#include "lepton/ParsedExpression.h"
+#include <utility>
+
+using namespace Lepton;
+using namespace std;
+#ifdef LEPTON_USE_JIT
+    using namespace asmjit;
+#endif
+
+CompiledExpression::CompiledExpression() : jitCode(NULL) {
+}
+
+CompiledExpression::CompiledExpression(const ParsedExpression& expression) : jitCode(NULL) {
+    ParsedExpression expr = expression.optimize(); // Just in case it wasn't already optimized.
+    vector<pair<ExpressionTreeNode, int> > temps;
+    compileExpression(expr.getRootNode(), temps);
+    int maxArguments = 1;
+    for (int i = 0; i < (int) operation.size(); i++)
+        if (operation[i]->getNumArguments() > maxArguments)
+            maxArguments = operation[i]->getNumArguments();
+    argValues.resize(maxArguments);
+#ifdef LEPTON_USE_JIT
+    generateJitCode();
+#endif
+}
+
+CompiledExpression::~CompiledExpression() {
+    for (int i = 0; i < (int) operation.size(); i++)
+        if (operation[i] != NULL)
+            delete operation[i];
+}
+
+CompiledExpression::CompiledExpression(const CompiledExpression& expression) : jitCode(NULL) {
+    *this = expression;
+}
+
+CompiledExpression& CompiledExpression::operator=(const CompiledExpression& expression) {
+    arguments = expression.arguments;
+    target = expression.target;
+    variableIndices = expression.variableIndices;
+    variableNames = expression.variableNames;
+    workspace.resize(expression.workspace.size());
+    argValues.resize(expression.argValues.size());
+    operation.resize(expression.operation.size());
+    for (int i = 0; i < (int) operation.size(); i++)
+        operation[i] = expression.operation[i]->clone();
+    setVariableLocations(variablePointers);
+    return *this;
+}
+
+void CompiledExpression::compileExpression(const ExpressionTreeNode& node, vector<pair<ExpressionTreeNode, int> >& temps) {
+    if (findTempIndex(node, temps) != -1)
+        return; // We have already processed a node identical to this one.
+
+    // Process the child nodes.
+
+    vector<int> args;
+    for (int i = 0; i < node.getChildren().size(); i++) {
+        compileExpression(node.getChildren()[i], temps);
+        args.push_back(findTempIndex(node.getChildren()[i], temps));
+    }
+
+    // Process this node.
+
+    if (node.getOperation().getId() == Operation::VARIABLE) {
+        variableIndices[node.getOperation().getName()] = (int) workspace.size();
+        variableNames.insert(node.getOperation().getName());
+    }
+    else {
+        int stepIndex = (int) arguments.size();
+        arguments.push_back(vector<int>());
+        target.push_back((int) workspace.size());
+        operation.push_back(node.getOperation().clone());
+        if (args.size() == 0)
+            arguments[stepIndex].push_back(0); // The value won't actually be used.  We just need something there.
+        else {
+            // If the arguments are sequential, we can just pass a pointer to the first one.
+
+            bool sequential = true;
+            for (int i = 1; i < args.size(); i++)
+                if (args[i] != args[i-1]+1)
+                    sequential = false;
+            if (sequential)
+                arguments[stepIndex].push_back(args[0]);
+            else
+                arguments[stepIndex] = args;
+        }
+    }
+    temps.push_back(make_pair(node, (int) workspace.size()));
+    workspace.push_back(0.0);
+}
+
+int CompiledExpression::findTempIndex(const ExpressionTreeNode& node, vector<pair<ExpressionTreeNode, int> >& temps) {
+    for (int i = 0; i < (int) temps.size(); i++)
+        if (temps[i].first == node)
+            return i;
+    return -1;
+}
+
+const set<string>& CompiledExpression::getVariables() const {
+    return variableNames;
+}
+
+double& CompiledExpression::getVariableReference(const string& name) {
+    map<string, double*>::iterator pointer = variablePointers.find(name);
+    if (pointer != variablePointers.end())
+        return *pointer->second;
+    map<string, int>::iterator index = variableIndices.find(name);
+    if (index == variableIndices.end())
+        throw Exception("getVariableReference: Unknown variable '"+name+"'");
+    return workspace[index->second];
+}
+
+void CompiledExpression::setVariableLocations(map<string, double*>& variableLocations) {
+    variablePointers = variableLocations;
+#ifdef LEPTON_USE_JIT
+    // Rebuild the JIT code.
+
+    if (workspace.size() > 0)
+        generateJitCode();
+#else
+    // Make a list of all variables we will need to copy before evaluating the expression.
+
+    variablesToCopy.clear();
+    for (map<string, int>::const_iterator iter = variableIndices.begin(); iter != variableIndices.end(); ++iter) {
+        map<string, double*>::iterator pointer = variablePointers.find(iter->first);
+        if (pointer != variablePointers.end())
+            variablesToCopy.push_back(make_pair(&workspace[iter->second], pointer->second));
+    }
+#endif
+}
+
+double CompiledExpression::evaluate() const {
+#ifdef LEPTON_USE_JIT
+    return jitCode();
+#else
+    for (int i = 0; i < variablesToCopy.size(); i++)
+        *variablesToCopy[i].first = *variablesToCopy[i].second;
+
+    // Loop over the operations and evaluate each one.
+
+    for (int step = 0; step < operation.size(); step++) {
+        const vector<int>& args = arguments[step];
+        if (args.size() == 1)
+            workspace[target[step]] = operation[step]->evaluate(&workspace[args[0]], dummyVariables);
+        else {
+            for (int i = 0; i < args.size(); i++)
+                argValues[i] = workspace[args[i]];
+            workspace[target[step]] = operation[step]->evaluate(&argValues[0], dummyVariables);
+        }
+    }
+    return workspace[workspace.size()-1];
+#endif
+}
+
+#ifdef LEPTON_USE_JIT
+static double evaluateOperation(Operation* op, double* args) {
+    static map<string, double> dummyVariables;
+    return op->evaluate(args, dummyVariables);
+}
+
+void CompiledExpression::generateJitCode() {
+    CodeHolder code;
+    code.init(runtime.getCodeInfo());
+    X86Compiler c(&code);
+    c.addFunc(FuncSignature0<double>());
+    vector<X86Xmm> workspaceVar(workspace.size());
+    for (int i = 0; i < (int) workspaceVar.size(); i++)
+        workspaceVar[i] = c.newXmmSd();
+    X86Gp argsPointer = c.newIntPtr();
+    c.mov(argsPointer, imm_ptr(&argValues[0]));
+
+    // Load the arguments into variables.
+
+    for (set<string>::const_iterator iter = variableNames.begin(); iter != variableNames.end(); ++iter) {
+        map<string, int>::iterator index = variableIndices.find(*iter);
+        X86Gp variablePointer = c.newIntPtr();
+        c.mov(variablePointer, imm_ptr(&getVariableReference(index->first)));
+        c.movsd(workspaceVar[index->second], x86::ptr(variablePointer, 0, 0));
+    }
+
+    // Make a list of all constants that will be needed for evaluation.
+
+    vector<int> operationConstantIndex(operation.size(), -1);
+    for (int step = 0; step < (int) operation.size(); step++) {
+        // Find the constant value (if any) used by this operation.
+
+        Operation& op = *operation[step];
+        double value;
+        if (op.getId() == Operation::CONSTANT)
+            value = dynamic_cast<Operation::Constant&>(op).getValue();
+        else if (op.getId() == Operation::ADD_CONSTANT)
+            value = dynamic_cast<Operation::AddConstant&>(op).getValue();
+        else if (op.getId() == Operation::MULTIPLY_CONSTANT)
+            value = dynamic_cast<Operation::MultiplyConstant&>(op).getValue();
+        else if (op.getId() == Operation::RECIPROCAL)
+            value = 1.0;
+        else if (op.getId() == Operation::STEP)
+            value = 1.0;
+        else if (op.getId() == Operation::DELTA)
+            value = 1.0;
+        else
+            continue;
+
+        // See if we already have a variable for this constant.
+
+        for (int i = 0; i < (int) constants.size(); i++)
+            if (value == constants[i]) {
+                operationConstantIndex[step] = i;
+                break;
+            }
+        if (operationConstantIndex[step] == -1) {
+            operationConstantIndex[step] = constants.size();
+            constants.push_back(value);
+        }
+    }
+
+    // Load constants into variables.
+
+    vector<X86Xmm> constantVar(constants.size());
+    if (constants.size() > 0) {
+        X86Gp constantsPointer = c.newIntPtr();
+        c.mov(constantsPointer, imm_ptr(&constants[0]));
+        for (int i = 0; i < (int) constants.size(); i++) {
+            constantVar[i] = c.newXmmSd();
+            c.movsd(constantVar[i], x86::ptr(constantsPointer, 8*i, 0));
+        }
+    }
+
+    // Evaluate the operations.
+
+    for (int step = 0; step < (int) operation.size(); step++) {
+        Operation& op = *operation[step];
+        vector<int> args = arguments[step];
+        if (args.size() == 1) {
+            // One or more sequential arguments.  Fill out the list.
+
+            for (int i = 1; i < op.getNumArguments(); i++)
+                args.push_back(args[0]+i);
+        }
+
+        // Generate instructions to execute this operation.
+
+        switch (op.getId()) {
+            case Operation::CONSTANT:
+                c.movsd(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);
+                break;
+            case Operation::ADD:
+                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                c.addsd(workspaceVar[target[step]], workspaceVar[args[1]]);
+                break;
+            case Operation::SUBTRACT:
+                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                c.subsd(workspaceVar[target[step]], workspaceVar[args[1]]);
+                break;
+            case Operation::MULTIPLY:
+                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                c.mulsd(workspaceVar[target[step]], workspaceVar[args[1]]);
+                break;
+            case Operation::DIVIDE:
+                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                c.divsd(workspaceVar[target[step]], workspaceVar[args[1]]);
+                break;
+            case Operation::POWER:
+                generateTwoArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], workspaceVar[args[1]], pow);
+                break;
+            case Operation::NEGATE:
+                c.xorps(workspaceVar[target[step]], workspaceVar[target[step]]);
+                c.subsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                break;
+            case Operation::SQRT:
+                c.sqrtsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                break;
+            case Operation::EXP:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], exp);
+                break;
+            case Operation::LOG:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], log);
+                break;
+            case Operation::SIN:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], sin);
+                break;
+            case Operation::COS:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], cos);
+                break;
+            case Operation::TAN:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], tan);
+                break;
+            case Operation::ASIN:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], asin);
+                break;
+            case Operation::ACOS:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], acos);
+                break;
+            case Operation::ATAN:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], atan);
+                break;
+            case Operation::ATAN2:
+                generateTwoArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], workspaceVar[args[1]], atan2);
+                break;
+            case Operation::SINH:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], sinh);
+                break;
+            case Operation::COSH:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], cosh);
+                break;
+            case Operation::TANH:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], tanh);
+                break;
+            case Operation::STEP:
+                c.xorps(workspaceVar[target[step]], workspaceVar[target[step]]);
+                c.cmpsd(workspaceVar[target[step]], workspaceVar[args[0]], imm(18)); // Comparison mode is _CMP_LE_OQ = 18
+                c.andps(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);
+                break;
+            case Operation::DELTA:
+                c.xorps(workspaceVar[target[step]], workspaceVar[target[step]]);
+                c.cmpsd(workspaceVar[target[step]], workspaceVar[args[0]], imm(16)); // Comparison mode is _CMP_EQ_OS = 16
+                c.andps(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);
+                break;
+            case Operation::SQUARE:
+                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                c.mulsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                break;
+            case Operation::CUBE:
+                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                c.mulsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                c.mulsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                break;
+            case Operation::RECIPROCAL:
+                c.movsd(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);
+                c.divsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                break;
+            case Operation::ADD_CONSTANT:
+                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                c.addsd(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);
+                break;
+            case Operation::MULTIPLY_CONSTANT:
+                c.movsd(workspaceVar[target[step]], workspaceVar[args[0]]);
+                c.mulsd(workspaceVar[target[step]], constantVar[operationConstantIndex[step]]);
+                break;
+            case Operation::ABS:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], fabs);
+                break;
+            case Operation::FLOOR:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], floor);
+                break;
+            case Operation::CEIL:
+                generateSingleArgCall(c, workspaceVar[target[step]], workspaceVar[args[0]], ceil);
+                break;
+            default:
+                // Just invoke evaluateOperation().
+
+                for (int i = 0; i < (int) args.size(); i++)
+                    c.movsd(x86::ptr(argsPointer, 8*i, 0), workspaceVar[args[i]]);
+                X86Gp fn = c.newIntPtr();
+                c.mov(fn, imm_ptr((void*) evaluateOperation));
+                CCFuncCall* call = c.call(fn, FuncSignature2<double, Operation*, double*>());
+                call->setArg(0, imm_ptr(&op));
+                call->setArg(1, imm_ptr(&argValues[0]));
+                call->setRet(0, workspaceVar[target[step]]);
+        }
+    }
+    c.ret(workspaceVar[workspace.size()-1]);
+    c.endFunc();
+    c.finalize();
+    runtime.add(&jitCode, &code);
+}
+
+void CompiledExpression::generateSingleArgCall(X86Compiler& c, X86Xmm& dest, X86Xmm& arg, double (*function)(double)) {
+    X86Gp fn = c.newIntPtr();
+    c.mov(fn, imm_ptr((void*) function));
+    CCFuncCall* call = c.call(fn, FuncSignature1<double, double>());
+    call->setArg(0, arg);
+    call->setRet(0, dest);
+}
+
+void CompiledExpression::generateTwoArgCall(X86Compiler& c, X86Xmm& dest, X86Xmm& arg1, X86Xmm& arg2, double (*function)(double, double)) {
+    X86Gp fn = c.newIntPtr();
+    c.mov(fn, imm_ptr((void*) function));
+    CCFuncCall* call = c.call(fn, FuncSignature2<double, double, double>());
+    call->setArg(0, arg1);
+    call->setArg(1, arg2);
+    call->setRet(0, dest);
+}
+#endif
index 65d3f0c..bbbae85 100644 (file)
@@ -6,7 +6,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-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -84,6 +84,11 @@ const Operation& ExpressionProgram::getOperation(int index) const {
     return *operations[index];
 }
 
+void ExpressionProgram::setOperation(int index, Operation* operation) {
+    delete operations[index];
+    operations[index] = operation;
+}
+
 int ExpressionProgram::getStackSize() const {
     return stackSize;
 }
index eadb20f..c30a8ce 100644 (file)
@@ -3,15 +3,15 @@
 
 /*
  * Up to version 11 (VC++ 2012), Microsoft does not support the
- * standard C99 erf() and erfc() functions so we have to fake them here. 
+ * standard C99 erf() and erfc() functions so we have to fake them here.
  * These were added in version 12 (VC++ 2013), which sets _MSC_VER=1800
  * (VC11 has _MSC_VER=1700).
  */
 
-#if defined(_MSC_VER) 
+#if defined(_MSC_VER)
 #define M_PI 3.14159265358979323846264338327950288
 
-#if _MSC_VER <= 1700 // 1700 is VC11, 1800 is VC12 
+#if _MSC_VER <= 1700 // 1700 is VC11, 1800 is VC12
 /***************************
 *   erf.cpp
 *   author:  Steve Strand
index 693dea2..78741c4 100644 (file)
@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -202,6 +202,16 @@ ExpressionTreeNode Operation::Atan::differentiate(const std::vector<ExpressionTr
                               childDerivs[0]);
 }
 
+ExpressionTreeNode Operation::Atan2::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
+    return ExpressionTreeNode(new Operation::Divide(),
+                              ExpressionTreeNode(new Operation::Subtract(),
+                                                 ExpressionTreeNode(new Operation::Multiply(), children[1], childDerivs[0]),
+                                                 ExpressionTreeNode(new Operation::Multiply(), children[0], childDerivs[1])),
+                              ExpressionTreeNode(new Operation::Add(),
+                                                 ExpressionTreeNode(new Operation::Square(), children[0]),
+                                                 ExpressionTreeNode(new Operation::Square(), children[1])));
+}
+
 ExpressionTreeNode Operation::Sinh::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
     return ExpressionTreeNode(new Operation::Multiply(),
                               ExpressionTreeNode(new Operation::Cosh(),
index 6effd06..fd3b091 100644 (file)
@@ -271,6 +271,16 @@ ExpressionTreeNode ParsedExpression::substituteSimplerExpression(const Expressio
                 return ExpressionTreeNode(new Operation::MultiplyConstant(-dynamic_cast<const Operation::MultiplyConstant*>(&node.getOperation())->getValue()), children[0].getChildren()[0]);
             break;
         }
+        case Operation::SQRT:
+        {
+            if (children[0].getOperation().getId() == Operation::SQUARE) // sqrt(square(x)) = abs(x)
+                return ExpressionTreeNode(new Operation::Abs(), children[0].getChildren()[0]);
+        }
+        case Operation::SQUARE:
+        {
+            if (children[0].getOperation().getId() == Operation::SQRT) // square(sqrt(x)) = x
+                return children[0].getChildren()[0];
+        }
         default:
         {
             // If operation ID is not one of the above,
index 6b19d73..e284add 100644 (file)
@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2009-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -66,7 +66,7 @@ private:
 
 string Parser::trim(const string& expression) {
     // Remove leading and trailing spaces.
-    
+
     int start, end;
     for (start = 0; start < (int) expression.size() && isspace(expression[start]); start++)
         ;
@@ -313,6 +313,7 @@ Operation* Parser::getFunctionOperation(const std::string& name, const map<strin
         opMap["asin"] = Operation::ASIN;
         opMap["acos"] = Operation::ACOS;
         opMap["atan"] = Operation::ATAN;
+        opMap["atan2"] = Operation::ATAN2;
         opMap["sinh"] = Operation::SINH;
         opMap["cosh"] = Operation::COSH;
         opMap["tanh"] = Operation::TANH;
@@ -368,6 +369,8 @@ Operation* Parser::getFunctionOperation(const std::string& name, const map<strin
             return new Operation::Acos();
         case Operation::ATAN:
             return new Operation::Atan();
+        case Operation::ATAN2:
+            return new Operation::Atan2();
         case Operation::SINH:
             return new Operation::Sinh();
         case Operation::COSH: