[NDTensors] Fix bugs in SortedSets and SmallVectors

[mtfishman]

Fixes a bug in SortedSets when the inputs were not already sorted:

julia> using NDTensors.SortedSets

julia> SortedSet([3, 2, 1])
3-element SortedSet{Int64, Vector{Int64}, Base.Order.ForwardOrdering}
 1
 2
 3

I also plan to fix some other bugs related to SmallVectors.

@emstoudenmire could you share the bug you came across in union?

[codecov-commenter]

Codecov Report

All modified lines are covered by tests ✅

Comparison is base (15decbd) 85.40% compared to head (68ff015) 67.22%.

❗ Current head 68ff015 differs from pull request most recent head b2573d9. Consider uploading reports for the commit b2573d9 to get more accurate results

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@             Coverage Diff             @@
##             main    #1211       +/-   ##
===========================================
- Coverage   85.40%   67.22%   -18.18%     
===========================================
  Files          88       87        -1     
  Lines        8426     8388       -38     
===========================================
- Hits         7196     5639     -1557     
- Misses       1230     2749     +1519     

see 40 files with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[emstoudenmire]

Sure thing: here is one of the bugs (I think it's a bug – shouldn't the duplicate element get overwritten/ignored?)

julia> s = SmallSet{5}([2,4,6,8])
4-element SmallSet{5, Int64, Base.Order.ForwardOrdering}
 2
 4
 6
 8

julia> union(s,[10])
5-element SmallSet{5, Int64, Base.Order.ForwardOrdering}
 2
 4
 6
 8
 10

julia> union(s,[2])
ERROR: ArgumentError: New length 7 must be ≤ the maximum length 5.
Stacktrace:
 [1] resize!
   @ NDTensors.SmallVectors ~/.julia/dev/ITensors/NDTensors/src/SmallVectors/src/msmallvector/msmallvector.jl:71 [inlined]
 [2] unionsortedunique(itr1::NDTensors.SmallVectors.SmallVector{…}, itr2::Vector{…}, order::Base.Order.ForwardOrdering)
   @ NDTensors.SmallVectors ~/.julia/dev/ITensors/NDTensors/src/SmallVectors/src/BaseExt/sortedunique.jl:95
 [3] union
   @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:258 [inlined]
 [4] union(inds::SmallSet{5, Int64, Base.Order.ForwardOrdering}, items::Vector{Int64})
   @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:263
 [5] top-level scope
   @ REPL[4]:1
Some type information was truncated. Use `show(err)` to see complete types.

[emstoudenmire]

Here is a different one: in this case the set is not at its full capacity yet, and if I union in an element that's already there, it hits an issue involving trying to call setindex! on a static vector

julia> s = SmallSet{5}([2,4,6])
3-element SmallSet{5, Int64, Base.Order.ForwardOrdering}
 2
 4
 6

julia> union(s,[2])
ERROR: setindex!(::StaticArraysCore.SVector{5, Int64}, value, ::Int) is not defined.
 Hint: Use `MArray` or `SizedArray` to create a mutable static array
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] setindex!(a::StaticArraysCore.SVector{5, Int64}, value::Int64, i::Int64)
    @ StaticArrays ~/.julia/packages/StaticArrays/dTwvg/src/indexing.jl:3
  [3] setindex!
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SmallVectors/src/subsmallvector/subsmallvector.jl:57 [inlined]
  [4] reverse!
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SmallVectors/src/abstractsmallvector/deque.jl:98 [inlined]
  [5] circshift!
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SmallVectors/src/abstractsmallvector/deque.jl:119 [inlined]
  [6] deleteat!
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SmallVectors/src/abstractsmallvector/deque.jl:138 [inlined]
  [7] uniquesorted(vec::NDTensors.SmallVectors.SmallVector{5, Int64}, order::Base.Order.ForwardOrdering)
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/BaseExt/sorted.jl:47
  [8] NDTensors.SortedSets.SortedIndices{…}(a::Vector{…}, order::Base.Order.ForwardOrdering; issorted::typeof(issorted), allunique::Function)
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:26 [inlined]
  [9] SortedIndices
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:19 [inlined]
 [10] #_#6
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:42 [inlined]
 [11] SortedIndices
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:39 [inlined]
 [12] #SortedIndices#7
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:48 [inlined]
 [13] SortedIndices
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:45 [inlined]
 [14] union
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:259 [inlined]
 [15] union(inds::SmallSet{5, Int64, Base.Order.ForwardOrdering}, items::Vector{Int64})
    @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:263
 [16] top-level scope
    @ REPL[8]:1
