Implement StructuredVector for storing scaling dimensions (#185)

* Implement `StructuredVector` for storing scaling dimensions

* fix include statement

* fix bug

* remove Adwait AI slop

I'm losing the battle against entropy

* small changes

* fix typo

Author: VictorVanthilt

docs/src/cft.md

@@ -26,7 +26,7 @@ The `CFTData` struct has two fields:
 - `scaling_dimensions`
 
 The `central_charge` can be either `missing` (when using the $[1, 1, 0]$ shape), or a number.
-The `scaling_dimensions` field is a `SectorVector` from TensorKit.jl.
+The `scaling_dimensions` field is a `StructuredVector`.
 
 The `scaling_dimensions` can be indexed like an `AbstractVector` (i.e. with scalars, slices, ...), or with sectors (e.g. `Z2Irrep(0)`), which will provide the scaling dimensions associated with that sector/charge.
 To check which sectors you can index the `scaling_dimensions` with you can use `keys(scaling_dimensions)`.

src/TNRKit.jl

@@ -52,6 +52,7 @@ include("schemes/correlationhotrg.jl")
 
 #Thermal TNR
 include("schemes/ttnr.jl")
+
 # Loop Methods
 include("schemes/looptnr.jl")
 include("schemes/symmetric_looptnr.jl")
@@ -124,6 +125,8 @@ export phi4_complex, phi4_complex_impϕ, phi4_complex_impϕdag, phi4_complex_imp
 include("utility/free_energy.jl")
 export free_energy
 
+include("utility/structuredvector.jl")
+
 include("utility/cft.jl")
 export CFTData, central_charge
 
@@ -145,4 +148,5 @@ include("utility/blocking.jl")
 export block_tensors
 
 include("utility/network_value.jl")
+
 end

src/models/clock.jl

@@ -148,5 +148,3 @@ function ZN_gauge_theory_dual(N::Int, β::Real; T::Type{<:Number} = ComplexF64)
             symmetric_space ⊗ symmetric_space ⊗ symmetric_space_dual ⊗ symmetric_space_dual,
     )
 end
-
-ZN_gaugetheory_dual(N::Int, β::Real; kwargs...) = ZN_gauge_theory_dual(N, β; kwargs...)

src/models/ising.jl

@@ -55,7 +55,7 @@ function ising_3D_free_energy_htse(β::Real; J::Real = 1.0, max_order::Int = 24)
     f -= series / β
     return f
 end
-ising_3d_free_energy_htse(; kwargs...) = ising_3D_free_energy_htse(ising_βc_3D; kwargs...)
+ising_3D_free_energy_htse(; kwargs...) = ising_3D_free_energy_htse(ising_βc_3D; kwargs...)
 
 """
     ising_bond_tensor(β::Real, T::Type{<:Number})

src/utility/cft.jl

@@ -1,44 +1,3 @@
-# Extensions on top of TensorKit.SectorVector
-Base.broadcasted(f, v::TensorKit.SectorVector) = TensorKit.SectorVector(broadcast(f, parent(v)), v.structure)
-Base.broadcasted(f, v::TensorKit.SectorVector, a) = TensorKit.SectorVector(broadcast(f, parent(v), a), v.structure)
-Base.broadcasted(f, a, v::TensorKit.SectorVector) = TensorKit.SectorVector(broadcast(f, a, parent(v)), v.structure)
-function Base.broadcasted(f, v1::TensorKit.SectorVector, v2::TensorKit.SectorVector)
-    if v1.structure != v2.structure
-        throw(ArgumentError("Cannot broadcast two SectorVectors with different structures"))
-    end
-    return TensorKit.SectorVector(broadcast(f, parent(v1), parent(v2)), v1.structure)
-end
-
-function Base.filter(f, v::TensorKit.SectorVector)
-    data = copy(parent(v))
-    structure = copy(v.structure)
-
-    kept_inds = findall(f, parent(v))
-    sectors = keys(structure)
-    for (i, sector) in enumerate(sectors)
-        structure[sector] = i == 1 ? (1:findlast(x -> x <= structure[sector].stop, kept_inds)) : ((structure[sectors[i - 1]].stop + 1):findlast(x -> x <= structure[sector].stop, kept_inds))
-    end
-    data = data[kept_inds]
-    return TensorKit.SectorVector(data, structure)
-end
-
-function Base.sort(v::TensorKit.SectorVector; kwargs...)
-    # sort within the sectors, then concatenate the data and update the structure
-    # Ideally this would sort the total data, but sectorvectors are only compatible with
-    # structures that contain unitranges, we cannot interweave the data of different sectors.
-    data = copy(parent(v))
-    newdata = similar(data)
-    structure = copy(v.structure)
-    sectors = keys(structure)
-    for sector in sectors
-        newdata[structure[sector]] = sort(data[structure[sector]]; kwargs...)
-    end
-    return TensorKit.SectorVector(newdata, structure)
-end
-
-Base.:*(a::Number, v::TensorKit.SectorVector) = scale(v, a)
-Base.:*(v::TensorKit.SectorVector, a::Number) = scale(v, a)
-
 """
     CFTData{E, I} where {E, I}
 
@@ -51,15 +10,15 @@ A struct to hold conformal data extracted from a TNR scheme.
 
 # Fields
     - `central_charge::Union{E, Missing}`: The central charge of the CFT. Will be `nothing` if not calculated.
