From ec3289b465e6ead6889127f2f8f3d1b1cc52d51a Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Tue, 30 Apr 2024 11:03:13 +0200 Subject: [PATCH 01/23] Small simplification. --- tools/gha_conda_docs.sh | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/tools/gha_conda_docs.sh b/tools/gha_conda_docs.sh index ff24a2df..8d410ecb 100644 --- a/tools/gha_conda_docs.sh +++ b/tools/gha_conda_docs.sh @@ -58,12 +58,7 @@ cd $HEYOKA_PY_PROJECT_DIR cd doc -# NOTE: run linkcheck only if we are on a pull request. -if [[ -n "${GITHUB_HEAD_REF}" ]]; then - make html linkcheck doctest -else - make html doctest -fi +make html doctest set +e set +x From 15aaf3356a7be01f6c95b7a241d81ab3efea14ce Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 23 May 2024 13:33:58 +0200 Subject: [PATCH 02/23] Anothe small doc addition. --- doc/notebooks/ensemble_mode.ipynb | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/doc/notebooks/ensemble_mode.ipynb b/doc/notebooks/ensemble_mode.ipynb index 0caffb2b..ae868f47 100644 --- a/doc/notebooks/ensemble_mode.ipynb +++ b/doc/notebooks/ensemble_mode.ipynb @@ -9,6 +9,10 @@ "\n", "# Ensemble propagations\n", "\n", + "```{versionadded} 0.17.0\n", + "\n", + "```\n", + "\n", "Starting with version 0.17.0, heyoka.py offers support for\n", "*ensemble propagations*. In ensemble mode, multiple distinct\n", "instances of the same ODE system are integrated in parallel,\n", From 4ff187454dfffe126142f3aac6fda61823a1aaf5 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Fri, 7 Jun 2024 11:14:15 +0200 Subject: [PATCH 03/23] Initial exposition of dfun. --- heyoka/expose_expression.cpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/heyoka/expose_expression.cpp b/heyoka/expose_expression.cpp index 8d5fb8ab..5624e244 100644 --- a/heyoka/expose_expression.cpp +++ b/heyoka/expose_expression.cpp @@ -490,6 +490,19 @@ void expose_expression(py::module_ &m) }, "y"_a.noconvert(), "x"_a.noconvert()); + // dfun(). + m.def( + "dfun", + [](std::string name, std::vector args, + std::optional>> didx) { + if (didx) { + return hey::dfun(std::move(name), std::move(args), std::move(*didx)); + } else { + return hey::dfun(std::move(name), std::move(args)); + } + }, + "name"_a, "args"_a, "didx"_a = py::none{}); + // Time. m.attr("time") = hey::time; From dfb48d0f007fe9ac11eb5cae37c1ff651f27f646 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Sun, 9 Jun 2024 21:56:22 +0200 Subject: [PATCH 04/23] Initial code for the exposition of variational ODEs. --- heyoka/CMakeLists.txt | 1 + heyoka/cfunc.cpp | 19 +- heyoka/common_utils.cpp | 64 +++++++ heyoka/common_utils.hpp | 16 ++ heyoka/core.cpp | 4 + heyoka/expose_batch_integrators.cpp | 267 ++++++++++++++++++++-------- heyoka/expose_expression.cpp | 88 ++------- heyoka/expose_var_ode_sys.cpp | 51 ++++++ heyoka/expose_var_ode_sys.hpp | 23 +++ heyoka/taylor_expose_integrator.cpp | 166 +++++++++++++++-- 10 files changed, 522 insertions(+), 177 deletions(-) create mode 100644 heyoka/expose_var_ode_sys.cpp create mode 100644 heyoka/expose_var_ode_sys.hpp diff --git a/heyoka/CMakeLists.txt b/heyoka/CMakeLists.txt index 24991218..34c6fc19 100644 --- a/heyoka/CMakeLists.txt +++ b/heyoka/CMakeLists.txt @@ -58,6 +58,7 @@ set(_HEYOKA_PY_CORE_SOURCES dtypes.cpp custom_casters.cpp expose_expression.cpp + expose_var_ode_sys.cpp expose_batch_integrators.cpp numpy_memory.cpp expose_models.cpp diff --git a/heyoka/cfunc.cpp b/heyoka/cfunc.cpp index 5b55ea67..8dd564bc 100644 --- a/heyoka/cfunc.cpp +++ b/heyoka/cfunc.cpp @@ -123,19 +123,6 @@ void expose_add_cfunc_impl(py::module &m, const char *suffix) // Fetch the dtype corresponding to T. const auto dt = get_dtype(); - // Small helper to facilitate the conversion of an iterable into - // a contiguous NumPy array of type dt. - const auto as_carray = [dt](const py::iterable &v) { - using namespace pybind11::literals; - - py::array ret = py::module_::import("numpy").attr("ascontiguousarray")(v, "dtype"_a = py::dtype(dt)); - - assert(ret.dtype().num() == dt); - assert(is_npy_array_carray(ret)); - - return ret; - }; - // Fetch the inputs array. // NOTE: this will either fetch the existing array, or convert // the input iterable to a py::array on the fly. @@ -162,7 +149,7 @@ void expose_add_cfunc_impl(py::module &m, const char *suffix) return std::move(v); } else { - return as_carray(v); + return as_carray(v, dt); } }, inputs_ob); @@ -261,7 +248,7 @@ void expose_add_cfunc_impl(py::module &m, const char *suffix) return std::move(v); } else { - return as_carray(v); + return as_carray(v, dt); } }, *pars_ob); @@ -328,7 +315,7 @@ void expose_add_cfunc_impl(py::module &m, const char *suffix) return std::move(v); } else { - return as_carray(v); + return as_carray(v, dt); } }, std::get(tm)); diff --git a/heyoka/common_utils.cpp b/heyoka/common_utils.cpp index 4524e98e..f89e9493 100644 --- a/heyoka/common_utils.cpp +++ b/heyoka/common_utils.cpp @@ -9,8 +9,13 @@ #include #include +#include #include #include +#include + +#include +#include #include #include @@ -33,6 +38,7 @@ #endif +#include #include #include "common_utils.hpp" @@ -108,4 +114,62 @@ bool with_pybind11_eh_impl() } // namespace detail +std::pair +dtens_t_it::operator()(const std::pair &p) const +{ + const auto &[sv_idx, ex] = p; + + return std::make_pair( + sparse_to_dense(sv_idx, boost::numeric_cast(dt->get_nargs())), ex); +} + +heyoka::dtens::v_idx_t dtens_t_it::sparse_to_dense(const heyoka::dtens::sv_idx_t &sv_idx, + heyoka::dtens::v_idx_t::size_type nargs) +{ + // Init the dense vector from the component index. + heyoka::dtens::v_idx_t ret{sv_idx.first}; + + // Transform the sparse index/order pairs into dense format. + // NOTE: no overflow check needed on ++idx because dtens ensures that + // the number of variables can be represented by std::uint32_t. + std::uint32_t idx = 0; + for (auto it = sv_idx.second.begin(); it != sv_idx.second.end(); ++idx) { + if (it->first == idx) { + // The current index shows up in the sparse vector, + // fetch the corresponding order and move to the next + // element of the sparse vector. + ret.push_back(it->second); + assert(it->second != 0u); + ++it; + } else { + // The current index does not show up in the sparse + // vector, set the order to zero. + ret.push_back(0); + } + } + + // Sanity check on the number of diff variables + // inferred from the sparse vector. + assert(ret.size() - 1u <= nargs); + + // Pad missing values at the end of ret. + ret.resize(boost::safe_numerics::safe(nargs) + 1); + + return ret; +} + +// Small helper to facilitate the conversion of an iterable into +// a contiguous NumPy array of type dt. +py::array as_carray(const py::iterable &v, int dt) +{ + using namespace pybind11::literals; + + py::array ret = py::module_::import("numpy").attr("ascontiguousarray")(v, "dtype"_a = py::dtype(dt)); + + assert(ret.dtype().num() == dt); + assert(is_npy_array_carray(ret)); + + return ret; +} + } // namespace heyoka_py diff --git a/heyoka/common_utils.hpp b/heyoka/common_utils.hpp index c5b4afa8..83f5e0bb 100644 --- a/heyoka/common_utils.hpp +++ b/heyoka/common_utils.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #if defined(__GLIBCXX__) @@ -24,6 +25,7 @@ #include +#include #include // NOTE: implementation of Py_SET_TYPE() for Python < 3.9. See: @@ -198,6 +200,20 @@ bool with_pybind11_eh(const F &f) } } +// Functor to transform on-the-fly the content of a dtens +// from sparse format into dense format. +struct dtens_t_it { + const heyoka::dtens *dt = nullptr; + + std::pair + operator()(const std::pair &) const; + + // Helper to implement the conversion from sparse to dense format. + static heyoka::dtens::v_idx_t sparse_to_dense(const heyoka::dtens::sv_idx_t &, heyoka::dtens::v_idx_t::size_type); +}; + +py::array as_carray(const py::iterable &, int); + } // namespace heyoka_py #endif diff --git a/heyoka/core.cpp b/heyoka/core.cpp index 2145e5e9..93858165 100644 --- a/heyoka/core.cpp +++ b/heyoka/core.cpp @@ -61,6 +61,7 @@ #include "expose_models.hpp" #include "expose_real.hpp" #include "expose_real128.hpp" +#include "expose_var_ode_sys.hpp" #include "logging.hpp" #include "pickle_wrappers.hpp" #include "setup_sympy.hpp" @@ -197,6 +198,9 @@ PYBIND11_MODULE(core, m) // Expression. heypy::expose_expression(m); + // var_ode_sys. + heypy::expose_var_ode_sys(m); + // Models. heypy::expose_models(m); diff --git a/heyoka/expose_batch_integrators.cpp b/heyoka/expose_batch_integrators.cpp index 42e87c15..8ff14cde 100644 --- a/heyoka/expose_batch_integrators.cpp +++ b/heyoka/expose_batch_integrators.cpp @@ -6,9 +6,12 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include +#include #include #include #include +#include #include #include #include @@ -31,6 +34,7 @@ #include #include #include +#include #include "common_utils.hpp" #include "custom_casters.hpp" @@ -43,6 +47,7 @@ namespace heyoka_py { namespace py = pybind11; +namespace hey = heyoka; namespace detail { @@ -50,10 +55,29 @@ namespace detail namespace { +// Helper to fetch the tstate from a variational integrator. +// Extracted here for re-use. +template +auto fetch_tstate(const py::object &o) +{ + const auto *ta = py::cast *>(o); + + assert(ta->get_tstate().size() % ta->get_batch_size() == 0u); + const auto n_orig_sv = boost::numeric_cast(ta->get_n_orig_sv()); + const auto bs = boost::numeric_cast(ta->get_batch_size()); + + auto ret + = py::array(py::dtype(get_dtype()), py::array::ShapeContainer{n_orig_sv, bs}, ta->get_tstate().data(), o); + + // Ensure the returned array is read-only. + ret.attr("flags").attr("writeable") = false; + + return ret; +} + template void expose_batch_integrator_impl(py::module_ &m, const std::string &suffix) { - namespace hey = heyoka; namespace kw = hey::kw; using namespace pybind11::literals; @@ -64,10 +88,11 @@ void expose_batch_integrator_impl(py::module_ &m, const std::string &suffix) using sys_t = std::vector>; // Implementation of the ctor. - auto tab_ctor_impl = [](const sys_t &sys, const py::iterable &state_ob, std::optional time_ob, - std::optional pars_ob, T tol, bool high_accuracy, bool compact_mode, - std::vector tes, std::vector ntes, bool parallel_mode, unsigned opt_level, - bool force_avx512, bool slp_vectorize, bool fast_math) { + auto tab_ctor_impl = [](std::variant vsys, const py::iterable &state_ob, + std::optional time_ob, std::optional pars_ob, T tol, + bool high_accuracy, bool compact_mode, std::vector tes, std::vector ntes, + bool parallel_mode, unsigned opt_level, bool force_avx512, bool slp_vectorize, + bool fast_math) { // Fetch the dtype corresponding to T. const auto dt = get_dtype(); @@ -115,77 +140,83 @@ void expose_batch_integrator_impl(py::module_ &m, const std::string &suffix) pars = py::cast>(pars_arr.attr("flatten")()); } - if (time_ob) { - // Times provided. - py::array time_arr = *time_ob; - if (time_arr.ndim() != 1 || boost::numeric_cast(time_arr.shape(0)) != batch_size) { - py_throw(PyExc_ValueError, - fmt::format("Invalid time vector passed to the constructor of a batch integrator: " - "the expected array shape is ({}), but the input array has either the wrong " - "number of dimensions or the wrong shape", - batch_size) - .c_str()); - } - auto time = py::cast>(time_arr); - - // NOTE: GIL release is fine here even if the events contain - // Python objects, as the event vectors are moved in - // upon construction and thus we should never end up calling - // into the interpreter. - py::gil_scoped_release release; - - return hey::taylor_adaptive_batch{sys, - std::move(state), - batch_size, - kw::time = std::move(time), - kw::tol = tol, - kw::high_accuracy = high_accuracy, - kw::compact_mode = compact_mode, - kw::pars = std::move(pars), - kw::t_events = std::move(tes), - kw::nt_events = std::move(ntes), - kw::parallel_mode = parallel_mode, - kw::opt_level = opt_level, - kw::force_avx512 = force_avx512, - kw::slp_vectorize = slp_vectorize, - kw::fast_math = fast_math}; - } else { - // Times not provided. - - // NOTE: GIL release is fine here even if the events contain - // Python objects, as the event vectors are moved in - // upon construction and thus we should never end up calling - // into the interpreter. - py::gil_scoped_release release; - - return hey::taylor_adaptive_batch{sys, - std::move(state), - batch_size, - kw::tol = tol, - kw::high_accuracy = high_accuracy, - kw::compact_mode = compact_mode, - kw::pars = std::move(pars), - kw::t_events = std::move(tes), - kw::nt_events = std::move(ntes), - kw::parallel_mode = parallel_mode, - kw::opt_level = opt_level, - kw::force_avx512 = force_avx512, - kw::slp_vectorize = slp_vectorize, - kw::fast_math = fast_math}; - } + return std::visit( + [&](auto &sys) { + if (time_ob) { + // Times provided. + py::array time_arr = *time_ob; + if (time_arr.ndim() != 1 || boost::numeric_cast(time_arr.shape(0)) != batch_size) { + py_throw( + PyExc_ValueError, + fmt::format("Invalid time vector passed to the constructor of a batch integrator: " + "the expected array shape is ({}), but the input array has either the wrong " + "number of dimensions or the wrong shape", + batch_size) + .c_str()); + } + auto time = py::cast>(time_arr); + + // NOTE: GIL release is fine here even if the events contain + // Python objects, as the event vectors are moved in + // upon construction and thus we should never end up calling + // into the interpreter. + py::gil_scoped_release release; + + return hey::taylor_adaptive_batch{std::move(sys), + std::move(state), + batch_size, + kw::time = std::move(time), + kw::tol = tol, + kw::high_accuracy = high_accuracy, + kw::compact_mode = compact_mode, + kw::pars = std::move(pars), + kw::t_events = std::move(tes), + kw::nt_events = std::move(ntes), + kw::parallel_mode = parallel_mode, + kw::opt_level = opt_level, + kw::force_avx512 = force_avx512, + kw::slp_vectorize = slp_vectorize, + kw::fast_math = fast_math}; + } else { + // Times not provided. + + // NOTE: GIL release is fine here even if the events contain + // Python objects, as the event vectors are moved in + // upon construction and thus we should never end up calling + // into the interpreter. + py::gil_scoped_release release; + + return hey::taylor_adaptive_batch{std::move(sys), + std::move(state), + batch_size, + kw::tol = tol, + kw::high_accuracy = high_accuracy, + kw::compact_mode = compact_mode, + kw::pars = std::move(pars), + kw::t_events = std::move(tes), + kw::nt_events = std::move(ntes), + kw::parallel_mode = parallel_mode, + kw::opt_level = opt_level, + kw::force_avx512 = force_avx512, + kw::slp_vectorize = slp_vectorize, + kw::fast_math = fast_math}; + } + }, + vsys); }; py::class_> tab_c(m, fmt::format("taylor_adaptive_batch_{}", suffix).c_str(), py::dynamic_attr{}); tab_c - .def(py::init([tab_ctor_impl](const sys_t &sys, const py::iterable &state, std::optional time, - std::optional pars, T tol, bool high_accuracy, bool compact_mode, - std::vector tes, std::vector ntes, bool parallel_mode, - unsigned opt_level, bool force_avx512, bool slp_vectorize, bool fast_math) { - return tab_ctor_impl(sys, state, std::move(time), std::move(pars), tol, high_accuracy, compact_mode, - std::move(tes), std::move(ntes), parallel_mode, opt_level, force_avx512, - slp_vectorize, fast_math); + .def(py::init([tab_ctor_impl](std::variant vsys, const py::iterable &state, + std::optional time, std::optional pars, T tol, + bool high_accuracy, bool compact_mode, std::vector tes, + std::vector ntes, bool parallel_mode, unsigned opt_level, + bool force_avx512, bool slp_vectorize, bool fast_math) { + return tab_ctor_impl(std::move(vsys), state, std::move(time), std::move(pars), tol, high_accuracy, + compact_mode, std::move(tes), std::move(ntes), parallel_mode, opt_level, + force_avx512, slp_vectorize, fast_math); }), "sys"_a, "state"_a, "time"_a = py::none{}, "pars"_a = py::none{}, "tol"_a.noconvert() = static_cast(0), "high_accuracy"_a = false, "compact_mode"_a = false, "t_events"_a = py::list{}, "nt_events"_a = py::list{}, @@ -509,6 +540,102 @@ void expose_batch_integrator_impl(py::module_ &m, const std::string &suffix) .def_property_readonly("compact_mode", &hey::taylor_adaptive_batch::get_compact_mode) .def_property_readonly("high_accuracy", &hey::taylor_adaptive_batch::get_high_accuracy) .def_property_readonly("with_events", &hey::taylor_adaptive_batch::with_events) + // Variational-specific bits. + .def_property_readonly("n_orig_sv", &hey::taylor_adaptive_batch::get_n_orig_sv) + .def_property_readonly("is_variational", &hey::taylor_adaptive_batch::is_variational) + .def_property_readonly("vargs", &hey::taylor_adaptive_batch::get_vargs) + .def_property_readonly("tstate", &fetch_tstate) + .def( + "get_vslice", + [](const hey::taylor_adaptive_batch &ta, std::uint32_t order, std::optional component) { + const auto ret = component ? ta.get_vslice(*component, order) : ta.get_vslice(order); + + return py::slice(boost::numeric_cast(ret.first), + boost::numeric_cast(ret.second), {}); + }, + "order"_a, "component"_a = py::none{}) + .def( + "get_mindex", + [](const hey::taylor_adaptive_batch &ta, std::uint32_t i) { + const auto &ret = ta.get_mindex(i); + + return dtens_t_it::sparse_to_dense( + ret, boost::numeric_cast(ta.get_vargs().size())); + }, + "i"_a) + .def("eval_taylor_map", + [](py::object &o, std::variant in) { + auto *ta = py::cast *>(o); + + // Fetch the dtype corresponding to T. + const auto dt = get_dtype(); + + // Fetch the inputs array. + // NOTE: this will either fetch the existing array, or convert + // the input iterable to a py::array on the fly. + const auto inputs = std::visit( + [&](U &v) { + if constexpr (std::same_as) { + // Check the dtype. + if (v.dtype().num() != dt) [[unlikely]] { + py_throw( + PyExc_TypeError, + fmt::format( + "Invalid dtype detected for the inputs of a Taylor map evaluation: the " + "expected dtype is '{}', but the dtype of the inputs array is '{}' instead", + str(py::dtype(dt)), str(v.dtype())) + .c_str()); + } + + // Check that the inputs array is a C-style contiguous array. + if (!is_npy_array_carray(v)) [[unlikely]] { + py_throw(PyExc_ValueError, + "Invalid inputs array detected in a Taylor map evaluation: the array is not " + "C-style contiguous, please " + "consider using numpy.ascontiguousarray() to turn it into one"); + } + + return std::move(v); + } else { + return as_carray(v, dt); + } + }, + in); + + // Validate the number of dimensions for the inputs. + // NOTE: this needs to be done regardless of the original type of in. + if (inputs.ndim() != 2) [[unlikely]] { + py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " + "of a Taylor map has {} dimensions, " + "but it must have 2 dimensions instead", + inputs.ndim()) + .c_str()); + } + + // Validate the shape for the inputs. + if (boost::numeric_cast(inputs.shape(0)) != ta->get_n_orig_sv()) [[unlikely]] { + py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " + "of a Taylor map has {} rows, " + "but it must have {} rows instead", + inputs.shape(0), ta->get_n_orig_sv()) + .c_str()); + } + + if (boost::numeric_cast(inputs.shape(1)) != ta->get_batch_size()) [[unlikely]] { + py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " + "of a Taylor map has {} columns, " + "but it must have {} columns instead", + inputs.shape(1), ta->get_batch_size()) + .c_str()); + } + + // Run the evaluation. + const auto *data_ptr = static_cast(inputs.data()); + ta->eval_taylor_map(std::ranges::subrange(data_ptr, data_ptr + inputs.size())); + + // Return the tstate. + return fetch_tstate(o); + }) // Event detection. .def_property_readonly("with_events", &hey::taylor_adaptive_batch::with_events) .def_property_readonly("te_cooldowns", &hey::taylor_adaptive_batch::get_te_cooldowns) diff --git a/heyoka/expose_expression.cpp b/heyoka/expose_expression.cpp index 8d5fb8ab..aa1075fe 100644 --- a/heyoka/expose_expression.cpp +++ b/heyoka/expose_expression.cpp @@ -24,7 +24,6 @@ #include #include -#include #include #include @@ -52,67 +51,12 @@ #include "common_utils.hpp" #include "custom_casters.hpp" #include "docstrings.hpp" +#include "expose_expression.hpp" #include "pickle_wrappers.hpp" namespace heyoka_py { -namespace detail -{ - -namespace -{ - -template -using uncvref_t = std::remove_cv_t>; - -// Functor to transform on-the-fly the content of a dtens -// from sparse format into dense format. -struct dtens_t_it { - const heyoka::dtens *dt = nullptr; - - std::pair - operator()(const std::pair &p) const - { - const auto &[sv_idx, ex] = p; - - // Init the dense vector from the component index. - heyoka::dtens::v_idx_t ret{sv_idx.first}; - - // Transform the sparse index/order pairs into dense format. - // NOTE: no overflow check needed on ++idx because dtens ensures that - // the number of variables can be represented by std::uint32_t. - std::uint32_t idx = 0; - for (auto it = sv_idx.second.begin(); it != sv_idx.second.end(); ++idx) { - if (it->first == idx) { - // The current index shows up in the sparse vector, - // fetch the corresponding order and move to the next - // element of the sparse vector. - ret.push_back(it->second); - assert(it->second != 0u); - ++it; - } else { - // The current index does not show up in the sparse - // vector, set the order to zero. - ret.push_back(0); - } - } - - // Sanity check on the number of diff variables - // inferred from the sparse vector. - assert(ret.size() - 1u <= dt->get_nargs()); - - // Pad missing values at the end of ret. - ret.resize(boost::safe_numerics::safe(dt->get_nargs()) + 1); - - return std::make_pair(std::move(ret), ex); - } -}; - -} // namespace - -} // namespace detail - namespace py = pybind11; // NOTE: regarding single-precision support: we only expose the expression ctor, @@ -381,8 +325,8 @@ void expose_expression(py::module_ &m) [](const mvf_arg &e, const mvf_arg &M) { return std::visit( [](const auto &a, const auto &b) -> hey::expression { - using tp1 = detail::uncvref_t; - using tp2 = detail::uncvref_t; + using tp1 = std::remove_cvref_t; + using tp2 = std::remove_cvref_t; if constexpr (!std::is_same_v && !std::is_same_v) { py_throw(PyExc_TypeError, "At least one of the arguments of kepE() must be an expression"); @@ -400,9 +344,9 @@ void expose_expression(py::module_ &m) [](const mvf_arg &h, const mvf_arg &k, const mvf_arg &lam) { return std::visit( [](const auto &a, const auto &b, const auto &c) -> hey::expression { - using tp1 = detail::uncvref_t; - using tp2 = detail::uncvref_t; - using tp3 = detail::uncvref_t; + using tp1 = std::remove_cvref_t; + using tp2 = std::remove_cvref_t; + using tp3 = std::remove_cvref_t; constexpr auto tp1_num = static_cast(!std::is_same_v); constexpr auto tp2_num = static_cast(!std::is_same_v); @@ -438,9 +382,9 @@ void expose_expression(py::module_ &m) [](const mvf_arg &s0, const mvf_arg &c0, const mvf_arg &DM) { return std::visit( [](const auto &a, const auto &b, const auto &c) -> hey::expression { - using tp1 = detail::uncvref_t; - using tp2 = detail::uncvref_t; - using tp3 = detail::uncvref_t; + using tp1 = std::remove_cvref_t; + using tp2 = std::remove_cvref_t; + using tp3 = std::remove_cvref_t; constexpr auto tp1_num = static_cast(!std::is_same_v); constexpr auto tp2_num = static_cast(!std::is_same_v); @@ -477,8 +421,8 @@ void expose_expression(py::module_ &m) [](const mvf_arg &y, const mvf_arg &x) { return std::visit( [](const auto &a, const auto &b) -> hey::expression { - using tp1 = detail::uncvref_t; - using tp2 = detail::uncvref_t; + using tp1 = std::remove_cvref_t; + using tp2 = std::remove_cvref_t; if constexpr (!std::is_same_v && !std::is_same_v) { py_throw(PyExc_TypeError, "At least one of the arguments of atan2() must be an expression"); @@ -553,7 +497,7 @@ void expose_expression(py::module_ &m) const auto s_idx = boost::numeric_cast::difference_type>(idx); - return detail::dtens_t_it{&dt}(dt.begin()[s_idx]); + return dtens_t_it{&dt}(dt.begin()[s_idx]); }); dtens_cl.def("__contains__", [](const hey::dtens &dt, const std::variant &v_idx_) { @@ -563,8 +507,8 @@ void expose_expression(py::module_ &m) dtens_cl.def( "__iter__", [](const hey::dtens &dt) { - auto t_begin = boost::iterators::make_transform_iterator(dt.begin(), detail::dtens_t_it{&dt}); - auto t_end = boost::iterators::make_transform_iterator(dt.end(), detail::dtens_t_it{&dt}); + auto t_begin = boost::iterators::make_transform_iterator(dt.begin(), dtens_t_it{&dt}); + auto t_end = boost::iterators::make_transform_iterator(dt.end(), dtens_t_it{&dt}); return py::make_key_iterator(t_begin, t_end); }, @@ -586,8 +530,8 @@ void expose_expression(py::module_ &m) [](const hey::dtens &dt, std::uint32_t order, std::optional component) { const auto sr = component ? dt.get_derivatives(*component, order) : dt.get_derivatives(order); - auto t_begin = boost::iterators::make_transform_iterator(sr.begin(), detail::dtens_t_it{&dt}); - auto t_end = boost::iterators::make_transform_iterator(sr.end(), detail::dtens_t_it{&dt}); + auto t_begin = boost::iterators::make_transform_iterator(sr.begin(), dtens_t_it{&dt}); + auto t_end = boost::iterators::make_transform_iterator(sr.end(), dtens_t_it{&dt}); return std::vector(t_begin, t_end); }, diff --git a/heyoka/expose_var_ode_sys.cpp b/heyoka/expose_var_ode_sys.cpp new file mode 100644 index 00000000..4b872500 --- /dev/null +++ b/heyoka/expose_var_ode_sys.cpp @@ -0,0 +1,51 @@ +// Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +// +// This file is part of the heyoka.py library. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include +#include +#include + +#include +#include + +#include + +#include "expose_var_ode_sys.hpp" + +namespace heyoka_py +{ + +namespace py = pybind11; + +void expose_var_ode_sys(py::module_ &m) +{ + namespace hey = heyoka; + // NOLINTNEXTLINE(google-build-using-namespace) + using namespace pybind11::literals; + + // var_args enum. + py::enum_(m, "var_args", py::arithmetic()) + .value("vars", hey::var_args::vars) + .value("params", hey::var_args::params) + .value("time", hey::var_args::time) + .value("all", hey::var_args::all) + .def("__or__", [](hey::var_args a, hey::var_args b) { return a | b; }); + + // var_ode_sys class. + py::class_ v_cl(m, "var_ode_sys", py::dynamic_attr{}); + v_cl.def(py::init([](const std::vector> &sys, + const std::variant> &args, + std::uint32_t order) { return hey::var_ode_sys(sys, args, order); }), + "sys"_a, "args"_a, "order"_a = static_cast(1)); + v_cl.def_property_readonly("sys", &hey::var_ode_sys::get_sys); + v_cl.def_property_readonly("vargs", &hey::var_ode_sys::get_vargs); + v_cl.def_property_readonly("n_orig_sv", &hey::var_ode_sys::get_n_orig_sv); + v_cl.def_property_readonly("order", &hey::var_ode_sys::get_order); +} + +} // namespace heyoka_py diff --git a/heyoka/expose_var_ode_sys.hpp b/heyoka/expose_var_ode_sys.hpp new file mode 100644 index 00000000..1b452567 --- /dev/null +++ b/heyoka/expose_var_ode_sys.hpp @@ -0,0 +1,23 @@ +// Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +// +// This file is part of the heyoka.py library. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef HEYOKA_PY_EXPOSE_VAR_ODE_SYS_HPP +#define HEYOKA_PY_EXPOSE_VAR_ODE_SYS_HPP + +#include + +namespace heyoka_py +{ + +namespace py = pybind11; + +void expose_var_ode_sys(py::module_ &); + +} // namespace heyoka_py + +#endif diff --git a/heyoka/taylor_expose_integrator.cpp b/heyoka/taylor_expose_integrator.cpp index d90cf4c7..aa7fb269 100644 --- a/heyoka/taylor_expose_integrator.cpp +++ b/heyoka/taylor_expose_integrator.cpp @@ -10,9 +10,11 @@ #include #include +#include #include -#include +#include #include +#include #include #include #include @@ -49,6 +51,7 @@ #include #include #include +#include #include "common_utils.hpp" #include "custom_casters.hpp" @@ -57,6 +60,12 @@ #include "step_cb_utils.hpp" #include "taylor_expose_integrator.hpp" +#if defined(HEYOKA_HAVE_REAL) + +#include "expose_real.hpp" + +#endif + namespace heyoka_py { @@ -78,6 +87,22 @@ constexpr bool default_cm = #endif ; +// Helper to fetch the tstate from a variational integrator. +// Extracted here for re-use. +template +auto fetch_tstate(const py::object &o) +{ + const auto *ta = py::cast *>(o); + auto ret = py::array(py::dtype(get_dtype()), + py::array::ShapeContainer{boost::numeric_cast(ta->get_tstate().size())}, + ta->get_tstate().data(), o); + + // Ensure the returned array is read-only. + ret.attr("flags").attr("writeable") = false; + + return ret; +} + // Implementation of the exposition of the scalar integrators. template void expose_taylor_integrator_impl(py::module &m, const std::string &suffix) @@ -91,30 +116,35 @@ void expose_taylor_integrator_impl(py::module &m, const std::string &suffix) using sys_t = std::vector>; py::class_> cl(m, (fmt::format("taylor_adaptive_{}", suffix)).c_str(), py::dynamic_attr{}); - cl.def(py::init([](const sys_t &sys, std::vector state, T time, std::vector pars, T tol, bool high_accuracy, - bool compact_mode, std::vector tes, std::vector ntes, bool parallel_mode, - unsigned opt_level, bool force_avx512, bool slp_vectorize, bool fast_math, long long prec) { + cl.def(py::init([](std::variant vsys, std::vector state, T time, std::vector pars, + T tol, bool high_accuracy, bool compact_mode, std::vector tes, std::vector ntes, + bool parallel_mode, unsigned opt_level, bool force_avx512, bool slp_vectorize, bool fast_math, + long long prec) { // NOTE: GIL release is fine here even if the events contain // Python objects, as the event vectors are moved in // upon construction and thus we should never end up calling // into the interpreter. py::gil_scoped_release release; - return hey::taylor_adaptive{sys, - std::move(state), - kw::time = time, - kw::tol = tol, - kw::high_accuracy = high_accuracy, - kw::compact_mode = compact_mode, - kw::pars = std::move(pars), - kw::t_events = std::move(tes), - kw::nt_events = std::move(ntes), - kw::parallel_mode = parallel_mode, - kw::opt_level = opt_level, - kw::force_avx512 = force_avx512, - kw::slp_vectorize = slp_vectorize, - kw::fast_math = fast_math, - kw::prec = prec}; + return std::visit( + [&](auto &sys) { + return hey::taylor_adaptive{std::move(sys), + std::move(state), + kw::time = time, + kw::tol = tol, + kw::high_accuracy = high_accuracy, + kw::compact_mode = compact_mode, + kw::pars = std::move(pars), + kw::t_events = std::move(tes), + kw::nt_events = std::move(ntes), + kw::parallel_mode = parallel_mode, + kw::opt_level = opt_level, + kw::force_avx512 = force_avx512, + kw::slp_vectorize = slp_vectorize, + kw::fast_math = fast_math, + kw::prec = prec}; + }, + vsys); }), "sys"_a, "state"_a.noconvert(), "time"_a.noconvert() = static_cast(0), "pars"_a.noconvert() = py::list{}, "tol"_a.noconvert() = static_cast(0), "high_accuracy"_a = false, "compact_mode"_a = default_cm, @@ -344,6 +374,104 @@ void expose_taylor_integrator_impl(py::module &m, const std::string &suffix) .def_property_readonly("compact_mode", &hey::taylor_adaptive::get_compact_mode) .def_property_readonly("high_accuracy", &hey::taylor_adaptive::get_high_accuracy) .def_property_readonly("with_events", &hey::taylor_adaptive::with_events) + // Variational-specific bits. + .def_property_readonly("n_orig_sv", &hey::taylor_adaptive::get_n_orig_sv) + .def_property_readonly("is_variational", &hey::taylor_adaptive::is_variational) + .def_property_readonly("vargs", &hey::taylor_adaptive::get_vargs) + .def_property_readonly("tstate", &fetch_tstate) + .def( + "get_vslice", + [](const hey::taylor_adaptive &ta, std::uint32_t order, std::optional component) { + const auto ret = component ? ta.get_vslice(*component, order) : ta.get_vslice(order); + + return py::slice(boost::numeric_cast(ret.first), + boost::numeric_cast(ret.second), {}); + }, + "order"_a, "component"_a = py::none{}) + .def( + "get_mindex", + [](const hey::taylor_adaptive &ta, std::uint32_t i) { + const auto &ret = ta.get_mindex(i); + + return dtens_t_it::sparse_to_dense( + ret, boost::numeric_cast(ta.get_vargs().size())); + }, + "i"_a) + .def("eval_taylor_map", + [](py::object &o, std::variant in) { + auto *ta = py::cast *>(o); + + // Fetch the dtype corresponding to T. + const auto dt = get_dtype(); + + // Fetch the inputs array. + // NOTE: this will either fetch the existing array, or convert + // the input iterable to a py::array on the fly. + const auto inputs = std::visit( + [&](U &v) { + if constexpr (std::same_as) { + // Check the dtype. + if (v.dtype().num() != dt) [[unlikely]] { + py_throw( + PyExc_TypeError, + fmt::format( + "Invalid dtype detected for the inputs of a Taylor map evaluation: the " + "expected dtype is '{}', but the dtype of the inputs array is '{}' instead", + str(py::dtype(dt)), str(v.dtype())) + .c_str()); + } + + // Check that the inputs array is a C-style contiguous array. + if (!is_npy_array_carray(v)) [[unlikely]] { + py_throw(PyExc_ValueError, + "Invalid inputs array detected in a Taylor map evaluation: the array is not " + "C-style contiguous, please " + "consider using numpy.ascontiguousarray() to turn it into one"); + } + + return std::move(v); + } else { + return as_carray(v, dt); + } + }, + in); + + // Validate the number of dimensions for the inputs. + // NOTE: this needs to be done regardless of the original type of in. + if (inputs.ndim() != 1) [[unlikely]] { + py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " + "of a Taylor map has {} dimensions, " + "but it must have 1 dimension instead", + inputs.ndim()) + .c_str()); + } + + // Validate the shape for the inputs. + if (boost::numeric_cast(inputs.shape(0)) != ta->get_n_orig_sv()) [[unlikely]] { + py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " + "of a Taylor map has {} elements, " + "but it must have {} elements instead", + inputs.shape(0), ta->get_n_orig_sv()) + .c_str()); + } + +#if defined(HEYOKA_HAVE_REAL) + + // Run the checks specific for mppp::real. + if constexpr (std::same_as) { + // Check that the inputs array contains values with the correct precision. + pyreal_check_array(inputs, ta->get_prec()); + } + +#endif + + // Run the evaluation. + const auto *data_ptr = static_cast(inputs.data()); + ta->eval_taylor_map(std::ranges::subrange(data_ptr, data_ptr + inputs.shape(0))); + + // Return the tstate. + return fetch_tstate(o); + }) // Event detection. .def_property_readonly("with_events", &hey::taylor_adaptive::with_events) .def_property_readonly("te_cooldowns", &hey::taylor_adaptive::get_te_cooldowns) From f61ba438a039cbfc370581d0fe8ec03bd6a24e0d Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Sun, 9 Jun 2024 22:15:31 +0200 Subject: [PATCH 05/23] More var_ode_sys exposition, testing. --- heyoka/CMakeLists.txt | 1 + heyoka/_test_model.py | 47 +++++++++++++++++++++++------------ heyoka/expose_var_ode_sys.cpp | 10 +++++++- heyoka/test.py | 2 ++ 4 files changed, 43 insertions(+), 17 deletions(-) diff --git a/heyoka/CMakeLists.txt b/heyoka/CMakeLists.txt index 34c6fc19..7f018229 100644 --- a/heyoka/CMakeLists.txt +++ b/heyoka/CMakeLists.txt @@ -30,6 +30,7 @@ set(HEYOKA_PY_PYTHON_FILES _test_vsop2013.py _test_elp2000.py _test_lagham.py + _test_var_ode_sys.py model/__init__.py callback/__init__.py ) diff --git a/heyoka/_test_model.py b/heyoka/_test_model.py index 8e6383be..cd4cb3bb 100644 --- a/heyoka/_test_model.py +++ b/heyoka/_test_model.py @@ -327,7 +327,11 @@ def test_nrlmsise00(self): "h", "lat", "lon", "f107", "f107a", "ap" ) nrlmsise00 = model.nrlmsise00_tn( - geodetic=[h, lat, lon], f107=f107, f107a=f107a, ap=ap, time_expr=time / 86400 + geodetic=[h, lat, lon], + f107=f107, + f107a=f107a, + ap=ap, + time_expr=time / 86400, ) nrlmsise00_cf = cfunc([nrlmsise00], vars=[h, lat, lon, f107, f107a, ap]) @@ -353,21 +357,32 @@ def test_nrlmsise00(self): def test_jb08(self): from . import model, make_vars, cfunc, time - h, lat, lon, f107a, f107, s107a, s107, m107a, m107, y107a, y107, dDstdT = ( - make_vars( - "h", - "lat", - "lon", - "f107a", - "f107", - "s107a", - "s107", - "m107a", - "m107", - "y107a", - "y107", - "dDstdT", - ) + ( + h, + lat, + lon, + f107a, + f107, + s107a, + s107, + m107a, + m107, + y107a, + y107, + dDstdT, + ) = make_vars( + "h", + "lat", + "lon", + "f107a", + "f107", + "s107a", + "s107", + "m107a", + "m107", + "y107a", + "y107", + "dDstdT", ) jb08 = model.jb08_tn( geodetic=[h, lat, lon], diff --git a/heyoka/expose_var_ode_sys.cpp b/heyoka/expose_var_ode_sys.cpp index 4b872500..ed3d231d 100644 --- a/heyoka/expose_var_ode_sys.cpp +++ b/heyoka/expose_var_ode_sys.cpp @@ -15,7 +15,9 @@ #include +#include "common_utils.hpp" #include "expose_var_ode_sys.hpp" +#include "pickle_wrappers.hpp" namespace heyoka_py { @@ -34,7 +36,8 @@ void expose_var_ode_sys(py::module_ &m) .value("params", hey::var_args::params) .value("time", hey::var_args::time) .value("all", hey::var_args::all) - .def("__or__", [](hey::var_args a, hey::var_args b) { return a | b; }); + .def("__or__", [](hey::var_args a, hey::var_args b) { return a | b; }) + .def("__and__", [](hey::var_args a, hey::var_args b) { return a & b; }); // var_ode_sys class. py::class_ v_cl(m, "var_ode_sys", py::dynamic_attr{}); @@ -46,6 +49,11 @@ void expose_var_ode_sys(py::module_ &m) v_cl.def_property_readonly("vargs", &hey::var_ode_sys::get_vargs); v_cl.def_property_readonly("n_orig_sv", &hey::var_ode_sys::get_n_orig_sv); v_cl.def_property_readonly("order", &hey::var_ode_sys::get_order); + // Copy/deepcopy. + v_cl.def("__copy__", copy_wrapper); + v_cl.def("__deepcopy__", deepcopy_wrapper, "memo"_a); + // Pickle support. + v_cl.def(py::pickle(&pickle_getstate_wrapper, &pickle_setstate_wrapper)); } } // namespace heyoka_py diff --git a/heyoka/test.py b/heyoka/test.py index 7f186fff..996b5cae 100644 --- a/heyoka/test.py +++ b/heyoka/test.py @@ -2079,6 +2079,7 @@ def run_test_suite(): _test_vsop2013, _test_elp2000, _test_lagham, + _test_var_ode_sys, ) import numpy as np from .model import nbody @@ -2098,6 +2099,7 @@ def run_test_suite(): suite.addTest( tl.loadTestsFromTestCase(_test_batch_integrator.batch_integrator_test_case) ) + suite.addTest(tl.loadTestsFromTestCase(_test_var_ode_sys.var_ode_sys_test_case)) suite.addTest(tl.loadTestsFromTestCase(_test_lagham.lagham_test_case)) suite.addTest(tl.loadTestsFromTestCase(_test_vsop2013.vsop2013_test_case)) suite.addTest(tl.loadTestsFromTestCase(_test_elp2000.elp2000_test_case)) From 47c61f2038f336a1397fbb27cbfffba916da5e43 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Sun, 9 Jun 2024 22:16:32 +0200 Subject: [PATCH 06/23] Missing file. --- heyoka/_test_var_ode_sys.py | 69 +++++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 heyoka/_test_var_ode_sys.py diff --git a/heyoka/_test_var_ode_sys.py b/heyoka/_test_var_ode_sys.py new file mode 100644 index 00000000..c44cee2e --- /dev/null +++ b/heyoka/_test_var_ode_sys.py @@ -0,0 +1,69 @@ +# Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +# +# This file is part of the heyoka.py library. +# +# This Source Code Form is subject to the terms of the Mozilla +# Public License v. 2.0. If a copy of the MPL was not distributed +# with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + + +import unittest as _ut + + +class var_ode_sys_test_case(_ut.TestCase): + def test_enum(self): + from . import var_args + + self.assertEqual(var_args.vars | var_args.time | var_args.params, var_args.all) + self.assertTrue((var_args.vars | var_args.time) & var_args.time) + + def test_basic(self): + from . import make_vars, var_ode_sys, var_args, sin, par, time + from copy import copy, deepcopy + from pickle import dumps, loads + + x, v = make_vars("x", "v") + + orig_sys = [(x, v), (v, -par[0] * sin(x) + time)] + + vsys = var_ode_sys(orig_sys, var_args.vars) + + self.assertEqual(orig_sys, vsys.sys[:2]) + self.assertEqual(vsys.vargs, [x, v]) + self.assertEqual(vsys.n_orig_sv, 2) + self.assertEqual(vsys.order, 1) + + vsys = var_ode_sys(sys=orig_sys, args=[v, time, x], order=2) + + self.assertEqual(orig_sys, vsys.sys[:2]) + self.assertEqual(vsys.vargs, [v, time, x]) + self.assertEqual(vsys.n_orig_sv, 2) + self.assertEqual(vsys.order, 2) + + vsys = var_ode_sys(orig_sys, var_args.vars | var_args.params) + + self.assertEqual(orig_sys, vsys.sys[:2]) + self.assertEqual(vsys.vargs, [x, v, par[0]]) + self.assertEqual(vsys.n_orig_sv, 2) + self.assertEqual(vsys.order, 1) + + vsys2 = copy(vsys) + + self.assertEqual(orig_sys, vsys2.sys[:2]) + self.assertEqual(vsys2.vargs, [x, v, par[0]]) + self.assertEqual(vsys2.n_orig_sv, 2) + self.assertEqual(vsys2.order, 1) + + vsys2 = deepcopy(vsys) + + self.assertEqual(orig_sys, vsys2.sys[:2]) + self.assertEqual(vsys2.vargs, [x, v, par[0]]) + self.assertEqual(vsys2.n_orig_sv, 2) + self.assertEqual(vsys2.order, 1) + + vsys2 = loads(dumps(vsys)) + + self.assertEqual(orig_sys, vsys2.sys[:2]) + self.assertEqual(vsys2.vargs, [x, v, par[0]]) + self.assertEqual(vsys2.n_orig_sv, 2) + self.assertEqual(vsys2.order, 1) From 010039ecbc128b2eda821794b48204fb00de894a Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Sun, 9 Jun 2024 23:08:59 +0200 Subject: [PATCH 07/23] Testing for the variational integrators. --- heyoka/CMakeLists.txt | 1 + heyoka/_test_var_integrator.py | 226 ++++++++++++++++++++++++++++ heyoka/expose_batch_integrators.cpp | 10 +- heyoka/test.py | 4 + 4 files changed, 236 insertions(+), 5 deletions(-) create mode 100644 heyoka/_test_var_integrator.py diff --git a/heyoka/CMakeLists.txt b/heyoka/CMakeLists.txt index 7f018229..699e7e49 100644 --- a/heyoka/CMakeLists.txt +++ b/heyoka/CMakeLists.txt @@ -31,6 +31,7 @@ set(HEYOKA_PY_PYTHON_FILES _test_elp2000.py _test_lagham.py _test_var_ode_sys.py + _test_var_integrator.py model/__init__.py callback/__init__.py ) diff --git a/heyoka/_test_var_integrator.py b/heyoka/_test_var_integrator.py new file mode 100644 index 00000000..6524b41a --- /dev/null +++ b/heyoka/_test_var_integrator.py @@ -0,0 +1,226 @@ +# Copyright 2020, 2021, 2022, 2023, 2024 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com) +# +# This file is part of the heyoka.py library. +# +# This Source Code Form is subject to the terms of the Mozilla +# Public License v. 2.0. If a copy of the MPL was not distributed +# with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + + +import unittest as _ut + + +class var_integrator_test_case(_ut.TestCase): + def test_scalar(self): + from . import ( + make_vars, + var_ode_sys, + var_args, + cos, + sin, + par, + time, + taylor_adaptive, + core, + ) + from sys import getrefcount + import numpy as np + + x, v = make_vars("x", "v") + + orig_sys = [(x, v), (v, cos(time) - par[0] * v - sin(x))] + + vsys = var_ode_sys(orig_sys, var_args.vars, order=2) + + ta = taylor_adaptive(vsys, [0.2, 0.3], pars=[0.4], time=0.5, compact_mode=True) + + self.assertTrue(ta.dim > 2) + self.assertEqual(ta.n_orig_sv, 2) + self.assertTrue(ta.is_variational) + self.assertEqual(ta.vargs, [x, v]) + + # Check that the refcount increases by when accessing tstate. + rc = getrefcount(ta) + ts = ta.tstate + self.assertEqual(getrefcount(ta), rc + 1) + self.assertEqual(ts.shape, (2,)) + self.assertTrue(np.all(ts == [0.0, 0.0])) + + self.assertEqual(ta.get_vslice(order=0), slice(0, 2, None)) + self.assertEqual(ta.get_vslice(order=0, component=1), slice(1, 2, None)) + + self.assertEqual(ta.get_vslice(order=1), slice(2, 6, None)) + self.assertEqual(ta.get_vslice(order=1, component=1), slice(4, 6, None)) + + self.assertEqual(ta.get_mindex(0), [0, 0, 0]) + self.assertEqual(ta.get_mindex(1), [1, 0, 0]) + self.assertEqual(ta.get_mindex(2), [0, 1, 0]) + self.assertEqual(ta.get_mindex(3), [0, 0, 1]) + self.assertEqual(ta.get_mindex(4), [1, 1, 0]) + self.assertEqual(ta.get_mindex(i=5), [1, 0, 1]) + + ta.propagate_until(3.0) + + ts2 = ta.eval_taylor_map([0.0, 0.0]) + self.assertEqual(getrefcount(ta), rc + 2) + self.assertTrue(np.shares_memory(ts, ts2)) + self.assertTrue(np.all(ts2 == ta.state[:2])) + + ts2 = ta.eval_taylor_map(np.array([0.0, 0.0])) + self.assertTrue(np.shares_memory(ts, ts2)) + self.assertTrue(np.all(ts2 == ta.state[:2])) + + with self.assertRaises(TypeError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0], dtype=np.int32)) + self.assertTrue( + "Invalid dtype detected for the inputs of a Taylor map evaluation:" + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0, 0.0, 0.0])[::2]) + self.assertTrue( + "Invalid inputs array detected in a Taylor map evaluation: the array is not C-style contiguous, please " + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([[0.0, 0.0], [0.0, 0.0]])) + self.assertTrue( + "The array of inputs provided for the evaluation of a Taylor map has 2 dimensions, " + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0, 0.0, 0.0])) + self.assertTrue( + "The array of inputs provided for the evaluation of a Taylor map has 4 elements, " + in str(cm.exception) + ) + + if not hasattr(core, "real"): + return + + from . import real + + prec = 14 + + ta = taylor_adaptive( + vsys, + [real(0.2, prec), real(0.3, prec)], + pars=[real(0.4, prec)], + time=real(0.5, prec), + compact_mode=True, + fp_type=real, + ) + + ta.propagate_until(real(3.0, prec)) + + with self.assertRaises(ValueError): + ta.eval_taylor_map(np.empty((2,), dtype=real)) + + with self.assertRaises(ValueError): + ta.eval_taylor_map(np.array([real(0.0, prec), real(0.0, prec - 1)])) + + ts2 = ta.eval_taylor_map(np.array([real(0.0, prec), real(0.0, prec)])) + self.assertTrue(np.all(ts2 == ta.state[:2])) + + def test_batch(self): + from . import ( + make_vars, + var_ode_sys, + var_args, + cos, + sin, + par, + time, + taylor_adaptive_batch, + core, + ) + from sys import getrefcount + import numpy as np + + x, v = make_vars("x", "v") + + orig_sys = [(x, v), (v, cos(time) - par[0] * v - sin(x))] + + vsys = var_ode_sys(orig_sys, var_args.vars, order=2) + + ta = taylor_adaptive_batch( + vsys, + [[0.2, 0.21], [0.3, 0.31]], + pars=[[0.4, 0.41]], + time=[0.5, 0.51], + compact_mode=True, + ) + + self.assertTrue(ta.dim > 2) + self.assertEqual(ta.n_orig_sv, 2) + self.assertTrue(ta.is_variational) + self.assertEqual(ta.vargs, [x, v]) + + # Check that the refcount increases by when accessing tstate. + rc = getrefcount(ta) + ts = ta.tstate + self.assertEqual(getrefcount(ta), rc + 1) + self.assertEqual(ts.shape, (2, 2)) + self.assertTrue(np.all(ts == [[0.0, 0.0], [0.0, 0.0]])) + + self.assertEqual(ta.get_vslice(order=0), slice(0, 2, None)) + self.assertEqual(ta.get_vslice(order=0, component=1), slice(1, 2, None)) + + self.assertEqual(ta.get_vslice(order=1), slice(2, 6, None)) + self.assertEqual(ta.get_vslice(order=1, component=1), slice(4, 6, None)) + + self.assertEqual(ta.get_mindex(0), [0, 0, 0]) + self.assertEqual(ta.get_mindex(1), [1, 0, 0]) + self.assertEqual(ta.get_mindex(2), [0, 1, 0]) + self.assertEqual(ta.get_mindex(3), [0, 0, 1]) + self.assertEqual(ta.get_mindex(4), [1, 1, 0]) + self.assertEqual(ta.get_mindex(i=5), [1, 0, 1]) + + ta.propagate_until(3.0) + + ts2 = ta.eval_taylor_map([[0.0, 0.0], [0.0, 0.0]]) + self.assertEqual(getrefcount(ta), rc + 2) + self.assertTrue(np.shares_memory(ts, ts2)) + self.assertTrue(np.all(ts2 == ta.state[:2])) + + ts2 = ta.eval_taylor_map(np.array([[0.0, 0.0], [0.0, 0.0]])) + self.assertTrue(np.shares_memory(ts, ts2)) + self.assertTrue(np.all(ts2 == ta.state[:2])) + + with self.assertRaises(TypeError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0], dtype=np.int32)) + self.assertTrue( + "Invalid dtype detected for the inputs of a Taylor map evaluation:" + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0, 0.0, 0.0])[::2]) + self.assertTrue( + "Invalid inputs array detected in a Taylor map evaluation: the array is not C-style contiguous, please " + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([0.0, 0.0])) + self.assertTrue( + "The array of inputs provided for the evaluation of a Taylor map has 1 dimension(s), " + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([[0.0, 0.0]])) + self.assertTrue( + "The array of inputs provided for the evaluation of a Taylor map has 1 row(s), but it must have 2 row(s) instead" + in str(cm.exception) + ) + + with self.assertRaises(ValueError) as cm: + ta.eval_taylor_map(np.array([[0.0], [0.0]])) + self.assertTrue( + "The array of inputs provided for the evaluation of a Taylor map has 1 column(s), but it must have 2 column(s) instead" + in str(cm.exception) + ) diff --git a/heyoka/expose_batch_integrators.cpp b/heyoka/expose_batch_integrators.cpp index 8ff14cde..2fcb6000 100644 --- a/heyoka/expose_batch_integrators.cpp +++ b/heyoka/expose_batch_integrators.cpp @@ -606,7 +606,7 @@ void expose_batch_integrator_impl(py::module_ &m, const std::string &suffix) // NOTE: this needs to be done regardless of the original type of in. if (inputs.ndim() != 2) [[unlikely]] { py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " - "of a Taylor map has {} dimensions, " + "of a Taylor map has {} dimension(s), " "but it must have 2 dimensions instead", inputs.ndim()) .c_str()); @@ -615,16 +615,16 @@ void expose_batch_integrator_impl(py::module_ &m, const std::string &suffix) // Validate the shape for the inputs. if (boost::numeric_cast(inputs.shape(0)) != ta->get_n_orig_sv()) [[unlikely]] { py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " - "of a Taylor map has {} rows, " - "but it must have {} rows instead", + "of a Taylor map has {} row(s), " + "but it must have {} row(s) instead", inputs.shape(0), ta->get_n_orig_sv()) .c_str()); } if (boost::numeric_cast(inputs.shape(1)) != ta->get_batch_size()) [[unlikely]] { py_throw(PyExc_ValueError, fmt::format("The array of inputs provided for the evaluation " - "of a Taylor map has {} columns, " - "but it must have {} columns instead", + "of a Taylor map has {} column(s), " + "but it must have {} column(s) instead", inputs.shape(1), ta->get_batch_size()) .c_str()); } diff --git a/heyoka/test.py b/heyoka/test.py index 996b5cae..c263dfc5 100644 --- a/heyoka/test.py +++ b/heyoka/test.py @@ -2080,6 +2080,7 @@ def run_test_suite(): _test_elp2000, _test_lagham, _test_var_ode_sys, + _test_var_integrator, ) import numpy as np from .model import nbody @@ -2100,6 +2101,9 @@ def run_test_suite(): tl.loadTestsFromTestCase(_test_batch_integrator.batch_integrator_test_case) ) suite.addTest(tl.loadTestsFromTestCase(_test_var_ode_sys.var_ode_sys_test_case)) + suite.addTest( + tl.loadTestsFromTestCase(_test_var_integrator.var_integrator_test_case) + ) suite.addTest(tl.loadTestsFromTestCase(_test_lagham.lagham_test_case)) suite.addTest(tl.loadTestsFromTestCase(_test_vsop2013.vsop2013_test_case)) suite.addTest(tl.loadTestsFromTestCase(_test_elp2000.elp2000_test_case)) From e94a08e1d2f1fdcd1feb227db650e947de1f128f Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Mon, 10 Jun 2024 08:24:04 +0200 Subject: [PATCH 08/23] Some testing for dfun(). --- heyoka/_test_expression.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/heyoka/_test_expression.py b/heyoka/_test_expression.py index 90b99cb4..7d677e1c 100644 --- a/heyoka/_test_expression.py +++ b/heyoka/_test_expression.py @@ -342,3 +342,14 @@ def test_relu_wrappers(self): self.assertEqual(leaky_relup(0.0)(x * y), relup(x * y)) self.assertEqual(leaky_relu(0.1)(x * y + y), relu(x * y + y, 0.1)) self.assertEqual(leaky_relup(0.1)(x * y + y), relup(x * y + y, 0.1)) + + def test_dfun(self): + from . import make_vars, dfun + + x, y = make_vars("x", "y") + + self.assertEqual(str(dfun("f", [x, y])), "(∂^0 f)") + self.assertEqual( + str(dfun(name="f", args=[x, y], didx=[(0, 3), (1, 5)])), + "(∂^8 f)/(∂a0^3 ∂a1^5)", + ) From 2b39559529c19605eeef092dab17f0198dbbf1f9 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Mon, 10 Jun 2024 08:51:00 +0200 Subject: [PATCH 09/23] Small doc tweak. --- doc/notebooks/Batch mode overview.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/notebooks/Batch mode overview.ipynb b/doc/notebooks/Batch mode overview.ipynb index e80155d3..3773df22 100644 --- a/doc/notebooks/Batch mode overview.ipynb +++ b/doc/notebooks/Batch mode overview.ipynb @@ -22,10 +22,10 @@ "It is important to emphasise that batch mode does not reduce\n", "the CPU time required to integrate a system of ODEs. Rather, as a fine-grained\n", "form of data parallelism, batch mode allows to integrate multiple ODE systems in parallel\n", - "at no additional cost, and it is thus most useful when the need arise\n", + "at (almost) no additional cost, and it is thus most useful when the need arise\n", "to integrate the same ODE system with different initial conditions and parameters.\n", "\n", - "Although batch mode can in principle be used with all floating-point types supported\n", + "Although batch mode can in principle be used with all the fundamental C++ floating-point types supported\n", "by heyoka.py, in practice at this time no CPU provides SIMD instructions for extended-precision\n", "datatypes. Thus, heyoka.py currently enables batch mode only for single and double precision\n", "computations. Specifically, in this tutorial we will consider the application of batch mode to\n", From 88797e8826bb0bed817bbad10059d211d18162fc Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Mon, 10 Jun 2024 08:51:18 +0200 Subject: [PATCH 10/23] Expose the vorder property. --- heyoka/expose_batch_integrators.cpp | 1 + heyoka/taylor_expose_integrator.cpp | 1 + 2 files changed, 2 insertions(+) diff --git a/heyoka/expose_batch_integrators.cpp b/heyoka/expose_batch_integrators.cpp index 2fcb6000..b4c28770 100644 --- a/heyoka/expose_batch_integrators.cpp +++ b/heyoka/expose_batch_integrators.cpp @@ -544,6 +544,7 @@ void expose_batch_integrator_impl(py::module_ &m, const std::string &suffix) .def_property_readonly("n_orig_sv", &hey::taylor_adaptive_batch::get_n_orig_sv) .def_property_readonly("is_variational", &hey::taylor_adaptive_batch::is_variational) .def_property_readonly("vargs", &hey::taylor_adaptive_batch::get_vargs) + .def_property_readonly("vorder", &hey::taylor_adaptive_batch::get_vorder) .def_property_readonly("tstate", &fetch_tstate) .def( "get_vslice", diff --git a/heyoka/taylor_expose_integrator.cpp b/heyoka/taylor_expose_integrator.cpp index aa7fb269..a92af47d 100644 --- a/heyoka/taylor_expose_integrator.cpp +++ b/heyoka/taylor_expose_integrator.cpp @@ -378,6 +378,7 @@ void expose_taylor_integrator_impl(py::module &m, const std::string &suffix) .def_property_readonly("n_orig_sv", &hey::taylor_adaptive::get_n_orig_sv) .def_property_readonly("is_variational", &hey::taylor_adaptive::is_variational) .def_property_readonly("vargs", &hey::taylor_adaptive::get_vargs) + .def_property_readonly("vorder", &hey::taylor_adaptive::get_vorder) .def_property_readonly("tstate", &fetch_tstate) .def( "get_vslice", From 659731d588ac0ca5c3ce28909d78064822f2b5c8 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Mon, 10 Jun 2024 08:52:20 +0200 Subject: [PATCH 11/23] Testing. --- heyoka/_test_var_integrator.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/heyoka/_test_var_integrator.py b/heyoka/_test_var_integrator.py index 6524b41a..bb6a7376 100644 --- a/heyoka/_test_var_integrator.py +++ b/heyoka/_test_var_integrator.py @@ -37,6 +37,7 @@ def test_scalar(self): self.assertTrue(ta.dim > 2) self.assertEqual(ta.n_orig_sv, 2) self.assertTrue(ta.is_variational) + self.assertEqual(ta.vorder, 2) self.assertEqual(ta.vargs, [x, v]) # Check that the refcount increases by when accessing tstate. @@ -157,6 +158,7 @@ def test_batch(self): self.assertTrue(ta.dim > 2) self.assertEqual(ta.n_orig_sv, 2) self.assertTrue(ta.is_variational) + self.assertEqual(ta.vorder, 2) self.assertEqual(ta.vargs, [x, v]) # Check that the refcount increases by when accessing tstate. From 874b6b703ef1405c551557a571068faea414b56b Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Mon, 10 Jun 2024 11:31:46 +0200 Subject: [PATCH 12/23] Initial work on the new tutorial. --- doc/notebooks/var_ode_sys.ipynb | 157 ++++++++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100644 doc/notebooks/var_ode_sys.ipynb diff --git a/doc/notebooks/var_ode_sys.ipynb b/doc/notebooks/var_ode_sys.ipynb new file mode 100644 index 00000000..531dee32 --- /dev/null +++ b/doc/notebooks/var_ode_sys.ipynb @@ -0,0 +1,157 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e916542d-c443-4403-bd35-e56122b2823b", + "metadata": {}, + "source": [ + "# Variational ODEs\n", + "\n", + "```{versionadded} 5.0.0\n", + "\n", + "```\n", + "\n", + "Consider a system of differential equations in the standard form\n", + "\n", + "$$\n", + "\\frac{d\\boldsymbol{x}}{dt} = f\\left(\\boldsymbol{x}, \\boldsymbol{\\alpha}, t\\right),\n", + "$$\n", + "\n", + "where $\\boldsymbol{x}$ is the vector of state variables and $\\boldsymbol{\\alpha}$ a vector of parameters. For a given set of initial conditions $\\boldsymbol{x}_0$ at time $t_0$, the solution of this system will be\n", + "\n", + "$$\n", + "\\boldsymbol{x} = \\boldsymbol{x}\\left(t, \\boldsymbol{x}_0, t_0, \\boldsymbol{\\alpha} \\right).\n", + "$$\n", + "\n", + "When solving numerically initial-value problems, it is often useful to compute not only the solution, but also its partial derivatives with respect to the initial conditions and/or the parameters. The derivatives with respect to the initial conditions, for instance, are needed for the computation of [chaos indicators](https://en.wikipedia.org/wiki/Lyapunov_exponent) and for [uncertainty propagation](https://en.wikipedia.org/wiki/Propagation_of_uncertainty), and they can also be used to propagate a small neighborhood in phase space around the initial conditions. The derivatives with respect to the parameters of the system are required when formulating optimisation and inversion problems such as orbit determination, trajectory optimisation and training of neural networks in [neural ODEs](./NeuralODEs.ipynb).\n", + "\n", + "There are two main methods for the computation of the partial derivatives. The first one is based on the application of automatic differentiation (AD) techniques directly to the numerical integration algorithm. This can be done either by replacing the algebra of floating-point numbers with the algebra of (generalised) [dual numbers](https://en.wikipedia.org/wiki/Dual_number) (aka truncated Taylor polynomials), or via [differentiable programming](https://en.wikipedia.org/wiki/Differentiable_programming) techniques. The former approach is used by libraries such as [pyaudi](https://github.com/darioizzo/audi), [desolver](https://github.com/Microno95/desolver) and [TaylorIntegration.jl](https://docs.sciml.ai/TaylorIntegration/stable/jet_transport/), while differentiable programming is popular in the [Julia programming language](https://en.wikipedia.org/wiki/Julia_(programming_language)) community.\n", + "\n", + "The second method is based on the formulation of the *variational equations*, that is, differential equations satisfied by the partial derivatives which are added to and solved together with the original ODEs. For instance, we can formulate differential equations for the first-order derivatives with respect to the initial conditions via elementary calculus:\n", + "\n", + "$$\n", + "\\frac{d}{dt}\\frac{\\partial \\boldsymbol{x}}{\\partial \\boldsymbol{x}_0} = \\frac{\\partial }{\\partial \\boldsymbol{x}_0} \\frac{d \\boldsymbol{x}}{dt} = \\frac{\\partial f}{\\partial \\boldsymbol{x}_0} = \\frac{\\partial f}{\\partial \\boldsymbol{x}} \\frac{\\partial \\boldsymbol{x}}{\\partial \\boldsymbol{x}_0}.\n", + "$$\n", + "\n", + "The variational ODE system then reads\n", + "\n", + "$$\n", + "\\begin{cases}\n", + "\\frac{d\\boldsymbol{x}}{dt} = f\\left(\\boldsymbol{x}, \\boldsymbol{\\alpha}, t\\right) \\\\\n", + "\\frac{d}{dt}\\frac{\\partial \\boldsymbol{x}}{\\partial \\boldsymbol{x}_0} = \\frac{\\partial f}{\\partial \\boldsymbol{x}} \\frac{\\partial \\boldsymbol{x}}{\\partial \\boldsymbol{x}_0}\n", + "\\end{cases},\n", + "$$\n", + "\n", + "and the original state vector $\\boldsymbol{x}$ has been extended to include the variational state variables $\\frac{\\partial \\boldsymbol{x}}{\\partial \\boldsymbol{x}_0}$.\n", + "\n", + "heyoka.py adopts the variational approach for the computation of the partial derivatives, supporting the formulation of variational ODEs at arbitrary differentiation orders and with respect to any combination of initial conditions, parameters and initial time. In this tutorial, we will explore this feature and show a couple of interesting use cases.\n", + "\n", + "Before beginning, however, let us point out for clarity (and for the benefit of the search engines indexing this page) that in the scientific literature there is a bewildering variety of different names and monikers used when discussing partial derivatives of ODEs and their applications. Here is a (partial) list:\n", + "\n", + "- in the astrodynamics community, the term *differential algebra* is often used to refer to the computation of partial derivatives via truncated Taylor polynomials (e.g., see [this paper](https://link.springer.com/article/10.1007/s10569-010-9283-5)). The term actually originates from the community of beam physics, where it has been used in the context of the theoretical modelling of particle accelerators since the 90s (e.g., see [this review](https://www.bmtdynamics.org/pub/papers/DAHAPE12/DAHAPE12.pdf));\n", + "- in the mathematical community, the term *jet transport* is sometimes used to refer to the propagation of a small neighborhood in phase space around the initial conditions via the Taylor series constructed form the partial derivatives (e.g., see [this paper](http://www.maia.ub.es/~angel/varis/granada09.pdf)). In heyoka.py, we refer to a similar idea as {ref}`Taylor map evaluation `;\n", + "- in the Julia programming language community, the term *local sensitivity analysis* refers to the computation of the partial derivatives via the variational equations, while *discrete sensitivity analysis* refers to the computation of the partial derivatives by directly differentiating the numerical method's steps (e.g., see [this review](https://arxiv.org/abs/1812.01892));\n", + "- in the space engineering community, the term *state transition tensors* is sometimes used to indicate the generalisations of the [state transition matrix](https://en.wikipedia.org/wiki/State-transition_matrix) (which in turn is built from the first-order partial derivatives) to higher differentiation orders.\n", + "\n", + "## Constructing a variational ODE system\n", + "\n", + "Let us begin with the construction of a simple ODE system:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "2c186e41-5526-4b7d-a11a-0ce38332c9ea", + "metadata": {}, + "outputs": [], + "source": [ + "import heyoka as hy\n", + "\n", + "# Create the symbolic variables x and v.\n", + "x, v = hy.make_vars(\"x\", \"v\")\n", + "\n", + "# Create an ODE system.\n", + "sys = [(x, v), (v, hy.cos(hy.time) - hy.par[0] * v - hy.sin(x))]" + ] + }, + { + "cell_type": "markdown", + "id": "1d9126e2-93b3-49a2-8fdb-33e31a161e35", + "metadata": {}, + "source": [ + "This is the forced damped pendulum system already considered in [another tutorial](<./Non-autonomous systems.ipynb>), where we have introduced the air friction coefficient as the [runtime parameter](<./ODEs with parameters.ipynb>) ``par[0]``.\n", + "\n", + "In order to formulate the variational equations in heyoka.py, we use the " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "ba980238-26c2-417c-817f-92fb2dbe3ec1", + "metadata": {}, + "outputs": [], + "source": [ + "vsys = hy.var_ode_sys(sys, hy.var_args.vars, order=2)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "e4715f11-f4c3-4953-8d6e-7d01cd14d32a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vsys.order" + ] + }, + { + "cell_type": "markdown", + "id": "d2ae2f2e-3b67-47be-a863-b17f44fb26c2", + "metadata": {}, + "source": [ + "(taylor_map)=\n", + "## Taylor map evaluation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3581bf58-ec4c-4d96-8f91-54f58d73c7ce", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 8c0a02c5160e7fc1079af25b48edb40dccdc4c2f Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Mon, 10 Jun 2024 13:41:28 +0200 Subject: [PATCH 13/23] More docs/tutorials WIP. --- doc/api_exsys.rst | 1 + doc/api_var_ode_sys.rst | 24 +++ doc/basic_tutorials.rst | 1 + doc/index.md | 1 + doc/notebooks/The variational equations.ipynb | 2 +- doc/notebooks/var_ode_sys.ipynb | 49 +++++- heyoka/__init__.py | 16 ++ heyoka/docstrings.cpp | 143 ++++++++++++++++++ heyoka/docstrings.hpp | 15 +- heyoka/expose_expression.cpp | 2 +- heyoka/expose_var_ode_sys.cpp | 24 +-- 11 files changed, 258 insertions(+), 20 deletions(-) create mode 100644 doc/api_var_ode_sys.rst diff --git a/doc/api_exsys.rst b/doc/api_exsys.rst index faed6fc3..6cc31324 100644 --- a/doc/api_exsys.rst +++ b/doc/api_exsys.rst @@ -34,6 +34,7 @@ Attributes :toctree: autosummary_generated par + time Enums ----- diff --git a/doc/api_var_ode_sys.rst b/doc/api_var_ode_sys.rst new file mode 100644 index 00000000..0aa44d2e --- /dev/null +++ b/doc/api_var_ode_sys.rst @@ -0,0 +1,24 @@ +.. _api_var_ode_sys: + +Variational ODE systems +======================= + +.. currentmodule:: heyoka + +Classes +------- + +.. autosummary:: + :toctree: autosummary_generated + :template: custom-class-template.rst + + var_ode_sys + +Enums +----- + +.. autosummary:: + :toctree: autosummary_generated + :template: custom-enum-template.rst + + var_args diff --git a/doc/basic_tutorials.rst b/doc/basic_tutorials.rst index 8f9a82ff..7f598ee5 100644 --- a/doc/basic_tutorials.rst +++ b/doc/basic_tutorials.rst @@ -34,3 +34,4 @@ Basic notebooks/Non-autonomous systems.ipynb notebooks/Dense output.ipynb notebooks/Event detection + notebooks/var_ode_sys.ipynb diff --git a/doc/index.md b/doc/index.md index c1d1cc90..0f3261fd 100644 --- a/doc/index.md +++ b/doc/index.md @@ -107,6 +107,7 @@ examples_others api_exsys api_integrators +api_var_ode_sys api_lagham api_model ``` diff --git a/doc/notebooks/The variational equations.ipynb b/doc/notebooks/The variational equations.ipynb index 63b0966e..2bcf22fe 100644 --- a/doc/notebooks/The variational equations.ipynb +++ b/doc/notebooks/The variational equations.ipynb @@ -292,7 +292,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.13" + "version": "3.10.13" } }, "nbformat": 4, diff --git a/doc/notebooks/var_ode_sys.ipynb b/doc/notebooks/var_ode_sys.ipynb index 531dee32..065c13d8 100644 --- a/doc/notebooks/var_ode_sys.ipynb +++ b/doc/notebooks/var_ode_sys.ipynb @@ -5,6 +5,8 @@ "id": "e916542d-c443-4403-bd35-e56122b2823b", "metadata": {}, "source": [ + "(var_ode_sys)=\n", + "\n", "# Variational ODEs\n", "\n", "```{versionadded} 5.0.0\n", @@ -55,12 +57,12 @@ "\n", "## Constructing a variational ODE system\n", "\n", - "Let us begin with the construction of a simple ODE system:" + "Let us begin with the definition of a simple ODE system:" ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 1, "id": "2c186e41-5526-4b7d-a11a-0ce38332c9ea", "metadata": {}, "outputs": [], @@ -81,7 +83,7 @@ "source": [ "This is the forced damped pendulum system already considered in [another tutorial](<./Non-autonomous systems.ipynb>), where we have introduced the air friction coefficient as the [runtime parameter](<./ODEs with parameters.ipynb>) ``par[0]``.\n", "\n", - "In order to formulate the variational equations in heyoka.py, we use the " + "We then proceed to create a {class}`~heyoka.var_ode_sys`" ] }, { @@ -96,7 +98,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 3, "id": "e4715f11-f4c3-4953-8d6e-7d01cd14d32a", "metadata": {}, "outputs": [ @@ -106,7 +108,7 @@ "2" ] }, - "execution_count": 7, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" } @@ -115,6 +117,41 @@ "vsys.order" ] }, + { + "cell_type": "code", + "execution_count": 12, + "id": "16e3bc0a-e7d5-499c-88bd-e17aaa6c1657", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[(x, v),\n", + " (v, ((cos(t) - (p0 * v)) - sin(x))),\n", + " (∂[(0, 1)]x, ∂[(0, 1)]v),\n", + " (∂[(1, 1)]x, ∂[(1, 1)]v),\n", + " (∂[(0, 1)]v, (-(∂[(0, 1)]x * cos(x)) - (∂[(0, 1)]v * p0))),\n", + " (∂[(1, 1)]v, (-(∂[(1, 1)]x * cos(x)) - (∂[(1, 1)]v * p0))),\n", + " (∂[(0, 2)]x, ∂[(0, 2)]v),\n", + " (∂[(0, 1), (1, 1)]x, ∂[(0, 1), (1, 1)]v),\n", + " (∂[(1, 2)]x, ∂[(1, 2)]v),\n", + " (∂[(0, 2)]v,\n", + " (-(∂[(0, 2)]v * p0) - (((∂[(0, 1)]x * -sin(x)) * ∂[(0, 1)]x) + (∂[(0, 2)]x * cos(x))))),\n", + " (∂[(0, 1), (1, 1)]v,\n", + " (-(∂[(0, 1), (1, 1)]v * p0) - (((∂[(0, 1)]x * -sin(x)) * ∂[(1, 1)]x) + (∂[(0, 1), (1, 1)]x * cos(x))))),\n", + " (∂[(1, 2)]v,\n", + " (-(∂[(1, 2)]v * p0) - (((∂[(1, 1)]x * -sin(x)) * ∂[(1, 1)]x) + (∂[(1, 2)]x * cos(x)))))]" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vsys.sys" + ] + }, { "cell_type": "markdown", "id": "d2ae2f2e-3b67-47be-a863-b17f44fb26c2", @@ -127,7 +164,7 @@ { "cell_type": "code", "execution_count": null, - "id": "3581bf58-ec4c-4d96-8f91-54f58d73c7ce", + "id": "40d0734e-25b3-4bc5-9f2b-13007790e3a0", "metadata": {}, "outputs": [], "source": [] diff --git a/heyoka/__init__.py b/heyoka/__init__.py index 4c724a15..72da6d02 100644 --- a/heyoka/__init__.py +++ b/heyoka/__init__.py @@ -382,3 +382,19 @@ def _create_par(): >>> p0 = par[0] # p0 will represent the parameter value at index 0 """ + +# Machinery for the time attribute. +def _create_time(): + from . import core + + return core._time + + +time = _create_time() +""" +Time placeholder. + +This global object is an :py:class:`~heyoka.expression` which is used to represent +time (i.e., the independent variable) in differential equations. + +""" diff --git a/heyoka/docstrings.cpp b/heyoka/docstrings.cpp index 43bf1a8c..6162293d 100644 --- a/heyoka/docstrings.cpp +++ b/heyoka/docstrings.cpp @@ -535,4 +535,147 @@ A few checks are run on the input arguments. Specifically, the number of geodesi )"; } +std::string var_ode_sys() +{ + return R"(Class to represent variational ODE systems. + +.. note:: + + A :ref:`tutorial ` explaining the use of this + class is available. + +)"; +} + +std::string var_ode_sys_init() +{ + return R"(__init__(self, sys: list[tuple[expression, expression]], args: var_args | list[expression], order: int = 1) + +Constructor. + +A variational ODE system is constructed from two mandatory arguments: the original ODE system +*sys* and an *args* parameter representing the quantities with respect to which the variational +equations will be formulated. + +If *args* is of type :class:`~heyoka.var_args`, then the variational equations will be constructed +with respect to arguments deduced from *sys*. E.g., if *sys* contains the two state variables +:math:`x` and :math:`y` and *args* is the *vars* enumerator of :class:`~heyoka.var_args`, then +the variational equations will be formuated with respect to the initial conditions of +:math:`x` and :math:`y`. Similarly, if *sys* contains two parameters ``par[0]`` and ``par[1]`` +and *args* is the *params* enumerator of :class:`~heyoka.var_args`, then +the variational equations will be formuated with respect to the two parameters. + +If *args* is a list of :class:`~heyoka.expression`, then the variational equations will be formulated +with respect to the quantities contained in the list. Specifically: + +- variable expression are used to request derivatives with respect to the intial conditions + for that state variable; +- parameter expressions are used to request derivatives with respect to those parameters; +- :attr:`heyoka.time` is used to request the derivative with respect to the initial integration time. + +Several checks are run on the input arguments. Specifically: + +- *sys* must be a valid ODE system, +- if *args* is of type :class:`~heyoka.var_args`, it must be either one of the valid enumerators + or a combination of the valid enumerators, +- if *args* is a list of expressions: + + - it must not be empty, + - it must consists only of variables, parameters or the :attr:`heyoka.time` placeholder, + - it must not contain duplicates, + - any variable expression must refer to a state variable in *sys*. + +Additionally, the differentiation order *order* must be at least 1. + +:param sys: the input ODE system. +:param args: the variational arguments. +:param order: the differentiation order. + +:raises ValueError: if one or more input arguments are malformed, as explained above. + +)"; +} + +std::string var_ode_sys_sys() +{ + return R"(The full system of equations (including partials). + +:rtype: list[tuple[expression, expression]] + +)"; +} + +std::string var_ode_sys_vargs() +{ + return R"(The list of variational arguments. + +:rtype: list[expression] + +)"; +} + +std::string var_ode_sys_n_orig_sv() +{ + return R"(The number of original state variables. + +:rtype: int + +)"; +} + +std::string var_ode_sys_order() +{ + return R"(The differentitation order. + +:rtype: int + +)"; +} + +std::string var_args() +{ + return R"(Enum for selecting variational arguments. + +Values of this enum can be used in the constructor of :class:`~heyoka.var_ode_sys()` to select +the arguments with respect to which the variational equations will be formulated. + +The enumerators can be combined with the logical OR ``|`` operator. + +Examples: + >>> from heyoka import var_args + >>> va = var_args.vars | var_args.params # Select differentiation wrt all initial conditions and parameters + >>> var_args.vars | var_args.params | var_args.time == var_args.all + True + +)"; +} + +std::string var_args_vars() +{ + return R"(Differentiate with respect to the initial conditions of all state variables. + +)"; +} + +std::string var_args_params() +{ + return R"(Differentiate with respect to all runtime parameters. + +)"; +} + +std::string var_args_time() +{ + return R"(Differentiate with respect to the initial integration time. + +)"; +} + +std::string var_args_all() +{ + return R"(Differentiate with respect to the initial conditions of all state variables, all runtime parameters and the initial integration time. + +)"; +} + } // namespace heyoka_py::docstrings diff --git a/heyoka/docstrings.hpp b/heyoka/docstrings.hpp index c455226b..67cb6c17 100644 --- a/heyoka/docstrings.hpp +++ b/heyoka/docstrings.hpp @@ -44,11 +44,24 @@ std::string diff_args_all(); std::string lagrangian(); std::string hamiltonian(); -// Models +// Models. std::string cart2geo(); std::string nrlmsise00_tn(); std::string jb08_tn(); +// var_ode_sys() and related. +std::string var_args(); +std::string var_args_vars(); +std::string var_args_params(); +std::string var_args_time(); +std::string var_args_all(); +std::string var_ode_sys(); +std::string var_ode_sys_init(); +std::string var_ode_sys_sys(); +std::string var_ode_sys_vargs(); +std::string var_ode_sys_n_orig_sv(); +std::string var_ode_sys_order(); + } // namespace heyoka_py::docstrings #endif diff --git a/heyoka/expose_expression.cpp b/heyoka/expose_expression.cpp index 6ce2811f..049c8649 100644 --- a/heyoka/expose_expression.cpp +++ b/heyoka/expose_expression.cpp @@ -448,7 +448,7 @@ void expose_expression(py::module_ &m) "name"_a, "args"_a, "didx"_a = py::none{}); // Time. - m.attr("time") = hey::time; + m.attr("_time") = hey::time; // pi. m.attr("pi") = hey::pi; diff --git a/heyoka/expose_var_ode_sys.cpp b/heyoka/expose_var_ode_sys.cpp index ed3d231d..60f823ba 100644 --- a/heyoka/expose_var_ode_sys.cpp +++ b/heyoka/expose_var_ode_sys.cpp @@ -16,6 +16,7 @@ #include #include "common_utils.hpp" +#include "docstrings.hpp" #include "expose_var_ode_sys.hpp" #include "pickle_wrappers.hpp" @@ -31,24 +32,25 @@ void expose_var_ode_sys(py::module_ &m) using namespace pybind11::literals; // var_args enum. - py::enum_(m, "var_args", py::arithmetic()) - .value("vars", hey::var_args::vars) - .value("params", hey::var_args::params) - .value("time", hey::var_args::time) - .value("all", hey::var_args::all) + py::enum_(m, "var_args", docstrings::var_args().c_str()) + .value("vars", hey::var_args::vars, docstrings::var_args_vars().c_str()) + .value("params", hey::var_args::params, docstrings::var_args_params().c_str()) + .value("time", hey::var_args::time, docstrings::var_args_time().c_str()) + .value("all", hey::var_args::all, docstrings::var_args_all().c_str()) .def("__or__", [](hey::var_args a, hey::var_args b) { return a | b; }) .def("__and__", [](hey::var_args a, hey::var_args b) { return a & b; }); // var_ode_sys class. - py::class_ v_cl(m, "var_ode_sys", py::dynamic_attr{}); + py::class_ v_cl(m, "var_ode_sys", py::dynamic_attr{}, docstrings::var_ode_sys().c_str()); v_cl.def(py::init([](const std::vector> &sys, const std::variant> &args, std::uint32_t order) { return hey::var_ode_sys(sys, args, order); }), - "sys"_a, "args"_a, "order"_a = static_cast(1)); - v_cl.def_property_readonly("sys", &hey::var_ode_sys::get_sys); - v_cl.def_property_readonly("vargs", &hey::var_ode_sys::get_vargs); - v_cl.def_property_readonly("n_orig_sv", &hey::var_ode_sys::get_n_orig_sv); - v_cl.def_property_readonly("order", &hey::var_ode_sys::get_order); + "sys"_a, "args"_a, "order"_a = static_cast(1), docstrings::var_ode_sys_init().c_str()); + v_cl.def_property_readonly("sys", &hey::var_ode_sys::get_sys, docstrings::var_ode_sys_sys().c_str()); + v_cl.def_property_readonly("vargs", &hey::var_ode_sys::get_vargs, docstrings::var_ode_sys_vargs().c_str()); + v_cl.def_property_readonly("n_orig_sv", &hey::var_ode_sys::get_n_orig_sv, + docstrings::var_ode_sys_n_orig_sv().c_str()); + v_cl.def_property_readonly("order", &hey::var_ode_sys::get_order, docstrings::var_ode_sys_order().c_str()); // Copy/deepcopy. v_cl.def("__copy__", copy_wrapper); v_cl.def("__deepcopy__", deepcopy_wrapper, "memo"_a); From ca3ec0706e3210743ee4d77674167fcd067a6952 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Mon, 10 Jun 2024 17:33:34 +0200 Subject: [PATCH 14/23] Typos. --- heyoka/docstrings.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/heyoka/docstrings.cpp b/heyoka/docstrings.cpp index 6162293d..128de087 100644 --- a/heyoka/docstrings.cpp +++ b/heyoka/docstrings.cpp @@ -560,10 +560,10 @@ equations will be formulated. If *args* is of type :class:`~heyoka.var_args`, then the variational equations will be constructed with respect to arguments deduced from *sys*. E.g., if *sys* contains the two state variables :math:`x` and :math:`y` and *args* is the *vars* enumerator of :class:`~heyoka.var_args`, then -the variational equations will be formuated with respect to the initial conditions of +the variational equations will be formulated with respect to the initial conditions of :math:`x` and :math:`y`. Similarly, if *sys* contains two parameters ``par[0]`` and ``par[1]`` and *args* is the *params* enumerator of :class:`~heyoka.var_args`, then -the variational equations will be formuated with respect to the two parameters. +the variational equations will be formulated with respect to the two parameters. If *args* is a list of :class:`~heyoka.expression`, then the variational equations will be formulated with respect to the quantities contained in the list. Specifically: @@ -575,9 +575,9 @@ with respect to the quantities contained in the list. Specifically: Several checks are run on the input arguments. Specifically: -- *sys* must be a valid ODE system, +- *sys* must be a valid ODE system; - if *args* is of type :class:`~heyoka.var_args`, it must be either one of the valid enumerators - or a combination of the valid enumerators, + or a combination of the valid enumerators; - if *args* is a list of expressions: - it must not be empty, From 9c7dce15fee77b4db0c2eed13634e312fd95e205 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Tue, 11 Jun 2024 10:19:51 +0200 Subject: [PATCH 15/23] Finish up the first version of the variational ODEs tutorial. --- .../Customising the adaptive integrator.ipynb | 4 +- doc/notebooks/var_ode_sys.ipynb | 633 +++++++++++++++++- 2 files changed, 615 insertions(+), 22 deletions(-) diff --git a/doc/notebooks/Customising the adaptive integrator.ipynb b/doc/notebooks/Customising the adaptive integrator.ipynb index 178cf959..e6a0a51e 100644 --- a/doc/notebooks/Customising the adaptive integrator.ipynb +++ b/doc/notebooks/Customising the adaptive integrator.ipynb @@ -128,6 +128,8 @@ "id": "national-investing", "metadata": {}, "source": [ + "(ta_compact_mode)=\n", + "\n", "Compact mode\n", "------------\n", "\n", @@ -290,7 +292,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.10.13" } }, "nbformat": 4, diff --git a/doc/notebooks/var_ode_sys.ipynb b/doc/notebooks/var_ode_sys.ipynb index 065c13d8..697a9dbd 100644 --- a/doc/notebooks/var_ode_sys.ipynb +++ b/doc/notebooks/var_ode_sys.ipynb @@ -83,7 +83,7 @@ "source": [ "This is the forced damped pendulum system already considered in [another tutorial](<./Non-autonomous systems.ipynb>), where we have introduced the air friction coefficient as the [runtime parameter](<./ODEs with parameters.ipynb>) ``par[0]``.\n", "\n", - "We then proceed to create a {class}`~heyoka.var_ode_sys`" + "We then proceed to create a {class}`~heyoka.var_ode_sys`:" ] }, { @@ -93,33 +93,41 @@ "metadata": {}, "outputs": [], "source": [ + "# Create the variational ODE system.\n", "vsys = hy.var_ode_sys(sys, hy.var_args.vars, order=2)" ] }, + { + "cell_type": "markdown", + "id": "bab348cb-ca8b-49f0-bc40-0c44ff091521", + "metadata": {}, + "source": [ + "Upon construction, {class}`~heyoka.var_ode_sys` formulates the variational equations for the input ODE system ``sys`` up to the specified differentiation ``order``. The second argument specifies with respect to which quantities the variational equations are formulated. In this case, we used the ``vars`` enumerator of the {class}`~heyoka.var_args` enum: this means that the variational equations will be formulated with respect to the initial conditions of the state variables. In a completely equivalent manner, we could have written instead:" + ] + }, { "cell_type": "code", "execution_count": 3, - "id": "e4715f11-f4c3-4953-8d6e-7d01cd14d32a", + "id": "220429e4-fd82-4495-873e-a08ae8c09a9a", "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "2" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "vsys.order" + "vsys = hy.var_ode_sys(sys, [x, v], order=2)" + ] + }, + { + "cell_type": "markdown", + "id": "6c86a6e8-20d2-4d78-8846-46ed6628e823", + "metadata": {}, + "source": [ + "In this case, instead of a {class}`~heyoka.var_args` enumerator, we passed an explicit list of state variables with respect to whose initial conditions we want to formulate the variational equations. Please refer to the documentation of {class}`~heyoka.var_ode_sys` for an exhaustive explanation of what can be passed as second argument to the constructor.\n", + "\n", + "Let us explore a bit the {class}`~heyoka.var_ode_sys` class. First of all, we can access the variational system of equations:" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 4, "id": "16e3bc0a-e7d5-499c-88bd-e17aaa6c1657", "metadata": {}, "outputs": [ @@ -143,7 +151,7 @@ " (-(∂[(1, 2)]v * p0) - (((∂[(1, 1)]x * -sin(x)) * ∂[(1, 1)]x) + (∂[(1, 2)]x * cos(x)))))]" ] }, - "execution_count": 12, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } @@ -152,22 +160,605 @@ "vsys.sys" ] }, + { + "cell_type": "markdown", + "id": "c969185b-7df4-4f95-9ab1-04af1205abdd", + "metadata": {}, + "source": [ + "The first two equations are from the original system of ODEs, while the remaining ones are the variational equations. The names of the variational variables begin with the $\\partial$ symbol, followed by a sparse multiindex encoding of the differentiation indices. For instance, the variational variable ``∂[(0, 1)]x`` is the first-order derivative of $x$ with respect to the first variational argument, that is, $\\frac{\\partial x}{\\partial x_0}$. Similarly, ``∂[(0, 1), (1, 1)]x`` is the second order derivative of $x$ with respect to both variational arguments, that is, $\\frac{\\partial^2 x}{\\partial x_0 \\partial y_0}$. The ordering of the variational equations follows the same scheme explained in the tutorial about [computing derivatives](<./computing_derivatives.ipynb>).\n", + "\n", + "We can also query other properties of ``vsys``. For instance, the differentiation order:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "e4715f11-f4c3-4953-8d6e-7d01cd14d32a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vsys.order" + ] + }, + { + "cell_type": "markdown", + "id": "27e8ac8e-1d81-4309-84d7-f18f16fc2d06", + "metadata": {}, + "source": [ + "The number of state variables in the original ODE system:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "2fad08b8-23d0-4572-a74c-edc42ff93c3b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vsys.n_orig_sv" + ] + }, + { + "cell_type": "markdown", + "id": "f5d76bce-6e58-4e87-b5ba-98fdb953dcc1", + "metadata": {}, + "source": [ + "And the list of arguments with respect to which the variational equations are formulated:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "9369b6f1-9b2c-4e7e-b70e-21679dd62fda", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[x, v]" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vsys.vargs" + ] + }, + { + "cell_type": "markdown", + "id": "276d5cc7-615c-44fe-b08d-48a9cee19fcb", + "metadata": {}, + "source": [ + "## Constructing a variational integrator\n", + "\n", + "After the construction of a variational ODE system, we are now ready to construct a variational integrator. We can do this by simply passing the variational ODE system (instead of the original, non-variational ODE system) as first input argument to the constructor:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "8c51dcf8-a089-42c2-9b84-5e2d5c22da24", + "metadata": {}, + "outputs": [], + "source": [ + "# Construct a variational integrator.\n", + "ta_var = hy.taylor_adaptive(vsys, [.2, .3], pars=[.4], compact_mode=True)" + ] + }, + { + "cell_type": "markdown", + "id": "5925b0ce-1b04-4174-86ef-593c7ffade74", + "metadata": {}, + "source": [ + "Note how we constructed the integrator with {ref}`compact mode ` enabled: the formulation of the variational equations, especially at high differentiation orders, greatly increases the size of the symbolic expressions that need to be just-in-time compiled during the creation of the integrator. By enabling compact mode, we keep the compilation time at manageable levels.\n", + "\n", + "Let us inspect the integrator:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "844c0ae8-9f32-4d09-a900-8098301863ff", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "C++ datatype : double\n", + "Tolerance : 2.220446049250313e-16\n", + "High accuracy : false\n", + "Compact mode : true\n", + "Taylor order : 20\n", + "Dimension : 12\n", + "Time : 0\n", + "State : [0.2, 0.3, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0]\n", + "Parameters : [0.4]\n", + "Variational order : 2" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta_var" + ] + }, + { + "cell_type": "markdown", + "id": "afc46400-dc3a-4e1d-bae6-aed5771d46a6", + "metadata": {}, + "source": [ + "The screen output informs us that ``ta_var`` is a variational integrator of order 2. We can also see that, although on construction we passed the initial conditions only for the ``x`` and ``v`` state variables, the integrator automatically set up appropriate initial conditions for the variational variables. Indeed, with respect to a non-variational integrator, the state vector has been augmented to store also the values of the variational variables:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "81cbba40-07fb-45b1-9622-29e9c966b273", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.2, 0.3, 1. , 0. , 0. , 1. , 0. , 0. , 0. , 0. , 0. , 0. ])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta_var.state" + ] + }, + { + "cell_type": "markdown", + "id": "8e5b8bd0-bf91-41f5-baf6-c1cb2d2e2e7a", + "metadata": {}, + "source": [ + "Alternatively, instead of relying on the integrator to automatically set up the initial conditions of the variational variables, we could also pass a full 12-elements vector of initial conditions in input.\n", + "\n", + "We are now ready to proceed to our first variational integration:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "afa9f717-bc57-4f19-9a25-198bff3cdd14", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.11757215, -1.24940656, -0.425478 , 0.41983649, -0.19171818,\n", + " -0.51871994, 0.27771857, 0.22392433, 0.60414865, -0.1122785 ,\n", + " -0.12689065, 0.12443791])" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Propagate until t=3.\n", + "ta_var.propagate_until(3.)\n", + "\n", + "# Print the full state vector.\n", + "ta_var.state" + ] + }, + { + "cell_type": "markdown", + "id": "a51db4d7-5c9f-446f-ae1f-eb9fc759c94d", + "metadata": {}, + "source": [ + "That is all fine an good, but how do we fetch the values of the derivatives we are interested in at the end of an integration? As mentioned earlier, the partial derivatives are ordered in the state vector following the same criterion explained in the tutorial about [computing derivatives](<./computing_derivatives.ipynb>): first by total order of differentiation, then by component (i.e., the derivatives of ``x`` precede the derivatives of ``v``) and finally by reverse lexicographic order with respect to the differentiation multiindices. However, navigating by hand this ordering scheme can be complicated, especially at high differentiation orders.\n", + "\n", + "Variational integrators provide a couple of methods that facilitate the task of locating specific derivatives in the state vector. The first helper is ``get_vslice()``. This method takes in input a differentiation order and, optionally, a component index, and returns a {class}`slice` into the state vector corresponding to the range of indices for the requested derivatives. Let us see a couple of examples:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "a78531a3-a1f9-4a60-bc21-5c37ab13f6ad", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "slice(6, 12, None)" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Fetch the range of all order-2 derivatives.\n", + "sl = ta_var.get_vslice(order=2)\n", + "sl" + ] + }, + { + "cell_type": "markdown", + "id": "bfd3b6b5-8cd6-4c92-be86-d1125cee17bf", + "metadata": {}, + "source": [ + "That is, the order-2 derivatives are between indices 6 and 12 in the state vector. We can use ``sl`` to index directly into the state vector:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "3fed6708-3e5e-44c2-871a-0237d005937e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.27771857, 0.22392433, 0.60414865, -0.1122785 , -0.12689065,\n", + " 0.12443791])" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Fetch all order-2 derivatives from the state vector.\n", + "ta_var.state[sl]" + ] + }, + { + "cell_type": "markdown", + "id": "455cbfba-964d-4cbe-b6c4-efdde7f85ad2", + "metadata": {}, + "source": [ + "If we are interested only in the order-2 derivatives of ``v``, we can pass the additional ``component`` keyword argument:" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "ff3ae54f-455d-44d3-b166-a0a607434d6a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([-0.1122785 , -0.12689065, 0.12443791])" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Fetch the range of the order-2 derivatives for v.\n", + "# component=1 means the second original state variable,\n", + "# i.e., v (component=0 would fetch the derivatives for x).\n", + "sl = ta_var.get_vslice(order=2, component=1)\n", + "\n", + "# Fetch the order-2 derivatives for v.\n", + "ta_var.state[sl]" + ] + }, + { + "cell_type": "markdown", + "id": "6fa1d1cd-840e-491c-aefb-535518c7bd2e", + "metadata": {}, + "source": [ + "Often fetching the values of the derivatives is not enough, and we also need to access the differentiation multiindices associated to each derivative. In order to do this, we can use the ``get_mindex()`` method, which takes in input a single index into the state vector and returns the corresponding differentiation multiindex.\n", + "\n", + "Let us see a couple of examples:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "5ae39909-0554-4def-b04b-4ed69ade36b7", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[0, 0, 0]" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta_var.get_mindex(i=0)" + ] + }, + { + "cell_type": "markdown", + "id": "1018d0c7-e7f3-463a-8e88-618d1600be31", + "metadata": {}, + "source": [ + "At ``i=0`` in the state vector we have the order-0 derivative of the first state variable - that is, a complicated way of saying that we have the current value of ``x``. Let us see a more useful example:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "064adc95-4ed0-4c6e-b927-57f8a607483c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Multiindex: [0, 2, 0], derivative value: 0.27771856625922825\n", + "Multiindex: [0, 1, 1], derivative value: 0.2239243290399873\n", + "Multiindex: [0, 0, 2], derivative value: 0.6041486469599123\n", + "Multiindex: [1, 2, 0], derivative value: -0.1122785039896587\n", + "Multiindex: [1, 1, 1], derivative value: -0.12689065229074553\n", + "Multiindex: [1, 0, 2], derivative value: 0.12443790781972988\n" + ] + } + ], + "source": [ + "# Fetch the range of all order-2 derivatives.\n", + "sl = ta_var.get_vslice(order=2)\n", + "\n", + "# Print the multiindices and associated values.\n", + "for idx in range(sl.start, sl.stop):\n", + " print(f\"Multiindex: {ta_var.get_mindex(idx)}, derivative value: {ta_var.state[idx]}\")" + ] + }, { "cell_type": "markdown", "id": "d2ae2f2e-3b67-47be-a863-b17f44fb26c2", "metadata": {}, "source": [ + "Recall that in a multiindex the first number refers to the component index (i.e., 0 for ``x`` and 1 for ``v``), while the remaining indices refer to the differentiation orders with respect to the variational arguments.\n", + "\n", "(taylor_map)=\n", - "## Taylor map evaluation" + "## Taylor map evaluation\n", + "\n", + "One of the most useful applications of the variational equations is the ability to compute how a small perturbation on the initial conditions and/or parameters of the system affects the current state of the system, and to do it quickly (i.e., without having to repeat the numerical integration with the updated initial conditions/parameters). This is accomplished by using the values of the partial derivatives to construct and evaluate the multivariate Taylor series of the solution around the original initial conditions/parameters of the system. This approach, when applied to perturbations on the initial conditions, is sometimes called *jet transport* in the literature. Here, more generally, we will call it *evaluation of the Taylor map*.\n", + "\n", + "Variational integrators provide a specific method called ``eval_taylor_map()`` to construct and evaluate the Taylor map. Let us see a simple example. We begin by re-creating from scratch our variational integrator:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "229a9a4e-bf5a-4b57-a599-e4571ebd7507", + "metadata": {}, + "outputs": [], + "source": [ + "ta_var = hy.taylor_adaptive(vsys, [.2, .3], pars=[.4], compact_mode=True)" + ] + }, + { + "cell_type": "markdown", + "id": "b09de854-aa19-4f08-a36d-66b926bd1d6b", + "metadata": {}, + "source": [ + "We define two small displacements on the state variables ``x`` and ``v``:" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "7c2d30d2-17cd-4ab4-889a-8bcb72bb9b11", + "metadata": {}, + "outputs": [], + "source": [ + "dx = 1e-4\n", + "dv = -2e-4" + ] + }, + { + "cell_type": "markdown", + "id": "eb9e6427-fb2c-4f11-91b2-85d9d45577b7", + "metadata": {}, + "source": [ + "And we create a non-variational integrator with displaced initial conditions with respect to the variational one:" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "1ab3b27e-3506-47e0-a0bb-35cd2076238f", + "metadata": {}, + "outputs": [], + "source": [ + "# Non-variational integrator with displaced\n", + "# initial conditions.\n", + "ta = hy.taylor_adaptive(sys, [.2 + dx, .3 + dv], pars=[.4], compact_mode=True)" + ] + }, + { + "cell_type": "markdown", + "id": "8fa77407-e47a-4338-b65e-f314a4a7afdd", + "metadata": {}, + "source": [ + "Next, we propagate both integrators up to $t=3$:" ] }, { "cell_type": "code", - "execution_count": null, - "id": "40d0734e-25b3-4bc5-9f2b-13007790e3a0", + "execution_count": 20, + "id": "42a8ca75-4eca-4da2-b33e-838085fd497f", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "ta_var.propagate_until(3.)\n", + "ta.propagate_until(3.);" + ] + }, + { + "cell_type": "markdown", + "id": "80f7bd5e-606b-425c-bfe8-e3b409a1975b", + "metadata": {}, + "source": [ + "Clearly, the values of ``x`` and ``v`` will differ in the two integrators due to the different initial conditions:" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "fa2c2e21-81f3-4e2a-a219-296b940bb3b9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Non-variational state: [ 0.11744564 -1.24932198]\n", + "Variational state : [ 0.11757215 -1.24940656]\n" + ] + } + ], + "source": [ + "print(f\"Non-variational state: {ta.state}\")\n", + "print(f\"Variational state : {ta_var.state[:2]}\")" + ] + }, + { + "cell_type": "markdown", + "id": "b102f4fd-63dc-45a6-81c1-1cf925aa8b15", + "metadata": {}, + "source": [ + "We can now use the ``eval_taylor_map()`` method on the variational integrator to compute the effect of the displacements ``dx`` and ``dv`` on the state of the system:" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "b7f778cc-db85-43aa-bf29-5445d6170958", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.11744564, -1.24932198])" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta_var.eval_taylor_map([dx, dv])" + ] + }, + { + "cell_type": "markdown", + "id": "a3f8b0d2-0446-4dd5-a8fb-0539863ed161", + "metadata": {}, + "source": [ + "``eval_taylor_map()`` takes in input a vector of displacements (one for each variational argument), and computes and evaluates the Taylor map, returning a reference to an internal array in ``ta_var`` storing the result of the evaluation (i.e., the updated values of the state variables). We can see that, in this specific case, the evaluation of the Taylor map reproduces accurately the state vector computed by the non-variational integrator with displaced initial conditions:" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "d561e6f0-2500-4faf-9caa-f8bf29ff99e4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Taylor map error: [6.58681443e-13 6.68798350e-13]\n" + ] + } + ], + "source": [ + "print(f'Taylor map error: {ta_var.eval_taylor_map([dx, dv]) - ta.state}')" + ] + }, + { + "cell_type": "markdown", + "id": "fb423c90-3458-4939-82ed-97326baf8be8", + "metadata": {}, + "source": [ + "Note that the Taylor map state vector can also be fetched via the ``tstate`` property of the integrator:" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "d1dc660b-aa56-47d0-bd83-d2f323b12dfd", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.11744564, -1.24932198])" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta_var.tstate" + ] + }, + { + "cell_type": "markdown", + "id": "0152b819-815d-4082-9e26-4da55c731d57", + "metadata": {}, + "source": [ + "```{warning}\n", + "Accessing the Taylor map state vector via ``tstate`` will **NOT** trigger any Taylor map evaluation, it will merely return a reference to the internal array storing the result of the evaluation. It is your responsibility to ensure that you called ``eval_taylor_map()`` before accessing this array via ``tstate``.\n", + "\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "fc4df22b-1edf-4700-987b-0b663999fb59", + "metadata": {}, + "source": [ + "## A note on computational efficiency\n", + "\n", + "{class}`~heyoka.var_ode_sys` uses internally the {func}`~heyoka.diff_tensors()` and {class}`~heyoka.dtens` API to formulate the variational equations. This means that the computation of the symbolic derivatives is performed in an efficient manner. For instance, reverse-mode symbolic automatic differentiation will be employed when computing the first-order variationals of ODE systems containing a large number of parameters (e.g., in [neural ODEs](./NeuralODEs.ipynb)).\n", + "\n", + "See the [computing derivatives](<./computing_derivatives.ipynb>) tutorial for a more in-depth discussion of how heyoka.py computes symbolic derivatives." + ] } ], "metadata": { From 3ded29752934be918ac0c0854110afc5a96b0d7d Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Tue, 11 Jun 2024 10:56:26 +0200 Subject: [PATCH 16/23] sympy fixes for time changes. --- heyoka/__init__.py | 1 + heyoka/_sympy_utils.py | 4 ++-- heyoka/_test_sympy.py | 15 ++++++++++++--- 3 files changed, 15 insertions(+), 5 deletions(-) diff --git a/heyoka/__init__.py b/heyoka/__init__.py index 72da6d02..a358d79e 100644 --- a/heyoka/__init__.py +++ b/heyoka/__init__.py @@ -383,6 +383,7 @@ def _create_par(): """ + # Machinery for the time attribute. def _create_time(): from . import core diff --git a/heyoka/_sympy_utils.py b/heyoka/_sympy_utils.py index baf33f6a..c3a20903 100644 --- a/heyoka/_sympy_utils.py +++ b/heyoka/_sympy_utils.py @@ -142,7 +142,7 @@ def _build_fmap(): if not _with_sympy: return None - from . import core, pi + from . import core, pi, time as htime retval = {} @@ -179,7 +179,7 @@ def mul_wrapper(*args): retval[_spy.Function("heyoka_kepE")] = core.kepE retval[_spy.Function("heyoka_kepF")] = core.kepF retval[_spy.Function("heyoka_kepDE")] = core.kepDE - retval[_spy.Function("heyoka_time")] = lambda: core.time + retval[_spy.Function("heyoka_time")] = lambda: htime return retval diff --git a/heyoka/_test_sympy.py b/heyoka/_test_sympy.py index 3c3d5863..fc7ae17b 100644 --- a/heyoka/_test_sympy.py +++ b/heyoka/_test_sympy.py @@ -192,7 +192,16 @@ def test_func_conversion(self): import sympy as spy - from . import core, make_vars, from_sympy, to_sympy, pi, sum as hsum, prod + from . import ( + core, + make_vars, + from_sympy, + to_sympy, + pi, + sum as hsum, + prod, + time as htime, + ) from .model import nbody @@ -320,8 +329,8 @@ def test_func_conversion(self): self.assertEqual(to_sympy(core.sigmoid(hx + hy)), 1.0 / (1.0 + spy.exp(-x - y))) - self.assertEqual(core.time, from_sympy(spy.Function("heyoka_time")())) - self.assertEqual(to_sympy(core.time), spy.Function("heyoka_time")()) + self.assertEqual(htime, from_sympy(spy.Function("heyoka_time")())) + self.assertEqual(to_sympy(htime), spy.Function("heyoka_time")()) with self.assertRaises(TypeError) as cm: from_sympy(abs(x)) From 44114aaaad00bb02591ec140570a1236b16e8f46 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Tue, 11 Jun 2024 17:41:25 +0200 Subject: [PATCH 17/23] Tutorial tweaks. --- doc/notebooks/var_ode_sys.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/notebooks/var_ode_sys.ipynb b/doc/notebooks/var_ode_sys.ipynb index 697a9dbd..347a92c3 100644 --- a/doc/notebooks/var_ode_sys.ipynb +++ b/doc/notebooks/var_ode_sys.ipynb @@ -27,7 +27,7 @@ "\n", "When solving numerically initial-value problems, it is often useful to compute not only the solution, but also its partial derivatives with respect to the initial conditions and/or the parameters. The derivatives with respect to the initial conditions, for instance, are needed for the computation of [chaos indicators](https://en.wikipedia.org/wiki/Lyapunov_exponent) and for [uncertainty propagation](https://en.wikipedia.org/wiki/Propagation_of_uncertainty), and they can also be used to propagate a small neighborhood in phase space around the initial conditions. The derivatives with respect to the parameters of the system are required when formulating optimisation and inversion problems such as orbit determination, trajectory optimisation and training of neural networks in [neural ODEs](./NeuralODEs.ipynb).\n", "\n", - "There are two main methods for the computation of the partial derivatives. The first one is based on the application of automatic differentiation (AD) techniques directly to the numerical integration algorithm. This can be done either by replacing the algebra of floating-point numbers with the algebra of (generalised) [dual numbers](https://en.wikipedia.org/wiki/Dual_number) (aka truncated Taylor polynomials), or via [differentiable programming](https://en.wikipedia.org/wiki/Differentiable_programming) techniques. The former approach is used by libraries such as [pyaudi](https://github.com/darioizzo/audi), [desolver](https://github.com/Microno95/desolver) and [TaylorIntegration.jl](https://docs.sciml.ai/TaylorIntegration/stable/jet_transport/), while differentiable programming is popular in the [Julia programming language](https://en.wikipedia.org/wiki/Julia_(programming_language)) community.\n", + "There are two main methods for the computation of the partial derivatives. The first one is based on the application of automatic differentiation (AD) techniques directly to the numerical integration algorithm. This can be done either by replacing the algebra of floating-point numbers with the algebra of (generalised) [dual numbers](https://en.wikipedia.org/wiki/Dual_number) (aka truncated Taylor polynomials), or via [differentiable programming](https://en.wikipedia.org/wiki/Differentiable_programming) techniques. The former approach is used by libraries such as [pyaudi](https://github.com/darioizzo/audi), [desolver](https://github.com/Microno95/desolver) and [TaylorIntegration.jl](https://docs.sciml.ai/TaylorIntegration/stable/jet_transport/), while differentiable programming is popular in the machine learning community with projects such as [PyTorch](https://pytorch.org/), [JAX](https://jax.readthedocs.io/en/latest/) and [TensorFlow](https://www.tensorflow.org/). Differentiable programming is also popular in the [Julia programming language](https://en.wikipedia.org/wiki/Julia_(programming_language)) community.\n", "\n", "The second method is based on the formulation of the *variational equations*, that is, differential equations satisfied by the partial derivatives which are added to and solved together with the original ODEs. For instance, we can formulate differential equations for the first-order derivatives with respect to the initial conditions via elementary calculus:\n", "\n", @@ -381,9 +381,9 @@ "id": "a51db4d7-5c9f-446f-ae1f-eb9fc759c94d", "metadata": {}, "source": [ - "That is all fine an good, but how do we fetch the values of the derivatives we are interested in at the end of an integration? As mentioned earlier, the partial derivatives are ordered in the state vector following the same criterion explained in the tutorial about [computing derivatives](<./computing_derivatives.ipynb>): first by total order of differentiation, then by component (i.e., the derivatives of ``x`` precede the derivatives of ``v``) and finally by reverse lexicographic order with respect to the differentiation multiindices. However, navigating by hand this ordering scheme can be complicated, especially at high differentiation orders.\n", + "That is all fine and good, but how do we fetch the values of the derivatives we are interested in at the end of an integration? As mentioned earlier, the partial derivatives are ordered in the state vector following the same criterion explained in the tutorial about [computing derivatives](<./computing_derivatives.ipynb>): first by total order of differentiation, then by component (i.e., the derivatives of ``x`` precede the derivatives of ``v``) and finally by reverse lexicographic order with respect to the differentiation multiindices. However, navigating by hand this ordering scheme can be complicated, especially at high differentiation orders.\n", "\n", - "Variational integrators provide a couple of methods that facilitate the task of locating specific derivatives in the state vector. The first helper is ``get_vslice()``. This method takes in input a differentiation order and, optionally, a component index, and returns a {class}`slice` into the state vector corresponding to the range of indices for the requested derivatives. Let us see a couple of examples:" + "Variational integrators provide a couple of methods that facilitate the task of locating specific derivatives in the state vector. The first helper is ``get_vslice()``. This method takes as input a differentiation order and, optionally, a component index, and returns a {class}`slice` into the state vector corresponding to the range of indices for the requested derivatives. Let us see a couple of examples:" ] }, { From 273f02bbf550e5d6295bdfc8c276d38f261c79f8 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Tue, 11 Jun 2024 17:47:01 +0200 Subject: [PATCH 18/23] Update changelog, readme. --- README.md | 2 ++ doc/changelog.rst | 5 +++++ doc/index.md | 2 ++ 3 files changed, 9 insertions(+) diff --git a/README.md b/README.md index 050fbfb1..c453c727 100644 --- a/README.md +++ b/README.md @@ -43,6 +43,8 @@ compilation via [LLVM](https://llvm.org/). Notable features include: * accurate and reliable event detection, * builtin support for analytical mechanics - bring your own Lagrangians/Hamiltonians and let heyoka.py formulate and solve the equations of motion, +* builtin support for high-order variational equations - compute not only the solution, + but also its partial derivatives, * builtin support for machine learning applications via neural network models, * the ability to maintain machine precision accuracy over tens of billions of timesteps, diff --git a/doc/changelog.rst b/doc/changelog.rst index b1ebf5b1..a86a84f0 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -9,6 +9,11 @@ Changelog New ~~~ +- Add support for variational ODE systems and Taylor map computation + (`#177 `__). +- Add thermonets: neural, differentiable, high-performance + models for the Earth's thermosphere density + (`#176 `__). - Add a vectorised implementation of ``diff()`` (`#173 `__). diff --git a/doc/index.md b/doc/index.md index 0f3261fd..acb9c2dc 100644 --- a/doc/index.md +++ b/doc/index.md @@ -16,6 +16,8 @@ compilation via [LLVM](https://llvm.org/). Notable features include: - accurate and reliable event detection, - builtin support for analytical mechanics - bring your own Lagrangians/Hamiltonians and let heyoka.py formulate and solve the equations of motion, +- builtin support for high-order variational equations - compute not only the solution, + but also its partial derivatives, - builtin support for machine learning applications via neural network models, - the ability to maintain machine precision accuracy over tens of billions of timesteps, From a87b0bebd3a7bf32964a95720cd629d80a7644cc Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Tue, 11 Jun 2024 17:50:46 +0200 Subject: [PATCH 19/23] Remove the now-outdated tutorial on variational equations. --- doc/examples_others.rst | 2 - doc/notebooks/The variational equations.ipynb | 300 ------------------ 2 files changed, 302 deletions(-) delete mode 100644 doc/notebooks/The variational equations.ipynb diff --git a/doc/examples_others.rst b/doc/examples_others.rst index d6b6d0dd..d12f4292 100644 --- a/doc/examples_others.rst +++ b/doc/examples_others.rst @@ -27,7 +27,5 @@ Others :maxdepth: 1 notebooks/ensemble_batch_perf.ipynb - notebooks/The variational equations notebooks/Optimal Control of the Lotka-Volterra equations notebooks/definite_integrals - diff --git a/doc/notebooks/The variational equations.ipynb b/doc/notebooks/The variational equations.ipynb deleted file mode 100644 index 2bcf22fe..00000000 --- a/doc/notebooks/The variational equations.ipynb +++ /dev/null @@ -1,300 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "excited-uganda", - "metadata": {}, - "source": [ - "The variational equations\n", - "=========================\n", - "\n", - "In this tutorial, we will show how it is possible to exploit heyoka.py's [expression system](<./The expression system.ipynb>) to implement first-order variational equations for a system of ODEs.\n", - "\n", - "\n", - "First, let's start with a brief recap on the variational equations. For simplicity, we will consider an autonomous ODE system in the variables $x$ and $y$ (the extension to more variables just complicates the notation):\n", - "\n", - "$$\n", - "\\begin{cases}\n", - "\\frac{dx}{dt} = f\\left(x, y \\right)\\\\\n", - "\\frac{dy}{dt} = g\\left(x, y \\right)\n", - "\\end{cases},\n", - "$$\n", - "\n", - "with initial conditions\n", - "\n", - "$$\n", - "\\begin{cases}\n", - "x\\left( 0 \\right) = x_0\\\\\n", - "y\\left( 0 \\right) = y_0\n", - "\\end{cases}.\n", - "$$\n", - "\n", - "The solutions of this ODE system will be the functions $x\\left( t; x_0, y_0 \\right)$ and $y\\left( t; x_0, y_0 \\right)$. Our objective is to determine the evolution in time of the first-order derivatives of these solutions with respect to the initial conditions $\\left( x_0, y_0 \\right)$:\n", - "\n", - "$$\n", - "\\begin{aligned}\n", - "\\left(\\frac{\\partial x}{\\partial x_0}\\right)\\left(t; x_0, y_0\\right),\\\\\n", - "\\left(\\frac{\\partial x}{\\partial y_0}\\right)\\left(t; x_0, y_0\\right),\\\\\n", - "\\left(\\frac{\\partial y}{\\partial x_0}\\right)\\left(t; x_0, y_0\\right),\\\\\n", - "\\left(\\frac{\\partial y}{\\partial y_0}\\right)\\left(t; x_0, y_0\\right).\\\\\n", - "\\end{aligned}\n", - "$$\n", - "\n", - "These functions tell us how the the solutions of the ODE system react to changes in the initial values $\\left( x_0, y_0 \\right)$. One way of determining numerically the evolution in time of the partial derivatives is to introduce them as additional state variables in the original ODE system, together with an additional set of differential equations (i.e., the variational equations). Let's see how this is done for, e.g., $\\frac{\\partial x}{\\partial x_0}$. We take the time derivative of $\\frac{\\partial x}{\\partial x_0}$,\n", - "\n", - "$$\n", - "\\frac{d \\left( \\frac{\\partial x}{\\partial x_0} \\right)}{dt},\n", - "$$\n", - "\n", - "and expand it, via elementary calculus rules, to\n", - "\n", - "$$\n", - "\\frac{d \\left( \\frac{\\partial x}{\\partial x_0} \\right)}{dt} = \\frac{\\partial \\left( \\frac{dx}{dt} \\right) }{\\partial x_0}=\\frac{\\partial f\\left(x, y\\right)}{\\partial x}\\frac{\\partial x}{\\partial x_0}+\\frac{\\partial f \\left(x, y\\right)}{\\partial y}\\frac{\\partial y}{\\partial x_0}.\n", - "$$\n", - "\n", - "We can write similar ODEs for the other partial derivatives, append them below the original ODE system and finally obtain an augmented ODE system containing the variational equations:\n", - "\n", - "$$\n", - "\\begin{cases}\n", - "\\frac{dx}{dt} = f\\left(x, y \\right)\\\\\n", - "\\frac{dy}{dt} = g\\left(x, y \\right)\\\\\n", - "\\frac{d \\left( \\frac{\\partial x}{\\partial x_0} \\right)}{dt} = \\frac{\\partial f\\left(x, y\\right)}{\\partial x}\\frac{\\partial x}{\\partial x_0}+\\frac{\\partial f\\left(x, y\\right)}{\\partial y}\\frac{\\partial y}{\\partial x_0}\\\\\n", - "\\frac{d \\left( \\frac{\\partial x}{\\partial y_0} \\right)}{dt} = \\frac{\\partial f\\left(x, y\\right)}{\\partial x}\\frac{\\partial x}{\\partial y_0}+\\frac{\\partial f\\left(x, y\\right)}{\\partial y}\\frac{\\partial y}{\\partial y_0}\\\\\n", - "\\frac{d \\left( \\frac{\\partial y}{\\partial x_0} \\right)}{dt} = \\frac{\\partial g\\left(x, y\\right)}{\\partial x}\\frac{\\partial x}{\\partial x_0}+\\frac{\\partial g\\left(x, y\\right)}{\\partial y}\\frac{\\partial y}{\\partial x_0}\\\\\n", - "\\frac{d \\left( \\frac{\\partial y}{\\partial y_0} \\right)}{dt} = \\frac{\\partial g\\left(x, y\\right)}{\\partial x}\\frac{\\partial x}{\\partial y_0}+\\frac{\\partial g\\left(x, y\\right)}{\\partial y}\\frac{\\partial y}{\\partial y_0}\n", - "\\end{cases}.\n", - "$$\n", - "\n", - "The introduction of the variational equations in an ODE system of $n$ equations results in additional $n^2$ equations and variables. The initial conditions for the new variables are:\n", - "\n", - "$$\n", - "\\begin{cases}\n", - "\\frac{\\partial x}{\\partial x_0} \\left( 0 \\right) = 1\\\\\n", - "\\frac{\\partial x}{\\partial y_0} \\left( 0 \\right) = 0\\\\\n", - "\\frac{\\partial y}{\\partial x_0} \\left( 0 \\right) = 0\\\\\n", - "\\frac{\\partial y}{\\partial y_0} \\left( 0 \\right) = 1\\\\\n", - "\\end{cases}.\n", - "$$\n", - "\n", - "Let us now try to implement the variational equations for the simple pendulum:" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "id": "polyphonic-argentina", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "Tolerance : 2.2204460492503131e-16\n", - "High accuracy : false\n", - "Compact mode : false\n", - "Taylor order : 20\n", - "Dimension : 6\n", - "Time : 0.0000000000000000\n", - "State : [0.050000000000000003, 0.025000000000000001, 1.0000000000000000, 0.0000000000000000, 0.0000000000000000, 1.0000000000000000]" - ] - }, - "execution_count": 1, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "import heyoka as hy\n", - "\n", - "# Create the state variables x and v.\n", - "x, v = hy.make_vars(\"x\", \"v\")\n", - "\n", - "# Create the symbolic variables for the variational equations.\n", - "x_x0, x_y0, y_x0, y_y0 = hy.make_vars(\"x_x0\", \"x_y0\", \"y_x0\", \"y_y0\")\n", - "\n", - "# Create the right-hand side of the ODE system.\n", - "f = v\n", - "g = -9.8 * hy.sin(x)\n", - "x_x0_p = hy.diff(f, x) * x_x0 + hy.diff(f, v) * y_x0\n", - "x_y0_p = hy.diff(f, x) * x_y0 + hy.diff(f, v) * y_y0\n", - "y_x0_p = hy.diff(g, x) * x_x0 + hy.diff(g, v) * y_x0\n", - "y_y0_p = hy.diff(g, x) * x_y0 + hy.diff(g, v) * y_y0\n", - "\n", - "# Create the integrator object.\n", - "ta = hy.taylor_adaptive(\n", - " # Definition of the ODE system.\n", - " sys = [\n", - " (x, f),\n", - " (v, g),\n", - " (x_x0, x_x0_p),\n", - " (x_y0, x_y0_p),\n", - " (y_x0, y_x0_p),\n", - " (y_y0, y_y0_p),\n", - " ],\n", - " # Initial conditions.\n", - " state = [0.05, 0.025, 1., 0., 0., 1.])\n", - "\n", - "ta" - ] - }, - { - "cell_type": "markdown", - "id": "unlike-fifty", - "metadata": {}, - "source": [ - "Here we used the ``diff()`` function to perform the symbolic differentiation of ``f`` and ``g`` with respect to the state variables.\n", - "\n", - "Let us now propagate the state of the system (including the variational equations) up to $t=50$:" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "recorded-shelf", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([ 0.03744787, 0.10667026, 0.80315068, -0.17679098, 1.82916222,\n", - " 0.84245788])" - ] - }, - "execution_count": 2, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "ta.propagate_until(50.)\n", - "ta.state" - ] - }, - { - "cell_type": "markdown", - "id": "blessed-arabic", - "metadata": {}, - "source": [ - "Next, we are going to create another integrator for the simple pendulum, this time *without* the variational equations and perturbing the initial value of $x$ by $10^{-8}$ with respect to the original initial conditions:" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "compound-installation", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "(,\n", - " 0.20208180443700685,\n", - " 0.2182277296610072,\n", - " 240,\n", - " None)" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "ta_dx = hy.taylor_adaptive(\n", - " # Definition of the ODE system.\n", - " sys = [\n", - " (x, f),\n", - " (v, g)\n", - " ],\n", - " # Initial conditions.\n", - " state = [0.05 + 1e-8, 0.025])\n", - "\n", - "ta_dx.propagate_until(50.)" - ] - }, - { - "cell_type": "markdown", - "id": "piano-thinking", - "metadata": {}, - "source": [ - "Because of the perturbation on the initial value of $x$, the value of $x$ at $t=50$ will differ in the two integrators:" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "express-pittsburgh", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "-8.031506786021492e-09" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "ta.state[0] - ta_dx.state[0]" - ] - }, - { - "cell_type": "markdown", - "id": "relative-briefing", - "metadata": {}, - "source": [ - "However, because in the variational integrator we now have the value of $\\frac{\\partial x}{\\partial x_0}$ at $t=50$, we can compute an approximation of the $\\Delta x$ induced by the perturbation on the initial state as\n", - "\n", - "$$\n", - "\\Delta x \\approx \\frac{\\partial x}{\\partial x_0}\\left( 50 \\right) \\cdot 10^{-8}.\n", - "$$\n", - "\n", - "Indeed:" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "confused-excellence", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "2.0816681711721685e-17" - ] - }, - "execution_count": 5, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "(ta.state[0] + ta.state[2]*1e-8) - ta_dx.state[0]" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.13" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} From 5300e720e3e87a5bc9b3266f42fb67ed93f12287 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Wed, 12 Jun 2024 08:13:31 +0200 Subject: [PATCH 20/23] Revert "Remove the now-outdated tutorial on variational equations." This reverts commit a87b0bebd3a7bf32964a95720cd629d80a7644cc. --- doc/examples_others.rst | 2 + doc/notebooks/The variational equations.ipynb | 300 ++++++++++++++++++ 2 files changed, 302 insertions(+) create mode 100644 doc/notebooks/The variational equations.ipynb diff --git a/doc/examples_others.rst b/doc/examples_others.rst index d12f4292..d6b6d0dd 100644 --- a/doc/examples_others.rst +++ b/doc/examples_others.rst @@ -27,5 +27,7 @@ Others :maxdepth: 1 notebooks/ensemble_batch_perf.ipynb + notebooks/The variational equations notebooks/Optimal Control of the Lotka-Volterra equations notebooks/definite_integrals + diff --git a/doc/notebooks/The variational equations.ipynb b/doc/notebooks/The variational equations.ipynb new file mode 100644 index 00000000..2bcf22fe --- /dev/null +++ b/doc/notebooks/The variational equations.ipynb @@ -0,0 +1,300 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "excited-uganda", + "metadata": {}, + "source": [ + "The variational equations\n", + "=========================\n", + "\n", + "In this tutorial, we will show how it is possible to exploit heyoka.py's [expression system](<./The expression system.ipynb>) to implement first-order variational equations for a system of ODEs.\n", + "\n", + "\n", + "First, let's start with a brief recap on the variational equations. For simplicity, we will consider an autonomous ODE system in the variables $x$ and $y$ (the extension to more variables just complicates the notation):\n", + "\n", + "$$\n", + "\\begin{cases}\n", + "\\frac{dx}{dt} = f\\left(x, y \\right)\\\\\n", + "\\frac{dy}{dt} = g\\left(x, y \\right)\n", + "\\end{cases},\n", + "$$\n", + "\n", + "with initial conditions\n", + "\n", + "$$\n", + "\\begin{cases}\n", + "x\\left( 0 \\right) = x_0\\\\\n", + "y\\left( 0 \\right) = y_0\n", + "\\end{cases}.\n", + "$$\n", + "\n", + "The solutions of this ODE system will be the functions $x\\left( t; x_0, y_0 \\right)$ and $y\\left( t; x_0, y_0 \\right)$. Our objective is to determine the evolution in time of the first-order derivatives of these solutions with respect to the initial conditions $\\left( x_0, y_0 \\right)$:\n", + "\n", + "$$\n", + "\\begin{aligned}\n", + "\\left(\\frac{\\partial x}{\\partial x_0}\\right)\\left(t; x_0, y_0\\right),\\\\\n", + "\\left(\\frac{\\partial x}{\\partial y_0}\\right)\\left(t; x_0, y_0\\right),\\\\\n", + "\\left(\\frac{\\partial y}{\\partial x_0}\\right)\\left(t; x_0, y_0\\right),\\\\\n", + "\\left(\\frac{\\partial y}{\\partial y_0}\\right)\\left(t; x_0, y_0\\right).\\\\\n", + "\\end{aligned}\n", + "$$\n", + "\n", + "These functions tell us how the the solutions of the ODE system react to changes in the initial values $\\left( x_0, y_0 \\right)$. One way of determining numerically the evolution in time of the partial derivatives is to introduce them as additional state variables in the original ODE system, together with an additional set of differential equations (i.e., the variational equations). Let's see how this is done for, e.g., $\\frac{\\partial x}{\\partial x_0}$. We take the time derivative of $\\frac{\\partial x}{\\partial x_0}$,\n", + "\n", + "$$\n", + "\\frac{d \\left( \\frac{\\partial x}{\\partial x_0} \\right)}{dt},\n", + "$$\n", + "\n", + "and expand it, via elementary calculus rules, to\n", + "\n", + "$$\n", + "\\frac{d \\left( \\frac{\\partial x}{\\partial x_0} \\right)}{dt} = \\frac{\\partial \\left( \\frac{dx}{dt} \\right) }{\\partial x_0}=\\frac{\\partial f\\left(x, y\\right)}{\\partial x}\\frac{\\partial x}{\\partial x_0}+\\frac{\\partial f \\left(x, y\\right)}{\\partial y}\\frac{\\partial y}{\\partial x_0}.\n", + "$$\n", + "\n", + "We can write similar ODEs for the other partial derivatives, append them below the original ODE system and finally obtain an augmented ODE system containing the variational equations:\n", + "\n", + "$$\n", + "\\begin{cases}\n", + "\\frac{dx}{dt} = f\\left(x, y \\right)\\\\\n", + "\\frac{dy}{dt} = g\\left(x, y \\right)\\\\\n", + "\\frac{d \\left( \\frac{\\partial x}{\\partial x_0} \\right)}{dt} = \\frac{\\partial f\\left(x, y\\right)}{\\partial x}\\frac{\\partial x}{\\partial x_0}+\\frac{\\partial f\\left(x, y\\right)}{\\partial y}\\frac{\\partial y}{\\partial x_0}\\\\\n", + "\\frac{d \\left( \\frac{\\partial x}{\\partial y_0} \\right)}{dt} = \\frac{\\partial f\\left(x, y\\right)}{\\partial x}\\frac{\\partial x}{\\partial y_0}+\\frac{\\partial f\\left(x, y\\right)}{\\partial y}\\frac{\\partial y}{\\partial y_0}\\\\\n", + "\\frac{d \\left( \\frac{\\partial y}{\\partial x_0} \\right)}{dt} = \\frac{\\partial g\\left(x, y\\right)}{\\partial x}\\frac{\\partial x}{\\partial x_0}+\\frac{\\partial g\\left(x, y\\right)}{\\partial y}\\frac{\\partial y}{\\partial x_0}\\\\\n", + "\\frac{d \\left( \\frac{\\partial y}{\\partial y_0} \\right)}{dt} = \\frac{\\partial g\\left(x, y\\right)}{\\partial x}\\frac{\\partial x}{\\partial y_0}+\\frac{\\partial g\\left(x, y\\right)}{\\partial y}\\frac{\\partial y}{\\partial y_0}\n", + "\\end{cases}.\n", + "$$\n", + "\n", + "The introduction of the variational equations in an ODE system of $n$ equations results in additional $n^2$ equations and variables. The initial conditions for the new variables are:\n", + "\n", + "$$\n", + "\\begin{cases}\n", + "\\frac{\\partial x}{\\partial x_0} \\left( 0 \\right) = 1\\\\\n", + "\\frac{\\partial x}{\\partial y_0} \\left( 0 \\right) = 0\\\\\n", + "\\frac{\\partial y}{\\partial x_0} \\left( 0 \\right) = 0\\\\\n", + "\\frac{\\partial y}{\\partial y_0} \\left( 0 \\right) = 1\\\\\n", + "\\end{cases}.\n", + "$$\n", + "\n", + "Let us now try to implement the variational equations for the simple pendulum:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "polyphonic-argentina", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Tolerance : 2.2204460492503131e-16\n", + "High accuracy : false\n", + "Compact mode : false\n", + "Taylor order : 20\n", + "Dimension : 6\n", + "Time : 0.0000000000000000\n", + "State : [0.050000000000000003, 0.025000000000000001, 1.0000000000000000, 0.0000000000000000, 0.0000000000000000, 1.0000000000000000]" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import heyoka as hy\n", + "\n", + "# Create the state variables x and v.\n", + "x, v = hy.make_vars(\"x\", \"v\")\n", + "\n", + "# Create the symbolic variables for the variational equations.\n", + "x_x0, x_y0, y_x0, y_y0 = hy.make_vars(\"x_x0\", \"x_y0\", \"y_x0\", \"y_y0\")\n", + "\n", + "# Create the right-hand side of the ODE system.\n", + "f = v\n", + "g = -9.8 * hy.sin(x)\n", + "x_x0_p = hy.diff(f, x) * x_x0 + hy.diff(f, v) * y_x0\n", + "x_y0_p = hy.diff(f, x) * x_y0 + hy.diff(f, v) * y_y0\n", + "y_x0_p = hy.diff(g, x) * x_x0 + hy.diff(g, v) * y_x0\n", + "y_y0_p = hy.diff(g, x) * x_y0 + hy.diff(g, v) * y_y0\n", + "\n", + "# Create the integrator object.\n", + "ta = hy.taylor_adaptive(\n", + " # Definition of the ODE system.\n", + " sys = [\n", + " (x, f),\n", + " (v, g),\n", + " (x_x0, x_x0_p),\n", + " (x_y0, x_y0_p),\n", + " (y_x0, y_x0_p),\n", + " (y_y0, y_y0_p),\n", + " ],\n", + " # Initial conditions.\n", + " state = [0.05, 0.025, 1., 0., 0., 1.])\n", + "\n", + "ta" + ] + }, + { + "cell_type": "markdown", + "id": "unlike-fifty", + "metadata": {}, + "source": [ + "Here we used the ``diff()`` function to perform the symbolic differentiation of ``f`` and ``g`` with respect to the state variables.\n", + "\n", + "Let us now propagate the state of the system (including the variational equations) up to $t=50$:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "recorded-shelf", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0.03744787, 0.10667026, 0.80315068, -0.17679098, 1.82916222,\n", + " 0.84245788])" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta.propagate_until(50.)\n", + "ta.state" + ] + }, + { + "cell_type": "markdown", + "id": "blessed-arabic", + "metadata": {}, + "source": [ + "Next, we are going to create another integrator for the simple pendulum, this time *without* the variational equations and perturbing the initial value of $x$ by $10^{-8}$ with respect to the original initial conditions:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "compound-installation", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(,\n", + " 0.20208180443700685,\n", + " 0.2182277296610072,\n", + " 240,\n", + " None)" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta_dx = hy.taylor_adaptive(\n", + " # Definition of the ODE system.\n", + " sys = [\n", + " (x, f),\n", + " (v, g)\n", + " ],\n", + " # Initial conditions.\n", + " state = [0.05 + 1e-8, 0.025])\n", + "\n", + "ta_dx.propagate_until(50.)" + ] + }, + { + "cell_type": "markdown", + "id": "piano-thinking", + "metadata": {}, + "source": [ + "Because of the perturbation on the initial value of $x$, the value of $x$ at $t=50$ will differ in the two integrators:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "express-pittsburgh", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "-8.031506786021492e-09" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ta.state[0] - ta_dx.state[0]" + ] + }, + { + "cell_type": "markdown", + "id": "relative-briefing", + "metadata": {}, + "source": [ + "However, because in the variational integrator we now have the value of $\\frac{\\partial x}{\\partial x_0}$ at $t=50$, we can compute an approximation of the $\\Delta x$ induced by the perturbation on the initial state as\n", + "\n", + "$$\n", + "\\Delta x \\approx \\frac{\\partial x}{\\partial x_0}\\left( 50 \\right) \\cdot 10^{-8}.\n", + "$$\n", + "\n", + "Indeed:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "confused-excellence", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2.0816681711721685e-17" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "(ta.state[0] + ta.state[2]*1e-8) - ta_dx.state[0]" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 7c88f37bf733bdb47a46a0c482925662b0074467 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Wed, 12 Jun 2024 08:19:51 +0200 Subject: [PATCH 21/23] Add deprecation notice to the old variational equations tutorial. --- doc/notebooks/The variational equations.ipynb | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/doc/notebooks/The variational equations.ipynb b/doc/notebooks/The variational equations.ipynb index 2bcf22fe..51f6c1b3 100644 --- a/doc/notebooks/The variational equations.ipynb +++ b/doc/notebooks/The variational equations.ipynb @@ -5,12 +5,17 @@ "id": "excited-uganda", "metadata": {}, "source": [ + "```{note}\n", + "\n", + "Starting with version 5.0.0, heyoka.py features builtin support for the formulation of the [variational equations](./var_ode_sys.ipynb). This tutorial is now deprecated and it will be removed in a future version.\n", + "\n", + "```\n", + "\n", "The variational equations\n", "=========================\n", "\n", "In this tutorial, we will show how it is possible to exploit heyoka.py's [expression system](<./The expression system.ipynb>) to implement first-order variational equations for a system of ODEs.\n", "\n", - "\n", "First, let's start with a brief recap on the variational equations. For simplicity, we will consider an autonomous ODE system in the variables $x$ and $y$ (the extension to more variables just complicates the notation):\n", "\n", "$$\n", From 768ebbd8429d0f40137c91f1717d88eee22fcfa7 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 13 Jun 2024 21:08:35 +0200 Subject: [PATCH 22/23] Disable linkcheck temporarily. --- tools/circleci_conda_heyoka_head_310.sh | 3 ++- tools/circleci_conda_heyoka_head_39.sh | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/tools/circleci_conda_heyoka_head_310.sh b/tools/circleci_conda_heyoka_head_310.sh index a92d60f0..499fde14 100644 --- a/tools/circleci_conda_heyoka_head_310.sh +++ b/tools/circleci_conda_heyoka_head_310.sh @@ -45,7 +45,8 @@ cd $HEYOKA_PY_PROJECT_DIR cd doc -make html linkcheck doctest +# make html linkcheck doctest +make html doctest set +e set +x diff --git a/tools/circleci_conda_heyoka_head_39.sh b/tools/circleci_conda_heyoka_head_39.sh index 8fa01aa9..c353beeb 100644 --- a/tools/circleci_conda_heyoka_head_39.sh +++ b/tools/circleci_conda_heyoka_head_39.sh @@ -45,7 +45,8 @@ cd $HEYOKA_PY_PROJECT_DIR cd doc -make html linkcheck doctest +#make html linkcheck doctest +make html doctest set +e set +x From a56eafed42ee764d8039116fe1ba9d05a26bc72d Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Thu, 13 Jun 2024 21:51:50 +0200 Subject: [PATCH 23/23] A couple of small test additions. --- heyoka/_test_var_integrator.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/heyoka/_test_var_integrator.py b/heyoka/_test_var_integrator.py index bb6a7376..d74bbdc0 100644 --- a/heyoka/_test_var_integrator.py +++ b/heyoka/_test_var_integrator.py @@ -47,6 +47,10 @@ def test_scalar(self): self.assertEqual(ts.shape, (2,)) self.assertTrue(np.all(ts == [0.0, 0.0])) + # Test that tstate cannot be written to. + with self.assertRaises(ValueError) as cm: + ta.tstate[0] = 0.5 + self.assertEqual(ta.get_vslice(order=0), slice(0, 2, None)) self.assertEqual(ta.get_vslice(order=0, component=1), slice(1, 2, None)) @@ -168,6 +172,10 @@ def test_batch(self): self.assertEqual(ts.shape, (2, 2)) self.assertTrue(np.all(ts == [[0.0, 0.0], [0.0, 0.0]])) + # Test that tstate cannot be written to. + with self.assertRaises(ValueError) as cm: + ta.tstate[0] = 0.5 + self.assertEqual(ta.get_vslice(order=0), slice(0, 2, None)) self.assertEqual(ta.get_vslice(order=0, component=1), slice(1, 2, None))