The assigned state matrix in ReferenceDensitySimulator is not persistent between rounds of run

[joshcombes]

Summary: If you initialize the quantum computer object with a density (e.g. qc.qam.wf_simulator.density = rho0 for some valid state rho0) the state is cleared after one call of qc.run.

For example suppose the program is

p = Program(I(0))
ro = p.declare('ro', 'BIT', 1)
p += MEASURE(0, ro[0])

Next load the pure state rho0= diag[0.0, 1.0] into the simulator and call run many times. You should get statistics like Pr(0) ~ 0.0 and Pr(1)~ 1.0. Instead you will get statistics like Pr(0) ~ 0.999 and Pr(1)~ 0.001 if you ran a thousand times. This happens because after the first run the reset method sets the state back to |0> rather than rho0.

I am arguing that this is incorrect behavior as the user has set the initial state it should reset to the state specified by the user until the user decides to change the initial state.

Details: There are some subtitles concerning if the actual behavior is reasonable / expected or not. So lets do a shallow belly flop (cf. deep dive) into this issue. The points made here are illustrated with working code below.

Let me state again that the desired statistics are approximately Pr(0) = 0.0 and Pr(1) = 1.0.

• qc.run behavior. The state is reset to |0> after every run. So if you put run(prog) into a loop the statistics will be Pr(0) = 1.0 and Pr(1) = 0.0.
• wrap_in_numshots_loop behavior. Now the statistics are correct i.e. Pr(0) = 0.0 and Pr(1) = 1.0. But if the program is run again then the statistics are incorrect.
• qc.run_and_measure behavior. This exhibits the same behavior as wrap_in_numshots_loop.
• operator_estimation. This feature perhaps best represents the expectation of a user. Suppose we want to estimate the expectation of Pauli X, Y, and Z of our program in a given state then we can use the tools in the module operator_estimation to do so. However because the QVM resets back to |0> the expectation values are incorrect.

Moreover, this behavior is, from a naive users perspective, inconsistent.

Code:

# make an abstract compiler
from pyquil import Program
from pyquil.api._qac import AbstractCompiler
class DummyCompiler(AbstractCompiler):
    def get_version_info(self):
        return {}
    def quil_to_native_quil(self, program: Program):
        return program
    def native_quil_to_executable(self, nq_program: Program):
        return nq_program

#make a quantum computer object
import networkx as nx
from pyquil.device import NxDevice
from pyquil.api import QVM, QuantumComputer
from pyquil.pyqvm import PyQVM
from pyquil.reference_simulator import ReferenceDensitySimulator

device = NxDevice(nx.complete_graph(1))
qc = QuantumComputer(name='testy!',
                     qam=PyQVM(n_qubits=1,quantum_simulator_type=ReferenceDensitySimulator),
                     device=device,
                     compiler=DummyCompiler())

# inspect the state of the QVM
print(qc.qam.wf_simulator.density)

# density matrix 
import numpy as np
rho1 = np.array([[ 0.0+0.j, 0.0+0.0j], [0.0-0.0j,  1.0+0.j]])

# load state into QVM
qc.qam.wf_simulator.density = rho1
# display QVM state
print(qc.qam.wf_simulator.density)

# Make some programs to run
from pyquil.gates import I, MEASURE
# run prog
prog = Program(I(0))
ro = prog.declare('ro', 'BIT', 1)
prog += MEASURE(0, ro[0])
# wrap in numshots prog
progwrap = prog.copy()
progwrap.wrap_in_numshots_loop(10)
print(prog)
# Run and measure style
progRAM = Program(I(0))
print(progRAM)

# run
qc.qam.wf_simulator.density = rho1
print(qc.run(prog))

for _ in range(0,4):
    print(qc.run(prog)) 
print('=========')
# compare with
qc.qam.wf_simulator.density = rho1
print(qc.run(prog))
for _ in range(0,4):
    qc.qam.wf_simulator.density = rho1
    print(qc.run(prog))

# wrap in numshots
qc.qam.wf_simulator.density = rho1
print(qc.qam.wf_simulator.density)
yep = qc.run(progwrap)
print(yep)
print(qc.qam.wf_simulator.density)
print(qc.run(progwrap))

# Run and measure
qc.qam.wf_simulator.density = rho1
print(qc.qam.wf_simulator.density)
print(qc.run_and_measure(progRAM,trials=40))
print(qc.qam.wf_simulator.density)
print(qc.run_and_measure(progRAM,trials=40))

# operator estimation
from forest.benchmarking.tomography import generate_state_tomography_experiment
experiment_1q = generate_state_tomography_experiment(program=progRAM, qubits=[0])
print(experiment_1q)
from pyquil.operator_estimation import measure_observables
qc.qam.wf_simulator.density = rho1
results = list(measure_observables(qc=qc, tomo_experiment=experiment_1q, n_shots=40000))
print(results)

from forest.benchmarking.tomography import iterative_mle_state_estimate

rho = iterative_mle_state_estimate(results=results, qubits=[0])
print('rho est', np.around(rho, decimals=2))
print('===')
print('rho true', np.around(rho1, decimals=2))

[joshcombes]

@kilimanjaro @kylegulshen, and @msohaibalam do you have thoughts?

[joshcombes]

closed by merging https://github.com/rigetti/pyquil/pull/920