refactor: transfer_spectrum and correlation_length default to sector-centered implementation

docs/src/changelog.md

@@ -33,6 +33,15 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   shorthand: a two-argument call is always the overlap, since the second argument cannot be
   disambiguated between a ket and an operator (undecidable for density matrices, where states
   and operators share a representation). ([#436](https://github.com/QuantumKitHub/MPSKit.jl/pull/436))
+- `transfer_spectrum` now computes the spectrum for all sectors of the transfer space at once,
+  returning a `TensorKit.SectorVector` that can be indexed per sector.
+  The `sector` keyword is removed; a specific selection of sectors (with per-sector counts) can be
+  requested by passing `howmany` as an `AbstractDict`/iterable of `sector => count` pairs.
+  The `below` state and the eigensolver algorithm are optional positional arguments, and the
+  algorithm can be resolved per sector. The Krylov dimension adapts to the number of values requested
+  in each sector, controlled by the new `oversampling` and `oversampling_factor` keywords.
+  Accordingly, `marek_gap` and `correlation_length` now return a `TensorKit.SectorDict` of
+  per-sector results by default; pass `sector = ...` to obtain a single sector's result as before.
 
 ### Deprecated
 

docs/src/examples/classic2d/1.hard-hexagon/index.md

@@ -12,7 +12,7 @@ EditURL = "../../../../../examples/classic2d/1.hard-hexagon/main.jl"
 
 Tensor networks are a natural way to do statistical mechanics on a lattice.
 As an example of this we will extract the central charge of the hard hexagon model.
-This model is known to have central charge 0.8, and has very peculiar non-local (anyonic) symmetries.
+This model is known to have central charge ``0.8``, and has very peculiar non-local (anyonic) symmetries.
 Because TensorKit supports anyonic symmetries, so does MPSKit.
 To follow the tutorial you need the following packages.
 
@@ -23,7 +23,8 @@ using MPSKit, MPSKitModels, TensorKit, Plots, Polynomials
 The [hard hexagon model](https://en.wikipedia.org/wiki/Hard_hexagon_model) is a 2-dimensional lattice model of a gas, where particles are allowed to be on the vertices of a triangular lattice, but no two particles may be adjacent.
 This can be encoded in a transfer matrix with a local MPO tensor using anyonic symmetries, and the resulting MPO has been implemented in MPSKitModels.
 
-In order to use these anyonic symmetries, we need to generalise the notion of the bond dimension and define how it interacts with the symmetry. Thus, we implement away of converting integers to symmetric spaces of the given dimension, which provides a crude guess for how the final MPS would distribute its Schmidt spectrum.
+In order to use these anyonic symmetries, we need to generalise the notion of the bond dimension and define how it interacts with the symmetry.
+Thus, we implement away of converting integers to symmetric spaces of the given dimension, which provides a crude guess for how the final MPS would distribute its Schmidt spectrum.
 
 ````julia
 mpo = hard_hexagon()
@@ -53,21 +54,24 @@ V = virtual_space(D)
 ) # use non-hermitian eigensolver
 F = real(expectation_value(ψ, mpo))
 S = real(first(entropy(ψ)))
-ξ = correlation_length(ψ)
+ξ = correlation_length(ψ; sector = leftunit(ψ))
 println("F = $F\tS = $S\tξ = $ξ")
 ````
 
 ````
-F = 0.8839037051703849	S = 1.280782962202701	ξ = 13.849682582739709
+F = 0.8839037051703852	S = 1.2807829621826905	ξ = 13.849682581482702
 
 ````
 
 ## The scaling hypothesis
 
-The dominant eigenvector is of course only an approximation. The finite bond dimension enforces a finite correlation length, which effectively introduces a length scale in the system. This can be exploited to formulate a scaling hypothesis [pollmann2009](@cite), which in turn allows to extract the central charge.
+The dominant eigenvector is of course only an approximation.
+The finite bond dimension enforces a finite correlation length, which effectively introduces a length scale in the system.
+This can be exploited to formulate a scaling hypothesis [pollmann2009](@cite), which in turn allows to extract the central charge.
 
-First we need to know the entropy and correlation length at a bunch of different bond dimensions. Our approach will be to re-use the previous approximated dominant eigenvector, and then expanding its bond dimension and re-running VUMPS.
-According to the scaling hypothesis we should have ``S \propto \frac{c}{6} log(ξ)``. Therefore we should find ``c`` using
+First we need to know the entropy and correlation length at a bunch of different bond dimensions.
+Our approach will be to re-use the previous approximated dominant eigenvector, and then expanding its bond dimension and re-running VUMPS.
+According to the scaling hypothesis we should have ``S ∝ \frac{c}{6} log(ξ)``. Therefore we should find ``c`` using
 
 ````julia
 function scaling_simulations(
@@ -77,16 +81,17 @@ function scaling_simulations(
     entropies = similar(Ds, Float64)
     correlations = similar(Ds, Float64)
     alg = VUMPS(; verbosity, tol, alg_eigsolve)
+    sector = unit(sectortype(mpo)) # dominant correlation functions are in the trivial sector
 
     ψ, envs, = leading_boundary(ψ₀, mpo, alg)
     entropies[1] = real(entropy(ψ)[1])
-    correlations[1] = correlation_length(ψ)
+    correlations[1] = correlation_length(ψ; sector)
 
     for (i, d) in enumerate(diff(Ds))
         ψ, envs = changebonds(ψ, mpo, OptimalExpand(; trscheme = truncrank(d)), envs)
         ψ, envs, = leading_boundary(ψ, mpo, alg, envs)
         entropies[i + 1] = real(entropy(ψ)[1])
-        correlations[i + 1] = correlation_length(ψ)
+        correlations[i + 1] = correlation_length(ψ; sector)
     end
     return entropies, correlations
 end
@@ -100,7 +105,7 @@ c = f.coeffs[2]
 ````
 
 ````
-0.8025237309498792
+0.802524639544328
 ````
 
 ````julia

docs/src/examples/classic2d/1.hard-hexagon/main.ipynb

@@ -10,7 +10,7 @@
     "\n",
     "Tensor networks are a natural way to do statistical mechanics on a lattice.\n",
     "As an example of this we will extract the central charge of the hard hexagon model.\n",
-    "This model is known to have central charge 0.8, and has very peculiar non-local (anyonic) symmetries.\n",
+    "This model is known to have central charge $0.8$, and has very peculiar non-local (anyonic) symmetries.\n",
     "Because TensorKit supports anyonic symmetries, so does MPSKit.\n",
     "To follow the tutorial you need the following packages."
    ]
@@ -31,7 +31,8 @@
     "The [hard hexagon model](https://en.wikipedia.org/wiki/Hard_hexagon_model) is a 2-dimensional lattice model of a gas, where particles are allowed to be on the vertices of a triangular lattice, but no two particles may be adjacent.\n",
     "This can be encoded in a transfer matrix with a local MPO tensor using anyonic symmetries, and the resulting MPO has been implemented in MPSKitModels.\n",
     "\n",
-    "In order to use these anyonic symmetries, we need to generalise the notion of the bond dimension and define how it interacts with the symmetry. Thus, we implement away of converting integers to symmetric spaces of the given dimension, which provides a crude guess for how the final MPS would distribute its Schmidt spectrum."
+    "In order to use these anyonic symmetries, we need to generalise the notion of the bond dimension and define how it interacts with the symmetry.\n",
+    "Thus, we implement away of converting integers to symmetric spaces of the given dimension, which provides a crude guess for how the final MPS would distribute its Schmidt spectrum."
    ]
   },
   {
@@ -77,7 +78,7 @@
     ") # use non-hermitian eigensolver\n",
     "F = real(expectation_value(ψ, mpo))\n",
     "S = real(first(entropy(ψ)))\n",
-    "ξ = correlation_length(ψ)\n",
+    "ξ = correlation_length(ψ; sector = leftunit(ψ))\n",
     "println(\"F = $F\\tS = $S\\tξ = $ξ\")"
    ]
   },
