Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2012-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "LeptonCall.h"
      23             : #include "OpenMP.h"
      24             : 
      25             : namespace PLMD {
      26             : 
      27        3082 : void LeptonCall::set(const std::string & func, const std::vector<std::string>& var, Action* action, const bool& a ) {
      28        3082 :   unsigned nth=OpenMP::getNumThreads();
      29        3082 :   expression.resize(nth);
      30        3082 :   expression_deriv.resize(var.size());
      31             :   // Resize the expression for the derivatives
      32        8067 :   for(unsigned i=0; i<expression_deriv.size(); ++i) {
      33        4985 :     expression_deriv[i].resize(OpenMP::getNumThreads());
      34             :   }
      35        3082 :   allow_extra_args=a;
      36        3082 :   nargs=var.size();
      37             : 
      38        3082 :   lepton_ref.resize(nth*nargs,nullptr);
      39        3082 :   lepton::ParsedExpression pe=lepton::Parser::parse(func).optimize(lepton::Constants());
      40             :   unsigned nt=0;
      41        3082 :   if( action ) {
      42        3082 :     action->log<<"  function as parsed by lepton: "<<pe<<"\n";
      43             :   }
      44        9246 :   for(auto & e : expression) {
      45        6164 :     e=pe.createCompiledExpression();
      46       16134 :     for(unsigned j=0; j<var.size(); ++j) {
      47             :       try {
      48        9970 :         lepton_ref[nt*var.size()+j]=&const_cast<lepton::CompiledExpression*>(&expression[nt])->getVariableReference(var[j]);
      49          88 :       } catch(const PLMD::lepton::Exception& exc) {
      50             : // this is necessary since in some cases lepton things a variable is not present even though it is present
      51             : // e.g. func=0*x
      52          88 :       }
      53             :     }
      54        6164 :     nt++;
      55             :   }
      56        8023 :   for(auto & p : expression[0].getVariables()) {
      57        4941 :     if(std::find(var.begin(),var.end(),p)==var.end()) {
      58           0 :       if( action ) {
      59           0 :         action->error("variable " + p + " is not defined");
      60             :       } else {
      61           0 :         plumed_merror("variable " + p + " is not defined in lepton function");
      62             :       }
      63             :     }
      64             :   }
      65        3082 :   if( action ) {
      66        3082 :     action->log<<"  derivatives as computed by lepton:\n";
      67             :   }
      68        3082 :   lepton_ref_deriv.resize(nth*nargs*nargs,nullptr);
      69        8067 :   for(unsigned i=0; i<var.size(); i++) {
      70        9970 :     lepton::ParsedExpression pe=lepton::Parser::parse(func).differentiate(var[i]).optimize(lepton::Constants());
      71             :     nt=0;
      72        4985 :     if( action ) {
      73        4985 :       action->log<<"    "<<pe<<"\n";
      74             :     }
      75       14955 :     for(auto & e : expression_deriv[i]) {
      76        9970 :       e=pe.createCompiledExpression();
      77       28896 :       for(unsigned j=0; j<var.size(); ++j) {
      78             :         try {
      79       18926 :           lepton_ref_deriv[i*OpenMP::getNumThreads()*var.size() + nt*var.size()+j]=&const_cast<lepton::CompiledExpression*>(&expression_deriv[i][nt])->getVariableReference(var[j]);
      80        6110 :         } catch(const PLMD::lepton::Exception& exc) {
      81             : // this is necessary since in some cases lepton things a variable is not present even though it is present
      82             : // e.g. func=0*x
      83        6110 :         }
      84             :       }
      85        9970 :       nt++;
      86             :     }
      87             :   }
      88        3082 : }
      89             : 
      90    22966063 : double LeptonCall::evaluate( const std::vector<double>& args ) const {
      91             :   plumed_dbg_assert( allow_extra_args || args.size()==nargs );
      92    22966063 :   const unsigned t=OpenMP::getThreadNum(), tbas=t*nargs;
      93    67582765 :   for(unsigned i=0; i<nargs; ++i) {
      94    44616702 :     if( lepton_ref[tbas+i] ) {
      95    44049736 :       *lepton_ref[tbas+i] = args[i];
      96             :     }
      97             :   }
      98    22966063 :   return expression[t].evaluate();
      99             : }
     100             : 
     101    32189400 : double LeptonCall::evaluateDeriv( const unsigned& ider, const std::vector<double>& args ) const {
     102             :   plumed_dbg_assert( allow_extra_args || args.size()==nargs );
     103             :   plumed_dbg_assert( ider<nargs );
     104    32189400 :   const unsigned t=OpenMP::getThreadNum(), dbas = ider*OpenMP::getNumThreads()*nargs + t*nargs;
     105   106819836 :   for(unsigned j=0; j<nargs; j++) {
     106    74630436 :     if(lepton_ref_deriv[dbas+j] ) {
     107    52954445 :       *lepton_ref_deriv[dbas+j] = args[j];
     108             :     }
     109             :   }
     110    32189400 :   return expression_deriv[ider][t].evaluate();
     111             : }
     112             : 
     113             : }