ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

import Lepton library with namespace and header changed to LMP_Lepton

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-12-21 14:18:39 -05:00
parent d8620fc34c
commit 517a2e5e26
24 changed files with 3975 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
lib/README
View File

@ -33,6 +33,8 @@ kim hooks to the KIM library, used by KIM package
from Ryan Elliott and Ellad Tadmor (U Minn)
kokkos Kokkos package for GPU and many-core acceleration
from Kokkos development team (Sandia)
lepton Lepton library for fast evaluation of mathematical
expressions from a string. Imported from OpenMM.
linalg set of BLAS and LAPACK routines needed by ATC package
from Axel Kohlmeyer (Temple U)
mdi hooks to the MDI library, used by MDI package

1
lib/lepton/Install.py Symbolic link
View File

@ -0,0 +1 @@
../Install.py

20
lib/lepton/LICENSE Normal file
View File

@ -0,0 +1,20 @@
Portions copyright (c) 2009-2019 Stanford University and the Authors.
Authors: Peter Eastman and OpenMM 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.

7
lib/lepton/Makefile.lammps.empty Normal file
View File

@ -0,0 +1,7 @@
# Settings that the LAMMPS build will import when this package library is used
# The default settings assume that HDF5 support is integrated into the standard
# distribution and search paths and thus only needs to link the HDF5 library.
lepton_SYSINC =
lepton_SYSLIB =
lepton_SYSPATH =

34
lib/lepton/Makefile.mpi Normal file
View File

@ -0,0 +1,34 @@
EXTRAMAKE=Makefile.lammps.empty
CC=mpicxx
# -DH5_NO_DEPRECATED_SYMBOLS is required here to ensure we are using
# the v1.8 API when HDF5 is configured to default to using the v1.6 API.
CXXFLAGS=-D_DEFAULT_SOURCE -O2 -Wall -fPIC -std=c++11
INC=-I include
AR=ar
ARFLAGS=rc
# need to build two libraries to not break compatibility and to support Install.py
LIB=liblepton.a
SRC=src/CompiledExpression.cpp src/ExpressionProgram.cpp src/ExpressionTreeNode.cpp src/Operation.cpp src/ParsedExpression.cpp src/Parser.cpp
OBJ=$(SRC:src/%.cpp=build/%.o)
all: $(LIB) Makefile.lammps
build:
mkdir -p build
build/%.o: src/%.cpp build
$(CXX) $(INC) $(CXXFLAGS) -c $< -o $@
Makefile.lammps:
cp $(EXTRAMAKE) $@
.PHONY: all lib clean
$(LIB) : $(OBJ)
$(AR) $(ARFLAGS) $@ $^
clean:
rm -f build/*.o $(LIB) *~

34
lib/lepton/Makefile.serial Normal file
View File

@ -0,0 +1,34 @@
EXTRAMAKE=Makefile.lammps.empty
CC=g++
# -DH5_NO_DEPRECATED_SYMBOLS is required here to ensure we are using
# the v1.8 API when HDF5 is configured to default to using the v1.6 API.
CXXFLAGS=-D_DEFAULT_SOURCE -O2 -Wall -fPIC -std=c++11
INC=-I include
AR=ar
ARFLAGS=rc
# need to build two libraries to not break compatibility and to support Install.py
LIB=liblepton.a
SRC=src/CompiledExpression.cpp src/ExpressionProgram.cpp src/ExpressionTreeNode.cpp src/Operation.cpp src/ParsedExpression.cpp src/Parser.cpp
OBJ=$(SRC:src/%.cpp=build/%.o)
all: $(LIB) Makefile.lammps
build:
mkdir -p build
build/%.o: src/%.cpp build
$(CXX) $(INC) $(CXXFLAGS) -c $< -o $@
Makefile.lammps:
cp $(EXTRAMAKE) $@
.PHONY: all lib clean
$(LIB) : $(OBJ)
$(AR) $(ARFLAGS) $@ $^
clean:
rm -f build/*.o $(LIB) *~

43
lib/lepton/README.md Normal file
View File

@ -0,0 +1,43 @@
This directory contains the lepton library from the OpenMM software
which allows to efficiently evaluate mathematical expressions from
strings. This library is used with the LEPTON package that support
force styles within LAMMPS that make use of this library.
You can type "make lib-lepton" from the src directory to see help on how
to build this library via make commands, or you can do the same thing
by typing "python Install.py" from within this directory, or you can
do it manually by following the instructions below.
---------------------
Lepton (short for “lightweight expression parser”) is a C++ library for
parsing, evaluating, differentiating, and analyzing mathematical
expressions. It takes expressions in the form of text strings, then
converts them to an internal representation suitable for evaluation or
analysis. Here are some of its major features:
- Support for a large number of standard mathematical functions and operations.
- Support for user defined custom functions.
- A variety of optimizations for automatically simplifying expressions.
- Computing analytic derivatives.
- Representing parsed expressions in two different forms (tree or program) suitable for
further analysis or processing.
Lepton was originally created for use in the [OpenMM project](https://openmm.org)
ch5md is developed by Pierre de Buyl and is released under the 3-clause BSD
license that can be found in the file LICENSE.
To use the h5md dump style in lammps, execute
make -f Makefile.h5cc
in this directory then
make yes-h5md
in the src directory of LAMMPS to rebuild LAMMPS.
Note that you must have the h5cc compiler installed to use
Makefile.h5cc. It should be part
If HDF5 is not in a standard system location, edit Makefile.lammps accordingly.
In the case of 2015 and more recent debian and ubuntu systems where concurrent
serial and mpi are possible, use the full platform depedent path, i.e.
`HDF5_PATH=/usr/lib/x86_64-linux-gnu/hdf5/serial`

43
lib/lepton/include/LMP_Lepton.h Normal file
View File

@ -0,0 +1,43 @@
#ifndef LMP_LEPTON_H_
#define LMP_LEPTON_H_
/* -------------------------------------------------------------------------- *
* 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) 2009 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/CustomFunction.h"
#include "lepton/ExpressionProgram.h"
#include "lepton/ExpressionTreeNode.h"
#include "lepton/Operation.h"
#include "lepton/ParsedExpression.h"
#include "lepton/Parser.h"
#endif /*LMP_LEPTON_H_*/

114
lib/lepton/include/lepton/CompiledExpression.h Normal file
View File

@ -0,0 +1,114 @@
#ifndef LEPTON_COMPILED_EXPRESSION_H_
#define LEPTON_COMPILED_EXPRESSION_H_
/* -------------------------------------------------------------------------- *
* 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 "ExpressionTreeNode.h"
#include "windowsIncludes.h"
#include <map>
#include <set>
#include <string>
#include <utility>
#include <vector>
#ifdef LEPTON_USE_JIT
#include "asmjit.h"
#endif
namespace LMP_Lepton {
class Operation;
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.
*/
class LEPTON_EXPORT CompiledExpression {
public:
CompiledExpression();
CompiledExpression(const CompiledExpression& expression);
~CompiledExpression();
CompiledExpression& operator=(const CompiledExpression& expression);
/**
* Get the names of all variables used by this expression.
*/
const std::set<std::string>& getVariables() const;
/**
* Get a reference to the memory location where the value of a particular variable is stored. This can be used
* to set the value of the variable before calling evaluate().
*/
double& getVariableReference(const std::string& name);
/**
* You can optionally specify the memory locations from which the values of variables should be read.
* This is useful, for example, when several expressions all use the same variable. You can then set
* the value of that variable in one place, and it will be seen by all of them.
*/
void setVariableLocations(std::map<std::string, double*>& variableLocations);
/**
* Evaluate the expression. The values of all variables should have been set before calling this.
*/
double evaluate() const;
private:
friend class ParsedExpression;
CompiledExpression(const ParsedExpression& expression);
void compileExpression(const ExpressionTreeNode& node, std::vector<std::pair<ExpressionTreeNode, int> >& temps);
int findTempIndex(const ExpressionTreeNode& node, std::vector<std::pair<ExpressionTreeNode, int> >& temps);
std::map<std::string, double*> variablePointers;
std::vector<std::pair<double*, double*> > variablesToCopy;
std::vector<std::vector<int> > arguments;
std::vector<int> target;
std::vector<Operation*> operation;
std::map<std::string, int> variableIndices;
std::set<std::string> variableNames;
mutable std::vector<double> workspace;
mutable std::vector<double> argValues;
std::map<std::string, double> dummyVariables;
double (*jitCode)();
#ifdef LEPTON_USE_JIT
void generateJitCode();
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
};
} // namespace LMP_Lepton
#endif /*LEPTON_COMPILED_EXPRESSION_H_*/

109
lib/lepton/include/lepton/CustomFunction.h Normal file
View File

