import numpy as np
import matplotlib.pylab as plt
import seaborn as sns

np.sign(2*np.random.rand(10)-1).astype('int')

class Ising2D():
    def __init__(self, N=10, seed=8848):
        
        self.N = N
        self.seed = seed
        np.random.seed(self.seed)
        # construct a lattice with random spins (+1, -1)
        self.Lattice = np.array([ np.sign(2*np.random.rand(self.N)-1).astype('int') for i in range(self.N)] )

    def Print_lattice(self):
        for row in self.Lattice:
            print (row)

    def Plot_lattice(self):
        fig,ax=plt.subplots(1,1, figsize=(6,6))
        sns.heatmap(self.Lattice, linecolor='white', cbar=True, ax=ax,\
                    square=True, cmap='viridis', xticklabels=[], yticklabels=[]);
        plt.show()
        
    def Energy(self, J=1):
        # one liner using list comprehension
        ene = np.sum([ -1.0*J*self.Lattice[i,j] *( self.Lattice[ (i+1)%self.N, j ] + self.Lattice[ (i-1)%self.N, j ] +\
                                            self.Lattice[ i, (j+1)%self.N ] + self.Lattice[ i, (j-1)%self.N ] ) \
                for i in range(self.N) for j in range(self.N) ])
        return ene
        
    def Magnetization(self):
        return np.sum(self.Lattice)

model = Ising2D(N=10, seed=1) # seed 0, N=12
model.Plot_lattice()
#model.Print_lattice()
Ene = model.Energy()
M = model.Magnetization()
print ( f"Energy of given configuration: {Ene} and Magnetization: {M}")

def calculate_energy(J, σ, h=[0]):
    energy = -sum(J_ij*σ[i]*σ[i+1] for i, J_ij in enumerate(J))
    energy += -sum(h_i*σ_i for h_i, σ_i in zip(h, σ))
    return energy

J = [1.0, -1.0]
σ = [+1, -1, +1]

calculate_energy(J, σ)

import itertools
for σ in itertools.product(*[{+1, -1}, {+1, -1}, {+1, -1}]):
#for σ in itertools.product(*[{+1,-1} for _ in range(3)]):
    print(f"Spin configuration:\t{σ} \tEnergy: \t {calculate_energy(J, σ)}")

J = [1., 1.]
σ = [+1, -1, +1]
h = [0.5, 0.5, 0.4]
print ( f"J={J} \tEnergy= {calculate_energy(J, σ, h)}")

J = [-1., -1.]
σ = [+1, -1, +1]
h = [0.5, 0.5, 0.4]
print ( f"J={J} \tEnergy= {calculate_energy(J, σ, h)}")

import itertools
J = [-1, -1]
h = [0, 0, 0]
minimum_energy = 0
for σ_c in itertools.product(*[{+1,-1} for _ in range(3)]):
    energy = calculate_energy(J, σ_c, h)
    if energy < minimum_energy:
        minimum_energy = energy
        σ = σ_c

σ, minimum_energy

#!pip3 install dimod

import dimod

J = {(0, 1): 1.0, (1, 2): -1.0}
h = {0:0, 1:0, 2:0}

model = dimod.BinaryQuadraticModel(h, J, 0.0, dimod.SPIN) # syntax (linear, quadratic, offset, vartype)

sampler = dimod.SimulatedAnnealingSampler()
response = sampler.sample(model, num_reads=10)

all_ene = [solution.energy for solution in response.data()]
print (f"All energies: {all_ene} with total count: {all_ene.count(-2.)}")

from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import execute
from qiskit import Aer
from qiskit.visualization import plot_state_qsphere
from qiskit_textbook.tools import array_to_latex


np.set_printoptions(precision=3, suppress=True)
backend = Aer.get_backend('statevector_simulator')
q = QuantumRegister(1)
c = ClassicalRegister(1)
circuit = QuantumCircuit(q, c)
display ( circuit.draw('mpl') )

def get_statevector(qc):
    return execute(qc, backend = Aer.get_backend('statevector_simulator') ).result().get_statevector()

