]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
1e59de90 TL |
14 | #include <cstdint> |
15 | #include <cmath> | |
7c673cae | 16 | #include <utility> |
7c673cae FG |
17 | #include <boost/math/tools/precision.hpp> |
18 | #include <boost/math/policies/policy.hpp> | |
7c673cae FG |
19 | |
20 | namespace boost{ namespace math{ namespace tools{ | |
21 | ||
22 | template <class F, class T> | |
1e59de90 TL |
23 | std::pair<T, T> brent_find_minima(F f, T min, T max, int bits, std::uintmax_t& max_iter) |
24 | noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>()))) | |
7c673cae FG |
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; | |
f67539c2 | 68 | // determine whether a parabolic step is acceptable or not: |
7c673cae FG |
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) | |
1e59de90 | 142 | noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>()))) |
7c673cae | 143 | { |
1e59de90 | 144 | std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)(); |
7c673cae FG |
145 | return brent_find_minima(f, min, max, digits, m); |
146 | } | |
147 | ||
148 | }}} // namespaces | |
149 | ||
150 | #endif | |
151 | ||
152 | ||
153 | ||
154 |