r"""
Advanced timedomain photonic circuits
======================================
*Authors: Fabian Laudenbach and Nicolas Quesada*
In a :doc:`previous tutorial `
we introduced the basics of timedomain photonic circuits
 the same technology underpinning Xanadu's `quantum computational advantage hardware `,
Borealis [[#advantage2022]_]  by translating
a simple standard Strawberry Fields program into a :class:`strawberryfields.TDMProgram`.
Note that the examples considered previously are relatively simple,
involving at most one spatial mode and two concurrent modes.
In this tutorial we will go into more depth and show how one
can take typical diagrams of experimental setups and from them write their
correct representation as timedomain programs.
.. note::
We encourage you to read the introductory :doc:`timedomain photonic circuits
tutorial ` before jumping into this one!
To achieve the goal of this tutorial we will propose a set of general rules
to translate diagrams of experimental setups into timedomain Strawberry Fields
programs. We will apply these rules to two concrete experimental
setups and then, after gaining some intuition, will justify their validity
and correctness.
Timedomain rules

Now we state the rules to translate a photonic circuit diagram into
a timedomain program:
1. Identify the number of spatial modes ``M`` in your circuit. This number
is equal to the number of wires appearing in the diagram.
2. Obtain the number of concurrent modes in each spatial mode.
This is the number of modes that are alive at any given time in each wire.
For each spatial mode the number of concurrent modes is one more than the number of
modes living in the delay lines in that wire (if there is any).
The number of concurrent modes is specified as a list ``N`` where each entry is an integer
greater than zero.
3. Once you know the number of concurrent modes ``n = sum(N)``, one
needs to assign them into the different parts of the
circuit. We start with the first spatial mode placing the first qumode
``q[0]`` at the end of the wire, where the detector is, and then we update its location
by tracing its path from right to left until we can place it on the right most timebin.
If there are more concurrent modes in the first spatial mode
because of a delay line we assign modes ``q[1]`` to ``q[delay1]`` until
every concurrent mode is accounted for in the first spatial mode. The indices of the
modes must always increase from right to left.
4. We now repeat step 3 for all the spatial modes. Assume we used ``delay1`` modes
for the first spatial mode. Then we place qumode ``q[delay1+1]`` at the end
of the second wire, slide back until the first concurrent mode we meet and
then start to assign more qumodes until all the concurrent modes of the second
spatial mode (wire) are accounted for. This we continue until all spatial modes
have been assigned.
5. Once the qumodes are assigned in a diagram read off which ones
enter into the different gates. This step is more easily visualized
once the previous steps are completed in a concrete example as we do
below.
One dimensional cluster state

As a first example we consider the experiment from Yoshikawa et al.
as described in `Figure 1 `_
of Ref. [[#yoshikawa2016]_]. Below we reproduce a simplified version of their
diagram:
.. image:: /tutorials/images/one_million.svg
:align: center
:width: 75%
:target: javascript:void(0);

Note that in this figure we removed some unnecessary details from theirs.
We also make explicit that the delay line, the loop in the figure, supports
only one time bin. Finally, we have added a rotation gate by :math:`\pi/2`
so that the squeezed states get entangled in the first beamsplitter.
This is necessary to compensate for the slightly different beamsplitter
convention we use here.
From the diagram above we can already start to apply the rules introduced
above.
We immediately conclude that there 2 spatial modes; we also easily see
that the first spatial mode supports 1 concurrent mode while the second
spatial mode supports 2 concurrent modes. This is because there is one
concurrent mode living in the delay line.
Thus we can already write
"""
N = [1, 2]
######################################################################
# This simple assignment completes steps 1 and 2.
# Our next step is to assign the labels for the first spatial
# mode. We first put ``q[0]`` at the end of the first spatial mode:
#
# .. image:: /tutorials/images/one_million_q.svg
# :align: center
# :width: 75%
# :target: javascript:void(0);
# 
#
# and then we slide them backwards until we hit a concurrent mode, represented
# by the red "pulse"
#
# .. image:: /tutorials/images/one_million_q1.svg
# :align: center
# :width: 75%
# :target: javascript:void(0);
# 
#
# This completes step 3.
# For step 4 we repeat this process obtaining
#
# .. image:: /tutorials/images/one_million_q2.svg
# :align: center
# :width: 75%
# :target: javascript:void(0);
# 
# With this final diagram in mind we can now move to step 5 and write
# the circuit.
import numpy as np
import strawberryfields as sf
sq_r = 4 # Set the value of the squeezing parameter
l = 100 # We will make a 100 temporal mode cluster state
# Measurement angles for the two spatial modes.
# Note the first half of the measurement are in the X quadrature
# and the second half is in the P quadrature
# This will be important when we process the samples below!
theta1 = [0] * int(l / 2) + [np.pi / 2] * int(l / 2)
theta2 = theta1
shots = 1
prog = sf.TDMProgram(N=[1, 2])
with prog.context(theta1, theta2) as (p, q):
sf.ops.Sgate(sq_r)  q[0]
sf.ops.Sgate(sq_r)  q[2]
sf.ops.Rgate(np.pi / 2)  q[0]
sf.ops.BSgate()  (q[0], q[2])
sf.ops.BSgate()  (q[0], q[1])
sf.ops.MeasureHomodyne(p[0])  q[0]
sf.ops.MeasureHomodyne(p[1])  q[1]
eng = sf.Engine("gaussian")
result = eng.run(prog, shots=shots)
samples = result.samples
print("The shape of our samples is", samples.shape)
######################################################################
# With these samples we can verify that we have generated an entangled
# state by computing the variances of the nullifiers of certain
# linear combinations of the quadrature operators. These nullifiers are
# given by (cf. Eq. 1 of Ref. [[#yoshikawa2016]_])
#
# .. math::
# X_k = x_{A,k}+x_{B,k}+x_{A,k+1}  x_{B,k+1}, \\
# P_k = p_{A,k}+p_{B,k}p_{A,k+1}  p_{B,k+1},
#
# where :math:`A` and :math:`B` refer to the two spatial modes and :math:`k` is an index labeling
# the different temporal modes.
#
X_A = samples[0][0][: l // 2] # X samples from first spatial mode
P_A = samples[0][0][l // 2 :] # P samples from first spatial mode
X_B = samples[0][1][: l // 2] # X samples from second spatial mode
P_B = samples[0][1][l // 2 :] # P samples from second spatial mode
######################################################################
# Having split the samples we easily calculate the nullifiers
ntot = len(X_A)  1
nX = np.array([X_A[i] + X_B[i] + X_A[i + 1]  X_B[i + 1] for i in range(ntot)])
nP = np.array([P_A[i] + P_B[i]  P_A[i + 1] + P_B[i + 1] for i in range(ntot)])
nXvar = np.var(nX)
nPvar = np.var(nP)
print("The variance of the X nullifier is", nXvar)
print("The variance of the P nullifier is", nPvar)
######################################################################
# We can compare these values with expected theoretical values
print(
"The expected value of the variance of the nullifiers is",
2 * sf.hbar * np.exp(2 * sq_r),
)
######################################################################
# Two dimensional cluster state
# 
# As a second example we consider the experiment from Asavanant et al.
# as described in `Figure 1.c `_
# of Ref. [[#asavanant2019]_]. Below we reproduce a simplified version of their
# diagram:
#
# .. image:: /tutorials/images/two_dimensional.svg
# :align: center
# :width: 75%
# :target: javascript:void(0);
# 
# As with the previous example we have simplified the diagram, added red pulses
# to make explicit how many concurrent modes are in each spatial mode
# and added two extra rotation gates to compensate for the way we define our
# beamsplitter. We have also assumed that the
# second longer loop supports 5 concurrent modes.
#
# We easily identify the number of current modes in each spatial mode.
# The first and last spatial modes have no delay lines so only have one
# concurrent mode. The second and third spatial modes have delay lines
# supporting 1 and 5 concurrent modes from which we easily deduce that
N = [1, 2, 6, 1]
######################################################################
# Now we start assigning qumodes to the concurrent modes.
# We start with the first spatial mode for which we put ``q[0]`` at the end
# of the wire and slide it back until we hit a concurrent mode to obtain
#
# .. image:: /tutorials/images/two_dimensional_1.svg
# :align: center
# :width: 75%
# :target: javascript:void(0);
# 
# Now we move to the second spatial mode. For which we assign as follows
#
# .. image:: /tutorials/images/two_dimensional_2.svg
# :align: center
# :width: 75%
# :target: javascript:void(0);
# 
# Next we look at the third spatial mode for which we need to assign
# 6 qumodes. The first one is easy
#
# .. image:: /tutorials/images/two_dimensional_3.svg
# :align: center
# :width: 75%
# :target: javascript:void(0);
# 
# Then we consecutively assign the rest of the modes to the loop
#
# .. image:: /tutorials/images/two_dimensional_4.svg
# :align: center
# :width: 75%
# :target: javascript:void(0);
# 
# And finally we assign the sixth qumode of this spatial mode
#
# .. image:: /tutorials/images/two_dimensional_5.svg
# :align: center
# :width: 75%
# :target: javascript:void(0);
# 
# We are left with the last qumode which is trivial to assign and we obtain
# a fully labeled circuit
#
# .. image:: /tutorials/images/two_dimensional_6.svg
# :align: center
# :width: 75%
# :target: javascript:void(0);
# 
# With this labeled circuit we can easily write the TDMProgram
sq_r = 5
n = 400 # number of timebins
theta_A = [0] * int(n / 2) + [np.pi / 2] * int(
n / 2
) # measurement angles for detector A
theta_B = theta_A # measurement angles for detector B
theta_C = theta_A
theta_D = theta_A
shots = 1
prog = sf.TDMProgram(N)
with prog.context(theta_A, theta_B, theta_C, theta_D) as (p, q):
sf.ops.Sgate(sq_r)  q[0]
sf.ops.Sgate(sq_r)  q[2]
sf.ops.Sgate(sq_r)  q[8]
sf.ops.Sgate(sq_r)  q[9]
sf.ops.Rgate(np.pi / 2)  q[0]
sf.ops.Rgate(np.pi / 2)  q[8]
sf.ops.BSgate()  (q[0], q[2])
sf.ops.BSgate()  (q[8], q[9])
sf.ops.BSgate()  (q[2], q[8])
sf.ops.BSgate()  (q[0], q[1])
sf.ops.BSgate()  (q[3], q[9])
sf.ops.MeasureHomodyne(p[0])  q[0]
sf.ops.MeasureHomodyne(p[1])  q[1]
sf.ops.MeasureHomodyne(p[2])  q[3]
sf.ops.MeasureHomodyne(p[3])  q[9]
eng = sf.Engine("gaussian")
result = eng.run(prog, shots=shots)
samples = result.samples
######################################################################
# To verify the correctness of the simulation we can again look at the
# variances of the nullifiers. For this twodimensional cluster
# the nullifiers are given by
#
# .. math::
# X_k = x_{A,k} + x_{B,k}  \frac{1}{\sqrt{2}}\left(x_{A,k+1}+x_{B,k+1}+x_{C,k+T}+x_{D,k+T} \right),\\
# P_k = p_{A,k} + p_{B,k} + \frac{1}{\sqrt{2}}\left(p_{A,k+1}+p_{B,k+1}+p_{C,k+T}+p_{D,k+T} \right),\\
# Y_k = x_{C,k}  x_{D,k}  \frac{1}{\sqrt{2}}\left(x_{A,k+1}+x_{B,k+1}x_{C,k+T}x_{D,k+T}\right),\\
# Q_k = p_{C,k}  p_{D,k} + \frac{1}{\sqrt{2}}\left(p_{A,k+1}+p_{B,k+1}p_{C,k+T}p_{D,k+T}\right),
#
# where we used :math:`A,B,C,D` to label the four spatial modes and the index
# :math:`k` to label temporal modes.
# The variances of the nullifiers should be given by :math:`2 \hbar e^{2 r}`
# where :math:`r` is the squeezing parameter of the input squeezed states.
# With this information we numerically compute the nullifiers
X_A = samples[0][0][: n // 2] # X samples from detector A
P_A = samples[0][0][n // 2 :] # P samples from detector A
X_B = samples[0][1][: n // 2] # X samples from detector B
P_B = samples[0][1][n // 2 :] # P samples from detector B
X_C = samples[0][2][: n // 2] # X samples from detector C
P_C = samples[0][2][n // 2 :] # P samples from detector C
X_D = samples[0][3][: n // 2] # X samples from detector D
P_D = samples[0][3][n // 2 :] # P samples from detector D
T = 5 # Delay in third spatial mode
c = np.sqrt(1 / 2)
ntot = len(X_A)  T  1
X = np.array(
[
X_A[k] + X_B[k]  c * (X_A[k + 1] + X_B[k + 1] + X_C[k + T] + X_D[k + T])
for k in range(ntot)
]
)
Y = np.array(
[
X_C[k]  X_D[k]  c * (X_A[k + 1] + X_B[k + 1]  X_C[k + T]  X_D[k + T])
for k in range(ntot)
]
)
P = np.array(
[
P_A[k] + P_B[k] + c * (P_A[k + 1] + P_B[k + 1] + P_C[k + T] + P_D[k + T])
for k in range(ntot)
]
)
Q = np.array(
[
P_C[k]  P_D[k] + c * (P_A[k + 1] + P_B[k + 1]  P_C[k + T]  P_D[k + T])
for k in range(ntot)
]
)
Xvar = np.var(X)
Yvar = np.var(Y)
Pvar = np.var(P)
Qvar = np.var(Q)
print("The variance of the X nullifier is", Xvar)
print("The variance of the Y nullifier is", Yvar)
print("The variance of the P nullifier is", Pvar)
print("The variance of the Q nullifier is", Qvar)
print(
"The expected value of the variances is", 2 * sf.hbar * np.exp(2 * sq_r),
)
######################################################################
# We obtain an excellent agreement between our numerics and the theory.
######################################################################
# Why does it work?
# 
# We conclude our tutorial by arguing why the rules we introduced earlier work.
# The first two rules should be fairly straightforward to understand:
# to simulate timedomain circuits we need to correctly account for
# how many spatial modes there are and also count correctly how many
# modes are alive in each wire at a given time.
# The third and fourth rules are a bit less obvious. Why do we insist
# that the qumodes with the smallest index are the ones closest to detectors?
# The reason is simple: it is because in our description of TDMPrograms
# we always shift the modes "downwards". Recall from the basic timedomain
# tutorial that every :class:`strawberryfields.TDMProgram` is equivalent to regular Strawberry Fields
# program where we apply a :func:`strawberryfields.shift_by` operation. This is why we always want
# to respect the chronology of the measurements. This is satisfied if and only if
# for every wire the qumodes are assigned consecutively and with indices decreasing
# from left to right.
######################################################################
# References
# 
#
# .. [#advantage2022]
#
# `Quantum computational advantage with a programmable photonic processor `__.
#
# .. [#yoshikawa2016]
#
# J. Yoshikawa, S. Yokoyama, T. Kaji, C. Sornphiphatphong, Y. Shiozawa, K. Makino and A. Furusawa.
# Generation of onemillionmode continuousvariable cluster state by unlimited timedomain multiplexing.
# APL Photonics, 2016. doi:10.1063/1.4962732 .
#
# .. [#asavanant2019]
#
# W.Asavanant, Y. Shiozawa, S. Yokoyama, B. Charoensombutamon, H. Emura, R. N. Alexander, S. Takeda, J. Yoshikawa, N. C. Menicucci, H. Yonezawa and A. Furusawa.
# Generation of timedomainmultiplexed twodimensional cluster state.
# Science, 2019. doi:10.1126/science.aay2645 .