BIG chunk of updates

- this commits a large number of updates to the BandGraphs code; many aspects have received corrections, very aggressive performance optimizations (e.g., work-arrays), as well as algorithmic improvements (e.g., checks of articulation point), and many changes to the types
- we also now do a much more general job of splitting degeneracies: in particular, we now allow splitting of degeneracies connected to non-nondegenerate nonmaximal irreps; this required multiset permutations.
- the main point of entry at the moment is in `test/scan-separable-irreps.jl`

BandGraphs/Project.toml

@@ -27,6 +27,7 @@ BandGraphsGraphMakieExt = "GraphMakie"
 
 [compat]
 GraphMakie = "0.5"
+JLD2 = "0.5.7"
 julia = "1.6"
 
 [extras]

BandGraphs/ext/BandGraphsGraphMakieExt.jl

@@ -218,15 +218,21 @@ function linearize_nonmaximal_y_coordinates(g_trail, ys::Vector)
         end
     end
     for is in values(track)
-        # a vertex assoc. w/ a nonmaximal irrep will only ever have two neighbors: it cannot
-        # have more than 2 irrep neighbors, as any connected maximal irrep must have equal
-        # or higher irrep dimensions; it cannot have less than 2 since it will always
+        # a vertex assoc. w/ a nonmaximal irrep will only ever have two neighbors (**): it
+        # cannot have more than 2 irrep neighbors, as any connected maximal irrep must have
+        # equal or higher irrep dimensions; it cannot have less than 2 since it will always
         # connect two different maximal irreps
+        # (**) This is true for a standard band graph, but not for a band graph that has
+        #      been "split" at a maximal vertex connected to a nonmaximal vertex of irrep
+        #      dimension higher than 1; ignoring this case, one could assert certain `only`
+        #      outcomes below, e.g., for `(in,out)neighbors(g_trail, i)`. To allow the
+        #      just-mentioned split-configuration, however, we just take averages in general
         avg_ys = Vector{Float64}(undef, length(is))
         for (idx_in_is, i) in enumerate(is)
-            in_i = only(inneighbors(g_trail, i))   # vertex index of ingoing maximal irrep
-            out_i = only(outneighbors(g_trail, i)) # vertex index of outgoing maximal irrep
-            avg_ys[idx_in_is] = (ys[in_i] + ys[out_i])/2
+            in_i = inneighbors(g_trail, i)   # vertex index of ingoing maximal irrep
+            out_i = outneighbors(g_trail, i) # vertex index of outgoing maximal irrep
+            avg_ys[idx_in_is] = (sum(@view ys[in_i])/length(in_i) + 
+                                 sum(@view ys[out_i])/length(out_i))/2
         end
         for (idx_in_is, i) in enumerate(is)
             y′ = avg_ys[idx_in_is]

BandGraphs/src/BandGraphs.jl

@@ -31,7 +31,12 @@ export
     permute_subgraphs,
     subduction_tables,
     plot_flattened_bandgraph,
+    Separability,
+    SeparabilityState,
+    is_separable,
+    is_feasible,
     findall_separable_vertices,
+    has_vertex,
     complete_split
 
 # ---------------------------------------------------------------------------------------- #
@@ -113,6 +118,8 @@ include("graphs.jl")
 include("subgraph-permutations.jl")
 include("graph_routing.jl")
 include("unfold.jl")
+include("multiset-permutation-indices.jl")
+include("count_connected_components.jl") # NB: remove when/if merged in Graphs.jl
 include("complete-split.jl")
 include("subsetsum.jl")
 include("separable.jl")

BandGraphs/src/count_connected_components.jl

