import numpy as np
from scipy.sparse import identity
from scipy import sparse
from qiskit import*

L = 4

gates_list = open('gates_list_'+str(L)+'.txt','r')

MCX_transpile = []
for gates in gates_list:
    MCX_transpile.append(gates.split(","))

len(MCX_transpile)

number_of_gates = len(MCX_transpile)
Seed = 4000
np.random.seed(Seed)
Noise = np.random.rand(number_of_gates)

def H_fixed(): # Hadamad gate acting on one qubit.
    
    return 1/np.sqrt(2)*np.matrix([[1,1],
                                   [1,-1]])

def Ry(theta):
    return np.matrix([[np.cos(theta/2), -np.sin(theta/2)],
                      [np.sin(theta/2), np.cos(theta/2)]])

def sigma_Z():
    return np.matrix([[1,0],
                      [0,-1]])

def sigma_X():
    return np.matrix([[0,1],
                      [1,0]])

def Rz(theta):
    return np.matrix([[np.exp(-1j*theta/2), 0],
                      [0, np.exp(1j*theta/2)]])

'''
For noise = 0; the following returns the fixed Hadamard gate.
'''
def Hadamard_variable(noise):
    return Ry(np.pi/2+noise)*sigma_Z()

def sigma_X_variable(noise):
    return Hadamard_variable(noise)*sigma_Z()*Hadamard_variable(noise)

def Cx(noise):
    Pi_0 = np.matrix([[1,0],[0,0]])
    final_matrix_1 = sparse.kron(Pi_0,np.identity(2))
    Pi_1 = np.matrix([[0,0],[0,1]])
    final_matrix_2 = sparse.kron(Pi_1,sigma_X_variable(noise))
    return final_matrix_1+final_matrix_2

def CU_gate(quantum_gate,control_qubit,target_qubit):
    Pi_0 = np.matrix([[1,0],[0,0]])
    Pi_1 = np.matrix([[0,0],[0,1]])
    '''
    List below will hold gates acting on one qubit. For example, for L = 8,
    gate U acting on control qubit c and target qubit t is
    I x I x ...x (Pi_0)_c x ...x (I)_t x ... x I + I x I x ...x (Pi_1)_c x ...x (U)_t x ... x I 
    ''' 
    qubits_list_1 = [] 
    for i in range(L):
        
        if i == control_qubit:
            
            qubits_list_1.append(Pi_0)
            
        else: # Other gates are identity operators.
            
            qubits_list_1.append(np.identity(2))
       
    qubits_list_2 = []
    for i in range(L):
        
        if i == control_qubit:
            
            qubits_list_2.append(Pi_1)
            
        elif i == target_qubit:
            
            qubits_list_2.append(quantum_gate)
            
        else: # Other gates are identity operators.
            
            qubits_list_2.append(np.identity(2))                    
    #qubits_list_1.reverse()    
    #qubits_list_2.reverse()    
    final_matrix_1 = sparse.csr_matrix(qubits_list_1[0]) # Initializes the final matrix.
    for g in range(1,len(qubits_list_1)):
        
        final_matrix_1 = sparse.kron(qubits_list_1[g],final_matrix_1) # kronecker product.
        
    final_matrix_2 = sparse.csr_matrix(qubits_list_2[0]) # Initializes the final matrix.
    for g in range(1,len(qubits_list_2)):
        
        final_matrix_2 = sparse.kron(qubits_list_2[g],final_matrix_2) # kronecker product.        
    return final_matrix_1+final_matrix_2   

def N_th_qubit_gate(quantum_gate,qubit_number):

    
    '''
    List below will hold gates acting on one qubit. For example, for L = 3,
    the Hadamard gate acting on the qubit 1 is given by = 1 x H x 1, where 
    x is the Kronecker product. Then, qubits_list = [1,H,1].

    ''' 
    qubits_list = [] 
    
    for i in range(L):
        
        if i == qubit_number: # qubit_number^th position in the list is the quantum_gate.
            
            qubits_list.append(quantum_gate)
            
        else: # Other gates are identity operators.
            
            qubits_list.append(np.identity(2))
            
    '''
    The following loop performs the Kronecker product.
    '''        
    
    final_matrix = sparse.csr_matrix(qubits_list[0]) # Initializes the final matrix.
    
    for i in range(1,len(qubits_list)):
        
        final_matrix = sparse.kron(final_matrix,qubits_list[i]) # kronecker product.
        
    return final_matrix    

U_0 = np.identity(2**L, dtype = complex);
Delta = 0.0
noise_counter = 0
for i in MCX_transpile:
    if i[0] == 'h':
        epsilon = Delta*Noise[noise_counter]
        noise_counter += 1
        qubit = int(i[2])
        U_0 = U_0*N_th_qubit_gate(Hadamard_variable(epsilon),qubit)
        #np.matmul(U_0,N_th_qubit_gate(Hadamard_variable(epsilon),qubit))
    elif i[0] == 'cx':
        epsilon = Delta*Noise[noise_counter]
        noise_counter += 1
        control_qubit = int(i[1])
        target_qubit = int(i[2])
        U_0 = U_0*CU_gate(sigma_X_variable(epsilon),control_qubit,target_qubit)
        #np.matmul(U_0,CU_gate(sigma_X_variable(epsilon),control_qubit,target_qubit)) 
    elif i[0] == 'rz':
        epsilon = Delta*Noise[noise_counter]
        noise_counter += 1
        qubit = int(i[2])
        angle = float(i[1])
        U_0 = U_0*N_th_qubit_gate(Rz(epsilon+angle),qubit)
