diff --git a/Project.toml b/Project.toml index f6be3a0..09603fd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,12 +1,13 @@ name = "ITensorBase" uuid = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" authors = ["ITensor developers and contributors"] -version = "0.1.2" +version = "0.1.3" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" DerivableInterfaces = "6c5e35bf-e59e-4898-b73c-732dcc4ba65f" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MapBroadcast = "ebd9b9da-f48d-417c-9660-449667d60261" NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde" UnallocatedArrays = "43c9e47c-e622-40fb-bf18-a09fc8c466b6" @@ -16,6 +17,7 @@ UnspecifiedTypes = "42b3faec-625b-4613-8ddc-352bf9672b8d" Accessors = "0.1.39" DerivableInterfaces = "0.3.7" FillArrays = "1.13.0" +LinearAlgebra = "1.10" MapBroadcast = "0.1.5" NamedDimsArrays = "0.3.0" UnallocatedArrays = "0.1.1" diff --git a/README.md b/README.md index a25436a..5216bc7 100644 --- a/README.md +++ b/README.md @@ -19,13 +19,8 @@ julia> Pkg.add("ITensorBase") ````julia using ITensorBase: ITensorBase, ITensor, Index -```` - -TODO: This should be `TensorAlgebra.qr`. - -````julia using LinearAlgebra: qr -using NamedDimsArrays: NamedDimsArray, aligndims, dimnames, name, unname +using NamedDimsArrays: aligndims, unname using Test: @test i = Index(2) j = Index(2) diff --git a/examples/README.jl b/examples/README.jl index 48de04f..70c1717 100644 --- a/examples/README.jl +++ b/examples/README.jl @@ -20,9 +20,8 @@ julia> Pkg.add("ITensorBase") # ## Examples using ITensorBase: ITensorBase, ITensor, Index -# TODO: This should be `TensorAlgebra.qr`. using LinearAlgebra: qr -using NamedDimsArrays: NamedDimsArray, aligndims, dimnames, name, unname +using NamedDimsArrays: aligndims, unname using Test: @test i = Index(2) j = Index(2) diff --git a/src/ITensorBase.jl b/src/ITensorBase.jl index f5b931b..ccecd09 100644 --- a/src/ITensorBase.jl +++ b/src/ITensorBase.jl @@ -1,5 +1,8 @@ module ITensorBase +export ITensor, Index + +using Accessors: @set using MapBroadcast: Mapped using NamedDimsArrays: NamedDimsArrays, @@ -14,16 +17,53 @@ using NamedDimsArrays: name, named, nameddimsindices, + setname, + setnameddimsindices, unname +const Tag = String +const TagSet = Set{Tag} + +tagset(tags::String) = Set(filter(!isempty, String.(strip.(split(tags, ","))))) +tagset(tags::TagSet) = tags + +function tagsstring(tags::TagSet) + str = "" + length(tags) == 0 && return str + tags_vec = collect(tags) + for n in 1:(length(tags_vec) - 1) + str *= "$(tags_vec[n])," + end + str *= "$(tags_vec[end])" + return str +end + @kwdef struct IndexName <: AbstractName id::UInt64 = rand(UInt64) + tags::TagSet = TagSet() plev::Int = 0 - tags::Set{String} = Set{String}() - namedtags::Dict{Symbol,String} = Dict{Symbol,String}() end NamedDimsArrays.randname(n::IndexName) = IndexName() +id(n::IndexName) = n.id +tags(n::IndexName) = n.tags +plev(n::IndexName) = n.plev + +settags(n::IndexName, tags) = @set n.tags = tags +addtags(n::IndexName, ts) = settags(n, tags(n) ∪ tagset(ts)) + +setprime(n::IndexName, plev) = @set n.plev = plev +prime(n::IndexName) = setprime(n, plev(n) + 1) + +function Base.show(io::IO, i::IndexName) + idstr = "id=$(id(i) % 1000)" + tagsstr = !isempty(tags(i)) ? "|\"$(tagsstring(tags(i)))\"" : "" + primestr = primestring(plev(i)) + str = "IndexName($(idstr)$(tagsstr))$(primestr)" + print(io, str) + return nothing +end + struct IndexVal{Value<:Integer} <: AbstractNamedInteger{Value,IndexName} value::Value name::IndexName @@ -41,7 +81,22 @@ struct Index{T,Value<:AbstractUnitRange{T}} <: AbstractNamedUnitRange{T,Value,In name::IndexName end -Index(length::Int) = Index(Base.OneTo(length), IndexName()) +function Index(length::Int; tags=TagSet(), kwargs...) + return Index(Base.OneTo(length), IndexName(; tags=tagset(tags), kwargs...)) +end +function Index(length::Int, tags::String; kwargs...) + return Index(Base.OneTo(length), IndexName(; kwargs..., tags=tagset(tags))) +end + +# TODO: Define for `NamedDimsArrays.NamedViewIndex`. +id(i::Index) = id(name(i)) +tags(i::Index) = tags(name(i)) +plev(i::Index) = plev(name(i)) + +# TODO: Define for `NamedDimsArrays.NamedViewIndex`. +addtags(i::Index, tags) = setname(i, addtags(name(i), tags)) +prime(i::Index) = setname(i, prime(name(i))) +Base.adjoint(i::Index) = prime(i) # Interface # TODO: Overload `Base.parent` instead. @@ -51,6 +106,29 @@ NamedDimsArrays.name(i::Index) = i.name # Constructor NamedDimsArrays.named(i::AbstractUnitRange, name::IndexName) = Index(i, name) +function primestring(plev) + if plev < 0 + return " (warning: prime level $plev is less than 0)" + end + if plev == 0 + return "" + elseif plev > 3 + return "'$plev" + else + return "'"^plev + end +end + +function Base.show(io::IO, i::Index) + lenstr = "length=$(dename(length(i)))" + idstr = "|id=$(id(i) % 1000)" + tagsstr = !isempty(tags(i)) ? "|\"$(tagsstring(tags(i)))\"" : "" + primestr = primestring(plev(i)) + str = "Index($(lenstr)$(idstr)$(tagsstr))$(primestr)" + print(io, str) + return nothing +end + struct NoncontiguousIndex{T,Value<:AbstractVector{T}} <: AbstractNamedVector{T,Value,IndexName} value::Value @@ -103,17 +181,30 @@ struct AllocatableArrayInterface <: AbstractAllocatableArrayInterface end unallocatable(a::AbstractITensor) = NamedDimsArray(a) -@interface ::AbstractAllocatableArrayInterface function Base.setindex!( - a::AbstractArray, value, I::Int... -) +function setindex_allocatable!(a::AbstractArray, value, I...) allocate!(specify_eltype!(a, typeof(value))) # TODO: Maybe use `@interface interface(a) a[I...] = value`? unallocatable(a)[I...] = value return a end +# TODO: Combine these by using `Base.to_indices`. +@interface ::AbstractAllocatableArrayInterface function Base.setindex!( + a::AbstractArray, value, I::Int... +) + setindex_allocatable!(a, value, I...) + return a +end +@interface ::AbstractAllocatableArrayInterface function Base.setindex!( + a::AbstractArray, value, I::AbstractNamedInteger... +) + setindex_allocatable!(a, value, I...) + return a +end + @derive AllocatableArrayInterface() (T=AbstractITensor,) begin Base.setindex!(::T, ::Any, ::Int...) + Base.setindex!(::T, ::Any, ::AbstractNamedInteger...) end mutable struct ITensor <: AbstractITensor @@ -127,9 +218,42 @@ using Accessors: @set setdenamed(a::ITensor, denamed) = (@set a.parent = denamed) setdenamed!(a::ITensor, denamed) = (a.parent = denamed) +function ITensor(elt::Type, I1::Index, I_rest::Index...) + I = (I1, I_rest...) + # TODO: Use `FillArrays.Zeros`. + return ITensor(zeros(elt, length.(dename.(I))...), I) +end + function ITensor(I1::Index, I_rest::Index...) I = (I1, I_rest...) return ITensor(Zeros{UnspecifiedZero}(length.(dename.(I))...), I) end +function ITensor() + return ITensor(Zeros{UnspecifiedZero}(), ()) +end + +inds(a::AbstractITensor) = nameddimsindices(a) +setinds(a::AbstractITensor, inds) = setnameddimsindices(a, inds) + +function uniqueinds(a1::AbstractITensor, a_rest::AbstractITensor...) + return setdiff(inds(a1), inds.(a_rest)...) +end +function uniqueind(a1::AbstractITensor, a_rest::AbstractITensor...) + return only(uniqueinds(a1, a_rest...)) +end + +function commoninds(a1::AbstractITensor, a_rest::AbstractITensor...) + return intersect(inds(a1), inds.(a_rest)...) +end +function commonind(a1::AbstractITensor, a_rest::AbstractITensor...) + return only(commoninds(a1, a_rest...)) +end + +# TODO: Use `replaceinds`/`mapinds`, based on +# `replacenameddimsindices`/`mapnameddimsindices`. +prime(a::AbstractITensor) = setinds(a, prime.(inds(a))) + +include("quirks.jl") + end diff --git a/src/quirks.jl b/src/quirks.jl new file mode 100644 index 0000000..f9d4e5c --- /dev/null +++ b/src/quirks.jl @@ -0,0 +1,31 @@ +# TODO: Define this properly. +dag(i::Index) = i +# TODO: Deprecate. +dim(i::Index) = dename(length(i)) +# TODO: Define this properly. +hasqns(i::Index) = false +# TODO: Deprecate. +itensor(parent::AbstractArray, nameddimsindices) = ITensor(parent, nameddimsindices) +function itensor(parent::AbstractArray, i1::Index, i_rest::Index...) + return ITensor(parent, (i1, i_rest...)) +end + +# This seems to be needed to get broadcasting working. +# TODO: Investigate this and see if we can get rid of it. +Base.Broadcast.extrude(a::AbstractITensor) = a + +# TODO: Generalize this. +# Maybe define it as `oneelement`, and base it on +# `FillArrays.OneElement` (https://juliaarrays.github.io/FillArrays.jl/stable/#FillArrays.OneElement). +function onehot(iv::Pair{<:Index,<:Int}) + a = ITensor(first(iv)) + a[last(iv)] = one(Bool) + return a +end + +using LinearAlgebra: svd +# TODO: Define this in `MatrixAlgebra.jl`/`TensorAlgebra.jl`. +function factorize(a::AbstractITensor, args...; kwargs...) + U, S, V = svd(a, args...; kwargs...) + return U, S * V +end