Merge pull request #3 from DuncDennis/feature/try_new_simulation_backend

New simulation backend, new systems, Lyapunov spectrum measure.

.github/workflows/test_and_coverage.yml

@@ -7,7 +7,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: ["3.8", "3.9", "3.10", "3.11"]
+        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
     name: Python ${{ matrix.python-version }} on ${{ matrix.os }}
     steps:
       - uses: actions/checkout@v3

.pre-commit-config.yaml

@@ -4,13 +4,14 @@ repos:
     hooks:
     -   id: trailing-whitespace
     -   id: check-added-large-files
+        args: ['--maxkb=1500']
 -   repo: https://github.com/psf/black
     rev: "23.1.0"
     hooks:
     -   id: black
         language_version: python3
 -   repo: https://github.com/charliermarsh/ruff-pre-commit
-    rev: 'v0.0.245'
+    rev: 'v0.1.0'
     hooks:
     -   id: ruff
 -   repo: https://github.com/pre-commit/mirrors-mypy

README.md

@@ -1,11 +1,19 @@
 # LorenzPy
 
+A Python package to simulate and measure chaotic dynamical systems.
+
 [![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black)
 [![Ruff](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/charliermarsh/ruff/main/assets/badge/v1.json)](https://github.com/charliermarsh/ruff)
 [![codecov](https://codecov.io/gh/DuncDennis/lorenzpy/branch/main/graph/badge.svg?token=ATWAEQHBYB)](https://codecov.io/gh/DuncDennis/lorenzpy)
 [![license: MIT](https://img.shields.io/badge/License-MIT-purple.svg)](LICENSE)
 [![Python versions](https://img.shields.io/badge/python-3.8+-blue.svg)](https://www.python.org/downloads/)
 
+------
+
+![Flow-Attractors](static/attractor_animation.gif)
+
+------
+
 ## ⚙️ Installation
 
 To install only the core functionality:
@@ -14,11 +22,13 @@ $ pip install lorenzpy
 ```
 
 To install with the additional plotting functionality.
-This also installs `plotly`.
+This also installs `plotly` and `matplotlib`. ⚠️ Plotting functionality not in a useful
+state.
 ```bash
 $ pip install lorenzpy[plot]
 ```
 
+
 ## ▶️ Usage
 
 LorenzPy can be used to simulate and measure chaotic dynamical systems.
@@ -40,7 +50,7 @@ data = l63_obj.simulate(5000)    # -> data.shape = (5000,3)
 iterator = l63_obj.iterate
 lle = lpy.measures.largest_lyapunov_exponent(
     iterator_func=iterator,
-    starting_point=l63_obj.default_starting_point,
+    starting_point=l63_obj.get_default_starting_pnt(),
     dt=l63_obj.dt
 )
 # -> lle = 0.905144329...
@@ -49,11 +59,33 @@ lle = lpy.measures.largest_lyapunov_exponent(
 The calculated largest Lyapunov exponent of *0.9051...* is very close to the literature
 value of *0.9056*[^SprottChaos].
 
+## 💫 Supported systems
+
+
+| Name                                  | Type                        | System Dimension |
+|:--------------------------------------|-----------------------------|:-----------------|
+| `Lorenz63`                            | autonomous dissipative flow | 3                |
+| `Roessler`                            | autonomous dissipative flow | 3                |
+| `ComplexButterfly`                    | autonomous dissipative flow | 3                |
+| `Chen`                                | autonomous dissipative flow | 3                |
+| `ChuaCircuit`                         | autonomous dissipative flow | 3                |
+| `Thomas`                              | autonomous dissipative flow | 3                |
+| `WindmiAttractor`                     | autonomous dissipative flow | 3                |
+| `Rucklidge`                     | autonomous dissipative flow | 3                |
+| `Halvorsen`                     | autonomous dissipative flow | 3                |
+| `DoubleScroll`                     | autonomous dissipative flow | 3                |
+| `Lorenz96`                            | autonomous dissipative flow | variable         |
+| `DoublePendulum`                      | conservative flow           | 4                |
+| `Logistic`                            | noninvertible map           | 1                |
+| `Henon`                               | dissipative map             | 2                |
+| `SimplestDrivenChaoticFlow`           | conservative flow           | 2 space + 1 time |
+| `KuramotoSivashinsky`                 | PDE                         | variable         |
+| `MackeyGlass`                         | delay differential equation | variable         |
 ## 📗 Documentation
 
 - The main documentation can be found here: https://duncdennis.github.io/lorenzpy/
     - ⚠️: The documentation is not in a useful state.
-##  ⚠️ Further notes:
+##  ⚠️ Further notes
 - So far the usefulness of this package is very limited.
 The authors main purpose to creating this package was to learn the full workflow to
 develop a Python package.

pyproject.toml

@@ -18,6 +18,7 @@ classifiers = [
 ]
 dependencies = [
     "numpy>=1.22.0",
+    "scipy>=1.10.0",
 ]
 
 [project.urls]
@@ -33,14 +34,15 @@ dev = [
     "pytest-cov==4.0",
     "black==23.1.0",
     "mypy==1.1.1",
-    "ruff==0.0.254",
+    "ruff==0.1.0",
     "mkdocs",  # add version?
     "mkdocstrings[python]",  # add version?
     "mkdocs-material",  # add version?
     "pre-commit==3.1.1",  # add version?
 ]
 plot = [
     "plotly>=5.11",
+    "matplotlib>=3.5"
 ]
 
 [tool.pytest.ini_options]
@@ -57,7 +59,7 @@ line-length = 88
 files = "src/lorenzpy/"
 
 [[tool.mypy.overrides]]
-module = ['plotly.*', 'numpy', 'pytest']   # ignore missing imports from the plotly package.
+module = ['plotly.*', 'numpy', 'pytest', "scipy.*", "matplotlib.*", "PIL"]   # ignore missing imports from the plotly package.
 ignore_missing_imports = true
 
 [tool.ruff]

src/lorenzpy/measures.py

@@ -74,3 +74,95 @@ def largest_lyapunov_exponent(
         )
     else:
         return float(np.average(log_divergence) / (dt * part_time_steps))
+
+
+def lyapunov_exponent_spectrum(
+    iterator_func: Callable[[np.ndarray], np.ndarray],
+    starting_point: np.ndarray,
+    deviation_scale: float = 1e-10,
+    steps: int = int(1e3),
+    part_time_steps: int = 15,
+    steps_skip: int = 50,
+    dt: float = 1.0,
+    m: int | None = None,
+    initial_pert_directions: np.ndarray | None = None,
+    return_convergence: bool = False,
+) -> np.ndarray:
+    """Calculate the spectrum of m largest lyapunov exponent given an iterator function.
+
+    A mixture of:
+    - The algorithm for the largest lyapunov exponent: Sprott, Julien Clinton, and
+        Julien C. Sprott. Chaos and time-series analysis. Vol. 69. Oxford: Oxford
+        university press, 2003.
+    - The algorithm for the spectrum given in 1902.09651 "LYAPUNOV EXPONENTS of the
+        KURAMOTO-SIVASHINSKY PDE".
+
+    Args:
+        iterator_func: Function to iterate the system to the next time step:
+                        x(i+1) = F(x(i))
+        starting_point: The starting_point of the main trajectory.
+        deviation_scale: The L2-norm of the initial perturbation.
+        steps: Number of renormalization steps.
+        part_time_steps: Time steps between renormalization steps.
+        steps_skip: Number of renormalization steps to perform, before tracking the log
+                    divergence.
+                Avoid transients by using steps_skip.
+        dt: Size of time step.
+        m: Number of Lyapunov exponents to compute. If None: take all (m = x_dim).
+        initial_pert_directions:
+            - If np.ndarray: The direction of the initial perturbations of shape
+                            (m, x_dim). Will be qr decomposed
+                             and the q matrix will be used.
+            - If None: The direction of the initial perturbation is assumed to be
+                        np.eye(x_dim)[:m, :].
+        return_convergence: If True, return the convergence of the largest LE; a
+                            numpy array of the shape (N, m).
+
+    Returns:
+        The Lyapunov exponent spectrum of largest m values. If return_convergence is
+        True: The convergence (2D N x m np.ndarray), else a (1D m-size np.ndarray),
+        which holds the last values in the convergence.
+    """
+    x_dim = starting_point.size
+
+    if m is None:
+        m = x_dim
+
+    if initial_pert_directions is None:
+        initial_perturbations = (
+            np.eye(x_dim)[:m, :] * deviation_scale
+        )  # each vector is of length deviation_scale
+    else:
+        q, _ = np.linalg.qr(initial_pert_directions)
+        initial_perturbations = q * deviation_scale
+
+    log_divergence_spec = np.zeros((steps, m))
+
+    x = starting_point
+    x_pert = starting_point + initial_perturbations
+
+    for i_n in range(steps + steps_skip):
+        for i_t in range(part_time_steps):
+            x = iterator_func(x)
+            for i_m in range(m):
+                x_pert[i_m, :] = iterator_func(x_pert[i_m, :])
+        dx = x_pert - x
+        q, r = np.linalg.qr(dx.T, mode="reduced")
+        x_pert = x + q.T * deviation_scale
+
+        if i_n >= steps_skip:
+            log_divergence_spec[i_n - steps_skip, :] = np.log(
+                np.abs(np.diag(r)) / deviation_scale
+            )
+
+    if return_convergence:
+        return np.array(
+            np.cumsum(log_divergence_spec, axis=0)
+            / (
+                np.repeat(np.arange(1, steps + 1)[:, np.newaxis], m, axis=1)
+                * dt
+                * part_time_steps
+            )
+        )
+    else:
+        return np.average(log_divergence_spec, axis=0) / (dt * part_time_steps)