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