@@ -87,10 +88,13 @@
    "source": [
     "## The scaling hypothesis\n",
     "\n",
-    "The dominant eigenvector is of course only an approximation. The finite bond dimension enforces a finite correlation length, which effectively introduces a length scale in the system. This can be exploited to formulate a scaling hypothesis [pollmann2009](@cite), which in turn allows to extract the central charge.\n",
+    "The dominant eigenvector is of course only an approximation.\n",
+    "The finite bond dimension enforces a finite correlation length, which effectively introduces a length scale in the system.\n",
+    "This can be exploited to formulate a scaling hypothesis [pollmann2009](@cite), which in turn allows to extract the central charge.\n",
     "\n",
-    "First we need to know the entropy and correlation length at a bunch of different bond dimensions. Our approach will be to re-use the previous approximated dominant eigenvector, and then expanding its bond dimension and re-running VUMPS.\n",
-    "According to the scaling hypothesis we should have $S \\propto \\frac{c}{6} log(ξ)$. Therefore we should find $c$ using"
+    "First we need to know the entropy and correlation length at a bunch of different bond dimensions.\n",
+    "Our approach will be to re-use the previous approximated dominant eigenvector, and then expanding its bond dimension and re-running VUMPS.\n",
+    "According to the scaling hypothesis we should have $S ∝ \\frac{c}{6} log(ξ)$. Therefore we should find $c$ using"
    ]
   },
   {
@@ -106,16 +110,17 @@
     "    entropies = similar(Ds, Float64)\n",
     "    correlations = similar(Ds, Float64)\n",
     "    alg = VUMPS(; verbosity, tol, alg_eigsolve)\n",
+    "    sector = unit(sectortype(mpo)) # dominant correlation functions are in the trivial sector\n",
     "\n",
     "    ψ, envs, = leading_boundary(ψ₀, mpo, alg)\n",
     "    entropies[1] = real(entropy(ψ)[1])\n",
-    "    correlations[1] = correlation_length(ψ)\n",
+    "    correlations[1] = correlation_length(ψ; sector)\n",
     "\n",
     "    for (i, d) in enumerate(diff(Ds))\n",
     "        ψ, envs = changebonds(ψ, mpo, OptimalExpand(; trscheme = truncrank(d)), envs)\n",
     "        ψ, envs, = leading_boundary(ψ, mpo, alg, envs)\n",
     "        entropies[i + 1] = real(entropy(ψ)[1])\n",
-    "        correlations[i + 1] = correlation_length(ψ)\n",
+    "        correlations[i + 1] = correlation_length(ψ; sector)\n",
     "    end\n",
     "    return entropies, correlations\n",
     "end\n",

docs/src/examples/quantum1d/5.haldane-spt/index.md

@@ -128,7 +128,7 @@ E_plus = expectation_value(ψ_plus, H)
 ````
 
 ````
--1.4014193313393009 - 5.2545708620134027e-17im
+-1.4014193313393009 - 3.851708855717825e-17im
 ````
 
 ````julia
@@ -139,7 +139,7 @@ E_minus = expectation_value(ψ_minus, H)
 ````
 
 ````
--1.401483973963084 - 2.875923408269285e-17im
+-1.4014839739630844 - 5.800167584873572e-17im
 ````
 
 ````julia
@@ -188,11 +188,12 @@ println("S_plus = $S_plus")
 ````
 
 ````
-S_minus + log(2) = 1.5486227235421324
-S_plus = 1.5450323530571919
+S_minus + log(2) = 1.548622723541372
+S_plus = 1.5450323530299226
 
 ````
 
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+

docs/src/examples/quantum1d/5.haldane-spt/main.ipynb

@@ -213,4 +213,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 3
-}
+}

