from matplotlib import pyplot
from qiskit import *
from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from qiskit import IBMQ, Aer, execute,assemble,QuantumCircuit, aqua
from qiskit.visualization import plot_histogram, plot_bloch_vector, plot_bloch_multivector
from qiskit.quantum_info import Statevector
from qiskit.extensions import *
provider = IBMQ.load_account()

from qiskit.quantum_info import random_unitary
from qiskit.providers.aer import QasmSimulator, StatevectorSimulator, UnitarySimulator
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import math
from math import pi, sqrt

from scipy.special import rel_entr
from random import seed
from random import random
import cmath


#Possible Bin
bins_list=[];
for i in range(76):
    bins_list.append((i)/75)
#Center of the Bean
bins_x=[]    
for i in range(75):
    bins_x.append(bins_list[1]+bins_list[i])
def P_harr(l,u,N):
    return (1-l)**(N-1)-(1-u)**(N-1)
#Harr historgram
P_harr_hist=[]
for i in range(75):
    P_harr_hist.append(P_harr(bins_list[i],bins_list[i+1],2))    
#Imaginary    
j=(-1)**(1/2)

backend = Aer.get_backend('qasm_simulator')
for x in range(1,20,1):
    print()
arr = []
for nsh in range(1,20,1):
    arr.append([])
    for lp in range(1,10,1):
        nshot=int(round(1000*(nsh)**1.5,0))
        nparam=1000*lp
        fidelity=[]    
        for x in range(nparam):
            qr = QuantumRegister(1)
            cr = ClassicalRegister(1)
            qc = QuantumCircuit(qr, cr)

            u13=UnitaryGate(random_unitary(2))
            qc.append(u13, [qr[0]] )
            u13=UnitaryGate(random_unitary(2))
            qc.append(u13, [qr[0]] )

            qc.measure(qr[0],cr[0])
            job = execute(qc, backend, shots=nshot)
            result = job.result()
            count =result.get_counts()

            if '0' in count and '1' in count:
                ratio=count['0']/nshot
            elif '0' in count and '1' not in count:
                ratio=count['0']/nshot
            else:
                ratio=0

            fidelity.append(ratio)   
            #Kullback Leibler divergence
        weights = np.ones_like(fidelity)/float(len(fidelity))            
        P_U_hist=np.histogram(fidelity, bins=bins_list, weights=weights, range=[0, 1])[0];
        #Kullback Leibler divergence
        print(nshot,'shots per simulation',nparam,'distribution size')
        arr[nsh-1].append(sum(rel_entr(P_U_hist, P_harr_hist)))


import sys
import numpy
numpy.set_printoptions(threshold=sys.maxsize)

np.savetxt("KLDivg_qasm.txt",qasm)

qasm

backend = QasmSimulator(method='density_matrix')

for x in range(1,20,1):
    print()
arr = []
for nsh in range(1,20,1):
    arr.append([])
    for lp in range(1,10,1):
        nshot=int(round(1000*(nsh)**1.5,0))
        nparam=1000*lp
        fidelity=[]    
        for x in range(nparam):
            qr = QuantumRegister(1)
            cr = ClassicalRegister(1)
            qc = QuantumCircuit(qr, cr)

            u13=UnitaryGate(random_unitary(2))
            qc.append(u13, [qr[0]] )
            u13=UnitaryGate(random_unitary(2))
            qc.append(u13, [qr[0]] )

            qc.measure(qr[0],cr[0])
            job = execute(qc, backend, shots=nshot)
            result = job.result()
            count =result.get_counts()

            if '0' in count and '1' in count:
                ratio=count['0']/nshot
            elif '0' in count and '1' not in count:
                ratio=count['0']/nshot
            else:
                ratio=0

            fidelity.append(ratio)   
            #Kullback Leibler divergence
        weights = np.ones_like(fidelity)/float(len(fidelity))            
        P_U_hist=np.histogram(fidelity, bins=bins_list, weights=weights, range=[0, 1])[0];
        #Kullback Leibler divergence
        
        arr[nsh-1].append(sum(rel_entr(P_U_hist, P_harr_hist)))
        print(nshot,'shots per simulation',nparam,'distribution size KL=',sum(rel_entr(P_U_hist, P_harr_hist)))
import sys
import numpy
density_mat=arr
numpy.set_printoptions(threshold=sys.maxsize)
np.savetxt("KLDivg_density_mat.txt",density_mat)

density_mat

backend = QasmSimulator(method='statevector')

for x in range(1,20,1):
    print()
arr = []
for nsh in range(1,20,1):
    arr.append([])
    for lp in range(1,10,1):
        nshot=int(round(1000*(nsh)**1.5,0))
        nparam=1000*lp
        fidelity=[]    
        for x in range(nparam):
            qr = QuantumRegister(1)
            cr = ClassicalRegister(1)
            qc = QuantumCircuit(qr, cr)

            u13=UnitaryGate(random_unitary(2))
            qc.append(u13, [qr[0]] )
            u13=UnitaryGate(random_unitary(2))
            qc.append(u13, [qr[0]] )

            qc.measure(qr[0],cr[0])
            job = execute(qc, backend, shots=nshot)
            result = job.result()
            count =result.get_counts()

            if '0' in count and '1' in count:
                ratio=count['0']/nshot
            elif '0' in count and '1' not in count:
                ratio=count['0']/nshot
            else:
                ratio=0

            fidelity.append(ratio)   
            #Kullback Leibler divergence
        weights = np.ones_like(fidelity)/float(len(fidelity))            
        P_U_hist=np.histogram(fidelity, bins=bins_list, weights=weights, range=[0, 1])[0];
        #Kullback Leibler divergence
        
        arr[nsh-1].append(sum(rel_entr(P_U_hist, P_harr_hist)))
        print(nshot,'shots per simulation',nparam,'distribution size KL=',sum(rel_entr(P_U_hist, P_harr_hist)))
import sys
import numpy
statevector1=arr
numpy.set_printoptions(threshold=sys.maxsize)
np.savetxt("KLDivg_statevector.txt",statevector1)

statevector1

backend = QasmSimulator(method='matrix_product_state')

for x in range(1,20,1):
    print()
arr = []
for nsh in range(1,20,1):
    arr.append([])
    for lp in range(1,10,1):
        nshot=int(round(1000*(nsh)**1.5,0))
        nparam=1000*lp
        fidelity=[]    
        for x in range(nparam):
            qr = QuantumRegister(1)
            cr = ClassicalRegister(1)
            qc = QuantumCircuit(qr, cr)

            u13=UnitaryGate(random_unitary(2))
            qc.append(u13, [qr[0]] )
            u13=UnitaryGate(random_unitary(2))
            qc.append(u13, [qr[0]] )

            qc.measure(qr[0],cr[0])
            job = execute(qc, backend, shots=nshot)
            result = job.result()
            count =result.get_counts()

            if '0' in count and '1' in count:
                ratio=count['0']/nshot
            elif '0' in count and '1' not in count:
                ratio=count['0']/nshot
            else:
                ratio=0

            fidelity.append(ratio)   
            #Kullback Leibler divergence
        weights = np.ones_like(fidelity)/float(len(fidelity))            
        P_U_hist=np.histogram(fidelity, bins=bins_list, weights=weights, range=[0, 1])[0];
        #Kullback Leibler divergence
        
        arr[nsh-1].append(sum(rel_entr(P_U_hist, P_harr_hist)))
        print(nshot,'shots per simulation',nparam,'distribution size KL=',sum(rel_entr(P_U_hist, P_harr_hist)))
import sys
import numpy
matrix_product1=arr
numpy.set_printoptions(threshold=sys.maxsize)
np.savetxt("KLDivg_matrix_product1.txt",matrix_product1)

matrix_product1

def plotdata(i,data):
    klll=np.transpose(data)
    x=[];
    y=[];
    for nsh in range(0,9,1):
        x.append(int(round(1000*(nsh)**1.5,0)))
        y.append(klll[nsh][i])
    return [x,y]

