Skip to content

Refactor to use PEPSKit.jl for certain computations. - #1

Draft
jack-dunham wants to merge 64 commits into
devfrom
pepskit-refactor
Draft

Refactor to use PEPSKit.jl for certain computations. #1
jack-dunham wants to merge 64 commits into
devfrom
pepskit-refactor

Conversation

@jack-dunham

Copy link
Copy Markdown
Owner

This is to make this package lighter and less bug-prone.

jack-dunham and others added 30 commits August 11, 2026 11:21
Design documents for retiring the parts of this package that PEPSKit 0.8.1
now implements: InfinitePEPO and SUWeight in place of PEPO/UnitCellBonds,
its SimpleUpdate for thermal annealing, FullEnvTruncation/ALSTruncation and
bondenv_ctm in place of the local FET and isometry code, SUGauge/BPGauge in
place of gauge/ and dev/vidalgauge.jl, and CTMRGEnv/leading_boundary in place
of contract/ and the TRGroups Renormalization glue.

The Lindblad update stays local: PEPSKit applies PEPO gates as exp(-Hdt/2) rho
exp(-Hdt/2) and asserts that real-time PEPO evolution is unimplemented, so a
vectorised Liouvillian acting jointly on ket and bra cannot use its _su_iter!.
The plan builds that one step on PEPSKit primitives instead.

Prerequisite recorded in both documents: PEPSKit needs TensorKit 0.16.5/0.17
while this package and TensorRenormalizationGroups are on 0.14.11, so nothing
can be outsourced before that migration.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Seven benchmarks compare the package against dense linear algebra rather than
against its own output, so a disagreement locates an error instead of a change.
Four agree to machine precision: one-local Lindblad dynamics against exp(t L)
(2.5e-15), trace/hermiticity/positivity of the one-site RDM (4.4e-16), and the
high-temperature limit of thermal annealing (2.1e-4, with beta = t).

Three disagree:

- The Hamiltonian in a Liouvillian enters as +i[H,rho] rather than -i[H,rho].
  Reversing the sign of H reproduces the package to 4.4e-16 for one spin with
  H = (h/2) Z at D=1 with no dissipation, and to 6.5e-6 for a two-body term.
  thermalvectorise is unaffected, so the error is in vectorise.
- pepositedata (src/pepo/pepo.jl:100-114) prepares rho^4 instead of rho: the
  input is written into all four directional tensors and contracted in a chain
  through the physical index. Matches rho^4 to 1.7e-16 for three mixed states.
  Pure states and ThermalState are unaffected, which is every existing example.
- Bond dimension above the rank the Trotter gate generates gives garbage, since
  simple update inverts the exact zeros in the padded bond weights: a one-local
  model exact at D=1,2 deviates by 2.0 at D=4 and 72.5 at D=8.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
tebdtensor built the vectorised superoperator as a matrix whose composite
index runs ket-slowest (row major), then moved it into the eight-leg array
with copyto!, which fills the fastest axis first. That reversed each group of
four legs, exchanging ket with bra and reversing the two sites, so evolution
ran with drho/dt = +i[H,rho] and with D[conj(L)] in place of D[L]. The site
reversal was invisible because every term in the tests and documentation is
symmetric under exchanging the two sites.

The ordering to obey is PEPSKit's, PEPOTensor : P (x) P' <- N (x) E (x) S (x) W,
with index 1 outgoing (ket, left-multiplied) and index 2 incoming (bra,
right-multiplied); trace_physicalspaces contracts O[p p; N E S W] and
_apply_sitegate applies gate * a to index 1 and a * gate to index 2. The
vectorisations already matched that convention, so only the conversion changes:
reshape into the reversed leg order, then permutedims back.

B6 (Larmor precession, D=1, no dissipation) moves from 1.13 to 4.4e-16 and its
Hamiltonian-reversed figure from 4.4e-16 to 1.13; B2 from 2.30 to 2.8e-6; a
jump operator with conj(L) != -L moves from matching D[conj(L)] to matching
D[L] at 3.9e-15. B1, B3 and B4 are unchanged, since their operators are
insensitive to the reversal - which is why this survived.

tebdtensor(model::Ising, ...) in heisenberg.jl reshapes a one-site
superoperator the same way and needs its own benchmark before being changed.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
tebdtensor(model::Ising, p1, p2, dt) reshaped `L * op * L` leg-1-fastest at
heisenberg.jl:207, the same reversal already fixed in the generic path, so the
documented Ising fast path also evolved with the ket and bra exchanged.

Both conversions now go through `vectorised_tensormap`, which reverses each
group of legs before wrapping the matrix, rather than repeating the permutation
at each call site.

B8 covers this path: at zero coupling the Ising gate reduces to a pure local
exponential, so it must build the same operator as the generic path, whose
absolute correctness B1 and B6 pin against dense linear algebra. B8 moves from
1.13 to 8.9e-16. With the coupling restored the two paths differ by 4.7e-4 at
dt = 0.05 and 1.2e-4 at dt = 0.025, a ratio of 3.99, which is the expected
second-order Trotter difference between exponentiating the sum and splitting
the local terms symmetrically around the interaction.

