Two-way Multiplexed protocol "Example"

examples/Project.toml

@@ -8,6 +8,7 @@ GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a"
+QuantumInterface = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5"
 QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
 QuantumSavory = "2de2e421-972c-4cb5-a0c3-999c85908079"
 QuantumSymbolics = "efa7fd63-0460-4890-beb7-be1bbdfbaeae"

examples/twoway_mtp/1_small_visualization.jl

@@ -0,0 +1,32 @@
+include("./setup.jl")
+
+L = 10^1        # Total network distance in km
+η_c = 0.9       # Coupling coefficient
+ϵ_g = 0.001     # Gate error rate
+n = 4           # Number of segments
+
+q = 32          # Number of qubits per interface
+T2 = 1.0        # T2 Dephasing time in seconds
+
+c = 2e5         # Speed of light in km/s
+l_att = 20      # Attenuation length in km
+
+l0 = L / n      # Internode distance in km
+p_ent = 0.5*η_c^2*exp(-l0/l_att)    # Entanglement generation probability
+ξ = 0.25ϵ_g     # Measurement error rate
+F = 1-1.25ϵ_g   # Initial bellpair fidelity
+
+t_comms = fill(l0/c, n)     # Internode communication time in seconds
+distil_sched = distil_scheds[(L, n, ϵ_g, η_c)]  # Distillation schedule
+
+
+net_param = NetworkParam(n, q; T2, F, p_ent, ϵ_g, ξ, t_comms, distil_sched)
+network = Network(net_param, rng=Xoshiro(1))
+
+# @time simulate!(network)
+
+sim = simulate_iter!(network)
+video_path = "./results/1_small_visualization.mp4"
+record(network.fig, video_path, 1:nsteps(net_param); framerate=2) do i
+    iterate(sim, nothing)
+end

examples/twoway_mtp/2_small_withdistil.jl

@@ -0,0 +1,33 @@
+include("./setup.jl")
+
+L = 10^1        # Total network distance in km
+η_c = 0.9       # Coupling coefficient
+ϵ_g = 0.001     # Gate error rate
+n = 4           # Number of segments
+
+q = 32          # Number of qubits per interface
+T2 = 1.0        # T2 Dephasing time in seconds
+
+c = 2e5         # Speed of light in km/s
+l_att = 20      # Attenuation length in km
+
+l0 = L / n      # Internode distance in km
+p_ent = 0.5*η_c^2*exp(-l0/l_att)    # Entanglement generation probability
+ξ = 0.25ϵ_g     # Measurement error rate
+F = 1-1.25ϵ_g   # Initial bellpair fidelity
+
+t_comms = fill(l0/c, n)     # Internode communication time in seconds
+# distil_sched = distil_scheds[(L, n, ϵ_g, η_c)]  # Distillation schedule
+distil_sched = fill(true, floor(Int, log2(n)))
+
+
+net_param = NetworkParam(n, q; T2, F, p_ent, ϵ_g, ξ, t_comms, distil_sched)
+network = Network(net_param, rng=Xoshiro(1))
+
+# @time simulate!(network)
+
+sim = simulate_iter!(network)
+video_path = "./results/2_small_withdistil.mp4"
+record(network.fig, video_path, 1:nsteps(net_param); framerate=2) do i
+    iterate(sim, nothing)
+end

examples/twoway_mtp/3_large_network.jl

