Python API Reference

Complete reference for the koopman_dmd Python module (PyO3 bindings)

Installation

The Python bindings are built with PyO3 and packaged using maturin. To build and install from source:

# Install the build tool
pip install maturin

# Navigate to the Python binding directory
cd koopman-dmd-py

# Build and install in the current environment
maturin develop --release

After installation, the module is available as koopman_dmd:

import koopman_dmd

DMD

The primary class for Dynamic Mode Decomposition. Supports standard DMD, centering, and Extended DMD via lifting functions.

Constructor

class DMD(rank=None, center=False, lifting=None, lifting_param=None)

from koopman_dmd import DMD

# Auto-rank, no centering
dmd = DMD()

# Explicit rank with centering
dmd = DMD(rank=4, center=True)

# Extended DMD with polynomial lifting (degree 3)
dmd = DMD(lifting="polynomial", lifting_param=3)

# Extended DMD with trigonometric lifting (2 harmonics)
dmd = DMD(lifting="trigonometric", lifting_param=2)

# Extended DMD with delay embedding (5 delays)
dmd = DMD(lifting="delay", lifting_param=5)

Methods

fit(data)

Fit the DMD model to a data matrix.

Returns: self -- the fitted DMD object, enabling method chaining.

import numpy as np
from koopman_dmd import DMD

data = np.random.randn(3, 100)  # 3 variables, 100 time steps
dmd = DMD(rank=2).fit(data)

predict(n_ahead, x0=None)

Predict future states using the spectral decomposition. Each time step is computed independently from eigenvalues, modes, and amplitudes (non-recursive). Suitable for long-range forecasting.

Returns: numpy.ndarray of shape (n_vars, n_ahead).

# Predict 50 steps ahead from the training initial condition
forecast = dmd.predict(50)

# Predict from a custom initial condition
x0 = np.array([1.0, 0.0, -0.5])
forecast = dmd.predict(50, x0=x0)

predict_matrix(n_ahead, x0=None)

Predict future states using recursive matrix multiplication. Each step is computed by applying the reconstructed linear operator to the previous state. Faster for short horizons.

Returns: numpy.ndarray of shape (n_vars, n_ahead).

forecast_mat = dmd.predict_matrix(50)

reconstruct(n_time, x0=None)

Reconstruct the fitted signal over a given number of time steps. Uses the spectral decomposition to reproduce the training dynamics. If centering was enabled, the stored mean is added back automatically.

Returns: numpy.ndarray of shape (n_vars, n_time).

fitted = dmd.reconstruct(100)

spectrum(dt)

Compute the continuous-time spectral properties of the decomposition. Converts discrete-time eigenvalues to frequencies, growth rates, and stability indicators.

Returns: list[dict] -- one dictionary per mode with keys:

spec = dmd.spectrum(dt=0.01)
for mode in spec:
    print(f"freq={mode['frequency']:.4f} Hz, growth={mode['growth_rate']:.4f}, {mode['stability']}")

stability(tol=1e-6)

Assess the overall stability of the fitted linear system by examining the eigenvalue magnitudes.

Returns: dict with keys:

stab = dmd.stability()
print(f"Spectral radius: {stab['spectral_radius']:.6f}")
print(f"Stable: {stab['is_stable']}")
print(f"Growing/Decaying/Neutral: {stab['n_growing']}/{stab['n_decaying']}/{stab['n_neutral']}")

error(data)

Compute reconstruction error metrics by comparing the DMD reconstruction to the original data.

Returns: dict with keys:

err = dmd.error(data)
print(f"RMSE:     {err['rmse']:.6f}")
print(f"Relative: {err['relative_error']:.6f}")
print(f"Max:      {err['max_error']:.6f}")

dominant_modes(criterion, n)

Extract the most significant modes according to a ranking criterion.

Returns: list[dict] -- one dictionary per mode, sorted by the chosen criterion in descending order.

top = dmd.dominant_modes(criterion="energy", n=3)
for m in top:
    print(m)

residual()

Compute per-mode residual norms, measuring how well each individual mode satisfies the DMD eigenvalue equation.

Returns: dict with keys:

res = dmd.residual()
print(f"Max residual:  {res['max_residual']:.2e}")
print(f"Mean residual: {res['mean_residual']:.2e}")

