#!/usr/bin/env python
# coding: utf-8

# # Direct Randomized Benchmarking

# This tutorial demonstrates how to generate ["Direct randomized benchmarking"](https://arxiv.org/abs/1807.07975) (DRB) circuits using pyGSTi (see the [RB analysis tutorial](RBAnalysis.ipynb) for RB data analysis functions). This is a recently proposed alternative to ["Clifford RB"](http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.106.180504), with the same core aim as Clifford RB - to estimate an error rate that provides a meaure of average gate performance. 
# 
# Clifford RB and Direct RB can be implemented (holistically) on a set of $n$ qubits whenever the $n$-qubit Clifford group can be generated by the gates in the device. But whereas Clifford RB has sequences of uniformly random $n$-qubit Cliffords (which must be compiled into the native gates of the device), a DRB circuit consists of:
# 
# 1. A circuit that generates a uniformly random $n$-qubit stabilizer state.
# 2. $m$ independently sampled layers of the native gates in the device, with these layers sampled according to a user-specified distribution $\Omega$ over all possible circuit layers.
# 3. A circuit that maps the ideal output of the preceeding circuit to a uniformly random computational basis state (or, if preferred, to the all-zeros state).
# 
# This protocol can be implemented on more qubits that Clifford RB, and has similar levels of reliability to Clifford RB (if $\Omega$ is chosen reasonably carefully). 
# 
# One important point to note is that the DRB error rate is $\Omega$-dependent. I.e., it quantifies gate performance over circuits that are sampled according to $\Omega$. This is analogous to the Clifford-compiler dependence of the Clifford RB error rate, but it is more easily controlled and understood. This tutorial will not provide comprehensive details on DRB; see ["Direct randomized benchmarking for multi-qubit devices"](https://arxiv.org/abs/1807.07975) for more information.

# In[1]:


from __future__ import print_function #python 2 & 3 compatibility
import pygsti
from pygsti.extras import rb


# ## Specifying the device to be benchmarked
# 
# To generate DRB circuits, you first need to specify the device to be benchmarked. This step ensures that the circuits returned will respect device connectivity, and contain only gates in the "native" gate-set of the device.
# 
# We do this using a `ProcessorSpec` object: see the [ProcessorSpec tutorial](../objects/advanced/ProcessorSpec.ipynb) on how to create these. Here we'll demonstrate creating DRB circuits for a device with:
# - Five qubits on a ring, labelled 'Q0', ..., 'Q4'
# - 1-qubit gates consisting of $\sigma_x$ and $\sigma_y$ rotations by $\pm \pi/2$, and an idle gate
# - Controlled-Z gates connecting adjacent qubits on the ring
# 
# Below, we generate the `ProcessorSpec` for this device:

# In[2]:


nQubits = 5 
qubit_labels = ['Q0','Q1','Q2','Q3','Q4'] 
gate_names = ['Gxpi2', 'Gxmpi2', 'Gypi2', 'Gympi2', 'Gcphase'] 
availability = {'Gcphase':[('Q0','Q1'), ('Q1','Q2'), ('Q2','Q3'), 
                           ('Q3','Q4'),('Q4','Q0')]}
pspec = pygsti.obj.ProcessorSpec(nQubits, gate_names, availability=availability, 
                                 qubit_labels=qubit_labels)


# ## Generating a Direct RB experiment
# 
# We can generate a set of DRB circuits using the `rb.sample.direct_rb_experiment()` function.
# 
# #### Essential parameters
# To sample a DRB experiment, it is necessary to specify:
# - The "DRB lengths" at which we will sample circuits
# - The number of circuits to sample at each length
# 
# To use our function, it is not *essential* to specify the DRB sampling distribution $\Omega$ (see above), but to make good use of DRB it is important to choose a sampling distribution that matches the information you're trying to obtain by doing DRB. We'll discuss this in more detail below
# 
# A "DRB length" ($m$) is the number of layers in the "core" of the DRB circuit, which consists of $\Omega$-distributed random circuit layers. So it does not include the stabilizer preparation and measurement circuits at either end of a DRB circuit. As with all RB samplers in pyGSTi, the minimal length is $m=0$.
# 
# Let's fix the DRB lengths to 0, 5, 10, 20, 25 and 30 and take the number of circuits at each length to be $k = 10$. (These are not recommendations for these parameters: these circuit lengths are potentially reasonable for 3-qubit DRB, but appropriate choices depend on the approximate quality of the gates; setting $k$ this low is probably not a good idea for an actual experiment, but it suffices here to demonstrate the method).

# In[3]:


lengths = [0,4,10,15,20,25,30]
k = 10


