import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

from qiskit import(QuantumCircuit, QuantumRegister, ClassicalRegister, transpile)
from qiskit.visualization import plot_histogram, plot_gate_map
from qiskit.quantum_info import Statevector
from qiskit.circuit.library import QFT
from qiskit_aer import AerSimulator

from qiskit_ibm_provider import IBMProvider
plt.style.use('default')


#IBMProvider.save_account('YOUR_TOKEN',overwrite=True) 
provider = IBMProvider(instance="ibm-q/open/main")

print(provider.backends(min_num_qubits=5, simulator=False, operational=True))

print(provider.backends(min_num_qubits=5, simulator=True, operational=True))

backend = provider.get_backend("ibm_lagos")

backend.status()

props = backend.properties()


def describe_qubit(qubit, properties):
    """Print a string describing some of reported properties of the given qubit."""

    # Conversion factors from standard SI units
    us = 1e6
    ns = 1e9
    GHz = 1e-9

    print("Qubit {0} has a \n"
          "  - T1 time of {1} microseconds\n"
          "  - T2 time of {2} microseconds\n"
          "  - U2 gate error of {3}\n"
          "  - U2 gate duration of {4} nanoseconds\n"
          "  - resonant frequency of {5} GHz".format(
              qubit,
              properties.t1(qubit) * us,
              properties.t2(qubit) * us,
              properties.gate_error('x', qubit),
              properties.gate_length('x', qubit) * ns,
              properties.frequency(qubit) * GHz))

describe_qubit(0, props)

plot_gate_map(backend)

simulator = AerSimulator()
backendSimulator =  AerSimulator.from_backend(backend)

qr = QuantumRegister(2, 'q')
cr = ClassicalRegister(2, 'c')
circuit = QuantumCircuit(qr, cr,name="bell state")
circuit.h(0)
circuit.cx(0,1)
circuit.draw('mpl')

vector=Statevector(circuit)
display(vector.draw('latex'))

circuit = QuantumCircuit(qr, cr,name="bell state")
circuit.h(0)
circuit.cx(0,1)
circuit.measure([0,1],[0,1])
circuit.draw('mpl')

counts=simulator.run(circuit,shots=1000).result().get_counts()

display(counts)
plot_histogram(counts)

circuit = QuantumCircuit(qr, cr,name="bell state")
circuit.h(0)
circuit.cx(0,1)
circuit.measure([0,1],[0,1])
circuit.draw('mpl')

layout=[0,1]
circuitTranspiled=transpile(circuit,backendSimulator,initial_layout=layout,optimization_level=1)
circuitTranspiled.draw('mpl')

counts=backendSimulator.run(circuitTranspiled,shots=1000).result().get_counts()

display(counts)
plot_histogram(counts)

layout=[0,6]
circuitTranspiled=transpile(circuit,backendSimulator,initial_layout=layout,optimization_level=1)
circuitTranspiled.draw('mpl')

counts=backendSimulator.run(circuitTranspiled,shots=1000).result().get_counts()

display(counts)
plot_histogram(counts)

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService(channel="ibm_quantum")
backend = service.backend("ibmq_qasm_simulator")

#backend = provider.get_backend("ibmq_qasm_simulator")

from qiskit.circuit.random import random_circuit

circuit = random_circuit(2, 2, seed=0, measure=True).decompose(reps=1)
circuit.draw("mpl")

from qiskit_ibm_runtime import Sampler

sampler = Sampler(backend=backend)

#from qiskit.primitives import Sampler

#sampler = Sampler()

job = sampler.run(circuit)

print(f">>> Job ID: {job.job_id()}")
print(f">>> Job Status: {job.status()}")

result = job.result()
print(f">>> {result}")
print(f"  > Quasi-probability distribution: {result.quasi_dists[0]}")
print(f"  > Metadata: {result.metadata[0]}")

from qiskit_ibm_runtime import Options

options = Options()
options.optimization_level = 1
options.execution.shots = 2048

sampler = Sampler(backend=backend, options=options)
job=sampler.run(circuit)

print(f">>> Job ID: {job.job_id()}")
print(f">>> Job Status: {job.status()}")

