Design Circuit with StrwaBerryFileds

Photonic Software Strawberry Fields
1

Hello!
I am trying to design the following circuits:
modernchip

In this circuit: first two modes corresponds to first qubit 0_a and 1_a and the other 2 modes corresponds to 0_b and 1_b. Blue boxes are my phase shofters. orange boxes are 50:50 beam siplitter
I started to code but I am stuck:

    import strawberryfields as sf
from strawberryfields.ops import Ket, BSgate, Interferometer
import numpy as np

cutoff_dim = 5  # (1+ total number of photons)
paths = 2
modes = 2 * paths #I created my modes and totally I have 4 modes

initial_state = np.zeros([cutoff_dim] * modes, dtype=np.complex)
# The ket below corresponds to a single horizontal photon in each of the modes
initial_state[0, 1, 0, 1] = 1
# Permutation matrix
X = np.array([[0, 1], [1, 0]])
prog = sf.Program(modes)
with prog.context as q:
    Ket(initial_state) | q  # Initial state preparation
    for i in range(paths):
        BSgate() | (q[2 * i], q[2 * i + 1])  # First layer of beamsplitters
        
eng = sf.Engine("fock", backend_options={"cutoff_dim": cutoff_dim})
result = eng.run(prog)
state = result.state

However after that I stuck: I should produce the state (|0_a 1b>+|1_a 0_b>)/sqrt(2) in correspondence of the green line, upon post-selection of one photon per qubit modesfor given an input state |0_a 0_b>

This is the first code with StrawberryField so I will be glas if someone can explain me briefly :slight_smile:
I took cutoff like 5 but I am not so sure that it is ok?
And I also do not know how to write code for after first beam splitters there is a place which s coming cross and there is no beamsplitter there and after that I need to do a produce the state (|0_a 1b>+|1_a 0_b>)/sqrt(2) in correspondence of the green line, upon post-selection of one photon per qubit modesfor given an input state |0_a 0_b>…

2

Hi @nisq — I see you have been reading our tutorial on Quantum information with polarized light . Could you explain what the green box in your diagram is?

Perhaps more generally, before going into coding could you explain to us what is that you are trying to do.

Thanks,

Nicolas

3

I apologize for confusion. It is nothing. There was a logo of company there and I wanted to close it sorry… But up to this green line, we will suppose to prodcue a Bell state like (|0_a 1b>+|1_a 0_b>)/sqrt(2).

(The idea is kind of optimization. Suppose that beam splitter is not perfect, or suppose that we have some noise in this circuit. In ideal case if we give our theoretical values to these phase shifters, we should find some unitary U, however because of the noise or imperfect beam splitter, we may have different result from U.Let’s say we will U’ and so that I will look for right values for my phase shifters for noisy cases for finding exact U. But before the optimization step, I am trying to simulate this circuit)

4

Hi @nisq!

If I understand correctly, you are trying to produce a state (|0_a 1b>+|1_a 0_b>)/sqrt(2). Do you know the beamsplitter values that reproduce this state?

If not, you can try to optimize your circuit so that it prepares the correct state. To do this, check out this tutorial (although you may use a different circuit to the one given in the tutorial).

5

I know this tutorial thanks a lot :slight_smile:
BUt My problem is preparing the circuit. I still could not achieve that.
For instance Ieven do not know cutoff value is correct or not.
I pasted my code in the previous answer and after that I am stuck.
Could you help me please?
And I the values of beam splitter is 50:50 . I already know the values for producing this state…(I will check this step but first I need to write code for this circuit)

6

Ok just for circuit, I wrote this code:
import strawberryfields as sf
from strawberryfields import ops

prog = sf.Program(4)

with prog.context as q:
    ops.BSgate() | (q[0], q[1])
    ops.BSgate() | (q[2], q[3])
    ops.Rgate(0.5719) | q[0]
    ops.BSgate() | (q[0], q[1])
    ops.Rgate(-1.9782) | q[1]
    ops.BSgate() | (q[0], q[1])
    ops.Rgate(-1.9782) | q[2]
    ops.BSgate() | (q[2], q[3])
    ops.Rgate(-1.9782) | q[2]
    ops.BSgate() | (q[2], q[3])
eng = sf.Engine("fock", backend_options={"cutoff_dim": 5})

prog.print()
#prog.draw_circuit(tex_dir='./circuit_tex', write_to_file=True)

However I am not so sure that… Can you please check it? (For now I used random values for gates and beamsplitters)I am a bit confused. For instance in this tutorial https://strawberryfields.readthedocs.io/en/stable/code/api/strawberryfields.ops.Interferometer.html and in the tutorial I have an option for rectangular interferometers…
What should I do ?..
this is confusing… even if my code is correct what can I do with this code…how can I simulate noises…

