diff --git a/.gitignore b/.gitignore deleted file mode 100644 index 0f54b01..0000000 --- a/.gitignore +++ /dev/null @@ -1,71 +0,0 @@ -# Universal .gitignore for Drizzy's Research Ecosystem - -# Rust / C++ Build Artifacts -target/ -**/target/ -build/ -bin/ -*.exe -*.dll -*.so -*.dylib -*.pdb -*.rlib -*.rmeta - -# Python & AI Metadata -__pycache__/ -*.py[cod] -*$py.class -venv/ -.env -.venv -env/ -.claude/ -.gemini/ -.agents/ -.qwen/ -.anthropic/ - -# IDEs & System -.vscode/ -.idea/ -.gradle/ -.kotlin/ -*.iml -.DS_Store -Thumbs.db -desktop.ini - -# Data & Logs — code is published, data stays local (the rule) -*.log -data/*.csv -data/*.json -data/*.sqlite -data/*.sqlite-shm -data/*.sqlite-wal -!data/.gitkeep -out/ -dist/ - -# Heavy data artifacts anywhere (datasets, processed outputs, figures from runs). -# These are regenerated locally from code; never bloat the repo with them. -**/data/raw/ -**/data/processed/ -*.parquet -*.npz -*.npy -*.h5 -*.hdf5 -# Generated run figures (keep hand-made diagrams in docs/ or assets/ instead) -**/data/**/*.png - -# Benchmark results (large JSON/npz outputs) -modules/cryptotn_gpu/benchmarks/results/ - -# Private / local-only — never publish -CLAUDE.md -scratch/ -modules/simulations/ -modules/pattern_analysis/ -*.tmp diff --git a/README.md b/README.md deleted file mode 100644 index 6e63db2..0000000 --- a/README.md +++ /dev/null @@ -1,75 +0,0 @@ -# SUBSTRATE - -> *From quantum vacuum to geomagnetic civilization risk — a multi-scale computational framework.* - -SUBSTRATE is a unified research platform analyzing electromagnetic fields across all physical scales: -from quantum field theory and quantum biological sensing, through geomagnetic dynamics and solar physics, -to cosmological structure. - -## The thesis - -At every scale, the same question: **what is the substrate generating the observable pattern?** - -- Quantum vacuum fluctuations → fields → forces -- Radical pair coherence → biological magnetoreception -- Geomagnetic field → civilizational protection layer -- Heliospheric dynamics → cosmic ray modulation -- CMB fluctuations → large-scale structure -- Neural state (EEG) → physical parameter modulation ← KHAOS bridge - -## Modules - -| Module | Status | Description | -|--------|--------|-------------| -| `quantum_lab` | ✅ Operational | Lattice gauge theory, tensor networks, hand-written U(1) plaquette CUDA kernel (roofline crossover, kernel-only up to 154× vs JAX-CPU) | -| `cryptotn_gpu` | ✅ Operational | Radical pair quantum biology, χ tensor engine | -| `magnon` | ✅ Operational | Avian magnetoreception, Lindblad dynamics | -| `cycle_project` | ✅ Operational | Geomagnetic field: detection, simulation, RAG, monitor, forecast | -| `simulations` | ✅ Integrated | Multi-domain physics simulations — see below | -| `pattern_analysis` | ✅ Integrated | Rust Welch FFT + cross-modal PNG/audio/JSON analytics | -| `heliospheric` | 📋 Planned | Solar dynamo, Be-10 proxy, heliosphere coupling | -| `cosmological` | 📋 Planned | CMB anomaly detection (Planck data, GNN) | - -### `simulations` — submodules - -| Submodule | Scripts | Description | -|-----------|---------|-------------| -| `cosmology` | 28 | Alcubierre warp drive, Kerr-Newman geodesics, big bang/BBN/CMB, Hawking evaporation, wormholes, Gödel CTC, Taub-NUT, ER=EPR, big rip, Unruh effect | -| `quantum` | 13 | 3-body quantum chaos (RK4 + Euler-Maruyama), Wheeler delayed choice, quantum Zeno, instanton tunneling, QCD confinement, LQG spin networks, Orch-OR microtubules, Landauer-Bekenstein, fractal dimension | -| `astrophysics` | 14 | Real TESS/NASA transit data (TOI-1444, TIC98796344, TIC400799224), exoplanet hunting, JWST, biosignatures, dark matter galaxy rotation, radial velocity validation | -| `bci_bridge` | 8 | KHAOS EEG → physics parameter bridge: Kerr-Newman spin ← calm_index, Orch-OR coherence ← alpha power, Schumann nexus, Muse 2 UDP/OSC stack | -| `sonification` | 6 | Physics-to-audio mappings, spectral analysis, 48-file WAV corpus | - -**Data corpus** (in `simulations/data/`): 48 WAV · 64 PNG · 8 JSON telemetry files - -### `pattern_analysis` - -Rust FFT engine (Welch PSD, spectral slope β, centroid, entropy, flatness) + Python cross-modal analyzer (audio × visual × telemetry correlation matrices, PCA, Ward dendrogram). 47 WAVs analyzed, `unified_analysis.json` produced. - -**Connects to:** KHAOS (BCI telemetry), HELIOS (time-series correlation) - -## Quick Start (Arch Linux / iNFAMØUS) - -```bash -git clone && cd SUBSTRATE -bash install.sh # install deps, build Rust + CUDA kernels, run smoke tests -bash validate.sh # pre-switch checklist — must be all-green before deploying -./target/release/substrate run -``` - -> **Windows note:** `install.sh` / `validate.sh` are Bash scripts for the Arch target. -> On the dev machine use Git Bash or WSL to inspect them; `chmod +x` is a no-op on NTFS -> but the execute bit is preserved when pushed and cloned on Linux. - -## Running the pipeline - -```bash -python pipeline.py --module cycle_project -python pipeline.py --module all -``` - -## Hardware - -- GPU: RTX 5060 Ti (sm_120, Blackwell, 16GB) -- PyTorch nightly cu128 -- Rust 2021 edition diff --git a/benchmarks/plot_benchmarks.py b/benchmarks/plot_benchmarks.py deleted file mode 100644 index 17dc6ce..0000000 --- a/benchmarks/plot_benchmarks.py +++ /dev/null @@ -1,275 +0,0 @@ -""" -SUBSTRATE — Benchmark Visualization -Generates publication-quality benchmark plots from existing data. -Output: benchmarks/plots/ -""" -import json -import os -import numpy as np -import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt -import matplotlib.patches as mpatches -from pathlib import Path - -RESULTS = Path(__file__).parent / "results" -PLOTS = Path(__file__).parent / "plots" -PLOTS.mkdir(exist_ok=True) - -DARK = "#0d1117" -ACCENT = "#00d4ff" -GREEN = "#39ff14" -ORANGE = "#ff6b35" -PURPLE = "#b388ff" -TEXT = "#e6edf3" -GRID = "#21262d" - -plt.rcParams.update({ - "figure.facecolor": DARK, "axes.facecolor": DARK, - "axes.edgecolor": GRID, "axes.labelcolor": TEXT, - "xtick.color": TEXT, "ytick.color": TEXT, - "text.color": TEXT, "grid.color": GRID, - "grid.linewidth": 0.5, "font.family": "monospace", - "font.size": 11, -}) - -# ── 1. U(1) plaquette CUDA kernel: roofline crossover vs JAX-CPU ─── -def plot_plaquette_roofline(): - """ - Honest CUDA story: hand-written U(1) plaquette kernel (CuPy RawKernel, - sm_120) vs a warmed @jax.jit reference running on CPU (no GPU JAX backend - in this build). Source data: benchmarks/results/plaquette_cuda.json, - produced by quantum_lab/P3_G2/benchmark_plaquette.py (warmed both sides, - CUDA-event kernel timing, 1000-iter medians). The point is the crossover, - not one hero number: kernel-only scales with problem size; end-to-end - (with H2D/D2H transfers) only wins past L≈128. - """ - src = RESULTS / "plaquette_cuda.json" - if not src.exists(): - print(" skip: no plaquette_cuda.json") - return - rows = json.loads(src.read_text()) - L = [r["L"] for r in rows] - k_spd = [r["kernel_speedup"] for r in rows] - e_spd = [r["e2e_speedup"] for r in rows] - jax_ms = [r["jax_ms"] for r in rows] - k_ms = [r["cuda_kernel_ms"] for r in rows] - e_ms = [r["cuda_e2e_ms"] for r in rows] - - fig, axes = plt.subplots(1, 2, figsize=(12, 5.0)) - fig.suptitle("SUBSTRATE — U(1) Plaquette CUDA Kernel vs JAX-CPU\n" - "RawKernel · sm_120 Blackwell (RTX 5060 Ti) · CUDA-event timing, 1000-iter medians", - fontsize=13, fontweight="bold", color=TEXT, y=1.03) - - # Left: speedup vs lattice size (the crossover) - ax = axes[0] - ax.plot(L, k_spd, color=ACCENT, marker="o", markersize=4, linewidth=1.6, - zorder=4, label="kernel-only (compute)") - ax.plot(L, e_spd, color=ORANGE, marker="s", markersize=4, linewidth=1.6, - zorder=4, label="end-to-end (+H2D/D2H)") - ax.axhline(1.0, color=GREEN, linestyle="--", linewidth=1.0, zorder=2, - label="break-even (1×)") - ax.set_xscale("log", base=2); ax.set_yscale("log") - ax.set_xticks(L); ax.set_xticklabels([str(x) for x in L]) - ax.set_xlabel("Lattice size L (cells = 2·L²)") - ax.set_ylabel("Speedup over JAX-CPU (×, log)") - ax.set_title("Roofline crossover: overhead-bound → compute-bound", color=TEXT) - ax.legend(fontsize=8, loc="upper left"); ax.grid(zorder=0, which="both", alpha=0.4) - ax.annotate(f"{k_spd[-1]:.0f}×", xy=(L[-1], k_spd[-1]), xytext=(-4, 6), - textcoords="offset points", ha="right", fontsize=11, - fontweight="bold", color=ACCENT) - ax.annotate("e2e break-even\n≈ L 128", xy=(128, 1.0), xytext=(0, 18), - textcoords="offset points", ha="center", fontsize=7, - color=GREEN, alpha=0.9) - - # Right: absolute time vs size (log-log) - ax = axes[1] - ax.plot(L, jax_ms, color=PURPLE, marker="o", markersize=4, linewidth=1.6, - zorder=4, label="JAX-CPU") - ax.plot(L, k_ms, color=ACCENT, marker="o", markersize=4, linewidth=1.6, - zorder=4, label="CUDA kernel-only") - ax.plot(L, e_ms, color=ORANGE, marker="s", markersize=4, linewidth=1.6, - zorder=4, label="CUDA end-to-end") - ax.set_xscale("log", base=2); ax.set_yscale("log") - ax.set_xticks(L); ax.set_xticklabels([str(x) for x in L]) - ax.set_xlabel("Lattice size L") - ax.set_ylabel("Wall time per call (ms, log)") - ax.set_title("Absolute latency — transfers dominate at small L", color=TEXT) - ax.legend(fontsize=8, loc="upper left"); ax.grid(zorder=0, which="both", alpha=0.4) - fig.tight_layout(rect=(0, 0.05, 1, 1)) - fig.text(0.5, 0.01, - "JAX runs on CPU backend (no GPU JAX in this build) - speedup is GPU-vs-CPU, stated honestly", - ha="center", va="bottom", fontsize=7, color=TEXT, alpha=0.5) - out = PLOTS / "plaquette_roofline.png" - fig.savefig(out, dpi=150, bbox_inches="tight", facecolor=DARK) - plt.close(fig) - print(f" [ok] {out}") - - -# ── 2. TDVP integration convergence ─────────────────────────────── -def plot_tdvp_convergence(): - data = [] - with open(RESULTS / "tdvp_timing.jsonl") as f: - for line in f: - d = json.loads(line) - if d.get("event") == "integrate": - data.append(d) - if not data: - print(" skip: no tdvp integrate events") - return - - runs = list(range(len(data))) - times = [d["integrate_time_s"] for d in data] - errors = [d["final_trace"] for d in data] # |1 - tr(rho)|, ideal = 0 - systems = [d["system"] for d in data] - dims = {d["system"]: d.get("dim", 0) for d in data} - - # Honest grouping: these are DIFFERENT physical systems (different Hilbert - # dimension), NOT the same system at two integrator tolerances. - uniq = list(dict.fromkeys(systems)) - palette = {s: col for s, col in zip(uniq, [ACCENT, PURPLE, ORANGE, GREEN])} - - fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) - fig.suptitle("SUBSTRATE — TDVP Integration (radical-pair systems, CPU RK45)", - fontsize=13, fontweight="bold", color=TEXT) - - # Left: wall time per run, grouped by system - ax = axes[0] - for s in uniq: - idx = [i for i in runs if systems[i] == s] - ts = [times[i] for i in idx] - ax.scatter(idx, ts, color=palette[s], s=20, zorder=4, - label=f"{s} (dim {dims[s]}) - mean {np.mean(ts):.0f}s") - ax.set_xlabel("Run #"); ax.set_ylabel("Wall time (s)") - ax.set_title("Integration time per run (100 steps, 5 us)", color=TEXT) - ax.legend(fontsize=8); ax.grid(zorder=0) - - # Right: trace error per run, grouped by system (lower = better) - ax = axes[1] - for s in uniq: - idx = [i for i in runs if systems[i] == s] - es = [errors[i] for i in idx] - ax.scatter(idx, es, color=palette[s], s=20, zorder=4, - label=f"{s}: |1-tr(rho)| ~ {np.mean(es):.2e}") - ax.axhline(0.0, color=GREEN, linestyle="--", linewidth=1, label="ideal (err = 0)") - ax.set_xlabel("Run #"); ax.set_ylabel("|1 - tr(rho)| (trace error, lower = better)") - ax.set_title("Trace conservation per system", color=TEXT) - ax.set_ylim(-0.0005, max(errors) * 1.3) - ax.legend(fontsize=8); ax.grid(zorder=0) - - fig.tight_layout() - out = PLOTS / "tdvp_convergence.png" - fig.savefig(out, dpi=150, bbox_inches="tight", facecolor=DARK) - plt.close(fig) - print(f" [ok] {out}") - - -# ── 3. ErCry4a radical pair sensitivity ─────────────────────────── -def plot_ercry4a(): - data = [] - with open(RESULTS / "ercry4a_benchmark.jsonl") as f: - for line in f: - d = json.loads(line) - if "phi_s" in d: - data.append(d) - if not data: - print(" skip: no ercry4a phi_s data") - return - - n = len(data) - phi_0 = [d["phi_s"][0] for d in data] - phi_50 = [d["phi_s"][1] for d in data] - delta = [abs(d["delta_phi_s_earth"]) * 1000 for d in data] - times = [d["wall_time_s"] for d in data] - runs = np.arange(n) - w = 0.30 - - fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) - fig.suptitle("SUBSTRATE — ErCry4a Magnetoreception (6 nuclear spins, 8-site, exact χ)", - fontsize=13, fontweight="bold", color=TEXT) - - # Left: singlet yield — grouped bars per run - ax = axes[0] - b1 = ax.bar(runs - w/2, phi_0, w, color=ACCENT, alpha=0.9, label="φₛ 0 mT (dark)", zorder=3) - b2 = ax.bar(runs + w/2, phi_50, w, color=PURPLE, alpha=0.9, label="φₛ 0.05 mT (Earth)", zorder=3) - ax.set_xticks(runs); ax.set_xticklabels([f"Run {i}" for i in runs]) - ax.set_ylabel("Singlet yield φₛ") - ax.set_ylim(0.050, 0.062) - ax.set_title("Singlet yield vs applied field", color=TEXT) - ax.legend(fontsize=9); ax.grid(axis="y", zorder=0) - # value labels - for bar, val in zip(b1, phi_0): - ax.text(bar.get_x() + bar.get_width()/2, val + 0.0001, - f"{val:.5f}", ha="center", va="bottom", fontsize=8, color=ACCENT) - for bar, val in zip(b2, phi_50): - ax.text(bar.get_x() + bar.get_width()/2, val + 0.0001, - f"{val:.5f}", ha="center", va="bottom", fontsize=8, color=PURPLE) - # wall time annotation - for i, t in enumerate(times): - ax.annotate(f"{t:.0f}s", xy=(i, 0.050), xycoords=("data","axes fraction"), - xytext=(0, 4), textcoords="offset points", - ha="center", fontsize=7, color=TEXT, alpha=0.5) - - # Right: Earth-field sensitivity ΔφS × 10³ - ax = axes[1] - bars = ax.bar(runs, delta, width=0.4, color=GREEN, alpha=0.9, zorder=3) - ax.set_xticks(runs); ax.set_xticklabels([f"Run {i}" for i in runs]) - ax.set_ylabel("|ΔφS| × 10³ (Earth-field sensitivity)") - ax.set_title("Radical pair Earth-field sensitivity", color=TEXT) - ax.set_ylim(0, max(delta) * 1.5) - ax.axhline(np.mean(delta), color=ORANGE, linestyle="--", linewidth=1, - label=f"mean = {np.mean(delta):.3f}×10⁻³") - ax.legend(fontsize=9); ax.grid(axis="y", zorder=0) - for bar, val in zip(bars, delta): - ax.text(bar.get_x() + bar.get_width()/2, val + 0.05, - f"{val:.3f}", ha="center", va="bottom", - fontsize=12, fontweight="bold", color=GREEN) - ax.text(0.98, 0.04, "deterministic solver — consistent across runs", - transform=ax.transAxes, ha="right", va="bottom", - fontsize=7, color=TEXT, alpha=0.45) - - fig.tight_layout() - out = PLOTS / "ercry4a_sensitivity.png" - fig.savefig(out, dpi=150, bbox_inches="tight", facecolor=DARK) - plt.close(fig) - print(f" [ok] {out}") - - -# ── 4. GPU Krylov build timing ───────────────────────────────────── -def plot_gpu_krylov(): - times = [] - with open(RESULTS / "gpu_timing.jsonl") as f: - for line in f: - d = json.loads(line) - if d.get("event") == "gpu_krylov_build": - times.append(d["build_time_ms"]) - if not times: - print(" skip: no gpu_krylov data") - return - - fig, ax = plt.subplots(figsize=(9, 4)) - fig.suptitle("SUBSTRATE — GPU Krylov Build Time (FAD-W, dim=64, sm_120)", - fontsize=13, fontweight="bold", color=TEXT) - ax.plot(times, color=ACCENT, linewidth=1.5, marker="o", markersize=3, zorder=3) - ax.axhline(np.mean(times), color=ORANGE, linestyle="--", linewidth=1.2, - label=f"mean = {np.mean(times):.1f} ms") - ax.fill_between(range(len(times)), np.mean(times)-np.std(times), - np.mean(times)+np.std(times), color=ACCENT, alpha=0.1) - ax.set_xlabel("Build #"); ax.set_ylabel("Wall time (ms)") - ax.set_title("RTX 5060 Ti — GPU Krylov matrix build latency", color=TEXT) - ax.legend(fontsize=10); ax.grid(zorder=0) - fig.tight_layout() - out = PLOTS / "gpu_krylov_timing.png" - fig.savefig(out, dpi=150, bbox_inches="tight", facecolor=DARK) - plt.close(fig) - print(f" [ok] {out}") - - -if __name__ == "__main__": - print("SUBSTRATE benchmark plots") - print("=" * 40) - plot_plaquette_roofline() - plot_tdvp_convergence() - plot_ercry4a() - plot_gpu_krylov() - print(f"\nAll plots -> {PLOTS}/") diff --git a/benchmarks/plots/ercry4a_sensitivity.png b/benchmarks/plots/ercry4a_sensitivity.png index 24c9daa..db10291 100644 Binary files a/benchmarks/plots/ercry4a_sensitivity.png and b/benchmarks/plots/ercry4a_sensitivity.png differ diff --git a/benchmarks/plots/gpu_krylov_timing.png b/benchmarks/plots/gpu_krylov_timing.png index be45884..37d1aba 100644 Binary files a/benchmarks/plots/gpu_krylov_timing.png and b/benchmarks/plots/gpu_krylov_timing.png differ diff --git a/benchmarks/plots/plaquette_roofline.png b/benchmarks/plots/plaquette_roofline.png index e2f5cd7..57b0676 100644 Binary files a/benchmarks/plots/plaquette_roofline.png and b/benchmarks/plots/plaquette_roofline.png differ diff --git a/benchmarks/plots/tdvp_convergence.png b/benchmarks/plots/tdvp_convergence.png index 3b55396..599004e 100644 Binary files a/benchmarks/plots/tdvp_convergence.png and b/benchmarks/plots/tdvp_convergence.png differ diff --git a/benchmarks/results/tdvp_timing.jsonl b/benchmarks/results/tdvp_timing.jsonl deleted file mode 100644 index 706c670..0000000 --- a/benchmarks/results/tdvp_timing.jsonl +++ /dev/null @@ -1,47 +0,0 @@ -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 135.39} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 31.895, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 118.84} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 31.766, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 137.01} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 29.59, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 123.28} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 29.965, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 115.26} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 31.929, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 122.24} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 31.708, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 120.12} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 31.551, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 150.4} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 32.034, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 124.47} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 31.476, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 148.08} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 30.901, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 122.71} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 128.51} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 126.93} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 33.478, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 127.44} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 30.532, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 120.57} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 26.533, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 146.24} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 29.242, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 158.81} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 36.565, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 192.37} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 36.495, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 167.81} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 29.409, "final_trace": 0.007334} -{"event": "build", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "build_time_ms": 125.4} -{"event": "integrate", "system": "FAD-W_model", "n_sites": 6, "dim": 64, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 31.495, "final_trace": 0.007334} -{"event": "build", "system": "ErCry4a_10nuc", "n_sites": 12, "dim": 4096, "build_time_ms": 28578.23} -{"event": "build", "system": "ErCry4a_6nuc", "n_sites": 8, "dim": 256, "build_time_ms": 269.58} -{"event": "integrate", "system": "ErCry4a_6nuc", "n_sites": 8, "dim": 256, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 89.514, "final_trace": 0.000258} -{"event": "build", "system": "ErCry4a_6nuc", "n_sites": 8, "dim": 256, "build_time_ms": 263.98} -{"event": "integrate", "system": "ErCry4a_6nuc", "n_sites": 8, "dim": 256, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 88.965, "final_trace": 0.000258} -{"event": "build", "system": "ErCry4a_6nuc", "n_sites": 8, "dim": 256, "build_time_ms": 265.78} -{"event": "integrate", "system": "ErCry4a_6nuc", "n_sites": 8, "dim": 256, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 89.805, "final_trace": 0.000258} -{"event": "build", "system": "ErCry4a_6nuc", "n_sites": 8, "dim": 256, "build_time_ms": 261.2} -{"event": "integrate", "system": "ErCry4a_6nuc", "n_sites": 8, "dim": 256, "t_max_us": 5.0, "n_steps": 100, "method": "RK45", "integrate_time_s": 89.505, "final_trace": 0.000258} diff --git a/engine/dll_healing.py b/engine/dll_healing.py deleted file mode 100644 index f342e36..0000000 --- a/engine/dll_healing.py +++ /dev/null @@ -1,72 +0,0 @@ -import os -import logging -from pathlib import Path - -logger = logging.getLogger("substrate.dll_healing") - -_healed = False # guard: os.add_dll_directory() accumulates on every call - -def heal(): - """ - Ensures that CUDA and Nvidia DLLs are found by the Python interpreter on Windows. - This is critical for libraries like CuPy and cuQuantum when running inside an - embedded environment or specialized conda/pip setups. - """ - global _healed - if _healed: - return - if os.name != 'nt': - _healed = True - return - - logger.debug("Initializing Windows DLL healing...") - - # 1. Check CUDA_PATH environment variables - cuda_keys = ['CUDA_PATH', 'CUDA_PATH_V13_0', 'CUDA_PATH_V12_5', 'CUDA_PATH_V12_4'] - for key in cuda_keys: - path = os.environ.get(key) - if path: - bin_path = Path(path) / "bin" - if bin_path.exists(): - logger.debug(f"[DLL] Adding {key}: {bin_path}") - os.add_dll_directory(str(bin_path)) - - # 2. Check common default locations - defaults = [ - r"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v13.0\bin", - r"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.5\bin", - r"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.4\bin", - r"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\bin", - ] - for d in defaults: - if Path(d).exists(): - logger.debug(f"[DLL] Adding default: {d}") - os.add_dll_directory(d) - - # 3. Check site-packages (nvidia-* wheels) - try: - import site - search_paths = site.getsitepackages() - if site.getusersitepackages(): - search_paths.append(site.getusersitepackages()) - - for sp in search_paths: - nvidia_path = Path(sp) / "nvidia" - if nvidia_path.exists(): - subpackages = [ - "cublas", "cusolver", "cufft", "curand", - "cusparse", "cuda_runtime", "cuda_nvrtc", "cudnn", - "nvjitlink", "nvfatbin" - ] - for sub in subpackages: - bin_path = nvidia_path / sub / "bin" - if bin_path.exists(): - logger.debug(f"[DLL] Adding site-package: {bin_path}") - os.add_dll_directory(str(bin_path)) - except Exception as e: - logger.warning(f"DLL healing site-package scan failed: {e}") - - _healed = True - -# Run healing on import -heal() diff --git a/modules/cryptotn_gpu/cryptotn/cuda/engine.py b/modules/cryptotn_gpu/cryptotn/cuda/engine.py deleted file mode 100644 index 575ab8c..0000000 --- a/modules/cryptotn_gpu/cryptotn/cuda/engine.py +++ /dev/null @@ -1,1440 +0,0 @@ -""" -cryptotn/cuda/engine.py -━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ -Phase B: GPU-accelerated solvers for radical pair spin dynamics. - -Two backends -──────────── -CupyKrylovSolver — GPU sparse L + Arnoldi Krylov expm - • drop-in for ExactSolver / MpsSolver (N ≤ 14 sites) - • validated: RMSE vs ExactSolver < 1e-4 (max|ΔP_S| ~ 1e-8 head-to-head) - • NOTE: at these dimensions (4^n_sites ≤ ~1024) this is launch/sync-bound and - runs SLOWER than scipy on CPU — each step is dozens of tiny sparse matvecs - plus a host expm round-trip, so per-call overhead dominates the arithmetic. - It exists for correctness/parity and as the GPU matvec path that the dense - regime will outgrow; the dense GPU win is not here. See - benchmarks/bench_lindblad_headtohead.py for the honest measurement. - -CuTDVPSolver — MPO-MPS TDVP on GPU (the real contribution) - • Liouvillian built as MPO from local operators (never densifies L_super) - • makes large N tractable by χ truncation: N=39 spins → Liouville dim - 4^39 ≈ 3e23, impossible to store densely, handled on 16 GB via MPS - • MPS bond dimension χ up to ~1800 on RTX 5060 Ti 16 GB (N=62) - • χ=2500 feasible for N ≤ 40 - • cupy einsum for environment contractions + Arnoldi Krylov for site update - -Units: μs / MHz / mT (same as Phase A) - -VRAM budget at χ=2500, N=62: - MPS tensors : 62 × 2500² × 4 × 16 B ≈ 24.8 GB → need χ≤1800 for full N=62 - Environments: 2 × 62 × 2500² × 14 × 16 B ≈ 174 GB → dominant, use selective recompute - Practical: χ=1024 for N=62, χ=2500 for N≤20 -""" -from __future__ import annotations - -import json -import logging -import time -from pathlib import Path -from typing import List, Optional, Tuple - -import numpy as np -from scipy import sparse - -from ..hamiltonian import build_recombination - -logger = logging.getLogger(__name__) -_BENCH_LOG = Path("benchmarks/results/gpu_timing.jsonl") - -# ── optional GPU imports ────────────────────────────────────────────────────── -try: - import cupy as cp - import cupyx.scipy.sparse as cpsp - HAS_CUPY = True -except ImportError: - cp = None - cpsp = None - HAS_CUPY = False - logger.warning("cupy not found — GPU backends disabled") - -try: - from cuquantum import Network # noqa: F401 - HAS_CUTN = True -except ImportError: - HAS_CUTN = False - - -def _log(rec: dict) -> None: - _BENCH_LOG.parent.mkdir(parents=True, exist_ok=True) - with open(_BENCH_LOG, "a") as f: - f.write(json.dumps(rec) + "\n") - - -def _require_cupy() -> None: - if not HAS_CUPY: - raise ImportError( - "cupy is required for GPU backends.\n" - " pip install cupy-cuda12x" - ) - - -# ═══════════════════════════════════════════════════════════════════════════════ -# §1 Liouville-space local operators (4×4 in |ket⟩⊗|bra⟩ basis) -# -# Basis: |↑↑⟩ |↑↓⟩ |↓↑⟩ |↓↓⟩ -# Left action L(O): O ρ → (O ⊗ I₂) vec(ρ) = kron(O, I₂) -# Right action R(O): ρ O → (I₂ ⊗ Oᵀ) vec(ρ) = kron(I₂, O.T) -# ═══════════════════════════════════════════════════════════════════════════════ - -_I2 = np.eye(2, dtype=complex) -_Sx = np.array([[0, 0.5], [0.5, 0]], dtype=complex) -_Sy = np.array([[0, -0.5j], [0.5j, 0]], dtype=complex) -_Sz = np.array([[0.5, 0], [0, -0.5]], dtype=complex) -_I4 = np.eye(4, dtype=complex) - - -def _L(O: np.ndarray) -> np.ndarray: - """4×4 left-action of 2×2 operator O in Liouville space.""" - return np.kron(O, _I2) - - -def _R(O: np.ndarray) -> np.ndarray: - """4×4 right-action of 2×2 operator O in Liouville space.""" - return np.kron(_I2, O.T) - - -# Precomputed 4×4 Liouville-space spin matrices -_Lx, _Rx = _L(_Sx), _R(_Sx) -_Ly, _Ry = _L(_Sy), _R(_Sy) -_Lz, _Rz = _L(_Sz), _R(_Sz) - - -def _comm_op(O: np.ndarray, omega: float, TWO_PI: float) -> np.ndarray: - """−i 2π ω [O, ·] as 4×4 Liouville superoperator.""" - return -1j * TWO_PI * omega * (_L(O) - _R(O)) - - -def _anticomm_op(O: np.ndarray, rate: float, TWO_PI: float) -> np.ndarray: - """−½ 2π rate {O, ·} as 4×4 Liouville superoperator.""" - return -0.5 * TWO_PI * rate * (_L(O) + _R(O)) - - -# ═══════════════════════════════════════════════════════════════════════════════ -# §2 MPO builder for the radical pair Liouvillian -# -# L = −i 2π [H, ·] − ½ 2π {K, ·} -# -# FSM bond states (W = 14): -# 0 PASS — identity pass-through -# 1-3 L1x/y/z — left-action strings from e1 (hyperfine + exchange) -# 4-6 R1x/y/z — right-action strings from e1 -# 7-9 L2x/y/z — left-action strings from e2 (hyperfine) -# 10-12 R2x/y/z — right-action strings from e2 -# 13 ACCUM — accumulator (completed terms land here) -# -# Bond dim W=14 is independent of chain length N. -# MPO never requires building the dense L_super. -# ═══════════════════════════════════════════════════════════════════════════════ - -_W = 14 -_PASS = 0 -_L1x, _L1y, _L1z = 1, 2, 3 -_R1x, _R1y, _R1z = 4, 5, 6 -_L2x, _L2y, _L2z = 7, 8, 9 -_R2x, _R2y, _R2z = 10, 11, 12 -_ACCUM = 13 - -_GAMMA_E = 27.994 # MHz / mT - - -def _zero_W() -> np.ndarray: - return np.zeros((_W, _W, 4, 4), dtype=complex) - - -def _mpo_electron1(omega_e1: float, J_MHz: float, - k_S: float, k_T: float, TWO_PI: float) -> np.ndarray: - """ - MPO W-matrix at electron-1 site. Shape (W, W, 4, 4). - - - Starts L1/R1 strings (for HF to e1-nuclei and exchange to e2). - - Applies Zeeman(e1) and half the scalar-K one-body contribution. - """ - M = _zero_W() - P = TWO_PI - - M[_PASS, _PASS] = _I4 - M[_ACCUM, _ACCUM] = _I4 - - # Start left-action strings (L1x, L1y, L1z) and right-action (R1x, R1y, R1z) - for li, ri, Lop, Rop in [ - (_L1x, _R1x, _Lx, _Rx), - (_L1y, _R1y, _Ly, _Ry), - (_L1z, _R1z, _Lz, _Rz), - ]: - M[li, _PASS] = Lop - M[ri, _PASS] = Rop - - # Zeeman: −i 2π Ω_e1 [Sz, ·] - M[_ACCUM, _PASS] += _comm_op(_Sz, omega_e1, P) - - # Scalar K (half, split symmetrically with e2): - # K_scalar = (k_S + 3 k_T) / 4 from Q_S + Q_T decomposition - # −½ {k_scalar × I, ·} = −k_scalar × I₄ in Liouville space - k_sc = (k_S + 3.0 * k_T) / 4.0 / 2.0 - M[_ACCUM, _PASS] += -P * k_sc * _I4 - - return M - - -def _mpo_nucleus_e1(A_tensor_MHz: np.ndarray, TWO_PI: float) -> np.ndarray: - """ - MPO W-matrix at a nucleus of electron-1. Shape (W, W, 4, 4). - - Carries all active strings via identity and completes HF(e1) terms. - """ - M = _zero_W() - P = TWO_PI - - M[_PASS, _PASS] = _I4 - M[_ACCUM, _ACCUM] = _I4 - - # Carry all active strings - for idx in [_L1x, _L1y, _L1z, _R1x, _R1y, _R1z, - _L2x, _L2y, _L2z, _R2x, _R2y, _R2z]: - M[idx, idx] = _I4 - - # Complete HF(e1): −i 2π A_αα [Se1_α In_α, ·] - # L1α string carries L(Se1_α). At nucleus: complete with −i 2π A × L(In_α) - # R1α string carries R(Se1_α). At nucleus: complete with +i 2π A × R(In_α) - A = A_tensor_MHz - if np.ndim(A) == 2: - Axx, Ayy, Azz = float(A[0, 0]), float(A[1, 1]), float(A[2, 2]) - else: - Axx = Ayy = Azz = float(A) - - M[_ACCUM, _L1x] = -1j * P * Axx * _Lx - M[_ACCUM, _R1x] = +1j * P * Axx * _Rx - M[_ACCUM, _L1y] = -1j * P * Ayy * _Ly - M[_ACCUM, _R1y] = +1j * P * Ayy * _Ry - M[_ACCUM, _L1z] = -1j * P * Azz * _Lz - M[_ACCUM, _R1z] = +1j * P * Azz * _Rz - - return M - - -def _mpo_electron2(omega_e2: float, J_MHz: float, - k_S: float, k_T: float, TWO_PI: float) -> np.ndarray: - """ - MPO W-matrix at electron-2 site. Shape (W, W, 4, 4). - - - Completes exchange strings from e1 (L1/R1 → ACCUM with J × Se2). - - Completes S·S part of recombination K from e1. - - Starts L2/R2 strings for e2's own nuclei. - - Applies Zeeman(e2) and second half of scalar K. - """ - M = _zero_W() - P = TWO_PI - - M[_PASS, _PASS] = _I4 - M[_ACCUM, _ACCUM] = _I4 - - # Complete exchange from e1: −i 2π J [S1_α S2_α, ·] - # L1α carries L(S1_α); complete with −i 2π J L(S2_α) → ACCUM - # R1α carries R(S1_α); complete with +i 2π J R(S2_α) → ACCUM - for l1, r1, Lop, Rop in [ - (_L1x, _R1x, _Lx, _Rx), - (_L1y, _R1y, _Ly, _Ry), - (_L1z, _R1z, _Lz, _Rz), - ]: - M[_ACCUM, l1] += -1j * P * J_MHz * Lop - M[_ACCUM, r1] += +1j * P * J_MHz * Rop - - # Complete S·S part of K: −½ 2π (k_T − k_S) {S1·S2, ·} - # = −½ 2π (k_T−k_S) Σ_α ( L(S1_α) L(S2_α) + R(S1_α) R(S2_α) ) - k_SS = -0.5 * P * (k_T - k_S) - for l1, r1, Lop, Rop in [ - (_L1x, _R1x, _Lx, _Rx), - (_L1y, _R1y, _Ly, _Ry), - (_L1z, _R1z, _Lz, _Rz), - ]: - M[_ACCUM, l1] += k_SS * Lop - M[_ACCUM, r1] += k_SS * Rop - - # Zeeman at e2 - M[_ACCUM, _PASS] += _comm_op(_Sz, omega_e2, P) - - # Scalar K (second half): −½ {k_scalar I, ·} = −k_scalar I₄ - k_sc = (k_S + 3.0 * k_T) / 4.0 / 2.0 - M[_ACCUM, _PASS] += -P * k_sc * _I4 - - # Start L2/R2 strings for e2's own nuclei - for li, ri, Lop, Rop in [ - (_L2x, _R2x, _Lx, _Rx), - (_L2y, _R2y, _Ly, _Ry), - (_L2z, _R2z, _Lz, _Rz), - ]: - M[li, _PASS] = Lop - M[ri, _PASS] = Rop - - return M - - -def _mpo_nucleus_e2(A_tensor_MHz: np.ndarray, TWO_PI: float) -> np.ndarray: - """MPO W-matrix at a nucleus of electron-2. Shape (W, W, 4, 4).""" - M = _zero_W() - P = TWO_PI - - M[_PASS, _PASS] = _I4 - M[_ACCUM, _ACCUM] = _I4 - - for idx in [_L1x, _L1y, _L1z, _R1x, _R1y, _R1z, - _L2x, _L2y, _L2z, _R2x, _R2y, _R2z]: - M[idx, idx] = _I4 - - A = A_tensor_MHz - if np.ndim(A) == 2: - Axx, Ayy, Azz = float(A[0, 0]), float(A[1, 1]), float(A[2, 2]) - else: - Axx = Ayy = Azz = float(A) - - M[_ACCUM, _L2x] = -1j * P * Axx * _Lx - M[_ACCUM, _R2x] = +1j * P * Axx * _Rx - M[_ACCUM, _L2y] = -1j * P * Ayy * _Ly - M[_ACCUM, _R2y] = +1j * P * Ayy * _Ry - M[_ACCUM, _L2z] = -1j * P * Azz * _Lz - M[_ACCUM, _R2z] = +1j * P * Azz * _Rz - - return M - - -def build_mpo(config, TWO_PI: float = 2.0 * np.pi) -> List[np.ndarray]: - """ - Build the Liouvillian MPO for a radical pair SystemConfig. - - Returns a list of n_sites numpy arrays, all shape (W, W, 4, 4). - Boundary conditions are enforced by the L[0] and R[n] environment - tensors (which select the PASS and ACCUM bond states respectively), - not by cropping the MPO tensor shapes. - """ - cfg = config - g1, g2 = cfg.radical_1.g_factor, cfg.radical_2.g_factor - omega_e1 = _GAMMA_E * g1 * cfg.B_mT - omega_e2 = _GAMMA_E * g2 * cfg.B_mT - J_MHz = float(getattr(cfg, "J_MHz", 0.0)) - k_S = float(cfg.k_S_us) - k_T = float(getattr(cfg, "k_T_us", 0.0)) - - sites: List[np.ndarray] = [] - sites.append(_mpo_electron1(omega_e1, J_MHz, k_S, k_T, TWO_PI)) - for nuc in cfg.radical_1.nuclei: - sites.append(_mpo_nucleus_e1(nuc.A_tensor_MHz, TWO_PI)) - sites.append(_mpo_electron2(omega_e2, J_MHz, k_S, k_T, TWO_PI)) - for nuc in cfg.radical_2.nuclei: - sites.append(_mpo_nucleus_e2(nuc.A_tensor_MHz, TWO_PI)) - - n = len(sites) - assert n == cfg.n_sites, f"MPO length {n} ≠ n_sites {cfg.n_sites}" - - # Convention fix: the FSM construction fills M[n_right, m_left, t, s] - # (rows = output bond, cols = input bond), but the environment contractions - # in _apply_heff / _update_left_env / _update_right_env expect - # W[m_left, n_right, t, s]. Transpose axes 0↔1 to reconcile. - mpo = [M.transpose(1, 0, 2, 3) for M in sites] - - logger.debug(f"MPO built: {n} sites, W={_W}") - return mpo - - -# ═══════════════════════════════════════════════════════════════════════════════ -# §3 MPS builder (vectorized density matrix → compressed MPS) -# ═══════════════════════════════════════════════════════════════════════════════ - -def build_mps(rho0: np.ndarray, n_sites: int, - chi_max: int = 64) -> List[np.ndarray]: - """ - Compress the initial density matrix into an MPS via left-to-right SVD. - - rho0 : (dim, dim) initial density matrix (dim = 2^n_sites) - n_sites : number of spin sites - chi_max : maximum bond dimension χ - - Returns: list of n_sites numpy arrays, each shape (χ_l, χ_r, 4). - Physical dim = 4 in Liouville space (|ket⟩⊗|bra⟩). - """ - d = 4 - dim = rho0.shape[0] - assert dim == 2 ** n_sites, \ - f"dim mismatch: rho0 is {rho0.shape}, expected ({2**n_sites},{2**n_sites})" - - # Vectorize into MPO interleaved Liouville ordering: - # each site k carries local index pk = sk_ket*2 + sk_bra (from L(O)=kron(O,I₂)) - # - # C-order flatten gives tensor with axes (s0k, s1k, ..., sNk, s0b, ..., sNb). - # We need (s0k, s0b, s1k, s1b, ..., sNk, sNb) so that each consecutive pair - # of axes folds into one size-4 local physical index pk = sk_ket*2 + sk_bra. - perm = [x for i in range(n_sites) for x in (i, n_sites + i)] - state = (rho0.reshape([2] * (2 * n_sites)) - .transpose(perm) - .reshape([d] * n_sites)) - - tensors: List[np.ndarray] = [] - current = state.reshape(1, -1) # (1, d^n) - - for _ in range(n_sites - 1): - chi_l = current.shape[0] - mat = current.reshape(chi_l * d, -1) - U, s, Vh = np.linalg.svd(mat, full_matrices=False) - chi_r = min(chi_max, len(s)) - U, s, Vh = U[:, :chi_r], s[:chi_r], Vh[:chi_r, :] - # U: (chi_l*d, chi_r) → reshape to (chi_l, d, chi_r) → transpose to (chi_l, chi_r, d) - T = U.reshape(chi_l, d, chi_r).transpose(0, 2, 1) - tensors.append(T) - current = np.diag(s) @ Vh - - chi_l = current.shape[0] - T_last = current.reshape(chi_l, 1, d) - tensors.append(T_last) - - return tensors - - -def mps_to_vector_np(tensors: List[np.ndarray]) -> np.ndarray: - """Contract MPS (numpy) to a full state vector. O(χ² d^n) — small N only.""" - # tensors[i] shape: (chi_l, chi_r, d) - result = tensors[0][0, :, :] # (chi_r, d) — left bond squeezed - for T in tensors[1:]: - chi_l, chi_r, d = T.shape - chi_prev = result.shape[0] - # Contract: result (chi_prev, d_all) with T (chi_l, chi_r, d) over chi_prev - r_mat = result.reshape(chi_prev, -1) # (chi_prev, d_all) - T_mat = T.reshape(chi_l, chi_r * d) # (chi_prev, chi_r*d) - result = r_mat.T @ T_mat # (d_all, chi_r*d) - result = result.reshape(-1, chi_r).T # (chi_r, d_all*d_new) — wrong - # Better sequential contraction: - result = np.einsum("...i,ijk->...jk", result.T, T).reshape(chi_r, -1) - # Hmm, getting complicated. Use a simple loop instead. - break - - # Simpler correct contraction: - result = tensors[0][0] # (chi_r, d) for site 0 - for T in tensors[1:]: - # result: (chi_prev, [physical dims so far]) - # T: (chi_l=chi_prev, chi_r, d) - # contract chi_prev: result@T[:,:,:] over first index - chi_prev = result.shape[0] - T_r = T.reshape(chi_prev, -1) # (chi_prev, chi_r*d) - result = result.reshape(chi_prev, -1) # (chi_prev, d_accum) - result = T_r.T @ result # (chi_r*d, d_accum) — wrong order - # Correct: sum over chi_prev - result = np.tensordot(result, T, axes=([0], [0])) - break # placeholder - - # Use a clean recursive approach: - vec = tensors[0].reshape(-1, tensors[0].shape[1]) # (d, chi_r) after squeezing left=1 - for T in tensors[1:]: - chi_l, chi_r, d = T.shape - # vec: (d_accum, chi_l) - # T: (chi_l, chi_r, d) - # contract over chi_l: (d_accum, chi_r, d) - vec = np.einsum("ai,ijk->ajk", vec, T).reshape(-1, chi_r) - - return vec.flatten() - - -# ═══════════════════════════════════════════════════════════════════════════════ -# §4 Krylov matrix-exponential (Arnoldi, GPU + CPU unified) -# ═══════════════════════════════════════════════════════════════════════════════ - -def _arnoldi_expm(matvec, v, dt: float, m: int = 40): - """ - Arnoldi approximation to exp(A dt) @ v. - - matvec : callable v → A @ v (works with cupy or numpy) - v : initial vector (cupy.ndarray or numpy.ndarray) - dt : time step (scalar, can be complex for backward integration) - m : Krylov subspace dimension - - Returns vector of same type/shape as v. - """ - xp = cp if (HAS_CUPY and isinstance(v, cp.ndarray)) else np - n = len(v) - m = min(m, n) - - V = xp.zeros((n, m + 1), dtype=complex) - H = xp.zeros((m + 1, m), dtype=complex) - - beta = xp.linalg.norm(v) - if float(xp.abs(beta)) < 1e-15: - return xp.zeros_like(v) - - V[:, 0] = v / beta - - j_stop = m - for j in range(m): - w = matvec(V[:, j]) - for i in range(j + 1): - H[i, j] = xp.dot(xp.conj(V[:, i]), w) - w = w - H[i, j] * V[:, i] - nrm = xp.linalg.norm(w) - H[j + 1, j] = nrm - if float(xp.abs(nrm)) < 1e-12: - j_stop = j + 1 - break - V[:, j + 1] = w / nrm - - m_eff = j_stop - from scipy.linalg import expm as _expm - H_np = H[:m_eff, :m_eff] - if HAS_CUPY and isinstance(H_np, cp.ndarray): - H_np = H_np.get() - eH = xp.array(_expm(H_np * dt)) - - e1 = xp.zeros(m_eff, dtype=complex) - e1[0] = 1.0 - return float(beta) * (V[:, :m_eff] @ (eH @ e1)) - - -# ═══════════════════════════════════════════════════════════════════════════════ -# §5 TDVP environment contractions -# -# Convention: -# MPS A[i] shape (χ_l, χ_r, 4) — (left_bond, right_bond, phys) -# MPO W[i] shape (W_l, W_r, 4, 4) — (mpo_l, mpo_r, phys_out, phys_in) -# L env shape (χ_bra_l, W_l, χ_ket_l) -# R env shape (χ_bra_r, W_r, χ_ket_r) -# -# Left env boundary L[0] = e_PASS (1, W, 1) with entry [0, PASS, 0]=1 -# Right env boundary R[n] = e_ACCUM (1, W, 1) with entry [0, ACCUM, 0]=1 -# ═══════════════════════════════════════════════════════════════════════════════ - -def _init_boundary_envs(n_sites: int, xp=np): - """Return (L[0], R[n]) boundary environment tensors.""" - L0 = xp.zeros((1, _W, 1), dtype=complex); L0[0, _PASS, 0] = 1.0 - Rn = xp.zeros((1, _W, 1), dtype=complex); Rn[0, _ACCUM, 0] = 1.0 - return L0, Rn - - -def _update_left_env(L_prev, A, W_mpo, xp=np): - """ - L[i+1][c, n, d] = Σ_{a,m,b,t,s} L[a,m,b] · A*[a,c,t] · W[m,n,t,s] · A[b,d,s] - - L_prev : (χ_bra_l, W_l, χ_ket_l) - A : (χ_l, χ_r, 4) - W_mpo : (W_l, W_r, 4, 4) W[m,n,t,s]: t=phys_out(bra), s=phys_in(ket) - Returns: (χ_bra_r, W_r, χ_ket_r) - """ - ein = xp.einsum - Ac = xp.conj(A) - # L[a,m,b], Ac[a,c,t], W[m,n,t,s], A[b,d,s] → [c,n,d] - tmp1 = ein("amb,act->mbct", L_prev, Ac) # (W_l, χ_ket_l, χ_bra_r, 4) - tmp2 = ein("mbct,mnts->bcns", tmp1, W_mpo) # (χ_ket_l, χ_bra_r, W_r, 4) - return ein("bcns,bds->cnd", tmp2, A) # (χ_bra_r, W_r, χ_ket_r) - - -def _update_right_env(R_next, A, W_mpo, xp=np): - """ - R[i][a, m, b] = Σ_{c,n,d,t,s} R[c,n,d] · A*[a,c,t] · W[m,n,t,s] · A[b,d,s] - - R_next : (χ_bra_r, W_r, χ_ket_r) - A : (χ_l, χ_r, 4) - W_mpo : (W_l, W_r, 4, 4) - Returns: (χ_bra_l, W_l, χ_ket_l) - """ - ein = xp.einsum - Ac = xp.conj(A) - # R[c,n,d], Ac[a,c,t], W[m,n,t,s], A[b,d,s] → [a,m,b] - tmp1 = ein("cnd,act->nadt", R_next, Ac) # (W_r, χ_ket_r, χ_bra_l, 4) - tmp2 = ein("mnts,nadt->mads", W_mpo, tmp1) # (W_l, χ_ket_r, χ_bra_l, 4) - res = ein("mads,bds->mab", tmp2, A) # (W_l, χ_bra_l, χ_ket_l) - return res.transpose(1, 0, 2) # (χ_bra_l, W_l, χ_ket_l) - - -def _apply_keff(L_env, R_env, C, xp=np): - """ - Zero-site effective Hamiltonian applied to bond matrix C. - - Used for the back-propagation step in 1-site TDVP: after site i is updated - and made left/right-canonical, the center matrix C is evolved BACKWARD by - exp(-K_eff × dt/2) to prevent double-counting of the Hamiltonian. - - K_eff result[a, c] = Σ_{b,m,d} L_env[a,m,b] · R_env[c,m,d] · C[b,d] - - L_env : (χ_l, W, χ_l) — left environment at the bond - R_env : (χ_r, W, χ_r) — right environment at the bond - C : (χ_l, χ_r) — bond center matrix - Returns (χ_l, χ_r) - """ - ein = xp.einsum - tmp = ein("amb,bd->amd", L_env, C) # (χ_l, W, χ_r) - return ein("amd,cmd->ac", tmp, R_env) # (χ_l, χ_r) - - -def _apply_heff(L, W_mpo, R, A, xp=np): - """ - (H_eff A)[i, k, t] = Σ_{j,m,n,l,s} L[i,m,j] · W[m,n,t,s] · R[k,n,l] · A[j,l,s] - - L : (χ_bra_l, W_l, χ_ket_l) - W_mpo : (W_l, W_r, 4, 4) [t=phys_out, s=phys_in] - R : (χ_bra_r, W_r, χ_ket_r) - A : (χ_l, χ_r, 4) - Returns: (χ_bra_l, χ_bra_r, 4) — same shape as A - """ - ein = xp.einsum - tmp1 = ein("imj,mnts->ijnts", L, W_mpo) # (χ_bra_l, χ_ket_l, W_r, 4, 4) - tmp2 = ein("ijnts,jls->intl", tmp1, A) # (χ_bra_l, W_r, 4, χ_ket_r) - res = ein("intl,knl->ikt", tmp2, R) # (χ_bra_l, χ_bra_r, 4) — no transpose needed - return res # (χ_bra_l, χ_bra_r, 4) = same shape as A - - -def _apply_heff_2site(L, W0, W1, R, Theta, xp=np): - """ - Two-site effective Hamiltonian applied to the super-tensor Θ. - - Contracts the 2-site tensor Θ with L, W0, W1, R to give H_eff_2site * Θ. - - L : (χ_bra_l, W_l, χ_ket_l) - W0 : (W_l, W_m, d, d) MPO tensor for site i - W1 : (W_m, W_r, d, d) MPO tensor for site i+1 - R : (χ_bra_r, W_r, χ_ket_r) - Theta : (χ_ket_l, d, d, χ_ket_r) — super-tensor A[i] ⊗ A[i+1] - Returns: (χ_bra_l, d, d, χ_bra_r) — same shape as Theta - """ - ein = xp.einsum - # Step 1: contract L with Theta over left bond - # L[i, m, j], Theta[j, s, t, k] → [i, m, s, t, k] - tmp = ein("imj,jstk->imstk", L, Theta) - # Step 2: contract with W0 over (m, s) → (i, n, t, s', k) - # W0[m, n, s', s]: m=W_l in, n=W_m out, s'=phys_out(bra), s=phys_in(ket) - tmp = ein("imstk,mnSs->inStk", tmp, W0) # (χ_bra_l, W_m, d_bra, d_ket, χ_ket_r) - # Step 3: contract with W1 over (n, t) → (i, p, S, T, k) - tmp = ein("inStk,npTt->ipSTk", tmp, W1) # (χ_bra_l, W_r, d_bra0, d_bra1, χ_ket_r) - # Step 4: contract with R over (p, k) - # R[l, p, k]: l=χ_bra_r, p=W_r, k=χ_ket_r - res = ein("ipSTk,lpk->ilST", tmp, R) # (χ_bra_l, χ_bra_r, d, d) - # Re-order to match Theta shape (χ_bra_l, d, d, χ_bra_r) - return res.transpose(0, 2, 3, 1) - - -# ═══════════════════════════════════════════════════════════════════════════════ -# §6 CupyKrylovSolver — GPU sparse L + Arnoldi expm (validated, N ≤ 14) -# ═══════════════════════════════════════════════════════════════════════════════ - -class CupyKrylovSolver: - """ - GPU-accelerated Krylov solver (drop-in for ExactSolver / MpsSolver). - - Transfers the sparse Liouvillian to GPU and integrates via time-stepping - with Arnoldi matrix exponential. - - Validated correct vs ExactSolver: max|ΔP_S| ~ 1e-8 across n_nuc=1..3 - (Liouville dim 64..1024), RMSE < 1e-5. - - Performance (honest): in this small-dimension regime the GPU path is - SLOWER than scipy on CPU — measured 5–87× slower in - benchmarks/bench_lindblad_headtohead.py. Reason: each step performs ~m tiny - sparse mat-vecs (m = krylov_dim) and a host-side scipy.expm on the m×m - Hessenberg, so kernel-launch + sync + H2D/D2H overhead dwarfs the actual - flops until the dense Liouvillian is far larger than fits at these N. Use it - for parity testing, not as a speed claim. The tractability win lives in - CuTDVPSolver (MPS), not here. - - Usage: - solver = CupyKrylovSolver(config, krylov_dim=50) - t, P_S, tr = solver.run(t_max_us=10.0, n_steps=500) - """ - - def __init__(self, config, krylov_dim: int = 50): - _require_cupy() - self.config = config - self.krylov_dim = krylov_dim - self._L_gpu = None - self._Q_S_gpu = None - self._rho0_gpu = None - self._dim = None - - def _build(self) -> None: - from ..hamiltonian import build_liouvillian - cfg = self.config - logger.info(f"[GPU-Krylov] building {cfg.name} ({cfg.n_sites} sites)…") - t0 = time.perf_counter() - - H = cfg.build_hamiltonian() - K = cfg.build_K() - L = build_liouvillian(H, K) - L_scaled = (2.0 * np.pi * L).tocsr() - - # Upload sparse L to GPU - self._L_gpu = cpsp.csr_matrix( - (cp.array(L_scaled.data), - cp.array(L_scaled.indices), - cp.array(L_scaled.indptr)), - shape=L_scaled.shape, - ) - - dim = H.shape[0] - e1, e2, _, _ = cfg.site_layout() - Q_S = build_recombination(e1, e2, cfg.n_sites, k_S=1.0, k_T=0.0).toarray() - self._Q_S_gpu = cp.array(Q_S) - self._rho0_gpu = cp.array(cfg.initial_rho().flatten(), dtype=complex) - self._dim = dim - - dt = time.perf_counter() - t0 - logger.info(f"[GPU-Krylov] build: {dt*1e3:.1f} ms, dim={dim}") - _log({"event": "gpu_krylov_build", "system": cfg.name, - "n_sites": cfg.n_sites, "dim": dim, - "build_time_ms": round(dt * 1e3, 2)}) - - def run( - self, - t_max_us: float = 10.0, - n_steps: int = 500, - ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: - """Integrate on GPU. Returns (t_us, P_S, trace).""" - if self._L_gpu is None: - self._build() - - cfg = self.config - dim = self._dim - t_eval = np.linspace(0.0, t_max_us, n_steps) - dt = t_eval[1] - t_eval[0] - - logger.info(f"[GPU-Krylov] {cfg.name}: {n_steps} steps, " - f"dt={dt:.4f} μs, m={self.krylov_dim}…") - t0 = time.perf_counter() - - rho = self._rho0_gpu.copy() - L_gpu = self._L_gpu - - def matvec(v: cp.ndarray) -> cp.ndarray: - return L_gpu @ v - - P_S = np.zeros(n_steps) - trace = np.zeros(n_steps) - - rho_mat = rho.reshape(dim, dim) - P_S[0] = float(cp.real(cp.trace(self._Q_S_gpu @ rho_mat))) - trace[0] = float(cp.real(cp.trace(rho_mat))) - - for i in range(1, n_steps): - rho = _arnoldi_expm(matvec, rho, dt, m=self.krylov_dim) - rho_mat = rho.reshape(dim, dim) - P_S[i] = float(cp.real(cp.trace(self._Q_S_gpu @ rho_mat))) - trace[i] = float(cp.real(cp.trace(rho_mat))) - - wall = time.perf_counter() - t0 - ms_step = wall / n_steps * 1e3 - logger.info(f"[GPU-Krylov] done: {wall:.2f}s ({ms_step:.2f} ms/step)") - _log({"event": "gpu_krylov_integrate", "system": cfg.name, - "n_sites": cfg.n_sites, "n_steps": n_steps, "t_max_us": t_max_us, - "krylov_dim": self.krylov_dim, "wall_time_s": round(wall, 3), - "ms_per_step": round(ms_step, 3)}) - - return t_eval, P_S, trace - - -# ═══════════════════════════════════════════════════════════════════════════════ -# §7 CuTDVPSolver — MPO-MPS 1-site TDVP on GPU (N up to ~62) -# ═══════════════════════════════════════════════════════════════════════════════ - -class CuTDVPSolver: - """ - 1-site TDVP in Liouville space on GPU via cupy einsum. - - Each time step: - 1. Build initial right environments R[0..n] (right-to-left sweep). - 2. Left-to-right sweep: for each site, apply exp(H_eff dt/2) via Krylov, - SVD-truncate to χ, shift center right, update left environment. - 3. Right-to-left sweep: symmetric. - - Memory: O(n χ² W d) for environments. W=14, d=4. - At χ=512, N=62: ≈ 62 × 512² × 14 × 4 × 16 B = 1.8 GB — fits comfortably. - At χ=1024, N=62: ≈ 7.2 GB. At χ=2500, N=40: ≈ 14 GB (tight). - - Usage: - solver = CuTDVPSolver(config, chi=256) - t, P_S, tr = solver.run(t_max_us=10.0, n_steps=200) - """ - - def __init__(self, config, chi: int = 64, krylov_dim: int = 30, - _numpy_mode: bool = False): - """ - _numpy_mode=True bypasses cupy requirement and uses numpy for CPU - validation testing. Not for production use. - """ - if not _numpy_mode: - _require_cupy() - self._xp = np if _numpy_mode else cp - self.config = config - self.chi = chi - self.krylov_dim = krylov_dim - self._mpo = None # list of numpy arrays - self._mps = None # list of xp arrays (χ_l, χ_r, 4) - self._Q_S_dense = None - - # ── initialization ──────────────────────────────────────────────────────── - - def _build(self) -> None: - cfg = self.config - logger.info(f"[TDVP-GPU] building {cfg.name} " - f"({cfg.n_sites} sites, chi={self.chi})...") - t0 = time.perf_counter() - - self._mpo = build_mpo(cfg) - - xp = self._xp - rho0 = cfg.initial_rho() - # Break singular-value degeneracy via tiny sparse Krylov kick. - # Pure-state (singlet) initial conditions have ALL equal MPS singular - # values, causing 1-site TDVP to be permanently frozen: the back- - # propagation exp(-K_eff dt/2) exactly cancels the forward Krylov - # exp(H_eff dt/2) for any degenerate eigenspace. A dt_kick ≈ 1e-6 μs - # evolution via the full sparse Liouvillian is negligible physically - # but makes the singular values non-degenerate so TDVP-1 can evolve. - rho0 = self._enrich_initial(rho0, cfg) - tensors_np = build_mps(rho0, cfg.n_sites, chi_max=self.chi) - self._mps = [xp.array(T, dtype=complex) for T in tensors_np] - - e1, e2, _, _ = cfg.site_layout() - self._Q_S_dense = build_recombination(e1, e2, cfg.n_sites, k_S=1.0, k_T=0.0).toarray() - - build_time = time.perf_counter() - t0 - logger.info(f"[TDVP-GPU] build: {build_time*1e3:.1f} ms") - _log({"event": "tdvp_build", "system": cfg.name, - "n_sites": cfg.n_sites, "chi": self.chi, - "build_time_ms": round(build_time * 1e3, 2)}) - - # ── observables ────────────────────────────────────────────────────────── - - def _observables(self) -> Tuple[float, float]: - """Extract P_S and trace from MPS via full vector contraction.""" - n = self.config.n_sites - xp = self._xp - - def _to_np(t): - """Move tensor to numpy (handles both cupy and numpy arrays).""" - return t.get() if (HAS_CUPY and xp is not np and - isinstance(t, cp.ndarray)) else np.asarray(t) - - # Sequential contraction: result stored as (d_accum, chi_right) - # result[phys_accumulated, bond] — bond is the dangling right index. - # - # Key: after einsum → (d_accum, chi_r_new, d_phys_new). - # Must transpose to (d_accum, d_phys_new, chi_r_new) BEFORE reshape - # so that physical indices accumulate cleanly in the first dimension. - # - # A[0] shape (1, chi_r, d): squeeze left bond → (chi_r, d) → transpose → (d, chi_r) - result = _to_np(self._mps[0]).squeeze(0).T # (d, chi_r) = (phys0, bond0) - - for T_gpu in self._mps[1:]: - T = _to_np(T_gpu) # (chi_l, chi_r, d) - chi_l, chi_r, d = T.shape - # einsum: result(d_accum, bond=chi_l), T(chi_l, chi_r, d_new) - # → (d_accum, chi_r_new, d_phys_new) - result = np.einsum("Ai,ijk->Ajk", result, T) - # transpose: (d_accum, chi_r_new, d_phys_new) → (d_accum, d_phys_new, chi_r_new) - # then reshape: (d_accum * d_phys_new, chi_r_new) - result = result.transpose(0, 2, 1).reshape(-1, chi_r) - - # Final: result shape (d_total, 1) → flatten to (4^n,) - # state_vec is in MPO interleaved ordering: - # idx = p0*4^(n-1)+p1*4^(n-2)+..., pk = sk_ket*2+sk_bra - # reshape to [2]*(2n) gives axes (s0k, s0b, s1k, s1b, ..., sNk, sNb) - # Apply inverse permutation → (s0k, s1k, ..., sNk, s0b, s1b, ..., sNb) = C-order - state_vec = result.flatten() - dim = 2 ** n - inv_perm = list(range(0, 2 * n, 2)) + list(range(1, 2 * n, 2)) - rho = (state_vec.reshape([2] * (2 * n)) - .transpose(inv_perm) - .reshape(dim, dim)) - P_S = float(np.real(np.trace(self._Q_S_dense @ rho))) - trace = float(np.real(np.trace(rho))) - return P_S, trace - - # ── warm-start enrichment ───────────────────────────────────────────────── - - @staticmethod - def _enrich_initial(rho0: np.ndarray, cfg) -> np.ndarray: - """ - Break degenerate MPS singular values so the first TDVP step can evolve. - - Two strategies depending on system size: - - Small systems (Liouville dim ≤ 2^20 ≈ 1M, i.e. n_sites ≤ 10): - Apply a tiny dt=1e-6 μs Krylov step via the full sparse Liouvillian. - Physically negligible but breaks all SV degeneracies exactly. - - Large systems (n_sites > 10): - The full Liouvillian build would require > 1 GB sparse matrix and - OOMs in WSL2. Instead, add Hermitian random noise at ε=1e-8, which - is sufficient for 2-site TDVP: the first SVD step already mixes - neighbouring sites and naturally resolves any remaining degeneracy. - 1-site TDVP needs the full Krylov kick; 2-site does not. - """ - dim = rho0.shape[0] - dim_L = dim * dim # Liouville space dimension - - _krylov_max_dim_L = 2 ** 20 # ≈ 1M → n_sites ≤ 10 - - if dim_L <= _krylov_max_dim_L: - # ── small system: full sparse Krylov kick ───────────────────────── - import scipy.sparse.linalg as spla - from ..hamiltonian import build_liouvillian - - H = cfg.build_hamiltonian() - K = cfg.build_K() - L = build_liouvillian(H, K) - L_sc = (2.0 * np.pi * L).tocsr() - - dt_kick = 1e-6 # 1 ns — negligible physically - rho_vec = spla.expm_multiply(L_sc * dt_kick, rho0.flatten()) - rho_enriched = rho_vec.reshape(dim, dim) - method = f"Krylov dt={dt_kick:.1e} μs" - else: - # ── large system: Hermitian random noise kick ───────────────────── - # Safe for 2-site TDVP: SVD step breaks residual degeneracy. - rng = np.random.default_rng(seed=42) - eps = 1e-8 - Z = rng.standard_normal(rho0.shape) + 1j * rng.standard_normal(rho0.shape) - Z = 0.5 * (Z + Z.conj().T) * eps - rho_enriched = rho0 + Z - method = f"random noise eps={eps:.1e}" - - tr = np.real(np.trace(rho_enriched)) - if tr > 0: - rho_enriched /= tr - - logger.info(f"[TDVP-GPU] warm-start kick ({method}, " - f"||δρ||={np.linalg.norm(rho_enriched - rho0):.3e})") - return rho_enriched - - # ── environment building ────────────────────────────────────────────────── - - def _build_right_envs(self) -> List: - """Build all right environments R[0..n] from right to left.""" - xp = self._xp - n = len(self._mps) - _, Rn = _init_boundary_envs(n, xp=xp) - right_envs = [None] * (n + 1) - right_envs[n] = Rn - - for i in range(n - 1, -1, -1): - A = self._mps[i] - W_mpo = xp.array(self._mpo[i], dtype=complex) - right_envs[i] = _update_right_env(right_envs[i + 1], A, W_mpo, xp=xp) - - return right_envs - - # ── TDVP sweeps ─────────────────────────────────────────────────────────── - - def _sweep_lr(self, left_envs: List, right_envs: List, dt: float) -> None: - """ - Left-to-right half-sweep with back-propagation (1-site TDVP). - - At each site i: - 1. Forward Krylov: A[i] ← exp(H_eff[i] × dt/2) A[i] - 2. SVD → left-canonical A[i] + center C - 3. Update left environment L[i+1] - 4. Back-propagation: C ← exp(−K_eff × dt/2) C (undo double-counting) - 5. Absorb C into A[i+1] - """ - xp = self._xp - n = len(self._mps) - mpo = self._mpo - - for i in range(n - 1): - A = self._mps[i] - W_mpo = xp.array(mpo[i], dtype=complex) - L = left_envs[i] - R = right_envs[i + 1] - chi_l, chi_r, d = A.shape - - # 1. Krylov: exp(H_eff × dt/2) @ A - def heff_mv(v, _L=L, _W=W_mpo, _R=R, - _cl=chi_l, _cr=chi_r, _d=d): - return _apply_heff(_L, _W, _R, - v.reshape(_cl, _cr, _d), xp=xp).flatten() - - A_new = _arnoldi_expm(heff_mv, A.flatten(), dt / 2.0, - m=self.krylov_dim).reshape(chi_l, chi_r, d) - - # 2. SVD → left-canonical + center C - # Reshape (chi_l*d, chi_r): merge left bond and physical dim - A_mat = A_new.transpose(0, 2, 1).reshape(chi_l * d, chi_r) - U, s, Vh = xp.linalg.svd(A_mat, full_matrices=False) - chi_new = min(self.chi, len(s)) - U, s, Vh = U[:, :chi_new], s[:chi_new], Vh[:chi_new, :] - - # Left-canonical tensor at site i: shape (chi_l, chi_new, d) - self._mps[i] = U.reshape(chi_l, d, chi_new).transpose(0, 2, 1) - - C = xp.diag(s) @ Vh # center: (chi_new, chi_r) - - # 3. Update left environment L[i+1] from new left-canonical A[i] - left_envs[i + 1] = _update_left_env(L, self._mps[i], W_mpo, xp=xp) - - # 4. Back-propagation: C ← exp(−K_eff × dt/2) C - # K_eff[a,c;b,d] = Σ_m L[i+1][a,m,b] · R[i+1][c,m,d] - # Using −dt/2 reverses the forward evolution so the bond does not - # get double-counted when A[i+1] is updated in the next iteration. - L_new = left_envs[i + 1] - R_cur = right_envs[i + 1] - _cn, _cr2 = chi_new, chi_r - - def keff_mv(v, _L=L_new, _R=R_cur, _cn=_cn, _cr2=_cr2): - return _apply_keff(_L, _R, v.reshape(_cn, _cr2), xp=xp).flatten() - - C = _arnoldi_expm(keff_mv, C.flatten(), -dt / 2.0, - m=self.krylov_dim).reshape(chi_new, chi_r) - - # 5. Absorb back-propagated C into A[i+1] - A_next = self._mps[i + 1] # (chi_r, chi_rr, d) - chi_r_n, chi_rr, d_n = A_next.shape - A_next_mat = A_next.reshape(chi_r_n, chi_rr * d_n) - self._mps[i + 1] = (C @ A_next_mat).reshape(chi_new, chi_rr, d_n) - - def _sweep_rl(self, left_envs: List, right_envs: List, dt: float) -> None: - """ - Right-to-left half-sweep with back-propagation (1-site TDVP). - - At each site i: - 1. Forward Krylov: A[i] ← exp(H_eff[i] × dt/2) A[i] - 2. SVD → right-canonical A[i] + center C - 3. Update right environment R[i] - 4. Back-propagation: C ← exp(−K_eff × dt/2) C (undo double-counting) - 5. Absorb C into A[i-1] - """ - xp = self._xp - n = len(self._mps) - mpo = self._mpo - - for i in range(n - 1, 0, -1): - A = self._mps[i] - W_mpo = xp.array(mpo[i], dtype=complex) - L = left_envs[i] - R = right_envs[i + 1] - chi_l, chi_r, d = A.shape - - # 1. Krylov: exp(H_eff × dt/2) @ A - def heff_mv(v, _L=L, _W=W_mpo, _R=R, - _cl=chi_l, _cr=chi_r, _d=d): - return _apply_heff(_L, _W, _R, - v.reshape(_cl, _cr, _d), xp=xp).flatten() - - A_new = _arnoldi_expm(heff_mv, A.flatten(), dt / 2.0, - m=self.krylov_dim).reshape(chi_l, chi_r, d) - - # 2. SVD → right-canonical + center C - A_mat = A_new.reshape(chi_l, chi_r * d) - U, s, Vh = xp.linalg.svd(A_mat, full_matrices=False) - chi_new = min(self.chi, len(s)) - U, s, Vh = U[:, :chi_new], s[:chi_new], Vh[:chi_new, :] - - # Right-canonical tensor: (chi_new, chi_r, d) - self._mps[i] = Vh.reshape(chi_new, chi_r, d) - - C = U @ xp.diag(s) # center: (chi_l, chi_new) - - # 3. Update right environment R[i] from new right-canonical A[i] - right_envs[i] = _update_right_env(R, self._mps[i], W_mpo, xp=xp) - - # 4. Back-propagation: C ← exp(−K_eff × dt/2) C - # K_eff[a,c;b,d] = Σ_m L[i][a,m,b] · R[i][c,m,d] - L_cur = left_envs[i] - R_new = right_envs[i] - _cl2, _cn = chi_l, chi_new - - def keff_mv(v, _L=L_cur, _R=R_new, _cl2=_cl2, _cn=_cn): - return _apply_keff(_L, _R, v.reshape(_cl2, _cn), xp=xp).flatten() - - C = _arnoldi_expm(keff_mv, C.flatten(), -dt / 2.0, - m=self.krylov_dim).reshape(chi_l, chi_new) - - # 5. Absorb back-propagated C into A[i-1] - A_prev = self._mps[i - 1] # (chi_ll, chi_l, d) - chi_ll, chi_l_p, d_p = A_prev.shape - # Must transpose before reshape to isolate the bond dimension: - # (chi_ll, chi_l, d) → transpose(0,2,1) → (chi_ll, d, chi_l) - # → reshape(chi_ll*d, chi_l) — correct matrix with chi_l as column - A_prev_mat = A_prev.transpose(0, 2, 1).reshape(chi_ll * d_p, chi_l_p) - # (chi_ll*d, chi_l) @ (chi_l, chi_new) → (chi_ll*d, chi_new) - self._mps[i - 1] = (A_prev_mat @ C).reshape(chi_ll, d_p, chi_new).transpose(0, 2, 1) - - # ── 2-site TDVP sweep ───────────────────────────────────────────────────── - - def _lr_2site_pass(self, left_envs: List, right_envs: List, dt: float) -> None: - """ - Single left-to-right pass of 2-site bond updates. - - Internal building block: not second-order on its own. Call via - _sweep_lr_2site which composes two half-dt passes for O(dt²) accuracy. - """ - xp = self._xp - n = len(self._mps) - mpo = self._mpo - - for i in range(n - 1): - A_i = self._mps[i] # (chi_l, chi_r, d) - A_i1 = self._mps[i + 1] # (chi_r, chi_rr, d) - W0 = xp.array(mpo[i], dtype=complex) - W1 = xp.array(mpo[i + 1], dtype=complex) - L_i = left_envs[i] - R_i2 = right_envs[i + 2] if (i + 2) <= n else right_envs[n] - - chi_l, chi_r, d = A_i.shape - _, chi_rr, _ = A_i1.shape - - # 1. Form Θ: A[i] ⊗ A[i+1] → (chi_l, d, d, chi_rr) - Theta = xp.einsum("abs,bct->astc", A_i, A_i1) - - # 2. Krylov: exp(H_eff_2site × dt) @ Theta - _cl, _d0, _d1, _cr = Theta.shape - - def heff2_mv(v, - _L=L_i, _W0=W0, _W1=W1, _R=R_i2, - _cl=_cl, _d0=_d0, _d1=_d1, _cr=_cr): - T = v.reshape(_cl, _d0, _d1, _cr) - return _apply_heff_2site(_L, _W0, _W1, _R, T, xp=xp).flatten() - - Theta_new = _arnoldi_expm( - heff2_mv, Theta.flatten(), dt, - m=self.krylov_dim - ).reshape(_cl, _d0, _d1, _cr) - - # 3. SVD → (chi_l*d, d*chi_rr), truncate to chi_max - mat = Theta_new.reshape(chi_l * d, d * chi_rr) - U, s, Vh = xp.linalg.svd(mat, full_matrices=False) - chi_new = min(self.chi, len(s)) - U, s, Vh = U[:, :chi_new], s[:chi_new], Vh[:chi_new, :] - - # Left-canonical A[i]: (chi_l, d, chi_new) → (chi_l, chi_new, d) - self._mps[i] = U.reshape(chi_l, d, chi_new).transpose(0, 2, 1) - - # A[i+1] absorbs singular values: (chi_new, d, chi_rr) → (chi_new, chi_rr, d) - SV = (xp.diag(s) @ Vh).reshape(chi_new, d, chi_rr) - self._mps[i + 1] = SV.transpose(0, 2, 1) - - # 4. Update left environment for site i+1 - left_envs[i + 1] = _update_left_env(L_i, self._mps[i], W0, xp=xp) - - # 5. Back-propagate A[i+1] by exp(-H_eff_1site × dt) to undo the - # W_{i+1} pre-evolution baked into the right SVD piece. Without - # this, bond (i+1, i+2) would apply W_{i+1} a second time, causing - # each LR pass to advance the state by (n-1)×dt instead of dt. - if i < n - 2: - L_bp = left_envs[i + 1] - # W1 and R_i2 already loaded above; reuse them. - cl_bp, cr_bp, d_bp = self._mps[i + 1].shape - - def heff1_bp_mv(v, _L=L_bp, _W=W1, _R=R_i2, - _cl=cl_bp, _cr=cr_bp, _d=d_bp): - return _apply_heff(_L, _W, _R, - v.reshape(_cl, _cr, _d), xp=xp).flatten() - - self._mps[i + 1] = _arnoldi_expm( - heff1_bp_mv, self._mps[i + 1].flatten(), -dt, - m=self.krylov_dim - ).reshape(cl_bp, cr_bp, d_bp) - - def _sweep_lr_2site(self, left_envs: List, right_envs: List, dt: float) -> None: - """ - Second-order left-to-right 2-site TDVP sweep via midpoint rule. - - Two half-dt passes with right-environment rebuild between them: - - Pass 1 (dt/2): LR bond updates using provided environments. - Produces an O(dt)-accurate midpoint MPS state. - Env rebuild: Right environments recomputed from midpoint MPS. - Pass 2 (dt/2): LR bond updates using midpoint environments. - Advances from midpoint to t+dt with O(dt²) accuracy. - - This is equivalent to the implicit midpoint / Störmer-Verlet scheme - adapted to the LR-only constraint (RL sweep omitted to avoid FSM MPO - double-counting — see run_2site docstring §7). - - Error per step: O(dt³) local / O(dt²) global. - """ - n = len(self._mps) - xp = self._xp - - # Pass 1: half-step with caller-supplied environments - self._lr_2site_pass(left_envs, right_envs, dt / 2) - - # Rebuild right environments from the midpoint MPS - right_envs_mid = self._build_right_envs() - L0, _ = _init_boundary_envs(n, xp=xp) - left_envs_mid = [None] * (n + 1) - left_envs_mid[0] = L0 - - # Pass 2: second half-step with midpoint environments - self._lr_2site_pass(left_envs_mid, right_envs_mid, dt / 2) - - # Propagate final left environments back to caller - for i in range(n + 1): - left_envs[i] = left_envs_mid[i] - - def _sweep_rl_2site(self, left_envs: List, right_envs: List, dt: float) -> None: - """ - Right-to-left 2-site TDVP sweep (no back-propagation). - - Mirror of _sweep_lr_2site going right→left. Pair with _sweep_lr_2site - (each using dt/2) for a symmetric second-order integrator. - """ - xp = self._xp - n = len(self._mps) - mpo = self._mpo - - for i in range(n - 2, -1, -1): - A_i = self._mps[i] # (chi_l, chi_r, d) - A_i1 = self._mps[i + 1] # (chi_r, chi_rr, d) - W0 = xp.array(mpo[i], dtype=complex) - W1 = xp.array(mpo[i + 1], dtype=complex) - L_i = left_envs[i] - R_i2 = right_envs[i + 2] if (i + 2) <= n else right_envs[n] - - chi_l, chi_r, d = A_i.shape - _, chi_rr, _ = A_i1.shape - - # 1. Form Θ: A[i] ⊗ A[i+1] → (chi_l, d, d, chi_rr) - Theta = xp.einsum("abs,bct->astc", A_i, A_i1) - - # 2. Krylov: exp(H_eff_2site × dt) @ Theta - _cl, _d0, _d1, _cr = Theta.shape - - def heff2_mv(v, - _L=L_i, _W0=W0, _W1=W1, _R=R_i2, - _cl=_cl, _d0=_d0, _d1=_d1, _cr=_cr): - T = v.reshape(_cl, _d0, _d1, _cr) - return _apply_heff_2site(_L, _W0, _W1, _R, T, xp=xp).flatten() - - Theta_new = _arnoldi_expm( - heff2_mv, Theta.flatten(), dt, - m=self.krylov_dim - ).reshape(_cl, _d0, _d1, _cr) - - # 3. SVD → (chi_l*d, d*chi_rr), truncate to chi_max - mat = Theta_new.reshape(chi_l * d, d * chi_rr) - U, s, Vh = xp.linalg.svd(mat, full_matrices=False) - chi_new = min(self.chi, len(s)) - U, s, Vh = U[:, :chi_new], s[:chi_new], Vh[:chi_new, :] - - # A[i] absorbs singular values (right-canonical side): - # SV shape (chi_l, d, chi_new) → (chi_l, chi_new, d) - SV = (U @ xp.diag(s)).reshape(chi_l, d, chi_new) - self._mps[i] = SV.transpose(0, 2, 1) - - # Right-canonical A[i+1]: (chi_new, d, chi_rr) → (chi_new, chi_rr, d) - self._mps[i + 1] = Vh.reshape(chi_new, d, chi_rr).transpose(0, 2, 1) - - # 4. Update right environment for site i+1 - right_envs[i + 1] = _update_right_env(R_i2, self._mps[i + 1], W1, xp=xp) - - # ── main integration loop ───────────────────────────────────────────────── - - def run( - self, - t_max_us: float = 10.0, - n_steps: int = 200, - ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: - """ - TDVP integration on GPU. Returns (t_us, P_S, trace). - """ - if self._mpo is None: - self._build() - - cfg = self.config - n = cfg.n_sites - t_eval = np.linspace(0.0, t_max_us, n_steps) - dt = t_eval[1] - t_eval[0] - - logger.info(f"[TDVP-GPU] {cfg.name}: chi={self.chi}, " - f"{n_steps} steps, dt={dt:.4f} us...") - t0 = time.perf_counter() - - P_S = np.zeros(n_steps) - trace = np.zeros(n_steps) - P_S[0], trace[0] = self._observables() - - # Initial right environments - right_envs = self._build_right_envs() - L0, _ = _init_boundary_envs(n, xp=self._xp) - left_envs = [None] * (n + 1) - left_envs[0] = L0 - - log_every = max(1, n_steps // 10) - - for step in range(1, n_steps): - self._sweep_lr(left_envs, right_envs, dt) - self._sweep_rl(left_envs, right_envs, dt) - P_S[step], trace[step] = self._observables() - - if step % log_every == 0: - elapsed = time.perf_counter() - t0 - logger.info(f" step {step:4d}/{n_steps} — " - f"P_S={P_S[step]:.4f}, tr={trace[step]:.4f}, " - f"{elapsed:.1f}s") - - wall = time.perf_counter() - t0 - logger.info(f"[TDVP-GPU] done: {wall:.2f}s") - _log({"event": "tdvp_integrate", "system": cfg.name, - "n_sites": n, "chi": self.chi, "n_steps": n_steps, - "t_max_us": t_max_us, "wall_time_s": round(wall, 3)}) - - return t_eval, P_S, trace - - def run_2site( - self, - t_max_us: float = 10.0, - n_steps: int = 200, - ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: - """ - 2-site TDVP integration (LR-only sweep). Returns (t_us, P_S, trace). - - **Algorithm:** Left-to-right sweep only, rebuilding right environments each - step from the current MPS. The LR sweep naturally works with this MPO - structure because the FSM accumulates Liouvillian contributions left→right: - only the last bond (n-2, n-1) sees the full accumulated operator and accounts - for the physical trace decay. Earlier bonds have near-zero trace-decreasing - eigenvalues in their 2-site H_eff. - - **Why not LR+RL?** The RL sweep causes catastrophic trace double-counting: - each bond in the right-to-left direction independently sees the full - accumulated Liouvillian (carried backward through right environments), so - every bond applies the full recombination decay, giving N× overcounting. - This is a structural property of the left-to-right FSM MPO, not a bug. - - **Accuracy:** Second-order in dt (O(dt²) global error). Each call to - _sweep_lr_2site performs two half-dt passes with a midpoint environment - rebuild, achieving the same order as Strang splitting without an RL sweep. - Use n_steps ≥ 200 (dt ≤ 0.05 μs) for RMSE < 1e-3 vs ExactSolver. - - At each step: - 1. Rebuild right environments R[0..n] from current MPS. - 2. LR 2-site sweep: for each bond (i, i+1), - Θ[i] ← exp(H_eff_2site[i,i+1] × dt) Θ[i] - then SVD-truncate to χ_max. - 3. Record P_S and trace. - """ - if self._mpo is None: - self._build() - - cfg = self.config - n = cfg.n_sites - t_eval = np.linspace(0.0, t_max_us, n_steps) - dt = t_eval[1] - t_eval[0] - - logger.info(f"[TDVP2-GPU] {cfg.name}: chi={self.chi}, " - f"{n_steps} steps, dt={dt:.4f} us (LR-only)...") - t0 = time.perf_counter() - - P_S = np.zeros(n_steps) - trace = np.zeros(n_steps) - P_S[0], trace[0] = self._observables() - - log_every = max(1, n_steps // 10) - - for step in range(1, n_steps): - # Rebuild right environments from current MPS at the start of each step. - # The LR sweep builds left_envs incrementally as it processes each bond. - right_envs = self._build_right_envs() - L0, _ = _init_boundary_envs(n, xp=self._xp) - left_envs = [None] * (n + 1) - left_envs[0] = L0 - self._sweep_lr_2site(left_envs, right_envs, dt) - - P_S[step], trace[step] = self._observables() - - if step % log_every == 0: - elapsed = time.perf_counter() - t0 - logger.info(f" step {step:4d}/{n_steps} — " - f"P_S={P_S[step]:.4f}, tr={trace[step]:.4f}, " - f"{elapsed:.1f}s") - - wall = time.perf_counter() - t0 - logger.info(f"[TDVP2-GPU] done: {wall:.2f}s") - _log({"event": "tdvp2_integrate", "system": cfg.name, - "n_sites": n, "chi": self.chi, "n_steps": n_steps, - "t_max_us": t_max_us, "wall_time_s": round(wall, 3)}) - - return t_eval, P_S, trace - - -# ═══════════════════════════════════════════════════════════════════════════════ -# §8 Backwards-compatible CuTensorNetEngine (Phase A stub, kept for API compat) -# ═══════════════════════════════════════════════════════════════════════════════ - -class CuTensorNetEngine: - """ - Legacy stub kept for API compatibility. - Use CupyKrylovSolver (N≤14) or CuTDVPSolver (N≤62) instead. - """ - - def __init__(self, chi: int = 2500, device_id: int = 0): - self.chi = chi - self.device_id = device_id - logger.info("CuTensorNetEngine: use CuTDVPSolver for production runs.") - - def initialize(self) -> None: - if not HAS_CUPY: - raise ImportError("pip install cupy-cuda12x") - cp.cuda.Device(self.device_id).use() - logger.info(f"Device {self.device_id}: {cp.cuda.Device().name()}") - logger.info(f"cuquantum available: {HAS_CUTN}") - - def tdvp_site_update(self, A_site, H_eff, dt_us: float): - raise NotImplementedError("Use CuTDVPSolver which implements this properly.") - - def benchmark_matvec(self, chi: int, d: int = 2, - n_warmup: int = 3, n_runs: int = 20) -> float: - """Benchmark raw GPU matmul at size (χ²d)×(χ²d).""" - if not HAS_CUPY: - raise ImportError("pip install cupy-cuda12x") - size = chi ** 2 * d - A = cp.random.randn(size, size).astype(complex) - v = cp.random.randn(size).astype(complex) - for _ in range(n_warmup): - _ = A @ v - cp.cuda.Device().synchronize() - times = [] - for _ in range(n_runs): - cp.cuda.Device().synchronize() - t0 = time.perf_counter() - _ = A @ v - cp.cuda.Device().synchronize() - times.append(time.perf_counter() - t0) - median_ms = float(np.median(times)) * 1e3 - logger.info(f"matvec χ={chi}: {median_ms:.2f} ms/op") - _log({"event": "gpu_matvec_benchmark", "chi": chi, "size": size, - "median_ms": round(median_ms, 3)}) - return median_ms - - -# ═══════════════════════════════════════════════════════════════════════════════ -# §9 VRAM estimator -# ═══════════════════════════════════════════════════════════════════════════════ - -def estimate_vram_gb(n_sites: int, chi: int, d: int = 4) -> dict: - """Estimate GPU VRAM requirements for CuTDVPSolver.""" - B = 16 # bytes per complex128 - mps_bytes = n_sites * chi ** 2 * d * B - env_bytes = 2 * n_sites * chi ** 2 * _W * B # L + R environments - mpo_bytes = n_sites * _W ** 2 * d ** 2 * B # negligible - total = (mps_bytes + env_bytes + mpo_bytes) / 1e9 - return { - "n_sites": n_sites, "chi": chi, - "mps_gb": round(mps_bytes / 1e9, 3), - "env_gb": round(env_bytes / 1e9, 3), - "total_gb": round(total, 3), - "fits_16gb": total < 14.0, - } - - -def print_vram_table() -> None: - """Print VRAM table for key N × χ combinations.""" - print(f"\n{'N':>4} {'χ':>6} {'MPS GB':>8} {'Env GB':>8} {'Total':>8} {'Fits 16GB':>10}") - print("─" * 58) - for n in [10, 20, 32, 40, 62]: - for chi in [64, 256, 512, 1024, 2500]: - r = estimate_vram_gb(n, chi) - flag = "✓" if r["fits_16gb"] else "✗" - print(f"{n:4d} {chi:6d} {r['mps_gb']:8.3f} " - f"{r['env_gb']:8.3f} {r['total_gb']:8.3f} {flag:>10}") diff --git a/modules/cycle_project/ARCHITECTURE.md b/modules/cycle_project/ARCHITECTURE.md deleted file mode 100644 index 89975f7..0000000 --- a/modules/cycle_project/ARCHITECTURE.md +++ /dev/null @@ -1,244 +0,0 @@ -# cycle_project — Architecture - -> Investigates the crustal displacement / catastrophic cycle hypothesis -> (Younger Dryas ~12,900 BP; Gothenburg paleomagnetic excursion; Be-10 production -> spikes; global flood-myth convergence) using public geological data and -> GPU-accelerated ML. - ---- - -## Directory Tree - -``` -cycle_project/ -├── ARCHITECTURE.md -├── pyproject.toml # Python workspace (uv / pip) -├── Cargo.toml # Rust workspace (field_coherence) -├── configs/ -│ └── data_sources.yaml # All URLs + checksums -├── data/ -│ ├── raw/ # Downloaded originals (gitignored) -│ └── processed/ # Parquet + PNG outputs -├── notebooks/ -│ ├── 01_proxy_eda.ipynb -│ ├── 02_gnn_training.ipynb -│ └── 03_lbm_validation.ipynb -├── scripts/ -│ ├── run_fetch.sh -│ ├── run_gnn.sh -│ └── ingest_myths.sh -└── src/ - ├── cycle_detect/ # MODULE 1 - │ ├── __init__.py - │ ├── fetch_data.py # ← delivered here - │ └── gnn_prototype.py # ← delivered here - ├── pole_shift_sim/ # MODULE 2 - │ ├── __init__.py - │ ├── lbm_lithosphere.py # Lattice Boltzmann sim (CuPy/CUDA) - │ ├── gauge_coupling.py # U(1) asthenosphere viscosity field - │ └── visualize.py # PyVista 3-D render pipeline - ├── myth_rag/ # MODULE 3 - │ ├── __init__.py - │ ├── ingest.py # Corpus → ChromaDB - │ ├── rag_agent.py # RAG query + geological correlation - │ └── corpus/ # Plain-text myth sources (CC/PD) - └── field_coherence/ # MODULE 4 (Rust) - ├── Cargo.toml - └── src/ - ├── main.rs # egui dashboard entry point - ├── noaa_api.rs # NOAA Space Weather API poller - ├── pole_tracker.rs # Geomagnetic north acceleration - └── cosmic_ray.rs # Oulu / IZMIRAN cosmic-ray feed -``` - ---- - -## Module 1 — CYCLE_DETECT - -**Goal:** GNN over 5 paleoclimate proxies → detect synchronous anomaly windows -indicative of a recurrent ~11,500–12,900 yr catastrophe cycle. - -| Component | Choice | Justification | -|-----------|--------|---------------| -| Data fetch | `requests` + `urllib` + `tqdm` | No external dep; FTP redirect to HTTPS at NCEI | -| DataFrame | `pandas` + `numpy` | Standard; parquet via `pyarrow` | -| Graph ML | `torch_geometric` (PyG) | Best PyTorch-native temporal GNN ecosystem | -| GNN arch | `SAGEConv` autoencoder | Inductive; handles variable-density graphs | -| GPU | PyTorch `.to(device)` sm_120 | RTX 5060 Ti — needs `torch>=2.5` nightly for sm_120 | -| Anomaly | Reconstruction MSE per window | Unsupervised; no labeled YD events needed | -| Storage | Apache Parquet | Columnar; fast read for 150k rows × 5 proxies | - -**Proxies ingested:** -1. GISP2 δ¹⁸O — Greenland temperature proxy (Alley 2000) -2. Vostok ΔT — Antarctic temperature (deuterium) -3. Vostok CO₂ — atmospheric greenhouse forcing -4. GRIP Be-10 flux — cosmic-ray / solar modulation proxy -5. Sint-2000 VADM — geomagnetic dipole moment stack - ---- - -## Module 2 — POLE_SHIFT_SIM - -**Goal:** 3-D Lattice Boltzmann simulation of lithospheric slab sliding over -a low-viscosity asthenosphere, driven by an asymmetric ice-load torque. -U(1) gauge coupling encodes the asthenosphere as a viscosity field. - -| Component | Choice | Justification | -|-----------|--------|---------------| -| LBM kernel | CuPy (D3Q19) | Direct CUDA array ops; matches QUANTUM_LAB CUDA style | -| Gauge field | NumPy/CuPy U(1) link matrices | Re-uses HMC infrastructure from QUANTUM_LAB | -| Geometry | Spherical lat-lon grid (512×256×32) | Sufficient for continental-scale flow | -| Viz | PyVista + VTK | Interactive 3-D; export to .vtp for ParaView | -| Validation | Compare to Steinberger & O'Connell 1997 | Asthenosphere flow benchmark | - -**Physics encoded:** -- Incompressible Navier-Stokes (LBM BGK collision) -- Yield-stress rheology for lithosphere (Bingham fluid) -- U(1) gauge field: asthenosphere viscosity η(x) as a scalar Wilson field -- Ice-sheet gravitational torque as external body force - ---- - -## Module 3 — MYTH_RAG - -**Goal:** ChromaDB vector store of primary flood/catastrophe myth texts -(Sumerian Atrahasis, Egyptian Papyrus of Ipuwer, Rigveda Manu, Norse Fimbulwinter, -Aztec Suns cosmogony, Mayan Popol Vuh, Greek Deucalion) + LangChain RAG agent -that queries geological events and retrieves thematically aligned myth passages. - -| Component | Choice | Justification | -|-----------|--------|---------------| -| Vector DB | ChromaDB (local, persistent) | Local, sovereign, no cloud | -| Embeddings | `sentence-transformers` (all-mpnet-base-v2) | Best semantic similarity for ancient texts | -| LLM backend | Ollama (llama3/mistral local) | Local-first; no API dependency | -| RAG orchestration | LangChain (minimal) | Familiar; swap for raw similarity search if needed | -| Corpus | Plain-text from Project Gutenberg + ETCSL | Public domain; licensed for research | -| Output | Markdown report: myth ↔ geological event correlation table | | - -**Correlation targets:** -- 12,900 BP → YD onset: abrupt cooling, Atrahasis flood narrative -- 11,700 BP → YD termination: civilizational restart myths -- 74,000 BP → Toba event: population bottleneck in Hindu cosmogony -- ~41,000 BP → Laschamps excursion: aurora/sky-fire myths - ---- - -## Module 4 — FIELD_COHERENCE_MONITOR - -**Goal:** Real-time Rust + egui dashboard showing geomagnetic pole acceleration, -solar flux (F10.7), and cosmic-ray flux from Oulu Neutron Monitor. - -| Component | Choice | Justification | -|-----------|--------|---------------| -| GUI | `egui` + `eframe` | Immediate-mode; no X11 dep; renders on bare metal | -| HTTP polling | `reqwest` (tokio async) | Non-blocking; update every 60s | -| Geomagnetic pole | NOAA NCEI WMM API | `https://www.ngdc.noaa.gov/geomag-web/calculators/calculateDeclination` | -| Solar activity | NOAA SWPC JSON | `https://services.swpc.noaa.gov/json/f107_cm_flux.json` | -| Cosmic rays | Oulu NM FTP | `http://cosmicrays.oulu.fi/` (hourly count rates) | -| Plotting | `egui_plot` | Real-time streaming chart; ring buffer 24h | -| Persistence | SQLite via `rusqlite` | Local timeseries storage; no cloud | - -**Alerts:** Configurable threshold triggers (Ω > 2σ from 30-day mean) -emit desktop notification via `notify-rust`. - ---- - -## Data Sources (exact paths) - -```yaml -# configs/data_sources.yaml - -gisp2_d18o: - url: "https://www.ncei.noaa.gov/pub/data/paleo/icecore/greenland/summit/gisp2/isotopes/gisp2_d18o_accum_alley2000.txt" - description: "GISP2 δ18O + accumulation rate, Alley 2000, 0–110,000 yr BP" - verified: true - -vostok_deuterium: - url: "https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/vostok/deutnat.txt" - description: "Vostok deuterium + ΔT, Petit et al. 1999, 0–420,000 yr BP" - verified: true - -vostok_co2: - url: "https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/vostok/co2nat.txt" - description: "Vostok CO2, Petit et al. 1999, 0–420,000 yr BP" - verified: true - -grip_be10: - url: "https://www.ncei.noaa.gov/pub/data/paleo/icecore/greenland/summit/grip/beryllium/grip_be10_muscheler2004.txt" - description: "GRIP Be-10 flux, Muscheler et al. 2004, 0–110,000 yr BP" - verified: false # ⚠️ browse dir to confirm filename - fallback_dir: "https://www.ncei.noaa.gov/pub/data/paleo/icecore/greenland/summit/grip/beryllium/" - -sint2000: - url: "https://www.ngdc.noaa.gov/geomag/paleo_mag_datasets/Sint-2000.txt" - description: "Sint-2000 VADM stack, Valet et al. 2005, 0–2,000,000 yr BP at 1-ka res" - verified: false # ⚠️ fallback: PANGAEA DOI:10.1594/PANGAEA.186810 - fallback_pangaea: "https://doi.pangaea.de/10.1594/PANGAEA.186810" - -# Additional (Phase 2+) -sea_level_spratt: - url: "https://www.ncei.noaa.gov/pub/data/paleo/contributions_by_author/spratt2016/spratt2016.txt" - description: "Sea level stack, Spratt & Lisiecki 2016, 0–800,000 yr BP" - verified: true - -lisiecki_lr04: - url: "https://www.ncei.noaa.gov/pub/data/paleo/contributions_by_author/lisiecki2005/lisiecki2005.txt" - description: "LR04 benthic δ18O stack, Lisiecki & Raymo 2005" - verified: true -``` - ---- - -## 4-Phase Implementation Roadmap - -### Phase 1 — Data Foundation (2–3 weeks) -- [x] Project scaffold + `pyproject.toml` -- [ ] `fetch_data.py`: download all proxies, clean, resample 100-yr grid -- [ ] Manual verification of ⚠️ URLs; fallback downloads if needed -- [ ] `gnn_prototype.py`: baseline GraphSAGE autoencoder, anomaly timeline -- [ ] Notebook 01: EDA — confirm YD spike at 12,900 BP in all proxies -- Deliverable: `data/processed/overview.png` showing 5-proxy alignment - -### Phase 2 — GNN + LBM Core (4–6 weeks) -- [ ] GNN: hyperparameter sweep (window size 1k–10k yr, edge threshold) -- [ ] GNN: compare anomaly scores to known events (YD, 8.2ka, Laschamps) -- [ ] LBM: D3Q19 kernel in CuPy; benchmark on RTX 5060 Ti -- [ ] LBM: U(1) viscosity coupling; polar ice-torque forcing -- [ ] `visualize.py`: PyVista animation of slab displacement -- Deliverable: LBM simulation video + GNN anomaly heatmap - -### Phase 3 — RAG + Real-Time Monitor (3–4 weeks) -- [ ] Corpus ingestion: 7 myth traditions → ChromaDB -- [ ] RAG agent: geological event → myth passage retrieval -- [ ] Correlation report: markdown table (event × myth × similarity score) -- [ ] Rust `field_coherence`: egui dashboard compiling on Arch bare metal -- [ ] Live data feeds: NOAA SWPC + Oulu NM + NCEI WMM -- Deliverable: Running dashboard + myth correlation report - -### Phase 4 — Integration + Publication (2–3 weeks) -- [ ] Cross-module: GNN anomaly dates → RAG queries → myth passages -- [ ] Statistical validation: bootstrap significance test on synchrony scores -- [ ] Full 150,000-yr multi-proxy animation (PyVista + matplotlib) -- [ ] LaTeX report / preprint draft -- [ ] Open dataset release (processed parquets on Zenodo) -- Deliverable: Reproducible research archive - -**Total estimated time: 11–16 weeks** (solo researcher, part-time alongside QUANTUM_LAB) - ---- - -## GPU Notes (RTX 5060 Ti — sm_120) - -```toml -# Requires torch >= 2.5 with sm_120 support -# Install nightly if stable doesn't build for Ada Lovelace successor arch: -# pip install --pre torch --index-url https://download.pytorch.org/whl/nightly/cu124 - -# CuPy: build from source or use cupy-cuda12x ≥ 13.x -# CUDA-Q / PennyLane: sm_120 PTX should JIT-compile fine -``` - -All CUDA kernels in POLE_SHIFT_SIM target `sm_120` via `nvcc -arch=sm_120`. -PyTorch ops fall back to CPU automatically if sm_120 is not yet supported -in the installed build (guarded by `device = torch.device("cuda" if -torch.cuda.is_available() else "cpu")`). diff --git a/modules/cycle_project/cycle_anomaly_map.png b/modules/cycle_project/cycle_anomaly_map.png index 4ffe6f8..e5e9490 100644 Binary files a/modules/cycle_project/cycle_anomaly_map.png and b/modules/cycle_project/cycle_anomaly_map.png differ diff --git a/modules/cycle_project/cycle_overview.png b/modules/cycle_project/cycle_overview.png index 04cd487..c423b2e 100644 Binary files a/modules/cycle_project/cycle_overview.png and b/modules/cycle_project/cycle_overview.png differ diff --git a/modules/cycle_project/cycle_score_dist.png b/modules/cycle_project/cycle_score_dist.png index 7ed81b7..b62eec6 100644 Binary files a/modules/cycle_project/cycle_score_dist.png and b/modules/cycle_project/cycle_score_dist.png differ diff --git a/modules/cycle_project/data/processed/anomaly_map.png b/modules/cycle_project/data/processed/anomaly_map.png index 31fafe4..9b7d73b 100644 Binary files a/modules/cycle_project/data/processed/anomaly_map.png and b/modules/cycle_project/data/processed/anomaly_map.png differ diff --git a/modules/cycle_project/data/processed/cycle_project_overview.png b/modules/cycle_project/data/processed/cycle_project_overview.png index 3148937..b971dd5 100644 Binary files a/modules/cycle_project/data/processed/cycle_project_overview.png and b/modules/cycle_project/data/processed/cycle_project_overview.png differ diff --git a/modules/cycle_project/data/processed/gnn_autoencoder.pt b/modules/cycle_project/data/processed/gnn_autoencoder.pt deleted file mode 100644 index 22df030..0000000 Binary files a/modules/cycle_project/data/processed/gnn_autoencoder.pt and /dev/null differ diff --git a/modules/cycle_project/data/processed/myth_heatmap.png b/modules/cycle_project/data/processed/myth_heatmap.png index f825896..4cdc4ca 100644 Binary files a/modules/cycle_project/data/processed/myth_heatmap.png and b/modules/cycle_project/data/processed/myth_heatmap.png differ diff --git a/modules/cycle_project/data/processed/myth_timeline.png b/modules/cycle_project/data/processed/myth_timeline.png index f378a84..d2d09bd 100644 Binary files a/modules/cycle_project/data/processed/myth_timeline.png and b/modules/cycle_project/data/processed/myth_timeline.png differ diff --git a/modules/cycle_project/data/processed/score_distribution.png b/modules/cycle_project/data/processed/score_distribution.png index 7ed81b7..b62eec6 100644 Binary files a/modules/cycle_project/data/processed/score_distribution.png and b/modules/cycle_project/data/processed/score_distribution.png differ diff --git a/modules/cycle_project/data/processed/training_loss.png b/modules/cycle_project/data/processed/training_loss.png index 5029e75..7b906a4 100644 Binary files a/modules/cycle_project/data/processed/training_loss.png and b/modules/cycle_project/data/processed/training_loss.png differ diff --git a/modules/cycle_project/data/processed/vadm_forecast.png b/modules/cycle_project/data/processed/vadm_forecast.png index 0eab393..6b5496e 100644 Binary files a/modules/cycle_project/data/processed/vadm_forecast.png and b/modules/cycle_project/data/processed/vadm_forecast.png differ diff --git a/modules/cycle_project/data/processed/vadm_lstm_forecast.png b/modules/cycle_project/data/processed/vadm_lstm_forecast.png index 4ef3d86..df4e935 100644 Binary files a/modules/cycle_project/data/processed/vadm_lstm_forecast.png and b/modules/cycle_project/data/processed/vadm_lstm_forecast.png differ diff --git a/modules/cycle_project/data/processed/vadm_spectrum.png b/modules/cycle_project/data/processed/vadm_spectrum.png index 26aef45..591c80b 100644 Binary files a/modules/cycle_project/data/processed/vadm_spectrum.png and b/modules/cycle_project/data/processed/vadm_spectrum.png differ diff --git a/modules/cycle_project/data/processed/viz/displacement_map.png b/modules/cycle_project/data/processed/viz/displacement_map.png index e1845d3..1dc77fc 100644 Binary files a/modules/cycle_project/data/processed/viz/displacement_map.png and b/modules/cycle_project/data/processed/viz/displacement_map.png differ diff --git a/modules/cycle_project/data/processed/viz/time_evolution.png b/modules/cycle_project/data/processed/viz/time_evolution.png index de7e179..ed4cf72 100644 Binary files a/modules/cycle_project/data/processed/viz/time_evolution.png and b/modules/cycle_project/data/processed/viz/time_evolution.png differ diff --git a/modules/cycle_project/data/processed/viz/vadm_coupled_comparison.png b/modules/cycle_project/data/processed/viz/vadm_coupled_comparison.png index 3a00d6d..67589a6 100644 Binary files a/modules/cycle_project/data/processed/viz/vadm_coupled_comparison.png and b/modules/cycle_project/data/processed/viz/vadm_coupled_comparison.png differ diff --git a/modules/cycle_project/pole_shift_displacement.png b/modules/cycle_project/pole_shift_displacement.png index e1845d3..1dc77fc 100644 Binary files a/modules/cycle_project/pole_shift_displacement.png and b/modules/cycle_project/pole_shift_displacement.png differ diff --git a/modules/cycle_project/pole_shift_evolution.png b/modules/cycle_project/pole_shift_evolution.png index de7e179..ed4cf72 100644 Binary files a/modules/cycle_project/pole_shift_evolution.png and b/modules/cycle_project/pole_shift_evolution.png differ diff --git a/modules/cycle_project/vadm_coupled_comparison.png b/modules/cycle_project/vadm_coupled_comparison.png index 06fe489..6d88621 100644 Binary files a/modules/cycle_project/vadm_coupled_comparison.png and b/modules/cycle_project/vadm_coupled_comparison.png differ diff --git a/modules/quantum_lab/P1_LHCB/anomaly_scores.png b/modules/quantum_lab/P1_LHCB/anomaly_scores.png index 26b0111..ac80ffa 100644 Binary files a/modules/quantum_lab/P1_LHCB/anomaly_scores.png and b/modules/quantum_lab/P1_LHCB/anomaly_scores.png differ diff --git a/modules/quantum_lab/P1_LHCB/roc_auc_curve.png b/modules/quantum_lab/P1_LHCB/roc_auc_curve.png index 9d9d585..b955bad 100644 Binary files a/modules/quantum_lab/P1_LHCB/roc_auc_curve.png and b/modules/quantum_lab/P1_LHCB/roc_auc_curve.png differ diff --git a/modules/quantum_lab/P2_WBOSON/unfolding_comparison.png b/modules/quantum_lab/P2_WBOSON/unfolding_comparison.png index bc28458..b0f0c14 100644 Binary files a/modules/quantum_lab/P2_WBOSON/unfolding_comparison.png and b/modules/quantum_lab/P2_WBOSON/unfolding_comparison.png differ diff --git a/modules/quantum_lab/P3_G2/cuda_accelerator.py b/modules/quantum_lab/P3_G2/cuda_accelerator.py deleted file mode 100644 index 77f834f..0000000 --- a/modules/quantum_lab/P3_G2/cuda_accelerator.py +++ /dev/null @@ -1,130 +0,0 @@ -""" -P3_G2/cuda_accelerator.py -Native C++/CUDA kernel for Front 3. -Computes the U(1) plaquette (Wilson action) and its block reduction (sum) -directly in the VRAM of the RTX 5060 Ti, bypassing the Python GIL. -""" -import numpy as np -try: - import cupy as cp -except ImportError: - cp = None - print("[CUDA] Warning: CuPy not detected. The CUDA accelerator will not run.") - -# ============================================================================== -# THE RAW C++ KERNEL (JIT-compiled at runtime by nvcc) -# ============================================================================== -# This kernel computes cos(P) at every site (x, y) of the lattice and then -# performs a block-level reduction using shared memory and __syncthreads() -# to guarantee a race-free sum. -# ============================================================================== - -cuda_source = r''' -extern "C" __global__ -void compute_plaquette_kernel(const float* theta, float* block_sums, int L) { - // Dynamic shared memory for the intra-block reduction - extern __shared__ float sdata[]; - - // Absolute lattice coordinates - int x = blockIdx.x * blockDim.x + threadIdx.x; - int y = blockIdx.y * blockDim.y + threadIdx.y; - - // Linear thread ID within the current block - int tid = threadIdx.y * blockDim.x + threadIdx.x; - - float val = 0.0f; - - // Stay inside the physical lattice bounds - if (x < L && y < L) { - // Map the 3D NumPy/CuPy tensor to contiguous 1D in C++: - // theta shape is (2, L, L) - // Channel 0: link in x (mu_0) - // Channel 1: link in y (mu_1) - - int offset_c0 = 0 * (L * L); - int offset_c1 = 1 * (L * L); - - int idx_t0_xy = offset_c0 + x * L + y; - int idx_t1_xy = offset_c1 + x * L + y; - - // Periodic boundary conditions (torus) - int xp1 = (x + 1) % L; - int yp1 = (y + 1) % L; - - int idx_t1_xp1_y = offset_c1 + xp1 * L + y; - int idx_t0_x_yp1 = offset_c0 + x * L + yp1; - - // P(x) = theta_0(x,y) + theta_1(x+1,y) - theta_0(x,y+1) - theta_1(x,y) - float t0 = theta[idx_t0_xy]; - float t1_shifted = theta[idx_t1_xp1_y]; - float t0_shifted = theta[idx_t0_x_yp1]; - float t1 = theta[idx_t1_xy]; - - float P = t0 + t1_shifted - t0_shifted - t1; - val = cosf(P); - } - - // Store this thread's value into the block's shared memory - sdata[tid] = val; - __syncthreads(); // BARRIER: wait until every thread has written its value - - // PARALLEL TREE REDUCTION - // Iteratively sum the halves of the block (assumes power-of-two block size) - int block_size = blockDim.x * blockDim.y; - for (unsigned int s = block_size / 2; s > 0; s >>= 1) { - if (tid < s) { - sdata[tid] += sdata[tid + s]; - } - __syncthreads(); // BARRIER: synchronize at every reduction step - } - - // The block's master thread (tid == 0) writes the partial sum to global memory - if (tid == 0) { - int block_id = blockIdx.y * gridDim.x + blockIdx.x; - block_sums[block_id] = sdata[0]; - } -} -''' - -# CuPy JIT compiler -if cp is not None: - plaquette_kernel = cp.RawKernel(cuda_source, 'compute_plaquette_kernel') - -def action_cuda(theta_np, beta, L): - """ - Python wrapper to invoke the C++ kernel for the Wilson action. - theta_np: ndarray (2, L, L) - """ - if cp is None: - raise RuntimeError("CuPy is not available.") - - # Copy host RAM to device VRAM (H2D) - theta_gpu = cp.asarray(theta_np, dtype=cp.float32) - - # Configure the CUDA thread topology - threads_x = 8 - threads_y = 8 - # Enough blocks to cover the full lattice L - blocks_x = (L + threads_x - 1) // threads_x - blocks_y = (L + threads_y - 1) // threads_y - - num_blocks = blocks_x * blocks_y - block_sums_gpu = cp.zeros(num_blocks, dtype=cp.float32) - - # 1 float (4 bytes) per thread in the block for shared memory - shared_mem_bytes = (threads_x * threads_y) * 4 - - # Launch the kernel - # signature: (grid, block, args, shared_mem) - plaquette_kernel((blocks_x, blocks_y), (threads_x, threads_y), - (theta_gpu, block_sums_gpu, L), - shared_mem=shared_mem_bytes) - - # Final reduction over the per-block partial sums (often a single block for small L) - total_sum = cp.sum(block_sums_gpu) - - # Action = -beta * sum(cos(P)) - action_val = -beta * total_sum - - # Bring the final scalar back to the CPU (D2H) - return action_val.get() diff --git a/modules/quantum_lab/P3_G2/main.py b/modules/quantum_lab/P3_G2/main.py deleted file mode 100644 index 56ade5c..0000000 --- a/modules/quantum_lab/P3_G2/main.py +++ /dev/null @@ -1,128 +0,0 @@ -""" -P3_G2/main.py -Phase 1 Orchestrator: U(1) Gauge Theory -> CNF -> Quimb -Calculates electromagnetic vacuum fluctuations (2D Lattice). -""" -import os -import numpy as np -import matplotlib.pyplot as plt - -from lattice_hmc import generate_u1_configs, plaquette -from cnf_flow import train_cnf -from tn_quimb import contract_lattice -from visualization_3d import plot_vacuum_3d - -def main(): - print("="*60) - print(" P3_G2: U(1) LATTICE | CNF | TENSOR NETWORKS ") - print("="*60) - - # Toy model parameters - L = 8 - beta = 1.0 # Inverse coupling (Temperature-like parameter) - - # ============================================================== - # 1. QUANTUM VACUUM GENERATION (Lattice QCD analog) - # ============================================================== - print("\n>>> PHASE A: Substrate Shadow Sampling (HMC)") - - # --- PRECISION CHECK: JAX vs RAW CUDA --- - # Numerical-agreement check only. For honest timing (warmed both sides, - # CUDA-event kernel timing, 1000-iter medians, size sweep, roofline - # crossover) run: python P3_G2/benchmark_plaquette.py - try: - from cuda_accelerator import action_cuda - from lattice_hmc import action as action_jax - import jax.numpy as jnp - - print("\n[CUDA] Verifying kernel against JAX reference (values only)...") - test_theta = np.random.uniform(-np.pi, np.pi, (2, L, L)).astype(np.float32) - - action_j = float(action_jax(jnp.array(test_theta), beta)) - action_c = float(action_cuda(test_theta, beta, L)) - - print(f" JAX Action : {action_j:.7f}") - print(f" CUDA Action : {action_c:.7f}") - diff = abs(action_j - action_c) - if diff < 1e-4: - print(f" [OK] Kernel matches JAX within fp32 tolerance (diff: {diff:.2e})") - else: - print(f" [WARNING] Divergence detected (diff: {diff:.2e})") - print(" Timing intentionally omitted here -> see benchmark_plaquette.py") - except Exception as e: - print(f"[CUDA] Precision check error: {e}") - - configs_hmc = generate_u1_configs(L=L, beta=beta, n_configs=200, n_steps=10, eps=0.1) - - # Extract primary physical observable: Plaquette Energy - import jax - plaqs = [np.mean(np.cos(jax.device_get(plaquette(c)))) for c in configs_hmc] - print(f"[Physics] Mean Plaquette Energy (HMC): {np.mean(plaqs):.4f} +/- {np.std(plaqs):.4f}") - - # ============================================================== - # 2. MACHINE LEARNING: CACHE LEARNING (Normalizing Flows) - # ============================================================== - print("\n>>> PHASE B: Training Generative Oracle (CNF)") - configs_ai = train_cnf(configs_hmc, epochs=50) - - plaqs_ai = [np.mean(np.cos(jax.device_get(plaquette(c)))) for c in configs_ai] - print(f"[Physics] Mean Plaquette Energy (AI): {np.mean(plaqs_ai):.4f} +/- {np.std(plaqs_ai):.4f}") - - # Visualization of Coherence Field assimilation - plt.figure(figsize=(9, 6)) - plt.hist(plaqs, bins=15, alpha=0.6, label="Quantum Mechanics (HMC)", density=True, color='blue') - plt.hist(plaqs_ai, bins=15, alpha=0.6, label="Artificial Intelligence (CNF)", density=True, color='orange') - plt.xlabel("Plaquette Energy $cos(P_{01})$") - plt.ylabel("Probability Density") - plt.title("U(1) Vacuum Generation: Physical Sampling vs AI Reconstruction") - plt.legend() - plt.grid(alpha=0.3) - - filename = "vacuum_shadow.png" - plt.savefig(filename) - print(f"[Main] Comparative shadow plot saved to '{filename}'.") - - print("\n[Main] Rendering 3D topology in graphics engine (PyVista)...") - try: - plot_vacuum_3d(configs_ai[-1], title="Vacuum Shadow - Generative Neural Network") - except Exception as e: - print(f"[Main] Error launching 3D visualization: {e}") - - # ============================================================== - # 3. EXACT CONTRACTION (Tensor Networks) - # ============================================================== - print("\n>>> PHASE C: Analytical Contraction of Tensor Network") - # Reduce L to 4 to ensure exact local contraction is viable within seconds - contract_lattice(L=4, beta=beta) - - # ============================================================== - # 4. SOVEREIGN TRACKING (Vacuum DevOps) - # ============================================================== - print("\n[Tracking] Logging metrics to central database...") - try: - import sys - # Add parent directory (QUANTUM_LAB) to path - sys.path.append(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) - from tracking import log_experiment - - log_experiment( - project="P3_G2", - parameters={"beta": beta, "lattice_size": f"{L}x{L}", "hmc_configs": 200, "cnf_epochs": 50}, - metrics={ - "plaquette_energy_hmc": float(np.mean(plaqs)), - "plaquette_energy_ai": float(np.mean(plaqs_ai)), - "logZ": 15.192149, - "free_energy": -0.949509 - }, - artifacts=["vacuum_shadow.png", "main.py"], - notes="Execution with native C++ CUDA acceleration and PyVista 3D rendering." - ) - except Exception as e: - print(f"[Tracking] Error connecting to SQLite DB: {e}") - - print("\n" + "="*60) - print(" PHASE 1 COMPLETED. SUBSTRATE CACHE HAS BEEN DECODED. ") - print("="*60) - -if __name__ == "__main__": - main() diff --git a/modules/quantum_lab/README.md b/modules/quantum_lab/README.md deleted file mode 100644 index cbf0d79..0000000 --- a/modules/quantum_lab/README.md +++ /dev/null @@ -1,186 +0,0 @@ -# QUANTUM-LAB - -Lattice field theory research facility. Three independent particle physics projects sharing a common infrastructure: CUDA-accelerated gauge computations, normalizing flow sampling, and a sovereign experiment tracking system. - -**Target hardware:** RTX 5060 Ti 16 GB (sm_120, Blackwell). No cloud. No notebooks-as-a-service. Bare metal. - ---- - -## results - -### P3 — U(1) gauge theory (CUDA reactor) - -Hand-written CUDA kernel (CuPy `RawKernel`, sm_120 Blackwell) for the U(1) -plaquette/Wilson action, benchmarked against a `@jax.jit` reference. Timing is -honest: both sides warmed up, JAX `block_until_ready()`, CUDA-event timing for -kernel-only GPU work, 1000-iteration medians. **JAX runs on CPU in this build** -(no GPU backend installed), so these are GPU-vs-CPU numbers — stated as such. - -| L | cells | JAX-CPU (ms) | CUDA kernel (ms) | CUDA e2e (ms) | kernel speedup | e2e speedup | max\|Δ\| | -|----:|------:|-------------:|-----------------:|--------------:|---------------:|------------:|--------:| -| 8 | 128 | 0.0272 | 0.00819 | 0.3767 | 3.3× | 0.07× | 0.0 | -| 16 | 512 | 0.0303 | 0.00816 | 0.4207 | 3.7× | 0.07× | 9.5e-07 | -| 32 | 2048 | 0.0367 | 0.00819 | 0.4054 | 4.5× | 0.09× | 4.8e-06 | -| 64 | 8192 | 0.1812 | 0.00822 | 0.4164 | 22.0× | 0.44× | 1.9e-05 | -| 128 | 32768 | 0.4183 | 0.00845 | 0.4470 | 49.5× | 0.94× | 1.9e-05 | -| 256 | 131072 | 0.6252 | 0.00870 | 0.4777 | 71.8× | 1.31× | 3.1e-05 | -| 512 | 524288 | 2.0204 | 0.01456 | 0.8433 | 138.8× | 2.40× | 4.1e-04 | - -**The honest story is the crossover, not a single hero number.** The CUDA kernel -itself is a stable 8–15 µs across all sizes; the kernel-only *speedup* rises from -~3× (small L) to roughly 139–157× at L=512 across runs (this table's run: 139×). -That top-end spread is JAX-CPU **baseline** noise (CPU wall-clock), not kernel -variance — the table above is one representative run. End-to-end (including H2D/D2H PCIe transfers) only -breaks even around **L≈128** — below that the workload is transfer/launch-bound -and the GPU loses. Numerical agreement with JAX is < 1e-3 (fp32 accumulation; -< 1e-4 up to L≈256, growing to ~4e-4 at L=512 as more terms are summed in fp32). -Reproduce: `python P3_G2/benchmark_plaquette.py` → `benchmark_plaquette_results.json`. - -HMC thermalization: 200 configs on 8×8 lattice, β=1.0, acceptance rate 0.98. -Plaquette energy: 0.4505 ± 0.0779. - -### P1 — LHCb anomaly detection - -GNN Autoencoder on CERN Open Data (b→sℓℓ decays). ROC-AUC: **0.765** on 4-feature latent space. Quantum Boltzmann Machine background model via PennyLane. - -### P2 — W-boson quantum unfolding - -Sequential QUBO unfolding of Jacobian peak distributions. Sliding window → Cirq/PennyLane solver. Benchmarked against classical SVD baseline. - ---- - -## architecture - -``` -QUANTUM_LAB/ -├── P1_LHCB/ -│ ├── main.py ← orchestrator (ingestion → GNN → QBM) -│ ├── ingestion.py ← uproot + graph construction -│ ├── gnn_autoencoder.py ← PyTorch Geometric autoencoder -│ └── qbm_pennylane.py ← Quantum Boltzmann Machine -│ -├── P2_WBOSON/ -│ ├── data_generator.py ← Jacobian peak toy data -│ ├── quantum_unfolder.py ← QUBO sliding window solver -│ └── benchmark.py ← SVD vs quantum comparison -│ -├── P3_G2/ -│ ├── main.py ← orchestrator (HMC → CNF → TN → CUDA) -│ ├── lattice_hmc.py ← JAX Hybrid Monte Carlo for U(1) -│ ├── cnf_flow.py ← Continuous Normalizing Flow (PyTorch) -│ ├── tn_quimb.py ← Tensor network contraction (quimb) -│ ├── cuda_accelerator.py ← CuPy RawKernel C++ plaquette action -│ └── visualization_3d.py ← PyVista/VTK bare-metal 3D rendering -│ -├── tracking.py ← SQLite sovereign experiment tracker -├── query.py ← CLI experiment query interface -├── plot_metrics.py ← Offline metric visualization -└── experiments.db ← Local experiment database -``` - ---- - -## stack - -| layer | technology | -|-------|-----------| -| gauge sampling | JAX (HMC with leapfrog integrator) | -| generative model | PyTorch + torchdiffeq (Neural ODE / CNF) | -| tensor networks | quimb + cotengra (PEPS, CTMRG, BMPS) | -| GPU acceleration | CuPy RawKernel (native C++ CUDA) | -| particle data | uproot (ROOT → NumPy) | -| quantum circuits | PennyLane + Cirq | -| 3D visualization | PyVista / VTK (bare-metal OpenGL) | -| experiment tracking | SQLite (no MLflow, no cloud) | - ---- - -## the CUDA kernel - -The plaquette action kernel maps each lattice site to a CUDA thread. The gauge -field U(1) action is computed as: - -``` -S = -β Σ Re(U_μ(n) · U_ν(n+μ) · U_μ(n+ν)* · U_ν(n)*) -``` - -Each thread computes one plaquette and accumulates into shared memory; a -block-level reduction (`__syncthreads()`, tree reduction over an 8×8 tile) -produces one partial sum per block, summed on-device with `cp.sum`. Keeping the -reduction in shared memory minimizes global-memory traffic, and the 8×8 block -maps cleanly to the 2D lattice for coalesced loads. - -```cpp -extern "C" __global__ void compute_plaquette_kernel( - const float* theta, // gauge angles, shape (2, L, L) flattened - float* block_sums, // one partial sum per block - int L) // lattice size -``` - -**Why kernel-only ≠ end-to-end.** The kernel itself is tiny (8–15 µs across all -sizes tested). At small L the runtime is dominated by the fixed cost of the H2D -copy of `theta` plus the D2H copy of the scalar result — that is why the -end-to-end column only beats JAX-CPU past L≈128. For a single action evaluation -this is the launch-overhead-bound regime; the win comes when the field already -lives on device (as it does inside an HMC trajectory, where the same buffer is -reused across leapfrog steps and no per-step transfer is paid). - ---- - -## sovereign tracking - -Every experiment run is logged to `experiments.db` (SQLite) with: -- Parameters (JSON): lattice size, β, epochs, bond dimension -- Metrics (JSON): plaquette energy, acceptance rate, loss curves -- Artifacts: paths to generated plots and configs -- Timestamp and run notes - -No MLflow. No localhost servers. No web dashboards. Query with `python query.py`. - ---- - -## run - -```bash -# full pipeline (HMC → CNF → TN → CUDA benchmark) -python P3_G2/main.py - -# LHCb anomaly detection -python P1_LHCB/main.py - -# W-boson unfolding benchmark -python P2_WBOSON/benchmark.py - -# query experiment history -python query.py --last 10 -``` - ---- - -## requirements - -``` -jax[cuda12] -torch -torchdiffeq -torch_geometric -quimb -cotengra -cupy-cuda12x -pennylane -cirq -uproot -pyvista -``` - ---- - -## related - -- [CryptoTN-GPU](https://github.com/QuantumDrizzy/CryptoTN-GPU) — GPU tensor networks for quantum biology -- [KHAOS](https://github.com/QuantumDrizzy/KHAOS) — BCI kernel with CUDA DSP -- [quantum-geo-metrology](https://github.com/QuantumDrizzy/quantum-geo-metrology) — geophysical quantum computing - ---- - -*Antonio Rodríguez (QuantumDrizzy) · research software engineer* diff --git a/modules/quantum_lab/reports/P3_G2_logZ.png b/modules/quantum_lab/reports/P3_G2_logZ.png index f50499d..27adc63 100644 Binary files a/modules/quantum_lab/reports/P3_G2_logZ.png and b/modules/quantum_lab/reports/P3_G2_logZ.png differ diff --git a/nexus/substrate_material_bridge.py b/nexus/substrate_material_bridge.py deleted file mode 100644 index 64bae7e..0000000 --- a/nexus/substrate_material_bridge.py +++ /dev/null @@ -1,238 +0,0 @@ -""" -nexus/substrate_material_bridge.py — SUBSTRATE → HELIOS material_states bridge -================================================================================ -Runs CryptoTN-GPU fast benchmarks (ErCry4a + FMO) and writes quantum material -outputs into HELIOS's material_states SQLite table. - -Mapping: - quantum_efficiency ← ErCry4a Φ_S(B=0) — singlet yield at zero field - electron_mobility ← ErCry4a ΔΦ_S — compass sensitivity (earth field) - lattice_stability ← 1 - FMO RMSE — FMO fidelity vs TENSO reference - temp_celsius ← last HELIOS telemetry ambient (default 25.0) - -Usage (standalone): - python nexus/substrate_material_bridge.py - -Usage (from pipeline.py): - from nexus.substrate_material_bridge import run_material_bridge - run_material_bridge() -""" - -from __future__ import annotations - -import json -import os -import sys -import sqlite3 -import time -import logging -from pathlib import Path - -import zmq - -logging.basicConfig( - level=logging.INFO, - format="%(asctime)s [SUBSTRATE.MATERIAL] %(message)s", -) -log = logging.getLogger(__name__) - -# ── paths ────────────────────────────────────────────────────────────────────── -_HERE = Path(__file__).parent # SUBSTRATE/nexus/ -_ROOT = _HERE.parent # SUBSTRATE/ -_CRYPTOTN = _ROOT / "modules" / "cryptotn_gpu" -_HELIOS_DB = Path(os.environ.get( - "HELIOS_DB", - r"C:\Users\Drizzy\Desktop\HELIOS\data\energy_bus.sqlite" -)) - -# add cryptotn_gpu to path so its imports resolve -if str(_CRYPTOTN) not in sys.path: - sys.path.insert(0, str(_CRYPTOTN)) - -# ── benchmark imports (deferred — only fail if module missing) ───────────────── -def _import_benchmarks(): - """Import benchmark runners. Returns (run_ercry4a, run_fmo) or raises.""" - benchmarks = _CRYPTOTN / "benchmarks" - if str(benchmarks) not in sys.path: - sys.path.insert(0, str(benchmarks)) - from bench_ercry4a import run_ercry4a # noqa: PLC0415 - from bench_fmo import run_fmo # noqa: PLC0415 - return run_ercry4a, run_fmo - - -# ── HELIOS DB helpers ────────────────────────────────────────────────────────── -def _get_ambient_temp(conn: sqlite3.Connection) -> float: - """Read latest ambient proxy from power_telemetry (placeholder: 25.0°C).""" - row = conn.execute( - "SELECT voltage FROM power_telemetry ORDER BY id DESC LIMIT 1" - ).fetchone() - # No real thermistor yet — return 25.0 as lab ambient - return 25.0 - - -def _write_material_state( - conn: sqlite3.Connection, - quantum_efficiency: float, - electron_mobility: float, - lattice_stability: float, - temp_celsius: float, -) -> int: - """Insert one row into material_states. Returns new row id.""" - cursor = conn.execute( - """ - INSERT INTO material_states - (quantum_efficiency, electron_mobility, lattice_stability, temp_celsius) - VALUES (?, ?, ?, ?) - """, - ( - round(quantum_efficiency, 6), - round(electron_mobility, 6), - round(lattice_stability, 6), - round(temp_celsius, 2), - ), - ) - conn.commit() - return cursor.lastrowid - - -# ── main bridge logic ────────────────────────────────────────────────────────── -def run_material_bridge( - n_nuc: int = 10, - fast: bool = True, -) -> dict | None: - """ - Run CryptoTN-GPU benchmarks and write results to HELIOS material_states. - - Parameters - ---------- - n_nuc : int - Nuclear spins for ErCry4a (10 = fast ExactSolver, ~1s). - fast : bool - If True: ErCry4a n_steps=100, FMO n_steps=100. ~2-3s total. - - Returns - ------- - dict with keys: quantum_efficiency, electron_mobility, lattice_stability, - temp_celsius, row_id, wall_time_s - """ - if not _HELIOS_DB.exists(): - log.error(f"HELIOS DB not found: {_HELIOS_DB}") - return None - - log.info(f"HELIOS DB: {_HELIOS_DB}") - log.info(f"CryptoTN-GPU path: {_CRYPTOTN}") - - try: - run_ercry4a, run_fmo = _import_benchmarks() - except ImportError as exc: - log.error(f"CryptoTN-GPU import failed: {exc}") - return None - - t0 = time.perf_counter() - - # ── 1. ErCry4a — singlet yield + compass sensitivity ────────────────────── - n_steps_cry = 100 if fast else 300 - log.info(f"Running ErCry4a ({n_nuc} nuclei, n_steps={n_steps_cry})...") - try: - ercry = run_ercry4a( - n_nuc=n_nuc, - B_field_values=[0.0, 0.05], # only need B=0 and earth-field - t_max_us=5.0, - n_steps=n_steps_cry, - chi=32, - use_mps=False, - ) - phi_s_zero = float(ercry["phi_s"][0]) # Φ_S at B=0 - delta_phi_s = float(ercry["delta_phi_s_earth"]) - except Exception as exc: - log.error(f"ErCry4a failed: {exc}") - phi_s_zero = 0.5 # neutral fallback (max entangled) - delta_phi_s = 0.0 - - # ── 2. FMO 77K — lattice stability via TENSO fidelity ──────────────────── - n_steps_fmo = 100 if fast else 200 - log.info(f"Running FMO 77K (n_steps={n_steps_fmo})...") - try: - fmo = run_fmo(T_K=77.0, t_max_fs=500.0, n_steps=n_steps_fmo) - rmse_tenso = float(fmo["rmse_vs_tenso"]) - except Exception as exc: - log.error(f"FMO failed: {exc}") - rmse_tenso = 0.0 - - # ── 3. Map to material_states schema ────────────────────────────────────── - quantum_efficiency = phi_s_zero # range ~[0.45, 0.55] - electron_mobility = delta_phi_s # range ~[0.001, 0.010] - lattice_stability = max(0.0, 1.0 - rmse_tenso * 10.0) # invert + scale - - # ── 4. Write to HELIOS DB ───────────────────────────────────────────────── - conn = sqlite3.connect(str(_HELIOS_DB)) - try: - temp_c = _get_ambient_temp(conn) - row_id = _write_material_state( - conn, - quantum_efficiency=quantum_efficiency, - electron_mobility=electron_mobility, - lattice_stability=lattice_stability, - temp_celsius=temp_c, - ) - finally: - conn.close() - - wall = time.perf_counter() - t0 - - result = { - "quantum_efficiency": round(quantum_efficiency, 6), - "electron_mobility": round(electron_mobility, 6), - "lattice_stability": round(lattice_stability, 6), - "temp_celsius": temp_c, - "row_id": row_id, - "wall_time_s": round(wall, 2), - "source": { - "phi_s_zero": round(phi_s_zero, 6), - "delta_phi_s": round(delta_phi_s, 6), - "rmse_vs_tenso": round(rmse_tenso, 6), - }, - } - - log.info( - f"material_states row {row_id} written — " - f"Φ_S={phi_s_zero:.5f} ΔΦ_S={delta_phi_s:.5f} " - f"RMSE={rmse_tenso:.5f} wall={wall:.2f}s" - ) - - # ── 5. Publish to nexus.quantum (port 5562) ─────────────────────────────── - try: - ctx = zmq.Context() - pub = ctx.socket(zmq.PUB) - pub.bind("tcp://*:5562") - pub.setsockopt(zmq.LINGER, 1000) # wait up to 1s on close - time.sleep(0.5) # warmup — give LOGOS time to detect new publisher - frame = json.dumps({ - "ts": time.time(), - "event": "substrate_result", - "payload": result, - }) - pub.send_string(f"nexus.quantum {frame}") # blocking send - time.sleep(0.2) # let subscribers drain before close - pub.close() - ctx.term() - log.info("[NEXUS.QUANTUM] substrate_result published on :5562") - except Exception as exc: - log.warning(f"[NEXUS.QUANTUM] publish failed (non-fatal): {exc}") - - return result - - -# ── CLI ──────────────────────────────────────────────────────────────────────── -if __name__ == "__main__": - import argparse, json - parser = argparse.ArgumentParser(description="SUBSTRATE → HELIOS material bridge") - parser.add_argument("--n-nuc", type=int, default=10, help="nuclear spins (default 10)") - parser.add_argument("--full", action="store_true", help="full benchmark (slower)") - args = parser.parse_args() - - result = run_material_bridge(n_nuc=args.n_nuc, fast=not args.full) - if result: - print(json.dumps(result, indent=2)) - else: - sys.exit(1) diff --git a/pipeline.py b/pipeline.py deleted file mode 100644 index 2696110..0000000 --- a/pipeline.py +++ /dev/null @@ -1,75 +0,0 @@ -""" -SUBSTRATE — Unified Research Pipeline -====================================== -Multi-scale electromagnetic field analysis: - quantum substrate → bio sensing → geomagnetic → heliospheric → cosmological - -Usage: - python pipeline.py --module all - python pipeline.py --module cycle_project - python pipeline.py --module cycle_project,cosmological -""" - -import argparse -import subprocess -import sys -from pathlib import Path - -ROOT = Path(__file__).parent -MODULES = ROOT / "modules" - -MODULE_RUNNERS = { - "cycle_project": MODULES / "cycle_project" / "src" / "pipeline.py", - # Add more as they are ready: - # "quantum_lab": MODULES / "quantum_lab" / "run.py", - # "magnon": MODULES / "magnon" / "run.py", -} - -# Modules with custom Python entry points (not subprocess runners) -INLINE_RUNNERS = { - "cryptotn_gpu": "nexus.substrate_material_bridge:run_material_bridge", -} - -def run_module(name: str): - # Check inline Python runners first - inline = INLINE_RUNNERS.get(name) - if inline: - print(f"\n{'='*60}\n SUBSTRATE: running {name} (inline)\n{'='*60}") - mod_path, func_name = inline.rsplit(":", 1) - # ensure SUBSTRATE root is in path - if str(ROOT) not in sys.path: - sys.path.insert(0, str(ROOT)) - try: - import importlib - mod = importlib.import_module(mod_path) - fn = getattr(mod, func_name) - result = fn() - if result: - import json - print(json.dumps(result, indent=2)) - except Exception as exc: - print(f"[SUBSTRATE] {name} failed: {exc}") - return - - runner = MODULE_RUNNERS.get(name) - if runner is None: - print(f"[SUBSTRATE] Module '{name}' not yet wired into pipeline.") - return - if not runner.exists(): - print(f"[SUBSTRATE] Runner not found: {runner}") - return - print(f"\n{'='*60}\n SUBSTRATE: running {name}\n{'='*60}") - result = subprocess.run([sys.executable, str(runner)], cwd=runner.parent.parent.parent) - if result.returncode != 0: - print(f"[SUBSTRATE] {name} exited with code {result.returncode}") - -if __name__ == "__main__": - parser = argparse.ArgumentParser(description="SUBSTRATE pipeline") - parser.add_argument("--module", default="all", help="Module(s) to run, comma-separated, or 'all'") - args = parser.parse_args() - - all_modules = list(MODULE_RUNNERS.keys()) + list(INLINE_RUNNERS.keys()) - targets = all_modules if args.module == "all" else args.module.split(",") - for mod in targets: - run_module(mod.strip()) - print("\n[SUBSTRATE] Done.") diff --git a/substrate.toml b/substrate.toml deleted file mode 100644 index 1c415cb..0000000 --- a/substrate.toml +++ /dev/null @@ -1,41 +0,0 @@ -[substrate] -name = "SUBSTRATE" -version = "0.2.0" - -[layers] -quantum = { enabled = true, weight = 1.5, params = { chi = 256, backend = "cupy" } } -geomagnetic = { enabled = true, weight = 2.0 } -magnon = { enabled = true, weight = 1.2 } -quantum_lab = { enabled = true, weight = 1.0, params = { chi = 512 } } -solar = { enabled = true, weight = 1.3 } -cosmological = { enabled = true, weight = 0.8 } -eeg = { enabled = true, weight = 1.8, params = { mode = "auto" } } -lunar = { enabled = true, weight = 1.4 } -radio = { enabled = true, weight = 1.6 } -seismic = { enabled = true, weight = 1.1 } - -# ── New modules (integrated 2026-05-28) ─────────────────────────── -[layers.simulations] -enabled = true -weight = 1.0 -description = "Multi-domain physics simulations: cosmology, quantum, astrophysics" -submodules = ["cosmology", "quantum", "astrophysics", "bci_bridge", "sonification"] -data_dir = "modules/simulations/data" - -[layers.simulations.bci_bridge] -enabled = true -weight = 1.8 -description = "KHAOS EEG → physics parameter bridge (mind_state.json IPC)" -mind_state_ttl_secs = 15 -connects_to = ["khaos", "helios"] - -[layers.pattern_analysis] -enabled = true -weight = 1.0 -description = "Cross-modal spectral analysis: Rust Welch FFT + Python PNG/JSON analytics" -rust_engine = "modules/pattern_analysis/rust_engine" -connects_to = ["simulations"] - -[output] -db_path = "data/substrate.db" -report_path = "data/processed/substrate_report.json" diff --git a/tools/exp0_logger.py b/tools/exp0_logger.py deleted file mode 100644 index 0bfa04b..0000000 --- a/tools/exp0_logger.py +++ /dev/null @@ -1,186 +0,0 @@ -#!/usr/bin/env python3 -"""Experiment 0 data logger -- persists (timestamp, dst_nT, singlet_yield) to disk. - -The SUBSTRATE GUI keeps only 240 frames in RAM. This daemon runs alongside the -motor, calls magnon_layer.run() on a fixed cadence, and appends each result to -a JSONL file that survives restarts. - -Usage: - python3 tools/exp0_logger.py # 60s cadence - python3 tools/exp0_logger.py --interval 30 - python3 tools/exp0_logger.py --once # single sample and exit - python3 tools/exp0_logger.py --export-csv # convert log to CSV for analysis - -Output: - ~/.cache/substrate/experiment_0/log.jsonl - -Each line: - {"t": 1778782109, "utc": "2026-05-14T18:08:29Z", - "dst_nT": -42.0, "dst_source": "noaa_dst_cache", - "singlet_yield": 0.7931, "fidelity": 0.9812, ...} -""" -from __future__ import annotations - -import argparse -import json -import multiprocessing -import os as _os -import signal -import sys -import time -from pathlib import Path - -_ROOT = Path(__file__).resolve().parent.parent -_LOG_DIR = Path.home() / ".cache" / "substrate" / "experiment_0" -_LOG_FILE = _LOG_DIR / "log.jsonl" -_MARKER = Path.home() / ".cache" / "substrate" / "geomagnetic" / "experiment_0_marker.json" - - -def _ensure_dirs(): - _LOG_DIR.mkdir(parents=True, exist_ok=True) - - -def _worker(queue: multiprocessing.Queue, root: str): - """Runs in a fresh subprocess — all C extensions and SDR buffers die with it.""" - import sys, os, time - sys.path.insert(0, str(root + "/engine")) - sys.path.insert(0, str(root + "/modules/magnon")) - - t0 = time.time() - utc = time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime(t0)) - try: - import magnon_layer # type: ignore - data = magnon_layer.run().get("data", {}) - queue.put({ - "t": int(t0), - "utc": utc, - "dst_nT": data.get("dst_nT", 0.0), - "dst_source": data.get("dst_source", "unknown"), - "singlet_yield": data.get("singlet_yield"), - "fidelity": data.get("fidelity"), - "noise_source": data.get("noise_source"), - "noise_power_ratio": data.get("noise_power_ratio"), - "b_earth_nT": data.get("b_earth_nT"), - "solve_ms": data.get("solve_ms"), - "error": data.get("error"), - }) - except Exception as exc: - queue.put({"t": int(t0), "utc": utc, "dst_nT": None, - "dst_source": "logger_error", "error": str(exc)}) - - -def _sample(): - """Spawn a subprocess for each sample — memory fully reclaimed on exit.""" - q = multiprocessing.Queue() - p = multiprocessing.Process(target=_worker, args=(q, str(_ROOT)), daemon=True) - p.start() - p.join(timeout=120) # 2 min hard timeout - if p.exitcode is None: - p.terminate() - t0 = int(time.time()) - return {"t": t0, "utc": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime(t0)), - "dst_nT": None, "dst_source": "logger_error", "error": "worker timeout"} - if not q.empty(): - return q.get_nowait() - t0 = int(time.time()) - return {"t": t0, "utc": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime(t0)), - "dst_nT": None, "dst_source": "logger_error", "error": "worker returned nothing"} - - -def _append(record): - _ensure_dirs() - with _LOG_FILE.open("a", encoding="utf-8") as f: - f.write(json.dumps(record) + "\n") - - -def _print_record(record): - dst = record.get("dst_nT") - ys = record.get("singlet_yield") - err = record.get("error") - if err: - print(f"[{record['utc']}] ERROR: {err}") - else: - print(f"[{record['utc']}] dst={dst:+.1f} nT [{record.get('dst_source','?')}]" - f" Y_s={ys:.8f} noise={record.get('noise_source','?')}" - f" {record.get('solve_ms',0):.1f}ms") - - -def run_daemon(interval): - _ensure_dirs() - print(f"Experiment 0 logger | interval={interval}s | {_LOG_FILE}") - print("Ctrl+C to stop.\n") - - def _stop(sig, frame): - print("\nLogger stopped.") - try: - m = json.loads(_MARKER.read_text()) - m["end_utc"] = time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()) - _MARKER.write_text(json.dumps(m, indent=2)) - except Exception: - pass - sys.exit(0) - - signal.signal(signal.SIGINT, _stop) - signal.signal(signal.SIGTERM, _stop) - - n = 0 - while True: - rec = _sample() - _append(rec) - _print_record(rec) - n += 1 - if n % 60 == 0: - lines = sum(1 for _ in _LOG_FILE.open()) - print(f" [{n} samples, {lines} total lines on disk]") - time.sleep(interval) - - -def run_once(): - rec = _sample() - print(json.dumps(rec, indent=2)) - _append(rec) - print(f"\nAppended to {_LOG_FILE}") - - -def export_csv(out_path=None): - import csv - if not _LOG_FILE.exists(): - print(f"No log file at {_LOG_FILE}") - return - out = Path(out_path) if out_path else _LOG_DIR / "exp0_data.csv" - fields = ["t", "utc", "dst_nT", "dst_source", "singlet_yield", - "fidelity", "noise_power_ratio", "noise_source", "solve_ms"] - total = kept = 0 - with _LOG_FILE.open() as fin, out.open("w", newline="") as fout: - w = csv.DictWriter(fout, fieldnames=fields, extrasaction="ignore") - w.writeheader() - for line in fin: - try: - row = json.loads(line) - total += 1 - if row.get("dst_source") == "dst_unavailable": - continue - if row.get("singlet_yield") is None: - continue - w.writerow({k: row.get(k, "") for k in fields}) - kept += 1 - except json.JSONDecodeError: - continue - print(f"Exported {kept}/{total} rows -> {out} ({total-kept} dropped)") - return out - - -if __name__ == "__main__": - p = argparse.ArgumentParser() - p.add_argument("--interval", type=int, default=60) - p.add_argument("--once", action="store_true") - p.add_argument("--export-csv", action="store_true") - p.add_argument("--out", default=None) - args = p.parse_args() - - if args.export_csv: - export_csv(args.out) - elif args.once: - run_once() - else: - run_daemon(args.interval) diff --git a/workbench/app.js b/workbench/app.js deleted file mode 100644 index 7e95867..0000000 --- a/workbench/app.js +++ /dev/null @@ -1,340 +0,0 @@ -/* app.js */ -const container = document.getElementById('canvas-3d-container'); -const scene = new THREE.Scene(); -const camera = new THREE.PerspectiveCamera(45, container.clientWidth / container.clientHeight, 0.1, 1000); -const renderer = new THREE.WebGLRenderer({ antialias: true, alpha: true }); -renderer.setSize(container.clientWidth, container.clientHeight); -container.appendChild(renderer.domElement); - -const geometry = new THREE.PlaneBufferGeometry(100, 100, 50, 50); -const material = new THREE.MeshPhongMaterial({ color: 0x8a8a8a, wireframe: true, transparent: true, opacity: 0.3, side: THREE.DoubleSide }); -const mesh = new THREE.Mesh(geometry, material); -mesh.rotation.x = -Math.PI / 2.5; -scene.add(mesh); - -const light = new THREE.DirectionalLight(0xffffff, 1); -light.position.set(0, 50, 50); -scene.add(light); -scene.add(new THREE.AmbientLight(0xffffff, 0.4)); -camera.position.set(0, 40, 80); -const controls = new THREE.OrbitControls(camera, renderer.domElement); -controls.enableDamping = true; -controls.dampingFactor = 0.05; -controls.screenSpacePanning = false; -controls.minDistance = 20; -controls.maxDistance = 200; -controls.maxPolarAngle = Math.PI / 2; - -let time = 0; -let activeViz = mesh; - -// --- CONFIGURACIÓN DE ÓRBITAS NEURO-CUÁNTICAS (BCI) --- -let bciCalm = 0.5; -let bciFocus = 0.5; -let orbitTime = 0; - -const maxTrailPoints = 150; -const bciOrbitGeom1 = new THREE.BufferGeometry(); -const bciOrbitGeom2 = new THREE.BufferGeometry(); -const trailPositions1 = new Float32Array(maxTrailPoints * 3); -const trailPositions2 = new Float32Array(maxTrailPoints * 3); - -bciOrbitGeom1.setAttribute('position', new THREE.BufferAttribute(trailPositions1, 3)); -bciOrbitGeom2.setAttribute('position', new THREE.BufferAttribute(trailPositions2, 3)); - -const bciOrbitMat1 = new THREE.LineBasicMaterial({ color: 0x00ffcc, linewidth: 2, transparent: true, opacity: 0.8 }); -const bciOrbitMat2 = new THREE.LineBasicMaterial({ color: 0xf72585, linewidth: 2, transparent: true, opacity: 0.8 }); - -const bciOrbit1 = new THREE.Line(bciOrbitGeom1, bciOrbitMat1); -const bciOrbit2 = new THREE.Line(bciOrbitGeom2, bciOrbitMat2); -scene.add(bciOrbit1); -scene.add(bciOrbit2); - -// Esferas para los cuerpos celestes/entrelazados BCI -const sphereGeom = new THREE.SphereGeometry(1.5, 16, 16); -const sphereMat1 = new THREE.MeshPhongMaterial({ color: 0x00ffcc, emissive: 0x002211, shininess: 100 }); -const sphereMat2 = new THREE.MeshPhongMaterial({ color: 0xf72585, emissive: 0x220011, shininess: 100 }); -const bciSphere1 = new THREE.Mesh(sphereGeom, sphereMat1); -const bciSphere2 = new THREE.Mesh(sphereGeom, sphereMat2); -scene.add(bciSphere1); -scene.add(bciSphere2); - -function animate() { - time += 0.01; - orbitTime += 0.015; - controls.update(); // Actualizar interacción de usuario - - // --- CÁLCULO DE ÓRBITAS ACOPLADAS AL BCI --- - const C = bciCalm; // 0.0 (estrés) a 1.0 (calma) - - // 1. Keplerian Lemniscate (Figura de 8, Calma Total) - const A = 24.0; - const den = 1.0 + Math.pow(Math.sin(orbitTime), 2); - const x8 = (A * Math.cos(orbitTime)) / den; - const y8 = (A * Math.sin(orbitTime) * Math.cos(orbitTime)) / den; - const z8 = 3.0 * Math.sin(2.0 * orbitTime); - - // 2. Inspiral / Caos post-Newtoniano (Estrés Cognitivo) - const tPeriod = 8.0; - const cycleTime = (orbitTime % tPeriod) / tPeriod; - const R = A * (0.05 + 0.95 * Math.pow(1.0 - cycleTime, 1.5)); - const freqSpiraling = 6.0 + 12.0 * cycleTime; - const xSpir = R * Math.cos(orbitTime * freqSpiraling); - const ySpir = R * Math.sin(orbitTime * freqSpiraling); - const zSpir = R * 0.25 * Math.sin(orbitTime * freqSpiraling); - - // Mezclar posiciones de órbita según calm_index - const p1x = C * x8 + (1 - C) * xSpir; - const p1y = C * y8 + (1 - C) * ySpir; - const p1z = C * z8 + (1 - C) * zSpir; - - const p2x = -p1x; - const p2y = -p1y; - const p2z = -p1z; - - // Actualizar esferas (mapeamos z a Y para que floten arriba de la cuadrícula) - bciSphere1.position.set(p1x, p1z, p1y); - bciSphere2.position.set(p2x, p2z, p2y); - - // Actualizar trails de órbita - const arr1 = bciOrbitGeom1.attributes.position.array; - const arr2 = bciOrbitGeom2.attributes.position.array; - - for (let i = maxTrailPoints - 1; i > 0; i--) { - arr1[i * 3] = arr1[(i - 1) * 3]; - arr1[i * 3 + 1] = arr1[(i - 1) * 3 + 1]; - arr1[i * 3 + 2] = arr1[(i - 1) * 3 + 2]; - - arr2[i * 3] = arr2[(i - 1) * 3]; - arr2[i * 3 + 1] = arr2[(i - 1) * 3 + 1]; - arr2[i * 3 + 2] = arr2[(i - 1) * 3 + 2]; - } - - arr1[0] = p1x; - arr1[1] = p1z; - arr1[2] = p1y; - - arr2[0] = p2x; - arr2[1] = p2z; - arr2[2] = p2y; - - bciOrbitGeom1.attributes.position.needsUpdate = true; - bciOrbitGeom2.attributes.position.needsUpdate = true; - - // --- DEFORMACIÓN DE LA MALLA RELATIVISTA DEL ESPACIOTIEMPO --- - if (activeViz && activeViz.geometry && activeViz.geometry.attributes.position) { - const positions = activeViz.geometry.attributes.position.array; - - if (activeViz === mesh) { - // Modulamos la deformación del plano: calma alta = ondas suaves; estrés alto = espigas caóticas - let amp = 1.5 + 8.5 * (1.0 - bciCalm); - let freq = 0.08 + 0.15 * bciFocus; - - for (let i = 0; i < positions.length; i += 3) { - const x = positions[i]; - const y = positions[i + 1]; - let z = Math.sin(x * freq + time) * Math.cos(y * freq + time) * amp; - - // Distancia a las esferas en el plano horizontal (x, y) - const d1 = Math.sqrt(Math.pow(x - p1x, 2) + Math.pow(y - p1y, 2)); - const d2 = Math.sqrt(Math.pow(x - p2x, 2) + Math.pow(y - p2y, 2)); - // Singularidades gravitatorias locales (pozos de potencial) - const gravityWell = -16.0 / (d1 + 1.8) - 16.0 / (d2 + 1.8); - - positions[i + 2] = z + gravityWell; - } - } else { - // Comportamiento por defecto para otros modos - for (let i = 0; i < positions.length; i += 3) { - const x = positions[i]; const y = positions[i + 1]; - let z = (activeViz.type === 'Points') ? Math.sin(x/5 + time)*2 : Math.sin(x/10+time)*Math.cos(y/10+time)*5 + Math.sin(x/5-time*0.5)*2; - positions[i+2] = z; - } - } - activeViz.geometry.attributes.position.needsUpdate = true; - } - renderer.render(scene, camera); - requestAnimationFrame(animate); -} -animate(); - -// --- INTERACTIVIDAD DE LABORATORIO --- - -// Listener para cambios de Tolerancia -document.querySelector('input[value="1.0e-8"]').addEventListener('change', function(e) { - console.log("Nueva tolerancia:", e.target.value); - if (window.eel) eel.update_solver_param('tolerance', e.target.value); -}); - -// Listener para cambio de Algoritmo -document.querySelector('select').addEventListener('change', function(e) { - console.log("Nuevo algoritmo:", e.target.value); - if (window.eel) eel.update_solver_param('algorithm', e.target.value); -}); - -// Botones de Play/Pause -document.getElementById('playBtn').addEventListener('click', () => { - document.getElementById('playBtn').style.color = 'var(--accent-green)'; - if (window.eel) eel.set_engine_state('RUNNING'); -}); - -document.getElementById('stopBtn').addEventListener('click', () => { - document.getElementById('playBtn').style.color = ''; - if (window.eel) eel.set_engine_state('STOPPED'); -}); - -let charts = []; - -function clearCharts() { - charts.forEach(c => c.destroy()); - charts = []; -} - -function initCharts(config) { - clearCharts(); - const ctx1 = document.getElementById('chart1').getContext('2d'); - const ctx2 = document.getElementById('chart2').getContext('2d'); - - charts.push(new Chart(ctx1, { - type: config.chart1Type || 'line', - data: { labels: Array.from({length: 30}, (_, i) => i), datasets: [{ label: config.chart1Label, borderColor: config.colorHex, data: Array.from({length: 30}, () => Math.random()), tension: 0.4, pointRadius: 0 }] }, - options: { responsive: true, maintainAspectRatio: false, scales: { y: { min: 0, max: 1 } }, plugins: { legend: { display: false } } } - })); - - charts.push(new Chart(ctx2, { - type: config.chart2Type || 'line', - data: { labels: Array.from({length: 10}, (_, i) => i), datasets: [{ label: config.chart2Label, backgroundColor: config.colorHex + '44', borderColor: config.colorHex, data: Array.from({length: 10}, () => Math.random()) }] }, - options: { responsive: true, maintainAspectRatio: false, plugins: { legend: { display: false } } } - })); -} - -// Configuración de Instrumentos -const instruments = { - "Quantum Layer (GPU)": { mode: "MESH", color: 0x0071e3, colorHex: "#0071e3", chart1Label: "Q-Strain", chart2Label: "Phase Shift" }, - "Geomagnetic Layer": { mode: "VECTORS", color: 0x27ae60, colorHex: "#27ae60", chart1Label: "Mag Vector", chart2Label: "Radial Var", chart2Type: "polarArea" }, - "Magnon Layer (Bio)": { mode: "SPIKES", color: 0x8e44ad, colorHex: "#8e44ad", chart1Label: "Bio-Potential", chart2Label: "Signal Frequency", chart2Type: "bar" }, - "Quantum Lab (TN)": { mode: "LATTICE", color: 0xd35400, colorHex: "#d35400", chart1Label: "Fidelity", chart2Label: "Bond Entropy", chart2Type: "bar" }, - "Solar Layer (Cycle)": { mode: "SPHERE", color: 0xf1c40f, colorHex: "#f1c40f", chart1Label: "Solar Flux", chart2Label: "Thermal Var" }, - "Cosmological Layer": { mode: "POINTS", color: 0x34495e, colorHex: "#34495e", chart1Label: "Expansion", chart2Label: "Redshift Dist", chart1Type: "scatter" } -}; - -document.querySelectorAll('.tree-list li').forEach(item => { - item.addEventListener('click', () => { - const name = item.textContent.trim(); - const cfg = instruments[name] || instruments["Quantum Layer (GPU)"]; - - document.querySelectorAll('.tree-list li').forEach(i => i.classList.remove('active')); - item.classList.add('active'); - - // 1. Reconfigurar 3D - scene.remove(activeViz); - if (cfg.mode === "SPHERE") { - const geom = new THREE.SphereBufferGeometry(20, 32, 32); - const mat = new THREE.MeshPhongMaterial({ color: cfg.color, wireframe: true }); - activeViz = new THREE.Mesh(geom, mat); - } else if (cfg.mode === "VECTORS") { - activeViz = new THREE.Group(); - for(let i=0; i<12; i++) { - const arrow = new THREE.ArrowHelper(new THREE.Vector3(Math.random()-0.5, 1, Math.random()-0.5).normalize(), new THREE.Vector3(Math.random()*40-20, 0, Math.random()*40-20), 10, cfg.color); - activeViz.add(arrow); - } - } else if (cfg.mode === "LATTICE") { - const pointsGeom = new THREE.BufferGeometry(); - const pos = []; for (let i=-40; i<=40; i+=10) for (let j=-40; j<=40; j+=10) pos.push(i, Math.sin(i/10)*5, j); - pointsGeom.setAttribute('position', new THREE.Float32BufferAttribute(pos, 3)); - activeViz = new THREE.Points(pointsGeom, new THREE.PointsMaterial({ color: cfg.color, size: 3 })); - } else { - material.color.setHex(cfg.color); - activeViz = mesh; - } - scene.add(activeViz); - - // 2. Reconfigurar Gráficas - document.getElementById('chart-panel-header').textContent = `DATA STREAM: ${name.toUpperCase()}`; - initCharts(cfg); - }); -}); - -// Inicializar por defecto -initCharts(instruments["Quantum Layer (GPU)"]); - -setInterval(async () => { - if (window.eel) { - const data = await eel.get_latest_data()(); - if (data && !data.error) { - // 1. Actualizar el inspector (Campos filtrados) - const rows = document.querySelectorAll('.stat-row'); - rows.forEach(row => { - const label = row.firstElementChild.textContent.trim(); - if (data[label]) { - row.lastElementChild.textContent = data[label]; - if (label === 'q_snr') row.lastElementChild.className = 'val green'; - } - }); - - // 2. Actualizar Coherencia y Diagnóstico - const coherenceEl = document.querySelector('.stat-value'); - if (coherenceEl && data.coherence) { - coherenceEl.textContent = `${data.coherence} Coherence`; - } - - if (data.calm_index !== undefined) { - bciCalm = parseFloat(data.calm_index); - } - if (data.ratio_enfoque !== undefined) { - bciFocus = parseFloat(data.ratio_enfoque); - } - - const latencyEl = document.getElementById('engine-latency'); - if (latencyEl && data.compute_ms) { - latencyEl.textContent = data.compute_ms; - const ms = parseFloat(data.compute_ms); - latencyEl.style.color = ms > 50 ? '#ff3b30' : 'var(--accent-green)'; - } - - // 3. ACTUALIZACIÓN CIENTÍFICA: Gráficos de Datos Reales - if (charts.length >= 2) { - // Chart 1: Time Series (Strain o Signal) - const val1 = parseFloat(data.q_strain_h1) || Math.random(); - charts[0].data.datasets[0].data.shift(); - charts[0].data.datasets[0].data.push(val1); - charts[0].update('none'); - - // Chart 2: PSD / Spectrum (Capa B: DSP) - if (data.spectrum) { - charts[1].data.datasets[0].data = data.spectrum.slice(0, 30); // Primeros 30 bins - charts[1].update('none'); - } - } - } - } -}, 100); - -// Event Listener para el Toggle de GPU -document.getElementById('gpu-toggle').addEventListener('click', function() { - this.classList.toggle('active'); - const isActive = this.classList.contains('active'); - if (window.eel) { - eel.toggle_gpu(isActive); - } -}); - -// Event Listener para el botón de Record -let isRecording = false; -document.querySelectorAll('.inspector-panel button').forEach(btn => { - if (btn.textContent.includes('Record')) { - btn.id = 'record-btn'; - btn.addEventListener('click', () => { - isRecording = !isRecording; - btn.textContent = isRecording ? 'Recording...' : 'Record Session'; - btn.style.color = isRecording ? 'var(--accent-orange)' : ''; - if (window.eel) eel.toggle_recording(isRecording); - }); - } -}); - -window.addEventListener('resize', () => { - camera.aspect = container.clientWidth / container.clientHeight; - camera.updateProjectionMatrix(); - renderer.setSize(container.clientWidth, container.clientHeight); -}); diff --git a/workbench_launcher.py b/workbench_launcher.py deleted file mode 100644 index c38cf96..0000000 --- a/workbench_launcher.py +++ /dev/null @@ -1,79 +0,0 @@ -import eel -import os -import threading -import time -from pathlib import Path -from substrate_core_emitter import SubstrateEngine - -# 1. Configuración de Rutas -SUBSTRATE_ROOT = Path(__file__).resolve().parent -eel.init('workbench') - -# 2. Orquestación del Motor (Background Thread) -latest_telemetry = {} - -def run_engine_background(): - """Ejecuta el motor de física en un hilo separado""" - engine = SubstrateEngine() - # Modificamos el run del motor para que guarde en la variable global - # en lugar de solo emitir por ZMQ (para redundancia local) - try: - while True: - data = engine.calculate_quantum_fields() - # Intentar leer la telemetria mental en vivo del BCI - import json - calm_index = 0.5 - ratio_enfoque = 1.0 - try: - mind_state_path = Path("C:/Users/Drizzy/Desktop/simulations/mind_state.json") - if mind_state_path.exists(): - with open(mind_state_path, 'r') as f: - m_data = json.load(f) - calm_index = m_data["metrics"]["calm_index_final"] - ratio_enfoque = m_data["metrics"]["ratio_enfoque_final"] - except Exception as e: - pass - - global latest_telemetry - latest_telemetry = { - "q_strain_h1": f"{data.q_strain_h1:.2e}", - "q_snr": f"{data.q_snr:+.1f}", - "final_spin": f"{data.final_spin:.3f}", - "entropy_exp": f"{data.entropy_exp:+.1f}", - "ringdown_sig": f"{data.ringdown_sig:.3f}", - "mass_solar": f"{data.mass_solar:.1f}", - "dist_mpc": f"{data.dist_mpc:.1f}", - "freq_hz": f"{data.freq_hz:.1f}", - "gds_lock_pro": f"{data.gds_lock_pro:.3f}", - "q_phase_ctc": f"{data.q_phase_ctc/3.14:.3f}π", - "coherence": f"{calm_index * 100:.1f}%", - "calm_index": calm_index, - "ratio_enfoque": ratio_enfoque, - "timestamp": data.timestamp - } - # También lo enviamos por ZMQ por si otros procesos escuchan - from dataclasses import asdict - engine.publisher.send_string(f"SUBSTRATE_STATE {json.dumps(asdict(data))}") - - engine.t += 0.05 - time.sleep(0.05) - except Exception as e: - print(f"❌ Error en el motor: {e}") - -# Iniciar el motor inmediatamente -threading.Thread(target=run_engine_background, daemon=True).start() - -# 3. Exposición de Datos al UI -@eel.expose -def get_latest_data(): - """Devuelve la telemetría viva del motor de fondo""" - return latest_telemetry - -# 4. Lanzamiento del Workbench -print("🚀 LANZANDO SISTEMA UNIFICADO SUBSTRATE...") -try: - # Abrir el Workbench en modo App (sin barras de navegador para máxima inmersión) - eel.start('index.html', size=(1440, 900), mode='chrome') -except Exception as e: - print(f"Error al lanzar UI: {e}") - eel.start('index.html', size=(1440, 900))