Beam Splitter Interferometer

Hello good morning dear colleagues. I’m trying to run code in Strawberry Fields and I’m facing problems in simulating and collecting results. The problem is specifically in the line “result =, shots=1)” which, when executed, the program enters a kind of infinite loop and does not stop executing unless told to stop. I would like to understand how I extract the results and reconstruct the output states.

# Define o número de modos
n_modes = 6
alpha = 2 + 0.0j
a = np.abs(alpha)

# Cria o programa
prog = sf.Program(n_modes)
eng = sf.Engine("bosonic")

with prog.context as q:
    # Prepara os estados nos modos
        Coherent(a, np.pi) | q[0]  # |A|
        Coherent(a, np.pi) | q[1]  # |B|
        Coherent(a, np.pi) | q[2]  # |C|
        Catstate(a, 0, 0, 'real', 1e-12, 2) | q[3]  # |cat|
        Catstate(a, 0, 0, 'real', 1e-12, 2) | q[4]  # |cat|
        Catstate(a, 0, 0, 'real', 1e-12, 2) | q[5]  # |cat|

    # Circuito 1
        BSgate(np.pi / 4, 0)  | (q[0], q[1])
        BSgate(np.pi / 4, 0)  | (q[2], q[3])
        BSgate(np.pi / 4, 0)  | (q[4], q[5])
        Rgate(np.pi) | q[0]

    # Medições heterodinas nos modos
        MeasureHeterodyne() | q[0]
        MeasureHeterodyne() | q[1]
        MeasureHeterodyne() | q[2]
        MeasureHeterodyne() | q[3]
        MeasureHeterodyne() | q[4]
        MeasureHeterodyne() | q[5]

result =, shots=1)

# Imprime os resultados das medições
heterodyne_samples = result.samples
print(f"Resultados das medições heterodinas: {heterodyne_samples}")
print(f"Resultados das medições heterodinas: {result}")

If you want help with diagnosing an error, please put the full error message below:

KeyboardInterrupt                         Traceback (most recent call last)
Cell In[5], line 1
----> 1 result =, shots=1).samples[:, 0]

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\, in, program, args, compile_options, **kwargs)
    565         if c.op.measurement_deps and eng_run_options["shots"] > 1:
    566             raise NotImplementedError(
    567                 "Feed-forwarding of measurements cannot be used together with multiple shots."
    568             )
--> 570 return super()._run(
    571     program_lst, args=args, compile_options=compile_options, **eng_run_options
    572 )

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\, in BaseEngine._run(self, program, args, compile_options, **kwargs)
    303 p.bind_params(args)
    304 p.lock()
--> 306 _, self.samples, self.samples_dict = self._run_program(p, **kwargs)
    307 self.run_progs.append(p)
    309 if isinstance(p, TDMProgram) and received_rolled:

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\, in BosonicEngine._run_program(self, prog, **kwargs)
    873 def _run_program(self, prog, **kwargs):
    874     # Custom Bosonic run code
--> 875     applied, samples_dict = self.backend.run_prog(prog, **kwargs)
    876     samples, samples_dict = self._combine_and_sort_samples(samples_dict)
    877     return applied, samples, samples_dict

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\backends\bosonicbackend\, in BosonicBackend.run_prog(self, prog, **kwargs)
    152 elif not isinstance(cmd.op, non_gauss_preps):
    153     try:
    154         # try to apply it to the backend and if op is a measurement, store outcome in values
--> 155         val = cmd.op.apply(cmd.reg, self, **kwargs)
    156         if val is not None:
    157             for i, r in enumerate(cmd.reg):
    158                 # Internally also store all the measurement outcomes

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\, in Measurement.apply(self, reg, backend, **kwargs)
    322 if _shots is None:
    323     return None
--> 325 values = super().apply(reg, backend, **kwargs)
    327 # store the results in the register reference objects
    328 for v, r in zip(np.transpose(values), reg):

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\, in Operation.apply(self, reg, backend, **kwargs)
    226 temp = [rr.ind for rr in reg]
    227 # call the child class specialized _apply method