Some type information was truncated. Use `show(err)` to see complete types.

[emstoudenmire]

Another one that might be related to the previous two. In this one I don't get an error, but and in a sense it works but the bug is that the length increases even though the number of elements does not:

julia> s = SmallSet{4}([1,2,3])
3-element SmallSet{4, Int64, Base.Order.ForwardOrdering}
 1
 2
 3

julia> union(s,[2])
4-element SmallSet{4, Int64, Base.Order.ForwardOrdering}
 1
 2
 3
 6101072064

[emstoudenmire]

One other question I have is, if the set has a non-trivial by function, will union use this function to determine equality for the purpose of overwriting (or ignoring, or throwing an error, whichever is the correct behavior) elements which are already there?

[mtfishman]

Ok thanks, I'll take a look at those union bugs.

Good question about union involving non-trivial by and elements that are equal according to by but not according to isequal. In Base Julia, this is the behavior for Set:

julia> struct Composite
         a::Symbol
         b::Float64
       end

julia> Base.isequal(x::Composite, y::Composite) = x.a == y.a

julia> Base.hash(x::Composite, h::UInt64) = hash(x.a, h)

julia> union(Set([Composite(:a, 2), Composite(:b, 3)]), [Composite(:a, 3)])
Set{Composite} with 2 elements:
  Composite(:a, 3.0)
  Composite(:b, 3.0)

so it replaces the element. Maybe we should follow that behavior.

[emstoudenmire]

Probably following what Julia does for sets is the best.

Do you think it would make sense to provide a keyword argument that specifies what to do for collisions? The default could be replace, but for Qn's I can think of cases where other behaviors would nicely implement some code patterns I need.

[mtfishman]

