From 67efbccfa66573312224f0e10bb5a9e7bf370796 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Mon, 4 Nov 2024 12:14:52 -0500 Subject: [PATCH 01/23] Change permutation in gray code min dist and add Brouwer early exit for information set --- src/Classical/linear_code.jl | 24 ++---- src/Classical/weight_dist.jl | 149 ++++++++++++++++++++++++----------- src/utils.jl | 45 +++++++++-- 3 files changed, 148 insertions(+), 70 deletions(-) diff --git a/src/Classical/linear_code.jl b/src/Classical/linear_code.jl index 29ad645..f25d71e 100644 --- a/src/Classical/linear_code.jl +++ b/src/Classical/linear_code.jl @@ -450,7 +450,7 @@ Return a random [n, k] linear code over F. - `k` - dimension of the code - `rng` - random number generator """ -function random_linear_code(F::CTFieldTypes, n::Int, k::Int; rng::AbstractRNG = Random.seed!()) +function random_linear_code(F::CTFieldTypes, n::Int, k::Int; rng::AbstractRNG = Random.seed!(), brute_force_WE = false) rand_mat = zero_matrix(F, k, n - k) for r in 1:nrows(rand_mat) for c in 1:ncols(rand_mat) @@ -458,7 +458,8 @@ function random_linear_code(F::CTFieldTypes, n::Int, k::Int; rng::AbstractRNG = end end full_mat = hcat(identity_matrix(F, k), rand_mat) - return LinearCode(full_mat) + parity = false + return LinearCode(full_mat, parity, brute_force_WE) end """ @@ -583,18 +584,7 @@ function show(io::IO, C::AbstractLinearCode) G = generator_matrix(C) nr, nc = size(G) println(io, "Generator matrix: $nr × $nc") - for i in 1:nr - print(io, "\t") - for j in 1:nc - if j != nc - print(io, "$(G[i, j]) ") - elseif j == nc && i != nr - println(io, "$(G[i, j])") - else - print(io, "$(G[i, j])") - end - end - end + println(io, "\t" * replace(repr(MIME("text/plain"), G), r"\n" => "\n\t")) end end # if !ismissing(C.weight_enum) @@ -640,8 +630,8 @@ end Return the syndrome of `v` with respect to `C`. """ -function syndrome(C::AbstractLinearCode, v::Union{CTMatrixTypes, Vector{Int}}) - w = isa(v, Vector{Int}) ? matrix(C.F, length(v), 1, v) : v +function syndrome(C::AbstractLinearCode, v::Union{CTMatrixTypes, Vector{Int}, Vector{fpFieldElem}, Vector{FpFieldElem}}) + w = isa(v, Union{Vector{Int}, Vector{fpFieldElem}, Vector{FpFieldElem}}) ? matrix(C.F, length(v), 1, v) : v H = parity_check_matrix(C) nc = ncols(H) (size(w) != (nc, 1) && size(w) != (1, nc)) && @@ -662,7 +652,7 @@ end Return whether or not `v` is a codeword of `C`. """ -in(v::Union{CTMatrixTypes, Vector{Int}}, C::AbstractLinearCode) = iszero(syndrome(C, v)) +in(v::Union{CTMatrixTypes, Vector{Int}, Vector{fpFieldElem}, Vector{FpFieldElem}}, C::AbstractLinearCode) = iszero(syndrome(C, v)) """ ⊆(C1::AbstractLinearCode, C2::AbstractLinearCode) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index f654c5b..59f8a8a 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -305,16 +305,23 @@ function information_sets(G::CTMatrixTypes, alg::Symbol = :Edmonds; permute::Boo rnk = nr if alg == :Brouwer for i in 0:Int(floor(nc / nr)) - 1 - rnk, Gi, Pi = _rref_col_swap(G, 1:nr, i * rnk + 1:nc) + start_ind = i * nr + 1 + rnk, Gi, Pi = _rref_col_swap(G, 1:nr, start_ind:nc) + if rnk < nr + println("rank ", rnk, " on mat of size ", size(Gi), " start ind ", start_ind) + break + end + if only_A - Ai = Gi[:, setdiff(1:nc, i * nr + 1:(i + 1) * nr)] + Ai = Gi[:, setdiff(1:nc, start_ind:(i + 1) * nr)] push!(gen_mats, Ai) push!(perms, Pi) push!(rnks, rnk) else if permute # permute identities to the front - pivots = collect(i * nr + 1:(i + 1) * nr) + pivots = collect(start_ind:(i + 1) * nr) + @assert det(Gi[:, pivots]) != 0 σ = [pivots; setdiff(1:nc, pivots)] Gi = Gi[:, σ] Pi = Pi[:, σ] @@ -413,17 +420,22 @@ end Return the minimum distance of `C` using a deterministic algorithm based on enumerating constant weight codewords of the binary reflected Gray code. If a word of minimum weight -is found before the lower and upper bounds cross, it is returned; otherwise, `missing` +is found before the lower and upper bounds cross, it is returned; otherwise, the zero vector is returned. """ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol = :auto, - verbose::Bool = false) + verbose::Bool = false, dbg = Dict()) ord_F = Int(order(C.F)) ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes.")) info_set_alg ∈ (:auto, :Brouwer, :Zimmermann, :White, :Chen, :Bouyuklieva, :Edmonds) || throw(ArgumentError("Unknown information set algorithm. Expected `:auto`, `:Brouwer`, `:Zimmermann`, `:White`, `:Chen`, `:Bouyuklieva`, or `:Edmonds`.")) - G = _remove_empty(generator_matrix(C, true), :cols) + generator_matrix(C, true) # ensure G_stand exists + if _has_empty_vec(C.G, :cols) + #TODO err string can instruct the user to construct a new code without 0 cols and tell them the function for that + throw(ArgumentError("Codes with standard form of generator matrix having 0 columns not supported")) + end + G = C.G #TODO remove local var G # generate if not pre-stored parity_check_matrix(C) @@ -441,10 +453,17 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol end # should have no need to permute to find better ranks because of Edmond's? A_mats, perms_mats, rnks = information_sets(G, info_set_alg, permute = true, only_A = false) + println("expect all k", rnks) + A_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Ai)') for Ai in A_mats] perms_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Pi)') for Pi in perms_mats] h = length(A_mats) + # println("Starting loop to refine upper bound. Initial upper bound ", C.u_bound, " num of mats is ", length(A_mats), " dimension ", size(A_mats[1])) rank_defs = zeros(Int, h) + + optimize_upper_bound = false # if this flag is true then we + dbg["exit_r"] = -1 + if verbose print("Generated $h information sets with ranks: ") for i in 1:h @@ -470,19 +489,29 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol end k, n = size(G) - found = zeros(UInt8, n) - perm = collect(1:n) - # while flag - for (j, g) in enumerate(A_mats) - # can make this faster with dots and views - w, i = _min_wt_col(g) - if w <= C.u_bound - found = g[:, i] - C.u_bound = w - perm = perms_mats[j] - end + found = zeros(Int, n) + if optimize_upper_bound + # C_orig_U_bound = C.u_bound #TODO restore this before returning if its smaller than C.u_bound at the end + C.u_bound = n + end + # initial_perm_ind will match the permutation we use for the 'found' vector if the found vector is nonzero. To simplify the code below we're going to choose an initial permutation arbitrarily. + initial_perm_ind = 1 + for (j, g) in enumerate(A_mats) # loop over the A_mats rather than the original G because it would add another case to deal with later + # can make this faster with dots and views + w, i = _min_wt_col(g) + if w <= C.u_bound + found = g[:, i] + println("setting found to the vector ") + println(found) + C.u_bound = w + initial_perm_ind = j + y = perms_mats[initial_perm_ind] * found #TODO shouldnt we multiply by inv(perm) because found is a vector from the permuted matrix + is_in = in(y, C) + @assert is_in end - # end + end + @assert length(A_mats) != 0 + verbose && println("Current upper bound: $(C.u_bound)") verbose && !iszero(found) && println("Found element matching upper bound.") @@ -507,8 +536,10 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol verbose && println("Upper bound: $(C.u_bound)") if C.l_bound >= C.u_bound C.d = C.u_bound - y = perm * found #TODO shouldnt we multiply by inv(perm) because found is a vector from the permuted matrix + y = perms_mats[initial_perm_ind] * found #TODO shouldnt we multiply by inv(perm) because found is a vector from the permuted matrix + @assert in(y, C) #TODO make into assertion ? iszero(C.H * transpose(y)) + dbg["exit_r"] = r return C.u_bound, y end @@ -526,50 +557,72 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol p = Int(characteristic(C.F)) uppers = [C.u_bound for _ in 1:num_thrds] founds = [found for _ in 1:num_thrds] - perms = [perm for _ in 1:num_thrds] - Threads.@threads for m in 1:num_thrds - c = zeros(Int, C.n) - prefix = digits(m - 1, base = 2, pad = power) - for u in GrayCode(C.k, r, prefix, mutate = true) - if flag[] - for i in 1:h - if r - rank_defs[i] > 0 - LinearAlgebra.mul!(c, A_mats[i], u) - w = r - @inbounds for j in 1:n - c[j] % p != 0 && (w += 1;) - end + exit_thread_indicator_vec = [initial_perm_ind for _ in 1:num_thrds] + use_sub_gray = false - if uppers[m] > w - uppers[m] = w - founds[m] .= c - perms[m] = i - verbose && println("Adjusting (local) upper bound: $w") - if C.l_bound == uppers[m] - Threads.atomic_cas!(flag, true, false) - break - else - r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds) - isnothing(r_term) && (r_term = k;) - verbose && println("Updated termination weight: $r_term") - end + c = zeros(Int, C.n) + if use_sub_gray + itrs = [SubsetGrayCode(C.k, r)] # no threads yet + vec = collect(1:C.k) # assumes 1 thread for now + else + itrs = [GrayCode(C.k, r, digits(m - 1, base = 2, pad = power), mutate = true) for m in 1:num_thrds] + end + # Threads.@threads for itr in itrs + first_itr = true + for m in 1:num_thrds + itr = itrs[m] + for u in itr + if use_sub_gray + if first_itr + first_itr = false + #TODO set vec + else + # its not necc to build vec explicitly but Ill do it for now for testing, at least + for ind in u + if ind != -1 + vec[ind] = vec[ind]==0 ? 1 : 0 end end end else - break + vec = u + end + + for i in 1:h + if r - rank_defs[i] > 0 + LinearAlgebra.mul!(c, A_mats[i], vec) + w = r + @inbounds for j in 1:n + c[j] % p != 0 && (w += 1;) + end + + if uppers[m] > w + uppers[m] = w + founds[m] = c + exit_thread_indicator_vec[m] = i + verbose && println("Adjusting (local) upper bound: $w") + if C.l_bound == uppers[m] + Threads.atomic_cas!(flag, true, false) + break + else + r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds) + isnothing(r_term) && (r_term = k;) + verbose && println("Updated termination weight: $r_term") + end + end + end end end end loc = argmin(uppers) C.u_bound = uppers[loc] found = founds[loc] - perm = perms[loc] + initial_perm_ind = exit_thread_indicator_vec[loc] end # at this point we are guaranteed to have found the answer C.d = C.u_bound - y = matrix(C.F, 1, n, found)[perm] + y = perms_mats[initial_perm_ind] * found # iszero(C.H * transpose(y)) return C.u_bound, y end diff --git a/src/utils.jl b/src/utils.jl index 737be13..d5ff391 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -372,6 +372,37 @@ end # return maxlen # end +function _has_empty_vec(A::Union{CTMatrixTypes, Matrix{<: Number}, BitMatrix, Matrix{Bool}}, + type::Symbol) + + type ∈ (:rows, :cols) || throw(ArgumentError("Unknown type in _remove_empty")) + + del = Vector{Int}() + if type == :rows + for r in axes(A, 1) + flag = true + for c in axes(A, 2) + if !iszero(A[r, c]) + flag = false + break + end + end + flag && append!(del, r) + end + elseif type == :cols + for c in axes(A, 2) + flag = true + for r in axes(A, 1) + if !iszero(A[r, c]) + flag = false + break + end + end + flag && append!(del, c) + end + end + return !isempty(del) +end function _remove_empty(A::Union{CTMatrixTypes, Matrix{<: Number}, BitMatrix, Matrix{Bool}}, type::Symbol) @@ -1962,11 +1993,15 @@ function _rand_invertible_matrix(F::CTFieldTypes, n::Integer) end function extended_binomial(x::UInt, y::UInt) - z = UInt(0) - if y <= x - z = binomial(big(x), big(y)) - end - return z + z = UInt(0) + if y <= x + z = binomial(big(x), big(y)) + end + return z +end + +function _value_distribution(vals) + return OrderedDict([(i, count(x -> (x == i), vals)) for i in collect(sort(unique(vals)))]) end # #= From 5ec69a4c8a83b87ad8b6fccf88761baadb38dca8 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Tue, 5 Nov 2024 19:57:07 -0500 Subject: [PATCH 02/23] a bug is living in here somewhere --- src/Classical/linear_code.jl | 46 ------ src/Classical/weight_dist.jl | 278 ++++++++++++++++++++++++++++++----- src/CodingTheory.jl | 2 +- src/iterators.jl | 4 +- test/iterators_test.jl | 2 +- 5 files changed, 245 insertions(+), 87 deletions(-) diff --git a/src/Classical/linear_code.jl b/src/Classical/linear_code.jl index 248c3cd..a0e061e 100644 --- a/src/Classical/linear_code.jl +++ b/src/Classical/linear_code.jl @@ -481,52 +481,6 @@ function random_information_set(C::AbstractLinearCode; rng::AbstractRNG = Random return on_sets(collect(1 : C.k), inv_perm) end -""" - random_linear_code(F::CTFieldTypes, n::Int, k::Int; rng::AbstractRNG = Random.seed!()) - -Return a random [n, k] linear code over F. - -# Arguments -- `F` - finite field -- `n` - length of the code -- `k` - dimension of the code -- `rng` - random number generator -""" -function random_linear_code(F::CTFieldTypes, n::Int, k::Int; rng::AbstractRNG = Random.seed!(), brute_force_WE = false) - rand_mat = zero_matrix(F, k, n - k) - for r in 1:nrows(rand_mat) - for c in 1:ncols(rand_mat) - rand_mat[r, c] = rand(rng, F) - end - end - full_mat = hcat(identity_matrix(F, k), rand_mat) - parity = false - return LinearCode(full_mat, parity, brute_force_WE) -end - -""" - random_linear_code(q::Int, n::Int, k::Int; rng::AbstractRNG = Random.seed!()) - -Return a random [n, k] linear code over GF(q). - -# Arguments -- `q` - a prime power -- `n` - length of the code -- `k` - dimension of the code -- `rng` - random number generator -""" -function random_linear_code(q::Int, n::Int, k::Int; rng::AbstractRNG = Random.seed!()) - is_pp, e, p = is_prime_power_with_data(q) - if e == 1 - field = Oscar.Nemo.Native.GF(p) - elseif e > 1 - field = GF(p, e, :x) - else - throw(ArgumentError("q must be a prime power")) - end - return random_linear_code(field, n, k, rng = rng) -end - function _standard_form(G::CTMatrixTypes) rnk, G_stand, P = _rref_col_swap(G, 1:nrows(G), 1:ncols(G)) F = base_ring(G_stand) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index aa9f6c8..72f152b 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -411,10 +411,6 @@ function information_set_lower_bound(r::Int, n::Int, k::Int, l::Int, rank_defs:: return lower end -# In the thesis, the h and GrayCode for loops are switched. This seems inefficient -# but does allow one to update the lower bound using a single matrix after every -# GrayCode loop has finished. I believe only enumerating the very expensive GrayCode -# once will overcome this. Could be tested though. """ Gray_code_minimum_distance(C::AbstractLinearCode; verbose::Bool = false) @@ -460,8 +456,9 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol # println("Starting loop to refine upper bound. Initial upper bound ", C.u_bound, " num of mats is ", length(A_mats), " dimension ", size(A_mats[1])) rank_defs = zeros(Int, h) - optimize_upper_bound = false # if this flag is true then we + optimize_upper_bound = false # true will be used in production. False makes it less likely that the algo will early exit at some trivial value dbg["exit_r"] = -1 + dbg["found_codewords"] = [] k, n = size(G) if info_set_alg == :Brouwer && rnks[h] != k @@ -506,7 +503,204 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol found = g[:, i] C.u_bound = w initial_perm_ind = j - y = perms_mats[initial_perm_ind] * found #TODO shouldnt we multiply by inv(perm) because found is a vector from the permuted matrix + y = perms_mats[initial_perm_ind] * found + push!(dbg["found_codewords"], y) + # println(dbg["found_codewords"]) + # is_in = in(y, C) + # @assert is_in + end + end + # @assert length(A_mats) != 0 + + verbose && println("Current upper bound: $(C.u_bound)") + verbose && !iszero(found) && println("Found element matching upper bound.") + + info_set_alg == :Chen ? (l = C.l;) : (l = 0;) + lower_bounds = [information_set_lower_bound(r, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k] + r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds) + isnothing(r_term) && (r_term = k;) + verbose && println("Predicted termination weight based on current upper bound: $r_term") + + num_thrds = Threads.nthreads() + verbose && println("Detected $num_thrds threads.") + power = Int(floor(log(2, num_thrds))) + + for r in 1:k + C.l_bound < lower_bounds[r] && (C.l_bound = lower_bounds[r];) + # an even code can't have have an odd minimum weight + # (!triply_even_flag && !doubly_even_flag && even_flag) && (C.l_bound += C.l_bound % 2;) + # (!triply_even_flag && doubly_even_flag) && (C.l_bound += 4 - C.l_bound % 4;) + # triply_even_flag && (C.l_bound += 8 - C.l_bound % 8;) + verbose && println("r: $r") + verbose && println("Lower bound: $(C.l_bound)") + verbose && println("Upper bound: $(C.u_bound)") + if C.l_bound >= C.u_bound + C.d = C.u_bound + y = perms_mats[initial_perm_ind] * found + # @assert in(y, C) + dbg["exit_r"] = r + return C.u_bound, y + end + + if verbose + count = 0 + for i in 1:h + r - rank_defs[i] ≤ 0 && (count += 1;) + end + count > 0 && println("$count of the original $h information sets no longer contribute to the lower bound") + end + + flag = Threads.Atomic{Bool}(true) + num_thrds = 1 + + p = Int(characteristic(C.F)) + uppers = [C.u_bound for _ in 1:num_thrds] + founds = [found for _ in 1:num_thrds] + exit_thread_indicator_vec = [initial_perm_ind for _ in 1:num_thrds] + + c = zeros(Int, C.n) + itrs = [GrayCode(C.k, r, digits(m - 1, base = 2, pad = power), mutate = true) for m in 1:num_thrds] + # Threads.@threads for itr in itrs + first_itr = true + for m in 1:num_thrds + itr = itrs[m] + # for u in itr #TODO its easier to test the iteration if the loops are in the order used by GW + for i in 1:h + # for i in 1:h #TODO + for u in itr + vec = u + if r - rank_defs[i] > 0 + LinearAlgebra.mul!(c, A_mats[i], vec) + + # for testing + F = Oscar.Nemo.Native.GF(2) + output_vec = F.(Vector(perms_mats[initial_perm_ind] * c)) + push!(dbg["found_codewords"], (r, deepcopy(vec), output_vec)) + + w = r + @inbounds for j in 1:n #TODO unless weve cropped the Gi down to the Ai shouldnt this be skipping the indexs corresp to the identity submatrix + c[j] % p != 0 && (w += 1;) + end + + if uppers[m] > w + uppers[m] = w + founds[m] = c + exit_thread_indicator_vec[m] = i + verbose && println("Adjusting (local) upper bound: $w") + if C.l_bound == uppers[m] + Threads.atomic_cas!(flag, true, false) + break + else + r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds) + isnothing(r_term) && (r_term = k;) + verbose && println("Updated termination weight: $r_term") + end + end + end + end + end + end + loc = argmin(uppers) + C.u_bound = uppers[loc] + found = founds[loc] + initial_perm_ind = exit_thread_indicator_vec[loc] + end + + # at this point we are guaranteed to have found the answer + C.d = C.u_bound + y = perms_mats[initial_perm_ind] * found + # iszero(C.H * transpose(y)) + return C.u_bound, y +end + +function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol = :auto, + verbose::Bool = false, dbg = Dict()) + + ord_F = Int(order(C.F)) + ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes.")) + info_set_alg ∈ (:auto, :Brouwer, :Zimmermann, :White, :Chen, :Bouyuklieva, :Edmonds) || throw(ArgumentError("Unknown information set algorithm. Expected `:auto`, `:Brouwer`, `:Zimmermann`, `:White`, `:Chen`, `:Bouyuklieva`, or `:Edmonds`.")) + + generator_matrix(C, true) # ensure G_stand exists + if _has_empty_vec(C.G, :cols) + #TODO err string can instruct the user to construct a new code without 0 cols and tell them the function for that + throw(ArgumentError("Codes with standard form of generator matrix having 0 columns not supported")) + end + G = C.G #TODO remove local var G + # generate if not pre-stored + parity_check_matrix(C) + + if info_set_alg == :auto + if typeof(C) <: AbstractCyclicCode + verbose && println("Detected a cyclic code, using Chen's adaption.") + info_set_alg = :Chen + # TODO: fix this case + elseif typeof(C) <: AbstractQuasiCyclicCode + verbose && println("Detected a quasi-cyclic code, using White's adaption.") + info_set_alg = :White + else + info_set_alg = :Edmonds + end + end + # should have no need to permute to find better ranks because of Edmond's? + A_mats, perms_mats, rnks = information_sets(G, info_set_alg, permute = true, only_A = false) + + A_mats_gf2 = [transpose(deepcopy(Ai)) for Ai in A_mats] + A_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Ai)') for Ai in A_mats] + perms_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Pi)') for Pi in perms_mats] + h = length(A_mats) + # println("Starting loop to refine upper bound. Initial upper bound ", C.u_bound, " num of mats is ", length(A_mats), " dimension ", size(A_mats[1])) + rank_defs = zeros(Int, h) + + optimize_upper_bound = false # true will be used in production. False makes it less likely that the algo will early exit at some trivial value + dbg["exit_r"] = -1 + dbg["found_codewords"] = [] + + k, n = size(G) + if info_set_alg == :Brouwer && rnks[h] != k + println("Rank of last matrix too small") + return + end + if verbose + print("Generated $h information sets with ranks: ") + for i in 1:h + i == h ? (println(rnks[i]);) : (print("$(rnks[i]), ");) + # will only be using the rank deficits here + # at the moment, the information sets are always disjoint so the relative + # rank is zero + # TODO huh? check this comment and setup properly + rank_defs[i] = C.k - rnks[i] + end + end + + even_flag = false + doubly_even_flag = false + triply_even_flag = false + ord_F == 2 && (even_flag = is_even(C);) + even_flag && (doubly_even_flag = is_doubly_even(C);) + doubly_even_flag && (triply_even_flag = is_triply_even(C);) + if verbose + triply_even_flag && println("Detected a triply even code.") + (!triply_even_flag && doubly_even_flag) && println("Detected a doubly even code.") + (!triply_even_flag && !doubly_even_flag && even_flag) && println("Detected an even code.") + end + + found = zeros(Int, n) + if optimize_upper_bound + # C_orig_U_bound = C.u_bound #TODO restore this before returning if its smaller than C.u_bound at the end + C.u_bound = n + end + # initial_perm_ind will match the permutation we use for the 'found' vector if the found vector is nonzero. To simplify the code below we're going to choose an initial permutation arbitrarily. + initial_perm_ind = 1 + for (j, g) in enumerate(A_mats) # loop over the A_mats rather than the original G because it would add another case to deal with later + # can make this faster with dots and views + w, i = _min_wt_col(g) + if w <= C.u_bound + found = g[:, i] + C.u_bound = w + initial_perm_ind = j + y = perms_mats[initial_perm_ind] * found + push!(dbg["found_codewords"], y) + println(dbg["found_codewords"]) is_in = in(y, C) @assert is_in end @@ -537,9 +731,8 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol verbose && println("Upper bound: $(C.u_bound)") if C.l_bound >= C.u_bound C.d = C.u_bound - y = perms_mats[initial_perm_ind] * found #TODO shouldnt we multiply by inv(perm) because found is a vector from the permuted matrix + y = perms_mats[initial_perm_ind] * found @assert in(y, C) - #TODO make into assertion ? iszero(C.H * transpose(y)) dbg["exit_r"] = r return C.u_bound, y end @@ -553,48 +746,60 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol end flag = Threads.Atomic{Bool}(true) - # num_thrds = 1 - # power = 0 + num_thrds = 1 + p = Int(characteristic(C.F)) uppers = [C.u_bound for _ in 1:num_thrds] founds = [found for _ in 1:num_thrds] exit_thread_indicator_vec = [initial_perm_ind for _ in 1:num_thrds] - use_sub_gray = false - c = zeros(Int, C.n) - if use_sub_gray - itrs = [SubsetGrayCode(C.k, r)] # no threads yet - vec = collect(1:C.k) # assumes 1 thread for now - else - itrs = [GrayCode(C.k, r, digits(m - 1, base = 2, pad = power), mutate = true) for m in 1:num_thrds] - end + itrs = [SubsetGrayCode(C.k, r)] # Threads.@threads for itr in itrs - first_itr = true for m in 1:num_thrds itr = itrs[m] - for u in itr - if use_sub_gray - if first_itr - first_itr = false - #TODO set vec - else + # for u in itr #TODO its easier to test the iteration if the loops are in the order used by GW + for i in 1:h + # for i in 1:h #TODO + vec = C.F.(vcat(fill(1, r), fill(0, C.k - r))) # initial vec corresponds to the subset {1,..,r} + c = zeros(C.F, C.n) + c_itr = zeros(C.F, C.n) + first_itr = true + for u in itr + A_col_inds = [] + if !first_itr # its not necc to build vec explicitly but Ill do it for now for testing, at least for ind in u if ind != -1 - vec[ind] = vec[ind]==0 ? 1 : 0 + vec[ind] = (vec[ind]==0) ? C.F(1) : C.F(0) + push!(A_col_inds, deepcopy(ind)) + if length(A_col_inds) == 0 + println("fail with u=$u") + end end end end - else - vec = u - end - - for i in 1:h if r - rank_defs[i] > 0 - LinearAlgebra.mul!(c, A_mats[i], vec) + c_itr = deepcopy(c) #TODO c is going to be removed then remove this line + for ind in A_col_inds + vec_to_add = A_mats_gf2[i][:, ind] + c_itr = deepcopy(c_itr + vec_to_add) + end + c = A_mats_gf2[i] * vec + if first_itr + first_itr = false + c_itr = c + end + # @assert (c == c_itr) "vec from mul not eq to vec from itr: $(repr(c)) != $(repr(c_itr))" + + # for testing + F = Oscar.Nemo.Native.GF(2) + output_vec = F.(Vector(perms_mats[initial_perm_ind] * c)) + println("pushing ", (r, deepcopy(vec), output_vec)) + push!(dbg["found_codewords"], (r, deepcopy(vec), output_vec)) + w = r - @inbounds for j in 1:n - c[j] % p != 0 && (w += 1;) + @inbounds for j in 1:n #TODO unless weve cropped the Gi down to the Ai shouldnt this be skipping the indexs corresp to the identity submatrix + c[j] != 0 && (w += 1;) end if uppers[m] > w @@ -627,7 +832,6 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol # iszero(C.H * transpose(y)) return C.u_bound, y end - """ words_of_weight(C::AbstractLinearCode, l_bound::Int, u_bound::Int; verbose::Bool = false) @@ -1415,10 +1619,10 @@ function minimum_distance(C::AbstractLinearCode; alg::Symbol = :trellis, sect::B for i in 1:length(HWE.polynomial)])) return C.d else - return Gray_code_minimum_distance(C, verbose = verbose) + return david_Gray_code_minimum_distance(C, verbose = verbose) end elseif alg == :Gray - return Gray_code_minimum_distance(C, verbose = verbose) + return david_Gray_code_minimum_distance(C, verbose = verbose) elseif alg == :trellis weight_enumerator_classical(syndrome_trellis(C, "primal", false), type = type) return C.d diff --git a/src/CodingTheory.jl b/src/CodingTheory.jl index b625534..1a47b51 100644 --- a/src/CodingTheory.jl +++ b/src/CodingTheory.jl @@ -378,7 +378,7 @@ export Trellis, vertices, edges, isisomorphic, isequal, loadbalancedecode, include("Classical/weight_dist.jl") export polynomial, type, CWE_to_HWE, weight_enumerator, MacWilliams_identity, weight_distribution, weight_plot, support, minimum_distance, Sterns_attack, - Gray_code_minimum_distance, minimum_words, words_of_weight + david_Gray_code_minimum_distance, minimum_words, words_of_weight ############################# # Quantum/weight_dist.jl diff --git a/src/iterators.jl b/src/iterators.jl index 949e865..0b4dfc7 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -177,10 +177,10 @@ end function _subset_unrank!(r::BigInt, n::UInt, T::Vector{UInt}) # Based on Algorithm 2.12 in kreher1999combinatorial k = length(T) - subset_size_str = "subset size k = $k must be smaller than the set size n = $n" + subset_size_str::String = "subset size k = $k must be smaller than the set size n = $n" k > n && throw(ArgumentError(subset_size_str)) bnd = binomial(n, k) - rank_size_str = "rank must be in [0, choose(n, k) - 1] = $bnd" + rank_size_str::String = "rank must be in [0, choose(n, k) - 1] = $bnd" r > bnd && throw(ArgumentError(rank_size_str)) x = 0 diff --git a/test/iterators_test.jl b/test/iterators_test.jl index 613f2f2..02ffa22 100644 --- a/test/iterators_test.jl +++ b/test/iterators_test.jl @@ -16,7 +16,7 @@ sort!(all_subsets_gray) all_subsets_hecke = CodingTheory.Oscar.subsets(collect(1:len), weight) sort!(all_subsets_hecke) - @assert all_subsets_gray == all_subsets_hecke + @test all_subsets_gray == all_subsets_hecke end function pushOrDel(a::Set{T}, b::T...) where T <: Any From 09679c6e40397eb5eeb4563199ea14dcdd01222d Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Wed, 6 Nov 2024 15:00:27 -0500 Subject: [PATCH 03/23] Appears to be working --- src/Classical/weight_dist.jl | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 72f152b..359f0c1 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -504,13 +504,8 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol C.u_bound = w initial_perm_ind = j y = perms_mats[initial_perm_ind] * found - push!(dbg["found_codewords"], y) - # println(dbg["found_codewords"]) - # is_in = in(y, C) - # @assert is_in end end - # @assert length(A_mats) != 0 verbose && println("Current upper bound: $(C.u_bound)") verbose && !iszero(found) && println("Found element matching upper bound.") @@ -699,13 +694,8 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S C.u_bound = w initial_perm_ind = j y = perms_mats[initial_perm_ind] * found - push!(dbg["found_codewords"], y) - println(dbg["found_codewords"]) - is_in = in(y, C) - @assert is_in end end - @assert length(A_mats) != 0 verbose && println("Current upper bound: $(C.u_bound)") verbose && !iszero(found) && println("Found element matching upper bound.") @@ -789,7 +779,6 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S first_itr = false c_itr = c end - # @assert (c == c_itr) "vec from mul not eq to vec from itr: $(repr(c)) != $(repr(c_itr))" # for testing F = Oscar.Nemo.Native.GF(2) From dadfaae65bc5a869191d99aa88b68740063fe483 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Fri, 8 Nov 2024 14:49:51 -0500 Subject: [PATCH 04/23] Use Blas --- src/Classical/weight_dist.jl | 99 ++++++++++++++++++++++-------------- src/utils.jl | 3 +- 2 files changed, 62 insertions(+), 40 deletions(-) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 359f0c1..9937ee2 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -7,7 +7,8 @@ ############################# # General ############################# - +using AllocCheck +using Profile """ polynomial(W::WeightEnumerator) @@ -519,6 +520,7 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol num_thrds = Threads.nthreads() verbose && println("Detected $num_thrds threads.") power = Int(floor(log(2, num_thrds))) + store_vecs = haskey(dbg, "found_codewords") for r in 1:k C.l_bound < lower_bounds[r] && (C.l_bound = lower_bounds[r];) @@ -567,10 +569,12 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol if r - rank_defs[i] > 0 LinearAlgebra.mul!(c, A_mats[i], vec) - # for testing - F = Oscar.Nemo.Native.GF(2) - output_vec = F.(Vector(perms_mats[initial_perm_ind] * c)) - push!(dbg["found_codewords"], (r, deepcopy(vec), output_vec)) + #TODO for testing + if store_vecs + F = Oscar.Nemo.Native.GF(2) + output_vec = F.(Vector(perms_mats[initial_perm_ind] * c)) + push!(dbg["found_codewords"], (r, deepcopy(vec), output_vec)) + end w = r @inbounds for j in 1:n #TODO unless weve cropped the Gi down to the Ai shouldnt this be skipping the indexs corresp to the identity submatrix @@ -639,16 +643,21 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S # should have no need to permute to find better ranks because of Edmond's? A_mats, perms_mats, rnks = information_sets(G, info_set_alg, permute = true, only_A = false) - A_mats_gf2 = [transpose(deepcopy(Ai)) for Ai in A_mats] A_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Ai)') for Ai in A_mats] perms_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Pi)') for Pi in perms_mats] h = length(A_mats) # println("Starting loop to refine upper bound. Initial upper bound ", C.u_bound, " num of mats is ", length(A_mats), " dimension ", size(A_mats[1])) rank_defs = zeros(Int, h) - optimize_upper_bound = false # true will be used in production. False makes it less likely that the algo will early exit at some trivial value + if length(keys(dbg)) > 0 + println("Debug mode ON") + end + optimize_uppr_bnd = false #TODO Set this from dbg dictionary. optimize_upper_bound controls the initial upper bound. False makes it less likely that the algo will exit early before any codewords are examined dbg["exit_r"] = -1 - dbg["found_codewords"] = [] + store_found_codewords = haskey(dbg, "found_codewords") + if store_found_codewords + verbose && println("Storing the codewords in the debug dictionary") + end k, n = size(G) if info_set_alg == :Brouwer && rnks[h] != k @@ -680,7 +689,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S end found = zeros(Int, n) - if optimize_upper_bound + if optimize_uppr_bnd # C_orig_U_bound = C.u_bound #TODO restore this before returning if its smaller than C.u_bound at the end C.u_bound = n end @@ -710,6 +719,8 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S verbose && println("Detected $num_thrds threads.") power = Int(floor(log(2, num_thrds))) + work_factor_up_to_log_field_size = 0 + expected_work_factor = 0 for r in 1:k C.l_bound < lower_bounds[r] && (C.l_bound = lower_bounds[r];) # an even code can't have have an odd minimum weight @@ -722,8 +733,9 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S if C.l_bound >= C.u_bound C.d = C.u_bound y = perms_mats[initial_perm_ind] * found - @assert in(y, C) + # @assert in(y, C) dbg["exit_r"] = r + println("work factor ", work_factor_up_to_log_field_size, " vs ", expected_work_factor) return C.u_bound, y end @@ -743,57 +755,66 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S founds = [found for _ in 1:num_thrds] exit_thread_indicator_vec = [initial_perm_ind for _ in 1:num_thrds] + # itrs_alloc = @allocated begin itrs = [SubsetGrayCode(C.k, r)] + # lens = [length(_) for _ in itrs] + # println(lens) + # expected_work_factor += h * sum(lens) + # end + # println("alloc itrs used ", itrs_alloc) + # Threads.@threads for itr in itrs for m in 1:num_thrds itr = itrs[m] # for u in itr #TODO its easier to test the iteration if the loops are in the order used by GW for i in 1:h # for i in 1:h #TODO - vec = C.F.(vcat(fill(1, r), fill(0, C.k - r))) # initial vec corresponds to the subset {1,..,r} - c = zeros(C.F, C.n) - c_itr = zeros(C.F, C.n) - first_itr = true + c_itr = zeros(Int, C.n) + is_first = true + curr_mat = A_mats[i] for u in itr - A_col_inds = [] - if !first_itr - # its not necc to build vec explicitly but Ill do it for now for testing, at least + work_factor_up_to_log_field_size += 1 + if r - rank_defs[i] > 0 for ind in u if ind != -1 - vec[ind] = (vec[ind]==0) ? C.F(1) : C.F(0) - push!(A_col_inds, deepcopy(ind)) - if length(A_col_inds) == 0 - println("fail with u=$u") + is_first = false + add_allocs = @allocated begin + axpy!(1, curr_mat[:, ind], c_itr) + end + println(add_allocs) + @inbounds @simd for j in 1:n + c_itr[j] %= p end end end - end - if r - rank_defs[i] > 0 - c_itr = deepcopy(c) #TODO c is going to be removed then remove this line - for ind in A_col_inds - vec_to_add = A_mats_gf2[i][:, ind] - c_itr = deepcopy(c_itr + vec_to_add) - end - c = A_mats_gf2[i] * vec - if first_itr - first_itr = false - c_itr = c + if is_first + # 2pm mul+reduce (mod 2) allocates ~450 bytes per iteration + vec = vcat(fill(1, r), fill(0, C.k - r)) # initial vec corresponds to the subset {1,..,r} + # c_itr = A_mats[i] * vec + LinearAlgebra.mul!(c_itr, curr_mat, vec) + @inbounds @simd for j in 1:n + c_itr[j] %= p + end + # @inbounds for j in 1:n + # c_itr[j] > 1 && (c_itr[j] %= 2) + # end end - # for testing - F = Oscar.Nemo.Native.GF(2) - output_vec = F.(Vector(perms_mats[initial_perm_ind] * c)) - println("pushing ", (r, deepcopy(vec), output_vec)) - push!(dbg["found_codewords"], (r, deepcopy(vec), output_vec)) + if store_found_codewords + output_vec = perms_mats[initial_perm_ind] * c_itr + push!(dbg["found_codewords"], (r, [], output_vec)) + end w = r @inbounds for j in 1:n #TODO unless weve cropped the Gi down to the Ai shouldnt this be skipping the indexs corresp to the identity submatrix - c[j] != 0 && (w += 1;) + # reduce c_itr and store + c_itr .%= p + c_itr[j] != 0 && (w += 1;) end if uppers[m] > w uppers[m] = w - founds[m] = c + founds[m] = c_itr exit_thread_indicator_vec[m] = i verbose && println("Adjusting (local) upper bound: $w") if C.l_bound == uppers[m] diff --git a/src/utils.jl b/src/utils.jl index 5688ad9..baa9e06 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -164,7 +164,8 @@ function _min_wt_row(A::Union{CTMatrixTypes, Matrix{S}, LinearAlgebra.Adjoint{S, return w, i end -function _min_wt_col(A::Union{CTMatrixTypes, Matrix{S}, LinearAlgebra.Adjoint{S, Matrix{S}}}) where S <: Integer +# function _min_wt_col(A::Union{CTMatrixTypes, Matrix{S}, LinearAlgebra.Adjoint{S, Matrix{S}}}) where S <: Integer +function _min_wt_col(A::Union{CTMatrixTypes, Matrix{S}, LinearAlgebra.Adjoint{S, Matrix{S}}, LinearAlgebra.Adjoint{Bool, BitMatrix}}) where S <: Integer w = size(A, 1) + 1 i = 0 for c in axes(A, 2) From ac2f74502549afdfc36335d7aea7224b0fe95135 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Tue, 12 Nov 2024 21:58:19 -0500 Subject: [PATCH 05/23] Remove allocs --- src/Classical/weight_dist.jl | 73 +++++++++++++++++++++--------------- 1 file changed, 42 insertions(+), 31 deletions(-) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 9937ee2..01421dd 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -308,8 +308,7 @@ function information_sets(G::CTMatrixTypes, alg::Symbol = :Edmonds; permute::Boo for i in 0:Int(floor(nc / nr)) - 1 start_ind = i * nr + 1 rnk, Gi, Pi = _rref_col_swap(G, 1:nr, start_ind:nc) - if rnk < nr - println("rank ", rnk, " on mat of size ", size(Gi), " start ind ", start_ind) + if rnk < nr # for Brouwer the Gi must all have full rank break end @@ -457,7 +456,7 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol # println("Starting loop to refine upper bound. Initial upper bound ", C.u_bound, " num of mats is ", length(A_mats), " dimension ", size(A_mats[1])) rank_defs = zeros(Int, h) - optimize_upper_bound = false # true will be used in production. False makes it less likely that the algo will early exit at some trivial value + optimize_upper_bound = true # true will be used in production. False makes it less likely that the algo will early exit at some trivial value dbg["exit_r"] = -1 dbg["found_codewords"] = [] @@ -536,7 +535,7 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol y = perms_mats[initial_perm_ind] * found # @assert in(y, C) dbg["exit_r"] = r - return C.u_bound, y + return C.u_bound, (C.F).(y) #TODO end if verbose @@ -609,7 +608,7 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol C.d = C.u_bound y = perms_mats[initial_perm_ind] * found # iszero(C.H * transpose(y)) - return C.u_bound, y + return C.u_bound, (C.F).(y) end function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol = :auto, @@ -617,6 +616,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S ord_F = Int(order(C.F)) ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes.")) + C.k < 2^16 || throw(DomainError("The given linear code has length k >= 2^16 which is not supported")) info_set_alg ∈ (:auto, :Brouwer, :Zimmermann, :White, :Chen, :Bouyuklieva, :Edmonds) || throw(ArgumentError("Unknown information set algorithm. Expected `:auto`, `:Brouwer`, `:Zimmermann`, `:White`, `:Chen`, `:Bouyuklieva`, or `:Edmonds`.")) generator_matrix(C, true) # ensure G_stand exists @@ -625,6 +625,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S throw(ArgumentError("Codes with standard form of generator matrix having 0 columns not supported")) end G = C.G #TODO remove local var G + # generate if not pre-stored parity_check_matrix(C) @@ -643,8 +644,8 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S # should have no need to permute to find better ranks because of Edmond's? A_mats, perms_mats, rnks = information_sets(G, info_set_alg, permute = true, only_A = false) - A_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Ai)') for Ai in A_mats] - perms_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Pi)') for Pi in perms_mats] + A_mats = [deepcopy(_Flint_matrix_to_Julia_T_matrix(Ai, UInt16)') for Ai in A_mats] + perms_mats = [deepcopy(_Flint_matrix_to_Julia_T_matrix(Pi, UInt16)') for Pi in perms_mats] h = length(A_mats) # println("Starting loop to refine upper bound. Initial upper bound ", C.u_bound, " num of mats is ", length(A_mats), " dimension ", size(A_mats[1])) rank_defs = zeros(Int, h) @@ -652,11 +653,22 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S if length(keys(dbg)) > 0 println("Debug mode ON") end - optimize_uppr_bnd = false #TODO Set this from dbg dictionary. optimize_upper_bound controls the initial upper bound. False makes it less likely that the algo will exit early before any codewords are examined - dbg["exit_r"] = -1 - store_found_codewords = haskey(dbg, "found_codewords") + dbg_key_unoptimized_uppr_bnd = "unoptimized_uppr_bnd" + dbg_key_found = "found_codewords" + dbg_key_exit_r = "exit_r" + + unoptimized_uppr_bnd = false + if haskey(dbg, dbg_key_unoptimized_uppr_bnd) + verbose && println("dbg Dict: Using unoptimized upper bound") + unoptimized_uppr_bnd = true + end + if haskey(dbg, dbg_key_exit_r) + verbose && println("dbg Dict: Storing the largest message weight reached @key=$dbg_key_found") + dbg["exit_r"] = -1 + end + store_found_codewords = haskey(dbg, dbg_key_found) if store_found_codewords - verbose && println("Storing the codewords in the debug dictionary") + verbose && println("dbg Dict: Storing the codewords in the debug dictionary, key=$dbg_key_found") end k, n = size(G) @@ -689,7 +701,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S end found = zeros(Int, n) - if optimize_uppr_bnd + if unoptimized_uppr_bnd # C_orig_U_bound = C.u_bound #TODO restore this before returning if its smaller than C.u_bound at the end C.u_bound = n end @@ -722,6 +734,9 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S work_factor_up_to_log_field_size = 0 expected_work_factor = 0 for r in 1:k + if r > 2^16 + verbose && println("reached an r with r>2^16. This is not implemented yet") + end C.l_bound < lower_bounds[r] && (C.l_bound = lower_bounds[r];) # an even code can't have have an odd minimum weight # (!triply_even_flag && !doubly_even_flag && even_flag) && (C.l_bound += C.l_bound % 2;) @@ -732,10 +747,10 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S verbose && println("Upper bound: $(C.u_bound)") if C.l_bound >= C.u_bound C.d = C.u_bound - y = perms_mats[initial_perm_ind] * found + y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) # @assert in(y, C) - dbg["exit_r"] = r - println("work factor ", work_factor_up_to_log_field_size, " vs ", expected_work_factor) + dbg[dbg_key_exit_r] = r + #TODO report work factors return C.u_bound, y end @@ -769,7 +784,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S # for u in itr #TODO its easier to test the iteration if the loops are in the order used by GW for i in 1:h # for i in 1:h #TODO - c_itr = zeros(Int, C.n) + c_itr = zeros(UInt16, C.n) is_first = true curr_mat = A_mats[i] for u in itr @@ -778,26 +793,26 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S for ind in u if ind != -1 is_first = false - add_allocs = @allocated begin - axpy!(1, curr_mat[:, ind], c_itr) + # 51.4 s (1.27% GC) to evaluate, + @simd for i in eachindex(c_itr) + @inbounds c_itr[i] += curr_mat[i, ind] end - println(add_allocs) @inbounds @simd for j in 1:n c_itr[j] %= p end + # @simd for i in eachindex(c_itr) + # @inbounds c_itr[i] = (c_itr[i] .+ curr_mat[i, ind]).%p + # println(c_itr[i]) + # c_itr[i] = .%p + # end end end if is_first - # 2pm mul+reduce (mod 2) allocates ~450 bytes per iteration vec = vcat(fill(1, r), fill(0, C.k - r)) # initial vec corresponds to the subset {1,..,r} - # c_itr = A_mats[i] * vec LinearAlgebra.mul!(c_itr, curr_mat, vec) @inbounds @simd for j in 1:n c_itr[j] %= p end - # @inbounds for j in 1:n - # c_itr[j] > 1 && (c_itr[j] %= 2) - # end end if store_found_codewords @@ -806,11 +821,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S end w = r - @inbounds for j in 1:n #TODO unless weve cropped the Gi down to the Ai shouldnt this be skipping the indexs corresp to the identity submatrix - # reduce c_itr and store - c_itr .%= p - c_itr[j] != 0 && (w += 1;) - end + w += sum(c_itr) if uppers[m] > w uppers[m] = w @@ -838,9 +849,9 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S # at this point we are guaranteed to have found the answer C.d = C.u_bound - y = perms_mats[initial_perm_ind] * found + y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) # iszero(C.H * transpose(y)) - return C.u_bound, y + return C.u_bound, y_Oscar end """ words_of_weight(C::AbstractLinearCode, l_bound::Int, u_bound::Int; verbose::Bool = false) From 4c8f5a9073544933c2d873249f242bcad4390560 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Thu, 14 Nov 2024 10:02:51 -0500 Subject: [PATCH 06/23] Add progress bar --- src/Classical/weight_dist.jl | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 01421dd..ce1e7ea 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -9,6 +9,7 @@ ############################# using AllocCheck using Profile +using ProgressMeter """ polynomial(W::WeightEnumerator) @@ -612,7 +613,7 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol end function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol = :auto, - verbose::Bool = false, dbg = Dict()) + verbose::Bool = false, dbg = Dict(), show_progress=true) ord_F = Int(order(C.F)) ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes.")) @@ -732,7 +733,12 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S power = Int(floor(log(2, num_thrds))) work_factor_up_to_log_field_size = 0 - expected_work_factor = 0 + predicted_work_factor = fld(n, k) * sum([binomial(k, i) for i in 1:r_term]) + verbose && println("Predicted work factor: $predicted_work_factor") + if show_progress + prog_bar = Progress(predicted_work_factor, dt=1.0, showspeed=true) # updates no faster than once every 1s + end + for r in 1:k if r > 2^16 verbose && println("reached an r with r>2^16. This is not implemented yet") @@ -751,6 +757,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S # @assert in(y, C) dbg[dbg_key_exit_r] = r #TODO report work factors + show_progress && ProgressMeter.finish!(prog_bar) return C.u_bound, y end @@ -789,30 +796,24 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S curr_mat = A_mats[i] for u in itr work_factor_up_to_log_field_size += 1 + show_progress && ProgressMeter.next!(prog_bar) if r - rank_defs[i] > 0 for ind in u + if ind != -1 is_first = false - # 51.4 s (1.27% GC) to evaluate, @simd for i in eachindex(c_itr) @inbounds c_itr[i] += curr_mat[i, ind] + # c_itr[i] %= p end - @inbounds @simd for j in 1:n - c_itr[j] %= p - end - # @simd for i in eachindex(c_itr) - # @inbounds c_itr[i] = (c_itr[i] .+ curr_mat[i, ind]).%p - # println(c_itr[i]) - # c_itr[i] = .%p - # end end end if is_first vec = vcat(fill(1, r), fill(0, C.k - r)) # initial vec corresponds to the subset {1,..,r} LinearAlgebra.mul!(c_itr, curr_mat, vec) - @inbounds @simd for j in 1:n - c_itr[j] %= p - end + # @inbounds @simd for j in 1:n + # c_itr[j] %= p + # end end if store_found_codewords @@ -820,12 +821,11 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S push!(dbg["found_codewords"], (r, [], output_vec)) end - w = r - w += sum(c_itr) + w = sum(c_itr) if uppers[m] > w uppers[m] = w - founds[m] = c_itr + founds[m] = c_itr #TODO reduce founds coeffs here and below exit_thread_indicator_vec[m] = i verbose && println("Adjusting (local) upper bound: $w") if C.l_bound == uppers[m] @@ -846,6 +846,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S found = founds[loc] initial_perm_ind = exit_thread_indicator_vec[loc] end + show_progress && ProgressMeter.finish!(prog_bar) # at this point we are guaranteed to have found the answer C.d = C.u_bound From 290e750cb5abae34d9ef1b24ffcaf37e82482e40 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Thu, 14 Nov 2024 13:58:10 -0500 Subject: [PATCH 07/23] Use xor and partial weight sum --- src/Classical/weight_dist.jl | 53 +++++++++++++++++++++--------------- 1 file changed, 31 insertions(+), 22 deletions(-) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index ce1e7ea..7b88312 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -741,7 +741,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S for r in 1:k if r > 2^16 - verbose && println("reached an r with r>2^16. This is not implemented yet") + verbose && println("Reached an r with r>2^16") end C.l_bound < lower_bounds[r] && (C.l_bound = lower_bounds[r];) # an even code can't have have an odd minimum weight @@ -756,7 +756,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) # @assert in(y, C) dbg[dbg_key_exit_r] = r - #TODO report work factors + #TODO report actual/expected work factors in dbg dict show_progress && ProgressMeter.finish!(prog_bar) return C.u_bound, y end @@ -786,7 +786,15 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S # println("alloc itrs used ", itrs_alloc) # Threads.@threads for itr in itrs + vec = vcat(fill(1, r), fill(0, C.k - r)) # initial vec corresponds to the subset {1,..,r} for m in 1:num_thrds + # if c_itr is a binary vector. If it were a random binary vector the expected + # weight of the subvector c_itr[1:2*uppers[m]] equals uppers[m]. + # So we use the heuristic that w(c_itr[1:2*uppers[m]+5]) is usually larger than uppers[m]. + # If not, we compute weight of all of c_itr + weight_sum_bound = min(2 * uppers[m] + 5, n) + + verbose && println("Computing weight up to $weight_sum_bound") itr = itrs[m] # for u in itr #TODO its easier to test the iteration if the loops are in the order used by GW for i in 1:h @@ -803,17 +811,15 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S if ind != -1 is_first = false @simd for i in eachindex(c_itr) - @inbounds c_itr[i] += curr_mat[i, ind] - # c_itr[i] %= p + @inbounds c_itr[i] = xor(c_itr[i], curr_mat[i, ind]) end end end if is_first - vec = vcat(fill(1, r), fill(0, C.k - r)) # initial vec corresponds to the subset {1,..,r} LinearAlgebra.mul!(c_itr, curr_mat, vec) - # @inbounds @simd for j in 1:n - # c_itr[j] %= p - # end + @inbounds @simd for j in 1:n + c_itr[j] %= p + end end if store_found_codewords @@ -821,20 +827,23 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S push!(dbg["found_codewords"], (r, [], output_vec)) end - w = sum(c_itr) - - if uppers[m] > w - uppers[m] = w - founds[m] = c_itr #TODO reduce founds coeffs here and below - exit_thread_indicator_vec[m] = i - verbose && println("Adjusting (local) upper bound: $w") - if C.l_bound == uppers[m] - Threads.atomic_cas!(flag, true, false) - break - else - r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds) - isnothing(r_term) && (r_term = k;) - verbose && println("Updated termination weight: $r_term") + partial_weight = sum(view(c_itr, 1:weight_sum_bound)) + + if uppers[m] > partial_weight + w = sum(c_itr) + if uppers[m] > w + uppers[m] = w + founds[m] = c_itr #TODO reduce founds coeffs here and below + exit_thread_indicator_vec[m] = i + verbose && println("Adjusting (local) upper bound: $w") + if C.l_bound == uppers[m] + Threads.atomic_cas!(flag, true, false) + break + else + r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds) + isnothing(r_term) && (r_term = k;) + verbose && println("Updated termination weight: $r_term") + end end end end From df95247cf88db0b5ed0f51dc05d7217fed705953 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Mon, 18 Nov 2024 11:37:46 -0500 Subject: [PATCH 08/23] Rename old and new mindist functions --- src/Classical/weight_dist.jl | 397 +++++++++++++++++++---------------- src/CodingTheory.jl | 2 +- 2 files changed, 215 insertions(+), 184 deletions(-) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 7b88312..8e613bd 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -305,6 +305,7 @@ function information_sets(G::CTMatrixTypes, alg::Symbol = :Edmonds; permute::Boo rnks = Vector{Int}() rnk = nr + #TODO Brouw and Zimm should use the same code for rref then discard the last matrix if alg == :Brouwer for i in 0:Int(floor(nc / nr)) - 1 start_ind = i * nr + 1 @@ -384,14 +385,16 @@ end function information_set_lower_bound(r::Int, n::Int, k::Int, l::Int, rank_defs::Vector{Int}, info_set_alg::Symbol; even::Bool = false, doubly_even::Bool = false, triply_even::Bool = false) + info_set_alg ∈ (:auto, :Brouwer, :Zimmermann, :White, :Chen, :Bouyuklieva, :Edmonds) || throw(ArgumentError("Unknown information set algorithm. Expected `:auto`, `:Brouwer`, `:Zimmermann`, `:White`, `:Chen`, `:Bouyuklieva`, or `:Edmonds`.")) + lower = 0 if info_set_alg == :Brouwer lower = r * length(rank_defs) elseif info_set_alg == :Zimmermann h = length(rank_defs) - lower = 0 + lower = count(x -> x != 0, rank_defs) for i in 1:h - lower += maximum([0, r - rank_defs[i]]) + lower += maximum([0, r - rank_defs[i]]) end elseif info_set_alg == :Chen lower = Int(ceil(n * r / k)) @@ -405,6 +408,9 @@ function information_set_lower_bound(r::Int, n::Int, k::Int, l::Int, rank_defs:: # elseif info_set_alg == :Edmonds # continue end + if lower == 0 + println("initial lower bound is set to 0") + end (!triply_even && !doubly_even && even) && (lower += lower % 2;) (!triply_even && doubly_even) && (lower += 4 - lower % 4;) @@ -413,19 +419,54 @@ function information_set_lower_bound(r::Int, n::Int, k::Int, l::Int, rank_defs:: end """ - Gray_code_minimum_distance(C::AbstractLinearCode; verbose::Bool = false) - + minimum_distance(C::AbstractLinearCode; alg::Symbol = :Zimmermann, v::Bool = false) Return the minimum distance of `C` using a deterministic algorithm based on enumerating constant weight codewords of the binary reflected Gray code. If a word of minimum weight is found before the lower and upper bounds cross, it is returned; otherwise, the zero vector is returned. """ -function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol = :auto, - verbose::Bool = false, dbg = Dict()) +function minimum_distance(C::AbstractLinearCode; alg::Symbol = :Zimmermann, v::Bool = false, show_progress=true) ord_F = Int(order(C.F)) ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes.")) - info_set_alg ∈ (:auto, :Brouwer, :Zimmermann, :White, :Chen, :Bouyuklieva, :Edmonds) || throw(ArgumentError("Unknown information set algorithm. Expected `:auto`, `:Brouwer`, `:Zimmermann`, `:White`, `:Chen`, `:Bouyuklieva`, or `:Edmonds`.")) + #TODO add :Bouyuklieva, and :Edmonds + + # We implement the following algorithms described in White's thesis: + # :Brouwer Algo 2.2 + # :Zimmermann Algo 2.4 + # :Chen Algo 2.6 + # :White Algo 3.1 + alg ∈ (:auto, :Brouwer, :Zimmermann, :White, :Chen) || throw(ArgumentError("Unknown information set algorithm. Expected `:auto`, `:Brouwer`, `:Zimmermann`, `:White`, `:Chen`")) + + if alg == :auto + if typeof(C) <: AbstractCyclicCode + v && println("Detected a cyclic code, using Chen's adaption.") + alg = :Chen + # TODO: fix this case + elseif typeof(C) <: AbstractQuasiCyclicCode + v && println("Detected a quasi-cyclic code, using White's adaption.") + alg = :White + else + v && println("Using Zimmermann's algorithm.") + alg = :Zimmermann + end + end + alg == :auto || throw(ErrorException("Could not determine minimum distance algo automatically")) + + if alg in (:Brouwer, :Zimmermann) || + return _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg = alg, verbose = v) + end + println("Warning: old enumeration algorithm selected. Performance will be slow") # TODO remove when all code updated + return _minimum_distance_enumeration_with_matrix_multiply(C::AbstractLinearCode; info_set_alg = alg) +end + +function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zimmermann, + verbose::Bool = false, dbg = Dict(), show_progress=true) + + ord_F = Int(order(C.F)) + ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes.")) + C.k < 2^16 || throw(DomainError("The given linear code has length k >= 2^16 which is not supported")) + info_set_alg ∈ (:Brouwer, :Zimmermann) || throw(ArgumentError("Unknown information set algorithm. Expected `:Brouwer`, `:Zimmermann`")) generator_matrix(C, true) # ensure G_stand exists if _has_empty_vec(C.G, :cols) @@ -433,33 +474,32 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol throw(ArgumentError("Codes with standard form of generator matrix having 0 columns not supported")) end G = C.G #TODO remove local var G + # generate if not pre-stored parity_check_matrix(C) - if info_set_alg == :auto - if typeof(C) <: AbstractCyclicCode - verbose && println("Detected a cyclic code, using Chen's adaption.") - info_set_alg = :Chen - # TODO: fix this case - elseif typeof(C) <: AbstractQuasiCyclicCode - verbose && println("Detected a quasi-cyclic code, using White's adaption.") - info_set_alg = :White - else - info_set_alg = :Edmonds - end - end - # should have no need to permute to find better ranks because of Edmond's? A_mats, perms_mats, rnks = information_sets(G, info_set_alg, permute = true, only_A = false) - A_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Ai)') for Ai in A_mats] - perms_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Pi)') for Pi in perms_mats] + A_mats = [deepcopy(_Flint_matrix_to_Julia_T_matrix(Ai, UInt16)') for Ai in A_mats] + perms_mats = [deepcopy(_Flint_matrix_to_Julia_T_matrix(Pi, UInt16)') for Pi in perms_mats] h = length(A_mats) # println("Starting loop to refine upper bound. Initial upper bound ", C.u_bound, " num of mats is ", length(A_mats), " dimension ", size(A_mats[1])) rank_defs = zeros(Int, h) - optimize_upper_bound = true # true will be used in production. False makes it less likely that the algo will early exit at some trivial value - dbg["exit_r"] = -1 - dbg["found_codewords"] = [] + if length(keys(dbg)) > 0 + println("Debug mode ON") + end + dbg_key_found = "found_codewords" + dbg_key_exit_r = "exit_r" + + if haskey(dbg, dbg_key_exit_r) + verbose && println("dbg Dict: Storing the largest message weight reached @key=$dbg_key_found") + dbg["exit_r"] = -1 + end + store_found_codewords = haskey(dbg, dbg_key_found) + if store_found_codewords + verbose && println("dbg Dict: Storing the codewords in the debug dictionary, key=$dbg_key_found") + end k, n = size(G) if info_set_alg == :Brouwer && rnks[h] != k @@ -490,39 +530,78 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol (!triply_even_flag && !doubly_even_flag && even_flag) && println("Detected an even code.") end - found = zeros(Int, n) - if optimize_upper_bound - # C_orig_U_bound = C.u_bound #TODO restore this before returning if its smaller than C.u_bound at the end - C.u_bound = n - end # initial_perm_ind will match the permutation we use for the 'found' vector if the found vector is nonzero. To simplify the code below we're going to choose an initial permutation arbitrarily. initial_perm_ind = 1 + # following loop is the r=1 case of the enumeration. We do this case here because we want to make a good guess at the terminating r before we start multiple threads for (j, g) in enumerate(A_mats) # loop over the A_mats rather than the original G because it would add another case to deal with later # can make this faster with dots and views + if store_found_codewords + for a in 1:size(g, 2) + output_vec = perms_mats[j] * g[:, a] + # @assert in(output_vec, C) + push!(dbg["found_codewords"], (1, [], output_vec)) + end + end + # println(dbg["found_codewords"]) w, i = _min_wt_col(g) if w <= C.u_bound found = g[:, i] C.u_bound = w - initial_perm_ind = j - y = perms_mats[initial_perm_ind] * found + y = perms_mats[j] * found end end verbose && println("Current upper bound: $(C.u_bound)") verbose && !iszero(found) && println("Found element matching upper bound.") + found = A_mats[1][:, 1] + + l = 0 + if verbose + _, _, b_rnks = information_sets(G, :Brouwer, permute = true, only_A = false) + b_h = length(b_rnks) + b_lower_bounds = [information_set_lower_bound(r+1, n, k, l, [k - 0 for i in 1:b_h], :Brouwer, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1] + b_r_term = findfirst(x -> x ≥ C.u_bound, b_lower_bounds) + + # _, _, z_rnks = information_sets(G, :Zimmermann, permute = true, only_A = false) + # z_h = length(b_rnks) + # z_lower_bounds = [information_set_lower_bound(r+1, n, k, l, [k - z_rnks[i] for i in 1:z_h], :Zimmermann, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1] + # z_r_term = findfirst(x -> x ≥ C.u_bound, z_lower_bounds) + # verbose && println("ranks: Brouwer $b_rnks Zimm $z_rnks") + # verbose && println("Predicted termination weight based on current upper bound: Brouwer $b_r_term Zimm $z_r_term") + end - info_set_alg == :Chen ? (l = C.l;) : (l = 0;) - lower_bounds = [information_set_lower_bound(r, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k] - r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds) - isnothing(r_term) && (r_term = k;) + #Note the r+1 here. + lower_bounds_for_prediction = [information_set_lower_bound(r+1, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1] + #TODO the next line is temporary while Zimmermann is broken. We use it because we want to compare our Brouwer's performance to Magma's Brouw-Zimm + r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds_for_prediction) + if isnothing(r_term) + raise(DomainError("invalid termination r")) + end verbose && println("Predicted termination weight based on current upper bound: $r_term") + #In the main loop we check if lower bound > upper bound before we enumerate and so the lower bounds for the loop use r not r+1 + lower_bounds = [information_set_lower_bound(r, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1] + if !store_found_codewords + lower_bounds = [lower_bounds[r]+(r+1) for r in 1:k-1] + end + num_thrds = Threads.nthreads() verbose && println("Detected $num_thrds threads.") power = Int(floor(log(2, num_thrds))) - store_vecs = haskey(dbg, "found_codewords") - for r in 1:k + work_factor_up_to_log_field_size = 0 + predicted_work_factor = fld(n, k) * sum([binomial(k, i) for i in 1:r_term]) + verbose && println("Predicted work factor: $predicted_work_factor") + if show_progress + prog_bar = Progress(predicted_work_factor, dt=1.0, showspeed=true) # updates no faster than once every 1s + end + weight_sum_bound = min(2 * C.u_bound + 5, n) + verbose && println("Codeword weights initially checked on first $weight_sum_bound entries") + + for r in 2:k + if r > 2^16 + verbose && println("Reached an r with r>2^16") + end C.l_bound < lower_bounds[r] && (C.l_bound = lower_bounds[r];) # an even code can't have have an odd minimum weight # (!triply_even_flag && !doubly_even_flag && even_flag) && (C.l_bound += C.l_bound % 2;) @@ -533,10 +612,12 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol verbose && println("Upper bound: $(C.u_bound)") if C.l_bound >= C.u_bound C.d = C.u_bound - y = perms_mats[initial_perm_ind] * found + y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) # @assert in(y, C) - dbg["exit_r"] = r - return C.u_bound, (C.F).(y) #TODO + dbg[dbg_key_exit_r] = r + #TODO report actual/expected work factors in dbg dict + show_progress && ProgressMeter.finish!(prog_bar) + return C.u_bound, y end if verbose @@ -555,44 +636,75 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol founds = [found for _ in 1:num_thrds] exit_thread_indicator_vec = [initial_perm_ind for _ in 1:num_thrds] - c = zeros(Int, C.n) - itrs = [GrayCode(C.k, r, digits(m - 1, base = 2, pad = power), mutate = true) for m in 1:num_thrds] + # itrs_alloc = @allocated begin + itrs = [SubsetGrayCode(C.k, r)] + # lens = [length(_) for _ in itrs] + # println(lens) + # expected_work_factor += h * sum(lens) + # end + # println("alloc itrs used ", itrs_alloc) + # Threads.@threads for itr in itrs - first_itr = true + vec = vcat(fill(1, r), fill(0, C.k - r)) # initial vec corresponds to the subset {1,..,r} for m in 1:num_thrds + #= + # c_itr is a binary vector. If it were a _random_ vector the expected + # weight of the subvector c_itr[1:2*uppers[m]] equals uppers[m]. + # So we use the heuristic that w(c_itr[1:2*uppers[m]+c]) for a small constant c is + # usually larger than uppers[m]. + # If not, we compute weight of all of c_itr + =# + weight_sum_bound = min(2 * uppers[m] + 5, n) + itr = itrs[m] # for u in itr #TODO its easier to test the iteration if the loops are in the order used by GW for i in 1:h # for i in 1:h #TODO + c_itr = zeros(UInt16, C.n) + is_first = true + curr_mat = A_mats[i] for u in itr - vec = u - if r - rank_defs[i] > 0 - LinearAlgebra.mul!(c, A_mats[i], vec) - - #TODO for testing - if store_vecs - F = Oscar.Nemo.Native.GF(2) - output_vec = F.(Vector(perms_mats[initial_perm_ind] * c)) - push!(dbg["found_codewords"], (r, deepcopy(vec), output_vec)) + work_factor_up_to_log_field_size += 1 + show_progress && ProgressMeter.next!(prog_bar) + if r - rank_defs[i] > 0 + for ind in u + + if ind != -1 + is_first = false + @simd for i in eachindex(c_itr) + @inbounds c_itr[i] = xor(c_itr[i], curr_mat[i, ind]) + end + end + end + if is_first + LinearAlgebra.mul!(c_itr, curr_mat, vec) + @inbounds @simd for j in 1:n + c_itr[j] %= p + end end - w = r - @inbounds for j in 1:n #TODO unless weve cropped the Gi down to the Ai shouldnt this be skipping the indexs corresp to the identity submatrix - c[j] % p != 0 && (w += 1;) + if store_found_codewords + output_vec = perms_mats[initial_perm_ind] * c_itr + push!(dbg["found_codewords"], (r, [], output_vec)) end - if uppers[m] > w - uppers[m] = w - founds[m] = c - exit_thread_indicator_vec[m] = i - verbose && println("Adjusting (local) upper bound: $w") - if C.l_bound == uppers[m] - Threads.atomic_cas!(flag, true, false) - break - else - r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds) - isnothing(r_term) && (r_term = k;) - verbose && println("Updated termination weight: $r_term") + partial_weight = sum(view(c_itr, 1:weight_sum_bound)) + + if uppers[m] > partial_weight + w = sum(c_itr) + if uppers[m] > w + uppers[m] = w + founds[m] = c_itr #TODO reduce founds coeffs here and below + exit_thread_indicator_vec[m] = i + verbose && println("Adjusting (local) upper bound: $w") + if C.l_bound == uppers[m] + Threads.atomic_cas!(flag, true, false) + break + else + r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds) + isnothing(r_term) && (r_term = k;) + verbose && println("Updated termination weight: $r_term") + end end end end @@ -604,20 +716,20 @@ function Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol found = founds[loc] initial_perm_ind = exit_thread_indicator_vec[loc] end + show_progress && ProgressMeter.finish!(prog_bar) # at this point we are guaranteed to have found the answer C.d = C.u_bound - y = perms_mats[initial_perm_ind] * found + y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) # iszero(C.H * transpose(y)) - return C.u_bound, (C.F).(y) + return C.u_bound, y_Oscar end -function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::Symbol = :auto, - verbose::Bool = false, dbg = Dict(), show_progress=true) +function _minimum_distance_enumeration_with_matrix_multiply(C::AbstractLinearCode; info_set_alg::Symbol = :auto, + verbose::Bool = false, dbg = Dict()) ord_F = Int(order(C.F)) ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes.")) - C.k < 2^16 || throw(DomainError("The given linear code has length k >= 2^16 which is not supported")) info_set_alg ∈ (:auto, :Brouwer, :Zimmermann, :White, :Chen, :Bouyuklieva, :Edmonds) || throw(ArgumentError("Unknown information set algorithm. Expected `:auto`, `:Brouwer`, `:Zimmermann`, `:White`, `:Chen`, `:Bouyuklieva`, or `:Edmonds`.")) generator_matrix(C, true) # ensure G_stand exists @@ -625,8 +737,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S #TODO err string can instruct the user to construct a new code without 0 cols and tell them the function for that throw(ArgumentError("Codes with standard form of generator matrix having 0 columns not supported")) end - G = C.G #TODO remove local var G - + G = C.G # generate if not pre-stored parity_check_matrix(C) @@ -644,39 +755,12 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S end # should have no need to permute to find better ranks because of Edmond's? A_mats, perms_mats, rnks = information_sets(G, info_set_alg, permute = true, only_A = false) - - A_mats = [deepcopy(_Flint_matrix_to_Julia_T_matrix(Ai, UInt16)') for Ai in A_mats] - perms_mats = [deepcopy(_Flint_matrix_to_Julia_T_matrix(Pi, UInt16)') for Pi in perms_mats] + A_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Ai)') for Ai in A_mats] + perms_mats = [deepcopy(_Flint_matrix_to_Julia_int_matrix(Pi)') for Pi in perms_mats] h = length(A_mats) - # println("Starting loop to refine upper bound. Initial upper bound ", C.u_bound, " num of mats is ", length(A_mats), " dimension ", size(A_mats[1])) rank_defs = zeros(Int, h) - if length(keys(dbg)) > 0 - println("Debug mode ON") - end - dbg_key_unoptimized_uppr_bnd = "unoptimized_uppr_bnd" - dbg_key_found = "found_codewords" - dbg_key_exit_r = "exit_r" - - unoptimized_uppr_bnd = false - if haskey(dbg, dbg_key_unoptimized_uppr_bnd) - verbose && println("dbg Dict: Using unoptimized upper bound") - unoptimized_uppr_bnd = true - end - if haskey(dbg, dbg_key_exit_r) - verbose && println("dbg Dict: Storing the largest message weight reached @key=$dbg_key_found") - dbg["exit_r"] = -1 - end - store_found_codewords = haskey(dbg, dbg_key_found) - if store_found_codewords - verbose && println("dbg Dict: Storing the codewords in the debug dictionary, key=$dbg_key_found") - end - k, n = size(G) - if info_set_alg == :Brouwer && rnks[h] != k - println("Rank of last matrix too small") - return - end if verbose print("Generated $h information sets with ranks: ") for i in 1:h @@ -701,13 +785,9 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S (!triply_even_flag && !doubly_even_flag && even_flag) && println("Detected an even code.") end - found = zeros(Int, n) - if unoptimized_uppr_bnd - # C_orig_U_bound = C.u_bound #TODO restore this before returning if its smaller than C.u_bound at the end - C.u_bound = n - end # initial_perm_ind will match the permutation we use for the 'found' vector if the found vector is nonzero. To simplify the code below we're going to choose an initial permutation arbitrarily. initial_perm_ind = 1 + found = A_mats[1][:, 1] for (j, g) in enumerate(A_mats) # loop over the A_mats rather than the original G because it would add another case to deal with later # can make this faster with dots and views w, i = _min_wt_col(g) @@ -732,17 +812,7 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S verbose && println("Detected $num_thrds threads.") power = Int(floor(log(2, num_thrds))) - work_factor_up_to_log_field_size = 0 - predicted_work_factor = fld(n, k) * sum([binomial(k, i) for i in 1:r_term]) - verbose && println("Predicted work factor: $predicted_work_factor") - if show_progress - prog_bar = Progress(predicted_work_factor, dt=1.0, showspeed=true) # updates no faster than once every 1s - end - for r in 1:k - if r > 2^16 - verbose && println("Reached an r with r>2^16") - end C.l_bound < lower_bounds[r] && (C.l_bound = lower_bounds[r];) # an even code can't have have an odd minimum weight # (!triply_even_flag && !doubly_even_flag && even_flag) && (C.l_bound += C.l_bound % 2;) @@ -753,12 +823,9 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S verbose && println("Upper bound: $(C.u_bound)") if C.l_bound >= C.u_bound C.d = C.u_bound - y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) + y = perms_mats[initial_perm_ind] * found # @assert in(y, C) - dbg[dbg_key_exit_r] = r - #TODO report actual/expected work factors in dbg dict - show_progress && ProgressMeter.finish!(prog_bar) - return C.u_bound, y + return C.u_bound, (C.F).(y) end if verbose @@ -770,70 +837,32 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S end flag = Threads.Atomic{Bool}(true) - num_thrds = 1 + verbose && println("Detected $num_thrds threads.") p = Int(characteristic(C.F)) uppers = [C.u_bound for _ in 1:num_thrds] founds = [found for _ in 1:num_thrds] exit_thread_indicator_vec = [initial_perm_ind for _ in 1:num_thrds] - # itrs_alloc = @allocated begin - itrs = [SubsetGrayCode(C.k, r)] - # lens = [length(_) for _ in itrs] - # println(lens) - # expected_work_factor += h * sum(lens) - # end - # println("alloc itrs used ", itrs_alloc) - - # Threads.@threads for itr in itrs - vec = vcat(fill(1, r), fill(0, C.k - r)) # initial vec corresponds to the subset {1,..,r} - for m in 1:num_thrds - # if c_itr is a binary vector. If it were a random binary vector the expected - # weight of the subvector c_itr[1:2*uppers[m]] equals uppers[m]. - # So we use the heuristic that w(c_itr[1:2*uppers[m]+5]) is usually larger than uppers[m]. - # If not, we compute weight of all of c_itr - weight_sum_bound = min(2 * uppers[m] + 5, n) - - verbose && println("Computing weight up to $weight_sum_bound") + c = zeros(Int, C.n) + itrs = [GrayCode(C.k, r, digits(m - 1, base = 2, pad = power), mutate = true) for m in 1:num_thrds] + Threads.@threads for m in 1:num_thrds itr = itrs[m] - # for u in itr #TODO its easier to test the iteration if the loops are in the order used by GW - for i in 1:h - # for i in 1:h #TODO - c_itr = zeros(UInt16, C.n) - is_first = true - curr_mat = A_mats[i] - for u in itr - work_factor_up_to_log_field_size += 1 - show_progress && ProgressMeter.next!(prog_bar) - if r - rank_defs[i] > 0 - for ind in u - - if ind != -1 - is_first = false - @simd for i in eachindex(c_itr) - @inbounds c_itr[i] = xor(c_itr[i], curr_mat[i, ind]) - end - end - end - if is_first - LinearAlgebra.mul!(c_itr, curr_mat, vec) - @inbounds @simd for j in 1:n - c_itr[j] %= p + for u in itr + if flag[] + for i in 1:h + vec = u + if r - rank_defs[i] > 0 + LinearAlgebra.mul!(c, A_mats[i], vec) + + w = r + @inbounds for j in 1:n + c[j] % p != 0 && (w += 1;) end - end - - if store_found_codewords - output_vec = perms_mats[initial_perm_ind] * c_itr - push!(dbg["found_codewords"], (r, [], output_vec)) - end - - partial_weight = sum(view(c_itr, 1:weight_sum_bound)) - if uppers[m] > partial_weight - w = sum(c_itr) - if uppers[m] > w + if uppers[m] > w uppers[m] = w - founds[m] = c_itr #TODO reduce founds coeffs here and below + founds[m] .= c exit_thread_indicator_vec[m] = i verbose && println("Adjusting (local) upper bound: $w") if C.l_bound == uppers[m] @@ -847,6 +876,8 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S end end end + else + break end end end @@ -855,14 +886,14 @@ function david_Gray_code_minimum_distance(C::AbstractLinearCode; info_set_alg::S found = founds[loc] initial_perm_ind = exit_thread_indicator_vec[loc] end - show_progress && ProgressMeter.finish!(prog_bar) # at this point we are guaranteed to have found the answer C.d = C.u_bound - y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) + y = perms_mats[initial_perm_ind] * found # iszero(C.H * transpose(y)) - return C.u_bound, y_Oscar + return C.u_bound, (C.F).(y) end + """ words_of_weight(C::AbstractLinearCode, l_bound::Int, u_bound::Int; verbose::Bool = false) @@ -1650,10 +1681,10 @@ function minimum_distance(C::AbstractLinearCode; alg::Symbol = :trellis, sect::B for i in 1:length(HWE.polynomial)])) return C.d else - return david_Gray_code_minimum_distance(C, verbose = verbose) + return _minimum_distance_BZ(C, verbose = verbose) end elseif alg == :Gray - return david_Gray_code_minimum_distance(C, verbose = verbose) + return _minimum_distance_BZ(C, verbose = verbose) elseif alg == :trellis weight_enumerator_classical(syndrome_trellis(C, "primal", false), type = type) return C.d diff --git a/src/CodingTheory.jl b/src/CodingTheory.jl index 1a47b51..a1be664 100644 --- a/src/CodingTheory.jl +++ b/src/CodingTheory.jl @@ -378,7 +378,7 @@ export Trellis, vertices, edges, isisomorphic, isequal, loadbalancedecode, include("Classical/weight_dist.jl") export polynomial, type, CWE_to_HWE, weight_enumerator, MacWilliams_identity, weight_distribution, weight_plot, support, minimum_distance, Sterns_attack, - david_Gray_code_minimum_distance, minimum_words, words_of_weight + _minimum_distance_BZ, minimum_words, words_of_weight ############################# # Quantum/weight_dist.jl From 234113cf56cea2cfcc8f37848089208ca1fd992b Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Mon, 18 Nov 2024 11:47:23 -0500 Subject: [PATCH 09/23] Add files for benchmarks --- benchmarks/iterators_bench.jl | 50 ++++++++++++++++++++ benchmarks/weight_dist_bench.jl | 81 +++++++++++++++++++++++++++++++++ 2 files changed, 131 insertions(+) create mode 100644 benchmarks/iterators_bench.jl create mode 100644 benchmarks/weight_dist_bench.jl diff --git a/benchmarks/iterators_bench.jl b/benchmarks/iterators_bench.jl new file mode 100644 index 0000000..39bab49 --- /dev/null +++ b/benchmarks/iterators_bench.jl @@ -0,0 +1,50 @@ +using CodingTheory +using Oscar +using LinearAlgebra +using BenchmarkTools + +function visit_all_subsets(mutate, g_len, v_weight, subset) + if subset + gi=CodingTheory.SubsetGrayCode(g_len, v_weight) + else + gi=CodingTheory.GrayCode(g_len, v_weight; mutate = mutate) + end + for subset in gi + # BenchmarkTools can elide simple computations so we need to do some nontrivial calculation here. + # In any realistic situation we need to look at least look at all entries of the vector. + # Ill model the task of looking at the entries using the 'all' function + Base.all(i->(i==0), subset) + end + return +end + +#= +SubsetGrayCode and GrayCode appear to be in the same ballpark. +GrayCode is faster, SubsetGrayCode uses less memory. + +###### +# n=25, k=12, mutate=true +###### +@benchmark visit_all_subsets(25, 12, false) +BenchmarkTools.Trial: 54 samples with 1 evaluation. + Range (min … max): 66.365 ms … 219.557 ms ┊ GC (min … max): 0.00% … 34.15% + Time (median): 92.912 ms ┊ GC (median): 26.24% + Time (mean ± σ): 93.526 ms ± 27.981 ms ┊ GC (mean ± σ): 19.78% ± 13.74% + + ▂▂ █▃ + ███▅▅▁▁▁▁▄▁▁█████▁▁▄▁▁▁▁▄▁▁▁▁▁▅▁▁▁▁▁▁▁▁▁▄▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄ ▁ + 66.4 ms Histogram: frequency by time 181 ms < + + Memory estimate: 158.70 MiB, allocs estimate: 5200322. + +BenchmarkTools.Trial: 65 samples with 1 evaluation. + Range (min … max): 51.159 ms … 103.931 ms ┊ GC (min … max): 0.00% … 50.88% + Time (median): 77.176 ms ┊ GC (median): 31.98% + Time (mean ± σ): 77.271 ms ± 6.824 ms ┊ GC (mean ± σ): 31.52% ± 7.35% + + █ ▁ + ▃▃▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▆████▅▆▅▁▃▃▁▁▄▁▁▃▁▁▁▁▃ ▁ + 51.2 ms Histogram: frequency by time 89.3 ms < + + Memory estimate: 238.05 MiB, allocs estimate: 5200308. +=# \ No newline at end of file diff --git a/benchmarks/weight_dist_bench.jl b/benchmarks/weight_dist_bench.jl new file mode 100644 index 0000000..76b7fb2 --- /dev/null +++ b/benchmarks/weight_dist_bench.jl @@ -0,0 +1,81 @@ +using Oscar, CodingTheory +using Combinatorics +using LinearAlgebra +using BenchmarkTools +using Random +using Dates +# using AllocCheck +# using Profile +CodingTheory.Random.seed!(0) + +function white_n70_k35(use_mine, terse_description) + F = Oscar.Nemo.Native.GF(2) + mat = matrix(F, [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 0 0 1 1 0 1 1 1 0 1 0 1 1 1 1 0 1 1 0 1 1 0 0 1 0 0 1 1 0 1; + 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 1 1 1 0 1 1 0 1 0 0 1 1 0 1 0 1 1 0 1 1 0 1 1 1 0 0 0 1 0; + 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 1 0 0 1 1 0 0 0 0 0 1 0 0 1 1 1 0 1 1 0 1 0 0 1 1 1 1 0 0; + 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 1 0 0 0 0 0 0 1 0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 0 0 1; + 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 1 0 0 0 0 1 0 0 0 1 0 1 1 0 1 0 0 0 1 1 0 1 0 0 0 0 0 0 1 0; + 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 1 1 0 0 0 0 1 1 0 1 0 1 0 1 0 1 0 0 0 1 1 0 0 1 0 0 1 0 1 1; + 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 1 1 0 0 1 0 0 1 1 1 1 1 0 0 1 1 1 0 1 1 1 0 0 0 0 1; + 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1 1 1 1 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 1 1 1 1 1 1 0; + 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 1 0 1 1 1 0 1 0 0 0 1 1 1 1 0 0 1 0 0 0 1 0 1 1 1 1 1 0 1 0 1; + 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 0 1 0 1 1 1 0 0 0 0 1 0 0 1 0 0 1 1 1 1 0 0 1 1 1 1 1 0; + 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 0 1 0 1 1 1 0 0 0 0 1 0 0 1 0 0 1 1 1 1 0 0 1 1 1 1 1; + 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 1 1 1 0 1 1 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 1 0 0; + 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 1 0 1 1 0 0 1 0 0 1 0 0 0 0 0 1 0 1 1; + 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 1 1 0 0 1 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 0 1 1 0 0 0 0 0 1; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 1 0 1 0 0 1 1 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 1 1; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 1 1 0 0 0 1 0 0 0 1 1 1 1 1 0 1 0 0 1 1 0 0 1 1 1 1 1; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 1 1 0 1 0 1 1 0 0 0 1 0 1 1 0 0 1 0 1 1 1 1 0 1 0 0 0 1 1 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 1 0 1 0 0 1 1 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 1 1 0 1 1 0 0 1 1 1 0 0 0 1 1 0 1 0 0 1 1 0 1 1 0 1 1 1 1 1; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 1 1 0 1 1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 1 1 0 0 1 1 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 1 1 0 1 0 1 1 0 0 0 1 0 1 1 0 0 1 0 1 1 1 1 0 1 0 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 0 0 0 0 0 0 0 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 0 1 1 0 1 1 1; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 1 1 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 0 0 0 0 0 1 0 0 1 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 1 1 0 1 1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 1 1 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 0 0 1 1 1 0 1 0 1 1 1 0 1 0 1 0 1 0 0 1 1 0 0 0 0 0 1 0 1 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 1 1 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 0 0 0 0 0 1 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 1 1 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 0 0 0 0 0 1; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 1 1 0 0 0 1 1 1 1 0 0 0 1 0 0 0 1 0 0 0 1 1; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 0 1 0 1 1 1 0 0 1 1 1 1 1 1 0 0 0 0 1 0 1 0 1 0 0 1 0 1 0 0 1 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 1 0 1 0 0 0 1 1 1 0 1 0 1 0 0 1 1 0 1 0 0 0 1 1 0 1 1 0 0 1 0 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 1 1 1 1 0 0 0 1 1 0 0 0 1 0 1 0 1 1 0 1 0 0 0 1 0 1 1 1 1 0 0 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 1 1 0 1 1 0 0 0 1 0 0 1 1 0 1 1 0 1 1 0 1 0 0 0 0 1 1 1 0 1 1 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 1 0 0 0 0 0 1 1 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 1 1 1 1 0 0; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 1 1 1 0 1 1 0 1 0 0 1 1 0 1 0 1 1 0 1 1 0 1 1 1 0 0 0 1 0 0 1 1; + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 1 1 0 1 1 0 0 0 1 0 0 1 1 0 1 1 0 1 1 0 1 0 0 0 0 1 1 1]) + + num_thrds = Threads.nthreads() + println("Benchmarking with $num_thrds threads.") + if use_mine + function_description = "_minimum_distance_BZ(codeN70K35,:Brouwer;false)" + terse_description ? println(function_description) : println(benchmark_description(function_description)) + return @benchmark CodingTheory._minimum_distance_BZ(codeN70K35, info_set_alg = :Brouwer; verbose = false) setup = (codeN70K35=LinearCode($mat, false, false)) + else + function_description = "_minimum_distance_enumeration_with_matrix_multiply(codeN70K35,:Brouwer;false)" + terse_description ? println(function_description) : println(benchmark_description(function_description)) + return @benchmark CodingTheory._minimum_distance_enumeration_with_matrix_multiply(codeN70K35, info_set_alg = :Brouwer; verbose = false) setup = (codeN70K35=LinearCode($mat, false, false)) + end +end + +function benchmark_description(function_description) + return "@" * Dates.format(now(), dateformat"MM:SS") * " benchmark " * function_description * + "\nMost recent git commit\n" * current_git_hash() * "\ndate: " * Dates.format(now(), dateformat"y-m-d : H:M:S") +end + +function current_git_hash() + try + # Run the Git command to get the latest commit hash + hash = readchomp(`git rev-parse --short HEAD`) + return hash + catch + # Return a default value if not in a Git repository or Git fails + return "no-git-repo" + end +end + +use_new=false +t = white_n70_k35(use_new, true) +# t = test_white_105(use_mine) +t \ No newline at end of file From e0593adb18f5d4bad93df1a6043191b241d26877 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Mon, 18 Nov 2024 12:51:57 -0500 Subject: [PATCH 10/23] Add utils file --- benchmarks/weight_dist_bench.jl | 19 ++----------------- 1 file changed, 2 insertions(+), 17 deletions(-) diff --git a/benchmarks/weight_dist_bench.jl b/benchmarks/weight_dist_bench.jl index 76b7fb2..61c2c34 100644 --- a/benchmarks/weight_dist_bench.jl +++ b/benchmarks/weight_dist_bench.jl @@ -6,6 +6,7 @@ using Random using Dates # using AllocCheck # using Profile +include("benchmark_utils.jl") CodingTheory.Random.seed!(0) function white_n70_k35(use_mine, terse_description) @@ -59,23 +60,7 @@ function white_n70_k35(use_mine, terse_description) end end -function benchmark_description(function_description) - return "@" * Dates.format(now(), dateformat"MM:SS") * " benchmark " * function_description * - "\nMost recent git commit\n" * current_git_hash() * "\ndate: " * Dates.format(now(), dateformat"y-m-d : H:M:S") -end - -function current_git_hash() - try - # Run the Git command to get the latest commit hash - hash = readchomp(`git rev-parse --short HEAD`) - return hash - catch - # Return a default value if not in a Git repository or Git fails - return "no-git-repo" - end -end - -use_new=false +use_new=true t = white_n70_k35(use_new, true) # t = test_white_105(use_mine) t \ No newline at end of file From 040cf7a12dd6145077809bd285d8dbf9917b28a4 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Mon, 18 Nov 2024 14:39:07 -0500 Subject: [PATCH 11/23] Remove old test code that made the upper/lower crossover too small --- benchmarks/weight_dist_bench.jl | 2 +- src/Classical/weight_dist.jl | 19 +++++++------------ 2 files changed, 8 insertions(+), 13 deletions(-) diff --git a/benchmarks/weight_dist_bench.jl b/benchmarks/weight_dist_bench.jl index 61c2c34..ba5759a 100644 --- a/benchmarks/weight_dist_bench.jl +++ b/benchmarks/weight_dist_bench.jl @@ -60,7 +60,7 @@ function white_n70_k35(use_mine, terse_description) end end -use_new=true +use_new=false t = white_n70_k35(use_new, true) # t = test_white_105(use_mine) t \ No newline at end of file diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 8e613bd..63906ad 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -572,7 +572,6 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim #Note the r+1 here. lower_bounds_for_prediction = [information_set_lower_bound(r+1, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1] - #TODO the next line is temporary while Zimmermann is broken. We use it because we want to compare our Brouwer's performance to Magma's Brouw-Zimm r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds_for_prediction) if isnothing(r_term) raise(DomainError("invalid termination r")) @@ -581,9 +580,6 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim #In the main loop we check if lower bound > upper bound before we enumerate and so the lower bounds for the loop use r not r+1 lower_bounds = [information_set_lower_bound(r, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1] - if !store_found_codewords - lower_bounds = [lower_bounds[r]+(r+1) for r in 1:k-1] - end num_thrds = Threads.nthreads() verbose && println("Detected $num_thrds threads.") @@ -595,6 +591,13 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim if show_progress prog_bar = Progress(predicted_work_factor, dt=1.0, showspeed=true) # updates no faster than once every 1s end + #= TODO rewrite + # c_itr is a binary vector. If it were a _random_ vector the expected + # weight of the subvector c_itr[1:2*uppers[m]] equals uppers[m]. + # So we use the heuristic that w(c_itr[1:2*uppers[m]+c]) for a small constant c is + # usually larger than uppers[m]. + # If not, we compute weight of all of c_itr + =# weight_sum_bound = min(2 * C.u_bound + 5, n) verbose && println("Codeword weights initially checked on first $weight_sum_bound entries") @@ -647,14 +650,6 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim # Threads.@threads for itr in itrs vec = vcat(fill(1, r), fill(0, C.k - r)) # initial vec corresponds to the subset {1,..,r} for m in 1:num_thrds - #= - # c_itr is a binary vector. If it were a _random_ vector the expected - # weight of the subvector c_itr[1:2*uppers[m]] equals uppers[m]. - # So we use the heuristic that w(c_itr[1:2*uppers[m]+c]) for a small constant c is - # usually larger than uppers[m]. - # If not, we compute weight of all of c_itr - =# - weight_sum_bound = min(2 * uppers[m] + 5, n) itr = itrs[m] # for u in itr #TODO its easier to test the iteration if the loops are in the order used by GW From 515b402170cc0422b07a1af223cd48f74a388acd Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Tue, 19 Nov 2024 10:20:31 -0500 Subject: [PATCH 12/23] Tweak flag --- src/Classical/weight_dist.jl | 15 +++++++++------ src/iterators.jl | 10 ---------- 2 files changed, 9 insertions(+), 16 deletions(-) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 63906ad..1ab715e 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -631,7 +631,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim count > 0 && println("$count of the original $h information sets no longer contribute to the lower bound") end - flag = Threads.Atomic{Bool}(true) + keep_going = Threads.Atomic{Bool}(true) num_thrds = 1 p = Int(characteristic(C.F)) @@ -652,13 +652,15 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim for m in 1:num_thrds itr = itrs[m] - # for u in itr #TODO its easier to test the iteration if the loops are in the order used by GW + # we loop over matrices first so that we can update the lower bound more often (see White Algo 7.1) for i in 1:h - # for i in 1:h #TODO c_itr = zeros(UInt16, C.n) is_first = true curr_mat = A_mats[i] for u in itr + if !keep_going + break + end work_factor_up_to_log_field_size += 1 show_progress && ProgressMeter.next!(prog_bar) if r - rank_defs[i] > 0 @@ -689,12 +691,11 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim w = sum(c_itr) if uppers[m] > w uppers[m] = w - founds[m] = c_itr #TODO reduce founds coeffs here and below + founds[m] = c_itr exit_thread_indicator_vec[m] = i verbose && println("Adjusting (local) upper bound: $w") if C.l_bound == uppers[m] - Threads.atomic_cas!(flag, true, false) - break + Threads.atomic_cas!(keep_going, true, false) else r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds) isnothing(r_term) && (r_term = k;) @@ -720,6 +721,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim return C.u_bound, y_Oscar end +#= function _minimum_distance_enumeration_with_matrix_multiply(C::AbstractLinearCode; info_set_alg::Symbol = :auto, verbose::Bool = false, dbg = Dict()) @@ -888,6 +890,7 @@ function _minimum_distance_enumeration_with_matrix_multiply(C::AbstractLinearCod # iszero(C.H * transpose(y)) return C.u_bound, (C.F).(y) end +=# """ words_of_weight(C::AbstractLinearCode, l_bound::Int, u_bound::Int; verbose::Bool = false) diff --git a/src/iterators.jl b/src/iterators.jl index 0b4dfc7..3281a58 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -122,16 +122,6 @@ end return (inds, (v, rank, inds)) end -@inline function rest(G::SubsetGrayCode, rank) - #TODO there will be a future PR for integrating this (its used for multithreaded weight calculation) - #= - _subset_unrank(rank, G.n, vec) - inds = Vector{Int}([-1,-1,-1]) - state = (vec, rank, inds) - return Base.rest(G, state) - =# -end - function _update_indices!(indices::Vector{Int}, x::Int, y::Int) _update_indices!(indices, x) _update_indices!(indices, y) From 70c2c4014812a5f2959b3d3cf9af3e739281b1c5 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Tue, 19 Nov 2024 16:16:15 -0500 Subject: [PATCH 13/23] Add support for multiple threads to iterator --- src/Classical/linear_code.jl | 2 +- src/Classical/weight_dist.jl | 10 ++++---- src/CodingTheory.jl | 4 ++-- src/iterators.jl | 41 ++++++++++++++++++++++++++------- src/utils.jl | 6 ++--- test/iterators_test.jl | 44 +++++++++++++++++++++++++++++++++--- 6 files changed, 85 insertions(+), 22 deletions(-) diff --git a/src/Classical/linear_code.jl b/src/Classical/linear_code.jl index a0e061e..1c1014f 100644 --- a/src/Classical/linear_code.jl +++ b/src/Classical/linear_code.jl @@ -344,7 +344,7 @@ minimum_distance_upper_bound(C::AbstractLinearCode) = C.u_bound Return `true` if code is maximum distance separable (MDS). """ function is_MDS(C::AbstractLinearCode) - ismissing(C.d) && minimum_distance(C) + ismissing(C.d) && minimum_distance_Gray(C) return C.d == Singleton_bound(C.n, C.k) end diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 1ab715e..e092358 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -419,13 +419,13 @@ function information_set_lower_bound(r::Int, n::Int, k::Int, l::Int, rank_defs:: end """ - minimum_distance(C::AbstractLinearCode; alg::Symbol = :Zimmermann, v::Bool = false) + minimum_distance_Gray(C::AbstractLinearCode; alg::Symbol = :Zimmermann, v::Bool = false) Return the minimum distance of `C` using a deterministic algorithm based on enumerating constant weight codewords of the binary reflected Gray code. If a word of minimum weight is found before the lower and upper bounds cross, it is returned; otherwise, the zero vector is returned. """ -function minimum_distance(C::AbstractLinearCode; alg::Symbol = :Zimmermann, v::Bool = false, show_progress=true) +function minimum_distance_Gray(C::AbstractLinearCode; alg::Symbol = :auto, v::Bool = false, show_progress=true) ord_F = Int(order(C.F)) ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes.")) @@ -454,7 +454,7 @@ function minimum_distance(C::AbstractLinearCode; alg::Symbol = :Zimmermann, v::B alg == :auto || throw(ErrorException("Could not determine minimum distance algo automatically")) if alg in (:Brouwer, :Zimmermann) || - return _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg = alg, verbose = v) + return _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg = alg, verbose = v, show_progress = show_progress) end println("Warning: old enumeration algorithm selected. Performance will be slow") # TODO remove when all code updated return _minimum_distance_enumeration_with_matrix_multiply(C::AbstractLinearCode; info_set_alg = alg) @@ -1679,10 +1679,10 @@ function minimum_distance(C::AbstractLinearCode; alg::Symbol = :trellis, sect::B for i in 1:length(HWE.polynomial)])) return C.d else - return _minimum_distance_BZ(C, verbose = verbose) + return minimum_distance_Gray(C, verbose = verbose) end elseif alg == :Gray - return _minimum_distance_BZ(C, verbose = verbose) + return minimum_distance_Gray(C, verbose = verbose) elseif alg == :trellis weight_enumerator_classical(syndrome_trellis(C, "primal", false), type = type) return C.d diff --git a/src/CodingTheory.jl b/src/CodingTheory.jl index 2455da6..0f21a34 100644 --- a/src/CodingTheory.jl +++ b/src/CodingTheory.jl @@ -377,8 +377,8 @@ export Trellis, vertices, edges, isisomorphic, isequal, loadbalancedecode, include("Classical/weight_dist.jl") export polynomial, type, CWE_to_HWE, weight_enumerator, MacWilliams_identity, - weight_distribution, weight_plot, support, minimum_distance, Sterns_attack, - _minimum_distance_BZ, minimum_words, words_of_weight + weight_distribution, weight_plot, support, minimum_distance_Gray, minimum_distance, + Sterns_attack, minimum_words, words_of_weight ############################# # Quantum/weight_dist.jl diff --git a/src/iterators.jl b/src/iterators.jl index 3281a58..1f92934 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -5,18 +5,43 @@ # LICENSE file in the root directory of this source tree. struct SubsetGrayCode + # details about the iterator using this struct are in the comments in the iterate() function below n::Int # length of codewords k::Int # weight of codewords len::Int # length of the iterator + init_rank::BigInt # iteration will start with the subset of this rank end function SubsetGrayCode(n::Int, k::Int) - if 0 <= k <= n - len = factorial(big(n)) ÷ (factorial(big(k)) * factorial(big(n - k))) - else - len = 0 + init_rank = big(1) + num_threads = 1 + return SubsetGrayCodeFromNumThreads(n, k, init_rank, num_threads) +end + +function SubsetGrayCodeFromNumThreads(n::Int, k::Int, init_rank::BigInt, num_threads::Int) + bin = extended_binomial(n, k) + if bin % num_threads != 0 + throw(ArgumentError("num_threads must divide binomial(n k)")) + end + len = fld(bin, num_threads) + return SubsetGrayCode(n, k, len, init_rank) +end + +function SubsetGrayCodesFromNumThreads(n::Int, k::Int, num_threads::Int) + bin = extended_binomial(n, k) + if bin % num_threads != 0 + throw(ArgumentError("num_threads must divide binomial(n k)")) + end + len = fld(bin, num_threads) + + itrs = fill(SubsetGrayCode(1,1,1,1), num_threads) + for i in 0:num_threads-1 + init_rank = 1 + i * len + + itrs[i+1] = SubsetGrayCode(n, k, len, init_rank) + @assert length(itrs[i+1]) == len end - return SubsetGrayCode(n, k, len) + return itrs end Base.IteratorEltype(::SubsetGrayCode) = Base.HasEltype() @@ -53,10 +78,10 @@ end Note that inds is part of the iterator's state inds only to prevent reallocation in each iteration. =# - rank = Int(1) - v = collect(1:G.k) + subset_vec = zeros(UInt, G.k) + CodingTheory._subset_unrank(G.init_rank, UInt(G.n), subset_vec) inds = fill(-1, 3) - (inds, (v, rank, inds)) + (inds, (subset_vec, G.init_rank, inds)) end @inline function Base.iterate(G::SubsetGrayCode, state) diff --git a/src/utils.jl b/src/utils.jl index baa9e06..5c91c38 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1993,10 +1993,10 @@ function _rand_invertible_matrix(F::CTFieldTypes, n::Integer) return A end -function extended_binomial(x::UInt, y::UInt) - z = UInt(0) +function extended_binomial(x::Union{Int, UInt}, y::Union{Int, UInt}) + z = big(0) if y <= x - z = binomial(big(x), big(y)) + z = binomial(big(x), big(y)) end return z end diff --git a/test/iterators_test.jl b/test/iterators_test.jl index 02ffa22..ec3566b 100644 --- a/test/iterators_test.jl +++ b/test/iterators_test.jl @@ -1,5 +1,5 @@ @testitem "iterators.jl" begin - @testset "SubsetGrayCode iterates over all weight k subsets of {1, .., n}" begin + @testset "SubsetGrayCode iterates over all weight k subsets of {1, .., n} (nthreads == 1)" begin len = 15 weight = 7 sgc = CodingTheory.SubsetGrayCode(len, weight) @@ -8,14 +8,52 @@ state = (subset, 1, fill(-1, 3)) for i in 1:length(sgc) all_subsets_gray[i] = deepcopy(subset) + next = iterate(sgc, state) + isnothing(next) && break + (_, state) = next + (subset, _, _) = state + end + sort!(all_subsets_gray) + all_subsets_hecke = CodingTheory.Oscar.subsets(collect(1:len), weight) + sort!(all_subsets_hecke) + @test length(all_subsets_gray) == 6435 + @test all_subsets_gray == all_subsets_hecke + end + + @testset "SubsetGrayCode iterates over all weight k subsets of {1, .., n} (nthreads > 1)" begin + len = 15 + weight = 7 + num_threads = 3 + itrs = CodingTheory.SubsetGrayCodesFromNumThreads(len, weight, num_threads) + @test length(itrs) == num_threads + initial_subset_vecs = fill(fill(0, weight), length(itrs)) + bin = extended_binomial(len, weight) + for j in eachindex(itrs) + itr = itrs[j] + @test length(itr) == fld(bin, num_threads) + subset_vec = zeros(UInt, itr.k) + CodingTheory._subset_unrank!(itr.init_rank, UInt(itr.n), subset_vec) + initial_subset_vecs[j] = subset_vec + end + + all_subsets_gray = fill(fill(0, weight), bin) + + for j in 0:length(itrs)-1 + sgc = itrs[j + 1] + subset_vec = initial_subset_vecs[j + 1] + state = (subset_vec, 1, fill(-1, 3)) + for i in 1:length(sgc) + all_subsets_gray[j * fld(bin, num_threads) + i] = deepcopy(subset_vec) next = iterate(sgc, state) - next === nothing && break + isnothing(next) && break (_, state) = next - (subset, _, _) = state + (subset_vec, _, _) = state + end end sort!(all_subsets_gray) all_subsets_hecke = CodingTheory.Oscar.subsets(collect(1:len), weight) sort!(all_subsets_hecke) + @test length(all_subsets_gray) == 6435 @test all_subsets_gray == all_subsets_hecke end From fb72608b4a16df12b57f5499d1b330d6d689bcdc Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Wed, 20 Nov 2024 14:08:31 -0500 Subject: [PATCH 14/23] Use index instead of rank --- src/iterators.jl | 29 +++++++++++++++-------------- test/iterators_test.jl | 16 ++++++++-------- 2 files changed, 23 insertions(+), 22 deletions(-) diff --git a/src/iterators.jl b/src/iterators.jl index 1f92934..324d5c8 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -8,17 +8,17 @@ struct SubsetGrayCode # details about the iterator using this struct are in the comments in the iterate() function below n::Int # length of codewords k::Int # weight of codewords - len::Int # length of the iterator + len::BigInt # length of the iterator init_rank::BigInt # iteration will start with the subset of this rank end function SubsetGrayCode(n::Int, k::Int) init_rank = big(1) num_threads = 1 - return SubsetGrayCodeFromNumThreads(n, k, init_rank, num_threads) + return _subset_gray_code_from_num_threads(n, k, init_rank, num_threads) end -function SubsetGrayCodeFromNumThreads(n::Int, k::Int, init_rank::BigInt, num_threads::Int) +function _subset_gray_code_from_num_threads(n::Int, k::Int, init_rank::BigInt, num_threads::Int) bin = extended_binomial(n, k) if bin % num_threads != 0 throw(ArgumentError("num_threads must divide binomial(n k)")) @@ -27,7 +27,9 @@ function SubsetGrayCodeFromNumThreads(n::Int, k::Int, init_rank::BigInt, num_thr return SubsetGrayCode(n, k, len, init_rank) end -function SubsetGrayCodesFromNumThreads(n::Int, k::Int, num_threads::Int) +function _subset_gray_codes_from_num_threads(n::Int, k::Int, num_threads::Int) + # This function splits a single iterator into several pieces. + # The intended usage is to do the iteration with multiple threads bin = extended_binomial(n, k) if bin % num_threads != 0 throw(ArgumentError("num_threads must divide binomial(n k)")) @@ -37,9 +39,7 @@ function SubsetGrayCodesFromNumThreads(n::Int, k::Int, num_threads::Int) itrs = fill(SubsetGrayCode(1,1,1,1), num_threads) for i in 0:num_threads-1 init_rank = 1 + i * len - itrs[i+1] = SubsetGrayCode(n, k, len, init_rank) - @assert length(itrs[i+1]) == len end return itrs end @@ -54,6 +54,7 @@ Base.in(v::Vector{Int}, G::SubsetGrayCode) = length(v) == G.k @inline function Base.isdone(G::SubsetGrayCode, state) (_, rank, _) = state + println("called") return rank == G.len end @@ -78,20 +79,20 @@ end Note that inds is part of the iterator's state inds only to prevent reallocation in each iteration. =# - subset_vec = zeros(UInt, G.k) - CodingTheory._subset_unrank(G.init_rank, UInt(G.n), subset_vec) + subset_vec = zeros(Int, G.k) + CodingTheory._subset_unrank!(G.init_rank, UInt(G.n), subset_vec) inds = fill(-1, 3) - (inds, (subset_vec, G.init_rank, inds)) + (inds, (subset_vec, 1, inds)) end @inline function Base.iterate(G::SubsetGrayCode, state) # Based on Algorithm 2.13 in kreher1999combinatorial - v, rank, inds = state + v, index, inds = state - if rank == G.len + if index == G.len return nothing end - rank += 1 + index += 1 @inbounds begin inds[1] = -1 @@ -144,7 +145,7 @@ end end end end - return (inds, (v, rank, inds)) + return (inds, (v, index, inds)) end function _update_indices!(indices::Vector{Int}, x::Int, y::Int) @@ -189,7 +190,7 @@ function _subset_rank(v::Vector{UInt}, k::UInt) return r end -function _subset_unrank!(r::BigInt, n::UInt, T::Vector{UInt}) +function _subset_unrank!(r::BigInt, n::UInt, T::Vector{Int}) # Based on Algorithm 2.12 in kreher1999combinatorial k = length(T) subset_size_str::String = "subset size k = $k must be smaller than the set size n = $n" diff --git a/test/iterators_test.jl b/test/iterators_test.jl index ec3566b..939cb39 100644 --- a/test/iterators_test.jl +++ b/test/iterators_test.jl @@ -1,5 +1,5 @@ @testitem "iterators.jl" begin - @testset "SubsetGrayCode iterates over all weight k subsets of {1, .., n} (nthreads == 1)" begin + @testset "SubsetGrayCode iterates over all weight k subsets of {1,..,n} (single iterator)" begin len = 15 weight = 7 sgc = CodingTheory.SubsetGrayCode(len, weight) @@ -20,18 +20,18 @@ @test all_subsets_gray == all_subsets_hecke end - @testset "SubsetGrayCode iterates over all weight k subsets of {1, .., n} (nthreads > 1)" begin + @testset "SubsetGrayCode iterates over all weight k subsets of {1,..,n} (split iterator)" begin len = 15 weight = 7 - num_threads = 3 - itrs = CodingTheory.SubsetGrayCodesFromNumThreads(len, weight, num_threads) + num_threads = 3 # split the iterator into 3 parts + itrs = CodingTheory._subset_gray_codes_from_num_threads(len, weight, num_threads) @test length(itrs) == num_threads initial_subset_vecs = fill(fill(0, weight), length(itrs)) bin = extended_binomial(len, weight) for j in eachindex(itrs) itr = itrs[j] @test length(itr) == fld(bin, num_threads) - subset_vec = zeros(UInt, itr.k) + subset_vec = zeros(Int, itr.k) CodingTheory._subset_unrank!(itr.init_rank, UInt(itr.n), subset_vec) initial_subset_vecs[j] = subset_vec end @@ -80,7 +80,7 @@ n1 = UInt(5) rank1 = CodingTheory._subset_rank(subset1, k1) @test rank1 == 0 - result1 = zeros(UInt, k1) + result1 = zeros(Int, k1) CodingTheory._subset_unrank!(rank1, n1, result1) @test result1 == subset1 @@ -89,13 +89,13 @@ n2 = UInt(5) rank2 = CodingTheory._subset_rank(subset2, k2) @test rank2 == 7 - result2 = zeros(UInt, k2) + result2 = zeros(Int, k2) CodingTheory._subset_unrank!(rank2, n2, result2) @test result2 == subset2 k3 = UInt(3) n3 = UInt(5) - result3 = zeros(UInt, k3) + result3 = zeros(Int, k3) results = Set() bin = binomial(n3, k3) for i::BigInt in collect(0: bin - 1) From d1a9c89d9eaaf68fb9912f3d6946605c03c16f62 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Thu, 23 Jan 2025 12:05:50 -0500 Subject: [PATCH 15/23] rename function --- src/Quantum/weight_dist.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Quantum/weight_dist.jl b/src/Quantum/weight_dist.jl index 08f58da..505b2d8 100644 --- a/src/Quantum/weight_dist.jl +++ b/src/Quantum/weight_dist.jl @@ -437,7 +437,7 @@ end Return the minimum distance of the stabilizer code if known, otherwise computes it. """ -function minimum_distance(S::AbstractStabilizerCode; alg::Symbol = :auto, verbose::Bool = false) +function minimum_distance_Gray(S::AbstractStabilizerCode; alg::Symbol = :auto, verbose::Bool = false) !ismissing(S.d) && return S.d # these should be different? weight? auto? BZ? @@ -577,7 +577,7 @@ end function is_pure(S::AbstractStabilizerCode) ismissing(S.pure) || return S.pure - minimum_distance(S) # this needs to force the weight enumerator approach + minimum_distance_Gray(S) # this needs to force the weight enumerator approach return S.pure end From ee5f03949fe72d4e0ed9b8b3e87648ddeba4540a Mon Sep 17 00:00:00 2001 From: Eric Sabo Date: Thu, 23 Jan 2025 11:04:56 -0500 Subject: [PATCH 16/23] Makie bug workaround --- Manifest.toml | 260 +++++++++-------- Project.toml | 32 ++- ext/JLD2Ext/JLD2Ext.jl | 4 - ext/JuMPExt/JuMPExt.jl | 5 - ext/MakieExt/Classical/Tanner.jl | 20 +- ext/MakieExt/LDPC/codes.jl | 132 ++++----- ext/MakieExt/MakieExt.jl | 27 +- ext/MakieExt/Quantum/weight_dist.jl | 20 +- src/Classical/TwistedReedSolomon.jl | 10 +- src/Classical/linear_code.jl | 44 +-- src/Classical/weight_dist.jl | 89 +++--- src/CodingTheory.jl | 2 + src/Convolutional/convolutional_code.jl | 366 ++++++++++++++++++++++++ src/Convolutional/types.jl | 43 +++ src/LDPC/types.jl | 2 +- src/Quantum/decoders/OTF.jl | 32 ++- src/iterators.jl | 9 +- test/LDPC/MP_decoders_test.jl | 10 +- test/Quantum/decoders/OFT_test.jl | 18 +- test/Quantum/product_codes_test.jl | 1 + 20 files changed, 779 insertions(+), 347 deletions(-) create mode 100644 src/Convolutional/convolutional_code.jl create mode 100644 src/Convolutional/types.jl diff --git a/Manifest.toml b/Manifest.toml index 157c29c..0dcfcda 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,8 +1,8 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.10.5" +julia_version = "1.11.3" manifest_format = "2.0" -project_hash = "163f71688893b1a2d676f6ba6e175eb2cf6dfafe" +project_hash = "bbb0e5646f1748a2c6f584642f98e12c63222479" [[deps.ASL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -11,10 +11,10 @@ uuid = "ae81ac8f-d209-56e5-92de-9978fef736f9" version = "0.1.3+0" [[deps.AbstractAlgebra]] -deps = ["InteractiveUtils", "LinearAlgebra", "MacroTools", "Preferences", "Random", "RandomExtensions", "SparseArrays", "Test"] -git-tree-sha1 = "f7cd148a76ce9668c3fe36f29fcc9201adf14265" +deps = ["LinearAlgebra", "MacroTools", "Preferences", "Random", "RandomExtensions", "SparseArrays", "Test"] +git-tree-sha1 = "4d02060414d46e141e79efbd1484901c0fb4ce24" uuid = "c3fe647b-3220-5bb0-a1ea-a7954cac585d" -version = "0.43.5" +version = "0.43.12" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] @@ -32,9 +32,9 @@ version = "1.5.0" [[deps.AlgebraicSolving]] deps = ["LinearAlgebra", "Logging", "Markdown", "Nemo", "Printf", "Random", "StaticArrays", "Test", "msolve_jll"] -git-tree-sha1 = "c68379967797de6ef0a9314f6c66d0803c455d48" +git-tree-sha1 = "406b5001ff9b0d62440d52b145b303e269e755a6" uuid = "66b61cbe-0446-4d5d-9090-1ff510639f9d" -version = "0.7.2" +version = "0.8.1" [[deps.AliasTables]] deps = ["PtrArrays", "Random"] @@ -44,7 +44,7 @@ version = "1.1.3" [[deps.ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" -version = "1.1.1" +version = "1.1.2" [[deps.ArnoldiMethod]] deps = ["LinearAlgebra", "Random", "StaticArrays"] @@ -54,6 +54,7 @@ version = "0.4.0" [[deps.Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +version = "1.11.0" [[deps.AutoHashEquals]] git-tree-sha1 = "4ec6b48702dacc5994a835c1189831755e4e76ef" @@ -62,6 +63,7 @@ version = "2.2.0" [[deps.Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +version = "1.11.0" [[deps.BinaryWrappers]] deps = ["JLLWrappers", "Scratch"] @@ -71,9 +73,9 @@ version = "0.1.3" [[deps.Bzip2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "9e2a6b69137e6969bab0152632dcb3bc108c8bdd" +git-tree-sha1 = "8873e196c2eb87962a2048b3b8e08946535864a1" uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" -version = "1.0.8+1" +version = "1.0.8+4" [[deps.Combinatorics]] git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860" @@ -97,9 +99,9 @@ version = "1.1.1+0" [[deps.CxxWrap]] deps = ["Libdl", "MacroTools", "libcxxwrap_julia_jll"] -git-tree-sha1 = "3345cb637ca1efb2ebf7f5145558522b92660d1f" +git-tree-sha1 = "9208ab6d13777dace5981b451171891faf59e6a2" uuid = "1f15a43c-97ca-5a2a-ae31-89f07a497df4" -version = "0.14.2" +version = "0.16.0" [[deps.DataAPI]] git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe" @@ -115,6 +117,7 @@ version = "0.18.20" [[deps.Dates]] deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +version = "1.11.0" [[deps.DecFP]] deps = ["DecFP_jll", "Printf", "Random", "SpecialFunctions"] @@ -131,12 +134,13 @@ version = "2.0.3+1" [[deps.Distributed]] deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" +version = "1.11.0" [[deps.Distributions]] deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"] -git-tree-sha1 = "d7477ecdafb813ddee2ae727afa94e9dcb5f3fb0" +git-tree-sha1 = "03aa5d44647eaec98e1920635cdfed5d5560a8b9" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" -version = "0.25.112" +version = "0.25.117" [deps.Distributions.extensions] DistributionsChainRulesCoreExt = "ChainRulesCore" @@ -169,16 +173,17 @@ version = "1.8.0" deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "4d81ed14783ec49ce9f2e168208a12ce1815aa25" uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" -version = "3.3.10+1" +version = "3.3.10+3" [[deps.FLINT_jll]] deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "MPFR_jll", "OpenBLAS32_jll"] -git-tree-sha1 = "4bde447e7bc6de73de36870fc942df02f4dd44b5" +git-tree-sha1 = "9327ebef3e04034b2aa1b9c59869f36717843ce5" uuid = "e134572f-a0d5-539d-bddf-3cad8db41a82" -version = "300.100.300+0" +version = "300.100.301+0" [[deps.FileWatching]] uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" +version = "1.11.0" [[deps.FillArrays]] deps = ["LinearAlgebra"] @@ -194,9 +199,9 @@ weakdeps = ["PDMats", "SparseArrays", "Statistics"] [[deps.GAP]] deps = ["AbstractAlgebra", "Artifacts", "Compat", "Downloads", "GAP_jll", "GAP_lib_jll", "GAP_pkg_juliainterface_jll", "InteractiveUtils", "Libdl", "MacroTools", "Markdown", "Ncurses_jll", "Pidfile", "Pkg", "REPL", "Random", "Scratch"] -git-tree-sha1 = "e251aa4e3693eed790e2ee77ec0126e128894a56" +git-tree-sha1 = "7edbab55aaed69062ea32c7030f2d477711cd98b" uuid = "c863536a-3901-11e9-33e7-d5cd0df7b904" -version = "0.12.0" +version = "0.12.3" [[deps.GAP_jll]] deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Readline_jll", "Zlib_jll"] @@ -218,14 +223,14 @@ version = "0.1200.0+0" [[deps.GLPK_jll]] deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "fe68622f32828aa92275895fdb324a85894a5b1b" +git-tree-sha1 = "6aa6294ba949ccfc380463bf50ff988b46de5bc7" uuid = "e8aa6df9-e6ca-548a-97ff-1f85fc5b8b98" -version = "5.0.1+0" +version = "5.0.1+1" [[deps.GMP_jll]] deps = ["Artifacts", "Libdl"] uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d" -version = "6.2.1+6" +version = "6.3.0+0" [[deps.Graphs]] deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"] @@ -235,9 +240,9 @@ version = "1.12.0" [[deps.Hecke]] deps = ["AbstractAlgebra", "Dates", "Distributed", "InteractiveUtils", "LazyArtifacts", "Libdl", "LinearAlgebra", "Markdown", "Nemo", "Pkg", "Printf", "Random", "RandomExtensions", "Serialization", "SparseArrays"] -git-tree-sha1 = "00a2e1d65c458948530c3c68ab519e6b7b80f053" +git-tree-sha1 = "ee9546ec43de9fa083094a40c999927b9c005b29" uuid = "3e1990a7-5d81-5526-99ce-9ba3ff248f21" -version = "0.34.4" +version = "0.34.9" weakdeps = ["GAP", "Polymake"] [deps.Hecke.extensions] @@ -246,15 +251,15 @@ weakdeps = ["GAP", "Polymake"] [[deps.Hwloc_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "dd3b49277ec2bb2c6b94eb1604d4d0616016f7a6" +git-tree-sha1 = "50aedf345a709ab75872f80a2779568dc0bb461b" uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8" -version = "2.11.2+0" +version = "2.11.2+3" [[deps.HypergeometricFunctions]] deps = ["LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"] -git-tree-sha1 = "7c4195be1649ae622304031ed46a2f4df989f1eb" +git-tree-sha1 = "b1c2585431c382e3fe5805874bda6aea90a95de9" uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" -version = "0.3.24" +version = "0.3.25" [[deps.Inflate]] git-tree-sha1 = "d1b1b796e47d94588b3757fe84fbf65a5ec4a80d" @@ -270,12 +275,13 @@ version = "2024.2.1+0" [[deps.InteractiveUtils]] deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +version = "1.11.0" [[deps.Ipopt_jll]] deps = ["ASL_jll", "Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "MUMPS_seq_jll", "SPRAL_jll", "libblastrampoline_jll"] -git-tree-sha1 = "a0950d209a055b3adb6d29ade5cbdf005a6bd290" +git-tree-sha1 = "4f55ad688c698a4f77d892a1cb673f7e8a30f178" uuid = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7" -version = "300.1400.1600+0" +version = "300.1400.1700+0" [[deps.IrrationalConstants]] git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" @@ -284,9 +290,9 @@ version = "0.2.2" [[deps.JLLWrappers]] deps = ["Artifacts", "Preferences"] -git-tree-sha1 = "f389674c99bfcde17dc57454011aa44d5a260a40" +git-tree-sha1 = "a007feb38b422fbdab534406aeca1b86823cb4d6" uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" -version = "1.6.0" +version = "1.7.0" [[deps.JSON]] deps = ["Dates", "Mmap", "Parsers", "Unicode"] @@ -296,9 +302,9 @@ version = "0.21.4" [[deps.JSON3]] deps = ["Dates", "Mmap", "Parsers", "PrecompileTools", "StructTypes", "UUIDs"] -git-tree-sha1 = "eb3edce0ed4fa32f75a0a11217433c31d56bd48b" +git-tree-sha1 = "1d322381ef7b087548321d3f878cb4c9bd8f8f9b" uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" -version = "1.14.0" +version = "1.14.1" [deps.JSON3.extensions] JSON3ArrowExt = ["ArrowTypes"] @@ -314,13 +320,14 @@ version = "18.1.7+0" [[deps.LZO_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "854a9c268c43b77b0a27f22d7fab8d33cdb3a731" +git-tree-sha1 = "1c602b1127f4751facb671441ca72715cc95938a" uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac" -version = "2.10.2+1" +version = "2.10.3+0" [[deps.LazyArtifacts]] deps = ["Artifacts", "Pkg"] uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" +version = "1.11.0" [[deps.LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] @@ -330,16 +337,17 @@ version = "0.6.4" [[deps.LibCURL_jll]] deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "8.4.0+0" +version = "8.6.0+0" [[deps.LibGit2]] deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" +version = "1.11.0" [[deps.LibGit2_jll]] deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" -version = "1.6.4+0" +version = "1.7.2+0" [[deps.LibSSH2_jll]] deps = ["Artifacts", "Libdl", "MbedTLS_jll"] @@ -348,16 +356,18 @@ version = "1.11.0+1" [[deps.Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +version = "1.11.0" [[deps.LinearAlgebra]] deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +version = "1.11.0" [[deps.LogExpFunctions]] deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "a2d09619db4e765091ee5c6ffe8872849de0feea" +git-tree-sha1 = "13ca9e2586b89836fd20cccf56e57e2b9ae7f38f" uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.28" +version = "0.3.29" [deps.LogExpFunctions.extensions] LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" @@ -371,18 +381,19 @@ version = "0.3.28" [[deps.Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" +version = "1.11.0" [[deps.Lz4_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "abf88ff67f4fd89839efcae2f4c39cbc4ecd0846" +git-tree-sha1 = "191686b1ac1ea9c89fc52e996ad15d1d241d1e33" uuid = "5ced341a-0733-55b8-9ab6-a4889d929147" -version = "1.10.0+1" +version = "1.10.1+0" [[deps.METIS_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "1fd0a97409e418b78c53fac671cf4622efdf0f21" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "2eefa8baa858871ae7770c98c3c2a7e46daba5b4" uuid = "d00139f3-1899-568f-a2f0-47f597d42d70" -version = "5.1.2+0" +version = "5.1.3+0" [[deps.MKL_jll]] deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "oneTBB_jll"] @@ -393,28 +404,28 @@ version = "2024.2.0+0" [[deps.MPFR_jll]] deps = ["Artifacts", "GMP_jll", "Libdl"] uuid = "3a97d323-0669-5f0c-9066-3539efd106a3" -version = "4.2.0+1" +version = "4.2.1+0" [[deps.MUMPS_seq_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "METIS_jll", "libblastrampoline_jll"] -git-tree-sha1 = "85047ac569761e3387717480a38a61d2a67df45c" +git-tree-sha1 = "0eab12f94948ca67908aec14b9f2ebefd17463fe" uuid = "d7ed1dd3-d0ae-5e8e-bfb4-87a502085b8d" -version = "500.700.300+0" +version = "500.700.301+0" [[deps.MacroTools]] -deps = ["Markdown", "Random"] -git-tree-sha1 = "2fa9ee3e63fd3a4f7a9a4f4744a52f4856de82df" +git-tree-sha1 = "72aebe0b5051e5143a079a4685a46da330a40472" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.13" +version = "0.5.15" [[deps.Markdown]] deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" +version = "1.11.0" [[deps.MbedTLS_jll]] deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.2+1" +version = "2.28.6+0" [[deps.Missings]] deps = ["DataAPI"] @@ -424,22 +435,23 @@ version = "1.2.0" [[deps.Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" +version = "1.11.0" [[deps.MongoC_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "OpenSSL_jll", "Zlib_jll", "Zstd_jll", "snappy_jll"] -git-tree-sha1 = "5a0f9f14a8186eae48608ff7922ac0c00ff52cdc" +git-tree-sha1 = "0ab1e13479581def03ea0b19733b8b830177889a" uuid = "90100e71-7732-535a-9be7-2e9affd1cfc1" -version = "1.25.1+0" +version = "1.28.1+0" [[deps.Mongoc]] deps = ["Dates", "DecFP", "MongoC_jll", "Serialization"] -git-tree-sha1 = "1f0d7579d1bacc1446751990d7e0c4c0bf0f3277" +git-tree-sha1 = "9b0ee6f5843a8df838d76794a5d50e77b8a3b08f" uuid = "4fe8b98c-fc19-5c23-8ec2-168ff83495f2" -version = "0.9.1" +version = "0.9.2" [[deps.MozillaCACerts_jll]] uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2023.1.10" +version = "2023.12.12" [[deps.Ncurses_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -449,9 +461,9 @@ version = "6.5.0+1" [[deps.Nemo]] deps = ["AbstractAlgebra", "FLINT_jll", "Libdl", "LinearAlgebra", "Pkg", "Random", "RandomExtensions", "SHA"] -git-tree-sha1 = "a644b4943424d9f1ffe7a7a1cdda617053bac5c0" +git-tree-sha1 = "c3960b4fe367d62818d03ba273ab8b00b96a34fa" uuid = "2edaba10-b0f1-5616-af89-8c11ac63239a" -version = "0.47.1" +version = "0.47.5" [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" @@ -465,14 +477,14 @@ version = "1.12.1+0" [[deps.OpenBLAS32_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] -git-tree-sha1 = "6065c4cff8fee6c6770b277af45d5082baacdba1" +git-tree-sha1 = "ece4587683695fe4c5f20e990da0ed7e83c351e7" uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2" -version = "0.3.24+0" +version = "0.3.29+0" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.23+4" +version = "0.3.27+1" [[deps.OpenLibm_jll]] deps = ["Artifacts", "Libdl"] @@ -481,34 +493,32 @@ version = "0.8.1+2" [[deps.OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "ad31332567b189f508a3ea8957a2640b1147ab00" +git-tree-sha1 = "7493f61f55a6cce7325f197443aa80d32554ba10" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "1.1.23+1" +version = "3.0.15+3" [[deps.OpenSpecFun_jll]] -deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1346c9208249809840c91b26703912dff463d335" uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" -version = "0.5.5+0" +version = "0.5.6+0" [[deps.OrderedCollections]] -git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5" +git-tree-sha1 = "12f1439c4f986bb868acda6ea33ebc78e19b95ad" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.6.3" +version = "1.7.0" [[deps.Oscar]] deps = ["AbstractAlgebra", "AlgebraicSolving", "Distributed", "GAP", "Hecke", "JSON", "JSON3", "LazyArtifacts", "Markdown", "Nemo", "Pkg", "Polymake", "Random", "RandomExtensions", "Serialization", "Singular", "TOPCOM_jll", "UUIDs", "cohomCalg_jll"] -git-tree-sha1 = "9aec1c8a631d474902d1274d6b4d683e4ae22357" -repo-rev = "master" -repo-url = "https://github.com/oscar-system/Oscar.jl.git" +git-tree-sha1 = "803604e2666bb0f007dced41968586e5134b9672" uuid = "f1435218-dba5-11e9-1e4d-f1a5fab5fc13" -version = "1.2.0-DEV" +version = "1.2.2" [[deps.PDMats]] deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "949347156c25054de2db3b166c52ac4728cbad65" +git-tree-sha1 = "966b85253e959ea89c53a9abebbf2e964fbf593b" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.31" +version = "0.11.32" [[deps.PPL_jll]] deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"] @@ -535,15 +545,19 @@ uuid = "fa939f87-e72e-5be4-a000-7fc836dbe307" version = "1.3.0" [[deps.Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] +deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.10.0" +version = "1.11.0" +weakdeps = ["REPL"] + + [deps.Pkg.extensions] + REPLExt = "REPL" [[deps.Polymake]] deps = ["AbstractAlgebra", "BinaryWrappers", "CxxWrap", "Downloads", "JSON", "Libdl", "Mongoc", "NetworkOptions", "Ninja_jll", "Perl_jll", "Pidfile", "Pkg", "REPL", "Scratch", "SparseArrays", "TOPCOM_jll", "lib4ti2_jll", "libpolymake_julia_jll", "polymake_jll", "polymake_oscarnumber_jll"] -git-tree-sha1 = "e065263f7e818a419bdc89979d7196222bb70efb" +git-tree-sha1 = "0a606d37843442f11b972046610eb94216b6089c" uuid = "d720cf60-89b5-51f5-aff5-213f193123e7" -version = "0.11.22" +version = "0.11.24" [[deps.PrecompileTools]] deps = ["Preferences"] @@ -560,6 +574,7 @@ version = "1.4.3" [[deps.Printf]] deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +version = "1.11.0" [[deps.ProgressMeter]] deps = ["Distributed", "Printf"] @@ -568,9 +583,9 @@ uuid = "92933f4c-e287-5a05-a399-4b506db050ca" version = "1.10.2" [[deps.PtrArrays]] -git-tree-sha1 = "77a42d78b6a92df47ab37e177b2deac405e1c88f" +git-tree-sha1 = "1d36ef11a9aaf1e8b74dacc6a731dd1de8fd493d" uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d" -version = "1.2.1" +version = "1.3.0" [[deps.QuadGK]] deps = ["DataStructures", "LinearAlgebra"] @@ -585,12 +600,14 @@ version = "2.11.1" Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" [[deps.REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] +deps = ["InteractiveUtils", "Markdown", "Sockets", "StyledStrings", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" +version = "1.11.0" [[deps.Random]] deps = ["SHA"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +version = "1.11.0" [[deps.RandomExtensions]] deps = ["Random", "SparseArrays"] @@ -645,10 +662,12 @@ version = "1.2.1" [[deps.Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" +version = "1.11.0" [[deps.SharedArrays]] deps = ["Distributed", "Mmap", "Random", "Serialization"] uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" +version = "1.11.0" [[deps.SimpleTraits]] deps = ["InteractiveUtils", "MacroTools"] @@ -658,18 +677,19 @@ version = "0.9.4" [[deps.Singular]] deps = ["AbstractAlgebra", "BinaryWrappers", "CxxWrap", "Libdl", "LinearAlgebra", "Nemo", "Pidfile", "Pkg", "Random", "RandomExtensions", "Singular_jll", "Statistics", "lib4ti2_jll", "libsingular_julia_jll"] -git-tree-sha1 = "bc00189bb74b6a3f774253064f945d09fb1668ab" +git-tree-sha1 = "8ec96ba2fc0137960c79d799e8c8f2b0296d5396" uuid = "bcd08a7b-43d2-5ff7-b6d4-c458787f915c" -version = "0.23.8" +version = "0.23.10" [[deps.Singular_jll]] deps = ["Artifacts", "FLINT_jll", "GMP_jll", "JLLWrappers", "Libdl", "MPFR_jll", "cddlib_jll"] -git-tree-sha1 = "6c85174749476dcd3f059a33884072d748033263" +git-tree-sha1 = "7038aa73c007d59620ea2a2db55698326a304266" uuid = "43d676ae-4934-50ba-8046-7a96366d613b" -version = "404.0.606+0" +version = "404.0.711+0" [[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" +version = "1.11.0" [[deps.SortingAlgorithms]] deps = ["DataStructures"] @@ -680,13 +700,13 @@ version = "1.2.1" [[deps.SparseArrays]] deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -version = "1.10.0" +version = "1.11.0" [[deps.SpecialFunctions]] deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "2f5d4697f21388cbe1ff299430dd169ef97d7e14" +git-tree-sha1 = "64cca0c26b4f31ba18f13f6c12af7c85f478cfde" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.4.0" +version = "2.5.0" [deps.SpecialFunctions.extensions] SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" @@ -696,9 +716,9 @@ version = "2.4.0" [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "eeafab08ae20c62c44c8399ccb9354a04b80db50" +git-tree-sha1 = "47091a0340a675c738b1304b58161f3b0839d454" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.9.7" +version = "1.9.10" [deps.StaticArrays.extensions] StaticArraysChainRulesCoreExt = "ChainRulesCore" @@ -714,9 +734,14 @@ uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" version = "1.4.3" [[deps.Statistics]] -deps = ["LinearAlgebra", "SparseArrays"] +deps = ["LinearAlgebra"] +git-tree-sha1 = "ae3bb1eb3bba077cd276bc5cfc337cc65c3075c0" uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -version = "1.10.0" +version = "1.11.1" +weakdeps = ["SparseArrays"] + + [deps.Statistics.extensions] + SparseArraysExt = ["SparseArrays"] [[deps.StatsAPI]] deps = ["LinearAlgebra"] @@ -725,10 +750,10 @@ uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" version = "1.7.0" [[deps.StatsBase]] -deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] -git-tree-sha1 = "5cf7606d6cef84b543b483848d4ae08ad9832b21" +deps = ["AliasTables", "DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] +git-tree-sha1 = "29321314c920c26684834965ec2ce0dacc9cf8e5" uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.34.3" +version = "0.34.4" [[deps.StatsFuns]] deps = ["HypergeometricFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] @@ -750,6 +775,10 @@ git-tree-sha1 = "159331b30e94d7b11379037feeb9b690950cace8" uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4" version = "1.11.0" +[[deps.StyledStrings]] +uuid = "f489334b-da3d-4c2e-b8f0-e476e12c162b" +version = "1.11.0" + [[deps.SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" @@ -757,7 +786,7 @@ uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" [[deps.SuiteSparse_jll]] deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" -version = "7.2.1+1" +version = "7.7.0+0" [[deps.TOML]] deps = ["Dates"] @@ -778,6 +807,7 @@ version = "1.10.0" [[deps.Test]] deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +version = "1.11.0" [[deps.TestItemRunner]] deps = ["Pkg", "TOML", "Test", "TestItems", "UUIDs"] @@ -793,9 +823,11 @@ version = "1.0.0" [[deps.UUIDs]] deps = ["Random", "SHA"] uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" +version = "1.11.0" [[deps.Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" +version = "1.11.0" [[deps.Zlib_jll]] deps = ["Libdl"] @@ -804,9 +836,9 @@ version = "1.2.13+1" [[deps.Zstd_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "555d1076590a6cc2fdee2ef1469451f872d8b41b" +git-tree-sha1 = "622cf78670d067c738667aaa96c553430b65e269" uuid = "3161d3a3-bdf6-5164-811a-617609db77b4" -version = "1.5.6+1" +version = "1.5.7+0" [[deps.bliss_jll]] deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"] @@ -845,21 +877,21 @@ version = "5.11.0+0" [[deps.libcxxwrap_julia_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "02d0a0a623248c709727088aaf722ab14f1463a5" +git-tree-sha1 = "b386562c39065094e2d1e1939ea289dbd448fe3e" uuid = "3eaa8342-bff7-56a5-9981-c04077f7cee7" -version = "0.11.2+1" +version = "0.13.3+0" [[deps.libpolymake_julia_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "FLINT_jll", "JLLWrappers", "Libdl", "TOPCOM_jll", "lib4ti2_jll", "libcxxwrap_julia_jll", "polymake_jll"] -git-tree-sha1 = "13e1cd24006a6c9f528125cc736fdf2981fba5c2" +git-tree-sha1 = "4885c6d17996ea3fa5492cfe8c68f0f52be82820" uuid = "4d8266f6-2b3b-57e3-ad7a-d431eaaac945" -version = "0.12.1+0" +version = "0.13.0+0" [[deps.libsingular_julia_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Singular_jll", "libcxxwrap_julia_jll"] -git-tree-sha1 = "436efe88c41337fd128062555fe8f8ea6c0e7309" +git-tree-sha1 = "8db9abc14f74d1e17d66466d5ee8a58581e3f3b8" uuid = "ae4fbd8f-ecdb-54f8-bbce-35570499b30e" -version = "0.45.6+0" +version = "0.46.1+0" [[deps.lrslib_jll]] deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"] @@ -869,9 +901,9 @@ version = "0.3.3+0" [[deps.msolve_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "FLINT_jll", "GMP_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl", "MPFR_jll"] -git-tree-sha1 = "7324fd011ebc4b9e43d9aa365c79797860459d0a" +git-tree-sha1 = "382191ffa0d6d873793b04659c3a99e6a332caca" uuid = "6d01cc9a-e8f6-580e-8c54-544227e08205" -version = "0.700.200+0" +version = "0.700.301+0" [[deps.nauty_jll]] deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"] @@ -882,7 +914,7 @@ version = "2.6.13+1" [[deps.nghttp2_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.52.0+1" +version = "1.59.0+0" [[deps.normaliz_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "FLINT_jll", "GMP_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl", "MPFR_jll", "nauty_jll"] @@ -903,18 +935,18 @@ version = "17.4.0+2" [[deps.polymake_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "FLINT_jll", "GMP_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl", "MPFR_jll", "MongoC_jll", "PPL_jll", "Perl_jll", "SCIP_jll", "bliss_jll", "boost_jll", "cddlib_jll", "lrslib_jll", "normaliz_jll"] -git-tree-sha1 = "bd0b5a3dd1f0d1309d6314b207d5e2bc51d42c29" +git-tree-sha1 = "bfd75c077ece1a66325ec61a6cb6b0b23420f4b2" uuid = "7c209550-9012-526c-9264-55ba7a78ba2c" -version = "400.1200.1+0" +version = "400.1300.2+0" [[deps.polymake_oscarnumber_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl", "libcxxwrap_julia_jll", "libpolymake_julia_jll", "polymake_jll"] -git-tree-sha1 = "e417d0fe884dc35f60c913a94c6ef3340f9a9d4f" +git-tree-sha1 = "0f1f590da4a2323defeecbf50a1d21448e7108fc" uuid = "10f31823-b687-53e6-9f29-edb9d4da9f9f" -version = "0.3.1+0" +version = "0.3.3+1" [[deps.snappy_jll]] deps = ["Artifacts", "JLLWrappers", "LZO_jll", "Libdl", "Lz4_jll", "Zlib_jll"] -git-tree-sha1 = "d34e37153e4c03205bde2f0f8ada6cced8c83978" +git-tree-sha1 = "a6c92533edec9755b3c1083b028da19a41f92e25" uuid = "fe1e1685-f7be-5f59-ac9f-4ca204017dfd" -version = "1.2.1+0" +version = "1.2.1+3" diff --git a/Project.toml b/Project.toml index 09723eb..359b24e 100644 --- a/Project.toml +++ b/Project.toml @@ -8,12 +8,12 @@ AutoHashEquals = "15f4f7f2-30c1-5605-9d31-71845cf9641f" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" -JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Oscar = "f1435218-dba5-11e9-1e4d-f1a5fab5fc13" -Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -37,24 +37,26 @@ WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008" [extensions] JLD2Ext = "JLD2" JuMPExt = ["JuMP", "GLPK"] -MakieExt = ["Makie", "NetworkLayout", "CairoMakie", "GraphMakie", "GLMakie", "WGLMakie", "GraphPlot"] +MakieExt = ["Makie"] [compat] -AutoHashEquals = "2.1.0" -CairoMakie = "0.10.11" -Colors = "0.12.10" +AutoHashEquals = "2.2.0" +CairoMakie = "0.13.1" +Colors = "0.13.0" Combinatorics = "1.0.2" -DataStructures = "0.18.15" -FFTW = "1.7.1" -GLMakie = "0.8.11" -GraphMakie = "0.5.6" -GraphPlot = "0.5.2" -Graphs = "1.9.0" +DataStructures = "0.18.20" +DocStringExtensions = "0.9.3" +FFTW = "1.8.0" +GLMakie = "0.11.2" +GraphMakie = "0.5.13" +GraphPlot = "0.6.0" +Graphs = "1.12.0" Makie = "0.19.11" -NetworkLayout = "0.4.6" -StatsBase = "0.34.2" +NetworkLayout = "0.4.9" +Oscar = "1.2.2" +StatsBase = "0.34.4" WGLMakie = "0.8.15" -julia = "≥ 1.9" +julia = "≥ 1.10" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/ext/JLD2Ext/JLD2Ext.jl b/ext/JLD2Ext/JLD2Ext.jl index 8d208b6..8e90cc3 100644 --- a/ext/JLD2Ext/JLD2Ext.jl +++ b/ext/JLD2Ext/JLD2Ext.jl @@ -1,13 +1,9 @@ module JLD2Ext import CodingTheory - import CodingTheory: TriangularColorCode488, TriangularColorCode666, StabilizerCode, set_logicals!, set_minimum_distance! #PlanarSurfaceCode3D, PlanarSurfaceCode3D_X, ToricCode3D, - import JLD2 - import JLD2: @load - import Oscar: GF, matrix include("Quantum/misc_known_codes.jl") diff --git a/ext/JuMPExt/JuMPExt.jl b/ext/JuMPExt/JuMPExt.jl index fce6f83..76295f2 100644 --- a/ext/JuMPExt/JuMPExt.jl +++ b/ext/JuMPExt/JuMPExt.jl @@ -1,15 +1,10 @@ module JuMPExt import CodingTheory - import CodingTheory: optimal_lambda, optimal_rho, optimal_lambda_and_rho, LP_decoder_LDPC, AbstractLinearCode, BinarySymmetricChannel - import JuMP - import JuMP: @variable, @constraint, @objective - import GLPK - import Oscar include("LDPC/decoders.jl") diff --git a/ext/MakieExt/Classical/Tanner.jl b/ext/MakieExt/Classical/Tanner.jl index 89a9cd4..c100446 100644 --- a/ext/MakieExt/Classical/Tanner.jl +++ b/ext/MakieExt/Classical/Tanner.jl @@ -20,34 +20,34 @@ function CodingTheory.Tanner_graph_plot(H::Union{CodingTheory.CTMatrixTypes, Mat # put in top right corner in order to get parents, children working A[1:nc, nc + 1:end] = transpose(M) - fig = Figure(); - ax = Axis(fig[1, 1], yreversed = true, xautolimitmargin = (0.15, 0.20), + fig = CairoMakie.Figure(); + ax = CairoMakie.Axis(fig[1, 1], yreversed = true, xautolimitmargin = (0.15, 0.20), yautolimitmargin = (0.15, 0.20)) - hidespines!(ax) - hidedecorations!(ax) + CairoMakie.hidespines!(ax) + CairoMakie.hidedecorations!(ax) left_x, left_y = zeros(nc), 1.:nc right_x, right_y = ones(nr) * nr, range(1, nc, nr) x = vcat(left_x, right_x) y = vcat(left_y, right_y) - points = CairoMakie.Point.(zip(x, y)) + points = CairoMakie.Point2f.(zip(x, y)) cols = (:aqua, :red, :orange, :green, :blue, :purple) - G = SimpleDiGraph(A) + G = Grphs.SimpleDiGraph(A) parents = [Grphs.inneighbors(G, i) for i in Grphs.vertices(G)] children = findall(x -> length(x) > 0, parents) for (i, v) in enumerate(children) for node in parents[v] - lines!(CairoMakie.Point(x[[node, v]]...), CairoMakie.Point(y[[node, v]]...), - color = cols[i % 6 + 1], linewidth = 5) + CairoMakie.lines!([CairoMakie.Point2f(x[[node, v]])...], [CairoMakie.Point2f(y[[node, + v]])...], color = cols[i % 6 + 1], linewidth = 5) end - text!(points[v], text = L"c_{%$i}", offset = (20, -15)) + CairoMakie.text!(points[v], text = L"c_{%$i}", offset = (20, -15)) end for (i, point) in enumerate(points[1:nc]) CairoMakie.scatter!(point, color = :black, marker = :circle, markersize = 25) - text!(point, text = L"v_{%$i}", offset = (-30, -10)) + CairoMakie.text!(point, text = L"v_{%$i}", offset = (-30, -10)) end for (i, point) in enumerate(points[nc + 1:end]) diff --git a/ext/MakieExt/LDPC/codes.jl b/ext/MakieExt/LDPC/codes.jl index 964676d..8c17f4e 100644 --- a/ext/MakieExt/LDPC/codes.jl +++ b/ext/MakieExt/LDPC/codes.jl @@ -43,81 +43,81 @@ function CodingTheory.degree_distributions_plot(C::AbstractLDPCCode) return f end -""" -$TYPEDSIGNATURES +# """ +# $TYPEDSIGNATURES -Return a bar graph and a dictionary of (length, count) pairs for unique short -cycles in the Tanner graph of `C`. An empty graph and dictionary are returned -when there are no cycles. +# Return a bar graph and a dictionary of (length, count) pairs for unique short +# cycles in the Tanner graph of `C`. An empty graph and dictionary are returned +# when there are no cycles. -# Note -- Short cycles are defined to be those with lengths between ``g`` and ``2g - 2``, - where ``g`` is the girth. -- Run `using Makie` to activate this extension. -""" -function CodingTheory.count_short_cycles_plot(C::AbstractLDPCCode) - if isempty(C.short_cycle_counts) || isempty(C.elementary_cycle_counts) - CodingTheory._count_cycles(C) - end +# # Note +# - Short cycles are defined to be those with lengths between ``g`` and ``2g - 2``, +# where ``g`` is the girth. +# - Run `using Makie` to activate this extension. +# """ +# function CodingTheory.count_short_cycles_plot(C::AbstractLDPCCode) +# if isempty(C.short_cycle_counts) || isempty(C.elementary_cycle_counts) +# CodingTheory._count_cycles(C) +# end - len = length(C.short_cycle_counts) - x_data = [0 for _ in 1:len] - y_data = [0 for _ in 1:len] - index = 1 - for (i, j) in C.short_cycle_counts - x_data[index] = i - y_data[index] = j - index += 1 - end +# len = length(C.short_cycle_counts) +# x_data = [0 for _ in 1:len] +# y_data = [0 for _ in 1:len] +# index = 1 +# for (i, j) in C.short_cycle_counts +# x_data[index] = i +# y_data[index] = j +# index += 1 +# end - fig = Figure(); - ax = Axis(fig[1, 1], xlabel = "Cycle Length", ylabel = "Occurrences", - title = "Short Cycle Counts") - barplot!(ax, x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data) - # fig = Plots.bar(x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data, - # legend = false, xlabel = "Cycle Length", ylabel = "Occurrences", - # title = "Short Cycle Counts") - display(fig) - return fig, C.short_cycle_counts -end +# fig = Figure(); +# ax = Axis(fig[1, 1], xlabel = "Cycle Length", ylabel = "Occurrences", +# title = "Short Cycle Counts") +# barplot!(ax, x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data) +# # fig = Plots.bar(x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data, +# # legend = false, xlabel = "Cycle Length", ylabel = "Occurrences", +# # title = "Short Cycle Counts") +# display(fig) +# return fig, C.short_cycle_counts +# end -""" -$TYPEDSIGNATURES +# """ +# $TYPEDSIGNATURES -Return a bar graph and a dictionary of (length, count) pairs for unique elementary -cycles in the Tanner graph of `C`. An empty graph and dictionary are returned -when there are no cycles. +# Return a bar graph and a dictionary of (length, count) pairs for unique elementary +# cycles in the Tanner graph of `C`. An empty graph and dictionary are returned +# when there are no cycles. -# Note -- Elementary cycles do not contain the same vertex twice and are unable to be - decomposed into a sequence of shorter cycles. -- Run `using Makie` to activate this extension. -""" -function CodingTheory.count_elementary_cycles_plot(C::AbstractLDPCCode) - if isempty(C.short_cycle_counts) || isempty(C.elementary_cycle_counts) - CodingTheory._count_cycles(C) - end +# # Note +# - Elementary cycles do not contain the same vertex twice and are unable to be +# decomposed into a sequence of shorter cycles. +# - Run `using Makie` to activate this extension. +# """ +# function CodingTheory.count_elementary_cycles_plot(C::AbstractLDPCCode) +# if isempty(C.short_cycle_counts) || isempty(C.elementary_cycle_counts) +# CodingTheory._count_cycles(C) +# end - len = length(C.elementary_cycle_counts) - x_data = [0 for _ in 1:len] - y_data = [0 for _ in 1:len] - index = 1 - for (i, j) in C.elementary_cycle_counts - x_data[index] = i - y_data[index] = j - index += 1 - end +# len = length(C.elementary_cycle_counts) +# x_data = [0 for _ in 1:len] +# y_data = [0 for _ in 1:len] +# index = 1 +# for (i, j) in C.elementary_cycle_counts +# x_data[index] = i +# y_data[index] = j +# index += 1 +# end - fig = Figure(); - ax = Axis(fig[1, 1], xlabel = "Cycle Length", ylabel = "Occurrences", - title = "Elementary Cycle Counts") - barplot!(ax, x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data) - # fig = Plots.bar(x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data, - # legend = false, xlabel = "Cycle Length", ylabel = "Occurrences", - # title = "Elementary Cycle Counts") - display(fig) - return fig, C.elementary_cycle_counts -end +# fig = Figure(); +# ax = Axis(fig[1, 1], xlabel = "Cycle Length", ylabel = "Occurrences", +# title = "Elementary Cycle Counts") +# barplot!(ax, x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data) +# # fig = Plots.bar(x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data, +# # legend = false, xlabel = "Cycle Length", ylabel = "Occurrences", +# # title = "Elementary Cycle Counts") +# display(fig) +# return fig, C.elementary_cycle_counts +# end """ $TYPEDSIGNATURES diff --git a/ext/MakieExt/MakieExt.jl b/ext/MakieExt/MakieExt.jl index 7c052d6..80d5d10 100644 --- a/ext/MakieExt/MakieExt.jl +++ b/ext/MakieExt/MakieExt.jl @@ -1,26 +1,13 @@ module MakieExt -import CodingTheory - -import CodingTheory: Tanner_graph_plot, weight_plot, EXIT_chart_plot, degree_distributions_plot, count_short_cycles_plot, count_elementary_cycles_plot, ACE_spectrum_plot, computation_graph, weight_plot_CSS_X, weight_plot_CSS_Z, weight_plot_CSS - -import Makie - -import NetworkLayout - -import CairoMakie - -import GraphMakie - -import GLMakie - -import WGLMakie - -import GraphPlot - -import Oscar - +import CodingTheory, Oscar +import CodingTheory: Tanner_graph_plot, weight_plot, EXIT_chart_plot, degree_distributions_plot, count_short_cycles_plot, count_elementary_cycles_plot, ACE_spectrum_plot, computation_graph, weight_plot_CSS_X, weight_plot_CSS_Z, weight_plot_CSS, AbstractLinearCode, AbstractLDPCCode, + LDPCEnsemble, AbstractClassicalNoiseChannel, AbstractStabilizerCode, AbstractStabilizerCodeCSS +# import Makie +using Makie, NetworkLayout, CairoMakie, GraphMakie, GLMakie, WGLMakie#, GraphPlot +import CairoMakie: @L_str #, Figure, Axis, hidespines!, hidedecorations!, lines!, text! import Graphs as Grphs +using DocStringExtensions # import CairoMakie: save include("Classical/Tanner.jl") diff --git a/ext/MakieExt/Quantum/weight_dist.jl b/ext/MakieExt/Quantum/weight_dist.jl index ed552d6..09ca6ef 100644 --- a/ext/MakieExt/Quantum/weight_dist.jl +++ b/ext/MakieExt/Quantum/weight_dist.jl @@ -8,19 +8,19 @@ # Weight Enumerators ############################# -""" -$TYPEDSIGNATURES +# """ +# $TYPEDSIGNATURES -Return a bar graph of the weight distribution related to `S`. +# Return a bar graph of the weight distribution related to `S`. -If `type` is `:stabilizer`, the weight distribution of the stabilizers are computed. -If `type` is `:normalizer`, the weight distrbution of the normalizer of the stabilizers -are computed. If `type` is `:quotient`, the weight distrbution of the normalizer mod the -stabilizers is computed. +# If `type` is `:stabilizer`, the weight distribution of the stabilizers are computed. +# If `type` is `:normalizer`, the weight distrbution of the normalizer of the stabilizers +# are computed. If `type` is `:quotient`, the weight distrbution of the normalizer mod the +# stabilizers is computed. -# Note -- Run `using Makie` to activate this extension. -""" +# # Note +# - Run `using Makie` to activate this extension. +# """ function CodingTheory.weight_plot(S::AbstractStabilizerCode; alg::Symbol = :auto, type::Symbol = :stabilizer) diff --git a/src/Classical/TwistedReedSolomon.jl b/src/Classical/TwistedReedSolomon.jl index cb7b80f..57fdf36 100644 --- a/src/Classical/TwistedReedSolomon.jl +++ b/src/Classical/TwistedReedSolomon.jl @@ -40,7 +40,7 @@ function TwistedReedSolomonCode(k::Int, α::Vector{T}, t::Vector{Int}, h::Vector G[i + 1, c] = g_i(α[c]) end end - display(G) + # display(G) t_dual = k .- h h_dual = (n - k) .- t @@ -55,11 +55,11 @@ function TwistedReedSolomonCode(k::Int, α::Vector{T}, t::Vector{Int}, h::Vector H[i + 1, c] = g_i(α[c]) end end - println(" ") - display(H) + # println(" ") + # display(H) - println(" ") - display(G * transpose(H)) + # println(" ") + # display(G * transpose(H)) C = LinearCode(G, H, false) return TwistedReedSolomonCode(C.F, C.n, C.k, C.d, C.l_bound, C.u_bound, C.G, C.H, C.G_stand, C.H_stand, C.P_stand, C.weight_enum, α, t, h, η, l) diff --git a/src/Classical/linear_code.jl b/src/Classical/linear_code.jl index 1c1014f..0f02922 100644 --- a/src/Classical/linear_code.jl +++ b/src/Classical/linear_code.jl @@ -9,7 +9,7 @@ ############################# """ - LinearCode(G::CTMatrixTypes, parity::Bool=false, brute_force_WE::Bool=true) + LinearCode(G::CTMatrixTypes, parity::Bool = false, brute_force_WE::Bool = true) Return the linear code constructed with generator matrix `G`. If the optional paramater `parity` is set to `true`, a linear code is built with `G` as the parity-check matrix. If the optional parameter @@ -23,10 +23,10 @@ function LinearCode(G::CTMatrixTypes, parity::Bool = false, brute_force_WE::Bool G_new = _remove_empty(G_new, :rows) C = if parity - H = kernel(G_new, side = :right) - rnk_H = rank(H) - # println(rnk_H) - # display(H) + H = kernel(G_new, side = :right) + rnk_H = rank(H) + # println(rnk_H) + # display(H) if ncols(H) == rnk_H H_tr = transpose(H) else @@ -189,45 +189,45 @@ end """ field(C::AbstractLinearCode) -Return the base ring of the generator matrix. +Return the base ring of the generator matrix of `C`. """ field(C::AbstractLinearCode) = C.F """ length(C::AbstractLinearCode) -Return the length of the code. +Return the length of `C`. """ length(C::AbstractLinearCode) = C.n """ dimension(C::AbstractLinearCode) -Return the dimension of the code. +Return the dimension of `C`. """ dimension(C::AbstractLinearCode) = C.k """ cardinality(C::AbstractLinearCode) -Return the cardinality of the code. +Return the cardinality of `C`. """ cardinality(C::AbstractLinearCode) = BigInt(order(C.F))^C.k """ rate(C::AbstractLinearCode) -Return the rate, ``R = k/n``, of the code. +Return the rate of `C`. """ rate(C::AbstractLinearCode) = C.k / C.n """ - generator_matrix(C::AbstractLinearCode, stand_form::Bool=false) + generator_matrix(C::AbstractLinearCode, stand_form::Bool = false) -Return the generator matrix of the code. If the optional parameter `stand_form` +Return the generator matrix of `C`. If the optional parameter `stand_form` is set to `true`, the standard form of the generator matrix is returned instead. """ -function generator_matrix(C::AbstractLinearCode, stand_form::Bool=false) +function generator_matrix(C::AbstractLinearCode, stand_form::Bool = false) if isa(C, QuasiCyclicCode) stand_form && !ismissing(C.G_stand) && (return C.G_stand;) if ismissing(C.G) @@ -259,13 +259,13 @@ function generator_matrix(C::AbstractLinearCode, stand_form::Bool=false) end """ - parity_check_matrix(C::AbstractLinearCode, stand_form::Bool=false) + parity_check_matrix(C::AbstractLinearCode, stand_form::Bool = false) -Return the parity-check matrix of the code. If the optional parameter +Return the parity-check matrix of `C`. If the optional parameter `stand_form` is set to `true`, the standard form of the parity-check matrix is returned instead. """ -function parity_check_matrix(C::AbstractLinearCode, stand_form::Bool=false) +function parity_check_matrix(C::AbstractLinearCode, stand_form::Bool = false) if isa(C, QuasiCyclicCode) if stand_form ismissing(C.H_stand) || (return C.H_stand;) @@ -312,7 +312,7 @@ standard_form_permutation(C::AbstractLinearCode) = C.P_stand """ relative_distance(C::AbstractLinearCode) -Return the relative minimum distance, ``\\delta = d / n`` of the code if ``d`` is known; +Return the relative minimum distance of `C` if ``d`` is known; otherwise return `missing`. """ relative_distance(C::AbstractLinearCode) = C.d / C.n @@ -653,7 +653,7 @@ Singleton_bound(C::AbstractLinearCode) = Singleton_bound(C.n, C.k) """ encode(C::AbstractLinearCode, v::Union{CTMatrixTypes, Vector{Int}}) -Return the encoding of `v` into `C` +Return the encoding of `v` into `C`. """ function encode(C::AbstractLinearCode, v::Union{CTMatrixTypes, Vector{Int}}) w = isa(v, Vector{Int}) ? matrix(C.F, 1, length(v), v) : v @@ -995,9 +995,9 @@ function is_triply_even(C::AbstractLinearCode) end """ - words(C::AbstractLinearCode, only_print::Bool=false) - codewords(C::AbstractLinearCode, only_print::Bool=false) - elements(C::AbstractLinearCode, only_print::Bool=false) + words(C::AbstractLinearCode, only_print::Bool = false) + codewords(C::AbstractLinearCode, only_print::Bool = false) + elements(C::AbstractLinearCode, only_print::Bool = false) Return the elements of `C`. @@ -1005,7 +1005,7 @@ Return the elements of `C`. * If `only_print` is `true`, the elements are only printed to the console and not returned. """ -function words(C::AbstractLinearCode, only_print::Bool=false) +function words(C::AbstractLinearCode, only_print::Bool = false) words = only_print ? nothing : Vector{typeof(C.G)}() G = ismissing(C.P_stand) ? generator_matrix(C, true) : generator_matrix(C, true) * C.P_stand E = base_ring(G) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index e092358..85e4857 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -7,8 +7,8 @@ ############################# # General ############################# -using AllocCheck -using Profile +# using AllocCheck +# using Profile using ProgressMeter """ polynomial(W::WeightEnumerator) @@ -581,26 +581,16 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim #In the main loop we check if lower bound > upper bound before we enumerate and so the lower bounds for the loop use r not r+1 lower_bounds = [information_set_lower_bound(r, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1] - num_thrds = Threads.nthreads() - verbose && println("Detected $num_thrds threads.") - power = Int(floor(log(2, num_thrds))) - work_factor_up_to_log_field_size = 0 predicted_work_factor = fld(n, k) * sum([binomial(k, i) for i in 1:r_term]) verbose && println("Predicted work factor: $predicted_work_factor") if show_progress prog_bar = Progress(predicted_work_factor, dt=1.0, showspeed=true) # updates no faster than once every 1s end - #= TODO rewrite - # c_itr is a binary vector. If it were a _random_ vector the expected - # weight of the subvector c_itr[1:2*uppers[m]] equals uppers[m]. - # So we use the heuristic that w(c_itr[1:2*uppers[m]+c]) for a small constant c is - # usually larger than uppers[m]. - # If not, we compute weight of all of c_itr - =# weight_sum_bound = min(2 * C.u_bound + 5, n) verbose && println("Codeword weights initially checked on first $weight_sum_bound entries") + num_thrds = Threads.nthreads() for r in 2:k if r > 2^16 verbose && println("Reached an r with r>2^16") @@ -632,52 +622,59 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end keep_going = Threads.Atomic{Bool}(true) - num_thrds = 1 p = Int(characteristic(C.F)) uppers = [C.u_bound for _ in 1:num_thrds] + println("uppers[1] $(uppers[1])") founds = [found for _ in 1:num_thrds] exit_thread_indicator_vec = [initial_perm_ind for _ in 1:num_thrds] - # itrs_alloc = @allocated begin - itrs = [SubsetGrayCode(C.k, r)] - # lens = [length(_) for _ in itrs] - # println(lens) - # expected_work_factor += h * sum(lens) + # if num_thrds == 1 + # itrs = [SubsetGrayCode(C.k, r, extended_binomial(C.k, r), 1)] + # else + # itrs = CodingTheory._subset_gray_codes_from_num_threads(C.k, r, num_thrds) # end - # println("alloc itrs used ", itrs_alloc) + itrs = CodingTheory._subset_gray_codes_from_num_threads(C.k, r, num_thrds) + num_thrds = length(itrs) - # Threads.@threads for itr in itrs - vec = vcat(fill(1, r), fill(0, C.k - r)) # initial vec corresponds to the subset {1,..,r} - for m in 1:num_thrds + vec = vcat(fill(1, 1), fill(0, C.k - r), fill(1, r-1)) # initial vec corresponds to the subset {1,..,r} + vec = UInt16.(vec) - itr = itrs[m] + # Threads.@threads for ind in 1:num_thrds + for ind in 1:num_thrds + m = Threads.threadid() # we loop over matrices first so that we can update the lower bound more often (see White Algo 7.1) for i in 1:h - c_itr = zeros(UInt16, C.n) + if keep_going[] == false + break + end + c_itr = zeros(UInt16, C.n) is_first = true curr_mat = A_mats[i] - for u in itr - if !keep_going + itr = itrs[m] + count = 0 + for u in itr + # count += 1 + if keep_going[] == false break end work_factor_up_to_log_field_size += 1 show_progress && ProgressMeter.next!(prog_bar) if r - rank_defs[i] > 0 - for ind in u - - if ind != -1 - is_first = false - @simd for i in eachindex(c_itr) - @inbounds c_itr[i] = xor(c_itr[i], curr_mat[i, ind]) - end - end - end if is_first LinearAlgebra.mul!(c_itr, curr_mat, vec) - @inbounds @simd for j in 1:n + @inbounds @simd for j in eachindex(c_itr) c_itr[j] %= p end + is_first = false + else + for ind in u + if ind != -1 + @simd for i in eachindex(c_itr) + @inbounds c_itr[i] = xor(c_itr[i], curr_mat[i, ind]) + end + end + end end if store_found_codewords @@ -689,12 +686,17 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim if uppers[m] > partial_weight w = sum(c_itr) + if w == 0 + verbose && println("WARNING found weight(c_itr)=0 at count ", count, "u is ", u) + continue + end if uppers[m] > w - uppers[m] = w + uppers[m] = w founds[m] = c_itr exit_thread_indicator_vec[m] = i verbose && println("Adjusting (local) upper bound: $w") if C.l_bound == uppers[m] + println("early exit") Threads.atomic_cas!(keep_going, true, false) else r_term = findfirst(x -> x ≥ C.u_bound, lower_bounds) @@ -706,11 +708,12 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end end end + loc = argmin(uppers) + println("old C.u_bound $(C.u_bound)") + C.u_bound = uppers[loc] + found = founds[loc] + initial_perm_ind = exit_thread_indicator_vec[loc] end - loc = argmin(uppers) - C.u_bound = uppers[loc] - found = founds[loc] - initial_perm_ind = exit_thread_indicator_vec[loc] end show_progress && ProgressMeter.finish!(prog_bar) @@ -721,7 +724,6 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim return C.u_bound, y_Oscar end -#= function _minimum_distance_enumeration_with_matrix_multiply(C::AbstractLinearCode; info_set_alg::Symbol = :auto, verbose::Bool = false, dbg = Dict()) @@ -890,7 +892,6 @@ function _minimum_distance_enumeration_with_matrix_multiply(C::AbstractLinearCod # iszero(C.H * transpose(y)) return C.u_bound, (C.F).(y) end -=# """ words_of_weight(C::AbstractLinearCode, l_bound::Int, u_bound::Int; verbose::Bool = false) diff --git a/src/CodingTheory.jl b/src/CodingTheory.jl index 0f21a34..84b70e2 100644 --- a/src/CodingTheory.jl +++ b/src/CodingTheory.jl @@ -19,6 +19,7 @@ using DataStructures using StatsBase using Distributions using ProgressMeter +using DocStringExtensions import LinearAlgebra: tr, Adjoint, transpose, kron, diagm import Oscar: dual, factor, transpose, order, polynomial, nrows, ncols, degree, @@ -52,6 +53,7 @@ const CTPolyRing = PolyRing{<:CTFieldElem} const CTPolyRingElem = PolyRingElem{<:CTFieldElem} const CTGroupAlgebra = GroupAlgebraElem{fpFieldElem, GroupAlgebra{fpFieldElem, FinGenAbGroup, FinGenAbGroupElem}} const CTChainComplex = Union{ComplexOfMorphisms{AbstractAlgebra.FPModule{fpFieldElem}}} # residue and group algebras later +const CTPolyMatrix = Union{AbstractAlgebra.Generic.MatSpaceElem{fpPolyRingElem}, AbstractAlgebra.Generic.MatSpaceElem{FqPolyRingElem}} include("Classical/types.jl") export AbstractCode, AbstractNonadditiveCode, AbstractNonlinearCode, AbstractAdditiveCode, diff --git a/src/Convolutional/convolutional_code.jl b/src/Convolutional/convolutional_code.jl new file mode 100644 index 0000000..64f75a0 --- /dev/null +++ b/src/Convolutional/convolutional_code.jl @@ -0,0 +1,366 @@ +# Copyright (c) 2024 Eric Sabo +# All rights reserved. +# +# This source code is licensed under the BSD-style license found in the +# LICENSE file in the root directory of this source tree. + +############################# + # constructors +############################# + +""" + ConvolutionalCode(G::CTPolyMatrix, parity::Bool = false) + +Return the convolutional code based on the matrix `G`. + +# Notes +- Currently only accepting matrices with polynomial entries; any code with rational entries is equivalent to a code with polynomial entries by multiplying by the lcm of the denominators. +- Currently only accepting full-rank matrices. +""" +function ConvolutionalCode(G::CTPolyMatrix, parity::Bool = false) + + # TODO check base ring on G is appriorate + D = gen(parent(G_new[1, 1])) + + G_new = deepcopy(G) + G_new = _remove_empty(G_new, :rows) + k = rank(G_new) + k == nrows(G_new) || throw(ArgumentError("The input matrix must have `rank(G)` rows.")) + + H = kernel(G_new, side = :right) + rnk_H = rank(H) + if ncols(H) == rnk_H + H_tr = transpose(H) + else + # remove empty columns for flint objects https://github.com/oscar-system/Oscar.jl/issues/1062 + nr = nrows(H) + H_tr = zero_matrix(base_ring(H), rnk_H, nr) + for r in 1:nr + for c in 1:rnk_H + !iszero(H[r, c]) && (H_tr[c, r] = H[r, c];) + end + end + end + H = H_tr + + if parity + temp = G_new + G_new = H + H = temp + end + + n = ncols(G_new) + vi = zeros(Int, k) + row_degs = zeros(Int, k) + for i in 1:k + row_deg = [degree(G_new[i, c]) for c in 1:n] + row_degs = sum(row_deg) + vi[i] = maximum(row_deg) + end + m = maximum(vi) + + mnrs = minors(G_new, k) + int_deg = maximum([degree(x) for x in mnrs]) + ext_deg = sum(row_degs) + return ConvolutionalCode(base_ring(G_new), n, k, missing, D, m, vi, mnrs, int_deg, ext_deg, G, H, missing) +end + +############################# + # getter functions +############################# + +""" + field(C::AbstractConvolutionalCode) + +Return the base ring of the generator matrix of `C`. +""" +field(C::AbstractConvolutionalCode) = C.F + +""" + length(C::AbstractConvolutionalCode) + +Return the length of `C`. +""" +length(C::AbstractConvolutionalCode) = C.n + +""" + dimension(C::AbstractConvolutionalCode) + +Return the dimension of the code. +""" +dimension(C::AbstractConvolutionalCode) = C.k + +""" + rate(C::AbstractConvolutionalCode) + +Return the rate of `C`. +""" +rate(C::AbstractConvolutionalCode) = C.k / C.n + +""" + parity_check_matrix(C::AbstractConvolutionalCode) + +Return the parity check matrix of `C`. +""" +parity_check_matrix(C::AbstractConvolutionalCode) = C.H + +""" + delay_operator(C::AbstractConvolutionalCode) + +Return the delay operator of `C`. +""" +delay_operator(C::AbstractConvolutionalCode) = C.D + +""" + constraint_lengths(C::AbstractConvolutionalCode) + +Return the constraint lengths of `C`. +""" +constraint_lengths(C::AbstractConvolutionalCode) = C.vi + +""" + overall_constraint_lengths(C::AbstractConvolutionalCode) + +Return the overall constraint length of `C`. +""" +overall_constraint_length(C::AbstractConvolutionalCode) = sum(C.vi) + +""" + memory(C::AbstractConvolutionalCode) + +Return the memory of `C`. +""" +memory(C::AbstractConvolutionalCode) = C.m + +""" + internal_degree(C::AbstractConvolutionalCode) + +Return the internal degree of `C`. +""" +internal_degree(C::AbstractConvolutionalCode) = C.int_deg + +""" + external_degree(C::AbstractConvolutionalCode) + +Return the external degree of `C`. +""" +external_degree(C::AbstractConvolutionalCode) = C.ext_deg + +############################# + # setter functions +############################# + +############################# + # general functions +############################# + +""" + is_realizable(C::AbstractConvolutionalCode) + +Return `true` if the generator matrix of `C` is realizable; otherwise, `false`. + +# Notes +- A generator matrix is realizable if it is a polynomial matrix or all the denominators of the rational functions are equal to one. Since only polynomial matrices are accepted at the moment, this function always returns true. +""" +is_realizable(C::AbstractConvolutionalCode) = true + +""" + is_polynomial(C::AbstractConvolutionalCode) + +Return `true` if the generator matrix of `C` is polynomial; otherwise, `false`. + +# Notes +- This would check all the denominators of the rational functions are equal to one; however, since only polynomial matrices are accepted at the moment, this function always returns true. +""" +is_polynomial(C::AbstractConvolutionalCode) = true + +""" + is_basic(C::AbstractConvolutionalCode) + +Return `true` if the generator matrix of `C` is basic; otherwise, `false`. +""" +function is_basic(C::AbstractConvolutionalCode) + SNF, _, _ = snf_with_transform(C.G) + # diag elements of SNF are invariant factors + # all invariant factors (SNF) are 1 + return all(is_one, [SNF[i, i] for i in 1:k]) +end + +# can always make minimal with row ops +""" + is_minimal(C::AbstractConvolutionalCode) + +Return `true` if the generator matrix of `C` is minimal; otherwise, `false`. +""" +function is_minimal(C::AbstractConvolutionalCode) + is_basic(C) || return false + G_h = zeros(UInt8, C.k, C.n) + for r in 1:C.k + for c in C.n + degree(C.G[r, c]) == C.vi[r] && (G_h[r, c] 1;) + end + end + return k == rank(G_h) +end + +""" + is_reduced(C::AbstractConvolutionalCode) + +Return `true` if the generator matrix of `C` is reduced; otherwise, `false`. +""" +is_reduced(C::AbstractConvolutionalCode) = internal_degree(C) == external_degree(C) + +""" + is_canonical(C::AbstractConvolutionalCode) + +Return `true` if the generator matrix of `C` is canonical; otherwise, `false`. +""" +is_canonical(C::AbstractConvolutionalCode) = is_basic(C) && is_reduced(C) + +""" + is_catastrophic(C::AbstractConvolutionalCode) + +Return `true` if the generator matrix of `C` is catastrophic; otherwise, `false`. +""" +is_catastrophic(C::AbstractConvolutionalCode) = !is_one(sum(ZZ.(coefficients(gcd(C.mnrs))))) + +""" + generator_matrix(C::AbstractConvolutionalCode, systematic::Bool = false) + +Return the generator matrix of `C` if `systematic` is `false` and a systematic generator matrix +otherwise. +""" +generator_matrix(C::AbstractConvolutionalCode, systematic::Bool = false) = + systematic ? (return hnf(fraction_field(base_ring(C.G)).(C.G));) : (return C.G;) + +function terminated_generator_matrix(C::AbstractConvolutionalCode, L::Int) + L ≥ 1 || throw(DomainError(L, "The number of rows must be at least one.")) + + G_mats = [zero_matrix(C.F, C.k, C.n) for _ in 1:C.m] + for r in 1:C.k + for c in 1:C.n + for i in 1:m + G_mats[m][r, c] = coeff(C.G[r, c], i - 1) + end + end + end + + # this matrix has L shifts and truncates at row L + G_L = zero_matrix(C.F, L * C.k, (m + L) * C.n) + for r in 1:L + for c in 1:m + G_L[(r - 1) * C.k + 1:r * C.k, (c - 1 + r - 1) * C.n + 1:(c + r - 1) * C.n] .= G_mats[c] + end + end + return G_L +end + +function truncated_generator_matrix(C::AbstractConvolutionalCode, L::Int) + L ≥ 1 || throw(DomainError(L, "The number of columns must be at least one.")) + + G_mats = [zero_matrix(C.F, C.k, C.n) for _ in 1:C.m] + for r in 1:C.k + for c in 1:C.n + for i in 1:m + G_mats[m][r, c] = coeff(C.G[r, c], i - 1) + end + end + end + + # this matrix truncates at column L + G_L = zeros(C.F, L * C.k, L * C.n) + for offset in 0:min(C.m - 1, L) + for i in 1:L - offset + rows = (1:C.k) .+ ((i - 1) * C.k) + cols = (1:C.n) .+ ((offset + i - 1) * C.n) + G_L[rows, cols] .= G_mats[offset + 1] + end + end + return G_L +end + +function tail_biting_generator_matrix(C::AbstractConvolutionalCode, L::Int) + L ≥ 1 || throw(DomainError(L, "The number of columns must be at least one.")) + + G_mats = [zero_matrix(C.F, C.k, C.n) for _ in 1:C.m] + for r in 1:C.k + for c in 1:C.n + for i in 1:m + G_mats[m][r, c] = coeff(C.G[r, c], i - 1) + end + end + end + + # this matrix has L shifts and truncates at row L + G_L = zero_matrix(C.F, L * C.k, L * C.n) + for r in 1:L + for c in 1:m + if c - 1 + r - 1 ≤ L + if c + r - 2 > L + cols = ((c - 1 + r - 1) % L) * C.n + 1:((c + r - 1) % L) * C.n + else + cols = (c - 1 + r - 1) * C.n + 1:(c + r - 1) * C.n + end + G_L[(r - 1) * C.k + 1:r * C.k, cols] .= G_mats[c] + end + end + end +end + +""" + dual(C::AbstractConvolutionalCode) + +Return the dual code of `C`. +""" +function dual(C::AbstractConvolutionalCode) + # TODO do better in the future, avoid the recomputation + return ConvolutionalCode(C.H) +end + +function show(io::IO, C::AbstractConvolutionalCode) + print(io, "($(C.n), $(C.k)") + !ismissing(C.d) && print(io, ", $(C.d)") + print(io, ")_$(order(C.F)) ") + println(io, "convolutional code") + + if get(io, :compact, true) + println(io, "Memory: $(C.m)") + println(io, "Constraint length: $(overall_constraint_length(C))") + if C.n ≤ 10 + println(io, "Generator matrix: $(C.k) × $(C.n)") + println(io, "\t" * replace(replace(replace(repr(MIME("text/plain"), C.G), + r"\n" => "\n\t"), r"\[" => ""), r"\]" => "")) + end + # if !ismissing(C.weight_enum) + # println(io, "\nComplete weight enumerator:") + # print(io, "\t", polynomial(C.weight_enum)) + # end + end +end + +""" + encode(C::AbstractConvolutionalCode, v::CTPolyMatrix) + +Return the encoding of `v` into `C`. +""" +function encode(C::AbstractConvolutionalCode, v::CTPolyMatrix) + (size(v) != (1, C.k) && size(v) != (C.k, 1)) && + throw(ArgumentError("Vector has incorrect dimension; expected length $(C.k), received: $(size(v)).")) + parent(v) == parent(C.G) || throw(ArgumentError("Vector must have the same parent as the generator matrix.")) + nrows(v) ≠ 1 || return v * C.G + return transpose(v) * C.G +end + +# free distance +# canonical_generator_matrix +# minimum_distance +# column_distance +# distance_profile +# row_distance +# free_distance +# also extended versions +# encode +# interleave::Bool +# scalar and polynomial versions +# syndrome + diff --git a/src/Convolutional/types.jl b/src/Convolutional/types.jl new file mode 100644 index 0000000..a065469 --- /dev/null +++ b/src/Convolutional/types.jl @@ -0,0 +1,43 @@ +# Copyright (c) 2024 Eric Sabo +# All rights reserved. +# +# This source code is licensed under the BSD-style license found in the +# LICENSE file in the root directory of this source tree. + +############################# + # abstract types +############################# + +abstract type AbstractCode end +abstract type AbstractConvolutionalCode <: AbstractCode end + +############################# + # concrete types +############################# + +############################# + # convolutional_code.jl +############################# + +struct PathEnumerator + polynomial::Union{ZZMPolyRingElem, Nemo.AbsSimpleNumFieldElem} + type::Symbol +end + +mutable struct ConvolutionalCode <: AbstractConvolutionalCode + F::CTFieldTypes # base field + n::Int # length + k::Int # dimension + d::Union{Int, Missing} # free distance + # l_bound::Int # lower bound on d + # u_bound::Int # upper bound on d + D::CTFieldElem # delay operator + m::Int # memory + vi::Vector{Int} # constraint lengths + mnrs::Union{Vector{fqPolyRingElem}, Vector{FqPolyRingElem}} + int_deg::Int # interal degree + ext_deg::Int # external degree + G::CTPolyMatrix + H::CTPolyMatrix + path_enum::Union{PathEnumerator, Missing} +end diff --git a/src/LDPC/types.jl b/src/LDPC/types.jl index aabbb46..2f6dbcc 100644 --- a/src/LDPC/types.jl +++ b/src/LDPC/types.jl @@ -49,7 +49,7 @@ mutable struct LDPCCode <: AbstractLDPCCode end ############################# - # LDPC/channels.jl + # LDPC/channels.jl ############################# struct BinaryErasureChannel <: AbstractBinaryErasureChannel diff --git a/src/Quantum/decoders/OTF.jl b/src/Quantum/decoders/OTF.jl index 9ea336d..3a06aaf 100644 --- a/src/Quantum/decoders/OTF.jl +++ b/src/Quantum/decoders/OTF.jl @@ -148,21 +148,31 @@ function _select_erased_columns(H::Matrix{UInt8}, ordered_indices::Vector{Int}, push!(output_indices, col) break end - # elseif count ≥ 1 - # _union_by_rank!(parents, depths, seen_roots_list[col][count], seen_roots_list[col][count + 1]) - # end - # count += 1 end - + if !flag - count = 0 - for row in var_adj_list[col] - if count ≥ 1 - _union_by_rank!(parents, depths, seen_roots_list[col][count], seen_roots_list[col][count + 1]) + # Find maximum of depths[see_roots_list[col][count]], also, see if there are two or more roots of maximum depth. In which case we should grow. + merged_root = -1 + merged_depth = 0 + growth = false + for row_root in seen_roots_list[col] + if depths[row_root] > merged_depth + merged_root = row_root + merged_depth = depths[row_root] + growth = false + elseif depths[row_root] == merged_depth + growth = true end - count += 1 end - # println(seen_roots_list) + + for row_root in seen_roots_list[col] + parents[row_root] = merged_root + end + + if growth + depths[merged_root] += 1 + end + end # println(parents) # println(depths) diff --git a/src/iterators.jl b/src/iterators.jl index 324d5c8..ecc2faf 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -21,7 +21,7 @@ end function _subset_gray_code_from_num_threads(n::Int, k::Int, init_rank::BigInt, num_threads::Int) bin = extended_binomial(n, k) if bin % num_threads != 0 - throw(ArgumentError("num_threads must divide binomial(n k)")) + throw(ArgumentError("num_threads=$num_threads, must divide binomial($n $k)")) end len = fld(bin, num_threads) return SubsetGrayCode(n, k, len, init_rank) @@ -31,13 +31,11 @@ function _subset_gray_codes_from_num_threads(n::Int, k::Int, num_threads::Int) # This function splits a single iterator into several pieces. # The intended usage is to do the iteration with multiple threads bin = extended_binomial(n, k) - if bin % num_threads != 0 - throw(ArgumentError("num_threads must divide binomial(n k)")) - end + i = 0 len = fld(bin, num_threads) itrs = fill(SubsetGrayCode(1,1,1,1), num_threads) - for i in 0:num_threads-1 + for i in 0:num_threads - 1 init_rank = 1 + i * len itrs[i+1] = SubsetGrayCode(n, k, len, init_rank) end @@ -54,7 +52,6 @@ Base.in(v::Vector{Int}, G::SubsetGrayCode) = length(v) == G.k @inline function Base.isdone(G::SubsetGrayCode, state) (_, rank, _) = state - println("called") return rank == G.len end diff --git a/test/LDPC/MP_decoders_test.jl b/test/LDPC/MP_decoders_test.jl index c0cdf1e..3faa7cc 100644 --- a/test/LDPC/MP_decoders_test.jl +++ b/test/LDPC/MP_decoders_test.jl @@ -38,11 +38,11 @@ @test flag == true && out == correct_v flag, out, iter, _ = sum_product(H, v, nm, schedule = :serial); @test flag == true && out == correct_v - flag, out, iter, _ = sum_product(H, v, nm, schedule = :layered); - @test flag == true && out == correct_v - flag, out, iter, _ = sum_product(H, v, nm, schedule = :layered, rand_sched = true); - @test flag == true && out == correct_v - flag, out, iter, _ = sum_product(H, v, nm, erasures = [rand(1:7)]); + # flag, out, iter, _ = sum_product(H, v, nm, schedule = :layered); + # @test flag == true && out == correct_v + # flag, out, iter, _ = sum_product(H, v, nm, schedule = :layered, rand_sched = true); + # @test flag == true && out == correct_v + # flag, out, iter, _ = sum_product(H, v, nm, erasures = [rand(1:7)]); @test flag == true && out == correct_v flag, out, iter, _ = min_sum(H, v, nm, erasures = [rand(1:7)]); @test flag == true && out == correct_v diff --git a/test/Quantum/decoders/OFT_test.jl b/test/Quantum/decoders/OFT_test.jl index 16b7d5e..a86c222 100644 --- a/test/Quantum/decoders/OFT_test.jl +++ b/test/Quantum/decoders/OFT_test.jl @@ -68,7 +68,7 @@ end # check it against the code in this library - for S in [SteaneCode(), Q15RM(), RotatedSurfaceCode(3), RotatedSurfaceCode(5), RotatedSurfaceCode(7)] + for S in [SteaneCode(), Q15RM(), RotatedSurfaceCode(3), RotatedSurfaceCode(5), RotatedSurfaceCode(7), GrossCode()] H_Int = CodingTheory._Flint_matrix_to_Julia_T_matrix(X_stabilizers(S), UInt8); ordering = shuffle(1:S.n); var_adj_list = [Int[] for _ in 1:size(H_Int, 2)]; @@ -81,14 +81,14 @@ end E = sort(CodingTheory._select_erased_columns(H_Int, ordering, var_adj_list)); T = findall(x -> x == true, get_indices(H_Int, ordering)); - # @test E == T - test = E == T - println(test) - if !test - println(ordering) - println(E) - println(T) - end + @test E == T + # test = E == T + # println(test) + # if !test + # println(ordering) + # println(E) + # println(T) + # end end end end diff --git a/test/Quantum/product_codes_test.jl b/test/Quantum/product_codes_test.jl index 68d3f1b..1ed9d91 100644 --- a/test/Quantum/product_codes_test.jl +++ b/test/Quantum/product_codes_test.jl @@ -297,6 +297,7 @@ @test weight_flag # [[400,16,6]] code from Table 1 of https://doi.org/10.1103/PhysRevResearch.2.043423 + F = Oscar.Nemo.Native.GF(2) H = matrix(F, 12, 16, [1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0; 0 0 1 0 0 0 1 1 0 0 0 0 1 0 0 0; From 41ff1cee4c08feb21c62b5303ae98efb71650cb4 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Tue, 28 Jan 2025 12:31:46 -0500 Subject: [PATCH 17/23] Add parallel --- src/Classical/weight_dist.jl | 46 +++++++++++++++++++----------------- 1 file changed, 24 insertions(+), 22 deletions(-) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 85e4857..ffb4511 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -461,7 +461,7 @@ function minimum_distance_Gray(C::AbstractLinearCode; alg::Symbol = :auto, v::Bo end function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zimmermann, - verbose::Bool = false, dbg = Dict(), show_progress=true) + verbose::Bool = false, dbg = Dict(), show_progress=false) ord_F = Int(order(C.F)) ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes.")) @@ -591,6 +591,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim verbose && println("Codeword weights initially checked on first $weight_sum_bound entries") num_thrds = Threads.nthreads() + verbose && println("number of threads ", num_thrds) for r in 2:k if r > 2^16 verbose && println("Reached an r with r>2^16") @@ -600,18 +601,17 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim # (!triply_even_flag && !doubly_even_flag && even_flag) && (C.l_bound += C.l_bound % 2;) # (!triply_even_flag && doubly_even_flag) && (C.l_bound += 4 - C.l_bound % 4;) # triply_even_flag && (C.l_bound += 8 - C.l_bound % 8;) - verbose && println("r: $r") - verbose && println("Lower bound: $(C.l_bound)") - verbose && println("Upper bound: $(C.u_bound)") if C.l_bound >= C.u_bound C.d = C.u_bound y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) # @assert in(y, C) dbg[dbg_key_exit_r] = r - #TODO report actual/expected work factors in dbg dict show_progress && ProgressMeter.finish!(prog_bar) return C.u_bound, y end + verbose && println("r: $r") + verbose && println("Lower bound: $(C.l_bound)") + verbose && println("Upper bound: $(C.u_bound)") if verbose count = 0 @@ -625,36 +625,38 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim p = Int(characteristic(C.F)) uppers = [C.u_bound for _ in 1:num_thrds] - println("uppers[1] $(uppers[1])") founds = [found for _ in 1:num_thrds] exit_thread_indicator_vec = [initial_perm_ind for _ in 1:num_thrds] - # if num_thrds == 1 - # itrs = [SubsetGrayCode(C.k, r, extended_binomial(C.k, r), 1)] - # else - # itrs = CodingTheory._subset_gray_codes_from_num_threads(C.k, r, num_thrds) - # end - itrs = CodingTheory._subset_gray_codes_from_num_threads(C.k, r, num_thrds) - num_thrds = length(itrs) - vec = vcat(fill(1, 1), fill(0, C.k - r), fill(1, r-1)) # initial vec corresponds to the subset {1,..,r} vec = UInt16.(vec) - # Threads.@threads for ind in 1:num_thrds - for ind in 1:num_thrds - m = Threads.threadid() + bin = extended_binomial(C.k, r) + # println("bin is ", bin) + len = fld(bin, num_thrds) + + # Threads.@threads for ind in 1:num_thrds + for ind in 1:num_thrds + # for ind in 1:num_thrds # we loop over matrices first so that we can update the lower bound more often (see White Algo 7.1) for i in 1:h + verbose && println("thread id=$(Threads.threadid())\t mat id=$(i)") if keep_going[] == false break end c_itr = zeros(UInt16, C.n) is_first = true curr_mat = A_mats[i] - itr = itrs[m] + init_rank = 1 + (ind - 1) * len + if ind == num_thrds + len = bin - (num_thrds - 1) * len + end count = 0 - for u in itr - # count += 1 + for u in SubsetGrayCode(C.k, r, len, init_rank) + m = Threads.threadid() + # println("Comparing sum with $(uppers[m]). Also C.u=$(C.u_bound)") + + count += 1 if keep_going[] == false break end @@ -708,8 +710,8 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end end end - loc = argmin(uppers) - println("old C.u_bound $(C.u_bound)") + loc = argmin(uppers) #TODO not threadsafe ? + # verbose && println("old C.u_bound $(C.u_bound)") C.u_bound = uppers[loc] found = founds[loc] initial_perm_ind = exit_thread_indicator_vec[loc] From f931d1f1dd6be1f9ee5f49c7f83b3118cdb0f67c Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Thu, 13 Feb 2025 11:56:05 -0500 Subject: [PATCH 18/23] Init vecs from unrank function --- Manifest.toml | 23 +++++- Project.toml | 8 +++ src/Classical/weight_dist.jl | 131 ++++++++++++++++++++++------------- src/CodingTheory.jl | 2 + src/iterators.jl | 59 +++++++++++++--- src/utils.jl | 1 + 6 files changed, 166 insertions(+), 58 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 0dcfcda..2e03a2c 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.11.3" manifest_format = "2.0" -project_hash = "bbb0e5646f1748a2c6f584642f98e12c63222479" +project_hash = "6b1804c68d878bc0c57ff542fcb5042f0144569a" [[deps.ASL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -65,6 +65,12 @@ version = "2.2.0" uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" version = "1.11.0" +[[deps.BenchmarkTools]] +deps = ["Compat", "JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"] +git-tree-sha1 = "e38fbc49a620f5d0b660d7f543db1009fe0f8336" +uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +version = "1.6.0" + [[deps.BinaryWrappers]] deps = ["JLLWrappers", "Scratch"] git-tree-sha1 = "7fea8f658689fa5062b23f4400eda888b7ae2aaa" @@ -427,6 +433,12 @@ deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" version = "2.28.6+0" +[[deps.Memoization]] +deps = ["MacroTools"] +git-tree-sha1 = "7dbf904fa6c4447bd1f1d316886bfbe29feacf45" +uuid = "6fafb56a-5788-4b4e-91ca-c0cea6611c73" +version = "0.2.2" + [[deps.Missings]] deps = ["DataAPI"] git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d" @@ -576,6 +588,10 @@ deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" version = "1.11.0" +[[deps.Profile]] +uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" +version = "1.11.0" + [[deps.ProgressMeter]] deps = ["Distributed", "Printf"] git-tree-sha1 = "8f6bc219586aef8baf0ff9a5fe16ee9c70cb65e4" @@ -820,6 +836,11 @@ git-tree-sha1 = "42fd9023fef18b9b78c8343a4e2f3813ffbcefcb" uuid = "1c621080-faea-4a02-84b6-bbd5e436b8fe" version = "1.0.0" +[[deps.ThreadSafeDicts]] +git-tree-sha1 = "300b753c0a786ea43fdafc26a4e50b87fb58cd50" +uuid = "4239201d-c60e-5e0a-9702-85d713665ba7" +version = "0.1.6" + [[deps.UUIDs]] deps = ["Random", "SHA"] uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" diff --git a/Project.toml b/Project.toml index 359b24e..69ae3dd 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,9 @@ authors = ["Eric Sabo ", "Ben Ide "] version = "0.1.0" [deps] +AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d" AutoHashEquals = "15f4f7f2-30c1-5605-9d31-71845cf9641f" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" @@ -13,6 +15,7 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Memoization = "6fafb56a-5788-4b4e-91ca-c0cea6611c73" Oscar = "f1435218-dba5-11e9-1e4d-f1a5fab5fc13" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -20,6 +23,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a" +ThreadSafeDicts = "4239201d-c60e-5e0a-9702-85d713665ba7" [weakdeps] CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" @@ -40,7 +44,9 @@ JuMPExt = ["JuMP", "GLPK"] MakieExt = ["Makie"] [compat] +AbstractAlgebra = "0.43.12" AutoHashEquals = "2.2.0" +BenchmarkTools = "1.6.0" CairoMakie = "0.13.1" Colors = "0.13.0" Combinatorics = "1.0.2" @@ -52,9 +58,11 @@ GraphMakie = "0.5.13" GraphPlot = "0.6.0" Graphs = "1.12.0" Makie = "0.19.11" +Memoization = "0.2.2" NetworkLayout = "0.4.9" Oscar = "1.2.2" StatsBase = "0.34.4" +ThreadSafeDicts = "0.1.6" WGLMakie = "0.8.15" julia = "≥ 1.10" diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index ffb4511..446b902 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -481,6 +481,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim A_mats, perms_mats, rnks = information_sets(G, info_set_alg, permute = true, only_A = false) A_mats = [deepcopy(_Flint_matrix_to_Julia_T_matrix(Ai, UInt16)') for Ai in A_mats] + # A_mats_trunc = () perms_mats = [deepcopy(_Flint_matrix_to_Julia_T_matrix(Pi, UInt16)') for Pi in perms_mats] h = length(A_mats) # println("Starting loop to refine upper bound. Initial upper bound ", C.u_bound, " num of mats is ", length(A_mats), " dimension ", size(A_mats[1])) @@ -502,6 +503,11 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end k, n = size(G) + A_mats_trunc = [Matrix{UInt16}(undef, k, n-k) for _ in 1:length(A_mats)] + for i in 1:size(A_mats, 1) + A_mats_trunc[i] = deepcopy(A_mats[i][k+1 : n, :]) + end + if info_set_alg == :Brouwer && rnks[h] != k println("Rank of last matrix too small") return @@ -537,7 +543,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim # can make this faster with dots and views if store_found_codewords for a in 1:size(g, 2) - output_vec = perms_mats[j] * g[:, a] + output_vec = Int.(perms_mats[j] * g[:, a]) # @assert in(output_vec, C) push!(dbg["found_codewords"], (1, [], output_vec)) end @@ -587,14 +593,14 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim if show_progress prog_bar = Progress(predicted_work_factor, dt=1.0, showspeed=true) # updates no faster than once every 1s end - weight_sum_bound = min(2 * C.u_bound + 5, n) + weight_sum_bound = min(2 * C.u_bound + 5, n-k) verbose && println("Codeword weights initially checked on first $weight_sum_bound entries") num_thrds = Threads.nthreads() verbose && println("number of threads ", num_thrds) for r in 2:k if r > 2^16 - verbose && println("Reached an r with r>2^16") + verbose && println("Reached an global_r with r>2^16") end C.l_bound < lower_bounds[r] && (C.l_bound = lower_bounds[r];) # an even code can't have have an odd minimum weight @@ -620,84 +626,112 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end count > 0 && println("$count of the original $h information sets no longer contribute to the lower bound") end - - keep_going = Threads.Atomic{Bool}(true) - p = Int(characteristic(C.F)) + uppers = [C.u_bound for _ in 1:num_thrds] founds = [found for _ in 1:num_thrds] exit_thread_indicator_vec = [initial_perm_ind for _ in 1:num_thrds] + keep_going = Threads.Atomic{Bool}(true) + # my_lock = ReentrantLock(); + # Threads.@threads for ind in 1:num_thrds - vec = vcat(fill(1, 1), fill(0, C.k - r), fill(1, r-1)) # initial vec corresponds to the subset {1,..,r} - vec = UInt16.(vec) + #TODO do the following 3 lines together for 1 thread testing + # num_thrds = 1 bin = extended_binomial(C.k, r) - # println("bin is ", bin) - len = fld(bin, num_thrds) - # Threads.@threads for ind in 1:num_thrds - for ind in 1:num_thrds + Threads.@threads for ind in 1:num_thrds + # num_thrds = 1 # for ind in 1:num_thrds + # local r = global_r + thrd_stop_msg = "Stopping current thread, main loop finished" + # a = @allocated begin + # end + # println("binomial() allocs $(a)") + + # set len + len = (ind == num_thrds) ? bin - (num_thrds - 1) * fld(bin, num_thrds) : fld(bin, num_thrds) + + # set first_vec + init_rank = 1 + (ind - 1) * fld(bin, num_thrds) + first_vec = zeros(Int, k) + if init_rank == 1 + for i in 1:r + first_vec[i] = 1 + end + else + CodingTheory._subset_unrank_to_vec!(big(init_rank), UInt64(r), first_vec) + # verbose && println("unranked with args rank=$(init_rank), r=$(r), got first_vec=$(first_vec)") + end + # we loop over matrices first so that we can update the lower bound more often (see White Algo 7.1) for i in 1:h - verbose && println("thread id=$(Threads.threadid())\t mat id=$(i)") + # verbose && println("mat id=$(i), thrd=$(num_thrds)") if keep_going[] == false + verbose && println(thrd_stop_msg) break end - c_itr = zeros(UInt16, C.n) + + c_itr = zeros(UInt16, C.n - C.k) is_first = true - curr_mat = A_mats[i] - init_rank = 1 + (ind - 1) * len - if ind == num_thrds - len = bin - (num_thrds - 1) * len - end - count = 0 - for u in SubsetGrayCode(C.k, r, len, init_rank) - m = Threads.threadid() - # println("Comparing sum with $(uppers[m]). Also C.u=$(C.u_bound)") - - count += 1 + curr_mat = A_mats_trunc[i] + # local count = 0 + # verbose && println("SubsetGrayCode with len=$(len), rank=$(init_rank)") + for u in SubsetGrayCode(k, r, len, init_rank) + # if (count + init_rank) > bin + # println("Error: current total rank=$(count + init_rank), length of all iterators=$(bin), index=$(ind), init_rank given to iterator is $(init_rank)") + # break + # end if keep_going[] == false + println(thrd_stop_msg) break end work_factor_up_to_log_field_size += 1 show_progress && ProgressMeter.next!(prog_bar) if r - rank_defs[i] > 0 if is_first - LinearAlgebra.mul!(c_itr, curr_mat, vec) + LinearAlgebra.mul!(c_itr, curr_mat, first_vec) @inbounds @simd for j in eachindex(c_itr) c_itr[j] %= p end is_first = false else - for ind in u - if ind != -1 + for ci in u + if ci != -1 @simd for i in eachindex(c_itr) - @inbounds c_itr[i] = xor(c_itr[i], curr_mat[i, ind]) + @inbounds c_itr[i] = xor(c_itr[i], curr_mat[i, ci]) end end end end - if store_found_codewords - output_vec = perms_mats[initial_perm_ind] * c_itr - push!(dbg["found_codewords"], (r, [], output_vec)) - end + # if store_found_codewords + # subset_vec_full = zeros(Int, k) + # verbose && println("unrank with init_rank=$(init_rank) and count=$(count)") + # CodingTheory._subset_unrank_to_vec!(big(init_rank + count), UInt64(r), subset_vec_full) + # output_vec = Int.(perms_mats[exit_thread_indicator_vec[ind]] * vcat(subset_vec_full, c_itr)) + # push!(dbg["found_codewords"], (r, [], Int.(output_vec))) + # end - partial_weight = sum(view(c_itr, 1:weight_sum_bound)) + partial_weight = r + sum(view(c_itr, 1:weight_sum_bound)) - if uppers[m] > partial_weight - w = sum(c_itr) + if uppers[ind] > partial_weight + w = r + sum(c_itr) if w == 0 - verbose && println("WARNING found weight(c_itr)=0 at count ", count, "u is ", u) + verbose && println("WARNING found weight(c_itr)=0 at count\t", count, "u is ", u, "c_itr_full\n", c_itr_full, "vec\n", first_vec_subset, "c_itr\n", c_itr) continue end - if uppers[m] > w - uppers[m] = w - founds[m] = c_itr - exit_thread_indicator_vec[m] = i - verbose && println("Adjusting (local) upper bound: $w") - if C.l_bound == uppers[m] + if uppers[ind] > w + subset_vec_full = zeros(Int, k) + CodingTheory._subset_unrank_to_vec!(big(init_rank + count), UInt64(r), subset_vec_full) + + uppers[ind] = w + founds[ind] = vcat(subset_vec_full, c_itr) + verbose && @assert size(founds[ind], 1) == C.n "found vector has length $(size(founds[ind], 1)) but should be n=$(C.n)" + exit_thread_indicator_vec[ind] = i + + verbose && println("Adjusting (local) upper bound: $w for c_itr=$(c_itr)") + if C.l_bound == uppers[ind] println("early exit") Threads.atomic_cas!(keep_going, true, false) else @@ -708,9 +742,10 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end end end + # count += 1 end end - loc = argmin(uppers) #TODO not threadsafe ? + loc = argmin(uppers) # verbose && println("old C.u_bound $(C.u_bound)") C.u_bound = uppers[loc] found = founds[loc] @@ -719,11 +754,11 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end show_progress && ProgressMeter.finish!(prog_bar) - # at this point we are guaranteed to have found the answer + # at this point we are guaranteed to have found the code's distance C.d = C.u_bound - y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) - # iszero(C.H * transpose(y)) - return C.u_bound, y_Oscar + y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) # weight(y) >= C.d with equality not being the typical case + verbose && @assert iszero(C.H * transpose(y)) + return C.u_bound, y end function _minimum_distance_enumeration_with_matrix_multiply(C::AbstractLinearCode; info_set_alg::Symbol = :auto, diff --git a/src/CodingTheory.jl b/src/CodingTheory.jl index 746e530..3120e09 100644 --- a/src/CodingTheory.jl +++ b/src/CodingTheory.jl @@ -20,6 +20,8 @@ using StatsBase using Distributions using ProgressMeter using DocStringExtensions +using Memoization +using ThreadSafeDicts import LinearAlgebra: tr, Adjoint, transpose, kron, diagm import Oscar: dual, factor, transpose, order, polynomial, nrows, ncols, degree, diff --git a/src/iterators.jl b/src/iterators.jl index ecc2faf..05da5d0 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -76,20 +76,34 @@ end Note that inds is part of the iterator's state inds only to prevent reallocation in each iteration. =# + # subset_vec = zeros(Int, G.k) + # CodingTheory._subset_unrank!(G.init_rank, UInt(G.n), subset_vec) + + # a = @allocated begyyin subset_vec = zeros(Int, G.k) - CodingTheory._subset_unrank!(G.init_rank, UInt(G.n), subset_vec) + if G.init_rank == 1 + subset_vec = Int.(collect(1 : G.k)) + # first_vec = vcat(ones(Int, r), zeros(Int, k - r)) + # for i in 1:G.k + # subset_vec[i] = 1 + # end + else + CodingTheory._subset_unrank!(G.init_rank, UInt(G.n), subset_vec) + end + # end; a > 0 && @show a + inds = fill(-1, 3) (inds, (subset_vec, 1, inds)) end @inline function Base.iterate(G::SubsetGrayCode, state) # Based on Algorithm 2.13 in kreher1999combinatorial - v, index, inds = state + v, rank, inds = state - if index == G.len + if rank == G.len return nothing end - index += 1 + rank += 1 @inbounds begin inds[1] = -1 @@ -142,7 +156,7 @@ end end end end - return (inds, (v, index, inds)) + return (inds, (v, rank, inds)) end function _update_indices!(indices::Vector{Int}, x::Int, y::Int) @@ -184,31 +198,58 @@ function _subset_rank(v::Vector{UInt}, k::UInt) if (k % 2) == 1 r = r - 1 end + r = r + 1 # we use 1-indexed ranks to follow the usual Julia convention return r end +function _subset_unrank_to_vec!(r::BigInt, k::UInt, vec::Vector{Int}) + # returns a {0,1} vector with nonzero entries corresponding to the subset returned by _subset_unrank! + n = length(vec) + subset_vec = Int.(zeros(k)) + CodingTheory._subset_unrank!(r, UInt(n), subset_vec) + for i in subset_vec + vec[i] = 1 + end +end + function _subset_unrank!(r::BigInt, n::UInt, T::Vector{Int}) # Based on Algorithm 2.12 in kreher1999combinatorial + r = r - 1 + k = length(T) subset_size_str::String = "subset size k = $k must be smaller than the set size n = $n" k > n && throw(ArgumentError(subset_size_str)) - bnd = binomial(n, k) - rank_size_str::String = "rank must be in [0, choose(n, k) - 1] = $bnd" - r > bnd && throw(ArgumentError(rank_size_str)) + # count_bins = 0 + # bnd = extended_binomial(n, k) + # # count_bins += 1 + # rank_size_str::String = "rank must be in [0, choose(n, k) - 1] = $bnd but is $(r)" + # # r > bnd && throw(ArgumentError(rank_size_str)) + # if r > bnd + # println(rank_size_str) + # println("bug exit _subset_unrank") + # return + # end x = 0 i = 0 - y = 0 + y = 0 x = n for i::UInt in k:-1:1 y = extended_binomial(x, i) + # println((Int16.(x),Int16.(i))) + # count_bins += 1 while y > r x = x - 1 y = extended_binomial(x, i) + # println((x,i)) + # count_bins += 1 end T[i] = x + 1 r = extended_binomial(x + 1, i) - r - 1 + # println((x+1,i)) + # count_bins += 1 end + # println("num calls to binomial(a,b) in 1 unrank = $(count_bins)") end struct GrayCode diff --git a/src/utils.jl b/src/utils.jl index 5c91c38..b2af121 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1993,6 +1993,7 @@ function _rand_invertible_matrix(F::CTFieldTypes, n::Integer) return A end +# @memoize ThreadSafeDict function extended_binomial(x::Union{Int, UInt}, y::Union{Int, UInt}) function extended_binomial(x::Union{Int, UInt}, y::Union{Int, UInt}) z = big(0) if y <= x From b836daa3a66a32e8f18e1b087ab1fe24ba189ab0 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Thu, 20 Feb 2025 15:21:04 -0500 Subject: [PATCH 19/23] Fix race and optimize --- src/Classical/weight_dist.jl | 71 ++++++++++-------------------------- src/iterators.jl | 14 +++---- src/utils.jl | 6 +-- 3 files changed, 27 insertions(+), 64 deletions(-) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 446b902..3f4ed1f 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -473,12 +473,10 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim #TODO err string can instruct the user to construct a new code without 0 cols and tell them the function for that throw(ArgumentError("Codes with standard form of generator matrix having 0 columns not supported")) end - G = C.G #TODO remove local var G - # generate if not pre-stored parity_check_matrix(C) - A_mats, perms_mats, rnks = information_sets(G, info_set_alg, permute = true, only_A = false) + A_mats, perms_mats, rnks = information_sets(C.G, info_set_alg, permute = true, only_A = false) A_mats = [deepcopy(_Flint_matrix_to_Julia_T_matrix(Ai, UInt16)') for Ai in A_mats] # A_mats_trunc = () @@ -499,10 +497,10 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end store_found_codewords = haskey(dbg, dbg_key_found) if store_found_codewords - verbose && println("dbg Dict: Storing the codewords in the debug dictionary, key=$dbg_key_found") + verbose && println("dbg Dict: Storing the codewords in the debug dictionary @key=$dbg_key_found") end - k, n = size(G) + k, n = size(C.G) A_mats_trunc = [Matrix{UInt16}(undef, k, n-k) for _ in 1:length(A_mats)] for i in 1:size(A_mats, 1) A_mats_trunc[i] = deepcopy(A_mats[i][k+1 : n, :]) @@ -544,11 +542,9 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim if store_found_codewords for a in 1:size(g, 2) output_vec = Int.(perms_mats[j] * g[:, a]) - # @assert in(output_vec, C) push!(dbg["found_codewords"], (1, [], output_vec)) end end - # println(dbg["found_codewords"]) w, i = _min_wt_col(g) if w <= C.u_bound found = g[:, i] @@ -563,7 +559,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim l = 0 if verbose - _, _, b_rnks = information_sets(G, :Brouwer, permute = true, only_A = false) + _, _, b_rnks = information_sets(C.G, :Brouwer, permute = true, only_A = false) b_h = length(b_rnks) b_lower_bounds = [information_set_lower_bound(r+1, n, k, l, [k - 0 for i in 1:b_h], :Brouwer, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1] b_r_term = findfirst(x -> x ≥ C.u_bound, b_lower_bounds) @@ -587,7 +583,6 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim #In the main loop we check if lower bound > upper bound before we enumerate and so the lower bounds for the loop use r not r+1 lower_bounds = [information_set_lower_bound(r, n, k, l, rank_defs, info_set_alg, even = even_flag, doubly_even = doubly_even_flag, triply_even = triply_even_flag) for r in 1:k-1] - work_factor_up_to_log_field_size = 0 predicted_work_factor = fld(n, k) * sum([binomial(k, i) for i in 1:r_term]) verbose && println("Predicted work factor: $predicted_work_factor") if show_progress @@ -597,10 +592,10 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim verbose && println("Codeword weights initially checked on first $weight_sum_bound entries") num_thrds = Threads.nthreads() - verbose && println("number of threads ", num_thrds) + verbose && println("Number of threads ", num_thrds) for r in 2:k if r > 2^16 - verbose && println("Reached an global_r with r>2^16") + verbose && println("Reached an r larger than 2^16") end C.l_bound < lower_bounds[r] && (C.l_bound = lower_bounds[r];) # an even code can't have have an odd minimum weight @@ -610,7 +605,6 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim if C.l_bound >= C.u_bound C.d = C.u_bound y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) - # @assert in(y, C) dbg[dbg_key_exit_r] = r show_progress && ProgressMeter.finish!(prog_bar) return C.u_bound, y @@ -632,27 +626,15 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim founds = [found for _ in 1:num_thrds] exit_thread_indicator_vec = [initial_perm_ind for _ in 1:num_thrds] keep_going = Threads.Atomic{Bool}(true) - # my_lock = ReentrantLock(); - # Threads.@threads for ind in 1:num_thrds - - #TODO do the following 3 lines together for 1 thread testing - # num_thrds = 1 bin = extended_binomial(C.k, r) + thrd_stop_msg = "Stopping current thread, main loop finished" + Threads.@threads for ind in 1:num_thrds - # num_thrds = 1 - # for ind in 1:num_thrds - # local r = global_r - thrd_stop_msg = "Stopping current thread, main loop finished" - # a = @allocated begin - # end - # println("binomial() allocs $(a)") - - # set len len = (ind == num_thrds) ? bin - (num_thrds - 1) * fld(bin, num_thrds) : fld(bin, num_thrds) - # set first_vec + # iteration begins with a single matrix multiplication of the generator matrix by first_vec init_rank = 1 + (ind - 1) * fld(bin, num_thrds) first_vec = zeros(Int, k) if init_rank == 1 @@ -660,13 +642,11 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim first_vec[i] = 1 end else - CodingTheory._subset_unrank_to_vec!(big(init_rank), UInt64(r), first_vec) - # verbose && println("unranked with args rank=$(init_rank), r=$(r), got first_vec=$(first_vec)") + CodingTheory._subset_unrank_to_vec!(init_rank, UInt64(r), first_vec) end - # we loop over matrices first so that we can update the lower bound more often (see White Algo 7.1) + # as in White Algo 7.1 we loop over matrices first so that we can update the lower bound more often for i in 1:h - # verbose && println("mat id=$(i), thrd=$(num_thrds)") if keep_going[] == false verbose && println(thrd_stop_msg) break @@ -675,18 +655,11 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim c_itr = zeros(UInt16, C.n - C.k) is_first = true curr_mat = A_mats_trunc[i] - # local count = 0 - # verbose && println("SubsetGrayCode with len=$(len), rank=$(init_rank)") for u in SubsetGrayCode(k, r, len, init_rank) - # if (count + init_rank) > bin - # println("Error: current total rank=$(count + init_rank), length of all iterators=$(bin), index=$(ind), init_rank given to iterator is $(init_rank)") - # break - # end if keep_going[] == false println(thrd_stop_msg) break end - work_factor_up_to_log_field_size += 1 show_progress && ProgressMeter.next!(prog_bar) if r - rank_defs[i] > 0 if is_first @@ -705,25 +678,21 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end end - # if store_found_codewords - # subset_vec_full = zeros(Int, k) - # verbose && println("unrank with init_rank=$(init_rank) and count=$(count)") - # CodingTheory._subset_unrank_to_vec!(big(init_rank + count), UInt64(r), subset_vec_full) - # output_vec = Int.(perms_mats[exit_thread_indicator_vec[ind]] * vcat(subset_vec_full, c_itr)) - # push!(dbg["found_codewords"], (r, [], Int.(output_vec))) - # end + if store_found_codewords + subset_vec_full = zeros(Int, k) + CodingTheory._subset_unrank_to_vec!(UInt128(init_rank + count), UInt64(r), subset_vec_full) + output_vec = Int.(perms_mats[exit_thread_indicator_vec[ind]] * vcat(subset_vec_full, c_itr)) + push!(dbg["found_codewords"], (r, [], Int.(output_vec))) + end partial_weight = r + sum(view(c_itr, 1:weight_sum_bound)) if uppers[ind] > partial_weight w = r + sum(c_itr) - if w == 0 - verbose && println("WARNING found weight(c_itr)=0 at count\t", count, "u is ", u, "c_itr_full\n", c_itr_full, "vec\n", first_vec_subset, "c_itr\n", c_itr) - continue - end + verbose && @assert w == 0 if uppers[ind] > w subset_vec_full = zeros(Int, k) - CodingTheory._subset_unrank_to_vec!(big(init_rank + count), UInt64(r), subset_vec_full) + CodingTheory._subset_unrank_to_vec!(UInt128(init_rank + count), UInt64(r), subset_vec_full) uppers[ind] = w founds[ind] = vcat(subset_vec_full, c_itr) @@ -742,11 +711,9 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end end end - # count += 1 end end loc = argmin(uppers) - # verbose && println("old C.u_bound $(C.u_bound)") C.u_bound = uppers[loc] found = founds[loc] initial_perm_ind = exit_thread_indicator_vec[loc] diff --git a/src/iterators.jl b/src/iterators.jl index 05da5d0..4c6d106 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -8,8 +8,8 @@ struct SubsetGrayCode # details about the iterator using this struct are in the comments in the iterate() function below n::Int # length of codewords k::Int # weight of codewords - len::BigInt # length of the iterator - init_rank::BigInt # iteration will start with the subset of this rank + len::UInt128 # length of the iterator + init_rank::UInt128 # iteration will start with the subset of this rank end function SubsetGrayCode(n::Int, k::Int) @@ -18,7 +18,7 @@ function SubsetGrayCode(n::Int, k::Int) return _subset_gray_code_from_num_threads(n, k, init_rank, num_threads) end -function _subset_gray_code_from_num_threads(n::Int, k::Int, init_rank::BigInt, num_threads::Int) +function _subset_gray_code_from_num_threads(n::Int, k::Int, init_rank::UInt128, num_threads::Int) bin = extended_binomial(n, k) if bin % num_threads != 0 throw(ArgumentError("num_threads=$num_threads, must divide binomial($n $k)")) @@ -189,8 +189,8 @@ function _subset_rank(v::Vector{UInt}, k::UInt) Based on Algorithm 2.11 in kreher1999combinatorial Results are undefined if the entries of v arent in {1,..,n} for n>=k =# - r = BigInt(0) - s = BigInt(1) + r = UInt128(0) + s = UInt128(1) for i in k:-1:1 r = r + extended_binomial(v[i], i) * s s = -s @@ -202,7 +202,7 @@ function _subset_rank(v::Vector{UInt}, k::UInt) return r end -function _subset_unrank_to_vec!(r::BigInt, k::UInt, vec::Vector{Int}) +function _subset_unrank_to_vec!(r::UInt128, k::UInt, vec::Vector{Int}) # returns a {0,1} vector with nonzero entries corresponding to the subset returned by _subset_unrank! n = length(vec) subset_vec = Int.(zeros(k)) @@ -212,7 +212,7 @@ function _subset_unrank_to_vec!(r::BigInt, k::UInt, vec::Vector{Int}) end end -function _subset_unrank!(r::BigInt, n::UInt, T::Vector{Int}) +function _subset_unrank!(r::UInt128, n::UInt, T::Vector{Int}) # Based on Algorithm 2.12 in kreher1999combinatorial r = r - 1 diff --git a/src/utils.jl b/src/utils.jl index b2af121..7ce029a 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1995,11 +1995,7 @@ end # @memoize ThreadSafeDict function extended_binomial(x::Union{Int, UInt}, y::Union{Int, UInt}) function extended_binomial(x::Union{Int, UInt}, y::Union{Int, UInt}) - z = big(0) - if y <= x - z = binomial(big(x), big(y)) - end - return z + return y <= x ? UInt128.(binomial(x, y)) : UInt128(0) end function _value_distribution(vals) From 7bad8f8503542303d0d77179d014e1bee32da937 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Thu, 20 Feb 2025 16:22:19 -0500 Subject: [PATCH 20/23] Clean --- benchmarks/weight_dist_bench.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/benchmarks/weight_dist_bench.jl b/benchmarks/weight_dist_bench.jl index ba5759a..a10f8c7 100644 --- a/benchmarks/weight_dist_bench.jl +++ b/benchmarks/weight_dist_bench.jl @@ -4,12 +4,10 @@ using LinearAlgebra using BenchmarkTools using Random using Dates -# using AllocCheck -# using Profile include("benchmark_utils.jl") CodingTheory.Random.seed!(0) -function white_n70_k35(use_mine, terse_description) +function white_n70_k35(use_new, terse_description) F = Oscar.Nemo.Native.GF(2) mat = matrix(F, [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 0 0 1 1 0 1 1 1 0 1 0 1 1 1 1 0 1 1 0 1 1 0 0 1 0 0 1 1 0 1; 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 1 1 1 0 1 1 0 1 0 0 1 1 0 1 0 1 1 0 1 1 0 1 1 1 0 0 0 1 0; @@ -49,7 +47,7 @@ function white_n70_k35(use_mine, terse_description) num_thrds = Threads.nthreads() println("Benchmarking with $num_thrds threads.") - if use_mine + if use_new function_description = "_minimum_distance_BZ(codeN70K35,:Brouwer;false)" terse_description ? println(function_description) : println(benchmark_description(function_description)) return @benchmark CodingTheory._minimum_distance_BZ(codeN70K35, info_set_alg = :Brouwer; verbose = false) setup = (codeN70K35=LinearCode($mat, false, false)) @@ -62,5 +60,4 @@ end use_new=false t = white_n70_k35(use_new, true) -# t = test_white_105(use_mine) t \ No newline at end of file From 1582bba8fdfe343e1f94d0df0a92c53daa1fa5e9 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Thu, 20 Feb 2025 16:59:44 -0500 Subject: [PATCH 21/23] Clean --- Manifest.toml | 23 +------ Project.toml | 8 --- ext/MakieExt/LDPC/codes.jl | 125 ------------------------------------- src/CodingTheory.jl | 2 - src/iterators.jl | 58 +---------------- src/utils.jl | 2 - 6 files changed, 2 insertions(+), 216 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 2e03a2c..0dcfcda 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.11.3" manifest_format = "2.0" -project_hash = "6b1804c68d878bc0c57ff542fcb5042f0144569a" +project_hash = "bbb0e5646f1748a2c6f584642f98e12c63222479" [[deps.ASL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -65,12 +65,6 @@ version = "2.2.0" uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" version = "1.11.0" -[[deps.BenchmarkTools]] -deps = ["Compat", "JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"] -git-tree-sha1 = "e38fbc49a620f5d0b660d7f543db1009fe0f8336" -uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -version = "1.6.0" - [[deps.BinaryWrappers]] deps = ["JLLWrappers", "Scratch"] git-tree-sha1 = "7fea8f658689fa5062b23f4400eda888b7ae2aaa" @@ -433,12 +427,6 @@ deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" version = "2.28.6+0" -[[deps.Memoization]] -deps = ["MacroTools"] -git-tree-sha1 = "7dbf904fa6c4447bd1f1d316886bfbe29feacf45" -uuid = "6fafb56a-5788-4b4e-91ca-c0cea6611c73" -version = "0.2.2" - [[deps.Missings]] deps = ["DataAPI"] git-tree-sha1 = "ec4f7fbeab05d7747bdf98eb74d130a2a2ed298d" @@ -588,10 +576,6 @@ deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" version = "1.11.0" -[[deps.Profile]] -uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" -version = "1.11.0" - [[deps.ProgressMeter]] deps = ["Distributed", "Printf"] git-tree-sha1 = "8f6bc219586aef8baf0ff9a5fe16ee9c70cb65e4" @@ -836,11 +820,6 @@ git-tree-sha1 = "42fd9023fef18b9b78c8343a4e2f3813ffbcefcb" uuid = "1c621080-faea-4a02-84b6-bbd5e436b8fe" version = "1.0.0" -[[deps.ThreadSafeDicts]] -git-tree-sha1 = "300b753c0a786ea43fdafc26a4e50b87fb58cd50" -uuid = "4239201d-c60e-5e0a-9702-85d713665ba7" -version = "0.1.6" - [[deps.UUIDs]] deps = ["Random", "SHA"] uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" diff --git a/Project.toml b/Project.toml index 69ae3dd..359b24e 100644 --- a/Project.toml +++ b/Project.toml @@ -4,9 +4,7 @@ authors = ["Eric Sabo ", "Ben Ide "] version = "0.1.0" [deps] -AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d" AutoHashEquals = "15f4f7f2-30c1-5605-9d31-71845cf9641f" -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" @@ -15,7 +13,6 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Memoization = "6fafb56a-5788-4b4e-91ca-c0cea6611c73" Oscar = "f1435218-dba5-11e9-1e4d-f1a5fab5fc13" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -23,7 +20,6 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a" -ThreadSafeDicts = "4239201d-c60e-5e0a-9702-85d713665ba7" [weakdeps] CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" @@ -44,9 +40,7 @@ JuMPExt = ["JuMP", "GLPK"] MakieExt = ["Makie"] [compat] -AbstractAlgebra = "0.43.12" AutoHashEquals = "2.2.0" -BenchmarkTools = "1.6.0" CairoMakie = "0.13.1" Colors = "0.13.0" Combinatorics = "1.0.2" @@ -58,11 +52,9 @@ GraphMakie = "0.5.13" GraphPlot = "0.6.0" Graphs = "1.12.0" Makie = "0.19.11" -Memoization = "0.2.2" NetworkLayout = "0.4.9" Oscar = "1.2.2" StatsBase = "0.34.4" -ThreadSafeDicts = "0.1.6" WGLMakie = "0.8.15" julia = "≥ 1.10" diff --git a/ext/MakieExt/LDPC/codes.jl b/ext/MakieExt/LDPC/codes.jl index 83ba8f1..cab1f92 100644 --- a/ext/MakieExt/LDPC/codes.jl +++ b/ext/MakieExt/LDPC/codes.jl @@ -43,131 +43,6 @@ function CodingTheory.degree_distributions_plot(C::AbstractLDPCCode) return f end -# """ -# $TYPEDSIGNATURES - -# Return a bar graph and a dictionary of (length, count) pairs for unique short -# cycles in the Tanner graph of `C`. An empty graph and dictionary are returned -# when there are no cycles. - -# # Note -# - Short cycles are defined to be those with lengths between ``g`` and ``2g - 2``, -# where ``g`` is the girth. -# - Run `using Makie` to activate this extension. -# """ -# function CodingTheory.count_short_cycles_plot(C::AbstractLDPCCode) -# if isempty(C.short_cycle_counts) || isempty(C.elementary_cycle_counts) -# CodingTheory._count_cycles(C) -# end - -# len = length(C.short_cycle_counts) -# x_data = [0 for _ in 1:len] -# y_data = [0 for _ in 1:len] -# index = 1 -# for (i, j) in C.short_cycle_counts -# x_data[index] = i -# y_data[index] = j -# index += 1 -# end - -# fig = Figure(); -# ax = Axis(fig[1, 1], xlabel = "Cycle Length", ylabel = "Occurrences", -# title = "Short Cycle Counts") -# barplot!(ax, x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data) -# # fig = Plots.bar(x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data, -# # legend = false, xlabel = "Cycle Length", ylabel = "Occurrences", -# # title = "Short Cycle Counts") -# display(fig) -# return fig, C.short_cycle_counts -# end - -# """ -# $TYPEDSIGNATURES - -# Return a bar graph and a dictionary of (length, count) pairs for unique elementary -# cycles in the Tanner graph of `C`. An empty graph and dictionary are returned -# when there are no cycles. - -# # Note -# - Elementary cycles do not contain the same vertex twice and are unable to be -# decomposed into a sequence of shorter cycles. -# - Run `using Makie` to activate this extension. -# """ -# function CodingTheory.count_elementary_cycles_plot(C::AbstractLDPCCode) -# if isempty(C.short_cycle_counts) || isempty(C.elementary_cycle_counts) -# CodingTheory._count_cycles(C) -# end - -# len = length(C.elementary_cycle_counts) -# x_data = [0 for _ in 1:len] -# y_data = [0 for _ in 1:len] -# index = 1 -# for (i, j) in C.elementary_cycle_counts -# x_data[index] = i -# y_data[index] = j -# index += 1 -# end - -# fig = Figure(); -# ax = Axis(fig[1, 1], xlabel = "Cycle Length", ylabel = "Occurrences", -# title = "Elementary Cycle Counts") -# barplot!(ax, x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data) -# # fig = Plots.bar(x_data, y_data, bar_width = 1, xticks = x_data, yticks = y_data, -# # legend = false, xlabel = "Cycle Length", ylabel = "Occurrences", -# # title = "Elementary Cycle Counts") -# display(fig) -# return fig, C.elementary_cycle_counts -# end - -""" -$TYPEDSIGNATURES - -Return an interactive figure and data for the ACE spectrum of the Tanner graph of `C`. - -# Note -- Run `using Makie` to activate this extension. -""" -function CodingTheory.ACE_spectrum_plot(C::AbstractLDPCCode) - # TODO: remove WGLMakie as a default use and only use for interactive plots - fig = Figure(); - ax = Axis(fig[1, 1], xlabel = "ACE", ylabel = "Occurrences", title = "ACE Spectrum") - sg = SliderGrid(fig[2, 1], (label = "Cycle Length", range = girth:2:2 * girth - 2, - startvalue = 4)) - - x_max = maximum(reduce(vcat, vs_ACEs)) - y_max = 0 - counts = ACE_spectrum(C) - X_data = Observable(Vector{Int}()) - Y_data = Observable(Vector{Int}()) - barplot!(ax, X_data, Y_data, bar_width = 1, xticks = X_data, yticks = Y_data) - for (k, l) in enumerate(girth:2:2 * girth - 2) - - ys = collect(values(counts[k])) - y_max = maximum([y_max; ys]) - - on(sg.sliders[1].value) do val - if to_value(val) == l - X_data.val = collect(keys(counts[k])) - Y_data.val = ys - notify(X_data) - notify(Y_data) - end - end - end - - GLMakie.limits!(0, x_max + 1, 0, y_max + 2) - display(fig) - return fig, counts -end - -mutable struct _ComputationGraphNode - id::Int - parent_id::Int - lvl::Int - vertex_number::Int - type::Symbol -end - # doesn't seem to be a point in making this dynamic with a slider, as it simply # continues in the same tree shape and no useful information is gained from watching it """ diff --git a/src/CodingTheory.jl b/src/CodingTheory.jl index 3120e09..746e530 100644 --- a/src/CodingTheory.jl +++ b/src/CodingTheory.jl @@ -20,8 +20,6 @@ using StatsBase using Distributions using ProgressMeter using DocStringExtensions -using Memoization -using ThreadSafeDicts import LinearAlgebra: tr, Adjoint, transpose, kron, diagm import Oscar: dual, factor, transpose, order, polynomial, nrows, ncols, degree, diff --git a/src/iterators.jl b/src/iterators.jl index 4c6d106..9c2705a 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -12,36 +12,6 @@ struct SubsetGrayCode init_rank::UInt128 # iteration will start with the subset of this rank end -function SubsetGrayCode(n::Int, k::Int) - init_rank = big(1) - num_threads = 1 - return _subset_gray_code_from_num_threads(n, k, init_rank, num_threads) -end - -function _subset_gray_code_from_num_threads(n::Int, k::Int, init_rank::UInt128, num_threads::Int) - bin = extended_binomial(n, k) - if bin % num_threads != 0 - throw(ArgumentError("num_threads=$num_threads, must divide binomial($n $k)")) - end - len = fld(bin, num_threads) - return SubsetGrayCode(n, k, len, init_rank) -end - -function _subset_gray_codes_from_num_threads(n::Int, k::Int, num_threads::Int) - # This function splits a single iterator into several pieces. - # The intended usage is to do the iteration with multiple threads - bin = extended_binomial(n, k) - i = 0 - len = fld(bin, num_threads) - - itrs = fill(SubsetGrayCode(1,1,1,1), num_threads) - for i in 0:num_threads - 1 - init_rank = 1 + i * len - itrs[i+1] = SubsetGrayCode(n, k, len, init_rank) - end - return itrs -end - Base.IteratorEltype(::SubsetGrayCode) = Base.HasEltype() Base.eltype(::SubsetGrayCode) = Array{Int, 1} Base.IteratorSize(::SubsetGrayCode) = Base.HasLength() @@ -76,21 +46,12 @@ end Note that inds is part of the iterator's state inds only to prevent reallocation in each iteration. =# - # subset_vec = zeros(Int, G.k) - # CodingTheory._subset_unrank!(G.init_rank, UInt(G.n), subset_vec) - - # a = @allocated begyyin subset_vec = zeros(Int, G.k) if G.init_rank == 1 subset_vec = Int.(collect(1 : G.k)) - # first_vec = vcat(ones(Int, r), zeros(Int, k - r)) - # for i in 1:G.k - # subset_vec[i] = 1 - # end else CodingTheory._subset_unrank!(G.init_rank, UInt(G.n), subset_vec) end - # end; a > 0 && @show a inds = fill(-1, 3) (inds, (subset_vec, 1, inds)) @@ -219,37 +180,20 @@ function _subset_unrank!(r::UInt128, n::UInt, T::Vector{Int}) k = length(T) subset_size_str::String = "subset size k = $k must be smaller than the set size n = $n" k > n && throw(ArgumentError(subset_size_str)) - # count_bins = 0 - # bnd = extended_binomial(n, k) - # # count_bins += 1 - # rank_size_str::String = "rank must be in [0, choose(n, k) - 1] = $bnd but is $(r)" - # # r > bnd && throw(ArgumentError(rank_size_str)) - # if r > bnd - # println(rank_size_str) - # println("bug exit _subset_unrank") - # return - # end - + x = 0 i = 0 y = 0 x = n for i::UInt in k:-1:1 y = extended_binomial(x, i) - # println((Int16.(x),Int16.(i))) - # count_bins += 1 while y > r x = x - 1 y = extended_binomial(x, i) - # println((x,i)) - # count_bins += 1 end T[i] = x + 1 r = extended_binomial(x + 1, i) - r - 1 - # println((x+1,i)) - # count_bins += 1 end - # println("num calls to binomial(a,b) in 1 unrank = $(count_bins)") end struct GrayCode diff --git a/src/utils.jl b/src/utils.jl index 7ce029a..97f8f85 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -164,7 +164,6 @@ function _min_wt_row(A::Union{CTMatrixTypes, Matrix{S}, LinearAlgebra.Adjoint{S, return w, i end -# function _min_wt_col(A::Union{CTMatrixTypes, Matrix{S}, LinearAlgebra.Adjoint{S, Matrix{S}}}) where S <: Integer function _min_wt_col(A::Union{CTMatrixTypes, Matrix{S}, LinearAlgebra.Adjoint{S, Matrix{S}}, LinearAlgebra.Adjoint{Bool, BitMatrix}}) where S <: Integer w = size(A, 1) + 1 i = 0 @@ -1993,7 +1992,6 @@ function _rand_invertible_matrix(F::CTFieldTypes, n::Integer) return A end -# @memoize ThreadSafeDict function extended_binomial(x::Union{Int, UInt}, y::Union{Int, UInt}) function extended_binomial(x::Union{Int, UInt}, y::Union{Int, UInt}) return y <= x ? UInt128.(binomial(x, y)) : UInt128(0) end From 445de92af47fe6779242a0360d82eec07ddf345a Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Mon, 24 Feb 2025 12:03:11 -0500 Subject: [PATCH 22/23] Remove stored found keywords --- src/Classical/weight_dist.jl | 69 ++++++++++++++---------------------- src/utils.jl | 1 + 2 files changed, 28 insertions(+), 42 deletions(-) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 3f4ed1f..52adedd 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -323,7 +323,6 @@ function information_sets(G::CTMatrixTypes, alg::Symbol = :Edmonds; permute::Boo if permute # permute identities to the front pivots = collect(start_ind:(i + 1) * nr) - @assert det(Gi[:, pivots]) != 0 σ = [pivots; setdiff(1:nc, pivots)] Gi = Gi[:, σ] Pi = Pi[:, σ] @@ -424,8 +423,12 @@ Return the minimum distance of `C` using a deterministic algorithm based on enum constant weight codewords of the binary reflected Gray code. If a word of minimum weight is found before the lower and upper bounds cross, it is returned; otherwise, the zero vector is returned. + +show_progress will display a progress meter for each iteration of a weight that takes longer than + a second """ -function minimum_distance_Gray(C::AbstractLinearCode; alg::Symbol = :auto, v::Bool = false, show_progress=true) +function minimum_distance_Gray(C::AbstractLinearCode; alg::Symbol = :auto, v::Bool = false, + show_progress = true) ord_F = Int(order(C.F)) ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes.")) @@ -451,9 +454,9 @@ function minimum_distance_Gray(C::AbstractLinearCode; alg::Symbol = :auto, v::Bo alg = :Zimmermann end end - alg == :auto || throw(ErrorException("Could not determine minimum distance algo automatically")) + alg == :auto && throw(ErrorException("Could not determine minimum distance algo automatically")) - if alg in (:Brouwer, :Zimmermann) || + if alg in (:Brouwer, :Zimmermann) return _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg = alg, verbose = v, show_progress = show_progress) end println("Warning: old enumeration algorithm selected. Performance will be slow") # TODO remove when all code updated @@ -462,6 +465,7 @@ end function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zimmermann, verbose::Bool = false, dbg = Dict(), show_progress=false) + dbg_key_exit_r = "exit_r" ord_F = Int(order(C.F)) ord_F == 2 || throw(ArgumentError("Currently only implemented for binary codes.")) @@ -488,16 +492,10 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim if length(keys(dbg)) > 0 println("Debug mode ON") end - dbg_key_found = "found_codewords" - dbg_key_exit_r = "exit_r" if haskey(dbg, dbg_key_exit_r) - verbose && println("dbg Dict: Storing the largest message weight reached @key=$dbg_key_found") - dbg["exit_r"] = -1 - end - store_found_codewords = haskey(dbg, dbg_key_found) - if store_found_codewords - verbose && println("dbg Dict: Storing the codewords in the debug dictionary @key=$dbg_key_found") + verbose && println("dbg Dict: largest message weight searched stored @key=$dbg_key_exit_r") + dbg[dbg_key_exit_r] = -1 end k, n = size(C.G) @@ -539,12 +537,6 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim # following loop is the r=1 case of the enumeration. We do this case here because we want to make a good guess at the terminating r before we start multiple threads for (j, g) in enumerate(A_mats) # loop over the A_mats rather than the original G because it would add another case to deal with later # can make this faster with dots and views - if store_found_codewords - for a in 1:size(g, 2) - output_vec = Int.(perms_mats[j] * g[:, a]) - push!(dbg["found_codewords"], (1, [], output_vec)) - end - end w, i = _min_wt_col(g) if w <= C.u_bound found = g[:, i] @@ -554,7 +546,6 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end verbose && println("Current upper bound: $(C.u_bound)") - verbose && !iszero(found) && println("Found element matching upper bound.") found = A_mats[1][:, 1] l = 0 @@ -595,7 +586,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim verbose && println("Number of threads ", num_thrds) for r in 2:k if r > 2^16 - verbose && println("Reached an r larger than 2^16") + verbose && println("Warning: Reached an r larger than 2^16") end C.l_bound < lower_bounds[r] && (C.l_bound = lower_bounds[r];) # an even code can't have have an odd minimum weight @@ -603,22 +594,19 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim # (!triply_even_flag && doubly_even_flag) && (C.l_bound += 4 - C.l_bound % 4;) # triply_even_flag && (C.l_bound += 8 - C.l_bound % 8;) if C.l_bound >= C.u_bound - C.d = C.u_bound - y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) - dbg[dbg_key_exit_r] = r - show_progress && ProgressMeter.finish!(prog_bar) - return C.u_bound, y + dbg[dbg_key_exit_r] = r-1 + break end verbose && println("r: $r") verbose && println("Lower bound: $(C.l_bound)") verbose && println("Upper bound: $(C.u_bound)") if verbose - count = 0 + i_count = 0 for i in 1:h - r - rank_defs[i] ≤ 0 && (count += 1;) + r - rank_defs[i] ≤ 0 && (i_count += 1;) end - count > 0 && println("$count of the original $h information sets no longer contribute to the lower bound") + i_count > 0 && println("$i_count of the original $h information sets no longer contribute to the lower bound") end p = Int(characteristic(C.F)) @@ -645,7 +633,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim CodingTheory._subset_unrank_to_vec!(init_rank, UInt64(r), first_vec) end - # as in White Algo 7.1 we loop over matrices first so that we can update the lower bound more often + # as in White Algo 7.1 we loop over matrices first for i in 1:h if keep_going[] == false verbose && println(thrd_stop_msg) @@ -655,6 +643,8 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim c_itr = zeros(UInt16, C.n - C.k) is_first = true curr_mat = A_mats_trunc[i] + count = UInt128(0) + for u in SubsetGrayCode(k, r, len, init_rank) if keep_going[] == false println(thrd_stop_msg) @@ -678,18 +668,11 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end end - if store_found_codewords - subset_vec_full = zeros(Int, k) - CodingTheory._subset_unrank_to_vec!(UInt128(init_rank + count), UInt64(r), subset_vec_full) - output_vec = Int.(perms_mats[exit_thread_indicator_vec[ind]] * vcat(subset_vec_full, c_itr)) - push!(dbg["found_codewords"], (r, [], Int.(output_vec))) - end - partial_weight = r + sum(view(c_itr, 1:weight_sum_bound)) if uppers[ind] > partial_weight w = r + sum(c_itr) - verbose && @assert w == 0 + verbose && @assert w != 0 if uppers[ind] > w subset_vec_full = zeros(Int, k) CodingTheory._subset_unrank_to_vec!(UInt128(init_rank + count), UInt64(r), subset_vec_full) @@ -699,7 +682,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim verbose && @assert size(founds[ind], 1) == C.n "found vector has length $(size(founds[ind], 1)) but should be n=$(C.n)" exit_thread_indicator_vec[ind] = i - verbose && println("Adjusting (local) upper bound: $w for c_itr=$(c_itr)") + println("Adjusting (local) upper bound: $w for c_itr=$(Int.(c_itr))") if C.l_bound == uppers[ind] println("early exit") Threads.atomic_cas!(keep_going, true, false) @@ -719,12 +702,15 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim initial_perm_ind = exit_thread_indicator_vec[loc] end end - show_progress && ProgressMeter.finish!(prog_bar) - # at this point we are guaranteed to have found the code's distance C.d = C.u_bound - y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) # weight(y) >= C.d with equality not being the typical case + y = matrix(C.F, 1, n, perms_mats[initial_perm_ind] * found) # weight(y) >= C.d, with equality not being the typical case verbose && @assert iszero(C.H * transpose(y)) + if dbg[dbg_key_exit_r] == -1 + dbg[dbg_key_exit_r] = r + end + show_progress && ProgressMeter.finish!(prog_bar) + verbose && println("Computation complete") return C.u_bound, y end @@ -827,7 +813,6 @@ function _minimum_distance_enumeration_with_matrix_multiply(C::AbstractLinearCod if C.l_bound >= C.u_bound C.d = C.u_bound y = perms_mats[initial_perm_ind] * found - # @assert in(y, C) return C.u_bound, (C.F).(y) end diff --git a/src/utils.jl b/src/utils.jl index 97f8f85..2b58ff1 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -403,6 +403,7 @@ function _has_empty_vec(A::Union{CTMatrixTypes, Matrix{<: Number}, BitMatrix, Ma end return !isempty(del) end + function _remove_empty(A::Union{CTMatrixTypes, Matrix{<: Number}, BitMatrix, Matrix{Bool}}, type::Symbol) From aef41b2de3a85d1219ebcdf8cf13a2a3cd5a0e37 Mon Sep 17 00:00:00 2001 From: "david.marquis" Date: Mon, 24 Feb 2025 16:08:20 -0500 Subject: [PATCH 23/23] Clean --- src/Classical/weight_dist.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 52adedd..5fc6093 100644 --- a/src/Classical/weight_dist.jl +++ b/src/Classical/weight_dist.jl @@ -694,6 +694,7 @@ function _minimum_distance_BZ(C::AbstractLinearCode; info_set_alg::Symbol = :Zim end end end + count = add!(count, count, 1) end end loc = argmin(uppers)