Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2018 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 "tools/ConjugateGradient.h"
      24             : #include "ActionWithGrid.h"
      25             : #include "EvaluateGridFunction.h"
      26             : #include "Interpolator.h"
      27             : 
      28             : //+PLUMEDOC GRIDANALYSIS FIND_GRID_MINIMUM
      29             : /*
      30             : Find the point with the lowest value of the function on the grid
      31             : 
      32             : This command takes a function on a grid as input and returns multiple scalars.  One of the returned scalars (the one with component `optval`) is the value of the input
      33             : function at its lowest point. The other scalars returned are the coordinates for which the function takes this particular value.
      34             : 
      35             : As the input grid has the values of the function evaluated over a grid of points we find this minimum value by finding the point on the grid where the function has its lowest value using
      36             : interpolation and conjugate gradients.  The coordinates of the grid point where the function's value is lowest will be used as an initial coordinate
      37             : for the optimisation.  We then interpolate the function to find the optimum.  You can use use the CGTOL keyword to control the conjugate gradient algorithm that is used to find the location of the
      38             : minimum in the interpolated function.
      39             : 
      40             : To print out the location on the grid where the function is minimised and the value of the function at that point you can use an input like the one shown below:
      41             : 
      42             : ```plumed
      43             : x: DISTANCE ATOMS=1,2
      44             : hA1: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
      45             : ff: CONVERT_TO_FES ARG=hA1 TEMP=300
      46             : min: FIND_GRID_MINIMUM ARG=ff
      47             : PRINT ARG=min.x_opt,min.optval STRIDE=0
      48             : ```
      49             : 
      50             : Notice that we set STRIDE=0 in the PRINT command here so the position of the minimum is only output once at the end of the simulation.
      51             : 
      52             : A more common usage of this action is illustrated in the following input, which demonstrates how you can use FIND_GRID_MINIMUM to shift a free energy surface so that the minimum in the output
      53             : fes has a value of zero.
      54             : 
      55             : ```plumed
      56             : x: DISTANCE ATOMS=1,2
      57             : hA1: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
      58             : ff: CONVERT_TO_FES ARG=hA1 TEMP=300
      59             : min: FIND_GRID_MINIMUM ARG=ff NOINTERPOL
      60             : sff: CUSTOM ARG=ff,min.optval FUNC=x+y PERIODIC=NO
      61             : DUMPGRID ARG=sff FILE=fes.dat
      62             : ```
      63             : 
      64             : Notice that the `NOINTERPOL` flag has been used here so the interpolation and conjugate gradient algorithm is not used. The output is simply the grid point at which the function takes its lowest value.
      65             : 
      66             : */
      67             : //+ENDPLUMEDOC
      68             : 
      69             : //+PLUMEDOC GRIDANALYSIS FIND_GRID_MAXIMUM
      70             : /*
      71             : Find the point with the highest value of the function on the grid
      72             : 
      73             : This command takes a function on a grid as input and returns multiple scalars.  One of the returned scalars (the one with component `optval`) is the value of the input
      74             : function at its highest point. The other scalars returned are the coordinates for which the function takes this particular value.
      75             : 
      76             : As the input grid has the values of the function evaluated over a grid of points we find this maximum value by finding the point on the grid where the function has its lowest value using
      77             : interpolation and conjugate gradients.  The coordinates of the grid point where the function's value is highest will be used as an initial coordinate
      78             : for the optimisation.  We then interpolate the function to find the optimum.  You can use use the CGTOL keyword to control the conjugate gradient algorithm that is used to find the location of the
      79             : maximum in the interpolated function.
      80             : 
      81             : To print out the location on the grid where the function is maximised and the value of the function at that point you can use an input like the one shown below:
      82             : 
      83             : ```plumed
      84             : x: DISTANCE ATOMS=1,2
      85             : hA1: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
      86             : max: FIND_GRID_MAXIMUM ARG=hA1
      87             : PRINT ARG=max.x_opt,max.optval STRIDE=0
      88             : ```
      89             : 
      90             : Notice that we set STRIDE=0 in the PRINT command here so the position of the maximum is only output once at the end of the simulation.
      91             : 
      92             : If for any reason you do not want to use interpolation and conjugate gradient to find the optimum you can use the `NOINTERPOL` flag as shown below.
      93             : 
      94             : ```plumed
      95             : x: DISTANCE ATOMS=1,2
      96             : hA1: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
      97             : max: FIND_GRID_MAXIMUM ARG=hA1 NOINTERPOL
      98             : PRINT ARG=max.x_opt,max.optval STRIDE=0
      99             : ```
     100             : 
     101             : The FIND_GRID_MAXIMUM action then outputs the coordinates of the grid point where the function is largest as well as the value of the function on that grid point.
     102             : 
     103             : */
     104             : //+ENDPLUMEDOC
     105             : 
     106             : namespace PLMD {
     107             : namespace gridtools {
     108             : 
     109             : class FindGridOptimum :
     110             :   public ActionWithValue,
     111             :   public ActionWithArguments {
     112             : private:
     113             :   bool domin;
     114             :   double cgtol;
     115             :   std::vector<double> spacing;
     116             :   std::unique_ptr<Interpolator> function;
     117             :   double calculateValueAndDerivatives( const std::vector<double>& pp, std::vector<double>& der );
     118             : public:
     119             :   static void registerKeywords( Keywords& keys );
     120             :   explicit FindGridOptimum(const ActionOptions&ao);
     121             :   unsigned getNumberOfDerivatives() override ;
     122             :   void calculate() override ;
     123           0 :   void apply() override {}
     124             : };
     125             : 
     126             : PLUMED_REGISTER_ACTION(FindGridOptimum,"FIND_GRID_MAXIMUM")
     127             : PLUMED_REGISTER_ACTION(FindGridOptimum,"FIND_GRID_MINIMUM")
     128             : 
     129          12 : void FindGridOptimum::registerKeywords( Keywords& keys ) {
     130          12 :   Action::registerKeywords( keys );
     131          12 :   ActionWithValue::registerKeywords( keys );
     132          12 :   ActionWithArguments::registerKeywords( keys );
     133          12 :   keys.remove("NUMERICAL_DERIVATIVES");
     134          24 :   keys.addInputKeyword("compulsory","ARG","grid","the label for the function on the grid that you would like to find the optimum in");
     135          12 :   keys.addFlag("NOINTERPOL",false,"do not interpolate the function when finding the optimum");
     136          12 :   keys.add("compulsory","CGTOL","1E-4","the tolerance for the conjugate gradient optimization");
     137          24 :   keys.addOutputComponent("optval","default","scalar","the value of the function at the optimum");
     138          24 :   keys.addOutputComponent("_opt","default","scalar","the values of the arguments of the function at the optimum can be referenced elsewhere in the input file "
     139             :                           "by using the names of the arguments followed by the string _opt");
     140          12 : }
     141             : 
     142           2 : FindGridOptimum::FindGridOptimum(const ActionOptions&ao):
     143             :   Action(ao),
     144             :   ActionWithValue(ao),
     145             :   ActionWithArguments(ao),
     146           2 :   cgtol(0) {
     147           2 :   if( getName()=="FIND_GRID_MAXIMUM" ) {
     148           0 :     domin=false;
     149           2 :   } else if( getName()=="FIND_GRID_MINIMUM" ) {
     150           2 :     domin=true;
     151             :   } else {
     152           0 :     plumed_error();
     153             :   }
     154             :   // Create value for this function
     155           2 :   ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
     156           2 :   if( !ag ) {
     157           0 :     error("input action should be a grid");
     158             :   }
     159           2 :   std::vector<std::string> argn( ag->getGridCoordinateNames() );
     160           2 :   std::vector<std::size_t> shape(0);
     161           5 :   for(unsigned i=0; i<argn.size(); ++i) {
     162           6 :     addComponent( argn[i] + "_opt", shape );
     163           6 :     componentIsNotPeriodic( argn[i] + "_opt" );
     164             :   }
     165           2 :   addComponent( "optval", shape );
     166           2 :   componentIsNotPeriodic( "optval" );
     167           2 :   bool nointerpol=false;
     168           2 :   parseFlag("NOINTERPOL",nointerpol);
     169           2 :   if( !nointerpol ) {
     170           2 :     parse("CGTOL",cgtol);
     171           4 :     function=Tools::make_unique<Interpolator>( getPntrToArgument(0), ag->getGridCoordinatesObject() );
     172             :   }
     173           2 : }
     174             : 
     175           0 : unsigned FindGridOptimum::getNumberOfDerivatives() {
     176           0 :   return 0;
     177             : }
     178             : 
     179         205 : double FindGridOptimum::calculateValueAndDerivatives( const std::vector<double>& pp, std::vector<double>& der ) {
     180         205 :   double val = function->splineInterpolation( pp, der );
     181             :   // We normalise the derivatives here and set them so that the linesearch is done over the cell that we know
     182             :   // in the grid that we know the minimum is inside
     183             :   double norm = 0;
     184         560 :   for(unsigned i=0; i<der.size(); ++i) {
     185         355 :     norm += der[i]*der[i];
     186             :   }
     187         205 :   norm = sqrt(norm);
     188         560 :   for(unsigned i=0; i<der.size(); ++i) {
     189         355 :     der[i] = spacing[i]*der[i] / norm;
     190             :   }
     191         205 :   if( domin ) {
     192             :     return val;
     193             :   }
     194             :   // If we are looking for a maximum
     195           0 :   for(unsigned i=0; i<der.size(); ++i) {
     196           0 :     der[i] = -der[i];
     197             :   }
     198           0 :   return -val;
     199             : }
     200             : 
     201           2 : void FindGridOptimum::calculate() {
     202           2 :   ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
     203           2 :   plumed_assert( ag );
     204           2 :   spacing = ag->getGridCoordinatesObject().getGridSpacing();
     205           2 :   const GridCoordinatesObject& ingrid = ag->getGridCoordinatesObject();
     206             :   Value* gval = getPntrToArgument(0);
     207           2 :   std::vector<double> optargs( gval->getRank() );
     208           2 :   std::vector<unsigned> gridind( gval->getRank() );
     209           2 :   double optval=gval->get( 0 );
     210           2 :   ingrid.getGridPointCoordinates( 0, gridind, optargs );
     211           2 :   unsigned nval = gval->getNumberOfValues();
     212             :   bool constant=true;
     213        2602 :   for(unsigned i=0; i<nval; ++i) {
     214        2600 :     double tval = gval->get( i );
     215        2600 :     if( domin && (tval<optval || std::isnan(optval)) ) {
     216             :       constant=false;
     217             :       optval=tval;
     218          33 :       ingrid.getGridPointCoordinates( i, gridind, optargs );
     219             :     }
     220        2600 :     if( !domin && (tval>optval || std::isnan(optval)) ) {
     221             :       constant=false;
     222             :       optval=tval;
     223           0 :       ingrid.getGridPointCoordinates( i, gridind, optargs );
     224             :     }
     225             :   }
     226             :   // This basically ensures we deal with cases where all points on the grid are infinity as isinf doesn't work on intel compiler
     227           2 :   if( constant ) {
     228           0 :     if( domin && gval->get(0)>=gval->get(1) ) {
     229             :       return;
     230           0 :     } else if( gval->get(0)<=gval->get(1) ) {
     231             :       return;
     232             :     }
     233             :   }
     234           2 :   if( std::isinf(optval) ) {
     235             :     return;
     236             :   }
     237             : 
     238           2 :   if( std::isnan(optval) ) {
     239           0 :     error("all values on grid are nans");
     240             :   }
     241             :   // Set the value of the arguments at the optimum before we optimise just in case the optimisation fails
     242           5 :   for(unsigned j=0; j<optargs.size(); ++j) {
     243           3 :     getPntrToComponent(j)->set( optargs[j] );
     244             :   }
     245             : 
     246             :   // And do conjugate gradient optimisation (because we can!!)
     247           2 :   if( cgtol>0 ) {
     248             :     ConjugateGradient<FindGridOptimum> myminimiser( this );
     249           2 :     int code = myminimiser.minimise( cgtol, optargs, &FindGridOptimum::calculateValueAndDerivatives );
     250           2 :     if( code>0 ) {
     251           0 :       warning("could not find optimum for function on input grid. Am thus just returning point in grid where value of function is lowest");
     252           0 :       getPntrToComponent(optargs.size())->set( optval );
     253           0 :       return;
     254             :     }
     255             :   }
     256             :   // And set the final value
     257           5 :   for(unsigned j=0; j<optargs.size(); ++j) {
     258           3 :     getPntrToComponent(j)->set( optargs[j] );
     259             :   }
     260           2 :   std::vector<double> optder( gval->getRank() );
     261           2 :   if( function ) {
     262           2 :     getPntrToComponent(optargs.size())->set( function->splineInterpolation( optargs, optder ) );
     263             :   } else {
     264           0 :     getPntrToComponent(optargs.size())->set( optval );
     265             :   }
     266           4 : }
     267             : 
     268             : }
     269             : }