You're reading the documentation of the v0.6. For the latest released version, please have a look at v0.11.

Boson Sampling

We are interested to simulate a boson sample with 14 photons and 60 modes, order of size comparable to what was done in Boson Sampling with 20 Input Photons and a 60-Mode Interferometer in a \(10^{14}\)-Dimensional Hilbert Space [1]

from IPython import display
from collections import Counter
from tabulate import tabulate
from tqdm.auto import tqdm
import gzip
import pickle
import time

import sympy as sp
import random

import perceval as pcvl
import perceval.lib.symb as symb

We define all the needed values below.

Perfect Boson sampling

n = 14       #number of photons at the input
m = 60       #number of modes
N = 5000000  #number of samplings

Generating a Haar random Unitary with Perceval

Unitary_60 = pcvl.Matrix.random_unitary(m) #creates a random unitary of dimension 60

A possible linear circuit realization of such matrix would be the following.

Here we define a 2-mode unitary circuit that we can use to decompose the 60 mode unitary

mzi = (symb.BS() // (0, symb.PS(phi=pcvl.Parameter("φ_a")))
       // symb.BS() // (1, symb.PS(phi=pcvl.Parameter("φ_b"))))
pcvl.pdisplay(mzi)

Let us decompose the unitary into a Reck’s type circuit [2] - this makes a huge circuit…

Linear_Circuit_60 = pcvl.Circuit.decomposition(Unitary_60, mzi,
                                               phase_shifter_fn=symb.PS,
                                               shape="triangle")

Running Simulation

Now we choose the way to perform the simulation with Perceval. The number of photons is within what we could simulate with a Naive backend (see here), however, the output space is far too big just to enumerate and store the states - so let us go with sampling using CliffordClifford2017 backend (see here).

Sampling_Backend = pcvl.BackendFactory().get_backend("CliffordClifford2017")

Select a random input:

#one can choose which mode he/she wants at input, or we can choose it randomly
def Generating_Input(n, m, modes = None):
    "This function randomly chooses an input with n photons in m modes."
    if modes == None :
        modes = sorted(random.sample(range(m),n))
    state = "|"
    for i in range(m):
        state = state + "0"*(1 - (i in modes)) +"1"*(i in modes)+ ","*(i < m-1)
    return pcvl.BasicState(state + ">")

input_state = Generating_Input(n, m)
print("The input state: ", input_state)

The input state:  |0,0,0,0,0,1,1,1,0,1,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,1,0,0,1,1,0,0,1,0,0,0,0,0>

Just to see that it outputs a statevectors of n photon(s) in m modes.

print("The sampled outputs are:")
for _ in range(10):
    print(Sampling_Backend(Unitary_60).sample(input_state))

The sampled outputs are:
|0,0,0,0,0,1,0,2,0,1,0,0,0,0,1,0,0,0,0,1,0,1,0,1,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0>
|0,0,0,0,1,1,0,1,0,1,0,0,0,0,0,0,0,0,0,2,2,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,1>
|0,0,0,1,0,0,1,0,0,1,0,0,1,0,1,0,0,0,0,0,1,0,0,1,0,0,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,2,0,0,0,0,0,0,0>
|0,0,1,0,0,0,0,1,0,0,0,0,0,1,0,2,0,0,0,1,0,0,0,0,0,1,2,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0>
|0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,1,2,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,1,0,0,0,1,0>
|2,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,1,1>
|0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,1,3,1,1,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0>
|1,0,0,1,0,0,0,2,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,2,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,1,0,0>
|2,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,1,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1,0,0,0,1,0,0,1,0,0,0>
|0,0,0,0,0,0,1,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,1,0,2,0,0,0,0,2,0,2,0,0,0,0,0,0,0,1,0,0,0>

We carry out the sampling, we do it N times, it will take some time, let us save the results to a file:

# if we want to launch parallel process
worker_id=1

#store the input and the unitary
with open("%dphotons_%dmodes_%dsamples-worker%s-unitary.pkl" %(n,m,N,worker_id), 'wb') as f:
    pickle.dump(Unitary_60, f)

with open("%dphotons_%dmodes_%dsamples-worker%s-inputstate.pkl" %(n,m,N,worker_id), 'w') as f:
    f.write(str(input_state)+"\n")


with gzip.open("%dphotons_%dmodes_%dsamples-worker%s-samples.txt.gz" %(n,m,N,worker_id), 'wb') as f:
    start = time.time()
    for i in range(N):
        f.write((str(Sampling_Backend(Unitary_60).sample(pcvl.BasicState(input_state)))+"\n").encode());
    end = time.time()
    f.write(str("==> %d\n" % (end-start)).encode())
f.close()

A little after (4 hours on a 3.1GHz Intel) - we do have 5M samples. We launched this on 32 threads for 2 days and collected 300 batches of 5M samples

Let us analyze the K-first mode bunching on these samples

import gzip

worker_id = 1
count = 0
bunching_distribution = Counter()

with gzip.open("%dphotons_%dmodes_%dsamples-worker%s-samples.txt.gz"%(n,m,N,worker_id), "rt") as f:
    for l in f:
        l = l.strip()
        if l.startswith("|") and l.endswith(">"):
            try:
                st = pcvl.BasicState(l)
                count+=1
                bunching_distribution[st.photon2mode(st.n-1)]+=1
            except Exception:
                pass
print(count, "samples")
print("Bunching Distribution:", "\t".join([str(bunching_distribution[k]) for k in range(m)]))

These numbers have been used on 300 samples for certification - see our article on Perceval for more details.

Boson sampling with non perfect sources

Let us explore around performing Boson sampling with a non perfect source. We declare a source with 90% brightness and purity.

source = pcvl.Source(brightness=0.90, purity=0.9)
QPU = pcvl.Processor({1:source,2:source, 3:source }, Linear_Circuit_60)

We can see what is the source distribution, so how likely a state at the input of the linear circuit will be.

pcvl.pdisplay(QPU.source_distribution, precision=1e-7, max_v=20)

[1]: Hui Wang, et al. Boson Sampling with 20 Input Photons and a 60-Mode Interferometer in a \(10^{14}\)-Dimensional Hilbert Space. Physical Review Letters, 123(25):250503, December 2019. Publisher: American Physical Society.

[2]: Michael Reck, Anton Zeilinger, Herbert J Bernstein, and Philip Bertani. Experimental realization of any discrete unitary operator. Physical review letters, 73(1):58, 1994.