@@ -0,0 +1,32 @@
+include("./setup.jl")
+
+L = 10^4        # Total network distance in km
+η_c = 0.9       # Coupling coefficient
+ϵ_g = 0.001     # Gate error rate
+n = 512           # Number of segments
+
+q = 1024          # Number of qubits per interface
+T2 = 1.0        # T2 Dephasing time in seconds
+
+c = 2e5         # Speed of light in km/s
+l_att = 20      # Attenuation length in km
+
+l0 = L / n      # Internode distance in km
+p_ent = 0.5*η_c^2*exp(-l0/l_att)    # Entanglement generation probability
+ξ = 0.25ϵ_g     # Measurement error rate
+F = 1-1.25ϵ_g   # Initial bellpair fidelity
+
+t_comms = fill(l0/c, n)     # Internode communication time in seconds
+distil_sched = distil_scheds[(L, n, ϵ_g, η_c)]  # Distillation schedule
+
+
+net_param = NetworkParam(n, q; T2, F, p_ent, ϵ_g, ξ, t_comms, distil_sched)
+network = Network(net_param, rng=Xoshiro(1))
+
+# @time simulate!(network)
+
+sim = simulate_iter!(network)
+video_path = "./results/3_large_network.mp4"
+record(network.fig, video_path, 1:nsteps(net_param); framerate=2) do i
+    iterate(sim, nothing)
+end

examples/twoway_mtp/README.md

