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 "core/ActionRegister.h"
      23             : #include "Function.h"
      24             : 
      25             : namespace PLMD {
      26             : namespace function {
      27             : 
      28             : //+PLUMEDOC FUNCTION PIECEWISE
      29             : /*
      30             : Compute a piece wise straight line through its arguments that passes through a set of ordered control points.
      31             : 
      32             : For variables less than the first
      33             : (greater than the last) point, the value of the first (last) point is used.
      34             : 
      35             : \f[
      36             : \frac{y_{i+1}-y_i}{x_{i+1}-x_i}(s-x_i)+y_i ;  if x_i<s<x_{i+1}
      37             : \f]
      38             : \f[
      39             : y_N ; if x>x_{N-1}
      40             : \f]
      41             : \f[
      42             : y_1 ; if x<x_0
      43             : \f]
      44             : 
      45             : Control points are passed using the POINT0=... POINT1=... syntax as in the example below
      46             : 
      47             : If one argument is supplied, it results in a scalar quantity.
      48             : If multiple arguments are supplied, it results
      49             : in a vector of values. Each value will be named as the name of the original
      50             : argument with suffix _pfunc.
      51             : 
      52             : \par Examples
      53             : 
      54             : \plumedfile
      55             : dist1: DISTANCE ATOMS=1,10
      56             : dist2: DISTANCE ATOMS=2,11
      57             : 
      58             : pw: PIECEWISE POINT0=1,10 POINT1=2,PI POINT2=3,10 ARG=dist1
      59             : ppww: PIECEWISE POINT0=1,10 POINT1=2,PI POINT2=3,10 ARG=dist1,dist2
      60             : PRINT ARG=pw,ppww.dist1_pfunc,ppww.dist2_pfunc
      61             : \endplumedfile
      62             : 
      63             : 
      64             : */
      65             : //+ENDPLUMEDOC
      66             : 
      67             : 
      68             : class Piecewise :
      69             :   public Function {
      70             :   std::vector<std::pair<double,double> > points;
      71             : public:
      72             :   explicit Piecewise(const ActionOptions&);
      73             :   void calculate() override;
      74             :   static void registerKeywords(Keywords& keys);
      75             : };
      76             : 
      77             : 
      78       13790 : PLUMED_REGISTER_ACTION(Piecewise,"PIECEWISE")
      79             : 
      80           7 : void Piecewise::registerKeywords(Keywords& keys) {
      81           7 :   Function::registerKeywords(keys);
      82           7 :   keys.use("ARG");
      83          14 :   keys.add("numbered","POINT","This keyword is used to specify the various points in the function above.");
      84          14 :   keys.reset_style("POINT","compulsory");
      85           7 :   componentsAreNotOptional(keys);
      86          14 :   keys.addOutputComponent("_pfunc","default","one or multiple instances of this quantity can be referenced elsewhere "
      87             :                           "in the input file.  These quantities will be named with the arguments of the "
      88             :                           "function followed by the character string _pfunc.  These quantities tell the "
      89             :                           "user the values of the piece wise functions of each of the arguments.");
      90           7 : }
      91             : 
      92           3 : Piecewise::Piecewise(const ActionOptions&ao):
      93             :   Action(ao),
      94           3 :   Function(ao) {
      95           9 :   for(int i=0;; i++) {
      96             :     std::vector<double> pp;
      97          24 :     if(!parseNumberedVector("POINT",i,pp) ) {
      98             :       break;
      99             :     }
     100           9 :     if(pp.size()!=2) {
     101           0 :       error("points should be in x,y format");
     102             :     }
     103           9 :     points.push_back(std::pair<double,double>(pp[0],pp[1]));
     104           9 :     if(i>0 && points[i].first<=points[i-1].first) {
     105           0 :       error("points abscissas should be monotonously increasing");
     106             :     }
     107           9 :   }
     108             : 
     109           6 :   for(unsigned i=0; i<getNumberOfArguments(); i++)
     110           4 :     if(getPntrToArgument(i)->isPeriodic()) {
     111           1 :       error("Cannot use PIECEWISE on periodic arguments");
     112             :     }
     113             : 
     114           2 :   if(getNumberOfArguments()==1) {
     115           1 :     addValueWithDerivatives();
     116           1 :     setNotPeriodic();
     117             :   } else {
     118           3 :     for(unsigned i=0; i<getNumberOfArguments(); i++) {
     119           3 :       addComponentWithDerivatives( getPntrToArgument(i)->getName()+"_pfunc" );
     120           2 :       getPntrToComponent(i)->setNotPeriodic();
     121             :     }
     122             :   }
     123           2 :   checkRead();
     124             : 
     125           2 :   log.printf("  on points:");
     126           8 :   for(unsigned i=0; i<points.size(); i++) {
     127           6 :     log.printf("   (%f,%f)",points[i].first,points[i].second);
     128             :   }
     129           2 :   log.printf("\n");
     130           4 : }
     131             : 
     132          10 : void Piecewise::calculate() {
     133          25 :   for(unsigned i=0; i<getNumberOfArguments(); i++) {
     134             :     double val=getArgument(i);
     135             :     unsigned p=0;
     136          37 :     for(; p<points.size(); p++) {
     137          33 :       if(val<points[p].first) {
     138             :         break;
     139             :       }
     140             :     }
     141             :     double f,d;
     142          15 :     if(p==0) {
     143           5 :       f=points[0].second;
     144             :       d=0.0;
     145          10 :     } else if(p==points.size()) {
     146           4 :       f=points[points.size()-1].second;
     147             :       d=0.0;
     148             :     } else {
     149           6 :       double m=(points[p].second-points[p-1].second) / (points[p].first-points[p-1].first);
     150           6 :       f=m*(val-points[p-1].first)+points[p-1].second;
     151             :       d=m;
     152             :     }
     153          15 :     if(getNumberOfArguments()==1) {
     154           5 :       setValue(f);
     155           5 :       setDerivative(i,d);
     156             :     } else {
     157          10 :       Value* v=getPntrToComponent(i);
     158          10 :       v->set(f);
     159             :       v->addDerivative(i,d);
     160             :     }
     161             :   }
     162          10 : }
     163             : 
     164             : }
     165             : }
     166             : 
     167             :