Slight optimizations to BPGates

src/BPGates.jl

@@ -437,6 +437,8 @@ const good_perm_qc = ( # From the appendix of Optimized Entanglement Purificatio
     (h*ph,p*hp)
 )
 
+# Does this only modify the second bell pair? No phase kickback?-- interesting.
+# phase kickback here contributes to a global phase
 const cnot_perm = (1, 2, 11, 12, 6, 5, 16, 15, 9, 10, 3, 4, 14, 13, 8, 7)
 
 function QuantumClifford.apply!(state::BellState, g::CNOTPerm) # TODO abstract away the permutation application as it is used by other gates too
@@ -571,7 +573,7 @@ end
 # TODO: Honestly this doesn't make that much sense to me. Think about why
 # this works.
 # Defines the probabilities of transitioning from one bell state to another
-function get_mixed_transition_probs(λ::Float64, cur_state_idx::Int)
+function get_mixed_transition_probs(λ, cur_state_idx)
     # 0 <= λ <= 1 (λ = 1 - e ^(-t / T1))
     mixed_state_tuple = (
         (0.5 * λ^2 - λ + 1, 0.5 * λ * (1 - λ), 0.5 * λ^2, 0.5 * λ * (1 - λ)),
@@ -591,7 +593,7 @@ struct MixedStateOp <: AbstractNoiseBellOp
     lambda::Float64
 end
 
-function select_rand_state(transition_probs::NTuple{4, Float64})
+function select_rand_state(transition_probs)
     # TODO: hardcoded for speed, BUT change this to a macro for readability.
     rand_prob_val = rand()
     new_phase_idx = 0

test/runtests.jl

@@ -16,19 +16,26 @@ end
 macro doset(descr)
     quote
         if doset($descr)
-            @safetestset $descr begin include("test_"*$descr*".jl") end
+            @safetestset $descr begin 
+            if $descr == "mixedstateop" || $descr == "paulizop"
+                include("test_kraus_mem_noise_helpers.jl")
+            end
+            include("test_"*$descr*".jl") end
         end
     end
 end
 
 println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREADS = $(Sys.CPU_THREADS)`...")
 
+# include("test_kraus_mem_noise_helpers.jl")
+
 @doset "quantumclifford"
 @doset "quantikz"
 get(ENV,"JET_TEST","")=="true" && @doset "jet"
 @doset "doctests"
 @doset "bpgates"
 @doset "mixedstateop"
+@doset "paulizop"
 
 using Aqua
 doset("aqua") && begin

test/test_kraus_mem_noise_helpers.jl

@@ -0,0 +1,259 @@
+module TestKrausNoiseHelpers
+
+using Test
+using Logging
+
+using BPGates: get_mixed_transition_probs, select_rand_state
+
+export test_bell_basis_probs, test_nonoise_bell_basis, test_almost_all_noise,
+test_random_lambdas, test_equal_mixed_probs_2, test_equal_mixed_probs_4,
+test_random_mixed_states, test_random_mixed_states_random_lambdas
+
+# using QuantumClifford: apply!
+
+# Get the bell states used for testing
+function get_bell_states()
+    return (1 / sqrt(2)) * [
+        [1.0; 0.0; 0.0; 1.0],
+        [0.0; 1.0; 1.0; 0.0],
+        [1.0; 0.0; 0.0; -1.0],
+        [0.0; 1.0; -1.0; 0.0]
+    ]
+end
+
+# Returns an index selected from the indices of a list, weighted by the probability at that index
+function get_weighted_index(prob_dist::Vector{Float64})
+    if !isapprox(sum(prob_dist), 1.0, atol=1e-9)
+        @error "get_weighted_index: prob_dist sums to $(sum(prob_dist)) and not 1"
+    end
+    rand_num = rand()
+    cur_sum = 0.0
+    for (idx, prob_val) in enumerate(prob_dist)
+        if rand_num < cur_sum + prob_val
+            return idx
+        end
+        cur_sum += prob_val
+    end
+end
+
+# Converts a bit to an int from left to right (LSB is the rightmost bit)
+function bit_to_int_testing(bits)
+    reduce(⊻,(bit<<(length(bits) - index) for (index,bit) in enumerate(bits))) + 1 # +1 so that we use julia indexing convenctions
+end
+
+# Returns the tensor product of a set of Kraus operators
+function get_tensored_kraus_ops(kraus_ops::Vector{Matrix{Complex}})
+    tensored_kraus_ops = []
+    for op1 in kraus_ops
+        for op2 in kraus_ops
+            push!(tensored_kraus_ops, kron(op1, op2))
+        end
+    end
+    return tensored_kraus_ops
+end
+
+# Returns the diagonal of a matrix
+function get_diagonal(matrix::Matrix{ComplexF64})::Vector{Float64}
+    diagonal_length = min(size(matrix)[1], size(matrix)[2])
+    if size(matrix)[1] != size(matrix)[2]
+        @error("get_diagonal: matrix is not square, rows: $(size(matrix)[1]) != cols: $(size(matrix)[2])")
+    end
+    diag_elements = []
+    for idx in 1:diagonal_length
+        push!(diag_elements, matrix[idx, idx])
+    end
+    return diag_elements
+end
+
+# Returns the bell state corresponding to the input phases
+function get_bell_state(phases::BitVector)
+    state_idx = bit_to_int_testing(phases)
+    bell_states = get_bell_states()
+    return bell_states[state_idx]
+end
+
+# Numerically computes the transition probabilities of one state to another,
+# given a mixed probability distribution of states and a set of kraus operators
+function get_numerical_trans_probs(mixed_state_probs::Vector{Float64}, kraus_ops::Vector{Matrix{Complex}})
+    bell_states = get_bell_states()
+    ρ_bell = zeros(ComplexF64, length(bell_states), length(bell_states))
+    for (bell_idx, bell_prob) in enumerate(mixed_state_probs)
+        ρ_bell += bell_prob * bell_states[bell_idx] * bell_states[bell_idx]'
+    end
+    @debug("get_numerical_trans_probs: ρ_bell: $ρ_bell")
+    tensored_ops = get_tensored_kraus_ops(kraus_ops)
+    # Apply the Kraus operators to the Bell state
+    ρ_new = zeros(ComplexF64, size(ρ_bell)[1], size(ρ_bell)[2])
+    for K in tensored_ops
+        ρ_new += K * ρ_bell * adjoint(K)
+    end
+    U_bell = hcat(get_bell_states()...)
+    @debug("get_numerical_trans_probs: ρ_new: $ρ_new")
+    @debug("get_numerical_trans_probs: U_bell: $U_bell")
+    ρ_Bell_basis = adjoint(U_bell) * ρ_new * U_bell
+    @debug("get_numerical_trans_probs: ρ_Bell_basis: $ρ_Bell_basis")
+    @debug("get_numerical_trans_probs: size(ρ_Bell_basis): $(size(ρ_Bell_basis))")
+    return get_diagonal(ρ_Bell_basis)
+end
+
+# Testing functions
+
+# Test our simulated noise for the pure states in the bell basis
+function test_bell_basis_probs(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    statedist_list = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
+    lambda = 1 - exp(-500e-9 / 260e-6)
+    num_samples = 1000000
+    kraus_ops = kraus_ops_factory(lambda)
+    for state_dist in statedist_list
+        @info "test_bell_basis_probs: state_dist: $state_dist"
+        numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+        simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+        # @info("test_A_state_probs: numerical_trans_probs: $numerical_trans_probs")
+        # @info("test_A_state_probs: simulated_trans_probs: $simulated_trans_probs")
+        @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+    end
+end
+
+# Test our simulated noise when there should be no noise applied for states in the bell basis
+function test_nonoise_bell_basis(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    statedist_list = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
+    lambda = 0.0
+    num_samples = 1000000
+    kraus_ops = kraus_ops_factory(lambda)
+    for state_dist in statedist_list
+        @info "test_nonoise_bell_basis: state_dist: $state_dist"
+        numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+        simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+        # @info("test_A_state_probs: numerical_trans_probs: $numerical_trans_probs")
+        # @info("test_A_state_probs: simulated_trans_probs: $simulated_trans_probs")
+        @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+    end
+end
+
+# Test our simulated noise with a noise factor close to 1 for states in the bell basis
+function test_almost_all_noise(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    statedist_list = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
+    lambda = 0.999999999
+    num_samples = 1000000
+    kraus_ops = kraus_ops_factory(lambda)
+    for state_dist in statedist_list
+        @info "test_almost_all_noise: state_dist: $state_dist"
+        numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+        simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+        @info("test_almost_all_noise: numerical_trans_probs: $numerical_trans_probs")
+        @info("test_almost_all_noise: simulated_trans_probs: $simulated_trans_probs")
+        @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+    end
+end
+
+# Test random noise factors for states in the bell basis
+function test_random_lambdas(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    statedist_list = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
+    num_samples = 1000000
+    num_random_lambdas = 5
+    for _ in 1:num_random_lambdas
+        lambda = rand()
+        num_samples = 1000000
+        kraus_ops = kraus_ops_factory(lambda)
+        for state_dist in statedist_list
+            @info "test_random_lambdas: state_dist: $state_dist"
+            @info "test_random_lambdas: lambda: $lambda"
+            numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+            simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+            @info("test_random_lambdas: numerical_trans_probs: $numerical_trans_probs")
+            @info("test_random_lambdas: simulated_trans_probs: $simulated_trans_probs")
+            @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+        end
+    end
+end
+
+# Test all mixed states composed of two states in the Bell basis with equal probability
+function test_equal_mixed_probs_2(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    bell_states = get_bell_states()
+    lambda = 1 - exp(-500e-9 / 260e-6)
+    kraus_ops = kraus_ops_factory(lambda)
+    num_samples = 1000000
+    for first_phases_idx in 1:length(bell_states)
+        for second_phases_idx in 1:length(bell_states)
+            if first_phases_idx == second_phases_idx
+                continue
+            end
+            state_dist = [0.0 for _ in 1:length(bell_states)]
+            for basis_idx in 1:length(bell_states)
+                if basis_idx == first_phases_idx || basis_idx == second_phases_idx
+                    state_dist[basis_idx] = 0.5
+                end
+            end
+            numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+            simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+            @info("test_equal_mixed_probs_2: numerical_trans_probs: $numerical_trans_probs")
+            @info("test_equal_mixed_probs_2: simulated_trans_probs: $simulated_trans_probs")
+            @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+        end
+    end
+end
+
+# Test a mixed state composed of all 4 bell basis states with equal probability
+function test_equal_mixed_probs_4(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    state_dist = [0.25 for _ in 1:length(get_bell_states())]
+    lambda = 1 - exp(-500e-9 / 260e-6)
+    num_samples = 1000000
+    kraus_ops = kraus_ops_factory(lambda)
+    @info "test_equal_mixed_probs_4: state_dist: $state_dist"
+    numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+    simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+    @info("test_equal_mixed_probs_4: numerical_trans_probs: $numerical_trans_probs")
+    @info("test_equal_mixed_probs_4: simulated_trans_probs: $simulated_trans_probs")
+    @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+end
+
+# Test random mixed states in the bell basis
+function test_random_mixed_states(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    num_random_states = 20
+    statedist_list = []
+    for _ in 1:num_random_states
+        state_dist = [rand() for _ in 1:length(get_bell_states())]
+        state_dist = state_dist / sum(state_dist)
+        push!(statedist_list, state_dist)
+    end
+    @debug "length(statedist_list): $(length(statedist_list))"
+    lambda = 1 - exp(-500e-9 / 260e-6)
+    num_samples = 1000000
+    kraus_ops = kraus_ops_factory(lambda)
+    for state_dist in statedist_list
+        @info "test_random_mixed_states: state_dist: $state_dist"
+        numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+        simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+        @info("test_random_mixed_states: numerical_trans_probs: $numerical_trans_probs")
+        @info("test_random_mixed_states: simulated_trans_probs: $simulated_trans_probs")
+        @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+    end
+end
+
+# Test random mixed states with random noise parameters in the Bell basis
+function test_random_mixed_states_random_lambdas(kraus_ops_factory::Function, get_simulated_trans_probs::Function)
+    num_random_states = 20
+    statedist_list = []
+    for _ in 1:num_random_states
+        state_dist = [rand() for _ in 1:length(get_bell_states())]
+        state_dist = state_dist / sum(state_dist)
+        push!(statedist_list, state_dist)
+    end
+    num_random_lambdas = 20
+    num_samples = 1000000
+    for _ in 1:num_random_lambdas
+        lambda = rand()
+        @info "test_random_mixed_states_random_lambdas: lambda: $lambda"
+        kraus_ops = kraus_ops_factory(lambda)
+        for state_dist in statedist_list
+            @info "test_random_mixed_states_random_lambdas: state_dist: $state_dist"
+            numerical_trans_probs = get_numerical_trans_probs(state_dist, kraus_ops)
+            simulated_trans_probs = get_simulated_trans_probs(state_dist, lambda, num_samples)
+            @info("test_random_mixed_states_random_lambdas: numerical_trans_probs: $numerical_trans_probs")
+            @info("test_random_mixed_states_random_lambdas: simulated_trans_probs: $simulated_trans_probs")
+            @test all(isapprox.(numerical_trans_probs, simulated_trans_probs, atol=(2 / sqrt(num_samples))))
+        end
+    end
+end
+
+end