root_solver.cpp

// ====================================================================
// This file is part of FlexibleSUSY.
//
// FlexibleSUSY is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// FlexibleSUSY is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FlexibleSUSY.  If not, see
// <http://www.gnu.org/licenses/>.
// ====================================================================

#include "root_solver.hpp"
#include "root_constraint.hpp"
#include "root_convergence_tester.hpp"
#include "root_initial_guesser.hpp"
#include "root_model.hpp"

#include "error.hpp"
#include "functors.hpp"
#include "logger.hpp"

#include <sstream>

#include <gsl/gsl_multiroots.h>

namespace flexiblesusy {

/**
 * Create empty two scale solver.
 * The RG running precision is set to the default value 0.001.
 */
RGFlow<Root>::RGFlow()
   : models()
   , iteration(0)
   , convergence_tester(NULL)
   , initial_guesser(NULL)
   , precision_goal(1e-4)
{
}

RGFlow<Root>::~RGFlow()
{
   delete_models();
}

/**
 * @brief Solves the boundary value problem.
 *
 * At first the initial_guess() is called.  Afterwards, the function
 * iteratively runs the tower up and down and imposes the boundary
 * conditions.  The iteration stops if either the maximum number of
 * iterations is reached or the precision goal is achieved (defined by
 * the convergence_tester).
 */
void RGFlow<Root>::solve()
{
   check_setup();

   const unsigned max_iterations = get_max_iterations();
   if (models.empty() || max_iterations == 0)
      return;

   initial_guess();

   const int status = find_root();

   if (status != GSL_SUCCESS)
      throw NoConvergenceError(max_iterations);

   VERBOSE_MSG("convergence reached after " << iteration << " iterations");
}

int RGFlow<Root>::func(const gsl_vector* x, void* params, gsl_vector* f)
{
   RGFlow<Root>* root_finder = static_cast<RGFlow<Root>*>(params);
   int status = GSL_SUCCESS;

   const std::size_t dimx = x->size;
   const std::size_t dimy = f->size;

   // vector for calculating function
   Eigen::ArrayXd xa(to_eigen_array(x));

   try {
      // calculate function
      const Eigen::ArrayXd func(root_finder->calculate_function(xa));

      to_gsl_vector(func, f);

   } catch (const Error& e) {

#ifdef ENABLE_VERBOSE
      WARNING("\tProblem in RGFlow<Root>::func: " << e.what());
#endif

      status = GSL_EDOM;
   }

   return status;
}

int RGFlow<Root>::find_root()
{
   const std::size_t dimension = get_dimension();

   int status;
   std::size_t iter = 0;
   gsl_multiroot_function f = {func, dimension, this};

   const gsl_multiroot_fsolver_type* solver_type
      = gsl_multiroot_fsolver_hybrid; // @todo make configurable

   gsl_multiroot_fsolver* solver
      = gsl_multiroot_fsolver_alloc(solver_type, dimension);

   if (!solver) {
      throw OutOfMemoryError(std::string("Cannot allocate gsl_multiroot_fsolver ") +
                             gsl_multiroot_fsolver_name(solver));
   }

#ifndef ENABLE_DEBUG
   gsl_set_error_handler_off();
#endif

   gsl_vector* x_init = gsl_vector_alloc(get_dimension());
   if (!x_init)
      throw OutOfMemoryError("GSL vector allocation failed in RGFlow<Root>()");

   starting_point(x_init);

   gsl_multiroot_fsolver_set(solver, &f, x_init);

#ifdef ENABLE_VERBOSE
   print_state(solver, iter); // @todo implement
#endif

   do {
      iter++;
      status = gsl_multiroot_fsolver_iterate(solver);

#ifdef ENABLE_VERBOSE
      print_state(solver, iter); // @todo implement
#endif

      if (status)   // check if solver is stuck
         break;

      status = gsl_multiroot_test_residual(solver->f, precision_goal);
   } while (status == GSL_CONTINUE && iter < get_max_iterations());

#ifdef ENABLE_VERBOSE
   printf("\tRoot_finder status = %s\n", gsl_strerror(status));
#endif

   set_model_parameters(to_eigen_array(solver->x));

   gsl_vector_free(x_init);
   gsl_multiroot_fsolver_free(solver);

   return status;
}

/**
 * Sanity checks the models and boundary condtitions.
 */
void RGFlow<Root>::check_setup() const
{
   for (size_t m = 0; m < models.size(); ++m) {
      const TModel* model = models[m];
      if (!model->model) {
         std::stringstream message;
         message << "RGFlow<Root>::Error: model pointer ["
                 << m << "] is NULL";
         throw SetupError(message.str());
      }
   }

   if (!convergence_tester) {
      throw SetupError("RGFlow<Root>::Error: convergence tester must "
                       "not be NULL");
   }
}

void RGFlow<Root>::clear_problems()
{
   VERBOSE_MSG("> clearing problems ...");

   const size_t number_of_models = models.size();
   for (size_t m = 0; m < number_of_models; ++m) {
      TModel* model = models[m];
      model->model->clear_problems();
   }
}

void RGFlow<Root>::delete_models()
{
   for_each(models.begin(), models.end(), Delete_object());
}

/**
 * Returns the maximum number of iterations set in the convergence
 * tester.
 */
unsigned int RGFlow<Root>::get_max_iterations() const
{
   return convergence_tester->max_iterations();
}

/**
 * Does the initial guess by calling the guess() method of the initial
 * guesser (if given).
 */
void RGFlow<Root>::initial_guess()
{
   if (initial_guesser)
      initial_guesser->guess();
}

/**
 * Returns the value returned by the accuracy_goal_reached() method of
 * the convergence tester.
 */
bool RGFlow<Root>::accuracy_goal_reached() const
{
   return convergence_tester->accuracy_goal_reached();
}

void RGFlow<Root>::add_model(Root_model* model,
                             const std::vector<Constraint<Root>*>& constraints)
{
   TModel* new_model = new TModel(model, constraints);

   for (Constraint_container::iterator it = new_model->constraints.begin(),
           end = new_model->constraints.end(); it != end; ++it)
      (*it)->set_model(model);

   models.push_back(new_model);
}

void RGFlow<Root>::set_convergence_tester(Convergence_tester<Root>* convergence_tester_)
{
   convergence_tester = convergence_tester_;
}

void RGFlow<Root>::set_initial_guesser(Initial_guesser<Root>* ig)
{
   initial_guesser = ig;
}

std::size_t RGFlow<Root>::get_dimension() const
{
   std::size_t dimension = 0;

   for (size_t m = 0; m < models.size(); ++m)
      dimension += models[m]->model->get_number_of_parameters();

   return dimension;
}

void RGFlow<Root>::starting_point(gsl_vector* x_init)
{
   std::size_t count = 0;

   for (std::size_t m = 0; m < models.size(); m++) {
      const Eigen::ArrayXd parameters(models[m]->model->get_parameters());

      for (std::size_t p = 0; p < parameters.rows(); p++, count++)
         gsl_vector_set(x_init, count, parameters(p));
   }

   assert(count == get_dimension());
}

Eigen::ArrayXd RGFlow<Root>::calculate_function(const Eigen::ArrayXd& x)
{
   Eigen::ArrayXd f(get_dimension());
   std::size_t offset = 0;

   for (std::size_t m = 0; m < models.size(); m++) {
      const Constraint_container constraints(models[m]->constraints);

      for (std::size_t c = 0; c < constraints.size(); c++) {
         const Eigen::ArrayXd func(constraints[c]->calculate_function(x));
         const std::size_t rows = func.rows();

         assert(func.rows() == constraints[c]->get_number_of_fixed_parameters());

         for (std::size_t i = 0; i < rows; i++)
            f(offset++) = func(i);
      }
   }

   assert(offset == get_dimension());

   return f;
}

void RGFlow<Root>::set_model_parameters(const Eigen::ArrayXd& p)
{
   std::size_t offset = 0;

   for (std::size_t m = 0; m < models.size(); m++) {
      const std::size_t dim = models[m]->model->get_number_of_parameters();
      models[m]->model->set_parameters(p.block(offset, 0, dim, 1));
      offset += dim;
   }

   assert(offset == get_dimension());
}

Eigen::ArrayXd RGFlow<Root>::to_eigen_array(const gsl_vector* v)
{
   const std::size_t dim = v->size;
   Eigen::ArrayXd xa(dim);

   for (std::size_t i = 0; i < dim; i++)
      xa(i) = gsl_vector_get(v, i);

   return xa;
}

gsl_vector* RGFlow<Root>::to_gsl_vector(const Eigen::ArrayXd& v)
{
   const std::size_t dim = v.rows();
   gsl_vector* result = gsl_vector_alloc(dim);

   for (std::size_t i = 0; i < dim; i++)
      gsl_vector_set(result, i, v(i));

   return result;
}

void RGFlow<Root>::to_gsl_vector(const Eigen::ArrayXd& v, gsl_vector* dst)
{
   const std::size_t dim = v.rows();

   assert(dim == dst->size);

   for (std::size_t i = 0; i < dim; i++)
      gsl_vector_set(dst, i, v(i));
}

} // namespace flexiblesusy