Add gaussian_elimination

Liozou · Jul 4, 2022 · 2c54e06 · 2c54e06
1 parent 072a612
commit 2c54e06
Show file tree

Hide file tree

Showing 3 changed files with 110 additions and 73 deletions.
diff --git a/docs/src/rings.md b/docs/src/rings.md
@@ -170,6 +170,7 @@ PeriodicGraphs.symdiff_cycles!
 PeriodicGraphs.symdiff_cycles
 PeriodicGraphs.IterativeGaussianElimination
 PeriodicGraphs.gaussian_elimination!
+PeriodicGraphs.gaussian_elimination
 PeriodicGraphs.intersect_cycles!
 PeriodicGraphs.intersect_cycles
 PeriodicGraphs.union_cycles!

diff --git a/src/algorithms/rings.jl b/src/algorithms/rings.jl
@@ -1187,26 +1187,7 @@ end
 IterativeGaussianElimination(ring, sizehint=ring[1]) = IterativeGaussianEliminationNone(ring, sizehint)
 IterativeGaussianElimination() = IterativeGaussianEliminationNone()
 
-# For an IterativeGaussianElimination of i rings of size at most l, symdiff_cycles cost at
-# most (l + lenr) + (2l + lenr) + ... + (il + lenr) = O(i^2*l +i*lenr)
-# If the size of the k-th ring is at most k*l, then the worst-case complexity is O(i^3*l)
-# Calling gaussian_elimination ν times thus leads to a worst-case complexity in O(ν^4*l)
-# The reasonning can probably be refined...
-"""
-    gaussian_elimination!(gauss::IterativeGaussianElimination, r::Vector{Int}) where T
-
-Test whether `r` can be expressed as a sum of vectors stored in `gauss`, and store `r` if
-not. "sum" refers to the symmetric difference of boolean vectors, represented in sparse
-format as the ordered list of non-zero indices.
-
-If `gauss isa IterativeGaussianEliminationLength`, return whether `r` can be expressed as a
-sum of strictly smaller vectors.
-
-Otherwise, return `true` when `r` is a sum of any previously encoutered vectors.
-If `gauss` isa `IterativeGaussianEliminationDecomposition`, query `retrieve_track(gauss)`
-to obtain the sorted list of indices of such previously encountered vectors.
-"""
-function gaussian_elimination!(gauss::IterativeGaussianElimination{T}, r::Vector{Int}) where T
+function _gaussian_elimination(gauss::IterativeGaussianElimination{T}, r::Vector{Int}) where T
     rings = gauss.rings
     shortcuts = gauss.shortcuts
     buffer1::Vector{Int} = gauss.buffer1
