diff --git a/Project.toml b/Project.toml
index 060484f18..5f6ecb111 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,7 +21,6 @@ OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
Preferences = "21216c6a-2e73-6563-6e65-726566657250"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
TensorKitManifolds = "11fa318c-39cb-4a83-b1ed-cdc7ba1e3684"
TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
@@ -30,9 +29,13 @@ VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
[weakdeps]
Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
[extensions]
MPSKitAdaptExt = "Adapt"
+MPSKitMakieExt = "Makie"
+MPSKitPlotsExt = "RecipesBase"
[compat]
Accessors = "0.1"
@@ -44,6 +47,7 @@ HalfIntegers = "1.6.0"
KrylovKit = "0.8.3, 0.9.2, 0.10"
LinearAlgebra = "1.6"
LoggingExtras = "1"
+Makie = "0.24.11"
MatrixAlgebraKit = "0.6.5"
OhMyThreads = "0.7, 0.8"
OptimKit = "0.4, 0.5"
diff --git a/docs/Project.toml b/docs/Project.toml
index 663381da8..6981f9dc2 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -8,7 +8,6 @@ DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index dff55f73b..d608cfead 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -129,6 +129,9 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
`expectation_value(::InfiniteMPS, ::InfiniteMPO)` forwards here. ([#508](https://github.com/QuantumKitHub/MPSKit.jl/pull/508))
- `InfiniteMultilineMPS`/`FiniteMultilineMPS` and `InfiniteMultilineMPO`/`FiniteMultilineMPO`
aliases for multiline objects with infinite or finite lines. ([#508](https://github.com/QuantumKitHub/MPSKit.jl/pull/508))
+- Makie.jl plotting functions `entanglement_plot` and `transfer_plot`, and their mutating
+ counterparts `entanglement_plot!` and `transfer_plot!`, through a new `MPSKitMakieExt` package
+ extension that is loaded with any Makie backend (e.g. `using CairoMakie`). ([#428](https://github.com/QuantumKitHub/MPSKit.jl/pull/428))
### Changed
@@ -200,6 +203,10 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
### Deprecated
- `DynamicTols.updatetol(alg, iter, ϵ)`; use `DynamicTols.adapt_solver(alg; iter, g_global = ϵ)` instead. ([#455](https://github.com/QuantumKitHub/MPSKit.jl/pull/455))
+- The Plots.jl recipes `entanglementplot` and `transferplot`, in favor of the Makie.jl based
+ `entanglement_plot` and `transfer_plot`. RecipesBase is now a weak dependency, so these are
+ only available once Plots.jl is loaded. ([#428](https://github.com/QuantumKitHub/MPSKit.jl/pull/428))
+
### Removed
- The custom `MPSTensor` constructors; construct `TensorMap`s with explicit `TensorMapSpace`s instead. ([#475](https://github.com/QuantumKitHub/MPSKit.jl/pull/475))
@@ -232,6 +239,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
matching tensors of a line, and accepts any row due to periodicity.
`checkbounds(Bool, ::PeriodicArray, I...)` now also accepts any index, matching `getindex`. ([#508](https://github.com/QuantumKitHub/MPSKit.jl/pull/508))
- `size`/`axes` for a `CView` over a `Multiline` with finite lines were missing. ([#508](https://github.com/QuantumKitHub/MPSKit.jl/pull/508))
+- `transferplot` shifted the angles the wrong way for a nonzero `thetaorigin`, moving eigenvalues
+ outside the plotted range. ([#428](https://github.com/QuantumKitHub/MPSKit.jl/pull/428))
- `excitations(::InfiniteMPO, ::QuasiparticleAnsatz, ::InfiniteQP, lenvs, renvs)` referenced `H_eff` before assigning. ([#489](https://github.com/QuantumKitHub/MPSKit.jl/pull/489))
- `Base.:+`/`-` on `FiniteMPS` returned a wrong state for near-parallel operands carried by
different tensor networks, e.g. `norm(E₀ * gs - H * gs)` coming out as `2 * norm(gs) * E₀`
diff --git a/docs/src/examples/dynamics/ising-dqpt/index.md b/docs/src/examples/dynamics/ising-dqpt/index.md
index c7d406254..e1dfed23b 100644
--- a/docs/src/examples/dynamics/ising-dqpt/index.md
+++ b/docs/src/examples/dynamics/ising-dqpt/index.md
@@ -111,7 +111,7 @@ dot(ψ₀, ψ₀)
````
````
-0.9999999999999999 + 6.62289120209154e-17im
+1.0000000000000002 - 4.9391993181912504e-17im
````
so the Loschmidt echo takes on the pleasant form
diff --git a/docs/src/examples/dynamics/ising-dqpt/main.ipynb b/docs/src/examples/dynamics/ising-dqpt/main.ipynb
index 7093b71a5..e89e2c865 100644
--- a/docs/src/examples/dynamics/ising-dqpt/main.ipynb
+++ b/docs/src/examples/dynamics/ising-dqpt/main.ipynb
@@ -390,7 +390,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Julia 1.13.1",
+ "display_name": "Julia 1.13.3",
"language": "julia",
"name": "julia-1.13"
},
@@ -398,7 +398,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
- "version": "1.13.1"
+ "version": "1.13.3"
}
},
"nbformat": 4,
diff --git a/docs/src/examples/dynamics/xy-finiteT/figure-1.png b/docs/src/examples/dynamics/xy-finiteT/figure-1.png
index be692709a..a3de820f7 100644
Binary files a/docs/src/examples/dynamics/xy-finiteT/figure-1.png and b/docs/src/examples/dynamics/xy-finiteT/figure-1.png differ
diff --git a/docs/src/examples/dynamics/xy-finiteT/figure-2.png b/docs/src/examples/dynamics/xy-finiteT/figure-2.png
index 508b94015..27d6ba5d7 100644
Binary files a/docs/src/examples/dynamics/xy-finiteT/figure-2.png and b/docs/src/examples/dynamics/xy-finiteT/figure-2.png differ
diff --git a/docs/src/examples/dynamics/xy-finiteT/figure-3.png b/docs/src/examples/dynamics/xy-finiteT/figure-3.png
index 9df578602..722da207f 100644
Binary files a/docs/src/examples/dynamics/xy-finiteT/figure-3.png and b/docs/src/examples/dynamics/xy-finiteT/figure-3.png differ
diff --git a/docs/src/examples/dynamics/xy-finiteT/figure-4.png b/docs/src/examples/dynamics/xy-finiteT/figure-4.png
index 701d842e0..22e977a66 100644
Binary files a/docs/src/examples/dynamics/xy-finiteT/figure-4.png and b/docs/src/examples/dynamics/xy-finiteT/figure-4.png differ
diff --git a/docs/src/examples/dynamics/xy-finiteT/figure-5.png b/docs/src/examples/dynamics/xy-finiteT/figure-5.png
index cfd09fdbc..d4838a8fa 100644
Binary files a/docs/src/examples/dynamics/xy-finiteT/figure-5.png and b/docs/src/examples/dynamics/xy-finiteT/figure-5.png differ
diff --git a/docs/src/examples/dynamics/xy-finiteT/figure-6.png b/docs/src/examples/dynamics/xy-finiteT/figure-6.png
index 37ddfbd45..8f29eaaab 100644
Binary files a/docs/src/examples/dynamics/xy-finiteT/figure-6.png and b/docs/src/examples/dynamics/xy-finiteT/figure-6.png differ
diff --git a/docs/src/examples/dynamics/xy-finiteT/index.md b/docs/src/examples/dynamics/xy-finiteT/index.md
index 04e0d6c1d..ccf256762 100644
--- a/docs/src/examples/dynamics/xy-finiteT/index.md
+++ b/docs/src/examples/dynamics/xy-finiteT/index.md
@@ -15,7 +15,7 @@ using .ExampleModels
using MPSKit: infinite_temperature_density_matrix
using QuadGK: quadgk
using SpecialFunctions: ellipe
-using Plots
+using CairoMakie
using LinearAlgebra
using BenchmarkFreeFermions
````
@@ -119,7 +119,7 @@ println("Exact (N=Inf):\t", groundstate_energy(J, Inf))
````
````
-Numerical: -0.15637803047011115
+Numerical: -0.1563780304701102
Exact (N=32): -0.15637803047254015
Exact (N=Inf): -0.15915494309189535
@@ -207,23 +207,39 @@ Z_taylor = map(Iterators.product(βs, expansion_orders)) do (β, expansion_order
end
F_taylor = -(1 ./ βs) .* Z_taylor
-p_taylor = let
- labels = reshape(
- map(expansion_orders) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, Z_analytic; label = "analytic", title = "Partition function",
- xlabel = "β", ylabel = "Z(β)"
- )
- plot!(p1, βs, Z_taylor; label = labels)
- p2 = plot(
- βs, F_analytic; label = "analytic", title = "Free energy",
- xlabel = "β", ylabel = "F(β)"
- )
- plot!(p2, βs, F_taylor; label = labels)
- plot(p1, p2)
+# Makie skips NaN, but picks wrong limits for infinite or (on a log scale) zero values
+nan_unless(f, ys) = [f(y) ? y : NaN for y in ys]
+
+# Plot the errors of several methods side by side, given as `label => (βs, ΔZ, ΔF)`.
+function plot_errors(errors; position = :lt, yscale = identity)
+ fig = Figure(; size = (800, 400))
+ axZ = Axis(fig[1, 1]; title = "Partition function error", xlabel = "β", ylabel = "ΔZ(β)", yscale)
+ axF = Axis(fig[1, 2]; title = "Free energy error", xlabel = "β", ylabel = "ΔF(β)", yscale)
+ valid = yscale === identity ? isfinite : y -> isfinite(y) && y > 0
+ for (label, (β, ΔZ, ΔF)) in errors
+ lines!(axZ, β, nan_unless(valid, ΔZ); label)
+ lines!(axF, β, nan_unless(valid, ΔF); label)
+ end
+ axislegend(axZ; position)
+ return fig
+end
+
+taylor_errors(Z, F, orders) = [
+ "Taylor N=$N" => (βs, abs.(Z[:, i] .- Z_analytic), abs.(F[:, i] .- F_analytic))
+ for (i, N) in enumerate(orders)
+]
+
+p_taylor = let fig = Figure(; size = (800, 400))
+ axZ = Axis(fig[1, 1]; title = "Partition function", xlabel = "β", ylabel = "Z(β)")
+ axF = Axis(fig[1, 2]; title = "Free energy", xlabel = "β", ylabel = "F(β)")
+ lines!(axZ, βs, Z_analytic; label = "analytic")
+ lines!(axF, βs, nan_unless(isfinite, F_analytic); label = "analytic")
+ for (i, N) in enumerate(expansion_orders)
+ lines!(axZ, βs, Z_taylor[:, i]; label = "Taylor N=$N")
+ lines!(axF, βs, nan_unless(isfinite, F_taylor[:, i]); label = "Taylor N=$N")
+ end
+ axislegend(axZ; position = :lt)
+ fig
end
````
@@ -261,23 +277,7 @@ expansion_orders = 2:3
Z_taylor = Z_taylor[:, 2:end]
F_taylor = F_taylor[:, 2:end]
-p_taylor_diff = let
- labels = reshape(
- map(expansion_orders) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, abs.(Z_taylor .- Z_analytic);
- label = labels, title = "Partition function error",
- xlabel = "β", ylabel = "ΔZ(β)", legend = :topleft
- )
- p2 = plot(
- βs, abs.(F_taylor .- F_analytic); label = labels,
- xlabel = "β", ylabel = "ΔF(β)", title = "Free energy error", legend = :topleft
- )
- plot(p1, p2)
-end
+p_taylor_diff = plot_errors(taylor_errors(Z_taylor, F_taylor, expansion_orders))
````
@@ -312,23 +312,7 @@ Z_taylor2 = map(Iterators.product(βs, expansion_orders)) do (β, expansion_orde
end
F_taylor2 = -(1 ./ βs) .* Z_taylor2
-p_taylor2_diff = let
- labels = reshape(
- map(expansion_orders[2:end]) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, abs.(Z_taylor2 .- Z_analytic);
- label = labels, title = "Partition function error",
- xlabel = "β", ylabel = "ΔZ(β)", legend = :topleft
- )
- p2 = plot(
- βs, abs.(F_taylor2 .- F_analytic); label = labels,
- xlabel = "β", ylabel = "ΔF(β)", title = "Free energy error", legend = :topleft
- )
- plot(p1, p2)
-end
+p_taylor2_diff = plot_errors(taylor_errors(Z_taylor2, F_taylor2, expansion_orders))
````
@@ -380,32 +364,12 @@ for i in 3:length(βs)
end
F_mpo_mul = -(1 ./ βs) .* Z_mpo_mul
-p_mpo_mul_diff = let
- labels = reshape(
- map(expansion_orders) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, abs.(Z_taylor2 .- Z_analytic);
- label = labels, title = "Partition function error",
- xlabel = "β", ylabel = "ΔZ(β)", legend = :bottomright, yscale = :log10
- )
- plot!(
- p1, βs, abs.(Z_mpo_mul .- Z_analytic);
- label = "MPO multiplication"
- )
- p2 = plot(
- βs, abs.(F_taylor2 .- F_analytic); label = labels,
- xlabel = "β", ylabel = "ΔF(β)", title = "Free energy error", legend = nothing,
- yscale = :log10
- )
- plot!(
- p2, βs, abs.(F_mpo_mul .- F_analytic);
- label = "MPO multiplication"
- )
- plot(p1, p2)
-end
+p_mpo_mul_diff = plot_errors(
+ [
+ taylor_errors(Z_taylor2, F_taylor2, expansion_orders);
+ "MPO multiplication" => (βs, abs.(Z_mpo_mul .- Z_analytic), abs.(F_mpo_mul .- F_analytic))
+ ]; position = :rb, yscale = log10
+)
````
@@ -464,28 +428,13 @@ for i in 2:length(βs_exp)
end
F_mpo_mul_exp = -(1 ./ βs_exp) .* Z_mpo_mul_exp
-p_mpo_mul_exp_diff = let
- labels = reshape(
- map(expansion_orders[2:end]) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, abs.(Z_taylor2 .- Z_analytic);
- label = labels, title = "Partition function error", xlabel = "β", ylabel = "ΔZ(β)",
- legend = :bottomright, yscale = :log10
- )
- plot!(p1, βs, abs.(Z_mpo_mul .- Z_analytic); label = "MPO multiplication")
- plot!(p1, βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp); label = "MPO multiplication exp")
-
- p2 = plot(
- βs, abs.(F_taylor2 .- F_analytic); label = labels, xlabel = "β", ylabel = "ΔF(β)",
- title = "Free energy error", legend = nothing, yscale = :log10
- )
- plot!(p2, βs, abs.(F_mpo_mul .- F_analytic); label = "MPO multiplication")
- plot!(p2, βs_exp, abs.(F_mpo_mul_exp .- F_analytic_exp); label = "MPO multiplication exp")
- plot(p1, p2)
-end
+p_mpo_mul_exp_diff = plot_errors(
+ [
+ taylor_errors(Z_taylor2, F_taylor2, expansion_orders);
+ "MPO multiplication" => (βs, abs.(Z_mpo_mul .- Z_analytic), abs.(F_mpo_mul .- F_analytic))
+ "MPO multiplication exp" => (βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp), abs.(F_mpo_mul_exp .- F_analytic_exp))
+ ]; position = :rb, yscale = log10
+)
````
@@ -531,31 +480,14 @@ for i in 2:length(βs)
end
F_tdvp = -(1 ./ βs) .* Z_tdvp
-p_mpo_mul_diff = let
- labels = reshape(
- map(expansion_orders) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, abs.(Z_taylor2 .- Z_analytic); label = labels,
- title = "Partition function error", xlabel = "β", ylabel = "ΔZ(β)",
- legend = :bottomright, yscale = :log10
- )
- plot!(p1, βs, abs.(Z_mpo_mul .- Z_analytic); label = "MPO multiplication")
- plot!(p1, βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp); label = "MPO multiplication exp")
- plot!(p1, βs, abs.(Z_tdvp .- Z_analytic); label = "TDVP")
-
- p2 = plot(
- βs, abs.(F_taylor2 .- F_analytic); label = labels, xlabel = "β", ylabel = "ΔF(β)",
- title = "Free energy error", legend = nothing, yscale = :log10
- )
- plot!(p2, βs, abs.(F_mpo_mul .- F_analytic); label = "MPO multiplication")
- plot!(p2, βs_exp, abs.(F_mpo_mul_exp .- F_analytic_exp); label = "MPO multiplication exp")
- plot!(p2, βs, abs.(F_tdvp .- F_analytic); label = "TDVP")
-
- plot(p1, p2)
-end
+p_mpo_mul_diff = plot_errors(
+ [
+ taylor_errors(Z_taylor2, F_taylor2, expansion_orders);
+ "MPO multiplication" => (βs, abs.(Z_mpo_mul .- Z_analytic), abs.(F_mpo_mul .- F_analytic))
+ "MPO multiplication exp" => (βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp), abs.(F_mpo_mul_exp .- F_analytic_exp))
+ "TDVP" => (βs, abs.(Z_tdvp .- Z_analytic), abs.(F_tdvp .- F_analytic))
+ ]; position = :rb, yscale = log10
+)
````
diff --git a/docs/src/examples/dynamics/xy-finiteT/main.ipynb b/docs/src/examples/dynamics/xy-finiteT/main.ipynb
index fed2a1e4f..c2e557d63 100644
--- a/docs/src/examples/dynamics/xy-finiteT/main.ipynb
+++ b/docs/src/examples/dynamics/xy-finiteT/main.ipynb
@@ -119,7 +119,7 @@
"using MPSKit: infinite_temperature_density_matrix\n",
"using QuadGK: quadgk\n",
"using SpecialFunctions: ellipe\n",
- "using Plots\n",
+ "using CairoMakie\n",
"using LinearAlgebra\n",
"using BenchmarkFreeFermions"
]
@@ -354,23 +354,39 @@
"end\n",
"F_taylor = -(1 ./ βs) .* Z_taylor\n",
"\n",
- "p_taylor = let\n",
- " labels = reshape(\n",
- " map(expansion_orders) do N\n",
- " return \"Taylor N=$N\"\n",
- " end, 1, :\n",
- " )\n",
- " p1 = plot(\n",
- " βs, Z_analytic; label = \"analytic\", title = \"Partition function\",\n",
- " xlabel = \"β\", ylabel = \"Z(β)\"\n",
- " )\n",
- " plot!(p1, βs, Z_taylor; label = labels)\n",
- " p2 = plot(\n",
- " βs, F_analytic; label = \"analytic\", title = \"Free energy\",\n",
- " xlabel = \"β\", ylabel = \"F(β)\"\n",
- " )\n",
- " plot!(p2, βs, F_taylor; label = labels)\n",
- " plot(p1, p2)\n",
+ "# Makie skips NaN, but picks wrong limits for infinite or (on a log scale) zero values\n",
+ "nan_unless(f, ys) = [f(y) ? y : NaN for y in ys]\n",
+ "\n",
+ "# Plot the errors of several methods side by side, given as `label => (βs, ΔZ, ΔF)`.\n",
+ "function plot_errors(errors; position = :lt, yscale = identity)\n",
+ " fig = Figure(; size = (800, 400))\n",
+ " axZ = Axis(fig[1, 1]; title = \"Partition function error\", xlabel = \"β\", ylabel = \"ΔZ(β)\", yscale)\n",
+ " axF = Axis(fig[1, 2]; title = \"Free energy error\", xlabel = \"β\", ylabel = \"ΔF(β)\", yscale)\n",
+ " valid = yscale === identity ? isfinite : y -> isfinite(y) && y > 0\n",
+ " for (label, (β, ΔZ, ΔF)) in errors\n",
+ " lines!(axZ, β, nan_unless(valid, ΔZ); label)\n",
+ " lines!(axF, β, nan_unless(valid, ΔF); label)\n",
+ " end\n",
+ " axislegend(axZ; position)\n",
+ " return fig\n",
+ "end\n",
+ "\n",
+ "taylor_errors(Z, F, orders) = [\n",
+ " \"Taylor N=$N\" => (βs, abs.(Z[:, i] .- Z_analytic), abs.(F[:, i] .- F_analytic))\n",
+ " for (i, N) in enumerate(orders)\n",
+ "]\n",
+ "\n",
+ "p_taylor = let fig = Figure(; size = (800, 400))\n",
+ " axZ = Axis(fig[1, 1]; title = \"Partition function\", xlabel = \"β\", ylabel = \"Z(β)\")\n",
+ " axF = Axis(fig[1, 2]; title = \"Free energy\", xlabel = \"β\", ylabel = \"F(β)\")\n",
+ " lines!(axZ, βs, Z_analytic; label = \"analytic\")\n",
+ " lines!(axF, βs, nan_unless(isfinite, F_analytic); label = \"analytic\")\n",
+ " for (i, N) in enumerate(expansion_orders)\n",
+ " lines!(axZ, βs, Z_taylor[:, i]; label = \"Taylor N=$N\")\n",
+ " lines!(axF, βs, nan_unless(isfinite, F_taylor[:, i]); label = \"Taylor N=$N\")\n",
+ " end\n",
+ " axislegend(axZ; position = :lt)\n",
+ " fig\n",
"end"
]
},
@@ -421,23 +437,7 @@
"Z_taylor = Z_taylor[:, 2:end]\n",
"F_taylor = F_taylor[:, 2:end]\n",
"\n",
- "p_taylor_diff = let\n",
- " labels = reshape(\n",
- " map(expansion_orders) do N\n",
- " return \"Taylor N=$N\"\n",
- " end, 1, :\n",
- " )\n",
- " p1 = plot(\n",
- " βs, abs.(Z_taylor .- Z_analytic);\n",
- " label = labels, title = \"Partition function error\",\n",
- " xlabel = \"β\", ylabel = \"ΔZ(β)\", legend = :topleft\n",
- " )\n",
- " p2 = plot(\n",
- " βs, abs.(F_taylor .- F_analytic); label = labels,\n",
- " xlabel = \"β\", ylabel = \"ΔF(β)\", title = \"Free energy error\", legend = :topleft\n",
- " )\n",
- " plot(p1, p2)\n",
- "end"
+ "p_taylor_diff = plot_errors(taylor_errors(Z_taylor, F_taylor, expansion_orders))"
]
},
{
@@ -480,23 +480,7 @@
"end\n",
"F_taylor2 = -(1 ./ βs) .* Z_taylor2\n",
"\n",
- "p_taylor2_diff = let\n",
- " labels = reshape(\n",
- " map(expansion_orders[2:end]) do N\n",
- " return \"Taylor N=$N\"\n",
- " end, 1, :\n",
- " )\n",
- " p1 = plot(\n",
- " βs, abs.(Z_taylor2 .- Z_analytic);\n",
- " label = labels, title = \"Partition function error\",\n",
- " xlabel = \"β\", ylabel = \"ΔZ(β)\", legend = :topleft\n",
- " )\n",
- " p2 = plot(\n",
- " βs, abs.(F_taylor2 .- F_analytic); label = labels,\n",
- " xlabel = \"β\", ylabel = \"ΔF(β)\", title = \"Free energy error\", legend = :topleft\n",
- " )\n",
- " plot(p1, p2)\n",
- "end"
+ "p_taylor2_diff = plot_errors(taylor_errors(Z_taylor2, F_taylor2, expansion_orders))"
]
},
{
@@ -557,32 +541,12 @@
"end\n",
"F_mpo_mul = -(1 ./ βs) .* Z_mpo_mul\n",
"\n",
- "p_mpo_mul_diff = let\n",
- " labels = reshape(\n",
- " map(expansion_orders) do N\n",
- " return \"Taylor N=$N\"\n",
- " end, 1, :\n",
- " )\n",
- " p1 = plot(\n",
- " βs, abs.(Z_taylor2 .- Z_analytic);\n",
- " label = labels, title = \"Partition function error\",\n",
- " xlabel = \"β\", ylabel = \"ΔZ(β)\", legend = :bottomright, yscale = :log10\n",
- " )\n",
- " plot!(\n",
- " p1, βs, abs.(Z_mpo_mul .- Z_analytic);\n",
- " label = \"MPO multiplication\"\n",
- " )\n",
- " p2 = plot(\n",
- " βs, abs.(F_taylor2 .- F_analytic); label = labels,\n",
- " xlabel = \"β\", ylabel = \"ΔF(β)\", title = \"Free energy error\", legend = nothing,\n",
- " yscale = :log10\n",
- " )\n",
- " plot!(\n",
- " p2, βs, abs.(F_mpo_mul .- F_analytic);\n",
- " label = \"MPO multiplication\"\n",
- " )\n",
- " plot(p1, p2)\n",
- "end"
+ "p_mpo_mul_diff = plot_errors(\n",
+ " [\n",
+ " taylor_errors(Z_taylor2, F_taylor2, expansion_orders);\n",
+ " \"MPO multiplication\" => (βs, abs.(Z_mpo_mul .- Z_analytic), abs.(F_mpo_mul .- F_analytic))\n",
+ " ]; position = :rb, yscale = log10\n",
+ ")"
]
},
{
@@ -654,28 +618,13 @@
"end\n",
"F_mpo_mul_exp = -(1 ./ βs_exp) .* Z_mpo_mul_exp\n",
"\n",
- "p_mpo_mul_exp_diff = let\n",
- " labels = reshape(\n",
- " map(expansion_orders[2:end]) do N\n",
- " return \"Taylor N=$N\"\n",
- " end, 1, :\n",
- " )\n",
- " p1 = plot(\n",
- " βs, abs.(Z_taylor2 .- Z_analytic);\n",
- " label = labels, title = \"Partition function error\", xlabel = \"β\", ylabel = \"ΔZ(β)\",\n",
- " legend = :bottomright, yscale = :log10\n",
- " )\n",
- " plot!(p1, βs, abs.(Z_mpo_mul .- Z_analytic); label = \"MPO multiplication\")\n",
- " plot!(p1, βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp); label = \"MPO multiplication exp\")\n",
- "\n",
- " p2 = plot(\n",
- " βs, abs.(F_taylor2 .- F_analytic); label = labels, xlabel = \"β\", ylabel = \"ΔF(β)\",\n",
- " title = \"Free energy error\", legend = nothing, yscale = :log10\n",
- " )\n",
- " plot!(p2, βs, abs.(F_mpo_mul .- F_analytic); label = \"MPO multiplication\")\n",
- " plot!(p2, βs_exp, abs.(F_mpo_mul_exp .- F_analytic_exp); label = \"MPO multiplication exp\")\n",
- " plot(p1, p2)\n",
- "end"
+ "p_mpo_mul_exp_diff = plot_errors(\n",
+ " [\n",
+ " taylor_errors(Z_taylor2, F_taylor2, expansion_orders);\n",
+ " \"MPO multiplication\" => (βs, abs.(Z_mpo_mul .- Z_analytic), abs.(F_mpo_mul .- F_analytic))\n",
+ " \"MPO multiplication exp\" => (βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp), abs.(F_mpo_mul_exp .- F_analytic_exp))\n",
+ " ]; position = :rb, yscale = log10\n",
+ ")"
]
},
{
@@ -735,31 +684,14 @@
"end\n",
"F_tdvp = -(1 ./ βs) .* Z_tdvp\n",
"\n",
- "p_mpo_mul_diff = let\n",
- " labels = reshape(\n",
- " map(expansion_orders) do N\n",
- " return \"Taylor N=$N\"\n",
- " end, 1, :\n",
- " )\n",
- " p1 = plot(\n",
- " βs, abs.(Z_taylor2 .- Z_analytic); label = labels,\n",
- " title = \"Partition function error\", xlabel = \"β\", ylabel = \"ΔZ(β)\",\n",
- " legend = :bottomright, yscale = :log10\n",
- " )\n",
- " plot!(p1, βs, abs.(Z_mpo_mul .- Z_analytic); label = \"MPO multiplication\")\n",
- " plot!(p1, βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp); label = \"MPO multiplication exp\")\n",
- " plot!(p1, βs, abs.(Z_tdvp .- Z_analytic); label = \"TDVP\")\n",
- "\n",
- " p2 = plot(\n",
- " βs, abs.(F_taylor2 .- F_analytic); label = labels, xlabel = \"β\", ylabel = \"ΔF(β)\",\n",
- " title = \"Free energy error\", legend = nothing, yscale = :log10\n",
- " )\n",
- " plot!(p2, βs, abs.(F_mpo_mul .- F_analytic); label = \"MPO multiplication\")\n",
- " plot!(p2, βs_exp, abs.(F_mpo_mul_exp .- F_analytic_exp); label = \"MPO multiplication exp\")\n",
- " plot!(p2, βs, abs.(F_tdvp .- F_analytic); label = \"TDVP\")\n",
- "\n",
- " plot(p1, p2)\n",
- "end"
+ "p_mpo_mul_diff = plot_errors(\n",
+ " [\n",
+ " taylor_errors(Z_taylor2, F_taylor2, expansion_orders);\n",
+ " \"MPO multiplication\" => (βs, abs.(Z_mpo_mul .- Z_analytic), abs.(F_mpo_mul .- F_analytic))\n",
+ " \"MPO multiplication exp\" => (βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp), abs.(F_mpo_mul_exp .- F_analytic_exp))\n",
+ " \"TDVP\" => (βs, abs.(Z_tdvp .- Z_analytic), abs.(F_tdvp .- F_analytic))\n",
+ " ]; position = :rb, yscale = log10\n",
+ ")"
]
},
{
@@ -784,7 +716,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Julia 1.13.1",
+ "display_name": "Julia 1.13.3",
"language": "julia",
"name": "julia-1.13"
},
@@ -792,7 +724,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
- "version": "1.13.1"
+ "version": "1.13.3"
}
},
"nbformat": 4,
diff --git a/docs/src/examples/excitations/haldane/figure-1.png b/docs/src/examples/excitations/haldane/figure-1.png
index e170413cc..2eedc44e0 100644
Binary files a/docs/src/examples/excitations/haldane/figure-1.png and b/docs/src/examples/excitations/haldane/figure-1.png differ
diff --git a/docs/src/examples/excitations/haldane/figure-2.png b/docs/src/examples/excitations/haldane/figure-2.png
index 87ebf95e3..d8ae5a954 100644
Binary files a/docs/src/examples/excitations/haldane/figure-2.png and b/docs/src/examples/excitations/haldane/figure-2.png differ
diff --git a/docs/src/examples/excitations/haldane/figure-3.png b/docs/src/examples/excitations/haldane/figure-3.png
index 3caef8e80..2c2149c8f 100644
Binary files a/docs/src/examples/excitations/haldane/figure-3.png and b/docs/src/examples/excitations/haldane/figure-3.png differ
diff --git a/docs/src/examples/excitations/haldane/index.md b/docs/src/examples/excitations/haldane/index.md
index 3c37b066a..104e2ab4e 100644
--- a/docs/src/examples/excitations/haldane/index.md
+++ b/docs/src/examples/excitations/haldane/index.md
@@ -12,7 +12,7 @@ In this tutorial we will calculate the Haldane gap (the energy gap in the ``S =
To follow the tutorial you need the following packages:
````julia
-using MPSKit, TensorKit, Plots, Polynomials
+using MPSKit, TensorKit, CairoMakie, Polynomials
using .ExampleModels
using TensorKitTensors.SpinOperators: S_exchange
````
@@ -62,22 +62,25 @@ En_2, st_2 = excitations(H, QuasiparticleAnsatz(), ψ, envs; sector = SU2Irrep(2
````
````
-0.7989253589480334
+0.7989253589480464
````
We can go even further and doublecheck the claim that ``S = 1`` is an edge excitation, by plotting the energy density.
````julia
-p_density = plot(; xaxis = "position", yaxis = "energy density")
excited_1 = convert(FiniteMPS, st_1[1])
excited_2 = convert(FiniteMPS, st_2[1])
SS = -S_exchange(ComplexF64, SU2Irrep; spin = 1)
e₀ = [real(expectation_value(ψ, (i, i + 1) => SS)) for i in 1:(L - 1)]
e₁ = [real(expectation_value(excited_1, (i, i + 1) => SS)) for i in 1:(L - 1)]
e₂ = [real(expectation_value(excited_2, (i, i + 1) => SS)) for i in 1:(L - 1)]
-plot!(p_density, e₀; label = "S = 0")
-plot!(p_density, e₁; label = "S = 1")
-plot!(p_density, e₂; label = "S = 2")
+fig = Figure()
+ax = Axis(fig[1, 1]; xlabel = "position", ylabel = "energy density")
+lines!(ax, e₀; label = "S = 0")
+lines!(ax, e₁; label = "S = 1")
+lines!(ax, e₂; label = "S = 2")
+axislegend(ax)
+fig
````
@@ -101,13 +104,16 @@ f = fit(Ls .^ (-2), ΔEs, 1)
````
````
-0.4517340158583773
+0.4517340158585122
````
````julia
-p_size_extrapolation = plot(; xaxis = "L^(-2)", yaxis = "ΔE", xlims = (0, 0.015))
-plot!(p_size_extrapolation, Ls .^ (-2), ΔEs; seriestype = :scatter, label = "numerical")
-plot!(p_size_extrapolation, x -> f(x); label = "fit")
+fig = Figure()
+ax = Axis(fig[1, 1]; xlabel = "L^(-2)", ylabel = "ΔE", limits = ((0, 0.015), nothing))
+scatter!(ax, Ls .^ (-2), ΔEs; label = "numerical")
+lines!(ax, range(0, 0.015, 100), x -> f(x); label = "fit")
+axislegend(ax; position = :lt)
+fig
````
@@ -136,12 +142,12 @@ println("minimum @k = $(kspace[idx]):\t ΔE = $(ΔE)")
````
````
-minimum @k = 3.141592653589793: ΔE = 0.41047924868496444
+minimum @k = 3.141592653589793: ΔE = 0.4104792486030011
````
````julia
-plot(kspace, real.(Es); xaxis = "momentum", yaxis = "ΔE", label = "S = 1")
+series(kspace, permutedims(real.(Es)); axis = (; xlabel = "momentum", ylabel = "ΔE", title = "S = 1"))
````
diff --git a/docs/src/examples/excitations/haldane/main.ipynb b/docs/src/examples/excitations/haldane/main.ipynb
index 5137d13ae..439983d75 100644
--- a/docs/src/examples/excitations/haldane/main.ipynb
+++ b/docs/src/examples/excitations/haldane/main.ipynb
@@ -16,7 +16,7 @@
"metadata": {},
"outputs": [],
"source": [
- "using MPSKit, TensorKit, Plots, Polynomials\n",
+ "using MPSKit, TensorKit, CairoMakie, Polynomials\n",
"# Model constructors for the MPSKit documentation and examples.\n",
"# Adapted from MPSKitModels.jl v0.4.7, src/models/{hamiltonians,transfermatrices}.jl:\n",
"# https://github.com/QuantumKitHub/MPSKitModels.jl/tree/v0.4.7/src/models\n",
@@ -199,16 +199,19 @@
"metadata": {},
"outputs": [],
"source": [
- "p_density = plot(; xaxis = \"position\", yaxis = \"energy density\")\n",
"excited_1 = convert(FiniteMPS, st_1[1])\n",
"excited_2 = convert(FiniteMPS, st_2[1])\n",
"SS = -S_exchange(ComplexF64, SU2Irrep; spin = 1)\n",
"e₀ = [real(expectation_value(ψ, (i, i + 1) => SS)) for i in 1:(L - 1)]\n",
"e₁ = [real(expectation_value(excited_1, (i, i + 1) => SS)) for i in 1:(L - 1)]\n",
"e₂ = [real(expectation_value(excited_2, (i, i + 1) => SS)) for i in 1:(L - 1)]\n",
- "plot!(p_density, e₀; label = \"S = 0\")\n",
- "plot!(p_density, e₁; label = \"S = 1\")\n",
- "plot!(p_density, e₂; label = \"S = 2\")"
+ "fig = Figure()\n",
+ "ax = Axis(fig[1, 1]; xlabel = \"position\", ylabel = \"energy density\")\n",
+ "lines!(ax, e₀; label = \"S = 0\")\n",
+ "lines!(ax, e₁; label = \"S = 1\")\n",
+ "lines!(ax, e₂; label = \"S = 2\")\n",
+ "axislegend(ax)\n",
+ "fig"
]
},
{
@@ -237,9 +240,12 @@
"\n",
"f = fit(Ls .^ (-2), ΔEs, 1)\n",
"ΔE_extrapolated = f.coeffs[1]\n",
- "p_size_extrapolation = plot(; xaxis = \"L^(-2)\", yaxis = \"ΔE\", xlims = (0, 0.015))\n",
- "plot!(p_size_extrapolation, Ls .^ (-2), ΔEs; seriestype = :scatter, label = \"numerical\")\n",
- "plot!(p_size_extrapolation, x -> f(x); label = \"fit\")"
+ "fig = Figure()\n",
+ "ax = Axis(fig[1, 1]; xlabel = \"L^(-2)\", ylabel = \"ΔE\", limits = ((0, 0.015), nothing))\n",
+ "scatter!(ax, Ls .^ (-2), ΔEs; label = \"numerical\")\n",
+ "lines!(ax, range(0, 0.015, 100), x -> f(x); label = \"fit\")\n",
+ "axislegend(ax; position = :lt)\n",
+ "fig"
]
},
{
@@ -273,7 +279,7 @@
"\n",
"ΔE, idx = findmin(real.(Es))\n",
"println(\"minimum @k = $(kspace[idx]):\\t ΔE = $(ΔE)\")\n",
- "plot(kspace, real.(Es); xaxis = \"momentum\", yaxis = \"ΔE\", label = \"S = 1\")"
+ "series(kspace, permutedims(real.(Es)); axis = (; xlabel = \"momentum\", ylabel = \"ΔE\", title = \"S = 1\"))"
]
},
{
@@ -288,7 +294,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Julia 1.13.1",
+ "display_name": "Julia 1.13.3",
"language": "julia",
"name": "julia-1.13"
},
@@ -296,7 +302,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
- "version": "1.13.1"
+ "version": "1.13.3"
}
},
"nbformat": 4,
diff --git a/docs/src/examples/groundstates/bose-hubbard/figure-1.png b/docs/src/examples/groundstates/bose-hubbard/figure-1.png
index 7bd25fdea..84f873d29 100644
Binary files a/docs/src/examples/groundstates/bose-hubbard/figure-1.png and b/docs/src/examples/groundstates/bose-hubbard/figure-1.png differ
diff --git a/docs/src/examples/groundstates/bose-hubbard/figure-2.png b/docs/src/examples/groundstates/bose-hubbard/figure-2.png
index b8ea6a15c..6eb23376c 100644
Binary files a/docs/src/examples/groundstates/bose-hubbard/figure-2.png and b/docs/src/examples/groundstates/bose-hubbard/figure-2.png differ
diff --git a/docs/src/examples/groundstates/bose-hubbard/figure-3.png b/docs/src/examples/groundstates/bose-hubbard/figure-3.png
index 680f6e6bb..823ed98ee 100644
Binary files a/docs/src/examples/groundstates/bose-hubbard/figure-3.png and b/docs/src/examples/groundstates/bose-hubbard/figure-3.png differ
diff --git a/docs/src/examples/groundstates/bose-hubbard/figure-4.png b/docs/src/examples/groundstates/bose-hubbard/figure-4.png
index 7f22fb7f5..1e023fd0f 100644
Binary files a/docs/src/examples/groundstates/bose-hubbard/figure-4.png and b/docs/src/examples/groundstates/bose-hubbard/figure-4.png differ
diff --git a/docs/src/examples/groundstates/bose-hubbard/figure-5.png b/docs/src/examples/groundstates/bose-hubbard/figure-5.png
index 394cdb16e..5b55141ff 100644
Binary files a/docs/src/examples/groundstates/bose-hubbard/figure-5.png and b/docs/src/examples/groundstates/bose-hubbard/figure-5.png differ
diff --git a/docs/src/examples/groundstates/bose-hubbard/figure-6.png b/docs/src/examples/groundstates/bose-hubbard/figure-6.png
index 9cd7a5d2e..32cb27604 100644
Binary files a/docs/src/examples/groundstates/bose-hubbard/figure-6.png and b/docs/src/examples/groundstates/bose-hubbard/figure-6.png differ
diff --git a/docs/src/examples/groundstates/bose-hubbard/index.md b/docs/src/examples/groundstates/bose-hubbard/index.md
index b73d62760..d478d4d0a 100644
--- a/docs/src/examples/groundstates/bose-hubbard/index.md
+++ b/docs/src/examples/groundstates/bose-hubbard/index.md
@@ -11,11 +11,10 @@ using Markdown
using TensorKitTensors.BosonOperators: b_min, b_plus_b_min, b_min_b_plus, b_num
using MPSKit, TensorKit
using .ExampleModels
-using Plots, LaTeXStrings
+using CairoMakie
-theme(:wong)
-default(fontfamily = "Computer Modern", label = nothing, dpi = 100, framestyle = :box)
+set_theme!(theme_latexfonts())
````
# 1D Bose-Hubbard model
@@ -133,9 +132,9 @@ println("Energy: ", expectation_value(ground_state, hamiltonian))
````
````
-[ Info: VUMPS init: obj = +3.565540128712e-01 err = 5.8200e-01
-[ Info: VUMPS conv 35: obj = -6.757777651150e-01 err = 9.9874991150e-07 time = 6.10 sec
-Energy: -0.6757777651149788 + 3.9753559577548215e-17im
+[ Info: VUMPS init: obj = +4.712416457085e-01 err = 5.7986e-01
+[ Info: VUMPS conv 67: obj = -6.757777651162e-01 err = 8.2673319905e-07 time = 1.30 min
+Energy: -0.675777765116182 + 5.4594185126295844e-17im
````
@@ -170,8 +169,12 @@ In case of the Hamiltonian, it is not necessary to specify the indices as it spa
We can now plot the correlation function $⟨â^{†}_i â_j⟩$.
````julia
-plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw = 2, xlabel = "Site index", ylabel = "Correlation function", yscale = :log10)
-hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black)
+fig, ax, _ = lines(
+ map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50);
+ linewidth = 2, axis = (; xlabel = "Site index", ylabel = "Correlation function", yscale = log10)
+)
+hlines!(ax, abs2(expectation_value(ground_state, (0,) => a_op)); linestyle = :dash, color = :black)
+fig
````
@@ -182,8 +185,13 @@ Let us now shift our parameters to probe other phases.
````julia
ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500)
-plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw = 2, xlabel = "Site index", ylabel = "Correlation function", yscale = :log10, xscale = :log10)
-hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black)
+fig, ax, _ = lines(
+ map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100);
+ linewidth = 2,
+ axis = (; xlabel = "Site index", ylabel = "Correlation function", xscale = log10, yscale = log10)
+)
+hlines!(ax, abs2(expectation_value(ground_state, (0,) => a_op)); linestyle = :dash, color = :black)
+fig
````
@@ -214,31 +222,26 @@ Threads.@threads for idx in eachindex(Ds)
two_point_correlation[idx, :] .= real.(expectation_value(states[idx], (1, i) => a_op' ⊗ a_op) for i in 1:npoints)
end
-p = plot(
- framestyle = :box, ylabel = "Correlation function " * L"\langle a_i^{\dagger}a_j \rangle",
- xlabel = "Distance " * L"|i-j|", xscale = :log10, yscale = :log10,
- xticks = ([10, 100], ["10", "100"]),
- yticks = ([0.25, 0.5, 1.0], ["0.25", "0.5", "1.0"])
-)
-
-plot!(
- p, 2:npoints, two_point_correlation[:, 2:end]',
- lab = "D = " .* string.(permutedims(Ds)), lw = 2
+fig = Figure()
+ax = Axis(
+ fig[1, 1];
+ xlabel = L"\text{Distance } |i-j|",
+ ylabel = L"\text{Correlation function } \langle a_i^{\dagger}a_j \rangle",
+ xscale = log10, yscale = log10, xticks = [10, 100], yticks = [0.25, 0.5, 1.0]
)
+for (D, correlation) in zip(Ds, eachrow(two_point_correlation))
+ lines!(ax, 2:npoints, correlation[2:end]; label = "D = $D", linewidth = 2)
+end
+axislegend(ax)
-scatter!(
- 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),
- subplot = 2,
- xticks = (20:10:50, string.(20:10:50)),
- yticks = ([50, 100], string.([50, 100])),
- xlabelfontsize = 8,
- ylabelfontsize = 8,
- ylims = [20, 130],
- xlims = [15, 60]
+inset = Axis(
+ fig[1, 1]; width = Relative(0.25), height = Relative(0.25), halign = 0.27, valign = 0.32,
+ xlabel = "Bond dimension", ylabel = "Correlation length", xlabelsize = 10, ylabelsize = 10,
+ xscale = log10, yscale = log10, xticks = 20:10:50, yticks = [50, 100],
+ limits = ((15, 60), (20, 130))
)
+scatter!(inset, Ds, map(ψ -> correlation_length(ψ; sector = leftunit(ψ)), states))
+fig
````
@@ -252,13 +255,13 @@ quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_o
````
7-element Vector{Float64}:
- 0.3097427801948516
- 0.2881477638457307
- 0.2702001643897378
- 0.2571272792016688
- 0.2468538581332228
- 0.23539756507406545
- 0.22799654612658096
+ 0.30974277570736086
+ 0.2881477298949796
+ 0.27020006032912913
+ 0.2571272839417027
+ 0.24685386116892385
+ 0.23539777044006097
+ 0.22799656077322028
````
We may now also visualize the momentum distribution function, which is obtained as the Fourier transform of the single-particle density matrix.
@@ -303,7 +306,13 @@ momentum_distribution = map(
)
)
momentum_distribution = vcat(momentum_distribution...)'
-plot(ks, momentum_distribution, lab = "D = " .* string.(permutedims(Ds)), lw = 1.5, xlabel = "Momentum k", ylabel = L"\langle n_k \rangle", ylim = [0, 50])
+fig = Figure()
+ax = Axis(fig[1, 1]; xlabel = "Momentum k", ylabel = L"\langle n_k \rangle", limits = (nothing, (0, 50)))
+for (D, n_k) in zip(Ds, eachcol(momentum_distribution))
+ lines!(ax, ks, n_k; label = "D = $D", linewidth = 1.5)
+end
+axislegend(ax)
+fig
````
@@ -358,7 +367,13 @@ function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ = 1.0e-4, npoints = 1
energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted))
end
- return plot(phis, energies, lw = 2, xlabel = "Phase twist per site" * L"(\phi)", ylabel = "Ground state energy", title = "t = $t | μ = $mu | D = $D | cutoff = $cutoff")
+ return lines(
+ phis, energies; linewidth = 2,
+ axis = (;
+ xlabel = L"\text{Phase twist per site } (\phi)", ylabel = "Ground state energy",
+ title = "t = $t | μ = $mu | D = $D | cutoff = $cutoff",
+ )
+ )
end
superfluid_stiffness_profile(0.2, 0.3, 5, 4) # superfluid
@@ -387,7 +402,12 @@ Threads.@threads for (i, j) in collect(Iterators.product(eachindex(mus), eachind
order_parameters[i, j] = abs(expectation_value(state, 0 => a_op))
end
-heatmap(ts, mus, order_parameters, xlabel = L"t/U", ylabel = L"\mu/U", title = L"\langle \hat{a}_i \rangle")
+fig, ax, hm = heatmap(
+ ts, mus, permutedims(order_parameters);
+ axis = (; xlabel = L"t/U", ylabel = L"\mu/U", title = L"\langle \hat{a}_i \rangle")
+)
+Colorbar(fig[1, 2], hm)
+fig
````
diff --git a/docs/src/examples/groundstates/bose-hubbard/main.ipynb b/docs/src/examples/groundstates/bose-hubbard/main.ipynb
index 35dfe1fe8..138a30a7c 100644
--- a/docs/src/examples/groundstates/bose-hubbard/main.ipynb
+++ b/docs/src/examples/groundstates/bose-hubbard/main.ipynb
@@ -115,11 +115,10 @@
"end\n",
"\n",
"using .ExampleModels\n",
- "using Plots, LaTeXStrings\n",
+ "using CairoMakie\n",
"\n",
"\n",
- "theme(:wong)\n",
- "default(fontfamily = \"Computer Modern\", label = nothing, dpi = 100, framestyle = :box)"
+ "set_theme!(theme_latexfonts())"
]
},
{
@@ -297,8 +296,12 @@
"metadata": {},
"outputs": [],
"source": [
- "plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw = 2, xlabel = \"Site index\", ylabel = \"Correlation function\", yscale = :log10)\n",
- "hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black)"
+ "fig, ax, _ = lines(\n",
+ " map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50);\n",
+ " linewidth = 2, axis = (; xlabel = \"Site index\", ylabel = \"Correlation function\", yscale = log10)\n",
+ ")\n",
+ "hlines!(ax, abs2(expectation_value(ground_state, (0,) => a_op)); linestyle = :dash, color = :black)\n",
+ "fig"
]
},
{
@@ -317,8 +320,13 @@
"source": [
"ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500)\n",
"\n",
- "plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw = 2, xlabel = \"Site index\", ylabel = \"Correlation function\", yscale = :log10, xscale = :log10)\n",
- "hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black)"
+ "fig, ax, _ = lines(\n",
+ " map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100);\n",
+ " linewidth = 2,\n",
+ " axis = (; xlabel = \"Site index\", ylabel = \"Correlation function\", xscale = log10, yscale = log10)\n",
+ ")\n",
+ "hlines!(ax, abs2(expectation_value(ground_state, (0,) => a_op)); linestyle = :dash, color = :black)\n",
+ "fig"
]
},
{
@@ -357,31 +365,26 @@
" two_point_correlation[idx, :] .= real.(expectation_value(states[idx], (1, i) => a_op' ⊗ a_op) for i in 1:npoints)\n",
"end\n",
"\n",
- "p = plot(\n",
- " framestyle = :box, ylabel = \"Correlation function \" * L\"\\langle a_i^{\\dagger}a_j \\rangle\",\n",
- " xlabel = \"Distance \" * L\"|i-j|\", xscale = :log10, yscale = :log10,\n",
- " xticks = ([10, 100], [\"10\", \"100\"]),\n",
- " yticks = ([0.25, 0.5, 1.0], [\"0.25\", \"0.5\", \"1.0\"])\n",
+ "fig = Figure()\n",
+ "ax = Axis(\n",
+ " fig[1, 1];\n",
+ " xlabel = L\"\\text{Distance } |i-j|\",\n",
+ " ylabel = L\"\\text{Correlation function } \\langle a_i^{\\dagger}a_j \\rangle\",\n",
+ " xscale = log10, yscale = log10, xticks = [10, 100], yticks = [0.25, 0.5, 1.0]\n",
")\n",
+ "for (D, correlation) in zip(Ds, eachrow(two_point_correlation))\n",
+ " lines!(ax, 2:npoints, correlation[2:end]; label = \"D = $D\", linewidth = 2)\n",
+ "end\n",
+ "axislegend(ax)\n",
"\n",
- "plot!(\n",
- " p, 2:npoints, two_point_correlation[:, 2:end]',\n",
- " lab = \"D = \" .* string.(permutedims(Ds)), lw = 2\n",
+ "inset = Axis(\n",
+ " fig[1, 1]; width = Relative(0.25), height = Relative(0.25), halign = 0.27, valign = 0.32,\n",
+ " xlabel = \"Bond dimension\", ylabel = \"Correlation length\", xlabelsize = 10, ylabelsize = 10,\n",
+ " xscale = log10, yscale = log10, xticks = 20:10:50, yticks = [50, 100],\n",
+ " limits = ((15, 60), (20, 130))\n",
")\n",
- "\n",
- "scatter!(\n",
- " p, Ds, map(ψ -> correlation_length(ψ; sector = leftunit(ψ)), states),\n",
- " ylabel = \"Correlation length\", xlabel = \"Bond dimension\",\n",
- " xscale = :log10, yscale = :log10,\n",
- " inset = bbox(0.2, 0.51, 0.25, 0.25),\n",
- " subplot = 2,\n",
- " xticks = (20:10:50, string.(20:10:50)),\n",
- " yticks = ([50, 100], string.([50, 100])),\n",
- " xlabelfontsize = 8,\n",
- " ylabelfontsize = 8,\n",
- " ylims = [20, 130],\n",
- " xlims = [15, 60]\n",
- ")"
+ "scatter!(inset, Ds, map(ψ -> correlation_length(ψ; sector = leftunit(ψ)), states))\n",
+ "fig"
]
},
{
@@ -453,7 +456,13 @@
" )\n",
")\n",
"momentum_distribution = vcat(momentum_distribution...)'\n",
- "plot(ks, momentum_distribution, lab = \"D = \" .* string.(permutedims(Ds)), lw = 1.5, xlabel = \"Momentum k\", ylabel = L\"\\langle n_k \\rangle\", ylim = [0, 50])"
+ "fig = Figure()\n",
+ "ax = Axis(fig[1, 1]; xlabel = \"Momentum k\", ylabel = L\"\\langle n_k \\rangle\", limits = (nothing, (0, 50)))\n",
+ "for (D, n_k) in zip(Ds, eachcol(momentum_distribution))\n",
+ " lines!(ax, ks, n_k; label = \"D = $D\", linewidth = 1.5)\n",
+ "end\n",
+ "axislegend(ax)\n",
+ "fig"
]
},
{
@@ -516,7 +525,13 @@
" energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted))\n",
" end\n",
"\n",
- " return plot(phis, energies, lw = 2, xlabel = \"Phase twist per site\" * L\"(\\phi)\", ylabel = \"Ground state energy\", title = \"t = $t | μ = $mu | D = $D | cutoff = $cutoff\")\n",
+ " return lines(\n",
+ " phis, energies; linewidth = 2,\n",
+ " axis = (;\n",
+ " xlabel = L\"\\text{Phase twist per site } (\\phi)\", ylabel = \"Ground state energy\",\n",
+ " title = \"t = $t | μ = $mu | D = $D | cutoff = $cutoff\",\n",
+ " )\n",
+ " )\n",
"end\n",
"\n",
"superfluid_stiffness_profile(0.2, 0.3, 5, 4) # superfluid\n",
@@ -553,7 +568,12 @@
" order_parameters[i, j] = abs(expectation_value(state, 0 => a_op))\n",
"end\n",
"\n",
- "heatmap(ts, mus, order_parameters, xlabel = L\"t/U\", ylabel = L\"\\mu/U\", title = L\"\\langle \\hat{a}_i \\rangle\")"
+ "fig, ax, hm = heatmap(\n",
+ " ts, mus, permutedims(order_parameters);\n",
+ " axis = (; xlabel = L\"t/U\", ylabel = L\"\\mu/U\", title = L\"\\langle \\hat{a}_i \\rangle\")\n",
+ ")\n",
+ "Colorbar(fig[1, 2], hm)\n",
+ "fig"
]
},
{
@@ -576,7 +596,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Julia 1.13.1",
+ "display_name": "Julia 1.13.3",
"language": "julia",
"name": "julia-1.13"
},
@@ -584,7 +604,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
- "version": "1.13.1"
+ "version": "1.13.3"
}
},
"nbformat": 4,
diff --git a/docs/src/examples/groundstates/haldane-spt/figure-1.png b/docs/src/examples/groundstates/haldane-spt/figure-1.png
index 9afbbea19..c4a912a1a 100644
Binary files a/docs/src/examples/groundstates/haldane-spt/figure-1.png and b/docs/src/examples/groundstates/haldane-spt/figure-1.png differ
diff --git a/docs/src/examples/groundstates/haldane-spt/figure-2.png b/docs/src/examples/groundstates/haldane-spt/figure-2.png
index 97219add5..9b3d17b61 100644
Binary files a/docs/src/examples/groundstates/haldane-spt/figure-2.png and b/docs/src/examples/groundstates/haldane-spt/figure-2.png differ
diff --git a/docs/src/examples/groundstates/haldane-spt/figure-3.png b/docs/src/examples/groundstates/haldane-spt/figure-3.png
index 42c179c9f..ff8bb1dda 100644
Binary files a/docs/src/examples/groundstates/haldane-spt/figure-3.png and b/docs/src/examples/groundstates/haldane-spt/figure-3.png differ
diff --git a/docs/src/examples/groundstates/haldane-spt/index.md b/docs/src/examples/groundstates/haldane-spt/index.md
index 97d3cae9d..efdb61d61 100644
--- a/docs/src/examples/groundstates/haldane-spt/index.md
+++ b/docs/src/examples/groundstates/haldane-spt/index.md
@@ -30,7 +30,7 @@ Consequently, the Hamiltonian also commutes with all elements of $SU(2)$.
````julia
using TensorKit
using MPSKit
-using Plots
+using CairoMakie
casimir(s::SU2Irrep) = s.j * (s.j + 1)
@@ -86,7 +86,9 @@ V_wrong = SU2Space(0 => 8, 1 // 2 => 8, 1 => 3, 3 // 2 => 3)
ψ = InfiniteMPS(ℋ, V_wrong)
ψ, environments, info = find_groundstate(ψ, H, VUMPS(; maxiter = 10))
sectors = SU2Irrep[0, 1 // 2, 1, 3 // 2]
-transferplot(ψ; sectors, title = "Transfer matrix spectrum", legend = :outertop)
+fig, ax, _ = transfer_plot(ψ; sectors, axis = (; title = "Transfer matrix spectrum"))
+Legend(fig[1, 2], ax)
+fig
````
@@ -109,7 +111,7 @@ E_plus = expectation_value(ψ_plus, H)
````
````
--1.4014193313393006 - 2.5835641861230988e-17im
+-1.4014193313393015 + 4.171495806356685e-17im
````
````julia
@@ -120,27 +122,24 @@ E_minus = expectation_value(ψ_minus, H)
````
````
--1.4014839739630836 + 1.5092213828325467e-16im
+-1.401483973963084 - 2.1341651321826966e-17im
````
````julia
-transferp_plus = transferplot(
- ψ_plus;
- sectors = SU2Irrep[0, 1, 2], title = "ψ_plus", legend = :outertop
-)
-transferp_minus = transferplot(
- ψ_minus;
- sectors = SU2Irrep[0, 1, 2], title = "ψ_minus", legend = :outertop
-)
-plot(transferp_plus, transferp_minus; layout = (1, 2), size = (800, 400))
+fig = Figure(; size = (800, 400))
+ax, _ = transfer_plot(fig[1, 1], ψ_plus; sectors = SU2Irrep[0, 1, 2], axis = (; title = "ψ_plus"))
+transfer_plot(fig[1, 2], ψ_minus; sectors = SU2Irrep[0, 1, 2], axis = (; title = "ψ_minus"))
+Legend(fig[2, 1:2], ax; orientation = :horizontal)
+fig
````
````julia
-entanglementp_plus = entanglementplot(ψ_plus; title = "ψ_plus", legend = :outertop)
-entanglementp_minus = entanglementplot(ψ_minus; title = "ψ_minus", legend = :outertop)
-plot(entanglementp_plus, entanglementp_minus; layout = (1, 2), size = (800, 400))
+fig = Figure(; size = (800, 400))
+entanglement_plot(fig[1, 1], ψ_plus; axis = (; title = "ψ_plus"))
+entanglement_plot(fig[1, 2], ψ_minus; axis = (; title = "ψ_minus"))
+fig
````
@@ -166,8 +165,8 @@ println("S_plus = $S_plus")
````
````
-S_minus + log(2) = 1.5486227235411327
-S_plus = 1.5450323530595034
+S_minus + log(2) = 1.548622723544442
+S_plus = 1.5450323530533225
````
diff --git a/docs/src/examples/groundstates/haldane-spt/main.ipynb b/docs/src/examples/groundstates/haldane-spt/main.ipynb
index 3f2a918be..f1fe9d20e 100644
--- a/docs/src/examples/groundstates/haldane-spt/main.ipynb
+++ b/docs/src/examples/groundstates/haldane-spt/main.ipynb
@@ -34,7 +34,7 @@
"source": [
"using TensorKit\n",
"using MPSKit\n",
- "using Plots\n",
+ "using CairoMakie\n",
"\n",
"casimir(s::SU2Irrep) = s.j * (s.j + 1)\n",
"\n",
@@ -90,7 +90,9 @@
"ψ = InfiniteMPS(ℋ, V_wrong)\n",
"ψ, environments, info = find_groundstate(ψ, H, VUMPS(; maxiter = 10))\n",
"sectors = SU2Irrep[0, 1 // 2, 1, 3 // 2]\n",
- "transferplot(ψ; sectors, title = \"Transfer matrix spectrum\", legend = :outertop)"
+ "fig, ax, _ = transfer_plot(ψ; sectors, axis = (; title = \"Transfer matrix spectrum\"))\n",
+ "Legend(fig[1, 2], ax)\n",
+ "fig"
]
},
{
@@ -122,18 +124,15 @@
"ψ_minus = InfiniteMPS(ℋ, V_minus)\n",
"ψ_minus, = find_groundstate(ψ_minus, H, VUMPS(; maxiter = 100))\n",
"E_minus = expectation_value(ψ_minus, H)\n",
- "transferp_plus = transferplot(\n",
- " ψ_plus;\n",
- " sectors = SU2Irrep[0, 1, 2], title = \"ψ_plus\", legend = :outertop\n",
- ")\n",
- "transferp_minus = transferplot(\n",
- " ψ_minus;\n",
- " sectors = SU2Irrep[0, 1, 2], title = \"ψ_minus\", legend = :outertop\n",
- ")\n",
- "plot(transferp_plus, transferp_minus; layout = (1, 2), size = (800, 400))\n",
- "entanglementp_plus = entanglementplot(ψ_plus; title = \"ψ_plus\", legend = :outertop)\n",
- "entanglementp_minus = entanglementplot(ψ_minus; title = \"ψ_minus\", legend = :outertop)\n",
- "plot(entanglementp_plus, entanglementp_minus; layout = (1, 2), size = (800, 400))"
+ "fig = Figure(; size = (800, 400))\n",
+ "ax, _ = transfer_plot(fig[1, 1], ψ_plus; sectors = SU2Irrep[0, 1, 2], axis = (; title = \"ψ_plus\"))\n",
+ "transfer_plot(fig[1, 2], ψ_minus; sectors = SU2Irrep[0, 1, 2], axis = (; title = \"ψ_minus\"))\n",
+ "Legend(fig[2, 1:2], ax; orientation = :horizontal)\n",
+ "fig\n",
+ "fig = Figure(; size = (800, 400))\n",
+ "entanglement_plot(fig[1, 1], ψ_plus; axis = (; title = \"ψ_plus\"))\n",
+ "entanglement_plot(fig[1, 2], ψ_minus; axis = (; title = \"ψ_minus\"))\n",
+ "fig"
]
},
{
@@ -178,7 +177,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Julia 1.13.1",
+ "display_name": "Julia 1.13.3",
"language": "julia",
"name": "julia-1.13"
},
@@ -186,7 +185,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
- "version": "1.13.1"
+ "version": "1.13.3"
}
},
"nbformat": 4,
diff --git a/docs/src/examples/groundstates/hubbard/figure-1.png b/docs/src/examples/groundstates/hubbard/figure-1.png
index 1689abf2b..4f1172e41 100644
Binary files a/docs/src/examples/groundstates/hubbard/figure-1.png and b/docs/src/examples/groundstates/hubbard/figure-1.png differ
diff --git a/docs/src/examples/groundstates/hubbard/figure-2.png b/docs/src/examples/groundstates/hubbard/figure-2.png
index 49b5daff9..ba3147753 100644
Binary files a/docs/src/examples/groundstates/hubbard/figure-2.png and b/docs/src/examples/groundstates/hubbard/figure-2.png differ
diff --git a/docs/src/examples/groundstates/hubbard/figure-3.png b/docs/src/examples/groundstates/hubbard/figure-3.png
index db0c2adce..b092541a9 100644
Binary files a/docs/src/examples/groundstates/hubbard/figure-3.png and b/docs/src/examples/groundstates/hubbard/figure-3.png differ
diff --git a/docs/src/examples/groundstates/hubbard/index.md b/docs/src/examples/groundstates/hubbard/index.md
index 6b234a1de..d5254248f 100644
--- a/docs/src/examples/groundstates/hubbard/index.md
+++ b/docs/src/examples/groundstates/hubbard/index.md
@@ -40,7 +40,7 @@ using MPSKit
using .ExampleModels
using SpecialFunctions: besselj0, besselj1
using QuadGK: quadgk
-using Plots
+using CairoMakie
using Interpolations
using Optim
````
@@ -121,7 +121,7 @@ Groundstate energy:
````
┌ Info: Groundstate energy:
-│ * numerical: -2.1899960609745746
+│ * numerical: -2.189996060974573
└ * analytic: -2.190038374277775
````
@@ -156,7 +156,7 @@ Groundstate energy:
````
┌ Info: Groundstate energy:
-│ * numerical: -2.1900153475144615
+│ * numerical: -2.1900153475144597
└ * analytic: -2.190038374277775
````
@@ -230,14 +230,19 @@ I_holon = sortperm(P_holon_analytic)
P_holon_analytic = P_holon_analytic[I_holon]
E_holon_analytic = E_holon_analytic[I_holon];
-p = let p_excitations = plot(; xaxis = "momentum", yaxis = "energy")
- scatter!(p_excitations, momenta, real(E_spinon); label = "spinon")
- plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
+# one point per momentum and excitation branch
+branches(momenta, E) = (repeat(momenta, size(E, 2)), vec(real(E)))
- scatter!(p_excitations, momenta, real(E_holon); label = "holon")
- plot!(p_excitations, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
+let fig = Figure()
+ ax = Axis(fig[1, 1]; xlabel = "momentum", ylabel = "energy")
+ scatter!(ax, branches(momenta, E_spinon)...; label = "spinon")
+ lines!(ax, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
- p_excitations
+ scatter!(ax, branches(momenta, E_holon)...; label = "holon")
+ lines!(ax, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
+
+ Legend(fig[1, 2], ax)
+ fig
end
````
@@ -250,14 +255,16 @@ Here, we can fix this shift by realizing that our choice of shifting the grounds
````julia
momenta_shifted = rem2pi.(momenta .- π / 2, RoundNearest)
-p = let p_excitations = plot(; xaxis = "momentum", yaxis = "energy", xlims = (-π, π))
- scatter!(p_excitations, momenta_shifted, real(E_spinon); label = "spinon")
- plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
+let fig = Figure()
+ ax = Axis(fig[1, 1]; xlabel = "momentum", ylabel = "energy", limits = ((-π, π), nothing))
+ scatter!(ax, branches(momenta_shifted, E_spinon)...; label = "spinon")
+ lines!(ax, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
- scatter!(p_excitations, momenta_shifted, real(E_holon); label = "holon")
- plot!(p_excitations, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
+ scatter!(ax, branches(momenta_shifted, E_holon)...; label = "holon")
+ lines!(ax, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
- p_excitations
+ Legend(fig[1, 2], ax)
+ fig
end
````
@@ -308,22 +315,22 @@ E_scattering_max = map(momenta_shifted) do p
return e
end;
-p = let p_excitations = plot(;
- xaxis = "momentum", yaxis = "energy", xlims = (-π, π), ylims = (-0.1, 5)
- )
- scatter!(p_excitations, momenta_shifted, real(E_spinon); label = "spinon")
- plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
+let fig = Figure()
+ ax = Axis(fig[1, 1]; xlabel = "momentum", ylabel = "energy", limits = ((-π, π), (-0.1, 5)))
+ scatter!(ax, branches(momenta_shifted, E_spinon)...; label = "spinon")
+ lines!(ax, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
- scatter!(p_excitations, momenta_shifted, real(E_holon); label = "holon")
- plot!(p_excitations, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
+ scatter!(ax, branches(momenta_shifted, E_holon)...; label = "holon")
+ lines!(ax, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
I = sortperm(momenta_shifted)
- plot!(
- p_excitations, momenta_shifted[I], E_scattering_min[I]; label = "scattering states",
- fillrange = E_scattering_max[I], fillalpha = 0.3, fillstyle = :x
+ band!(
+ ax, momenta_shifted[I], E_scattering_min[I], E_scattering_max[I];
+ label = "scattering states", color = (:gray, 0.3)
)
- p_excitations
+ Legend(fig[1, 2], ax)
+ fig
end
````
diff --git a/docs/src/examples/groundstates/hubbard/main.ipynb b/docs/src/examples/groundstates/hubbard/main.ipynb
index 858a5f677..df7c961c8 100644
--- a/docs/src/examples/groundstates/hubbard/main.ipynb
+++ b/docs/src/examples/groundstates/hubbard/main.ipynb
@@ -154,7 +154,7 @@
"using .ExampleModels\n",
"using SpecialFunctions: besselj0, besselj1\n",
"using QuadGK: quadgk\n",
- "using Plots\n",
+ "using CairoMakie\n",
"using Interpolations\n",
"using Optim"
]
@@ -376,14 +376,19 @@
"P_holon_analytic = P_holon_analytic[I_holon]\n",
"E_holon_analytic = E_holon_analytic[I_holon];\n",
"\n",
- "p = let p_excitations = plot(; xaxis = \"momentum\", yaxis = \"energy\")\n",
- " scatter!(p_excitations, momenta, real(E_spinon); label = \"spinon\")\n",
- " plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label = \"spinon (analytic)\")\n",
+ "# one point per momentum and excitation branch\n",
+ "branches(momenta, E) = (repeat(momenta, size(E, 2)), vec(real(E)))\n",
"\n",
- " scatter!(p_excitations, momenta, real(E_holon); label = \"holon\")\n",
- " plot!(p_excitations, P_holon_analytic, E_holon_analytic; label = \"holon (analytic)\")\n",
+ "let fig = Figure()\n",
+ " ax = Axis(fig[1, 1]; xlabel = \"momentum\", ylabel = \"energy\")\n",
+ " scatter!(ax, branches(momenta, E_spinon)...; label = \"spinon\")\n",
+ " lines!(ax, P_spinon_analytic, E_spinon_analytic; label = \"spinon (analytic)\")\n",
"\n",
- " p_excitations\n",
+ " scatter!(ax, branches(momenta, E_holon)...; label = \"holon\")\n",
+ " lines!(ax, P_holon_analytic, E_holon_analytic; label = \"holon (analytic)\")\n",
+ "\n",
+ " Legend(fig[1, 2], ax)\n",
+ " fig\n",
"end"
]
},
@@ -404,14 +409,16 @@
"outputs": [],
"source": [
"momenta_shifted = rem2pi.(momenta .- π / 2, RoundNearest)\n",
- "p = let p_excitations = plot(; xaxis = \"momentum\", yaxis = \"energy\", xlims = (-π, π))\n",
- " scatter!(p_excitations, momenta_shifted, real(E_spinon); label = \"spinon\")\n",
- " plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label = \"spinon (analytic)\")\n",
+ "let fig = Figure()\n",
+ " ax = Axis(fig[1, 1]; xlabel = \"momentum\", ylabel = \"energy\", limits = ((-π, π), nothing))\n",
+ " scatter!(ax, branches(momenta_shifted, E_spinon)...; label = \"spinon\")\n",
+ " lines!(ax, P_spinon_analytic, E_spinon_analytic; label = \"spinon (analytic)\")\n",
"\n",
- " scatter!(p_excitations, momenta_shifted, real(E_holon); label = \"holon\")\n",
- " plot!(p_excitations, P_holon_analytic, E_holon_analytic; label = \"holon (analytic)\")\n",
+ " scatter!(ax, branches(momenta_shifted, E_holon)...; label = \"holon\")\n",
+ " lines!(ax, P_holon_analytic, E_holon_analytic; label = \"holon (analytic)\")\n",
"\n",
- " p_excitations\n",
+ " Legend(fig[1, 2], ax)\n",
+ " fig\n",
"end"
]
},
@@ -470,22 +477,22 @@
" return e\n",
"end;\n",
"\n",
- "p = let p_excitations = plot(;\n",
- " xaxis = \"momentum\", yaxis = \"energy\", xlims = (-π, π), ylims = (-0.1, 5)\n",
- " )\n",
- " scatter!(p_excitations, momenta_shifted, real(E_spinon); label = \"spinon\")\n",
- " plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label = \"spinon (analytic)\")\n",
+ "let fig = Figure()\n",
+ " ax = Axis(fig[1, 1]; xlabel = \"momentum\", ylabel = \"energy\", limits = ((-π, π), (-0.1, 5)))\n",
+ " scatter!(ax, branches(momenta_shifted, E_spinon)...; label = \"spinon\")\n",
+ " lines!(ax, P_spinon_analytic, E_spinon_analytic; label = \"spinon (analytic)\")\n",
"\n",
- " scatter!(p_excitations, momenta_shifted, real(E_holon); label = \"holon\")\n",
- " plot!(p_excitations, P_holon_analytic, E_holon_analytic; label = \"holon (analytic)\")\n",
+ " scatter!(ax, branches(momenta_shifted, E_holon)...; label = \"holon\")\n",
+ " lines!(ax, P_holon_analytic, E_holon_analytic; label = \"holon (analytic)\")\n",
"\n",
" I = sortperm(momenta_shifted)\n",
- " plot!(\n",
- " p_excitations, momenta_shifted[I], E_scattering_min[I]; label = \"scattering states\",\n",
- " fillrange = E_scattering_max[I], fillalpha = 0.3, fillstyle = :x\n",
+ " band!(\n",
+ " ax, momenta_shifted[I], E_scattering_min[I], E_scattering_max[I];\n",
+ " label = \"scattering states\", color = (:gray, 0.3)\n",
" )\n",
"\n",
- " p_excitations\n",
+ " Legend(fig[1, 2], ax)\n",
+ " fig\n",
"end"
]
},
@@ -501,7 +508,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Julia 1.13.1",
+ "display_name": "Julia 1.13.3",
"language": "julia",
"name": "julia-1.13"
},
@@ -509,7 +516,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
- "version": "1.13.1"
+ "version": "1.13.3"
}
},
"nbformat": 4,
diff --git a/docs/src/examples/groundstates/ising-cft/figure-1.png b/docs/src/examples/groundstates/ising-cft/figure-1.png
index ec004e5a0..671babf1f 100644
Binary files a/docs/src/examples/groundstates/ising-cft/figure-1.png and b/docs/src/examples/groundstates/ising-cft/figure-1.png differ
diff --git a/docs/src/examples/groundstates/ising-cft/figure-2.png b/docs/src/examples/groundstates/ising-cft/figure-2.png
index 53fb8d8f2..f221b7dc0 100644
Binary files a/docs/src/examples/groundstates/ising-cft/figure-2.png and b/docs/src/examples/groundstates/ising-cft/figure-2.png differ
diff --git a/docs/src/examples/groundstates/ising-cft/figure-3.png b/docs/src/examples/groundstates/ising-cft/figure-3.png
index 92dbcc3b6..878f7677a 100644
Binary files a/docs/src/examples/groundstates/ising-cft/figure-3.png and b/docs/src/examples/groundstates/ising-cft/figure-3.png differ
diff --git a/docs/src/examples/groundstates/ising-cft/index.md b/docs/src/examples/groundstates/ising-cft/index.md
index b7acb9d8d..8f701f4b1 100644
--- a/docs/src/examples/groundstates/ising-cft/index.md
+++ b/docs/src/examples/groundstates/ising-cft/index.md
@@ -12,7 +12,7 @@ This tutorial is meant to show the finite size CFT spectrum for the quantum Isin
We do this by first employing an exact diagonalization technique, and then extending the analysis to larger system sizes through the use of MPS techniques.
````julia
-using MPSKit, TensorKit, Plots, KrylovKit
+using MPSKit, TensorKit, CairoMakie, KrylovKit
using .ExampleModels
using LinearAlgebra: eigvals, diagm, Hermitian
````
@@ -51,10 +51,7 @@ Thus, optimizing the middle tensor is equivalent to optimizing a state in the en
````julia
energies, states = exact_diagonalization(H; num = 18, alg = Lanczos(; krylovdim = 200));
-plot(
- real.(energies);
- seriestype = :scatter, legend = false, ylabel = "energy", xlabel = "#eigenvalue"
-)
+scatter(real.(energies); axis = (; xlabel = "#eigenvalue", ylabel = "energy"))
````
@@ -119,24 +116,24 @@ append!(momenta, fix_degeneracies(states[17:18]))
````
18-element Vector{Float64}:
- -2.1278136533158578e-17
- 2.3050344087760185e-17
- -5.0324217423779645e-17
- 0.5235987755982991
- -0.5235987755982991
- -1.0471975511965976
- 1.0471975511965974
- -0.523598775598299
+ 1.4801423271942085e-18
+ -5.394357855450486e-18
+ -6.467597289392403e-18
0.5235987755982988
- 1.0471975511965979
- -1.0471975511965976
- -1.9939344511686805e-17
- -1.5707963267948966
- 1.5707963267948961
+ -0.523598775598299
+ 1.0471975511965983
+ -1.047197551196597
+ -0.5235987755982986
+ 0.523598775598299
+ -1.047197551196598
+ 1.0471975511965974
+ 1.413981352586968e-18
+ -1.570796326794897
+ 1.5707963267948963
-1.047197551196598
- 1.047197551196598
- -1.5707963267948968
- 1.5707963267948966
+ 1.0471975511965976
+ 1.5707963267948968
+ -1.5707963267948966
````
We can compute the scaling dimensions $Δ_n$ of the operators in the CFT from the energy gap of the corresponding excitations as $E_n - E_0 = \frac{2π v}{L} Δ_n$, where $v = 2$.
@@ -147,14 +144,12 @@ v = 2.0
Δ = real.(energies[1:18] .- energies[1]) ./ (2π * v / L)
S = momenta ./ (2π / L)
-p = plot(
- S, real.(Δ);
- seriestype = :scatter, xlabel = "conformal spin (S)", ylabel = "scaling dimension (Δ)",
- legend = false
+fig, ax, _ = scatter(
+ S, real.(Δ); axis = (; xlabel = "conformal spin (S)", ylabel = "scaling dimension (Δ)")
)
-vline!(p, -3:3; color = "gray", linestyle = :dash)
-hline!(p, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = "gray", linestyle = :dash)
-p
+vlines!(ax, -3:3; color = :gray, linestyle = :dash)
+hlines!(ax, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = :gray, linestyle = :dash)
+fig
````
@@ -193,14 +188,12 @@ v = 2.0
Δ_mps = real.(energies_mps[1:18] .- energies_mps[1]) ./ (2π * v / L_mps)
S_mps = momenta_mps ./ (2π / L_mps)
-p_mps = plot(
- S_mps, real.(Δ_mps);
- seriestype = :scatter, xlabel = "conformal spin (S)",
- ylabel = "scaling dimension (Δ)", legend = false
+fig, ax, _ = scatter(
+ S_mps, real.(Δ_mps); axis = (; xlabel = "conformal spin (S)", ylabel = "scaling dimension (Δ)")
)
-vline!(p_mps, -3:3; color = "gray", linestyle = :dash)
-hline!(p_mps, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = "gray", linestyle = :dash)
-p_mps
+vlines!(ax, -3:3; color = :gray, linestyle = :dash)
+hlines!(ax, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = :gray, linestyle = :dash)
+fig
````
diff --git a/docs/src/examples/groundstates/ising-cft/main.ipynb b/docs/src/examples/groundstates/ising-cft/main.ipynb
index 81e2702a7..e7f662b9b 100644
--- a/docs/src/examples/groundstates/ising-cft/main.ipynb
+++ b/docs/src/examples/groundstates/ising-cft/main.ipynb
@@ -16,7 +16,7 @@
"metadata": {},
"outputs": [],
"source": [
- "using MPSKit, TensorKit, Plots, KrylovKit\n",
+ "using MPSKit, TensorKit, CairoMakie, KrylovKit\n",
"# Model constructors for the MPSKit documentation and examples.\n",
"# Adapted from MPSKitModels.jl v0.4.7, src/models/{hamiltonians,transfermatrices}.jl:\n",
"# https://github.com/QuantumKitHub/MPSKitModels.jl/tree/v0.4.7/src/models\n",
@@ -160,10 +160,7 @@
"outputs": [],
"source": [
"energies, states = exact_diagonalization(H; num = 18, alg = Lanczos(; krylovdim = 200));\n",
- "plot(\n",
- " real.(energies);\n",
- " seriestype = :scatter, legend = false, ylabel = \"energy\", xlabel = \"#eigenvalue\"\n",
- ")"
+ "scatter(real.(energies); axis = (; xlabel = \"#eigenvalue\", ylabel = \"energy\"))"
]
},
{
@@ -264,14 +261,12 @@
"Δ = real.(energies[1:18] .- energies[1]) ./ (2π * v / L)\n",
"S = momenta ./ (2π / L)\n",
"\n",
- "p = plot(\n",
- " S, real.(Δ);\n",
- " seriestype = :scatter, xlabel = \"conformal spin (S)\", ylabel = \"scaling dimension (Δ)\",\n",
- " legend = false\n",
+ "fig, ax, _ = scatter(\n",
+ " S, real.(Δ); axis = (; xlabel = \"conformal spin (S)\", ylabel = \"scaling dimension (Δ)\")\n",
")\n",
- "vline!(p, -3:3; color = \"gray\", linestyle = :dash)\n",
- "hline!(p, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = \"gray\", linestyle = :dash)\n",
- "p"
+ "vlines!(ax, -3:3; color = :gray, linestyle = :dash)\n",
+ "hlines!(ax, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = :gray, linestyle = :dash)\n",
+ "fig"
]
},
{
@@ -328,14 +323,12 @@
"Δ_mps = real.(energies_mps[1:18] .- energies_mps[1]) ./ (2π * v / L_mps)\n",
"S_mps = momenta_mps ./ (2π / L_mps)\n",
"\n",
- "p_mps = plot(\n",
- " S_mps, real.(Δ_mps);\n",
- " seriestype = :scatter, xlabel = \"conformal spin (S)\",\n",
- " ylabel = \"scaling dimension (Δ)\", legend = false\n",
+ "fig, ax, _ = scatter(\n",
+ " S_mps, real.(Δ_mps); axis = (; xlabel = \"conformal spin (S)\", ylabel = \"scaling dimension (Δ)\")\n",
")\n",
- "vline!(p_mps, -3:3; color = \"gray\", linestyle = :dash)\n",
- "hline!(p_mps, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = \"gray\", linestyle = :dash)\n",
- "p_mps"
+ "vlines!(ax, -3:3; color = :gray, linestyle = :dash)\n",
+ "hlines!(ax, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = :gray, linestyle = :dash)\n",
+ "fig"
]
},
{
@@ -350,7 +343,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Julia 1.13.1",
+ "display_name": "Julia 1.13.3",
"language": "julia",
"name": "julia-1.13"
},
@@ -358,7 +351,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
- "version": "1.13.1"
+ "version": "1.13.3"
}
},
"nbformat": 4,
diff --git a/docs/src/examples/groundstates/xxz-heisenberg/figure-1.png b/docs/src/examples/groundstates/xxz-heisenberg/figure-1.png
index a69545d1b..61e38e9be 100644
Binary files a/docs/src/examples/groundstates/xxz-heisenberg/figure-1.png and b/docs/src/examples/groundstates/xxz-heisenberg/figure-1.png differ
diff --git a/docs/src/examples/groundstates/xxz-heisenberg/figure-2.png b/docs/src/examples/groundstates/xxz-heisenberg/figure-2.png
index c19958f0d..e5b2b3be3 100644
Binary files a/docs/src/examples/groundstates/xxz-heisenberg/figure-2.png and b/docs/src/examples/groundstates/xxz-heisenberg/figure-2.png differ
diff --git a/docs/src/examples/groundstates/xxz-heisenberg/index.md b/docs/src/examples/groundstates/xxz-heisenberg/index.md
index 00f23316e..8222e8961 100644
--- a/docs/src/examples/groundstates/xxz-heisenberg/index.md
+++ b/docs/src/examples/groundstates/xxz-heisenberg/index.md
@@ -12,7 +12,7 @@ In this file we will give step by step instructions on how to analyze the spin 1
The necessary packages to follow this tutorial are:
````julia
-using MPSKit, TensorKit, Plots
+using MPSKit, TensorKit, CairoMakie
using .ExampleModels
````
@@ -70,7 +70,7 @@ info
````
AlgorithmInfo:
converged = false after 200 iterations
- galerkin = 0.370633
+ galerkin = 0.37968
````
@@ -86,7 +86,7 @@ info
````
AlgorithmInfo:
converged = false after 20 iterations
- gradientnorm = 0.00747596
+ gradientnorm = 0.00594684
````
@@ -94,7 +94,7 @@ Convergence is quite slow and even fails after sufficiently many iterations.
To understand why, we can look at the transfer matrix spectrum.
````julia
-transferplot(groundstate, groundstate)
+transfer_plot(groundstate, groundstate)
````
@@ -137,7 +137,7 @@ groundstate, envs, info = find_groundstate(
````
````
-┌ Warning: VUMPS cancel 100: obj = -8.862417624752e-01 err = 4.2010390011e-06 time = 2.81 sec
+┌ Warning: VUMPS cancel 100: obj = -8.862417624752e-01 err = 4.2010394685e-06 time = 2.25 sec
└ @ MPSKit src/algorithms/groundstate/vumps.jl:96
````
@@ -149,7 +149,7 @@ The reason behind this becomes more obvious at higher bond dimensions:
groundstate, envs, info = find_groundstate(
state, H2, IDMRG2(; trunc = truncrank(50), maxiter = 20, tol = 1.0e-12, verbosity = 1)
);
-entanglementplot(groundstate)
+entanglement_plot(groundstate)
````
diff --git a/docs/src/examples/groundstates/xxz-heisenberg/main.ipynb b/docs/src/examples/groundstates/xxz-heisenberg/main.ipynb
index e2ca4974b..020584d6d 100644
--- a/docs/src/examples/groundstates/xxz-heisenberg/main.ipynb
+++ b/docs/src/examples/groundstates/xxz-heisenberg/main.ipynb
@@ -16,7 +16,7 @@
"metadata": {},
"outputs": [],
"source": [
- "using MPSKit, TensorKit, Plots\n",
+ "using MPSKit, TensorKit, CairoMakie\n",
"# Model constructors for the MPSKit documentation and examples.\n",
"# Adapted from MPSKitModels.jl v0.4.7, src/models/{hamiltonians,transfermatrices}.jl:\n",
"# https://github.com/QuantumKitHub/MPSKitModels.jl/tree/v0.4.7/src/models\n",
@@ -229,7 +229,7 @@
"metadata": {},
"outputs": [],
"source": [
- "transferplot(groundstate, groundstate)"
+ "transfer_plot(groundstate, groundstate)"
]
},
{
@@ -299,7 +299,7 @@
"groundstate, envs, info = find_groundstate(\n",
" state, H2, IDMRG2(; trunc = truncrank(50), maxiter = 20, tol = 1.0e-12, verbosity = 1)\n",
");\n",
- "entanglementplot(groundstate)"
+ "entanglement_plot(groundstate)"
]
},
{
@@ -389,7 +389,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Julia 1.13.1",
+ "display_name": "Julia 1.13.3",
"language": "julia",
"name": "julia-1.13"
},
@@ -397,7 +397,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
- "version": "1.13.1"
+ "version": "1.13.3"
}
},
"nbformat": 4,
diff --git a/docs/src/examples/statmech/hard-hexagon/figure-1.png b/docs/src/examples/statmech/hard-hexagon/figure-1.png
index 7174423ee..5a0f6b62d 100644
Binary files a/docs/src/examples/statmech/hard-hexagon/figure-1.png and b/docs/src/examples/statmech/hard-hexagon/figure-1.png differ
diff --git a/docs/src/examples/statmech/hard-hexagon/index.md b/docs/src/examples/statmech/hard-hexagon/index.md
index 66c3d3207..d82d445db 100644
--- a/docs/src/examples/statmech/hard-hexagon/index.md
+++ b/docs/src/examples/statmech/hard-hexagon/index.md
@@ -17,7 +17,7 @@ Because TensorKit supports anyonic symmetries, so does MPSKit.
To follow the tutorial you need the following packages.
````julia
-using MPSKit, TensorKit, Plots, Polynomials
+using MPSKit, TensorKit, CairoMakie, Polynomials
using .ExampleModels
````
@@ -60,7 +60,7 @@ println("F = $F\tS = $S\tξ = $ξ")
````
````
-F = 0.12340715308061596 S = 1.2807829607025263 ξ = 13.849682470052896
+F = 0.12340715308061621 S = 1.2807829607025245 ξ = 13.849682470052645
````
@@ -107,13 +107,17 @@ c = f.coeffs[2]
````
````
-0.8025235781819059
+0.8025235781819732
````
````julia
-p = plot(; xlabel = "logarithmic correlation length", ylabel = "entanglement entropy")
-p = plot(log.(ξs), Ss; seriestype = :scatter, label = nothing)
-plot!(p, ξ -> f(ξ) / 6; label = "fit")
+fig, ax, _ = scatter(
+ log.(ξs), Ss;
+ axis = (; xlabel = "logarithmic correlation length", ylabel = "entanglement entropy")
+)
+lines!(ax, range(extrema(log.(ξs))..., 100), ξ -> f(ξ) / 6; label = "fit")
+axislegend(ax; position = :lt)
+fig
````
diff --git a/docs/src/examples/statmech/hard-hexagon/main.ipynb b/docs/src/examples/statmech/hard-hexagon/main.ipynb
index 0791de51e..7d7b0dcb5 100644
--- a/docs/src/examples/statmech/hard-hexagon/main.ipynb
+++ b/docs/src/examples/statmech/hard-hexagon/main.ipynb
@@ -21,7 +21,7 @@
"metadata": {},
"outputs": [],
"source": [
- "using MPSKit, TensorKit, Plots, Polynomials\n",
+ "using MPSKit, TensorKit, CairoMakie, Polynomials\n",
"# Model constructors for the MPSKit documentation and examples.\n",
"# Adapted from MPSKitModels.jl v0.4.7, src/models/{hamiltonians,transfermatrices}.jl:\n",
"# https://github.com/QuantumKitHub/MPSKitModels.jl/tree/v0.4.7/src/models\n",
@@ -238,9 +238,13 @@
"\n",
"f = fit(log.(ξs), 6 * Ss, 1)\n",
"c = f.coeffs[2]\n",
- "p = plot(; xlabel = \"logarithmic correlation length\", ylabel = \"entanglement entropy\")\n",
- "p = plot(log.(ξs), Ss; seriestype = :scatter, label = nothing)\n",
- "plot!(p, ξ -> f(ξ) / 6; label = \"fit\")"
+ "fig, ax, _ = scatter(\n",
+ " log.(ξs), Ss;\n",
+ " axis = (; xlabel = \"logarithmic correlation length\", ylabel = \"entanglement entropy\")\n",
+ ")\n",
+ "lines!(ax, range(extrema(log.(ξs))..., 100), ξ -> f(ξ) / 6; label = \"fit\")\n",
+ "axislegend(ax; position = :lt)\n",
+ "fig"
]
},
{
@@ -255,7 +259,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Julia 1.13.1",
+ "display_name": "Julia 1.13.3",
"language": "julia",
"name": "julia-1.13"
},
@@ -263,7 +267,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
- "version": "1.13.1"
+ "version": "1.13.3"
}
},
"nbformat": 4,
diff --git a/examples/Cache.toml b/examples/Cache.toml
index da0152820..6f2cc32b2 100644
--- a/examples/Cache.toml
+++ b/examples/Cache.toml
@@ -1,16 +1,16 @@
[excitations]
-haldane = "6303c13abd414504a3a2a74f787ca409e8b4259222b875dc2790960d26eff603"
+haldane = "2e79dd49111097b0541bb6039e2c4e23f615b61d59611da6fd633df229a20b5f"
[statmech]
-hard-hexagon = "1532649b29e1d4584cbf823a35ef66ed7c5aea40f58021a0e0a59588b2857bee"
+hard-hexagon = "68269e430555ac77e0ed1d1c3b44e553ffa3bfb6a37112eb8c2aff350cd81e1c"
[dynamics]
-ising-dqpt = "ed3256a3aa9bf17a6e717e0758c560e85b5c151d0041fd6951d396f211d31121"
-xy-finiteT = "bb3866372b40c75b32f1f6ea540ed398fa0ca914ebd72de17688677d7ba05226"
+ising-dqpt = "cd6d0545172ada42735a7847d7da35be46bdf7a1b038366c7fa8046585269c92"
+xy-finiteT = "11e794f7e8700da83bbeb75da8c697a47957c9469b830496974c7df442d854ab"
[groundstates]
-bose-hubbard = "547847f63135707738efd75cf589b666a15839a7758ad3e2147a84ca5ddb7695"
-haldane-spt = "d61d0ef25c51e5fcbf00499b51295260fbbf818714a4eeede2f20a0af21d1a86"
-xxz-heisenberg = "a1f4d49cc4518c4773a0422a23a1f5998bd2d9b5016dda7e64b19fbc18c2c3c0"
-ising-cft = "78fd4f151783014b824a68fcd08b3865b5e268a619924c48ba11829038b81374"
-hubbard = "5b0dae1f6c0de939b7cbfdeaf16fb7eb945b4a8574acbbb375d4e4bb9aa1586f"
+bose-hubbard = "5411a7d968e9fb0071528f81946536046d12820f39dc358aaae1e1a3be9580bf"
+haldane-spt = "4ad3562ff82691ac31e0851c30bf5c7975d7bbcb06494f93ba168084ef912d42"
+xxz-heisenberg = "aad20bd24766dce3ee103e815723d1d5c15951a6f69557660525308a94b5e229"
+ising-cft = "9715c2348b7bf056bf1742d4f0940faa3eb3d30857d439e027993fb8fbb98e0d"
+hubbard = "66faf9d13d42abae30248d288f1e628c5cc5664e0ee69c6faf34b6b06d5ef75c"
diff --git a/examples/Project.toml b/examples/Project.toml
index 827d34e2e..262cac131 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -1,13 +1,12 @@
[deps]
BenchmarkFreeFermions = "5b68a00d-ca4d-4b58-ab0c-dc082ade624e"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
-LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
Optim = "429524aa-4258-5aef-a3af-852621145aeb"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
diff --git a/examples/dynamics/xy-finiteT/main.jl b/examples/dynamics/xy-finiteT/main.jl
index ea457a35f..6cf0269c1 100644
--- a/examples/dynamics/xy-finiteT/main.jl
+++ b/examples/dynamics/xy-finiteT/main.jl
@@ -7,7 +7,7 @@ using .ExampleModels
using MPSKit: infinite_temperature_density_matrix
using QuadGK: quadgk
using SpecialFunctions: ellipe
-using Plots
+using CairoMakie
using LinearAlgebra
using BenchmarkFreeFermions
@@ -180,23 +180,39 @@ Z_taylor = map(Iterators.product(βs, expansion_orders)) do (β, expansion_order
end
F_taylor = -(1 ./ βs) .* Z_taylor
-p_taylor = let
- labels = reshape(
- map(expansion_orders) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, Z_analytic; label = "analytic", title = "Partition function",
- xlabel = "β", ylabel = "Z(β)"
- )
- plot!(p1, βs, Z_taylor; label = labels)
- p2 = plot(
- βs, F_analytic; label = "analytic", title = "Free energy",
- xlabel = "β", ylabel = "F(β)"
- )
- plot!(p2, βs, F_taylor; label = labels)
- plot(p1, p2)
+## Makie skips NaN, but picks wrong limits for infinite or (on a log scale) zero values
+nan_unless(f, ys) = [f(y) ? y : NaN for y in ys]
+
+## Plot the errors of several methods side by side, given as `label => (βs, ΔZ, ΔF)`.
+function plot_errors(errors; position = :lt, yscale = identity)
+ fig = Figure(; size = (800, 400))
+ axZ = Axis(fig[1, 1]; title = "Partition function error", xlabel = "β", ylabel = "ΔZ(β)", yscale)
+ axF = Axis(fig[1, 2]; title = "Free energy error", xlabel = "β", ylabel = "ΔF(β)", yscale)
+ valid = yscale === identity ? isfinite : y -> isfinite(y) && y > 0
+ for (label, (β, ΔZ, ΔF)) in errors
+ lines!(axZ, β, nan_unless(valid, ΔZ); label)
+ lines!(axF, β, nan_unless(valid, ΔF); label)
+ end
+ axislegend(axZ; position)
+ return fig
+end
+
+taylor_errors(Z, F, orders) = [
+ "Taylor N=$N" => (βs, abs.(Z[:, i] .- Z_analytic), abs.(F[:, i] .- F_analytic))
+ for (i, N) in enumerate(orders)
+]
+
+p_taylor = let fig = Figure(; size = (800, 400))
+ axZ = Axis(fig[1, 1]; title = "Partition function", xlabel = "β", ylabel = "Z(β)")
+ axF = Axis(fig[1, 2]; title = "Free energy", xlabel = "β", ylabel = "F(β)")
+ lines!(axZ, βs, Z_analytic; label = "analytic")
+ lines!(axF, βs, nan_unless(isfinite, F_analytic); label = "analytic")
+ for (i, N) in enumerate(expansion_orders)
+ lines!(axZ, βs, Z_taylor[:, i]; label = "Taylor N=$N")
+ lines!(axF, βs, nan_unless(isfinite, F_taylor[:, i]); label = "Taylor N=$N")
+ end
+ axislegend(axZ; position = :lt)
+ fig
end
md"""
@@ -234,23 +250,7 @@ expansion_orders = 2:3
Z_taylor = Z_taylor[:, 2:end]
F_taylor = F_taylor[:, 2:end]
-p_taylor_diff = let
- labels = reshape(
- map(expansion_orders) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, abs.(Z_taylor .- Z_analytic);
- label = labels, title = "Partition function error",
- xlabel = "β", ylabel = "ΔZ(β)", legend = :topleft
- )
- p2 = plot(
- βs, abs.(F_taylor .- F_analytic); label = labels,
- xlabel = "β", ylabel = "ΔF(β)", title = "Free energy error", legend = :topleft
- )
- plot(p1, p2)
-end
+p_taylor_diff = plot_errors(taylor_errors(Z_taylor, F_taylor, expansion_orders))
md"""
We can now clearly see that, somewhat unsurprisingly, the error increases the larger $β$ becomes.
@@ -284,23 +284,7 @@ Z_taylor2 = map(Iterators.product(βs, expansion_orders)) do (β, expansion_orde
end
F_taylor2 = -(1 ./ βs) .* Z_taylor2
-p_taylor2_diff = let
- labels = reshape(
- map(expansion_orders[2:end]) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, abs.(Z_taylor2 .- Z_analytic);
- label = labels, title = "Partition function error",
- xlabel = "β", ylabel = "ΔZ(β)", legend = :topleft
- )
- p2 = plot(
- βs, abs.(F_taylor2 .- F_analytic); label = labels,
- xlabel = "β", ylabel = "ΔF(β)", title = "Free energy error", legend = :topleft
- )
- plot(p1, p2)
-end
+p_taylor2_diff = plot_errors(taylor_errors(Z_taylor2, F_taylor2, expansion_orders))
md"""
### MPO multiplication approach (linear)
@@ -350,32 +334,12 @@ for i in 3:length(βs)
end
F_mpo_mul = -(1 ./ βs) .* Z_mpo_mul
-p_mpo_mul_diff = let
- labels = reshape(
- map(expansion_orders) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, abs.(Z_taylor2 .- Z_analytic);
- label = labels, title = "Partition function error",
- xlabel = "β", ylabel = "ΔZ(β)", legend = :bottomright, yscale = :log10
- )
- plot!(
- p1, βs, abs.(Z_mpo_mul .- Z_analytic);
- label = "MPO multiplication"
- )
- p2 = plot(
- βs, abs.(F_taylor2 .- F_analytic); label = labels,
- xlabel = "β", ylabel = "ΔF(β)", title = "Free energy error", legend = nothing,
- yscale = :log10
- )
- plot!(
- p2, βs, abs.(F_mpo_mul .- F_analytic);
- label = "MPO multiplication"
- )
- plot(p1, p2)
-end
+p_mpo_mul_diff = plot_errors(
+ [
+ taylor_errors(Z_taylor2, F_taylor2, expansion_orders);
+ "MPO multiplication" => (βs, abs.(Z_mpo_mul .- Z_analytic), abs.(F_mpo_mul .- F_analytic))
+ ]; position = :rb, yscale = log10
+)
md"""
This approach clearly improves the accuracy of the results, indicating that we can indeed compute partition functions at larger $β$ values.
@@ -434,28 +398,13 @@ for i in 2:length(βs_exp)
end
F_mpo_mul_exp = -(1 ./ βs_exp) .* Z_mpo_mul_exp
-p_mpo_mul_exp_diff = let
- labels = reshape(
- map(expansion_orders[2:end]) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, abs.(Z_taylor2 .- Z_analytic);
- label = labels, title = "Partition function error", xlabel = "β", ylabel = "ΔZ(β)",
- legend = :bottomright, yscale = :log10
- )
- plot!(p1, βs, abs.(Z_mpo_mul .- Z_analytic); label = "MPO multiplication")
- plot!(p1, βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp); label = "MPO multiplication exp")
-
- p2 = plot(
- βs, abs.(F_taylor2 .- F_analytic); label = labels, xlabel = "β", ylabel = "ΔF(β)",
- title = "Free energy error", legend = nothing, yscale = :log10
- )
- plot!(p2, βs, abs.(F_mpo_mul .- F_analytic); label = "MPO multiplication")
- plot!(p2, βs_exp, abs.(F_mpo_mul_exp .- F_analytic_exp); label = "MPO multiplication exp")
- plot(p1, p2)
-end
+p_mpo_mul_exp_diff = plot_errors(
+ [
+ taylor_errors(Z_taylor2, F_taylor2, expansion_orders);
+ "MPO multiplication" => (βs, abs.(Z_mpo_mul .- Z_analytic), abs.(F_mpo_mul .- F_analytic))
+ "MPO multiplication exp" => (βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp), abs.(F_mpo_mul_exp .- F_analytic_exp))
+ ]; position = :rb, yscale = log10
+)
md"""
Clearly, the exponential approach allows us to reach larger $β$ values much quicker, but there is again a trade-off.
@@ -501,31 +450,14 @@ for i in 2:length(βs)
end
F_tdvp = -(1 ./ βs) .* Z_tdvp
-p_mpo_mul_diff = let
- labels = reshape(
- map(expansion_orders) do N
- return "Taylor N=$N"
- end, 1, :
- )
- p1 = plot(
- βs, abs.(Z_taylor2 .- Z_analytic); label = labels,
- title = "Partition function error", xlabel = "β", ylabel = "ΔZ(β)",
- legend = :bottomright, yscale = :log10
- )
- plot!(p1, βs, abs.(Z_mpo_mul .- Z_analytic); label = "MPO multiplication")
- plot!(p1, βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp); label = "MPO multiplication exp")
- plot!(p1, βs, abs.(Z_tdvp .- Z_analytic); label = "TDVP")
-
- p2 = plot(
- βs, abs.(F_taylor2 .- F_analytic); label = labels, xlabel = "β", ylabel = "ΔF(β)",
- title = "Free energy error", legend = nothing, yscale = :log10
- )
- plot!(p2, βs, abs.(F_mpo_mul .- F_analytic); label = "MPO multiplication")
- plot!(p2, βs_exp, abs.(F_mpo_mul_exp .- F_analytic_exp); label = "MPO multiplication exp")
- plot!(p2, βs, abs.(F_tdvp .- F_analytic); label = "TDVP")
-
- plot(p1, p2)
-end
+p_mpo_mul_diff = plot_errors(
+ [
+ taylor_errors(Z_taylor2, F_taylor2, expansion_orders);
+ "MPO multiplication" => (βs, abs.(Z_mpo_mul .- Z_analytic), abs.(F_mpo_mul .- F_analytic))
+ "MPO multiplication exp" => (βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp), abs.(F_mpo_mul_exp .- F_analytic_exp))
+ "TDVP" => (βs, abs.(Z_tdvp .- Z_analytic), abs.(F_tdvp .- F_analytic))
+ ]; position = :rb, yscale = log10
+)
md"""
!!! note
diff --git a/examples/excitations/haldane/main.jl b/examples/excitations/haldane/main.jl
index 49c2a6486..4023ca35e 100644
--- a/examples/excitations/haldane/main.jl
+++ b/examples/excitations/haldane/main.jl
@@ -5,7 +5,7 @@ In this tutorial we will calculate the Haldane gap (the energy gap in the ``S =
To follow the tutorial you need the following packages:
"""
-using MPSKit, TensorKit, Plots, Polynomials
+using MPSKit, TensorKit, CairoMakie, Polynomials
include(joinpath(@__DIR__, "..", "..", "..", "docs", "src", "assets", "models.jl")) #hide
using .ExampleModels
using TensorKitTensors.SpinOperators: S_exchange
@@ -54,16 +54,19 @@ md"""
We can go even further and doublecheck the claim that ``S = 1`` is an edge excitation, by plotting the energy density.
"""
-p_density = plot(; xaxis = "position", yaxis = "energy density")
excited_1 = convert(FiniteMPS, st_1[1])
excited_2 = convert(FiniteMPS, st_2[1])
SS = -S_exchange(ComplexF64, SU2Irrep; spin = 1)
e₀ = [real(expectation_value(ψ, (i, i + 1) => SS)) for i in 1:(L - 1)]
e₁ = [real(expectation_value(excited_1, (i, i + 1) => SS)) for i in 1:(L - 1)]
e₂ = [real(expectation_value(excited_2, (i, i + 1) => SS)) for i in 1:(L - 1)]
-plot!(p_density, e₀; label = "S = 0")
-plot!(p_density, e₁; label = "S = 1")
-plot!(p_density, e₂; label = "S = 2")
+fig = Figure()
+ax = Axis(fig[1, 1]; xlabel = "position", ylabel = "energy density")
+lines!(ax, e₀; label = "S = 0")
+lines!(ax, e₁; label = "S = 1")
+lines!(ax, e₂; label = "S = 2")
+axislegend(ax)
+fig
md"""
Finally, we can obtain a value for the Haldane gap by extrapolating our results for different system sizes.
@@ -85,9 +88,12 @@ f = fit(Ls .^ (-2), ΔEs, 1)
#+
-p_size_extrapolation = plot(; xaxis = "L^(-2)", yaxis = "ΔE", xlims = (0, 0.015))
-plot!(p_size_extrapolation, Ls .^ (-2), ΔEs; seriestype = :scatter, label = "numerical")
-plot!(p_size_extrapolation, x -> f(x); label = "fit")
+fig = Figure()
+ax = Axis(fig[1, 1]; xlabel = "L^(-2)", ylabel = "ΔE", limits = ((0, 0.015), nothing))
+scatter!(ax, Ls .^ (-2), ΔEs; label = "numerical")
+lines!(ax, range(0, 0.015, 100), x -> f(x); label = "fit")
+axislegend(ax; position = :lt)
+fig
md"""
## Thermodynamic limit
@@ -114,4 +120,4 @@ println("minimum @k = $(kspace[idx]):\t ΔE = $(ΔE)")
#+
-plot(kspace, real.(Es); xaxis = "momentum", yaxis = "ΔE", label = "S = 1")
+series(kspace, permutedims(real.(Es)); axis = (; xlabel = "momentum", ylabel = "ΔE", title = "S = 1"))
diff --git a/examples/figure_externalization.jl b/examples/figure_externalization.jl
index 54db91a87..5cfc09459 100644
--- a/examples/figure_externalization.jl
+++ b/examples/figure_externalization.jl
@@ -17,7 +17,7 @@ using Base64
# A `@raw html` fenced block whose body is a single `![]()
` tag.
const _RAW_IMG_BLOCK = r"```@raw html[ \t]*\r?\n[ \t]*(]*>)[ \t]*\r?\n```"
# A base64 image data URI inside such a tag: capture the MIME subtype and the payload.
-const _DATA_URI = r"data:image/([A-Za-z0-9.+-]+);base64,([A-Za-z0-9+/=]+)"
+const _DATA_URI = r"data:image/([A-Za-z0-9.+-]+);base64,\s*([A-Za-z0-9+/=]+)"
_ext_for(fmt) = fmt == "svg+xml" ? "svg" : fmt
diff --git a/examples/groundstates/bose-hubbard/main.jl b/examples/groundstates/bose-hubbard/main.jl
index a06049764..e47405007 100644
--- a/examples/groundstates/bose-hubbard/main.jl
+++ b/examples/groundstates/bose-hubbard/main.jl
@@ -3,14 +3,13 @@ using TensorKitTensors.BosonOperators: b_min, b_plus_b_min, b_min_b_plus, b_num
using MPSKit, TensorKit
include(joinpath(@__DIR__, "..", "..", "..", "docs", "src", "assets", "models.jl")) #hide
using .ExampleModels
-using Plots, LaTeXStrings
+using CairoMakie
#src # for reproducibility:
#src using Random
#src Random.seed!(123)
-theme(:wong)
-default(fontfamily = "Computer Modern", label = nothing, dpi = 100, framestyle = :box)
+set_theme!(theme_latexfonts())
md"""
# 1D Bose-Hubbard model
@@ -131,8 +130,12 @@ In case of the Hamiltonian, it is not necessary to specify the indices as it spa
We can now plot the correlation function $⟨â^{†}_i â_j⟩$.
"""
-plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw = 2, xlabel = "Site index", ylabel = "Correlation function", yscale = :log10)
-hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black)
+fig, ax, _ = lines(
+ map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50);
+ linewidth = 2, axis = (; xlabel = "Site index", ylabel = "Correlation function", yscale = log10)
+)
+hlines!(ax, abs2(expectation_value(ground_state, (0,) => a_op)); linestyle = :dash, color = :black)
+fig
md"""
We see that the correlations drop off exponentially, indicating the existence of a gapped Mott insulating phase.
@@ -141,8 +144,13 @@ Let us now shift our parameters to probe other phases.
ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500)
-plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw = 2, xlabel = "Site index", ylabel = "Correlation function", yscale = :log10, xscale = :log10)
-hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black)
+fig, ax, _ = lines(
+ map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100);
+ linewidth = 2,
+ axis = (; xlabel = "Site index", ylabel = "Correlation function", xscale = log10, yscale = log10)
+)
+hlines!(ax, abs2(expectation_value(ground_state, (0,) => a_op)); linestyle = :dash, color = :black)
+fig
md"""
In this case, the correlation function drops off algebraically and eventually saturates as $\lim_{i → ∞}⟨â_i^{†} â_j⟩ ≈ ⟨â_i^{†}⟩ ⟨â_j⟩ = |⟨a_i⟩|^2 ≠ 0$.
@@ -171,31 +179,26 @@ Threads.@threads for idx in eachindex(Ds)
two_point_correlation[idx, :] .= real.(expectation_value(states[idx], (1, i) => a_op' ⊗ a_op) for i in 1:npoints)
end
-p = plot(
- framestyle = :box, ylabel = "Correlation function " * L"\langle a_i^{\dagger}a_j \rangle",
- xlabel = "Distance " * L"|i-j|", xscale = :log10, yscale = :log10,
- xticks = ([10, 100], ["10", "100"]),
- yticks = ([0.25, 0.5, 1.0], ["0.25", "0.5", "1.0"])
-)
-
-plot!(
- p, 2:npoints, two_point_correlation[:, 2:end]',
- lab = "D = " .* string.(permutedims(Ds)), lw = 2
+fig = Figure()
+ax = Axis(
+ fig[1, 1];
+ xlabel = L"\text{Distance } |i-j|",
+ ylabel = L"\text{Correlation function } \langle a_i^{\dagger}a_j \rangle",
+ xscale = log10, yscale = log10, xticks = [10, 100], yticks = [0.25, 0.5, 1.0]
)
+for (D, correlation) in zip(Ds, eachrow(two_point_correlation))
+ lines!(ax, 2:npoints, correlation[2:end]; label = "D = $D", linewidth = 2)
+end
+axislegend(ax)
-scatter!(
- 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),
- subplot = 2,
- xticks = (20:10:50, string.(20:10:50)),
- yticks = ([50, 100], string.([50, 100])),
- xlabelfontsize = 8,
- ylabelfontsize = 8,
- ylims = [20, 130],
- xlims = [15, 60]
+inset = Axis(
+ fig[1, 1]; width = Relative(0.25), height = Relative(0.25), halign = 0.27, valign = 0.32,
+ xlabel = "Bond dimension", ylabel = "Correlation length", xlabelsize = 10, ylabelsize = 10,
+ xscale = log10, yscale = log10, xticks = 20:10:50, yticks = [50, 100],
+ limits = ((15, 60), (20, 130))
)
+scatter!(inset, Ds, map(ψ -> correlation_length(ψ; sector = leftunit(ψ)), states))
+fig
md"""
This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, forming a Bose-Einstein condensate which introduces erroneous long-distance behaviour of correlation functions.
@@ -247,7 +250,13 @@ momentum_distribution = map(
)
)
momentum_distribution = vcat(momentum_distribution...)'
-plot(ks, momentum_distribution, lab = "D = " .* string.(permutedims(Ds)), lw = 1.5, xlabel = "Momentum k", ylabel = L"\langle n_k \rangle", ylim = [0, 50])
+fig = Figure()
+ax = Axis(fig[1, 1]; xlabel = "Momentum k", ylabel = L"\langle n_k \rangle", limits = (nothing, (0, 50)))
+for (D, n_k) in zip(Ds, eachcol(momentum_distribution))
+ lines!(ax, ks, n_k; label = "D = $D", linewidth = 1.5)
+end
+axislegend(ax)
+fig
md"""
We see that the density seems to peak around $k=0$, this time seemingly becoming more prominent as $D → ∞$ which seems to suggest again that there is a condensate.
@@ -300,7 +309,13 @@ function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ = 1.0e-4, npoints = 1
energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted))
end
- return plot(phis, energies, lw = 2, xlabel = "Phase twist per site" * L"(\phi)", ylabel = "Ground state energy", title = "t = $t | μ = $mu | D = $D | cutoff = $cutoff")
+ return lines(
+ phis, energies; linewidth = 2,
+ axis = (;
+ xlabel = L"\text{Phase twist per site } (\phi)", ylabel = "Ground state energy",
+ title = "t = $t | μ = $mu | D = $D | cutoff = $cutoff",
+ )
+ )
end
superfluid_stiffness_profile(0.2, 0.3, 5, 4) # superfluid
@@ -327,7 +342,12 @@ Threads.@threads for (i, j) in collect(Iterators.product(eachindex(mus), eachind
order_parameters[i, j] = abs(expectation_value(state, 0 => a_op))
end
-heatmap(ts, mus, order_parameters, xlabel = L"t/U", ylabel = L"\mu/U", title = L"\langle \hat{a}_i \rangle")
+fig, ax, hm = heatmap(
+ ts, mus, permutedims(order_parameters);
+ axis = (; xlabel = L"t/U", ylabel = L"\mu/U", title = L"\langle \hat{a}_i \rangle")
+)
+Colorbar(fig[1, 2], hm)
+fig
md"""
Although the bond dimension here is quite low, we already see the deformation of the Mott insulator lobes to give way to the well known BKT transition that happens at commensurate density.
diff --git a/examples/groundstates/haldane-spt/main.jl b/examples/groundstates/haldane-spt/main.jl
index d566a35f5..d359baceb 100644
--- a/examples/groundstates/haldane-spt/main.jl
+++ b/examples/groundstates/haldane-spt/main.jl
@@ -23,7 +23,7 @@ Consequently, the Hamiltonian also commutes with all elements of $SU(2)$.
using TensorKit
using MPSKit
-using Plots
+using CairoMakie
casimir(s::SU2Irrep) = s.j * (s.j + 1)
@@ -69,7 +69,9 @@ V_wrong = SU2Space(0 => 8, 1 // 2 => 8, 1 => 3, 3 // 2 => 3)
ψ = InfiniteMPS(ℋ, V_wrong)
ψ, environments, info = find_groundstate(ψ, H, VUMPS(; maxiter = 10))
sectors = SU2Irrep[0, 1 // 2, 1, 3 // 2]
-transferplot(ψ; sectors, title = "Transfer matrix spectrum", legend = :outertop)
+fig, ax, _ = transfer_plot(ψ; sectors, axis = (; title = "Transfer matrix spectrum"))
+Legend(fig[1, 2], ax)
+fig
md"""
Nevertheless, using the symmetry, this can be remedied rather easily, by imposing the ground state to belong to a single class, and comparing the results.
@@ -97,21 +99,18 @@ E_minus = expectation_value(ψ_minus, H)
#+
-transferp_plus = transferplot(
- ψ_plus;
- sectors = SU2Irrep[0, 1, 2], title = "ψ_plus", legend = :outertop
-)
-transferp_minus = transferplot(
- ψ_minus;
- sectors = SU2Irrep[0, 1, 2], title = "ψ_minus", legend = :outertop
-)
-plot(transferp_plus, transferp_minus; layout = (1, 2), size = (800, 400))
+fig = Figure(; size = (800, 400))
+ax, _ = transfer_plot(fig[1, 1], ψ_plus; sectors = SU2Irrep[0, 1, 2], axis = (; title = "ψ_plus"))
+transfer_plot(fig[1, 2], ψ_minus; sectors = SU2Irrep[0, 1, 2], axis = (; title = "ψ_minus"))
+Legend(fig[2, 1:2], ax; orientation = :horizontal)
+fig
#+
-entanglementp_plus = entanglementplot(ψ_plus; title = "ψ_plus", legend = :outertop)
-entanglementp_minus = entanglementplot(ψ_minus; title = "ψ_minus", legend = :outertop)
-plot(entanglementp_plus, entanglementp_minus; layout = (1, 2), size = (800, 400))
+fig = Figure(; size = (800, 400))
+entanglement_plot(fig[1, 1], ψ_plus; axis = (; title = "ψ_plus"))
+entanglement_plot(fig[1, 2], ψ_minus; axis = (; title = "ψ_minus"))
+fig
md"""
diff --git a/examples/groundstates/hubbard/main.jl b/examples/groundstates/hubbard/main.jl
index ef0f9f398..885a4201e 100644
--- a/examples/groundstates/hubbard/main.jl
+++ b/examples/groundstates/hubbard/main.jl
@@ -32,7 +32,7 @@ include(joinpath(@__DIR__, "..", "..", "..", "docs", "src", "assets", "models.jl
using .ExampleModels
using SpecialFunctions: besselj0, besselj1
using QuadGK: quadgk
-using Plots
+using CairoMakie
using Interpolations
using Optim
@@ -202,14 +202,19 @@ I_holon = sortperm(P_holon_analytic)
P_holon_analytic = P_holon_analytic[I_holon]
E_holon_analytic = E_holon_analytic[I_holon];
-p = let p_excitations = plot(; xaxis = "momentum", yaxis = "energy")
- scatter!(p_excitations, momenta, real(E_spinon); label = "spinon")
- plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
+## one point per momentum and excitation branch
+branches(momenta, E) = (repeat(momenta, size(E, 2)), vec(real(E)))
- scatter!(p_excitations, momenta, real(E_holon); label = "holon")
- plot!(p_excitations, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
+let fig = Figure()
+ ax = Axis(fig[1, 1]; xlabel = "momentum", ylabel = "energy")
+ scatter!(ax, branches(momenta, E_spinon)...; label = "spinon")
+ lines!(ax, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
- p_excitations
+ scatter!(ax, branches(momenta, E_holon)...; label = "holon")
+ lines!(ax, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
+
+ Legend(fig[1, 2], ax)
+ fig
end
md"""
@@ -220,14 +225,16 @@ Here, we can fix this shift by realizing that our choice of shifting the grounds
"""
momenta_shifted = rem2pi.(momenta .- π / 2, RoundNearest)
-p = let p_excitations = plot(; xaxis = "momentum", yaxis = "energy", xlims = (-π, π))
- scatter!(p_excitations, momenta_shifted, real(E_spinon); label = "spinon")
- plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
+let fig = Figure()
+ ax = Axis(fig[1, 1]; xlabel = "momentum", ylabel = "energy", limits = ((-π, π), nothing))
+ scatter!(ax, branches(momenta_shifted, E_spinon)...; label = "spinon")
+ lines!(ax, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
- scatter!(p_excitations, momenta_shifted, real(E_holon); label = "holon")
- plot!(p_excitations, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
+ scatter!(ax, branches(momenta_shifted, E_holon)...; label = "holon")
+ lines!(ax, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
- p_excitations
+ Legend(fig[1, 2], ax)
+ fig
end
md"""
@@ -276,20 +283,20 @@ E_scattering_max = map(momenta_shifted) do p
return e
end;
-p = let p_excitations = plot(;
- xaxis = "momentum", yaxis = "energy", xlims = (-π, π), ylims = (-0.1, 5)
- )
- scatter!(p_excitations, momenta_shifted, real(E_spinon); label = "spinon")
- plot!(p_excitations, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
+let fig = Figure()
+ ax = Axis(fig[1, 1]; xlabel = "momentum", ylabel = "energy", limits = ((-π, π), (-0.1, 5)))
+ scatter!(ax, branches(momenta_shifted, E_spinon)...; label = "spinon")
+ lines!(ax, P_spinon_analytic, E_spinon_analytic; label = "spinon (analytic)")
- scatter!(p_excitations, momenta_shifted, real(E_holon); label = "holon")
- plot!(p_excitations, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
+ scatter!(ax, branches(momenta_shifted, E_holon)...; label = "holon")
+ lines!(ax, P_holon_analytic, E_holon_analytic; label = "holon (analytic)")
I = sortperm(momenta_shifted)
- plot!(
- p_excitations, momenta_shifted[I], E_scattering_min[I]; label = "scattering states",
- fillrange = E_scattering_max[I], fillalpha = 0.3, fillstyle = :x
+ band!(
+ ax, momenta_shifted[I], E_scattering_min[I], E_scattering_max[I];
+ label = "scattering states", color = (:gray, 0.3)
)
- p_excitations
+ Legend(fig[1, 2], ax)
+ fig
end
diff --git a/examples/groundstates/ising-cft/main.jl b/examples/groundstates/ising-cft/main.jl
index b2acd743c..2cfa46cc0 100644
--- a/examples/groundstates/ising-cft/main.jl
+++ b/examples/groundstates/ising-cft/main.jl
@@ -5,7 +5,7 @@ This tutorial is meant to show the finite size CFT spectrum for the quantum Isin
We do this by first employing an exact diagonalization technique, and then extending the analysis to larger system sizes through the use of MPS techniques.
"""
-using MPSKit, TensorKit, Plots, KrylovKit
+using MPSKit, TensorKit, CairoMakie, KrylovKit
include(joinpath(@__DIR__, "..", "..", "..", "docs", "src", "assets", "models.jl")) #hide
using .ExampleModels
using LinearAlgebra: eigvals, diagm, Hermitian
@@ -25,10 +25,7 @@ Thus, optimizing the middle tensor is equivalent to optimizing a state in the en
"""
energies, states = exact_diagonalization(H; num = 18, alg = Lanczos(; krylovdim = 200));
-plot(
- real.(energies);
- seriestype = :scatter, legend = false, ylabel = "energy", xlabel = "#eigenvalue"
-)
+scatter(real.(energies); axis = (; xlabel = "#eigenvalue", ylabel = "energy"))
md"""
!!! note "Krylov dimension"
@@ -95,14 +92,12 @@ v = 2.0
Δ = real.(energies[1:18] .- energies[1]) ./ (2π * v / L)
S = momenta ./ (2π / L)
-p = plot(
- S, real.(Δ);
- seriestype = :scatter, xlabel = "conformal spin (S)", ylabel = "scaling dimension (Δ)",
- legend = false
+fig, ax, _ = scatter(
+ S, real.(Δ); axis = (; xlabel = "conformal spin (S)", ylabel = "scaling dimension (Δ)")
)
-vline!(p, -3:3; color = "gray", linestyle = :dash)
-hline!(p, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = "gray", linestyle = :dash)
-p
+vlines!(ax, -3:3; color = :gray, linestyle = :dash)
+hlines!(ax, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = :gray, linestyle = :dash)
+fig
md"""
## Finite bond dimension
@@ -139,11 +134,9 @@ v = 2.0
Δ_mps = real.(energies_mps[1:18] .- energies_mps[1]) ./ (2π * v / L_mps)
S_mps = momenta_mps ./ (2π / L_mps)
-p_mps = plot(
- S_mps, real.(Δ_mps);
- seriestype = :scatter, xlabel = "conformal spin (S)",
- ylabel = "scaling dimension (Δ)", legend = false
+fig, ax, _ = scatter(
+ S_mps, real.(Δ_mps); axis = (; xlabel = "conformal spin (S)", ylabel = "scaling dimension (Δ)")
)
-vline!(p_mps, -3:3; color = "gray", linestyle = :dash)
-hline!(p_mps, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = "gray", linestyle = :dash)
-p_mps
+vlines!(ax, -3:3; color = :gray, linestyle = :dash)
+hlines!(ax, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color = :gray, linestyle = :dash)
+fig
diff --git a/examples/groundstates/xxz-heisenberg/main.jl b/examples/groundstates/xxz-heisenberg/main.jl
index 1ca9716c5..f6513aeb5 100644
--- a/examples/groundstates/xxz-heisenberg/main.jl
+++ b/examples/groundstates/xxz-heisenberg/main.jl
@@ -5,7 +5,7 @@ In this file we will give step by step instructions on how to analyze the spin 1
The necessary packages to follow this tutorial are:
"""
-using MPSKit, TensorKit, Plots
+using MPSKit, TensorKit, CairoMakie
include(joinpath(@__DIR__, "..", "..", "..", "docs", "src", "assets", "models.jl")) #hide
using .ExampleModels
@@ -52,7 +52,7 @@ Convergence is quite slow and even fails after sufficiently many iterations.
To understand why, we can look at the transfer matrix spectrum.
"""
-transferplot(groundstate, groundstate)
+transfer_plot(groundstate, groundstate)
md"""
We can clearly see multiple eigenvalues close to the unit circle.
@@ -89,7 +89,7 @@ The reason behind this becomes more obvious at higher bond dimensions:
groundstate, envs, info = find_groundstate(
state, H2, IDMRG2(; trunc = truncrank(50), maxiter = 20, tol = 1.0e-12, verbosity = 1)
);
-entanglementplot(groundstate)
+entanglement_plot(groundstate)
md"""
We see that some eigenvalues clearly belong to a group, and are almost degenerate.
diff --git a/examples/make.jl b/examples/make.jl
index e925df768..415ba6f6c 100644
--- a/examples/make.jl
+++ b/examples/make.jl
@@ -144,12 +144,12 @@ function strip_precompilation_output(content::AbstractString)
return join(kept, '\n')
end
-# All examples run in this one process, and Plots keeps `theme`/`default` settings as global
-# state, so an example that styles its plots would otherwise restyle every later example.
-const PLOTS_PKGID = Base.PkgId(Base.UUID("91a5bcdd-55d7-5caf-9e0b-520d859cae80"), "Plots")
+# All examples run in this one process, and Makie keeps the theme as global state,
+# so an example that styles its plots would otherwise restyle every later example.
+const MAKIE_PKGID = Base.PkgId(Base.UUID("ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"), "Makie")
function reset_plot_defaults()
- plots = get(Base.loaded_modules, PLOTS_PKGID, nothing)
- isnothing(plots) || Base.invokelatest(plots.reset_defaults)
+ makie = get(Base.loaded_modules, MAKIE_PKGID, nothing)
+ isnothing(makie) || Base.invokelatest(makie.set_theme!)
return nothing
end
diff --git a/examples/statmech/hard-hexagon/main.jl b/examples/statmech/hard-hexagon/main.jl
index 5972d088f..f935952ef 100644
--- a/examples/statmech/hard-hexagon/main.jl
+++ b/examples/statmech/hard-hexagon/main.jl
@@ -10,7 +10,7 @@ Because TensorKit supports anyonic symmetries, so does MPSKit.
To follow the tutorial you need the following packages.
"""
-using MPSKit, TensorKit, Plots, Polynomials
+using MPSKit, TensorKit, CairoMakie, Polynomials
include(joinpath(@__DIR__, "..", "..", "..", "docs", "src", "assets", "models.jl")) #hide
using .ExampleModels
@@ -99,6 +99,10 @@ c = f.coeffs[2]
#+
-p = plot(; xlabel = "logarithmic correlation length", ylabel = "entanglement entropy")
-p = plot(log.(ξs), Ss; seriestype = :scatter, label = nothing)
-plot!(p, ξ -> f(ξ) / 6; label = "fit")
+fig, ax, _ = scatter(
+ log.(ξs), Ss;
+ axis = (; xlabel = "logarithmic correlation length", ylabel = "entanglement entropy")
+)
+lines!(ax, range(extrema(log.(ξs))..., 100), ξ -> f(ξ) / 6; label = "fit")
+axislegend(ax; position = :lt)
+fig
diff --git a/examples/windowmps.jl b/examples/windowmps.jl
index cb3ffdccb..56b9a3e2b 100644
--- a/examples/windowmps.jl
+++ b/examples/windowmps.jl
@@ -1,4 +1,4 @@
-using MPSKit, TensorKit, Plots
+using MPSKit, TensorKit, CairoMakie
include(joinpath(@__DIR__, "..", "docs", "src", "assets", "models.jl"))
using .ExampleModels
diff --git a/ext/MPSKitMakieExt.jl b/ext/MPSKitMakieExt.jl
new file mode 100644
index 000000000..283cd7cab
--- /dev/null
+++ b/ext/MPSKitMakieExt.jl
@@ -0,0 +1,112 @@
+module MPSKitMakieExt
+
+using Makie
+using MPSKit: MPSKit, AbstractMPS, _firstspace
+using TensorKit: dim
+import MPSKit: entanglement_plot, entanglement_plot!, transfer_plot, transfer_plot!
+
+# Makie lowercases the type name to get the function name, hence the underscores
+@recipe Entanglement_Plot (mps::AbstractMPS,) begin
+ site = 0
+ expand_symmetry = false
+ sortby = maximum
+ sector_margin = 1 // 10
+ sector_formatter = string
+ colors = Makie.wong_colors()
+ markersize = @inherit markersize
+end
+
+function entanglement_data(mps, site, expand_symmetry, sortby, sector_margin)
+ data = MPSKit.entanglement_plot_data(mps, site; expand_symmetry, sortby, sector_margin)
+ return (; data.sectors, points = map((x, y) -> Point2d.(x, y), data.xs, data.ys), data.dropped)
+end
+
+function Makie.plot!(p::Entanglement_Plot)
+ data = lift(entanglement_data, p, p.mps, p.site, p.expand_symmetry, p.sortby, p.sector_margin)
+ data[].dropped && @warn "Entanglement spectrum contains vanishing Schmidt values. These are omitted from the plot."
+ add_sector_scatters!(p, data)
+ return p
+end
+
+function Makie.preferred_axis_attributes(::Type{Axis}, p::Entanglement_Plot)
+ mps, site = p.mps[], p.site[]
+ sectors, points = entanglement_data(mps, site, p.expand_symmetry[], p.sortby[], p.sector_margin[])
+ n = length(sectors)
+ smallest = minimum(last, Iterators.flatten(points); init = 1.0)
+ # limits must precede `yscale`: an empty axis has limits (0, 10), invalid on a log scale
+ return (;
+ limits = ((1, n + 1), (exp10(floor(log10(smallest))), 1.1)),
+ title = "Entanglement Spectrum",
+ xlabel = "χ = $(dim(_firstspace(mps.C[site])))",
+ ylabel = "Schmidt values",
+ xticks = ((1:n) .+ 0.5, p.sector_formatter[].(sectors)),
+ xticklabelrotation = π / 4,
+ xgridvisible = false,
+ yscale = log10,
+ )
+end
+
+Makie.get_plots(p::Entanglement_Plot) = p.plots
+
+@recipe Transfer_Plot (above::AbstractMPS, below::AbstractMPS) begin
+ sectors = nothing
+ transferkwargs = (;)
+ thetaorigin = 0.0
+ sector_formatter = string
+ colors = Makie.wong_colors()
+ markersize = @inherit markersize
+end
+
+Makie.convert_arguments(::Type{<:Transfer_Plot}, above::AbstractMPS) = (above, above)
+
+function transfer_data(above, below, sectors, transferkwargs, thetaorigin)
+ # Makie stores NamedTuple attributes as nested `Attributes` of observables
+ transferkwargs = (; (k => to_value(v) for (k, v) in pairs(transferkwargs))...)
+ data = MPSKit.transfer_plot_data(above, below; sectors, transferkwargs, thetaorigin)
+ return (; data.sectors, points = map((θ, r) -> Point2d.(θ, r), data.θs, data.rs))
+end
+
+function Makie.plot!(p::Transfer_Plot)
+ data = lift(transfer_data, p, p.above, p.below, p.sectors, p.transferkwargs, p.thetaorigin)
+ add_sector_scatters!(p, data)
+ return p
+end
+
+function Makie.preferred_axis_attributes(::Type{Axis}, p::Transfer_Plot)
+ θ₀ = float(p.thetaorigin[])
+ ks = ceil(Int, 4θ₀ / π):floor(Int, 4(θ₀ + 2π) / π)
+ return (;
+ title = "Transfer Spectrum",
+ xlabel = "θ",
+ ylabel = "r",
+ xticks = (ks .* (π / 4), pilabel.(ks .// 4)),
+ limits = ((θ₀ - 0.1, θ₀ + 2π + 0.1), (nothing, 1.05)),
+ )
+end
+
+Makie.get_plots(p::Transfer_Plot) = p.plots
+
+# one labelled scatter per sector so `axislegend` lists them; the sector count is fixed at creation
+function add_sector_scatters!(p, data)
+ sectors = data[].sectors
+ isempty(sectors) && return scatter!(p, Point2d[]; visible = false)
+ colors = p.colors[]
+ for (i, sector) in enumerate(sectors)
+ scatter!(
+ p, lift(d -> d.points[i], p, data);
+ color = colors[mod1(i, length(colors))],
+ markersize = p.markersize,
+ label = p.sector_formatter[](sector),
+ )
+ end
+ return p
+end
+
+function pilabel(x::Rational)
+ n, d = numerator(x), denominator(x)
+ iszero(n) && return "0"
+ s = (abs(n) == 1 ? (n < 0 ? "-" : "") : string(n)) * "π"
+ return d == 1 ? s : "$s/$d"
+end
+
+end
diff --git a/ext/MPSKitPlotsExt.jl b/ext/MPSKitPlotsExt.jl
new file mode 100644
index 000000000..0b24d4465
--- /dev/null
+++ b/ext/MPSKitPlotsExt.jl
@@ -0,0 +1,87 @@
+module MPSKitPlotsExt
+
+using RecipesBase
+using MPSKit, TensorKit
+import MPSKit: entanglementplot, transferplot
+
+@userplot EntanglementPlot
+
+@recipe function f(
+ h::EntanglementPlot; site = 0, expand_symmetry = false, sortby = maximum,
+ sector_margin = 1 // 10, sector_formatter = string
+ )
+ Base.depwarn(
+ "`entanglementplot` (Plots.jl) is deprecated, use `entanglement_plot` (Makie.jl) instead.",
+ :entanglementplot
+ )
+ mps = h.args[1]
+ (; sectors, xs, ys, dropped) = MPSKit.entanglement_plot_data(
+ mps, site; expand_symmetry, sortby, sector_margin
+ )
+ dropped && @warn "Entanglement spectrum contains vanishing Schmidt values. These are omitted from the plot."
+
+ for (x, y, sector) in zip(xs, ys, sectors)
+ @series begin
+ seriestype := :scatter
+ label := sector_formatter(sector)
+ return x, y
+ end
+ end
+
+ title --> "Entanglement Spectrum"
+ legend --> false
+ grid --> :xy
+ widen --> true
+
+ xguide --> "χ = $(dim(MPSKit._firstspace(mps.C[site])))"
+ xticks --> (1:length(sectors), sector_formatter.(sectors))
+ xtickfonthalign --> :center
+ xtick_direction --> :out
+ xrotation --> 45
+ xlims --> (1, length(sectors) + 1)
+
+ ylims --> (-Inf, 1 + 1.0e-1)
+ yscale --> :log10
+ label := nothing
+
+ return []
+end
+
+@userplot TransferPlot
+
+@recipe function f(
+ h::TransferPlot; sectors = nothing, transferkwargs = (;), thetaorigin = 0,
+ sector_formatter = string
+ )
+ Base.depwarn(
+ "`transferplot` (Plots.jl) is deprecated, use `transfer_plot` (Makie.jl) instead.",
+ :transferplot
+ )
+ data = MPSKit.transfer_plot_data(h.args...; sectors, transferkwargs, thetaorigin)
+
+ for (θ, r, sector) in zip(data.θs, data.rs, data.sectors)
+ @series begin
+ yguide --> "r"
+ ylims --> (-Inf, 1.05)
+
+ xguide --> "θ"
+ xlims --> (thetaorigin, thetaorigin + 2pi)
+ xticks --> range(0, 2pi; length = 7)
+ xformatter --> x -> "$(rationalize(x / π, tol = 0.05))π"
+ xwiden --> true
+ seriestype := :scatter
+ markershape --> :auto
+ label := sector_formatter(sector)
+ return θ, r
+ end
+ end
+
+ title --> "Transfer Spectrum"
+ legend --> false
+ grid --> :xy
+ framestyle --> :zerolines
+
+ return nothing
+end
+
+end
diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index 468e15a46..e7b2c253b 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -53,6 +53,7 @@ export correlation_length, marek_gap, transfer_spectrum
export entropy, entanglement_spectrum
export open_boundary_conditions, periodic_boundary_conditions
export entanglementplot, transferplot
+export entanglement_plot, entanglement_plot!, transfer_plot, transfer_plot!
export r_LL, l_LL, r_RR, l_RR, r_RL, r_LR, l_RL, l_LR # TODO: rename
# unexported
@@ -80,7 +81,6 @@ using KrylovKit: KrylovAlgorithm
using OptimKit
using Base.Threads
using Base.Iterators
-using RecipesBase
using VectorInterface
using Accessors
using HalfIntegers
diff --git a/src/utility/plotting.jl b/src/utility/plotting.jl
index b58b912cc..992484747 100644
--- a/src/utility/plotting.jl
+++ b/src/utility/plotting.jl
@@ -1,7 +1,11 @@
"""
- entanglementplot(state; site = 0[, kwargs...])
+ entanglement_plot(state; site = 0, kwargs...)
+ entanglement_plot(gridposition, state; kwargs...)
+ entanglement_plot!([axis], state; kwargs...)
-Plot the entanglement spectrum (see [`entanglement_spectrum`](@ref)) of a given MPS `state`.
+Plot the entanglement spectrum (see [`entanglement_spectrum`](@ref)) of a given MPS `state`
+using [Makie.jl](https://github.com/MakieOrg/Makie.jl). The Schmidt values of each sector
+are drawn as a group of points, on a logarithmic scale.
# Arguments
@@ -15,83 +19,90 @@ Plot the entanglement spectrum (see [`entanglement_spectrum`](@ref)) of a given
- `sortby = maximum`: the method of sorting the sectors.
- `sector_margin = 1 // 10`: the amount of whitespace between sectors.
- `sector_formatter = string`: how to convert sectors to strings.
-- `kwargs...`: other kwargs are passed on to the plotting backend.
+- `colors = Makie.wong_colors()`: the colors that are cycled through for the sectors.
+
+All other keywords are handled by Makie as usual, e.g. `axis = (; title = "custom title")`
+or `figure = (; size = (800, 600))` for the non-mutating version. The mutating version does
+not modify the axis it plots into, unless `use_axis_hints = true` is passed.
+Each sector is labelled, such that `axislegend(ax)` adds a legend.
!!! note
- You will need to manually import [Plots.jl](https://github.com/JuliaPlots/Plots.jl) to
- be able to use this function. MPSKit.jl defines its plots based on
- [RecipesBase.jl](https://github.com/JuliaPlots/Plots.jl/tree/v2/RecipesBase), but the
- user still has to add `using Plots` to be able to actually produce the plots.
+ This function requires a Makie backend to be loaded, e.g. `using CairoMakie`.
+"""
+function entanglement_plot end
+function entanglement_plot! end
"""
-function entanglementplot end
-@userplot EntanglementPlot
+ transfer_plot(above, below = above; sectors = nothing, transferkwargs = (;), kwargs...)
+ transfer_plot(gridposition, above, below = above; kwargs...)
+ transfer_plot!([axis], above, below = above; kwargs...)
-@recipe function f(
- h::EntanglementPlot; site = 0, expand_symmetry = false, sortby = maximum,
- sector_margin = 1 // 10, sector_formatter = string
- )
- mps = h.args[1]
+Plot the partial transfer matrix spectrum (see [`transfer_spectrum`](@ref)) of two
+`InfiniteMPS`s using [Makie.jl](https://github.com/MakieOrg/Makie.jl), as the moduli
+against the angles of the eigenvalues.
- spectra = entanglement_spectrum(mps, site)
- sectors = []
- spectrum = []
- for (c, b) in pairs(spectra)
- if expand_symmetry # Duplicate entries according to the quantum dimension.
- b′ = repeat(b, dim(c))
- sort!(b′; rev = true)
- push!(spectrum, b′)
- else
- push!(spectrum, b)
- end
- push!(sectors, c)
- end
+# Arguments
- if length(spectrum) > 1
- order = sortperm(spectrum; by = sortby, rev = true)
- spectrum = spectrum[order]
- sectors = sectors[order]
- end
+- `above::InfiniteMPS`: above mps for [`transfer_spectrum`](@ref).
+- `below::InfiniteMPS = above`: below mps for [`transfer_spectrum`](@ref).
- for (i, (partial_spectrum, sector)) in enumerate(zip(spectrum, sectors))
- @series begin
- seriestype := :scatter
- label := sector_formatter(sector)
- n_spectrum = length(partial_spectrum)
-
- # Put single dot in the middle, or a linear range with padding.
- if n_spectrum == 1
- x = [i + 1 // 2]
- else
- x = range(i + sector_margin, i + 1 - sector_margin; length = n_spectrum)
- end
- return x, partial_spectrum
- end
- end
+# Keyword Arguments
+
+- `sectors = nothing`: restrict the spectrum to the given sectors; by default all sectors of
+ the transfer space are included.
+- `transferkwargs = (;)`: kwargs for call to [`transfer_spectrum`](@ref).
+- `thetaorigin = 0`: origin of the angle range.
+- `sector_formatter = string`: how to convert sectors to strings.
+- `colors = Makie.wong_colors()`: the colors that are cycled through for the sectors.
+
+All other keywords are handled by Makie as usual, e.g. `axis = (; title = "custom title")`
+or `figure = (; size = (800, 600))` for the non-mutating version. The mutating version does
+not modify the axis it plots into, unless `use_axis_hints = true` is passed.
+Each sector is labelled, such that `axislegend(ax)` adds a legend.
- title --> "Entanglement Spectrum"
- legend --> false
- grid --> :xy
- widen --> true
+!!! note
+ This function requires a Makie backend to be loaded, e.g. `using CairoMakie`.
+"""
+function transfer_plot end
+function transfer_plot! end
- xguide --> "χ = $(dim(_firstspace(mps.C[site])))"
- xticks --> (1:length(sectors), sector_formatter.(sectors))
- xtickfonthalign --> :center
- xtick_direction --> :out
- xrotation --> 45
- xlims --> (1, length(sectors) + 1)
+"""
+ entanglementplot(state; site = 0[, kwargs...])
- ylims --> (-Inf, 1 + 1.0e-1)
- yscale --> :log10
- label := nothing
+Plot the entanglement spectrum (see [`entanglement_spectrum`](@ref)) of a given MPS `state`
+using [Plots.jl](https://github.com/JuliaPlots/Plots.jl).
- return []
-end
+!!! warning "Deprecated"
+ The Plots.jl recipes are deprecated, use the Makie.jl based [`entanglement_plot`](@ref) instead.
+
+# Arguments
+
+- `state`: the MPS for which to compute the entanglement spectrum.
+
+# Keyword Arguments
+
+- `site::Int = 0`: MPS index for multisite unit cells. The spectrum is computed for the bond
+ between `site` and `site + 1`.
+- `expand_symmetry = false`: add quantum dimension degeneracies.
+- `sortby = maximum`: the method of sorting the sectors.
+- `sector_margin = 1 // 10`: the amount of whitespace between sectors.
+- `sector_formatter = string`: how to convert sectors to strings.
+- `kwargs...`: other kwargs are passed on to the plotting backend.
+
+!!! note
+ You will need to manually import [Plots.jl](https://github.com/JuliaPlots/Plots.jl) to
+ be able to use this function.
+"""
+function entanglementplot end
"""
transferplot(above, below = above; sectors = nothing, transferkwargs = (;)[, kwargs...])
-Plot the partial transfer matrix spectrum of two InfiniteMPS's.
+Plot the partial transfer matrix spectrum of two InfiniteMPS's using
+[Plots.jl](https://github.com/JuliaPlots/Plots.jl).
+
+!!! warning "Deprecated"
+ The Plots.jl recipes are deprecated, use the Makie.jl based [`transfer_plot`](@ref) instead.
# Arguments
@@ -109,50 +120,44 @@ Plot the partial transfer matrix spectrum of two InfiniteMPS's.
!!! note
You will need to manually import [Plots.jl](https://github.com/JuliaPlots/Plots.jl) to
- be able to use this function. MPSKit.jl defines its plots based on
- [RecipesBase.jl](https://github.com/JuliaPlots/Plots.jl/tree/v2/RecipesBase), but the
- user still has to add `using Plots` to be able to actually produce the plots.
-
+ be able to use this function.
"""
function transferplot end
-@userplot TransferPlot
-@recipe function f(
- h::TransferPlot; sectors = nothing, transferkwargs = (;), thetaorigin = 0,
- sector_formatter = string
+# Backend-independent data for the plotting extensions, one entry per sector.
+function entanglement_plot_data(
+ mps, site = 0; expand_symmetry = false, sortby = maximum, sector_margin = 1 // 10
+ )
+ spectra = entanglement_spectrum(mps, site)
+ sectors = collect(keys(spectra))
+ ys = map(sectors) do c
+ b = collect(Float64, spectra[c])
+ return expand_symmetry ? sort!(repeat(b, dim(c)); rev = true) : b
+ end
+ order = sortperm(ys; by = sortby, rev = true)
+ sectors, ys = sectors[order], ys[order]
+ dropped = any(y -> any(<=(0), y), ys)
+ ys = [filter(>(0), y) for y in ys]
+ margin = float(sector_margin)
+ xs = map(enumerate(ys)) do (i, y)
+ n = length(y)
+ return n == 1 ? [i + 0.5] : collect(range(i + margin, i + 1 - margin; length = n))
+ end
+ return (; sectors, xs, ys, dropped)
+end
+
+function transfer_plot_data(
+ above, below = above; sectors = nothing, transferkwargs = (;), thetaorigin = 0
)
- below = length(h.args) == 1 ? h.args[1] : h.args[2]
kwargs = (; transferkwargs...)
if sectors !== nothing && get(kwargs, :howmany, 20) isa Int
# restrict the computation to the requested sectors
howmany = Dict(c => get(kwargs, :howmany, 20) for c in sectors)
kwargs = (; kwargs..., howmany)
end
- spectra = transfer_spectrum(h.args[1], below; kwargs...)
-
- for (sector, spectrum) in pairs(spectra)
- sectors === nothing || sector in sectors || continue
-
- @series begin
- yguide --> "r"
- ylims --> (-Inf, 1.05)
-
- xguide --> "θ"
- xlims --> (thetaorigin, thetaorigin + 2pi)
- xticks --> range(0, 2pi; length = 7)
- xformatter --> x -> "$(rationalize(x / π, tol = 0.05))π"
- xwiden --> true
- seriestype := :scatter
- markershape --> :auto
- label := sector_formatter(sector)
- return mod2pi.(angle.(spectrum) .+ thetaorigin) .- thetaorigin, abs.(spectrum)
- end
- end
-
- title --> "Transfer Spectrum"
- legend --> false
- grid --> :xy
- framestyle --> :zerolines
-
- return nothing
+ spectra = transfer_spectrum(above, below; kwargs...)
+ plotted = [c for c in keys(spectra) if sectors === nothing || c in sectors]
+ θs = [mod2pi.(angle.(spectra[c]) .- thetaorigin) .+ thetaorigin for c in plotted]
+ rs = [abs.(spectra[c]) for c in plotted]
+ return (; sectors = plotted, θs, rs)
end
diff --git a/test/Project.toml b/test/Project.toml
index 5597b8c84..7520013a3 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -8,6 +8,7 @@ BlockTensorKit = "5f87ffc2-9cf1-4a46-8172-465d160bd8cd"
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
@@ -27,6 +28,7 @@ MPSKit = {path = ".."}
[compat]
AMDGPU = "2"
Aqua = "0.8.9"
+CairoMakie = "0.15"
CUDA = "6"
Combinatorics = "1"
ParallelTestRunner = "2"
diff --git a/test/internals/makie.jl b/test/internals/makie.jl
new file mode 100644
index 000000000..7b52dd134
--- /dev/null
+++ b/test/internals/makie.jl
@@ -0,0 +1,49 @@
+println("
+-----------------------------------
+| Plot tests with Makie.jl |
+-----------------------------------
+")
+
+using .TestSetup
+using Test, TestExtras
+using MPSKit
+using TensorKit
+using TensorKit: ℙ
+using CairoMakie
+using Plots: Plots
+
+@testset "plot tests" begin
+ ψ = InfiniteMPS([ℙ^2], [ℙ^5])
+ _, ax, _ = entanglement_plot(ψ)
+ @test ax.title[] == "Entanglement Spectrum"
+ _, ax, _ = transfer_plot(ψ, ψ; transferkwargs = (; howmany = 3), axis = (; title = "custom"))
+ @test ax.title[] == "custom"
+
+ # mutating plots only style the axis on request
+ fig = Figure()
+ ax = Axis(fig[1, 1]; title = "mine")
+ @test transfer_plot!(ax, ψ) isa Makie.Plot
+ @test ax.title[] == "mine"
+ ax = Axis(fig[1, 2])
+ entanglement_plot!(ax, ψ; use_axis_hints = true)
+ @test ax.title[] == "Entanglement Spectrum"
+end
+
+@testset "graded plots" begin
+ ψ = InfiniteMPS([Z2Space(0 => 1, 1 => 1)], [Z2Space(0 => 4, 1 => 4)])
+ @test entanglement_plot(ψ; site = 1, expand_symmetry = true) isa Makie.FigureAxisPlot
+
+ _, ax, p = transfer_plot(ψ)
+ @test length(p.plots) == 2
+ @test axislegend(ax) isa Legend
+ _, _, p = transfer_plot(ψ; sectors = [unit(sectortype(ψ))])
+ @test length(p.plots) == 1
+ θs = MPSKit.transfer_plot_data(ψ; thetaorigin = -π).θs
+ @test all(θ -> -π <= θ < π, reduce(vcat, θs))
+end
+
+@testset "coexistence with Plots.jl" begin
+ ψ = InfiniteMPS([ℙ^2], [ℙ^5])
+ @test entanglement_plot(ψ) isa Makie.FigureAxisPlot
+ @test (@test_deprecated MPSKit.entanglementplot(ψ)) isa Plots.Plot
+end
diff --git a/test/internals/plots.jl b/test/internals/plots.jl
index 412d69d93..a9a253f2d 100644
--- a/test/internals/plots.jl
+++ b/test/internals/plots.jl
@@ -1,7 +1,7 @@
println("
----------------------
-| Plot tests |
----------------------
+-----------------------------------
+| Plot tests with Plots.jl |
+-----------------------------------
")
using .TestSetup
@@ -13,6 +13,6 @@ using Plots
@testset "plot tests" begin
ψ = InfiniteMPS([ℙ^2], [ℙ^5])
- @test transferplot(ψ) isa Plots.Plot
- @test entanglementplot(ψ) isa Plots.Plot
+ @test (@test_deprecated transferplot(ψ)) isa Plots.Plot
+ @test (@test_deprecated entanglementplot(ψ)) isa Plots.Plot
end