Why Release Your Code Part II: Reproducing Scientific Results
It has been a while since we have discussed the topic of publishing your source code, and the issue of reproducibility in numerical simulations. I thought it be insightful to consider an example where we try to reproduce published results using QuTiP. Here we will try to reproduce Figure #1 from Qian et. al., "Quantum Signatures of the Optomechanical Instability", Phys. Rev. Lett. 109, 253601 (2012).
This is an ideal example for several reasons: I) It is an interesting topic. II) This paper relies almost entirely on numerical results to support their conclusion. III) It is a perfect project for applying QuTiP.
In essence, the goal in this work is to find steady state solutions for an optomechanical system where the Wigner function of the mechanical oscillator, after performing a partial trace, has a sufficient amount of negative area, such that the state can be considered non-classical. Therefore, we will be making extensive use of the steadystate solver in QuTiP. Before we discuss how to do this, let us just present the final results:
First of all yes, it does not look as good, but hey, I am also not submitting it to PRL. Second, you will notice that the subfigure on the far right looks a bit different from the results in the original PRL figure. In particular, there is an extra region of stability near Delta/wm=1 for lower coupling strengths, and the graph is quite "blocky" compared to the smooth regions presented by the original authors. Why this is the case will be the focus of the discussion that follows.
Since this is an optomechanical system, it consists of two harmonic oscillators. Therefore to begin, we we need to specify how many Fock states we should consider for both the cavity (Nc) and mechanical (Nm) resonators. Here we pick Nc=6 and Nm=45. Originally I did Nc=4 and got the wrong answer. After talking with the authors of the paper, I adjusted my values to the same values that they used. All of the couplings, and environmental decay rates can be obtained directly from the paper (in units of the mechanical frequency wm). To find the steady-state of the system in QuTiP all we need to do is define the Hamiltonian and collapse operators and call:
rho_ss=steadystate(H,c_op_list)
where c_op_list is the list of collapse operators, and rho_ss is the steady state density matrix. To generate subplots (a)-(h), all we need to do is calculate this steady-state for each value of detuning and coupling strength given in the paper. Things like the Fano factor (part f) and Wigner functions can be calculated directly from this state.
The most difficult part to reproduce was part (i) from the figure. To generate this figure, I initially broke up the region of interest, detuning from -1 to +1 and coupling from g0/kappa=0.8 to 1.2, into a grid of points and calculated the Wigner function at each of these locations. From these Wigner functions, we can extract the minimum value, and the ratio of the number of negative Wigner values to number of positive values. These two numbers give the criteria for non-classicality used in the paper. Per their definition, for a state to be nonclassical, the magnitude of the minimal Wigner value should greater than 5% of the maximum Wigner value, and the number of negative Wigner elements should be >= 3% of the number of positive elements.
The first of these criteria is a valid numerical benchmark, however physically, having just one point that is negative doesn't mean much for the entire state. Therefore, the authors supplemented this with the much more physically relevant condition that the total negative Wigner area should be an appreciable part of the total. The problem is that, as it stands, this criteria is subject to numerical errors. This is because values that should be identically zero, sometimes pickup small positive or negative components at the end of a numerical calculation. Therefore little values such at +/- 1e-15 can distort the ratio of negative to positive Wigner terms. To get around this, I set a threshold of tol=1e-8 such that any Wigner value -1e-8 <= w <= +1e-8 is set to zero. This is quite a liberal bound for numerical errors but works none the less.
After setting up the calculation, the results did not agree with those of the original paper. After contacting the lead author, I found out why. First, in generating part (i), the authors change the power of the input drive called alpha_L. However, in the caption of the figure, and the text of the paper, the authors say this parameter is fixed. Therefore, the reader would have no idea that they vary this parameter. Furthermore, the equation that they use to vary the laser power is not in the paper, nor in the supplementary material. The authors plan to fix this in an errata. In addition, the authors do not worry about the numerical errors from small non-zero terms. To generate part (i) the authors calculated the states at several points and then, by eye, determined whether they satisfy the non-classicallty constrains. Subfigure (i) plots the boundaries of the regions found via this method. Having fixed our simulation to account for the changing input power alpha_L, to obtain sufficiently similar results to the original paper, it was found that our tolerance value for zeros must be set to tol=1e-4. This is equivalent to saying that, in an experiment, one can only measure values with magnitude greater than this threshold.
The final result is presented in the QuTiP generated figure above. The final thing to consider is the extra region on non-classicality near Delta/wm=+1. It could be that the authors missed this section using their search method. An indeed, using their parameters this seems to be the case. However, this extra region of non-classicallty is actually a numerical artifact resulting from not considering enough basis states for the cavity resonator at these larger driving powers (again this fact is missing in the paper). Increasing the number of basis states to Nc=8 and Nm=60 makes this region disappear.
What is the end message? Well, what originally was planned to be a single days worth of calculations turned out to be a weeks worth of emailing with the authors and adjusting parameters to try to reproduce their results. Had they included their source code with the paper, this could have been avoided. The author stated that they did not have any room to discuss computational details in the manuscript, however I find this answer pretty weak given that PRL has supplementary material, which the authors used, made explicitly for details such as these. After discussing with the authors, I now understand what they were trying to show with their calculations, and was actually quite impressed with the knowledge of the main author with respect to numerical simulations. Unfortunately, a lot of this didn't come through in the paper. Hopefully when the errata get published this will all become clear.
QuTiP Script for Generating Part (i) from Original Plot:
"""
This script regenerates the data for Figure #1 part (i) from:
Qian et. al. PRL 109, 253601 (2012).
This script is meant to be used in collaboration with the:
opto_figure.py
script that plots all of Figure #1. This file MUST BE run first.
"""
from scipy import *
from qutip import *
#parameters (all rates in units of wm)
#-------------------------------------
Nc=6 #Number of cavity states
Nm=45 #Number of mechanical states
kappa=0.3 #Cavity damping rate
gamma=0.00147 #Mech damping rate
tol=1e-4 #sets limit for what is considered zero (measureable).
#operators
idc=qeye(Nc)
idm=qeye(Nm)
a=tensor(destroy(Nc),idm)
b=tensor(idc,destroy(Nm))
#collapse operators
cc=sqrt(kappa)*a
cm=sqrt(gamma)*b
c_op_list=[cc,cm]
#number of steps in detuning (delta) direction
steps=21
delta=linspace(-1.0,1.0,steps)
#range and number of steps for Wigner function
xvec=linspace(-10,10,500)
#steady state solver function for parfor
#---------------------------------------
def func(args):
x=args[0];g0=args[1]
#change laser driving strength, this is omitted in the paper :(
alpha=sqrt(0.0125/(g0**2))
#-------------------------
H=(-delta[x]+g0*(b.dag()+b))*(a.dag()*a)+b.dag()*b+alpha*(a.dag()+a)
rho_ss=steadystate(H,c_op_list)
rho_b=ptrace(rho_ss,1)
W=wigner(rho_b,xvec,xvec)
neg_inds=where(W<-tol) #find all negative values
pos_inds=where(W>tol) #find all positive values
zero_inds=where(abs(W)<tol) #find all zero values
min_val=W.min() #get minimum Wigner value
num_neg=len(W[neg_inds[0],neg_inds[1]]) #number of neg values
num_pos=len(W[pos_inds[0],pos_inds[1]]) #number of pos values
#determine if state is nonclassical according to definition
if min_val<=-0.05*W.max() and num_neg>=0.03*num_pos:
nonclass=1 #yes
else:
nonclass=0 #no
return nonclass
#number of steps in the coupling (G) direction
steps2=21
G=linspace(0.8,2.0,steps2)*kappa
non_class=zeros((steps2,steps))
#loop over the array of coupling values
for kk in range(steps2):
args=array([[jj,G[kk]] for jj in range(steps)])
results=parfor(func,args)
non_class[kk,:]=results
print kk, " done" #gives some feedback as to progress of calculation
#use QuTiP's store to file function
file_data_store("opto_data_%s_%s.dat" % (steps2,steps),non_class,numtype='real')
QuTiP Script for Generating the Entire Figure (After running the above script):
"""
This script regenerates Figure #1 from:
Qian et. al. PRL 109, 253601 (2012).
To use this script, you MUST first run its companion script:
opto_regions.py
that generates the required data file:
opto_data_21_21.dat
corresponding the the nonclassically tests for part (i)
in the original PRL Figure.
"""
from scipy import *
from qutip import *
from pylab import *
#parameters (all rates in units of wm)
#-------------------------------------
Nc=6 #Number of cavity states
Nm=45 #Number of mechanical states
alpha=0.311 #Coherent state amplitude
g0=0.36 #Coupling strength
kappa=0.3 #Cavity damping rate
gamma=0.00147 #Mech damping rate
xvec=linspace(-10,10,500) #Wigner area and sampling
#operators
idc=qeye(Nc)
idm=qeye(Nm)
a=tensor(destroy(Nc),idm)
b=tensor(idc,destroy(Nm))
#Hamiltonian terms that do not depend on detuning or coupling
H_pre=b.dag()*b+alpha*(a.dag()+a)
#collapse operators
cc=sqrt(kappa)*a
cm=sqrt(gamma)*b
c_op_list=[cc,cm]
#array of points we specifically want
elems=[-0.8,-0.6,-0.25,0.0]
#setup array of detunings to iterate over
delta=linspace_with(-1.0,0.2,20,elems)
#get indicies of elems in delta
inds=array([i for i in range(len(delta)) if delta[i] in elems])
#steady state solver routine
#---------------------------
def func(x):
H=(-delta[x]+g0*(b.dag()+b))*(a.dag()*a)+H_pre
rho_ss=steadystate(H,c_op_list)
numa=expect(a.dag()*a,rho_ss)
rho_b=ptrace(rho_ss,1)
numb=expect(num(Nm),rho_b)
numb2=expect(num(Nm)**2,rho_b)
fano=(numb2-numb**2)/numb
print delta[x]," done"
return numa,numb,fano,rho_b
#run steadystate in parallel over delta array
results=parfor(func,range(len(delta)))
#Make plot of results
#--------------------
#setup a plot grid using gridspec so we can
#adjust the spacing between subplots
import matplotlib.gridspec as gridspec
gs = gridspec.GridSpec(3,9)
#set subplot spacings
gs.update(wspace=0.75,hspace=0.40)
#add all wigner function plots by iterating over inds
for kk in range(4):
ind=inds[kk]
W=wigner(results[3][ind],xvec,xvec)
wmap=wigner_cmap(W)
if kk!=3:
ax=subplot(gs[kk,0])
contourf(xvec,xvec,W,50,cmap=wmap)
ax.xaxis.set_major_locator(MaxNLocator(3))
ax.yaxis.set_major_locator(MaxNLocator(3))
else:
ax=subplot(gs[0:2,1:3])
ax.xaxis.set_major_locator(MaxNLocator(3))
ax.yaxis.set_major_locator(MaxNLocator(3))
wig_den=contourf(xvec,xvec,W,50,cmap=wmap)
title("Wigner Density")
cbar=colorbar(wig_den)
cbar.set_ticks(linspace(round(W.min(),2),round(W.max(),2),6))
#add plot for phonon occupation number
ax=subplot(gs[2:3,1:3])
ax.xaxis.set_major_locator(MaxNLocator(4))
ax.yaxis.set_major_locator(MaxNLocator(4))
plot(delta,results[1],'g',lw=1.5)
plot(delta[inds],results[1][inds],'ko')
xlabel(r"Detuning $\Delta/\omega_{m}$")
ylabel(r"Phonon Number")
xlim([-1.0,0.2])
#add plot for Fano factor
ax=subplot(gs[0:1,3:7])
ax.xaxis.set_major_locator(MaxNLocator(4))
ax.yaxis.set_major_locator(MaxNLocator(4))
plot(delta,results[2],'r',lw=1.5)
plot(delta[inds],results[2][inds],'ko')
xlabel(r"Detuning $\Delta/\omega_{m}$")
ylabel(r"Fano Factor")
axhline(y=1,lw=1)
xlim([-1.0,0.2])
#add plot for P(n) at g0=0.36
ax=subplot(gs[1:,3:5])
ax.xaxis.set_major_locator(MaxNLocator(4))
ax.yaxis.set_major_locator(MaxNLocator(4))
plot(range(Nm),results[3][inds[-1]].diag(),'k',lw=1.5)
plot(range(Nm),results[3][inds[-2]].diag(),'r--',lw=1.5)
ax.xaxis.set_major_locator(MaxNLocator(5,integer=True))
ax.yaxis.set_major_locator(MaxNLocator(4))
xlabel(r"$n$")
ylabel(r"$P(n)$")
xlim([0,Nm-1])
#generate steady state solutions at larger coupling strength g0=0.6
#required from reproducing part (h) of original plot
g0=0.6 #stronger coupling
alpha=0.186 #laser driving strength
#use delta = -0.5 and get steady state
delta=-0.5
H05=(-delta+g0*(b.dag()+b))*(a.dag()*a)+b.dag()*b+alpha*(a.dag()+a)
rho05=steadystate(H05,c_op_list)
#mechanical DM at delta=-0.5
rhob5=ptrace(rho05,1)
#change detuning to delta = -0.25 and get steady state
delta=-0.25
H025=(-delta+g0*(b.dag()+b))*(a.dag()*a)+b.dag()*b+alpha*(a.dag()+a)
rho025=steadystate(H025,c_op_list)
#mechanical DM at delta=-0.25
rhob25=ptrace(rho025,1)
#add new results to plot
ax=subplot(gs[1:,5:7])
ax.xaxis.set_major_locator(MaxNLocator(5,integer=True))
ax.yaxis.set_major_locator(MaxNLocator(4))
plot(range(Nm),rhob5.diag(),'r--',lw=1.5)
plot(range(Nm),rhob25.diag(),'k',lw=1.5)
xlim([0,Nm-1])
xlabel(r"$n$")
ylabel(r"$P(n)$")
#Plot data generated from parameter space plot
#of non-classicality of the mechanical resonator.
#get data from file
data=file_data_read('opto_data_21_21.dat')
delta_from_data=linspace(-1.0,1.0,len(data[0,:]))
coupling_from_data=linspace(0.8,2.0,len(data[:,0]))
ax=subplot(gs[0:,7:])
ax.set_xlim(0.8,2.0)
ax.set_ylim(-1.0,1.0)
#add markers for points A,B,C,D
plot(1.2*ones(4),[0,-0.25,-0.6,-0.8],'ro',markersize=10)
text(1.25, 0, "D",ha='left', va='bottom',fontsize=18, color='g')
text(1.25, -0.25, "C",ha='left', va='bottom',fontsize=18, color='g')
text(1.25, -0.6, "B",ha='left', va='bottom',fontsize=18, color='g')
text(1.25, -0.8, "A",ha='left', va='bottom',fontsize=18, color='g')
#plot data
im=imshow(flipud(data.T),cmap=mpl.cm.binary)
im.set_extent([0.8,2.0,-1,1])
ax.xaxis.set_major_locator(MaxNLocator(4))
ax.yaxis.set_major_locator(MaxNLocator(4))
xlabel(r"Coupling Strength $g_{0}/\kappa$",fontsize=12)
ylabel(r"Detuning $\Delta/\omega_{m}$",fontsize=12)
#see what we get. Figure needs to be expanded manually (sorry)
show()
