[WIP] Rewrite based on NamedDimsArrays.jl

Project.toml

@@ -1,20 +1,30 @@
 name = "ITensorMPS"
 uuid = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 authors = ["Matthew Fishman <[email protected]>", "Miles Stoudenmire <[email protected]>"]
-version = "0.3.18"
+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"
+GradedArrays = "bc96ca6e-b7c8-4bb6-888e-c93f838762c2"
 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"
+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"
+SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
+TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 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"
@@ -31,22 +41,32 @@ ITensorMPSPackageCompilerExt = "PackageCompiler"
 ITensorMPSZygoteRulesExt = ["ChainRulesCore", "ZygoteRules"]
 
 [compat]
-Adapt = "4.1.0"
+Adapt = "4.1"
+BackendSelection = "0.1"
 ChainRulesCore = "1.10"
 Compat = "4.16.0"
+DiagonalArrays = "0.3"
+FillArrays = "1.13"
+GradedArrays = "0.4"
 HDF5 = "0.14, 0.15, 0.16, 0.17"
-ITensors = "0.8, 0.9"
-IsApprox = "2.0.0"
+ITensors = "0.10"
+IsApprox = "2.0"
 KrylovKit = "0.8.1, 0.9"
 LinearAlgebra = "1.10"
-NDTensors = "0.4"
+NamedDimsArrays = "0.7"
 Observers = "0.2"
 PackageCompiler = "1, 2"
 Printf = "1.10"
+QuantumOperatorAlgebra = "0.1"
+QuantumOperatorDefinitions = "0.2"
 Random = "1.10"
-SerializedElementArrays = "0.1.0"
+SerializedElementArrays = "0.1"
+SparseArraysBase = "0.5"
+TensorAlgebra = "0.3"
 Test = "1.10"
-TupleTools = "1.5.0"
+TupleTools = "1.5"
+TypeParameterAccessors = "0.3"
+VectorInterface = "0.5"
 ZygoteRules = "0.2.2"
 julia = "1.10"
 

examples/dmrg/1d_heisenberg.jl

@@ -2,37 +2,62 @@ using ITensors, ITensorMPS
 using Printf
 using Random
 
-Random.seed!(1234)
+using Adapt: adapt
+using JLArrays: JLArray
+using Metal: MtlArray
+using SerializedArrays: SerializedArray
 
-let
-  N = 100
+dev = adapt(SerializedArray)
 
-  # 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)
+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
 
-  # 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 = 10
+
+# 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 = dev(MPO(os, sites))
+
+# Create an initial random matrix product state
+psi0 = dev(random_mps(sites; maxdim=10))
+# psi0 = random_mps(j -> isodd(j) ? "↑" : "↓", sites; maxdim=10)
+# psi0 = MPS(j -> isodd(j) ? "↑" : "↓", sites)
+
+# 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)

examples/dmrg/1d_heisenberg_conserve_spin.jl

@@ -1,39 +1,44 @@
+using BlockSparseArrays
+using GradedArrays
 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)
 
+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)
+
+# Create an initial random matrix product state
+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
+H = MPO(os, s)
+
+# 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)

examples/dmrg/Project.toml

@@ -1,6 +1,18 @@
 [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"
+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"
+TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"

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...

ext/ITensorMPSPackageCompilerExt/compile.jl

@@ -1,4 +1,4 @@
-using NDTensors: @Algorithm_str
+using BackendSelection: @Algorithm_str
 using ITensors: ITensors
 using PackageCompiler: PackageCompiler
 

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")
@@ -20,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")