result = job.result()

print(f">>> Metadata: {result.metadata[0]}")

from qiskit_aer.noise import NoiseModel
device = provider.get_backend("ibm_lagos")
noise_model = NoiseModel.from_backend(device)
coupling_map = device.configuration().coupling_map
basis_gates=device.configuration().basis_gates


options_with_nm = Options()
options_with_nm.simulator = {
    "noise_model": noise_model,
    "basis_gates": basis_gates,
    "coupling_map": coupling_map,
    "seed_simulator": 42
}

# Set number of shots, optimization_level and resilience_level
options_with_nm.execution.shots = 1000
options_with_nm.optimization_level = 0 # no optimization
options_with_nm.resilience_level = 1 # M3 for Sampler and T-REx for Estimator

sampler = Sampler(backend=backend, options=options_with_nm)
job = sampler.run(circuit)

print(f">>> Job ID: {job.job_id()}")
print(f">>> Job Status: {job.status()}")

result=job.result()
print(f"  > Quasi-probability distribution: {result.quasi_dists[0]}")
plot_histogram(result.quasi_dists[0])

circuitList=[]

for t in np.linspace(0,10000,101):
    circuit = QuantumCircuit(1, 1)
    circuit.x(0)
    #circuit.h(0) 
    circuit.barrier()
    for i in range(int(t)):
        circuit.id(0)
        circuit.barrier()
    #circuit.h(0)
    circuit.measure(0,0)
    circuitList.append(circuit)
   
circuitListTranspiled=transpile(circuitList,backendSimulator,initial_layout=[0],optimization_level=0)

countsList=backendSimulator.run(circuitListTranspiled,shots=1000).result().get_counts()

probabiltyVec=[]
for counts in countsList:
    probabiltyVec.append(counts.get('1')/1000)

timeVec=np.linspace(0,10000,101)* props.gate_length('id', 0)   *10**6
plt.plot(timeVec,probabiltyVec)
plt.xlabel('Time [$\mu$s]')
plt.ylabel('$p_1$')
plt.show()

def func(x, amp, tau,c):
    return amp* np.exp(-x/tau)+c

popt, pcov = curve_fit(func, timeVec, probabiltyVec) # scipy's curve fitting tool
print(popt)

plt.plot(timeVec,probabiltyVec)
plt.plot(timeVec, func(timeVec, *popt), label="Fitted Curve") 
plt.xlabel('Time [$\mu$s]')
plt.ylabel('$p_1$')
plt.legend()
plt.show()

transpile(QFT(4,do_swaps=True,inverse=False,insert_barriers=True),simulator).draw('mpl')

backend = provider.get_backend("ibm_brisbane")
plot_gate_map(backend,figsize=(16,10))


backendSimulator =  AerSimulator.from_backend(backend)

circuitQFT=QFT(5,do_swaps=True,inverse=False)
circuitQFT.draw('mpl')

nQubit=8

qr = QuantumRegister(nQubit, 'q')
circuit = QuantumCircuit(qr,name=" QFT calculation")

#circuit.x(4)
for i in range(nQubit):
    circuit.h(i)
circuit.barrier()
circuit.append(QFT(nQubit,do_swaps=True,inverse=True),range(nQubit))

circuit.draw('mpl')

vector=Statevector(circuit)
display( vector.draw('latex') )

circuit.measure_all()
circuit.draw('mpl')

layout=[0,34,4,8,100,45,12,5]
circuitTranspiled=transpile(circuit,backendSimulator,initial_layout=layout,optimization_level=1)
#circuitTranspiled.draw('mpl')

# NON ESEGUIRE!!!!!
counts=backendSimulator.run(circuitTranspiled,shots=1000).result().get_counts()

layout=[0,1,2,3,4,5,6,7]
circuitTranspiled=transpile(circuit,backendSimulator,initial_layout=layout,optimization_level=1)
#circuitTranspiled.draw('mpl')

counts=backendSimulator.run(circuitTranspiled,shots=1000).result().get_counts()

#display(counts)
plot_histogram(counts,figsize=(26,6))