circuit.z(q[0])
state = get_statevector(circuit)
array_to_latex ( state, pretext="\\text{ State : }\n" )
plot_state_qsphere(state.data)

circuit = QuantumCircuit(q, c)
circuit.z(q[0])
circuit.x(q[0])

state = get_statevector(circuit)
array_to_latex ( state, pretext="\\text{ State : }\n" )
plot_state_qsphere(state.data)

def calculate_energy_expectation(state, hamiltonian):
    return np.squeeze(np.dot(np.matrix.getH(state),  np.dot(hamiltonian, state) ))
    #return float((state.T.conj() @ hamiltonian @ state).real)
    #return float((state.T.conj() @ hamiltonian @ state).real)

# built-in operators
from qiskit.quantum_info import Operator, Pauli

array_to_latex ( Operator(Pauli(label='X')).data, pretext="\\text{ Pauli-X operator : }\n" )
array_to_latex ( Operator(Pauli(label='XZ')).data, pretext="\\text{ XZ operator : }\n" )

PauliZ = np.array([[1, 0], [0, -1]])
I2 = np.eye(2)
IZ = np.kron(I2, PauliZ);
ZI = np.kron(PauliZ, I2);
ZZ = np.kron(PauliZ, PauliZ);
H = -ZZ + -0.5*(ZI+IZ);
ψ = np.kron([[1], [0]], [[1], [0]]) ;
ψ_H_ψ = calculate_energy_expectation(ψ, H);

array_to_latex ( IZ, pretext="\\text{ IZ = }\n" )
array_to_latex ( ZI, pretext="\\text{ ZI = }\n" )
array_to_latex ( ZZ, pretext="\\text{ ZZ = }\n" )
array_to_latex ( H, pretext="\\text{ H = }\n" )
array_to_latex ( ψ, pretext="\\text{ ψ = }\n" )
 
print (f"Energy expectetation value: {ψ_H_ψ}")


circuit = QuantumCircuit(q, c)
circuit.x(q[0]); 
circuit.z(q[0])
state = get_statevector(circuit)

op = Operator(circuit)
array_to_latex ( op.data, pretext="\\text{ Pauli-X, then Pauli-Z, Circuit as operator : }\n" )
array_to_latex ( state, pretext="\\text{ Statevector : }\n" )

circuit = QuantumCircuit(q, c)
circuit.z(q[0]);
circuit.x(q[0])
state = get_statevector(circuit)

op = Operator(circuit)
array_to_latex ( op.data, pretext="\\text{ Pauli-Z, then Pauli-X, Circuit as operator : }\n" )
array_to_latex ( state, pretext="\\text{ Statevector: }\n" )


X = np.array([[0, 1], [1, 0]])
IX = np.kron(I2, X)
XI = np.kron(X, I2)
H = H - XI - IX
H

ψ = np.kron([[1], [0]], [[1], [0]])
calculate_energy_expectation(ψ, H)

eigenvalues, eigenvectors = np.linalg.eigh(H)
ψ = eigenvectors[:, 0:1]
energy_expectation_value = calculate_energy_expectation(ψ, H)
energy_expectation_value

from qiskit.providers.aer.noise import NoiseModel
from qiskit.providers.aer.noise.errors import pauli_error, depolarizing_error

from qiskit import ignis
from qiskit.visualization import plot_histogram


# Construct a 5% single-qubit Bit-flip error
p_error = 0.05
bit_flip = pauli_error([('X', p_error), ('I', 1 - p_error)])
phase_flip = pauli_error([('Z', p_error), ('I', 1 - p_error)])
print(bit_flip)
print(phase_flip)

def get_noise(p):
    # flip the bit with a given p
    error_meas = pauli_error([('X',p), ('I', 1 - p)])
    noise_model = NoiseModel()
    noise_model.add_all_qubit_quantum_error(error_meas, "measure") 
    # measurement error is applied to measurements
    return noise_model

noise_model = get_noise(0.01)
noise_model

qc = QuantumCircuit(2,2)
qc.x(0)
qc.h(1)
qc.measure(0, 1)
display( qc.draw('mpl') )

backend = Aer.get_backend('qasm_simulator')

noises = np.array([0.01, 0.1, .5])
noises_label = ['noise='+str(n) for n in noises]

all_counts =[]
for noise in noises:
    
    noise_model = get_noise(noise)
    job = execute(qc, backend=backend, noise_model=noise_model, shots=10000)
    counts = job.result().get_counts()
    
    all_counts.append(counts)

plot_histogram(all_counts, legend=noises_label, figsize=(6, 4))



noise_model = get_noise(0.01)

qc = QuantumCircuit(2,2)
for state in ['00','01','10','11']:
    if state[0]=='1':
        qc.x(1)
    if state[1]=='1':
        qc.x(0)
    qc.measure(qc.qregs[0],qc.cregs[0])
    
    backend = Aer.get_backend('qasm_simulator')
    job = execute(qc, backend=backend, noise_model=noise_model, shots=10000)
    counts = job.result().get_counts()
    print ( f"Pure State: {state}, mixed states counts: {counts}" )

def Get_matrix():
    M = np.zeros((2**2, 2**2), dtype=float)
    for i, state in enumerate(['00','01','10','11']):
        qc = QuantumCircuit(2,2)
        if state[0]=='1': qc.x(1)
        if state[1]=='1': qc.x(0)

        qc.measure(qc.qregs[0],qc.cregs[0])
        
        backend = Aer.get_backend('qasm_simulator')
        job = execute(qc, backend=backend, noise_model=noise_model, shots=10000)
        counts = job.result().get_counts()
        
        M[i] = np.array([counts.get(s, 0) for s in ['00','01','10','11'] ])/10000
        
    return M

M = Get_matrix()
array_to_latex(M, pretext="M=")

C_ideal = np.array([0, 5000, 5000, 0] )
C_noisy = np.dot(M, C_ideal)
print (f"C_ideal: {C_ideal}")
print (f"C_Noisy: {C_noisy}")

C_mitigated = np.dot(np.linalg.inv(M), C_noisy)
print (f"C_mitigated: {C_mitigated}")

from qiskit.ignis.mitigation.measurement import (complete_meas_cal, CompleteMeasFitter)

qr = QuantumRegister(2)
meas_calibs, state_labels = complete_meas_cal(qr=qr, circlabel='mcal')

print ("The basis states: ", state_labels)

for circuit in meas_calibs:
    print ("Circuit:",circuit.name)
    display(circuit.draw('mpl'))


# Execute the calibration circuits without noise
backend = Aer.get_backend('qasm_simulator')
job = execute(meas_calibs, backend=backend, shots=1000)
cal_results = job.result()
#cal_results

meas_fitter = CompleteMeasFitter(cal_results, state_labels, circlabel='mcal')
print(meas_fitter.cal_matrix)

noise_model = get_noise(0.1)

backend = Aer.get_backend('qasm_simulator')
job = execute(meas_calibs, backend=backend, shots=1000, noise_model=noise_model)
cal_results = job.result()

meas_fitter = CompleteMeasFitter(cal_results, state_labels, circlabel='mcal')
print(meas_fitter.cal_matrix)

qc = QuantumCircuit(2,2)
qc.h(0)
qc.cx(0,1)
display ( qc.draw('mpl') )
qc.measure(qc.qregs[0],qc.cregs[0])

job = execute(qc, backend=backend, shots=10000, noise_model=noise_model)
results = job.result()

noisy_counts = results.get_counts()
print("Noisy Counts:",noisy_counts)

# Get the filter object
meas_filter = meas_fitter.filter

# Results with mitigation
mitigated_results = meas_filter.apply(results)
mitigated_counts = mitigated_results.get_counts(0)
print ("mitigated counts:\n", mitigated_counts)

print ( noisy_counts)
print ( mitigated_counts )

%config InlineBackend.figure_format = 'svg' # Makes the images look nice
plot_histogram([noisy_counts, mitigated_counts], legend=['noisy', 'mitigated'], figsize=(6, 4))





