Skip to content

Commit

Permalink
test only running model parameters for convergence
Browse files Browse the repository at this point in the history 
The GSL root solvers sometimes don't change the parameters from one
iteration to the other.  This might lead to an early stop of the
algorithm, which has not yet converged to a good root.
  • Loading branch information
Alexander Voigt committed Jul 13, 2015
1 parent 812de0d commit fd2f63f
Show file tree
Hide file tree
Showing 3 changed files with 17 additions and 114 deletions.
97 changes: 14 additions & 83 deletions src/root_convergence_tester_drbar.hpp
View file
Expand Up @@ -21,6 +21,7 @@

#include "root_convergence_tester.hpp"
#include "convergence_tester_drbar.hpp"
#include "functors.hpp"
#include "logger.hpp"
#include "numerics2.hpp"

Expand All @@ -36,7 +37,7 @@ template <template<class Method> class Model>
class Convergence_tester_DRbar<Model<Root> > :
public Convergence_tester<Root> {
public:
Convergence_tester_DRbar(Model<Root>*, double);
Convergence_tester_DRbar(double);
virtual ~Convergence_tester_DRbar();

virtual bool accuracy_goal_reached(const Eigen::ArrayXd&);
Expand All @@ -46,36 +47,22 @@ class Convergence_tester_DRbar<Model<Root> > :
void set_max_iterations(std::size_t); ///< set maximum number of iterations

protected:
std::size_t get_iteration() const; ///< get current iteration number
const Model<Root>& get_model() const; ///< get model
const Model<Root>& get_last_iteration_model() const; ///< get model state during last iteration
virtual double max_rel_diff() const = 0; ///< maximum relative difference to last iteration
virtual double rel_scale_difference() const; ///< relative scale difference
virtual double scale_difference() const; ///< absolute scale difference
virtual bool scale_has_changed() const; ///< returns true if scale has changed
virtual double max_rel_diff() const; ///< maximum relative difference to last iteration

private:
Model<Root>* model; ///< pointer to model
Model<Root> last_iteration_model; ///< model state at last iteration
std::size_t it_count; ///< iteration
std::size_t max_it; ///< maximum number of iterations
double accuracy_goal; ///< accuracy goal
double current_accuracy; ///< current accuracy
};

template <template<class Method> class Model>
Convergence_tester_DRbar<Model<Root> >::Convergence_tester_DRbar
(Model<Root>* model_, double accuracy_goal_)
(double accuracy_goal_)
: Convergence_tester<Root>()
, model(model_)
, last_iteration_model()
, it_count(0)
, max_it(static_cast<int>(-log(accuracy_goal_) / log(10.0) * 10))
, accuracy_goal(accuracy_goal_)
, current_accuracy(std::numeric_limits<double>::infinity())
{
assert(model && "Error: Convergence_tester_DRbar<Model<Root>>: "
"model pointer must not be zero!");
}

template <template<class Method> class Model>
Expand All @@ -84,34 +71,15 @@ Convergence_tester_DRbar<Model<Root> >::~Convergence_tester_DRbar()
}

template <template<class Method> class Model>
bool Convergence_tester_DRbar<Model<Root> >::accuracy_goal_reached(const Eigen::ArrayXd&)
bool Convergence_tester_DRbar<Model<Root> >::accuracy_goal_reached(const Eigen::ArrayXd& x)
{
bool precision_reached;
if (it_count == 0) {
// this is the first run => no comparison possible => assume
// that accuracy goal has not been reached
precision_reached = false;
} else {
const double scale_accuracy_goal = accuracy_goal * 16*M_PI*M_PI;
if (rel_scale_difference() < scale_accuracy_goal) {
current_accuracy = max_rel_diff();
precision_reached = current_accuracy < accuracy_goal;
VERBOSE_MSG("Convergence_tester_DRbar: current accuracy = "
<< current_accuracy
<< ", accuracy goal = " << accuracy_goal);
} else {
precision_reached = false;
VERBOSE_MSG("scale has changed by " << scale_difference()
<< " GeV (" << rel_scale_difference()
<< "), skipping parameter comparison");
}
}

// save old model parameters
last_iteration_model = *model;
++it_count;

return precision_reached;
current_accuracy = x.cwiseAbs().unaryExpr(Chop<double>(1e-10)).maxCoeff();

VERBOSE_MSG("Convergence_tester_DRbar: current accuracy = "
<< current_accuracy
<< ", accuracy goal = " << accuracy_goal);

return current_accuracy < accuracy_goal;
}

