Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-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_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        2187 : 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        2187 :   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        2187 : Brent1DRootSearch<FCLASS>::Brent1DRootSearch( const FCLASS& pf, const double& t ):
      64        2187 :   bracketed(false),
      65        2187 :   tol(t),
      66        2187 :   ITMAX(100),
      67        2187 :   EPS(3.0E-8),
      68        2187 :   EXPAND(1.6),
      69        2187 :   ax(0), bx(0),
      70        2187 :   fa(0), fb(0),
      71        2187 :   myclass_func(pf) {
      72        2187 : }
      73             : 
      74             : template <class FCLASS>
      75        2187 : void Brent1DRootSearch<FCLASS>::bracket( const double& a, const double& b, eng_pointer eng ) {
      76        2187 :   plumed_assert( a!=b );
      77        2187 :   ax=a;
      78        2187 :   bx=b;
      79        2187 :   fa=(myclass_func.*eng)(a);
      80        2187 :   fb=(myclass_func.*eng)(b);
      81        2187 :   if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
      82           0 :     plumed_merror("input points do not bracket root");
      83             :   }
      84        2187 :   bracketed=true;
      85        2187 : }
      86             : 
      87             : template <class FCLASS>
      88        2187 : double Brent1DRootSearch<FCLASS>::search( eng_pointer eng ) {
      89             :   plumed_dbg_assert( bracketed );
      90             : 
      91        2187 :   double cx=bx, d, e, min1, min2, fc=fb, p, q, r, s, tol1, xm;
      92       13534 :   for(unsigned iter=0; iter<ITMAX; iter++) {
      93       13534 :     if ( (fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0) ) {
      94        9920 :       cx=ax;
      95        9920 :       fc=fa;
      96        9920 :       e=d=bx-ax;
      97             :     }
      98       13534 :     if( std::fabs(fc) < std::fabs(fb) ) {
      99        3796 :       ax=bx;
     100        3796 :       bx=cx;
     101             :       cx=ax;
     102        3796 :       fa=fb;
     103        3796 :       fb=fc;
     104             :       fc=fa;
     105             :     }
     106       13534 :     tol1=2*EPS*std::fabs(bx)+0.5*tol;
     107       13534 :     xm=0.5*(cx-bx);
     108       13534 :     if( std::fabs(xm) <= tol1 || fb == 0.0 ) {
     109        2187 :       return bx;
     110             :     }
     111       11347 :     if( std::fabs(e) >= tol1 && std::fabs(fa) > std::fabs(fb) ) {
     112       11219 :       s=fb/fa;
     113       11219 :       if( ax==cx ) {
     114        7754 :         p=2.0*xm*s;
     115        7754 :         q=1.0-s;
     116             :       } else {
     117        3465 :         q=fa/fc;
     118        3465 :         r=fb/fc;
     119        3465 :         p=s*(2.0*xm*q*(q-r)-(bx-ax)*(r-1.0));
     120        3465 :         q=(q-1.0)*(r-1.0)*(s-1.0);
     121             :       }
     122       11219 :       if (p > 0.0) {
     123        5741 :         q = -q;
     124             :       }
     125       11219 :       p=std::fabs(p);
     126       11219 :       min1=3.0*xm*q-std::fabs(tol1*q);
     127       11219 :       min2=std::fabs(e*q);
     128       21784 :       if (2.0*p < (min1 < min2 ? min1 : min2)) {
     129             :         e=d;
     130       10537 :         d=p/q;
     131             :       } else {
     132             :         d=xm;
     133             :         e=d;
     134             :       }
     135             :     } else {
     136             :       d=xm;
     137             :       e=d;
     138             :     }
     139       11347 :     ax=bx;
     140       11347 :     fa=fb;
     141       11347 :     if( std::fabs(d) > tol1 ) {
     142        9300 :       bx+=d;
     143        2047 :     } else if(xm<0 ) {
     144        1016 :       bx -= std::fabs(tol1);  // SIGN(tol1,xm);
     145             :     } else {
     146        1031 :       bx += tol1;
     147             :     }
     148       11347 :     fb = (myclass_func.*eng)(bx);
     149             :   }
     150             : 
     151           0 :   plumed_merror("Too many interactions in zbrent");
     152             : }
     153             : 
     154             : }
     155             : #endif