@ -0,0 +1,109 @@
#ifndef LEPTON_CUSTOM_FUNCTION_H_
#define LEPTON_CUSTOM_FUNCTION_H_
/* -------------------------------------------------------------------------- *
* 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) 2009 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 "windowsIncludes.h"
namespace LMP_Lepton {
/**
* This class is the interface for defining your own function that may be included in expressions.
* To use it, create a concrete subclass that implements all of the virtual methods for each new function
* you want to define. Then when you call Parser::parse() to parse an expression, pass a map of
* function names to CustomFunction objects.
*/
class LEPTON_EXPORT CustomFunction {
public:
virtual ~CustomFunction() {
}
/**
* Get the number of arguments this function expects.
*/
virtual int getNumArguments() const = 0;
/**
* Evaluate the function.
*
* @param arguments the array of argument values
*/
virtual double evaluate(const double* arguments) const = 0;
/**
* Evaluate a derivative of the function.
*
* @param arguments the array of argument values
* @param derivOrder an array specifying the number of times the function has been differentiated
* with respect to each of its arguments. For example, the array {0, 2} indicates
* a second derivative with respect to the second argument.
*/
virtual double evaluateDerivative(const double* arguments, const int* derivOrder) const = 0;
/**
* Create a new duplicate of this object on the heap using the "new" operator.
*/
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 LMP_Lepton
#endif /*LEPTON_CUSTOM_FUNCTION_H_*/

59
lib/lepton/include/lepton/Exception.h Normal file
View File

@ -0,0 +1,59 @@
#ifndef LEPTON_EXCEPTION_H_
#define LEPTON_EXCEPTION_H_
/* -------------------------------------------------------------------------- *
* 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) 2009 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 <exception>
#include <string>
namespace LMP_Lepton {
/**
* This class is used for all exceptions thrown by Lepton.
*/
class Exception : public std::exception {
public:
Exception(const std::string& message) : message(message) {
}
~Exception() throw() {
}
const char* what() const throw() {
return message.c_str();
}
private:
std::string message;
};
} // namespace LMP_Lepton
#endif /*LEPTON_EXCEPTION_H_*/

103
lib/lepton/include/lepton/ExpressionProgram.h Normal file
View File

@ -0,0 +1,103 @@
#ifndef LEPTON_EXPRESSION_PROGRAM_H_
#define LEPTON_EXPRESSION_PROGRAM_H_
/* -------------------------------------------------------------------------- *
* 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) 2009-2018 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 "ExpressionTreeNode.h"
#include "windowsIncludes.h"
#include <map>
#include <string>
#include <vector>
namespace LMP_Lepton {
class ParsedExpression;
/**
* An ExpressionProgram is a linear sequence of Operations for evaluating an expression. The evaluation
* is done with a stack. The arguments to each Operation are first taken off the stack in order, then it is
* evaluated and the result is pushed back onto the stack. At the end, the stack contains a single value,
* which is the value of the expression.
*
* An ExpressionProgram is created by calling createProgram() on a ParsedExpression.
*/
class LEPTON_EXPORT ExpressionProgram {
public:
ExpressionProgram();
ExpressionProgram(const ExpressionProgram& program);
~ExpressionProgram();
ExpressionProgram& operator=(const ExpressionProgram& program);
/**
* Get the number of Operations that make up this program.
*/
int getNumOperations() const;
/**
* 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.
*/
int getStackSize() const;
/**
* Evaluate the expression. If the expression involves any variables, this method will throw an exception.
*/
double evaluate() const;
/**
* Evaluate the expression.
*
* @param variables a map specifying the values of all variables that appear in the expression. If any
* variable appears in the expression but is not included in this map, an exception
* will be thrown.
*/
double evaluate(const std::map<std::string, double>& variables) const;
private:
friend class ParsedExpression;
ExpressionProgram(const ParsedExpression& expression);
void buildProgram(const ExpressionTreeNode& node);
std::vector<Operation*> operations;
int maxArgs, stackSize;
};
} // namespace LMP_Lepton
#endif /*LEPTON_EXPRESSION_PROGRAM_H_*/

105
lib/lepton/include/lepton/ExpressionTreeNode.h Normal file
View File

@ -0,0 +1,105 @@
#ifndef LEPTON_EXPRESSION_TREE_NODE_H_
#define LEPTON_EXPRESSION_TREE_NODE_H_
/* -------------------------------------------------------------------------- *
* 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) 2009 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 "windowsIncludes.h"
#include <string>
#include <vector>
namespace LMP_Lepton {
class Operation;
/**
* This class represents a node in the abstract syntax tree representation of an expression.
* Each node is defined by an Operation and a set of children. When the expression is
* evaluated, each child is first evaluated in order, then the resulting values are passed
* as the arguments to the Operation's evaluate() method.
*/
class LEPTON_EXPORT ExpressionTreeNode {
public:
/**
* Create a new ExpressionTreeNode.
*
* @param operation the operation for this node. The ExpressionTreeNode takes over ownership
* of this object, and deletes it when the node is itself deleted.
* @param children the children of this node
*/
ExpressionTreeNode(Operation* operation, const std::vector<ExpressionTreeNode>& children);
/**
* Create a new ExpressionTreeNode with two children.
*
* @param operation the operation for this node. The ExpressionTreeNode takes over ownership
* of this object, and deletes it when the node is itself deleted.
* @param child1 the first child of this node
* @param child2 the second child of this node
*/
ExpressionTreeNode(Operation* operation, const ExpressionTreeNode& child1, const ExpressionTreeNode& child2);
/**
* Create a new ExpressionTreeNode with one child.
*
* @param operation the operation for this node. The ExpressionTreeNode takes over ownership
* of this object, and deletes it when the node is itself deleted.
* @param child the child of this node
*/
ExpressionTreeNode(Operation* operation, const ExpressionTreeNode& child);
/**
* Create a new ExpressionTreeNode with no children.
*
* @param operation the operation for this node. The ExpressionTreeNode takes over ownership
* of this object, and deletes it when the node is itself deleted.
*/
ExpressionTreeNode(Operation* operation);
ExpressionTreeNode(const ExpressionTreeNode& node);
ExpressionTreeNode();
~ExpressionTreeNode();
bool operator==(const ExpressionTreeNode& node) const;
bool operator!=(const ExpressionTreeNode& node) const;
ExpressionTreeNode& operator=(const ExpressionTreeNode& node);
/**
* Get the Operation performed by this node.
*/
const Operation& getOperation() const;
/**
* Get this node's child nodes.
*/
const std::vector<ExpressionTreeNode>& getChildren() const;
private:
Operation* operation;
std::vector<ExpressionTreeNode> children;
};
} // namespace LMP_Lepton
#endif /*LEPTON_EXPRESSION_TREE_NODE_H_*/

1193
lib/lepton/include/lepton/Operation.h Normal file
View File

File diff suppressed because it is too large Load Diff

131
lib/lepton/include/lepton/ParsedExpression.h Normal file
View File

