Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-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             : #ifndef __PLUMED_tools_Minimise1DBrent_h
      23             : #define __PLUMED_tools_Minimise1DBrent_h
      24             : 
      25             : #include "Tools.h"
      26             : 
      27             : #include <vector>
      28             : #include <string>
      29             : 
      30             : namespace PLMD {
      31             : 
      32             : /// A class for doing parabolic interpolation and minimisation of
      33             : /// 1D functions using Brent's method.
      34             : template <class FCLASS>
      35        1192 : class Minimise1DBrent {
      36             : private:
      37             : /// Has the minimum been bracketed
      38             :   bool bracketed;
      39             : /// Has the function been minimised
      40             :   bool minimised;
      41             : /// The tolerance for the line minimiser
      42             :   double tol;
      43             : /// The default ration by which successive intervals are magnified
      44             :   const double GOLD;
      45             : /// The maximum magnification allowed for a parabolic fit step
      46             :   const double GLIMIT;
      47             : /// Use to prevent any possible division by zero
      48             :   const double TINY;
      49             : /// Maximum number of interactions in line minimiser
      50             :   const unsigned ITMAX;
      51             : /// The value of the golden ratio
      52             :   const double CGOLD;
      53             : /// A small number that protects against trying to achieve fractional
      54             : /// accuracy for a minimum that happens to be exactly zero
      55             :   const double ZEPS;
      56             : /// This is the type specifier for the function to minimise
      57             :   typedef double(FCLASS::*eng_pointer)( const double& val );
      58             : /// Three points bracketting the minimum and the corresponding function values
      59             :   double ax,bx,cx,fa,fb,fc,fmin;
      60             : /// The class containing the function we are trying to minimise
      61             :   FCLASS myclass_func;
      62             : public:
      63        1192 :   explicit Minimise1DBrent( const FCLASS& pf,  const double& t=3.0E-8 );
      64             : /// Bracket the minium
      65             :   void bracket( const double& ax, const double& xx, eng_pointer eng );
      66             : /// Find the minimum between two brackets
      67             :   double minimise( eng_pointer eng );
      68             : /// Return the value of the function at the minimum
      69             :   double getMinimumValue() const ;
      70             : };
      71             : 
      72             : template <class FCLASS>
      73        1192 : Minimise1DBrent<FCLASS>::Minimise1DBrent( const FCLASS& pf, const double& t ):
      74        1192 :   bracketed(false),
      75        1192 :   minimised(false),
      76        1192 :   tol(t),
      77        1192 :   GOLD(1.618034),
      78        1192 :   GLIMIT(100.0),
      79        1192 :   TINY(1.0E-20),
      80        1192 :   ITMAX(100),
      81        1192 :   CGOLD(0.3819660),
      82        1192 :   ZEPS(epsilon*1.0E-3),
      83        1192 :   ax(0),bx(0),cx(0),
      84        1192 :   fa(0),fb(0),fc(0),
      85        1192 :   fmin(0),
      86        1192 :   myclass_func(pf) {
      87        1192 : }
      88             : 
      89             : template <class FCLASS>
      90        1192 : void Minimise1DBrent<FCLASS>::bracket( const double& a, const double& b, eng_pointer eng ) {
      91        1192 :   ax=a;
      92        1192 :   bx=b;
      93             :   double fu;
      94        1192 :   fa=(myclass_func.*eng)(ax);
      95        1192 :   fb=(myclass_func.*eng)(bx);
      96        1192 :   if( fb>fa ) {
      97             :     double tmp;
      98         155 :     tmp=ax;
      99         155 :     ax=bx;
     100         155 :     bx=tmp;
     101             :     tmp=fa;
     102         155 :     fa=fb;
     103         155 :     fb=tmp;
     104             :   }
     105        1192 :   cx=bx+GOLD*(bx-ax);
     106        1192 :   fc=(myclass_func.*eng)(cx);
     107        2638 :   while ( fb > fc ) {
     108        1521 :     double r=(bx-ax)*(fb-fc);
     109        1521 :     double q=(bx-cx)*(fb-fa);
     110        1521 :     double u=bx-((bx-cx)*q-(bx-ax)*r)/(2.0*(std::fabs(q-r)>TINY?std::fabs(q-r):TINY)*(q-r>=0?1:-1));
     111        1521 :     double ulim=bx+GLIMIT*(cx-bx);
     112        1521 :     if((bx-u)*(u-cx) > 0.0 ) {
     113         151 :       fu=(myclass_func.*eng)(u);
     114         151 :       if( fu < fc ) {
     115          67 :         ax=bx;
     116          67 :         bx=u;
     117          67 :         fa=fb;
     118          67 :         fb=fu;
     119          67 :         bracketed=true;
     120          75 :         return;
     121          84 :       } else if( fu > fb ) {
     122           8 :         cx=u;
     123           8 :         fc=fu;
     124           8 :         bracketed=true;
     125           8 :         return;
     126             :       }
     127          76 :       u=cx+GOLD*(cx-bx);
     128          76 :       fu=(myclass_func.*eng)(u);
     129        1370 :     } else if((cx-u)*(u-ulim) > 0.0 ) {
     130        1052 :       fu=(myclass_func.*eng)(u);
     131        1052 :       if( fu<fc ) {
     132        1041 :         bx=cx;
     133        1041 :         cx=u;
     134        1041 :         u+=GOLD*(u-bx);
     135        1041 :         fb=fc;
     136        1041 :         fc=fu;
     137        1041 :         fu=(myclass_func.*eng)(u);
     138             :       }
     139         318 :     } else if( (u-ulim)*(ulim-cx) >= 0.0 ) {
     140         283 :       u=ulim;
     141         283 :       fu=(myclass_func.*eng)(u);
     142             :     } else {
     143          35 :       u=cx+GOLD*(cx-bx);
     144          35 :       fu=(myclass_func.*eng)(u);
     145             :     }
     146        1446 :     ax=bx;
     147        1446 :     bx=cx;
     148        1446 :     cx=u;
     149        1446 :     fa=fb;
     150        1446 :     fb=fc;
     151        1446 :     fc=fu;
     152             :   }
     153        1117 :   bracketed=true;
     154             : }
     155             : 
     156             : template <class FCLASS>
     157        1192 : double Minimise1DBrent<FCLASS>::minimise( eng_pointer eng ) {
     158             :   plumed_dbg_assert( bracketed );
     159             : 
     160             :   double a,b,d=0.0,etemp,fu,fv,fw,fx;
     161             :   double p,q,r,tol1,tol2,u,v,w,x,xm;
     162             :   double e=0.0;
     163             : 
     164        1192 :   a=(ax < cx ? ax : cx );
     165        1192 :   b=(ax >= cx ? ax : cx );
     166        1192 :   x=w=v=bx;
     167        1192 :   fw=fv=fx=(myclass_func.*eng)(x);
     168       18942 :   for(unsigned iter=0; iter<ITMAX; ++iter) {
     169       18942 :     xm=0.5*(a+b);
     170       18942 :     tol2=2.0*(tol1=tol*std::fabs(x)+ZEPS);
     171       18942 :     if( std::fabs(x-xm) <= (tol2-0.5*(b-a))) {
     172        1192 :       fmin=fx;
     173        1192 :       minimised=true;
     174        1192 :       return x;
     175             :     }
     176       17750 :     if( std::fabs(e) > tol1 ) {
     177       16377 :       r=(x-w)*(fx-fv);
     178       16377 :       q=(x-v)*(fx-fw);
     179       16377 :       p=(x-v)*q-(x-w)*r;
     180       16377 :       q=2.0*(q-r);
     181       16377 :       if( q > 0.0 ) {
     182        7818 :         p = -p;
     183             :       }
     184       16377 :       q=std::fabs(q);
     185             :       etemp=e;
     186             :       e=d;
     187       16377 :       if( std::fabs(p) >= std::fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x) ) {
     188        7607 :         d = CGOLD*(e=(x >= xm ? a-x : b-x ));
     189             :       } else {
     190        8770 :         d=p/q;
     191        8770 :         u=x+d;
     192        8770 :         if(u-a < tol2 || b-u < tol2 ) {
     193         705 :           d=(xm-x>=0?std::fabs(tol1):-std::fabs(tol1));
     194             :         }
     195             :       }
     196             :     } else {
     197        1373 :       d=CGOLD*(e=( x >= xm ? a-x : b-x ));
     198             :     }
     199       17750 :     if( std::fabs(d)>=tol1) {
     200       16896 :       u=x+d;
     201             :     } else {
     202         854 :       u=x+(d>=0?std::fabs(tol1):-std::fabs(tol1));
     203             :     }
     204       17750 :     fu=(myclass_func.*eng)(u);
     205       17750 :     if( fu <= fx ) {
     206        5911 :       if( u >= x ) {
     207             :         a=x;
     208             :       } else {
     209             :         b=x;
     210             :       }
     211             :       v=w;
     212             :       fv=fw;
     213             :       w=x;
     214             :       fw=fx;
     215        5911 :       x=u;
     216             :       fx=fu;
     217             :     } else {
     218       11839 :       if( u < x ) {
     219             :         a=u;
     220             :       } else {
     221             :         b=u;
     222             :       }
     223       11839 :       if( fu <=fw || w==x ) {
     224             :         v=w;
     225             :         w=u;
     226             :         fv=fw;
     227             :         fw=fu;
     228        7526 :       } else if( fu <= fv || v==x || v==w ) {
     229             :         v=u;
     230             :         fv=fu;
     231             :       }
     232             :     }
     233             :   }
     234           0 :   plumed_merror("Too many interactions in brent");
     235             : }
     236             : 
     237             : template <class FCLASS>
     238        1192 : double Minimise1DBrent<FCLASS>::getMinimumValue() const {
     239             :   plumed_dbg_assert( minimised );
     240        1192 :   return fmin;
     241             : }
     242             : 
     243             : }
     244             : #endif