From d0198d122015651489d65a717362fad8cb2b7bf7 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Sun, 22 Oct 2023 16:11:32 +0200 Subject: [PATCH] Expose/test kepF and kepDE. --- CMakeLists.txt | 4 +- doc/changelog.rst | 9 +- heyoka/CMakeLists.txt | 2 + heyoka/_sympy_utils.py | 2 + heyoka/_test_celmec.py | 165 +++++++++++++++++ heyoka/_test_mp.py | 8 +- heyoka/_test_sympy.py | 325 ++++++++++++++++++++++++++++++++ heyoka/expose_expression.cpp | 191 +++++++++++++------ heyoka/setup_sympy.cpp | 10 +- heyoka/test.py | 347 +---------------------------------- 10 files changed, 657 insertions(+), 406 deletions(-) create mode 100644 heyoka/_test_celmec.py create mode 100644 heyoka/_test_sympy.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 1cfa4d8d..65447515 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,7 +11,7 @@ if(NOT CMAKE_BUILD_TYPE) FORCE) endif() -project(heyoka.py VERSION 3.0.0 LANGUAGES CXX C) +project(heyoka.py VERSION 3.1.0 LANGUAGES CXX C) list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" "${CMAKE_CURRENT_SOURCE_DIR}/cmake/yacma") @@ -118,7 +118,7 @@ find_package(fmt REQUIRED CONFIG) message(STATUS "fmt version: ${fmt_VERSION}") # heyoka. -find_package(heyoka 3.0.0 REQUIRED CONFIG) +find_package(heyoka 3.1.0 REQUIRED CONFIG) # Python. diff --git a/doc/changelog.rst b/doc/changelog.rst index 72bc6bca..37d4cfb8 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -3,9 +3,16 @@ Changelog ========= -3.0.1 (unreleased) +3.1.0 (unreleased) ------------------ +Changes +~~~~~~~ + +- heyoka.py now requires version 3.1.0 of the + heyoka C++ library + (`#140 `__). + Fix ~~~ diff --git a/heyoka/CMakeLists.txt b/heyoka/CMakeLists.txt index 736c4b3b..0fc87252 100644 --- a/heyoka/CMakeLists.txt +++ b/heyoka/CMakeLists.txt @@ -25,6 +25,8 @@ set(HEYOKA_PY_PYTHON_FILES _test_batch_integrator.py _test_ensemble.py _test_memcache.py + _test_celmec.py + _test_sympy.py model/__init__.py ) diff --git a/heyoka/_sympy_utils.py b/heyoka/_sympy_utils.py index ec0d7c8c..6be58131 100644 --- a/heyoka/_sympy_utils.py +++ b/heyoka/_sympy_utils.py @@ -177,6 +177,8 @@ def mul_wrapper(*args): retval[_spy.Mul] = mul_wrapper 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 return retval diff --git a/heyoka/_test_celmec.py b/heyoka/_test_celmec.py new file mode 100644 index 00000000..86573102 --- /dev/null +++ b/heyoka/_test_celmec.py @@ -0,0 +1,165 @@ +# Copyright 2020, 2021, 2022, 2023 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 kepE_test_case(_ut.TestCase): + def test_expr(self): + from . import kepE, diff, make_vars, sin, cos, core + from .core import _ppc_arch + import numpy as np + + x, y = make_vars("x", "y") + + # Try a few overloads. + kepE(x, y) + kepE(0.1, y) + kepE(e=x, M=0.2) + + self.assertEqual( + diff(kepE(x, y), x), sin(kepE(x, y)) / (1.0 - x * cos(kepE(x, y))) + ) + self.assertEqual(diff(kepE(x, y), y), 1.0 / (1.0 - x * cos(kepE(x, y)))) + + if not _ppc_arch: + self.assertEqual( + diff(kepE(x, np.longdouble("1.1")), x), + sin(kepE(x, np.longdouble("1.1"))) + / (1.0 - x * cos(kepE(x, np.longdouble("1.1")))), + ) + self.assertEqual( + diff(kepE(np.longdouble("1.1"), y), y), + 1.0 / (1.0 - np.longdouble("1.1") * cos(kepE(np.longdouble("1.1"), y))), + ) + + kepE(np.longdouble(0.1), y) + kepE(x, np.longdouble(0.2)) + + with self.assertRaises(TypeError) as cm: + kepE(0.1, np.longdouble("1.1")) + self.assertTrue( + "At least one of the arguments of kepE() must be an expression" + in str(cm.exception) + ) + + with self.assertRaises(TypeError) as cm: + kepE(0.1, 0.2) + self.assertTrue( + "At least one of the arguments of kepE() must be an expression" + in str(cm.exception) + ) + + if not hasattr(core, "real128"): + return + + from .core import real128 + + self.assertEqual( + diff(kepE(x, real128("1.1")), x), + sin(kepE(x, real128("1.1"))) / (1.0 - x * cos(kepE(x, real128("1.1")))), + ) + self.assertEqual( + diff(kepE(real128("1.1"), y), y), + 1.0 / (1.0 - real128("1.1") * cos(kepE(real128("1.1"), y))), + ) + + +class kepF_test_case(_ut.TestCase): + def test_expr(self): + from . import kepF, make_vars, core + from .core import _ppc_arch + import numpy as np + + x, y, z = make_vars("x", "y", "z") + + # Try a few overloads. + kepF(x, y, z) + kepF(0.1, y, z) + kepF(h=0.1, k=0.2, lam=z) + + if not _ppc_arch: + kepF(x, y, np.longdouble("1.1")) + kepF(x, np.longdouble(".1"), np.longdouble("1.1")) + + with self.assertRaises(TypeError) as cm: + kepF(x, 0.1, np.longdouble("1.1")) + self.assertTrue( + "The numerical arguments of kepF() must be all of the same type" + in str(cm.exception) + ) + + with self.assertRaises(TypeError) as cm: + kepF(0.1, 0.2, 0.3) + self.assertTrue( + "At least one of the arguments of kepF() must be an expression" + in str(cm.exception) + ) + + if not hasattr(core, "real128"): + return + + from .core import real128 + + kepF(real128(0.1), y, z) + kepF(real128(0.1), real128(0.2), z) + + with self.assertRaises(TypeError) as cm: + kepF(x, 0.1, real128("1.1")) + self.assertTrue( + "The numerical arguments of kepF() must be all of the same type" + in str(cm.exception) + ) + + +class kepDE_test_case(_ut.TestCase): + def test_expr(self): + from . import kepDE, make_vars, core + from .core import _ppc_arch + import numpy as np + + x, y, z = make_vars("x", "y", "z") + + # Try a few overloads. + kepDE(x, y, z) + kepDE(0.1, y, z) + kepDE(s0=0.1, c0=0.2, DM=z) + + if not _ppc_arch: + kepDE(x, y, np.longdouble("1.1")) + kepDE(x, np.longdouble(".1"), np.longdouble("1.1")) + + with self.assertRaises(TypeError) as cm: + kepDE(x, 0.1, np.longdouble("1.1")) + self.assertTrue( + "The numerical arguments of kepDE() must be all of the same type" + in str(cm.exception) + ) + + with self.assertRaises(TypeError) as cm: + kepDE(0.1, 0.2, 0.3) + self.assertTrue( + "At least one of the arguments of kepDE() must be an expression" + in str(cm.exception) + ) + + if not hasattr(core, "real128"): + return + + from .core import real128 + + kepDE(real128(0.1), y, z) + kepDE(real128(0.1), real128(0.2), z) + + with self.assertRaises(TypeError) as cm: + kepDE(x, 0.1, real128("1.1")) + self.assertTrue( + "The numerical arguments of kepDE() must be all of the same type" + in str(cm.exception) + ) diff --git a/heyoka/_test_mp.py b/heyoka/_test_mp.py index bd47585d..f95af169 100644 --- a/heyoka/_test_mp.py +++ b/heyoka/_test_mp.py @@ -368,7 +368,7 @@ def test_expression(self): if not hasattr(core, "real"): return - from . import expression as ex, real, kepE, atan2 + from . import expression as ex, real, kepE, kepF, kepDE, atan2 self.assertEqual( str(ex(real("1.1", 128))), "1.100000000000000000000000000000000000001" @@ -415,6 +415,12 @@ def test_expression(self): kepE(ex("x"), real("1.1", 128)) kepE(real("1.1", 128), ex("x")) + kepF(ex("x"), real("1.1", 128), ex("y")) + kepF(real("1.1", 128), ex("y"), ex("x")) + + kepDE(ex("x"), real("1.1", 128), ex("y")) + kepDE(real("1.1", 128), ex("y"), ex("x")) + atan2(ex("x"), real("1.1", 128)) atan2(real("1.1", 128), ex("x")) diff --git a/heyoka/_test_sympy.py b/heyoka/_test_sympy.py new file mode 100644 index 00000000..8f6115b3 --- /dev/null +++ b/heyoka/_test_sympy.py @@ -0,0 +1,325 @@ +# Copyright 2020, 2021, 2022, 2023 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 sympy_test_case(_ut.TestCase): + def test_basic(self): + try: + import sympy + except ImportError: + return + + from . import from_sympy, make_vars, sum as hsum + + with self.assertRaises(TypeError) as cm: + from_sympy(3.5) + self.assertTrue( + "The 'ex' parameter must be a sympy expression but it is of type" + in str(cm.exception) + ) + + with self.assertRaises(TypeError) as cm: + from_sympy(sympy.Symbol("x"), []) + self.assertTrue( + "The 's_dict' parameter must be a dict but it is of type" + in str(cm.exception) + ) + + with self.assertRaises(TypeError) as cm: + from_sympy(sympy.Symbol("x"), {3.5: 3.5}) + self.assertTrue( + "The keys in 's_dict' must all be sympy expressions" in str(cm.exception) + ) + + with self.assertRaises(TypeError) as cm: + from_sympy(sympy.Symbol("x"), {sympy.Symbol("x"): 3.5}) + self.assertTrue( + "The values in 's_dict' must all be heyoka expressions" in str(cm.exception) + ) + + # Test the s_dict functionality of from_sympy(). + x, y = sympy.symbols("x y", real=True) + hx, hy, hz = make_vars("x", "y", "z") + + self.assertEqual( + from_sympy((x - y) * (x + y), s_dict={x - y: hz}), hsum([hx, hy]) * hz + ) + + def test_number_conversion(self): + try: + import sympy + except ImportError: + return + + from . import to_sympy, from_sympy, expression, core + from .core import _ppc_arch + from sympy import Float, Rational, Integer + from mpmath import workprec + import numpy as np + + with self.assertRaises(ValueError) as cm: + from_sympy(Rational(3, 5)) + self.assertTrue( + "Cannot convert from sympy a rational number whose denominator is not a power of 2" + in str(cm.exception) + ) + + # From integer. + self.assertEqual(from_sympy(Integer(-42)), expression(-42.0)) + + # From rational. + self.assertEqual(from_sympy(Rational(42, -2)), expression(-21.0)) + + # Double precision. + with workprec(53): + self.assertEqual(to_sympy(expression(1.1)), Float(1.1)) + self.assertEqual(from_sympy(Float(1.1)), expression(1.1)) + + self.assertEqual( + to_sympy(expression((2**40 + 1) / (2**128))), + Rational(2**40 + 1, 2**128), + ) + self.assertEqual( + from_sympy(Rational(2**40 + 1, 2**128)), + expression((2**40 + 1) / (2**128)), + ) + + # Long double precision. + if not _ppc_arch: + with workprec(np.finfo(np.longdouble).nmant + 1): + self.assertEqual( + to_sympy(expression(np.longdouble("1.1"))), + Float("1.1", precision=np.finfo(np.longdouble).nmant + 1), + ) + + # NOTE: on platforms where long double is not wider than + # double (e.g., MSVC), conversion from sympy will produce a double + # and these tests will fail. + if np.finfo(np.longdouble).nmant > np.finfo(float).nmant: + self.assertEqual( + from_sympy(Float("1.1")), expression(np.longdouble("1.1")) + ) + + expo = np.finfo(np.longdouble).nmant - 10 + self.assertEqual( + to_sympy( + expression( + np.longdouble(2**expo + 1) / np.longdouble(2**128) + ) + ), + Rational(2**expo + 1, 2**128), + ) + self.assertEqual( + from_sympy(Rational(2**expo + 1, 2**128)), + expression( + np.longdouble(2**expo + 1) / np.longdouble(2**128) + ), + ) + + # Too high precision. + if not hasattr(core, "real"): + with self.assertRaises(ValueError) as cm: + from_sympy(Integer(2**500 + 1)) + self.assertTrue("the required precision" in str(cm.exception)) + + if not hasattr(core, "real128") or _ppc_arch: + return + + from .core import real128 + + # Quad precision. + with workprec(113): + self.assertEqual( + to_sympy(expression(real128("1.1"))), Float("1.1", precision=113) + ) + self.assertEqual(from_sympy(Float("1.1")), expression(real128("1.1"))) + + expo = 100 + self.assertEqual( + to_sympy(expression(real128(2**expo + 1) / real128(2**128))), + Rational(2**expo + 1, 2**128), + ) + self.assertEqual( + from_sympy(Rational(2**expo + 1, 2**128)), + expression(real128(2**expo + 1) / real128(2**128)), + ) + + def test_sympar_conversion(self): + try: + import sympy + except ImportError: + return + + from . import to_sympy, from_sympy, expression, par + from sympy import Symbol + + self.assertEqual(Symbol("x", real=True), to_sympy(expression("x"))) + self.assertEqual(Symbol("par[0]", real=True), to_sympy(par[0])) + self.assertEqual(Symbol("par[9]", real=True), to_sympy(par[9])) + self.assertEqual(Symbol("par[123]", real=True), to_sympy(par[123])) + self.assertEqual( + Symbol("par[-123]", real=True), to_sympy(expression("par[-123]")) + ) + self.assertEqual(Symbol("par[]", real=True), to_sympy(expression("par[]"))) + + self.assertEqual(from_sympy(Symbol("x")), expression("x")) + self.assertEqual(from_sympy(Symbol("par[0]")), par[0]) + self.assertEqual(from_sympy(Symbol("par[9]")), par[9]) + self.assertEqual(from_sympy(Symbol("par[123]")), par[123]) + self.assertEqual(from_sympy(Symbol("par[-123]")), expression("par[-123]")) + self.assertEqual(from_sympy(Symbol("par[]")), expression("par[]")) + + def test_func_conversion(self): + try: + import sympy + except ImportError: + return + + import sympy as spy + + from . import core, make_vars, from_sympy, to_sympy, pi, sum as hsum + + from .model import nbody + + x, y, z, a, b, c = spy.symbols("x y z a b c", real=True) + hx, hy, hz, ha, hb, hc = make_vars("x", "y", "z", "a", "b", "c") + + self.assertEqual(core.acos(hx), from_sympy(spy.acos(x))) + self.assertEqual(to_sympy(core.acos(hx)), spy.acos(x)) + + self.assertEqual(core.acosh(hx), from_sympy(spy.acosh(x))) + self.assertEqual(to_sympy(core.acosh(hx)), spy.acosh(x)) + + self.assertEqual(core.asin(hx), from_sympy(spy.asin(x))) + self.assertEqual(to_sympy(core.asin(hx)), spy.asin(x)) + + self.assertEqual(core.asinh(hx), from_sympy(spy.asinh(x))) + self.assertEqual(to_sympy(core.asinh(hx)), spy.asinh(x)) + + self.assertEqual(core.atan(hx), from_sympy(spy.atan(x))) + self.assertEqual(to_sympy(core.atan(hx)), spy.atan(x)) + + self.assertEqual(core.atan2(hy, hx), from_sympy(spy.atan2(y, x))) + self.assertEqual(to_sympy(core.atan2(hy, hx)), spy.atan2(y, x)) + + self.assertEqual(core.atanh(hx), from_sympy(spy.atanh(x))) + self.assertEqual(to_sympy(core.atanh(hx)), spy.atanh(x)) + + self.assertEqual(core.cos(hx), from_sympy(spy.cos(x))) + self.assertEqual(to_sympy(core.cos(hx)), spy.cos(x)) + + self.assertEqual(core.cosh(hx), from_sympy(spy.cosh(x))) + self.assertEqual(to_sympy(core.cosh(hx)), spy.cosh(x)) + + self.assertEqual(core.erf(hx), from_sympy(spy.erf(x))) + self.assertEqual(to_sympy(core.erf(hx)), spy.erf(x)) + + self.assertEqual(core.exp(hx), from_sympy(spy.exp(x))) + self.assertEqual(to_sympy(core.exp(hx)), spy.exp(x)) + + self.assertEqual(core.log(hx), from_sympy(spy.log(x))) + self.assertEqual(to_sympy(core.log(hx)), spy.log(x)) + + self.assertEqual(core.sin(hx), from_sympy(spy.sin(x))) + self.assertEqual(to_sympy(core.sin(hx)), spy.sin(x)) + + self.assertEqual(core.sinh(hx), from_sympy(spy.sinh(x))) + self.assertEqual(to_sympy(core.sinh(hx)), spy.sinh(x)) + + self.assertEqual(core.sqrt(hx), from_sympy(spy.sqrt(x))) + self.assertEqual(to_sympy(core.sqrt(hx)), spy.sqrt(x)) + + self.assertEqual(core.tan(hx), from_sympy(spy.tan(x))) + self.assertEqual(to_sympy(core.tan(hx)), spy.tan(x)) + + self.assertEqual(core.tanh(hx), from_sympy(spy.tanh(x))) + self.assertEqual(to_sympy(core.tanh(hx)), spy.tanh(x)) + + self.assertEqual(hx**3.5, from_sympy(x**3.5)) + self.assertEqual(to_sympy(hx**3.5), x**3.5) + + self.assertEqual(hsum([hx, hy, hz]), from_sympy(x + y + z)) + self.assertEqual(to_sympy(hx + hy + hz), x + y + z) + self.assertEqual(to_sympy(hsum([hx, hy, hz])), x + y + z) + self.assertEqual(to_sympy(hsum([hx])), x) + self.assertEqual(to_sympy(hsum([])), 0.0) + self.assertEqual( + hsum([ha, hb, hc, hx, hy, hz]), from_sympy(x + y + z + a + b + c) + ) + self.assertEqual(to_sympy(ha + hb + hc + hx + hy + hz), x + y + z + a + b + c) + self.assertEqual( + to_sympy(hsum([ha, hb, hc, hx, hy, hz])), x + y + z + a + b + c + ) + + self.assertEqual(hx * hy * hz, from_sympy(x * y * z)) + self.assertEqual(to_sympy(hx * hy * hz), x * y * z) + self.assertEqual( + (ha * hb) * (hc * hx) * (hy * hz), from_sympy(x * y * z * a * b * c) + ) + self.assertEqual(to_sympy(ha * hb * hc * hx * hy * hz), x * y * z * a * b * c) + + self.assertEqual(hsum([hx, -1.0 * hy, -1.0 * hz]), from_sympy(x - y - z)) + self.assertEqual(to_sympy(hx - hy - hz), x - y - z) + + # Run a test in the vector form as well. + self.assertEqual(to_sympy([hx - hy - hz, hx * hy * hz]), [x - y - z, x * y * z]) + + self.assertEqual(hx * hz**-1.0, from_sympy(x / z)) + self.assertEqual(to_sympy(hx / hz), x / z) + + self.assertEqual( + core.kepE(hx, hy), from_sympy(spy.Function("heyoka_kepE")(x, y)) + ) + self.assertEqual(to_sympy(core.kepE(hx, hy)), spy.Function("heyoka_kepE")(x, y)) + + self.assertEqual( + core.kepF(hx, hy, hz), from_sympy(spy.Function("heyoka_kepF")(x, y, z)) + ) + self.assertEqual( + to_sympy(core.kepF(hx, hy, hz)), spy.Function("heyoka_kepF")(x, y, z) + ) + + self.assertEqual( + core.kepDE(hx, hy, hz), from_sympy(spy.Function("heyoka_kepDE")(x, y, z)) + ) + self.assertEqual( + to_sympy(core.kepDE(hx, hy, hz)), spy.Function("heyoka_kepDE")(x, y, z) + ) + + self.assertEqual(-1.0 * hx, from_sympy(-x)) + self.assertEqual(to_sympy(-hx), -x) + + 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")()) + + with self.assertRaises(TypeError) as cm: + from_sympy(abs(x)) + self.assertTrue("Unable to convert the sympy object" in str(cm.exception)) + + # Test caching behaviour. + foo = hx + hy + bar = foo / (foo * hz + 1.0) + bar_spy = to_sympy(bar) + self.assertEqual( + id(bar_spy.args[1]), id(bar_spy.args[0].args[0].args[1].args[1]) + ) + + # pi constant. + self.assertEqual(to_sympy(pi), spy.pi) + self.assertEqual(from_sympy(spy.pi), pi) + self.assertEqual(to_sympy(from_sympy(spy.pi)), spy.pi) + + # nbody helper. + [to_sympy(_[1]) for _ in nbody(2)] + [to_sympy(_[1]) for _ in nbody(4)] + [to_sympy(_[1]) for _ in nbody(10)] diff --git a/heyoka/expose_expression.cpp b/heyoka/expose_expression.cpp index aafd8e08..2978b561 100644 --- a/heyoka/expose_expression.cpp +++ b/heyoka/expose_expression.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -50,6 +51,19 @@ namespace heyoka_py { +namespace detail +{ + +namespace +{ + +template +using uncvref_t = std::remove_cv_t>; + +} // namespace + +} // namespace detail + namespace py = pybind11; void expose_expression(py::module_ &m) @@ -263,9 +277,10 @@ void expose_expression(py::module_ &m) // Prod. m.def("prod", &hey::prod, "terms"_a); + // NOTE: need explicit casts for sqrt and exp due to the presence of overloads for number. m.def("sqrt", static_cast(&hey::sqrt)); + m.def("exp", static_cast(&hey::exp)); m.def("log", &hey::log); - m.def("exp", [](hey::expression e) { return hey::exp(std::move(e)); }); m.def("sin", &hey::sin); m.def("cos", &hey::cos); m.def("tan", &hey::tan); @@ -281,76 +296,134 @@ void expose_expression(py::module_ &m) m.def("sigmoid", &hey::sigmoid); m.def("erf", &hey::erf); - // kepE(). - m.def( - "kepE", [](hey::expression e, hey::expression M) { return hey::kepE(std::move(e), std::move(M)); }, "e"_a, - "M"_a); - m.def( - "kepE", [](double e, hey::expression M) { return hey::kepE(e, std::move(M)); }, "e"_a.noconvert(), "M"_a); - m.def( - "kepE", [](long double e, hey::expression M) { return hey::kepE(e, std::move(M)); }, "e"_a.noconvert(), "M"_a); + // NOTE: when exposing multivariate functions, we want to be able to pass + // in numerical arguments for convenience. Thus, we expose such functions taking + // in input a union of expression and supported numerical types. + using mvf_arg = std::variant; + // kepE(). m.def( - "kepE", [](hey::expression e, double M) { return hey::kepE(std::move(e), M); }, "e"_a, "M"_a.noconvert()); - m.def( - "kepE", [](hey::expression e, long double M) { return hey::kepE(std::move(e), M); }, "e"_a, "M"_a.noconvert()); -#if defined(HEYOKA_HAVE_REAL128) - m.def( - "kepE", [](hey::expression e, mppp::real128 M) { return hey::kepE(std::move(e), M); }, "e"_a, - "M"_a.noconvert()); -#endif -#if defined(HEYOKA_HAVE_REAL) - m.def( - "kepE", [](hey::expression e, mppp::real M) { return hey::kepE(std::move(e), std::move(M)); }, "e"_a, - "M"_a.noconvert()); -#endif + "kepE", + [](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; + + 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"); + } else { + return hey::kepE(a, b); + } + }, + e, M); + }, + "e"_a, "M"_a); - // atan2(). + // kepF(). m.def( - "atan2", [](hey::expression y, hey::expression x) { return hey::atan2(std::move(y), std::move(x)); }, "y"_a, - "x"_a); + "kepF", + [](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; + + constexpr auto tp1_num = static_cast(!std::is_same_v); + constexpr auto tp2_num = static_cast(!std::is_same_v); + constexpr auto tp3_num = static_cast(!std::is_same_v); + + constexpr auto n_num = tp1_num + tp2_num + tp3_num; + + if constexpr (n_num == 3) { + py_throw(PyExc_TypeError, "At least one of the arguments of kepF() must be an expression"); + } else if constexpr (n_num == 2) { + constexpr auto flag = tp1_num + (tp2_num << 1) + (tp3_num << 2); + + if constexpr (flag == 6 && std::is_same_v) { + return hey::kepF(a, b, c); + } else if constexpr (flag == 5 && std::is_same_v) { + return hey::kepF(a, b, c); + } else if constexpr (flag == 3 && std::is_same_v) { + return hey::kepF(a, b, c); + } else { + py_throw(PyExc_TypeError, "The numerical arguments of kepF() must be all of the same type"); + } + } else { + return hey::kepF(a, b, c); + } + }, + h, k, lam); + }, + "h"_a, "k"_a, "lam"_a); + // kepDE(). m.def( - "atan2", [](double y, hey::expression x) { return hey::atan2(y, std::move(x)); }, "y"_a.noconvert(), "x"_a); - m.def( - "atan2", [](long double y, hey::expression x) { return hey::atan2(y, std::move(x)); }, "y"_a.noconvert(), - "x"_a); -#if defined(HEYOKA_HAVE_REAL128) - m.def( - "atan2", [](mppp::real128 y, hey::expression x) { return hey::atan2(y, std::move(x)); }, "y"_a.noconvert(), - "x"_a); -#endif -#if defined(HEYOKA_HAVE_REAL) - m.def( - "atan2", [](mppp::real y, hey::expression x) { return hey::atan2(std::move(y), std::move(x)); }, - "y"_a.noconvert(), "x"_a); -#endif + "kepDE", + [](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; + + constexpr auto tp1_num = static_cast(!std::is_same_v); + constexpr auto tp2_num = static_cast(!std::is_same_v); + constexpr auto tp3_num = static_cast(!std::is_same_v); + + constexpr auto n_num = tp1_num + tp2_num + tp3_num; + + if constexpr (n_num == 3) { + py_throw(PyExc_TypeError, "At least one of the arguments of kepDE() must be an expression"); + } else if constexpr (n_num == 2) { + constexpr auto flag = tp1_num + (tp2_num << 1) + (tp3_num << 2); + + if constexpr (flag == 6 && std::is_same_v) { + return hey::kepDE(a, b, c); + } else if constexpr (flag == 5 && std::is_same_v) { + return hey::kepDE(a, b, c); + } else if constexpr (flag == 3 && std::is_same_v) { + return hey::kepDE(a, b, c); + } else { + py_throw(PyExc_TypeError, + "The numerical arguments of kepDE() must be all of the same type"); + } + } else { + return hey::kepDE(a, b, c); + } + }, + s0, c0, DM); + }, + "s0"_a, "c0"_a, "DM"_a); + // atan2(). m.def( - "atan2", [](hey::expression y, double x) { return hey::atan2(std::move(y), x); }, "y"_a, "x"_a.noconvert()); - m.def( - "atan2", [](hey::expression y, long double x) { return hey::atan2(std::move(y), x); }, "y"_a, - "x"_a.noconvert()); -#if defined(HEYOKA_HAVE_REAL128) - m.def( - "atan2", [](hey::expression y, mppp::real128 x) { return hey::atan2(std::move(y), x); }, "y"_a, - "x"_a.noconvert()); -#endif -#if defined(HEYOKA_HAVE_REAL) - m.def( - "atan2", [](hey::expression y, mppp::real x) { return hey::atan2(std::move(y), std::move(x)); }, "y"_a, - "x"_a.noconvert()); -#endif + "atan2", + [](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; + + 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"); + } else { + return hey::atan2(a, b); + } + }, + y, x); + }, + "y"_a, "x"_a); // Time. m.attr("time") = hey::time; diff --git a/heyoka/setup_sympy.cpp b/heyoka/setup_sympy.cpp index facbf850..ad515349 100644 --- a/heyoka/setup_sympy.cpp +++ b/heyoka/setup_sympy.cpp @@ -297,11 +297,17 @@ void setup_sympy(py::module &m) } }; - // kepE. - // NOTE: this will remain an unevaluated binary function. + // kepE, kepF, kepDE. + // NOTE: these will remain unevaluated functions. auto sympy_kepE = py::object(detail::spy->attr("Function")("heyoka_kepE")); detail::fmap[typeid(hy::detail::kepE_impl)] = sympy_kepE; + auto sympy_kepF = py::object(detail::spy->attr("Function")("heyoka_kepF")); + detail::fmap[typeid(hy::detail::kepF_impl)] = sympy_kepF; + + auto sympy_kepDE = py::object(detail::spy->attr("Function")("heyoka_kepDE")); + detail::fmap[typeid(hy::detail::kepDE_impl)] = sympy_kepDE; + // sigmoid. detail::fmap[typeid(hy::detail::sigmoid_impl)] = [](std::unordered_map &func_map, const hy::func &f) { diff --git a/heyoka/test.py b/heyoka/test.py index 20e378c0..8ad5e73f 100644 --- a/heyoka/test.py +++ b/heyoka/test.py @@ -1873,345 +1873,6 @@ def cb3(ta, mr, d_sgn): ) -class kepE_test_case(_ut.TestCase): - def test_expr(self): - from . import kepE, diff, make_vars, sin, cos, core - from .core import _ppc_arch - import numpy as np - - x, y = make_vars("x", "y") - self.assertEqual( - diff(kepE(x, y), x), sin(kepE(x, y)) / (1.0 - x * cos(kepE(x, y))) - ) - self.assertEqual(diff(kepE(x, y), y), 1.0 / (1.0 - x * cos(kepE(x, y)))) - - if not _ppc_arch: - self.assertEqual( - diff(kepE(x, np.longdouble("1.1")), x), - sin(kepE(x, np.longdouble("1.1"))) - / (1.0 - x * cos(kepE(x, np.longdouble("1.1")))), - ) - self.assertEqual( - diff(kepE(np.longdouble("1.1"), y), y), - 1.0 / (1.0 - np.longdouble("1.1") * cos(kepE(np.longdouble("1.1"), y))), - ) - - if not hasattr(core, "real128"): - return - - from .core import real128 - - self.assertEqual( - diff(kepE(x, real128("1.1")), x), - sin(kepE(x, real128("1.1"))) / (1.0 - x * cos(kepE(x, real128("1.1")))), - ) - self.assertEqual( - diff(kepE(real128("1.1"), y), y), - 1.0 / (1.0 - real128("1.1") * cos(kepE(real128("1.1"), y))), - ) - - -class sympy_test_case(_ut.TestCase): - def test_basic(self): - try: - import sympy - except ImportError: - return - - from . import from_sympy, make_vars, sum as hsum - - with self.assertRaises(TypeError) as cm: - from_sympy(3.5) - self.assertTrue( - "The 'ex' parameter must be a sympy expression but it is of type" - in str(cm.exception) - ) - - with self.assertRaises(TypeError) as cm: - from_sympy(sympy.Symbol("x"), []) - self.assertTrue( - "The 's_dict' parameter must be a dict but it is of type" - in str(cm.exception) - ) - - with self.assertRaises(TypeError) as cm: - from_sympy(sympy.Symbol("x"), {3.5: 3.5}) - self.assertTrue( - "The keys in 's_dict' must all be sympy expressions" in str(cm.exception) - ) - - with self.assertRaises(TypeError) as cm: - from_sympy(sympy.Symbol("x"), {sympy.Symbol("x"): 3.5}) - self.assertTrue( - "The values in 's_dict' must all be heyoka expressions" in str(cm.exception) - ) - - # Test the s_dict functionality of from_sympy(). - x, y = sympy.symbols("x y", real=True) - hx, hy, hz = make_vars("x", "y", "z") - - self.assertEqual( - from_sympy((x - y) * (x + y), s_dict={x - y: hz}), hsum([hx, hy]) * hz - ) - - def test_number_conversion(self): - try: - import sympy - except ImportError: - return - - from . import to_sympy, from_sympy, expression, core - from .core import _ppc_arch - from sympy import Float, Rational, Integer - from mpmath import workprec - import numpy as np - - with self.assertRaises(ValueError) as cm: - from_sympy(Rational(3, 5)) - self.assertTrue( - "Cannot convert from sympy a rational number whose denominator is not a power of 2" - in str(cm.exception) - ) - - # From integer. - self.assertEqual(from_sympy(Integer(-42)), expression(-42.0)) - - # From rational. - self.assertEqual(from_sympy(Rational(42, -2)), expression(-21.0)) - - # Double precision. - with workprec(53): - self.assertEqual(to_sympy(expression(1.1)), Float(1.1)) - self.assertEqual(from_sympy(Float(1.1)), expression(1.1)) - - self.assertEqual( - to_sympy(expression((2**40 + 1) / (2**128))), - Rational(2**40 + 1, 2**128), - ) - self.assertEqual( - from_sympy(Rational(2**40 + 1, 2**128)), - expression((2**40 + 1) / (2**128)), - ) - - # Long double precision. - if not _ppc_arch: - with workprec(np.finfo(np.longdouble).nmant + 1): - self.assertEqual( - to_sympy(expression(np.longdouble("1.1"))), - Float("1.1", precision=np.finfo(np.longdouble).nmant + 1), - ) - - # NOTE: on platforms where long double is not wider than - # double (e.g., MSVC), conversion from sympy will produce a double - # and these tests will fail. - if np.finfo(np.longdouble).nmant > np.finfo(float).nmant: - self.assertEqual( - from_sympy(Float("1.1")), expression(np.longdouble("1.1")) - ) - - expo = np.finfo(np.longdouble).nmant - 10 - self.assertEqual( - to_sympy( - expression( - np.longdouble(2**expo + 1) / np.longdouble(2**128) - ) - ), - Rational(2**expo + 1, 2**128), - ) - self.assertEqual( - from_sympy(Rational(2**expo + 1, 2**128)), - expression( - np.longdouble(2**expo + 1) / np.longdouble(2**128) - ), - ) - - # Too high precision. - if not hasattr(core, "real"): - with self.assertRaises(ValueError) as cm: - from_sympy(Integer(2**500 + 1)) - self.assertTrue("the required precision" in str(cm.exception)) - - if not hasattr(core, "real128") or _ppc_arch: - return - - from .core import real128 - - # Quad precision. - with workprec(113): - self.assertEqual( - to_sympy(expression(real128("1.1"))), Float("1.1", precision=113) - ) - self.assertEqual(from_sympy(Float("1.1")), expression(real128("1.1"))) - - expo = 100 - self.assertEqual( - to_sympy(expression(real128(2**expo + 1) / real128(2**128))), - Rational(2**expo + 1, 2**128), - ) - self.assertEqual( - from_sympy(Rational(2**expo + 1, 2**128)), - expression(real128(2**expo + 1) / real128(2**128)), - ) - - def test_sympar_conversion(self): - try: - import sympy - except ImportError: - return - - from . import to_sympy, from_sympy, expression, par - from sympy import Symbol - - self.assertEqual(Symbol("x", real=True), to_sympy(expression("x"))) - self.assertEqual(Symbol("par[0]", real=True), to_sympy(par[0])) - self.assertEqual(Symbol("par[9]", real=True), to_sympy(par[9])) - self.assertEqual(Symbol("par[123]", real=True), to_sympy(par[123])) - self.assertEqual( - Symbol("par[-123]", real=True), to_sympy(expression("par[-123]")) - ) - self.assertEqual(Symbol("par[]", real=True), to_sympy(expression("par[]"))) - - self.assertEqual(from_sympy(Symbol("x")), expression("x")) - self.assertEqual(from_sympy(Symbol("par[0]")), par[0]) - self.assertEqual(from_sympy(Symbol("par[9]")), par[9]) - self.assertEqual(from_sympy(Symbol("par[123]")), par[123]) - self.assertEqual(from_sympy(Symbol("par[-123]")), expression("par[-123]")) - self.assertEqual(from_sympy(Symbol("par[]")), expression("par[]")) - - def test_func_conversion(self): - try: - import sympy - except ImportError: - return - - import sympy as spy - - from . import core, make_vars, from_sympy, to_sympy, pi, sum as hsum - - from .model import nbody - - x, y, z, a, b, c = spy.symbols("x y z a b c", real=True) - hx, hy, hz, ha, hb, hc = make_vars("x", "y", "z", "a", "b", "c") - - self.assertEqual(core.acos(hx), from_sympy(spy.acos(x))) - self.assertEqual(to_sympy(core.acos(hx)), spy.acos(x)) - - self.assertEqual(core.acosh(hx), from_sympy(spy.acosh(x))) - self.assertEqual(to_sympy(core.acosh(hx)), spy.acosh(x)) - - self.assertEqual(core.asin(hx), from_sympy(spy.asin(x))) - self.assertEqual(to_sympy(core.asin(hx)), spy.asin(x)) - - self.assertEqual(core.asinh(hx), from_sympy(spy.asinh(x))) - self.assertEqual(to_sympy(core.asinh(hx)), spy.asinh(x)) - - self.assertEqual(core.atan(hx), from_sympy(spy.atan(x))) - self.assertEqual(to_sympy(core.atan(hx)), spy.atan(x)) - - self.assertEqual(core.atan2(hy, hx), from_sympy(spy.atan2(y, x))) - self.assertEqual(to_sympy(core.atan2(hy, hx)), spy.atan2(y, x)) - - self.assertEqual(core.atanh(hx), from_sympy(spy.atanh(x))) - self.assertEqual(to_sympy(core.atanh(hx)), spy.atanh(x)) - - self.assertEqual(core.cos(hx), from_sympy(spy.cos(x))) - self.assertEqual(to_sympy(core.cos(hx)), spy.cos(x)) - - self.assertEqual(core.cosh(hx), from_sympy(spy.cosh(x))) - self.assertEqual(to_sympy(core.cosh(hx)), spy.cosh(x)) - - self.assertEqual(core.erf(hx), from_sympy(spy.erf(x))) - self.assertEqual(to_sympy(core.erf(hx)), spy.erf(x)) - - self.assertEqual(core.exp(hx), from_sympy(spy.exp(x))) - self.assertEqual(to_sympy(core.exp(hx)), spy.exp(x)) - - self.assertEqual(core.log(hx), from_sympy(spy.log(x))) - self.assertEqual(to_sympy(core.log(hx)), spy.log(x)) - - self.assertEqual(core.sin(hx), from_sympy(spy.sin(x))) - self.assertEqual(to_sympy(core.sin(hx)), spy.sin(x)) - - self.assertEqual(core.sinh(hx), from_sympy(spy.sinh(x))) - self.assertEqual(to_sympy(core.sinh(hx)), spy.sinh(x)) - - self.assertEqual(core.sqrt(hx), from_sympy(spy.sqrt(x))) - self.assertEqual(to_sympy(core.sqrt(hx)), spy.sqrt(x)) - - self.assertEqual(core.tan(hx), from_sympy(spy.tan(x))) - self.assertEqual(to_sympy(core.tan(hx)), spy.tan(x)) - - self.assertEqual(core.tanh(hx), from_sympy(spy.tanh(x))) - self.assertEqual(to_sympy(core.tanh(hx)), spy.tanh(x)) - - self.assertEqual(hx**3.5, from_sympy(x**3.5)) - self.assertEqual(to_sympy(hx**3.5), x**3.5) - - self.assertEqual(hsum([hx, hy, hz]), from_sympy(x + y + z)) - self.assertEqual(to_sympy(hx + hy + hz), x + y + z) - self.assertEqual(to_sympy(hsum([hx, hy, hz])), x + y + z) - self.assertEqual(to_sympy(hsum([hx])), x) - self.assertEqual(to_sympy(hsum([])), 0.0) - self.assertEqual( - hsum([ha, hb, hc, hx, hy, hz]), from_sympy(x + y + z + a + b + c) - ) - self.assertEqual(to_sympy(ha + hb + hc + hx + hy + hz), x + y + z + a + b + c) - self.assertEqual( - to_sympy(hsum([ha, hb, hc, hx, hy, hz])), x + y + z + a + b + c - ) - - self.assertEqual(hx * hy * hz, from_sympy(x * y * z)) - self.assertEqual(to_sympy(hx * hy * hz), x * y * z) - self.assertEqual( - (ha * hb) * (hc * hx) * (hy * hz), from_sympy(x * y * z * a * b * c) - ) - self.assertEqual(to_sympy(ha * hb * hc * hx * hy * hz), x * y * z * a * b * c) - - self.assertEqual(hsum([hx, -1.0 * hy, -1.0 * hz]), from_sympy(x - y - z)) - self.assertEqual(to_sympy(hx - hy - hz), x - y - z) - - # Run a test in the vector form as well. - self.assertEqual(to_sympy([hx - hy - hz, hx * hy * hz]), [x - y - z, x * y * z]) - - self.assertEqual(hx * hz**-1.0, from_sympy(x / z)) - self.assertEqual(to_sympy(hx / hz), x / z) - - self.assertEqual( - core.kepE(hx, hy), from_sympy(spy.Function("heyoka_kepE")(x, y)) - ) - self.assertEqual(to_sympy(core.kepE(hx, hy)), spy.Function("heyoka_kepE")(x, y)) - - self.assertEqual(-1.0 * hx, from_sympy(-x)) - self.assertEqual(to_sympy(-hx), -x) - - 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")()) - - with self.assertRaises(TypeError) as cm: - from_sympy(abs(x)) - self.assertTrue("Unable to convert the sympy object" in str(cm.exception)) - - # Test caching behaviour. - foo = hx + hy - bar = foo / (foo * hz + 1.0) - bar_spy = to_sympy(bar) - self.assertEqual( - id(bar_spy.args[1]), id(bar_spy.args[0].args[0].args[1].args[1]) - ) - - # pi constant. - self.assertEqual(to_sympy(pi), spy.pi) - self.assertEqual(from_sympy(spy.pi), pi) - self.assertEqual(to_sympy(from_sympy(spy.pi)), spy.pi) - - # nbody helper. - [to_sympy(_[1]) for _ in nbody(2)] - [to_sympy(_[1]) for _ in nbody(4)] - [to_sympy(_[1]) for _ in nbody(10)] - - class llvm_state_test_case(_ut.TestCase): def test_copy(self): from . import make_vars, sin, taylor_adaptive @@ -2913,6 +2574,8 @@ def run_test_suite(): _test_batch_integrator, _test_ensemble, _test_memcache, + _test_celmec, + _test_sympy, ) import numpy as np from .model import nbody @@ -2947,8 +2610,10 @@ def run_test_suite(): suite.addTest(tl.loadTestsFromTestCase(llvm_state_test_case)) suite.addTest(tl.loadTestsFromTestCase(event_classes_test_case)) suite.addTest(tl.loadTestsFromTestCase(event_detection_test_case)) - suite.addTest(tl.loadTestsFromTestCase(kepE_test_case)) - suite.addTest(tl.loadTestsFromTestCase(sympy_test_case)) + suite.addTest(tl.loadTestsFromTestCase(_test_celmec.kepE_test_case)) + suite.addTest(tl.loadTestsFromTestCase(_test_celmec.kepF_test_case)) + suite.addTest(tl.loadTestsFromTestCase(_test_celmec.kepDE_test_case)) + suite.addTest(tl.loadTestsFromTestCase(_test_sympy.sympy_test_case)) suite.addTest(tl.loadTestsFromTestCase(_test_memcache.memcache_test_case)) test_result = _ut.TextTestRunner(verbosity=2).run(suite)