From d9ddc9013af89924c4d575fb007b0e380b62def6 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 6 Aug 2012 11:21:29 +0200
Subject: [PATCH] Move minimisation code to a separate header
---
dune/tectonic/ellipticenergy.hh | 170 ---------------------------
dune/tectonic/minimisation.hh | 182 +++++++++++++++++++++++++++++
dune/tectonic/myblockproblem.hh | 1 +
src/test-gradient-method-helper.hh | 1 +
4 files changed, 184 insertions(+), 170 deletions(-)
create mode 100644 dune/tectonic/minimisation.hh
diff --git a/dune/tectonic/ellipticenergy.hh b/dune/tectonic/ellipticenergy.hh
index c30626f7..274d8c4c 100644
--- a/dune/tectonic/ellipticenergy.hh
+++ b/dune/tectonic/ellipticenergy.hh
@@ -6,11 +6,8 @@
#include <dune/common/stdstreams.hh>
#include <dune/fufem/interval.hh>
-#include <dune/tnnmg/problem-classes/bisection.hh>
-#include "mydirectionalconvexfunction.hh"
#include "localnonlinearity.hh"
-#include "circularconvexfunction.hh"
namespace Dune {
template <int dim> class EllipticEnergy {
@@ -146,172 +143,5 @@ template <int dim> class EllipticEnergy {
y.axpy(-dx / xx, x);
}
};
-
-template <class Functional>
-void descentMinimisation(Functional const &J,
- typename Functional::SmallVector &x,
- typename Functional::SmallVector const &descDir,
- Bisection const &bisection) {
- typedef typename Functional::SmallVector SmallVector;
- typedef typename Functional::NonlinearityType LocalNonlinearityType;
-
- // {{{ Construct a restriction of J to the line x + t * descDir
-
- /* We have
-
- 1/2 <A(u+xv),u+xv>-<b,u+xv> = 1/2 <Av,v> x^2 - <b-Au,v> x + <1/2 Au-b,u>
-
- since A is symmetric.
- */
- SmallVector tmp = J.b; // b
- J.A.mmv(x, tmp); // b-Au
- double const JRestb = tmp * descDir; // <b-Au,v>
-
- J.A.mv(descDir, tmp); // Av
- double const JRestA = tmp * descDir; // <Av,v>
-
- MyDirectionalConvexFunction<LocalNonlinearityType> const JRest(
- JRestA, JRestb, *J.phi, x, descDir);
- // }}}
-
- { // Debug
- Interval<double> D;
- JRest.subDiff(0, D);
-
- dverb
- << "## Directional derivative (as per subdifferential of restriction): "
- << D[1] << " (coordinates of the restriction)" << std::endl;
- /*
- It is possible that this differs quite a lot from the
- directional derivative computed in the descentDirection()
- method:
-
- If phi is nonsmooth at x, so that the directional
- derivatives jump, and |x| is computed to be too small or too
- large globally or locally, the locally computed
- subdifferential and the globally computed subdifferential
- will no longer coincide!
-
- The assertion D[1] <= 0 may thus fail.
- */
- }
-
- int count;
- double const stepsize = bisection.minimize(JRest, 0.0, 0.0, count);
- dverb << "Number of iterations in the bisection method: " << count
- << std::endl;
- ;
-
- x.axpy(stepsize, descDir);
-}
-
-template <class Functional>
-void tangentialMinimisation(Functional const &J,
- typename Functional::SmallVector &x,
- typename Functional::SmallVector const &descDir,
- Bisection const &bisection) {
- typedef typename Functional::SmallMatrix SmallMatrix;
- typedef typename Functional::SmallVector SmallVector;
-
- // We completely ignore the nonlinearity here -- when restricted
- // to a circle, it just enters as a constant!
- CircularConvexFunction<SmallMatrix, SmallVector> const JRest(J.A, J.b, x,
- descDir);
-
- int count;
- double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
- dverb << "Number of iterations in the bisection method: " << count
- << std::endl;
- ;
-
- // Since x is used in the computation of the rhs, do not write to it directly
- SmallVector tmp;
- JRest.cartesian(stepsize, tmp);
- x = tmp;
-}
-
-template <class Functional>
-void minimise(Functional const &J, typename Functional::SmallVector &x,
- size_t steps, Bisection const &bisection) {
- typedef typename Functional::SmallVector SmallVector;
-
- for (size_t step = 0; step < steps; ++step) {
- SmallVector descDir;
- bool const linesearchp = J.descentDirection(x, descDir);
-
- if (descDir == SmallVector(0.0))
- return;
-
- if (linesearchp) {
- descentMinimisation(J, x, descDir, bisection);
- } else {
- Bisection slowBisection(bisection);
- slowBisection.setFastQuadratic(false);
-
- tangentialMinimisation(J, x, descDir, slowBisection);
- }
- }
-}
-
-// NOTE: only works for the 2D case
-template <class Functional>
-void minimise2(Functional const &J, typename Functional::SmallVector &x,
- size_t steps, Bisection const &bisection) {
- typedef typename Functional::SmallVector SmallVector;
- minimisationInitialiser(J, bisection, x);
-
- { // Make a single step where we already know we're not differentiable
- SmallVector descDir;
- if (x == SmallVector(0))
- J.descentAtZero(descDir);
- else
- J.descentDirectionNew(x, descDir);
- descentMinimisation(J, x, descDir, bisection);
- }
-
- // From here on, we're in a C1 region and stay there.
- for (size_t i = 1; i < steps; ++i) {
- SmallVector descDir;
- J.descentDirectionNew(x, descDir);
- descentMinimisation(J, x, descDir, bisection);
- }
-}
-
-template <class Functional>
-void minimisationInitialiser(Functional const &J, Bisection const &bisection,
- typename Functional::SmallVector &startingPoint) {
- typedef typename Functional::SmallVector SmallVector;
-
- std::vector<double> const kinks = { 5, 10,
- 15 }; // FIXME: We happen to know these
-
- SmallVector x_old(0);
- double Jx_old = J(x_old);
-
- for (auto &kink : kinks) {
- SmallVector x1 = { kink, 0 }; // Random vector that has the right norm
- SmallVector const descDir1 = { x1[1], -x1[0] };
- tangentialMinimisation(J, x1, descDir1, bisection);
- double const Jx1 = J(x1);
-
- SmallVector x2(0);
- x2.axpy(-1, x1);
- SmallVector const descDir2 = { x2[1], -x2[0] };
- tangentialMinimisation(J, x2, descDir2, bisection);
- double const Jx2 = J(x2);
-
- double const Jx = (Jx2 < Jx1 ? Jx2 : Jx1);
- SmallVector const x = (Jx2 < Jx1 ? x2 : x1);
-
- if (Jx >= Jx_old) {
- startingPoint = x_old;
- return;
- }
-
- Jx_old = Jx;
- x_old = x;
- }
- startingPoint = x_old;
-}
}
#endif
diff --git a/dune/tectonic/minimisation.hh b/dune/tectonic/minimisation.hh
new file mode 100644
index 00000000..113c101a
--- /dev/null
+++ b/dune/tectonic/minimisation.hh
@@ -0,0 +1,182 @@
+#ifndef MINIMISATION_HH
+#define MINIMISATION_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/stdstreams.hh>
+
+#include <dune/fufem/interval.hh>
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include "circularconvexfunction.hh"
+#include "mydirectionalconvexfunction.hh"
+
+namespace Dune {
+template <class Functional>
+void descentMinimisation(Functional const &J,
+ typename Functional::SmallVector &x,
+ typename Functional::SmallVector const &descDir,
+ Bisection const &bisection) {
+ typedef typename Functional::SmallVector SmallVector;
+ typedef typename Functional::NonlinearityType LocalNonlinearityType;
+
+ // {{{ Construct a restriction of J to the line x + t * descDir
+
+ /* We have
+
+ 1/2 <A(u+xv),u+xv>-<b,u+xv> = 1/2 <Av,v> x^2 - <b-Au,v> x + <1/2 Au-b,u>
+
+ since A is symmetric.
+ */
+ SmallVector tmp = J.b; // b
+ J.A.mmv(x, tmp); // b-Au
+ double const JRestb = tmp * descDir; // <b-Au,v>
+
+ J.A.mv(descDir, tmp); // Av
+ double const JRestA = tmp * descDir; // <Av,v>
+
+ MyDirectionalConvexFunction<LocalNonlinearityType> const JRest(
+ JRestA, JRestb, *J.phi, x, descDir);
+ // }}}
+
+ { // Debug
+ Interval<double> D;
+ JRest.subDiff(0, D);
+
+ dverb
+ << "## Directional derivative (as per subdifferential of restriction): "
+ << D[1] << " (coordinates of the restriction)" << std::endl;
+ /*
+ It is possible that this differs quite a lot from the
+ directional derivative computed in the descentDirection()
+ method:
+
+ If phi is nonsmooth at x, so that the directional
+ derivatives jump, and |x| is computed to be too small or too
+ large globally or locally, the locally computed
+ subdifferential and the globally computed subdifferential
+ will no longer coincide!
+
+ The assertion D[1] <= 0 may thus fail.
+ */
+ }
+
+ int count;
+ double const stepsize = bisection.minimize(JRest, 0.0, 0.0, count);
+ dverb << "Number of iterations in the bisection method: " << count
+ << std::endl;
+ ;
+
+ x.axpy(stepsize, descDir);
+}
+
+template <class Functional>
+void tangentialMinimisation(Functional const &J,
+ typename Functional::SmallVector &x,
+ typename Functional::SmallVector const &descDir,
+ Bisection const &bisection) {
+ typedef typename Functional::SmallMatrix SmallMatrix;
+ typedef typename Functional::SmallVector SmallVector;
+
+ // We completely ignore the nonlinearity here -- when restricted
+ // to a circle, it just enters as a constant!
+ CircularConvexFunction<SmallMatrix, SmallVector> const JRest(J.A, J.b, x,
+ descDir);
+
+ int count;
+ double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
+ dverb << "Number of iterations in the bisection method: " << count
+ << std::endl;
+ ;
+
+ // Since x is used in the computation of the rhs, do not write to it directly
+ SmallVector tmp;
+ JRest.cartesian(stepsize, tmp);
+ x = tmp;
+}
+
+template <class Functional>
+void minimise(Functional const &J, typename Functional::SmallVector &x,
+ size_t steps, Bisection const &bisection) {
+ typedef typename Functional::SmallVector SmallVector;
+
+ for (size_t step = 0; step < steps; ++step) {
+ SmallVector descDir;
+ bool const linesearchp = J.descentDirection(x, descDir);
+
+ if (descDir == SmallVector(0.0))
+ return;
+
+ if (linesearchp) {
+ descentMinimisation(J, x, descDir, bisection);
+ } else {
+ Bisection slowBisection(bisection);
+ slowBisection.setFastQuadratic(false);
+
+ tangentialMinimisation(J, x, descDir, slowBisection);
+ }
+ }
+}
+
+// NOTE: only works for the 2D case
+template <class Functional>
+void minimise2(Functional const &J, typename Functional::SmallVector &x,
+ size_t steps, Bisection const &bisection) {
+ typedef typename Functional::SmallVector SmallVector;
+ minimisationInitialiser(J, bisection, x);
+
+ { // Make a single step where we already know we're not differentiable
+ SmallVector descDir;
+ if (x == SmallVector(0))
+ J.descentAtZero(descDir);
+ else
+ J.descentDirectionNew(x, descDir);
+ descentMinimisation(J, x, descDir, bisection);
+ }
+
+ // From here on, we're in a C1 region and stay there.
+ for (size_t i = 1; i < steps; ++i) {
+ SmallVector descDir;
+ J.descentDirectionNew(x, descDir);
+ descentMinimisation(J, x, descDir, bisection);
+ }
+}
+
+template <class Functional>
+void minimisationInitialiser(Functional const &J, Bisection const &bisection,
+ typename Functional::SmallVector &startingPoint) {
+ typedef typename Functional::SmallVector SmallVector;
+
+ std::vector<double> const kinks = { 5, 10,
+ 15 }; // FIXME: We happen to know these
+
+ SmallVector x_old(0);
+ double Jx_old = J(x_old);
+
+ for (auto &kink : kinks) {
+ SmallVector x1 = { kink, 0 }; // Random vector that has the right norm
+ SmallVector const descDir1 = { x1[1], -x1[0] };
+ tangentialMinimisation(J, x1, descDir1, bisection);
+ double const Jx1 = J(x1);
+
+ SmallVector x2(0);
+ x2.axpy(-1, x1);
+ SmallVector const descDir2 = { x2[1], -x2[0] };
+ tangentialMinimisation(J, x2, descDir2, bisection);
+ double const Jx2 = J(x2);
+
+ double const Jx = (Jx2 < Jx1 ? Jx2 : Jx1);
+ SmallVector const x = (Jx2 < Jx1 ? x2 : x1);
+
+ if (Jx >= Jx_old) {
+ startingPoint = x_old;
+ return;
+ }
+
+ Jx_old = Jx;
+ x_old = x;
+ }
+ startingPoint = x_old;
+}
+}
+#endif
diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index dc447c98..fb070099 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -12,6 +12,7 @@
#include "globalnonlinearity.hh"
#include "localnonlinearity.hh"
+#include "minimisation.hh"
#include "mydirectionalconvexfunction.hh"
#include "ellipticenergy.hh"
diff --git a/src/test-gradient-method-helper.hh b/src/test-gradient-method-helper.hh
index f3016d5a..a9479ebb 100644
--- a/src/test-gradient-method-helper.hh
+++ b/src/test-gradient-method-helper.hh
@@ -6,6 +6,7 @@
#include <dune/tnnmg/problem-classes/bisection.hh>
#include <dune/tectonic/ellipticenergy.hh>
+#include <dune/tectonic/minimisation.hh>
template <int dim>
double functionTester(Dune::EllipticEnergy<dim> J,
--
GitLab