Using lepton in MATHEVAL

Notice that Lepton is always available but should be triggered with
export PLUMED_USE_LEPTON=yes

Also notice that to behave as matheval it requires defining a pi variable.

src/function/Makefile

-USE=core reference tools
+USE=core reference tools lepton
 
 #generic makefile
 include ../maketools/make.module

src/function/Matheval.cpp

@@ -26,6 +26,8 @@
 #include <matheval.h>
 #endif
 
+#include "lepton/Lepton.h"
+
 using namespace std;
 
 namespace PLMD {
@@ -170,6 +172,7 @@ progression (S) and distance (Z) variables \cite perez2015atp.
 class Matheval :
   public Function
 {
+  const bool use_lepton;
   void* evaluator;
   vector<void*> evaluator_deriv;
   vector<string> var;
@@ -183,7 +186,6 @@ public:
   static void registerKeywords(Keywords& keys);
 };
 
-#ifdef __PLUMED_HAS_MATHEVAL
 PLUMED_REGISTER_ACTION(Matheval,"MATHEVAL")
 
 void Matheval::registerKeywords(Keywords& keys) {
@@ -196,6 +198,7 @@ void Matheval::registerKeywords(Keywords& keys) {
 Matheval::Matheval(const ActionOptions&ao):
   Action(ao),
   Function(ao),
+  use_lepton(std::getenv("PLUMED_USE_LEPTON")),
   evaluator_deriv(getNumberOfArguments()),
   values(getNumberOfArguments()),
   names(getNumberOfArguments())
@@ -215,56 +218,119 @@ Matheval::Matheval(const ActionOptions&ao):
   addValueWithDerivatives();
   checkRead();
 
-  evaluator=evaluator_create(const_cast<char*>(func.c_str()));
+  if(use_lepton) {
+    func=func+" ; pi=3.14159265358979";
+    log<<"  WARNING: you are using lepton as a replacement for libmatheval\n";
+    evaluator=new lepton::CompiledExpression(lepton::Parser::parse(func).optimize().createCompiledExpression());
+  } else {
+#ifdef __PLUMED_HAS_MATHEVAL
+    evaluator=evaluator_create(const_cast<char*>(func.c_str()));
+#else
+    error("MATHEVAL not available, please export PLUMED_USE_LEPTON=yes");
+#endif
+  }
 
   if(!evaluator) error("There was some problem in parsing matheval formula "+func);
 
-  char **check_names;
-  int    check_count;
-  evaluator_get_variables(evaluator,&check_names,&check_count);
-  if(check_count!=int(getNumberOfArguments())) {
-    string sc;
-    Tools::convert(check_count,sc);
-    error("Your function string contains "+sc+" arguments. This should be equal to the number of ARGs");
-  }
-  for(unsigned i=0; i<getNumberOfArguments(); i++) {
-    bool found=false;
-    for(unsigned j=0; j<getNumberOfArguments(); j++) {
-      if(var[i]==check_names[j])found=true;
+  if(use_lepton) {
+    for(unsigned i=0; i<getNumberOfArguments(); i++)
+      evaluator_deriv[i]=new
+      lepton::CompiledExpression(lepton::Parser::parse(func).differentiate(var[i]).optimize().createCompiledExpression());
+  } else {
+#ifdef __PLUMED_HAS_MATHEVAL
+    char **check_names;
+    int    check_count;
+    evaluator_get_variables(evaluator,&check_names,&check_count);
+    if(check_count!=int(getNumberOfArguments())) {
+      string sc;
+      Tools::convert(check_count,sc);
+      error("Your function string contains "+sc+" arguments. This should be equal to the number of ARGs");
     }
-    if(!found)
-      error("Variable "+var[i]+" cannot be found in your function string");
+    for(unsigned i=0; i<getNumberOfArguments(); i++) {
+      bool found=false;
+      for(unsigned j=0; j<getNumberOfArguments(); j++) {
+        if(var[i]==check_names[j])found=true;
+      }
+      if(!found)
+        error("Variable "+var[i]+" cannot be found in your function string");
+    }
+    for(unsigned i=0; i<getNumberOfArguments(); i++)
+      evaluator_deriv[i]=evaluator_derivative(evaluator,const_cast<char*>(var[i].c_str()));
+#else
+    error("MATHEVAL not available, please export PLUMED_USE_LEPTON=yes");
+#endif
   }
 
-  for(unsigned i=0; i<getNumberOfArguments(); i++)
-    evaluator_deriv[i]=evaluator_derivative(evaluator,const_cast<char*>(var[i].c_str()));
-
-
   log.printf("  with function : %s\n",func.c_str());
   log.printf("  with variables :");
   for(unsigned i=0; i<var.size(); i++) log.printf(" %s",var[i].c_str());
   log.printf("\n");
-  log.printf("  function as parsed by matheval: %s\n", evaluator_get_string(evaluator));
-  log.printf("  derivatives as computed by matheval:\n");
-  for(unsigned i=0; i<var.size(); i++) log.printf("    %s\n",evaluator_get_string(evaluator_deriv[i]));
+  if(use_lepton) {
+    log<<"  function as parsed by lepton: "<<lepton::Parser::parse(func).optimize()<<"\n";
+    log<<"  derivatives as computed by lepton:\n";
+    for(unsigned i=0; i<var.size(); i++) log<<"    "<<lepton::Parser::parse(func).differentiate(var[i]).optimize()<<"\n";
+  } else {
+#ifdef __PLUMED_HAS_MATHEVAL
+    log.printf("  function as parsed by matheval: %s\n", evaluator_get_string(evaluator));
+    log.printf("  derivatives as computed by matheval:\n");
+    for(unsigned i=0; i<var.size(); i++) log.printf("    %s\n",evaluator_get_string(evaluator_deriv[i]));
+#else
+    error("MATHEVAL not available, please export PLUMED_USE_LEPTON=yes");
+#endif
+  }
 }
 
 void Matheval::calculate() {
-  for(unsigned i=0; i<getNumberOfArguments(); i++) values[i]=getArgument(i);
-  for(unsigned i=0; i<getNumberOfArguments(); i++) names[i]=const_cast<char*>(var[i].c_str());
-  setValue(evaluator_evaluate(evaluator,names.size(),&names[0],&values[0]));
-
-  for(unsigned i=0; i<getNumberOfArguments(); i++) {
-    setDerivative(i,evaluator_evaluate(evaluator_deriv[i],names.size(),&names[0],&values[0]));
+  if(use_lepton) {
+    for(unsigned i=0; i<getNumberOfArguments(); i++) {
+      try {
+        static_cast<lepton::CompiledExpression*>(evaluator)->getVariableReference(var[i])=getArgument(i);
+      } catch(PLMD::lepton::Exception& exc) {
+// this is necessary since in some cases lepton things a variable is not present even though it is present
+// e.g. func=0*x
+      }
+    }
+    setValue(static_cast<lepton::CompiledExpression*>(evaluator)->evaluate());
+    for(unsigned i=0; i<getNumberOfArguments(); i++) {
+      for(unsigned j=0; j<getNumberOfArguments(); j++) {
+        try {
+          static_cast<lepton::CompiledExpression*>(evaluator_deriv[i])->getVariableReference(var[j])=getArgument(j);
+        } catch(PLMD::lepton::Exception& exc) {
+// this is necessary since in some cases lepton things a variable is not present even though it is present
+// e.g. func=0*x
+        }
+      }
+      setDerivative(i,static_cast<lepton::CompiledExpression*>(evaluator_deriv[i])->evaluate());
+    }
+  } else {
+#ifdef __PLUMED_HAS_MATHEVAL
+    for(unsigned i=0; i<getNumberOfArguments(); i++) values[i]=getArgument(i);
+    for(unsigned i=0; i<getNumberOfArguments(); i++) names[i]=const_cast<char*>(var[i].c_str());
+    setValue(evaluator_evaluate(evaluator,names.size(),&names[0],&values[0]));
+    for(unsigned i=0; i<getNumberOfArguments(); i++) {
+      setDerivative(i,evaluator_evaluate(evaluator_deriv[i],names.size(),&names[0],&values[0]));
+    }
+#else
+    error("MATHEVAL not available, please export PLUMED_USE_LEPTON=yes");
+#endif
   }
 }
 
 Matheval::~Matheval() {
-  evaluator_destroy(evaluator);
-  for(unsigned i=0; i<evaluator_deriv.size(); i++)evaluator_destroy(evaluator_deriv[i]);
-}
-
+  if(use_lepton) {
+    for(unsigned i=0; i<evaluator_deriv.size(); i++) {
+      delete static_cast<lepton::CompiledExpression*>(evaluator_deriv[i]);
+    }
+    delete static_cast<lepton::CompiledExpression*>(evaluator);
+  } else {
+#ifdef __PLUMED_HAS_MATHEVAL
+    evaluator_destroy(evaluator);
+    for(unsigned i=0; i<evaluator_deriv.size(); i++)evaluator_destroy(evaluator_deriv[i]);
+#else
+    error("MATHEVAL not available, please export PLUMED_USE_LEPTON=yes");
 #endif
+  }
+}
 
 }
 }