@ -0,0 +1,131 @@
#ifndef LEPTON_PARSED_EXPRESSION_H_
#define LEPTON_PARSED_EXPRESSION_H_
/* -------------------------------------------------------------------------- *
* 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) 2009=2013 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 "ExpressionTreeNode.h"
#include "windowsIncludes.h"
#include <map>
#include <string>
namespace LMP_Lepton {
class CompiledExpression;
class ExpressionProgram;
/**
* This class represents the result of parsing an expression. It provides methods for working with the
* expression in various ways, such as evaluating it, getting the tree representation of the expresson, etc.
*/
class LEPTON_EXPORT ParsedExpression {
public:
/**
* Create an uninitialized ParsedExpression. This exists so that ParsedExpressions can be put in STL containers.
* Doing anything with it will produce an exception.
*/
ParsedExpression();
/**
* Create a ParsedExpression. Normally you will not call this directly. Instead, use the Parser class
* to parse expression.
*/
ParsedExpression(const ExpressionTreeNode& rootNode);
/**
* Get the root node of the expression's abstract syntax tree.
*/
const ExpressionTreeNode& getRootNode() const;
/**
* Evaluate the expression. If the expression involves any variables, this method will throw an exception.
*/
double evaluate() const;
/**
* Evaluate the expression.
*
* @param variables a map specifying the values of all variables that appear in the expression. If any
* variable appears in the expression but is not included in this map, an exception
* will be thrown.
*/
double evaluate(const std::map<std::string, double>& variables) const;
/**
* Create a new ParsedExpression which produces the same result as this one, but is faster to evaluate.
*/
ParsedExpression optimize() const;
/**
* Create a new ParsedExpression which produces the same result as this one, but is faster to evaluate.
*
* @param variables a map specifying values for a subset of variables that appear in the expression.
* All occurrences of these variables in the expression are replaced with the values
* specified.
*/
ParsedExpression optimize(const std::map<std::string, double>& variables) const;
/**
* Create a new ParsedExpression which is the analytic derivative of this expression with respect to a
* particular variable.
*
* @param variable the variable with respect to which the derivate should be taken
*/
ParsedExpression differentiate(const std::string& variable) const;
/**
* Create an ExpressionProgram that represents the same calculation as this expression.
*/
ExpressionProgram createProgram() const;
/**
* Create a CompiledExpression that represents the same calculation as this expression.
*/
CompiledExpression createCompiledExpression() const;
/**
* Create a new ParsedExpression which is identical to this one, except that the names of some
* variables have been changed.
*
* @param replacements a map whose keys are the names of variables, and whose values are the
* new names to replace them with
*/
ParsedExpression renameVariables(const std::map<std::string, std::string>& replacements) const;
private:
static double evaluate(const ExpressionTreeNode& node, const std::map<std::string, double>& variables);
static ExpressionTreeNode preevaluateVariables(const ExpressionTreeNode& node, const std::map<std::string, double>& variables);
static ExpressionTreeNode precalculateConstantSubexpressions(const ExpressionTreeNode& node);
static ExpressionTreeNode substituteSimplerExpression(const ExpressionTreeNode& node);
static ExpressionTreeNode differentiate(const ExpressionTreeNode& node, const std::string& variable);
static bool isConstant(const ExpressionTreeNode& node);
static double getConstantValue(const ExpressionTreeNode& node);
static ExpressionTreeNode renameNodeVariables(const ExpressionTreeNode& node, const std::map<std::string, std::string>& replacements);
ExpressionTreeNode rootNode;
};
LEPTON_EXPORT std::ostream& operator<<(std::ostream& out, const ExpressionTreeNode& node);
LEPTON_EXPORT std::ostream& operator<<(std::ostream& out, const ParsedExpression& exp);
} // namespace LMP_Lepton
#endif /*LEPTON_PARSED_EXPRESSION_H_*/

77
lib/lepton/include/lepton/Parser.h Normal file
View File

@ -0,0 +1,77 @@
#ifndef LEPTON_PARSER_H_
#define LEPTON_PARSER_H_
/* -------------------------------------------------------------------------- *
* 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) 2009 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 "windowsIncludes.h"
#include <map>
#include <string>
#include <vector>
namespace LMP_Lepton {
class CustomFunction;
class ExpressionTreeNode;
class Operation;
class ParsedExpression;
class ParseToken;
/**
* This class provides the main interface for parsing expressions.
*/
class LEPTON_EXPORT Parser {
public:
/**
* Parse a mathematical expression and return a representation of it as an abstract syntax tree.
*/
static ParsedExpression parse(const std::string& expression);
/**
* Parse a mathematical expression and return a representation of it as an abstract syntax tree.
*
* @param customFunctions a map specifying user defined functions that may appear in the expression.
* The key are function names, and the values are corresponding CustomFunction objects.
*/
static ParsedExpression parse(const std::string& expression, const std::map<std::string, CustomFunction*>& customFunctions);
private:
static std::string trim(const std::string& expression);
static std::vector<ParseToken> tokenize(const std::string& expression);
static ParseToken getNextToken(const std::string& expression, int start);
static ExpressionTreeNode parsePrecedence(const std::vector<ParseToken>& tokens, int& pos, const std::map<std::string, CustomFunction*>& customFunctions,
const std::map<std::string, ExpressionTreeNode>& subexpressionDefs, int precedence);
static Operation* getOperatorOperation(const std::string& name);
static Operation* getFunctionOperation(const std::string& name, const std::map<std::string, CustomFunction*>& customFunctions);
};
} // namespace LMP_Lepton
#endif /*LEPTON_PARSER_H_*/

41
lib/lepton/include/lepton/windowsIncludes.h Normal file
View File

@ -0,0 +1,41 @@
#ifndef LEPTON_WINDOW_INCLUDE_H_
#define LEPTON_WINDOW_INCLUDE_H_
/*
* Shared libraries are messy in Visual Studio. We have to distinguish three
* cases:
* (1) this header is being used to build the Lepton shared library
* (dllexport)
* (2) this header is being used by a *client* of the Lepton shared
* library (dllimport)
* (3) we are building the Lepton static library, or the client is
* being compiled with the expectation of linking with the
* Lepton static library (nothing special needed)
* In the CMake script for building this library, we define one of the symbols
* Lepton_BUILDING_{SHARED|STATIC}_LIBRARY
* Client code normally has no special symbol defined, in which case we'll
* assume it wants to use the shared library. However, if the client defines
* the symbol LEPTON_USE_STATIC_LIBRARIES we'll suppress the dllimport so
* that the client code can be linked with static libraries. Note that
* the client symbol is not library dependent, while the library symbols
* affect only the Lepton library, meaning that other libraries can
* be clients of this one. However, we are assuming all-static or all-shared.
*/
#ifdef _MSC_VER
// We don't want to hear about how sprintf is "unsafe".
#pragma warning(disable:4996)
// Keep MS VC++ quiet about lack of dll export of private members.
#pragma warning(disable:4251)
#if defined(LEPTON_BUILDING_SHARED_LIBRARY)
#define LEPTON_EXPORT __declspec(dllexport)
#elif defined(LEPTON_BUILDING_STATIC_LIBRARY) || defined(LEPTON_USE_STATIC_LIBRARIES)
#define LEPTON_EXPORT
#else
#define LEPTON_EXPORT __declspec(dllimport) // i.e., a client of a shared library
#endif
#else
#define LEPTON_EXPORT // Linux, Mac
#endif
#endif // LEPTON_WINDOW_INCLUDE_H_

418
lib/lepton/src/CompiledExpression.cpp Normal file
View File

