bisection-simpler-example2-gradient.cc

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

// FIXME: is there a bette way to do this?
//#define DUNE_MINIMAL_DEBUG_LEVEL 2
#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 <limits>

template <int dimension>
class SampleFunctional : public SmallFunctional<dimension> {
public:
  typedef typename SmallFunctional<dimension>::SmallMatrix SmallMatrix;
  typedef typename SmallFunctional<dimension>::SmallVector SmallVector;

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

  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; // <1/2 Av - b,v>
  }

  SmallVector d(const SmallVector v) const {
    SmallVector y;
    A_.mv(v, y); // y = Av
    y -= b_;     // y = Av - b
    return y;
  }

  SmallMatrix d2(const SmallVector v) const { return A_; }

  // extending the interface here
  double directionalDerivative(const SmallVector x,
                               const SmallVector dir) const {
    return d(x) * dir;
  }

  SmallVector minimise(const SmallVector x, int iterations) const {
    SmallVector descDir = d(x);
    descDir *= -1; // The negative gradient

    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 = descDir;
    tmp *= r;
    while (directionalDerivative(x + tmp, descDir) < 0) {
      l = r;
      r *= 2;
      tmp *= 2;
      Dune::dverb << "Widened interval!" << std::endl;
    }
    Dune::dverb << "Interval now [" << l << "," << r << "]" << std::endl;

    // Debugging
    {
      SmallVector tmpl = tmp;
      tmpl *= l;
      SmallVector tmpr = tmp;
      tmpr *= r;
      assert(directionalDerivative(x + tmpl, descDir) < 0);
      assert(directionalDerivative(x + tmpr, descDir) > 0);
    }

    double m = l / 2 + r / 2;
    SmallVector middle;
    for (size_t 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 derivative = directionalDerivative(x + middle, descDir);

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

private:
  SmallMatrix A_;
  SmallVector b_;
};

int main() {
  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);

  SampleFunctional::SmallVector start = b;
  start *= 17;
  SampleFunctional::SmallVector correction = J.minimise(start, 20);
  assert(J(start + correction) <= J(start));
  std::cout << J(start + correction) << std::endl;

  return 0;
}