diff --git a/AGENTS.md b/AGENTS.md index 8d5b878..abcff19 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -11,7 +11,7 @@ Reihenfolge: §1 Working agreements > §2 Conventions > §3 Don't > §4 When stu ## Project - **Was:** Adaptiver RG-QEC-Simulator als **Diagnostik-/Verifikations-Harness mit Konvergenz-Guards** (KEIN Frontier-Threshold-Tool). Spec/Theorie + lauffähiger Phase-1-MVP (seit PR#4). - **Inhalt:** gehärtete Kernel-Spec v1.0 + Proof-Block v1.1 / KernelSpec v1.2 (`spec/`) + Phase-1-MVP-Code (`src/adaptiverg_qec/`, `tests/`, CI). -- **MVP-Stand:** real = A-Kernel-MCMC + Foster-Lyapunov-Drift-Guard + MCRG-Map/Jacobian + **skalarer Swendsen-MCRG-Schätzer** (sample-geschätzte R̂ `T̂=⟨S'S⟩_c/⟨S'S'⟩_c` vs `tanh(2K)` validiert; 1D-Ising-Instanz). **Phase-3a (NEU): T̂ aus dem korrelierten A-Kernel-MCMC (statt exakt-i.i.d.) mit autokorrelations-bewussten Fehlerbalken** — `autocorr.py`: FFT-ρ (Wiener-Khinchin), τ_int + Wolff-g-Windowing, Binning-Cross-Check, Block-Jackknife fürs Verhältnis; belegt N_effσ_iid (`results/phase3a-akernel-autocorr.json`). **Phase-3b (NEU): Multi-Operator-Swendsen-MATRIX auf 2D-Ising** — `ising2d.py` (vektorisierter Checkerboard-Metropolis + Majority-Rule-Blocking b=2) + `mcrg_matrix.py` (gerade Operatoren S₁/S₂/S₃, connected-corr-Matrizen A,B, `T=A·B⁻¹` via `np.linalg.solve`, Eigenwert-Exponenten, Block-Jackknife). y_t=0.97±0.01 vs Onsager-Orakel y_t=1 — **ehrlich GROB** (single-spin + 1 RG-Stufe → kritisches Slowing-Down; Plausibilität, KEIN Frontier-Wert), `results/phase3b-swendsen-matrix.json`. **Phase-4 (NEU): Wolff-Cluster + Multi-RG + y_h** — `wolff2d.py` (Wolff-Single-Cluster, `P_add=1−e^{−2K}`, vektorisierter BFS, rejection-free; Energie vs exakte L=4-Enum |err|<0.008; `τ_int(Wolff)≪τ_int(Metropolis)`, ×12–16 @ L=32) + `mcrg_multirg.py` (iterierte Majority-Stufen L=32→16→8→4: gerader y_t konvergiert, bester `|y_t−1|≈0.006` vs 3b 0.035; ungerade Operatoren O₁=M, O₂=3-Spin-L-Cluster → ungerade Swendsen-Matrix `T_h=A·B⁻¹` → y_h, bester Iterationswert `|y_h−15/8|≈0.002` (Minimum über Iterationen, NICHT tiefste Stufe; tiefste-Iter ≈0.003) vs Onsager y_h=15/8=1.875). Gates G24–G29, `results/phase4-wolff-multirg.json`. **Ehrlich:** Rest-finite-Size bleibt, KEINE L→∞-FSS, KEIN Frontier-Wert. **Phase-5 (NEU): CLT-Varianz σ²_g + R̂-Multichain + Run-Manifest** — `clt.py` (MCMC-CLT `σ²_g=2·τ_int·Var(g)` via Γ-Methode + UNABHÄNGIGER OBM-Schätzer Flegal-Jones; gegen geschlossene AR(1)-Form `σ²_g=σ²_ε/(1−φ)²` validiert, φ=0.8 Orakel 25.0 → Γ rel.0.010/OBM rel.0.020; Coverage beidseitig: korrekt 0.937 vs iid-falsch 0.490), `rhat.py` (rank-normalized split-R̂ + folded-R̂ + bulk/tail-ESS nach Vehtari et al. 2021 doi:10.1214/20-BA1221; A-Kernel M=4 → R̂=1.0003 converged; beidseitig: Mittel-Drift R̂=1.62 bulk / Skalen-Drift R̂=1.27 folded>bulk), `manifest.py` + CLI `phase5 --from-manifest` (Seeds/Parameter/Versionen/git-SHA/Plattform; byte-identische Reproduktion via SHA-256 result_hash; beidseitig: Round-trip==Hash, geänderter Seed→anderer Hash). Gates G33–G38, `results/phase5-clt-rhat-manifest.json`. **Phase-6 (NEU): SNIS + Surrogate-DA + Checkpoint/Restart-Lockfile — schliesst die drei dokumentierten Phase-4/5-Lücken.** `snis.py` (Self-Normalized Importance Sampling auf der offenen 1D-Ising-Kette; χ²-Divergenz GESCHLOSSEN `[cosh(2K_t−K_p)cosh(K_p)/cosh²(K_t)]^{L−1}−1` → Orakel für ESS `1/(1+χ²)`, führenden O(1/N)-Bias `(1+χ²)(tanh K_t−tanh(2K_t−K_p))/N` und MSE-Bound `4(1+χ²)/N` (Agapiou et al. 2017); ESS-Kollaps-Guard beidseitig), `surrogate.py` (Delayed-Acceptance-Metropolis nach Christen & Fox 2005, Surrogat `β̃=β(1+γ)`; γ=0 BIT-IDENTISCH zum Metropolis-A-Kernel, miskalibriertes Surrogat bleibt exakt vs Transfer-Matrix-Orakel; Surrogate-Drift-Guard feuert/hält beidseitig), `checkpoint.py` (Philox-State-Serialisierung + gemeinsamer `advance_chain`/`postprocess_multichain`-Code-Pfad → Interrupt+Resume ergibt BYTE-IDENTISCHEN `result_hash`; O_EXCL-Lockfile + SHA-256-Integritäts-Hash fail-closed). Gates G39–G45, `results/phase6-snis-surrogate-checkpoint.json`, CLI `phase6`. Zusätzlich Audit-Härtung: zell-eigene Seeds in QEC-Sweeps (vorher rangkorrelierte Zellen), Jeffreys-regularisierte `std_err` (Null-Ereignis-Zellen falsifizierbar), `n_sigma` auf SEM statt Einzel-Seed-std (√8 strenger), R̂ erkennt konstante Ketten mit verschiedenen Mitteln (vorher „converged"), Manifest-`__post_init__`-Validierung, rg_map-dtype-Fix. Selftest gesamt **45/45 [PASS]**. Keine neuen Deps. Offen (ehrlich, NICHT erledigt): Defensive Mixture, SNIS auf 2D/RBIM-Targets, DA×Diminishing-Adaptation-Kombination, MMD-Drift; bekannte dokumentierte Limitationen: Majority-Tie-Break nicht Z2-äquivariant (`ising2d.majority_block_b2`), RBIM-Scan nutzt gemeinsame Seeds über p (common random numbers), τ_int-Clamp ≥0.5 (konservativ bei Anti-Korrelation). +- **MVP-Stand:** real = A-Kernel-MCMC + Foster-Lyapunov-Drift-Guard + MCRG-Map/Jacobian + **skalarer Swendsen-MCRG-Schätzer** (sample-geschätzte R̂ `T̂=⟨S'S⟩_c/⟨S'S'⟩_c` vs `tanh(2K)` validiert; 1D-Ising-Instanz). **Phase-3a (NEU): T̂ aus dem korrelierten A-Kernel-MCMC (statt exakt-i.i.d.) mit autokorrelations-bewussten Fehlerbalken** — `autocorr.py`: FFT-ρ (Wiener-Khinchin), τ_int + Wolff-g-Windowing, Binning-Cross-Check, Block-Jackknife fürs Verhältnis; belegt N_effσ_iid (`results/phase3a-akernel-autocorr.json`). **Phase-3b (NEU): Multi-Operator-Swendsen-MATRIX auf 2D-Ising** — `ising2d.py` (vektorisierter Checkerboard-Metropolis + Majority-Rule-Blocking b=2) + `mcrg_matrix.py` (gerade Operatoren S₁/S₂/S₃, connected-corr-Matrizen A,B, `T=A·B⁻¹` via `np.linalg.solve`, Eigenwert-Exponenten, Block-Jackknife). y_t=0.97±0.01 vs Onsager-Orakel y_t=1 — **ehrlich GROB** (single-spin + 1 RG-Stufe → kritisches Slowing-Down; Plausibilität, KEIN Frontier-Wert), `results/phase3b-swendsen-matrix.json`. **Phase-4 (NEU): Wolff-Cluster + Multi-RG + y_h** — `wolff2d.py` (Wolff-Single-Cluster, `P_add=1−e^{−2K}`, vektorisierter BFS, rejection-free; Energie vs exakte L=4-Enum |err|<0.008; `τ_int(Wolff)≪τ_int(Metropolis)`, ×12–16 @ L=32) + `mcrg_multirg.py` (iterierte Majority-Stufen L=32→16→8→4: gerader y_t konvergiert, bester `|y_t−1|≈0.006` vs 3b 0.035; ungerade Operatoren O₁=M, O₂=3-Spin-L-Cluster → ungerade Swendsen-Matrix `T_h=A·B⁻¹` → y_h, bester Iterationswert `|y_h−15/8|≈0.002` (Minimum über Iterationen, NICHT tiefste Stufe; tiefste-Iter ≈0.003) vs Onsager y_h=15/8=1.875). Gates G24–G29, `results/phase4-wolff-multirg.json`. **Ehrlich:** Rest-finite-Size bleibt, KEINE L→∞-FSS, KEIN Frontier-Wert. **Phase-5 (NEU): CLT-Varianz σ²_g + R̂-Multichain + Run-Manifest** — `clt.py` (MCMC-CLT `σ²_g=2·τ_int·Var(g)` via Γ-Methode + UNABHÄNGIGER OBM-Schätzer Flegal-Jones; gegen geschlossene AR(1)-Form `σ²_g=σ²_ε/(1−φ)²` validiert, φ=0.8 Orakel 25.0 → Γ rel.0.010/OBM rel.0.020; Coverage beidseitig: korrekt 0.937 vs iid-falsch 0.490), `rhat.py` (rank-normalized split-R̂ + folded-R̂ + bulk/tail-ESS nach Vehtari et al. 2021 doi:10.1214/20-BA1221; A-Kernel M=4 → R̂=1.0003 converged; beidseitig: Mittel-Drift R̂=1.62 bulk / Skalen-Drift R̂=1.27 folded>bulk), `manifest.py` + CLI `phase5 --from-manifest` (Seeds/Parameter/Versionen/git-SHA/Plattform; byte-identische Reproduktion via SHA-256 result_hash; beidseitig: Round-trip==Hash, geänderter Seed→anderer Hash). Gates G33–G38, `results/phase5-clt-rhat-manifest.json`. **Phase-6 (NEU): SNIS + Surrogate-DA + Checkpoint/Restart-Lockfile — schliesst die drei dokumentierten Phase-4/5-Lücken.** `snis.py` (Self-Normalized Importance Sampling auf der offenen 1D-Ising-Kette; χ²-Divergenz GESCHLOSSEN `[cosh(2K_t−K_p)cosh(K_p)/cosh²(K_t)]^{L−1}−1` → Orakel für ESS `1/(1+χ²)`, führenden O(1/N)-Bias `(1+χ²)(tanh K_t−tanh(2K_t−K_p))/N` und MSE-Bound `4(1+χ²)/N` (Agapiou et al. 2017); ESS-Kollaps-Guard beidseitig), `surrogate.py` (Delayed-Acceptance-Metropolis nach Christen & Fox 2005, Surrogat `β̃=β(1+γ)`; γ=0 BIT-IDENTISCH zum Metropolis-A-Kernel, miskalibriertes Surrogat bleibt exakt vs Transfer-Matrix-Orakel; Surrogate-Drift-Guard feuert/hält beidseitig), `checkpoint.py` (Philox-State-Serialisierung + gemeinsamer `advance_chain`/`postprocess_multichain`-Code-Pfad → Interrupt+Resume ergibt BYTE-IDENTISCHEN `result_hash`; O_EXCL-Lockfile + SHA-256-Integritäts-Hash fail-closed). Gates G39–G45, `results/phase6-snis-surrogate-checkpoint.json`, CLI `phase6`. Zusätzlich Audit-Härtung: zell-eigene Seeds in QEC-Sweeps (vorher rangkorrelierte Zellen), Jeffreys-regularisierte `std_err` (Null-Ereignis-Zellen falsifizierbar), `n_sigma` auf SEM statt Einzel-Seed-std (√8 strenger), R̂ erkennt konstante Ketten mit verschiedenen Mitteln (vorher „converged"), Manifest-`__post_init__`-Validierung, rg_map-dtype-Fix. Selftest gesamt **45/45 [PASS]**. Keine neuen Deps. Offen (ehrlich, NICHT erledigt): Defensive Mixture, SNIS auf 2D/RBIM-Targets, DA×Diminishing-Adaptation-Kombination, MMD-Drift; bekannte dokumentierte Limitationen: τ_int-Clamp ≥0.5 (konservativ bei Anti-Korrelation). **Z2/Seed-Härtung (NEU):** `ising2d.majority_block_b2` löst 2+2-Ties jetzt durch Auswahl eines der vier Original-Spins des Blocks → `B(-s) == -B(s)` gilt konstruktiv exakt statt nur im Mittel; die dadurch verschobenen Phase-3b/4-Referenzwerte (G32, `test_fix3_central_values_unchanged`) wurden neu erhoben, die externen Onsager-Gates G22/G27/G28 blieben PASS. `rbim_scan.py` macht die Stream-Politik explizit (`SeedSequence` + `spawn(2)`, Default `independent`; CRN nur als benannte Option); die historische `rbim_nishimori`-Baseline behält bewusst ihre additiven Seeds. Offen: Regeneration der `results/`-Artefakte. ## Working agreements 1. **Kein Overclaim „implementiert".** `src/` enthält den Phase-1-MVP; README/SOURCES/Status nennen ehrlich, was MVP-real vs. gestubbt ist (stochastische R̂, SNIS, Multichain offen). Keine Komponente als „validiert" behaupten ohne lauffähigen Code + Gate-Log in `results/`. diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index e72532d..1c934b8 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -141,13 +141,28 @@ Duplikat-p-Guard; Surface-Threshold-Fenster verbreitert (0.09–0.115, 80k Shots `crossing_found`-Flag statt stillem NaN. **Bekannte, BEWUSST nicht in diesem Inkrement gefixte Limitationen (dokumentiert, Follow-up):** -- `ising2d.majority_block_b2`-Tie-Break ist deterministisch, aber nicht Z2-äquivariant - (Tie-Break flippt nicht unter s→−s); in-Repo-Aufrufer übergeben `config_index`, der Effekt - ist im Rahmen der ausgewiesenen Grobheit enthalten. Ein äquivarianter Fix ändert alle - committeten Phase-3b/4-Baselines (inkl. G32-Goldwerte) und gehört in ein eigenes Inkrement. -- `rbim_nishimori`-Scan nutzt denselben `base_seed` über alle p (common random numbers über - die Kurve; Disorder-Realisierungen über p genestet) — die p*-Lokalisierung bleibt auf - Plausibilitäts-Niveau, wie ausgewiesen. +- ~~`ising2d.majority_block_b2`-Tie-Break nicht Z2-äquivariant~~ — **erledigt.** Der Tie-Break + wählt jetzt einen der vier Original-Spins des 2x2-Blocks, womit `B(-s) == -B(s)` konstruktiv + exakt gilt (vorher nur im Mittel 50/50). Wie hier vorhergesagt hat der Fix die committeten + Phase-3b/4-Referenzwerte verschoben: die G32-/`test_fix3_central_values_unchanged`-Werte + wurden unter dem äquivarianten Tie-Break neu erhoben (Herleitung im jeweiligen Docstring), + die externen Onsager-Gates G22/G27/G28 blieben unverändert PASS. **Evidenz nachgezogen:** + die drei Artefakte, die die Blocking-Regel berühren, wurden aus diesem Commit neu erzeugt — + `results/phase3b-swendsen-matrix.json` (`python -m adaptiverg_qec.mcrg_matrix`, 18 s), + `results/phase4-wolff-multirg.json` (`python -m adaptiverg_qec.mcrg_multirg`, 135 s) und der + Gate-Log `results/selftest.json` + (`python -m adaptiverg_qec.cli --selftest --json results/selftest.json`, 414 s, 45/45 PASS). + Nur diese drei sind betroffen: `majority_block_b2` wird ausschliesslich von `mcrg_matrix` und + `mcrg_multirg` aufgerufen. Ein Gate-Log, der eine Transformation beschreibt, die es in diesem + Commit nicht mehr gibt, wäre irreführend — auch wenn ihn kein Test liest. +- `rbim_nishimori`-Scan (historische Baseline, bewusst unverändert) nutzt weiterhin + arithmetische Seeds `base_seed + d` / `base_seed + 10000 + d` — also implizite common random + numbers über p. Der neue Pfad `rbim_scan.py` macht die Wahl explizit: Default + `seed_policy="independent"` mischt `p`/`L`/Replikat über `SeedSequence` und trennt Bond- und + Thermal-Strom per `spawn(2)`; CRN ist nur noch als ausdrückliche Option erreichbar. + **Offen an der Baseline:** die additive Ableitung lässt Bond- und MCMC-Seeds ab + `n_disorder >= 10001` exakt überlappen (`base_seed + 10000` tritt in beiden Familien auf); + `n_disorder` ist nur gegen `< 1` geprüft. - `autocorr.integrated_autocorr_time` klemmt τ_int ≥ 0.5 (für anti-korrelierte Reihen bewusst konservativ; jetzt im Code dokumentiert). - G26 vergleicht τ_int in Update-Einheiten (1 Wolff-Cluster vs 1 Metropolis-Sweep), nicht diff --git a/pyproject.toml b/pyproject.toml index 9a75e38..4391f0f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,21 +1,22 @@ [build-system] -requires = ["setuptools>=68"] +requires = ["setuptools>=77.0.3"] build-backend = "setuptools.build_meta" [project] name = "adaptiverg-qec" version = "0.1.0.dev0" -description = "AdaptiveRG-QEC Phase-1 MVP: Diagnostik-/Verifikations-Harness mit Konvergenz-Guards (MCMC + minimaler MCRG)." +description = "AdaptiveRG-QEC diagnostics and verification harness for MCMC, MCRG, and quantum-error-correction experiments." readme = "README.md" requires-python = ">=3.12" -license = { text = "Apache-2.0" } +license = "Apache-2.0" +license-files = ["LICENSE"] authors = [{ name = "Coworker Research" }] keywords = ["qec", "mcmc", "renormalization-group", "foster-lyapunov", "diagnostics"] classifiers = [ "Development Status :: 3 - Alpha", "Intended Audience :: Science/Research", - "License :: OSI Approved :: Apache Software License", "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.12", "Topic :: Scientific/Engineering :: Physics", ] dependencies = [ @@ -25,14 +26,18 @@ dependencies = [ [project.optional-dependencies] dev = ["pytest>=7.4", "ruff>=0.5"] -# Inkrement-3: ECHTES MWPM-Decoding hinter optional-dependency-Gate. -# Ohne dieses Extra importiert surface_decoder.py nicht und seine Tests SKIPPEN -# (sie failen NICHT hart). Kein Eigen-Decoder: Stim (Sampling) + PyMatching 2 (MWPM). +# Real MWPM decoding behind an optional dependency gate. Keep decoder +# dependencies out of the core install so the analytical/MCMC harness remains +# lightweight and independently testable. surface = ["stim>=1.13", "pymatching>=2.1"] [project.scripts] adaptiverg-qec = "adaptiverg_qec.cli:main" +[project.urls] +Repository = "https://github.com/marcohost33-maker/qec-engine" +Issues = "https://github.com/marcohost33-maker/qec-engine/issues" + [tool.setuptools.packages.find] where = ["src"] @@ -48,7 +53,8 @@ src = ["src", "tests"] [tool.ruff.lint] select = ["E", "F", "I", "B", "UP", "SIM"] -# B008: scipy/numpy style; allow defaults. We keep strict otherwise. +# B904 is intentionally ignored for compatibility with existing exception +# chaining style. Keep this comment synchronized with the configured rule. ignore = ["B904"] [tool.ruff.lint.per-file-ignores] diff --git a/results/phase3b-swendsen-matrix.json b/results/phase3b-swendsen-matrix.json index f7a2e8f..a5ba5ac 100644 --- a/results/phase3b-swendsen-matrix.json +++ b/results/phase3b-swendsen-matrix.json @@ -16,68 +16,68 @@ "rows": [ { "seed": 0, - "y_t": 0.9620743073104707, - "y_t_error_jackknife": 0.01282754451200649, - "abs_error_vs_oracle": 0.03792569268952928, - "n_sigma": 2.956582427293945, - "lambda_max": 1.9481088736086845, - "cond_B": 87.26739243966706, - "tau_int_max": 4.292362206104164, + "y_t": 0.9721966879350383, + "y_t_error_jackknife": 0.013094867322879159, + "abs_error_vs_oracle": 0.02780331206496167, + "n_sigma": 2.123222128137498, + "lambda_max": 1.9618254526442591, + "cond_B": 89.95801868899628, + "tau_int_max": 4.1479082343668985, "block_size": 9, "n_blocks": 666, "n_samples": 6000, "eigenvalues_abs": [ - 1.9481088736086845, - 0.20520006968096904 + 1.9618254526442591, + 0.21590292842653014 ], "exponents": [ - 0.9620743073104707, - -2.2848968740380515 + 0.9721966879350383, + -2.211545283218558 ] }, { "seed": 1, - "y_t": 0.9747667403326046, - "y_t_error_jackknife": 0.012729263847867515, - "abs_error_vs_oracle": 0.025233259667395425, - "n_sigma": 1.9823031377908515, - "lambda_max": 1.9653234114638973, - "cond_B": 88.49498716449236, - "tau_int_max": 4.754650650677398, + "y_t": 0.9609243556959352, + "y_t_error_jackknife": 0.013187918611818722, + "abs_error_vs_oracle": 0.03907564430406485, + "n_sigma": 2.962987978182252, + "lambda_max": 1.9465566825425507, + "cond_B": 88.81785823291553, + "tau_int_max": 4.691073292390602, "block_size": 10, "n_blocks": 600, "n_samples": 6000, "eigenvalues_abs": [ - 1.9653234114638973, - 0.18004396177200732 + 1.9465566825425507, + 0.22382184954666418 ], "exponents": [ - 0.9747667403326046, - -2.4735788789618196 + 0.9609243556959352, + -2.1595772154364883 ] }, { "seed": 2, - "y_t": 0.9586219810979993, - "y_t_error_jackknife": 0.013977778172939664, - "abs_error_vs_oracle": 0.04137801890200066, - "n_sigma": 2.960271538870648, - "lambda_max": 1.943452680453265, - "cond_B": 89.74556191382386, - "tau_int_max": 3.5879848820483846, + "y_t": 0.9482400730709694, + "y_t_error_jackknife": 0.013477340065490635, + "abs_error_vs_oracle": 0.05175992692903064, + "n_sigma": 3.8405150183576935, + "lambda_max": 1.9295174256779086, + "cond_B": 93.55665165239814, + "tau_int_max": 3.552955853157631, "block_size": 8, "n_blocks": 750, "n_samples": 6000, "eigenvalues_abs": [ - 1.943452680453265, - 0.24657952643829528 + 1.9295174256779086, + 0.23995835909788088 ], "exponents": [ - 0.9586219810979993, - -2.0198750775154 + 0.9482400730709694, + -2.0591440237836443 ] } ], - "y_t_mean": 0.9651543429136916, - "y_t_multiseed_spread": 0.008501663200437526 + "y_t_mean": 0.9604537055673142, + "y_t_multiseed_spread": 0.011985240197593486 } \ No newline at end of file diff --git a/results/phase4-wolff-multirg.json b/results/phase4-wolff-multirg.json index a3c7dd0..3166c08 100644 --- a/results/phase4-wolff-multirg.json +++ b/results/phase4-wolff-multirg.json @@ -23,37 +23,37 @@ { "seed": 0, "y_t_per_iter": [ - 0.9300108484561757, - 0.9956698050798066, - 1.0219696993830758 + 0.9322035709243033, + 1.012672564728223, + 1.005350865628758 ], "y_t_err_per_iter": [ - 0.023375888974070447, - 0.022268876815738393, - 0.022553185950075737 + 0.021794240268322725, + 0.020336709170486834, + 0.02173754991078718 ], "y_t_abs_err_per_iter": [ - 0.06998915154382435, - 0.004330194920193375, - 0.021969699383075803 + 0.06779642907569672, + 0.01267256472822309, + 0.0053508656287579726 ], - "y_t_best_abs_err": 0.004330194920193375, + "y_t_best_abs_err": 0.0053508656287579726, "y_h_per_iter": [ - 1.8809880939402193, - 1.8728283632124236, - 1.8710143051843333 + 1.8807404284210785, + 1.8730156253034003, + 1.872500085001765 ], "y_h_err_per_iter": [ - 0.0007798049668828549, - 0.0013884662786459575, - 0.002523636665580784 + 0.0007746333092357208, + 0.0013699893997804746, + 0.0025225884859981272 ], "y_h_abs_err_per_iter": [ - 0.0059880939402192634, - 0.00217163678757637, - 0.003985694815666685 + 0.0057404284210784695, + 0.0019843746965997333, + 0.0024999149982349866 ], - "y_h_best_abs_err": 0.00217163678757637, + "y_h_best_abs_err": 0.0019843746965997333, "tau_int_wolff": 2.8048948088192778, "tau_int_metropolis": 45.868154037227725, "tau_ratio_metro_over_wolff": 16.352896334296386, @@ -62,37 +62,37 @@ { "seed": 1, "y_t_per_iter": [ - 0.9511674785369428, - 1.0018521587307372, - 0.9772596907030568 + 0.9283771581485563, + 0.9886259053795388, + 0.9852681391818852 ], "y_t_err_per_iter": [ - 0.02121857283050567, - 0.021384623734982365, - 0.027062828232580285 + 0.02109262529857761, + 0.022542493947599187, + 0.023903640747491645 ], "y_t_abs_err_per_iter": [ - 0.048832521463057166, - 0.0018521587307371856, - 0.02274030929694315 + 0.07162284185144374, + 0.01137409462046124, + 0.01473186081811484 ], - "y_t_best_abs_err": 0.0018521587307371856, + "y_t_best_abs_err": 0.01137409462046124, "y_h_per_iter": [ - 1.8818738251638143, - 1.8773650370592705, - 1.8797913359995964 + 1.8816035697812679, + 1.8773449589771227, + 1.8792501762475406 ], "y_h_err_per_iter": [ - 0.0007371943566163144, - 0.0013432812124068113, - 0.0023584694388396154 + 0.0007231482535378979, + 0.0013925884366575037, + 0.0023015727701970394 ], "y_h_abs_err_per_iter": [ - 0.0068738251638142955, - 0.0023650370592704917, - 0.004791335999596358 + 0.0066035697812678595, + 0.002344958977122724, + 0.004250176247540649 ], - "y_h_best_abs_err": 0.0023650370592704917, + "y_h_best_abs_err": 0.002344958977122724, "tau_int_wolff": 2.2252076610591995, "tau_int_metropolis": 26.997274792427625, "tau_ratio_metro_over_wolff": 12.1324743145891, @@ -101,44 +101,44 @@ { "seed": 2, "y_t_per_iter": [ - 0.9839328834061714, - 1.0122269043974723, - 0.9856284941274613 + 0.9830995508507281, + 1.0104677513182485, + 1.0044551296508368 ], "y_t_err_per_iter": [ - 0.020402384021764607, - 0.021701359113449857, - 0.023390115584646282 + 0.020645081764230314, + 0.021203789863875146, + 0.02080542949305262 ], "y_t_abs_err_per_iter": [ - 0.016067116593828645, - 0.012226904397472316, - 0.014371505872538659 + 0.016900449149271912, + 0.010467751318248508, + 0.004455129650836831 ], - "y_t_best_abs_err": 0.012226904397472316, + "y_t_best_abs_err": 0.004455129650836831, "y_h_per_iter": [ - 1.881775212139035, - 1.877433948460453, - 1.874656697234665 + 1.8809373380169399, + 1.876009294189666, + 1.8759663640446882 ], "y_h_err_per_iter": [ - 0.0007649463551989142, - 0.0013849452387626887, - 0.002429899789878017 + 0.0007748360631930655, + 0.0013650956350285649, + 0.002447779772077367 ], "y_h_abs_err_per_iter": [ - 0.006775212139034981, - 0.0024339484604529726, - 0.00034330276533500736 + 0.005937338016939853, + 0.001009294189665999, + 0.0009663640446881949 ], - "y_h_best_abs_err": 0.00034330276533500736, + "y_h_best_abs_err": 0.0009663640446881949, "tau_int_wolff": 2.7994646805179157, "tau_int_metropolis": 32.91937002274189, "tau_ratio_metro_over_wolff": 11.759166047650092, "mean_cluster_fraction": 0.4456201171875 } ], - "y_t_best_abs_err_mean": 0.006136419349467626, - "y_h_best_abs_err_mean": 0.0016266588707272895, + "y_t_best_abs_err_mean": 0.007060029966685348, + "y_h_best_abs_err_mean": 0.0017652325728035507, "tau_ratio_mean": 13.41484556551186 } \ No newline at end of file diff --git a/results/selftest.json b/results/selftest.json index 2198e70..adce976 100644 --- a/results/selftest.json +++ b/results/selftest.json @@ -4,7 +4,7 @@ "n_pass": 45, "n_total": 45, "all_pass": true, - "elapsed_s": 169.919, + "elapsed_s": 411.083, "gates": [ { "gate": "G1 Analytic-Oracle (MCMC vs Transfer-Matrix)", @@ -24,7 +24,7 @@ { "gate": "G4 Jacobian-Consistency (CS==FD==analytic)", "pass": true, - "detail": "max|CS-analytic|=1.11e-16 max|CS-FD|=2.12e-11" + "detail": "max|CS-analytic|=1.11e-16 max|CS-FD|=1.18e-11" }, { "gate": "G5 RG-Fixpoint (R(0)=0, iter->0)", @@ -104,17 +104,17 @@ { "gate": "G20 corr-matrices A,B symmetric/PSD/conditioned", "pass": true, - "detail": "B symmetric(max|B-B.T|=0.0e+00) PSD(min_eig=1.16e+01)=True cond(B)=87.1 A_finite=True" + "detail": "B symmetric(max|B-B.T|=0.0e+00) PSD(min_eig=1.12e+01)=True cond(B)=89.9 A_finite=True" }, { "gate": "G21 T=A.B^-1 linear-solve consistency (resid~eps)", "pass": true, - "detail": "max|T@B - A|=2.27e-13 rel=1.72e-16(<1e-9) (linear-solve, no explicit inverse)" + "detail": "max|T@B - A|=4.55e-13 rel=3.45e-16(<1e-9) (linear-solve, no explicit inverse)" }, { "gate": "G22 y_t matrix vs Onsager y_t=1 (honest band)", "pass": true, - "detail": "y_t per seed=[0.962 0.975 0.959] mean=0.9652 oracle=1.0 |err|=0.035(<0.2) multiseed_spread=0.009 jackknife_sigma~0.013 lambda_max~1.948 tau~4.3 (COARSE by design)" + "detail": "y_t per seed=[0.972 0.961 0.948] mean=0.9605 oracle=1.0 |err|=0.040(<0.2) multiseed_spread=0.012 jackknife_sigma~0.013 lambda_max~1.962 tau~4.1 (COARSE by design)" }, { "gate": "G23 Swendsen-matrix reproducibility (seed)", @@ -139,12 +139,12 @@ { "gate": "G27 y_t multi-RG converges (beats Phase-3b)", "pass": true, - "detail": "y_t per iter=[0.9300 0.9957 1.0220] oracle=1.0 |err|=[0.0700 0.0043 0.0220] best|err|=0.0043(<0.035 Phase-3b) improved-vs-iter1=True L=32" + "detail": "y_t per iter=[0.9322 1.0127 1.0054] oracle=1.0 |err|=[0.0678 0.0127 0.0054] best|err|=0.0054(<0.035 Phase-3b) improved-vs-iter1=True L=32" }, { "gate": "G28 y_h matrix vs Onsager 15/8 (best iter)", "pass": true, - "detail": "y_h per iter=[1.8810 1.8728 1.8710] oracle=15/8=1.8750 |err|=[0.0060 0.0022 0.0040] best|err|=0.0022(<0.05) sigma_jk~0.0008 L=32 (best=min over iters, NOT deepest; deepest|err|=0.0040)" + "detail": "y_h per iter=[1.8807 1.8730 1.8725] oracle=15/8=1.8750 |err|=[0.0057 0.0020 0.0025] best|err|=0.0020(<0.05) sigma_jk~0.0008 L=32 (best=min over iters, NOT deepest; deepest|err|=0.0025)" }, { "gate": "G29 Phase-4 reproducibility (seed)", diff --git a/src/adaptiverg_qec/cli.py b/src/adaptiverg_qec/cli.py index b24a639..4967265 100644 --- a/src/adaptiverg_qec/cli.py +++ b/src/adaptiverg_qec/cli.py @@ -978,16 +978,45 @@ def _g32_jackknife_block_per_iter() -> tuple[bool, str]: """Codex-Fix 3: Jackknife-Blockgroesse PRO ITERATION (nicht global Level 0). Verifiziert (a) per-iter Blockgroessen werden gemeldet, (b) die y_t/y_h- - ZENTRALWERTE bleiben byte-identisch zur Baseline (nur Fehlerbalken aendern), - (c) per-iter Blockgroessen koennen zwischen Stufen variieren. + ZENTRALWERTE reproduzieren die committete Referenz byte-genau (nur + Fehlerbalken haengen an der Blockgroesse), (c) per-iter Blockgroessen + koennen zwischen Stufen variieren. + + REFERENZ-HERKUNFT (wichtig -- die Zahlen sind KEIN frei nachgezogener + Messwert). Die urspruengliche Referenz stammte aus der Zeit des ALTEN, + NICHT Z2-aequivarianten Majority-Tie-Breaks. Der Tie-Break in + `ising2d.majority_block_b2` waehlt seither einen der vier Original-Spins + des 2x2-Blocks, womit `B(-s) == -B(s)` konstruktiv exakt gilt (vorher nur + im Mittel 50/50). Die Blocking-Abbildung selbst ist damit bewusst eine + andere; y_t/y_h sind Funktionale der geblockten Konfigurationen und MUESSEN + sich verschieben. Die alte Referenz haette behauptet, eine absichtliche + Aenderung der RG-Abbildung sei wirkungslos -- sie war ab dem Fix falsch. + + Dass die NEUE Referenz die richtige ist, haengt nicht daran, dass sie + herauskam, sondern an drei unabhaengigen Belegen: + 1. Symmetrie-Orakel: die erzeugende Abbildung ist jetzt exakt + Z2-aequivariant (erschoepfend geprueft, tests/test_ising2d.py). + Exakte Spin-Flip-Symmetrie ist eine Forderung an eine legitime + Ising-Blocking-Regel, keine freie Wahl. + 2. Externes Physik-Orakel (Onsager, unabhaengig von diesem Code): + y_t=1 und y_h=15/8 werden mit den neuen Zentralwerten weiterhin in + den ausgewiesenen ehrlichen Banden getroffen -- G22/G27/G28 bleiben + unveraendert PASS (best |y_t-1|=0.0054 < 0.035, best |y_h-15/8|=0.0020 + < 0.05); der ungerade Sektor wird sogar leicht besser. + 3. Struktur-Invariante: dass die Zentralwerte ueberhaupt nicht an der + Jackknife-Blockgroesse haengen -- die eigentliche Fix-3-Behauptung -- + wird snapshot-frei in + tests/test_mcrg_multirg.py::test_jackknife_block_size_moves_only_error_bars + geprueft (per-iter-Default vs global fixe Blockgroesse). """ v = mcrg_multirg.validate_multirg_2d( L=32, n_op_even=2, n_op_odd=2, n_levels=3, n_records=3000, burn_in=400, seed=0 ) bsz_t = np.asarray(v.multirg.block_size_per_iter) bsz_h = np.asarray(v.multirg_odd.block_size_per_iter) - base_yt = np.array([0.93001085, 0.99566981, 1.0219697]) - base_yh = np.array([1.88098809, 1.87282836, 1.87101431]) + # Referenz unter dem Z2-aequivarianten Tie-Break (siehe Docstring). + base_yt = np.array([0.93220357, 1.01267256, 1.00535087]) + base_yh = np.array([1.88074043, 1.87301563, 1.87250009]) yt_same = np.allclose(v.multirg.y_t_per_iter, base_yt, rtol=0, atol=1e-7) yh_same = np.allclose(v.multirg_odd.y_h_per_iter, base_yh, rtol=0, atol=1e-7) sized = bsz_t.shape[0] == v.multirg.n_iters and bsz_h.shape[0] == v.multirg_odd.n_iters diff --git a/src/adaptiverg_qec/ising2d.py b/src/adaptiverg_qec/ising2d.py index ff9347d..0f02223 100644 --- a/src/adaptiverg_qec/ising2d.py +++ b/src/adaptiverg_qec/ising2d.py @@ -1,56 +1,28 @@ -"""Phase-3b: vektorisierter 2D-Ising-Sampler + Majority-Rule-Blocking. - -KONTEXT (Phase-3b). Phase-2/3a war SKALAR (1D-Ising, T_c=0 -> kein nicht-trivialer -Fixpunkt -> nur die Schaetz-Maschinerie validierbar, KEIN echter Exponent). Das -2D-Ising-Modell auf dem Quadratgitter hat dagegen einen NICHT-TRIVIALEN Fixpunkt -bei T_c = 2 / ln(1+sqrt(2)) ~ 2.2692 J (Onsager 1944), also echte messbare -RG-Exponenten -- die Voraussetzung fuer die Multi-Operator-Swendsen-MATRIX -(mcrg_matrix.py). - ------------------------------------------------------------------------------- -KONVENTION ------------------------------------------------------------------------------- -Spin s_ij in {+1,-1} auf dem L x L Quadratgitter mit PERIODISCHEN Randbedingungen. -Hamiltonian (ferromagnetisch, kein Feld): - H(s) = - J * sum_ s_i s_j (Summe ueber NN-Bonds). -Boltzmann-Ziel pi ~ exp(-H/T) = exp(K * sum_ s_i s_j) mit der dimensionslosen -Kopplung K = J/T. Wir setzen J=1, parametrisieren also direkt ueber K (bzw. T=1/K). -T_c entspricht K_c = ln(1+sqrt(2))/2 ~ 0.4407. - ------------------------------------------------------------------------------- -SAMPLER: Checkerboard-(Schachbrett-)Metropolis -- VEKTORISIERT ------------------------------------------------------------------------------- -Single-Spin-Metropolis, aber OHNE Python-Doppelschleife ueber Spins: das Gitter -wird in zwei interpenetrierende Untergitter (schwarz/weiss wie ein Schachbrett) -zerlegt. Innerhalb eines Untergitters sind alle Spins UNABHAENGIG (ihre Nachbarn -liegen komplett im anderen Untergitter), also kann ein ganzes Untergitter in EINEM -vektorisierten numpy-Schritt aktualisiert werden. Ein "Sweep" = beide Untergitter -einmal. Das ist exakte Single-Spin-Metropolis-Dynamik (detailed balance pro Spin), -nur datenparallel ausgewertet. - -EHRLICHE SCOPE-GRENZE (AGENTS.md Sec.1): Single-Spin-Metropolis nahe T_c leidet -unter KRITISCHEM SLOWING-DOWN (tau_int ~ L^z, z~2.17). Bei kleinem L + 1 RG-Stufe -ist die Exponenten-Praezision daher nur GROB; ein Cluster-Algorithmus (Wolff/ -Swendsen-Wang) + mehrere RG-Iterationen waeren fuer Hochpraezision noetig (= Phase-4). -Hier wird die MATRIX-MASCHINERIE auf Plausibilitaets-Niveau validiert, KEIN -Frontier-Wert. Das ist bewusst und ehrlich. - ------------------------------------------------------------------------------- -BLOCKING: Majority-Rule (Kadanoff), b=2 ------------------------------------------------------------------------------- -2x2-Bloecke -> ein Block-Spin = Vorzeichen der Block-Summe (Mehrheitsregel). -Bei einem GLEICHSTAND (Summe 0, zwei +1 und zwei -1) gibt es keine Mehrheit; die -Wahl muss UNBIASED sein (sonst bricht die +/- -Symmetrie und verfaelscht den -geraden/ungeraden Sektor). Wir brechen Gleichstaende deterministisch+reproduzierbar -und im Mittel unbiased ueber einen pro-Block- und pro-Konfigurations-Hash (Splitmix- -artig aus Block-Index + globalem Seed): das Vorzeichen ist eine feste Funktion von -(config_index, block_row, block_col, seed), je zur Haelfte +1/-1. Dokumentiert & -getestet (test_ising2d: Tie-Break im Mittel symmetrisch). - -Referenzen (Methode, NICHT Dependency): -- L. Onsager, Phys. Rev. 65 (1944) 117 (exakte 2D-Loesung, T_c, Exponenten). -- L. P. Kadanoff, Physics 2 (1966) 263 (Block-Spin / Majority-Rule). -- R. H. Swendsen, Phys. Rev. Lett. 42 (1979) 859 (MCRG-Matrix). +"""Phase-3b: vectorized 2D Ising sampler and majority-rule blocking. + +The module provides the homogeneous 2D Ising reference model used by the MCRG +experiments. The blocking rule is deliberately deterministic for reproducible +runs, but its tie rule must also preserve the exact global spin-flip symmetry. +In particular, for every configuration ``s`` and fixed tie-breaking metadata, + + majority_block_b2(-s) == -majority_block_b2(s) + +must hold. This is stronger than merely obtaining a 50/50 tie distribution in +aggregate and prevents artificial coupling of even and odd RG sectors. + +Conventions +----------- +Spins ``s_ij`` are in ``{+1, -1}`` on an ``L x L`` square lattice with periodic +boundary conditions. The Hamiltonian is ``H = -J sum_ s_i s_j`` and the +Boltzmann target is proportional to ``exp(K sum_ s_i s_j)`` with ``K=J/T``. +For ``J=1`` the critical coupling is +``K_c = log(1 + sqrt(2))/2``. + +References +---------- +- L. Onsager, Phys. Rev. 65 (1944) 117. +- L. P. Kadanoff, Physics 2 (1966) 263. +- R. H. Swendsen, Phys. Rev. Lett. 42 (1979) 859. """ from __future__ import annotations @@ -71,16 +43,15 @@ "exact_energy_per_spin_2x2", ] -# Onsager: K_c = ln(1+sqrt(2))/2 ; T_c = 1/K_c = 2/ln(1+sqrt(2)). -KC_2D: float = 0.5 * math.log(1.0 + math.sqrt(2.0)) # ~ 0.4406868 -TC_2D: float = 1.0 / KC_2D # ~ 2.2691853 +KC_2D: float = 0.5 * math.log(1.0 + math.sqrt(2.0)) +TC_2D: float = 1.0 / KC_2D @dataclass(frozen=True) class Ising2DChain: - """Resultat eines 2D-Ising-Metropolis-Laufs (Konfigurations-Trajektorie). + """Result of a 2D-Ising Metropolis run. - configs: (n_records, L, L) in {+1,-1}, je ein aufgezeichneter Sweep nach Burn-in. + ``configs`` has shape ``(n_records, L, L)`` and values in ``{+1, -1}``. """ configs: np.ndarray @@ -92,10 +63,20 @@ class Ising2DChain: def __post_init__(self) -> None: if self.configs.ndim != 3: raise ValueError(f"configs must be 3D (n,L,L), got ndim={self.configs.ndim}") + if self.configs.shape[1:] != (self.L, self.L): + raise ValueError( + f"configs trailing shape must be ({self.L},{self.L}), got {self.configs.shape[1:]}" + ) + if not np.all((self.configs == 1) | (self.configs == -1)): + raise ValueError("configs must contain only +/-1 spins") + if not np.isfinite(self.K) or self.K <= 0.0: + raise ValueError(f"K must be finite and > 0, got {self.K}") + if not np.isfinite(self.acceptance) or not (0.0 <= self.acceptance <= 1.0): + raise ValueError(f"acceptance must be in [0,1], got {self.acceptance}") def _neighbor_field(s: np.ndarray) -> np.ndarray: - """Summe der 4 NN-Spins je Gitterplatz (periodisch), vektorisiert.""" + """Return the sum of the four periodic nearest-neighbour spins.""" return ( np.roll(s, 1, axis=0) + np.roll(s, -1, axis=0) @@ -105,7 +86,7 @@ def _neighbor_field(s: np.ndarray) -> np.ndarray: def _checkerboard_masks(L: int) -> tuple[np.ndarray, np.ndarray]: - """Schwarz/Weiss-Untergitter-Masken (Schachbrett): (i+j) gerade/ungerade.""" + """Return even/odd checkerboard sub-lattice masks.""" ii, jj = np.indices((L, L)) even = (ii + jj) % 2 == 0 return even, ~even @@ -120,39 +101,27 @@ def checkerboard_metropolis( seed: int, record_every: int = 1, ) -> Ising2DChain: - """Vektorisierter Checkerboard-Metropolis fuer das 2D-Ising-Modell. - - Ein Sweep = beide Untergitter (schwarz, weiss) je einmal vektorisiert - aktualisiert. Metropolis-Akzeptanz: ein Flip s->-s aendert die lokale Energie - um dE = +2*K*s*nbsum (in pi-Konvention pi~exp(K*sum bonds), also Flip-Faktor - exp(-dlogpi) mit dlogpi = -2*K*s*nbsum). Akzeptiere mit min(1, exp(-2*K*s*nbsum)). - - Args: - K: Kopplung J/T (J=1). Endlich, > 0 (ferromagnetisch). - L: Gitterkante (>=4, gerade -> b=2-Blocking ergibt L/2 sauber). - n_sweeps: Anzahl Sweeps NACH burn_in, die fuer Records zaehlen. - burn_in: Anzahl Sweeps Equilibrierung (verworfen). - seed: PCG64-Seed (Reproduzierbarkeit). - record_every: zeichne jeden record_every-ten Sweep auf (Ausduennung gegen - Autokorrelation -- der Record-Abstand reduziert tau_int der Reihe). - - Returns: - Ising2DChain mit configs (n_records, L, L) in {+1,-1}. + """Run vectorized checkerboard single-spin Metropolis dynamics. + + A sweep updates both independent checkerboard sub-lattices once. The flip + log-probability ratio is ``-2*K*s*sum_neighbours``. """ if not np.isfinite(K): raise ValueError(f"K must be finite, got {K}") if K <= 0.0: raise ValueError(f"K must be > 0 (ferromagnetic), got {K}") + if not isinstance(L, (int, np.integer)): + raise TypeError(f"L must be an integer, got {type(L).__name__}") if L < 4: raise ValueError(f"L must be >= 4, got {L}") if L % 2 != 0: raise ValueError(f"L must be even for b=2 majority blocking, got {L}") - if n_sweeps < 1: - raise ValueError(f"n_sweeps must be >= 1, got {n_sweeps}") - if burn_in < 0: - raise ValueError(f"burn_in must be >= 0, got {burn_in}") - if record_every < 1: - raise ValueError(f"record_every must be >= 1, got {record_every}") + if not isinstance(n_sweeps, (int, np.integer)) or n_sweeps < 1: + raise ValueError(f"n_sweeps must be an integer >= 1, got {n_sweeps}") + if not isinstance(burn_in, (int, np.integer)) or burn_in < 0: + raise ValueError(f"burn_in must be an integer >= 0, got {burn_in}") + if not isinstance(record_every, (int, np.integer)) or record_every < 1: + raise ValueError(f"record_every must be an integer >= 1, got {record_every}") rng = np.random.default_rng(seed) s = np.where(rng.random((L, L)) < 0.5, 1.0, -1.0) @@ -166,11 +135,9 @@ def checkerboard_metropolis( for sweep in range(total): for mask in (even, odd): nb = _neighbor_field(s) - # dlogpi for a flip = -2*K*s*nb; accept with prob min(1, exp(dlogpi)). dlogpi = -2.0 * K * s * nb - accept_prob = np.exp(np.minimum(dlogpi, 0.0)) # min(1, exp(dlogpi)) - r = rng.random((L, L)) - flip = mask & (r < accept_prob) + accept_prob = np.exp(np.minimum(dlogpi, 0.0)) + flip = mask & (rng.random((L, L)) < accept_prob) if sweep >= burn_in: n_accept += int(flip.sum()) n_attempt += int(mask.sum()) @@ -179,25 +146,23 @@ def checkerboard_metropolis( records.append(s.copy()) configs = np.asarray(records, dtype=np.int8) - acceptance = (n_accept / n_attempt) if n_attempt > 0 else 0.0 + acceptance = n_accept / n_attempt if n_attempt else 0.0 return Ising2DChain( - configs=configs, K=float(K), L=int(L), acceptance=acceptance, seed=int(seed) + configs=configs, + K=float(K), + L=int(L), + acceptance=float(acceptance), + seed=int(seed), ) def energy_per_spin(s: np.ndarray) -> np.ndarray: - """Energie pro Spin E/N = -(1/N) sum_ s_i s_j (J=1), periodisch. - - Args: - s: (...,L,L) in {+1,-1}. - - Returns: - E/N je Konfiguration: Shape (...,) (Skalar fuer 2D-Eingabe). - """ + """Return periodic nearest-neighbour Ising energy per spin, with ``J=1``.""" s = np.asarray(s, dtype=np.float64) + if s.ndim < 2: + raise ValueError(f"s must have at least two dimensions, got ndim={s.ndim}") if s.shape[-1] < 2 or s.shape[-2] < 2: raise ValueError(f"need L>=2 in both dims, got {s.shape[-2:]}") - # Jeder Bond einmal: rechts + unten (periodisch). right = s * np.roll(s, -1, axis=-1) down = s * np.roll(s, -1, axis=-2) bonds = right.sum(axis=(-1, -2)) + down.sum(axis=(-1, -2)) @@ -206,84 +171,93 @@ def energy_per_spin(s: np.ndarray) -> np.ndarray: def magnetization_per_spin(s: np.ndarray) -> np.ndarray: - """Magnetisierung pro Spin m = (1/N) sum_ij s_ij (vorzeichenbehaftet).""" + """Return signed magnetization per spin.""" s = np.asarray(s, dtype=np.float64) + if s.ndim < 2: + raise ValueError(f"s must have at least two dimensions, got ndim={s.ndim}") n = s.shape[-1] * s.shape[-2] return s.sum(axis=(-1, -2)) / n def exact_energy_per_spin_2x2(K: float) -> float: - """Exakte E/N der 2x2-periodischen 2D-Ising-Plaquette (Lehrbuch-Orakel). - - Das 2x2-Gitter mit periodischen Randbedingungen hat DOPPELTE Bonds (jeder der - 4 Spins ist mit jedem Nachbar zweimal verbunden via wrap). Statt eine - Konvention zu erraten, enumerieren wir die 2^4=16 Zustaende mit GENAU der - energy_per_spin-Bond-Konvention (right+down, periodisch) -> exaktes, - selbstkonsistentes Orakel fuer den Sampler-Test. - - Args: - K: Kopplung. - - Returns: - exakt fuer L=2. - """ + """Return exact ```` for the periodic 2x2 lattice by enumeration.""" if not np.isfinite(K) or K <= 0.0: raise ValueError(f"K must be finite and > 0, got {K}") states = np.arange(16, dtype=np.int64) bits = ((states[:, None] >> np.arange(4)[None, :]) & 1).astype(np.int8) s = (1 - 2 * bits).reshape(16, 2, 2).astype(np.float64) - e = energy_per_spin(s) # (16,) E/N per state - # Boltzmann-Gewicht pi ~ exp(K * sum bonds) = exp(-K*N*e), N=4. - w = np.exp(-K * 4.0 * e) - w = w / w.sum() + e = energy_per_spin(s) + logw = -K * 4.0 * e + logw -= np.max(logw) + w = np.exp(logw) + w /= w.sum() return float((w * e).sum()) -def majority_block_b2(s: np.ndarray, *, config_index: int = 0, seed: int = 0) -> np.ndarray: - """Majority-Rule-Blocking b=2: L x L -> L/2 x L/2 Block-Spins. +def _splitmix64_grid(config_index: int, br: np.ndarray, bc: np.ndarray, seed: int) -> np.ndarray: + """Return a deterministic uint64 hash for each block coordinate.""" + mask = (1 << 64) - 1 + ci_u = np.uint64(int(config_index) & mask) + seed_u = np.uint64(int(seed) & mask) + with np.errstate(over="ignore"): + h = ( + (ci_u * np.uint64(0x9E3779B97F4A7C15)) + ^ (br.astype(np.uint64) * np.uint64(0xBF58476D1CE4E5B9)) + ^ (bc.astype(np.uint64) * np.uint64(0x94D049BB133111EB)) + ^ seed_u + ) + h ^= h >> np.uint64(30) + h *= np.uint64(0xBF58476D1CE4E5B9) + h ^= h >> np.uint64(27) + h *= np.uint64(0x94D049BB133111EB) + h ^= h >> np.uint64(31) + return h + - Jeder 2x2-Block -> Vorzeichen der 4er-Summe. Gleichstand (Summe 0) wird - UNBIASED + reproduzierbar gebrochen: ein deterministischer Splitmix-Hash aus - (config_index, block_row, block_col, seed) liefert je zur Haelfte +1/-1 - (im Mittel symmetrisch -> +/- -Symmetrie bleibt erhalten). +def majority_block_b2(s: np.ndarray, *, config_index: int = 0, seed: int = 0) -> np.ndarray: + """Apply a ``2x2 -> 1`` majority-rule block-spin transformation. - Args: - s: (L,L) in {+1,-1} (eine Konfiguration). L gerade. - config_index: laufender Index der Konfiguration (entkoppelt Tie-Breaks - ueber die Trajektorie -> kein systematischer Bias). - seed: globaler Seed fuer die Tie-Break-Hash-Funktion. + Non-tied blocks use the sign of the four-spin sum. A two-vs-two tie is + resolved by selecting one of the *four input spins* using a deterministic + SplitMix64-style block hash. Selecting an input spin, instead of generating + an independent ``+/-1`` hash bit, enforces exact global-spin-flip + equivariance while retaining reproducibility and an unbiased selector: - Returns: - Block-Spins (L/2, L/2) in {+1,-1}. + ``majority_block_b2(-s, ...) == -majority_block_b2(s, ...)``. """ - s = np.asarray(s, dtype=np.int64) - if s.ndim != 2: - raise ValueError(f"s must be 2D (L,L), got ndim={s.ndim}") - L = s.shape[0] - if s.shape[1] != L: - raise ValueError(f"s must be square, got {s.shape}") + # Die Spin-Pruefung laeuft auf der UNVERAENDERTEN Eingabe. Ein vorgezogener + # Cast nach int64 wuerde genau die Werte unsichtbar machen, gegen die sie + # schuetzt: 1.5 und -1.5 wuerden zu 1 und -1 abgeschnitten und danach als + # gueltige Spins durchgehen, NaN/inf wuerden ueber einen undefinierten Cast + # (RuntimeWarning) laufen statt ueber die dokumentierte Fehlermeldung. Erst + # pruefen, dann casten. + s_in = np.asarray(s) + if s_in.ndim != 2: + raise ValueError(f"s must be 2D (L,L), got ndim={s_in.ndim}") + L = s_in.shape[0] + if s_in.shape[1] != L: + raise ValueError(f"s must be square, got {s_in.shape}") if L < 2 or L % 2 != 0: raise ValueError(f"L must be even and >= 2, got {L}") + if not np.all((s_in == 1) | (s_in == -1)): + raise ValueError("s must contain only +/-1 spins") + s = s_in.astype(np.int64) + lb = L // 2 - block_sum = s.reshape(lb, 2, lb, 2).sum(axis=(1, 3)) # (lb, lb), in {-4,-2,0,2,4} + block_view = s.reshape(lb, 2, lb, 2) + block_sum = block_view.sum(axis=(1, 3)) - # Unbiased deterministic tie-break for block_sum == 0. - # Splitmix64-artiger Hash; das niederwertige Bit entscheidet +/- (50/50). - # uint64-Arithmetik ist BEWUSST modular (2^64 wraparound) -- numpy warnt sonst - # vor "overflow"; wir unterdruecken die Warnung lokal (kein Bug, gewollte Mod-Arith). br, bc = np.indices((lb, lb)) - with np.errstate(over="ignore"): - h = ( - (np.uint64(config_index) * np.uint64(0x9E3779B97F4A7C15)) - ^ (br.astype(np.uint64) * np.uint64(0xBF58476D1CE4E5B9)) - ^ (bc.astype(np.uint64) * np.uint64(0x94D049BB133111EB)) - ^ np.uint64(seed) - ) - h = h ^ (h >> np.uint64(30)) - h = h * np.uint64(0xBF58476D1CE4E5B9) - h = h ^ (h >> np.uint64(27)) - tie = np.where((h & np.uint64(1)) == np.uint64(1), 1, -1) + h = _splitmix64_grid(config_index, br, bc, seed) + selector = (h & np.uint64(3)).astype(np.intp) + + # Convert (block-row, intra-row, block-col, intra-col) to + # (block-row, block-col, flat-intra-block-index) and select the same + # physical input position for s and -s. The selected spin therefore flips + # exactly under the global Z2 transformation. + blocks = block_view.transpose(0, 2, 1, 3).reshape(lb, lb, 4) + tie_spin = np.take_along_axis(blocks, selector[..., None], axis=2)[..., 0] out = np.sign(block_sum).astype(np.int64) - out = np.where(block_sum == 0, tie, out) + out = np.where(block_sum == 0, tie_spin, out) return out.astype(np.int8) diff --git a/src/adaptiverg_qec/rbim_scan.py b/src/adaptiverg_qec/rbim_scan.py new file mode 100644 index 0000000..ba17c5b --- /dev/null +++ b/src/adaptiverg_qec/rbim_scan.py @@ -0,0 +1,187 @@ +"""Reproducible RBIM scan orchestration with explicit random-stream policy. + +The historical :mod:`rbim_nishimori` baseline used arithmetic seeds of the form +``base_seed + d`` at every p point. Reusing those seeds can be a deliberate +common-random-numbers (CRN) design, but it must not happen implicitly because it +correlates scan points. + +This module keeps the historical baseline untouched and provides a versionable +scan path with two explicit policies: + +``independent`` (default) + Mix ``p``, lattice size, replicate id and the root seed through NumPy + ``SeedSequence``. Different p points therefore receive independent (with + very high probability) disorder and thermal streams. + +``common_random_numbers`` + Deliberately omit ``p`` from the stream identity. The same underlying + streams are then reused across p values. This is useful only when CRN is an + intentional variance-reduction design and the induced covariance is handled + in the analysis. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from typing import Literal + +import numpy as np + +from .rbim_nishimori import ( + DisorderResult, + nishimori_beta, + rbim_wolff_sample, + sample_bonds, +) + +SeedPolicy = Literal["independent", "common_random_numbers"] + + +@dataclass(frozen=True) +class StreamSeeds: + """Concrete child seeds used for one disorder replicate.""" + + bond_seed: int + thermal_seed: int + policy: SeedPolicy + p: float + L: int + replicate: int + + +def _float64_words(value: float) -> tuple[int, int]: + """Encode a finite float64 exactly as two uint32 words for SeedSequence.""" + if not np.isfinite(value): + raise ValueError(f"value must be finite, got {value}") + bits = int(np.asarray(value, dtype=np.float64).view(np.uint64)) + return bits & 0xFFFFFFFF, bits >> 32 + + +def derive_stream_seeds( + *, + base_seed: int, + p: float, + L: int, + replicate: int, + policy: SeedPolicy = "independent", +) -> StreamSeeds: + """Derive reproducible disorder/thermal child streams for one scan cell. + + ``SeedSequence`` hashes the complete stream identity instead of relying on + neighbouring integer seeds. Child streams are created with ``spawn(2)`` so + bond generation and thermal sampling never share a generator state. + """ + if not isinstance(base_seed, (int, np.integer)) or int(base_seed) < 0: + raise ValueError(f"base_seed must be a non-negative integer, got {base_seed}") + if not np.isfinite(p) or not (0.0 < p < 0.5): + raise ValueError(f"p must be in (0, 0.5), got {p}") + if not isinstance(L, (int, np.integer)) or int(L) < 4: + raise ValueError(f"L must be an integer >= 4, got {L}") + if not isinstance(replicate, (int, np.integer)) or int(replicate) < 0: + raise ValueError(f"replicate must be a non-negative integer, got {replicate}") + if policy not in ("independent", "common_random_numbers"): + raise ValueError(f"unknown seed policy: {policy}") + + p_lo, p_hi = _float64_words(p) + policy_tag = 0x494E4450 if policy == "independent" else 0x43524E00 + entropy = [int(base_seed), int(L), int(replicate), policy_tag] + if policy == "independent": + entropy.extend((p_lo, p_hi)) + + root = np.random.SeedSequence(entropy) + bond_ss, thermal_ss = root.spawn(2) + bond_seed = int(bond_ss.generate_state(1, dtype=np.uint64)[0]) + thermal_seed = int(thermal_ss.generate_state(1, dtype=np.uint64)[0]) + return StreamSeeds( + bond_seed=bond_seed, + thermal_seed=thermal_seed, + policy=policy, + p=float(p), + L=int(L), + replicate=int(replicate), + ) + + +def nishimori_scan_seeded( + p: float, + L: int, + *, + n_disorder: int, + n_records: int, + burn_in: int, + base_seed: int, + seed_policy: SeedPolicy = "independent", + n_skip: int = 1, + sweeps_per_step: int = 2, + aligned_start: bool = True, +) -> DisorderResult: + """Run one Nishimori scan point using an explicit random-stream policy.""" + if not isinstance(n_disorder, (int, np.integer)) or n_disorder < 1: + raise ValueError(f"n_disorder must be an integer >= 1, got {n_disorder}") + beta = nishimori_beta(p) + n = L * L + per_real_absm: list[float] = [] + per_real_m2: list[float] = [] + per_real_m4: list[float] = [] + cfracs: list[float] = [] + + for replicate in range(n_disorder): + seeds = derive_stream_seeds( + base_seed=base_seed, + p=p, + L=L, + replicate=replicate, + policy=seed_policy, + ) + bonds = sample_bonds(p, L, seed=seeds.bond_seed) + configs, cf = rbim_wolff_sample( + bonds, + beta, + n_records=n_records, + burn_in=burn_in, + seed=seeds.thermal_seed, + n_skip=n_skip, + sweeps_per_step=sweeps_per_step, + aligned_start=aligned_start, + ) + m = configs.reshape(configs.shape[0], -1).sum(axis=1) / n + per_real_absm.append(float(np.mean(np.abs(m)))) + per_real_m2.append(float(np.mean(m * m))) + per_real_m4.append(float(np.mean(m**4))) + cfracs.append(cf) + + absm_arr = np.asarray(per_real_absm, dtype=np.float64) + m2 = float(np.mean(per_real_m2)) + m4 = float(np.mean(per_real_m4)) + binder = 1.0 - m4 / (3.0 * m2 * m2) if m2 > 0.0 else float("nan") + sem = float(np.std(absm_arr, ddof=1) / math.sqrt(len(absm_arr))) if len(absm_arr) > 1 else 0.0 + return DisorderResult( + p=float(p), + beta=float(beta), + L=int(L), + abs_m=float(np.mean(absm_arr)), + abs_m_err=sem, + m2=m2, + m4=m4, + binder=binder, + mean_cluster_frac=float(np.mean(cfracs)), + n_disorder=int(n_disorder), + n_records=int(n_records), + ) + + +def nishimori_scan_grid( + ps: tuple[float, ...] | list[float], + L: int, + **kwargs, +) -> list[DisorderResult]: + """Run an ordered p-grid with the same explicit seed policy at every point.""" + ps_arr = np.asarray(ps, dtype=np.float64) + if ps_arr.ndim != 1 or len(ps_arr) < 2: + raise ValueError("ps must be a one-dimensional sequence with at least two points") + if not np.all(np.isfinite(ps_arr)) or np.any((ps_arr <= 0.0) | (ps_arr >= 0.5)): + raise ValueError("all p values must be finite and in (0, 0.5)") + if np.any(np.diff(ps_arr) <= 0.0): + raise ValueError("ps must be strictly increasing") + return [nishimori_scan_seeded(float(p), L, **kwargs) for p in ps_arr] diff --git a/tests/test_ising2d.py b/tests/test_ising2d.py index 65e49b4..3493abd 100644 --- a/tests/test_ising2d.py +++ b/tests/test_ising2d.py @@ -1,4 +1,4 @@ -"""Tests fuer den 2D-Ising-Sampler + Majority-Rule-Blocking (Phase-3b).""" +"""Tests for the 2D-Ising sampler and majority-rule blocking.""" from __future__ import annotations @@ -9,14 +9,12 @@ def test_kc_tc_constants() -> None: - """K_c = ln(1+sqrt(2))/2, T_c = 1/K_c (Onsager).""" assert pytest.approx(0.5 * np.log(1 + np.sqrt(2))) == i2.KC_2D assert pytest.approx(1.0 / i2.KC_2D) == i2.TC_2D assert pytest.approx(2.2691853, abs=1e-5) == i2.TC_2D def test_energy_per_spin_ground_state() -> None: - """Voll ausgerichtetes Gitter -> E/N = -2 (alle 2N Bonds = +1).""" s = np.ones((8, 8)) assert float(i2.energy_per_spin(s)) == pytest.approx(-2.0) assert float(i2.energy_per_spin(-s)) == pytest.approx(-2.0) @@ -29,14 +27,16 @@ def test_magnetization_per_spin() -> None: def test_metropolis_vs_exact_L4() -> None: - """Sampler reproduziert die exakte L=4-Enumeration der Energie (Orakel).""" + """Sampler reproduces the exact L=4 energy enumeration.""" L, n = 4, 16 states = np.arange(1 << n, dtype=np.int64) bits = ((states[:, None] >> np.arange(n)[None, :]) & 1).astype(np.int8) s = (1 - 2 * bits).reshape(-1, L, L).astype(np.float64) e_all = i2.energy_per_spin(s) for K in (0.25, i2.KC_2D, 0.55): - w = np.exp(-K * n * e_all) + logw = -K * n * e_all + logw -= logw.max() + w = np.exp(logw) w /= w.sum() e_exact = float((w * e_all).sum()) ch = i2.checkerboard_metropolis(K, L, n_sweeps=20000, burn_in=3000, seed=7, record_every=2) @@ -45,7 +45,6 @@ def test_metropolis_vs_exact_L4() -> None: def test_metropolis_reproducible() -> None: - """Gleicher Seed -> bit-identische Trajektorie; anderer Seed -> verschieden.""" kw = dict(n_sweeps=500, burn_in=100, record_every=1) a = i2.checkerboard_metropolis(0.4, 16, seed=5, **kw) b = i2.checkerboard_metropolis(0.4, 16, seed=5, **kw) @@ -55,7 +54,6 @@ def test_metropolis_reproducible() -> None: def test_metropolis_ordered_above_kc() -> None: - """Bei K > K_c (kalt) -> grosse Magnetisierung; bei K < K_c -> klein.""" cold = i2.checkerboard_metropolis(0.6, 16, n_sweeps=3000, burn_in=1500, seed=1, record_every=5) hot = i2.checkerboard_metropolis(0.25, 16, n_sweeps=3000, burn_in=1500, seed=1, record_every=5) m_cold = float(np.abs(i2.magnetization_per_spin(cold.configs)).mean()) @@ -66,7 +64,6 @@ def test_metropolis_ordered_above_kc() -> None: def test_majority_block_shape_and_rule() -> None: - """b=2-Blocking halbiert jede Dimension; Mehrheit korrekt.""" s = np.ones((8, 8), dtype=np.int64) b = i2.majority_block_b2(s) assert b.shape == (4, 4) @@ -76,10 +73,10 @@ def test_majority_block_shape_and_rule() -> None: def test_majority_tie_break_unbiased() -> None: - """Gleichstand-Bloecke (Summe 0) werden im Mittel symmetrisch +/- gebrochen.""" + """Tie selector is balanced in aggregate across block/config hashes.""" L = 32 ii, jj = np.indices((L, L)) - s = np.where((ii + jj) % 2 == 0, 1, -1) # jeder 2x2-Block summiert zu 0 + s = np.where((ii + jj) % 2 == 0, 1, -1) plus = total = 0 for ci in range(200): b = i2.majority_block_b2(s, config_index=ci, seed=42) @@ -97,6 +94,39 @@ def test_majority_tie_break_reproducible() -> None: assert not np.array_equal(a, c) +def test_majority_block_exact_z2_equivariance() -> None: + """Global spin flip must commute with blocking, including every tie block.""" + rng = np.random.default_rng(2026) + + # Explicit all-tie checkerboard stresses the formerly defective path. + ii, jj = np.indices((16, 16)) + tied = np.where((ii + jj) % 2 == 0, 1, -1).astype(np.int8) + for seed in (0, 1, 7, 2**63 + 5): + for config_index in (0, 3, 99): + b = i2.majority_block_b2(tied, config_index=config_index, seed=seed) + b_flip = i2.majority_block_b2(-tied, config_index=config_index, seed=seed) + assert np.array_equal(b_flip, -b) + + # Also exercise mixed majority/tie configurations. + for config_index in range(20): + s = rng.choice(np.array([-1, 1], dtype=np.int8), size=(16, 16)) + b = i2.majority_block_b2(s, config_index=config_index, seed=42) + b_flip = i2.majority_block_b2(-s, config_index=config_index, seed=42) + assert np.array_equal(b_flip, -b) + + +def test_majority_tie_output_is_an_input_spin() -> None: + """On ties the selected coarse spin must come from the corresponding 2x2 block.""" + rng = np.random.default_rng(11) + s = rng.choice(np.array([-1, 1], dtype=np.int8), size=(12, 12)) + b = i2.majority_block_b2(s, config_index=5, seed=17) + lb = s.shape[0] // 2 + blocks = s.reshape(lb, 2, lb, 2).transpose(0, 2, 1, 3).reshape(lb, lb, 4) + sums = blocks.sum(axis=2) + for r, c in zip(*np.where(sums == 0), strict=True): + assert b[r, c] in blocks[r, c] + + @pytest.mark.parametrize( "bad", [ @@ -109,8 +139,182 @@ def test_majority_tie_break_reproducible() -> None: lambda: i2.exact_energy_per_spin_2x2(-1.0), lambda: i2.majority_block_b2(np.ones((3, 3))), lambda: i2.majority_block_b2(np.ones((4, 6))), + lambda: i2.majority_block_b2(np.array([[1, 0], [-1, 1]])), ], ) def test_edge_inputs_raise(bad) -> None: with pytest.raises((ValueError, TypeError)): bad() + + +def _all_2x2_configs() -> np.ndarray: + """Alle 2^4 = 16 moeglichen 2x2-Bloecke in {+1,-1}.""" + states = np.arange(16, dtype=np.int64) + bits = ((states[:, None] >> np.arange(4)[None, :]) & 1).astype(np.int8) + return (1 - 2 * bits).reshape(16, 2, 2) + + +def test_majority_block_z2_equivariance_exhaustive_2x2() -> None: + """Z2-Aequivarianz erschoepfend auf der Block-Ebene. + + ``majority_block_b2`` wirkt blockweise und unabhaengig je 2x2-Block; die + einzige Kopplung an die Gittergroesse ist der Hash ueber (block_row, + block_col). Daher ist ein Sweep ueber ALLE 16 moeglichen Blockinhalte, + gekreuzt mit allen vier erreichbaren Selektorwerten, ein vollstaendiger + Nachweis der Eigenschaft auf Blockebene -- keine Stichprobe. + + Der Tie-Pfad ist dabei nicht Beiwerk: 6 der 16 Blockinhalte (C(4,2)) sind + 2+2-Ties, also genau der Pfad, auf dem die alte Hash-Bit-Regel + ``B(-s) != -B(s)`` lieferte. + """ + blocks = _all_2x2_configs() + seen_selectors: set[int] = set() + n_tie_checked = 0 + n_checked = 0 + + for seed in (0, 1, 17, 2**63 + 5): + for config_index in range(64): + selector = int( + i2._splitmix64_grid( + config_index, + np.zeros((1, 1), dtype=np.int64), + np.zeros((1, 1), dtype=np.int64), + seed, + )[0, 0] + & np.uint64(3) + ) + seen_selectors.add(selector) + for s in blocks: + b = i2.majority_block_b2(s, config_index=config_index, seed=seed) + b_flip = i2.majority_block_b2(-s, config_index=config_index, seed=seed) + assert np.array_equal(b_flip, -b), (s, config_index, seed) + n_checked += 1 + if s.sum() == 0: + n_tie_checked += 1 + + # Nicht-Vakuitaet: alle vier Selektor-Slots und alle 6 Tie-Muster wurden + # wirklich durchlaufen -- sonst wuerde ein Sweep, der den Tie-Pfad nie + # trifft, ebenfalls bestehen und nichts beweisen. + assert seen_selectors == {0, 1, 2, 3}, seen_selectors + assert n_checked == 16 * 64 * 4 + assert n_tie_checked == 6 * 64 * 4 + + +def test_majority_block_z2_equivariance_randomized_large() -> None: + """Dieselbe Invariante auf grossen Gittern, deterministisch geseedet. + + Deckt die Hash-Indizierung ueber (block_row, block_col) ab, die der + erschoepfende 1x1-Block-Sweep nicht beruehrt. Fester ``Generator``, damit + ein Fehlschlag reproduzierbar ist. + """ + rng = np.random.default_rng(20260910) + tie_blocks_seen = 0 + for config_index in range(40): + for L in (8, 16, 32): + s = rng.choice(np.array([-1, 1], dtype=np.int8), size=(L, L)) + b = i2.majority_block_b2(s, config_index=config_index, seed=4711) + b_flip = i2.majority_block_b2(-s, config_index=config_index, seed=4711) + assert np.array_equal(b_flip, -b) + lb = L // 2 + tie_blocks_seen += int((s.reshape(lb, 2, lb, 2).sum(axis=(1, 3)) == 0).sum()) + # Zufallskonfigurationen erzeugen reichlich Ties; ohne sie waere der Sweep + # blind fuer genau den reparierten Pfad. + assert tie_blocks_seen > 1000, tie_blocks_seen + + +def test_majority_block_non_tie_blocks_ignore_the_hash() -> None: + """Positiv-Kontrolle: der Tie-Pfad darf Nicht-Tie-Bloecke NICHT anfassen. + + Eine Tie-Regel, die alles ueberschreibt (oder eine Implementierung, die + jeden Block als Tie behandelt), wuerde jeden Aequivarianz-Test bestehen und + waere trotzdem falsch: Z2-Aequivarianz allein ist von ``B(s) = s[0,0]`` + ebenfalls erfuellt. Hier wird daher festgehalten, dass Bloecke mit echter + Mehrheit unabhaengig von config_index/seed das Vorzeichen der Blocksumme + liefern. + """ + blocks = _all_2x2_configs() + non_tie = np.array([blk for blk in blocks if blk.sum() != 0]) + assert non_tie.shape[0] == 10 # 16 - 6 Ties + + for blk in non_tie: + expected = np.sign(blk.sum()) + for config_index in (0, 5, 12345): + for seed in (0, 99, 2**63 + 5): + out = i2.majority_block_b2(blk, config_index=config_index, seed=seed) + assert out.shape == (1, 1) + assert out[0, 0] == expected, (blk, config_index, seed) + + # Und auf einem grossen Gitter ohne jeden Tie: Ergebnis == reines Mehrheits- + # Vorzeichen, hash-unabhaengig. + rng = np.random.default_rng(7) + # Jeder 2x2-Block einheitlich +1 oder -1 -> Blocksumme immer +/-4, nie ein Tie, + # aber das Gitter ist echt variiert (kein triviales Eins-Gitter). + coarse = rng.choice(np.array([-1, 1], dtype=np.int8), size=(8, 8)) + s = np.kron(coarse, np.ones((2, 2), dtype=np.int8)) + block_sum = s.reshape(8, 2, 8, 2).sum(axis=(1, 3)) + assert np.array_equal(np.sign(block_sum).astype(np.int8), coarse) + assert np.all(block_sum != 0) + a = i2.majority_block_b2(s, config_index=1, seed=1) + b = i2.majority_block_b2(s, config_index=2, seed=2) + assert np.array_equal(a, b) + assert np.array_equal(a, np.sign(block_sum).astype(np.int8)) + + +@pytest.mark.parametrize( + "bad, warum", + [ + (np.array([[1.5, -1.5], [1.5, -1.5]]), "1.5/-1.5 wuerde ein Cast zu 1/-1 abschneiden"), + (np.array([[0.5, 0.5], [0.5, 0.5]]), "0.5 -> 0"), + (np.array([[0.0, 0.0], [0.0, 0.0]]), "0 ist kein Spin"), + (np.array([[2.0, -2.0], [2.0, -2.0]]), "Betrag != 1"), + (np.array([[0.9999999999, -1.0], [1.0, -1.0]]), "knapp neben +1 -> 0"), + (np.array([[np.nan, 1.0], [1.0, -1.0]]), "NaN"), + (np.array([[np.inf, 1.0], [1.0, -1.0]]), "inf"), + ], +) +def test_majority_block_rejects_non_spin_values_before_casting(bad, warum) -> None: + """Die ±1-Pruefung muss VOR dem int-Cast greifen. + + Laeuft der Cast zuerst, kann die Pruefung nicht mehr sehen, wogegen sie + schuetzt: ``1.5`` und ``-1.5`` werden zu ``1`` und ``-1`` abgeschnitten und + danach als gueltige Spins akzeptiert -- die Funktion blockt dann still + veraenderte Daten, statt den dokumentierten ±1-Fehler zu werfen. + + Breit gefangen und der TYP geprueft: stirbt der Aufruf an einer anderen + Ausnahme (etwa einem Cast-Fehler), waere der Test sonst nicht einzuordnen. + """ + with pytest.raises(Exception) as exc: + i2.majority_block_b2(bad) + assert isinstance(exc.value, ValueError), ( + f"erwartet ValueError ({warum}), kam {type(exc.value).__name__}: {exc.value}" + ) + assert "+/-1" in str(exc.value), f"erwartet die dokumentierte ±1-Meldung, kam: {exc.value}" + + +def test_majority_block_still_accepts_valid_spins_in_any_container() -> None: + """Positiv-Kontrolle: gueltige ±1-Daten muessen weiter akzeptiert werden. + + Ein Validator, der jede float-Eingabe ablehnt, bestuende jeden Negativtest + und waere trotzdem falsch. Float-Spins (1.0/-1.0), int8-Arrays und rohe + Listen sind gueltige Eingaben und muessen dasselbe Ergebnis liefern. + """ + ref = np.array([[1, -1, -1, 1], [1, -1, 1, 1], [-1, -1, 1, -1], [1, 1, -1, -1]]) + erwartet = i2.majority_block_b2(ref, config_index=3, seed=17) + for variante in ( + ref.astype(np.float64), + ref.astype(np.float32), + ref.astype(np.int8), + ref.tolist(), + ): + wie = getattr(variante, "dtype", type(variante).__name__) + # Eine Ablehnung ist hier ein FEHLSCHLAG, keine Ausnahme, die den Test + # abstuerzen laesst: ein abgestuerzter Test ist nicht einzuordnen und + # taugt nicht als Beleg. Darum in eine Zusicherung uebersetzen. + try: + out = i2.majority_block_b2(variante, config_index=3, seed=17) + except Exception as exc: + raise AssertionError( + f"gueltige +/-1-Daten als {wie} wurden abgelehnt: {type(exc).__name__}: {exc}" + ) from exc + assert np.array_equal(out, erwartet), f"abweichendes Ergebnis fuer {wie}" + assert out.dtype == np.int8, f"erwartet int8 fuer {wie}, kam {out.dtype}" diff --git a/tests/test_mcrg_multirg.py b/tests/test_mcrg_multirg.py index b42a9a5..4e6e9b0 100644 --- a/tests/test_mcrg_multirg.py +++ b/tests/test_mcrg_multirg.py @@ -148,22 +148,36 @@ def test_block_size_per_iter_exposed() -> None: def test_fix3_central_values_unchanged() -> None: - """Codex-Fix 3 aendert NUR Fehlerbalken, NICHT die y_t/y_h-Zentralwerte. - - Der Punktschaetzer nutzt weiter das gemeinsame Level-0-Fenster `keep`; nur - die Jackknife-Partition wird pro Iteration gewaehlt. Verankert die EXAKTEN - G27/G28-Zentralwerte (validate_multirg_2d, seed=0, L=32, n_op=2, burn_in=400) - gegen einen Regress. Die Baseline ist tool-gemessen VOR Fix 3 und wird nach - Fix 3 byte-identisch reproduziert (nur die Fehlerbalken duerfen sich aendern). + """Reproduzierbarkeits-Anker fuer die G27/G28-Zentralwerte. + + Verankert die EXAKTEN Zentralwerte (validate_multirg_2d, seed=0, L=32, + n_op=2, burn_in=400) gegen einen unbeabsichtigten Regress. + + REFERENZ-HERKUNFT. Die urspruengliche Referenz war unter dem ALTEN, nicht + Z2-aequivarianten Majority-Tie-Break gemessen. Seit + `ising2d.majority_block_b2` bei einem 2+2-Tie einen der vier Original-Spins + waehlt, gilt `B(-s) == -B(s)` exakt statt nur im Mittel -- die + Blocking-Abbildung ist bewusst eine andere, und y_t/y_h als Funktionale der + geblockten Konfigurationen verschieben sich zwangslaeufig. Die alte + Erwartung haette ab dem Fix behauptet, eine absichtliche Aenderung der + RG-Abbildung sei wirkungslos. + + Die neue Referenz ist nicht deshalb richtig, weil sie herauskam, sondern + weil (1) die erzeugende Abbildung jetzt exakt Z2-aequivariant ist + (erschoepfend belegt in test_ising2d.py), (2) das externe Onsager-Orakel + weiter getroffen wird (G22/G27/G28 unveraendert PASS) und (3) die + eigentliche Fix-3-Behauptung -- Jackknife-Blockgroesse bewegt nur + Fehlerbalken -- snapshot-frei in + test_jackknife_block_size_moves_only_error_bars geprueft wird. """ v = mcrg_multirg.validate_multirg_2d( L=32, n_op_even=2, n_op_odd=2, n_levels=3, n_records=3000, burn_in=400, seed=0 ) np.testing.assert_allclose( - v.multirg.y_t_per_iter, [0.93001085, 0.99566981, 1.0219697], rtol=0, atol=1e-7 + v.multirg.y_t_per_iter, [0.93220357, 1.01267256, 1.00535087], rtol=0, atol=1e-7 ) np.testing.assert_allclose( - v.multirg_odd.y_h_per_iter, [1.88098809, 1.87282836, 1.87101431], rtol=0, atol=1e-7 + v.multirg_odd.y_h_per_iter, [1.88074043, 1.87301563, 1.87250009], rtol=0, atol=1e-7 ) # Alle Fehlerbalken endlich + nicht-negativ. assert np.all(np.isfinite(v.multirg.y_t_err_per_iter)) and np.all( @@ -179,3 +193,48 @@ def test_explicit_block_size_respected() -> None: chain = _wolff_chain(L=32, n=3000, seed=0) res = mcrg_multirg.multi_rg_y_t(chain, n_op=2, n_levels=3, block_size=10) assert np.all(np.asarray(res.block_size_per_iter) == 10) + + +def test_jackknife_block_size_moves_only_error_bars() -> None: + """Fix-3-Kernbehauptung, snapshot-frei: die Jackknife-Blockgroesse darf NUR + die Fehlerbalken bewegen, nie die Zentralwerte. + + Warum das die richtige Formulierung ist: der Zentralwert wird in + ``multi_rg_y_t``/``multi_rg_y_h`` auf dem vollen getrimmten Fenster berechnet, + BEVOR und unabhaengig davon die Jackknife-Schleife die Blockgroesse benutzt + (mcrg_multirg: ``y_iter[p] = _y_t_even_from_levels(S_lo, S_hi, ...)`` vor der + ``for j in range(nb)``-Schleife). Ein Regressionstest gegen fest verdrahtete + Zahlen prueft das nur mittelbar und bricht bei jeder absichtlichen + Physik-Aenderung; diese Fassung prueft die Invariante direkt und ist + unabhaengig von der Tie-Break-Konvention. + + ``block_size=50`` teilt ``n_records=800`` ohne Rest (800//50 = 16 Bloecke), + also ist das getrimmte Fenster ``keep`` in beiden Laeufen identisch — die + Zentralwerte MUESSEN daher bit-identisch sein. + """ + n_records, global_block = 800, 50 + chain = _wolff_chain(16, n_records, seed=0) + + per_iter_t = mcrg_multirg.multi_rg_y_t(chain, n_op=2, n_levels=3, seed=0) + global_t = mcrg_multirg.multi_rg_y_t(chain, n_op=2, n_levels=3, seed=0, block_size=global_block) + per_iter_h = mcrg_multirg.multi_rg_y_h(chain, n_op=2, n_levels=3, seed=0) + global_h = mcrg_multirg.multi_rg_y_h(chain, n_op=2, n_levels=3, seed=0, block_size=global_block) + + # (b) Zentralwerte bit-identisch -- die eigentliche Invariante. + assert np.array_equal(per_iter_t.y_t_per_iter, global_t.y_t_per_iter) + assert np.array_equal(per_iter_h.y_h_per_iter, global_h.y_h_per_iter) + + # Nicht-Vakuitaet: die beiden Politiken muessen sich ueberhaupt unterscheiden, + # sonst wuerde die Gleichheit oben nichts bedeuten. + assert not np.array_equal(per_iter_t.block_size_per_iter, global_t.block_size_per_iter) + assert not np.array_equal(per_iter_h.block_size_per_iter, global_h.block_size_per_iter) + assert not np.allclose(per_iter_t.y_t_err_per_iter, global_t.y_t_err_per_iter) + assert not np.allclose(per_iter_h.y_h_err_per_iter, global_h.y_h_err_per_iter) + + # (a)/(c) per-iter Blockgroessen werden gemeldet und sind gueltig. + assert per_iter_t.block_size_per_iter.shape[0] == per_iter_t.n_iters + assert per_iter_h.block_size_per_iter.shape[0] == per_iter_h.n_iters + assert np.all(per_iter_t.block_size_per_iter >= 1) + assert np.all(per_iter_h.block_size_per_iter >= 1) + assert np.all(np.isfinite(per_iter_t.y_t_err_per_iter)) + assert np.all(np.isfinite(per_iter_h.y_h_err_per_iter)) diff --git a/tests/test_rbim_scan.py b/tests/test_rbim_scan.py new file mode 100644 index 0000000..6c47816 --- /dev/null +++ b/tests/test_rbim_scan.py @@ -0,0 +1,91 @@ +"""Tests for explicit RBIM random-stream policies.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from adaptiverg_qec.rbim_scan import derive_stream_seeds, nishimori_scan_grid + + +def test_independent_policy_is_reproducible_and_p_specific() -> None: + a = derive_stream_seeds(base_seed=2026, p=0.10, L=8, replicate=3) + b = derive_stream_seeds(base_seed=2026, p=0.10, L=8, replicate=3) + c = derive_stream_seeds(base_seed=2026, p=0.11, L=8, replicate=3) + assert a == b + assert (a.bond_seed, a.thermal_seed) != (c.bond_seed, c.thermal_seed) + assert a.bond_seed != a.thermal_seed + + +def test_common_random_numbers_policy_is_explicitly_p_shared() -> None: + a = derive_stream_seeds( + base_seed=2026, + p=0.10, + L=8, + replicate=3, + policy="common_random_numbers", + ) + b = derive_stream_seeds( + base_seed=2026, + p=0.11, + L=8, + replicate=3, + policy="common_random_numbers", + ) + assert (a.bond_seed, a.thermal_seed) == (b.bond_seed, b.thermal_seed) + + +def test_stream_identity_changes_with_lattice_and_replicate() -> None: + base = derive_stream_seeds(base_seed=17, p=0.1, L=8, replicate=0) + other_l = derive_stream_seeds(base_seed=17, p=0.1, L=10, replicate=0) + other_rep = derive_stream_seeds(base_seed=17, p=0.1, L=8, replicate=1) + assert base.bond_seed != other_l.bond_seed + assert base.bond_seed != other_rep.bond_seed + + +def test_many_child_streams_have_no_deterministic_collisions() -> None: + seeds = { + derive_stream_seeds(base_seed=99, p=0.1, L=8, replicate=r).bond_seed for r in range(1000) + } + assert len(seeds) == 1000 + + +@pytest.mark.parametrize( + "kwargs", + [ + {"base_seed": -1, "p": 0.1, "L": 8, "replicate": 0}, + {"base_seed": 1, "p": 0.0, "L": 8, "replicate": 0}, + {"base_seed": 1, "p": np.nan, "L": 8, "replicate": 0}, + {"base_seed": 1, "p": 0.1, "L": 2, "replicate": 0}, + {"base_seed": 1, "p": 0.1, "L": 8, "replicate": -1}, + {"base_seed": 1, "p": 0.1, "L": 8, "replicate": 0, "policy": "implicit"}, + ], +) +def test_invalid_seed_identity_rejected(kwargs) -> None: + with pytest.raises(ValueError): + derive_stream_seeds(**kwargs) + + +def test_seeded_scan_grid_reproducible_smoke() -> None: + kwargs = dict( + n_disorder=2, + n_records=8, + burn_in=8, + base_seed=123, + n_skip=1, + sweeps_per_step=1, + aligned_start=True, + ) + a = nishimori_scan_grid([0.08, 0.12], 4, **kwargs) + b = nishimori_scan_grid([0.08, 0.12], 4, **kwargs) + assert [(r.abs_m, r.binder) for r in a] == [(r.abs_m, r.binder) for r in b] + + +def test_scan_grid_requires_strictly_increasing_valid_ps() -> None: + kwargs = dict(n_disorder=1, n_records=2, burn_in=0, base_seed=1) + with pytest.raises(ValueError): + nishimori_scan_grid([0.1, 0.1], 4, **kwargs) + with pytest.raises(ValueError): + nishimori_scan_grid([0.2, 0.1], 4, **kwargs) + with pytest.raises(ValueError): + nishimori_scan_grid([0.0, 0.1], 4, **kwargs)