Added new lyapunov spectrum measure with tests.

README.md

@@ -40,7 +40,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...

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)

tests/test_measures.py

@@ -1,19 +1,35 @@
 """Testing simulations."""
+
 import numpy as np
+import pytest
 
-from lorenzpy import measures
+from lorenzpy import measures, simulations
 
 
-def test_largest_lyapunov_one_d_linear_time_dependent():
-    """Testing measures.largest_lyapunov_exponent."""
+@pytest.fixture
+def get_demo_exponential_decay_system():
+    """Fixture to get a demo iterator func and starting point.
+
+    Returns the iterator_func, starting_point, dt and a for exponential decay.
+    """
     dt = 0.1
     a = -0.1
 
     def iterator_func(x):
         return x + dt * a * x
 
-    return_convergence = False
     starting_point = np.array(1)
+
+    return iterator_func, starting_point, dt, a
+
+
+def test_largest_lyapunov_one_d_linear_time_dependent(
+    get_demo_exponential_decay_system,
+):
+    """Testing measures.largest_lyapunov_exponent."""
+    iterator_func, starting_point, dt, a = get_demo_exponential_decay_system
+
+    return_convergence = False
     deviation_scale = 1e-10
     steps = 10
     steps_skip = 1
@@ -33,16 +49,13 @@ def iterator_func(x):
     np.testing.assert_almost_equal(actual, desired, 2)
 
 
-def test_largest_lyapunov_one_d_linear_time_dependent_return_conv():
+def test_largest_lyapunov_one_d_linear_time_dependent_return_conv(
+    get_demo_exponential_decay_system,
+):
     """Testing measures.largest_lyapunov_exponent with convergence."""
-    dt = 0.1
-    a = -0.1
-
-    def iterator_func(x):
-        return x + dt * a * x
+    iterator_func, starting_point, dt, a = get_demo_exponential_decay_system
 
     return_convergence = True
-    starting_point = np.array(1)
     deviation_scale = 1e-10
     steps = 10
     steps_skip = 1
@@ -65,3 +78,93 @@ def iterator_func(x):
         * steps
     )
     np.testing.assert_almost_equal(actual, desired, 2)
+
+
+def test_lyapunov_spectrum_lorenz63():
+    """Testing lyapunov_exponent_spectrum with initial_pert_directions."""
+    dt = 0.05
+    lorenz63_obj = simulations.Lorenz63(dt=dt)
+    iterator_func = lorenz63_obj.iterate
+    starting_point = lorenz63_obj.get_default_starting_pnt()
+
+    m = 3
+    deviation_scale = 1e-10
+    steps = 1000
+    part_time_steps = 15
+    steps_skip = 50
+    initial_pert_directions = None
+    return_convergence = False
+
+    lyap_spec_conv = measures.lyapunov_exponent_spectrum(
+        iterator_func=iterator_func,
+        starting_point=starting_point,
+        deviation_scale=deviation_scale,
+        steps=steps,
+        part_time_steps=part_time_steps,
+        steps_skip=steps_skip,
+        dt=dt,
+        m=m,
+        initial_pert_directions=initial_pert_directions,
+        return_convergence=return_convergence,
+    )
+
+    expected_lyapunovs = np.array([9.05e-01, 1.70e-03, -1.44e01])
+    np.testing.assert_almost_equal(lyap_spec_conv, expected_lyapunovs, decimal=1)
+
+
+def test_lyapunov_spectrum_vs_largest_le_lorenz63():
+    """Testing lyapunov_exponent_spectrum with initial_pert_directions."""
+    dt = 0.05
+    lorenz63_obj = simulations.Lorenz63(dt=dt)
+    iterator_func = lorenz63_obj.iterate
+    starting_point = lorenz63_obj.get_default_starting_pnt()
+
+    m = 1
+    deviation_scale = 1e-10
+    steps = 500
+    part_time_steps = 10
+    steps_skip = 50
+    return_convergence = True
+
+    lyap_spec_conv = measures.lyapunov_exponent_spectrum(
+        iterator_func=iterator_func,
+        starting_point=starting_point,
+        deviation_scale=deviation_scale,
+        steps=steps,
+        part_time_steps=part_time_steps,
+        steps_skip=steps_skip,
+        dt=dt,
+        m=m,
+        initial_pert_directions=None,
+        return_convergence=return_convergence,
+    )
+
+    lyap_largest_conv = measures.largest_lyapunov_exponent(
+        iterator_func=iterator_func,
+        starting_point=starting_point,
+        deviation_scale=deviation_scale,
+        steps=steps,
+        part_time_steps=part_time_steps,
+        steps_skip=steps_skip,
+        dt=dt,
+        initial_pert_direction=None,
+        return_convergence=return_convergence,
+    )
+
+    np.testing.assert_almost_equal(lyap_spec_conv, lyap_largest_conv[:, np.newaxis], 3)
+
+
+def test_lyapunov_spectrum_custom_pert(get_demo_exponential_decay_system):
+    """Test if it works if custom initial_pert_directions is given."""
+    iterator_func, starting_point, dt, a = get_demo_exponential_decay_system
+
+    lle = measures.lyapunov_exponent_spectrum(
+        iterator_func=iterator_func,
+        starting_point=starting_point,
+        dt=dt,
+        steps=10,
+        part_time_steps=1,
+        steps_skip=0,
+        initial_pert_directions=np.array([[1]]),
+    )
+    np.testing.assert_almost_equal(lle, a, decimal=2)