@@ -0,0 +1,37 @@
+# Two-way Multiplexed Protocol
+
+This example demonstrates a high-performance simulation of a multiplexed two-way quantum repeater architecture based on the protocol outlined in [Mantri et al., 2024](https://arxiv.org/abs/2409.06152). It models entanglement generation, DEJMPS purification on a schedule, and entanglement swapping under realistic noise, leveraging QuantumSavory.jl’s flexible architecture.
+
+## Overview
+This simulation models a linear quantum network consisting of:
+- A sender (Alice), a receiver (Bob), and one or more intermediate Repeaters
+- Multiplexed imperfect Bell pair generation on each link
+- A realistic noise model involving T2 decoherence, gate errors and measurement infidelity
+
+The experiment aims to evaluate:
+- End-to-end fidelity after purification and swapping
+- Secret Key Rate (SKR) across different configurations
+
+![SimGIF](https://pouch.jumpshare.com/preview/DC3f7WLV8MZBij8oZZJzsy3wfS5RCLIgpxfLik0SPj12KPPdBe-5SZOvvPgGz_iPxHIA3sErGJ_1XUOG1nWkWrem2COdyPx78xsPzwhcFZA)
+
+## Running the simulation
+Clone the repository and install the necessary dependencies:
+```sh
+git clone https://github.com/QuantumSavory/QuantumSavory.jl.git
+julia --project=examples -e "using Pkg; Pkg.instantiate()"
+```
+
+Run one of the example simulations:
+```sh
+julia --project=examples ./examples/twoway_mtp/n_example.jl
+```
+
+## References
+- Mantri, P., Goodenough, K., & Towsley, D. (2024, September 10). Comparing one- and two-way quantum repeater architectures. arXiv.org. https://arxiv.org/abs/2409.06152
+- Main project repository and Experiment results: http://github.com/sagnikpal2004/QNet-MTP
+
+## To Do List:
+- Upstream [`RGate`](baseops/RGate.jl) to QuantumSymbolics.jl - in progress  [QuantumSymbolics.jl:Pull#95](https://github.com/QuantumSavory/QuantumSymbolics.jl/pull/95)
+- [`DEJMPSProtocol`](noisyops/CircuitZoo.jl) is also missing from QuantumSavory.jl - in progress [QuantumSavory.jl:Issue#237](https://github.com/QuantumSavory/QuantumSavory.jl/issues/237)
+- Modify measurement error model to integrate into the projectors themselves instead of using a random number to flip the measurement using [`project_traceout!`](noisyops/traceout.jl) (Do we need a `project_traceout!` that can work with POVM as well and not just Kets?)
+- Figure out a better way to do noisy [`apply!`](noisyops/apply.jl)

examples/twoway_mtp/baseops/RGate.jl

@@ -0,0 +1,25 @@
+import QuantumSymbolics: Metadata
+
+QuantumSymbolics.@withmetadata struct RGate <: QuantumSymbolics.AbstractSingleQubitGate
+    dir::Symbol
+    θ::Float64
+end
+QuantumSymbolics.symbollabel(g::RGate) = "R$(g.dir)($(g.θ))"
+QuantumSymbolics.ishermitian(::RGate) = true
+QuantumSymbolics.isunitary(::RGate) = true
+
+QuantumSymbolics.express_nolookup(gate::RGate, ::QuantumSymbolics.QuantumOpticsRepr) = QuantumOptics.Operator(
+    QuantumInterface.SpinBasis(1//2),
+    if gate.dir == :x
+        [cos(gate.θ/2) -im*sin(gate.θ/2); -im*sin(gate.θ/2) cos(gate.θ/2)]
+    elseif gate.dir == :y
+        [cos(gate.θ/2) -sin(gate.θ/2); sin(gate.θ/2) cos(gate.θ/2)]
+    elseif gate.dir == :z
+        [exp(-im*gate.θ/2) 0; 0 exp(im*gate.θ/2)]
+    end
+)
+
+
+Rx(θ::Float64) = RGate(:x, θ)
+Ry(θ::Float64) = RGate(:y, θ)
+Rz(θ::Float64) = RGate(:z, θ)

examples/twoway_mtp/baseops/uptotime.jl

@@ -0,0 +1,19 @@
+"""Updates the time of a register"""
+function uptotime!(reg::QuantumSavory.Register, t::Float64)
+    for i in 1:length(reg.traits)
+        QuantumSavory.uptotime!(reg[i], t)
+    end
+end
+
+
+"""Updates the time of the network"""
+function uptotime!(N::Network, t::Float64)
+    for node in N.nodes
+        if !isnothing(node.left)
+            uptotime!(node.left, t)
+        end
+        if !isnothing(node.right)
+            uptotime!(node.right, t)
+        end
+    end
+end

examples/twoway_mtp/noisyops/CircuitZoo.jl

@@ -0,0 +1,78 @@
+import QuantumSavory
+include("./apply.jl")
+include("./traceout.jl")
+include("../baseops/RGate.jl")
+
+
+_b2 = SpinBasis(1//2)
+_l0 = spinup(_b2)
+_l1 = spindown(_b2)
+_l⁺ = (_l0 + _l1)/√2
+_l⁻ = (_l0 - _l1)/√2
+_l00 = projector(_l0)
+_l11 = projector(_l1)
+_l⁺⁺ = projector(_l⁺)
+_l⁻⁻ = projector(_l⁻)
+
+noisy_zbasis(ξ) = [(1-ξ)*_l00 + ξ*_l11, (1-ξ)*_l11 + ξ*_l00]
+noisy_xbasis(ξ) = [(1-ξ)*_l⁺⁺ + ξ*_l⁻⁻, (1-ξ)*_l⁻⁻ + ξ*_l⁺⁺]
+
+
+struct EntanglementSwap <: QuantumSavory.CircuitZoo.AbstractCircuit
+    ϵ_g::Float64
+    ξ::Float64
+    rng::AbstractRNG
+
+    function EntanglementSwap(ϵ_g::Float64, ξ::Float64, rng::AbstractRNG=Random.GLOBAL_RNG)
+        @assert 0 <= ϵ_g <= 1   "ϵ_g must be in [0, 1]"
+        @assert 0 <=  ξ  <= 1   "ξ must be in [0, 1]"
+
+        new(ϵ_g, ξ, rng)
+    end
+end
+function (circuit::EntanglementSwap)(localL, remoteL, localR, remoteR)
+    apply!((localL, localR), QuantumSavory.CNOT; ϵ_g=circuit.ϵ_g)
+    xmeas = project_traceout!(localL, noisy_xbasis(circuit.ξ))
+    zmeas = project_traceout!(localR, noisy_zbasis(circuit.ξ))
+    if xmeas==2
+        QuantumSavory.apply!(remoteL, QuantumSavory.Z)
+    end
+    if zmeas==2
+        QuantumSavory.apply!(remoteR, QuantumSavory.X)
+    end
+    return xmeas, zmeas
+end
+inputqubits(::EntanglementSwap) = 4
+
+
+struct DEJMPSProtocol <: QuantumSavory.CircuitZoo.AbstractCircuit
+    ϵ_g::Float64
+    ξ::Float64
+
+    function DEJMPSProtocol(ϵ_g::Float64, ξ::Float64)
+        @assert 0 <= ϵ_g <= 1   "ϵ_g must be in [0, 1]"
+        @assert 0 <=  ξ  <= 1   "ξ must be in [0, 1]"
+
+        new(ϵ_g, ξ)
+    end
+end
+function (circuit::DEJMPSProtocol)(purifiedL, purifiedR, sacrificedL, sacrificedR)
+    QuantumSavory.apply!(purifiedL, Rx(π/2))
+    QuantumSavory.apply!(sacrificedL, Rx(π/2))
+    QuantumSavory.apply!(purifiedR, Rx(-π/2))
+    QuantumSavory.apply!(sacrificedR, Rx(-π/2))
+
+    apply!([purifiedL, sacrificedL], QuantumSavory.CNOT; ϵ_g=circuit.ϵ_g)
+    apply!([purifiedR, sacrificedR], QuantumSavory.CNOT; ϵ_g=circuit.ϵ_g)
+
+    measa = project_traceout!(sacrificedL, noisy_zbasis(circuit.ξ))
+    measb = project_traceout!(sacrificedR, noisy_zbasis(circuit.ξ))
+
+    success = measa == measb
+    if !success
+        QuantumSavory.traceout!(purifiedL)
+        QuantumSavory.traceout!(purifiedR)
+    end
+    return success
+end
+inputqubits(::DEJMPSProtocol) = 4

examples/twoway_mtp/noisyops/apply.jl

@@ -0,0 +1,87 @@
+function noise(state::QuantumOptics.Operator, indices)
+    mixed_state = QuantumSymbolics.express(QuantumSavory.IdentityOp(QuantumOptics.basis(state)) / length(QuantumOptics.basis(state)))
+
+    if !isa(QuantumOptics.basis(state), QuantumInterface.CompositeBasis)
+        return mixed_state
+    elseif length(indices) == length(QuantumOptics.basis(state).bases)
+        return mixed_state
+    else
+        mixed_basis = QuantumInterface.CompositeBasis(QuantumOptics.basis(state).bases[indices])
+        mixed_state = QuantumSymbolics.express(QuantumSavory.IdentityOp(mixed_basis) / length(mixed_basis))
+        traced_state = QuantumOptics.ptrace(state, indices)
+
+        bit_order = vcat([i for i in 1:length(QuantumOptics.basis(state).bases) if !(i in indices)], indices)
+        perm = [findfirst(==(x), bit_order) for x in 1:length(bit_order)]
+
+        perfect_state = QuantumOptics.:⊗(traced_state, mixed_state)
+        noisy_state = QuantumOptics.permutesystems(perfect_state, perm)
+
+        return noisy_state
+    end
+end
+
+function apply!(state::QuantumOptics.Operator, indices, operation::QuantumOptics.Operator; ϵ_g::Float64=0.0)
+    op = QuantumOpticsBase.is_apply_shortcircuit(state, indices, operation) ? operation : QuantumOptics.embed(QuantumOptics.basis(state), indices, operation)
+    state.data = ((1-ϵ_g)*(op*state*op') + ϵ_g*noise(state, indices)).data
+    return state
+end
+function apply!(state::QuantumOptics.Ket, indices, operation::QuantumOptics.Operator; ϵ_g::Float64=0.0)
+    ϵ_g > 0.0 && return apply!(QuantumOptics.dm(state), indices, operation; ϵ_g)
+
+    op = QuantumOpticsBase.is_apply_shortcircuit(state, indices, operation) ? operation : QuantumOptics.embed(QuantumOptics.basis(state), indices, operation)
+    state.data = (op*state).data
+    return state
+end
+apply!(state::QuantumOptics.Ket, indices, operation::T; ϵ_g::Float64=0.0) where {T<:QuantumInterface.AbstractSuperOperator} = apply!(QuantumInterface.dm(state), indices, operation; ϵ_g)
+# function apply!(state::QuantumOptics.Operator, indices, operation::T; ϵ_g::Float64=0.0) where {T<:QuantumInterface.AbstractSuperOperator}
+#     if QuantumOpticsBase.is_apply_shortcircuit(state, indices, operation)
+#         state.data = (operation*state).data
+#         return state
+#     else
+#         error("`apply!` does not yet support QuantumOptics.embedding superoperators acting only on a subsystem of the given state")
+#     end
+# end
+
+
+function apply!(state, indices::Base.AbstractVecOrTuple{Int}, operation::QuantumSymbolics.Symbolic{QuantumInterface.AbstractOperator}; ϵ_g::Float64=0.0)
+    repr = QuantumSavory.default_repr(state)
+    apply!(state, indices, QuantumSymbolics.express(operation, repr, QuantumSymbolics.UseAsOperation()); ϵ_g)
+end
+function apply!(state, indices::Base.AbstractVecOrTuple{Int}, operation::QuantumSymbolics.Symbolic{QuantumInterface.AbstractSuperOperator}; ϵ_g::Float64=0.0)
+    repr = QuantumSavory.default_repr(state)
+    apply!(state, indices, QuantumSymbolics.express(operation, repr, QuantumSymbolics.UseAsOperation()); ϵ_g)
+end
+function apply!(regs::Vector{QuantumSavory.Register}, indices::Base.AbstractVecOrTuple{Int}, operation::QuantumSymbolics.Symbolic{QuantumInterface.AbstractOperator}; time=nothing, ϵ_g::Float64=0.0)
+    invoke(apply!, Tuple{Vector{QuantumSavory.Register}, Base.AbstractVecOrTuple{Int}, Any}, regs, indices, operation; time, ϵ_g)
+end
+function apply!(regs::Vector{QuantumSavory.Register}, indices::Base.AbstractVecOrTuple{Int}, operation::QuantumSymbolics.Symbolic{QuantumInterface.AbstractSuperOperator}; time=nothing, ϵ_g::Float64=0.0)
+    invoke(apply!, Tuple{Vector{QuantumSavory.Register}, Base.AbstractVecOrTuple{Int}, Any}, regs, indices, operation; time, ϵ_g)
+end
+
+
+"""
+Apply a given operation on the given set of register slots.
+
+`apply!([regA, regB], [slot1, slot2], Gates.CNOT)` would apply a CNOT gate
+on the content of the given registers at the given slots.
+The appropriate representation of the gate is used,
+depending on the formalism under which a quantum state is stored in the given registers.
+The Hilbert spaces of the registers are automatically joined if necessary.
+"""
+function apply!(regs::Vector{QuantumSavory.Register}, indices::Base.AbstractVecOrTuple{Int}, operation; time=nothing, ϵ_g::Float64=0.0)
+    max_time = maximum((r.accesstimes[i] for (r,i) in zip(regs,indices)))
+    if !isnothing(time)
+        time<max_time && error("The simulation was commanded to apply $(operation) at time t=$(time) although the current simulation time is higher at t=$(max_time). Consider using locks around the offending operations.")
+        max_time = time
+    end
+    QuantumSavory.uptotime!(regs, indices, max_time)
+    QuantumSavory.subsystemcompose(regs,indices)
+    state = regs[1].staterefs[indices[1]].state[]
+    state_indices = [r.stateindices[i] for (r,i) in zip(regs, indices)]
+    state = apply!(state, state_indices, operation; ϵ_g)
+    regs[1].staterefs[indices[1]].state[] = state
+    return regs, max_time
+end
+apply!(refs::Vector{QuantumSavory.RegRef}, operation; time=nothing, ϵ_g::Float64=0.0) = apply!([r.reg for r in refs], [r.idx for r in refs], operation; time, ϵ_g)
+apply!(refs::NTuple{N,QuantumSavory.RegRef}, operation; time=nothing, ϵ_g::Float64=0.0) where {N} = apply!([r.reg for r in refs], [r.idx for r in refs], operation; time, ϵ_g) # TODO temporary array allocated here
+apply!(ref::QuantumSavory.RegRef, operation; time=nothing, ϵ_g::Float64=0.0) = apply!([ref.reg], [ref.idx], operation; time, ϵ_g)