@ -0,0 +1,418 @@
/* -------------------------------------------------------------------------- *
* 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 LMP_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

110
lib/lepton/src/ExpressionProgram.cpp Normal file
View File

@ -0,0 +1,110 @@
/* -------------------------------------------------------------------------- *
* 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) 2009-2018 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/ExpressionProgram.h"
#include "lepton/Operation.h"
#include "lepton/ParsedExpression.h"
using namespace LMP_Lepton;
using namespace std;
ExpressionProgram::ExpressionProgram() : maxArgs(0), stackSize(0) {
}
ExpressionProgram::ExpressionProgram(const ParsedExpression& expression) : maxArgs(0), stackSize(0) {
buildProgram(expression.getRootNode());
int currentStackSize = 0;
for (int i = 0; i < (int) operations.size(); i++) {
int args = operations[i]->getNumArguments();
if (args > maxArgs)
maxArgs = args;
currentStackSize += 1-args;
if (currentStackSize > stackSize)
stackSize = currentStackSize;
}
}
ExpressionProgram::~ExpressionProgram() {
for (int i = 0; i < (int) operations.size(); i++)
delete operations[i];
}
ExpressionProgram::ExpressionProgram(const ExpressionProgram& program) {
*this = program;
}
ExpressionProgram& ExpressionProgram::operator=(const ExpressionProgram& program) {
maxArgs = program.maxArgs;
stackSize = program.stackSize;
operations.resize(program.operations.size());
for (int i = 0; i < (int) operations.size(); i++)
operations[i] = program.operations[i]->clone();
return *this;
}
void ExpressionProgram::buildProgram(const ExpressionTreeNode& node) {
for (int i = (int) node.getChildren().size()-1; i >= 0; i--)
buildProgram(node.getChildren()[i]);
operations.push_back(node.getOperation().clone());
}
int ExpressionProgram::getNumOperations() const {
return (int) operations.size();
}
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;
}
double ExpressionProgram::evaluate() const {
return evaluate(map<string, double>());
}
double ExpressionProgram::evaluate(const std::map<std::string, double>& variables) const {
vector<double> stack(stackSize+1);
int stackPointer = stackSize;
for (int i = 0; i < (int) operations.size(); i++) {
int numArgs = operations[i]->getNumArguments();
double result = operations[i]->evaluate(&stack[stackPointer], variables);
stackPointer += numArgs-1;
stack[stackPointer] = result;
}
return stack[stackSize-1];
}

107
lib/lepton/src/ExpressionTreeNode.cpp Normal file
View File

@ -0,0 +1,107 @@
/* -------------------------------------------------------------------------- *
* 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) 2009-2015 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/ExpressionTreeNode.h"
#include "lepton/Exception.h"
#include "lepton/Operation.h"
using namespace LMP_Lepton;
using namespace std;
ExpressionTreeNode::ExpressionTreeNode(Operation* operation, const vector<ExpressionTreeNode>& children) : operation(operation), children(children) {
if (operation->getNumArguments() != children.size())
throw Exception("wrong number of arguments to function: "+operation->getName());
}
ExpressionTreeNode::ExpressionTreeNode(Operation* operation, const ExpressionTreeNode& child1, const ExpressionTreeNode& child2) : operation(operation) {
children.push_back(child1);
children.push_back(child2);
if (operation->getNumArguments() != children.size())
throw Exception("wrong number of arguments to function: "+operation->getName());
}
ExpressionTreeNode::ExpressionTreeNode(Operation* operation, const ExpressionTreeNode& child) : operation(operation) {
children.push_back(child);
if (operation->getNumArguments() != children.size())
throw Exception("wrong number of arguments to function: "+operation->getName());
}
ExpressionTreeNode::ExpressionTreeNode(Operation* operation) : operation(operation) {
if (operation->getNumArguments() != children.size())
throw Exception("wrong number of arguments to function: "+operation->getName());
}
ExpressionTreeNode::ExpressionTreeNode(const ExpressionTreeNode& node) : operation(node.operation == NULL ? NULL : node.operation->clone()), children(node.getChildren()) {
}
ExpressionTreeNode::ExpressionTreeNode() : operation(NULL) {
}
ExpressionTreeNode::~ExpressionTreeNode() {
if (operation != NULL)
delete operation;
}
bool ExpressionTreeNode::operator!=(const ExpressionTreeNode& node) const {
if (node.getOperation() != getOperation())
return true;
if (getOperation().isSymmetric() && getChildren().size() == 2) {
if (getChildren()[0] == node.getChildren()[0] && getChildren()[1] == node.getChildren()[1])
return false;
if (getChildren()[0] == node.getChildren()[1] && getChildren()[1] == node.getChildren()[0])
return false;
return true;
}
for (int i = 0; i < (int) getChildren().size(); i++)
if (getChildren()[i] != node.getChildren()[i])
return true;
return false;
}
bool ExpressionTreeNode::operator==(const ExpressionTreeNode& node) const {
return !(*this != node);
}
ExpressionTreeNode& ExpressionTreeNode::operator=(const ExpressionTreeNode& node) {
if (operation != NULL)
delete operation;
operation = node.getOperation().clone();
children = node.getChildren();
return *this;
}
const Operation& ExpressionTreeNode::getOperation() const {
return *operation;
}
const vector<ExpressionTreeNode>& ExpressionTreeNode::getChildren() const {
return children;
}

91
lib/lepton/src/MSVC_erfc.h Normal file
View File

@ -0,0 +1,91 @@
#ifndef LEPTON_MSVC_ERFC_H_
#define LEPTON_MSVC_ERFC_H_
/*
* Up to version 11 (VC++ 2012), Microsoft does not support the
* 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) || defined(__MINGW32__)
#if !defined(M_PI)
#define M_PI 3.14159265358979323846264338327950288
#endif
#endif
#if defined(_MSC_VER)
#if _MSC_VER <= 1700 // 1700 is VC11, 1800 is VC12
/***************************
* erf.cpp
* author: Steve Strand
* written: 29-Jan-04
***************************/
#include <cmath>
static const double rel_error= 1E-12; //calculate 12 significant figures
//you can adjust rel_error to trade off between accuracy and speed
//but don't ask for > 15 figures (assuming usual 52 bit mantissa in a double)
static double erfc(double x);
static double erf(double x)
//erf(x) = 2/sqrt(pi)*integral(exp(-t^2),t,0,x)
// = 2/sqrt(pi)*[x - x^3/3 + x^5/5*2! - x^7/7*3! + ...]
// = 1-erfc(x)
{
static const double two_sqrtpi= 1.128379167095512574; // 2/sqrt(pi)
if (fabs(x) > 2.2) {
return 1.0 - erfc(x); //use continued fraction when fabs(x) > 2.2
}
double sum= x, term= x, xsqr= x*x;
int j= 1;
do {
term*= xsqr/j;
sum-= term/(2*j+1);
++j;
term*= xsqr/j;
sum+= term/(2*j+1);
++j;
} while (fabs(term)/sum > rel_error);
return two_sqrtpi*sum;
}
static double erfc(double x)
//erfc(x) = 2/sqrt(pi)*integral(exp(-t^2),t,x,inf)
// = exp(-x^2)/sqrt(pi) * [1/x+ (1/2)/x+ (2/2)/x+ (3/2)/x+ (4/2)/x+ ...]
// = 1-erf(x)
//expression inside [] is a continued fraction so '+' means add to denominator only
{
static const double one_sqrtpi= 0.564189583547756287; // 1/sqrt(pi)
if (fabs(x) < 2.2) {
return 1.0 - erf(x); //use series when fabs(x) < 2.2
}
// Don't look for x==0 here!
if (x < 0) { //continued fraction only valid for x>0
return 2.0 - erfc(-x);
}
double a=1, b=x; //last two convergent numerators
double c=x, d=x*x+0.5; //last two convergent denominators
double q1, q2= b/d; //last two convergents (a/c and b/d)
double n= 1.0, t;
do {
t= a*n+b*x;
a= b;
b= t;
t= c*n+d*x;
c= d;
d= t;
n+= 0.5;
q1= q2;
q2= b/d;
} while (fabs(q1-q2)/q2 > rel_error);
return one_sqrtpi*exp(-x*x)*q2;
}
#endif // _MSC_VER <= 1700
#endif // _MSC_VER
#endif // LEPTON_MSVC_ERFC_H_

345
lib/lepton/src/Operation.cpp Normal file
View File

