from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, execute, Aer
from qiskit.tools.visualization import plot_histogram, plot_bloch_multivector
#import numpy as np

q = QuantumRegister(2)
c = ClassicalRegister(2)

# Define Quantum circuit with 2 quantum register and 2 classical register
bell = QuantumCircuit(q,c)

# Create a Bell state
bell.x(q[0])
bell.x(q[1])
bell.h(q[0])
bell.cx(q[0], q[1])

# Draw circuit
bell.draw(output = 'mpl')

backend_statevector = Aer.get_backend('statevector_simulator')
job = execute(bell, backend_statevector)
result = job.result()
outputstate = result.get_statevector(bell, decimals = 3)
print(outputstate)
#plot_bloch_multivector(job.result().get_statevector(bell))

bell.measure(q, c)
bell.draw('mpl')

backend = Aer.get_backend('qasm_simulator')
job = execute(bell, backend)
result = job.result()
result.get_counts(bell)
plot_histogram(job.result().get_counts(bell))

import numpy as np
zero_ket = np.array([[1], [0]])
one_ket = np.array([[0], [1]])
ψ = (zero_ket + one_ket)/np.sqrt(2)
print("Density matrix of the equal superposition")
print(ψ @ ψ.T.conj())
print("Density matrix of the equally mixed state of |0><0| and |1><1|")
print((zero_ket @ zero_ket.T.conj()+one_ket @ one_ket.T.conj())/2)

def mixed_state(pure_state, visibility):
    #density_matrix = pure_state @ pure_state.T.conj()
    density_matrix = np.dot( pure_state, np.matrix.getH(pure_state) ) 
    maximally_mixed_state = np.eye(4)/2**2
    return visibility*density_matrix + (1-visibility)*maximally_mixed_state

ϕ = np.array([[1],[0],[0],[1]])/np.sqrt(2)
#print (array_to_latex(ϕ, pretext="$\phi$ = "))
print(f"Maximum visibility is a pure state: \n{mixed_state(ϕ, 1.0)}\n")
print(f"The state is still entangled with visibility 0.8:\n {mixed_state(ϕ, 0.8)}\n")
print(f"Entanglement is lost by 0.6:\n{mixed_state(ϕ, 0.6)}\n")
print(f"Barely any coherence remains by 0.2:\n{mixed_state(ϕ, 0.2)}\n")

import matplotlib.pyplot as plt

temperatures = [.5, 5, 2000]
energies = np.linspace(0, 20, 100)
fig, ax = plt.subplots()
for i, T in enumerate(temperatures):
    probabilities = np.exp(-energies/T)
    Z = probabilities.sum()
    probabilities /= Z
    ax.plot(energies, probabilities, linewidth=3, label = "$T_" + str(i+1)+"$="+str(T))
ax.set_xlim(0, 20)
ax.set_ylim(0, 1.2*probabilities.max())
ax.set_xticks([])
ax.set_yticks([])
ax.set_xlabel('Energy')
ax.set_ylabel('Probability')
ax.legend()

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, execute, Aer

cq = QuantumRegister(2,'code\_qubit\_') 
lq = QuantumRegister(1,'ancilla\_qubit\_') 
sb = ClassicalRegister(1,'syndrome\_bit')  
qc = QuantumCircuit(cq,lq,sb)
qc.cx(cq[0],lq[0])
qc.cx(cq[1],lq[0])
qc.measure(lq,sb)
qc.draw('mpl')

qc_init = QuantumCircuit(cq)

(qc_init+qc).draw('mpl')

counts = execute( qc_init+qc, Aer.get_backend('qasm_simulator')).result().get_counts()
print('Results:',counts)

qc_init = QuantumCircuit(cq)
qc_init.x(cq)

(qc_init+qc).draw('mpl')

counts = execute(qc_init+qc, Aer.get_backend('qasm_simulator')).result().get_counts()
print('Results:',counts)

qc_init = QuantumCircuit(cq)
qc_init.h(cq[0])
qc_init.cx(cq[0],cq[1])

(qc_init+qc).draw('mpl')

counts = execute(qc_init+qc, Aer.get_backend('qasm_simulator')).result().get_counts()
print('Results:',counts)

qc_init = QuantumCircuit(cq)
qc_init.h(cq[0])
qc_init.cx(cq[0],cq[1])
qc_init.x(cq[0])

(qc_init+qc).draw('mpl')

counts = execute(qc_init+qc, Aer.get_backend('qasm_simulator')).result().get_counts()
print('Results:',counts)