test-gradient-method.cc

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

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

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

#include "samplefunctional.hh"

#include <cassert>

template <int dim>
double functionTester(Dune::SampleFunctional<dim> J,
                      typename Dune::SampleFunctional<dim>::SmallVector &start,
                      size_t runs) {
  typename Dune::SampleFunctional<dim>::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 testIdentity() {
  int const dim = 2;
  typedef Dune::SampleFunctional<dim> Functional;

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

  Dune::Identity f;
  Functional J(A, b, Dune::MyNonlinearity<dim>(f));

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

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

  // Something random
  start[0] = 279;
  start[1] = -96;

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

  start[0] = 0;
  start[1] = 0;

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

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

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

  Dune::SampleFunction f;
  Functional J(A, b, Dune::MyNonlinearity<dim>(f));

  Functional::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))
  */
  Functional::SmallVector error;
  error[0] = -(102 - 1 + 2 / sqrt(5));
  error[1] = -(161.5 - 2 + 4 / sqrt(5));
  Functional::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);

  // Something random
  start[0] = 279;
  start[1] = -96;

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

  start[0] = 0;
  start[1] = 0;

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

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

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

  Dune::SampleFunction f;
  Functional J(A, b, Dune::MyNonlinearity<dim>(f));

  Functional::SmallVector start;
  Functional::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);

    Functional::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);

    Functional::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> Functional;

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

  Dune::TrivialFunction f;
  Functional J(A, b, Dune::MyNonlinearity<dim>(f));

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

  /*
    j(x)
    = Ax - b
    = 17*(6, 9.5) - (1, 2)
    = (102 - 1, 161.5 - 2)
  */
  Functional::SmallVector error;
  error[0] = -101;
  error[1] = -159.5;
  Functional::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;

  // Something random
  start[0] = 279;
  start[1] = -96;

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

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

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

  Dune::HorribleFunction f;
  Functional J(A, b, Dune::MyNonlinearity<dim>(f));

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

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

  // Something random
  start[0] = 279;
  start[1] = -96;

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

  start[0] = 0;
  start[1] = 0;

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

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

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

  Dune::HorribleFunctionLogarithmic f;
  Functional J(A, b, Dune::MyNonlinearity<dim>(f));

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

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

  // Something random
  start[0] = 279;
  start[1] = -96;

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

  start[0] = 0;
  start[1] = 0;

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

void testSampleFunction3D() {
  int const dim = 3;
  typedef Dune::SampleFunctional<dim> Functional;

  Functional::SmallMatrix A(0);
  A[0][0] = 3;
  A[0][1] = 1.5;
  A[1][0] = 1.5;
  A[1][1] = 4;
  A[1][2] = 1.5;
  A[2][1] = 1.5;
  A[2][2] = 5;
  Functional::SmallVector b;
  b[0] = 1;
  b[1] = 2;
  b[2] = 3;

  Dune::SampleFunction f;
  Functional J(A, b, Dune::MyNonlinearity<dim>(f));

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

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

  // Something random
  start[0] = 279;
  start[1] = -96;
  start[2] = -15;

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

  start[0] = 0;
  start[1] = 0;
  start[2] = 0;

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

// Checks if reaching the minimum in one step leads to problems
void testSampleFunction2() {
  int const dim = 2;
  typedef Dune::SampleFunctional<dim> Functional;

  Functional::SmallMatrix A;
  A[0][0] = 1;
  A[0][1] = 0;
  A[1][0] = 0;
  A[1][1] = 1;
  Functional::SmallVector b;
  b[0] = 1;
  b[1] = 1;

  Dune::SampleFunction f;
  Functional J(A, b, Dune::MyNonlinearity<dim>(f));

  Functional::SmallVector start = b;

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

  // Something random
  start[0] = 279;
  start[1] = -96;

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

  start[0] = 0;
  start[1] = 0;

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

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