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           0 : LeptonCall::LeptonCall(const LeptonCall& other):
      28           0 :   nargs(other.nargs),
      29           0 :   allow_extra_args(other.allow_extra_args),
      30           0 :   varnames(other.varnames),
      31           0 :   expression(other.expression),
      32           0 :   expression_deriv(other.expression_deriv) {
      33           0 :   linkVariables();
      34           0 : }
      35             : 
      36           0 : LeptonCall& LeptonCall::operator=(const LeptonCall& other) {
      37           0 :   if( this!=&other ) {
      38           0 :     nargs=other.nargs;
      39           0 :     allow_extra_args=other.allow_extra_args;
      40           0 :     varnames=other.varnames;
      41           0 :     expression=other.expression;
      42           0 :     expression_deriv=other.expression_deriv;
      43           0 :     linkVariables();
      44             :   }
      45           0 :   return *this;
      46             : }
      47             : 
      48        3492 : void LeptonCall::linkVariables() {
      49             :   const std::size_t nth=expression.size();
      50        3492 :   const std::size_t na=nargs;
      51        3492 :   lepton_ref.assign(nth*na,nullptr);
      52       10476 :   for(std::size_t t=0; t<nth; ++t) {
      53       18022 :     for(std::size_t j=0; j<na; ++j) {
      54             :       try {
      55       11038 :         lepton_ref[t*na+j]=&expression[t].getVariableReference(varnames[j]);
      56         154 :       } catch(const PLMD::lepton::Exception& exc) {
      57             : // this is necessary since in some cases lepton things a variable is not present even though it is present
      58             : // e.g. func=0*x
      59         154 :       }
      60             :     }
      61             :   }
      62        3492 :   lepton_ref_deriv.assign(nth*na*na,nullptr);
      63        9011 :   for(std::size_t i=0; i<expression_deriv.size(); i++) {
      64        5519 :     plumed_assert( expression_deriv[i].size()==nth );
      65       16557 :     for(std::size_t t=0; t<nth; ++t) {
      66       31776 :       for(std::size_t j=0; j<na; ++j) {
      67             :         try {
      68       20738 :           lepton_ref_deriv[i*nth*na + t*na+j]=&expression_deriv[i][t].getVariableReference(varnames[j]);
      69        7396 :         } catch(const PLMD::lepton::Exception& exc) {
      70        7396 :         }
      71             :       }
      72             :     }
      73             :   }
      74        3492 : }
      75             : 
      76        3492 : void LeptonCall::set(const std::string & func,
      77             :                      const std::vector<std::string>& var,
      78             :                      Action* action,
      79             :                      const bool extraArgs ) {
      80        3492 :   allow_extra_args=extraArgs;
      81        3492 :   varnames=var;
      82        3492 :   nargs=var.size();
      83        3492 :   const unsigned nth=OpenMP::getNumThreads();
      84        3492 :   lepton::ParsedExpression pe=lepton::Parser::parse(func).optimize(lepton::Constants());
      85        3492 :   if( action ) {
      86        3383 :     action->log<<"  function as parsed by lepton: "<<pe<<"\n";
      87             :   }
      88        3492 :   expression.assign(nth,pe.createCompiledExpression());
      89        8934 :   for(auto & p : expression[0].getVariables()) {
      90        5442 :     if(std::find(var.begin(),var.end(),p)==var.end()) {
      91           0 :       if( action ) {
      92           0 :         action->error("variable " + p + " is not defined");
      93             :       } else {
      94           0 :         plumed_merror("variable " + p + " is not defined in lepton function");
      95             :       }
      96             :     }
      97             :   }
      98        3492 :   if( action ) {
      99        3383 :     action->log<<"  derivatives as computed by lepton:\n";
     100             :   }
     101        3492 :   expression_deriv.resize(var.size());
     102        9011 :   for(unsigned i=0; i<var.size(); i++) {
     103        5519 :     lepton::ParsedExpression pevar=lepton::Parser::parse(func).differentiate(var[i]).optimize(lepton::Constants());
     104        5519 :     if( action ) {
     105        5301 :       action->log<<"    "<<pevar<<"\n";
     106             :     }
     107       11038 :     expression_deriv[i].assign(nth,pevar.createCompiledExpression());
     108             :   }
     109        3492 :   linkVariables();
     110        3492 : }
     111             : 
     112      551477 : double LeptonCall::evaluate( const std::vector<double>& args ) const {
     113      551477 :   return evaluate( View<const double>( args.data(), args.size() ) );
     114             : }
     115             : 
     116    41922378 : double LeptonCall::evaluate(const View<const double> args ) const {
     117             :   plumed_dbg_assert( allow_extra_args || args.size()==nargs );
     118    41922378 :   const std::size_t t=OpenMP::getThreadNum(), tbas=t*nargs;
     119             :   plumed_dbg_assert( t<expression.size() );
     120   123617773 :   for(unsigned i=0; i<nargs; ++i) {
     121    81695395 :     if( lepton_ref[tbas+i] ) {
     122    81001259 :       *lepton_ref[tbas+i] = args[i];
     123             :     }
     124             :   }
     125    41922378 :   return expression[t].evaluate();
     126             : }
     127             : 
     128      296550 : double LeptonCall::evaluateDeriv( const unsigned ider, const std::vector<double>& args ) const {
     129      296550 :   return evaluateDeriv( ider, View<const double>( args.data(), args.size() ) );
     130             : }
     131             : 
     132    30057796 : double LeptonCall::evaluateDeriv( const unsigned ider, const View<const double> args ) const {
     133             :   plumed_dbg_assert( allow_extra_args || args.size()==nargs );
     134             :   plumed_dbg_assert( ider<nargs );
     135    30057796 :   const std::size_t t=OpenMP::getThreadNum();
     136    30057796 :   const std::size_t dbas = static_cast<std::size_t>(ider)*expression.size()*nargs + t*nargs;
     137    95965896 :   for(unsigned j=0; j<nargs; j++) {
     138    65908100 :     if(lepton_ref_deriv[dbas+j] ) {
     139    43373050 :       *lepton_ref_deriv[dbas+j] = args[j];
     140             :     }
     141             :   }
     142    30057796 :   return expression_deriv[ider][t].evaluate();
     143             : }
     144             : 
     145             : }