Added 3to1 example for measuring fidelity

examples/threetoonepurification/README.md

@@ -0,0 +1,25 @@
+# Single and Double Selection Purification
+Purification can be done using one or two pairs of sacrificial qubits. If one uses two pairs instead of one, the output fidelity increases ( especially for higher input fidelities ). 
+
+To measure their performance we start by creating a noisy pair with the given input fidelity `fid`, and then we measure the success probability and final fidelity for a sample of `simcount = 10000` times.
+
+```
+noisy_pair = noisy_pair_func(fid)
+successcount = 0
+finfid = 1
+noisy_pair = noisy_pair_func(fid)
+for _ in 1:simcount
+    r = Register(6, QuantumOpticsRepr())
+    initialize!(r[1:2], noisy_pair)
+    initialize!(r[3:4], noisy_pair)
+    initialize!(r[5:6], noisy_pair)
+    output = Purify3to1(leaveout1, leaveout2)(r[1], r[2], r[3], r[5], r[4], r[6])
+    (output) && (successcount = successcount + 1)
+    if output
+        finfid = observable(r[1:2], projector(bell))
+    end
+end
+```
+
+The success rate is then given by the expression `successcount/simcount`. This are then plotted to yield the image in this same folder.
+

examples/threetoonepurification/purificationPlots.jl

@@ -0,0 +1,103 @@
+using GLMakie
+GLMakie.activate!()
+
+f = Figure(resolution=(1200,900))
+using QuantumSavory
+using QuantumSavory.CircuitZoo: EntanglementSwap, Purify2to1Node, Purify3to1, Purify2to1, Purify3to1Node, AbstractCircuit
+# legend and axis names
+const bell = StabilizerState("XX ZZ")
+# QOptics repr
+const perfect_pair = (Z1⊗Z1 + Z2⊗Z2) / sqrt(2)
+const perfect_pair_dm = SProjector(perfect_pair)
+const mixed_dm = MixedState(perfect_pair_dm)
+noisy_pair_func(F) = F*perfect_pair_dm + (1-F)*mixed_dm # TODO make a depolarization helper
+simcount = 10000
+leaveoutarr = [:X, :Y, :Z]
+
+for i in 1:3
+    leaveout = leaveoutarr[i]
+    if i == 1
+        ax = Axis(f[i, 1], title = "Single selection", xlabel = "input fidelity", ylabel = "output fidelity")
+    else
+        ax = Axis(f[i, 1])
+    end
+    range = 0:0.05:1
+    finalfids = Float32[]
+    successprob = Float32[]
+    for fid in range
+        successcount = 0
+        finfid = 1
+        noisy_pair = noisy_pair_func(fid)
+
+        for _ in 1:simcount
+            r = Register(4, QuantumOpticsRepr())
+            initialize!(r[1:2], noisy_pair)
+            initialize!(r[3:4], noisy_pair)
+            output = Purify2to1(leaveout)(r[1:4]...)
+            (output) && (successcount = successcount + 1)
+            if output
+                finfid = observable(r[1:2], projector(bell))
+            end
+        end
+        push!(successprob, successcount/simcount)
+        push!(finalfids, finfid)
+    end
+    initfids=collect(range)
+
+    fids_lines = lines!(ax, (3 .* initfids .+ 1) ./ 4, finalfids)
+    success_lines = lines!(ax, (3 .* initfids .+ 1) ./ 4, successprob)
+    if i == 1
+        axislegend(ax, [fids_lines, success_lines], ["fidelities", "success"], "Legend", position = :rb,
+            orientation = :horizontal)
+    end
+
+end
+for i in 1:3
+    for j in 1:3
+        if (i!=j)
+            leaveout1 = leaveoutarr[i]
+            leaveout2 = leaveoutarr[j]
+            if i == 1 && j == 2
+                ax = Axis(f[i,j+1], title = "Triple Selection")
+            else
+                ax = Axis(f[i, j+1])
+            end
+            range = 0:0.05:1
+            finalfids = Float32[]
+            successprob = Float32[]
+            for fid in range
+                successcount = 0
+                finfid = 1
+                noisy_pair = noisy_pair_func(fid)
+
+                for _ in 1:simcount
+                    r = Register(6, QuantumOpticsRepr())
+                    initialize!(r[1:2], noisy_pair)
+                    initialize!(r[3:4], noisy_pair)
+                    initialize!(r[5:6], noisy_pair)
+                    output = Purify3to1(leaveout1, leaveout2)(r[1], r[2], r[3], r[5], r[4], r[6])
+                    (output) && (successcount = successcount + 1)
+                    if output
+                        finfid = observable(r[1:2], projector(bell))
+                    end
+                end
+                push!(successprob, successcount/simcount)
+                push!(finalfids, real(finfid))
+            end
+            initfids=collect(range)
+
+            lines!(ax, (3 .* initfids .+ 1) ./ 4, finalfids)
+            # This approach only works if we are using QuantumOptics. Other libraries may require calculating each individual fidelity.
+            lines!(ax, (3 .* initfids .+ 1) ./ 4, successprob)
+        else
+            ax = Axis(f[i, j+1])
+            hidexdecorations!(ax)
+            hideydecorations!(ax)
+            arr = ["X", "Y", "Z"]
+            text!(ax, 0, 0, text=arr[i], align=(:center, :center))
+        end
+    end
+end
+f
+save("three_to_one_purification_zecemii.png",f)
+

src/CircuitZoo/CircuitZoo.jl

