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