git.sesse.net Git - casparcg/blob - dependencies/boost/boost/math/tools/minima.hpp

   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 {
  25    BOOST_MATH_STD_USING
  26    bits = (std::min)(policies::digits<T, policies::policy<> >() / 2, bits);
  27    T tolerance = static_cast<T>(ldexp(1.0, 1-bits));
  28    T x;  // minima so far
  29    T w;  // second best point
  30    T v;  // previous value of w
  31    T u;  // most recent evaluation point
  32    T delta;  // The distance moved in the last step
  33    T delta2; // The distance moved in the step before last
  34    T fu, fv, fw, fx;  // function evaluations at u, v, w, x
  35    T mid; // midpoint of min and max
  36    T fract1, fract2;  // minimal relative movement in x
  37
  38    static const T golden = 0.3819660f;  // golden ratio, don't need too much precision here!
  39
  40    x = w = v = max;
  41    fw = fv = fx = f(x);
  42    delta2 = delta = 0;
  43
  44    uintmax_t count = max_iter;
  45
  46    do{
  47       // get midpoint
  48       mid = (min + max) / 2;
  49       // work out if we're done already:
  50       fract1 = tolerance * fabs(x) + tolerance / 4;
  51       fract2 = 2 * fract1;
  52       if(fabs(x - mid) <= (fract2 - (max - min) / 2))
  53          break;
  54
  55       if(fabs(delta2) > fract1)
  56       {
  57          // try and construct a parabolic fit:
  58          T r = (x - w) * (fx - fv);
  59          T q = (x - v) * (fx - fw);
  60          T p = (x - v) * q - (x - w) * r;
  61          q = 2 * (q - r);
  62          if(q > 0)
  63             p = -p;
  64          q = fabs(q);
  65          T td = delta2;
  66          delta2 = delta;
  67          // determine whether a parabolic step is acceptible or not:
  68          if((fabs(p) >= fabs(q * td / 2)) || (p <= q * (min - x)) || (p >= q * (max - x)))
  69          {
  70             // nope, try golden section instead
  71             delta2 = (x >= mid) ? min - x : max - x;
  72             delta = golden * delta2;
  73          }
  74          else
  75          {
  76             // whew, parabolic fit:
  77             delta = p / q;
  78             u = x + delta;
  79             if(((u - min) < fract2) || ((max- u) < fract2))
  80                delta = (mid - x) < 0 ? (T)-fabs(fract1) : (T)fabs(fract1);
  81          }
  82       }
  83       else
  84       {
  85          // golden section:
  86          delta2 = (x >= mid) ? min - x : max - x;
  87          delta = golden * delta2;
  88       }
  89       // update current position:
  90       u = (fabs(delta) >= fract1) ? T(x + delta) : (delta > 0 ? T(x + fabs(fract1)) : T(x - fabs(fract1)));
  91       fu = f(u);
  92       if(fu <= fx)
  93       {
  94          // good new point is an improvement!
  95          // update brackets:
  96          if(u >= x)
  97             min = x;
  98          else
  99             max = x;
 100          // update control points:
 101          v = w;
 102          w = x;
 103          x = u;
 104          fv = fw;
 105          fw = fx;
 106          fx = fu;
 107       }
 108       else
 109       {
 110          // Oh dear, point u is worse than what we have already,
 111          // even so it *must* be better than one of our endpoints:
 112          if(u < x)
 113             min = u;
 114          else
 115             max = u;
 116          if((fu <= fw) || (w == x))
 117          {
 118             // however it is at least second best:
 119             v = w;
 120             w = u;
 121             fv = fw;
 122             fw = fu;
 123          }
 124          else if((fu <= fv) || (v == x) || (v == w))
 125          {
 126             // third best:
 127             v = u;
 128             fv = fu;
 129          }
 130       }
 131
 132    }while(--count);
 133
 134    max_iter -= count;
 135
 136    return std::make_pair(x, fx);
 137 }
 138
 139 template <class F, class T>
 140 inline std::pair<T, T> brent_find_minima(F f, T min, T max, int digits)
 141 {
 142    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 143    return brent_find_minima(f, min, max, digits, m);
 144 }
 145
 146 }}} // namespaces
 147
 148 #endif
 149
 150
 151
 152