diff --git a/lib/colvars/README b/lib/colvars/README index 087528748b..65821d68c7 100644 --- a/lib/colvars/README +++ b/lib/colvars/README @@ -50,6 +50,9 @@ settings in Makefile.common should work. Note: some Colvars functions use the Lepton mathematical expression parser, which is here included (no additional steps required). For more details, see: https://simtk.org/projects/lepton +Starting from 2019-06-02, Lepton requires C++11. Please see: + https://colvars.github.io/README-c++11.html + https://lammps.sandia.gov/doc/Build_settings.html#cxx11 ## Documentation diff --git a/lib/colvars/lepton/include/lepton/CompiledExpression.h b/lib/colvars/lepton/include/lepton/CompiledExpression.h index 67442e0cf5..c7e393e93b 100644 --- a/lib/colvars/lepton/include/lepton/CompiledExpression.h +++ b/lib/colvars/lepton/include/lepton/CompiledExpression.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) 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 workspace; mutable std::vector argValues; std::map 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 constants; asmjit::JitRuntime runtime; #endif diff --git a/lib/colvars/lepton/include/lepton/CustomFunction.h b/lib/colvars/lepton/include/lepton/CustomFunction.h index 5c5586105f..fbb0ddd52a 100644 --- a/lib/colvars/lepton/include/lepton/CustomFunction.h +++ b/lib/colvars/lepton/include/lepton/CustomFunction.h @@ -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_*/ diff --git a/lib/colvars/lepton/include/lepton/ExpressionProgram.h b/lib/colvars/lepton/include/lepton/ExpressionProgram.h index 94d37f471d..a49a9094d0 100644 --- a/lib/colvars/lepton/include/lepton/ExpressionProgram.h +++ b/lib/colvars/lepton/include/lepton/ExpressionProgram.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 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. diff --git a/lib/colvars/lepton/include/lepton/Operation.h b/lib/colvars/lepton/include/lepton/Operation.h index f7a8b78163..1ddde0b8c0 100644 --- a/lib/colvars/lepton/include/lepton/Operation.h +++ b/lib/colvars/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-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& 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& children, const std::vector& 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& variables) const { + return std::atan2(args[0], args[1]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& 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& 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) { diff --git a/lib/colvars/lepton/src/CompiledExpression.cpp b/lib/colvars/lepton/src/CompiledExpression.cpp index 302f294ee2..1ad348b47d 100644 --- a/lib/colvars/lepton/src/CompiledExpression.cpp +++ b/lib/colvars/lepton/src/CompiledExpression.cpp @@ -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) 2013-2016 Stanford University and the Authors. * + * Portions copyright (c) 2013-2019 Stanford University and the Authors. * * Authors: Peter Eastman * * Contributors: * * * @@ -84,17 +84,17 @@ CompiledExpression& CompiledExpression::operator=(const CompiledExpression& expr void CompiledExpression::compileExpression(const ExpressionTreeNode& node, vector >& temps) { if (findTempIndex(node, temps) != -1) return; // We have already processed a node identical to this one. - + // Process the child nodes. - + vector 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()); @@ -108,7 +108,7 @@ void CompiledExpression::compileExpression(const ExpressionTreeNode& node, vecto 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) @@ -148,12 +148,12 @@ void CompiledExpression::setVariableLocations(map& variableLoca 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::const_iterator iter = variableIndices.begin(); iter != variableIndices.end(); ++iter) { map::iterator pointer = variablePointers.find(iter->first); @@ -165,13 +165,13 @@ void CompiledExpression::setVariableLocations(map& variableLoca double CompiledExpression::evaluate() const { #ifdef LEPTON_USE_JIT - return ((double (*)()) jitCode)(); + 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& args = arguments[step]; if (args.size() == 1) @@ -188,34 +188,36 @@ double CompiledExpression::evaluate() const { #ifdef LEPTON_USE_JIT static double evaluateOperation(Operation* op, double* args) { - map* dummyVariables = NULL; - return op->evaluate(args, *dummyVariables); + static map dummyVariables; + return op->evaluate(args, dummyVariables); } void CompiledExpression::generateJitCode() { - X86Compiler c(&runtime); - c.addFunc(kFuncConvHost, FuncBuilder0()); - vector workspaceVar(workspace.size()); + CodeHolder code; + code.init(runtime.getCodeInfo()); + X86Compiler c(&code); + c.addFunc(FuncSignature0()); + vector workspaceVar(workspace.size()); for (int i = 0; i < (int) workspaceVar.size(); i++) - workspaceVar[i] = c.newXmmVar(kX86VarTypeXmmSd); - X86GpVar argsPointer(c); + workspaceVar[i] = c.newXmmSd(); + X86Gp argsPointer = c.newIntPtr(); c.mov(argsPointer, imm_ptr(&argValues[0])); - + // Load the arguments into variables. - + for (set::const_iterator iter = variableNames.begin(); iter != variableNames.end(); ++iter) { map::iterator index = variableIndices.find(*iter); - X86GpVar variablePointer(c); + 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 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) @@ -232,9 +234,9 @@ void CompiledExpression::generateJitCode() { 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; @@ -245,33 +247,33 @@ void CompiledExpression::generateJitCode() { constants.push_back(value); } } - + // Load constants into variables. - - vector constantVar(constants.size()); + + vector constantVar(constants.size()); if (constants.size() > 0) { - X86GpVar constantsPointer(c); + X86Gp constantsPointer = c.newIntPtr(); c.mov(constantsPointer, imm_ptr(&constants[0])); for (int i = 0; i < (int) constants.size(); i++) { - constantVar[i] = c.newXmmVar(kX86VarTypeXmmSd); + 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 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]]); @@ -292,6 +294,9 @@ void CompiledExpression::generateJitCode() { 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]]); @@ -323,6 +328,9 @@ void CompiledExpression::generateJitCode() { 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; @@ -374,12 +382,12 @@ void CompiledExpression::generateJitCode() { 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]]); - X86GpVar fn(c, kVarTypeIntPtr); + X86Gp fn = c.newIntPtr(); c.mov(fn, imm_ptr((void*) evaluateOperation)); - X86CallNode* call = c.call(fn, kFuncConvHost, FuncBuilder2()); + CCFuncCall* call = c.call(fn, FuncSignature2()); call->setArg(0, imm_ptr(&op)); call->setArg(1, imm_ptr(&argValues[0])); call->setRet(0, workspaceVar[target[step]]); @@ -387,14 +395,24 @@ void CompiledExpression::generateJitCode() { } c.ret(workspaceVar[workspace.size()-1]); c.endFunc(); - jitCode = c.make(); + c.finalize(); + runtime.add(&jitCode, &code); } -void CompiledExpression::generateSingleArgCall(X86Compiler& c, X86XmmVar& dest, X86XmmVar& arg, double (*function)(double)) { - X86GpVar fn(c, kVarTypeIntPtr); +void CompiledExpression::generateSingleArgCall(X86Compiler& c, X86Xmm& dest, X86Xmm& arg, double (*function)(double)) { + X86Gp fn = c.newIntPtr(); c.mov(fn, imm_ptr((void*) function)); - X86CallNode* call = c.call(fn, kFuncConvHost, FuncBuilder1()); + CCFuncCall* call = c.call(fn, FuncSignature1()); 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()); + call->setArg(0, arg1); + call->setArg(1, arg2); + call->setRet(0, dest); +} #endif diff --git a/lib/colvars/lepton/src/ExpressionProgram.cpp b/lib/colvars/lepton/src/ExpressionProgram.cpp index 65d3f0c79a..bbbae8533f 100644 --- a/lib/colvars/lepton/src/ExpressionProgram.cpp +++ b/lib/colvars/lepton/src/ExpressionProgram.cpp @@ -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; } diff --git a/lib/colvars/lepton/src/MSVC_erfc.h b/lib/colvars/lepton/src/MSVC_erfc.h index eadb20fdf8..c30a8ce542 100644 --- a/lib/colvars/lepton/src/MSVC_erfc.h +++ b/lib/colvars/lepton/src/MSVC_erfc.h @@ -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 diff --git a/lib/colvars/lepton/src/Operation.cpp b/lib/colvars/lepton/src/Operation.cpp index 693dea2ede..78741c4814 100644 --- a/lib/colvars/lepton/src/Operation.cpp +++ b/lib/colvars/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-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& children, const std::vector& 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& children, const std::vector& childDerivs, const std::string& variable) const { return ExpressionTreeNode(new Operation::Multiply(), ExpressionTreeNode(new Operation::Cosh(), diff --git a/lib/colvars/lepton/src/ParsedExpression.cpp b/lib/colvars/lepton/src/ParsedExpression.cpp index 6effd06007..fd3b091d3c 100644 --- a/lib/colvars/lepton/src/ParsedExpression.cpp +++ b/lib/colvars/lepton/src/ParsedExpression.cpp @@ -271,6 +271,16 @@ ExpressionTreeNode ParsedExpression::substituteSimplerExpression(const Expressio return ExpressionTreeNode(new Operation::MultiplyConstant(-dynamic_cast(&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, diff --git a/lib/colvars/lepton/src/Parser.cpp b/lib/colvars/lepton/src/Parser.cpp index 6b19d7370d..e284add258 100644 --- a/lib/colvars/lepton/src/Parser.cpp +++ b/lib/colvars/lepton/src/Parser.cpp @@ -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