Line data Source code
1 0 : // 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 "Parallel/Printf.hpp"
17 : #include "Utilities/ErrorHandling/Error.hpp"
18 : #include "Utilities/ErrorHandling/Exceptions.hpp"
19 : #include "Utilities/Gsl.hpp"
20 : #include "Utilities/Requires.hpp"
21 : #include "Utilities/TypeTraits/CreateIsCallable.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 : CREATE_IS_CALLABLE_V(jacobian)
128 : } // namespace gsl_multiroot_detail
129 :
130 : /*!
131 : * \ingroup NumericalAlgorithmsGroup
132 : * \brief The different options for the verbosity of gsl_multiroot.
133 : */
134 1 : enum class Verbosity {
135 : /// Do not print anything
136 : Silent,
137 : /// Print only "success" or "failed" on termination
138 : Quiet,
139 : /// Print final functions values on termination
140 : Verbose,
141 : /// Print function values on every iteration
142 : Debug
143 : };
144 :
145 : std::ostream& operator<<(std::ostream& /*os*/,
146 : const Verbosity& /*verbosity*/) noexcept;
147 :
148 : /*!
149 : * \ingroup NumericalAlgorithmsGroup
150 : * \brief The different options for the rootfinding method of gsl_multiroot.
151 : *
152 : * This enum is for setting the method used the rootfinder.
153 : * The precise method used by the gsl rootfinder depends on whether or not the
154 : * function passed to it has a callable `jacobian` member function. In the
155 : * case where it doesn't, the jacobian is approximated with a finite
156 : * difference. For example, if the Method specified is Hybrid, gsl will use
157 : * the gsl_multiroot_fdfsolver_hybridj method in the case where a `jacobian`
158 : * is provided, and gsl_multiroot_fsolver_hybrid in the case where one isn't.
159 : * See
160 : * [GSL's documentation for multidimensional
161 : * rootfinding](https://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html)
162 : * for information on the different methods.
163 : * \note gsl does not provide a finite difference version for the modified
164 : * Newton method (gsl_multiroot_fdfsolver_gnewton). In the case where a
165 : * jacobian is not provided the method used will be a non-modified Newton
166 : * method.
167 : *
168 : */
169 1 : enum class Method {
170 : /// Hybrid of Newton's method along with following the gradient direction.
171 : /// \note Sometimes Hybrids works only with the Absolute stopping condition.
172 : Hybrids,
173 : /// "Unscaled version of Hybrids that uses a spherical trust region," see
174 : /// GSL documentation for more details.
175 : Hybrid,
176 : /// If an analytic jacobian is provided, gsl uses a modification of Newton's
177 : /// method to improve global convergence. Uses vanilla Newton's method if no
178 : /// jacobian is provided.
179 : Newton
180 : };
181 :
182 : std::ostream& operator<<(std::ostream& /*os*/,
183 : const Method& /*method*/) noexcept;
184 :
185 : /*!
186 : * \ingroup NumericalAlgorithmsGroup
187 : * \brief The different options for the convergence criterion of gsl_multiroot.
188 : *
189 : * See
190 : * [GSL's documentation for multidimensional
191 : * rootfinding](https://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html)
192 : * for information on the different stopping conditions.
193 : */
194 1 : enum class StoppingCondition {
195 : /// See GSL documentation for gsl_multiroot_test_delta.
196 : AbsoluteAndRelative,
197 : /// See GSL documentation for gsl_multiroot_test_residual.
198 : Absolute
199 : };
200 :
201 : std::ostream& operator<<(std::ostream& /*os*/,
202 : const StoppingCondition& /*condition*/) noexcept;
203 :
204 : namespace gsl_multiroot_detail {
205 : template <typename SolverType, typename SolverAlloc, typename SolverSet,
206 : typename SolverIterate, typename SolverFree, size_t Dim,
207 : typename Function>
208 : std::array<double, Dim> gsl_multiroot_impl(
209 : Function& f, const std::array<double, Dim>& initial_guess,
210 : const double absolute_tolerance, const size_t maximum_iterations,
211 : const double relative_tolerance, const Verbosity verbosity,
212 : const double maximum_absolute_tolerance, const Method method,
213 : const SolverType solver_type, const StoppingCondition condition,
214 : const SolverAlloc solver_alloc, const SolverSet solver_set,
215 : const SolverIterate solver_iterate, const SolverFree solver_free) {
216 : // Check for valid stopping condition:
217 : if (UNLIKELY(condition != StoppingCondition::AbsoluteAndRelative and
218 : condition != StoppingCondition::Absolute)) {
219 : ERROR(
220 : "Invalid stopping condition. Has to be either AbsoluteAndRelative"
221 : "or Absolute.");
222 : }
223 :
224 : // Supply gsl_root with the initial guess:
225 : gsl_vector* const gsl_root = gsl_vector_alloc(Dim);
226 : gsl_vector_set_with_std_array(gsl_root, initial_guess);
227 : auto* const solver = solver_alloc(solver_type, Dim);
228 : solver_set(solver, &f, gsl_root);
229 :
230 : // Take iterations:
231 : int status;
232 : size_t iteration_number = 0;
233 : do {
234 : if (UNLIKELY(verbosity == Verbosity::Debug)) {
235 : print_state<Dim>(iteration_number, solver, iteration_number == 0);
236 : }
237 : iteration_number++;
238 : status = solver_iterate(solver);
239 : // Check if solver is stuck
240 : if (UNLIKELY(status == GSL_ENOPROG)) {
241 : if (UNLIKELY(verbosity == Verbosity::Debug)) {
242 : Parallel::printf(
243 : "The iteration is not making any progress, preventing the "
244 : "algorithm from continuing.");
245 : }
246 : break;
247 : }
248 : if (condition == StoppingCondition::AbsoluteAndRelative) {
249 : status = gsl_multiroot_test_delta(solver->dx, solver->x,
250 : absolute_tolerance, relative_tolerance);
251 : } else { // condition is StoppingCondition::Absolute
252 : // NOTE: Sometimes hybridsj works only with the test_residual condition
253 : status = gsl_multiroot_test_residual(solver->f, absolute_tolerance);
254 : }
255 : } while (status == GSL_CONTINUE and iteration_number < maximum_iterations);
256 : if (UNLIKELY(verbosity == Verbosity::Verbose or
257 : verbosity == Verbosity::Debug)) {
258 : Parallel::printf("Finished iterating:\n");
259 : print_state<Dim>(iteration_number, solver, verbosity == Verbosity::Verbose);
260 : }
261 : bool success = (status == GSL_SUCCESS);
262 : if (UNLIKELY(verbosity != Verbosity::Silent)) {
263 : Parallel::printf("\n");
264 : if (not success) {
265 : const std::string ascii_divider = std::string(70, '#');
266 : const std::string failure_message =
267 : ascii_divider + "\n\t\tWARNING: Root Finding FAILED\n" +
268 : ascii_divider;
269 : Parallel::printf("%s\n", failure_message);
270 : } else {
271 : Parallel::printf("Root finder converged.\n");
272 : }
273 : }
274 : // If maximum_absolute_tolerance is given, return success = true
275 : // as long as maximum_absolute_tolerance is achieved even if the
276 : // root finder doesn't converge.
277 : bool success_with_tolerance = true;
278 : bool failed_root_is_forgiven = false;
279 : if (not success and maximum_absolute_tolerance > 0.0) {
280 : for (size_t i = 0; i < Dim; ++i) {
281 : if (fabs(gsl_vector_get(solver->f, i)) > maximum_absolute_tolerance) {
282 : success_with_tolerance = false;
283 : }
284 : }
285 : failed_root_is_forgiven = success_with_tolerance;
286 : }
287 : if (success_with_tolerance and maximum_absolute_tolerance > 0.0) {
288 : success = true;
289 : }
290 : if (UNLIKELY(failed_root_is_forgiven and (verbosity == Verbosity::Verbose or
291 : verbosity == Verbosity::Debug))) {
292 : Parallel::printf(
293 : "The failed root was forgiven as each component was found to be under "
294 : "maximum_absolute_tolerance %f",
295 : maximum_absolute_tolerance);
296 : }
297 :
298 : if (UNLIKELY(not success)) {
299 : std::stringstream error_message;
300 : error_message << "The root find failed and was not forgiven. An exception "
301 : "has been thrown.\n"
302 : << "The gsl error returned is: " << gsl_strerror(status)
303 : << "\n"
304 : << "Verbosity: " << verbosity << "\n"
305 : << "Method: " << method << "\n"
306 : << "StoppingCondition: " << condition << "\n"
307 : << "Maximum absolute tolerance: "
308 : << maximum_absolute_tolerance << "\n"
309 : << "Absolute tolerance: " << absolute_tolerance << "\n"
310 : << "Relative tolerance: " << relative_tolerance << "\n"
311 : << "Maximum number of iterations: " << maximum_iterations
312 : << "\n"
313 : << "Number of iterations reached: " << iteration_number
314 : << "\n"
315 : << "The last value of f in the root solver is:\n";
316 : for (size_t i = 0; i < Dim; i++) {
317 : error_message << gsl_vector_get(solver->f, i) << "\n";
318 : }
319 : error_message << "The last value of x in the root solver is:\n";
320 : for (size_t i = 0; i < Dim; i++) {
321 : error_message << gsl_vector_get(solver->x, i) << "\n";
322 : }
323 : error_message << "The last value of dx in the root solver is:\n";
324 : for (size_t i = 0; i < Dim; i++) {
325 : error_message << gsl_vector_get(solver->dx, i) << "\n";
326 : }
327 :
328 : if (UNLIKELY(verbosity == Verbosity::Debug)) {
329 : Parallel::printf("Error: %s\n", gsl_strerror(status));
330 : if (iteration_number >= maximum_iterations) {
331 : Parallel::printf(
332 : "The number of iterations (%d) has reached the maximum number of "
333 : "iterations (%d)\n",
334 : iteration_number, maximum_iterations);
335 : } else {
336 : Parallel::printf(
337 : "The number of iterations (%d) failed to reach the maximum number of "
338 : "iterations (%d)\n",
339 : iteration_number, maximum_iterations);
340 : }
341 : }
342 : throw convergence_error(error_message.str());
343 : }
344 :
345 : // Store the converged root in result
346 : std::array<double, Dim> result = gsl_to_std_array<Dim>(solver->x);
347 : solver_free(solver);
348 : gsl_vector_free(gsl_root);
349 :
350 : return result;
351 : }
352 :
353 : void print_rootfinding_parameters(Method method, double absolute_tolerance,
354 : double relative_tolerance,
355 : double maximum_absolute_tolerance,
356 : StoppingCondition condition) noexcept;
357 : } // namespace gsl_multiroot_detail
358 :
359 : // @{
360 : /*!
361 : * \ingroup NumericalAlgorithmsGroup
362 : * \brief A multidimensional root finder supporting Newton and Hybrid
363 : * methods, as well as modified methods based on these.
364 : *
365 : * This root finder accepts function objects with and without a callable
366 : * `jacobian` member function. This member function both accepts and
367 : * returns a `std::array<double, Dim>`, the dimension of the domain and range
368 : * of the function the root find is being performed on. Whether the jacobian
369 : * is provided determines the details of the implementation of the
370 : * root-finding method that is selected by the user using the Method enum.
371 : * That is, whether the jacobian is computed analytically via the `jacobian`
372 : * member function, or whether the jacobian is computed numerically via a
373 : * finite difference approximation.
374 : * \note GSL does not provide a finite difference version of its modified
375 : * Newton method, so the unmodified one is used instead when the user
376 : * uses the Method::Newton method.
377 : *
378 : * The user can select one of two possible criteria for convergence,
379 : * StoppingCondition::Absolute, where the sum of the absolute values of the
380 : * components of the residual vector f are compared against the value
381 : * provided to `absolute_tolerance`, and
382 : * StoppingCondition::AbsoluteAndRelative, where the size of the most recent
383 : * step taken in the root-finding iteration is compared against
384 : * `absolute_tolerance` + `relative_tolerance` * |x_i|, for each component.
385 : * In either case, a `maximum_absolute_tolerance` may be specified if the user
386 : * anticipates that the convergence criterion specified with StoppingCondition
387 : * will be too strict for a few points out of a population of points found with
388 : * a sequence of root finds.
389 : *
390 : * See
391 : * [GSL's documentation for multidimensional
392 : * rootfinding](https://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html)
393 : * for reference.
394 : *
395 : * \param func Function whose root is to be found.
396 : * \param initial_guess Contains initial guess.
397 : * \param absolute_tolerance The absolute tolerance.
398 : * \param maximum_iterations The maximum number of iterations.
399 : * \param relative_tolerance The relative tolerance.
400 : * \param verbosity Whether to print diagnostic messages.
401 : * \param maximum_absolute_tolerance Acceptable absolute tolerance when
402 : * root finder doesn't converge.
403 : * You may wish to use this if there
404 : * are only a few "problematic" points where
405 : * it is difficult to do a precise root find.
406 : * \param method The method to use. See the documentation for the Method enum.
407 : * \param condition The convergence condition to use. See the documentation
408 : * for the StoppingCondition enum.
409 : */
410 : template <size_t Dim, typename Function,
411 : Requires<gsl_multiroot_detail::is_jacobian_callable_v<
412 : Function, std::array<double, Dim>>> = nullptr>
413 1 : std::array<double, Dim> gsl_multiroot(
414 : const Function& func, const std::array<double, Dim>& initial_guess,
415 : const double absolute_tolerance, const size_t maximum_iterations,
416 : const double relative_tolerance = 0.0,
417 : const Verbosity verbosity = Verbosity::Silent,
418 : const double maximum_absolute_tolerance = 0.0,
419 : const Method method = Method::Newton,
420 : const StoppingCondition condition = StoppingCondition::Absolute) {
421 : gsl_multiroot_function_fdf gsl_func = {
422 : &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_f<Dim, Function>,
423 : &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_df<Dim, Function>,
424 : &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_fdf<Dim,
425 : Function>,
426 : Dim, const_cast<Function*>(&func)}; //NOLINT
427 :
428 : // Set up method for solver:
429 : const gsl_multiroot_fdfsolver_type* solver_type;
430 : if (method == Method::Newton) {
431 : solver_type = gsl_multiroot_fdfsolver_gnewton;
432 : } else if (method == Method::Hybrids) {
433 : solver_type = gsl_multiroot_fdfsolver_hybridsj;
434 : } else if (method == Method::Hybrid) {
435 : solver_type = gsl_multiroot_fdfsolver_hybridj;
436 : } else {
437 : ERROR(
438 : "Invalid method. Has to be one of Newton, Hybrids or "
439 : "Hybrid.");
440 : }
441 : // Print initial parameters
442 : if (UNLIKELY(verbosity == Verbosity::Verbose or
443 : verbosity == Verbosity::Debug)) {
444 : gsl_multiroot_detail::print_rootfinding_parameters(
445 : method, absolute_tolerance, relative_tolerance,
446 : maximum_absolute_tolerance, condition);
447 : }
448 : return gsl_multiroot_detail::gsl_multiroot_impl(
449 : gsl_func, initial_guess, absolute_tolerance, maximum_iterations,
450 : relative_tolerance, verbosity, maximum_absolute_tolerance, method,
451 : solver_type, condition, &gsl_multiroot_fdfsolver_alloc,
452 : &gsl_multiroot_fdfsolver_set, &gsl_multiroot_fdfsolver_iterate,
453 : &gsl_multiroot_fdfsolver_free);
454 : }
455 :
456 : template <size_t Dim, typename Function,
457 : Requires<not gsl_multiroot_detail::is_jacobian_callable_v<
458 : Function, std::array<double, Dim>>> = nullptr>
459 : std::array<double, Dim> gsl_multiroot(
460 : const Function& func, const std::array<double, Dim>& initial_guess,
461 : const double absolute_tolerance, const size_t maximum_iterations,
462 : const double relative_tolerance = 0.0,
463 : const Verbosity verbosity = Verbosity::Silent,
464 : const double maximum_absolute_tolerance = 0.0,
465 : const Method method = Method::Newton,
466 : const StoppingCondition condition = StoppingCondition::Absolute) {
467 : gsl_multiroot_function gsl_func = {
468 : &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_f<Dim, Function>,
469 : Dim, const_cast<Function*>(&func)}; // NOLINT
470 :
471 : // Set up method for solver:
472 : const gsl_multiroot_fsolver_type* solver_type;
473 : if (method == Method::Newton) {
474 : solver_type = gsl_multiroot_fsolver_dnewton;
475 : } else if (method == Method::Hybrids) {
476 : solver_type = gsl_multiroot_fsolver_hybrids;
477 : } else if (method == Method::Hybrid) {
478 : solver_type = gsl_multiroot_fsolver_hybrid;
479 : } else {
480 : ERROR(
481 : "Invalid method. Has to be one of Newton, Hybrids or "
482 : "Hybrid.");
483 : }
484 : // Print initial parameters
485 : if (UNLIKELY(verbosity == Verbosity::Verbose or
486 : verbosity == Verbosity::Debug)) {
487 : gsl_multiroot_detail::print_rootfinding_parameters(
488 : method, absolute_tolerance, relative_tolerance,
489 : maximum_absolute_tolerance, condition);
490 : }
491 : return gsl_multiroot_detail::gsl_multiroot_impl(
492 : gsl_func, initial_guess, absolute_tolerance, maximum_iterations,
493 : relative_tolerance, verbosity, maximum_absolute_tolerance, method,
494 : solver_type, condition, &gsl_multiroot_fsolver_alloc,
495 : &gsl_multiroot_fsolver_set, &gsl_multiroot_fsolver_iterate,
496 : &gsl_multiroot_fsolver_free);
497 : }
498 : // @}
499 : } // namespace RootFinder
|
