ellipticenergy.hh

#ifndef ELLIPTIC_ENERGY_HH
#define ELLIPTIC_ENERGY_HH

#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/stdstreams.hh>

#include <dune/fufem/interval.hh>
#include <dune/tnnmg/problem-classes/bisection.hh>

#include "mydirectionalconvexfunction.hh"
#include "localnonlinearity.hh"
#include "circularconvexfunction.hh"

namespace Dune {
template <int dim> class EllipticEnergy {
public:
  typedef FieldVector<double, dim> SmallVector;
  typedef FieldMatrix<double, dim, dim> SmallMatrix;

  typedef LocalNonlinearity<dim> NonlinearityType;

  EllipticEnergy(SmallMatrix const &A, SmallVector const &b,
                 shared_ptr<NonlinearityType const> phi, int ignore = dim)
      : A(A), b(b), phi(phi), ignore(ignore) {}

  double operator()(SmallVector const &v) const {
    SmallVector y(0);
    A.usmv(0.5, v, y);        //  1/2 Av
    y -= b;                   //  1/2 Av - b
    return y * v + (*phi)(v); // <1/2 Av - b,v> + H(|v|)
  }

  void descentAtZero(SmallVector &ret) const {
    SmallVector const zero(0);
    // If there is a direction of descent, this is it
    SmallVector d;
    smoothGradient(zero, d);
    d *= -1;

    Interval<double> D;
    phi->directionalSubDiff(zero, d, D);
    double const nonlinearDecline = D[1];
    double const smoothDecline = -(d * d);
    double const combinedDecline = smoothDecline + nonlinearDecline;

    if (combinedDecline < 0) {
      ret = d;
    } else {
      ret = 0;
    }
  }

  // returns false if the direction is tangential
  void descentDirectionNew(SmallVector const &x, SmallVector &ret) const {
    SmallVector pg;
    upperGradient(x, pg);
    SmallVector mg;
    lowerGradient(x, mg);
    double const pgx = pg * x;
    double const mgx = mg * x;
    if (pgx >= 0 && mgx >= 0) {
      ret = pg;
      ret *= -1;
    } else if (pgx <= 0 && mgx <= 0) {
      ret = mg;
      ret *= -1;
    } else {
      assert(false);
    }
  }

  bool descentDirection(SmallVector const &x, SmallVector &ret) const {
    // Check the squared norm rather than each component because
    // complementaryProjection() divides by it
    if (x.two_norm2() == 0.0) {
      descentAtZero(ret);
      return true;
    }

    SmallVector pg;
    upperGradient(x, pg);
    SmallVector mg;
    lowerGradient(x, mg);
    double const pgx = pg * x;
    double const mgx = mg * x;
    if (pgx >= 0 && mgx >= 0) {
      ret = pg;
      ret *= -1;
      return true;
    } else if (pgx <= 0 && mgx <= 0) {
      ret = mg;
      ret *= -1;
      return true;
    } else {
      // Includes the case that pg points in direction x and mg
      // points in direction -x. The projection will then be zero.
      SmallVector d;
      smoothGradient(x, d);
      complementaryProjection(d, x, ret);
      ret *= -1;
      return false;
    }
  }

  SmallMatrix const &A;
  SmallVector const &b;
  shared_ptr<NonlinearityType const> const phi;
  int const ignore; // Dimension that should be ignored; goes from 0
                    // to dim-1; the special value dim means that no
                    // dimension should be ignored

private:
  // Gradient of the smooth part
  void smoothGradient(SmallVector const &x, SmallVector &y) const {
    A.mv(x, y); // y = Av
    y -= b;     // y = Av - b
    if (ignore != dim)
      y[ignore] = 0;
  }

  void upperGradient(SmallVector const &x, SmallVector &y) const {
    phi->upperGradient(x, y);
    SmallVector z;
    smoothGradient(x, z);
    y += z;
    if (ignore != dim)
      y[ignore] = 0;
  }

  void lowerGradient(SmallVector const &x, SmallVector &y) const {
    phi->lowerGradient(x, y);
    SmallVector z;
    smoothGradient(x, z);
    y += z;
    if (ignore != dim)
      y[ignore] = 0;
  }

  // y = (id - P)(d) where P is the projection onto the line t*x
  void complementaryProjection(SmallVector const &d, SmallVector const &x,
                               SmallVector &y) const {
    double const dx = d * x;
    double const xx = x.two_norm2();
    y = d;
    y.axpy(-dx / xx, x);
  }
};

template <class Functional>
void descentMinimisation(Functional const &J,
                         typename Functional::SmallVector &x,
                         typename Functional::SmallVector const &descDir,
                         Bisection const &bisection) {
  typedef typename Functional::SmallVector SmallVector;
  typedef typename Functional::NonlinearityType LocalNonlinearityType;

  // {{{ Construct a restriction of J to the line x + t * descDir

  /* We have

     1/2 <A(u+xv),u+xv>-<b,u+xv> = 1/2 <Av,v> x^2 - <b-Au,v> x + <1/2 Au-b,u>

     since A is symmetric.
  */
  SmallVector tmp = J.b;               //  b
  J.A.mmv(x, tmp);                     //  b-Au
  double const JRestb = tmp * descDir; // <b-Au,v>

  J.A.mv(descDir, tmp);                //  Av
  double const JRestA = tmp * descDir; // <Av,v>

  MyDirectionalConvexFunction<LocalNonlinearityType> const JRest(
      JRestA, JRestb, *J.phi, x, descDir);
  // }}}

  { // Debug
    Interval<double> D;
    JRest.subDiff(0, D);

    dverb
        << "## Directional derivative (as per subdifferential of restriction): "
        << D[1] << " (coordinates of the restriction)" << std::endl;
    /*
      It is possible that this differs quite a lot from the
      directional derivative computed in the descentDirection()
      method:

      If phi is nonsmooth at x, so that the directional
      derivatives jump, and |x| is computed to be too small or too
      large globally or locally, the locally computed
      subdifferential and the globally computed subdifferential
      will no longer coincide!

      The assertion D[1] <= 0 may thus fail.
    */
  }

  int count;
  double const stepsize = bisection.minimize(JRest, 0.0, 0.0, count);
  dverb << "Number of iterations in the bisection method: " << count
        << std::endl;
  ;

  x.axpy(stepsize, descDir);
}

template <class Functional>
void tangentialMinimisation(Functional const &J,
                            typename Functional::SmallVector &x,
                            typename Functional::SmallVector const &descDir,
                            Bisection const &bisection) {
  typedef typename Functional::SmallMatrix SmallMatrix;
  typedef typename Functional::SmallVector SmallVector;

  // We completely ignore the nonlinearity here -- when restricted
  // to a circle, it just enters as a constant!
  CircularConvexFunction<SmallMatrix, SmallVector> const JRest(J.A, J.b, x,
                                                               descDir);

  int count;
  double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
  dverb << "Number of iterations in the bisection method: " << count
        << std::endl;
  ;

  // Since x is used in the computation of the rhs, do not write to it directly
  SmallVector tmp;
  JRest.cartesian(stepsize, tmp);
  x = tmp;
}

template <class Functional>
void minimise(Functional const &J, typename Functional::SmallVector &x,
              size_t steps, Bisection const &bisection) {
  typedef typename Functional::SmallVector SmallVector;

  for (size_t step = 0; step < steps; ++step) {
    SmallVector descDir;
    bool const linesearchp = J.descentDirection(x, descDir);

    if (descDir == SmallVector(0.0))
      return;

    if (linesearchp) {
      descentMinimisation(J, x, descDir, bisection);
    } else {
      Bisection slowBisection(bisection);
      slowBisection.setFastQuadratic(false);

      tangentialMinimisation(J, x, descDir, slowBisection);
    }
  }
}

// NOTE: only works for the 2D case
template <class Functional>
void minimise2(Functional const &J, typename Functional::SmallVector &x,
               size_t steps, Bisection const &bisection) {
  typedef typename Functional::SmallVector SmallVector;
  minimisationInitialiser(J, bisection, x);

  { // Make a single step where we already know we're not differentiable
    SmallVector descDir;
    if (x == SmallVector(0))
      J.descentAtZero(descDir);
    else
      J.descentDirectionNew(x, descDir);
    descentMinimisation(J, x, descDir, bisection);
  }

  // From here on, we're in a C1 region and stay there.
  for (size_t i = 1; i < steps; ++i) {
    SmallVector descDir;
    J.descentDirectionNew(x, descDir);
    descentMinimisation(J, x, descDir, bisection);
  }
}

template <class Functional>
void minimisationInitialiser(Functional const &J, Bisection const &bisection,
                             typename Functional::SmallVector &startingPoint) {
  typedef typename Functional::SmallVector SmallVector;

  std::vector<double> const kinks = { 5, 10,
                                      15 }; // FIXME: We happen to know these

  SmallVector x_old(0);
  double Jx_old = J(x_old);

  for (auto &kink : kinks) {
    SmallVector x1 = { kink, 0 }; // Random vector that has the right norm
    SmallVector const descDir1 = { x1[1], -x1[0] };
    tangentialMinimisation(J, x1, descDir1, bisection);
    double const Jx1 = J(x1);

    SmallVector x2 = { -x1[0], -x1[1] };
    SmallVector const descDir2 = { x2[1], -x2[0] };
    tangentialMinimisation(J, x2, descDir2, bisection);
    double const Jx2 = J(x2);

    double const Jx = (Jx2 < Jx1 ? Jx2 : Jx1);
    SmallVector const x = (Jx2 < Jx1 ? x2 : x1);

    if (Jx >= Jx_old) {
      startingPoint = x_old;
      return;
    }

    Jx_old = Jx;
    x_old = x;
  }
  startingPoint = x_old;
}
}
#endif