From 39a4b4b9cc37cd73e827e0d0f2d7f778137ad5d3 Mon Sep 17 00:00:00 2001 From: Ben Corbett <32752943+corbett5@users.noreply.github.com> Date: Wed, 21 Feb 2024 11:30:57 -0800 Subject: [PATCH] Moved code to new repo and setup packaging. (#2) --- Project.toml | 7 + docs/make.jl | 29 ++- examples/electrons.jl | 122 ++++++++++++ src/ITensorMPOConstruction.jl | 26 ++- src/MPOConstruction.jl | 274 +++++++++++++++++++++++++++ src/OpIDSum.jl | 336 ++++++++++++++++++++++++++++++++++ src/graph.jl | 176 ++++++++++++++++++ test/runtests.jl | 238 +++++++++++++++++++++++- 8 files changed, 1190 insertions(+), 18 deletions(-) create mode 100644 examples/electrons.jl create mode 100644 src/MPOConstruction.jl create mode 100644 src/OpIDSum.jl create mode 100644 src/graph.jl diff --git a/Project.toml b/Project.toml index 923e2ae..caf50c6 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,13 @@ uuid = "f1f99f3b-4b00-4d2f-a9c2-ce76efdc8f31" authors = ["Ben Corbett and contributors"] version = "0.1.0" +[deps] +ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Memoize = "c03570c3-d221-55d1-a50c-7939bbd78826" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" + [compat] julia = "1.10" diff --git a/docs/make.jl b/docs/make.jl index d6ea6e3..40ff161 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,23 +1,20 @@ using ITensorMPOConstruction using Documenter -DocMeta.setdocmeta!(ITensorMPOConstruction, :DocTestSetup, :(using ITensorMPOConstruction); recursive=true) +DocMeta.setdocmeta!( + ITensorMPOConstruction, :DocTestSetup, :(using ITensorMPOConstruction); recursive=true +) makedocs(; - modules=[ITensorMPOConstruction], - authors="Ben Corbett and contributors", - sitename="ITensorMPOConstruction.jl", - format=Documenter.HTML(; - canonical="https://ITensor.github.io/ITensorMPOConstruction.jl", - edit_link="main", - assets=String[], - ), - pages=[ - "Home" => "index.md", - ], + modules=[ITensorMPOConstruction], + authors="Ben Corbett and contributors", + sitename="ITensorMPOConstruction.jl", + format=Documenter.HTML(; + canonical="https://ITensor.github.io/ITensorMPOConstruction.jl", + edit_link="main", + assets=String[], + ), + pages=["Home" => "index.md"], ) -deploydocs(; - repo="github.com/ITensor/ITensorMPOConstruction.jl", - devbranch="main", -) +deploydocs(; repo="github.com/ITensor/ITensorMPOConstruction.jl", devbranch="main") diff --git a/examples/electrons.jl b/examples/electrons.jl new file mode 100644 index 0000000..a4f176d --- /dev/null +++ b/examples/electrons.jl @@ -0,0 +1,122 @@ +using ITensorMPOConstruction +using ITensors +using TimerOutputs + +function electron_op_cache_vec(sites::Vector{<:Index}) + operatorNames = [ + "I", + "Cdn", + "Cup", + "Cdagdn", + "Cdagup", + "Ndn", + "Nup", + "Cup * Cdn", + "Cup * Cdagdn", + "Cup * Ndn", + "Cdagup * Cdn", + "Cdagup * Cdagdn", + "Cdagup * Ndn", + "Nup * Cdn", + "Nup * Cdagdn", + "Nup * Ndn", + ] + + return [ + [OpInfo(ITensors.Op(name, n), sites[n]) for name in operatorNames] for + n in eachindex(sites) + ] +end + +function Fermi_Hubbard_momentum_space(N::Int, t::Real=1.0, U::Real=4.0)::MPO + ## Create the OpSum as per usual. + os = OpSum{Float64}() + for k in 1:N + epsilon = cospi(2 * k / N) + os .+= -t * epsilon, "Nup", k + os .+= -t * epsilon, "Ndn", k + end + + for p in 1:N + for q in 1:N + for k in 1:N + os .+= U / N, "Cdagup", mod1(p - k, N), "Cdagdn", mod1(q + k, N), "Cdn", q, "Cup", p + end + end + end + + sites = siteinds("Electron", N; conserve_qns=true) + + ## The only additional step required is to provide an operator basis in which to represent the OpSum. + opCacheVec = electron_op_cache_vec(sites) + return MPO_new(os, sites; basisOpCacheVec=opCacheVec) +end + +function electronic_structure_example(N::Int, complexBasisFunctions::Bool)::MPO + coeff() = rand() + 1im * complexBasisFunctions * rand() + + os = complexBasisFunctions ? OpSum{ComplexF64}() : OpSum{Float64}() + for a in 1:N + for b in a:N + epsilon_ab = coeff() + os .+= epsilon_ab, "Cdagup", a, "Cup", b + os .+= epsilon_ab, "Cdagdn", a, "Cdn", b + + a == b && continue + os .+= conj(epsilon_ab), "Cdagup", b, "Cup", a + os .+= conj(epsilon_ab), "Cdagdn", b, "Cdn", a + end + end + + for j in 1:N + for s_j in ("dn", "up") + for k in 1:N + s_k = s_j + (s_k, k) <= (s_j, j) && continue + + for l in 1:N + for s_l in ("dn", "up") + (s_l, l) <= (s_j, j) && continue + + for m in 1:N + s_m = s_l + (s_m, m) <= (s_k, k) && continue + + value = coeff() + os .+= value, "Cdag$s_j", j, "Cdag$s_l", l, "C$s_m", m, "C$s_k", k + os .+= conj(value), "Cdag$s_k", k, "Cdag$s_m", m, "C$s_l", l, "C$s_j", j + end + end + end + end + end + end + + sites = siteinds("Electron", N; conserve_qns=true) + + ## The only additional step required is to provide an operator basis in which to represent the OpSum. + opCacheVec = electron_op_cache_vec(sites) + return MPO_new(os, sites; basisOpCacheVec=opCacheVec) +end + +let + N = 50 + + # Ensure compilation + Fermi_Hubbard_momentum_space(10) + + println("Constructing the Fermi-Hubbard momentum space MPO for $N sites") + @time Fermi_Hubbard_momentum_space(N) +end + +let + N = 25 + + # Ensure compilation + electronic_structure_example(5, false) + + println("Constructing the Electronic Structure MPO for $N orbitals") + H = electronic_structure_example(N, false) +end + +nothing; diff --git a/src/ITensorMPOConstruction.jl b/src/ITensorMPOConstruction.jl index 3ced44e..aa91d0b 100644 --- a/src/ITensorMPOConstruction.jl +++ b/src/ITensorMPOConstruction.jl @@ -1,5 +1,29 @@ module ITensorMPOConstruction -# Write your package code here. +##################################### +# External packages +# +using ITensors +using LinearAlgebra +using SparseArrays +using Memoize +using TimerOutputs + +##################################### +# MPO Construction +# +include("OpIDSum.jl") +include("Graph.jl") +include("MPOConstruction.jl") + +##################################### +# Exports +# + +# OpIDSum.jl +export OpInfo, OpCacheVec, OpID, OpIDSum + +# MPOConstruction.jl +export MPO_new end diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl new file mode 100644 index 0000000..7495b99 --- /dev/null +++ b/src/MPOConstruction.jl @@ -0,0 +1,274 @@ +BlockSparseMatrix{C} = Dict{Tuple{Int,Int},Matrix{C}} + +@timeit function sparse_qr( + A::SparseMatrixCSC, tol::Real +)::Tuple{SparseMatrixCSC,SparseMatrixCSC,Vector{Int},Vector{Int},Int} + if tol < 0 + ret = qr(A) + else + ret = qr(A; tol=tol) + end + + return sparse(ret.Q), ret.R, ret.prow, ret.pcol, rank(ret) +end + +function for_non_zeros(f::Function, A::SparseMatrixCSC, maxRow::Int, maxCol::Int)::Nothing + @assert maxRow <= size(A, 1) + @assert maxCol <= size(A, 2) + + rows = rowvals(A) + vals = nonzeros(A) + for col in 1:maxCol + for idx in nzrange(A, col) + row = rows[idx] + if row <= maxRow + value = vals[idx] + f(value, row, col) + end + end + end +end + +function for_non_zeros_batch(f::Function, A::SparseMatrixCSC, maxCol::Int)::Nothing + @assert maxCol <= size(A, 2) "$maxCol, $(size(A, 2))" + + rows = rowvals(A) + vals = nonzeros(A) + for col in 1:maxCol + range = nzrange(A, col) + isempty(range) && continue + f((@view vals[range]), (@view rows[range]), col) + end +end + +function for_non_zeros_batch(f::Function, A::SparseMatrixCSC, cols::Vector{Int})::Nothing + rows = rowvals(A) + vals = nonzeros(A) + for col in cols + @assert 0 < col <= size(A, 2) + range = nzrange(A, col) + isempty(range) && continue + f((@view vals[range]), (@view rows[range]), col) + end +end + +@timeit function at_site!( + graphs::Dict{QN,MPOGraph{C}}, + n::Int, + sites::Vector{<:Index}, + tol::Real, + opCacheVec::OpCacheVec, +)::Tuple{Dict{QN,MPOGraph{C}},BlockSparseMatrix{C},Index} where {C} + hasQNs = hasqns(sites) + nextGraphs = Dict{QN,MPOGraph{C}}() + + matrix = BlockSparseMatrix{C}() + + qi = Vector{Pair{QN,Int}}() + outgoingLinkOffset = 0 + + mIdentityToID = Vector{Int}() + + for (inQN, g) in graphs + W = sparse_edge_weights(g) + Q, R, prow, pcol, rank = sparse_qr(W, tol) + + W = nothing + + #= + If we are at the last site, then R will be a 1x1 matrix containing an overall scaling factor + that we need to account for. + =# + if n == length(sites) + Q *= only(R) + end + + if hasQNs + append!(qi, [inQN => rank]) + end + + # Form the local transformation tensor. + # At the 0th dummy site we don't need to do this. + @timeit "constructing the sparse matrix" if n != 0 + ## TODO: This is wastefull + localMatrices = [ + my_matrix([lv.opOnSite], opCacheVec[n]; needsJWString=lv.needsJWString) for + lv in g.leftVertexValues + ] + + for_non_zeros(Q, left_size(g), rank) do weight, i, m + rightLink = m + outgoingLinkOffset + lv = left_value(g, prow[i]) + + matrixElement = get!(matrix, (lv.link, rightLink)) do + return zeros(C, dim(sites[n]), dim(sites[n])) + end + + matrixElement .+= weight * localMatrices[prow[i]] + end + end + + ## Free up some memory + Q = prow = nothing + + # Connect this output \tilde{A}_m to \tilde{B}_m = (R \vec{B})_m + # If we are processing the last site then there's no need to do this. + @timeit "constructing the next graphs" if n != length(sites) + nextGraphOfZeroFlux = get!(nextGraphs, inQN) do + return MPOGraph{C}() + end + + resize!(mIdentityToID, rank) + mIdentityToID .= 0 + for_non_zeros_batch(R, right_size(g)) do weights, ms, j + j = pcol[j] + + rv = right_value(g, j) + if !isempty(rv.ops) && rv.ops[end].n == n + 1 + onsite = pop!(rv.ops) + flux = opCacheVec[n + 1][onsite.id].qnFlux + newIsFermionic = xor(rv.isFermionic, opCacheVec[n + 1][onsite.id].isFermionic) + + rv = RightVertex(rv.ops, newIsFermionic) + + nextGraph = get!(nextGraphs, inQN + flux) do + return MPOGraph{C}() + end + + add_edges!(nextGraph, rv, rank, ms, outgoingLinkOffset, onsite, weights) + else + add_edges_id!( + nextGraphOfZeroFlux, + rv, + rank, + ms, + outgoingLinkOffset, + OpID(1, n + 1), + weights, + mIdentityToID, + ) + end + end + + if num_edges(nextGraphOfZeroFlux) == 0 + delete!(nextGraphs, inQN) + end + + for (_, nextGraph) in nextGraphs + empty!(nextGraph.leftVertexIDs) + end + end + + empty!(g) + outgoingLinkOffset += rank + end + + for (_, nextGraph) in nextGraphs + empty!(nextGraph.rightVertexIDs) + end + + if hasQNs + outgoingLink = Index(qi...; tags="Link,l=$n", dir=ITensors.Out) + else + outgoingLink = Index(outgoingLinkOffset; tags="Link,l=$n") + end + + return nextGraphs, matrix, outgoingLink +end + +@timeit function svdMPO_new( + ValType::Type{<:Number}, + os::OpIDSum{C}, + opCacheVec::OpCacheVec, + sites::Vector{<:Index}; + tol::Real=-1, +)::MPO where {C} + # TODO: This should be fixed. + @assert !ITensors.using_auto_fermion() + + N = length(sites) + + graphs = Dict{QN,MPOGraph{C}}(QN() => MPOGraph{C}(os)) + + H = MPO(sites) + + llinks = Vector{Index}(undef, N + 1) + for n in 0:N + graphs, symbolicMatrix, llinks[n + 1] = at_site!(graphs, n, sites, tol, opCacheVec) + + # For the 0th iteration we only care about constructing the graphs for the next site. + n == 0 && continue + + # Constructing the tensor from an array is much faster than setting the components of the ITensor directly. + # Especially for sparse tensors, and works for both sparse and dense tensors. + @timeit "creating ITensor" let + tensor = zeros( + ValType, + dim(dag(llinks[n])), + dim(llinks[n + 1]), + dim(prime(sites[n])), + dim(dag(sites[n])), + ) + + for ((leftLink, rightLink), localOpMatrix) in symbolicMatrix + tensor[leftLink, rightLink, :, :] = localOpMatrix + end + + H[n] = itensor( + tensor, + dag(llinks[n]), + llinks[n + 1], + prime(sites[n]), + dag(sites[n]); + tol=1e-10, + checkflux=false, + ) + end + end + + # Remove the dummy link indices from the left and right. + L = ITensor(llinks[1]) + L[end] = 1.0 + + R = ITensor(dag(llinks[N + 1])) + R[1] = 1.0 + + H[1] *= L + H[N] *= R + + return H +end + +function svdMPO_new(os::OpIDSum{C}, opCacheVec::OpCacheVec, sites; kwargs...)::MPO where {C} + # Function barrier to improve type stability + ValType = determine_val_type(os) + return svdMPO_new(ValType, os, opCacheVec, sites; kwargs...) +end + +function MPO_new( + os::OpIDSum, + sites::Vector{<:Index}, + opCacheVec::OpCacheVec; + basisOpCacheVec::Union{Nothing,OpCacheVec}=nothing, + kwargs..., +)::MPO + os, opCacheVec = prepare_opID_sum!(os, sites, opCacheVec, basisOpCacheVec) + + return svdMPO_new(os, opCacheVec, sites; kwargs...) +end + +function MPO_new(os::OpSum, sites::Vector{<:Index}; kwargs...)::MPO + opIDSum, opCacheVec = op_sum_to_opID_sum(os, sites) + return MPO_new(opIDSum, sites, opCacheVec; kwargs...) +end + +function sparsity(mpo::MPO)::Float64 + numEntries = 0 + numZeros = 0 + for tensor in mpo + numEntries += prod(size(tensor)) + numZeros += prod(size(tensor)) - ITensors.nnz(tensor) + end + + return numZeros / numEntries +end diff --git a/src/OpIDSum.jl b/src/OpIDSum.jl new file mode 100644 index 0000000..8a3202c --- /dev/null +++ b/src/OpIDSum.jl @@ -0,0 +1,336 @@ +struct OpInfo + matrix::Matrix + isFermionic::Bool + qnFlux::QN +end + +function OpInfo(localOp::Op, site::Index) + tensor = op(site, ITensors.which_op(localOp); ITensors.params(localOp)...) + isFermionic = has_fermion_string(ITensors.name(localOp), site) + qnFlux = flux(tensor) + + return OpInfo(copy(array(tensor)), isFermionic, isnothing(qnFlux) ? QN() : qnFlux) +end + +OpCacheVec = Vector{Vector{OpInfo}} + +struct OpID + id::Int16 + n::Int16 +end + +struct OpIDSum{C} + terms::Vector{OpID} + offsets::Vector{Int} + scalars::Vector{C} +end + +function OpIDSum{C}()::OpIDSum{C} where {C} + return OpIDSum(OpID[], Int[1], C[]) +end + +function Base.length(os::OpIDSum)::Int + return length(os.scalars) +end + +function Base.eachindex(os::OpIDSum)::UnitRange{Int} + return 1:length(os) +end + +function Base.getindex(os::OpIDSum, i::Integer) + return os.scalars[i], @view os.terms[os.offsets[i]:(os.offsets[i + 1] - 1)] +end + +function Base.push!(os::OpIDSum{C}, scalar::C, ops)::OpIDSum{C} where {C} + append!(os.terms, ops) + push!(os.offsets, os.offsets[end] + length(ops)) + push!(os.scalars, scalar) + + return os +end + +function Base.push!(os::OpIDSum{C}, scalar::C, ops::OpID...)::OpIDSum{C} where {C} + return push!(os, scalar, ops) +end + +function Base.push!(os::OpIDSum{C}, scalar::C, op::OpID)::OpIDSum{C} where {C} + push!(os.terms, op) + push!(os.offsets, os.offsets[end] + 1) + push!(os.scalars, scalar) + + return os +end + +function set_scalar!(os::OpIDSum{C}, i::Integer, scalar::C)::Nothing where {C} + os.scalars[i] = scalar + return nothing +end + +function add_to_scalar!(os::OpIDSum{C}, i::Integer, scalar::C)::Nothing where {C} + os.scalars[i] += scalar + return nothing +end + +# TODO: Define as `C`. Rename `coefficient_type`. +function determine_val_type(os::OpIDSum{C}) where {C} + for i in eachindex(os) + scalar, ops = os[i] + (!isreal(scalar)) && return ComplexF64 + end + + return Float64 +end + +@timeit function sort_each_term!( + os::OpIDSum, opCacheVec::OpCacheVec, sites::Vector{<:Index} +)::Nothing + isless_site(o1::OpID, o2::OpID) = o1.n < o2.n + + perm = Vector{Int}() + for j in eachindex(os) + scalar, ops = os[j] + Nt = length(ops) + + # Sort operators by site order, + # and keep the permutation used, perm, for analysis below + resize!(perm, Nt) + sortperm!(perm, ops; alg=InsertionSort, lt=isless_site) + ops .= ops[perm] + + # Identify fermionic operators, + # zeroing perm for bosonic operators, + parity = +1 + for (i, op) in enumerate(ops) + op.n > length(sites) && error( + "The OpSum contains an operator acting on site $(op.n) that extends beyond the number of sites $(length(sites)).", + ) + + if opCacheVec[op.n][op.id].isFermionic + parity = -parity + else + # Ignore bosonic operators in perm + # by zeroing corresponding entries + perm[i] = 0 + end + end + + if parity == -1 + error("Parity-odd fermionic terms not yet supported by AutoMPO") + end + + # Keep only fermionic op positions (non-zero entries) + filter!(!iszero, perm) + + # and account for anti-commuting, fermionic operators + # during above sort; put resulting sign into coef + set_scalar!(os, j, ITensors.parity_sign(perm) * scalar) + end +end + +@timeit function merge_terms!(os::OpIDSum)::Nothing + uniqueTermLocations = Dict{ + SubArray{OpID,1,Vector{OpID},Tuple{UnitRange{Int64}},true},Int + }() + for i in eachindex(os) + scalar, ops = os[i] + + loc = get!(uniqueTermLocations, ops, i) + loc == i && continue + + add_to_scalar!(os, loc, scalar) + set_scalar!(os, i, 0.0) + end + + return nothing +end + +function convert_to_basis(op::Matrix, basisCache::Vector{OpInfo})::Tuple{ComplexF64,Int} + function check_ratios_are_equal(a1::Number, b1::Number, ai::Number, bi::Number)::Bool + # TODO: Add a tolerance param here + return abs(a1 * bi - b1 * ai) <= 1e-10 + end + + for i in eachindex(basisCache) + basisOp = basisCache[i].matrix + + j = findfirst(val -> val != 0, basisOp) + + if all( + check_ratios_are_equal(basisOp[j], op[j], basisOp[k], op[k]) for + k in CartesianIndices(op) + ) + return op[j] / basisOp[j], i + end + end + + return 0, 0 +end + +function for_equal_sites(f::Function, ops::AbstractVector{OpID})::Nothing + i = 1 + while i <= length(ops) + j = i + while j <= length(ops) + (j == length(ops) || ops[i].n != ops[j + 1].n) && break + j += 1 + end + + f(i, j) + i = j + 1 + end + + return nothing +end + +@timeit function rewrite_in_operator_basis( + os::OpIDSum{C}, opCacheVec::OpCacheVec, basisOpCacheVec::OpCacheVec +)::OpIDSum{C} where {C} + @memoize Dict function convert_to_basis_memoized( + ops::AbstractVector{OpID} + )::Tuple{ComplexF64,Int} + n = ops[1].n + localOp = my_matrix(ops, opCacheVec[n]) + return convert_to_basis(localOp, basisOpCacheVec[n]) + end + + newOpIdSum = OpIDSum{C}() + + opsInBasis = Vector{OpID}() + for i in eachindex(os) + scalar, ops = os[i] + + empty!(opsInBasis) + + for_equal_sites(ops) do a, b + coeff, basisID = convert_to_basis_memoized(@view ops[a:b]) + + if basisID == 0 + error( + "The following operator product cannot be simplified into a single basis operator.\n" * + "\tOperator product = $(ops[a:b])\n", + ) + end + + scalar *= coeff + push!(opsInBasis, OpID(basisID, ops[a].n)) + end + + push!(newOpIdSum, C(scalar), opsInBasis) + end + + return newOpIdSum +end + +@timeit function op_sum_to_opID_sum( + os::OpSum{C}, sites::Vector{<:Index} +)::Tuple{OpIDSum{C},OpCacheVec} where {C} + N = length(sites) + + opsOnEachSite = [Dict{Op,Int}(Op("I", n) => 1) for n in 1:N] + + opID_sum = OpIDSum{C}() + + opIDTerm = Vector{OpID}() + for (i, term) in enumerate(os) + resize!(opIDTerm, 0) + for op in ITensors.terms(term) + op.which_op == "I" && continue + + n = ITensors.site(op) + opID = get!(opsOnEachSite[n], op, length(opsOnEachSite[n]) + 1) + push!(opIDTerm, OpID(opID, n)) + end + + push!(opID_sum, ITensors.coefficient(term), opIDTerm) + end + + opCacheVec = OpCacheVec() + for (n, opsOnSite) in enumerate(opsOnEachSite) + opCache = Vector{OpInfo}(undef, length(opsOnSite)) + + for (op, id) in opsOnSite + @assert ITensors.site(op) == n + opCache[id] = OpInfo(op, sites[n]) + end + + push!(opCacheVec, opCache) + end + + return opID_sum, opCacheVec +end + +@timeit function check_for_errors(os::OpIDSum, opCacheVec::OpCacheVec)::Nothing + for i in eachindex(os) + _, ops = os[i] + + flux = QN() + fermionParity = 0 + for j in eachindex(ops) + opj = ops[j] + flux += opCacheVec[opj.n][opj.id].qnFlux + fermionParity += opCacheVec[opj.n][opj.id].isFermionic + + if j < length(ops) + ops[j].n > ops[j + 1].n && error("The operators are not sorted by site.") + ops[j].n == ops[j + 1].n && + error("A site has more than one operator acting on it in a term.") + end + end + + flux != QN() && error("The term does not have zero flux.") + mod(fermionParity, 2) != 0 && error("Odd parity fermion terms not supported.") + end +end + +@timeit function prepare_opID_sum!( + os::OpIDSum, + sites::Vector{<:Index}, + opCacheVec::OpCacheVec, + basisOpCacheVec::Union{Nothing,OpCacheVec}, +)::Tuple{OpIDSum,OpCacheVec} + sort_each_term!(os, opCacheVec, sites) + merge_terms!(os) + + if !isnothing(basisOpCacheVec) + os, opCacheVec = rewrite_in_operator_basis(os, opCacheVec, basisOpCacheVec), + basisOpCacheVec + merge_terms!(os) + end + + check_for_errors(os, opCacheVec) + + # This is to account for the dummy site at the end + push!(opCacheVec, OpInfo[]) + + return os, opCacheVec +end + +function my_matrix( + term::AbstractVector{OpID}, opCache::Vector{OpInfo}; needsJWString::Bool=false +) + @assert all(op.n == term[1].n for op in term) + + if isempty(term) + localMatrix = Matrix{Float64}(I, size(opCache[begin].matrix)...) + else + localMatrix = prod(opCache[op.id].matrix for op in term) + end + + if needsJWString + # localMatrix can be returned directly from the opCache, and in this case we need to copy it + # before modification. + localMatrix = copy(localMatrix) + + # This is a weird way of applying the Jordan-Wigner string it but op("F", sites[n]) is very slow. + if size(localMatrix, 1) == 2 + localMatrix[:, 2] *= -1 + elseif size(localMatrix, 1) == 4 + localMatrix[:, 2] *= -1 + localMatrix[:, 3] *= -1 + else + error("Unknown fermionic site.") + end + end + + return localMatrix +end diff --git a/src/graph.jl b/src/graph.jl new file mode 100644 index 0000000..24c271a --- /dev/null +++ b/src/graph.jl @@ -0,0 +1,176 @@ +struct LeftVertex + link::Int32 + opOnSite::OpID + needsJWString::Bool +end + +function Base.hash(v::LeftVertex, h::UInt) + return hash((v.link, v.opOnSite), h) +end + +function Base.:(==)(v1::LeftVertex, v2::LeftVertex)::Bool + return (v1.link, v1.opOnSite) == (v2.link, v2.opOnSite) +end + +struct RightVertex + ops::Vector{OpID} + isFermionic::Bool +end + +function Base.hash(v::RightVertex, h::UInt) + return hash(v.ops, h) +end + +function Base.:(==)(v1::RightVertex, v2::RightVertex)::Bool + return v1.ops == v2.ops +end + +struct MPOGraph{C} + edgeLeftVertex::Vector{Int} + edgeRightVertex::Vector{Int} + edgeWeight::Vector{C} + leftVertexIDs::Dict{LeftVertex,Int} + rightVertexIDs::Dict{RightVertex,Int} + leftVertexValues::Vector{LeftVertex} + rightVertexValues::Vector{RightVertex} +end + +function MPOGraph{C}() where {C} + return MPOGraph( + Vector{Int}(), + Vector{Int}(), + Vector{C}(), + Dict{LeftVertex,Int}(), + Dict{RightVertex,Int}(), + Vector{LeftVertex}(), + Vector{RightVertex}(), + ) +end + +""" +Construct the graph of the sum split about a dummy site 0. + +Every term in the sum must be unique and sorted by site, have 0 flux, and even fermionic parity. +""" +@timeit function MPOGraph{C}(os::OpIDSum) where {C} + g = MPOGraph{C}() + lv = LeftVertex(0, OpID(0, 0), false) + + push!(g.leftVertexValues, lv) + + for i in eachindex(os) + scalar, ops = os[i] + + scalar == 0 && continue + + rv = RightVertex(reverse(ops), false) + push!(g.rightVertexValues, rv) + + push!(g.edgeLeftVertex, 1) + push!(g.edgeRightVertex, length(g.rightVertexValues)) + push!(g.edgeWeight, scalar) + end + + return g +end + +function Base.empty!(g::MPOGraph)::Nothing + empty!(g.edgeLeftVertex) + empty!(g.edgeRightVertex) + empty!(g.edgeWeight) + empty!(g.leftVertexValues) + empty!(g.rightVertexValues) + + return nothing +end + +function left_size(g::MPOGraph)::Int + return length(g.leftVertexValues) +end + +function right_size(g::MPOGraph)::Int + return length(g.rightVertexValues) +end + +function num_edges(g::MPOGraph)::Int + return length(g.edgeLeftVertex) +end + +function left_value(g::MPOGraph{C}, leftID::Int)::LeftVertex where {C} + return g.leftVertexValues[leftID] +end + +function right_value(g::MPOGraph{C}, rightID::Int)::RightVertex where {C} + return g.rightVertexValues[rightID] +end + +function add_edges!( + g::MPOGraph{C}, + rv::RightVertex, + rank::Int, + ms::AbstractVector{Int}, + mOffset::Int, + onsiteOp::OpID, + weights::AbstractVector{C}, +)::Nothing where {C} + rightID = get!(g.rightVertexIDs, rv) do + push!(g.rightVertexValues, rv) + return length(g.rightVertexValues) + end + + for i in 1:length(ms) + ms[i] > rank && return nothing + + lv = LeftVertex(ms[i] + mOffset, onsiteOp, rv.isFermionic) + + leftID = get!(g.leftVertexIDs, lv) do + push!(g.leftVertexValues, lv) + return length(g.leftVertexValues) + end + + push!(g.edgeLeftVertex, leftID) + push!(g.edgeRightVertex, rightID) + push!(g.edgeWeight, weights[i]) + end + + return nothing +end + +function add_edges_id!( + g::MPOGraph{C}, + rv::RightVertex, + rank::Int, + ms::AbstractVector{Int}, + mOffset::Int, + onsiteOp::OpID, + weights::AbstractVector{C}, + fooBar::Vector{Int}, +)::Nothing where {C} + rightID = get!(g.rightVertexIDs, rv) do + push!(g.rightVertexValues, rv) + return length(g.rightVertexValues) + end + + for i in 1:length(ms) + m = ms[i] + m > rank && return nothing + + if fooBar[m] == 0 + push!(g.leftVertexValues, LeftVertex(m + mOffset, onsiteOp, rv.isFermionic)) + fooBar[m] = length(g.leftVertexValues) + end + + push!(g.edgeLeftVertex, fooBar[m]) + push!(g.edgeRightVertex, rightID) + push!(g.edgeWeight, weights[i]) + end + + return nothing +end + +@timeit function sparse_edge_weights(g::MPOGraph{C})::SparseMatrixCSC{C,Int} where {C} + @assert length(g.edgeLeftVertex) == length(g.edgeRightVertex) + @assert length(g.edgeLeftVertex) == length(g.edgeWeight) + + return sparse(g.edgeLeftVertex, g.edgeRightVertex, g.edgeWeight) +end diff --git a/test/runtests.jl b/test/runtests.jl index 6d1375b..528abd4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,242 @@ using ITensorMPOConstruction +using ITensors using Test +function compare_MPOs(A::MPO, B::MPO; tol::Real=1e-7)::Nothing + diff = add(A, -1 * B; alg="directsum") + relativeDiff = sqrt(abs(dot(diff, diff))) / norm(B) + @test relativeDiff < tol + return nothing +end + +function test_from_OpSum( + os::OpSum, sites::Vector{<:Index}, basisOpCacheVec::Union{Nothing,OpCacheVec}, tol::Real +)::Tuple{MPO,MPO} + mpo = MPO_new(os, sites; tol=tol, basisOpCacheVec=basisOpCacheVec) + + if tol < 0 + mpoFromITensor = MPO(os, sites) + else + mpoFromITensor = MPO(os, sites; cutoff=tol) + end + + @test all(linkdims(mpo) .<= linkdims(mpoFromITensor)) + + compare_MPOs(mpo, mpoFromITensor) + + return mpo, mpoFromITensor +end + +function random_complex()::ComplexF64 + return 2 * rand(ComplexF64) - ComplexF64(1, 1) +end + +function test_IXYZ(N::Int64, tol::Real) + β = zeros(ComplexF64, N, 4) + for i in eachindex(β) + β[i] = random_complex() + end + + localOps = ["I", "X", "Y", "Z"] + + os = OpSum{ComplexF64}() + for string in CartesianIndex(ones(Int64, N)...):CartesianIndex((4 * ones(Int64, N))...) + string = Tuple(string) + + weight = 1 + args = [] + for i in eachindex(string) + whichOp = string[i] + weight *= β[i, whichOp] + append!(args, (localOps[whichOp], i)) + end + + if weight != 0 + os .+= weight, args... + end + end + + sites = siteinds("Qubit", N) + algMPO, iTensorMPO = test_from_OpSum(os, sites, nothing, tol) + + exact = MPO(sites) + + llinks = Vector{Index}(undef, N + 1) + llinks[1] = Index(1; tags="Link,l=0") + + for (n, site) in enumerate(sites) + localOp = β[n, 1] * op(site, "I") + localOp += β[n, 2] * op(site, "X") + localOp += β[n, 3] * op(site, "Y") + localOp += β[n, 4] * op(site, "Z") + + llinks[n + 1] = Index(1; tags="Link,l=$n") + exact[n] = ITensor( + ComplexF64, dag(llinks[n]), llinks[n + 1], prime(sites[n]), dag(sites[n]) + ) + exact[n][1, 1, :, :] = array(localOp, prime(sites[n]), dag(sites[n])) + end + + L = ITensor(llinks[1]) + L[end] = 1.0 + + R = ITensor(dag(llinks[N + 1])) + R[1] = 1.0 + + exact[1] *= L + exact[N] *= R + + compare_MPOs(algMPO, exact) + compare_MPOs(iTensorMPO, exact) + + return nothing +end + +function all_ops_of_weight_or_less(N::Integer, weight::Integer)::Vector + localOps = "X", "Y", "Z" + + if weight == 1 + ops = [] + for i in 1:N + for op_i in localOps + push!(ops, [op_i, i]) + end + end + else + ops = all_ops_of_weight_or_less(N, weight - 1) + for i in eachindex(ops) + string = ops[i] + if length(ops[i]) != 2 * (weight - 1) + continue + end + + lastSite = string[end] + for j in (lastSite + 1):N + for op in localOps + newOp = copy(string) + append!(newOp, [op, j]) + push!(ops, newOp) + end + end + end + end + + return ops +end + +function test_random_operator(N::Integer, maxWeight::Integer, tol::Real)::Nothing + ops = all_ops_of_weight_or_less(N, maxWeight) + + os = OpSum() + for op in ops + os .+= random_complex(), op... + end + + sites = siteinds("Qubit", N) + test_from_OpSum(os, sites, nothing, tol) + + return nothing +end + +function test_Fermi_Hubbard(N::Int, tol::Real)::Nothing + t, U = 1, 4 + sites = siteinds("Electron", N; conserve_qns=true) + + os = OpSum{Float64}() + + for k in 1:N + epsilon = cospi(2 * k / N) + os .+= -t * epsilon, "Nup", k + os .+= -t * epsilon, "Ndn", k + end + + for q in 1:N + for p in 1:N + for k in 1:N + os .+= U / N, "Cdagup", mod1(q - k, N), "Cdagdn", mod1(p + k, N), "Cdn", p, "Cup", q + end + end + end + + operatorNames = [ + "I", + "Cdn", + "Cup", + "Cdagdn", + "Cdagup", + "Ndn", + "Nup", + "Cup * Cdn", + "Cup * Cdagdn", + "Cup * Ndn", + "Cdagup * Cdn", + "Cdagup * Cdagdn", + "Cdagup * Ndn", + "Nup * Cdn", + "Nup * Cdagdn", + "Nup * Ndn", + ] + + opCacheVec = [ + [OpInfo(ITensors.Op(name, n), sites[n]) for name in operatorNames] for + n in eachindex(sites) + ] + + test_from_OpSum(os, sites, opCacheVec, tol) + return nothing +end + +ITensors.op(::OpName"00", ::SiteType"Qubit") = [1 0; 0 0] + +ITensors.op(::OpName"01", ::SiteType"Qubit") = [0 1; 0 0] + +ITensors.op(::OpName"10", ::SiteType"Qubit") = [0 0; 1 0] + +ITensors.op(::OpName"11", ::SiteType"Qubit") = [0 0; 0 1] + +function test_qft(N::Int64, applyR::Bool, tol::Real) + function bits_to_int(bits) + N = length(bits) + return sum(bitj << (N - j) for (j, bitj) in enumerate(bits)) + end + + os = OpSum{ComplexF64}() + for qBits in CartesianIndex(0 * ones(Int64, N)...):CartesianIndex((ones(Int64, N))...) + for qPrimeBits in + CartesianIndex(0 * ones(Int64, N)...):CartesianIndex((ones(Int64, N))...) + qBits = Tuple(qBits) + qPrimeBits = Tuple(qPrimeBits) + + args = Vector{Op}(undef, N) + for j in 1:N + args[j] = Op("$(qBits[j])$(qPrimeBits[j])", j) + end + + if applyR + q = bits_to_int(reverse(qBits)) + else + q = bits_to_int(qBits) + end + + qPrime = bits_to_int(qPrimeBits) + + weight = 1 / sqrt(2^N) * exp(1im * 2 * π * q * qPrime / 2^N) + add!(os, weight * Prod(args)) + end + end + + sites = siteinds("Qubit", N) + return algMPO, iTensorMPO = test_from_OpSum(os, sites, nothing, tol) +end + @testset "ITensorMPOConstruction.jl" begin - # Write your tests here. + test_IXYZ(5, -1) + + test_random_operator(8, 4, -1) + + test_qft(6, false, -1) + + test_qft(6, true, -1) + + test_Fermi_Hubbard(12, -1) end