Dev/decoherence

CHANGELOG

@@ -0,0 +1,3 @@
+- Created an abstract noise operation class `AbstractNoiseBellOp`
+- Implemented `PauliZOp` for phase damping errors and `MixedStateOp` for amplitude damping errors
+- Added tests for `PauliZOp` and `MixedStateOp` for comparison to exact density matrix computation

Project.toml

@@ -17,4 +17,4 @@ BPGatesQuantikzExt = "Quantikz"
 Quantikz = "1.3.1"
 QuantumClifford = "0.8, 0.9"
 Random = "1"
-julia = "1.9"
+julia = "1.9"

src/BPGates.jl

@@ -9,7 +9,8 @@
     BellSinglePermutation, BellDoublePermutation, BellPauliPermutation,
     BellMeasure, bellmeasure!,
     BellGate, CNOTPerm, GoodSingleQubitPerm,
-    PauliNoiseOp, PauliNoiseBellGate, NoisyBellMeasure, NoisyBellMeasureNoisyReset
+    PauliNoiseOp, PauliZOp, PauliNoiseBellGate, NoisyBellMeasure, NoisyBellMeasureNoisyReset, AbstractNoiseBellOp,
+    MixedStateOp, AsymmDepOp
 
 function int_to_bit(int,digits)
     int = int - 1 # -1 so that we use julia indexing conventions
@@ -168,6 +169,7 @@
     sidx::Int
     function BellMeasure(p,s)
         1 <= p <= 3 || throw(ArgumentError("The basis measurement index needs to be between 1 and 3"))
+        # why is there no check here for the measured bell pair? ensuring that the index is >= 1?
         new(p,s)
     end
 end
@@ -214,9 +216,13 @@
 
 """The permutations realized by [`BellPauliPermutation`](@ref)."""
 const pauli_perm_tuple = (
+    # 3 is actually 2; 2 is actually 3
     (1, 2, 3, 4),
+    # Z gate?
     (3, 4, 1, 2), ## X flip
+    # X gate?
     (2, 1, 4, 3), ## Z flip
+    # bit and phase flip
     (4, 3, 2, 1)  ## Y flip
 )
 
@@ -226,9 +232,16 @@
 function QuantumClifford.apply!(state::BellState, op::BellPauliPermutation) # TODO abstract away the permutation application as it is used by other gates too
     phase = state.phases
     @inbounds phase_idx = bit_to_int(phase[op.sidx*2-1],phase[op.sidx*2])
+    # println("QuantumClifford.apply!, BellPauliPermutation, op.sidx: $(op.sidx)")
+    # println("QuantumClifford.apply!, BellPauliPermutation, phase_idx: $phase_idx")
+    # println("QuantumClifford.apply!, BellPauliPermutation, op.pidx: $(op.pidx)")
     @inbounds perm = pauli_perm_tuple[op.pidx]
     @inbounds permuted_idx = perm[phase_idx]
+    # println("QuantumClifford.apply!, BellPauliPermutation, permuted_idx: $permuted_idx")
+    # this makes no sense. Make bit 2 go before bit1.
     bit1, bit2 = int_to_bit(permuted_idx,Val(2))
+    # println("QuantumClifford.apply!, BellPauliPermutation, bit1: $bit1")
+    # println("QuantumClifford.apply!, BellPauliPermutation, bit2: $bit2")
     @inbounds phase[op.sidx*2-1] = bit1
     @inbounds phase[op.sidx*2] = bit2
     return state
@@ -321,6 +334,7 @@
 function bellmeasure!(state::BellState, op::BellMeasure) # TODO document which index corresponds to which measurement
     phase = state.phases
     result = measure_tuple[bit_to_int(phase[op.sidx*2-1],phase[op.sidx*2])][op.midx]
+    # TODO: introduce latency here?
     phase[op.sidx*2-1:op.sidx*2] .= 0 # reset the measured pair to 00
     return state, result
 end
