minima.hpp 4.2 KB

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