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/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)) {
297 : Parallel::printf("Finished iterating:\n");
298 : print_state<Dim>(iteration_number, *solver,
299 : verbosity >= Verbosity::Verbose);
300 : }
301 : bool success = (status == GSL_SUCCESS);
302 : if (UNLIKELY(verbosity > Verbosity::Silent)) {
303 : Parallel::printf("\n");
304 : if (not success) {
305 : const std::string ascii_divider = std::string(70, '#');
306 : const std::string failure_message =
307 : ascii_divider + "\n\t\tWARNING: Root Finding FAILED\n" +
308 : ascii_divider;
309 : Parallel::printf("%s\n", failure_message);
310 : } else {
311 : Parallel::printf("Root finder converged.\n");
312 : }
313 : }
314 : // If maximum_absolute_tolerance is given, return success = true
315 : // as long as maximum_absolute_tolerance is achieved even if the
316 : // root finder doesn't converge.
317 : bool success_with_tolerance = true;
318 : bool failed_root_is_forgiven = false;
319 : if (not success and maximum_absolute_tolerance > 0.0) {
320 : for (size_t i = 0; i < Dim; ++i) {
321 : if (fabs(gsl_vector_get(solver->f, i)) > maximum_absolute_tolerance) {
322 : success_with_tolerance = false;
323 : }
324 : }
325 : failed_root_is_forgiven = success_with_tolerance;
326 : }
327 : if (success_with_tolerance and maximum_absolute_tolerance > 0.0) {
328 : success = true;
329 : }
330 : if (UNLIKELY(failed_root_is_forgiven and verbosity >= Verbosity::Verbose)) {
331 : Parallel::printf(
332 : "The failed root was forgiven as each component was found to be under "
333 : "maximum_absolute_tolerance %f",
334 : maximum_absolute_tolerance);
335 : }
336 :
337 : if (UNLIKELY(not success)) {
338 : std::stringstream error_message;
339 : error_message << "The root find failed and was not forgiven. An exception "
340 : "has been thrown.\n"
341 : << "The gsl error returned is: " << gsl_strerror(status)
342 : << "\n"
343 : << "Verbosity: " << verbosity << "\n"
344 : << "Method: " << method << "\n"
345 : << "StoppingCondition: " << condition << "\n"
346 : << "Maximum absolute tolerance: "
347 : << maximum_absolute_tolerance << "\n"
348 : << "Maximum number of iterations: " << maximum_iterations
349 : << "\n"
350 : << "Number of iterations reached: " << iteration_number
351 : << "\n"
352 : << "The last value of f in the root solver is:\n";
353 : for (size_t i = 0; i < Dim; i++) {
354 : error_message << gsl_vector_get(solver->f, i) << "\n";
355 : }
356 : error_message << "The last value of x in the root solver is:\n";
357 : for (size_t i = 0; i < Dim; i++) {
358 : error_message << gsl_vector_get(solver->x, i) << "\n";
359 : }
360 : error_message << "The last value of dx in the root solver is:\n";
361 : for (size_t i = 0; i < Dim; i++) {
362 : error_message << gsl_vector_get(solver->dx, i) << "\n";
363 : }
364 :
365 : if (UNLIKELY(verbosity >= Verbosity::Debug)) {
366 : Parallel::printf("Error: %s\n", gsl_strerror(status));
367 : if (iteration_number >= maximum_iterations) {
368 : Parallel::printf(
369 : "The number of iterations (%d) has reached the maximum number of "
370 : "iterations (%d)\n",
371 : iteration_number, maximum_iterations);
372 : } else {
373 : Parallel::printf(
374 : "The number of iterations (%d) failed to reach the maximum number "
375 : "of iterations (%d)\n",
376 : iteration_number, maximum_iterations);
377 : }
378 : }
379 : throw convergence_error(error_message.str());
380 : }
381 :
382 : return gsl_to_std_array<Dim>(solver->x);
383 : }
384 :
385 : void print_rootfinding_parameters(Method method,
386 : double maximum_absolute_tolerance,
387 : const StoppingCondition& condition);
388 : } // namespace gsl_multiroot_detail
389 :
390 : /// @{
391 : /*!
392 : * \ingroup NumericalAlgorithmsGroup
393 : * \brief A multidimensional root finder supporting Newton and Hybrid
394 : * methods, as well as modified methods based on these.
395 : *
396 : * This root finder accepts function objects with and without a
397 : * callable `jacobian` member function. The call operator both
398 : * accepts and returns a `std::array<double, Dim>`, of the appropriate
399 : * dimension for the domain and range the function the root find is
400 : * being performed on.
401 : *
402 : * If a `jacobian` function is provided, it must accept a
403 : * `std::array<double, Dim>` and return a
404 : * `std::array<std::array<double, Dim>, Dim>` representing the
405 : * derivative of the call operator, with
406 : * \begin{equation}
407 : * \text{jacobian[i][j]} = \frac{\partial f_i}{x_j}.
408 : * \end{equation}
409 : *
410 : * Whether the jacobian is provided determines the details of the
411 : * implementation of the root-finding method that is selected by the
412 : * user using the Method enum. That is, whether the jacobian is
413 : * computed analytically via the `jacobian` member function, or
414 : * whether the jacobian is computed numerically via a finite
415 : * difference approximation.
416 : *
417 : * \note GSL does not provide a finite difference version of its modified
418 : * Newton method, so the unmodified one is used instead when the user
419 : * uses the Method::Newton method.
420 : *
421 : * The user can select one of two possible criteria for convergence,
422 : * StoppingCondition::Residual, where the sum of the absolute values of the
423 : * components of the residual vector f are compared against the value
424 : * provided to `absolute_tolerance`, and
425 : * StoppingCondition::Convergence, where the size of the most recent
426 : * step taken in the root-finding iteration is compared against
427 : * `absolute_tolerance` + `relative_tolerance` * |x_i|, for each component.
428 : * In either case, a `maximum_absolute_tolerance` may be specified if the user
429 : * anticipates that the convergence criterion specified with StoppingCondition
430 : * will be too strict for a few points out of a population of points found with
431 : * a sequence of root finds.
432 : *
433 : * See
434 : * [GSL's documentation for multidimensional
435 : * rootfinding](https://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html)
436 : * for reference.
437 : *
438 : * \param func Function whose root is to be found.
439 : * \param initial_guess Contains initial guess.
440 : * \param maximum_iterations The maximum number of iterations.
441 : * \param verbosity Whether to print diagnostic messages.
442 : * \param maximum_absolute_tolerance Acceptable absolute tolerance when
443 : * root finder doesn't converge.
444 : * You may wish to use this if there
445 : * are only a few "problematic" points where
446 : * it is difficult to do a precise root find.
447 : * \param method The method to use. See the documentation for the Method enum.
448 : * \param condition The convergence condition to use. See the documentation
449 : * for the StoppingCondition enum.
450 : */
451 : template <size_t Dim, typename Function,
452 : Requires<gsl_multiroot_detail::is_jacobian_callable_v<
453 : Function, std::array<double, Dim>>> = nullptr>
454 1 : std::array<double, Dim> gsl_multiroot(
455 : const Function& func, const std::array<double, Dim>& initial_guess,
456 : const StoppingCondition& condition, const size_t maximum_iterations,
457 : const Verbosity verbosity = Verbosity::Silent,
458 : const double maximum_absolute_tolerance = 0.0,
459 : const Method method = Method::Newton) {
460 : gsl_multiroot_function_fdf gsl_func = {
461 : &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_f<Dim, Function>,
462 : &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_df<Dim, Function>,
463 : &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_fdf<Dim,
464 : Function>,
465 : Dim, const_cast<Function*>(&func)}; //NOLINT
466 :
467 : // Set up method for solver:
468 : const gsl_multiroot_fdfsolver_type* solver_type;
469 : if (method == Method::Newton) {
470 : solver_type = gsl_multiroot_fdfsolver_gnewton;
471 : } else if (method == Method::Hybrids) {
472 : solver_type = gsl_multiroot_fdfsolver_hybridsj;
473 : } else if (method == Method::Hybrid) {
474 : solver_type = gsl_multiroot_fdfsolver_hybridj;
475 : } else {
476 : ERROR(
477 : "Invalid method. Has to be one of Newton, Hybrids or "
478 : "Hybrid.");
479 : }
480 : // Print initial parameters
481 : if (UNLIKELY(verbosity >= Verbosity::Verbose)) {
482 : gsl_multiroot_detail::print_rootfinding_parameters(
483 : method, maximum_absolute_tolerance, condition);
484 : }
485 : return gsl_multiroot_detail::gsl_multiroot_impl(
486 : gsl_func, initial_guess, condition, maximum_iterations, verbosity,
487 : maximum_absolute_tolerance, method, solver_type,
488 : &gsl_multiroot_fdfsolver_alloc, &gsl_multiroot_fdfsolver_set,
489 : &gsl_multiroot_fdfsolver_iterate, &gsl_multiroot_fdfsolver_free);
490 : }
491 :
492 : template <size_t Dim, typename Function,
493 : Requires<not gsl_multiroot_detail::is_jacobian_callable_v<
494 : Function, std::array<double, Dim>>> = nullptr>
495 : std::array<double, Dim> gsl_multiroot(
496 : const Function& func, const std::array<double, Dim>& initial_guess,
497 : const StoppingCondition& condition, const size_t maximum_iterations,
498 : const Verbosity verbosity = Verbosity::Silent,
499 : const double maximum_absolute_tolerance = 0.0,
500 : const Method method = Method::Newton) {
501 : gsl_multiroot_function gsl_func = {
502 : &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_f<Dim, Function>,
503 : Dim, const_cast<Function*>(&func)}; // NOLINT
504 :
505 : // Set up method for solver:
506 : const gsl_multiroot_fsolver_type* solver_type;
507 : if (method == Method::Newton) {
508 : solver_type = gsl_multiroot_fsolver_dnewton;
509 : } else if (method == Method::Hybrids) {
510 : solver_type = gsl_multiroot_fsolver_hybrids;
511 : } else if (method == Method::Hybrid) {
512 : solver_type = gsl_multiroot_fsolver_hybrid;
513 : } else {
514 : ERROR(
515 : "Invalid method. Has to be one of Newton, Hybrids or "
516 : "Hybrid.");
517 : }
518 : // Print initial parameters
519 : if (UNLIKELY(verbosity >= Verbosity::Verbose)) {
520 : gsl_multiroot_detail::print_rootfinding_parameters(
521 : method, maximum_absolute_tolerance, condition);
522 : }
523 : return gsl_multiroot_detail::gsl_multiroot_impl(
524 : gsl_func, initial_guess, condition, maximum_iterations, verbosity,
525 : maximum_absolute_tolerance, method, solver_type,
526 : &gsl_multiroot_fsolver_alloc, &gsl_multiroot_fsolver_set,
527 : &gsl_multiroot_fsolver_iterate, &gsl_multiroot_fsolver_free);
528 : }
529 : /// @}
530 : } // namespace RootFinder
|