docs/src/examples/quantum1d/8.bose-hubbard/index.md

@@ -154,9 +154,9 @@ println("Energy: ", expectation_value(ground_state, hamiltonian))
 ````
 
 ````
-[ Info: VUMPS init:	obj = +4.962471958690e-01	err = 5.8876e-01
-[ Info: VUMPS conv 50:	obj = -6.757777651150e-01	err = 9.9482617748e-07	time = 6.72 sec
-Energy: -0.6757777651149999 - 6.246015340721835e-17im
+[ Info: VUMPS init:	obj = +3.450895756898e-01	err = 6.6224e-01
+[ Info: VUMPS conv 72:	obj = -6.756981551609e-01	err = 9.5394492271e-07	time = 2.91 sec
+Energy: -0.6756981551608794 - 1.967768893425733e-17im
 
 ````
 
@@ -257,7 +257,7 @@ plot!(
 )
 
 scatter!(
-    p, Ds, correlation_length.(states),
+    p, Ds, map(ψ -> correlation_length(ψ; sector = leftunit(ψ)), states),
     ylabel = "Correlation length", xlabel = "Bond dimension",
     xscale = :log10, yscale = :log10,
     inset = bbox(0.2, 0.51, 0.25, 0.25),
@@ -285,13 +285,13 @@ quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_o
 
 ````
 7-element Vector{Float64}:
- 0.30974277207656425
- 0.28814775930068737
- 0.2702000951730164
- 0.25712728156142256
- 0.2468538617129866
- 0.2353979140629328
- 0.2279965552408022
+ 0.31098779070601257
+ 0.2881478589434881
+ 0.2702000230423913
+ 0.25712728516508654
+ 0.24685385948652017
+ 0.23539753899204166
+ 0.22799654088348645
 ````
 
 We may now also visualize the momentum distribution function, which is obtained as the
@@ -466,3 +466,4 @@ using what we have learnt in this tutorial.
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+