Sf.plot_quad producing different plots for the same circuit

I am using sf.plot_quad for the same circuit. However, I get different y-axis values depending on the range of x-axis values over which I produce the plots. Below is a short code snippet that should hopefully allow you to easily reproduce this on your end. The circuit is simple enough (3 Gaussian gates on a single qumode) that the correct output state can be computed by hand. In case this is useful, I have done this calculation and I can confirm that the correct output state is being computed. In other words, results.state.means() is returning the correct values.

#
# code snippet begins
#

import strawberryfields as sf
from strawberryfields import ops
import numpy as np

# set hbar
sf.hbar = 1
# initialize simulation engine
eng = sf.Engine('gaussian')
# initialize qumode(s)
prog = sf.Program(1)
with prog.context as q:
    # create circuit
    ops.Vac | (q[0])
    ops.Pgate(-6.419) | (q[0])
    ops.Zgate(12.84) | (q[0])
    ops.Rgate(-0.0866) | (q[0])
results = eng.run(prog)
### print <x> and <p>
means = results.state.means()
print("<x>,<p>={0}".format(means))

# plot 1
# to the naked eye, this plot appears correct
# <x> and <p> both appear correct here
xvec = np.arange(0, 25, 0.01)
pvec = np.arange(0, 25, 0.1)
sf.plot_quad(state=results.state, modes=[0],
             xvec=xvec, pvec=pvec, renderer="notebook")

# plot 2
# this plot appears incorrect
# <x> appears incorrect here
# also curves appear not to integrate to 1
xvec = np.arange(1, 2.5, 0.01)
pvec = np.arange(1, 2.5, 0.01)
sf.plot_quad(state=results.state, modes=[0],
             xvec=xvec, pvec=pvec, renderer="notebook")

Hi @Juan_Adame! Thanks for your question.

Here’s the plots I generated from your code, and I can see what you mean. Not only are the peaks scaled differently, they’re in completely different places!

I’ll continue to investigate and update you when I find something more.

I think I’ve figured it out.

The values for the x and p quadratures are created by integrating over the Wigner function. To get the x quadrature, you have to integrate over the provided momentum values.

But the mean of the momentum is at 12. If you cut the integration over p off at only 2.5, you are missing the bulk of the information.

You can cut down on the size of the x values to just around the peak, and cut down on the resolution in the momentum direction, but you need to keep the momentum values covering the bulk of its behaviour.

Hope that helps :slight_smile:

Hi Christina. Thanks for your quick reply!

Let me try to rephrase your findings and you can tell me if my interpretation is correct or not.

When I use sf.plot_quad, I should make sure to choose xvec and pvec such that they cover the bulk of both the x and p quadratures? And I suppose a good place to start is to provide a lot of values centered around their mean(s)?

Exactly. The pvec in the second case was not covering the bulk of the p quadrature.

Alongside the mean, you can use the covariance to estimate how widely to space the data points.

>>> results.state.cov()
array([[ 0.10100814, -1.38627303],
       [-1.38627303, 21.50077236]])

This gives us a standard deviation (square root of covariance) for momentum of about 4.63 . You’d probably want to double that to figure out how much higher and lower than the mean you want to go to start with.