-    - `scaling_dimensions::TensorKit.SectorVector{E, I}`: The scaling dimensions of the CFT, organized in a `TensorKit.SectorVector` where the sectors correspond to different spin sectors (or other quantum numbers) and the data contains the scaling dimensions within those sectors
+    - `scaling_dimensions::StructuredVector{E, K, V, A}`: The scaling dimensions of the CFT, organized in a `StructuredVector` where the sectors correspond to different spin sectors (or other quantum numbers) and the data contains the scaling dimensions within those sectors
 
 """
-struct CFTData{E, I}
+struct CFTData{E, K, V, A <: AbstractVector{E}}
     "Central charge of the CFT. Will be `nothing` if not calculated."
     central_charge::Union{E, Missing}
 
     "Scaling dimensions of the CFT."
-    scaling_dimensions::TensorKit.SectorVector{E, I}
+    scaling_dimensions::StructuredVector{E, K, V, A}
 end
 
 function Base.show(io::IO, data::CFTData)
@@ -71,7 +30,8 @@ end
 
 function CFTData(T::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], kwargs...) where {E, S}
     if shape == [1, 1, 0] # trivial implementation
-        return CFTData(missing, _scaling_dimensions(T))
+        scaling_dimensions = _scaling_dimensions(T)
+        return CFTData(missing, StructuredVector(scaling_dimensions, Dict([Trivial => collect(eachindex(scaling_dimensions))])))
     else
         CFTData(T, T; shape, kwargs...)
     end
@@ -199,23 +159,23 @@ function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
     norm_const_0 = spec_sector[one(I)][1]
     central_charge = 6 / pi / (Imτ - area / 4) * log(norm_const_0)
 
-    # Construct a SectorVector from the data of the different sectors
+    # Construct a StructuredVector from the data of the different sectors
     data = ComplexF64[]
-    structure = TensorKit.SectorDict{sectortype(xspace), UnitRange{Int}}()
+    structure = Dict{eltype(sectors(fuse(xspace))), Vector{Int}}()
     last_index = 1
     for charge in sectors(fuse(xspace))
         DeltaS = -1 / (2 * pi * Imτ) * log.(spec_sector[charge] / norm_const_0)
         if !(relative_shift ≈ 0)
             push!(data, (real.(DeltaS) + imag.(DeltaS) / relative_shift * im)...)
-            structure[charge] = last_index:(last_index + length(DeltaS) - 1)
+            structure[charge] = [last_index:(last_index + length(DeltaS) - 1)...]
         else
             push!(data, real.(DeltaS)...)
-            structure[charge] = last_index:(last_index + length(DeltaS) - 1)
+            structure[charge] = [last_index:(last_index + length(DeltaS) - 1)...]
         end
         last_index += length(DeltaS)
     end
 
-    sv = TensorKit.SectorVector(data, structure)
+    sv = StructuredVector(data, structure)
     sv = sort(sv; by = real)
     sv = filter(x -> real(x) ≤ 1.0e16, sv)
 

src/utility/structuredvector.jl

@@ -0,0 +1,45 @@
+struct StructuredVector{E, K, V, A <: AbstractVector{E}} <: AbstractVector{E}
+    data::A
+    structure::Dict{K, V}
+end
+
+@inline Base.getindex(v::StructuredVector, i::Int) = getindex(parent(v), i)
+@inline Base.getindex(v::StructuredVector{E, K}, keys::K) where {E, K} = parent(v)[v.structure[keys]]
+@inline Base.setindex!(v::StructuredVector, val, i::Int) = setindex!(parent(v), val, i)
+
+Base.size(v::StructuredVector, args...) = size(parent(v), args...)
+Base.size(v::StructuredVector) = size(parent(v))
+Base.copy(v::StructuredVector) = StructuredVector(copy(v.data), v.structure)
+Base.parent(v::StructuredVector) = v.data
+
+function Base.sort(v::StructuredVector; kwargs...)
+    p = sortperm(v.data; kwargs...)
+    inv_p = invperm(p)
+    newdict = Dict(k => sort(inv_p[v.structure[k]]) for k in keys(v.structure))
+    return StructuredVector(v.data[p], newdict)
+end
+
+function Base.filter(f, v::StructuredVector)
+    kept_inds = findall(f, parent(v))
+    data = parent(v)[kept_inds]
+    old_to_new = Dict(old_ind => new_ind for (new_ind, old_ind) in enumerate(kept_inds))
+
+    new_structure = Dict{keytype(v.structure), Vector{Int}}()
+    for (sector, inds) in v.structure
+        new_structure[sector] = [old_to_new[ind] for ind in inds if haskey(old_to_new, ind)]
+    end
+
+    return StructuredVector(data, new_structure)
+end
+
+function Base.show(io::IO, v::StructuredVector)
+    println(io, "StructuredVector with keys: ", keys(v.structure))
+    return print(io, "Data: ", v.data)
+end
+
+Base.:*(v::StructuredVector, x::Number) = StructuredVector(v.data .* x, v.structure)
+Base.:*(x::Number, v::StructuredVector) = StructuredVector(x .* v.data, v.structure)
+Base.:/(v::StructuredVector, x::Number) = StructuredVector(v.data ./ x, v.structure)
+Base.:/(x::Number, v::StructuredVector) = StructuredVector(x ./ v.data, v.structure)
+
+Base.keys(v::StructuredVector) = keys(v.structure)

test/schemes/schemes.jl

@@ -9,7 +9,7 @@ println("---------------------")
 
 T = classical_ising()
 T_3D = classical_ising_3D()
-const f_benchmark3D = ising_3D_free_energy_htse(ising_βc_3D)
+const f_benchmark3D = ising_3D_free_energy_htse()
 
 function cft_finalize!(scheme)
     finalize!(scheme)