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