05_PerformanceChecks.Rmd

---
jupyter:
  jupytext:
    formats: ipynb,Rmd
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
      jupytext_version: 1.3.2
  kernelspec:
    display_name: Python 3
    language: python
    name: python3
---

```{python}
# Visualization
import matplotlib.pyplot as plt
import matplotlib
from tqdm.notebook import tqdm
from time import time
from copy import deepcopy
# %config InlineBackend.figure_format = 'svg' # Makes the images look nice

# Numpy
import numpy as np
from numpy import linalg as LA

#importing cirq
import cirq

# importing Qiskit
from qiskit import QuantumCircuit, execute, Aer
from qiskit.circuit import library as lb
from qiskit.ignis.verification import get_ghz_simple 

# Tensor networks
import quimb as quimb
from ncon import ncon

#User-defined functions in support files
from MPS_QFT.helper import print_state, right_contract, left_contract
from MPS_QFT.helper import to_full_MPS, to_dense, to_approx_MPS

from MPS_QFT.manual import apply_two_qubit_gate_full, max_bond_dimension, apply_two_qubit_gate, apply_one_qubit_gate

from MPS_QFT.gates import CPHASE, cphase_swap_qiskit, cphase_swap_quimb

from MPS_QFT.circuit import qft_circuit_swap_full, qft_circuit_swap_approx
from MPS_QFT.circuit import qft_circuit_qiskit, circ_data, MPS_circ 
```

```{python}
# Plotting LateX figures
matplotlib.use("pgf")
matplotlib.rcParams.update({ 
    "pgf.texsystem": "pdflatex",
    'font.family': 'serif',
    'text.usetex': False,
    'pgf.rcfonts': False,
    "pgf.preamble": [ r"\usepackage[utf8]{inputenc}" ]
})

#Font size configuration
SMALL_SIZE = 8
MEDIUM_SIZE = 10
BIGGER_SIZE = 11
BIGGEST_SIZE = 12

#All sizes are customizable here
plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
plt.rc('axes', titlesize=BIGGER_SIZE)    # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE)   # legend fontsize
plt.rc('figure', titlesize=BIGGEST_SIZE) # fontsize of the figure title
#plt.rcParams['axes.facecolor'] = 'white'
```

```{python}
def get_figsize(wf=0.5, hf=(5.**0.5-1.0)/2.0, ):
    """Parameters:
      - wf [float]:  width fraction in columnwidth units
      - hf [float]:  height fraction in columnwidth units.
                     Set by default to golden ratio.
      - columnwidth [float]: width of the column in latex. Get this from LaTeX 
                             using \showthe\columnwidth
    Returns:  [fig_width,fig_height]: that should be given to matplotlib
    """
    columnwidth = 510.0 #! The width of the Latex paper should be put here [OK]
    
    fig_width_pt = columnwidth*wf 
    inches_per_pt = 1.0/72.27               # Convert pt to inch
    fig_width = fig_width_pt*inches_per_pt  # width in inches
    fig_height = fig_width*hf      # height in inches
    return [fig_width, fig_height]
```

## Quimb-qiskit interface

```{python}
# Controlled NOT
CNOT = quimb.controlled('not')

# Hadamard
H = quimb.gen.operators.hadamard()

# SWAP
SWAP = quimb.gen.operators.swap()

# Dictionary of the gates
gates = {'h': H, 
        'cx': CNOT,
        'cp': CPHASE,
        'swap': SWAP}
```

## Time performance tests


### Quimb - $\chi=2$

```{python}
#run QFT on different number of qubits, store execution times
#quimb implementation with bond dimension chi=2
N_quimb = np.arange(4, 17, dtype=int)

times_quimb_approx = np.zeros( (len(N_quimb), 2) )

for i, n in enumerate(tqdm(N_quimb)):
    times = []
    #GHZ initial state
    state = quimb.tensor.tensor_gen.MPS_ghz_state(n) #do not measure time to generate state!
    for _ in range(50):
        start = time()
        
        qc = QuantumCircuit(n)
        qft_circuit_qiskit(qc, n)
        psi0 = MPS_circ(qc, gates, init_state=state, chi=2)
        psi_dense = psi0.to_dense()
        
        times.append( time()-start )
    
    times_quimb_approx[i, :] = ([np.mean(times), np.std(times) ])

```

