jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Upgrade Lepton library to version 2019-06-02

This commit is contained in:
Giacomo Fiorin 2019-10-14 14:59:34 -04:00
parent 8f431b0fb8
commit e04a18fc4f
11 changed files with 175 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
lib/colvars/README
View File

@ -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

11
lib/colvars/lepton/include/lepton/CompiledExpression.h
View 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

32
lib/colvars/lepton/include/lepton/CustomFunction.h
View 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_*/

10
lib/colvars/lepton/include/lepton/ExpressionProgram.h
View 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.

34
lib/colvars/lepton/include/lepton/Operation.h
View 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) {

100
lib/colvars/lepton/src/CompiledExpression.cpp
View 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) 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<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());
@ -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<string, double*>& 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<string, int>::const_iterator iter = variableIndices.begin(); iter != variableIndices.end(); ++iter) {
map<string, double*>::iterator pointer = variablePointers.find(iter->first);
@ -165,13 +165,13 @@ void CompiledExpression::setVariableLocations(map<string, double*>& 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<int>& 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<string, double>* dummyVariables = NULL;
return op->evaluate(args, *dummyVariables);
static map<string, double> dummyVariables;
return op->evaluate(args, dummyVariables);
}
void CompiledExpression::generateJitCode() {
X86Compiler c(&runtime);
c.addFunc(kFuncConvHost, FuncBuilder0<double>());
vector<X86XmmVar> workspaceVar(workspace.size());
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.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<string>::const_iterator iter = variableNames.begin(); iter != variableNames.end(); ++iter) {
map<string, int>::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<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)
@ -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<X86XmmVar> constantVar(constants.size());
vector<X86Xmm> 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<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]]);
@ -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<double, Operation*, double*>());
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]]);
@ -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<double, double>());
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

7
lib/colvars/lepton/src/ExpressionProgram.cpp
View 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;
}

6
lib/colvars/lepton/src/MSVC_erfc.h
View 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

12
lib/colvars/lepton/src/Operation.cpp
View 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(),

10
lib/colvars/lepton/src/ParsedExpression.cpp
View 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,

7
lib/colvars/lepton/src/Parser.cpp
View 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: