tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add Lepton library within lib/colvars

Browse Source
This commit is contained in:
Giacomo Fiorin
2017-10-09 14:42:47 -04:00
parent 88a33edb50
commit d90aad887e
21 changed files with 3941 additions and 30 deletions
Download Patch File Download Diff File Expand all files Collapse all files

25
lib/colvars/Makefile.common
View File

@ -12,7 +12,7 @@
ifeq ($(COLVARS_DEBUG),)
COLVARS_DEBUG_INCFLAGS =
else
COLVARS_DEBUG_INCFLAGS= -DCOLVARS_DEBUG
COLVARS_DEBUG_INCFLAGS = -DCOLVARS_DEBUG
endif
COLVARS_INCFLAGS = $(COLVARS_DEBUG_INCFLAGS) $(COLVARS_PYTHON_INCFLAGS)
@ -21,6 +21,7 @@ COLVARS_INCFLAGS = $(COLVARS_DEBUG_INCFLAGS) $(COLVARS_PYTHON_INCFLAGS)
.SUFFIXES:
.SUFFIXES: .cpp .o
COLVARS_SRCS = \
colvaratoms.cpp \
colvarbias_abf.cpp \
@ -45,21 +46,33 @@ COLVARS_SRCS = \
colvartypes.cpp \
colvarvalue.cpp
COLVARS_OBJS = $(COLVARS_SRCS:.cpp=.o)
LEPTON_SRCS = \
lepton/src/CompiledExpression.cpp lepton/src/ExpressionTreeNode.cpp \
lepton/src/ParsedExpression.cpp lepton/src/ExpressionProgram.cpp \
lepton/src/Operation.cpp lepton/src/Parser.cpp
.cpp.o:
$(CXX) $(CXXFLAGS) $(COLVARS_INCFLAGS) -c $<
LEPTON_OBJS = \
lepton/src/CompiledExpression.o lepton/src/ExpressionTreeNode.o \
lepton/src/ParsedExpression.o lepton/src/ExpressionProgram.o \
lepton/src/Operation.o lepton/src/Parser.o
COLVARS_OBJS = $(COLVARS_SRCS:.cpp=.o) $(LEPTON_OBJS)
%.o: %.cpp
$(CXX) $(CXXFLAGS) $(COLVARS_INCFLAGS) \
-Ilepton/include -DLEPTON -c -o $@ $<
$(COLVARS_LIB): Makefile.deps $(COLVARS_OBJS)
$(AR) $(ARFLAGS) $(COLVARS_LIB) $(COLVARS_OBJS)
$(AR) $(ARFLAGS) $(COLVARS_LIB) $(COLVARS_OBJS) $(LEPTON_OBJS)
Makefile.deps: $(COLVARS_SRCS)
@echo > $@
@for src in $^ ; do \
obj=`basename $$src .cpp`.o ; \
$(CXX) -MM $(COLVARS_INCFLAGS) \
$(CXX) -MM $(COLVARS_INCFLAGS) -Ilepton/include -DLEPTON \
-MT '$$(COLVARS_OBJ_DIR)'$$obj $$src >> $@ ; \
done
include Makefile.deps
include Makefile.lepton.deps # Hand-generated

206
lib/colvars/Makefile.deps
View File