### Quimb - full

```{python}
#run QFT on different number of qubits, store execution times
#quimb implementation with maximum choice of bond dimension
times_quimb_full = np.zeros( (len(N_quimb), 2) )

for i, n in enumerate(tqdm(N_quimb)):
    times = []
    #GHZ initial state
    state = quimb.tensor.tensor_gen.MPS_ghz_state(n) 
    for _ in range(50):
        start = time()
        
        #choose the maximum bond dimension possible
        chi = 2**(np.floor(n/2)) 
        
        qc = QuantumCircuit(n)
        qft_circuit_qiskit(qc, n)
        psi0 = MPS_circ(qc, gates, init_state=state, chi=chi)
        psi_dense = psi0.to_dense()
        
        times.append( time()-start )
    
    times_quimb_full[i, :] = ([np.mean(times), np.std(times) ])

```

### Manual - with bond dimension $\chi=2$

```{python}
#run QFT on different number of qubits, store execution times
#manual implementation with bond dimension chi fixed to 2

times_manual_chi = np.zeros( (len(N_quimb), 2) )

for i,n in enumerate(tqdm(N_quimb)):
    times = []
    
    #GHZ initial state
    state = quimb.tensor.tensor_gen.MPS_ghz_state(n) 
    mps = [s.data for s in state]
    state = [mps[0]] + [s.transpose(0,2,1) for s in mps[1:-1]] + [mps[-1]]
    
    for _ in range(50):
        start = time()
    
        result = qft_circuit_swap_approx(state, n, chi=2)
        result = to_dense(result).flatten()
        
        times.append( time()-start )
        
    times_manual_chi[i,:] = ( [np.mean(times), np.std(times) ]) 
    

```

### Manual - full

```{python}
#run QFT on different number of qubits, store execution times
#manual implementation with maximum choice of bond dimension
N_manual = np.arange(4, 17, dtype=int)

times_manual_full = np.zeros( (len(N_manual), 2) )

for i,n in enumerate(tqdm(N_manual)):
    times = []
    
    #GHZ initial state
    state = quimb.tensor.tensor_gen.MPS_ghz_state(n) 
    mps = [s.data for s in state]
    state = [mps[0]] + [s.transpose(0,2,1) for s in mps[1:-1]] + [mps[-1]]
    
    for _ in range(50):
        start = time()
        
        result = qft_circuit_swap_full(state, n)
        result = to_dense(result).flatten()
        
        times.append( time()-start )
        
    times_manual_full[i,:] = ( [np.mean(times), np.std(times) ]) 
    
```

#### Save/load data

```{python}
save = True
if save:
    np.save("times_quimb_approx", times_quimb_approx)
    np.save("times_quimb_full", times_quimb_full)
    np.save("times_manual_chi", times_manual_chi)
    np.save("times_manual_full", times_manual_full)

# times_quimb_approx = np.load(".\\data\\times_quimb_approx.npy")
# times_quimb_full = np.load(".\\data\\times_quimb_full.npy")
# times_manual_chi = np.load(".\\data\\times_manual_chi.npy")
# times_manual_full = np.load(".\\data\\times_manual_full.npy")

# N_quimb = np.arange(4, 17, dtype=int)
# N_manual = np.arange(4, 17, dtype=int)

```

### Fit

```{python}
#Data to be fit is in xs and ys
from scipy.optimize import curve_fit
from scipy import stats

def straight_line(x, a, b):
    return a + b*x

fit_log_xs = np.log10(N_quimb) #Move to log-log space

#mask
mask = fit_log_xs > np.log10(6)
fit_log_xs = fit_log_xs[mask]

fit_log_ys_q_approx = np.log10(times_quimb_approx[mask,0])
err_log_q_approx = np.log10(times_quimb_approx[mask,1])

fit_log_ys_q = np.log10(times_quimb_full[mask,0])
err_log_q = np.log10(times_quimb_full[mask,1])

fit_log_ys_man = np.log10(times_manual_chi[mask,0])
err_log_man = np.log10(times_manual_chi[mask,1])

popt_q_approx, pcov_q_approx = curve_fit(straight_line, fit_log_xs, fit_log_ys_q_approx, sigma=err_log_q_approx)
popt_q, pcov_q = curve_fit(straight_line, fit_log_xs, fit_log_ys_q, sigma= err_log_q)
popt_man, pcov_man = curve_fit(straight_line, fit_log_xs, fit_log_ys_man, sigma = err_log_man)

#b
b_q_approx = popt_q_approx[-1]
b_q = popt_q[-1]
b_man = popt_man[-1]

#errors on b
err_b_approx = np.sqrt(np.diag(pcov_q_approx))[-1]
err_b_q = np.sqrt(np.diag(pcov_q))[-1]
err_b_man = np.sqrt(np.diag(pcov_man))[-1]

x_fit = np.linspace(N_quimb[0], N_quimb[-1], 100) #put the plot range here

```