Properties

print("Eigenvalues:", dmd.eigenvalues)
print("Modes shape:", dmd.modes.shape)
print("Amplitudes: ", dmd.amplitudes)
print("Rank:       ", dmd.rank)

HankelDMD

DMD with time-delay embedding via Hankel matrices. Captures higher-order dynamics from scalar or low-dimensional measurements by augmenting the state space with lagged copies of the data.

Constructor

class HankelDMD(delays=None, rank=None, dt=1.0)

from koopman_dmd import HankelDMD

hdmd = HankelDMD(delays=10, rank=4, dt=0.01)

Methods

fit(data)

Fit the Hankel-DMD model. Internally constructs the Hankel matrix from the input data, then performs standard DMD on the augmented state.

Returns: self

predict(n_ahead)

Predict future states. Output is projected back to the original observation dimension.

Returns: numpy.ndarray of shape (n_vars, n_ahead).

reconstruct()

Reconstruct the training data from the Hankel-DMD model.

Returns: numpy.ndarray of shape (n_vars, n_time).

Properties

import numpy as np
from koopman_dmd import HankelDMD

# Scalar time series with delay embedding
t = np.linspace(0, 10, 500)
signal = np.sin(2.0 * np.pi * t) + 0.5 * np.sin(6.0 * np.pi * t)
data = signal.reshape(1, -1)

hdmd = HankelDMD(delays=20, rank=4).fit(data)

print("Eigenvalues:", hdmd.eigenvalues)
print("Rank:", hdmd.rank)

recon = hdmd.reconstruct()
forecast = hdmd.predict(100)

GLA

Generalized Laplace Analysis for extracting Koopman eigenvalues and eigenfunctions from trajectory data using spectral methods. GLA solves an optimization problem to find the eigenvalues that best explain the observed dynamics.

Constructor

class GLA(n_eigenvalues=2, tol=1e-6, max_iter=None)

from koopman_dmd import GLA

gla = GLA(n_eigenvalues=4, tol=1e-8)

Methods

fit(data)

Fit the GLA model to trajectory data.

Returns: self

predict(n_ahead)

Predict future states using the extracted Koopman eigenvalues.

Returns: numpy.ndarray

reconstruct()

Reconstruct the training data from the GLA model.

Returns: numpy.ndarray

Properties

import numpy as np
from koopman_dmd import GLA

t = np.linspace(0, 20, 1000)
signal = np.cos(1.5 * t) + 0.3 * np.cos(3.7 * t)
data = signal.reshape(1, -1)

gla = GLA(n_eigenvalues=4, tol=1e-8).fit(data)

print("Eigenvalues:", gla.eigenvalues)
print("Convergence rate:", gla.convergence_rate)

recon = gla.reconstruct()
forecast = gla.predict(200)

Functions

generate_trajectory

generate_trajectory(ic, map_name, n_steps, params=None)

Generate a trajectory from one of the built-in dynamical maps.

Returns: numpy.ndarray of shape (n_vars, n_steps).

Map parameters:

from koopman_dmd import generate_trajectory

# Standard map trajectory
traj = generate_trajectory(
    ic=[0.1, 0.2],
    map_name="standard",
    n_steps=10000,
    params={"epsilon": 0.12}
)

# Henon map with custom parameters
traj = generate_trajectory(
    ic=[0.0, 0.0],
    map_name="henon",
    n_steps=5000,
    params={"a": 1.4, "b": 0.3}
)

# Logistic map
traj = generate_trajectory(
    ic=[0.5],
    map_name="logistic",
    n_steps=1000,
    params={"r": 3.9}
)

harmonic_time_average

harmonic_time_average(ic, map_name, observable, omega, n_iter, params=None)

Compute the harmonic time average of an observable along a trajectory at a specified frequency. The harmonic time average isolates the Fourier component of the observable at frequency omega, providing a finite-time approximation of the corresponding Koopman eigenfunction.

Returns: dict with keys:

from koopman_dmd import harmonic_time_average

result = harmonic_time_average(
    ic=[0.1, 0.2],
    map_name="standard",
    observable="identity",
    omega=0.3,
    n_iter=10000,
    params={"epsilon": 0.12}
)