m = np.around(U_0,2)
m = m/m[0,0]
m.real

def Grover_operator(Delta):
    noise_counter = 0
    '''
    First the U_x will be created.
    '''
    U_x = np.identity(2**L, dtype = complex);
    for i in range(L):
        epsilon = Delta*Noise[noise_counter]
        noise_counter += 1
        U_x = U_x*N_th_qubit_gate(Hadamard_variable(epsilon),i)
        #np.matmul(U_x,N_th_qubit_gate(Hadamard_variable(epsilon),i))
    for i in range(L):
        epsilon = Delta*Noise[noise_counter]
        noise_counter += 1
        U_x = U_x*N_th_qubit_gate(sigma_X_variable(epsilon),i)
        #np.matmul(U_x,N_th_qubit_gate(sigma_X_variable(epsilon),i))
    
    # mcx gate of U_x
    for i in MCX_transpile:
        if i[0] == 'h':
            epsilon = Delta*Noise[noise_counter]
            noise_counter += 1
            qubit = int(i[2])
            U_x = U_x*N_th_qubit_gate(Hadamard_variable(epsilon),qubit)
            #np.matmul(U_x,N_th_qubit_gate(Hadamard_variable(epsilon),qubit))
        elif i[0] == 'cx':
            epsilon = Delta*Noise[noise_counter]
            noise_counter += 1
            control_qubit = int(i[1])
            target_qubit = int(i[2])
            U_x = U_x*CU_gate(sigma_X_variable(epsilon),control_qubit,target_qubit)
            #np.matmul(U_x,CU_gate(sigma_X_variable(epsilon),control_qubit,target_qubit)) 
        elif i[0] == 'rz':
            epsilon = Delta*Noise[noise_counter]
            noise_counter += 1
            qubit = int(i[2])
            angle = float(i[1])
            U_x = U_x*N_th_qubit_gate(Rz(epsilon+angle),qubit)
            #np.matmul(U_x,N_th_qubit_gate(Rz_variable(epsilon+angle),qubit)) 
            
        for i in range(L):
            epsilon = Delta*Noise[noise_counter]
            noise_counter += 1
            U_x = U_x*N_th_qubit_gate(sigma_X_variable(epsilon),i)
            #np.matmul(U_x,N_th_qubit_gate(sigma_X_variable(epsilon),i))     
        for i in range(L):
            epsilon = Delta*Noise[noise_counter]
            noise_counter += 1
            U_x = U_x*N_th_qubit_gate(Hadamard_variable(epsilon),i)
            #np.matmul(U_x,N_th_qubit_gate(Hadamard_variable(epsilon),i))
     
    '''
    First the U_0 will be created.
    '''       
    U_0 = np.identity(2**L, dtype = complex);
    for i in range(L):
        epsilon = Delta*Noise[noise_counter]
        noise_counter += 1
        U_0 = U_0*N_th_qubit_gate(sigma_X_variable(epsilon),i)
        #np.matmul(U_0,N_th_qubit_gate(sigma_X_variable(epsilon),i))
        
    epsilon = Delta*Noise[noise_counter]
    noise_counter += 1
    U_0 = U_0*N_th_qubit_gate(Hadamard_variable(epsilon),L)
    #np.matmul(U_0,N_th_qubit_gate(Hadamard_variable(epsilon),L))
    
    # mcx gate of U_0
    for i in MCX_transpile:
        if i[0] == 'h':
            epsilon = Delta*Noise[noise_counter]
            noise_counter += 1
            qubit = int(i[2])
            U_0 = U_0*N_th_qubit_gate(Hadamard_variable(epsilon),qubit)
            #np.matmul(U_0,N_th_qubit_gate(Hadamard_variable(epsilon),qubit))
        elif i[0] == 'cx':
            epsilon = Delta*Noise[noise_counter]
            noise_counter += 1
            control_qubit = int(i[1])
            target_qubit = int(i[2])
            U_0 = U_0*CU_gate(sigma_X_variable(epsilon),control_qubit,target_qubit)
            #np.matmul(U_0,CU_gate(sigma_X_variable(epsilon),control_qubit,target_qubit)) 
        elif i[0] == 'rz':
            epsilon = Delta*Noise[noise_counter]
            noise_counter += 1
            qubit = int(i[2])
            angle = float(i[1])
            U_0 = U_0*N_th_qubit_gate(Rz(epsilon+angle),qubit)
            #np.matmul(U_0,N_th_qubit_gate(Rz_variable(epsilon+angle),qubit))     

        for i in range(L):
            epsilon = Delta*Noise[noise_counter]
            noise_counter += 1
            U_0 = U_0*N_th_qubit_gate(sigma_X_variable(epsilon),i)
            #np.matmul(U_0,N_th_qubit_gate(sigma_X_variable(epsilon),i))

        epsilon = Delta*Noise[noise_counter]
        noise_counter += 1
        U_0 = U_0*N_th_qubit_gate(Hadamard_variable(epsilon),L)
        #np.matmul(U_0,N_th_qubit_gate(Hadamard_variable(epsilon),L))  
        
    return U_0