Add photon number and threshold detection to the Bosonic backend

[josh146]

Currently, the Bosonic backend does not support the MeasureFock and MeasureThreshold measurement operations. It would be great to add support for these two operations to the Bosonic backend.

This would entail writing functions to calculate the different probabilities up to a given cutoff, using numpy to sample from this distribution and writing update rules for updating the post measurement state.

[JeyKelly]

Hi @josh146, I've been looking at this for the past week, getting more familiarized with the software and code. I've managed to get a partially working version of measure_fock so far, by recycling the gaussian version of it, and by using the methods of the bosonic backend. It's partially working under the condition that I have a state which requires only one Gaussian state to be expressed and I let go of extra dimension added to the get_covmat_xp() and get_mean_xp() matrices. Specifically, the output of get_covmat_xp() has an extra dimension ([n, 2, 2]) compared with the scovmatxp() method to take in account the different Gaussians needed to express the state. Mu is also in the format [n, 2] instead of [1] in the Gaussian format. This leads me to think I can simply iterate and take in account the different weights which are linked to the different Gaussian modes of the bosonic state expression. However, when debugging, I see that self.state().weights() gives me outrageous numbers for any Fock state I generate with more than one photon (I'm getting a two dimensional weight of 400, -399 for single photon), so I'm not sure if I'm misunderstanding the model, or there is some kind of normalization that needs to be added.

This leads me to ask - could you help me understand why is the mu array output from BosonicBackend.circuit.mean gives a [n, 2] array, why the mean is in the format [n, 1, 2] (as opposed to a [2, 1] sized mean in the gaussian backend), and why the weight seems so big in my tests? Also, maybe the previous question would help me understand this better, but do you think sampling each Gaussian modes of the bosonic backend and summing their weighted (normalized?) contribution renders an accurate photon number at the end of the calculation?

[nquesada]

Hi @JeyKelly --- The outrageous weights from the Fock states are to be expected. Can I suggest, for ease of testing that you try a cat state?

[nquesada]

As for your second question, regarding the shape of the means, have you looked into this tutorial: https://strawberryfields.ai/photonics/demos/run_intro_bosonic.html ? For a cat state, for which one need 4 Gaussians, one will get a vector of 4 means (2-dimensional vectors). Does that help?

[nquesada]

Also, in a first stab at this, I'd suggest you consider first implementing the threshold detectors and then once you gain more understanding go with the ones for Fock measurements. Two things you might find useful are the following: A threshold measurement has two output, click and no click. If there is no click, you just need to project the measured mode into vacuum, and you don't need to change the number of weights/covs/means. If on the other hand you get a click, you need to double the number of weights/covs/means.

[JeyKelly]

As for your second question, regarding the shape of the means, have you looked into this tutorial: https://strawberryfields.ai/photonics/demos/run_intro_bosonic.html ? For a cat state, for which one need 4 Gaussians, one will get a vector of 4 means (2-dimensional vectors). Does that help?

Hi @nquesada, I did take a look at the tutorial and I have a hard time understanding why there are two seperate entries for the mean per mode. If I understand correctly, in the cat state example, the \ket{\alpha}\bra{\alpha}, \ket{-\alpha}\bra{\alpha}, \ket{\alpha}\bra{-\alpha} and \ket{-\alpha}\bra{-\alpha} have their contributions shown in the state.means() array. Why is a 1x2 matrix used to describe one mode?

[nquesada]

I think I understand your question better! So, in the Gaussian backend the internal representation is complex with \alpha \propto x + i p. In the Bosonic backend we separate the two quadratures x and p and hence we get two dimensions instead of one. Now, there is an extra complication in that Gaussian states would have only purely real means, but because we allow for things like \ket{\alpha} \bra{-\alpha} (note the minus sign, which implies this is not a hermitian operator!) then the nx2 tensor of means can be complex. Does that make more sense now?

[JeyKelly]

I think I understand your question better! So, in the Gaussian backend the internal representation is complex with \alpha \propto x + i p. In the Bosonic backend we separate the two quadratures x and p and hence we get two dimensions instead of one. Now, there is an extra complication in that Gaussian states would have only purely real means, but because we allow for things like \ket{\alpha} \bra{-\alpha} (note the minus sign, which implies this is not a hermitian operator!) then the nx2 tensor of means can be complex. Does that make more sense now?

ooh, yes it does! Thank you

[JeyKelly]

@nquesada I have managed to apply measurements of both the Fock and threshold method giving what sounds to me like reasonable outcomes (consistent with other backends, except for Fock states that have odd weights), but I'm unsure how to update the state post-measurement. I've seen that certain measurement methods of other backends refer to their self.circuit methods instead of the backend methods that are being pointed to in the ops.py file (i.e self.circuit.measure_fock in the measure_fock of the Fock backend), and these usually update their state array. I wasn't too familiar with the way private methods are declared in python (putting the "" prefix seems to be the way to do this, if I understand), and in this specific case it seems to me that the "_state" method isn't defined in the circuit.py files. In fact, I can't find it declared anywhere, even when I try to look at the call stack to see where the function comes from.

I imagine that I could simply use a circuit class method and update self._state according to measurement outcomes (i.e if I measure no click in the threshold measurement, I have a vacuum state in my channel), however I wouldn't understand what is going on in the back. Do you think you could help understand of how the state in the channel is manipulated in the program?

Also, regarding this:

If on the other hand you get a click, you need to double the number of weights/covs/means.

I can imagine that if I get no click, I simply project the mode on vacuum, however I'm not too sure what happens when it does click. By checking the section IV C) of this : https://arxiv.org/pdf/2103.05530.pdf, it seems to me like I should be expecting the outcome density matrix of a click to be the sum of a thermal and vacuum Gaussian state (edit: for a Fock state). How would this translate to doubling the weights/covs/means like you've mentioned?

[josh146]

@JeyKelly, while waiting for @nquesada to reply, perhaps you could create a work-in-progress PR with your currently modified branch? That way, we can more closely see your working code, and you can tag us for questions 🙂

When making the PR, simply prepend the title with [WIP].

[sduquemesa]

Hi All! I'm closing this issue as all the related pull requests are closed. If this is not yet solved feel free to reopen this issue.