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=false;
      39             : /// The tolerance for the line minimiser
      40             :   double tol;
      41             : /// Maximum number of interactions in line minimiser
      42             :   static constexpr unsigned ITMAX=100;
      43             : /// A small number that protects against trying to achieve fractional
      44             : /// accuracy for a minimum that happens to be exactly zero
      45             :   static constexpr double EPS=3.0e-8;
      46             : /// The factor by which to expand the range when bracketing
      47             :   static constexpr double EXPAND=1.6;
      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=0.0, bx=0.0, fa=0.0, fb=0.0;
      52             : /// The class containing the function we are trying to minimise
      53             :   FCLASS myclass_func;
      54             : public:
      55             :   explicit Brent1DRootSearch( const FCLASS& pf,  double t=3.0E-8 );
      56             : /// Bracket the minium
      57             :   void bracket( double ax, 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 :   tol(t),
      65        2187 :   myclass_func(pf) {
      66        2187 : }
      67             : 
      68             : template <class FCLASS>
      69        2187 : void Brent1DRootSearch<FCLASS>::bracket( const double a, const double b,  eng_pointer eng ) {
      70        2187 :   plumed_assert( a!=b );
      71        2187 :   ax=a;
      72        2187 :   bx=b;
      73        2187 :   fa=(myclass_func.*eng)(a);
      74        2187 :   fb=(myclass_func.*eng)(b);
      75        2187 :   if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
      76           0 :     plumed_merror("input points do not bracket root");
      77             :   }
      78        2187 :   bracketed=true;
      79        2187 : }
      80             : 
      81             : template <class FCLASS>
      82        2187 : double Brent1DRootSearch<FCLASS>::search( eng_pointer eng ) {
      83             :   plumed_dbg_assert( bracketed );
      84             :   // starting with these parameters:
      85             :   // double cx=bx;
      86             :   // double fc=fb;
      87             :   //by definition this is true :
      88             :   // if ( (fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0) ) {
      89             :   // so we initialize the variable/registers
      90             :   // with the body of the first if statement
      91        2187 :   double cx=ax;
      92        2187 :   double fc=fa;
      93        2187 :   double d = bx -ax;
      94             :   double e=d;
      95             :   double min1, min2, p, q, r, s, tol1, xm;
      96       13540 :   for(unsigned iter=0; iter<ITMAX; iter++) {
      97       13540 :     if ( (fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0) ) {
      98        7745 :       cx=ax;
      99        7745 :       fc=fa;
     100        7745 :       e=d=bx-ax;
     101             :     }
     102       13540 :     if( std::fabs(fc) < std::fabs(fb) ) {
     103        3802 :       ax=bx;
     104        3802 :       bx=cx;
     105             :       cx=ax;
     106        3802 :       fa=fb;
     107        3802 :       fb=fc;
     108             :       fc=fa;
     109             :     }
     110       13540 :     tol1=2*EPS*std::fabs(bx)+0.5*tol;
     111       13540 :     xm=0.5*(cx-bx);
     112       13540 :     if( std::fabs(xm) <= tol1 || fb == 0.0 ) {
     113        2187 :       return bx;
     114             :     }
     115       11353 :     if( std::fabs(e) >= tol1 && std::fabs(fa) > std::fabs(fb) ) {
     116       11225 :       s=fb/fa;
     117       11225 :       if( ax==cx ) {
     118        7760 :         p=2.0*xm*s;
     119        7760 :         q=1.0-s;
     120             :       } else {
     121        3465 :         q=fa/fc;
     122        3465 :         r=fb/fc;
     123        3465 :         p=s*(2.0*xm*q*(q-r)-(bx-ax)*(r-1.0));
     124        3465 :         q=(q-1.0)*(r-1.0)*(s-1.0);
     125             :       }
     126       11225 :       if (p > 0.0) {
     127        5744 :         q = -q;
     128             :       }
     129       11225 :       p=std::fabs(p);
     130       11225 :       min1=3.0*xm*q-std::fabs(tol1*q);
     131       11225 :       min2=std::fabs(e*q);
     132       21796 :       if (2.0*p < (min1 < min2 ? min1 : min2)) {
     133             :         e=d;
     134       10543 :         d=p/q;
     135             :       } else {
     136             :         d=xm;
     137             :         e=d;
     138             :       }
     139             :     } else {
     140             :       d=xm;
     141             :       e=d;
     142             :     }
     143       11353 :     ax=bx;
     144       11353 :     fa=fb;
     145       11353 :     if( std::fabs(d) > tol1 ) {
     146        9300 :       bx+=d;
     147        2053 :     } else if(xm<0 ) {
     148        1013 :       bx -= std::fabs(tol1);  // SIGN(tol1,xm);
     149             :     } else {
     150        1040 :       bx += tol1;
     151             :     }
     152       11353 :     fb = (myclass_func.*eng)(bx);
     153             :   }
     154             : 
     155           0 :   plumed_merror("Too many interactions in zbrent");
     156             : }
     157             : 
     158             : }
     159             : #endif