|
| 1 | +#BPOTS Tests |
| 2 | + |
| 3 | +using Test |
| 4 | +using SparseArrays |
| 5 | +using Random |
| 6 | +using QuantumClifford.ECC |
| 7 | +import LDPCDecoders: BPOTSDecoder, BeliefPropagationDecoder, decode!, reset! |
| 8 | + |
| 9 | + |
| 10 | +#= |
| 11 | +Tests syndrome matching |
| 12 | +create_cycle_matrix(n)- creates a cycle graph with known trapping set properties |
| 13 | +generate_random_syndrome(H)- generates random error patterns (computes their syndromes using parity check matrices) |
| 14 | +Remove the excess debug statements and make them more concise |
| 15 | +=# |
| 16 | +@testset "Comprehensive BP-OTS Decoder Tests" begin |
| 17 | + # Utility functions |
| 18 | + function create_cycle_matrix(n::Int) |
| 19 | + I_idx = Int[] |
| 20 | + J_idx = Int[] |
| 21 | + for j in 1:n |
| 22 | + push!(I_idx, j) |
| 23 | + push!(J_idx, mod(j,n)+1) |
| 24 | + push!(I_idx, j) |
| 25 | + push!(J_idx, j) |
| 26 | + end |
| 27 | + V = fill(true, length(I_idx)) |
| 28 | + return sparse(I_idx, J_idx, V, n, n) |
| 29 | + end |
| 30 | + |
| 31 | + function generate_random_syndrome(H::SparseMatrixCSC{Bool}) |
| 32 | + m, n = size(H) |
| 33 | + error = rand(Bool, n) |
| 34 | + return Vector{Bool}(mod.(H * error, 2)) |
| 35 | + end |
| 36 | + |
| 37 | + function test_decoder_performance(decoder_type, H, noise_levels, max_iterations) |
| 38 | + results = [] |
| 39 | + for noise in noise_levels |
| 40 | + success_count = 0 |
| 41 | + total_runs = 100 |
| 42 | + |
| 43 | + for _ in 1:total_runs |
| 44 | + syndrome = generate_random_syndrome(H) |
| 45 | + decoder = decoder_type(H, noise, max_iterations) |
| 46 | + reset!(decoder) |
| 47 | + result, converged = decode!(decoder, syndrome) |
| 48 | + decoded_syndrome = Bool.(mod.(H * result, 2)) |
| 49 | + |
| 50 | + if decoded_syndrome == syndrome |
| 51 | + success_count += 1 |
| 52 | + end |
| 53 | + end |
| 54 | + |
| 55 | + push!(results, (noise=noise, success_rate=success_count/total_runs)) |
| 56 | + end |
| 57 | + return results |
| 58 | + end |
| 59 | + |
| 60 | + @testset "Trapping Set Resistance" begin |
| 61 | + # Test with cycle matrices of different sizes |
| 62 | + for n in [4, 8, 16] |
| 63 | + H = create_cycle_matrix(n) |
| 64 | + |
| 65 | + # Test BP-OTS with different configurations |
| 66 | + parameter_sets = [ |
| 67 | + (T=3, C=1.0), |
| 68 | + (T=5, C=2.0), |
| 69 | + (T=9, C=3.0) |
| 70 | + ] |
| 71 | + |
| 72 | + for (T, C) in parameter_sets |
| 73 | + bpots_decoder = BPOTSDecoder(H, 0.01, 100; T=T, C=C) |
| 74 | + reset!(bpots_decoder) |
| 75 | + |
| 76 | + # Creates a weight-2 error pattern (known to form trapping sets) |
| 77 | + error = zeros(Bool, n) |
| 78 | + error[1:2] .= true |
| 79 | + syndrome = Vector{Bool}(mod.(H * error, 2)) |
| 80 | + |
| 81 | + result, converged = decode!(bpots_decoder, syndrome) |
| 82 | + decoded_syndrome = Bool.(mod.(H * result, 2)) |
| 83 | + |
| 84 | + @test decoded_syndrome == syndrome |
| 85 | + println("Cycle Matrix n=$n, T=$T, C=$C: Decoded successfully") |
| 86 | + end |
| 87 | + end |
| 88 | + end |
| 89 | + |
| 90 | + @testset "Syndrome Matching Robustness" begin |
| 91 | + # Test BPOTS on different types of matrices |
| 92 | + test_matrices = [ |
| 93 | + sparse([1,2,3,4], [1,2,3,4], true, 4, 4), # Diagonal matrix |
| 94 | + create_cycle_matrix(4), # Cycle matrix |
| 95 | + sparse([1,1,2,2,3,3], [1,2,2,3,3,4], true, 4, 4) # Irregular matrix |
| 96 | + ] |
| 97 | + |
| 98 | + noise_levels = [0.01, 0.1, 0.5] |
| 99 | + |
| 100 | + for H in test_matrices |
| 101 | + # Test BP-OTS performance |
| 102 | + performance_results = test_decoder_performance( |
| 103 | + (H, noise, max_iter) -> BPOTSDecoder(H, noise, max_iter; T=9, C=3.0), |
| 104 | + H, noise_levels, 100 |
| 105 | + ) |
| 106 | + |
| 107 | + println("Performance for matrix:") |
| 108 | + for result in performance_results |
| 109 | + println(" Noise: $(result.noise), Success Rate: $(result.success_rate)") |
| 110 | + @test result.success_rate > 0.5 # Expect more than 50% success |
| 111 | + end |
| 112 | + end |
| 113 | + end |
| 114 | + |
| 115 | + @testset "Parameter Sensitivity" begin |
| 116 | + H = create_cycle_matrix(4) |
| 117 | + |
| 118 | + # Test T (threshold) values |
| 119 | + for T in [3, 5, 9, 15] |
| 120 | + bpots = BPOTSDecoder(H, 0.01, 100; T=T, C=3.0) |
| 121 | + reset!(bpots) |
| 122 | + |
| 123 | + syndrome = generate_random_syndrome(H) |
| 124 | + result, converged = decode!(bpots, syndrome) |
| 125 | + decoded_syndrome = Bool.(mod.(H * result, 2)) |
| 126 | + |
| 127 | + println("T=$T - Converged: $converged, Syndrome match: $(decoded_syndrome == syndrome)") |
| 128 | + @test decoded_syndrome == syndrome |
| 129 | + end |
| 130 | + |
| 131 | + # Test C values |
| 132 | + for C in [1.0, 2.0, 5.0, 10.0] |
| 133 | + bpots = BPOTSDecoder(H, 0.01, 100; T=9, C=C) |
| 134 | + reset!(bpots) |
| 135 | + |
| 136 | + syndrome = generate_random_syndrome(H) |
| 137 | + result, converged = decode!(bpots, syndrome) |
| 138 | + decoded_syndrome = Bool.(mod.(H * result, 2)) |
| 139 | + |
| 140 | + println("C=$C - Converged: $converged, Syndrome match: $(decoded_syndrome == syndrome)") |
| 141 | + @test decoded_syndrome == syndrome |
| 142 | + end |
| 143 | + end |
| 144 | + |
| 145 | + """ |
| 146 | + generate_toric_code_matrix(; d::Int=3) |
| 147 | + Generate a parity check matrix for a toric code with distance `d`. |
| 148 | + Just use toric tht exists |
| 149 | + parity_checks_x(Toric(n,n)) |
| 150 | + parity_checks_z(Toric(n,n)) |
| 151 | + Can implement surface as well (edit(Toric)) |
| 152 | + """ |
| 153 | + function generate_toric_code_matrix(; d::Int=3) |
| 154 | + # Number of physical qubits is 2 * d^2 |
| 155 | + n_qubits = 2 * d^2 |
| 156 | + |
| 157 | + # Number of stabilizers (checks) is d^2 for X-type and d^2 for Z-type |
| 158 | + n_checks = 2 * d^2 |
| 159 | + |
| 160 | + # Prepare for sparse matrix construction |
| 161 | + I_idx = Int[] # Row indices |
| 162 | + J_idx = Int[] # Column indices |
| 163 | + V = Bool[] # Values |
| 164 | + |
| 165 | + # Create X-type stabilizers (plaquettes) |
| 166 | + for i in 0:(d-1) |
| 167 | + for j in 0:(d-1) |
| 168 | + check_idx = i*d + j + 1 |
| 169 | + |
| 170 | + # Four edges of the plaquette |
| 171 | + qubit_h1 = i*d + j + 1 |
| 172 | + qubit_h2 = i*d + ((j+1) % d) + 1 |
| 173 | + qubit_v1 = d^2 + i*d + j + 1 |
| 174 | + qubit_v2 = d^2 + ((i+1) % d)*d + j + 1 |
| 175 | + |
| 176 | + # Add entries to sparse matrix |
| 177 | + for qubit in [qubit_h1, qubit_h2, qubit_v1, qubit_v2] |
| 178 | + push!(I_idx, check_idx) |
| 179 | + push!(J_idx, qubit) |
| 180 | + push!(V, true) |
| 181 | + end |
| 182 | + end |
| 183 | + end |
| 184 | + |
| 185 | + # Create Z-type stabilizers (stars) |
| 186 | + for i in 0:(d-1) |
| 187 | + for j in 0:(d-1) |
| 188 | + check_idx = d^2 + i*d + j + 1 |
| 189 | + |
| 190 | + # Four edges meeting at a vertex |
| 191 | + qubit_h1 = i*d + j + 1 |
| 192 | + qubit_h2 = i*d + ((j-1+d) % d) + 1 |
| 193 | + qubit_v1 = d^2 + i*d + j + 1 |
| 194 | + qubit_v2 = d^2 + ((i-1+d) % d)*d + j + 1 |
| 195 | + |
| 196 | + # Add entries to sparse matrix |
| 197 | + for qubit in [qubit_h1, qubit_h2, qubit_v1, qubit_v2] |
| 198 | + push!(I_idx, check_idx) |
| 199 | + push!(J_idx, qubit) |
| 200 | + push!(V, true) |
| 201 | + end |
| 202 | + end |
| 203 | + end |
| 204 | + |
| 205 | + return sparse(I_idx, J_idx, V, n_checks, n_qubits) |
| 206 | + end |
| 207 | + |
| 208 | + """ |
| 209 | + generate_surface_code_matrix(; d::Int=3) |
| 210 | + Generate a parity check matrix for a surface code with distance `d`. |
| 211 | + """ |
| 212 | + function generate_surface_code_matrix(; d::Int=3) |
| 213 | + if d < 2 |
| 214 | + error("Surface code distance must be at least 2") |
| 215 | + end |
| 216 | + |
| 217 | + # For a simple surface code, use a lattice with open boundaries |
| 218 | + # Number of qubits: (d-1)*d horizontal + d*(d-1) vertical |
| 219 | + n_qubits = 2 * d * (d-1) |
| 220 | + |
| 221 | + # Number of checks: (d-1)² plaquettes + d² vertices - 1 (dependent check) |
| 222 | + n_checks = (d-1)^2 + d^2 - 1 |
| 223 | + |
| 224 | + # Prepare for sparse matrix construction |
| 225 | + I_idx = Int[] # Row indices |
| 226 | + J_idx = Int[] # Column indices |
| 227 | + V = Bool[] # Values |
| 228 | + |
| 229 | + # Horizontal edge at (row, col) |
| 230 | + h_edge(row, col) = row*(d-1) + col + 1 |
| 231 | + |
| 232 | + # Vertical edge at (row, col) |
| 233 | + v_edge(row, col) = (d-1)*d + row*d + col + 1 |
| 234 | + |
| 235 | + # Create plaquette operators (X-type stabilizers) |
| 236 | + check_idx = 1 |
| 237 | + for row in 0:(d-2) |
| 238 | + for col in 0:(d-2) |
| 239 | + edges = [ |
| 240 | + h_edge(row, col), # Top |
| 241 | + h_edge(row+1, col), # Bottom |
| 242 | + v_edge(row, col), # Left |
| 243 | + v_edge(row, col+1) # Right |
| 244 | + ] |
| 245 | + |
| 246 | + for edge in edges |
| 247 | + push!(I_idx, check_idx) |
| 248 | + push!(J_idx, edge) |
| 249 | + push!(V, true) |
| 250 | + end |
| 251 | + check_idx += 1 |
| 252 | + end |
| 253 | + end |
| 254 | + |
| 255 | + # Create vertex operators (Z-type stabilizers) |
| 256 | + # Skip the last vertex as it's dependent on others |
| 257 | + for row in 0:d-1 |
| 258 | + for col in 0:d-1 |
| 259 | + # Skip the last vertex |
| 260 | + if row == d-1 && col == d-1 |
| 261 | + continue |
| 262 | + end |
| 263 | + |
| 264 | + # Collect edges connected to this vertex |
| 265 | + edges = Int[] |
| 266 | + |
| 267 | + # Horizontal edges |
| 268 | + if col > 0 |
| 269 | + push!(edges, h_edge(row, col-1)) # Left |
| 270 | + end |
| 271 | + if col < d-1 |
| 272 | + push!(edges, h_edge(row, col)) # Right |
| 273 | + end |
| 274 | + |
| 275 | + # Vertical edges |
| 276 | + if row > 0 |
| 277 | + push!(edges, v_edge(row-1, col)) # Top |
| 278 | + end |
| 279 | + if row < d-1 |
| 280 | + push!(edges, v_edge(row, col)) # Bottom |
| 281 | + end |
| 282 | + |
| 283 | + # Add entries to sparse matrix |
| 284 | + for edge in edges |
| 285 | + push!(I_idx, check_idx) |
| 286 | + push!(J_idx, edge) |
| 287 | + push!(V, true) |
| 288 | + end |
| 289 | + check_idx += 1 |
| 290 | + end |
| 291 | + end |
| 292 | + |
| 293 | + return sparse(I_idx, J_idx, V, n_checks, n_qubits) |
| 294 | + end |
| 295 | + |
| 296 | + @testset "Toric and Surface Code Performance" begin |
| 297 | + # Generate matrices with a smaller distance for testing |
| 298 | + d_test = 3 |
| 299 | + println("Generating toric code matrix (d=$d_test)...") |
| 300 | + #H_toric = generate_toric_code_matrix(d=d_test) |
| 301 | + H_toric = parity_checks_x(Toric(d_test, d_test)) |
| 302 | + |
| 303 | + println("Generating surface code matrix (d=$d_test)...") |
| 304 | + H_surface = generate_surface_code_matrix(d=d_test) |
| 305 | + |
| 306 | + # Test with reduced noise levels and iterations for faster testing |
| 307 | + noise_levels = [0.01, 0.05, 0.1] |
| 308 | + max_iter = 50 |
| 309 | + |
| 310 | + for (matrix_idx, H) in enumerate([H_toric, H_surface]) |
| 311 | + matrix_name = matrix_idx == 1 ? "Toric Code" : "Surface Code" |
| 312 | + println("\nTesting $matrix_name ($(size(H)[1])×$(size(H)[2]))...") |
| 313 | + |
| 314 | + performance_results = test_decoder_performance( |
| 315 | + (H, noise, max_iter) -> BPOTSDecoder(H, noise, max_iter; T=9, C=3.0), |
| 316 | + H, noise_levels, max_iter |
| 317 | + ) |
| 318 | + |
| 319 | + for result in performance_results |
| 320 | + println(" $matrix_name @ Noise $(result.noise): Success Rate $(result.success_rate)") |
| 321 | + @test result.success_rate > 0.5 |
| 322 | + end |
| 323 | + end |
| 324 | + end |
| 325 | + |
| 326 | +end |
0 commit comments