7

Hi @nisq!

Good question. When using backends like the "fock" backend, the state is a tensor that records the amplitudes of photon number states. There is a countably infinite number of such states, so we must introduce a cutoff and only record amplitudes up to that amount. Ideally, we want to set this number to be high to reduce inaccuracy in the simulation, but this comes at a cost of increased simulation time. So the choice of cutoff is a compromise and there is no single correct answer. As a rule of thumb, try to get an idea of the typical number of photons in your system, and you will want to set the cutoff to be suitably larger than that. For example, if we expect say 1 or 2 photons typically, a cutoff of 3 or 4 should be good. For your circuit, it looks like you plan to have 4 (?) photons coming in and being passed through passive gates that do not generate or destroy photons. Hence, a cutoff of 5 should be great for this simulation.

The circuit that you shared looks along the right lines! Are you experiencing an error when you try to run? The next step would be to run the program using eng.run() and access the output state of the circuit. You could train that state using the tutorial.

You have manually coded up a circuit based on elementary gates like beamsplitters. The Interferometer gate can be used instead to create a fixed array of beamsplitters that realizes a general interferometer. Whether you use the pre-made interferometer or construct your own circuit depends on the problem you’re considering.

If you want to add noise, you can add noise channels such as the loss channel.

1 Like
8

Many thanks :slight_smile:

Here is my result. I have 2 disconnected component but I could not understand …

output
output1295×200 36.8 KB

I will try to use your suggestion for training and in the code I will the circuit component with my circuit component

And let’s make one point clear: In my curcuit I have 2 qubit and each qubit includes 2 modes. This leans that my modes are 4 right. I mean I al using modes in its right meaning

Ah the warning is coming from q[1] and q[2].
I mean the warning coming from the part which is just before the green line

sktr

But how can I code in this are. There is no beam splitter and nothing in this point on my chip…

9

Hi @nisq,

Yes, it looks like there is no connection between the [0, 1] and [2, 3] modes. This is not necessarily a problem, it just means that your four modes could be simulated separately as two systems of two modes.

Looking at your diagram, the [0, 1] and [2, 3] pairs are correlated due to the beamsplitter to the left of the green box.

Yes, I believe you are using a dual-rail type encoding where a qubit state is given by a photon being in one of two modes.

I’m not quite sure what you mean, but from the diagram it looks to me that there is a 50:50 beamsplitter just to the left of the green box, between modes 1 and 2.

10

But there should be a connection for producing a Bell state I guess…
And according to my code, do I really realize my circuit? Because for now I do not think so we will see cross line. This seems to me like a normal(logical gate)circuit…

And in my optimization scenario, we assume that: we care in the sense that each beam-splitter is a component that might not work as ideally as we expect. But it’s not something that we can control directly. So the parameters we’re interested in tuning are just the phi’s in this instance.
SO at the end of the circuit, our aim get a maximally violated Bell state.
And in the tutorial https://strawberryfields.ai/photonics/demos/run_state_learner.html I do not think so our aim is the same.
Or maybe I al wrong and these are the same purposes??
My aim is that:

In case of beam splitter is not perfect or because of other noises, I will need to chance the angle of my phase shifters for finding right unitary function. For instance in the ideal case I will have U(a1,a2,a3,a4) but because of the noises I may find different U for the same angles like U’(a1,a2,a3,a4) so for finding exact U, I need the change my angle like U(a1’,a2’,a3,a4) or U(a1,a2’,a3’,a4)…) so do you think the paper(which you gave me link) will work for me if I update the circuit and parameter parts? For instance in the link they have 7 different parameters, and in my case I have 6 beam splitter and 4 phase shifters so in total I will have 16 parameters?? I think I should think my whole circle as 1 layer

yeah maybe it is work actually…

But at which point I will be able to choose my right angles and at which point I will be able to deal with noises…

I am really sorry for my questions… This is really new for me…

11

Hi @nisq

I spent some time to recreate your circuit and to train it to produce a target_state. Check out the following code:

import strawberryfields as sf
import tensorflow as tf
from strawberryfields.ops import Ket, BSgate, Rgate
import numpy as np

cutoff_dim = 5
paths = 2
modes = 2 * paths

initial_state = np.zeros([cutoff_dim] * modes, dtype=np.complex)
initial_state[0, 1, 0, 1] = 1

prog = sf.Program(modes)

bs_theta_names = [f"bs_theta_{str(i) + str(j)}" for i in range(2) for j in range(3)]
bs_phi_names = [f"bs_phi_{str(i) + str(j)}" for i in range(2) for j in range(3)]
r_names = [f"r_{i}" for i in range(4)]

