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 "FunctionSetup.h"
      24             : #include "FunctionShortcut.h"
      25             : #include "FunctionOfScalar.h"
      26             : 
      27             : namespace PLMD {
      28             : namespace function {
      29             : 
      30             : //+PLUMEDOC FUNCTION PIECEWISE
      31             : /*
      32             : Compute a piece wise straight line through its arguments that passes through a set of ordered control points.
      33             : 
      34             : This action can be used to calculate a piecewise linear function of an input argument such as the one given below:
      35             : 
      36             : $$
      37             : f(x) = \begin{cases}
      38             : 10 & \textrm{if} \quad x<1 \\
      39             : 10 + \frac{\pi - 10}{2-1}(x-1) & \textrm{if} \quad 1 \le x < 2 \\
      40             : \pi + \frac{10 - \pi}{3-2}(x-2) & \textrm{if} \quad 2 \le x \le 3 \\
      41             : 10 & \textrm{otherwise}
      42             : \end{cases}
      43             : $$
      44             : 
      45             : The example shown below illustrates how one can use PLUMED to evaluate the function described above for the distance
      46             : between atom 1 and atom 2.
      47             : 
      48             : ```plumed
      49             : dist1: DISTANCE ATOMS=1,10
      50             : pw: PIECEWISE POINT0=1,10 POINT1=2,PI POINT2=3,10 ARG=dist1
      51             : PRINT ARG=pw FILE=colvar
      52             : ```
      53             : 
      54             : As you can see from the example above, the control points for the picewise function are passed using the `POINT0=...` `POINT1=...` syntax.
      55             : You can specify as many of these control points as you want.  These control points then serce as the $x_i$ and $y_i$ values in the following expression.
      56             : 
      57             : For variables less than the first
      58             : (greater than the last) point, the value of the first (last) point is used.
      59             : 
      60             : $$
      61             : \frac{y_{i+1}-y_i}{x_{i+1}-x_i}(s-x_i)+y_i ;  if x_i<s<x_{i+1}
      62             : $$
      63             : 
      64             : If the input value of $s$ is smaller than the lowest specified $x_i$ value then this action outputs the $y_i$ value that corresponds to the smallest of the input $x_i$ values.
      65             : Similarly, if the input value of $s$ is larger than the highest specified $x_i$ value the $y_i$ value that corresponds to the largest of the input $x_i$ values is output.
      66             : 
      67             : ## Using multiple scalars in input
      68             : 
      69             : The following example illustrates what happens when multiple scalar arguments are passed to this action:
      70             : 
      71             : ```plumed
      72             : dist1: DISTANCE ATOMS=1,10
      73             : dist2: DISTANCE ATOMS=2,11
      74             : 
      75             : ppww: PIECEWISE POINT0=1,10 POINT1=2,PI POINT2=3,10 ARG=dist1,dist2
      76             : PRINT ARG=ppww.dist1_pfunc,ppww.dist2_pfunc
      77             : ```
      78             : 
      79             : In essence the piecewise function is applied to each of the input arguments in turn.  Hence, in the example above the PIECEWISE command outputs two values.  The first of
      80             : these `ppww.dist1_pfunc` is the result that is obtained when the piecewise function is applied to the argument `dist1`.  The second is then the result that is obtained when the
      81             : piecewise function is applied to the argument `dist2`.
      82             : 
      83             : ## Non rank zero arguments
      84             : 
      85             : This argument currently cannot accept non-rank zero arguments.  However, it would be extremely straightforward to add functionality to ensure that if a PIECEWISE command receives
      86             : a vector, matrix or function on grid in input it will output a vector, matrix or function on grid that is obtained by applying the piecewise function elementwise to each of the
      87             : elements of the input vector, matrix or function.
      88             : 
      89             : */
      90             : //+ENDPLUMEDOC
      91             : 
      92           4 : class Piecewise {
      93             : public:
      94             :   std::vector<std::pair<double,double> > points;
      95             :   static void registerKeywords(Keywords& keys);
      96             :   static void read( Piecewise& func, ActionWithArguments* action, FunctionOptions& options );
      97             :   static void calc( const Piecewise& func, bool noderiv, View<const double> args, FunctionOutput& funcout );
      98             : };
      99             : 
     100             : 
     101             : typedef FunctionShortcut<Piecewise> PiecewiseShortcut;
     102             : PLUMED_REGISTER_ACTION(PiecewiseShortcut,"PIECEWISE")
     103             : typedef FunctionOfScalar<Piecewise> ScalarPiecewise;
     104             : PLUMED_REGISTER_ACTION(ScalarPiecewise,"PIECEWISE_SCALAR")
     105             : 
     106          16 : void Piecewise::registerKeywords(Keywords& keys) {
     107          16 :   keys.add("numbered","POINT","This keyword is used to specify the various points in the function above.");
     108          32 :   keys.reset_style("POINT","compulsory");
     109          32 :   keys.addOutputComponent("_pfunc","default","scalar","one or multiple instances of this quantity can be referenced elsewhere "
     110             :                           "in the input file.  These quantities will be named with the arguments of the "
     111             :                           "function followed by the character string _pfunc.  These quantities tell the "
     112             :                           "user the values of the piece wise functions of each of the arguments.");
     113          16 : }
     114             : 
     115           3 : void Piecewise::read( Piecewise& func, ActionWithArguments* action, FunctionOptions& options ) {
     116           9 :   for(int i=0;; i++) {
     117          12 :     std::vector<double> pp;
     118          24 :     if(!action->parseNumberedVector("POINT",i,pp) ) {
     119             :       break;
     120             :     }
     121           9 :     if(pp.size()!=2) {
     122           0 :       action->error("points should be in x,y format");
     123             :     }
     124           9 :     func.points.push_back(std::pair<double,double>(pp[0],pp[1]));
     125           9 :     if(i>0 && func.points[i].first<=func.points[i-1].first) {
     126           0 :       action->error("points abscissas should be monotonously increasing");
     127             :     }
     128          12 :   }
     129             : 
     130           3 :   if( func.points.empty() ) {
     131           0 :     action->error("no points found: use POINT0=x,y POINT1=x,y ...");
     132             :   }
     133           6 :   for(unsigned i=0; i<action->getNumberOfArguments(); i++) {
     134           4 :     if(action->getPntrToArgument(i)->isPeriodic()) {
     135           1 :       action->error("Cannot use PIECEWISE on periodic arguments");
     136             :     }
     137             :   }
     138           2 :   action->log.printf("  on points:");
     139           8 :   for(unsigned i=0; i<func.points.size(); i++) {
     140           6 :     action->log.printf("   (%f,%f)",func.points[i].first,func.points[i].second);
     141             :   }
     142           2 :   action->log.printf("\n");
     143           2 : }
     144             : 
     145          10 : void Piecewise::calc( const Piecewise& func, bool noderiv, const View<const double> args, FunctionOutput& funcout ) {
     146          25 :   for(unsigned i=0; i<args.size(); i++) {
     147             :     unsigned p=0;
     148          37 :     for(; p<func.points.size(); p++) {
     149          33 :       if(args[i]<func.points[p].first) {
     150             :         break;
     151             :       }
     152             :     }
     153          15 :     if(p==0) {
     154           5 :       funcout.values[i]=func.points[0].second;
     155           5 :       if( !noderiv ) {
     156           5 :         funcout.derivs[i][i]=0.0;
     157             :       }
     158          10 :     } else if(p==func.points.size()) {
     159           4 :       funcout.values[i]=func.points[func.points.size()-1].second;
     160           4 :       if( !noderiv ) {
     161           4 :         funcout.derivs[i][i]=0.0;
     162             :       }
     163             :     } else {
     164           6 :       double m=(func.points[p].second-func.points[p-1].second) / (func.points[p].first-func.points[p-1].first);
     165           6 :       funcout.values[i]=m*(args[i]-func.points[p-1].first)+func.points[p-1].second;
     166           6 :       if( !noderiv ) {
     167           6 :         funcout.derivs[i][i]=m;
     168             :       }
     169             :     }
     170             :   }
     171          10 : }
     172             : 
     173             : }
     174             : }
     175             : 
     176             :