@ -4,73 +4,231 @@ $(COLVARS_OBJ_DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h \
colvarparse.h colvaratoms.h colvardeps.h
$(COLVARS_OBJ_DIR)colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h colvar.h \
colvarparse.h colvardeps.h colvarbias_abf.h colvarbias.h colvargrid.h
colvarparse.h colvardeps.h lepton/include/Lepton.h \
lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarbias_abf.h colvarbias.h colvargrid.h
$(COLVARS_OBJ_DIR)colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarbias_alb.h colvar.h colvarparse.h colvardeps.h colvarbias.h
colvarbias_alb.h colvar.h colvarparse.h colvardeps.h \
lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarbias.h
$(COLVARS_OBJ_DIR)colvarbias.o: colvarbias.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h \
colvar.h colvarparse.h colvardeps.h
colvar.h colvarparse.h colvardeps.h lepton/include/Lepton.h \
lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h
$(COLVARS_OBJ_DIR)colvarbias_histogram.o: colvarbias_histogram.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvar.h colvarparse.h colvardeps.h colvarbias_histogram.h \
colvarbias.h colvargrid.h
colvarvalue.h colvar.h colvarparse.h colvardeps.h \
lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarbias_histogram.h colvarbias.h colvargrid.h
$(COLVARS_OBJ_DIR)colvarbias_meta.o: colvarbias_meta.cpp colvar.h \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvardeps.h colvarbias_meta.h colvarbias.h \
colvargrid.h
colvarvalue.h colvarparse.h colvardeps.h lepton/include/Lepton.h \
lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarbias_meta.h colvarbias.h colvargrid.h
$(COLVARS_OBJ_DIR)colvarbias_restraint.o: colvarbias_restraint.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarbias_restraint.h colvarbias.h colvar.h colvarparse.h \
colvardeps.h
colvardeps.h lepton/include/Lepton.h \
lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h
$(COLVARS_OBJ_DIR)colvarcomp_angles.o: colvarcomp_angles.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvar.h colvarparse.h colvardeps.h colvarcomp.h \
colvaratoms.h
colvarvalue.h colvar.h colvarparse.h colvardeps.h \
lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarcomp.h colvaratoms.h
$(COLVARS_OBJ_DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvaratoms.h colvardeps.h colvar.h \
lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarcomp.h
$(COLVARS_OBJ_DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h colvar.h \
colvarparse.h colvardeps.h colvarcomp.h colvaratoms.h
colvarparse.h colvardeps.h lepton/include/Lepton.h \
lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarcomp.h colvaratoms.h
$(COLVARS_OBJ_DIR)colvarcomp_distances.o: colvarcomp_distances.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvar.h colvardeps.h colvarcomp.h \
colvaratoms.h
colvarvalue.h colvarparse.h colvar.h colvardeps.h \
lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarcomp.h colvaratoms.h
$(COLVARS_OBJ_DIR)colvarcomp_protein.o: colvarcomp_protein.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvar.h colvardeps.h colvarcomp.h \
colvaratoms.h
colvarvalue.h colvarparse.h colvar.h colvardeps.h \
lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarcomp.h colvaratoms.h
$(COLVARS_OBJ_DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h colvar.h colvardeps.h colvarcomp.h \
colvaratoms.h
colvarvalue.h colvarparse.h colvar.h colvardeps.h \
lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarcomp.h colvaratoms.h
$(COLVARS_OBJ_DIR)colvar.o: colvar.cpp colvarmodule.h colvars_version.h \
colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \
colvardeps.h colvarcomp.h colvaratoms.h colvarscript.h colvarbias.h
colvardeps.h lepton/include/Lepton.h \
lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarcomp.h colvaratoms.h colvarscript.h colvarbias.h
$(COLVARS_OBJ_DIR)colvardeps.o: colvardeps.cpp colvardeps.h \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarparse.h
$(COLVARS_OBJ_DIR)colvargrid.o: colvargrid.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarparse.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \
colvargrid.h
colvarparse.h colvar.h colvardeps.h lepton/include/Lepton.h \
lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarcomp.h colvaratoms.h colvargrid.h
$(COLVARS_OBJ_DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarparse.h colvar.h colvardeps.h colvarbias.h colvarbias_abf.h \
colvargrid.h colvarbias_alb.h colvarbias_histogram.h colvarbias_meta.h \
colvarbias_restraint.h colvarscript.h colvaratoms.h
colvarparse.h colvar.h colvardeps.h lepton/include/Lepton.h \
lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvarbias.h colvarbias_abf.h colvargrid.h colvarbias_alb.h \
colvarbias_histogram.h colvarbias_meta.h colvarbias_restraint.h \
colvarscript.h colvaratoms.h
$(COLVARS_OBJ_DIR)colvarparse.o: colvarparse.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarparse.h
$(COLVARS_OBJ_DIR)colvarproxy.o: colvarproxy.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarscript.h colvarbias.h colvar.h colvarparse.h colvardeps.h \
lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
colvaratoms.h
$(COLVARS_OBJ_DIR)colvarscript.o: colvarscript.cpp colvarscript.h \
colvarmodule.h colvars_version.h colvartypes.h colvarproxy.h \
colvarvalue.h colvarbias.h colvar.h colvarparse.h colvardeps.h
colvarvalue.h colvarbias.h colvar.h colvarparse.h colvardeps.h \
lepton/include/Lepton.h lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h
$(COLVARS_OBJ_DIR)colvartypes.o: colvartypes.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarproxy.h colvarvalue.h \
colvarparse.h

40
lib/colvars/Makefile.lepton.deps Normal file
View File

@ -0,0 +1,40 @@
lepton/src/CompiledExpression.o: lepton/src/CompiledExpression.cpp \
lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h
lepton/src/ExpressionProgram.o: lepton/src/ExpressionProgram.cpp \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h
lepton/src/ExpressionTreeNode.o: lepton/src/ExpressionTreeNode.cpp \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/Exception.h lepton/include/lepton/Operation.h \
lepton/include/lepton/CustomFunction.h lepton/include/lepton/Exception.h
lepton/src/Operation.o: lepton/src/Operation.cpp \
lepton/include/lepton/Operation.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h lepton/include/lepton/Exception.h \
lepton/include/lepton/ExpressionTreeNode.h lepton/src/MSVC_erfc.h
lepton/src/ParsedExpression.o: lepton/src/ParsedExpression.cpp \
lepton/include/lepton/ParsedExpression.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CompiledExpression.h \
lepton/include/lepton/ExpressionProgram.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h
lepton/src/Parser.o: lepton/src/Parser.cpp \
lepton/include/lepton/Parser.h lepton/include/lepton/windowsIncludes.h \
lepton/include/lepton/CustomFunction.h lepton/include/lepton/Exception.h \
lepton/include/lepton/ExpressionTreeNode.h \
lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
lepton/include/lepton/Exception.h \
lepton/include/lepton/ParsedExpression.h \
lepton/include/lepton/ExpressionTreeNode.h

4
lib/colvars/README
View File

@ -47,6 +47,10 @@ correct for your system, else the LAMMPS build will likely fail.
If you want to set a debug flag recognized by the library, the
settings in Makefile.common should work.
Note: some Colvars functions use the Lepton mathematical expression parser,
which is here included (no additional steps required). For more details, see:
https://simtk.org/projects/lepton
## Documentation

43
lib/colvars/lepton/include/Lepton.h Normal file
View File

@ -0,0 +1,43 @@
#ifndef LEPTON_H_
#define 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 /*LEPTON_H_*/

113
lib/colvars/lepton/include/lepton/CompiledExpression.h Normal file
View File

@ -0,0 +1,113 @@
#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-2016 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 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;
void* jitCode;
#ifdef LEPTON_USE_JIT
void generateJitCode();
void generateSingleArgCall(asmjit::X86Compiler& c, asmjit::X86XmmVar& dest, asmjit::X86XmmVar& arg, double (*function)(double));
std::vector<double> constants;
asmjit::JitRuntime runtime;
#endif
};
} // namespace Lepton
#endif /*LEPTON_COMPILED_EXPRESSION_H_*/

77
lib/colvars/lepton/include/lepton/CustomFunction.h Normal file
View File

@ -0,0 +1,77 @@
#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 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;
};
} // namespace Lepton
#endif /*LEPTON_CUSTOM_FUNCTION_H_*/

59
lib/colvars/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 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 Lepton
#endif /*LEPTON_EXCEPTION_H_*/

95
lib/colvars/lepton/include/lepton/ExpressionProgram.h Normal file
View File

@ -0,0 +1,95 @@
#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 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 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;
/**
* 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 Lepton
#endif /*LEPTON_EXPRESSION_PROGRAM_H_*/

105
lib/colvars/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 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 Lepton
#endif /*LEPTON_EXPRESSION_TREE_NODE_H_*/

1165
lib/colvars/lepton/include/lepton/Operation.h Normal file
View File

File diff suppressed because it is too large Load Diff

130
lib/colvars/lepton/include/lepton/ParsedExpression.h Normal file
View File

@ -0,0 +1,130 @@
#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 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 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 Lepton
#endif /*LEPTON_PARSED_EXPRESSION_H_*/

77
lib/colvars/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 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 Lepton
#endif /*LEPTON_PARSER_H_*/

41
lib/colvars/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_

400
lib/colvars/lepton/src/CompiledExpression.cpp Normal file
View File

@ -0,0 +1,400 @@
/* -------------------------------------------------------------------------- *
* 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-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "lepton/CompiledExpression.h"
#include "lepton/Operation.h"
#include "lepton/ParsedExpression.h"
#include <utility>
using namespace Lepton;
using namespace std;
#ifdef LEPTON_USE_JIT
using namespace asmjit;
#endif
CompiledExpression::CompiledExpression() : jitCode(NULL) {
}
CompiledExpression::CompiledExpression(const ParsedExpression& expression) : jitCode(NULL) {
ParsedExpression expr = expression.optimize(); // Just in case it wasn't already optimized.
vector<pair<ExpressionTreeNode, int> > temps;
compileExpression(expr.getRootNode(), temps);
int maxArguments = 1;
for (int i = 0; i < (int) operation.size(); i++)
if (operation[i]->getNumArguments() > maxArguments)
maxArguments = operation[i]->getNumArguments();
argValues.resize(maxArguments);
#ifdef LEPTON_USE_JIT
generateJitCode();
#endif
}
CompiledExpression::~CompiledExpression() {
for (int i = 0; i < (int) operation.size(); i++)
if (operation[i] != NULL)
delete operation[i];
}
CompiledExpression::CompiledExpression(const CompiledExpression& expression) : jitCode(NULL) {
*this = expression;
}
CompiledExpression& CompiledExpression::operator=(const CompiledExpression& expression) {
arguments = expression.arguments;
target = expression.target;
variableIndices = expression.variableIndices;
variableNames = expression.variableNames;
workspace.resize(expression.workspace.size());
argValues.resize(expression.argValues.size());
operation.resize(expression.operation.size());
for (int i = 0; i < (int) operation.size(); i++)
operation[i] = expression.operation[i]->clone();
setVariableLocations(variablePointers);
return *this;
}
void CompiledExpression::compileExpression(const ExpressionTreeNode& node, vector<pair<ExpressionTreeNode, int> >& temps) {
if (findTempIndex(node, temps) != -1)
return; // We have already processed a node identical to this one.
// Process the child nodes.
vector<int> args;
for (int i = 0; i < node.getChildren().size(); i++) {
compileExpression(node.getChildren()[i], temps);
args.push_back(findTempIndex(node.getChildren()[i], temps));
}
// Process this node.
if (node.getOperation().getId() == Operation::VARIABLE) {
variableIndices[node.getOperation().getName()] = (int) workspace.size();
variableNames.insert(node.getOperation().getName());
}
else {
int stepIndex = (int) arguments.size();
arguments.push_back(vector<int>());
target.push_back((int) workspace.size());
operation.push_back(node.getOperation().clone());
if (args.size() == 0)
arguments[stepIndex].push_back(0); // The value won't actually be used. We just need something there.
else {
// If the arguments are sequential, we can just pass a pointer to the first one.
bool sequential = true;
for (int i = 1; i < args.size(); i++)
if (args[i] != args[i-1]+1)
sequential = false;
if (sequential)
arguments[stepIndex].push_back(args[0]);
else
arguments[stepIndex] = args;
}
}
temps.push_back(make_pair(node, (int) workspace.size()));
workspace.push_back(0.0);
}
int CompiledExpression::findTempIndex(const ExpressionTreeNode& node, vector<pair<ExpressionTreeNode, int> >& temps) {
for (int i = 0; i < (int) temps.size(); i++)
if (temps[i].first == node)
return i;
return -1;
}
const set<string>& CompiledExpression::getVariables() const {
return variableNames;
}
double& CompiledExpression::getVariableReference(const string& name) {
map<string, double*>::iterator pointer = variablePointers.find(name);
if (pointer != variablePointers.end())
return *pointer->second;
map<string, int>::iterator index = variableIndices.find(name);
if (index == variableIndices.end())
throw Exception("getVariableReference: Unknown variable '"+name+"'");
return workspace[index->second];
}
void CompiledExpression::setVariableLocations(map<string, double*>& variableLocations) {
variablePointers = variableLocations;
#ifdef LEPTON_USE_JIT
// Rebuild the JIT code.
if (workspace.size() > 0)
generateJitCode();
#else
// Make a list of all variables we will need to copy before evaluating the expression.
variablesToCopy.clear();
for (map<string, int>::const_iterator iter = variableIndices.begin(); iter != variableIndices.end(); ++iter) {
map<string, double*>::iterator pointer = variablePointers.find(iter->first);
if (pointer != variablePointers.end())
variablesToCopy.push_back(make_pair(&workspace[iter->second], pointer->second));
}
#endif
}
double CompiledExpression::evaluate() const {
#ifdef LEPTON_USE_JIT
return ((double (*)()) 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) {
map<string, double>* dummyVariables = NULL;
return op->evaluate(args, *dummyVariables);
}
void CompiledExpression::generateJitCode() {
X86Compiler c(&runtime);
c.addFunc(kFuncConvHost, FuncBuilder0<double>());
vector<X86XmmVar> workspaceVar(workspace.size());
for (int i = 0; i < (int) workspaceVar.size(); i++)
workspaceVar[i] = c.newXmmVar(kX86VarTypeXmmSd);
X86GpVar argsPointer(c);
c.mov(argsPointer, imm_ptr(&argValues[0]));
// Load the arguments into variables.
for (set<string>::const_iterator iter = variableNames.begin(); iter != variableNames.end(); ++iter) {
map<string, int>::iterator index = variableIndices.find(*iter);
X86GpVar variablePointer(c);
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<X86XmmVar> constantVar(constants.size());
if (constants.size() > 0) {
X86GpVar constantsPointer(c);
c.mov(constantsPointer, imm_ptr(&constants[0]));
for (int i = 0; i < (int) constants.size(); i++) {
constantVar[i] = c.newXmmVar(kX86VarTypeXmmSd);
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::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::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]]);
X86GpVar fn(c, kVarTypeIntPtr);
c.mov(fn, imm_ptr((void*) evaluateOperation));
X86CallNode* call = c.call(fn, kFuncConvHost, FuncBuilder2<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();
jitCode = c.make();
}
void CompiledExpression::generateSingleArgCall(X86Compiler& c, X86XmmVar& dest, X86XmmVar& arg, double (*function)(double)) {
X86GpVar fn(c, kVarTypeIntPtr);
c.mov(fn, imm_ptr((void*) function));
X86CallNode* call = c.call(fn, kFuncConvHost, FuncBuilder1<double, double>());
call->setArg(0, arg);
call->setRet(0, dest);
}
#endif

105
lib/colvars/lepton/src/ExpressionProgram.cpp Normal file
View File

@ -0,0 +1,105 @@
/* -------------------------------------------------------------------------- *
* 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 "lepton/ExpressionProgram.h"
#include "lepton/Operation.h"
#include "lepton/ParsedExpression.h"
using namespace 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];
}
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/colvars/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 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;
}

87
lib/colvars/lepton/src/MSVC_erfc.h Normal file
View File

@ -0,0 +1,87 @@
#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)
#define M_PI 3.14159265358979323846264338327950288
#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_

335
lib/colvars/lepton/src/Operation.cpp Normal file
View File

@ -0,0 +1,335 @@
/* -------------------------------------------------------------------------- *
* 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/Operation.h"
#include "lepton/ExpressionTreeNode.h"
#include "MSVC_erfc.h"
using namespace 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::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);
}

351
lib/colvars/lepton/src/ParsedExpression.cpp Normal file
View File

@ -0,0 +1,351 @@
/* -------------------------------------------------------------------------- *
* 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 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]);
switch (node.getOperation().getId()) {
case Operation::ADD:
{
double first = getConstantValue(children[0]);
double second = getConstantValue(children[1]);
if (first == 0.0) // Add 0
return children[1];
if (second == 0.0) // Add 0
return children[0];
if (first == first) // Add a constant
return ExpressionTreeNode(new Operation::AddConstant(first), children[1]);
if (second == second) // 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
double first = getConstantValue(children[0]);
if (first == 0.0) // Subtract from 0
return ExpressionTreeNode(new Operation::Negate(), children[1]);
double second = getConstantValue(children[1]);
if (second == 0.0) // Subtract 0
return children[0];
if (second == second) // 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:
{
double first = getConstantValue(children[0]);
double second = getConstantValue(children[1]);
if (first == 0.0 || second == 0.0) // Multiply by 0
return ExpressionTreeNode(new Operation::Constant(0.0));
if (first == 1.0) // Multiply by 1
return children[1];
if (second == 1.0) // Multiply by 1
return children[0];
if (children[0].getOperation().getId() == Operation::CONSTANT) { // 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 (children[1].getOperation().getId() == Operation::CONSTANT) { // 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
double numerator = getConstantValue(children[0]);
if (numerator == 0.0) // 0 divided by something
return ExpressionTreeNode(new Operation::Constant(0.0));
if (numerator == 1.0) // 1 divided by something
return ExpressionTreeNode(new Operation::Reciprocal(), children[1]);
double denominator = getConstantValue(children[1]);
if (denominator == 1.0) // Divide by 1
return children[0];
if (children[1].getOperation().getId() == Operation::CONSTANT) {
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()/denominator), children[0].getChildren()[0]);
return ExpressionTreeNode(new Operation::MultiplyConstant(1.0/denominator), 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:
{
double base = getConstantValue(children[0]);
if (base == 0.0) // 0 to any power is 0
return ExpressionTreeNode(new Operation::Constant(0.0));
if (base == 1.0) // 1 to any power is 1
return ExpressionTreeNode(new Operation::Constant(1.0));
double exponent = getConstantValue(children[1]);
if (exponent == 0.0) // x^0 = 1
return ExpressionTreeNode(new Operation::Constant(1.0));
if (exponent == 1.0) // x^1 = x
return children[0];
if (exponent == -1.0) // x^-1 = recip(x)
return ExpressionTreeNode(new Operation::Reciprocal(), children[0]);
if (exponent == 2.0) // x^2 = square(x)
return ExpressionTreeNode(new Operation::Square(), children[0]);
if (exponent == 3.0) // x^3 = cube(x)
return ExpressionTreeNode(new Operation::Cube(), children[0]);
if (exponent == 0.5) // x^0.5 = sqrt(x)
return ExpressionTreeNode(new Operation::Sqrt(), children[0]);
if (exponent == exponent) // Constant power
return ExpressionTreeNode(new Operation::PowerConstant(exponent), 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 (children[0].getOperation().getId() == Operation::CONSTANT) // Negate a constant
return ExpressionTreeNode(new Operation::Constant(-getConstantValue(children[0])));
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 (children[0].getOperation().getId() == Operation::CONSTANT) // 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;
}
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);
}
double ParsedExpression::getConstantValue(const ExpressionTreeNode& node) {
if (node.getOperation().getId() == Operation::CONSTANT)
return dynamic_cast<const Operation::Constant&>(node.getOperation()).getValue();
return numeric_limits<double>::quiet_NaN();
}
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& 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& Lepton::operator<<(ostream& out, const ParsedExpression& exp) {
out << exp.getRootNode();
return out;
}

406
lib/colvars/lepton/src/Parser.cpp Normal file
View File

@ -0,0 +1,406 @@
/* -------------------------------------------------------------------------- *
* 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/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 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 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["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::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");
}
}