Random failure with `measure_fock` method

[say4n]

Hi folks, while working with a very simple quantum circuit with the fock backend, I came across an issue where the script crashes and this stacktrace shows up sporadically.

  File "/Users/Sayan/.pyenv/versions/3.9.7/lib/python3.9/site-packages/strawberryfields/ops.py", line 1173, in _apply
    samples = backend.measure_fock(reg, shots=shots, select=self.select, **kwargs)
  File "/Users/Sayan/.pyenv/versions/3.9.7/lib/python3.9/site-packages/strawberryfields/backends/fockbackend/backend.py", line 295, in measure_fock
    return self.circuit.measure_fock(self._remap_modes(modes), select=select)
  File "/Users/Sayan/.pyenv/versions/3.9.7/lib/python3.9/site-packages/strawberryfields/backends/fockbackend/circuit.py", line 683, in measure_fock
    i = np.random.choice(list(range(len(dist))), p=dist / sum(dist))
  File "mtrand.pyx", line 928, in numpy.random.mtrand.RandomState.choice
ValueError: probabilities contain NaN

The quantum circuit that I am simulating is

with program.context as q:
  # Prepare the coherent states for the layer. Appropriately scales
  # the amplitudes for each of the layers.
  for m in range(config.NUM_MODES):
      sf.ops.Coherent(amplitudes[layer_number] * params[f"input_codeword_arg_{m}"]) | q[m]
  
  # Displace each of the modes by using the displacement magnitudes
  # generated by the ML backend.
  for m in range(config.NUM_MODES):
      sf.ops.Dgate(params[f"displacement_magnitudes_for_each_mode_arg_{m}"]) | q[m]
  
  # Perform measurements.
  sf.ops.MeasureFock() | q

Any pointers as to what might be the issue here?

Other relevant details

Platform: macOS Monetery 12.1 Python Version: 3.9.7 Strawberry Fields Version: 0.21.0

[say4n]

Folks I did a little further digging and it turns out that this crash is being caused by the vector dist to constitute of all zeros. This crash makes sense because we normalize the probabilities in the next lines if the sum is not 1 by dividing the values by the sum of the values.

But, the question now is why is the probability distribution a zero vector?

> /Users/Sayan/.pyenv/versions/3.9.7/lib/python3.9/site-packages/strawberryfields/backends/fockbackend/circuit.py(684)measure_fock()
-> if sum(dist) != 1:
(Pdb) p dist
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
(Pdb)

Could this be a side-effect of the fixes introduced as a response to #354 (and #364)?

[josh146]

Thanks for digging into this a bit further @say4n! When you are debugging, is the line where the crash is occurring at the bugfix?

https://github.com/XanaduAI/strawberryfields/blob/4de1d4676660a55fe2fef296b365af00d62b6e8e/strawberryfields/backends/fockbackend/circuit.py#L767-L770

[say4n]

@josh146 happy to help!

The crash occurs in these lines: https://github.com/XanaduAI/strawberryfields/blob/4de1d4676660a55fe2fef296b365af00d62b6e8e/strawberryfields/backends/fockbackend/circuit.py#L682-L689

Particulary, it crashes because of p=dist / sum(dist) in L#687. Maybe L#682 is making the values all 0?

[say4n]

Hi folks, any updates on triaging the issue?

[thisac]

Hi @say4n. Could you provide a minimal non-working example of your code? Specifically, I'd like to know which parameter values that you're using in the circuit, how many modes you are using (or does it fail no matter how many modes the circuit has?), the cutoff-value that you're passing to the Fock backend, and any other backend or run options that you might have. This could help me replicating the issue and solve it.

[say4n]

Hi @thisac I haven't had the time to prepare a MWE for the code yet, but I found that replacing the line https://github.com/XanaduAI/strawberryfields/blob/4de1d4676660a55fe2fef296b365af00d62b6e8e/strawberryfields/backends/fockbackend/circuit.py#L681-L682

to the following seems to have solved my crash

- dist = dist * ~np.isclose(dist, 0.0)
+ dist = np.abs(dist)

Once, I have a little more time at hand I will craft the MWE to reproduce the crash and post it here.

[thisac]

Thank you @say4n. I can see why that could solve the issue with tiny negative values, although I'm still a bit curious why your dist, i.e., the reduced state, only consists of values very close to zero in the first place. Would you expect this to be the case?

Once, I have a little more time at hand I will craft the MWE to reproduce the crash and post it here.

That would be fantastic. Thank you!