diff --git a/src/Makefile.am b/src/Makefile.am
index 39fc4c205f221a5f0028753592680c160e5a65f0..b1eb1ad10590ef42da8c0fb37c9b672f4c2c100a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -2,25 +2,13 @@
SUBDIRS =
noinst_PROGRAMS = \
- bisection-example \
- bisection-example-flexible \
- bisection-simpler-example \
- bisection-simpler-example2 \
- bisection-simpler-example2-gradient
+ bisection-example-flexible
-bisection_example_SOURCES = bisection-example.cc
bisection_example_flexible_SOURCES = \
bisection-example-flexible.cc \
properscalarincreasingconvexfunction.hh
-bisection_simpler_example_SOURCES = bisection-simpler-example.cc
-bisection_simpler_example2_SOURCES = bisection-simpler-example2.cc
-bisection_simpler_example2_gradient_SOURCES = bisection-simpler-example2-gradient.cc
check-am:
- ./bisection-simpler-example
- ./bisection-simpler-example2
- ./bisection-simpler-example2-gradient
- ./bisection-example
./bisection-example-flexible
AM_CXXFLAGS = -Wall -Wextra
diff --git a/src/bisection-example.cc b/src/bisection-example.cc
deleted file mode 100644
index 37c345bacd80d821b4b1ead13e2bc62e87b59db9..0000000000000000000000000000000000000000
--- a/src/bisection-example.cc
+++ /dev/null
@@ -1,207 +0,0 @@
-/* -*- 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>
-
-template <int dimension> class SampleFunctional {
-public:
- typedef Dune::FieldVector<double, dimension> SmallVector;
- typedef Dune::FieldMatrix<double, dimension, dimension> SmallMatrix;
-
- // FIXME: hardcoded function H
- 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 + H(v.two_norm()); // <1/2 Av - b,v> + H(|v|)
- }
-
- // |dir|-times the directional derivative wrt dir/|dir|.
- double pseudoDirectionalDerivative(const SmallVector x,
- const SmallVector dir) const {
- if (x == SmallVector(0.0))
- return HPrimePlus(0) * dir.two_norm();
-
- if (x * dir > 0)
- return PlusGrad(x) * dir;
- else
- return MinusGrad(x) * dir;
- }
-
- 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);
-
- 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 (directionalDerivative(tmp, descDir) >= 0)
- break;
-
- l = r;
- r *= 2;
- Dune::dverb << "Widened interval!" << std::endl;
- }
- Dune::dverb << "Interval now [" << l << "," << r << "]" << std::endl;
-
- // Debugging
- {
- SmallVector tmpl = x;
- tmpl.axpy(l, descDir);
- SmallVector tmpr = x;
- tmpr.axpy(r, descDir);
- assert(directionalDerivative(tmpl, descDir) < 0);
- assert(directionalDerivative(tmpr, descDir) > 0);
- }
-
- double m = l / 2 + r / 2;
- SmallVector middle;
- 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 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_;
-
- double H(double s) const { return (s < 1) ? s : (2 * s - 1); }
-
- double HPrimeMinus(double s) const { return (s <= 1) ? 1 : 2; }
-
- double HPrimePlus(double s) const { return (s < 1) ? 1 : 2; }
-
- // 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(HPrimePlus(x.two_norm()) / x.two_norm(), x);
- return y;
- }
-
- SmallVector MinusGrad(const SmallVector x) const {
- SmallVector y = SmoothGrad(x);
- y.axpy(HPrimeMinus(x.two_norm()) / x.two_norm(), x);
- return y;
- }
-
- SmallVector ModifiedGradient(const SmallVector x) const {
- if (x == SmallVector(0.0))
- // TODO
- DUNE_THROW(Dune::Exception, "The case x = 0 is not yet handled.");
-
- 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;
- }
-};
-
-int main() {
- try {
- 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.pseudoDirectionalDerivative(b, b) == 10 + 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.254644) < 1e-8);
- std::cout << J(start + correction) << std::endl;
-
- return 0;
- }
- catch (Dune::Exception &e) {
- Dune::derr << "Dune reported error: " << e << std::endl;
- }
-}
diff --git a/src/bisection-simpler-example.cc b/src/bisection-simpler-example.cc
deleted file mode 100644
index 5c12666376e7152554706d6d3c24e5639b576c71..0000000000000000000000000000000000000000
--- a/src/bisection-simpler-example.cc
+++ /dev/null
@@ -1,58 +0,0 @@
-/* -*- mode:c++; mode: flymake -*- */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/fufem/interval.hh>
-
-#include <dune/tnnmg/nonlinearities/smallfunctional.hh>
-#include <dune/tnnmg/problem-classes/bisection.hh>
-
-#include <cassert>
-#include <limits>
-
-class SampleFunctional : public SmallFunctional<2> {
-public:
- 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) const { return A_; }
-
-private:
- SmallMatrix A_;
- SmallVector b_;
-};
-
-int main() {
- 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);
-
- assert(J(b) == 2.5);
-
- assert(b.two_norm() == sqrt(5));
- assert(b.two_norm2() == 5);
- return 0;
-}
diff --git a/src/bisection-simpler-example2-gradient.cc b/src/bisection-simpler-example2-gradient.cc
deleted file mode 100644
index 90d09a43e6b9e3045445192ea1fcbdcbad0c23d1..0000000000000000000000000000000000000000
--- a/src/bisection-simpler-example2-gradient.cc
+++ /dev/null
@@ -1,135 +0,0 @@
-/* -*- mode:c++; mode: flymake -*- */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#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 <iostream>
-
-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) 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, unsigned 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;
- while (true) {
- tmp = x;
- tmp.axpy(r, descDir);
- if (directionalDerivative(tmp, descDir) >= 0)
- break;
-
- l = r;
- r *= 2;
- Dune::dverb << "Widened interval!" << std::endl;
- }
- Dune::dverb << "Interval now [" << l << "," << r << "]" << std::endl;
-
- // Debugging
- {
- SmallVector tmpl = x;
- tmpl.axpy(l, descDir);
- SmallVector tmpr = x;
- tmpr.axpy(r, descDir);
- assert(directionalDerivative(tmpl, descDir) < 0);
- assert(directionalDerivative(tmpr, descDir) > 0);
- }
-
- double m = l / 2 + r / 2;
- SmallVector middle;
- 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 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));
- assert(std::abs(J(start + correction) + 0.833333) < 1e-6);
- std::cout << J(start + correction) << std::endl;
-
- return 0;
-}
diff --git a/src/bisection-simpler-example2.cc b/src/bisection-simpler-example2.cc
deleted file mode 100644
index a34baf5d2635663990da3aba88cfcc847cc23106..0000000000000000000000000000000000000000
--- a/src/bisection-simpler-example2.cc
+++ /dev/null
@@ -1,64 +0,0 @@
-/* -*- mode:c++; mode: flymake -*- */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#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) const { return A_; }
-
-private:
- SmallMatrix A_;
- SmallVector b_;
-};
-
-int main() {
- int const dim = 2;
-
- SampleFunctional<dim>::SmallMatrix A;
- A[0][0] = 3;
- A[0][1] = 0;
- A[1][0] = 0;
- A[1][1] = 3;
- SampleFunctional<dim>::SmallVector b;
- b[0] = 1;
- b[1] = 2;
-
- SampleFunctional<dim> J(A, b);
-
- assert(J(b) == 2.5);
-
- assert(b.two_norm() == sqrt(5));
- assert(b.two_norm2() == 5);
- return 0;
-}