# #### Optional parameters
# 
# The RB samplers in `pyGSTi` allow the user to benchmark a subset of the qubits, by specifying a **`subsetQs`** list. This then means that a `ProcessorSpec` can be specified for an entire device even if you only wish to benchmark some subset of it. If this is not specified it is assumed that you want RB circuits for holistically benchmarking the entire device.
# 
# This set of qubits must be connected (otherwise it is not possible to generate a uniformly random $n$-qubit stabilizer state over these $n$ qubits, which is the first part of a DRB circuit).
# 
# Let's demonstrate generating circuits to benchmark 3 of the qubits:

# In[4]:


subsetQs = ['Q0','Q1','Q2']


# Another important optional parameter is **`randomizeout`**. This specificies whether the perfect output of the circuits should be the input state (assumed to be $0,0,0...$ herein, although any computational basis state can be accounted for) or a random computational basis state. There are many good reasons to instead set this to True, so we'll do that here.

# In[5]:


randomizeout = True


# Another useful parameter is **`citerations`**, which is the number of iterations used in the randomized compilers that construct the stabilizer state preparation and measurement circuits. Increasing this will reduce the average depth of these subcircuits. Note that, because these circuits are not included in the DRB length, reducing their depth effectively reduces the SPAM error in the DRB analysis and improves the estimate of the DRB number.  This contrasts with Clifford RB, for which the benchmarking score *is* compilation dependent.  
# 
# But while more iterations is better from an experimental perspective, any increase will cause the circuit generation computation to take longer to run. For this notebook, we'll leave it at the default value. For the experiments presented in  ["Direct randomized benchmarking for multi-qubit devices"](https://arxiv.org/abs/1807.07975), we increased it to 200.

# In[6]:


citerations = 20


# Next, we'll specify the DRB sampler. There are a few circuit layer samplers built into `pyGSTi`. This includes all the DRB samplers used in the experiments and simulations of ["Direct randomized benchmarking for multi-qubit devices"](https://arxiv.org/abs/1807.07975).
# 
# For all the available options, you can investigate all of the functions beginning `rb.sample.circuit_layer_by_` and/or take a look at `rb.sample.random_circuit()`. Here we'll over-view the simplest option that is valid for any device, as well the most flexible (and likely most useful) option. 
# 
# If an in-built sampler is to be used, it is specified by the setting the optional argument **`sampler`** to the relevant string. Let's set this to `'Qelimination'`, which is the default. Note that this is **not** a sampler that we particularly recommend, but it works with all device connectivities **and** doesn't require any user-input **and** it will result in reliable DRB in most circumstances. So it is a reasonable option for a first attempt at DRB.

# In[7]:


sampler = 'Qelimination'


# This sampler picks a circuit layer in the following way (for more information see the `rb.sample.circuit_layer_by_Qelimination()` docstring). Until all the qubits have a gate acting on them in the layer it repeats the following steps:
# 1. Pick a qubit $q_1$ uniformly at random from the qubits not yet used in the layer.
# 2. Pick another qubit $q_2$, uniformly at random, that is not yet used in the layer **and** that is connected to this qubit, if there are any such qubits.
# 3. If such a qubit has been found, apply a two-qubit gate to this pair of qubits ($q_1$ and $q_2$) with the probability $p$ specified by the user.
# 4. If a two-qubit gate has not been chosen for $q_1$, pick a uniformly random 1-qubit gate to apply to the qubit.
# 
# Each of the samplers have some user-specifiable arguments, which we set via the list **`samplerargs`**. Here, there is only one variable in this sampler: the probability $p$ that appears in step 3. 
# 
# Let's set $p=0.5$ (this is actually the default, but we'll include it explicitly here).

# In[8]:


samplerargs = [0.5]


# To sample a set of DRB circuits using this DRB specification, we simply pass all of these arguments to the `rb.sample.direct_rb_experiment()` function:

# In[9]:


exp_dict = rb.sample.direct_rb_experiment(pspec, lengths, k, subsetQs=subsetQs, sampler=sampler,
                                         samplerargs=samplerargs, randomizeout=randomizeout)


# And that's it!
# 
# Before discussing what's in the output, we'll go throught this again with a different, more flexible, sampler.
# 
# ## A flexible circuit layer sampler: the "compatible two-qubit gates" sampler
# 
# The sampler used above is very simple to specify, but the properties of the layers it samples are fairly opaque (e.g., the expected number of 2-qubit gates depends on device connectivity, as does how often each 2-qubit gate is used). There is an in-built sampler that we have found to be very useful for DRB experiments, the "compatible two-qubit gates" sampler, specified by setting:

