| // (C) Copyright John Maddock 2006. |
| // Use, modification and distribution are subject to the |
| // Boost Software License, Version 1.0. (See accompanying file |
| // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
| |
| |
| #ifndef BOOST_MATH_TOOLS_MINIMA_HPP |
| #define BOOST_MATH_TOOLS_MINIMA_HPP |
| |
| #ifdef _MSC_VER |
| #pragma once |
| #endif |
| |
| #include <utility> |
| #include <boost/config/no_tr1/cmath.hpp> |
| #include <boost/math/tools/precision.hpp> |
| #include <boost/math/policies/policy.hpp> |
| #include <boost/cstdint.hpp> |
| |
| namespace boost{ namespace math{ namespace tools{ |
| |
| template <class F, class T> |
| std::pair<T, T> brent_find_minima(F f, T min, T max, int bits, boost::uintmax_t& max_iter) |
| { |
| BOOST_MATH_STD_USING |
| bits = (std::min)(policies::digits<T, policies::policy<> >() / 2, bits); |
| T tolerance = static_cast<T>(ldexp(1.0, 1-bits)); |
| T x; // minima so far |
| T w; // second best point |
| T v; // previous value of w |
| T u; // most recent evaluation point |
| T delta; // The distance moved in the last step |
| T delta2; // The distance moved in the step before last |
| T fu, fv, fw, fx; // function evaluations at u, v, w, x |
| T mid; // midpoint of min and max |
| T fract1, fract2; // minimal relative movement in x |
| |
| static const T golden = 0.3819660f; // golden ratio, don't need too much precision here! |
| |
| x = w = v = max; |
| fw = fv = fx = f(x); |
| delta2 = delta = 0; |
| |
| uintmax_t count = max_iter; |
| |
| do{ |
| // get midpoint |
| mid = (min + max) / 2; |
| // work out if we're done already: |
| fract1 = tolerance * fabs(x) + tolerance / 4; |
| fract2 = 2 * fract1; |
| if(fabs(x - mid) <= (fract2 - (max - min) / 2)) |
| break; |
| |
| if(fabs(delta2) > fract1) |
| { |
| // try and construct a parabolic fit: |
| T r = (x - w) * (fx - fv); |
| T q = (x - v) * (fx - fw); |
| T p = (x - v) * q - (x - w) * r; |
| q = 2 * (q - r); |
| if(q > 0) |
| p = -p; |
| q = fabs(q); |
| T td = delta2; |
| delta2 = delta; |
| // determine whether a parabolic step is acceptible or not: |
| if((fabs(p) >= fabs(q * td / 2)) || (p <= q * (min - x)) || (p >= q * (max - x))) |
| { |
| // nope, try golden section instead |
| delta2 = (x >= mid) ? min - x : max - x; |
| delta = golden * delta2; |
| } |
| else |
| { |
| // whew, parabolic fit: |
| delta = p / q; |
| u = x + delta; |
| if(((u - min) < fract2) || ((max- u) < fract2)) |
| delta = (mid - x) < 0 ? (T)-fabs(fract1) : (T)fabs(fract1); |
| } |
| } |
| else |
| { |
| // golden section: |
| delta2 = (x >= mid) ? min - x : max - x; |
| delta = golden * delta2; |
| } |
| // update current position: |
| u = (fabs(delta) >= fract1) ? T(x + delta) : (delta > 0 ? T(x + fabs(fract1)) : T(x - fabs(fract1))); |
| fu = f(u); |
| if(fu <= fx) |
| { |
| // good new point is an improvement! |
| // update brackets: |
| if(u >= x) |
| min = x; |
| else |
| max = x; |
| // update control points: |
| v = w; |
| w = x; |
| x = u; |
| fv = fw; |
| fw = fx; |
| fx = fu; |
| } |
| else |
| { |
| // Oh dear, point u is worse than what we have already, |
| // even so it *must* be better than one of our endpoints: |
| if(u < x) |
| min = u; |
| else |
| max = u; |
| if((fu <= fw) || (w == x)) |
| { |
| // however it is at least second best: |
| v = w; |
| w = u; |
| fv = fw; |
| fw = fu; |
| } |
| else if((fu <= fv) || (v == x) || (v == w)) |
| { |
| // third best: |
| v = u; |
| fv = fu; |
| } |
| } |
| |
| }while(--count); |
| |
| max_iter -= count; |
| |
| return std::make_pair(x, fx); |
| } |
| |
| template <class F, class T> |
| inline std::pair<T, T> brent_find_minima(F f, T min, T max, int digits) |
| { |
| boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)(); |
| return brent_find_minima(f, min, max, digits, m); |
| } |
| |
| }}} // namespaces |
| |
| #endif |
| |
| |
| |
| |