#loading dataset
fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(20, 5))
datalist=[density_mat,statevector1,matrix_product1,qasm]
datalistn=['density_mat','statevector','matrix_product','qasm']
indx=0;
for indx in range(4):
    data=datalist[indx];
    axes[indx].plot(plotdata(0,data)[0],plotdata(0,data)[1], color='red', label='npram=1000');
    axes[indx].plot(plotdata(1,data)[0],plotdata(1,data)[1], color='blue', label='npram=2000');
    axes[indx].plot(plotdata(2,data)[0],plotdata(2,data)[1], color='black', label='npram=3000');
    axes[indx].plot(plotdata(3,data)[0],plotdata(3,data)[1], color='green', label='npram=4000');
    axes[indx].plot(plotdata(4,data)[0],plotdata(4,data)[1], color='purple',  label='npram=5000');
    axes[indx].plot(plotdata(5,data)[0],plotdata(5,data)[1], color='gray',  label='npram=6000');
    axes[indx].plot(plotdata(6,data)[0],plotdata(6,data)[1], color='black',  label='npram=7000');
    axes[indx].plot(plotdata(7,data)[0],plotdata(7,data)[1], color='yellow',  label='npram=8000');
    axes[indx].plot(plotdata(8,data)[0],plotdata(8,data)[1], color='pink',  label='npram=9000');

    axes[indx].set_ylim([0.003, 0.05])
    axes[indx].legend(loc='upper right')
    axes[indx].set_title(datalistn[indx])
    axes[indx].set_yscale('log',base=10) 
    axes[indx].set_ylabel('Expressibility')

from matplotlib import pyplot as plt
fig.savefig('ExpressibilitybySimulator.png')

fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(20, 5))
indx=0
for indx in range(4):
    data=density_mat;
    axes[indx].plot(plotdata(indx,data)[0],plotdata(indx,data)[1], color='red', label='density_mat');
    data=statevector1;
    axes[indx].plot(plotdata(indx,data)[0],plotdata(indx,data)[1], color='blue', label='statevector');
    data=matrix_product1;
    axes[indx].plot(plotdata(indx,data)[0],plotdata(indx,data)[1], color='green', label='matrix_product1');
    data=qasm;
    axes[indx].plot(plotdata(indx,data)[0],plotdata(indx,data)[1], color='pink', label='qasm');
    axes[indx].set_title('nparam='+str((indx+1)*1000))
    axes[indx].set_yscale('log',base=10) 
    axes[indx].set_ylim([0.003, 0.05])
    axes[indx].set_ylabel('Expressibility')
    axes[indx].legend(['density_mat','statevector','matrix_product1','qasm'])
    axes[indx].set_xlabel('iteration')
fig.savefig('ExpressibilitybySimulatornparam1.png')

fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(20, 5))
indx=0
for indx1 in range(4,8,1):
    indx=indx1-4
    data=density_mat;
    axes[indx].plot(plotdata(indx1,data)[0],plotdata(indx1,data)[1], color='red', label='density_mat');
    data=statevector1;
    axes[indx].plot(plotdata(indx1,data)[0],plotdata(indx1,data)[1], color='blue', label='statevector');
    data=matrix_product1;
    axes[indx].plot(plotdata(indx1,data)[0],plotdata(indx1,data)[1], color='green', label='matrix_product1');
    data=qasm;
    axes[indx].plot(plotdata(indx1,data)[0],plotdata(indx1,data)[1], color='pink', label='qasm');
    axes[indx].set_title('nparam='+str((indx1+1)*1000))
    axes[indx].set_yscale('log',base=10) 
    axes[indx].set_ylabel('Expressibility')
    axes[indx].set_ylim([0.003, 0.05])
    axes[indx].legend(['density_mat','statevector','matrix_product1','qasm'])
    axes[indx].set_xlabel('iteration')
fig.savefig('ExpressibilitybySimulatornparam2.png')