# In[10]:


sampler = 'co2Qgates'


# This sampler may seem rather complicated at first - because it's not as simple to specify as the `Qelimination` sampler - but it actually creates much more easily understood circuits. This sampler requires the user to specify sets of compatible 2-qubit gates, meaning 2-qubit gates that can applied in parallel. We specifying this as a list of lists of `Label` objects (see the [Ciruit tutorial](../objects/Circuit.ipynb) for more on `Label` objects), so let's import the `Label` object:

# In[11]:


from pygsti.baseobjs import Label as L


# In this example, we are benchmarking 3 qubits for a device containing 5 qubits with ring connectivity. So we can easily write down all of the possible compatible 2-qubit gate lists over these 3 qubits. There are only 3 of them: a list containing no 2-qubit gates, and two lists containing only one 2-qubit gate:

# In[12]:


C2QGs1 = [] #  A list containing no 2-qubit gates is an acceptable set of compatible 2-qubit gates.
C2QGs2 = [L('Gcphase',('Q0','Q1')),] # A controlled-Z between Q0 and Q1
C2QGs3 = [L('Gcphase',('Q1','Q2')),] # A controlled-Z between Q1 and Q2.


# Note that we often wouldn't want to start by writting down all possible sets of compatible 2-qubit gates - there can be a lot of them. That'll hopefully become clear below. 
# 
# Let's continue with this example, as it is particularly easy to follow. We put all of these possible sets of compatible 2-qubit gates into a list **`co2Qgates`**, we also pick a probability distribution over this list **`co2Qgatesprob`**, and we pick a probability **`twoQprob`** between 0 and 1.

# In[13]:


co2Qgates = [C2QGs1,C2QGs2,C2QGs3]
co2Qgatesprob = [0.5,0.25,0.25]
twoQprob = 1


# The sampler then picks a layer as follows:
# 1. Sample a list from `co2Qgates` according to the distribution `co2Qgatesprob`.
# 2. Consider each gate in this list, and add it to the layer with probability `twoQprob`.
# 3. For every qubit that doesn't yet have a gate acting on it in the layer, independently and uniformly at random, sample a 1-qubit gate to act on that qubit (sampled from the "native" 1-qubit gates in the device).
# 
# So with the example above there is a 50% probability of no 2-qubit gates in a layer, a 50% chance that there is one 2-qubit gate in the layer, there is no probability of more than one 2-qubit gate in the layer (which here is impossible anyway), and each of the two possible 2-qubit gates is equally likely to appear in a layer.
# 
# To clarify this method, note that there is more than one way to achieve the same sampling here. Instead, we could have set `co2Qgatesprob = [0,0.5,0.5]` and `twoQprob = 0.5`.
# 
# To use these sampler parameters, we put them (in this order) into the samplerargs list:

# In[14]:


samplerargs = [co2Qgates,co2Qgatesprob,twoQprob]


# And then we run exactly the same function as before:

# In[15]:


exp_dict = rb.sample.direct_rb_experiment(pspec, lengths, k, subsetQs=subsetQs, sampler=sampler,
                                         samplerargs=samplerargs, randomizeout=randomizeout)


# We've found that it's fairly useful to implement DRB with sampling that picks a single uniformly random 2-qubit gate with some probability and implements 1-qubit gates on all the other qubits (this is useful for fairly few-qubit DRB). And it's a bit inconvenient to specify this in the method above. So, there is also an option for "nested" sets of compatible two-qubit gates. That is, `co2Qgates` can be a list where some or all of the elements are not lists containing compatable two-qubit gates, but are instead lists of lists of compatible two-qubit gates. 
# 
# An element of `co2Qgates` is sampled according to the `co2Qgatesprob` distribution (which defaults to the uniform distribution if not specified). If the chosen element is just a list of `Labels` (i.e., a list of compatible 2-qubit gates), the algorithm proceeds as above. But if the chosen element is a list of lists of `Labels`, the sampler picks one of these sublists uniformly at random; this sublist should be a list of compatible 2-qubit gates.
# 
# This may sound complicated, so below we show how to redo the previous example in this format. 

# In[16]:


co2Qgates = [C2QGs1,[C2QGs2,C2QGs3]]
co2Qgatesprob = [0.5,0.5] # This doesn't need to be specified, as the uniform dist is the default.
twoQprob = 1 # This also doesn't need to be specifed, as this value is the default.
# We leave the latter two values of this list, because we are using the default values.
samplerargs = [co2Qgates,]


