Improve quasiparticle type stability (#272)

* Quasiparticle type stability improvements

Test type stability excitations

* Add Multiline VectorInterface support

* Add multiline QP constructor

* utility `eltype` in type domain

* More utility functions

* Replace `@constinferred` with `@inferred`

* Refactor QP to improve type stability (again)

* various fixes

* restore `effective_excitation_hamiltonian`

* attempt to fix type stability for lts

Author: lkdvos

src/algorithms/excitation/quasiparticleexcitation.jl

@@ -150,7 +150,8 @@ end
 function environments(exci::FiniteQP,
                       H::FiniteMPOHamiltonian,
                       lenvs=environments(exci.left_gs, H),
-                      renvs=exci.trivial ? lenvs : environments(exci.right_gs, H))
+                      renvs=exci.trivial ? lenvs : environments(exci.right_gs, H);
+                      kwargs...)
     AL = exci.left_gs.AL
     AR = exci.right_gs.AR
 

src/environments/qp_envs.jl

@@ -214,7 +214,8 @@ Utility
 
 Base.size(ψ::InfiniteMPS, args...) = size(ψ.AL, args...)
 Base.length(ψ::InfiniteMPS) = length(ψ.AL)
-Base.eltype(ψ::InfiniteMPS) = eltype(ψ.AL)
+Base.eltype(ψ::InfiniteMPS) = eltype(typeof(ψ))
+Base.eltype(::Type{<:InfiniteMPS{A}}) where {A} = A
 Base.copy(ψ::InfiniteMPS) = InfiniteMPS(copy(ψ.AL), copy(ψ.AR), copy(ψ.C), copy(ψ.AC))
 function Base.copy!(ψ::InfiniteMPS, ϕ::InfiniteMPS)
     copy!.(ψ.AL, ϕ.AL)

src/states/infinitemps.jl

@@ -74,6 +74,10 @@ bond_type(::Type{Multiline{S}}) where {S} = bond_type(S)
 site_type(st::Multiline) = site_type(typeof(st))
 bond_type(st::Multiline) = bond_type(typeof(st))
 VectorInterface.scalartype(::Multiline{T}) where {T} = scalartype(T)
+TensorKit.sectortype(t::Multiline) = sectortype(typeof(t))
+TensorKit.sectortype(::Type{Multiline{T}}) where {T} = sectortype(T)
+TensorKit.spacetype(t::Multiline) = spacetype(typeof(t))
+TensorKit.spacetype(::Type{Multiline{T}}) where {T} = spacetype(T)
 
 function TensorKit.dot(a::MultilineMPS, b::MultilineMPS; kwargs...)
     return sum(dot.(parent(a), parent(b); kwargs...))

src/states/multilinemps.jl

@@ -26,22 +26,40 @@ struct RightGaugedQP{S,T1,T2,E<:Number}
     momentum::E
 end
 
+function leftgaugedqptype(::Type{S}, ::Type{E}) where {S,E<:Number}
+    T1 = eltype(S)
+    T2 = tensormaptype(spacetype(T1), 1, 2, storagetype(T1))
+    return LeftGaugedQP{S,T1,T2,E}
+end
+
 #constructors
 function LeftGaugedQP(datfun, left_gs, right_gs=left_gs;
                       sector=one(sectortype(left_gs)), momentum=0.0)
     # find the left null spaces for the TNS
     excitation_space = Vect[typeof(sector)](sector => 1)