print(f"Magnitude: {result['magnitude']:.6f}")
print(f"Phase:     {result['phase']:.6f}")
print(f"Real:      {result['real']:.6f}")
print(f"Imag:      {result['imag']:.6f}")

mesochronic_compute

mesochronic_compute(map_name, x_range, y_range, resolution, observable, omega, n_iter, params=None)

Compute the harmonic time average over a grid of initial conditions, producing matrices suitable for mesochronic analysis and visualization. This function evaluates the HTA at each grid point in the specified rectangular domain.

Returns: dict with keys:

from koopman_dmd import mesochronic_compute

result = mesochronic_compute(
    map_name="standard",
    x_range=(0.0, 1.0),
    y_range=(0.0, 1.0),
    resolution=200,
    observable="identity",
    omega=0.0,
    n_iter=5000,
    params={"epsilon": 0.12}
)

# Visualize with matplotlib
import matplotlib.pyplot as plt
plt.pcolormesh(
    result["x_coords"],
    result["y_coords"],
    result["hta_matrix"],
    cmap="viridis"
)
plt.colorbar(label="HTA magnitude")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Mesochronic plot -- Standard Map")
plt.show()

classify_phase_space

classify_phase_space(ic, map_name, omega, n_iter, params=None)

Classify a single initial condition as belonging to a regular (quasiperiodic) or chaotic region of phase space, based on the convergence behavior of the harmonic time average.

Returns: dict with keys:

from koopman_dmd import classify_phase_space

result = classify_phase_space(
    ic=[0.1, 0.2],
    map_name="standard",
    omega=0.0,
    n_iter=50000,
    params={"epsilon": 0.12}
)

print(f"Classification: {result['classification']}")
print(f"HTA magnitude:  {result['hta_magnitude']:.6f}")

hta_convergence

hta_convergence(ic, map_name, observable, omega, max_iter, step, params=None)

Compute the harmonic time average at increasing iteration counts to study its convergence behavior. Returns the HTA value at evenly spaced checkpoints, which can be used to assess whether the average has stabilized or is still fluctuating (indicating chaotic dynamics).

Returns: dict -- contains convergence data (iteration counts and corresponding HTA values at each checkpoint).

from koopman_dmd import hta_convergence

conv = hta_convergence(
    ic=[0.1, 0.2],
    map_name="standard",
    observable="identity",
    omega=0.3,
    max_iter=50000,
    step=500,
    params={"epsilon": 0.12}
)

Complete Examples

DMD: Fit, Predict, and Analyze

import numpy as np
from koopman_dmd import DMD

# Generate a multi-component signal
n = 200
dt = 0.02
t = np.arange(n) * dt
x1 = np.sin(2.0 * np.pi * 1.5 * t)
x2 = np.cos(2.0 * np.pi * 1.5 * t) * np.exp(-0.1 * t)
x3 = np.sin(2.0 * np.pi * 3.0 * t) * 0.5
data = np.vstack([x1, x2, x3])  # shape (3, 200)

# Fit the DMD model
dmd = DMD(rank=4, center=True).fit(data)

# Inspect eigenvalues and modes
print("Rank:       ", dmd.rank)
print("Eigenvalues:", dmd.eigenvalues)
print("Modes shape:", dmd.modes.shape)
print("Amplitudes: ", dmd.amplitudes)

# Spectral analysis
spec = dmd.spectrum(dt=dt)
for i, mode in enumerate(spec):
    print(f"Mode {i}: freq={mode['frequency']:.2f} Hz, "
          f"growth={mode['growth_rate']:.4f}, {mode['stability']}")

# Stability check
stab = dmd.stability()
print(f"System stable: {stab['is_stable']}")

# Reconstruction error
err = dmd.error(data)
print(f"RMSE:           {err['rmse']:.6f}")
print(f"Relative error: {err['relative_error']:.6f}")

# Dominant modes
top = dmd.dominant_modes(criterion="energy", n=2)
print("Top 2 modes by energy:", top)

# Mode residuals
res = dmd.residual()
print(f"Max residual:  {res['max_residual']:.2e}")

# Predict 100 steps ahead
forecast = dmd.predict(100)
print(f"Forecast shape: {forecast.shape}")

# Reconstruct the training window
fitted = dmd.reconstruct(n)
print(f"Reconstruction shape: {fitted.shape}")

