ciq_tutorial.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0a97c59f",
   "metadata": {},
   "source": [
    "# Introduction to Perceval and Remote Quantum Computing with the Quandela Cloud\n",
    "\n",
    "In this tutorial, we will cover some basic code to get our hands on Perceval.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "4ad79fd5",
   "metadata": {},
   "source": [
    "## I. Introduction\n",
    "\n",
    "### 1. Perceval installation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7537273d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import perceval as pcvl\n",
    "pcvl.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28abc3e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from perceval.components.unitary_components import PS, BS, PERM\n",
    "from perceval.algorithm import Sampler, Analyzer\n",
    "import numpy as np"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "8817daed",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 2. BasicStates\n",
    "\n",
    "In Linear Optical (LO) circuits, photons can have many discrete degrees of freedom, which can be used as computational modes: frequency, polarisation, position, ...\n",
    "\n",
    "We represent these degrees of freedom with Fock states. If we have $n$ photons over $m$ modes, the Fock state $|s_1,s_2,...,s_m\\rangle$ means we have $s_i$ photons in the $i^{th}$ mode, i.e. $\\sum_{i=1}^m s_i =n$.\n",
    "\n",
    "In Perceval, we will use the module `pcvl.BasicState`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "435b102b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Difference syntaxes of BasicStates (list, string, etc.)\n",
    "bs1 = pcvl.BasicState('|0,2,0,1>')\n",
    "bs2 = pcvl.BasicState([0, 2, 0, 1])\n",
    "\n",
    "if bs1==bs2:print('Those are the same states')\n",
    "\n",
    "## You can iterate on modes\n",
    "for i, photon_count in enumerate(bs2):\n",
    "    print(f'There are {photon_count} photons in mode {i}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "5fe29a50",
   "metadata": {},
   "source": [
    "### 3. LO-components\n",
    "\n",
    "LO components are the elementary blocks which will act on our Fock states.\n",
    "\n",
    "It's important to know all the possible components that can be found in Perceval and understand their effects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3dad8b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Permutation\n",
    "perm=PERM([2,0,1])\n",
    "\n",
    "print(perm.name)\n",
    "print(perm.describe())\n",
    "pcvl.pdisplay(perm.definition())\n",
    "pcvl.pdisplay(perm)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "92d6b1db",
   "metadata": {},
   "source": [
    "The output matrix from above is $\\left[\\begin{matrix}0 & 1 & 0\\\\0 & 0 & 1\\\\1 & 0 & 0\\end{matrix}\\right]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d0cf3b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Phase shifter\n",
    "ps = PS(phi=np.pi)\n",
    "\n",
    "print(ps.name)\n",
    "print(ps.describe())  \n",
    "pcvl.pdisplay(ps.definition()) \n",
    "pcvl.pdisplay(ps)  # A pdisplay call on a circuit/processor needs to be the last line of a cell"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "17b45c73",
   "metadata": {},
   "source": [
    "The output matrix from above is $\\left[\\begin{matrix}e^{i \\phi}\\end{matrix}\\right]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d71c153b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Beam splitters\n",
    "bs_rx = BS.Rx()  # By default, a beam splitter follows the Rx gate convention, so bs=BS() has the same matrix\n",
    "\n",
    "# But other conventions of beam splitter phase exist too:\n",
    "bs_h = BS.H() \n",
    "bs_ry = BS.Ry()\n",
    "\n",
    "## Check the difference in the unitary definition:\n",
    "print('BS.Rx() unitary matrix')\n",
    "pcvl.pdisplay(bs_rx.definition())\n",
    "print('BS.H() unitary matrix')\n",
    "pcvl.pdisplay(bs_h.definition())\n",
    "print('BS.Ry() unitary matrix')\n",
    "pcvl.pdisplay(bs_ry.definition())\n",
    "print('BS displays its convention as a small label:')\n",
    "pcvl.pdisplay(bs_ry)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "38a6f440",
   "metadata": {},
   "source": [
    "The output matrices from above are:\n",
    "\n",
    "$\\left[\\begin{matrix}e^{i \\left(\\phi_{tl} + \\phi_{tr}\\right)} \\cos{\\left(\\frac{\\theta}{2} \\right)} & i e^{i \\left(\\phi_{bl} + \\phi_{tr}\\right)} \\sin{\\left(\\frac{\\theta}{2} \\right)}\\\\i e^{i \\left(\\phi_{br} + \\phi_{tl}\\right)} \\sin{\\left(\\frac{\\theta}{2} \\right)} & e^{i \\left(\\phi_{bl} + \\phi_{br}\\right)} \\cos{\\left(\\frac{\\theta}{2} \\right)}\\end{matrix}\\right]$\n",
    "\n",
    "$\\left[\\begin{matrix}e^{i \\left(\\phi_{tl} + \\phi_{tr}\\right)} \\cos{\\left(\\frac{\\theta}{2} \\right)} & e^{i \\left(\\phi_{bl} + \\phi_{tr}\\right)} \\sin{\\left(\\frac{\\theta}{2} \\right)}\\\\e^{i \\left(\\phi_{br} + \\phi_{tl}\\right)} \\sin{\\left(\\frac{\\theta}{2} \\right)} & - e^{i \\left(\\phi_{bl} + \\phi_{br}\\right)} \\cos{\\left(\\frac{\\theta}{2} \\right)}\\end{matrix}\\right]$\n",
    "\n",
    "$\\left[\\begin{matrix}e^{i \\left(\\phi_{tl} + \\phi_{tr}\\right)} \\cos{\\left(\\frac{\\theta}{2} \\right)} & - e^{i \\left(\\phi_{bl} + \\phi_{tr}\\right)} \\sin{\\left(\\frac{\\theta}{2} \\right)}\\\\e^{i \\left(\\phi_{br} + \\phi_{tl}\\right)} \\sin{\\left(\\frac{\\theta}{2} \\right)} & e^{i \\left(\\phi_{bl} + \\phi_{br}\\right)} \\cos{\\left(\\frac{\\theta}{2} \\right)}\\end{matrix}\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd1e84bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can ask for the symbolic matrix value of your component with the attribute U\n",
    "my_ps = PS(phi=np.pi/8)\n",
    "pcvl.pdisplay(my_ps.U)\n",
    "# And for the numerical value with the method compute_unitary\n",
    "pcvl.pdisplay(my_ps.compute_unitary())\n",
    "print('')\n",
    "\n",
    "# If you do it for a beam splitter, you can see that by default theta=pi/2, and the phi's are 0\n",
    "print('A default beam splitter:')\n",
    "pcvl.pdisplay(BS().compute_unitary())  # this is a balanced beam splitter\n",
    "print('')\n",
    "\n",
    "# To control the value of the parameters of a component, several choices are possible:\n",
    "#  - by setting a numerical value during the creation of the component\n",
    "print('A beam splitter with a numerical value for theta:')\n",
    "bs_rx = BS.Rx(theta=10)\n",
    "pcvl.pdisplay(bs_rx.U)\n",
    "pcvl.pdisplay(bs_rx.compute_unitary())\n",
    "print('')\n",
    "\n",
    "#  - by using the syntax pcvl.P to create a symbolic variable \n",
    "#    (note that you cannot compute the numerical value of your component anymore)\n",
    "print('A phase shifter with a symbolic value for phi:')\n",
    "ps = PS(phi=pcvl.P('\\psi'))\n",
    "pcvl.pdisplay(ps.U)\n",
    "print('')\n",
    "\n",
    "#  - you can still modify the value of a symbolic variable after its creation\n",
    "#    This is not true for a numerical variable!\n",
    "print('A beam splitter with a symbolic variable...')\n",
    "bs_rx = BS(theta=pcvl.P('toto'))\n",
    "pcvl.pdisplay(bs_rx.U)\n",
    "bs_rx.assign({'toto':10})\n",
    "print('... set to a numerical value')\n",
    "pcvl.pdisplay(bs_rx.compute_unitary())\n",
    "print('')\n",
    "\n",
    "# To check which parameters can be modified, you can call the method get_parameters\n",
    "# You can also directly change the output of get_parameters to change the values of the parameters\n",
    "bs_rx = BS(theta=pcvl.P('toto'), phi_tl = pcvl.P('tata'), phi_tr = -1)\n",
    "parameters = bs_rx.get_parameters()\n",
    "parameters[0].set_value(np.pi)\n",
    "print('Modified parameters...')\n",
    "for param in parameters:\n",
    "    print('    ', param)\n",
    "print('... and successfully modified beam splitter:')\n",
    "pcvl.pdisplay(bs_rx.U)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7cb0eead",
   "metadata": {},
   "source": [
    "_The above output in markdown:_\n",
    "\n",
    "$\\left[\\begin{matrix}e^{0.392699081698724 i}\\end{matrix}\\right]$\n",
    "\n",
    "$\\left[\\begin{matrix}0.92388 + 0.382683 i\\end{matrix}\\right]$\n",
    "\n",
    "A default beam splitter:\n",
    "\n",
    "$\\left[\\begin{matrix}\\frac{\\sqrt{2}}{2} & \\frac{\\sqrt{2} i}{2}\\\\\\frac{\\sqrt{2} i}{2} & \\frac{\\sqrt{2}}{2}\\end{matrix}\\right]$\n",
    "\n",
    "A beam splitter with a numerical value for theta:\n",
    "\n",
    "$\\left[\\begin{matrix}\\cos{\\left(5 \\right)} & i \\sin{\\left(5 \\right)}\\\\i \\sin{\\left(5 \\right)} & \\cos{\\left(5 \\right)}\\end{matrix}\\right]$\n",
    "\n",
    "$\\left[\\begin{matrix}0.283662 & - 0.958924 i\\\\- 0.958924 i & 0.283662\\end{matrix}\\right]$\n",
    "\n",
    "A phase shifter with a symbolic value for phi:\n",
    "\n",
    "$\\left[\\begin{matrix}e^{i \\psi}\\end{matrix}\\right]$\n",
    "\n",
    "A beam splitter with a symbolic variable...\n",
    "\n",
    "$\\left[\\begin{matrix}\\cos{\\left(\\frac{toto}{2} \\right)} & i \\sin{\\left(\\frac{toto}{2} \\right)}\\\\i \\sin{\\left(\\frac{toto}{2} \\right)} & \\cos{\\left(\\frac{toto}{2} \\right)}\\end{matrix}\\right]$\n",
    "\n",
    "... set to a numerical value\n",
    "\n",
    "$\\left[\\begin{matrix}0.283662 & - 0.958924 i\\\\- 0.958924 i & 0.283662\\end{matrix}\\right]$\n",
    "\n",
    "Modified parameters...\n",
    "\n",
    "     Parameter(name='toto', value=3.141592653589793, min_v=0.0, max_v=12.566370614359172)\n",
    "\n",
    "     Parameter(name='tata', value=None, min_v=0.0, max_v=6.283185307179586)\n",
    "\n",
    "... and successfully modified beam splitter:\n",
    "\n",
    "$\\left[\\begin{matrix}6.12323399573677 \\cdot 10^{-17} e^{i \\left(tata + 5.28318530717959\\right)} & 1.0 i e^{5.28318530717959 i}\\\\1.0 i e^{i tata} & 6.12323399573677 \\cdot 10^{-17}\\end{matrix}\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b04b1194",
   "metadata": {},
   "outputs": [],
   "source": [
    "## to understand the conventions, you can note that a BS.Rx with the 4 phases phi (top left/right and bottom left/right) can be represented like that \n",
    "bs_rx_circuit=pcvl.Circuit(2) // (0,PS(phi=pcvl.P('phi_tl'))) // (1,PS(phi=pcvl.P('phi_bl'))) // BS(theta=pcvl.P('theta')) // (0,PS(phi=pcvl.P('phi_tr'))) // (1,PS(phi=pcvl.P('phi_br')))\n",
    "\n",
    "pcvl.pdisplay(bs_rx_circuit.U)\n",
    "\n",
    "# we can check it's the same as bs_rx.definition()\n",
    "pcvl.pdisplay(bs_rx_circuit)\n",
    "\n",
    "## For this cell, we needed the syntax to builds circuits... Good transition!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "94a55663",
   "metadata": {},
   "source": [
    "The output matrix from above is $\\left[\\begin{matrix}e^{i \\left(\\phi_{tl} + \\phi_{tr}\\right)} \\cos{\\left(\\frac{\\theta}{2} \\right)} & i e^{i \\left(\\phi_{bl} + \\phi_{tr}\\right)} \\sin{\\left(\\frac{\\theta}{2} \\right)}\\\\i e^{i \\left(\\phi_{br} + \\phi_{tl}\\right)} \\sin{\\left(\\frac{\\theta}{2} \\right)} & e^{i \\left(\\phi_{bl} + \\phi_{br}\\right)} \\cos{\\left(\\frac{\\theta}{2} \\right)}\\end{matrix}\\right]$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "bb4cd4d2",
   "metadata": {},
   "source": [
    "## II. LO-circuits\n",
    "\n",
    "From the LO-components, we can build a LO-circuit, i.e. a sequence of those components acting on our different modes.\n",
    "\n",
    "### 1. Syntax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ab8b90e",
   "metadata": {},
   "outputs": [],
   "source": [
    "circuit = pcvl.Circuit(3)  # Create a 3-mode circuit\n",
    "\n",
    "circuit.add(0, BS())\n",
    "circuit.add(0, PS(phi=np.pi/2)).add(1, PS(phi=pcvl.P('phi'))).add(1, BS())\n",
    "\n",
    "# Equivalent syntax:\n",
    "# circuit // BS() // PS(phi=np.pi/2) // (1, PS(phi=pcvl.P('phi'))) // (1, BS())\n",
    "\n",
    "pcvl.pdisplay(circuit.U)\n",
    "pcvl.pdisplay(circuit)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7f824d7c",
   "metadata": {},
   "source": [
    "The output matrix from above is $\\left[\\begin{matrix}\\frac{\\sqrt{2} e^{1.5707963267949 i}}{2} & \\frac{\\sqrt{2} i e^{1.5707963267949 i}}{2} & 0\\\\\\frac{i e^{i \\phi}}{2} & \\frac{e^{i \\phi}}{2} & \\frac{\\sqrt{2} i}{2}\\\\- \\frac{e^{i \\phi}}{2} & \\frac{i e^{i \\phi}}{2} & \\frac{\\sqrt{2}}{2}\\end{matrix}\\right]$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0efeef18",
   "metadata": {},
   "source": [
    "The syntax ``pcvl.P('phi')`` allows you to use parameters in the circuit, where you can assign a value or not. The behavior of the parameters of a circuit is similar to the case of the components.\n",
    "\n",
    "For instance, you can use :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1230c948",
   "metadata": {},
   "outputs": [],
   "source": [
    "params=circuit.get_parameters()\n",
    "print(params) # list of the parameters\n",
    "\n",
    "# the value is None, but we can change that with\n",
    "params[0].set_value(np.pi)\n",
    "pcvl.pdisplay(circuit)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "75e4440d",
   "metadata": {},
   "source": [
    "### 2. Mach-Zehnder interferometers\n",
    "\n",
    "The beamsplitter's angle $\\theta$ can also be defined as a parameter.\n",
    "\n",
    "However, as the reflexivity depends on the mirror, it's hard to have adaptibility on the angle. \n",
    "Therefore, in practice, we use a [Mach-Zehnder Interferometer](https://en.wikipedia.org/wiki/Mach%E2%80%93Zehnder_interferometer). \n",
    "\n",
    "The beamsplitter with a parameterised $\\theta$ is therefore implemented with a parameterised phase shifter $\\phi$ between two fixed beamsplitters.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8598d2e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO: build a circuit implementing an MZI\n",
    "\n",
    "## TODO: Check that the parameterised phase allows you to change the reflexivity of your MZI\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "9d8e718b",
   "metadata": {},
   "source": [
    "### 3. Universal circuits\n",
    "\n",
    "An operation on the modes of our circuit can also be expressed as a unitary.\n",
    "\n",
    "For three modes, the unitary $U=\\begin{pmatrix}\n",
    "a_{1,1} & a_{1,2} & a_{1,3}\\\\\n",
    "a_{2,1} & a_{2,2} & a_{2,3} \\\\ \n",
    "a_{3,1} & a_{3,2} & a_{3,3}\n",
    "\\end{pmatrix}$ performs the following operation on the Fock state basis:\n",
    "\n",
    "$$\\begin{array}{rcl}\n",
    "|1,0,0\\rangle &  \\mapsto&  a_{1,1}|1,0,0\\rangle + a_{1,2}|0,1,0\\rangle + a_{1,3}|0,0,1\\rangle\\\\\n",
    "|0,1,0\\rangle &  \\mapsto&  a_{2,1}|1,0,0\\rangle + a_{2,2}|0,1,0\\rangle + a_{2,3}|0,0,1\\rangle\\\\\n",
    "|0,0,1\\rangle &  \\mapsto&  a_{3,1}|1,0,0\\rangle + a_{3,2}|0,1,0\\rangle + a_{3,3}|0,0,1\\rangle\n",
    "\\end{array}$$\n",
    "\n",
    "In 1994, [Reck et al](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.73.58) showed that any $U$ on these modes can be implemented as an LO-circuit.\n",
    "\n",
    "This decomposition can easily be done in Perceval using beam splitters and phase shifters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a043bb05",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Generate a random unitary on n modes\n",
    "n=3\n",
    "U=pcvl.Matrix.random_unitary(n)\n",
    "\n",
    "## Decompose the unitary according to Reck et al\n",
    "circuit_u=pcvl.Circuit.decomposition(U,BS(theta=pcvl.P('theta'),phi_tr=pcvl.P('phi')),phase_shifter_fn=PS)\n",
    "\n",
    "pcvl.pdisplay(circuit_u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The error between the two unitaries is',np.linalg.norm(U-circuit_u.compute_unitary()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20b31e3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO: decompose your unitary with only phase shifters and balanced beamsplitters (MZIs).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d77fa5a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO: check the norm of the difference to be sure it has worked well \n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "b6990985",
   "metadata": {},
   "source": [
    "### 4. Black box\n",
    "\n",
    "To improve readibility, the circuit can be constructed in multiple steps which are then combined as black boxes. This will also help when we'll need generic operations.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf442c25",
   "metadata": {},
   "outputs": [],
   "source": [
    "pre_MZI = (pcvl.Circuit(4, name='Bell State Preparation')\n",
    "           .add(0, BS())\n",
    "           .add(2, BS())\n",
    "           .add(1, PERM([1, 0])))\n",
    "\n",
    "upper_MZI = (pcvl.Circuit(2, name='upper MZI')\n",
    "             .add(0, PS(phi=pcvl.P('phi_0')))\n",
    "             .add(0, BS())\n",
    "             .add(0, PS(phi=pcvl.P('phi_2')))\n",
    "             .add(0, BS()))\n",
    "\n",
    "lower_MZI = (pcvl.Circuit(2, name='lower MZI')\n",
    "             .add(0, PS(phi=pcvl.P('phi_1')))\n",
    "             .add(0, BS())\n",
    "             .add(0, PS(phi=pcvl.P('phi_3')))\n",
    "             .add(0, BS()))\n",
    "\n",
    "chip = (pcvl.Circuit(4)\n",
    "              .add(0, pre_MZI)\n",
    "              .add(0, upper_MZI, merge=False) # merge is False to show boxes, True shows sub-circuits\n",
    "              .add(2, lower_MZI, merge=False))\n",
    "\n",
    "pcvl.pdisplay(chip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74c28a38",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## You can still display the inside of black boxes with:\n",
    "pcvl.pdisplay(chip, recursive=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1bd8e2f0",
   "metadata": {},
   "source": [
    "## III. Simulation\n",
    "\n",
    "Up until this point, we have focused on creating circuits.\n",
    "It's time to learn how to sample from them or describe their output distribution.\n",
    "\n",
    "### 1. Computing probabilities\n",
    "\n",
    "We will take the [Hong-Ou-Mandel](https://en.wikipedia.org/wiki/Hong%E2%80%93Ou%E2%80%93Mandel_effect) experiment as an example.\n",
    "\n",
    "Making two indistinguishable photons, one in each mode, enter one balanced beamsplitter $BS=\\frac{1}{\\sqrt{2}} \\left[\\begin{matrix}1 & 1\\\\1& -1\\end{matrix}\\right]$, we expect the outcome to be bunched:\n",
    "\n",
    "$$|1,1\\rangle \\mapsto \\frac{|2,0\\rangle - |0,2\\rangle}{\\sqrt{2}}  $$\n",
    "\n",
    "We will show how to verify this in the next steps using the Naive backend to recover the full probability distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73e9230b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO: build the circuit with the convention above\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2818570a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Syntax to compute the amplitudes\n",
    "backend = pcvl.BackendFactory.get_backend('Naive') \n",
    "\n",
    "simulator = backend(circuit)\n",
    "print(simulator.probampli(pcvl.BasicState([1,1]), pcvl.BasicState([2,0]))) # note that it's the amplitude! \n",
    "print(simulator.probampli(pcvl.BasicState([1,1]), pcvl.BasicState([0,2])))\n",
    "print(simulator.prob(pcvl.BasicState([1,1]), pcvl.BasicState([2,0]))) # note that it's the probability!\n",
    "print(simulator.prob(pcvl.BasicState([1,1]), pcvl.BasicState([0,2])))\n",
    "\n",
    "## We can also use the Analyser module to compute a table of probabilities\n",
    "## The Analyzer works with a Processor. A Processor simulates a photon source plugged into a circuit.\n",
    "## The main syntax is:\n",
    "## >>> p = pcvl.Processor(backend_name, circuit, source)\n",
    "p = pcvl.Processor('Naive', BS())\n",
    "analyzer = Analyzer(p, [pcvl.BasicState([1,1])], '*')\n",
    "pcvl.pdisplay(analyzer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcf74970",
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO:  Choose a random unitary 3x3 U, and output the  probablity table given the input |1,1,0>.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd69cf36",
   "metadata": {},
   "source": [
    "### 2. Sampling\n",
    "\n",
    "Although it's crucial to compute the output distribution, it's not what we can expect from a photonic chip. Indeed, realistically, we only can obtain a single sample from the distribution each time we run the circuit. This can be done using the backend SLOS.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab11d846",
   "metadata": {},
   "outputs": [],
   "source": [
    "p = pcvl.Processor('SLOS', BS())\n",
    "p.with_input(pcvl.BasicState([1,1]))\n",
    "\n",
    "# The sampler holds 'probs', 'sample_count' and 'samples' calls. You can use the one that fits your needs!\n",
    "sampler = Sampler(p)  \n",
    "\n",
    "sample_count = sampler.sample_count(1000)\n",
    "print(sample_count['results'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecbaaa5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO: implement the code to sample from the 3x3 Unitary of earlier\n",
    "\n",
    "## Question: how many states do we have for 3 modes and 2 photons?\n",
    "## Answer: \n",
    "\n",
    "## Question: how many states do we have for m modes and n photons? \n",
    "## Answer: "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "58b02992",
   "metadata": {},
   "source": [
    "### 3. Variational algorithm\n",
    "\n",
    "In variational algorithms, the samples from a quantum circuit allow us to approximate an expectation value, which is then used to determine the value of a loss function. This loss function is chosen such that minimising it yields a solution to a given problem. By changing the values of the parameters in our quantum circuit, we can search for this minimum.\n",
    "\n",
    "We won't go into the details of variational algorithms. However, it may be useful to see how to perform an optimisation with Perceval.\n",
    "\n",
    "We will use the library [scipy.optimise](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize).\n",
    "\n",
    "The following code solves the problem of finding an LO-circuit which, given an input Fock state $|1,1,1,1\\rangle$, maximizes the probability of outputting $|4,0,0,0\\rangle$.\n",
    "\n",
    "The solution below works for an arbitrary $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eaa4cab5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize\n",
    "\n",
    "## Data\n",
    "n=4\n",
    "input=pcvl.BasicState([1]*n)\n",
    "output_to_max=pcvl.BasicState([n]+[0]*(n-1))\n",
    "backend = pcvl.BackendFactory.get_backend('Naive') \n",
    "\n",
    "## TODO: implement a generic circuit of size n with parameters. \n",
    "## Code the loss function to maximise the desired output. \n",
    "## Launch the optimisation procedure. \n",
    "## Output the probability and circuit obtained.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3a78cc87",
   "metadata": {},
   "source": [
    "## IV. Remote computation on a real QPU\n",
    "\n",
    "Perceval can also connect to Quandela's photonic quantum processor on the cloud."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "348fd89d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "from tqdm.notebook import tqdm"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "af07f3d3",
   "metadata": {},
   "source": [
    "First, define your circuit and input state using Perceval, as usual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be3d343f",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_state = pcvl.BasicState([1, 1])\n",
    "\n",
    "c = pcvl.Circuit(2)\n",
    "c.add(0, pcvl.BS())\n",
    "c.add(0, pcvl.PS(phi = np.pi/4))\n",
    "c.add(0, pcvl.BS())\n",
    "\n",
    "pcvl.pdisplay(c)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "339a932a",
   "metadata": {},
   "source": [
    "Now, create an account on [cloud.quandela.com](https://cloud.quandela.com/) and login to see which QPU and simulators are available, as well as their specifications. Once you have chosen, all you have to do is to copy the machine’s name. You can now define a RemoteProcessor using the name of the machine and your token. Before using your token, don’t forget to give it the right to be used on the machine you want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4893a8fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use your key here to let the system know who you are\n",
    "token_qcloud = 'YOUR_API_KEY'\n",
    "remote_simulator = pcvl.RemoteProcessor('sim:ascella', token_qcloud)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d14482b2",
   "metadata": {},
   "source": [
    "You can now remotely access the simulator specs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6baf028",
   "metadata": {},
   "outputs": [],
   "source": [
    "specs = remote_simulator.specs\n",
    "pcvl.pdisplay(specs['specific_circuit'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3337d69",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(specs['constraints'])\n",
    "print(specs['parameters'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6f47cda8",
   "metadata": {},
   "source": [
    "Now we have to specify which parameters we want to give to compute. For specific parameters, we have to use a special `set_parameter` function (or `set_parameters`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff239bed",
   "metadata": {},
   "outputs": [],
   "source": [
    "remote_simulator = pcvl.RemoteProcessor('sim:ascella', token_qcloud)\n",
    "remote_simulator.set_circuit(c)\n",
    "remote_simulator.with_input(input_state)\n",
    "remote_simulator.min_detected_photons_filter(1)\n",
    "\n",
    "remote_simulator.set_parameters({\n",
    "    'HOM': .9461,\n",
    "    'transmittance': .0765,\n",
    "    'g2': .0086\n",
    "})\n",
    "\n",
    "nsample = 10000"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1c2404b1",
   "metadata": {},
   "source": [
    "We can now use the `Sampler` with our `RemoteProcessor`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef641191",
   "metadata": {},
   "outputs": [],
   "source": [
    "sampler = Sampler(remote_simulator)\n",
    "\n",
    "remote_job = sampler.sample_count\n",
    "remote_job.name = 'My CiQ sampling job' # All jobs created by this sampler instance will have this custom name on the cloud\n",
    "remote_job.execute_async(nsample)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "944995c1",
   "metadata": {},
   "source": [
    "The job has now been sent to a distant computer. As it is an async computation, we can do other things locally before the results arrive. In our case, we will just wait for the end of the computation. If you go to the cloud website again, you can see the job and its completion status."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05ecc8c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "previous_prog = 0\n",
    "with tqdm(total=1, bar_format='{desc}{percentage:3.0f}%|{bar}|') as tq:\n",
    "    tq.set_description(f'Get {nsample} samples from {remote_simulator.name}')\n",
    "    while not remote_job.is_complete:\n",
    "        tq.update(remote_job.status.progress/100-previous_prog)\n",
    "        previous_prog = remote_job.status.progress/100\n",
    "        time.sleep(1)\n",
    "    tq.update(1-previous_prog)\n",
    "    tq.close()\n",
    "\n",
    "print(f'Job status = {remote_job.status()}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "39ef9118",
   "metadata": {},
   "source": [
    "Once the previous cell has stopped, the job is finished. We can now retrieve the results for some computation. Here, the computation should be relatively fast (unless the simulator is unavailable or there are many requests on it), so we can use the job object we created before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee7b2f4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = remote_job.get_results()\n",
    "print(results['results'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "596cbf5d",
   "metadata": {},
   "source": [
    "You can also run the same sampling on the QPU:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6e75f27",
   "metadata": {},
   "outputs": [],
   "source": [
    "remote_qpu = pcvl.RemoteProcessor('qpu:ascella', token_qcloud)\n",
    "remote_qpu.set_circuit(c)\n",
    "remote_qpu.with_input(input_state)\n",
    "remote_qpu.min_detected_photons_filter(1)\n",
    "\n",
    "nsample = 10000\n",
    "\n",
    "sampler_on_qpu = Sampler(remote_qpu)\n",
    "\n",
    "remote_job = sampler_on_qpu.sample_count\n",
    "remote_job.name = 'QPU CiQ sampling job'\n",
    "remote_job.execute_async(nsample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26442e6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "previous_prog = 0\n",
    "with tqdm(total=1, bar_format='{desc}{percentage:3.0f}%|{bar}|') as tq:\n",
    "    tq.set_description(f'Get {nsample} samples from {remote_qpu.name}')\n",
    "    while not remote_job.is_complete:\n",
    "        tq.update(remote_job.status.progress/100-previous_prog)\n",
    "        previous_prog = remote_job.status.progress/100\n",
    "        time.sleep(1)\n",
    "    tq.update(1-previous_prog)\n",
    "    tq.close()\n",
    "\n",
    "print(f'Job status = {remote_job.status()}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61a34809",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = remote_job.get_results()\n",
    "print(results['results'])\n",
    "## may print >10k counts if your job takes less than 1 second\n",
    "## because collection rate is 1Hz, but photon clock rate is >>1Hz"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "bdc9f6d6cea172960c76841e513a9eff3b3efd9cc47388e8d51859232db82bd4"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}