pyquil icon indicating copy to clipboard operation
pyquil copied to clipboard
verified publisher icon rigetti

Parameters are incorrectly assumed to be 2𝜋 periodic

Open bramathon opened this issue 3 years ago • 2 comments

Pre-Report Checklist

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

Issue Description

The headline summary is that the XY gate works as expected when provided a numerical value, but not when provided a parameter.

before-andy-fix

This behaviour can be traced to an incorrect assumption about parameters deep in the stack, and may affect other gates (although we haven't yet found any).

To reproduce:

Code Snippet

import numpy as np
from pyquil.api import get_qc
from pyquil.quil import Program
from pyquil.gates import RX, XY, MEASURE

qc = get_qc("Aspen-M-3")
qvm = get_qc("Aspen-M-3", as_qvm=True)

program = Program()
theta = program.declare("theta", "REAL", 1)
ro = program.declare("ro", "BIT", 2)
program += RX(np.pi/2, 112)
program += XY(theta, 112, 125)
program += RX(np.pi/2, 112)
program += MEASURE(112, ro[0])
program += MEASURE(125, ro[1])
program.wrap_in_numshots_loop(512)
exe = qc.compiler.native_quil_to_executable(program)
qvm_exe = qvm.compiler.native_quil_to_executable(program)

results = []
qvm_results = []
x = np.linspace(-2*np.pi, 2*np.pi, num=24)
for theta in x:
    # QPU
    exe.write_memory(region_name="theta", value=[theta])
    res = qc.run(exe).readout_data["ro"]
    states, counts = np.unique(res, return_counts=True, axis=0)
    hist = {"".join(str(s) for s in state): c for state, c in zip(states, counts)}
    # print(f"QPU {hist}")
    results.append(hist)
    
    # QVM
    qvm_exe.write_memory(region_name="theta", value=[theta])
    sim_res = qvm.run(qvm_exe).readout_data["ro"]
    states, counts = np.unique(sim_res, return_counts=True, axis=0)
    hist = {"".join(str(s) for s in state): c for state, c in zip(states, counts)}
    # print(f"QVM {hist}")
    qvm_results.append(hist)

import matplotlib.pyplot as plt

y = [res["10"] if "10" in res else 0 for res in results]
y_qvm = [res["10"] if "10" in res else 0 for res in qvm_results]
plt.plot(x, y, "-o", label="QPU")
plt.plot(x, y_qvm, "-o", label="QVM")
plt.legend(loc="upper left")
plt.title("XY(theta). QPU (theta), QVM(theta)")

How does this happen?

First, we observe that the behaviour of the QPU is consistent with the angle being mod 2𝜋'd. In other words, we can reproduce the same behaviour of the simply by replacing:

qvm_exe.write_memory(region_name="theta", value=[theta])

with

qvm_exe.write_memory(region_name="theta", value=[theta % 2*np.pi])

The XY(theta) is implemented as a 0.5*%theta phase shift on the xy frame in the defcalibration. When the calibration is expanded in the program instruction sequence, we thus end up with something like:

DECLARE theta REAL[1]
DECLARE ro BIT[2]
RX(pi/2) 112
XY(theta) 112 125
RX(pi/2) 112

becoming,

DECLARE __P2 REAL[1]
DECLARE theta REAL[1]
DECLARE ro BIT[2]
FENCE 112
NONBLOCKING PULSE 112 "rf" drag_gaussian(duration: 2.4000000000000003e-08, fwhm: 6.000000000000001e-09, t0: 1.2000000000000002e-08, anh: -190000000.0, alpha: -2.4831333182584574, scale: 0.312478348136509, phase: 0.0, detuning: -64046.20547445084)
FENCE 112
FENCE
SHIFT-PHASE 112 125 "xy" -1.9162803766959502
NONBLOCKING PULSE 112 125 "xy" q112_q125_xy/sqrtiSWAP
SHIFT-PHASE 112 125 "xy" 1.9162803766959502
SHIFT-PHASE 112 125 "xy" -0.7901340066069987 + -0.5*__P2[0]
NONBLOCKING PULSE 112 125 "xy" q112_q125_xy/sqrtiSWAP
SHIFT-PHASE 112 125 "xy" -1*(-0.7901340066069987 + -0.5*__P2[0])
DELAY 112 "rf" 2.4e-07
SHIFT-PHASE 112 "rf" 3.8186114768472033
DELAY 125 "rf" 2.4e-07
SHIFT-PHASE 125 "rf" 5.57043817144071
SHIFT-PHASE 125 "rf" 0.5*__P2[0]
SHIFT-PHASE 112 "rf" -0.5*__P2[0]
FENCE
SHIFT-PHASE 111 112 "xy" 1.9093057384236016
SHIFT-PHASE 112 113 "xy" 1.9093057384236016
SHIFT-PHASE 112 125 "xy" 1.9093057384236016
SHIFT-PHASE 112 125 "xy" -2.785219085720355
SHIFT-PHASE 124 125 "xy" -2.785219085720355
SHIFT-PHASE 125 126 "xy" -2.785219085720355
SHIFT-PHASE 112 125 "xy" -0.25*__P2[0]
SHIFT-PHASE 124 125 "xy" -0.25*__P2[0]
SHIFT-PHASE 125 126 "xy" -0.25*__P2[0]
SHIFT-PHASE 111 112 "xy" -0.25*__P2[0]
SHIFT-PHASE 112 113 "xy" -0.25*__P2[0]
SHIFT-PHASE 112 125 "xy" -0.25*__P2[0]
FENCE 112
NONBLOCKING PULSE 112 "rf" drag_gaussian(duration: 2.4000000000000003e-08, fwhm: 6.000000000000001e-09, t0: 1.2000000000000002e-08, anh: -190000000.0, alpha: -2.4831333182584574, scale: 0.312478348136509, phase: 0.0, detuning: -64046.20547445084)
FENCE 112

Here, the parameters have been recalculated by rewrite_arithmetic and the calibrations have been expanded. So theta is now __P2[0]. Additionally, in the recalculation table, the parameter is divided by 2𝜋, giving the mapping:

{ParameterAref(index=0, name='__P2'): Div(<MRef theta[0]>,6.283185307179586)}

This is what is submitted to the API.

Somewhere behind the API, __P2 is being modulo'd by 1.0.

For many instructions that is the correct behaviour, but in this case the full range of the parameter is 8𝜋. To be specific:

SHIFT-PHASE 112 125 "xy" -0.25*__P2[0] has a domain of 0 to 8𝜋, but the scaling of the parameter and subsequent modulo treat it as if the domain is 0 to 2𝜋.

Why doesn't this affect numbers?

For real numbers, the parameter expression is evaluated prior to modulo. So the output of -0.25*__P2[0] is periodic in 2𝜋.

Why do XY gates work at all?

The RZs defcals include numerous xy frame updates, so we might wonder how anything works. It seems that the frame updates in RZ may indeed be off by a sign, but without a physical effect.

Possible solutions

  1. It may be necessary for the stack to be more aware of the periodicity of parameters
  2. Redefine theta to 4*theta in the XY gate to ensure the shift-phase instructions are all, in fact 2pi periodic
  3. Evaluate shift-phase expressions inline on the client side
  4. ???

We can see that the definition of theta to 4theta is effective, however, this would break all existing programs using XY.

after-andy-fix

Environment Context

Aspen devices

bramathon avatar Dec 07 '22 21:12 bramathon

Snippet to modify an XY defcal:

calibration_program = qc.compiler.get_calibration_program()
defcal = calibration_program.get_calibration(XY(np.pi, 112, 125))

new_instrs = []
for inst in defcal.instrs:
    if isinstance(inst, ShiftPhase) and isinstance(inst.phase, Expression):
        if isinstance(inst.phase, Add):
            new_instrs.append(ShiftPhase(inst.frame, Add(inst.phase.op1, Mul(4.0*inst.phase.op2.op1, inst.phase.op2.op2))))
        elif isinstance(inst.phase, Mul) and isinstance(inst.phase.op2, Parameter):
            new_instrs.append(ShiftPhase(inst.frame, Mul(4.0*inst.phase.op1, inst.phase.op2)))
        elif isinstance(inst.phase, Mul) and isinstance(inst.phase.op2, Add) and isinstance(inst.phase.op2.op2, Mul):
            new_instrs.append(ShiftPhase(inst.frame, Mul(inst.phase.op1, Add(inst.phase.op2.op1, Mul(4.0*inst.phase.op2.op2.op1, inst.phase.op2.op2.op2)))))
    else:
        new_instrs.append(inst)
                    
defcal.instrs = new_instrs

bramathon avatar Dec 07 '22 21:12 bramathon

Note that the plot titles are incorrect. The QVM is operating here on theta, not theta % 2pi

bramathon avatar Dec 08 '22 01:12 bramathon