@@ -0,0 +1,92 @@
+# Copy of https://github.com/JuliaGraphs/Graphs.jl/pull/407, until it is merged:
+# implements `count_connected_components(g, [label, search_queue])`.
+# TODO: remove when merged
+
+"""
+    _connected_components!(label, g, [search_queue])
+
+Fill `label` with the `id` of the connected component in the undirected graph
+`g` to which it belongs. Return a vector representing the component assigned
+to each vertex. The component value is the smallest vertex ID in the component.
+
+A `search_queue`, an empty `Vector{eltype(edgetype(g))}`, can be provided to reduce
+allocations if `_connected_components!` is intended to be called multiple times sequentially.
+If not provided, it is automatically instantiated.
+
+### Performance
+This algorithm is linear in the number of edges of the graph.
+"""
+function _connected_components!(
+    label::AbstractVector, g::AbstractGraph{T}, search_queue::Vector{T}=Vector{T}()
+) where {T}
+    isempty(search_queue) || error("provided `search_queue` is not empty")
+    for u in vertices(g)
+        label[u] != zero(T) && continue
+        label[u] = u
+        push!(search_queue, u)
+        while !isempty(search_queue)
+            src = popfirst!(search_queue)
+            for vertex in all_neighbors(g, src)
+                if label[vertex] == zero(T)
+                    push!(search_queue, vertex)
+                    label[vertex] = u
+                end
+            end
+        end
+    end
+    return label
+end
+
+"""
+    count_connected_components( g, [label, search_queue]; reset_label::Bool=false)
+
+Return the number of connected components in `g`.
+
+Equivalent to `length(connected_components(g))` but uses fewer allocations by not
+materializing the component vectors explicitly. Additionally, mutated work-arrays `label`
+and `search_queue` can be provided to reduce allocations further (see
+[`_connected_components!`](@ref)).
+
+## Keyword arguments
+- `reset_label :: Bool` (default, `false`): if `true`, `label` is reset to zero before
+  returning.
+
+## Example
+```
+julia> using Graphs
+
+julia> g = Graph(Edge.([1=>2, 2=>3, 3=>1, 4=>5, 5=>6, 6=>4, 7=>8]));
+
+length> connected_components(g)
+3-element Vector{Vector{Int64}}:
+ [1, 2, 3]
+ [4, 5, 6]
+ [7, 8]
+
+julia> count_connected_components(g)
+3
+```
+"""
+function count_connected_components(
+    g::AbstractGraph{T},
+    label::AbstractVector=zeros(T, nv(g)),
+    search_queue::Vector{T}=Vector{T}();
+    reset_label::Bool=false
+) where T
+    _connected_components!(label, g, search_queue)
+    c = count_unique(label)
+    reset_label && fill!(label, zero(eltype(label)))
+    return c
+end
+
+function count_unique(label::Vector{T}) where T
+    seen = Set{T}()
+    c = 0
+    for l in label
+        if l ∉ seen
+            push!(seen, l)
+            c += 1
+        end
+    end
+    return c
+end

BandGraphs/src/graphs.jl

@@ -3,7 +3,8 @@
 # While the graph associated with `A` in principle could be constructed from the `Graph(A)`
 # this loses a lot of relevant information. Instead, we construct it directly from
 # the subgraph and partition structures, and return additional metadata on vertices & edges
