--- /dev/null
+New: Add line minimization functions.
+<br>
+(Denis Davydov 2018/09/13)
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+
+#ifndef dealii_line_minimization_h
+#define dealii_line_minimization_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/numbers.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/numerics/history.h>
+
+#include <boost/optional.hpp>
+
+#include <errno.h>
+#include <sys/stat.h>
+
+#include <fstream>
+#include <string>
+
+DEAL_II_NAMESPACE_OPEN
+
+using namespace dealii;
+
+/**
+ * A namespace for various algorithms related to minimization over line.
+ */
+namespace LineMinimization
+{
+ /**
+ * Given $x\_low$ and $x\_hi$ together with values of function
+ * $f(x\_low)$ and $f(x\_hi)$ and the gradient $g(x\_low)$, return the local
+ * minimizer of the quadratic interpolation function.
+ *
+ * The return type is optional to fit with similar function which may
+ * not have solution for given parameters.
+ */
+ template <typename NumberType>
+ 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);
+
+ /**
+ * Given $x\_low$ and $x\_hi$ together with values of function
+ * $f(x\_low)$ and $f(x\_hi)$) and its gradients ($g(x\_low)*g(x\_hi) < 0$) at
+ * those points, return the local minimizer of the cubic interpolation
+ * function. That is the location where the cubic interpolation function
+ * attains its minimum.
+ *
+ * The return type is optional as the real-valued solution might not exist.
+ */
+ template <typename NumberType>
+ 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);
+
+ /**
+ * Find the minimizer of a cubic polynomial that goes through the
+ * points $f\_low=f(x\_low)$, $f\_hi=f(x\_hi)$ and $f\_rec(x\_rec)$
+ * and has derivatve $g\_low$ at $x\_low$.
+ *
+ * The return type is optional as the real-valued solution might not exist.
+ */
+ template <typename NumberType>
+ 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);
+
+ /**
+ * Return the minimizer of a polynomial using function values @p f_low @p f_hi @p f_rec[0] at three points
+ * @p x_low @p x_hi @p x_rec[0] and derivatives at two points @p g_low and @p g_hi. The returned point
+ * should be within the bounds @p bounds .
+ *
+ * This function will first try the cubic fit (see cubic_fit ). If it's
+ * unsuccessfull or not
+ * within the provided @p bounds the quadratic fit will be performed (see quadratic_fit ). The function will
+ * fallback to bisection if quadratic fit fails.
+ */
+ template <typename NumberType>
+ 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);
+
+ /**
+ * Same as above but doing cubic fit with three points (see
+ * cubic_fit_three_points ).
+ */
+ template <typename NumberType>
+ 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);
+
+
+ /**
+ * Perform a line search in $(0,max]$ with strong Wolfe conditions
+ * \f[
+ * f(\alpha) \le f(0) + \alpha \mu f'(0) \\
+ * |f'(\alpha)| \le \eta |f'(0)|
+ * \f]
+ * using one dimensional
+ * functions @p func and a function @p interpolate to choose a new point
+ * from the interval based on the function values and derivatives at its ends.
+ * @p a1 is a trial estimate of the first step.
+ * Interpolation can be done using poly_fit or poly_fit_three_points .
+ *
+ * The function implements Algorithms
+ * 2.6.2 and 2.6.4 on pages 34-35 in Fletcher, 2013, Practical methods of
+ * optimization. these are minor variations of the Algorithm 3.5 and 3.6 on
+ * pages 60-61 in Nocedal and Wright, Numerical optimization.
+ * It consists of a bracketing phase and a zoom phase, where @p interpolate is used.
+ *
+ * The function returns the step size and the number of times function @p func was called.
+ *
+ * @param func a one dimensional function which returns value and derivative
+ * at the given point.
+ * @param f0 function value the origin
+ * @param g0 function derivative the origin
+ * @param interpolate a function which determines how interpolation is done
+ * during the zoom phase. It takes values and derivatives at the current
+ * interval/bracket ($f\_low$, $f\_hi$) as well as up to 5 values and
+ * derivatives at previous steps. The returned value is to be provided within
+ * the given bounds.
+ * @param a1 initial trial step for bracketing phase
+ * @param eta a parameter in the second Wolfe condition (curvature condition)
+ * @param mu a parameter in the first Wolfe condition (sufficient decrease)
+ * @param a_max maximum allowed step size
+ * @param max_evaluations maximum allowed number of function evaluations
+ * @param debug_output a flag to do extra debug output into deallog static
+ * object
+ */
+ template <typename NumberType>
+ 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);
+
+ // ------------------- inline and template functions ----------------
+#ifndef DOXYGEN
+
+ template <typename NumberType>
+ boost::optional<NumberType>
+ quadratic_fit(const NumberType x1,
+ const NumberType f1,
+ const NumberType g1,
+ const NumberType x2,
+ const NumberType f2)
+ {
+ Assert(x1 != x2, ExcMessage("Point are the same"));
+ const NumberType denom = (2. * g1 * x2 - 2. * g1 * x1 - 2. * f2 + 2. * f1);
+ if (denom == 0)
+ return boost::none;
+ else
+ return (g1 * (x2 * x2 - x1 * x1) + 2. * (f1 - f2) * x1) / denom;
+ }
+
+ template <typename NumberType>
+ boost::optional<NumberType>
+ cubic_fit(const NumberType x1,
+ const NumberType f1,
+ const NumberType g1,
+ const NumberType x2,
+ const NumberType f2,
+ const NumberType g2)
+ {
+ Assert(x1 != x2, ExcMessage("Points are the same"));
+ const NumberType beta1 = g1 + g2 - 3. * (f1 - f2) / (x1 - x2);
+ const NumberType s = beta1 * beta1 - g1 * g2;
+ if (s < 0)
+ return boost::none;
+
+ const NumberType beta2 = std::sqrt(s);
+ const NumberType denom =
+ x1 < x2 ? g2 - g1 + 2. * beta2 : g1 - g2 + 2. * beta2;
+ if (denom == 0.)
+ return boost::none;
+
+ return x1 < x2 ? x2 - (x2 - x1) * (g2 + beta2 - beta1) / denom :
+ x1 - (x1 - x2) * (g1 + beta2 - beta1) / denom;
+ }
+
+
+
+ template <typename NumberType>
+ boost::optional<NumberType>
+ cubic_fit_three_points(const NumberType x1,
+ const NumberType f1,
+ const NumberType g1,
+ const NumberType x2,
+ const NumberType f2,
+ const NumberType x3,
+ const NumberType f3)
+ {
+ Assert(x1 != x2, ExcMessage("Points are the same"));
+ Assert(x1 != x3, ExcMessage("Points are the same"));
+ // f(x) = A *(x-x1)^3 + B*(x-x1)^2 + C*(x-x1) + D
+ // =>
+ // D = f1
+ // C = g1
+
+ // the rest is a system of 2 equations:
+
+ const NumberType x2_shift = x2 - x1;
+ const NumberType x3_shift = x3 - x1;
+ const NumberType r1 = f2 - f1 - g1 * x2_shift;
+ const NumberType r2 = f3 - f1 - g1 * x3_shift;
+ const NumberType denom =
+ std::pow(x2_shift * x3_shift, 2) * (x2_shift - x3_shift);
+ if (denom == 0.)
+ return boost::none;
+
+ const NumberType A =
+ (r1 * std::pow(x3_shift, 2) - r2 * std::pow(x2_shift, 2)) / denom;
+ const NumberType B =
+ (r2 * std::pow(x2_shift, 3) - r1 * std::pow(x3_shift, 3)) / denom;
+ const NumberType &C = g1;
+
+ // now get the minimizer:
+ const NumberType radical = B * B - A * C * 3;
+ if (radical < 0)
+ return boost::none;
+
+ return x1 + (-B + std::sqrt(radical)) / (A * 3);
+ }
+
+
+ template <typename NumberType>
+ NumberType
+ poly_fit(const NumberType x1,
+ const NumberType f1,
+ const NumberType g1,
+ const NumberType x2,
+ const NumberType f2,
+ const NumberType g2,
+ const FiniteSizeHistory<NumberType> &,
+ const FiniteSizeHistory<NumberType> &,
+ const FiniteSizeHistory<NumberType> &,
+ const std::pair<NumberType, NumberType> bounds)
+ {
+ Assert(bounds.first < bounds.second, ExcMessage("Incorrect bounds"));
+
+ // Similar to scipy implementation but we fit based on two points
+ // with their gradients and do bisection on bounds.
+ // https://github.com/scipy/scipy/blob/v1.0.0/scipy/optimize/linesearch.py#L555-L563
+
+ // First try cubic interpolation
+ boost::optional<NumberType> res = cubic_fit(x1, f1, g1, x2, f2, g2);
+ if (res && *res >= bounds.first && *res <= bounds.second)
+ return *res;
+
+ // cubic either fails or outside of safe region, do quadratic:
+ res = quadratic_fit(x1, f1, g1, x2, f2);
+ if (res && *res >= bounds.first && *res <= bounds.second)
+ return *res;
+
+ // quadratic either failed or outside of safe region. Do bisection
+ // on safe region
+ return (bounds.first + bounds.second) * 0.5;
+ }
+
+
+
+ template <typename NumberType>
+ NumberType
+ poly_fit_three_points(const NumberType x1,
+ const NumberType f1,
+ const NumberType g1,
+ const NumberType x2,
+ const NumberType f2,
+ const NumberType g2,
+ const FiniteSizeHistory<NumberType> &x_rec,
+ const FiniteSizeHistory<NumberType> &f_rec,
+ const FiniteSizeHistory<NumberType> & /*g_rec*/,
+ const std::pair<NumberType, NumberType> bounds)
+ {
+ Assert(bounds.first < bounds.second, ExcMessage("Incorrect bounds"));
+ AssertDimension(x_rec.size(), f_rec.size());
+
+ // Same as scipy implementation where cubic fit is using 3 points
+ // https://github.com/scipy/scipy/blob/v1.0.0/scipy/optimize/linesearch.py#L555-L563
+
+ // First try cubic interpolation after first iteration
+ boost::optional<NumberType> res =
+ x_rec.size() > 0 ?
+ cubic_fit_three_points(x1, f1, g1, x2, f2, x_rec[0], f_rec[0]) :
+ boost::none;
+ if (res && *res >= bounds.first && *res <= bounds.second)
+ return *res;
+
+ // cubic either fails or outside of safe region, do quadratic:
+ res = quadratic_fit(x1, f1, g1, x2, f2);
+ if (res && *res >= bounds.first && *res <= bounds.second)
+ return *res;
+
+ // quadratic either failed or outside of safe region. Do bisection
+ // on safe region
+ return (bounds.first + bounds.second) * 0.5;
+ }
+
+
+ template <typename NumberType>
+ 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)> &choose,
+ const NumberType a1,
+ const NumberType eta,
+ const NumberType mu,
+ const NumberType a_max,
+ const unsigned int max_evaluations,
+ const bool debug_output)
+ {
+ // Note that scipy use dcsrch() from Minpack2 Fortran lib for line search
+ Assert(mu < 0.5 && mu > 0, ExcMessage("mu is not in (0,1/2)."));
+ Assert(eta < 1. && eta > mu, ExcMessage("eta is not in (mu,1)."));
+ Assert(a_max > 0, ExcMessage("max is not positive."));
+ Assert(a1 > 0 && a1 <= a_max, ExcMessage("a1 is not in (0,max]."));
+ Assert(g0 < 0, ExcMessage("Initial slope is not negative"));
+
+ // Growth parameter for bracketing phase:
+ // 1 < tau1
+ const NumberType tau1 = 9.;
+ // shrink parameters for sectioning phase to prevent ai from being
+ // arbitrary close to the extremes of the interval.
+ // 0 < tau2 < tau3 <= 1/2
+ // tau2 <= eta is advisable
+ const NumberType tau2 = 0.1; // bound for closeness to a_lo
+ const NumberType tau3 = 0.5; // bound for closeness to a_hi
+
+ const NumberType g0_abs = std::abs(g0);
+ const NumberType f_min = f0 + a_max * mu * g0;
+
+ // return True if the first Wolfe condition (sufficient decrease) is
+ // satisfied
+ const auto w1 = [&](const NumberType a, const NumberType f) {
+ return f <= f0 + a * mu * g0;
+ };
+
+ // return True if the second Wolfe condition (curvature condition) is
+ // satisfied
+ const auto w2 = [&](const NumberType g) {
+ return std::abs(g) <= eta * g0_abs;
+ };
+
+ // Bracketing phase (Algorithm 2.6.2): look for a non-trivial interval
+ // which is known to contain an interval of acceptable points.
+ // We adopt notation of Noceal.
+ const NumberType x = std::numeric_limits<NumberType>::signaling_NaN();
+ NumberType a_lo = x, f_lo = x, g_lo = x;
+ NumberType a_hi = x, f_hi = x, g_hi = x;
+ NumberType ai = x, fi = x, gi = x;
+
+ // count function calls in i:
+ unsigned int i = 0;
+ {
+ NumberType f_prev, g_prev, a_prev;
+ ai = a1;
+ f_prev = f0;
+ g_prev = g0;
+ a_prev = 0;
+
+ while (i < max_evaluations)
+ {
+ const auto fgi = func(ai);
+ fi = fgi.first;
+ gi = fgi.second;
+ i++;
+
+ if (debug_output)
+ deallog << "Bracketing phase: " << i << std::endl
+ << ai << " " << fi << " " << gi << std::endl;
+
+ // first check if we can stop bracketing or the whole line search:
+ if (fi <= f_min || ai == a_max)
+ return std::make_pair(ai, i);
+
+ if (!w1(ai, fi) ||
+ (fi >= f_prev && i > 1)) // violate first Wolfe or not descending
+ {
+ a_lo = a_prev;
+ f_lo = f_prev;
+ g_lo = g_prev;
+
+ a_hi = ai;
+ f_hi = fi;
+ g_hi = gi;
+ break; // end bracketing
+ }
+
+ if (w2(gi)) // satisfies both Wolfe conditions
+ {
+ Assert(w1(ai, fi), ExcInternalError());
+ return std::make_pair(ai, i);
+ }
+
+ if (gi >= 0) // not descending
+ {
+ a_lo = ai;
+ f_lo = fi;
+ g_lo = gi;
+
+ a_hi = a_prev;
+ f_hi = f_prev;
+ g_hi = g_prev;
+ break; // end bracketing
+ }
+
+ // extrapolation step with the bounds
+ const auto bounds =
+ std::make_pair(2. * ai - a_prev,
+ std::min(a_max, ai + tau1 * (ai - a_prev)));
+
+ a_prev = ai;
+ f_prev = fi;
+ g_prev = gi;
+
+ // NOTE: Fletcher's 2.6.2 includes optional extrapolation, we
+ // simply take the upper bound
+ // Scipy increases by factor of two:
+ // https://github.com/scipy/scipy/blob/v1.0.0/scipy/optimize/linesearch.py#L447
+ ai = bounds.second;
+ }
+ }
+
+ AssertThrow(
+ i < max_evaluations,
+ ExcMessage(
+ "Could not find the initial bracket within the given number of iterations."));
+
+ // Check properties of the bracket (Theorem 3.2 in More and Thuente, 94
+ // and Eq. 2.6.3 in Fletcher 2013
+
+ // FIXME: these conditions are actually violated for Fig3 and a1=10^3 in
+ // More and Thorenton, 94.
+
+ /*
+ Assert((f_lo < f_hi) && w1(a_lo, f_lo), ExcInternalError());
+ Assert(((a_hi - a_lo) * g_lo < 0) && !w2(g_lo), ExcInternalError());
+ Assert((w1(a_hi, f_hi) || f_hi >= f_lo), ExcInternalError());
+ */
+
+ // keep short history of last points to improve interpolation
+ FiniteSizeHistory<NumberType> a_rec(5), f_rec(5), g_rec(5);
+ // if neither a_lo nor a_hi are zero:
+ if (std::abs(a_lo) > std::numeric_limits<NumberType>::epsilon() &&
+ std::abs(a_hi) > std::numeric_limits<NumberType>::epsilon())
+ {
+ a_rec.add(0);
+ f_rec.add(f0);
+ g_rec.add(g0);
+ }
+
+ // Now sectioning phase: we allow both [a_lo, a_hi] and [a_hi, a_lo]
+ while (i < max_evaluations)
+ {
+ const NumberType a_lo_safe = a_lo + tau2 * (a_hi - a_lo);
+ const NumberType a_hi_safe = a_hi - tau3 * (a_hi - a_lo);
+ const auto bounds = std::minmax(a_lo_safe, a_hi_safe);
+
+ ai = choose(
+ a_lo, f_lo, g_lo, a_hi, f_hi, g_hi, a_rec, f_rec, g_rec, bounds);
+
+ const std::pair<NumberType, NumberType> fgi = func(ai);
+ fi = fgi.first;
+ gi = fgi.second;
+ i++;
+
+ if (debug_output)
+ deallog << "Sectioning phase: " << i << std::endl
+ << a_lo << " " << f_lo << " " << g_lo << std::endl
+ << a_hi << " " << f_hi << " " << g_hi << std::endl
+ << ai << " " << fi << " " << gi << std::endl;
+
+ if (!w1(ai, fi) || fi >= f_lo)
+ // take [a_lo, ai]
+ {
+ a_rec.add(a_hi);
+ f_rec.add(f_hi);
+ g_rec.add(g_hi);
+
+ a_hi = ai;
+ f_hi = fi;
+ g_hi = gi;
+ }
+ else
+ {
+ if (w2(gi)) // satisfies both wolf
+ {
+ Assert(w1(ai, fi), ExcInternalError());
+ return std::make_pair(ai, i);
+ }
+
+ if (gi * (a_hi - a_lo) >= 0)
+ // take [ai, a_lo]
+ {
+ a_rec.add(a_hi);
+ f_rec.add(f_hi);
+ g_rec.add(g_hi);
+
+ a_hi = a_lo;
+ f_hi = f_lo;
+ g_hi = g_lo;
+ }
+ else
+ // take [ai, a_hi]
+ {
+ a_rec.add(a_lo);
+ f_rec.add(f_lo);
+ g_rec.add(g_lo);
+ }
+
+ a_lo = ai;
+ f_lo = fi;
+ g_lo = gi;
+ }
+ }
+
+ // if we got here, we could not find the solution
+ AssertThrow(
+ false,
+ ExcMessage(
+ "Could not could complete the sectioning phase within the given number of iterations."));
+ return std::make_pair(std::numeric_limits<NumberType>::signaling_NaN(), i);
+ }
+
+#endif
+
+} // namespace LineMinimization
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // dealii_line_minimization_h
--- /dev/null
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+INCLUDE(../setup_testsubproject.cmake)
+PROJECT(testsuite CXX)
+DEAL_II_PICKUP_TESTS()
+
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+
+// check minimization of the cubic fit based on f(x1), f(x2) and f'(x1)
+// and f'(x2)
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/optimization/line_minimization.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+void
+test()
+{
+ // test 1:
+ {
+ auto f = [](double x) {
+ return std::pow(x, 4) - 20. * std::pow(x, 3) + 0.1 * x;
+ };
+ auto g = [](double x) {
+ return 4. * std::pow(x, 3) - 60. * std::pow(x, 2) + 0.1;
+ };
+
+ const double x1 = 5;
+ const double x2 = 17;
+ const double f1 = f(x1);
+ const double f2 = f(x2);
+ const double g1 = g(x1);
+ const double g2 = g(x2);
+ const double res = *LineMinimization::cubic_fit(x1, f1, g1, x2, f2, g2);
+ const double res2 = *LineMinimization::cubic_fit(x2, f2, g2, x1, f1, g1);
+ deallog << x1 << " " << x2 << std::endl
+ << f1 << " " << f2 << std::endl
+ << g1 << " " << g2 << std::endl
+ << res << std::endl
+ << res2 << std::endl;
+ }
+}
+
+
+int
+main(int argc, char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile, /*do not print job id*/ false);
+ deallog.depth_console(0);
+
+ test();
+}
--- /dev/null
+DEAL::5.00000 17.0000
+DEAL::-1874.50 -14737.3
+DEAL::-999.900 2312.10
+DEAL::14.6115
+DEAL::14.6115
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+
+// check minimization of the cubic fit based on f(x1), f(x2) and f'(x1)
+// and f'(x2)
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/optimization/line_minimization.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+void
+test()
+{
+ // test 1:
+ {
+ auto f = [](double x) {
+ return std::pow(x, 4) - 20. * std::pow(x, 3) + 0.1 * x;
+ };
+ auto g = [](double x) {
+ return 4. * std::pow(x, 3) - 60. * std::pow(x, 2) + 0.1;
+ };
+
+ const double x1 = 17;
+ const double x2 = 10;
+ const double x3 = 5;
+ const double f1 = f(x1);
+ const double f2 = f(x2);
+ const double f3 = f(x3);
+ const double g1 = g(x1);
+ const double res =
+ *LineMinimization::cubic_fit_three_points(x1, f1, g1, x2, f2, x3, f3);
+ deallog << x1 << " " << f1 << " " << g1 << std::endl
+ << x2 << " " << f2 << std::endl
+ << x3 << " " << f3 << std::endl
+ << res << std::endl;
+ }
+}
+
+
+int
+main(int argc, char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile, /*do not print job id*/ false);
+ deallog.depth_console(0);
+
+ test();
+}
--- /dev/null
+DEAL::17.0000 -14737.3 2312.10
+DEAL::10.0000 -9999.00
+DEAL::5.00000 -1874.50
+DEAL::14.8441
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+
+// check line minimization with strong Wolfe conditions
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/optimization/line_minimization.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+using namespace dealii;
+
+/*
+ * MWE in Maxima
+
+func(x):=100 * x^4 + (1-x)^2;
+gfunc(x):=''(diff(func(x),x));
+eta : 0.1;
+mu : 0.01;
+w1(x) := func(x) - func(0) - x * mu * gfunc(0);
+w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0));
+plot2d([func(x),gfunc(x),w1(x), w2(x)], [x,0.1,0.2]);
+w1(0.159668);
+w2(0.159668);
+bfloat(solve(gfunc(x)=0)[3]);
+
+ */
+
+
+void
+test()
+{
+ // test 1:
+ {
+ const double min_x = 0.161262023139589;
+ auto func = [](const double x) {
+ const double f = 100. * std::pow(x, 4) + std::pow(1. - x, 2);
+ const double g = 400. * std::pow(x, 3) - 2. * (1. - x);
+ return std::make_pair(f, g);
+ };
+
+ const auto fg0 = func(0);
+ const auto res =
+ LineMinimization::line_search<double>(func,
+ fg0.first,
+ fg0.second,
+ LineMinimization::poly_fit<double>,
+ 0.1,
+ 0.1,
+ 0.01,
+ 100,
+ 20,
+ true);
+ deallog << "Solution: " << res.first << std::endl
+ << "Distance: " << std::fabs(res.first - min_x) << std::endl;
+ }
+}
+
+
+int
+main(int argc, char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile, /*do not print job id*/ false);
+ deallog.depth_console(0);
+
+ test();
+}
--- /dev/null
+DEAL::Bracketing phase: 1
+DEAL::0.100000 0.820000 -1.40000
+DEAL::Bracketing phase: 2
+DEAL::1.00000 100.000 400.000
+DEAL::Sectioning phase: 3
+DEAL::0.100000 0.820000 -1.40000
+DEAL::1.00000 100.000 400.000
+DEAL::0.333333 1.67901 13.4815
+DEAL::Sectioning phase: 4
+DEAL::0.100000 0.820000 -1.40000
+DEAL::0.333333 1.67901 13.4815
+DEAL::0.167641 0.771802 0.219785
+DEAL::Sectioning phase: 5
+DEAL::0.167641 0.771802 0.219785
+DEAL::0.100000 0.820000 -1.40000
+DEAL::0.159668 0.771152 -0.0524435
+DEAL::Solution: 0.159668
+DEAL::Distance: 0.00159407
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+
+// check line minimization with strong Wolfe conditions
+// similar to line_minimization.cc but a different function and
+// different initial steps
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/optimization/line_minimization.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+using namespace dealii;
+
+/*
+ * MWE in Maxima
+
+func(x):=-3*x/(x^2+2)-0.03*x;
+gfunc(x):=''(diff(func(x),x));
+mu:0.025;
+eta : 0.1;
+a_max : 30;
+f_min : func(0) + a_max * mu * gfunc(0);
+w1(x) := func(x) - func(0) - x * mu * gfunc(0);
+w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0));
+w(x) := signum( signum(w1(x)) + signum(w2(x)) + 1 );
+pline(x):=func(0)+x*mu*gfunc(0);
+plot2d([func(x),pline(x),w(x)], [x,0,20]);
+bfloat(solve(gfunc(x)=0)[4]);
+
+ */
+
+
+void
+test()
+{
+ // test 1:
+ {
+ const double min_x = 1.474531468108294;
+ auto func = [](const double x) {
+ const double f = (-3. * x) / (x * x + 2.) - 0.03 * x;
+ const double g =
+ -3. / (x * x + 2.) + (6. * x * x) / std::pow(x * x + 2., 2) - 0.03;
+ return std::make_pair(f, g);
+ };
+
+ const auto fg0 = func(0);
+
+ {
+ deallog << "Case 1:" << std::endl;
+ // First, overshoot and get to solution immediately
+ const auto res = LineMinimization::line_search<double>(
+ func,
+ fg0.first,
+ fg0.second,
+ LineMinimization::poly_fit<double>,
+ 13,
+ 0.1,
+ 0.025,
+ 30,
+ 20,
+ true);
+ deallog << "Solution: " << res.first << std::endl
+ << "Distance: " << std::fabs(res.first - min_x) << std::endl;
+ }
+
+ {
+ deallog << "Case 2:" << std::endl;
+ // Now a small step to converge where needed:
+ const auto res = LineMinimization::line_search<double>(
+ func,
+ fg0.first,
+ fg0.second,
+ LineMinimization::poly_fit<double>,
+ 0.1,
+ 0.1,
+ 0.025,
+ 30,
+ 20,
+ true);
+ deallog << "Solution: " << res.first << std::endl
+ << "Distance: " << std::fabs(res.first - min_x) << std::endl;
+ }
+
+ {
+ deallog << "Case 3:" << std::endl;
+ // Now do a big step so that next one in bracketing satisfies both Wolf:
+ // at the termination point the derivative is alos negative!
+ // Also that point satisfies both Wolfe conditions as well, but
+ // we are interested in another segment, which contains local
+ // minimizer
+ const auto res = LineMinimization::line_search<double>(
+ func,
+ fg0.first,
+ fg0.second,
+ LineMinimization::poly_fit<double>,
+ 1,
+ 0.1,
+ 0.025,
+ 30,
+ 20,
+ true);
+ deallog << "Solution: " << res.first << std::endl
+ << "Distance: " << std::fabs(res.first - min_x) << std::endl;
+ }
+ }
+}
+
+
+int
+main(int argc, char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile, /*do not print job id*/ false);
+ deallog.depth_console(0);
+
+ test();
+}
--- /dev/null
+DEAL::Case 1:
+DEAL::Bracketing phase: 1
+DEAL::13.0000 -0.618070 -0.0128665
+DEAL::Solution: 13.0000
+DEAL::Distance: 11.5255
+DEAL::Case 2:
+DEAL::Bracketing phase: 1
+DEAL::0.100000 -0.152254 -1.50769
+DEAL::Bracketing phase: 2
+DEAL::1.00000 -1.03000 -0.363333
+DEAL::Bracketing phase: 3
+DEAL::9.10000 -0.594896 0.00370484
+DEAL::Sectioning phase: 4
+DEAL::1.00000 -1.03000 -0.363333
+DEAL::9.10000 -0.594896 0.00370484
+DEAL::3.09290 -0.895024 0.139676
+DEAL::Sectioning phase: 5
+DEAL::1.00000 -1.03000 -0.363333
+DEAL::3.09290 -0.895024 0.139676
+DEAL::1.60613 -1.10031 0.0529144
+DEAL::Solution: 1.60613
+DEAL::Distance: 0.131601
+DEAL::Case 3:
+DEAL::Bracketing phase: 1
+DEAL::1.00000 -1.03000 -0.363333
+DEAL::Bracketing phase: 2
+DEAL::10.0000 -0.594118 -0.00174164
+DEAL::Sectioning phase: 3
+DEAL::1.00000 -1.03000 -0.363333
+DEAL::10.0000 -0.594118 -0.00174164
+DEAL::3.36365 -0.858821 0.127629
+DEAL::Sectioning phase: 4
+DEAL::1.00000 -1.03000 -0.363333
+DEAL::3.36365 -0.858821 0.127629
+DEAL::1.65166 -1.09756 0.0676998
+DEAL::Solution: 1.65166
+DEAL::Distance: 0.177131
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+
+// check line minimization with strong Wolfe conditions
+// use functions from More Thuente 1994, Line Search Algorithms with
+// Guaranteed Sufficient Decrease
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/optimization/line_minimization.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+using namespace dealii;
+
+/*
+ * MWE in Maxima
+
+
+Case 1 (compare to Table I) Function 5.1 Figure 3:
+
+b:2;
+func(x):=-x/(x^2+b);
+gfunc(x):=''(diff(func(x),x));
+mu:0.001;
+eta : 0.1;
+w1(x) := func(x) - func(0) - x * mu * gfunc(0);
+w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0));
+w(x) := signum( signum(w1(x)) + signum(w2(x)) + 1 );
+pline(x):=func(0)+x*mu*gfunc(0);
+plot2d([func(x),pline(x),w(x)], [x,0,16]);
+
+Case 2, Function 5.2 Figure 4 (compared to Table II we have more iterations, up
+to 8 vs 16):
+
+b:0.004;
+func(x):=(x+b)^5-2*(x+b)^4;
+gfunc(x):=''(diff(func(x),x));
+mu:0.1;
+eta : 0.1;
+w1(x) := func(x) - func(0) - x * mu * gfunc(0);
+w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0));
+w(x) := signum( signum(w1(x)) + signum(w2(x)) + 1 );
+pline(x):=func(0)+x*mu*gfunc(0);
+plot2d([func(x),pline(x),w(x)], [x,0,2]);
+
+Case 3-5 Function 5.4 Figure 6 (Table IV-V-VI):
+
+b1:0.001;
+b2:0.001;
+
+b1:0.01;
+b2:0.001;
+
+b1:0.001;
+b2:0.01;
+
+g(x):=(1+x^2)^(1/2)-x;
+func(x):=g(b1)*((1-x)^2+b2^2)^(1/2) + g(b2)*(x^2+b1^2)^(1/2);
+gfunc(x):=''(diff(func(x),x));
+plot2d([func(x)], [x,0,1]);
+
+ */
+
+
+void
+test()
+{
+ const std::vector<double> values = {{1e-3, 1e-1, 1e+1, 1e+3}};
+ double f0, g0, fi, gi;
+
+ {
+ deallog << "Table 1:" << std::endl;
+ const double b = 2;
+ auto func = [&](const double &x) {
+ const double f = -x / (x * x + b);
+ const double g = 2. * x * x / std::pow(x * x + 2., 2) - 1. / (x * x + 2.);
+ return std::make_pair(f, g);
+ };
+
+ const auto fg0 = func(0);
+
+ for (auto a1 : values)
+ {
+ const auto res = LineMinimization::line_search<double>(
+ func,
+ fg0.first,
+ fg0.second,
+ LineMinimization::poly_fit<double>,
+ a1,
+ 0.1,
+ 0.001);
+
+ const auto fgi = func(res.first);
+ deallog << res.second << " " << res.first << " " << fgi.second
+ << std::endl;
+ }
+ }
+
+ {
+ deallog << "Table 2:" << std::endl;
+ const double b = 0.004;
+ auto func = [&](const double x) {
+ const double f = std::pow(x + b, 5) - 2. * std::pow(x + b, 4);
+ const double g = 5. * std::pow(x + b, 4) - 8. * std::pow(x + b, 3);
+ return std::make_pair(f, g);
+ };
+
+ const auto fg0 = func(0);
+
+ for (auto a1 : values)
+ {
+ const auto res = LineMinimization::line_search<double>(
+ func,
+ fg0.first,
+ fg0.second,
+ LineMinimization::poly_fit<double>,
+ a1,
+ 0.100001,
+ 0.1);
+
+ const auto fgi = func(res.first);
+ deallog << res.second << " " << res.first << " " << fgi.second
+ << std::endl;
+ }
+ }
+
+ {
+ const std::vector<std::pair<double, double>> params = {
+ {{0.001, 0.001}, {0.01, 0.001}, {0.001, 0.01}}};
+
+ unsigned int ind = 4;
+ for (auto p : params)
+ {
+ deallog << "Table " << ind++ << ":" << std::endl;
+ const double b1 = p.first;
+ const double b2 = p.second;
+
+ const double gb1 = std::sqrt(1. + b1 * b1) - b1;
+ const double gb2 = std::sqrt(1. + b2 * b2) - b2;
+
+ auto func = [&](const double x) {
+ const double f = gb1 * std::sqrt(std::pow(1. - x, 2) + b2 * b2) +
+ gb2 * std::sqrt(x * x + b1 * b1);
+ const double g =
+ gb2 * x / sqrt(x * x + b1 * b1) -
+ gb1 * (1. - x) / std::sqrt(std::pow(1 - x, 2) + b2 * b2);
+ return std::make_pair(f, g);
+ };
+
+ const auto fg0 = func(0);
+
+ for (auto a1 : values)
+ {
+ const auto res = LineMinimization::line_search<double>(
+ func,
+ fg0.first,
+ fg0.second,
+ LineMinimization::poly_fit<double>,
+ a1,
+ 0.00100001,
+ 0.001);
+
+ const auto fgi = func(res.first);
+ deallog << res.second << " " << res.first << " " << fgi.second
+ << std::endl;
+ }
+ }
+ }
+}
+
+
+int
+main(int argc, char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile, /*do not print job id*/ false);
+ deallog.depth_console(0);
+
+ test();
+}
--- /dev/null
+DEAL::Table 1:
+DEAL::7 1.47453 0.00999934
+DEAL::5 1.51075 0.0153974
+DEAL::1 10.0000 0.00941945
+DEAL::4 37.0539 0.000725162
+DEAL::Table 2:
+DEAL::19 1.59600 -2.81091e-10
+DEAL::16 1.59600 -9.10774e-11
+DEAL::13 1.59600 3.23723e-10
+DEAL::16 1.59600 2.57199e-09
+DEAL::Table 4:
+DEAL::3 0.0910000 -5.97089e-05
+DEAL::1 0.100000 -4.93296e-05
+DEAL::3 0.371819 -2.34721e-06
+DEAL::5 0.545285 7.35858e-07
+DEAL::Table 5:
+DEAL::5 0.0743487 3.55960e-05
+DEAL::3 0.0744514 5.98435e-05
+DEAL::10 0.0733164 -0.000213775
+DEAL::13 0.0767649 0.000581084
+DEAL::Table 6:
+DEAL::12 0.925705 -2.28542e-05
+DEAL::11 0.923466 -0.000531186
+DEAL::12 0.925750 -1.23021e-05
+DEAL::9 0.924924 -0.000205208
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+
+// same as line_minimization_03 but use three points cubic fit.
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/optimization/line_minimization.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+using namespace dealii;
+
+/*
+ * MWE in Maxima
+
+
+Case 1 (compare to Table I) Function 5.1 Figure 3:
+
+b:2;
+func(x):=-x/(x^2+b);
+gfunc(x):=''(diff(func(x),x));
+mu:0.001;
+eta : 0.1;
+w1(x) := func(x) - func(0) - x * mu * gfunc(0);
+w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0));
+w(x) := signum( signum(w1(x)) + signum(w2(x)) + 1 );
+pline(x):=func(0)+x*mu*gfunc(0);
+plot2d([func(x),pline(x),w(x)], [x,0,16]);
+
+Case 2, Function 5.2 Figure 4 (compared to Table II we have more iterations, up
+to 8 vs 16):
+
+b:0.004;
+func(x):=(x+b)^5-2*(x+b)^4;
+gfunc(x):=''(diff(func(x),x));
+mu:0.1;
+eta : 0.1;
+w1(x) := func(x) - func(0) - x * mu * gfunc(0);
+w2(x) := abs(gfunc(x)) - eta * abs(gfunc(0));
+w(x) := signum( signum(w1(x)) + signum(w2(x)) + 1 );
+pline(x):=func(0)+x*mu*gfunc(0);
+plot2d([func(x),pline(x),w(x)], [x,0,2]);
+
+Case 3-5 Function 5.4 Figure 6 (Table IV-V-VI):
+
+b1:0.001;
+b2:0.001;
+
+b1:0.01;
+b2:0.001;
+
+b1:0.001;
+b2:0.01;
+
+g(x):=(1+x^2)^(1/2)-x;
+func(x):=g(b1)*((1-x)^2+b2^2)^(1/2) + g(b2)*(x^2+b1^2)^(1/2);
+gfunc(x):=''(diff(func(x),x));
+plot2d([func(x)], [x,0,1]);
+
+ */
+
+
+void
+test()
+{
+ const std::vector<double> values = {{1e-3, 1e-1, 1e+1, 1e+3}};
+ double f0, g0, fi, gi;
+
+ {
+ deallog << "Table 1:" << std::endl;
+ const double b = 2;
+ auto func = [&](const double x) {
+ const double f = -x / (x * x + b);
+ const double g = 2. * x * x / std::pow(x * x + 2., 2) - 1. / (x * x + 2.);
+ return std::make_pair(f, g);
+ };
+
+ const auto fg0 = func(0);
+
+ for (auto a1 : values)
+ {
+ const auto res = LineMinimization::line_search<double>(
+ func,
+ fg0.first,
+ fg0.second,
+ LineMinimization::poly_fit_three_points<double>,
+ a1,
+ 0.1,
+ 0.001);
+
+ const auto fgi = func(res.first);
+ deallog << res.second << " " << res.first << " " << fgi.second
+ << std::endl;
+ }
+ }
+
+ {
+ deallog << "Table 2:" << std::endl;
+ const double b = 0.004;
+ auto func = [&](const double x) {
+ const double f = std::pow(x + b, 5) - 2. * std::pow(x + b, 4);
+ const double g = 5. * std::pow(x + b, 4) - 8. * std::pow(x + b, 3);
+ return std::make_pair(f, g);
+ };
+
+ const auto fg0 = func(0);
+
+ for (auto a1 : values)
+ {
+ const auto res = LineMinimization::line_search<double>(
+ func,
+ fg0.first,
+ fg0.second,
+ LineMinimization::poly_fit_three_points<double>,
+ a1,
+ 0.100001,
+ 0.1);
+
+ const auto fgi = func(res.first);
+ deallog << res.second << " " << res.first << " " << fgi.second
+ << std::endl;
+ }
+ }
+
+ {
+ const std::vector<std::pair<double, double>> params = {
+ {{0.001, 0.001}, {0.01, 0.001}, {0.001, 0.01}}};
+
+ unsigned int ind = 4;
+ for (auto p : params)
+ {
+ deallog << "Table " << ind++ << ":" << std::endl;
+ const double b1 = p.first;
+ const double b2 = p.second;
+
+ const double gb1 = std::sqrt(1. + b1 * b1) - b1;
+ const double gb2 = std::sqrt(1. + b2 * b2) - b2;
+
+ auto func = [&](const double x) {
+ const double f = gb1 * std::sqrt(std::pow(1. - x, 2) + b2 * b2) +
+ gb2 * std::sqrt(x * x + b1 * b1);
+ const double g =
+ gb2 * x / sqrt(x * x + b1 * b1) -
+ gb1 * (1. - x) / std::sqrt(std::pow(1 - x, 2) + b2 * b2);
+ return std::make_pair(f, g);
+ };
+
+ const auto fg0 = func(0);
+
+ for (auto a1 : values)
+ {
+ const auto res = LineMinimization::line_search<double>(
+ func,
+ fg0.first,
+ fg0.second,
+ LineMinimization::poly_fit_three_points<double>,
+ a1,
+ 0.00100001,
+ 0.001);
+
+ const auto fgi = func(res.first);
+ deallog << res.second << " " << res.first << " " << fgi.second
+ << std::endl;
+ }
+ }
+ }
+}
+
+
+int
+main(int argc, char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile, /*do not print job id*/ false);
+ deallog.depth_console(0);
+
+ test();
+}
--- /dev/null
+DEAL::Table 1:
+DEAL::8 1.49141 0.0125702
+DEAL::6 1.49995 0.0138338
+DEAL::1 10.0000 0.00941945
+DEAL::5 25.6769 0.00150302
+DEAL::Table 2:
+DEAL::15 1.59600 2.94690e-10
+DEAL::19 1.59600 9.73245e-10
+DEAL::11 1.59600 -2.83316e-08
+DEAL::19 1.59600 -7.43938e-12
+DEAL::Table 4:
+DEAL::3 0.0910000 -5.97089e-05
+DEAL::1 0.100000 -4.93296e-05
+DEAL::3 0.431189 -1.14275e-06
+DEAL::5 0.693314 4.27147e-06
+DEAL::Table 5:
+DEAL::6 0.0741039 -2.26077e-05
+DEAL::3 0.0754771 0.000296756
+DEAL::11 0.0734116 -0.000190322
+DEAL::14 0.0723147 -0.000465817
+DEAL::Table 6:
+DEAL::13 0.924531 -0.000294927
+DEAL::11 0.925274 -0.000124134
+DEAL::6 0.926593 0.000191539
+DEAL::10 0.928406 0.000653661
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+
+// check minimization of the quadratic fit based on f(x1), f(x2) and f'(x1)
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/optimization/line_minimization.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+using namespace dealii;
+
+
+void
+test()
+{
+ // test 1:
+ {
+ auto f = [](double x) {
+ return std::pow(x, 4) - 20. * std::pow(x, 3) + 0.1 * x;
+ };
+ auto g = [](double x) {
+ return 4. * std::pow(x, 3) - 60. * std::pow(x, 2) + 0.1;
+ };
+
+ const double x1 = 10;
+ const double x2 = 17;
+ const double f1 = f(x1);
+ const double f2 = f(x2);
+ const double g1 = g(x1);
+ const double g2 = g(x2);
+ const double res = *LineMinimization::quadratic_fit(x1, f1, g1, x2, f2);
+ deallog << x1 << " " << x2 << std::endl
+ << f1 << " " << f2 << std::endl
+ << g1 << " " << g2 << std::endl
+ << res << std::endl;
+ }
+}
+
+
+int
+main(int argc, char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile, /*do not print job id*/ false);
+ deallog.depth_console(0);
+
+ test();
+}
--- /dev/null
+DEAL::10.0000 17.0000
+DEAL::-9999.00 -14737.3
+DEAL::-1999.90 2312.10
+DEAL::15.2907
projects / dealii.git / commitdiff