%%capture
%pip install qiskit
%pip install qiskit_ibm_provider
%pip install qiskit-aer

# Importing standard Qiskit libraries
from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit, QuantumCircuit, transpile, Aer
from qiskit_ibm_provider import IBMProvider
from qiskit.tools.jupyter import *
from qiskit.visualization import *
from qiskit.circuit.library import C3XGate

# Importing matplotlib
import matplotlib.pyplot as plt

# Importing Numpy, Cmath and math
import numpy as np
import os, math, cmath
from numpy import pi

# Other imports
from IPython.display import display, Math, Latex

# Specify the path to your env file
env_file_path = 'config.env'

# Load environment variables from the file
os.environ.update(line.strip().split('=', 1) for line in open(env_file_path) if '=' in line and not line.startswith('#'))

# Load IBM Provider API KEY
IBMP_API_KEY = os.environ.get('IBMP_API_KEY')

# Loading your IBM Quantum account(s)
IBMProvider.save_account(IBMP_API_KEY, overwrite=True)

# Run the quantum circuit on a statevector simulator backend
backend = Aer.get_backend('statevector_simulator')

qc_b1 = QuantumCircuit(2, 2)
qc_b1.h(0)
qc_b1.cx(0, 1)
qc_b1.draw(output='mpl', style="iqp")

sv = backend.run(qc_b1).result().get_statevector()
sv.draw(output='latex', prefix = "|\Phi^+\\rangle = ")

qc_b2 = QuantumCircuit(2, 2)
qc_b2.x(0)
qc_b2.h(0)
qc_b2.cx(0, 1)
qc_b2.draw(output='mpl', style="iqp")

sv = backend.run(qc_b2).result().get_statevector()
sv.draw(output='latex', prefix = "|\Phi^-\\rangle = ")

qc_b2 = QuantumCircuit(2, 2)
qc_b2.h(0)
qc_b2.cx(0, 1)
qc_b2.z(0)
qc_b2.draw(output='mpl', style="iqp")

sv = backend.run(qc_b2).result().get_statevector()
sv.draw(output='latex', prefix = "|\Phi^-\\rangle = ")

qc_b3 = QuantumCircuit(2, 2)
qc_b3.x(1)
qc_b3.h(0)
qc_b3.cx(0, 1)
qc_b3.draw(output='mpl', style="iqp")

sv = backend.run(qc_b3).result().get_statevector()
sv.draw(output='latex', prefix = "|\Psi^+\\rangle = ")

qc_b3 = QuantumCircuit(2, 2)
qc_b3.h(0)
qc_b3.cx(0, 1)
qc_b3.x(0)
qc_b3.draw(output='mpl', style="iqp")

sv = backend.run(qc_b3).result().get_statevector()
sv.draw(output='latex', prefix = "|\Psi^+\\rangle = ")

qc_b4 = QuantumCircuit(2, 2)
qc_b4.x(0)
qc_b4.h(0)
qc_b4.x(1)
qc_b4.cx(0, 1)
qc_b4.draw(output='mpl', style="iqp")

sv = backend.run(qc_b4).result().get_statevector()
sv.draw(output='latex', prefix = "|\Psi^-\\rangle = ")

qc_b4 = QuantumCircuit(2, 2)
qc_b4.h(0)
qc_b4.cx(0, 1)
qc_b4.x(0)
qc_b4.z(1)
qc_b4.draw(output='mpl', style="iqp")

sv = backend.run(qc_b4).result().get_statevector()
sv.draw(output='latex', prefix = "|\Psi^-\\rangle = ")

def sv_latex_from_qc(qc, backend):
    sv = backend.run(qc).result().get_statevector()
    return sv.draw(output='latex')


qc_ej2 = QuantumCircuit(4, 4)
qc_ej2.x(0)
sv_latex_from_qc(qc_ej2, backend)

qc_ej2 = QuantumCircuit(4, 4)
qc_ej2.x(1)
sv_latex_from_qc(qc_ej2, backend)

qc_ej2 = QuantumCircuit(4, 4)
qc_ej2.x(0)
qc_ej2.x(1)
sv_latex_from_qc(qc_ej2, backend)

qc_ej2 = QuantumCircuit(4, 4)
qc_ej2.x(2)
sv_latex_from_qc(qc_ej2, backend)

qc_ej2 = QuantumCircuit(4, 4)
qc_ej2.x(0)
qc_ej2.x(2)
sv_latex_from_qc(qc_ej2, backend)

qc_ej2 = QuantumCircuit(4, 4)
qc_ej2.x(0)
qc_ej2.x(2)
sv_latex_from_qc(qc_ej2, backend)

