Python Qiskit and R in RStudio – first experiment

How to use Qiskit (the Python package) to define and run a quantum circuit on an IBMQ computer, then plot the results in R – all in RStudio. Here’s the repo. Or click the button at the top right hand side of this page to download the Markdown file.

I had already installed Anaconda and run:

pip install qiskit[visualization]

R libraries:

library(reticulate)
library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.3      ✔ forcats 0.5.2 ── Conflicts ──────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

{reticulate} is used to get R to communicate with Python.

Now the Qiskit code to get going in a Python chunk:

from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit, IBMQ, assemble, transpile, execute
from numpy import pi
from qiskit.tools import job_monitor

Here’s a simple system with three qubits, two entangled.

\[|\psi\rangle = \frac{|000\rangle + |011\rangle + |100\rangle + |111\rangle}{2}\]

Numbering the qubits right to left, the first two (call them \(q[0]\) and \(q[1]\)) are entangled, the third, \(q[2]\), is unentangled with the first two. We are going to measure \(q[0]\) first, which should cause \(q[1]\) to collapse to the same state but leave \(q[2]\) free to choose at random between 1 and 0. I’ve fiddled around with the sums for this system over here.

Back to Python, to define the circuit:

qreg_q = QuantumRegister(3, 'q')
creg_c = ClassicalRegister(3, 'c')
circuit = QuantumCircuit(qreg_q, creg_c)

circuit.reset(qreg_q[0])
circuit.reset(qreg_q[1])
circuit.reset(qreg_q[2])
circuit.h(qreg_q[0])
circuit.h(qreg_q[2])
circuit.cx(qreg_q[0], qreg_q[1])
circuit.measure(qreg_q[0], creg_c[0])
circuit.measure(qreg_q[1], creg_c[1])
circuit.measure(qreg_q[2], creg_c[2])

Draw the circuit:

circuit.draw('mpl').show()

This will pop up in a new window; I can’t work out how to get it to appear in the Markdown other than by saving it and inserting a png…

Enter your TOP SECRET API token here:

IBMQ.save_account('1337', overwrite=True)

Load up the account:

provider = IBMQ.load_account()

See what backends I can access (I put this on one line as the output appears to be printed before print is called otherwise):

for backend in provider.backends(): print(backend)

ibmq_qasm_simulator
ibmq_lima
ibmq_belem
ibmq_quito
simulator_statevector
simulator_mps
simulator_extended_stabilizer
simulator_stabilizer
ibmq_manila
ibm_nairobi
ibm_oslo

I chose ibmq_lima:

lima = provider.get_backend('ibmq_lima')

(And, note to self – look up how to query queue length.)

Send off the job!

job = execute(circuit,
              backend = lima,
              shots = 1024)
job.job_id()

If RStudio crashes for any reason (and it did), then the job can be retrieved with the ID that would have been printed above before it crashed (or log into your IBM Quantum web account and find the ID there):

job = lima.retrieve_job('635aea27333aeb047289f3ba')

Use this command to query the status of the job:

job.status()

<JobStatus.DONE: 'job has successfully run'>

There’s another command that busy-waits until the job is done. Trying to interrupt it is what crashed RStudio, so I wouldn’t recommend it!

Once the job is done, the results:

res = job.result()

Here are the frequencies of each outcome in Python:

vals = res.get_counts(circuit)
vals

{'000': 266, '001': 21, '010': 10, '011': 247, '100': 248, '101': 13, '110': 6, '111': 213}

We can also see them in R using py$ as a prefix on variable names:

py$vals

{'000': 266, '001': 21, '010': 10, '011': 247, '100': 248, '101': 13, '110': 6, '111': 213}

Now how do I transform that (it’s a dictionary of some description) into a data frame…? I couldn’t see an easy way in R (I didn’t think very long). It was easier at the Python end:

import pandas
df = pandas.DataFrame(list(vals.items()),
                      columns = ['Outcome', 'Freq'])
df

  Outcome  Freq
0     000   266
1     001    21
2     010    10
3     011   247
4     100   248
5     101    13
6     110     6
7     111   213

Back in R:

py$df

And now we can ask ggplot to make a pic!

py$df |> 
  ggplot(aes(x = Outcome, y = Freq, fill = Outcome)) +
  geom_col() +
  theme(legend.position="none")