Metadata-Version: 2.3
Name: pywynda
Version: 0.1.0
Summary: A zero-optimization, online streaming adaptive observer for nonlinear dynamical systems
Author: Arpit Singh
Author-email: Arpit Singh <iarpitsingh.bms@gmail.com>
License: MIT
Classifier: Programming Language :: Python :: 3
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: OS Independent
Classifier: Topic :: Scientific/Engineering :: Mathematics
Requires-Dist: numpy>=2.5.0
Requires-Dist: scipy>=1.18.0
Requires-Dist: matplotlib>=3.8.0
Description-Content-Type: text/markdown

﻿<div align="center">

<h1>
  <img src="https://img.shields.io/badge/pywynda-v0.1.0-4A90E2?style=for-the-badge&labelColor=1a1a2e" alt="pywynda"/>
</h1>

**Wide-Array of Nonlinear Dynamics Approximation**

*A zero-optimization, online streaming adaptive observer for real-time identification of nonlinear dynamical systems with time-varying parameters.*

[![Python](https://img.shields.io/badge/Python-3.13%2B-blue?logo=python&logoColor=white)](https://www.python.org/)
[![NumPy](https://img.shields.io/badge/NumPy-2.5%2B-013243?logo=numpy&logoColor=white)](https://numpy.org/)
[![SciPy](https://img.shields.io/badge/SciPy-1.18%2B-8CAAE6?logo=scipy&logoColor=white)](https://scipy.org/)
[![Tests](https://img.shields.io/badge/pytest-26%20passed-brightgreen?logo=pytest&logoColor=white)](tests/)
[![License](https://img.shields.io/badge/License-MIT-green)](LICENSE)
[![Reference](https://img.shields.io/badge/DOI-10.1016%2Fj.mex.2024.102625-orange)](https://doi.org/10.1016/j.mex.2024.102625)

</div>

---

## Overview

**pywynda** is a Python implementation of the WyNDA (Wide-Array of Nonlinear Dynamics Approximation) framework, a simultaneous state-and-parameter observer derived from the theory published in:

> Hasan, A. (2024). *WyNDA: A Method to Discover Mathematical Models of Nonlinear Dynamical Systems.* **MethodsX**, 12, 102625. https://doi.org/10.1016/j.mex.2024.102625

The framework solves the **nonlinear system identification problem** in a strictly online, recursive fashion — consuming one measurement sample at a time with no batch history, no gradient descent, and no optimization solver of any kind. Its core engine is a **dual-gain Kalman-type recursion** that simultaneously corrects state estimates and parameter estimates at every time step, with exponential forgetting factors that make the observer adaptive to structural changes in the system dynamics.

---

## Three Core Design Constraints

### 1. Zero-Optimization Constraint

> **No `scikit-learn`. No `scipy.optimize`. No neural network training loop. No gradient.**

Every matrix operation reduces to a single `numpy.linalg.solve` call per time step. The observer is a pure **matrix-algebra recursion** — not a regression, not a minimization, and not a learned model. This makes it:

- Deterministic and reproducible without random seeds
- O(n³) per step in parameter dimension, with no iterative convergence
- Deployable on hardware that cannot run an optimization runtime

### 2. Strict Online Streaming

> **One sample in → corrected state + updated parameters out.**

The `.step(y, Psi)` method maintains all internal state between calls. There is no buffer, no window, no queue. The observer is fully compatible with real-time sensor streams, embedded control loops, and hardware-in-the-loop pipelines.

### 3. Non-Stationary Adaptation Edge

> **Exponential forgetting discounts stale history. No retraining required.**

Two scalar forgetting factors — $\lambda_x \in (0,1)$ for state estimation and $\lambda_\theta \in (0,1)$ for parameter estimation — inflate the respective covariance matrices at each step, ensuring that recent data always dominates the gain computation. When system parameters change abruptly (component failure, structural drift, regime transition), the observer automatically re-identifies the new dynamics within a time window of approximately $\tau \approx \Delta t \,/\,(1 - \lambda_\theta)$.

---

## Mathematical Engine

The observer solves the following discrete-time system identification problem at each step $k$.

### System Model

| Equation | Role |
|---|---|
| $y(k) = x(k) + v(k)$ | Measurement with noise $v(k)$ |
| $x(k+1) = x(k) + \Psi(y(k), u(k))\,\theta$ | Euler-discretised nonlinear dynamics |

where $\Psi \in \mathbb{R}^{n \times p}$ is the **block-diagonal regressor library** evaluated at the current measurement $y(k)$ and control input $u(k)$, and $\theta \in \mathbb{R}^p$ is the unknown parameter vector.

### Gain Computations

$$K_{x}(k) = P_{x}(k|k-1)\bigl(P_{x}(k|k-1) + R_{x}(k)\bigr)^{-1}$$

$$\Omega(k) = \Gamma(k|k-1)\,P_{\theta}(k|k-1)\,\Gamma^{\top}(k|k-1) + R_{\theta}(k)$$

$$K_{\theta}(k) = P_{\theta}(k|k-1)\,\Gamma^{\top}(k|k-1)\,\Omega(k)^{-1}$$

where $\Gamma(k|k-1)$ is the coupling matrix that links the parameter subspace to the innovation signal, and $R_x$, $R_\theta$ are measurement noise covariances acting as natural Tikhonov-style regularisers.

### State and Parameter Correction

$$\bar{x}(k|k) = \bar{x}(k|k-1) + \bigl(K_{x}(k) + \Gamma(k|k-1)\,K_{\theta}(k)\bigr)\bigl(y(k) - \bar{x}(k|k-1)\bigr)$$

$$\bar{\theta}(k|k) = \bar{\theta}(k|k-1) - K_{\theta}(k)\bigl(y(k) - \bar{x}(k|k-1)\bigr)$$

### One-Step-Ahead Predictions and Covariance Transitions

$$\bar{x}(k+1|k) = \bar{x}(k|k) + \Psi(y(k), u(k))\,\bar{\theta}(k|k)$$

$$P_{x}(k+1|k) = \lambda_{x}^{-1}\bigl(I - K_{x}(k)\bigr)P_{x}(k|k-1)$$

$$P_{\theta}(k+1|k) = \lambda_{\theta}^{-1}\bigl(I - K_{\theta}(k)\,\Gamma(k|k-1)\bigr)P_{\theta}(k|k-1)$$

$$\Gamma(k+1|k) = \bigl(I - K_{x}(k)\bigr)\Gamma(k|k-1) - \Psi(y(k), u(k))$$

The $\lambda^{-1}$ inflation in the covariance transitions is the mechanism of **exponential forgetting** — each step slightly enlarges $P_x$ and $P_\theta$, increasing the observer gain and responsiveness to new data.

### Theorem 1 (MethodsX 2024): Exponential Error Decay

Under the **Persistency of Excitation (PE)** condition — that the regressor $\Psi$ spans the full parameter space over a finite time window — both error norms decay exponentially:

$$\|\varepsilon_x(k)\| \leq C_x\,e^{-\alpha_x k}, \qquad \|\varepsilon_\theta(k)\| \leq C_\theta\,e^{-\alpha_\theta k}$$

This guarantee is validated empirically across all five benchmark examples in this repository.

---

## Repository Layout

```
pywynda/
├── src/
│   └── pywynda/
│       ├── __init__.py          # Public API surface
│       ├── core.py              # WyNDAObserver — the recursive gain engine
│       ├── utils.py             # Trajectory generation, PE checks, error metrics
│       └── py.typed             # PEP 561 marker (fully typed)
│
├── tests/
│   └── test_observer.py         # 26 unit tests across 8 validation categories
│
├── examples/
│   ├── ex1_msd.py               # Linear MSD — over-parameterised library
│   ├── ex2_nonlinear.py         # Lorenz, Rossler, Lotka-Volterra, Van der Pol
│   ├── ex3_control_pe.py        # PE vs. no-PE comparison on forced oscillator
│   ├── ex4_maglev.py            # Magnetic levitation — 3-state physical system
│   ├── ex5_sindy_comparison.py  # WyNDA vs. Batch SINDy: non-stationary drift
│   ├── run_all.py               # Sequential manifest runner for all examples
│   └── plots/                   # Auto-generated PNG convergence/comparison plots
│
├── pyproject.toml
├── uv.lock
└── README.md
```

---

## Installation

pywynda uses the [`uv`](https://docs.astral.sh/uv/) toolchain for environment management.

```bash
# Clone the repository
git clone https://github.com/scarwizz/pywynda.git
cd pywynda

# Create and activate the virtual environment
uv venv
# Windows:
.venv\Scripts\activate
# Linux / macOS:
source .venv/bin/activate

# Install all dependencies (numpy, scipy, matplotlib, pytest)
uv sync

# Verify installation with the full test suite
uv run pytest tests/ -v
```

> **Runtime dependencies:** `numpy >= 2.5`, `scipy >= 1.18`, `matplotlib >= 3.8`
> **No** `scikit-learn`, `torch`, `tensorflow`, `jax`, or any optimization framework.

---

## Quickstart

```python
import numpy as np
from pywynda import WyNDAObserver

# ── System dimensions ──────────────────────────────────────────────────────
n_states = 2   # [position, velocity]
n_params = 6   # 3 basis functions per state equation (block-diagonal library)
dt       = 0.001

# ── Initialise the observer ────────────────────────────────────────────────
obs = WyNDAObserver(
    n_states   = n_states,
    n_params   = n_params,
    lambda_x   = 0.999,          # state-estimation forgetting factor
    lambda_theta = 0.998,        # parameter-estimation forgetting factor (tau ~ 0.5 s)
    P_x0       = 100.0 * np.eye(n_states),   # initial state covariance
    P_theta0   = 500.0 * np.eye(n_params),   # initial parameter covariance
    R_x        = 0.05**2 * np.eye(n_states), # measurement noise covariance
    R_theta    = 5e-3  * np.eye(n_states),   # parameter noise covariance (n x n)
    x0_hat     = np.zeros(n_states),         # initial state estimate
    theta0_hat = np.zeros(n_params),         # initial parameter guess
)

# ── Build the block-diagonal regressor Psi (n x p) at measurement y ───────
def build_psi(y: np.ndarray) -> np.ndarray:
    x1, x2 = y[0], y[1]
    phi = dt * np.array([x1, x2, x1**2])     # 3-entry basis per equation
    Psi = np.zeros((n_states, n_params))
    Psi[0, :3] = phi                          # x1 equation
    Psi[1, 3:] = phi                          # x2 equation
    return Psi

# ── Streaming loop — one sample at a time ─────────────────────────────────
# (replace with your real sensor stream)
for k, y_k in enumerate(your_measurement_sequence):
    state = obs.step(y=y_k, Psi=build_psi(y_k))

    x_hat     = state.x_hat      # corrected state estimate, shape (n,)
    theta_hat = state.theta_hat  # updated parameter vector, shape (p,)
    inno      = state.innovation # innovation signal y(k) - x_hat(k|k-1)

# ── Known-physics feed-forward (optional) ─────────────────────────────────
# Inject analytically-known contributions (e.g. gravity, kinematic coupling)
# without contaminating the parametric subspace:
obs.inject_feedforward(delta_x=np.array([dt * float(y_k[1]), dt * 9.81]))
```

---

## Benchmark Examples

Run all five examples sequentially and regenerate all plots:

```bash
uv run python examples/run_all.py
```

Or run individual examples:

```bash
uv run python examples/ex1_msd.py          # ~3 s
uv run python examples/ex2_nonlinear.py    # ~25 s (4 chaotic systems)
uv run python examples/ex3_control_pe.py   # ~5 s
uv run python examples/ex4_maglev.py       # ~3 s
uv run python examples/ex5_sindy_comparison.py  # ~10 s
```

### Example Suite Summary

| Example | System | Key Result |
|---|---|---|
| **Ex 1** | Mass-Spring-Damper (over-parameterised, 12 basis) | $k/m$ and $b/m$ recovered; spurious coefficients collapse toward zero |
| **Ex 2a** | Lorenz strange attractor | $\hat{\sigma}=9.78$, $\hat{\rho}=27.79$, $\hat{\beta}=2.98$; $\varepsilon_\theta = 0.53$ at $t=2$ s |
| **Ex 2b** | Rössler system | Expected PE-deficit failure on $c$-parameter — documented in §4.2 of reference |
| **Ex 2c** | Lotka-Volterra predator-prey | All 4 parameters within 5 % of truth; $\varepsilon_\theta(\text{final}) = 0.034$ |
| **Ex 2d** | Van der Pol oscillator | $\hat{\mu} = 1.196$ (0.3 % error); dual-coefficient consistency confirmed |
| **Ex 3** | Forced oscillator — PE vs. no-PE | Empirical proof of Theorem 1 necessity condition |
| **Ex 4** | Magnetic levitation (3-state) | $\hat{L}=3.002$ H, $\hat{R}=6.004\ \Omega$, $\hat{C}=1.497$ F (sub-1 % error in 0.2 s) |
| **Ex 5** | MSD with sudden parameter drift | **WyNDA vs. Batch SINDy** — see dedicated section below |

---

## Benchmark Showcase: WyNDA vs. Batch SINDy under Non-Stationary Dynamics

Example 5 (`ex5_sindy_comparison.py`) is a rigorous head-to-head evaluation on the most challenging scenario in data-driven system identification: **sudden parameter drift** mid-operation.

### Experiment Design

A Mass-Spring-Damper system is subjected to a **hard step change in spring stiffness** at $t = 2.5$ s, simulating a sudden structural failure or component degradation:

$$k(t) = \begin{cases} 84.0 \text{ N/m} & 0.0 \leq t < 2.5 \text{ s} \\ 40.0 \text{ N/m} & 2.5 \leq t \leq 5.0 \text{ s} \end{cases}$$

Multi-frequency PE forcing ($3\sin(2\pi \cdot 1.5\, t) + 2\sin(2\pi \cdot 3.7\, t)$) ensures the PE lower bound is satisfied throughout. White Gaussian noise ($\sigma = 0.05$) is applied to all state measurements.

Both algorithms run on the **exact same trajectory** (identical random seed):

| Parameter | WyNDA | Batch SINDy |
|---|---|---|
| Processing mode | Online streaming (1 sample at a time) | Batch (re-fit every 200 steps over all history) |
| Forgetting | $\lambda_\theta = 0.998$ ($\tau \approx 0.5$ s) | None — accumulates all data from $t=0$ |
| Solver | `np.linalg.solve` (one call/step) | `np.linalg.solve` on $\Theta^\top\Theta + \alpha I$ |
| Sklearn / optimize | **No** | **No** (zero-optimization compliant) |

### Performance Telemetry

| Metric | WyNDA | Batch SINDy | Advantage |
|---|---|---|---|
| Full-run $k$ RMSE | 20.4 N/m | 35.3 N/m | WyNDA **1.7×** lower |
| **Post-drift $k$ RMSE** | **23.3 N/m** | **41.3 N/m** | **WyNDA 1.8× lower** |
| $\hat{k}$ at $t = 5.0$ s | **41.5 N/m** ✓ | **79.4 N/m** ✗ | WyNDA converges; SINDy does not |
| Detection latency (5 N/m tol.) | **1040 ms** | **> T_end (never)** | SINDy structurally fails |

### Why SINDy Fails Here

SINDy's batch regression accumulates all historical data $[0 \ldots k]$ into the design matrix $\Theta$. After the drift at $t = 2.5$ s, the regression is contaminated by 2 500 steps of pre-drift dynamics where $k = 84$ N/m. The normal equations mix two incompatible physical regimes and converge to a **weighted mean** of the two spring constants — stuck near 79 N/m at end-of-run — regardless of how many post-drift samples arrive.

WyNDA's forgetting factor $\lambda_\theta = 0.998$ means data from $3\tau = 1.5$ s ago carries only $\lambda_\theta^{1500} \approx 5\%$ weight in the gain computation. The observer effectively erases the pre-drift regime within two forgetting time-constants and re-identifies $k = 40$ N/m cleanly.

> **This demonstrates WyNDA's definitive architectural edge:** global batch regression is structurally incompatible with non-stationary dynamical systems. Online streaming with forgetting is not a heuristic optimisation — it is the mathematically correct formulation for this class of problems.

---

## Test Suite

The full validation suite confirms all theoretical guarantees from MethodsX 12 (2024) 102625.

```
========================= test session starts ==========================
platform win32 -- Python 3.13.11, pytest-9.1.1, pluggy-1.6.0
rootdir: C:\...\pywynda
configfile: pyproject.toml
collected 26 items

tests/test_observer.py::TestOutputShapes::test_observer_state_shapes         PASSED [  3%]
tests/test_observer.py::TestOutputShapes::test_observer_history_shapes        PASSED [  7%]
tests/test_observer.py::TestOutputShapes::test_psi_sequence_shape             PASSED [ 11%]
tests/test_observer.py::TestConvergence::test_state_estimation_convergence    PASSED [ 15%]
tests/test_observer.py::TestConvergence::test_parameter_estimation_convergence PASSED [ 19%]
tests/test_observer.py::TestConvergence::test_final_parameter_accuracy         PASSED [ 23%]
tests/test_observer.py::TestConvergence::test_state_error_last_quarter_small  PASSED [ 26%]
tests/test_observer.py::TestConvergence::test_exponential_decay_envelope      PASSED [ 30%]
tests/test_observer.py::TestPersistenceOfExcitation::test_chirp_satisfies_pe  PASSED [ 34%]
tests/test_observer.py::TestPersistenceOfExcitation::test_white_noise_satisfies_pe PASSED [ 38%]
tests/test_observer.py::TestPersistenceOfExcitation::test_gram_matrix_positive_definite PASSED [ 42%]
tests/test_observer.py::TestDeterminism::test_same_seed_same_results          PASSED [ 46%]
tests/test_observer.py::TestDeterminism::test_different_seeds_different_results PASSED [ 50%]
tests/test_observer.py::TestForgettingFactor::test_lower_lambda_faster_initial_theta_convergence PASSED [ 53%]
tests/test_observer.py::TestPredictionIntegrity::test_x_pred_finite           PASSED [ 57%]
tests/test_observer.py::TestPredictionIntegrity::test_theta_pred_finite       PASSED [ 61%]
tests/test_observer.py::TestPredictionIntegrity::test_innovations_bounded     PASSED [ 65%]
tests/test_observer.py::TestPredictionIntegrity::test_covariances_decrease    PASSED [ 69%]
tests/test_observer.py::TestInputValidation::test_invalid_lambda_x_zero       PASSED [ 73%]
tests/test_observer.py::TestInputValidation::test_invalid_lambda_x_too_large  PASSED [ 76%]
tests/test_observer.py::TestInputValidation::test_invalid_lambda_theta_negative PASSED [ 80%]
tests/test_observer.py::TestInputValidation::test_wrong_P_x_shape             PASSED [ 84%]
tests/test_observer.py::TestInputValidation::test_run_mismatched_steps        PASSED [ 88%]
tests/test_observer.py::TestRunAPI::test_run_produces_all_fields              PASSED [ 92%]
tests/test_observer.py::TestRunAPI::test_run_step_consistency                 PASSED [ 96%]
tests/test_observer.py::TestRunAPI::test_multiple_input_types_converge        PASSED [100%]

========================= 26 passed in 4.77s ===========================
```

### Test Category Coverage

| # | Category | Scenarios | Validates |
|---|---|---|---|
| 1 | **Shape Tests** | 3 | All output matrices carry correct dimensions |
| 2 | **Convergence Tests** | 5 | Exponential decay of $\varepsilon_x$ and $\varepsilon_\theta$ (Theorem 1) |
| 3 | **Persistency of Excitation** | 3 | PE lower bound satisfied for chirp and white-noise inputs |
| 4 | **Determinism** | 2 | Same seed → bit-identical trajectories |
| 5 | **Forgetting Factor Sensitivity** | 1 | Lower $\lambda$ gives faster initial convergence |
| 6 | **Prediction Integrity** | 4 | $\hat{x}(k+1|k)$ stays finite; covariances decrease |
| 7 | **Input Validation** | 5 | Constructor and `.step()` raise on bad arguments |
| 8 | **Run() API** | 3 | `ObserverHistory` consistency with step-by-step calls |

---

## Public API Reference

### `WyNDAObserver`

```python
class WyNDAObserver:
    def __init__(
        self,
        n_states:     int,
        n_params:     int,
        lambda_x:     float,          # state forgetting factor, in (0, 1)
        lambda_theta: float,          # parameter forgetting factor, in (0, 1)
        P_x0:         NDArray,        # initial state covariance, shape (n, n)
        P_theta0:     NDArray,        # initial parameter covariance, shape (p, p)
        R_x:          NDArray,        # measurement noise covariance, shape (n, n)
        R_theta:      NDArray,        # parameter noise covariance, shape (n, n)
        x0_hat:       NDArray,        # initial state estimate, shape (n,)
        theta0_hat:   NDArray,        # initial parameter estimate, shape (p,)
    ) -> None: ...

    def step(self, y: NDArray, Psi: NDArray) -> ObserverState:
        """Process one measurement. Returns the corrected state and parameters."""

    def inject_feedforward(self, delta_x: NDArray) -> None:
        """Add a known non-parametric dynamics correction to the one-step-ahead prior.
        Use for gravity, kinematic coupling, or other analytically-known contributions."""

    def run(self, Y: NDArray, Psi_seq: NDArray) -> ObserverHistory:
        """Process a full measurement sequence and return stacked history."""
```

### `ObserverState` (returned by `.step()`)

| Field | Shape | Description |
|---|---|---|
| `x_hat` | `(n,)` | Corrected state estimate $\bar{x}(k|k)$ |
| `theta_hat` | `(p,)` | Updated parameter vector $\bar{\theta}(k|k)$ |
| `x_pred` | `(n,)` | One-step-ahead prediction $\bar{x}(k+1|k)$ |
| `theta_pred` | `(p,)` | One-step-ahead parameter prediction |
| `P_x` | `(n,n)` | State covariance $P_x(k+1|k)$ |
| `P_theta` | `(p,p)` | Parameter covariance $P_\theta(k+1|k)$ |
| `Gamma` | `(n,p)` | Coupling matrix $\Gamma(k+1|k)$ |
| `K_x` | `(n,n)` | State Kalman gain $K_x(k)$ |
| `K_theta` | `(p,n)` | Parameter Kalman gain $K_\theta(k)$ |
| `innovation` | `(n,)` | Innovation signal $y(k) - \bar{x}(k|k-1)$ |

---

## Design Principles

### Regressor Library `Psi` — Block-Diagonal Structure

The regressor $\Psi \in \mathbb{R}^{n \times p}$ must be **block-diagonal** to decouple each state equation's parameter subspace. For an $n$-state system with $m$ basis functions per equation:

```
Psi = [ phi(y, u) |   0     ]   Row 0: x1-equation parameters
      [   0       | phi(y,u)]   Row 1: x2-equation parameters
                   p = n * m
```

This structure ensures the parameter correction $K_\theta \cdot \text{innovation}$ updates only the relevant equation coefficients at each step.

### Forgetting Factor Selection Guide

| Regime | $\lambda_\theta$ recommendation | Effective memory window |
|---|---|---|
| Stationary system | $0.999$ | ~1 000 steps |
| Slow drift | $0.998$ | ~500 steps |
| Sudden step-change | $0.995$ – $0.998$ | 200–500 steps |
| Rapid adaptation needed | $0.990$ | ~100 steps |

> Setting $\lambda_\theta$ too low causes **covariance wind-up** during stable phases; too high causes **slow drift tracking**. The sweet spot for step-change detection is $\lambda_\theta \approx 0.998$, giving a forgetting time-constant $\tau = \Delta t / (1 - \lambda_\theta)$.

### Solver Compliance

Every linear system in the recursion is solved via `numpy.linalg.solve`, never `numpy.linalg.inv`. This follows numerical best practice: solving $A x = b$ directly via LU factorisation is both faster and numerically stabler than computing $A^{-1} b$.

---

## Reference

```bibtex
@article{hasan2024wynda,
  title   = {{WyNDA}: A Method to Discover Mathematical Models of
              Nonlinear Dynamical Systems},
  author  = {Hasan, Agus},
  journal = {MethodsX},
  volume  = {12},
  pages   = {102625},
  year    = {2024},
  doi     = {10.1016/j.mex.2024.102625},
  url     = {https://doi.org/10.1016/j.mex.2024.102625}
}
```

---

## License

MIT License — see [`LICENSE`](LICENSE) for details.

---

<div align="center">
<sub>Built with NumPy. Zero optimization dependencies. Strictly online. Reference implementation of MethodsX 12 (2024) 102625.</sub>
</div>
