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 GRIDCALC FIND_GRID_MINIMUM
      29             : /*
      30             : Find the point with the lowest value of the function on the grid
      31             : 
      32             : \par Examples
      33             : 
      34             : */
      35             : //+ENDPLUMEDOC
      36             : 
      37             : //+PLUMEDOC GRIDCALC FIND_GRID_MAXIMUM
      38             : /*
      39             : Find the point with the highest value of the function on the grid
      40             : 
      41             : \par Examples
      42             : 
      43             : */
      44             : //+ENDPLUMEDOC
      45             : 
      46             : namespace PLMD {
      47             : namespace gridtools {
      48             : 
      49             : class FindGridOptimum : public ActionWithGrid {
      50             : private:
      51             :   bool domin;
      52             :   double cgtol;
      53             :   std::unique_ptr<Interpolator> function;
      54             : public:
      55             :   static void registerKeywords( Keywords& keys );
      56             :   explicit FindGridOptimum(const ActionOptions&ao);
      57           0 :   void setupOnFirstStep( const bool incalc ) override { plumed_error(); }
      58             :   unsigned getNumberOfDerivatives() override ;
      59             :   void calculate() override ;
      60             :   std::vector<std::string> getGridCoordinateNames() const override ;
      61             :   const GridCoordinatesObject& getGridCoordinatesObject() const override ;
      62             :   double calculateValueAndDerivatives( const std::vector<double>& pp, std::vector<double>& der );
      63           0 :   void performTask( const unsigned& current, MultiValue& myvals ) const override { plumed_error(); }
      64             : };
      65             : 
      66             : PLUMED_REGISTER_ACTION(FindGridOptimum,"FIND_GRID_MAXIMUM")
      67             : PLUMED_REGISTER_ACTION(FindGridOptimum,"FIND_GRID_MINIMUM")
      68             : 
      69           8 : void FindGridOptimum::registerKeywords( Keywords& keys ) {
      70           8 :   ActionWithGrid::registerKeywords( keys );
      71          16 :   keys.addInputKeyword("compulsory","ARG","grid","the label for the function on the grid that you would like to find the optimum in");
      72          16 :   keys.addFlag("NOINTERPOL",false,"do not interpolate the function when finding the optimum");
      73          16 :   keys.add("compulsory","CGTOL","1E-4","the tolerance for the conjugate gradient optimization");
      74          16 :   keys.addOutputComponent("optval","default","scalar","the value of the function at the optimum");
      75          16 :   keys.addOutputComponent("_opt","default","scalar","the values of the arguments of the function at the optimum can be referenced elsewhere in the input file "
      76             :                           "by using the names of the arguments followed by the string _opt");
      77           8 : }
      78             : 
      79           2 : FindGridOptimum::FindGridOptimum(const ActionOptions&ao):
      80             :   Action(ao),
      81             :   ActionWithGrid(ao),
      82           2 :   cgtol(0)
      83             : {
      84           2 :   if( getName()=="FIND_GRID_MAXIMUM" ) domin=false;
      85           2 :   else if( getName()=="FIND_GRID_MINIMUM" ) domin=true;
      86           0 :   else plumed_error();
      87             :   // Create value for this function
      88           2 :   std::vector<std::string> argn( getGridCoordinateNames() );
      89           5 :   std::vector<unsigned> shape(0); for(unsigned i=0; i<argn.size(); ++i) addComponent( argn[i] + "_opt", shape );
      90           4 :   addComponent( "optval", shape ); componentIsNotPeriodic( "optval" );
      91           2 :   bool nointerpol=false; parseFlag("NOINTERPOL",nointerpol);
      92           6 :   if( !nointerpol ) { parse("CGTOL",cgtol); function=Tools::make_unique<Interpolator>( getPntrToArgument(0), getGridCoordinatesObject() ); }
      93           2 : }
      94             : 
      95           0 : unsigned FindGridOptimum::getNumberOfDerivatives() {
      96           0 :   return 0;
      97             : }
      98             : 
      99         240 : const GridCoordinatesObject& FindGridOptimum::getGridCoordinatesObject() const {
     100         240 :   ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
     101         240 :   plumed_assert( ag ); return ag->getGridCoordinatesObject();
     102             : }
     103             : 
     104           2 : std::vector<std::string> FindGridOptimum::getGridCoordinateNames() const {
     105           2 :   ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
     106           2 :   plumed_assert( ag ); return ag->getGridCoordinateNames();
     107             : }
     108             : 
     109         236 : double FindGridOptimum::calculateValueAndDerivatives( const std::vector<double>& pp, std::vector<double>& der ) {
     110         236 :   double val = function->splineInterpolation( pp, der );
     111             :   // We normalise the derivatives here and set them so that the linesearch is done over the cell that we know
     112             :   // in the grid that we know the minimum is inside
     113         236 :   std::vector<double> spacing( getGridCoordinatesObject().getGridSpacing() );
     114         623 :   double norm = 0; for(unsigned i=0; i<der.size(); ++i) norm += der[i]*der[i];
     115         623 :   norm = sqrt(norm); for(unsigned i=0; i<der.size(); ++i) der[i] = spacing[i]*der[i] / norm;
     116         236 :   if( domin ) return val;
     117             :   // If we are looking for a maximum
     118           0 :   for(unsigned i=0; i<der.size(); ++i) der[i] = -der[i];
     119           0 :   return -val;
     120             : }
     121             : 
     122           2 : void FindGridOptimum::calculate() {
     123           2 :   const GridCoordinatesObject& ingrid = getGridCoordinatesObject(); Value* gval = getPntrToArgument(0);
     124           2 :   std::vector<double> optargs( gval->getRank() ); std::vector<unsigned> gridind( gval->getRank() );
     125           2 :   double optval=gval->get( 0 ); ingrid.getGridPointCoordinates( 0, gridind, optargs );
     126           2 :   unsigned nval = gval->getNumberOfValues(); bool constant=true;
     127        2602 :   for(unsigned i=0; i<nval; ++i) {
     128        2600 :     double tval = gval->get( i );
     129        2600 :     if( domin && (tval<optval || std::isnan(optval)) ) { constant=false; optval=tval; ingrid.getGridPointCoordinates( i, gridind, optargs ); }
     130        2600 :     if( !domin && (tval>optval || std::isnan(optval)) ) { constant=false; optval=tval; ingrid.getGridPointCoordinates( i, gridind, optargs ); }
     131             :   }
     132             :   // This basically ensures we deal with cases where all points on the grid are infinity as isinf doesn't work on intel compiler
     133           2 :   if( constant ) {
     134           0 :     if( domin && gval->get(0)>=gval->get(1) ) return;
     135           0 :     else if( gval->get(0)<=gval->get(1) ) return;
     136             :   }
     137           2 :   if( std::isinf(optval) ) { return; }
     138             : 
     139           2 :   if( std::isnan(optval) ) error("all values on grid are nans");
     140             :   // And do conjugate gradient optimisation (because we can!!)
     141           2 :   if( cgtol>0 ) {
     142             :     ConjugateGradient<FindGridOptimum> myminimiser( this );
     143           2 :     myminimiser.minimise( cgtol, optargs, &FindGridOptimum::calculateValueAndDerivatives );
     144             :   }
     145             :   // And set the final value
     146           5 :   for(unsigned j=0; j<optargs.size(); ++j) getPntrToComponent(j)->set( optargs[j] );
     147           2 :   std::vector<double> optder( gval->getRank() );
     148           2 :   getPntrToComponent(optargs.size())->set( function->splineInterpolation( optargs, optder ) );
     149             : }
     150             : 
     151             : }
     152             : }