]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add line minimization functions
authorDenis Davydov <davydden@gmail.com>
Thu, 13 Sep 2018 10:14:53 +0000 (12:14 +0200)
committerDenis Davydov <davydden@gmail.com>
Mon, 17 Sep 2018 15:35:32 +0000 (17:35 +0200)
17 files changed:
doc/news/changes/minor/20180913DenisDavydov [new file with mode: 0644]
include/deal.II/optimization/line_minimization.h [new file with mode: 0644]
tests/optimization/CMakeLists.txt [new file with mode: 0644]
tests/optimization/cubic_fit.cc [new file with mode: 0644]
tests/optimization/cubic_fit.output [new file with mode: 0644]
tests/optimization/cubic_fit_three_points.cc [new file with mode: 0644]
tests/optimization/cubic_fit_three_points.output [new file with mode: 0644]
tests/optimization/line_minimization.cc [new file with mode: 0644]
tests/optimization/line_minimization.output [new file with mode: 0644]
tests/optimization/line_minimization_02.cc [new file with mode: 0644]
tests/optimization/line_minimization_02.output [new file with mode: 0644]
tests/optimization/line_minimization_03.cc [new file with mode: 0644]
tests/optimization/line_minimization_03.output [new file with mode: 0644]
tests/optimization/line_minimization_03b.cc [new file with mode: 0644]
tests/optimization/line_minimization_03b.output [new file with mode: 0644]
tests/optimization/quadratic_fit.cc [new file with mode: 0644]
tests/optimization/quadratic_fit.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20180913DenisDavydov b/doc/news/changes/minor/20180913DenisDavydov
new file mode 100644 (file)
index 0000000..31324ec
--- /dev/null
@@ -0,0 +1,3 @@
+New: Add line minimization functions.
+<br>
+(Denis Davydov 2018/09/13)
diff --git a/include/deal.II/optimization/line_minimization.h b/include/deal.II/optimization/line_minimization.h
new file mode 100644 (file)
index 0000000..160913b
--- /dev/null
@@ -0,0 +1,606 @@
+//-----------------------------------------------------------
+//
+//    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
diff --git a/tests/optimization/CMakeLists.txt b/tests/optimization/CMakeLists.txt
new file mode 100644 (file)
index 0000000..3772fac
--- /dev/null
@@ -0,0 +1,5 @@
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+INCLUDE(../setup_testsubproject.cmake)
+PROJECT(testsuite CXX)
+DEAL_II_PICKUP_TESTS()
+
diff --git a/tests/optimization/cubic_fit.cc b/tests/optimization/cubic_fit.cc
new file mode 100644 (file)
index 0000000..742fa70
--- /dev/null
@@ -0,0 +1,67 @@
+//-----------------------------------------------------------
+//
+//    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();
+}
diff --git a/tests/optimization/cubic_fit.output b/tests/optimization/cubic_fit.output
new file mode 100644 (file)
index 0000000..b0ad696
--- /dev/null
@@ -0,0 +1,5 @@
+DEAL::5.00000 17.0000
+DEAL::-1874.50 -14737.3
+DEAL::-999.900 2312.10
+DEAL::14.6115
+DEAL::14.6115
diff --git a/tests/optimization/cubic_fit_three_points.cc b/tests/optimization/cubic_fit_three_points.cc
new file mode 100644 (file)
index 0000000..0f970ec
--- /dev/null
@@ -0,0 +1,67 @@
+//-----------------------------------------------------------
+//
+//    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();
+}
diff --git a/tests/optimization/cubic_fit_three_points.output b/tests/optimization/cubic_fit_three_points.output
new file mode 100644 (file)
index 0000000..fa97885
--- /dev/null
@@ -0,0 +1,4 @@
+DEAL::17.0000 -14737.3 2312.10
+DEAL::10.0000 -9999.00
+DEAL::5.00000 -1874.50
+DEAL::14.8441
diff --git a/tests/optimization/line_minimization.cc b/tests/optimization/line_minimization.cc
new file mode 100644 (file)
index 0000000..8dfb886
--- /dev/null
@@ -0,0 +1,83 @@
+//-----------------------------------------------------------
+//
+//    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();
+}
diff --git a/tests/optimization/line_minimization.output b/tests/optimization/line_minimization.output
new file mode 100644 (file)
index 0000000..e1a7ffb
--- /dev/null
@@ -0,0 +1,18 @@
+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
diff --git a/tests/optimization/line_minimization_02.cc b/tests/optimization/line_minimization_02.cc
new file mode 100644 (file)
index 0000000..cc49a4f
--- /dev/null
@@ -0,0 +1,133 @@
+//-----------------------------------------------------------
+//
+//    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();
+}
diff --git a/tests/optimization/line_minimization_02.output b/tests/optimization/line_minimization_02.output
new file mode 100644 (file)
index 0000000..5a8e71a
--- /dev/null
@@ -0,0 +1,37 @@
+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
diff --git a/tests/optimization/line_minimization_03.cc b/tests/optimization/line_minimization_03.cc
new file mode 100644 (file)
index 0000000..7774bf6
--- /dev/null
@@ -0,0 +1,195 @@
+//-----------------------------------------------------------
+//
+//    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();
+}
diff --git a/tests/optimization/line_minimization_03.output b/tests/optimization/line_minimization_03.output
new file mode 100644 (file)
index 0000000..688f24c
--- /dev/null
@@ -0,0 +1,25 @@
+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
diff --git a/tests/optimization/line_minimization_03b.cc b/tests/optimization/line_minimization_03b.cc
new file mode 100644 (file)
index 0000000..8e19d3b
--- /dev/null
@@ -0,0 +1,193 @@
+//-----------------------------------------------------------
+//
+//    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();
+}
diff --git a/tests/optimization/line_minimization_03b.output b/tests/optimization/line_minimization_03b.output
new file mode 100644 (file)
index 0000000..43e82b9
--- /dev/null
@@ -0,0 +1,25 @@
+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
diff --git a/tests/optimization/quadratic_fit.cc b/tests/optimization/quadratic_fit.cc
new file mode 100644 (file)
index 0000000..8183bb8
--- /dev/null
@@ -0,0 +1,64 @@
+//-----------------------------------------------------------
+//
+//    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();
+}
diff --git a/tests/optimization/quadratic_fit.output b/tests/optimization/quadratic_fit.output
new file mode 100644 (file)
index 0000000..8537015
--- /dev/null
@@ -0,0 +1,4 @@
+DEAL::10.0000 17.0000
+DEAL::-9999.00 -14737.3
+DEAL::-1999.90 2312.10
+DEAL::15.2907

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.