r"""
.. _borealis_advanced:
Operating Borealis --- Advanced
===============================
.. raw:: html
Borealis will only be available to the public until June 2, 2023.
Public access to this device will not be available after this date.
*Authors: Fabian Laudenbach and Theodor Isacsson*
.. related::
run_advanced_time_domain Advanced time-domain photonic circuits
run_time_domain Time-domain photonic circuits
tutorial_borealis_beginner Operating Borealis --- Beginner
tutorial_borealis_quickstart Borealis Quickstart
.. note::
If you're in a hurry, check out our :doc:`Borealis Quickstart `.
In the :doc:`Borealis beginner tutorial ` we
described how Xanadu's quantum-advantage machine `Borealis
`__ can be programmed in a fairly
straightforward way, using convenience functions to make the user-defined
circuit compatible with Borealis. In this advanced tutorial, we will again
define and submit an instance of Gaussian Boson Sampling (GBS), but this time we
will instead define the circuit and its parameters manually, without using any
of the utility functions. This will allow us to better understand the subtleties
of time-domain programs and to exercise more control over the quantum circuit we
create. In our `publication in Nature `__ [[#advantage2022]_]
you can read more about how we can successfully sample GBS instances orders
of magnitudes faster than any known classical algorithm could using exact methods.
Let's have a look at the schematic below. It looks very much like the one shown
in the beginner's tutorial, only it contains more details:
- There is now a phase gate located in each of the three delay lines.
- After the three cascaded loops, the one spatial mode branches out into a
binary tree of 16 spatial modes, each one terminating with a measurement by
its own photon-number resolving detector.
.. image:: /tutorials/images/borealis_schematic_advanced.png
:width: 100%
:align: center
This advanced tutorial will direct more focus on these additional details and
help us get a more comprehensive understanding of Borealis and its gate
arguments.
"""
######################################################################
# Loading the device specification
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# We start off by importing Strawberry Fields and NumPy, and loading the
# Borealis :class:`device ` object which contains relevant and up-to-date
# information about the hardware.
#
import strawberryfields as sf
import numpy as np
eng = sf.RemoteEngine("borealis")
device = eng.device
######################################################################
# .. note::
#
# This demo is about remotely submitting jobs to our quantum-advantage device
# Borealis. Borealis is now deprecated, but you can download a dummy :download:`device
# specification <../tutorials/borealis_device_specification.json>` and :download:`hardware certificate <../tutorials/borealis_device_certificate.json>`.
# Load these files and use them to create a :class:`device ` object like what is shown
# below. This object will allow you to go through most steps of this tutorial,
# although you won't be able to run anything on the Borealis hardware
# device, nor receive any samples.
#
# .. code::
#
# import json
#
# device_spec = json.load(open("borealis_device_spec.json"))
# device_cert = json.load(open("borealis_device_certificate.json"))
# device = sf.Device(spec=device_spec, cert=device_cert)
#
######################################################################
# Defining a GBS circuit
# ~~~~~~~~~~~~~~~~~~~~~~
#
# As in the beginner's tutorial, here we will define a GBS circuit similar to the one in the
# beginner's tutorial which we can later submit to the hardware and/or to local
# simulations. This time, we want to define the circuit and all of the gate
# parameters manually, relying less on the handy convenience functions that we
# used in the beginner's tutorial.
#
# In the beginner's tutorial matching the squeezing values to the closest one
# supported by the hardware was done behind the scenes. Here, we will instead
# do it manually by using one of the allowed squeezing values from the device
# certificate and broadcasting it over a list that is as long as the number of modes
# that we wish to use.
modes = 216
r0 = device.certificate["squeezing_parameters_mean"]["high"]
r = [r0] * modes
######################################################################
# Next, we will define our GBS circuit by randomizing the phase-gate and
# beamsplitter arguments, very much like we did in the beginner's tutorial. We
# will also have to make sure that the loops are completely filled with computational
# modes at the beginning of the circuit by opening up the loops in the
# initial time bins by setting them to :math:`T=1` (:math:`\alpha=0`).
min_phi, max_phi = -np.pi / 2, np.pi / 2
min_T, max_T = 0.4, 0.6
# rotation-gate parameters
phi_0 = np.random.uniform(low=min_phi, high=max_phi, size=modes).tolist()
phi_1 = np.random.uniform(low=min_phi, high=max_phi, size=modes).tolist()
phi_2 = np.random.uniform(low=min_phi, high=max_phi, size=modes).tolist()
# beamsplitter parameters
T_0 = np.random.uniform(low=min_T, high=max_T, size=modes)
T_1 = np.random.uniform(low=min_T, high=max_T, size=modes)
T_2 = np.random.uniform(low=min_T, high=max_T, size=modes)
alpha_0 = np.arccos(np.sqrt(T_0)).tolist()
alpha_1 = np.arccos(np.sqrt(T_1)).tolist()
alpha_2 = np.arccos(np.sqrt(T_2)).tolist()
# the travel time per delay line in time bins
delay_0, delay_1, delay_2 = 1, 6, 36
# set the first beamsplitter arguments to `T=1` (`alpha=0`) to fill the
# loops with pulses
alpha_0[:delay_0] = [0.0] * delay_0
alpha_1[:delay_1] = [0.0] * delay_1
alpha_2[:delay_2] = [0.0] * delay_2
######################################################################
# Here comes a tricky part. Although these gate sequences are of length 216
# (the number of squeezed-light modes in our GBS instance) the actual
# time-domain program needs to be longer than that. Why is that the case?
#
# We just imposed that the first computational mode is fully coupled into the
# first loop. This means that it will arrive at the second loop after exactly
# one time bin (i.e. one roundtrip in the first delay line). We also imposed
# that the first mode should be completely coupled into the second loop,
# which means it will arrive at the third loop no earlier than 7 time bins after the
# program has started (1 + 6 time bins delay). Therefore, the gates associated
# to the second and third loop need some idle time bins, awaiting the first
# computational mode.
#
# To do this, we need to append 0-valued parameters (i.e. full transmission) to the
# gates in the first two loops. This will make sure that the non-trivial random
# gates are applied to the computational modes exclusively.
#
# second loop needs to await the first mode for 1 time bin (delay_0 = 1)
phi_1 = [0] * delay_0 + phi_1
alpha_1 = [0] * delay_0 + alpha_1
# second loop needs to await the first mode for 7 time bins (delay_0 + delay_1 = 1 + 6)
phi_2 = [0] * (delay_0 + delay_1) + phi_2
alpha_2 = [0] * (delay_0 + delay_1) + alpha_2
######################################################################
# Toward the end of the time-domain circuit, we need to make sure that all
# computational modes in the loops are entirely coupled out again, avoiding
# interfering with the vacuum modes that will follow the pulse train of 216
# squeezed-light modes. This means that we need to set the beamsplitters to full
# transmission again, and, since all gate sequences should have the
# same length, the gates of the first and second loops need to remain idle until
# the very last computational mode has exited the third loop.
#
# This requires us to add 0-valued parameters corresponding to the number of
# modes required to exit through the current and all following loops; i.e., 43
# (1 + 6 + 36) for the first loop, 42 (6 + 36) for the second loop, and 36 for
# the third loop.
#
# empty first loop (delay_0) and await the emptying of second and third loop
# (delay_1 + delay_2)
r += [0] * (delay_0 + delay_1 + delay_2)
phi_0 += [0] * (delay_0 + delay_1 + delay_2)
alpha_0 += [0] * (delay_0 + delay_1 + delay_2)
# empty second loop (delay_1) and await the emptying of third loop (delay_2)
phi_1 += [0] * (delay_1 + delay_2)
alpha_1 += [0] * (delay_1 + delay_2)
# empty third loop
phi_2 += [0] * delay_2
alpha_2 += [0] * delay_2
######################################################################
# Unfortunately, all of this padding with idle time bins before and after the
# random gates is necessary to avoid interfering computational modes with
# vacuum modes. This simply reflects the reality that we need a total unitary that
# fits the size of the total input state: 216 computational modes and 43 vacuum
# modes (1, 6, 36 of them in the first, second, and third loop, respectively).
# Our total unitary thus needs to be of size 259 :math:`\times` 259.
#
# The image below provides an overview of the necessary vacuum mode padding for
# all the three loops, where each row corresponds to a gate parameter array with
# length 259.
#
######################################################################
# .. image:: /tutorials/images/borealis_padding.png
# :width: 100%
# :align: center
# :target: javascript:void(0);
#
######################################################################
# The reward for taking all of these precautions is that the unitary that we
# desire to implement is now a block matrix within the total unitary, acting on the 216
# computational modes exclusively. Moreover, the length of the argument lists
# submitted to the seven quantum gates is now exactly 259, matching the size of
# the total unitary. Let's collect all of the gate sequences in a single nested list that we
# can submit to the :class:`TDMProgram `. The list, aptly named ``gate_args_list``,
# contains all gate parameters in the order in which the corresponding gates
# will be applied.
gate_args_list = [r, phi_0, alpha_0, phi_1, alpha_1, phi_2, alpha_2]
[len(gate_args) for gate_args in (gate_args_list)]
######################################################################
# .. rst-class:: sphx-glr-script-out
#
# Out:
#
# .. code-block:: none
#
# [259, 259, 259, 259, 259, 259, 259]
#
######################################################################
# As mentioned in the beginner's tutorial, each loop applies a static phase rotation
# to the modes passing through it. Currently, this phase rotation cannot be programmed. All
# we can do is to *measure* this phase and take it into account whenever we run
# a program on the hardware. This is why the three loop offsets ``"loop_phases"``
# are part of the device certificate. Even though we have no direct control over
# them, we need to deal with them in order to obtain a correct representation of
# the circuit that is applied. This can be done in two different ways:
#
# - We include these offsets explicitly into our :class:`TDMProgram `. In this case,
# no hidden compensation takes place and the modulators apply the exact
# phases defined by the user.
#
# - We can stay agnostic about the static offsets and use our phase modulators
# to actively compensate for them in the background. This is done
# automatically when the loop offsets are not found in the :class:`TDMProgram `
# during compilation. In this case, the Strawberry Fields compiler
# will absorb the three roundtrip phases into the instructions for the phase
# modulators such that they *effectively* apply the user-defined phase-gate
# arguments. There is a caveat to this approach though, since ~50% of all the
# loop offsets *cannot* be compensated for due to the limited range of our
# phase modulators. In these cases, those phases are shifted by pi.
#
# While in the beginner's tutorial we went for the second option, we will go with the
# first option this time. In other words, we will acknowledge the roundtrip phases
# as part of our circuit and will explicitly add them to our time-domain
# program. In order to do so, let's first obtain them from the certificate.
# the intrinsic roundtrip phase of each delay line
phis_loop = device.certificate["loop_phases"]
######################################################################
# Executing the circuit
# ~~~~~~~~~~~~~~~~~~~~~
#
# Submit to hardware
# <<<<<<<<<<<<<<<<<<
#
# From here on, we can proceed very much like we did in the beginner's tutorial.
# Only this time, our :class:`TDMProgram ` explicitly includes rotation-gates with the
# loop-roundtrip phase shifts applied. The results retrieved from running this program on
# the Borealis device will essentially be the same (albeit with differently randomized
# parameters) as with the program we ran in the beginner's tutorial, except for the lack of automatic
# cropping, which we will do manually instead.
from strawberryfields.tdm import get_mode_indices
from strawberryfields.ops import Sgate, Rgate, BSgate, MeasureFock
delays = [1, 6, 36]
vac_modes = sum(delays)
n, N = get_mode_indices(delays)
prog = sf.TDMProgram(N)
with prog.context(*gate_args_list) as (p, q):
Sgate(p[0]) | q[n[0]]
for i in range(len(delays)):
Rgate(p[2 * i + 1]) | q[n[i]]
BSgate(p[2 * i + 2], np.pi / 2) | (q[n[i + 1]], q[n[i]])
Rgate(phis_loop[i]) | q[n[i]]
MeasureFock() | q[0]
shots = 250_000
eng = sf.RemoteEngine("borealis")
results = eng.run(prog, shots=shots)
######################################################################
# Running the circuit remotely on Borealis will return a :class:`Result ` object,
# which contains the samples that are obtained by the hardware. The samples
# returned by any :class:`TDMProgram ` are of shape ``(shots, spatial_modes,
# temporal_modes)``, so in our case whe should get a three-dimensional array
# with shape ``(shots, 1, 259)``.
#
# Since the first 43 modes that arrived at the detectors are vacuum
# modes that live inside the loops before the program starts, we can crop the
# samples array and focus on the 216 *computational* modes. Recall from the beginner's
# tutorial that this can also be achieved by passing the ``crop=True`` flag to
# the :meth:`eng.run() ` call above.
samples = results.samples[..., 43:]
#######################################################################
# Submit to realistic simulation
# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
#
# Now, let us switch to a local simulation. In order to make it more realistic, we
# need to take optical loss into account. In the beginner's tutorial, realistic loss
# was added automatically by the compiler in the backend. Here, we instead want to insert the various
# loss components manually into our circuit. The daily calibration routine that
# runs on Borealis decomposes the total loss suffered by our pulses into three
# categories:
#
# - ``eta_glob``: The total loss experienced by the different modes depends on
# which loops they have travelled through and in which detector they ended up.
# So it will be highly inhomogeneous within the 216 computational modes.
# ``eta_glob``, on the other hand, describes their *common* efficiency. It is
# experienced by a mode that bypasses all three loops and ends up in our most
# efficient detector. In other words ``1 - eta_glob`` is the minimal loss that
# every mode will suffer. On top of that, most modes will suffer additional
# losses imposed by the loops and a less efficient detector.
# - ``eta_loop``: Each loop will add an additional loss to each mode that passes through it.
# ``eta_loop`` is a list containing the respective efficiencies of the three loops.
# - ``eta_ch_rel``: In the experiment, two consecutive pulses are only 167 ns
# apart, which exceeds the rate that one single photon-number resolving (PNR)
# detector can presently handle. So although our batches of 216 squeezed-light modes
# travel towards the detection module in one single fibre, they need to be
# demultiplexed into 16 spatial modes, each one of them probed by its own PNR
# detector. You can find a schematic of our demultiplexing-and-detection
# module illustrated in the top figure of this tutorial.
#
# We have 16 PNRs operating at 6/16 MHz. The demultiplexing pattern is
# repetitive, meaning that detector :math:`i` will detect all temporal modes
# whose index satisfies :math:`i + 16n`, where :math:`n` is an integer. Since
# the 16 spatial modes and their respective detectors do not share the exact
# same efficiency, ``eta_ch_rel`` keeps track of the loss suffered by each one
# of our modes individually.
#
# All of these loss components can be obtained from the :class:`device ` object.
eta_glob = device.certificate["common_efficiency"]
etas_loop = device.certificate["loop_efficiencies"]
etas_ch_rel = device.certificate["relative_channel_efficiencies"]
######################################################################
# Note that the detector-channel efficiencies are cyclic, so mode ``i`` sees the
# same efficiency as mode ``i + 16``, thus necessitating us to repeat the 16
# efficiency values for the length of the program. We will also add the channel
# efficiencies to our list of gate arguments that we submit to the
# :class:`TDMProgram `. This is because the relative detection efficiency experienced
# by the modes needs to be updated at each time bin and is therefore a dynamic
# TDM gate, very much like the phase modulators and beamsplitters.
prog_length = len(gate_args_list[0])
reps = int(np.ceil(prog_length / 16))
etas_ch_rel = np.tile(etas_ch_rel, reps)[:prog_length]
etas_ch_rel = etas_ch_rel.tolist()
gate_args_list += [etas_ch_rel]
######################################################################
# With that, we are now able to redefine the circuit that we submitted to
# Borealis for simulations under realistic losses.
from strawberryfields.ops import LossChannel
prog_sim = sf.TDMProgram(N)
with prog_sim.context(*gate_args_list) as (p, q):
Sgate(p[0]) | q[n[0]]
LossChannel(eta_glob) | q[n[0]]
for i in range(len(delays)):
Rgate(p[2 * i + 1]) | q[n[i]]
BSgate(p[2 * i + 2], np.pi / 2) | (q[n[i + 1]], q[n[i]])
Rgate(phis_loop[i]) | q[n[i]]
LossChannel(etas_loop[i]) | q[n[i]]
LossChannel(p[7]) | q[0]
MeasureFock() | q[0]
######################################################################
# As we have now inserted loss explicitly to the circuit, we can skip any
# Borealis-specific compiling steps and submit our program directly to the
# Gaussian backend.
eng_sim = sf.Engine(backend="gaussian")
results_sim = eng_sim.run(prog_sim, shots=None, space_unroll=True)
######################################################################
# The covariance matrix returned by :meth:`results_sim.state.cov() ` describes our 216
# computational modes as well as the 43 vacuum modes that were stored in the
# loops before executing the program. In the beginner's tutorial, we made sure
# that only the computational modes are considered in the covariance matrix by
# adding ``crop=True`` to the :meth:`run() ` method. We could have done the same
# thing here, but instead let's crop out the vacuum modes manually using the
# package *The Walrus*.
from thewalrus.quantum import reduced_gaussian
# total covariance matrix and mean vector describe 259 modes
cov_tot = results_sim.state.cov()
mu_tot = np.zeros(len(cov_tot))
# cropped covariance matrix and mean vector describe 216 modes
mu, cov = reduced_gaussian(mu_tot, cov_tot, range(vac_modes, prog_length))
######################################################################
# From here you can use the methods described in the beginner's tutorial to compare
# your experimental results to the ones obtained by the simulation.
#
# Note that in the above :meth:`eng_sim.run() ` call we added the flag
# ``space_unroll=True``. This feature was already used and very briefly
# explained in the beginner's tutorial. Let's go into more details here. The
# circuit above will be rolled-up, using symbolic placeholder variables for the
# gate parameters. The circuit must be unrolled with each mode only being
# measured once, a process we call *space-unrolling*.
#
# Space-unrolling might seem slightly confusing, but there is a straight-forward
# way to understand it. Instead of having several modes going through the same path,
# we interpret the TDM loop interactions as spatially separate. This way, each
# measurement is made on a different spatial mode, as can be seen in the
# right-hand side of the image below.
#
# The left-hand side depicts how an actual single loop TDM circuit would work,
# while the right-hand side is how Strawberry Fields simulates it. Since each
# mode can only be measured once (i.e., modes cannot be reused when using Fock
# measurements, opposite to the `Entanglement Synthesizer `__
# , which uses homodyne measurements), the circuit needs to be unrolled. To get
# an understanding for how this space-unrolling works, try to follow the
# interactions happening in the left-hand side single-loop setup, and try to
# spot them in the right-hand side circuit diagram. They are equal in practice!
#
######################################################################
# .. image:: /tutorials/images/borealis_space_unrolling.png
# :width: 100%
# :align: center
# :target: javascript:void(0);
#
######################################################################
# Calculating the transfer matrix
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# The circuit above can be described as a combination of initial state
# preparation and an interferometer over the 216 modes used in the
# circuit. This interferometer setup can be described by a so-called transfer
# matrix. It describes the quantum circuit applied to any initial state. In the
# absence of loss, the transfer matrix is equivalent to the unitary. By
# having access to the tranfer matrix, it's possible to recreate the same circuit
# using different initial states. Instead of starting with a squeezed state, as
# we did above, we could, for example, apply the transfer matrix to a coherent state or
# GKP states.
#
# There is a way to quickly calculate the transfer matrix of a circuit — as
# long as it only contains beamsplitters, rotation gates and loss channels
# — without having to execute it. This is done by using the "passive"
# compiler to compile the circuit into a single transfer gate. Thus, we
# remove the squeezing from the circuit and simply compile the program
# using the "passive" compiler, extracting the parameter from the
# resulting :meth:`PassiveChannel `.
prog_passive = sf.TDMProgram(N)
with prog_passive.context(*gate_args_list, etas_ch_rel) as (p, q):
LossChannel(eta_glob) | q[n[0]]
for i in range(len(delays)):
Rgate(p[2 * i + 1]) | q[n[i]]
BSgate(p[2 * i + 2], np.pi / 2) | (q[n[i + 1]], q[n[i]])
Rgate(phis_loop[i]) | q[n[i]]
LossChannel(etas_loop[i]) | q[n[i]]
LossChannel(p[7]) | q[0]
prog_passive.space_unroll()
prog_passive = prog_passive.compile(compiler="passive")
######################################################################
# Keep in mind that the circuit now only consists of a single :meth:`PassiveChannel `,
# which has a single matrix parameter. This parameter is the circuit's transfer
# matrix. The transfer matrix then allows us to experiment using the same setup,
# but with different state preparations.
transfer_matrix_tot = prog_passive.circuit[0].op.p[0]
# crop out the vacuum modes
transfer_matrix = transfer_matrix_tot[vac_modes:prog_length, vac_modes:prog_length]
######################################################################
# A matrix plot will give us a good intuition on the connectivity of our
# temporal modes.
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(18, 8))
ax[0].matshow(transfer_matrix.real)
ax[1].matshow(transfer_matrix.imag)
plt.show()
######################################################################
# .. image:: /tutorials/images/borealis_transfer_matrix.png
# :width: 100%
# :align: center
# :target: javascript:void(0);
#
######################################################################
# Conclusion
# ~~~~~~~~~~
#
# This tutorial covered how to submit jobs to our temporal-division multiplexing
# device, Borealis. It's rather exciting to be able to demonstrate the
# concept of quantum computational advantage by examples running on actual
# hardware. As you now have seen, there are quite a few details to keep track
# of, but as long as the circuit follows the correct layout and the parameters
# are correctly declared, anyone with Xanadu cloud access can submit jobs to be
# excecuted on the Borealis hardware, and get back samples that would be difficult, if
# not outright impossible, to retrieve using a local simulator.
#
# We hope that you've enjoyed this demonstration, and that you're as excited as
# we are about the possibilities. If you wish to dig deeper into the science
# behind all of this, we recommend that you read our paper [[#advantage2022]_].
#
######################################################################
# References
# ~~~~~~~~~~
#
# .. [#advantage2022]
#
# Madsen, L.S., Laudenbach, F., Askarani, M.F. et al.
# "Quantum computational advantage with a programmable photonic processor",
# `Nature 606, 75-81 `__, 2022.
#
# About the authors
# <<<<<<<<<<<<<<<<<
#############################################################################################
# .. bio:: Theodor Isacsson
# :photo: ../_static/Theodor.jpg
#
# Theodor Isacsson is a quantum software developer at Xanadu. He works on the photonics software stack, leading the efforts to expand and maintain packages such as Strawberry Fields and The Walrus.
#
# .. bio:: Fabian Laudenbach
# :photo: ../_static/Fabian.png
#
# Fabian Laudenbach is an experimental physicist at Xanadu. He is one of the lead authors of the publication showcasing quantum computational advantage with Borealis.