Extended DMD with Polynomial Lifting

import numpy as np
from koopman_dmd import DMD

# Nonlinear signal: sin^2(t)
n = 300
dt = 0.05
t = np.arange(n) * dt
signal = np.sin(t) ** 2
data = signal.reshape(1, -1)

# Standard DMD (will have high error on nonlinear signal)
dmd_std = DMD().fit(data)
err_std = dmd_std.error(data)
print(f"Standard DMD RMSE: {err_std['rmse']:.6f}")

# Extended DMD with polynomial lifting (degree 2)
dmd_ext = DMD(lifting="polynomial", lifting_param=2).fit(data)
err_ext = dmd_ext.error(data)
print(f"Extended DMD RMSE: {err_ext['rmse']:.6f}")

# Predictions are automatically back-projected to original space
forecast = dmd_ext.predict(100)
print(f"Forecast shape: {forecast.shape}")  # (1, 100)

Hankel-DMD for Scalar Time Series

import numpy as np
from koopman_dmd import HankelDMD

# Two-frequency signal from scalar measurement
t = np.linspace(0, 10, 500)
signal = np.sin(2.0 * np.pi * t) + 0.5 * np.sin(6.0 * np.pi * t)
data = signal.reshape(1, -1)

# Hankel-DMD with time-delay embedding
hdmd = HankelDMD(delays=20, rank=4, dt=t[1] - t[0]).fit(data)

print("Eigenvalues:", hdmd.eigenvalues)
print("Rank:", hdmd.rank)

# Reconstruct and predict
recon = hdmd.reconstruct()
forecast = hdmd.predict(100)

print(f"Reconstruction shape: {recon.shape}")
print(f"Forecast shape:       {forecast.shape}")

GLA for Spectral Extraction

import numpy as np
from koopman_dmd import GLA

# Signal with two distinct frequencies
t = np.linspace(0, 20, 1000)
signal = np.cos(1.5 * t) + 0.3 * np.cos(3.7 * t)
data = signal.reshape(1, -1)

# Extract 4 Koopman eigenvalues via GLA
gla = GLA(n_eigenvalues=4, tol=1e-8).fit(data)

print("Koopman eigenvalues:")
for ev in gla.eigenvalues:
    print(f"  {ev}")
print(f"Convergence rate: {gla.convergence_rate:.2e}")

# Reconstruct and forecast
recon = gla.reconstruct()
forecast = gla.predict(200)
print(f"Forecast shape: {forecast.shape}")

Trajectory Generation and HTA Analysis

import numpy as np
from koopman_dmd import (
    generate_trajectory,
    harmonic_time_average,
    classify_phase_space,
    hta_convergence,
    mesochronic_compute,
)

# Generate a trajectory on the standard map
traj = generate_trajectory(
    ic=[0.1, 0.2],
    map_name="standard",
    n_steps=10000,
    params={"epsilon": 0.12}
)
print(f"Trajectory shape: {traj.shape}")  # (2, 10000)

# Compute harmonic time average
hta = harmonic_time_average(
    ic=[0.1, 0.2],
    map_name="standard",
    observable="identity",
    omega=0.3,
    n_iter=10000,
    params={"epsilon": 0.12}
)
print(f"HTA magnitude: {hta['magnitude']:.6f}")
print(f"HTA phase:     {hta['phase']:.6f}")

# Classify the initial condition
cls = classify_phase_space(
    ic=[0.1, 0.2],
    map_name="standard",
    omega=0.0,
    n_iter=50000,
    params={"epsilon": 0.12}
)
print(f"Region: {cls['classification']}")

# Study convergence of the HTA
conv = hta_convergence(
    ic=[0.1, 0.2],
    map_name="standard",
    observable="identity",
    omega=0.3,
    max_iter=50000,
    step=500,
    params={"epsilon": 0.12}
)

# Compute mesochronic plot data
meso = mesochronic_compute(
    map_name="standard",
    x_range=(0.0, 1.0),
    y_range=(0.0, 1.0),
    resolution=100,
    observable="identity",
    omega=0.0,
    n_iter=5000,
    params={"epsilon": 0.12}
)
print(f"HTA matrix shape:   {meso['hta_matrix'].shape}")
print(f"Phase matrix shape: {meso['phase_matrix'].shape}")