Implementing Cutoff protocol, related synchronization communication t…

…ags, and general Improvements to protocol stability (#120)


Co-authored-by: Abhishek Bhatt <[email protected]>
Co-authored-by: Stefan Krastanov <[email protected]>
Co-authored-by: Stefan Krastanov <[email protected]>

CHANGELOG.md

@@ -1,10 +1,18 @@
 # News
 
-## v0.5.0 - 2024-08-11
 
-- `observable` now takes a default value as a kwarg, i.e., you need to make the substitution `observable(regs, obs, 0.0; time)` ↦ `observable(regs, obs; something=0.0, time)`
+## v0.5.0 - 2024-09-05
+
+- Develop `CutoffProt` to deal with deadlocks in a simulation
+- Expand `SwapperProt` with `agelimit` to permit cutoff policies (with `CutoffProt`)
+- Tutorial and interactive examples for entanglement distribution on a grid with local-only knowledge
+- **(breaking)** `observable` now takes a default value as a kwarg, i.e., you need to make the substitution `observable(regs, obs, 0.0; time)` ↦ `observable(regs, obs; something=0.0, time)`
 - Bump QuantumSymbolics and QuantumOpticsBase compat bound and bump julia compat to 1.10.
 
+## v0.4.2 - 2024-08-13
+
+- Incorrect breaking release. It should have been 0.5 (see above).
+
 ## v0.4.1 - 2024-06-05
 
 - Significant improvements to the performance of `query`.

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumSavory"
 uuid = "2de2e421-972c-4cb5-a0c3-999c85908079"
 authors = ["Stefan Krastanov <[email protected]>"]
-version = "0.4.2"
+version = "0.5"
 
 [deps]
 Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"

docs/make.jl

@@ -4,6 +4,7 @@ push!(LOAD_PATH,"../src/")
 using Documenter
 using DocumenterCitations
 using QuantumSavory
+using QuantumSavory.ProtocolZoo # TODO is this the correct place to place this to ensure cross_references work
 
 DocMeta.setdocmeta!(QuantumSavory, :DocTestSetup, :(using QuantumSavory); recursive=true)
 

docs/src/howto/repeatergrid/repeatergrid.md

@@ -0,0 +1,173 @@
+# [Entanglement Generation On A Repeater Grid](@id Entanglement-Generation-On-A-Repeater-Grid)
+
+```@meta
+DocTestSetup = quote
+    using QuantumSavory
+end
+```
+
+This section provides a detailed walkthrough of how QuantumSavory.jl can be used to simulate entanglement generation on a network of repeaters where each repeater relies only on local knowledge of the network.
+This is only one of many ways in which such a network can be set up, focusing on one particular "no global knowledge" approach.
+
+A complete implementation of the simulation described here is available in the `examples/repeatergrid` folder of the `QuantumSavory` repository.
+
+For this example, we consider a square grid topology in which each node is connected to its nearest neighbors.
+The registers act as repeater nodes. The nodes on the diagonal corners are Alice and Bob, the two special nodes that the network is trying to entangle through generating link-level entanglement at each edge and performing appropriate swaps at each node.
+
+The goal is to establish entanglement between Alice and Bob by routing entanglement through any of the possible paths (horizontal or vertical) formed by local entanglement links and then swapping those links by performing entanglement swaps.
+
+This employs functionality from the `ProtocolZoo` module of QuantumSavory to run the following Quantum Networking protocols:
+
+- [`EntanglerProt`](@ref): Entangler protocol to produce link-level entanglement at each edge in the network
+
+- [`SwapperProt`](@ref): Swapper protocol runs at each node except at the Alice and Bob nodes, to perform swaps. The swaps are performed only if a query deems them useful for propagating entanglement closer and closer to Alice and Bob.
+
+- [`EntanglementTracker`](@ref): Entanglement Tracker protocol to keep track of and update the local link state (classical knowledge) by listening for "entanglement update" messages generated by the other protocols (`SwapperProt` & `CutoffProt` specifically).
+
+- [`CutoffProt`](@ref): As the simulation progresses, the unused entangled pairs generated by `EntanglerProt` need to be discarded due to the loss of fidelity under noise as they might not be suitable for networking applications beyond a certain cutoff interval of time from their instantiation. The `CutoffProt` is instantiated with a `retention_time` parameter that discards such qubits in each node.
+
+- [`EntanglementConsumer`](@ref): This protocol runs between the end nodes and consumes the final entanglement generated as a result of all of the above protocols, which is supposed to represent the qubits being consumed in a networking application.
+
+All of the above protocols rely on the query and tagging functionality as described in the [Tagging and Querying](@ref tagging-and-querying) section.
+
+Other than that, `ConcurrentSim` and `ResumableFunctions` are used in the backend to run the discrete event simulation. `Graphs` helps with some functionality needed for the `RegisterNet` datastructure that forms the grid. `GLMakie` and `NetworkLayout` are used for visualization along with the visualization functionality implemented in `QuantumSavory` itself.
+
+# Custom Predicate And Choosing function
+
+```julia
+function check_nodes(net, c_node, node; low=true)
+    n = Int(sqrt(size(net.graph)[1])) # grid size
+    c_x = c_node%n == 0 ? c_node ÷ n : (c_node ÷ n) + 1
+    c_y = c_node - n*(c_x-1)
+    x = node%n == 0 ? node ÷ n : (node ÷ n) + 1
+    y = node - n*(x-1)
+    return low ? (c_x - x) >= 0 && (c_y - y) >= 0 : (c_x - x) <= 0 && (c_y - y) <= 0
+end
+```
+
+The Swapper Protocol is initialized with a custom predicate function which is then placed in a call to `queryall` inside the Swapper to pick the nodes that are suitable to perform a swap with. The criteria for "suitability" is described below.
+
+This predicate function encodes most of the "logic" a local node will be performing.
+
+The custom predicate function shown above is parametrized with `net` and `c_node` along with the keyword argument `low`, when initializing the Swapper Protocol. This predicate function `Int->Bool` selects the target remote nodes for which a swap is appropriate. The arguments are:
+
+- `net`: The network of register nodes representing the graph structure, an instance of `RegisterNet`.
+
+- `c_node`: The node in which the Swapper protocol would be running.
+
+- `node`: As the [`queryall`](@ref) function goes through all the nodes linked with the current node, the custom predicate filters them depending on whether the node is suitable for a swap or not.
+
+- `low`: The nodes in the grid are numbered as consecutive integers starting from 1. If the Swapper is running at some node n, we want a link closest to Alice and another closest to Bob to perform a swap. We communicate whether we are looking for nodes of the first kind or the latter with the `low` keyword.
+
+Out of all the links at some node, the suitable ones are picked by computing the difference between the coordinates of the current node with the coordinates of the candidate node. A `low` node should have both of the `x` and `y` coordinate difference positive and vice versa for a non-`low` node.
+
+As the Swapper gets a list of suitable candidates for a swap in each direction, the one with the furthest distance from the current node is chosen by summing the x distance and y-distance.
+
+```julia
+function choose_node(net, node, arr; low=true)
+    grid_size = Int(sqrt(size(net.graph)[1]))
+    return low ? argmax((distance.(grid_size, node, arr))) : argmin((distance.(grid_size, node, arr)))
+end
+
+function distance(n, a, b)
+    x1 = a%n == 0 ? a ÷ n : (a ÷ n) + 1
+    x2 = b%n == 0 ? b ÷ n : (b ÷ n) + 1
+    y1 = a - n*(x1-1)
+    y2 = b - n*(x2-1)
+
+    return x1 - x2 + y1 - y2
+end
+```
+
+# Simulation and Visualization
+
+```julia
+n = 6 # the size of the square grid network (n × n)
+regsize = 8 # the size of the quantum registers at each node
+
+graph = grid([n,n])
+
+net = RegisterNet(graph, [Register(regsize, fill(5.0, regsize)) for i in 1:n^2])
+
+sim = get_time_tracker(net)
+
+# each edge is capable of generating raw link-level entanglement
+for (;src, dst) in edges(net)
+    eprot = EntanglerProt(sim, net, src, dst; rounds=-1, randomize=true)
+    @process eprot()
+end
+
+# each node except the corners on one of the diagonals is capable of swapping entanglement
+for i in 2:(n^2 - 1)
+    l(x) = check_nodes(net, i, x)
+    h(x) = check_nodes(net, i, x; low=false)
+    cL(arr) = choose_node(net, i, arr)
+    cH(arr) = choose_node(net, i, arr; low=false)
+    swapper = SwapperProt(sim, net, i; nodeL = l, nodeH = h, chooseL = cL, chooseH = cH, rounds=-1)
+    @process swapper()
+end
+
+for v in vertices(net)
+    tracker = EntanglementTracker(sim, net, v)
+    @process tracker()
+end
+
+# Entanglement usage/consumption by the network end nodes
+
+consumer = EntanglementConsumer(sim, net, 1, n^2)
+@process consumer()
+
+# decoherence protocol runs at each node to free up slots that haven't been used past the retention time
+for v in vertices(net)
+    decprot = DecoherenceProt(sim, net, v)
+    @process decprot()
+end
+```
+
+We set up the simulation to run with a 6x6 grid of nodes above. Here, each node has 8 qubit slots.
+Each vertical and horizontal edge runs an entanglement generation protocol.
+Each node in the network runs an entanglement tracker protocol and all of the nodes except the nodes that we're trying to connect,
+i.e., Alice's and Bob's nodes which are at the diagonal ends of the grid run the swapper protocol.
+The code that runs and visualizes this simulation is shown below
+
+```julia
+fig = Figure(;size=(600, 600))
+
+# the network part of the visualization
+layout = SquareGrid(cols=:auto, dx=10.0, dy=-10.0)(graph) # provided by NetworkLayout, meant to simplify plotting of graphs in 2D
+_, ax, _, obs = registernetplot_axis(fig[1:2,1], net;registercoords=layout)
+
+# the performance log part of the visualization
+entlog = Observable(consumer.log) # Observables are used by Makie to update the visualization in real-time in an automated reactive way
+ts = @lift [e[1] for e in $entlog]  # TODO this needs a better interface, something less cluncky, maybe also a whole Makie recipe
+tzzs = @lift [Point2f(e[1],e[2]) for e in $entlog]
+txxs = @lift [Point2f(e[1],e[3]) for e in $entlog]
+Δts = @lift length($ts)>1 ? $ts[2:end] .- $ts[1:end-1] : [0.0]
+entlogaxis = Axis(fig[1,2], xlabel="Time", ylabel="Entanglement", title="Entanglement Successes")
+ylims!(entlogaxis, (-1.04,1.04))
+stem!(entlogaxis, tzzs)
+histaxis = Axis(fig[2,2], xlabel="ΔTime", title="Histogram of Time to Successes")
+hist!(histaxis, Δts)
+
+display(fig)
+
+step_ts = range(0, 200, step=0.1)
+record(fig, "grid_sim6x6hv.mp4", step_ts; framerate=10, visible=true) do t
+    run(sim, t)
+    notify.((obs,entlog))
+    ylims!(entlogaxis, (-1.04,1.04))
+    xlims!(entlogaxis, max(0,t-50), 1+t)
+    autolimits!(histaxis)
+end
+
+```
+
+# Result
+
+```@repl
+include("../../../../examples/repeatergrid/1a_async_interactive_visualization.jl") # hide
+```
+
+```@raw html
+<video src="../grid_sim6x6hv.mp4" autoplay loop muted></video>
+```

docs/src/symbolics.md

@@ -202,5 +202,4 @@ express(MixedState(X1)/2+SProjector(Z1)/2, CliffordRepr())
 ```
 
 !!! warning "Stabilizer state expressions"
-
     The state written as $\frac{|Z₁⟩⊗|Z₁⟩+|Z₂⟩⊗|Z₂⟩}{√2}$ is a well known stabilizer state, namely a Bell state. However, automatically expressing it as a stabilizer is a prohibitively expensive computational operation in general. We do not perform that computation automatically. If you want to ensure that states you define can be automatically converted to tableaux for Clifford simulations, avoid using summation of kets. On the other hand, in all of our Clifford Monte-Carlo simulations, `⊗` is fully supported, as well as [`SProjector`](@ref), [`MixedState`](@ref), [`StabilizerState`](@ref), and summation of density matrices.

docs/src/visualizations.md

@@ -18,6 +18,7 @@ The [`registernetplot_axis`](@ref) function can be used to draw a given set of r
 
 ```@example vis
 using GLMakie
+GLMakie.activate!()
 using QuantumSavory
 
 # create a network of qubit registers
@@ -57,17 +58,21 @@ initialize!(network[1,1]) # hide
 initialize!(network[2,3], X₁) # hide
 initialize!((network[3,1],network[4,2]), X₁⊗Z₂) # hide
 apply!((network[2,3],network[3,1]), CNOT) # hide
-fig = Figure(size=(400,400)) # hide
-_, _, plt, obs = registernetplot_axis(fig[1,1],network) # hide
+fig = Figure(size=(700,400)) # hide
+_, ax, plt, obs = registernetplot_axis(fig[1,1],network) # hide
+fig
+```
+
+```@example vis
 QuantumSavory.showmetadata(fig,ax,plt,1,1)
 fig
 ```
 
 And here with some extra tag metadata.
 
 ```@example vis
-tag!(net[2,3], :specialplace, 1, 2)
-tag!(net[2,3], :otherdata, 3, 4)
+tag!(network[2,3], :specialplace, 1, 2)
+tag!(network[2,3], :otherdata, 3, 4)
 QuantumSavory.showmetadata(fig,ax,plt,2,3)
 fig
 ```

examples/congestionchain/1_visualization.jl

@@ -45,13 +45,16 @@ scatter!(ax_fidZZ,ts,fidZZ,label="ZZ",color=(c2,0.1))
 display(fig)
 
 step_ts = range(0, 1000, step=0.1)
-for t in step_ts
+
+record(fig, "congestionchain.mp4", step_ts; framerate=50, visible=true) do t
     run(sim, t)
     ax.title = "t=$(t)"
-    notify(obs)
-    notify(ts)
-    autolimits!(ax_fidXX)
-    autolimits!(ax_fidZZ)
+    if length(ts[])>2 # to avoid failing autolimits on empty plots
+        notify(obs)
+        notify(ts)
+        autolimits!(ax_fidXX)
+        autolimits!(ax_fidZZ)
+    end
 end
 
 ##

examples/repeatergrid/1a_async_interactive_visualization.jl

@@ -0,0 +1,79 @@
+using GLMakie
+# TODO significant code duplication with the other examples
+
+include("setup.jl")
+
+sim, net, graph, consumer, params... = prepare_simulation()
+
+fig = Figure(;size=(1200, 1100))
+
+# the network part of the visualization
+layout = SquareGrid(cols=:auto, dx=30.0, dy=-30.0)(graph) # provided by NetworkLayout, meant to simplify plotting of graphs in 2D
+_, ax, _, obs = registernetplot_axis(fig[1:2,1], net;registercoords=layout)
+
+# the performance log part of the visualization
+entlog = Observable(consumer.log) # Observables are used by Makie to update the visualization in real-time in an automated reactive way
+ts = @lift [e[1] for e in $entlog]  # TODO this needs a better interface, something less cluncky, maybe also a whole Makie recipe
+tzzs = @lift [Point2f(e[1],e[2]) for e in $entlog]
+txxs = @lift [Point2f(e[1],e[3]) for e in $entlog]
+Δts = @lift length($ts)>1 ? $ts[2:end] .- $ts[1:end-1] : [0.0]
+entlogaxis = Axis(fig[1,2], xlabel="Time", ylabel="Entanglement", title="Entanglement Successes")
+ylims!(entlogaxis, (-1.04,1.04))
+stem!(entlogaxis, txxs)
+histaxis = Axis(fig[2,2], xlabel="ΔTime", title="Histogram of Time to Successes")
+hist!(histaxis, Δts)
+
+avg_fids = @lift cumsum([e[3] for e in $entlog])./cumsum(ones(length($entlog))) #avg fidelity per unit time
+fid_info = @lift [Point2f(t,f) for (t,f) in zip($ts, $avg_fids)]
+fid_axis = Axis(fig[3,1], xlabel="Time", ylabel="Avg. Fidelity", title="Time evolution of Average Fidelity")
+lines!(fid_axis, fid_info)
+
+num_epr = @lift cumsum(ones(length($entlog)))./($ts) #avg number of pairs per unit time
+num_epr_info = @lift [Point2f(t,n) for (t,n) in zip($ts, $num_epr)]
+num_epr_axis = Axis(fig[3,2], xlabel="Time", title="Avg. Number of Entangled Pairs between Alice and Bob")
+lines!(num_epr_axis, num_epr_info)
+
+
+#  sliders
+sg = SliderGrid( # TODO significant code duplication with the other examples
+    fig[4,:],
+    (label="Probability of success of Entanglement generation at each attempt",
+        range=0.001:0.05:1.0, format="{:.3f}", startvalue=0.001),
+    (label="Local busy time for swapper",
+        range=0.001:0.5:10.0, format="{:.3f}", startvalue=0.001),
+    (label="Wait time after failure to lock qubits for a swap",
+        range=0.1:0.05:1.0, format="{:.2f}", startvalue=0.1),
+    (label="Retention time for an unused qubit",
+        range=0.1:0.1:10.0, format="{:.2f}", startvalue=5.0),
+    (label="Time before a qubit's retention time runs out (for `agelimit`)",
+        range=0.1:0.5:10.0, format="{:.2f}", startvalue=0.5),
+    (label="Period of time between subsequent queries at the consumer",
+        range=0.001:0.05:1.0, format="{:.3f}", startvalue=0.001),
+    (label="Period of time between subsequent queries at the DecoherenceProtocol",
+        range=0.001:0.05:1.0, format="{:.3f}", startvalue=0.001),
+
+    width = 600,
+    tellheight = false)
+
+for (param, slider) in zip(params, sg.sliders)
+    on(slider.value) do val
+        param[] = val
+    end
+end
+
+
+display(fig)
+
+step_ts = range(0, 50, step=0.1)
+record(fig, "grid_sim6x6hv.mp4", step_ts; framerate=10, visible=true) do t
+    run(sim, t)
+    notify.((obs,entlog))
+    notify.(params)
+    ylims!(entlogaxis, (-1.04,1.04))
+    xlims!(entlogaxis, max(0,t-50), 1+t)
+    ylims!(fid_axis, (0, 1.04))
+    xlims!(fid_axis, max(0, t-50), 1+t)
+    autolimits!(histaxis)
+    ylims!(num_epr_axis, (0, 4))
+    xlims!(num_epr_axis, max(0, t-50), 1+t)
+end