Skip to content

Commit c322570

Browse files
lkdvoslkdvos
authored
Restore "expert-mode" Hamiltonian constructors (#257)
* Implement FiniteMPOHamiltonian(W_matrices) constructor * Add (limited) testset * fix typo * fix unbound types * make types function more generically * Add InfiniteMPOHamiltonian version * Add test * fix typo
1 parent 941b747 commit c322570

File tree

3 files changed

+280
-2
lines changed

3 files changed

+280
-2
lines changed

‎src/operators/mpohamiltonian.jl‎

Lines changed: 239 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -54,7 +54,245 @@ function InfiniteMPOHamiltonian(Ws::AbstractVector{O}) where {O<:MPOTensor}
5454
return InfiniteMPOHamiltonian{O}(Ws)
5555
end
5656

57-
# TODO: consider if we want MPOHamiltonian(x::AbstractArray{<:Any,3}) constructor
57+
"""
58+
FiniteMPOHamiltonian(Ws::Vector{<:Matrix})
59+
60+
Create a `FiniteMPOHamiltonian` from a vector of matrices, such that `Ws[i][j, k]` represents
61+
the the operator at site `i`, left level `j` and right level `k`. Here, the entries can be
62+
either `MPOTensor`, `Missing` or `Number`.
63+
"""
64+
function FiniteMPOHamiltonian(Ws::Vector{<:Matrix})
65+
T = promote_type(_split_mpoham_types.(Ws)...)
66+
return FiniteMPOHamiltonian{T}(Ws)
67+
end
68+
function FiniteMPOHamiltonian{O}(W_mats::Vector{<:Matrix}) where {O<:MPOTensor}
69+
T = scalartype(O)
70+
L = length(W_mats)
71+
# initialize sumspaces
72+
S = spacetype(O)
73+
Vspaces = Vector{SumSpace{S}}(undef, L + 1)
74+
Pspaces = Vector{S}(undef, L)
75+
76+
# left end
77+
nlvls = size(W_mats[1], 1)
78+
@assert nlvls == 1 "left boundary should have a single level"
79+
Vspaces[1] = SumSpace(oneunit(S))
80+
# right end
81+
nlvls = size(W_mats[end], 2)
82+
@assert nlvls == 1 "right boundary should have a single level"
83+
Vspaces[end] = SumSpace(oneunit(S))
84+
85+
# start filling spaces
86+
# note that we assume that the FSA does not contain "dead ends", as this would mess with the
87+
# ability to deduce spaces
88+
for (site, W_mat) in enumerate(W_mats)
89+
# physical space
90+
operator_id = findfirst(x -> x isa O, W_mat)
91+
@assert !isnothing(operator_id) "could not determine physical space at site $site"
92+
Pspaces[site] = physicalspace(W_mat[operator_id])
93+
94+
Vs_left = Vspaces[site]
95+
if site == L
96+
Vs_right = Vspaces[site + 1]
97+
else
98+
# start by assuming trivial spaces everywhere -- replace everything that we know
99+
# assume spacecheck errors will happen when filling the BlockTensors
100+
nlvls = size(W_mat, 2)
101+
Vs_right = SumSpace(fill(oneunit(S), nlvls))
102+
end
103+
104+
for I in eachindex(IndexCartesian(), W_mat)
105+
Welem = W_mat[I]
106+
ismissing(Welem) && continue
107+
row, col = I.I
108+
if Welem isa MPOTensor
109+
V_left = left_virtualspace(Welem)
110+
@assert Vs_left[row] == V_left "incompatible space between sites $(site-1) and $site at level $row"
111+
V_right = right_virtualspace(Welem)
112+
Vs_right[col] = V_right
113+
elseif !iszero(Welem) # Welem isa Number
114+
V_left = V_right = Vs_left[row]
115+
Vs_right[col] = V_right
116+
end
117+
end
118+
119+
Vspaces[site + 1] = Vs_right
120+
end
121+
122+
# instantiate tensors
123+
Ws = map(enumerate(W_mats)) do (site, W_mat)
124+
W = jordanmpotensortype(S, T)(undef,
125+
Vspaces[site] ⊗ Pspaces[site] ←
126+
Pspaces[site] ⊗ Vspaces[site + 1])
127+
for (I, v) in enumerate(W_mat)
128+
ismissing(v) && continue
129+
if v isa MPOTensor
130+
W[I] = v
131+
elseif !iszero(v)
132+
τ = BraidingTensor{T}(eachspace(W)[I])
133+
W[I] = isone(v) ? τ : τ * v
134+
end
135+
end
136+
return W
137+
end
138+
139+
return FiniteMPOHamiltonian(Ws)
140+
end
141+
142+
"""
143+
InfiniteMPOHamiltonian(Ws::Vector{<:Matrix})
144+
145+
Create a `InfiniteMPOHamiltonian` from a vector of matrices, such that `Ws[i][j, k]` represents
146+
the the operator at site `i`, left level `j` and right level `k`. Here, the entries can be
147+
either `MPOTensor`, `Missing` or `Number`.
148+
"""
149+
function InfiniteMPOHamiltonian(Ws::Vector{<:Matrix})
150+
T = promote_type(_split_mpoham_types.(Ws)...)
151+
return InfiniteMPOHamiltonian{T}(Ws)
152+
end
153+
function InfiniteMPOHamiltonian{O}(W_mats::Vector{<:Matrix}) where {O<:MPOTensor}
154+
# InfiniteMPOHamiltonian only works for square matrices:
155+
for W_mat in W_mats
156+
size(W_mat, 1) == size(W_mat, 2) ||
157+
throw(ArgumentError("matrices should be square"))
158+
end
159+
allequal(Base.Fix2(size, 1), W_mats) ||
160+
throw(ArgumentError("matrices should have the same size"))
161+
nlvls = size(W_mats[1], 1)
162+
163+
T = scalartype(O)
164+
L = length(W_mats)
165+
# initialize sumspaces
166+
S = spacetype(O)
167+
168+
# physical spaces
169+
Pspaces = map(W_mats) do W_mat
170+
operator_id = findfirst(x -> x isa O, W_mat)
171+
@assert !isnothing(operator_id) "could not determine physical space"
172+
return physicalspace(W_mat[operator_id])
173+
end
174+
175+
# virtual spaces:
176+
# note that we assume that the FSA does not contain "dead ends", as this would mess with the
177+
# ability to deduce spaces.
178+
# also assume spacecheck errors will happen when filling the BlockTensors
179+
MissingS = Union{Missing,S}
180+
Vspaces = PeriodicArray([Vector{MissingS}(missing, nlvls) for _ in 1:L])
181+
for V in Vspaces
182+
V[1] = V[end] = oneunit(S)
183+
end
184+
185+
haschanged = true
186+
while haschanged
187+
haschanged = false
188+
# sweep left-to-right-to-left
189+
for site in vcat(1:length(W_mats), reverse(1:(length(W_mats) - 1)))
190+
W_mat = W_mats[site]
191+
Vs_left = Vspaces[site]
192+
Vs_right = Vspaces[site + 1]
193+
194+
for I in eachindex(IndexCartesian(), W_mat)
195+
Welem = W_mat[I]
196+
ismissing(Welem) && continue
197+
row, col = I.I
198+
if Welem isa MPOTensor
199+
V_left = left_virtualspace(Welem)
200+
if ismissing(Vs_left[row])
201+
Vs_left[row] = V_left
202+
haschanged = true
203+
else
204+
@assert Vs_left[row] == V_left "incompatible space between sites $(site-1) and $site at level $row"
205+
end
206+
207+
V_right = right_virtualspace(Welem)
208+
if ismissing(Vs_right[col])
209+
Vs_right[col] = V_right
210+
haschanged = true
211+
else
212+
@assert Vs_right[col] == V_right "incompatible space between sites $(site) and $(site+1) at level $col"
213+
end
214+
elseif !iszero(Welem) # Welem isa Number
215+
if ismissing(Vs_left[row]) && !ismissing(Vs_right[col])
216+
Vs_left[row] = Vs_right[col]
217+
haschanged = true
218+
elseif !ismissing(Vs_left[row]) && ismissing(Vs_right[col])
219+
Vs_right[col] = Vs_left[row]
220+
haschanged = true
221+
else
222+
@assert Vs_left[row] == Vs_right[col] "incompatible space between sites $(site-1) and $site at level $row"
223+
end
224+
end
225+
end
226+
227+
Vspaces[site] = Vs_left
228+
Vspaces[site + 1] = Vs_right
229+
end
230+
end
231+
232+
foreach(Base.Fix2(replace!, missing => oneunit(S)), Vspaces)
233+
Vsumspaces = map(Vspaces) do V
234+
return SumSpace(collect(S, V))
235+
end
236+
237+
# instantiate tensors
238+
Ws = map(enumerate(W_mats)) do (site, W_mat)
239+
W = jordanmpotensortype(S, T)(undef,
240+
Vsumspaces[site] ⊗ Pspaces[site] ←
241+
Pspaces[site] ⊗ Vsumspaces[site + 1])
242+
for (I, v) in enumerate(W_mat)
243+
ismissing(v) && continue
244+
if v isa MPOTensor
245+
W[I] = v
246+
elseif !iszero(v)
247+
τ = BraidingTensor{T}(eachspace(W)[I])
248+
W[I] = isone(v) ? τ : τ * v
249+
end
250+
end
251+
return W
252+
end
253+
254+
return InfiniteMPOHamiltonian(Ws)
255+
end
256+
257+
function _split_mpoham_types(W::Matrix)::Type{<:MPOTensor}
258+
# attempt to deduce from eltype -- hopefully type-stable
259+
T = eltype(W)
260+
if T <: MPOTensor
261+
return T
262+
elseif T <: Union{Missing,Number,MPOTensor}
263+
Ts = collect(DataType, Base.uniontypes(T))
264+
# find MPO type
265+
iTO = findall(x -> x <: MPOTensor, Ts)
266+
@assert !isempty(iTO) "should not happen"
267+
TO = promote_type(Ts[iTO]...)
268+
# check scalar type
269+
iTE = findall(x -> x <: Number, Ts)
270+
if !isempty(iTE)
271+
all(i -> Ts[i] <: scalartype(TO), iTE) ||
272+
throw(ArgumentError("scalar type should be a subtype of the tensor scalar type"))
273+
end
274+
return TO
275+
end
276+
277+
# didn't work, so we check all types
278+
TO = Base.Bottom # mpotensor type
279+
TE = Base.Bottom # scalar type
280+
for x in W
281+
Tx = typeof(x)
282+
if Tx <: MPOTensor
283+
TO = promote_type(TO, Tx)
284+
elseif Tx <: Number
285+
TE = promote_type(TE, Tx)
286+
else
287+
Tx === Missing || throw(ArgumentError("invalid type $Tx in matrix"))
288+
end
289+
end
290+
TO === Base.Bottom && throw(ArgumentError("no MPOTensor found in matrix"))
291+
TE <: scalartype(TO) ||
292+
throw(ArgumentError("scalar type should be a subtype of the tensor scalar type"))
293+
294+
return TO
295+
end
58296

59297
"""
60298
instantiate_operator(lattice::AbstractArray{<:VectorSpace}, O::Pair)

‎src/states/finitemps.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -132,7 +132,7 @@ _not_missing_type(::Type{Missing}) = throw(ArgumentError("Only missing type pres
132132
function _not_missing_type(::Type{T}) where {T}
133133
if T isa Union
134134
return (!(T.a === Missing) && !(T.b === Missing)) ? T :
135-
!(T.a === Missing) ? _not_missing_type(T.b) : _not_missing_type(T.b)
135+
!(T.a === Missing) ? _not_missing_type(T.a) : _not_missing_type(T.b)
136136
else
137137
return T
138138
end

‎test/operators.jl‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,46 @@ vspaces = (ℙ^10, Rep[U₁]((0 => 20)), Rep[SU₂](1 // 2 => 10, 3 // 2 => 5, 5
7171
end
7272
end
7373

74+
@testset "MPOHamiltonian constructors" begin
75+
P = ℂ^2
76+
T = Float64
77+
78+
H1 = randn(T, P ← P)
79+
H1 += H1'
80+
D = FiniteMPO(H1)[1]
81+
82+
H2 = randn(T, P^2 ← P^2)
83+
H2 += H2'
84+
C, B = FiniteMPO(H2)[1:2]
85+
86+
Elt = Union{Missing,typeof(D),scalartype(D)}
87+
Wmid = Elt[1.0 C D; 0.0 0.0 B; 0.0 0.0 1.0]
88+
Wleft = Wmid[1:1, :]
89+
Wright = Wmid[:, end:end]
90+
91+
# Finite
92+
Ws = [Wleft, Wmid, Wmid, Wright]
93+
H = FiniteMPOHamiltonian(fill(P, 4),
94+
[(i,) => H1 for i in 1:4]...,
95+
[(i, i + 1) => H2 for i in 1:3]...)
96+
H′ = FiniteMPOHamiltonian(Ws)
97+
@test H ≈ H′
98+
99+
H′ = FiniteMPOHamiltonian(map(Base.Fix1(collect, Any), Ws)) # without type info
100+
@test H ≈ H′
101+
102+
# Infinite
103+
Ws = [Wmid]
104+
H = InfiniteMPOHamiltonian(fill(P, 1),
105+
[(i,) => H1 for i in 1:1]...,
106+
[(i, i + 1) => H2 for i in 1:1]...)
107+
H′ = InfiniteMPOHamiltonian(Ws)
108+
@test all(parent(H) .≈ parent(H′))
109+
110+
H′ = InfiniteMPOHamiltonian(map(Base.Fix1(collect, Any), Ws)) # without type info
111+
@test all(parent(H) .≈ parent(H′))
112+
end
113+
74114
@testset "Finite MPOHamiltonian" begin
75115
L = 3
76116
T = ComplexF64

0 commit comments

Comments
 (0)