bisection-example-flexible.cc

/* -*- mode:c++; mode: flymake -*- */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <dune/common/exceptions.hh>
#include <dune/common/stdstreams.hh>

#include <dune/fufem/interval.hh>

#include <dune/tnnmg/nonlinearities/smallfunctional.hh>
#include <dune/tnnmg/problem-classes/bisection.hh>

#include <cassert>
#include <cstdlib>
#include <limits>

#include "properscalarincreasingconvexfunction.hh"

template <int dimension, class Function = TrivialFunction>
class SampleFunctional {
public:
  typedef Dune::FieldVector<double, dimension> SmallVector;
  typedef Dune::FieldMatrix<double, dimension, dimension> SmallMatrix;

  SampleFunctional(SmallMatrix A, SmallVector b) : A_(A), b_(b), func_() {}

  double operator()(const SmallVector v) const {
    SmallVector y;
    A_.mv(v, y);                        // y = Av
    y /= 2;                             // y = 1/2 Av
    y -= b_;                            // y = 1/2 Av - b
    return y * v + func_(v.two_norm()); // <1/2 Av - b,v> + H(|v|)
  }

  double directionalDerivative(const SmallVector x,
                               const SmallVector dir) const {
    double norm = dir.two_norm();

    if (norm == 0)
      DUNE_THROW(Dune::Exception, "Directional derivatives cannot be computed "
                                  "w.r.t. the zero-direction.");

    SmallVector tmp = dir;
    tmp /= norm;

    return pseudoDirectionalDerivative(x, tmp);
  }

  SmallVector minimise(const SmallVector x, unsigned int iterations) const {
    SmallVector descDir = ModifiedGradient(x);
    if (descDir == SmallVector(0.0))
      return SmallVector(0.0);

    Dune::dverb << "Starting at x with J(x) = " << operator()(x) << std::endl;
    Dune::dverb << "Minimizing in direction w with dJ(x,w) = "
                << directionalDerivative(x, descDir) << std::endl;

    double l = 0;
    double r = 1;
    SmallVector tmp;
    while (true) {
      tmp = x;
      tmp.axpy(r, descDir);
      if (pseudoDirectionalDerivative(tmp, descDir) >= 0)
        break;

      l = r;
      r *= 2;
      Dune::dverb << "Widened interval!" << std::endl;
    }
    Dune::dverb << "Interval now [" << l << "," << r << "]" << std::endl;

#ifndef NDEBUG
    {
      SmallVector tmpl = x;
      tmpl.axpy(l, descDir);
      SmallVector tmpr = x;
      tmpr.axpy(r, descDir);
      assert(directionalDerivative(tmpl, descDir) < 0);
      assert(directionalDerivative(tmpr, descDir) > 0);
    }
#endif

    double m = l / 2 + r / 2;
    SmallVector middle = SmallVector(0.0);
    for (unsigned int count = 0; count < iterations; ++count) {
      Dune::dverb << "now at m = " << m << std::endl;
      Dune::dverb << "Value of J here: " << operator()(x + middle) << std::endl;

      middle = descDir;
      middle *= m;

      double pseudoDerivative =
          pseudoDirectionalDerivative(x + middle, descDir);

      if (pseudoDerivative < 0) {
        l = m;
        m = (m + r) / 2;
      } else if (pseudoDerivative > 0) {
        r = m;
        m = (l + m) / 2;
      } else
        break;
    }
    return middle;
  }

private:
  SmallMatrix A_;
  SmallVector b_;

  Function func_;

  // Gradient of the smooth part
  SmallVector SmoothGrad(const SmallVector x) const {
    SmallVector y;
    A_.mv(x, y); // y = Av
    y -= b_;     // y = Av - b
    return y;
  }

  SmallVector PlusGrad(const SmallVector x) const {
    SmallVector y = SmoothGrad(x);
    y.axpy(func_.rightDifferential(x.two_norm()) / x.two_norm(), x);
    return y;
  }