```{python}
fig, ax = plt.subplots(1, 2, figsize=(6.68, 2.35)) #the first parameter is the width fraction. .45 = half-page (spans one column), .95 = full-page (spans both columns) (accounting for a .05 margin)

# create a color palette
palette = plt.get_cmap('Set2')

ax[0].set_ylim(1e-4, 5)
ax[1].set_ylim(1e-4, 5)



#Comparison of "approximated" computations of QFT
ax[0].plot(N_quimb, times_quimb_approx[:, 0], color=palette(1), linestyle=':', label='quimb')
ax[0].fill_between(N_quimb, times_quimb_approx[:,0]-times_quimb_approx[:,1], 
                   times_quimb_approx[:,0]+times_quimb_approx[:,1],
                   color=palette(1), alpha=0.5)

ax[0].plot(N_quimb, times_manual_chi[:, 0], color=palette(0), linestyle='-.', label='manual')
ax[0].fill_between(N_quimb, times_manual_chi[:,0]-times_manual_chi[:,1], 
                   times_manual_chi[:,0]+times_manual_chi[:,1],
                   color=palette(0), alpha=0.4)
#Fit
ax[0].plot(x_fit, 10**(straight_line(np.log10(x_fit), *popt_q_approx)), '--r', 
           linewidth=0.8, label=f'$b=({np.round(b_q_approx, 1)}\\pm{np.round(err_b_approx, 2)})$') 
ax[0].plot(x_fit, 10**(straight_line(np.log10(x_fit), *popt_man)), '--g', 
           linewidth=0.8, label=f'$b=({np.round(b_man, 1)}\\pm{np.round(err_b_man, 2)})$') 

ax[0].set_ylabel('Time [s]')
ax[0].set_xlabel('Number of qubits')
ax[0].set_title('Time scaling for QFT - $\chi=2$')
ax[0].legend(loc='best')
ax[0].loglog()

#Comparison of "full" computations of QFT
ax[1].plot(N_quimb, times_quimb_full[:, 0], color=palette(1), linestyle=':', label='quimb')
ax[1].fill_between(N_quimb, times_quimb_full[:,0]-times_quimb_full[:,1], 
                   times_quimb_full[:,0]+times_quimb_full[:,1],
                   color=palette(1), alpha=0.5)

ax[1].plot(N_manual, times_manual_full[:, 0], color=palette(0), linestyle='-.', label='manual')
ax[1].fill_between(N_manual, times_manual_full[:,0]-times_manual_full[:,1], 
                   times_manual_full[:,0]+times_manual_full[:,1],
                   color=palette(0), alpha=0.4)

#Fit
ax[1].plot(x_fit, 10**(straight_line(np.log10(x_fit), *popt_q)), '--r', 
           linewidth=0.8, label=f'$b=({np.round(b_q, 1)}\\pm{np.round(err_b_q, 2)})$') 


#ax[1].set_ylabel('Time [s]')
ax[1].set_xlabel('Number of qubits')
ax[1].set_title('Time scaling for QFT - $\chi=2^{\\left \\lfloor n/2 \\right \\rfloor}$')
ax[1].legend(loc='lower right')
ax[1].loglog()

ticks = [4,5,6,7,8,9,10,16]
ax[0].set_xticks(ticks)
ax[1].set_xticks(ticks)
ax[0].set_xticklabels(['${}$'.format(l) for l in ticks])
ax[1].set_xticklabels(['${}$'.format(l) for l in ticks])

ax[0].get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax[1].get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())

ax[0].set_xlim(5,16)
ax[1].set_xlim(5,16)

plt.show()
fig.tight_layout()
plt.savefig("time_scaling.pgf") 
```