@@ -418,6 +432,8 @@
     (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
@@ -471,7 +487,10 @@
 
 """A wrapper for `BellGate` that implements Pauli noise in addition to the gate."""
 struct PauliNoiseBellGate{G} <: BellOp where {G<:BellOp} # TODO make it work with the QuantumClifford noise ops
+    # to do the above TODO, probably have to make it match some interface?
     g::G
+    # probability that the state changes in these bases after the noise is applied
+    # probability that the qubit at the specified index changes to an X gate.
     px::Float64
     py::Float64
     pz::Float64
@@ -518,6 +537,112 @@
     return state
 end
 
+# Type for abstracting away noisy operations on bell states.
+abstract type AbstractNoiseBellOp <: BellOp end
+
+"""
+PauliZOp(idx, pz) causes qubit at idx `idx` to have a Z flip with probabilty `pz`.
+"""
+struct PauliZOp <: AbstractNoiseBellOp
+    idx::Int
+    # assert that 0 <= pz <= 1? not strictly necessary. it wouldn't make sense, but it
+    # also wouldn't break the code. natural floor and ceiling here.
+    pz::Float64
+end
+
+function QuantumClifford.apply!(state::BellState, g::PauliZOp)
+    i = g.idx
+    r = rand()
+    if r<g.pz
+        # I believe this actually is a Z flip.
+        apply!(state, BellPauliPermutation(2, i))
+    end
+    return state
+end
+
+# TODO: continue from here
+# TODO: assert that the values add to 1
+# This is NOT at all 'extensible' (though, it is 'fast')
+# 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(λ, 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 - λ)),
+        (0.5 * λ, 1 - λ, 0.5 * λ, 0.0),
+        (0.5 * λ^2, 0.5 * λ * (1 - λ), 0.5 * λ^2 - λ + 1, 0.5 * λ * (1 - λ)),
+        (0.5 * λ, 0.0, 0.5 * λ, 1 - λ)
+    )
+    return mixed_state_tuple[cur_state_idx]
+end
+
+"""
+MixedStateOp(idx, lambda) causes Bell state at idx `idx` to change to another
+Bell state determined by `mixed_state_tuple` and `lambda`.
+"""
+struct MixedStateOp <: AbstractNoiseBellOp
+    idx::Int
+    lambda::Float64
+end
+
+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
+    if rand_prob_val < transition_probs[1]
+        new_phase_idx = 1
+    elseif rand_prob_val < transition_probs[1] + transition_probs[2]
+        new_phase_idx = 2
+    elseif rand_prob_val < transition_probs[1] + transition_probs[2] + transition_probs[3]
+        new_phase_idx = 3
+    else
+        new_phase_idx = 4
+    end
+    return new_phase_idx
+end
+
+function QuantumClifford.apply!(state::BellState, g::MixedStateOp)
+    state_idx = g.idx
+    phase = state.phases
+    @inbounds phase_idx = bit_to_int(phase[state_idx * 2 - 1],phase[state_idx * 2])
+    transition_probs = get_mixed_transition_probs(g.lambda, phase_idx)
+    new_phase_idx = select_rand_state(transition_probs)
+    # this should never happen; can remove for speed
+    @assert new_phase_idx != 0
+    state_bit1, state_bit2 = int_to_bit(new_phase_idx, Val(2))
+    @inbounds phase[state_idx * 2 - 1] = state_bit1
+    @inbounds phase[state_idx * 2] = state_bit2
+    return state
+end
+
+struct AsymmDepOp <: AbstractNoiseBellOp
+    idx::Int
+    px::Float64
+    py::Float64
+    pz::Float64
+end
+
+function QuantumClifford.apply!(state::BellState, g::AsymmDepOp)
+    i = g.idx
+    # TODO repetition with ...NoisyReset and PauliNoise...
+    # ^ ?
+    r = rand()
+    # 2: Z, 3: X, 4: Y
+    # println("QuantumClifford.apply!, AsymmDepOp")
+    # print("QuantumClifford.apply!, AsymmDepOp, px_bell: $px_bell")
+    # print("QuantumClifford.apply!, AsymmDepOp, py_bell: $py_bell")
+    # print("QuantumClifford.apply!, AsymmDepOp, pz_bell: $pz_bell")
+    if r<g.px
+        apply!(state, BellPauliPermutation(3, i))
+    elseif r<g.px+g.pz
+        apply!(state, BellPauliPermutation(2, i))
+    elseif r<g.px+g.pz+g.py
+        apply!(state, BellPauliPermutation(4, i))
+    end
+    return state
+end
+
 """A wrapper for [`BellMeasure`](@ref) that implements measurement noise."""
 struct NoisyBellMeasure <: BellOp # TODO make it work with the QuantumClifford noise ops
     m::BellMeasure
@@ -527,7 +652,9 @@
 """A wrapper for [`BellMeasure`](@ref) that implements measurement noise and Pauli noise after the reset."""
 struct NoisyBellMeasureNoisyReset <: BellOp # TODO make it work with the QuantumClifford noise ops
     m::BellMeasure
+    # this p is the probability that the incorrect measurement is returned?
     p::Float64
+    # probability of resetting to a different Bell state after measurement
     px::Float64
     py::Float64
     pz::Float64
