16 #ifndef dealii_line_minimization_h 17 #define dealii_line_minimization_h 19 #include <deal.II/base/config.h> 21 #include <deal.II/base/exceptions.h> 22 #include <deal.II/base/logstream.h> 23 #include <deal.II/base/numbers.h> 24 #include <deal.II/base/utilities.h> 26 #include <deal.II/numerics/history.h> 28 #include <boost/optional.hpp> 37 DEAL_II_NAMESPACE_OPEN
52 template <
typename NumberType>
53 boost::optional<NumberType>
55 const NumberType f_low,
56 const NumberType g_low,
57 const NumberType x_hi,
58 const NumberType f_hi);
69 template <
typename NumberType>
70 boost::optional<NumberType>
72 const NumberType f_low,
73 const NumberType g_low,
74 const NumberType x_hi,
75 const NumberType f_hi,
76 const NumberType g_hi);
85 template <
typename NumberType>
86 boost::optional<NumberType>
88 const NumberType f_low,
89 const NumberType g_low,
90 const NumberType x_hi,
91 const NumberType f_hi,
92 const NumberType x_rec,
93 const NumberType f_rec);
106 template <
typename NumberType>
109 const NumberType f_low,
110 const NumberType g_low,
111 const NumberType x_hi,
112 const NumberType f_hi,
113 const NumberType g_hi,
117 const std::pair<NumberType, NumberType> bounds);
123 template <
typename NumberType>
126 const NumberType f_low,
127 const NumberType g_low,
128 const NumberType x_hi,
129 const NumberType f_hi,
130 const NumberType g_hi,
134 const std::pair<NumberType, NumberType> bounds);
323 template <
typename NumberType>
324 std::pair<NumberType, unsigned int>
326 const std::function<std::pair<NumberType, NumberType>(
const NumberType x)>
331 NumberType(
const NumberType x_low,
332 const NumberType f_low,
333 const NumberType g_low,
334 const NumberType x_hi,
335 const NumberType f_hi,
336 const NumberType g_hi,
340 const std::pair<NumberType, NumberType> bounds)> &interpolate,
342 const NumberType eta = 0.9,
343 const NumberType
mu = 0.01,
344 const NumberType a_max = std::numeric_limits<NumberType>::max(),
345 const unsigned int max_evaluations = 20,
346 const bool debug_output =
false);
355 template <
typename NumberType>
356 boost::optional<NumberType>
364 const NumberType denom = (2. * g1 * x2 - 2. * g1 * x1 - 2. * f2 + 2. * f1);
368 return (g1 * (x2 * x2 - x1 * x1) + 2. * (f1 - f2) * x1) / denom;
373 template <
typename NumberType>
374 boost::optional<NumberType>
383 const NumberType beta1 = g1 + g2 - 3. * (f1 - f2) / (x1 - x2);
384 const NumberType s = beta1 * beta1 - g1 * g2;
388 const NumberType beta2 = std::sqrt(s);
389 const NumberType denom =
390 x1 < x2 ? g2 - g1 + 2. * beta2 : g1 - g2 + 2. * beta2;
394 return x1 < x2 ? x2 - (x2 - x1) * (g2 + beta2 - beta1) / denom :
395 x1 - (x1 - x2) * (g1 + beta2 - beta1) / denom;
400 template <
typename NumberType>
401 boost::optional<NumberType>
419 const NumberType x2_shift = x2 - x1;
420 const NumberType x3_shift = x3 - x1;
421 const NumberType r1 = f2 - f1 - g1 * x2_shift;
422 const NumberType r2 = f3 - f1 - g1 * x3_shift;
423 const NumberType denom =
424 std::pow(x2_shift * x3_shift, 2) * (x2_shift - x3_shift);
429 (r1 * std::pow(x3_shift, 2) - r2 * std::pow(x2_shift, 2)) / denom;
431 (r2 * std::pow(x2_shift, 3) - r1 * std::pow(x3_shift, 3)) / denom;
432 const NumberType &C = g1;
435 const NumberType radical = B * B - A * C * 3;
439 return x1 + (-B + std::sqrt(radical)) / (A * 3);
444 template <
typename NumberType>
455 const std::pair<NumberType, NumberType> bounds)
464 boost::optional<NumberType> res =
cubic_fit(x1, f1, g1, x2, f2, g2);
465 if (res && *res >= bounds.first && *res <= bounds.second)
470 if (res && *res >= bounds.first && *res <= bounds.second)
475 return (bounds.first + bounds.second) * 0.5;
480 template <
typename NumberType>
491 const std::pair<NumberType, NumberType> bounds)
500 boost::optional<NumberType> res =
504 if (res && *res >= bounds.first && *res <= bounds.second)
509 if (res && *res >= bounds.first && *res <= bounds.second)
514 return (bounds.first + bounds.second) * 0.5;
519 template <
typename NumberType>
520 std::pair<NumberType, unsigned int>
522 const std::function<std::pair<NumberType, NumberType>(
const NumberType x)>
527 NumberType(
const NumberType x_low,
528 const NumberType f_low,
529 const NumberType g_low,
530 const NumberType x_hi,
531 const NumberType f_hi,
532 const NumberType g_hi,
536 const std::pair<NumberType, NumberType> bounds)> &choose,
538 const NumberType eta,
540 const NumberType a_max,
541 const unsigned int max_evaluations,
542 const bool debug_output)
553 const NumberType tau1 = 9.;
558 const NumberType tau2 = 0.1;
559 const NumberType tau3 = 0.5;
561 const NumberType g0_abs = std::abs(g0);
562 const NumberType f_min = f0 + a_max *
mu * g0;
566 const auto w1 = [&](
const NumberType a,
const NumberType f) {
567 return f <= f0 + a *
mu * g0;
572 const auto w2 = [&](
const NumberType g) {
573 return std::abs(g) <= eta * g0_abs;
579 const NumberType x = std::numeric_limits<NumberType>::signaling_NaN();
580 NumberType a_lo = x, f_lo = x, g_lo = x;
581 NumberType a_hi = x, f_hi = x, g_hi = x;
582 NumberType ai = x, fi = x, gi = x;
587 NumberType f_prev, g_prev, a_prev;
593 while (i < max_evaluations)
595 const auto fgi = func(ai);
601 deallog <<
"Bracketing phase: " << i << std::endl
602 << ai <<
" " << fi <<
" " << gi <<
" " << w1(ai, fi) <<
" " 603 << w2(gi) <<
" " << f_min << std::endl;
606 if (fi <= f_min || ai == a_max)
609 deallog <<
"Reached the maximum step size." << std::endl;
610 return std::make_pair(ai, i);
614 (fi >= f_prev && i > 1))
629 deallog <<
"Satisfied both Wolfe conditions during Bracketing." 633 return std::make_pair(ai, i);
650 std::make_pair(2. * ai - a_prev,
651 std::min(a_max, ai + tau1 * (ai - a_prev)));
668 "Could not find the initial bracket within the given number of iterations."));
685 if (std::abs(a_lo) > std::numeric_limits<NumberType>::epsilon() &&
686 std::abs(a_hi) > std::numeric_limits<NumberType>::epsilon())
694 while (i < max_evaluations)
696 const NumberType a_lo_safe = a_lo + tau2 * (a_hi - a_lo);
697 const NumberType a_hi_safe = a_hi - tau3 * (a_hi - a_lo);
698 const auto bounds = std::minmax(a_lo_safe, a_hi_safe);
701 a_lo, f_lo, g_lo, a_hi, f_hi, g_hi, a_rec, f_rec, g_rec, bounds);
703 const std::pair<NumberType, NumberType> fgi = func(ai);
709 deallog <<
"Sectioning phase: " << i << std::endl
710 << a_lo <<
" " << f_lo <<
" " << g_lo <<
" " << w1(a_lo, f_lo)
711 <<
" " << w2(g_lo) << std::endl
712 << a_hi <<
" " << f_hi <<
" " << g_hi <<
" " << w1(a_hi, f_hi)
713 <<
" " << w2(g_hi) << std::endl
714 << ai <<
" " << fi <<
" " << gi <<
" " << w1(ai, fi) <<
" " 715 << w2(gi) << std::endl;
717 if (!w1(ai, fi) || fi >= f_lo)
733 deallog <<
"Satisfied both Wolfe conditions." << std::endl;
735 return std::make_pair(ai, i);
738 if (gi * (a_hi - a_lo) >= 0)
767 "Could not could complete the sectioning phase within the given number of iterations."));
768 return std::make_pair(std::numeric_limits<NumberType>::signaling_NaN(), i);
775 DEAL_II_NAMESPACE_CLOSE
777 #endif // dealii_line_minimization_h #define AssertDimension(dim1, dim2)
boost::optional< NumberType > cubic_fit_three_points(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi, const NumberType x_rec, const NumberType f_rec)
#define AssertThrow(cond, exc)
boost::optional< NumberType > cubic_fit(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi, const NumberType g_hi)
static ::ExceptionBase & ExcMessage(std::string arg1)
NumberType poly_fit_three_points(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi, const NumberType g_hi, const FiniteSizeHistory< NumberType > &x_rec, const FiniteSizeHistory< NumberType > &f_rec, const FiniteSizeHistory< NumberType > &g_rec, const std::pair< NumberType, NumberType > bounds)
#define Assert(cond, exc)
NumberType poly_fit(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi, const NumberType g_hi, const FiniteSizeHistory< NumberType > &x_rec, const FiniteSizeHistory< NumberType > &f_rec, const FiniteSizeHistory< NumberType > &g_rec, const std::pair< NumberType, NumberType > bounds)
void add(const T &element)
boost::optional< NumberType > quadratic_fit(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi)
std::pair< NumberType, unsigned int > line_search(const std::function< std::pair< NumberType, NumberType >(const NumberType x)> &func, const NumberType f0, const NumberType g0, const std::function< NumberType(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi, const NumberType g_hi, const FiniteSizeHistory< NumberType > &x_rec, const FiniteSizeHistory< NumberType > &f_rec, const FiniteSizeHistory< NumberType > &g_rec, const std::pair< NumberType, NumberType > bounds)> &interpolate, const NumberType a1, const NumberType eta=0.9, const NumberType mu=0.01, const NumberType a_max=std::numeric_limits< NumberType >::max(), const unsigned int max_evaluations=20, const bool debug_output=false)
static ::ExceptionBase & ExcInternalError()