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