--> 228 return self._apply(temp, backend, **kwargs)

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\, in MeasureHeterodyne._apply(self, reg, backend, shots, **kwargs)
   1317 def _apply(self, reg, backend, shots=1, **kwargs):
-> 1318     return backend.measure_heterodyne(*reg, shots=shots,, **kwargs)

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\backends\bosonicbackend\, in BosonicBackend.measure_heterodyne(self, mode, shots, select, **kwargs)
    787 def measure_heterodyne(self, mode, shots=1, select=None, **kwargs):
    788     if select is None:
--> 789         res = 0.5 * self.circuit.heterodyne(mode, shots=shots)
    790         return np.array([res[:, 0] + 1j * res[:, 1]]).T
    792     res = select

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\backends\bosonicbackend\, in BosonicModes.heterodyne(self, mode, shots)
    883 r"""Performs a heterodyne measurement on a mode, simulated by a generaldyne
    884 onto a coherent state.
    891     array: heterodyne outcome
    892 """
    893 covmat = self.hbar * np.eye(2) / 2
--> 894 return self.measure_dyne(covmat, [mode], shots=shots)

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\backends\bosonicbackend\, in BosonicModes.measure_dyne(self, covmat, modes, shots)
    747 diff_sample = peak_sample - means_quad
    748 # Calculate arguments for the Gaussian functions used to calculate
    749 # the exact probability distribution at the sampled point
--> 750 exp_arg = np.einsum(
    751     "...j,...jk,...k",
    752     (diff_sample),
    753     np.linalg.inv(covs_quad),
    754     (diff_sample),
    755 )
    757 # Make a copy to calculate the exponential arguments of the
    758 # upper bounding function at the point
    759 ub_exp_arg = np.copy(exp_arg)

File C:\ProgramData\anaconda3\Lib\site-packages\numpy\core\, in einsum(out, optimize, *operands, **kwargs)
   1369     if specified_out:
   1370         kwargs['out'] = out
-> 1371     return c_einsum(*operands, **kwargs)
   1373 # Check the kwargs to avoid a more cryptic error later, without having to
   1374 # repeat default values here
   1375 valid_einsum_kwargs = ['dtype', 'order', 'casting']


Anaconda 3
strawberryfields 0.23

Hi @Antonio_Aguiar ,

It’s taking a long time for me too, I’m not sure why. Have you tried with a smaller example? Maybe 2 modes instead of 6?

Hello, thanks for responding. I’ve already done tests and with up to 5 modes it runs in seconds. From the moment I add a mode and another cat state it simply does not finish its execution.

Hi @Antonio_Aguiar , it looks like the issue is adding the cat state, not the additional mode. If you use 6 modes but coherent states on the first 5 it runs instantly. I’d say you need to change the ampl_cutoff to something higher like 1e-6 or something like this. Maybe this will be too high but at least you’ll have better chances of finishing the simulation. It seems that unfortunately cat states are very computationally expensive.

1 Like

Hello again, good morning. I did what I suggested and the simulation finished in a matter of minutes. Then I tried to run it and more than an hour passed and it didn’t finish running. So I tried restarting everything and this time it finished running after almost an hour. Do you have any ideas on how to lower the simulation time? Is there anything I can do to have predictability regarding the simulation time?

Hi @Antonio_Aguiar ,

The very random simulation times sound like a cache issue, although I’m not sure and I wouldn’t know how to solve it.

Maybe you can go to the source code of Catstate and run it line by line to understand where the bottleneck is.

I’m sorry I don’t have a better answer.

Alternatively, cat states are available in Strawberry Fields as a NumPy array via sf.utils.cat_state. This isn’t something you use as part of a program, you just run it as-is. I’m not sure if it can help you but I’m sharing it just in case.

Hello, good morning Catalina Albornoz. I’m sad to hear this, but I understand the limitations of a project in development. Now I need to redefine my route and choose the most suitable path.

