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 ); keys.use("ARG"); 71 16 : keys.addFlag("NOINTERPOL",false,"do not interpolate the function when finding the optimum"); 72 16 : keys.add("compulsory","CGTOL","1E-4","the tolerance for the conjugate gradient optimization"); 73 16 : keys.addOutputComponent("optval","default","the value of the function at the optimum"); 74 16 : keys.addOutputComponent("_opt","default","the values of the arguments of the function at the optimum can be referenced elsewhere in the input file " 75 : "by using the names of the arguments followed by the string _opt"); 76 8 : } 77 : 78 2 : FindGridOptimum::FindGridOptimum(const ActionOptions&ao): 79 : Action(ao), 80 : ActionWithGrid(ao), 81 2 : cgtol(0) 82 : { 83 2 : if( getName()=="FIND_GRID_MAXIMUM" ) domin=false; 84 2 : else if( getName()=="FIND_GRID_MINIMUM" ) domin=true; 85 0 : else plumed_error(); 86 : // Create value for this function 87 2 : std::vector<std::string> argn( getGridCoordinateNames() ); 88 5 : std::vector<unsigned> shape(0); for(unsigned i=0; i<argn.size(); ++i) addComponent( argn[i] + "_opt", shape ); 89 4 : addComponent( "optval", shape ); componentIsNotPeriodic( "optval" ); 90 2 : bool nointerpol=false; parseFlag("NOINTERPOL",nointerpol); 91 6 : if( !nointerpol ) { parse("CGTOL",cgtol); function=Tools::make_unique<Interpolator>( getPntrToArgument(0), getGridCoordinatesObject() ); } 92 2 : } 93 : 94 0 : unsigned FindGridOptimum::getNumberOfDerivatives() { 95 0 : return 0; 96 : } 97 : 98 240 : const GridCoordinatesObject& FindGridOptimum::getGridCoordinatesObject() const { 99 240 : ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() ); 100 240 : plumed_assert( ag ); return ag->getGridCoordinatesObject(); 101 : } 102 : 103 2 : std::vector<std::string> FindGridOptimum::getGridCoordinateNames() const { 104 2 : ActionWithGrid* ag=ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() ); 105 2 : plumed_assert( ag ); return ag->getGridCoordinateNames(); 106 : } 107 : 108 236 : double FindGridOptimum::calculateValueAndDerivatives( const std::vector<double>& pp, std::vector<double>& der ) { 109 236 : double val = function->splineInterpolation( pp, der ); 110 : // We normalise the derivatives here and set them so that the linesearch is done over the cell that we know 111 : // in the grid that we know the minimum is inside 112 236 : std::vector<double> spacing( getGridCoordinatesObject().getGridSpacing() ); 113 623 : double norm = 0; for(unsigned i=0; i<der.size(); ++i) norm += der[i]*der[i]; 114 623 : norm = sqrt(norm); for(unsigned i=0; i<der.size(); ++i) der[i] = spacing[i]*der[i] / norm; 115 236 : if( domin ) return val; 116 : // If we are looking for a maximum 117 0 : for(unsigned i=0; i<der.size(); ++i) der[i] = -der[i]; 118 0 : return -val; 119 : } 120 : 121 2 : void FindGridOptimum::calculate() { 122 2 : const GridCoordinatesObject& ingrid = getGridCoordinatesObject(); Value* gval = getPntrToArgument(0); 123 2 : std::vector<double> optargs( gval->getRank() ); std::vector<unsigned> gridind( gval->getRank() ); 124 2 : double optval=gval->get( 0 ); ingrid.getGridPointCoordinates( 0, gridind, optargs ); 125 2 : unsigned nval = gval->getNumberOfValues(); bool constant=true; 126 2602 : for(unsigned i=0; i<nval; ++i) { 127 2600 : double tval = gval->get( i ); 128 2600 : if( domin && (tval<optval || std::isnan(optval)) ) { constant=false; optval=tval; ingrid.getGridPointCoordinates( i, gridind, optargs ); } 129 2600 : if( !domin && (tval>optval || std::isnan(optval)) ) { constant=false; optval=tval; ingrid.getGridPointCoordinates( i, gridind, optargs ); } 130 : } 131 : // This basically ensures we deal with cases where all points on the grid are infinity as isinf doesn't work on intel compiler 132 2 : if( constant ) { 133 0 : if( domin && gval->get(0)>=gval->get(1) ) return; 134 0 : else if( gval->get(0)<=gval->get(1) ) return; 135 : } 136 2 : if( std::isinf(optval) ) { return; } 137 : 138 2 : if( std::isnan(optval) ) error("all values on grid are nans"); 139 : // And do conjugate gradient optimisation (because we can!!) 140 2 : if( cgtol>0 ) { 141 : ConjugateGradient<FindGridOptimum> myminimiser( this ); 142 2 : myminimiser.minimise( cgtol, optargs, &FindGridOptimum::calculateValueAndDerivatives ); 143 : } 144 : // And set the final value 145 5 : for(unsigned j=0; j<optargs.size(); ++j) getPntrToComponent(j)->set( optargs[j] ); 146 2 : std::vector<double> optder( gval->getRank() ); 147 2 : getPntrToComponent(optargs.size())->set( function->splineInterpolation( optargs, optder ) ); 148 : } 149 : 150 : } 151 : }