B1 through B7 are unchanged. Two conversions keep the defect, both unreachable:
tebdtensor(model::Ising, p, method, dt) at heisenberg.jl:222 returns
TensorMap(data, ...) with `data` never assigned, and the AbstractMatrix method
at dynamicalmap.jl:401 is reached only through `timestepper`, which has no live
callers.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
B1 through B8 all use terms invariant under exchanging the two sites of a bond,
so none of them constrained the site order that the leg-reversal fix also
corrected. B9 uses a X_j Y_k and reproduces the exact generator to 6.6e-5 while
deviating by 1.60 from the reference for the exchanged term a Y_j X_k.

Getting that power took two attempts, both recorded in the docstrings. On a
uniform product state the two orderings are indistinguishable: summing a
one-site observable over the four bonds places the site in the first factor
twice and the second factor twice, so exchanging the factors permutes the sum.
A checkerboard does not help either, since both neighbours along an axis hold
the same state. Exchanging the factors of every bond term is a lattice
inversion, so it is only observable when the state breaks inversion symmetry;
B9 therefore uses a three-period stripe, where a site has B on one side and C
on the other.

This also bounds the severity of that half of the earlier fix: the ket-bra
exchange was a physics error, whereas the site reversal was a relabelling for
every model the package has been run on.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
TRGroups stays on TensorKit 0.14, so raising the bound breaks it while it is
still a dependency: task 0.3 and Phase 1 merge into one stretch with no working
build. That makes deletion-first ordering necessary, since porting code the plan
later deletes would be work done blind and then discarded, and it makes the
benchmark suite the exit criterion.

Dependency resolution checked against the local PEPSKit clone: TensorKit 0.17.1,
MPSKit 0.13.13, MatrixAlgebraKit 0.6.9, KrylovKit 0.10.4, TensorOperations 5.7.0,
Zygote 0.7.12 on Julia 1.12.5. KrylovKit moving 0.9.5 to 0.10.4 is a second
migration surface the plan had not recorded.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
PEPSKit has no minimal canonical form. Its gauges are SUGauge, which fixes the
gauge by trivial simple update, and BPGauge, which diagonalises the local belief
propagation messages; the third gauge_fix, in algorithms/ctmrg/gaugefix.jl,
fixes the CTMRG environment's gauge for fixed-point differentiation. None of
these is the gradient flow that minimalcanonicalform! implements, so the earlier
mapping of MCF onto BPGauge was wrong and is corrected.

The file moves to src/dev/mcf.jl and carries the MCF struct with it so it is
self-contained for a later rewrite onto InfinitePEPO. Two defects to repair at
that point: it throws for unit cells larger than 1x1, and fixgauge!(x, ::MCF)
takes two arguments where its only callers in truncation/isometries.jl pass
three, so selecting SimpleUpdate(gauge = MCF()) raises a MethodError instead of
gauging anything. Its one working entry point is IterativeSU with gauge == -1.

The package still loads and MCF still resolves; nothing else changes.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Records the rule that no code leaves the repository: whatever PEPSKit does not
provide moves to src/dev with its include commented out, and whatever PEPSKit
does provide moves there too, ceasing to be loaded only once its replacement
demonstrably works.

The unit of work is a subsystem rather than a file, because most files mix the
two categories - pepo/weights.jl holds setweight! alongside bondentropy, and
truncation/truncation.jl holds the FET and SVDU structs alongside NTU - so files
are parked whole as each replacement lands.

Also flags two mappings in the API table as unverified, so they cannot drive a
parking decision yet. WTG is not a state gauge: fixgauge!(iso, metric, ::WTG)
gauges the truncation isometries, so its analogue is bondenv/gaugefix.jl rather
than SUGauge. And the product-state constructors have no PEPSKit builder beyond
infinite_temperature_density_matrix, which covers only the identity.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Three of its statements had gone stale. The TRGroups migration is declined, so
the stretch with no working build is a decision rather than an open question.
Phase 5 parks what it does not rebuild instead of deleting it, and the NTU call
no longer reads as delete-or-port. Phase completion is the exact benchmark suite
rather than recorded trajectories, which would have preserved the leg reversal
that Phase 0 found. Benchmark counts updated from seven to nine.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Swaps the dependency set (TensorKit 0.17.1, PEPSKit 0.8.1, KrylovKit 0.10.4,
no TensorRenormalizationGroups) and replaces the state and observable layers.
Everything else is unloaded rather than removed: the includes are commented out
in the module, awaiting their replacements, and no file is deleted.

src/state.jl builds a product-state density-matrix PEPO on PEPSKit's PEPOTensor
convention, P (x) P' <- N (x) E (x) S (x) W, index 1 outgoing (ket) and index 2
incoming (bra). This drops the rho^4 defect rather than patching it: the local
density matrix is written straight into the tensor, so a mixed initial state is
prepared as requested. bloch now rejects Bloch vectors longer than one instead
of silently returning a non-positive matrix.