@ -0,0 +1,345 @@
/* -------------------------------------------------------------------------- *
* 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) 2009-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/Operation.h"
#include "lepton/ExpressionTreeNode.h"
#include "MSVC_erfc.h"
using namespace LMP_Lepton;
using namespace std;
double Operation::Erf::evaluate(double* args, const map<string, double>& variables) const {
return erf(args[0]);
}
double Operation::Erfc::evaluate(double* args, const map<string, double>& variables) const {
return erfc(args[0]);
}
ExpressionTreeNode Operation::Constant::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Constant(0.0));
}
ExpressionTreeNode Operation::Variable::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
if (variable == name)
return ExpressionTreeNode(new Operation::Constant(1.0));
return ExpressionTreeNode(new Operation::Constant(0.0));
}
ExpressionTreeNode Operation::Custom::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
if (function->getNumArguments() == 0)
return ExpressionTreeNode(new Operation::Constant(0.0));
ExpressionTreeNode result = ExpressionTreeNode(new Operation::Multiply(), ExpressionTreeNode(new Operation::Custom(*this, 0), children), childDerivs[0]);
for (int i = 1; i < getNumArguments(); i++) {
result = ExpressionTreeNode(new Operation::Add(),
result,
ExpressionTreeNode(new Operation::Multiply(), ExpressionTreeNode(new Operation::Custom(*this, i), children), childDerivs[i]));
}
return result;
}
ExpressionTreeNode Operation::Add::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Add(), childDerivs[0], childDerivs[1]);
}
ExpressionTreeNode Operation::Subtract::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Subtract(), childDerivs[0], childDerivs[1]);
}
ExpressionTreeNode Operation::Multiply::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Add(),
ExpressionTreeNode(new Operation::Multiply(), children[0], childDerivs[1]),
ExpressionTreeNode(new Operation::Multiply(), children[1], childDerivs[0]));
}
ExpressionTreeNode Operation::Divide::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::Square(), children[1]));
}
ExpressionTreeNode Operation::Power::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Add(),
ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Multiply(),
children[1],
ExpressionTreeNode(new Operation::Power(),
children[0], ExpressionTreeNode(new Operation::AddConstant(-1.0), children[1]))),
childDerivs[0]),
ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Log(), children[0]),
ExpressionTreeNode(new Operation::Power(), children[0], children[1])),
childDerivs[1]));
}
ExpressionTreeNode Operation::Negate::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Negate(), childDerivs[0]);
}
ExpressionTreeNode Operation::Sqrt::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::MultiplyConstant(0.5),
ExpressionTreeNode(new Operation::Reciprocal(),
ExpressionTreeNode(new Operation::Sqrt(), children[0]))),
childDerivs[0]);
}
ExpressionTreeNode Operation::Exp::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Exp(), children[0]),
childDerivs[0]);
}
ExpressionTreeNode Operation::Log::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Reciprocal(), children[0]),
childDerivs[0]);
}
ExpressionTreeNode Operation::Sin::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Cos(), children[0]),
childDerivs[0]);
}
ExpressionTreeNode Operation::Cos::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Negate(),
ExpressionTreeNode(new Operation::Sin(), children[0])),
childDerivs[0]);
}
ExpressionTreeNode Operation::Sec::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Sec(), children[0]),
ExpressionTreeNode(new Operation::Tan(), children[0])),
childDerivs[0]);
}
ExpressionTreeNode Operation::Csc::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Negate(),
ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Csc(), children[0]),
ExpressionTreeNode(new Operation::Cot(), children[0]))),
childDerivs[0]);
}
ExpressionTreeNode Operation::Tan::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Square(),
ExpressionTreeNode(new Operation::Sec(), children[0])),
childDerivs[0]);
}
ExpressionTreeNode Operation::Cot::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Negate(),
ExpressionTreeNode(new Operation::Square(),
ExpressionTreeNode(new Operation::Csc(), children[0]))),
childDerivs[0]);
}
ExpressionTreeNode Operation::Asin::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Reciprocal(),
ExpressionTreeNode(new Operation::Sqrt(),
ExpressionTreeNode(new Operation::Subtract(),
ExpressionTreeNode(new Operation::Constant(1.0)),
ExpressionTreeNode(new Operation::Square(), children[0])))),
childDerivs[0]);
}
ExpressionTreeNode Operation::Acos::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Negate(),
ExpressionTreeNode(new Operation::Reciprocal(),
ExpressionTreeNode(new Operation::Sqrt(),
ExpressionTreeNode(new Operation::Subtract(),
ExpressionTreeNode(new Operation::Constant(1.0)),
ExpressionTreeNode(new Operation::Square(), children[0]))))),
childDerivs[0]);
}
ExpressionTreeNode Operation::Atan::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Reciprocal(),
ExpressionTreeNode(new Operation::AddConstant(1.0),
ExpressionTreeNode(new Operation::Square(), children[0]))),
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(),
children[0]),
childDerivs[0]);
}
ExpressionTreeNode Operation::Cosh::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Sinh(),
children[0]),
childDerivs[0]);
}
ExpressionTreeNode Operation::Tanh::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Subtract(),
ExpressionTreeNode(new Operation::Constant(1.0)),
ExpressionTreeNode(new Operation::Square(),
ExpressionTreeNode(new Operation::Tanh(), children[0]))),
childDerivs[0]);
}
ExpressionTreeNode Operation::Erf::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Constant(2.0/sqrt(M_PI))),
ExpressionTreeNode(new Operation::Exp(),
ExpressionTreeNode(new Operation::Negate(),
ExpressionTreeNode(new Operation::Square(), children[0])))),
childDerivs[0]);
}
ExpressionTreeNode Operation::Erfc::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Constant(-2.0/sqrt(M_PI))),
ExpressionTreeNode(new Operation::Exp(),
ExpressionTreeNode(new Operation::Negate(),
ExpressionTreeNode(new Operation::Square(), children[0])))),
childDerivs[0]);
}
ExpressionTreeNode Operation::Step::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Constant(0.0));
}
ExpressionTreeNode Operation::Delta::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Constant(0.0));
}
ExpressionTreeNode Operation::Square::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::MultiplyConstant(2.0),
children[0]),
childDerivs[0]);
}
ExpressionTreeNode Operation::Cube::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::MultiplyConstant(3.0),
ExpressionTreeNode(new Operation::Square(), children[0])),
childDerivs[0]);
}
ExpressionTreeNode Operation::Reciprocal::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::Negate(),
ExpressionTreeNode(new Operation::Reciprocal(),
ExpressionTreeNode(new Operation::Square(), children[0]))),
childDerivs[0]);
}
ExpressionTreeNode Operation::AddConstant::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return childDerivs[0];
}
ExpressionTreeNode Operation::MultiplyConstant::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::MultiplyConstant(value),
childDerivs[0]);
}
ExpressionTreeNode Operation::PowerConstant::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Multiply(),
ExpressionTreeNode(new Operation::MultiplyConstant(value),
ExpressionTreeNode(new Operation::PowerConstant(value-1),
children[0])),
childDerivs[0]);
}
ExpressionTreeNode Operation::Min::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
ExpressionTreeNode step(new Operation::Step(),
ExpressionTreeNode(new Operation::Subtract(), children[0], children[1]));
return ExpressionTreeNode(new Operation::Subtract(),
ExpressionTreeNode(new Operation::Multiply(), childDerivs[1], step),
ExpressionTreeNode(new Operation::Multiply(), childDerivs[0],
ExpressionTreeNode(new Operation::AddConstant(-1), step)));
}
ExpressionTreeNode Operation::Max::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
ExpressionTreeNode step(new Operation::Step(),
ExpressionTreeNode(new Operation::Subtract(), children[0], children[1]));
return ExpressionTreeNode(new Operation::Subtract(),
ExpressionTreeNode(new Operation::Multiply(), childDerivs[0], step),
ExpressionTreeNode(new Operation::Multiply(), childDerivs[1],
ExpressionTreeNode(new Operation::AddConstant(-1), step)));
}
ExpressionTreeNode Operation::Abs::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
ExpressionTreeNode step(new Operation::Step(), children[0]);
return ExpressionTreeNode(new Operation::Multiply(),
childDerivs[0],
ExpressionTreeNode(new Operation::AddConstant(-1),
ExpressionTreeNode(new Operation::MultiplyConstant(2), step)));
}
ExpressionTreeNode Operation::Floor::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Constant(0.0));
}
ExpressionTreeNode Operation::Ceil::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
return ExpressionTreeNode(new Operation::Constant(0.0));
}
ExpressionTreeNode Operation::Select::differentiate(const std::vector<ExpressionTreeNode>& children, const std::vector<ExpressionTreeNode>& childDerivs, const std::string& variable) const {
vector<ExpressionTreeNode> derivChildren;
derivChildren.push_back(children[0]);
derivChildren.push_back(childDerivs[1]);
derivChildren.push_back(childDerivs[2]);
return ExpressionTreeNode(new Operation::Select(), derivChildren);
}

379
lib/lepton/src/ParsedExpression.cpp Normal file
View File

@ -0,0 +1,379 @@
/* -------------------------------------------------------------------------- *
* 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) 2009 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/ParsedExpression.h"
#include "lepton/CompiledExpression.h"
#include "lepton/ExpressionProgram.h"
#include "lepton/Operation.h"
#include <limits>
#include <vector>
using namespace LMP_Lepton;
using namespace std;
ParsedExpression::ParsedExpression() : rootNode(ExpressionTreeNode()) {
}
ParsedExpression::ParsedExpression(const ExpressionTreeNode& rootNode) : rootNode(rootNode) {
}
const ExpressionTreeNode& ParsedExpression::getRootNode() const {
if (&rootNode.getOperation() == NULL)
throw Exception("Illegal call to an initialized ParsedExpression");
return rootNode;
}
double ParsedExpression::evaluate() const {
return evaluate(getRootNode(), map<string, double>());
}
double ParsedExpression::evaluate(const map<string, double>& variables) const {
return evaluate(getRootNode(), variables);
}
double ParsedExpression::evaluate(const ExpressionTreeNode& node, const map<string, double>& variables) {
int numArgs = (int) node.getChildren().size();
vector<double> args(max(numArgs, 1));
for (int i = 0; i < numArgs; i++)
args[i] = evaluate(node.getChildren()[i], variables);
return node.getOperation().evaluate(&args[0], variables);
}
ParsedExpression ParsedExpression::optimize() const {
ExpressionTreeNode result = precalculateConstantSubexpressions(getRootNode());
while (true) {
ExpressionTreeNode simplified = substituteSimplerExpression(result);
if (simplified == result)
break;
result = simplified;
}
return ParsedExpression(result);
}
ParsedExpression ParsedExpression::optimize(const map<string, double>& variables) const {
ExpressionTreeNode result = preevaluateVariables(getRootNode(), variables);
result = precalculateConstantSubexpressions(result);
while (true) {
ExpressionTreeNode simplified = substituteSimplerExpression(result);
if (simplified == result)
break;
result = simplified;
}
return ParsedExpression(result);
}
ExpressionTreeNode ParsedExpression::preevaluateVariables(const ExpressionTreeNode& node, const map<string, double>& variables) {
if (node.getOperation().getId() == Operation::VARIABLE) {
const Operation::Variable& var = dynamic_cast<const Operation::Variable&>(node.getOperation());
map<string, double>::const_iterator iter = variables.find(var.getName());
if (iter == variables.end())
return node;
return ExpressionTreeNode(new Operation::Constant(iter->second));
}
vector<ExpressionTreeNode> children(node.getChildren().size());
for (int i = 0; i < (int) children.size(); i++)
children[i] = preevaluateVariables(node.getChildren()[i], variables);
return ExpressionTreeNode(node.getOperation().clone(), children);
}
ExpressionTreeNode ParsedExpression::precalculateConstantSubexpressions(const ExpressionTreeNode& node) {
vector<ExpressionTreeNode> children(node.getChildren().size());
for (int i = 0; i < (int) children.size(); i++)
children[i] = precalculateConstantSubexpressions(node.getChildren()[i]);
ExpressionTreeNode result = ExpressionTreeNode(node.getOperation().clone(), children);
if (node.getOperation().getId() == Operation::VARIABLE || node.getOperation().getId() == Operation::CUSTOM)
return result;
for (int i = 0; i < (int) children.size(); i++)
if (children[i].getOperation().getId() != Operation::CONSTANT)
return result;
return ExpressionTreeNode(new Operation::Constant(evaluate(result, map<string, double>())));
}
ExpressionTreeNode ParsedExpression::substituteSimplerExpression(const ExpressionTreeNode& node) {
vector<ExpressionTreeNode> children(node.getChildren().size());
for (int i = 0; i < (int) children.size(); i++)
children[i] = substituteSimplerExpression(node.getChildren()[i]);
// Collect some info on constant expressions in children
bool first_const = children.size() > 0 && isConstant(children[0]); // is first child constant?
bool second_const = children.size() > 1 && isConstant(children[1]); ; // is second child constant?
double first, second; // if yes, value of first and second child
if (first_const)
first = getConstantValue(children[0]);
if (second_const)
second = getConstantValue(children[1]);
switch (node.getOperation().getId()) {
case Operation::ADD:
{
if (first_const) {
if (first == 0.0) { // Add 0
return children[1];
} else { // Add a constant
return ExpressionTreeNode(new Operation::AddConstant(first), children[1]);
}
}
if (second_const) {
if (second == 0.0) { // Add 0
return children[0];
} else { // Add a constant
return ExpressionTreeNode(new Operation::AddConstant(second), children[0]);
}
}
if (children[1].getOperation().getId() == Operation::NEGATE) // a+(-b) = a-b
return ExpressionTreeNode(new Operation::Subtract(), children[0], children[1].getChildren()[0]);
if (children[0].getOperation().getId() == Operation::NEGATE) // (-a)+b = b-a
return ExpressionTreeNode(new Operation::Subtract(), children[1], children[0].getChildren()[0]);
break;
}
case Operation::SUBTRACT:
{
if (children[0] == children[1])
return ExpressionTreeNode(new Operation::Constant(0.0)); // Subtracting anything from itself is 0
if (first_const) {
if (first == 0.0) // Subtract from 0
return ExpressionTreeNode(new Operation::Negate(), children[1]);
}
if (second_const) {
if (second == 0.0) { // Subtract 0
return children[0];
} else { // Subtract a constant
return ExpressionTreeNode(new Operation::AddConstant(-second), children[0]);
}
}
if (children[1].getOperation().getId() == Operation::NEGATE) // a-(-b) = a+b
return ExpressionTreeNode(new Operation::Add(), children[0], children[1].getChildren()[0]);
break;
}
case Operation::MULTIPLY:
{
if ((first_const && first == 0.0) || (second_const && second == 0.0)) // Multiply by 0
return ExpressionTreeNode(new Operation::Constant(0.0));
if (first_const && first == 1.0) // Multiply by 1
return children[1];
if (second_const && second == 1.0) // Multiply by 1
return children[0];
if (first_const) { // Multiply by a constant
if (children[1].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Combine two multiplies into a single one
return ExpressionTreeNode(new Operation::MultiplyConstant(first*dynamic_cast<const Operation::MultiplyConstant*>(&children[1].getOperation())->getValue()), children[1].getChildren()[0]);
return ExpressionTreeNode(new Operation::MultiplyConstant(first), children[1]);
}
if (second_const) { // Multiply by a constant
if (children[0].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Combine two multiplies into a single one
return ExpressionTreeNode(new Operation::MultiplyConstant(second*dynamic_cast<const Operation::MultiplyConstant*>(&children[0].getOperation())->getValue()), children[0].getChildren()[0]);
return ExpressionTreeNode(new Operation::MultiplyConstant(second), children[0]);
}
if (children[0].getOperation().getId() == Operation::NEGATE && children[1].getOperation().getId() == Operation::NEGATE) // The two negations cancel
return ExpressionTreeNode(new Operation::Multiply(), children[0].getChildren()[0], children[1].getChildren()[0]);
if (children[0].getOperation().getId() == Operation::NEGATE && children[1].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Negate the constant
return ExpressionTreeNode(new Operation::Multiply(), children[0].getChildren()[0], ExpressionTreeNode(new Operation::MultiplyConstant(-dynamic_cast<const Operation::MultiplyConstant*>(&children[1].getOperation())->getValue()), children[1].getChildren()[0]));
if (children[1].getOperation().getId() == Operation::NEGATE && children[0].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Negate the constant
return ExpressionTreeNode(new Operation::Multiply(), ExpressionTreeNode(new Operation::MultiplyConstant(-dynamic_cast<const Operation::MultiplyConstant*>(&children[0].getOperation())->getValue()), children[0].getChildren()[0]), children[1].getChildren()[0]);
if (children[0].getOperation().getId() == Operation::NEGATE) // Pull the negation out so it can possibly be optimized further
return ExpressionTreeNode(new Operation::Negate(), ExpressionTreeNode(new Operation::Multiply(), children[0].getChildren()[0], children[1]));
if (children[1].getOperation().getId() == Operation::NEGATE) // Pull the negation out so it can possibly be optimized further
return ExpressionTreeNode(new Operation::Negate(), ExpressionTreeNode(new Operation::Multiply(), children[0], children[1].getChildren()[0]));
if (children[1].getOperation().getId() == Operation::RECIPROCAL) // a*(1/b) = a/b
return ExpressionTreeNode(new Operation::Divide(), children[0], children[1].getChildren()[0]);
if (children[0].getOperation().getId() == Operation::RECIPROCAL) // (1/a)*b = b/a
return ExpressionTreeNode(new Operation::Divide(), children[1], children[0].getChildren()[0]);
if (children[0] == children[1])
return ExpressionTreeNode(new Operation::Square(), children[0]); // x*x = square(x)
if (children[0].getOperation().getId() == Operation::SQUARE && children[0].getChildren()[0] == children[1])
return ExpressionTreeNode(new Operation::Cube(), children[1]); // x*x*x = cube(x)
if (children[1].getOperation().getId() == Operation::SQUARE && children[1].getChildren()[0] == children[0])
return ExpressionTreeNode(new Operation::Cube(), children[0]); // x*x*x = cube(x)
break;
}
case Operation::DIVIDE:
{
if (children[0] == children[1])
return ExpressionTreeNode(new Operation::Constant(1.0)); // Dividing anything from itself is 0
if (first_const && first == 0.0) // 0 divided by something
return ExpressionTreeNode(new Operation::Constant(0.0));
if (first_const && first == 1.0) // 1 divided by something
return ExpressionTreeNode(new Operation::Reciprocal(), children[1]);
if (second_const && second == 1.0) // Divide by 1
return children[0];
if (second_const) {
if (children[0].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Combine a multiply and a divide into one multiply
return ExpressionTreeNode(new Operation::MultiplyConstant(dynamic_cast<const Operation::MultiplyConstant*>(&children[0].getOperation())->getValue()/second), children[0].getChildren()[0]);
return ExpressionTreeNode(new Operation::MultiplyConstant(1.0/second), children[0]); // Replace a divide with a multiply
}
if (children[0].getOperation().getId() == Operation::NEGATE && children[1].getOperation().getId() == Operation::NEGATE) // The two negations cancel
return ExpressionTreeNode(new Operation::Divide(), children[0].getChildren()[0], children[1].getChildren()[0]);
if (children[1].getOperation().getId() == Operation::NEGATE && children[0].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Negate the constant
return ExpressionTreeNode(new Operation::Divide(), ExpressionTreeNode(new Operation::MultiplyConstant(-dynamic_cast<const Operation::MultiplyConstant*>(&children[0].getOperation())->getValue()), children[0].getChildren()[0]), children[1].getChildren()[0]);
if (children[0].getOperation().getId() == Operation::NEGATE) // Pull the negation out so it can possibly be optimized further
return ExpressionTreeNode(new Operation::Negate(), ExpressionTreeNode(new Operation::Divide(), children[0].getChildren()[0], children[1]));
if (children[1].getOperation().getId() == Operation::NEGATE) // Pull the negation out so it can possibly be optimized further
return ExpressionTreeNode(new Operation::Negate(), ExpressionTreeNode(new Operation::Divide(), children[0], children[1].getChildren()[0]));
if (children[1].getOperation().getId() == Operation::RECIPROCAL) // a/(1/b) = a*b
return ExpressionTreeNode(new Operation::Multiply(), children[0], children[1].getChildren()[0]);
break;
}
case Operation::POWER:
{
if (first_const && first == 0.0) // 0 to any power is 0
return ExpressionTreeNode(new Operation::Constant(0.0));
if (first_const && first == 1.0) // 1 to any power is 1
return ExpressionTreeNode(new Operation::Constant(1.0));
if (second_const) { // Constant exponent
if (second == 0.0) // x^0 = 1
return ExpressionTreeNode(new Operation::Constant(1.0));
if (second == 1.0) // x^1 = x
return children[0];
if (second == -1.0) // x^-1 = recip(x)
return ExpressionTreeNode(new Operation::Reciprocal(), children[0]);
if (second == 2.0) // x^2 = square(x)
return ExpressionTreeNode(new Operation::Square(), children[0]);
if (second == 3.0) // x^3 = cube(x)
return ExpressionTreeNode(new Operation::Cube(), children[0]);
if (second == 0.5) // x^0.5 = sqrt(x)
return ExpressionTreeNode(new Operation::Sqrt(), children[0]);
// Constant power
return ExpressionTreeNode(new Operation::PowerConstant(second), children[0]);
}
break;
}
case Operation::NEGATE:
{
if (children[0].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Combine a multiply and a negate into a single multiply
return ExpressionTreeNode(new Operation::MultiplyConstant(-dynamic_cast<const Operation::MultiplyConstant*>(&children[0].getOperation())->getValue()), children[0].getChildren()[0]);
if (first_const) // Negate a constant
return ExpressionTreeNode(new Operation::Constant(-first));
if (children[0].getOperation().getId() == Operation::NEGATE) // The two negations cancel
return children[0].getChildren()[0];
break;
}
case Operation::MULTIPLY_CONSTANT:
{
if (children[0].getOperation().getId() == Operation::MULTIPLY_CONSTANT) // Combine two multiplies into a single one
return ExpressionTreeNode(new Operation::MultiplyConstant(dynamic_cast<const Operation::MultiplyConstant*>(&node.getOperation())->getValue()*dynamic_cast<const Operation::MultiplyConstant*>(&children[0].getOperation())->getValue()), children[0].getChildren()[0]);
if (first_const) // Multiply two constants
return ExpressionTreeNode(new Operation::Constant(dynamic_cast<const Operation::MultiplyConstant*>(&node.getOperation())->getValue()*getConstantValue(children[0])));
if (children[0].getOperation().getId() == Operation::NEGATE) // Combine a multiply and a negate into a single multiply
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,
// we don't substitute a simpler expression.
break;
}
}
return ExpressionTreeNode(node.getOperation().clone(), children);
}
ParsedExpression ParsedExpression::differentiate(const string& variable) const {
return differentiate(getRootNode(), variable);
}
ExpressionTreeNode ParsedExpression::differentiate(const ExpressionTreeNode& node, const string& variable) {
vector<ExpressionTreeNode> childDerivs(node.getChildren().size());
for (int i = 0; i < (int) childDerivs.size(); i++)
childDerivs[i] = differentiate(node.getChildren()[i], variable);
return node.getOperation().differentiate(node.getChildren(),childDerivs, variable);
}
bool ParsedExpression::isConstant(const ExpressionTreeNode& node) {
return (node.getOperation().getId() == Operation::CONSTANT);
}
double ParsedExpression::getConstantValue(const ExpressionTreeNode& node) {
if (node.getOperation().getId() != Operation::CONSTANT) {
throw Exception("getConstantValue called on a non-constant ExpressionNode");
}
return dynamic_cast<const Operation::Constant&>(node.getOperation()).getValue();
}
ExpressionProgram ParsedExpression::createProgram() const {
return ExpressionProgram(*this);
}
CompiledExpression ParsedExpression::createCompiledExpression() const {
return CompiledExpression(*this);
}
ParsedExpression ParsedExpression::renameVariables(const map<string, string>& replacements) const {
return ParsedExpression(renameNodeVariables(getRootNode(), replacements));
}
ExpressionTreeNode ParsedExpression::renameNodeVariables(const ExpressionTreeNode& node, const map<string, string>& replacements) {
if (node.getOperation().getId() == Operation::VARIABLE) {
map<string, string>::const_iterator replace = replacements.find(node.getOperation().getName());
if (replace != replacements.end())
return ExpressionTreeNode(new Operation::Variable(replace->second));
}
vector<ExpressionTreeNode> children;
for (int i = 0; i < (int) node.getChildren().size(); i++)
children.push_back(renameNodeVariables(node.getChildren()[i], replacements));
return ExpressionTreeNode(node.getOperation().clone(), children);
}
ostream& LMP_Lepton::operator<<(ostream& out, const ExpressionTreeNode& node) {
if (node.getOperation().isInfixOperator() && node.getChildren().size() == 2) {
out << "(" << node.getChildren()[0] << ")" << node.getOperation().getName() << "(" << node.getChildren()[1] << ")";
}
else if (node.getOperation().isInfixOperator() && node.getChildren().size() == 1) {
out << "(" << node.getChildren()[0] << ")" << node.getOperation().getName();
}
else {
out << node.getOperation().getName();
if (node.getChildren().size() > 0) {
out << "(";
for (int i = 0; i < (int) node.getChildren().size(); i++) {
if (i > 0)
out << ", ";
out << node.getChildren()[i];
}
out << ")";
}
}
return out;
}
ostream& LMP_Lepton::operator<<(ostream& out, const ParsedExpression& exp) {
out << exp.getRootNode();
return out;
}

409
lib/lepton/src/Parser.cpp Normal file
View File

@ -0,0 +1,409 @@
/* -------------------------------------------------------------------------- *
* 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) 2009-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/Parser.h"
#include "lepton/CustomFunction.h"
#include "lepton/Exception.h"
#include "lepton/ExpressionTreeNode.h"
#include "lepton/Operation.h"
#include "lepton/ParsedExpression.h"
#include <cctype>
#include <iostream>
using namespace LMP_Lepton;
using namespace std;
static const string Digits = "0123456789";
static const string Operators = "+-*/^";
static const bool LeftAssociative[] = {true, true, true, true, false};
static const int Precedence[] = {0, 0, 1, 1, 3};
static const Operation::Id OperationId[] = {Operation::ADD, Operation::SUBTRACT, Operation::MULTIPLY, Operation::DIVIDE, Operation::POWER};
class LMP_Lepton::ParseToken {
public:
enum Type {Number, Operator, Variable, Function, LeftParen, RightParen, Comma, Whitespace};
ParseToken(string text, Type type) : text(text), type(type) {
}
const string& getText() const {
return text;
}
Type getType() const {
return type;
}
private:
string text;
Type type;
};
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++)
;
for (end = (int) expression.size()-1; end > start && isspace(expression[end]); end--)
;
if (start == end && isspace(expression[end]))
return "";
return expression.substr(start, end-start+1);
}
ParseToken Parser::getNextToken(const string& expression, int start) {
char c = expression[start];
if (c == '(')
return ParseToken("(", ParseToken::LeftParen);
if (c == ')')
return ParseToken(")", ParseToken::RightParen);
if (c == ',')
return ParseToken(",", ParseToken::Comma);
if (Operators.find(c) != string::npos)
return ParseToken(string(1, c), ParseToken::Operator);
if (isspace(c)) {
// White space
for (int pos = start+1; pos < (int) expression.size(); pos++) {
if (!isspace(expression[pos]))
return ParseToken(expression.substr(start, pos-start), ParseToken::Whitespace);
}
return ParseToken(expression.substr(start, string::npos), ParseToken::Whitespace);
}
if (c == '.' || Digits.find(c) != string::npos) {
// A number
bool foundDecimal = (c == '.');
bool foundExp = false;
int pos;
for (pos = start+1; pos < (int) expression.size(); pos++) {
c = expression[pos];
if (Digits.find(c) != string::npos)
continue;
if (c == '.' && !foundDecimal) {
foundDecimal = true;
continue;
}
if ((c == 'e' || c == 'E') && !foundExp) {
foundExp = true;
if (pos < (int) expression.size()-1 && (expression[pos+1] == '-' || expression[pos+1] == '+'))
pos++;
continue;
}
break;
}
return ParseToken(expression.substr(start, pos-start), ParseToken::Number);
}
// A variable, function, or left parenthesis
for (int pos = start; pos < (int) expression.size(); pos++) {
c = expression[pos];
if (c == '(')
return ParseToken(expression.substr(start, pos-start+1), ParseToken::Function);
if (Operators.find(c) != string::npos || c == ',' || c == ')' || isspace(c))
return ParseToken(expression.substr(start, pos-start), ParseToken::Variable);
}
return ParseToken(expression.substr(start, string::npos), ParseToken::Variable);
}
vector<ParseToken> Parser::tokenize(const string& expression) {
vector<ParseToken> tokens;
int pos = 0;
while (pos < (int) expression.size()) {
ParseToken token = getNextToken(expression, pos);
if (token.getType() != ParseToken::Whitespace)
tokens.push_back(token);
pos += (int) token.getText().size();
}
return tokens;
}
ParsedExpression Parser::parse(const string& expression) {
return parse(expression, map<string, CustomFunction*>());
}
ParsedExpression Parser::parse(const string& expression, const map<string, CustomFunction*>& customFunctions) {
try {
// First split the expression into subexpressions.
string primaryExpression = expression;
vector<string> subexpressions;
while (true) {
string::size_type pos = primaryExpression.find_last_of(';');
if (pos == string::npos)
break;
string sub = trim(primaryExpression.substr(pos+1));
if (sub.size() > 0)
subexpressions.push_back(sub);
primaryExpression = primaryExpression.substr(0, pos);
}
// Parse the subexpressions.
map<string, ExpressionTreeNode> subexpDefs;
for (int i = 0; i < (int) subexpressions.size(); i++) {
string::size_type equalsPos = subexpressions[i].find('=');
if (equalsPos == string::npos)
throw Exception("subexpression does not specify a name");
string name = trim(subexpressions[i].substr(0, equalsPos));
if (name.size() == 0)
throw Exception("subexpression does not specify a name");
vector<ParseToken> tokens = tokenize(subexpressions[i].substr(equalsPos+1));
int pos = 0;
subexpDefs[name] = parsePrecedence(tokens, pos, customFunctions, subexpDefs, 0);
if (pos != tokens.size())
throw Exception("unexpected text at end of subexpression: "+tokens[pos].getText());
}
// Now parse the primary expression.
vector<ParseToken> tokens = tokenize(primaryExpression);
int pos = 0;
ExpressionTreeNode result = parsePrecedence(tokens, pos, customFunctions, subexpDefs, 0);
if (pos != tokens.size())
throw Exception("unexpected text at end of expression: "+tokens[pos].getText());
return ParsedExpression(result);
}
catch (Exception& ex) {
throw Exception("Parse error in expression \""+expression+"\": "+ex.what());
}
}
ExpressionTreeNode Parser::parsePrecedence(const vector<ParseToken>& tokens, int& pos, const map<string, CustomFunction*>& customFunctions,
const map<string, ExpressionTreeNode>& subexpressionDefs, int precedence) {
if (pos == tokens.size())
throw Exception("unexpected end of expression");
// Parse the next value (number, variable, function, parenthesized expression)
ParseToken token = tokens[pos];
ExpressionTreeNode result;
if (token.getType() == ParseToken::Number) {
double value;
stringstream(token.getText()) >> value;
result = ExpressionTreeNode(new Operation::Constant(value));
pos++;
}
else if (token.getType() == ParseToken::Variable) {
map<string, ExpressionTreeNode>::const_iterator subexp = subexpressionDefs.find(token.getText());
if (subexp == subexpressionDefs.end()) {
Operation* op = new Operation::Variable(token.getText());
result = ExpressionTreeNode(op);
}
else
result = subexp->second;
pos++;
}
else if (token.getType() == ParseToken::LeftParen) {
pos++;
result = parsePrecedence(tokens, pos, customFunctions, subexpressionDefs, 0);
if (pos == tokens.size() || tokens[pos].getType() != ParseToken::RightParen)
throw Exception("unbalanced parentheses");
pos++;
}
else if (token.getType() == ParseToken::Function) {
pos++;
vector<ExpressionTreeNode> args;
bool moreArgs;
do {
args.push_back(parsePrecedence(tokens, pos, customFunctions, subexpressionDefs, 0));
moreArgs = (pos < (int) tokens.size() && tokens[pos].getType() == ParseToken::Comma);
if (moreArgs)
pos++;
} while (moreArgs);
if (pos == tokens.size() || tokens[pos].getType() != ParseToken::RightParen)
throw Exception("unbalanced parentheses");
pos++;
Operation* op = getFunctionOperation(token.getText(), customFunctions);
try {
result = ExpressionTreeNode(op, args);
}
catch (...) {
delete op;
throw;
}
}
else if (token.getType() == ParseToken::Operator && token.getText() == "-") {
pos++;
ExpressionTreeNode toNegate = parsePrecedence(tokens, pos, customFunctions, subexpressionDefs, 2);
result = ExpressionTreeNode(new Operation::Negate(), toNegate);
}
else
throw Exception("unexpected token: "+token.getText());
// Now deal with the next binary operator.
while (pos < (int) tokens.size() && tokens[pos].getType() == ParseToken::Operator) {
token = tokens[pos];
int opIndex = (int) Operators.find(token.getText());
int opPrecedence = Precedence[opIndex];
if (opPrecedence < precedence)
return result;
pos++;
ExpressionTreeNode arg = parsePrecedence(tokens, pos, customFunctions, subexpressionDefs, LeftAssociative[opIndex] ? opPrecedence+1 : opPrecedence);
Operation* op = getOperatorOperation(token.getText());
try {
result = ExpressionTreeNode(op, result, arg);
}
catch (...) {
delete op;
throw;
}
}
return result;
}
Operation* Parser::getOperatorOperation(const std::string& name) {
switch (OperationId[Operators.find(name)]) {
case Operation::ADD:
return new Operation::Add();
case Operation::SUBTRACT:
return new Operation::Subtract();
case Operation::MULTIPLY:
return new Operation::Multiply();
case Operation::DIVIDE:
return new Operation::Divide();
case Operation::POWER:
return new Operation::Power();
default:
throw Exception("unknown operator");
}
}
Operation* Parser::getFunctionOperation(const std::string& name, const map<string, CustomFunction*>& customFunctions) {
static map<string, Operation::Id> opMap;
if (opMap.size() == 0) {
opMap["sqrt"] = Operation::SQRT;
opMap["exp"] = Operation::EXP;
opMap["log"] = Operation::LOG;
opMap["sin"] = Operation::SIN;
opMap["cos"] = Operation::COS;
opMap["sec"] = Operation::SEC;
opMap["csc"] = Operation::CSC;
opMap["tan"] = Operation::TAN;
opMap["cot"] = Operation::COT;
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;
opMap["erf"] = Operation::ERF;
opMap["erfc"] = Operation::ERFC;
opMap["step"] = Operation::STEP;
opMap["delta"] = Operation::DELTA;
opMap["square"] = Operation::SQUARE;
opMap["cube"] = Operation::CUBE;
opMap["recip"] = Operation::RECIPROCAL;
opMap["min"] = Operation::MIN;
opMap["max"] = Operation::MAX;
opMap["abs"] = Operation::ABS;
opMap["floor"] = Operation::FLOOR;
opMap["ceil"] = Operation::CEIL;
opMap["select"] = Operation::SELECT;
}
string trimmed = name.substr(0, name.size()-1);
// First check custom functions.
map<string, CustomFunction*>::const_iterator custom = customFunctions.find(trimmed);
if (custom != customFunctions.end())
return new Operation::Custom(trimmed, custom->second->clone());
// Now try standard functions.
map<string, Operation::Id>::const_iterator iter = opMap.find(trimmed);
if (iter == opMap.end())
throw Exception("unknown function: "+trimmed);
switch (iter->second) {
case Operation::SQRT:
return new Operation::Sqrt();
case Operation::EXP:
return new Operation::Exp();
case Operation::LOG:
return new Operation::Log();
case Operation::SIN:
return new Operation::Sin();
case Operation::COS:
return new Operation::Cos();
case Operation::SEC:
return new Operation::Sec();
case Operation::CSC:
return new Operation::Csc();
case Operation::TAN:
return new Operation::Tan();
case Operation::COT:
return new Operation::Cot();
case Operation::ASIN:
return new Operation::Asin();
case Operation::ACOS:
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:
return new Operation::Cosh();
case Operation::TANH:
return new Operation::Tanh();
case Operation::ERF:
return new Operation::Erf();
case Operation::ERFC:
return new Operation::Erfc();
case Operation::STEP:
return new Operation::Step();
case Operation::DELTA:
return new Operation::Delta();
case Operation::SQUARE:
return new Operation::Square();
case Operation::CUBE:
return new Operation::Cube();
case Operation::RECIPROCAL:
return new Operation::Reciprocal();
case Operation::MIN:
return new Operation::Min();
case Operation::MAX:
return new Operation::Max();
case Operation::ABS:
return new Operation::Abs();
case Operation::FLOOR:
return new Operation::Floor();
case Operation::CEIL:
return new Operation::Ceil();
case Operation::SELECT:
return new Operation::Select();
default:
throw Exception("unknown function");
}
}