prog_params = prog.params(*bs_theta_names, *bs_phi_names, *r_names)

with prog.context as q:
    Ket(initial_state) | q  # Initial state preparation
    
    BSgate(prog_params[0], prog_params[6]) | (q[0], q[1])
    BSgate(prog_params[3], prog_params[9]) | (q[2], q[3])
    
    BSgate() | (q[1], q[2])
    
    Rgate(prog_params[12]) | q[0]
    Rgate(prog_params[13]) | q[2]
    
    BSgate(prog_params[1], prog_params[7]) | (q[0], q[1])
    BSgate(prog_params[4], prog_params[10]) | (q[2], q[3])
    
    Rgate(prog_params[14]) | q[0]
    Rgate(prog_params[15]) | q[2]
    
    BSgate(prog_params[2], prog_params[8]) | (q[0], q[1])
    BSgate(prog_params[5], prog_params[11]) | (q[2], q[3])

params_bs = tf.Variable(np.random.random((2, 2, 3)))
params_r = tf.Variable(np.random.random(4))

params = [params_bs, params_r]
    
target_state = np.zeros([cutoff_dim] * modes, dtype=np.complex)
target_state[1, 0, 1, 0] = 1
target_state = tf.constant(target_state, dtype=tf.complex64)

def overlap(params):
    mapping_theta = {f"bs_theta_{str(i) + str(j)}": params_bs[0, i, j] for i in range(2) for j in range(3)}
    mapping_phi = {f"bs_phi_{str(i) + str(j)}": params_bs[1, i, j] for i in range(2) for j in range(3)}
    mapping_r = {f"r_{str(i)}": params_r[i] for i in range(4)}
    mapping = {**mapping_theta, **mapping_phi, **mapping_r}

    eng = sf.Engine(backend="tf", backend_options={"cutoff_dim": cutoff_dim})
    result = eng.run(prog, args=mapping)

    ket = result.state.ket()
    return tf.abs(tf.reduce_sum(tf.math.conj(ket) * target_state) - 1)

reps = 1000
opt = tf.keras.optimizers.Adam()

print("Optimizing:")

for i in range(reps):
    opt.minimize(lambda: overlap(params), params)
    
    if i % 100 == 0:
        print(f"Reached step {i} with cost function {overlap(params)}")

mapping_theta = {f"bs_theta_{str(i) + str(j)}": params_bs[0, i, j] for i in range(2) for j in range(3)}
mapping_phi = {f"bs_phi_{str(i) + str(j)}": params_bs[1, i, j] for i in range(2) for j in range(3)}
mapping_r = {f"r_{str(i)}": params_r[i] for i in range(4)}
mapping = {**mapping_theta, **mapping_phi, **mapping_r}

eng = sf.Engine(backend="tf", backend_options={"cutoff_dim": cutoff_dim})
result = eng.run(prog, args=mapping)

ket = result.state.ket()
print(f"Amplitude of target state: {ket[1, 0, 1, 0]}")

Hopefully this will help with your goal. As far as I understand, your next step would be to update the target_state and introduce noise into your circuit using the LossChannel. You would then find the beamsplitter parameters that get you close to the target state. However, note that introducing noise will make your state mixed, so you will need to tweak the code to support the calculations for density matrices.

Good luck!

1 Like
12

oh thanks !
I will check the code and if I have question I will ask you.
Many thanks! you spent your times for me
so sorry
One quick question tf is coming tensorflow I guess right? because we have sf but sf has not keras so it should be tensorflow tf I guess?

13

Thanks @nisq! Good point, yes tf is TensorFlow. I added the following line to the above code:

import tensorflow as tf
14

First, I really thank you for your helps:

Here I have a few questions

1)Now we do not have 50:50 beam splitters right? we filled our parameters randomly. ? And for the phase shifters, angles were choosen randomly right

  1. I could not understand what target_state means. It means which state we want to have at the end of the code right?? For instance I want to have maximally violated bell test at the end so should I change this line?

target_state[1, 0, 1, 0] = 1

And why this is [1,0,1,0]?

  1. My circuit part is actually just this part

     prog_params = prog.params(*bs_theta_names, *bs_phi_names, *r_names)

    print(“prog_params”,prog_params)
    with prog.context as q:
    Ket(initial_state) | q # Initial state preparation

     BSgate(prog_params[0], prog_params[6]) | (q[0], q[1])
 BSgate(prog_params[3], prog_params[9]) | (q[2], q[3])
 
 BSgate() | (q[1], q[2])
 
 Rgate(prog_params[12]) | q[0]
 Rgate(prog_params[13]) | q[2]
 
 BSgate(prog_params[1], prog_params[7]) | (q[0], q[1])
 BSgate(prog_params[4], prog_params[10]) | (q[2], q[3])
 
 Rgate(prog_params[14]) | q[0]
 Rgate(prog_params[15]) | q[2]
 
 BSgate(prog_params[2], prog_params[8]) | (q[0], q[1])
 BSgate(prog_params[5], prog_params[11]) | (q[2], q[3])

