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 897 : 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 897 : 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 897 : Minimise1DBrent<FCLASS>::Minimise1DBrent( const FCLASS& pf, const double& t ):
74 897 : bracketed(false),
75 897 : minimised(false),
76 897 : tol(t),
77 897 : GOLD(1.618034),
78 897 : GLIMIT(100.0),
79 897 : TINY(1.0E-20),
80 897 : ITMAX(100),
81 897 : CGOLD(0.3819660),
82 897 : ZEPS(epsilon*1.0E-3),
83 897 : ax(0),bx(0),cx(0),
84 897 : fa(0),fb(0),fc(0),
85 897 : fmin(0),
86 897 : myclass_func(pf)
87 : {
88 897 : }
89 :
90 : template <class FCLASS>
91 897 : void Minimise1DBrent<FCLASS>::bracket( const double& a, const double& b, eng_pointer eng ) {
92 897 : ax=a; bx=b; double fu;
93 897 : fa=(myclass_func.*eng)(ax); fb=(myclass_func.*eng)(bx);
94 897 : if( fb>fa ) {
95 : double tmp;
96 809 : tmp=ax; ax=bx; bx=tmp;
97 809 : tmp=fa; fa=fb; fb=tmp;
98 : }
99 897 : cx=bx+GOLD*(bx-ax);
100 897 : fc=(myclass_func.*eng)(cx);
101 1113 : while ( fb > fc ) {
102 220 : double r=(bx-ax)*(fb-fc);
103 220 : double q=(bx-cx)*(fb-fa);
104 220 : 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));
105 220 : double ulim=bx+GLIMIT*(cx-bx);
106 220 : if((bx-u)*(u-cx) > 0.0 ) {
107 5 : fu=(myclass_func.*eng)(u);
108 5 : if( fu < fc ) {
109 4 : ax=bx; bx=u; fa=fb; fb=fu; bracketed=true; return;
110 1 : } else if( fu > fb ) {
111 0 : cx=u; fc=fu; bracketed=true; return;
112 : }
113 1 : u=cx+GOLD*(cx-bx); fu=(myclass_func.*eng)(u);
114 215 : } else if((cx-u)*(u-ulim) > 0.0 ) {
115 81 : fu=(myclass_func.*eng)(u);
116 81 : if( fu<fc ) {
117 80 : bx=cx; cx=u; u+=GOLD*(u-bx);
118 80 : fb=fc; fc=fu; fu=(myclass_func.*eng)(u);
119 : }
120 134 : } else if( (u-ulim)*(ulim-cx) >= 0.0 ) {
121 124 : u=ulim;
122 124 : fu=(myclass_func.*eng)(u);
123 : } else {
124 10 : u=cx+GOLD*(cx-bx);
125 10 : fu=(myclass_func.*eng)(u);
126 : }
127 216 : ax=bx; bx=cx; cx=u;
128 216 : fa=fb; fb=fc; fc=fu;
129 : }
130 893 : bracketed=true;
131 : }
132 :
133 : template <class FCLASS>
134 897 : double Minimise1DBrent<FCLASS>::minimise( eng_pointer eng ) {
135 : plumed_dbg_assert( bracketed );
136 :
137 : double a,b,d=0.0,etemp,fu,fv,fw,fx;
138 : double p,q,r,tol1,tol2,u,v,w,x,xm;
139 : double e=0.0;
140 :
141 897 : a=(ax < cx ? ax : cx );
142 897 : b=(ax >= cx ? ax : cx );
143 897 : x=w=v=bx;
144 897 : fw=fv=fx=(myclass_func.*eng)(x);
145 14768 : for(unsigned iter=0; iter<ITMAX; ++iter) {
146 14768 : xm=0.5*(a+b);
147 14768 : tol2=2.0*(tol1=tol*std::fabs(x)+ZEPS);
148 14768 : if( std::fabs(x-xm) <= (tol2-0.5*(b-a))) {
149 897 : fmin=fx; minimised=true; return x;
150 : }
151 13871 : if( std::fabs(e) > tol1 ) {
152 12810 : r=(x-w)*(fx-fv);
153 12810 : q=(x-v)*(fx-fw);
154 12810 : p=(x-v)*q-(x-w)*r;
155 12810 : q=2.0*(q-r);
156 12810 : if( q > 0.0 ) p = -p;
157 12810 : q=std::fabs(q);
158 : etemp=e;
159 : e=d;
160 12810 : if( std::fabs(p) >= std::fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x) ) {
161 5079 : d = CGOLD*(e=(x >= xm ? a-x : b-x ));
162 : } else {
163 7731 : d=p/q; u=x+d;
164 7731 : if(u-a < tol2 || b-u < tol2 ) d=(xm-x>=0?std::fabs(tol1):-std::fabs(tol1));
165 : }
166 : } else {
167 1061 : d=CGOLD*(e=( x >= xm ? a-x : b-x ));
168 : }
169 13871 : if( std::fabs(d)>=tol1) u=x+d; else u=x+(d>=0?std::fabs(tol1):-std::fabs(tol1));
170 13871 : fu=(myclass_func.*eng)(u);
171 13871 : if( fu <= fx ) {
172 5411 : if( u >= x ) a=x; else b=x;
173 : v=w; fv=fw;
174 : w=x; fw=fx;
175 5411 : x=u; fx=fu;
176 : } else {
177 8460 : if( u < x ) a=u; else b=u;
178 8460 : if( fu <=fw || w==x ) {
179 : v=w; w=u; fv=fw; fw=fu;
180 5096 : } else if( fu <= fv || v==x || v==w ) {
181 : v=u; fv=fu;
182 : }
183 : }
184 : }
185 0 : plumed_merror("Too many interactions in brent");
186 : }
187 :
188 : template <class FCLASS>
189 897 : double Minimise1DBrent<FCLASS>::getMinimumValue() const {
190 : plumed_dbg_assert( minimised );
191 897 : return fmin;
192 : }
193 :
194 : }
195 : #endif
|
