Dev/decoherence

CHANGELOG

@@ -0,0 +1,2 @@
+- Created an abstract noise operation class `AbstractNoiseBellOp`
+- Implemented `PauliZOp` for phase damping errors and `MixedStateOp` for amplitude damping errors

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
 
 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
@@ -215,8 +217,11 @@
 """The permutations realized by [`BellPauliPermutation`](@ref)."""
 const pauli_perm_tuple = (
     (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
 )
 
@@ -260,6 +265,7 @@
 
 """The permutations realized by [`BellDoublePermutation`](@ref) as Clifford operations."""
 const double_perm_qc = ( # TODO switch to symbolic gates
+    # Q: are these the stabilizer states for each qubit?
     C"X_ _Z Z_ _X",
     C"_X X_ _Z Z_",
     C"XX _X Z_ ZZ",
@@ -319,8 +325,11 @@
 ```
 """
 function bellmeasure!(state::BellState, op::BellMeasure) # TODO document which index corresponds to which measurement
+    # does this apply both coincidence AND anti-coincidence measurements?
+    # i assume that it's hard coded into measure_tuple
     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
@@ -379,6 +388,7 @@
     idx2::Int
     function CNOTPerm(s1,s2,i1,i2)
         (1 <= s1 <= 6 && 1 <= s2 <= 6) || throw(ArgumentError("The permutation index needs to be between 1 and 6."))
+        # no other checks on the BP indices? maybe that's the job of the other component
         (i1 > 0 && i2 > 0) || throw(ArgumentError("The Bell pair indices have to be positive integers."))
         i1 != i2 || throw(ArgumentError("The gate has to act on two different Bell pairs, i.e. idx1!=idx2."))
         new(s1,s2,i1,i2)
@@ -409,6 +419,7 @@
 const p = tPhase
 const hp = h*p
 const ph = p*h
+# TODO: What is this????
 const good_perm_qc = ( # From the appendix of Optimized Entanglement Purification, but be careful with index notation being different
     (tId1,tId1), # TODO switch to symbolic gates
     (h*ph*ph,h*hp*hp*hp*hp),
@@ -422,6 +433,7 @@
 
 function QuantumClifford.apply!(state::BellState, g::CNOTPerm) # TODO abstract away the permutation application as it is used by other gates too
     phase = state.phases
+    # why can you do inbounds here?
     @inbounds phase_idxa = bit_to_int(phase[g.idx1*2-1],phase[g.idx1*2])
     @inbounds phase_idxb = bit_to_int(phase[g.idx2*2-1],phase[g.idx2*2])
     if phase_idxa==phase_idxb==1
@@ -471,7 +483,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
@@ -507,6 +522,7 @@
 function QuantumClifford.apply!(state::BellState, g::PauliNoiseOp)
     i = g.idx
     # TODO repetition with ...NoisyReset and PauliNoise...
+    # ^ ?
     r = rand()
     if r<g.px
         apply!(state, BellPauliPermutation(2, i))
@@ -518,6 +534,80 @@
     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(λ::Float64, cur_state_idx::Int)
+    # 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.5 * λ^2, 0.5 * λ * (1 - λ), 0.5 * λ^2 - λ + 1, 0.5 * λ * (1 - λ)),
+        (0.5 * λ, 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 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])
+    rand_prob_val = rand()
+    transition_probs = get_mixed_transition_probs(g.lambda, phase_idx)
+    # TODO: hardcoded for speed, BUT change this to a macro for readability.
+    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
+    # 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
+
 """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 +617,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 +630,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

@@ -27,6 +27,7 @@ println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREA
 @doset "quantikz"
 get(ENV,"JET_TEST","")=="true" && @doset "jet"
 @doset "doctests"
+@doset "bpgates"
 
 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