template <template<class Method> class Model>
Expand All @@ -127,23 +95,9 @@ double Convergence_tester_DRbar<Model<Root> >::get_current_accuracy() const
}

template <template<class Method> class Model>
std::size_t Convergence_tester_DRbar<Model<Root> >::get_iteration() const
double Convergence_tester_DRbar<Model<Root> >::max_rel_diff() const
{
return it_count;
}

template <template<class Method> class Model>
const Model<Root>&
Convergence_tester_DRbar<Model<Root> >::get_model() const
{
return *model;
}

template <template<class Method> class Model>
const Model<Root>&
Convergence_tester_DRbar<Model<Root> >::get_last_iteration_model() const
{
return last_iteration_model;
return std::numeric_limits<double>::infinity();
}

template <template<class Method> class Model>
Expand All @@ -160,29 +114,6 @@ std::size_t Convergence_tester_DRbar<Model<Root> >::max_iterations()
return max_it;
}

template <template<class Method> class Model>
bool Convergence_tester_DRbar<Model<Root> >::scale_has_changed() const
{
return !is_zero(scale_difference());
}

template <template<class Method> class Model>
double Convergence_tester_DRbar<Model<Root> >::scale_difference() const
{
return model->get_scale() - last_iteration_model.get_scale();
}

template <template<class Method> class Model>
double Convergence_tester_DRbar<Model<Root> >::rel_scale_difference()
const
{
const double diff = scale_difference();
const double last_scale = last_iteration_model.get_scale();
if (!is_zero(last_scale))
return diff / last_scale;
return std::numeric_limits<double>::infinity();
}

} // namespace flexiblesusy

#endif
31 changes: 3 additions & 28 deletions templates/root_convergence_tester.cpp.in
View file
Expand Up @@ -19,42 +19,17 @@
// File generated at @DateAndTime@

#include "@ModelName@_root_convergence_tester.hpp"
#include <cmath>
#include <algorithm>
#include "wrappers.hpp"

namespace flexiblesusy {

#define OLD(p) ol.get_##p()
#define NEW(p) ne.get_##p()

#define OLD1(p,i) ol.get_##p()(i)
#define NEW1(p,i) ne.get_##p()(i)

#define OLD2(p,i,j) ol.get_##p(i,j)
#define NEW2(p,i,j) ne.get_##p(i,j)

#define OLD3(p,i,j,k) ol.get_##p(i,j,k)
#define NEW3(p,i,j,k) ne.get_##p(i,j,k)

#define OLD4(p,i,j,k,l) ol.get_##p(i,j,k,l)
#define NEW4(p,i,j,k,l) ne.get_##p(i,j,k,l)

@ModelName@_convergence_tester<Root>::@ModelName@_convergence_tester(@ModelName@<Root>* model, double accuracy_goal)
: Convergence_tester_DRbar<@ModelName@<Root> >(model, accuracy_goal)
@ModelName@_convergence_tester<Root>::@ModelName@_convergence_tester(
@ModelName@<Root>*, double accuracy_goal)
: Convergence_tester_DRbar<@ModelName@<Root> >(accuracy_goal)
{
}

@ModelName@_convergence_tester<Root>::~@ModelName@_convergence_tester()
{
}

double @ModelName@_convergence_tester<Root>::max_rel_diff() const
{
const @ModelName@<Root>& ol = get_last_iteration_model();
const @ModelName@<Root>& ne = get_model();

@compareFunction@
}

} // namespace flexiblesusy
3 changes: 0 additions & 3 deletions templates/root_convergence_tester.hpp.in
View file
Expand Up @@ -34,9 +34,6 @@ class @ModelName@_convergence_tester<Root> : public Convergence_tester_DRbar<@Mo
public:
@ModelName@_convergence_tester(@ModelName@<Root>*, double);
virtual ~@ModelName@_convergence_tester();

protected:
virtual double max_rel_diff() const;
};

} // namespace flexiblesusy
Expand Down

0 comments on commit fd2f63f

Please sign in to comment.