math/tools/minima.hpp

0001 //  (C) Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006
0007 #ifndef BOOST_MATH_TOOLS_MINIMA_HPP
0008 #define BOOST_MATH_TOOLS_MINIMA_HPP
0009
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #endif
0013
0014 #include <cstdint>
0015 #include <cmath>
0016 #include <utility>
0017 #include <boost/math/tools/precision.hpp>
0018 #include <boost/math/policies/policy.hpp>
0019
0020 namespace boost{ namespace math{ namespace tools{
0021
0022 template <class F, class T>
0023 std::pair<T, T> brent_find_minima(F f, T min, T max, int bits, std::uintmax_t& max_iter)
0024    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
0025 {
0026    BOOST_MATH_STD_USING
0027    bits = (std::min)(policies::digits<T, policies::policy<> >() / 2, bits);
0028    T tolerance = static_cast<T>(ldexp(1.0, 1-bits));
0029    T x;  // minima so far
0030    T w;  // second best point
0031    T v;  // previous value of w
0032    T u;  // most recent evaluation point
0033    T delta;  // The distance moved in the last step
0034    T delta2; // The distance moved in the step before last
0035    T fu, fv, fw, fx;  // function evaluations at u, v, w, x
0036    T mid; // midpoint of min and max
0037    T fract1, fract2;  // minimal relative movement in x
0038
0039    static const T golden = 0.3819660f;  // golden ratio, don't need too much precision here!
0040
0041    x = w = v = max;
0042    fw = fv = fx = f(x);
0043    delta2 = delta = 0;
0044
0045    uintmax_t count = max_iter;
0046
0047    do{
0048       // get midpoint
0049       mid = (min + max) / 2;
0050       // work out if we're done already:
0051       fract1 = tolerance * fabs(x) + tolerance / 4;
0052       fract2 = 2 * fract1;
0053       if(fabs(x - mid) <= (fract2 - (max - min) / 2))
0054          break;
0055
0056       if(fabs(delta2) > fract1)
0057       {
0058          // try and construct a parabolic fit:
0059          T r = (x - w) * (fx - fv);
0060          T q = (x - v) * (fx - fw);
0061          T p = (x - v) * q - (x - w) * r;
0062          q = 2 * (q - r);
0063          if(q > 0)
0064             p = -p;
0065          q = fabs(q);
0066          T td = delta2;
0067          delta2 = delta;
0068          // determine whether a parabolic step is acceptable or not:
0069          if((fabs(p) >= fabs(q * td / 2)) || (p <= q * (min - x)) || (p >= q * (max - x)))
0070          {
0071             // nope, try golden section instead
0072             delta2 = (x >= mid) ? min - x : max - x;
0073             delta = golden * delta2;
0074          }
0075          else
0076          {
0077             // whew, parabolic fit:
0078             delta = p / q;
0079             u = x + delta;
0080             if(((u - min) < fract2) || ((max- u) < fract2))
0081                delta = (mid - x) < 0 ? (T)-fabs(fract1) : (T)fabs(fract1);
0082          }
0083       }
0084       else
0085       {
0086          // golden section:
0087          delta2 = (x >= mid) ? min - x : max - x;
0088          delta = golden * delta2;
0089       }
0090       // update current position:
0091       u = (fabs(delta) >= fract1) ? T(x + delta) : (delta > 0 ? T(x + fabs(fract1)) : T(x - fabs(fract1)));
0092       fu = f(u);
0093       if(fu <= fx)
0094       {
0095          // good new point is an improvement!
0096          // update brackets:
0097          if(u >= x)
0098             min = x;
0099          else
0100             max = x;
0101          // update control points:
0102          v = w;
0103          w = x;
0104          x = u;
0105          fv = fw;
0106          fw = fx;
0107          fx = fu;
0108       }
0109       else
0110       {
0111          // Oh dear, point u is worse than what we have already,
0112          // even so it *must* be better than one of our endpoints:
0113          if(u < x)
0114             min = u;
0115          else
0116             max = u;
0117          if((fu <= fw) || (w == x))
0118          {
0119             // however it is at least second best:
0120             v = w;
0121             w = u;
0122             fv = fw;
0123             fw = fu;
0124          }
0125          else if((fu <= fv) || (v == x) || (v == w))
0126          {
0127             // third best:
0128             v = u;
0129             fv = fu;
0130          }
0131       }
0132
0133    }while(--count);
0134
0135    max_iter -= count;
0136
0137    return std::make_pair(x, fx);
0138 }
0139
0140 template <class F, class T>
0141 inline std::pair<T, T> brent_find_minima(F f, T min, T max, int digits)
0142    noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
0143 {
0144    std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)();
0145    return brent_find_minima(f, min, max, digits, m);
0146 }
0147
0148 }}} // namespaces
0149
0150 #endif
0151
0152
0153
0154