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

Merge pull request #686 from giacomofiorin/lepton-library

Browse Source 
Add Lepton library within lib/colvars
This commit is contained in:
Steve Plimpton
2017-10-13 09:48:42 -06:00
committed by GitHub
parent d3ef4bd594 d90aad887e
commit 18ca2292c2
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");
}
}