## Testing bond dimension $\chi$


### Different $\chi$ for manual method

```{python}
chis = np.array([2**i for i in range(2,6)])

Ns = np.arange(4, 17, dtype=int)

time_chis = np.zeros( (len(chis), len(Ns), 2) )

for i, x in enumerate(tqdm(chis)):
    for j, N in enumerate(Ns):
        times = []
        
        #GHZ initial state
        state = quimb.tensor.tensor_gen.MPS_ghz_state(N) 
        mps = [s.data for s in state]
        state = [mps[0]] + [s.transpose(0,2,1) for s in mps[1:-1]] + [mps[-1]]
    
        for _ in range(50):
            start = time()

            result = qft_circuit_swap_approx(state, N, chi=x)
            #result = to_dense(result).flatten()
            times.append( time()-start )

        time_chis[i, j, :] = ([np.mean(times), np.std(times) ])

```

### Different $\chi$ for quimb implementation

```{python}
time_chis_quimb = np.zeros( (len(chis), len(Ns), 2) )

for i, x in enumerate(tqdm(chis)):
    for j, n in enumerate(Ns):
        times = []
        #GHZ initial state
        state = quimb.tensor.tensor_gen.MPS_ghz_state(n) 
        for _ in range(50):
            start = time()
        
            qc = QuantumCircuit(n)
            qft_circuit_qiskit(qc, n)
            psi0 = MPS_circ(qc, gates, init_state=state, chi=x)
            #psi_dense = psi0.to_dense()

            times.append( time()-start )

        time_chis_quimb[i, j, :] = ([np.mean(times), np.std(times) ])
```

#### Save/load data

```{python}
np.save(".\\data\\time_chis_quimb", time_chis_quimb)
np.save(".\\data\\time_chis", time_chis)

time_chis = np.load(".\\data\\time_chis.npy")
time_chis_quimb = np.load(".\\data\\time_chis_quimb.npy")

chis = np.array([2**i for i in range(2,6)])

Ns = np.arange(4, 17, dtype=int)

N_list = np.array( [8, 16] )

idx_N = N_list - 4

chi_x = np.array([2**i for i in range(1,6)])
```

```{python}
y_man = time_chis.transpose(1,0,2)
y_q   = time_chis_quimb.transpose(1,0,2)

y_man = np.append(times_manual_chi.reshape(1,13,2).transpose(1,0,2), y_man, axis=1)
y_q = np.append(times_quimb_approx.reshape(1,13,2).transpose(1,0,2), y_q, axis=1)
```

```{python}
fig, ax = plt.subplots(1, 2, figsize=(6.68, 2.20)) #the first parameter is the width fraction. .45 = half-page (spans one column), .95 = full-page (spans both columns) (accounting for a .05 margin)

#Manual implementation

for i in range(len(N_list)):
    ax[0].plot(chi_x, y_man[idx_N[i], :, 0], label='$n='+str(N_list[i])+'$', 
               color=palette(i+2))#, alpha=0.2*(len(chis)-i))
    ax[0].fill_between(chi_x, y_man[idx_N[i], :, 0]-y_man[idx_N[i], :, 1], 
                       y_man[idx_N[i], :, 0]+y_man[idx_N[i], :, 1],
                       color=palette(i+2), alpha=0.1*(len(chis)-i))

ax[0].set_ylabel('Time [s]')
ax[0].set_xlabel('Bond dimension $\chi$')
ax[0].set_title('Manual QFT circuit')
ax[0].legend(loc='best')
ax[0].loglog()

#Quimb implementation

for i in range(len(N_list)):
    ax[1].plot(chi_x, y_q[idx_N[i], :, 0], label='$n='+str(N_list[i])+'$',
               color=palette(i+2))
    ax[1].fill_between(chi_x, y_q[idx_N[i], :, 0]-y_q[idx_N[i], :, 1], 
                   y_q[idx_N[i], :, 0]+y_q[idx_N[i], :, 1],
                   color=palette(i+2), alpha=0.1*(len(chis)-i))

ax[1].set_xlabel('Bond dimension $\chi$')
ax[1].set_title('quimb QFT circuit')
ax[1].legend(loc='best', )
ax[1].loglog()

plt.show()
fig.tight_layout()
plt.savefig("chi_scaling.pgf") 

```

