From d90aad887eb6c5885100a6c79147e5ab8defb969 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Mon, 9 Oct 2017 14:42:47 -0400 Subject: [PATCH] Add Lepton library within lib/colvars --- lib/colvars/Makefile.common | 25 +- lib/colvars/Makefile.deps | 206 ++- lib/colvars/Makefile.lepton.deps | 40 + lib/colvars/README | 4 + lib/colvars/lepton/include/Lepton.h | 43 + .../include/lepton/CompiledExpression.h | 113 ++ .../lepton/include/lepton/CustomFunction.h | 77 ++ lib/colvars/lepton/include/lepton/Exception.h | 59 + .../lepton/include/lepton/ExpressionProgram.h | 95 ++ .../include/lepton/ExpressionTreeNode.h | 105 ++ lib/colvars/lepton/include/lepton/Operation.h | 1165 +++++++++++++++++ .../lepton/include/lepton/ParsedExpression.h | 130 ++ lib/colvars/lepton/include/lepton/Parser.h | 77 ++ .../lepton/include/lepton/windowsIncludes.h | 41 + lib/colvars/lepton/src/CompiledExpression.cpp | 400 ++++++ lib/colvars/lepton/src/ExpressionProgram.cpp | 105 ++ lib/colvars/lepton/src/ExpressionTreeNode.cpp | 107 ++ lib/colvars/lepton/src/MSVC_erfc.h | 87 ++ lib/colvars/lepton/src/Operation.cpp | 335 +++++ lib/colvars/lepton/src/ParsedExpression.cpp | 351 +++++ lib/colvars/lepton/src/Parser.cpp | 406 ++++++ 21 files changed, 3941 insertions(+), 30 deletions(-) create mode 100644 lib/colvars/Makefile.lepton.deps create mode 100644 lib/colvars/lepton/include/Lepton.h create mode 100644 lib/colvars/lepton/include/lepton/CompiledExpression.h create mode 100644 lib/colvars/lepton/include/lepton/CustomFunction.h create mode 100644 lib/colvars/lepton/include/lepton/Exception.h create mode 100644 lib/colvars/lepton/include/lepton/ExpressionProgram.h create mode 100644 lib/colvars/lepton/include/lepton/ExpressionTreeNode.h create mode 100644 lib/colvars/lepton/include/lepton/Operation.h create mode 100644 lib/colvars/lepton/include/lepton/ParsedExpression.h create mode 100644 lib/colvars/lepton/include/lepton/Parser.h create mode 100644 lib/colvars/lepton/include/lepton/windowsIncludes.h create mode 100644 lib/colvars/lepton/src/CompiledExpression.cpp create mode 100644 lib/colvars/lepton/src/ExpressionProgram.cpp create mode 100644 lib/colvars/lepton/src/ExpressionTreeNode.cpp create mode 100644 lib/colvars/lepton/src/MSVC_erfc.h create mode 100644 lib/colvars/lepton/src/Operation.cpp create mode 100644 lib/colvars/lepton/src/ParsedExpression.cpp create mode 100644 lib/colvars/lepton/src/Parser.cpp diff --git a/lib/colvars/Makefile.common b/lib/colvars/Makefile.common index f47403f771..5e2f0f268e 100644 --- a/lib/colvars/Makefile.common +++ b/lib/colvars/Makefile.common @@ -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 diff --git a/lib/colvars/Makefile.deps b/lib/colvars/Makefile.deps index f463da5f86..fec7b2be12 100644 --- a/lib/colvars/Makefile.deps +++ b/lib/colvars/Makefile.deps @@ -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 diff --git a/lib/colvars/Makefile.lepton.deps b/lib/colvars/Makefile.lepton.deps new file mode 100644 index 0000000000..93c3912384 --- /dev/null +++ b/lib/colvars/Makefile.lepton.deps @@ -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 diff --git a/lib/colvars/README b/lib/colvars/README index 5df9612dfa..087528748b 100644 --- a/lib/colvars/README +++ b/lib/colvars/README @@ -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 diff --git a/lib/colvars/lepton/include/Lepton.h b/lib/colvars/lepton/include/Lepton.h new file mode 100644 index 0000000000..22edcb3ff9 --- /dev/null +++ b/lib/colvars/lepton/include/Lepton.h @@ -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_*/ diff --git a/lib/colvars/lepton/include/lepton/CompiledExpression.h b/lib/colvars/lepton/include/lepton/CompiledExpression.h new file mode 100644 index 0000000000..67442e0cf5 --- /dev/null +++ b/lib/colvars/lepton/include/lepton/CompiledExpression.h @@ -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 +#include +#include +#include +#include +#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& 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& 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 >& temps); + int findTempIndex(const ExpressionTreeNode& node, std::vector >& temps); + std::map variablePointers; + std::vector > variablesToCopy; + std::vector > arguments; + std::vector target; + std::vector operation; + std::map variableIndices; + std::set variableNames; + mutable std::vector workspace; + mutable std::vector argValues; + std::map 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 constants; + asmjit::JitRuntime runtime; +#endif +}; + +} // namespace Lepton + +#endif /*LEPTON_COMPILED_EXPRESSION_H_*/ diff --git a/lib/colvars/lepton/include/lepton/CustomFunction.h b/lib/colvars/lepton/include/lepton/CustomFunction.h new file mode 100644 index 0000000000..5c5586105f --- /dev/null +++ b/lib/colvars/lepton/include/lepton/CustomFunction.h @@ -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_*/ diff --git a/lib/colvars/lepton/include/lepton/Exception.h b/lib/colvars/lepton/include/lepton/Exception.h new file mode 100644 index 0000000000..5ad55714d1 --- /dev/null +++ b/lib/colvars/lepton/include/lepton/Exception.h @@ -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 +#include + +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_*/ diff --git a/lib/colvars/lepton/include/lepton/ExpressionProgram.h b/lib/colvars/lepton/include/lepton/ExpressionProgram.h new file mode 100644 index 0000000000..94d37f471d --- /dev/null +++ b/lib/colvars/lepton/include/lepton/ExpressionProgram.h @@ -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 +#include +#include + +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& variables) const; +private: + friend class ParsedExpression; + ExpressionProgram(const ParsedExpression& expression); + void buildProgram(const ExpressionTreeNode& node); + std::vector operations; + int maxArgs, stackSize; +}; + +} // namespace Lepton + +#endif /*LEPTON_EXPRESSION_PROGRAM_H_*/ diff --git a/lib/colvars/lepton/include/lepton/ExpressionTreeNode.h b/lib/colvars/lepton/include/lepton/ExpressionTreeNode.h new file mode 100644 index 0000000000..bf3a9a0902 --- /dev/null +++ b/lib/colvars/lepton/include/lepton/ExpressionTreeNode.h @@ -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 +#include + +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& 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& getChildren() const; +private: + Operation* operation; + std::vector children; +}; + +} // namespace Lepton + +#endif /*LEPTON_EXPRESSION_TREE_NODE_H_*/ diff --git a/lib/colvars/lepton/include/lepton/Operation.h b/lib/colvars/lepton/include/lepton/Operation.h new file mode 100644 index 0000000000..f7a8b78163 --- /dev/null +++ b/lib/colvars/lepton/include/lepton/Operation.h @@ -0,0 +1,1165 @@ +#ifndef LEPTON_OPERATION_H_ +#define LEPTON_OPERATION_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-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 "windowsIncludes.h" +#include "CustomFunction.h" +#include "Exception.h" +#include +#include +#include +#include +#include +#include + +namespace Lepton { + +class ExpressionTreeNode; + +/** + * An Operation represents a single step in the evaluation of an expression, such as a function, + * an operator, or a constant value. Each Operation takes some number of values as arguments + * and produces a single value. + * + * This is an abstract class with subclasses for specific operations. + */ + +class LEPTON_EXPORT Operation { +public: + virtual ~Operation() { + } + /** + * This enumeration lists all Operation subclasses. This is provided so that switch statements + * can be used when processing or analyzing parsed expressions. + */ + enum Id {CONSTANT, VARIABLE, CUSTOM, ADD, SUBTRACT, MULTIPLY, DIVIDE, POWER, NEGATE, SQRT, EXP, LOG, + SIN, COS, SEC, CSC, TAN, COT, ASIN, ACOS, ATAN, SINH, COSH, TANH, ERF, ERFC, STEP, DELTA, SQUARE, CUBE, RECIPROCAL, + ADD_CONSTANT, MULTIPLY_CONSTANT, POWER_CONSTANT, MIN, MAX, ABS, FLOOR, CEIL, SELECT}; + /** + * Get the name of this Operation. + */ + virtual std::string getName() const = 0; + /** + * Get this Operation's ID. + */ + virtual Id getId() const = 0; + /** + * Get the number of arguments this operation expects. + */ + virtual int getNumArguments() const = 0; + /** + * Create a clone of this Operation. + */ + virtual Operation* clone() const = 0; + /** + * Perform the computation represented by this operation. + * + * @param args the array of arguments + * @param variables a map containing the values of all variables + * @return the result of performing the computation. + */ + virtual double evaluate(double* args, const std::map& variables) const = 0; + /** + * Return an ExpressionTreeNode which represents the analytic derivative of this Operation with respect to a variable. + * + * @param children the child nodes + * @param childDerivs the derivatives of the child nodes with respect to the variable + * @param variable the variable with respect to which the derivate should be taken + */ + virtual ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const = 0; + /** + * Get whether this operation should be displayed with infix notation. + */ + virtual bool isInfixOperator() const { + return false; + } + /** + * Get whether this is a symmetric binary operation, such that exchanging its arguments + * does not affect the result. + */ + virtual bool isSymmetric() const { + return false; + } + virtual bool operator!=(const Operation& op) const { + return op.getId() != getId(); + } + virtual bool operator==(const Operation& op) const { + return !(*this != op); + } + class Constant; + class Variable; + class Custom; + class Add; + class Subtract; + class Multiply; + class Divide; + class Power; + class Negate; + class Sqrt; + class Exp; + class Log; + class Sin; + class Cos; + class Sec; + class Csc; + class Tan; + class Cot; + class Asin; + class Acos; + class Atan; + class Sinh; + class Cosh; + class Tanh; + class Erf; + class Erfc; + class Step; + class Delta; + class Square; + class Cube; + class Reciprocal; + class AddConstant; + class MultiplyConstant; + class PowerConstant; + class Min; + class Max; + class Abs; + class Floor; + class Ceil; + class Select; +}; + +class LEPTON_EXPORT Operation::Constant : public Operation { +public: + Constant(double value) : value(value) { + } + std::string getName() const { + std::stringstream name; + name << value; + return name.str(); + } + Id getId() const { + return CONSTANT; + } + int getNumArguments() const { + return 0; + } + Operation* clone() const { + return new Constant(value); + } + double evaluate(double* args, const std::map& variables) const { + return value; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; + double getValue() const { + return value; + } + bool operator!=(const Operation& op) const { + const Constant* o = dynamic_cast(&op); + return (o == NULL || o->value != value); + } +private: + double value; +}; + +class LEPTON_EXPORT Operation::Variable : public Operation { +public: + Variable(const std::string& name) : name(name) { + } + std::string getName() const { + return name; + } + Id getId() const { + return VARIABLE; + } + int getNumArguments() const { + return 0; + } + Operation* clone() const { + return new Variable(name); + } + double evaluate(double* args, const std::map& variables) const { + std::map::const_iterator iter = variables.find(name); + if (iter == variables.end()) + throw Exception("No value specified for variable "+name); + return iter->second; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; + bool operator!=(const Operation& op) const { + const Variable* o = dynamic_cast(&op); + return (o == NULL || o->name != name); + } +private: + std::string name; +}; + +class LEPTON_EXPORT Operation::Custom : public Operation { +public: + Custom(const std::string& name, CustomFunction* function) : name(name), function(function), isDerivative(false), derivOrder(function->getNumArguments(), 0) { + } + Custom(const Custom& base, int derivIndex) : name(base.name), function(base.function->clone()), isDerivative(true), derivOrder(base.derivOrder) { + derivOrder[derivIndex]++; + } + ~Custom() { + delete function; + } + std::string getName() const { + return name; + } + Id getId() const { + return CUSTOM; + } + int getNumArguments() const { + return function->getNumArguments(); + } + Operation* clone() const { + Custom* clone = new Custom(name, function->clone()); + clone->isDerivative = isDerivative; + clone->derivOrder = derivOrder; + return clone; + } + double evaluate(double* args, const std::map& variables) const { + if (isDerivative) + return function->evaluateDerivative(args, &derivOrder[0]); + return function->evaluate(args); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; + const std::vector& getDerivOrder() const { + return derivOrder; + } + bool operator!=(const Operation& op) const { + const Custom* o = dynamic_cast(&op); + return (o == NULL || o->name != name || o->isDerivative != isDerivative || o->derivOrder != derivOrder); + } +private: + std::string name; + CustomFunction* function; + bool isDerivative; + std::vector derivOrder; +}; + +class LEPTON_EXPORT Operation::Add : public Operation { +public: + Add() { + } + std::string getName() const { + return "+"; + } + Id getId() const { + return ADD; + } + int getNumArguments() const { + return 2; + } + Operation* clone() const { + return new Add(); + } + double evaluate(double* args, const std::map& variables) const { + return args[0]+args[1]; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; + bool isInfixOperator() const { + return true; + } + bool isSymmetric() const { + return true; + } +}; + +class LEPTON_EXPORT Operation::Subtract : public Operation { +public: + Subtract() { + } + std::string getName() const { + return "-"; + } + Id getId() const { + return SUBTRACT; + } + int getNumArguments() const { + return 2; + } + Operation* clone() const { + return new Subtract(); + } + double evaluate(double* args, const std::map& variables) const { + return args[0]-args[1]; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; + bool isInfixOperator() const { + return true; + } +}; + +class LEPTON_EXPORT Operation::Multiply : public Operation { +public: + Multiply() { + } + std::string getName() const { + return "*"; + } + Id getId() const { + return MULTIPLY; + } + int getNumArguments() const { + return 2; + } + Operation* clone() const { + return new Multiply(); + } + double evaluate(double* args, const std::map& variables) const { + return args[0]*args[1]; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; + bool isInfixOperator() const { + return true; + } + bool isSymmetric() const { + return true; + } +}; + +class LEPTON_EXPORT Operation::Divide : public Operation { +public: + Divide() { + } + std::string getName() const { + return "/"; + } + Id getId() const { + return DIVIDE; + } + int getNumArguments() const { + return 2; + } + Operation* clone() const { + return new Divide(); + } + double evaluate(double* args, const std::map& variables) const { + return args[0]/args[1]; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; + bool isInfixOperator() const { + return true; + } +}; + +class LEPTON_EXPORT Operation::Power : public Operation { +public: + Power() { + } + std::string getName() const { + return "^"; + } + Id getId() const { + return POWER; + } + int getNumArguments() const { + return 2; + } + Operation* clone() const { + return new Power(); + } + double evaluate(double* args, const std::map& variables) const { + return std::pow(args[0], args[1]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; + bool isInfixOperator() const { + return true; + } +}; + +class LEPTON_EXPORT Operation::Negate : public Operation { +public: + Negate() { + } + std::string getName() const { + return "-"; + } + Id getId() const { + return NEGATE; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Negate(); + } + double evaluate(double* args, const std::map& variables) const { + return -args[0]; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Sqrt : public Operation { +public: + Sqrt() { + } + std::string getName() const { + return "sqrt"; + } + Id getId() const { + return SQRT; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Sqrt(); + } + double evaluate(double* args, const std::map& variables) const { + return std::sqrt(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Exp : public Operation { +public: + Exp() { + } + std::string getName() const { + return "exp"; + } + Id getId() const { + return EXP; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Exp(); + } + double evaluate(double* args, const std::map& variables) const { + return std::exp(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Log : public Operation { +public: + Log() { + } + std::string getName() const { + return "log"; + } + Id getId() const { + return LOG; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Log(); + } + double evaluate(double* args, const std::map& variables) const { + return std::log(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Sin : public Operation { +public: + Sin() { + } + std::string getName() const { + return "sin"; + } + Id getId() const { + return SIN; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Sin(); + } + double evaluate(double* args, const std::map& variables) const { + return std::sin(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Cos : public Operation { +public: + Cos() { + } + std::string getName() const { + return "cos"; + } + Id getId() const { + return COS; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Cos(); + } + double evaluate(double* args, const std::map& variables) const { + return std::cos(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Sec : public Operation { +public: + Sec() { + } + std::string getName() const { + return "sec"; + } + Id getId() const { + return SEC; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Sec(); + } + double evaluate(double* args, const std::map& variables) const { + return 1.0/std::cos(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Csc : public Operation { +public: + Csc() { + } + std::string getName() const { + return "csc"; + } + Id getId() const { + return CSC; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Csc(); + } + double evaluate(double* args, const std::map& variables) const { + return 1.0/std::sin(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Tan : public Operation { +public: + Tan() { + } + std::string getName() const { + return "tan"; + } + Id getId() const { + return TAN; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Tan(); + } + double evaluate(double* args, const std::map& variables) const { + return std::tan(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Cot : public Operation { +public: + Cot() { + } + std::string getName() const { + return "cot"; + } + Id getId() const { + return COT; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Cot(); + } + double evaluate(double* args, const std::map& variables) const { + return 1.0/std::tan(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Asin : public Operation { +public: + Asin() { + } + std::string getName() const { + return "asin"; + } + Id getId() const { + return ASIN; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Asin(); + } + double evaluate(double* args, const std::map& variables) const { + return std::asin(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Acos : public Operation { +public: + Acos() { + } + std::string getName() const { + return "acos"; + } + Id getId() const { + return ACOS; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Acos(); + } + double evaluate(double* args, const std::map& variables) const { + return std::acos(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Atan : public Operation { +public: + Atan() { + } + std::string getName() const { + return "atan"; + } + Id getId() const { + return ATAN; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Atan(); + } + double evaluate(double* args, const std::map& variables) const { + return std::atan(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Sinh : public Operation { +public: + Sinh() { + } + std::string getName() const { + return "sinh"; + } + Id getId() const { + return SINH; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Sinh(); + } + double evaluate(double* args, const std::map& variables) const { + return std::sinh(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Cosh : public Operation { +public: + Cosh() { + } + std::string getName() const { + return "cosh"; + } + Id getId() const { + return COSH; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Cosh(); + } + double evaluate(double* args, const std::map& variables) const { + return std::cosh(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Tanh : public Operation { +public: + Tanh() { + } + std::string getName() const { + return "tanh"; + } + Id getId() const { + return TANH; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Tanh(); + } + double evaluate(double* args, const std::map& variables) const { + return std::tanh(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Erf : public Operation { +public: + Erf() { + } + std::string getName() const { + return "erf"; + } + Id getId() const { + return ERF; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Erf(); + } + double evaluate(double* args, const std::map& variables) const; + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Erfc : public Operation { +public: + Erfc() { + } + std::string getName() const { + return "erfc"; + } + Id getId() const { + return ERFC; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Erfc(); + } + double evaluate(double* args, const std::map& variables) const; + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Step : public Operation { +public: + Step() { + } + std::string getName() const { + return "step"; + } + Id getId() const { + return STEP; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Step(); + } + double evaluate(double* args, const std::map& variables) const { + return (args[0] >= 0.0 ? 1.0 : 0.0); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Delta : public Operation { +public: + Delta() { + } + std::string getName() const { + return "delta"; + } + Id getId() const { + return DELTA; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Delta(); + } + double evaluate(double* args, const std::map& variables) const { + return (args[0] == 0.0 ? 1.0 : 0.0); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Square : public Operation { +public: + Square() { + } + std::string getName() const { + return "square"; + } + Id getId() const { + return SQUARE; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Square(); + } + double evaluate(double* args, const std::map& variables) const { + return args[0]*args[0]; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Cube : public Operation { +public: + Cube() { + } + std::string getName() const { + return "cube"; + } + Id getId() const { + return CUBE; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Cube(); + } + double evaluate(double* args, const std::map& variables) const { + return args[0]*args[0]*args[0]; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Reciprocal : public Operation { +public: + Reciprocal() { + } + std::string getName() const { + return "recip"; + } + Id getId() const { + return RECIPROCAL; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Reciprocal(); + } + double evaluate(double* args, const std::map& variables) const { + return 1.0/args[0]; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::AddConstant : public Operation { +public: + AddConstant(double value) : value(value) { + } + std::string getName() const { + std::stringstream name; + name << value << "+"; + return name.str(); + } + Id getId() const { + return ADD_CONSTANT; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new AddConstant(value); + } + double evaluate(double* args, const std::map& variables) const { + return args[0]+value; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; + double getValue() const { + return value; + } + bool operator!=(const Operation& op) const { + const AddConstant* o = dynamic_cast(&op); + return (o == NULL || o->value != value); + } +private: + double value; +}; + +class LEPTON_EXPORT Operation::MultiplyConstant : public Operation { +public: + MultiplyConstant(double value) : value(value) { + } + std::string getName() const { + std::stringstream name; + name << value << "*"; + return name.str(); + } + Id getId() const { + return MULTIPLY_CONSTANT; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new MultiplyConstant(value); + } + double evaluate(double* args, const std::map& variables) const { + return args[0]*value; + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; + double getValue() const { + return value; + } + bool operator!=(const Operation& op) const { + const MultiplyConstant* o = dynamic_cast(&op); + return (o == NULL || o->value != value); + } +private: + double value; +}; + +class LEPTON_EXPORT Operation::PowerConstant : public Operation { +public: + PowerConstant(double value) : value(value) { + intValue = (int) value; + isIntPower = (intValue == value); + } + std::string getName() const { + std::stringstream name; + name << "^" << value; + return name.str(); + } + Id getId() const { + return POWER_CONSTANT; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new PowerConstant(value); + } + double evaluate(double* args, const std::map& variables) const { + if (isIntPower) { + // Integer powers can be computed much more quickly by repeated multiplication. + + int exponent = intValue; + double base = args[0]; + if (exponent < 0) { + exponent = -exponent; + base = 1.0/base; + } + double result = 1.0; + while (exponent != 0) { + if ((exponent&1) == 1) + result *= base; + base *= base; + exponent = exponent>>1; + } + return result; + } + else + return std::pow(args[0], value); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; + double getValue() const { + return value; + } + bool operator!=(const Operation& op) const { + const PowerConstant* o = dynamic_cast(&op); + return (o == NULL || o->value != value); + } + bool isInfixOperator() const { + return true; + } +private: + double value; + int intValue; + bool isIntPower; +}; + +class LEPTON_EXPORT Operation::Min : public Operation { +public: + Min() { + } + std::string getName() const { + return "min"; + } + Id getId() const { + return MIN; + } + int getNumArguments() const { + return 2; + } + Operation* clone() const { + return new Min(); + } + double evaluate(double* args, const std::map& variables) const { + // parens around (std::min) are workaround for horrible microsoft max/min macro trouble + return (std::min)(args[0], args[1]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Max : public Operation { +public: + Max() { + } + std::string getName() const { + return "max"; + } + Id getId() const { + return MAX; + } + int getNumArguments() const { + return 2; + } + Operation* clone() const { + return new Max(); + } + double evaluate(double* args, const std::map& variables) const { + // parens around (std::min) are workaround for horrible microsoft max/min macro trouble + return (std::max)(args[0], args[1]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Abs : public Operation { +public: + Abs() { + } + std::string getName() const { + return "abs"; + } + Id getId() const { + return ABS; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Abs(); + } + double evaluate(double* args, const std::map& variables) const { + return std::abs(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Floor : public Operation { +public: + + Floor() { + } + std::string getName() const { + return "floor"; + } + Id getId() const { + return FLOOR; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Floor(); + } + double evaluate(double* args, const std::map& variables) const { + return std::floor(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Ceil : public Operation { +public: + Ceil() { + } + std::string getName() const { + return "ceil"; + } + Id getId() const { + return CEIL; + } + int getNumArguments() const { + return 1; + } + Operation* clone() const { + return new Ceil(); + } + double evaluate(double* args, const std::map& variables) const { + return std::ceil(args[0]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +class LEPTON_EXPORT Operation::Select : public Operation { +public: + Select() { + } + std::string getName() const { + return "select"; + } + Id getId() const { + return SELECT; + } + int getNumArguments() const { + return 3; + } + Operation* clone() const { + return new Select(); + } + double evaluate(double* args, const std::map& variables) const { + return (args[0] != 0.0 ? args[1] : args[2]); + } + ExpressionTreeNode differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const; +}; + +} // namespace Lepton + +#endif /*LEPTON_OPERATION_H_*/ diff --git a/lib/colvars/lepton/include/lepton/ParsedExpression.h b/lib/colvars/lepton/include/lepton/ParsedExpression.h new file mode 100644 index 0000000000..d88b3d5829 --- /dev/null +++ b/lib/colvars/lepton/include/lepton/ParsedExpression.h @@ -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 +#include + +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& 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& 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& replacements) const; +private: + static double evaluate(const ExpressionTreeNode& node, const std::map& variables); + static ExpressionTreeNode preevaluateVariables(const ExpressionTreeNode& node, const std::map& 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& 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_*/ diff --git a/lib/colvars/lepton/include/lepton/Parser.h b/lib/colvars/lepton/include/lepton/Parser.h new file mode 100644 index 0000000000..63d5988d5f --- /dev/null +++ b/lib/colvars/lepton/include/lepton/Parser.h @@ -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 +#include +#include + +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& customFunctions); +private: + static std::string trim(const std::string& expression); + static std::vector tokenize(const std::string& expression); + static ParseToken getNextToken(const std::string& expression, int start); + static ExpressionTreeNode parsePrecedence(const std::vector& tokens, int& pos, const std::map& customFunctions, + const std::map& subexpressionDefs, int precedence); + static Operation* getOperatorOperation(const std::string& name); + static Operation* getFunctionOperation(const std::string& name, const std::map& customFunctions); +}; + +} // namespace Lepton + +#endif /*LEPTON_PARSER_H_*/ diff --git a/lib/colvars/lepton/include/lepton/windowsIncludes.h b/lib/colvars/lepton/include/lepton/windowsIncludes.h new file mode 100644 index 0000000000..798229850e --- /dev/null +++ b/lib/colvars/lepton/include/lepton/windowsIncludes.h @@ -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_ diff --git a/lib/colvars/lepton/src/CompiledExpression.cpp b/lib/colvars/lepton/src/CompiledExpression.cpp new file mode 100644 index 0000000000..302f294ee2 --- /dev/null +++ b/lib/colvars/lepton/src/CompiledExpression.cpp @@ -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 + +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 > 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 >& temps) { + if (findTempIndex(node, temps) != -1) + return; // We have already processed a node identical to this one. + + // Process the child nodes. + + vector args; + for (int i = 0; i < node.getChildren().size(); i++) { + compileExpression(node.getChildren()[i], temps); + args.push_back(findTempIndex(node.getChildren()[i], temps)); + } + + // Process this node. + + if (node.getOperation().getId() == Operation::VARIABLE) { + variableIndices[node.getOperation().getName()] = (int) workspace.size(); + variableNames.insert(node.getOperation().getName()); + } + else { + int stepIndex = (int) arguments.size(); + arguments.push_back(vector()); + 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 >& temps) { + for (int i = 0; i < (int) temps.size(); i++) + if (temps[i].first == node) + return i; + return -1; +} + +const set& CompiledExpression::getVariables() const { + return variableNames; +} + +double& CompiledExpression::getVariableReference(const string& name) { + map::iterator pointer = variablePointers.find(name); + if (pointer != variablePointers.end()) + return *pointer->second; + map::iterator index = variableIndices.find(name); + if (index == variableIndices.end()) + throw Exception("getVariableReference: Unknown variable '"+name+"'"); + return workspace[index->second]; +} + +void CompiledExpression::setVariableLocations(map& 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::const_iterator iter = variableIndices.begin(); iter != variableIndices.end(); ++iter) { + map::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& 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* dummyVariables = NULL; + return op->evaluate(args, *dummyVariables); +} + +void CompiledExpression::generateJitCode() { + X86Compiler c(&runtime); + c.addFunc(kFuncConvHost, FuncBuilder0()); + vector 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::const_iterator iter = variableNames.begin(); iter != variableNames.end(); ++iter) { + map::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 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(op).getValue(); + else if (op.getId() == Operation::ADD_CONSTANT) + value = dynamic_cast(op).getValue(); + else if (op.getId() == Operation::MULTIPLY_CONSTANT) + value = dynamic_cast(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 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 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()); + 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()); + call->setArg(0, arg); + call->setRet(0, dest); +} +#endif diff --git a/lib/colvars/lepton/src/ExpressionProgram.cpp b/lib/colvars/lepton/src/ExpressionProgram.cpp new file mode 100644 index 0000000000..65d3f0c79a --- /dev/null +++ b/lib/colvars/lepton/src/ExpressionProgram.cpp @@ -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()); +} + +double ExpressionProgram::evaluate(const std::map& variables) const { + vector 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]; +} diff --git a/lib/colvars/lepton/src/ExpressionTreeNode.cpp b/lib/colvars/lepton/src/ExpressionTreeNode.cpp new file mode 100644 index 0000000000..8519878262 --- /dev/null +++ b/lib/colvars/lepton/src/ExpressionTreeNode.cpp @@ -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& 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::getChildren() const { + return children; +} diff --git a/lib/colvars/lepton/src/MSVC_erfc.h b/lib/colvars/lepton/src/MSVC_erfc.h new file mode 100644 index 0000000000..eadb20fdf8 --- /dev/null +++ b/lib/colvars/lepton/src/MSVC_erfc.h @@ -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 + +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_ diff --git a/lib/colvars/lepton/src/Operation.cpp b/lib/colvars/lepton/src/Operation.cpp new file mode 100644 index 0000000000..693dea2ede --- /dev/null +++ b/lib/colvars/lepton/src/Operation.cpp @@ -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& variables) const { + return erf(args[0]); +} + +double Operation::Erfc::evaluate(double* args, const map& variables) const { + return erfc(args[0]); +} + +ExpressionTreeNode Operation::Constant::differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const { + return ExpressionTreeNode(new Operation::Constant(0.0)); +} + +ExpressionTreeNode Operation::Variable::differentiate(const std::vector& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& childDerivs, const std::string& variable) const { + return ExpressionTreeNode(new Operation::Add(), childDerivs[0], childDerivs[1]); +} + +ExpressionTreeNode Operation::Subtract::differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const { + return ExpressionTreeNode(new Operation::Subtract(), childDerivs[0], childDerivs[1]); +} + +ExpressionTreeNode Operation::Multiply::differentiate(const std::vector& children, const std::vector& 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& children, const std::vector& childDerivs, const std::string& variable) const { + return ExpressionTreeNode(new Operation::Divide(), + ExpressionTreeNode(new Operation::Subtract(), + ExpressionTreeNode(new Operation::Multiply(), children[1], childDerivs[0]), + ExpressionTreeNode(new Operation::Multiply(), children[0], childDerivs[1])), + ExpressionTreeNode(new Operation::Square(), children[1])); +} + +ExpressionTreeNode Operation::Power::differentiate(const std::vector& children, const std::vector& 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& children, const std::vector& childDerivs, const std::string& variable) const { + return ExpressionTreeNode(new Operation::Negate(), childDerivs[0]); +} + +ExpressionTreeNode Operation::Sqrt::differentiate(const std::vector& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& childDerivs, const std::string& variable) const { + return ExpressionTreeNode(new Operation::Constant(0.0)); +} + +ExpressionTreeNode Operation::Delta::differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const { + return ExpressionTreeNode(new Operation::Constant(0.0)); +} + +ExpressionTreeNode Operation::Square::differentiate(const std::vector& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& childDerivs, const std::string& variable) const { + return childDerivs[0]; +} + +ExpressionTreeNode Operation::MultiplyConstant::differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const { + return ExpressionTreeNode(new Operation::MultiplyConstant(value), + childDerivs[0]); +} + +ExpressionTreeNode Operation::PowerConstant::differentiate(const std::vector& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& 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& children, const std::vector& childDerivs, const std::string& variable) const { + return ExpressionTreeNode(new Operation::Constant(0.0)); +} + +ExpressionTreeNode Operation::Ceil::differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const { + return ExpressionTreeNode(new Operation::Constant(0.0)); +} + +ExpressionTreeNode Operation::Select::differentiate(const std::vector& children, const std::vector& childDerivs, const std::string& variable) const { + vector derivChildren; + derivChildren.push_back(children[0]); + derivChildren.push_back(childDerivs[1]); + derivChildren.push_back(childDerivs[2]); + return ExpressionTreeNode(new Operation::Select(), derivChildren); +} diff --git a/lib/colvars/lepton/src/ParsedExpression.cpp b/lib/colvars/lepton/src/ParsedExpression.cpp new file mode 100644 index 0000000000..6effd06007 --- /dev/null +++ b/lib/colvars/lepton/src/ParsedExpression.cpp @@ -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 +#include + +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()); +} + +double ParsedExpression::evaluate(const map& variables) const { + return evaluate(getRootNode(), variables); +} + +double ParsedExpression::evaluate(const ExpressionTreeNode& node, const map& variables) { + int numArgs = (int) node.getChildren().size(); + vector 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& 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& variables) { + if (node.getOperation().getId() == Operation::VARIABLE) { + const Operation::Variable& var = dynamic_cast(node.getOperation()); + map::const_iterator iter = variables.find(var.getName()); + if (iter == variables.end()) + return node; + return ExpressionTreeNode(new Operation::Constant(iter->second)); + } + vector 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 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()))); +} + +ExpressionTreeNode ParsedExpression::substituteSimplerExpression(const ExpressionTreeNode& node) { + vector 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(&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(&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(&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(&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(&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(&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(&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(&node.getOperation())->getValue()*dynamic_cast(&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(&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(&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 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(node.getOperation()).getValue(); + return numeric_limits::quiet_NaN(); +} + +ExpressionProgram ParsedExpression::createProgram() const { + return ExpressionProgram(*this); +} + +CompiledExpression ParsedExpression::createCompiledExpression() const { + return CompiledExpression(*this); +} + +ParsedExpression ParsedExpression::renameVariables(const map& replacements) const { + return ParsedExpression(renameNodeVariables(getRootNode(), replacements)); +} + +ExpressionTreeNode ParsedExpression::renameNodeVariables(const ExpressionTreeNode& node, const map& replacements) { + if (node.getOperation().getId() == Operation::VARIABLE) { + map::const_iterator replace = replacements.find(node.getOperation().getName()); + if (replace != replacements.end()) + return ExpressionTreeNode(new Operation::Variable(replace->second)); + } + vector 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; +} diff --git a/lib/colvars/lepton/src/Parser.cpp b/lib/colvars/lepton/src/Parser.cpp new file mode 100644 index 0000000000..6b19d7370d --- /dev/null +++ b/lib/colvars/lepton/src/Parser.cpp @@ -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 +#include + +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 Parser::tokenize(const string& expression) { + vector 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()); +} + +ParsedExpression Parser::parse(const string& expression, const map& customFunctions) { + try { + // First split the expression into subexpressions. + + string primaryExpression = expression; + vector 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 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 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 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& tokens, int& pos, const map& customFunctions, + const map& 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::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 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& customFunctions) { + + static map 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::const_iterator custom = customFunctions.find(trimmed); + if (custom != customFunctions.end()) + return new Operation::Custom(trimmed, custom->second->clone()); + + // Now try standard functions. + + map::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"); + } +}