From fb4c69da26e8a50f1e56ff1395afe25cad869826 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 10 Jan 2025 17:42:37 -0500 Subject: [PATCH 01/33] Get package to compile --- Project.toml | 8 +- .../abstractmps.jl | 2 +- src/ITensorMPS.jl | 2 +- src/abstractmps.jl | 62 ++- src/abstractprojmpo/diskprojmpo.jl | 4 +- src/dmrg.jl | 135 ++--- src/imports.jl | 39 +- src/mpo.jl | 3 +- src/mps.jl | 4 +- src/opsum_to_mpo/opsum_to_mpo.jl | 2 +- src/opsum_to_mpo/opsum_to_mpo_generic.jl | 4 +- src/opsum_to_mpo/opsum_to_mpo_qn.jl | 522 +++++++++--------- src/opsum_to_mpo/qnmatelem.jl | 60 +- src/solvers/expand.jl | 2 +- test/Project.toml | 1 - test/base/Project.toml | 3 +- 16 files changed, 429 insertions(+), 424 deletions(-) diff --git a/Project.toml b/Project.toml index 34654fb..5dc149a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorMPS" uuid = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" authors = ["Matthew Fishman ", "Miles Stoudenmire "] -version = "0.3.3" +version = "0.4.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" @@ -10,11 +10,11 @@ ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" IsApprox = "28f27b66-4bd8-47e7-9110-e2746eb8bed7" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SerializedElementArrays = "d3ce8812-9567-47e9-a7b5-65a6d70a3065" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" +TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" [weakdeps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" @@ -35,11 +35,10 @@ Adapt = "4.1.0" ChainRulesCore = "1.10" Compat = "4.16.0" HDF5 = "0.14, 0.15, 0.16, 0.17" -ITensors = "0.7" +ITensors = "0.8.0" IsApprox = "2.0.0" KrylovKit = "0.8.1" LinearAlgebra = "1.10" -NDTensors = "0.3.46" Observers = "0.2" PackageCompiler = "1, 2" Printf = "1.10" @@ -47,6 +46,7 @@ Random = "1.10" SerializedElementArrays = "0.1.0" Test = "1.10" TupleTools = "1.5.0" +TypeParameterAccessors = "0.2.1" ZygoteRules = "0.2.2" julia = "1.10" diff --git a/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl b/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl index 53c6980..a15d13f 100644 --- a/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl +++ b/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl @@ -3,7 +3,7 @@ using ChainRulesCore: ChainRulesCore, HasReverseMode, NoTangent, RuleConfig, rru using ITensors: ITensors, ITensor, dag, hassameinds, inds, itensor, mapprime, replaceprime, swapprime using ITensorMPS: ITensorMPS, MPO, MPS, apply, inner, siteinds -using NDTensors: datatype +## using NDTensors: datatype function ChainRulesCore.rrule( ::Type{T}, x::Vector{<:ITensor}; kwargs... diff --git a/src/ITensorMPS.jl b/src/ITensorMPS.jl index e90ec72..30add93 100644 --- a/src/ITensorMPS.jl +++ b/src/ITensorMPS.jl @@ -43,5 +43,5 @@ include("solvers/reducedlinearproblem.jl") include("solvers/linsolve.jl") include("solvers/expand.jl") include("lib/Experimental/src/Experimental.jl") -include("lib/ITensorMPSNamedDimsArraysExt/src/ITensorMPSNamedDimsArraysExt.jl") +## include("lib/ITensorMPSNamedDimsArraysExt/src/ITensorMPSNamedDimsArraysExt.jl") end diff --git a/src/abstractmps.jl b/src/abstractmps.jl index bd96f45..9e823c6 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -1,10 +1,14 @@ +# TODO: Define in `BackendSelection.jl`. +using ITensors: @Algorithm_str + using IsApprox: Approx, IsApprox -using ITensors: ITensors -using NDTensors: NDTensors, using_auto_fermion, scalartype, tensor +using ITensors: ITensors, ITensor, Index, commonind, commoninds, uniqueind, uniqueinds +## using NDTensors: NDTensors, using_auto_fermion, scalartype, tensor using ITensors.Ops: Prod -using ITensors.QuantumNumbers: QuantumNumbers, removeqn +## using ITensors.QuantumNumbers: QuantumNumbers, removeqn using ITensors.SiteTypes: SiteTypes, siteinds -using ITensors.TagSets: TagSets +## using ITensors.TagSets: TagSets +using LinearAlgebra: LinearAlgebra abstract type AbstractMPS end @@ -53,9 +57,9 @@ have type `ComplexF64`, return `ComplexF64`. """ promote_itensor_eltype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m) -NDTensors.scalartype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m) -NDTensors.scalartype(m::Array{ITensor}) = LinearAlgebra.promote_leaf_eltypes(m) -NDTensors.scalartype(m::Array{<:Array{ITensor}}) = LinearAlgebra.promote_leaf_eltypes(m) +## NDTensors.scalartype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m) +## NDTensors.scalartype(m::Array{ITensor}) = LinearAlgebra.promote_leaf_eltypes(m) +## NDTensors.scalartype(m::Array{<:Array{ITensor}}) = LinearAlgebra.promote_leaf_eltypes(m) """ eltype(m::MPS) @@ -711,9 +715,9 @@ for (fname, fname!) in [ (:(ITensors.noprime), :(ITensors.noprime!)), (:(ITensors.swapprime), :(ITensors.swapprime!)), (:(ITensors.replaceprime), :(ITensors.replaceprime!)), - (:(TagSets.addtags), :(ITensors.addtags!)), - (:(TagSets.removetags), :(ITensors.removetags!)), - (:(TagSets.replacetags), :(ITensors.replacetags!)), + (:(ITensors.addtags), :(ITensors.addtags!)), + (:(ITensors.removetags), :(ITensors.removetags!)), + (:(ITensors.replacetags), :(ITensors.replacetags!)), (:(ITensors.settags), :(ITensors.settags!)), ] @eval begin @@ -851,13 +855,13 @@ function hassamenuminds(::typeof(siteinds), M1::AbstractMPS, M2::AbstractMPS) end for (fname, fname!) in [ - (:(NDTensors.sim), :(sim!)), + (:(ITensors.sim), :(sim!)), (:(ITensors.prime), :(ITensors.prime!)), (:(ITensors.setprime), :(ITensors.setprime!)), (:(ITensors.noprime), :(ITensors.noprime!)), - (:(TagSets.addtags), :(ITensors.addtags!)), - (:(TagSets.removetags), :(ITensors.removetags!)), - (:(TagSets.replacetags), :(ITensors.replacetags!)), + (:(ITensors.addtags), :(ITensors.addtags!)), + (:(ITensors.removetags), :(ITensors.removetags!)), + (:(ITensors.replacetags), :(ITensors.replacetags!)), (:(ITensors.settags), :(ITensors.settags!)), ] @eval begin @@ -1602,11 +1606,11 @@ out-of-place with `orthogonalize`. """ function orthogonalize!(M::AbstractMPS, j::Int; maxdim=nothing, normalize=nothing) # TODO: Delete `maxdim` and `normalize` keyword arguments. - @debug_check begin - if !(1 <= j <= length(M)) - error("Input j=$j to `orthogonalize!` out of range (valid range = 1:$(length(M)))") - end - end + ## @debug_check begin + ## if !(1 <= j <= length(M)) + ## error("Input j=$j to `orthogonalize!` out of range (valid range = 1:$(length(M)))") + ## end + ## end while leftlim(M) < (j - 1) (leftlim(M) < 0) && setleftlim!(M, 0) b = leftlim(M) + 1 @@ -1772,12 +1776,12 @@ end _isodd_fermionic_parity(s::Index, ::Integer) = false -function _isodd_fermionic_parity(s::QNIndex, n::Integer) - qn_n = qn(space(s)[n]) - fermionic_qn_pos = findfirst(q -> isfermionic(q), qn_n) - isnothing(fermionic_qn_pos) && return false - return isodd(val(qn_n[fermionic_qn_pos])) -end +## function _isodd_fermionic_parity(s::QNIndex, n::Integer) +## qn_n = qn(space(s)[n]) +## fermionic_qn_pos = findfirst(q -> isfermionic(q), qn_n) +## isnothing(fermionic_qn_pos) && return false +## return isodd(val(qn_n[fermionic_qn_pos])) +## end function _fermionic_swap(s1::Index, s2::Index) T = ITensor(QN(), s1', s2', dag(s1), dag(s2)) @@ -1911,7 +1915,7 @@ replacesites!(ψ::AbstractMPS, args...; kwargs...) = setindex!(ψ, args...; kwar replacesites(ψ::AbstractMPS, args...; kwargs...) = setindex!(copy(ψ), args...; kwargs...) _number_inds(s::Index) = 1 -_number_inds(s::IndexSet) = length(s) +_number_inds(s::Union{Vector{<:Index},Tuple{Vararg{Index}}}) = length(s) _number_inds(sites) = sum(_number_inds(s) for s in sites) """ @@ -2356,9 +2360,9 @@ function splitblocks(::typeof(linkinds), M::AbstractMPS; tol=0) end removeqns(M::AbstractMPS) = map(removeqns, M; set_limits=false) -function QuantumNumbers.removeqn(M::AbstractMPS, qn_name::String) - return map(m -> removeqn(m, qn_name), M; set_limits=false) -end +## function QuantumNumbers.removeqn(M::AbstractMPS, qn_name::String) +## return map(m -> removeqn(m, qn_name), M; set_limits=false) +## end # # Broadcasting diff --git a/src/abstractprojmpo/diskprojmpo.jl b/src/abstractprojmpo/diskprojmpo.jl index dd9d1df..7630a15 100644 --- a/src/abstractprojmpo/diskprojmpo.jl +++ b/src/abstractprojmpo/diskprojmpo.jl @@ -1,3 +1,5 @@ +using ITensors: OneITensor +using SerializedElementArrays: SerializedElementVector """ A DiskProjMPO computes and stores the projection of an @@ -26,7 +28,7 @@ mutable struct DiskProjMPO <: AbstractProjMPO rpos::Int nsite::Int H::MPO - LR::DiskVector{ITensor} + LR::SerializedElementVector{ITensor} Lcache::Union{ITensor,OneITensor} lposcache::Union{Int,Nothing} Rcache::Union{ITensor,OneITensor} diff --git a/src/dmrg.jl b/src/dmrg.jl index 39ca7a2..3e0e65c 100644 --- a/src/dmrg.jl +++ b/src/dmrg.jl @@ -1,6 +1,6 @@ using Adapt: adapt using KrylovKit: eigsolve -using NDTensors: scalartype, timer +## using NDTensors: scalartype, timer using Printf: @printf using TupleTools: TupleTools @@ -55,7 +55,8 @@ function dmrg(H::MPO, Ms::Vector{MPS}, psi0::MPS, sweeps::Sweeps; weight=true, k return dmrg(PMM, psi0, sweeps; kwargs...) end -using NDTensors.TypeParameterAccessors: unwrap_array_type +using TypeParameterAccessors: unwrap_array_type + """ dmrg(H::MPO, psi0::MPS; kwargs...) dmrg(H::MPO, psi0::MPS, sweeps::Sweeps; kwargs...) @@ -179,12 +180,12 @@ function dmrg( ) end - @debug_check begin - # Debug level checks - # Enable with ITensors.enable_debug_checks() - checkflux(psi0) - checkflux(PH) - end + ## @debug_check begin + ## # Debug level checks + ## # Enable with ITensors.enable_debug_checks() + ## checkflux(psi0) + ## checkflux(PH) + ## end psi = copy(psi0) N = length(psi) @@ -217,36 +218,36 @@ function dmrg( end for (b, ha) in sweepnext(N) - @debug_check begin - checkflux(psi) - checkflux(PH) - end - - @timeit_debug timer "dmrg: position!" begin - PH = position!(PH, psi, b) - end - - @debug_check begin - checkflux(psi) - checkflux(PH) - end - - @timeit_debug timer "dmrg: psi[b]*psi[b+1]" begin - phi = psi[b] * psi[b + 1] - end - - @timeit_debug timer "dmrg: eigsolve" begin - vals, vecs = eigsolve( - PH, - phi, - 1, - eigsolve_which_eigenvalue; - ishermitian, - tol=eigsolve_tol, - krylovdim=eigsolve_krylovdim, - maxiter=eigsolve_maxiter, - ) - end + ## @debug_check begin + ## checkflux(psi) + ## checkflux(PH) + ## end + + ## @timeit_debug timer "dmrg: position!" begin + PH = position!(PH, psi, b) + ## end + + ## @debug_check begin + ## checkflux(psi) + ## checkflux(PH) + ## end + + ## @timeit_debug timer "dmrg: psi[b]*psi[b+1]" begin + phi = psi[b] * psi[b + 1] + ## end + + ## @timeit_debug timer "dmrg: eigsolve" begin + vals, vecs = eigsolve( + PH, + phi, + 1, + eigsolve_which_eigenvalue; + ishermitian, + tol=eigsolve_tol, + krylovdim=eigsolve_krylovdim, + maxiter=eigsolve_maxiter, + ) + ## end energy = vals[1] ## Right now there is a conversion problem in CUDA.jl where `UnifiedMemory` Arrays are being converted @@ -263,40 +264,40 @@ function dmrg( drho = nothing if noise(sweeps, sw) > 0 - @timeit_debug timer "dmrg: noiseterm" begin - # Use noise term when determining new MPS basis. - # This is used to preserve the element type of the MPS. - elt = real(scalartype(psi)) - drho = elt(noise(sweeps, sw)) * noiseterm(PH, phi, ortho) - end + ## @timeit_debug timer "dmrg: noiseterm" begin + # Use noise term when determining new MPS basis. + # This is used to preserve the element type of the MPS. + elt = real(scalartype(psi)) + drho = elt(noise(sweeps, sw)) * noiseterm(PH, phi, ortho) + ## end end - @debug_check begin - checkflux(phi) - end + ## @debug_check begin + ## checkflux(phi) + ## end - @timeit_debug timer "dmrg: replacebond!" begin - spec = replacebond!( - PH, - psi, - b, - phi; - maxdim=maxdim(sweeps, sw), - mindim=mindim(sweeps, sw), - cutoff=cutoff(sweeps, sw), - eigen_perturbation=drho, - ortho, - normalize=true, - which_decomp, - svd_alg, - ) - end + ## @timeit_debug timer "dmrg: replacebond!" begin + spec = replacebond!( + PH, + psi, + b, + phi; + maxdim=maxdim(sweeps, sw), + mindim=mindim(sweeps, sw), + cutoff=cutoff(sweeps, sw), + eigen_perturbation=drho, + ortho, + normalize=true, + which_decomp, + svd_alg, + ) + ## end maxtruncerr = max(maxtruncerr, spec.truncerr) - @debug_check begin - checkflux(psi) - checkflux(PH) - end + ## @debug_check begin + ## checkflux(psi) + ## checkflux(PH) + ## end if outputlevel >= 2 @printf("Sweep %d, half %d, bond (%d,%d) energy=%s\n", sw, ha, b, b + 1, energy) diff --git a/src/imports.jl b/src/imports.jl index d3a4b93..e72aba8 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -98,28 +98,27 @@ import Base.Broadcast: broadcastable, instantiate -import ..ITensors.NDTensors: - Algorithm, - @Algorithm_str, - EmptyNumber, - _Tuple, - _NTuple, - blas_get_num_threads, - datatype, - dense, - diagind, - disable_auto_fermion, - double_precision, - eachblock, - eachdiagblock, - enable_auto_fermion, - fill!!, - randn!!, - permutedims, - permutedims! +## import ..ITensors.NDTensors: +## Algorithm, +## @Algorithm_str, +## EmptyNumber, +## _Tuple, +## _NTuple, +## blas_get_num_threads, +## datatype, +## dense, +## diagind, +## disable_auto_fermion, +## double_precision, +## eachblock, +## eachdiagblock, +## enable_auto_fermion, +## fill!!, +## randn!!, +## permutedims, +## permutedims! import ..ITensors: - AbstractRNG, Apply, apply, argument, diff --git a/src/mpo.jl b/src/mpo.jl index 9acd561..4e19724 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -1,6 +1,7 @@ using Adapt: adapt using LinearAlgebra: dot -using Random: Random +using Random: Random, AbstractRNG +using ITensors: @ts_str, Algorithm, outer using ITensors.Ops: OpSum using ITensors.SiteTypes: SiteTypes, siteind, siteinds diff --git a/src/mps.jl b/src/mps.jl index de71318..a1029a6 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -1,6 +1,6 @@ using Adapt: adapt -using NDTensors: using_auto_fermion -using Random: Random +## using NDTensors: using_auto_fermion +using Random: Random, AbstractRNG using ITensors.SiteTypes: SiteTypes, siteind, siteinds, state """ diff --git a/src/opsum_to_mpo/opsum_to_mpo.jl b/src/opsum_to_mpo/opsum_to_mpo.jl index 114eec3..9142e02 100644 --- a/src/opsum_to_mpo/opsum_to_mpo.jl +++ b/src/opsum_to_mpo/opsum_to_mpo.jl @@ -1,4 +1,4 @@ -using NDTensors: using_auto_fermion +## using NDTensors: using_auto_fermion # `ValType::Type{<:Number}` is used instead of `ValType::Type` for efficiency, possibly due to increased method specialization. # See https://github.com/ITensor/ITensors.jl/pull/1183. diff --git a/src/opsum_to_mpo/opsum_to_mpo_generic.jl b/src/opsum_to_mpo/opsum_to_mpo_generic.jl index 915fd0e..009a699 100644 --- a/src/opsum_to_mpo/opsum_to_mpo_generic.jl +++ b/src/opsum_to_mpo/opsum_to_mpo_generic.jl @@ -1,4 +1,4 @@ -using NDTensors: using_auto_fermion +## using NDTensors: using_auto_fermion using ITensors.Ops: Ops, Op, OpSum, Scaled, Sum, coefficient using ITensors.SiteTypes: has_fermion_string, op @@ -110,7 +110,7 @@ function Base.copyto!(os, bc::Broadcast.Broadcasted{OpSumAddTermStyle,<:Any,type end # XXX: Create a new function name for this. -isempty(op_qn::Pair{Vector{Op},QN}) = isempty(op_qn.first) +## isempty(op_qn::Pair{Vector{Op},QN}) = isempty(op_qn.first) # the key type is Prod{Op} for the dense case # and is Pair{Prod{Op},QN} for the QN conserving case diff --git a/src/opsum_to_mpo/opsum_to_mpo_qn.jl b/src/opsum_to_mpo/opsum_to_mpo_qn.jl index 25f75b3..a1f31a2 100644 --- a/src/opsum_to_mpo/opsum_to_mpo_qn.jl +++ b/src/opsum_to_mpo/opsum_to_mpo_qn.jl @@ -1,261 +1,261 @@ -using NDTensors: using_auto_fermion - -# `ValType::Type{<:Number}` is used instead of `ValType::Type` for efficiency, possibly due to increased method specialization. -# See https://github.com/ITensor/ITensors.jl/pull/1183. -function qn_svdMPO( - ValType::Type{<:Number}, os::OpSum{C}, sites; mindim=1, maxdim=typemax(Int), cutoff=1e-15 -)::MPO where {C} - N = length(sites) - - # Specifying the element type with `Dict{QN,Matrix{ValType}}[...]` improves type inference and therefore efficiency. - # See https://github.com/ITensor/ITensors.jl/pull/1183. - Vs = Dict{QN,Matrix{ValType}}[Dict{QN,Matrix{ValType}}() for n in 1:(N + 1)] - sparse_MPO = [QNMatElem{Scaled{C,Prod{Op}}}[] for n in 1:N] - - function crosses_bond(t::Scaled{C,Prod{Op}}, n::Int) - return (only(site(t[1])) <= n <= only(site(t[end]))) - end - - # A cache of the ITensor operators on a certain site - # of a certain type - op_cache = Dict{Pair{String,Int},ITensor}() - function calcQN(term::Vector{Op}) - q = QN() - for st in term - op_tensor = get(op_cache, which_op(st) => only(site(st)), nothing) - if op_tensor === nothing - op_tensor = op(sites[only(site(st))], which_op(st); params(st)...) - op_cache[which_op(st) => only(site(st))] = op_tensor - end - q -= flux(op_tensor) - end - return q - end - - Hflux = -calcQN(terms(first(terms(os)))) - - rightmap = Dict{Pair{Vector{Op},QN},Int}() - next_rightmap = Dict{Pair{Vector{Op},QN},Int}() - - for n in 1:N - h_sparse = Dict{QN,Vector{MatElem{ValType}}}() - - leftmap = Dict{Pair{Vector{Op},QN},Int}() - for term in os - crosses_bond(term, n) || continue - - left = filter(t -> (only(site(t)) < n), terms(term)) - onsite = filter(t -> (only(site(t)) == n), terms(term)) - right = filter(t -> (only(site(t)) > n), terms(term)) - - lqn = calcQN(left) - sqn = calcQN(onsite) - - bond_row = -1 - bond_col = -1 - if !isempty(left) - bond_row = posInLink!(leftmap, left => lqn) - bond_col = posInLink!(rightmap, vcat(onsite, right) => lqn) - bond_coef = convert(ValType, coefficient(term)) - q_h_sparse = get!(h_sparse, lqn, MatElem{ValType}[]) - push!(q_h_sparse, MatElem(bond_row, bond_col, bond_coef)) - end - - rqn = sqn + lqn - A_row = bond_col - A_col = posInLink!(next_rightmap, right => rqn) - site_coef = one(C) - if A_row == -1 - site_coef = coefficient(term) - end - if isempty(onsite) - if !using_auto_fermion() && isfermionic(right, sites) - push!(onsite, Op("F", n)) - else - push!(onsite, Op("Id", n)) - end - end - el = QNMatElem(lqn, rqn, A_row, A_col, site_coef * Prod(onsite)) - push!(sparse_MPO[n], el) - end - remove_dups!(sparse_MPO[n]) - - if n > 1 && !isempty(h_sparse) - for (q, mat) in h_sparse - h = toMatrix(mat) - U, S, V = svd(h) - P = S .^ 2 - truncate!(P; maxdim, cutoff, mindim) - tdim = length(P) - Vs[n][q] = Matrix{ValType}(V[:, 1:tdim]) - end - end - - rightmap = next_rightmap - next_rightmap = Dict{Pair{Vector{Op},QN},Int}() - end - - # - # Make MPO link indices - # - llinks = Vector{QNIndex}(undef, N + 1) - # Set dir=In for fermionic ordering, avoid arrow sign - # : - linkdir = using_auto_fermion() ? ITensors.In : ITensors.Out - llinks[1] = Index([QN() => 1, Hflux => 1]; tags="Link,l=0", dir=linkdir) - for n in 1:N - qi = Vector{Pair{QN,Int}}() - push!(qi, QN() => 1) - for (q, Vq) in Vs[n + 1] - cols = size(Vq, 2) - if using_auto_fermion() # - push!(qi, (-q) => cols) - else - push!(qi, q => cols) - end - end - push!(qi, Hflux => 1) - llinks[n + 1] = Index(qi...; tags="Link,l=$n", dir=linkdir) - end - - H = MPO(N) - - # Find location where block of Index i - # matches QN q, but *not* 1 or dim(i) - # which are special ending/starting states - function qnblock(i::Index, q::QN) - for b in 2:(nblocks(i) - 1) - flux(i, Block(b)) == q && return b - end - return error("Could not find block of QNIndex with matching QN") - end - qnblockdim(i::Index, q::QN) = blockdim(i, qnblock(i, q)) - - for n in 1:N - ll = llinks[n] - rl = llinks[n + 1] - - begin_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}() - cont_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}() - end_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}() - onsite_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}() - - for el in sparse_MPO[n] - t = el.val - (abs(coefficient(t)) > eps()) || continue - A_row = el.row - A_col = el.col - ct = convert(ValType, coefficient(t)) - - ldim = (A_row == -1) ? 1 : qnblockdim(ll, el.rowqn) - rdim = (A_col == -1) ? 1 : qnblockdim(rl, el.colqn) - zero_mat() = zeros(ValType, ldim, rdim) - - if A_row == -1 && A_col == -1 - # Onsite term - M = get!(onsite_block, (el.rowqn, terms(t)), zeros(ValType, 1, 1)) - M[1, 1] += ct - elseif A_row == -1 - # Operator beginning a term on site n - M = get!(begin_block, (el.rowqn, terms(t)), zero_mat()) - VR = Vs[n + 1][el.colqn] - for c in 1:size(VR, 2) - M[1, c] += ct * VR[A_col, c] - end - elseif A_col == -1 - # Operator ending a term on site n - M = get!(end_block, (el.rowqn, terms(t)), zero_mat()) - VL = Vs[n][el.rowqn] - for r in 1:size(VL, 2) - M[r, 1] += ct * conj(VL[A_row, r]) - end - else - # Operator continuing a term on site n - M = get!(cont_block, (el.rowqn, terms(t)), zero_mat()) - VL = Vs[n][el.rowqn] - VR = Vs[n + 1][el.colqn] - for r in 1:size(VL, 2), c in 1:size(VR, 2) - M[r, c] += ct * conj(VL[A_row, r]) * VR[A_col, c] - end - end - end - - H[n] = ITensor() - - # Helper functions to compute block locations - # of various blocks within the onsite blocks, - # begin blocks, etc. - loc_onsite(rq, cq) = Block(nblocks(ll), 1) - loc_begin(rq, cq) = Block(nblocks(ll), qnblock(rl, cq)) - loc_cont(rq, cq) = Block(qnblock(ll, rq), qnblock(rl, cq)) - loc_end(rq, cq) = Block(qnblock(ll, rq), 1) - - for (loc, block) in ( - (loc_onsite, onsite_block), - (loc_begin, begin_block), - (loc_end, end_block), - (loc_cont, cont_block), - ) - for (q_op, M) in block - op_prod = q_op[2] - Op = computeSiteProd(sites, Prod(op_prod)) - (nnzblocks(Op) == 0) && continue - - rq = q_op[1] - sq = flux(Op) - cq = rq - if !isnothing(sq) - # By convention, if `Op` has no blocks it has a flux - # of `nothing`, catch this case - cq -= sq - - if using_auto_fermion() - # : - # MPO is defined with Index order - # of (rl,s[n]',s[n],cl) where rl = row link, cl = col link - # so compute sign that would result by permuting cl from - # second position to last position: - if fparity(sq) == 1 && fparity(cq) == 1 - Op .*= -1 - end - end - end - - b = loc(rq, cq) - T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl)) - T[b] .= M - - H[n] += (itensor(T) * Op) - end - end - - # Put in ending identity operator - Id = op("Id", sites[n]) - b = Block(1, 1) - T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl)) - T[b] = 1 - H[n] += (itensor(T) * Id) - - # Put in starting identity operator - b = Block(nblocks(ll), nblocks(rl)) - T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl)) - T[b] = 1 - H[n] += (itensor(T) * Id) - end # for n in 1:N - - L = ITensor(llinks[1]) - L[llinks[1] => end] = 1.0 - H[1] *= L - - R = ITensor(dag(llinks[N + 1])) - R[dag(llinks[N + 1]) => 1] = 1.0 - H[N] *= R - - return H -end #qn_svdMPO - -function qn_svdMPO(os::OpSum{C}, sites; kwargs...)::MPO where {C} - # Function barrier to improve type stability - ValType = determineValType(terms(os)) - return qn_svdMPO(ValType, os, sites; kwargs...) -end +## ## using NDTensors: using_auto_fermion +## +## # `ValType::Type{<:Number}` is used instead of `ValType::Type` for efficiency, possibly due to increased method specialization. +## # See https://github.com/ITensor/ITensors.jl/pull/1183. +## function qn_svdMPO( +## ValType::Type{<:Number}, os::OpSum{C}, sites; mindim=1, maxdim=typemax(Int), cutoff=1e-15 +## )::MPO where {C} +## N = length(sites) +## +## # Specifying the element type with `Dict{QN,Matrix{ValType}}[...]` improves type inference and therefore efficiency. +## # See https://github.com/ITensor/ITensors.jl/pull/1183. +## Vs = Dict{QN,Matrix{ValType}}[Dict{QN,Matrix{ValType}}() for n in 1:(N + 1)] +## sparse_MPO = [QNMatElem{Scaled{C,Prod{Op}}}[] for n in 1:N] +## +## function crosses_bond(t::Scaled{C,Prod{Op}}, n::Int) +## return (only(site(t[1])) <= n <= only(site(t[end]))) +## end +## +## # A cache of the ITensor operators on a certain site +## # of a certain type +## op_cache = Dict{Pair{String,Int},ITensor}() +## function calcQN(term::Vector{Op}) +## q = QN() +## for st in term +## op_tensor = get(op_cache, which_op(st) => only(site(st)), nothing) +## if op_tensor === nothing +## op_tensor = op(sites[only(site(st))], which_op(st); params(st)...) +## op_cache[which_op(st) => only(site(st))] = op_tensor +## end +## q -= flux(op_tensor) +## end +## return q +## end +## +## Hflux = -calcQN(terms(first(terms(os)))) +## +## rightmap = Dict{Pair{Vector{Op},QN},Int}() +## next_rightmap = Dict{Pair{Vector{Op},QN},Int}() +## +## for n in 1:N +## h_sparse = Dict{QN,Vector{MatElem{ValType}}}() +## +## leftmap = Dict{Pair{Vector{Op},QN},Int}() +## for term in os +## crosses_bond(term, n) || continue +## +## left = filter(t -> (only(site(t)) < n), terms(term)) +## onsite = filter(t -> (only(site(t)) == n), terms(term)) +## right = filter(t -> (only(site(t)) > n), terms(term)) +## +## lqn = calcQN(left) +## sqn = calcQN(onsite) +## +## bond_row = -1 +## bond_col = -1 +## if !isempty(left) +## bond_row = posInLink!(leftmap, left => lqn) +## bond_col = posInLink!(rightmap, vcat(onsite, right) => lqn) +## bond_coef = convert(ValType, coefficient(term)) +## q_h_sparse = get!(h_sparse, lqn, MatElem{ValType}[]) +## push!(q_h_sparse, MatElem(bond_row, bond_col, bond_coef)) +## end +## +## rqn = sqn + lqn +## A_row = bond_col +## A_col = posInLink!(next_rightmap, right => rqn) +## site_coef = one(C) +## if A_row == -1 +## site_coef = coefficient(term) +## end +## if isempty(onsite) +## if !using_auto_fermion() && isfermionic(right, sites) +## push!(onsite, Op("F", n)) +## else +## push!(onsite, Op("Id", n)) +## end +## end +## el = QNMatElem(lqn, rqn, A_row, A_col, site_coef * Prod(onsite)) +## push!(sparse_MPO[n], el) +## end +## remove_dups!(sparse_MPO[n]) +## +## if n > 1 && !isempty(h_sparse) +## for (q, mat) in h_sparse +## h = toMatrix(mat) +## U, S, V = svd(h) +## P = S .^ 2 +## truncate!(P; maxdim, cutoff, mindim) +## tdim = length(P) +## Vs[n][q] = Matrix{ValType}(V[:, 1:tdim]) +## end +## end +## +## rightmap = next_rightmap +## next_rightmap = Dict{Pair{Vector{Op},QN},Int}() +## end +## +## # +## # Make MPO link indices +## # +## llinks = Vector{QNIndex}(undef, N + 1) +## # Set dir=In for fermionic ordering, avoid arrow sign +## # : +## linkdir = using_auto_fermion() ? ITensors.In : ITensors.Out +## llinks[1] = Index([QN() => 1, Hflux => 1]; tags="Link,l=0", dir=linkdir) +## for n in 1:N +## qi = Vector{Pair{QN,Int}}() +## push!(qi, QN() => 1) +## for (q, Vq) in Vs[n + 1] +## cols = size(Vq, 2) +## if using_auto_fermion() # +## push!(qi, (-q) => cols) +## else +## push!(qi, q => cols) +## end +## end +## push!(qi, Hflux => 1) +## llinks[n + 1] = Index(qi...; tags="Link,l=$n", dir=linkdir) +## end +## +## H = MPO(N) +## +## # Find location where block of Index i +## # matches QN q, but *not* 1 or dim(i) +## # which are special ending/starting states +## function qnblock(i::Index, q::QN) +## for b in 2:(nblocks(i) - 1) +## flux(i, Block(b)) == q && return b +## end +## return error("Could not find block of QNIndex with matching QN") +## end +## qnblockdim(i::Index, q::QN) = blockdim(i, qnblock(i, q)) +## +## for n in 1:N +## ll = llinks[n] +## rl = llinks[n + 1] +## +## begin_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}() +## cont_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}() +## end_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}() +## onsite_block = Dict{Tuple{QN,Vector{Op}},Matrix{ValType}}() +## +## for el in sparse_MPO[n] +## t = el.val +## (abs(coefficient(t)) > eps()) || continue +## A_row = el.row +## A_col = el.col +## ct = convert(ValType, coefficient(t)) +## +## ldim = (A_row == -1) ? 1 : qnblockdim(ll, el.rowqn) +## rdim = (A_col == -1) ? 1 : qnblockdim(rl, el.colqn) +## zero_mat() = zeros(ValType, ldim, rdim) +## +## if A_row == -1 && A_col == -1 +## # Onsite term +## M = get!(onsite_block, (el.rowqn, terms(t)), zeros(ValType, 1, 1)) +## M[1, 1] += ct +## elseif A_row == -1 +## # Operator beginning a term on site n +## M = get!(begin_block, (el.rowqn, terms(t)), zero_mat()) +## VR = Vs[n + 1][el.colqn] +## for c in 1:size(VR, 2) +## M[1, c] += ct * VR[A_col, c] +## end +## elseif A_col == -1 +## # Operator ending a term on site n +## M = get!(end_block, (el.rowqn, terms(t)), zero_mat()) +## VL = Vs[n][el.rowqn] +## for r in 1:size(VL, 2) +## M[r, 1] += ct * conj(VL[A_row, r]) +## end +## else +## # Operator continuing a term on site n +## M = get!(cont_block, (el.rowqn, terms(t)), zero_mat()) +## VL = Vs[n][el.rowqn] +## VR = Vs[n + 1][el.colqn] +## for r in 1:size(VL, 2), c in 1:size(VR, 2) +## M[r, c] += ct * conj(VL[A_row, r]) * VR[A_col, c] +## end +## end +## end +## +## H[n] = ITensor() +## +## # Helper functions to compute block locations +## # of various blocks within the onsite blocks, +## # begin blocks, etc. +## loc_onsite(rq, cq) = Block(nblocks(ll), 1) +## loc_begin(rq, cq) = Block(nblocks(ll), qnblock(rl, cq)) +## loc_cont(rq, cq) = Block(qnblock(ll, rq), qnblock(rl, cq)) +## loc_end(rq, cq) = Block(qnblock(ll, rq), 1) +## +## for (loc, block) in ( +## (loc_onsite, onsite_block), +## (loc_begin, begin_block), +## (loc_end, end_block), +## (loc_cont, cont_block), +## ) +## for (q_op, M) in block +## op_prod = q_op[2] +## Op = computeSiteProd(sites, Prod(op_prod)) +## (nnzblocks(Op) == 0) && continue +## +## rq = q_op[1] +## sq = flux(Op) +## cq = rq +## if !isnothing(sq) +## # By convention, if `Op` has no blocks it has a flux +## # of `nothing`, catch this case +## cq -= sq +## +## if using_auto_fermion() +## # : +## # MPO is defined with Index order +## # of (rl,s[n]',s[n],cl) where rl = row link, cl = col link +## # so compute sign that would result by permuting cl from +## # second position to last position: +## if fparity(sq) == 1 && fparity(cq) == 1 +## Op .*= -1 +## end +## end +## end +## +## b = loc(rq, cq) +## T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl)) +## T[b] .= M +## +## H[n] += (itensor(T) * Op) +## end +## end +## +## # Put in ending identity operator +## Id = op("Id", sites[n]) +## b = Block(1, 1) +## T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl)) +## T[b] = 1 +## H[n] += (itensor(T) * Id) +## +## # Put in starting identity operator +## b = Block(nblocks(ll), nblocks(rl)) +## T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl)) +## T[b] = 1 +## H[n] += (itensor(T) * Id) +## end # for n in 1:N +## +## L = ITensor(llinks[1]) +## L[llinks[1] => end] = 1.0 +## H[1] *= L +## +## R = ITensor(dag(llinks[N + 1])) +## R[dag(llinks[N + 1]) => 1] = 1.0 +## H[N] *= R +## +## return H +## end #qn_svdMPO +## +## function qn_svdMPO(os::OpSum{C}, sites; kwargs...)::MPO where {C} +## # Function barrier to improve type stability +## ValType = determineValType(terms(os)) +## return qn_svdMPO(ValType, os, sites; kwargs...) +## end diff --git a/src/opsum_to_mpo/qnmatelem.jl b/src/opsum_to_mpo/qnmatelem.jl index 7ec55c4..9245187 100644 --- a/src/opsum_to_mpo/qnmatelem.jl +++ b/src/opsum_to_mpo/qnmatelem.jl @@ -1,30 +1,30 @@ -struct QNMatElem{T} - rowqn::QN - colqn::QN - row::Int - col::Int - val::T -end - -function Base.:(==)(m1::QNMatElem{T}, m2::QNMatElem{T})::Bool where {T} - return ( - m1.row == m2.row && - m1.col == m2.col && - m1.val == m2.val && - m1.rowqn == m2.rowqn && - m1.colqn == m2.colqn - ) -end - -function Base.isless(m1::QNMatElem{T}, m2::QNMatElem{T})::Bool where {T} - if m1.rowqn != m2.rowqn - return m1.rowqn < m2.rowqn - elseif m1.colqn != m2.colqn - return m1.colqn < m2.colqn - elseif m1.row != m2.row - return m1.row < m2.row - elseif m1.col != m2.col - return m1.col < m2.col - end - return m1.val < m2.val -end +## struct QNMatElem{T} +## rowqn::QN +## colqn::QN +## row::Int +## col::Int +## val::T +## end +## +## function Base.:(==)(m1::QNMatElem{T}, m2::QNMatElem{T})::Bool where {T} +## return ( +## m1.row == m2.row && +## m1.col == m2.col && +## m1.val == m2.val && +## m1.rowqn == m2.rowqn && +## m1.colqn == m2.colqn +## ) +## end +## +## function Base.isless(m1::QNMatElem{T}, m2::QNMatElem{T})::Bool where {T} +## if m1.rowqn != m2.rowqn +## return m1.rowqn < m2.rowqn +## elseif m1.colqn != m2.colqn +## return m1.colqn < m2.colqn +## elseif m1.row != m2.row +## return m1.row < m2.row +## elseif m1.col != m2.col +## return m1.col < m2.col +## end +## return m1.val < m2.val +## end diff --git a/src/solvers/expand.jl b/src/solvers/expand.jl index c970e45..7291390 100644 --- a/src/solvers/expand.jl +++ b/src/solvers/expand.jl @@ -15,7 +15,7 @@ using ITensors: scalartype, uniqueinds using LinearAlgebra: normalize, svd, tr -using NDTensors: unwrap_array_type +using TypeParameterAccessors: unwrap_array_type # Possible improvements: # - Allow a maxdim argument to be passed to `expand`. diff --git a/test/Project.toml b/test/Project.toml index a8c85de..e176add 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -8,7 +8,6 @@ ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf" Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0" OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" diff --git a/test/base/Project.toml b/test/base/Project.toml index e84dc5a..e42d08c 100644 --- a/test/base/Project.toml +++ b/test/base/Project.toml @@ -1,9 +1,8 @@ [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" +ITensors = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" -ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" From ca7418d2cd0e460caaa58fe736e6a8464ef54205 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Jan 2025 19:51:59 -0500 Subject: [PATCH 02/33] Work towards getting non-QN DMRG working --- examples/dmrg/Project.toml | 2 + src/ITensorMPS.jl | 1 - src/abstractmps.jl | 5 +- src/abstractprojmpo/projmposum.jl | 4 +- src/mpo.jl | 2 +- src/opsum_to_mpo/opsum_to_mpo.jl | 92 ++++++++++++++++++++++++ src/opsum_to_mpo/opsum_to_mpo_generic.jl | 30 ++++++-- src/solvers/timedependentsum.jl | 5 +- 8 files changed, 129 insertions(+), 12 deletions(-) diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index e6b9177..7b391c3 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -1,6 +1,8 @@ [deps] +ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" diff --git a/src/ITensorMPS.jl b/src/ITensorMPS.jl index 30add93..1c42b83 100644 --- a/src/ITensorMPS.jl +++ b/src/ITensorMPS.jl @@ -43,5 +43,4 @@ include("solvers/reducedlinearproblem.jl") include("solvers/linsolve.jl") include("solvers/expand.jl") include("lib/Experimental/src/Experimental.jl") -## include("lib/ITensorMPSNamedDimsArraysExt/src/ITensorMPSNamedDimsArraysExt.jl") end diff --git a/src/abstractmps.jl b/src/abstractmps.jl index 9e823c6..014bdd1 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -2,7 +2,8 @@ using ITensors: @Algorithm_str using IsApprox: Approx, IsApprox -using ITensors: ITensors, ITensor, Index, commonind, commoninds, uniqueind, uniqueinds +using ITensors: + ITensors, ITensor, Index, commonind, commoninds, dag, hasqns, uniqueind, uniqueinds ## using NDTensors: NDTensors, using_auto_fermion, scalartype, tensor using ITensors.Ops: Prod ## using ITensors.QuantumNumbers: QuantumNumbers, removeqn @@ -2293,7 +2294,7 @@ end Return true if the MPS or MPO has tensors which carry quantum numbers. """ -hasqns(M::AbstractMPS) = any(hasqns, data(M)) +ITensors.hasqns(M::AbstractMPS) = any(hasqns, data(M)) # Trait type version of hasqns # Note this is not inferrable, so hasqns would be preferred diff --git a/src/abstractprojmpo/projmposum.jl b/src/abstractprojmpo/projmposum.jl index 3d5d5d9..571c44d 100644 --- a/src/abstractprojmpo/projmposum.jl +++ b/src/abstractprojmpo/projmposum.jl @@ -1,8 +1,8 @@ -using Compat: allequal +using ITensors.Ops: Ops abstract type AbstractSum end -terms(sum::AbstractSum) = sum.terms +Ops.terms(sum::AbstractSum) = sum.terms function set_terms(sum::AbstractSum, terms) return error("Please implement `set_terms` for the `AbstractSum` type `$(typeof(sum))`.") diff --git a/src/mpo.jl b/src/mpo.jl index 4e19724..828ee0b 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -1,7 +1,7 @@ using Adapt: adapt using LinearAlgebra: dot using Random: Random, AbstractRNG -using ITensors: @ts_str, Algorithm, outer +using ITensors: @ts_str, Algorithm, dag, outer using ITensors.Ops: OpSum using ITensors.SiteTypes: SiteTypes, siteind, siteinds diff --git a/src/opsum_to_mpo/opsum_to_mpo.jl b/src/opsum_to_mpo/opsum_to_mpo.jl index 9142e02..15e855c 100644 --- a/src/opsum_to_mpo/opsum_to_mpo.jl +++ b/src/opsum_to_mpo/opsum_to_mpo.jl @@ -1,4 +1,96 @@ ## using NDTensors: using_auto_fermion +using ITensors: dim + +replace_nothing(::Nothing, y) = y +replace_nothing(x, y) = x + +default_mindim(a) = true +default_use_absolute_cutoff(a) = false +default_use_relative_cutoff(a) = true + +function truncate!( + P::AbstractVector; + mindim=nothing, + maxdim=nothing, + cutoff=nothing, + use_absolute_cutoff=nothing, + use_relative_cutoff=nothing, +) + mindim = replace_nothing(mindim, default_mindim(P)) + maxdim = replace_nothing(maxdim, length(P)) + cutoff = replace_nothing(cutoff, typemin(eltype(P))) + use_absolute_cutoff = replace_nothing(use_absolute_cutoff, default_use_absolute_cutoff(P)) + use_relative_cutoff = replace_nothing(use_relative_cutoff, default_use_relative_cutoff(P)) + + origm = length(P) + docut = zero(eltype(P)) + + #if P[1] <= 0.0 + # P[1] = 0.0 + # resize!(P, 1) + # return 0.0, 0.0 + #end + + if origm == 1 + docut = abs(P[1]) / 2 + return zero(eltype(P)), docut + end + + s = sign(P[1]) + s < 0 && (P .*= s) + + #Zero out any negative weight + for n in origm:-1:1 + (P[n] >= zero(eltype(P))) && break + P[n] = zero(eltype(P)) + end + + n = origm + truncerr = zero(eltype(P)) + while n > maxdim + truncerr += P[n] + n -= 1 + end + + if use_absolute_cutoff + #Test if individual prob. weights fall below cutoff + #rather than using *sum* of discarded weights + while P[n] <= cutoff && n > mindim + truncerr += P[n] + n -= 1 + end + else + scale = one(eltype(P)) + if use_relative_cutoff + scale = sum(P) + (scale == zero(eltype(P))) && (scale = one(eltype(P))) + end + + #Continue truncating until *sum* of discarded probability + #weight reaches cutoff reached (or m==mindim) + while (truncerr + P[n] <= cutoff * scale) && (n > mindim) + truncerr += P[n] + n -= 1 + end + + truncerr /= scale + end + + if n < 1 + n = 1 + end + + if n < origm + docut = (P[n] + P[n + 1]) / 2 + if abs(P[n] - P[n + 1]) < eltype(P)(1e-3) * P[n] + docut += eltype(P)(1e-3) * P[n] + end + end + + s < 0 && (P .*= s) + resize!(P, n) + return truncerr, docut +end # `ValType::Type{<:Number}` is used instead of `ValType::Type` for efficiency, possibly due to increased method specialization. # See https://github.com/ITensor/ITensors.jl/pull/1183. diff --git a/src/opsum_to_mpo/opsum_to_mpo_generic.jl b/src/opsum_to_mpo/opsum_to_mpo_generic.jl index 009a699..1d9a0cb 100644 --- a/src/opsum_to_mpo/opsum_to_mpo_generic.jl +++ b/src/opsum_to_mpo/opsum_to_mpo_generic.jl @@ -1,7 +1,28 @@ ## using NDTensors: using_auto_fermion -using ITensors.Ops: Ops, Op, OpSum, Scaled, Sum, coefficient +using ITensors.Ops: + Ops, Op, OpSum, Scaled, Sum, argument, coefficient, site, sites, terms, which_op using ITensors.SiteTypes: has_fermion_string, op +""" + parity_sign(P) + +Given an array or tuple of integers representing +a permutation or a subset of a permutation, +compute the parity sign defined as -1 for a +permutation consisting of an odd number of swaps +and +1 for an even number of swaps. This +implementation uses an O(n^2) algorithm and is +intended for small permutations only. +""" +function parity_sign(P)::Int + L = length(P) + s = +1 + for i in 1:L, j in (i + 1):L + s *= sign(P[j] - P[i]) + end + return s +end + # TODO: Deprecate. const AutoMPO = OpSum @@ -167,7 +188,7 @@ function sorteachterm(os::OpSum, sites) os = copy(os) for (j, t) in enumerate(os) - if maximum(ITensors.sites(t)) > length(sites) + if maximum(Ops.sites(t)) > length(sites) error( "The OpSum contains a term $t that extends beyond the number of sites $(length(sites)).", ) @@ -188,7 +209,8 @@ function sorteachterm(os::OpSum, sites) t_parity = +1 for n in reverse(1:Nt) site_n = only(site(t[n])) - if !using_auto_fermion() && (t_parity == -1) && (site_n < prevsite) + ## if !using_auto_fermion() && (t_parity == -1) && (site_n < prevsite) + if (t_parity == -1) && (site_n < prevsite) # Insert local piece of Jordan-Wigner string emanating # from fermionic operators to the right # (Remaining F operators will be put in by svdMPO) @@ -212,7 +234,7 @@ function sorteachterm(os::OpSum, sites) filter!(!iszero, perm) # and account for anti-commuting, fermionic operators # during above sort; put resulting sign into coef - t *= ITensors.parity_sign(perm) + t *= parity_sign(perm) terms(os)[j] = t end diff --git a/src/solvers/timedependentsum.jl b/src/solvers/timedependentsum.jl index f2998de..d556a7d 100644 --- a/src/solvers/timedependentsum.jl +++ b/src/solvers/timedependentsum.jl @@ -1,4 +1,5 @@ using ITensors: ITensor, inds, permute +using ITensors.Ops: Ops # Represents a time-dependent sum of terms: # @@ -10,7 +11,7 @@ struct TimeDependentSum{Coefficients,Terms} end coefficients(expr::TimeDependentSum) = expr.coefficients -terms(expr::TimeDependentSum) = expr.terms +Ops.terms(expr::TimeDependentSum) = expr.terms function Base.copy(expr::TimeDependentSum) return TimeDependentSum(coefficients(expr), copy.(terms(expr))) end @@ -51,7 +52,7 @@ struct ScaledSum{Coefficients,Terms} end coefficients(expr::ScaledSum) = expr.coefficients -terms(expr::ScaledSum) = expr.terms +Ops.terms(expr::ScaledSum) = expr.terms # Apply the scaled sum of terms: # From 7ed11a49ec1cfef14f516fc28a5b19069f74f014 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sun, 12 Jan 2025 09:42:06 -0500 Subject: [PATCH 03/33] Start fixing MPO constructor --- src/opsum_to_mpo/opsum_to_mpo.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/opsum_to_mpo/opsum_to_mpo.jl b/src/opsum_to_mpo/opsum_to_mpo.jl index 15e855c..31e6836 100644 --- a/src/opsum_to_mpo/opsum_to_mpo.jl +++ b/src/opsum_to_mpo/opsum_to_mpo.jl @@ -1,5 +1,5 @@ ## using NDTensors: using_auto_fermion -using ITensors: dim +using ITensors: dim, inds, itensor replace_nothing(::Nothing, y) = y replace_nothing(x, y) = x @@ -207,6 +207,10 @@ function svdMPO( end T = itensor(M, ll, rl) + + @show inds(T) + @show inds(computeSiteProd(sites, argument(t))) + H[n] += T * computeSiteProd(sites, argument(t)) end From 0ac4690e8f1f1c81ef38975e2a0a52a7a3afbea3 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Thu, 16 Jan 2025 20:27:27 -0500 Subject: [PATCH 04/33] Progress getting DMRG working --- examples/dmrg/1d_heisenberg.jl | 59 +++++++++++++----------- examples/dmrg/Project.toml | 1 - src/abstractmps.jl | 16 +++++-- src/dmrg.jl | 31 +++++++------ src/mps.jl | 27 +++++------ src/opsum_to_mpo/opsum_to_mpo.jl | 3 -- src/opsum_to_mpo/opsum_to_mpo_generic.jl | 7 +-- 7 files changed, 80 insertions(+), 64 deletions(-) diff --git a/examples/dmrg/1d_heisenberg.jl b/examples/dmrg/1d_heisenberg.jl index 4718721..7195bf5 100644 --- a/examples/dmrg/1d_heisenberg.jl +++ b/examples/dmrg/1d_heisenberg.jl @@ -2,37 +2,42 @@ using ITensors, ITensorMPS using Printf using Random -Random.seed!(1234) - -let - N = 100 - - # Create N spin-one degrees of freedom - sites = siteinds("S=1", N) - # Alternatively can make spin-half sites instead - #sites = siteinds("S=1/2", N) - - # Input operator terms which define a Hamiltonian +function heisenberg(N) os = OpSum() for j in 1:(N - 1) os += "Sz", j, "Sz", j + 1 os += 0.5, "S+", j, "S-", j + 1 os += 0.5, "S-", j, "S+", j + 1 end - # Convert these terms to an MPO tensor network - H = MPO(os, sites) - - # Create an initial random matrix product state - psi0 = random_mps(sites; linkdims=10) - - # Plan to do 5 DMRG sweeps: - nsweeps = 5 - # Set maximum MPS bond dimensions for each sweep - maxdim = [10, 20, 100, 100, 200] - # Set maximum truncation error allowed when adapting bond dimensions - cutoff = [1E-11] - - # Run the DMRG algorithm, returning energy and optimized MPS - energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff) - @printf("Final energy = %.12f\n", energy) + return os end + +Random.seed!(1234) + +N = 100 + +# Create N spin-one degrees of freedom +sites = siteinds("S=1", N) +# Alternatively can make spin-half sites instead +#sites = siteinds("S=1/2", N) + +os = heisenberg(N) + +# Input operator terms which define a Hamiltonian +# Convert these terms to an MPO tensor network +H = MPO(os, sites) + +# Create an initial random matrix product state +# psi0 = random_mps(sites; linkdims=10) +psi0 = MPS(sites, j -> isodd(j) ? "↑" : "↓") + +# Plan to do 5 DMRG sweeps: +nsweeps = 5 +# Set maximum MPS bond dimensions for each sweep +maxdim = [10, 20, 100, 100, 200] +# Set maximum truncation error allowed when adapting bond dimensions +cutoff = [1E-11] + +# Run the DMRG algorithm, returning energy and optimized MPS +energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff) +@printf("Final energy = %.12f\n", energy) diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index 7b391c3..68b25a2 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -3,6 +3,5 @@ ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" diff --git a/src/abstractmps.jl b/src/abstractmps.jl index 014bdd1..cf0d7e8 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -3,7 +3,17 @@ using ITensors: @Algorithm_str using IsApprox: Approx, IsApprox using ITensors: - ITensors, ITensor, Index, commonind, commoninds, dag, hasqns, uniqueind, uniqueinds + ITensors, + ITensor, + Index, + commonind, + commoninds, + dag, + factorize, + hasqns, + tags, + uniqueind, + uniqueinds ## using NDTensors: NDTensors, using_auto_fermion, scalartype, tensor using ITensors.Ops: Prod ## using ITensors.QuantumNumbers: QuantumNumbers, removeqn @@ -79,7 +89,7 @@ Base.imag(ψ::AbstractMPS) = imag.(ψ) Base.conj(ψ::AbstractMPS) = conj.(ψ) function convert_leaf_eltype(eltype::Type, ψ::AbstractMPS) - return map(ψᵢ -> convert_leaf_eltype(eltype, ψᵢ), ψ; set_limits=false) + return map(ψᵢ -> eltype.(ψᵢ), ψ; set_limits=false) end """ @@ -760,7 +770,7 @@ function check_hascommoninds(::typeof(siteinds), A::AbstractMPS, B::AbstractMPS) ) end for n in 1:N - !hascommoninds(siteinds(A, n), siteinds(B, n)) && error( + isdisjoint(siteinds(A, n), siteinds(B, n)) && error( "$(typeof(A)) A and $(typeof(B)) B must share site indices. On site $n, A has site indices $(siteinds(A, n)) while B has site indices $(siteinds(B, n)).", ) end diff --git a/src/dmrg.jl b/src/dmrg.jl index 3e0e65c..5114837 100644 --- a/src/dmrg.jl +++ b/src/dmrg.jl @@ -1,29 +1,32 @@ using Adapt: adapt +using ITensors: plev using KrylovKit: eigsolve ## using NDTensors: scalartype, timer using Printf: @printf using TupleTools: TupleTools -function permute( - M::AbstractMPS, ::Tuple{typeof(linkind),typeof(siteinds),typeof(linkind)} -)::typeof(M) - M̃ = typeof(M)(length(M)) - for n in 1:length(M) - lₙ₋₁ = linkind(M, n - 1) - lₙ = linkind(M, n) - s⃗ₙ = TupleTools.sort(Tuple(siteinds(M, n)); by=plev) - M̃[n] = ITensors.permute(M[n], filter(!isnothing, (lₙ₋₁, s⃗ₙ..., lₙ))) - end - set_ortho_lims!(M̃, ortho_lims(M)) - return M̃ -end +## TODO: Add this back. +## function permute( +## M::AbstractMPS, ::Tuple{typeof(linkind),typeof(siteinds),typeof(linkind)} +## )::typeof(M) +## M̃ = typeof(M)(length(M)) +## for n in 1:length(M) +## lₙ₋₁ = linkind(M, n - 1) +## lₙ = linkind(M, n) +## s⃗ₙ = TupleTools.sort(Tuple(siteinds(M, n)); by=plev) +## M̃[n] = ITensors.permute(M[n], filter(!isnothing, (lₙ₋₁, s⃗ₙ..., lₙ))) +## end +## set_ortho_lims!(M̃, ortho_lims(M)) +## return M̃ +## end function dmrg(H::MPO, psi0::MPS, sweeps::Sweeps; kwargs...) check_hascommoninds(siteinds, H, psi0) check_hascommoninds(siteinds, H, psi0') # Permute the indices to have a better memory layout # and minimize permutations - H = permute(H, (linkind, siteinds, linkind)) + ## TODO: Add this back. + ## H = permute(H, (linkind, siteinds, linkind)) PH = ProjMPO(H) return dmrg(PH, psi0, sweeps; kwargs...) end diff --git a/src/mps.jl b/src/mps.jl index a1029a6..6b3112e 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -420,19 +420,20 @@ function MPS(eltype::Type{<:Number}, sites::Vector{<:Index}, states_) states = [state(sites[j], states_[j]) for j in 1:N] - if hasqns(states[1]) - lflux = QN() - for j in 1:(N - 1) - lflux += flux(states[j]) - end - links = Vector{QNIndex}(undef, N - 1) - for j in (N - 1):-1:1 - links[j] = dag(Index(lflux => 1; tags="Link,l=$j")) - lflux -= flux(states[j]) - end - else - links = [Index(1; tags="Link,l=$n") for n in 1:N] - end + ## TODO: Add this back. + ## if hasqns(states[1]) + ## lflux = QN() + ## for j in 1:(N - 1) + ## lflux += flux(states[j]) + ## end + ## links = Vector{QNIndex}(undef, N - 1) + ## for j in (N - 1):-1:1 + ## links[j] = dag(Index(lflux => 1; tags="Link,l=$j")) + ## lflux -= flux(states[j]) + ## end + ## else + links = [Index(1; tags="Link,l=$n") for n in 1:N] + ## end M[1] = ITensor(sites[1], links[1]) M[1] += states[1] * state(links[1], 1) diff --git a/src/opsum_to_mpo/opsum_to_mpo.jl b/src/opsum_to_mpo/opsum_to_mpo.jl index 31e6836..97b7e8e 100644 --- a/src/opsum_to_mpo/opsum_to_mpo.jl +++ b/src/opsum_to_mpo/opsum_to_mpo.jl @@ -208,9 +208,6 @@ function svdMPO( T = itensor(M, ll, rl) - @show inds(T) - @show inds(computeSiteProd(sites, argument(t))) - H[n] += T * computeSiteProd(sites, argument(t)) end diff --git a/src/opsum_to_mpo/opsum_to_mpo_generic.jl b/src/opsum_to_mpo/opsum_to_mpo_generic.jl index 1d9a0cb..dc0275f 100644 --- a/src/opsum_to_mpo/opsum_to_mpo_generic.jl +++ b/src/opsum_to_mpo/opsum_to_mpo_generic.jl @@ -320,9 +320,10 @@ function MPO(os::OpSum, sites::Vector{<:Index}; splitblocks=true, kwargs...)::MP return qn_svdMPO(os, sites; kwargs...) end M = svdMPO(os, sites; kwargs...) - if splitblocks - M = ITensors.splitblocks(linkinds, M) - end + ## TODO: Add this back. + ## if splitblocks + ## M = ITensors.splitblocks(linkinds, M) + ## end return M end From e3d740dbd751cebeec1a37ceada34adfec334ff7 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 17 Jan 2025 11:00:26 -0500 Subject: [PATCH 05/33] More fixes --- Project.toml | 4 +- examples/dmrg/Project.toml | 1 - src/abstractprojmpo/abstractprojmpo.jl | 3 ++ src/dmrg.jl | 68 +++++++++++++------------- src/mps.jl | 40 ++++++++------- 5 files changed, 61 insertions(+), 55 deletions(-) diff --git a/Project.toml b/Project.toml index 5dc149a..6ade918 100644 --- a/Project.toml +++ b/Project.toml @@ -15,6 +15,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SerializedElementArrays = "d3ce8812-9567-47e9-a7b5-65a6d70a3065" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" +VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" [weakdeps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" @@ -37,7 +38,7 @@ Compat = "4.16.0" HDF5 = "0.14, 0.15, 0.16, 0.17" ITensors = "0.8.0" IsApprox = "2.0.0" -KrylovKit = "0.8.1" +KrylovKit = "0.8.1, 0.9" LinearAlgebra = "1.10" Observers = "0.2" PackageCompiler = "1, 2" @@ -47,6 +48,7 @@ SerializedElementArrays = "0.1.0" Test = "1.10" TupleTools = "1.5.0" TypeParameterAccessors = "0.2.1" +VectorInterface = "0.5.0" ZygoteRules = "0.2.2" julia = "1.10" diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index 68b25a2..e6b9177 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -1,5 +1,4 @@ [deps] -ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/abstractprojmpo/abstractprojmpo.jl b/src/abstractprojmpo/abstractprojmpo.jl index 05aae2d..393bbb4 100644 --- a/src/abstractprojmpo/abstractprojmpo.jl +++ b/src/abstractprojmpo/abstractprojmpo.jl @@ -1,3 +1,6 @@ +# TODO: Deprecate, use `ndims` instead. +using ITensors: order + abstract type AbstractProjMPO end copy(::AbstractProjMPO) = error("Not implemented") diff --git a/src/dmrg.jl b/src/dmrg.jl index 5114837..c150caf 100644 --- a/src/dmrg.jl +++ b/src/dmrg.jl @@ -1,32 +1,34 @@ using Adapt: adapt -using ITensors: plev +using ITensors: ITensors, plev using KrylovKit: eigsolve -## using NDTensors: scalartype, timer using Printf: @printf using TupleTools: TupleTools - -## TODO: Add this back. -## function permute( -## M::AbstractMPS, ::Tuple{typeof(linkind),typeof(siteinds),typeof(linkind)} -## )::typeof(M) -## M̃ = typeof(M)(length(M)) -## for n in 1:length(M) -## lₙ₋₁ = linkind(M, n - 1) -## lₙ = linkind(M, n) -## s⃗ₙ = TupleTools.sort(Tuple(siteinds(M, n)); by=plev) -## M̃[n] = ITensors.permute(M[n], filter(!isnothing, (lₙ₋₁, s⃗ₙ..., lₙ))) -## end -## set_ortho_lims!(M̃, ortho_lims(M)) -## return M̃ -## end +using VectorInterface: scalartype + +## TODO: Add this back? +## using NDTensors: timer + +function permute( + M::AbstractMPS, ::Tuple{typeof(linkind),typeof(siteinds),typeof(linkind)} +)::typeof(M) + M̃ = typeof(M)(length(M)) + for n in 1:length(M) + lₙ₋₁ = linkind(M, n - 1) + lₙ = linkind(M, n) + s⃗ₙ = TupleTools.sort(Tuple(siteinds(M, n)); by=plev) + # TODO: Use `NamedDimsArrays.aligndims` instead. + M̃[n] = ITensors.permute(M[n], filter(!isnothing, (lₙ₋₁, s⃗ₙ..., lₙ))) + end + set_ortho_lims!(M̃, ortho_lims(M)) + return M̃ +end function dmrg(H::MPO, psi0::MPS, sweeps::Sweeps; kwargs...) check_hascommoninds(siteinds, H, psi0) check_hascommoninds(siteinds, H, psi0') # Permute the indices to have a better memory layout # and minimize permutations - ## TODO: Add this back. - ## H = permute(H, (linkind, siteinds, linkind)) + H = permute(H, (linkind, siteinds, linkind)) PH = ProjMPO(H) return dmrg(PH, psi0, sweeps; kwargs...) end @@ -221,25 +223,22 @@ function dmrg( end for (b, ha) in sweepnext(N) + ## TODO: Add back `@debug_check`, `checkflux`. ## @debug_check begin ## checkflux(psi) ## checkflux(PH) ## end - ## @timeit_debug timer "dmrg: position!" begin PH = position!(PH, psi, b) - ## end + ## TODO: Add back `@debug_check`, `checkflux`. ## @debug_check begin ## checkflux(psi) ## checkflux(PH) ## end - ## @timeit_debug timer "dmrg: psi[b]*psi[b+1]" begin phi = psi[b] * psi[b + 1] - ## end - ## @timeit_debug timer "dmrg: eigsolve" begin vals, vecs = eigsolve( PH, phi, @@ -250,36 +249,37 @@ function dmrg( krylovdim=eigsolve_krylovdim, maxiter=eigsolve_maxiter, ) - ## end energy = vals[1] + + ## TODO: Bring back some version of this logic. ## Right now there is a conversion problem in CUDA.jl where `UnifiedMemory` Arrays are being converted ## into `DeviceMemory`. This conversion line is here temporarily to fix that problem when it arises ## Adapt is only called when using CUDA backend. CPU will work as implemented previously. ## TODO this might be the only place we really need iscu if its not fixed. - phi = if NDTensors.iscu(phi) && NDTensors.iscu(vecs[1]) - adapt(ITensors.set_eltype(unwrap_array_type(phi), eltype(vecs[1])), vecs[1]) - else - vecs[1] - end + ## phi = if NDTensors.iscu(phi) && NDTensors.iscu(vecs[1]) + ## adapt(ITensors.set_eltype(unwrap_array_type(phi), eltype(vecs[1])), vecs[1]) + ## else + ## vecs[1] + ## end + + phi = vecs[1] ortho = ha == 1 ? "left" : "right" drho = nothing if noise(sweeps, sw) > 0 - ## @timeit_debug timer "dmrg: noiseterm" begin # Use noise term when determining new MPS basis. # This is used to preserve the element type of the MPS. elt = real(scalartype(psi)) drho = elt(noise(sweeps, sw)) * noiseterm(PH, phi, ortho) - ## end end + ## TODO: Add back `@debug_check`, `checkflux`. ## @debug_check begin ## checkflux(phi) ## end - ## @timeit_debug timer "dmrg: replacebond!" begin spec = replacebond!( PH, psi, @@ -294,9 +294,9 @@ function dmrg( which_decomp, svd_alg, ) - ## end maxtruncerr = max(maxtruncerr, spec.truncerr) + ## TODO: Add back `@debug_check`, `checkflux`. ## @debug_check begin ## checkflux(psi) ## checkflux(PH) diff --git a/src/mps.jl b/src/mps.jl index 6b3112e..a5721a3 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -1,7 +1,10 @@ using Adapt: adapt -## using NDTensors: using_auto_fermion -using Random: Random, AbstractRNG +using ITensors: hasqns using ITensors.SiteTypes: SiteTypes, siteind, siteinds, state +using Random: Random, AbstractRNG + +## TODO: Add this back. +## using NDTensors: using_auto_fermion """ MPS @@ -420,20 +423,19 @@ function MPS(eltype::Type{<:Number}, sites::Vector{<:Index}, states_) states = [state(sites[j], states_[j]) for j in 1:N] - ## TODO: Add this back. - ## if hasqns(states[1]) - ## lflux = QN() - ## for j in 1:(N - 1) - ## lflux += flux(states[j]) - ## end - ## links = Vector{QNIndex}(undef, N - 1) - ## for j in (N - 1):-1:1 - ## links[j] = dag(Index(lflux => 1; tags="Link,l=$j")) - ## lflux -= flux(states[j]) - ## end - ## else - links = [Index(1; tags="Link,l=$n") for n in 1:N] - ## end + if hasqns(states[1]) + lflux = QN() + for j in 1:(N - 1) + lflux += flux(states[j]) + end + links = Vector{QNIndex}(undef, N - 1) + for j in (N - 1):-1:1 + links[j] = dag(Index(lflux => 1; tags="Link,l=$j")) + lflux -= flux(states[j]) + end + else + links = [Index(1; tags="Link,l=$n") for n in 1:N] + end M[1] = ITensor(sites[1], links[1]) M[1] += states[1] * state(links[1], 1) @@ -551,9 +553,9 @@ function replacebond!( use_relative_cutoff=nothing, min_blockdim=nothing, ) - normalize = NDTensors.replace_nothing(normalize, false) - swapsites = NDTensors.replace_nothing(swapsites, false) - ortho = NDTensors.replace_nothing(ortho, "left") + normalize = replace_nothing(normalize, false) + swapsites = replace_nothing(swapsites, false) + ortho = replace_nothing(ortho, "left") indsMb = inds(M[b]) if swapsites From 01bb8ae742e70963da0305e878e736fbc1fb7df2 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 17 Jan 2025 17:28:06 -0500 Subject: [PATCH 06/33] Overload LinearAlgebra.norm --- src/abstractmps.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/abstractmps.jl b/src/abstractmps.jl index cf0d7e8..ed217d9 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -1244,7 +1244,7 @@ the full inner product of the MPS/MPO with itself. See also [`lognorm`](@ref). """ -function norm(M::AbstractMPS) +function LinearAlgebra.norm(M::AbstractMPS) if isortho(M) return norm(M[orthocenter(M)]) end From 924854c8fa3d31181228eb3fca8dd33e813af6fb Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 18 Jan 2025 13:29:41 -0500 Subject: [PATCH 07/33] Fix inner --- examples/dmrg/1d_heisenberg.jl | 2 +- src/abstractmps.jl | 32 +++++++++++++++++--------------- src/mpo.jl | 2 +- 3 files changed, 19 insertions(+), 17 deletions(-) diff --git a/examples/dmrg/1d_heisenberg.jl b/examples/dmrg/1d_heisenberg.jl index 7195bf5..6e63465 100644 --- a/examples/dmrg/1d_heisenberg.jl +++ b/examples/dmrg/1d_heisenberg.jl @@ -14,7 +14,7 @@ end Random.seed!(1234) -N = 100 +N = 10 # Create N spin-one degrees of freedom sites = siteinds("S=1", N) diff --git a/src/abstractmps.jl b/src/abstractmps.jl index ed217d9..7c0e9f9 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -11,6 +11,8 @@ using ITensors: dag, factorize, hasqns, + replaceinds, + sim, tags, uniqueind, uniqueinds @@ -702,12 +704,12 @@ function replaceinds!(::typeof(linkinds), M::AbstractMPS, l̃s::Vector{<:Index}) return M end -function replaceinds(::typeof(linkinds), M::AbstractMPS, l̃s::Vector{<:Index}) +function ITensors.replaceinds(::typeof(linkinds), M::AbstractMPS, l̃s::Vector{<:Index}) return replaceinds!(linkinds, copy(M), l̃s) end # TODO: change kwarg from `set_limits` to `preserve_ortho` -function map!(f::Function, M::AbstractMPS; set_limits::Bool=true) +function Base.map!(f::Function, M::AbstractMPS; set_limits::Bool=true) for i in eachindex(M) M[i, set_limits=set_limits] = f(M[i]) end @@ -715,7 +717,7 @@ function map!(f::Function, M::AbstractMPS; set_limits::Bool=true) end # TODO: change kwarg from `set_limits` to `preserve_ortho` -function map(f::Function, M::AbstractMPS; set_limits::Bool=true) +function Base.map(f::Function, M::AbstractMPS; set_limits::Bool=true) return map!(f, copy(M); set_limits=set_limits) end @@ -777,44 +779,44 @@ function check_hascommoninds(::typeof(siteinds), A::AbstractMPS, B::AbstractMPS) return nothing end -function map!(f::Function, ::typeof(linkinds), M::AbstractMPS) +function Base.map!(f::Function, ::typeof(linkinds), M::AbstractMPS) for i in eachindex(M)[1:(end - 1)] l = linkinds(M, i) if !isempty(l) - l̃ = f(l) + l̃ = f.(l) @preserve_ortho M begin - M[i] = replaceinds(M[i], l, l̃) - M[i + 1] = replaceinds(M[i + 1], l, l̃) + M[i] = replaceinds(M[i], (l .=> l̃)...) + M[i + 1] = replaceinds(M[i + 1], (l .=> l̃)...) end end end return M end -map(f::Function, ::typeof(linkinds), M::AbstractMPS) = map!(f, linkinds, copy(M)) +Base.map(f::Function, ::typeof(linkinds), M::AbstractMPS) = map!(f, linkinds, copy(M)) -function map!(f::Function, ::typeof(siteinds), M::AbstractMPS) +function Base.map!(f::Function, ::typeof(siteinds), M::AbstractMPS) for i in eachindex(M) s = siteinds(M, i) if !isempty(s) @preserve_ortho M begin - M[i] = replaceinds(M[i], s, f(s)) + M[i] = replaceinds(M[i], s, f.(s)) end end end return M end -map(f::Function, ::typeof(siteinds), M::AbstractMPS) = map!(f, siteinds, copy(M)) +Base.map(f::Function, ::typeof(siteinds), M::AbstractMPS) = map!(f, siteinds, copy(M)) -function map!( +function Base.map!( f::Function, ::typeof(siteinds), ::typeof(commoninds), M1::AbstractMPS, M2::AbstractMPS ) length(M1) != length(M2) && error("MPOs/MPSs must be the same length") for i in eachindex(M1) s = siteinds(commoninds, M1, M2, i) if !isempty(s) - s̃ = f(s) + s̃ = f.(s) @preserve_ortho (M1, M2) begin M1[i] = replaceinds(M1[i], s .=> s̃) M2[i] = replaceinds(M2[i], s .=> s̃) @@ -824,7 +826,7 @@ function map!( return M1, M2 end -function map!( +function Base.map!( f::Function, ::typeof(siteinds), ::typeof(uniqueinds), M1::AbstractMPS, M2::AbstractMPS ) length(M1) != length(M2) && error("MPOs/MPSs must be the same length") @@ -839,7 +841,7 @@ function map!( return M1 end -function map( +function Base.map( f::Function, ffilter::typeof(siteinds), fsubset::Union{typeof(commoninds),typeof(uniqueinds)}, diff --git a/src/mpo.jl b/src/mpo.jl index 828ee0b..23dc5c9 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -283,7 +283,7 @@ function hassameinds(::typeof(siteinds), ψ::MPS, Hϕ::Tuple{MPO,MPS}) N = length(ψ) @assert N == length(Hϕ[1]) == length(Hϕ[1]) for n in 1:N - !hassameinds(siteinds(Hϕ, n), siteinds(ψ, n)) && return false + !issetequal(siteinds(Hϕ, n), siteinds(ψ, n)) && return false end return true end From a16f115711b5764dfeffbb8ba4747784e81a7122 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 18 Jan 2025 16:59:22 -0500 Subject: [PATCH 08/33] More fixes for inner, dmrg --- examples/dmrg/1d_heisenberg.jl | 4 +++- src/abstractmps.jl | 4 ++-- src/dmrg.jl | 1 + 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/examples/dmrg/1d_heisenberg.jl b/examples/dmrg/1d_heisenberg.jl index 6e63465..11938be 100644 --- a/examples/dmrg/1d_heisenberg.jl +++ b/examples/dmrg/1d_heisenberg.jl @@ -34,10 +34,12 @@ psi0 = MPS(sites, j -> isodd(j) ? "↑" : "↓") # Plan to do 5 DMRG sweeps: nsweeps = 5 # Set maximum MPS bond dimensions for each sweep -maxdim = [10, 20, 100, 100, 200] +maxdim = [10] # Set maximum truncation error allowed when adapting bond dimensions cutoff = [1E-11] # Run the DMRG algorithm, returning energy and optimized MPS energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff) +@show inner(psi', H, psi) +@show inner(psi, psi) @printf("Final energy = %.12f\n", energy) diff --git a/src/abstractmps.jl b/src/abstractmps.jl index 7c0e9f9..a615e1f 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -854,7 +854,7 @@ end function hassameinds(::typeof(siteinds), M1::AbstractMPS, M2::AbstractMPS) length(M1) ≠ length(M2) && return false for n in 1:length(M1) - !hassameinds(siteinds(all, M1, n), siteinds(all, M2, n)) && return false + !issetequal(siteinds(all, M1, n), siteinds(all, M2, n)) && return false end return true end @@ -1106,7 +1106,7 @@ function deprecate_make_inds_match!( make_inds_match = false end if !hassameinds(siteinds, M1dag, M2) && make_inds_match - ITensors.warn_once(inner_mps_mpo_mps_deprecation_warning(), :inner_mps_mps) + #ITensors.warn_once(inner_mps_mpo_mps_deprecation_warning(), :inner_mps_mps) replace_siteinds!(M1dag, siteindsM2) end return M1dag, M2 diff --git a/src/dmrg.jl b/src/dmrg.jl index c150caf..e881e49 100644 --- a/src/dmrg.jl +++ b/src/dmrg.jl @@ -248,6 +248,7 @@ function dmrg( tol=eigsolve_tol, krylovdim=eigsolve_krylovdim, maxiter=eigsolve_maxiter, + verbosity=eigsolve_verbosity, ) energy = vals[1] From e3232655fb60564df3ef0f41fe1c3b36a8eb777d Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sun, 19 Jan 2025 21:19:57 -0500 Subject: [PATCH 09/33] Make use of new quantum operator libraries --- Project.toml | 4 +++ src/abstractmps.jl | 34 +++++++++++++++--------- src/abstractprojmpo/projmposum.jl | 4 +-- src/exports.jl | 1 - src/imports.jl | 6 ++--- src/mpo.jl | 14 +++++----- src/mps.jl | 13 ++++++--- src/opsum_to_mpo/opsum_to_mpo_generic.jl | 24 ++++++++++++----- src/solvers/timedependentsum.jl | 6 ++--- 9 files changed, 68 insertions(+), 38 deletions(-) diff --git a/Project.toml b/Project.toml index 6ade918..40c66a0 100644 --- a/Project.toml +++ b/Project.toml @@ -6,11 +6,13 @@ version = "0.4.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" +ITensorQuantumOperatorDefinitions = "fd9b415b-e710-4e2a-b407-cba023081494" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" IsApprox = "28f27b66-4bd8-47e7-9110-e2746eb8bed7" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +QuantumOperatorAlgebra = "c7a6d0f7-daa6-4368-ba67-89ed64127c3b" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SerializedElementArrays = "d3ce8812-9567-47e9-a7b5-65a6d70a3065" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" @@ -36,6 +38,7 @@ Adapt = "4.1.0" ChainRulesCore = "1.10" Compat = "4.16.0" HDF5 = "0.14, 0.15, 0.16, 0.17" +ITensorQuantumOperatorDefinitions = "0.1.1" ITensors = "0.8.0" IsApprox = "2.0.0" KrylovKit = "0.8.1, 0.9" @@ -43,6 +46,7 @@ LinearAlgebra = "1.10" Observers = "0.2" PackageCompiler = "1, 2" Printf = "1.10" +QuantumOperatorAlgebra = "0.1.0" Random = "1.10" SerializedElementArrays = "0.1.0" Test = "1.10" diff --git a/src/abstractmps.jl b/src/abstractmps.jl index a615e1f..9eccdac 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -17,9 +17,9 @@ using ITensors: uniqueind, uniqueinds ## using NDTensors: NDTensors, using_auto_fermion, scalartype, tensor -using ITensors.Ops: Prod +using QuantumOperatorAlgebra: Prod ## using ITensors.QuantumNumbers: QuantumNumbers, removeqn -using ITensors.SiteTypes: SiteTypes, siteinds +using ITensorQuantumOperatorDefinitions: ITensorQuantumOperatorDefinitions, siteinds ## using ITensors.TagSets: TagSets using LinearAlgebra: LinearAlgebra @@ -485,7 +485,7 @@ end Get the site index (or indices) of MPO `A` that is unique to `A` (not shared with MPS/MPO `B`). """ -function SiteTypes.siteinds( +function ITensorQuantumOperatorDefinitions.siteinds( f::Union{typeof(uniqueinds),typeof(uniqueind)}, A::AbstractMPS, B::AbstractMPS, @@ -505,7 +505,7 @@ end Get the site indices of MPO `A` that are unique to `A` (not shared with MPS/MPO `B`), as a `Vector{<:Index}`. """ -function SiteTypes.siteinds( +function ITensorQuantumOperatorDefinitions.siteinds( f::Union{typeof(uniqueinds),typeof(uniqueind)}, A::AbstractMPS, B::AbstractMPS; kwargs... ) return [siteinds(f, A, B, j; kwargs...) for j in eachindex(A)] @@ -517,7 +517,7 @@ end Get the site index (or indices) of the `j`th MPO tensor of `A` that is shared with MPS/MPO `B`. """ -function SiteTypes.siteinds( +function ITensorQuantumOperatorDefinitions.siteinds( f::Union{typeof(commoninds),typeof(commonind)}, A::AbstractMPS, B::AbstractMPS, @@ -533,7 +533,7 @@ end Get a vector of the site index (or indices) of MPO `A` that is shared with MPS/MPO `B`. """ -function SiteTypes.siteinds( +function ITensorQuantumOperatorDefinitions.siteinds( f::Union{typeof(commoninds),typeof(commonind)}, A::AbstractMPS, B::AbstractMPS; kwargs... ) return [siteinds(f, A, B, j) for j in eachindex(A)] @@ -639,7 +639,9 @@ Return the first site Index found on the MPS or MPO You can choose different filters, like prime level and tags, with the `kwargs`. """ -function SiteTypes.siteind(::typeof(first), M::AbstractMPS, j::Integer; kwargs...) +function ITensorQuantumOperatorDefinitions.siteind( + ::typeof(first), M::AbstractMPS, j::Integer; kwargs... +) N = length(M) (N == 1) && return firstind(M[1]; kwargs...) if j == 1 @@ -661,7 +663,7 @@ at the site `j` as an IndexSet. Optionally filter prime tags and prime levels with keyword arguments like `plev` and `tags`. """ -function SiteTypes.siteinds(M::AbstractMPS, j::Integer; kwargs...) +function ITensorQuantumOperatorDefinitions.siteinds(M::AbstractMPS, j::Integer; kwargs...) N = length(M) (N == 1) && return inds(M[1]; kwargs...) if j == 1 @@ -674,19 +676,27 @@ function SiteTypes.siteinds(M::AbstractMPS, j::Integer; kwargs...) return si end -function SiteTypes.siteinds(::typeof(all), ψ::AbstractMPS, n::Integer; kwargs...) +function ITensorQuantumOperatorDefinitions.siteinds( + ::typeof(all), ψ::AbstractMPS, n::Integer; kwargs... +) return siteinds(ψ, n; kwargs...) end -function SiteTypes.siteinds(::typeof(first), ψ::AbstractMPS; kwargs...) +function ITensorQuantumOperatorDefinitions.siteinds( + ::typeof(first), ψ::AbstractMPS; kwargs... +) return [siteind(first, ψ, j; kwargs...) for j in 1:length(ψ)] end -function SiteTypes.siteinds(::typeof(only), ψ::AbstractMPS; kwargs...) +function ITensorQuantumOperatorDefinitions.siteinds( + ::typeof(only), ψ::AbstractMPS; kwargs... +) return [siteind(only, ψ, j; kwargs...) for j in 1:length(ψ)] end -function SiteTypes.siteinds(::typeof(all), ψ::AbstractMPS; kwargs...) +function ITensorQuantumOperatorDefinitions.siteinds( + ::typeof(all), ψ::AbstractMPS; kwargs... +) return [siteinds(ψ, j; kwargs...) for j in 1:length(ψ)] end diff --git a/src/abstractprojmpo/projmposum.jl b/src/abstractprojmpo/projmposum.jl index 571c44d..2b3372b 100644 --- a/src/abstractprojmpo/projmposum.jl +++ b/src/abstractprojmpo/projmposum.jl @@ -1,8 +1,8 @@ -using ITensors.Ops: Ops +using QuantumOperatorAlgebra: QuantumOperatorAlgebra abstract type AbstractSum end -Ops.terms(sum::AbstractSum) = sum.terms +QuantumOperatorAlgebra.terms(sum::AbstractSum) = sum.terms function set_terms(sum::AbstractSum, terms) return error("Please implement `set_terms` for the `AbstractSum` type `$(typeof(sum))`.") diff --git a/src/exports.jl b/src/exports.jl index 749a73c..a620fd4 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -10,7 +10,6 @@ export Apply, Op, OpName, - Ops, Prod, Scaled, SiteType, diff --git a/src/imports.jl b/src/imports.jl index e72aba8..6779a7b 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -1,7 +1,7 @@ # Primarily used to import names into the `ITensorMPS` # module from submodules or from `ITensors` so they can # be reexported. -using ITensors.SiteTypes: +using ITensorQuantumOperatorDefinitions: @OpName_str, @SiteType_str, @StateName_str, @@ -13,7 +13,7 @@ using ITensors.SiteTypes: TagType, ValName, ops -using ITensors.Ops: Trotter +using QuantumOperatorAlgebra: Trotter import Base: # types @@ -168,6 +168,6 @@ import ..ITensors: truncate!, which_op -import ..ITensors.Ops: params +import QuantumOperatorAlgebra: params import SerializedElementArrays: disk diff --git a/src/mpo.jl b/src/mpo.jl index 23dc5c9..fe476f7 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -2,8 +2,9 @@ using Adapt: adapt using LinearAlgebra: dot using Random: Random, AbstractRNG using ITensors: @ts_str, Algorithm, dag, outer -using ITensors.Ops: OpSum -using ITensors.SiteTypes: SiteTypes, siteind, siteinds +using QuantumOperatorAlgebra: OpSum +using ITensorQuantumOperatorDefinitions: + ITensorQuantumOperatorDefinitions, siteind, siteinds """ MPO @@ -242,7 +243,8 @@ end Get the first site Index of the MPO found, by default with prime level 0. """ -SiteTypes.siteind(M::MPO, j::Int; kwargs...) = siteind(first, M, j; plev=0, kwargs...) +ITensorQuantumOperatorDefinitions.siteind(M::MPO, j::Int; kwargs...) = + siteind(first, M, j; plev=0, kwargs...) # TODO: make this return the site indices that would have # been used to create the MPO? I.e.: @@ -252,9 +254,9 @@ SiteTypes.siteind(M::MPO, j::Int; kwargs...) = siteind(first, M, j; plev=0, kwar Get a Vector of IndexSets of all the site indices of M. """ -SiteTypes.siteinds(M::MPO; kwargs...) = siteinds(all, M; kwargs...) +ITensorQuantumOperatorDefinitions.siteinds(M::MPO; kwargs...) = siteinds(all, M; kwargs...) -function SiteTypes.siteinds(Mψ::Tuple{MPO,MPS}, n::Int; kwargs...) +function ITensorQuantumOperatorDefinitions.siteinds(Mψ::Tuple{MPO,MPS}, n::Int; kwargs...) return siteinds(uniqueinds, Mψ[1], Mψ[2], n; kwargs...) end @@ -265,7 +267,7 @@ function nsites(Mψ::Tuple{MPO,MPS}) return N end -function SiteTypes.siteinds(Mψ::Tuple{MPO,MPS}; kwargs...) +function ITensorQuantumOperatorDefinitions.siteinds(Mψ::Tuple{MPO,MPS}; kwargs...) return [siteinds(Mψ, n; kwargs...) for n in 1:nsites(Mψ)] end diff --git a/src/mps.jl b/src/mps.jl index a5721a3..a56852c 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -1,6 +1,7 @@ using Adapt: adapt using ITensors: hasqns -using ITensors.SiteTypes: SiteTypes, siteind, siteinds, state +using ITensorQuantumOperatorDefinitions: + ITensorQuantumOperatorDefinitions, siteind, siteinds, state using Random: Random, AbstractRNG ## TODO: Add this back. @@ -486,14 +487,17 @@ MPS(sites::Vector{<:Index}, states) = MPS(Float64, sites, states) Get the first site Index of the MPS. Return `nothing` if none is found. """ -SiteTypes.siteind(M::MPS, j::Int; kwargs...) = siteind(first, M, j; kwargs...) +ITensorQuantumOperatorDefinitions.siteind(M::MPS, j::Int; kwargs...) = + siteind(first, M, j; kwargs...) """ siteind(::typeof(only), M::MPS, j::Int; kwargs...) Get the only site Index of the MPS. Return `nothing` if none is found. """ -function SiteTypes.siteind(::typeof(only), M::MPS, j::Int; kwargs...) +function ITensorQuantumOperatorDefinitions.siteind( + ::typeof(only), M::MPS, j::Int; kwargs... +) is = siteinds(M, j; kwargs...) if isempty(is) return nothing @@ -515,7 +519,8 @@ Get a vector of the only site Index found on each tensor of the MPS. Errors if m Get a vector of the all site Indices found on each tensor of the MPS. Returns a Vector of IndexSets. """ -SiteTypes.siteinds(M::MPS; kwargs...) = siteinds(first, M; kwargs...) +ITensorQuantumOperatorDefinitions.siteinds(M::MPS; kwargs...) = + siteinds(first, M; kwargs...) function replace_siteinds!(M::MPS, sites) for j in eachindex(M) diff --git a/src/opsum_to_mpo/opsum_to_mpo_generic.jl b/src/opsum_to_mpo/opsum_to_mpo_generic.jl index dc0275f..bf776b9 100644 --- a/src/opsum_to_mpo/opsum_to_mpo_generic.jl +++ b/src/opsum_to_mpo/opsum_to_mpo_generic.jl @@ -1,7 +1,17 @@ ## using NDTensors: using_auto_fermion -using ITensors.Ops: - Ops, Op, OpSum, Scaled, Sum, argument, coefficient, site, sites, terms, which_op -using ITensors.SiteTypes: has_fermion_string, op +using QuantumOperatorAlgebra: + QuantumOperatorAlgebra, + Op, + OpSum, + Scaled, + Sum, + argument, + coefficient, + site, + sites, + terms, + which_op +using ITensorQuantumOperatorDefinitions: has_fermion_string, op """ parity_sign(P) @@ -85,10 +95,10 @@ end add!(os::OpSum, o::Op) = add!(os, Prod{Op}() * o) add!(os::OpSum, o::Scaled{C,Op}) where {C} = add!(os, Prod{Op}() * o) add!(os::OpSum, o::Prod{Op}) = add!(os, one(Float64) * o) -add!(os::OpSum, o::Tuple) = add!(os, Ops.op_term(o)) +add!(os::OpSum, o::Tuple) = add!(os, QuantumOperatorAlgebra.op_term(o)) add!(os::OpSum, a1::String, args...) = add!(os, (a1, args...)) add!(os::OpSum, a1::Number, args...) = add!(os, (a1, args...)) -subtract!(os::OpSum, o::Tuple) = add!(os, -Ops.op_term(o)) +subtract!(os::OpSum, o::Tuple) = add!(os, -QuantumOperatorAlgebra.op_term(o)) function isfermionic(t::Vector{Op}, sites) p = +1 @@ -188,7 +198,7 @@ function sorteachterm(os::OpSum, sites) os = copy(os) for (j, t) in enumerate(os) - if maximum(Ops.sites(t)) > length(sites) + if maximum(QuantumOperatorAlgebra.sites(t)) > length(sites) error( "The OpSum contains a term $t that extends beyond the number of sites $(length(sites)).", ) @@ -362,7 +372,7 @@ function MPO( return MPO(eltype, OpSum{C}() + o, s; kwargs...) end -# Like `Ops.OpSumLike` but without `OpSum` included. +# Like `QuantumOperatorAlgebra.OpSumLike` but without `OpSum` included. const OpSumLikeWithoutOpSum{C} = Union{ Op,Scaled{C,Op},Sum{Op},Prod{Op},Scaled{C,Prod{Op}},Sum{Scaled{C,Op}} } diff --git a/src/solvers/timedependentsum.jl b/src/solvers/timedependentsum.jl index d556a7d..9acf431 100644 --- a/src/solvers/timedependentsum.jl +++ b/src/solvers/timedependentsum.jl @@ -1,5 +1,5 @@ using ITensors: ITensor, inds, permute -using ITensors.Ops: Ops +using QuantumOperatorAlgebra: QuantumOperatorAlgebra # Represents a time-dependent sum of terms: # @@ -11,7 +11,7 @@ struct TimeDependentSum{Coefficients,Terms} end coefficients(expr::TimeDependentSum) = expr.coefficients -Ops.terms(expr::TimeDependentSum) = expr.terms +QuantumOperatorAlgebra.terms(expr::TimeDependentSum) = expr.terms function Base.copy(expr::TimeDependentSum) return TimeDependentSum(coefficients(expr), copy.(terms(expr))) end @@ -52,7 +52,7 @@ struct ScaledSum{Coefficients,Terms} end coefficients(expr::ScaledSum) = expr.coefficients -Ops.terms(expr::ScaledSum) = expr.terms +QuantumOperatorAlgebra.terms(expr::ScaledSum) = expr.terms # Apply the scaled sum of terms: # From f550c6d85a06574b59ba32f053cb7db9d7e04957 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sun, 19 Jan 2025 22:22:09 -0500 Subject: [PATCH 10/33] Fix op call --- src/opsum_to_mpo/opsum_to_mpo_generic.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/opsum_to_mpo/opsum_to_mpo_generic.jl b/src/opsum_to_mpo/opsum_to_mpo_generic.jl index bf776b9..79a5f5a 100644 --- a/src/opsum_to_mpo/opsum_to_mpo_generic.jl +++ b/src/opsum_to_mpo/opsum_to_mpo_generic.jl @@ -165,10 +165,10 @@ end function computeSiteProd(sites, ops::Prod{Op})::ITensor i = only(site(ops[1])) - T = op(sites[i], which_op(ops[1]); params(ops[1])...) + T = op(which_op(ops[1]), sites[i]; params(ops[1])...) for j in 2:length(ops) (only(site(ops[j])) != i) && error("Mismatch of site number in computeSiteProd") - opj = op(sites[i], which_op(ops[j]); params(ops[j])...) + opj = op(which_op(ops[j]), sites[i]; params(ops[j])...) T = product(T, opj) end return T From 6350e39633b376195432e4f3efc22acba0133e9d Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 25 Jan 2025 21:26:41 -0500 Subject: [PATCH 11/33] Simplifications --- Project.toml | 6 +- examples/dmrg/Project.toml | 1 + src/ITensorMPS.jl | 1 + src/abstractmps.jl | 65 +++++----- src/abstractprojmpo/abstractprojmpo.jl | 5 +- src/deprecated.jl | 1 - src/imports.jl | 114 ++++++++---------- .../src/ITensorMPSNamedDimsArraysExt.jl | 3 - .../src/to_nameddimsarray.jl | 6 - src/mpo.jl | 28 +++-- src/mps.jl | 27 ++--- src/opsum_to_mpo/opsum_to_mpo.jl | 4 +- src/opsum_to_mpo/opsum_to_mpo_generic.jl | 2 +- src/siteinds.jl | 8 ++ src/solvers/ITensorsExtensions.jl | 10 +- src/solvers/alternating_update.jl | 5 +- src/solvers/dmrg.jl | 3 +- src/solvers/dmrg_x.jl | 5 +- src/solvers/expand.jl | 2 +- src/solvers/sweep_update.jl | 3 +- src/solvers/timedependentsum.jl | 5 +- 21 files changed, 149 insertions(+), 155 deletions(-) delete mode 100644 src/lib/ITensorMPSNamedDimsArraysExt/src/ITensorMPSNamedDimsArraysExt.jl delete mode 100644 src/lib/ITensorMPSNamedDimsArraysExt/src/to_nameddimsarray.jl create mode 100644 src/siteinds.jl diff --git a/Project.toml b/Project.toml index 40c66a0..7185dc6 100644 --- a/Project.toml +++ b/Project.toml @@ -6,13 +6,14 @@ version = "0.4.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" -ITensorQuantumOperatorDefinitions = "fd9b415b-e710-4e2a-b407-cba023081494" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" IsApprox = "28f27b66-4bd8-47e7-9110-e2746eb8bed7" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuantumOperatorAlgebra = "c7a6d0f7-daa6-4368-ba67-89ed64127c3b" +QuantumOperatorDefinitions = "826dd319-6fd5-459a-a990-3a4f214664bf" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SerializedElementArrays = "d3ce8812-9567-47e9-a7b5-65a6d70a3065" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" @@ -38,15 +39,16 @@ Adapt = "4.1.0" ChainRulesCore = "1.10" Compat = "4.16.0" HDF5 = "0.14, 0.15, 0.16, 0.17" -ITensorQuantumOperatorDefinitions = "0.1.1" ITensors = "0.8.0" IsApprox = "2.0.0" KrylovKit = "0.8.1, 0.9" LinearAlgebra = "1.10" +NamedDimsArrays = "0.3.12" Observers = "0.2" PackageCompiler = "1, 2" Printf = "1.10" QuantumOperatorAlgebra = "0.1.0" +QuantumOperatorDefinitions = "0.1.2" Random = "1.10" SerializedElementArrays = "0.1.0" Test = "1.10" diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index e6b9177..c88c435 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -2,5 +2,6 @@ ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +QuantumOperatorDefinitions = "826dd319-6fd5-459a-a990-3a4f214664bf" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" diff --git a/src/ITensorMPS.jl b/src/ITensorMPS.jl index 1c42b83..288085e 100644 --- a/src/ITensorMPS.jl +++ b/src/ITensorMPS.jl @@ -2,6 +2,7 @@ module ITensorMPS using ITensors include("imports.jl") include("exports.jl") +include("siteinds.jl") include("abstractmps.jl") include("mps.jl") include("mpo.jl") diff --git a/src/abstractmps.jl b/src/abstractmps.jl index 9eccdac..c18ab7e 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -16,12 +16,11 @@ using ITensors: tags, uniqueind, uniqueinds -## using NDTensors: NDTensors, using_auto_fermion, scalartype, tensor using QuantumOperatorAlgebra: Prod ## using ITensors.QuantumNumbers: QuantumNumbers, removeqn -using ITensorQuantumOperatorDefinitions: ITensorQuantumOperatorDefinitions, siteinds ## using ITensors.TagSets: TagSets using LinearAlgebra: LinearAlgebra +using VectorInterface: VectorInterface, scalartype abstract type AbstractMPS end @@ -70,9 +69,11 @@ have type `ComplexF64`, return `ComplexF64`. """ promote_itensor_eltype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m) -## NDTensors.scalartype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m) -## NDTensors.scalartype(m::Array{ITensor}) = LinearAlgebra.promote_leaf_eltypes(m) -## NDTensors.scalartype(m::Array{<:Array{ITensor}}) = LinearAlgebra.promote_leaf_eltypes(m) +VectorInterface.scalartype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m) + +## TODO: Are these definitions needed? +## VectorInterface.scalartype(m::Array{ITensor}) = LinearAlgebra.promote_leaf_eltypes(m) +## VectorInterface.scalartype(m::Array{<:Array{ITensor}}) = LinearAlgebra.promote_leaf_eltypes(m) """ eltype(m::MPS) @@ -485,7 +486,7 @@ end Get the site index (or indices) of MPO `A` that is unique to `A` (not shared with MPS/MPO `B`). """ -function ITensorQuantumOperatorDefinitions.siteinds( +function siteinds( f::Union{typeof(uniqueinds),typeof(uniqueind)}, A::AbstractMPS, B::AbstractMPS, @@ -505,7 +506,7 @@ end Get the site indices of MPO `A` that are unique to `A` (not shared with MPS/MPO `B`), as a `Vector{<:Index}`. """ -function ITensorQuantumOperatorDefinitions.siteinds( +function siteinds( f::Union{typeof(uniqueinds),typeof(uniqueind)}, A::AbstractMPS, B::AbstractMPS; kwargs... ) return [siteinds(f, A, B, j; kwargs...) for j in eachindex(A)] @@ -517,7 +518,7 @@ end Get the site index (or indices) of the `j`th MPO tensor of `A` that is shared with MPS/MPO `B`. """ -function ITensorQuantumOperatorDefinitions.siteinds( +function siteinds( f::Union{typeof(commoninds),typeof(commonind)}, A::AbstractMPS, B::AbstractMPS, @@ -533,7 +534,7 @@ end Get a vector of the site index (or indices) of MPO `A` that is shared with MPS/MPO `B`. """ -function ITensorQuantumOperatorDefinitions.siteinds( +function siteinds( f::Union{typeof(commoninds),typeof(commonind)}, A::AbstractMPS, B::AbstractMPS; kwargs... ) return [siteinds(f, A, B, j) for j in eachindex(A)] @@ -639,9 +640,7 @@ Return the first site Index found on the MPS or MPO You can choose different filters, like prime level and tags, with the `kwargs`. """ -function ITensorQuantumOperatorDefinitions.siteind( - ::typeof(first), M::AbstractMPS, j::Integer; kwargs... -) +function siteind(::typeof(first), M::AbstractMPS, j::Integer; kwargs...) N = length(M) (N == 1) && return firstind(M[1]; kwargs...) if j == 1 @@ -663,7 +662,7 @@ at the site `j` as an IndexSet. Optionally filter prime tags and prime levels with keyword arguments like `plev` and `tags`. """ -function ITensorQuantumOperatorDefinitions.siteinds(M::AbstractMPS, j::Integer; kwargs...) +function siteinds(M::AbstractMPS, j::Integer; kwargs...) N = length(M) (N == 1) && return inds(M[1]; kwargs...) if j == 1 @@ -676,27 +675,19 @@ function ITensorQuantumOperatorDefinitions.siteinds(M::AbstractMPS, j::Integer; return si end -function ITensorQuantumOperatorDefinitions.siteinds( - ::typeof(all), ψ::AbstractMPS, n::Integer; kwargs... -) +function siteinds(::typeof(all), ψ::AbstractMPS, n::Integer; kwargs...) return siteinds(ψ, n; kwargs...) end -function ITensorQuantumOperatorDefinitions.siteinds( - ::typeof(first), ψ::AbstractMPS; kwargs... -) +function siteinds(::typeof(first), ψ::AbstractMPS; kwargs...) return [siteind(first, ψ, j; kwargs...) for j in 1:length(ψ)] end -function ITensorQuantumOperatorDefinitions.siteinds( - ::typeof(only), ψ::AbstractMPS; kwargs... -) +function siteinds(::typeof(only), ψ::AbstractMPS; kwargs...) return [siteind(only, ψ, j; kwargs...) for j in 1:length(ψ)] end -function ITensorQuantumOperatorDefinitions.siteinds( - ::typeof(all), ψ::AbstractMPS; kwargs... -) +function siteinds(::typeof(all), ψ::AbstractMPS; kwargs...) return [siteinds(ψ, j; kwargs...) for j in 1:length(ψ)] end @@ -738,10 +729,11 @@ for (fname, fname!) in [ (:(ITensors.noprime), :(ITensors.noprime!)), (:(ITensors.swapprime), :(ITensors.swapprime!)), (:(ITensors.replaceprime), :(ITensors.replaceprime!)), - (:(ITensors.addtags), :(ITensors.addtags!)), - (:(ITensors.removetags), :(ITensors.removetags!)), - (:(ITensors.replacetags), :(ITensors.replacetags!)), - (:(ITensors.settags), :(ITensors.settags!)), + # TODO: Delete these, use `settag` instead. + # (:(ITensors.addtags), :(ITensors.addtags!)), + # (:(ITensors.removetags), :(ITensors.removetags!)), + # (:(ITensors.replacetags), :(ITensors.replacetags!)), + # (:(ITensors.settags), :(ITensors.settags!)), ] @eval begin """ @@ -882,10 +874,11 @@ for (fname, fname!) in [ (:(ITensors.prime), :(ITensors.prime!)), (:(ITensors.setprime), :(ITensors.setprime!)), (:(ITensors.noprime), :(ITensors.noprime!)), - (:(ITensors.addtags), :(ITensors.addtags!)), - (:(ITensors.removetags), :(ITensors.removetags!)), - (:(ITensors.replacetags), :(ITensors.replacetags!)), - (:(ITensors.settags), :(ITensors.settags!)), + # TODO: Delete these, use `settag` instead. + # (:(ITensors.addtags), :(ITensors.addtags!)), + # (:(ITensors.removetags), :(ITensors.removetags!)), + # (:(ITensors.replacetags), :(ITensors.replacetags!)), + # (:(ITensors.settags), :(ITensors.settags!)), ] @eval begin """ @@ -1332,7 +1325,7 @@ of the orthogonality center to avoid numerical overflow in the case of diverging See also [`normalize!`](@ref), [`norm`](@ref), [`lognorm`](@ref). """ -function normalize(M::AbstractMPS; (lognorm!)=[]) +function LinearAlgebra.normalize(M::AbstractMPS; (lognorm!)=[]) return normalize!(deepcopy_ortho_center(M); (lognorm!)=lognorm!) end @@ -1604,7 +1597,7 @@ truncation. - `cutoff::Real`: singular value truncation cutoff - `maxdim::Int`: maximum MPS/MPO bond dimension """ -function sum(ψ⃗::Vector{T}; kwargs...) where {T<:AbstractMPS} +function Base.sum(ψ⃗::Vector{T}; kwargs...) where {T<:AbstractMPS} iszero(length(ψ⃗)) && return T() isone(length(ψ⃗)) && return copy(only(ψ⃗)) return +(ψ⃗...; kwargs...) @@ -2432,7 +2425,7 @@ function Base.show(io::IO, M::AbstractMPS) println(io, "#undef") else A = M[i] - if order(A) != 0 + if ndims(A) != 0 println(io, "[$i] $(inds(A))") else println(io, "[$i] ITensor()") diff --git a/src/abstractprojmpo/abstractprojmpo.jl b/src/abstractprojmpo/abstractprojmpo.jl index 393bbb4..1502d43 100644 --- a/src/abstractprojmpo/abstractprojmpo.jl +++ b/src/abstractprojmpo/abstractprojmpo.jl @@ -1,5 +1,4 @@ -# TODO: Deprecate, use `ndims` instead. -using ITensors: order +using ITensors: noprime abstract type AbstractProjMPO end @@ -72,7 +71,7 @@ shorthand for `product(P,v)`. """ function product(P::AbstractProjMPO, v::ITensor)::ITensor Pv = contract(P, v) - if order(Pv) != order(v) + if ndims(Pv) != ndims(v) error( string( "The order of the ProjMPO-ITensor product P*v is not equal to the order of the ITensor v, ", diff --git a/src/deprecated.jl b/src/deprecated.jl index 35c775b..94875b2 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -8,7 +8,6 @@ @deprecate error_mpoprod(args...; kwargs...) error_contract(args...; kwargs...) @deprecate error_mul(args...; kwargs...) error_contract(args...; kwargs...) @deprecate multMPO(args...; kwargs...) contract(args...; kwargs...) -@deprecate sum(A::AbstractMPS, B::AbstractMPS; kwargs...) add(A, B; kwargs...) @deprecate multmpo(args...; kwargs...) contract(args...; kwargs...) @deprecate set_leftlim!(args...; kwargs...) ITensors.setleftlim!(args...; kwargs...) @deprecate set_rightlim!(args...; kwargs...) ITensors.setrightlim!(args...; kwargs...) diff --git a/src/imports.jl b/src/imports.jl index 6779a7b..820a85d 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -1,18 +1,10 @@ # Primarily used to import names into the `ITensorMPS` # module from submodules or from `ITensors` so they can # be reexported. -using ITensorQuantumOperatorDefinitions: - @OpName_str, - @SiteType_str, - @StateName_str, - @TagType_str, - @ValName_str, - OpName, - SiteType, - StateName, - TagType, - ValName, - ops +# TODO: Delete these and expect users to load +# `QuantumOperatorDefinitions`? +using QuantumOperatorDefinitions: + @OpName_str, @SiteType_str, @StateName_str, OpName, SiteType, StateName using QuantumOperatorAlgebra: Trotter import Base: @@ -118,55 +110,55 @@ import Base.Broadcast: ## permutedims, ## permutedims! -import ..ITensors: - Apply, - apply, - argument, - Broadcasted, - @Algorithm_str, - checkflux, - convert_leaf_eltype, - commontags, - @debug_check, - dag, - data, - DefaultArrayStyle, - DiskVector, - flux, - hascommoninds, - hasqns, - hassameinds, - inner, - isfermionic, - maxdim, - mindim, - noprime, - noprime!, - norm, - normalize, - outer, - OneITensor, - permute, - prime, - prime!, - product, - QNIndex, - replaceinds, - replaceprime, - replacetags, - setprime, - sim, - site, - splitblocks, - store, - Style, - sum, - swapprime, - symmetrystyle, - terms, - @timeit_debug, - truncate!, - which_op +## import ..ITensors: +## Apply, +## apply, +## argument, +## Broadcasted, +## @Algorithm_str, +## checkflux, +## convert_leaf_eltype, +## commontags, +## @debug_check, +## dag, +## data, +## DefaultArrayStyle, +## DiskVector, +## flux, +## hascommoninds, +## hasqns, +## hassameinds, +## inner, +## isfermionic, +## maxdim, +## mindim, +## noprime, +## noprime!, +## norm, +## normalize, +## outer, +## OneITensor, +## permute, +## prime, +## prime!, +## product, +## QNIndex, +## replaceinds, +## replaceprime, +## replacetags, +## setprime, +## sim, +## site, +## splitblocks, +## store, +## Style, +## sum, +## swapprime, +## symmetrystyle, +## terms, +## @timeit_debug, +## truncate!, +## which_op import QuantumOperatorAlgebra: params diff --git a/src/lib/ITensorMPSNamedDimsArraysExt/src/ITensorMPSNamedDimsArraysExt.jl b/src/lib/ITensorMPSNamedDimsArraysExt/src/ITensorMPSNamedDimsArraysExt.jl deleted file mode 100644 index 97229dd..0000000 --- a/src/lib/ITensorMPSNamedDimsArraysExt/src/ITensorMPSNamedDimsArraysExt.jl +++ /dev/null @@ -1,3 +0,0 @@ -module ITensorMPSNamedDimsArraysExt -include("to_nameddimsarray.jl") -end diff --git a/src/lib/ITensorMPSNamedDimsArraysExt/src/to_nameddimsarray.jl b/src/lib/ITensorMPSNamedDimsArraysExt/src/to_nameddimsarray.jl deleted file mode 100644 index f6596cb..0000000 --- a/src/lib/ITensorMPSNamedDimsArraysExt/src/to_nameddimsarray.jl +++ /dev/null @@ -1,6 +0,0 @@ -using ..ITensorMPS: AbstractMPS -using ITensors.ITensorsNamedDimsArraysExt: ITensorsNamedDimsArraysExt, to_nameddimsarray - -function ITensorsNamedDimsArraysExt.to_nameddimsarray(x::AbstractMPS) - return to_nameddimsarray.(x) -end diff --git a/src/mpo.jl b/src/mpo.jl index fe476f7..d0fb2a7 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -1,10 +1,12 @@ using Adapt: adapt using LinearAlgebra: dot using Random: Random, AbstractRNG -using ITensors: @ts_str, Algorithm, dag, outer +using ITensors: dag +# TODO: Remove from ITensors? +using ITensors: @ts_str, Algorithm +# TODO: Add this back. +# using ITensors: outer using QuantumOperatorAlgebra: OpSum -using ITensorQuantumOperatorDefinitions: - ITensorQuantumOperatorDefinitions, siteind, siteinds """ MPO @@ -44,7 +46,7 @@ function MPO(::Type{ElT}, sites::Vector{<:Index}) where {ElT<:Number} return MPO(v) end space_ii = all(hasqns, sites) ? [QN() => 1] : 1 - l = [Index(space_ii, "Link,l=$ii") for ii in 1:(N - 1)] + l = [Index(space_ii; tags=Dict("l" => "$ii")) for ii in 1:(N - 1)] for ii in eachindex(sites) s = sites[ii] if ii == 1 @@ -148,6 +150,8 @@ function outer_mps_mps_deprecation_warning() return "Calling `outer(ψ::MPS, ϕ::MPS)` for MPS `ψ` and `ϕ` with shared indices is deprecated. Currently, we automatically prime `ψ` to make sure the site indices don't clash, but that will no longer be the case in ITensors v0.4. To upgrade your code, call `outer(ψ', ϕ)`. Although the new interface seems less convenient, it will allow `outer` to accept more general outer products going forward, such as outer products where some indices are shared (a batched outer product) or outer products of MPS between site indices that aren't just related by a single prime level." end +function outer end + function deprecate_make_inds_unmatch(::typeof(outer), ψ::MPS, ϕ::MPS; kw...) if hassameinds(siteinds, ψ, ϕ) ITensors.warn_once(outer_mps_mps_deprecation_warning(), :outer_mps_mps) @@ -243,8 +247,7 @@ end Get the first site Index of the MPO found, by default with prime level 0. """ -ITensorQuantumOperatorDefinitions.siteind(M::MPO, j::Int; kwargs...) = - siteind(first, M, j; plev=0, kwargs...) +siteind(M::MPO, j::Int; kwargs...) = siteind(first, M, j; plev=0, kwargs...) # TODO: make this return the site indices that would have # been used to create the MPO? I.e.: @@ -254,9 +257,9 @@ ITensorQuantumOperatorDefinitions.siteind(M::MPO, j::Int; kwargs...) = Get a Vector of IndexSets of all the site indices of M. """ -ITensorQuantumOperatorDefinitions.siteinds(M::MPO; kwargs...) = siteinds(all, M; kwargs...) +siteinds(M::MPO; kwargs...) = siteinds(all, M; kwargs...) -function ITensorQuantumOperatorDefinitions.siteinds(Mψ::Tuple{MPO,MPS}, n::Int; kwargs...) +function siteinds(Mψ::Tuple{MPO,MPS}, n::Int; kwargs...) return siteinds(uniqueinds, Mψ[1], Mψ[2], n; kwargs...) end @@ -267,7 +270,7 @@ function nsites(Mψ::Tuple{MPO,MPS}) return N end -function ITensorQuantumOperatorDefinitions.siteinds(Mψ::Tuple{MPO,MPS}; kwargs...) +function siteinds(Mψ::Tuple{MPO,MPS}; kwargs...) return [siteinds(Mψ, n; kwargs...) for n in 1:nsites(Mψ)] end @@ -452,9 +455,10 @@ Same as [`dot`](@ref). """ inner(y::MPS, A::MPO, x::MPS; kwargs...) = dot(y, A, x; kwargs...) -function inner(y::MPS, Ax::Apply{Tuple{MPO,MPS}}) - return inner(y', Ax.args[1], Ax.args[2]) -end +## TODO: Add this back? +## function inner(y::MPS, Ax::Apply{Tuple{MPO,MPS}}) +## return inner(y', Ax.args[1], Ax.args[2]) +## end """ loginner(y::MPS, A::MPO, x::MPS) diff --git a/src/mps.jl b/src/mps.jl index a56852c..8e03d47 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -1,7 +1,6 @@ using Adapt: adapt using ITensors: hasqns -using ITensorQuantumOperatorDefinitions: - ITensorQuantumOperatorDefinitions, siteind, siteinds, state +using QuantumOperatorDefinitions: QuantumOperatorDefinitions, state using Random: Random, AbstractRNG ## TODO: Add this back. @@ -418,11 +417,11 @@ function MPS(eltype::Type{<:Number}, sites::Vector{<:Index}, states_) M = MPS(N) if N == 1 - M[1] = state(sites[1], states_[1]) + M[1] = state(states_[1], sites[1]) return convert_leaf_eltype(eltype, M) end - states = [state(sites[j], states_[j]) for j in 1:N] + states = [state(states_[j], sites[j]) for j in 1:N] if hasqns(states[1]) lflux = QN() @@ -431,21 +430,21 @@ function MPS(eltype::Type{<:Number}, sites::Vector{<:Index}, states_) end links = Vector{QNIndex}(undef, N - 1) for j in (N - 1):-1:1 - links[j] = dag(Index(lflux => 1; tags="Link,l=$j")) + links[j] = dag(Index(lflux => 1; tags=Dict("l" => "$j"))) lflux -= flux(states[j]) end else - links = [Index(1; tags="Link,l=$n") for n in 1:N] + links = [Index(1; tags=Dict("l" => "$n")) for n in 1:N] end M[1] = ITensor(sites[1], links[1]) - M[1] += states[1] * state(links[1], 1) + M[1] += states[1] * state(1, links[1]) for n in 2:(N - 1) M[n] = ITensor(dag(links[n - 1]), sites[n], links[n]) - M[n] += state(dag(links[n - 1]), 1) * states[n] * state(links[n], 1) + M[n] += state(1, dag(links[n - 1])) * states[n] * state(1, links[n]) end M[N] = ITensor(dag(links[N - 1]), sites[N]) - M[N] += state(dag(links[N - 1]), 1) * states[N] + M[N] += state(1, dag(links[N - 1])) * states[N] return convert_leaf_eltype(eltype, M) end @@ -487,17 +486,14 @@ MPS(sites::Vector{<:Index}, states) = MPS(Float64, sites, states) Get the first site Index of the MPS. Return `nothing` if none is found. """ -ITensorQuantumOperatorDefinitions.siteind(M::MPS, j::Int; kwargs...) = - siteind(first, M, j; kwargs...) +siteind(M::MPS, j::Int; kwargs...) = siteind(first, M, j; kwargs...) """ siteind(::typeof(only), M::MPS, j::Int; kwargs...) Get the only site Index of the MPS. Return `nothing` if none is found. """ -function ITensorQuantumOperatorDefinitions.siteind( - ::typeof(only), M::MPS, j::Int; kwargs... -) +function siteind(::typeof(only), M::MPS, j::Int; kwargs...) is = siteinds(M, j; kwargs...) if isempty(is) return nothing @@ -519,8 +515,7 @@ Get a vector of the only site Index found on each tensor of the MPS. Errors if m Get a vector of the all site Indices found on each tensor of the MPS. Returns a Vector of IndexSets. """ -ITensorQuantumOperatorDefinitions.siteinds(M::MPS; kwargs...) = - siteinds(first, M; kwargs...) +siteinds(M::MPS; kwargs...) = siteinds(first, M; kwargs...) function replace_siteinds!(M::MPS, sites) for j in eachindex(M) diff --git a/src/opsum_to_mpo/opsum_to_mpo.jl b/src/opsum_to_mpo/opsum_to_mpo.jl index 97b7e8e..2d0f7d0 100644 --- a/src/opsum_to_mpo/opsum_to_mpo.jl +++ b/src/opsum_to_mpo/opsum_to_mpo.jl @@ -159,7 +159,7 @@ function svdMPO( end llinks = Vector{Index{Int}}(undef, N + 1) - llinks[1] = Index(2, "Link,l=0") + llinks[1] = Index(2; tags=Dict("l" => "0")) H = MPO(sites) @@ -171,7 +171,7 @@ function svdMPO( VR = Vs[n] tdim = isempty(rightmaps[n]) ? 0 : size(VR, 2) - llinks[n + 1] = Index(2 + tdim, "Link,l=$n") + llinks[n + 1] = Index(2 + tdim; tags=Dict("l" => "$n")) ll = llinks[n] rl = llinks[n + 1] diff --git a/src/opsum_to_mpo/opsum_to_mpo_generic.jl b/src/opsum_to_mpo/opsum_to_mpo_generic.jl index 79a5f5a..97c94ed 100644 --- a/src/opsum_to_mpo/opsum_to_mpo_generic.jl +++ b/src/opsum_to_mpo/opsum_to_mpo_generic.jl @@ -11,7 +11,7 @@ using QuantumOperatorAlgebra: sites, terms, which_op -using ITensorQuantumOperatorDefinitions: has_fermion_string, op +using QuantumOperatorDefinitions: has_fermion_string, op """ parity_sign(P) diff --git a/src/siteinds.jl b/src/siteinds.jl new file mode 100644 index 0000000..eb2c58f --- /dev/null +++ b/src/siteinds.jl @@ -0,0 +1,8 @@ +using ITensors: Index, settag +using QuantumOperatorDefinitions: SiteType + +# TODO: Move to `ITensorSites.jl`? +# TODO: Add support for symmetry sectors. +function siteinds(t::String, n::Int) + return [settag(Index(length(SiteType(t))), "sitetype", t) for _ in 1:n] +end diff --git a/src/solvers/ITensorsExtensions.jl b/src/solvers/ITensorsExtensions.jl index 90ed8b1..2782697 100644 --- a/src/solvers/ITensorsExtensions.jl +++ b/src/solvers/ITensorsExtensions.jl @@ -1,9 +1,13 @@ module ITensorsExtensions -using ITensors: ITensor, array, inds, itensor + +using ITensors: ITensor, inds +using NamedDimsArrays: unname + function to_vec(x::ITensor) function to_itensor(x_vec) - return itensor(x_vec, inds(x)) + return ITensor(x_vec, inds(x)) end - return vec(array(x)), to_itensor + return vec(unname(x)), to_itensor end + end diff --git a/src/solvers/alternating_update.jl b/src/solvers/alternating_update.jl index 75ecc4b..330dfca 100644 --- a/src/solvers/alternating_update.jl +++ b/src/solvers/alternating_update.jl @@ -1,4 +1,5 @@ -using ITensors: ITensors, permute +using ITensors: ITensors +using VectorInterface: scalartype function _extend_sweeps_param(param, nsweeps) if param isa Number @@ -39,7 +40,7 @@ function alternating_update( normalize=default_normalize(), maxdim=default_maxdim(), mindim=default_mindim(), - cutoff=default_cutoff(ITensors.scalartype(init)), + cutoff=default_cutoff(scalartype(init)), noise=default_noise(), ) reduced_operator = ITensorMPS.reduced_operator(operator) diff --git a/src/solvers/dmrg.jl b/src/solvers/dmrg.jl index 90fc9df..d6be0ff 100644 --- a/src/solvers/dmrg.jl +++ b/src/solvers/dmrg.jl @@ -1,5 +1,6 @@ using ITensors: ITensors using KrylovKit: eigsolve +using VectorInterface: scalartype function eigsolve_updater( operator, @@ -7,7 +8,7 @@ function eigsolve_updater( internal_kwargs, which_eigval=:SR, ishermitian=true, - tol=10^2 * eps(real(ITensors.scalartype(init))), + tol=10^2 * eps(real(scalartype(init))), krylovdim=3, maxiter=1, verbosity=0, diff --git a/src/solvers/dmrg_x.jl b/src/solvers/dmrg_x.jl index 48201e6..92c5f33 100644 --- a/src/solvers/dmrg_x.jl +++ b/src/solvers/dmrg_x.jl @@ -1,12 +1,13 @@ -using ITensors: array, contract, dag, uniqueind, onehot +using ITensors: contract, dag, uniqueind, onehot using LinearAlgebra: eigen +using NamedDimsArrays: unname function eigen_updater(operator, state; internal_kwargs) contracted_operator = contract(operator, ITensor(true)) d, u = eigen(contracted_operator; ishermitian=true) u_ind = uniqueind(u, contracted_operator) u′_ind = uniqueind(d, u) - max_overlap, max_index = findmax(abs, array(state * dag(u))) + max_overlap, max_index = findmax(abs, unname(state * dag(u))) u_max = u * dag(onehot(eltype(u), u_ind => max_index)) d_max = d[u′_ind => max_index, u_ind => max_index] return u_max, (; eigval=d_max) diff --git a/src/solvers/expand.jl b/src/solvers/expand.jl index 7291390..15fbc6d 100644 --- a/src/solvers/expand.jl +++ b/src/solvers/expand.jl @@ -12,10 +12,10 @@ using ITensors: directsum, hasqns, prime, - scalartype, uniqueinds using LinearAlgebra: normalize, svd, tr using TypeParameterAccessors: unwrap_array_type +using VectorInterface: scalartype # Possible improvements: # - Allow a maxdim argument to be passed to `expand`. diff --git a/src/solvers/sweep_update.jl b/src/solvers/sweep_update.jl index 46ffa3e..3a3754c 100644 --- a/src/solvers/sweep_update.jl +++ b/src/solvers/sweep_update.jl @@ -1,6 +1,7 @@ using ITensors: ITensors, uniqueinds using LinearAlgebra: norm, normalize!, svd using Printf: @printf +using VectorInterface: scalartype function sweep_update( order::TDVPOrder, @@ -76,7 +77,7 @@ function sub_sweep_update( outputlevel=default_outputlevel(), maxdim=default_maxdim(), mindim=default_mindim(), - cutoff=default_cutoff(ITensors.scalartype(state)), + cutoff=default_cutoff(scalartype(state)), noise=default_noise(), ) reduced_operator = copy(reduced_operator) diff --git a/src/solvers/timedependentsum.jl b/src/solvers/timedependentsum.jl index 9acf431..ed8dd12 100644 --- a/src/solvers/timedependentsum.jl +++ b/src/solvers/timedependentsum.jl @@ -1,4 +1,5 @@ -using ITensors: ITensor, inds, permute +using ITensors: ITensor, inds +using NamedDimsArrays: aligndims using QuantumOperatorAlgebra: QuantumOperatorAlgebra # Represents a time-dependent sum of terms: @@ -66,4 +67,4 @@ function scaledsum_apply(expr, x) end end (expr::ScaledSum)(x) = scaledsum_apply(expr, x) -(expr::ScaledSum)(x::ITensor) = permute(scaledsum_apply(expr, x), inds(x)) +(expr::ScaledSum)(x::ITensor) = aligndims(scaledsum_apply(expr, x), inds(x)) From 549b05d7177a4649d7181128afae416c81855e4e Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 28 Jan 2025 16:02:44 -0500 Subject: [PATCH 12/33] Simplifications --- Project.toml | 4 +++ examples/dmrg/Project.toml | 2 +- src/abstractmps.jl | 44 +++++++++++++------------- src/abstractprojmpo/abstractprojmpo.jl | 13 ++++---- src/abstractprojmpo/diskprojmpo.jl | 19 +++++------ src/dmrg.jl | 14 ++++---- src/exports.jl | 2 +- src/mpo.jl | 18 ++++++----- src/mps.jl | 8 ++--- src/opsum_to_mpo/opsum_to_mpo.jl | 6 ++-- src/opsum_to_mpo/opsum_to_mpo_qn.jl | 6 ++-- src/solvers/contract.jl | 2 +- src/solvers/expand.jl | 2 +- src/solvers/reducedconstantterm.jl | 4 +-- 14 files changed, 76 insertions(+), 68 deletions(-) diff --git a/Project.toml b/Project.toml index 7185dc6..2a125e3 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,8 @@ version = "0.4.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" +DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77" +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" IsApprox = "28f27b66-4bd8-47e7-9110-e2746eb8bed7" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" @@ -38,6 +40,8 @@ ITensorMPSZygoteRulesExt = ["ChainRulesCore", "ZygoteRules"] Adapt = "4.1.0" ChainRulesCore = "1.10" Compat = "4.16.0" +DiagonalArrays = "0.2.3" +FillArrays = "1.13.0" HDF5 = "0.14, 0.15, 0.16, 0.17" ITensors = "0.8.0" IsApprox = "2.0.0" diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index c88c435..68b25a2 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -1,7 +1,7 @@ [deps] +ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -QuantumOperatorDefinitions = "826dd319-6fd5-459a-a990-3a4f214664bf" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" diff --git a/src/abstractmps.jl b/src/abstractmps.jl index c18ab7e..baa3b2b 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -619,7 +619,7 @@ findsites(ψ::AbstractMPS, s::Index) = findsites(ψ, IndexSet(s)) Return the first site of the MPS or MPO that has the site index `i`. """ -findfirstsiteind(ψ::AbstractMPS, s::Index) = findfirst(hasind(s), ψ) +findfirstsiteind(ψ::AbstractMPS, s::Index) = findfirst(ψᵢ -> s ∈ inds(ψᵢ), ψ) # TODO: depracate in favor of findsite. """ @@ -629,7 +629,7 @@ findfirstsiteind(ψ::AbstractMPS, s::Index) = findfirst(hasind(s), ψ) Return the first site of the MPS or MPO that has the site indices `is`. """ -findfirstsiteinds(ψ::AbstractMPS, s) = findfirst(hasinds(s), ψ) +findfirstsiteinds(ψ::AbstractMPS, s) = findfirst(ψᵢ -> s ⊆ inds(ψᵢ), ψ) """ siteind(::typeof(first), M::Union{MPS,MPO}, j::Integer; kwargs...) @@ -724,16 +724,16 @@ end for (fname, fname!) in [ (:(ITensors.dag), :(dag!)), - (:(ITensors.prime), :(ITensors.prime!)), - (:(ITensors.setprime), :(ITensors.setprime!)), - (:(ITensors.noprime), :(ITensors.noprime!)), - (:(ITensors.swapprime), :(ITensors.swapprime!)), - (:(ITensors.replaceprime), :(ITensors.replaceprime!)), + (:(ITensors.prime), :(prime!)), + (:(ITensors.setprime), :(setprime!)), + (:(ITensors.noprime), :(noprime!)), + (:(ITensors.swapprime), :(swapprime!)), + (:(ITensors.replaceprime), :(replaceprime!)), # TODO: Delete these, use `settag` instead. - # (:(ITensors.addtags), :(ITensors.addtags!)), - # (:(ITensors.removetags), :(ITensors.removetags!)), - # (:(ITensors.replacetags), :(ITensors.replacetags!)), - # (:(ITensors.settags), :(ITensors.settags!)), + # (:(ITensors.addtags), :(addtags!)), + # (:(ITensors.removetags), :(removetags!)), + # (:(ITensors.replacetags), :(replacetags!)), + # (:(ITensors.settags), :(settags!)), ] @eval begin """ @@ -871,14 +871,14 @@ end for (fname, fname!) in [ (:(ITensors.sim), :(sim!)), - (:(ITensors.prime), :(ITensors.prime!)), - (:(ITensors.setprime), :(ITensors.setprime!)), - (:(ITensors.noprime), :(ITensors.noprime!)), + (:(ITensors.prime), :(prime!)), + (:(ITensors.setprime), :(setprime!)), + (:(ITensors.noprime), :(noprime!)), # TODO: Delete these, use `settag` instead. - # (:(ITensors.addtags), :(ITensors.addtags!)), - # (:(ITensors.removetags), :(ITensors.removetags!)), - # (:(ITensors.replacetags), :(ITensors.replacetags!)), - # (:(ITensors.settags), :(ITensors.settags!)), + # (:(ITensors.addtags), :(addtags!)), + # (:(ITensors.removetags), :(removetags!)), + # (:(ITensors.replacetags), :(replacetags!)), + # (:(ITensors.settags), :(settags!)), ] @eval begin """ @@ -1955,9 +1955,9 @@ function (::Type{MPST})( ) where {MPST<:AbstractMPS} N = length(sites) for s in sites - @assert hasinds(A, s) + @assert s ⊆ inds(A) end - @assert isnothing(leftinds) || hasinds(A, leftinds) + @assert isnothing(leftinds) || leftinds ⊆ inds(A) @assert 1 ≤ orthocenter ≤ N @@ -1986,7 +1986,7 @@ function (::Type{MPST})( end function (::Type{MPST})(A::AbstractArray, sites; kwargs...) where {MPST<:AbstractMPS} - return MPST(itensor(A, sites...), sites; kwargs...) + return MPST(ITensor(A, sites...), sites; kwargs...) end """ @@ -2220,7 +2220,7 @@ function ITensors.op(::OpName"CX", ::SiteType"S=1/2", s1::Index, s2::Index) 0 1 0 0 0 0 0 1 0 0 1 0] - return itensor(mat, s2', s1', s2, s1) + return ITensor(mat, s2', s1', s2, s1) end os = [("CX", 1, 3), ("σz", 3)] diff --git a/src/abstractprojmpo/abstractprojmpo.jl b/src/abstractprojmpo/abstractprojmpo.jl index 1502d43..260f8c3 100644 --- a/src/abstractprojmpo/abstractprojmpo.jl +++ b/src/abstractprojmpo/abstractprojmpo.jl @@ -1,3 +1,4 @@ +using FillArrays: Ones using ITensors: noprime abstract type AbstractProjMPO end @@ -25,18 +26,18 @@ the length of the MPO used to construct it """ Base.length(P::AbstractProjMPO) = length(P.H) -function lproj(P::AbstractProjMPO)::Union{ITensor,OneITensor} - (P.lpos <= 0) && return OneITensor() +function lproj(P::AbstractProjMPO) + (P.lpos <= 0) && return ITensor(Ones{Bool}()) return P.LR[P.lpos] end -function rproj(P::AbstractProjMPO)::Union{ITensor,OneITensor} - (P.rpos >= length(P) + 1) && return OneITensor() +function rproj(P::AbstractProjMPO) + (P.rpos >= length(P) + 1) && return ITensor(Ones{Bool}()) return P.LR[P.rpos] end -function ITensors.contract(P::AbstractProjMPO, v::ITensor)::ITensor - itensor_map = Union{ITensor,OneITensor}[lproj(P)] +function ITensors.contract(P::AbstractProjMPO, v::ITensor) + itensor_map = [lproj(P)] append!(itensor_map, P.H[site_range(P)]) push!(itensor_map, rproj(P)) diff --git a/src/abstractprojmpo/diskprojmpo.jl b/src/abstractprojmpo/diskprojmpo.jl index 7630a15..f70d3b2 100644 --- a/src/abstractprojmpo/diskprojmpo.jl +++ b/src/abstractprojmpo/diskprojmpo.jl @@ -1,4 +1,5 @@ -using ITensors: OneITensor +using FillArrays: Ones +using ITensors: ITensor using SerializedElementArrays: SerializedElementVector """ @@ -29,9 +30,9 @@ mutable struct DiskProjMPO <: AbstractProjMPO nsite::Int H::MPO LR::SerializedElementVector{ITensor} - Lcache::Union{ITensor,OneITensor} + Lcache::ITensor lposcache::Union{Int,Nothing} - Rcache::Union{ITensor,OneITensor} + Rcache::ITensor rposcache::Union{Int,Nothing} end @@ -61,9 +62,9 @@ function DiskProjMPO(H::MPO) 2, H, disk(Vector{ITensor}(undef, length(H))), - OneITensor, + ITensor, nothing, - OneITensor, + ITensor, nothing, ) end @@ -86,8 +87,8 @@ disk(pm::DiskProjMPO; kwargs...) = pm # Special overload of lproj which uses the cached # version of the left projected MPO, and if the # cache doesn't exist it loads it from disk. -function lproj(P::DiskProjMPO)::Union{ITensor,OneITensor} - (P.lpos <= 0) && return OneITensor() +function lproj(P::DiskProjMPO) + (P.lpos <= 0) && return ITensor(Ones{Bool}()) if (P.lpos ≠ P.lposcache) || (P.lpos == 1) # Need to update the cache P.Lcache = P.LR[P.lpos] @@ -99,8 +100,8 @@ end # Special overload of rproj which uses the cached # version of the right projected MPO, and if the # cache doesn't exist it loads it from disk. -function rproj(P::DiskProjMPO)::Union{ITensor,OneITensor} - (P.rpos >= length(P) + 1) && return OneITensor() +function rproj(P::DiskProjMPO) + (P.rpos >= length(P) + 1) && return ITensor(Ones{Bool}()) if (P.rpos ≠ P.rposcache) || (P.rpos == length(P)) # Need to update the cache P.Rcache = P.LR[P.rpos] diff --git a/src/dmrg.jl b/src/dmrg.jl index e881e49..ad03190 100644 --- a/src/dmrg.jl +++ b/src/dmrg.jl @@ -1,6 +1,7 @@ using Adapt: adapt using ITensors: ITensors, plev using KrylovKit: eigsolve +using NamedDimsArrays: NamedDimsArrays, aligndims using Printf: @printf using TupleTools: TupleTools using VectorInterface: scalartype @@ -8,7 +9,7 @@ using VectorInterface: scalartype ## TODO: Add this back? ## using NDTensors: timer -function permute( +function NamedDimsArrays.aligndims( M::AbstractMPS, ::Tuple{typeof(linkind),typeof(siteinds),typeof(linkind)} )::typeof(M) M̃ = typeof(M)(length(M)) @@ -16,8 +17,7 @@ function permute( lₙ₋₁ = linkind(M, n - 1) lₙ = linkind(M, n) s⃗ₙ = TupleTools.sort(Tuple(siteinds(M, n)); by=plev) - # TODO: Use `NamedDimsArrays.aligndims` instead. - M̃[n] = ITensors.permute(M[n], filter(!isnothing, (lₙ₋₁, s⃗ₙ..., lₙ))) + M̃[n] = aligndims(M[n], filter(!isnothing, (lₙ₋₁, s⃗ₙ..., lₙ))) end set_ortho_lims!(M̃, ortho_lims(M)) return M̃ @@ -28,7 +28,7 @@ function dmrg(H::MPO, psi0::MPS, sweeps::Sweeps; kwargs...) check_hascommoninds(siteinds, H, psi0') # Permute the indices to have a better memory layout # and minimize permutations - H = permute(H, (linkind, siteinds, linkind)) + H = aligndims(H, (linkind, siteinds, linkind)) PH = ProjMPO(H) return dmrg(PH, psi0, sweeps; kwargs...) end @@ -38,7 +38,7 @@ function dmrg(Hs::Vector{MPO}, psi0::MPS, sweeps::Sweeps; kwargs...) check_hascommoninds(siteinds, H, psi0) check_hascommoninds(siteinds, H, psi0') end - Hs .= permute.(Hs, Ref((linkind, siteinds, linkind))) + Hs .= aligndims.(Hs, Ref((linkind, siteinds, linkind))) PHS = ProjMPOSum(Hs) return dmrg(PHS, psi0, sweeps; kwargs...) end @@ -49,8 +49,8 @@ function dmrg(H::MPO, Ms::Vector{MPS}, psi0::MPS, sweeps::Sweeps; weight=true, k for M in Ms check_hascommoninds(siteinds, M, psi0) end - H = permute(H, (linkind, siteinds, linkind)) - Ms .= permute.(Ms, Ref((linkind, siteinds, linkind))) + H = aligndims(H, (linkind, siteinds, linkind)) + Ms .= aligndims.(Ms, Ref((linkind, siteinds, linkind))) if weight <= 0 error( "weight parameter should be > 0.0 in call to excited-state dmrg (value passed was weight=$weight)", diff --git a/src/exports.jl b/src/exports.jl index a620fd4..417a90b 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -7,7 +7,7 @@ export @StateName_str, @TagType_str, @ValName_str, - Apply, + ## Apply, Op, OpName, Prod, diff --git a/src/mpo.jl b/src/mpo.jl index d0fb2a7..cf8768b 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -1,12 +1,13 @@ using Adapt: adapt -using LinearAlgebra: dot -using Random: Random, AbstractRNG +using DiagonalArrays: δ using ITensors: dag +using LinearAlgebra: dot # TODO: Remove from ITensors? -using ITensors: @ts_str, Algorithm -# TODO: Add this back. +using ITensors: Algorithm +# TODO: Add this back? # using ITensors: outer using QuantumOperatorAlgebra: OpSum +using Random: Random, AbstractRNG """ MPO @@ -547,7 +548,7 @@ function logdot(M1::MPO, M2::MPO; make_inds_match::Bool=false, kwargs...) return _log_or_not_dot(M1, M2, true; make_inds_match=make_inds_match) end -function LinearAlgebra.tr(M::MPO; plev::Pair{Int,Int}=0 => 1, tags::Pair=ts"" => ts"") +function LinearAlgebra.tr(M::MPO; plev::Pair{Int,Int}=0 => 1, tags::Pair="" => "") N = length(M) # # TODO: choose whether to contract or trace @@ -615,9 +616,10 @@ end (A::MPO)(ψ::MPS; kwargs...) = apply(A, ψ; kwargs...) -function Apply(A::MPO, ψ::MPS; kwargs...) - return ITensors.LazyApply.Applied(apply, (A, ψ), NamedTuple(kwargs)) -end +## TODO: Decide what to do about this. +## function Apply(A::MPO, ψ::MPS; kwargs...) +## return ITensors.LazyApply.Applied(apply, (A, ψ), NamedTuple(kwargs)) +## end function ITensors.contract(A::MPO, ψ::MPS; alg=nothing, method=alg, kwargs...) # TODO: Delete `method` since it is deprecated. diff --git a/src/mps.jl b/src/mps.jl index 8e03d47..3e0e648 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -99,7 +99,7 @@ function randomU(rng::AbstractRNG, eltype::Type{<:Number}, s1::Index, s2::Index) mdim = dim(s1) * dim(s2) RM = randn(rng, eltype, mdim, mdim) Q, _ = NDTensors.qr_positive(RM) - G = itensor(Q, dag(s1), dag(s2), s1', s2') + G = ITensor(Q, dag(s1), dag(s2), s1', s2') else M = random_itensor(rng, eltype, QN(), s1', s2', dag(s1), dag(s2)) U, S, V = svd(M, (s1', s2')) @@ -185,7 +185,7 @@ function randomCircuitMPS( chi = min(linkdims[N - 1], d) l[N - 1] = Index(chi, "Link,l=$(N-1)") O = NDTensors.random_unitary(rng, eltype, chi, d) - M[N] = itensor(O, l[N - 1], sites[N]) + M[N] = ITensor(O, l[N - 1], sites[N]) for j in (N - 1):-1:2 chi *= dim(sites[j]) @@ -193,13 +193,13 @@ function randomCircuitMPS( l[j - 1] = Index(chi, "Link,l=$(j-1)") O = NDTensors.random_unitary(rng, eltype, chi, dim(sites[j]) * dim(l[j])) T = reshape(O, (chi, dim(sites[j]), dim(l[j]))) - M[j] = itensor(T, l[j - 1], sites[j], l[j]) + M[j] = ITensor(T, l[j - 1], sites[j], l[j]) end O = NDTensors.random_unitary(rng, eltype, 1, dim(sites[1]) * dim(l[1])) l0 = Index(1, "Link,l=0") T = reshape(O, (1, dim(sites[1]), dim(l[1]))) - M[1] = itensor(T, l0, sites[1], l[1]) + M[1] = ITensor(T, l0, sites[1], l[1]) M[1] *= onehot(eltype, l0 => 1) M.llim = 0 diff --git a/src/opsum_to_mpo/opsum_to_mpo.jl b/src/opsum_to_mpo/opsum_to_mpo.jl index 2d0f7d0..f76f00a 100644 --- a/src/opsum_to_mpo/opsum_to_mpo.jl +++ b/src/opsum_to_mpo/opsum_to_mpo.jl @@ -1,5 +1,5 @@ ## using NDTensors: using_auto_fermion -using ITensors: dim, inds, itensor +using ITensors: dim, inds replace_nothing(::Nothing, y) = y replace_nothing(x, y) = x @@ -206,7 +206,7 @@ function svdMPO( end end - T = itensor(M, ll, rl) + T = ITensor(M, ll, rl) H[n] += T * computeSiteProd(sites, argument(t)) end @@ -218,7 +218,7 @@ function svdMPO( idM = zeros(ValType, dim(ll), dim(rl)) idM[1, 1] = 1.0 idM[end, end] = 1.0 - T = itensor(idM, ll, rl) + T = ITensor(idM, ll, rl) H[n] += T * computeSiteProd(sites, Prod([Op("Id", n)])) end diff --git a/src/opsum_to_mpo/opsum_to_mpo_qn.jl b/src/opsum_to_mpo/opsum_to_mpo_qn.jl index a1f31a2..1f9b98c 100644 --- a/src/opsum_to_mpo/opsum_to_mpo_qn.jl +++ b/src/opsum_to_mpo/opsum_to_mpo_qn.jl @@ -225,7 +225,7 @@ ## T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl)) ## T[b] .= M ## -## H[n] += (itensor(T) * Op) +## H[n] += (ITensor(T) * Op) ## end ## end ## @@ -234,13 +234,13 @@ ## b = Block(1, 1) ## T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl)) ## T[b] = 1 -## H[n] += (itensor(T) * Id) +## H[n] += (ITensor(T) * Id) ## ## # Put in starting identity operator ## b = Block(nblocks(ll), nblocks(rl)) ## T = ITensors.NDTensors.BlockSparseTensor(ValType, [b], (dag(ll), rl)) ## T[b] = 1 -## H[n] += (itensor(T) * Id) +## H[n] += (ITensor(T) * Id) ## end # for n in 1:N ## ## L = ITensor(llinks[1]) diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl index 1919fdd..bda56e6 100644 --- a/src/solvers/contract.jl +++ b/src/solvers/contract.jl @@ -1,4 +1,4 @@ -using ITensors: ITensors, Index, ITensor, @Algorithm_str, commoninds, contract, hasind, sim +using ITensors: ITensors, ITensor, Index, @Algorithm_str, contract function contract_operator_state_updater(operator, init; internal_kwargs) # TODO: Use `contract(operator)`. diff --git a/src/solvers/expand.jl b/src/solvers/expand.jl index 15fbc6d..a85e113 100644 --- a/src/solvers/expand.jl +++ b/src/solvers/expand.jl @@ -1,11 +1,11 @@ using Adapt: adapt +using DiagonalArrays: δ using ITensors: ITensors, Algorithm, Index, ITensor, @Algorithm_str, - δ, commonind, dag, denseblocks, diff --git a/src/solvers/reducedconstantterm.jl b/src/solvers/reducedconstantterm.jl index 52c882b..ff254ad 100644 --- a/src/solvers/reducedconstantterm.jl +++ b/src/solvers/reducedconstantterm.jl @@ -105,7 +105,7 @@ function ITensorMPS.makeR!(reduced_state::ReducedConstantTerm, basis::MPS, posit end function ITensors.contract(reduced_state::ReducedConstantTerm, v::ITensor) - reduced_state_tensors = Union{ITensor,OneITensor}[lproj(reduced_state)] + reduced_state_tensors = [lproj(reduced_state)] append!( reduced_state_tensors, [prime(t, "Link") for t in reduced_state.state[site_range(reduced_state)]], @@ -129,7 +129,7 @@ end # Contract the reduced constant term down to a since ITensor. function ITensors.contract(reduced_state::ReducedConstantTerm) - reduced_state_tensors = Union{ITensor,OneITensor}[lproj(reduced_state)] + reduced_state_tensors = [lproj(reduced_state)] append!( reduced_state_tensors, [dag(prime(t, "Link")) for t in reduced_state.state[site_range(reduced_state)]], From 79d4b7e9dec01cfbe0ff16635889cf077bdd648f Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 28 Jan 2025 17:49:48 -0500 Subject: [PATCH 13/33] Use BackendSelection --- Project.toml | 6 +- examples/dmrg/Project.toml | 1 - ext/ITensorMPSPackageCompilerExt/compile.jl | 2 +- src/abstractmps.jl | 4 +- src/imports.jl | 74 ++------------------- src/mpo.jl | 3 +- src/solvers/contract.jl | 4 +- src/solvers/expand.jl | 3 +- src/solvers/tdvp.jl | 2 +- 9 files changed, 17 insertions(+), 82 deletions(-) diff --git a/Project.toml b/Project.toml index 2a125e3..e54e632 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.4.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +BackendSelection = "680c2d7c-f67a-4cc9-ae9c-da132b1447a5" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" @@ -18,6 +19,7 @@ QuantumOperatorAlgebra = "c7a6d0f7-daa6-4368-ba67-89ed64127c3b" QuantumOperatorDefinitions = "826dd319-6fd5-459a-a990-3a4f214664bf" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SerializedElementArrays = "d3ce8812-9567-47e9-a7b5-65a6d70a3065" +TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" @@ -38,6 +40,7 @@ ITensorMPSZygoteRulesExt = ["ChainRulesCore", "ZygoteRules"] [compat] Adapt = "4.1.0" +BackendSelection = "0.1.0" ChainRulesCore = "1.10" Compat = "4.16.0" DiagonalArrays = "0.2.3" @@ -47,7 +50,7 @@ ITensors = "0.8.0" IsApprox = "2.0.0" KrylovKit = "0.8.1, 0.9" LinearAlgebra = "1.10" -NamedDimsArrays = "0.3.12" +NamedDimsArrays = "0.4" Observers = "0.2" PackageCompiler = "1, 2" Printf = "1.10" @@ -55,6 +58,7 @@ QuantumOperatorAlgebra = "0.1.0" QuantumOperatorDefinitions = "0.1.2" Random = "1.10" SerializedElementArrays = "0.1.0" +TensorAlgebra = "0.1.7" Test = "1.10" TupleTools = "1.5.0" TypeParameterAccessors = "0.2.1" diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index 68b25a2..e6b9177 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -1,5 +1,4 @@ [deps] -ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/ext/ITensorMPSPackageCompilerExt/compile.jl b/ext/ITensorMPSPackageCompilerExt/compile.jl index c14d0e9..58f3452 100644 --- a/ext/ITensorMPSPackageCompilerExt/compile.jl +++ b/ext/ITensorMPSPackageCompilerExt/compile.jl @@ -1,4 +1,4 @@ -using NDTensors: @Algorithm_str +using BackendSelection: @Algorithm_str using ITensors: ITensors using PackageCompiler: PackageCompiler diff --git a/src/abstractmps.jl b/src/abstractmps.jl index baa3b2b..7463bd6 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -1,6 +1,4 @@ -# TODO: Define in `BackendSelection.jl`. -using ITensors: @Algorithm_str - +using BackendSelection: @Algorithm_str, Algorithm using IsApprox: Approx, IsApprox using ITensors: ITensors, diff --git a/src/imports.jl b/src/imports.jl index 820a85d..1083896 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -7,6 +7,7 @@ using QuantumOperatorDefinitions: @OpName_str, @SiteType_str, @StateName_str, OpName, SiteType, StateName using QuantumOperatorAlgebra: Trotter +# TODO: Delete this, overload explicitly. import Base: # types Array, @@ -77,6 +78,7 @@ import Base: # macros @propagate_inbounds +# TODO: Delete this, overload explicitly. import Base.Broadcast: # types AbstractArrayStyle, @@ -90,76 +92,8 @@ import Base.Broadcast: broadcastable, instantiate -## import ..ITensors.NDTensors: -## Algorithm, -## @Algorithm_str, -## EmptyNumber, -## _Tuple, -## _NTuple, -## blas_get_num_threads, -## datatype, -## dense, -## diagind, -## disable_auto_fermion, -## double_precision, -## eachblock, -## eachdiagblock, -## enable_auto_fermion, -## fill!!, -## randn!!, -## permutedims, -## permutedims! - -## import ..ITensors: -## Apply, -## apply, -## argument, -## Broadcasted, -## @Algorithm_str, -## checkflux, -## convert_leaf_eltype, -## commontags, -## @debug_check, -## dag, -## data, -## DefaultArrayStyle, -## DiskVector, -## flux, -## hascommoninds, -## hasqns, -## hassameinds, -## inner, -## isfermionic, -## maxdim, -## mindim, -## noprime, -## noprime!, -## norm, -## normalize, -## outer, -## OneITensor, -## permute, -## prime, -## prime!, -## product, -## QNIndex, -## replaceinds, -## replaceprime, -## replacetags, -## setprime, -## sim, -## site, -## splitblocks, -## store, -## Style, -## sum, -## swapprime, -## symmetrystyle, -## terms, -## @timeit_debug, -## truncate!, -## which_op - +# TODO: Delete this, overload explicitly. import QuantumOperatorAlgebra: params +# TODO: Delete this, overload explicitly. import SerializedElementArrays: disk diff --git a/src/mpo.jl b/src/mpo.jl index cf8768b..775fdc9 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -1,9 +1,8 @@ using Adapt: adapt +using BackendSelection: @Algorithm_str, Algorithm using DiagonalArrays: δ using ITensors: dag using LinearAlgebra: dot -# TODO: Remove from ITensors? -using ITensors: Algorithm # TODO: Add this back? # using ITensors: outer using QuantumOperatorAlgebra: OpSum diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl index bda56e6..664f307 100644 --- a/src/solvers/contract.jl +++ b/src/solvers/contract.jl @@ -1,4 +1,6 @@ -using ITensors: ITensors, ITensor, Index, @Algorithm_str, contract +using BackendSelection: @Algorithm_str +using ITensors: ITensors, ITensor, Index +using TensorAlgebra: contract function contract_operator_state_updater(operator, init; internal_kwargs) # TODO: Use `contract(operator)`. diff --git a/src/solvers/expand.jl b/src/solvers/expand.jl index a85e113..68029d8 100644 --- a/src/solvers/expand.jl +++ b/src/solvers/expand.jl @@ -1,11 +1,10 @@ using Adapt: adapt +using BackendSelection: @Algorithm_str, Algorithm using DiagonalArrays: δ using ITensors: ITensors, - Algorithm, Index, ITensor, - @Algorithm_str, commonind, dag, denseblocks, diff --git a/src/solvers/tdvp.jl b/src/solvers/tdvp.jl index e8ed9d6..08b2ce4 100644 --- a/src/solvers/tdvp.jl +++ b/src/solvers/tdvp.jl @@ -1,4 +1,4 @@ -using ITensors: Algorithm, @Algorithm_str +using BackendSelection: @Algorithm_str, Algorithm using KrylovKit: exponentiate function exponentiate_updater(operator, init; internal_kwargs, kwargs...) From 58a3a0b80ff80bb89c5df63e515c0c71d72529f0 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 28 Jan 2025 18:06:24 -0500 Subject: [PATCH 14/33] Delete deprecations --- src/ITensorMPS.jl | 1 - src/deprecated.jl | 105 ---------------------------------------------- 2 files changed, 106 deletions(-) delete mode 100644 src/deprecated.jl diff --git a/src/ITensorMPS.jl b/src/ITensorMPS.jl index 288085e..e044803 100644 --- a/src/ITensorMPS.jl +++ b/src/ITensorMPS.jl @@ -21,7 +21,6 @@ include("opsum_to_mpo/qnmatelem.jl") include("opsum_to_mpo/opsum_to_mpo_generic.jl") include("opsum_to_mpo/opsum_to_mpo.jl") include("opsum_to_mpo/opsum_to_mpo_qn.jl") -include("deprecated.jl") include("defaults.jl") include("update_observer.jl") include("lattices/lattices.jl") diff --git a/src/deprecated.jl b/src/deprecated.jl deleted file mode 100644 index 94875b2..0000000 --- a/src/deprecated.jl +++ /dev/null @@ -1,105 +0,0 @@ -# mps/abstractmps.jl -@deprecate orthoCenter(args...; kwargs...) orthocenter(args...; kwargs...) -@deprecate store(m::AbstractMPS) data(m) false -@deprecate replacesites!(args...; kwargs...) ITensors.replace_siteinds!(args...; kwargs...) -@deprecate applyMPO(args...; kwargs...) contract(args...; kwargs...) -@deprecate applympo(args...; kwargs...) contract(args...; kwargs...) -@deprecate errorMPOprod(args...; kwargs...) error_contract(args...; kwargs...) -@deprecate error_mpoprod(args...; kwargs...) error_contract(args...; kwargs...) -@deprecate error_mul(args...; kwargs...) error_contract(args...; kwargs...) -@deprecate multMPO(args...; kwargs...) contract(args...; kwargs...) -@deprecate multmpo(args...; kwargs...) contract(args...; kwargs...) -@deprecate set_leftlim!(args...; kwargs...) ITensors.setleftlim!(args...; kwargs...) -@deprecate set_rightlim!(args...; kwargs...) ITensors.setrightlim!(args...; kwargs...) -@deprecate tensors(args...; kwargs...) ITensors.data(args...; kwargs...) -@deprecate primelinks!(args...; kwargs...) ITensors.prime_linkinds!(args...; kwargs...) -@deprecate simlinks!(args...; kwargs...) ITensors.sim_linkinds!(args...; kwargs...) -@deprecate mul(A::AbstractMPS, B::AbstractMPS; kwargs...) contract(A, B; kwargs...) - -# mps/mpo.jl -@deprecate toMPO(args...; kwargs...) MPO(args...; kwargs...) -@deprecate MPO(A::MPS; kwargs...) outer(A', A; kwargs...) -@deprecate randomMPO(args...; kwargs...) random_mpo(args...; kwargs...) - -# mps/mps.jl -@deprecate randomMPS(args...; kwargs...) random_mps(args...; kwargs...) - -# Deprecated syntax for specifying link dimensions. -@deprecate randomMPS(elt::Type{<:Number}, sites::Vector{<:Index}, state, linkdims::Integer) random_mps( - elt, sites, state; linkdims -) -@deprecate randomMPS(elt::Type{<:Number}, sites::Vector{<:Index}, linkdims::Integer) random_mps( - elt, sites; linkdims -) -@deprecate randomMPS(sites::Vector{<:Index}, state, linkdims::Integer) random_mps( - sites, state; linkdims -) -@deprecate randomMPS(sites::Vector{<:Index}, linkdims::Integer) random_mps(sites; linkdims) - -# Pass throughs of old name to new name: - -unique_siteind(A::AbstractMPS, B::AbstractMPS, j::Integer) = siteinds(uniqueind, A, B, j) - -unique_siteinds(A::AbstractMPS, B::AbstractMPS, j::Integer) = siteinds(uniqueinds, A, B, j) - -unique_siteinds(A::AbstractMPS, B::AbstractMPS) = siteinds(uniqueind, A, B) - -common_siteind(A::AbstractMPS, B::AbstractMPS, j::Integer) = siteinds(commonind, A, B, j) - -common_siteinds(A::AbstractMPS, B::AbstractMPS, j::Integer) = siteinds(commoninds, A, B, j) - -common_siteinds(A::AbstractMPS, B::AbstractMPS) = siteinds(commoninds, A, B) - -firstsiteind(M::AbstractMPS, j::Integer; kwargs...) = siteind(first, M, j; kwargs...) - -map_linkinds!(f::Function, M::AbstractMPS) = map!(f, linkinds, M) - -map_linkinds(f::Function, M::AbstractMPS) = map(f, linkinds, M) - -function map_common_siteinds!(f::Function, M1::AbstractMPS, M2::AbstractMPS) - return map!(f, siteinds, commoninds, M1, M2) -end - -function map_common_siteinds(f::Function, M1::AbstractMPS, M2::AbstractMPS) - return map(f, siteinds, commoninds, M1, M2) -end - -function map_unique_siteinds!(f::Function, M1::AbstractMPS, M2::AbstractMPS) - return map!(f, siteinds, uniqueinds, M1, M2) -end - -function map_unique_siteinds(f::Function, M1::AbstractMPS, M2::AbstractMPS) - return map(f, siteinds, uniqueinds, M1, M2) -end - -for fname in - (:sim, :prime, :setprime, :noprime, :addtags, :removetags, :replacetags, :settags) - @eval begin - function $(Symbol(fname, :_linkinds))(M::AbstractMPS, args...; kwargs...) - return map(i -> $fname(i, args...; kwargs...), linkinds, M) - end - function $(Symbol(fname, :_linkinds!))(M::AbstractMPS, args...; kwargs...) - return map!(i -> $fname(i, args...; kwargs...), linkinds, M) - end - function $(Symbol(fname, :_common_siteinds))( - M1::AbstractMPS, M2::AbstractMPS, args...; kwargs... - ) - return map(i -> $fname(i, args...; kwargs...), siteinds, commoninds, M1, M2) - end - function $(Symbol(fname, :_common_siteinds!))( - M1::AbstractMPS, M2::AbstractMPS, args...; kwargs... - ) - return map!(i -> $fname(i, args...; kwargs...), siteinds, commoninds, M1, M2) - end - function $(Symbol(fname, :_unique_siteinds))( - M1::AbstractMPS, M2::AbstractMPS, args...; kwargs... - ) - return map(i -> $fname(i, args...; kwargs...), siteinds, uniqueinds, M1, M2) - end - function $(Symbol(fname, :_unique_siteinds!))( - M1::AbstractMPS, M2::AbstractMPS, args...; kwargs... - ) - return map!(i -> $fname(i, args...; kwargs...), siteinds, uniqueinds, M1, M2) - end - end -end From f738f5837a022289022eb62ac3c744e0c4ecd968 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 31 Jan 2025 12:42:37 -0500 Subject: [PATCH 15/33] Improve siteinds, simplify example --- examples/dmrg/1d_heisenberg_conserve_spin.jl | 56 +++++++++++--------- examples/dmrg/Project.toml | 1 + src/siteinds.jl | 7 ++- 3 files changed, 34 insertions(+), 30 deletions(-) diff --git a/examples/dmrg/1d_heisenberg_conserve_spin.jl b/examples/dmrg/1d_heisenberg_conserve_spin.jl index f8ec1fb..9421c5e 100644 --- a/examples/dmrg/1d_heisenberg_conserve_spin.jl +++ b/examples/dmrg/1d_heisenberg_conserve_spin.jl @@ -1,39 +1,43 @@ using ITensors, ITensorMPS -using LinearAlgebra using Printf using Random -using Strided - -Random.seed!(1234) -BLAS.set_num_threads(1) -Strided.set_num_threads(1) -ITensors.enable_threaded_blocksparse() -#ITensors.disable_threaded_blocksparse() - -let - N = 100 - - sites = siteinds("S=1", N; conserve_qns=true) +using SymmetrySectors +function heisenberg(N) os = OpSum() for j in 1:(N - 1) os += 0.5, "S+", j, "S-", j + 1 os += 0.5, "S-", j, "S+", j + 1 os += "Sz", j, "Sz", j + 1 end - H = MPO(os, sites) + return os +end + +Random.seed!(1234) - state = [isodd(n) ? "Up" : "Dn" for n in 1:N] - psi0 = random_mps(sites, state; linkdims=10) +N = 10 - # Plan to do 5 DMRG sweeps: - nsweeps = 5 - # Set maximum MPS bond dimensions for each sweep - maxdim = [10, 20, 100, 100, 200] - # Set maximum truncation error allowed when adapting bond dimensions - cutoff = [1E-11] +s = siteinds("S=1/2", N; gradings=("Sz",)) - # Run the DMRG algorithm, returning energy and optimized MPS - energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff) - @printf("Final energy = %.12f\n", energy) -end +os = heisenberg(N) + +# Input operator terms which define a Hamiltonian +# Convert these terms to an MPO tensor network +H = MPO(os, s) + +# Create an initial random matrix product state +# psi0 = random_mps(s; linkdims=10) +psi0 = MPS(s, j -> isodd(j) ? "↑" : "↓") + +# Plan to do 5 DMRG sweeps: +nsweeps = 5 +# Set maximum MPS bond dimensions for each sweep +maxdim = [10] +# Set maximum truncation error allowed when adapting bond dimensions +cutoff = [1E-11] + +# Run the DMRG algorithm, returning energy and optimized MPS +energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff) +@show inner(psi', H, psi) +@show inner(psi, psi) +@printf("Final energy = %.12f\n", energy) diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index e6b9177..39065e0 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -4,3 +4,4 @@ ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" +SymmetrySectors = "f8a8ad64-adbc-4fce-92f7-ffe2bb36a86e" diff --git a/src/siteinds.jl b/src/siteinds.jl index eb2c58f..1e5e864 100644 --- a/src/siteinds.jl +++ b/src/siteinds.jl @@ -1,8 +1,7 @@ -using ITensors: Index, settag +using ITensors: ITensors, Index, settag using QuantumOperatorDefinitions: SiteType # TODO: Move to `ITensorSites.jl`? -# TODO: Add support for symmetry sectors. -function siteinds(t::String, n::Int) - return [settag(Index(length(SiteType(t))), "sitetype", t) for _ in 1:n] +function siteinds(t::String, n::Int; kwargs...) + return [settag(Index(SiteType(t; kwargs...)), "n", string(j)) for j in 1:n] end From d896610f10202d550026bf39742d3c786545dcb3 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sun, 2 Feb 2025 08:43:52 -0500 Subject: [PATCH 16/33] Simplified siteinds --- src/siteinds.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/siteinds.jl b/src/siteinds.jl index 1e5e864..9877979 100644 --- a/src/siteinds.jl +++ b/src/siteinds.jl @@ -1,7 +1,7 @@ using ITensors: ITensors, Index, settag -using QuantumOperatorDefinitions: SiteType +using QuantumOperatorDefinitions: sites # TODO: Move to `ITensorSites.jl`? function siteinds(t::String, n::Int; kwargs...) - return [settag(Index(SiteType(t; kwargs...)), "n", string(j)) for j in 1:n] + return sites(Index, t, n; kwargs...) end From 1846b5832a82a30d5c1654304465df95acc50292 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sun, 2 Feb 2025 22:03:38 -0500 Subject: [PATCH 17/33] Replace `onehot` with new `SparseArraysBase.oneelement` --- Project.toml | 2 ++ src/mps.jl | 3 ++- src/solvers/dmrg_x.jl | 6 ++++-- test/base/test_autompo.jl | 43 ++++++++++++++++++++------------------- test/base/test_mps.jl | 5 +++-- 5 files changed, 33 insertions(+), 26 deletions(-) diff --git a/Project.toml b/Project.toml index e54e632..c3883d1 100644 --- a/Project.toml +++ b/Project.toml @@ -19,6 +19,7 @@ QuantumOperatorAlgebra = "c7a6d0f7-daa6-4368-ba67-89ed64127c3b" QuantumOperatorDefinitions = "826dd319-6fd5-459a-a990-3a4f214664bf" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SerializedElementArrays = "d3ce8812-9567-47e9-a7b5-65a6d70a3065" +SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208" TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" @@ -58,6 +59,7 @@ QuantumOperatorAlgebra = "0.1.0" QuantumOperatorDefinitions = "0.1.2" Random = "1.10" SerializedElementArrays = "0.1.0" +SparseArraysBase = "0.2.11" TensorAlgebra = "0.1.7" Test = "1.10" TupleTools = "1.5.0" diff --git a/src/mps.jl b/src/mps.jl index 3e0e648..75ec057 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -2,6 +2,7 @@ using Adapt: adapt using ITensors: hasqns using QuantumOperatorDefinitions: QuantumOperatorDefinitions, state using Random: Random, AbstractRNG +using SparseArraysBase: oneelement ## TODO: Add this back. ## using NDTensors: using_auto_fermion @@ -200,7 +201,7 @@ function randomCircuitMPS( l0 = Index(1, "Link,l=0") T = reshape(O, (1, dim(sites[1]), dim(l[1]))) M[1] = ITensor(T, l0, sites[1], l[1]) - M[1] *= onehot(eltype, l0 => 1) + M[1] *= oneelement(eltype, l0 => 1) M.llim = 0 M.rlim = 2 diff --git a/src/solvers/dmrg_x.jl b/src/solvers/dmrg_x.jl index 92c5f33..6f36b81 100644 --- a/src/solvers/dmrg_x.jl +++ b/src/solvers/dmrg_x.jl @@ -1,6 +1,8 @@ -using ITensors: contract, dag, uniqueind, onehot +using ITensors: dag, uniqueind using LinearAlgebra: eigen using NamedDimsArrays: unname +using SparseArraysBase: oneelement +using TensorAlgebra: contract function eigen_updater(operator, state; internal_kwargs) contracted_operator = contract(operator, ITensor(true)) @@ -8,7 +10,7 @@ function eigen_updater(operator, state; internal_kwargs) u_ind = uniqueind(u, contracted_operator) u′_ind = uniqueind(d, u) max_overlap, max_index = findmax(abs, unname(state * dag(u))) - u_max = u * dag(onehot(eltype(u), u_ind => max_index)) + u_max = u * dag(oneelement(eltype(u), u_ind => max_index)) d_max = d[u′_ind => max_index, u_ind => max_index] return u_max, (; eigval=d_max) end diff --git a/test/base/test_autompo.jl b/test/base/test_autompo.jl index c4ebeec..cacdc39 100644 --- a/test/base/test_autompo.jl +++ b/test/base/test_autompo.jl @@ -1,6 +1,7 @@ @eval module $(gensym()) using ITensorMPS, ITensors, Test, Random, JLD2 using NDTensors: scalartype +using SparseArraysBase: oneelement include(joinpath(@__DIR__, "utils", "util.jl")) @@ -98,26 +99,26 @@ function NNheisenbergMPO(sites, J1::Float64, J2::Float64)::MPO ll = dag(link[n]) rl = link[n + 1] H[n] = ITensor(ll, dag(s), s', rl) - H[n] += onehot(ll => 1) * onehot(rl => 1) * op(sites, "Id", n) - H[n] += onehot(ll => 8) * onehot(rl => 8) * op(sites, "Id", n) - - H[n] += onehot(ll => 2) * onehot(rl => 1) * op(sites, "S-", n) - H[n] += onehot(ll => 5) * onehot(rl => 2) * op(sites, "Id", n) - H[n] += onehot(ll => 8) * onehot(rl => 2) * op(sites, "S+", n) * J1 / 2 - H[n] += onehot(ll => 8) * onehot(rl => 5) * op(sites, "S+", n) * J2 / 2 - - H[n] += onehot(ll => 3) * onehot(rl => 1) * op(sites, "S+", n) - H[n] += onehot(ll => 6) * onehot(rl => 3) * op(sites, "Id", n) - H[n] += onehot(ll => 8) * onehot(rl => 3) * op(sites, "S-", n) * J1 / 2 - H[n] += onehot(ll => 8) * onehot(rl => 6) * op(sites, "S-", n) * J2 / 2 - - H[n] += onehot(ll => 4) * onehot(rl => 1) * op(sites, "Sz", n) - H[n] += onehot(ll => 7) * onehot(rl => 4) * op(sites, "Id", n) - H[n] += onehot(ll => 8) * onehot(rl => 4) * op(sites, "Sz", n) * J1 - H[n] += onehot(ll => 8) * onehot(rl => 7) * op(sites, "Sz", n) * J2 + H[n] += oneelement(ll => 1) * oneelement(rl => 1) * op(sites, "Id", n) + H[n] += oneelement(ll => 8) * oneelement(rl => 8) * op(sites, "Id", n) + + H[n] += oneelement(ll => 2) * oneelement(rl => 1) * op(sites, "S-", n) + H[n] += oneelement(ll => 5) * oneelement(rl => 2) * op(sites, "Id", n) + H[n] += oneelement(ll => 8) * oneelement(rl => 2) * op(sites, "S+", n) * J1 / 2 + H[n] += oneelement(ll => 8) * oneelement(rl => 5) * op(sites, "S+", n) * J2 / 2 + + H[n] += oneelement(ll => 3) * oneelement(rl => 1) * op(sites, "S+", n) + H[n] += oneelement(ll => 6) * oneelement(rl => 3) * op(sites, "Id", n) + H[n] += oneelement(ll => 8) * oneelement(rl => 3) * op(sites, "S-", n) * J1 / 2 + H[n] += oneelement(ll => 8) * oneelement(rl => 6) * op(sites, "S-", n) * J2 / 2 + + H[n] += oneelement(ll => 4) * oneelement(rl => 1) * op(sites, "Sz", n) + H[n] += oneelement(ll => 7) * oneelement(rl => 4) * op(sites, "Id", n) + H[n] += oneelement(ll => 8) * oneelement(rl => 4) * op(sites, "Sz", n) * J1 + H[n] += oneelement(ll => 8) * oneelement(rl => 7) * op(sites, "Sz", n) * J2 end - H[1] *= onehot(link[1] => 8) - H[N] *= onehot(dag(link[N + 1]) => 1) + H[1] *= oneelement(link[1] => 8) + H[N] *= oneelement(dag(link[N + 1]) => 1) return H end @@ -1062,8 +1063,8 @@ end q = flux(op(which_op, sites[j])) links = [Index([n < j ? q => 1 : QN() => 1], "Link,l=$n") for n in 1:N] for n in 1:(N - 1) - M[n] *= onehot(links[n] => 1) - M[n + 1] *= onehot(dag(links[n]) => 1) + M[n] *= oneelement(links[n] => 1) + M[n + 1] *= oneelement(dag(links[n]) => 1) end return M end diff --git a/test/base/test_mps.jl b/test/base/test_mps.jl index 65eed0e..229ce89 100644 --- a/test/base/test_mps.jl +++ b/test/base/test_mps.jl @@ -4,6 +4,7 @@ using ITensorMPS using ITensors using LinearAlgebra: diag using Random +using SparseArraysBase: oneelement using Test Random.seed!(1234) @@ -1425,8 +1426,8 @@ end CSWAP = [op("CSWAP", s, n, m, k) for n in 1:N, m in 1:N, k in 1:N] CCCNOT = [op("CCCNOT", s, n, m, k, l) for n in 1:N, m in 1:N, k in 1:N, l in 1:N] - v0 = [onehot(s[n] => "0") for n in 1:N] - v1 = [onehot(s[n] => "1") for n in 1:N] + v0 = [oneelement(s[n] => "0") for n in 1:N] + v1 = [oneelement(s[n] => "1") for n in 1:N] # Single qubit @test product(I[1], v0[1]) ≈ v0[1] From df4a99f9ab639361e13d858685b7f24a944b3437 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 11 Feb 2025 21:46:07 -0500 Subject: [PATCH 18/33] Updates for dependent packages --- Project.toml | 2 ++ src/abstractmps.jl | 4 ++-- src/abstractprojmpo/abstractprojmpo.jl | 1 + src/abstractprojmpo/projmps.jl | 1 + src/mpo.jl | 2 +- src/mps.jl | 1 + src/siteinds.jl | 6 ++++-- src/solvers/dmrg_x.jl | 3 ++- src/solvers/expand.jl | 12 ++---------- src/solvers/reducedconstantterm.jl | 3 ++- 10 files changed, 18 insertions(+), 17 deletions(-) diff --git a/Project.toml b/Project.toml index c3883d1..61fe906 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ BackendSelection = "680c2d7c-f67a-4cc9-ae9c-da132b1447a5" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" +GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" IsApprox = "28f27b66-4bd8-47e7-9110-e2746eb8bed7" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" @@ -46,6 +47,7 @@ ChainRulesCore = "1.10" Compat = "4.16.0" DiagonalArrays = "0.2.3" FillArrays = "1.13.0" +GradedUnitRanges = "0.1.5" HDF5 = "0.14, 0.15, 0.16, 0.17" ITensors = "0.8.0" IsApprox = "2.0.0" diff --git a/src/abstractmps.jl b/src/abstractmps.jl index 7463bd6..196a316 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -1,4 +1,5 @@ using BackendSelection: @Algorithm_str, Algorithm +using GradedUnitRanges: GradedUnitRanges, dag using IsApprox: Approx, IsApprox using ITensors: ITensors, @@ -6,7 +7,6 @@ using ITensors: Index, commonind, commoninds, - dag, factorize, hasqns, replaceinds, @@ -721,7 +721,7 @@ function Base.map(f::Function, M::AbstractMPS; set_limits::Bool=true) end for (fname, fname!) in [ - (:(ITensors.dag), :(dag!)), + (:(GradedUnitRanges.dag), :(dag!)), (:(ITensors.prime), :(prime!)), (:(ITensors.setprime), :(setprime!)), (:(ITensors.noprime), :(noprime!)), diff --git a/src/abstractprojmpo/abstractprojmpo.jl b/src/abstractprojmpo/abstractprojmpo.jl index 260f8c3..38d9532 100644 --- a/src/abstractprojmpo/abstractprojmpo.jl +++ b/src/abstractprojmpo/abstractprojmpo.jl @@ -1,4 +1,5 @@ using FillArrays: Ones +using GradedUnitRanges: dag using ITensors: noprime abstract type AbstractProjMPO end diff --git a/src/abstractprojmpo/projmps.jl b/src/abstractprojmpo/projmps.jl index 0ad60a5..c29db74 100644 --- a/src/abstractprojmpo/projmps.jl +++ b/src/abstractprojmpo/projmps.jl @@ -1,3 +1,4 @@ +using GradedUnitRanges: dag mutable struct ProjMPS lpos::Int diff --git a/src/mpo.jl b/src/mpo.jl index 775fdc9..5959275 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -1,7 +1,7 @@ using Adapt: adapt using BackendSelection: @Algorithm_str, Algorithm using DiagonalArrays: δ -using ITensors: dag +using GradedUnitRanges: dag using LinearAlgebra: dot # TODO: Add this back? # using ITensors: outer diff --git a/src/mps.jl b/src/mps.jl index 75ec057..c798342 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -1,4 +1,5 @@ using Adapt: adapt +using GradedUnitRanges: dag using ITensors: hasqns using QuantumOperatorDefinitions: QuantumOperatorDefinitions, state using Random: Random, AbstractRNG diff --git a/src/siteinds.jl b/src/siteinds.jl index 9877979..819946b 100644 --- a/src/siteinds.jl +++ b/src/siteinds.jl @@ -1,7 +1,9 @@ using ITensors: ITensors, Index, settag -using QuantumOperatorDefinitions: sites +using QuantumOperatorDefinitions: QuantumOperatorDefinitions # TODO: Move to `ITensorSites.jl`? function siteinds(t::String, n::Int; kwargs...) - return sites(Index, t, n; kwargs...) + # TODO: QuantumOperatorAlgebra.jl defines its own `QuantumOperatorAlgebra.sites` + # with a different meaning, decide if we should keep both. + return QuantumOperatorDefinitions.sites(Index, t, n; kwargs...) end diff --git a/src/solvers/dmrg_x.jl b/src/solvers/dmrg_x.jl index 6f36b81..0ed1f49 100644 --- a/src/solvers/dmrg_x.jl +++ b/src/solvers/dmrg_x.jl @@ -1,4 +1,5 @@ -using ITensors: dag, uniqueind +using GradedUnitRanges: dag +using ITensors: uniqueind using LinearAlgebra: eigen using NamedDimsArrays: unname using SparseArraysBase: oneelement diff --git a/src/solvers/expand.jl b/src/solvers/expand.jl index 68029d8..1fca947 100644 --- a/src/solvers/expand.jl +++ b/src/solvers/expand.jl @@ -1,17 +1,9 @@ using Adapt: adapt using BackendSelection: @Algorithm_str, Algorithm using DiagonalArrays: δ +using GradedUnitRanges: dag using ITensors: - ITensors, - Index, - ITensor, - commonind, - dag, - denseblocks, - directsum, - hasqns, - prime, - uniqueinds + ITensors, Index, ITensor, commonind, denseblocks, directsum, hasqns, prime, uniqueinds using LinearAlgebra: normalize, svd, tr using TypeParameterAccessors: unwrap_array_type using VectorInterface: scalartype diff --git a/src/solvers/reducedconstantterm.jl b/src/solvers/reducedconstantterm.jl index ff254ad..c3e206f 100644 --- a/src/solvers/reducedconstantterm.jl +++ b/src/solvers/reducedconstantterm.jl @@ -1,4 +1,5 @@ -using ITensors: ITensors, ITensor, dag, dim, prime +using GradedUnitRanges: dag +using ITensors: ITensors, ITensor, dim, prime """ Holds the following data where basis From 9a9399d19f5ea128b5bb24679c31a087663754d5 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Wed, 12 Feb 2025 11:07:10 -0500 Subject: [PATCH 19/33] Use settag --- src/mpo.jl | 6 +++--- src/mps.jl | 32 ++++++++++++++++---------------- src/opsum_to_mpo/opsum_to_mpo.jl | 6 +++--- 3 files changed, 22 insertions(+), 22 deletions(-) diff --git a/src/mpo.jl b/src/mpo.jl index 5959275..5dfbda6 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -46,7 +46,7 @@ function MPO(::Type{ElT}, sites::Vector{<:Index}) where {ElT<:Number} return MPO(v) end space_ii = all(hasqns, sites) ? [QN() => 1] : 1 - l = [Index(space_ii; tags=Dict("l" => "$ii")) for ii in 1:(N - 1)] + l = [settag(Index(space_ii), "l", "$ii") for ii in 1:(N - 1)] for ii in eachindex(sites) s = sites[ii] if ii == 1 @@ -558,10 +558,10 @@ function LinearAlgebra.tr(M::MPO; plev::Pair{Int,Int}=0 => 1, tags::Pair="" => " # # So tracing first is better if d > √χ. # - L = tr(M[1]; plev=plev, tags=tags) + L = tr(M[1]; plev, tags) for j in 2:N L *= M[j] - L = tr(L; plev=plev, tags=tags) + L = tr(L; plev, tags) end return L end diff --git a/src/mps.jl b/src/mps.jl index c798342..bc332e2 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -1,5 +1,5 @@ using Adapt: adapt -using GradedUnitRanges: dag +using GradedUnitRanges: dual using ITensors: hasqns using QuantumOperatorDefinitions: QuantumOperatorDefinitions, state using Random: Random, AbstractRNG @@ -82,9 +82,9 @@ function MPS( if ii == 1 v[ii] = ITensor(T, l[ii], s) elseif ii == N - v[ii] = ITensor(T, dag(l[ii - 1]), s) + v[ii] = ITensor(T, dual(l[ii - 1]), s) else - v[ii] = ITensor(T, dag(l[ii - 1]), s, l[ii]) + v[ii] = ITensor(T, dual(l[ii - 1]), s, l[ii]) end end return MPS(v) @@ -101,9 +101,9 @@ function randomU(rng::AbstractRNG, eltype::Type{<:Number}, s1::Index, s2::Index) mdim = dim(s1) * dim(s2) RM = randn(rng, eltype, mdim, mdim) Q, _ = NDTensors.qr_positive(RM) - G = ITensor(Q, dag(s1), dag(s2), s1', s2') + G = ITensor(Q, dual(s1), dual(s2), s1', s2') else - M = random_itensor(rng, eltype, QN(), s1', s2', dag(s1), dag(s2)) + M = random_itensor(rng, eltype, QN(), s1', s2', dual(s1), dual(s2)) U, S, V = svd(M, (s1', s2')) u = commonind(U, S) v = commonind(S, V) @@ -354,7 +354,7 @@ function MPS(::Type{T}, ivals::Vector{<:Pair{<:Index}}) where {T<:Number} end links = Vector{QNIndex}(undef, N - 1) for j in (N - 1):-1:1 - links[j] = dag(Index(lflux => 1; tags="Link,l=$j")) + links[j] = dual(Index(lflux => 1; tags="Link,l=$j")) lflux -= qn(ivals[j]) end else @@ -365,10 +365,10 @@ function MPS(::Type{T}, ivals::Vector{<:Pair{<:Index}}) where {T<:Number} M[1][ivals[1], links[1] => 1] = one(T) for n in 2:(N - 1) s = ind(ivals[n]) - M[n] = ITensor(T, dag(links[n - 1]), s, links[n]) + M[n] = ITensor(T, dual(links[n - 1]), s, links[n]) M[n][links[n - 1] => 1, ivals[n], links[n] => 1] = one(T) end - M[N] = ITensor(T, dag(links[N - 1]), ind(ivals[N])) + M[N] = ITensor(T, dual(links[N - 1]), ind(ivals[N])) M[N][links[N - 1] => 1, ivals[N]] = one(T) return M @@ -432,21 +432,21 @@ function MPS(eltype::Type{<:Number}, sites::Vector{<:Index}, states_) end links = Vector{QNIndex}(undef, N - 1) for j in (N - 1):-1:1 - links[j] = dag(Index(lflux => 1; tags=Dict("l" => "$j"))) + links[j] = dual(settag(Index(lflux => 1), "l", "$j")) lflux -= flux(states[j]) end else - links = [Index(1; tags=Dict("l" => "$n")) for n in 1:N] + links = [settag(Index(1), "l", "$n") for n in 1:N] end M[1] = ITensor(sites[1], links[1]) M[1] += states[1] * state(1, links[1]) for n in 2:(N - 1) - M[n] = ITensor(dag(links[n - 1]), sites[n], links[n]) - M[n] += state(1, dag(links[n - 1])) * states[n] * state(1, links[n]) + M[n] = ITensor(dual(links[n - 1]), sites[n], links[n]) + M[n] += state(1, dual(links[n - 1])) * states[n] * state(1, links[n]) end - M[N] = ITensor(dag(links[N - 1]), sites[N]) - M[N] += state(1, dag(links[N - 1])) * states[N] + M[N] = ITensor(dual(links[N - 1]), sites[N]) + M[N] += state(1, dual(links[N - 1])) * states[N] return convert_leaf_eltype(eltype, M) end @@ -807,7 +807,7 @@ function correlation_matrix( L = ITensor(1.0) else lind = commonind(psi[start_site], psi[start_site - 1]) - L = delta(dag(lind), lind') + L = delta(dual(lind), lind') end pL = start_site - 1 @@ -901,7 +901,7 @@ function correlation_matrix( Li21 *= oᵢ * dag(psi[pL21])' else sᵢ = siteind(psi, pL21) - Li21 *= prime(dag(si[pL21]), !sᵢ) + Li21 *= prime(dual(si[pL21]), !sᵢ) end Li21 *= psi[pL21] end diff --git a/src/opsum_to_mpo/opsum_to_mpo.jl b/src/opsum_to_mpo/opsum_to_mpo.jl index f76f00a..916d3e7 100644 --- a/src/opsum_to_mpo/opsum_to_mpo.jl +++ b/src/opsum_to_mpo/opsum_to_mpo.jl @@ -1,5 +1,5 @@ ## using NDTensors: using_auto_fermion -using ITensors: dim, inds +using ITensors: dim, inds, settag replace_nothing(::Nothing, y) = y replace_nothing(x, y) = x @@ -159,7 +159,7 @@ function svdMPO( end llinks = Vector{Index{Int}}(undef, N + 1) - llinks[1] = Index(2; tags=Dict("l" => "0")) + llinks[1] = settag(Index(2), "l", "0") H = MPO(sites) @@ -171,7 +171,7 @@ function svdMPO( VR = Vs[n] tdim = isempty(rightmaps[n]) ? 0 : size(VR, 2) - llinks[n + 1] = Index(2 + tdim; tags=Dict("l" => "$n")) + llinks[n + 1] = settag(Index(2 + tdim), "l", "$n") ll = llinks[n] rl = llinks[n + 1] From b278cb646a817bbeae0fe0a70754c1dec54870d5 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Thu, 13 Feb 2025 13:35:06 -0500 Subject: [PATCH 20/33] Add back non-QN random_mps --- examples/dmrg/1d_heisenberg.jl | 4 ++-- src/mps.jl | 25 +++++++++++++++++++------ 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/examples/dmrg/1d_heisenberg.jl b/examples/dmrg/1d_heisenberg.jl index 11938be..e22df2e 100644 --- a/examples/dmrg/1d_heisenberg.jl +++ b/examples/dmrg/1d_heisenberg.jl @@ -28,8 +28,8 @@ os = heisenberg(N) H = MPO(os, sites) # Create an initial random matrix product state -# psi0 = random_mps(sites; linkdims=10) -psi0 = MPS(sites, j -> isodd(j) ? "↑" : "↓") +psi0 = random_mps(sites; linkdims=10) +# psi0 = MPS(sites, j -> isodd(j) ? "↑" : "↓") # Plan to do 5 DMRG sweeps: nsweeps = 5 diff --git a/src/mps.jl b/src/mps.jl index bc332e2..f7047d9 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -1,6 +1,7 @@ using Adapt: adapt using GradedUnitRanges: dual using ITensors: hasqns +using LinearAlgebra: qr using QuantumOperatorDefinitions: QuantumOperatorDefinitions, state using Random: Random, AbstractRNG using SparseArraysBase: oneelement @@ -165,6 +166,18 @@ function randomCircuitMPS( return randomCircuitMPS(Random.default_rng(), eltype, sites, linkdims; kwargs...) end +function random_unitary(rng::AbstractRNG, elt::Type, n::Int, m::Int) + if n < m + return copy(random_unitary(rng, elt, m, n)') + end + F = qr(randn(rng, elt, n, m)) + Q = Matrix(F.Q) + for c in 1:size(Q, 2) + Q[:, c] .*= sign(F.factors[c, c]) + end + return Q +end + function randomCircuitMPS( rng::AbstractRNG, eltype::Type{<:Number}, @@ -185,21 +198,21 @@ function randomCircuitMPS( d = dim(sites[N]) chi = min(linkdims[N - 1], d) - l[N - 1] = Index(chi, "Link,l=$(N-1)") - O = NDTensors.random_unitary(rng, eltype, chi, d) + l[N - 1] = settag(Index(chi), "l", "$(N-1)") + O = random_unitary(rng, eltype, chi, d) M[N] = ITensor(O, l[N - 1], sites[N]) for j in (N - 1):-1:2 chi *= dim(sites[j]) chi = min(linkdims[j - 1], chi) - l[j - 1] = Index(chi, "Link,l=$(j-1)") - O = NDTensors.random_unitary(rng, eltype, chi, dim(sites[j]) * dim(l[j])) + l[j - 1] = settag(Index(chi), "l", "$(j-1)") + O = random_unitary(rng, eltype, chi, dim(sites[j]) * dim(l[j])) T = reshape(O, (chi, dim(sites[j]), dim(l[j]))) M[j] = ITensor(T, l[j - 1], sites[j], l[j]) end - O = NDTensors.random_unitary(rng, eltype, 1, dim(sites[1]) * dim(l[1])) - l0 = Index(1, "Link,l=0") + O = random_unitary(rng, eltype, 1, dim(sites[1]) * dim(l[1])) + l0 = settag(Index(1), "l", "0") T = reshape(O, (1, dim(sites[1]), dim(l[1]))) M[1] = ITensor(T, l0, sites[1], l[1]) M[1] *= oneelement(eltype, l0 => 1) From 775f79708a10166649c528c4d40a2716e10538e3 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 14 Feb 2025 09:08:32 -0500 Subject: [PATCH 21/33] Get more of random_mps working --- Project.toml | 2 +- examples/dmrg/1d_heisenberg.jl | 1 + examples/dmrg/1d_heisenberg_conserve_spin.jl | 9 ++++--- examples/dmrg/Project.toml | 3 +++ src/abstractmps.jl | 8 +++--- src/mpo.jl | 2 +- src/mps.jl | 27 +++++++++++++++----- 7 files changed, 36 insertions(+), 16 deletions(-) diff --git a/Project.toml b/Project.toml index 61fe906..93c23ce 100644 --- a/Project.toml +++ b/Project.toml @@ -65,7 +65,7 @@ SparseArraysBase = "0.2.11" TensorAlgebra = "0.1.7" Test = "1.10" TupleTools = "1.5.0" -TypeParameterAccessors = "0.2.1" +TypeParameterAccessors = "0.3" VectorInterface = "0.5.0" ZygoteRules = "0.2.2" julia = "1.10" diff --git a/examples/dmrg/1d_heisenberg.jl b/examples/dmrg/1d_heisenberg.jl index e22df2e..3bbe7c0 100644 --- a/examples/dmrg/1d_heisenberg.jl +++ b/examples/dmrg/1d_heisenberg.jl @@ -29,6 +29,7 @@ H = MPO(os, sites) # Create an initial random matrix product state psi0 = random_mps(sites; linkdims=10) +# psi0 = random_mps(sites, j -> isodd(j) ? "↑" : "↓"; linkdims=10) # psi0 = MPS(sites, j -> isodd(j) ? "↑" : "↓") # Plan to do 5 DMRG sweeps: diff --git a/examples/dmrg/1d_heisenberg_conserve_spin.jl b/examples/dmrg/1d_heisenberg_conserve_spin.jl index 9421c5e..ccf54d0 100644 --- a/examples/dmrg/1d_heisenberg_conserve_spin.jl +++ b/examples/dmrg/1d_heisenberg_conserve_spin.jl @@ -1,3 +1,4 @@ +using BlockSparseArrays using ITensors, ITensorMPS using Printf using Random @@ -21,14 +22,14 @@ s = siteinds("S=1/2", N; gradings=("Sz",)) os = heisenberg(N) +# Create an initial random matrix product state +psi0 = random_mps(s, j -> isodd(j) ? "↑" : "↓"; linkdims=10) +# psi0 = MPS(s, j -> isodd(j) ? "↑" : "↓") + # Input operator terms which define a Hamiltonian # Convert these terms to an MPO tensor network H = MPO(os, s) -# Create an initial random matrix product state -# psi0 = random_mps(s; linkdims=10) -psi0 = MPS(s, j -> isodd(j) ? "↑" : "↓") - # Plan to do 5 DMRG sweeps: nsweeps = 5 # Set maximum MPS bond dimensions for each sweep diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index 39065e0..a5c923c 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -1,7 +1,10 @@ [deps] +BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" +GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" SymmetrySectors = "f8a8ad64-adbc-4fce-92f7-ffe2bb36a86e" +TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" diff --git a/src/abstractmps.jl b/src/abstractmps.jl index 196a316..b93eb9c 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -818,8 +818,8 @@ function Base.map!( if !isempty(s) s̃ = f.(s) @preserve_ortho (M1, M2) begin - M1[i] = replaceinds(M1[i], s .=> s̃) - M2[i] = replaceinds(M2[i], s .=> s̃) + M1[i] = replaceinds(M1[i], (s .=> s̃)...) + M2[i] = replaceinds(M2[i], (s .=> s̃)...) end end end @@ -834,7 +834,7 @@ function Base.map!( s = siteinds(uniqueinds, M1, M2, i) if !isempty(s) @preserve_ortho M1 begin - M1[i] = replaceinds(M1[i], s .=> f(s)) + M1[i] = replaceinds(M1[i], (s .=> f(s))...) end end end @@ -956,7 +956,7 @@ for (fname, fname!) in [ args...; kwargs..., ) - return map(i -> $fname(i, args...; kwargs...), ffilter, fsubset, M1, M2) + return map(i -> $fname.(i, args...; kwargs...), ffilter, fsubset, M1, M2) end function $(fname!)( diff --git a/src/mpo.jl b/src/mpo.jl index 5dfbda6..87bac8d 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -622,7 +622,7 @@ end function ITensors.contract(A::MPO, ψ::MPS; alg=nothing, method=alg, kwargs...) # TODO: Delete `method` since it is deprecated. - alg = NDTensors.replace_nothing(method, "densitymatrix") + alg = replace_nothing(method, "densitymatrix") # Keyword argument deprecations # TODO: Delete these. diff --git a/src/mps.jl b/src/mps.jl index f7047d9..47a7a1c 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -101,7 +101,7 @@ function randomU(rng::AbstractRNG, eltype::Type{<:Number}, s1::Index, s2::Index) if !hasqns(s1) && !hasqns(s2) mdim = dim(s1) * dim(s2) RM = randn(rng, eltype, mdim, mdim) - Q, _ = NDTensors.qr_positive(RM) + Q, _ = qr_positive(RM) G = ITensor(Q, dual(s1), dual(s2), s1', s2') else M = random_itensor(rng, eltype, QN(), s1', s2', dual(s1), dual(s2)) @@ -140,6 +140,11 @@ function randomizeMPS!( s1 = sites[b] s2 = sites[b + db] G = randomU(rng, eltype, s1, s2) + + @show inds(G) + @show inds(M[b]) + @show inds(M[b + db]) + T = noprime(G * M[b] * M[b + db]) rinds = uniqueinds(M[b], M[b + db]) @@ -166,15 +171,25 @@ function randomCircuitMPS( return randomCircuitMPS(Random.default_rng(), eltype, sites, linkdims; kwargs...) end +# TODO: Move this to MatrixAlgebraKit.jl/TensorAlgebra.jl? +using LinearAlgebra: Diagonal, diag +function nonzero_sign(x) + iszero(x) && return one(x) + return sign(x) +end +function qr_positive(M::AbstractMatrix) + Q′, R = qr(M) + Q = convert(typeof(R), Q′) + signs = nonzero_sign.(diag(R)) + Q = Q * Diagonal(signs) + R = Diagonal(conj.(signs)) * R + return Q, R +end function random_unitary(rng::AbstractRNG, elt::Type, n::Int, m::Int) if n < m return copy(random_unitary(rng, elt, m, n)') end - F = qr(randn(rng, elt, n, m)) - Q = Matrix(F.Q) - for c in 1:size(Q, 2) - Q[:, c] .*= sign(F.factors[c, c]) - end + Q, _ = qr_positive(randn(rng, elt, n, m)) return Q end From 4b9759ac28898eb3a9fda4ff36edd6a546e2a4b4 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 14 Feb 2025 09:56:11 -0500 Subject: [PATCH 22/33] Fix shape issue in randomU --- examples/dmrg/1d_heisenberg.jl | 4 ++-- src/mps.jl | 11 ++++------- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/examples/dmrg/1d_heisenberg.jl b/examples/dmrg/1d_heisenberg.jl index 3bbe7c0..86bd7fd 100644 --- a/examples/dmrg/1d_heisenberg.jl +++ b/examples/dmrg/1d_heisenberg.jl @@ -28,8 +28,8 @@ os = heisenberg(N) H = MPO(os, sites) # Create an initial random matrix product state -psi0 = random_mps(sites; linkdims=10) -# psi0 = random_mps(sites, j -> isodd(j) ? "↑" : "↓"; linkdims=10) +# psi0 = random_mps(sites; linkdims=10) +psi0 = random_mps(sites, j -> isodd(j) ? "↑" : "↓"; linkdims=10) # psi0 = MPS(sites, j -> isodd(j) ? "↑" : "↓") # Plan to do 5 DMRG sweeps: diff --git a/src/mps.jl b/src/mps.jl index 47a7a1c..9d6f316 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -99,10 +99,12 @@ end function randomU(rng::AbstractRNG, eltype::Type{<:Number}, s1::Index, s2::Index) if !hasqns(s1) && !hasqns(s2) - mdim = dim(s1) * dim(s2) + mdim = Int(length(s1)) * Int(length(s2)) RM = randn(rng, eltype, mdim, mdim) Q, _ = qr_positive(RM) - G = ITensor(Q, dual(s1), dual(s2), s1', s2') + inds = (dual(s1), dual(s2), s1', s2') + dims = Int.(length.(inds)) + G = ITensor(reshape(Q, dims), inds) else M = random_itensor(rng, eltype, QN(), s1', s2', dual(s1), dual(s2)) U, S, V = svd(M, (s1', s2')) @@ -140,11 +142,6 @@ function randomizeMPS!( s1 = sites[b] s2 = sites[b + db] G = randomU(rng, eltype, s1, s2) - - @show inds(G) - @show inds(M[b]) - @show inds(M[b + db]) - T = noprime(G * M[b] * M[b + db]) rinds = uniqueinds(M[b], M[b + db]) From dd63d7f31e729a434595c0f2a0cb01bdc57cf96d Mon Sep 17 00:00:00 2001 From: mtfishman Date: Wed, 19 Feb 2025 07:52:30 -0500 Subject: [PATCH 23/33] Fix example --- examples/dmrg/1d_heisenberg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/dmrg/1d_heisenberg.jl b/examples/dmrg/1d_heisenberg.jl index 86bd7fd..3bbe7c0 100644 --- a/examples/dmrg/1d_heisenberg.jl +++ b/examples/dmrg/1d_heisenberg.jl @@ -28,8 +28,8 @@ os = heisenberg(N) H = MPO(os, sites) # Create an initial random matrix product state -# psi0 = random_mps(sites; linkdims=10) -psi0 = random_mps(sites, j -> isodd(j) ? "↑" : "↓"; linkdims=10) +psi0 = random_mps(sites; linkdims=10) +# psi0 = random_mps(sites, j -> isodd(j) ? "↑" : "↓"; linkdims=10) # psi0 = MPS(sites, j -> isodd(j) ? "↑" : "↓") # Plan to do 5 DMRG sweeps: From cfb04a3b1889763fa0159e099be8c8bf73fd373d Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 21 Feb 2025 14:35:01 -0500 Subject: [PATCH 24/33] Start rewriting random_mps --- Project.toml | 2 +- examples/dmrg/1d_heisenberg.jl | 10 +- src/abstractmps.jl | 24 +- src/mps.jl | 415 ++------------------------------- 4 files changed, 45 insertions(+), 406 deletions(-) diff --git a/Project.toml b/Project.toml index 93c23ce..734f75a 100644 --- a/Project.toml +++ b/Project.toml @@ -49,7 +49,7 @@ DiagonalArrays = "0.2.3" FillArrays = "1.13.0" GradedUnitRanges = "0.1.5" HDF5 = "0.14, 0.15, 0.16, 0.17" -ITensors = "0.8.0" +ITensors = "0.9.0" IsApprox = "2.0.0" KrylovKit = "0.8.1, 0.9" LinearAlgebra = "1.10" diff --git a/examples/dmrg/1d_heisenberg.jl b/examples/dmrg/1d_heisenberg.jl index 3bbe7c0..191cd6c 100644 --- a/examples/dmrg/1d_heisenberg.jl +++ b/examples/dmrg/1d_heisenberg.jl @@ -17,9 +17,9 @@ Random.seed!(1234) N = 10 # Create N spin-one degrees of freedom -sites = siteinds("S=1", N) +#sites = siteinds("S=1", N) # Alternatively can make spin-half sites instead -#sites = siteinds("S=1/2", N) +sites = siteinds("S=1/2", N) os = heisenberg(N) @@ -28,9 +28,9 @@ os = heisenberg(N) H = MPO(os, sites) # Create an initial random matrix product state -psi0 = random_mps(sites; linkdims=10) -# psi0 = random_mps(sites, j -> isodd(j) ? "↑" : "↓"; linkdims=10) -# psi0 = MPS(sites, j -> isodd(j) ? "↑" : "↓") +psi0 = random_mps(sites; maxdim=10) +# psi0 = random_mps(j -> isodd(j) ? "↑" : "↓", sites; linkdims=10) +# psi0 = MPS(j -> isodd(j) ? "↑" : "↓", sites) # Plan to do 5 DMRG sweeps: nsweeps = 5 diff --git a/src/abstractmps.jl b/src/abstractmps.jl index b93eb9c..ff40ad4 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -367,7 +367,11 @@ If there is no link Index, return `nothing`. function linkind(M::AbstractMPS, j::Integer) N = length(M) (j ≥ length(M) || j < 1) && return nothing - return commonind(M[j], M[j + 1]) + is = commoninds(M[j], M[j + 1]) + if iszero(length(is)) + return nothing + end + return only(is) end """ @@ -600,7 +604,7 @@ findsites(M, (s[4]', s[4])) == [4] findsites(M, (s[4]', s[3])) == [3, 4] ``` """ -findsites(ψ::AbstractMPS, is) = findall(hascommoninds(is), ψ) +findsites(ψ::AbstractMPS, is) = findall(ψᵢ -> !isdisjoint(inds(ψᵢ), inds(is)), ψ) findsites(ψ::AbstractMPS, s::Index) = findsites(ψ, IndexSet(s)) @@ -1633,7 +1637,7 @@ function orthogonalize!(M::AbstractMPS, j::Int; maxdim=nothing, normalize=nothin if !isnothing(lb) ltags = tags(lb) else - ltags = TagSet("Link,l=$b") + ltags = Dict("l" => "$b") end L, R = factorize(M[b], linds; tags=ltags, maxdim) M[b] = L @@ -1654,7 +1658,7 @@ function orthogonalize!(M::AbstractMPS, j::Int; maxdim=nothing, normalize=nothin if !isnothing(lb) ltags = tags(lb) else - ltags = TagSet("Link,l=$b") + ltags = Dict("l" => "$b") end L, R = factorize(M[b + 1], rinds; tags=ltags, maxdim) M[b + 1] = L @@ -2121,7 +2125,7 @@ to false. - `move_sites_back::Bool = true`: after the ITensors are applied to the MPS or MPO, move the sites of the MPS or MPO back to their original locations. """ -function product( +function apply( o::ITensor, ψ::AbstractMPS, ns=findsites(ψ, o); @@ -2147,7 +2151,7 @@ function product( for n in 2:N ϕ *= ψ[ns′[n]] end - ϕ = product(o, ϕ; apply_dag=apply_dag) + ϕ = apply(o, ϕ; apply_dag=apply_dag) ψ[ns′[1]:ns′[end], kwargs...] = ϕ if move_sites_back # Move the sites back to their original positions @@ -2254,7 +2258,7 @@ expτH = ops(os, s) ψτ = apply(expτH, ψ0) ``` """ -function product( +function apply( As::Vector{ITensor}, ψ::AbstractMPS; move_sites_back_between_gates::Bool=true, @@ -2263,7 +2267,7 @@ function product( ) Aψ = ψ for A in As - Aψ = product(A, Aψ; move_sites_back=move_sites_back_between_gates, kwargs...) + Aψ = apply(A, Aψ; move_sites_back=move_sites_back_between_gates, kwargs...) end if !move_sites_back_between_gates && move_sites_back s = siteinds(Aψ) @@ -2287,8 +2291,8 @@ end # # # U|ψ⟩ = Z₁X₁|ψ⟩ # apply(U, -function product(o::Prod{ITensor}, ψ::AbstractMPS; kwargs...) - return product(reverse(terms(o)), ψ; kwargs...) +function apply(o::Prod{ITensor}, ψ::AbstractMPS; kwargs...) + return apply(reverse(terms(o)), ψ; kwargs...) end function (o::Prod{ITensor})(ψ::AbstractMPS; kwargs...) diff --git a/src/mps.jl b/src/mps.jl index 9d6f316..d0fef6c 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -50,358 +50,29 @@ function MPS(N::Int; ortho_lims::UnitRange=1:N) return MPS(Vector{ITensor}(undef, N); ortho_lims=ortho_lims) end -""" - MPS([::Type{ElT} = Float64, ]sites; linkdims=1) - -Construct an MPS filled with Empty ITensors of type `ElT` from a collection of indices. - -Optionally specify the link dimension with the keyword argument `linkdims`, which by default is 1. - -In the future we may generalize `linkdims` to allow specifying each individual link dimension as a vector, -and additionally allow specifying quantum numbers. -""" -function MPS( - ::Type{T}, sites::Vector{<:Index}; linkdims::Union{Integer,Vector{<:Integer}}=1 -) where {T<:Number} - _linkdims = _fill_linkdims(linkdims, sites) - N = length(sites) - v = Vector{ITensor}(undef, N) - if N == 1 - v[1] = ITensor(T, sites[1]) - return MPS(v) - end - - spaces = if hasqns(sites) - [[QN() => _linkdims[j]] for j in 1:(N - 1)] - else - [_linkdims[j] for j in 1:(N - 1)] - end - - l = [Index(spaces[ii], "Link,l=$ii") for ii in 1:(N - 1)] - for ii in eachindex(sites) - s = sites[ii] - if ii == 1 - v[ii] = ITensor(T, l[ii], s) - elseif ii == N - v[ii] = ITensor(T, dual(l[ii - 1]), s) - else - v[ii] = ITensor(T, dual(l[ii - 1]), s, l[ii]) - end - end - return MPS(v) -end - -MPS(sites::Vector{<:Index}, args...; kwargs...) = MPS(Float64, sites, args...; kwargs...) - -function randomU(eltype::Type{<:Number}, s1::Index, s2::Index) - return randomU(Random.default_rng(), eltype, s1, s2) -end - -function randomU(rng::AbstractRNG, eltype::Type{<:Number}, s1::Index, s2::Index) - if !hasqns(s1) && !hasqns(s2) - mdim = Int(length(s1)) * Int(length(s2)) - RM = randn(rng, eltype, mdim, mdim) - Q, _ = qr_positive(RM) - inds = (dual(s1), dual(s2), s1', s2') - dims = Int.(length.(inds)) - G = ITensor(reshape(Q, dims), inds) - else - M = random_itensor(rng, eltype, QN(), s1', s2', dual(s1), dual(s2)) - U, S, V = svd(M, (s1', s2')) - u = commonind(U, S) - v = commonind(S, V) - replaceind!(U, u, v) - G = U * V - end - return G -end - -function randomizeMPS!(eltype::Type{<:Number}, M::MPS, sites::Vector{<:Index}, linkdims=1) - return randomizeMPS!(Random.default_rng(), eltype, M, sites, linkdims) -end - -function randomizeMPS!( - rng::AbstractRNG, eltype::Type{<:Number}, M::MPS, sites::Vector{<:Index}, linkdims=1 -) - _linkdims = _fill_linkdims(linkdims, sites) - if isone(length(sites)) - randn!(rng, M[1]) - normalize!(M) - return M - end - N = length(sites) - c = div(N, 2) - max_pass = 100 - for pass in 1:max_pass, half in 1:2 - if half == 1 - (db, brange) = (+1, 1:1:(N - 1)) - else - (db, brange) = (-1, N:-1:2) - end - for b in brange - s1 = sites[b] - s2 = sites[b + db] - G = randomU(rng, eltype, s1, s2) - T = noprime(G * M[b] * M[b + db]) - rinds = uniqueinds(M[b], M[b + db]) - - b_dim = half == 1 ? b : b + db - U, S, V = svd(T, rinds; maxdim=_linkdims[b_dim], utags="Link,l=$(b-1)") - M[b] = U - M[b + db] = S * V - M[b + db] /= norm(M[b + db]) - end - if half == 2 && dim(commonind(M[c], M[c + 1])) >= _linkdims[c] - break - end - end - setleftlim!(M, 0) - setrightlim!(M, 2) - if dim(commonind(M[c], M[c + 1])) < _linkdims[c] - @warn "MPS center bond dimension is less than requested (you requested $(_linkdims[c]), but in practice it is $(dim(commonind(M[c], M[c + 1]))). This is likely due to technicalities of truncating quantum number sectors." - end -end - -function randomCircuitMPS( - eltype::Type{<:Number}, sites::Vector{<:Index}, linkdims::Vector{<:Integer}; kwargs... -) - return randomCircuitMPS(Random.default_rng(), eltype, sites, linkdims; kwargs...) -end - -# TODO: Move this to MatrixAlgebraKit.jl/TensorAlgebra.jl? -using LinearAlgebra: Diagonal, diag -function nonzero_sign(x) - iszero(x) && return one(x) - return sign(x) -end -function qr_positive(M::AbstractMatrix) - Q′, R = qr(M) - Q = convert(typeof(R), Q′) - signs = nonzero_sign.(diag(R)) - Q = Q * Diagonal(signs) - R = Diagonal(conj.(signs)) * R - return Q, R -end -function random_unitary(rng::AbstractRNG, elt::Type, n::Int, m::Int) - if n < m - return copy(random_unitary(rng, elt, m, n)') - end - Q, _ = qr_positive(randn(rng, elt, n, m)) - return Q -end - -function randomCircuitMPS( - rng::AbstractRNG, - eltype::Type{<:Number}, - sites::Vector{<:Index}, - linkdims::Vector{<:Integer}; - kwargs..., -) - N = length(sites) - M = MPS(N) - - if N == 1 - M[1] = ITensor(randn(rng, eltype, dim(sites[1])), sites[1]) - M[1] /= norm(M[1]) - return M - end - - l = Vector{Index}(undef, N) - - d = dim(sites[N]) - chi = min(linkdims[N - 1], d) - l[N - 1] = settag(Index(chi), "l", "$(N-1)") - O = random_unitary(rng, eltype, chi, d) - M[N] = ITensor(O, l[N - 1], sites[N]) - - for j in (N - 1):-1:2 - chi *= dim(sites[j]) - chi = min(linkdims[j - 1], chi) - l[j - 1] = settag(Index(chi), "l", "$(j-1)") - O = random_unitary(rng, eltype, chi, dim(sites[j]) * dim(l[j])) - T = reshape(O, (chi, dim(sites[j]), dim(l[j]))) - M[j] = ITensor(T, l[j - 1], sites[j], l[j]) - end - - O = random_unitary(rng, eltype, 1, dim(sites[1]) * dim(l[1])) - l0 = settag(Index(1), "l", "0") - T = reshape(O, (1, dim(sites[1]), dim(l[1]))) - M[1] = ITensor(T, l0, sites[1], l[1]) - M[1] *= oneelement(eltype, l0 => 1) - - M.llim = 0 - M.rlim = 2 - - return M -end - -function randomCircuitMPS(sites::Vector{<:Index}, linkdims::Vector{<:Integer}; kwargs...) - return randomCircuitMPS(Random.default_rng(), sites, linkdims; kwargs...) -end - -function randomCircuitMPS( - rng::AbstractRNG, sites::Vector{<:Index}, linkdims::Vector{<:Integer}; kwargs... -) - return randomCircuitMPS(rng, Float64, sites, linkdims; kwargs...) -end - -function _fill_linkdims(linkdims::Vector{<:Integer}, sites::Vector{<:Index}) - @assert length(linkdims) == length(sites) - 1 - return linkdims -end - -function _fill_linkdims(linkdims::Integer, sites::Vector{<:Index}) - return fill(linkdims, length(sites) - 1) -end - -""" - random_mps(eltype::Type{<:Number}, sites::Vector{<:Index}; linkdims=1) - -Construct a random MPS with link dimension `linkdims` of -type `eltype`. - -`linkdims` can also accept a `Vector{Int}` with -`length(linkdims) == length(sites) - 1` for constructing an -MPS with non-uniform bond dimension. -""" -function random_mps( - ::Type{ElT}, sites::Vector{<:Index}; linkdims::Union{Integer,Vector{<:Integer}}=1 -) where {ElT<:Number} - return random_mps(Random.default_rng(), ElT, sites; linkdims) -end - -function random_mps( - rng::AbstractRNG, - ::Type{ElT}, - sites::Vector{<:Index}; - linkdims::Union{Integer,Vector{<:Integer}}=1, -) where {ElT<:Number} - _linkdims = _fill_linkdims(linkdims, sites) - if any(hasqns, sites) - error("initial state required to use random_mps with QNs") +function apply_random_staircase_circuit(rng::AbstractRNG, elt::Type, state::MPS; depth, kwargs...) + n = length(state) + layer = [(j - 1, j) for j in reverse(2:n)] + s = siteinds(state) + gate_layers = mapreduce(vcat, 1:depth) do _ + # TODO: Pass `elt` and `rng` as kwargs to `op`, to be + # used as parameters in `OpName` by `QuantumOperaterDefinitions.jl`. + return map(((i, j),) -> op("RandomUnitary", (s[i], s[j])), layer) end - - # For non-QN-conserving MPS, instantiate - # the random MPS directly as a circuit: - return randomCircuitMPS(rng, ElT, sites, _linkdims) -end - -""" - random_mps(sites::Vector{<:Index}; linkdims=1) - random_mps(eltype::Type{<:Number}, sites::Vector{<:Index}; linkdims=1) - -Construct a random MPS with link dimension `linkdims` which by -default has element type `Float64`. - -`linkdims` can also accept a `Vector{Int}` with -`length(linkdims) == length(sites) - 1` for constructing an -MPS with non-uniform bond dimension. -""" -function random_mps(sites::Vector{<:Index}; linkdims::Union{Integer,Vector{<:Integer}}=1) - return random_mps(Random.default_rng(), sites; linkdims) + # TODO: Use `apply`, for some reason that can't be found right now. + return apply(gate_layers, state; kwargs...) end -function random_mps( - rng::AbstractRNG, sites::Vector{<:Index}; linkdims::Union{Integer,Vector{<:Integer}}=1 -) - return random_mps(rng, Float64, sites; linkdims) -end - -function random_mps( - sites::Vector{<:Index}, state; linkdims::Union{Integer,Vector{<:Integer}}=1 -) - return random_mps(Random.default_rng(), sites, state; linkdims) +function random_mps(rng::AbstractRNG, elt::Type, state, sites::Vector{<:Index}; maxdim, kwargs...) + x = MPS(elt, state, sites) + depth = ceil(Int, log(minimum(Int ∘ length, sites), maxdim)) + return apply_random_staircase_circuit(rng, elt, x; depth, maxdim, kwargs...) end - -function random_mps( - rng::AbstractRNG, - sites::Vector{<:Index}, - state; - linkdims::Union{Integer,Vector{<:Integer}}=1, -) - return random_mps(rng, Float64, sites, state; linkdims) -end - -function random_mps( - eltype::Type{<:Number}, - sites::Vector{<:Index}, - state; - linkdims::Union{Integer,Vector{<:Integer}}=1, -) - return random_mps(Random.default_rng(), eltype, sites, state; linkdims) -end - -function random_mps( - rng::AbstractRNG, - eltype::Type{<:Number}, - sites::Vector{<:Index}, - state; - linkdims::Union{Integer,Vector{<:Integer}}=1, -)::MPS - M = MPS(eltype, sites, state) - if any(>(1), linkdims) - randomizeMPS!(rng, eltype, M, sites, linkdims) - end - return M -end - -@doc """ - random_mps(sites::Vector{<:Index}, state; linkdims=1) - -Construct a real, random MPS with link dimension `linkdims`, -made by randomizing an initial product state specified by -`state`. This version of `random_mps` is necessary when creating -QN-conserving random MPS (consisting of QNITensors). The initial -`state` array provided determines the total QN of the resulting -random MPS. -""" random_mps(::Vector{<:Index}, ::Any) - -""" - MPS(::Type{T<:Number}, ivals::Vector{<:Pair{<:Index}}) - -Construct a product state MPS with element type `T` and -nonzero values determined from the input IndexVals. -""" -function MPS(::Type{T}, ivals::Vector{<:Pair{<:Index}}) where {T<:Number} - N = length(ivals) - M = MPS(N) - - if N == 1 - M[1] = ITensor(T, ind(ivals[1])) - M[1][ivals[1]] = one(T) - return M - end - - if hasqns(ind(ivals[1])) - lflux = QN() - for j in 1:(N - 1) - lflux += qn(ivals[j]) - end - links = Vector{QNIndex}(undef, N - 1) - for j in (N - 1):-1:1 - links[j] = dual(Index(lflux => 1; tags="Link,l=$j")) - lflux -= qn(ivals[j]) - end - else - links = [Index(1, "Link,l=$n") for n in 1:(N - 1)] - end - - M[1] = ITensor(T, ind(ivals[1]), links[1]) - M[1][ivals[1], links[1] => 1] = one(T) - for n in 2:(N - 1) - s = ind(ivals[n]) - M[n] = ITensor(T, dual(links[n - 1]), s, links[n]) - M[n][links[n - 1] => 1, ivals[n], links[n] => 1] = one(T) - end - M[N] = ITensor(T, dual(links[N - 1]), ind(ivals[N])) - M[N][links[N - 1] => 1, ivals[N]] = one(T) - - return M +function random_mps(sites; kwargs...) + state = fill("0", length(sites)) + return random_mps(Random.default_rng(), Float64, state, sites; kwargs...) end -# For backwards compatibility -const productMPS = MPS - """ MPS(ivals::Vector{<:Pair{<:Index}}) @@ -436,59 +107,23 @@ psi = MPS(ComplexF64, sites, states) phi = MPS(sites, "Up") ``` """ -function MPS(eltype::Type{<:Number}, sites::Vector{<:Index}, states_) - if length(sites) != length(states_) - throw(DimensionMismatch("Number of sites and and initial vals don't match")) - end - N = length(states_) - M = MPS(N) - - if N == 1 - M[1] = state(states_[1], sites[1]) - return convert_leaf_eltype(eltype, M) - end - - states = [state(states_[j], sites[j]) for j in 1:N] - - if hasqns(states[1]) - lflux = QN() - for j in 1:(N - 1) - lflux += flux(states[j]) - end - links = Vector{QNIndex}(undef, N - 1) - for j in (N - 1):-1:1 - links[j] = dual(settag(Index(lflux => 1), "l", "$j")) - lflux -= flux(states[j]) - end - else - links = [settag(Index(1), "l", "$n") for n in 1:N] - end - - M[1] = ITensor(sites[1], links[1]) - M[1] += states[1] * state(1, links[1]) - for n in 2:(N - 1) - M[n] = ITensor(dual(links[n - 1]), sites[n], links[n]) - M[n] += state(1, dual(links[n - 1])) * states[n] * state(1, links[n]) - end - M[N] = ITensor(dual(links[N - 1]), sites[N]) - M[N] += state(1, dual(links[N - 1])) * states[N] - - return convert_leaf_eltype(eltype, M) +function MPS(elt::Type{<:Number}, states_, sites::Vector{<:Index}) + return MPS([state(elt, states_[j], sites[j]) for j in 1:length(sites)]) end function MPS( - ::Type{T}, sites::Vector{<:Index}, state::Union{String,Integer} + ::Type{T}, state::Union{String,Integer}, sites::Vector{<:Index} ) where {T<:Number} - return MPS(T, sites, fill(state, length(sites))) + return MPS(T, fill(state, length(sites)), sites) end -function MPS(::Type{T}, sites::Vector{<:Index}, states::Function) where {T<:Number} +function MPS(::Type{T}, states::Function, sites::Vector{<:Index}) where {T<:Number} states_vec = [states(n) for n in 1:length(sites)] - return MPS(T, sites, states_vec) + return MPS(T, states_vec, sites) end """ - MPS(sites::Vector{<:Index},states) + MPS(states, sites::Vector{<:Index}) Construct a product state MPS having site indices `sites`, and which corresponds to the initial @@ -506,7 +141,7 @@ states = [isodd(n) ? "Up" : "Dn" for n in 1:N] psi = MPS(sites, states) ``` """ -MPS(sites::Vector{<:Index}, states) = MPS(Float64, sites, states) +MPS(state, sites::Vector{<:Index}) = MPS(Float64, state, sites) """ siteind(M::MPS, j::Int; kwargs...) From 4f7f0b24e9f334869535aa3e54763c54e9990c65 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 21 Feb 2025 18:12:42 -0500 Subject: [PATCH 25/33] Start adding apply/product for ITensors --- src/abstractmps.jl | 47 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/src/abstractmps.jl b/src/abstractmps.jl index ff40ad4..c97e43e 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -2101,6 +2101,53 @@ function movesites(ψ::AbstractMPS, ns, ns′; kwargs...) return ψ end +# TODO: Fix this, move to `ITensorBase.jl`, rewrite in terms of +# a new `ITensorMap` type to generalize beyond prime levels. +using ITensors: setprime +function apply(A::ITensor, B::ITensor; apply_dag::Bool=false) + commonindsAB = filter(i -> iszero(plev(i)), commoninds(A, B)) + isempty(commonindsAB) && error("In product, must have common indices with prime level 0.") + common_paired_indsA = filter( + i -> ((i ∈ commonindsAB) && (setprime(i, 1) ∈ inds(A))), inds(A) + ) + common_paired_indsB = filter( + i -> ((i ∈ commonindsAB) && (setprime(i, 1) ∈ inds(B))), inds(B) + ) + + if !isempty(common_paired_indsA) + commoninds_pairs = unioninds(common_paired_indsA, common_paired_indsA') + elseif !isempty(common_paired_indsB) + commoninds_pairs = unioninds(common_paired_indsB, common_paired_indsB') + else + # vector-vector product + apply_dag && error("apply_dag not supported for vector-vector product") + return A * B + end + danglings_indsA = uniqueinds(A, commoninds_pairs) + danglings_indsB = uniqueinds(B, commoninds_pairs) + danglings_inds = unioninds(danglings_indsA, danglings_indsB) + if hassameinds(common_paired_indsA, common_paired_indsB) + # matrix-matrix product + A′ = prime(A; inds=!danglings_inds) + AB = mapprime(A′ * B, 2 => 1; inds=!danglings_inds) + if apply_dag + AB′ = prime(AB; inds=!danglings_inds) + Adag = swapprime(dag(A), 0 => 1; inds=!danglings_inds) + return mapprime(AB′ * Adag, 2 => 1; inds=!danglings_inds) + end + return AB + elseif isempty(common_paired_indsA) && !isempty(common_paired_indsB) + # vector-matrix product + apply_dag && error("apply_dag not supported for matrix-vector product") + A′ = prime(A; inds=!danglings_inds) + return A′ * B + elseif !isempty(common_paired_indsA) && isempty(common_paired_indsB) + # matrix-vector product + apply_dag && error("apply_dag not supported for vector-matrix product") + return replaceprime(A * B, 1 => 0; inds=!danglings_inds) + end +end + """ apply(o::ITensor, ψ::Union{MPS, MPO}, [ns::Vector{Int}]; kwargs...) product([...]) From c8eec73aa80a61f754063bf6338f5d0b44ee07f2 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 22 Feb 2025 10:51:07 -0500 Subject: [PATCH 26/33] Start updating apply --- src/abstractmps.jl | 117 ++++++++++++++++++++++++++------------------- 1 file changed, 68 insertions(+), 49 deletions(-) diff --git a/src/abstractmps.jl b/src/abstractmps.jl index c97e43e..bed08bf 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -1872,49 +1872,50 @@ function setindex!( #s = collect(Iterators.flatten(sites)) indsA = filter(x -> !isnothing(x), [lind, Iterators.flatten(sites)..., rind]) - @assert hassameinds(A, indsA) + @assert issetequal(inds(A), indsA) # For MPO case, restrict to 0 prime level #sites = filter(hasplev(0), sites) - if !isnothing(perm) - sites0 = sites - sites = sites0[[perm...]] - # Check if the site indices - # are fermionic - if !using_auto_fermion() && any(ITensors.anyfermionic, sites) - if length(sites) == 2 && ψ isa MPS - if all(ITensors.allfermionic, sites) - s0 = Index.(sites0) - - # TODO: the Fermionic swap is could be diagonal, - # if we combine the site indices - #C = combiner(s0[1], s0[2]) - #c = combinedind(C) - #AC = A * C - #AC = noprime(AC * _fermionic_swap(c)) - #A = AC * dag(C) - - FSWAP = adapt(datatype(A), _fermionic_swap(s0[1], s0[2])) - A = noprime(A * FSWAP) - end - elseif ψ isa MPO - @warn "In setindex!(MPO, ::ITensor, ::UnitRange), " * - "fermionic signs are only not handled properly for non-trivial " * - "permutations of sites. Please inform the developers of ITensors " * - "if you require this feature (otherwise, fermionic signs can be " * - "put in manually with fermionic swap gates)." - else - @warn "In setindex!(::Union{MPS, MPO}, ::ITensor, ::UnitRange), " * - "fermionic signs are only handled properly for permutations involving 2 sites. " * - "The original sites are $sites0, with a permutation $perm. " * - "To have the fermion sign handled correctly, we recommend performing your permutation " * - "pairwise." - end - end - end + ## TODO: Bring this package to handle fermions. + ## if !isnothing(perm) + ## sites0 = sites + ## sites = sites0[[perm...]] + ## # Check if the site indices + ## # are fermionic + ## if !using_auto_fermion() && any(ITensors.anyfermionic, sites) + ## if length(sites) == 2 && ψ isa MPS + ## if all(ITensors.allfermionic, sites) + ## s0 = Index.(sites0) + + ## # TODO: the Fermionic swap is could be diagonal, + ## # if we combine the site indices + ## #C = combiner(s0[1], s0[2]) + ## #c = combinedind(C) + ## #AC = A * C + ## #AC = noprime(AC * _fermionic_swap(c)) + ## #A = AC * dag(C) + + ## FSWAP = adapt(datatype(A), _fermionic_swap(s0[1], s0[2])) + ## A = noprime(A * FSWAP) + ## end + ## elseif ψ isa MPO + ## @warn "In setindex!(MPO, ::ITensor, ::UnitRange), " * + ## "fermionic signs are only not handled properly for non-trivial " * + ## "permutations of sites. Please inform the developers of ITensors " * + ## "if you require this feature (otherwise, fermionic signs can be " * + ## "put in manually with fermionic swap gates)." + ## else + ## @warn "In setindex!(::Union{MPS, MPO}, ::ITensor, ::UnitRange), " * + ## "fermionic signs are only handled properly for permutations involving 2 sites. " * + ## "The original sites are $sites0, with a permutation $perm. " * + ## "To have the fermion sign handled correctly, we recommend performing your permutation " * + ## "pairwise." + ## end + ## end + ## end - ψA = MPST(A, sites; leftinds=lind, orthocenter=orthocenter - first(r) + 1, kwargs...) + ψA = MPST(A, sites; leftinds=(lind,), orthocenter=orthocenter - first(r) + 1, kwargs...) #@assert prod(ψA) ≈ A ψ[firstsite:lastsite] = ψA @@ -1959,6 +1960,7 @@ function (::Type{MPST})( for s in sites @assert s ⊆ inds(A) end + @assert isnothing(leftinds) || leftinds ⊆ inds(A) @assert 1 ≤ orthocenter ≤ N @@ -1970,9 +1972,9 @@ function (::Type{MPST})( # 1:orthocenter and reverse(orthocenter:N) # so the orthogonality center is set correctly. for n in 1:(N - 1) - Lis = IndexSet(sites[n]) + Lis = (sites[n],) if !isnothing(l) - Lis = unioninds(Lis, l) + Lis = Lis ∪ l end L, R = factorize(Ã, Lis; kwargs..., tags="Link,n=$n", ortho="left") l = commonind(L, R) @@ -2115,21 +2117,31 @@ function apply(A::ITensor, B::ITensor; apply_dag::Bool=false) ) if !isempty(common_paired_indsA) - commoninds_pairs = unioninds(common_paired_indsA, common_paired_indsA') + ## Old version: `commoninds_pairs = unioninds(common_paired_indsA, common_paired_indsA')` + commoninds_pairs = common_paired_indsA ∪ prime.(common_paired_indsA) elseif !isempty(common_paired_indsB) - commoninds_pairs = unioninds(common_paired_indsB, common_paired_indsB') + ## Old version: `commoninds_pairs = unioninds(common_paired_indsB, common_paired_indsB')` + commoninds_pairs = common_paired_indsB ∪ prime.(common_paired_indsB) else # vector-vector product apply_dag && error("apply_dag not supported for vector-vector product") return A * B end - danglings_indsA = uniqueinds(A, commoninds_pairs) - danglings_indsB = uniqueinds(B, commoninds_pairs) - danglings_inds = unioninds(danglings_indsA, danglings_indsB) - if hassameinds(common_paired_indsA, common_paired_indsB) + danglings_indsA = setdiff(inds(A), commoninds_pairs) + danglings_indsB = setdiff(inds(B), commoninds_pairs) + danglings_inds = danglings_indsA ∪ danglings_indsB + if issetequal(common_paired_indsA, common_paired_indsB) # matrix-matrix product - A′ = prime(A; inds=!danglings_inds) - AB = mapprime(A′ * B, 2 => 1; inds=!danglings_inds) + ## Old version: `prime(A; inds=!danglings_inds)` + A′ = replaceinds(i -> i ∉ danglings_inds ? i' : i, A) + A′B = A′ * B + ## Old version: `AB = mapprime(A′B, 2 => 1; inds=!danglings_inds)` + A′ = replaceinds(A) do i + if i ∉ danglings_inds && plev(i) == 2 + return setprime(i, 1) + end + return i + end if apply_dag AB′ = prime(AB; inds=!danglings_inds) Adag = swapprime(dag(A), 0 => 1; inds=!danglings_inds) @@ -2144,7 +2156,14 @@ function apply(A::ITensor, B::ITensor; apply_dag::Bool=false) elseif !isempty(common_paired_indsA) && isempty(common_paired_indsB) # matrix-vector product apply_dag && error("apply_dag not supported for vector-matrix product") - return replaceprime(A * B, 1 => 0; inds=!danglings_inds) + ## Old version: `replaceprime(A * B, 1 => 0; inds=!danglings_inds)` + AB = A * B + return replaceinds(AB) do i + if i ∉ danglings_inds && plev(i) == 1 + return setprime(i, 0) + end + return i + end end end From 419a1c3c7019c5b91af1647053ccf82573d15558 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Mon, 24 Feb 2025 10:28:41 -0500 Subject: [PATCH 27/33] Gate application and new random_mps implementation --- examples/dmrg/1d_heisenberg.jl | 6 +-- examples/dmrg/1d_heisenberg_conserve_spin.jl | 4 +- src/abstractmps.jl | 39 +++++++------------- src/mps.jl | 17 ++++++++- 4 files changed, 34 insertions(+), 32 deletions(-) diff --git a/examples/dmrg/1d_heisenberg.jl b/examples/dmrg/1d_heisenberg.jl index 191cd6c..5cc6c15 100644 --- a/examples/dmrg/1d_heisenberg.jl +++ b/examples/dmrg/1d_heisenberg.jl @@ -17,9 +17,9 @@ Random.seed!(1234) N = 10 # Create N spin-one degrees of freedom -#sites = siteinds("S=1", N) +sites = siteinds("S=1", N) # Alternatively can make spin-half sites instead -sites = siteinds("S=1/2", N) +# sites = siteinds("S=1/2", N) os = heisenberg(N) @@ -29,7 +29,7 @@ H = MPO(os, sites) # Create an initial random matrix product state psi0 = random_mps(sites; maxdim=10) -# psi0 = random_mps(j -> isodd(j) ? "↑" : "↓", sites; linkdims=10) +# psi0 = random_mps(j -> isodd(j) ? "↑" : "↓", sites; maxdim=10) # psi0 = MPS(j -> isodd(j) ? "↑" : "↓", sites) # Plan to do 5 DMRG sweeps: diff --git a/examples/dmrg/1d_heisenberg_conserve_spin.jl b/examples/dmrg/1d_heisenberg_conserve_spin.jl index ccf54d0..b854053 100644 --- a/examples/dmrg/1d_heisenberg_conserve_spin.jl +++ b/examples/dmrg/1d_heisenberg_conserve_spin.jl @@ -23,8 +23,8 @@ s = siteinds("S=1/2", N; gradings=("Sz",)) os = heisenberg(N) # Create an initial random matrix product state -psi0 = random_mps(s, j -> isodd(j) ? "↑" : "↓"; linkdims=10) -# psi0 = MPS(s, j -> isodd(j) ? "↑" : "↓") +psi0 = random_mps(j -> isodd(j) ? "↑" : "↓", s; maxdim=10) +# psi0 = MPS(j -> isodd(j) ? "↑" : "↓", s) # Input operator terms which define a Hamiltonian # Convert these terms to an MPO tensor network diff --git a/src/abstractmps.jl b/src/abstractmps.jl index bed08bf..edf0d02 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -383,7 +383,7 @@ MPS or MPO tensor on site j to site j+1. """ function linkinds(M::AbstractMPS, j::Integer) N = length(M) - (j ≥ length(M) || j < 1) && return IndexSet() + (j ≥ length(M) || j < 1) && return [] return commoninds(M[j], M[j + 1]) end @@ -1637,7 +1637,7 @@ function orthogonalize!(M::AbstractMPS, j::Int; maxdim=nothing, normalize=nothin if !isnothing(lb) ltags = tags(lb) else - ltags = Dict("l" => "$b") + ltags = ["l" => "$b"] end L, R = factorize(M[b], linds; tags=ltags, maxdim) M[b] = L @@ -1647,9 +1647,7 @@ function orthogonalize!(M::AbstractMPS, j::Int; maxdim=nothing, normalize=nothin setrightlim!(M, leftlim(M) + 2) end end - N = length(M) - while rightlim(M) > (j + 1) (rightlim(M) > (N + 1)) && setrightlim!(M, N + 1) b = rightlim(M) - 2 @@ -1658,12 +1656,11 @@ function orthogonalize!(M::AbstractMPS, j::Int; maxdim=nothing, normalize=nothin if !isnothing(lb) ltags = tags(lb) else - ltags = Dict("l" => "$b") + ltags = ["l" => "$b"] end L, R = factorize(M[b + 1], rinds; tags=ltags, maxdim) M[b + 1] = L M[b] *= R - setrightlim!(M, b + 1) if leftlim(M) > rightlim(M) - 2 setleftlim!(M, rightlim(M) - 2) @@ -1865,14 +1862,11 @@ function setindex!( # Check that A has the proper common # indices with ψ - lind = linkind(ψ, firstsite - 1) - rind = linkind(ψ, lastsite) - + linds = linkinds(ψ, firstsite - 1) + rinds = linkinds(ψ, lastsite) sites = [siteinds(ψ, j) for j in firstsite:lastsite] - #s = collect(Iterators.flatten(sites)) - indsA = filter(x -> !isnothing(x), [lind, Iterators.flatten(sites)..., rind]) - @assert issetequal(inds(A), indsA) + @assert issetequal(inds(A), [linds; reduce(vcat, sites); rinds]) # For MPO case, restrict to 0 prime level #sites = filter(hasplev(0), sites) @@ -1915,7 +1909,7 @@ function setindex!( ## end ## end - ψA = MPST(A, sites; leftinds=(lind,), orthocenter=orthocenter - first(r) + 1, kwargs...) + ψA = MPST(A, sites; leftinds=linds, orthocenter=orthocenter - first(r) + 1, kwargs...) #@assert prod(ψA) ≈ A ψ[firstsite:lastsite] = ψA @@ -1946,7 +1940,7 @@ by site according to the site indices `sites`. # Keywords -- `leftinds = nothing`: optional left dangling indices. Indices that are not +- `leftinds = []`: optional left dangling indices. Indices that are not in `sites` and `leftinds` will be dangling off of the right side of the MPS/MPO. - `orthocenter::Integer = length(sites)`: the desired final orthogonality center of the output MPS/MPO. @@ -1954,15 +1948,13 @@ by site according to the site indices `sites`. - `maxdim`: the maximum link dimension. """ function (::Type{MPST})( - A::ITensor, sites; leftinds=nothing, orthocenter::Integer=length(sites), kwargs... + A::ITensor, sites; leftinds=[], orthocenter::Integer=length(sites), kwargs... ) where {MPST<:AbstractMPS} N = length(sites) for s in sites @assert s ⊆ inds(A) end - - @assert isnothing(leftinds) || leftinds ⊆ inds(A) - + @assert leftinds ⊆ inds(A) @assert 1 ≤ orthocenter ≤ N ψ = Vector{ITensor}(undef, N) @@ -1972,12 +1964,9 @@ function (::Type{MPST})( # 1:orthocenter and reverse(orthocenter:N) # so the orthogonality center is set correctly. for n in 1:(N - 1) - Lis = (sites[n],) - if !isnothing(l) - Lis = Lis ∪ l - end - L, R = factorize(Ã, Lis; kwargs..., tags="Link,n=$n", ortho="left") - l = commonind(L, R) + Lis = sites[n] ∪ l + L, R = factorize(Ã, Lis; kwargs..., tags=["site" => "$n"], ortho="left") + l = commoninds(L, R) ψ[n] = L Ã = R end @@ -2136,7 +2125,7 @@ function apply(A::ITensor, B::ITensor; apply_dag::Bool=false) A′ = replaceinds(i -> i ∉ danglings_inds ? i' : i, A) A′B = A′ * B ## Old version: `AB = mapprime(A′B, 2 => 1; inds=!danglings_inds)` - A′ = replaceinds(A) do i + AB = replaceinds(A′B) do i if i ∉ danglings_inds && plev(i) == 2 return setprime(i, 1) end diff --git a/src/mps.jl b/src/mps.jl index d0fef6c..083c7a0 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -50,7 +50,9 @@ function MPS(N::Int; ortho_lims::UnitRange=1:N) return MPS(Vector{ITensor}(undef, N); ortho_lims=ortho_lims) end -function apply_random_staircase_circuit(rng::AbstractRNG, elt::Type, state::MPS; depth, kwargs...) +function apply_random_staircase_circuit( + rng::AbstractRNG, elt::Type, state::MPS; depth, kwargs... +) n = length(state) layer = [(j - 1, j) for j in reverse(2:n)] s = siteinds(state) @@ -63,11 +65,22 @@ function apply_random_staircase_circuit(rng::AbstractRNG, elt::Type, state::MPS; return apply(gate_layers, state; kwargs...) end -function random_mps(rng::AbstractRNG, elt::Type, state, sites::Vector{<:Index}; maxdim, kwargs...) +function random_mps( + rng::AbstractRNG, elt::Type, state, sites::Vector{<:Index}; maxdim, kwargs... +) x = MPS(elt, state, sites) depth = ceil(Int, log(minimum(Int ∘ length, sites), maxdim)) return apply_random_staircase_circuit(rng, elt, x; depth, maxdim, kwargs...) end +function random_mps(rng::AbstractRNG, state, sites::Vector{<:Index}; kwargs...) + return random_mps(rng, Float64, state, sites; kwargs...) +end +function random_mps(elt::Type, state, sites::Vector{<:Index}; kwargs...) + return random_mps(Random.default_rng(), elt, state, sites; kwargs...) +end +function random_mps(state, sites::Vector{<:Index}; kwargs...) + return random_mps(Random.default_rng(), Float64, state, sites; kwargs...) +end function random_mps(sites; kwargs...) state = fill("0", length(sites)) return random_mps(Random.default_rng(), Float64, state, sites; kwargs...) From e536fd94707f0fd95ad0e93554da40d02dd30eb1 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Mon, 10 Mar 2025 11:50:03 -0400 Subject: [PATCH 28/33] Bump package versions --- Project.toml | 26 +++++++++++++------------- src/mps.jl | 2 +- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/Project.toml b/Project.toml index 734f75a..890de3a 100644 --- a/Project.toml +++ b/Project.toml @@ -41,32 +41,32 @@ ITensorMPSPackageCompilerExt = "PackageCompiler" ITensorMPSZygoteRulesExt = ["ChainRulesCore", "ZygoteRules"] [compat] -Adapt = "4.1.0" -BackendSelection = "0.1.0" +Adapt = "4.1" +BackendSelection = "0.1" ChainRulesCore = "1.10" Compat = "4.16.0" -DiagonalArrays = "0.2.3" -FillArrays = "1.13.0" +DiagonalArrays = "0.3" +FillArrays = "1.13" GradedUnitRanges = "0.1.5" HDF5 = "0.14, 0.15, 0.16, 0.17" -ITensors = "0.9.0" -IsApprox = "2.0.0" +ITensors = "0.9" +IsApprox = "2.0" KrylovKit = "0.8.1, 0.9" LinearAlgebra = "1.10" -NamedDimsArrays = "0.4" +NamedDimsArrays = "0.5" Observers = "0.2" PackageCompiler = "1, 2" Printf = "1.10" -QuantumOperatorAlgebra = "0.1.0" +QuantumOperatorAlgebra = "0.1" QuantumOperatorDefinitions = "0.1.2" Random = "1.10" -SerializedElementArrays = "0.1.0" -SparseArraysBase = "0.2.11" -TensorAlgebra = "0.1.7" +SerializedElementArrays = "0.1" +SparseArraysBase = "0.5" +TensorAlgebra = "0.2" Test = "1.10" -TupleTools = "1.5.0" +TupleTools = "1.5" TypeParameterAccessors = "0.3" -VectorInterface = "0.5.0" +VectorInterface = "0.5" ZygoteRules = "0.2.2" julia = "1.10" diff --git a/src/mps.jl b/src/mps.jl index 083c7a0..a2d37bd 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -59,7 +59,7 @@ function apply_random_staircase_circuit( gate_layers = mapreduce(vcat, 1:depth) do _ # TODO: Pass `elt` and `rng` as kwargs to `op`, to be # used as parameters in `OpName` by `QuantumOperaterDefinitions.jl`. - return map(((i, j),) -> op("RandomUnitary", (s[i], s[j])), layer) + return map(((i, j),) -> op("RandomUnitary", (s[i], s[j]); rng, eltype=elt), layer) end # TODO: Use `apply`, for some reason that can't be found right now. return apply(gate_layers, state; kwargs...) From f0b87172f5cf0d1a5c8e6371370c649e8605e5a3 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Wed, 9 Apr 2025 10:30:38 -0400 Subject: [PATCH 29/33] Remove stale example deps --- examples/dmrg/Project.toml | 3 --- 1 file changed, 3 deletions(-) diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index a5c923c..83d9e46 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -1,10 +1,7 @@ [deps] BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" -GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" -SymmetrySectors = "f8a8ad64-adbc-4fce-92f7-ffe2bb36a86e" -TypeParameterAccessors = "7e5a90cf-f82e-492e-a09b-e3e26432c138" From 7359f70fb78460b3d468dc972fc21d090640cf05 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 11 Apr 2025 15:49:54 -0400 Subject: [PATCH 30/33] Update for dependent packages --- Project.toml | 10 +++++----- examples/dmrg/Project.toml | 2 +- src/abstractmps.jl | 16 ++++++++-------- src/abstractprojmpo/abstractprojmpo.jl | 2 +- src/abstractprojmpo/projmps.jl | 2 +- src/lib/Experimental/src/dmrg.jl | 2 +- src/mpo.jl | 2 +- src/mps.jl | 2 +- src/solvers/alternating_update.jl | 2 +- src/solvers/dmrg_x.jl | 2 +- src/solvers/expand.jl | 2 +- src/solvers/reducedconstantterm.jl | 2 +- src/solvers/sweep_update.jl | 4 ++-- src/solvers/tdvp.jl | 2 +- test/base/test_mpo.jl | 4 ++-- test/base/test_mps.jl | 2 +- test/base/test_qnmpo.jl | 1 + test/base/test_solvers/test_tdvp.jl | 2 +- .../test_optimization.jl | 2 +- 19 files changed, 32 insertions(+), 31 deletions(-) diff --git a/Project.toml b/Project.toml index 890de3a..7cc2d5d 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,7 @@ BackendSelection = "680c2d7c-f67a-4cc9-ae9c-da132b1447a5" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" -GradedUnitRanges = "e2de450a-8a67-46c7-b59c-01d5a3d041c5" +GradedArrays = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" IsApprox = "28f27b66-4bd8-47e7-9110-e2746eb8bed7" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" @@ -47,22 +47,22 @@ ChainRulesCore = "1.10" Compat = "4.16.0" DiagonalArrays = "0.3" FillArrays = "1.13" -GradedUnitRanges = "0.1.5" +GradedArrays = "0.3.1" HDF5 = "0.14, 0.15, 0.16, 0.17" ITensors = "0.9" IsApprox = "2.0" KrylovKit = "0.8.1, 0.9" LinearAlgebra = "1.10" -NamedDimsArrays = "0.5" +NamedDimsArrays = "0.7" Observers = "0.2" PackageCompiler = "1, 2" Printf = "1.10" QuantumOperatorAlgebra = "0.1" -QuantumOperatorDefinitions = "0.1.2" +QuantumOperatorDefinitions = "0.2" Random = "1.10" SerializedElementArrays = "0.1" SparseArraysBase = "0.5" -TensorAlgebra = "0.2" +TensorAlgebra = "0.3" Test = "1.10" TupleTools = "1.5" TypeParameterAccessors = "0.3" diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index 83d9e46..0b09591 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -1,7 +1,7 @@ [deps] -BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" diff --git a/src/abstractmps.jl b/src/abstractmps.jl index 3825b28..fa421c5 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -1,5 +1,5 @@ using BackendSelection: @Algorithm_str, Algorithm -using GradedUnitRanges: GradedUnitRanges, dag +using GradedArrays: GradedUnitRanges, dag using IsApprox: Approx, IsApprox using ITensors: ITensors, @@ -714,7 +714,7 @@ end # TODO: change kwarg from `set_limits` to `preserve_ortho` function Base.map!(f::Function, M::AbstractMPS; set_limits::Bool=true) for i in eachindex(M) - M[i, set_limits=set_limits] = f(M[i]) + M[i, set_limits = set_limits] = f(M[i]) end return M end @@ -1328,7 +1328,7 @@ of the orthogonality center to avoid numerical overflow in the case of diverging See also [`normalize!`](@ref), [`norm`](@ref), [`lognorm`](@ref). """ function LinearAlgebra.normalize(M::AbstractMPS; (lognorm!)=[]) - return normalize!(deepcopy_ortho_center(M); (lognorm!)=lognorm!) + return normalize!(deepcopy_ortho_center(M); (lognorm!)=(lognorm!)) end """ @@ -2014,7 +2014,7 @@ function swapbondsites(ψ::AbstractMPS, b::Integer; ortho="right", kwargs...) elseif rightlim(ψ) > b + 2 ψ = orthogonalize(ψ, b + 1) end - ψ[b:(b + 1), orthocenter=orthocenter, perm=[2, 1], kwargs...] = ψ[b] * ψ[b + 1] + ψ[b:(b + 1), orthocenter = orthocenter, perm = [2, 1], kwargs...] = ψ[b] * ψ[b + 1] return ψ end @@ -2147,15 +2147,15 @@ function apply(A::ITensor, B::ITensor; apply_dag::Bool=false) return i end if apply_dag - AB′ = prime(AB; inds=!danglings_inds) - Adag = swapprime(dag(A), 0 => 1; inds=!danglings_inds) - return mapprime(AB′ * Adag, 2 => 1; inds=!danglings_inds) + AB′ = prime(AB; inds=(!danglings_inds)) + Adag = swapprime(dag(A), 0 => 1; inds=(!danglings_inds)) + return mapprime(AB′ * Adag, 2 => 1; inds=(!danglings_inds)) end return AB elseif isempty(common_paired_indsA) && !isempty(common_paired_indsB) # vector-matrix product apply_dag && error("apply_dag not supported for matrix-vector product") - A′ = prime(A; inds=!danglings_inds) + A′ = prime(A; inds=(!danglings_inds)) return A′ * B elseif !isempty(common_paired_indsA) && isempty(common_paired_indsB) # matrix-vector product diff --git a/src/abstractprojmpo/abstractprojmpo.jl b/src/abstractprojmpo/abstractprojmpo.jl index 38d9532..8dd5a87 100644 --- a/src/abstractprojmpo/abstractprojmpo.jl +++ b/src/abstractprojmpo/abstractprojmpo.jl @@ -1,5 +1,5 @@ using FillArrays: Ones -using GradedUnitRanges: dag +using GradedArrays: dag using ITensors: noprime abstract type AbstractProjMPO end diff --git a/src/abstractprojmpo/projmps.jl b/src/abstractprojmpo/projmps.jl index c29db74..5ff2785 100644 --- a/src/abstractprojmpo/projmps.jl +++ b/src/abstractprojmpo/projmps.jl @@ -1,4 +1,4 @@ -using GradedUnitRanges: dag +using GradedArrays: dag mutable struct ProjMPS lpos::Int diff --git a/src/lib/Experimental/src/dmrg.jl b/src/lib/Experimental/src/dmrg.jl index 8796936..6739c13 100644 --- a/src/lib/Experimental/src/dmrg.jl +++ b/src/lib/Experimental/src/dmrg.jl @@ -10,7 +10,7 @@ function dmrg( operator, init::MPS; updater=eigsolve_updater, (observer!)=default_observer(), kwargs... ) info_ref! = Ref{Any}() - info_observer! = values_observer(; info=info_ref!) + info_observer! = values_observer(; info=(info_ref!)) observer! = compose_observers(observer!, info_observer!) state = alternating_update(operator, init; updater, observer!, kwargs...) return info_ref![].eigval, state diff --git a/src/mpo.jl b/src/mpo.jl index 87bac8d..fb87e65 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -1,7 +1,7 @@ using Adapt: adapt using BackendSelection: @Algorithm_str, Algorithm using DiagonalArrays: δ -using GradedUnitRanges: dag +using GradedArrays: dag using LinearAlgebra: dot # TODO: Add this back? # using ITensors: outer diff --git a/src/mps.jl b/src/mps.jl index a2d37bd..40b1d9d 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -1,5 +1,5 @@ using Adapt: adapt -using GradedUnitRanges: dual +using GradedArrays: dual using ITensors: hasqns using LinearAlgebra: qr using QuantumOperatorDefinitions: QuantumOperatorDefinitions, state diff --git a/src/solvers/alternating_update.jl b/src/solvers/alternating_update.jl index 330dfca..79c607f 100644 --- a/src/solvers/alternating_update.jl +++ b/src/solvers/alternating_update.jl @@ -105,7 +105,7 @@ function alternating_update( flush(stdout) end isdone = checkdone(; - state, sweep, outputlevel, observer=observer!, sweep_observer=sweep_observer! + state, sweep, outputlevel, observer=(observer!), sweep_observer=(sweep_observer!) ) isdone && break end diff --git a/src/solvers/dmrg_x.jl b/src/solvers/dmrg_x.jl index 0ed1f49..8d9972a 100644 --- a/src/solvers/dmrg_x.jl +++ b/src/solvers/dmrg_x.jl @@ -1,4 +1,4 @@ -using GradedUnitRanges: dag +using GradedArrays: dag using ITensors: uniqueind using LinearAlgebra: eigen using NamedDimsArrays: unname diff --git a/src/solvers/expand.jl b/src/solvers/expand.jl index 1fca947..85d5bc6 100644 --- a/src/solvers/expand.jl +++ b/src/solvers/expand.jl @@ -1,7 +1,7 @@ using Adapt: adapt using BackendSelection: @Algorithm_str, Algorithm using DiagonalArrays: δ -using GradedUnitRanges: dag +using GradedArrays: dag using ITensors: ITensors, Index, ITensor, commonind, denseblocks, directsum, hasqns, prime, uniqueinds using LinearAlgebra: normalize, svd, tr diff --git a/src/solvers/reducedconstantterm.jl b/src/solvers/reducedconstantterm.jl index c3e206f..9b8916b 100644 --- a/src/solvers/reducedconstantterm.jl +++ b/src/solvers/reducedconstantterm.jl @@ -1,4 +1,4 @@ -using GradedUnitRanges: dag +using GradedArrays: dag using ITensors: ITensors, ITensor, dim, prime """ diff --git a/src/solvers/sweep_update.jl b/src/solvers/sweep_update.jl index 3a3754c..f4d33b4 100644 --- a/src/solvers/sweep_update.jl +++ b/src/solvers/sweep_update.jl @@ -303,7 +303,7 @@ function region_update!( end set_nsite!(reduced_operator, nsite - 1) position!(reduced_operator, state, b1) - internal_kwargs = (; current_time, time_step=-time_step, outputlevel) + internal_kwargs = (; current_time, time_step=(-time_step), outputlevel) bond_reduced_state, info = updater( reduced_operator, bond_reduced_state; internal_kwargs, updater_kwargs... ) @@ -438,7 +438,7 @@ function region_update!( bond_reduced_state = state[b1] set_nsite!(reduced_operator, nsite - 1) position!(reduced_operator, state, b1) - internal_kwargs = (; current_time, time_step=-time_step, outputlevel) + internal_kwargs = (; current_time, time_step=(-time_step), outputlevel) bond_reduced_state, info = updater( reduced_operator, bond_reduced_state; internal_kwargs, updater_kwargs... ) diff --git a/src/solvers/tdvp.jl b/src/solvers/tdvp.jl index 41a1e48..a02e9b1 100644 --- a/src/solvers/tdvp.jl +++ b/src/solvers/tdvp.jl @@ -79,7 +79,7 @@ function tdvp( nsweeps=nothing, nsteps=nsweeps, (step_observer!)=default_sweep_observer(), - (sweep_observer!)=step_observer!, + (sweep_observer!)=(step_observer!), kwargs..., ) time_step, nsteps = time_step_and_nsteps(t, time_step, nsteps) diff --git a/test/base/test_mpo.jl b/test/base/test_mpo.jl index 047e738..9116a26 100644 --- a/test/base/test_mpo.jl +++ b/test/base/test_mpo.jl @@ -562,7 +562,7 @@ end ψ = orthogonalize(ψ0, 2) A = prod(ITensors.data(ψ)[2:(N - 1)]) randn!(A) - ψ[2:(N - 1), orthocenter=3] = A + ψ[2:(N - 1), orthocenter = 3] = A @test prod(ψ) ≈ ψ[1] * A * ψ[N] @test maxlinkdim(ψ) == 4 @test ITensorMPS.orthocenter(ψ) == 3 @@ -728,7 +728,7 @@ end s = siteinds(H2; plev=1) C = combiner.(s; tags="X") H2 .*= C - H2H2 = prime(H2; tags=!ts"X") * dag(H2) + H2H2 = prime(H2; tags=(!ts"X")) * dag(H2) @test @disable_warn_order prod(HH) ≈ prod(H2H2) end diff --git a/test/base/test_mps.jl b/test/base/test_mps.jl index 2acd673..f9d6c8b 100644 --- a/test/base/test_mps.jl +++ b/test/base/test_mps.jl @@ -1335,7 +1335,7 @@ end ψ = orthogonalize(ψ0, 2) A = prod(ITensorMPS.data(ψ)[2:(N - 1)]) randn!(A) - ψ[2:(N - 1), orthocenter=3] = A + ψ[2:(N - 1), orthocenter = 3] = A @test prod(ψ) ≈ ψ[1] * A * ψ[N] @test maxlinkdim(ψ) == 4 @test ITensorMPS.orthocenter(ψ) == 3 diff --git a/test/base/test_qnmpo.jl b/test/base/test_qnmpo.jl index 038f081..2e1eb78 100644 --- a/test/base/test_qnmpo.jl +++ b/test/base/test_qnmpo.jl @@ -337,6 +337,7 @@ test_combos = [(make_heisenberg_opsum, "S=1/2"), (make_hubbard_opsum, "Electron" @testset "QR/QL MPO tensors with complex block structures, H=$(test_combo[1])" for test_combo in test_combos + N, NNN = 10, 7 #10 lattice site, up 7th neight interactions sites = siteinds(test_combo[2], N; conserve_qns=true) H = test_combo[1](sites, NNN) diff --git a/test/base/test_solvers/test_tdvp.jl b/test/base/test_solvers/test_tdvp.jl index df6f4f3..b50c247 100644 --- a/test/base/test_solvers/test_tdvp.jl +++ b/test/base/test_solvers/test_tdvp.jl @@ -43,7 +43,7 @@ const elts = (Float32, Float64, Complex{Float32}, Complex{Float64}) total_time = elt(0.1) * im # matching total time and time step directions @test_throws "signs agree" tdvp(H, -total_time, ψ0; time_step) - @test_throws "signs agree" tdvp(H, total_time, ψ0; time_step=-time_step) + @test_throws "signs agree" tdvp(H, total_time, ψ0; time_step=(-time_step)) # total time is an integer of time steps @test_throws "integer" tdvp(H, total_time * 1.5, ψ0; time_step) end diff --git a/test/ext/ITensorMPSChainRulesCoreExt/test_optimization.jl b/test/ext/ITensorMPSChainRulesCoreExt/test_optimization.jl index 3835f95..b2965a4 100644 --- a/test/ext/ITensorMPSChainRulesCoreExt/test_optimization.jl +++ b/test/ext/ITensorMPSChainRulesCoreExt/test_optimization.jl @@ -124,7 +124,7 @@ include(joinpath(@__DIR__, "utils", "circuit.jl")) end @testset "State preparation (MPS)" begin - for gate in ["Ry"] #="Rx", =# + for gate in ["Ry"] nsites = 4 # Number of sites nlayers = 2 # Layers of gates in the ansatz gradtol = 1e-3 # Tolerance for stopping gradient descent From 35f04d858e3f9400015f77e458bd4094adf32e69 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 11 Apr 2025 16:02:06 -0400 Subject: [PATCH 31/33] Fix --- src/abstractmps.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/abstractmps.jl b/src/abstractmps.jl index e129a2d..be7ccb3 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -1925,7 +1925,7 @@ function setindex!( ## end # use the first link index if present, otherwise use the default tag - linktags = TagSet[ + linktags = [ (b=linkind(ψ, i); isnothing(b) ? defaultlinkindtags(i) : tags(b)) for i in firstsite:(lastsite - 1) ] @@ -1933,7 +1933,7 @@ function setindex!( ψA = MPST( A, sites; - leftinds=lind, + leftinds=linds, orthocenter=orthocenter - first(r) + 1, tags=linktags, kwargs..., From 5e79373343cb2f068195d66088811d782aa28ea3 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Mon, 14 Apr 2025 14:59:09 -0400 Subject: [PATCH 32/33] Format --- examples/dmrg/1d_heisenberg.jl | 21 +++++++++++++++++++-- examples/dmrg/Project.toml | 7 +++++++ src/dmrg.jl | 2 +- src/mps.jl | 5 +++-- 4 files changed, 30 insertions(+), 5 deletions(-) diff --git a/examples/dmrg/1d_heisenberg.jl b/examples/dmrg/1d_heisenberg.jl index 5cc6c15..8e7febe 100644 --- a/examples/dmrg/1d_heisenberg.jl +++ b/examples/dmrg/1d_heisenberg.jl @@ -2,6 +2,23 @@ using ITensors, ITensorMPS using Printf using Random +using Adapt: adapt +using JLArrays: JLArray +using Metal: MtlArray +using SerializedArrays: SerializedArray + +dev = adapt(SerializedArray) + +using ConstructionBase: constructorof +using GPUArraysCore: AbstractGPUArray +using TensorAlgebra: TensorAlgebra, factorize +function TensorAlgebra.factorize( + a::AbstractGPUArray, labels, codomain_labels, domain_labels; kwargs... +) + x, y = factorize(Array(a), labels, codomain_labels, domain_labels; kwargs...) + return constructorof(typeof(a)).((x, y)) +end + function heisenberg(N) os = OpSum() for j in 1:(N - 1) @@ -25,10 +42,10 @@ os = heisenberg(N) # Input operator terms which define a Hamiltonian # Convert these terms to an MPO tensor network -H = MPO(os, sites) +H = dev(MPO(os, sites)) # Create an initial random matrix product state -psi0 = random_mps(sites; maxdim=10) +psi0 = dev(random_mps(sites; maxdim=10)) # psi0 = random_mps(j -> isodd(j) ? "↑" : "↓", sites; maxdim=10) # psi0 = MPS(j -> isodd(j) ? "↑" : "↓", sites) diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index 0b09591..ec860b4 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -1,7 +1,14 @@ [deps] +ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" +ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" +JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4" +Metal = "dde4c033-4e86-420c-a63e-0dd931031962" NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SerializedArrays = "621c0da3-e96e-4f80-bd06-5ae31cdfcb39" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" diff --git a/src/dmrg.jl b/src/dmrg.jl index ad03190..61ecd5f 100644 --- a/src/dmrg.jl +++ b/src/dmrg.jl @@ -295,7 +295,7 @@ function dmrg( which_decomp, svd_alg, ) - maxtruncerr = max(maxtruncerr, spec.truncerr) + # maxtruncerr = max(maxtruncerr, spec.truncerr) ## TODO: Add back `@debug_check`, `checkflux`. ## @debug_check begin diff --git a/src/mps.jl b/src/mps.jl index 40b1d9d..ac7c7d1 100644 --- a/src/mps.jl +++ b/src/mps.jl @@ -238,7 +238,8 @@ function replacebond!( sbp1 = siteind(M, b + 1) indsMb = replaceind(indsMb, sb, sbp1) end - L, R, spec = factorize( + indsMb = indsMb ∩ inds(phi) + L, R = factorize( phi, indsMb; mindim, @@ -268,7 +269,7 @@ function replacebond!( "In replacebond!, got ortho = $ortho, only currently supports `left` and `right`." ) end - return spec + return nothing end """ From 98370a6f5c24e0be51b421f656afb0fe787d4914 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Tue, 13 May 2025 15:43:38 -0400 Subject: [PATCH 33/33] Update dependencies --- Project.toml | 4 ++-- examples/dmrg/1d_heisenberg_conserve_spin.jl | 2 +- examples/dmrg/Project.toml | 4 ++++ src/abstractmps.jl | 4 ++-- 4 files changed, 9 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 7cc2d5d..20c6724 100644 --- a/Project.toml +++ b/Project.toml @@ -47,9 +47,9 @@ ChainRulesCore = "1.10" Compat = "4.16.0" DiagonalArrays = "0.3" FillArrays = "1.13" -GradedArrays = "0.3.1" +GradedArrays = "0.4" HDF5 = "0.14, 0.15, 0.16, 0.17" -ITensors = "0.9" +ITensors = "0.10" IsApprox = "2.0" KrylovKit = "0.8.1, 0.9" LinearAlgebra = "1.10" diff --git a/examples/dmrg/1d_heisenberg_conserve_spin.jl b/examples/dmrg/1d_heisenberg_conserve_spin.jl index b854053..6b2ec08 100644 --- a/examples/dmrg/1d_heisenberg_conserve_spin.jl +++ b/examples/dmrg/1d_heisenberg_conserve_spin.jl @@ -1,8 +1,8 @@ using BlockSparseArrays +using GradedArrays using ITensors, ITensorMPS using Printf using Random -using SymmetrySectors function heisenberg(N) os = OpSum() diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index ec860b4..0ca1a4d 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -1,6 +1,9 @@ [deps] +Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +BlockSparseArrays = "2c9a651f-6452-4ace-a6ac-809f4280fbb4" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" +GradedArrays = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2" ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" @@ -12,3 +15,4 @@ NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SerializedArrays = "621c0da3-e96e-4f80-bd06-5ae31cdfcb39" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" +TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" diff --git a/src/abstractmps.jl b/src/abstractmps.jl index be7ccb3..bcbc091 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -1,5 +1,5 @@ using BackendSelection: @Algorithm_str, Algorithm -using GradedArrays: GradedUnitRanges, dag +using GradedArrays: GradedArrays, dag using IsApprox: Approx, IsApprox using ITensors: ITensors, @@ -725,7 +725,7 @@ function Base.map(f::Function, M::AbstractMPS; set_limits::Bool=true) end for (fname, fname!) in [ - (:(GradedUnitRanges.dag), :(dag!)), + (:(GradedArrays.dag), :(dag!)), (:(ITensors.prime), :(prime!)), (:(ITensors.setprime), :(setprime!)), (:(ITensors.noprime), :(noprime!)),