nQubit=8

qr = QuantumRegister(nQubit, 'q')
circuit = QuantumCircuit(qr,name=" QFT calculation")

circuit.x(3)
circuit.append(QFT(nQubit,do_swaps=False,inverse=False),range(nQubit))
circuit.barrier()
circuit.append(QFT(nQubit,do_swaps=False,inverse=True),range(nQubit))
circuit.barrier()
circuit.append(QFT(nQubit,do_swaps=False,inverse=False),range(nQubit))
circuit.barrier()
circuit.append(QFT(nQubit,do_swaps=False,inverse=True),range(nQubit))

circuit.measure_all()

layout=[0,1,2,3,4,5,6,7]
circuitTranspiled=transpile(circuit,backendSimulator,initial_layout=layout,optimization_level=1)
counts=backendSimulator.run(circuitTranspiled,shots=1000).result().get_counts()
plot_histogram(counts,figsize=(40,6))

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService(channel="ibm_quantum")
backend = service.backend("ibmq_qasm_simulator")

from qiskit.circuit.random import random_circuit
circuit = random_circuit(2, 2, seed=0).decompose(reps=1)
display(circuit.draw("mpl"))

from qiskit.quantum_info import SparsePauliOp
observable = SparsePauliOp("XZ")
print(f">>> Observable: {observable.paulis}")

from qiskit_ibm_runtime import Estimator

estimator = Estimator(backend=backend)

job = estimator.run(circuit, observable)

print(f">>> Job ID: {job.job_id()}")
print(f">>> Job Status: {job.status()}")

result = job.result()
print(f">>> {result}")
print(f"  > Expectation value: {result.values[0]}")

circuits = (
    random_circuit(2, 2, seed=0).decompose(reps=1),
    random_circuit(2, 2, seed=1).decompose(reps=1),
)
observables = (
    SparsePauliOp("XZ"),
    SparsePauliOp("IY"),
)

job = estimator.run(circuits, observables)


print(f">>> Job ID: {job.job_id()}")
print(f">>> Job Status: {job.status()}")

result = job.result()
[display(cir.draw("mpl")) for cir in circuits]
print(f">>> Observables: {[obs.paulis for obs in observables]}")
print(f">>> Expectation values: {result.values.tolist()}")

from qiskit.circuit.library import RealAmplitudes

circuit = RealAmplitudes(num_qubits=2, reps=2).decompose(reps=1)
observable = SparsePauliOp("ZI")
parameter_values = [0, 1, 2, 3, 4, 5]

job = estimator.run(circuit, observable, parameter_values)

print(f">>> Job ID: {job.job_id()}")
print(f">>> Job Status: {job.status()}")


result = job.result()
display(circuit.draw("mpl"))
print(f">>> Observable: {observable.paulis}")
print(f">>> Parameter values: {parameter_values}")
print(f">>> Expectation value: {result.values[0]}")

# General imports
import time
import numpy as np

# Pre-defined ansatz circuit and operator class for Hamiltonian
from qiskit.circuit.library import EfficientSU2
from qiskit.quantum_info import SparsePauliOp

# The IBM Qiskit Runtime / Primitives
#from qiskit_ibm_runtime import QiskitRuntimeService
#from qiskit_ibm_runtime import Estimator, Session
from qiskit.primitives import Estimator

# SciPy minimizer routine
from scipy.optimize import minimize

# Plotting functions
import matplotlib.pyplot as plt



hamiltonian = SparsePauliOp.from_list(
    [("YZ", 0.3980), ("ZI", -0.3980), ("ZZ", -0.0113), ("XX", 0.1810)]
)

ansatz = EfficientSU2(hamiltonian.num_qubits)
ansatz.draw("mpl")

transpile(ansatz,simulator).draw('mpl')

num_params = ansatz.num_parameters
num_params

