Code review fixes: correctness, robustness, packaging, ITensorMPS 0.4

Project.toml

@@ -3,19 +3,15 @@ uuid = "9fc92dc4-a6d3-4a85-b5fe-6a206e4c3402"
 version = "0.1.0"
 
 [deps]
-EDKit = "f6058819-fd6f-45d3-a3c8-4753d8c54d1a"
 ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [compat]
-EDKit = "0.4.7, 0.5"
-ITensorMPS = "0.3.44"
+ITensorMPS = "0.3.44, 0.4"
 ITensors = "0.9.25"
 KrylovKit = "0.10.2"
-Plots = "1.41.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/make.jl

@@ -1,6 +1,10 @@
 using Documenter
 using MPSToolkit
 
+# Make the package available inside `jldoctest` blocks so docstring examples are executed
+# and verified during the docs build (`makedocs` runs doctests by default).
+DocMeta.setdocmeta!(MPSToolkit, :DocTestSetup, :(using MPSToolkit); recursive=true)
+
 const PRETTY_URLS = get(ENV, "CI", "false") == "true"
 
 makedocs(

src/MPSToolkit.jl

@@ -24,8 +24,8 @@ include("models/spinhalf.jl")
 include("operator_space/helpers.jl")
 include("operator_space/dmt.jl")
 include("operator_space/daoe.jl")
-include("chebyshev/reconstruction.jl")
 include("chebyshev/types.jl")
+include("chebyshev/reconstruction.jl")
 include("chebyshev/moments.jl")
 
 """
@@ -137,6 +137,8 @@ Namespace for Chebyshev moments, kernels, reconstruction, and energy-window cuto
 module Chebyshev
 using ..MPSToolkit:
   ChebyshevRescaling,
+  chebyshev_rescaling,
+  rescale_hamiltonian,
   SpectralFunction,
   chebyshev_moments,
   energy_cutoff!,
@@ -145,6 +147,8 @@ using ..MPSToolkit:
   reconstruct_chebyshev,
   spectral_function
 export ChebyshevRescaling,
+  chebyshev_rescaling,
+  rescale_hamiltonian,
   SpectralFunction,
   chebyshev_moments,
   energy_cutoff!,
@@ -162,6 +166,6 @@ export tebd_evolve!, dmt_evolve!, tdvp_evolve!, local_gates_from_hamiltonians, t
 export pauli_matrices, pauli_basis, pauli_components
 export spinhalf_matrices, spinhalf_xyz_bond_hamiltonian, spinhalf_tfim_bond_hamiltonian
 export pauli_siteinds, pauli_basis_state, pauli_total_sz_state, pauli_gate, pauli_gate_from_hamiltonian, pauli_lindblad_generator, pauli_gate_from_lindbladian, DMTOptions, dmt_step!, dmt_evolve!, pauli_daoe_projector, fdaoe_projector
-export ChebyshevRescaling, SpectralFunction, chebyshev_moments, energy_cutoff!, jackson_damping, jackson_kernel, reconstruct_chebyshev, spectral_function
+export ChebyshevRescaling, chebyshev_rescaling, rescale_hamiltonian, SpectralFunction, chebyshev_moments, energy_cutoff!, jackson_damping, jackson_kernel, reconstruct_chebyshev, spectral_function
 
 end

src/bases/pauli.jl

@@ -11,11 +11,14 @@ Return the dense spin-1/2 Pauli matrices.
   included, and `(X, Y, Z)` otherwise.
 
 # Examples
-```julia
-julia> pauli_matrices().Z
-2x2 Matrix{ComplexF64}:
- 1+0im   0+0im
- 0+0im  -1+0im
+```jldoctest
+julia> Z = pauli_matrices().Z;
+
+julia> Z == [1 0; 0 -1]
+true
+
+julia> keys(pauli_matrices(; include_identity=false))
+(:X, :Y, :Z)
 ```
 """
 function pauli_matrices(; include_identity::Bool=true)
@@ -60,6 +63,13 @@ Decompose a `2 x 2` operator into Pauli-basis coefficients.
 # Returns
 - A named tuple of Hilbert-Schmidt coefficients, normalized so that
   `operator == sum(coeffs[name] * pauli_matrices()[name] for name in keys(coeffs))`.
+
+# Notes
+- These coefficients are taken against the **unnormalized** basis `{I, X, Y, Z}`. They are
+  **not** the operator-space state amplitudes produced by [`pauli_basis_state`](@ref) and the
+  other operator-space helpers, which use the normalized basis `σ / √2`. The two conventions
+  differ by a factor of `(√2)^L` over `L` sites, so do not feed `pauli_components` output
+  directly into the operator-space state constructors.
 """
 function pauli_components(operator::AbstractMatrix; include_identity::Bool=true)
   size(operator) == (2, 2) || throw(ArgumentError("Pauli decomposition requires a 2 x 2 operator"))

src/chebyshev/moments.jl

@@ -195,6 +195,15 @@ the Hamiltonian has already been rescaled into the target Chebyshev window.
 
 # Returns
 - The mutated `psi`.
+
+# Notes
+- This filters each site against its **local** effective problem (the single-site effective
+  Hamiltonian from the MPO environment), not the global Hamiltonian spectrum. Out-of-window
+  weight that is delocalized across entangled bonds is therefore only weakly removed, so a
+  strongly entangling state may retain such weight even after many sweeps.
+- `tol` is compared against a heuristic per-sweep convergence estimate (the accumulated
+  relative removed weight), not a calibrated RMS energy error; treat it as a stopping monitor
+  rather than an error bound.
 """
 function energy_cutoff!(
   psi::MPS,
@@ -227,7 +236,8 @@ end
 """
     _energy_cutoff_sweep!(projector, psi; krylovdim, window)
 
-Run one bidirectional energy-cutoff sweep and return its RMS error estimate.
+Run one bidirectional energy-cutoff sweep and return a heuristic convergence estimate (the
+accumulated relative removed weight across the site visits, not a calibrated RMS error).
 """
 function _energy_cutoff_sweep!(projector::ProjMPO, psi::MPS; krylovdim::Integer, window::Real)
   nsites = length(psi)
@@ -286,6 +296,10 @@ Project a local Krylov basis state onto the eigenmodes whose energies lie within
 - `(projected_coefficients, removed_weight)`.
 """
 function _project_energy_window(eigenvalues::AbstractVector, eigenvectors::AbstractMatrix; window::Real=1.0)
+  # The input local state is e_1 in this Krylov basis because `_lanczos_tridiagonal` seeds
+  # `basis[1] = normalize(tensor)`. That invariant is load-bearing here: `coefficients[1] = 1`
+  # represents the (normalized) input state, and `vector[1]` is each eigenmode's overlap with
+  # it. If the Lanczos seed ordering ever changes, this projection must be revisited.
   coefficients = zeros(eltype(eigenvectors), size(eigenvectors, 1))
   coefficients[1] = one(eltype(coefficients))
   removed_weight = 0.0

src/chebyshev/types.jl

@@ -7,6 +7,12 @@ coordinate `x = (ω - center) / halfwidth`.
 # Fields
 - `center`: Physical frequency mapped to `x = 0`.
 - `halfwidth`: Positive scale factor mapping the target energy window to `[-1, 1]`.
+
+# Notes
+- Choose `halfwidth` with a small safety margin (e.g. `(E_max - E_min) / (2 - ε)`) so that
+  every eigenvalue maps strictly inside `(-1, 1)`. Spectral weight that maps to exactly `±1`
+  lands on the band edge, where the Chebyshev measure `1 / (π √(1 - x²))` diverges and the
+  reconstructed [`SpectralFunction`](@ref) returns `0.0`, silently dropping that weight.
 """
 struct ChebyshevRescaling
   center::Float64
@@ -20,6 +26,57 @@ struct ChebyshevRescaling
   end
 end
 
+"""
+    chebyshev_rescaling(emin, emax; margin=0.05)
+
+Build the [`ChebyshevRescaling`](@ref) that maps the spectral window `[emin, emax]` into the
+Chebyshev coordinate so the extreme eigenvalues land at `±(1 - margin)`, leaving a safety
+buffer so no spectral weight reaches the band edge `±1`.
+
+# Arguments
+- `emin`, `emax`: Lower and upper bounds (or estimates) of the Hamiltonian spectrum.
+
+# Keyword Arguments
+- `margin`: Fractional safety buffer in `[0, 1)`. The extreme eigenvalues map to `±(1 - margin)`.
+
+# Returns
+- A [`ChebyshevRescaling`](@ref) with `center = (emax + emin) / 2` and
+  `halfwidth = (emax - emin) / (2 * (1 - margin))`.
+"""
+function chebyshev_rescaling(emin::Real, emax::Real; margin::Real=0.05)
+  emax > emin || throw(ArgumentError("chebyshev_rescaling requires emax > emin"))
+  0 <= margin < 1 || throw(ArgumentError("chebyshev_rescaling requires 0 <= margin < 1"))
+  center = (emax + emin) / 2
+  halfwidth = (emax - emin) / (2 * (1 - margin))
+  return ChebyshevRescaling(center, halfwidth)
+end
+
+"""
+    rescale_hamiltonian(H::MPO, emin, emax; margin=0.05, cutoff=1e-14)
+
+Rescale an MPO Hamiltonian into the Chebyshev band.
+
+# Arguments
+- `H`: Hamiltonian `MPO`.
+- `emin`, `emax`: Spectral bounds (or estimates) of `H`.
+
+# Keyword Arguments
+- `margin`: Safety buffer forwarded to [`chebyshev_rescaling`](@ref).
+- `cutoff`: Truncation cutoff used when subtracting the identity shift.
+
+# Returns
+- `(H_rescaled, rescaling)` where `H_rescaled = (H - center) / halfwidth` has its spectrum
+  strictly inside `(-1, 1)`, and `rescaling::`[`ChebyshevRescaling`](@ref) inverts the map for
+  [`spectral_function`](@ref). Feed `H_rescaled` to [`chebyshev_moments`](@ref) and reuse
+  `rescaling.center` / `rescaling.halfwidth` when reconstructing the physical spectral function.
+"""
+function rescale_hamiltonian(H::MPO, emin::Real, emax::Real; margin::Real=0.05, cutoff::Real=1e-14)
+  rescaling = chebyshev_rescaling(emin, emax; margin=margin)
+  identity_mpo = MPO(firstsiteinds(H), "Id")
+  shifted = add(H, -rescaling.center * identity_mpo; cutoff=cutoff)
+  return shifted / rescaling.halfwidth, rescaling
+end
+
 """
     SpectralFunction(moments, kernel, rescaling)
 
@@ -65,6 +122,11 @@ end
     (spectrum::SpectralFunction)(ω)
 
 Evaluate a reconstructed spectral function at one physical frequency `ω`.
+
+# Notes
+- Frequencies that map outside the open interval `(-1, 1)` (i.e. `|x| ≥ 1`, including the band
+  edges `x = ±1`) return `0.0`. Keep the spectrum strictly inside the band through the
+  rescaling (see [`ChebyshevRescaling`](@ref)).
 """
 function (spectrum::SpectralFunction)(ω::Real)
   x = (ω - spectrum.rescaling.center) / spectrum.rescaling.halfwidth

src/evolution/tdvp.jl

@@ -11,7 +11,9 @@ Run one TDVP evolution call on a finite OBC `MPS`.
 - The mutated `psi`.
 
 # Notes
-- The effective TDVP step count is `evo.nsteps` if present, otherwise `evo.nsweeps`.
+- The effective TDVP step count is `evo.nsteps` if present, otherwise `evo.nsweeps`. When both
+  are `nothing`, the step count is left to `tdvp`, which derives it from `evo.t` and
+  `evo.time_step`.
 - The underlying `tdvp` call returns a new state, so this wrapper copies the result back
   into the original storage to preserve the package-wide in-place API convention.
 """

src/evolution/types.jl

@@ -193,7 +193,10 @@ Construct a [`TDVPEvolution`](@ref) for finite OBC `MPS` states.
 - `updater_backend`: Backend string forwarded to `tdvp`.
 - `updater`: Optional custom updater callback.
 - `normalize`: Whether to normalize the state after evolution.
-- `solver_kwargs`: Additional keyword arguments forwarded to `tdvp`.
+- `solver_kwargs`: Additional keyword arguments forwarded to `tdvp` (e.g. `maxdim`, `cutoff`).
+  Must not contain a key that already has a dedicated argument above (`time_step`, `nsteps`,
+  `reverse_step`, `updater_backend`, `updater`, `normalize`); passing one throws, because it
+  would otherwise be forwarded last and silently override the dedicated argument.
 - `schedule`: Optional metadata used by higher-level workflows.
 
 # Returns
@@ -214,6 +217,10 @@ function TDVPEvolution(
 )
   isnothing(nsteps) || nsteps >= 1 || throw(ArgumentError("TDVPEvolution requires nsteps >= 1 when provided"))
   isnothing(nsweeps) || nsweeps >= 1 || throw(ArgumentError("TDVPEvolution requires nsweeps >= 1 when provided"))
+  reserved_solver_keys = (:time_step, :nsteps, :reverse_step, :updater_backend, :updater, :normalize)
+  conflicting_keys = intersect(keys(solver_kwargs), reserved_solver_keys)
+  isempty(conflicting_keys) ||
+    throw(ArgumentError("TDVPEvolution solver_kwargs must not contain reserved keys $(collect(conflicting_keys)); these are forwarded to `tdvp` last and would silently override the dedicated TDVPEvolution keyword arguments. Set them through the dedicated keyword arguments instead."))
   return TDVPEvolution(
     generator,
     t,

src/observables/entanglement.jl

@@ -17,28 +17,28 @@ function _diagonal_values(tensor::ITensor)
 end
 
 """
-    _entropy_from_values(values)
-
-Compute the von Neumann entropy from Schmidt values or equivalent amplitudes.
+    _entropy_from_probabilities(probabilities)
 
-# Arguments
-- `values`: Schmidt coefficients, amplitudes, or any vector whose squared magnitudes should
-  be interpreted as probabilities.
-
-# Returns
-- The Shannon/von-Neumann entropy `-∑ p log(p)` after normalization.
+Compute the von Neumann entropy `-∑ p log(p)` from a vector of (sub)normalized probabilities.
 
 # Notes
-- Zero-probability entries are skipped safely.
+- The probabilities are renormalized to sum to one; zero entries are skipped safely.
 """
-function _entropy_from_values(values)
-  probabilities = abs2.(values)
+function _entropy_from_probabilities(probabilities)
   total = sum(probabilities)
   iszero(total) && return 0.0
   normalized = probabilities ./ total
   return -sum(p -> iszero(p) ? 0.0 : p * log(p), normalized)
 end
 
+"""
+    _entropy_from_values(values)
+
+Compute the von Neumann entropy from Schmidt values or amplitudes, interpreting the squared
+magnitudes as probabilities.
+"""
+_entropy_from_values(values) = _entropy_from_probabilities(abs2.(values))
+
 function _validate_entanglement_bond(psi::MPS, bond_index::Int)
   1 <= bond_index < length(psi) || throw(ArgumentError("entanglement bond must lie in 1:length(psi)-1"))
   return bond_index
@@ -57,13 +57,11 @@ Return the bond entropy for a finite `MPS`.
 - The von Neumann entropy associated with the Schmidt values across the selected cut.
 
 # Notes
-- This method delegates spectrum extraction to [`entanglement_spectrum`](@ref) and then
-  converts the normalized Schmidt probabilities back into amplitudes before evaluating the
-  entropy.
+- This method delegates spectrum extraction to [`entanglement_spectrum`](@ref), which already
+  returns normalized Schmidt probabilities, and evaluates the entropy directly from them.
 """
 function bond_entropy(psi::MPS, bond::Union{Nothing,Int})
-  spectrum = entanglement_spectrum(psi, bond)
-  return _entropy_from_values(sqrt.(spectrum))
+  return _entropy_from_probabilities(entanglement_spectrum(psi, bond))
 end
 
 """
@@ -80,13 +78,13 @@ Return the normalized Schmidt probabilities for a finite `MPS`.
 
 # Notes
 - Out-of-range bonds throw an `ArgumentError`.
-- The state is copied before orthogonalization, so the input `psi` is not mutated.
+- The non-mutating `orthogonalize` returns a new state, so the input `psi` is not mutated.
 """
 function entanglement_spectrum(psi::MPS, bond::Union{Nothing,Int})
   bond_index = isnothing(bond) ? max(1, length(psi) ÷ 2) : bond
   _validate_entanglement_bond(psi, bond_index)
 
-  centered = orthogonalize(copy(psi), bond_index)
+  centered = orthogonalize(psi, bond_index)
   left_inds = uniqueinds(centered[bond_index], centered[bond_index + 1])
   _, singular_values, _ = svd(centered[bond_index], left_inds)
   values = abs2.(_diagonal_values(singular_values))