  SmallVector MinusGrad(const SmallVector x) const {
    SmallVector y = SmoothGrad(x);
    y.axpy(func_.leftDifferential(x.two_norm()) / x.two_norm(), x);
    return y;
  }

  // |dir|-times the directional derivative wrt dir/|dir|.  If only
  // the sign of the directionalDerivative matters, this saves the
  // cost of normalising.
  double pseudoDirectionalDerivative(const SmallVector x,
                                     const SmallVector dir) const {
    if (x == SmallVector(0.0))
      return func_.rightDifferential(0) * dir.two_norm();

    if (x * dir > 0)
      return PlusGrad(x) * dir;
    else
      return MinusGrad(x) * dir;
  }

  SmallVector ModifiedGradient(const SmallVector x) const {
    if (x == SmallVector(0.0)) {
      SmallVector d = SmoothGrad(x);
      // Decline of the smooth part in the negative gradient direction
      double smoothDecline = -(d * d);
      double nonlinearDecline =
          func_.rightDifferential(0.0) * d.two_norm(); // TODO: is this correct?
      double combinedDecline = smoothDecline + nonlinearDecline;

      return (combinedDecline < 0) ? d : SmallVector(0.0);
    }

    SmallVector const pg = PlusGrad(x);
    SmallVector const mg = MinusGrad(x);
    SmallVector ret;
    // TODO: collinearity checks suck
    if (pg * x == pg.two_norm() * x.two_norm() &&
        -(mg * x) == mg.two_norm() * x.two_norm()) {
      return SmallVector(0);
    } else if (pg * x >= 0 && mg * x >= 0) {
      ret = pg;
    } else if (pg * x <= 0 && mg * x <= 0) {
      ret = mg;
    } else {
      ret = project(SmoothGrad(x), x);
    }
    ret *= -1;
    return ret;
  }

  // No normalising is done!
  SmallVector project(const SmallVector z, const SmallVector x) const {
    SmallVector y = z;
    y.axpy(-(z * x) / x.two_norm2(), x);
    return y;
  }
};

void testSampleFunction() {
  int const dim = 2;
  typedef SampleFunctional<dim, SampleFunction> SampleFunctional;

  SampleFunctional::SmallMatrix A;
  A[0][0] = 3;
  A[0][1] = 0;
  A[1][0] = 0;
  A[1][1] = 3;
  SampleFunctional::SmallVector b;
  b[0] = 1;
  b[1] = 2;

  SampleFunctional J(A, b);

  std::cout << J.directionalDerivative(b, b) << std::endl;
  assert(J.directionalDerivative(b, b) == 2 * sqrt(5) + 2);

  SampleFunctional::SmallVector start = b;
  start *= 17;
  SampleFunctional::SmallVector correction = J.minimise(start, 20);
  assert(J(start + correction) <= J(start));
  assert(std::abs(J(start + correction) + 0.254644) < 1e-8);
  std::cout << J(start + correction) << std::endl;
}

void testTrivialFunction() {
  int const dim = 2;
  typedef SampleFunctional<dim> SampleFunctional;

  SampleFunctional::SmallMatrix A;
  A[0][0] = 3;
  A[0][1] = 0;
  A[1][0] = 0;
  A[1][1] = 3;
  SampleFunctional::SmallVector b;
  b[0] = 1;
  b[1] = 2;

  SampleFunctional J(A, b);

  std::cout << J.directionalDerivative(b, b) << std::endl;
  assert(J.directionalDerivative(b, b) == 2 * sqrt(5));

  SampleFunctional::SmallVector start = b;
  start *= 17;
  SampleFunctional::SmallVector correction = J.minimise(start, 20);
  assert(J(start + correction) <= J(start));
  assert(std::abs(J(start + correction) + 0.83333333) < 1e-8);
  std::cout << J(start + correction) << std::endl;
}

int main() {
  try {
    testSampleFunction();
    testTrivialFunction();
    return 0;
  }
  catch (Dune::Exception &e) {
    Dune::derr << "Dune reported error: " << e << std::endl;
  }
}