From 8150a3be0508dc115f19f35a8f4f546851ca655c Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Wed, 29 Mar 2023 15:27:43 -0600 Subject: [PATCH 01/14] Import Loic MpoMatrix->InfiniteMPO code Setup unit tests to confirm InfiniteSum{MPO} == l*InfiniteMPO*r (l* *r means terminated) --- src/infinitempo.jl | 2 + src/infinitempomatrix.jl | 382 ++++++++++++++++++++++++++++++++++- src/models/models.jl | 8 + test/test_iMPOConversions.jl | 115 +++++++++++ 4 files changed, 506 insertions(+), 1 deletion(-) create mode 100644 test/test_iMPOConversions.jl diff --git a/src/infinitempo.jl b/src/infinitempo.jl index 1c2c75f..db9c44f 100644 --- a/src/infinitempo.jl +++ b/src/infinitempo.jl @@ -11,3 +11,5 @@ mutable struct InfiniteMPO <: AbstractInfiniteMPS end translator(mpo::InfiniteMPO) = mpo.data.translator + +InfiniteMPO(data::CelledVector{ITensor}) = InfiniteMPO(data, 0, size(data, 1), false) diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index 32cd2f6..d72c501 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -43,4 +43,384 @@ function InfiniteMPOMatrix(data::Vector{Matrix{ITensor}}, translator::Function) return InfiniteMPOMatrix(CelledVector(data, translator), 0, size(data)[1], false) end -#nrange(H::InfiniteMPOMatrix) = [size(H[j])[1] - 1 for j in 1:nsites(H)] +function matrixITensorToITensor( + H::Matrix{ITensor}, com_inds, left_dir, right_dir; kwargs... +) + init_all = get(kwargs, :init_all, true) + init_left_first = get(kwargs, :init_left_first, init_all) + init_right_first = get(kwargs, :init_right_first, init_all) + init_left_last = get(kwargs, :init_left_last, init_all) + init_right_last = get(kwargs, :init_right_last, init_all) + #TODO: fix the rev + rev_all = get(kwargs, :rev_all, false) + rev_left = get(kwargs, :rev_left, rev_all) + rev_right = get(kwargs, :rev_right, rev_all) + #TODO some fixing when rev and init are mixed up + + lx, ly = size(H) + #Generate in order the leftbasis + left_basis = valtype(com_inds)[] + for j in 1:lx + for k in 1:ly + temp_ind = filter(x -> dir(x) == left_dir, uniqueinds(H[j, k], com_inds)) + if length(temp_ind) == 1 + append!(left_basis, temp_ind) + break + end + end + end + if init_left_first + left_basis = [ + Index( + QN() => 1; + dir=left_dir, + tags=length(left_basis) > 0 ? tags(left_basis[1]) : "left_link", + ), + left_basis..., + ] #Dummy index for the first line + end + init_left_last && append!( + left_basis, + [ + Index( + QN() => 1; + dir=left_dir, + tags=length(left_basis) > 0 ? tags(left_basis[1]) : "left_link", + ), + ], + ) #Dummy index for the last line + + right_basis = valtype(com_inds)[] + for k in 1:ly + for j in 1:lx + temp_ind = filter(x -> dir(x) == right_dir, uniqueinds(H[j, k], com_inds)) + if length(temp_ind) == 1 + append!(right_basis, temp_ind) + break + end + end + end + if init_right_first + right_basis = [ + Index( + QN() => 1; + dir=right_dir, + tags=length(right_basis) > 0 ? tags(right_basis[1]) : "right_link", + ), + right_basis..., + ] #Dummy index for the first column + end + init_right_last && append!( + right_basis, + [ + Index( + QN() => 1; + dir=right_dir, + tags=length(right_basis) > 0 ? tags(right_basis[1]) : "right_link", + ), + ], + ) #Dummy index for the last column + + left_block = Vector{Pair{QN,Int64}}() + dic_inv_left_ind = Dict{Tuple{UInt64,Int64},Int64}() + for index in left_basis + for (n, qp) in enumerate(index.space) + append!(left_block, [qp]) + dic_inv_left_ind[index.id, n] = length(left_block) + end + end + new_left_index = if length(left_basis) == 1 + left_basis[1] + else + Index(left_block; dir=left_dir, tags=tags(left_basis[1])) + end + + right_block = Vector{Pair{QN,Int64}}() + dic_inv_right_ind = Dict{Tuple{UInt64,Int64},Int64}() + for index in right_basis + for (n, qp) in enumerate(index.space) + append!(right_block, [qp]) + dic_inv_right_ind[index.id, n] = length(right_block) + end + end + new_right_index = if length(right_basis) == 1 + right_basis[1] + else + Index(right_block; dir=right_dir, tags=tags(right_basis[1])) + end + #Determine the non-zero blocks, not efficient in memory for now TODO: improve memory use + temp_block = Block{4}[] + elements = [] + dummy_left = ITensor(1, Index(QN() => 1)) + dum_left_ind = only(inds(dummy_left)) + dummy_right = ITensor(1, Index(QN() => 1)) + dum_right_ind = only(inds(dummy_right)) + for x in 1:lx + for y in 1:ly + isempty(H[x, y]) && continue + tli = filter(x -> dir(x) == dir(left_basis[1]), commoninds(H[x, y], left_basis)) + #This first part find which structure the local Ham has and ensure H[x, y] is properly ordered + case = 0 + if !isempty(tli) + li = only(tli) + tri = filter(x -> dir(x) == dir(right_basis[1]), commoninds(H[x, y], right_basis)) + if !isempty(tri) + ri = only(tri) + T = permute(H[x, y], com_inds..., li, ri; allow_alias=true) + case = 3 #This is the default case, both legs exists + else + #println(x, y) + !(y == 1 || (x == lx && y == ly)) && error("Incompatible leg") + T = permute( + H[x, y] * dummy_right, com_inds..., li, dum_right_ind; allow_alias=true + ) + case = 1 + end + else + !(x == lx || (x == 1 && y == 1)) && error("Incompatible leg") + tri = filter(x -> dir(x) == dir(right_basis[1]), commoninds(H[x, y], right_basis)) + if !isempty(tri) + ri = only(tri) + T = permute(H[x, y] * dummy_left, com_inds..., dum_left_ind, ri; allow_alias=true) + case = 2 + else + !((x==1 && y == 1) || (x==lx && y == 1) ||(x == lx && y == ly)) && error("Incompatible leg") + T = permute( + H[x, y] * dummy_left * dummy_right, + com_inds..., + dum_left_ind, + dum_right_ind; + allow_alias=true, + ) + end + end + for (n, b) in enumerate(eachnzblock(T)) + #TODO not completely ok for attribution to 1 and end when stuff is missing + norm(T[b]) == 0 && continue + if case == 0 + if x == 1 && y == 1 + append!( + temp_block, + [ + Block( + b[1], + b[2], + dic_inv_left_ind[left_basis[1].id, 1], + dic_inv_right_ind[right_basis[1].id, 1], + ), + ], + ) + elseif x==lx && y == 1 + append!( + temp_block, + [ + Block( + b[1], + b[2], + dic_inv_left_ind[left_basis[end].id, 1], + dic_inv_right_ind[right_basis[1].id, 1], + ), + ], + ) + elseif x == lx && y == ly + append!( + temp_block, + [ + Block( + b[1], + b[2], + dic_inv_left_ind[left_basis[end].id, 1], + dic_inv_right_ind[right_basis[end].id, 1], + ), + ], + ) + else + error("Something went wrong at $((x, y))") + end + elseif case == 1 + append!( + temp_block, + [ + Block( + b[1], + b[2], + dic_inv_left_ind[li.id, b[3]], + dic_inv_right_ind[right_basis[1].id, 1], + ), + ], + ) + elseif case == 2 + append!( + temp_block, + [ + Block( + b[1], + b[2], + dic_inv_left_ind[left_basis[end].id, 1], + dic_inv_right_ind[ri.id, b[4]], + ), + ], + ) + elseif case == 3 #Default case + append!( + temp_block, + [ + Block( + b[1], b[2], dic_inv_left_ind[li.id, b[3]], dic_inv_right_ind[ri.id, b[4]] + ), + ], + ) + else + println("Not treated case") + end + append!(elements, [T[b]]) + end + end + end + Hf = ITensors.BlockSparseTensor( + eltype(elements[1]), undef, temp_block, (com_inds..., new_left_index, new_right_index) + ) + for (n, b) in enumerate(temp_block) + ITensors.blockview(Hf, b) .= elements[n] + end + return itensor(Hf), new_left_index, new_right_index +end + +function matrixITensorToITensor(H::Vector{ITensor}, com_inds; rev=false, kwargs...) + rev && error("not yet implemented") + init_all = get(kwargs, :init_all, true) + init_first = get(kwargs, :init_first, init_all) + init_last = get(kwargs, :init_last, init_all) + + lx = length(H) + #Generate in order the leftbasis + left_basis = valtype(com_inds)[] #Dummy index for the first line + for j in 1:lx + append!(left_basis, uniqueinds(H[j], com_inds)) + end + left_dir = dir(left_basis[1]) + if init_first + left_basis = [Index(QN() => 1; dir=left_dir), left_basis...] + end + init_last && append!(left_basis, [Index(QN() => 1; dir=left_dir)]) #Dummy index for the last line + + left_block = Vector{Pair{QN,Int64}}() + dic_inv_left_ind = Dict{Tuple{UInt64,Int64},Int64}() + for index in left_basis + for (n, qp) in enumerate(index.space) + append!(left_block, [qp]) + dic_inv_left_ind[index.id, n] = length(left_block) + end + end + new_left_index = Index(left_block; dir=left_dir, tags="left_link") + + #Determine the non-zero blocks, not efficient in memory for now TODO: improve memory use + temp_block = Block{3}[] + elements = [] + for x in 1:lx + isempty(H[x]) && continue + tli = commoninds(H[x], left_basis) + #This first part find which structure the local Ham has and ensure H[x, y] is properly ordered + case = 0 + if !isempty(tli) + li = only(tli) + T = permute(H[x], com_inds..., li; allow_alias=true) + case = 1 #This is the default case, the leg exists + else + !(x == lx || x == 1) && error("Incompatible leg") + T = permute(H[x], com_inds...; allow_alias=true) + end + for (n, b) in enumerate(eachnzblock(T)) + norm(T[b]) == 0 && continue + if case == 0 + if x == 1 + append!(temp_block, [Block(b[1], b[2], dic_inv_left_ind[left_basis[1].id, 1])]) + elseif x == lx + append!(temp_block, [Block(b[1], b[2], dic_inv_left_ind[left_basis[end].id, 1])]) + else + error("Something went wrong") + end + elseif case == 1 #Default case + append!(temp_block, [Block(b[1], b[2], dic_inv_left_ind[li.id, b[3]])]) + else + println("Not treated case") + end + append!(elements, [T[b]]) + end + end + Hf = ITensors.BlockSparseTensor( + eltype(elements[1]), undef, temp_block, (com_inds..., new_left_index) + ) + for (n, b) in enumerate(temp_block) + ITensors.blockview(Hf, b) .= elements[n] + end + return itensor(Hf), new_left_index +end + +function matrixITensorToITensor(H::Matrix{ITensor}, idleft, idright; kwargs...) + lx, ly = size(H) + s = commoninds(H[1, 1], H[end, end]) + right_ind = uniqueinds(H[end, end], s) + for j in reverse(1:(ly - 1)) + length(right_ind) == 1 && break + right_ind = uniqueinds(H[end, j], s) + end + length(right_ind) != 1 && error("Not able to isolate the right index") + left_ind = uniqueinds(H[end, end], s) + for j in 2:lx + left_ind = uniqueinds(H[j, 1], s) + length(left_ind) == 1 && break + end + length(left_ind) != 1 && error("Not able to isolate the left index") + dir_left_ind = dir(only(left_ind)) + dir_right_ind = dir(only(right_ind)) + #return matrixITensorToITensor(H, s, left_tag, right_tag; dir_left_ind, dir_right_ind) + if length(collect(idleft)) == 1 || length(collect(idright)) == 1 + return matrixITensorToITensor(H[idleft, idright], s; kwargs...) + end + return matrixITensorToITensor( + H[idleft, idright], s, dir_left_ind, dir_right_ind; kwargs... + ) +end + +function matrixITensorToITensor(H::Matrix{ITensor}; kwargs...) + lx, ly = size(H) + s = commoninds(H[1, 1], H[end, end]) + right_ind = uniqueinds(H[end, end], s) + for j in reverse(1:(ly - 1)) + length(right_ind) == 1 && break + right_ind = uniqueinds(H[end, j], s) + end + length(right_ind) != 1 && error("Not able to isolate the right index") + left_ind = uniqueinds(H[end, end], s) + for j in 2:lx + left_ind = uniqueinds(H[j, 1], s) + length(left_ind) == 1 && break + end + length(left_ind) != 1 && error("Not able to isolate the left index") + dir_left_ind = dir(only(left_ind)) + dir_right_ind = dir(only(right_ind)) + #return matrixITensorToITensor(H, s, left_tag, right_tag; dir_left_ind, dir_right_ind) + return matrixITensorToITensor(H, s, dir_left_ind, dir_right_ind; kwargs...) +end + +function InfiniteMPO(H::InfiniteMPOMatrix) + temp = matrixITensorToITensor.(H) + new_H = CelledVector([x[1] for x in temp], translator(H)) + lis = CelledVector([x[2] for x in temp], translator(H)) + ris = CelledVector([x[3] for x in temp], translator(H)) + #retags the right_links + s = [commoninds(H[j][1, 1], H[j][end, end])[1] for j in 1:nsites(H)] + for j in 1:nsites(H) + newTag = "Link,c=$(getcell(s[j])),n=$(getsite(s[j]))" + temp = replacetags(ris[j], tags(ris[j]), newTag) + new_H[j] = replaceinds(new_H[j], [ris[j]], [temp]) + ris[j] = replacetags(ris[j], tags(ris[j]), newTag) + end + # joining the indexes + for j in 1:nsites(H) + temp = δ(dag(ris[j]), dag(lis[j + 1])) + new_H[j + 1] *= temp + end + return InfiniteMPO(new_H) +end diff --git a/src/models/models.jl b/src/models/models.jl index 75d5085..271c6ff 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -82,6 +82,14 @@ end import ITensors: op op(::OpName"Zero", ::SiteType, s::Index) = ITensor(s', dag(s)) +function InfiniteMPO(model::Model, s::CelledVector; kwargs...) + return InfiniteMPO(model, s, translator(s); kwargs...) +end + +function InfiniteMPO(model::Model, s::CelledVector, translator::Function; kwargs...) + InfiniteMPO(InfiniteMPOMatrix(model,s,translator;kwargs...)) +end + function InfiniteMPOMatrix(model::Model, s::CelledVector; kwargs...) return InfiniteMPOMatrix(model, s, translator(s); kwargs...) end diff --git a/test/test_iMPOConversions.jl b/test/test_iMPOConversions.jl new file mode 100644 index 0000000..7748424 --- /dev/null +++ b/test/test_iMPOConversions.jl @@ -0,0 +1,115 @@ +using ITensors +using ITensorInfiniteMPS +using Test + +# +# InfiniteMPO has dangling links at the end of the chain. We contract these on the outside +# with l,r terminating vectors, to make a finite lattice MPO. +# +function terminate(h::InfiniteMPO)::MPO + Ncell=nsites(h) + # left termination vector + il1=commonind(h[1],h[2]) + il0,=noncommoninds(h[1],il1,tags="Link") + l=ITensor(0.0,il0) + l[il0=>dim(il0)]=1.0 #assuming lower reg form in h + # right termination vector + iln=commonind(h[Ncell-1],h[Ncell]) + ilnp,=noncommoninds(h[Ncell],iln,tags="Link") + r=ITensor(0.0,ilnp) + r[ilnp=>1]=1.0 #assuming lower reg form in h + # build up a finite MPO + hf=MPO(Ncell) + hf[1]=dag(l)*h[1] #left terminate + hf[Ncell]=h[Ncell]*dag(r) #right terminate + for n in 2:Ncell-1 + hf[n]=h[n] #fill in the bulk. + end + return hf +end +# +# Terminate and then call expect +# for inf ψ and finite h, which is already supported in src/infinitecanonicalmps.jl +# +function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteMPO) + return expect(ψ,terminate(h)) #defer to src/infinitecanonicalmps.jl +end + +#H = ΣⱼΣn (½ S⁺ⱼS⁻ⱼ₊n + ½ S⁻ⱼS⁺ⱼ₊n + SᶻⱼSᶻⱼ₊n) +function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenbergNNN";NNN::Int64) + opsum = OpSum() + for n=1:NNN + J=1.0/n + opsum += J*0.5, "S+", 1, "S-", 1+n + opsum += J*0.5, "S-", 1, "S+", 1+n + opsum += J, "Sz", 1, "Sz", 1+n + end + return opsum +end + +function ITensorInfiniteMPS.unit_cell_terms(::Model"hubbardNNN";NNN::Int64) + U::Float64=0.25 + t::Float64=1.0 + V::Float64=0.5 + opsum = OpSum() + opsum += (U, "Nupdn", 1) + for n=1:NNN + tj,Vj=t/n,V/n + opsum += -tj, "Cdagup", 1, "Cup", 1 + n + opsum += -tj, "Cdagup", 1 + n, "Cup", 1 + opsum += -tj, "Cdagdn", 1, "Cdn", 1 + n + opsum += -tj, "Cdagdn", 1 + n, "Cdn", 1 + opsum += Vj, "Ntot" , 1, "Ntot", 1 + n + end + return opsum +end +# +# InfiniteMPO has dangling links at the end of the chain. We contract these outside +# with l,r terminating vectors, to make a finite lattice MPO. And then call expect +# for inf ψ and finite h, whihc is alread supported in src/infinitecanonicalmps.jl +# +function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteMPO) + Ncell=nsites(h) + # left termination vector + il1=commonind(h[1],h[2]) + il0,=noncommoninds(h[1],il1,tags="Link") + l=ITensor(0.0,il0) + l[il0=>dim(il0)]=1.0 #assuming lower reg form in h + # right termination vector + iln=commonind(h[Ncell-1],h[Ncell]) + ilnp,=noncommoninds(h[Ncell],iln,tags="Link") + r=ITensor(0.0,ilnp) + r[ilnp=>1]=1.0 #assuming lower reg form in h + # build up a finite MPO + hf=MPO(Ncell) + hf[1]=dag(l)*h[1] #left terminate + hf[Ncell]=h[Ncell]*dag(r) #right terminate + for n in 2:Ncell-1 + hf[n]=h[n] #fill in the bulk. + end + return expect(ψ,hf) #defer to src/infinitecanonicalmps.jl +end + +@testset verbose=true "InfiniteMPOMatrix -> InfiniteMPO" begin + ferro(n) = "↑" + antiferro(n) = isodd(n) ? "↑" : "↓" + + models=[(Model"heisenbergNNN"(),"S=1/2"),(Model"hubbardNNN"(),"Electron")] + @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af" for (model,site) in models, Ncell in 2:6, NNN in 1:Ncell-1, Af in [true,false] + if isodd(Ncell) && Af #skip test since Af state does fit inside odd cells. + continue + end + initstate(n) = Af ? antiferro(n) : ferro(n) + model_kwargs=(NNN=NNN,) + s = infsiteinds(site, Ncell; initstate,conserve_qns=true); + ψ = InfMPS(s, initstate) + Hi=InfiniteMPO(model, s;model_kwargs...) + Hs=InfiniteSum{MPO}(model, s;model_kwargs...) + Es=expect(ψ,Hs) + Ei=expect(ψ,Hi) + #@show Es Ei + @test sum(Es[1:Ncell-NNN]) ≈ Ei atol=1e-14 + end + +end +nothing \ No newline at end of file From a3524aefd8e272f40a758df30efefa6d51521e50 Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Thu, 30 Mar 2023 13:06:44 -0600 Subject: [PATCH 02/14] Support dense conversions. --- src/ITensorInfiniteMPS.jl | 1 + src/infinitempomatrix.jl | 453 ++++++++++++++++++++++------------- src/subtensor.jl | 400 +++++++++++++++++++++++++++++++ test/test_iMPOConversions.jl | 12 +- 4 files changed, 692 insertions(+), 174 deletions(-) create mode 100644 src/subtensor.jl diff --git a/src/ITensorInfiniteMPS.jl b/src/ITensorInfiniteMPS.jl index d0ef847..0d9eef0 100644 --- a/src/ITensorInfiniteMPS.jl +++ b/src/ITensorInfiniteMPS.jl @@ -52,6 +52,7 @@ include("vumps_generic.jl") include("vumps_localham.jl") include("vumps_nonlocalham.jl") include("vumps_mpo.jl") +include("subtensor.jl") export Cell, CelledVector, diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index d72c501..866da09 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -44,7 +44,7 @@ function InfiniteMPOMatrix(data::Vector{Matrix{ITensor}}, translator::Function) end function matrixITensorToITensor( - H::Matrix{ITensor}, com_inds, left_dir, right_dir; kwargs... + H::Matrix{ITensor}, com_inds, left_dir, right_dir,left_inds,right_inds; kwargs... ) init_all = get(kwargs, :init_all, true) init_left_first = get(kwargs, :init_left_first, init_all) @@ -56,23 +56,37 @@ function matrixITensorToITensor( rev_left = get(kwargs, :rev_left, rev_all) rev_right = get(kwargs, :rev_right, rev_all) #TODO some fixing when rev and init are mixed up + qns=hasqns(H[1,1]) + sample_index=com_inds[1] + tspace=ITensors.trivial_space(sample_index) + # spaceT=typeof(space(sample_index)) + # @show qns spaceT ITensors.trivial_space(sample_index) + # @show com_inds valtype(com_inds) + # indexT=valtype(com_inds) + # @show ITensors.trivial_space(sample_index) + left_inds=reverse(left_inds) + right_inds=reverse(right_inds) lx, ly = size(H) #Generate in order the leftbasis - left_basis = valtype(com_inds)[] - for j in 1:lx - for k in 1:ly - temp_ind = filter(x -> dir(x) == left_dir, uniqueinds(H[j, k], com_inds)) - if length(temp_ind) == 1 - append!(left_basis, temp_ind) - break - end - end - end + # left_basis = valtype(com_inds)[] + # for j in 1:lx + # for k in 1:ly + # temp_ind = filter(x -> dir(x) == left_dir, uniqueinds(H[j, k], com_inds)) + # if length(temp_ind) == 1 + # append!(left_basis, temp_ind) + # break + # end + # end + # end + # #@show left_basis left_inds + # @assert all(left_basis .== left_inds) + left_basis=left_inds + if init_left_first left_basis = [ Index( - QN() => 1; + tspace; dir=left_dir, tags=length(left_basis) > 0 ? tags(left_basis[1]) : "left_link", ), @@ -83,27 +97,31 @@ function matrixITensorToITensor( left_basis, [ Index( - QN() => 1; + tspace; dir=left_dir, tags=length(left_basis) > 0 ? tags(left_basis[1]) : "left_link", ), ], ) #Dummy index for the last line - right_basis = valtype(com_inds)[] - for k in 1:ly - for j in 1:lx - temp_ind = filter(x -> dir(x) == right_dir, uniqueinds(H[j, k], com_inds)) - if length(temp_ind) == 1 - append!(right_basis, temp_ind) - break - end - end - end + # right_basis = valtype(com_inds)[] + # for k in 1:ly + # for j in 1:lx + # temp_ind = filter(x -> dir(x) == right_dir, uniqueinds(H[j, k], com_inds)) + # if length(temp_ind) == 1 + # append!(right_basis, temp_ind) + # break + # end + # end + # end + # #@show right_basis right_inds + # @assert all(right_basis.==right_inds) + right_basis=right_inds + if init_right_first right_basis = [ Index( - QN() => 1; + tspace; dir=right_dir, tags=length(right_basis) > 0 ? tags(right_basis[1]) : "right_link", ), @@ -114,176 +132,232 @@ function matrixITensorToITensor( right_basis, [ Index( - QN() => 1; + tspace; dir=right_dir, tags=length(right_basis) > 0 ? tags(right_basis[1]) : "right_link", ), ], ) #Dummy index for the last column - left_block = Vector{Pair{QN,Int64}}() - dic_inv_left_ind = Dict{Tuple{UInt64,Int64},Int64}() - for index in left_basis - for (n, qp) in enumerate(index.space) - append!(left_block, [qp]) - dic_inv_left_ind[index.id, n] = length(left_block) + #Ts=typeof(left_basis[1][1]) + #@show left_basis dim.(left_basis) + if qns + left_block = Vector{Pair{QN,Int64}}() + dic_inv_left_ind = Dict{Tuple{UInt64,Int64},Int64}() + for index in left_basis + for (n, qp) in enumerate(index.space) + append!(left_block, [qp]) + dic_inv_left_ind[index.id, n] = length(left_block) + end end + else + left_block=sum(dim.(left_basis)) end new_left_index = if length(left_basis) == 1 left_basis[1] else Index(left_block; dir=left_dir, tags=tags(left_basis[1])) end + #@show new_left_index - right_block = Vector{Pair{QN,Int64}}() - dic_inv_right_ind = Dict{Tuple{UInt64,Int64},Int64}() - for index in right_basis - for (n, qp) in enumerate(index.space) - append!(right_block, [qp]) - dic_inv_right_ind[index.id, n] = length(right_block) + if qns + right_block = Vector{Pair{QN,Int64}}() + dic_inv_right_ind = Dict{Tuple{UInt64,Int64},Int64}() + for index in right_basis + for (n, qp) in enumerate(index.space) + append!(right_block, [qp]) + dic_inv_right_ind[index.id, n] = length(right_block) + end end + else + right_block=sum(dim.(right_basis)) end new_right_index = if length(right_basis) == 1 right_basis[1] else Index(right_block; dir=right_dir, tags=tags(right_basis[1])) end - #Determine the non-zero blocks, not efficient in memory for now TODO: improve memory use - temp_block = Block{4}[] - elements = [] - dummy_left = ITensor(1, Index(QN() => 1)) - dum_left_ind = only(inds(dummy_left)) - dummy_right = ITensor(1, Index(QN() => 1)) - dum_right_ind = only(inds(dummy_right)) + #@show left_basis right_basis + #@show new_left_index new_right_index + Hf1=ITensor(0.0,new_left_index,new_right_index',com_inds) + @assert length(left_basis)==lx + @assert length(right_basis)==ly + xf=1 #coordinates into Hf1 for x in 1:lx + ilb=left_basis[x] + Dx=dim(ilb) + yf=1 for y in 1:ly - isempty(H[x, y]) && continue - tli = filter(x -> dir(x) == dir(left_basis[1]), commoninds(H[x, y], left_basis)) - #This first part find which structure the local Ham has and ensure H[x, y] is properly ordered - case = 0 - if !isempty(tli) - li = only(tli) - tri = filter(x -> dir(x) == dir(right_basis[1]), commoninds(H[x, y], right_basis)) - if !isempty(tri) - ri = only(tri) - T = permute(H[x, y], com_inds..., li, ri; allow_alias=true) - case = 3 #This is the default case, both legs exists - else - #println(x, y) - !(y == 1 || (x == lx && y == ly)) && error("Incompatible leg") - T = permute( - H[x, y] * dummy_right, com_inds..., li, dum_right_ind; allow_alias=true - ) - case = 1 + W=copy(H[x,y]) + irb=right_basis[y] + Dy=dim(irb) + if !isempty(W) + if !hasind(W,ilb) + @assert Dx==1 + W*=onehot(eltype(W), ilb => 1) end - else - !(x == lx || (x == 1 && y == 1)) && error("Incompatible leg") - tri = filter(x -> dir(x) == dir(right_basis[1]), commoninds(H[x, y], right_basis)) - if !isempty(tri) - ri = only(tri) - T = permute(H[x, y] * dummy_left, com_inds..., dum_left_ind, ri; allow_alias=true) - case = 2 - else - !((x==1 && y == 1) || (x==lx && y == 1) ||(x == lx && y == ly)) && error("Incompatible leg") - T = permute( - H[x, y] * dummy_left * dummy_right, - com_inds..., - dum_left_ind, - dum_right_ind; - allow_alias=true, - ) - end - end - for (n, b) in enumerate(eachnzblock(T)) - #TODO not completely ok for attribution to 1 and end when stuff is missing - norm(T[b]) == 0 && continue - if case == 0 - if x == 1 && y == 1 - append!( - temp_block, - [ - Block( - b[1], - b[2], - dic_inv_left_ind[left_basis[1].id, 1], - dic_inv_right_ind[right_basis[1].id, 1], - ), - ], - ) - elseif x==lx && y == 1 - append!( - temp_block, - [ - Block( - b[1], - b[2], - dic_inv_left_ind[left_basis[end].id, 1], - dic_inv_right_ind[right_basis[1].id, 1], - ), - ], - ) - elseif x == lx && y == ly - append!( - temp_block, - [ - Block( - b[1], - b[2], - dic_inv_left_ind[left_basis[end].id, 1], - dic_inv_right_ind[right_basis[end].id, 1], - ), - ], - ) - else - error("Something went wrong at $((x, y))") - end - elseif case == 1 - append!( - temp_block, - [ - Block( - b[1], - b[2], - dic_inv_left_ind[li.id, b[3]], - dic_inv_right_ind[right_basis[1].id, 1], - ), - ], - ) - elseif case == 2 - append!( - temp_block, - [ - Block( - b[1], - b[2], - dic_inv_left_ind[left_basis[end].id, 1], - dic_inv_right_ind[ri.id, b[4]], - ), - ], - ) - elseif case == 3 #Default case - append!( - temp_block, - [ - Block( - b[1], b[2], dic_inv_left_ind[li.id, b[3]], dic_inv_right_ind[ri.id, b[4]] - ), - ], - ) - else - println("Not treated case") + if !hasind(W,irb) + @assert Dy==1 + W*=onehot(eltype(W), irb => 1) end - append!(elements, [T[b]]) - end - end - end - Hf = ITensors.BlockSparseTensor( - eltype(elements[1]), undef, temp_block, (com_inds..., new_left_index, new_right_index) - ) - for (n, b) in enumerate(temp_block) - ITensors.blockview(Hf, b) .= elements[n] - end - return itensor(Hf), new_left_index, new_right_index + prime!(W,irb) + #@show inds(W,tags="Link") ilb irb + @assert dim(inds(W,tags=tags(ilb),plev=0)[1])==Dx + @assert dim(inds(W,tags=tags(irb),plev=1))[1]==Dy + W=replacetags(W,tags(ilb),tags(new_left_index),plev=0) + W=replacetags(W,tags(irb),tags(new_right_index),plev=1) + Hf1[new_left_index=>xf:xf+Dx-1,new_right_index'=>yf:yf+Dy-1]=W + + end #isempty + yf+=Dy + end #for y + xf+=Dx + end #for x + + #Determine the non-zero blocks, not efficient in memory for now TODO: improve memory use + # temp_block = Block{4}[] + # elements = [] + # dummy_left = ITensor(1, Index(tspace)) + # dum_left_ind = only(inds(dummy_left)) + # dummy_right = ITensor(1, Index(tspace)) + # dum_right_ind = only(inds(dummy_right)) + # for x in 1:lx + # for y in 1:ly + # isempty(H[x, y]) && continue + # tli = filter(x -> dir(x) == dir(left_basis[1]), commoninds(H[x, y], left_basis)) + # #This first part find which structure the local Ham has and ensure H[x, y] is properly ordered + # case = 0 + # if !isempty(tli) + # li = only(tli) + # tri = filter(x -> dir(x) == dir(right_basis[1]), commoninds(H[x, y], right_basis)) + # if !isempty(tri) + # ri = only(tri) + # T = permute(H[x, y], com_inds..., li, ri; allow_alias=true) + # case = 3 #This is the default case, both legs exists + # else + # #println(x, y) + # !(y == 1 || (x == lx && y == ly)) && error("Incompatible leg") + # T = permute( + # H[x, y] * dummy_right, com_inds..., li, dum_right_ind; allow_alias=true + # ) + # case = 1 + # end + # else + # !(x == lx || (x == 1 && y == 1)) && error("Incompatible leg") + # tri = filter(x -> dir(x) == dir(right_basis[1]), commoninds(H[x, y], right_basis)) + # if !isempty(tri) + # ri = only(tri) + # T = permute(H[x, y] * dummy_left, com_inds..., dum_left_ind, ri; allow_alias=true) + # case = 2 + # else + # !((x==1 && y == 1) || (x==lx && y == 1) ||(x == lx && y == ly)) && error("Incompatible leg") + # T = permute( + # H[x, y] * dummy_left * dummy_right, + # com_inds..., + # dum_left_ind, + # dum_right_ind; + # allow_alias=true, + # ) + # end + # end + # for (n, b) in enumerate(eachnzblock(T)) + # #TODO not completely ok for attribution to 1 and end when stuff is missing + # norm(T[b]) == 0 && continue + # if case == 0 + # if x == 1 && y == 1 + # append!( + # temp_block, + # [ + # Block( + # b[1], + # b[2], + # dic_inv_left_ind[left_basis[1].id, 1], + # dic_inv_right_ind[right_basis[1].id, 1], + # ), + # ], + # ) + # elseif x==lx && y == 1 + # append!( + # temp_block, + # [ + # Block( + # b[1], + # b[2], + # dic_inv_left_ind[left_basis[end].id, 1], + # dic_inv_right_ind[right_basis[1].id, 1], + # ), + # ], + # ) + # elseif x == lx && y == ly + # append!( + # temp_block, + # [ + # Block( + # b[1], + # b[2], + # dic_inv_left_ind[left_basis[end].id, 1], + # dic_inv_right_ind[right_basis[end].id, 1], + # ), + # ], + # ) + # else + # error("Something went wrong at $((x, y))") + # end + # elseif case == 1 + # append!( + # temp_block, + # [ + # Block( + # b[1], + # b[2], + # dic_inv_left_ind[li.id, b[3]], + # dic_inv_right_ind[right_basis[1].id, 1], + # ), + # ], + # ) + # elseif case == 2 + # append!( + # temp_block, + # [ + # Block( + # b[1], + # b[2], + # dic_inv_left_ind[left_basis[end].id, 1], + # dic_inv_right_ind[ri.id, b[4]], + # ), + # ], + # ) + # elseif case == 3 #Default case + # append!( + # temp_block, + # [ + # Block( + # b[1], b[2], dic_inv_left_ind[li.id, b[3]], dic_inv_right_ind[ri.id, b[4]] + # ), + # ], + # ) + # else + # println("Not treated case") + # end + # append!(elements, [T[b]]) + # end + + + # end + # end + # Hf = ITensors.BlockSparseTensor( + # eltype(elements[1]), undef, temp_block, (com_inds..., new_left_index, new_right_index) + # ) + # for (n, b) in enumerate(temp_block) + # ITensors.blockview(Hf, b) .= elements[n] + # end + # iHf=itensor(Hf) + # #@show dense(slice(iHf,new_left_index=>1,new_right_index=>1)) + iHf1=itensor(noprime(Hf1,tags="Link")) + #@show dense(slice(iHf1,new_left_index=>1,new_right_index=>1)) + #@show Hf1[new_left_index=>1,new_right_index=>1] + #@show nzblocks(Hf) + return iHf1, new_left_index, new_right_index end function matrixITensorToITensor(H::Vector{ITensor}, com_inds; rev=false, kwargs...) @@ -383,6 +457,43 @@ function matrixITensorToITensor(H::Matrix{ITensor}, idleft, idright; kwargs...) ) end +# +# Find all the right and left pointing link indices. Index IDs don't line up from one +# tensor to the next, so we have to rely on tags to decide what to do. +# +function find_links(H::Matrix{ITensor}) + lx, ly = size(H) + @assert lx==ly + + Ms=map(n->H[lx-n+1,ly-n],1:lx-1) #Get the sub diagonal into an array. + + indexT=typeof(inds(Ms[1],tags="Link")[1]) + left_inds,right_inds=indexT[],indexT[] + @assert length(inds(Ms[1],tags="Link"))==1 + ir,=inds(Ms[1],tags="Link") + push!(right_inds,ir) + for n in 2:length(Ms) + il,=inds(Ms[n],tags=tags(ir)) + #@show ir noncommonind(Ms[n],ir,tags="Link") + ir=noncommonind(Ms[n],il,tags="Link") + push!(left_inds,il) + if !isnothing(ir) + push!(right_inds,ir) + end + end + return left_inds,right_inds +end + +# +# H should have the form below. Only link indices are shown. +# +# I 0 0 0 0 +# M--l=n 0 0 0 0 +# 0 l=n--M--l=n-1 ... 0 0 0 +# : : : : : +# 0 0 ... l=2--M--l=1 0 0 +# 0 0 ... 0 l=1--M I +# function matrixITensorToITensor(H::Matrix{ITensor}; kwargs...) lx, ly = size(H) s = commoninds(H[1, 1], H[end, end]) @@ -401,7 +512,9 @@ function matrixITensorToITensor(H::Matrix{ITensor}; kwargs...) dir_left_ind = dir(only(left_ind)) dir_right_ind = dir(only(right_ind)) #return matrixITensorToITensor(H, s, left_tag, right_tag; dir_left_ind, dir_right_ind) - return matrixITensorToITensor(H, s, dir_left_ind, dir_right_ind; kwargs...) + #@show s right_ind left_ind dir_left_ind dir_right_ind + left_inds,right_inds=find_links(H) + return matrixITensorToITensor(H, s, dir_left_ind, dir_right_ind,left_inds,right_inds; kwargs...) end function InfiniteMPO(H::InfiniteMPOMatrix) diff --git a/src/subtensor.jl b/src/subtensor.jl new file mode 100644 index 0000000..fa9edff --- /dev/null +++ b/src/subtensor.jl @@ -0,0 +1,400 @@ +import ITensors: dim, dims, DenseTensor, QNIndex, eachindval, eachval, getindex, setindex! +import NDTensors: getperm, permute, BlockDim, blockstart, blockend +import Base.range + + +struct IndexRange + index::Index + range::UnitRange{Int64} + function IndexRange(i::Index,r::UnitRange) + return new(i,r) + end +end +irPair{T}=Pair{<:Index{T},UnitRange{Int64}} where {T} +irPairU{T}=Union{irPair{T},IndexRange} where {T} +IndexRange(ir::irPair{T}) where {T} =IndexRange(ir.first,ir.second) +IndexRange(ir::IndexRange)=IndexRange(ir.index,ir.range) + +start(ir::IndexRange)=range(ir).start +range(ir::IndexRange)=ir.range +range(i::Index)=1:dim(i) +ranges(irs::Tuple) = ntuple(i -> range(irs[i]), Val(length(irs))) +indices(irs::Tuple{Vararg{IndexRange}}) = map((ir)->ir.index ,irs) +indranges(ips::Tuple{Vararg{irPairU}}) = map((ip)->IndexRange(ip) ,ips) + +dim(ir::IndexRange)=dim(range(ir)) +dim(r::UnitRange{Int64})=r.stop-r.start+1 +dims(irs::Tuple{Vararg{IndexRange}})=map((ir)->dim(ir),irs) +redim(ip::irPair{T}) where {T} = redim(IndexRange(ip)) +redim(ir::IndexRange) = redim(ir.index,dim(ir),start(ir)-1) +redim(irs::Tuple{Vararg{IndexRange}}) = map((ir)->redim(ir) ,irs) + + +eachval(ir::IndexRange) = range(ir) +eachval(irs::Tuple{Vararg{IndexRange}}) = CartesianIndices(ranges(irs)) + +eachindval(irs::Tuple{Vararg{IndexRange}}) = (indices(irs).=> Tuple(ns) for ns in eachval(irs)) + + +#-------------------------------------------------------------------------------------- +# +# NDTensor level code which distinguishes between Dense and BlockSparse storage +# +function in_range(block_start::NTuple{N, Int64},block_end::NTuple{N, Int64},rs::UnitRange{Int64}...) where {N} + ret=true + for i in eachindex(rs) + if block_start[i]>rs[i].stop || block_end[i]=rs[i].start #in range? + istart=Base.max(dest_block_start[i],rs[i].start) + istop=Base.min(dest_block_end[i],rs[i].stop) + rs1[i]=istart-dest_block_start[i]+1:istop-dest_block_start[i]+1 + end + #@show rs1 + return Tuple(rs1) +end +# +# the range r will in general start at some index start(r) >1. This function +# counts how many QN blocks are between 1 and start(r) +# +function get_offset_block_count(in::QNIndex,r::UnitRange{Int64}) + # @show in r + rs=r.start-1 + if rs==0 + return 0,0 + end + @assert rs>0 + nb=0 + qns=space(in) + for n in eachindex(qns) + rs-=qns[n].second #dim of space + nb+=1 #increment block counts + if rs==0 + break + elseif rs<0 + nb-=1 + #@show "mid block slice" + #@show qns[n] rs nb + # @error("Slicing mid block is not supported yet.") + end + end + #@show nb rs + return nb,rs +end + +using StaticArrays +function get_offset_block_counts(inds::NTuple{N,IndsT},rs::UnitRange{Int64}...) where {N,IndsT} + # @show inds rs + dbs=StaticArrays.MVector{N,UInt}(undef) + shifts=StaticArrays.MVector{N,Int}(undef) + for i in eachindex(inds) + dbi,shift=get_offset_block_count(inds[i],rs[i]) + dbs[i]=dbi + shifts[i]=shift + end + return Block{N}(dbs),shifts +end + +function get_subtensor(T::BlockSparseTensor{ElT,N},new_inds,rs::UnitRange{Int64}...) where {ElT,N} + Ds = Vector{DenseTensor{ElT,N}}() + bs = Vector{Block{N}}() + dbs,=get_offset_block_counts(inds(T),rs...) + for (jj, b) in enumerate(eachnzblock(T)) + blockT = blockview(T, b) + if in_range(blockstart(T,b),blockend(T,b),rs...) + rs1=fix_ranges(blockstart(T,b),blockend(T,b),rs...) + #@show "In of range" b dims(blockT) rs rs1 NDTensors.blockstart(T,b) + push!(Ds,blockT[rs1...]) + bc=CartesianIndex(b)-CartesianIndex(dbs) + b=bc #Decrement block numbers by the number of skipped blocks. + push!(bs,b) + end + end + if length(Ds)==0 + return BlockSparseTensor(new_inds) + end + # + # JR: All attempts at building the new indices here at the NDTensors level failed. + # The only thing I could make work was to pass the new indices down from the ITensors + # level and use those. + # + #@show bs[1] inds(T) new_inds space(new_inds[1]) Block(bs[1][1]) + #bd=blockdim(new_inds[1],bs[1][1]) + #@show bs rs new_inds + T_sub = BlockSparseTensor(ElT, undef, bs, new_inds) + for ib in eachindex(Ds) + blockT_sub = bs[ib] + blockview(T_sub, blockT_sub) .= Ds[ib] + end + return T_sub +end + +function blockrange(T::Tensor{<:Number,N},tb::Block{N})::NTuple{N,UnitRange{Int64}} where {N} + bs=blockstart(T,tb) + be=blockend(T,tb) + return ntuple(i->bs[i]:be[i],N) +end + +function fix_ranges(dest_range::NTuple{N,UnitRange{Int64}},src_range::NTuple{N,UnitRange{Int64}},rs::UnitRange{Int64}...) where {N} + rs1=Vector{UnitRange{Int64}}(undef,N) + @assert length(rs)==N + #@show rs dest_range src_range + for i in eachindex(rs1) + @assert dest_range[i].start<=rs[i].stop || dest_range[i].stop>=rs[i].start #in range? + @assert src_range[i].start<=rs[i].stop || src_range[i].stop>=rs[i].start #in range? + ds_start=Base.max(dest_range[i].start,src_range[i].start) + dsrc = src_range[i].stop-src_range[i].start + istart=Base.max(ds_start,rs[i].start) -dest_range[i].start+1 + #istop=Base.min(ds_end,rs[i].stop) + #@show i ds_start dsrc istart + rs1[i]=istart:istart+dsrc + #@show rs1[i] + end + #@show rs1 + return Tuple(rs1) +end + +function set_subtensor(T::BlockSparseTensor{ElT,N},A::BlockSparseTensor{ElT,N},rs::UnitRange{Int64}...) where {ElT,N} +# @mpoc_assert nzblocks(T)==nzblocks(A) + # dbs,shifts=get_offset_block_counts(inds(T),rs...) + # dbs_ci=CartesianIndex(dbs) + insert=false + #rsa=ntuple(i->rs[i].start-dbs_ci[i]:rs[i].stop-dbs_ci[i],N) + rsa=ntuple(i->1:dim(inds(A)[i]),N) + #@show rs dbs_ci rsa + #@show inds(A) inds(T) + for ab in eachnzblock(A) + #@show CartesianIndex(dbs) CartesianIndex(ab) + blockA = blockview(A, ab) + #@show ab blockA blockstart(A,ab) blockend(A,ab) + + if in_range(blockstart(A,ab),blockend(A,ab),rsa...) + iA=blockstart(A,ab) + iT=ntuple(i->iA[i]+rs[i].start-1,N) + #iA=ntuple(i->iA[i]+dbs_ci[i],N) + #@show "inrange" + index_within_block,tb=blockindex(T,Tuple(iT)...) + blockT = blockview(T, tb) + #@show iA iT index_within_block tb blockT + if blockT==nothing + insertblock!(T,tb) + blockT = blockview(T, tb) + #@show "insert missing block" + #@show iA iT ab tb blockA blockT + insert=true + # index_within_block,tb=blockindex(T,Tuple(it)...) + # @show index_within_block tb + end + + # rs1=fix_ranges(blockrange(T,tb),blockrange(A,ab),rs...) + # @show rs1 + dA=[dims(blockA)...] + dT=[blockend(T,tb)...]-[index_within_block...]+fill(1,N) + rs1=ntuple(i->index_within_block[i]:index_within_block[i]+Base.min(dA[i],dT[i])-1,N) + #@show rs1 dA dT + if length(findall(dA.>dT))==0 + blockT[rs1...]=blockA #Dense assignment for each block + #@show blockT + else + rsa1=ntuple(i->rsa[i].start:rsa[i].start+Base.min(dA[i],dT[i])-1,N) + #@show rsa1 blockA[rsa1...] + blockT[rs1...]=blockA[rsa1...] #partial block Dense assignment for each block + @error "Incomplete bloc transfer." + @assert false + #@show blockT + end + end + # @show "-----block-------------" + end + if insert + #@show "blocks inserted" + #@assert false + end + #@show "---------------------------------" +end +function set_subtensor(T::DiagBlockSparseTensor{ElT,N},A::DiagBlockSparseTensor{ElT,N},rs::UnitRange{Int64}...) where {ElT,N} + dbs,=get_offset_block_counts(inds(T),rs...) + dbs_ci=CartesianIndex(dbs) + rsa=ntuple(i->rs[i].start-dbs_ci[i]:rs[i].stop-dbs_ci[i],N) + for ab in eachnzblock(A) + if in_range(blockstart(A,ab),blockend(A,ab),rsa...) + it=blockstart(A,ab) + it=ntuple(i->it[i]+dbs_ci[i],N) + index_within_block,tb=blockindex(T,Tuple(it)...) + blockT = blockview(T, tb) + blockA = blockview(A, ab) + rs1=fix_ranges(blockrange(T,tb),blockrange(A,ab),rs...) + #@show rs1 ab tb + + blockT[rs1...]=blockA #Diag assignment for each block + end + end +end + +function set_subtensor(T::DiagTensor{ElT},A::DiagTensor{ElT},rs::UnitRange{Int64}...) where {ElT} + if !all(y->y==rs[1],rs) + @error("set_subtensor(DiagTensor): All ranges must be the same, rs=$(rs).") + end + N=length(rs) + #only assign along the diagonal. + for i in rs[1] + is=ntuple(i1 -> i , N) + js=ntuple(j1 -> i-rs[1].start+1 , N) + T[is...]=A[js...] + end +end + +setindex!(T::BlockSparseTensor{ElT,N},A::BlockSparseTensor{ElT,N},irs::Vararg{UnitRange{Int64},N}) where {ElT,N} = set_subtensor(T,A,irs...) +setindex!(T::DiagTensor{ElT,N},A::DiagTensor{ElT,N},irs::Vararg{UnitRange{Int64},N}) where {ElT,N} = set_subtensor(T,A,irs...) +setindex!(T::DiagBlockSparseTensor{ElT,N},A::DiagBlockSparseTensor{ElT,N},irs::Vararg{UnitRange{Int64},N}) where {ElT,N} = set_subtensor(T,A,irs...) + + +#------------------------------------------------------------------------------------ +# +# ITensor level wrappers which allows us to handle the indices in a different manner +# depending on dense/block-sparse +# +function get_subtensor_wrapper(T::DenseTensor{ElT,N},new_inds,rs::UnitRange{Int64}...) where {ElT,N} + return ITensor(T[rs...],new_inds) +end + +function get_subtensor_wrapper(T::BlockSparseTensor{ElT,N},new_inds,rs::UnitRange{Int64}...) where {ElT,N} + return ITensor(get_subtensor(T,new_inds,rs...)) +end + + +function permute(indsT::T,irs::IndexRange...) where T<:(Tuple{Vararg{T, N}} where {N, T}) + ispec=indices(irs) #indices caller specified ranges for + inot=Tuple(noncommoninds(indsT,ispec)) #indices not specified by caller + isort=ispec...,inot... #all indices sorted so user specified ones are first. + isort_sub=redim(irs)...,inot... #all indices for subtensor + p=getperm(indsT, ntuple(n -> isort[n], length(isort))) + #@show p + return permute(isort_sub,p),permute((ranges(irs)...,ranges(inot)...),p) +end +# +# Use NDTensors T[3:4,1:3,1:6...] syntax to extract the subtensor. +# +function get_subtensor_ND(T::ITensor,irs::IndexRange...) + isub,rsub=permute(inds(T),irs...) #get indices and ranges for the subtensor + return get_subtensor_wrapper(tensor(T),isub,rsub...) #virtual dispatch based on Dense or BlockSparse +end + +function match_tagplev(i1::Index,i2s::Index...) + for i in eachindex(i2s) + #@show i tags(i1) tags(i2s[i]) plev(i1) plev(i2s[i]) + if tags(i1)==tags(i2s[i]) && plev(i1)==plev(i2s[i]) + return i + end + end + @error("match_tagplev: unable to find tag/plev for $i1 in Index set $i2s.") + return nothing +end + +function getperm_tagplev(s1, s2) + N = length(s1) + r = Vector{Int}(undef, N) + return map!(i -> match_tagplev(s1[i], s2...), r, 1:length(s1)) +end +# +# +# Permute is1 to be in the order of is2. +# BUT: match based on tags and plevs instread of IDs +# +function permute_tagplev(is1,is2) + length(is1) != length(is2) && throw( + ArgumentError( + "length of first index set, $(length(is1)) does not match length of second index set, $(length(is2))", + ), + ) + perm=getperm_tagplev(is1,is2) + #@show perm is1 is2 + return is1[invperm(perm)] +end + +# This operation is non trivial because we need to +# 1) Establish that the non-requested (not in irs) indices are identical (same ID) between T & A +# 2) Establish that requested (in irs) indices have the same tags and +# prime levels (dims are different so they cannot possibly have the same IDs) +# 3) Use a combination of tags/primes (for irs indices) or IDs (for non irs indices) to permut +# the indices of A prior to conversion to a Tensor. +# +function set_subtensor_ND(T::ITensor, A::ITensor,irs::IndexRange...) + ireqT=indices(irs) #indices caller requested ranges for + inotT=Tuple(noncommoninds(inds(T),ireqT)) #indices not requested by caller + ireqA=Tuple(noncommoninds(inds(A),inotT)) #Get the requested indices for a + inotA=Tuple(noncommoninds(inds(A),ireqA)) + if length(ireqT)!=length(ireqA) || inotA!=inotT + @show inotA inotT ireqT ireqA inotA!=inotT length(ireqT) length(ireqA) + @error("subtensor assign, incompatable indices\ndestination inds=$(inds(T)),\n source inds=$(inds(A)).") + @assert(false) + end + isortT=ireqT...,inotT... #all indices sorted so user specified ones are first. + p=getperm(inds(T), ntuple(n -> isortT[n], length(isortT))) # find p such that isort[p]==inds(T) + rsortT=permute((ranges(irs)...,ranges(inotT)...),p) #sorted ranges for T + ireqAp=permute_tagplev(ireqA,ireqT) #match based on tags & plev, NOT IDs since dims are different. + isortA=permute((ireqAp...,inotA...),p) #inotA is the same inotT, using inotA here for a less confusing read. + Ap=ITensors.permute(A,isortA...;allow_alias=true) + tensor(T)[rsortT...]=tensor(Ap) +end +function set_subtensor_ND(T::ITensor, v::Number,irs::IndexRange...) + ireqT=indices(irs) #indices caller requested ranges for + inotT=Tuple(noncommoninds(inds(T),ireqT)) #indices not requested by caller + isortT=ireqT...,inotT... #all indices sorted so user specified ones are first. + p=getperm(inds(T), ntuple(n -> isortT[n], length(isortT))) # find p such that isort[p]==inds(T) + rsortT=permute((ranges(irs)...,ranges(inotT)...),p) #sorted ranges for T + tensor(T)[rsortT...].=v +end + +getindex(T::ITensor, irs::Vararg{IndexRange,N}) where {N} = get_subtensor_ND(T,irs...) +getindex(T::ITensor, irs::Vararg{irPairU,N}) where {N} = get_subtensor_ND(T,indranges(irs)...) +setindex!(T::ITensor, A::ITensor,irs::Vararg{IndexRange,N}) where {N} = set_subtensor_ND(T,A,irs...) +setindex!(T::ITensor, A::ITensor,irs::Vararg{irPairU,N}) where {N} = set_subtensor_ND(T,A,indranges(irs)...) +setindex!(T::ITensor, v::Number,irs::Vararg{IndexRange,N}) where {N} = set_subtensor_ND(T,v,irs...) +setindex!(T::ITensor, v::Number,irs::Vararg{irPairU,N}) where {N} = set_subtensor_ND(T,v,indranges(irs)...) +#= +#-------------------------------------------------------------------------------- +# +# Some overloads of similar so we can easily create new ITensors from a template. +# +function similar(T::DenseTensor{ElT},is::Index...) where {ElT} + return ITensor(DenseTensor(ElT, undef, is)) +end +function similar(T::BlockSparseTensor{ElT},is::Index...) where {ElT} + return ITensor(BlockSparseTensor(ElT, undef, nzblocks(T), is)) +end + +function similar(T::DiagTensor{ElT},is::Index...) where {ElT} + ds=[dims(is)...] + if !all(y->y==ds[1],ds) + @error("similar(DiagTensor): All indices must have the same dimensions, dims(is)=$(dims(is)).") + end + N=dim(is[1]) + return ITensor(Diag(ElT, N),is) +end +function similar(T::DiagBlockSparseTensor{ElT},is::Index...) where {ElT} + ds=[dims(is)...] + if !all(y->y==ds[1],ds) + @error("similar(DiagTensor): All indices must have the same dimensions, dims(is)=$(dims(is)).") + end + return ITensor(DiagBlockSparseTensor(ElT, undef, nzblocks(T), is)) +end + + +function similar(T::ITensor,is::Index...) + return similar(tensor(T),is...) +end + =# + diff --git a/test/test_iMPOConversions.jl b/test/test_iMPOConversions.jl index 7748424..3e28d3f 100644 --- a/test/test_iMPOConversions.jl +++ b/test/test_iMPOConversions.jl @@ -94,14 +94,18 @@ end ferro(n) = "↑" antiferro(n) = isodd(n) ? "↑" : "↓" - models=[(Model"heisenbergNNN"(),"S=1/2"),(Model"hubbardNNN"(),"Electron")] - @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af" for (model,site) in models, Ncell in 2:6, NNN in 1:Ncell-1, Af in [true,false] - if isodd(Ncell) && Af #skip test since Af state does fit inside odd cells. + models=[ + (Model"heisenbergNNN"(),"S=1/2"), + (Model"hubbardNNN"(),"Electron") + ] + @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af, qns=$qns" for (model,site) in models, qns=[false,true], Ncell in 2:6, NNN in 1:Ncell-1, Af in [true,false] +# @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af" for (model,site) in models, qns in [true], Ncell in 3:3, NNN in 2:2, Af in [false] + if isodd(Ncell) && Af #skip test since Af state does fit inside odd cells. continue end initstate(n) = Af ? antiferro(n) : ferro(n) model_kwargs=(NNN=NNN,) - s = infsiteinds(site, Ncell; initstate,conserve_qns=true); + s = infsiteinds(site, Ncell; initstate,conserve_qns=qns); ψ = InfMPS(s, initstate) Hi=InfiniteMPO(model, s;model_kwargs...) Hs=InfiniteSum{MPO}(model, s;model_kwargs...) From fb6ab84150d1347776ceb4bf2dfb949e810859db Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Fri, 31 Mar 2023 08:50:46 -0600 Subject: [PATCH 03/14] Add FQHE tests. Reduce code volume. --- src/infinitempomatrix.jl | 421 +++-------------------------------- test/test_iMPOConversions.jl | 81 ++++--- 2 files changed, 86 insertions(+), 416 deletions(-) diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index 866da09..285e020 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -44,142 +44,48 @@ function InfiniteMPOMatrix(data::Vector{Matrix{ITensor}}, translator::Function) end function matrixITensorToITensor( - H::Matrix{ITensor}, com_inds, left_dir, right_dir,left_inds,right_inds; kwargs... + H::Matrix{ITensor}; kwargs... ) - init_all = get(kwargs, :init_all, true) - init_left_first = get(kwargs, :init_left_first, init_all) - init_right_first = get(kwargs, :init_right_first, init_all) - init_left_last = get(kwargs, :init_left_last, init_all) - init_right_last = get(kwargs, :init_right_last, init_all) - #TODO: fix the rev - rev_all = get(kwargs, :rev_all, false) - rev_left = get(kwargs, :rev_left, rev_all) - rev_right = get(kwargs, :rev_right, rev_all) - #TODO some fixing when rev and init are mixed up + com_inds = commoninds(H[1, 1], H[end, end]) + left_inds,right_inds=find_links(H) qns=hasqns(H[1,1]) sample_index=com_inds[1] tspace=ITensors.trivial_space(sample_index) - # spaceT=typeof(space(sample_index)) - # @show qns spaceT ITensors.trivial_space(sample_index) - # @show com_inds valtype(com_inds) - # indexT=valtype(com_inds) - # @show ITensors.trivial_space(sample_index) + left_inds=reverse(left_inds) right_inds=reverse(right_inds) + @assert length(left_inds)>0 + @assert length(right_inds)>0 + + left_dir=dir(left_inds[1]) + right_dir=dir(right_inds[1]) + # Should we check that all dirs are consistent? lx, ly = size(H) - #Generate in order the leftbasis - # left_basis = valtype(com_inds)[] - # for j in 1:lx - # for k in 1:ly - # temp_ind = filter(x -> dir(x) == left_dir, uniqueinds(H[j, k], com_inds)) - # if length(temp_ind) == 1 - # append!(left_basis, temp_ind) - # break - # end - # end - # end - # #@show left_basis left_inds - # @assert all(left_basis .== left_inds) - left_basis=left_inds - - if init_left_first - left_basis = [ - Index( - tspace; - dir=left_dir, - tags=length(left_basis) > 0 ? tags(left_basis[1]) : "left_link", - ), - left_basis..., - ] #Dummy index for the first line - end - init_left_last && append!( - left_basis, - [ - Index( - tspace; - dir=left_dir, - tags=length(left_basis) > 0 ? tags(left_basis[1]) : "left_link", - ), - ], - ) #Dummy index for the last line - - # right_basis = valtype(com_inds)[] - # for k in 1:ly - # for j in 1:lx - # temp_ind = filter(x -> dir(x) == right_dir, uniqueinds(H[j, k], com_inds)) - # if length(temp_ind) == 1 - # append!(right_basis, temp_ind) - # break - # end - # end - # end - # #@show right_basis right_inds - # @assert all(right_basis.==right_inds) - right_basis=right_inds + left_basis = [ + Index(tspace;dir=left_dir,tags= tags(left_inds[1]) ), + left_inds..., + Index(tspace;dir=left_dir,tags= tags(left_inds[1]) ) + ] + + right_basis = [ + Index(tspace;dir=right_dir,tags=tags(right_inds[1]) ), + right_inds..., + Index(tspace;dir=right_dir,tags=tags(right_inds[1]) ) + ] - if init_right_first - right_basis = [ - Index( - tspace; - dir=right_dir, - tags=length(right_basis) > 0 ? tags(right_basis[1]) : "right_link", - ), - right_basis..., - ] #Dummy index for the first column - end - init_right_last && append!( - right_basis, - [ - Index( - tspace; - dir=right_dir, - tags=length(right_basis) > 0 ? tags(right_basis[1]) : "right_link", - ), - ], - ) #Dummy index for the last column - - #Ts=typeof(left_basis[1][1]) - #@show left_basis dim.(left_basis) + # if qns - left_block = Vector{Pair{QN,Int64}}() - dic_inv_left_ind = Dict{Tuple{UInt64,Int64},Int64}() - for index in left_basis - for (n, qp) in enumerate(index.space) - append!(left_block, [qp]) - dic_inv_left_ind[index.id, n] = length(left_block) - end - end + left_block = vcat(space.(left_basis)...) #like directsum() + right_block = vcat(space.(right_basis)...) else left_block=sum(dim.(left_basis)) - end - new_left_index = if length(left_basis) == 1 - left_basis[1] - else - Index(left_block; dir=left_dir, tags=tags(left_basis[1])) - end - #@show new_left_index - - if qns - right_block = Vector{Pair{QN,Int64}}() - dic_inv_right_ind = Dict{Tuple{UInt64,Int64},Int64}() - for index in right_basis - for (n, qp) in enumerate(index.space) - append!(right_block, [qp]) - dic_inv_right_ind[index.id, n] = length(right_block) - end - end - else right_block=sum(dim.(right_basis)) end - new_right_index = if length(right_basis) == 1 - right_basis[1] - else - Index(right_block; dir=right_dir, tags=tags(right_basis[1])) - end - #@show left_basis right_basis - #@show new_left_index new_right_index - Hf1=ITensor(0.0,new_left_index,new_right_index',com_inds) + new_left_index=Index(left_block; dir=left_dir, tags=tags(left_basis[1])) + new_right_index=Index(right_block; dir=right_dir, tags=tags(right_basis[1])) + + Hf=ITensor(0.0,new_left_index,new_right_index',com_inds) @assert length(left_basis)==lx @assert length(right_basis)==ly xf=1 #coordinates into Hf1 @@ -201,260 +107,19 @@ function matrixITensorToITensor( W*=onehot(eltype(W), irb => 1) end prime!(W,irb) - #@show inds(W,tags="Link") ilb irb @assert dim(inds(W,tags=tags(ilb),plev=0)[1])==Dx @assert dim(inds(W,tags=tags(irb),plev=1))[1]==Dy W=replacetags(W,tags(ilb),tags(new_left_index),plev=0) W=replacetags(W,tags(irb),tags(new_right_index),plev=1) - Hf1[new_left_index=>xf:xf+Dx-1,new_right_index'=>yf:yf+Dy-1]=W + Hf[new_left_index=>xf:xf+Dx-1,new_right_index'=>yf:yf+Dy-1]=W end #isempty yf+=Dy end #for y xf+=Dx end #for x - - #Determine the non-zero blocks, not efficient in memory for now TODO: improve memory use - # temp_block = Block{4}[] - # elements = [] - # dummy_left = ITensor(1, Index(tspace)) - # dum_left_ind = only(inds(dummy_left)) - # dummy_right = ITensor(1, Index(tspace)) - # dum_right_ind = only(inds(dummy_right)) - # for x in 1:lx - # for y in 1:ly - # isempty(H[x, y]) && continue - # tli = filter(x -> dir(x) == dir(left_basis[1]), commoninds(H[x, y], left_basis)) - # #This first part find which structure the local Ham has and ensure H[x, y] is properly ordered - # case = 0 - # if !isempty(tli) - # li = only(tli) - # tri = filter(x -> dir(x) == dir(right_basis[1]), commoninds(H[x, y], right_basis)) - # if !isempty(tri) - # ri = only(tri) - # T = permute(H[x, y], com_inds..., li, ri; allow_alias=true) - # case = 3 #This is the default case, both legs exists - # else - # #println(x, y) - # !(y == 1 || (x == lx && y == ly)) && error("Incompatible leg") - # T = permute( - # H[x, y] * dummy_right, com_inds..., li, dum_right_ind; allow_alias=true - # ) - # case = 1 - # end - # else - # !(x == lx || (x == 1 && y == 1)) && error("Incompatible leg") - # tri = filter(x -> dir(x) == dir(right_basis[1]), commoninds(H[x, y], right_basis)) - # if !isempty(tri) - # ri = only(tri) - # T = permute(H[x, y] * dummy_left, com_inds..., dum_left_ind, ri; allow_alias=true) - # case = 2 - # else - # !((x==1 && y == 1) || (x==lx && y == 1) ||(x == lx && y == ly)) && error("Incompatible leg") - # T = permute( - # H[x, y] * dummy_left * dummy_right, - # com_inds..., - # dum_left_ind, - # dum_right_ind; - # allow_alias=true, - # ) - # end - # end - # for (n, b) in enumerate(eachnzblock(T)) - # #TODO not completely ok for attribution to 1 and end when stuff is missing - # norm(T[b]) == 0 && continue - # if case == 0 - # if x == 1 && y == 1 - # append!( - # temp_block, - # [ - # Block( - # b[1], - # b[2], - # dic_inv_left_ind[left_basis[1].id, 1], - # dic_inv_right_ind[right_basis[1].id, 1], - # ), - # ], - # ) - # elseif x==lx && y == 1 - # append!( - # temp_block, - # [ - # Block( - # b[1], - # b[2], - # dic_inv_left_ind[left_basis[end].id, 1], - # dic_inv_right_ind[right_basis[1].id, 1], - # ), - # ], - # ) - # elseif x == lx && y == ly - # append!( - # temp_block, - # [ - # Block( - # b[1], - # b[2], - # dic_inv_left_ind[left_basis[end].id, 1], - # dic_inv_right_ind[right_basis[end].id, 1], - # ), - # ], - # ) - # else - # error("Something went wrong at $((x, y))") - # end - # elseif case == 1 - # append!( - # temp_block, - # [ - # Block( - # b[1], - # b[2], - # dic_inv_left_ind[li.id, b[3]], - # dic_inv_right_ind[right_basis[1].id, 1], - # ), - # ], - # ) - # elseif case == 2 - # append!( - # temp_block, - # [ - # Block( - # b[1], - # b[2], - # dic_inv_left_ind[left_basis[end].id, 1], - # dic_inv_right_ind[ri.id, b[4]], - # ), - # ], - # ) - # elseif case == 3 #Default case - # append!( - # temp_block, - # [ - # Block( - # b[1], b[2], dic_inv_left_ind[li.id, b[3]], dic_inv_right_ind[ri.id, b[4]] - # ), - # ], - # ) - # else - # println("Not treated case") - # end - # append!(elements, [T[b]]) - # end - - - # end - # end - # Hf = ITensors.BlockSparseTensor( - # eltype(elements[1]), undef, temp_block, (com_inds..., new_left_index, new_right_index) - # ) - # for (n, b) in enumerate(temp_block) - # ITensors.blockview(Hf, b) .= elements[n] - # end - # iHf=itensor(Hf) - # #@show dense(slice(iHf,new_left_index=>1,new_right_index=>1)) - iHf1=itensor(noprime(Hf1,tags="Link")) - #@show dense(slice(iHf1,new_left_index=>1,new_right_index=>1)) - #@show Hf1[new_left_index=>1,new_right_index=>1] - #@show nzblocks(Hf) - return iHf1, new_left_index, new_right_index -end - -function matrixITensorToITensor(H::Vector{ITensor}, com_inds; rev=false, kwargs...) - rev && error("not yet implemented") - init_all = get(kwargs, :init_all, true) - init_first = get(kwargs, :init_first, init_all) - init_last = get(kwargs, :init_last, init_all) - - lx = length(H) - #Generate in order the leftbasis - left_basis = valtype(com_inds)[] #Dummy index for the first line - for j in 1:lx - append!(left_basis, uniqueinds(H[j], com_inds)) - end - left_dir = dir(left_basis[1]) - if init_first - left_basis = [Index(QN() => 1; dir=left_dir), left_basis...] - end - init_last && append!(left_basis, [Index(QN() => 1; dir=left_dir)]) #Dummy index for the last line - - left_block = Vector{Pair{QN,Int64}}() - dic_inv_left_ind = Dict{Tuple{UInt64,Int64},Int64}() - for index in left_basis - for (n, qp) in enumerate(index.space) - append!(left_block, [qp]) - dic_inv_left_ind[index.id, n] = length(left_block) - end - end - new_left_index = Index(left_block; dir=left_dir, tags="left_link") - - #Determine the non-zero blocks, not efficient in memory for now TODO: improve memory use - temp_block = Block{3}[] - elements = [] - for x in 1:lx - isempty(H[x]) && continue - tli = commoninds(H[x], left_basis) - #This first part find which structure the local Ham has and ensure H[x, y] is properly ordered - case = 0 - if !isempty(tli) - li = only(tli) - T = permute(H[x], com_inds..., li; allow_alias=true) - case = 1 #This is the default case, the leg exists - else - !(x == lx || x == 1) && error("Incompatible leg") - T = permute(H[x], com_inds...; allow_alias=true) - end - for (n, b) in enumerate(eachnzblock(T)) - norm(T[b]) == 0 && continue - if case == 0 - if x == 1 - append!(temp_block, [Block(b[1], b[2], dic_inv_left_ind[left_basis[1].id, 1])]) - elseif x == lx - append!(temp_block, [Block(b[1], b[2], dic_inv_left_ind[left_basis[end].id, 1])]) - else - error("Something went wrong") - end - elseif case == 1 #Default case - append!(temp_block, [Block(b[1], b[2], dic_inv_left_ind[li.id, b[3]])]) - else - println("Not treated case") - end - append!(elements, [T[b]]) - end - end - Hf = ITensors.BlockSparseTensor( - eltype(elements[1]), undef, temp_block, (com_inds..., new_left_index) - ) - for (n, b) in enumerate(temp_block) - ITensors.blockview(Hf, b) .= elements[n] - end - return itensor(Hf), new_left_index -end - -function matrixITensorToITensor(H::Matrix{ITensor}, idleft, idright; kwargs...) - lx, ly = size(H) - s = commoninds(H[1, 1], H[end, end]) - right_ind = uniqueinds(H[end, end], s) - for j in reverse(1:(ly - 1)) - length(right_ind) == 1 && break - right_ind = uniqueinds(H[end, j], s) - end - length(right_ind) != 1 && error("Not able to isolate the right index") - left_ind = uniqueinds(H[end, end], s) - for j in 2:lx - left_ind = uniqueinds(H[j, 1], s) - length(left_ind) == 1 && break - end - length(left_ind) != 1 && error("Not able to isolate the left index") - dir_left_ind = dir(only(left_ind)) - dir_right_ind = dir(only(right_ind)) - #return matrixITensorToITensor(H, s, left_tag, right_tag; dir_left_ind, dir_right_ind) - if length(collect(idleft)) == 1 || length(collect(idright)) == 1 - return matrixITensorToITensor(H[idleft, idright], s; kwargs...) - end - return matrixITensorToITensor( - H[idleft, idright], s, dir_left_ind, dir_right_ind; kwargs... - ) + Hf=noprime(Hf,tags="Link") + return itensor(Hf), new_left_index, new_right_index end # @@ -494,28 +159,6 @@ end # 0 0 ... l=2--M--l=1 0 0 # 0 0 ... 0 l=1--M I # -function matrixITensorToITensor(H::Matrix{ITensor}; kwargs...) - lx, ly = size(H) - s = commoninds(H[1, 1], H[end, end]) - right_ind = uniqueinds(H[end, end], s) - for j in reverse(1:(ly - 1)) - length(right_ind) == 1 && break - right_ind = uniqueinds(H[end, j], s) - end - length(right_ind) != 1 && error("Not able to isolate the right index") - left_ind = uniqueinds(H[end, end], s) - for j in 2:lx - left_ind = uniqueinds(H[j, 1], s) - length(left_ind) == 1 && break - end - length(left_ind) != 1 && error("Not able to isolate the left index") - dir_left_ind = dir(only(left_ind)) - dir_right_ind = dir(only(right_ind)) - #return matrixITensorToITensor(H, s, left_tag, right_tag; dir_left_ind, dir_right_ind) - #@show s right_ind left_ind dir_left_ind dir_right_ind - left_inds,right_inds=find_links(H) - return matrixITensorToITensor(H, s, dir_left_ind, dir_right_ind,left_inds,right_inds; kwargs...) -end function InfiniteMPO(H::InfiniteMPOMatrix) temp = matrixITensorToITensor.(H) @@ -525,7 +168,7 @@ function InfiniteMPO(H::InfiniteMPOMatrix) #retags the right_links s = [commoninds(H[j][1, 1], H[j][end, end])[1] for j in 1:nsites(H)] for j in 1:nsites(H) - newTag = "Link,c=$(getcell(s[j])),n=$(getsite(s[j]))" + newTag = "Link,c=$(getcell(s[j])),l=$(getsite(s[j]))" temp = replacetags(ris[j], tags(ris[j]), newTag) new_H[j] = replaceinds(new_H[j], [ris[j]], [temp]) ris[j] = replacetags(ris[j], tags(ris[j]), newTag) diff --git a/test/test_iMPOConversions.jl b/test/test_iMPOConversions.jl index 3e28d3f..56ebaff 100644 --- a/test/test_iMPOConversions.jl +++ b/test/test_iMPOConversions.jl @@ -1,7 +1,6 @@ using ITensors using ITensorInfiniteMPS using Test - # # InfiniteMPO has dangling links at the end of the chain. We contract these on the outside # with l,r terminating vectors, to make a finite lattice MPO. @@ -63,31 +62,34 @@ function ITensorInfiniteMPS.unit_cell_terms(::Model"hubbardNNN";NNN::Int64) end return opsum end -# -# InfiniteMPO has dangling links at the end of the chain. We contract these outside -# with l,r terminating vectors, to make a finite lattice MPO. And then call expect -# for inf ψ and finite h, whihc is alread supported in src/infinitecanonicalmps.jl -# -function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteMPO) - Ncell=nsites(h) - # left termination vector - il1=commonind(h[1],h[2]) - il0,=noncommoninds(h[1],il1,tags="Link") - l=ITensor(0.0,il0) - l[il0=>dim(il0)]=1.0 #assuming lower reg form in h - # right termination vector - iln=commonind(h[Ncell-1],h[Ncell]) - ilnp,=noncommoninds(h[Ncell],iln,tags="Link") - r=ITensor(0.0,ilnp) - r[ilnp=>1]=1.0 #assuming lower reg form in h - # build up a finite MPO - hf=MPO(Ncell) - hf[1]=dag(l)*h[1] #left terminate - hf[Ncell]=h[Ncell]*dag(r) #right terminate - for n in 2:Ncell-1 - hf[n]=h[n] #fill in the bulk. + +function ITensors.space(::SiteType"FermionK", pos::Int; p=1, q=1, conserve_momentum=true) + if !conserve_momentum + return [QN("Nf", -p) => 1, QN("Nf", q - p) => 1] + else + return [ + QN(("Nf", -p), ("NfMom", -p * pos)) => 1, + QN(("Nf", q - p), ("NfMom", (q - p) * pos)) => 1, + ] end - return expect(ψ,hf) #defer to src/infinitecanonicalmps.jl +end + +function fermion_momentum_translator(i::Index, n::Integer; N=6) + ts = tags(i) + translated_ts = ITensorInfiniteMPS.translatecelltags(ts, n) + new_i = replacetags(i, ts => translated_ts) + for j in 1:length(new_i.space) + ch = new_i.space[j][1][1].val + mom = new_i.space[j][1][2].val + new_i.space[j] = Pair( + QN(("Nf", ch), ("NfMom", mom + n * N * ch)), new_i.space[j][2] + ) + end + return new_i + end + +function ITensors.op!(Op::ITensor, opname::OpName, ::SiteType"FermionK", s::Index...) + return ITensors.op!(Op, opname, SiteType("Fermion"), s...) end @testset verbose=true "InfiniteMPOMatrix -> InfiniteMPO" begin @@ -96,11 +98,13 @@ end models=[ (Model"heisenbergNNN"(),"S=1/2"), - (Model"hubbardNNN"(),"Electron") + (Model"hubbardNNN"(),"Electron"), + ] @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af, qns=$qns" for (model,site) in models, qns=[false,true], Ncell in 2:6, NNN in 1:Ncell-1, Af in [true,false] # @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af" for (model,site) in models, qns in [true], Ncell in 3:3, NNN in 2:2, Af in [false] - if isodd(Ncell) && Af #skip test since Af state does fit inside odd cells. + + if isodd(Ncell) && Af #skip test since Af state does fit inside odd cells. continue end initstate(n) = Af ? antiferro(n) : ferro(n) @@ -115,5 +119,28 @@ end @test sum(Es[1:Ncell-NNN]) ≈ Ei atol=1e-14 end + @testset "FQHE Hamitonian" begin + N = 6 + model = Model"fqhe_2b_pot"() + model_params = (Vs= [1.0, 0.0, 1.0, 0.0, 0.1], Ly = 3.0, prec = 1e-8) + trf=fermion_momentum_translator + function initstate(n) + mod1(n, 3) == 1 && return 2 + return 1 + end + p = 1 + q = 3 + conserve_momentum = true + s = infsiteinds("FermionK", N; translator = trf, initstate, conserve_momentum, p, q); + ψ = InfMPS(s, initstate) + Hs=InfiniteSum{MPO}(model, s; model_params...); + Hi=InfiniteMPO(model, s,trf;model_params...) + #@show inds(Hi[1]) + Es=expect(ψ,Hs) + Ei=expect(ψ,Hi) + @show Es Ei + @test Es[1] ≈ Ei atol=1e-14 + end + end nothing \ No newline at end of file From c1ad2f97e050b4b40b71546e344330116ce17d0d Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Fri, 31 Mar 2023 14:52:35 -0600 Subject: [PATCH 04/14] Get directsum working --- src/infinitempomatrix.jl | 49 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 44 insertions(+), 5 deletions(-) diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index 285e020..7f8ba45 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -47,13 +47,13 @@ function matrixITensorToITensor( H::Matrix{ITensor}; kwargs... ) com_inds = commoninds(H[1, 1], H[end, end]) - left_inds,right_inds=find_links(H) + left_inds1,right_inds1,Ms=find_links(H) qns=hasqns(H[1,1]) sample_index=com_inds[1] tspace=ITensors.trivial_space(sample_index) - left_inds=reverse(left_inds) - right_inds=reverse(right_inds) + left_inds=reverse(left_inds1) + right_inds=reverse(right_inds1) @assert length(left_inds)>0 @assert length(right_inds)>0 @@ -119,7 +119,45 @@ function matrixITensorToITensor( xf+=Dx end #for x Hf=noprime(Hf,tags="Link") - return itensor(Hf), new_left_index, new_right_index + + #@show left_inds1 right_inds1 inds(Ms[1]) + ils,irs=left_inds1,right_inds1 + Ncell=lx + @assert length(Ms)==Ncell-1 + i0=Index(tspace;dir=left_dir,tags="Link,l=0") + in=Index(tspace;dir=right_dir,tags="Link,l=$Ncell") + T=eltype(Ms[1]) + Ms[1]*=onehot(T, i0 => 1) + Ms[Ncell-1]*=onehot(T, in => 1) + + Hf1=H[end,end]*onehot(T, left_basis[end] => 1)*onehot(T, in => 1) + ih= left_basis[end],in + ilm=ils[Ncell-2],in + Hf1,ih=directsum(Hf1=>left_basis[end],Ms[Ncell-1]=>ils[Ncell-2]) + ih=ih,in + # @show ih + #@show inds(Hf1,tags="Link") ih + for i in length(Ms)-1:-1:2 + @assert hasind(Ms[i],ils[i-1]) + @assert hasind(Ms[i],irs[i]) + ilm=ils[i-1],irs[i] + # #@show dim.(ih) dir.(ih) dim.(ilm) dir.(ilm) + Hf1,ih=directsum(Hf1=>ih,Ms[i]=>ilm,tags=["Link,l=0","Link,l=$Ncell"]) + end + ilm=i0,irs[1] + Hf1,ih=directsum(Hf1=>ih,Ms[1]=>ilm,tags=["Link,l=1","Link,l=2"]) + # @show ih i0 + M=H[1,1]*onehot(T, dag(i0) => 1)*onehot(T, ih[1] => dim(ih[1])) + # @show inds(M,tags="Link") inds(Hf1,tags="Link") right_basis[1] right_basis[2] + #@show inds(M,tags="Link") dir.(inds(M,tags="Link")) dir.(inds(Hf1,tags="Link")) + Hf1,ih1=directsum(Hf1=>ih[2],M=>i0,tags="Link,l=2") + + ih=inds(Hf1,tags="Link") + #@show ih + #@show new_left_index new_right_index ih + # #@show ih typeof(Hf1) + return Hf1,ih[1],ih[2] +# return itensor(Hf), new_left_index, new_right_index end # @@ -131,6 +169,7 @@ function find_links(H::Matrix{ITensor}) @assert lx==ly Ms=map(n->H[lx-n+1,ly-n],1:lx-1) #Get the sub diagonal into an array. + #Ms=map(n->H[n+1,n],1:lx-1) #Get the sub diagonal into an array. indexT=typeof(inds(Ms[1],tags="Link")[1]) left_inds,right_inds=indexT[],indexT[] @@ -146,7 +185,7 @@ function find_links(H::Matrix{ITensor}) push!(right_inds,ir) end end - return left_inds,right_inds + return left_inds,right_inds,Ms end # From 42bce1b3985fd008b2db14af499b4753217c2c93 Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Fri, 31 Mar 2023 17:49:11 -0600 Subject: [PATCH 05/14] Stash --- src/infinitempomatrix.jl | 153 ++++++++++++++++----------------------- 1 file changed, 61 insertions(+), 92 deletions(-) diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index 7f8ba45..bc8adb0 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -43,120 +43,89 @@ function InfiniteMPOMatrix(data::Vector{Matrix{ITensor}}, translator::Function) return InfiniteMPOMatrix(CelledVector(data, translator), 0, size(data)[1], false) end +# +# H should have the form below. Only link indices are shown. +# +# I 0 0 0 0 +# M-->l=n 0 0 0 0 +# 0 l=n-->M-->l=n-1 ... 0 0 0 +# : : : : : +# 0 0 ... l=2-->M-->l=1 0 0 +# 0 0 ... 0 l=1-->M I +# +# We need to capture the corner I's and the sub-diagonal Ms, and paste them together into on ITensor. +# In addition we need make all elements of H into order(4) tensors by adding dummy Dw=1 indices. +# function matrixITensorToITensor( - H::Matrix{ITensor}; kwargs... + Hm::Matrix{ITensor}; kwargs... ) - com_inds = commoninds(H[1, 1], H[end, end]) - left_inds1,right_inds1,Ms=find_links(H) - qns=hasqns(H[1,1]) - sample_index=com_inds[1] - tspace=ITensors.trivial_space(sample_index) - - left_inds=reverse(left_inds1) - right_inds=reverse(right_inds1) - @assert length(left_inds)>0 - @assert length(right_inds)>0 - - left_dir=dir(left_inds[1]) - right_dir=dir(right_inds[1]) - # Should we check that all dirs are consistent? - - lx, ly = size(H) - left_basis = [ - Index(tspace;dir=left_dir,tags= tags(left_inds[1]) ), - left_inds..., - Index(tspace;dir=left_dir,tags= tags(left_inds[1]) ) - ] - - right_basis = [ - Index(tspace;dir=right_dir,tags=tags(right_inds[1]) ), - right_inds..., - Index(tspace;dir=right_dir,tags=tags(right_inds[1]) ) - ] - - # - if qns - left_block = vcat(space.(left_basis)...) #like directsum() - right_block = vcat(space.(right_basis)...) - else - left_block=sum(dim.(left_basis)) - right_block=sum(dim.(right_basis)) - end - new_left_index=Index(left_block; dir=left_dir, tags=tags(left_basis[1])) - new_right_index=Index(right_block; dir=right_dir, tags=tags(right_basis[1])) + indexT=typeof(inds(Hm[1,1])[1]) + T=eltype(Hm[1,1]) - Hf=ITensor(0.0,new_left_index,new_right_index',com_inds) - @assert length(left_basis)==lx - @assert length(right_basis)==ly - xf=1 #coordinates into Hf1 - for x in 1:lx - ilb=left_basis[x] - Dx=dim(ilb) - yf=1 - for y in 1:ly - W=copy(H[x,y]) - irb=right_basis[y] - Dy=dim(irb) - if !isempty(W) - if !hasind(W,ilb) - @assert Dx==1 - W*=onehot(eltype(W), ilb => 1) - end - if !hasind(W,irb) - @assert Dy==1 - W*=onehot(eltype(W), irb => 1) - end - prime!(W,irb) - @assert dim(inds(W,tags=tags(ilb),plev=0)[1])==Dx - @assert dim(inds(W,tags=tags(irb),plev=1))[1]==Dy - W=replacetags(W,tags(ilb),tags(new_left_index),plev=0) - W=replacetags(W,tags(irb),tags(new_right_index),plev=1) - Hf[new_left_index=>xf:xf+Dx-1,new_right_index'=>yf:yf+Dy-1]=W - - end #isempty - yf+=Dy - end #for y - xf+=Dx - end #for x - Hf=noprime(Hf,tags="Link") - - #@show left_inds1 right_inds1 inds(Ms[1]) - ils,irs=left_inds1,right_inds1 - Ncell=lx - @assert length(Ms)==Ncell-1 + lx, ly = size(Hm) + @assert lx==ly + # + # Extract the sub diagonal + # + Ms=map(n->Hm[lx-n+1,ly-n],1:lx-1) #Get the sub diagonal into an array. + #Ms=map(n->H[n+1,n],1:lx-1) #Get the sub diagonal into an array, reverse order. + # + # Create list of all left and right link indices on the Ms. In orefr to support dense storage we + # can't QN directions to decide left/right. Instead we look for common tags among the neighbourinf Ms. + # + N=length(Ms) + @assert N==lx-1 + ils,irs=indexT[],indexT[] + @assert length(inds(Ms[1],tags="Link"))==1 + ir,=inds(Ms[1],tags="Link") + push!(irs,ir) + for n in 2:N + il,=inds(Ms[n],tags=tags(ir)) + ir=noncommonind(Ms[n],il,tags="Link") + push!(ils,il) + if !isnothing(ir) # ir==nothing on the last M + push!(irs,ir) + end + end + # + # Make some dummy indices + # + left_dir=dir(ils[1]) + right_dir=dir(irs[1]) + tspace=ITensors.trivial_space(ils[1]) i0=Index(tspace;dir=left_dir,tags="Link,l=0") - in=Index(tspace;dir=right_dir,tags="Link,l=$Ncell") - T=eltype(Ms[1]) + in=Index(tspace;dir=right_dir,tags="Link,l=$N") + inp=prime(dag(in)) + Ms[1]*=onehot(T, i0 => 1) - Ms[Ncell-1]*=onehot(T, in => 1) + Ms[N]*=onehot(T, in => 1) + @show inds(Ms[1]) inds(Ms[N]) - Hf1=H[end,end]*onehot(T, left_basis[end] => 1)*onehot(T, in => 1) - ih= left_basis[end],in - ilm=ils[Ncell-2],in - Hf1,ih=directsum(Hf1=>left_basis[end],Ms[Ncell-1]=>ils[Ncell-2]) + H=Hm[1,1]*onehot(T, inp => 1)*onehot(T, in => 1) #I op in the top left of Hm + H,ih=directsum(H=>inp,Ms[N]=>ils[N-1]) ih=ih,in # @show ih #@show inds(Hf1,tags="Link") ih - for i in length(Ms)-1:-1:2 + for i in N-1:-1:2 @assert hasind(Ms[i],ils[i-1]) @assert hasind(Ms[i],irs[i]) ilm=ils[i-1],irs[i] # #@show dim.(ih) dir.(ih) dim.(ilm) dir.(ilm) - Hf1,ih=directsum(Hf1=>ih,Ms[i]=>ilm,tags=["Link,l=0","Link,l=$Ncell"]) + H,ih=directsum(H=>ih,Ms[i]=>ilm,tags=["Link,l=0","Link,l=$N"]) end ilm=i0,irs[1] - Hf1,ih=directsum(Hf1=>ih,Ms[1]=>ilm,tags=["Link,l=1","Link,l=2"]) + H,ih=directsum(H=>ih,Ms[1]=>ilm,tags=["Link,l=1","Link,l=2"]) # @show ih i0 - M=H[1,1]*onehot(T, dag(i0) => 1)*onehot(T, ih[1] => dim(ih[1])) + M=Hm[N+1,N+1]*onehot(T, dag(i0) => 1)*onehot(T, ih[1] => dim(ih[1])) # @show inds(M,tags="Link") inds(Hf1,tags="Link") right_basis[1] right_basis[2] #@show inds(M,tags="Link") dir.(inds(M,tags="Link")) dir.(inds(Hf1,tags="Link")) - Hf1,ih1=directsum(Hf1=>ih[2],M=>i0,tags="Link,l=2") + H,ih1=directsum(H=>ih[2],M=>i0,tags="Link,l=2") - ih=inds(Hf1,tags="Link") + ih=inds(H,tags="Link") #@show ih #@show new_left_index new_right_index ih # #@show ih typeof(Hf1) - return Hf1,ih[1],ih[2] + return H,ih[1],ih[2] # return itensor(Hf), new_left_index, new_right_index end From 7102772c5923e88a6e8a6223a5821395e70b5819 Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Fri, 31 Mar 2023 17:58:10 -0600 Subject: [PATCH 06/14] clean up direct sum code --- src/infinitempomatrix.jl | 41 +++++++++++++--------------------------- 1 file changed, 13 insertions(+), 28 deletions(-) diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index bc8adb0..d76710d 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -61,7 +61,6 @@ function matrixITensorToITensor( ) indexT=typeof(inds(Hm[1,1])[1]) T=eltype(Hm[1,1]) - lx, ly = size(Hm) @assert lx==ly # @@ -88,45 +87,31 @@ function matrixITensorToITensor( end end # - # Make some dummy indices + # Make some dummy indices. Details of the Link indices do not matter from here on. # left_dir=dir(ils[1]) right_dir=dir(irs[1]) tspace=ITensors.trivial_space(ils[1]) - i0=Index(tspace;dir=left_dir,tags="Link,l=0") - in=Index(tspace;dir=right_dir,tags="Link,l=$N") + i0=Index(tspace;dir=left_dir,tags="Link") + in=Index(tspace;dir=right_dir,tags="Link") inp=prime(dag(in)) + ils=[i0,ils...] Ms[1]*=onehot(T, i0 => 1) Ms[N]*=onehot(T, in => 1) - @show inds(Ms[1]) inds(Ms[N]) - H=Hm[1,1]*onehot(T, inp => 1)*onehot(T, in => 1) #I op in the top left of Hm - H,ih=directsum(H=>inp,Ms[N]=>ils[N-1]) + H=Hm[1,1]*onehot(T, inp => 1)*onehot(T, in => 1) #Bootstrap with the I op in the top left of Hm + H,ih=directsum(H=>inp,Ms[N]=>ils[N]) #1-D directsum to put M[N] directly below the I op ih=ih,in - # @show ih - #@show inds(Hf1,tags="Link") ih - for i in N-1:-1:2 - @assert hasind(Ms[i],ils[i-1]) - @assert hasind(Ms[i],irs[i]) - ilm=ils[i-1],irs[i] - # #@show dim.(ih) dir.(ih) dim.(ilm) dir.(ilm) - H,ih=directsum(H=>ih,Ms[i]=>ilm,tags=["Link,l=0","Link,l=$N"]) - end - ilm=i0,irs[1] - H,ih=directsum(H=>ih,Ms[1]=>ilm,tags=["Link,l=1","Link,l=2"]) - # @show ih i0 - M=Hm[N+1,N+1]*onehot(T, dag(i0) => 1)*onehot(T, ih[1] => dim(ih[1])) - # @show inds(M,tags="Link") inds(Hf1,tags="Link") right_basis[1] right_basis[2] - #@show inds(M,tags="Link") dir.(inds(M,tags="Link")) dir.(inds(Hf1,tags="Link")) - H,ih1=directsum(H=>ih[2],M=>i0,tags="Link,l=2") + for i in N-1:-1:1 + ilm=ils[i],irs[i] + H,ih=directsum(H=>ih,Ms[i]=>ilm,tags=["Link","Link"]) #2-D directsums to place new blocks below and to the right. + end + IN=Hm[N+1,N+1]*onehot(T, dag(i0) => 1)*onehot(T, ih[1] => dim(ih[1])) #I op in the bottom left of Hm + H,_=directsum(H=>ih[2],IN=>i0,tags="Link") #1-D directsum to the put I in the bottom row, to the right of M[1] - ih=inds(H,tags="Link") - #@show ih - #@show new_left_index new_right_index ih - # #@show ih typeof(Hf1) + ih=inds(H,tags="Link") return H,ih[1],ih[2] -# return itensor(Hf), new_left_index, new_right_index end # From 83a3941c6fb42f870fbcc6f0a5f7a3fd7a31dcee Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Sat, 1 Apr 2023 08:51:44 -0600 Subject: [PATCH 07/14] Finish clean of matrixITensorToITensor --- src/infinitempomatrix.jl | 112 +++++++++++++-------------------------- 1 file changed, 36 insertions(+), 76 deletions(-) diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index d76710d..965c1c1 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -44,7 +44,7 @@ function InfiniteMPOMatrix(data::Vector{Matrix{ITensor}}, translator::Function) end # -# H should have the form below. Only link indices are shown. +# Hm should have the form below. Only link indices are shown. # # I 0 0 0 0 # M-->l=n 0 0 0 0 @@ -54,11 +54,9 @@ end # 0 0 ... 0 l=1-->M I # # We need to capture the corner I's and the sub-diagonal Ms, and paste them together into on ITensor. -# In addition we need make all elements of H into order(4) tensors by adding dummy Dw=1 indices. +# In addition we need make all elements of Hm into order(4) tensors by adding dummy Dw=1 indices. # -function matrixITensorToITensor( - Hm::Matrix{ITensor}; kwargs... -) +function matrixITensorToITensor(Hm::Matrix{ITensor}) indexT=typeof(inds(Hm[1,1])[1]) T=eltype(Hm[1,1]) lx, ly = size(Hm) @@ -68,91 +66,53 @@ function matrixITensorToITensor( # Ms=map(n->Hm[lx-n+1,ly-n],1:lx-1) #Get the sub diagonal into an array. #Ms=map(n->H[n+1,n],1:lx-1) #Get the sub diagonal into an array, reverse order. - # - # Create list of all left and right link indices on the Ms. In orefr to support dense storage we - # can't QN directions to decide left/right. Instead we look for common tags among the neighbourinf Ms. - # N=length(Ms) @assert N==lx-1 - ils,irs=indexT[],indexT[] + # + # We need a dummy link index, but it has to have the right direction. + # Ms[1] should only have one link index, and it should be to the right. + # @assert length(inds(Ms[1],tags="Link"))==1 ir,=inds(Ms[1],tags="Link") - push!(irs,ir) - for n in 2:N - il,=inds(Ms[n],tags=tags(ir)) - ir=noncommonind(Ms[n],il,tags="Link") + left_dir=dir(dag(ir)) + tspace=ITensors.trivial_space(ir) + il0=Index(tspace;dir=left_dir,tags="Link,l=0") #Dw=1 left dummy index. + # + # Convert edge Ms to order 4 ITensors using the dummy index. + # + Ms[1]*=onehot(T, il0 => 1) + Ms[N]*=onehot(T, dag(il0) => 1) #We don't need distinct index IDs for this. directsum makes all new index anyway. + # + # Create list of all left and right link indices on the Ms. In order to support dense storage we + # can't use QN directions to decide left/right. Instead we look for common tags among the neighbouring Ms. + # Right now the code assumes plev=0 for all Ms. + # + ils,irs=indexT[],indexT[] #left the right link indices. + ir=il0 #set up recursion below. ir from the previous M, should have the same tags as il on the next M. + for n in 1:N + il,=inds(Ms[n],tags=tags(ir)) #Find the left index + ir=noncommonind(Ms[n],il,tags="Link") #New right index by elimination push!(ils,il) - if !isnothing(ir) # ir==nothing on the last M - push!(irs,ir) - end + push!(irs,ir) end # - # Make some dummy indices. Details of the Link indices do not matter from here on. + # Now direct sum everything non-empty in Hm, including those I ops in the corners. + # Order is critical. + # Details of the Link index tags like "l=n" do not matter from here on. # - left_dir=dir(ils[1]) - right_dir=dir(irs[1]) - tspace=ITensors.trivial_space(ils[1]) - i0=Index(tspace;dir=left_dir,tags="Link") - in=Index(tspace;dir=right_dir,tags="Link") - inp=prime(dag(in)) - ils=[i0,ils...] - - Ms[1]*=onehot(T, i0 => 1) - Ms[N]*=onehot(T, in => 1) - - H=Hm[1,1]*onehot(T, inp => 1)*onehot(T, in => 1) #Bootstrap with the I op in the top left of Hm - H,ih=directsum(H=>inp,Ms[N]=>ils[N]) #1-D directsum to put M[N] directly below the I op - ih=ih,in - for i in N-1:-1:1 - ilm=ils[i],irs[i] - H,ih=directsum(H=>ih,Ms[i]=>ilm,tags=["Link","Link"]) #2-D directsums to place new blocks below and to the right. + H=Hm[1,1]*onehot(T, il0' => 1)*onehot(T, dag(il0) => 1) #Bootstrap with the I op in the top left of Hm + H,ih=directsum(H=>il0',Ms[N]=>ils[N]) #1-D directsum to put M[N] directly below the I op + ih=ih,dag(il0) #setup recusion. + for i in N-1:-1:1 #2-D directsums to place new blocks below and to the right. + H,ih=directsum(H=>ih,Ms[i]=>(ils[i],irs[i]),tags=["Link","Link"]) end - IN=Hm[N+1,N+1]*onehot(T, dag(i0) => 1)*onehot(T, ih[1] => dim(ih[1])) #I op in the bottom left of Hm - H,_=directsum(H=>ih[2],IN=>i0,tags="Link") #1-D directsum to the put I in the bottom row, to the right of M[1] + IN=Hm[N+1,N+1]*onehot(T, dag(il0) => 1)*onehot(T, ih[1] => dim(ih[1])) #I op in the bottom left of Hm + H,_=directsum(H=>ih[2],IN=>il0,tags="Link") #1-D directsum to the put I in the bottom row, to the right of M[1] ih=inds(H,tags="Link") return H,ih[1],ih[2] end -# -# Find all the right and left pointing link indices. Index IDs don't line up from one -# tensor to the next, so we have to rely on tags to decide what to do. -# -function find_links(H::Matrix{ITensor}) - lx, ly = size(H) - @assert lx==ly - - Ms=map(n->H[lx-n+1,ly-n],1:lx-1) #Get the sub diagonal into an array. - #Ms=map(n->H[n+1,n],1:lx-1) #Get the sub diagonal into an array. - - indexT=typeof(inds(Ms[1],tags="Link")[1]) - left_inds,right_inds=indexT[],indexT[] - @assert length(inds(Ms[1],tags="Link"))==1 - ir,=inds(Ms[1],tags="Link") - push!(right_inds,ir) - for n in 2:length(Ms) - il,=inds(Ms[n],tags=tags(ir)) - #@show ir noncommonind(Ms[n],ir,tags="Link") - ir=noncommonind(Ms[n],il,tags="Link") - push!(left_inds,il) - if !isnothing(ir) - push!(right_inds,ir) - end - end - return left_inds,right_inds,Ms -end - -# -# H should have the form below. Only link indices are shown. -# -# I 0 0 0 0 -# M--l=n 0 0 0 0 -# 0 l=n--M--l=n-1 ... 0 0 0 -# : : : : : -# 0 0 ... l=2--M--l=1 0 0 -# 0 0 ... 0 l=1--M I -# - function InfiniteMPO(H::InfiniteMPOMatrix) temp = matrixITensorToITensor.(H) new_H = CelledVector([x[1] for x in temp], translator(H)) From af4b4e290cf4e9f58ff154e701be8bce4d57ecb3 Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Sat, 1 Apr 2023 09:23:39 -0600 Subject: [PATCH 08/14] Git rid of subtensor code --- src/subtensor.jl | 400 ----------------------------------------------- 1 file changed, 400 deletions(-) delete mode 100644 src/subtensor.jl diff --git a/src/subtensor.jl b/src/subtensor.jl deleted file mode 100644 index fa9edff..0000000 --- a/src/subtensor.jl +++ /dev/null @@ -1,400 +0,0 @@ -import ITensors: dim, dims, DenseTensor, QNIndex, eachindval, eachval, getindex, setindex! -import NDTensors: getperm, permute, BlockDim, blockstart, blockend -import Base.range - - -struct IndexRange - index::Index - range::UnitRange{Int64} - function IndexRange(i::Index,r::UnitRange) - return new(i,r) - end -end -irPair{T}=Pair{<:Index{T},UnitRange{Int64}} where {T} -irPairU{T}=Union{irPair{T},IndexRange} where {T} -IndexRange(ir::irPair{T}) where {T} =IndexRange(ir.first,ir.second) -IndexRange(ir::IndexRange)=IndexRange(ir.index,ir.range) - -start(ir::IndexRange)=range(ir).start -range(ir::IndexRange)=ir.range -range(i::Index)=1:dim(i) -ranges(irs::Tuple) = ntuple(i -> range(irs[i]), Val(length(irs))) -indices(irs::Tuple{Vararg{IndexRange}}) = map((ir)->ir.index ,irs) -indranges(ips::Tuple{Vararg{irPairU}}) = map((ip)->IndexRange(ip) ,ips) - -dim(ir::IndexRange)=dim(range(ir)) -dim(r::UnitRange{Int64})=r.stop-r.start+1 -dims(irs::Tuple{Vararg{IndexRange}})=map((ir)->dim(ir),irs) -redim(ip::irPair{T}) where {T} = redim(IndexRange(ip)) -redim(ir::IndexRange) = redim(ir.index,dim(ir),start(ir)-1) -redim(irs::Tuple{Vararg{IndexRange}}) = map((ir)->redim(ir) ,irs) - - -eachval(ir::IndexRange) = range(ir) -eachval(irs::Tuple{Vararg{IndexRange}}) = CartesianIndices(ranges(irs)) - -eachindval(irs::Tuple{Vararg{IndexRange}}) = (indices(irs).=> Tuple(ns) for ns in eachval(irs)) - - -#-------------------------------------------------------------------------------------- -# -# NDTensor level code which distinguishes between Dense and BlockSparse storage -# -function in_range(block_start::NTuple{N, Int64},block_end::NTuple{N, Int64},rs::UnitRange{Int64}...) where {N} - ret=true - for i in eachindex(rs) - if block_start[i]>rs[i].stop || block_end[i]=rs[i].start #in range? - istart=Base.max(dest_block_start[i],rs[i].start) - istop=Base.min(dest_block_end[i],rs[i].stop) - rs1[i]=istart-dest_block_start[i]+1:istop-dest_block_start[i]+1 - end - #@show rs1 - return Tuple(rs1) -end -# -# the range r will in general start at some index start(r) >1. This function -# counts how many QN blocks are between 1 and start(r) -# -function get_offset_block_count(in::QNIndex,r::UnitRange{Int64}) - # @show in r - rs=r.start-1 - if rs==0 - return 0,0 - end - @assert rs>0 - nb=0 - qns=space(in) - for n in eachindex(qns) - rs-=qns[n].second #dim of space - nb+=1 #increment block counts - if rs==0 - break - elseif rs<0 - nb-=1 - #@show "mid block slice" - #@show qns[n] rs nb - # @error("Slicing mid block is not supported yet.") - end - end - #@show nb rs - return nb,rs -end - -using StaticArrays -function get_offset_block_counts(inds::NTuple{N,IndsT},rs::UnitRange{Int64}...) where {N,IndsT} - # @show inds rs - dbs=StaticArrays.MVector{N,UInt}(undef) - shifts=StaticArrays.MVector{N,Int}(undef) - for i in eachindex(inds) - dbi,shift=get_offset_block_count(inds[i],rs[i]) - dbs[i]=dbi - shifts[i]=shift - end - return Block{N}(dbs),shifts -end - -function get_subtensor(T::BlockSparseTensor{ElT,N},new_inds,rs::UnitRange{Int64}...) where {ElT,N} - Ds = Vector{DenseTensor{ElT,N}}() - bs = Vector{Block{N}}() - dbs,=get_offset_block_counts(inds(T),rs...) - for (jj, b) in enumerate(eachnzblock(T)) - blockT = blockview(T, b) - if in_range(blockstart(T,b),blockend(T,b),rs...) - rs1=fix_ranges(blockstart(T,b),blockend(T,b),rs...) - #@show "In of range" b dims(blockT) rs rs1 NDTensors.blockstart(T,b) - push!(Ds,blockT[rs1...]) - bc=CartesianIndex(b)-CartesianIndex(dbs) - b=bc #Decrement block numbers by the number of skipped blocks. - push!(bs,b) - end - end - if length(Ds)==0 - return BlockSparseTensor(new_inds) - end - # - # JR: All attempts at building the new indices here at the NDTensors level failed. - # The only thing I could make work was to pass the new indices down from the ITensors - # level and use those. - # - #@show bs[1] inds(T) new_inds space(new_inds[1]) Block(bs[1][1]) - #bd=blockdim(new_inds[1],bs[1][1]) - #@show bs rs new_inds - T_sub = BlockSparseTensor(ElT, undef, bs, new_inds) - for ib in eachindex(Ds) - blockT_sub = bs[ib] - blockview(T_sub, blockT_sub) .= Ds[ib] - end - return T_sub -end - -function blockrange(T::Tensor{<:Number,N},tb::Block{N})::NTuple{N,UnitRange{Int64}} where {N} - bs=blockstart(T,tb) - be=blockend(T,tb) - return ntuple(i->bs[i]:be[i],N) -end - -function fix_ranges(dest_range::NTuple{N,UnitRange{Int64}},src_range::NTuple{N,UnitRange{Int64}},rs::UnitRange{Int64}...) where {N} - rs1=Vector{UnitRange{Int64}}(undef,N) - @assert length(rs)==N - #@show rs dest_range src_range - for i in eachindex(rs1) - @assert dest_range[i].start<=rs[i].stop || dest_range[i].stop>=rs[i].start #in range? - @assert src_range[i].start<=rs[i].stop || src_range[i].stop>=rs[i].start #in range? - ds_start=Base.max(dest_range[i].start,src_range[i].start) - dsrc = src_range[i].stop-src_range[i].start - istart=Base.max(ds_start,rs[i].start) -dest_range[i].start+1 - #istop=Base.min(ds_end,rs[i].stop) - #@show i ds_start dsrc istart - rs1[i]=istart:istart+dsrc - #@show rs1[i] - end - #@show rs1 - return Tuple(rs1) -end - -function set_subtensor(T::BlockSparseTensor{ElT,N},A::BlockSparseTensor{ElT,N},rs::UnitRange{Int64}...) where {ElT,N} -# @mpoc_assert nzblocks(T)==nzblocks(A) - # dbs,shifts=get_offset_block_counts(inds(T),rs...) - # dbs_ci=CartesianIndex(dbs) - insert=false - #rsa=ntuple(i->rs[i].start-dbs_ci[i]:rs[i].stop-dbs_ci[i],N) - rsa=ntuple(i->1:dim(inds(A)[i]),N) - #@show rs dbs_ci rsa - #@show inds(A) inds(T) - for ab in eachnzblock(A) - #@show CartesianIndex(dbs) CartesianIndex(ab) - blockA = blockview(A, ab) - #@show ab blockA blockstart(A,ab) blockend(A,ab) - - if in_range(blockstart(A,ab),blockend(A,ab),rsa...) - iA=blockstart(A,ab) - iT=ntuple(i->iA[i]+rs[i].start-1,N) - #iA=ntuple(i->iA[i]+dbs_ci[i],N) - #@show "inrange" - index_within_block,tb=blockindex(T,Tuple(iT)...) - blockT = blockview(T, tb) - #@show iA iT index_within_block tb blockT - if blockT==nothing - insertblock!(T,tb) - blockT = blockview(T, tb) - #@show "insert missing block" - #@show iA iT ab tb blockA blockT - insert=true - # index_within_block,tb=blockindex(T,Tuple(it)...) - # @show index_within_block tb - end - - # rs1=fix_ranges(blockrange(T,tb),blockrange(A,ab),rs...) - # @show rs1 - dA=[dims(blockA)...] - dT=[blockend(T,tb)...]-[index_within_block...]+fill(1,N) - rs1=ntuple(i->index_within_block[i]:index_within_block[i]+Base.min(dA[i],dT[i])-1,N) - #@show rs1 dA dT - if length(findall(dA.>dT))==0 - blockT[rs1...]=blockA #Dense assignment for each block - #@show blockT - else - rsa1=ntuple(i->rsa[i].start:rsa[i].start+Base.min(dA[i],dT[i])-1,N) - #@show rsa1 blockA[rsa1...] - blockT[rs1...]=blockA[rsa1...] #partial block Dense assignment for each block - @error "Incomplete bloc transfer." - @assert false - #@show blockT - end - end - # @show "-----block-------------" - end - if insert - #@show "blocks inserted" - #@assert false - end - #@show "---------------------------------" -end -function set_subtensor(T::DiagBlockSparseTensor{ElT,N},A::DiagBlockSparseTensor{ElT,N},rs::UnitRange{Int64}...) where {ElT,N} - dbs,=get_offset_block_counts(inds(T),rs...) - dbs_ci=CartesianIndex(dbs) - rsa=ntuple(i->rs[i].start-dbs_ci[i]:rs[i].stop-dbs_ci[i],N) - for ab in eachnzblock(A) - if in_range(blockstart(A,ab),blockend(A,ab),rsa...) - it=blockstart(A,ab) - it=ntuple(i->it[i]+dbs_ci[i],N) - index_within_block,tb=blockindex(T,Tuple(it)...) - blockT = blockview(T, tb) - blockA = blockview(A, ab) - rs1=fix_ranges(blockrange(T,tb),blockrange(A,ab),rs...) - #@show rs1 ab tb - - blockT[rs1...]=blockA #Diag assignment for each block - end - end -end - -function set_subtensor(T::DiagTensor{ElT},A::DiagTensor{ElT},rs::UnitRange{Int64}...) where {ElT} - if !all(y->y==rs[1],rs) - @error("set_subtensor(DiagTensor): All ranges must be the same, rs=$(rs).") - end - N=length(rs) - #only assign along the diagonal. - for i in rs[1] - is=ntuple(i1 -> i , N) - js=ntuple(j1 -> i-rs[1].start+1 , N) - T[is...]=A[js...] - end -end - -setindex!(T::BlockSparseTensor{ElT,N},A::BlockSparseTensor{ElT,N},irs::Vararg{UnitRange{Int64},N}) where {ElT,N} = set_subtensor(T,A,irs...) -setindex!(T::DiagTensor{ElT,N},A::DiagTensor{ElT,N},irs::Vararg{UnitRange{Int64},N}) where {ElT,N} = set_subtensor(T,A,irs...) -setindex!(T::DiagBlockSparseTensor{ElT,N},A::DiagBlockSparseTensor{ElT,N},irs::Vararg{UnitRange{Int64},N}) where {ElT,N} = set_subtensor(T,A,irs...) - - -#------------------------------------------------------------------------------------ -# -# ITensor level wrappers which allows us to handle the indices in a different manner -# depending on dense/block-sparse -# -function get_subtensor_wrapper(T::DenseTensor{ElT,N},new_inds,rs::UnitRange{Int64}...) where {ElT,N} - return ITensor(T[rs...],new_inds) -end - -function get_subtensor_wrapper(T::BlockSparseTensor{ElT,N},new_inds,rs::UnitRange{Int64}...) where {ElT,N} - return ITensor(get_subtensor(T,new_inds,rs...)) -end - - -function permute(indsT::T,irs::IndexRange...) where T<:(Tuple{Vararg{T, N}} where {N, T}) - ispec=indices(irs) #indices caller specified ranges for - inot=Tuple(noncommoninds(indsT,ispec)) #indices not specified by caller - isort=ispec...,inot... #all indices sorted so user specified ones are first. - isort_sub=redim(irs)...,inot... #all indices for subtensor - p=getperm(indsT, ntuple(n -> isort[n], length(isort))) - #@show p - return permute(isort_sub,p),permute((ranges(irs)...,ranges(inot)...),p) -end -# -# Use NDTensors T[3:4,1:3,1:6...] syntax to extract the subtensor. -# -function get_subtensor_ND(T::ITensor,irs::IndexRange...) - isub,rsub=permute(inds(T),irs...) #get indices and ranges for the subtensor - return get_subtensor_wrapper(tensor(T),isub,rsub...) #virtual dispatch based on Dense or BlockSparse -end - -function match_tagplev(i1::Index,i2s::Index...) - for i in eachindex(i2s) - #@show i tags(i1) tags(i2s[i]) plev(i1) plev(i2s[i]) - if tags(i1)==tags(i2s[i]) && plev(i1)==plev(i2s[i]) - return i - end - end - @error("match_tagplev: unable to find tag/plev for $i1 in Index set $i2s.") - return nothing -end - -function getperm_tagplev(s1, s2) - N = length(s1) - r = Vector{Int}(undef, N) - return map!(i -> match_tagplev(s1[i], s2...), r, 1:length(s1)) -end -# -# -# Permute is1 to be in the order of is2. -# BUT: match based on tags and plevs instread of IDs -# -function permute_tagplev(is1,is2) - length(is1) != length(is2) && throw( - ArgumentError( - "length of first index set, $(length(is1)) does not match length of second index set, $(length(is2))", - ), - ) - perm=getperm_tagplev(is1,is2) - #@show perm is1 is2 - return is1[invperm(perm)] -end - -# This operation is non trivial because we need to -# 1) Establish that the non-requested (not in irs) indices are identical (same ID) between T & A -# 2) Establish that requested (in irs) indices have the same tags and -# prime levels (dims are different so they cannot possibly have the same IDs) -# 3) Use a combination of tags/primes (for irs indices) or IDs (for non irs indices) to permut -# the indices of A prior to conversion to a Tensor. -# -function set_subtensor_ND(T::ITensor, A::ITensor,irs::IndexRange...) - ireqT=indices(irs) #indices caller requested ranges for - inotT=Tuple(noncommoninds(inds(T),ireqT)) #indices not requested by caller - ireqA=Tuple(noncommoninds(inds(A),inotT)) #Get the requested indices for a - inotA=Tuple(noncommoninds(inds(A),ireqA)) - if length(ireqT)!=length(ireqA) || inotA!=inotT - @show inotA inotT ireqT ireqA inotA!=inotT length(ireqT) length(ireqA) - @error("subtensor assign, incompatable indices\ndestination inds=$(inds(T)),\n source inds=$(inds(A)).") - @assert(false) - end - isortT=ireqT...,inotT... #all indices sorted so user specified ones are first. - p=getperm(inds(T), ntuple(n -> isortT[n], length(isortT))) # find p such that isort[p]==inds(T) - rsortT=permute((ranges(irs)...,ranges(inotT)...),p) #sorted ranges for T - ireqAp=permute_tagplev(ireqA,ireqT) #match based on tags & plev, NOT IDs since dims are different. - isortA=permute((ireqAp...,inotA...),p) #inotA is the same inotT, using inotA here for a less confusing read. - Ap=ITensors.permute(A,isortA...;allow_alias=true) - tensor(T)[rsortT...]=tensor(Ap) -end -function set_subtensor_ND(T::ITensor, v::Number,irs::IndexRange...) - ireqT=indices(irs) #indices caller requested ranges for - inotT=Tuple(noncommoninds(inds(T),ireqT)) #indices not requested by caller - isortT=ireqT...,inotT... #all indices sorted so user specified ones are first. - p=getperm(inds(T), ntuple(n -> isortT[n], length(isortT))) # find p such that isort[p]==inds(T) - rsortT=permute((ranges(irs)...,ranges(inotT)...),p) #sorted ranges for T - tensor(T)[rsortT...].=v -end - -getindex(T::ITensor, irs::Vararg{IndexRange,N}) where {N} = get_subtensor_ND(T,irs...) -getindex(T::ITensor, irs::Vararg{irPairU,N}) where {N} = get_subtensor_ND(T,indranges(irs)...) -setindex!(T::ITensor, A::ITensor,irs::Vararg{IndexRange,N}) where {N} = set_subtensor_ND(T,A,irs...) -setindex!(T::ITensor, A::ITensor,irs::Vararg{irPairU,N}) where {N} = set_subtensor_ND(T,A,indranges(irs)...) -setindex!(T::ITensor, v::Number,irs::Vararg{IndexRange,N}) where {N} = set_subtensor_ND(T,v,irs...) -setindex!(T::ITensor, v::Number,irs::Vararg{irPairU,N}) where {N} = set_subtensor_ND(T,v,indranges(irs)...) -#= -#-------------------------------------------------------------------------------- -# -# Some overloads of similar so we can easily create new ITensors from a template. -# -function similar(T::DenseTensor{ElT},is::Index...) where {ElT} - return ITensor(DenseTensor(ElT, undef, is)) -end -function similar(T::BlockSparseTensor{ElT},is::Index...) where {ElT} - return ITensor(BlockSparseTensor(ElT, undef, nzblocks(T), is)) -end - -function similar(T::DiagTensor{ElT},is::Index...) where {ElT} - ds=[dims(is)...] - if !all(y->y==ds[1],ds) - @error("similar(DiagTensor): All indices must have the same dimensions, dims(is)=$(dims(is)).") - end - N=dim(is[1]) - return ITensor(Diag(ElT, N),is) -end -function similar(T::DiagBlockSparseTensor{ElT},is::Index...) where {ElT} - ds=[dims(is)...] - if !all(y->y==ds[1],ds) - @error("similar(DiagTensor): All indices must have the same dimensions, dims(is)=$(dims(is)).") - end - return ITensor(DiagBlockSparseTensor(ElT, undef, nzblocks(T), is)) -end - - -function similar(T::ITensor,is::Index...) - return similar(tensor(T),is...) -end - =# - From 56de067ba36f7c415f55d53f4ba808cd87b85fc7 Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Sat, 1 Apr 2023 10:09:15 -0600 Subject: [PATCH 09/14] FInish cleanup of MPOMatrix-->InfiniteMPO conversion --- src/infinitempomatrix.jl | 75 ++++++++++++++++++++-------------------- 1 file changed, 37 insertions(+), 38 deletions(-) diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index 965c1c1..e22f940 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -7,6 +7,7 @@ mutable struct InfiniteMPOMatrix <: AbstractInfiniteMPS end translator(mpo::InfiniteMPOMatrix) = mpo.data.translator +data(mpo::InfiniteMPOMatrix) = mpo.data # TODO better printing? function Base.show(io::IO, M::InfiniteMPOMatrix) @@ -53,10 +54,12 @@ end # 0 0 ... l=2-->M-->l=1 0 0 # 0 0 ... 0 l=1-->M I # -# We need to capture the corner I's and the sub-diagonal Ms, and paste them together into on ITensor. -# In addition we need make all elements of Hm into order(4) tensors by adding dummy Dw=1 indices. +# We need to capture the corner I's and the sub-diagonal Ms, and paste them together into an ITensor. +# This is facilutated by making all elements of Hm into order(4) tensors by adding dummy Dw=1 indices. +# We ITensors.directsum() to join all the blocks into the final ITensor. +# This code should be 100% dense/blocks-sparse agnostic. # -function matrixITensorToITensor(Hm::Matrix{ITensor}) +function matrixITensorToITensor(Hm::Matrix{ITensor})::ITensor indexT=typeof(inds(Hm[1,1])[1]) T=eltype(Hm[1,1]) lx, ly = size(Hm) @@ -65,18 +68,15 @@ function matrixITensorToITensor(Hm::Matrix{ITensor}) # Extract the sub diagonal # Ms=map(n->Hm[lx-n+1,ly-n],1:lx-1) #Get the sub diagonal into an array. - #Ms=map(n->H[n+1,n],1:lx-1) #Get the sub diagonal into an array, reverse order. N=length(Ms) @assert N==lx-1 # # We need a dummy link index, but it has to have the right direction. - # Ms[1] should only have one link index, and it should be to the right. + # Ms[1] should only have one link index which should be pointing to the right. # @assert length(inds(Ms[1],tags="Link"))==1 ir,=inds(Ms[1],tags="Link") - left_dir=dir(dag(ir)) - tspace=ITensors.trivial_space(ir) - il0=Index(tspace;dir=left_dir,tags="Link,l=0") #Dw=1 left dummy index. + il0=Index(ITensors.trivial_space(ir);dir=dir(dag(ir)),tags="Link,l=0") #Dw=1 left dummy index. # # Convert edge Ms to order 4 ITensors using the dummy index. # @@ -87,13 +87,12 @@ function matrixITensorToITensor(Hm::Matrix{ITensor}) # can't use QN directions to decide left/right. Instead we look for common tags among the neighbouring Ms. # Right now the code assumes plev=0 for all Ms. # - ils,irs=indexT[],indexT[] #left the right link indices. - ir=il0 #set up recursion below. ir from the previous M, should have the same tags as il on the next M. + is=NTuple{2,indexT}[] + ir=il0 #set up recursion below. ir is from the previous M, should have the same tags as il on the next M. for n in 1:N il,=inds(Ms[n],tags=tags(ir)) #Find the left index - ir=noncommonind(Ms[n],il,tags="Link") #New right index by elimination - push!(ils,il) - push!(irs,ir) + ir=noncommonind(Ms[n],il,tags="Link") #Find the new right index by elimination + push!(is,(il,ir)) end # # Now direct sum everything non-empty in Hm, including those I ops in the corners. @@ -101,35 +100,35 @@ function matrixITensorToITensor(Hm::Matrix{ITensor}) # Details of the Link index tags like "l=n" do not matter from here on. # H=Hm[1,1]*onehot(T, il0' => 1)*onehot(T, dag(il0) => 1) #Bootstrap with the I op in the top left of Hm - H,ih=directsum(H=>il0',Ms[N]=>ils[N]) #1-D directsum to put M[N] directly below the I op + H,ih=directsum(H=>il0',Ms[N]=>is[N][1]) #1-D directsum to put M[N] directly below the I op ih=ih,dag(il0) #setup recusion. for i in N-1:-1:1 #2-D directsums to place new blocks below and to the right. - H,ih=directsum(H=>ih,Ms[i]=>(ils[i],irs[i]),tags=["Link","Link"]) + H,ih=directsum(H=>ih,Ms[i]=>is[i],tags=["Link,left","Link,right"]) end IN=Hm[N+1,N+1]*onehot(T, dag(il0) => 1)*onehot(T, ih[1] => dim(ih[1])) #I op in the bottom left of Hm - H,_=directsum(H=>ih[2],IN=>il0,tags="Link") #1-D directsum to the put I in the bottom row, to the right of M[1] - - ih=inds(H,tags="Link") - return H,ih[1],ih[2] + H,_=directsum(H=>ih[2],IN=>il0,tags="Link,right") #1-D directsum to the put I in the bottom row, to the right of M[1] + return H end - -function InfiniteMPO(H::InfiniteMPOMatrix) - temp = matrixITensorToITensor.(H) - new_H = CelledVector([x[1] for x in temp], translator(H)) - lis = CelledVector([x[2] for x in temp], translator(H)) - ris = CelledVector([x[3] for x in temp], translator(H)) - #retags the right_links - s = [commoninds(H[j][1, 1], H[j][end, end])[1] for j in 1:nsites(H)] - for j in 1:nsites(H) - newTag = "Link,c=$(getcell(s[j])),l=$(getsite(s[j]))" - temp = replacetags(ris[j], tags(ris[j]), newTag) - new_H[j] = replaceinds(new_H[j], [ris[j]], [temp]) - ris[j] = replacetags(ris[j], tags(ris[j]), newTag) - end - # joining the indexes - for j in 1:nsites(H) - temp = δ(dag(ris[j]), dag(lis[j + 1])) - new_H[j + 1] *= temp +# +# Hm is the InfiniteMPOMatrix +# Hs is an array of ITensors, one for the site in the unit cell. +# Hi is a CelledVector of ITensors. +# +function InfiniteMPO(Hm::InfiniteMPOMatrix) + Hs = matrixITensorToITensor.(data(Hm)) + Hi = CelledVector([H for H in Hs], translator(Hm)) + lis = CelledVector([inds(H,tags="left")[1] for H in Hs], translator(Hm)) + ris = CelledVector([inds(H,tags="right")[1] for H in Hs], translator(Hm)) + site_inds = [commoninds(Hm[j][1, 1], Hm[j][end, end])[1] for j in 1:nsites(Hm)] + # + # Create new tags with proper cell and link numbers. Also daisy chain + # all the indices so right index at j = dag(left index at j+1) + # + for j in 1:nsites(Hm) + newTag = "Link,c=$(getcell(site_inds[j])),l=$(getsite(site_inds[j]))" + ir = replacetags(ris[j], tags(ris[j]), newTag) #new right index + Hi[j] = replaceinds(Hi[j], [ris[j]], [ir]) + Hi[j + 1] = replaceinds(Hi[j + 1],[lis[j + 1]],[dag(ir)]) end - return InfiniteMPO(new_H) + return InfiniteMPO(Hi) end From 8fed6b5658b52f2a88ab818c63b013c5a95a1a30 Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Sat, 1 Apr 2023 18:30:09 -0600 Subject: [PATCH 10/14] Stop including subtensor.jl --- src/ITensorInfiniteMPS.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ITensorInfiniteMPS.jl b/src/ITensorInfiniteMPS.jl index 0d9eef0..d0ef847 100644 --- a/src/ITensorInfiniteMPS.jl +++ b/src/ITensorInfiniteMPS.jl @@ -52,7 +52,6 @@ include("vumps_generic.jl") include("vumps_localham.jl") include("vumps_nonlocalham.jl") include("vumps_mpo.jl") -include("subtensor.jl") export Cell, CelledVector, From 99ec250b8934cf191f2993108551dc4f8277ca77 Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Sun, 2 Apr 2023 07:47:54 -0600 Subject: [PATCH 11/14] Format --- src/infinitempomatrix.jl | 54 ++++----- src/models/models.jl | 2 +- test/test_iMPOConversions.jl | 221 +++++++++++++++++------------------ 3 files changed, 137 insertions(+), 140 deletions(-) diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index e22f940..cf21e2b 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -60,53 +60,53 @@ end # This code should be 100% dense/blocks-sparse agnostic. # function matrixITensorToITensor(Hm::Matrix{ITensor})::ITensor - indexT=typeof(inds(Hm[1,1])[1]) - T=eltype(Hm[1,1]) + indexT = typeof(inds(Hm[1, 1])[1]) + T = eltype(Hm[1, 1]) lx, ly = size(Hm) - @assert lx==ly + @assert lx == ly # # Extract the sub diagonal # - Ms=map(n->Hm[lx-n+1,ly-n],1:lx-1) #Get the sub diagonal into an array. - N=length(Ms) - @assert N==lx-1 + Ms = map(n -> Hm[lx - n + 1, ly - n], 1:(lx - 1)) #Get the sub diagonal into an array. + N = length(Ms) + @assert N == lx - 1 # # We need a dummy link index, but it has to have the right direction. # Ms[1] should only have one link index which should be pointing to the right. # - @assert length(inds(Ms[1],tags="Link"))==1 - ir,=inds(Ms[1],tags="Link") - il0=Index(ITensors.trivial_space(ir);dir=dir(dag(ir)),tags="Link,l=0") #Dw=1 left dummy index. + @assert length(inds(Ms[1]; tags="Link")) == 1 + ir, = inds(Ms[1]; tags="Link") + il0 = Index(ITensors.trivial_space(ir); dir=dir(dag(ir)), tags="Link,l=0") #Dw=1 left dummy index. # # Convert edge Ms to order 4 ITensors using the dummy index. # - Ms[1]*=onehot(T, il0 => 1) - Ms[N]*=onehot(T, dag(il0) => 1) #We don't need distinct index IDs for this. directsum makes all new index anyway. + Ms[1] *= onehot(T, il0 => 1) + Ms[N] *= onehot(T, dag(il0) => 1) #We don't need distinct index IDs for this. directsum makes all new index anyway. # # Create list of all left and right link indices on the Ms. In order to support dense storage we # can't use QN directions to decide left/right. Instead we look for common tags among the neighbouring Ms. # Right now the code assumes plev=0 for all Ms. # - is=NTuple{2,indexT}[] - ir=il0 #set up recursion below. ir is from the previous M, should have the same tags as il on the next M. + is = NTuple{2,indexT}[] + ir = il0 #set up recursion below. ir is from the previous M, should have the same tags as il on the next M. for n in 1:N - il,=inds(Ms[n],tags=tags(ir)) #Find the left index - ir=noncommonind(Ms[n],il,tags="Link") #Find the new right index by elimination - push!(is,(il,ir)) - end + il, = inds(Ms[n]; tags=tags(ir)) #Find the left index + ir = noncommonind(Ms[n], il; tags="Link") #Find the new right index by elimination + push!(is, (il, ir)) + end # # Now direct sum everything non-empty in Hm, including those I ops in the corners. # Order is critical. # Details of the Link index tags like "l=n" do not matter from here on. # - H=Hm[1,1]*onehot(T, il0' => 1)*onehot(T, dag(il0) => 1) #Bootstrap with the I op in the top left of Hm - H,ih=directsum(H=>il0',Ms[N]=>is[N][1]) #1-D directsum to put M[N] directly below the I op - ih=ih,dag(il0) #setup recusion. - for i in N-1:-1:1 #2-D directsums to place new blocks below and to the right. - H,ih=directsum(H=>ih,Ms[i]=>is[i],tags=["Link,left","Link,right"]) + H = Hm[1, 1] * onehot(T, il0' => 1) * onehot(T, dag(il0) => 1) #Bootstrap with the I op in the top left of Hm + H, ih = directsum(H => il0', Ms[N] => is[N][1]) #1-D directsum to put M[N] directly below the I op + ih = ih, dag(il0) #setup recusion. + for i in (N - 1):-1:1 #2-D directsums to place new blocks below and to the right. + H, ih = directsum(H => ih, Ms[i] => is[i]; tags=["Link,left", "Link,right"]) end - IN=Hm[N+1,N+1]*onehot(T, dag(il0) => 1)*onehot(T, ih[1] => dim(ih[1])) #I op in the bottom left of Hm - H,_=directsum(H=>ih[2],IN=>il0,tags="Link,right") #1-D directsum to the put I in the bottom row, to the right of M[1] + IN = Hm[N + 1, N + 1] * onehot(T, dag(il0) => 1) * onehot(T, ih[1] => dim(ih[1])) #I op in the bottom left of Hm + H, _ = directsum(H => ih[2], IN => il0; tags="Link,right") #1-D directsum to the put I in the bottom row, to the right of M[1] return H end # @@ -117,8 +117,8 @@ end function InfiniteMPO(Hm::InfiniteMPOMatrix) Hs = matrixITensorToITensor.(data(Hm)) Hi = CelledVector([H for H in Hs], translator(Hm)) - lis = CelledVector([inds(H,tags="left")[1] for H in Hs], translator(Hm)) - ris = CelledVector([inds(H,tags="right")[1] for H in Hs], translator(Hm)) + lis = CelledVector([inds(H; tags="left")[1] for H in Hs], translator(Hm)) + ris = CelledVector([inds(H; tags="right")[1] for H in Hs], translator(Hm)) site_inds = [commoninds(Hm[j][1, 1], Hm[j][end, end])[1] for j in 1:nsites(Hm)] # # Create new tags with proper cell and link numbers. Also daisy chain @@ -128,7 +128,7 @@ function InfiniteMPO(Hm::InfiniteMPOMatrix) newTag = "Link,c=$(getcell(site_inds[j])),l=$(getsite(site_inds[j]))" ir = replacetags(ris[j], tags(ris[j]), newTag) #new right index Hi[j] = replaceinds(Hi[j], [ris[j]], [ir]) - Hi[j + 1] = replaceinds(Hi[j + 1],[lis[j + 1]],[dag(ir)]) + Hi[j + 1] = replaceinds(Hi[j + 1], [lis[j + 1]], [dag(ir)]) end return InfiniteMPO(Hi) end diff --git a/src/models/models.jl b/src/models/models.jl index 271c6ff..6646d03 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -87,7 +87,7 @@ function InfiniteMPO(model::Model, s::CelledVector; kwargs...) end function InfiniteMPO(model::Model, s::CelledVector, translator::Function; kwargs...) - InfiniteMPO(InfiniteMPOMatrix(model,s,translator;kwargs...)) + return InfiniteMPO(InfiniteMPOMatrix(model, s, translator; kwargs...)) end function InfiniteMPOMatrix(model::Model, s::CelledVector; kwargs...) diff --git a/test/test_iMPOConversions.jl b/test/test_iMPOConversions.jl index 56ebaff..c57936c 100644 --- a/test/test_iMPOConversions.jl +++ b/test/test_iMPOConversions.jl @@ -6,141 +6,138 @@ using Test # with l,r terminating vectors, to make a finite lattice MPO. # function terminate(h::InfiniteMPO)::MPO - Ncell=nsites(h) - # left termination vector - il1=commonind(h[1],h[2]) - il0,=noncommoninds(h[1],il1,tags="Link") - l=ITensor(0.0,il0) - l[il0=>dim(il0)]=1.0 #assuming lower reg form in h - # right termination vector - iln=commonind(h[Ncell-1],h[Ncell]) - ilnp,=noncommoninds(h[Ncell],iln,tags="Link") - r=ITensor(0.0,ilnp) - r[ilnp=>1]=1.0 #assuming lower reg form in h - # build up a finite MPO - hf=MPO(Ncell) - hf[1]=dag(l)*h[1] #left terminate - hf[Ncell]=h[Ncell]*dag(r) #right terminate - for n in 2:Ncell-1 - hf[n]=h[n] #fill in the bulk. - end - return hf + Ncell = nsites(h) + # left termination vector + il1 = commonind(h[1], h[2]) + il0, = noncommoninds(h[1], il1; tags="Link") + l = ITensor(0.0, il0) + l[il0 => dim(il0)] = 1.0 #assuming lower reg form in h + # right termination vector + iln = commonind(h[Ncell - 1], h[Ncell]) + ilnp, = noncommoninds(h[Ncell], iln; tags="Link") + r = ITensor(0.0, ilnp) + r[ilnp => 1] = 1.0 #assuming lower reg form in h + # build up a finite MPO + hf = MPO(Ncell) + hf[1] = dag(l) * h[1] #left terminate + hf[Ncell] = h[Ncell] * dag(r) #right terminate + for n in 2:(Ncell - 1) + hf[n] = h[n] #fill in the bulk. + end + return hf end # # Terminate and then call expect # for inf ψ and finite h, which is already supported in src/infinitecanonicalmps.jl # function ITensors.expect(ψ::InfiniteCanonicalMPS, h::InfiniteMPO) - return expect(ψ,terminate(h)) #defer to src/infinitecanonicalmps.jl + return expect(ψ, terminate(h)) #defer to src/infinitecanonicalmps.jl end #H = ΣⱼΣn (½ S⁺ⱼS⁻ⱼ₊n + ½ S⁻ⱼS⁺ⱼ₊n + SᶻⱼSᶻⱼ₊n) -function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenbergNNN";NNN::Int64) - opsum = OpSum() - for n=1:NNN - J=1.0/n - opsum += J*0.5, "S+", 1, "S-", 1+n - opsum += J*0.5, "S-", 1, "S+", 1+n - opsum += J, "Sz", 1, "Sz", 1+n - end - return opsum +function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenbergNNN"; NNN::Int64) + opsum = OpSum() + for n in 1:NNN + J = 1.0 / n + opsum += J * 0.5, "S+", 1, "S-", 1 + n + opsum += J * 0.5, "S-", 1, "S+", 1 + n + opsum += J, "Sz", 1, "Sz", 1 + n + end + return opsum end -function ITensorInfiniteMPS.unit_cell_terms(::Model"hubbardNNN";NNN::Int64) - U::Float64=0.25 - t::Float64=1.0 - V::Float64=0.5 - opsum = OpSum() - opsum += (U, "Nupdn", 1) - for n=1:NNN - tj,Vj=t/n,V/n - opsum += -tj, "Cdagup", 1, "Cup", 1 + n - opsum += -tj, "Cdagup", 1 + n, "Cup", 1 - opsum += -tj, "Cdagdn", 1, "Cdn", 1 + n - opsum += -tj, "Cdagdn", 1 + n, "Cdn", 1 - opsum += Vj, "Ntot" , 1, "Ntot", 1 + n - end - return opsum +function ITensorInfiniteMPS.unit_cell_terms(::Model"hubbardNNN"; NNN::Int64) + U::Float64 = 0.25 + t::Float64 = 1.0 + V::Float64 = 0.5 + opsum = OpSum() + opsum += (U, "Nupdn", 1) + for n in 1:NNN + tj, Vj = t / n, V / n + opsum += -tj, "Cdagup", 1, "Cup", 1 + n + opsum += -tj, "Cdagup", 1 + n, "Cup", 1 + opsum += -tj, "Cdagdn", 1, "Cdn", 1 + n + opsum += -tj, "Cdagdn", 1 + n, "Cdn", 1 + opsum += Vj, "Ntot", 1, "Ntot", 1 + n + end + return opsum end function ITensors.space(::SiteType"FermionK", pos::Int; p=1, q=1, conserve_momentum=true) - if !conserve_momentum - return [QN("Nf", -p) => 1, QN("Nf", q - p) => 1] - else - return [ - QN(("Nf", -p), ("NfMom", -p * pos)) => 1, - QN(("Nf", q - p), ("NfMom", (q - p) * pos)) => 1, - ] - end + if !conserve_momentum + return [QN("Nf", -p) => 1, QN("Nf", q - p) => 1] + else + return [ + QN(("Nf", -p), ("NfMom", -p * pos)) => 1, + QN(("Nf", q - p), ("NfMom", (q - p) * pos)) => 1, + ] + end end function fermion_momentum_translator(i::Index, n::Integer; N=6) - ts = tags(i) - translated_ts = ITensorInfiniteMPS.translatecelltags(ts, n) - new_i = replacetags(i, ts => translated_ts) - for j in 1:length(new_i.space) - ch = new_i.space[j][1][1].val - mom = new_i.space[j][1][2].val - new_i.space[j] = Pair( - QN(("Nf", ch), ("NfMom", mom + n * N * ch)), new_i.space[j][2] - ) - end - return new_i - end + ts = tags(i) + translated_ts = ITensorInfiniteMPS.translatecelltags(ts, n) + new_i = replacetags(i, ts => translated_ts) + for j in 1:length(new_i.space) + ch = new_i.space[j][1][1].val + mom = new_i.space[j][1][2].val + new_i.space[j] = Pair(QN(("Nf", ch), ("NfMom", mom + n * N * ch)), new_i.space[j][2]) + end + return new_i +end function ITensors.op!(Op::ITensor, opname::OpName, ::SiteType"FermionK", s::Index...) - return ITensors.op!(Op, opname, SiteType("Fermion"), s...) + return ITensors.op!(Op, opname, SiteType("Fermion"), s...) end -@testset verbose=true "InfiniteMPOMatrix -> InfiniteMPO" begin - ferro(n) = "↑" - antiferro(n) = isodd(n) ? "↑" : "↓" +@testset verbose = true "InfiniteMPOMatrix -> InfiniteMPO" begin + ferro(n) = "↑" + antiferro(n) = isodd(n) ? "↑" : "↓" - models=[ - (Model"heisenbergNNN"(),"S=1/2"), - (Model"hubbardNNN"(),"Electron"), - - ] - @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af, qns=$qns" for (model,site) in models, qns=[false,true], Ncell in 2:6, NNN in 1:Ncell-1, Af in [true,false] -# @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af" for (model,site) in models, qns in [true], Ncell in 3:3, NNN in 2:2, Af in [false] - - if isodd(Ncell) && Af #skip test since Af state does fit inside odd cells. - continue - end - initstate(n) = Af ? antiferro(n) : ferro(n) - model_kwargs=(NNN=NNN,) - s = infsiteinds(site, Ncell; initstate,conserve_qns=qns); - ψ = InfMPS(s, initstate) - Hi=InfiniteMPO(model, s;model_kwargs...) - Hs=InfiniteSum{MPO}(model, s;model_kwargs...) - Es=expect(ψ,Hs) - Ei=expect(ψ,Hi) - #@show Es Ei - @test sum(Es[1:Ncell-NNN]) ≈ Ei atol=1e-14 - end + models = [(Model"heisenbergNNN"(), "S=1/2"), (Model"hubbardNNN"(), "Electron")] + @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af, qns=$qns" for (model, site) in + models, + qns in [false, true], + Ncell in 2:6, + NNN in 1:(Ncell - 1), + Af in [true, false] + # @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af" for (model,site) in models, qns in [true], Ncell in 3:3, NNN in 2:2, Af in [false] - @testset "FQHE Hamitonian" begin - N = 6 - model = Model"fqhe_2b_pot"() - model_params = (Vs= [1.0, 0.0, 1.0, 0.0, 0.1], Ly = 3.0, prec = 1e-8) - trf=fermion_momentum_translator - function initstate(n) - mod1(n, 3) == 1 && return 2 - return 1 - end - p = 1 - q = 3 - conserve_momentum = true - s = infsiteinds("FermionK", N; translator = trf, initstate, conserve_momentum, p, q); - ψ = InfMPS(s, initstate) - Hs=InfiniteSum{MPO}(model, s; model_params...); - Hi=InfiniteMPO(model, s,trf;model_params...) - #@show inds(Hi[1]) - Es=expect(ψ,Hs) - Ei=expect(ψ,Hi) - @show Es Ei - @test Es[1] ≈ Ei atol=1e-14 + if isodd(Ncell) && Af #skip test since Af state does fit inside odd cells. + continue end + initstate(n) = Af ? antiferro(n) : ferro(n) + model_kwargs = (NNN=NNN,) + s = infsiteinds(site, Ncell; initstate, conserve_qns=qns) + ψ = InfMPS(s, initstate) + Hi = InfiniteMPO(model, s; model_kwargs...) + Hs = InfiniteSum{MPO}(model, s; model_kwargs...) + Es = expect(ψ, Hs) + Ei = expect(ψ, Hi) + #@show Es Ei + @test sum(Es[1:(Ncell - NNN)]) ≈ Ei atol = 1e-14 + end + @testset "FQHE Hamitonian" begin + N = 6 + model = Model"fqhe_2b_pot"() + model_params = (Vs=[1.0, 0.0, 1.0, 0.0, 0.1], Ly=3.0, prec=1e-8) + trf = fermion_momentum_translator + function initstate(n) + mod1(n, 3) == 1 && return 2 + return 1 + end + p = 1 + q = 3 + conserve_momentum = true + s = infsiteinds("FermionK", N; translator=trf, initstate, conserve_momentum, p, q) + ψ = InfMPS(s, initstate) + Hs = InfiniteSum{MPO}(model, s; model_params...) + Hi = InfiniteMPO(model, s, trf; model_params...) + Es = expect(ψ, Hs) + Ei = expect(ψ, Hi) + #@show Es Ei + @test Es[1] ≈ Ei atol = 1e-14 + end end -nothing \ No newline at end of file +nothing From fa2f76902fea20e7586154dab063a773cb0ff44a Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Sun, 2 Apr 2023 10:28:58 -0600 Subject: [PATCH 12/14] Rename function matrixITensorToITensor--> cat_to_itensor --- src/infinitempomatrix.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index cf21e2b..3d48da0 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -59,7 +59,7 @@ end # We ITensors.directsum() to join all the blocks into the final ITensor. # This code should be 100% dense/blocks-sparse agnostic. # -function matrixITensorToITensor(Hm::Matrix{ITensor})::ITensor +function cat_to_itensor(Hm::Matrix{ITensor})::ITensor indexT = typeof(inds(Hm[1, 1])[1]) T = eltype(Hm[1, 1]) lx, ly = size(Hm) @@ -115,7 +115,7 @@ end # Hi is a CelledVector of ITensors. # function InfiniteMPO(Hm::InfiniteMPOMatrix) - Hs = matrixITensorToITensor.(data(Hm)) + Hs = cat_to_itensor.(data(Hm)) Hi = CelledVector([H for H in Hs], translator(Hm)) lis = CelledVector([inds(H; tags="left")[1] for H in Hs], translator(Hm)) ris = CelledVector([inds(H; tags="right")[1] for H in Hs], translator(Hm)) From 4d65e7c828a125ea23fee2944e1fcc44b96aa5d6 Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Sun, 2 Apr 2023 10:34:43 -0600 Subject: [PATCH 13/14] Find index hunting syntax inds(H;tags="Link") -->inds(H,tags="Link") --- src/infinitempomatrix.jl | 18 +++++++++--------- test/test_iMPOConversions.jl | 4 ++-- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index 3d48da0..4be036b 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -74,9 +74,9 @@ function cat_to_itensor(Hm::Matrix{ITensor})::ITensor # We need a dummy link index, but it has to have the right direction. # Ms[1] should only have one link index which should be pointing to the right. # - @assert length(inds(Ms[1]; tags="Link")) == 1 - ir, = inds(Ms[1]; tags="Link") - il0 = Index(ITensors.trivial_space(ir); dir=dir(dag(ir)), tags="Link,l=0") #Dw=1 left dummy index. + @assert length(inds(Ms[1], tags="Link")) == 1 + ir, = inds(Ms[1], tags="Link") + il0 = Index(ITensors.trivial_space(ir), dir=dir(dag(ir)), tags="Link,l=0") #Dw=1 left dummy index. # # Convert edge Ms to order 4 ITensors using the dummy index. # @@ -90,8 +90,8 @@ function cat_to_itensor(Hm::Matrix{ITensor})::ITensor is = NTuple{2,indexT}[] ir = il0 #set up recursion below. ir is from the previous M, should have the same tags as il on the next M. for n in 1:N - il, = inds(Ms[n]; tags=tags(ir)) #Find the left index - ir = noncommonind(Ms[n], il; tags="Link") #Find the new right index by elimination + il, = inds(Ms[n], tags=tags(ir)) #Find the left index + ir = noncommonind(Ms[n], il, tags="Link") #Find the new right index by elimination push!(is, (il, ir)) end # @@ -103,10 +103,10 @@ function cat_to_itensor(Hm::Matrix{ITensor})::ITensor H, ih = directsum(H => il0', Ms[N] => is[N][1]) #1-D directsum to put M[N] directly below the I op ih = ih, dag(il0) #setup recusion. for i in (N - 1):-1:1 #2-D directsums to place new blocks below and to the right. - H, ih = directsum(H => ih, Ms[i] => is[i]; tags=["Link,left", "Link,right"]) + H, ih = directsum(H => ih, Ms[i] => is[i], tags=["Link,left", "Link,right"]) end IN = Hm[N + 1, N + 1] * onehot(T, dag(il0) => 1) * onehot(T, ih[1] => dim(ih[1])) #I op in the bottom left of Hm - H, _ = directsum(H => ih[2], IN => il0; tags="Link,right") #1-D directsum to the put I in the bottom row, to the right of M[1] + H, _ = directsum(H => ih[2], IN => il0, tags="Link,right") #1-D directsum to the put I in the bottom row, to the right of M[1] return H end # @@ -117,8 +117,8 @@ end function InfiniteMPO(Hm::InfiniteMPOMatrix) Hs = cat_to_itensor.(data(Hm)) Hi = CelledVector([H for H in Hs], translator(Hm)) - lis = CelledVector([inds(H; tags="left")[1] for H in Hs], translator(Hm)) - ris = CelledVector([inds(H; tags="right")[1] for H in Hs], translator(Hm)) + lis = CelledVector([inds(H, tags="left")[1] for H in Hs], translator(Hm)) + ris = CelledVector([inds(H, tags="right")[1] for H in Hs], translator(Hm)) site_inds = [commoninds(Hm[j][1, 1], Hm[j][end, end])[1] for j in 1:nsites(Hm)] # # Create new tags with proper cell and link numbers. Also daisy chain diff --git a/test/test_iMPOConversions.jl b/test/test_iMPOConversions.jl index c57936c..cf88085 100644 --- a/test/test_iMPOConversions.jl +++ b/test/test_iMPOConversions.jl @@ -9,12 +9,12 @@ function terminate(h::InfiniteMPO)::MPO Ncell = nsites(h) # left termination vector il1 = commonind(h[1], h[2]) - il0, = noncommoninds(h[1], il1; tags="Link") + il0, = noncommoninds(h[1], il1, tags="Link") l = ITensor(0.0, il0) l[il0 => dim(il0)] = 1.0 #assuming lower reg form in h # right termination vector iln = commonind(h[Ncell - 1], h[Ncell]) - ilnp, = noncommoninds(h[Ncell], iln; tags="Link") + ilnp, = noncommoninds(h[Ncell], iln, tags="Link") r = ITensor(0.0, ilnp) r[ilnp => 1] = 1.0 #assuming lower reg form in h # build up a finite MPO From 426f9c37ad84edc040b4307c4d94c74cae92c793 Mon Sep 17 00:00:00 2001 From: Jan Reimers Date: Mon, 3 Apr 2023 09:22:47 -0600 Subject: [PATCH 14/14] Stop filtering inds using tags="link" Use network connectivity and set operations. --- src/infinitempomatrix.jl | 96 ++++++++++++++++++++++-------------- test/test_iMPOConversions.jl | 12 ++--- 2 files changed, 63 insertions(+), 45 deletions(-) diff --git a/src/infinitempomatrix.jl b/src/infinitempomatrix.jl index 4be036b..29ec966 100644 --- a/src/infinitempomatrix.jl +++ b/src/infinitempomatrix.jl @@ -59,9 +59,7 @@ end # We ITensors.directsum() to join all the blocks into the final ITensor. # This code should be 100% dense/blocks-sparse agnostic. # -function cat_to_itensor(Hm::Matrix{ITensor})::ITensor - indexT = typeof(inds(Hm[1, 1])[1]) - T = eltype(Hm[1, 1]) +function cat_to_itensor(Hm::Matrix{ITensor}, inds_array)::Tuple{ITensor,Index,Index} lx, ly = size(Hm) @assert lx == ly # @@ -71,54 +69,76 @@ function cat_to_itensor(Hm::Matrix{ITensor})::ITensor N = length(Ms) @assert N == lx - 1 # - # We need a dummy link index, but it has to have the right direction. - # Ms[1] should only have one link index which should be pointing to the right. - # - @assert length(inds(Ms[1], tags="Link")) == 1 - ir, = inds(Ms[1], tags="Link") - il0 = Index(ITensors.trivial_space(ir), dir=dir(dag(ir)), tags="Link,l=0") #Dw=1 left dummy index. - # # Convert edge Ms to order 4 ITensors using the dummy index. # + il0 = inds_array[1][1] + T = eltype(Ms[1]) Ms[1] *= onehot(T, il0 => 1) Ms[N] *= onehot(T, dag(il0) => 1) #We don't need distinct index IDs for this. directsum makes all new index anyway. # - # Create list of all left and right link indices on the Ms. In order to support dense storage we - # can't use QN directions to decide left/right. Instead we look for common tags among the neighbouring Ms. - # Right now the code assumes plev=0 for all Ms. + # Start direct sums to build the single ITensor. Order is critical to get the desired form. # - is = NTuple{2,indexT}[] - ir = il0 #set up recursion below. ir is from the previous M, should have the same tags as il on the next M. - for n in 1:N - il, = inds(Ms[n], tags=tags(ir)) #Find the left index - ir = noncommonind(Ms[n], il, tags="Link") #Find the new right index by elimination - push!(is, (il, ir)) + H = Hm[1, 1] * onehot(T, il0' => 1) * onehot(T, dag(il0) => 1) #Bootstrap with the I op in the top left of Hm + H, il = directsum(H => il0', Ms[N] => inds_array[N][1]) #1-D directsum to put M[N] directly below the I op + ir = dag(il0) #setup recusion. + for i in (N - 1):-1:1 #2-D directsums to place new blocks below and to the right. + H, (il, ir) = directsum( + H => (il, ir), Ms[i] => inds_array[i]; tags=["Link,left", "Link,right"] + ) end + IN = Hm[N + 1, N + 1] * onehot(T, dag(il0) => 1) * onehot(T, il => dim(il)) #I op in the bottom left of Hm + H, ir = directsum(H => ir, IN => il0; tags="Link,right") #1-D directsum to the put I in the bottom row, to the right of M[1] + + return H, il, ir +end + +function find_all_links(Hms::InfiniteMPOMatrix) + is = inds(Hms[1][1, 1]) #site inds + ir = noncommonind(Hms[1][2, 1], is) #This op should have one link index pointing to the right. # - # Now direct sum everything non-empty in Hm, including those I ops in the corners. - # Order is critical. - # Details of the Link index tags like "l=n" do not matter from here on. + # Set up return array of 2-tuples # - H = Hm[1, 1] * onehot(T, il0' => 1) * onehot(T, dag(il0) => 1) #Bootstrap with the I op in the top left of Hm - H, ih = directsum(H => il0', Ms[N] => is[N][1]) #1-D directsum to put M[N] directly below the I op - ih = ih, dag(il0) #setup recusion. - for i in (N - 1):-1:1 #2-D directsums to place new blocks below and to the right. - H, ih = directsum(H => ih, Ms[i] => is[i], tags=["Link,left", "Link,right"]) + indexT = typeof(ir) + TupleT = NTuple{2,indexT} + inds_array = Vector{TupleT}[] + # + # Make a dummy index + # + il0 = Index(ITensors.trivial_space(ir); dir=dir(dag(ir)), tags="Link,l=0") + # + # Loop over sites and nonzero matrix elements which are linked into the next + # and previous iMPOMatrix. + # + for n in 1:nsites(Hms) + Hm = Hms[n] + inds_n = TupleT[] + lx, ly = size(Hm) + @assert lx == ly + for iy in (ly - 1):-1:1 + ix = iy + 1 + il = ix < lx ? commonind(Hm[ix, iy], Hms[n - 1][ix + 1, iy + 1]) : dag(il0) + ir = iy > 1 ? commonind(Hm[ix, iy], Hms[n + 1][ix - 1, iy - 1]) : il0 + push!(inds_n, (il, ir)) + end + push!(inds_array, inds_n) end - IN = Hm[N + 1, N + 1] * onehot(T, dag(il0) => 1) * onehot(T, ih[1] => dim(ih[1])) #I op in the bottom left of Hm - H, _ = directsum(H => ih[2], IN => il0, tags="Link,right") #1-D directsum to the put I in the bottom row, to the right of M[1] - return H + return inds_array end + # # Hm is the InfiniteMPOMatrix -# Hs is an array of ITensors, one for the site in the unit cell. +# Hlrs is an array of {ITensor,Index,Index}s, one for each site in the unit cell. # Hi is a CelledVector of ITensors. # function InfiniteMPO(Hm::InfiniteMPOMatrix) - Hs = cat_to_itensor.(data(Hm)) - Hi = CelledVector([H for H in Hs], translator(Hm)) - lis = CelledVector([inds(H, tags="left")[1] for H in Hs], translator(Hm)) - ris = CelledVector([inds(H, tags="right")[1] for H in Hs], translator(Hm)) + inds_array = find_all_links(Hm) + Hlrs = cat_to_itensor.(Hm, inds_array) #return an array of {ITensor,Index,Index} + # + # Unpack the array of tuples into three arrays. And also get an array site indices. + # + Hi = CelledVector([Hlr[1] for Hlr in Hlrs], translator(Hm)) + ils = CelledVector([Hlr[2] for Hlr in Hlrs], translator(Hm)) + irs = CelledVector([Hlr[3] for Hlr in Hlrs], translator(Hm)) site_inds = [commoninds(Hm[j][1, 1], Hm[j][end, end])[1] for j in 1:nsites(Hm)] # # Create new tags with proper cell and link numbers. Also daisy chain @@ -126,9 +146,9 @@ function InfiniteMPO(Hm::InfiniteMPOMatrix) # for j in 1:nsites(Hm) newTag = "Link,c=$(getcell(site_inds[j])),l=$(getsite(site_inds[j]))" - ir = replacetags(ris[j], tags(ris[j]), newTag) #new right index - Hi[j] = replaceinds(Hi[j], [ris[j]], [ir]) - Hi[j + 1] = replaceinds(Hi[j + 1], [lis[j + 1]], [dag(ir)]) + ir = replacetags(irs[j], tags(irs[j]), newTag) #new right index + Hi[j] = replaceinds(Hi[j], [irs[j]], [ir]) + Hi[j + 1] = replaceinds(Hi[j + 1], [ils[j + 1]], [dag(ir)]) end return InfiniteMPO(Hi) end diff --git a/test/test_iMPOConversions.jl b/test/test_iMPOConversions.jl index cf88085..2c73e83 100644 --- a/test/test_iMPOConversions.jl +++ b/test/test_iMPOConversions.jl @@ -8,15 +8,13 @@ using Test function terminate(h::InfiniteMPO)::MPO Ncell = nsites(h) # left termination vector - il1 = commonind(h[1], h[2]) - il0, = noncommoninds(h[1], il1, tags="Link") + il0 = commonind(h[1], h[0]) l = ITensor(0.0, il0) l[il0 => dim(il0)] = 1.0 #assuming lower reg form in h # right termination vector - iln = commonind(h[Ncell - 1], h[Ncell]) - ilnp, = noncommoninds(h[Ncell], iln, tags="Link") - r = ITensor(0.0, ilnp) - r[ilnp => 1] = 1.0 #assuming lower reg form in h + iln = commonind(h[Ncell], h[Ncell + 1]) + r = ITensor(0.0, iln) + r[iln => 1] = 1.0 #assuming lower reg form in h # build up a finite MPO hf = MPO(Ncell) hf[1] = dag(l) * h[1] #left terminate @@ -75,6 +73,7 @@ function ITensors.space(::SiteType"FermionK", pos::Int; p=1, q=1, conserve_momen end function fermion_momentum_translator(i::Index, n::Integer; N=6) + #@show n ts = tags(i) translated_ts = ITensorInfiniteMPS.translatecelltags(ts, n) new_i = replacetags(i, ts => translated_ts) @@ -101,7 +100,6 @@ end Ncell in 2:6, NNN in 1:(Ncell - 1), Af in [true, false] - # @testset "H=$model, Ncell=$Ncell, NNN=$NNN, Antiferro=$Af" for (model,site) in models, qns in [true], Ncell in 3:3, NNN in 2:2, Af in [false] if isodd(Ncell) && Af #skip test since Af state does fit inside odd cells. continue