1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2016-2021 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 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>
Brent1DRootSearch(const FCLASS & pf,const double & t)63 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 myclass_func(pf)
72 {
73 }
74
75 template <class FCLASS>
bracket(const double & a,const double & b,eng_pointer eng)76 void Brent1DRootSearch<FCLASS>::bracket( const double& a, const double& b, eng_pointer eng ) {
77 plumed_assert( a!=b ); ax=a; bx=b; fa=(myclass_func.*eng)(a); fb=(myclass_func.*eng)(b);
78 if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) plumed_merror("input points do not bracket root");
79 bracketed=true;
80 }
81
82 template <class FCLASS>
search(eng_pointer eng)83 double Brent1DRootSearch<FCLASS>::search( eng_pointer eng ) {
84 plumed_dbg_assert( bracketed );
85
86 double cx=bx, d, e, min1, min2, fc=fb, p, q, r, s, tol1, xm;
87 for(unsigned iter=0; iter<ITMAX; iter++) {
88 if ( (fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0) ) { cx=ax; fc=fa; e=d=bx-ax; }
89 if( fabs(fc) < fabs(fb) ) { ax=bx; bx=cx; cx=ax; fa=fb; fb=fc; fc=fa; }
90 tol1=2*EPS*fabs(bx)+0.5*tol; xm=0.5*(cx-bx);
91 if( fabs(xm) <= tol1 || fb == 0.0 ) return bx;
92 if( fabs(e) >= tol1 && fabs(fa) > fabs(fb) ) {
93 s=fb/fa;
94 if( ax==cx ) {
95 p=2.0*xm*s; q=1.0-s;
96 } else {
97 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 if (p > 0.0) q = -q;
100 p=fabs(p); min1=3.0*xm*q-fabs(tol1*q); min2=fabs(e*q);
101 if (2.0*p < (min1 < min2 ? min1 : min2)) {
102 e=d; d=p/q;
103 } else {
104 d=xm; e=d;
105 }
106 } else {
107 d=xm; e=d;
108 }
109 ax=bx; fa=fb;
110 if( fabs(d) > tol1 ) bx+=d;
111 else if(xm<0 ) bx -= fabs(tol1); // SIGN(tol1,xm);
112 else bx += tol1;
113 fb = (myclass_func.*eng)(bx);
114 }
115
116 plumed_merror("Too many interactions in zbrent");
117 }
118
119 }
120 #endif
121