and the rest of the code is tranining/optimization right?
In this case beam splitters are 50:50 or we just created the circuit and we will fill the parameters randomly?
4)Why are we using complex type like np.complex or tf.complex64
5)in this code scripts:

params_bs = tf.Variable(np.random.random((2, 2, 3)))
params_r = tf.Variable(np.random.random(4))

params_r for pahse shifters and since we have 4 pahse shifters we gave 4 as a parameter. Right?
But I could not understand for beam splitter why do we write (2,2,3) what is the meaning of 3
6)and according to this code, once we created circuit, then we are filling all the parameters randomly in the circuit and then we tried to achieve our target state right??? so is it possible to right bell states as a target state or at the end of the code can I have bell test measurements?
7) when i print initial state all matrices are zero. And when I print target state still all matrices are zero. I know in the code we wrote np.zeros but then we also wrote [1,01,0]=1 and is it normal to see zero matrix
And finally I could not manage to install cudo so at the beginning I get some errors/warning but code is working. Do you think that it will be a problem even if code works or it is not so important…

Thanks for the patience and your times :slight_smile:

15

Hey @nisq,

Now we do not have 50:50 beam splitters right? we filled our parameters randomly. ? And for the phase shifters, angles were choosen randomly right

Right, most of the beamsplitters have controllable parameters which we are hoping to change. These parameters are randomly initialized.

I could not understand what target_state means. It means which state we want to have at the end of the code right?? For instance I want to have maximally violated bell test at the end so should I change this line?

Yes, the target state is a tensor of shape (5, 5, 5, 5) which corresponds to the state we want the circuit to prepare. This tensor is a complex vector holding the state amplitudes. You will need to update it to train on your target state, I just set target_state[1, 0, 1, 0] = 1 as a simple example of a case where we’re just swapping the mode of the two photons (remember, initial_state[0, 1, 0, 1] = 1).

My circuit part is actually just this part

The circuit is everything contained within with prog.context as q:. The code before that is setup, and after that is optimization. Most of the beamsplitters aren’t 50:50 and are controllable, with the exception of BSgate() | (q[1], q[2]) (which I made 50:50 based on your diagram, but feel free to update).

Why are we using complex type like np.complex or tf.complex64

Ah, there it is because we are loading a numpy array. For example, np.zeros(..., dtype=np.complex). When we make TF tensors or convert NumPy arrays to TF tensors, we can specify a TF dtype, e.g.,
tf.constant(..., dtype=tf.complex64).

params_r for pahse shifters and since we have 4 pahse shifters we gave 4 as a parameter. Right?
But I could not understand for beam splitter why do we write (2,2,3) what is the meaning of 3

Yes there are 4 phase shifters based on the diagram. The (2,2,3) of params_bs could be made a bit tidier, but the first 2 is for selecting thetas vs phis. For example, mapping_theta = {f"bs_theta_{str(i) + str(j)}": params_bs[0, i, j] for i in range(2) for j in range(3)}. The second 2 is for selecting whether it is a beamsplitter on the top 2 rows or bottom 2. The final 3 is because there are overall 3 beamsplitters on each pair of rows.

and according to this code, once we created circuit, then we are filling all the parameters randomly in the circuit and then we tried to achieve our target state right??? so is it possible to right bell states as a target state or at the end of the code can I have bell test measurements?

Yes the circuit is just training to create a specified target state, so your next step can be to work out how to fill in target_state according to your objective. Note that the circuit doesn’t create or destroy photons, so your target state should have as many photons as the input state (currently 2).

when i print initial state all matrices are zero. And when I print target state still all matrices are zero. I know in the code we wrote np.zeros but then we also wrote [1,01,0]=1 and is it normal to see zero matrix

initial_state has shape (5, 5, 5, 5). There will be quite a few matrices that look all zero and only one of them with a single nonzero element. You could try using np.unique() to convince yourself that the 1 is there.

And finally I could not manage to install cudo so at the beginning I get some errors/warning but code is working. Do you think that it will be a problem even if code works or it is not so important…

It should be ok, if the model is training successfully with the code above.

Thanks!

1 Like
16

Thank you very much !

17

No problem! :100: :rocket:

1 Like