Good question, I was wondering about that. As we discussed, at least when fusing abelian QNs it is very closely related to union except it involves modular arithmetic of irreps as well. I think maybe that more sophisticated behavior should just be handled by another function, however. I would feel uncomfortable having union have different behavior for SortedSet as opposed to Base.Set (if you decided to try out Base.Set as a QN storage, suddenly your code wouldn't work anymore), and we can't tack on new behavior for union of Base.Set.

Maybe best to make up a new name for that more general functionality, maybe it could be called mapunion. For example it may have an interface like:

julia> mapunion(+, [1, 2], [2, 3]; default=0)
3-element Vector{Int64}:
 1
 4
 3

default defines what happens when one of the sets has a unique element. Then fusing abelian QNs would be:

mapunion(+, QN(("N", 1), ("Sz", 2)), QN(("N", 2)); default=Sector()) == QN(("N", 3), ("Sz", 2))

[mtfishman]

@emstoudenmire I believe I've fixed the issues you brought up, could you try out this branch and see if they are fixed?

[emstoudenmire]

Looks like all the cases I listed above are fixed for integers. Unfortunately I am running into an issue when I try to put Sector objects into SmallSet. I can give you the definition of Sector for testing if you want, but probably for these purposes it is just some struct holding an InlineStrings.String7 and an integer (technically it's a HalfInteger).

The error I get seems to be related to some code paths that is perhaps ending up calling an overly generic constructor of SortedSet, based on the error message and looking at the code. Ultimately this leads to calling collect on a Sector, which isn't defined:

julia> a1 = Sector("A",1)
Sector("A",1,SU{2})

julia> a2 = Sector("A",2)
Sector("A",2,SU{2})

julia> b = Sector("B",3)
Sector("B",3,SU{2})

julia> s = SmallSet{4}([a1,b];by=ITensorFusion.name)
2-element SmallSet{4, Sector, Base.Order.By{typeof(ITensorFusion.name), Base.Order.ForwardOrdering}}
 Sector("A",1,SU{2})
 Sector("B",3,SU{2})

julia> union(s,a2)
ERROR: MethodError: no method matching iterate(::Sector)

Closest candidates are:
  iterate(::ITensors.TruncEigen)
   @ ITensors ~/.julia/dev/ITensors/src/tensor_operations/matrix_decomposition.jl:196
  iterate(::ITensors.TruncEigen, ::Val{:V})
   @ ITensors ~/.julia/dev/ITensors/src/tensor_operations/matrix_decomposition.jl:197
  iterate(::ITensors.TruncEigen, ::Val{:spec})
   @ ITensors ~/.julia/dev/ITensors/src/tensor_operations/matrix_decomposition.jl:198
  ...

Stacktrace:
 [1] copyto!(dest::Vector{Any}, src::Sector)
   @ Base ./abstractarray.jl:938
 [2] _collect(cont::UnitRange{Int64}, itr::Sector, ::Base.HasEltype, isz::Base.HasLength)
   @ Base ./array.jl:763
 [3] collect(itr::Sector)
   @ Base ./array.jl:757
 [4] SortedSet{T, Data, Order}(a::Data, order::Order) where {T, Data<:(AbstractArray{T}), Order<:Base.Order.Ordering}
   @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:56 [inlined]
 [5] union(set::SmallSet{4, Sector, Base.Order.By{typeof(ITensorFusion.name), Base.Order.ForwardOrdering}}, items::Sector)
   @ NDTensors.SortedSets ~/.julia/dev/ITensors/NDTensors/src/SortedSets/src/sortedset.jl:277
 [6] top-level scope
   @ REPL[35]:1

[emstoudenmire]

Separately from that issue, I really like your idea of having a "map" variation of union, which would make it clearer that the user need to provide a function to handle collisions or reductions. (It's a lot better than my first idea, which could result in having an over-designed union function that takes too many keyword arguments.)

[mtfishman]

That's not really a SortedSet/SmallSet issue per se, just how generic code deals with treating objects as either iterable or not. A priori Julia doesn't know if it should treat Sector as a single object that you want to put into a set, or itself is a collection of more primitive objects. You would see the same issue using a Base.Set of Sector objects and union.

Two solutions are:

1. Call it as union(s, [a2]) or union(s, (a2,)), i.e. make it the only object in an iterable container like a Vector or Tuple.
2. Define iteration for Sector in a way such that it is treated as a trivial collection of itself, i.e. define Base.iterate(s::Sector, done=false) = done ? nothing : (s, true) and also Base.length(::Sector) = 1. Likely it makes sense to usually treat Sector objects as "scalar" objects, i.e. you wouldn't need to iterate over the data in a Sector object, just access them individually like ask for the name or irrep. That would make union(s, a2) work directly.

[mtfishman]

Additionally, something to consider to avoid having to pass by as a keyword argument in the SortedSet/SmallSet constructor is define comparison functions for Sector objects based on the name.

You could define things like:

Base.isless(s1::Sector, s2::Sector) = isless(name(s1), name(s2))
Base.:(==)(s1::Sector, s2::Sector) = name(s1) == name(s2)
Base.hash(s::Sector, h::UInt) = hash(name(s), hash(:Sector, h))

Then for all other purposes (using in a hash table like a Base.Dict or Base.Set, sorting, etc.) only the name will be used for comparison. I'm not sure if that's desirable to make those global definitions, but it would make it so you don't have to specify the by keyword argument in SortedSet. Note that this could be true for Sector, but then QN objects themselves could have other logic for comparison which make use of the sector names and irreps.

[emstoudenmire]

Ah ok, that was just my mistake. I had known you had to call union as union(s,[1]) for integers, but just forgot the brackets when testing for Sector entries. It's one of those things where it would be nice if the error message could be nicer, but those are extremely hard error messages to automatically determine (e.g. because like you said, how would Julia know if a Sector is meant to be a collection of things you are trying to pass to union?).

I'll just go with the solution of calling union(q,[Sector("A",1)]) .

Next I'll see how far I get with the existing set tools (union, intersect, setdiff) and let you know if something like mapunion is actually needed.

[emstoudenmire]

I just tried intersect with Sector data too and it works properly. I think all the known bugs in SortedSets are now fixed to the best of my knowledge, so this PR may be ready to merge.