@@ -1220,6 +1201,7 @@ function gaussian_elimination!(gauss::IterativeGaussianElimination{T}, r::Vector
         empty!(track)
     end
     r1 = r[1]
+    maxlen = 0
 
     idx::Int32 = r1 > lenshort ? zero(Int32) : shortcuts[r1]
     if !iszero(idx)
@@ -1230,7 +1212,7 @@ function gaussian_elimination!(gauss::IterativeGaussianElimination{T}, r::Vector
             push!(track, idx)
         end
         symdiff_cycles!(buffer1, r, ridx)
-        isempty(buffer1) && @goto notindependentreturn
+        isempty(buffer1) && return true, r1, maxlen, buffer1
         r1 = buffer1[1]
         idx = r1 > lenshort ? zero(Int32) : shortcuts[r1]
     else
@@ -1244,30 +1226,76 @@ function gaussian_elimination!(gauss::IterativeGaussianElimination{T}, r::Vector
             push!(track, idx)
         end
         symdiff_cycles!(buffer2, buffer1, ridx)
-        isempty(buffer2) && @goto notindependentreturn
+        isempty(buffer2) && return true, r1, maxlen, buffer1
         r1 = buffer2[1]
         idx = r1 > lenshort ? zero(Int32) : shortcuts[r1]
         buffer2, buffer1 = buffer1, buffer2
     end
 
+    return false, r1, maxlen, buffer1
+end
+
+"""
+    gaussian_elimination(gauss::IterativeGaussianElimination, r::Vector{Int})
+
+Test whether `r` can be expressed as a sum of vectors stored in `gauss`.
+
+See [`PeriodicGraphs.gaussian_elimination!`](@ref) to store `r` in `gauss` if not, and for
+more details dependending on the type of `gauss`.
+"""
+gaussian_elimination(gauss::IterativeGaussianElimination, r::Vector{Int}) = first(_gaussian_elimination(gauss, r))
+
+# For an IterativeGaussianElimination of i rings of size at most l, symdiff_cycles cost at
+# most (l + lenr) + (2l + lenr) + ... + (il + lenr) = O(i^2*l +i*lenr)
+# If the size of the k-th ring is at most k*l, then the worst-case complexity is O(i^3*l)
+# Calling gaussian_elimination ν times thus leads to a worst-case complexity in O(ν^4*l)
+# The reasonning can probably be refined...
+"""
+    gaussian_elimination!(gauss::IterativeGaussianElimination, r::Vector{Int})
+
+Test whether `r` can be expressed as a sum of vectors stored in `gauss`, and store `r` if
+not. "sum" refers to the symmetric difference of boolean vectors, represented in sparse
+format as the ordered list of non-zero indices.
+
+If `gauss isa IterativeGaussianEliminationLength`, return whether `r` can be expressed as a
+sum of strictly smaller vectors.
+
+Otherwise, return `true` when `r` is a sum of any previously encoutered vectors.
+If `gauss` isa `IterativeGaussianEliminationDecomposition`, query `retrieve_track(gauss)`
+to obtain the sorted list of indices of such previously encountered vectors.
+
+See also [`gaussian_elimination`](@ref) to test `r` without storing it.
+"""
+function gaussian_elimination!(gauss::IterativeGaussianElimination{T}, r::Vector{Int}) where T
+    rings = gauss.rings
+    shortcuts = gauss.shortcuts
+    if gauss isa IterativeGaussianEliminationLength
+        len = length(r) % UInt8
+    elseif gauss isa IterativeGaussianEliminationDecomposition
+        track::Vector{Int32} = first(gauss.track)
+    end
+
+    notindependent, r1, maxlen, buffer1 = _gaussian_elimination(gauss, r)
+
+    if notindependent
+        gauss isa IterativeGaussianEliminationLength && return maxlen < len
+        if gauss isa IterativeGaussianEliminationDecomposition
+            push!(rings, buffer1) # dummy, used to keep track of dependent rings
+            push!(last(gauss.track), track) # also dummy
+        end
+        return true
+    end
+
     push!(rings, copy(buffer1))
     r1 > length(shortcuts) && append!(shortcuts, zero(Int32) for _ in 1:(r1-length(shortcuts)))
     shortcuts[r1] = length(rings) % Int32
     if gauss isa IterativeGaussianEliminationLength
         # the ring was not a sum of strictly smaller ones (since it's not a sum of previous ones at all)
-        push!(lengths, len)
+        push!(gauss.track, len)
     elseif gauss isa IterativeGaussianEliminationDecomposition
         push!(last(gauss.track), copy(sort!(track)))
     end
     return false # the new ring was independent from the other ones
-
-    @label notindependentreturn
-    gauss isa IterativeGaussianEliminationLength && return maxlen < len
-    if gauss isa IterativeGaussianEliminationDecomposition
-        push!(rings, buffer2) # dummy, used to keep track of dependent rings
-        push!(last(gauss.track), track) # also dummy
-    end
-    return true
 end
 
 """

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -736,6 +736,56 @@ end
     end
 end
 
+@testset "Iterative gaussian elimination" begin
+    gausslengths = PeriodicGraphs.IterativeGaussianEliminationLength([3, 4, 7])
+    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [3, 6])
+    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [4, 7])
+    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [1, 3, 5])
+    @test length(gausslengths.track) == length(gausslengths.rings) == 4
+    @test PeriodicGraphs.gaussian_elimination(gausslengths, [1, 4, 5, 7])
+    @test PeriodicGraphs.gaussian_elimination!(gausslengths, [1, 4, 5, 7])
+    @test length(gausslengths.track) == length(gausslengths.rings) == 4
+    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [3, 5, 6, 7])
+    @test !PeriodicGraphs.gaussian_elimination(gausslengths, [2, 3])
+    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [2, 3])
+    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [1, 2, 3])
+    @test length(gausslengths.track) == length(gausslengths.rings) == 7
+
+    gaussnone = PeriodicGraphs.IterativeGaussianElimination([3, 4, 7])
+    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [3, 6])
+    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [4, 7])
+    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [1, 3, 5])
+    @test length(gaussnone.rings) == 4
+    @test PeriodicGraphs.gaussian_elimination(gaussnone, [1, 4, 5, 7])
+    @test PeriodicGraphs.gaussian_elimination!(gaussnone, [1, 4, 5, 7])
+    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [3, 5, 6, 7])
+    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [2, 3])
+    @test !PeriodicGraphs.gaussian_elimination(gaussnone, [1, 2, 3])
+    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [1, 2, 3])
+    @test PeriodicGraphs.gaussian_elimination!(gaussnone, [1, 2, 3])
+    @test gaussnone.shortcuts == Int32[4, 6, 1, 2, 5, 3, 7]
+    length(gaussnone.rings) == 7
+
+    gausstrack = PeriodicGraphs.IterativeGaussianEliminationDecomposition([3, 4, 7])
+    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [3, 6])
+    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [4, 7])
+    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [1, 3, 5])
+    @test length(gausstrack.rings) == 4
+    @test PeriodicGraphs.gaussian_elimination(gausstrack, [1, 4, 5, 7])
+    @test PeriodicGraphs.gaussian_elimination!(gausstrack, [1, 4, 5, 7])
+    @test PeriodicGraphs.retrieve_track!(gausstrack) == Int32[5, 4, 1]
+    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [3, 5, 6, 7])
+    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [2, 3])
+    @test !PeriodicGraphs.gaussian_elimination(gausstrack, [1, 2, 3])
+    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [1, 2, 3])
+    @test PeriodicGraphs.gaussian_elimination!(gausstrack, [1, 2, 3])
+    @test PeriodicGraphs.retrieve_track!(gausstrack) == Int32[9, 8]
+    @test gausstrack.shortcuts == Int32[4, 7, 1, 2, 6, 3, 8]
+    @test length(gausstrack.rings) == 9
+
+    @test PeriodicGraphs.IterativeGaussianElimination().track == PeriodicGraphs.IterativeGaussianElimination{Nothing}().track == nothing
+end
+
 const cpi = PeriodicGraph("2 1 2 0 0 1 3 0 0 1 4 0 0 1 5 0 0 2 3 0 -1 2 5 1 -1 3 4 1 0 4 5 0 -1");
 const nab = PeriodicGraph("3 1 2 0 -1 0 1 3 0 -1 -1 1 4 -1 0 0 1 5 -1 0 -1 2 3 0 0 -1 2 4 -1 0 0 2 5 -1 0 -1 3 4 0 0 0 3 5 0 0 0 4 5 0 0 0");
 const cai = PeriodicGraph("3 1 1 1 0 0 1 2 0 0 0 1 3 0 0 0 1 4 0 0 0 1 5 0 0 0 2 2 1 0 0 2 5 0 -1 0 2 6 0 0 0 3 3 1 0 0 3 4 0 0 1 3 6 0 0 1 4 4 1 0 0 4 6 0 1 0 5 5 1 0 0 5 6 0 1 1 6 6 1 0 0");
@@ -875,48 +925,6 @@ end
     @test PeriodicGraphs.union_cycles([3,4], [4,5,6]) == PeriodicGraphs.union_cycles([4,5,6], [3,4]) == [3,4,5,6]
     @test PeriodicGraphs.union_cycles([4,6], [3,5,6]) == PeriodicGraphs.union_cycles([3,5,6], [3,4,6]) == [3,4,5,6]
 
-    gausslengths = PeriodicGraphs.IterativeGaussianEliminationLength([3, 4, 7])
-    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [3, 6])
-    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [4, 7])
-    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [1, 3, 5])
-    @test length(gausslengths.track) == length(gausslengths.rings) == 4
-    @test PeriodicGraphs.gaussian_elimination!(gausslengths, [1, 4, 5, 7])
-    @test length(gausslengths.track) == length(gausslengths.rings) == 4
-    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [3, 5, 6, 7])
-    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [2, 3])
-    @test !PeriodicGraphs.gaussian_elimination!(gausslengths, [1, 2, 3])
-    @test length(gausslengths.track) == length(gausslengths.rings) == 7
-
-    gaussnone = PeriodicGraphs.IterativeGaussianElimination([3, 4, 7])
-    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [3, 6])
-    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [4, 7])
-    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [1, 3, 5])
-    @test length(gaussnone.rings) == 4
-    @test PeriodicGraphs.gaussian_elimination!(gaussnone, [1, 4, 5, 7])
-    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [3, 5, 6, 7])
-    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [2, 3])
-    @test !PeriodicGraphs.gaussian_elimination!(gaussnone, [1, 2, 3])
-    @test PeriodicGraphs.gaussian_elimination!(gaussnone, [1, 2, 3])
-    @test gaussnone.shortcuts == Int32[4, 6, 1, 2, 5, 3, 7]
-    length(gaussnone.rings) == 7
-
-    gausstrack = PeriodicGraphs.IterativeGaussianEliminationDecomposition([3, 4, 7])
-    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [3, 6])
-    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [4, 7])
-    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [1, 3, 5])
-    @test length(gausstrack.rings) == 4
-    @test PeriodicGraphs.gaussian_elimination!(gausstrack, [1, 4, 5, 7])
-    @test PeriodicGraphs.retrieve_track!(gausstrack) == Int32[5, 4, 1]
-    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [3, 5, 6, 7])
-    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [2, 3])
-    @test !PeriodicGraphs.gaussian_elimination!(gausstrack, [1, 2, 3])
-    @test PeriodicGraphs.gaussian_elimination!(gausstrack, [1, 2, 3])
-    @test PeriodicGraphs.retrieve_track!(gausstrack) == Int32[9, 8]
-    @test gausstrack.shortcuts == Int32[4, 7, 1, 2, 6, 3, 8]
-    @test length(gausstrack.rings) == 9
-
-    @test PeriodicGraphs.IterativeGaussianElimination().track == PeriodicGraphs.IterativeGaussianElimination{Nothing}().track == nothing
-
     # keep track of the limitations
     @test_throws ErrorException rings(lta, 63)
     very_high_degree = PeriodicGraph{0}(128)