Commit 9ef511b2 authored by Carsten Gräser's avatar Carsten Gräser
Browse files

Remove old bisection implementation

This was only available using a compile time switch.
Since the old style TNNMG implementation was ajusted
to the new bisection interface, this is no longer needed.
parent d04b5f3f
Pipeline #32281 passed with stage
in 5 minutes and 8 seconds
......@@ -18,20 +18,13 @@ class Bisection
* \param safety acceptance factor for inexact minimization
*/
Bisection(double acceptError = 0.0, double acceptFactor = 1.0, double requiredResidual = 1e-12,
#ifdef USE_OLD_TNNMG
bool fastQuadratic = true,
#endif
double safety = 1e-14) :
#ifdef USE_OLD_TNNMG
fastQuadratic_(fastQuadratic),
#endif
acceptFactor_(acceptFactor),
acceptError_(acceptError),
requiredResidual_(requiredResidual),
safety_(safety)
{};
#ifndef USE_OLD_TNNMG
/** \brief minimize convex functional
*
* Computes an approximate minimizer of the convex functional
......@@ -211,218 +204,6 @@ class Bisection
return x;
}
#else // The following is for the old tnnmg implementation
/** \brief minimize convex functional
*
* Computes an approximate minimizer of the convex functional
* @f[ f(x) = \frac 1 2 ax^2 - rx + \phi(x) @f]
* using bisection. If \f$ x^* \f$ is the exact minimizer
* and \f$ x^{old} \f$ is the old value the returned
* approximation @f$ x @f$ satisfies
* \f[
* (x - x^{old}) \in [ factor , 1] (x^*- x^{old}).
* \f]
* This guarantees enough descent for every fixed \f$ factor \f$ and hence convergence of the Gauss-Seidel method.
*
* \param J The convex functional
* \param x initial value
* \param x_old old value, needed for inexact minimization
* \param [out] count return the number of times the subdifferential of J was evaluated
* \return approximate minimizer
*/
template <class Functional>
double minimize(const Functional& J, double x, double x_old, int& count, int verbosity=0) const
{
size_t n=0;
Dune::Solvers::Interval<double> I;
Dune::Solvers::Interval<double> DJ;
Dune::Solvers::Interval<double> dom;
count = 0;
J.domain(dom);
double A = J.quadraticPart();
double b = J.linearPart();
// try to minimize the quadratic part exactly
if (fastQuadratic_)
{
double z = dom.projectIn(b/A);
J.subDiff(z, DJ);
++count;
if (DJ.containsZero(safety_))
return z;
}
// try if projected initial guess is sufficient
x = dom.projectIn(x);
J.subDiff(x, DJ);
++count;
if (DJ.containsZero(safety_))
{
if (verbosity > 0)
std::cout << "Bisection: initial iterate " << x << " accepted, DJ = " << DJ << std::endl;
return x;
}
// compute initial interval
// if quadratic part is strictly positive we can compute one bound from the other
if (DJ[0] > 0.0)
{
I[1] = x;
if (A > 0.0)
I[0] = dom.projectFromBelow(x - DJ[0]/A);
else
{
I[0] = dom.projectFromBelow(I[1] - DJ[0]);
J.subDiff(I[0], DJ);
++count;
while (DJ[0] > safety_)
{
I[0] = I[1] - 2.0*(I[1]-I[0]);
if (I[0] < dom[0])
{
I[0] = dom[0];
break;
}
J.subDiff(I[0], DJ);
++count;
}
}
}
else
{
I[0] = x;
if (A > 0.0)
I[1] = dom.projectFromAbove(x - DJ[1]/A);
else
{
I[1] = dom.projectFromAbove(I[0] - DJ[1]);
J.subDiff(I[1], DJ);
++count;
while (DJ[1] < -safety_)
{
I[1] = I[0] - 2.0*(I[0]-I[1]);
if (I[1] > dom[1])
{
I[1] = dom[1];
break;
}
J.subDiff(I[1], DJ);
++count;
}
}
}
if (verbosity>0)
{
std::cout << " # | I[0] | x | I[1] | dJ(I[0]) | dJ(x) | dJ(I[1]) |" << std::endl;
std::cout << "----|--------------|--------------|--------------|--------------|--------------|--------------|" << std::endl;
}
// compute midpoint
x = (I[0] + I[1]) / 2.0;
// We remember the value of x from the last loop in order
// to check if the iteration becomes stationary. This
// is necessary to guarantee termination because it
// might happen, that I[0] and I[1] are successive
// floating point numbers, while the values at both
// interval boundaries do not match the normal
// termination criterion.
//
// Checking ((I[0] < x) and (x < I[1])) does not work
// since this is true in the above case if x is taken
// from the register of a floating point unit that
// uses a higher accuracy internally
//
// The above does e.g. happen with the x87 floating
// point unit that uses 80 bit internally. In this
// specific case the problem could be avoided by
// changing the x87 rounding mode with
//
// int mode = 0x27F;
// asm ("fldcw %0" : : "m" (*&mode));
//
// In the general case one can avoid the problem at least
// for the equality check using a binary comparison
// that circumvents the floating point unit.
double xLast = I[0];
// iterate while mid point does not coincide with boundary numerically
while (not(binaryEqual(x, xLast)))
{
++n;
// evaluate subdifferential
J.subDiff(x, DJ);
++count;
if (verbosity>0)
{
Dune::Solvers::Interval<double> DJ0, DJ1;
J.subDiff(I[0],DJ0);
J.subDiff(I[1],DJ1);
const std::streamsize p = std::cout.precision();
std::cout << std::setw(4) << n << "|"
<< std::setprecision(8)
<< std::setw(14) << I[0] << "|"
<< std::setw(14) << x << "|"
<< std::setw(14) << I[1] << "|"
<< std::setw(14) << DJ0[0] << "|"
<< std::setw(14) << DJ[0] << "|"
<< std::setw(14) << DJ1[0] << "|" << std::endl;
std::cout << std::setprecision(p);
}
// stop if at least 'factor' of full minimization step is realized
// if DJ[0]<= 0 we have
//
// x_old <= I[0] <= x <= x* <= I[1]
//
// and hence the criterion
//
// alpha|x*-x_old| <= alpha|I[1]-x_old| <= |x-x_old| <= |x*-x_old|
//
// similar for DJ[1] >= 0
if ((I[0]>=x_old) and (DJ[0]<=0.0) and (std::abs((x-x_old)/(I[1]-x_old)) >= acceptFactor_) and (-requiredResidual_ <= DJ[1]))
break;
if ((I[1]<=x_old) and (DJ[1]>=0.0) and (std::abs((x-x_old)/(I[0]-x_old)) >= acceptFactor_) and (DJ[0] <= requiredResidual_))
break;
// stop if error is less than acceptError_
if (std::abs(I[1]-I[0]) <= acceptError_)
break;
// stop if 0 in DJ numerically
if (DJ.containsZero(safety_))
break;
// adjust bounds and compute new mid point
if (DJ[0] > 0.0)
I[1] = x;
else
I[0] = x;
xLast = x;
x = (I[0] + I[1]) / 2.0;
}
return x;
}
void
setFastQuadratic(bool fastQuadratic) {
fastQuadratic_ = fastQuadratic;
}
#endif
private:
......@@ -487,9 +268,6 @@ class Bisection
return true;
}
#ifdef USE_OLD_TNNMG
bool fastQuadratic_;
#endif
double acceptFactor_;
double acceptError_;
double requiredResidual_;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment