Documentation: "Adding Decoherence Noise" Example Needs Changes to Work

Pre-Report Checklist

• I am running the latest versions of pyQuil and the Forest SDK
• I checked to make sure that this bug has not already been reported

Issue Description

The example at:

https://pyquil-docs.rigetti.com/en/stable/noise.html?highlight=Adding%20Decoherence%20Noise#adding-decoherence-noise

requires two changes to work correctly.

How to Reproduce

The code snippet below contains the changes needed to make the example work.

The changes are:

 ####CORRECTION1 - Replace following commented line with the following line
        #noisy = add_decoherence_noise(prog, T1=t1).inst([
        noisy = add_decoherence_noise(prog, T1=t1, T2=2*t1).inst([
            Declare("ro", "BIT", 2),
            MEASURE(0, ("ro", 0)),
            MEASURE(1, ("ro", 1)),
        ])

and

  ####CORRECTION2 - Add "noisy.wrap_in_numshots_loop(1000)"
        noisy.wrap_in_numshots_loop(1000)

If CORRECTION1 is not made errors will result.

If CORRECTION2 is not made the output will be incorrect.

Code Snippet

# general imports
import matplotlib.pyplot as plt
%matplotlib inline
import string
import time
import numpy as np
import random as rd

import json
import os
from numpy import pi

from pyquil import Program, get_qc
from pyquil.api import QVM
from pyquil.gates import RX, MEASURE, H, S, CNOT, S, RZ, RY, CZ, Declare
from pyquil.quil import Pragma
from pyquil.noise import add_decoherence_noise

def get_compiled_prog(theta):
    return Program([
        RZ(-pi/2, 0),
        RX(-pi/2, 0),
        RZ(-pi/2, 1),
        RX( pi/2, 1),
        CZ(1, 0),
        RZ(-pi/2, 1),
        RX(-pi/2, 1),
        RZ(theta, 1),
        RX( pi/2, 1),
        CZ(1, 0),
        RX( pi/2, 0),
        RZ( pi/2, 0),
        RZ(-pi/2, 1),
        RX( pi/2, 1),
        RZ(-pi/2, 1),
    ])

import numpy as np
from pyquil import get_qc
qc = get_qc("2q-qvm")

t1s = np.logspace(-6, -5, num=3)
thetas = np.linspace(-pi, pi, num=20)
t1s * 1e6 # us

from pyquil.noise import add_decoherence_noise
records = []
for theta in thetas:
    for t1 in t1s:
        prog = get_compiled_prog(theta)
        ####CORRECTION1 - Replace following commented line with the following line
        #noisy = add_decoherence_noise(prog, T1=t1).inst([
        noisy = add_decoherence_noise(prog, T1=t1, T2=2*t1).inst([
            Declare("ro", "BIT", 2),
            MEASURE(0, ("ro", 0)),
            MEASURE(1, ("ro", 1)),
        ])
        ####CORRECTION2 - Add "noisy.wrap_in_numshots_loop(1000)"
        noisy.wrap_in_numshots_loop(1000)
        bitstrings = qc.run(noisy).readout_data.get("ro")

        # Expectation of Z0 and Z1
        z0, z1 = 1 - 2*np.mean(bitstrings, axis=0)

        # Expectation of ZZ by computing the parity of each pair
        zz = 1 - (np.sum(bitstrings, axis=1) % 2).mean() * 2

        record = {
            'z0': z0,
            'z1': z1,
            'zz': zz,
            'theta': theta,
            't1': t1,
        }
        records += [record]

%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
sns.set(style='ticks', palette='colorblind')

import pandas as pd
df_all = pd.DataFrame(records)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12,4))

for t1 in t1s:
    df = df_all.query('t1 == @t1')

    ax1.plot(df['theta'], df['z0'], 'o-')
    ax2.plot(df['theta'], df['z1'], 'o-')
    ax3.plot(df['theta'], df['zz'], 'o-', label='T1 = {:.0f} us'.format(t1*1e6))

ax3.legend(loc='best')

ax1.set_ylabel('Z0')
ax2.set_ylabel('Z1')
ax3.set_ylabel('ZZ')
ax2.set_xlabel(r'$\theta$')
fig.tight_layout()               

Error Output

First error seen if CORRECTION1 is not made:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-28-8002bb82ebb6> in <module>
      6        ####CORRECTION1 - UNCOMMENT TO FIX
      7        #noisy = add_decoherence_noise(prog, T1=t1, T2=2*t1).inst([
----> 8        noisy = add_decoherence_noise(prog, T1=t1).inst([
      9            Declare("ro", "BIT", 2),
     10            MEASURE(0, ("ro", 0)),

/opt/conda/lib/python3.8/site-packages/pyquil/noise.py in add_decoherence_noise(prog, T1, T2, gate_time_1q, gate_time_2q, ro_fidelity)
    626     """
    627     gates = _get_program_gates(prog)
--> 628     noise_model = _decoherence_noise_model(
    629         gates,
    630         T1=T1,

/opt/conda/lib/python3.8/site-packages/pyquil/noise.py in _decoherence_noise_model(gates, T1, T2, gate_time_1q, gate_time_2q, ro_fidelity)
    461         ro_fidelity = {q: ro_fidelity for q in all_qubits}
    462 
--> 463     noisy_identities_1q = {
    464         q: damping_after_dephasing(T1.get(q, INFINITY), T2.get(q, INFINITY), gate_time_1q) for q in all_qubits
    465     }

/opt/conda/lib/python3.8/site-packages/pyquil/noise.py in <dictcomp>(.0)
    462 
    463     noisy_identities_1q = {
--> 464         q: damping_after_dephasing(T1.get(q, INFINITY), T2.get(q, INFINITY), gate_time_1q) for q in all_qubits
    465     }
    466     noisy_identities_2q = {

/opt/conda/lib/python3.8/site-packages/pyquil/noise.py in damping_after_dephasing(T1, T2, gate_time)
    344         if T1 != INFINITY:
    345             if T2 > 2 * T1:
--> 346                 raise ValueError("T2 is upper bounded by 2 * T1")
    347             gamma_phi -= float(gate_time) / float(2 * T1)
    348 

ValueError: T2 is upper bounded by 2 * T1 

Environment Context

Operating System:

Python Version (python -V):

Quilc Version (quilc --version):

QVM Version (qvm --version):

Python Environment Details (pip freeze or conda list):

Copy and paste the output of `pip freeze` or `conda list` here.

Example of correct output:

Open the linked pyQuil noise documentation and locate the “Adding Decoherence Noise” example. Run the snippet to verify both reported corrections, then update the example so it executes without the T2 error and produces the expected output.