Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2019 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_Brent1DRootSearch_h
      23             : #define __PLUMED_tools_Brent1DRootSearch_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        1516 : class Brent1DRootSearch {
      36             : private:
      37             : /// Has the minimum been bracketed
      38             :   bool bracketed;
      39             : /// The tolerance for the line minimiser
      40             :   double tol;
      41             : /// Maximum number of interactions in line minimiser
      42             :   const unsigned ITMAX;
      43             : /// A small number that protects against trying to achieve fractional
      44             : /// accuracy for a minimum that happens to be exactly zero
      45             :   const double EPS;
      46             : /// The factor by which to expand the range when bracketing
      47             :   const double EXPAND;
      48             : /// This is the type specifier for the function to minimise
      49             :   typedef double(FCLASS::*eng_pointer)( const double& val );
      50             : /// Three points bracketting the minimum and the corresponding function values
      51             :   double ax,bx,fa,fb;
      52             : /// The class containing the function we are trying to minimise
      53             :   FCLASS myclass_func;
      54             : public:
      55             :   explicit Brent1DRootSearch( const FCLASS& pf,  const double& t=3.0E-8 );
      56             : /// Bracket the minium
      57             :   void bracket( const double& ax, const double& xx, eng_pointer eng );
      58             : /// Find the minimum between two brackets
      59             :   double search( eng_pointer eng );
      60             : };
      61             : 
      62             : template <class FCLASS>
      63        1516 : Brent1DRootSearch<FCLASS>::Brent1DRootSearch( const FCLASS& pf, const double& t ):
      64             :   bracketed(false),
      65             :   tol(t),
      66             :   ITMAX(100),
      67             :   EPS(3.0E-8),
      68             :   EXPAND(1.6),
      69             :   ax(0), bx(0),
      70             :   fa(0), fb(0),
      71        1516 :   myclass_func(pf)
      72             : {
      73        1516 : }
      74             : 
      75             : template <class FCLASS>
      76        1516 : void Brent1DRootSearch<FCLASS>::bracket( const double& a, const double& b, eng_pointer eng ) {
      77        1516 :   plumed_assert( a!=b ); ax=a; bx=b; fa=(myclass_func.*eng)(a); fb=(myclass_func.*eng)(b);
      78        1516 :   if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) plumed_merror("input points do not bracket root");
      79        1516 :   bracketed=true;
      80        1516 : }
      81             : 
      82             : template <class FCLASS>
      83        1516 : double Brent1DRootSearch<FCLASS>::search( eng_pointer eng ) {
      84             :   plumed_dbg_assert( bracketed );
      85             : 
      86        1516 :   double cx=bx, d, e, min1, min2, fc=fb, p, q, r, s, tol1, xm;
      87       16344 :   for(unsigned iter=0; iter<ITMAX; iter++) {
      88        8930 :     if ( (fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0) ) { cx=ax; fc=fa; e=d=bx-ax; }
      89        8930 :     if( fabs(fc) < fabs(fb) ) { ax=bx; bx=cx; cx=ax; fa=fb; fb=fc; fc=fa; }
      90        8930 :     tol1=2*EPS*fabs(bx)+0.5*tol; xm=0.5*(cx-bx);
      91        8930 :     if( fabs(xm) <= tol1 || fb == 0.0 ) return bx;
      92        7414 :     if( fabs(e) >= tol1 && fabs(fa) > fabs(fb) ) {
      93        7350 :       s=fb/fa;
      94        7350 :       if( ax==cx ) {
      95        5192 :         p=2.0*xm*s; q=1.0-s;
      96             :       } else {
      97        2158 :         q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(bx-ax)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0);
      98             :       }
      99        7350 :       if (p > 0.0) q = -q;
     100        7350 :       p=fabs(p); min1=3.0*xm*q-fabs(tol1*q); min2=fabs(e*q);
     101        7350 :       if (2.0*p < (min1 < min2 ? min1 : min2)) {
     102        7091 :         e=d; d=p/q;
     103             :       } else {
     104             :         d=xm; e=d;
     105             :       }
     106             :     } else {
     107             :       d=xm; e=d;
     108             :     }
     109        7414 :     ax=bx; fa=fb;
     110        7414 :     if( fabs(d) > tol1 ) bx+=d;
     111        1427 :     else if(xm<0 ) bx -= fabs(tol1); // SIGN(tol1,xm);
     112         698 :     else bx += tol1;
     113        7414 :     fb = (myclass_func.*eng)(bx);
     114             :   }
     115             : 
     116           0 :   plumed_merror("Too many interactions in zbrent");
     117             : }
     118             : 
     119             : }
     120             : #endif