src/observables.jl wraps the boundary. A mixed-state reduced density matrix is
contracted against the environment of the *traced* network, since
InfiniteSquareNetwork(pepo) sandwiches the PEPO with its adjoint and yields the
Hilbert-Schmidt norm instead; trace_physicalspaces gives the partition function
whose boundary is wanted. partialtrace then defers to reduced_densitymatrix,
which normalises to unit trace and accepts non-adjacent indices.

test/state_test.jl covers it: three mixed states are prepared to 1e-10 in every
Pauli component where the old constructor deviated by 0.476, padding the bond
space leaves the state unchanged, and the tensor spaces match the documented
convention. 32 tests pass.

PEPSKit's SUWeight(pepo) already builds trivial weights, so an initialweights
helper written here was dropped in favour of it. Trivial weights are a gauge
choice for a product state, whose true spectrum (1, 0, ...) is not invertible.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The term algebra needed no migration at all: abstractterm.jl, terms/*.jl,
hamiltonian.jl, liouvillian.jl and ising.jl are matrix-level code with no
reference to TensorKit or TRGroups, so they are simply loaded again. Only
trotterlayer.jl and equation.jl touch either, and both belong to the lattice
layer that the evolution rebuild replaces; they stay unloaded.

src/propagator.jl builds the two-site object fresh. A Liouvillian acts jointly
on the ket and bra of a site, which PEPSKit's own two-site operators cannot
express, so each site's physical pair is fused into one leg and the propagator
becomes an ordinary (2, 2) map on fuse(P (x) P') - the shape PEPSKit requires.
The vectorised superoperator's legs are reordered per site before fusing, using
the same vectorised_tensormap helper that fixed the earlier leg reversal.

Named propagator rather than gate throughout: this package simulates open
system dynamics, and the package's own vocabulary is stepper and propagator,
with dynamicalmap.jl for the evolution layer. The plan's Phase 4 design was
renamed to match.

test/propagator_test.jl, 9 tests: the propagator's spectrum matches
exp(dt * eigvals(generator)) for both the Lindblad and thermal generators,
which tests the generator and time step without depending on any index
convention; the Lindblad generator's spectrum matches one built independently
from -i[H,rho] + sum(L rho L' - {L'L, rho}/2); the trace functional is left
invariant by the Lindblad propagator and not by the thermal one; dt = 0 gives
the identity; and the spaces are (F (x) F) <- (F (x) F).

Three of my own errors, all caught: the module dropped Accessors and TupleTools
when trimmed, vectorisedidentity built its index in the dual space through an
ambiguous @tensor form and now returns a plain vector, and the spectrum test
needed a rounded sort key since conjugate pairs degenerate at 1e-15 were paired
inconsistently between the two lists.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The propagator no longer fuses each site's physical pair: it is an AbsTen{4,4}
on (ket1, bra1, ket2, bra2), which is the object built before fusing, so the
change removes a step rather than adding one.

Keeping the legs separate costs exactly four mirrored functions, in
src/reducedbond.jl. PEPSKit's bond_tensor_first/bond_tensor_last split off one
physical leg, selected by gate_ax, because their two-site operators act on the
ket and the bra in separate passes; a Liouvillian acts on both at once, so the
reduced tensor is a GenericMPSTensor{S,3} carrying (bond, ket, bra) rather than
(bond, phys). Cluster extraction, weight absorption, bond rotation and the
truncations are reused unchanged.

test/reducedbond_test.jl, 12 tests: splitting is invertible at D = 1, 2, 4 for
both ends of a bond, the reduced tensors carry both physical legs with the
expected spaces, and a propagator built at dt = 0 leaves a random pair unchanged
through the whole split, apply, truncate, recombine path.

Two bugs of mine, both caught. applypropagator permuted the SVD's right factor
to (chi, bond, ket) instead of (bond, ket, bra), so the recombined pair had
mismatched spaces. And the trace-preservation test compared two tensors living
in mutually dual spaces; it now works in dense linear algebra where the index
layout is written out.

MatrixAlgebraKit becomes a direct dependency, since TruncationStrategy appears
in a signature here.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
bondstep! mirrors PEPSKit's _su_iter! with only the reduced bond split and the
propagator application replaced: cluster extraction, weight absorption, bond
rotation, the SVD truncation and the write-back are all theirs. trotterstep!
sweeps every bond of the unit cell along both axes and returns the largest
truncation error.

Verified against dense linear algebra rather than against the old code: a
one-local Lindblad model at D = 1 on a 2x2 cell reproduces exp(t L) to better
than 1e-8 over ten steps of dt = 0.02. That case is exact for two independent
reasons - no interaction means the state stays a product state so nothing is
truncated, and every bond generator commutes with every other so the sweep
carries no Trotter error - which makes it a sharp test of the whole path: state
construction, propagator, split, application, truncation, weight handling,
rotations and the observable contraction. This is benchmark B1 restored on the
new stack. An interacting sweep at D = 4 preserves unit trace, hermiticity and
positivity of the one-site reduced density matrix, which is B3.

Two facts the sweep depends on, both documented at the call site. Each site
belongs to four bonds and each bond carries half of the one-local part of the
propagator, so a sweep applies that part twice per site and the model handed to
bondpropagator must halve it; the old code did the same thing through
changeaxis(op, axis; a = 0.5). And truncspace rejects dual spaces, so the bond
space taken from the weight is flipped when dual, as PEPSKit does in
_get_fixedspacetrunc.

58 tests pass across the four new files.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
trotterise turns a model into a sequence of (propagator, axis) pairs, one per
layer per axis, halving the one-local part of any layer that is not restricted
to a single axis: a full step sweeps both axes and each sweep applies that part
once per site. This is where the old changeaxis(op, axis; a = 0.5) bookkeeping
now lives, so the sweep functions no longer carry it.

The sweep is split accordingly: axissweep! updates every bond along one axis and
timestep! runs a whole decomposition.

Simulation holds the state, its weights and the decomposition; simulate! runs the
steps with the callback schedule, wall-clock budget and convergence break of the
original driver. Convergence is weightdistance, the largest relative change per
unit time in the bond weights, with a changed bond space counting as Inf.

Model and TrotterLayer stay parked: they carry the lattice-broadcasting methods
typed on TRGroups geometries, which belong to the FSM and cluster work. A layer
here is a bare term and a model is a term or a vector of them, optionally paired
with an axis as `term => axis`.

test/simulation_test.jl restores two benchmarks through the public API. B1: a
one-local Lindblad trajectory matches exp(t L) to better than 1e-8 over ten
steps. B6: Larmor precession under H = (h/2) Z from a state along +x gives
<Y> = +sin(omega t) to 1e-8, with the reversed convention excluded by 0.1 rather
than merely disfavoured. Also covered: the callback schedule under maxshots, the
accessors, and the early stop on breaktol.

74 tests pass across five files.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
expval now takes one operator per site. Its ordering had to be pinned rather
than assumed: the dense reshape runs its composite index with the last site
slowest, so the operator for the first site is the last kron factor, which a
product state with distinct states per site exposed immediately.

B2, B4, B7 and B9 move to test/benchmarks_test.jl and run through Simulation.
Each keeps its discriminating comparison, so passing means the right answer
rather than a loose bound: B2 excludes the reversed Hamiltonian sign by more
than 1 while agreeing to 1e-3, B4's measured slopes fix beta = t rather than
2t, and B9 excludes the exchanged site order by more than 1.

B7 now agrees at every D from 1 to 8, where the old code deviated by 2.02 at
D = 4 and 72.5 at D = 8. That defect is resolved by construction: PEPSKit guards
the weight inversion with sdiag_pow(s, pow; tol), so nothing here had to fix it.
Likewise the rho^4 constructor defect is gone because productpepo writes the
density matrix straight into the tensor.

The plan now records where each benchmark lives and that all three defects found
in task 0.1 are resolved, only the first by an actual fix.

85 tests pass across six files.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Adds correlationlength, correlator, negativity and PurifiedDensityMatrix.

A purified state rho = X'X uses the boundary of the two-layer network that
InfiniteSquareNetwork builds for a PEPO, and reads observables through
reduced_densitymatrix with the PEPO as both ket and bra; DensityMatrix instead
uses the traced network, since a Lindblad state is already a density matrix.
Verified by putting sqrt(rho) on every site, which must return rho.

negativity is (norm(rho^T_B, 1) - 1)/2, taken densely: the partial transpose is
a permutation of the reduced density matrix's four indices. Checked against a
product state (0), a Bell state (1/2), and Werner states either side of the
p = 1/3 separability threshold.

correlationlength returns the largest boundary correlation length over rows and
columns; PEPSKit's correlation_length gives the per-row and per-column values.

Three of my expectations were wrong and are now stated correctly. A product
state's correlation length is not zero but ~1/46: the second transfer eigenvalue
is ~1e-20 rather than 0, and -1/log|lambda2/lambda1| turns that into a finite
noise floor, which also fluctuates with the random CTMRG initialisation. And the
connected ZZ correlator of the tutorial model stays at 1e-13 while the
transverse channel reaches 1.5e-3, because (J/4) sum ZZ commutes with every Z_i,
so the Z sector closes and no ZZ correlation is generated from a Z-diagonal
initial state. The test now asserts both halves, which pins a symmetry of the
model rather than merely the presence of correlations.

pepo_test.jl, densitymatrix_test.jl and trotterise_test.jl move to test/legacy
with the _test suffix dropped so ReTestItems no longer discovers them, alongside
the old-API benchmark script and the trajectory recorder, with a README saying
what replaced each. Nothing is deleted, and `Pkg.test` is meaningful again:
runtests(TimeEvolutionPEPO) discovers 10 test items and 107 tests pass.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
PowOp and GaussOp were loaded but non-functional: both call
approx_as_sum_of_exp, which lives in longrangeheisenberg.jl, and that file was
not included. Constructing any power-law term raised UndefVarError. fsm.jl comes
with it, since it is what turns a long-range term into an operator PEPO and its
state count is that operator's bond dimension.

Neither file needed porting to load: longrangeheisenberg.jl touches TensorKit
only inside its two fsm* builders, and fsm.jl only inside fsmpepo, none of which
run at load time. Both are placed after ising.jl, since fsm.jl names Hamiltonian
and Liouvillian in method signatures.

107 tests still pass.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
gaugefix! wraps PEPSKit's SUGauge and BPGauge, writing the result back into the
state and weights a Simulation holds, since gauge_fix returns new objects.
bpenvironment computes the belief propagation fixed point BPGauge needs, and
bondenvironment the CTMRG boundary of the PEPO sandwiched with its adjoint.

test/gauge_test.jl, 7 tests: both algorithms leave one-site observables
unchanged, which is what a gauge transformation must do while the weights move;
SUGauge is idempotent, reaching the same weights when applied twice; and a
product state is already in the fixed-point gauge. Note BPGauge fuses the two
physical legs internally and unfuses on the way out - PEPSKit's implementation
detail, not a change to how the state is stored here.

The environment-based truncation does not work and is parked as
test/envtruncation_wip.jl, undiscovered by ReTestItems. bondenv_ctm contracts a
physical index between the two layers of the orthogonal remainder, which exists
only because PEPSKit's split leaves one physical leg there; ours leaves none.
A standalone probe appeared to accept a four-leg remainder, but that was the
@tensor macro mis-binding indices in a case where the dimensions coincided, so
the environment must be mirrored for a four-leg remainder instead.

applypropagator gains an untruncated method, and bondstep!, axissweep! and
timestep! take truncalg and env keywords. That path is inert until the
environment is mirrored; the default SVD truncation is unchanged.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Two blockers from the previous attempt are resolved. bondenv_ctm does accept a
remainder with no physical leg: _prepare_site_tensor(t::PEPSOrth) pads it with a
trivial unit leg, so no mirror is needed, and my earlier claim that it had to be
rewritten was wrong. Vertical bonds additionally need the environment rotated as
_bond_rotation rotated the tensors, with the bond position carried through that
rotation, which _rotatedbond now does.

What remains is not a convention problem but a design one: the environment is
computed once per time step while each bondstep! rewrites two site tensors and a
bond space, so by the second axis the state and the environment disagree and the
contraction fails outright on a dual mismatch. A full-environment update must
recompute the environment as the sweep proceeds, which is precisely the cost
that separates it from simple update.

The test stays parked as test/envtruncation_wip.jl with that diagnosis recorded.
The truncalg and env keywords remain inert, and the default SVD path is
unchanged.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
`timestep!` takes `truncalg` (`FullEnvTruncation` or `ALSTruncation`) and `envalg`
(a `BondEnvironment`), and recomputes the CTMRG boundary before each axis sweep,
warm-started from the previous one when the state's virtual spaces still match it.
`simulate!` passes the same three keywords through.

`bondstep!` now restores the arrow the bond carried before the update. The SVD in
`applypropagator` picks the arrow itself, and `_bond_rotation` maps the rotated
frame's east bond onto a north-south bond carrying the opposite one, so the first
update of every vertical bond used to flip the north and south legs it wrote back.
Nothing in simple update reads those arrows, but `bondenv_ctm` contracts the state
against an environment built earlier in the sweep, and after a flip that
contraction failed with `SpaceMismatch: ((C^1)' x (C^1)') != ((C^1)' x C^1)`.

test/envtruncation_test.jl runs again: the D = 1 one-local model reproduces
exp(t L) to 1e-8 through the environment path, full-environment truncation is at
least as faithful as a bare SVD at the same bond dimension, and ALS truncation
preserves trace, hermiticity and positivity.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Replaces the arrow-restoring flip of d0992cc. `axissweep!` now takes `envalg` and
solves for a CTMRG boundary immediately before each `bondstep!`, warm-started from
the boundary of the previous bond when the state's virtual spaces still match it,
so the environment always describes the state the bond is about to be updated in.
The state's arrows are left as the SVD and `_bond_rotation` produce them.

Cost against one boundary per axis: test/envtruncation_test.jl goes from 4.3 s to
6.8 s and from 9.9 GB to 18.6 GB of allocation on a 2x2 unit cell, where per-bond
means four boundary solves per axis sweep instead of one.

Known consequence, unchanged from before d0992cc: a vertical bond update flips the
arrow on the bond it writes back, so every bond's weight changes space on its first
update and `weightdistance` reports `Inf` for the first time step of a simulation.
Later steps report a finite convergence measure.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
`axissweep!` passes a vertical bond as [(r+1, c), (r, c)], the orientation
`PEPSKit._trivial_gates` uses, and swaps the propagator's two sites to match.
The other order makes `_nn_bondrev` report `rev = true`, which sends the pair
through `rot180`/`rotl90`; those relabel the four virtual legs without touching
their arrows, so the north-south legs land in the slots where the SVD in
`applypropagator` writes the opposite arrow, and every bond came back with its
arrow flipped on first update.

With the arrows now invariant across a sweep the boundary is solved for once per
axis again, warm-started from the previous axis, replacing the per-bond solve of
b648d74: test/envtruncation_test.jl drops from 6.8 s to 4.4 s and from 18.6 GB
to 10.1 GB of allocation. `weightdistance` also reports a finite convergence
measure for the first time step now (9.99 against `Inf` before), since bond
weights no longer change space.

B9 covers the propagator swap: an asymmetric two-site term evolves the wrong way
round if the sites are reordered without it.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Every page described the pre-refactor architecture: `PEPO`, `Simulation(pepo,
model, method)`, TRGroups unit cells and boundary algorithms, and extension
points (`updatetensors` A/B/C, `environment!`, `newtruncators`) that no longer
exist. Rewritten against the current API — `productpepo`, `InfinitePEPO`,
`SUWeight`, `Simulation(state, model; timestep)`, `DensityMatrix(state, bonddim)`
— with `interfaces.md` now documenting the `trunc`/`truncalg`/`envalg` keywords
and PEPSKit's `bond_truncate` as the way to add a truncation scheme.

man/tensormap.md had the virtual index order as east-south-west-north; PEPSKit's
PEPOTensor is north-east-south-west. Its `permutebasis` doctest is replaced by an
explicit row/column permutation, since that function was not carried over.

The `@docs` blocks listed 17 names with no binding left (`TEBD`, `Model`,
`TrotterLayer`, `NTU`, `FET`, `EVB`, `DZM`, `AbstractPEPO`, `expandspace`, …).
Pruned to what exists, and docstrings added for `quantumstate`, `bondweights`,
`isexpandable`, `isaxial`, `oneterms` and `dissipator`, which the library pages
reference but which had none. A dangling `[`TEBD`](@ref)` in the `Ising`
docstring is now prose.

The docs environment could not have built on CI: the workflow installed
TensorRenormalizationGroups from git, which the package no longer depends on.
Removed, and docs/Project.toml + Manifest.toml resolved against registered
PEPSKit v0.8.1 with no absolute paths. `julia --project=docs docs/make.jl`
completes with no errors.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Site-anchored terms, per-endpoint one-local placement, push!/distribute!
separated, expandable terms deferred to the long-range work.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
jack-dunham and others added 30 commits August 18, 2026 12:48
Tasks 5 and 6 rewrite src/model/{trotterlayer,model}.jl in place instead of
git mv-ing the parked versions to src/dev/, and Task 7 adds a bondpropagator
method beside its sibling in src/propagator.jl instead of extracting a helper
into another file. Relocating code makes the diff unreadable; rearranging files
is a separate change.
Ising was exported and documented but generator threw for it, so the type could
not be evolved. It is expanded to the equivalent Hamiltonian or Liouvillian; the
low-rank exponentiation the parked equation.jl provided is not restored here.
trotterise([LocalOp(X,X), LocalOp(Y,Y)], ...) threw on scalelocal, so the
multi-layer form only accepted terms already wrapped in a Hamiltonian. Each
layer is now promoted, dissipative terms to a Liouvillian. A layer holding a
long-range term now names the type it cannot exponentiate.
productpepo(rho::AbstractMatrix, bond) selected the unit-cell method, so a single
density matrix was broadcast over its scalar entries. A matrix of numbers is now
one local state; a matrix of matrices, tuples or ThermalStates is a unit cell.
sparsetensormap.jl is superseded — TensorKit 0.17 covers what it provided, and
its only remaining callers (evolution/abstractevolution.jl, truncation/
truncation.jl) are themselves unloaded and slated for deletion. Its commented-out
include goes with it.

model/terms/toricop.jl is unfinished work rather than superseded work: the
PlaquetteVertex term was never wired to a generator or a propagator, so it moves
to src/dev to be finished later.

Also drops the word "parked" from the documents in favour of naming which
operation happened - deleted, moved to src/dev, or left in the tree unloaded -
since the three differ in whether the code can be recovered or run.
space(bondpropagator(...)) == space(bondpropagator(...)) compares
P ⊗ P ⊗ P' ⊗ P', which is fixed by P alone and independent of the model, so
the check could not fail regardless of correctness. Compares the propagators
directly instead, which also exercises the exponentiation.
promote_type on two Hamiltonians already yields Hamiltonian, so no rule is needed
there; Hamiltonian with Liouvillian yields AbstractSumOfTerms, so the Liouvillian
promote_rule is required rather than optional.
data[r, c] is anchored at site (r, c): its two-local part is the coupling to the
neighbour in the +axis direction, its one-local part belongs to that site alone.
Indexing wraps periodically so a 1x1 layer is uniform. Replaces the superseded
TRGroups-based version in place; git history holds the old one.
isaxial(::Dissipator) was unconditional false, so composing it with
isaxial(::Liouvillian) made Liouvillian(H) with no jump operators
non-axial too, spuriously blocking a plain Hamiltonian from pinning
to an axis. isaxial(D) = isempty(terms(D)) fixes the empty case only.
push! always appends a layer; distribute! spreads a term over the existing
layers with a 1/nlayers prefactor and throws on an empty model. The superseded
version triggered distribution by argument type, which made the behaviour
depend on whether the term was wrapped. Replaces the superseded version in place.

Deviates from the brief's draft implementation in two ways, both required to
make it work: `layers` is `Vector{TrotterLayer}`, not `Vector{TrotterLayer{T}}`,
because TrotterLayer's own constructor always narrows to the concrete shape of
its terms (e.g. Hamiltonian{Tuple{LocalOp{2}}}), which is never a subtype of
TrotterLayer{Hamiltonian} since struct type parameters are invariant; storing
into a Vector{TrotterLayer{Hamiltonian}} then fails with a TypeError on every
non-empty push. `modelkind` maps a term's concrete type to its family
(Hamiltonian, Liouvillian or Ising) for the push!/distribute! checks, since
`eltype(layer.data) <: T` compared two different concrete shapes of the same
family and rejected valid pushes.
MODEL_SPEC.md and MODEL_PLAN.md were working documents for this piece of work.
They remain in git history and in the untracked execution workspace.
…tating

distribute! added the term to every layer via t + share without checking it
against the model's type parameter, so distribute!(model, Dissipator(Z)) on a
Model{Hamiltonian} silently turned every layer into a Liouvillian. Extracts
push!'s family check into _checkfamily, reused by both, and calls it in
distribute! before share is computed or any layer is touched. Also trims the
modelkind comment to two lines.
Design and planning documents are working artifacts; carrying them in the repo
makes every diff that touches them noisy. They remain in git history, and local
copies are kept in the gitignored dev/design/.
Model layer, tasks 1-6 of 10:
- generator places one-local terms per bond endpoint, so a bond propagator can
  represent a field that differs between its two sites
- Ising reaches a generator; bare terms promote to Trotter layers; a bare density
  matrix is one site in productpepo
- TrotterLayer carries a term per unit-cell site; Model is an ordered list of them
  with push!/distribute! separated

Also: environment-based bond truncation with one CTMRG boundary per axis sweep,
vertical bonds swept south to north so bond arrows stay invariant, docs rewritten
for the PEPSKit API, SparseTensorMap deleted and the toric-code term moved to
src/dev.

trotterise does not yet accept a Model and axissweep! still applies one propagator
to every bond, so Model and TrotterLayer are exported but not reachable from
Simulation. Tasks 7-10 close that.
generator is defined in src/propagator.jl, but the Ising methods were added to
src/model/ising.jl, splitting the function across two files. They move back
beside the other generator methods. _equation stays in ising.jl, where it reads
Ising's fields alongside scalelocal, which reconstructs the coupling the same way.
bondenvironment, BondEnvironment and _matchesnetwork were introduced in
src/gauge.jl and relocated to src/reducedbond.jl when the environment truncation
was wired up, which made that commit's diff read as a deletion plus an unrelated
addition. They go back; truncatebond and the reduced bond split stay where they
were introduced.
Return relocated functionality to its original file: generator's Ising methods to
src/propagator.jl, and bondenvironment to src/gauge.jl.
The pretty-printer was dropped in 55cd7da when the simulation driver was rewritten,
so a Simulation printed with Julia's default struct display. Restored with the
METHOD section replaced by the model and the Trotter layer count, since there is
no method field any more.
Seven files added during the refactor are folded back into the pre-refactor files
whose functionality they replaced, so the diff for each subsystem reads as a change
to one file rather than a deletion plus an unrelated addition:

  state.jl              -> pepo/pepo.jl
  observables.jl        -> densitymatrix.jl
  propagator.jl         -> model/equation.jl
  reducedbond.jl        -> truncation/isometries.jl
  simpleupdate.jl       -> truncation/simpleupdate.jl
  trotterdecomposition.jl -> evolution/dynamicalmap.jl
  gauge.jl              -> gauge/abstractgauge.jl

generator/bondpropagator and gaugefix! had no pre-refactor equivalent; they go to
the nearest subsystem, model/equation.jl and gauge/abstractgauge.jl. The include
order is unchanged. These names now describe the architecture they replaced rather
than their contents, which the module header and CLAUDE.md record.
Restore the Simulation show method, and fold the seven source files added during
the refactor into the pre-refactor files whose functionality they replaced.
Each layer becomes a Matrix of propagators indexed as its term matrix is; the
bond anchored at (r, c) takes its coupling and first-slot one-local part from
data[r, c] and its second-slot part from the neighbour in the +axis direction.
A uniform layer yields a 1x1 matrix.

Also fixes a stale assertion in propagator_test.jl that checked trotterise's
pre-Task-7 return shape.
axissweep! takes a Matrix of propagators and indexes it with mod1, so a 1x1
matrix serves any unit cell and a full-size one varies per bond. Simulation
rejects a layer whose size is neither 1x1 nor the state's unit cell rather
than tiling it silently.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
At D = 1 with no coupling each site evolves independently, so both sublattices
are checked against dense linear algebra. This pins the one-local prefactor and
which kron factor of the propagator is the first site, neither of which a
uniform model can distinguish.

Brief used nsteps = 5 (t = 0.1), giving norm(want_active - want_idle) = 0.076,
which fails its own > 0.1 discrimination bound. Raised nsteps to 15 (t = 0.3,
separation 0.208); the simulation still matches exp(tL) to ~1e-14.
Documents the Model/TrotterLayer/push!/distribute! API in the manual and
library reference, and corrects the pipeline description in interfaces.md
now that axissweep! takes a propagator matrix. Also drops a stale @ref to
push! (undocumented) that broke the cross-reference build.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
- Simulation docstring now names Model as the general form of `model`.
- timestep!/axissweep! docstrings and a loop variable consistently say
  `propagators` where the value is a per-layer matrix, not a bare tensor.
- bondpropagator's 3-argument method delegates to the 4-argument one
  instead of duplicating its body, mirroring generator's own delegation.
- distribute! now builds every replacement layer before writing any of
  them back, so a later layer's axial-check throw leaves the model with
  its original layers rather than a partially replaced set. Added a test
  reproducing the partial mutation and covering the fix.
- Removed Model's incorrect eltype method; iterate already yields
  TrotterLayer and nothing in the tree calls eltype(::Model).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
On-site dissipation has no direction to be inconsistent with an axis, so a
dissipative Liouvillian layer can now be pinned to one, matching the existing
behavior of a dissipative Ising layer. Updates the two model_test.jl tests that
asserted the old refusal, and the "later layer rejects the axial merge" test
that used Dissipator as its non-axial example (now uses GaussOp2d over a
Model{Hamiltonian}, since GaussOp2d does not compose with Model{Liouvillian}).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Model layer, tasks 7-10 of 10 plus the pre-merge fix wave:
- trotterise returns a Matrix of propagators per Trotter layer, one per unit-cell
  site; axissweep! indexes it per bond; Simulation rejects a layer whose matrix is
  neither 1x1 nor the state's unit cell
- staggered-field acceptance test against dense exp(tL): at D = 1 with no coupling
  each site evolves independently, so both sublattices are checked directly. This
  is the only test that can detect the one-local prefactor or which propagator slot
  is the first site
- documentation for Model, TrotterLayer and site-dependent terms
- a Dissipator is now axial, so a dissipative layer can be pinned to one axis; a
  dissipative Liouvillian and an equivalent Ising now agree
- distribute! builds every replacement layer before assigning any, so a rejected
  term leaves the model untouched

Full suite 170/170 across 16 test items.
…c/dev

Every file under src/ outside src/dev/ is now loaded: 24 files, 24 includes.
The block of commented-out includes is gone.

Deleted, superseded by PEPSKit or by the new code (21 files, ~4200 lines):
  abstractlattice.jl, edgelattice.jl      InfinitePEPO indexes its own unit cell
  contract/{ctmrg,vumps}.jl               PEPSKit's leading_boundary and CTMRGEnv
  pepo/*.jl (7 files)                     InfinitePEPO, SUWeight, InfiniteSquareNetwork
  evolution/{abstractevolution,timeevolution}.jl   the TEBD method zoo
  updatestyle.jl, model/abstractaxis.jl   the UpdateStyle trait zoo
  model/{hamiltonianterm,lindbladterm,longrangeterm}.jl   superseded by model/terms/
  truncation/truncation.jl                the new sweep and PEPSKit's truncations
  utils.jl, tensorkitutils.jl             both encode the pre-PEPSKit leg order --
    _DIR_TO_DOMIND is east-south-west-north and westbond/southbond read the
    codomain, so their accessors name the wrong legs for a PEPSKit PEPOTensor

Moved to src/dev/, unfinished work with no replacement (9 files, ~3400 lines):
  truncation/clu/{clu,metric}.jl   the NTU cluster metric
  truncation/iterative.jl          IterativeSU
  gauge/wtg.jl                     WTG
  model/{mpo,pepso,heisenberg,axialheisenberg}.jl   the long-range PEPO builders
    the deferred finite-state-machine work needs
  model.jl                         model recipes written against the removed ham/lin API

No live file depends on any deleted definition. Fifteen names appeared to be
shared, but each is either also defined in a loaded file or occurs only in a
comment, a docstring, or a parameter name: verified by checking that adjacent,
bonds, ham, metric, permutedom, unpack, tensors and AbsTen resolve nowhere in the
loaded module, then reading every occurrence.

Full suite 170/170 across 16 test items.
Delete the 21 superseded source files, move the 9 unfinished ones to src/dev/.
Every file under src/ outside src/dev/ is now loaded and the commented-out
include block is gone. Full suite 170/170.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant