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
|
#ifndef MNBRAK_HPP
#define MNBRAK_HPP
#define OPT_HEADER
//#include "optimizer.hpp"
#include "bas_util.hpp"
namespace opt_utilities
{
template <typename T>
void mnbrak(T& ax,T& bx,T& cx,T& fa,T& fb,T& fc,func_obj<T,T>& func)
{
const T GOLD=1.618034;
const T GLIMIT=100;
const T TINY=std::numeric_limits<T>::epsilon();
T ulim,u,r,q,fu;
fa=func.eval(ax);
fb=func.eval(bx);
if(fb>fa)
{
//shft(dum,ax,bx,dum);
//shft(dum,fb,fa,dum);
std::swap(ax,bx);
std::swap(fa,fb);
}
cx=bx+GOLD*(bx-ax);
fc=func.eval(cx);
while(fb>fc)
{
r=(bx-ax)*(fb-fc);
q=(bx-cx)*(fb-fa);
u=bx-T((bx-cx)*q-(bx-ax)*r)/
T(T(2.)*sign(T(max(T(tabs(T(q-r))),T(TINY))),T(q-r)));
ulim=bx+GLIMIT*(cx-bx);
if((bx-u)*(u-cx)>0.)
{
fu=func.eval(u);
if(fu<fc)
{
ax=bx;
bx=u;
fa=fb;
fb=fu;
return;
}
else if(fu>fb)
{
cx=u;
fc=fu;
return;
}
u=cx+GOLD*(cx-bx);
fu=func.eval(u);
}
else if((cx-u)*(u-ulim)>0.)
{
fu=func.eval(u);
if(fu<fc)
{
shft3(bx,cx,u,T(cx+GOLD*(cx-bx)));
shft3(fb,fc,fu,func.eval(u));
}
}
else if((u-ulim)*(ulim-cx)>=0)
{
u=ulim;
fu=func.eval(u);
}
else
{
u=cx+GOLD*(cx-bx);
fu=func.eval(u);
}
shft3(ax,bx,cx,u);
shft3(fa,fb,fc,fu);
}
}
}
#endif
|
