Chris@16
|
1 // (C) Copyright John Maddock 2006.
|
Chris@16
|
2 // Use, modification and distribution are subject to the
|
Chris@16
|
3 // Boost Software License, Version 1.0. (See accompanying file
|
Chris@16
|
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
Chris@16
|
5
|
Chris@16
|
6
|
Chris@16
|
7 #ifndef BOOST_MATH_TOOLS_MINIMA_HPP
|
Chris@16
|
8 #define BOOST_MATH_TOOLS_MINIMA_HPP
|
Chris@16
|
9
|
Chris@16
|
10 #ifdef _MSC_VER
|
Chris@16
|
11 #pragma once
|
Chris@16
|
12 #endif
|
Chris@16
|
13
|
Chris@16
|
14 #include <utility>
|
Chris@16
|
15 #include <boost/config/no_tr1/cmath.hpp>
|
Chris@16
|
16 #include <boost/math/tools/precision.hpp>
|
Chris@16
|
17 #include <boost/math/policies/policy.hpp>
|
Chris@16
|
18 #include <boost/cstdint.hpp>
|
Chris@16
|
19
|
Chris@16
|
20 namespace boost{ namespace math{ namespace tools{
|
Chris@16
|
21
|
Chris@16
|
22 template <class F, class T>
|
Chris@16
|
23 std::pair<T, T> brent_find_minima(F f, T min, T max, int bits, boost::uintmax_t& max_iter)
|
Chris@16
|
24 {
|
Chris@16
|
25 BOOST_MATH_STD_USING
|
Chris@16
|
26 bits = (std::min)(policies::digits<T, policies::policy<> >() / 2, bits);
|
Chris@16
|
27 T tolerance = static_cast<T>(ldexp(1.0, 1-bits));
|
Chris@16
|
28 T x; // minima so far
|
Chris@16
|
29 T w; // second best point
|
Chris@16
|
30 T v; // previous value of w
|
Chris@16
|
31 T u; // most recent evaluation point
|
Chris@16
|
32 T delta; // The distance moved in the last step
|
Chris@16
|
33 T delta2; // The distance moved in the step before last
|
Chris@16
|
34 T fu, fv, fw, fx; // function evaluations at u, v, w, x
|
Chris@16
|
35 T mid; // midpoint of min and max
|
Chris@16
|
36 T fract1, fract2; // minimal relative movement in x
|
Chris@16
|
37
|
Chris@16
|
38 static const T golden = 0.3819660f; // golden ratio, don't need too much precision here!
|
Chris@16
|
39
|
Chris@16
|
40 x = w = v = max;
|
Chris@16
|
41 fw = fv = fx = f(x);
|
Chris@16
|
42 delta2 = delta = 0;
|
Chris@16
|
43
|
Chris@16
|
44 uintmax_t count = max_iter;
|
Chris@16
|
45
|
Chris@16
|
46 do{
|
Chris@16
|
47 // get midpoint
|
Chris@16
|
48 mid = (min + max) / 2;
|
Chris@16
|
49 // work out if we're done already:
|
Chris@16
|
50 fract1 = tolerance * fabs(x) + tolerance / 4;
|
Chris@16
|
51 fract2 = 2 * fract1;
|
Chris@16
|
52 if(fabs(x - mid) <= (fract2 - (max - min) / 2))
|
Chris@16
|
53 break;
|
Chris@16
|
54
|
Chris@16
|
55 if(fabs(delta2) > fract1)
|
Chris@16
|
56 {
|
Chris@16
|
57 // try and construct a parabolic fit:
|
Chris@16
|
58 T r = (x - w) * (fx - fv);
|
Chris@16
|
59 T q = (x - v) * (fx - fw);
|
Chris@16
|
60 T p = (x - v) * q - (x - w) * r;
|
Chris@16
|
61 q = 2 * (q - r);
|
Chris@16
|
62 if(q > 0)
|
Chris@16
|
63 p = -p;
|
Chris@16
|
64 q = fabs(q);
|
Chris@16
|
65 T td = delta2;
|
Chris@16
|
66 delta2 = delta;
|
Chris@16
|
67 // determine whether a parabolic step is acceptible or not:
|
Chris@16
|
68 if((fabs(p) >= fabs(q * td / 2)) || (p <= q * (min - x)) || (p >= q * (max - x)))
|
Chris@16
|
69 {
|
Chris@16
|
70 // nope, try golden section instead
|
Chris@16
|
71 delta2 = (x >= mid) ? min - x : max - x;
|
Chris@16
|
72 delta = golden * delta2;
|
Chris@16
|
73 }
|
Chris@16
|
74 else
|
Chris@16
|
75 {
|
Chris@16
|
76 // whew, parabolic fit:
|
Chris@16
|
77 delta = p / q;
|
Chris@16
|
78 u = x + delta;
|
Chris@16
|
79 if(((u - min) < fract2) || ((max- u) < fract2))
|
Chris@16
|
80 delta = (mid - x) < 0 ? (T)-fabs(fract1) : (T)fabs(fract1);
|
Chris@16
|
81 }
|
Chris@16
|
82 }
|
Chris@16
|
83 else
|
Chris@16
|
84 {
|
Chris@16
|
85 // golden section:
|
Chris@16
|
86 delta2 = (x >= mid) ? min - x : max - x;
|
Chris@16
|
87 delta = golden * delta2;
|
Chris@16
|
88 }
|
Chris@16
|
89 // update current position:
|
Chris@16
|
90 u = (fabs(delta) >= fract1) ? T(x + delta) : (delta > 0 ? T(x + fabs(fract1)) : T(x - fabs(fract1)));
|
Chris@16
|
91 fu = f(u);
|
Chris@16
|
92 if(fu <= fx)
|
Chris@16
|
93 {
|
Chris@16
|
94 // good new point is an improvement!
|
Chris@16
|
95 // update brackets:
|
Chris@16
|
96 if(u >= x)
|
Chris@16
|
97 min = x;
|
Chris@16
|
98 else
|
Chris@16
|
99 max = x;
|
Chris@16
|
100 // update control points:
|
Chris@16
|
101 v = w;
|
Chris@16
|
102 w = x;
|
Chris@16
|
103 x = u;
|
Chris@16
|
104 fv = fw;
|
Chris@16
|
105 fw = fx;
|
Chris@16
|
106 fx = fu;
|
Chris@16
|
107 }
|
Chris@16
|
108 else
|
Chris@16
|
109 {
|
Chris@16
|
110 // Oh dear, point u is worse than what we have already,
|
Chris@16
|
111 // even so it *must* be better than one of our endpoints:
|
Chris@16
|
112 if(u < x)
|
Chris@16
|
113 min = u;
|
Chris@16
|
114 else
|
Chris@16
|
115 max = u;
|
Chris@16
|
116 if((fu <= fw) || (w == x))
|
Chris@16
|
117 {
|
Chris@16
|
118 // however it is at least second best:
|
Chris@16
|
119 v = w;
|
Chris@16
|
120 w = u;
|
Chris@16
|
121 fv = fw;
|
Chris@16
|
122 fw = fu;
|
Chris@16
|
123 }
|
Chris@16
|
124 else if((fu <= fv) || (v == x) || (v == w))
|
Chris@16
|
125 {
|
Chris@16
|
126 // third best:
|
Chris@16
|
127 v = u;
|
Chris@16
|
128 fv = fu;
|
Chris@16
|
129 }
|
Chris@16
|
130 }
|
Chris@16
|
131
|
Chris@16
|
132 }while(--count);
|
Chris@16
|
133
|
Chris@16
|
134 max_iter -= count;
|
Chris@16
|
135
|
Chris@16
|
136 return std::make_pair(x, fx);
|
Chris@16
|
137 }
|
Chris@16
|
138
|
Chris@16
|
139 template <class F, class T>
|
Chris@16
|
140 inline std::pair<T, T> brent_find_minima(F f, T min, T max, int digits)
|
Chris@16
|
141 {
|
Chris@16
|
142 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
|
Chris@16
|
143 return brent_find_minima(f, min, max, digits, m);
|
Chris@16
|
144 }
|
Chris@16
|
145
|
Chris@16
|
146 }}} // namespaces
|
Chris@16
|
147
|
Chris@16
|
148 #endif
|
Chris@16
|
149
|
Chris@16
|
150
|
Chris@16
|
151
|
Chris@16
|
152
|