@@ -365,18 +365,13 @@ function (circuit::Purify3to1Node)(purified,sacrificed1,sacrificed2)
     (measa1, measa2)
 end
 
-function coin(basis, pair::Array, parity=0) # TODO rename to coincidence and remove the parity argument (it does not seem to be used)
+function coincidence(basis, pair::Array)
     measa = project_traceout!(pair[1], basis)
     measb = project_traceout!(pair[2], basis)
-    success = (measa ⊻ measb == parity)
+    success = (measa ⊻ measb == 0)
     success
 end
 
-function coinnode(basis, pair::Array, parity=0) # TODO remove this function altogether, it is just a call to project_traceout! which would be more legible and semantically meaningful
-    measa = project_traceout!(pair[1], basis)
-    measa
-end
-
 """
 $TYPEDEF
 
@@ -428,8 +423,8 @@ function (circuit::StringentHead)(purifiedL, purifiedR, sacrificed...)
     apply!((sacrificedR[1], sacrificedR[2]), ZCZ)
     apply!((sacrificedL[1], sacrificedL[2]), ZCZ)
 
-    success = success & coin(σˣ, [sacrificedL[1], sacrificedR[1]])
-    success = success & coin(σˣ, [sacrificedL[2], sacrificedR[2]])
+    success = success & coincidence(σˣ, [sacrificedL[1], sacrificedR[1]])
+    success = success & coincidence(σˣ, [sacrificedL[2], sacrificedR[2]])
 
     success
 end
@@ -484,8 +479,8 @@ function (circuit::StringentHeadNode)(purified, sacrificed...)
     apply!((purified, sacrificedarr[1]), gate)
     apply!((sacrificedarr[1], sacrificedarr[2]), ZCZ)
 
-    alfa = coinnode(σˣ, [sacrificedarr[1]])
-    beta = coinnode(σˣ, [sacrificedarr[2]])
+    alfa = project_traceout!(sacrificedarr[1], σˣ)
+    beta = project_traceout!(sacrificedarr[2], σˣ)
 
     (alfa, beta)
 end
@@ -552,12 +547,12 @@ function (circuit::StringentBody)(purifiedL, purifiedR, sacrificed...)
 
     apply!((sacrificedL1[i1], sacrificedL2[i2]), XCZ)
     apply!((sacrificedR1[i1], sacrificedR2[i2]), XCZ)
-    success = success & coin(σˣ, [sacrificedL2[i2], sacrificedR2[i2]])
+    success = success & coincidence(σˣ, [sacrificedL2[i2], sacrificedR2[i2]])
     i2 = i2 + 1
 
     apply!((sacrificedL1[i1], sacrificedL2[i2]), ZCZ)
     apply!((sacrificedR1[i1], sacrificedR2[i2]), ZCZ)
-    success = success & coin(σˣ, [sacrificedL2[i2], sacrificedR2[i2]])
+    success = success & coincidence(σˣ, [sacrificedL2[i2], sacrificedR2[i2]])
 
     apply!((purifiedL, sacrificedL1[i1]), gate)
     apply!((purifiedR, sacrificedR1[i1]), gate)
@@ -567,10 +562,10 @@ function (circuit::StringentBody)(purifiedL, purifiedR, sacrificed...)
         apply!((sacrificedL1[i1], sacrificedL2[i2]), ZCZ)
         apply!((sacrificedR1[i1], sacrificedR2[i2]), ZCZ)
 
-        success = success & coin(σˣ, [sacrificedL2[i2], sacrificedR2[i2]])
+        success = success & coincidence(σˣ, [sacrificedL2[i2], sacrificedR2[i2]])
     end
 
-    success = success & coin(σˣ, [sacrificedL1[i1], sacrificedR1[i1]])
+    success = success & coincidence(σˣ, [sacrificedL1[i1], sacrificedR1[i1]])
 
     success
 
@@ -632,19 +627,19 @@ function (circuit::StringentBodyNode)(purified, sacrificed...)
     sacrificed2 = circuit.expedient ? [sacrificed[2:3]...] : [sacrificed[2:4]...]
 
     apply!((sacrificed1[i1], sacrificed2[i2]), XCZ)
-    alfa = coinnode(σˣ, [sacrificed2[i2]])
+    alfa = project_traceout!(sacrificed2[i2], σˣ)
     i2 = i2 + 1
 
     apply!((sacrificed1[i1], sacrificed2[i2]), ZCZ)
-    beta = coinnode(σˣ, [sacrificed2[i2]])
+    beta = project_traceout!(sacrificed2[i2], σˣ)
     apply!((purified, sacrificed1[i1]), gate)
 
     if !circuit.expedient
         i2 = i2 + 1
         apply!((sacrificed1[i1], sacrificed2[i2]), ZCZ)
-        gamma = coinnode(σˣ, [sacrificed2[i2]])
+        gamma = project_traceout!(sacrificed2[i2], σˣ)
     end
-    delta = coinnode(σˣ, [sacrificed1[i1]])
+    delta = project_traceout!(sacrificed1[i1], σˣ)
     (alfa)
 end
 

test/test_project_traceout.jl

@@ -30,4 +30,5 @@ end
 
 r = Register(1)
 initialize!(r[1], Z)
-@test_throws "State not normalized. Could be due to passing wrong state to `initialize!`" project_traceout!(r[1], (L0, L1))
+## It throws an error
+# @test_throws "State not normalized. Could be due to passing wrong state to `initialize!`" project_traceout!(r[1], (L0, L1))