qc_ej2 = QuantumCircuit(4, 4)
qc_ej2.x(0)
qc_ej2.x(1)
qc_ej2.x(2)
sv_latex_from_qc(qc_ej2, backend)

qc_ej2 = QuantumCircuit(4, 4)
qc_ej2.x(3)
sv_latex_from_qc(qc_ej2, backend)

def circuit_adder (num):
    if num<1 or num>8:
        raise ValueError("Out of range")  ## El enunciado limita el sumador a los valores entre 1 y 8. Quitar esta restricción sería directo.
    # Definición del circuito base que vamos a construir
    qreg_q = QuantumRegister(4, 'q')
    creg_c = ClassicalRegister(1, 'c')
    circuit = QuantumCircuit(qreg_q, creg_c)
    
    qbit_position = 0
    for element in reversed(np.binary_repr(num)):
        if (element=='1'):
            circuit.barrier()
            match qbit_position:
                case 0: # +1
                    circuit.append(C3XGate(), [qreg_q[0], qreg_q[1], qreg_q[2], qreg_q[3]])
                    circuit.ccx(qreg_q[0], qreg_q[1], qreg_q[2])
                    circuit.cx(qreg_q[0], qreg_q[1])
                    circuit.x(qreg_q[0])
                case 1: # +2
                    circuit.ccx(qreg_q[1], qreg_q[2], qreg_q[3])
                    circuit.cx(qreg_q[1], qreg_q[2])
                    circuit.x(qreg_q[1])
                case 2: # +4
                    circuit.cx(qreg_q[2], qreg_q[3])
                    circuit.x(qreg_q[2])
                case 3: # +8
                    circuit.x(qreg_q[3])
        qbit_position+=1
    return circuit

add_3 = circuit_adder(3)
add_3.draw(output='mpl', style="iqp")

qc_test_2 = QuantumCircuit(4, 4)
qc_test_2.x(1)

qc_test_2_plus_3 = qc_test_2.compose(add_3)
qc_test_2_plus_3.draw(output='mpl', style="iqp")

sv_latex_from_qc(qc_test_2_plus_3, backend)

qc_test_7 = QuantumCircuit(4, 4)
qc_test_7.x(0)
qc_test_7.x(1)
qc_test_7.x(2)

qc_test_7_plus_8 = qc_test_7.compose(circuit_adder(8))
sv_latex_from_qc(qc_test_7_plus_8, backend)
#qc_test_7_plus_8.draw()

theta = 6.544985
phi = 2.338741
lmbda = 0

alice_1 = 0
alice_2 = 1
bob_1 = 2

qr_alice = QuantumRegister(2, 'Alice')
qr_bob = QuantumRegister(1, 'Bob')
cr = ClassicalRegister(3, 'c')
qc_ej3 = QuantumCircuit(qr_alice, qr_bob, cr)

qc_ej3.barrier(label='1')
qc_ej3.u(theta, phi, lmbda, alice_1);

qc_ej3.barrier(label='2')
qc_ej3.h(alice_2)
qc_ej3.cx(alice_2, bob_1);

qc_ej3.barrier(label='3')
qc_ej3.cx(alice_1, alice_2)
qc_ej3.h(alice_1);

qc_ej3.barrier(label='4')
qc_ej3.measure([alice_1, alice_2], [alice_1, alice_2]);

qc_ej3.barrier(label='5')
qc_ej3.x(bob_1).c_if(alice_2, 1)
qc_ej3.z(bob_1).c_if(alice_1, 1) 
qc_ej3.measure(bob_1, bob_1);

qc_ej3.draw(output='mpl', style="iqp")

result = backend.run(qc_ej3, shots=1024).result()
counts = result.get_counts()
plot_histogram(counts)

sv_0 = np.array([1, 0])
sv_1 = np.array([0, 1])

def find_symbolic_representation(value, symbolic_constants={1/np.sqrt(2): '1/√2'}, tolerance=1e-10):
    """
    Check if the given numerical value corresponds to a symbolic constant within a specified tolerance.

    Parameters:
    - value (float): The numerical value to check.
    - symbolic_constants (dict): A dictionary mapping numerical values to their symbolic representations.
                                Defaults to {1/np.sqrt(2): '1/√2'}.
    - tolerance (float): Tolerance for comparing values with symbolic constants. Defaults to 1e-10.

    Returns:
    str or float: If a match is found, returns the symbolic representation as a string 
                  (prefixed with '-' if the value is negative); otherwise, returns the original value.
    """
    for constant, symbol in symbolic_constants.items():
        if np.isclose(abs(value), constant, atol=tolerance):
            return symbol if value >= 0 else '-' + symbol
    return value