@@ -538,8 +665,10 @@
     state, result⊻(rand()<op.p) ? continue_stat : failure_stat
 end
 
+# This function is like a 'continuation' of the above function
 function QuantumClifford.applywstatus!(state::BellState, op::NoisyBellMeasureNoisyReset)
     state, result = bellmeasure!(state, op.m)
+    # pretty high probability of flipping to a correct measurement?
     cont = result⊻(rand()<op.p)
     cont && apply!(state, PauliNoiseOp(op.m.sidx,op.px,op.py,op.pz))
     state, cont ? continue_stat : failure_stat

test/Project.toml

@@ -6,6 +6,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Quantikz = "b0d11df0-eea3-4d79-b4a5-421488cbf74b"
 QuantumClifford = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

test/runtests.jl

@@ -16,17 +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_bpgates.jl

@@ -0,0 +1,108 @@
+using Test
+using Logging
+
+using BPGates
+
+using QuantumClifford
+
+function test_pauliz_constructor()
+
+    test_pauliz = PauliZOp(1, 0.5)
+
+    @test test_pauliz isa PauliZOp
+
+    @test test_pauliz.idx == 1
+
+    @test test_pauliz.pz == 0.5
+end
+
+function test_pauliz_application_guaranteed()
+    n_bellpairs = 3
+    changed_bp_idx = 2
+    test_state = BellState(n_bellpairs)
+    test_guaranteed_pauliz = PauliZOp(changed_bp_idx, 1)
+
+    QuantumClifford.apply!(test_state, test_guaranteed_pauliz)
+
+    @test test_state.phases[2*changed_bp_idx - 1] == 0
+    @test test_state.phases[2*changed_bp_idx] == 1
+end
+
+function test_pauliz_application_guaranteed_none()
+    n_bellpairs = 3
+    changed_bp_idx = 2
+    test_state = BellState(n_bellpairs)
+    test_guaranteed_pauliz = PauliZOp(changed_bp_idx, 0)
+
+    # TODO: do I have to import QuantumClifford to use it when testing?
+    QuantumClifford.apply!(test_state, test_guaranteed_pauliz)
+
+    @test test_state.phases[2*changed_bp_idx - 1] == 0
+    @test test_state.phases[2*changed_bp_idx] == 0
+end
+
+function test_mixed_state_op_constructor()
+    mixed_state_op = MixedStateOp(1, 0.5)
+
+    @test mixed_state_op isa MixedStateOp
+
+    @test mixed_state_op.idx == 1
+
+    @test mixed_state_op.lambda == 0.5
+end
+
+function test_apply_mixed_state_op()
+    n_bellpairs = 1
+    changed_bp_idx = 1
+    test_state = BellState(n_bellpairs)
+    mixed_state_op = MixedStateOp(changed_bp_idx, 0.0)
+
+    QuantumClifford.apply!(test_state, mixed_state_op)
+
+    @test test_state.phases[2 * changed_bp_idx - 1] == 0
+    @test test_state.phases[2 * changed_bp_idx] == 0
+end
+
+function test_apply_mixed_state_op_diff_bellstate()
+    n_bellpairs = 1
+    changed_bp_idx = 1
+    test_state = BellState((0, 1))
+    mixed_state_op = MixedStateOp(changed_bp_idx, 0.0)
+
+    QuantumClifford.apply!(test_state, mixed_state_op)
+
+    @test test_state.phases[2 * changed_bp_idx - 1] == 0
+    @test test_state.phases[2 * changed_bp_idx] == 1
+end
+
+function test_apply_both_memory_errors()
+    n_bellpairs = 1
+    changed_bp_idx = 1
+
+    test_state = BellState(n_bellpairs)
+
+    noise_ops = Vector{AbstractNoiseBellOp}()
+
+    push!(noise_ops, PauliZOp(changed_bp_idx, 0))
+    push!(noise_ops, MixedStateOp(changed_bp_idx, 0.0))
+
+    for noise_op in noise_ops
+        QuantumClifford.apply!(test_state, noise_op)
+    end
+
+    @test test_state.phases[2 * changed_bp_idx - 1] == 0
+    @test test_state.phases[2 * changed_bp_idx] == 0
+end
+
+@testset "BPGates.jl, PauliZOp tests" begin
+    test_pauliz_constructor()
+    test_pauliz_application_guaranteed()
+    test_pauliz_application_guaranteed_none()
+end
+
+@testset "BPGates.jl, MixedStateOp tests" begin
+    test_mixed_state_op_constructor()
+    test_apply_mixed_state_op()
+    test_apply_mixed_state_op_diff_bellstate()
+    test_apply_both_memory_errors()
+end