Regarding simulations on real machines, where would I run circuits with cat states at the input?
I made a request to have access to the X8, does it take a long time for it to be used?

In Xanadu Cloud’s Quantun Devices, is “simulon_gaussian” a real machine or a virtual simulation environment in the cloud?

If I choose to work only with coherent states on my machine, does it make a difference between running on the bosonic and gaussian backend? Is Gaussian faster?

Hi @Antonio_Aguiar ,

You should receive a response to your access request to X8 within about 1 week. If it takes longer than that please let me know so I can check with our system admin (also please make sure to check your spam folder).

In our X8 tutorial you’ll see under device details that only a few gates are supported. Catstate initializations are not supported here.

simulon_gaussian is indeed a simulator.

The Gaussian backend is efficient in terms of memory scaling but it’s not guaranteed to always be faster. I would recommend that you test them both for your specific use case and choose the one that works better.

Hi @CatalinaAlbornoz Catalina Albornoz. I was about to run my example in X8 and the following error message came up. Do you have any idea what I may have done or not done?

2024-06-17 15:38:59,218 - INFO - Compiling program for device X8_01 using compiler Xunitary.

CircuitError Traceback (most recent call last)
Cell In[6], line 1
----> 1 results =, shots=10)

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\, in, program, compile_options, recompile, **kwargs)
672 def run(
673 self, program: Program, *, compile_options=None, recompile=False, **kwargs
674 ) → Optional[Result]:
675 “”“Runs a blocking job.
677 In the blocking mode, the engine blocks until the job is completed, failed, or
704 FailedJobError: if the remote job fails on the server side (“cancelled” or “failed”)
705 “””
→ 706 job = self.run_async(
707 program, compile_options=compile_options, recompile=recompile, **kwargs
708 )
709 try:
710 while True:
711 # TODO: needed to refresh connection; remove once xcc.Connection
712 # is able to refresh config info dynamically

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\, in RemoteEngine.run_async(self, program, compile_options, recompile, **kwargs)
800 msg = f"Compiling program for device {} using compiler {compiler_name}."
→ 802 program = program.compile(device=device, **compile_options)
804 elif recompile:
805 # recompiling program
806 if compile_options:

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\, in Program.compile(self, device, compiler, **kwargs)
723 compiler = _get_compiler(compiler)
724 target = compiler.short_name
→ 726 seq = compiler.decompose(self.circuit)
728 if kwargs.get(“warn_connected”, True):
729 DAG = pu.list_to_DAG(seq)

File ~\AppData\Roaming\Python\Python311\site-packages\strawberryfields\compilers\, in Compiler.decompose(self, seq)
290 compiled.append(cmd)
292 else:
→ 293 raise pu.CircuitError(
294 “The operation {} cannot be used with the compiler ‘{}’.”.format(
295, self.short_name
296 )
297 )
299 return compiled

CircuitError: The operation Coherent cannot be used with the compiler ‘Xunitary’.

What changes do I need to make to the circuit code in question? I’m going to replace the cat states with coherent states. Is the unitary interferometer matrix mandatory in all simulations?

Hi @Antonio_Aguiar ,

In our X8 tutorial you’ll see under device details that only a few gates are supported. Coherent initializations are not supported here. Only S2gate, BSgate, MZgate, Rgate, and Interferometer operations are supported.

Hi @CatalinaAlbornoz, I noticed the limitations of the X8 regarding the operations allowed on its hardware. I would like to know if it is possible to compile the above circuit on ports compatible with your architecture.
I would also like to know if it is possible that if I have a unitary transfer matrix on hand, I can decompose it into ports compatible with the X8 architecture. For example, the unitary matrix of the three-qubit Toffoli gate, could I decompose it according to the ports compatible with X8?

Hi @Antonio_Aguiar ,
No, in general you can’t do this. X8 is not universal meaning that you can’t decompose any unitary into the gates that run on it. This applies both to general CV gates such as Catstate and to DV gates such as the Toffoli.