test-gradient-method.cc

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

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

#include <dune/common/exceptions.hh>

#include "samplefunctional.hh"

#include <cassert>

template <int dim, class Function>
double functionTester(
    Dune::SampleFunctional<dim, Function> J,
    typename Dune::SampleFunctional<dim, Function>::SmallVector &start,
    size_t runs) {
  typename Dune::SampleFunctional<dim, Function>::SmallVector correction;
  std::cout << "Old value: J(...) = " << J(start) << std::endl;
  for (size_t i = 1; i <= runs; ++i) {
    Dune::minimise(J, start, correction);
    start += correction;
    if (i != runs)
      std::cout << "New value: J(...) = " << J(start) << std::endl;
  }
  double const final = J(start);
  std::cout << "Final value J(...) = " << final << std::endl;
  return final;
}

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

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

  SampleFunctional J(A, b);

  SampleFunctional::SmallVector start = b;
  start *= 17;

  /*
    j(x)
    = Ax - b + 2/|x| x
    = 17*(6, 9.5) - (1, 2) + 2/sqrt(5) (1, 2)
    = (102 - 1 + 2/sqrt(5), 161.5 - 2 + 4/sqrt(5))
  */
  SampleFunctional::SmallVector error;
  error[0] = -(102 - 1 + 2 / sqrt(5));
  error[1] = -(161.5 - 2 + 4 / sqrt(5));
  SampleFunctional::SmallVector returned;
  J.descentDirection(start, returned);
  error -= returned;
  assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?

  double const ret1 = functionTester(J, start, 6);

  start[0] = 279;
  start[1] = -96;

  double const ret2 = functionTester(J, start, 10);
  assert(std::abs(ret1 - ret2) < 1e-5);
}

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

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

  SampleFunctional J(A, b);

  SampleFunctional::SmallVector start;
  SampleFunctional::SmallVector error;
  /*
    for x = b/|b|:

    j_(x)
    = Ax - b + 1/|x| x
    = 1/sqrt(5) (6, 9.5) - (1, 2) + 1/sqrt(5) (1, 2)
    = (7/sqrt(5) - 1, 11.5/sqrt(5) - 2)

    j+(x)
    = Ax - b + 2/|x| x
    = 1/sqrt(5) (6, 9.5) - (1, 2) + 2/sqrt(5) (1, 2)
    = (8/sqrt(5) - 1, 13.5/sqrt(5) - 2)
  */
  {
    start = b;
    start /= (b.two_norm() + 1e-12);
    assert(start.two_norm() < 1);

    SampleFunctional::SmallVector returned;
    J.descentDirection(start, returned);
    error[0] = -(7 / sqrt(5) - 1);
    error[1] = -(11.5 / sqrt(5) - 2);
    error -= returned;
    assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?

    functionTester(J, start, 6);
  }
  std::cout << std::endl;
  {
    start = b;
    start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1;
    assert(start.two_norm() > 1);

    SampleFunctional::SmallVector returned;
    J.descentDirection(start, returned);
    error[0] = -(8 / sqrt(5) - 1);
    error[1] = -(13.5 / sqrt(5) - 2);
    error -= returned;
    assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?

    functionTester(J, start, 6);
  }
}

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

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

  SampleFunctional J(A, b);

  SampleFunctional::SmallVector start = b;
  start *= 17;

  /*
    j(x)
    = Ax - b
    = 17*(6, 9.5) - (1, 2)
    = (102 - 1, 161.5 - 2)
  */
  SampleFunctional::SmallVector error;
  error[0] = -101;
  error[1] = -159.5;
  SampleFunctional::SmallVector returned;
  J.descentDirection(start, returned);
  error -= returned;
  assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?

  double const ret1 = functionTester(J, start, 6);
  std::cout << std::endl;

  start[0] = 279;
  start[1] = -96;

  double const ret2 = functionTester(J, start, 16);
  assert(std::abs(ret1 - ret2) < 1e-5);
}

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