# If you want to play around with the various in-built circuit samplers, take a look at the `rb.sample.random_circuit()` function, and the `rb.sample.circuit_layer_by_` functions, which have reasonable docstrings. The  `rb.sample.random_circuit()` function is what is used to sample the "core" circuit in DRB. Note that this function and, in turn, `rb.sample.direct_rb_experiment()`  can also be used with user-defined circuit layer samplers. To use this functionality the `sampler` argument is a user-defined function, and this function needs to take the `ProcessorSpec` and `subsetQs` arguments as the first two inputs and return a circuit layer.

# ## Be careful when picking the sampling distribution!
# 
# Direct RB is not "reliable" with a completely arbitrary sampling distribution, i.e., the observed DRB decay may not be a single exponential and/or the DRB error rate may not be directly related to the error rate of the gates via the formula derived in ["Direct randomized benchmarking for multi-qubit devices"](https://arxiv.org/abs/1807.07975). This is obvious if we consider that one possible sampler would be to deterministically idle every qubit in every layer, and clearly there would be no guarantee of an expoential decay in this case (e.g., in this case large, coherent $\sigma_z$-rotation errors would cause oscillations in the DRB decay curve).
# 
# DRB *is* reliable when the sampling distribution will "scramble" errors fairly quickly, meaning that:
# 1. Unitary errors are quickly converted to stochastic errors (local randomization by 1-qubit gates is sufficient for this).
# 2. Low-weight stochastic errors are (on average) quickly converted to high-weight stochastic errors. (Entangling gates need to be sufficiently often for this to occur).
# 
# For more information on this see ["Direct randomized benchmarking for multi-qubit devices"](https://arxiv.org/abs/1807.07975).

# ## What's in the output?
# 
# The output of `rb.sample.direct_rb_experiment()` has exactly the same format as with the Clifford RB experiment generation function, and this output was explained in detail in the previous Clifford RB tutorial. So see that tutorial for full details; here we only give a brief overview.
# 
# The returned dictionary contains a full specification for the DRB circuits. This dictionary contains 4 keys:

# In[17]:


print(exp_dict.keys())


# - The 'spec' key gives the specification used to generate the circuits, including all of the sampling parameters.
# - The 'circuits' key returns a dictionary containing the sampled circuits.
# - The 'idealout' key returns the error-free output of each circuit, to be compared to observed outputs with "success" declared iff the output is the bit-string specified by 'idealout'.
# - 'qubitordering' specifies the ordering of the qubits in the 'idealout' bit-strings.

# Each of these circuits can be converted to OpenQasm of Quil using the methods shown in the tutorial introducing the `Circuit` object. We can also write them to file using the same method as demonstrated in the [Clifford RB tutorial](CliffordRB.ipynb):

# In[18]:


circuitlist = [exp_dict['circuits'][m,i] for m in lengths for i in range(k)]
pygsti.io.write_circuit_list("../tutorial_files/DirectRBCircuits.txt",circuitlist,
                                "Direct RB circuits")


# Before we finish, let's have a quick look at a couple of the sampled circuits. If we compare the first one of these to a Clifford RB circuit at length 0 (see the end of the Clifford RB tutorial) we'll see that **typically** it contains many fewer 2-qubit gates. The average "cost" of the length 0 circuits is the predominent factor in whether or not the RB method is feasable on a given device: these circuits need to be implementable with a high enough success probability, on average, for an exponential decay to be observable with reasonable amounts of data.

# In[19]:


print("The first circuit sampled at Direct RB length 0:")
print('')
# print(exp_dict['circuits'][0,0])
circuit_string = str(exp_dict['circuits'][0,0])
circuit_string = circuit_string.split('\n')
for bar in range(int(len(circuit_string[0])/80)+1):
    for ind in range(len(circuit_string)):
        print(circuit_string[ind][80*bar:80*(bar+1)])
print("The circuit size is: ", exp_dict['circuits'][0,0].size())
print("The circuit depth is: ", exp_dict['circuits'][0,0].depth())
print("The circuit multi-qubit-gate count is: ", exp_dict['circuits'][0,0].multiQgate_count())


# In[20]:


print("The first circuit sampled at Direct RB length 30:")
print('')
circuit_string = str(exp_dict['circuits'][30,0])
circuit_string = circuit_string.split('\n')
for bar in range(int(len(circuit_string[0])/80)+1):
    for ind in range(len(circuit_string)):
        print(circuit_string[ind][80*bar:80*(bar+1)])
print("The circuit size is: ", exp_dict['circuits'][30,0].size())
print("The circuit depth is: ", exp_dict['circuits'][30,0].depth())
print("The circuit multi-qubit-gate count is: ", exp_dict['circuits'][30,0].multiQgate_count())


# In[ ]: