From 01d5ca769531cace52bfc3e6cfafe83519a6a0b0 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 15 Jul 2025 10:10:51 +0300 Subject: [PATCH 1/5] Fix latex in fhn example --- docs/index.rst | 2 +- examples/double_fhn_example.py | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/index.rst b/docs/index.rst index 009b925..f8767f0 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -40,7 +40,7 @@ As with SciPy’s ODE, this documentation assumes that the differential equation .. math:: - \dot{y} = f(t,y) + \dot{y} = f(t,y). .. _example: diff --git a/examples/double_fhn_example.py b/examples/double_fhn_example.py index eb12361..b18f737 100644 --- a/examples/double_fhn_example.py +++ b/examples/double_fhn_example.py @@ -5,13 +5,13 @@ .. math:: - f(y) = \\left( - \\begin{matrix} - y_0 ( a-y_0 ) ( y_0-1) - y_1 + k (y_2 - y_0) \\\\ - b_1 y_0 - c y_1 \\\\ - y_2 ( a-y_2 ) ( y_2-1 ) - y_3 + k (y_0 - y_2)\\\\ + f(y) = + \begin{bmatrix} + y_0 ( a-y_0 ) ( y_0-1) - y_1 + k (y_2 - y_0) \\ + b_1 y_0 - c y_1 \\ + y_2 ( a-y_2 ) ( y_2-1 ) - y_3 + k (y_0 - y_2)\\ b_2 y_2 - c y_3 - \\end{matrix} \\right), + \end{bmatrix}, and :math:`a = -0.025794`, :math:`b_1 = 0.0065`, :math:`b_2 = 0.0135`, :math:`c = 0.02`, and :math:`k = 0.128`. Then the following code integrates the above for 100000 time units (after discarding 2000 time units of transients), with :math:`y(t=0) = (1,2,3,4)`, and writes the results to :code:`timeseries.dat`: From e8740ffa724ca6878782c870d90c825475640919 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 15 Jul 2025 10:19:05 +0300 Subject: [PATCH 2/5] Fix line numbers in Lotka-Volterra example --- examples/lotka_volterra.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/examples/lotka_volterra.py b/examples/lotka_volterra.py index 5a1dee7..7b543c2 100644 --- a/examples/lotka_volterra.py +++ b/examples/lotka_volterra.py @@ -16,21 +16,21 @@ .. literalinclude:: ../examples/lotka_volterra.py :start-after: example-st\u0061rt - :lines: 1-2 + :lines: 1-3 :dedent: 1 … and defining the control parameters: .. literalinclude:: ../examples/lotka_volterra.py :start-after: example-st\u0061rt - :lines: 4-7 + :lines: 5-8 :dedent: 1 The `y` that we imported from `jitcode` has to be used for the dynamical variables. However, to make our code use the same notation as the above equation, we can rename them: .. literalinclude:: ../examples/lotka_volterra.py :start-after: example-st\u0061rt - :lines: 9 + :lines: 10 :dedent: 1 We here might as well have written `R,B = y(0),y(1)`, but the above is more handy for larger systems. Note that the following code is written such that the order of our dynamical variables does not matter. @@ -39,7 +39,7 @@ .. literalinclude:: ../examples/lotka_volterra.py :start-after: example-st\u0061rt - :lines: 11-14 + :lines: 12-15 :dedent: 1 Note that there are other ways to define the derivative, e.g., used in the previous and following example. @@ -48,35 +48,35 @@ .. literalinclude:: ../examples/lotka_volterra.py :start-after: example-st\u0061rt - :lines: 16 + :lines: 17 :dedent: 1 Before we start integrating, we have to choose an integrator and define initial conditions. We here choose the 5th-order Dormand–Prince method. Moreover the initial :math:`B` shall be 0.5, the initial :math:`R` shall be 0.2, and the starting time should be :math:`t=0`. .. literalinclude:: ../examples/lotka_volterra.py :start-after: example-st\u0061rt - :lines: 17-19 + :lines: 18-20 :dedent: 1 We then define an array of time points where we want to observe the system, and a dictionary of empty lists that shall be filled with our results. .. literalinclude:: ../examples/lotka_volterra.py :start-after: example-st\u0061rt - :lines: 21-22 + :lines: 22-23 :dedent: 1 Finally, we perform the actual integration by calling `ODE.integrate` for each of our `times`. After this, we use `ODE.y_dict` to access the current state as a dictionary and appending its values to the respective lists. .. literalinclude:: ../examples/lotka_volterra.py :start-after: example-st\u0061rt - :lines: 23-26 + :lines: 24-27 :dedent: 1 We can now plot or analyse our data, but that’s beyond the scope of JiTCODE. So we just save it: .. literalinclude:: ../examples/lotka_volterra.py :start-after: example-st\u0061rt - :lines: 28 + :lines: 29 :dedent: 1 Taking everything together, our code is: @@ -101,9 +101,9 @@ R,B = dynvars = [ y(i) for i in range(2) ] lotka_volterra_diff = { - B: γ*B - φ*R*B, - R: -ω*R + ν*R*B, - } + B: γ*B - φ*R*B, + R: -ω*R + ν*R*B, + } ODE = jitcode(lotka_volterra_diff) ODE.set_integrator("dopri5") From b4104176331138e028bcea8eab434ecc126829dd Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 15 Jul 2025 10:47:06 +0300 Subject: [PATCH 3/5] Fix line numbers in Lyapunov Roesslers example --- docs/index.rst | 2 +- examples/SW_of_Roesslers.py | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/index.rst b/docs/index.rst index f8767f0..ed316e3 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -114,7 +114,7 @@ Estimates for the Lyapunov vectors are returned as well. .. automodule:: double_fhn_lyapunov .. literalinclude:: ../examples/double_fhn_lyapunov.py - :emphasize-lines: 20-21,26,28-36 + :emphasize-lines: 22-23,27,29-37 :start-after: example-start :dedent: 1 :linenos: diff --git a/examples/SW_of_Roesslers.py b/examples/SW_of_Roesslers.py index 75608f6..4e49764 100644 --- a/examples/SW_of_Roesslers.py +++ b/examples/SW_of_Roesslers.py @@ -20,17 +20,17 @@ :linenos: :dedent: 1 :lines: 63- - :emphasize-lines: 9, 25-26, 34, 31, 43 + :emphasize-lines: 10, 26-27, 35, 32, 44 Explanation of selected features and choices: -* The values of :math:`ω` are initialised globally (line 9). We cannot just define a function here, because the parameter is used twice for each oscillator. Moreover, if we were trying to calculate Lyapunov exponents or the Jacobian, the generator function would be called multiple times, and thus the value of the parameter would not be consistent (which would be desastrous). +* The values of :math:`ω` are initialised globally (line 10). We cannot just define a function here, because the parameter is used twice for each oscillator. Moreover, if we were trying to calculate Lyapunov exponents or the Jacobian, the generator function would be called multiple times, and thus the value of the parameter would not be consistent (which would be desastrous). -* Since we need :math:`\\sum_{j=0}^N x_j` to calculate the derivative of :math:`z` for every oscillator, it is prudent to only calculate this once. Therefore we define a helper symbol for this in lines 25 and 26, which we employ in line 34. (See the arguments of `jitcode` for details.) +* Since we need :math:`\\sum_{j=0}^N x_j` to calculate the derivative of :math:`z` for every oscillator, it is prudent to only calculate this once. Therefore we define a helper symbol for this in lines 26 and 27, which we employ in line 35. (See the arguments of `jitcode` for details.) -* In line 31, we implement :math:`\\sin(t)`. For this purpose we had to import `t` in line 1. Also, we need to use `symengine.sin` – in contrast to `math.sin` or `numpy.sin`. +* In line 32, we implement :math:`\\sin(t)`. For this purpose we had to import `t` in line 3. Also, we need to use `symengine.sin` – in contrast to `math.sin` or `numpy.sin`. -* As this is a large system, we use a generator function instead of a list to define :math:`f` (lines 28-35) and have the code automatically be split into chunks of 150 lines, corresponding to the equations of fifty oscillators, in line 43. (For this system, any reasonably sized multiple of 3 is a good choice for `chunk_size`; for other systems, the precise choice of the value may be more relevant.) See `large_systems` for the rationale. +* As this is a large system, we use a generator function instead of a list to define :math:`f` (lines 29-36) and have the code automatically be split into chunks of 150 lines, corresponding to the equations of fifty oscillators, in line 44. (For this system, any reasonably sized multiple of 3 is a good choice for `chunk_size`; for other systems, the precise choice of the value may be more relevant.) See `large_systems` for the rationale. """ if __name__ == "__main__": From c3782890978a568714798eb15c054dd454542dca Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 15 Jul 2025 10:47:20 +0300 Subject: [PATCH 4/5] Fix line numbers in Lyapunov FHN example --- examples/double_fhn_transversal_lyap.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/double_fhn_transversal_lyap.py b/examples/double_fhn_transversal_lyap.py index dedeb43..9ad8140 100644 --- a/examples/double_fhn_transversal_lyap.py +++ b/examples/double_fhn_transversal_lyap.py @@ -4,12 +4,12 @@ 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): .. literalinclude:: ../examples/double_fhn_transversal_lyap.py - :emphasize-lines: 6, 12, 14, 17, 19-20 + :emphasize-lines: 7, 13, 15, 18, 20-22 :start-after: example-st\u0061rt :dedent: 1 :linenos: -Note that the initial state (line 17) is reduced in dimensionality and there is only one component for each synchronisation group. +Note that the initial state (line 18) is reduced in dimensionality and there is only one component for each synchronisation group. """ if __name__ == "__main__": From c13fd2851ae7f3634d1791b9b21399476272f5e0 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 15 Jul 2025 11:06:00 +0300 Subject: [PATCH 5/5] Fix links and references in docs --- docs/conf.py | 7 +++++++ docs/index.rst | 16 +++++++++------- examples/SW_of_Roesslers.py | 4 ++-- 3 files changed, 18 insertions(+), 9 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 58cfd96..2131f0f 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -28,6 +28,7 @@ class GroupHandler_mock: needs_sphinx = "1.3" extensions = [ + "sphinx.ext.intersphinx", "sphinx.ext.autodoc", "sphinx.ext.autosummary", "sphinx.ext.mathjax", @@ -61,6 +62,12 @@ class GroupHandler_mock: toc_object_entries_show_parents = "hide" +intersphinx_mapping = { + "numpy": ("https://numpy.org/doc/stable", None), + "python": ("https://docs.python.org/3", None), + "scipy": ("https://docs.scipy.org/doc/scipy", None), +} + def on_missing_reference(app, env, node, contnode): if node["reftype"] == "any": return contnode diff --git a/docs/index.rst b/docs/index.rst index ed316e3..62714ad 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -7,7 +7,7 @@ The JiTCODE module Overview -------- -JiTCODE (just-in-time compilation for ordinary differential equations) is an extension of SciPy’s `ODE`_ (`scipy.integrate.ode`) or `Solve IVP`_ (`scipy.integrate.solve_ivp`). +JiTCODE (just-in-time compilation for ordinary differential equations) is an extension of SciPy’s `ODE`_ (:class:`scipy.integrate.ode`) or `Solve IVP`_ (:func:`scipy.integrate.solve_ivp`). Where the latter take a Python function as an argument, JiTCODE takes an iterable (or generator function or dictionary) of symbolic expressions, which it translates to C code, compiles on the fly, and uses as the function to feed into SciPy’s ODE or Solve IVP. Symbolic expressions are mostly handled by `SymEngine`_, `SymPy`_’s compiled-backend-to-be (see `SymPy vs. SymEngine`_ for details). @@ -69,9 +69,11 @@ An explicit example Details of the building process ------------------------------- +.. currentmodule:: jitcode + To generate the functions needed by SciPy’s ODE or Solve IVP, JiTCODE executes a series of distinct processing steps, each of which can be invoked through a command and tweaked as needed. These commands execute previous steps as needed, i.e., if they have not been performed yet. -If you are happy with the default options, however, you do not need to bother with this and can just use the commands at the very end of the chain, namely `set_integrator`, `set_initial_value`, or `save_compiled`. +If you are happy with the default options, however, you do not need to bother with this and can just use the commands at the very end of the chain, namely :meth:`~_jitcode.jitcode.set_integrator`, :meth:`~_jitcode.jitcode.set_initial_value`, or :meth:`~_jitcode.jitcode.save_compiled`. The following diagram details which command calls which other command when needed: @@ -102,12 +104,12 @@ A more complicated example Calculating Lyapunov exponents with `jitcode_lyap` -------------------------------------------------- -`jitcode_lyap` is a simple extension of `jitcode` that almost automatically handles calculating Lyapunov exponents by evolving tangent vectors [BGGS80]_. -It works just like `jitcode`, except that it generates and integrates additional differential equations for the tangent vectors. -After every call of `integrate`, the tangent vectors are orthonormalised, and the “local” Lyapunov exponents for this integration step are returned alongside with the system’s state. +:class:`~_jitcode.jitcode_lyap` is a simple extension of :class:`~_jitcode.jitcode` that almost automatically handles calculating Lyapunov exponents by evolving tangent vectors [BGGS80]_. +It works just like :class:`~_jitcode.jitcode`, except that it generates and integrates additional differential equations for the tangent vectors. +After every call of :meth:`~_jitcode.jitcode_lyap.integrate`, the tangent vectors are orthonormalised, and the “local” Lyapunov exponents for this integration step are returned alongside with the system’s state. These can then be further processed to obtain the Lyapunov exponents. The tangent vectors are initialised with random vectors, and you have to take care of the preiterations that the tangent vectors require to align themselves. -You also have to take care not to `integrate` for too long to avoid a numerical overflow in the tangent vectors. +You also have to take care not to :meth:`~_jitcode.jitcode_lyap.integrate` for too long to avoid a numerical overflow in the tangent vectors. Estimates for the Lyapunov vectors are returned as well. @@ -124,7 +126,7 @@ Estimates for the Lyapunov vectors are returned as well. Calculating transversal Lyapunov exponents with `jitcode_transversal_lyap` -------------------------------------------------------------------------- -`jitcode_transversal_lyap` is a variant of `jitcode_lyap` that calculates Lyapunov exponents transversal to a user-defined synchronisation manifold. +:class:`~_jitcode.jitcode_transversal_lyap` is a variant of :class:`~_jitcode.jitcode_lyap` that calculates Lyapunov exponents transversal to a user-defined synchronisation manifold. It automatically conflates the differential equations of a group of synchronised components into one equation on the synchronisation manifold. Moreover, it transforms the equations for the tangent vectors such that the tangent vector is automatically orthogonal to the synchronisation manifold. If you are interested in the mathematical details, please read `the accompanying paper`_. diff --git a/examples/SW_of_Roesslers.py b/examples/SW_of_Roesslers.py index 4e49764..699a4ae 100644 --- a/examples/SW_of_Roesslers.py +++ b/examples/SW_of_Roesslers.py @@ -12,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` 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)`. +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 :any:`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: @@ -28,7 +28,7 @@ * Since we need :math:`\\sum_{j=0}^N x_j` to calculate the derivative of :math:`z` for every oscillator, it is prudent to only calculate this once. Therefore we define a helper symbol for this in lines 26 and 27, which we employ in line 35. (See the arguments of `jitcode` for details.) -* In line 32, we implement :math:`\\sin(t)`. For this purpose we had to import `t` in line 3. Also, we need to use `symengine.sin` – in contrast to `math.sin` or `numpy.sin`. +* In line 32, we implement :math:`\\sin(t)`. For this purpose we had to import `t` in line 3. Also, we need to use `symengine.sin` – in contrast to :func:`math.sin` or :data:`numpy.sin`. * As this is a large system, we use a generator function instead of a list to define :math:`f` (lines 29-36) and have the code automatically be split into chunks of 150 lines, corresponding to the equations of fifty oscillators, in line 44. (For this system, any reasonably sized multiple of 3 is a good choice for `chunk_size`; for other systems, the precise choice of the value may be more relevant.) See `large_systems` for the rationale. """