Merge branch 'dev' into fixdev

Author: Fe-r-oz

src/Quantum/weight_dist.jl

@@ -657,124 +657,49 @@ function random_information_set_minimum_distance_bound!(::HasGauges, ::IsCSS, ::
     # this is a CSS subsystem code
 
     n = S.n
-    if dressed
-        if which == :full
-            verbose && println("Bounding the full dressed distance")
-            stabs = _Flint_matrix_to_Julia_T_matrix(stabilizers(S), Int)
-            _rref_no_col_swap_binary!(stabs)
-            stabs = _remove_empty(stabs, :rows)
-            gauges = _Flint_matrix_to_Julia_T_matrix(gauges_matrix(S), Int)
-            logs = _Flint_matrix_to_Julia_T_matrix(logicals_matrix(S), Int)
-            operators_to_reduce = vcat(stabs, gauges, logs)
-            check_against = permutedims(logs[:, [n + 1:2n; 1:n]])
-            curr_l_bound = S.l_bound_dressed
-            verbose && println("Starting lower bound: $curr_l_bound")
-
-            # this is done in the constructor but the logical is not stored at the time
-            # so must redo here
-            mat = _rref_no_col_swap_binary(operators_to_reduce)
-            anti = mat * check_against
-            curr_u_bound, index = findmin(row_wts_symplectic(mat[findall(!iszero(anti[i, :]) for i in axes(anti, 1)), :]))
-            found = operators_to_reduce[index, :]
-            verbose && println("Starting upper bound: $curr_u_bound")
-        elseif which == :X
-            verbose && println("Bounding the dressed X-distance")
-            stabs = _Flint_matrix_to_Julia_T_matrix(stabilizers(S), Int)
-            _rref_no_col_swap_binary!(stabs)
-            stabs = _remove_empty(stabs, :rows)
-            gauges = _Flint_matrix_to_Julia_T_matrix(reduce(vcat, [log[1][:, 1:n] for log in S.gauge_ops]), Int)
-            logs = _Flint_matrix_to_Julia_T_matrix(reduce(vcat, [log[1][:, 1:n] for log in S.logicals]), Int)
-            operators_to_reduce = vcat(stabs, gauges, logs)
-            check_against = permutedims(reduce(vcat, [log[2][:, n + 1:end] for log in S.logicals])[:, [n + 1:2n; 1:n]])
-            curr_l_bound = S.l_bound_dx_dressed
-            verbose && println("Starting lower bound: $curr_l_bound")
-
-            # this is done in the constructor but the logical is not stored at the time
-            # so must redo here
-            mat = _rref_no_col_swap_binary(operators_to_reduce)
-            anti = mat * check_against
-            curr_u_bound, index = findmin(row_wts_symplectic(mat[findall(!iszero(anti[i, :]) for i in axes(anti, 1)), :]))
-            found = operators_to_reduce[index, :]
-            verbose && println("Starting upper bound: $curr_u_bound")
+    if which == :full
+        upperx, foundx = random_information_set_minimum_distance_bound!(HasGauges(), IsCSS(), HasLogicals(), S, :X, dressed, max_iters, verbose)
+        upperz, foundz = random_information_set_minimum_distance_bound!(HasGauges(), IsCSS(), HasLogicals(), S, :Z, dressed, max_iters, verbose)
+        if upperx <= upperz
+            if dressed
+                S.u_bound_dressed = min(upperx, S.u_bound_dressed)
+            else
+                S.u_bound_bare = min(upperx, S.u_bound_bare)
+            end
+            return upperx, foundx
         else
-            verbose && println("Bounding the dressed Z-distance")
-            stabs = _Flint_matrix_to_Julia_T_matrix(stabilizers(S), Int)
-            _rref_no_col_swap_binary!(stabs)
-            stabs = _remove_empty(stabs, :rows)
-            gauges = _Flint_matrix_to_Julia_T_matrix(reduce(vcat, [log[2][:, n + 1:end] for log in S.gauge_ops]), Int)
-            logs = _Flint_matrix_to_Julia_T_matrix(reduce(vcat, [log[2][:, n + 1:end] for log in S.logicals]), Int)
-            operators_to_reduce = vcat(stabs, gauges, logs)
-            check_against = permutedims(reduce(vcat, [log[1][:, 1:n] for log in S.logicals]))
-            curr_l_bound = S.l_bound_dz_dressed
-            verbose && println("Starting lower bound: $curr_l_bound")
-
-            # this is done in the constructor but the logical is not stored at the time
-            # so must redo here
-            mat = _rref_no_col_swap_binary(operators_to_reduce)
-            anti = mat * check_against
-            curr_u_bound, index = findmin(row_wts_symplectic(mat[findall(!iszero(anti[i, :]) for i in axes(anti, 1)), :]))
-            found = operators_to_reduce[index, :]
-            verbose && println("Starting upper bound: $curr_u_bound")
+            if dressed
+                S.u_bound_dressed = min(upperz, S.u_bound_dressed)
+            else
+                S.u_bound_bare = min(upperz, S.u_bound_bare)
+            end
+            return upperz, foundz
         end
     else
-        if which == :full
-            verbose && println("Bounding the full bare distance")
-            stabs = _Flint_matrix_to_Julia_T_matrix(stabilizers(S), Int)
-            _rref_no_col_swap_binary!(stabs)
-            stabs = _remove_empty(stabs, :rows)
-            logs = _Flint_matrix_to_Julia_T_matrix(logicals_matrix(S), Int)
-            operators_to_reduce = vcat(stabs, logs)
-            check_against = permutedims(logs[:, [n + 1:2n; 1:n]])
-            curr_l_bound = S.l_bound_bare
-            verbose && println("Starting lower bound: $curr_l_bound")
-
-            # this is done in the constructor but the logical is not stored at the time
-            # so must redo here
-            mat = _rref_no_col_swap_binary(operators_to_reduce)
-            anti = mat * check_against
-            curr_u_bound, index = findmin(row_wts_symplectic(mat[findall(!iszero(anti[i, :]) for i in axes(anti, 1)), :]))
-            found = operators_to_reduce[index, :]
-            verbose && println("Starting upper bound: $curr_u_bound")
-        elseif which == :X
-            verbose && println("Bounding the bare X-distance")
-            stabs = _Flint_matrix_to_Julia_T_matrix(stabilizers(S), Int)
-            _rref_no_col_swap_binary!(stabs)
-            stabs = _remove_empty(stabs, :rows)
-            logs = _Flint_matrix_to_Julia_T_matrix(reduce(vcat, [log[1][:, 1:n] for log in S.logicals]), Int)
-            operators_to_reduce = vcat(stabs, logs)
-            check_against = permutedims(reduce(vcat, [log[2][:, n + 1:end] for log in S.logicals])[:, [n + 1:2n; 1:n]])
-            curr_l_bound = S.l_bound_dx_bare
-            verbose && println("Starting lower bound: $curr_l_bound")
-
-            # this is done in the constructor but the logical is not stored at the time
-            # so must redo here
-            mat = _rref_no_col_swap_binary(operators_to_reduce)
-            anti = mat * check_against
-            curr_u_bound, index = findmin(row_wts_symplectic(mat[findall(!iszero(anti[i, :]) for i in axes(anti, 1)), :]))
-            found = operators_to_reduce[index, :]
-            verbose && println("Starting upper bound: $curr_u_bound")
+        stabs = _Flint_matrix_to_Julia_T_matrix(stabilizers(S)[:, (which == :X ? (1:n) : (n + 1:2n))], UInt8)
+        if dressed
+            gauges = _Flint_matrix_to_Julia_T_matrix(gauge_operators_matrix(S)[:, (which == :X ? (1:n) : (n + 1:2n))], UInt8)
+            stabs = [stabs; gauges]
+        end
+        _rref_no_col_swap_binary!(stabs)
+        stabs = _remove_empty(stabs, :rows)
+        logs = _Flint_matrix_to_Julia_T_matrix(logicals_matrix(S)[:, (which == :X ? (1:n) : (n + 1:2n))], UInt8)
+        logs = _remove_empty(logs, :rows)
+        operators_to_reduce = vcat(stabs, logs)
+        check_against = _Flint_matrix_to_Julia_T_matrix(logicals_matrix(S)[:, (which == :X ? (n + 1:2n) : (1:n))], UInt8)
+        check_against = permutedims(_remove_empty(check_against, :rows))
+        curr_l_bound = if dressed
+            which == :X ? S.l_bound_dx_dressed : S.l_bound_dz_dressed
         else
-            verbose && println("Bounding the bare Z-distance")
-            stabs = _Flint_matrix_to_Julia_T_matrix(stabilizers(S), Int)
-            _rref_no_col_swap_binary!(stabs)
-            stabs = _remove_empty(stabs, :rows)
-            logs = _Flint_matrix_to_Julia_T_matrix(reduce(vcat, [log[2][:, n + 1:end] for log in S.logicals]), Int)
-            operators_to_reduce = vcat(stabs, logs)
-            check_against = permutedims(reduce(vcat, [log[1][:, 1:n] for log in S.logicals]))
-            curr_l_bound = S.l_bound_dx_bare
-            verbose && println("Starting lower bound: $curr_l_bound")
-
-            # this is done in the constructor but the logical is not stored at the time
-            # so must redo here
-            mat = _rref_no_col_swap_binary(operators_to_reduce)
-            anti = mat * check_against
-            curr_u_bound, index = findmin(row_wts_symplectic(mat[findall(!iszero(anti[i, :]) for i in axes(anti, 1)), :]))
-            found = operators_to_reduce[index, :]
-            verbose && println("Starting upper bound: $curr_u_bound")
+            which == :X ? S.l_bound_dx_bare : S.l_bound_dz_bare
         end
+        verbose && println("Starting lower bound: $curr_l_bound")
+        curr_u_bound, index = findmin(count(!iszero, logs[i, :]) for i in 1:size(logs, 1))
+        found = logs[index, :]
+        verbose && println("Starting upper bound: $curr_u_bound")
     end
 
-    uppers, founds = _RIS_bound_loop_symp!(operators_to_reduce, check_against, curr_l_bound,
+    uppers, founds = _RIS_bound_loop!(operators_to_reduce, check_against, curr_l_bound,
         curr_u_bound, found, max_iters, n, verbose)
     loc = argmin(uppers)
     if dressed
@@ -795,7 +720,14 @@ function random_information_set_minimum_distance_bound!(::HasGauges, ::IsCSS, ::
         end
     end
     verbose && println("Ending $max_iters iterations with an upper bound of $(uppers[loc])")
-    return uppers[loc], founds[loc]
+    flint_mat_found = if which == :full
+        matrix(field(S), permutedims(founds[loc]))
+    elseif which == :X
+        matrix(field(S), [permutedims(founds[loc]) zeros(Int, 1, n)])
+    else
+        matrix(field(S), [zeros(Int, 1, n) permutedims(founds[loc])])
+    end
+    return uppers[loc], flint_mat_found
 end
 
 function random_information_set_minimum_distance_bound!(::HasNoGauges, ::IsCSS, ::HasLogicals,
@@ -979,7 +911,8 @@ function _RIS_bound_loop_symp!(operators_to_reduce, check_against, curr_l_bound:
                 shuffle!(perm)
                 _col_permutation_symp!(perm_ops, orig_ops, perm)
                 # modifying this in place is not thread safe (apparently)
-                perm_ops = _rref_no_col_swap_binary(perm_ops)
+                # perm_ops = _rref_no_col_swap_binary(perm_ops)
+                _rref_no_col_swap_binary!(perm_ops)
                 _col_permutation_symp!(ops, perm_ops, invperm(perm))
                 LinearAlgebra.mul!(log_test, ops, check_against)
                 for i in axes(log_test, 1)

src/utils.jl

@@ -425,31 +425,31 @@ function _remove_empty(A::Union{CTMatrixTypes, Matrix{<: Number}, BitMatrix, Mat
     end
 end
 
-function _rref_no_col_swap(M::CTMatrixTypes, row_range::UnitRange{Int}, col_range::UnitRange{Int})
+function _rref_no_col_swap(M::CTMatrixTypes, row_range::AbstractUnitRange{Int}, col_range::AbstractUnitRange{Int})
     A = deepcopy(M)
     _rref_no_col_swap!(A, row_range, col_range)
     return A
 end
-_rref_no_col_swap(M::CTMatrixTypes, row_range::Base.OneTo{Int}, col_range::Base.OneTo{Int}) = _rref_no_col_swap(M, 1:row_range.stop, 1:col_range.stop)
+_rref_no_col_swap(M::CTMatrixTypes, row_range::AbstractUnitRange{Int}, col_range::AbstractUnitRange{Int}) = _rref_no_col_swap(M, 1:row_range.stop, 1:col_range.stop)
 _rref_no_col_swap(M::CTMatrixTypes) = _rref_no_col_swap(M, axes(M, 1), axes(M, 2))
 
 function _rref_no_col_swap_binary(A::Union{BitMatrix, Matrix{Bool}, Matrix{<: Integer}},
-    row_range::UnitRange{Int} = 1:size(A, 1), col_range::UnitRange{Int} = 1:size(A, 2))
+    row_range::AbstractUnitRange{Int} = 1:size(A, 1), col_range::AbstractUnitRange{Int} = 1:size(A, 2))
 
     B = deepcopy(A)
     _rref_no_col_swap_binary!(B, row_range, col_range)
     return B
 end
 
-function _rref_no_col_swap!(A::CTMatrixTypes, row_range::UnitRange{Int}, col_range::UnitRange{Int})
+function _rref_no_col_swap!(A::CTMatrixTypes, row_range::AbstractUnitRange{Int}, col_range::AbstractUnitRange{Int})
     # don't do anything to A if the range is empty
     isempty(row_range) && return nothing
     isempty(col_range) && return nothing
 
-    i = row_range.start
-    j = col_range.start
-    nr = row_range.stop
-    nc = col_range.stop
+    i = first(row_range)
+    j = first(col_range)
+    nr = last(row_range)
+    nc = last(col_range)
     if Int(order(base_ring(A))) != 2
         while i <= nr && j <= nc
             # find first pivot
@@ -523,14 +523,14 @@ function _rref_no_col_swap!(A::CTMatrixTypes, row_range::UnitRange{Int}, col_ran
 end
 
 function _rref_no_col_swap_binary!(A::Union{BitMatrix, Matrix{Bool}, Matrix{<: Integer}},
-    row_range::UnitRange{Int} = 1:size(A, 1), col_range::UnitRange{Int} = 1:size(A, 2))
+    row_range::AbstractUnitRange{Int} = 1:size(A, 1), col_range::AbstractUnitRange{Int} = 1:size(A, 2))
 
     isempty(row_range) && return nothing
     isempty(col_range) && return nothing
-    i = row_range.start
-    j = col_range.start
-    nr = row_range.stop
-    nc = col_range.stop
+    i = first(row_range)
+    j = first(col_range)
+    nr = last(row_range)
+    nc = last(col_range)
     while i <= nr && j <= nc
         # find first pivot
         ind = 0
@@ -565,15 +565,15 @@ function _rref_no_col_swap_binary!(A::Union{BitMatrix, Matrix{Bool}, Matrix{<: I
     return nothing
 end
 
-function _rref_col_swap(M::CTMatrixTypes, row_range::UnitRange{Int}, col_range::UnitRange{Int})
+function _rref_col_swap(M::CTMatrixTypes, row_range::AbstractUnitRange{Int}, col_range::AbstractUnitRange{Int})
     A = deepcopy(M)
     rnk, P = _rref_col_swap!(A, row_range, col_range)
     return rnk, A, P
 end
-_rref_col_swap(M::CTMatrixTypes, row_range::Base.OneTo{Int}, col_range::Base.OneTo{Int}) = _rref_col_swap(M, 1:row_range.stop, 1:col_range.stop)
+_rref_col_swap(M::CTMatrixTypes, row_range::AbstractUnitRange{Int}, col_range::AbstractUnitRange{Int}) = _rref_col_swap(M, 1:row_range.stop, 1:col_range.stop)
 _rref_col_swap(M::CTMatrixTypes) = _rref_col_swap(M, axes(M, 1), axes(M, 2))
 
-function _rref_col_swap!(A::CTMatrixTypes, row_range::UnitRange{Int}, col_range::UnitRange{Int})
+function _rref_col_swap!(A::CTMatrixTypes, row_range::AbstractUnitRange{Int}, col_range::AbstractUnitRange{Int})
     # don't do anything to A if the range is empty, return rank 0 and missing permutation matrix
     isempty(row_range) && return 0, missing
     isempty(col_range) && return 0, missing
@@ -583,10 +583,10 @@ function _rref_col_swap!(A::CTMatrixTypes, row_range::UnitRange{Int}, col_range:
     nc_A = ncols(A)
 
     rnk = 0
-    i = row_range.start
-    j = col_range.start
-    nr = row_range.stop
-    nc = col_range.stop
+    i = first(row_range)
+    j = first(col_range)
+    nr = last(row_range)
+    nc = last(col_range)
     if Int(order(base_ring(A))) != 2
         while i <= nr && j <= nc
             # find first pivot
@@ -626,7 +626,7 @@ function _rref_col_swap!(A::CTMatrixTypes, row_range::UnitRange{Int}, col_range:
                 ind != i && swap_rows!(A, ind, i)
 
                 # eliminate
-                for k = row_range.start:nr
+                for k = first(row_range):nr
                     if k != i
                         # do a manual loop here to reduce allocations
                         d = A[k, j]
@@ -674,7 +674,7 @@ function _rref_col_swap!(A::CTMatrixTypes, row_range::UnitRange{Int}, col_range:
                 ind != i && swap_rows!(A, ind, i)
 
                 # eliminate
-                for k = row_range.start:nr
+                for k = first(row_range):nr
                     if k != i
                         if isone(A[k, j])
                             # do a manual loop here to reduce allocations
@@ -693,7 +693,7 @@ function _rref_col_swap!(A::CTMatrixTypes, row_range::UnitRange{Int}, col_range:
     return rnk, P
 end
 
-function _rref_symp_col_swap!(A::CTMatrixTypes, row_range::UnitRange{Int}, col_range::UnitRange{Int})
+function _rref_symp_col_swap!(A::CTMatrixTypes, row_range::AbstractUnitRange{Int}, col_range::AbstractUnitRange{Int})
     # don't do anything to A if the range is empty, return rank 0 and missing permutation matrix
     isempty(row_range) && return 0, missing
     isempty(col_range) && return 0, missing
@@ -703,10 +703,10 @@ function _rref_symp_col_swap!(A::CTMatrixTypes, row_range::UnitRange{Int}, col_r
     nc_A = ncols(A)
 
     rnk = 0
-    i = row_range.start
-    j = col_range.start
-    nr = row_range.stop
-    nc = col_range.stop
+    i = first(row_range)
+    j = first(col_range)
+    nr = last(row_range)
+    nc = last(col_range)
     if Int(order(base_ring(A))) != 2
         while i <= nr && j <= nc
             # find first pivot
@@ -750,7 +750,7 @@ function _rref_symp_col_swap!(A::CTMatrixTypes, row_range::UnitRange{Int}, col_r
                 ind != i && swap_rows!(A, ind, i)
 
                 # eliminate
-                for k = row_range.start:nr
+                for k = first(row_range):nr
                     if k != i
                         # do a manual loop here to reduce allocations
                         d = A[k, j]
@@ -802,7 +802,7 @@ function _rref_symp_col_swap!(A::CTMatrixTypes, row_range::UnitRange{Int}, col_r
                 ind != i && swap_rows!(A, ind, i)
 
                 # eliminate
-                for k = row_range.start:nr
+                for k = first(row_range):nr
                     if k != i
                         if isone(A[k, j])
                             # do a manual loop here to reduce allocations
@@ -824,10 +824,10 @@ _rref_symp_col_swap!(A::CTMatrixTypes) = _rref_symp_col_swap!(A, axes(A, 1), axe
 _rref_symp_col_swap(A::CTMatrixTypes) = (B = deepcopy(A); return _rref_symp_col_swap!(B);)
 
 function _col_permutation!(X::Matrix{T}, A::Matrix{T}, p::AbstractVector{Int}) where T
-    length(p) == size(A, 2) || error()
-    size(X) == size(A) || error()
-    for i in axes(X, 1)
-        for j in axes(X, 2)
+    length(p) == size(A, 2) || throw(ArgumentError("`p` should have length `size(A, 2)`."))
+    size(X) == size(A) || throw(ArgumentError("`X` and `A` should have the same shape."))
+    for j in axes(X, 2)
+        for i in axes(X, 1)
             X[i, j] = A[i, p[j]]
         end
     end
@@ -836,13 +836,18 @@ end
 
 function _col_permutation_symp!(X::Matrix{T}, A::Matrix{T}, p::AbstractVector{Int}) where T
     n = length(p)
-    2n == size(A, 2) || error()
-    size(X) == size(A) || error()
-    for i in axes(X, 1)
-        for j in axes(X, 2)
+    2n == size(A, 2) || throw(ArgumentError("`p` should have length `size(A, 2)/2`."))
+    size(X) == size(A) || throw(ArgumentError("`X` and `A` should have the same shape."))
+    for j in 1:n
+        for i in axes(X, 1)
             X[i, j] = A[i, p[mod1(j, n)]]
         end
     end
+    for j in 1:n
+        for i in axes(X, 1)
+            X[i, j + n] = A[i, p[n] + n]
+        end
+    end
     return nothing
 end