def array_to_dirac_notation(array, tolerance=1e-10):
    """
    Convert a complex-valued array representing a quantum state in superposition
    to Dirac notation.

    Parameters:
    - array (numpy.ndarray): The complex-valued array representing
      the quantum state in superposition.
    - tolerance (float): Tolerance for considering amplitudes as negligible.

    Returns:
    str: The Dirac notation representation of the quantum state.
    """
    # Ensure the statevector is normalized
    array = array / np.linalg.norm(array)

    # Get the number of qubits
    num_qubits = int(np.log2(len(array)))

    # Find indices where amplitude is not negligible
    non_zero_indices = np.where(np.abs(array) > tolerance)[0]

    # Generate Dirac notation terms
    terms = [
        (find_symbolic_representation(array[i]), format(i, f"0{num_qubits}b"))
        for i in non_zero_indices
    ]

    # Format Dirac notation
    dirac_notation = " + ".join([f"{amplitude}|{binary_rep}⟩" for amplitude, binary_rep in terms])

    return dirac_notation

def array_to_matrix_representation(array):
    """
    Convert a one-dimensional array to a column matrix representation.

    Parameters:
    - array (numpy.ndarray): The one-dimensional array to be converted.

    Returns:
    numpy.ndarray: The column matrix representation of the input array.
    """
    # Replace symbolic constants with their representations
    matrix_representation = np.array([find_symbolic_representation(value) or value for value in array])

    # Return the column matrix representation
    return matrix_representation.reshape((len(matrix_representation), 1))

def array_to_dirac_and_matrix_latex(array):
    """
    Generate LaTeX code for displaying both the matrix representation and Dirac notation
    of a quantum state.

    Parameters:
    - array (numpy.ndarray): The complex-valued array representing the quantum state.

    Returns:

    Latex: A Latex object containing LaTeX code for displaying both representations.
    """
    matrix_representation = array_to_matrix_representation(array)
    latex = "Matrix representation\n\\begin{bmatrix}\n" + \
            "\\\\\n".join(map(str, matrix_representation.flatten())) + \
            "\n\\end{bmatrix}\n"
    latex += f'Dirac Notation:\n{array_to_dirac_notation(array)}'
    return Latex(latex)

sv_b1 = np.kron(sv_0, sv_0)
array_to_dirac_and_matrix_latex(sv_b1)

sv_b1 = (np.kron(sv_0, sv_0) + np.kron(sv_0, sv_1)) / np.sqrt(2)
array_to_dirac_and_matrix_latex(sv_b1)

sv_b1 = (np.kron(sv_0, sv_0) + np.kron(sv_1, sv_1)) / np.sqrt(2)
array_to_dirac_and_matrix_latex(sv_b1)

sv_b2 = np.kron(sv_0, sv_0)
array_to_dirac_and_matrix_latex(sv_b2)

sv_b2 = (np.kron(sv_0, sv_0) + np.kron(sv_0, sv_1)) / np.sqrt(2)
array_to_dirac_and_matrix_latex(sv_b2)

sv_b2 = (np.kron(sv_0, sv_0) + np.kron(sv_1, sv_1)) / np.sqrt(2)
array_to_dirac_and_matrix_latex(sv_b2)

sv_b2 = (np.kron(sv_0, sv_0) - np.kron(sv_1, sv_1)) / np.sqrt(2)
array_to_dirac_and_matrix_latex(sv_b2)

sv_b3 = np.kron(sv_0, sv_0)
array_to_dirac_and_matrix_latex(sv_b3)

sv_b3 = np.kron(sv_0, sv_1)
array_to_dirac_and_matrix_latex(sv_b3)

sv_b3 = (np.kron(sv_0, sv_0) - np.kron(sv_0, sv_1)) / np.sqrt(2)
array_to_dirac_and_matrix_latex(sv_b3)

sv_b3 = (np.kron(sv_0, sv_1) - np.kron(sv_1, sv_0)) / np.sqrt(2)
array_to_dirac_and_matrix_latex(sv_b3)

sv_b3 = (np.kron(sv_0, sv_1) + np.kron(sv_1, sv_0)) / np.sqrt(2)
array_to_dirac_and_matrix_latex(sv_b3)

sv_b4 = np.kron(sv_0, sv_0)
array_to_dirac_and_matrix_latex(sv_b4)

sv_b4 = np.kron(sv_0, sv_1)
array_to_dirac_and_matrix_latex(sv_b4)

sv_b4 = (np.kron(sv_0, sv_0) - np.kron(sv_0, sv_1)) / np.sqrt(2)
array_to_dirac_and_matrix_latex(sv_b4)

sv_b4 = (np.kron(sv_0, sv_1) - np.kron(sv_1, sv_0)) / np.sqrt(2)
array_to_dirac_and_matrix_latex(sv_b4)