r"""
.. role:: html(raw)
   :format: html

.. _boson_tutorial:

Boson sampling and the permanent
================================

    "The extended Church-Turing thesis states that *'Any algorithmic process can be simulated
    efficiently using a probabilistic Turing machine'.* This ad-hoc modification of the
    Church-Turing thesis should leave you feeling rather queasy." --- Nielsen and Chuang
    [[1]_].

In this tutorial, we will walk through the application of the boson sampling.

Background theory
-----------------

Introduced by Aaronson and Arkhipov [[2]_], boson sampling presented a slight
deviation from the general approach in quantum computation. Rather than presenting a theoretical
model of universal quantum computation (i.e., a framework that enables quantum simulation of any
arbitrary Hamiltonian [[1]_]), boson sampling-based devices are instead an example
of an **intermediate quantum computer**, designed to experimentally implement a computation that
is thought to be intractable classically [[3]_].

Boson sampling proposes the following `quantum linear optics
<https://en.wikipedia.org/wiki/Linear_optical_quantum_computing>`_ scheme. An array of single-photon
sources is set up, designed to simultaneously emit single photon states into a multimode linear
`interferometer <https://en.wikipedia.org/wiki/Interferometry>`_; the results are then generated by
sampling from the probability of single photon measurements from the output of the linear
interferometer.

For example, consider :math:`N` single photon Fock states,
:math:`\ket{\psi}=\ket{m_1,m_2,\dots,m_N}`, composed of :math:`b=\sum_i m_i` photons, incident on an
:math:`N`-mode linear interferometer, which performs the following linear transformation of the
input mode creation and annihilations operators:

.. math::

    &\hat{a}^\dagger_{out_k} = \sum_{j=0}^N U_{kj}\hat{a}^\dagger_{in_j}\\ &\hat{a}_{out_k} =
    \sum_{j=0}^N U_{jk}^\dagger\hat{a}_{in_j}


Here, the unitary :math:`U` completely describes the interferometer. Thus, the probability of
detecting :math:`n_j` photons at the :math:`j\text{th}` output mode is given by

.. math:: \left|\braketT{n_1,n_2,\dots,n_N}{W}{\psi}\right|^2

where :math:`W` represents the action of :math:`U` on the Fock basis (:math:`W` is simply a
homomorphism of :math:`U`). The remarkable nature of the boson sampling problem to challenge the
extended Church-Turing thesis lies in the fact that [[2]_]

.. math::

    \left|\braketT{n_1,n_2,\dots,n_N}{W}{\psi}\right|^2 =
    \frac{\left|\text{Per}(U_{st})\right|^2}{m_1!m_2!\cdots m_N!n_1!n_2!\cdots n_N!}

i.e., the sampled single photon probability distribution is proportional to the **permanent** of
:math:`U_{st}`, a submatrix of the interferometer unitary, dependent upon the input and output Fock
states.

.. note::

    The permanent of a matrix, defined by

    .. math:: \text{Per}(A) = \sum_{\sigma=S_N}\prod_{i=1}^N A_{i\sigma(i)}

    where :math:`S_N` is the set of all permutations of :math:`N` elements, is calculated in a
    similar fashion to the determinant, but unlike the determinant, the signatures of the
    permutations are not taken into account - every permutation is taken as a positive quantity.

    In graph theory, the permanent calculates the number of perfect `matchings
    <https://en.wikipedia.org/wiki/Matching_(graph_theory)>`_ in a `bipartite graph
    <https://en.wikipedia.org/wiki/Bipartite_graph>`_ with adjacency matrix :math:`A`.

Whilst the determinant can be calculated efficiently on classical computers, computing the permanent
belongs to the computational complexity class of #P-Hard problems [[4]_], which are
strongly believed to be classically hard to calculate. (Surprisingly, even calculating the permanent
in an *approximate* manner is a member of #P and intractable classically).

This implies that simulating boson sampling cannot be done efficiently on a classical computer,
providing a potential challenge to the extended Church-Turing thesis, and demonstrating the power of
(non-universal) quantum computation.

Circuit construction and simulation
-----------------------------------

In quantum linear optics, the multimode linear interferometer is commonly decomposed into two-mode
beamsplitters (:class:`~strawberryfields.ops.BSgate`) and single-mode phase shifters
(:class:`~strawberryfields.ops.Rgate`)
[[5]_], allowing for a straightforward translation into a CV quantum circuit.

For example, in the case of a 4 mode interferometer, with arbitrary :math:`4\times 4` unitary
:math:`U`, the quantum photonics circuit is given by

.. image:: /tutorials/images/boson_sampling_ex.svg
    :align: center
    :width: 70%
    :target: javascript:void(0);

:html:`<br>`

In the above, the detectors perform Fock state measurements, and the parameters of the beamsplitters
and the rotation gates determines the unitary :math:`U`. Note that, in order to allow for arbitrary
linear unitaries for :math:`m` imput modes, we must have a minimum of :math:`m+1` columns in the
beamsplitter array [[6]_].

Simulating this circuit using Strawberry Fields is easy; we can simply read off the gates from left
to right, and convert it into the Blackbird circuit language.

To begin, we create the boson sampling quantum program using Strawberry Fields:
"""
import numpy as np

