Merge pull request #4 from DuncDennis/feature/get_ready_for_v002

Feature/get ready for v002

Author: DuncDennis

CHANGELOG.md

@@ -2,6 +2,15 @@
 
 <!--next-version-placeholder-->
 
+### v0.0.2 (21.10.2023)
+- Added more Systems.
+- Added new Lyapunov spectrum measure.
+- Added new simulation backend.
+- Added example folder.
+- Added static folder to generate animation in readme.
+- Changed plotly dependency to matplotlib.
+- More smaller improvements...
+
 ### v0.0.1 (16.03.2023)
 - First release of `lorenzpy`
 

README.md

@@ -22,8 +22,7 @@ $ pip install lorenzpy
 ```
 
 To install with the additional plotting functionality.
-This also installs `plotly` and `matplotlib`. ⚠️ Plotting functionality not in a useful
-state.
+This also installs `matplotlib`. ⚠️ Plotting functionality not in a useful state.
 ```bash
 $ pip install lorenzpy[plot]
 ```
@@ -59,6 +58,8 @@ 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].
 
+For more examples see the [examples folder](examples/README.md).
+
 ## 💫 Supported systems
 
 

docs/reference.md

@@ -3,6 +3,7 @@
 
 ### simulations module:
 ::: lorenzpy.simulations
+::: lorenzpy.simulations.solvers
 
 ### measures module:
 ::: lorenzpy.measures

examples/README.md

@@ -0,0 +1,8 @@
+# Examples:
+
+### Contents
+- ``double_pendulum_lyapunov_spectrum.py``: Plot the lyapunov spectrum of the double pendulum.
+- ``lyapunov_spectra_of_3d_autonomous_flows.py``: Plot the Lyapunov spectrum of all 3D autonomous dissipative flow systems.
+
+⚠️ Not many examples here yet, and examples might be flawed.
+

examples/double_pendulum_lyapunov_spectrum.py

@@ -0,0 +1,68 @@
+"""Lyapunov Spectrum of single system."""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from lorenzpy import measures as meas
+from lorenzpy import simulations as sims
+
+sys_obj = sims.DoublePendulum(dt=0.1)
+dt = sys_obj.dt
+
+# Calculate exponents:
+m = 4
+deviation_scale = 1e-10
+steps = 1000
+part_time_steps = 15
+steps_skip = 0
+
+iterator_func = sys_obj.iterate
+starting_point = sys_obj.get_default_starting_pnt()
+
+le_spectrum = meas.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=True,
+)
+
+fig, ax = plt.subplots(
+    1, 1, figsize=(6, 6), layout="constrained", sharex=True, sharey=False
+)
+
+# x and y titles:
+fig.supxlabel("Number of renormalization steps")
+fig.supylabel("Lyapunov exponent convergence")
+
+x = np.arange(1, steps + 1)
+ax.plot(
+    x,
+    le_spectrum,
+    linewidth=1,
+)
+ax.grid(True)
+
+final_les = np.round(le_spectrum[-1, :], 4).tolist()
+final_les = [str(x) for x in final_les]
+le_string = "\n".join(final_les)
+le_string = "Final LEs: \n" + le_string
+x_position = 0.1  # X-coordinate of the upper-left corner for each subplot
+y_position = 0.5
+ax.text(
+    x_position,
+    y_position,
+    le_string,
+    fontsize=10,
+    bbox=dict(facecolor="white", edgecolor="black", boxstyle="round"),
+    verticalalignment="center",
+    horizontalalignment="left",
+    transform=ax.transAxes,
+)
+
+plt.show()

examples/lyapunov_spectra_of_3d_autonomous_flows.py

@@ -0,0 +1,91 @@
+"""Calculate and plot the Lyapunov Spectra of all three-dimensional chaotic flows."""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from lorenzpy import measures as meas
+from lorenzpy import simulations as sims
+
+systems = [
+    "Lorenz63",
+    "Chen",
+    "ChuaCircuit",
+    "ComplexButterfly",
+    "DoubleScroll",
+    "Halvorsen",
+    "Roessler",
+    "Rucklidge",
+    "Thomas",
+    "WindmiAttractor",
+]
+
+# Calculate exponents:
+m = 3
+deviation_scale = 1e-10
+steps = 1000
+part_time_steps = 15
+steps_skip = 50
+
+solver = "rk4"
+# solver = sims.solvers.create_scipy_ivp_solver(method="RK45")
+
+lyap_dict = {}
+for i_sys, system in enumerate(systems):
+    print(system)
+    sys_obj = getattr(sims, system)(solver=solver)
+    iterator_func = sys_obj.iterate
+    starting_point = sys_obj.get_default_starting_pnt()
+    dt = sys_obj.dt
+
+    lyap_dict[system] = meas.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=True,
+    )
+
+fig, axs = plt.subplots(
+    2, 5, figsize=(15, 8), layout="constrained", sharex=True, sharey=False
+)
+
+# x and y titles:
+fig.supxlabel("Number of renormalization steps")
+fig.supylabel("Lyapunov exponent convergence")
+
+axs = axs.flatten()
+x = np.arange(1, steps + 1)
+for i_ax, ax in enumerate(axs):
+    system = systems[i_ax]
+    le_spectrum = lyap_dict[system]
+    ax.title.set_text(system)
+    ax.plot(
+        x,
+        le_spectrum,
+        linewidth=1,
+    )
+    ax.grid(True)
+
+    final_les = np.round(le_spectrum[-1, :], 4).tolist()
+    final_les = [str(x) for x in final_les]
+    le_string = "\n".join(final_les)
+    le_string = "Final LEs: \n" + le_string
+    x_position = 0.1  # X-coordinate of the upper-left corner for each subplot
+    y_position = 0.5
+    ax.text(
+        x_position,
+        y_position,
+        le_string,
+        fontsize=10,
+        bbox=dict(facecolor="white", edgecolor="black", boxstyle="round"),
+        verticalalignment="center",
+        horizontalalignment="left",
+        transform=ax.transAxes,
+    )
+
+plt.show()

pyproject.toml

@@ -1,7 +1,7 @@
 [project]
 name = "lorenzpy"
 readme = "README.md"
-version = "0.0.1"
+version = "0.0.2"
 description = "A Python package to simulate and measure chaotic dynamical systems."
 authors = [
   {name = "Dennis Duncan", email = "[email protected]"},
@@ -41,7 +41,6 @@ dev = [
     "pre-commit==3.1.1",  # add version?
 ]
 plot = [
-    "plotly>=5.11",
     "matplotlib>=3.5"
 ]
 
@@ -59,7 +58,7 @@ line-length = 88
 files = "src/lorenzpy/"
 
 [[tool.mypy.overrides]]
-module = ['plotly.*', 'numpy', 'pytest', "scipy.*", "matplotlib.*", "PIL"]   # ignore missing imports from the plotly package.
+module = ['numpy', 'pytest', "scipy.*", "matplotlib.*", "PIL", "mpl_toolkits.*"]
 ignore_missing_imports = true
 
 [tool.ruff]

src/lorenzpy/__init__.py

@@ -8,4 +8,4 @@
 
 from . import measures, simulations
 
-__version__ = "0.0.1"
+__version__ = "0.0.2"

src/lorenzpy/plot/__init__.py

@@ -1 +1,2 @@
 """Submodule to plot chaotic dynamics systems."""
+from .plot import create_3d_line_plot

src/lorenzpy/plot/plot.py

@@ -0,0 +1,39 @@
+"""Plot the data of dynamical systems."""
+
+import matplotlib.pyplot as plt
+import numpy as np
+from mpl_toolkits.mplot3d import Axes3D
+
+
+def create_3d_line_plot(data: np.ndarray, ax: "Axes3D" = None, **kwargs) -> "Axes3D":
+    """Create a three-dimensional line plot of data.
+
+    Args:
+        data (np.ndarray): A NumPy array containing 3D data points with shape (n, 3).
+        ax (Axes3D, optional): A Matplotlib 3D axis to use for plotting.
+        If not provided, a new 3D plot will be created.
+        **kwargs: Additional keyword arguments to pass to ax.plot.
+
+    Returns:
+        Axes3D: The Matplotlib 3D axis used for the plot.
+
+    Example:
+        >>> data = np.random.rand(100, 3)  # Replace this with your own 3D data
+        >>> create_3d_line_plot(data, color='b', linestyle='--')
+        >>> plt.show()
+    """
+    if ax is None:
+        fig = plt.figure()
+        ax = fig.add_subplot(111, projection="3d")
+
+    x = data[:, 0]
+    y = data[:, 1]
+    z = data[:, 2]
+
+    ax.plot(x, y, z, **kwargs)  # Use plot for a line plot with kwargs
+
+    ax.set_xlabel("X Axis")
+    ax.set_ylabel("Y Axis")
+    ax.set_zlabel("Z Axis")
+
+    return ax

src/lorenzpy/plot/plot_3d.py

@@ -28,24 +28,8 @@
     >>> data = sims.Lorenz63().simulate(1000)
     >>> data.shape
     (1000, 3)
-
-
-TODO <below>
-    - Probably for each concrete simulation class + public methods. Compare with sklearn
-    - Find out which functionality is missing. E.g. Raising error when wrong values are
-      parsed.
-    - Check where to add proper tests and how to add them efficiently. Fixtures?
-      Parametrization?
-    - Implement all the other dynamical systems.
-    - Check if the names of files and functions make sense?
-    - Add functionality to add your own dynamical system? As my base-classes are
-      protected this is maybe not so easy? -> Make ABC public?
-    - Think about adding NARMA? Maybe I need a random number generator framework.
-    - Check if I can further reduce code duplication. Maybe regarding solvers.
-    - Check for proper doc-generation. It seems that the methods of inhereted members
-      is not implemented yet. See:
-      https://github.com/mkdocstrings/mkdocstrings/issues/78
 """
+
 from . import solvers
 from .autonomous_flows import (
     Chen,

src/lorenzpy/simulations/__init__.py

@@ -95,7 +95,7 @@ class ComplexButterfly(_BaseSimFlow):
     def __init__(
         self,
         a: float = 0.55,
-        dt: float = 0.05,
+        dt: float = 0.1,
         solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
     ):
         """Initialize the ComplexButterfly simulation object.

src/lorenzpy/simulations/autonomous_flows.py

@@ -228,13 +228,12 @@ def simulate(
         if starting_point.size == self.history_steps + 1:
             initial_history = starting_point
         elif starting_point.size == 1:
-            initial_history = np.repeat(starting_point, self.history_steps)
+            initial_history = np.repeat(starting_point, self.history_steps + 1)
         else:
             raise ValueError("Wrong size of starting point.")
 
         traj_w_hist = np.zeros((self.history_steps + time_steps, 1))
-        traj_w_hist[: self.history_steps, :] = initial_history[:, np.newaxis]
-        traj_w_hist[self.history_steps, :] = starting_point
+        traj_w_hist[: self.history_steps + 1, :] = initial_history[:, np.newaxis]
 
         for t in range(1, time_steps + transient):
             t_shifted = t + self.history_steps