## Testing maximum number of Qubits


### Manual implementation - $\chi=2$

```{python}
#run QFT on different number of qubits, store execution times
#manual implementation with bond dimension chi fixed to 2
N_man_tot = np.arange(4, 26, dtype=int)

times_manual_tot = np.zeros( (len(N_man_tot), 2) )

for i,n in enumerate(tqdm(N_man_tot)):
    times = []
    
    #GHZ initial state
    state = quimb.tensor.tensor_gen.MPS_ghz_state(n) 
    mps = [s.data for s in state]
    state = [mps[0]] + [s.transpose(0,2,1) for s in mps[1:-1]] + [mps[-1]]
        
    for _ in range(10):
        start = time()

        result = qft_circuit_swap_approx(state, n, chi=2)
        #result = to_dense(result).flatten()
        
        times.append( time()-start )
        
    times_manual_tot[i,:] = ( [np.mean(times), np.std(times) ]) 
```

### Quimb - authomatic $\chi$

```{python}
#run QFT on different number of qubits, store execution times
#quimb implementation with authomatic choice of bond dimension
N_tot =  np.arange(4, 51, dtype=int)

times_quimb_auto = np.zeros( (len(N_tot), 2) )

for i, n in enumerate(tqdm(N_tot)):
    times = []
    
    #GHZ initial state
    state = quimb.tensor.tensor_gen.MPS_ghz_state(n) 
    
    for _ in range(10):
        start = time()    
        
        qc = QuantumCircuit(n)
        qft_circuit_qiskit(qc, n)
        psi0 = MPS_circ(qc, gates, init_state=state)
        #psi_dense = psi0.to_dense()
        
        times.append( time()-start )
    
    times_quimb_auto[i, :] = ([np.mean(times), np.std(times) ])

```

### Qiskit

```{python}
#run QFT on different number of qubits, store execution times
#qiskit implementation 
N_qiskit = np.arange(4, 26, dtype=int)

times_qiskit = np.zeros( (len(N_qiskit), 2) )

for i, n in enumerate(tqdm(N_qiskit)):
    times = []
    
    #GHZ initial state
    state = get_ghz_simple(n, measure=False)
    for _ in range(10):
        start = time()
        
        #QFT circuit
        qc = QuantumCircuit(n)
        qft_circuit_qiskit(qc, n)
        #combine the two circuits
        qc = state + qc
        # Select the QasmSimulator from the Aer provider
        simulator = Aer.get_backend('qasm_simulator')
        # Execute 
        result = execute(qc, simulator, shots=1).result()
        
        times.append( time()-start )
    
    times_qiskit[i, :] = ([np.mean(times), np.std(times) ])

```

### Cirq

```{python}
def cphase_swap_cirq(ctrl, target, phase):
    yield cirq.CZ(ctrl, target) ** phase
    yield cirq.SWAP(ctrl, target)

def qft_circuit_cirq(qubits, circuit=[]):
    """
    Build a circuit implementing the QFT algorithm on the given @qubits. 
    The order of @qubits is preserved by SWAP operations.
    Implemented using only local operations, i.e. gates acting on neighbouring qubits.
    Adapted from: https://github.com/quantumlib/Cirq/blob/master/examples/quantum_fourier_transform.py and extended to
    n generic qubits through recursion.
    """
    n = len(qubits)
    assert n > 0, "Number of qubits must be > 0"
    
    if (n == 1):
        circuit.append(cirq.H(qubits[0]))
        return cirq.Circuit(circuit, strategy=cirq.InsertStrategy.EARLIEST)
    else:
        circuit.append(cirq.H(qubits[0]))
        circuit.extend(cphase_swap_cirq(qubits[i], qubits[i+1], 1/2**(i+1)) for i in range(n-1))
        return qft_circuit_cirq(qubits[:n-1], circuit)

```

