From bc4756f768b6338dde09e7ddfb5b9de927f81698 Mon Sep 17 00:00:00 2001 From: tbridel Date: Thu, 1 Dec 2022 22:49:02 +0100 Subject: [PATCH 1/8] stump for rk merson implementation --- CMakeLists.txt | 8 +++++ setup.py | 2 +- src/RK_merson.h | 21 +++++++++++++ src/RK_merson_wrapper.cpp | 64 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 94 insertions(+), 1 deletion(-) create mode 100644 src/RK_merson.h create mode 100644 src/RK_merson_wrapper.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 44b80ad..1bcca88 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -41,12 +41,20 @@ add_library(solve_ivp SHARED set_target_properties(solve_ivp PROPERTIES PREFIX "lib") target_link_libraries(solve_ivp dop853_mod) +# RK 56 Merson +add_library(rk_merson SHARED + ${CMAKE_CURRENT_SOURCE_DIR}/src/RK_merson_wrapper.cpp +) +set_target_properties(rk_merson PROPERTIES PREFIX "lib") + if (SKBUILD) install(TARGETS lsoda DESTINATION numbalsoda) install(TARGETS dop853 DESTINATION numbalsoda) install(TARGETS solve_ivp DESTINATION numbalsoda) + install(TARGETS rk_merson DESTINATION numbalsoda) else() install(TARGETS lsoda DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/) install(TARGETS dop853 DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/) install(TARGETS solve_ivp DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/) + install(TARGETS rk_merson DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/) endif() diff --git a/setup.py b/setup.py index bad23ce..e6b9d63 100644 --- a/setup.py +++ b/setup.py @@ -8,7 +8,7 @@ setup( name="numbalsoda", packages=['numbalsoda'], - version='0.3.5', + version='0.3.6.dev0', license='MIT', install_requires=['numpy','numba'], author = 'Nicholas Wogan', diff --git a/src/RK_merson.h b/src/RK_merson.h new file mode 100644 index 0000000..12859ee --- /dev/null +++ b/src/RK_merson.h @@ -0,0 +1,21 @@ +/** + * @file RK_merson.h + * @author Thibault Bridel-Bertomeu (re-cae.com) + * @brief Merson Runge Kutta 5(6) ODE integrator utilities. + * @version 1.0 + * @date 2022-12-01 + */ + +#ifndef RK_MERSON_H +#define RK_MERSON_H + +#include + +#define RK_MERSON_NSTEPS 5 + +extern const std::array rk_merson_nodes; +extern const std::array rk_merson_weights; +extern const std::array rk_merson_starweights; +extern const std::array, RK_MERSON_NSTEPS> rk_merson_matrix; + +#endif diff --git a/src/RK_merson_wrapper.cpp b/src/RK_merson_wrapper.cpp new file mode 100644 index 0000000..b18d154 --- /dev/null +++ b/src/RK_merson_wrapper.cpp @@ -0,0 +1,64 @@ +/** + * @file RK_merson_wrapper.cpp + * @author Thibault Bridel-Bertomeu (re-cae.com) + * @brief Merson Runge Kutta 5(6) ODE integrator. + * @version 1.0 + * @date 2022-12-01 + */ + +#include "RK_merson.h" + +#include +#include + +const std::array rk_merson_nodes = {0.0, (1.0/3.0), (1.0/3.0), 0.5, 1.0}; +const std::array rk_merson_weights = {(1.0/6.0), 0.0, 0.0, 2.0*(1.0/3.0), (1.0/6.0)}; +const std::array rk_merson_starweights = {0.5, 0.0, -3.0*0.5, 2.0, 0.0}; +const std::array, RK_MERSON_NSTEPS> rk_merson_matrix = { + 0.0, 0.0, 0.0, 0.0, 0.0, + (1.0/3.0), 0.0, 0.0, 0.0, 0.0, + (1.0/6.0), (1.0/6.0), 0.0, 0.0, 0.0, + (1.0/8.0), 0.0, 3.0*(1.0/8.0), 0.0, 0.0, + 0.5, 0.0, -3.0*0.5, 2.0, 0.0, +}; + +extern "C" +{ + +/** + * @brief ODE solver based on Merson's Runge-Kutta 5(6) integrator. + * @details This method implements Merson's Runge-Kutta 5(6) method with adaptive time-stepping + * to solve an ODE. + * + * @param rhs Pointer to the method to compute the right hand side of the ODE. + * @param neq Number of equations. + * @param u0 Initial solution/guess. + * @param data Extra data to be passed down to the RHS. + * @param nt Size of the teval vector - see below. + * @param teval Vector of times whereat to solve the ODE (optional). + * @param usol Vector holding the solution at a given instant. + * @param rtol Relative tolerance for the integration. + * @param atol Absolute tolerance for the integration. + * @param mxstep Maximum number of time steps taken. + * @param success Boolean indicating whether integration was successful or errored. + */ +void rk_merson_wrapper( + void (*rhs)(double t, double *u, double *du, void* data), + int neq, + double* u0, + void* data, + int nt, + double* teval, + double* usol, + double rtol, + double atol, + int mxstep, + int* success +) +{ + + + +} + +} From c48c45d3cb082d1e64d5b0ddf69272198babedff Mon Sep 17 00:00:00 2001 From: tbridel Date: Sat, 3 Dec 2022 15:19:39 +0100 Subject: [PATCH 2/8] Work in progess writing the cpp rk45 --- CMakeLists.txt | 15 +- src/RK45.cpp | 148 ++++++++++++++++++ src/RK45.h | 307 ++++++++++++++++++++++++++++++++++++++ src/RK45_wrapper.cpp | 67 +++++++++ src/RK_merson.h | 21 --- src/RK_merson_wrapper.cpp | 64 -------- 6 files changed, 530 insertions(+), 92 deletions(-) create mode 100644 src/RK45.cpp create mode 100644 src/RK45.h create mode 100644 src/RK45_wrapper.cpp delete mode 100644 src/RK_merson.h delete mode 100644 src/RK_merson_wrapper.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 1bcca88..41c9655 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,7 +11,7 @@ endif() set(CMAKE_POSITION_INDEPENDENT_CODE ON) set(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS ON) -set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD 17) # lsoda add_library(lsoda SHARED @@ -41,20 +41,21 @@ add_library(solve_ivp SHARED set_target_properties(solve_ivp PROPERTIES PREFIX "lib") target_link_libraries(solve_ivp dop853_mod) -# RK 56 Merson -add_library(rk_merson SHARED - ${CMAKE_CURRENT_SOURCE_DIR}/src/RK_merson_wrapper.cpp +# RKF 45 +add_library(rk45 SHARED + ${CMAKE_CURRENT_SOURCE_DIR}/src/RK45.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/RK45_wrapper.cpp ) -set_target_properties(rk_merson PROPERTIES PREFIX "lib") +set_target_properties(rk45 PROPERTIES PREFIX "lib") if (SKBUILD) install(TARGETS lsoda DESTINATION numbalsoda) install(TARGETS dop853 DESTINATION numbalsoda) install(TARGETS solve_ivp DESTINATION numbalsoda) - install(TARGETS rk_merson DESTINATION numbalsoda) + install(TARGETS rk45 DESTINATION numbalsoda) else() install(TARGETS lsoda DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/) install(TARGETS dop853 DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/) install(TARGETS solve_ivp DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/) - install(TARGETS rk_merson DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/) + install(TARGETS rk45 DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/) endif() diff --git a/src/RK45.cpp b/src/RK45.cpp new file mode 100644 index 0000000..5f26a01 --- /dev/null +++ b/src/RK45.cpp @@ -0,0 +1,148 @@ +/** + * @file RK45.cpp + * @author Thibault Bridel-Bertomeu (re-cae.com) + * @brief Runge-Kutta-Felhberg method of order 5(4); class definition; + * @version 1.0 + * @date 2022-12-01 + */ + +#include +#include +#include +#include + +#include "RK45.h" + +/** + * @brief Construct a new RK45::RK45 object + * + * @param fun Callable, right-hand side of the system. + * @param t0 Initial time. + * @param y0 Initial state. + * @param t_bound Boundary time, the integration won't continue beyond it. + * @param max_step Maximum allowed step size; default is infinity, i.e. the step size + * is not bounded and determined solely by the solver. + * @param rtol Relative tolerance. + * @param atol Absolute tolerance. + * @param first_step Initial step size; default is "None" which means that the algorithm + * should choose. + */ +RK45::RK45(RK45_ODE_SYSTEM_TYPE rhs_handle, double t0, std::vector &y0, double tf, + int largest_step, + double _rtol, + double _atol, + std::optional first_step) +{ + _m_status = "running"; + m_fun = rhs_handle; + + _m_t_old = std::nullopt; + m_t = t0; + m_t_bound = tf; + m_direction = (m_t_bound != t0) ? sgn(m_t_bound - t0) : 1; + + if (_rtol < 100 * std::numeric_limits::epsilon()) { + std::cout << "Warning: '_rtol' is too small. Setting to max(rtol, 100*EPS)." << std::endl; + m_rtol = std::max(_rtol, 100 * std::numeric_limits::epsilon()); + } + if (_atol < 0) { + throw std::invalid_argument("'_atol' must be positive."); + } + m_atol = _atol; + + m_y = y0; + _m_y_old = std::nullopt; + m_n_eqns = m_y.size(); + m_fun(m_t, m_y.data(), m_f.data(), NULL); + + + if (largest_step <= int(0)) + throw std::invalid_argument("'largest_step' must be positive."); + m_max_step = largest_step; + if (first_step.has_value()) { + if (first_step.value() <= 0) + throw std::invalid_argument("'first_step' must be positive."); + if (first_step.value() > std::abs(tf - t0)) + throw std::invalid_argument("'first_step' exceeds bounds."); + m_h_abs = first_step.value(); + } else { + m_h_abs = select_initial_step(m_fun, m_t, m_y, m_f, m_direction, _m_error_estimation_order, m_rtol, m_atol); + } + _m_h_previous = std::nullopt; + + _m_K.reserve(_m_n_stages + 1); + for (auto i = 0; i < _m_n_stages + 1; i++) + _m_K[i].reserve(m_n_eqns); +} + +/** + * @brief Estimate the current integration error. + * + * @param h Current step size. + * @return std::vector Current integration error per equation. + */ +std::vector RK45::_estimate_error(double h) +{ + std::vector dot_K_E(m_n_eqns); + dgemv("T", _m_n_stages + 1, m_n_eqns, 1.0, _m_K, _m_n_stages + 1, _m_E, 1, 0.0, dot_K_E, 1); + + for (auto i = 0; i < dot_K_E.size(); i++) { + dot_K_E[i] *= h; + } + + return dot_K_E; +} + +/** + * @brief Norm of the current integration error. + * + * @param h Current step size. + * @param scale Scaling factor for the current integration error per equation. + * @return double Norm of the current integration error, each equation rescaled by the + * ad hoc factor based on user specified tolerances. + */ +double RK45::_estimate_error_norm(double h, std::vector scale) +{ + std::vector itg_error = _estimate_error(h); + double norm = double(0); + for (auto i = 0; i < scale.size(); i++) { + norm += (itg_error[i] / scale[i]) * (itg_error[i] / scale[i]); + } + + return std::sqrt(norm); +} + +/** + * @brief Computations for one integration step. + * + * @return bool Whether the step was successful. + */ +bool RK45::step() +{ + double min_step = double(10) * std::numeric_limits::epsilon(); + double h_abs; + if (m_h_abs > m_max_step) { + h_abs = m_max_step; + } else if (m_h_abs < min_step) { + h_abs = min_step; + } else { + h_abs = m_h_abs; + } + + bool step_accepted = false; + bool step_rejected = false; + + while (not step_accepted) { + if (h_abs < min_step) { + return false; + } + + double h = h_abs * m_direction; + double t_new = m_t + dt; + if (m_direction * (t_new - m_t_bound) > double(0)) { + t_new = m_t_bound; + } + h = t_new - m_t; + h_abs = std::abs(h); + } +} diff --git a/src/RK45.h b/src/RK45.h new file mode 100644 index 0000000..d1a6e26 --- /dev/null +++ b/src/RK45.h @@ -0,0 +1,307 @@ +/** + * @file RK45.h + * @author Thibault Bridel-Bertomeu (re-cae.com) + * @brief Runge-Kutta-Felhberg method of order 5(4); class declaration. + * @version 1.0 + * @date 2022-12-01 + */ + +#pragma once + +#include +#include +#include +#include +#include +#include + +/** + * @brief Multiply steps computed from asymptotic behavior of errors by this. + */ +#define SAFETY 0.9 + +/** + * @brief Minimum allowed decrease in a step size. + */ +#define MIN_FACTOR 0.2 + +/** + * @brief Maximum allowed increased in a step size. + */ +#define MAX_FACTOR 10.0 + +/** + * @brief Type definition of RK45 ode system. + */ +typedef void (*RK45_ODE_SYSTEM_TYPE)(double t, double *y, double *dydt, void*); + +/** + * @brief Runge-Kutta-Fehlberg 5(4) explicit integration method. + * @details This uses the Dormand-Prince pair of formulas [1]. The error is controlled + * assuming accuracy of the fourth-order method accuracy, but steps are taken + * using the fifth-order accurate formula (local extrapolation is done). + * A quartic interpolation polynomial is used for the dense output [2]. + * + * References: + * 1. J. R. Dormand, P. J. Prince, "A family of embedded Runge-Kutta + * formulae", Journal of Computational and Applied Mathematics, Vol. 6, + * No. 1, pp. 19-26, 1980. + * 2. L. W. Shampine, "Some Practical Runge-Kutta Formulas", Mathematics + * of Computation,, Vol. 46, No. 173, pp. 135-150, 1986. + */ +class RK45 { + +public: + RK45(RK45_ODE_SYSTEM_TYPE fun, double t0, std::vector &y0, double tf, + int largest_step=std::numeric_limits::max(), + double _rtol=1e-3, + double _atol=1e-6, + std::optional first_step=std::nullopt); + ~RK45(); + + RK45_ODE_SYSTEM_TYPE m_fun; + int m_n_eqns; + int m_max_step; + double m_t_bound; + double m_rtol; + double m_atol; + + std::vector m_y; + std::vector m_f; + double m_h_abs; + double m_t; + int m_direction; + + bool step(); + +private: + static int const _m_order = 5; + static int const _m_error_estimation_order = 4; + static constexpr double const _m_error_exponent = -1.0 / (_m_error_estimation_order + 1.0); + static int const _m_n_stages = 6; + + std::vector const _m_C {0.0, 1.0/5.0, 3.0/10.0, 4.0/5.0, 8.0/9.0, 1.0}; // size _n_stages + std::vector> const _m_A { + { 0.0, 0.0, 0.0, 0.0, 0.0 }, + { 1.0/5.0, 0.0, 0.0, 0.0, 0.0 }, + { 3.0/40.0, 9.0/40.0, 0.0, 0.0, 0.0 }, + { 44.0/45.0, -56.0/15.0, 32.0/9.0, 0.0, 0.0 }, + { 19372.0/6561.0, -25360.0/2187.0, 64448.0/6561.0, -212.0/729.0, 0.0 }, + { 9017.0/3168.0, -355.0/33.0, 46732.0/5247.0, 49.0/176.0, -5103.0/1865.0 }, + }; // size _n_stages x _n_stages + std::vector const _m_B {35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0}; // size _n_stages + std::vector const _m_E {-71.0/57600.0, 0.0, 71.0/16695.0, -71.0/1920.0, 17253.0/339200.0, -22.0/525.0, 1.0/40.0}; // size _n_stages + 1 + + std::optional _m_h_previous; + std::optional _m_t_old; + std::optional> _m_y_old; + std::vector> _m_K; + + std::string _m_status; + + std::vector _estimate_error(double h); + double _estimate_error_norm(double h, std::vector scale); +}; + + +/** + * @brief Implements signum (-1, 0, 1). + * + * @tparam T Type of the number whose sign to evaluate. + * @param val Number whose sign to evaluate. + * @return int -1 if the number is negative, 1 if the number is positive, 0 elsehow. + */ +template int sgn(T val) +{ + return (T(0) < val) - (val < T(0)); +} + + +/** + * @brief Empirically select a good initial step. + * @details The algorithm is described in [1]. + * + * References: + * - E. Hairer, S. P. Norsett, G. Wanner, "Solving Ordinary Differential + * Equations I: Nonstiff problems", Sec. II.4. + * + * @param fun Callable, right-hand side of the system. + * @param t0 Initial value of the independant variable. + * @param y0 Initial value of the dependant variable. + * @param f0 Initial value of the derivatives. + * @param direction Integration direction. + * @param order Error estimation order, i.e. the error controlled by the + * algorithm is proportional to 'step_size ** (order + 1)'. + * @param rtol Desired relative tolerance. + * @param atol Desired absolute tolerance. + * @return double Absolute value of the suggested initial step. + */ +double select_initial_step( + RK45_ODE_SYSTEM_TYPE fun, + double t0, + std::vector &y0, + std::vector &f0, + int direction, + int order, + double rtol, + double atol +) +{ + int n_eqns = y0.size(); + if (n_eqns == 0) + return std::numeric_limits::max(); + + std::vector scale(n_eqns); + double d0 = 0.0, d1 = 0.0; + for (auto i = 0; i < n_eqns; i++) { + scale[i] = atol + std::abs(y0[i]) * rtol; + d0 += (y0[i] / scale[i]) * (y0[i] / scale[i]); + d1 += (f0[i] / scale[i]) * (f0[i] / scale[i]); + } + d0 = std::sqrt(d0); + d1 = std::sqrt(d1); + + double h0; + if ((d0 < 1e-5) || (d1 < 1e-5)) { + h0 = 1e-6; + } else { + h0 = 0.01 * d0 / d1; + } + + std::vector y1(n_eqns), f1(n_eqns); + for (auto i = 0; i < n_eqns; i++) { + y1[i] = y0[i] + h0 * direction * f0[i]; + } + fun(t0 + h0 * direction, y1.data(), f1.data(), NULL); + double d2 = 0.0; + for (auto i = 0; i < n_eqns; i++) { + d2 += (f1[i] - f0[i]) / scale[i]; + } + d2 = std::sqrt(d2) / h0; + + double h1; + if ((d1 <= 1e-15) && (d2 <= 1e-15)) { + h1 = std::max(1e-06, h0 * 1e-03); + } else { + h1 = std::pow((0.01 / std::max(d1, d2)), 1.0 / (order + 1)); + } + + return std::min(100 * h0, h1); +} + + +/** + * @brief Performs y = alpha*A*x + beta*y or y = alpha*A^T*x + beta*y + * + * @param trans Specifies the operation to be performed. If equals "N" then + * y = alpha * A * x + beta * y, if equals "T" then y = alpha * A^T * x + beta * y. + * @param m Number of rows of the matrix A. + * @param n Number of columns of the matrix A. + * @param alpha Scalar alpha. + * @param a m x n matrix of coefficients A. + * @param lda First dimension of A. + * @param x Vector X. + * @param incx Increment for the elements of X. + * @param beta Scalar beta. + * @param y Vector Y; will be overwritten by the updated vector Y. + * @param incy Increment for the elements of Y. + */ +void dgemv(const std::string trans, const size_t m, const size_t n, + const double alpha, const std::vector> &a, + const int lda, const std::vector &x, const int incx, + const double beta, std::vector &y, const int incy) +{ + + int lenx, leny; + if (trans == "N") { + lenx = n; + leny = m; + } else { + lenx = m; + leny = n; + } + + int kx, ky; + if (incx > 0) { + kx = 0; + } else { + kx = - (lenx - 1) * incx; + } + if (incy > 0) { + ky = 0; + } else { + ky = - (leny - 1) * incy; + } + + if (std::abs(beta - double(1)) > 0) { + if (incy == 1) { + if (std::abs(beta) <= double(0)) { + std::fill(y.begin(), y.end(), double(0)); + } else { + for (auto i = 0; i < leny; i++) + y[i] *= beta; + } + } else { + int iy = ky; + if (std::abs(beta) <= double(0)) { + for (auto i = 0; i < leny; i++) { + y[iy] = double(0); + iy += incy; + } + } else { + for (auto i = 0; i < leny; i++) { + y[iy] *= beta; + iy += incy; + } + } + } + } + + // Form y = alpha * A * x + beta * y + if (trans == "N") { + int jx = kx; + if (incy == 1) { + for (auto j = 0; j < n; j++) { + double temp = alpha * x[jx]; + for (auto i = 0; i < m; i++) { + y[i] = y[i] + temp * a[i][j]; + } + jx += incx; + } + } else { + for (auto j = 0; j < n; j++) { + double temp = alpha * x[jx]; + int iy = ky; + for (auto i = 0; i < m; i++) { + y[iy] = y[iy] + temp * a[i][j]; + iy += incy; + } + jx += incx; + } + } + // Form y = alpha * A^T * x + beta * y + } else { + int jy = ky; + if (incx == 1) { + for (auto j = 1; j < n; j++) { + double temp = double(0); + for(auto i = 0; i < m; i++) { + temp += a[i][j] * x[i]; + } + y[jy] += alpha * temp; + jy += incy; + } + } else { + for (auto j = 1; j < n; j++) { + double temp = double(0); + int ix = kx; + for (auto i = 0; i < m; i++) { + temp += a[i][j] * x[ix]; + ix += incx; + } + y[jy] += alpha * temp; + jy += incy; + } + } + } +} diff --git a/src/RK45_wrapper.cpp b/src/RK45_wrapper.cpp new file mode 100644 index 0000000..fc825cb --- /dev/null +++ b/src/RK45_wrapper.cpp @@ -0,0 +1,67 @@ +/** + * @file RK45_wrapper.cpp + * @author Thibault Bridel-Bertomeu (re-cae.com) + * @brief Runge-Kutta-Felhberg method of order 5(4) with adaptive timestepping. + * @version 1.0 + * @date 2022-12-01 + */ + +#include +#include +#include +#include + +struct RK45 { + static const int n_stages = 6; + const std::array rk_merson_nodes = {0.0, (1.0/3.0), (1.0/3.0), 0.5, 1.0}; + const std::array rk_merson_weights = {(1.0/6.0), 0.0, 0.0, 2.0*(1.0/3.0), (1.0/6.0)}; + const std::array rk_merson_starweights = {0.5, 0.0, -3.0*0.5, 2.0, 0.0}; + const std::array, n_stages> rk_merson_matrix = { + 0.0, 0.0, 0.0, 0.0, 0.0, + (1.0/3.0), 0.0, 0.0, 0.0, 0.0, + (1.0/6.0), (1.0/6.0), 0.0, 0.0, 0.0, + (1.0/8.0), 0.0, (3.0/8.0), 0.0, 0.0, + 0.5, 0.0, -1.5, 2.0, 0.0 + }; +}; + +extern "C" +{ + +/** + * @brief ODE solver based on Merson's Runge-Kutta 5(6) integrator. + * @details This method implements Merson's Runge-Kutta 5(6) method with adaptive time-stepping + * to solve an ODE. + * + * @param rhs Pointer to the method to compute the right hand side of the ODE. + * @param neq Number of equations and unknowns + * @param data Extra data to be passed down to the RHS. + * @param t0 Initial time. + * @param u0 Initial solution/guess. + * @param tf Final time + * @param itf Maximum number of iterations + * @param rtol Relative tolerance for the integration. + * @param atol Absolute tolerance for the integration. + * @param dt0 Initial time step, optional + * @param usol Vector holding the solution at a given instant. + * @param success Boolean indicating whether integration was successful or errored. + */ +void rk_merson_wrapper( + void (*rhs)(double t, double *u, double *du, void* data), + int neq, + void* data, + double t0, + double* u0, + double tf, + int itf, + double rtol, + double atol, + double dt0, + double* usol, + int* success +) +{ + +} // end void rk_merson_wrapper(...) + +} // end extern "C" diff --git a/src/RK_merson.h b/src/RK_merson.h deleted file mode 100644 index 12859ee..0000000 --- a/src/RK_merson.h +++ /dev/null @@ -1,21 +0,0 @@ -/** - * @file RK_merson.h - * @author Thibault Bridel-Bertomeu (re-cae.com) - * @brief Merson Runge Kutta 5(6) ODE integrator utilities. - * @version 1.0 - * @date 2022-12-01 - */ - -#ifndef RK_MERSON_H -#define RK_MERSON_H - -#include - -#define RK_MERSON_NSTEPS 5 - -extern const std::array rk_merson_nodes; -extern const std::array rk_merson_weights; -extern const std::array rk_merson_starweights; -extern const std::array, RK_MERSON_NSTEPS> rk_merson_matrix; - -#endif diff --git a/src/RK_merson_wrapper.cpp b/src/RK_merson_wrapper.cpp deleted file mode 100644 index b18d154..0000000 --- a/src/RK_merson_wrapper.cpp +++ /dev/null @@ -1,64 +0,0 @@ -/** - * @file RK_merson_wrapper.cpp - * @author Thibault Bridel-Bertomeu (re-cae.com) - * @brief Merson Runge Kutta 5(6) ODE integrator. - * @version 1.0 - * @date 2022-12-01 - */ - -#include "RK_merson.h" - -#include -#include - -const std::array rk_merson_nodes = {0.0, (1.0/3.0), (1.0/3.0), 0.5, 1.0}; -const std::array rk_merson_weights = {(1.0/6.0), 0.0, 0.0, 2.0*(1.0/3.0), (1.0/6.0)}; -const std::array rk_merson_starweights = {0.5, 0.0, -3.0*0.5, 2.0, 0.0}; -const std::array, RK_MERSON_NSTEPS> rk_merson_matrix = { - 0.0, 0.0, 0.0, 0.0, 0.0, - (1.0/3.0), 0.0, 0.0, 0.0, 0.0, - (1.0/6.0), (1.0/6.0), 0.0, 0.0, 0.0, - (1.0/8.0), 0.0, 3.0*(1.0/8.0), 0.0, 0.0, - 0.5, 0.0, -3.0*0.5, 2.0, 0.0, -}; - -extern "C" -{ - -/** - * @brief ODE solver based on Merson's Runge-Kutta 5(6) integrator. - * @details This method implements Merson's Runge-Kutta 5(6) method with adaptive time-stepping - * to solve an ODE. - * - * @param rhs Pointer to the method to compute the right hand side of the ODE. - * @param neq Number of equations. - * @param u0 Initial solution/guess. - * @param data Extra data to be passed down to the RHS. - * @param nt Size of the teval vector - see below. - * @param teval Vector of times whereat to solve the ODE (optional). - * @param usol Vector holding the solution at a given instant. - * @param rtol Relative tolerance for the integration. - * @param atol Absolute tolerance for the integration. - * @param mxstep Maximum number of time steps taken. - * @param success Boolean indicating whether integration was successful or errored. - */ -void rk_merson_wrapper( - void (*rhs)(double t, double *u, double *du, void* data), - int neq, - double* u0, - void* data, - int nt, - double* teval, - double* usol, - double rtol, - double atol, - int mxstep, - int* success -) -{ - - - -} - -} From dc5fbec661db39d81eb129f67cea12f6a81a1e7c Mon Sep 17 00:00:00 2001 From: tbridel Date: Mon, 5 Dec 2022 14:59:59 +0100 Subject: [PATCH 3/8] First functional version of RK45 --- numbalsoda/__init__.py | 3 +- numbalsoda/driver_rk45.py | 52 ++++++ src/RK45.cpp | 333 ++++++++++++++++++++++++++++++++++++-- src/RK45.h | 202 +---------------------- src/RK45_wrapper.cpp | 89 ++++++---- tests/test_rk45.py | 40 +++++ 6 files changed, 471 insertions(+), 248 deletions(-) create mode 100644 numbalsoda/driver_rk45.py create mode 100644 tests/test_rk45.py diff --git a/numbalsoda/__init__.py b/numbalsoda/__init__.py index 1a3fe85..5ebadf6 100644 --- a/numbalsoda/__init__.py +++ b/numbalsoda/__init__.py @@ -1,3 +1,4 @@ from .driver import lsoda_sig, lsoda, address_as_void_pointer from .driver_dop853 import dop853 -from .driver_solve_ivp import solve_ivp, ODEResult \ No newline at end of file +from .driver_solve_ivp import solve_ivp, ODEResult +from .driver_rk45 import rk45_sig, rk45 diff --git a/numbalsoda/driver_rk45.py b/numbalsoda/driver_rk45.py new file mode 100644 index 0000000..3f52963 --- /dev/null +++ b/numbalsoda/driver_rk45.py @@ -0,0 +1,52 @@ +import ctypes as ct +import numba as nb +from numba import njit, types +import numpy as np +import os +import platform + +rk45_sig = types.void(types.double, + types.CPointer(types.double), + types.CPointer(types.double), + types.CPointer(types.double)) + +rootdir = os.path.dirname(os.path.realpath(__file__))+'/' + +if platform.uname()[0] == "Windows": + name = "librk45.dll" +elif platform.uname()[0] == "Linux": + name = "librk45.so" +else: + name = "librk45.dylib" +librk45 = ct.CDLL(rootdir+name) +rk45_wrapper = librk45.rk45_wrapper +rk45_wrapper.argtypes = [ct.c_void_p, ct.c_int, ct.c_void_p, ct.c_void_p,\ + ct.c_double, ct.c_double, ct.c_double, ct.c_int,\ + ct.c_void_p, ct.c_double, ct.c_double, ct.c_double, ct.c_void_p,\ + ct.c_void_p, ct.c_void_p] +rk45_wrapper.restype = None + +@njit +def rk45(funcptr, u0, dt0, t0, tf, itf, data = np.array([0.0], np.float64), rtol = 1.0e-3, atol = 1.0e-6, mxstep = 10000.0): + neq = len(u0) + usol = np.empty((itf+1, neq),dtype=np.float64) + success = np.array([999], np.int32) + actual_final_time = np.array([-1.0], dtype=np.float64) + actual_final_iteration = np.array([-1], dtype=np.int32) + rk45_wrapper(funcptr, neq, u0.ctypes.data, data.ctypes.data, dt0, t0, tf, itf, + usol.ctypes.data, rtol, atol, mxstep, success.ctypes.data, + actual_final_time.ctypes.data, actual_final_iteration.ctypes.data) + success_ = True + if success != 1: + success_ = False + return usol[:actual_final_iteration[0],:], actual_final_time[0], success_ + +@nb.extending.intrinsic +def address_as_void_pointer(typingctx, src): + """ returns a void pointer from a given memory address """ + from numba import types + from numba.core import cgutils + sig = types.voidptr(src) + def codegen(cgctx, builder, sig, args): + return builder.inttoptr(args[0], cgutils.voidptr_t) + return sig, codegen diff --git a/src/RK45.cpp b/src/RK45.cpp index 5f26a01..00bba70 100644 --- a/src/RK45.cpp +++ b/src/RK45.cpp @@ -14,7 +14,261 @@ #include "RK45.h" /** - * @brief Construct a new RK45::RK45 object + * @brief Empirically select a good initial step. + * @details The algorithm is described in [1]. + * + * References: + * - E. Hairer, S. P. Norsett, G. Wanner, "Solving Ordinary Differential + * Equations I: Nonstiff problems", Sec. II.4. + * + * @param fun Callable, right-hand side of the system. + * @param t0 Initial value of the independant variable. + * @param y0 Initial value of the dependant variable. + * @param f0 Initial value of the derivatives. + * @param direction Integration direction. + * @param order Error estimation order, i.e. the error controlled by the + * algorithm is proportional to 'step_size ** (order + 1)'. + * @param rtol Desired relative tolerance. + * @param atol Desired absolute tolerance. + * @return double Absolute value of the suggested initial step. + */ +double select_initial_step( + RK45_ODE_SYSTEM_TYPE fun, + double t0, + std::vector &y0, + std::vector &f0, + int direction, + int order, + double rtol, + double atol +) +{ + int n_eqns = y0.size(); + if (n_eqns == 0) + return std::numeric_limits::max(); + + std::vector scale(n_eqns); + double d0 = 0.0, d1 = 0.0; + for (auto i = 0; i < n_eqns; i++) { + scale[i] = atol + std::abs(y0[i]) * rtol; + d0 += (y0[i] / scale[i]) * (y0[i] / scale[i]); + d1 += (f0[i] / scale[i]) * (f0[i] / scale[i]); + } + d0 = std::sqrt(d0); + d1 = std::sqrt(d1); + + double h0; + if ((d0 < 1e-5) || (d1 < 1e-5)) { + h0 = 1e-6; + } else { + h0 = 0.01 * d0 / d1; + } + + std::vector y1(n_eqns), f1(n_eqns); + for (auto i = 0; i < n_eqns; i++) { + y1[i] = y0[i] + h0 * direction * f0[i]; + } + fun(t0 + h0 * direction, y1.data(), f1.data(), NULL); + double d2 = 0.0; + for (auto i = 0; i < n_eqns; i++) { + d2 += (f1[i] - f0[i]) / scale[i]; + } + d2 = std::sqrt(d2) / h0; + + double h1; + if ((d1 <= 1e-15) && (d2 <= 1e-15)) { + h1 = std::max(1e-06, h0 * 1e-03); + } else { + h1 = std::pow((0.01 / std::max(d1, d2)), 1.0 / (order + 1)); + } + + return std::min(100 * h0, h1); +} + +/** + * @brief Perform a single Runge-Kutta step. + * @details This method is integrator-agnostic and is valid for all explicit Runge-Kutta + * built with a Butcher tableau. + * + * @param fun Callable, right-hand side of the system. + * @param t Current value of the independent variable. + * @param y Current value of the dependent variable. + * @param f Current value of the derivative. + * @param h Step to use. + * @param A Coefficients for combining previous RK stages to compute the next stage. + * @param B Coefficients for combining RK stages for computing the final prediction. + * @param C Coefficients for incrementing time for consecutive RK stages. + * @param K Storage array for putting RK stages here. Stages are stored in rows. The last + * row is a linear combination of the previous rows with coefficients. + * @return std::array, 2> Solution at t + h and derivative thereof. + */ +std::array, 2> rk_step( + RK45_ODE_SYSTEM_TYPE fun, + double const t, + std::vector const& y, + std::vector const& f, + double const h, + std::vector> const& A, + std::vector const& B, + std::vector const& C, + std::vector>& K +) +{ + int const n_eqns = y.size(); + int const n_stages = K.size() - 1; + + for (auto i = 0; i < n_eqns; i++) { + K[0][i] = f[i]; + } + for (auto s = 1; s < n_stages; s++) { + std::vector dy(n_eqns, double(0)); + for (auto iss = 0; iss < s; iss++) { + for (auto i = 0; i < n_eqns; i++) { + dy[i] += K[iss][i] * A[s][iss] * h; + } + } + for (auto i = 0; i < n_eqns; i++) { + dy[i] += y[i]; + } + fun(t + C[s] * h, dy.data(), K[s].data(), NULL); + } + + std::vector y_new(n_eqns, double(0)), f_new(n_eqns); + for (auto s = 0; s < n_stages; s++) { + for (auto i = 0; i < n_eqns; i++) { + y_new[i] += K[s][i] * B[s] * h; + } + } + for (auto i = 0; i < n_eqns; i++) { + y_new[i] += y[i]; + } + fun(t + h, y_new.data(), f_new.data(), NULL); + for (auto i = 0; i < n_eqns; i++) { + K[n_stages][i] = f_new[i]; + } + + std::array, 2> result{y_new, f_new}; + return result; +} + +/** + * @brief Performs y = alpha*A*x + beta*y or y = alpha*A^T*x + beta*y + * + * @param trans Specifies the operation to be performed. If equals "N" then + * y = alpha * A * x + beta * y, if equals "T" then y = alpha * A^T * x + beta * y. + * @param m Number of rows of the matrix A. + * @param n Number of columns of the matrix A. + * @param alpha Scalar alpha. + * @param a m x n matrix of coefficients A. + * @param lda First dimension of A. + * @param x Vector X. + * @param incx Increment for the elements of X. + * @param beta Scalar beta. + * @param y Vector Y; will be overwritten by the updated vector Y. + * @param incy Increment for the elements of Y. + */ +void dgemv(const std::string trans, const size_t m, const size_t n, + const double alpha, const std::vector> &a, + const int lda, const std::vector &x, const int incx, + const double beta, std::vector &y, const int incy) +{ + + int lenx, leny; + if (trans == "N") { + lenx = n; + leny = m; + } else { + lenx = m; + leny = n; + } + + int kx, ky; + if (incx > 0) { + kx = 0; + } else { + kx = - (lenx - 1) * incx; + } + if (incy > 0) { + ky = 0; + } else { + ky = - (leny - 1) * incy; + } + + if (std::abs(beta - double(1)) > 0) { + if (incy == 1) { + if (std::abs(beta) <= double(0)) { + std::fill(y.begin(), y.end(), double(0)); + } else { + for (auto i = 0; i < leny; i++) + y[i] *= beta; + } + } else { + int iy = ky; + if (std::abs(beta) <= double(0)) { + for (auto i = 0; i < leny; i++) { + y[iy] = double(0); + iy += incy; + } + } else { + for (auto i = 0; i < leny; i++) { + y[iy] *= beta; + iy += incy; + } + } + } + } + + // Form y = alpha * A * x + beta * y + if (trans == "N") { + int jx = kx; + if (incy == 1) { + for (auto j = 0; j < n; j++) { + double temp = alpha * x[jx]; + for (auto i = 0; i < m; i++) { + y[i] = y[i] + temp * a[i][j]; + } + jx += incx; + } + } else { + for (auto j = 0; j < n; j++) { + double temp = alpha * x[jx]; + int iy = ky; + for (auto i = 0; i < m; i++) { + y[iy] = y[iy] + temp * a[i][j]; + iy += incy; + } + jx += incx; + } + } + // Form y = alpha * A^T * x + beta * y + } else { + int jy = ky; + if (incx == 1) { + for (auto j = 1; j < n; j++) { + double temp = double(0); + for(auto i = 0; i < m; i++) { + temp += a[i][j] * x[i]; + } + y[jy] += alpha * temp; + jy += incy; + } + } else { + for (auto j = 1; j < n; j++) { + double temp = double(0); + int ix = kx; + for (auto i = 0; i < m; i++) { + temp += a[i][j] * x[ix]; + ix += incx; + } + y[jy] += alpha * temp; + jy += incy; + } + } + } +} + +/** + * @brief Construct a new RK45::RK45 object. * * @param fun Callable, right-hand side of the system. * @param t0 Initial time. @@ -28,15 +282,14 @@ * should choose. */ RK45::RK45(RK45_ODE_SYSTEM_TYPE rhs_handle, double t0, std::vector &y0, double tf, - int largest_step, + double largest_step, double _rtol, double _atol, - std::optional first_step) + double first_step) { _m_status = "running"; m_fun = rhs_handle; - _m_t_old = std::nullopt; m_t = t0; m_t_bound = tf; m_direction = (m_t_bound != t0) ? sgn(m_t_bound - t0) : 1; @@ -51,30 +304,32 @@ RK45::RK45(RK45_ODE_SYSTEM_TYPE rhs_handle, double t0, std::vector &y0, m_atol = _atol; m_y = y0; - _m_y_old = std::nullopt; m_n_eqns = m_y.size(); + m_f.resize(m_n_eqns); m_fun(m_t, m_y.data(), m_f.data(), NULL); - - if (largest_step <= int(0)) + if (largest_step <= double(0)) throw std::invalid_argument("'largest_step' must be positive."); m_max_step = largest_step; - if (first_step.has_value()) { - if (first_step.value() <= 0) - throw std::invalid_argument("'first_step' must be positive."); - if (first_step.value() > std::abs(tf - t0)) + if (first_step > 0) { + if (first_step > std::abs(tf - t0)) throw std::invalid_argument("'first_step' exceeds bounds."); - m_h_abs = first_step.value(); + m_h_abs = first_step; } else { m_h_abs = select_initial_step(m_fun, m_t, m_y, m_f, m_direction, _m_error_estimation_order, m_rtol, m_atol); } - _m_h_previous = std::nullopt; + _m_h_previous = -1.0; - _m_K.reserve(_m_n_stages + 1); + _m_K.resize(_m_n_stages + 1); for (auto i = 0; i < _m_n_stages + 1; i++) - _m_K[i].reserve(m_n_eqns); + _m_K[i].resize(m_n_eqns); } +/** + * @brief Destroy the RK45::RK45 object. + */ +RK45::~RK45() {}; + /** * @brief Estimate the current integration error. * @@ -132,17 +387,61 @@ bool RK45::step() bool step_accepted = false; bool step_rejected = false; + double h, t_new; + std::vector y_new(m_n_eqns, double(0)), f_new(m_n_eqns, double(0)); while (not step_accepted) { if (h_abs < min_step) { return false; } - double h = h_abs * m_direction; - double t_new = m_t + dt; + h = h_abs * m_direction; + t_new = m_t + h; if (m_direction * (t_new - m_t_bound) > double(0)) { t_new = m_t_bound; } h = t_new - m_t; h_abs = std::abs(h); + + auto y_and_f_new = rk_step(m_fun, m_t, m_y, m_f, h, _m_A, _m_B, _m_C, _m_K); + y_new = y_and_f_new[0]; + f_new = y_and_f_new[1]; + + std::vector scale(m_n_eqns); + for (auto i = 0; i < m_n_eqns; i++) { + scale[i] = m_atol + std::max(std::abs(m_y[i]), std::abs(y_new[i])) * m_rtol; + } + double error_norm = _estimate_error_norm(h, scale); + + double factor = double(0); + if (error_norm < 1) { + if (std::abs(error_norm) < std::numeric_limits::epsilon()) { + factor = MAX_FACTOR; + } else { + factor = std::min(MAX_FACTOR, SAFETY * std::pow(error_norm, _m_error_exponent)); + } + + if (step_rejected) { + factor = std::min(double(1), factor); + } + + h_abs *= factor; + + step_accepted = true; + } else { + h_abs *= std::max(MIN_FACTOR, SAFETY * std::pow(error_norm, _m_error_exponent)); + step_rejected = true; + } } + + _m_h_previous = h; + m_t = t_new; + for (auto i = 0; i < m_n_eqns; i++) { + m_y[i] = y_new[i]; + } + m_h_abs = h_abs; + for (auto i = 0; i < m_n_eqns; i++) { + m_f[i] = f_new[i]; + } + + return true; } diff --git a/src/RK45.h b/src/RK45.h index d1a6e26..88a10ab 100644 --- a/src/RK45.h +++ b/src/RK45.h @@ -53,15 +53,15 @@ class RK45 { public: RK45(RK45_ODE_SYSTEM_TYPE fun, double t0, std::vector &y0, double tf, - int largest_step=std::numeric_limits::max(), + double largest_step=std::numeric_limits::max(), double _rtol=1e-3, double _atol=1e-6, - std::optional first_step=std::nullopt); + double first_step=double(-1)); ~RK45(); RK45_ODE_SYSTEM_TYPE m_fun; int m_n_eqns; - int m_max_step; + double m_max_step; double m_t_bound; double m_rtol; double m_atol; @@ -87,14 +87,12 @@ class RK45 { { 3.0/40.0, 9.0/40.0, 0.0, 0.0, 0.0 }, { 44.0/45.0, -56.0/15.0, 32.0/9.0, 0.0, 0.0 }, { 19372.0/6561.0, -25360.0/2187.0, 64448.0/6561.0, -212.0/729.0, 0.0 }, - { 9017.0/3168.0, -355.0/33.0, 46732.0/5247.0, 49.0/176.0, -5103.0/1865.0 }, + { 9017.0/3168.0, -355.0/33.0, 46732.0/5247.0, 49.0/176.0, -5103.0/18656.0 }, }; // size _n_stages x _n_stages std::vector const _m_B {35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0}; // size _n_stages std::vector const _m_E {-71.0/57600.0, 0.0, 71.0/16695.0, -71.0/1920.0, 17253.0/339200.0, -22.0/525.0, 1.0/40.0}; // size _n_stages + 1 - std::optional _m_h_previous; - std::optional _m_t_old; - std::optional> _m_y_old; + double _m_h_previous; std::vector> _m_K; std::string _m_status; @@ -115,193 +113,3 @@ template int sgn(T val) { return (T(0) < val) - (val < T(0)); } - - -/** - * @brief Empirically select a good initial step. - * @details The algorithm is described in [1]. - * - * References: - * - E. Hairer, S. P. Norsett, G. Wanner, "Solving Ordinary Differential - * Equations I: Nonstiff problems", Sec. II.4. - * - * @param fun Callable, right-hand side of the system. - * @param t0 Initial value of the independant variable. - * @param y0 Initial value of the dependant variable. - * @param f0 Initial value of the derivatives. - * @param direction Integration direction. - * @param order Error estimation order, i.e. the error controlled by the - * algorithm is proportional to 'step_size ** (order + 1)'. - * @param rtol Desired relative tolerance. - * @param atol Desired absolute tolerance. - * @return double Absolute value of the suggested initial step. - */ -double select_initial_step( - RK45_ODE_SYSTEM_TYPE fun, - double t0, - std::vector &y0, - std::vector &f0, - int direction, - int order, - double rtol, - double atol -) -{ - int n_eqns = y0.size(); - if (n_eqns == 0) - return std::numeric_limits::max(); - - std::vector scale(n_eqns); - double d0 = 0.0, d1 = 0.0; - for (auto i = 0; i < n_eqns; i++) { - scale[i] = atol + std::abs(y0[i]) * rtol; - d0 += (y0[i] / scale[i]) * (y0[i] / scale[i]); - d1 += (f0[i] / scale[i]) * (f0[i] / scale[i]); - } - d0 = std::sqrt(d0); - d1 = std::sqrt(d1); - - double h0; - if ((d0 < 1e-5) || (d1 < 1e-5)) { - h0 = 1e-6; - } else { - h0 = 0.01 * d0 / d1; - } - - std::vector y1(n_eqns), f1(n_eqns); - for (auto i = 0; i < n_eqns; i++) { - y1[i] = y0[i] + h0 * direction * f0[i]; - } - fun(t0 + h0 * direction, y1.data(), f1.data(), NULL); - double d2 = 0.0; - for (auto i = 0; i < n_eqns; i++) { - d2 += (f1[i] - f0[i]) / scale[i]; - } - d2 = std::sqrt(d2) / h0; - - double h1; - if ((d1 <= 1e-15) && (d2 <= 1e-15)) { - h1 = std::max(1e-06, h0 * 1e-03); - } else { - h1 = std::pow((0.01 / std::max(d1, d2)), 1.0 / (order + 1)); - } - - return std::min(100 * h0, h1); -} - - -/** - * @brief Performs y = alpha*A*x + beta*y or y = alpha*A^T*x + beta*y - * - * @param trans Specifies the operation to be performed. If equals "N" then - * y = alpha * A * x + beta * y, if equals "T" then y = alpha * A^T * x + beta * y. - * @param m Number of rows of the matrix A. - * @param n Number of columns of the matrix A. - * @param alpha Scalar alpha. - * @param a m x n matrix of coefficients A. - * @param lda First dimension of A. - * @param x Vector X. - * @param incx Increment for the elements of X. - * @param beta Scalar beta. - * @param y Vector Y; will be overwritten by the updated vector Y. - * @param incy Increment for the elements of Y. - */ -void dgemv(const std::string trans, const size_t m, const size_t n, - const double alpha, const std::vector> &a, - const int lda, const std::vector &x, const int incx, - const double beta, std::vector &y, const int incy) -{ - - int lenx, leny; - if (trans == "N") { - lenx = n; - leny = m; - } else { - lenx = m; - leny = n; - } - - int kx, ky; - if (incx > 0) { - kx = 0; - } else { - kx = - (lenx - 1) * incx; - } - if (incy > 0) { - ky = 0; - } else { - ky = - (leny - 1) * incy; - } - - if (std::abs(beta - double(1)) > 0) { - if (incy == 1) { - if (std::abs(beta) <= double(0)) { - std::fill(y.begin(), y.end(), double(0)); - } else { - for (auto i = 0; i < leny; i++) - y[i] *= beta; - } - } else { - int iy = ky; - if (std::abs(beta) <= double(0)) { - for (auto i = 0; i < leny; i++) { - y[iy] = double(0); - iy += incy; - } - } else { - for (auto i = 0; i < leny; i++) { - y[iy] *= beta; - iy += incy; - } - } - } - } - - // Form y = alpha * A * x + beta * y - if (trans == "N") { - int jx = kx; - if (incy == 1) { - for (auto j = 0; j < n; j++) { - double temp = alpha * x[jx]; - for (auto i = 0; i < m; i++) { - y[i] = y[i] + temp * a[i][j]; - } - jx += incx; - } - } else { - for (auto j = 0; j < n; j++) { - double temp = alpha * x[jx]; - int iy = ky; - for (auto i = 0; i < m; i++) { - y[iy] = y[iy] + temp * a[i][j]; - iy += incy; - } - jx += incx; - } - } - // Form y = alpha * A^T * x + beta * y - } else { - int jy = ky; - if (incx == 1) { - for (auto j = 1; j < n; j++) { - double temp = double(0); - for(auto i = 0; i < m; i++) { - temp += a[i][j] * x[i]; - } - y[jy] += alpha * temp; - jy += incy; - } - } else { - for (auto j = 1; j < n; j++) { - double temp = double(0); - int ix = kx; - for (auto i = 0; i < m; i++) { - temp += a[i][j] * x[ix]; - ix += incx; - } - y[jy] += alpha * temp; - jy += incy; - } - } - } -} diff --git a/src/RK45_wrapper.cpp b/src/RK45_wrapper.cpp index fc825cb..6f1a3ab 100644 --- a/src/RK45_wrapper.cpp +++ b/src/RK45_wrapper.cpp @@ -9,59 +9,82 @@ #include #include #include +#include #include -struct RK45 { - static const int n_stages = 6; - const std::array rk_merson_nodes = {0.0, (1.0/3.0), (1.0/3.0), 0.5, 1.0}; - const std::array rk_merson_weights = {(1.0/6.0), 0.0, 0.0, 2.0*(1.0/3.0), (1.0/6.0)}; - const std::array rk_merson_starweights = {0.5, 0.0, -3.0*0.5, 2.0, 0.0}; - const std::array, n_stages> rk_merson_matrix = { - 0.0, 0.0, 0.0, 0.0, 0.0, - (1.0/3.0), 0.0, 0.0, 0.0, 0.0, - (1.0/6.0), (1.0/6.0), 0.0, 0.0, 0.0, - (1.0/8.0), 0.0, (3.0/8.0), 0.0, 0.0, - 0.5, 0.0, -1.5, 2.0, 0.0 - }; -}; +#include "RK45.h" extern "C" { /** - * @brief ODE solver based on Merson's Runge-Kutta 5(6) integrator. - * @details This method implements Merson's Runge-Kutta 5(6) method with adaptive time-stepping - * to solve an ODE. + * @brief ODE solver based on Runge-Kutta-Fehlberg 5(4) method. + * @details This method implements the Runge-Kutta-Fehlberg 5(4) time integrator with adaptive + * timestepping to solve an ODE. * * @param rhs Pointer to the method to compute the right hand side of the ODE. * @param neq Number of equations and unknowns + * @param u0 Initial solution/guess. * @param data Extra data to be passed down to the RHS. + * @param dt0 Initial time step, optional * @param t0 Initial time. - * @param u0 Initial solution/guess. * @param tf Final time * @param itf Maximum number of iterations + * @param usol Vector holding the solution at a given instant. * @param rtol Relative tolerance for the integration. * @param atol Absolute tolerance for the integration. - * @param dt0 Initial time step, optional - * @param usol Vector holding the solution at a given instant. + * @param mxstep Maximum allowed step size. * @param success Boolean indicating whether integration was successful or errored. */ -void rk_merson_wrapper( - void (*rhs)(double t, double *u, double *du, void* data), - int neq, - void* data, - double t0, - double* u0, - double tf, - int itf, - double rtol, - double atol, - double dt0, - double* usol, - int* success +void rk45_wrapper( + void (*rhs)(double t, double *u, double *du, void* data), + int neq, + double* u0, + void* data, + double dt0, + double t0, + double tf, + int itf, + double* usol, + double rtol, + double atol, + double mxstep, + int* success, + double* actual_final_time, + int* actual_final_iteration ) { + std::vector y0(neq, double(0)); + for (auto i = 0; i < neq; i++) { + y0[i] = u0[i]; + } + + RK45 rk45(rhs, t0, y0, tf, mxstep, rtol, atol, dt0); + + int itnum = 1; + while (itnum < itf + 1) { + bool successful_step = rk45.step(); + + if (not successful_step) { + *success = 0; + return; + } + + for (auto i = 0; i < neq; i++) { + usol[i + (itnum - 1) * neq] = rk45.m_y[i]; + } + + if (rk45.m_direction * (rk45.m_t - rk45.m_t_bound) >= 0) { + break; + } + + itnum += 1; + } + + *actual_final_time = rk45.m_t; + *actual_final_iteration = itnum; + *success = 1; -} // end void rk_merson_wrapper(...) +} // end void rk45_wrapper(...) } // end extern "C" diff --git a/tests/test_rk45.py b/tests/test_rk45.py new file mode 100644 index 0000000..7de8192 --- /dev/null +++ b/tests/test_rk45.py @@ -0,0 +1,40 @@ +from numbalsoda import rk45_sig, rk45 +from numba import cfunc +import numpy as np + +@cfunc(rk45_sig) +def f(t, u, du, p): + du[0] = u[0]-u[0]*u[1] + du[1] = u[0]*u[1]-u[1] + +funcptr = f.address + +def test_rk45(): + u0 = np.array([5.,0.8]) + data = np.array([1.0]) + dt0 = -1.0 + t0 = 0.0 + tf = 10.0 + itf = 1000000 + + usol, actual_tf, success = rk45( + funcptr, + u0, + dt0, + t0, + tf, + itf, + rtol=1.e-8, + atol=1.e-8 + ) + + print(actual_tf) + print(usol.shape) + print(usol[-20:,:]) + + assert success + assert np.isclose(usol[-1,0], 0.38583246250193476) + assert np.isclose(usol[-1,1], 4.602012234037773) + +if __name__ == "__main__": + test_rk45() From d8fad64a1eda98222dd5e521399f9e1367b066c9 Mon Sep 17 00:00:00 2001 From: tbridel Date: Mon, 5 Dec 2022 15:44:30 +0100 Subject: [PATCH 4/8] Put the RK45 test back with the others. --- tests/test_numbalsoda.py | 19 ++++++++++++++++++- tests/test_rk45.py | 40 ---------------------------------------- 2 files changed, 18 insertions(+), 41 deletions(-) delete mode 100644 tests/test_rk45.py diff --git a/tests/test_numbalsoda.py b/tests/test_numbalsoda.py index 68a1417..0ac7b05 100644 --- a/tests/test_numbalsoda.py +++ b/tests/test_numbalsoda.py @@ -1,4 +1,4 @@ -from numbalsoda import lsoda_sig, lsoda, dop853, solve_ivp +from numbalsoda import lsoda_sig, lsoda, rk45, dop853, solve_ivp from numba import njit, cfunc import numpy as np @@ -30,6 +30,22 @@ def test_dop853(): assert np.isclose(usol[-1,0], 0.38583246250193476) assert np.isclose(usol[-1,1], 4.602012234037773) +def test_rk45(): + u0 = np.array([5.,0.8]) + data = np.array([1.0]) + dt0 = -1.0 + t0 = 0.0 + tf = 10.0 + itf = 1000000 + + usol, _, success = rk45( + funcptr, u0, dt0, t0, tf, itf, rtol=1.e-8, atol=1.e-8 + ) + + assert success + assert np.isclose(usol[-1,0], 0.38583246250193476) + assert np.isclose(usol[-1,1], 4.602012234037773) + def test_solve_ivp_1(): u0 = np.array([5.,0.8]) data = np.array([1.0]) @@ -58,5 +74,6 @@ def test_solve_ivp_2(): if __name__ == "__main__": test_lsoda() test_dop853() + test_rk45() test_solve_ivp_1() test_solve_ivp_2() \ No newline at end of file diff --git a/tests/test_rk45.py b/tests/test_rk45.py deleted file mode 100644 index 7de8192..0000000 --- a/tests/test_rk45.py +++ /dev/null @@ -1,40 +0,0 @@ -from numbalsoda import rk45_sig, rk45 -from numba import cfunc -import numpy as np - -@cfunc(rk45_sig) -def f(t, u, du, p): - du[0] = u[0]-u[0]*u[1] - du[1] = u[0]*u[1]-u[1] - -funcptr = f.address - -def test_rk45(): - u0 = np.array([5.,0.8]) - data = np.array([1.0]) - dt0 = -1.0 - t0 = 0.0 - tf = 10.0 - itf = 1000000 - - usol, actual_tf, success = rk45( - funcptr, - u0, - dt0, - t0, - tf, - itf, - rtol=1.e-8, - atol=1.e-8 - ) - - print(actual_tf) - print(usol.shape) - print(usol[-20:,:]) - - assert success - assert np.isclose(usol[-1,0], 0.38583246250193476) - assert np.isclose(usol[-1,1], 4.602012234037773) - -if __name__ == "__main__": - test_rk45() From 3e854126e5db3cc85b7963888f22d03e4e04a72a Mon Sep 17 00:00:00 2001 From: tbridel Date: Tue, 6 Dec 2022 07:49:37 +0100 Subject: [PATCH 5/8] Fix RK45 to pass low tolerance test. Return also the vector of evaluation times for consistency. --- .gitignore | 3 + comparison2scipy.ipynb | 154 ++++++++++++++++---------------------- numbalsoda/driver_rk45.py | 11 +-- src/RK45.cpp | 1 + src/RK45_wrapper.cpp | 10 ++- tests/test_numbalsoda.py | 7 +- 6 files changed, 89 insertions(+), 97 deletions(-) diff --git a/.gitignore b/.gitignore index b0d46e5..d7e60f5 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,6 @@ Photo.run build *.a *.dylib +_skbuild +numbalsoda.egg-info +.vscode diff --git a/comparison2scipy.ipynb b/comparison2scipy.ipynb index b115853..4ba9019 100644 --- a/comparison2scipy.ipynb +++ b/comparison2scipy.ipynb @@ -2,11 +2,11 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "from numbalsoda import lsoda_sig, lsoda\n", + "from numbalsoda import lsoda_sig, lsoda, rk45\n", "from numba import njit, cfunc\n", "from scipy.integrate import solve_ivp \n", "import numpy as np\n", @@ -15,7 +15,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -31,19 +31,22 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "funcptr = f.address\n", - "u0 = np.array([5.,0.8])\n", + "u0 = np.array([5., 0.8])\n", "data = np.array([1.0])\n", - "t_eval = np.linspace(0.0,50.0,1000)" + "t0 = 0.0\n", + "tf = 50.0\n", + "itf = 1000\n", + "t_eval = np.linspace(t0, tf, itf)" ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -59,22 +62,9 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAFHCAYAAAABXlMnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACGB0lEQVR4nO29eZgkV3Xm/bsRudRe1dVdvXertbT2FQQCxA5GBlmAbDBgwGCw8YYXbMOMbGNj42Hs4RsG2xhjZrAxxoBBgFkMyAIkARYSau1bt9bel9q3zMo17vfHjciMzIqIjMilKiI63uepJ7uzMrPiZtx7zznvec+5QkpJggQJEiRIEAZo630BCRIkSJAggYXEKCVIkCBBgtAgMUoJEiRIkCA0SIxSggQJEiQIDRKjlCBBggQJQoPEKCVIkCBBgtAg1es/sGnTJrlnz55e/5kECRIkSBAR3H333dNSygmn3/XcKO3Zs4d9+/b1+s8kSJAgQYKIQAhxyO13CX2XIEGCBAlCg8QoJUiQIEGC0CAxSgkSJEiQIDRIjFKCBAkSJAgNEqOUIEGCBAlCg8QoJUiQIEGC0KClURJCvF0IIR1+fm0tLjBBggQJEpw+CFKn9FJgxfb/p7p8LQkSJEiQ4DRHEKN0l5RyuWdXkiBBggQJTntEIqd0YmGFlVK18cnDd0JuZn0uKEGCBAkS9ARBjNKTQoiKEOKAEOJXe3ZFTShXDZ77P7/Pr3zG1qqoUoR/fAX83bPX6jLWDeWqwa/+yz6+89CJ+pOVEnz6Z+ChL6/fha0RlosVPv+Tw1SqRv3JShEe+TpUy+t3YWuEqiF5+PhC45Mr83Dvv4KU63JNCRL0En6M0gng/cBbgeuAO4FPCCHe08sLs/DAUbUgf/TEdP3Jqf3qMT8d+4V5y/5Jbnr4FH/1nQP1J5/8Phz8IXxzTW7BuuKzdxzihq88yKdvP1h/8v4vwBffCj/8yLpd11rhvTfez7V/8yOemrIx51/5Ffjab8CJ+9btutYK77vxfn7nC/c2PnnTH8G3/9v6XNAaompIPvKfBzg6l68/KSXs+ydYnly/C+sxWholKeVNUsq/kFL+p5Ty21LKXwS+CPyxEMLx/UKIdwkh9gkh9k1NTXV0gRm9/idq3vLMk/UXrMx19Plhx9PTOQAKZRt9OW/2MqwU1+GK1hZPTqrN2PoeAJh+TD3OPL4OV7S2+Pp9xwF45MRi/cmTD6nHI3etwxWtHQrlKl/cd5Sv3Xe8Pv+lhB9/DO78BORn1/cCe4wHjs7zN99/gnd+2sYSnXoYvvm78E+vXLfr6jXazSndCIwDe5x+KaX8pJTySinllRMTjt3JfeOSnaP8+WsuAmAub9I1K7bJOHewo88POw7PKi9puVipP2l5SdUSGIbDu+IDa/ynFm0GePbpxscYI206ZU9O2oyyZuqT5l0bLccC9gjBmgcs2WjsmN//+47MAzCTs819iyWaeWLtL2iN0KnQYU24s01DWQCml82bY4+OctMO74gPTiwUAFgqVKga5te9fEo9SkNRmDHG1JK655NLhfqT1pjtG1QMsVQos2JGCHP5knpSSsiZTknMx398vn7PZ5bN8edszMvC4TW+orXFsbmV1U/aDXFMc6rtGqWfA6aBNXHVNg5mANvEXJmv/zLm9N1MrlT79+KKOQnthnjx+Bpf0dpi0jRKlnECoGAm/mNO30wv1+/9gnXvi4tQMTfrxXgbpZMLdaM0a62DvE1xu3B0ja9obXHKnPMzuVI9dWFniWKaV2pZpySE+DLwE+ABQAfeYP78tpRyTbijoT51mTUKa2UO0gNQzkNhfi0uYd0wawvdF1bKbBjMQHEJhGZGSvGVxRfK1do9b6AvLaeksgLlFUj3r/3FrQGa7z3Q6JBZEXNMMbVcH3/tu7A7IjF3Sk4tKqMspTJMW0b6mhzyWRjdsT4X10P4iZQOAO8Avgx8CbgQ+EUp5d/28sLsGMoqo5QvmRtTYQHGdqt/xzxSmsuV2bNxAIB5a2MqLcHITvXvwoLLO6OPxYIa74aBNLliBWkpLQsLkB1R/47xxmQxA6P9aeYt+q5kqvD6x2N970FFR5mU2qJq+WT7/Y752p+2GeWlgrX3zddfENO570d994dSyvOklANSyn4p5TOllP+yFhdnYSCjjFLO8pZLOcgOQ99orCdmqWKwXKxwxsZBQOUY1C9yMGoZpfn1ubg1QL6o8ilbRvowJCq/UimqCGn8TPWilXguTKhTVmdNDNroO9Moje5Q9z7GJRFz+RITQ1kyKY2c5ZAWTUM8dkas5z7A4kqFLSMqn56zs0T94+a/4zn3I9HRYTCrA5CzujqU84q+6xuLtVGyIoWtI30A9a4WxeW6UbKH8zGDRdltNse/XKgo6hLUpgSx9RahTtltH+snb937kjn+0V1gVNRaiCnmciU2DKYZzOg1B4XyiqKuhzbHeu1LKVlcKbN9TFHTdaM0X3fIYjr3I2GU+tM6QtgjpTxkBqF/Q6wnphWyW96SpcSitAyDm0DPxJrCse73lmE1/uWibRMe3aUeY+otgnJKdE2wcTBjM0qmNHzEzCXE2CmZy5fZMJBhIJOqR0rlAqT6VbQQ47VfKBuUqkbNKC1Ze1/ZxpLEdO5HwigJIRjMpMgVmyKl2Bsl5SlP2CMlw1BGKTOk6MsYUxjWRrzFipSKFeWQQD3BG1NvEZRTMpRNMZBJ1R2SGn0X/5zicrHCcF+KwawtUqqYwpaYr32LJdk+quZ+zSGvFBVDlB6IrUMSCaMEisKr3Zhy/rSYmLVIyYwU8qVqPVLIDKrJGfNNCWDziMP4a5FCvI3ScF+K/rROqWKoOrVSs1GaX7fr6zXyxQoDGWWU85ZRLhdsa39+Xa+vl7BTt2AzSrXxj8fWIYuMUepP6xQqFoVxutB3amJaOZVaoh/MiTkW64VpLcRxs06tUK6qnAKosacHIR/v+z/Sl2Ygo3KqK2WbUR7eph5j7JTkSlWGslakZHNIU33q/hcXY1tAahmlbaPKKC3XIsUCpLIwsCG2DllkjFI2pdf7X5VzptBhVC3KmCqQFs1IaeNgBk2Y9J1VOKlnzEhpft2ur9ewhC11o2TUjVJ6QCkwi4tub488Fs1Iqc80SvlSRXWIB5Xoh9g6JVJKcsUKAxmdgUyqXqdWKUC6TzmkEFujbBXKW/nkYsWk7qvFek4tiZTWF9m0RrFiKM/IqJib0pD6d7XU+gMiCIu+U95yStFXVTNSSvXVjXJMYUVKGwdtC9OKFNL96v6X4nvupKLv0gyklVEqlAy1KQu9LguO6f0vVQ0qhmQwm6Ivrau1D2ax9EDdKMWUKbEipQ0DGTK6phwyyyFNZaFvJLYOWWSMUl9Kp9h8YzJD6t/FeG5MFn031KcW5krZ5imnMqcFfZdNabWSgAb6Lj2g7n9M7z1Y9F2KfitSKpsOmOWQQGwjZUvUNJjRyegaJcsoVQpq/Nlh9f+YbsxWpDTSnyab0tTct/a+dD9khmM79yNjlLJpTeWUapuyzShZtRsxw1KhwmBGR9cE2ZRGqSLrkZKeVQszxpFCrlSpecqA8pYbIqV4j98SOvSl1TJV9G1ROSR6SuXUYhopWVHyQDZFJqU1RUr9sXdIF1YsliRF1ooUaw55n8kSxHPfi45RSmkqUrKoOj2jbgzEdmIuFcoM96UByKQ0SlWjLnRIZZXYo1qqG+qYIVesMpjVa0apwVtM9cU6UpJS1u5/RlfjL1dNp0RXdGac6UurHGAomzIdMlvxrLUpQ2zHv1ws05/WSekafWmNop0lsM/9GObTo2OULPVdLVLIqE0ZYjsxLU8Z1GGH5UqzUTIpjJiOf7lYYTCTos/sf1YoG3W1lZ6OtbeYK1UxJAz3pWr930rW/U+ZRikzVC+mjRksYcNARq87ZGAKHQbqcz+mTslyscqg2fOzz9r7aspb0yhLm5MWI0THKFmRUgN9F+9NealQqXVIT6eEWph2+i7mRjlv0ncpXSOlCRUpGaZR0tKxjpSsfOJwX7pmlMrVZqM0GNvxW82Xh7KpxpxSeUVtyqfB3B8yc6n1vc+KlPpjbZQjY5RqCpwafZeOPX2XK1VqHdIzumZuSjahQ8zHv8pbLBtQNaXBtUgpnmO3cirDfSnSugCoz/8afTcc20ipllPKqEjRkKgzhSpmm6GY03e5YrXWiHpVpJTK2sYfP6YgMkZJeQvVpkgh3hMzX6wyaE7MtK41Jjsbxh/PjSlfbPQWCxUrUhKg6cpbrBTqhipGsHIqAxmdrEXf1SIlVbdFZjCWmxLU1XdDWTt9Wa1HSmkzUoqpQ6ZYAjX3+9JaU41evIUeETJKVk7JllOI8Y0BFSkNmBMzkzIpjOrpoz7MmW1mQBnlStXMKelK/BFnb9HalAcyKdK6Sd/Vckqqw0ecc0oWfTeQVZJwgFKpqPIoqX7QNHP8MV37trlfaxzQrL6DWI4/MkYpk9IoVyXSXqcU4xsDylu2IqU6fdekvoPYbkzLxTp9mU4JpT4zKiqfBLF2SqxNeTCr1yMFK6eo2yKlGI4d6m11BjP1SKlcsJUDgJlTjJ9DAkroMmhjCcrVJkl4klNaf6Q1xatXyzZJeCqrNqiYGiWrzQo4REp6NtY5JSmlMsrmwkzrpgKrWlY1OhBrpyRXcoiUmoUOMc4p5UsVNKGoq5pRrhklM1KMcU4xbypPQc39clWqZqxQV99BLFmC1HpfgF+kzYlZLRfURVveYjaeCqyqISlWjAb6qjFSygDx3ZSLFdVmZsAeKVYMlVOqRUrx9RatBqQNkZJF39kjpZJZqyLEel1qT2CVAwghajm1ctE0SilbpBRTo6wipfraV/feVqcklLMWx7kfnUjJ9BYrZVtOBWLLK9vpG7BFSs3FoxDL8VvqqyHbwqwYUkVKWnOkFD9v0R4p1XIqFn1nzykhY3n6bN6mvLTGXy2YBihtG38MN2V7M1qATEo0OaRJTikUyJiyWMMKYWveYlyNUn1TAou+ko0dLdL96mjoGHqLdvUZQFo3F6ZRqQsd4pxTsheP6vZIqdSovoNYjr9Z5ANQKduUtxDb4ml7M1qAlGayJA3qO4sliN/4I2OUUpa3VLZ1dIDY0ne5YmOkVGu1YtE3QqifmHqLyw6RUsnqEr8qUorf+HOlKpmURlrX0DRBSjONckObofgWj+dsOZVsSq2BGktid0piOPfzxWaHzMwpVYqAUOtfTykaMzFK6weLvjPsQgc4bSIlS33YUDwJ9bxCzFA3yvbxN0dKMc4plVQzXgv1vEIzfUc8779N5FKLlCpWNw+bUxLLsTfO/Vo3l4rZ98/KH8Z0/BEySiZ9V8upxFsWW9uUbfRVqVbRbjdK8ZyYVk7Frr5bJQmPc07JVtEPNqekuXgWYjv/B20OGUDV6mZiGaWYRko5mxweVE6tXqOXqb8wpuOPkFFyiZRienxBLVKqJXt1qobEsEuCwYyU4pdTao6UajkluyTcUiHFcGHaK/rBjJTKFaU+XEXfxe/+50sOQodm+i47rKIHo7oel9gz5GyFw6DuvSHBsM99SCKl9YZllGSlKdkZ20ihKVJKqUhRlouN3lI2nod9LdciRbvQo0kSLkRsF2au1BgpZVNaPZ/aHCnFMlKsONB3VqRk0bfxbMqab4qUag55pVQfO8T2oL8IGSWLvispxZnlMcSUvqslO5u8RaPsFCnFcfxNOSWrTsveZghiuzDzxeZISTSWA0Csex82ttmxIiWrxZiNvoPY3f9cUzlIbe+rVurUJcRWfRgZo5SxR0oNkcKQ8p6tCComaI6ULG9Rns45pYqVU4o/hdEcKWVSmnJIoLF4FmK3KRuGJF+uNohcwKSvwCZ0iKf6sFaj2JRTk5VSI32X5JTWFyndfmPsm3I8efVV6rsGo9wUKcVwYi4XK6R1UZMDp1PCJVKKZ07NSX1X65BvLxyH2I1/pVxFSptDVqOvbGdpQWwjpeUaS6LGn9JsRtlO32Xj2dEiMkbJCmFVnYZ9UxpQjzG7OflShZQmal5SQ07NHinFtP9Z3kbfgFqYq3JKEFujlCtWa9QtmN5ycz41lVFRU8wonGZJdC1SqjQLHeIpic+vyqea+eRqpWnvS4zSusLalHFSn0Hsbo6SBNc95VoI7xQplHNgGGt9iT3FcrFaK5wFW53SabIwnSIlw37Ao4UYjr+W6G8SOhjNdUpxXfsmS9KfbqLuq2V1jpgFe+/DGCF6RmnVphxPCsM6CtyCNTEbOhpAffwx63+2WhJtHV3RPP74CT0MQ3VIb1bfrcqpQCzzCsu2U2cBUppACDAM8zBHe50SxG/tF5VDopknI9RYkmb6LjOI6n24sg5X2TtEyCjZ6TunSCleC1Mlum2Rkt0or5qYxG78y030XVrXqBrSOVKM2aa0Um6MFMC8/0ZTTgUUhVWO1/itfKoVKQshyOiayieDraNHPOe+6vvXOPcBD4c8XuOPkFGy3xg7fRHPEF5Jgh0iJaOyWoEDsRt/znbAHzTdf62ZvovXoqwVTzYZZWnRV3pTpBizSClna0ZrIaNryiGB+NN3xWoTdWvmlOzdTCAxSusNa1MW1dJqTh3iNzGbIqXaptwsiY6pt2g/4A9s6kOjvNool/OxqupvzqmAmv/SGuMq+jJuc7+xGS9AShdq7kM9WkjH0yjlS40sQZ0lqazOKUHsxh8Zo5Qy+VVRLblESnHblOu9v6ApUnKk7+I1MZeLjeOv07cu449RTq1lpLQqUozXvW8uHAfQNa1ulKzx6ylVSByztZ8rNjpk6drad6CuIXb3PzJGyboxwjVSiNeNyTdJgq1NWThFChC78eea6MvTaWFaOZVmp6RGX63KqcVrU15uakYMav6LapPQAWIbKQ76pa4hdvc/MkbJCmFF86YU0xA+51Q8CYqmOg3oO/tx0NC8MONtlOvNeO30pY2+ivmmnHeIFFO6AGlFis2y6HiNP9fEElgskVr7trFbdVoxyykGNkpCiB1CiGUhhBRCDPXiopxQ25Rk06Yc0xA+33R0gW5OTM2peBRitTDLVYNSxWgwynWnJP45tbxDol/XNDR5ehil5WKVjK7VKWsgrWkqUtLS9fOEIJZCl3yz8jbl4pDHcO1De5HSh4E1nwW6WasgmjdliN3ClFKaIbwtUjJbjQjpUDwKsRp/87EVYDkl8rRYmDkH+i6tCzV2WH3/Y3Z8Q3ONGqhIadW9h1jm1FZR17pbPjl+ax8CGiUhxAuAnwb+v95cjjfSuoaQ1cacCsTOKBUrBoZ0oC8AcRrQd83NWEFtyhpm5XrMF2advrJHSsIWKcXcKDexBGD2f2tmSSB2OTXlkFZXd4jHK1KKz/ghgFESQujA3wJ/Dkz37Io8kNE1tGb6BmIXwtebsdYnpsUra7Jp/KmsOuguVpuSQ6SU0khjSYKbuoQDFOPT/6128miDJFpTDgnEXhacK66OlGqRoqNRis/YixWDqiEdJeFCNuWUUn3qGJ8YjR+CRUq/BvQBf9eja2mJlC7MGxNv+s6peDBl0ldaM30nROwojGUHo5TRNVJYm3K8I4V8qYIm6ucIAaQ1QUqY4487fdukPgM1/1fNfYitQ+okctKa6bva2o/P+MGnURJCbAQ+CPyelJYExvP17xJC7BNC7Juamur0GmtIWxMz5t5SbWJmG4UOOmbTVUejHJ+J2XzyJqhIsWaU4r4pF6sMZlIIW0Jf10U9UnSkL+MUKTaqz8CkL5s3ZYjd2q85pA7lEEK65dTis/bBf6T0P4A7pZTf8vNiKeUnpZRXSimvnJiYaP/qmpCpeUtO9F2MJqZDTiGt2zZlrZHaiNvCrEdKjQWEaVw6GkCsFma+VGmQg4MSujgb5ThGilVn+q6ZvoLYzX2nGjUrp+ScuojX+AFSrV4ghLgIeAfwQiHEmPm0eYgRo0KIqpRyTdrUpnWBVna7MTHalBxyCrpm85RXeUvxmph5hzYzalN2GL+eUfMhRuPPlaqOkYLu5JTEMVIsrY6UUpqZT3aKFCoF80iTlttZ6OHokGn2nFJilAD2Amngxw6/Owp8Cvjlbl6UG1K6hn465JScIiXNJacC8YsUm44uAJVPrOVUVvHq8br/+aJDpKQLW6QYbwWWOuBw9fgVde+w9kF1StdH1+gKe4d8aXU+VdOEckpcc2rxmfvgzyj9CHhJ03M/Dfw34FXAU92+KDeonNLpEMI3njwJamKmhZVTahp/dgiWTqzV5fUc1nHQQ01tltJO9BXEbmHmmhpygiqePV3ou+Y6HTAjJUfq3jb+vugbJUt5aR3wZyGjSQTS2SjHaO2DD6MkpZwGbrU/J4TYY/7zh1LKNXPR0rqHt2AVEDZv2BFErtaQsnEsfZpplJzouxi1GrHUZ33pesoz1RApxpy+LVXZOJhpeK4xUoxv78dK1aBYMVbTd7owHVKHtQ+xGb8TdQ3Qr1trv2nuZ4dgJh5jtxCZ3ndgKnBwoe8gdhOzeWH26V7qu3iMHeodwu3qs5TulVOLlwIpV2w85A3qQhcp9NVtdiA24885KE+hhfIWYjf+Zoc0azmkp0FOqS2jJKX8tJRSrGWUBJARUlX1u07MeNwc9xDeLVKKGX3lQN+k3eqUIHYLM19qPOQN6vSdbJ77qYz6PmIy/rpD1jx+N5YkXms/X3RxSDWHbiYQu7UPEYuUsm4hbAxD+IGMjqaJhudr9J1jTm0ZpFyjK+wtmtusQHOdUrw7euSKq3NKaU3l1GTzpgSxMspOdTpgUve4dHOBWI1fiNUOab/ulk+N19qHiBmljPAIYSE2G5M6dXZ1uq8ewjt5ixLKa6LM7zkcE9265qw+g1htylJKxzodSxIunXKmMfKWczWRS7NTorWIlGK09tOrHdJMLVJy2vvis/YhYkYp6xrCxi+Eb96UALLCg76D2IzfqaI/pQlSIv51WsWKQaWp9xko+jJNFSncIqWYbMoO5QBgCR2M2FP3qnB6tUPa75pTildOESJnlNzUV/HalJcduiQDZFzHHzNvsVh1iJRc6nQgVpGCVdHfrL7STfpyVU4JlAIrJurLnMv4VU7RS+gQj/uvWkw5OKSaRzkExGbtQ1SNkmutQjxuTL5UWUVfgD2nFPOF6XCeTmObHRcFkmGs0RX2Dk7NeKEuCTeEg1GKUaTodGwHqEhZp+q+KcfEKOcdatTAzhLFe+1DxIxSRpweN8Ytp5Rx6hIN8Ru/Q05J00R9/G68eiX6vLpTM16o12nFPafk1CEeVE4xJSurhR6pjGo1FROHVLEEq++xu/I2XmsfomaUNA/6BmJzY9xySu7JzniF8LlidRV9Ax5Cj2x8vOVlr0iJikekFP2xg3PfRzAjJWGoOq1mxCxSdIyUXB2yeK19iJpROk3Ud/lWkVKMhR5VQ7JSrq7alMF2/2PMq7tV9Kc1JXSIO31XM8rp1UY5TQXDKacWo0gxV3J2yDKeylti4ZBZiKZRas4pxKxT9HKx4jIxW+XUoj/+nMumDC2q2iEe47daTDl2CTcwHCOF+GzKbjV6Vk4x/pFixdEhc82nZ+PFEkHkjJJLpBCzTtHWwmyGe04tRpGCy6YMp0dOrd4l2qlLuAd9V46J0MOFJbDUl6dDpNhMXYKftR+P8UPEjFLajb6D2HiLpYpBuSqdJ6ZX8SjEYvxO58lYyJ4GOUW3Oh1LEl5125QByvleX17PkXPJp6Z0DT3mRskqnPaMlFzXfvQdUguRMkqu9B3UvcWIw00SC5B2C+HTA4CIxcJ0y6mAx/2P0cL0rNPxkoRDLO6/dRR8M9I1oxxf+rJUVYXTng5pM0uQ6gOhxWL8FiJllNJuxaMQG29p2aUhI9iFHk70ZTwW5rJLpACnh9Ajb/Y+sx/bAfVIyTWnBPEwyh6Rkjd9F/2x16lrD4e0ee+rrf3oj99CpIxSCrUpO1a1p+NhlNzqVADSbrJQMBfmUi8vbU2Qczjgz0KmZVV79O+/dRS6/dgOsCThVSoxj5TyJeecSloz0IR0py9jMPacy5E14JFTgsQorScsT7mKW61C9G9MvUuyg7fkFsJDbBamW6IfPCLFmlGKg1F2FrnUJOGOcz9GkVLJhb4z576zUYoHS+DLIY3x2rcQKaNkTcyK04G5MbkxtYnpSd/Fd/xuFf2gnBIDDbSmaZvKgtBjMX63OhVdV212nCOlGEWKbkbZnPtVpy0rJupDL4fUtZsJxGbtW4iUUUq1mpgxuDFuvc8A5+OwLcTEW8x5GKW0qFJxihRilFPLFyvOUbJZp+PKEkA8IiVXSbTFksRXfejlkHqzJPGY+xYiZZTS5tEFzpFSPG6MV/Foa/ouDpuSmexNO3iLbpJoiM/4XdrM6OYhf95GKdrzX0rpeMAj2CIltzZDEPnxt++QxmPuW4iUUbKEDhXX/lfLkT+BsbYpOy5MF/UZxCpSdKroBzV+x00ZYjR+56ML1HlSVXeHDCI//mLFoOoqiW7hkELkN+aa0MGrHMJp7cfo6BKInFFqIXSQBlQKa3xV3UXeQ4GTooqBWJ1TgVhFik6LEtT4HTcliI9Rchm/ZtXpxJi+86Kv6tS9l9Aj2vffckgdnRKsAy6TnFKokKoJHeI/Mfsd6Cu1KXtFCtHelMC9QziovIJjlAyxMcp5l+JRMO+/0/hjIvTwpK88135cjLIldHA3yo7lMDGZ+xYiZZQs+qosXYQOEPmJ6UVfpahSli3oq8jTl87qKzhNjHLJWegAKqdYcbr/MRF6eOdTLfou/g6pUz413cooxyB1YSFSRknHhwIn6hOztPoocAuuOQVQ4zcqUC318Op6D7eGlGBtyi7jjwGvLqVU6jOXSElvdf8jbpTrkmgHoUerchCI/PjdOqSDLVJ0mv/WIZfl6B9yCREzStaNKTlGSvHwlvKliiOnDJCSVSputyw243en71Ki4jH+6PPqxYqBIZ0T3WAa5RiPv97Nw119503fRXz8Lh3SQUWKhhSUWW2w4iL0sBAtoyTjzyvniu4T06KvqoZDmB6b8XvQd7JK2TVSiAF95dEhHSkVfRtj9WG9GbFDTkVa9F2MHVKXvn8AOgZldCpVp7WfGKV1g4afnFLEJ2bJa2KqTbniVLkek/HnSs4HHIIaf8ucUoR5dUt95uiUGGpTds0pZocjT18uFz3Ud9bajzF95xUppahQRadSje/atxApo1SbmK68KpG/MW4V7WBOTKm18JaiPn7vSNFT6BFxXr3eId5hjIZHpACxyCktF8oADPU55ZQ8hA4xUh+6UvemQ1b2ZEmiPX4LkTJKOlUMKajI+PKqbg0poR4pVGJK3xmGVJGSw6YEtKCvom+U8x7Fk1TVhu3okEEs6LtcrSGpgyRcWg5pnNWHVUeRB4AuK5TRqSb0XbigS5Xo9t6Uoz0x8x45FTUxU84hfDYGm3K5ipTOiW5QkaJ3pESkF2ateNJp/K3ouxgYpeVihYyukU2tHqPuVQ4CsYgU88WK69zXMaiiU/ai7iNO31qIllEyK/odcyrpAfUY8YXpJQnXqVJFcxE6mEYpwhPTSvQPZR1aqQC6rFJqaZSie/+9Ev21SMlLfRnhsQMsF9zzqZqMv1HOe+SUrEjJkbqPgUNqR6SMkmZGSmWnG6NpkT/oz6pT8Y6UWvHK0TVKSwUP9Rkqr+DpKUOkx++V6PcXKUVb6OGZT7XouxirD3Me5SA6FSpSp+wodEiM0rpBlx6SaIh8CF+qGlRcGlKCNf6UC68c/UihHil5eYspl0hxWD1G+P57HXCIoSIl70gx+kIPr3sPXkYp+pFivuidU4pzPtmOSBklTSpZpKO3AJH3lvIeDRkBNGkpcOJJX/oxSmphxlMWW88pOdF36rtxLByHWHjLXkbJKgcpGfHMKZUqBqWq4b72MdTcd9r7Un0gtEjfezsiZZR0U33lyKtC5L2lnEdDRjDpSzdJuKYrwxThhbnkccAfKKPsXkAYfaOUL1XQBGRTDsvSsIxSnIUe7vSdFnOhh2c+kbpD5pi6iIn60EKkjJImK1RlfOm7nFdOASun5iL0gMgvTCtSGnaRhOtS8epxrdNaNvveCeFQ8lCj71rl1KI9frdIyRd9GfGxg3ONFqi1X3ZjCcA0Sku9urw1RbSMkmEl+mO6KXvlFGgRKUDkx7/cMlLyuP8xiBTyRXflZb1O6TQ1ShZ96UbfRbyjhTX3h92oa8N0SGO69u1oaZSEEK8TQtwuhJgRQhSEEAeEEH8shMisxQU2XIusmq024nlj8l45Beo5NcdkJ0Q+hF9ukVOqRYpO9z/dD4hIj9/r2IoafeeaU7GEHtH1lnNeRtkPfRlh9eFyoYVDZpaDxDWfbofLDGjARuAW4MPAPPBs4APAVuDdvbowJ2heChRQm3IMvCU3SXht/J4TM7rjzxUrpDThnFOhPn7HhRkDXt3r2Ar/OaVojr/WzcPVKFv0nQO1CY3qw8xAby6yh2hF3wlDFY7H1SG1o6VRklL+Q9NTtwghRoDfFEL8lpRr55oIw6KvXDbl7DAUo+spWhNzpM+5eNSi7zxzaoXFXl1ez6GKJ11yKth5dY/xR/z+u+XTLPqu6KU+g8huTLVuHi3G726Uba12ImyU3Og7TZbNfLLH3F8+2avLW1O0m1OaAdacvtNk2Zu+6huB4mKEQ3izIaXLxBRG2Ww1Ek+hx7LHUeigcoqqTs3NKYm2t7hU8Er0q02rYMSz91+uRT7RGn/RbfxZk76MqFNi0XeekVIrliTCLJEdfug7AIQQOpAFngH8NvD33YiSFhcXmZycpFwut37xBe9B7K1yZV+eRx99dPXvx38KrnkWPPqI0u1HDBcOlPm/r97GicNPMpNOs3nzZkZGRmq/16Tqku26KWeinuwtu2/KUtZySo6yWIg8r75cdG9GW+vo0CpSiuqm3CKfSLVMFc25GTPYjFI0mYLW+dQW+fSIO2R2+DZKQA5llAA+A7y30z++uLjIqVOn2LFjB/39/a60TQ3TOsvFKstDZ7J1tM/hCqdh4Qhs2Qv6mgdyHePEwgoDyyUu2D7CysoKx44dA6gZJmFUzGSnV6QYzU0JrES3d6LftU4LIs+rLxUqrvRNjb5z25TT/aClInv/WxVOU5v7blGy6bxFdPy1FlsuOUVhsgTuyuNoz307goQTzwNeAPw+8BrgY24vFEK8SwixTwixb2pqyvUDJycn2bFjBwMDA60NEoAEiUDisilZ0ZHbjQs5qoZEEwIhBAMDA+zYsYPJycna71UI7yELzQ4rTzGi41eRgnM+zdqUK6RalAREM1KUUpo5JZfx1+grlyUrRP3+RxCt1GfKKLm0mILo03emHF7TnPdBYZR9lINEV31oh2+jJKW8R0r5IynlR1D03a8LIc52ee0npZRXSimvnJiYcP3McrlMf39/gMuVDQ+roJlettm8MWowDNBtk7K/v7+B1hRSeYuuBXTZYZQCKZoe07JH6/4afdWqo0dEN6VC2aBqyJb0nWtOBSIt9PFD31WES0cDiL5R8sonosphXBuygk19mO/NBa4h2k283GM+ntnpBfiKkCxIiRTC1SYhzAVrRNQoSYndUWr+bkSrArqIUxg5z4r++smjcVRfLhW9RS51oYPHesmORnb8VuG4V0eHqlc3E2vuR1R96plPBDDK3uUwNaMcTabAjnaN0tXm49PduhB/kIDAVV9Ri5SiSV9VDdkQKa2C0aLViDUxo7owC+69z+r0nYf60FJfRhAWfdVSEl71WLLZ4Ujfe/Cm7wzhUQ7RZzlk0Rz/kpdDRj2n5O6QRXv8dvjp6PAdIcQfCCFeKYR4hRDiz4D/DfyblPLJ3l+iDVKaOSUXWDmlkNF3H//4x7n22mvZuHEjQghuvfVWx9dVpcopOUJKs6OFy8m7EOlISUqreNK795nrybugxl8pQKXUo6vsHVrSVzZJuKtTFmWjbHYz8WozVCXlTl+lskrcFMG5D6ocxMsoUS27N2SFyEeKdviJlO4C3g58CfgicB1wA/DW3l2WG9QNcc3lhZS++8xnPsPs7CzXXHON5+sM6REp1boke9B3EfYWV8pVDOmdUwBUQ9YYGmVLfdXKKLmePAyRFjrkihV0TdCXdtmSjDJVr0gJIk3fevb9QzUOqAoP+rK29hd6cHVrCz8dHd4PvH8NrqU1WkVKmhUphYu+u/3229E0jYceeojPf/7zrq8zDOmqvrE2Ze9IKbq1Gn7oG2ghdLAvzMGN3b7EnmKpRfFk/Th0ZZRTTnqHiG/KgxndPcdcLVMVHjVqEGmjvFxonVMyhFedUrSFHnZErMJUKumrW6gkNEBb00jpxS9+Ma973esanrv11lsRQvDQQw8BoGmtv2YpJVUDdLeUkmHblGJYq7Hc4tgKe06pZU4twuN3azFVF3p4UFjZ6NaptYoUVE7JQxIO8R5/tYwhPNoMnWb0XXggTaGD12s0LXQ5JT+QEiQekZJpaD37X8VgU3ZvSGoTOsSQV2/VYsoafxWP4uHsMFRLUC704hJ7Cq8D/oCa0MHVIENkjVK9Rs3t3lcBqYyS2/gjTN03I1pGCZVV8qwPE3oki0er5qB0D/oCaFHVHl31XauzlKzzdDyFDn3RjxTd6UvllJS9nJK+UfUY0fF70letIgWILH3XMp9qRslS81CeRpglaUaQNkNrgj/7xsM8ctxlYpWWqTBDRRynL+1SRFjOK4ovdSLw375w+wh/et1Fgd/XDRjmZHOPlOr0jSuFoemRLSC1Tt119xZ9SMIjLItdKlbIpjQyLsd2WE6JgfBBXy7CkHvRehjRmr4rY2it6LtoSuJbNWOt3XuRdnfINB3S0T4lwELkIqXWEO4dH0KMlpGSuSlLr6p2ML3F6Clwls3iUT91Si1zahFcmMsFD/oGzER3ChA+kt3RG//iStk9nwZQrSCFRz4NItv7callOYBplDQP5S2Y44/e2m9G6CIl10hFSjhxH7PaOHP6Rs6eGHJ+3exTUCnC5gt6d5E29PX1USo11sXMzs4G/pyWkVLVCuE96CuIrAJrqVXxqCX0kB6y4Ahvyp7HVoDKqWjq93GVxHsbZVPo0NIhW1J7RZBOMeuM1oXT6veIlDtLAOr+R9Aha0YEIyXROqe0hpLwnTt3sn///obnbr755sCfY601d/WdzSjFUIHU0ihV7eozl/Gn+8wCyugtzNY5FRUpAN5tliCSG9NiocxIv0ekZJSRmg9JuFFWBdQRQr1w2lt5GVeHtBmhi5RcYVmiVh6Qpq+pJPz666/nU5/6FO95z3u49tprueWWW7jpppsaXrNv3z4OHjzIkSNHALjtttuYnp5mz549XHnllYCfnJJJ32keBXQQWV59caVMNqWRdSzAwZ8kHiLrLbZqyIlRQZqRUtyakpYqBoWy4X5sB0C1gqENup8lBo2RYjpIo+f1RevC6fra9zTKfSNQiD59F6FIqX4zPM8WFKYkfI1auF977bV86EMf4sYbb+T666/n0KFDfPSjH214zcc+9jFe//rX83u/93sAfOADH+D1r389H/tY/fSP1jklNXGVLDR+Ve0tPWU/ve8gsuNf8jq2AmqRAuBR1W+p76JllJdMOXzLSKnl3I8mfemnQzqA1NItHNJoOmTNiE6kVDNKLeqUrFZD0qj/u8e44YYbuOGGGxqesxvOT3/603z605/2/Ay/OSXRir6LaLJ30UdOATA3plbJ7ugtzOVimWEr0nFCtWIzSvHKqbWkbkHdf83jLC2I7PhrNWotji2RrYQOEXXImhGdSMlG33kGQSFtNdQKVSkRQrg3ZLUiJc1DFgpmTilaixL8qK9Mb1H3kVOLoLfYus1MBWo5JZfxp7KgZyO3MS1akVKL++9LEg6Ru//1GjXvs8TQWqkPRyO59psRHaPUECl50XfhbMraCoYh3ak7qPHKeBXQgVqYpeXIjX+pUGlJ3wAILR27qn6ror9VnU4tUmqV7I7Ypuw/UkpTrkp3+j6iObWlQoWMVz61aq19H22Wyvn66yOK6Bglu9DBM1KK5umzVVkP8pxfYE1MH5ESRO5Y8MVC2R9956tWI1qbcrFiUK56nDoLavy6MtqtuxpEa1NeXPGRU6rWjbLr8CNqlBYLZUY9HTJz7utpb4c0wh1N7IiOUQqaU4pYpNA6UjLHE1MKY3Gl0rJ4Eqxkb7xazdQihRbqM1oJHSCSRtlvpCQ1NT9cI+WItllaWPFnlERLhzSaObVmRMco1WxSq5xSNE+frXodWwF1+k73UasBkVuYS4UyI/3e9BWA0P3k1JbWTH3ZDSz6VJ/RShIOkaQv/Y2/gjDXduvi6WjJopVD5t33DwDd4+gKiHRHEzuiY5QC55QqPb+ibqIqJakWR6EDir5r5SlDpDamQrlKsWL4EjrQSujQN6IckgjRlzX6ylMSbouUYmeUKggBQ24d4kHRdxZ92UroEbFaHT+FwwBCz3irD2PSKTw6Rsm3+i6iOaVW9F3VSvS3koWaFEaEFqZF33h6i76FDtGLFBd85lTqOaV4CR0WV9RR4K2YAmHRd17j7x+L1NwHH/SdvRzED0sSsfvfjOgYJXvxrNfLhDmkiOWUWtN3ZqSUahUpWEZpvmvX1mv4Kp60+n/p6dhRGJZR8s4rVBGtIgVQ9z9C9x5M5aVXlCilyinpymnxzKn2jcHKfFevr9doWQ5h2Oe+l0MSzZxaM6JjlBrUd170nTDPVIqOUTKkxJAS3YdRapns7B9TjxFamIu+Et1lEDop3ePkWYhkVwNr/C1zan6EDv1jauwRmv+tlZdqLForoQOYkdJ89y6ux5BSslio+MqnanqLhqwJfbfWaFTfebYa0lKRWpSW5+dplKr1RL93/6sx9Rihhekrp2LSV2ldtE70Q6QonEVfkVKlFin5u//RGf9SoVWkUM8nQrwipVypStWQLei7ek7JVzlIhO69E6JjlEwjJMy8i/eR6DrIcAgdTpw4wXvf+14uu+wyhoaG2LVrF29729s4fvx47TWGH6Nkq+r2XJSpDKQHIrUw65LgVon+NLomvMffv0E9Rmj8Cytl+tIexZMA1QpC87Ep1yLlue5dYI+xuNIiUrDV6EEroxwt+tK3yAVaU9fpPkj1RWr8ToiOUcJG3+HnSPRwREp33303X/3qV3nTm97EN77xDT784Q9z55138rznPY/lZaUQq0VKnnVKJn2XSnsnekF5ixGamHVJcIuNSU+R0j2Og4f6phyl8bfKKYBK9Os+6KsIRspLxXJrhwQ198GHUY5QpOBXDg8Wfedj7UfIIXNCdBqyyjp9Zz5h+3cTNB3KJeffrTGe//zns3//flKp+lf9jGc8g/POO48vf/nLvO1tb6t3CPcRKWlamkq1xdiitjB9eYtl0BR9Vyj72JQjFCm0VF+Bou9SPjo6RDBSbFmnY9TVZ+DDKBcWwTBatEgJBxbyPqhbi75LpakaRe8PjFhOzQnhv2tNEH4iJU1fM0n4i1/8Yl73utc1PHfrrbcihOChhx5ibGyswSABnHvuuQwMDDA5OQkEzCml0t6eIkTOW1oqVNA1wUDGm75CT5PSNG9ePZWB9GCkxr+w0qJOBaBaDkbfRWRjklKyVGgRKZlzX/MbKSEjU0BbE7n4iBS1Vr3/QDklEZr7ToiOUbI6NPjNKRlrd6ZSUDzwwAPk83kuvPBCwKdRsofwXpsyRM5bstRXolVDWi3VWugAkRx/60ipWtuUfdF3EdmYcqUqhvShPAR1qjCt2iyNqceIjL/e9691Tk1LqfHHSejhhPDRd9/+73DywdXPm8ccD+n9nFUBPaO7n0JbLUG1CJkhXCk+J2y9BF75l21dtl8YhsHv/M7vsHfvXl7xilcAwXJKWquOBmBOzAe6cblrAl85FVN9l9I0700JIrcwF1bK7N3scZYSmMWjKTTRok4pYpGSVaPmHSlZc9+H0KE2/qhESv6UlwBaqn6elqsmpn8DnHqom5e45gifUWoJQYs4yWasPPJO64QbbriBH//4x9x2222k0yYdISWaEN7Fs1UVKaR0PXaR0lKrA/6gllNK6cJ7UwaTwohQTilf9s6pQN0o65p3sjvdr1rtRMQoL674PLYCakIP70ghWsXjVuF0q2NLAHQzUixXDfrSLlapfyxSc98J4TNKbpFKbhoWjpAbOZdD82XO3TLsfmNW5mDuIEycrxZpSPDxj3+cD3/4w3z+85/nqquuqj3fspsD1CTRqVaSaFCRQmm5oTVNmLHYqk4FzJxSipQmWiuQ+sdg9qmuXV8vYRiSpWLFF32HpsZfjRF9aW3KY/0Z9xfVikfjR18urqhztFK6RybFihTTPjp69G+I1Np3QnRySlh1Sj4ueQ2Pr+jr66NUalTDzc7Ornrdl7/8ZX7rt36L//W//hdveMMbGn7Xsu8d1BpypnTNR6Qwph4jQmGoRL/fSElrvSlHiL5bKlaQsoUkGGo5tZQmfNK38926xJ5iPq/WzthAN4UORMYoK4fMx9xHoFs5RT85tYisfSdExyg1ScK9OzqsXVPWnTt3sn///obnbr755ob/33rrrbz5zW/m3e9+N3/wB3+w6jOqRosWQ2Ae8mYm+v3kVCBCG1OZDQMenjLU6Es1fh9GOSIUhq8D7qB2yF+6VZ0WRCpSmjcl0Z5GqZZTCdDRIiJz35fy0nRI0+Ye4SunGJH574Tw0XeusCIl+/9coK1dpHT99dfzqU99ive85z1ce+213HLLLdx000213z/66KO89rWv5fzzz+cNb3gDd9xxR+13ExMTnH322VQN6R2+Qz2npGmqP6UX5Rchb1FKyXy+zKjXpgTqXvqRhIMaf2UFKkV1nEGI4asZK9hyij7p26UT3bnAHmN+xYqUvOg7yyhlAMN7/JlB1WYsIpHCos9yACufCD7oO4iMUXZC9CIlk77z7uhg2to1OFPp2muv5UMf+hA33ngj119/PYcOHeKjH/1o7fd33nknCwsL3H///Vx99dU897nPrf188IMfBJTQoXWkVE/0g88QPgITM1+qUqoarSMlo74pt6QvIzR+f33vqoA0c4parCTxc/kyKU0w6FmjVm9ICi0k4UJEqtXQYqsO6VDLJ6cDrf0kUloDBOx9B2tWQHvDDTdwww03NDxn0YsXX3wxb3/72z3f74++qye6QXlLroKdCEVK8+amvKFVpFQtQ3ZY0VcthQ6WtzgHw1u6cJW9w0KQ3mearoyyL0l8NCKF+XyZsYFM6xo1QE9lgFJroxyhnNriSpkLt414v8ik7nVf9J059yOw9t0QoUjJfKzxd62Or9BC0//OC1JKf0IHW+83aNFqJkLe0lxO0TejXuorqEeKmh9J+Jh6jMDCrNWp+Ej0K/rS5/iLC5GY//P5knc+Cerqs5rQIT45NV8iHxt1Dz6OLoHIGGUnRMcomTVH9gokT4SoKasXDD9976BBEg74O1MpAgtzwXekVKkZ5YrRotVKX3R49Xqk1LpOx6LvfEVKEIm8ihK5+FAeAqmUVacTj0ipVDFYLlYYb0ldN9J3LQ95hEg4pG6IjlGSEoTwFSgBkTlTyVeLIbBJws2J6Xl8RRZS/ZFYmHN5H4luaIiUID6y4IWVMppoVTxZP7bEX/HwmHqMwPjnV8qto2QzUtTTPtrsQGQaEtfk8IM+xt/Akng4JXoaMsORuPduaGmUhBCvF0J8XQhxTAixLIS4WwjxprW4uEYEjJRCdKaSFyzjkmpllMyJmTZD+LjIgufyAXJKetqfUbbnlEIOSw7vmVOp0Xfm0R1+1HcQCafEH33X2Put5dzvG4vEvbfmvu9IydwjfAldInDv3eAnUvo9YBl4D/Bq4Bbgc0KI3+rlha1C8yF/rUIlLbUm6rtOETRS0v1EChAZCmPB9BZbS8KthenDKNcojPkuXGFvMZcvsaGVp2yLlNKa8JdTgUg4JYHoO7+R0sC4GnvImZJZM5/qa/xayp8kHCJVp+cEP+q766SU07b/f18IsR1lrP620wuQUnp7ifVXujdgdUJE6LuKh1FqMLzNkvBWE3NgPBITcy5fZiCje5+6CjYKwwevrumQjYYseDZX8uEp109e1TU/XdLH1WPI73+hXGWlXG1N3VqRUtrqEt5q7m9UpwoUFtQ6CCks+q6lU1LLp/qQhEPkDvlsRstIqckgWbgX2NzpH0+n06ysrPh7sTTpO3PvbjUv1fEVldAeX2HBapnjRN+trKzUmrZaxaNpP7wyqIWZn+nqtfYCvro5QEObIfCxMPtHQ78pA8zlymwY9FE4DLX737J4eGCjesyvbncVJtT63vmk79JmIXRL+q42/nDP/1nLKLU0yiXQMzWWwF9D4vkuXOH6oF2hw/OARzr945s3b+bYsWPk8/nWdJwZKfnOKmlWAW24o6WKIRE0RkpSSvL5PMeOHWPzZtP2V8uqTsVPrQJEyCiVWnczADNSyvhrtQIqWgj5pgxqYxr3k+gGW52Sj00JQn//ayIXP+UAqEhJCD/0XTSMkq8WS6DGr2dsLIEP+nYl/HPfDYGLZ4UQLwNeA7yj0z8+MqKKxo4fP065XPZ+cW4aqmWq04JTCwWK02lOeSmWSjk1KWcfCXW33Pl8iZVSlf1Ljd3M0+k0W7ZsqX1HNUm4n0Q/mEZpNvTHQqucil+j5LPVCkTCKEspmcuV/HWzAP91WnpKUTg5J5IjPJj3LXIxGx6bOcXW1HU0jNJsrsRARnc/7cCCOffrHR18rn0ZMOUREgQySkKIPcDngK9JKT/t8bp3Ae8C2L17t+dnjoyM1DdeL/zbW2D6CWbf/gN+5oM384HrLuTtl5/p/vonvgtf/Xl4x3/C7ktbf/464Tf/9R72n1zke7//DO8X1rpEW5uyDwpDVlURpeU5hxDzK2W2jbU4XkRKk8KoS8Jb0neDm2Dm8S5dZW+wVKxQMaQPo2SThPupUwI1/pBvypZRailyMYtnrWihpdAjIkZpLu/DIQE199MDAdb+JrVfFBfrop8IwbcLLYQYB74NHAbe4vVaKeUnpZRXSimvnJiY6PASTRjVxkR3TEL42ZwP+gZWJzv9eou5cI9/Pl9mzM9ZQsgGCqM1hbMp9PSd1c3CV6Ibave/ZaQEkYgU5/3WqNUiJd2f0CMiObW5XBCWIOOvzRAohwRCHym7wZdREkIMAN8EMsC1UspcT6/KCbVWG/EySnP5UutFCfXiWdNbarkpD4Z//IYhmffjLVr0ldklHPwku8fVYWflQheutDewJMHjLYUO9Y4OabOjRUtYFE6I4bvvoZlTQQjSutZ67mcGVPF4iOc+KOWpv0jJrFH002IMlEMGoR+/G/wUz6aALwF7gVdKKSd7flVOqFX0B6CvIPQ3xpckGIJ1CYdIjH+pWMGQftRXtpyCH0k41L3FfHi9xXpOxW9OSdWptZz7oIxyiMcOyiHL6Br9fnIqmpojuuajIS1EwigHou/sQoeWa9+UwUc0UvKTU/o48Crgd4BxIcRzbL+7V0pZ7MmVNaNaacwptNqU0gOQ6gv1piyl9Fc8CXVJuF9ZaASMkn/6xoqUMv5arUCjtzi6s5PL7BnqkZLPnJKZ7G6Z6AaTvpwJdbJ7wTxHq2Wdou1o77Qf+g5MoxzeuQ8mfdfKIYNapFgvHPfrkIV7/G7wY5ReYT7+tcPvzgQOdu1qvGBUIJVB0wSaH1moEEoWHGJpZK5UpVyVrekbqEvCazmV6EdKNUmsn0POwDTKQXNq4fUW54IUT4JZEuCjTgnU+KslRWFmhzu80t5gernERl/51FLNKOl+DjmE0OfUKlWDxULFn0NqO+DReq8nBsLPEnihpVGSUu5Zg+toDaMM2iCA2f8r+iF8LdEdJIT3uylnBk1ePbwTcyanguyNQz4T3UEk4RHwFmdzJVKaYNirtAFs48/6q1OCRqckpEZpJldk05CPk4GtnBKYknCfa3/+cIdX2DvU82nt0Hc+cmrpgVA7ZF4IbwFLM5pCeH8KpHCH8L7pGzDHn/VPX0HojfL0shp/y42pgb6LT07NErm0pq8so5TxV6cENgVWeMc/s1xq7ZBALVIAYqM+rFPXfliSShN952fth78kwA3RMUrmyaugkp1xCOFn/dI3sLpOJwZGecY0Si03pgb1nUlf+jlTR+ih9hZVOYBP6hZqkaJvhwRCfv+LbBz0ESlV65GS7zqtgY2q/5tFfYYMszmzQ7jvtZ+qpS58O+QhnvteiJBRqntLaT1ACB/iRembvqsVjwao04HQj396uchARmcg45O+sqsvW21MmhZ6ozyX8ysJrkdKVqK/ZVsuS4EV0vGvlKrkSlWfkVI9p5ROaZT8bsoQ2v6Hc3773kFt7UOA1MXgplBT916IjlGy0XeBQviV+dD2v6vRdy0lwbbiUb+SeAi9UZpZLvqnb0Btyn6Lh8GkMMK7MH31vYMGo6Sb97+lTxLySMnKJ27yc/+NSm3tZ3WNciX6kaLFErRWXhqqM0stp+Z370vou97DLB4FK4T3m+yVoe2YO5cvoWuCYa+jsKEh0R9oUx7cFO6cQq7kn76BRqGDXwonxOOf91sO0CD0sO5/i/FnR1RtT0iNco269XX/S7U6pXRK+GdJILQb8/SyT5GPrUYNVKTkyyEN+dr3QnSMUkMI67OAzupqkFufet9WmFlWxXNay1NnbYnuoJtycaG+qYcMU0s+1Vf2TTlITm0wvJFi1ZD+C6ft9J1fBZYQoY6UfSsvoVHkpGuUAhmlcBrlqaUio/1pf+eIgS2n5rdObRzKOSj7PBooRIimUfIbwg5tUY/L4TRK08tFJoaDRQrB6Ktw5xVmciV/9E0DfedTEg6hpu9mcyUMCZtHgtz/+vh9U1gh9ZZ9Ky9htVHyM/ZBs+dm5Nd+3SEB5ZC3FPlAvVYpgmKHCBmlMpiHfKX9KpAGzbOIQjoxp5aCT8yMuSn5W5jhHb9hRgq+PGW7+s5vqxWwHd8Rvpzi1JKKFCb8RopCa+h/5ovCGtwEualOLrNn8E1fQa3FFkAm5VPkNLhJfWchnPtgsQRBHLJ66sK30AFC65R5ITpGqVKsV3X7jpSsTflUDy+sfUwtFf1vSgC6qmlJ68IfhRHiSHF+pUzVkMHoOy0ofbcJlVMMnwJrytyUfTslpqecSaklW/TjlAxvheWTbV9jLzGzXPKnvITG8es+zlMCdfL04ERo176KlPpav9BopO/SvkVe4S8ed0M0jJKUZlW3WsBKFunjxvSNqveEcGJKKZlqg74Dc2H62pQsoxS+8c/UPOVg9JWVU2vZZglC7ZTUIiW/99+2KYPPSGlos3JIWp7qvPbwrbyE2rEtoDZlX2OH+vhDiOllv9R13SED/NepWZHScjgjZS9EwyjZEt2gZJG+NiUhVLQQQgpjYaVMuSrb4pVVrUYQ+i583rIVKWzy280Cggsdhraqx6UQjn/JkkT7jZTq9BXgM1LeCpWCOuwtZPCtvISGSMl3TgnU2g+hQ7JSqrJcrPjPp0G9HMZvQ9phc+6HcO23QsSMkpqYvg76sjC0OZQTM5in3GSU/BYPZwaUNDiE3mK9m0Mw9V0goUPII6XBjM5gq753sGpTBihXAgh9lsI3ft+RAjTklHyr78A0SuGb+9OBqNtm+s6nJDwzCJnhUN77VoiGUaqYm5Jd6BDxiRk40Q0N9J2vnAKo8YcwUrDou0AUhnn6phA+hQ7DIY6U/FK30KA+q0dKPsQbITbKvlsMQcP4s36FDlB3SENGX9byiW2sfd8NeUHR90mk1CM43BhfbXYgvJFSzVsKtimDpUDyO/5wGuXp5RKa8HGWEjScJwRWp2gf47e8xTDe/6VCAKNUakh0A5T8REo1Cidc4w+kvIRVkvBAc98oh07oEjifCE1z369R3hpKh6wVImaUTKFDIPpui9Lqh6wxYz1S8qHAaTZKukap4lPmHFKjPJMrMj6oIp+WWJXs9Xn6KihvMYQL03c5ADQYpWygnFI4I6X5lTIVv8pLaBB6WMeh+3JKQzr+6eUA+cQm9Z3vFmsQ2rnfChEzSvVOwYFCeGTo9PpTy0UyusZIv5+cQpO3lApglIe3hm5RAkwu+uzmAA73PwCFMbQ1lJGi73IAWBUpgM/i2b4x5ciFbGM6uVAAYOuoD4cMGpsxp3y2WQJbSUS45r/lkPpuRgsN6jtfymMw53746MtWiJZRStXpq0DJTgjlxJwYzrY+SwccI6VARrm0DMXlNq+0Nzi5WPC/KVUb6TvfslgIJa9erFRZLFTaipQCqe8s9WnIjPKpJWWUtvjpZgGr6pTAb6QYzjq96eUiGwbqoh1POLQZCsQSlPNQXGrzStcH0TBKFeVZ2CkM/7LQcHY1mFoqsilIohsahB6BhA4Quv5/pxYLbB3xa5RKylM2DbjvNlNg8urhckisFjvBhA5N6rugyf4Q4dSCZZR83H8pG7qEW0bZV6QYVvpuqRSMuoQGSbj/fHo4c4qtEA2j1Fw82pZRCteNCUbfNKnvAimQwicLLlcNppdL/jYlaPCUIWCye3iLakwZIm8xUKIbGuuU9AAdHSCU9O3JRWWUNvvpaGCtfa2RvvQVKWVHINUXuvFPLbdDXduFHgEiJQgdfdsKETFKVqSkbmQgozQYTqM0vVwKtilBk9AhuvTlpLkp+88p1D1lCNAlHmwFtOEZfyCRC5gtthrpu0hHSotFNg5mamPxhLX2bSwB+KzTEiKUXR1OLhTYFmTuQ9MJCUmktP5wkEQX/S7KzABkR0PlLZSrBjO5gHUq0CQJj65RqiW6A9F3NqMUiL4LX1eLU2ak0A59F6ghL6iNKT9Tr/ULAU4tFvxHybUaRfX6miQ+orLoqiE5tVhg21iAuQ82+k4Lpr6DUI3fDyJilKyciplTMiOFlkdCWxjZDovHe3RxwXFqsYCUsN13on91CO97UxrYqPIxIZqY1qbcGX0XgL6CUI3/xMIKuibaou/S7URKEKqc4qkgIpeKmiu1tR90/MNbYelE0EvsGWaWi1QMydbRfn9vaFLfBer9Z6kvQ+SQ+UE0jFKz0CGtDsby7S2N7oCFo724srZwwowUto0FnJjtFM9qGgxvC5VRDiwJtkmiISiFEb5I8cR8gS3DWX81WuCsPvPrlIxsV4+L4dmYVaTk0yDXjJIVKQU0SqM7YeFYaGTRtbXvO1JspC8DzX0hzFql8Mx9P4iGUarRV2ZOKfDC3AGLx3pxZW3h+Lw6DdI3r+zQ/8p3ohtCN/5TiwUyKY0NA+nWLwa1MaXq35UepE6tf4N6b4iM8omFgn+HBJrUdxZ95XNjGt2pHheOBLnEnqFUaUPkAqvUh4HWfjkXmq4OJ4I6ZE1GOVCNJiiHNESRoh9ExChZkVJT/y+/E3N0p+oUbnkd64yatxSEvhOaOiMGyAQJ4cH0FsOxKYFSX20Z8VmjBeq+peqedTqILFaI0BnlEwsr/u89NERKQohgQpeRHeoxJOO32mv5zie6REqBWBIIzfhPLAR0SCvNQo8A+VQI3dz3g4gYpdX0FQSYmCFbmCfmVxjOphju8xkpNOVUAqkPQS3MxePgV7HWY5xcCFCjBKsipUCtVsA0yuGgb6WUnFgosL3NSAmC5hVGVf+/kIz/ZJAaJbAJHZrVh37zyVakGI61f3JBsQTjfo5sAWWUtHTNIQ1UOA51+jIka98PImKUGotHa7Ua5YC8ekgmpqJvgnjKzZtSwBB+dJcybCE5VyqQ+gpWR0q6zyOhLYztCs2mPJsrUawYbURKdQcmkFMiRKiM8mRQkUtTpJQJ0mYJbJFSOMZ/wpSDB2MJ6t9V2uz76VvkNbpLMU0ha7PmhWgYpUojfZdNBw3hTW8pLJHSQoFtftU34LgpVQyJ4ZfCGgnPwpRScmqxGDBSalyYgegrUAtz6WQoZNF16jbo/e/EKQmP0McqnPUvcmmShKeCSsK3KPVpSBzS9liCukNWP3nZ59of26Ue58ND37dCNIxSU5fwtoQOEJqFeWJhhe1BIiVb8SS0w6tbFMb6j39xpcJKuep/U4JVkVJfWqdQ9tklHczxS1haf7FDYJGLUQVZ7ZC+3Rkah+zY3Ap96YAiF2i/zZKmm+rTcIz/xGLAfGKTQ2ZJ4n0LnUImdPGDiBmlRqGD7xuTGVAqrBBMzGKlyvRyKZin3OQtBTq+AGwTc/3Hf2QuD8DODUHH37gwA6kPQ2SUrUjBN31bVkaMdFOkGCRSGjGFPuWC//f0CEfnVti5YSAYfQWr6LvARjkEc98wpIqUAq/9ukPSvlFa/7nvF9ExSlq61pAzsPoO1MIMwcQ8GVR5B+bErE/kQMcXgDLI6YFQTMyjNaM04P9Nzd5iWg9olMJDYRyfL5DWBZv8nrratClDm5EShMIpOzqfD+iQWONvFDoEL4lY/7k/kytRrso21n799X1mjaZvpqBvLFRCFz+IhlGqlBwjhaLfg+7AVKCt/6I8Pq+MUiD1VbnQ6CkHVSDVZNHrPzGPzinPf1cgo7Q6UgxE34WIvj2xsMLW0T40v4WzFTNSsie728kpQSjGf2R2JXiUDLXx96UCbsoQGvVpYDk4rKKurXy6b6NcE7qsv0PmF9EwSk2J/qw5MYNFSuFI9rY3MQurNiVoh8JY//Efmc0znE35O9zQQlOk1Bc0Ukr3qca8IViYJ+YDilwsyi1df0+g88QgNJHSYqHMwko5mEPSVA4SeFMGxZJUS+uuQDs8q1iCXeNBxt8099syyolR6j6aNuXAdUoAY7uhMA+FxS5fXDAcmVVGKVCk5Dr+gN5iGIzS3Ao7xwPkFMCRVy9VDP/qQwiPUZ4LSl9ZkYJdEi+CO2Sw7uM/ZkbJwajbxkgpm9IQAoqBhS6sO31rGaXdQYySg8gHAhqlEJVE+EE0jVI7kcL4mepx7uluXllgHJrNsXWkrza5fKFcaPCU+0yjVPBbpwUwdobq/2YlztcJR4NuykZVHYfdkFNqwykJgVEqlKucXCxwxvig/zfVNmXb/Q8aKaaySho9f8j/e3qAozWjFMQoW5JwtTELIRR9G2T8G/aox3Ve+0dm82wayjCYDcISrBb5QMBIcXSn6hRfyvt/zzoiGkapvKIS9SbaEjpsMI3S7PpOzMMzeXZvDOApgcordOotWeOfW7+NSUrJkdmVgPmkxjYr0CaFMbYb5g+va2POo3N5pIQzgtx/B/VdXyqgJB5g/CyYPRjsPV3Gkdk2lZdaqtbRANooCbCM0nqv/dl8MOoOHHJK7dB3ptAnIhRehIzSavoqkLcQEm/p0GyeM9qamPWF3J9RE3MlyMQcP0s9zj4V7G93EbO5EivlakeJbmgzr7BhjzLu63iExaEZk74JYpQc1Hf9GT3YvQfTKK3fvQcVKfWndf8tdkDlglKN+dfARjkzoGqV1nvtz+SDUXeg5r/uJPKKnkPuF76MkhDiHCHEPwgh7hdCVIUQt/b4uhrRJInOthMp9Y2os4XW8cbkSxWmlorBPGVYZZTrkUIb9OU6bkxHLOVdUE4dOo+UNp6tHmef9P+eLsMySoGcEgf1XV9aC3bvQW1MS8fXlcI5Opdn13h/8Hyi3mjE2h7/Os79ctXg+PxKG0apCzmlEMz9IPAbKV0EvAp4zPxZW7hESoFyCqAm5jp6S7VE58YAOQVYxSv3Z9T4A3nL/RtUc851HP+RmvpoHSKlEESKh2fzDGVTwSIFB/VdX1qnUAoaKVn07cFg7+sijpiFs4HQNPehDfoOzEhx/eb+8fkVDBnQIQOHOqWAfT+hvvZn4mWUviGl3CWlfD3wcC8vyBHllUZJrN7GpgRqYa4jr96WpwwOyc42vCUh1p3Csbo57AjUIbsx0Q3tymJ3qQLsdVyYh2Zy7G5HeQir2ywFqdGDulFeJ6dESsnhmRy7glC3YNYoNhrxbFoPJnQAGN+jTmAt5YK9r0uwHNLga7+5RlPN/UA1mkLA+NnrTt/6hS+jJKVc36qzykoDfWedKRPoxoCKlBaPrltjTitSCETfVStgVFZ5yhBQFgvrTmE8PZVj01DW/5Ed0L1ISdNVXmkdx39oNh+cunVQ3/WndcpVSSUIU7DO9O3kUpFcqcrZm4eCvdEpUgpaPA02oc/BYO/rEuosSXcipcD05cazY0ffrS+aOhqAySu3Q2FIQ6mw1gGHZvKM9KUYGwhA39RyCnVvqS2hAyhvef5I/SiQNcZT0znOnghKXTrklNrh1WFdI0XDkBydXQm+KTmp76yNKYhR7t+gftZp/E9OLgNw1qaARqmpQzqYbabaufewbhTe4Zk8GV1jy3CQI0sqqhmvA0sS2CEfP0uVRITkoFMvRMMoVRol4QADmVTwTXnD+tYqKU+53U25wzolMI1ydd2M8pNTy5w10YanDM61Gm15i0+tiyz85GKBUtUIVqMEjuq7KBrlJ6cVbXZWYKekscUUWJFSG3Mf1m3tPzWdY/fGAf/tpcA29xsPeNREO2v/bOWQr2NO0S96YpSEEO8SQuwTQuybmurCwXLl1SH8QEYnHzhSMr2ldcorHJzOtae8gwZPOaVrpHXRXqQE67IwZ3Ml5vPlDiKl1Zty8JziWVDOr4ss/KC5Ke9pp0ZN6A1ttqzxrwSd/+tI3z41tUx/Wg92lhAotWCTQ9pWTq1/g2pOuk5r/8nJZc4J7JCtnvtCiPaFHhCJvFJPjJKU8pNSyiullFdOTEx09mFGVfV/SjcmSPszevBFObRZTcyp/Z1dUxtYKVU5Mpfn3C3Dwd7oEClAuwWUpjR0+olg7+sCnpxS9M3Z7UZK+ur2/e0vzLXfmB47tQTAOVsCjt/BIasb5YDj33i2onDW4QiLp6ZynLlpMFikAMqJyDQ6MkoSHnDsABPnwfTai4eLlSqHZvPsDXrvHah7aOPoFqjLwiOgwAs/feeyKfen2yggFAImzoepA126OP94cmoZKWFvO4leWL0xZdowSjWj/Giw93UBT7VrlBwk0XVevc2FOf14sPd1AY9PLjPan2ZiyOeRFRYqq/Op/bVIKeD4J85XFM46bMxPTS8Hp+5AqeVWGSU9OH0FyihNPrrm9O3B6TxVQ3JO0LVv1ZRlGt/XVqQ0MG5GimvvkAZF+I1SbVNqDOH726HvADafrzblNZ6Yj08qT3lv0EjJYVOGNgsIhYDNF8Lk2keKT07lyKQ0dgSVBJeUMbNvTHUFUsD7P7ob0oPrEik/fmqZc7cMBZODw6rCcbALHQKOf/MF6nGNx18oVzk6txI8nwgqUnKg7wJHiQATF8DKrDrwcA3xhCnyCGyUyqZ8vWn8bUVKoO7/5No7pEHht6PDgBDidUKI1wE7gAnr/0KIgCR5QFRW51RA5ZQC03egvMWVuTWfmI+dWiatizYkwc4hfH+6zfFvPn9dvMWnppY5c+MgelD6xqorsRmltiMlTVPjP7W2pXZSSh6bXOKczQEdElA5RYd7D23klMbPVn3kJh8Jfh0d4NCM6vkXOJ8IKlpojpRMoYMMOoc3n68e19goPz65hBBtsAQOcx/ajJSgbpTWsf+jH/iNlDYDXzJ/ngNcaPv/5t5cmgkr0d/kLQ5kUuTLleCfN2FOzDX2GB4/tcyZmwZrZyH5hoP6DtpM9oLyFosLsHQi+Hs7wJNTOc7e3CZ9AyrCMWHllALTt7Au3uKMKfIITN2Cuv+rouQ21XepDGzcu+aRshUpBN6UpVSRcnOk0K7QZcKMFNdh/Ls2DAQ7GQBsRqnxewvcJd3C5gvV2l88Hvy9awi/xbMHpZTC5edgT6/QQX0GamEG5tShbpTWOK/0+ORScOoOWoy/zU0Z1nRjzpcqHJzJtRkp5FQ+Ta+3+9c0odSXxTacks0XqcPeltcuUrZEDoFFLqDoq1VCh06M8vlrHik9emIRXRPB6atKAZCOkQK0URIwvBWyo2ueU31icjn42MFmlFaXw7Q39y9Uj2t8/4Mi/DmlinNORdF3bdyY4a2qD9QaTsxCucrh2Xx7nrJDTgWsSKlNXhnWlMI4cHIJKeGi7SPB31zKrfKUAQazKXLt3P+aUV47Cs+KFAKrr0Dd/+zqRDe0sSmDihbmD61pu51HTyxy9sRgG5GClehvnPs1+rIdodPm89c0Uqoakqemc22ufWf6bjCbYrkto2TN/cQodQZX+k4nX64G55UtBd4aTswnJpXyri1PuWgZpcb39qXa6GgBMLgJBjat6cR85IQ67ffCbW0apczqBT2Y0VkuthMpWN7i2jklj51aYqQvxebhgMo7UPffQX0FbQgdwOaUrB1T8MiJRS5o5967JPoHs2r8bW3ME2srdHpyaplSxWg/SoZV938o26bIa2BcHeFxKjFKncGFvurP6EjZBq8MsOViOPUQGGvT0u+R42pTPn9rGxOzpKgfso3vbetMHQtrnFd55Pgiw32pYOcoWSgtr/IUwYyU2tmUhjZD//iaGuWHjy9y/raR4Mo7MCOlpnvfrtAB1txbnsuVOLFQaNMhsTblRqM0ZJ7c2tb933KREjqtUV7loWMLAFyyczT4my2WZJVRbnPug6m+TYxSZ9h+Ofz8Z+qFnyY6WpjbL4fi4pp1Nnjw2ALD2RR7grYYAuUpa6lVCqzBbIp8O/QVwLbLlAJtjXrgPXpikQvb3pRXq6+gg4UphNqYTj4U/L1toFI1eOT4IpfuaGNTAigurfKUrbmfa2fuj5+lRCMn7m/vegLiUTNK7ixSWk1fQZuR0vYr1OOJ+4K/tw08dGyRvrQWXOQBiiUQuuPab2vsoJySqQOqr15IEX6jNLwVLnwN9I81PD1gNiXNtxMtWBPz+L0dXpw/PHBsgYt2jASvZgczUhhSm6kNw9kUS4U2J9b2K1Subg08pqoh2X9yiQvbySeBY/EkKG+5rZwSKKfk1ENr0i3+8cllihWjPU/ZUp815ZQ0TTCUTbHczv3XdDX+Y/cEf28beKQTo+SSUxnqxChtuVht9MfvC/7eNvDQ8QUu3DYSvBQCTIds9dofzKQoVoxgXeItbL9CdchZhwJ6vwi/UXJBf0ZNzLZrlVJ9a2KUylWDR08scunOsfY+oLiavgG1MIsVI9jpuxbW0CgfmsmRL1Xbo2/APaeUTZFrJ6cEsP0Zqvv0GogdHjxq0jftREqVojq2xGH8Q9kUS4U2I93tV8DJB9fEKD9yYpGJ4SwT7eTTWtB3bRnlzICKFtZg7huG5JHji1zcbpRcWl41dqjn1Nqa/zueoR6P7mvvmtYA0TVKndB3elp5TGvgLT12aolSxehgYq6mbwCG+zrwFsfPUi1H1sBbfuh4B54yuC/MjN4+hbHjmerx2N3tvT8AHjy2wFC71K2VU8iu/u6G+zqgcNbQW37g6AIXtxslu9B3Q+bc7yhSPn5vz8UOB2dyLBcrXLy93bXvzhIALLcz/g1nqua0azD320VkjdJgxuLVO1iYJ+7vudjB8pTbzymspm8AhsyD8tryFoVQ418Db/GeQ3P0p/X2RB7gujAHs23WagCM7VYKxDUwyg8cW+DidqnboiVycbr/HRgly1vu8fgX8mWemFzmmWdsaO8DWkRKHdHX+WnVnLaHsByyth1Sh2a0UM+ptZ1T3fHMNaNv20FkjZJ1emlHE7O01PMGhQ8eW2C4LxW8vZCF0mpJMNgWZrEDCmfykZ53jL770ByX7xojFbSThYVy3p2+K1UxjDa83drC7K23aFG3bVF3YKtRc4qU0yy2O/ctb7nHTsk9R+YAeMbuNo2SJYluipSyKY2UJtpXoK0RfX3v4Tn60lp79Wlg1uh5REqdMAVTj9bLTUKGyBqlkX51YxZX2tyUd16pHo/c2aUrcsa+g2pTbkt5Bq6R0khfh97ijmeofMXJB9p7vw/kihUeObHYvqdsGB4URgdCFzAX5oF6NNIDPHhsgVLF4PJdbY7f2jQc7v9wNsVyuzklK1LusVG+99AcmoDLdo219wEFxTLQ10j/CSE6k0VvuRj0bM/X/l0HZ7li14bgrcUseJRDQJuREqi5L401U2AGRWSN0mi/ipQW2jVKm85VFM6h/+riVTViNlfiwKklrjpzvP0PKS2vKpyFOq/eFn0HsOs56rGH47//6DxVQ/LMPZ14yqvbzEAXFubOZ6rPPnpXe+/3gZ88PQvAs9u9/7VIyVno0rZDArD7uaosID/b/me0wN2H57hg20jtXgXGyryKFGwHHFoYyqZYavfep7LKKe3h3F8qlHnk+CLP6mTtF5ccRU51oUOHOdUeG+V2EVmjZNF3i514i2c8r6cT866DasFfddbG9j+kuOScU+g0hB+aUC1nDv6o/WtrgbsPmvRNu5FCYV499o2t+tVgpkOjtOsqVf/19A/be78P3PnUDGdPDLanPAPPnFJHQgeAPc8HJBz+cfuf4YGqIbnv8Hz71B2o+9/v/P7hvg4iJYAzrlaRQmGx/c/wwD2H5zEkPHtPB0ZpZc5x/PW13yZLMLip52u/E0TWKOmaYDibaj9SArUw5w/D/JHuXZgNdz41SzalcWk7NSqg1EEOxZNgz6l1OP7Dd/SsiHbfoTn2bh5idGC1p+sLK8qoeS/MNjem7LCShh/sjVGqGpJ9B+c6c0g8ckpDfSnypWp7tSqgvOVUX882pkeOL5IrVdunbkFFSk31iRY6KiAF2HO1orB6FC385OkZdE1wxe6x9j7AMFyNUscsAcCZL1AOyRqUBQRFZI0SwEh/msWVDr0l6Fm0dOfTM1yxe6x2/k9glFdAVh1DeEsS3jaFAcoolZZ7wi0XK1V+8vQszz27g03ZopYcFqZl6ObzHRjUM1+oVEg9yCs9emKRpWKlM+rW8uId778af9u1Wqks7Hp2z4zyD59QXdifd04H939lzjFKBtovHraw89mgpXtmlO96eo6Ld4y2T12WlpTRdJj7I30dpi5Azf1yPpTS8EgbpeG+DiOlzReqSd+DhbmwUuaRE4s8+8xOFqW5KQ+s/gxLgdTRwrSMcg/Gv+/gHCvlKi/cO9H+h1iR0sDqjX3DQAaAuXwHnt6ZL1BG/1D3Kaw7npoBOsgnAeRnFMXYtzrSHu5UfQmw5wWq3VIP8ko/fGyaC7aNsHm4r/WL3VCYd42URvrTna39zICKFntglJaLFe49MsdzzuqQugNHo5RJaQxnU8zmOpj7Z1wNCHj6B+1/Ro8QaaM02p9uP6cE6iTSs14Mj3+364V0P3hsCinhReduav9DctPq0cEoCSEY7kt1Nv6hCWWYn/he+5/hgh88NkVaF51FSh4Lc3zQNEqdLMxdV4Gegadva/8zXHDLgUn2bh5i22gbTWgt5GfUvXdQblqRckdMwZkvBGTXx58vVdh3aJYX7u1g7oOi71wipY2DGWY6ufeg1v7xeyA309nnNOFHj09Trkpecl4H5596zH2ADYOZzhyygXHYeklP5n6niLRRUvRdh/mQc38alk92ncL6/v5JNgyk25cDg9qUQCUmHTA+mGEu1+n4r1Hc8sp8Z5/ThNsem+KZZ2xon74Az4U52p9GCJjthL5L9yuxy2M3tf8ZDlgqlLnzqVleekGHhzJbRskBGwa7ECnuuFJ1TD/wnfY/wwF3PjVLuSp54bkdRMngGSmND2ZYKlTaa7Nl4dxrFEX2xM3tf4YDbj0wyXBfqsN8Wmuj1FGkBHD2S1VOzfpbIUGkjdJoN4zS3lcAAh7r3sKsVA1uOTDJS87b3F4jRgt5d/oOYONQlunlYvufD3Deq1S90hPf7exzbJhcLLD/5FLnm9LKrErGp1dHG7omGO1PM9/Jpgxw/s/AzOMw9Vhnn2PDDx+fpmJIXnb+ls4+yMMobRpSir6O7r+eUhvz4zd1tWv0rQcm6UtrnW3KlaLKebgYJcsod3T/t10OQ1vhwLfb/4wmSCm55cAkL9w70X59ErQ0SuMD6c4cElBz36jA4901yp0i0kZppC/NfKdGaXCjSvh20Sjde2Se+Xy5C56yO30HsGko07lR2vFMVa/VxYV508MnAXj5BR1uyi7qIwvjA13wFs97pXo88B+dfY4N3330FGMDaZ7RrvLKgqdRUpvy9HIXxr8y1zUVmmFIvv3QSV587ubgJ83aYUXuHvQd0BmFp2nKKD/xva6p0B4+vsipxSIvPq9Th8wHfdcpS7LjmTC0BfZ/s7PP6TIibZQmhrPkS9XOpJGgJubxe7vWC+s/Hz5JShOdRwr5GRCax8LMds6ra7qiMB+/uWsL8xsPnODcLUPtnbZpx8q8opdcMNYNb3F0p/KY93fHKJWrBrfsn+TF506031rJgodRGulLk9IEM506JWe/VOXVDnyrs88xse/QHJNLRV516bbOPmjFXXkJ9ZxiV5yS0lLXxD7fevAEuiZ4yfkdOqQ1ozTm+OvxgQ5zSqCM8nmvUka5x+3GgiDSRsk6XnpyqcOFedHPqscHv9ThFan6lK/ff5wXnzdRk262jfyM2pQ159u0cSjDfL5Mud1aFQsXvhqKC13h1k8tFrjr4CzXXrK9489S4/eIlLrhLYKiMY7e1ZV6tdsOTDGXL3PdZR2O36iqjcnFKGmaYHwww0ynkVJ2WBmmh76s/maH+I8HjpNNabys00156YR6HHY2bl2JlADOeglkR7uy9g1D8rX7jvP8czbV6NW2kZtWnTxSzp+zYTBDvlSl0G6bLQvn/4wqC3my+2KndhFpo7RlRMlNJxc7tPLjZyol1v3/1rEK7/Ynpzm1WOT6K3Z2dk2gJqbLpgQqpwQdKtBAbUqDE3D/5zv7HJSnKCVce+nWjj+LxeMw4u5xb+gGfQdw6evV4/1f6PijvnrvMcYHM11I8i+oJLzH/d80lGUm16FDBnDZG5UR6FCJVakairo7b6IzgQvAommUXO5/LVLqNFJM98FFr4VHvt5xg9J7Ds9xbH6F117RBYds8RiMuH9O1yLFs14Eg5vhvs919jldRKSN0uYRtSmf6jRSArj0DapzbocNSr967zGG+1K8rNN8EiijNOi+uW0a7FJeQU/DJT+vVFgd1KxIKfnivqNcuG2EczZ3SN0ZhtooPRbm5hEl9Ki20yncjg17VM3Off/akVOyWChz86OnuO7SbZ0luQGWVF7OTXkJKlKe6vTeA5z7SlULdV9nTsmtB6aYXCpy/RU7Or+mFpHS2EBGqS+74ZRc9iZ1dlOHuZV/v+8YfWmNn7qwWw6Z+9yfMB3SU5065HoaLnuDyqkvT3X2WV1CpI3SluEuRUoAF12vuPV7PtP2RyyslPnOQye59pJtnSV5ax94ROU8XLCxGwosC5e9EYxyRzTGPYfnefTEIm9+zu7Oryc/o06HHXHf4LaP9VMxJJNLXbj/l78Z5p6GQ7e3/RFfu+84pYrB9c/oQpS8YFKJY+7f5aahLNPdcMjSfYrCfvQbHZUG/Oudh9g8nOVlnQpcQBmlvjFH5SUo9eWW4T6OL3Th3u9+DoydAfd+tu2PWCqU+fd7j/PKi7fVWmB1hMUTnnN/xwb1vRyf78bcf4tS4T3wb51/VhcQaaM00p8im9I6zymBKia7+HUqjG1Tt/+FnxwmX6ry1uee0fn1VMsqhB/b5fqS7WPKKB+bX+n87227VKlx7vxE27mFf73jEEPZFK+9vAue8qIpOvHwFneMqYV5bK4L47/w1Sq3cNf/bevthiH5px89zaU7R7ms3V6HdswfVo+j7vd/22gfpxYL7fe/s+PKX4LKCtzzz229/chsnlsfm+KNz9rVeZQI5qbsTYPt3NDPkdl8539LCHjWO5XY4UR7TMmNdx9luVjh7c/b0/n1VCuqdtKDuraM0rH5Lox/8/mq7dJd/68recVOEWmjJIRg80iWk93wlgCe+xuqNuLu4AuzUjX459sP8pyzxrmo3eOP7Vg8rnIKHp7yttF+0rrg0EwXJibAc98Ns0+1JQ+fXCzwzQdPcP0VOzrPJ4AaP3h7i5ZR6oZRzgzCs94Bj3xNfQcBcetjkzw1neOdzz+z/bOz7Fg4oiL3Ifeo44yNA1QMyYluzP9tlykK845PtKXC/PTtB9GE4A3P7kKUDCpSGvamwXZu6OdoNxwSgGe8TR2TccfHA7/VMCT/fPtBrtg91v7ZUXYsn1Jr38Moj/SlGc6muuOQATzv3YopePQb3fm8DhBpowSwe3yAQ93wlkC13djzArjzHwJLJL/5wAmOLxR45/PP6s61WPSNh6esa4KdGwa64y0CXPBqGN0Nt/9N4NzKJ257iqoh+eUXnNmda/FhlLZ30ygBXPVrqtfc7R8L9DYpJZ+47Sm2jvTxqks6lEJbmDepWxflJcDucXXOVNeckuf9Niwdh4e/EuhtU0tF/vXOQ7z28h01R6FjzB/2pK4Bdo0PcLJbkWL/GFzxFnjwxsClId966AQHZ/L80tXdmvvH1KPH3AcVLXVt7p//MzB+FvzXX3e95VpQRN4onT0xxFOTy8hufZEvep9amAFonFLF4H/ffIALto10LoW1MN86pwBqYR7ullHSU3D1b6tCysf/0/fbJhcL/Oudh/jZK3ZwxsbVB/K1hZknID3gKfQYzKYYG0h3z1sc3qqS3vf+C8wd9P22Wx+b4idPz/LrLz67O9QVmJuyu0MCsHvjAED37v85L4ctl8AtH1IdFXzi//7wKUoVg3e/9JzuXEd+VhWObzrX82U7N/RT7VakCCpaEEKN3ycqVYOP/OdjnLtliGu75ZBMm91FNnp/nzvGuhgparpySo7f09VGAm1dyrr+9S7grE2DLBUrTHUj2Q+qSeU5L4cf/H++k76fu/MQR2ZX+O+vPB+tk7ZCdsw8AUJv6S2eMT7AoZlcd/4mwDPfDuNnw81/4rv1zEdufoyqIbu3KQFMPgoT53tGCqDu/+OTnUl5G/Di/66ipe/+ma+XG4bkf33nALvHB3hTt6grKVXro3HvqHvrSB8ZXePQbJfuv6bBT/0ZzB+Cuz7l6y2HZ/J8+vaDvPbyHZy5qUsOydQB9bjpPM+X7RpXRvnp6S6Nf2w3PPtdKq986mFfb7nx7qM8NZ3j919xXmctxeyYfFQd175hj+fLzpoY5KnpXHciRVCR4sa98J9/3LMz1vwg8kbp7M3qALSnprq4Mb/8A6pO5Ht/3vKlk4sFPnLzYzzv7I2dd0W24+SDalN2KZ6zcNbEIIuFSufSUAt6Wo1/ar9KfLbA3Ydm+cJdR3jn88/sXpQEamFuvqDlyy7YNsKjJxa7FymPbFce48NfgYOtz9n6lzsO8eiJRf7gmvPIpLq0nBaOqPm39RLPl+ma4MxNgzx2sovnQZ3zMlVQettf1mXpHvizbzxMWhP8t1ee371rmDaN0oR3pHTB1hEAHjnRxdNjX/D70DcC33qvKkvwwPRykb/8zn6uPGMDr7iwC4pDC1MHVJSoeSt4L9g2QqlidM8o62l4xV8oh/jOf+jOZ7aByBulc0yjtL+bE3PrJfDc34R9n/I8LltKyfu/9hDFisFfvPbi7iS4LZx8oOWmBNROtX3g6EL3/vYF16lGtd/7M8+kf6Fc5Q+/8hDbRvv47Zft7d7fz81AbtK3UVoqVLrHrYOiMDfsga/9hmdB5ZHZPH/1nf288NwJruu0rY4dJx9Sjz7u/yU7R3nw2EL3jDLAqz6s6Ltv/K5nfuHr9x/ne/sn+d2Xn1srZO8KJvdDql/lNz2wYTDDjrF+Hj7exbU/MK425kP/pda/Bz74zUfIFSv8z5+9pLtrf2q/UsS1wAXbemCUz70G9l4D3/9gV5sUB0HkjdLWkT62jfZx9+H57n7wS/4INpwJX/21+rlGTfjCXUe46eFTvOenzuWsidVHVreN5SmlPvKxKV24bRRdEzx4dL57f18IuO6v1cmcX3mXa37hg998hAOnlvjQz17SHcWdheP3qEcf479gmyrS7erGlBmE1/49zB2C7/w3x425VDH43X+7DwHd35ROPgAIddZVC1y6c5Tp5VL38ioAm/bCS/8YHvs23P1Pji85PJPnD7/yIM88YwNvv3pP9/42wJE7VHlCC+oW4KLtIzx8vIsOGcAVb1VdTm7+E1ca74v7jvC1+47z7pfsZW+nPR7tWDyuIuVtl7V86dkTQ2R0jUe6OfeFgFf/jaoP++q71qUnXuSNkhCCZ56xgbuenu2ut5gZgNf/E+Sm4EtvX8Wx3n1ojj/92sO88NwJfuUFXVLcWThknoa581ktX9qf0Tl3yzB3H+7ymSgj29XkPHoX/MfvrdqYv3jXEf71zsP86gvP6uwwMyc8/QMlh9757JYvvWj7KH1pjR8/2d2D2jjjeYrKufezq6gMKSUf+MbD3H1ojr963aXdU5xZePqHsOViyLZ2dC7bOQao+dhVPOc3VG71W+9bVVC8WCjzrn/Zhybgr994effEHaCOpj9xv/r+feCK3Rt4airXnQJ6C0LAaz4O2RH4wi+s6nLywNF53v/vD3H1ORv5zZec3b2/C3XKeM8LWr40k9K4fNcY//Wks9PcNoa3wqv/VjWp/sbvrLkaL/JGCeCqM8c5uVjgyW7mlQC2X6EihoM/hC//ci3x/8jxRd7x6bvYNtbH37zx8u4lOC088V3V9mXHM329/IXnbuInT892dgqtEy56LbzwvWpj/t6f1ybnzY+c4oavPsgL9m7iD67xTka3hadvUwY5M9DypX1pnavO3MgPH+9Bi5SX/BGcdy3cdIOSCpv4Pzc/xufuPMyvvehsfubSLvQ5s6O4pCKFvS/39fKLd4wyPpjhe4+e6u51aDr83Kdgwxnw+TfB8fsAWClV+dXP3M0Tk8t87Beewc4Nre9RIBz8karR8WmUXmT2GLz1sS7f/5Ft8IbPqsjlsz+ncnzAY6eWeNs//oRNQ1n+5o1XdN4JvhlP36rWvg+WANTaf+jYYufd4ptxwXVq/j/wBRUxmmv/gaPz3XX+HRALo/SKi7YiBHzzgePd//DL36Q45kf+Hb70Nu5+/Ahv/n93MJjR+ew7r2JsINPdv1cpqZNQz36pkmj7wE9dsIVyVXLL/snuXgvAi/9QKfJ+9BH4zg38+90H+fXP3s3F20f4xFue2V0vGWD6CeUp732F/0s8b4Inp3I8dqqLCX9Q9NHPfhJ2Pxe+8isYd/0T//Nbj/I333+CNz5rF+/rhUE+8B3V8uWcn/L1cl0TvPT8zXx//yTFSper8fvH4C1fURHDZ17D8v5beeun7uSOp2f4/15/WedNZ53w4I2qM/wZV/t6+QXbhtk60sdND7UWZQTGrmfBz/+LEh3986t54JFHeeMn7yCta3zuV66qtfnqGsoFVby695qWIgcLLzZZim/3YvwvfC9c+U64/W+Q33ovf/e9R3n1x/6Lbz5wovt/y4ZYGKUtI30858yNfGnf0c6OR3bD834L45q/xNj/LQY/+you7pvi8+96Tk2S2lU8+nVFGV7+Ft9vuWL3BnaN9/MvPz7U/evRNLj2/1B51rvgzr9n47+/mZft0vjsL1/V3TyShXs+raTwl73R91tec/kOMrrGZ+/owfizQ/DmL1E640Vo//G7nHH7DfzSs7fyoesv6Z783467/p/KZe5+ru+3vPbyHSwWKnztvh44ZRvOgLd/k0J2nL4v/CyXH/88f/emy3ltN5quNmN5Up1rddH1kPLn7AkheN0zd3LLgcnuFZHbcd5PI9/wWcqTj7H1336aq9MH+MK7ntNdpamFR/5dRWSX/4Lvt1y0fYQLt43w2TsOdT+CEQKu/d8UnvUbiLv+L8+87Zf4xUv6eMVFXVQaOsCXURJCXCiE+J4QIi+EOC6E+HMhRBc6jnYP73rRWRybX+Fzd3Z/Y9p/cpE3P3gFbyu9j536HJ8pvoczHv44lLuo+AIlKLjlQ6pW4OyX+n6brgneefWZ7Ds015No6QdPzPCqx36G95bfxXNTB/jEwq8y/PDnut8na/4I/OT/wcU/27LFjB3jgxlee8V2vvCTIxzsljzWRKVq8MUHZnnuoV/lH4zX8AupW/jTI+9Ee/K73efaD3xHUXdX/ZqvJL+Fq8/ZyAXbRvjb7z/OSqm792QhX+aD/5Xjqqk/4g7tGfyx/hleddc76grBbuK2v1JNeJ/zm4He9ubn7Cala3z4pgNdv6T9Jxd56w83cG3+T5HpAf628EecdccfKYVoN1HKq7W/5RI480W+3yaE4Jeu3sP+k0tdd0oMQ/L1B07wovtfxu9XfoNnpp7iz46+k+zh7hyI6IaWM18IsQH4LiCB1wB/Dvw+4K+6cI3w4nMneMHeTfzld/Zz96H2j1+wIKXkzqdmePfn7uGVf/1DHjq2wLWvfTOD79mHOO+V8P2/gI9eqtpydKPlu2HAt/4AZp+En/7LQJsSwJuu2s3ezUO878sPcLgLbWcqVYPvPHSSN37yx/ziP/6EQtngVb/4XlK/cTtiy0Xwjd+Gjz0L7v40FLqg/ikuw43vUCftvvSPA7/9919xHtmUxrs/fw8LK53n1paLFT5352Gu+egPeN+ND7B7YoSXvfvj8ItfU9TKv74OPvUKRbd048Te2afg6+9WBaNXviPQW4UQ/Ol1F6oC7q880PlRHqjWTX/57f288MO38I//9TSvuvJcLn3vt5QAYPoAfOJq+Le3qMR8i3oeX3joyypKfNYvw6ZgRdjbRvv59RedzdfvP86/dCFallLyk6dn+a3P38sr//qH3H90nre+5pVsfu9dqj/kPZ+Bj14CN/2Rops7hVFV62n+MFzzPwKv/Z99xk4u2znK+7/2UFeUePlShS/uO8K1f/sjfvvz97JxMMsv/fp/J/3rP0Jsv6Jlp4lOIVqFfEKIG4D3AWdIKRfN594HfADYaj3nhiuvvFLu27evO1fbApOLBd7wyTs4sbDCe685n1949m76M/4Dunypwt2H5rjtwBTf2z/J09M5hrMp3vycM/i1F53VmD86+F+qwPDpHyi66ZyXwdkvU4dmTZyvQl+/mDoAN/2hEji84A/gZe8PMOo6Hj+1xOs+8WOEgD++9kJec/n2QDmfqaUi9x6e45YDk9z8yCTTy0V2jPXzS1fv4a3PPYNsyvwuDUOdPfODDyv5cqofzvtpFd2d+UJ1DIDf8UupvsNvvw+mH1eKxwtf08bo4fv7T/Guz9zNrvEB/vS6C3nh3gnfFJuUkqNzK/z4qRluOzDFLQcmyZeqnL91mN99+V6uuWhrXfZdKaqN6b/+BhYOqxzIBdfBWS+GM54PwwHojUpJHRdy8/vV5vTOm1sWjbrh7255gg/fdIBnnznOn/zMhVy8w39j4FLF4NETi9z59Aw3P3KKfYfmEMA1F23l3S89p7HJcH4W7vh7pUosLqhOCBe8Wt373c9Vxad+sXRK9Vq84+NKbfmLX1NHaQREpWrwK5/Zxy0HpnjTs3fzOy/by9ZR/5+zsFLmgaPzjWu/L8UvXLWbX3/R2Y1rf+qA6vjy0I1KlLHjSlXfc+aLYMczVBGqH0ipyh/+8/2qLuplf6IUn23g6Fyen/v721kqVHjfNefxxmfv9n18jjX37zo4yy0Hprh1/yRLxQrnbB7i3S85h1dftr3rVLUQ4m4p5ZWOv/NhlH4AHJdSvtH23G7gEPBqKaVnW9m1NEqgqqx/74v384PHphjI6Dz3LEVtbB3tY6Q/TUoTlKsGxYrBbK7EyYUCJxZWeOzUMgdnckippJbP3jPOay7fzrWXbmMg45E7mTqg2pI88u/1fmmZIVX4uelcRUMNbVFnw2i6+imvqNqn2afg6D449aA6+vjlf6o8xQ5qXp6ezvE7X7iXB44uMNqf5nlnb2Tv5iG2jvYzkNHJpjQKlSq5YpXFQpmjcyscmc3z1FSuVoA6lE3xovMmuO7S7bz8gs3uCiMp1fXf/3kVMeRM6rB/XBnmTXvV2Ic2K0WR0NTYisvqtdOPw6Efq419ZKeSoJ/zsrbHDnDXwVl+9wv3cWx+hS0jWZ5z1kbOmRhi03C2lgMrlqsUKgZTiwWOLxQ4Opfn0RNLtQhLnQm0mZ+/cheX7xpzr0GqluHJW9Q5NI/dBCVTaDG8TY194151//s31I91lxIK86pbwtR+tRmtzKmN7fpPqPd1gBvvPsqfff1hlooVzpoY5MozNnDGxkE2DGToS2vomiBXrJIvqdZcR2bzHJld4cCppVo+9sJtI7zioi28/spd3nL3Ug4e/aYa/8EfQbUICGWkNl+gHgcn1E8qq36HVDnTxRNK0HLkTpBVFR2+4i9UjVibKFcN/urb+/mn2w9iSMnlu8a4ePsou8cHGO5L0Z/RKVUMChWDpUKZE/MFjs+v8OTUMgdNdiGja1x11jjXXbqd6y7b7u3ULp6AB7+oorwTD6ixaWl1DyfOU/NgcEKdHqzpav5Ximr8cwfV2GeeUGvjmv8JV7y57bGDOvDvPf92H7c/OUN/Wueqs8Y5d8sw20b7GMymSOuCQtmgUK4ymytxbH6FY3Pq3s/n1dzfNJTlJedN8Pord/GsPRu6W39nQ6dGaRL4uJTyA03P54APSCk/7PX+tTZKoCz/XQfn+Oq9x7j70CxPTuVcKY2hbIotI1n2bh7mgm0jXLprlOecuTFQhFXD3CHl9Z98ECYfURNueVItOidkR2HHFSrCuOwXYKg7aibDkNz62CTfevAkP3l6lqNzedwYnQ0DaXaND3DGxkEu2znK5bvGuGTnaD0q8gspVSNJa/xTBxQVmZtGMb8OGNoKu56tvMyLftaXBNwPShWD/3jwODc9dIoHjy24dnsQQhmg7WP9nL91hIt3jHDFrg1csG04+GKsVlTUePBHqkXS9GPq/hfmnV8vNNU1YtdVauznvDwwbeOGhZUyX73nKLccmOLh44uuh0BmdI2dG/rZNT7AuVuGuHzXBp5xxhjbRtuouyoX4OhP4PAdavxTB9SZWAWXwtb0gCoOPvMFStQTkLLzwqGZHP9+73F+8PgUj51aYqng3MNxtD/N9rF+zhgf4JKdo1y8Y5Qrz9jQnoAnP6tKR47fW7//S6fUibZOGJyAbZfDea+Ei39OKR27AIt6/OYDJ7jr4CxPTeccxV+aMBsPjPWzd/MQl+wc5bKdY1y4baQ3Ap4mdGqUysB7pZQfbXr+KPAZKeUfer1/PYxSM0oVg/l8iYWVMlUpSesaGV1jw2CmO6dEesEwYGVWLU6jquS+qayalNnhjqIivyhVDGZyRQplg2KlSjalM5jRGepLeUeB3YBRVafIFhYBqf6fGVTjb4OmaQcrpSrzKyWWzc2pL62TTWtsGMh0X9LejEpJRULW5iyEuu+DE75lv50iV6ywVKiwUq5SNSRD2RSDWZ3BTKr3G1ClqByTagnlnAgVOazR3JdSsliokCuq8Wd0jb60zmBW7/3cBxVN5mcVzWdFUoMTvtWFnaJqSBZXyiwXK1QMSV9aoy+lM9yX6n6NVQB4GSW/d8XJcgmX5xFCvAt4F8Du3V3qnNwBMimNzSN9bO5mfy6/0DQY3KR+1gmZlNae99sNaLqi74a63PUhAPozOv2ZfujC2YuBkcqoHFOQPFOXMZhN9Ua+7wepLIz2QD7uE0IIRvvTjPb7zPN0G5nBjijJTqFrgg2DGTYMro0R7Ab8mMo5YMzh+VFg3ukNUspPSimvlFJeOTHRgwK7BAkSJEgQS/gxSvuBhpa1QohdwKD5uwQJEiRIkKAr8GOUvg1cI4Swt8J9A7AC3NaTq0qQIEGCBKcl/BilTwBF4CtCiJeb+aIPAB9pVaOUIEGCBAkSBEHL7KeUck4I8TLgY8A3UHmk/4MyTAkSJEiQIEHX4EuSI6V8BPDfjC1BggQJEiRoA7HoEp4gQYIECeKBxCglSJAgQYLQIDFKCRIkSJAgNEiMUoIECRIkCA1a9r7r+A8IMYXqKN4JNgHTXbicOCL5btyRfDfeSL4fdyTfjTu68d2cIaV0bPfTc6PUDQgh9rk17zvdkXw37ki+G28k3487ku/GHb3+bhL6LkGCBAkShAaJUUqQIEGCBKFBVIzSJ9f7AkKM5LtxR/LdeCP5ftyRfDfu6Ol3E4mcUoIECRIkOD0QlUgpQYIECRKcBgi1URJCXCiE+J4QIi+EOC6E+HMhxNqcIR0iCCHOEUL8gxDifiFEVQhxq8NrhBDiD4UQR4QQK0KIHwghLl/7q107CCFeL4T4uhDimBBiWQhxtxDiTU2vOe2+FwAhxOuEELcLIWaEEAUhxAEhxB8LITK215yW300zhBA7zPkjhRBDtudPu+9HCPF283to/vk122t6+r2E1igJITYA30Uduf4a4M+B3wf+bD2va51wEfAq4DHzxwn/HXg/8FfAdcAy8F0hxNY1ucL1we+hxvke4NXALcDnhBC/ZXvN6fi9AGxEfR+/DLwS+Efgj4CP2F5zun43zfgwauzNOJ2/n5cCz7X9fMX2u95+L1LKUP4AN6COYh+xPfc+IG9/7nT4ATTbv28Ebm36fR+wAPyJ7blBYAr4i/W+/h5+L5scnvsc8PTp/L14fF//A3X0jEi+m9qYXwDMAn+AcoCHTue5A7zd/j04/L7n30toIyWUd3eTbDxI8AtAP/Ci9bmk9YGU0mjxkucBI8AXbe/Joc6/emUPL21dIaV0qiq/F9hs/vu0/F48MANY9N1p/92YqYC/RbEwzXPptP9+XNDz7yXMRul8YL/9CSnlYVSkdP66XFF4cT5QBR5vev5RTr/v6nnAI+a/T/vvRQihCyEGhBDPB34b+Hup3NvT/rsBfg3l+f+dw+9O9+/nSSFExcxF/qrt+Z5/L74O+VsnbEBRDc2YM3+XoI4NwLKUstr0/BwwIITISClL63BdawrzhOTXAO8wn0q+F8gBWfPfnwHea/77tP5uhBAbgQ8Cb5FSloUQzS85Xb+fE6h80U8AHXgT8AkhxICU8v+wBt9LmI0SKG6zGcLl+dMdbt+V2+9iBSHEHlQ+6WtSyk/bfnVafy+oyHEAeDbwJ8DHgN8wf3c6fzf/A7hTSvktj9ecdt+PlPIm4CbbU98WQmSBPxZC/LX1Moe3du17CbNRmgPGHJ4fxTmCOp0xBwwLIfQmD2YMyEspy+tzWWsDIcQ48G3gMPAW269O6+8FQEp5j/nPHwkhpoF/FkL8b07j70YIcREqmn6hEGLMfHrAfBwVQlQ5jb8fB9wI/DywhzX4XsKcU9pPE0cphNiFUnrsd3zH6Yv9qFD7nKbnV+Xl4gYhxADwTVQC/1oz6WrhtP1eXGAZqDM5vb+bvUAa+DFqk52jnlc6ihI/nM7fjxska/C9hNkofRu4RggxbHvuDcAKcNv6XFJocTuwCLzeesLcrK9DfY+xhBAiBXwJtcm8Uko52fSS0/J78cDV5uPTnN7fzY+AlzT9/JX5u1eh6pZO5++nGT+HUiceYg2+lzDTd59AqYW+IoT4K+As4APAR5pk4rGHedNfZf53BzAihHid+f9vSSnzQoi/BN4vhJhDeSy/h3I6/nbNL3jt8HHU9/I7wLgQ4jm2390rpSycpt8LQojvoIrPH0appa5GFZ//m5TySfM1p+V3Y5YS3Gp/zsxJAvxQSrlsPnfafT9CiC+jRA4PoCKiN5g/v22WpvR+Ta13sVaLQq4Lge+joqMTKLWMvt7XtQ7fwx5U6Oz0s8d8jUBV7B81v68fAles97X3+Hs5mHwvrt/NB4GHUNX28yjq7reAtO01p+V34/J9vZ2motHT8fsBPgQcQJXerAB3A29tek1Pv5ekS3iCBAkSJAgNwpxTSpAgQYIEpxkSo5QgQYIECUKDxCglSJAgQYLQIDFKCRIkSJAgNEiMUoIECRIkCA0So5QgQYIECUKDxCglSNBjCCF+Xgjx9vW+jgQJooCkTilBgh5DCHEj6pTcF6/3tSRIEHYkkVKCBAkSJAgNEqOUIEEPIYT4NKqh5YuEENL8+cD6XlWCBOFFmBuyJkgQB3wQ2I06b8Y6XO/oul1NggQhR2KUEiToIaSUTwohZgFNSnnHel9PggRhR0LfJUiQIEGC0CAxSgkSJEiQIDRIjFKCBAkSJAgNEqOUIEHvUQL61vsiEiSIAhKjlCBB77EfuEQI8VohxJVCiO3rfUEJEoQViVFKkKD3+Djwn8A/AncB71rfy0mQILxI2gwlSJAgQYLQIImUEiRIkCBBaJAYpQQJEiRIEBokRilBggQJEoQGiVFKkCBBggShQWKUEiRIkCBBaJAYpQQJEiRIEBokRilBggQJEoQGiVFKkCBBggShQWKUEiRIkCBBaPD/A+Ek/l6m//qNAAAAAElFTkSuQmCC\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "plt.rcParams.update({'font.size': 15})\n", "fig,ax = plt.subplots(1,1,figsize=[7,5])\n", @@ -96,30 +86,42 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "True\n" - ] - } - ], + "outputs": [], "source": [ - "# Check that numbalsoda and scipy match.\n", + "# Check that numbalsoda and scipy match for the hybrid LSODA solver.\n", "\n", "# numbalsoda\n", "usol, success = lsoda(funcptr, u0, t_eval, data, rtol = 1e-9, atol = 1e-30)\n", "\n", "# scipy\n", - "t_span = [min(t_eval),max(t_eval)]\n", + "t_span = [t0, tf]\n", "sol = solve_ivp(f_scipy, t_span, u0, t_eval = t_eval,\\\n", " rtol = 1e-9, atol = 1e-30, method='LSODA')\n", "print(np.all(np.isclose(sol.y.T,usol)))" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Check that numbalsoda and scipy match for the RK45 explicit solver.\n", + "\n", + "# numbalsoda\n", + "tsol, usol, _, success = rk45(\n", + " funcptr, u0, -1.0, t0, tf, 200, data, rtol = 1e-9, atol = 1e-30\n", + ")\n", + "\n", + "# scipy\n", + "t_span = [t0, tf]\n", + "sol = solve_ivp(f_scipy, t_span, u0, t_eval=tsol,\\\n", + " rtol = 1e-9, atol = 1e-30, method='RK45')\n", + "print(np.all(np.isclose(sol.y.T,usol)))" + ] + }, { "cell_type": "code", "execution_count": null, @@ -129,26 +131,28 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "152 µs ± 256 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)\n", - "11.1 ms ± 55 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n", - "\n", - "scipy took 73 times longer than numbalsoda\n" - ] - } - ], + "outputs": [], "source": [ - "t_nb = %timeit -o usol, success = lsoda(funcptr, u0, t_eval, data, rtol = 1e-6, atol = 1e-30)\n", + "t_nb_lsoda = %timeit -o usol, success = lsoda(funcptr, u0, t_eval, data, rtol = 1e-6, atol = 1e-30)\n", "t_sp = %timeit -o sol = solve_ivp(f_scipy, t_span, u0, t_eval = t_eval,\\\n", " rtol = 1e-6, atol = 1e-30, method='LSODA')\n", "\n", - "print(\"\\nscipy took \"+'%i'%(t_sp.average/t_nb.average)+\" times longer than numbalsoda\")" + "print(\"\\nscipy took \"+'%i'%(t_sp.average/t_nb_lsoda.average)+\" times longer than numbalsoda (LSODA integration)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "t_nb_rk45 = %timeit -o tsol, usol, _, success = rk45(funcptr, u0, -1.0, t0, tf, 200, data, rtol = 1e-6, atol = 1e-30)\n", + "t_sp = %timeit -o sol = solve_ivp(f_scipy, t_span, u0, t_eval = tsol,\\\n", + " rtol = 1e-6, atol = 1e-30, method='RK45')\n", + "\n", + "print(\"\\nscipy took \"+'%i'%(t_sp.average/t_nb_rk45.average)+\" times longer than numbalsoda (RK45 integration)\")" ] }, { @@ -160,26 +164,9 @@ }, { "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[5. , 0.8 ],\n", - " [5.02801788, 0.97855571],\n", - " [5.00657991, 1.19693864],\n", - " ...,\n", - " [0.1390893 , 0.1149374 ],\n", - " [0.1454087 , 0.1101119 ],\n", - " [0.15204821, 0.10552706]])" - ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - } - ], + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "# numbalsoda within jit compiled function works\n", "@njit\n", @@ -191,23 +178,9 @@ }, { "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "ename": "TypingError", - "evalue": "Failed in nopython mode pipeline (step: nopython frontend)\n\u001b[1mUntyped global name 'solve_ivp':\u001b[0m \u001b[1m\u001b[1mCannot determine Numba type of \u001b[0m\n\u001b[1m\nFile \"../../../../../../../var/folders/sf/43vm953d201c4jw4yg22hnhc0000gn/T/ipykernel_58698/4157945421.py\", line 4:\u001b[0m\n\u001b[1m\u001b[0m\n\u001b[0m", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mTypingError\u001b[0m Traceback (most recent call last)", - "Input \u001b[0;32mIn [9]\u001b[0m, in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 4\u001b[0m sol \u001b[38;5;241m=\u001b[39m solve_ivp(f_scipy, t_span, u0, t_eval \u001b[38;5;241m=\u001b[39m t_eval, method\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mLSODA\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m sol\n\u001b[0;32m----> 6\u001b[0m \u001b[43mtest_sp\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", - "File \u001b[0;32m~/miniforge3/lib/python3.9/site-packages/numba/core/dispatcher.py:468\u001b[0m, in \u001b[0;36m_DispatcherBase._compile_for_args\u001b[0;34m(self, *args, **kws)\u001b[0m\n\u001b[1;32m 464\u001b[0m msg \u001b[38;5;241m=\u001b[39m (\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mstr\u001b[39m(e)\u001b[38;5;241m.\u001b[39mrstrip()\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m \u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mThis error may have been caused \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 465\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mby the following argument(s):\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;132;01m{\u001b[39;00margs_str\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 466\u001b[0m e\u001b[38;5;241m.\u001b[39mpatch_message(msg)\n\u001b[0;32m--> 468\u001b[0m \u001b[43merror_rewrite\u001b[49m\u001b[43m(\u001b[49m\u001b[43me\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mtyping\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 469\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m errors\u001b[38;5;241m.\u001b[39mUnsupportedError \u001b[38;5;28;01mas\u001b[39;00m e:\n\u001b[1;32m 470\u001b[0m \u001b[38;5;66;03m# Something unsupported is present in the user code, add help info\u001b[39;00m\n\u001b[1;32m 471\u001b[0m error_rewrite(e, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124munsupported_error\u001b[39m\u001b[38;5;124m'\u001b[39m)\n", - "File \u001b[0;32m~/miniforge3/lib/python3.9/site-packages/numba/core/dispatcher.py:409\u001b[0m, in \u001b[0;36m_DispatcherBase._compile_for_args..error_rewrite\u001b[0;34m(e, issue_type)\u001b[0m\n\u001b[1;32m 407\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m e\n\u001b[1;32m 408\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 409\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m e\u001b[38;5;241m.\u001b[39mwith_traceback(\u001b[38;5;28;01mNone\u001b[39;00m)\n", - "\u001b[0;31mTypingError\u001b[0m: Failed in nopython mode pipeline (step: nopython frontend)\n\u001b[1mUntyped global name 'solve_ivp':\u001b[0m \u001b[1m\u001b[1mCannot determine Numba type of \u001b[0m\n\u001b[1m\nFile \"../../../../../../../var/folders/sf/43vm953d201c4jw4yg22hnhc0000gn/T/ipykernel_58698/4157945421.py\", line 4:\u001b[0m\n\u001b[1m\u001b[0m\n\u001b[0m" - ] - } - ], + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "# scipy within jit compiled function does not work\n", "@njit\n", @@ -227,7 +200,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3.10.8 64-bit", "language": "python", "name": "python3" }, @@ -241,7 +214,12 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.7" + "version": "3.10.8" + }, + "vscode": { + "interpreter": { + "hash": "b0fa6594d8f4cbf19f97940f81e996739fb7646882a419484c72d19e05852a7e" + } } }, "nbformat": 4, diff --git a/numbalsoda/driver_rk45.py b/numbalsoda/driver_rk45.py index 3f52963..390106b 100644 --- a/numbalsoda/driver_rk45.py +++ b/numbalsoda/driver_rk45.py @@ -22,24 +22,25 @@ rk45_wrapper = librk45.rk45_wrapper rk45_wrapper.argtypes = [ct.c_void_p, ct.c_int, ct.c_void_p, ct.c_void_p,\ ct.c_double, ct.c_double, ct.c_double, ct.c_int,\ - ct.c_void_p, ct.c_double, ct.c_double, ct.c_double, ct.c_void_p,\ - ct.c_void_p, ct.c_void_p] + ct.c_void_p, ct.c_void_p, ct.c_double, ct.c_double, ct.c_double,\ + ct.c_void_p, ct.c_void_p, ct.c_void_p] rk45_wrapper.restype = None @njit def rk45(funcptr, u0, dt0, t0, tf, itf, data = np.array([0.0], np.float64), rtol = 1.0e-3, atol = 1.0e-6, mxstep = 10000.0): neq = len(u0) - usol = np.empty((itf+1, neq),dtype=np.float64) + usol = np.empty((itf + 1, neq),dtype=np.float64) + tsol = np.empty(itf + 1, dtype=np.float64) success = np.array([999], np.int32) actual_final_time = np.array([-1.0], dtype=np.float64) actual_final_iteration = np.array([-1], dtype=np.int32) rk45_wrapper(funcptr, neq, u0.ctypes.data, data.ctypes.data, dt0, t0, tf, itf, - usol.ctypes.data, rtol, atol, mxstep, success.ctypes.data, + usol.ctypes.data, tsol.ctypes.data, rtol, atol, mxstep, success.ctypes.data, actual_final_time.ctypes.data, actual_final_iteration.ctypes.data) success_ = True if success != 1: success_ = False - return usol[:actual_final_iteration[0],:], actual_final_time[0], success_ + return tsol[:actual_final_iteration[0]], usol[:actual_final_iteration[0],:], actual_final_time[0], success_ @nb.extending.intrinsic def address_as_void_pointer(typingctx, src): diff --git a/src/RK45.cpp b/src/RK45.cpp index 00bba70..d9bd0ae 100644 --- a/src/RK45.cpp +++ b/src/RK45.cpp @@ -294,6 +294,7 @@ RK45::RK45(RK45_ODE_SYSTEM_TYPE rhs_handle, double t0, std::vector &y0, m_t_bound = tf; m_direction = (m_t_bound != t0) ? sgn(m_t_bound - t0) : 1; + m_rtol = _rtol; if (_rtol < 100 * std::numeric_limits::epsilon()) { std::cout << "Warning: '_rtol' is too small. Setting to max(rtol, 100*EPS)." << std::endl; m_rtol = std::max(_rtol, 100 * std::numeric_limits::epsilon()); diff --git a/src/RK45_wrapper.cpp b/src/RK45_wrapper.cpp index 6f1a3ab..fdd8641 100644 --- a/src/RK45_wrapper.cpp +++ b/src/RK45_wrapper.cpp @@ -31,6 +31,7 @@ extern "C" * @param tf Final time * @param itf Maximum number of iterations * @param usol Vector holding the solution at a given instant. + * @param tsol Vector holding the instants at which the solution has been evaluated. * @param rtol Relative tolerance for the integration. * @param atol Absolute tolerance for the integration. * @param mxstep Maximum allowed step size. @@ -46,6 +47,7 @@ void rk45_wrapper( double tf, int itf, double* usol, + double* tsol, double rtol, double atol, double mxstep, @@ -61,6 +63,11 @@ void rk45_wrapper( RK45 rk45(rhs, t0, y0, tf, mxstep, rtol, atol, dt0); + for (auto i = 0; i < neq; i++) { + usol[i + (0) * neq] = rk45.m_y[i]; + } + tsol[(0)] = rk45.m_t; + int itnum = 1; while (itnum < itf + 1) { bool successful_step = rk45.step(); @@ -71,8 +78,9 @@ void rk45_wrapper( } for (auto i = 0; i < neq; i++) { - usol[i + (itnum - 1) * neq] = rk45.m_y[i]; + usol[i + itnum * neq] = rk45.m_y[i]; } + tsol[itnum] = rk45.m_t; if (rk45.m_direction * (rk45.m_t - rk45.m_t_bound) >= 0) { break; diff --git a/tests/test_numbalsoda.py b/tests/test_numbalsoda.py index 0ac7b05..8239f41 100644 --- a/tests/test_numbalsoda.py +++ b/tests/test_numbalsoda.py @@ -36,13 +36,14 @@ def test_rk45(): dt0 = -1.0 t0 = 0.0 tf = 10.0 - itf = 1000000 + itf = 200 - usol, _, success = rk45( + _, usol, actual_tf, success = rk45( funcptr, u0, dt0, t0, tf, itf, rtol=1.e-8, atol=1.e-8 ) assert success + assert np.isclose(actual_tf, tf) assert np.isclose(usol[-1,0], 0.38583246250193476) assert np.isclose(usol[-1,1], 4.602012234037773) @@ -76,4 +77,4 @@ def test_solve_ivp_2(): test_dop853() test_rk45() test_solve_ivp_1() - test_solve_ivp_2() \ No newline at end of file + test_solve_ivp_2() From c8e85a36a9efdff7a9e66d602d5e725ff99a0289 Mon Sep 17 00:00:00 2001 From: tbridel Date: Tue, 6 Dec 2022 08:14:15 +0100 Subject: [PATCH 6/8] Fix passing rhs params to RK45 wrapper and stepping method. --- benchmark/benchmark_lorenz.py | 17 +++++++++++++---- benchmark/benchmark_robber.py | 2 +- src/RK45.cpp | 27 +++++++++++++++++---------- src/RK45.h | 6 +++--- src/RK45_wrapper.cpp | 4 ++-- 5 files changed, 36 insertions(+), 20 deletions(-) diff --git a/benchmark/benchmark_lorenz.py b/benchmark/benchmark_lorenz.py index fdaac63..4ffb0e2 100644 --- a/benchmark/benchmark_lorenz.py +++ b/benchmark/benchmark_lorenz.py @@ -1,5 +1,5 @@ import numpy as np -from numbalsoda import lsoda_sig, lsoda, dop853 +from numbalsoda import lsoda_sig, lsoda, dop853, rk45 from scipy.integrate import solve_ivp import timeit import numba as nb @@ -24,14 +24,16 @@ def f_sp(t, u, sigma, rho, beta): u0 = np.array([1.0,0.0,0.0]) tspan = np.array([0., 100.]) -t_eval = np.linspace(tspan[0], tspan[1], 1001) args = np.array([10.0,28.0,8./3.]) args_tuple = tuple(args) rtol = 1.0e-8 atol = 1.0e-8 +tsol_nb_rk, usol_nb_rk, tf_rk, success = rk45(funcptr, u0, -1.0, tspan[0], tspan[1], 30000, data=args, rtol=rtol, atol=atol) +assert np.isclose(tf_rk, tspan[1]) +t_eval = tsol_nb_rk +usol_nb, success = lsoda(funcptr, u0, t_eval, args, rtol=rtol, atol=atol) usol_sp = solve_ivp(f_sp, tspan, u0, t_eval = t_eval, args=args_tuple, rtol=rtol, atol=atol, method='LSODA') -usol_nb, success = lsoda(funcptr, u0, t_eval, args, rtol=rtol, atol=atol) @nb.njit(boundscheck=False) def time_nb(): @@ -41,6 +43,10 @@ def time_nb(): def time_dop853(): usol_nb, success = dop853(funcptr, u0, t_eval, args, rtol=rtol, atol=atol) +@nb.njit(boundscheck=False) +def time_rk45(): + _, _, _, _ = rk45(funcptr, u0, -1.0, tspan[0], tspan[1], 30000, data=args, rtol=rtol, atol=atol) + def time_sp_LSODA(): usol_sp = solve_ivp(f_sp, tspan, u0, t_eval = t_eval, args=args_tuple, rtol=rtol, atol=atol, method='LSODA') @@ -53,17 +59,20 @@ def time_sp_DOP853(): # time time_nb() time_dop853() +time_rk45() time_sp_LSODA() time_sp_RK45() time_sp_DOP853() -iters = 10 +iters = 100 t_nb = timeit.Timer(time_nb).timeit(number=iters)/iters*1e3 t_dop853 = timeit.Timer(time_dop853).timeit(number=iters)/iters*1e3 +t_rk45 = timeit.Timer(time_rk45).timeit(number=iters)/iters*1e3 t_sp_LSODA = timeit.Timer(time_sp_LSODA).timeit(number=iters)/iters*1e3 t_sp_RK45 = timeit.Timer(time_sp_RK45).timeit(number=iters)/iters*1e3 t_sp_DOP853 = timeit.Timer(time_sp_DOP853).timeit(number=iters)/iters*1e3 print("numbalsoda lsoda time =",'%.3f'%t_nb,'ms') print("numbalsoda dop853 time =",'%.3f'%t_dop853,'ms') +print("numbalsoda rk45 time =",'%.3f'%t_rk45,'ms') print("Scipy LSODA time =",'%.3f'%t_sp_LSODA,'ms') print("Scipy RK45 time =",'%.3f'%t_sp_RK45,'ms') print("Scipy DOP853 time =",'%.3f'%t_sp_DOP853,'ms') diff --git a/benchmark/benchmark_robber.py b/benchmark/benchmark_robber.py index 07f9b80..d6c2933 100644 --- a/benchmark/benchmark_robber.py +++ b/benchmark/benchmark_robber.py @@ -1,4 +1,4 @@ -from numbalsoda import lsoda_sig, lsoda +from numbalsoda import lsoda_sig, lsoda, rk45 import numba as nb import numpy as np from scipy.integrate import solve_ivp diff --git a/src/RK45.cpp b/src/RK45.cpp index d9bd0ae..35f449a 100644 --- a/src/RK45.cpp +++ b/src/RK45.cpp @@ -30,6 +30,7 @@ * algorithm is proportional to 'step_size ** (order + 1)'. * @param rtol Desired relative tolerance. * @param atol Desired absolute tolerance. + * @param rhs_args Extra parameters to pass to the RHS function. * @return double Absolute value of the suggested initial step. */ double select_initial_step( @@ -40,7 +41,8 @@ double select_initial_step( int direction, int order, double rtol, - double atol + double atol, + void *const rhs_args ) { int n_eqns = y0.size(); @@ -68,7 +70,7 @@ double select_initial_step( for (auto i = 0; i < n_eqns; i++) { y1[i] = y0[i] + h0 * direction * f0[i]; } - fun(t0 + h0 * direction, y1.data(), f1.data(), NULL); + fun(t0 + h0 * direction, y1.data(), f1.data(), rhs_args); double d2 = 0.0; for (auto i = 0; i < n_eqns; i++) { d2 += (f1[i] - f0[i]) / scale[i]; @@ -100,6 +102,7 @@ double select_initial_step( * @param C Coefficients for incrementing time for consecutive RK stages. * @param K Storage array for putting RK stages here. Stages are stored in rows. The last * row is a linear combination of the previous rows with coefficients. + * @param rhs_args Extra parameters to pass to the RHS function. * @return std::array, 2> Solution at t + h and derivative thereof. */ std::array, 2> rk_step( @@ -111,7 +114,8 @@ std::array, 2> rk_step( std::vector> const& A, std::vector const& B, std::vector const& C, - std::vector>& K + std::vector>& K, + void *const rhs_args ) { int const n_eqns = y.size(); @@ -130,7 +134,7 @@ std::array, 2> rk_step( for (auto i = 0; i < n_eqns; i++) { dy[i] += y[i]; } - fun(t + C[s] * h, dy.data(), K[s].data(), NULL); + fun(t + C[s] * h, dy.data(), K[s].data(), rhs_args); } std::vector y_new(n_eqns, double(0)), f_new(n_eqns); @@ -142,7 +146,7 @@ std::array, 2> rk_step( for (auto i = 0; i < n_eqns; i++) { y_new[i] += y[i]; } - fun(t + h, y_new.data(), f_new.data(), NULL); + fun(t + h, y_new.data(), f_new.data(), rhs_args); for (auto i = 0; i < n_eqns; i++) { K[n_stages][i] = f_new[i]; } @@ -280,12 +284,14 @@ void dgemv(const std::string trans, const size_t m, const size_t n, * @param atol Absolute tolerance. * @param first_step Initial step size; default is "None" which means that the algorithm * should choose. + * @param rhs_args Extra parameters to pass to the RHS function. */ RK45::RK45(RK45_ODE_SYSTEM_TYPE rhs_handle, double t0, std::vector &y0, double tf, double largest_step, double _rtol, double _atol, - double first_step) + double first_step, + void *const rhs_args) { _m_status = "running"; m_fun = rhs_handle; @@ -307,7 +313,7 @@ RK45::RK45(RK45_ODE_SYSTEM_TYPE rhs_handle, double t0, std::vector &y0, m_y = y0; m_n_eqns = m_y.size(); m_f.resize(m_n_eqns); - m_fun(m_t, m_y.data(), m_f.data(), NULL); + m_fun(m_t, m_y.data(), m_f.data(), rhs_args); if (largest_step <= double(0)) throw std::invalid_argument("'largest_step' must be positive."); @@ -317,7 +323,7 @@ RK45::RK45(RK45_ODE_SYSTEM_TYPE rhs_handle, double t0, std::vector &y0, throw std::invalid_argument("'first_step' exceeds bounds."); m_h_abs = first_step; } else { - m_h_abs = select_initial_step(m_fun, m_t, m_y, m_f, m_direction, _m_error_estimation_order, m_rtol, m_atol); + m_h_abs = select_initial_step(m_fun, m_t, m_y, m_f, m_direction, _m_error_estimation_order, m_rtol, m_atol, rhs_args); } _m_h_previous = -1.0; @@ -371,9 +377,10 @@ double RK45::_estimate_error_norm(double h, std::vector scale) /** * @brief Computations for one integration step. * + * @param rhs_args Extra parameters to pass to the RHS function. * @return bool Whether the step was successful. */ -bool RK45::step() +bool RK45::step(void *const rhs_args) { double min_step = double(10) * std::numeric_limits::epsilon(); double h_abs; @@ -403,7 +410,7 @@ bool RK45::step() h = t_new - m_t; h_abs = std::abs(h); - auto y_and_f_new = rk_step(m_fun, m_t, m_y, m_f, h, _m_A, _m_B, _m_C, _m_K); + auto y_and_f_new = rk_step(m_fun, m_t, m_y, m_f, h, _m_A, _m_B, _m_C, _m_K, rhs_args); y_new = y_and_f_new[0]; f_new = y_and_f_new[1]; diff --git a/src/RK45.h b/src/RK45.h index 88a10ab..6e2635a 100644 --- a/src/RK45.h +++ b/src/RK45.h @@ -56,7 +56,8 @@ class RK45 { double largest_step=std::numeric_limits::max(), double _rtol=1e-3, double _atol=1e-6, - double first_step=double(-1)); + double first_step=double(-1), + void *const rhs_args=NULL); ~RK45(); RK45_ODE_SYSTEM_TYPE m_fun; @@ -72,7 +73,7 @@ class RK45 { double m_t; int m_direction; - bool step(); + bool step(void *const rhs_args); private: static int const _m_order = 5; @@ -94,7 +95,6 @@ class RK45 { double _m_h_previous; std::vector> _m_K; - std::string _m_status; std::vector _estimate_error(double h); diff --git a/src/RK45_wrapper.cpp b/src/RK45_wrapper.cpp index fdd8641..b9b8b99 100644 --- a/src/RK45_wrapper.cpp +++ b/src/RK45_wrapper.cpp @@ -61,7 +61,7 @@ void rk45_wrapper( y0[i] = u0[i]; } - RK45 rk45(rhs, t0, y0, tf, mxstep, rtol, atol, dt0); + RK45 rk45(rhs, t0, y0, tf, mxstep, rtol, atol, dt0, data); for (auto i = 0; i < neq; i++) { usol[i + (0) * neq] = rk45.m_y[i]; @@ -70,7 +70,7 @@ void rk45_wrapper( int itnum = 1; while (itnum < itf + 1) { - bool successful_step = rk45.step(); + bool successful_step = rk45.step(data); if (not successful_step) { *success = 0; From 15063819398e88563aeb2f31cd0b62bfc1b5ca4a Mon Sep 17 00:00:00 2001 From: tbridel Date: Tue, 6 Dec 2022 08:17:22 +0100 Subject: [PATCH 7/8] Roll back testing rk45 on implicit robber problem --- benchmark/benchmark_robber.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmark/benchmark_robber.py b/benchmark/benchmark_robber.py index d6c2933..07f9b80 100644 --- a/benchmark/benchmark_robber.py +++ b/benchmark/benchmark_robber.py @@ -1,4 +1,4 @@ -from numbalsoda import lsoda_sig, lsoda, rk45 +from numbalsoda import lsoda_sig, lsoda import numba as nb import numpy as np from scipy.integrate import solve_ivp From 739d93022c55869534bf6bbe3d7487fa33eaf022 Mon Sep 17 00:00:00 2001 From: tbridel Date: Tue, 6 Dec 2022 08:21:03 +0100 Subject: [PATCH 8/8] Update documentation and benchmarks. --- README.md | 7 +++++-- benchmark/README.md | 1 + 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index a174822..4923c61 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ `numbalsoda` also wraps the `dop853` explicit Runge-Kutta method from this repository: https://github.com/jacobwilliams/dop853 -This package is very similar to `scipy.integrate.solve_ivp` ([see here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)), when you set `method = 'LSODA'` or `method = DOP853`. But, `scipy.integrate.solve_ivp` invokes the python interpreter every time step which can be slow. Also, `scipy.integrate.solve_ivp` can not be used within numba jit-compiled python functions. In contrast, `numbalsoda` never invokes the python interpreter during integration and can be used within a numba compiled function which makes `numbalsoda` a lot faster than scipy for most problems, and achieves similar performance to Julia's DifferentialEquations.jl in some cases (see `benchmark` folder). +This package is very similar to `scipy.integrate.solve_ivp` ([see here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)), when you set `method = 'LSODA'`, `method = DOP853` or `method = RK45`. But, `scipy.integrate.solve_ivp` invokes the python interpreter every time step which can be slow. Also, `scipy.integrate.solve_ivp` can not be used within numba jit-compiled python functions. In contrast, `numbalsoda` never invokes the python interpreter during integration and can be used within a numba compiled function which makes `numbalsoda` a lot faster than scipy for most problems, and achieves similar performance to Julia's DifferentialEquations.jl in some cases (see `benchmark` folder). ## Installation Conda: @@ -19,7 +19,7 @@ python -m pip install numbalsoda ## Basic usage ```python -from numbalsoda import lsoda_sig, lsoda, dop853 +from numbalsoda import lsoda_sig, lsoda, dop853, rk45 from numba import njit, cfunc import numpy as np @@ -39,6 +39,9 @@ usol, success = lsoda(funcptr, u0, t_eval, data = data) # integrate with dop853 method usol1, success1 = dop853(funcptr, u0, t_eval, data = data) +# integrate with rk45 method (automatic adaptive time stepping) +tsol2, usol2, tf2, success2 = rk45(funcptr, u0, -1.0, np.amin(t_eval), np.amax(t_eval), 100000, data = data) + # usol = solution # success = True/False ``` diff --git a/benchmark/README.md b/benchmark/README.md index 3b3bdf4..4c4d6d2 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -12,6 +12,7 @@ Also, Scipy might be faster than numbalsoda for stiff problems with > ~500 ODEs | -------------------------- | ---------- | ---------------------- | | numbalsoda lsoda | 2.594 ms | 1x | | numbalsoda dop853 | 0.700 ms | 0.27x | +| numbalsoda rk45 | 5.785 ms | 2.23x | Scipy LSODA | 109.349 ms | 42x | | Scipy RK45 | 292.826 ms | 113x | | Scipy DOP853 | 180.997 ms | 70x |