1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
|
#ifndef BRENT_HPP
#define BRENT_HPP
#include <iostream>
#include "bas_util.hpp"
//#include "optimizer.hpp"
namespace opt_utilities
{
template<typename T>
T brent(T ax,T bx,T cx,func_obj<T,T>& f,T tol,T& xmin)
{
const int ITMAX=100;
const T CGOLD=0.3819660;
const T ZEPS=std::numeric_limits<T>::epsilon()*1.e-3;
int iter;
T a=0,b=0,d(0),etemp=0,fu=0,fv=0,fw=0,fx=0,p=0,q=0
,r=0,tol1=0,tol2=0,u=0,v=0,w=0,x=0,xm=0;
T e=0.;
a=(ax<cx?ax:cx);
b=(ax>cx?ax:cx);
x=w=v=bx;
fw=fv=fx=f.eval(x);
for(iter=0;iter<ITMAX;++iter)
{
xm=.5*(a+b);
tol2=2.*(tol1=tol*tabs(x)+ZEPS);
if(tabs(T(x-xm))<=(tol2-.5*(b-a)))
{
xmin=x;
return fx;
}
if(tabs(e)>tol1)
{
r=(x-w)*(fx-fv);
q=(x-v)*(fx-fw);
p=(x-v)*q-(x-w)*r;
q=2.*(q-r);
if(q>0.)
{
p=-p;
}
q=tabs(q);
etemp=e;
e=d;
if(tabs(p)>=tabs(T(T(.5)*p*etemp))||p<=q*(a-x)||p>=q*(b-x))
{
d=CGOLD*(e=(x>=xm?a-x:b-x));
}
else
{
d=p/q;
u=x+d;
if(u-a<tol2||b-u<tol2)
{
d=sign(tol1,T(xm-x));
}
}
}
else
{
d=CGOLD*(e=(x>=xm?a-x:b-x));
}
u=(tabs(d)>=tol1?x+d:x+sign(tol1,d));
fu=f.eval(u);
if(fu<=fx)
{
if(u>=x)
{
a=x;
}
else
{
b=x;
}
shft3(v,w,x,u);
shft3(fv,fw,fx,fu);
}
else
{
if(u<x)
{
a=u;
}
else
{
b=u;
}
if(fu<=fw||w==x)
{
v=w;
w=u;
fv=fw;
fw=fu;
}
else if(fu<=fv||v==x||v==w)
{
v=u;
fv=fu;
}
}
}
std::cerr<<"Too many iterations in brent"<<std::endl;
xmin=x;
return fx;
}
}
#endif
|
