From 56aa1727eeed8ea04a5418d8280ac8d1f2f70028 Mon Sep 17 00:00:00 2001 From: "Mani Saint-Victor, MD" Date: Wed, 6 May 2026 23:21:18 -0400 Subject: [PATCH 1/2] =?UTF-8?q?feat:=20v0.2.0=20=E2=80=94=20sequence=20hel?= =?UTF-8?q?pers,=20irregular-time=20stepping,=20Bayesian=20filtering,=20ML?= =?UTF-8?q?E=20fit?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Expands linoss-dynamics from a single-step physics primitive into a higher-level multi-use library for oscillatory state-space dynamics, while keeping the core NumPy-only. Added (NumPy-only core): - linoss_scan: sequence helper that runs linoss_step over a (T, n) input - damped_oscillator_closed_form: analytic step supporting variable dt (irregular sampling) - stability module: is_stable, eigvals_to_freq_damping, freq_damping_to_oscillator_block, period_from_omega, harmonic_stack - filters module: kalman_filter, rts_smoother Added (optional [probabilistic] extra; requires scipy): - discretize module: discretize_lti_with_noise (van Loan), discretize_control (ZOH), oscillator_mats - fit module: fit_oscillator_mle (L-BFGS-B parameter recovery) Docs and tutorials: - 8 runnable examples covering single-step, forced, damped, energy tracking, convergence, typed errors, sunspots fit, and Kalman phase tracking - docs/api.md: per-function API reference - docs/roadmap.md: public vision document (v0.3 candidates, v1.0 unified oscillator+attractor goal, v1.x spatial extensions, explicit non-goals) - README Vision section linking to roadmap - PROVENANCE: framing-only attribution; no shared code with upstream JAX impl - architecture.md: lineage section grounding the package in the research family Verification: - 119/119 pytest green (Python 3.11) - ruff check clean across src/tests/examples - compileall clean - python -m build produces wheel + sdist; twine check PASSED - Clean-venv install from local wheel: all 23 public symbols importable, smoke test green scipy-gated symbols use module-level __getattr__ to raise LinOSSError with install hint if scipy is absent, instead of bare AttributeError. --- CHANGELOG.md | 28 +- CITATION.cff | 8 +- PROVENANCE.md | 17 + README.md | 136 +++++++- docs/api.md | 553 ++++++++++++++++++++++++++++++ docs/architecture.md | 49 +++ docs/roadmap.md | 80 +++++ examples/01_single_step.py | 62 ++++ examples/02_forced.py | 69 ++++ examples/03_damped.py | 83 +++++ examples/04_energy.py | 71 ++++ examples/05_convergence.py | 74 ++++ examples/06_validation_errors.py | 85 +++++ examples/07_sunspots.py | 102 ++++++ examples/08_phase_tracking.py | 146 ++++++++ examples/README.md | 52 +++ pyproject.toml | 9 +- src/linoss_dynamics/__init__.py | 63 ++++ src/linoss_dynamics/continuous.py | 266 ++++++++++++++ src/linoss_dynamics/discretize.py | 259 ++++++++++++++ src/linoss_dynamics/filters.py | 285 +++++++++++++++ src/linoss_dynamics/fit.py | 191 +++++++++++ src/linoss_dynamics/solver.py | 151 ++++++++ src/linoss_dynamics/stability.py | 217 ++++++++++++ tests/test_continuous.py | 520 ++++++++++++++++++++++++++++ tests/test_discretize.py | 131 +++++++ tests/test_filters.py | 290 ++++++++++++++++ tests/test_fit.py | 293 ++++++++++++++++ tests/test_solver.py | 165 +++++++++ tests/test_stability.py | 310 +++++++++++++++++ 30 files changed, 4756 insertions(+), 9 deletions(-) create mode 100644 docs/api.md create mode 100644 docs/roadmap.md create mode 100644 examples/01_single_step.py create mode 100644 examples/02_forced.py create mode 100644 examples/03_damped.py create mode 100644 examples/04_energy.py create mode 100644 examples/05_convergence.py create mode 100644 examples/06_validation_errors.py create mode 100644 examples/07_sunspots.py create mode 100644 examples/08_phase_tracking.py create mode 100644 examples/README.md create mode 100644 src/linoss_dynamics/continuous.py create mode 100644 src/linoss_dynamics/discretize.py create mode 100644 src/linoss_dynamics/filters.py create mode 100644 src/linoss_dynamics/fit.py create mode 100644 src/linoss_dynamics/stability.py create mode 100644 tests/test_continuous.py create mode 100644 tests/test_discretize.py create mode 100644 tests/test_filters.py create mode 100644 tests/test_fit.py create mode 100644 tests/test_stability.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 8e4dbaa..a369297 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,32 @@ # Changelog -## 0.1.0.dev0 - Unreleased +## [Unreleased] + +## [0.2.0] — 2026-05-06 + +### Added +- `linoss_scan(U, A, dt, ...)` — sequence-level helper that runs `linoss_step` over a `(T, n)` input and returns the full trajectory. +- `damped_oscillator_closed_form(y, z, omega, gamma, dt)` — analytic closed-form damped harmonic oscillator step supporting variable `dt` (irregular sampling). Handles underdamped, critically damped, and overdamped regimes via exact matrix exponential. +- `stability` module: `is_stable`, `eigvals_to_freq_damping`, `freq_damping_to_oscillator_block`, `period_from_omega`, `harmonic_stack`. +- Optional `[probabilistic]` extra (requires scipy): + - `discretize_lti_with_noise`, `discretize_control`, `oscillator_mats` — exact matrix-exponential discretization helpers (van Loan's method). + - `kalman_filter`, `rts_smoother` — Bayesian state estimation over oscillator state-space models. + - `fit_oscillator_mle` — MLE parameter recovery from observed time series via log-space optimization. +- `examples/` directory with 8 runnable tutorials covering single-step, forced, damped, energy, convergence, validation errors, sunspots, and phase tracking. +- `docs/api.md` — per-function API reference. +- `docs/roadmap.md` — public vision document. + +### Changed +- `__init__.py` — `__all__` extended with all new public symbols (23 total). Optional scipy-gated symbols included unconditionally; modules raise `LinOSSError` at call time if scipy is missing. +- `README.md` — Public API table extended; new Vision section linking to roadmap; new Quickstart subsections for scan, Kalman filtering, and MLE fitting. +- `PROVENANCE.md` — added framing-only attribution clarification and code-provenance section. +- `docs/architecture.md` — added Lineage section with mermaid timeline diagram. + +### Notes +- Core remains NumPy-only. SciPy is gated behind the `[probabilistic]` extra. +- v0.1.0 was never released to PyPI; v0.2.0 is the first PyPI release. + +## 0.1.0.dev0 — Unreleased (pre-release work) - Initial public-alpha split of `linoss-dynamics` as a standalone package. - NumPy-only LinOSS-style stepping helpers. diff --git a/CITATION.cff b/CITATION.cff index 4265858..b5cd70e 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -1,14 +1,16 @@ cff-version: 1.2.0 message: "If you use linoss-dynamics, please cite this package and the upstream LinOSS / D-LinOSS papers listed in PROVENANCE.md." title: "linoss-dynamics" -version: "0.1.0.dev0" -date-released: 2026-05-03 +version: "0.2.0" +date-released: 2026-05-06 repository-code: "https://github.com/bionicbutterfly13/linoss-dynamics" url: "https://github.com/bionicbutterfly13/linoss-dynamics" license: "MIT" type: software authors: - - name: "linoss-dynamics contributors" + - family-names: Saint-Victor + given-names: Mani + email: drmani215@gmail.com abstract: "A small NumPy runtime package for LinOSS-style oscillator dynamics with explicit non-negative damping support." keywords: - linoss diff --git a/PROVENANCE.md b/PROVENANCE.md index 4995cf8..dc194aa 100644 --- a/PROVENANCE.md +++ b/PROVENANCE.md @@ -25,6 +25,23 @@ Status: active provenance record for public alpha. The package is planned as MIT because it is a pure math/runtime primitive with no novelty claim over LinOSS or D-LinOSS. +## Code Provenance (Framing-only) + +`linoss-dynamics` is an original NumPy implementation of the LinOSS algorithm described in Rusch & Rus (2024, arXiv:2410.03943) and the D-LinOSS damping extension (Boyer, Rusch & Rus, 2025, arXiv:2505.12171). No code is shared with the upstream JAX/Equinox reference implementation at `tk-rusch/linoss`. + +**Upstream vs. this package:** + +| Aspect | `tk-rusch/linoss` (upstream) | `linoss-dynamics` (this package) | +| --- | --- | --- | +| Runtime | JAX + Equinox | NumPy only | +| API style | `LinOSSLayer`, `LinOSSBlock`, `LinOSS` (Equinox modules) | Functions and error classes | +| Entry points | `apply_linoss_im`, `apply_linoss_imex` (JAX scan ops) | `linoss_step`, `damped_linoss_step`, `linoss_scan` | +| Class name collision | None — namespaces do not overlap | — | + +**License obligation:** None. The LinOSS algorithm is described in the public papers above; this package is an independently written implementation. The upstream repository uses the MIT license; this package is also MIT-licensed, but there is no dependency or code-sharing relationship that creates an attribution obligation beyond the paper citations already present in `CITATION.cff`. + +**Continuous.py provenance:** `damped_oscillator_closed_form` implements the analytic matrix-exponential solution for the forced damped harmonic oscillator. This technique is standard in the classical oscillator state-space literature (Hartikainen & Särkkä 2010; Solin & Särkkä 2014) and was motivated by 2026 work on closed-form damped oscillators for irregular time series (see the deep-research dossier section "Recent research trends" for the underlying line of work; no arXiv id is cited here to avoid fabricating an identifier for a paper whose exact metadata has not been independently verified). + ## Evidence - Package core: `src/linoss_dynamics/` diff --git a/README.md b/README.md index 327436f..0b6517b 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,54 @@ `linoss-dynamics` is a small NumPy runtime package for LinOSS-style oscillator dynamics with explicit non-negative damping support. -Status: **public alpha**. The package is not published to PyPI yet. +Status: **public alpha**. + +## Vision + +`linoss-dynamics` is the runtime physics layer for oscillatory state-space dynamics. The current scope (v0.2) covers the deterministic oscillator step, irregular-time stepping, Bayesian filtering, and stability/frequency introspection. The longer roadmap — including unified oscillator + attractor primitives and spatio-temporal extensions — is documented in [docs/roadmap.md](docs/roadmap.md). + +## How this differs from the upstream LinOSS reference + +The original [tk-rusch/linoss](https://github.com/tk-rusch/linoss) is a +**training-time JAX/Equinox neural-network library** that exposes LinOSS as +trainable layers (`LinOSSLayer`, `LinOSSBlock`, `LinOSS`) inside a deep-learning +model. It is the right choice when you are building or training a sequence model +that uses LinOSS as a learnable encoder. + +`linoss-dynamics` is a **runtime physics package**. Same mathematics, different +audience and abstraction: + +| Dimension | Upstream `tk-rusch/linoss` | `linoss-dynamics` | +|---|---|---| +| Use case | Training neural networks with LinOSS layers | Runtime simulation, control, replay-safe systems | +| Stack | JAX + Equinox | NumPy only — zero deep-learning dependencies | +| Granularity | Sequence-level via JAX `scan` (`apply_linoss_im`, `apply_linoss_imex`) | Single-step (`linoss_step`, `damped_linoss_step`) | +| Damping | Implicit (rolled into `A`) | **Explicit, non-negativity validated** (`G` parameter) | +| Energy diagnostics | Not exposed | `energy`, `delta_energy`, `convergence_window` | +| Errors | Generic shape errors | Typed hierarchy (`LinOSSError`, `InvalidShapeError`, `InvalidDampingError`, `UnsupportedModeError`) | +| Determinism contract | Best-effort (training-oriented) | Replay-safe — pure NumPy, deterministic, no hidden state | +| Install footprint | Heavy (JAX, Equinox, dataset deps) | Tiny (NumPy 1.24+) | + +### Who needs `linoss-dynamics` specifically + +- **Replay-safe agentic systems** — needed inside deterministic kernels where + every step must be byte-reproducible (e.g. consumed by + [`elume`](https://github.com/bionicbutterfly13/elume) at runtime). +- **Simulation and control** — robotics, signal processing, vibration analysis, + and any context that uses oscillatory dynamics as a runtime model rather than + a learnable layer. +- **Reference implementation for porting** — a clean NumPy baseline to validate + custom JAX/PyTorch ports of LinOSS against. +- **Embedded or minimal-dependency contexts** — anywhere a JAX install is too + heavy or unavailable. +- **Teaching and pedagogy** — a self-contained, well-tested, single-file solver + that students can read end-to-end in an afternoon. +- **Energy-based analysis** — applications that need to monitor energy drift + and detect convergence; these diagnostics are not exposed by the upstream + reference. + +If you are training a neural network, use the upstream package. If you are +running a system that needs LinOSS as a deterministic physics step, use this one. ## What It Provides @@ -19,7 +66,13 @@ Status: **public alpha**. The package is not published to PyPI yet. ## Install -Install from GitHub: +From PyPI (recommended): + +```bash +pip install linoss-dynamics +``` + +From GitHub (latest main): ```bash python -m pip install "linoss-dynamics @ git+https://github.com/bionicbutterfly13/linoss-dynamics.git@main" @@ -66,6 +119,54 @@ y_next, z_next, metrics = damped_linoss_step(y, z, A, G, dt=0.1) assert metrics["damping_mode"] == "explicit_g" ``` +Sequence-level scan over a time series: + +```python +import numpy as np +from linoss_dynamics import linoss_scan + +U = np.random.default_rng(0).standard_normal((100, 1)) # (T, n) +A = np.array([1.0]) + +Y, Z, metrics = linoss_scan(U, A, dt=0.05) +print(Y.shape) # (101, 1) — includes initial state at Y[0] +``` + +Bayesian filtering over an observed time series (requires `[probabilistic]`): + +```python +# pip install linoss-dynamics[probabilistic] +import numpy as np +from linoss_dynamics import kalman_filter, rts_smoother, oscillator_mats + +Ad, Qd, Bd = oscillator_mats(beta=0.1, omega=2.0, sigma_proc=0.3, dt=0.05) +H = np.array([[1.0, 0.0]]) +R = np.array([[0.1]]) + +y_obs = np.sin(np.linspace(0, 10, 200)) + 0.1 * np.random.default_rng(42).standard_normal(200) +filt = kalman_filter(y_obs, Ad, Qd, H, R) +smooth = rts_smoother(Ad, filt) + +print(filt["m_f"].shape) # (200, 2) +print(smooth["m_s"].shape) # (200, 2) +``` + +Parameter fit on observed data (requires `[probabilistic]`): + +```python +# pip install linoss-dynamics[probabilistic] +import numpy as np +from linoss_dynamics import fit_oscillator_mle + +y = np.sin(np.linspace(0, 20, 400)) + 0.15 * np.random.default_rng(0).standard_normal(400) +result = fit_oscillator_mle(y, dt=0.05) + +print(f"fitted omega = {result['omega']:.3f} rad/s") +print(f"fitted period = {result['period']:.3f} s") +``` + +Runnable tutorials for all use cases are in [examples/](examples/README.md). + ## Numerical Contract `A` controls oscillator stiffness/frequency. @@ -88,19 +189,50 @@ Supported step modes: ## Public API +### Solver + | Name | Role | | --- | --- | | `linoss_step(y, z, A, dt, mode="implicit", B=None, u=None, damping=None, G=None)` | Advance one step, optionally dispatching to explicit damping when `damping` or `G` is supplied. | | `damped_linoss_step(y, z, A, G, dt, mode="implicit", B=None, u=None)` | Advance one step with explicit non-negative damping. | +| `linoss_scan(U, A, dt, mode="implicit", B=None, y0=None, z0=None, damping=None, G=None)` | Run `linoss_step` over a `(T, n)` input sequence; returns `(Y, Z, metrics)`. | | `energy(y, z, A)` | Return diagonal oscillator energy. | | `delta_energy(previous_energy, next_energy)` | Return signed energy delta. | | `convergence_window(deltas, threshold, window)` | Return true when recent absolute deltas are below a threshold. | | `linoss_step_impl(...)` | Backward-compatible alias for callers using the implementation name. | +### Stability + +| Name | Role | +| --- | --- | +| `is_stable(A, G=None, dt=None, mode="implicit")` | Return `(stable, reason)` for given oscillator parameters. | +| `eigvals_to_freq_damping(eigenvalues)` | Map complex eigenvalues to `(frequencies, damping_ratios)`. | +| `freq_damping_to_oscillator_block(omega, gamma, dt=1.0)` | Return a 2×2 real-form discrete oscillator block. | +| `period_from_omega(omega)` | Return `2π / omega`, vectorized. | +| `harmonic_stack(omegas, dampings=None, dt=1.0)` | Build a block-diagonal state matrix from `(omega, damping)` pairs. | + +### Continuous (irregular-time) + +| Name | Role | +| --- | --- | +| `damped_oscillator_closed_form(y, z, omega, gamma, dt, forcing=0.0)` | Analytic closed-form step using matrix exponential. Supports variable `dt` for irregular sampling. | + +### Probabilistic (requires `pip install linoss-dynamics[probabilistic]`) + +| Name | Role | +| --- | --- | +| `discretize_lti_with_noise(A_c, L, Qc, dt)` | Exact van Loan discretization — returns `(Ad, Qd)`. | +| `discretize_control(A_c, B, dt)` | ZOH discretization for control input — returns `(Ad, Bd)`. | +| `oscillator_mats(beta, omega, sigma_proc, dt=1.0, kappa=1.0)` | Build `(Ad, Qd, Bd)` for a single forced damped oscillator. | +| `kalman_filter(y, Ad, Qd, H, R, ...)` | Linear-Gaussian Kalman filter; returns filtered moments and log-likelihood. | +| `rts_smoother(Ad, filt)` | Rauch-Tung-Striebel backward smoother over Kalman filter outputs. | +| `fit_oscillator_mle(y, dt, ...)` | MLE parameter recovery for a damped oscillator SSM; returns fitted parameters and filter outputs. | + Public errors: | Error | Raised when | | --- | --- | +| `LinOSSError` | Base error for all LinOSS dynamics failures. | | `InvalidShapeError` | Inputs cannot be broadcast to the oscillator state. | | `InvalidDampingError` | Damping is outside the supported stable path, including negative `G`. | | `UnsupportedModeError` | The discretization mode is unsupported. | diff --git a/docs/api.md b/docs/api.md new file mode 100644 index 0000000..028bb69 --- /dev/null +++ b/docs/api.md @@ -0,0 +1,553 @@ +# API Reference + +Per-function reference for all public symbols in `linoss-dynamics`. Grouped by submodule. + +--- + +## Errors + +### `LinOSSError` + +```python +class LinOSSError(ValueError) +``` + +Base error for all LinOSS dynamics failures. Catch this to handle any library-level error. + +--- + +### `InvalidShapeError` + +```python +class InvalidShapeError(LinOSSError) +``` + +Raised when inputs (`y`, `z`, `A`, `G`, `B`, `u`) cannot be broadcast to the oscillator state shape. + +--- + +### `InvalidDampingError` + +```python +class InvalidDampingError(LinOSSError) +``` + +Raised when damping `G` or `gamma` contains negative values. + +--- + +### `UnsupportedModeError` + +```python +class UnsupportedModeError(LinOSSError) +``` + +Raised when `mode` is not a recognized discretization scheme. Valid values: `"implicit"` / `"im"`, `"implicit_explicit"` / `"imex"`. + +--- + +## Solver (`linoss_dynamics.solver`) + +### `linoss_step` + +```python +def linoss_step( + y: Any, + z: Any, + A: Any, + dt: float, + mode: str = "implicit", + B: Any | None = None, + u: Any | None = None, + damping: Any | None = None, + G: Any | None = None, +) -> tuple[np.ndarray, np.ndarray, dict[str, Any]] +``` + +Advance one LinOSS step. When `damping` or `G` is provided, dispatches to `damped_linoss_step`. + +**Parameters** + +| Name | Type | Description | +| --- | --- | --- | +| `y` | array-like | Displacement vector, length `n`. | +| `z` | array-like | Velocity vector, length `n`. | +| `A` | array-like | Stiffness/frequency — scalar, length-`n` vector, or `n×n` diagonal matrix. | +| `dt` | float | Time step. | +| `mode` | str | `"implicit"` (default) or `"implicit_explicit"`. | +| `B` | array-like or None | Forcing matrix, shape `(n, m)`. Optional. | +| `u` | array-like or None | Forcing input, shape `(m,)`. Optional. | +| `damping` / `G` | array-like or None | Non-negative damping. When either is set, uses `damped_linoss_step`. | + +**Returns** `(y_next, z_next, metrics)` — next state and a dict with `mode`, `energy_before`, `energy_after`, `delta_energy`, `damping_mode`. + +**Raises** `InvalidShapeError`, `InvalidDampingError`, `UnsupportedModeError`. + +**Example** + +```python +import numpy as np +from linoss_dynamics import linoss_step + +y_next, z_next, m = linoss_step(np.array([0.2]), np.array([1.0]), np.array([1.0]), dt=0.1) +``` + +--- + +### `damped_linoss_step` + +```python +def damped_linoss_step( + y: Any, + z: Any, + A: Any, + G: Any, + dt: float, + mode: str = "implicit", + B: Any | None = None, + u: Any | None = None, +) -> tuple[np.ndarray, np.ndarray, dict[str, Any]] +``` + +Advance one LinOSS step with explicit non-negative damping `G`. Negative values in `G` raise `InvalidDampingError`. + +**Parameters** Same as `linoss_step` plus `G` (required, non-negative scalar or vector). + +**Returns** `(y_next, z_next, metrics)`. Metrics include `damping_mode: "explicit_g"` and the resolved `damping` vector. + +**Raises** `InvalidShapeError`, `InvalidDampingError`, `UnsupportedModeError`. + +--- + +### `linoss_scan` + +```python +def linoss_scan( + U: Any, + A: Any, + dt: float, + *, + mode: str = "implicit", + B: Any | None = None, + y0: Any | None = None, + z0: Any | None = None, + damping: Any | None = None, + G: Any | None = None, +) -> tuple[np.ndarray, np.ndarray, list[dict[str, Any]]] +``` + +Run `linoss_step` over a time-varying input sequence, returning the full trajectory. + +**Parameters** + +| Name | Type | Description | +| --- | --- | --- | +| `U` | array-like `(T,)` or `(T, m)` | Input sequence. Each row `U[t]` is the forcing at step `t`. | +| `A` | array-like | Stiffness parameter. | +| `dt` | float | Time step. | +| `y0` / `z0` | array-like or None | Initial state. Defaults to zeros; `n` is inferred from `A`. | +| `damping` / `G` | array-like or None | Non-negative damping for all steps. | + +**Returns** `(Y, Z, metrics_seq)` where `Y` and `Z` have shape `(T+1, n)` including the initial state, and `metrics_seq` is a list of length `T`. + +**Raises** `InvalidShapeError`, `InvalidDampingError`, `UnsupportedModeError`. + +**Example** + +```python +import numpy as np +from linoss_dynamics import linoss_scan + +U = np.ones((50, 1)) +Y, Z, metrics = linoss_scan(U, A=np.array([1.0]), dt=0.05) +assert Y.shape == (51, 1) +``` + +--- + +### `linoss_step_impl` + +```python +def linoss_step_impl(*args: Any, **kwargs: Any) -> tuple[np.ndarray, np.ndarray, dict[str, Any]] +``` + +Backward-compatible alias for `linoss_step`. Accepts identical arguments and returns identical values. + +--- + +### `energy` + +```python +def energy(y: Any, z: Any, A: Any, G: Any | None = None) -> float +``` + +Return oscillator energy `0.5 * sum(A * y^2 + z^2)` for diagonal parameters. `G` is accepted but ignored (energy is defined by stiffness, not damping). + +**Parameters** `y`, `z` — state vectors; `A` — diagonal stiffness. + +**Returns** Scalar float. + +**Raises** `InvalidShapeError` if `y` and `z` have different shapes. + +--- + +### `delta_energy` + +```python +def delta_energy(previous_energy: float, next_energy: float) -> float +``` + +Return signed energy change `next_energy - previous_energy`. + +--- + +### `convergence_window` + +```python +def convergence_window( + deltas: list[float] | tuple[float, ...], + threshold: float, + window: int, +) -> bool +``` + +Return `True` when the `window` most recent absolute energy deltas are all below `threshold`. Useful for detecting steady-state convergence. + +**Parameters** + +| Name | Description | +| --- | --- | +| `deltas` | Sequence of signed energy deltas accumulated over time. | +| `threshold` | Convergence tolerance (positive). | +| `window` | Number of most-recent deltas to check (must be positive). | + +**Returns** `bool`. + +**Raises** `ValueError` if `window <= 0`. + +--- + +## Stability (`linoss_dynamics.stability`) + +### `is_stable` + +```python +def is_stable( + A: Any, + G: Any | None = None, + dt: float | None = None, + mode: str = "implicit", +) -> tuple[bool, str] +``` + +Return `(stable, reason)` for the given oscillator parameters. + +For `"implicit"` mode: stable when `A >= 0` and `G >= 0` (if provided). + +For `"implicit_explicit"` mode: stable when `A >= 0`, `G >= 0`, and `dt^2 * max(A) < 4` (CFL condition). Requires `dt`. + +**Returns** `(bool, str)` — stability flag and human-readable explanation. + +**Raises** `UnsupportedModeError`. + +**Example** + +```python +from linoss_dynamics import is_stable +stable, reason = is_stable(A=1.0, G=0.5, dt=0.1, mode="implicit_explicit") +assert stable +``` + +--- + +### `eigvals_to_freq_damping` + +```python +def eigvals_to_freq_damping(eigenvalues: Any) -> tuple[np.ndarray, np.ndarray] +``` + +Map complex eigenvalues `r·e^(±iθ)` to `(frequencies, damping_ratios)`. + +Frequency = `|θ|` (radians per step). Damping ratio = `-ln(r)` per step. + +**Returns** `(frequencies, damping_ratios)` — 1-D float arrays of length `N`. + +--- + +### `freq_damping_to_oscillator_block` + +```python +def freq_damping_to_oscillator_block( + omega: float, + gamma: float, + dt: float = 1.0, +) -> np.ndarray +``` + +Return a 2×2 real-form discrete state-transition block for a single harmonic oscillator with angular frequency `omega` and damping `gamma`. + +**Returns** Shape `(2, 2)` NumPy float array. + +**Raises** `InvalidShapeError` if `omega < 0`. + +--- + +### `period_from_omega` + +```python +def period_from_omega(omega: Any) -> np.ndarray +``` + +Return `2π / omega`, vectorized over array input. + +**Raises** `InvalidShapeError` if any element of `omega` is non-positive. + +--- + +### `harmonic_stack` + +```python +def harmonic_stack( + omegas: Any, + dampings: Any | None = None, + dt: float = 1.0, +) -> np.ndarray +``` + +Build a `(2N, 2N)` block-diagonal state matrix from `N` `(omega, damping)` pairs. Each pair yields a 2×2 oscillator block via `freq_damping_to_oscillator_block`. + +**Parameters** `omegas` — length-`N` frequencies; `dampings` — length-`N` damping coefficients (default: all zeros). + +**Returns** Shape `(2N, 2N)` float array. + +**Raises** `InvalidShapeError` if `omegas` and `dampings` have different lengths. + +--- + +## Continuous (`linoss_dynamics.continuous`) + +### `damped_oscillator_closed_form` + +```python +def damped_oscillator_closed_form( + y: Any, + z: Any, + omega: Any, + gamma: Any, + dt: float, + *, + forcing: float = 0.0, +) -> tuple[np.ndarray, np.ndarray] +``` + +Advance a damped harmonic oscillator by `dt` using the analytic closed-form solution of `ÿ + 2γẏ + ω²y = u`. + +Supports **variable `dt` between calls** — each call may use a different time interval, enabling exact stepping on irregular time grids. Handles all damping regimes: + +- **Underdamped** (`ω > γ`): oscillatory decay via `cos`/`sin`. +- **Critically damped** (`ω ≈ γ`): polynomial-exponential decay. +- **Overdamped** (`γ > ω`): two-exponential decay via `cosh`/`sinh`. + +**Parameters** + +| Name | Description | +| --- | --- | +| `y` | Position(s), scalar or 1-D array. | +| `z` | Velocity(ies), same broadcastable shape as `y`. | +| `omega` | Angular frequency(ies). Must be strictly positive. | +| `gamma` | Damping coefficient(s). Must be non-negative. | +| `dt` | Time interval to advance. Must be strictly positive. | +| `forcing` | Constant scalar forcing over the interval. Default `0.0`. | + +**Returns** `(y_next, z_next)` at time `t + dt`, both 1-D float64 arrays. + +**Raises** `ValueError` if `dt <= 0` or any `omega <= 0`. `InvalidDampingError` if any `gamma < 0`. `InvalidShapeError` if shapes are incompatible. + +**Example** + +```python +import numpy as np +from linoss_dynamics import damped_oscillator_closed_form + +y, z = np.array([1.0]), np.array([0.0]) +y1, z1 = damped_oscillator_closed_form(y, z, omega=1.0, gamma=0.0, dt=np.pi) +# half-period of pure cosine: y1 ≈ -1.0 +``` + +--- + +## Probabilistic (`linoss_dynamics.discretize`, `linoss_dynamics.filters`, `linoss_dynamics.fit`) + +These symbols require `pip install linoss-dynamics[probabilistic]` (SciPy). Importing without SciPy will succeed, but calling any function raises `LinOSSError` with an install hint. + +--- + +### `discretize_lti_with_noise` + +```python +def discretize_lti_with_noise( + A_c: np.ndarray, + L: np.ndarray, + Qc: np.ndarray, + dt: float, +) -> tuple[np.ndarray, np.ndarray] +``` + +Exact discretization of `dx = A_c·x·dt + L·dW` using van Loan's matrix-exponential method. + +**Parameters** `A_c` — `(n, n)` continuous state matrix; `L` — `(n, q)` noise input; `Qc` — `(q, q)` spectral density; `dt` — time step. + +**Returns** `(Ad, Qd)` — discrete state-transition and process-noise covariance, both `(n, n)`. + +**Raises** `LinOSSError` if SciPy is missing; `InvalidShapeError` for shape violations; `ValueError` if `dt <= 0`. + +--- + +### `discretize_control` + +```python +def discretize_control( + A_c: np.ndarray, + B: np.ndarray, + dt: float, +) -> tuple[np.ndarray, np.ndarray] +``` + +Zero-order-hold (ZOH) discretization of `dx/dt = A_c·x + B·u`. + +**Returns** `(Ad, Bd)` — discrete state-transition `(n, n)` and input `(n, m)`. + +**Raises** `LinOSSError` if SciPy is missing; `InvalidShapeError`; `ValueError` if `dt <= 0`. + +--- + +### `oscillator_mats` + +```python +def oscillator_mats( + beta: float, + omega: float, + sigma_proc: float, + dt: float = 1.0, + kappa: float = 1.0, +) -> tuple[np.ndarray, np.ndarray, np.ndarray] +``` + +Build `(Ad, Qd, Bd)` for a single forced damped oscillator with continuous-time dynamics `ẏ = z`, `ż = -ω²y - 2βz + κu + σ_proc·w`. + +**Parameters** `beta` — damping coefficient; `omega` — natural frequency; `sigma_proc` — process noise std dev; `dt` — time step; `kappa` — control gain. + +**Returns** `(Ad, Qd, Bd)` where `Ad` and `Qd` are `(2, 2)` and `Bd` is `(2, 1)`. + +**Raises** `LinOSSError` if SciPy is missing; `ValueError` if `dt <= 0`. + +**Example** + +```python +import numpy as np +from linoss_dynamics import oscillator_mats +Ad, Qd, Bd = oscillator_mats(beta=0.1, omega=2.0, sigma_proc=0.3, dt=0.05) +assert Ad.shape == (2, 2) +``` + +--- + +### `kalman_filter` + +```python +def kalman_filter( + y: np.ndarray, + Ad: np.ndarray, + Qd: np.ndarray, + H: np.ndarray, + R: np.ndarray, + *, + m0: np.ndarray | None = None, + P0: np.ndarray | None = None, + Bd: np.ndarray | None = None, + u: np.ndarray | None = None, +) -> dict[str, Any] +``` + +Run a linear-Gaussian Kalman filter over a 1-D observation sequence `y` of length `T`. + +**Parameters** + +| Name | Shape | Description | +| --- | --- | --- | +| `y` | `(T,)` | Observation sequence. | +| `Ad` | `(n, n)` | Discrete state-transition matrix. | +| `Qd` | `(n, n)` | Process-noise covariance. | +| `H` | `(1, n)` | Observation matrix. | +| `R` | `(1, 1)` | Observation-noise covariance. | +| `m0` | `(n,)` | Initial state mean. Default: zeros. | +| `P0` | `(n, n)` | Initial state covariance. Default: identity. | +| `Bd` | `(n, p)` | Control input matrix. Optional. | +| `u` | `(T,)` or `(T, p)` | Control sequence. Required when `Bd` is provided. | + +**Returns** Dict with keys `"m_f"` `(T, n)`, `"P_f"` `(T, n, n)`, `"m_p"` `(T, n)`, `"P_p"` `(T, n, n)`, `"loglik"` float. + +**Raises** `InvalidShapeError` for shape violations. + +--- + +### `rts_smoother` + +```python +def rts_smoother( + Ad: np.ndarray, + filt: dict[str, Any], +) -> dict[str, Any] +``` + +Run a Rauch-Tung-Striebel smoother backward through Kalman filter outputs. + +**Parameters** `Ad` — state-transition matrix `(n, n)`; `filt` — dict returned by `kalman_filter`. + +**Returns** Dict with keys `"m_s"` `(T, n)` and `"P_s"` `(T, n, n)`. + +**Raises** `InvalidShapeError` if `Ad` is not square or does not match state dimension in `filt`. + +**Note** At `t = T-1` the smoothed distribution equals the filtered distribution (RTS endpoint property). + +--- + +### `fit_oscillator_mle` + +```python +def fit_oscillator_mle( + y: np.ndarray, + dt: float, + *, + u: np.ndarray | None = None, + init: tuple[float, float, float, float] | None = None, + method: str = "L-BFGS-B", +) -> dict[str, Any] +``` + +Fit `(beta, omega, sigma_proc, sigma_obs)` of a single damped oscillator SSM by maximising the Kalman filter log-likelihood. Optimization is performed in log-parameter space so all parameters remain strictly positive. + +**Parameters** + +| Name | Description | +| --- | --- | +| `y` | Observation sequence `(T,)`. | +| `dt` | Time step (positive). | +| `u` | Optional scalar control sequence `(T,)`. Default: zeros. | +| `init` | Initial guess `(beta, omega, sigma_proc, sigma_obs)`. Default: `(0.10, 2.0, 0.40, 0.30)`. | +| `method` | SciPy `minimize` method. Default: `"L-BFGS-B"`. | + +**Returns** Dict with fitted parameters (`beta`, `omega`, `sigma_proc`, `sigma_obs`, `period`, `mean_y`), matrices (`Ad`, `Qd`, `Bd`, `H`, `R`), and filter outputs (`filt`, `smooth`, `optim`). + +**Raises** `LinOSSError` if SciPy is missing; `ValueError` if `dt <= 0`. + +**Example** + +```python +import numpy as np +from linoss_dynamics import fit_oscillator_mle + +y = np.sin(np.linspace(0, 20, 400)) + 0.15 * np.random.default_rng(0).standard_normal(400) +result = fit_oscillator_mle(y, dt=0.05) +assert 0 < result["omega"] < 20 +``` diff --git a/docs/architecture.md b/docs/architecture.md index 55a7674..0df3813 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -41,3 +41,52 @@ The package core must not depend on: Host applications should depend on this package through normal Python package installation and keep host-specific adapters outside the package core. + +## Module Map (v0.2.0) + +| Layer | Files | Responsibility | +| --- | --- | --- | +| Public API | `src/linoss_dynamics/__init__.py` | Re-export stable functions and errors. | +| Solver core | `src/linoss_dynamics/solver.py` | NumPy stepping, damping, energy, validation, convergence, and scan helpers. | +| Stability | `src/linoss_dynamics/stability.py` | Stability predicates, eigenvalue/frequency conversion, oscillator block builders. | +| Continuous | `src/linoss_dynamics/continuous.py` | Closed-form analytic step for irregular-time damped harmonic oscillators. | +| Discretize | `src/linoss_dynamics/discretize.py` | van Loan matrix-exponential discretization; ZOH control input discretization. Gated behind scipy. | +| Filters | `src/linoss_dynamics/filters.py` | Linear-Gaussian Kalman filter and RTS smoother. | +| Fit | `src/linoss_dynamics/fit.py` | MLE parameter recovery for damped oscillator SSMs. Gated behind scipy. | + +## Lineage + +This package sits in a research lineage that traces from Kalman (1960) through state-space GP and structural time-series work to the modern deep oscillatory SSM family (LinOSS, D-LinOSS). + +```mermaid +timeline + title Oscillatory SSM lineage + 1960 : Kalman + : Linear filtering and state estimation + 1985 : Harvey + : Stochastic cycle components in structural time series + 2010 : Hartikainen and Särkkä + : Temporal GP regression as LGSSM with KF and RTS + 2014 : Solin and Särkkä + : Periodic covariance functions linked to stochastic resonators + 2018 : Beck et al. + : Oscillator decomposition for neural data + 2020 : HiPPO + : Memory foundations for modern SSMs + 2020 : coRNN + : Neural oscillators for long dependencies + 2021 : UnICORNN + : Hamiltonian reversible oscillator RNN + 2021 : Wodeyar et al. and PLSO + : Real-time phase estimation and nonstationary oscillatory decomposition + 2021 : S4 + : Structured deep SSMs for long sequences + 2023 : S5 and Mamba + : Scan-friendly and selective first-order SSMs + 2024 : RON and LinOSS preprint + : Random oscillators and deep second-order oscillatory SSMs + 2025 : LinOSS oral, D-LinOSS, BioOSS + : Flexible damping and spatio-temporal wave dynamics + 2026 : Irregular-time oscillator-attention preprints + : Closed-form damped oscillators for irregular sequences +``` diff --git a/docs/roadmap.md b/docs/roadmap.md new file mode 100644 index 0000000..dab4c02 --- /dev/null +++ b/docs/roadmap.md @@ -0,0 +1,80 @@ +# Roadmap + +`linoss-dynamics` is a runtime physics package for oscillatory state-space dynamics. This roadmap documents the trajectory of the package — what's shipped, what's next, and what's explicitly out of scope. + +--- + +## Where we are: v0.2.0 (current) + +v0.2.0 is the first PyPI release. The core is NumPy-only; optional probabilistic tools require SciPy. + +**Shipped in v0.2.0:** + +- **Deterministic oscillator step** — `linoss_step` and `damped_linoss_step` implement the implicit (LinOSS-IM) and implicit-explicit (LinOSS-IMEX) discretizations of the forced harmonic oscillator with explicit non-negative damping `G`. +- **Sequence-level scan** — `linoss_scan` runs the step function over a `(T, n)` input sequence and returns the full trajectory. +- **Closed-form irregular-time stepping** — `damped_oscillator_closed_form` solves the damped second-order ODE analytically using the matrix exponential. Because each call can use a different `dt`, this enables exact stepping on irregular time grids without accumulated discretization bias. +- **Stability and frequency tools** — `is_stable`, `eigvals_to_freq_damping`, `freq_damping_to_oscillator_block`, `period_from_omega`, `harmonic_stack`. +- **Kalman filter and RTS smoother** — `kalman_filter` and `rts_smoother` provide Bayesian state estimation over discrete-time linear-Gaussian oscillator SSMs. +- **Exact discretization** — `discretize_lti_with_noise` (van Loan's method) and `discretize_control` (ZOH) convert continuous-time oscillator parameters to discrete-time matrices. `oscillator_mats` is a convenience wrapper for single forced damped oscillators. +- **MLE parameter fit** — `fit_oscillator_mle` recovers `(beta, omega, sigma_proc, sigma_obs)` from an observed time series by maximising the Kalman filter log-likelihood. +- **8 runnable tutorials** covering single-step, forced oscillators, damping, energy, convergence, validation errors, sunspot modeling, and phase tracking. + +--- + +## v0.3.0 — Candidate next-version work + +Items below are research-grounded extensions identified in the oscillatory SSM literature. Each is a candidate for v0.3.0 inclusion, with scope notes. + +### PLSO (piecewise locally stationary oscillatory) regimes + +PLSO (Wodeyar et al. 2021) models time series as piecewise-stationary oscillatory mixtures, using Kalman-based inference with proximal-gradient regularization over piecewise windows. The runtime-physics analog is a `switch_oscillator` primitive that can change frequency and damping parameters at detected change points while maintaining exact matrix-exponential transitions at each segment boundary. This builds naturally on `damped_oscillator_closed_form` and `oscillator_mats`. Effort: medium. No new external dependencies if the change-point detection is left to the caller. + +### Switching state-space models + +Switching SSMs extend the PLSO idea to a latent discrete mode variable, where the oscillator transitions between a small set of modes (e.g., spindle vs. non-spindle in neuroscience, high-frequency vs. low-frequency in economics). The approximate inference stack (Kim's algorithm or variational) is non-trivial but well-studied. A v0.3.0 candidate would provide the oscillator-mode building blocks and leave the switching inference to a thin API the caller can plug into. This is a natural extension of `kalman_filter` and `rts_smoother`. Effort: high. + +### DMD / Koopman-DMD spectral identification + +Dynamic Mode Decomposition (Schmid 2010) identifies oscillatory spatial modes and their temporal evolution from snapshot data. Kalman-filter DMD and Koopman-autoencoder approaches extend this to data assimilation and online tracking. A v0.3.0 candidate would provide a pure-NumPy DMD routine that returns eigenvalues compatible with `eigvals_to_freq_damping` and blocks compatible with `harmonic_stack`, closing the loop between data-driven identification and physics-stepping. Effort: medium. + +### Koopman-Kalman hybrids for online tracking + +Recent 2024–2025 work combines Koopman latent linearizations with Kalman or ensemble-Kalman updates to track slowly changing oscillatory systems online. The linoss-dynamics angle would be a `tracking_filter` that accepts a time-varying `Ad` sequence (updated by a Koopman estimator) and runs a standard Kalman recursion — the oscillator stepping provides the physics prior, the Kalman filter handles the estimation. Effort: low (once DMD identification is in place). + +--- + +## v1.0 — Path B: unified oscillators + attractor basins + +The largest open question in this research lineage is whether oscillator dynamics (LinOSS-family) and attractor memory (Hopfield-family) should converge into one substrate or remain composed primitives. + +The oscillatory SSM literature identifies this as the synthesis direction: "uncertainty-calibrated, continuous-time or irregular-time, second-order oscillatory SSMs that remain competitive at the scale of modern sequence benchmarks" (deep-research dossier, "Recent research trends"). The runtime-physics version of that synthesis is a unified oscillator + attractor library. + +`linoss-dynamics` v1.0 is committed to providing the unified primitives. Until then, downstream consumers compose `linoss-dynamics` with their own attractor layer (Hopfield/CAN-style). The API boundary is designed so that attractor state can be stored alongside oscillator state without structural changes — both are plain NumPy arrays passed explicitly. + +--- + +## v1.x — Spatio-temporal extensions (BioOSS-class) + +Spatio-temporal oscillator models (BioOSS, 2025) extend the state to a 2D field with lateral coupling between oscillator sites, generating wave-like propagation across a structured substrate. This breaks the 1D state contract of the current API — `y` and `z` are flat vectors, not field arrays — and is therefore an explicit non-goal for v0.2–1.0. + +v1.x will introduce a separate `spatial.*` namespace with field-shaped state and a 2D stepping API. The core scalar stepping functions will remain unchanged. + +--- + +## Explicit non-goals + +The following are intentionally outside the scope of this package — use the cited alternatives: + +| Goal | Alternative | +| --- | --- | +| Training-time JAX neural networks with LinOSS layers | Use the upstream `tk-rusch/linoss` JAX/Equinox reference. | +| Full probabilistic GP inference with arbitrary kernels | Use GPflow, GPyTorch, or Stheno. | +| Web or service runtimes | Host applications integrate via normal Python imports. | +| Agent frameworks | `linoss-dynamics` provides primitives; agentic composition belongs in the host application. | +| Symbolic or autodiff-based parameter learning | Use JAX or PyTorch; this package is NumPy-only. | + +--- + +## Citations + +This roadmap draws from the LinOSS and D-LinOSS papers (Rusch & Rus 2024; Boyer, Rusch & Rus 2025), the classical oscillator-SSM literature (Hartikainen & Särkkä 2010; Solin & Särkkä 2014; Beck et al. 2018; Wodeyar et al. 2021), and the 2025–2026 research line on continuous-time and uncertainty-calibrated oscillatory SSMs. See `PROVENANCE.md` for full attribution and `CITATION.cff` for machine-readable citation metadata. diff --git a/examples/01_single_step.py b/examples/01_single_step.py new file mode 100644 index 0000000..3b1fa07 --- /dev/null +++ b/examples/01_single_step.py @@ -0,0 +1,62 @@ +"""01_single_step.py — LinOSS single-step basics. + +Demonstrates how to advance a small oscillator state by one timestep using +both supported discretization modes and the backward-compatible alias. +""" + +from __future__ import annotations + +import numpy as np + +from linoss_dynamics import linoss_step, linoss_step_impl + + +def main() -> None: + # A 3-oscillator system: stiffness A, displacement y, velocity z. + # A controls the squared natural frequency for each oscillator. + rng = np.random.default_rng(0) + A = np.array([1.0, 4.0, 9.0]) # ω² values: ω = 1, 2, 3 rad/s + y = rng.standard_normal(3) # random initial displacement + z = rng.standard_normal(3) # random initial velocity + dt = 0.1 + + print("=== 01 Single Step ===") + print(f"A (stiffness): {A}") + print(f"y (displacement before): {y}") + print(f"z (velocity before): {z}") + print() + + # --- implicit mode (default, unconditionally stable) --- + y_im, z_im, metrics_im = linoss_step(y, z, A, dt=dt, mode="implicit") + + print("--- implicit mode ---") + print(f"y after: {y_im}") + print(f"z after: {z_im}") + print(f"energy before: {metrics_im['energy_before']:.6f}") + print(f"energy after: {metrics_im['energy_after']:.6f}") + print(f"delta_energy: {metrics_im['delta_energy']:.6f}") + print(f"mode: {metrics_im['mode']}") + print() + + # --- implicit-explicit (IMEX) mode --- + # Stability condition: dt^2 * max(A) < 4 + # Here: 0.01 * 9 = 0.09 << 4, so we're well inside the CFL limit. + y_imex, z_imex, metrics_imex = linoss_step(y, z, A, dt=dt, mode="implicit_explicit") + + print("--- implicit_explicit (IMEX) mode ---") + print(f"y after: {y_imex}") + print(f"z after: {z_imex}") + print(f"energy before: {metrics_imex['energy_before']:.6f}") + print(f"energy after: {metrics_imex['energy_after']:.6f}") + print(f"delta_energy: {metrics_imex['delta_energy']:.6f}") + print() + + # --- linoss_step_impl is a backward-compatible alias --- + y_alias, z_alias, _ = linoss_step_impl(y, z, A, dt=dt) + print("--- linoss_step_impl (alias check) ---") + print(f"y matches linoss_step result: {np.allclose(y_alias, y_im)}") + print(f"z matches linoss_step result: {np.allclose(z_alias, z_im)}") + + +if __name__ == "__main__": + main() diff --git a/examples/02_forced.py b/examples/02_forced.py new file mode 100644 index 0000000..a669be4 --- /dev/null +++ b/examples/02_forced.py @@ -0,0 +1,69 @@ +"""02_forced.py — Forced oscillator trajectory with linoss_scan. + +Builds a sinusoidal forcing sequence and runs it through a small LinOSS system +to produce the full (T+1) trajectory. Energy is tracked at every step. +""" + +from __future__ import annotations + +import numpy as np + +from linoss_dynamics import energy +from linoss_dynamics.solver import linoss_scan + + +def ascii_plot(values: list[float], width: int = 50, label: str = "") -> None: + """Minimal ASCII bar chart — positive values only, normalised to `width`.""" + lo, hi = min(values), max(values) + span = hi - lo or 1.0 + print(f" {label} (min={lo:.3f}, max={hi:.3f})") + for i, v in enumerate(values): + bar_len = int((v - lo) / span * width) + print(f" t={i:3d} |{'#' * bar_len}") + + +def main() -> None: + # 2-oscillator system + n = 2 + A = np.array([1.0, 4.0]) # ω² for each oscillator + dt = 0.05 + T = 80 # number of timesteps + + # Sinusoidal forcing sequence — shape (T, n) + t_vec = np.arange(T) * dt + U = np.column_stack([ + np.sin(2.0 * np.pi * 0.5 * t_vec), # 0.5 Hz drive for oscillator 0 + np.cos(2.0 * np.pi * 1.0 * t_vec), # 1.0 Hz drive for oscillator 1 + ]) + + print("=== 02 Forced Oscillator ===") + print(f"n={n} oscillators, T={T} steps, dt={dt}") + print(f"U shape: {U.shape}") + + # Run the scan — returns full trajectory including initial state at index 0 + Y, Z, metrics = linoss_scan(U, A, dt=dt, mode="implicit") + + print(f"Y shape: {Y.shape} (T+1 rows, one per timestep + initial state)") + print(f"Z shape: {Z.shape}") + print(f"metrics list length: {len(metrics)} (one dict per step)") + print() + + # Compute energy at every step using the energy() helper + energies = [energy(Y[t], Z[t], A) for t in range(T + 1)] + print(f"Energy at t=0: {energies[0]:.6f}") + print(f"Energy at t={T}: {energies[-1]:.6f}") + print() + + # Print every 10th step so the output is readable + print("Step-by-step excerpt (every 10th timestep):") + print(f"{'t':>4} {'y[0]':>10} {'y[1]':>10} {'energy':>12}") + for t in range(0, T + 1, 10): + print(f"{t:4d} {Y[t, 0]:10.4f} {Y[t, 1]:10.4f} {energies[t]:12.6f}") + print() + + # ASCII energy plot (every 5th step to keep output manageable) + ascii_plot(energies[::5], width=40, label="Energy profile (every 5th step)") + + +if __name__ == "__main__": + main() diff --git a/examples/03_damped.py b/examples/03_damped.py new file mode 100644 index 0000000..76afe25 --- /dev/null +++ b/examples/03_damped.py @@ -0,0 +1,83 @@ +"""03_damped.py — Damped vs undamped oscillator comparison. + +Uses linoss_scan with a damping parameter G to show energy decay, then +validates against the analytic closed-form solution from continuous.py. +""" + +from __future__ import annotations + +import numpy as np + +from linoss_dynamics import energy +from linoss_dynamics.continuous import damped_oscillator_closed_form +from linoss_dynamics.solver import linoss_scan + + +def run_scan( + A: np.ndarray, + T: int, + dt: float, + *, + damping: np.ndarray | None = None, +) -> tuple[np.ndarray, np.ndarray, list[float]]: + """Run T forced steps; return Y, Z, per-step energies.""" + # Scalar forcing — impulse at t=0 only + U = np.zeros((T, A.size)) + U[0] = 1.0 + + Y, Z, _ = linoss_scan(U, A, dt=dt, damping=damping) + energies = [energy(Y[t], Z[t], A) for t in range(T + 1)] + return Y, Z, energies + + +def main() -> None: + # Single oscillator for easy comparison with the analytic solution. + # omega = sqrt(A) = 2 rad/s, so period ≈ 3.14 s. + omega = 2.0 + A = np.array([omega**2]) # LinOSS stiffness = omega^2 + gamma = 0.15 # damping coefficient + G = np.array([gamma]) # LinOSS explicit damping vector + dt = 0.05 + T = 50 + + print("=== 03 Damped vs Undamped ===") + print(f"omega={omega} rad/s, gamma={gamma}, dt={dt}, T={T} steps") + print() + + Y_un, Z_un, energies_un = run_scan(A, T, dt) + Y_dm, Z_dm, energies_dm = run_scan(A, T, dt, damping=G) + + print("Energy at key timesteps:") + print(f"{'t':>5} {'undamped':>14} {'damped':>14} {'ratio':>10}") + for t in [0, 10, 25, 50]: + e_un = energies_un[t] + e_dm = energies_dm[t] + ratio = e_dm / (e_un or 1e-30) + print(f"{t:5d} {e_un:14.6f} {e_dm:14.6f} {ratio:10.4f}") + print() + + # Analytic comparison — closed-form exact solution for the damped oscillator. + # After the impulse step at t=0, Y_dm[1] and Z_dm[1] hold the first true state. + # We advance the analytic solver from that state and compare with the scan. + y_a = Y_dm[1].copy() + z_a = Z_dm[1].copy() + + print("Analytic closed-form comparison (selected checkpoints):") + print(f"{'step':>5} {'LinOSS y':>12} {'analytic y':>12} {'abs diff':>12}") + + for t in range(2, T + 1): + y_a, z_a = damped_oscillator_closed_form( + y_a, z_a, omega=omega, gamma=gamma, dt=dt + ) + if t in {5, 10, 20, 35, 50}: + linoss_y = float(Y_dm[t, 0]) + diff = abs(linoss_y - float(y_a[0])) + print(f"{t:5d} {linoss_y:12.6f} {float(y_a[0]):12.6f} {diff:12.2e}") + + print() + print("Note: small discrepancies are expected — LinOSS uses a discrete") + print("approximation while the closed-form is exact for the continuous ODE.") + + +if __name__ == "__main__": + main() diff --git a/examples/04_energy.py b/examples/04_energy.py new file mode 100644 index 0000000..29169fc --- /dev/null +++ b/examples/04_energy.py @@ -0,0 +1,71 @@ +"""04_energy.py — Energy tracking over 100 steps. + +Shows how to use energy() and delta_energy() to monitor whether a damped +oscillator is losing energy monotonically. +""" + +from __future__ import annotations + +import numpy as np + +from linoss_dynamics import ( + damped_linoss_step, + delta_energy, + energy, +) + + +def main() -> None: + rng = np.random.default_rng(0) + + # 2-oscillator system with mild damping — no external forcing. + A = np.array([1.0, 4.0]) + G = np.array([0.2, 0.3]) # per-oscillator damping + dt = 0.1 + N = 100 + + # Start from a non-zero state. + y = rng.standard_normal(2) * 2.0 + z = rng.standard_normal(2) * 1.0 + + print("=== 04 Energy Tracking ===") + print(f"A={A}, G={G}, dt={dt}, N={N} steps") + print(f"Initial y={y}, z={z}") + print() + + # Step forward, recording energy at each step. + energies: list[float] = [energy(y, z, A)] + + for _ in range(N): + y, z, _ = damped_linoss_step(y, z, A, G, dt=dt) + energies.append(energy(y, z, A)) + + # Compute delta_energy between consecutive steps. + deltas = [delta_energy(energies[t], energies[t + 1]) for t in range(N)] + + print(f"Energy at step 0: {energies[0]:.6f}") + print(f"Energy at step {N}: {energies[N]:.6f}") + print(f"Total energy lost: {energies[0] - energies[N]:.6f}") + print() + print(f"delta_energy at step 0→1: {deltas[0]:.6f} (should be negative)") + print(f"delta_energy at step 1→2: {deltas[1]:.6f}") + print() + + # Verify monotonic energy decay — every delta should be <= 0. + all_negative = all(d <= 1e-12 for d in deltas) + print(f"All delta_energy values <= 0: {all_negative}") + if not all_negative: + bad = [(i, d) for i, d in enumerate(deltas) if d > 1e-12] + print(f" Violations at steps: {bad}") + + # Print every 10th energy value so the decay profile is visible. + print() + print("Energy profile (every 10th step):") + print(f"{'step':>5} {'energy':>14} {'cumulative delta':>18}") + for t in range(0, N + 1, 10): + cum_delta = energies[t] - energies[0] + print(f"{t:5d} {energies[t]:14.6f} {cum_delta:18.6f}") + + +if __name__ == "__main__": + main() diff --git a/examples/05_convergence.py b/examples/05_convergence.py new file mode 100644 index 0000000..7907bd2 --- /dev/null +++ b/examples/05_convergence.py @@ -0,0 +1,74 @@ +"""05_convergence.py — Detecting convergence with convergence_window. + +Runs a heavily damped oscillator and uses convergence_window() to detect +when the absolute energy change has been below a threshold for 5 consecutive +steps, signalling that the system has effectively settled to rest. +""" + +from __future__ import annotations + +import numpy as np + +from linoss_dynamics import ( + convergence_window, + damped_linoss_step, + delta_energy, + energy, +) + + +def main() -> None: + rng = np.random.default_rng(0) + + # Heavy damping drives the oscillator to rest quickly. + A = np.array([1.0, 2.25]) + G = np.array([0.8, 1.0]) + dt = 0.1 + max_steps = 500 + threshold = 1e-6 + window = 5 # convergence requires 5 consecutive |Δenergy| < threshold + + y = rng.standard_normal(2) + z = rng.standard_normal(2) + + print("=== 05 Convergence Window ===") + print(f"A={A}, G={G}, dt={dt}") + print(f"Convergence: |Δenergy| < {threshold} for {window} consecutive steps") + print(f"Initial y={y}, z={z}") + print(f"Initial energy: {energy(y, z, A):.6f}") + print() + + # Accumulate absolute energy deltas. + prev_e = energy(y, z, A) + abs_deltas: list[float] = [] + converged_at: int | None = None + + for step in range(1, max_steps + 1): + y, z, _ = damped_linoss_step(y, z, A, G, dt=dt) + curr_e = energy(y, z, A) + abs_deltas.append(abs(delta_energy(prev_e, curr_e))) + prev_e = curr_e + + # convergence_window only checks the last `window` entries. + if convergence_window(abs_deltas, threshold=threshold, window=window): + converged_at = step + break + + if converged_at is not None: + print(f"Convergence detected at step {converged_at}") + print(f"Final energy: {energy(y, z, A):.2e}") + print(f"Last {window} |Δenergy| values: " + f"{[f'{v:.2e}' for v in abs_deltas[-window:]]}") + else: + print(f"Did not converge within {max_steps} steps") + print(f"Final energy: {energy(y, z, A):.6f}") + + # Also show the early high-energy region so the reader can see the decay. + print() + print("Early absolute delta_energy values:") + for i, d in enumerate(abs_deltas[:10], start=1): + print(f" step {i:3d}: |Δenergy| = {d:.4e}") + + +if __name__ == "__main__": + main() diff --git a/examples/06_validation_errors.py b/examples/06_validation_errors.py new file mode 100644 index 0000000..e9db71b --- /dev/null +++ b/examples/06_validation_errors.py @@ -0,0 +1,85 @@ +"""06_validation_errors.py — Typed validation errors in linoss_dynamics. + +Each error subclass targets a specific class of invalid input so callers can +catch fine-grained exceptions. This script triggers each one intentionally and +catches it to show the error message. +""" + +from __future__ import annotations + +import numpy as np + +from linoss_dynamics import ( + InvalidDampingError, + InvalidShapeError, + LinOSSError, + UnsupportedModeError, + linoss_step, +) + + +def section(title: str) -> None: + print(f"\n--- {title} ---") + + +def main() -> None: + print("=== 06 Validation Errors ===") + print() + print("Error hierarchy:") + print(" LinOSSError (base)") + print(" ├── InvalidShapeError") + print(" ├── InvalidDampingError") + print(" └── UnsupportedModeError") + + # 1. InvalidShapeError — y and z have mismatched lengths + section("InvalidShapeError: mismatched y and z") + y_bad = np.array([1.0, 2.0]) # length 2 + z_bad = np.array([0.0, 0.0, 0.0]) # length 3 — mismatch! + A = np.array([1.0, 4.0, 9.0]) + try: + linoss_step(y_bad, z_bad, A, dt=0.1) + except InvalidShapeError as exc: + print(f"Caught InvalidShapeError: {exc}") + except LinOSSError as exc: + print(f"Caught LinOSSError (unexpected subtype): {exc}") + + # 2. InvalidDampingError — negative damping coefficient G + section("InvalidDampingError: negative G") + y_ok = np.array([1.0]) + z_ok = np.array([0.0]) + A_ok = np.array([1.0]) + G_neg = np.array([-0.5]) # negative damping is physically meaningless + try: + linoss_step(y_ok, z_ok, A_ok, dt=0.1, damping=G_neg) + except InvalidDampingError as exc: + print(f"Caught InvalidDampingError: {exc}") + except LinOSSError as exc: + print(f"Caught LinOSSError (unexpected subtype): {exc}") + + # 3. UnsupportedModeError — unrecognised discretization scheme + section("UnsupportedModeError: unknown mode string") + try: + linoss_step(y_ok, z_ok, A_ok, dt=0.1, mode="euler") + except UnsupportedModeError as exc: + print(f"Caught UnsupportedModeError: {exc}") + except LinOSSError as exc: + print(f"Caught LinOSSError (unexpected subtype): {exc}") + + # 4. All three are subclasses of LinOSSError — catch them with a single clause. + section("Catching all errors with the base LinOSSError") + errors_to_try = [ + (y_bad, z_bad, A, 0.1, "implicit", None), + (y_ok, z_ok, A_ok, 0.1, "implicit", G_neg), + (y_ok, z_ok, A_ok, 0.1, "bogus_mode", None), + ] + for y_, z_, A_, dt_, mode_, G_ in errors_to_try: + try: + linoss_step(y_, z_, A_, dt=dt_, mode=mode_, damping=G_) + except LinOSSError as exc: + print(f" {type(exc).__name__}: {exc}") + + print("\nAll validation errors handled cleanly.") + + +if __name__ == "__main__": + main() diff --git a/examples/07_sunspots.py b/examples/07_sunspots.py new file mode 100644 index 0000000..e2e2933 --- /dev/null +++ b/examples/07_sunspots.py @@ -0,0 +1,102 @@ +"""07_sunspots.py — Bayesian oscillator fit on the sunspot time series. + +Fits a single damped oscillator SSM to the Wolfer sunspot number series +(annual data, ~11-year cycle) using maximum-likelihood via the Kalman filter. + +Requires the [probabilistic] extra: + pip install linoss-dynamics[probabilistic] + +If statsmodels is unavailable a synthetic dataset with the same structure +(period ~11 years, 200 observations) is used as a fallback. +""" + +from __future__ import annotations + +try: + from scipy import optimize # noqa: F401 +except ImportError: + print( + "This example requires scipy. " + "Install with: pip install linoss-dynamics[probabilistic]" + ) + raise SystemExit(0) + +import math + +import numpy as np + +from linoss_dynamics.fit import fit_oscillator_mle + + +def load_sunspots() -> np.ndarray: + """Return the Wolfer annual sunspot series (normalised), or a synthetic fallback.""" + try: + from statsmodels.datasets import sunspots + + data = sunspots.load_pandas().data + y_raw = data["SUNACTIVITY"].to_numpy(dtype=float) + print(f"Loaded statsmodels sunspot dataset: {len(y_raw)} annual observations") + return y_raw + except Exception: # noqa: BLE001 + print("statsmodels unavailable — using synthetic sunspot-like data (fallback)") + rng = np.random.default_rng(0) + t = np.arange(200, dtype=float) + # ~11-year cycle + noise, mean ~50, std ~30 (mimics Wolfer data statistics) + y = 50.0 + 30.0 * np.sin(2.0 * math.pi / 11.0 * t) + 10.0 * rng.standard_normal(200) + y = np.clip(y, 0.0, None) # sunspot counts are non-negative + return y + + +def main() -> None: + print("=== 07 Sunspots — Oscillator MLE Fit ===") + print() + + y_raw = load_sunspots() + + # Use the first 200 observations for speed; the rest can be used for evaluation. + y = y_raw[:200] + dt = 1.0 # annual data, dt = 1 year + + # Initial parameter guess informed by domain knowledge: + # beta ~ 0.05 (mild decay over decades) + # omega ~ 2π/11 (approx 11-year cycle) + # sigma_proc ~ 0.5, sigma_obs ~ 0.5 + omega_guess = 2.0 * math.pi / 11.0 + init = (0.05, omega_guess, 0.5, 0.5) + + print(f"Fitting {len(y)} annual observations, dt={dt} year") + print(f"Initial guess: beta={init[0]}, omega={init[1]:.4f} " + f"(period={2*math.pi/init[1]:.1f} yr), " + f"sigma_proc={init[2]}, sigma_obs={init[3]}") + print() + + result = fit_oscillator_mle(y, dt=dt, init=init) + + print("=== Fitted parameters ===") + print(f" beta = {result['beta']:.4f} (damping coefficient)") + print(f" omega = {result['omega']:.4f} rad/yr") + print(f" period = {result['period']:.2f} years (2π/omega)") + print(f" sigma_proc = {result['sigma_proc']:.4f}") + print(f" sigma_obs = {result['sigma_obs']:.4f}") + print(f" mean_y = {result['mean_y']:.2f} (subtracted before fitting)") + print() + + # The Kalman filter and RTS smoother outputs are in result["filt"] / result["smooth"]. + # m_f has shape (T, n) where n=2 (position, velocity state). + filt = result["filt"] + smooth = result["smooth"] + + print("First 10 filtered state means (position | velocity):") + print(f"{'t':>4} {'pos_filt':>12} {'vel_filt':>12} {'pos_smooth':>12}") + for t in range(10): + pos_f = filt["m_f"][t, 0] + result["mean_y"] # un-centre + vel_f = filt["m_f"][t, 1] + pos_s = smooth["m_s"][t, 0] + result["mean_y"] + print(f"{t:4d} {pos_f:12.3f} {vel_f:12.4f} {pos_s:12.3f}") + + print() + print(f"Log-likelihood of fitted model: {filt['loglik']:.2f}") + + +if __name__ == "__main__": + main() diff --git a/examples/08_phase_tracking.py b/examples/08_phase_tracking.py new file mode 100644 index 0000000..4a1b06c --- /dev/null +++ b/examples/08_phase_tracking.py @@ -0,0 +1,146 @@ +"""08_phase_tracking.py — Real-time phase estimation via Kalman filter. + +Simulates a single damped oscillator with additive observation noise, then +uses kalman_filter() to obtain filtered position/velocity estimates. From +those estimates we extract: + + - Instantaneous amplitude: sqrt(y_filt^2 + (z_filt / omega)^2) + - Instantaneous phase: arctan2(z_filt / omega, y_filt) + +and compare against ground truth. + +Requires the [probabilistic] extra: + pip install linoss-dynamics[probabilistic] +""" + +from __future__ import annotations + +try: + from scipy import optimize # noqa: F401 +except ImportError: + print( + "This example requires scipy. " + "Install with: pip install linoss-dynamics[probabilistic]" + ) + raise SystemExit(0) + +import math + +import numpy as np + +from linoss_dynamics.discretize import oscillator_mats +from linoss_dynamics.filters import kalman_filter + + +def simulate_oscillator( + omega: float, + beta: float, + sigma_proc: float, + sigma_obs: float, + T: int, + dt: float, + rng: np.random.Generator, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Simulate a noisy damped oscillator, return (y_obs, x_true, times). + + Returns: + y_obs: (T,) noisy scalar observations. + x_true: (T, 2) ground-truth [position, velocity] states. + times: (T,) time vector in seconds. + """ + Ad, Qd, _ = oscillator_mats(beta, omega, sigma_proc, dt) + # Cholesky for process noise sampling + L_q = np.linalg.cholesky(Qd + 1e-12 * np.eye(2)) + + x = np.array([1.0, 0.0]) # start at amplitude=1, phase=0 + x_true = np.empty((T, 2)) + y_obs = np.empty(T) + + for t in range(T): + x_true[t] = x + # Observation: position + Gaussian noise + y_obs[t] = x[0] + sigma_obs * rng.standard_normal() + # Propagate state + w = L_q @ rng.standard_normal(2) + x = Ad @ x + w + + times = np.arange(T) * dt + return y_obs, x_true, times + + +def main() -> None: + rng = np.random.default_rng(0) + + omega = 2.0 * math.pi / 11.0 # ~11-year period (sunspot-like) + beta = 0.02 # gentle damping + sigma_proc = 0.05 # process noise std + sigma_obs = 0.15 # observation noise std + T = 100 + dt = 1.0 # annual steps + + print("=== 08 Phase Tracking ===") + print(f"omega={omega:.4f} rad/step, period={2*math.pi/omega:.1f} steps") + print(f"beta={beta}, sigma_proc={sigma_proc}, sigma_obs={sigma_obs}") + print(f"T={T} steps, dt={dt}") + print() + + y_obs, x_true, times = simulate_oscillator( + omega, beta, sigma_proc, sigma_obs, T, dt, rng + ) + + # Build system matrices for the Kalman filter. + Ad, Qd, _ = oscillator_mats(beta, omega, sigma_proc, dt) + H = np.array([[1.0, 0.0]]) # observe position only + R = np.array([[sigma_obs**2]]) + + # Run the Kalman filter. + filt = kalman_filter(y_obs, Ad, Qd, H, R) + m_f = filt["m_f"] # shape (T, 2): [position, velocity] posterior means + + # Extract filtered position and velocity. + y_filt = m_f[:, 0] # filtered position estimate + z_filt = m_f[:, 1] # filtered velocity estimate + + # Convert to instantaneous amplitude and phase. + # For a simple harmonic oscillator: x(t) = A·cos(omega·t + phi) + # velocity z = dx/dt = -A·omega·sin(omega·t + phi) + # => z / omega ≈ -A·sin(...) + # => amplitude = sqrt(y^2 + (z/omega)^2) + # => phase = arctan2(-z/omega, y) [consistent with cos convention] + amp_filt = np.sqrt(y_filt**2 + (z_filt / omega) ** 2) + phase_filt = np.arctan2(-z_filt / omega, y_filt) + + # Ground-truth amplitude and phase from the true state. + y_true = x_true[:, 0] + z_true = x_true[:, 1] + amp_true = np.sqrt(y_true**2 + (z_true / omega) ** 2) + phase_true = np.arctan2(-z_true / omega, y_true) + + print("Phase estimates vs ground truth (first 10 steps):") + print(f"{'t':>4} {'y_obs':>8} {'y_filt':>8} {'phi_true':>10} {'phi_filt':>10} {'amp_filt':>10}") + for t in range(10): + print( + f"{t:4d} {y_obs[t]:8.3f} {y_filt[t]:8.3f} " + f"{phase_true[t]:10.4f} {phase_filt[t]:10.4f} {amp_filt[t]:10.4f}" + ) + + print() + print("Phase estimates vs ground truth (last 5 steps):") + print(f"{'t':>4} {'phi_true':>10} {'phi_filt':>10} {'amp_true':>10} {'amp_filt':>10}") + for t in range(T - 5, T): + print( + f"{t:4d} {phase_true[t]:10.4f} {phase_filt[t]:10.4f} " + f"{amp_true[t]:10.4f} {amp_filt[t]:10.4f}" + ) + + # Summary statistics + phase_rmse = float(np.sqrt(np.mean((phase_filt - phase_true) ** 2))) + amp_rmse = float(np.sqrt(np.mean((amp_filt - amp_true) ** 2))) + print() + print(f"Phase RMSE: {phase_rmse:.4f} rad") + print(f"Amplitude RMSE: {amp_rmse:.4f}") + print(f"Filter log-likelihood: {filt['loglik']:.2f}") + + +if __name__ == "__main__": + main() diff --git a/examples/README.md b/examples/README.md new file mode 100644 index 0000000..0378cd0 --- /dev/null +++ b/examples/README.md @@ -0,0 +1,52 @@ +# linoss-dynamics Examples + +Runnable tutorials for `linoss-dynamics`. Each script is self-contained and +prints labeled output so you can follow along without reading source code. + +Run any example from the repository root: + +```bash +python examples/