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

# # Simulation of quantum circuits on the QLM: introduction<br><br>
# 
# This notebook aims at introducing **how simulations are run in the QLM**. It is valid for both noisy and ideal circuit simulators. <br><br>
# 
# - When simulating a quantum circuit, one can aim for different kinds of results:<br><br>
# 
#     - You might want to get the **full amplitude vector** back, e.g for detailed analysis and debugging of your circuit implementation.<br><br>
#     - On the contrary, you may be interested in **strictly emulating** the behavior of a quantum computer, and getting a list of **measurement results** as a result of your computation.<br><br>
#     - Finally, perhaps you are only interested in the average value of an **observable** at the output of your circuit. It might be useful when dealing, for instance, when dealing with hybrid variational approaches (QAOA, VQE...).<br><br>
#     
# 
# - **Table of contents of this notebook**
#     - [Introduction of the process, with default options](#default)
#     - [Working with a subset of qubits](#subset)
#     - [Getting a list of measurement results](#emulation) 
#     - [Grouping measurement results by value: aggregate_data](#aggregate)
#     - [Getting the average value of an observable](#observable)
#     - [A little hack: directly getting the table of amplitudes as a numpy array](#wavefunction)
#     
# Within this notebook, we only work with a simple Bell-state-creation quantum circuit, and simulate it with our generic simulator **PyLinalg** (based on linear algebra).<br><br>
# 
# For more details about **how to write observables** see [this notebook](../basic/observables.ipynb) and the [Sphinx Documentation](../../doc/index.html)

# <a id="default"></a>
# ## Overall process
# 
# <br><br>A simulation is started by sending a **job** to a **qpu** (quantum processing unit) via its **submit** method. 
# 
# In our case, a qpu is a simulator. The job is created from a quantum circuit. Simulation options are specified at the creation of this job.<br><br>
# 
# The following snippet example of the process:

# In[ ]:


from qat.lang.AQASM import Program, H, CNOT
from qat.qpus import PyLinalg

# qpu creation
qpu = PyLinalg() 

# program creation and gate applications
my_prog = Program()
qbits = my_prog.qalloc(2)
my_prog.apply(H,qbits[0])
my_prog.apply(CNOT,qbits)

# converting into a circuit
circ = my_prog.to_circ()

job = circ.to_job() # specify simulation options. Here: default.

result = qpu.submit(job)

for state in result:
    print(state)


# ## What happened here ?
# 
# "result" contains **all states with non-zero amplitude**. <br><br> 
# 
# The job was created with default arguments. Which in particular that the "nbshots" argument was equal to 0 (see docstring below).
# 
# - Concerning nbshots, the convention regarding the behavior depending on its value is:<br><br>
#     - **if nbshots = 0** then the qpu returns the **best it can do**. In the case of LinAlg the best it can do is just to return the **full state distribution**. Because states are, by default, filtered with an **amplitude threshold** (see docstring below), the result here consists in the two possible outputs, for a Bell state. It might be the case, for instance with an actual quantum chip, that "the best the qpu can do" is to return a list of measurement result, as long as reasonably possible.<br><br>
#     - **if nbshots > 0** then the qpu returns a **list of samples**, obtained by measuring the output probability distribution of the circuit. As we will see later, in that case, the format of the output also depends on the "aggregate_data" argument.<br><br>
#     
# As you can see, it is also at job creation that one can specify the **subset of qubits** of interest. Apart from the number of qubits in the returned samples, as we will see, the behavior is strictly identical to the default case where all qubits are measured. 

# In[ ]:


help(circ.to_job)


# #### We are now going to play around with all these arguments, in order to illustrate the behaviors they lead to.
# 
# 

# <a id="subset"></a>
# ## Focusing on a subset of qubits

# In[ ]:


job = circ.to_job(qubits=[0]) # focusing on the first qubit. 

result = qpu.submit(job)

for state in result:
    print(state) # 1-qubit states, because we focus on one qubit.


# Or, equivalently, but more general (e.g when working with **several registers** of qubits):

# In[ ]:


job = circ.to_job(qubits=[qbits[0].index]) # focusing on the first qubit. the index contains the 
                                           # index of the qubit within the entire set of qubits,    
                                           # composed of all registers.
result = qpu.submit(job)

for state in result:
    print(state)


# <a id="emulation"></a>
# ## Strict emulation: getting a list of samples from the output probability distribution.
# 
# By "strict emulation", we mean that, even though here we are using a classical simulator to process our quantum circuits, the result looks **exactly** like the **raw output** of an **actual quantum computer**.

# In[ ]:


job = circ.to_job(nbshots=10, aggregate_data=False) # see below for aggregate_data.

result = qpu.submit(job)

for state in result: # 10 results are printed, as required.
    print(state) 


# <a id="aggregate"></a>
# ## Grouping results by value: aggregate_data option (defaults to True)<br><br>
# 
# When this option is active (which is the default), all samples containing a given state are "aggregated" together. The result object then contains one unique sample per possible output. It comes with an **empirical estimation of the probability of that output**, and an "error" associated to that estimation. 

# In[ ]:


job = circ.to_job(nbshots=10) # see below for aggregate_data.

result = qpu.submit(job)

for state in result: # 10 results are printed, as required.
    print(state) 


# **Notice how the error decreases when taking a greater number of samples:** TODECIDE: to we remove this or do we implement the error estimation for linalg ?

# In[ ]:


job = circ.to_job(nbshots=1000) # see below for aggregate_data.

result = qpu.submit(job)

for state in result: # 10 results are printed, as required.
    print(state) 


# <a id="observable"></a>
# ## Average value of an observable
# 
# Often, for instance when dealing with hybrid variational algorithms, we are actually interested in the **average value of an observable**, and not measurement values.
# 
# Here, we only give a simple example of the "ZZ" observable, which is always equal to 1 for a bell pair. For more advanced used of observables, and details regarding their construction, we refer the reader to [Sphinx Documentation](../../doc/index.html) and [this notebook](../basic/observables.ipynb).

# In[ ]:


from qat.core import Observable, Term

obs = Observable(2., pauli_terms=[Term(1., "ZZ", [0, 1])])
job = circ.to_job("OBS", observable=obs)

result = qpu.submit(job)

print("Should be 1, as we are working with a Bell pair: ", result.value)


# <a id="wavefunction"></a>
# ## Getting the tables of amplitude directly, as a numpy array: wavefunction

# In[ ]:


from qat.core.simutil import wavefunction

qpu = PyLinalg() 

wf = wavefunction(circ, qpu)

print(wf) # simple and efficient.


# In[ ]: