from qiskit import QuantumCircuit
from qiskit.primitives import Sampler, BackendSampler
from qiskit.circuit.library import *
from qiskit.circuit.random import random_circuit
from qiskit.visualization import plot_histogram, array_to_latex, plot_bloch_multivector
from qiskit.quantum_info import Statevector

from qiskit_ibm_runtime import Sampler, QiskitRuntimeService, Options, Session
from qiskit.providers.fake_provider import FakeManila
from qiskit_aer.noise import NoiseModel

import numpy as np

print('X (1q gate)')
xgate = XGate()
array_to_latex(xgate.to_matrix())

qc_x = QuantumCircuit(1)
qc_x.x(0)
qc_x.draw('mpl', scale=1.5)

state = Statevector(qc_x)
plot_bloch_multivector(state)

qc_x = QuantumCircuit(1,1)
qc_x.x(0)
qc_x.measure([0],[0])
qc_x.draw('mpl', scale=1.5)

# load the service and set the backend to the simulator
service = QiskitRuntimeService(channel="ibm_quantum")
backend = "ibmq_qasm_simulator"

# Make a noise model
fake_backend = FakeManila()
noise_model = NoiseModel.from_backend(fake_backend)

# Set options to include the noise model
options = Options()
options.simulator = {
    "noise_model": noise_model,
    "basis_gates": fake_backend.configuration().basis_gates,
    "coupling_map": fake_backend.configuration().coupling_map,
    "seed_simulator": 42
}

# Set number of shots, optimization_level and resilience_level
options.execution.shots = 1000
options.optimization_level = 0
options.resilience_level = 0

with Session(service=service, backend=backend):
    sampler = Sampler(options=options)
    job = sampler.run(qc_x)


result_prob = job.result().quasi_dists[0].binary_probabilities()
plot_histogram(result_prob)

print('H (1q gate)')
hgate = HGate()
array_to_latex(hgate.to_matrix())

qc_h = QuantumCircuit(1,1)
qc_h.h(0)
qc_h.draw('mpl', scale=1.5)

state = Statevector(qc_h)
plot_bloch_multivector(state)

qc_h.measure([0],[0])

with Session(service=service, backend=backend):
    sampler = Sampler(options=options)
    job = sampler.run(qc_h)

result_prob = job.result().quasi_dists[0].binary_probabilities()
plot_histogram(result_prob)

options.resilience_level = 1  # add in M3 measurement mitigation

with Session(service=service, backend=backend):
    sampler = Sampler(options=options)
    job = sampler.run(qc_h)

result_prob = job.result().quasi_dists[0].binary_probabilities()
plot_histogram(result_prob)

qc_cx = QuantumCircuit(2)
qc_cx.cx(0,1)
qc_cx.draw('mpl')

print('CX (2q gate)')
cxgate = CXGate()
array_to_latex(cxgate.to_matrix())

qc2 = QuantumCircuit(2, 2)
qc2.h(0)
qc2.h(1)
qc2.draw('mpl', scale=1.5)

state = Statevector(qc2)
plot_bloch_multivector(state)

qc2.measure([0,1],[0,1])

options.resilience_level = 0  # no error mitigation

with Session(service=service, backend=backend):
    sampler = Sampler(options=options)
    job = sampler.run(qc2)

result_prob = job.result().quasi_dists[0].binary_probabilities()
plot_histogram(result_prob)

options.resilience_level = 1  # add in M3 measurement mitigation

with Session(service=service, backend=backend):
    sampler = Sampler(options=options)
    job = sampler.run(qc2)

result_prob = job.result().quasi_dists[0].binary_probabilities()
plot_histogram(result_prob)

qc_bell = QuantumCircuit(2, 2)
qc_bell.h(0)
qc_bell.cx(0,1)
qc_bell.draw('mpl')

state = Statevector(qc_bell)
plot_bloch_multivector(state)

qc_bell.measure([0,1],[0,1]) # add the measurement
options.resilience_level = 0  # no mitigation

with Session(service=service, backend=backend):
    sampler = Sampler(options=options)
    job = sampler.run(qc_bell)

result_prob = job.result().quasi_dists[0].binary_probabilities()
plot_histogram(result_prob)

options.resilience_level = 1  # add in M3 measurement mitigation

with Session(service=service, backend=backend):
    sampler = Sampler(options=options)
    job = sampler.run(qc_bell)

result_prob = job.result().quasi_dists[0].binary_probabilities()
plot_histogram(result_prob)

print(result_prob)

least_busy_device = service.least_busy(filters=lambda b: b.configuration().simulator==False)
least_busy_device

hw_result = None
with Session(service, backend=least_busy_device) as session:
    
    options = Options(resilience_level=0)
    sampler0 = Sampler(options=options)
    options = Options(resilience_level=1)
    sampler1 = Sampler(options=options)
    job0 = sampler0.run(circuits=[qc_bell], shots=8000)
    job1 = sampler1.run(circuits=[qc_bell], shots=8000)
    
    # You can see that the results are quasi probabilities, not counts
    hw_results = [job0.result().quasi_dists[0].binary_probabilities(), 
                  job1.result().quasi_dists[0].binary_probabilities()]

print('Results resilience level 0 ', hw_results[0])
print('Results resilience level 1 ', hw_results[1])

plot_histogram(hw_results, legend=['resilience 0', 'resilience 1'])