-    VLs = [adjoint(rightnull(adjoint(v))) for v in left_gs.AL]
-    Xs = [TensorMap{scalartype(left_gs)}(undef, _lastspace(VLs[loc])',
-                                         excitation_space' *
-                                         right_virtualspace(right_gs, loc))
-          for loc in 1:length(left_gs)]
-    fill_data!.(Xs, datfun)
+    VLs = convert(Vector, map(leftnull, left_gs.AL))
+    Xs = map(enumerate(VLs)) do (loc, vl)
+        x = similar(vl,
+                    right_virtualspace(vl) ←
+                    excitation_space' ⊗ right_virtualspace(right_gs, loc))
+        fill_data!(x, datfun)
+        return x
+    end
     left_gs isa InfiniteMPS ||
         momentum == zero(momentum) ||
         @warn "momentum is ignored for finite quasiparticles"
     return LeftGaugedQP(left_gs, right_gs, VLs, Xs, momentum)
 end
+function LeftGaugedQP(datfun, left_gs::MultilineMPS, right_gs::MultilineMPS=left_gs;
+                      sector=one(sectortype(left_gs)), momentum=0.0)
+    # not sure why this is needed for type stability
+    Tresult = leftgaugedqptype(eltype(parent(left_gs)), typeof(momentum))
+    qp_rows = Vector{Tresult}(undef, size(left_gs, 1))
+    for row in eachindex(qp_rows)
+        qp_rows[row] = LeftGaugedQP(datfun, left_gs[row], right_gs[row]; sector, momentum)
+    end
+    return Multiline(qp_rows)
+end
 
 function RightGaugedQP(datfun, left_gs, right_gs=left_gs;
                        sector=one(sectortype(left_gs)), momentum=0.0)
@@ -59,13 +77,11 @@ function RightGaugedQP(datfun, left_gs, right_gs=left_gs;
 end
 
 #gauge dependent code
-function Base.similar(v::LeftGaugedQP, T=scalartype(v))
-    return LeftGaugedQP(v.left_gs, v.right_gs, v.VLs, map(e -> similar(e, T), v.Xs),
-                        v.momentum)
+function Base.similar(v::LeftGaugedQP, ::Type{T}=scalartype(v)) where {T<:Number}
+    return LeftGaugedQP(v.left_gs, v.right_gs, v.VLs, similar.(v.Xs, T), v.momentum)
 end
-function Base.similar(v::RightGaugedQP, T=scalartype(v))
-    return RightGaugedQP(v.left_gs, v.right_gs, map(e -> similar(e, T), v.Xs), v.VRs,
-                         v.momentum)
+function Base.similar(v::RightGaugedQP, ::Type{T}=scalartype(v)) where {T<:Number}
+    return RightGaugedQP(v.left_gs, v.right_gs, similar.(v.Xs, T), v.VRs, v.momentum)
 end
 
 Base.getindex(v::LeftGaugedQP, i::Int) = v.VLs[mod1(i, end)] * v.Xs[mod1(i, end)];

src/states/quasiparticle_state.jl

@@ -51,3 +51,49 @@ function Base.repeat(A::Multiline, rows::Int, cols::Int)
     outer = repeat(inner, rows)
     return Multiline(outer)
 end
+
+# VectorInterface
+# ---------------
+VectorInterface.scalartype(::Type{Multiline{T}}) where {T} = scalartype(T)
+
+function VectorInterface.zerovector(x::Multiline, ::Type{S}) where {S<:Number}
+    return Multiline(zerovector.(parent(x), S))
+end
+VectorInterface.zerovector!(x::Multiline) = (zerovector!.(parent(x)); x)
+
+function VectorInterface.scale(x::Multiline, α::Number)
+    return scale!(zerovector(x, VectorInterface.promote_scale(x, α)), x, α)
+end
+
+function VectorInterface.scale!(x::Multiline, α::Number)
+    scale!.(parent(x), α)
+    return x
+end
+VectorInterface.scale!!(x::Multiline, α::Number) = scale!(x, α)
+
+function VectorInterface.scale!(x::Multiline, x′::Multiline, α::Number)
+    scale!.(parent(x), parent(x′), α)
+    return x
+end
+
+VectorInterface.scale!!(x::Multiline, x′::Multiline, α::Number) = scale!(x, x′, α)
+
+function VectorInterface.add(x::Multiline, y::Multiline, α::Number, β::Number)
+    z = zerovector(x, VectorInterface.promote_add(x, y, α, β))
+    return add!(scale!(z, x, β), y, α)
+end
+
+function VectorInterface.add!(x::Multiline, y::Multiline, α::Number, β::Number)
+    add!.(parent(x), parent(y), α, β)
+    return x
+end
+
+VectorInterface.add!!(x::Multiline, y::Multiline, α::Number, β::Number) = add!(x, y, α, β)
+
+function VectorInterface.inner(x::Multiline, y::Multiline)
+    T = VectorInterface.promote_inner(x, y)
+    init = zero(T)
+    return sum(splat(inner), zip(parent(parent(x)), parent(parent(y))); init)
+end
+
+LinearAlgebra.norm(x::Multiline) = sqrt(real(inner(x, x)))

src/utility/multiline.jl

@@ -435,7 +435,8 @@ end
         ψ = InfiniteMPS([ℙ^3, ℙ^3], [ℙ^48, ℙ^48])
         ψ, envs, _ = find_groundstate(ψ, H; maxiter=400, verbosity=verbosity_conv,
                                       tol=1e-10)
-        energies, ϕs = excitations(H, QuasiparticleAnsatz(), Float64(pi), ψ, envs)
+        energies, ϕs = @inferred excitations(H, QuasiparticleAnsatz(), Float64(pi), ψ,
+                                             envs)
         @test energies[1] ≈ 0.41047925 atol = 1e-4
         @test variance(ϕs[1], H) < 1e-8
     end
@@ -445,8 +446,9 @@ end
         ψ, envs, _ = leading_boundary(ψ, H,
                                       VUMPS(; maxiter=400, verbosity=verbosity_conv,
                                             tol=1e-10))
-        energies, ϕs = excitations(H, QuasiparticleAnsatz(), [0.0, Float64(pi / 2)], ψ,
-                                   envs; verbosity=0)
+        energies, ϕs = @inferred excitations(H, QuasiparticleAnsatz(),
+                                             [0.0, Float64(pi / 2)], ψ,
+                                             envs; verbosity=0)
         @test abs(energies[1]) > abs(energies[2]) # has a minimum at pi/2
     end
 
@@ -455,7 +457,7 @@ end
         H_inf = force_planar(transverse_field_ising())
         ψ_inf = InfiniteMPS([ℙ^2], [ℙ^10])
         ψ_inf, envs, _ = find_groundstate(ψ_inf, H_inf; maxiter=400, verbosity, tol=1e-9)
-        energies, ϕs = excitations(H_inf, QuasiparticleAnsatz(), 0.0, ψ_inf, envs)
+        energies, ϕs = @inferred excitations(H_inf, QuasiparticleAnsatz(), 0.0, ψ_inf, envs)
         inf_en = energies[1]
 
         fin_en = map([20, 10]) do len
@@ -464,20 +466,20 @@ end
             ψ, envs, = find_groundstate(ψ, H; verbosity)
 
             # find energy with quasiparticle ansatz
-            energies_QP, ϕs = excitations(H, QuasiparticleAnsatz(), ψ, envs)
+            energies_QP, ϕs = @inferred excitations(H, QuasiparticleAnsatz(), ψ, envs)
             @test variance(ϕs[1], H) < 1e-6
 
             # find energy with normal dmrg
             for gsalg in (DMRG(; verbosity, tol=1e-6),
                           DMRG2(; verbosity, tol=1e-6, trscheme=truncbelow(1e-4)))
-                energies_dm, _ = excitations(H, FiniteExcited(; gsalg), ψ)
+                energies_dm, _ = @inferred excitations(H, FiniteExcited(; gsalg), ψ)
                 @test energies_dm[1] ≈ energies_QP[1] + expectation_value(ψ, H, envs) atol = 1e-4
             end
 
             # find energy with Chepiga ansatz
-            energies_ch, _ = excitations(H, ChepigaAnsatz(), ψ, envs)
+            energies_ch, _ = @inferred excitations(H, ChepigaAnsatz(), ψ, envs)
             @test energies_ch[1] ≈ energies_QP[1] + expectation_value(ψ, H, envs) atol = 1e-4
-            energies_ch2, _ = excitations(H, ChepigaAnsatz2(), ψ, envs)
+            energies_ch2, _ = @inferred excitations(H, ChepigaAnsatz2(), ψ, envs)
             @test energies_ch2[1] ≈ energies_QP[1] + expectation_value(ψ, H, envs) atol = 1e-4
             return energies_QP[1]
         end