def cost_func(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    energy = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
    return energy

x0 = 2 * np.pi * np.random.random(num_params)

estimator = Estimator(options={"shots": int(1e4)})
res = minimize(cost_func, x0, args=(ansatz, hamiltonian, estimator), method="cobyla")

display(res)

def build_callback(ansatz, hamiltonian, estimator, callback_dict):
    """Return callback function that uses Estimator instance,
    and stores intermediate values into a dictionary.

    Parameters:
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance
        callback_dict (dict): Mutable dict for storing values

    Returns:
        Callable: Callback function object
    """

    def callback(current_vector):
        """Callback function storing previous solution vector,
        computing the intermediate cost value, and displaying number
        of completed iterations and average time per iteration.

        Values are stored in pre-defined 'callback_dict' dictionary.

        Parameters:
            current_vector (ndarray): Current vector of parameters
                                      returned by optimizer
        """
        # Keep track of the number of iterations
        callback_dict["iters"] += 1
        # Set the prev_vector to the latest one
        callback_dict["prev_vector"] = current_vector
        # Compute the value of the cost function at the current vector
        callback_dict["cost_history"].append(
            estimator.run(ansatz, hamiltonian, parameter_values=current_vector).result().values[0]
        )
        # Grab the current time
        current_time = time.perf_counter()
        # Find the total time of the execute (after the 1st iteration)
        if callback_dict["iters"] > 1:
            callback_dict["_total_time"] += current_time - callback_dict["_prev_time"]
        # Set the previous time to the current time
        callback_dict["_prev_time"] = current_time
        # Compute the average time per iteration and round it
        time_str = (
            round(callback_dict["_total_time"] / (callback_dict["iters"] - 1), 2)
            if callback_dict["_total_time"]
            else "-"
        )
        # Print to screen on single line
        print(
            "Iters. done: {} [Avg. time per iter: {}]".format(callback_dict["iters"], time_str),
            end="\r",
            flush=True,
        )

    return callback

callback_dict = {
    "prev_vector": None,
    "iters": 0,
    "cost_history": [],
    "_total_time": 0,
    "_prev_time": None,}

estimator = Estimator(options={"shots": int(1e4)})
callback = build_callback(ansatz, hamiltonian, estimator, callback_dict)
res = minimize(cost_func,x0,args=(ansatz, hamiltonian, estimator),method="cobyla",callback=callback)

all(callback_dict["prev_vector"] == res.x)

callback_dict["iters"] == res.nfev

fig, ax = plt.subplots()
ax.plot(range(callback_dict["iters"]), callback_dict["cost_history"])
ax.set_xlabel("Iterations")
ax.set_ylabel("Cost");

# General imports
import time
import numpy as np

# The IBM Qiskit Runtime / Primitives
from qiskit.primitives import Sampler

# SciPy minimizer routine
from scipy.optimize import minimize

# Datasets
from sklearn.datasets import make_moons

# Plotting functions
import matplotlib.pyplot as plt


# Set random seeds
np.random.seed(42)

# number of samples
N=500

# Let's use the Moons dataset
X, Y = make_moons(n_samples=N, noise=0.1)

c = ["#1f77b4" if y_ == 0 else "#ff7f0e" for y_ in Y]  # colours for each class
#plt.axis("off")
plt.scatter(X[:, 0], X[:, 1], c=c)
plt.show()

display(X,Y)


from qiskit.circuit import ParameterVector

# data re-uploading classifier circuit

L=5

x=ParameterVector('x',2)
thetas = ParameterVector('θ',3*L)

qr = QuantumRegister(1, name="q")
cr = ClassicalRegister(1, name='c')
ansatz = QuantumCircuit(qr, cr)
for k in range(L):
    ansatz.u(x[0],x[1],0.0,qr)
    ansatz.u(thetas[3*k+0], thetas[3*k+1], thetas[3*k+2], qr) # su(2) rotation
ansatz.measure(qr, cr)
ansatz.draw('mpl')

num_params = ansatz.num_parameters
display(num_params)

params0 = 2 * np.pi * np.random.random(3*L)

circuit=ansatz.assign_parameters({x:X[0], thetas:params0})
circuit.draw('mpl')

sampler=Sampler()

def loss(params, ansatz, x , y , sampler):
    """Return the loss for a pattern 

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        x (numpyArray): Array of the pattern features
        y (integer): Class of the pattern 0 or 1
        sampler (Estimator): Sampler primitive instance

    Returns:
        float: loss function (cost function)
    """
    parameter_values=np.concatenate((x,params))
    job = sampler.run(ansatz, parameter_values=parameter_values)
    loss = 1 - job.result().quasi_dists[0].get(y)
    return loss

def cost_func( params, ansatz, trainingX , trainingY , sampler):
    N=len(trainingY)
    return sum(loss(params,ansatz,trainingX[k],trainingY[k],sampler) for k in range(N))/N

def predict(params,x,ansatz,sampler):
    if loss(params,ansatz,x,0,sampler) < 0.5:
        return 0
    else:
        return 1

cost_func(params0,ansatz,X,Y,sampler)

Ypred = [predict(params0,x,ansatz,sampler) for x in X]
c = ["#1f77b4" if y == 0 else "#ff7f0e" for y in Y]  # colours for each class
cpred = ["#1f77b4" if y == 0 else "#ff7f0e" for y in Ypred]  # colours for each class
plt.figure(figsize=(15,5))
plt.subplot(1,2,1)
plt.axis("off")
plt.title("true")
plt.scatter(X[:, 0], X[:, 1], c=c)
plt.subplot(1,2,2)
plt.axis("off")
plt.title("predicted")
plt.scatter(X[:, 0], X[:, 1], c=cpred)
plt.show()

trainingX=X[0:100]
trainingY=Y[0:100]
res = minimize(cost_func, params0, args=(ansatz, trainingX , trainingY , sampler), method="cobyla")

res.fun


paramsOpt=res.x
Ypred = [predict(paramsOpt,x,ansatz,sampler) for x in X[200:]]
c = ["#1f77b4" if y == 0 else "#ff7f0e" for y in Y[200:]]  # colours for each class
cpred = ["#1f77b4" if y == 0 else "#ff7f0e" for y in Ypred]  # colours for each class
plt.figure(figsize=(15,5))
plt.subplot(1,2,1)
plt.axis("off")
plt.title("true")
plt.scatter(X[200:, 0], X[200:, 1], c=c)
plt.subplot(1,2,2)
plt.axis("off")
plt.title("predicted")
plt.scatter(X[200:, 0], X[200:, 1], c=cpred)
plt.show()

from qutip import basis,Bloch, expect, sigmax,sigmay,sigmaz
#import qutip as qt
L=5
x=ParameterVector('x',2)
thetas = ParameterVector('θ',3*L)
qr = QuantumRegister(1, name="q")
cr = ClassicalRegister(1, name='c')
qc = QuantumCircuit(qr, cr)
for k in range(L):
    qc.u(x[0],x[1],0.0,qr)
    qc.u(thetas[3*k+0], thetas[3*k+1], thetas[3*k+2], qr) # su(2) rotation


def blochPoints(X,Y,params,ansatz):
    xp1=[]
    yp1=[]
    zp1=[]

    xp0=[]
    yp0=[]
    zp0=[]

    for i in range(len(X)):
        circuit=ansatz.assign_parameters({x:X[i], thetas:params})
        state=Statevector(circuit)
  
        pnt=state[0]*basis(2,0)+state[1]*basis(2,1)
        if Y[i]==0:
            xp0.append(expect(sigmax(),pnt))
            yp0.append(expect(sigmay(),pnt))
            zp0.append(expect(sigmaz(),pnt))
        else:
            xp1.append(expect(sigmax(),pnt))
            yp1.append(expect(sigmay(),pnt))
            zp1.append(expect(sigmaz(),pnt))
    return [[xp0,yp0,zp0],[xp1,yp1,zp1]]


b1=Bloch()
b2=Bloch()

b1.add_points(blochPoints(X[200:300],Y[200:300],params0,qc)[0])
b1.add_points(blochPoints(X[200:300],Y[200:300],params0,qc)[1])

b2.add_points(blochPoints(X[200:300],Y[200:300],paramsOpt,qc)[0])
b2.add_points(blochPoints(X[200:300],Y[200:300],paramsOpt,qc)[1])

b1.show()
b2.show()


import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright