Skip to content
Snippets Groups Projects
Commit e56be5fb authored by Elias Pipping's avatar Elias Pipping Committed by Elias Pipping
Browse files

Nuke obsolete files

parent a1823d67
No related branches found
No related tags found
No related merge requests found
......@@ -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
......
/* -*- 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;
}
}
/* -*- 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;
}
/* -*- 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;
}
/* -*- 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;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment