diff --git a/.editorconfig b/.editorconfig index 8d3905d..e269037 100644 --- a/.editorconfig +++ b/.editorconfig @@ -7,6 +7,7 @@ root = true [*] end_of_line = lf insert_final_newline = true +trim_trailing_whitespace = false # Matches multiple files with brace expansion notation # Set default charset diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1df0059..b29faa9 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -14,6 +14,20 @@ concurrency: cancel-in-progress: true jobs: + typos: + name: Spelling (typos) + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: crate-ci/typos@master + + ruff: + name: Linting (ruff) + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: chartboost/ruff-action@v1 + test: name: Unittest (${{ matrix.os }}-${{ matrix.compiler }}-py${{ matrix.python-version }}) runs-on: ${{ matrix.os }} diff --git a/docs/conf.py b/docs/conf.py index a305e49..431b03f 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,42 +1,48 @@ -import sys import os +import sys from unittest.mock import MagicMock as Mock + from setuptools_scm import get_version + MOCK_MODULES = [ - 'numpy', 'numpy.testing', 'numpy.random', - 'scipy', 'scipy.integrate', 'scipy.integrate._ode', 'scipy.stats', 'scipy.integrate._ivp.ivp', - 'symengine', 'symengine.printing', 'symengine.lib.symengine_wrapper', - 'jitcxde_common.helpers','jitcxde_common.numerical','jitcxde_common.symbolic','jitcxde_common.transversal', + "numpy", "numpy.testing", "numpy.random", + "scipy", "scipy.integrate", "scipy.integrate._ode", "scipy.stats", "scipy.integrate._ivp.ivp", + "symengine", "symengine.printing", "symengine.lib.symengine_wrapper", + "jitcxde_common.helpers","jitcxde_common.numerical","jitcxde_common.symbolic","jitcxde_common.transversal", ] sys.modules.update((mod_name, Mock()) for mod_name in MOCK_MODULES) -class ode_mock(object): pass -sys.modules['scipy.integrate'] = Mock(ode=ode_mock) -class GroupHandler_mock(object): pass -sys.modules['jitcxde_common.transversal'] = Mock(GroupHandler=GroupHandler_mock) +class ode_mock: + pass + +class GroupHandler_mock: + pass + +sys.modules["scipy.integrate"] = Mock(ode=ode_mock) +sys.modules["jitcxde_common.transversal"] = Mock(GroupHandler=GroupHandler_mock) sys.path.insert(0,os.path.abspath("../examples")) sys.path.insert(0,os.path.abspath("../jitcode")) -needs_sphinx = '1.3' +needs_sphinx = "1.3" extensions = [ - 'sphinx.ext.autodoc', - 'sphinx.ext.autosummary', - 'sphinx.ext.mathjax', - 'numpydoc', - 'sphinx.ext.graphviz' + "sphinx.ext.autodoc", + "sphinx.ext.autosummary", + "sphinx.ext.mathjax", + "sphinx.ext.graphviz" + "numpydoc", ] -source_suffix = '.rst' +source_suffix = ".rst" -master_doc = 'index' +master_doc = "index" -project = u'JiTCODE' -copyright = u'2016, Gerrit Ansmann' +project = "JiTCODE" +copyright = "2016, Gerrit Ansmann" -release = version = get_version(root='..', relative_to=__file__) +release = version = get_version(root="..", relative_to=__file__) default_role = "any" @@ -44,22 +50,22 @@ class GroupHandler_mock(object): pass add_module_names = False -html_theme = 'nature' -pygments_style = 'colorful' -htmlhelp_basename = 'JiTCODEdoc' +html_theme = "nature" +pygments_style = "colorful" +htmlhelp_basename = "JiTCODEdoc" numpydoc_show_class_members = False -autodoc_member_order = 'bysource' +autodoc_member_order = "bysource" graphviz_output_format = "svg" -toc_object_entries_show_parents = 'hide' +toc_object_entries_show_parents = "hide" def on_missing_reference(app, env, node, contnode): - if node['reftype'] == 'any': + if node["reftype"] == "any": return contnode else: return None def setup(app): - app.connect('missing-reference', on_missing_reference) + app.connect("missing-reference", on_missing_reference) diff --git a/examples/SW_of_Roesslers.py b/examples/SW_of_Roesslers.py index eac9a08..1401ffe 100644 --- a/examples/SW_of_Roesslers.py +++ b/examples/SW_of_Roesslers.py @@ -1,5 +1,4 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- """ This example showcases several advanced features of JiTCODE that are relevant for an efficient integration of more complex systems as well as how to deal with some special situations. Therefore it is pretty bizarre from a dynamical perspective. @@ -13,7 +12,7 @@ \\dot{z}_i &= b + z_i (x_i -c) + k \\sum_{j=0}^N (z_j-z_i) \\end{alignedat} -The control parameters shall be :math:`a = 0.165`, :math:`b = 0.2`, :math:`c = 10.0`, and :math:`k = 0.01`. The (frequency) parameter :math:`ω_i` shall be picked randomly from the uniform distribution on :math:`[0.8,1.0]` for each :math:`i`. :math:`A∈ℝ^{N×N}` shall be the adjacency matrix of a one-dimensional small-world network (which shall be provided by a function `small_world_network` in the following example code). So, the :math:`x` compenents are coupled diffusively with a small-world coupling topology, while the :math:`z` components are coupled diffusively to their mean field, with the coupling term being modulated with :math:`\\sin(t)`. +The control parameters shall be :math:`a = 0.165`, :math:`b = 0.2`, :math:`c = 10.0`, and :math:`k = 0.01`. The (frequency) parameter :math:`ω_i` shall be picked randomly from the uniform distribution on :math:`[0.8,1.0]` for each :math:`i`. :math:`A∈ℝ^{N×N}` shall be the adjacency matrix of a one-dimensional small-world network (which shall be provided by a function `small_world_network` in the following example code). So, the :math:`x` components are coupled diffusively with a small-world coupling topology, while the :math:`z` components are coupled diffusively to their mean field, with the coupling term being modulated with :math:`\\sin(t)`. Without further ado, here is the example code (`complete running example `_); highlighted lines will be commented upon below: @@ -47,10 +46,10 @@ def small_world_network(number_of_nodes, nearest_neighbours, rewiring_probabilit # rewiring for i in range(n): for j in range(i): - if A[i,j] and (np.random.random() < rewiring_probability): + if A[i,j] and (rng.random() < rewiring_probability): A[j,i] = A[i,j] = False while True: - i_new,j_new = np.random.randint(0,n,2) + i_new,j_new = rng.integers(0,n,2) if A[i_new,j_new] or i_new==j_new: continue else: @@ -60,15 +59,17 @@ def small_world_network(number_of_nodes, nearest_neighbours, rewiring_probabilit return A # example-start - from jitcode import jitcode, y, t import numpy as np import symengine + + from jitcode import jitcode, t, y # parameters # ---------- + rng = np.random.default_rng() N = 500 - ω = np.random.uniform(0.8,1.0,N) + ω = rng.uniform(0.8,1.0,N) a = 0.165 b = 0.2 c = 10.0 @@ -99,11 +100,11 @@ def f(): # integrate # --------- - initial_state = np.random.random(3*N) + initial_state = rng.random(3*N) ODE = jitcode(f, helpers=helpers, n=3*N) ODE.generate_f_C(simplify=False, do_cse=False, chunk_size=150) - ODE.set_integrator('dopri5') + ODE.set_integrator("dopri5") ODE.set_initial_value(initial_state,0.0) # data structure: x[0], v[0], z[0], x[1], …, x[N], v[N], z[N] diff --git a/examples/benchmark.py b/examples/benchmark.py index 38efbf6..ac29330 100644 --- a/examples/benchmark.py +++ b/examples/benchmark.py @@ -1,16 +1,18 @@ -import numpy as np -from numpy.random import choice, uniform from time import process_time -from scipy.integrate import ode, solve_ivp, odeint + +import numpy as np +from scipy.integrate import ode, odeint, solve_ivp from scipy.integrate._ivp.ivp import METHODS -from jitcode import jitcode, y from symengine import sin +from jitcode import jitcode, y + + solver_ode = "dopri5" solver_ivp = "RK45" # Context manager for timing -class timer(object): +class timer: def __init__(self,name): self.name = name @@ -20,105 +22,105 @@ def __enter__(self): def __exit__(self,*args): end = process_time() duration = end-self.start - print("%s took %.5f s" % (self.name,duration)) + print(f"{self.name} took {duration:.5f}s") # The actual test def test_scenario(name,fun,initial,times,rtol,atol): print(40*"-",name,40*"-",sep="\n") - with timer("ode (%s)"%solver_ode): - I = ode(fun) - I.set_integrator(solver_ode,rtol=rtol,atol=atol,nsteps=10**8) - I.set_initial_value(initial,0.0) - result = np.vstack([I.integrate(time) for time in times]) - assert I.successful() + with timer(f"ode ({solver_ode})"): + solver = ode(fun) + solver.set_integrator(solver_ode,rtol=rtol,atol=atol,nsteps=10**8) + solver.set_initial_value(initial,0.0) + _result = np.vstack([solver.integrate(time) for time in times]) + assert solver.successful() inv_fun = lambda y,t: fun(t,y) with timer("odeint with suboptimal function (LSODA)"): - result = odeint( + _result = odeint( func=inv_fun, - y0=initial, t=[0.0]+list(times), + y0=initial, t=[0.0, *times], rtol=rtol, atol=atol, mxstep=10**8 ) - with timer("solve_ivp (%s) without result"%solver_ivp): - I = solve_ivp( + with timer(f"solve_ivp ({solver_ivp}) without result"): + solver = solve_ivp( fun, t_span=(0,times[-1]), y0=initial, method=solver_ivp, rtol=rtol, atol=atol ) - assert I.status != -1 + assert solver.status != -1 - with timer("solve_ivp (%s)"%solver_ivp): - I = solve_ivp( + with timer(f"solve_ivp ({solver_ivp})"): + solver = solve_ivp( fun, t_span=(0,times[-1]), t_eval=times, y0=initial, method=solver_ivp, rtol=rtol, atol=atol ) - result = I.y - assert I.status != -1 + _result = solver.y + assert solver.status != -1 - with timer("solve_ivp (%s) with dense_output"%solver_ivp): - I = solve_ivp( + with timer(f"solve_ivp ({solver_ivp}) with dense_output"): + solver = solve_ivp( fun, t_span=(0,times[-1]), y0=initial, method=solver_ivp, rtol=rtol, atol=atol, dense_output=True ) - result = np.vstack([I.sol(time) for time in times]) - assert I.status != -1 + _result = np.vstack([solver.sol(time) for time in times]) + assert solver.status != -1 - with timer("%s with dense output"%solver_ivp): - I = METHODS[solver_ivp]( + with timer(f"{solver_ivp} with dense output"): + solver = METHODS[solver_ivp]( fun=fun, y0=initial, t0=0.0, t_bound=times[-1], rtol=rtol, atol=atol ) def solutions(): for time in times: - while I.t < time: - I.step() - yield I.dense_output()(time) - result = np.vstack(list(solutions())) - assert I.status != "failed" + while solver.t < time: + solver.step() + yield solver.dense_output()(time) + _result = np.vstack(list(solutions())) + assert solver.status != "failed" - with timer("%s with manual resetting"%solver_ivp): - I = METHODS[solver_ivp]( + with timer(f"{solver_ivp} with manual resetting"): + solver = METHODS[solver_ivp]( fun=fun, y0=initial, t0=0.0, t_bound=times[-1], rtol=rtol, atol=atol ) def solutions(): for time in times: - I.t_bound = time - I.status = "running" - while I.status == "running": - I.step() - yield I.y - result = np.vstack(list(solutions())) - assert I.status != "failed" + solver.t_bound = time + solver.status = "running" + while solver.status == "running": + solver.step() + yield solver.y + _result = np.vstack(list(solutions())) + assert solver.status != "failed" - with timer("%s with reinitialising"%solver_ivp): + with timer(f"{solver_ivp} with reinitialising"): def solutions(): current_time = 0.0 state = initial for time in times: - I = METHODS[solver_ivp]( + solver = METHODS[solver_ivp]( fun=fun, y0=state, t0=current_time, t_bound=time, rtol=rtol, atol=atol ) - while I.status == "running": - I.step() - assert I.status != "failed" + while solver.status == "running": + solver.step() + assert solver.status != "failed" current_time = time - state = I.y + state = solver.y yield state - result = np.vstack(list(solutions())) + _result = np.vstack(list(solutions())) # Using compiled functions to make things faster def get_compiled_function(f): @@ -142,9 +144,12 @@ def get_compiled_function(f): ) +rng = np.random.default_rng() + n, c, q = 100, 3.0, 0.2 -A = choice( [1,0], size=(n,n), p=[q,1-q] ) -omega = uniform(-0.5,0.5,n) +A = rng.choice( [1,0], size=(n,n), p=[q,1-q] ) +omega = rng.uniform(-0.5,0.5,n) + def kuramotos_f(): for i in range(n): coupling_sum = sum( @@ -157,10 +162,8 @@ def kuramotos_f(): test_scenario( name = "random network of Kuramoto oscillators", fun = get_compiled_function(kuramotos_f), - initial = uniform(0,2*np.pi,n), + initial = rng.uniform(0,2*np.pi,n), times = range(1,10001,10), rtol = 1e-13, atol = 1e-6, ) - - diff --git a/examples/double_fhn_example.py b/examples/double_fhn_example.py index 9b7c5df..eb12361 100644 --- a/examples/double_fhn_example.py +++ b/examples/double_fhn_example.py @@ -1,7 +1,6 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- -""" +r""" Suppose our differential equation is :math:`\dot{y} = f(y)` with :math:`y∈ℝ^4`, .. math:: @@ -20,8 +19,9 @@ if __name__ == "__main__": # example-start - from jitcode import jitcode, y import numpy as np + + from jitcode import jitcode, y a = -0.025794 b1 = 0.0065 diff --git a/examples/double_fhn_lyapunov.py b/examples/double_fhn_lyapunov.py index c73555c..4f5caf0 100644 --- a/examples/double_fhn_lyapunov.py +++ b/examples/double_fhn_lyapunov.py @@ -1,5 +1,4 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- """ For instance, we can calculate and print the Lyapunov exponents for the system from `example` as follows (changes highlighted): @@ -7,9 +6,10 @@ if __name__ == "__main__": # example-start - from jitcode import jitcode_lyap, y - from scipy.stats import sem import numpy as np + from scipy.stats import sem + + from jitcode import jitcode_lyap, y a = -0.025794 b1 = 0.0065 @@ -42,4 +42,4 @@ for i in range(n): lyap = np.average(lyaps[1000:,i]) stderr = sem(lyaps[1000:,i]) # Note that this only an estimate - print("%i. Lyapunov exponent: % .4f ± %.4f" % (i+1,lyap,stderr)) + print(f"{i+1}. Lyapunov exponent: {lyap:.4f} ± {stderr:.4f}") diff --git a/examples/double_fhn_restricted_lyap.py b/examples/double_fhn_restricted_lyap.py index f76d555..6f88efc 100644 --- a/examples/double_fhn_restricted_lyap.py +++ b/examples/double_fhn_restricted_lyap.py @@ -1,10 +1,13 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- -from jitcode import jitcode_restricted_lyap, y import numpy as np from scipy.stats import sem +from jitcode import jitcode_restricted_lyap, y + + +rng = np.random.default_rng() + a = -0.025794 b1 = 0.01 b2 = 0.01 @@ -18,7 +21,7 @@ b2*y(2) - c*y(3) ] -initial_state = np.random.random(4) +initial_state = rng.random(4) vectors = [ np.array([1.,0.,1.,0.]), diff --git a/examples/double_fhn_transversal_lyap.py b/examples/double_fhn_transversal_lyap.py index 0afb4e1..dedeb43 100644 --- a/examples/double_fhn_transversal_lyap.py +++ b/examples/double_fhn_transversal_lyap.py @@ -1,5 +1,4 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- """ For instance, let’s interpret the system from `example` as two oscillators (which is what it is), one consisting of the first and second and one of the third and fourth component. Furthermore, let’s change the control parameters a bit to make the two oscillators identical. We can then calculate the transversal Lyapunov exponents to the synchronisation manifold as follows (important changes are highlighted): @@ -15,9 +14,10 @@ if __name__ == "__main__": # example-start - from jitcode import jitcode_transversal_lyap, y - from scipy.stats import sem import numpy as np + from scipy.stats import sem + + from jitcode import jitcode_transversal_lyap, y a = -0.025794 b = 0.01 @@ -45,5 +45,4 @@ lyap = np.average(lyaps[1000:]) stderr = sem(lyaps[1000:]) # Note that this only an estimate - print("transversal Lyapunov exponent: % .4f ± %.4f" % (lyap,stderr)) - + print(f"transversal Lyapunov exponent: {lyap:.4f} ± {stderr:.4f}") diff --git a/examples/kuramoto_network.py b/examples/kuramoto_network.py index c34479a..8a9e194 100755 --- a/examples/kuramoto_network.py +++ b/examples/kuramoto_network.py @@ -1,15 +1,18 @@ -from numpy import pi, vstack, savetxt -from jitcode import jitcode, y +import numpy as np from symengine import sin -from numpy.random import choice, uniform + +from jitcode import jitcode, y + + +rng = np.random.default_rng() n = 100 c = 3.0 q = 0.2 -A = choice( [1,0], size=(n,n), p=[q,1-q] ) +A = rng.choice( [1,0], size=(n,n), p=[q,1-q] ) -omega = uniform(-0.5,0.5,n) +omega = rng.uniform(-0.5,0.5,n) # Sorting does not effect the qualitative dynamics, but only the resulting plots: omega.sort() @@ -22,14 +25,14 @@ def kuramotos_f(): ) yield omega[i] + c/(n-1)*coupling_sum -I = jitcode(kuramotos_f,n=n) -I.set_integrator("dop853",atol=1e-6,rtol=0) +solver = jitcode(kuramotos_f,n=n) +solver.set_integrator("dop853",atol=1e-6,rtol=0) -initial_state = uniform(0,2*pi,n) -I.set_initial_value(initial_state,time=0.0) +initial_state = rng.uniform(0,2*np.pi,n) +solver.set_initial_value(initial_state,time=0.0) times = range(0,2001) -data = vstack ([ I.integrate(time) for time in times ]) -data %= 2*pi -savetxt("kuramotos.dat",data) +data = np.vstack ([ solver.integrate(time) for time in times ]) +data %= 2*np.pi +np.savetxt("kuramotos.dat",data) diff --git a/examples/lotka_volterra.py b/examples/lotka_volterra.py index 357562e..5a1dee7 100644 --- a/examples/lotka_volterra.py +++ b/examples/lotka_volterra.py @@ -1,5 +1,4 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- """ Suppose, we want to implement the Lotka–Volterra model, which is described by the following equations: @@ -90,8 +89,9 @@ if __name__ == "__main__": # example-start - from jitcode import y, jitcode import numpy as np + + from jitcode import jitcode, y γ = 0.6 φ = 1.0 diff --git a/jitcode/__init__.py b/jitcode/__init__.py index b77aed5..cab0ab8 100644 --- a/jitcode/__init__.py +++ b/jitcode/__init__.py @@ -1,12 +1,9 @@ -from ._jitcode import ( - t, y, - jitcode, jitcode_lyap, jitcode_restricted_lyap, jitcode_transversal_lyap, - UnsuccessfulIntegration, - test, - ) +from ._jitcode import jitcode, jitcode_lyap, jitcode_restricted_lyap, jitcode_transversal_lyap, t, test, y # noqa: F401 +from .integrator_tools import UnsuccessfulIntegration # noqa: F401 + try: - from .version import version as __version__ + from .version import version as __version__ # noqa: F401 except ImportError: from warnings import warn - warn('Failed to find (autogenerated) version.py. Do not worry about this unless you really need to know the version.') + warn("Failed to find (autogenerated) version.py. Do not worry about this unless you really need to know the version.", stacklevel=1) diff --git a/jitcode/_jitcode.py b/jitcode/_jitcode.py index 874d3e1..d8efc81 100644 --- a/jitcode/_jitcode.py +++ b/jitcode/_jitcode.py @@ -1,22 +1,20 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- -from warnings import warn -from types import FunctionType, BuiltinFunctionType -from inspect import signature from itertools import count +from types import BuiltinFunctionType, FunctionType +from warnings import warn -from numpy import hstack, log import numpy as np import symengine -from jitcxde_common import jitcxde, checker -from jitcxde_common.helpers import sympify_helpers, sort_helpers, find_dependent_helpers -from jitcxde_common.numerical import random_direction, orthonormalise_qr +from jitcxde_common import checker, jitcxde +from jitcxde_common.helpers import find_dependent_helpers, sort_helpers, sympify_helpers +from jitcxde_common.numerical import orthonormalise_qr, random_direction from jitcxde_common.symbolic import collect_arguments, ordered_subs, replace_function from jitcxde_common.transversal import GroupHandler -from jitcode.integrator_tools import empty_integrator, IVP_wrapper, IVP_wrapper_no_interpolation, ODE_wrapper, integrator_info, UnsuccessfulIntegration +from jitcode.integrator_tools import IVP_wrapper, IVP_wrapper_no_interpolation, ODE_wrapper, empty_integrator, integrator_info + #: the symbol for the state that must be used to define the differential equation. It is a function and the integer argument denotes the component. You may just as well define an analogous function directly with SymEngine or SymPy, but using this function is the best way to get the most of future versions of JiTCODE, in particular avoiding incompatibilities. You can import a SymPy variant from the submodule `sympy_symbols` instead (see `SymPy vs. SymEngine`_ for details). y = symengine.Function("y") @@ -72,7 +70,7 @@ class jitcode(jitcxde): * A SymEngine function object used in `f_sym` to represent the function call. If you want to use any JiTCODE features that need the derivative, this must have a properly defined `f_diff` method with the derivative being another callback function (or constant). * The Python function to be called. This function will receive the state array (`y`) as the first argument. All further arguments are whatever you use as arguments of the SymEngine function in `f_sym`. These can be any expression that you might use in the definition of the derivative and contain, e.g., dynamical variables, time, control parameters, and helpers. The only restriction is that the arguments are floats (and not vectors or similar). The return value must also be a float (or something castable to float). It is your responsibility to ensure that this function adheres to these criteria, is deterministic and sufficiently smooth with respect its arguments; expect nasty errors otherwise. - * The number of arguments, **excluding** the state array as mandatory first argument. This means if you have a variadic Python function, you cannot just call it with different numbers of arguments in `f_sym`, but you have to define separate callbacks for each of numer of arguments. + * The number of arguments, **excluding** the state array as mandatory first argument. This means if you have a variadic Python function, you cannot just call it with different numbers of arguments in `f_sym`, but you have to define separate callbacks for each of number of arguments. See `this example `_ (for JiTCDDE) for how to use this. @@ -137,11 +135,11 @@ def _check_valid_arguments(self): for argument in collect_arguments(entry,y): self._check_assert( argument[0] >= 0, - "y is called with a negative argument (%i) in equation %i." % (argument[0],i), + f"y is called with a negative argument ({argument[0]}) in equation {i}.", ) self._check_assert( argument[0] < self.n, - "y is called with an argument (%i) higher than the system’s dimension (%i) in equation %i." % (argument[0],self.n,i) + f"y is called with an argument ({argument[0]}) higher than the system’s dimension ({self.n}) in equation {i}." ) @checker @@ -152,7 +150,7 @@ def _check_valid_symbols(self): for symbol in entry.atoms(symengine.Symbol): self._check_assert( symbol in valid_symbols, - "Invalid symbol (%s) in equation %i." % (symbol.name,i), + f"Invalid symbol ({symbol.name}) in equation {i}.", ) @property @@ -245,7 +243,7 @@ def generate_f_C(self, simplify=None, do_cse=False, chunk_size=100): (set_dy(i,entry) for i,entry in enumerate(f_sym_wc)), name = "f", chunk_size = chunk_size, - arguments = arguments+[("dY", "PyArrayObject *__restrict const")] + arguments = [*arguments, ("dY", "PyArrayObject *__restrict const")] ) self._f_C_source = True @@ -321,7 +319,7 @@ def generate_jac_C(self, do_cse=False, chunk_size=100, sparse=True): ), name = "jac", chunk_size = chunk_size, - arguments = arguments+[("dfdY", "PyArrayObject *__restrict const")] + arguments = [*arguments, ("dfdY", "PyArrayObject *__restrict const")] ) self._jac_C_source = True @@ -355,7 +353,7 @@ def generate_helpers_C(self, chunk_size=100): (set_helper(i, helper[1].subs(self.general_subs)) for i,helper in enumerate(self.helpers)), name = "general_helpers", chunk_size = chunk_size, - arguments = self._default_arguments() + [("general_helper","double *__restrict const")], + arguments = [*self._default_arguments(), ("general_helper","double *__restrict const")], omp = False, ) @@ -424,12 +422,12 @@ def _prepare_lambdas(self): if not hasattr(self,"_lambda_subs") or not hasattr(self,"_lambda_args"): if self.helpers: - warn("Lambdification handles helpers by plugging them in. This may be very inefficient") + warn("Lambdification handles helpers by plugging them in. This may be very inefficient", stacklevel=3) self._lambda_subs = list(reversed(self.helpers)) self._lambda_args = [t] for i in range(self.n): - symbol = symengine.Symbol("dummy_argument_%i"%i) + symbol = symengine.Symbol(f"dummy_argument_{i}") self._lambda_args.append(symbol) self._lambda_subs.append((y(i),symbol)) self._lambda_args.extend(self.control_pars) @@ -582,7 +580,7 @@ def set_integrator(self,name,nsteps=10**6,interpolate=True,**integrator_params): The `solve_ivp` methods are usually slightly faster for large differential equations, but they come with a massive overhead that makes them considerably slower for small differential equations. Implicit solvers are slower than explicit ones, except for stiff problems. If you don’t know what to choose, start with `"dopri5"`. nsteps: integer - Same as the respective parameter of the `ode` solvers, but with a higher default value to avoid annoying errors when getting rid of transients. + Same as the respective parameter of the `ode` solvers, but with a higher default value to avoid annoying errors when getting rid of transients. interpolate: boolean Whether the sampled solutions for `solve_ivp` solvers shall be obtained using interpolation. If your sampling step is small, this may make things faster; otherwise it depends. This may also make the results slightly less accurate. @@ -657,11 +655,11 @@ def set_parameters(self,*args): self.initialise() def set_f_params(self, *args): - warn("This function has been replaced by `set_parameters`") + warn("This function has been replaced by `set_parameters`", stacklevel=2) self.set_parameters(*args) def set_jac_params(self, *args): - warn("This function has been replaced by `set_parameters`") + warn("This function has been replaced by `set_parameters`", stacklevel=2) self.set_parameters(*args) def integrate(self,*args,**kwargs): @@ -689,7 +687,7 @@ def __init__( self, f_sym=(), n_lyap=-1, simplify=None, **kwargs ): f_basic = self._handle_input(f_sym,n_basic=True) self._n_lyap = n_lyap if (0<=n_lyap<=self.n_basic) else self.n_basic if self._n_lyap>10: - warn("You are about to calculate %i Lyapunov exponents. This is very likely more than you need and may lead to severe difficulties with compilation and integration. Unless you really know what you are doing, consider how many Lyapunov exponents you actually need and set the parameter `n_lyap` accordingly." % self._n_lyap) + warn(f"You are about to calculate {self._n_lyap} Lyapunov exponents. This is very likely more than you need and may lead to severe difficulties with compilation and integration. Unless you really know what you are doing, consider how many Lyapunov exponents you actually need and set the parameter `n_lyap` accordingly.", stacklevel=2) if simplify is None: simplify = self.n_basic<=10 @@ -715,7 +713,7 @@ def f_lyap(): expression = expression.simplify(ratio=1.0) yield expression - super(jitcode_lyap, self).__init__( + super().__init__( f_lyap, helpers = helpers, n = self.n_basic*(self._n_lyap+1), @@ -730,7 +728,7 @@ def set_initial_value(self, y, time=0.0): for _ in range(self._n_lyap): new_y.append(random_direction(self.n_basic)) - super(jitcode_lyap, self).set_initial_value(hstack(new_y), time) + super().set_initial_value(np.hstack(new_y), time) def set_integrator(self,name,interpolate=None,**kwargs): """ @@ -739,15 +737,15 @@ def set_integrator(self,name,interpolate=None,**kwargs): if interpolate is None: interpolate = name in ["RK45","Radau"] if name == "LSODA": - warn("Using LSODA for Lyapunov exponents is discouraged since interpolation errors may accumulate.") - super(jitcode_lyap,self).set_integrator(name,interpolate,**kwargs) + warn("Using LSODA for Lyapunov exponents is discouraged since interpolation errors may accumulate.", stacklevel=2) + super().set_integrator(name,interpolate,**kwargs) def norms(self): n = self.n_basic tangent_vectors = self._y[n:].reshape(self._n_lyap,n) tangent_vectors,norms = orthonormalise_qr(tangent_vectors) if not np.all(np.isfinite(norms)): - warn("Norms of perturbation vectors for Lyapunov exponents out of numerical bounds. You probably waited too long before renormalising and should call integrate with smaller intervals between steps (as renormalisations happen once with every call of integrate).") + warn("Norms of perturbation vectors for Lyapunov exponents out of numerical bounds. You probably waited too long before renormalising and should call integrate with smaller intervals between steps (as renormalisations happen once with every call of integrate).", stacklevel=2) return norms, tangent_vectors def integrate(self, *args, **kwargs): @@ -767,12 +765,12 @@ def integrate(self, *args, **kwargs): """ old_t = self.t - super(jitcode_lyap, self).integrate(*args, **kwargs) + super().integrate(*args, **kwargs) delta_t = self.t-old_t norms, tangent_vectors = self.norms() - lyaps = log(norms) / delta_t + lyaps = np.log(norms) / delta_t self._y[self.n_basic:] = tangent_vectors.flatten() - super(jitcode_lyap, self).set_initial_value(self._y, self.t) + super().set_initial_value(self._y, self.t) return self._y[:self.n_basic], lyaps, tangent_vectors @@ -844,7 +842,7 @@ def finalise(entry): def f_lyap(): for entry in self.iterate(tangent_vector_f()): - if type(entry)==int: # i.e., if main index + if type(entry) is int: # i.e., if main index if average_dynamics: group = groups[entry] yield sum( finalise(f_list[i]) for i in group )/len(group) @@ -855,7 +853,7 @@ def f_lyap(): helpers = ((helper[0],finalise(helper[1])) for helper in helpers) - super(jitcode_transversal_lyap, self).__init__( + super().__init__( f_lyap, helpers = helpers, n = self.n, @@ -873,7 +871,7 @@ def set_initial_value(self, y, time=0.0): new_y = np.empty(self.n) new_y[self.main_indices] = y new_y[self.tangent_indices] = random_direction(len(self.tangent_indices)) - super(jitcode_transversal_lyap, self).set_initial_value(new_y, time) + super().set_initial_value(new_y, time) def set_integrator(self,name,interpolate=False,**kwargs): """ @@ -882,15 +880,15 @@ def set_integrator(self,name,interpolate=False,**kwargs): if interpolate is None: interpolate = name in ["RK45","Radau"] if name == "LSODA": - warn("Using LSODA for Lyapunov exponents is discouraged since interpolation errors may accumulate.") - super(jitcode_transversal_lyap,self).set_integrator(name,interpolate,**kwargs) + warn("Using LSODA for Lyapunov exponents is discouraged since interpolation errors may accumulate.", stacklevel=2) + super().set_integrator(name,interpolate,**kwargs) def norm(self): tangent_vector = self._y[self.tangent_indices] norm = np.linalg.norm(tangent_vector) tangent_vector /= norm if not np.isfinite(norm): - warn("Norm of perturbation vector for Lyapunov exponent out of numerical bounds. You probably waited too long before renormalising and should call integrate with smaller intervals between steps (as renormalisations happen once with every call of integrate).") + warn("Norm of perturbation vector for Lyapunov exponent out of numerical bounds. You probably waited too long before renormalising and should call integrate with smaller intervals between steps (as renormalisations happen once with every call of integrate).", stacklevel=2) self._y[self.tangent_indices] = tangent_vector return norm @@ -908,11 +906,11 @@ def integrate(self, *args, **kwargs): """ old_t = self.t - super(jitcode_transversal_lyap, self).integrate(*args, **kwargs) + super().integrate(*args, **kwargs) delta_t = self.t-old_t norm = self.norm() - lyap = log(norm) / delta_t - super(jitcode_transversal_lyap, self).set_initial_value(self._y, self.t) + lyap = np.log(norm) / delta_t + super().set_initial_value(self._y, self.t) return self._y[self.main_indices], lyap @@ -933,7 +931,7 @@ class jitcode_restricted_lyap(jitcode_lyap): def __init__(self, f_sym=(), vectors=(), **kwargs): kwargs["n_lyap"] = 1 - super(jitcode_restricted_lyap,self).__init__(f_sym,**kwargs) + super().__init__(f_sym,**kwargs) self.vectors = [ vector/np.linalg.norm(vector) for vector in vectors ] def norms(self): @@ -944,7 +942,7 @@ def norms(self): norm = np.linalg.norm(tangent_vector) tangent_vector /= norm if not np.isfinite(norm): - warn("Norm of perturbation vector for Lyapunov exponent out of numerical bounds. You probably waited too long before renormalising and should call integrate with smaller intervals between steps (as renormalisations happen once with every call of integrate).") + warn("Norm of perturbation vector for Lyapunov exponent out of numerical bounds. You probably waited too long before renormalising and should call integrate with smaller intervals between steps (as renormalisations happen once with every call of integrate).", stacklevel=2) return norm, tangent_vector def test(omp=True,sympy=True): @@ -965,4 +963,3 @@ def test(omp=True,sympy=True): ODE.set_integrator("dopri5") ODE.set_initial_value([1,2]) ODE.integrate(0.1) - diff --git a/jitcode/integrator_tools.py b/jitcode/integrator_tools.py index 858c253..7b2fa57 100644 --- a/jitcode/integrator_tools.py +++ b/jitcode/integrator_tools.py @@ -1,9 +1,10 @@ from inspect import signature +import numpy as np from scipy.integrate import ode -from scipy.integrate._ode import find_integrator from scipy.integrate._ivp.ivp import METHODS as ivp_methods -from numpy import inf +from scipy.integrate._ode import find_integrator + class UnsuccessfulIntegration(Exception): """ @@ -15,7 +16,7 @@ def integrator_info(name): """ Finds out the integrator from a given name, what backend it uses, and whether it can use a Jacobian. """ - if name == 'zvode': + if name == "zvode": raise NotImplementedError("JiTCODE does not natively support complex numbers yet.") if name in ivp_methods.keys(): @@ -35,7 +36,7 @@ def integrator_info(name): "integrator": integrator } -class IVP_wrapper(object): +class IVP_wrapper: """ This is a wrapper around the integrators from scipy.integrate.solve_ivp making them work like scipy.integrate.ode (mostly). """ @@ -56,7 +57,7 @@ def __init__( # Dictionary to be passed as arguments to the integrator and store stuff self.kwargs = { - "t_bound": inf, + "t_bound": np.inf, "vectorized": False, "fun": self.f } @@ -132,7 +133,7 @@ class ODE_wrapper(ode): """ def integrate(self,t,step=False,relax=False): if t>self.t or step or relax: - result = super(ODE_wrapper,self).integrate(t,step,relax) + result = super().integrate(t,step,relax) if self.successful(): return result else: @@ -149,9 +150,9 @@ def params(self): def set_params(self,*args): raise NotImplementedError("This method should not be called anymore") -class empty_integrator(object): +class empty_integrator: """ - This is a dummy class that mimicks some basic properties of scipy.integrate.ode or the above wrappers, respectively. It exists to store states and parameters and to raise exceptions in the same interface. + This is a dummy class that mimics some basic properties of scipy.integrate.ode or the above wrappers, respectively. It exists to store states and parameters and to raise exceptions in the same interface. """ def __init__(self): @@ -181,4 +182,3 @@ def integrate(self,*args,**kwargs): def successful(self): raise RuntimeError("You must call set_integrator first.") - diff --git a/jitcode/sympy_symbols.py b/jitcode/sympy_symbols.py index b97c692..65d9688 100644 --- a/jitcode/sympy_symbols.py +++ b/jitcode/sympy_symbols.py @@ -1,5 +1,6 @@ import sympy + #: the symbol for the state that must be used to define the differential equation. It is a function and the integer argument denotes the component. This one is different from the one you can import from jitcode directly by being defined via SymPy and thus being better suited for some symbolic processing techniques that are not available in SymEngine yet. y = sympy.Function("y",real=True) diff --git a/pyproject.toml b/pyproject.toml index 2b9a3d4..1976822 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -54,3 +54,47 @@ jitcode = [ [tool.setuptools_scm] write_to = "jitcode/version.py" + +[tool.ruff] +target-version = "py37" +line-length = 320 + +[tool.ruff.lint] +select = [ + "A", # flake8-builtins + "B", # flake8-bugbear + "C4", # flake8-comprehensions + "E", # flake8 + "F", # flake8 + "I", # flake8-isort + "NPY", # numpy + "Q", # flake8-quotes + "RUF", # ruff + "UP", # pyupgrade + "W", # flake8 +] +ignore = [ + "A001", # shadowing-builtin-variable + "A005", # shadowing-builtin-module + "C409", # incorrectly-parenthesized-tuple-in-subscript + "E203", # whitespace-before-punctuation + "E402", # module-import-not-at-top-of-file + "E501", # line-too-long + "E731", # assign-to-lambda + "RUF001", # ambiguous-unicode-character-string + "RUF002", # ambiguous-unicode-character-docstring + "RUF003", # ambiguous-unicode-character-comment + "W191", # indentation-contains-tabs + "W293", # blank-line-with-whitespace +] + +[tool.ruff.lint.flake8-quotes] +docstring-quotes = "double" +inline-quotes = "double" +multiline-quotes = "double" + +[tool.ruff.lint.isort] +combine-as-imports = true +known-local-folder = [ "jitcode", "scenarios" ] +known-first-party = [ "jitcxde_common" ] +lines-after-imports = 2 diff --git a/tests/scenarios.py b/tests/scenarios.py index bb4e177..7744079 100644 --- a/tests/scenarios.py +++ b/tests/scenarios.py @@ -1,10 +1,12 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- + +from random import shuffle import numpy as np -from symengine import symbols, Function +from symengine import Function, symbols + from jitcode import y -from random import shuffle + n = 4 @@ -27,9 +29,9 @@ ]) y1 = np.array([ - 0.0011789485114731, + 0.0011789485114731, -0.0021947158873226, - 0.0195744683782066, + 0.0195744683782066, -0.0057801623466600, ]) @@ -59,8 +61,7 @@ # with generator def f_generator(): - for entry in f: - yield entry + yield from f with_generator = { "f_sym":f_generator, "n":n } @@ -132,5 +133,3 @@ def first_component(Y): (call_c_times,c_times,1) ], } - - diff --git a/tests/test_generation.py b/tests/test_generation.py index 6c8c774..8df0529 100644 --- a/tests/test_generation.py +++ b/tests/test_generation.py @@ -1,23 +1,18 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- """ -Tests that the code-generation and compilaton or lambdification, respectively, of the derivative and Jacobian works as intended in all kinds of scenarios. +Tests that the code-generation and compilation or lambdification, respectively, of the derivative and Jacobian works as intended in all kinds of scenarios. """ import unittest -from numpy.random import random +import numpy as np from numpy.testing import assert_allclose -from jitcode import jitcode, jitcode_lyap, y +from jitcode import jitcode, jitcode_lyap from jitcode._jitcode import _is_C, _is_lambda +from scenarios import callback, f_of_y0, jac_of_y0, n, params_args, vanilla, with_dictionary, with_generator, with_helpers, with_params, y0 -from scenarios import ( - y0, f_of_y0, jac_of_y0, - vanilla, with_params, with_helpers, with_generator, with_dictionary, callback, - n, params_args - ) class TestBasic(unittest.TestCase): init_params = () @@ -60,7 +55,7 @@ def tearDown(self): self.assertIsNotNone(self.ODE.f) self.ODE.set_parameters(*self.init_params) assert_allclose( self.ODE.f(0.0,y0), f_of_y0, rtol=1e-5 ) - if not self.ODE.jac is None: + if self.ODE.jac is not None: assert_allclose( self.ODE.jac(0.0,y0), jac_of_y0, rtol=1e-5) class TestHelpers(TestBasic): @@ -69,7 +64,8 @@ def setUp(self): class FurtherHelpersTests(unittest.TestCase): def test_identity_of_jacs(self): - x = random(n) + rng = np.random.default_rng() + x = rng.random(n) def evaluate(scenario): ODE = jitcode(**scenario) @@ -87,7 +83,8 @@ def evaluate(scenario): ) def test_identity_of_lyaps(self): - x = random((n+1)*n) + rng = np.random.default_rng() + x = rng.random((n+1)*n) def evaluate(scenario): ODE = jitcode_lyap(**scenario,n_lyap=n) @@ -130,4 +127,3 @@ def test_lambdas(self): if __name__ == "__main__": unittest.main(buffer=True) - diff --git a/tests/test_integrator_tools.py b/tests/test_integrator_tools.py index 3030d12..0fcecb2 100644 --- a/tests/test_integrator_tools.py +++ b/tests/test_integrator_tools.py @@ -1,16 +1,16 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- import unittest import numpy as np from numpy.testing import assert_allclose -from numpy.random import random from jitcode import jitcode -from jitcode.integrator_tools import empty_integrator, IVP_wrapper, IVP_wrapper_no_interpolation, ODE_wrapper, UnsuccessfulIntegration +from jitcode.integrator_tools import IVP_wrapper, IVP_wrapper_no_interpolation, ODE_wrapper, UnsuccessfulIntegration, empty_integrator +from scenarios import n, vanilla, y0, y1 -from scenarios import y0, y1, vanilla, n + +rng = np.random.default_rng() # Generating compiled functions @@ -30,9 +30,9 @@ # ----------------------------- -class TestSkeleton(object): +class TestSkeleton: """ - This class exists to be inherited by a test that adds self.initialise to intialise self.integrator. + This class exists to be inherited by a test that adds self.initialise to initialise self.integrator. """ def control_result(self): @@ -46,13 +46,13 @@ def test_no_params(self): def test_initial_twice(self): self.initialise(f,jac,rtol=1e-5) - self.integrator.set_initial_value(random(n)) + self.integrator.set_initial_value(rng.random(n)) self.integrator.set_initial_value(y0) self.control_result() def test_zero_integration(self): self.initialise(f,jac) - initial = random(n) + initial = rng.random(n) self.integrator.set_initial_value(initial) assert_allclose(initial,self.integrator.integrate(0)) @@ -143,7 +143,7 @@ def setUp(self): def test_t(self): with self.assertRaises(RuntimeError): - self.integrator.t + _ = self.integrator.t def test_set_integrator(self): with self.assertRaises(RuntimeError): @@ -154,11 +154,10 @@ def test_integrate(self): self.integrator.integrate(2.3) def test_set_initial(self): - initial = random(5) + initial = rng.random(5) self.integrator.set_initial_value(initial,1.2) assert np.all( self.integrator._y == initial ) assert self.integrator.t == 1.2 if __name__ == "__main__": unittest.main(buffer=True) - diff --git a/tests/test_jitcode.py b/tests/test_jitcode.py index d035386..8f5f622 100644 --- a/tests/test_jitcode.py +++ b/tests/test_jitcode.py @@ -1,5 +1,4 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- import unittest @@ -8,12 +7,10 @@ from scipy.stats import sem as standard_error from symengine import symbols -from jitcode import jitcode, y, jitcode_lyap, UnsuccessfulIntegration, test +from jitcode import jitcode, jitcode_lyap, y from jitcode._jitcode import _is_C, _is_lambda +from scenarios import f_of_y0, jac_of_y0, lyaps, n, vanilla, y0, y1 -from scenarios import ( - y0, f_of_y0, jac_of_y0, y1, lyaps, vanilla, n, - ) class TestOrders(unittest.TestCase): """ @@ -39,7 +36,7 @@ def test_initialise_first(self): self.ODE.set_integrator("dop853") self.assertTrue(_is_C(self.ODE.f)) - def test_initalise_with_dict(self): + def test_initialise_with_dict(self): self.ODE = jitcode(**vanilla) initial_value = {y(i):y0[i] for i in range(n)} self.ODE.set_initial_value(initial_value) @@ -91,7 +88,7 @@ def tearDown(self): self.ODE.check() self.assertIsNotNone(self.ODE.f) assert_allclose( self.ODE.f(0.0,y0), f_of_y0, rtol=1e-5 ) - if not self.ODE.jac is None: + if self.ODE.jac is not None: assert_allclose( self.ODE.jac(0.0,y0), jac_of_y0, rtol=1e-5) assert_allclose( self.ODE.integrate(1.0), y1, rtol=1e-4 ) for i in reversed(range(n)): @@ -174,7 +171,7 @@ def test_duplicate_error(self): def test_wrong_n(self): with self.assertRaises(ValueError): - ODE = jitcode(**vanilla,n=2*n) + _ODE = jitcode(**vanilla,n=2*n) def test_dimension_mismatch(self): ODE = jitcode(**vanilla) @@ -209,4 +206,3 @@ def test_no_interpolation_LSODA(self): if __name__ == "__main__": unittest.main(buffer=True) - diff --git a/tests/test_order.py b/tests/test_order.py index 682a064..cd24206 100644 --- a/tests/test_order.py +++ b/tests/test_order.py @@ -1,13 +1,16 @@ -from jitcode import jitcode -from symengine import Symbol from itertools import permutations -p = Symbol('p') +from symengine import Symbol + +from jitcode import jitcode + + +p = Symbol("p") f = [1/p] for integrator in ["dopri5","RK45"]: - def set_integrator(ODE): + def set_integrator(ODE, integrator=integrator): ODE.set_integrator(integrator) def set_parameters(ODE): @@ -34,7 +37,7 @@ def set_initial_value(ODE): f_2 = [1] for integrator in ["dopri5","RK45"]: - def set_integrator(ODE): + def set_integrator(ODE, integrator=integrator): ODE.set_integrator(integrator) def set_initial_value(ODE): diff --git a/tests/test_sympy_input.py b/tests/test_sympy_input.py index 263dede..421042e 100644 --- a/tests/test_sympy_input.py +++ b/tests/test_sympy_input.py @@ -2,10 +2,12 @@ Tests whether things works independent of where symbols are imported from. """ +import symengine +import sympy + import jitcode import jitcode.sympy_symbols -import sympy -import symengine + symengine_manually = [ symengine.Symbol("t",real=True), diff --git a/tests/test_transversal_lyap.py b/tests/test_transversal_lyap.py index f6bd61b..d2e5554 100644 --- a/tests/test_transversal_lyap.py +++ b/tests/test_transversal_lyap.py @@ -1,17 +1,22 @@ #!/usr/bin/python3 -# -*- coding: utf-8 -*- """ Integration test of jitcode_restricted_lyap and jitcode_transversal_lyap by comparing their results to each other for a synchronised scenario. """ from itertools import combinations -from jitcode import jitcode_restricted_lyap, y, jitcode_transversal_lyap -from jitcxde_common import DEFAULT_COMPILE_ARGS + import numpy as np from scipy.stats import sem from symengine import Symbol +from jitcxde_common import DEFAULT_COMPILE_ARGS + +from jitcode import jitcode_restricted_lyap, jitcode_transversal_lyap, y + + +rng = np.random.default_rng() + a = -0.025794 b = 0.01 c = 0.02 @@ -80,7 +85,7 @@ ) # Simplification or compiler optimisation would lead to trajectories diverging from the synchronisation manifold due to numerical noise. ODE1.generate_f_C(simplify=False) - ODE1.compile_C( extra_compile_args = DEFAULT_COMPILE_ARGS + ["-O2"] ) + ODE1.compile_C( extra_compile_args = [*DEFAULT_COMPILE_ARGS, "-O2"] ) ODE1.set_integrator("dopri5") ODE2 = jitcode_transversal_lyap( @@ -97,23 +102,23 @@ ODE2.set_parameters(coupling["k"]) if coupling["sign"]<0: - initial_state = np.random.random(n) + initial_state = rng.random(n) else: - single = np.random.random(2) + single = rng.random(2) initial_state = np.empty(n) for j,group in enumerate(scenario["groups"]): for i in group: initial_state[i] = single[j] ODE1.set_initial_value(initial_state,0.0) - ODE2.set_initial_value(np.random.random(2),0.0) + ODE2.set_initial_value(rng.random(2),0.0) times = range(100,100000,100) lyaps1 = np.hstack([ODE1.integrate(time)[1] for time in times]) lyaps2 = np.hstack([ODE2.integrate(time)[1] for time in times]) # Check that we are still on the synchronisation manifold: - message = f"The dynamics left the synchronisation manifold when {scenario['name']} with coupling {coupling}. If this fails, this is a problem with the test and not with what is tested or any software involved.\n\nSpecifically, this test only works when the backend (Symengine plus compiler) implents certain computations completely symmetrically. This needs not and cannot be reasonably controlled (and no, turning off compiler optimisation doesn’t necessarily help as it often restores symmetries broken by Symengine). It’s only something exploited by this test to make it work in the first place." + message = f"The dynamics left the synchronisation manifold when {scenario['name']} with coupling {coupling}. If this fails, this is a problem with the test and not with what is tested or any software involved.\n\nSpecifically, this test only works when the backend (Symengine plus compiler) implements certain computations completely symmetrically. This needs not and cannot be reasonably controlled (and no, turning off compiler optimisation doesn’t necessarily help as it often restores symmetries broken by Symengine). It’s only something exploited by this test to make it work in the first place." for group in scenario["groups"]: for i,j in combinations(group,2): assert ODE1.y[i]==ODE1.y[j], message @@ -126,7 +131,7 @@ sign2 = np.sign(Lyap2) if abs(Lyap2)>margin2 else 0 assert sign1==coupling["sign"] assert sign2==coupling["sign"] - assert abs(Lyap1-Lyap2)