"""Weak Koopman latent-dynamics toy verification. Run with: python3 weak_koopman_experiment.py Outputs results.json with analytic predictions and measured values. """ import json from pathlib import Path import numpy as np SEED = 1425 rng = np.random.default_rng(SEED) def psi_and_derivative(t, a, b): """Compactly supported Hann test function and its physical derivative.""" u = (t - a) / (b - a) inside = (u >= 0) & (u <= 1) psi = np.zeros_like(t, dtype=float) dpsi = np.zeros_like(t, dtype=float) psi[inside] = np.sin(np.pi * u[inside]) ** 2 dpsi[inside] = (np.pi / (b - a)) * np.sin(2 * np.pi * u[inside]) return psi, dpsi def weak_weights(n, dt, width): t = np.arange(n) * dt psi, dpsi = psi_and_derivative(t, 0.0, width * dt) # B = - integral psi'(t) z(t) dt; sign is irrelevant for variance. return psi, -dt * dpsi def analytic_variances(sigma, dt, width): _, wprime = weak_weights(width, dt, width) weak = sigma**2 * np.sum(wprime**2) fd = 2.0 * sigma**2 / dt**2 return weak, fd def identity_check(): # Smooth oscillator observable z=[sin(t), cos(t)] and exact generator. dt, n = 0.002, 1501 t = np.arange(n) * dt z = np.column_stack([np.sin(t), np.cos(t)]) A = np.array([[0.0, 1.0], [-1.0, 0.0]]) psi, wp = weak_weights(n, dt, n - 1) G = np.sum((dt * psi)[:, None] * z, axis=0) B = np.sum(wp[:, None] * z, axis=0) # B_l = integral psi * (A z)_l = (G @ A.T)_l. rhs = G @ A.T return float(np.max(np.abs(B - rhs)) / (np.max(np.abs(rhs)) + 1e-12)) def variance_sweep(): # Prediction 1: empirical weak noise variance / analytic variance ~= 1, # and FD variance / analytic variance ~= 1, over sigma. dt, width, reps = 0.01, 2.0, 12000 sigmas = [0.02, 0.05, 0.10, 0.20] _, ww = weak_weights(int(round(width / dt)) + 1, dt, int(round(width / dt))) # use exactly the samples represented by the returned weights rows_sigma = [] for sigma in sigmas: weak_samples = sigma * (rng.standard_normal((reps, len(ww))) @ ww) fd_samples = sigma * (rng.standard_normal((reps, 2))) @ np.array([-1.0, 1.0]) / dt weak_emp = float(np.var(weak_samples, ddof=1)) fd_emp = float(np.var(fd_samples, ddof=1)) weak_pred = float(sigma**2 * np.sum(ww**2)) fd_pred = float(2 * sigma**2 / dt**2) rows_sigma.append({"sigma": sigma, "weak_emp": weak_emp, "weak_pred": weak_pred, "fd_emp": fd_emp, "fd_pred": fd_pred, "weak_ratio": weak_emp / weak_pred, "fd_ratio": fd_emp / fd_pred}) # Prediction 2: at fixed physical window, weak variance is proportional to # sigma^2 and approximately dt, while pointwise FD is proportional dt^-2. # Prediction 3: with fixed dt, widening a Hann test function gives W^-3. rows_dt, rows_width = [], [] sigma, physical_width = 0.10, 2.0 for dt_i in [0.005, 0.01, 0.02, 0.04, 0.08]: m = int(round(physical_width / dt_i)) weak_pred, fd_pred = analytic_variances(sigma, dt_i, m) _, ww_i = weak_weights(m + 1, dt_i, m) weak_emp = float(np.var(sigma * (rng.standard_normal((reps, len(ww_i))) @ ww_i), ddof=1)) fd_emp = float(np.var(sigma * (rng.standard_normal((reps, 2))) @ np.array([-1., 1.]) / dt_i, ddof=1)) rows_dt.append({"dt": dt_i, "weak_emp": weak_emp, "weak_pred": weak_pred, "fd_emp": fd_emp, "fd_pred": fd_pred}) for m in [50, 75, 100, 150, 200, 300, 400]: dt_i = 0.01 weak_pred, fd_pred = analytic_variances(sigma, dt_i, m) _, ww_i = weak_weights(m + 1, dt_i, m) weak_emp = float(np.var(sigma * (rng.standard_normal((reps, len(ww_i))) @ ww_i), ddof=1)) rows_width.append({"width": m * dt_i, "weak_emp": weak_emp, "weak_pred": weak_pred, "fd_pred": fd_pred}) # Log-log slopes quantify the parameter scalings. dt_arr = np.array([r["dt"] for r in rows_dt]) weak_dt_slope = float(np.polyfit(np.log(dt_arr), np.log([r["weak_emp"] for r in rows_dt]), 1)[0]) fd_dt_slope = float(np.polyfit(np.log(dt_arr), np.log([r["fd_emp"] for r in rows_dt]), 1)[0]) w_arr = np.array([r["width"] for r in rows_width]) weak_w_slope = float(np.polyfit(np.log(w_arr), np.log([r["weak_emp"] for r in rows_width]), 1)[0]) return {"sigma": rows_sigma, "dt": rows_dt, "width": rows_width, "slopes": {"weak_vs_dt": weak_dt_slope, "fd_vs_dt": fd_dt_slope, "weak_vs_window": weak_w_slope}} def generator_mini_experiment(): """Compare FD and weak regression for the exact 2D oscillator generator.""" dt, n, trajectories, sigma = 0.02, 101, 300, 0.15 t = np.arange(n) * dt A_true = np.array([[0., 1.], [-1., 0.]]) errs_fd, errs_weak = [], [] width = (n - 1) * dt psi, wp = weak_weights(n, dt, n - 1) # Multiple random phase trajectories make a small EDMD-like regression. for _ in range(trajectories): phase = rng.uniform(0, 2*np.pi) clean = np.column_stack([np.sin(t + phase), np.cos(t + phase)]) noisy = clean + sigma * rng.standard_normal(clean.shape) # Standard pointwise derivative matching, central differences. dz = (noisy[2:] - noisy[:-2]) / (2*dt) Zmid = noisy[1:-1] Afd_T = np.linalg.lstsq(Zmid, dz, rcond=None)[0] Afd = Afd_T.T # One weak equation per trajectory is insufficient for 2D A, so use # several overlapping windows as rows in G and B. Gs, Bs = [], [] win = 51 for start in range(0, n - win + 1, 10): zwin = noisy[start:start+win] ps, ws = weak_weights(win, dt, win - 1) Gs.append(np.sum((dt * ps)[:, None] * zwin, axis=0)) Bs.append(np.sum(ws[:, None] * zwin, axis=0)) G = np.asarray(Gs); B = np.asarray(Bs) Aw_T = np.linalg.solve(G.T @ G + 1e-5*np.eye(2), G.T @ B) Aw = Aw_T.T errs_fd.append(np.linalg.norm(Afd - A_true)) errs_weak.append(np.linalg.norm(Aw - A_true)) return {"sigma": sigma, "fd_generator_error": float(np.mean(errs_fd)), "weak_generator_error": float(np.mean(errs_weak)), "weak_over_fd": float(np.mean(errs_weak) / np.mean(errs_fd))} def main(): result = {"seed": SEED, "identity_relative_error": identity_check(), "variance_sweeps": variance_sweep(), "generator_mini_experiment": generator_mini_experiment()} Path("results.json").write_text(json.dumps(result, indent=2)) print(json.dumps(result, indent=2)) if __name__ == "__main__": main()