GslMultiRoot.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <cmath>
8 #include <cstddef>
9 #include <gsl/gsl_errno.h>
10 #include <gsl/gsl_matrix_double.h>
11 #include <gsl/gsl_multiroots.h>
12 #include <gsl/gsl_vector_double.h>
13 #include <ostream>
14 #include <string>
15 
16 #include "ErrorHandling/Error.hpp"
17 #include "ErrorHandling/Exceptions.hpp"
18 #include "Parallel/Printf.hpp"
19 #include "Utilities/Gsl.hpp"
20 #include "Utilities/Requires.hpp"
21 #include "Utilities/TypeTraits.hpp"
22 
23 namespace RootFinder {
24 namespace gsl_multiroot_detail {
25 template <size_t Dim, typename Solver>
26 void print_state(const size_t iteration_number, const Solver* const solver,
27  const bool print_header = false) noexcept {
28  if (print_header) {
29  Parallel::printf("Iter\t");
30  for (size_t i = 0; i < Dim; ++i) {
31  Parallel::printf(" x[%u]\t", i);
32  }
33  for (size_t i = 0; i < Dim; ++i) {
34  Parallel::printf(" f[%u]\t", i);
35  }
36  Parallel::printf("\n");
37  }
38 
39  Parallel::printf("%u\t", iteration_number);
40  for (size_t i = 0; i < Dim; ++i) {
41  Parallel::printf("%3.4f ", gsl_vector_get(solver->x, i));
42  }
43  for (size_t i = 0; i < Dim; ++i) {
44  Parallel::printf("%1.3e ", gsl_vector_get(solver->f, i));
45  }
46  Parallel::printf("\n");
47 }
48 
49 template <size_t Dim>
50 std::array<double, Dim> gsl_to_std_array(const gsl_vector* const x) noexcept {
51  std::array<double, Dim> input_as_std_array{};
52  for (size_t i = 0; i < Dim; i++) {
53  gsl::at(input_as_std_array, i) = gsl_vector_get(x, i);
54  }
55  return input_as_std_array;
56 }
57 
58 template <size_t Dim>
59 void gsl_vector_set_with_std_array(
60  gsl_vector* const func,
61  const std::array<double, Dim>& result_as_std_array) noexcept {
62  for (size_t i = 0; i < Dim; i++) {
63  gsl_vector_set(func, i, gsl::at(result_as_std_array, i));
64  }
65 }
66 
67 template <size_t Dim>
68 void gsl_matrix_set_with_std_array(
69  gsl_matrix* const matrix,
70  const std::array<std::array<double, Dim>, Dim>& matrix_array) noexcept {
71  for (size_t i = 0; i < Dim; i++) {
72  for (size_t j = 0; j < Dim; j++) {
73  gsl_matrix_set(matrix, i, j, gsl::at(gsl::at(matrix_array, i), j));
74  }
75  }
76 }
77 
78 // The gsl_multiroot_function_fdf expects its functions to be of the form
79 // int (* f) (const gsl_vector * x, void * params, gsl_vector * f).
80 // However, we would like to be able to perform rootfinding on functions
81 // of the form std::array<double, Dim> f(const std::array<double, Dim>& x).
82 // So, we pass the function wrapper below to gsl_multiroot_function_fdf.
83 // In the gsl documentation the third parameter is refered to as "void* params",
84 // referring to the parameters that select out a particular function out of a
85 // family of possible functions described in a class, but we instead pass the
86 // pointer to the entire function object itself here. The type of the function
87 // object is passed through with the Function template parameter.
88 template <size_t Dim, typename Function>
89 int gsl_multirootfunctionfdf_wrapper_f(const gsl_vector* const x,
90  void* const untyped_function_object,
91  gsl_vector* const f_of_x) noexcept {
92  const auto function_object =
93  static_cast<const Function*>(untyped_function_object);
94  gsl_vector_set_with_std_array(
95  f_of_x, function_object->operator()(gsl_to_std_array<Dim>(x)));
96  return GSL_SUCCESS;
97 }
98 
99 template <size_t Dim, typename Function>
100 int gsl_multirootfunctionfdf_wrapper_df(const gsl_vector* const x,
101  void* const untyped_function_object,
102  gsl_matrix* const jacobian) noexcept {
103  const auto function_object =
104  static_cast<const Function*>(untyped_function_object);
105  gsl_matrix_set_with_std_array(
106  jacobian, function_object->jacobian(gsl_to_std_array<Dim>(x)));
107 
108  return GSL_SUCCESS;
109 }
110 
111 template <size_t Dim, typename Function>
112 int gsl_multirootfunctionfdf_wrapper_fdf(const gsl_vector* const x,
113  void* const untyped_function_object,
114  gsl_vector* const f_of_x,
115  gsl_matrix* const jacobian) noexcept {
116  const auto function_object =
117  static_cast<const Function*>(untyped_function_object);
118  const std::array<double, Dim> x_as_std_array = gsl_to_std_array<Dim>(x);
119  gsl_vector_set_with_std_array(f_of_x,
120  function_object->operator()(x_as_std_array));
121  gsl_matrix_set_with_std_array(jacobian,
122  function_object->jacobian(x_as_std_array));
123  return GSL_SUCCESS;
124 }
125 
126 CREATE_IS_CALLABLE(jacobian)
127 } // namespace gsl_multiroot_detail
128 
129 /*!
130  * \ingroup NumericalAlgorithmsGroup
131  * \brief The different options for the verbosity of gsl_multiroot.
132  */
133 enum class Verbosity {
134  /// Do not print anything
135  Silent,
136  /// Print only "success" or "failed" on termination
137  Quiet,
138  /// Print final functions values on termination
139  Verbose,
140  /// Print function values on every iteration
141  Debug
142 };
143 
144 std::ostream& operator<<(std::ostream& /*os*/,
145  const Verbosity& /*verbosity*/) noexcept;
146 
147 /*!
148  * \ingroup NumericalAlgorithmsGroup
149  * \brief The different options for the rootfinding method of gsl_multiroot.
150  *
151  * This enum is for setting the method used the rootfinder.
152  * The precise method used by the gsl rootfinder depends on whether or not the
153  * function passed to it has a callable `jacobian` member function. In the
154  * case where it doesn't, the jacobian is approximated with a finite
155  * difference. For example, if the Method specified is Hybrid, gsl will use
156  * the gsl_multiroot_fdfsolver_hybridj method in the case where a `jacobian`
157  * is provided, and gsl_multiroot_fsolver_hybrid in the case where one isn't.
158  * See
159  * [GSL's documentation for multidimensional
160  * rootfinding](https://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html)
161  * for information on the different methods.
162  * \note gsl does not provide a finite difference version for the modified
163  * Newton method (gsl_multiroot_fdfsolver_gnewton). In the case where a
164  * jacobian is not provided the method used will be a non-modified Newton
165  * method.
166  *
167  */
168 enum class Method {
169  /// Hybrid of Newton's method along with following the gradient direction.
170  /// \note Sometimes Hybrids works only with the Absolute stopping condition.
171  Hybrids,
172  /// "Unscaled version of Hybrids that uses a spherical trust region," see
173  /// GSL documentation for more details.
174  Hybrid,
175  /// If an analytic jacobian is provided, gsl uses a modification of Newton's
176  /// method to improve global convergence. Uses vanilla Newton's method if no
177  /// jacobian is provided.
178  Newton
179 };
180 
181 std::ostream& operator<<(std::ostream& /*os*/,
182  const Method& /*method*/) noexcept;
183 
184 /*!
185  * \ingroup NumericalAlgorithmsGroup
186  * \brief The different options for the convergence criterion of gsl_multiroot.
187  *
188  * See
189  * [GSL's documentation for multidimensional
190  * rootfinding](https://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html)
191  * for information on the different stopping conditions.
192  */
193 enum class StoppingCondition {
194  /// See GSL documentation for gsl_multiroot_test_delta.
196  /// See GSL documentation for gsl_multiroot_test_residual.
197  Absolute
198 };
199 
200 std::ostream& operator<<(std::ostream& /*os*/,
201  const StoppingCondition& /*condition*/) noexcept;
202 
203 namespace gsl_multiroot_detail {
204 template <typename SolverType, typename SolverAlloc, typename SolverSet,
205  typename SolverIterate, typename SolverFree, size_t Dim,
206  typename Function>
207 std::array<double, Dim> gsl_multiroot_impl(
208  Function& f, const std::array<double, Dim>& initial_guess,
209  const double absolute_tolerance, const size_t maximum_iterations,
210  const double relative_tolerance, const Verbosity verbosity,
211  const double maximum_absolute_tolerance, const Method method,
212  const SolverType solver_type, const StoppingCondition condition,
213  const SolverAlloc solver_alloc, const SolverSet solver_set,
214  const SolverIterate solver_iterate, const SolverFree solver_free) {
215  // Check for valid stopping condition:
217  condition != StoppingCondition::Absolute)) {
218  ERROR(
219  "Invalid stopping condition. Has to be either AbsoluteAndRelative"
220  "or Absolute.");
221  }
222 
223  // Supply gsl_root with the initial guess:
224  gsl_vector* const gsl_root = gsl_vector_alloc(Dim);
225  gsl_vector_set_with_std_array(gsl_root, initial_guess);
226  auto* const solver = solver_alloc(solver_type, Dim);
227  solver_set(solver, &f, gsl_root);
228 
229  // Take iterations:
230  int status;
231  size_t iteration_number = 0;
232  do {
233  if (UNLIKELY(verbosity == Verbosity::Debug)) {
234  print_state<Dim>(iteration_number, solver, iteration_number == 0);
235  }
236  iteration_number++;
237  status = solver_iterate(solver);
238  // Check if solver is stuck
239  if (UNLIKELY(status == GSL_ENOPROG)) {
240  if (UNLIKELY(verbosity == Verbosity::Debug)) {
242  "The iteration is not making any progress, preventing the "
243  "algorithm from continuing.");
244  }
245  break;
246  }
247  if (condition == StoppingCondition::AbsoluteAndRelative) {
248  status = gsl_multiroot_test_delta(solver->dx, solver->x,
249  absolute_tolerance, relative_tolerance);
250  } else { // condition is StoppingCondition::Absolute
251  // NOTE: Sometimes hybridsj works only with the test_residual condition
252  status = gsl_multiroot_test_residual(solver->f, absolute_tolerance);
253  }
254  } while (status == GSL_CONTINUE and iteration_number < maximum_iterations);
255  if (UNLIKELY(verbosity == Verbosity::Verbose or
256  verbosity == Verbosity::Debug)) {
257  Parallel::printf("Finished iterating:\n");
258  print_state<Dim>(iteration_number, solver, verbosity == Verbosity::Verbose);
259  }
260  bool success = (status == GSL_SUCCESS);
261  if (UNLIKELY(verbosity != Verbosity::Silent)) {
262  Parallel::printf("\n");
263  if (not success) {
264  const std::string ascii_divider = std::string(70, '#');
265  const std::string failure_message =
266  ascii_divider + "\n\t\tWARNING: Root Finding FAILED\n" +
267  ascii_divider;
268  Parallel::printf("%s\n", failure_message);
269  } else {
270  Parallel::printf("Root finder converged.\n");
271  }
272  }
273  // If maximum_absolute_tolerance is given, return success = true
274  // as long as maximum_absolute_tolerance is achieved even if the
275  // root finder doesn't converge.
276  bool success_with_tolerance = true;
277  bool failed_root_is_forgiven = false;
278  if (not success and maximum_absolute_tolerance > 0.0) {
279  for (size_t i = 0; i < Dim; ++i) {
280  if (fabs(gsl_vector_get(solver->f, i)) > maximum_absolute_tolerance) {
281  success_with_tolerance = false;
282  }
283  }
284  failed_root_is_forgiven = success_with_tolerance;
285  }
286  if (success_with_tolerance and maximum_absolute_tolerance > 0.0) {
287  success = true;
288  }
289  if (UNLIKELY(failed_root_is_forgiven and (verbosity == Verbosity::Verbose or
290  verbosity == Verbosity::Debug))) {
292  "The failed root was forgiven as each component was found to be under "
293  "maximum_absolute_tolerance %f",
294  maximum_absolute_tolerance);
295  }
296 
297  if (UNLIKELY(not success)) {
298  std::stringstream error_message;
299  error_message << "The root find failed and was not forgiven. An exception "
300  "has been thrown.\n"
301  << "The gsl error returned is: " << gsl_strerror(status)
302  << "\n"
303  << "Verbosity: " << verbosity << "\n"
304  << "Method: " << method << "\n"
305  << "StoppingCondition: " << condition << "\n"
306  << "Maximum absolute tolerance: "
307  << maximum_absolute_tolerance << "\n"
308  << "Absolute tolerance: " << absolute_tolerance << "\n"
309  << "Relative tolerance: " << relative_tolerance << "\n"
310  << "Maximum number of iterations: " << maximum_iterations
311  << "\n"
312  << "Number of iterations reached: " << iteration_number
313  << "\n"
314  << "The last value of f in the root solver is:\n";
315  for (size_t i = 0; i < Dim; i++) {
316  error_message << gsl_vector_get(solver->f, i) << "\n";
317  }
318  error_message << "The last value of x in the root solver is:\n";
319  for (size_t i = 0; i < Dim; i++) {
320  error_message << gsl_vector_get(solver->x, i) << "\n";
321  }
322  error_message << "The last value of dx in the root solver is:\n";
323  for (size_t i = 0; i < Dim; i++) {
324  error_message << gsl_vector_get(solver->dx, i) << "\n";
325  }
326 
327  if (UNLIKELY(verbosity == Verbosity::Debug)) {
328  Parallel::printf("Error: %s\n", gsl_strerror(status));
329  if (iteration_number >= maximum_iterations) {
331  "The number of iterations (%d) has reached the maximum number of "
332  "iterations (%d)\n",
333  iteration_number, maximum_iterations);
334  } else {
336  "The number of iterations (%d) failed to reach the maximum number of "
337  "iterations (%d)\n",
338  iteration_number, maximum_iterations);
339  }
340  }
341  throw convergence_error(error_message.str());
342  }
343 
344  // Store the converged root in result
345  std::array<double, Dim> result = gsl_to_std_array<Dim>(solver->x);
346  solver_free(solver);
347  gsl_vector_free(gsl_root);
348 
349  return result;
350 }
351 
352 void print_rootfinding_parameters(Method method, double absolute_tolerance,
353  double relative_tolerance,
354  double maximum_absolute_tolerance,
355  StoppingCondition condition) noexcept;
356 } // namespace gsl_multiroot_detail
357 
358 // @{
359 /*!
360  * \ingroup NumericalAlgorithmsGroup
361  * \brief A multidimensional root finder supporting Newton and Hybrid
362  * methods, as well as modified methods based on these.
363  *
364  * This root finder accepts function objects with and without a callable
365  * `jacobian` member function. This member function both accepts and
366  * returns a `std::array<double, Dim>`, the dimension of the domain and range
367  * of the function the root find is being performed on. Whether the jacobian
368  * is provided determines the details of the implementation of the
369  * root-finding method that is selected by the user using the Method enum.
370  * That is, whether the jacobian is computed analytically via the `jacobian`
371  * member function, or whether the jacobian is computed numerically via a
372  * finite difference approximation.
373  * \note GSL does not provide a finite difference version of its modified
374  * Newton method, so the unmodified one is used instead when the user
375  * uses the Method::Newton method.
376  *
377  * The user can select one of two possible criteria for convergence,
378  * StoppingCondition::Absolute, where the sum of the absolute values of the
379  * components of the residual vector f are compared against the value
380  * provided to `absolute_tolerance`, and
381  * StoppingCondition::AbsoluteAndRelative, where the size of the most recent
382  * step taken in the root-finding iteration is compared against
383  * `absolute_tolerance` + `relative_tolerance` * |x_i|, for each component.
384  * In either case, a `maximum_absolute_tolerance` may be specified if the user
385  * anticipates that the convergence criterion specified with StoppingCondition
386  * will be too strict for a few points out of a population of points found with
387  * a sequence of root finds.
388  *
389  * See
390  * [GSL's documentation for multidimensional
391  * rootfinding](https://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html)
392  * for reference.
393  *
394  * \param func Function whose root is to be found.
395  * \param initial_guess Contains initial guess.
396  * \param absolute_tolerance The absolute tolerance.
397  * \param maximum_iterations The maximum number of iterations.
398  * \param relative_tolerance The relative tolerance.
399  * \param verbosity Whether to print diagnostic messages.
400  * \param maximum_absolute_tolerance Acceptable absolute tolerance when
401  * root finder doesn't converge.
402  * You may wish to use this if there
403  * are only a few "problematic" points where
404  * it is difficult to do a precise root find.
405  * \param method The method to use. See the documentation for the Method enum.
406  * \param condition The convergence condition to use. See the documentation
407  * for the StoppingCondition enum.
408  */
409 template <size_t Dim, typename Function,
410  Requires<gsl_multiroot_detail::is_jacobian_callable_v<
411  Function, std::array<double, Dim>>> = nullptr>
412 std::array<double, Dim> gsl_multiroot(
413  const Function& func, const std::array<double, Dim>& initial_guess,
414  const double absolute_tolerance, const size_t maximum_iterations,
415  const double relative_tolerance = 0.0,
416  const Verbosity verbosity = Verbosity::Silent,
417  const double maximum_absolute_tolerance = 0.0,
418  const Method method = Method::Newton,
420  gsl_multiroot_function_fdf gsl_func = {
421  &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_f<Dim, Function>,
422  &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_df<Dim, Function>,
423  &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_fdf<Dim,
424  Function>,
425  Dim, const_cast<Function*>(&func)}; //NOLINT
426 
427  // Set up method for solver:
428  const gsl_multiroot_fdfsolver_type* solver_type;
429  if (method == Method::Newton) {
430  solver_type = gsl_multiroot_fdfsolver_gnewton;
431  } else if (method == Method::Hybrids) {
432  solver_type = gsl_multiroot_fdfsolver_hybridsj;
433  } else if (method == Method::Hybrid) {
434  solver_type = gsl_multiroot_fdfsolver_hybridj;
435  } else {
436  ERROR(
437  "Invalid method. Has to be one of Newton, Hybrids or "
438  "Hybrid.");
439  }
440  // Print initial parameters
441  if (UNLIKELY(verbosity == Verbosity::Verbose or
442  verbosity == Verbosity::Debug)) {
443  gsl_multiroot_detail::print_rootfinding_parameters(
444  method, absolute_tolerance, relative_tolerance,
445  maximum_absolute_tolerance, condition);
446  }
447  return gsl_multiroot_detail::gsl_multiroot_impl(
448  gsl_func, initial_guess, absolute_tolerance, maximum_iterations,
449  relative_tolerance, verbosity, maximum_absolute_tolerance, method,
450  solver_type, condition, &gsl_multiroot_fdfsolver_alloc,
451  &gsl_multiroot_fdfsolver_set, &gsl_multiroot_fdfsolver_iterate,
452  &gsl_multiroot_fdfsolver_free);
453 }
454 
455 template <size_t Dim, typename Function,
456  Requires<not gsl_multiroot_detail::is_jacobian_callable_v<
457  Function, std::array<double, Dim>>> = nullptr>
458 std::array<double, Dim> gsl_multiroot(
459  const Function& func, const std::array<double, Dim>& initial_guess,
460  const double absolute_tolerance, const size_t maximum_iterations,
461  const double relative_tolerance = 0.0,
462  const Verbosity verbosity = Verbosity::Silent,
463  const double maximum_absolute_tolerance = 0.0,
464  const Method method = Method::Newton,
466  gsl_multiroot_function gsl_func = {
467  &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_f<Dim, Function>,
468  Dim, const_cast<Function*>(&func)}; // NOLINT
469 
470  // Set up method for solver:
471  const gsl_multiroot_fsolver_type* solver_type;
472  if (method == Method::Newton) {
473  solver_type = gsl_multiroot_fsolver_dnewton;
474  } else if (method == Method::Hybrids) {
475  solver_type = gsl_multiroot_fsolver_hybrids;
476  } else if (method == Method::Hybrid) {
477  solver_type = gsl_multiroot_fsolver_hybrid;
478  } else {
479  ERROR(
480  "Invalid method. Has to be one of Newton, Hybrids or "
481  "Hybrid.");
482  }
483  // Print initial parameters
484  if (UNLIKELY(verbosity == Verbosity::Verbose or
485  verbosity == Verbosity::Debug)) {
486  gsl_multiroot_detail::print_rootfinding_parameters(
487  method, absolute_tolerance, relative_tolerance,
488  maximum_absolute_tolerance, condition);
489  }
490  return gsl_multiroot_detail::gsl_multiroot_impl(
491  gsl_func, initial_guess, absolute_tolerance, maximum_iterations,
492  relative_tolerance, verbosity, maximum_absolute_tolerance, method,
493  solver_type, condition, &gsl_multiroot_fsolver_alloc,
494  &gsl_multiroot_fsolver_set, &gsl_multiroot_fsolver_iterate,
495  &gsl_multiroot_fsolver_free);
496 }
497 // @}
498 } // namespace RootFinder
If an analytic jacobian is provided, gsl uses a modification of Newton&#39;s method to improve global con...
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:35
#define CREATE_IS_CALLABLE(METHOD_NAME)
Generate a type trait to check if a class has a member function that can be invoked with arguments of...
Definition: TypeTraits.hpp:870
Definition: GslMultiRoot.cpp:10
#define UNLIKELY(x)
Definition: Gsl.hpp:72
See GSL documentation for gsl_multiroot_test_residual.
StoppingCondition
The different options for the convergence criterion of gsl_multiroot.
Definition: GslMultiRoot.hpp:193
Defines the type alias Requires.
Print final functions values on termination.
Method
The different options for the rootfinding method of gsl_multiroot.
Definition: GslMultiRoot.hpp:168
Hybrid of Newton&#39;s method along with following the gradient direction.
Do not print anything.
Exception indicating convergence failure.
Definition: Exceptions.hpp:10
"Unscaled version of Hybrids that uses a spherical trust region," see GSL documentation for more deta...
Defines Parallel::printf for writing to stdout.
std::array< double, Dim > gsl_multiroot(const Function &func, const std::array< double, Dim > &initial_guess, const double absolute_tolerance, const size_t maximum_iterations, const double relative_tolerance=0.0, const Verbosity verbosity=Verbosity::Silent, const double maximum_absolute_tolerance=0.0, const Method method=Method::Newton, const StoppingCondition condition=StoppingCondition::Absolute)
A multidimensional root finder supporting Newton and Hybrid methods, as well as modified methods base...
Definition: GslMultiRoot.hpp:412
See GSL documentation for gsl_multiroot_test_delta.
Defines functions and classes from the GSL.
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t ...
Definition: Requires.hpp:67
Print only "success" or "failed" on termination.
Print function values on every iteration.
Defines macro ERROR.
Defines type traits, some of which are future STL type_traits header.
void printf(const std::string &format, Args &&... args)
Print an atomic message to stdout with C printf usage.
Definition: Printf.hpp:100
Verbosity
The different options for the verbosity of gsl_multiroot.
Definition: GslMultiRoot.hpp:133
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid...
Definition: Gsl.hpp:124