diff --git a/AGENTS.md b/AGENTS.md index 07c6010..9722257 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -48,6 +48,41 @@ At the start of a new task: - `pepsy.sampling`: `MpsSampler`, `VecSampler`, `PepsBpSampler`, and result dataclasses. - `pepsy._internal`: private formatting and small utility helpers only. +## Upstream Tensor-Network Substrate + +- Treat `autoray`, `cotengra`, `cotengrust`, and `quimb` as first-class + dependencies and prefer their APIs over local reimplementations wherever they + cover the needed tensor-network, contraction-planning, and backend-dispatch + behavior. +- Use `quimb` tensor-network objects and methods for Tensor/TensorNetwork, + MPS/MPO/PEPS/PEPO construction, contraction, boundary, gate, geometry, + sampling, and optimization workflows when practical. Keep Pepsy wrappers + thin and focused on Pepsy conventions, compatibility, and stable public APIs. +- Use `cotengra` for contraction path/tree optimization, hyper-optimization, + slicing, subtree reconfiguration, compressed contraction planning, and + reusable optimizer objects. Forward cotengra-compatible options rather than + inventing parallel path-search configuration inside Pepsy. +- Use `cotengrust` through cotengra's public acceleration paths for Rust-backed + greedy, random-greedy, optimal, and subtree-reconfiguration path-finding. + Keep cotengra as the main optimizer interface unless a task specifically + needs a lower-level cotengrust primitive. +- Use `autoray` for backend inference and backend-agnostic array operations, + including creation, conversion, `einsum`/`tensordot`, reshaping, conjugation, + linalg, and registered backend-specific gradient/linalg fixes. Preserve + compatibility with NumPy, Torch, JAX, CuPy, and other autoray-compatible + arrays when optional dependencies are installed. +- Preserve compatibility with the supported dependency ranges in + `pyproject.toml`; handle version-specific upstream behavior with feature + detection and focused regression tests. +- Do not vendor or copy upstream internals. If Pepsy needs a workaround for an + upstream behavior, isolate it behind a small adapter, document why it exists, + and test it against the closest public upstream API. +- `symmray` is the planned symmetric/block-sparse integration path. Design new + tensor, gate, boundary, and optimizer code so Symmray-style arrays can flow + through quimb/autoray paths where possible, keep `symmray` optional unless it + becomes a declared dependency, and use `pytest.importorskip("symmray")` for + Symmray-specific tests. + ## Core Workflows - Boundary contraction flow: `build_bra_ket(ket, bra?) -> BdyMPS(...) -> contract_boundary(...)`. @@ -115,6 +150,14 @@ NUMBA_CACHE_DIR=/tmp/numba_cache MPLCONFIGDIR=/tmp/mplconfig PYTHONPYCACHEPREFIX ## Examples and Notebooks - Keep examples lightweight and deterministic where practical. +- Write examples compositionally: expose small setup/build/evolve/measure steps + that compose public APIs instead of copying internal implementation details. +- For Pepsy examples, prefer Pepsy public functionality wherever possible; drop + to `quimb`, `autoray`, or `cotengra` only when the example specifically needs + lower-level tensor-network, backend, or contraction-planning control. +- For Gaugy examples, prefer Gaugy public functionality first, then Pepsy public + functionality, then `quimb`; use `autoray`/`cotengra` only for explicit + backend or contraction-optimizer concerns. - Do not modify generated notebook outputs unless explicitly requested. - When changing public examples, verify that imports use current namespaces and public API names. - If a notebook or helper still shows an old flat import path, migrate it deliberately instead of copying that pattern into new code. diff --git a/PLAN.md b/PLAN.md index c936059..e860184 100644 --- a/PLAN.md +++ b/PLAN.md @@ -1,10 +1,10 @@ -# PLAN.md — Belief Propagation, quimb integration, and Symmetric Tensor Networks +# PLAN.md — BP, quimb, Symmetric TN, and Differentiable Truncations Status: draft / living document -Last updated: 2026-06-24 +Last updated: 2026-06-28 Owners: pepsy maintainers + coding agents -This document plans three related workstreams to extend `pepsy`: +This document plans four related workstreams to extend `pepsy`: 1. **BP** — Belief Propagation (BP) contraction and gauging, plus the *loop series / loop cluster expansion* corrections that turn BP from an @@ -15,10 +15,15 @@ This document plans three related workstreams to extend `pepsy`: 3. **Symmetric TN** — block-sparse abelian-symmetric and fermionic tensors via [`symmray`](https://github.com/jcmgray/symmray), so PEPS/MPS workflows can exploit U(1)/Z2 symmetries and fermionic statistics. +4. **Differentiable truncations** — implement stable, efficient AD rules for + truncated SVD/eigendecomposition/QR-style decompositions, with the + Francuz-Schuch-Vanhecke missing truncated-spectrum term and a clear split + between "full SVD then truncate" and true low-rank `svds` paths. -These three are deliberately grouped: `quimb` is the substrate, BP is the new -algorithmic layer, and `symmray` provides the symmetric array backend that both -can run on top of (all three share the same author/ecosystem, so they compose). +These are deliberately grouped: `quimb` is the substrate, BP is the new +algorithmic layer, `symmray` provides the symmetric array backend that both can +run on top of, and differentiable truncations make gradient-based PEPS/MPS/MPO +optimization reliable across all of those paths. Companion folders: - `history/` — session-to-session hand-off notes for agents (read on start, @@ -67,6 +72,31 @@ Companion folders: `ham_tfim_from_edges`, `ham_heisenberg_from_edges`, `ham_fermi_hubbard_*`. Reference: Gao et al., *Phys. Rev. Research* 7, 023193 (2025). +### Stable AD for truncated decompositions +- **Stable and efficient differentiation of tensor network algorithms** — + Francuz, Schuch, Vanhecke, *Phys. Rev. Research* 7, 013237 (2025). + DOI: https://doi.org/10.1103/PhysRevResearch.7.013237 · + arXiv: https://arxiv.org/abs/2311.11894 + - Core message for `pepsy`: differentiating a true truncated SVD is not the + same as differentiating a full SVD and then slicing. The true truncated + pullback contains an extra term from the discarded spectrum; older AD + treatments effectively assumed that discarded spectrum was zero. + - If the forward pass already computes a full SVD and only then truncates, + ordinary full-SVD AD plus separate slicing is a valid baseline. The paper's + efficient rule matters when the forward pass only computes the kept + subspace (`svds`/iterative/randomized low-rank SVD), or when memory and + large matrix dimensions make full SVD undesirable. + - Gauge-dependent singular-vector losses are ill-conditioned near + degeneracy. Tests must use gauge/phase-invariant objectives, or explicitly + remove the gauge component of the cotangents. +- **TensorKit / MatrixAlgebraKit reference implementation** + - Current TensorKit moved factorization pullbacks into + `src/factorizations/pullbacks.jl` and delegates dense block pullbacks to + MatrixAlgebraKit, rather than keeping the old + `ext/TensorKitChainRulesCoreExt/factorizations.jl` entry point. + - Useful design cue for `pepsy`: keep tensor wrappers thin and put the AD + rules at the array/decomposition layer. + --- ## 1. Where this lands in pepsy @@ -77,7 +107,7 @@ Current relevant architecture (do not break these public entry points): - `pepsy.tensors`: `OneDMap`, product/identity constructors, contraction optimizers, observables, `tn_norm`, `tn_fidelity`. - `pepsy.backends`: backend inference/conversion + torch/JAX/CuPy linalg - registration. + registration, including custom linalg AD rules. - `pepsy.operators`, `pepsy.solvers`, `pepsy.optimizers`, `pepsy.sampling`. BP is conceptually a *peer* of boundary-MPS contraction: another way to @@ -96,6 +126,28 @@ Symmetric tensors are a *backend/array* concern, not an algorithm: validate symmetry tags, and bridge pepsy lattice tags (`X{i}`, `Y{j}`, `I…`, `k…`/`b…` legs) to symmray index `dual`/charge conventions. +Differentiable truncations are also a *backend/decomposition* concern: +- The existing `pepsy.backends.linalg_torch.SVD` registers a stabilized full + SVD via `autoray` as `torch.linalg.svd`. The preferred public hook is + `pepsy.tensors.reg_rel_svd_torch()`; `reg_complex_svd_torch()` and + `register_torch_linalg(mode="complex")` are compatibility paths to the same + relative-regularized rule. It uses Townsend's rectangular full-SVD backward + structure, Lorentzian relative broadening for singular-value denominators, + the complex SVD gauge correction, and a CPU SciPy `gesvd` forward fallback + when `torch.linalg.svd` fails. This helps full-SVD AD and robustness, but by + itself does not implement the truncated-SVD pullback from the paper. +- Installed `quimb` routes tensor splits through `quimb.tensor.decomp.array_split`. + The default `svd_truncated` currently computes `xp.linalg.svd(x)` and then + slices/trims the factors in `_trim_and_renorm_svd_result`. Thus: + - For dense full-SVD paths, `pepsy` can continue to register full-SVD AD and + treat truncation as an independent slice. + - For efficient low-rank paths, `pepsy` needs a new split driver or adapter + that exposes a custom VJP/JVP for `svds`-style truncated factors. +- Primary call sites that should benefit without scattered rewrites: + `Tensor.split`/`tensor_network_gate_inds`, `compress_between`, + `compress_all_`, MPS/MPO `left_compress`, PEPS warmstarts, and boundary + fitting canonicalization/compression. + Keep the public API lazy and centralized in `src/pepsy/__init__.py` (`_SYMBOL_MODULES`, `_MODULE_EXPORTS`, `__all__`) per AGENTS.md. @@ -220,7 +272,109 @@ Goal: let PEPS/MPS workflows carry abelian symmetry and fermionic statistics. --- -## 5. Milestones (suggested order) +## 5. Workstream D — Differentiable truncated SVD / QR / eig + +Goal: make the decomposition layer usable for gradient optimization when +truncation is real, not just a post-hoc slice of a full decomposition. + +### D0. Audit and terminology +- Name the modes explicitly: + - `full_svd_slice`: compute full/compact SVD, then slice. Existing full SVD + AD is acceptable here and should remain the default baseline. + - `truncated_svd_ad`: compute only the kept subspace and use the paper's + truncated pullback. This is the efficient path Bram described in email. + - `stable_qr`: QR/LQ backward with phase/gauge stabilization, useful for + canonicalization but not a substitute for the missing SVD truncation term. +- Document which `quimb` methods map to each mode: + - `method="svd"` in current `quimb`: full SVD then trim. + - `method="svds"`, `"isvd"`, `"rsvd"`: low-rank candidates, but their AD + behavior is not automatically correct just because the forward is fast. + - `method="qr"`/`"lq"`: no singular-value truncation; useful for stable + canonicalization and gauge fixing. + +### D1. Centralize decomposition policy +- Add a small `pepsy.backends.decompositions` or `pepsy.backends.truncation` + module that owns: + - registration of full torch/JAX SVD/QR rules; + - optional registration of a `quimb` split driver such as + `method="svd:pepsy-trunc-ad"`; + - a context manager/config flag for choosing `full_svd_slice` vs + `truncated_svd_ad`. +- Prefer registering a new `quimb.tensor.decomp.register_split_driver(...)` + driver over changing every optimizer call site. Existing pepsy calls already + pass `max_bond`, `cutoff`, and `cutoff_mode` through Quimb. +- Keep the public optimizer APIs stable. Add keyword plumbing only where users + already pass decomposition/compression options; otherwise default to current + behavior. + +### D2. Implement the torch prototype first +- Start with torch because `PepsOptimizer` already has + `register_torch_svd=True` and `pepsy.backends.linalg_torch` already contains + custom SVD/QR autograd functions. +- Implement a custom `torch.autograd.Function` for true truncated SVD: + - forward accepts matrix `A`, target rank/cutoff policy, and method selector; + - returns `U_k`, `S_k`, `Vh_k` and enough metadata for backward; + - backward computes the paper's truncated-SVD VJP, including the discarded + spectrum contribution; + - for the first version, allow an internal full SVD in backward to verify the + formula, then replace the expensive part with Sylvester/linear-solve forms + when the tests pass. +- Handle adaptive cutoff carefully: + - the selected rank is piecewise constant and non-differentiable at threshold + crossings; + - stop gradients through the rank decision and test away from exact cutoff + crossings; + - record selected rank and discarded norm in `info` for diagnostics. + +### D3. JAX and backend follow-up +- After torch validation, implement the same primitive as a JAX `custom_vjp`. +- Do not promise NumPy gradients. NumPy remains the numerical/reference path. +- Symmray/block-sparse support should reuse the same array-level primitive per + dense block, following the TensorKit/MatrixAlgebraKit design rather than + writing tensor-object-specific AD rules. + +### D4. Integrate with Quimb without fighting Quimb +- Register a new split driver that mirrors `quimb.tensor.decomp.svd_truncated` + but calls the pepsy AD primitive instead of plain `xp.linalg.svd`. +- Respect Quimb's absorb modes (`both`, `left`, `right`, `None`) by applying + absorption after the custom primitive, just like Quimb's `_do_absorb`. +- Keep `info["error"]` semantics compatible with Quimb's current SVD driver. +- For call sites that cannot select `method=...`, add narrowly scoped adapter + kwargs in pepsy optimizers, e.g. `decomp_method="auto"` or + `compression_method="svd:pepsy-trunc-ad"`, defaulting to existing Quimb + behavior. + +### D5. Numerical tests +- Unit-level gradient checks: + - small complex matrices with fixed rank `k`; compare VJP to high-order + finite differences of gauge/phase-invariant scalar losses; + - repeated or near-repeated singular values; verify no blow-up beyond the + intended regularization; + - adaptive cutoff away from threshold crossings; verify selected-rank + diagnostics and gradients. +- Integration checks: + - MPS/MPO `mode="svd"` gate update with torch tensors and `max_bond < full`; + - PEPS warmstart split path; + - boundary-MPS fitting on a tiny PEPS; + - compare optimization steps against the current full-SVD baseline on small + examples where full SVD is feasible. +- Tests should avoid objectives depending on arbitrary singular-vector phases. + Use reconstructed factors, projectors/subspaces, norms, or physical + observables. + +### D6. Acceptance criteria +- Existing tests pass with default behavior unchanged. +- New focused tests pass under torch when optional torch is installed. +- Users can opt into the new path without touching Quimb internals directly. +- Docs explain when to use: + - current full-SVD AD: small matrices or when Quimb already computes all + singular vectors; + - true truncated AD: large matrices with small retained rank; + - finite differences: debugging/reference only. + +--- + +## 6. Milestones (suggested order) 1. **M0 — scaffolding (this PR):** `PLAN.md`, `history/`, `learning/`. 2. **M1 — quimb BP wrapper:** BP fixed point + norm via quimb, pepsy adapter + @@ -233,14 +387,20 @@ Goal: let PEPS/MPS workflows carry abelian symmetry and fermionic statistics. PEPS norm test. (C1–C3, C5) 6. **M5 — fermionic support:** fermionic constructors, phase handling, Hubbard energy test. (C4) -7. **M6 — docs:** `docs/howto/` + `docs/tutorials/` pages; update `docs/api/` +7. **M6 — AD decomposition prototype:** torch true-truncated-SVD primitive, + Quimb split driver registration, fixed-rank gradient tests, and one MPS/MPO + integration smoke. (D1–D5) +8. **M7 — AD decomposition integration:** adaptive-rank diagnostics, JAX + custom-VJP if needed, PEPS/boundary integration tests, docs explaining the + full-SVD vs true-truncated distinction. (D3–D6) +9. **M8 — docs:** `docs/howto/` + `docs/tutorials/` pages; update `docs/api/` and `tests/test_public_api.py` for any new public symbols. Each milestone is independently shippable and testable. --- -## 6. Cross-cutting requirements (from AGENTS.md) +## 7. Cross-cutting requirements (from AGENTS.md) - Use new package namespaces; never add old flat modules (`pepsy.core`, `pepsy.gates`, …). @@ -254,13 +414,14 @@ Each milestone is independently shippable and testable. - API/layout: `pytest -q tests/test_public_api.py tests/test_package_layout.py` - Boundary/contraction: `pytest -q tests/test_prepare_boundary_inputs.py` - Core/observables: `pytest -q tests/test_core_seed.py` + - Decomposition AD: `pytest -q tests/test_linalg_ad.py` - New BP: `pytest -q tests/test_bp.py` - New symmetry: `pytest -q tests/test_symmetry.py` - Docs: `sphinx-build -W -b html docs docs/_build/html` --- -## 7. Open questions / risks +## 8. Open questions / risks - BP fixed point may not exist / be unique for ground states of local Hamiltonians; loop expansion diverges for (near-)degenerate fixed points @@ -270,10 +431,13 @@ Each milestone is independently shippable and testable. - Exact division of labor between pepsy and quimb for BP (wrap vs. own) — decide in M1 after the quimb API audit; record the decision in `history/`. - symmray API is young (v0.2.x); pin a version and isolate behind the bridge. +- The efficient truncated-SVD VJP is subtle for complex tensors, adaptive + cutoffs, and degenerate singular values. Start with fixed-rank torch tests + and only wire into optimization once finite-difference checks are boring. --- -## 8. How agents should use this file +## 9. How agents should use this file 1. On session start: read this `PLAN.md` and the latest entry in `history/`. 2. Pick the current milestone; do the smallest shippable slice. diff --git a/README.md b/README.md index f80adc0..b2efbe5 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,8 @@ Current package version: `0.2.0` (from `pyproject.toml` / `pepsy.__version__`). - `optimizers/`: sweep, global, energy, MPS, MPO, and PEPS optimizers - `sampling/`: `MpsSampler` and related sampling utilities - `_internal/`: private formatting and utility helpers -- `examples/`: runnable examples (e.g. `MpsMagnetization/` 2D ITF Trotter MPS evolution) +- `examples/`: runnable examples, including `MpsMagnetization/` and direct + fermionic Symmray Fermi-Hubbard starters under `pepsy_examples/Fermi_Hubbard/` - `docs/`: Sphinx documentation source - `tests/`: package tests @@ -47,6 +48,37 @@ res = pepsy.contract_boundary(norm=norm, bdy=bdy, direction="y", n_iter=2) print(pepsy.__version__, res.cost) ``` +## Symmetric Fermionic States + +Pepsy includes optional Symmray-backed symmetric tensor-network wrappers. For +spinful Fermi-Hubbard work, `model="fermi_hubbard"` uses total particle-number +`U1`, while `model="fermi_hubbard_u1u1"` uses spin-resolved `U1U1` charges +`(N_up, N_down)`. + +For direct fermionic Fermi-Hubbard work, the main Pepsy/Symmray methods +reference is Gao et al., "Fermionic tensor network contraction for arbitrary +geometries", Phys. Rev. Research 7, 023193 (2025), +https://doi.org/10.1103/PhysRevResearch.7.023193. It motivates keeping +fermionic parity and leg-order metadata in Symmray arrays while letting quimb +choose graph-level contraction orders. + +```python +import pepsy as py + +psi = py.SymMPS.for_model( + "fermi_hubbard_u1u1", + 16, + bond_dim=4, + site_charge=py.site_charge_from_occupations([(1, 0), (0, 1)] * 8), +) + +assert psi.overall_charge() == (8, 8) + +ordering = psi.fermionic_ordering() +assert ordering["enabled"] +assert ordering["methods_reference"]["doi"] == "10.1103/PhysRevResearch.7.023193" +``` + ## Documentation Build docs locally: diff --git a/docs/api/tensors/core.md b/docs/api/tensors/core.md index 4447759..c59bc24 100644 --- a/docs/api/tensors/core.md +++ b/docs/api/tensors/core.md @@ -1,5 +1,10 @@ # `pepsy.tensors.core` +`reg_rel_svd_torch` is the preferred torch SVD registration for tensor-network +autodiff. It installs the relative-regularized SVD backward rule used by +`reg_complex_svd_torch`, and its CPU forward path falls back to SciPy `gesvd` +if `torch.linalg.svd` fails. + ```{eval-rst} .. automodule:: pepsy.tensors.core :members: diff --git a/docs/api/tensors/symmetric.md b/docs/api/tensors/symmetric.md index 7d1a811..e2c6f28 100644 --- a/docs/api/tensors/symmetric.md +++ b/docs/api/tensors/symmetric.md @@ -4,14 +4,16 @@ Physical sectors are charge maps: ``{charge: sector_size}``. For example, ``{0: 1, 1: 1}`` is a two-state U(1) or Z2 local space, while -``{0: 1, 1: 2, 2: 1}`` is the spinful Fermi-Hubbard local space -``|0>, |up>, |down>, |up down>``. +``{0: 1, 1: 2, 2: 1}`` is the spinful Fermi-Hubbard local space with total +particle-number U(1). For spin-resolved particle-number sectors use +``model="fermi_hubbard_u1u1"``, whose local charges are ``(n_up, n_down)``. ```python import pepsy as py spinless = py.default_physical_sectors("U1", 2) spinful = py.default_physical_sectors(model="fermi_hubbard") +spinful_spin_resolved = py.default_physical_sectors(model="fermi_hubbard_u1u1") psi = py.SymMPS.random( 6, @@ -23,8 +25,9 @@ psi = py.SymMPS.random( ``` The local ``site_charge`` pattern fixes the global sector represented by the -state. For U(1), ``psi.overall_charge()`` is the sum of local charges. For Z2, -``psi.overall_parity()`` is the same sum modulo two. +state. For U(1), ``psi.overall_charge()`` is the sum of local charges. For +``U1U1``, it is a pair such as ``(N_up, N_down)``. For Z2, +``psi.overall_parity()`` is the charge sum modulo two. ```python even = py.site_charge_uniform(0) @@ -43,6 +46,44 @@ peps.site_charges() peps.overall_parity() ``` +```python +half_filled_4x4 = py.site_charge_from_occupations( + [(1, 0), (0, 1)] * 8, +) + +fh = py.SymMPS.for_model( + "fermi_hubbard_u1u1", + 16, + bond_dim=4, + site_charge=half_filled_4x4, +) + +fh.overall_charge() # (8, 8) +``` + +For direct fermionic Fermi-Hubbard examples, use Gao et al., "Fermionic tensor +network contraction for arbitrary geometries", Phys. Rev. Research 7, 023193 +(2025), https://doi.org/10.1103/PhysRevResearch.7.023193 as the primary +methods reference. Pepsy should keep fermionic parity, additional Abelian +symmetries such as ``U1`` or ``U1U1``, and leg-order metadata in Symmray-backed +arrays, while relying on quimb/cotengra graph optimizers for contraction +ordering. + +Use ``state.fermionic_ordering()`` when a workflow needs the package-level +record of the graph and local order data carried by a symmetric state: + +```python +ordering = fh.fermionic_ordering() + +ordering["enabled"] # True for direct fermionic states +ordering["site_order"] # the site labels in tensor-network order +ordering["edge_order"] # the stored graph edge order +ordering["edges"][0]["index_directions"] +``` + +The same record is available as ``summary["fermionic_ordering"]`` from +``symmray_mps_summary`` and ``symmray_peps_summary``. + ## Time evolution Hamiltonians produce a canonical bundled gate stream, so the same stream can be @@ -56,6 +97,15 @@ gates = ham.gate_stream(0.01) psi.time_evolve_mps_optimizer(0.01, hamiltonian=ham, chi=16, mode="mpo") ``` +For Symmray-backed MPS gate streams, ``mode="swap"`` and ``mode="svd"`` use +quimb's block-aware auto-swap split path for nonlocal 1D gate streams such as a +row-major square lattice. ``mode="mpo"`` uses its usual sub-MPO compression for +nearest-neighbor gates and falls back to the same Symmray auto-swap path for +nonlocal gates, because the current quimb/Symmray sub-MPO path mixes in dense +helper tensors. ``mode="exact"`` is useful as a small-system reference, and +``mode="dmrg"`` works when the initial symmetric MPS was built with enough +block-sparse bond capacity, for example ``bond_dim >= chi``. + ```python peps = py.SymPEPS.for_model("itf", 4, 4, bond_dim=2) ham = peps.build_hamiltonian(jx=-1.0, hz=-0.5) @@ -103,6 +153,93 @@ peps.measure( ) ``` +## Inspecting Symmray blocks + +Use ``symmray_block_summary`` to inspect the charge sectors and stored block +dimensions of any Symmray-backed local operator or tensor. Use +``draw_symmray_blocks`` for a lightweight ``quimb.schematic`` drawing of the +same information, including the stored-vs-dense entry count. + +For whole MPS states, use ``symmray_mps_summary`` and ``draw_symmray_mps``. For +PEPS, use ``symmray_peps_summary`` and ``draw_symmray_peps``. These expose the +scientific structure that is usually hidden in a dense drawing: each site +tensor's block count, physical charge sectors, virtual-bond sector maps, and +aggregate block-sparse storage density. The default schematics follow the +compact quimb style with tensor nodes, physical legs, and bond arrows; labels +and diagnostics are opt-in. + +In the detailed drawing mode, ``T_i`` is the site tensor, ``q`` is that tensor's +Symmray charge, ``B`` is the number of stored block sectors, ``e_i`` is a +virtual bond, and ``q_e``/``q_p`` are the virtual/physical charge-sector maps. +Bond labels include the two local index orientations, for example +``out->in``, so the charge-flow convention is visible on the same line as the +bond dimension. Diagnostics include both ``charge_total`` and ``Q_total``; for +``Z2`` states ``Q_total`` is reduced modulo two. Colored block tiles are +available with ``show_blocks=True`` for a focused block-sector view, but the +overview diagrams leave them off by default and show ``B`` instead. +The MPS/PEPS summary dictionaries also expose ``charge_total`` and ``Q_total`` +alongside the legacy ``total_charge`` key. + +```python +psi = py.SymMPS.for_model( + "itf", + 4, + bond_dim=2, + site_charge=py.site_charge_from_occupations([0] * 4), +) +rz_gate = psi.operator_from_dense(py.rz(0.1), charge=0, sites=1) + +summary = py.symmray_block_summary(rz_gate) +summary["blocks"] + +gate_drawing = py.draw_symmray_blocks(rz_gate, title="Z2 RZ gate") +display(gate_drawing.fig) + +mps_summary = py.symmray_mps_summary(psi.tn) +mps_summary["bonds"] + +mps_drawing = py.draw_symmray_mps(psi.tn, title="Symmray ITF MPS") +display(mps_drawing.fig) + +detailed_mps_drawing = py.draw_symmray_mps( + psi.tn, + title="Symmray ITF MPS with dimensions", + show_bond_labels=True, + show_phys_labels=True, + show_diagnostics=True, +) +display(detailed_mps_drawing.fig) + +peps = py.SymPEPS.for_model("itf", 3, 3, bond_dim=2) +peps_summary = py.symmray_peps_summary(peps) +peps_summary["bonds"] + +peps_drawing = py.draw_symmray_peps( + peps, + title="Symmray ITF PEPS", + show_bond_labels=True, + show_phys_labels=True, + show_diagnostics=True, +) +display(peps_drawing.fig) +``` + +At the end of a Symmray ``MpsOptimizer`` notebook, pass the optimized chain +directly: + +```python +opt = py.MpsOptimizer(psi.tn.copy(), gates, chi=8, mode="mpo") +opt.run(progbar=False) + +py.draw_symmray_mps( + opt.p, + title="Final Symmray ITF MPS", + center="middle", + show_bond_labels=True, + show_diagnostics=True, +) +``` + ```{eval-rst} .. automodule:: pepsy.tensors.symmetric :members: diff --git a/examples/MpsMagnetization/helper.py b/examples/MpsMagnetization/helper.py index 2e41e6a..7f45dff 100644 --- a/examples/MpsMagnetization/helper.py +++ b/examples/MpsMagnetization/helper.py @@ -13,22 +13,29 @@ def build_lattice(Lx, Ly, coupling_j, field_h, cyclic=True, lattice="square", mo Returns ------- - dict with keys: mpo_H, edges_1d, sites, mapper, res (full build_itf_lattice output). + dict with keys: mpo_H, edges_2d, edges_1d, sites, mapper, res + (full build_itf_lattice output), and the 1D <-> 2D maps. """ - mapper = OneDMap(Lx, Ly, mode="hilbert") + mapper = OneDMap(Lx, Ly, mode=mode) res = py.ham_tn.build_itf_lattice( L_x=Lx, L_y=Ly, lattice=lattice, cyclic=cyclic, J=coupling_j, field=field_h, mapper=mapper, ) mpo_H = res["mpo"] + edges_2d = res["edges"] edges_1d = res["edges_1d"] - sites = sorted({(site,) for edge in edges_1d for site in edge}) + one_d_to_two_d = res["one_d_to_lattice"] + two_d_to_one_d = res["lattice_to_one_d"] + sites = tuple((site,) for site in sorted(one_d_to_two_d)) return { "mpo_H": mpo_H, + "edges_2d": edges_2d, "edges_1d": edges_1d, "sites": sites, "mapper": mapper, + "one_d_to_two_d": one_d_to_two_d, + "two_d_to_one_d": two_d_to_one_d, "res": res, } diff --git a/examples/MpsMagnetization/mps_simulator.ipynb b/examples/MpsMagnetization/mps_simulator.ipynb index fcbd7a0..799d7f9 100644 --- a/examples/MpsMagnetization/mps_simulator.ipynb +++ b/examples/MpsMagnetization/mps_simulator.ipynb @@ -13,14 +13,20 @@ "source": [ "## MPS Magnetization — 2D ITF Trotter Evolution\n", "\n", - "Time-evolve a product state under the 2D transverse-field Ising model\n", - "using pepsy's `MpsOptimizer` with different modes and backends.\n", + "Time-evolve a product state under the 2D transverse-field Ising model using one clear MPS optimizer path.\n", "\n", - "**Parameters** (matching `compare_exact_vs_color.py`):\n", - "- 4×4 square lattice, PBC\n", - "- H = J ΣZZ + h ΣX, J = −1, h = 2\n", + "This notebook keeps the geometry and mapping visible:\n", + "\n", + "1. build the 2D square-lattice edges,\n", + "2. map those 2D edges to 1D MPS sites with one `OneDMap`,\n", + "3. build an explicit 2nd-order Trotter gate stream,\n", + "4. evolve with `MpsOptimizer(mode=\"dmrg\")`.\n", + "\n", + "**Parameters**:\n", + "- 4x4 square lattice, PBC\n", + "- H = J ΣZZ + h ΣX, J = -1, h = 2\n", "- 2nd-order Trotter with dt = 0.25\n", - "- Initial state: intermediate-temperature quench (arXiv:2503.20870)" + "- Initial state: intermediate-temperature quench (arXiv:2503.20870)\n" ] }, { @@ -36,14 +42,14 @@ "import pepsy as py\n", "from pepsy import tensors\n", "from helper import (\n", - " build_lattice, build_initial_state, build_trotter_gates,\n", + " build_lattice, build_initial_state,\n", " build_mpo_z_sq, measure_energy, measure_z, measure_z_sq,\n", ")\n", "\n", "try:\n", " import torch\n", "except Exception:\n", - " torch = None" + " torch = None\n" ] }, { @@ -83,7 +89,9 @@ "id": "d4df0e16", "metadata": {}, "source": [ - "### Lattice & Hamiltonian" + "### Lattice, Mapping, And Hamiltonian\n", + "\n", + "`build_lattice(...)` uses Pepsy's ITF lattice builder. For `lattice=\"square\"`, Pepsy delegates the 2D edge generation to `qtn.edges_2d_square(...)`, then maps every 2D coordinate to the 1D MPS chain with the single `mapper` below.\n" ] }, { @@ -93,7 +101,7 @@ "metadata": {}, "outputs": [], "source": [ - "# ── Physical parameters ──────────────────────────────────────────────────────\n", + "# -- Physical parameters -----------------------------------------------------\n", "Lx, Ly = 4, 4\n", "L = Lx * Ly\n", "coupling_j = -1.0 # J\n", @@ -101,20 +109,51 @@ "cyclic = True\n", "lattice = \"square\"\n", "dt = 0.25 # Trotter step size\n", - "\n", - "# ── Build lattice ────────────────────────────────────────────────────────────\n", - "# mode: \"snake\" | \"snake-row-major\" | \"row-major\" | \"col-major\" |\n", - "# \"hilbert\" | \"hilbert-row-major\" | \"diag\"\n", - "lat = build_lattice(Lx, Ly, coupling_j, field_h, cyclic=cyclic, lattice=lattice, mode=\"hilbert\")\n", + "mapping_mode = \"hilbert\"\n", + "\n", + "# -- Build lattice, MPO Hamiltonian, and one shared 2D <-> 1D mapper ----------\n", + "lat = build_lattice(\n", + " Lx,\n", + " Ly,\n", + " coupling_j,\n", + " field_h,\n", + " cyclic=cyclic,\n", + " lattice=lattice,\n", + " mode=mapping_mode,\n", + ")\n", "mpo_H = lat[\"mpo_H\"]\n", + "edges_2d = lat[\"edges_2d\"]\n", "edges_1d = lat[\"edges_1d\"]\n", "sites = lat[\"sites\"]\n", "mapper = lat[\"mapper\"]\n", + "one_d_to_two_d = lat[\"one_d_to_two_d\"]\n", + "two_d_to_one_d = lat[\"two_d_to_one_d\"]\n", + "\n", + "# Same square-lattice edges, checked directly against quimb's generator.\n", + "if lattice == \"square\":\n", + " edges_2d_direct = tuple(\n", + " tuple(tuple(site) for site in edge)\n", + " for edge in qtn.edges_2d_square(Lx, Ly, cyclic=cyclic)\n", + " )\n", + " assert {frozenset(edge) for edge in edges_2d} == {\n", + " frozenset(edge) for edge in edges_2d_direct\n", + " }\n", "\n", - "mapper.show(title=\"Hilbert-curve 2D → 1D mapping\")\n", - "print(f\"Lattice: {Lx}×{Ly}, L={L}, {'PBC' if cyclic else 'OBC'}\")\n", - "print(f\"H = J·ΣZZ + h·ΣX, J={coupling_j}, h={field_h}\")\n", - "print(f\"Trotter dt = {dt}\")" + "# Same mapping used for Hamiltonian MPO and for the Trotter gate locations.\n", + "edges_1d_from_map = tuple(\n", + " (two_d_to_one_d[tuple(site_a)], two_d_to_one_d[tuple(site_b)])\n", + " for site_a, site_b in edges_2d\n", + ")\n", + "assert edges_1d == edges_1d_from_map\n", + "assert sites == tuple((site,) for site in range(L))\n", + "\n", + "mapper.show(title=f\"{mapping_mode} 2D -> 1D mapping\")\n", + "print(f\"Lattice: {Lx}x{Ly}, L={L}, {'PBC' if cyclic else 'OBC'}\")\n", + "print(f\"H = J*ΣZZ + h*ΣX, J={coupling_j}, h={field_h}\")\n", + "print(f\"Mapping mode: {mapping_mode}\")\n", + "print(f\"2D edges: {len(edges_2d)} | mapped 1D edges: {len(edges_1d)}\")\n", + "print(f\"First 2D edge {edges_2d[0]} -> 1D edge {edges_1d[0]}\")\n", + "print(f\"Trotter dt = {dt}\")\n" ] }, { @@ -151,9 +190,15 @@ "id": "8bc80571", "metadata": {}, "source": [ - "### 2nd-order Trotter gates\n", + "### Explicit 2nd-Order Trotter Gates\n", + "\n", + "For one Trotter step we use Strang splitting:\n", "\n", - "Strang splitting: RX(h·dt/2) → RZZ(2J·dt) → RX(h·dt/2)" + "1. apply the X half-step on every site,\n", + "2. apply the ZZ full-step on every mapped lattice edge,\n", + "3. apply the X half-step again.\n", + "\n", + "Pepsy's `rx(theta)` convention is `exp(-i theta X / 2)`, so the half-step `exp(-i h dt X / 2)` is `py.rx(h * dt)`. Likewise, `py.rzz(theta)` is `exp(-i theta ZZ / 2)`, so the full bond step `exp(-i J dt ZZ)` is `py.rzz(2 * J * dt)`.\n" ] }, { @@ -163,9 +208,22 @@ "metadata": {}, "outputs": [], "source": [ - "# ── Build 2nd-order Trotter gate stream ──────────────────────────────────────\n", - "gates_trotter = build_trotter_gates(sites, edges_1d, field_h, coupling_j, dt, to_backend)\n", - "print(f\"Trotter gates: {len(gates_trotter)} total\")" + "# -- Build one explicit 2nd-order Trotter gate stream ------------------------\n", + "rx_half = to_backend(py.rx(field_h * dt))\n", + "rzz_full = to_backend(py.rzz(2.0 * coupling_j * dt))\n", + "\n", + "gates_x_left = [(rx_half, site) for site in sites]\n", + "gates_zz = [(rzz_full, edge) for edge in edges_1d]\n", + "gates_x_right = [(rx_half, site) for site in sites]\n", + "\n", + "gates_trotter = gates_x_left + gates_zz + gates_x_right\n", + "\n", + "print(f\"X half-step gates: {len(gates_x_left)}\")\n", + "print(f\"ZZ edge gates: {len(gates_zz)}\")\n", + "print(f\"X half-step gates: {len(gates_x_right)}\")\n", + "print(f\"Total gates: {len(gates_trotter)}\")\n", + "print(f\"Example one-site gate location: {gates_x_left[0][1]}\")\n", + "print(f\"Example two-site gate location: {gates_zz[0][1]}\")\n" ] }, { @@ -173,7 +231,9 @@ "id": "dd3f1a09", "metadata": {}, "source": [ - "### MPS evolution parameters" + "### MPS Evolution Parameters\n", + "\n", + "Keep one optimizer mode in focus here: `dmrg`. The other MPS modes are useful diagnostics, but comparing modes is a separate experiment from checking the lattice mapping and gate stream.\n" ] }, { @@ -183,105 +243,17 @@ "metadata": {}, "outputs": [], "source": [ - "# ── Evolution parameters ─────────────────────────────────────────────────────\n", + "# -- Evolution parameters ----------------------------------------------------\n", "n_steps = 20 # number of Trotter steps\n", "chi = 120 # MPS bond dimension\n", + "mode = \"dmrg\"\n", "\n", - "mode = \"dmrg\" # MPS optimization mode: \"mpo\", \"dmrg\", \"svd\", \"swap\", \"exact\"\n", - "\n", - "# DMRG-specific (only used when mode=\"dmrg\"):\n", + "# DMRG-specific controls:\n", "n_iter = 5 # inner FIT iterations per 2q gate\n", "k_2q_batch = 1 # batch this many 2q gates into one FIT window\n", "\n", "print(f\"n_steps = {n_steps}, T = {n_steps * dt}, chi = {chi}\")\n", - "print(f\"Mode: {mode}\")" - ] - }, - { - "cell_type": "markdown", - "id": "nonunitary-check-heading", - "metadata": {}, - "source": [ - "### Non-unitary MPS gate-stream check\n", - "\n", - "A compact CPU check for the `non_unitary=True` path before the heavier time evolution below. It compares `dmrg` and `mpo`, separates represented norm from raw tensor-data norm, and enables the true normalized-overlap diagnostic explicitly.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "nonunitary-check-code", - "metadata": {}, - "outputs": [], - "source": [ - "# -- Non-unitary stream smoke check -----------------------------------------\n", - "H = np.array([[1.0, 1.0], [1.0, -1.0]], dtype=complex) / np.sqrt(2.0)\n", - "non_unitary_filter = np.diag([1.0, 0.5, 0.5, 2.0]).astype(complex)\n", - "gates_non_unitary = [\n", - " (H, (0,)),\n", - " (H, (1,)),\n", - " (non_unitary_filter, (0, 1)),\n", - "]\n", - "\n", - "def run_non_unitary_check(mode_check):\n", - " psi_check = qtn.MPS_computational_state(\"0000\", dtype=\"complex128\")\n", - " engine_check = py.MpsOptimizer(\n", - " psi_check, gates_non_unitary, chi=1, mode=mode_check,\n", - " contraction_opt=\"auto-hq\",\n", - " )\n", - " engine_check.run(\n", - " progbar=False, cutoff=1e-12, n_iter=2,\n", - " non_unitary=True, normalize_every=1,\n", - " track_infidelity=True,\n", - " )\n", - " sample = engine_check.get_infidelity_samples()[-1]\n", - " data_state = engine_check.p.copy()\n", - " data_state.exponent = 0.0\n", - " return {\n", - " \"mode\": mode_check,\n", - " \"represented_norm\": float(np.real(engine_check.p.norm())),\n", - " \"data_norm\": float(np.real(data_state.norm())),\n", - " \"exponent\": float(np.real(getattr(engine_check.p, \"exponent\", 0.0))),\n", - " \"geometric_fidelity\": sample[\"geometric_fidelity\"],\n", - " \"loss_fidelity\": engine_check.losses[-1],\n", - " \"local_infidelity\": sample[\"local_infidelity\"],\n", - " \"infidelity\": engine_check.get_infidelities()[-1],\n", - " \"normalizations\": len(engine_check.get_normalizations()),\n", - " \"state\": engine_check.p.copy(),\n", - " }\n", - "\n", - "non_unitary_checks = [run_non_unitary_check(m) for m in (\"dmrg\", \"mpo\")]\n", - "\n", - "# With one compressed update, the tracked true metric equals the final target infidelity.\n", - "target = qtn.MPS_computational_state(\"0000\", dtype=\"complex128\")\n", - "for gate_i, where_i in gates_non_unitary[:2]:\n", - " py.gate(target, gate_i, where_i, contract=True, cutoff=1e-12, inplace=True)\n", - "target = py.gate(\n", - " target, non_unitary_filter, (0, 1),\n", - " contract=\"reduce-split\", cutoff=1e-12, cutoff_mode=\"rel\", inplace=False,\n", - ")\n", - "for row in non_unitary_checks:\n", - " target_fidelity = py.tn_fidelity(row[\"state\"], target, contraction_opt=\"auto-hq\")\n", - " row[\"target_infidelity\"] = 1.0 - float(np.real(target_fidelity))\n", - " assert np.isclose(row[\"infidelity\"], row[\"target_infidelity\"])\n", - " assert np.isclose(row[\"geometric_fidelity\"], row[\"loss_fidelity\"])\n", - "\n", - "dmrg_check = next(row for row in non_unitary_checks if row[\"mode\"] == \"dmrg\")\n", - "mpo_check = next(row for row in non_unitary_checks if row[\"mode\"] == \"mpo\")\n", - "assert np.isclose(dmrg_check[\"infidelity\"], mpo_check[\"infidelity\"])\n", - "assert np.isclose(dmrg_check[\"represented_norm\"], mpo_check[\"represented_norm\"])\n", - "\n", - "for row in non_unitary_checks:\n", - " assert np.isclose(row[\"data_norm\"], 1.0)\n", - " assert 0.0 <= row[\"infidelity\"] <= 1.0\n", - " print(\n", - " f\"{row['mode']}: represented norm={row['represented_norm']:.12f}, \"\n", - " f\"data norm={row['data_norm']:.12f}, exponent={row['exponent']:.6f}, \"\n", - " f\"geometric fidelity={row['geometric_fidelity']:.6e}, \"\n", - " f\"local infidelity={row['local_infidelity']:.6e}, \"\n", - " f\"cumulative infidelity={row['infidelity']:.6e}, \"\n", - " f\"normalizations={row['normalizations']}\"\n", - " )\n" + "print(f\"Mode: {mode}\")\n" ] }, { @@ -336,12 +308,16 @@ "energy = [measure_energy(engine.p.copy(), mpo_H, L, optimizer)]\n", "z_sq = [measure_z_sq(engine.p.copy(), mpo_z_sq_offdiag, diagonal_shift, optimizer)]\n", "z_mag = [measure_z(engine.p.copy(), mpo_z, optimizer)]\n", - "norm_proxy_per_step = []\n", + "norm_proxy_per_step = []\n", "\n", "pbar = tqdm(range(n_steps), desc=mode)\n", "for step in pbar:\n", - " engine.run(progbar=False, mode=mode, fidelity_samples=10,\n", - " n_iter=n_iter, k_2q_batch=k_2q_batch)\n", + " engine.run(\n", + " progbar=False,\n", + " fidelity_samples=10,\n", + " n_iter=n_iter,\n", + " k_2q_batch=k_2q_batch,\n", + " )\n", "\n", " t_now = (step + 1) * dt\n", " psi_t = engine.p.copy()\n", @@ -351,16 +327,19 @@ " z_mag.append(measure_z(psi_t, mpo_z, optimizer))\n", "\n", " loss = engine.losses[-1]\n", - " norm_proxy_per_step.append(loss)\n", - " pbar.set_postfix({\"norm\": f\"{loss:.6f}\", \"E/L\": f\"{energy[-1]:.4f}\",\n", - " \"Z²\": f\"{z_sq[-1]:.4f}\", \"Z\": f\"{z_mag[-1]:.4f}\"})\n", + " norm_proxy_per_step.append(loss)\n", + " pbar.set_postfix({\n", + " \"norm\": f\"{loss:.6f}\",\n", + " \"E/L\": f\"{energy[-1]:.4f}\",\n", + " \"Z^2\": f\"{z_sq[-1]:.4f}\",\n", + " \"Z\": f\"{z_mag[-1]:.4f}\",\n", + " })\n", "\n", "times = np.array(times)\n", "energy = np.array(energy)\n", "z_sq = np.array(z_sq)\n", "z_mag = np.array(z_mag)\n", - "\n", - "norm_proxy_per_step = np.array(norm_proxy_per_step)\n", + "norm_proxy_per_step = np.array(norm_proxy_per_step)\n", "print(f\"\\nDone. {len(times)} measurement points.\")\n" ] }, @@ -431,24 +410,24 @@ "metadata": {}, "outputs": [], "source": [ - "# ── Norm proxy per Trotter step ──────────────────────────────────────────────\n", - "step_indices = np.arange(1, n_steps + 1)\n", - "\n", - "fig, ax = plt.subplots(figsize=(9, 3.5))\n", - "ax.plot(step_indices, norm_proxy_per_step, \">-\", lw=1.2, markersize=4,\n", - " color=\"#7b2cbf\", alpha=0.85)\n", - "ax.axhline(norm_proxy_per_step.mean(), ls=\"--\", lw=0.8, color=\"gray\", alpha=0.6,\n", - " label=f\"mean = {norm_proxy_per_step.mean():.6f}\")\n", - "\n", - "# Auto y-limits to show variations clearly\n", - "ymin = norm_proxy_per_step.min() - 0.005\n", - "ymax = min(norm_proxy_per_step.max() + 0.005, 1.0)\n", - "ax.set_ylim(ymin, ymax)\n", - "\n", - "ax.set_xlabel(\"Trotter step\")\n", - "ax.set_ylabel(\"Norm proxy\")\n", - "ax.set_title(f\"Norm proxy per step — {mode}, $\\\\chi$={chi}, $dt$={dt}\",\n", - " fontweight=\"bold\")\n", + "# ── Norm proxy per Trotter step ──────────────────────────────────────────────\n", + "step_indices = np.arange(1, n_steps + 1)\n", + "\n", + "fig, ax = plt.subplots(figsize=(9, 3.5))\n", + "ax.plot(step_indices, norm_proxy_per_step, \">-\", lw=1.2, markersize=4,\n", + " color=\"#7b2cbf\", alpha=0.85)\n", + "ax.axhline(norm_proxy_per_step.mean(), ls=\"--\", lw=0.8, color=\"gray\", alpha=0.6,\n", + " label=f\"mean = {norm_proxy_per_step.mean():.6f}\")\n", + "\n", + "# Auto y-limits to show variations clearly\n", + "ymin = norm_proxy_per_step.min() - 0.005\n", + "ymax = min(norm_proxy_per_step.max() + 0.005, 1.0)\n", + "ax.set_ylim(ymin, ymax)\n", + "\n", + "ax.set_xlabel(\"Trotter step\")\n", + "ax.set_ylabel(\"Norm proxy\")\n", + "ax.set_title(f\"Norm proxy per step — {mode}, $\\\\chi$={chi}, $dt$={dt}\",\n", + " fontweight=\"bold\")\n", "ax.grid(True, alpha=0.25, linestyle=\"--\")\n", "ax.set_xlim(0.5, n_steps + 0.5)\n", "ax.spines[\"top\"].set_visible(False)\n", @@ -457,8 +436,8 @@ "plt.tight_layout()\n", "plt.show()\n", "\n", - "print(f\"Norm proxy per step: min={norm_proxy_per_step.min():.8f}, \"\n", - " f\"max={norm_proxy_per_step.max():.8f}, mean={norm_proxy_per_step.mean():.8f}\")" + "print(f\"Norm proxy per step: min={norm_proxy_per_step.min():.8f}, \"\n", + " f\"max={norm_proxy_per_step.max():.8f}, mean={norm_proxy_per_step.mean():.8f}\")" ] }, { diff --git a/examples/pepsy_examples/Fermi_Hubbard/README.md b/examples/pepsy_examples/Fermi_Hubbard/README.md new file mode 100644 index 0000000..0b2ac3d --- /dev/null +++ b/examples/pepsy_examples/Fermi_Hubbard/README.md @@ -0,0 +1,54 @@ +# Direct Fermionic Fermi-Hubbard Examples + +This folder contains Pepsy/Symmray examples for Fermi-Hubbard simulations that +work directly in fermionic tensor-network space. These examples do not use the +fermion-to-qubit encoding from arXiv:2511.02125. + +## Main methods reference + +Use Gao et al., "Fermionic tensor network contraction for arbitrary +geometries", Phys. Rev. Research 7, 023193 (2025), +https://doi.org/10.1103/PhysRevResearch.7.023193 as the main Pepsy/Symmray +reference for direct fermionic Fermi-Hubbard tensor networks. + +Implementation cues from that paper: + +- Store parity, Abelian charge labels, and leg order in the fermionic array + object, not in a separate fermion-to-qubit mapping layer. +- Prefer local edge/order metadata for arbitrary graphs; square-lattice PEPS, + snake MPS embeddings, checkerboards, bilayers, and random graphs should all + flow through the same fermionic tensor path when possible. +- Let quimb/cotengra optimize the tensor-network graph contraction order; the + fermionic signs are handled by Symmray's swap/parity algebra. +- Keep dense fallbacks as small-system references only, and verify that gates, + observables, and evolved states remain Symmray/FermionicArray-backed. + +## Current notebook + +- `half_filled_4x4_direct_fermions.ipynb` + +The notebook builds both MPS and PEPS versions of a 4 x 4 spinful +Fermi-Hubbard model with: + +- symmetry: `U1U1` +- local sectors: `(0, 0)`, `(0, 1)`, `(1, 0)`, `(1, 1)` +- half-filled charge: `(N_up, N_down) = (8, 8)` +- parameters: `t = 1`, `U/t = 8` +- no Jordan-Wigner, compact, or Octagon fermion-to-qubit mapping + +The energy and imaginary-time cells are off by default so the notebook opens +quickly. Turn them on after the construction and block summaries look right. + +## Paper targets to test next + +For the doped checkerboard case in arXiv:2511.02125: + +- lattice: 6 x 6 split into 2 x 2 plaquettes +- direct-fermion charge: `(15, 15)` +- doping: 30 particles on 36 sites, i.e. six holes from half filling +- weak-coupling point: `t_prime = 0`, `U/t = 2` +- exact shifted energy target: `/N - U/4 = -1.27367` +- d-wave average target: `0.108` theory, `0.079 +/- 0.005` experiment + +The next implementation step is adding weighted checkerboard Hubbard edges and +fermionic singlet-pair observables so those numbers can be tested directly. diff --git a/examples/pepsy_examples/Fermi_Hubbard/half_filled_4x4_direct_fermions.ipynb b/examples/pepsy_examples/Fermi_Hubbard/half_filled_4x4_direct_fermions.ipynb new file mode 100644 index 0000000..793be57 --- /dev/null +++ b/examples/pepsy_examples/Fermi_Hubbard/half_filled_4x4_direct_fermions.ipynb @@ -0,0 +1,386 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 4x4 Fermi-Hubbard Starter: Direct Fermionic Pepsy/Symmray\n", + "\n", + "This notebook is a small starting point for direct fermionic Fermi-Hubbard simulations in Pepsy/Symmray. The main methods reference is Gao et al., *Fermionic tensor network contraction for arbitrary geometries*, Phys. Rev. Research 7, 023193 (2025), https://doi.org/10.1103/PhysRevResearch.7.023193. The square-lattice physics settings are taken from arXiv:2511.02125, whose hardware workflow uses a fermion-to-qubit encoding; this notebook does **not** use that mapping. It works directly with fermionic Symmray tensor-network objects through Pepsy.\n", + "\n", + "The first target is intentionally modest: build a 4x4 spinful fermionic PEPS at half filling, construct the Hubbard Hamiltonian with `t = 1` and `U/t = 8`, inspect the block-sparse structure, and optionally run one tiny imaginary-time step." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "import os\n", + "\n", + "os.environ.setdefault(\"NUMBA_CACHE_DIR\", \"/tmp/numba_cache_pepsy\")\n", + "os.environ.setdefault(\"MPLCONFIGDIR\", \"/tmp/mplconfig_pepsy\")\n", + "os.makedirs(os.environ[\"NUMBA_CACHE_DIR\"], exist_ok=True)\n", + "os.makedirs(os.environ[\"MPLCONFIGDIR\"], exist_ok=True)\n", + "\n", + "import numpy as np\n", + "\n", + "import pepsy as py\n", + "\n", + "try:\n", + " import symmray as sr\n", + "except ImportError as exc:\n", + " raise ImportError(\n", + " \"This example needs the optional dependency `symmray`. \"\n", + " \"Install it in the Pepsy environment before running the notebook.\"\n", + " ) from exc\n", + "\n", + "print(\"pepsy\", py.__version__)\n", + "print(\"symmray\", getattr(sr, \"__version__\", \"installed\"))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Model Settings\n", + "\n", + "The paper's half-filled Hubbard section uses a periodic 6x6 square lattice and mixed periodic boundary-condition averaging. For a first direct-fermion tensor-network smoke test, this notebook starts with a 4x4 lattice. `CYCLIC = False` keeps the first contraction cheaper; set it to `True` once the small open-boundary workflow runs cleanly." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Lx = 4\n", + "Ly = 4\n", + "t = 1.0\n", + "U = 8.0\n", + "D = 2\n", + "CYCLIC = False\n", + "SEED = 251102125\n", + "\n", + "num_sites = Lx * Ly\n", + "num_spin_orbitals = 2 * num_sites\n", + "target_up = num_sites // 2\n", + "target_down = num_sites // 2\n", + "target_particles = target_up + target_down\n", + "target_charge = (target_up, target_down)\n", + "\n", + "phys_sectors = py.default_physical_sectors(model=\"fermi_hubbard_u1u1\")\n", + "half_filled_charges = {\n", + " (x, y): (1, 0) if (x + y) % 2 == 0 else (0, 1)\n", + " for x in range(Lx)\n", + " for y in range(Ly)\n", + "}\n", + "site_charge = py.site_charge_from_occupations(half_filled_charges)\n", + "\n", + "print(\"local sectors:\", phys_sectors)\n", + "print(\"sites:\", num_sites)\n", + "print(\"spin orbitals:\", num_spin_orbitals)\n", + "print(\"target total particle number:\", target_particles)\n", + "print(\"target (N_up, N_down):\", target_charge)\n", + "print(\"filling N / (2 sites):\", target_particles / num_spin_orbitals)\n", + "print(\"boundary:\", \"periodic\" if CYCLIC else \"open smoke test\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Half Filling And Symmetry Sector\n", + "\n", + "For the spinful Hubbard model, each site has two spin orbitals. Half filling means `N_particles = N_sites`, not `2 * N_sites`. Here that is 16 particles on 16 lattice sites, or 16 particles in 32 spin orbitals.\n", + "\n", + "The symmetry used below is spin-resolved `U1U1`, with local charges `(n_up, n_down)`. The notebook fixes `N_up = 8` and `N_down = 8`, so the total particle number is `16` on `16` lattice sites. The simple alternating charge pattern is only a sector choice for initialization; it is not a fermion-to-qubit mapping." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Build A Direct Fermionic PEPS\n", + "\n", + "Each lattice site has the spinful Hubbard local basis with charge sectors `(0, 0)`, `(0, 1)`, `(1, 0)`, and `(1, 1)`. The `site_charge` map fixes the global `U1U1` sector to `(8, 8)`, i.e. 16 fermions on this 4x4 lattice." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "psi = py.SymPEPS.random(\n", + " Lx,\n", + " Ly,\n", + " symmetry=\"U1U1\",\n", + " fermionic=True,\n", + " phys_dim=phys_sectors,\n", + " site_charge=site_charge,\n", + " bond_dim=D,\n", + " cyclic=CYCLIC,\n", + " seed=SEED,\n", + " dtype=\"complex128\",\n", + " contraction_opt=\"auto-hq\",\n", + ")\n", + "\n", + "ham = py.SymHamiltonian.from_edges(\n", + " \"fermi_hubbard_u1u1\",\n", + " \"U1U1\",\n", + " psi.edges,\n", + " t=t,\n", + " U=U,\n", + " mu=0.0,\n", + ")\n", + "\n", + "summary = py.symmray_peps_summary(psi)\n", + "\n", + "print(\"sites:\", psi.num_sites)\n", + "print(\"edges / local Hamiltonian terms:\", len(ham.terms))\n", + "print(\"fermionic tensors:\", psi.fermionic)\n", + "print(\"total U1U1 charge:\", psi.overall_charge())\n", + "print(\"max PEPS bond dim:\", summary[\"max_bond_dim\"])\n", + "print(\"stored / dense entries:\", summary[\"total_stored_size\"], \"/\", summary[\"total_dense_size\"])\n", + "print(\"block storage density:\", f\"{summary['density']:.3f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Same Half-Filled Sector As An MPS\n", + "\n", + "This builds the same direct-fermion half-filled 4x4 problem as an MPS. The MPS uses a snake ordering for the sites; the Hamiltonian still contains the 2D square-lattice edges, so some gates are nonlocal in the 1D order." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "snake_sites = []\n", + "for x in range(Lx):\n", + " ys = range(Ly) if x % 2 == 0 else range(Ly - 1, -1, -1)\n", + " for y in ys:\n", + " snake_sites.append((x, y))\n", + "\n", + "coord_to_mps = {site: i for i, site in enumerate(snake_sites)}\n", + "mps_edges = tuple(\n", + " (coord_to_mps[a], coord_to_mps[b])\n", + " for a, b in psi.edges\n", + ")\n", + "\n", + "psi_mps = py.SymMPS.random(\n", + " num_sites,\n", + " symmetry=\"U1U1\",\n", + " fermionic=True,\n", + " phys_dim=phys_sectors,\n", + " site_charge=py.site_charge_from_occupations([\n", + " half_filled_charges[site]\n", + " for site in snake_sites\n", + " ]),\n", + " bond_dim=D,\n", + " seed=SEED,\n", + " dtype=\"complex128\",\n", + ")\n", + "ham_mps = py.SymHamiltonian.from_edges(\n", + " \"fermi_hubbard_u1u1\",\n", + " \"U1U1\",\n", + " mps_edges,\n", + " t=t,\n", + " U=U,\n", + " mu=0.0,\n", + ")\n", + "mps_summary = py.symmray_mps_summary(psi_mps)\n", + "\n", + "print(\"MPS sites:\", psi_mps.num_sites)\n", + "print(\"MPS total U1U1 charge:\", psi_mps.overall_charge())\n", + "print(\"MPS Hamiltonian terms:\", len(ham_mps.terms))\n", + "print(\"MPS max bond dim:\", mps_summary[\"max_bond_dim\"])\n", + "print(\"first snake sites:\", snake_sites[:6])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Inspect One Fermionic Hubbard Term\n", + "\n", + "This checks that the Hamiltonian terms are Symmray block-sparse arrays, not qubit Pauli operators." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "edge, term = next(iter(ham.terms.items()))\n", + "term_summary = py.symmray_block_summary(term)\n", + "\n", + "print(\"example edge:\", edge)\n", + "print(\"term type:\", type(term))\n", + "print(\"term shape:\", term_summary[\"shape\"])\n", + "print(\"number of stored blocks:\", term_summary[\"num_blocks\"])\n", + "print(\"stored / dense entries:\", term_summary[\"stored_size\"], \"/\", term_summary[\"dense_size\"])\n", + "term_summary[\"blocks\"][:8]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Optional Energy Smoke Test\n", + "\n", + "For a 4x4 PEPS this exact local-energy contraction can still be nontrivial. Keep `COMPUTE_ENERGY = False` while editing quickly, then turn it on for the first numerical check." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "COMPUTE_ENERGY = False\n", + "\n", + "if COMPUTE_ENERGY:\n", + " e0 = psi.energy_density(ham)\n", + " print(\"initial /N:\", np.real_if_close(e0))\n", + " print(\"initial /N - U/4:\", np.real_if_close(e0 - U / 4))\n", + "else:\n", + " print(\"Set COMPUTE_ENERGY = True to contract the initial energy density.\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Optional Tiny Imaginary-Time Step\n", + "\n", + "The paper uses optimized preparation circuits before probing eta pairing. As a direct-fermion TN starting point, a small imaginary-time step is a simpler first sanity check. It is not meant to reproduce the paper's optimized state yet." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "RUN_ONE_IMAGINARY_TIME_STEP = False\n", + "\n", + "if RUN_ONE_IMAGINARY_TIME_STEP:\n", + " psi_it = psi.ground_state(\n", + " dt=0.01,\n", + " steps=1,\n", + " hamiltonian=ham,\n", + " order=1,\n", + " max_bond=4,\n", + " cutoff=1.0e-10,\n", + " method=\"gate\",\n", + " inplace=False,\n", + " )\n", + " print(\"after one step max bond:\", psi_it.tn.max_bond())\n", + " print(\"after one step norm:\", np.real_if_close(psi_it.norm()))\n", + "else:\n", + " psi_it = None\n", + " print(\"Set RUN_ONE_IMAGINARY_TIME_STEP = True to run one tiny projection step.\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Paper Light-Pulse Constants\n", + "\n", + "These constants are included here so the next notebook iteration can add Peierls phases directly to fermionic hopping terms, still without any fermion-to-qubit mapping." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "omega = 4 * math.pi / 3\n", + "tau = math.pi / (2 * omega)\n", + "\n", + "def vector_potential(s):\n", + " return math.pi * (1 - math.cos(omega * s)) / 2\n", + "\n", + "print(\"omega:\", omega)\n", + "print(\"tau:\", tau)\n", + "print(\"pulse end time pi/omega:\", math.pi / omega)\n", + "print(\"A(0):\", vector_potential(0.0))\n", + "print(\"A(pi/omega):\", vector_potential(math.pi / omega))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Doped Checkerboard Targets From The Paper\n", + "\n", + "For the next direct-fermion notebook, the doped checkerboard case should use a 6x6 lattice with spin-resolved charge `(15, 15)`. That is 30 fermions on 36 sites, i.e. six holes relative to half filling. The paper gives two useful comparison numbers at the weak-coupling checkerboard point `t' = 0`, `U/t = 2`: exact `/N - U/4 = -1.27367`, and theoretical d-wave average `0.108`. The experimental d-wave average is `0.079 +/- 0.005`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "doped_Lx = 6\n", + "doped_Ly = 6\n", + "doped_sites = doped_Lx * doped_Ly\n", + "doped_target_charge = (15, 15)\n", + "doped_holes = doped_sites - sum(doped_target_charge)\n", + "checkerboard_reference_energy_shifted = -1.27367\n", + "checkerboard_theory_dwave_average = 0.108\n", + "checkerboard_experiment_dwave_average = (0.079, 0.005)\n", + "\n", + "print(\"doped checkerboard target charge:\", doped_target_charge)\n", + "print(\"holes relative to half filling:\", doped_holes)\n", + "print(\"reference /N - U/4:\", checkerboard_reference_energy_shifted)\n", + "print(\"theory d-wave average:\", checkerboard_theory_dwave_average)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Next Direct-Fermion Steps\n", + "\n", + "- Add weighted or phased Fermi-Hubbard terms so the Peierls light pulse acts directly on hopping edges.\n", + "- Add onsite eta-pair correlators `Delta_i^dag Delta_j + h.c.` as fermionic Symmray observables.\n", + "- Repeat the 4x4 run with `CYCLIC = True`, then compare against a small exact reference.\n", + "- Only after the 4x4 checks are stable, move the same workflow toward the paper's 6x6 half-filled setting." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/pepsy_examples/plan.md b/examples/pepsy_examples/plan.md new file mode 100644 index 0000000..4ce2ee3 --- /dev/null +++ b/examples/pepsy_examples/plan.md @@ -0,0 +1,272 @@ +# Pepsy/Symmray Fermi-Hubbard Simulation Plan + +## Goal + +Build Pepsy examples that use fermionic Symmray tensor networks for direct +Fermi-Hubbard simulations. The main methods reference for the Pepsy/Symmray +path is Gao et al., "Fermionic tensor network contraction for arbitrary +geometries", Phys. Rev. Research 7, 023193 (2025). + +Use arXiv:2511.02125, "Superconducting pairing correlations on a trapped-ion +quantum computer", as the physics benchmark source for square-lattice, +checkerboard, and bilayer settings. + +The examples should prioritize reproducible classical checks over reproducing +the trapped-ion circuit encoding. Use Pepsy's native fermionic tensor-network +path whenever possible so that the examples test the physics directly, without +the Octagon fermion-to-qubit overhead used in the hardware experiment. + +Paper sources: + +- Main direct-fermion TN reference: + https://doi.org/10.1103/PhysRevResearch.7.023193 +- Open arXiv text for the methods reference: + https://arxiv.org/abs/2410.02215 +- Physics benchmark paper: +- https://arxiv.org/pdf/2511.02125 +- Version checked while writing this plan: arXiv v3, dated 2026-02-18. + +## Implementation Notes From PRR 7, 023193 + +- Keep the fermions native: parity, Abelian charge sectors, and leg-order + metadata should stay with Symmray fermionic arrays rather than being encoded + through Jordan-Wigner, compact, or Octagon qubit mappings. +- Prefer graph-level geometry metadata. The PRR paper supports both globally + ordered and locally ordered conventions, with local ordering most natural for + arbitrary graphs; Pepsy examples should record site maps, edge orientations, + and any MPS snake embedding explicitly. +- Contraction planning can remain graph based. quimb/cotengra should choose + exact or approximate contraction trees from the tensor-network graph, while + Symmray handles fermionic swap/parity signs during transposition, + contraction, and decomposition. +- The PRR benchmarks use fermionic PEPS for half-filled Hubbard models on 3D + diamond lattices and random 3-regular graphs. Those are not the immediate + square-lattice targets below, but they are useful future regression tests for + Pepsy arbitrary-graph support. +- The paper's cluster-energy idea suggests an eventual Pepsy smoke benchmark: + compare full approximate contraction against small-radius graph clusters + with simple-update gauges on a fixed fermionic PEPS. + +## Existing Pepsy Pieces To Reuse + +- `pepsy.SymMPS` and `pepsy.SymPEPS` for Symmray-backed MPS/PEPS states. +- `pepsy.SymHamiltonian.from_edges("fermi_hubbard", ...)` for local + Symmray Fermi-Hubbard terms. +- `pepsy.default_physical_sectors(model="fermi_hubbard_u1u1")`, which gives + spin-resolved spinful local occupation sectors + `{(0, 0): 1, (0, 1): 1, (1, 0): 1, (1, 1): 1}`. +- `pepsy.site_charge_from_occupations(...)` to fix `U1U1` sectors such as + `(N_up, N_down)`. +- `state.time_evolve(...)`, `state.apply_gates(...)`, and + `state.measure(...)` for gate application and local-observable checks. +- `pepsy.symmray_mps_summary(...)` and `pepsy.symmray_peps_summary(...)` for + block-size diagnostics. + +Use spin-resolved `U1U1` for the paper-style direct-fermion Hubbard examples. +The older `model="fermi_hubbard"` preset remains useful when only total +particle-number `U1` is desired. + +## Target Paper Settings + +### 1. Half-filled square-lattice Hubbard and light-induced eta pairing + +Paper setting: + +- Lattice: periodic 6 x 6 square lattice. +- Model: nearest-neighbor spinful Fermi-Hubbard, `t = 1`. +- Direct-fermion symmetry target: spin-resolved `U1U1` charge `(18, 18)` for + the paper's 6 x 6 half-filled setting. +- Boundary treatment: average over periodic/anti-periodic choices in x and y + where feasible; otherwise start with periodic and document the difference. +- Benchmark limit: `U = 0`, 36 fermions, track imbalance `I_A = n_A - n_Abar` + under hopping evolution with Trotter step `tau = 0.5`. +- Interacting low-energy setting: `U/t = 8`. +- Best shallow state-preparation depth in the paper: + `(D_Heisenberg, D_Hubbard) = (1, 1)`. +- Paper reference values for sanity checks: + energy density `/N - U/4` is about `-2.36` raw, `-2.43` mitigated, and + periodic-boundary ground-state benchmark is about `-2.5278`. +- Light pulse: + `A(s) = pi * (1 - cos(omega * s)) / 2`, aligned with the lattice diagonal, + `omega = 4*pi/3 ~= U/2`. +- Pulse evolution: evolve to `pi/omega` using two second-order Trotter steps + with `tau = pi/(2*omega) = 0.375`. +- Relaxation check: apply two additional Hubbard Trotter steps at `tau = 0.375` + with the field off. +- Observables: + onsite doublon-pair correlations + `P_eta(x, y) = (1/N) sum_i `, + staggered average `P_eta_stag`, spin-spin correlations, energy density, and + the free-fermion imbalance benchmark. + +Example target: + +- `half_filled_eta_pairing.py` +- First make a 4 x 4 exact or high-bond reference. +- Then run 6 x 6 SymPEPS with a small sweep over PEPS bond dimension and + boundary/environment `chi`. +- Save compact JSON/CSV summaries: energy, norm drift, max bond, block density, + `P_eta(x, y)`, and `P_eta_stag`. + +### 2. Doped checkerboard Hubbard and d-wave bond pairing + +Paper setting: + +- Lattice: 6 x 6 square lattice split into 2 x 2 plaquettes. +- Doping: 1/6 hole doping. +- Direct-fermion symmetry target: spin-resolved `U1U1` charge `(15, 15)`, + i.e. 30 particles on 36 sites and six holes relative to half filling. +- Initial weak-coupling point: decoupled plaquettes, `t_prime = 0`, `U/t = 2`. +- Strong intra-plaquette hopping: `t`. +- Weak inter-plaquette hopping: `t_prime`. +- Effective preparation model: ferromagnetic 3 x 3 XXZ model with + anisotropy approximately `delta = -1` and `sum_i Z_i = -3`. +- Plaquette injection: + `|0> -> |s>` in the quarter-filled 2 x 2 plaquette sector, + `|1> -> |d>` in the half-filled 2 x 2 plaquette sector. +- Paper reference values: + exact weak-coupling energy density `/N - U/4 = -1.27367`; + d-wave average is about `0.108` in theory and `0.079 +/- 0.005` in the + experiment. +- Adiabatic target check: + one second-order Trotter step of size `tau = 0.3` toward + `U/t = 8`, `t_prime/t = 1/2`. +- Observables: + strong-bond singlet pair correlations + `P_b(x, y) = (4/N) sum_ `, + d-wave signed average, and energy density. + +Example target: + +- `checkerboard_dwave_pairing.py` +- Implement edge grouping for strong and weak bonds. +- Build exact 2 x 2 plaquette states `|s>` and `|d>` in dense form, then + convert to fermionic Symmray-compatible tensors or initialize a product PEPS + over plaquettes. +- Validate 2 x 2 plaquette energies before the 6 x 6 run. +- Run the weak-coupling state first, then the one-step ramp to the target + checkerboard parameters. + +### 3. Bilayer Hubbard nickelate-inspired s-wave rung pairing + +Paper setting: + +- Lattice: two coupled 4 x 4 square Hubbard layers, i.e. 4 x 4 x 2 sites. +- Filling: two quarter-filled layers in the bilayer construction. +- Model: + `H_bilayer = H_Hubbard^A + H_Hubbard^B + H_exchange`, + where + `H_exchange = J sum_i (S_iA . S_iB - n_iA n_iB / 4)`. +- Perturbative limit: `t/J -> 0`, `U >= 0`; ground subspace is rung singlets + and hole pairs. +- Effective preparation model: ferromagnetic 4 x 4 XXZ with + `delta = -2/3`. +- Rung injection: + `|0> -> |vac>` and + `|1> -> (c_Aup^dag c_Bdown^dag - c_Adown^dag c_Bup^dag) |vac> / sqrt(2)`. +- Target ramp: + `(t/J, U/J) = (0, 0) -> (0.7, 7)` with `J = 1`, so target + `t/J = 0.7`, `U/t = 10`. +- Trotter settings: `M = 1, 2` steps, `tau = 0.4`. +- Observables: + rung-rung singlet pairing + `P_r(x, y) = (1/N) sum_i `, + optional slanted-pair correlations, and target-Hamiltonian energy density. + +Example target: + +- `bilayer_rung_pairing.py` +- Represent bilayer sites as `(layer, x, y)` or flatten them with a reversible + map stored in the output metadata. +- Add bilayer Hamiltonian term construction for intra-layer hopping, onsite + `U`, and inter-layer exchange. +- Validate on 2 x 2 x 2 or 3 x 3 x 2 before running the paper's 4 x 4 x 2 + setting. + +## Implementation Phases + +1. Infrastructure and metadata + - Add a small `settings.py` with lattice maps, boundary-condition helpers, + output paths, and deterministic seeds. + - Add `observables.py` for number, spin, onsite eta-pair, bond-singlet, and + rung-singlet observables. + - Add `hamiltonians.py` for weighted Fermi-Hubbard and bilayer exchange + term dictionaries that can be wrapped in `pepsy.SymHamiltonian`. + +2. Small exact checks + - 2 x 2 spinful Hubbard plaquette at `U/t = 2`; verify `|s>` and `|d>`. + - 4 x 4 or smaller free-fermion imbalance at `U = 0`. + - 2 x 2 x 2 bilayer rung-singlet injection and exchange energy. + - Use dense NumPy/quimb exact diagonalization only at these small sizes. + +3. Symmray tensor-network checks + - Recreate the same small systems with `SymMPS` or `SymPEPS`. + - Assert that tensors and gates remain Symmray/FermionicArray-backed after + gate application. + - Track `overall_charge`, `max_bond`, block density, norm, and energy. + +4. Paper-size exploratory runs + - Half-filled 6 x 6: PEPS bond dimensions such as `D = 2, 4, 6` and + boundary `chi` sweeps. + - Checkerboard 6 x 6: weak-coupling product/plaquette state first, then + one Trotter ramp step. + - Bilayer 4 x 4 x 2: start with MPS snake-order references, then compare to + PEPS or layered PEPS if Pepsy supports the needed geometry cleanly. + +5. Result comparison + - Produce tables with paper reference values beside Pepsy/Symmray values. + - Plot or save `P_eta`, `P_b`, and `P_r` versus displacement. + - Report truncation sensitivity by `D`, `chi`, and cutoff, not just one + number. + +## Pepsy Gaps To Resolve Before Full Reproduction + +- Fermionic graph metadata: keep site ordering, edge orientation, graph + distance, and any local/global ordering convention explicit in outputs and + helpers, especially for non-square, bilayer, and snake-MPS embeddings. +- Weighted Hubbard edges: the checkerboard model needs per-edge `t` values, + and the light pulse needs time-dependent Peierls phases. +- Bilayer exchange: the current uniform Fermi-Hubbard builder is not enough for + `S_iA . S_iB - n_iA n_iB / 4`; add a local exchange-term builder. +- Pairing observables: add tested dense-to-Symmray constructors for charge-zero + two-site/four-site pair correlators, including fermionic signs. +- Geometry: decide whether bilayer examples should use a flattened MPS, + a custom PEPS-like graph through Symmray-from-edges, or a Pepsy lattice + wrapper extension. +- Boundary conditions: implement periodic/anti-periodic sign choices and + mixed-periodic averaging explicitly in metadata. + +## Validation Checklist + +- Run focused API and Symmray checks: + `pytest -q tests/test_symmetric_tensors.py tests/test_public_api.py` +- For new Hamiltonian/observable helpers, add example-local unit tests or + promote them to `tests/` if they become package APIs. +- For each example, write a small mode that finishes quickly: + `--size smoke`, `--bond-dim 2`, `--chi 8`, and fixed seeds. +- Store generated numeric outputs under an ignored results directory, not in + the notebook or source tree. +- Document any mismatch from the paper caused by boundary conditions, + finite-bond truncation, or using native fermionic TN instead of the + hardware compact encoding. + +## Suggested File Layout + +```text +examples/pepsy_examples/ + plan.md + Fermi_Hubbard/ + README.md + half_filled_4x4_direct_fermions.ipynb + fermihubbard_symmray/ # future script/module form + README.md + settings.py + hamiltonians.py + observables.py + half_filled_eta_pairing.py + checkerboard_dwave_pairing.py + bilayer_rung_pairing.py + tests/ + test_small_hubbard_checks.py + test_observables.py +``` diff --git a/history/2026-06-28-relative-svd-gesvd-fallback.md b/history/2026-06-28-relative-svd-gesvd-fallback.md new file mode 100644 index 0000000..85d6432 --- /dev/null +++ b/history/2026-06-28-relative-svd-gesvd-fallback.md @@ -0,0 +1,42 @@ +# 2026-06-28 — Relative torch SVD registration and `gesvd` fallback + +- Milestone: M6 — AD decomposition prototype / full-SVD hardening +- Branch / commit: `main` and `develop` / `acf207e` + +## What changed +- Added `pepsy.tensors.reg_rel_svd_torch()` as the preferred torch full-SVD + registration. +- Kept `reg_complex_svd_torch()` and `register_torch_linalg(mode="complex")` + as compatibility paths to the same relative-regularized rule. +- Documented the backward rule as Townsend rectangular full-SVD terms with + Lorentzian relative broadening and the complex SVD gauge correction. +- Added a CPU SciPy `gesvd` forward fallback for the relative/complex SVD path + when `torch.linalg.svd` fails, including batched real and complex inputs. + +## Why +- The first registration pass stabilized the backward rule but still used plain + CPU `torch.linalg.svd` in forward, so ill-conditioned high-chi bonds could + still crash on `gesdd` failures. +- The real-mode registration already had a forward fallback; the relative and + complex registration paths now have the same robustness. + +## How it was validated +- `pytest -q tests/test_core_seed.py::test_reg_rel_svd_torch_forward_falls_back_to_scipy_gesvd` -> `2 passed`. +- SVD-focused tests in `tests/test_core_seed.py` -> `9 passed`. +- `pytest -q tests/test_public_api.py tests/test_package_layout.py tests/test_core_seed.py` -> `470 passed`. +- `python -m pyflakes src/pepsy/backends/linalg_torch.py src/pepsy/backends/linalg.py tests/test_core_seed.py` -> passed. + +## Decisions / findings +- This remains a full-SVD autodiff shim. It is not the true truncated-SVD + pullback planned in `PLAN.md`. +- For real float64 forward-only runs, `register_torch_linalg(mode="real")` + remains appropriate. For complex/autodiff full-SVD paths, use + `reg_rel_svd_torch()` or `register_torch_linalg(mode="complex")`. + +## Next step (do this first next time) +- Keep true truncated-SVD AD work separate from this full-SVD robustness layer. + Add a dedicated split driver or custom VJP before changing Quimb truncation + behavior. + +## Open questions / blockers +- None for the full-SVD fallback. True truncated-SVD AD remains open. diff --git a/pyproject.toml b/pyproject.toml index 8f434ac..60e1c0d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,6 +13,7 @@ dependencies = [ "numpy", "quimb", "cotengra", + "cotengrust", "autoray", "tqdm", ] diff --git a/src/pepsy/__init__.py b/src/pepsy/__init__.py index 3d4f056..beec2ab 100644 --- a/src/pepsy/__init__.py +++ b/src/pepsy/__init__.py @@ -114,6 +114,9 @@ "build_optimizer": ".tensors", "contract_hypercompressed_tn": ".tensors", "default_physical_sectors": ".tensors", + "draw_symmray_blocks": ".tensors", + "draw_symmray_mps": ".tensors", + "draw_symmray_peps": ".tensors", "expec_mpo": ".tensors", "haar_random_state": ".tensors", "hrps_to_mps": ".tensors", @@ -132,6 +135,9 @@ "site_charge_from_map": ".tensors", "site_charge_from_occupations": ".tensors", "site_charge_uniform": ".tensors", + "symmray_block_summary": ".tensors", + "symmray_mps_summary": ".tensors", + "symmray_peps_summary": ".tensors", "symm_operator_from_dense": ".tensors", "tn_fidelity": ".tensors", "tn_norm": ".tensors", @@ -229,6 +235,9 @@ "SymMPS", "SymPEPS", "default_physical_sectors", + "draw_symmray_blocks", + "draw_symmray_mps", + "draw_symmray_peps", "expec_mpo", "haar_random_state", "hrps_to_mps", @@ -247,6 +256,9 @@ "site_charge_from_map", "site_charge_from_occupations", "site_charge_uniform", + "symmray_block_summary", + "symmray_mps_summary", + "symmray_peps_summary", "symm_operator_from_dense", "tn_fidelity", "tn_norm", @@ -345,6 +357,9 @@ def __getattr__(name): SymMPS, SymPEPS, default_physical_sectors, + draw_symmray_blocks, + draw_symmray_mps, + draw_symmray_peps, expec_mpo, haar_random_state, hrps_to_mps, @@ -363,6 +378,9 @@ def __getattr__(name): site_charge_from_map, site_charge_from_occupations, site_charge_uniform, + symmray_block_summary, + symmray_mps_summary, + symmray_peps_summary, symm_operator_from_dense, tn_fidelity, tn_norm, diff --git a/src/pepsy/backends/linalg.py b/src/pepsy/backends/linalg.py index ec595fa..2a615f8 100644 --- a/src/pepsy/backends/linalg.py +++ b/src/pepsy/backends/linalg.py @@ -6,6 +6,7 @@ import autoray as ar import jax.numpy as jnp +import numpy as np import scipy.linalg import torch from jax import custom_vjp @@ -62,6 +63,50 @@ def safe_inverse_2(x, eps): return x.clamp_min(eps).reciprocal() +def _scipy_gesvd(A): + """Compute a thin SVD through SciPy's more robust ``gesvd`` driver.""" + A_np = A.detach().cpu().numpy() + batch_shape = A_np.shape[:-2] + m, n = A_np.shape[-2:] + k = min(m, n) + + if batch_shape: + flat = A_np.reshape((-1,) + A_np.shape[-2:]) + if flat.shape[0] == 0: + U_np = np.empty(batch_shape + (m, k), dtype=A_np.dtype) + S_np = np.empty(batch_shape + (k,), dtype=A_np.real.dtype) + Vh_np = np.empty(batch_shape + (k, n), dtype=A_np.dtype) + else: + parts = [ + scipy.linalg.svd( + mat, + full_matrices=False, + lapack_driver="gesvd", + ) + for mat in flat + ] + U_np = np.stack([part[0] for part in parts]).reshape( + batch_shape + parts[0][0].shape + ) + S_np = np.stack([part[1] for part in parts]).reshape( + batch_shape + parts[0][1].shape + ) + Vh_np = np.stack([part[2] for part in parts]).reshape( + batch_shape + parts[0][2].shape + ) + else: + U_np, S_np, Vh_np = scipy.linalg.svd( + A_np, + full_matrices=False, + lapack_driver="gesvd", + ) + + U = torch.from_numpy(np.asarray(U_np)).to(device=A.device, dtype=A.dtype) + S = torch.from_numpy(np.asarray(S_np)).to(device=A.device, dtype=A.real.dtype) + Vh = torch.from_numpy(np.asarray(Vh_np)).to(device=A.device, dtype=A.dtype) + return U, S, Vh + + class SVD(torch.autograd.Function): """Custom torch SVD with stabilized backward for near-degenerate spectra. @@ -71,10 +116,13 @@ class SVD(torch.autograd.Function): @staticmethod def forward(ctx, A): - if A.is_cuda: - U, S, Vh = torch.linalg.svd(A, full_matrices=False, driver='gesvd') - else: - U, S, Vh = torch.linalg.svd(A, full_matrices=False) + try: + if A.is_cuda: + U, S, Vh = torch.linalg.svd(A, full_matrices=False, driver='gesvd') + else: + U, S, Vh = torch.linalg.svd(A, full_matrices=False) + except Exception: + U, S, Vh = _scipy_gesvd(A) # A = U @ diag(S) @ Vh ctx.save_for_backward(U, S, Vh) @@ -457,9 +505,13 @@ def reg_complex_svd_jax(): +def reg_rel_svd_torch(): + """Register the relative-regularized torch SVD rule in autoray.""" + ar.register_function('torch', 'linalg.svd', SVD.apply) + def reg_complex_svd_torch(): """Register the complex torch SVD autograd implementation in autoray.""" - ar.register_function('torch', 'linalg.svd', SVD.apply) + reg_rel_svd_torch() def reg_real_svd_torch(): """Register the real torch SVD autograd implementation in autoray.""" diff --git a/src/pepsy/backends/linalg_torch.py b/src/pepsy/backends/linalg_torch.py index 4787c25..0882b4a 100644 --- a/src/pepsy/backends/linalg_torch.py +++ b/src/pepsy/backends/linalg_torch.py @@ -1,6 +1,7 @@ """Torch-side linalg registrations with stabilized autodiff rules.""" import autoray as ar +import numpy as np import torch try: @@ -32,15 +33,75 @@ def safe_inverse_2(x, eps): return x.clamp_min(eps).reciprocal() +def _scipy_gesvd(A, exc): + """Compute a thin SVD through SciPy's more robust ``gesvd`` driver.""" + if scipy_linalg is None: + raise RuntimeError( + "torch.linalg.svd failed and scipy is unavailable for gesvd fallback." + ) from exc + + A_np = A.detach().cpu().numpy() + batch_shape = A_np.shape[:-2] + m, n = A_np.shape[-2:] + k = min(m, n) + + if batch_shape: + flat = A_np.reshape((-1,) + A_np.shape[-2:]) + if flat.shape[0] == 0: + U_np = np.empty(batch_shape + (m, k), dtype=A_np.dtype) + S_np = np.empty(batch_shape + (k,), dtype=A_np.real.dtype) + Vh_np = np.empty(batch_shape + (k, n), dtype=A_np.dtype) + else: + parts = [ + scipy_linalg.svd( + mat, + full_matrices=False, + lapack_driver="gesvd", + ) + for mat in flat + ] + U_np = np.stack([part[0] for part in parts]).reshape( + batch_shape + parts[0][0].shape + ) + S_np = np.stack([part[1] for part in parts]).reshape( + batch_shape + parts[0][1].shape + ) + Vh_np = np.stack([part[2] for part in parts]).reshape( + batch_shape + parts[0][2].shape + ) + else: + U_np, S_np, Vh_np = scipy_linalg.svd( + A_np, + full_matrices=False, + lapack_driver="gesvd", + ) + + U = torch.from_numpy(np.asarray(U_np)).to(device=A.device, dtype=A.dtype) + S = torch.from_numpy(np.asarray(S_np)).to(device=A.device, dtype=A.real.dtype) + Vh = torch.from_numpy(np.asarray(Vh_np)).to(device=A.device, dtype=A.dtype) + return U, S, Vh + + class SVD(torch.autograd.Function): - """Custom torch SVD with stabilized backward for near-degenerate spectra.""" + """Torch SVD with a relative-regularized reverse-mode rule. + + The rectangular real-SVD terms follow Townsend's reverse update, with the + singular-gap and inverse-singular-value reciprocals regularized as + ``x / (x**2 + eps)``. The scale-aware ``eps`` keeps the stabilizer relative + to the current singular spectrum, which is the Lorentzian broadening used in + differentiable tensor-network SVDs. Complex inputs additionally include the + phase/gauge term from the complex-valued SVD backward formula. + """ @staticmethod def forward(ctx, A): - if A.is_cuda: - U, S, Vh = torch.linalg.svd(A, full_matrices=False, driver="gesvd") - else: - U, S, Vh = torch.linalg.svd(A, full_matrices=False) + try: + if A.is_cuda: + U, S, Vh = torch.linalg.svd(A, full_matrices=False, driver="gesvd") + else: + U, S, Vh = torch.linalg.svd(A, full_matrices=False) + except Exception as exc: # pragma: no cover - backend failure dependent + U, S, Vh = _scipy_gesvd(A, exc) ctx.save_for_backward(U, S, Vh) return U, S, Vh @@ -85,6 +146,8 @@ def backward(ctx, gu, gsigma, gvh): eps_scale=sigma_scale, ) + # Townsend's F+/F- terms, written as 1/(s_j - s_i) and + # 1/(s_i + s_j), with relative Lorentzian broadening. F = sigma.unsqueeze(-2) - sigma.unsqueeze(-1) F = safe_inverse( F, @@ -137,6 +200,7 @@ def backward(ctx, gu, gsigma, gvh): dA = u_term + sigma_term + v_term if (u.is_complex() or v.is_complex()) and gu is not None: + # Complex-SVD gauge correction from arXiv:1909.02659. phase_diag = (uh @ gu).diagonal(0, -2, -1) L = 1j * phase_diag.imag * sigma_inv imag_term = (u * L.unsqueeze(-2)) @ vh @@ -287,11 +351,20 @@ def backward(ctx, dQ, dR): return dA -def reg_complex_svd_torch(): - """Register the complex torch SVD autograd implementation in autoray.""" +def reg_rel_svd_torch(): + """Register the relative-regularized torch SVD rule in autoray.""" ar.register_function("torch", "linalg.svd", SVD.apply) +def reg_complex_svd_torch(): + """Register the complex torch SVD autograd implementation in autoray. + + This compatibility name installs the same relative-regularized SVD rule as + :func:`reg_rel_svd_torch`. + """ + reg_rel_svd_torch() + + def reg_real_svd_torch(): """Register the real torch SVD autograd implementation in autoray.""" ar.register_function("torch", "linalg.svd", SVD_real.apply) diff --git a/src/pepsy/operators/gates.py b/src/pepsy/operators/gates.py index d40228d..45cecf3 100644 --- a/src/pepsy/operators/gates.py +++ b/src/pepsy/operators/gates.py @@ -1207,6 +1207,117 @@ def _rectangular_swap_gate(dim_a, dim_b, *, dtype="complex128"): return swap_gate +def _is_symmray_array(value): + """Return whether ``value`` looks like a Symmray block-sparse array.""" + return hasattr(value, "blocks") and hasattr(value, "indices") + + +def _symmray_sample_data_from_tn(tn): + """Return representative Symmray array data from a tensor network.""" + tensor_map = getattr(tn, "tensor_map", None) + if tensor_map: + tensors = tensor_map.values() + else: + try: + tensors = tuple(tn) + except TypeError: + tensors = () + + for tensor in tensors: + data = getattr(tensor, "data", None) + if _is_symmray_array(data): + return data + return None + + +def _symmray_index_map_for_tn_ind(tn, ind): + """Return the Symmray charge map for a live tensor-network index.""" + tensor_ids = getattr(tn, "ind_map", {}).get(ind, ()) + tensors = [] + tensor_map = getattr(tn, "tensor_map", {}) + for tid in tensor_ids: + try: + tensors.append(tensor_map[tid]) + except KeyError: + continue + if not tensors: + try: + tensors = tuple(tn) + except TypeError: + tensors = () + + for tensor in tensors: + data = getattr(tensor, "data", None) + if not _is_symmray_array(data) or ind not in getattr(tensor, "inds", ()): + continue + axis = tensor.inds.index(ind) + return dict(data.indices[axis].chargemap) + + raise ValueError(f"Could not infer Symmray charge map for index {ind!r}.") + + +def _symmray_dense_index_map_from_chargemap(chargemap): + """Expand ``{charge: size}`` into Symmray's ``{dense_index: charge}``.""" + index_map = {} + dense_index = 0 + for charge, size in dict(chargemap).items(): + for _ in range(int(size)): + index_map[dense_index] = charge + dense_index += 1 + return index_map + + +def _symmray_dense_gate_from_site_maps(tn, gate, output_sites, input_sites, ind_id): + """Convert a dense site gate to a Symmray array using live site sectors.""" + sample = _symmray_sample_data_from_tn(tn) + if sample is None: + return None + + index_maps = [] + for site in tuple(output_sites) + tuple(input_sites): + ind = _physical_index_for_site(tn, site, ind_id) + index_maps.append( + _symmray_dense_index_map_from_chargemap( + _symmray_index_map_for_tn_ind(tn, ind) + ) + ) + + nout = len(tuple(output_sites)) + nin = len(tuple(input_sites)) + array_cls = type(sample) + kwargs = {} + if array_cls.__name__ in {"AbelianArray", "FermionicArray"}: + symmetry = getattr(sample, "symmetry", None) + if symmetry is not None: + kwargs["symmetry"] = symmetry + + return array_cls.from_dense( + np.asarray(gate), + index_maps=tuple(index_maps), + duals=(False,) * nout + (True,) * nin, + charge=0, + **kwargs, + ) + + +def _symmray_swap_gate_for_site_pair(tn, site_a, site_b, *, ind_id=None, dtype="complex128"): + """Return a Symmray block-sparse SWAP for two live physical site legs.""" + sample = _symmray_sample_data_from_tn(tn) + if sample is None: + return None + + dim_a = _physical_index_size_for_site(tn, site_a, ind_id) + dim_b = _physical_index_size_for_site(tn, site_b, ind_id) + dense_swap = _rectangular_swap_gate(dim_a, dim_b, dtype=dtype) + return _symmray_dense_gate_from_site_maps( + tn, + dense_swap, + output_sites=(site_b, site_a), + input_sites=(site_a, site_b), + ind_id=ind_id, + ) + + def _convert_internal_gate_to_backend(gate, inferred_converter): """Best-effort conversion for internally generated exact gates.""" if inferred_converter is None: @@ -1227,6 +1338,16 @@ def _swap_gate_for_site_pair( inferred_converter=None, ): """Return a SWAP tensor matching the sites' current physical dimensions.""" + symmray_swap = _symmray_swap_gate_for_site_pair( + tn, + site_a, + site_b, + ind_id=ind_id, + dtype=dtype, + ) + if symmray_swap is not None: + return symmray_swap + dim_a = _physical_index_size_for_site(tn, site_a, ind_id) dim_b = _physical_index_size_for_site(tn, site_b, ind_id) swap_gate = _rectangular_swap_gate(dim_a, dim_b, dtype=dtype) @@ -1633,8 +1754,10 @@ def gate(tn, gates, where=None, which=None, **kwargs): when ``contract`` is ``"split"`` or ``"reduce-split"``. For other contract modes, the gate is applied directly to the requested endpoints. By default, 2D/3D routing uses ``sequence="auto"`` to choose a shortest - route with the smallest current virtual-bond bottleneck. Pass an explicit - deterministic sequence to force a particular route. + route with the smallest current virtual-bond bottleneck, and + ``contract="reduce-split"`` for quimb's reduced two-site split path. Pass + an explicit deterministic sequence or contract mode to force a particular + route or split strategy. The efficient routed pattern is:: gate( @@ -1730,7 +1853,7 @@ def gate(tn, gates, where=None, which=None, **kwargs): if arity == 2: opts_local = dict(opts) - opts_local.setdefault("contract", "split") + opts_local.setdefault("contract", "reduce-split") if which_payload is not None: opts_local["ind_id"] = _ind_id_from_which(which_payload, 2) elif (which_default is not None) and ("ind_id" not in opts_local): @@ -1748,7 +1871,7 @@ def gate(tn, gates, where=None, which=None, **kwargs): if arity == 3: opts_local = dict(opts) - opts_local.setdefault("contract", "split") + opts_local.setdefault("contract", "reduce-split") if which_payload is not None: opts_local["ind_id"] = _ind_id_from_which(which_payload, 3) elif (which_default is not None) and ("ind_id" not in opts_local): @@ -2203,7 +2326,7 @@ def _apply_gate_2d( bond_dim=None, max_bond=None, bra=False, - contract="split", + contract="reduce-split", tags=None, dtype="complex128", cutoff=1.0e-12, @@ -2478,7 +2601,7 @@ def _apply_gate_3d( bond_dim=None, max_bond=None, bra=False, - contract="split", + contract="reduce-split", tags=None, dtype="complex128", cutoff=1.0e-12, @@ -2692,7 +2815,7 @@ def build_pepo_from_gates( dtype="complex128", max_bond=16, sequence="auto", - contract="split", + contract="reduce-split", ind_id="k{},{}", ): """Build a PEPO from gate-style input on top of a PEPO identity. @@ -2725,8 +2848,9 @@ def build_pepo_from_gates( 2D SWAP-path preference for long-range two-site gates. Defaults to ``"auto"`` for the same lower-bond smart routing used by :func:`gate`. - contract : str, default="split" - Gate contraction mode. + contract : str, default="reduce-split" + Gate contraction mode. The default uses quimb's reduced two-site split + path, which is usually cheaper than ``"split"`` for PEPO/PEPS tensors. ind_id : str, default="k{},{}" Physical index format used for PEPO ket-family indices. @@ -2806,7 +2930,7 @@ def build_mpo_from_gates( mpo_=None, dtype="complex128", max_bond=16, - contract="split", + contract="reduce-split", ind_id="k{}", ): """Build an MPO from gate-style input on top of an MPO identity. @@ -2835,8 +2959,9 @@ def build_mpo_from_gates( max_bond : int, default=16 Per-gate local split truncation cap and fallback construction compression cap. - contract : str, default="split" - Gate contraction mode. + contract : str, default="reduce-split" + Gate contraction mode. The default uses quimb's reduced two-site split + path, which is usually cheaper than ``"split"`` for MPO tensors. ind_id : str, default="k{}" Physical index format used for MPO ket-family indices. diff --git a/src/pepsy/operators/hamiltonians.py b/src/pepsy/operators/hamiltonians.py index 2d38009..c9b4020 100644 --- a/src/pepsy/operators/hamiltonians.py +++ b/src/pepsy/operators/hamiltonians.py @@ -1178,17 +1178,17 @@ def _show_mpo_schematic_2d( ) try: - from quimb import schematic + from matplotlib import colormaps except ImportError as exc: # pragma: no cover - optional dependency raise ImportError( - "Schematic plotting requires quimb.schematic to be available." + "Schematic plotting requires matplotlib to be available." ) from exc try: - from matplotlib import colormaps + from quimb import schematic except ImportError as exc: # pragma: no cover - optional dependency raise ImportError( - "Schematic plotting requires matplotlib to be available." + "Schematic plotting requires quimb.schematic to be available." ) from exc positions = ( diff --git a/src/pepsy/optimizers/mps/optimizer.py b/src/pepsy/optimizers/mps/optimizer.py index e5e4840..28c3d7c 100644 --- a/src/pepsy/optimizers/mps/optimizer.py +++ b/src/pepsy/optimizers/mps/optimizer.py @@ -194,6 +194,82 @@ def _infer_gate_dims(gate, where): return None return dims_in + @staticmethod + def _is_symmray_array(value): + """Return whether ``value`` looks like a Symmray block-sparse array.""" + return hasattr(value, "blocks") and hasattr(value, "indices") + + @classmethod + def _has_symmray_data(cls, tn): + """Return whether any tensor data in ``tn`` is Symmray-backed.""" + return any( + cls._is_symmray_array(tensor.data) + for tensor in getattr(tn, "tensors", ()) + ) + + @staticmethod + def _is_nearest_neighbor_1d(where): + """Return whether an integer two-site location is adjacent in MPS order.""" + if len(where) != 2: + return True + site0, site1 = where + if not isinstance(site0, Integral) or not isinstance(site1, Integral): + return True + return abs(int(site0) - int(site1)) == 1 + + def _validate_symmray_mode_support(self): + """Fail early for Symmray/MPS mode combinations with known bad paths.""" + if not self._has_symmray_data(self.p): + return + + if self.mode == "dmrg" and self.p.max_bond() < self.chi: + raise ValueError( + "Symmray MPS data cannot be automatically expanded for " + "mode='dmrg' because quimb's bond-dimension expansion calls " + "dense-style pad on the Symmray backend. Construct the SymMPS " + "with bond_dim >= chi, or run with chi <= p.max_bond()." + ) + + def _apply_symmray_auto_swap_gate( + self, + p, + gate, + where, + *, + cutoff, + cutoff_mode, + max_bond=None, + info=None, + ): + """Apply a Symmray two-site gate through quimb's block-aware swaps.""" + compress_opts = { + "cutoff": cutoff, + "cutoff_mode": cutoff_mode, + } + if max_bond is not None: + compress_opts["max_bond"] = max_bond + p.gate_with_auto_swap_( + gate, + where, + info=self.info_c if info is None else info, + swap_back=True, + **compress_opts, + ) + return p + + def _build_symmray_auto_swap_target(self, p, gate, where, cutoff, cutoff_mode): + """Build an un-chi-capped target using Symmray-aware swap routing.""" + p_target = p.copy() + self._apply_symmray_auto_swap_gate( + p_target, + gate, + where, + cutoff=cutoff, + cutoff_mode=cutoff_mode, + info={}, + ) + return p_target + def _apply_gate(self, p, gate, where, **kwargs): """Apply a gate using this optimizer's physical-index convention.""" kwargs.setdefault("ind_id", self.ind_id) @@ -366,6 +442,7 @@ def run( # pylint: disable=too-many-arguments,too-many-positional-arguments where_seq = list(self.where) if not G_seq: return self.p + self._validate_symmray_mode_support() non_unitary = bool(non_unitary) if not non_unitary: @@ -1212,32 +1289,55 @@ def _run_mpo( # pylint: disable=too-many-locals raise ValueError("Each gate location must have one or two sites.") two_qubit_count += 1 xmin, xmax = sorted(where) + use_symmray_auto_swap = ( + self._has_symmray_data(p) + and not self._is_nearest_neighbor_1d(where) + ) if track_norm_infidelity: self.canonize_mps(p, (xmin, xmax)) if track_norm_infidelity or track_infidelity: - p_target = self._build_norm_target( - p, - gate, - where, - cutoff, - cutoff_mode, - ) + if use_symmray_auto_swap: + p_target = self._build_symmray_auto_swap_target( + p, + gate, + where, + cutoff, + cutoff_mode, + ) + else: + p_target = self._build_norm_target( + p, + gate, + where, + cutoff, + cutoff_mode, + ) else: p_target = None if track_norm_infidelity: target_norm = self._canonical_span_norm(p_target, (xmin, xmax)) else: target_norm = None - p.gate_nonlocal_( - gate, - where, - dims=self._infer_gate_dims(gate, where), - max_bond=self.chi, - info=self.info_c, - method="direct", - cutoff=cutoff, - cutoff_mode=cutoff_mode, - ) + if use_symmray_auto_swap: + self._apply_symmray_auto_swap_gate( + p, + gate, + where, + cutoff=cutoff, + cutoff_mode=cutoff_mode, + max_bond=self.chi, + ) + else: + p.gate_nonlocal_( + gate, + where, + dims=self._infer_gate_dims(gate, where), + max_bond=self.chi, + info=self.info_c, + method="direct", + cutoff=cutoff, + cutoff_mode=cutoff_mode, + ) idx += 1 advanced = 1 last_where = (xmin, xmax) @@ -1478,7 +1578,9 @@ def _run_svd( # pylint: disable=too-many-locals """Apply gates with local SVD compression for nonlocal 2-site gates. Two-site gates are applied with ``contract="reduce-split"`` then - compressed on the local span to ``max_bond=self.chi``. + compressed on the local span to ``max_bond=self.chi``. Symmray-backed + MPS data use quimb's block-aware auto-swap split path instead, since + ``reduce-split`` loses Symmray fusion metadata in current quimb/Symmray. """ p = self.p two_qubit_count = 0 @@ -1532,35 +1634,61 @@ def _run_svd( # pylint: disable=too-many-locals raise ValueError("Each gate location must have one or two sites.") two_qubit_count += 1 - compress_opts = {"cutoff": cutoff} + compress_opts = {"cutoff": cutoff, "cutoff_mode": cutoff_mode} xmin, xmax = sorted(where) + use_symmray_auto_swap = self._has_symmray_data(p) if track_norm_infidelity: self.canonize_mps(p, (xmin, xmax)) - self._apply_gate( - p, - gate, - where, - contract="reduce-split", - cutoff=cutoff, - cutoff_mode=cutoff_mode, - inplace=True, - ) - target_norm = ( - self._canonical_span_norm(p, (xmin, xmax)) - if track_norm_infidelity - else None - ) - p_target = p.copy() if track_infidelity else None - self.canonize_mps(p, (xmin, xmax)) - - for i in range(xmax, xmin, -1): - p.right_canonize_site(i, bra=None) - p.left_compress( - start=xmin, - stop=xmax, - max_bond=self.chi, - **compress_opts, - ) + if use_symmray_auto_swap: + if track_norm_infidelity or track_infidelity: + p_target = self._build_symmray_auto_swap_target( + p, + gate, + where, + cutoff, + cutoff_mode, + ) + else: + p_target = None + target_norm = ( + self._canonical_span_norm(p_target, (xmin, xmax)) + if track_norm_infidelity + else None + ) + self._apply_symmray_auto_swap_gate( + p, + gate, + where, + cutoff=cutoff, + cutoff_mode=cutoff_mode, + max_bond=self.chi, + ) + else: + self._apply_gate( + p, + gate, + where, + contract="reduce-split", + cutoff=cutoff, + cutoff_mode=cutoff_mode, + inplace=True, + ) + target_norm = ( + self._canonical_span_norm(p, (xmin, xmax)) + if track_norm_infidelity + else None + ) + p_target = p.copy() if track_infidelity else None + self.canonize_mps(p, (xmin, xmax)) + + for i in range(xmax, xmin, -1): + p.right_canonize_site(i, bra=None) + p.left_compress( + start=xmin, + stop=xmax, + max_bond=self.chi, + **compress_opts, + ) idx += 1 advanced = 1 diff --git a/src/pepsy/optimizers/sweep/optimizer.py b/src/pepsy/optimizers/sweep/optimizer.py index 3ebc35b..c930c4c 100644 --- a/src/pepsy/optimizers/sweep/optimizer.py +++ b/src/pepsy/optimizers/sweep/optimizer.py @@ -569,8 +569,11 @@ def _maybe_store_best_state(self, loss_value): if loss_value < 0.0: return if loss_value < float(self.best_loss): + state = getattr(self, "state", None) + if state is None: + return self.best_loss = loss_value - self.best_state = self.state.copy() + self.best_state = state.copy() def _ensure_boundary_chi(self, chi): """Retune both stored boundary objects to at least ``chi``. diff --git a/src/pepsy/tensors/README.md b/src/pepsy/tensors/README.md index 3a6e7a1..bb789cd 100644 --- a/src/pepsy/tensors/README.md +++ b/src/pepsy/tensors/README.md @@ -47,7 +47,10 @@ Backend helpers manage package-wide defaults and optional linalg shims: - `set_default_array_backend(...)` / `get_default_array_backend()` - `set_default_grad_backend(...)` / `get_default_grad_backend()` - `reset_default_backends()` -- torch and JAX complex-SVD/stop-gradient registrations +- torch and JAX linalg/stop-gradient registrations. For torch SVD, + `reg_rel_svd_torch()` is the preferred full-SVD autodiff shim; it installs + the relative-regularized backward rule also used by `reg_complex_svd_torch()` + and falls back to SciPy `gesvd` on CPU forward-driver failures. ## Tag and index conventions @@ -67,6 +70,25 @@ these conventions for shape inference and layer construction. helpers. Symmray remains optional. Code and tests that depend on it should import lazily or use `pytest.importorskip("symmray")`. +For spinful Fermi-Hubbard states, the named model presets are: + +- `fermi_hubbard`: total particle-number `U1` sectors. +- `fermi_hubbard_u1u1`: spin-resolved `U1U1` sectors with charges + `(N_up, N_down)`. + +Use Gao et al., "Fermionic tensor network contraction for arbitrary +geometries", Phys. Rev. Research 7, 023193 (2025), +https://doi.org/10.1103/PhysRevResearch.7.023193 as the main methods +reference for Pepsy/Symmray Fermi-Hubbard examples. The relevant design cue is +to preserve Symmray fermionic parity, symmetry, and leg-order metadata through +gate application, measurement, boundary contraction, and any future arbitrary +graph lattice wrappers. + +`SymMPS.fermionic_ordering()` and `SymPEPS.fermionic_ordering()` expose the +package-level record of site order, edge order, local index directions, and the +methods reference. The same record is included in `symmray_mps_summary(...)` +and `symmray_peps_summary(...)` under the `fermionic_ordering` key. + When changing symmetric behavior, check both dense compatibility and Symmray routing through PEPS optimizers and boundary contraction paths. diff --git a/src/pepsy/tensors/__init__.py b/src/pepsy/tensors/__init__.py index ac05765..e44f5d8 100644 --- a/src/pepsy/tensors/__init__.py +++ b/src/pepsy/tensors/__init__.py @@ -9,11 +9,17 @@ SymMPS, SymPEPS, default_physical_sectors, + draw_symmray_blocks, + draw_symmray_mps, + draw_symmray_peps, sector_index_map, site_charge_alternating, site_charge_from_map, site_charge_from_occupations, site_charge_uniform, + symmray_block_summary, + symmray_mps_summary, + symmray_peps_summary, symm_operator_from_dense, ) from .core import ( @@ -43,6 +49,7 @@ random_haar_qubit, reg_complex_svd_jax, reg_complex_svd_torch, + reg_rel_svd_torch, reg_stop_gradient_torch, register_torch_linalg, reset_default_backends, @@ -70,6 +77,9 @@ "contract_hypercompressed_tn", "expec_mpo", "default_physical_sectors", + "draw_symmray_blocks", + "draw_symmray_mps", + "draw_symmray_peps", "get_default_array_backend", "get_default_grad_backend", "haar_random_state", @@ -86,6 +96,7 @@ "random_haar_qubit", "reg_complex_svd_jax", "reg_complex_svd_torch", + "reg_rel_svd_torch", "reg_stop_gradient_torch", "register_torch_linalg", "reset_default_backends", @@ -97,6 +108,9 @@ "site_charge_from_occupations", "site_charge_uniform", "stop_grad", + "symmray_block_summary", + "symmray_mps_summary", + "symmray_peps_summary", "symm_operator_from_dense", "tn_fidelity", "tn_norm", diff --git a/src/pepsy/tensors/core.py b/src/pepsy/tensors/core.py index 2a0c2f3..23e7081 100644 --- a/src/pepsy/tensors/core.py +++ b/src/pepsy/tensors/core.py @@ -23,6 +23,7 @@ "backend_cupy", "backend_jax", "register_torch_linalg", + "reg_rel_svd_torch", "reg_complex_svd_torch", "reg_complex_svd_jax", "reg_stop_gradient_torch", @@ -486,17 +487,17 @@ def _show_schematic_2d( figsize=None, ): try: - from quimb import schematic + from matplotlib import colormaps except ImportError as exc: # pragma: no cover - optional dependency raise ImportError( - "Schematic plotting requires quimb.schematic to be available." + "Schematic plotting requires matplotlib to be available." ) from exc try: - from matplotlib import colormaps + from quimb import schematic except ImportError as exc: # pragma: no cover - optional dependency raise ImportError( - "Schematic plotting requires matplotlib to be available." + "Schematic plotting requires quimb.schematic to be available." ) from exc if title is None: @@ -910,7 +911,7 @@ def register_torch_linalg(mode="complex"): from ..backends import linalg_torch as lr # pylint: disable=import-outside-toplevel if mode == "complex": - lr.reg_complex_svd_torch() + lr.reg_rel_svd_torch() lr.reg_complex_qr_torch() return if mode == "real": @@ -920,8 +921,29 @@ def register_torch_linalg(mode="complex"): raise ValueError("mode must be 'complex' or 'real'") +def reg_rel_svd_torch(): + """Register torch SVD with a stable relative-regularized backward rule. + + The registered autoray ``torch`` SVD uses Townsend's rectangular SVD + reverse-mode update, Lorentzian broadening of singular-value denominators + from differentiable tensor-network practice, and the complex phase/gauge + correction for complex-valued SVDs. + """ + if torch is None: # pragma: no cover - exercised in no-torch CI + raise ImportError( + "reg_rel_svd_torch requires optional dependency 'torch'. " + "Install it with: pip install pepsy[torch] (or pip install torch)." + ) + from ..backends import linalg_torch as lr # pylint: disable=import-outside-toplevel + + lr.reg_rel_svd_torch() + + def reg_complex_svd_torch(): - """Register complex torch SVD autograd rule in autoray.""" + """Register complex torch SVD autograd rule in autoray. + + Compatibility wrapper for :func:`reg_rel_svd_torch`. + """ if torch is None: # pragma: no cover - exercised in no-torch CI raise ImportError( "reg_complex_svd_torch requires optional dependency 'torch'. " @@ -929,7 +951,7 @@ def reg_complex_svd_torch(): ) from ..backends import linalg_torch as lr # pylint: disable=import-outside-toplevel - lr.reg_complex_svd_torch() + lr.reg_rel_svd_torch() def reg_complex_svd_jax(): @@ -976,10 +998,22 @@ def stop_grad(x): return x +def _ensure_cotengrust(): + """Import cotengrust so cotengra can use accelerated pathfinders.""" + try: + import cotengrust as ctgr # pylint: disable=import-outside-toplevel + except ImportError as exc: # pragma: no cover - exercised by packaging/env failures + raise ImportError( + "Pepsy requires 'cotengrust' for accelerated cotengra path search. " + "Install it with: pip install cotengrust" + ) from exc + return ctgr + + def build_optimizer( progbar: bool = False, alpha: int = 64, - max_time="rate:7e8", + max_time="rate:7e8", max_repeats: int = 2**8, parallel="auto", optlib: str = "cmaes", @@ -993,6 +1027,10 @@ def build_optimizer( ): """Build a reusable cotengra contraction optimizer. + ``cotengrust`` is imported up front so cotengra's ``accel="auto"`` + pathfinders can use the Rust implementations when constructing greedy, + random-greedy, optimal, and reconfiguration paths. + Parameters ---------- progbar : bool, optional @@ -1022,6 +1060,8 @@ def build_optimizer( slicing_reconf_opts : dict | None, optional Options for interleaved slicing and reconfiguration. """ + _ensure_cotengrust() + # cotengra expects directory to be str, True, or None — not False. if directory is False: directory = None @@ -1050,6 +1090,7 @@ def build_optimizer( return ctg.ReusableHyperOptimizer(**kwargs) + def build_compressed_optimizer( progbar=True, chi=4, @@ -1059,12 +1100,17 @@ def build_compressed_optimizer( ): """Build and return a reusable cotengra compressed optimizer. + ``cotengrust`` is imported up front so cotengra can use accelerated + contraction-ordering primitives in any supported compressed path searches. + Parameters ---------- directory : None, True, or str, optional Passed directly to cotengra. ``None`` disables caching; ``True`` auto-generates a directory in the current working directory. """ + _ensure_cotengrust() + copt = ctg.ReusableHyperCompressedOptimizer( chi, max_repeats=max_repeats, @@ -1639,17 +1685,19 @@ def expec_mpo(mpo, mps, *, contraction_opt=None): mps_n = mps.copy() norm_ = mps_n.normalize() L = mps.L + divisor = 1.0 else: mps_n = mps.copy() L = len(mps.outer_inds()) norm_ = tn_norm(mps_n, contraction_opt=contraction_opt) - + divisor = norm_ + if norm_ == 0.0: raise ValueError("Cannot compute normalized expectation for a zero-norm state.") mps_h = mps_n.H mps_h.reindex_({f"k{i}": f"b{i}" for i in range(L)}) - return (mps_h | mpo | mps_n).contract(all, optimize=contraction_opt) / norm_ + return (mps_h | mpo | mps_n).contract(all, optimize=contraction_opt) / divisor def ps_to_peps(Lx: int, Ly: int, dtype: str = "complex128", theta: float = 0.0, cyclic: bool = False, chi: int = 1, rand_strength: float = 0.0): diff --git a/src/pepsy/tensors/symmetric.py b/src/pepsy/tensors/symmetric.py index a45982c..5fb7adf 100644 --- a/src/pepsy/tensors/symmetric.py +++ b/src/pepsy/tensors/symmetric.py @@ -12,16 +12,28 @@ __all__ = ["SymGateStream", "SymHamiltonian", "SymMPS", "SymPEPS"] __all__ += [ "default_physical_sectors", + "draw_symmray_blocks", + "draw_symmray_mps", + "draw_symmray_peps", "sector_index_map", "site_charge_alternating", "site_charge_from_map", "site_charge_from_occupations", "site_charge_uniform", + "symmray_block_summary", + "symmray_mps_summary", + "symmray_peps_summary", "symm_operator_from_dense", ] _SYMMRAY_AUTORAY_REGISTERED = False +_FERMIONIC_TN_METHODS_REFERENCE = { + "title": "Fermionic tensor network contraction for arbitrary geometries", + "doi": "10.1103/PhysRevResearch.7.023193", + "url": "https://doi.org/10.1103/PhysRevResearch.7.023193", +} + def _require_symmray(): """Import symmray with a clear optional-dependency message.""" @@ -73,6 +85,10 @@ def _allclose(a, b, rtol=1e-5, atol=1e-8, **kwargs): "fermi-hubbard": "fermi_hubbard", "hubbard": "fermi_hubbard", "fh": "fermi_hubbard", + "fermi_hubbard_u1u1": "fermi_hubbard_u1u1", + "fermi-hubbard-u1u1": "fermi_hubbard_u1u1", + "hubbard_u1u1": "fermi_hubbard_u1u1", + "fh_u1u1": "fermi_hubbard_u1u1", "spinless_fermi_hubbard": "fermi_hubbard_spinless", "spinless-fermi-hubbard": "fermi_hubbard_spinless", "fermi_hubbard_spinless": "fermi_hubbard_spinless", @@ -85,6 +101,7 @@ def _allclose(a, b, rtol=1e-5, atol=1e-8, **kwargs): "tfim": {"symmetry": "Z2", "fermionic": False, "phys_dim": 2}, "heisenberg": {"symmetry": "U1", "fermionic": False, "phys_dim": 2}, "fermi_hubbard": {"symmetry": "U1", "fermionic": True, "phys_dim": 4}, + "fermi_hubbard_u1u1": {"symmetry": "U1U1", "fermionic": True, "phys_dim": 4}, "fermi_hubbard_spinless": {"symmetry": "U1", "fermionic": True, "phys_dim": 2}, } @@ -133,6 +150,1603 @@ def sector_index_map(sectors): return out +def _require_symmray_array(value, *, name="value"): + if _is_symmray_array(getattr(value, "data", None)): + return value.data + if not _is_symmray_array(value): + raise TypeError(f"{name} must be a Symmray block-sparse array.") + return value + + +def _mapping_items(mapping): + return sorted(mapping.items(), key=lambda item: repr(item[0])) + + +def _as_tuple(value): + if isinstance(value, tuple): + return value + if isinstance(value, list): + return tuple(value) + return (value,) + + +def _format_compact_mapping(mapping, *, max_items=4): + items = _mapping_items(mapping) + pieces = [f"{_format_charge(key)}:{val}" for key, val in items[:max_items]] + if len(items) > max_items: + pieces.append("...") + return "{" + ", ".join(pieces) + "}" + + +def _format_charge(charge): + if isinstance(charge, tuple): + return "(" + ", ".join(str(x) for x in charge) + ")" + return str(charge) + + +def _format_sector(sector): + sector = _as_tuple(sector) + return "(" + ", ".join(_format_charge(x) for x in sector) + ")" + + +def _format_shape(shape): + return "x".join(str(int(dim)) for dim in shape) or "scalar" + + +def _add_charges(a, b): + if isinstance(a, tuple) or isinstance(b, tuple): + a_t = a if isinstance(a, tuple) else (a,) * len(b) + b_t = b if isinstance(b, tuple) else (b,) * len(a) + return tuple(x + y for x, y in zip(a_t, b_t)) + return a + b + + +def _sum_charges(charges): + charges = [charge for charge in charges if charge is not None] + if not charges: + return None + total = charges[0] + for charge in charges[1:]: + total = _add_charges(total, charge) + return total + + +def _mod_charge(charge, mod): + if charge is None: + return None + if isinstance(charge, tuple): + return tuple(int(x) % mod for x in charge) + return int(charge) % mod + + +def _resolve_total_charge(source, tensors): + if hasattr(source, "overall_charge"): + try: + return source.overall_charge() + except Exception: # pragma: no cover - summary should stay best-effort + pass + return _sum_charges(tensor.get("charge") for tensor in tensors) + + +def _resolve_q_total(symmetry, total_charge): + if total_charge is None: + return None + if str(symmetry).upper() == "Z2": + return _mod_charge(total_charge, 2) + return total_charge + + +def _charge_summary_text(summary): + pieces = [] + symmetry = summary.get("symmetry") + if symmetry is not None: + pieces.append(f"sym={symmetry}") + if summary.get("fermionic"): + pieces.append("fermionic") + total_charge = summary.get("charge_total", summary.get("total_charge")) + if total_charge is not None: + pieces.append(f"charge_total={_format_charge(total_charge)}") + q_total = summary.get("Q_total") + if q_total is None: + q_total = _resolve_q_total(symmetry, total_charge) + if str(symmetry).upper() == "Z2": + pieces.append(f"Q_total={_format_charge(q_total)} mod 2") + else: + pieces.append(f"Q_total={_format_charge(q_total)}") + return " | ".join(pieces) + + +def _flow_text(*directions): + return "->".join(str(direction) for direction in directions if direction is not None) + + +def _flow_math(*directions): + parts = [str(direction) for direction in directions if direction is not None] + if not parts: + return "" + return r"\to".join(rf"\mathrm{{{part}}}" for part in parts) + + +def _shape_size(shape): + return int(np.prod(tuple(shape), dtype=int)) if shape else 1 + + +def _summarize_symmray_index(index, axis): + chargemap = { + charge: int(size) + for charge, size in _mapping_items(dict(getattr(index, "chargemap", {}))) + } + dual = bool(getattr(index, "dual", False)) + return { + "axis": int(axis), + "dual": dual, + "direction": "in" if dual else "out", + "chargemap": chargemap, + "sectors": [ + {"charge": charge, "size": size} + for charge, size in _mapping_items(chargemap) + ], + "dim": int(sum(chargemap.values())), + "num_sectors": len(chargemap), + } + + +def symmray_block_summary(array): + """Return leg and block metadata for a Symmray block-sparse array. + + Parameters + ---------- + array : symmray array or quimb.Tensor + Block-sparse array exposing ``indices`` and ``get_sector_block_pairs``. + + Returns + ------- + dict + Summary with ``shape``, optional ``charge``, per-leg ``indices``, + present block-sector records under ``blocks``, and dense/stored size + counts useful for diagnostics. + """ + array = _require_symmray_array(array, name="array") + + indices = [ + _summarize_symmray_index(index, axis) + for axis, index in enumerate(getattr(array, "indices", ())) + ] + + blocks = [] + for sector, block in array.get_sector_block_pairs(): + shape = tuple(int(dim) for dim in getattr(block, "shape", ())) + size = int(getattr(block, "size", _shape_size(shape))) + blocks.append( + { + "sector": _as_tuple(sector), + "shape": shape, + "size": size, + "dtype": str(getattr(block, "dtype", "")), + } + ) + + shape = tuple(int(dim) for dim in getattr(array, "shape", ())) + dense_size = _shape_size(shape) + stored_size = int(sum(block["size"] for block in blocks)) + return { + "shape": shape, + "charge": getattr(array, "charge", None), + "indices": indices, + "blocks": blocks, + "num_blocks": len(blocks), + "dense_size": dense_size, + "stored_size": stored_size, + "density": stored_size / dense_size if dense_size else 0.0, + } + + +def _as_mps_tensor_network(value): + if hasattr(value, "tn"): + value = value.tn + elif hasattr(value, "p"): + value = value.p + if not hasattr(value, "sites") or not hasattr(value, "site_ind"): + raise TypeError("mps must be a SymMPS, MpsOptimizer, or quimb MPS object.") + return value + + +def _mps_sites(tn): + sites = tuple(getattr(tn, "sites", ())) + if not sites and hasattr(tn, "gen_sites_present"): + sites = tuple(tn.gen_sites_present()) + if not sites: + raise ValueError("mps does not expose any site labels.") + return sites + + +def _mps_site_tensor(tn, site): + try: + return tn[site] + except Exception as exc: # pragma: no cover - defensive for quimb variants + if hasattr(tn, "site_tag"): + try: + return tn[tn.site_tag(site)] + except Exception: + pass + raise ValueError(f"Could not resolve MPS tensor for site {site!r}.") from exc + + +def _shared_virtual_ind(tensor_a, tensor_b, physical_inds): + physical_inds = set(physical_inds) + for ind in tensor_a.inds: + if ind in tensor_b.inds and ind not in physical_inds: + return ind + return None + + +def _shared_virtual_inds(tensor_a, tensor_b, physical_inds): + physical_inds = set(physical_inds) + return tuple( + ind + for ind in tensor_a.inds + if ind in tensor_b.inds and ind not in physical_inds + ) + + +def _is_fermionic_array_data(data): + if data is None: + return False + if bool(getattr(data, "fermionic", False)): + return True + return "fermionic" in type(data).__name__.lower() + + +def _infer_fermionic(source, tensors): + fermionic = getattr(source, "fermionic", None) + if fermionic is not None: + return bool(fermionic) + return any(_is_fermionic_array_data(getattr(tensor, "data", None)) for tensor in tensors) + + +def _source_edges(source, bonds): + edges = getattr(source, "edges", None) + if edges is not None: + return _as_edges(edges) + return tuple(tuple(bond["between"]) for bond in bonds) + + +def _edge_lookup(edges): + order_by_pair = {} + edge_by_pair = {} + for position, edge in enumerate(edges): + edge = tuple(edge) + reverse_edge = tuple(reversed(edge)) + order_by_pair.setdefault(edge, int(position)) + order_by_pair.setdefault(reverse_edge, int(position)) + edge_by_pair.setdefault(edge, edge) + edge_by_pair.setdefault(reverse_edge, edge) + return order_by_pair, edge_by_pair + + +def _bond_endpoint_directions(bond): + if "left_site" in bond: + return { + bond["left_site"]: bond.get("left_direction"), + bond["right_site"]: bond.get("right_direction"), + } + return { + bond["site_a"]: bond.get("site_a_direction"), + bond["site_b"]: bond.get("site_b_direction"), + } + + +def _fermionic_edge_record(bond, order_by_pair, edge_by_pair): + between = tuple(bond["between"]) + edge = edge_by_pair.get(between, between) + directions = _bond_endpoint_directions(bond) + return { + "position": int(bond["position"]), + "edge_order": order_by_pair.get(between), + "edge": edge, + "between": between, + "ind": bond["ind"], + "lattice_direction": bond.get("direction"), + "index_directions": tuple( + {"site": site, "direction": directions.get(site)} + for site in edge + ), + } + + +def _fermionic_ordering_summary(source, *, network_kind, sites, bonds, fermionic): + edges = _source_edges(source, bonds) + order_by_pair, edge_by_pair = _edge_lookup(edges) + return { + "enabled": bool(fermionic), + "network_kind": network_kind, + "methods_reference": dict(_FERMIONIC_TN_METHODS_REFERENCE), + "site_order": tuple(sites), + "edge_order": edges, + "edges": tuple( + _fermionic_edge_record(bond, order_by_pair, edge_by_pair) + for bond in bonds + ), + } + + +def _resolve_mps_position(tensors, value, *, name="position"): + if value is None: + return None + if value == "middle": + return tensors[len(tensors) // 2]["position"] if tensors else None + for tensor in tensors: + if value == tensor["position"] or value == tensor["site"]: + return tensor["position"] + raise ValueError(f"{name}={value!r} does not identify a shown MPS site.") + + +def symmray_mps_summary(mps): + """Return site, bond, and block metadata for a Symmray-backed MPS. + + Parameters + ---------- + mps : SymMPS, MpsOptimizer, or quimb MatrixProductState + Object whose site tensors store Symmray block-sparse arrays. + + Returns + ------- + dict + Summary with per-site tensor block counts, physical-sector maps, + nearest-neighbor bond-sector maps, and aggregate storage diagnostics. + """ + source = mps + tn = _as_mps_tensor_network(mps) + sites = _mps_sites(tn) + + tensors = [] + site_tensors = [] + index_maps = [] + physical_inds = [] + + for position, site in enumerate(sites): + tensor = _mps_site_tensor(tn, site) + array_summary = symmray_block_summary(tensor) + site_ind = tn.site_ind(site) + + index_by_ind = {} + indices = [] + for ind, index_summary in zip(tensor.inds, array_summary["indices"]): + entry = dict(index_summary) + entry["ind"] = ind + index_by_ind[ind] = entry + indices.append(entry) + + physical = index_by_ind.get(site_ind) + if physical is None: + raise ValueError( + f"MPS tensor for site {site!r} does not expose physical index {site_ind!r}." + ) + + tensors.append( + { + "position": int(position), + "site": site, + "site_tag": tn.site_tag(site) if hasattr(tn, "site_tag") else None, + "site_ind": site_ind, + "inds": tuple(tensor.inds), + "shape": array_summary["shape"], + "charge": array_summary["charge"], + "indices": indices, + "physical": physical, + "left_bond": None, + "right_bond": None, + "blocks": array_summary["blocks"], + "num_blocks": array_summary["num_blocks"], + "dense_size": array_summary["dense_size"], + "stored_size": array_summary["stored_size"], + "density": array_summary["density"], + } + ) + site_tensors.append(tensor) + index_maps.append(index_by_ind) + physical_inds.append(site_ind) + + bonds = [] + for position, (site_l, site_r) in enumerate(zip(sites[:-1], sites[1:])): + ind = _shared_virtual_ind( + site_tensors[position], + site_tensors[position + 1], + (physical_inds[position], physical_inds[position + 1]), + ) + if ind is None: + continue + left_index = index_maps[position].get(ind) + right_index = index_maps[position + 1].get(ind) + index_summary = left_index or right_index + bond = { + "position": int(position), + "left_position": int(position), + "right_position": int(position + 1), + "left_site": site_l, + "right_site": site_r, + "between": (site_l, site_r), + "ind": ind, + "left_direction": left_index["direction"] if left_index is not None else None, + "right_direction": right_index["direction"] if right_index is not None else None, + "dim": index_summary["dim"], + "num_sectors": index_summary["num_sectors"], + "chargemap": index_summary["chargemap"], + "sectors": index_summary["sectors"], + } + bonds.append(bond) + tensors[position]["right_bond"] = bond + tensors[position + 1]["left_bond"] = bond + + total_dense_size = int(sum(tensor["dense_size"] for tensor in tensors)) + total_stored_size = int(sum(tensor["stored_size"] for tensor in tensors)) + total_charge = _resolve_total_charge(source, tensors) + symmetry = getattr(source, "symmetry", None) + fermionic = _infer_fermionic(source, site_tensors) + q_total = _resolve_q_total(symmetry, total_charge) + return { + "num_sites": len(sites), + "sites": sites, + "tensors": tensors, + "bonds": bonds, + "symmetry": symmetry, + "fermionic": fermionic, + "fermionic_ordering": _fermionic_ordering_summary( + source, + network_kind="mps", + sites=sites, + bonds=bonds, + fermionic=fermionic, + ), + "total_charge": total_charge, + "charge_total": total_charge, + "Q_total": q_total, + "total_parity": _mod_charge(total_charge, 2), + "max_bond_dim": max((bond["dim"] for bond in bonds), default=1), + "max_bond_sectors": max((bond["num_sectors"] for bond in bonds), default=0), + "total_dense_size": total_dense_size, + "total_stored_size": total_stored_size, + "density": total_stored_size / total_dense_size if total_dense_size else 0.0, + } + + +def draw_symmray_blocks( + array, + *, + ax=None, + title=None, + max_blocks=12, + show_leg_chargemaps=True, + figsize=None, + return_summary=False, +): + """Draw a lightweight sector schematic for a Symmray array. + + The diagram uses :mod:`quimb.schematic` to show array legs, their charge + maps, and the present block sectors with block shapes. + """ + summary = symmray_block_summary(array) + blocks = summary["blocks"] + max_blocks = int(max_blocks) + if max_blocks < 1: + raise ValueError("max_blocks must be >= 1.") + shown_blocks = blocks[:max_blocks] + + try: + from quimb import schematic # pylint: disable=import-outside-toplevel + except ImportError as exc: # pragma: no cover - quimb is a declared dep + raise ImportError("draw_symmray_blocks requires quimb.schematic.") from exc + + try: + from matplotlib.patches import Rectangle # pylint: disable=import-outside-toplevel + except ImportError as exc: # pragma: no cover - plotting optionality + raise ImportError("draw_symmray_blocks requires matplotlib.") from exc + + rank = max(1, len(summary["indices"])) + block_count = max(1, len(shown_blocks)) + if figsize is None: + figsize = ( + max(5.5, 1.15 * max(rank, block_count) + 2.2), + 4.2, + ) + + presets = { + "array": { + "color": schematic.get_color("blue"), + "alpha": 0.28, + "linewidth": 1.6, + }, + "leg": { + "color": (0.32, 0.35, 0.38, 1.0), + "linewidth": 1.5, + }, + "block": { + "alpha": 0.62, + "linewidth": 1.2, + }, + } + drawing = schematic.Drawing(presets=presets, ax=ax, figsize=figsize) + + drawing.ax.add_patch( + Rectangle( + (-1.05, -0.38), + 2.10, + 0.76, + facecolor=schematic.get_color("blue"), + edgecolor=(0.12, 0.25, 0.34, 1.0), + alpha=0.70, + linewidth=1.3, + zorder=3, + ) + ) + center_label = ( + f"shape {_format_shape(summary['shape'])}" + f"\n{summary['num_blocks']} blocks" + f"\n{summary['stored_size']}/{summary['dense_size']} stored" + ) + if summary["charge"] is not None: + center_label += f"\ncharge {summary['charge']}" + drawing.text( + (0.0, 0.0), + center_label, + fontsize=10, + ha="center", + va="center", + color="white", + zorder=4, + ) + if title is not None: + drawing.ax.set_title(str(title)) + + if summary["indices"]: + if len(summary["indices"]) == 1: + xs = [0.0] + else: + xs = np.linspace(-1.8, 1.8, len(summary["indices"])) + for x_pos, index in zip(xs, summary["indices"]): + drawing.line((0.0, 0.35), (float(x_pos), 0.95), preset="leg") + label = f"axis {index['axis']} ({index['direction']})\ndim {index['dim']}" + if show_leg_chargemaps: + label += "\n" + _format_compact_mapping(index["chargemap"]) + drawing.text((float(x_pos), 1.08), label, fontsize=8, ha="center", va="bottom") + + if shown_blocks: + if len(shown_blocks) == 1: + block_xs = [0.0] + else: + block_xs = np.linspace(-1.9, 1.9, len(shown_blocks)) + for x_pos, block in zip(block_xs, shown_blocks): + color = schematic.hash_to_color(str(block["sector"])) + drawing.ax.add_patch( + Rectangle( + (float(x_pos) - 0.42, -1.28), + 0.84, + 0.56, + facecolor=color, + edgecolor=(0.16, 0.18, 0.21, 0.75), + alpha=0.78, + linewidth=0.6, + zorder=3, + ) + ) + drawing.text( + (float(x_pos), -0.91), + _format_sector(block["sector"]), + fontsize=8, + ha="center", + ) + drawing.text( + (float(x_pos), -1.43), + _format_shape(block["shape"]), + fontsize=8, + ha="center", + ) + drawing.text((0.0, -0.55), "present blocks", fontsize=9, ha="center") + if len(blocks) > len(shown_blocks): + drawing.text( + (2.45, -1.0), + f"+{len(blocks) - len(shown_blocks)} more", + fontsize=8, + ha="left", + ) + else: + drawing.text((0.0, -1.0), "no stored blocks", fontsize=9, ha="center") + + drawing.ax.set_xlim(-2.55, 2.55) + drawing.ax.set_ylim(-1.62, 1.58) + drawing.ax.axis("off") + if return_summary: + return drawing, summary + return drawing + + +def draw_symmray_mps( + mps, + *, + ax=None, + title=None, + max_sites=None, + center="middle", + highlight_pair=True, + show_regions=True, + show_arrows=True, + show_leg_chargemaps=True, + show_bond_labels=False, + show_phys_labels=False, + show_tensor_labels=True, + show_diagnostics=False, + show_blocks=False, + show_block_labels=False, + max_blocks_per_site=4, + node_shape="circle", + node_radius=0.24, + figsize=None, + return_summary=False, +): + """Draw a block-aware MPS schematic for a Symmray-backed MPS. + + The diagram uses :mod:`quimb.schematic` in the same style as quimb's manual + tensor-network schematics: tensor nodes, virtual bonds, physical legs, + optional canonical-flow arrows, and optional left/center/right region + highlighting. Symmray-specific dimensions remain available through + ``return_summary=True`` and can be overlaid with the label options. + """ + summary = symmray_mps_summary(mps) + max_blocks_per_site = int(max_blocks_per_site) + if max_blocks_per_site < 1: + raise ValueError("max_blocks_per_site must be >= 1.") + if max_sites is None: + shown_tensors = list(summary["tensors"]) + else: + max_sites = int(max_sites) + if max_sites < 1: + raise ValueError("max_sites must be >= 1.") + shown_tensors = list(summary["tensors"][:max_sites]) + + shown_positions = {tensor["position"] for tensor in shown_tensors} + center_position = _resolve_mps_position(shown_tensors, center, name="center") + pair_right_position = None + if highlight_pair and center_position is not None: + candidate = center_position + 1 + if candidate in shown_positions: + pair_right_position = candidate + shown_bonds = [ + bond + for bond in summary["bonds"] + if bond["left_position"] in shown_positions + and bond["right_position"] in shown_positions + ] + + try: + from quimb import schematic # pylint: disable=import-outside-toplevel + except ImportError as exc: # pragma: no cover - quimb is a declared dep + raise ImportError("draw_symmray_mps requires quimb.schematic.") from exc + + try: + from matplotlib.patches import Rectangle # pylint: disable=import-outside-toplevel + except ImportError as exc: # pragma: no cover - plotting optionality + raise ImportError("draw_symmray_mps requires matplotlib.") from exc + + n_show = len(shown_tensors) + detailed_labels = bool(show_leg_chargemaps and (show_bond_labels or show_phys_labels)) + if detailed_labels: + spacing = 1.72 if n_show <= 12 else 1.32 + else: + spacing = 1.0 if n_show > 12 else 1.12 + if figsize is None: + if show_bond_labels or show_phys_labels or show_blocks: + height = 5.25 if detailed_labels else 3.9 + if show_block_labels: + height += 0.35 + if show_diagnostics: + height += 0.15 + figsize = (max(8.0, spacing * max(1, n_show - 1) + 4.0), height) + else: + figsize = (max(6.5, 0.95 * n_show + 1.5), 2.9) + + presets = { + "bond": { + "color": (0.12, 0.14, 0.16, 1.0), + "linewidth": 3.0, + }, + "phys": { + "color": (0.12, 0.14, 0.16, 1.0), + "linewidth": 1.55, + }, + "center": { + "color": schematic.get_color("orange"), + "hatch": "/////", + "linewidth": 1.3, + "edgecolor": (0.55, 0.38, 0.00, 1.0), + }, + "left": { + "color": schematic.get_color("bluedark"), + "linewidth": 1.3, + "edgecolor": (0.02, 0.22, 0.34, 1.0), + }, + "right": { + "color": schematic.get_color("blue"), + "linewidth": 1.3, + "edgecolor": (0.05, 0.34, 0.50, 1.0), + }, + "pair": { + "facecolor": (0.20, 0.80, 0.50, 0.34), + "edgecolor": (0.18, 0.50, 0.30, 0.72), + "linestyle": ":", + "linewidth": 1.5, + }, + "region": { + "facecolor": (0.83, 0.83, 0.83, 0.42), + "edgecolor": (0.42, 0.42, 0.42, 0.72), + "linestyle": ":", + "linewidth": 1.4, + }, + } + drawing = schematic.Drawing(presets=presets, ax=ax, figsize=figsize) + + x_by_position = { + tensor["position"]: i * spacing + for i, tensor in enumerate(shown_tensors) + } + y0 = 0.0 + phys_y = -1.05 if (show_phys_labels and show_blocks) else -0.68 + block_y = -0.50 if (show_phys_labels and show_blocks) else -0.46 + max_dim = max((bond["dim"] for bond in shown_bonds), default=1) + + left_region = [] + right_region = [] + if center_position is not None: + right_region_start = ( + pair_right_position + 1 + if pair_right_position is not None + else center_position + 1 + ) + left_region = [ + (x_by_position[tensor["position"]], y0) + for tensor in shown_tensors + if tensor["position"] < center_position + ] + right_region = [ + (x_by_position[tensor["position"]], y0) + for tensor in shown_tensors + if tensor["position"] >= right_region_start + ] + if show_regions and left_region: + drawing.patch_around(left_region, radius=0.46, preset="region", zorder=0) + drawing.text( + (left_region[-1][0] - 0.25 * spacing, 1.36 if detailed_labels else 0.90), + "LEFT", + fontsize=12, + color=(0.18, 0.19, 0.21, 1.0), + ha="center", + ) + if show_regions and right_region: + drawing.patch_around(right_region, radius=0.46, preset="region", zorder=0) + drawing.text( + (right_region[0][0] + 0.35 * spacing, 1.36 if detailed_labels else 0.90), + "RIGHT", + fontsize=12, + color=(0.18, 0.19, 0.21, 1.0), + ha="center", + ) + if show_regions and pair_right_position is not None: + drawing.patch_around_circles( + (x_by_position[center_position], y0), + node_radius + 0.04, + (x_by_position[pair_right_position], y0), + node_radius + 0.04, + padding=0.22, + preset="pair", + zorder=0, + ) + + for bond in shown_bonds: + x0 = x_by_position[bond["left_position"]] + x1 = x_by_position[bond["right_position"]] + width = 2.4 + 1.2 * (bond["dim"] / max_dim) + drawing.line((x0, y0), (x1, y0), preset="bond", linewidth=width, zorder=1) + + if show_arrows: + if center_position is None: + drawing.arrowhead((x0, y0), (x1, y0), preset="bond", center=0.58, width=0.10) + elif bond["right_position"] <= center_position: + drawing.arrowhead((x0, y0), (x1, y0), preset="bond", center=0.58, width=0.10) + elif ( + pair_right_position is not None + and bond["left_position"] >= pair_right_position + ) or ( + pair_right_position is None + and bond["left_position"] >= center_position + ): + drawing.arrowhead((x1, y0), (x0, y0), preset="bond", center=0.58, width=0.10) + + mid = 0.5 * (x0 + x1) + if show_bond_labels: + flow = _flow_math(bond.get("left_direction"), bond.get("right_direction")) + if flow: + label = rf"$e_{{{bond['position']}}}: {flow}, \chi={bond['dim']}$" + else: + label = rf"$e_{{{bond['position']}}}: \chi={bond['dim']}$" + if show_leg_chargemaps: + label += "\n" + rf"$q_e:$ {_format_compact_mapping(bond['chargemap'], max_items=4)}" + label_y = 0.28 + (0.18 if detailed_labels and bond["position"] % 2 else 0.0) + drawing.text( + (mid, label_y), + label, + fontsize=6.6, + ha="center", + va="bottom", + color=(0.18, 0.20, 0.23, 1.0), + zorder=5, + ) + + show_block_labels = bool(show_block_labels) + for tensor in shown_tensors: + x_pos = x_by_position[tensor["position"]] + position = tensor["position"] + if position == center_position: + preset = "center" + elif pair_right_position is not None and position == pair_right_position: + preset = "right" + elif center_position is not None and position < center_position: + preset = "left" + else: + preset = "right" + + if node_shape == "circle": + drawing.circle((x_pos, y0), radius=node_radius, preset=preset, zorder=3) + elif node_shape == "cube": + drawing.cube((x_pos, y0, 0.0), preset=preset, zorder=3) + else: + raise ValueError("node_shape must be 'circle' or 'cube'.") + + drawing.line( + (x_pos, y0), + (x_pos, phys_y), + preset="phys", + zorder=1, + ) + if show_arrows and position != center_position: + drawing.arrowhead( + (x_pos, phys_y), + (x_pos, y0), + preset="phys", + center=0.55, + width=0.09, + ) + drawing.dot( + (x_pos, phys_y), + color=(0.12, 0.14, 0.16, 1.0), + radius=0.028, + zorder=2, + ) + + if show_tensor_labels: + label_lines = [rf"$T_{{{tensor['site']}}}$"] + if tensor["charge"] is not None: + label_lines.append(rf"$q={_format_charge(tensor['charge'])}$") + if not show_blocks: + label_lines.append(rf"$B={tensor['num_blocks']}$") + label = "\n".join(label_lines) + label_color = ( + schematic.get_color("orange") + if position == center_position + else (0.06, 0.20, 0.30, 1.0) + ) + tensor_label_y = 0.76 if detailed_labels else 0.54 + drawing.text( + (x_pos, tensor_label_y), + label, + fontsize=8.5, + ha="center", + va="bottom", + color=label_color, + zorder=5, + ) + + physical = tensor["physical"] + if show_phys_labels: + phys_label = ( + rf"$p_{{{tensor['site']}}}: \mathrm{{{physical['direction']}}}, d={physical['dim']}$" + ) + if show_leg_chargemaps: + phys_label += "\n" + rf"$q_p:$ {_format_compact_mapping(physical['chargemap'], max_items=4)}" + drawing.text( + (x_pos, phys_y - 0.12), + phys_label, + fontsize=6.4, + ha="center", + va="top", + color=(0.18, 0.20, 0.23, 1.0), + ) + + if show_blocks and tensor["blocks"]: + blocks = tensor["blocks"][:max_blocks_per_site] + block_width = 0.18 + block_height = 0.14 + gap = 0.035 + total_width = len(blocks) * block_width + (len(blocks) - 1) * gap + start = x_pos - 0.5 * total_width + for block_pos, block in enumerate(blocks): + bx = start + block_pos * (block_width + gap) + color = schematic.hash_to_color(str(block["sector"])) + drawing.ax.add_patch( + Rectangle( + (bx, block_y), + block_width, + block_height, + facecolor=color, + edgecolor=(0.16, 0.18, 0.21, 0.75), + alpha=0.86, + linewidth=0.65, + zorder=4, + ) + ) + if show_block_labels: + drawing.ax.text( + bx + 0.5 * block_width, + block_y - 0.045, + _format_sector(block["sector"]), + fontsize=5.2, + ha="center", + va="top", + rotation=45, + color=(0.15, 0.17, 0.20, 1.0), + ) + drawing.ax.text( + x_pos, + block_y + block_height + 0.025, + rf"$B={tensor['num_blocks']}$", + fontsize=6.5, + ha="center", + va="bottom", + color=(0.15, 0.17, 0.20, 1.0), + bbox={ + "boxstyle": "round,pad=0.06", + "facecolor": (1.0, 1.0, 1.0, 0.78), + "edgecolor": (1.0, 1.0, 1.0, 0.0), + "linewidth": 0.0, + }, + zorder=5, + ) + + last_x = x_by_position[shown_tensors[-1]["position"]] if shown_tensors else 0.0 + right_pad = 0.75 + if show_diagnostics: + charge_line = _charge_summary_text(summary) + diagnostic_lines = [f"sites {summary['num_sites']}"] + if charge_line: + diagnostic_lines.append(charge_line) + diagnostic_lines += [ + f"max bond {summary['max_bond_dim']}", + f"bond sectors {summary['max_bond_sectors']}", + f"stored {summary['total_stored_size']}/{summary['total_dense_size']}", + f"density {summary['density']:.3f}", + ] + diagnostic = "\n".join(diagnostic_lines) + if show_blocks: + diagnostic += "\ncolored tiles: stored blocks" + if show_arrows: + diagnostic += "\narrows/labels: charge in/out flow" + if len(shown_tensors) < summary["num_sites"]: + diagnostic += f"\n+{summary['num_sites'] - len(shown_tensors)} sites hidden" + summary_x = last_x + 0.82 + drawing.ax.text( + summary_x, + 0.52, + diagnostic, + fontsize=8, + ha="left", + va="top", + color=(0.15, 0.17, 0.20, 1.0), + bbox={ + "boxstyle": "round,pad=0.22", + "facecolor": (1.0, 1.0, 1.0, 0.92), + "edgecolor": (0.68, 0.70, 0.74, 1.0), + "linewidth": 0.8, + }, + ) + right_pad = 1.80 + elif len(shown_tensors) < summary["num_sites"]: + drawing.text( + (last_x + 0.55, y0), + f"+{summary['num_sites'] - len(shown_tensors)}", + fontsize=9, + ha="left", + va="center", + color=(0.18, 0.20, 0.23, 1.0), + ) + right_pad = 1.10 + + if title is None: + title = "Symmray MPS block structure" + drawing.ax.set_title(title) + charge_text = _charge_summary_text(summary) + if charge_text and not show_diagnostics: + drawing.text( + (-0.55, 1.12), + charge_text, + fontsize=8, + ha="left", + va="center", + color=(0.18, 0.20, 0.23, 1.0), + ) + drawing.ax.set_xlim(-0.65, last_x + right_pad) + y_min = -1.55 if (show_phys_labels and show_blocks) else (-1.18 if (show_phys_labels or show_block_labels) else -0.96) + y_max = 1.78 if detailed_labels else (1.35 if (show_regions or show_leg_chargemaps) else 0.92) + drawing.ax.set_ylim(y_min, y_max) + drawing.ax.axis("off") + if return_summary: + return drawing, summary + return drawing + + +def _as_peps_tensor_network(value): + if hasattr(value, "tn"): + value = value.tn + elif hasattr(value, "state"): + value = value.state + if not hasattr(value, "sites") or not hasattr(value, "site_ind"): + raise TypeError("peps must be a SymPEPS, PepsOptimizer, or quimb PEPS object.") + if not hasattr(value, "Lx") or not hasattr(value, "Ly"): + raise TypeError("peps must expose PEPS lattice dimensions Lx and Ly.") + return value + + +def _peps_sites(tn): + sites = tuple(getattr(tn, "sites", ())) + if not sites and hasattr(tn, "gen_site_coos"): + sites = tuple(tn.gen_site_coos()) + if not sites: + sites = tuple((i, j) for i in range(int(tn.Lx)) for j in range(int(tn.Ly))) + if not all(isinstance(site, tuple) and len(site) == 2 for site in sites): + raise ValueError("peps sites must be two-dimensional coordinate tuples.") + return tuple(sorted(tuple(int(x) for x in site) for site in sites)) + + +def _peps_site_tensor(tn, site): + try: + return tn[site] + except Exception as exc: # pragma: no cover - defensive for quimb variants + if hasattr(tn, "site_tag"): + try: + return tn[tn.site_tag(site)] + except Exception: + pass + raise ValueError(f"Could not resolve PEPS tensor for site {site!r}.") from exc + + +def _peps_site_tag(tn, site): + if hasattr(tn, "site_tag"): + return tn.site_tag(site) + return None + + +def _peps_site_ind(tn, site): + if hasattr(tn, "site_ind"): + return tn.site_ind(site) + return _format_site_ind(site, getattr(tn, "site_ind_id", "k{},{}")) + + +def _peps_relative_direction(site_a, site_b): + dx = int(site_b[0]) - int(site_a[0]) + dy = int(site_b[1]) - int(site_a[1]) + if dx == 0 and dy == 1: + return "right" + if dx == 0 and dy == -1: + return "left" + if dx == 1 and dy == 0: + return "down" + if dx == -1 and dy == 0: + return "up" + return "other" + + +def _opposite_direction(direction): + return { + "right": "left", + "left": "right", + "down": "up", + "up": "down", + }.get(direction, "other") + + +def symmray_peps_summary(peps): + """Return site, bond, and block metadata for a Symmray-backed PEPS. + + Parameters + ---------- + peps : SymPEPS, PepsOptimizer, or quimb PEPS + Object whose site tensors store Symmray block-sparse arrays. + + Returns + ------- + dict + Summary with per-site tensor block counts, physical-sector maps, + virtual-bond sector maps, lattice dimensions, and aggregate storage + diagnostics. + """ + source = peps + tn = _as_peps_tensor_network(peps) + sites = _peps_sites(tn) + site_set = set(sites) + + tensors = [] + tensors_by_site = {} + site_tensors = {} + index_maps = {} + physical_inds = {} + + for position, site in enumerate(sites): + tensor = _peps_site_tensor(tn, site) + array_summary = symmray_block_summary(tensor) + site_ind = _peps_site_ind(tn, site) + + index_by_ind = {} + indices = [] + for ind, index_summary in zip(tensor.inds, array_summary["indices"]): + entry = dict(index_summary) + entry["ind"] = ind + index_by_ind[ind] = entry + indices.append(entry) + + physical = index_by_ind.get(site_ind) + if physical is None: + raise ValueError( + f"PEPS tensor for site {site!r} does not expose physical index {site_ind!r}." + ) + + entry = { + "position": int(position), + "site": site, + "site_tag": _peps_site_tag(tn, site), + "site_ind": site_ind, + "inds": tuple(tensor.inds), + "shape": array_summary["shape"], + "charge": array_summary["charge"], + "indices": indices, + "physical": physical, + "bonds": {}, + "blocks": array_summary["blocks"], + "num_blocks": array_summary["num_blocks"], + "dense_size": array_summary["dense_size"], + "stored_size": array_summary["stored_size"], + "density": array_summary["density"], + } + tensors.append(entry) + tensors_by_site[site] = entry + site_tensors[site] = tensor + index_maps[site] = index_by_ind + physical_inds[site] = site_ind + + bonds = [] + for left_pos, site_a in enumerate(sites): + for site_b in sites[left_pos + 1 :]: + shared_inds = _shared_virtual_inds( + site_tensors[site_a], + site_tensors[site_b], + (physical_inds[site_a], physical_inds[site_b]), + ) + for ind in shared_inds: + index_a = index_maps[site_a].get(ind) + index_b = index_maps[site_b].get(ind) + index_summary = index_a or index_b + direction = _peps_relative_direction(site_a, site_b) + bond = { + "position": len(bonds), + "site_a": site_a, + "site_b": site_b, + "between": (site_a, site_b), + "ind": ind, + "direction": direction, + "site_a_direction": index_a["direction"] if index_a is not None else None, + "site_b_direction": index_b["direction"] if index_b is not None else None, + "dim": index_summary["dim"], + "num_sectors": index_summary["num_sectors"], + "chargemap": index_summary["chargemap"], + "sectors": index_summary["sectors"], + } + bonds.append(bond) + tensors_by_site[site_a]["bonds"][direction] = bond + tensors_by_site[site_b]["bonds"][_opposite_direction(direction)] = bond + + total_dense_size = int(sum(tensor["dense_size"] for tensor in tensors)) + total_stored_size = int(sum(tensor["stored_size"] for tensor in tensors)) + Lx = int(getattr(tn, "Lx", max(site[0] for site in site_set) + 1)) + Ly = int(getattr(tn, "Ly", max(site[1] for site in site_set) + 1)) + total_charge = _resolve_total_charge(source, tensors) + symmetry = getattr(source, "symmetry", None) + fermionic = _infer_fermionic(source, site_tensors.values()) + q_total = _resolve_q_total(symmetry, total_charge) + return { + "Lx": Lx, + "Ly": Ly, + "num_sites": len(sites), + "sites": sites, + "tensors": tensors, + "bonds": bonds, + "symmetry": symmetry, + "fermionic": fermionic, + "fermionic_ordering": _fermionic_ordering_summary( + source, + network_kind="peps", + sites=sites, + bonds=bonds, + fermionic=fermionic, + ), + "total_charge": total_charge, + "charge_total": total_charge, + "Q_total": q_total, + "total_parity": _mod_charge(total_charge, 2), + "max_bond_dim": max((bond["dim"] for bond in bonds), default=1), + "max_bond_sectors": max((bond["num_sectors"] for bond in bonds), default=0), + "total_dense_size": total_dense_size, + "total_stored_size": total_stored_size, + "density": total_stored_size / total_dense_size if total_dense_size else 0.0, + } + + +def _resolve_peps_center(tensors, value, *, Lx, Ly, name="center"): + if value is None: + return None + sites = {tensor["site"] for tensor in tensors} + if value == "middle": + target_x = (int(Lx) - 1) / 2 + target_y = (int(Ly) - 1) / 2 + return min( + sites, + key=lambda site: ( + abs(site[0] - target_x) + abs(site[1] - target_y), + site[0], + site[1], + ), + ) + if value in sites: + return tuple(value) + raise ValueError(f"{name}={value!r} does not identify a shown PEPS site.") + + +def _peps_site_distance(site, center): + return abs(int(site[0]) - int(center[0])) + abs(int(site[1]) - int(center[1])) + + +def draw_symmray_peps( + peps, + *, + ax=None, + title=None, + max_sites=None, + center="middle", + show_region=True, + show_arrows=True, + show_leg_chargemaps=True, + show_bond_labels=False, + show_phys_labels=False, + show_tensor_labels=True, + show_diagnostics=False, + show_blocks=False, + show_block_labels=False, + max_blocks_per_site=4, + node_shape="circle", + node_radius=0.22, + figsize=None, + return_summary=False, +): + """Draw a block-aware PEPS schematic for a Symmray-backed PEPS. + + The schematic follows the compact :mod:`quimb.schematic` style with a PEPS + lattice, virtual-bond arrows, physical legs, and optional Symmray block and + dimension labels. + """ + summary = symmray_peps_summary(peps) + max_blocks_per_site = int(max_blocks_per_site) + if max_blocks_per_site < 1: + raise ValueError("max_blocks_per_site must be >= 1.") + if max_sites is None: + shown_tensors = list(summary["tensors"]) + else: + max_sites = int(max_sites) + if max_sites < 1: + raise ValueError("max_sites must be >= 1.") + shown_tensors = list(summary["tensors"][:max_sites]) + + shown_sites = {tensor["site"] for tensor in shown_tensors} + center_site = _resolve_peps_center( + shown_tensors, + center, + Lx=summary["Lx"], + Ly=summary["Ly"], + name="center", + ) + shown_bonds = [ + bond + for bond in summary["bonds"] + if bond["site_a"] in shown_sites and bond["site_b"] in shown_sites + ] + + try: + from quimb import schematic # pylint: disable=import-outside-toplevel + except ImportError as exc: # pragma: no cover - quimb is a declared dep + raise ImportError("draw_symmray_peps requires quimb.schematic.") from exc + + try: + from matplotlib.patches import Rectangle # pylint: disable=import-outside-toplevel + except ImportError as exc: # pragma: no cover - plotting optionality + raise ImportError("draw_symmray_peps requires matplotlib.") from exc + + spacing = 1.16 + if figsize is None: + width = max(4.8, 1.15 * max(1, summary["Ly"]) + 1.8) + height = max(4.4, 1.20 * max(1, summary["Lx"]) + 1.9) + if show_bond_labels or show_phys_labels or show_blocks: + width += 0.85 if show_leg_chargemaps else 0.55 + height += 0.55 if show_leg_chargemaps else 0.35 + if show_diagnostics: + width += 1.35 + figsize = (width, height) + + presets = { + "bond": { + "color": (0.12, 0.14, 0.16, 1.0), + "linewidth": 3.0, + }, + "phys": { + "color": (0.12, 0.14, 0.16, 1.0), + "linewidth": 1.45, + }, + "center": { + "color": schematic.get_color("orange"), + "hatch": "/////", + "linewidth": 1.2, + "edgecolor": (0.55, 0.38, 0.00, 1.0), + }, + "site_a": { + "color": schematic.get_color("bluedark"), + "linewidth": 1.2, + "edgecolor": (0.02, 0.22, 0.34, 1.0), + }, + "site_b": { + "color": schematic.get_color("blue"), + "linewidth": 1.2, + "edgecolor": (0.05, 0.34, 0.50, 1.0), + }, + "region": { + "facecolor": (0.83, 0.83, 0.83, 0.38), + "edgecolor": (0.42, 0.42, 0.42, 0.72), + "linestyle": ":", + "linewidth": 1.35, + }, + } + drawing = schematic.Drawing(presets=presets, ax=ax, figsize=figsize) + + xy_by_site = { + tensor["site"]: (tensor["site"][1] * spacing, -tensor["site"][0] * spacing) + for tensor in shown_tensors + } + max_dim = max((bond["dim"] for bond in shown_bonds), default=1) + + if show_region and shown_tensors: + drawing.patch_around(list(xy_by_site.values()), radius=0.42, preset="region", zorder=0) + + for bond in shown_bonds: + xy_a = xy_by_site[bond["site_a"]] + xy_b = xy_by_site[bond["site_b"]] + width = 2.25 + 1.05 * (bond["dim"] / max_dim) + drawing.line(xy_a, xy_b, preset="bond", linewidth=width, zorder=1) + + if show_arrows: + start, stop = xy_a, xy_b + if center_site is not None: + dist_a = _peps_site_distance(bond["site_a"], center_site) + dist_b = _peps_site_distance(bond["site_b"], center_site) + if dist_a < dist_b: + start, stop = xy_b, xy_a + drawing.arrowhead(start, stop, preset="bond", center=0.58, width=0.10) + + if show_bond_labels: + mid = (0.5 * (xy_a[0] + xy_b[0]), 0.5 * (xy_a[1] + xy_b[1])) + offset = (0.0, 0.16) if bond["direction"] in {"left", "right"} else (0.17, 0.0) + flow = _flow_math(bond.get("site_a_direction"), bond.get("site_b_direction")) + if flow: + label = rf"$e_{{{bond['position']}}}: {flow}, \chi={bond['dim']}$" + else: + label = rf"$e_{{{bond['position']}}}: \chi={bond['dim']}$" + if show_leg_chargemaps: + label += "\n" + rf"$q_e:$ {_format_compact_mapping(bond['chargemap'], max_items=4)}" + drawing.text( + (mid[0] + offset[0], mid[1] + offset[1]), + label, + fontsize=6.2, + ha="center" if offset[0] == 0.0 else "left", + va="bottom" if offset[1] > 0.0 else "center", + color=(0.18, 0.20, 0.23, 1.0), + zorder=5, + ) + + show_block_labels = bool(show_block_labels) + for tensor in shown_tensors: + site = tensor["site"] + x_pos, y_pos = xy_by_site[site] + preset = "center" if site == center_site else ("site_a" if sum(site) % 2 == 0 else "site_b") + + if node_shape == "circle": + drawing.circle((x_pos, y_pos), radius=node_radius, preset=preset, zorder=3) + elif node_shape == "cube": + drawing.cube((x_pos, y_pos, 0.0), preset=preset, zorder=3) + else: + raise ValueError("node_shape must be 'circle' or 'cube'.") + + phys_xy = (x_pos - 0.30, y_pos - 0.42) + drawing.line((x_pos, y_pos), phys_xy, preset="phys", zorder=1) + if show_arrows: + drawing.arrowhead(phys_xy, (x_pos, y_pos), preset="phys", center=0.56, width=0.08) + drawing.dot( + phys_xy, + color=(0.12, 0.14, 0.16, 1.0), + radius=0.025, + zorder=2, + ) + + if show_tensor_labels: + label_lines = [rf"$T_{{({site[0]},{site[1]})}}$"] + if tensor["charge"] is not None: + label_lines.append(rf"$q={_format_charge(tensor['charge'])}$") + if not show_blocks: + label_lines.append(rf"$B={tensor['num_blocks']}$") + label = "\n".join(label_lines) + label_color = ( + schematic.get_color("orange") + if site == center_site + else (0.06, 0.20, 0.30, 1.0) + ) + drawing.text( + (x_pos, y_pos + 0.34), + label, + fontsize=7.5, + ha="center", + va="bottom", + color=label_color, + zorder=5, + ) + + physical = tensor["physical"] + if show_phys_labels: + phys_label = ( + rf"$p_{{({site[0]},{site[1]})}}: \mathrm{{{physical['direction']}}}, d={physical['dim']}$" + ) + if show_leg_chargemaps: + phys_label += "\n" + rf"$q_p:$ {_format_compact_mapping(physical['chargemap'], max_items=4)}" + drawing.text( + (phys_xy[0] - 0.05, phys_xy[1] - 0.10), + phys_label, + fontsize=5.8, + ha="right", + va="top", + color=(0.18, 0.20, 0.23, 1.0), + ) + + if show_blocks and tensor["blocks"]: + blocks = tensor["blocks"][:max_blocks_per_site] + block_width = 0.15 + block_height = 0.12 + gap = 0.032 + total_width = len(blocks) * block_width + (len(blocks) - 1) * gap + start = x_pos - 0.5 * total_width + block_y = y_pos - 0.34 + for block_pos, block in enumerate(blocks): + bx = start + block_pos * (block_width + gap) + color = schematic.hash_to_color(str(block["sector"])) + drawing.ax.add_patch( + Rectangle( + (bx, block_y), + block_width, + block_height, + facecolor=color, + edgecolor=(0.16, 0.18, 0.21, 0.75), + alpha=0.86, + linewidth=0.60, + zorder=4, + ) + ) + if show_block_labels: + drawing.ax.text( + bx + 0.5 * block_width, + block_y - 0.04, + _format_sector(block["sector"]), + fontsize=4.8, + ha="center", + va="top", + rotation=45, + color=(0.15, 0.17, 0.20, 1.0), + ) + drawing.ax.text( + x_pos, + block_y + block_height + 0.022, + rf"$B={tensor['num_blocks']}$", + fontsize=5.8, + ha="center", + va="bottom", + color=(0.15, 0.17, 0.20, 1.0), + bbox={ + "boxstyle": "round,pad=0.05", + "facecolor": (1.0, 1.0, 1.0, 0.78), + "edgecolor": (1.0, 1.0, 1.0, 0.0), + "linewidth": 0.0, + }, + zorder=5, + ) + + xs = [xy[0] for xy in xy_by_site.values()] or [0.0] + ys = [xy[1] for xy in xy_by_site.values()] or [0.0] + right_pad = 0.78 + if show_diagnostics: + charge_line = _charge_summary_text(summary) + diagnostic_lines = [f"sites {summary['num_sites']}"] + if charge_line: + diagnostic_lines.append(charge_line) + diagnostic_lines += [ + f"max bond {summary['max_bond_dim']}", + f"bond sectors {summary['max_bond_sectors']}", + f"stored {summary['total_stored_size']}/{summary['total_dense_size']}", + f"density {summary['density']:.3f}", + ] + diagnostic = "\n".join(diagnostic_lines) + if show_blocks: + diagnostic += "\ncolored tiles: stored blocks" + if show_arrows: + diagnostic += "\narrows/labels: charge in/out flow" + if len(shown_tensors) < summary["num_sites"]: + diagnostic += f"\n+{summary['num_sites'] - len(shown_tensors)} sites hidden" + drawing.ax.text( + max(xs) + 0.82, + max(ys) + 0.46, + diagnostic, + fontsize=8, + ha="left", + va="top", + color=(0.15, 0.17, 0.20, 1.0), + bbox={ + "boxstyle": "round,pad=0.22", + "facecolor": (1.0, 1.0, 1.0, 0.92), + "edgecolor": (0.68, 0.70, 0.74, 1.0), + "linewidth": 0.8, + }, + ) + right_pad = 1.85 + elif len(shown_tensors) < summary["num_sites"]: + drawing.text( + (max(xs) + 0.55, min(ys)), + f"+{summary['num_sites'] - len(shown_tensors)}", + fontsize=9, + ha="left", + va="center", + color=(0.18, 0.20, 0.23, 1.0), + ) + right_pad = 1.10 + + if title is None: + title = "Symmray PEPS block structure" + drawing.ax.set_title(title) + charge_text = _charge_summary_text(summary) + if charge_text and not show_diagnostics: + drawing.text( + (min(xs) - 0.68, max(ys) + 0.58), + charge_text, + fontsize=8, + ha="left", + va="center", + color=(0.18, 0.20, 0.23, 1.0), + ) + drawing.ax.set_xlim(min(xs) - 0.82, max(xs) + right_pad) + y_min = min(ys) - (0.92 if show_phys_labels else 0.76) + y_max = max(ys) + 0.76 + drawing.ax.set_ylim(y_min, y_max) + drawing.ax.axis("off") + if return_summary: + return drawing, summary + return drawing + + def _site_parity(site): if isinstance(site, (tuple, list)): return sum(int(x) for x in site) % 2 @@ -360,7 +1974,7 @@ def _hamiltonian_from_edges(model, symmetry, edges, *, flat=False, **params): return sr.ham_tfim_from_edges(symmetry, edges, flat=flat, **params) if model == "heisenberg": return sr.ham_heisenberg_from_edges(symmetry, edges, flat=flat, **params) - if model == "fermi_hubbard": + if model in {"fermi_hubbard", "fermi_hubbard_u1u1"}: return sr.ham_fermi_hubbard_from_edges(symmetry, edges, flat=flat, **params) if model == "fermi_hubbard_spinless": return sr.ham_fermi_hubbard_spinless_from_edges(symmetry, edges, flat=flat, **params) @@ -524,6 +2138,21 @@ def overall_parity(self): """Return the configured total Z2 parity, i.e. charge sum modulo 2.""" return self.overall_charge(mod=2) + def fermionic_ordering(self): + """Return site and edge ordering metadata for this Symmray state. + + The metadata records the site order, stored edge order, local bond index + directions, and the direct-fermion tensor-network reference used by the + Symmray summaries. For fermionic states, this is the package-level + record of the graph/order data that Symmray uses for parity-aware + contractions. + """ + if self.site_ind_id == "k{}": + return symmray_mps_summary(self)["fermionic_ordering"] + if self.site_ind_id == "k{},{}": + return symmray_peps_summary(self)["fermionic_ordering"] + raise ValueError("Unsupported symmetric state site index convention.") + def operator_from_dense(self, array, *, charge=0, sectors=None, sites=None): """Convert a dense local observable/operator to this state's symmetry.""" sectors_use = self.phys_sectors if sectors is None else sectors @@ -689,7 +2318,7 @@ def apply_gates( self, gates, *, - contract="split", + contract="auto", max_bond=None, cutoff=1e-10, normalize=False, @@ -702,6 +2331,7 @@ def apply_gates( """Apply a bundled local gate stream to this state.""" target = self if inplace else self.copy() method = str(method).strip().lower() + contract_auto = contract is None or str(contract).strip().lower() == "auto" if max_bond is not None: compress_opts.setdefault("max_bond", max_bond) if cutoff is not None: @@ -711,7 +2341,8 @@ def apply_gates( from ..operators import gate as pepsy_gate opts = dict(compress_opts) - opts.setdefault("contract", contract) + if not contract_auto: + opts.setdefault("contract", contract) opts.update({} if gate_kwargs is None else dict(gate_kwargs)) target.network = pepsy_gate( target.network, @@ -752,7 +2383,7 @@ def apply_gates( target.network, gate, inds, - contract=contract, + contract="split" if contract_auto else contract, tags=[], info=None, inplace=True, @@ -775,7 +2406,7 @@ def time_evolve( max_bond=None, cutoff=1e-10, normalize=None, - contract="split", + contract="auto", inplace=True, method="direct", gauges=None, diff --git a/tests/test_core_seed.py b/tests/test_core_seed.py index 70a31d3..689d5ef 100644 --- a/tests/test_core_seed.py +++ b/tests/test_core_seed.py @@ -11,6 +11,46 @@ import quimb.tensor as qtn +def test_build_optimizer_loads_cotengrust_before_cotengra(monkeypatch): + """build_optimizer should make cotengrust available to cotengra.""" + events = [] + + def fake_ensure_cotengrust(): + events.append("cotengrust") + return object() + + class DummyOpt: + def __init__(self, **kwargs): + _ = kwargs + events.append("cotengra") + + monkeypatch.setattr(core, "_ensure_cotengrust", fake_ensure_cotengrust) + monkeypatch.setattr(core.ctg, "ReusableHyperOptimizer", DummyOpt) + + out = core.build_optimizer( + progbar=False, + parallel=False, + optlib="random", + directory=None, + max_repeats=1, + max_time="rate:1e2", + ) + + assert isinstance(out, DummyOpt) + assert events == ["cotengrust", "cotengra"] + + +def test_expec_mpo_mps_normalized_copy_does_not_divide_twice(): + """MPS expectations should be invariant to the input state's scale.""" + mps = qtn.MPS_rand_state(5, bond_dim=2, phys_dim=2, dtype="complex128", seed=26) + mps[0].modify(data=0.2 * mps[0].data) + mpo = core.id_to_mpo(5, dtype="complex128") + + measured = core.expec_mpo(mpo, mps, contraction_opt="auto-hq") + + assert abs(measured - 1.0) < 1e-12 + + def test_build_optimizer_constructs_without_seed(monkeypatch): """build_optimizer should configure cotengra optimizer without seed kwarg.""" captured = {} @@ -84,6 +124,34 @@ def __init__(self, *args, **kwargs): assert "seed" not in captured +def test_build_compressed_optimizer_loads_cotengrust_before_cotengra(monkeypatch): + """build_compressed_optimizer should make cotengrust available to cotengra.""" + events = [] + + def fake_ensure_cotengrust(): + events.append("cotengrust") + return object() + + class DummyCOpt: + def __init__(self, *args, **kwargs): + _ = args, kwargs + events.append("cotengra") + + monkeypatch.setattr(core, "_ensure_cotengrust", fake_ensure_cotengrust) + monkeypatch.setattr(core.ctg, "ReusableHyperCompressedOptimizer", DummyCOpt) + + out = core.build_compressed_optimizer( + progbar=False, + chi=4, + directory=None, + max_repeats=1, + max_time="rate:1e2", + ) + + assert isinstance(out, DummyCOpt) + assert events == ["cotengrust", "cotengra"] + + def test_tn_fidelity_uses_default_optimizer_settings(monkeypatch): """tn_fidelity should call build_optimizer with progbar disabled.""" captured = {} @@ -226,6 +294,7 @@ def test_default_backend_setters_roundtrip(): def test_core_all_exports_backend_jax(): """core.__all__ should include the JAX backend caster for parity.""" assert "backend_jax" in core.__all__ + assert "reg_rel_svd_torch" in core.__all__ assert "reg_complex_svd_torch" in core.__all__ assert "reg_complex_svd_jax" in core.__all__ @@ -807,6 +876,48 @@ def test_complex_svd_torch_backward_supports_batched_autoray(): assert torch.isfinite(A.grad).all() +def test_reg_rel_svd_torch_registers_stable_autoray_rule(): + """Relative SVD registration should stabilize degenerate spectra.""" + torch = pytest.importorskip("torch") + + core.reg_rel_svd_torch() + + A = torch.eye(3, dtype=torch.complex128, requires_grad=True) + U, S, Vh = ar.do("linalg.svd", A) + loss = S.sum() + (U[0, 0] * Vh.conj().T[0, 0].conj()).real + loss.backward() + + assert A.grad is not None + assert torch.isfinite(A.grad).all() + + +@pytest.mark.parametrize("dtype", ["float64", "complex128"]) +def test_reg_rel_svd_torch_forward_falls_back_to_scipy_gesvd(monkeypatch, dtype): + """Relative SVD should recover when torch's CPU SVD driver fails.""" + torch = pytest.importorskip("torch") + pytest.importorskip("scipy.linalg") + from pepsy.backends import linalg_torch as lrt # pylint: disable=import-outside-toplevel + + torch_dtype = getattr(torch, dtype) + + def fail_svd(*_args, **_kwargs): + raise RuntimeError("forced gesdd failure") + + monkeypatch.setattr(lrt.torch.linalg, "svd", fail_svd) + + torch.manual_seed(151617) + A = torch.randn(2, 4, 3, dtype=torch_dtype, requires_grad=True) + U, S, Vh = lrt.SVD.apply(A) + recon = (U * S.unsqueeze(-2)) @ Vh + + torch.testing.assert_close(recon, A.detach(), rtol=1e-10, atol=1e-10) + loss = recon.real.square().sum() + S.sum() + loss.backward() + + assert A.grad is not None + assert torch.isfinite(A.grad).all() + + def test_complex_svd_torch_backward_regularizes_degenerate_spectra(): """Custom complex SVD backward should stay finite at degenerate spectra.""" torch = pytest.importorskip("torch") diff --git a/tests/test_gate.py b/tests/test_gate.py index 5f2304c..cadd409 100644 --- a/tests/test_gate.py +++ b/tests/test_gate.py @@ -611,6 +611,38 @@ def _fake_bulk_helper(tn_i, G_arg, where=None, **kwargs): assert has_inplace_kw is (dim == "1d") +def test_gate_2d_3d_default_contract_is_reduce_split(monkeypatch): + """PEPS-like direct gate routing should default to quimb's reduce-split.""" + calls = [] + + def _fake_gate_2d(tn, G_arg, where=None, **kwargs): + calls.append(("2d", where, kwargs.get("contract"))) + return tn + + def _fake_gate_3d(tn, G_arg, where=None, **kwargs): + calls.append(("3d", where, kwargs.get("contract"))) + return tn + + monkeypatch.setattr("pepsy.operators.gates._apply_gate_2d", _fake_gate_2d) + monkeypatch.setattr("pepsy.operators.gates._apply_gate_3d", _fake_gate_3d) + + class _Dummy3DTN: # pylint: disable=too-few-public-methods + Lx = 1 + Ly = 1 + Lz = 2 + + gate_op = np.eye(4, dtype=np.complex128).reshape(2, 2, 2, 2) + peps = ps_to_peps(1, 2, dtype="complex128") + tn3d = _Dummy3DTN() + + assert apply_gate(peps, gate_op, ((0, 0), (0, 1))) is peps + assert apply_gate(tn3d, gate_op, ((0, 0, 0), (0, 0, 1))) is tn3d + assert calls == [ + ("2d", ((0, 0), (0, 1)), "reduce-split"), + ("3d", ((0, 0, 0), (0, 0, 1)), "reduce-split"), + ] + + @pytest.mark.parametrize("dim", ("1d", "2d", "3d")) def test_gate_sequence_dispatch_inplace_false_copies(monkeypatch, dim): """Bundled streams should copy TN first when ``inplace=False``.""" @@ -1071,6 +1103,7 @@ def _fake_gate(tn, gate, where, **kwargs): where, kwargs = calls[0] assert where == (0, 2) assert kwargs["max_bond"] == 5 + assert kwargs["contract"] == "reduce-split" assert "bond_dim" not in kwargs assert kwargs["ind_id"] == "k{}" @@ -1123,6 +1156,7 @@ def _fake_gate(tn, gate, where, **kwargs): where, kwargs = calls[0] assert where == ((0, 0), (1, 2)) assert kwargs["max_bond"] == 6 + assert kwargs["contract"] == "reduce-split" assert "bond_dim" not in kwargs assert kwargs["sequence"] == "auto" assert kwargs["ind_id"] == "k{},{}" diff --git a/tests/test_ham.py b/tests/test_ham.py index e8ed4c3..58d7382 100644 --- a/tests/test_ham.py +++ b/tests/test_ham.py @@ -1,5 +1,7 @@ """Regression tests for :mod:`pepsy.operators.hamiltonians` builders.""" +import builtins + import numpy as np import pytest import quimb @@ -536,6 +538,40 @@ def test_map_builder_instance_show_returns_schematic_drawing(): assert drawing.ax.get_title() == "Instance Mapper" +def test_map_builder_show_reports_missing_matplotlib(monkeypatch): + """show() should report the missing plotting package directly.""" + real_import = builtins.__import__ + + def import_without_matplotlib(name, *args, **kwargs): + if name == "matplotlib" or name.startswith("matplotlib."): + raise ImportError("No module named 'matplotlib'") + return real_import(name, *args, **kwargs) + + monkeypatch.setattr(builtins, "__import__", import_without_matplotlib) + + with pytest.raises(ImportError, match="matplotlib"): + OneDMap.show(2, 2, mode="snake") + + +def test_mpo_schematic_reports_missing_matplotlib(monkeypatch): + """MPO schematic rendering should report the missing plotting package directly.""" + real_import = builtins.__import__ + + def import_without_matplotlib(name, *args, **kwargs): + if name == "matplotlib" or name.startswith("matplotlib."): + raise ImportError("No module named 'matplotlib'") + return real_import(name, *args, **kwargs) + + class DummyMPO: + L = 4 + + monkeypatch.setattr(builtins, "__import__", import_without_matplotlib) + + builder = py.ham_tn(Lx=2, Ly=2) + with pytest.raises(ImportError, match="matplotlib"): + builder._show_mpo_schematic_2d(DummyMPO(), [((0, 0), (1, 0))]) + + def test_map_builder_show_rejects_3d(): """show() should fail clearly for 3D maps until a schematic 3D view exists.""" with pytest.raises(NotImplementedError, match="only available for 2D lattices"): diff --git a/tests/test_optimize_mps.py b/tests/test_optimize_mps.py index ff5ec65..7773b4c 100644 --- a/tests/test_optimize_mps.py +++ b/tests/test_optimize_mps.py @@ -183,6 +183,28 @@ def test_mps_optimizer_unitary_default_still_samples_norm_proxy(): assert len(opt.get_fidelities()) > 1 +def test_mps_optimizer_svd_forwards_cutoff_mode_to_final_compression(monkeypatch): + """SVD mode should honor explicit cutoff_mode in its chi compression pass.""" + p0 = qtn.MPS_computational_state("0000", dtype="complex128") + gates = [(qu.CNOT(), (0, 3))] + + opt = py.MpsOptimizer(p0.copy(), gates=gates, chi=2, mode="svd") + calls = [] + original_left_compress = opt.p.left_compress + + def _recording_left_compress(*args, **kwargs): + calls.append(dict(kwargs)) + return original_left_compress(*args, **kwargs) + + monkeypatch.setattr(opt.p, "left_compress", _recording_left_compress) + + opt.run(progbar=False, cutoff=1.0e-9, cutoff_mode="rsum2") + + assert calls + assert calls[-1]["cutoff"] == pytest.approx(1.0e-9) + assert calls[-1]["cutoff_mode"] == "rsum2" + + def test_mps_optimizer_rejects_negative_fidelity_samples(): """Negative fidelity_samples should fail clearly.""" p0 = qtn.MPS_computational_state("0000", dtype="complex128") diff --git a/tests/test_package_layout.py b/tests/test_package_layout.py index d3df532..fedb943 100644 --- a/tests/test_package_layout.py +++ b/tests/test_package_layout.py @@ -28,6 +28,7 @@ ps_to_3dpeps, ps_to_peps, reg_complex_svd_torch, + reg_rel_svd_torch, site_charge_from_occupations, ) @@ -58,6 +59,7 @@ def test_new_namespace_imports_resolve(): assert callable(haar_random_state) assert callable(ps_to_peps) assert callable(ps_to_3dpeps) + assert callable(reg_rel_svd_torch) assert callable(reg_complex_svd_torch) assert callable(site_charge_from_occupations) diff --git a/tests/test_public_api.py b/tests/test_public_api.py index 654820f..f81f909 100644 --- a/tests/test_public_api.py +++ b/tests/test_public_api.py @@ -27,10 +27,11 @@ def test_package_version_available(): "FDSolver", "MpsOptimizer", "MpoOptimizer", "PepsOptimizer", "PEPSSampleResult", "PepsBpSampler", "MpsSampler", "MpsSampleResult", "VecSampler", "gate", "tn_fidelity", "tn_norm", "SymGateStream", "SymHamiltonian", "SymMPS", "SymPEPS", - "default_physical_sectors", "sector_index_map", + "default_physical_sectors", "draw_symmray_blocks", "draw_symmray_mps", "draw_symmray_peps", + "sector_index_map", "site_charge_alternating", "site_charge_from_map", "site_charge_from_occupations", "site_charge_uniform", - "symm_operator_from_dense", + "symmray_block_summary", "symmray_mps_summary", "symmray_peps_summary", "symm_operator_from_dense", ] _EXPECTED_NOT_IN_ALL = [ @@ -41,7 +42,7 @@ def test_package_version_available(): "gates_tn_1d", "gates_tn_2d", "gates_tn_3d", "apply_2d_gate", "apply_2d_gates", "apply_2dtn_", "gate_2d", "gate_to_pepo", "gate_1d", "canonize_mps", "apply_gates_", "expec_TN_1D", "peps_I", - "reg_complex_svd_torch", "reg_complex_svd_jax", + "reg_rel_svd_torch", "reg_complex_svd_torch", "reg_complex_svd_jax", "reg_stop_gradient_torch", "stop_grad", "MPSOptimizer", "MPOOptimizer", "boundary_metrics", "boundary_states", "boundary_sweeps", "core", "fit", @@ -75,10 +76,11 @@ def test_internal_symbol_not_exported(name): "id_to_mpo", "id_to_pepo", "ps_to_pepo", "ps_to_mpo", "SweepOptimizer", "FDSolver", "MpsOptimizer", "MpoOptimizer", "PepsOptimizer", "PEPSSampleResult", "PepsBpSampler", "tn_fidelity", "tn_norm", "SymGateStream", "SymHamiltonian", "SymMPS", "SymPEPS", - "default_physical_sectors", "sector_index_map", + "default_physical_sectors", "draw_symmray_blocks", "draw_symmray_mps", "draw_symmray_peps", + "sector_index_map", "site_charge_alternating", "site_charge_from_map", "site_charge_from_occupations", "site_charge_uniform", - "symm_operator_from_dense", + "symmray_block_summary", "symmray_mps_summary", "symmray_peps_summary", "symm_operator_from_dense", ] _BLOCKED_NAMES = _EXPECTED_NOT_IN_ALL @@ -130,6 +132,8 @@ def test_optional_linalg_registrations_resolve(): if has_torch: import torch + assert callable(pepsy.tensors.core.reg_rel_svd_torch) + assert callable(pepsy.tensors.reg_rel_svd_torch) assert callable(pepsy.tensors.core.reg_complex_svd_torch) assert callable(pepsy.tensors.reg_complex_svd_torch) x = torch.tensor([1.0], dtype=torch.float64, requires_grad=True) diff --git a/tests/test_symmetric_tensors.py b/tests/test_symmetric_tensors.py index 6b60e39..f502bb2 100644 --- a/tests/test_symmetric_tensors.py +++ b/tests/test_symmetric_tensors.py @@ -11,11 +11,17 @@ SymMPS, SymPEPS, default_physical_sectors, + draw_symmray_blocks, + draw_symmray_mps, + draw_symmray_peps, sector_index_map, site_charge_alternating, site_charge_from_map, site_charge_from_occupations, site_charge_uniform, + symmray_block_summary, + symmray_mps_summary, + symmray_peps_summary, symm_operator_from_dense, ) @@ -23,10 +29,339 @@ sr = pytest.importorskip("symmray") +def test_symmray_block_summary_and_schematic_for_z2_gate(): + """Symmray gate helpers should expose and draw block-sector structure.""" + state = SymMPS.random( + 4, + symmetry="Z2", + phys_dim={0: 1, 1: 1}, + site_charge=site_charge_from_occupations([0] * 4), + bond_dim=2, + seed=1, + dtype="complex128", + ) + rz_gate = state.operator_from_dense(pepsy.rz(0.1), charge=0, sites=1) + + summary = symmray_block_summary(rz_gate) + + assert summary["shape"] == (2, 2) + assert summary["num_blocks"] == 2 + assert summary["stored_size"] == 2 + assert summary["dense_size"] == 4 + assert [block["sector"] for block in summary["blocks"]] == [(0, 0), (1, 1)] + assert [block["shape"] for block in summary["blocks"]] == [(1, 1), (1, 1)] + assert summary["indices"][0]["direction"] == "out" + assert summary["indices"][1]["direction"] == "in" + + pytest.importorskip("matplotlib") + drawing, drawn_summary = draw_symmray_blocks( + rz_gate, + title="Z2 RZ gate", + return_summary=True, + ) + assert drawing is not None + assert drawn_summary["blocks"] == summary["blocks"] + + +def test_symmray_mps_summary_and_schematic_for_z2_chain(): + """Symmray MPS drawings should expose site blocks and bond-sector metadata.""" + state = SymMPS.random( + 4, + symmetry="Z2", + phys_dim={0: 1, 1: 1}, + site_charge=site_charge_from_occupations([0] * 4), + bond_dim=2, + seed=2, + dtype="complex128", + ) + + summary = symmray_mps_summary(state.tn) + + assert summary["num_sites"] == 4 + assert summary["max_bond_dim"] == 2 + assert summary["max_bond_sectors"] == 2 + assert summary["total_stored_size"] == 12 + assert summary["total_dense_size"] == 24 + assert summary["charge_total"] == summary["total_charge"] + assert summary["Q_total"] == summary["total_parity"] + assert summary["tensors"][0]["site"] == 0 + assert summary["tensors"][0]["physical"]["chargemap"] == {0: 1, 1: 1} + assert summary["tensors"][0]["physical"]["direction"] in {"in", "out"} + assert summary["tensors"][1]["left_bond"]["between"] == (0, 1) + assert summary["bonds"][0]["chargemap"] == {0: 1, 1: 1} + assert summary["bonds"][0]["left_direction"] in {"in", "out"} + assert summary["bonds"][0]["right_direction"] in {"in", "out"} + + pytest.importorskip("matplotlib") + drawing, drawn_summary = draw_symmray_mps( + state.tn, + title="Z2 MPS", + return_summary=True, + ) + assert hasattr(drawing, "fig") + assert hasattr(drawing, "ax") + assert drawn_summary["bonds"] == summary["bonds"] + + +def test_symmray_peps_summary_and_schematic_for_z2_grid(): + """Symmray PEPS drawings should expose grid bonds and block-sector metadata.""" + state = SymPEPS.random( + 2, + 2, + symmetry="Z2", + phys_dim={0: 1, 1: 1}, + site_charge=site_charge_from_occupations({(i, j): 0 for i in range(2) for j in range(2)}), + bond_dim=2, + seed=3, + dtype="complex128", + ) + + summary = symmray_peps_summary(state) + + assert summary["Lx"] == 2 + assert summary["Ly"] == 2 + assert summary["num_sites"] == 4 + assert len(summary["bonds"]) == 4 + assert summary["max_bond_dim"] == 2 + assert summary["max_bond_sectors"] == 2 + assert summary["total_stored_size"] == 16 + assert summary["total_dense_size"] == 32 + assert summary["charge_total"] == summary["total_charge"] + assert summary["Q_total"] == summary["total_parity"] + assert summary["tensors"][0]["site"] == (0, 0) + assert summary["tensors"][0]["physical"]["chargemap"] == {0: 1, 1: 1} + assert summary["tensors"][0]["physical"]["direction"] in {"in", "out"} + assert summary["tensors"][0]["bonds"]["right"]["between"] == ((0, 0), (0, 1)) + assert summary["tensors"][0]["bonds"]["down"]["between"] == ((0, 0), (1, 0)) + assert summary["bonds"][0]["site_a_direction"] in {"in", "out"} + assert summary["bonds"][0]["site_b_direction"] in {"in", "out"} + + pytest.importorskip("matplotlib") + drawing, drawn_summary = draw_symmray_peps( + state, + title="Z2 PEPS", + return_summary=True, + ) + assert hasattr(drawing, "fig") + assert hasattr(drawing, "ax") + assert drawn_summary["bonds"] == summary["bonds"] + + +def _square_lattice_edges(Lx, Ly): + """Return nearest-neighbor square-lattice edges in row-major MPS order.""" + edges = [] + + def site(x, y): + return x * Ly + y + + for x in range(Lx): + for y in range(Ly): + if x + 1 < Lx: + edges.append((site(x, y), site(x + 1, y))) + if y + 1 < Ly: + edges.append((site(x, y), site(x, y + 1))) + return tuple(edges) + + +def _square_lattice_coordinate_edges(Lx, Ly): + """Return nearest-neighbor square-lattice edges as PEPS coordinates.""" + edges = [] + for x in range(Lx): + for y in range(Ly): + if y + 1 < Ly: + edges.append(((x, y), (x, y + 1))) + if x + 1 < Lx: + edges.append(((x, y), (x + 1, y))) + return tuple(edges) + + +def _xy_u1_hamiltonian(edges): + """Build the U(1)-symmetric XY Hamiltonian from an explicit dense term.""" + sx = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=np.complex128) + sy = np.array([[0.0, -1.0j], [1.0j, 0.0]], dtype=np.complex128) + xy_term_dense = 0.5 * (np.kron(sx, sx) + np.kron(sy, sy)) + xy_term = symm_operator_from_dense( + xy_term_dense, + {0: 1, 1: 1}, + symmetry="U1", + charge=0, + sites=2, + ) + return SymHamiltonian( + model="xy", + symmetry="U1", + edges=tuple(edges), + terms={edge: xy_term for edge in edges}, + parameters={"J": 1.0}, + ) + + +def _all_tensor_data_symmray(tn): + """Return whether every tensor stores Symmray block-sparse data.""" + return all( + hasattr(tensor.data, "blocks") and hasattr(tensor.data, "indices") + for tensor in tn.tensors + ) + + +def _finite_double_layer_norm(tn): + """Return whether the direct double-layer norm contraction is finite.""" + norm = (tn.H & tn).contract(all, optimize="auto-hq") + return np.isfinite(np.real(norm)) + + +def _raw_mps_norm(mps): + """Return an MPS norm with PEPSY's stored exponent removed.""" + raw = mps.copy() + if hasattr(raw, "exponent"): + raw.exponent = 0.0 + return raw.norm() + + +def _build_3x3_symmray_mps_case(name): + """Build an explicit 3x3 state, Hamiltonian, and gate stream.""" + edges = _square_lattice_edges(3, 3) + + if name == "itf_z2": + state = SymMPS.random( + 9, + symmetry="Z2", + phys_dim={0: 1, 1: 1}, + site_charge=site_charge_from_occupations([0] * 9), + bond_dim=4, + seed=41, + dtype="complex128", + ) + hamiltonian = SymHamiltonian.from_edges( + "itf", + "Z2", + edges, + jx=-1.0, + hz=-0.5, + ) + gates = hamiltonian.gate_stream(0.001, imaginary=False) + return state, hamiltonian, gates, False, 0 + + if name == "xy_u1": + occupations = [1, 0, 1, 0, 1, 0, 1, 0, 1] + state = SymMPS.random( + 9, + symmetry="U1", + phys_dim={0: 1, 1: 1}, + site_charge=site_charge_from_occupations(occupations), + bond_dim=4, + seed=42, + dtype="complex128", + ) + hamiltonian = _xy_u1_hamiltonian(edges) + gates = hamiltonian.gate_stream(0.001, imaginary=False) + return state, hamiltonian, gates, False, sum(occupations) + + if name == "fermi_hubbard_u1": + occupations = [1] * 9 + state = SymMPS.random( + 9, + symmetry="U1", + phys_dim={0: 1, 1: 2, 2: 1}, + fermionic=True, + site_charge=site_charge_from_occupations(occupations), + bond_dim=4, + seed=43, + dtype="complex128", + ) + hamiltonian = SymHamiltonian.from_edges( + "fermi_hubbard", + "U1", + edges, + t=1.0, + U=2.0, + mu=0.1, + ) + gates = hamiltonian.gate_stream(0.0005, imaginary=True) + return state, hamiltonian, gates, True, sum(occupations) + + raise ValueError(f"Unknown symmetric MPS case {name!r}.") + + +def _build_3x3_symmray_peps_case(name, *, edges=None): + """Build an explicit 3x3 PEPS state, Hamiltonian, and gate stream.""" + edges = _square_lattice_coordinate_edges(3, 3) if edges is None else tuple(edges) + + if name == "itf_z2": + charges = {(i, j): 0 for i in range(3) for j in range(3)} + state = SymPEPS.random( + 3, + 3, + symmetry="Z2", + phys_dim={0: 1, 1: 1}, + site_charge=site_charge_from_occupations(charges), + bond_dim=2, + seed=51, + dtype="complex128", + ) + hamiltonian = SymHamiltonian.from_edges( + "itf", + "Z2", + edges, + jx=-1.0, + hz=-0.5, + ) + gates = hamiltonian.gate_stream(0.001, imaginary=False) + return state, hamiltonian, gates, 0 + + if name == "xy_u1": + charges = {(i, j): (i + j) % 2 for i in range(3) for j in range(3)} + state = SymPEPS.random( + 3, + 3, + symmetry="U1", + phys_dim={0: 1, 1: 1}, + site_charge=site_charge_from_occupations(charges), + bond_dim=2, + seed=52, + dtype="complex128", + ) + hamiltonian = _xy_u1_hamiltonian(edges) + gates = hamiltonian.gate_stream(0.001, imaginary=False) + return state, hamiltonian, gates, sum(charges.values()) + + if name == "fermi_hubbard_u1": + charges = {(i, j): 1 for i in range(3) for j in range(3)} + state = SymPEPS.random( + 3, + 3, + symmetry="U1", + phys_dim={0: 1, 1: 2, 2: 1}, + fermionic=True, + site_charge=site_charge_from_occupations(charges), + bond_dim=2, + seed=53, + dtype="complex128", + ) + hamiltonian = SymHamiltonian.from_edges( + "fermi_hubbard", + "U1", + edges, + t=1.0, + U=2.0, + mu=0.1, + ) + gates = hamiltonian.gate_stream(0.0005, imaginary=True) + return state, hamiltonian, gates, sum(charges.values()) + + raise ValueError(f"Unknown symmetric PEPS case {name!r}.") + + def test_sector_and_charge_helpers_make_total_charge_explicit(): """Physical sectors and local charges should be easy to inspect.""" assert default_physical_sectors("U1", 4) == {0: 1, 1: 2, 2: 1} assert default_physical_sectors(model="fermi_hubbard") == {0: 1, 1: 2, 2: 1} + assert default_physical_sectors(model="fermi_hubbard_u1u1") == { + (0, 0): 1, + (0, 1): 1, + (1, 0): 1, + (1, 1): 1, + } assert sector_index_map({0: 1, 1: 2, 2: 1}) == {0: 0, 1: 1, 2: 1, 3: 2} occupations = [1, 0, 1, 0] @@ -138,6 +473,95 @@ def test_symmps_fermi_hubbard_defaults_to_fermionic_u1(): assert evolved.norm() == pytest.approx(1.0) +def test_fermi_hubbard_u1u1_preset_uses_spin_resolved_fermionic_tensors(): + """U1U1 preset should expose spin-resolved spinful Hubbard sectors.""" + sectors = default_physical_sectors(model="fermi_hubbard_u1u1") + mps = SymMPS.for_model( + "fermi_hubbard_u1u1", + 4, + bond_dim=2, + site_charge=site_charge_from_occupations([(1, 0), (0, 1), (1, 0), (0, 1)]), + seed=3, + dtype="complex128", + ) + ham_mps = mps.build_hamiltonian(t=1.0, U=4.0, mu=0.0) + evolved_mps = mps.time_evolve( + 0.001, + steps=1, + hamiltonian=ham_mps, + imaginary=True, + max_bond=4, + inplace=False, + ) + + assert mps.symmetry == "U1U1" + assert mps.fermionic + assert mps.phys_sectors == sectors + assert mps.overall_charge() == (2, 2) + assert all(type(term).__name__ == "U1U1FermionicArray" for term in ham_mps.terms.values()) + mps_ordering = mps.fermionic_ordering() + assert mps_ordering["enabled"] is True + assert mps_ordering["network_kind"] == "mps" + assert mps_ordering["methods_reference"]["doi"] == "10.1103/PhysRevResearch.7.023193" + assert mps_ordering["site_order"] == (0, 1, 2, 3) + assert mps_ordering["edge_order"] == mps.edges + assert mps_ordering["edges"][0]["edge"] == (0, 1) + assert mps_ordering["edges"][0]["edge_order"] == 0 + assert mps_ordering["edges"][0]["index_directions"][0]["site"] == 0 + assert mps_ordering["edges"][0]["index_directions"][0]["direction"] in {"in", "out"} + assert evolved_mps.overall_charge() == (2, 2) + assert evolved_mps.tn.max_bond() <= 4 + assert evolved_mps.norm() == pytest.approx(1.0) + + peps_charges = { + (0, 0): (1, 0), + (0, 1): (0, 1), + (1, 0): (1, 0), + (1, 1): (0, 1), + } + peps = SymPEPS.for_model( + "fermi_hubbard_u1u1", + 2, + 2, + bond_dim=2, + site_charge=site_charge_from_occupations(peps_charges), + seed=4, + dtype="complex128", + ) + ham_peps = peps.build_hamiltonian(t=1.0, U=4.0, mu=0.0) + evolved_peps = peps.time_evolve( + 0.001, + steps=1, + hamiltonian=ham_peps, + imaginary=True, + max_bond=4, + method="gate", + inplace=False, + ) + + assert peps.symmetry == "U1U1" + assert peps.fermionic + assert peps.phys_sectors == sectors + assert peps.overall_charge() == (2, 2) + assert all(type(term).__name__ == "U1U1FermionicArray" for term in ham_peps.terms.values()) + peps_ordering = peps.fermionic_ordering() + assert peps_ordering["enabled"] is True + assert peps_ordering["network_kind"] == "peps" + assert peps_ordering["methods_reference"]["doi"] == "10.1103/PhysRevResearch.7.023193" + assert peps_ordering["site_order"] == ((0, 0), (0, 1), (1, 0), (1, 1)) + assert peps_ordering["edge_order"] == peps.edges + peps_first_edge = peps.edges[0] + peps_first_record = next( + record for record in peps_ordering["edges"] + if record["edge"] == peps_first_edge + ) + assert peps_first_record["edge_order"] == 0 + assert tuple(item["site"] for item in peps_first_record["index_directions"]) == peps_first_edge + assert peps_first_record["index_directions"][0]["direction"] in {"in", "out"} + assert evolved_peps.overall_charge() == (2, 2) + assert evolved_peps.tn.max_bond() <= 4 + + def test_symmps_gate_stream_runs_mps_optimizer_mpo_heisenberg(): """Symmray U(1) gates should run through MpsOptimizer(mode='mpo').""" state = SymMPS.for_model( @@ -218,6 +642,146 @@ def test_symmps_mps_optimizer_coerces_dense_hamiltonian_terms(): assert np.isfinite(np.real(evolved.norm())) +@pytest.mark.parametrize("case_name", ["itf_z2", "xy_u1", "fermi_hubbard_u1"]) +@pytest.mark.parametrize("mode", ["dmrg", "mpo", "swap", "svd", "exact"]) +def test_symmps_mps_optimizer_3x3_streams_cover_supported_modes(case_name, mode): + """Explicit 3x3 Symmray streams should run through supported MPS modes.""" + state, hamiltonian, gates, non_unitary, expected_charge = _build_3x3_symmray_mps_case( + case_name + ) + + assert len(hamiltonian.edges) == 12 + assert len(gates) == 12 + assert state.L == 9 + assert state.overall_charge() == expected_charge + valid_gate_shapes = {(2, 2, 2, 2), (4, 4, 4, 4)} + assert all(term.shape in valid_gate_shapes for term in hamiltonian.terms.values()) + assert all(gate.shape in valid_gate_shapes for gate, _ in gates) + + opt = pepsy.MpsOptimizer(state.tn.copy(), gates, chi=4, mode=mode) + run_kwargs = { + "progbar": False, + "cutoff": 1.0e-10, + "fidelity_samples": 0, + "n_iter": 4, + } + if non_unitary and mode != "exact": + run_kwargs.update( + { + "non_unitary": True, + "normalize_every": 1, + "normalize_final": True, + } + ) + + out = opt.run(**run_kwargs) + + assert _all_tensor_data_symmray(out) + assert _finite_double_layer_norm(out) + + if mode == "exact": + assert len(out.tensors) == 1 + else: + assert out.L == 9 + assert out.max_bond() <= 4 + + if non_unitary and mode != "exact": + assert len(opt.get_normalizations()) == len(gates) + assert _raw_mps_norm(out) == pytest.approx(1.0) + else: + assert opt.get_normalizations() == [] + + +def test_symmps_mps_optimizer_symmray_dmrg_expansion_caveat_is_explicit(): + """DMRG should fail clearly when Symmray bond expansion would be required.""" + state = SymMPS.random( + 4, + symmetry="Z2", + phys_dim={0: 1, 1: 1}, + site_charge=site_charge_from_occupations([0] * 4), + bond_dim=2, + seed=46, + dtype="complex128", + ) + nearest_ham = SymHamiltonian.from_edges( + "itf", + "Z2", + [(0, 1)], + jx=-1.0, + hz=-0.5, + ) + nonlocal_ham = SymHamiltonian.from_edges( + "itf", + "Z2", + [(0, 2)], + jx=-1.0, + hz=-0.5, + ) + + nearest_gates = nearest_ham.gate_stream(0.001) + nonlocal_gates = nonlocal_ham.gate_stream(0.001) + + for mode, gates in [ + ("mpo", nearest_gates), + ("mpo", nonlocal_gates), + ("svd", nonlocal_gates), + ]: + out = pepsy.MpsOptimizer( + state.tn.copy(), + gates, + chi=2, + mode=mode, + ).run(progbar=False, cutoff=1.0e-10, fidelity_samples=0) + assert out.L == 4 + assert out.max_bond() <= 2 + + with pytest.raises(ValueError, match="bond_dim >= chi"): + pepsy.MpsOptimizer(state.tn.copy(), nearest_gates, chi=4, mode="dmrg").run( + progbar=False + ) + + +@pytest.mark.parametrize("mode", ["mpo", "svd"]) +def test_symmps_mps_optimizer_symmray_auto_swap_tracks_infidelity(mode): + """Symmray auto-swap fallbacks should still support true infidelity samples.""" + state = SymMPS.random( + 4, + symmetry="Z2", + phys_dim={0: 1, 1: 1}, + site_charge=site_charge_from_occupations([0] * 4), + bond_dim=2, + seed=47, + dtype="complex128", + ) + hamiltonian = SymHamiltonian.from_edges( + "itf", + "Z2", + [(0, 2)], + jx=-1.0, + hz=-0.5, + ) + + opt = pepsy.MpsOptimizer( + state.tn.copy(), + hamiltonian.gate_stream(0.001), + chi=2, + mode=mode, + ) + out = opt.run( + progbar=False, + cutoff=1.0e-10, + fidelity_samples=0, + track_infidelity=True, + ) + + samples = opt.get_infidelity_samples() + assert out.L == 4 + assert out.max_bond() <= 2 + assert len(samples) == 1 + assert 0.0 <= samples[0]["fidelity"] <= 1.0 + assert 0.0 <= opt.get_infidelities()[-1] <= 1.0 + + def test_sympeps_tfim_builds_z2_terms_and_step(): """SymPEPS should build Z2 TFIM terms on a square grid.""" state = SymPEPS.for_model( @@ -388,6 +952,107 @@ def test_sympeps_gate_stream_runs_pepsy_gate_and_gate_simple(): assert np.isfinite(np.real(out_simple.norm())) +@pytest.mark.parametrize("case_name", ["itf_z2", "xy_u1", "fermi_hubbard_u1"]) +@pytest.mark.parametrize("method", ["gate", "simple"]) +def test_sympeps_gate_wrappers_3x3_streams_cover_symmetries(case_name, method): + """PEPSY gate wrappers should handle 3x3 Symmray PEPS gate streams.""" + state, hamiltonian, gates, expected_charge = _build_3x3_symmray_peps_case( + case_name + ) + + assert len(hamiltonian.edges) == 12 + assert len(gates) == 12 + assert state.Lx == 3 + assert state.Ly == 3 + assert state.overall_charge() == expected_charge + + if method == "gate": + out = gate( + state.tn.copy(), + gates, + max_bond=4, + cutoff=1.0e-10, + inplace=False, + ) + else: + gauges = {} + out = gate_simple( + state.tn.copy(), + gates, + gauges=gauges, + max_bond=4, + cutoff=1.0e-10, + inplace=False, + ) + assert len(gauges) == len(gates) + + assert out.Lx == 3 + assert out.Ly == 3 + assert out.max_bond() <= 4 + assert _all_tensor_data_symmray(out) + + +@pytest.mark.parametrize("case_name", ["itf_z2", "xy_u1", "fermi_hubbard_u1"]) +@pytest.mark.parametrize("method", ["gate", "simple"]) +def test_sympeps_gate_wrappers_route_nonlocal_symmray_swaps(case_name, method): + """Internal routed SWAPs should be Symmray arrays for nonlocal PEPS gates.""" + nonlocal_edge = (((0, 0), (2, 2)),) + state, _, gates, _ = _build_3x3_symmray_peps_case( + case_name, + edges=nonlocal_edge, + ) + + if method == "gate": + out = gate( + state.tn.copy(), + gates, + max_bond=4, + cutoff=1.0e-10, + inplace=False, + ) + else: + gauges = {} + out = gate_simple( + state.tn.copy(), + gates, + gauges=gauges, + max_bond=4, + cutoff=1.0e-10, + inplace=False, + ) + assert len(gauges) > 0 + + assert out.max_bond() <= 4 + assert _all_tensor_data_symmray(out) + + +def test_sympeps_gate_method_preserves_pepsy_gate_contract_default(monkeypatch): + """SymPEPS method='gate' should not override pepsy.gate's default.""" + state = SymPEPS.for_model( + "heisenberg", + 2, + 2, + bond_dim=2, + seed=12, + dtype="complex128", + ) + calls = [] + + def _fake_gate(tn, gates, **kwargs): + calls.append((gates, kwargs.copy())) + return tn + + monkeypatch.setattr("pepsy.operators.gate", _fake_gate) + + out = state.copy().apply_gates( + ((np.eye(2, dtype=np.complex128), ((0, 0),)),), + method="gate", + ) + + assert out.tn is not None + assert "contract" not in calls[0][1] + + def test_sympeps_raw_pepsy_gate_functions_accept_symmray_streams(): """The plain gate functions should accept a SymGateStream directly.""" state = SymPEPS.for_model("itf", 2, 2, bond_dim=2, seed=10, dtype="complex128")