#ifndef MINIMISATION_HH
#define MINIMISATION_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 "circularconvexfunction.hh"
#include "mydirectionalconvexfunction.hh"

namespace Dune {
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);
    }
  }
}
}
#endif