-function assemble_graph(subgraphs, partitions)
+function assemble_graph(subgraphs, partitions;
+            may_have_multiple_nonmax_to_max_edges::Bool=true)
     g = MetaGraph(Graph();
             label_type       = Tuple{String, Int},
             vertex_data_type = @NamedTuple{lgir::LGIrrep, maximal::Bool},
@@ -34,18 +35,35 @@ function assemble_graph(subgraphs, partitions)
             i′ = max_iridxs[i] # global column index in overall graph
             weight = a[i]
             add_edge!(g, label_for(g, i′), label_for(g, j′), (; weight=weight))
+            if may_have_multiple_nonmax_to_max_edges
+                # for a "standard" band graph, there will only ever be one nonzero element
+                # per column, since each nonmaximal vertex connects to exactly one maximal
+                # vertex; however, when we artificially split a maximal vertex, and this
+                # vertex is connected to a non-unit irdim nonmaximal node, this may no
+                # longer apply; so we allow this if `may_have_multiple_nonmax_to_max_edges`
+                # is true
+                while (i = findnext(≠(0), a, i+1); !isnothing(i))
+                    i′ = max_iridxs[something(i)]
+                    weight = a[something(i)]
+                    add_edge!(g, label_for(g, i′), label_for(g, j′), (; weight=weight))
+                end
+            end
         end
     end
     return g
 end
-assemble_graph(bandg::BandGraph) = assemble_graph(bandg.subgraphs, bandg.partitions)
+assemble_graph(bandg::BandGraph; kws...) = assemble_graph(bandg.subgraphs, bandg.partitions; kws...)
 
 
-function assemble_simple_graph(subgraphs, partitions)
+function assemble_simple_graph!(
+        g :: Graph,
+        subgraphs :: AbstractVector{<:SubGraph};
+        reset_edges::Bool = true,
+        may_have_multiple_nonmax_to_max_edges :: Bool = true)
     # does not have the right "weights" (and hence, not the right adjacency matrix either),
     # since a simple graph must have unit weights, but can still be used for connectivity 
     # analysis for instance
-    g = Graph(last(partitions[end].iridxs))
+    reset_edges && reset_edges!(g) # allows reuse of graph without reallocation
 
     # add edges to graph
     for subgraph in subgraphs
@@ -55,11 +73,28 @@ function assemble_simple_graph(subgraphs, partitions)
             i = something(findfirst(≠(0), a)) # local row index in block/subgraph
             i′ = max_iridxs[i] # global column index in overall graph
             add_edge!(g, i′, j′)
+            if may_have_multiple_nonmax_to_max_edges
+                while (i = findnext(≠(0), a, i+1); !isnothing(i))
+                    i′ = max_iridxs[something(i)]
+                    add_edge!(g, i′, j′)
+                end
+            end
         end
     end
     return g
 end
-assemble_simple_graph(bandg::BandGraph) = assemble_simple_graph(bandg.subgraphs, bandg.partitions)
+function assemble_simple_graph!(g::Graph, bandg::BandGraph; kws...)
+    assemble_simple_graph!(g, bandg.subgraphs; kws...)
+end
+reset_edges!(g::Graph) = (foreach(l->resize!(l, 0), g.fadjlist); g.ne = 0; g)
+function assemble_simple_graph(subgraphs :: AbstractVector{<:SubGraph}, partitions; kws...)
+    g = Graph(last(partitions[end].iridxs))
+    assemble_simple_graph!(g, subgraphs; reset_edges=false, kws...)
+end
+function assemble_simple_graph(bandg::BandGraph; kws...)
+    g = Graph(nv(bandg))
+    assemble_simple_graph!(g, bandg.subgraphs; reset_edges=false, kws...)
+end
 
 # ---------------------------------------------------------------------------------------- #
 

BandGraphs/src/multiset-permutation-indices.jl

@@ -0,0 +1,83 @@
+struct MultiSetPermutationIndices
+    f::Vector{Int}
+end
+Base.eltype(::Type{MultiSetPermutationIndices}) = Vector{Int}
+Base.length(c::MultiSetPermutationIndices) = _multinomial(c.f)
+
+# copy of https://github.com/JuliaMath/Combinatorics.jl/pull/170, until merged & public
+function _multinomial(k)
+    s = zero(eltype(k))
+    result = one(eltype(k))
+    @inbounds for i in k
+        s += i
+        result *= binomial(s, i)
+    end
+    result
+end
+
+"""
+    multiset_permutation_indices(f :: Vector{Int})
+
+Return a generator over all index permutations corresponding to the multiset permutations of
+a multiset `f = {f[1]*1, f[2]*2, ...}`, i.e. of a multiset with `f[1]` elements of `1`,
+`f[2]` elements of `2`, and so on.
+
+More explicitly, `f[i]` indicates the frequency of the `i`th element of the underlying
+multiset `a`, which can consequently be considered as the sequence or sorted multiset:
+`a` = [`1` repeated `f[1]` times, `2` repeated `f[2]` times, ...].
+
+The indices `inds` of the unique permutations of `a` are returned. Compared to the iterator
+`multiset_permutations(a, length(a))` - which returns the actual permutations, not the
+indices - from Combinatorics.jl, this means that
+`a[ind] == collect(multiset_permutations(a), length(a))[ind]` for all `ind in inds`.
+
+## Implementation
+The implementation here is adapted from the implementation in `multiset_permutations` in
+Combinatorics.jl, which in turn appears to follow the Algorithm L (p. 319) from Knuth's
+TAOCP Vol. 4A.
+"""
+multiset_permutation_indices(f::Vector{Int}) = MultiSetPermutationIndices(f)
+
+function canonical_multiset(f::Vector{Int})
+    # build the sequence `a` = [`1` repeated `f[1]` times, `2` repeated `f[2]` times, ...].
+    a = Vector{Int}(undef, sum(f))
+    N = 1
+    for (i, n) in enumerate(f)
+        @inbounds a[N:N+n-1] .= i
+        N += n
+    end
+    return a
+end
+
+function Base.iterate(
+        p::MultiSetPermutationIndices, 
+        composite_state = (a=canonical_multiset(p.f); (a, collect(eachindex(a))))
+        )
+    state, indices = composite_state
+    if !isempty(state) && state[1] > length(state)
+        return nothing
+    end
+    return nextpermutation!(state, indices)
+end
+
+function nextpermutation!(state::Vector{Int}, indices::Vector{Int})
+    i = length(state) - 1
+    indices_next = copy(indices)
+    while i ≥ 1 && state[i] ≥ state[i+1]
+        i -= 1
+    end
+    if i > 0
+        j = length(state)
+        while j > i && state[i] ≥ state[j]
+            j -= 1
+        end
+        state[i], state[j] = state[j], state[i]
+        indices_next[i], indices_next[j] = indices_next[j], indices_next[i]
+        reverse!(state, i+1)
+        reverse!(indices_next, i+1)
+    else
+        state[1] = length(state) + 1 # finished; no more permutations
+    end
+
+    return indices, (state, indices_next)
+end