# set the random seed
np.random.seed(42)

# import Strawberry Fields
import strawberryfields as sf
from strawberryfields.ops import *

# initialize a 4 mode program
boson_sampling = sf.Program(4)

with boson_sampling.context as q:
    # prepare the input fock states
    Fock(1) | q[0]
    Fock(1) | q[1]
    Vac     | q[2]
    Fock(1) | q[3]

    # rotation gates
    Rgate(0.5719)  | q[0]
    Rgate(-1.9782) | q[1]
    Rgate(2.0603)  | q[2]
    Rgate(0.0644)  | q[3]

    # beamsplitter array
    BSgate(0.7804, 0.8578)  | (q[0], q[1])
    BSgate(0.06406, 0.5165) | (q[2], q[3])
    BSgate(0.473, 0.1176)   | (q[1], q[2])
    BSgate(0.563, 0.1517)   | (q[0], q[1])
    BSgate(0.1323, 0.9946)  | (q[2], q[3])
    BSgate(0.311, 0.3231)   | (q[1], q[2])
    BSgate(0.4348, 0.0798)  | (q[0], q[1])
    BSgate(0.4368, 0.6157)  | (q[2], q[3])

######################################################################
# A couple of things to note in this particular example:
#
# ..
#
# 1. We initialize the input state :math:`\ket{\psi} = \ket{1,1,0,1}`, by creating a single photon
#    Fock state in modes 0, 1 and 3. This is done via the :class:`~strawberryfields.ops.Fock`
#    operator. Mode 2 is initialized as a vacuum state using the ``Vac`` operator (a shortcut to
#    :class:`~strawberryfields.ops.Vacuum`). This is *optional* - modes are automatically created in
#    the vacuum state upon engine initialization.
#
# ..
#
# 2. Next we apply the rotation gates, :class:`~strawberryfields.ops.Rgate`, to each mode. The
#    resulting rotation in the phase space occurs in the anticlockwise direction, with angle
#    :math:`\phi`.
#
# ..
#
# 3. Finally, we apply the array of beamsplitters, using the :class:`~strawberryfields.ops.BSgate`
#    operator, with arguments ``(theta,phi)``.
#
#    * The transmission amplitude is then given by :math:`t=\cos\theta`
#    * The reflection amplitude is given by :math:`r=e^{i\phi}\sin{\theta}`
#
# ..
#
# 4. The rotation gate and beamsplitter parameters have been chosen at random, leading to an
#    arbitrary unitary :math:`U` acting on the input modes annihilation and creation operators.
#    After leaving the beamsplitter array, we will denote the output state by :math:`\ket{\psi'}`.
#
# ..
#
# 5. We are **not** performing Fock measurements at the output; this is to ensure the state is
#    preserved, so we can extract the joint Fock state probabilities after the beamsplitter array.
#
#    If we wish to simulate Fock measurements, we can additionally include
#
#    .. code-block:: python
#
#        MeasureFock() | q
#
#    after the beamsplitter array. After constructing the circuit and running the engine, the values of
#    the Fock state measurements will be available within the ``samples`` attribute of the
#    :class:`~strawberryfields.Result` object returned by the engine.
#
# Now that the program is created, we can initialize an engine, and execute the program
# on one of the built-in Strawberry Fields state simulators.
#
# As only Fock backends support the non-Gaussian initial states used in boson sampling,
# we will use the NumPy ``'fock'`` backend, along with a Fock state
# truncation/cutoff dimension of 7 (i.e., all information of quantum amplitudes of Fock states
# :math:`\ket{n}`, :math:`n\geq 7`, is discarded).