```{python}
#run QFT on different number of qubits, store execution times
#qiskit implementation 
N_cirq = np.arange(4, 26, dtype=int)

times_cirq = np.zeros( (len(N_cirq), 2) )

for i, n in enumerate(tqdm(N_cirq)):
    times = []
    
    #GHZ initial state
    state = np.zeros(2**n)
    state[0] = 1
    state[-1] = 1
    state = state / np.sqrt(2)
        
    for _ in range(10):
        start = time()
        
        #QFT circuit
        qubits = cirq.LineQubit.range(n)
        qc = qft_circuit_cirq(qubits, [])

        simulator = cirq.Simulator()
        result = simulator.simulate(qc, initial_state=state)

        times.append( time()-start )
    
    times_cirq[i, :] = ([np.mean(times), np.std(times) ])

```

#### Save/load data

```{python}
# np.save("times_quimb_auto", times_quimb_auto)
# np.save("times_manual_tot", times_manual_tot)
# np.save("times_cirq", times_cirq)
# np.save("times_qiskit", times_qiskit)

times_quimb_auto = np.load(".\\data\\times_quimb_auto.npy")
times_manual_tot = np.load(".\\data\\times_manual_tot.npy")
times_cirq = np.load(".\\data\\times_cirq.npy")
times_qiskit = np.load(".\\data\\times_qiskit.npy")

N_tot =  np.arange(4, 51, dtype=int)
N_man_tot = np.arange(4, 26, dtype=int)
N_cirq = np.arange(4, 26, dtype=int)
N_qiskit = np.arange(4, 26, dtype=int)

```

```{python}
#Data to be fit is in xs and ys
from scipy.optimize import curve_fit
from scipy import stats

def straight_line(x, a, b):
    return a + b*x

fit_log_xs = np.log10(N_tot) #Move to log-log space

#mask
mask = fit_log_xs > np.log10(6)
fit_log_xs = fit_log_xs[mask]

fit_q_tot = np.log10(times_quimb_auto[mask,0])

popt_q_tot, pcov_tot = curve_fit(straight_line, fit_log_xs, fit_q_tot)

#b
b_q_tot = popt_q_tot[-1]
err_b_tot = np.sqrt(np.diag(pcov_tot))[-1]

x_fit = np.linspace(N_tot[0], N_tot[-1], 100) #put the plot range here

```

```{python}
fig, ax = plt.subplots(figsize=get_figsize(0.95)) #the first parameter is the width fraction. .45 = half-page (spans one column), .95 = full-page (spans both columns) (accounting for a .05 margin)

#quimb implementation
plt.plot(N_tot, times_quimb_auto[:, 0], color=palette(1), linestyle=':', label='quimb')
plt.fill_between(N_tot, times_quimb_auto[:,0]-times_quimb_auto[:,1], 
                   times_quimb_auto[:,0]+times_quimb_auto[:,1],
                   color=palette(1), alpha=0.6)

#manual implementation
plt.plot(N_man_tot, times_manual_tot[:, 0], color=palette(0), linestyle='-.', label='manual, $\chi=2$')
plt.fill_between(N_man_tot, times_manual_tot[:,0]-times_manual_tot[:,1], 
                   times_manual_tot[:,0]+times_manual_tot[:,1],
                   color=palette(0), alpha=0.6)

#cirq
plt.plot(N_cirq, times_cirq[:, 0], color='plum', label='cirq', alpha=0.8)
plt.fill_between(N_cirq, times_cirq[:,0]-times_cirq[:,1], 
                   times_cirq[:,0]+times_cirq[:,1],
                   color='thistle', alpha=0.6)

#qiskit
plt.plot(N_qiskit, times_qiskit[:, 0], color='cadetblue', label='qiskit', alpha=0.8)
plt.fill_between(N_qiskit, times_qiskit[:,0]-times_qiskit[:,1], 
                   times_qiskit[:,0]+times_qiskit[:,1],
                   color='paleturquoise', alpha=0.6)

#Fit
plt.plot(x_fit, 10**(straight_line(np.log10(x_fit), *popt_q_tot)), '--r', 
           linewidth=0.8, label=f'$b=({np.round(b_q_tot, 2)}\\pm{np.round(err_b_tot, 2)})$') 


plt.ylabel('Time [s]')
plt.xlabel('Number of qubits')
plt.title('Time scaling for QFT circuit')
plt.legend(loc='best')
plt.loglog()

plt.show()

fig.tight_layout()
plt.savefig("max_n_qubits.pgf") #Save as pdf in the "report_plots" folder

```

```{python}

```

```{python}

```