eng = sf.Engine(backend="fock", backend_options={"cutoff_dim": 7})

######################################################################
# We can now execute the program with the engine:

results = eng.run(boson_sampling)

######################################################################
# The state method :meth:`~strawberryfields.backends.BaseFockState.all_fock_probs` returns a
# :math:`\overbrace{D\times D\times\cdots\times D}^{\text{num modes}}` array, where :math:`D` is the
# Fock basis cutoff truncation, containing the marginal Fock state probabilities - in this example,
# it would have dimension :math:`7\times 7\times 7\times 7`. We can then access the probability of
# measuring a particular output state by indexing. For example,
#
# .. math:: \texttt{probs[1, 1, 1, 0]} = \left|\braketD{1,1,1,0}{\psi'}\right|^2

# extract the joint Fock probabilities
probs = results.state.all_fock_probs()

# print the joint Fock state probabilities
print(probs[1, 1, 0, 1])
print(probs[2, 0, 0, 1])


######################################################################
# Calculating the unitary
# -----------------------
#
# Let's explore these results further. To do so, however, we'll need to compute
# the unitary that is being applied by our beamsplitters and rotation gates.
#
# By hand
# ~~~~~~~
#
# First, import some useful libraries, such as NumPy, as well
# as the ``multi_dot`` and ``block_diag`` functions from NumPy and SciPy respectively:

import numpy as np
from numpy.linalg import multi_dot
from scipy.linalg import block_diag

######################################################################
# To start with, let's generate the matrix representing the unitary transformation of the input mode
# annihilation and creation operators. The rotation gates simply act as follows,
#
# .. math:: R(\phi)\hat{a} = \hat{a}e^{i\phi},
#
# and thus the column of rotation gates has the following block diagonal form:
#
# .. math::
#     U_\phi = \left[\begin{matrix}
#             e^{i\phi_1} & 0 & \cdots\\
#             0 & e^{i\phi_2}  & \cdots \\
#             \vdots & \vdots & \ddots \\
#         \end{matrix}\right]
#
# Generating this in NumPy:

Uphase = np.diag([np.exp(0.5719*1j),np.exp(-1.9782*1j),np.exp(2.0603*1j),np.exp(0.0644*1j)])

######################################################################
# A single beamsplitter, acting on two input modes :math:`(\hat{a}_1,\hat{a}_2)`, instead acts as
# follows:
#
# .. math::
#
#     BS(\theta, \phi) (\hat{a}_1,\hat{a}_2) = \left[\begin{matrix}
#             t & -r^*\\
#             r & t\\
#         \end{matrix}\right] \left[\begin{matrix}
#             \hat{a}_1\\
#             \hat{a}_2
#         \end{matrix}\right]
#
# where :math:`t=\cos(\theta)` and :math:`r=e^{i\phi}\sin(\theta)`. Again, like the rotation gate,
# they combine as block diagonal matrices.
#
# First of all, we need to convert the ``BSgate`` arguments, ``theta`` and ``phi`` (reproduced below
# for convenience),

BSargs = [
    (0.7804, 0.8578),
    (0.06406, 0.5165),
    (0.473, 0.1176),
    (0.563, 0.1517),
    (0.1323, 0.9946),
    (0.311, 0.3231),
    (0.4348, 0.0798),
    (0.4368, 0.6157)
]

######################################################################
# into transmission and reflection amplitudes :math:`t` and :math:`r`:

t_r_amplitudes = [(np.cos(q), np.exp(p*1j)*np.sin(q)) for q,p in BSargs]

######################################################################
# Next, we generate the individual beamsplitter unitary transformations,

BSunitaries = [np.array([[t, -np.conj(r)], [r, t]]) for t,r in t_r_amplitudes]

######################################################################
# before using the SciPy function ``scipy.linalg.block_diag`` to calculate the resulting
# :math:`U_{BS_i}`, i.e., the unitary corresponding to each column of 'independent' beamsplitters in
# the above circuit diagram:

UBS1 = block_diag(*BSunitaries[0:2])
UBS2 = block_diag([[1]], BSunitaries[2], [[1]])
UBS3 = block_diag(*BSunitaries[3:5])
UBS4 = block_diag([[1]], BSunitaries[5], [[1]])
UBS5 = block_diag(*BSunitaries[6:8])

######################################################################
# Finally, we combine the unitaries to form a single :math:`4\times 4` unitary via matrix
# multiplication; :math:`U = U_{BS_5}U_{BS_4}U_{BS_3}U_{BS_2}U_{BS_1}U_{\phi}`. Since
# ``numpy.dot`` only supports matrix multiplication of two arrays, we instead use
# ``numpy.linalg.multi_dot``:

U = multi_dot([UBS5, UBS4, UBS3, UBS2, UBS1, Uphase])
print(np.round(U, 4))

######################################################################
# We find that
#
# .. math::
#     U = \left[\begin{matrix}
#         0.2195-0.2565i & 0.6111+0.5242i & -0.1027+0.4745i & -0.0273+0.0373i\\
#         0.4513+0.6026i & 0.4570+0.0123i & 0.1316-0.4504i & 0.0353-0.0532i\\
#         0.0387+0.4927i & -0.0192-0.3218i & -0.2408+0.5244i & -0.4584+0.3296i\\
#         -0.1566+0.2246i & 0.1100-0.1638i & -0.4212+0.1836i & 0.8188+0.068i
#     \end{matrix}\right]
#
# The Gaussian compiler
# ~~~~~~~~~~~~~~~~~~~~~
#
# While we have done it by hand, Strawberry Fields also supports a Gaussian unitary
# compiler, that allows us to compile our program into a single Gaussian unitary.
#
# Note that we must create a new program without the initial Fock states,
# as the Gaussian unitary compiler only works with *Gaussian operations*.
# Compiling the program:

prog_unitary = sf.Program(4)
prog_unitary.circuit = boson_sampling.circuit[4:]
prog_compiled = prog_unitary.compile(compiler="gaussian_unitary")


######################################################################
# Printing ``prog_compiled``, we see it now consists of a single
# :class:`~strawberryfields.ops.GaussianTransform` operation, consisting of a single
# symplectic matrix:

prog_compiled.print()

######################################################################
# We can easily extract this symplectic matrix, and rewrite it as a unitary
# matrix:

S = prog_compiled.circuit[0].op.p[0]
U = S[:4, :4] + 1j*S[4:, :4]
print(U)


######################################################################
# which agrees with the result above.
#
# The Interferometer operation
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# Strawberry Fields supports the :class:`~strawberryfields.ops.Interferometer` operation, which
# allows numeric unitary matrices to be directly embedded in programs and decomposed into
# beamsplitters and rotation gates:

boson_sampling = sf.Program(4)

with boson_sampling.context as q:
    # prepare the input fock states
    Fock(1) | q[0]
    Fock(1) | q[1]
    Vac     | q[2]
    Fock(1) | q[3]

    Interferometer(U) | q

######################################################################
# Compiling this for the Fock backend, and printing the result:

boson_sampling.compile(compiler="fock").print()

######################################################################
# Comparing to the permanent
# --------------------------
#
# Now that we have the interferometer unitary transformation :math:`U`, as well as the 'experimental'
# results, let's compare the two, and see if the boson sampling result,
#
# .. math::
#
#     \left|\braketD{n_1,n_2,\dots,n_N}{\psi'}\right|^2 =
#     \frac{\left|\text{Perm}(U_{st})\right|^2}{m_1!m_2!\cdots m_N!n_1!n_2!\cdots n_N!},
#
# holds.
#
# For this example, we'll consider the output state :math:`\ket{2,0,0,1}`. Extracting
# :math:`\left|\braketT{2,0,0,1}{\psi'}\right|^2` from the output data, we see that

print(probs[2,0,0,1])

######################################################################
# Before we can calculate the right-hand-side of the equation, we need a method of calculating the
# permanent. Since the permanent is classically hard to compute, it is not provided in either NumPy
# *or* SciPy, so we will use `The Walrus <https://the-walrus.readthedocs.io>`_ library, installed
# alongside Strawberry Fields:

from thewalrus import perm

######################################################################
# Finally, how do we determine the submatrix :math:`U_{st}`? This isn't too hard [[7]_];
# first, we calculate the submatrix :math:`U_s` by taking :math:`m_k` copies of the :math:`k\text{th}`
# **columns** of :math:`U`, where :math:`m_k` is the photon number of the :math:`k\text{th}` input
# state.
#
# Since our input state is :math:`\ket{\psi}=\ket{1,1,0,1}`, we simply take single copies of the
# first, second, and fourth columns:

U[:,[0,1,3]]

######################################################################
# Next, we take :math:`n_k` copies of the :math:`k\text{th}` **row** of :math:`U_s`, where
# :math:`n_k` is the photon number of the :math:`k\text{th}` **output** state that is measured.
# Here, our measurement is :math:`\ket{2,0,0,1}\bra{2,0,0,1}` so we take two copies of the first
# row, and one copy of the last row:

U[:,[0,1,3]][[0,0,3]]

######################################################################
# Now, using the permanent function we imported above, we can take the absolute value squared of the
# permanent. Finally, we divide by the product of the input and output state number factorials. Since
# :math:`0!=1!=1`, we only need to take into account the case :math:`2!=2`:

print(np.abs(perm(U[:, [0,1,3]][[0,0,3]]))**2 / 2)

######################################################################
# Comparing this to the result from Strawberry Fields, we can see that they only differ at the **17th
# decimal point**. Calculating the exact percentage difference,

BS = np.abs(perm(U[:, [0,1,3]][[0,0,3]]))**2 / 2
SF = probs[2,0,0,1]

print(100*np.abs(BS-SF)/BS)

######################################################################
# They agree with almost negligable error! This is due to the high accuracy of the Fock backend,
# despite the Fock state truncation/cutoff.
#
# This close result stands for any other output Fock state measurement that preserves the photon
# number, for example:
#
# :math:`\ket{3,0,0,0}\bra{3,0,0,0}`:

print(probs[3,0,0,0])
print(np.abs(perm(U[:, [0,1,3]][[0,0,0]]))**2 / 6)

######################################################################
# :math:`\ket{1,1,0,1}\bra{1,1,0,1}`:

print(probs[1,1,0,1])
print(np.abs(perm(U[:, [0,1,3]][[0,1,3]]))**2 / 1)

######################################################################
# .. note::
#
#     Although returning only an **approximation** of the Fock state joint probability, and thus
#     only approximating the various submatrix permanents, the Fock backends are still computing a
#     classically intractable problem.
#
#     This is because approximating the matrix permanent remains a countably hard classical problem.
#
# References
# ----------
#
# .. [1] M.A. Nielsen and I.L. Chuang. Quantum Computation and Quantum Information. Cambridge
#        University Press, 2010. ISBN 9780511992773. URL:
#        https://books.google.ca/books?id=JRz3jgEACAAJ.
#
# .. [2] Scott Aaronson and Alex Arkhipov. The computational complexity of linear optics. Theory of
#        Computing, 9(1):143–252, 2013. doi:10.4086/toc.2013.v009a004.
#
# .. [3] Max Tillmann, Borivoje Dakić, René Heilmann, Stefan Nolte, Alexander Szameit, and Philip
#        Walther. Experimental boson sampling. Nature Photonics, 7(7):540–544, May 2013.
#        doi:10.1038/nphoton.2013.102.
#
# .. [4] L.G. Valiant. The complexity of computing the permanent. Theoretical Computer Science,
#        8(2):189–201, 1979. doi:10.1016/0304-3975(79)90044-6.
#
# .. [5] Michael Reck, Anton Zeilinger, Herbert J. Bernstein, and Philip Bertani. Experimental
#        realization of any discrete unitary operator. Physical Review Letters, 73(1):58–61, Jul 1994.
#        doi:10.1103/physrevlett.73.58.
#
# .. [6] William R Clements, Peter C Humphreys, Benjamin J Metcalf, W Steven Kolthammer, and Ian A
#        Walsmley. Optimal design for universal multiport interferometers. Optica, 3(12):1460–1465, 2016.
#        doi:10.1364/OPTICA.3.001460.
#
# .. [7] Max Tillmann, Borivoje Dakić, René Heilmann, Stefan Nolte, Alexander Szameit, and Philip Walther.
#        Experimental boson sampling. Nature Photonics, 7(7):540–544, May 2013.
#        doi:10.1038/nphoton.2013.102.
