Skip to content

Instantly share code, notes, and snippets.

@marekyggdrasil marekyggdrasil/README.md
Last active May 22, 2020

Embed
What would you like to do?

This is code repository for Quantum Adiabatic Optimization article.

Start by running instance.py that will build list-based definiton of the graph from its set-based definition, this will let you keep consistent indices between all the runs.

The run.py is the adiabatic optimization code.

Energy spectrum can be examined by running eigenenergies.py with single command line argument being one of HMVC, HMVC_, HMIS or HMIS_. I don't suggest running all the diagonalizations in single run, graphs are big and you are likely to get a segmentation fault error. This is why it is good to keep consistent indices and store your graph using lists and not sets.

You can plot figures using plot.py.

import numpy as np
import sys
import pickle
# Hamiltonians
from hamiltonians import HMVC, HMVC_, HMIS, HMIS_
# utilities
from utils import countMISViolations, extractSpinConfiguration, makePairs
# a function that diagonalizes the Hamiltonian and prints `top` best solutions
# displays spin configuration of each solution, its energy and number violations
# returns the spin configurations a list of lists
def diagonalizeHamiltonian(H, N, top, pairs, sel_mis=1):
evals, ekets = H.eigenstates(eigvals=top)
print('{0:3} {1:10} {2:20} {3:10} {4:10} {5:10}'.format(
'#', 'Energy', '|psi>', 'Size MIS', 'Size VC', 'Violations MIS'))
for i in range(top):
e = evals[i]
ket = ekets[i]
sol = extractSpinConfiguration(N, ket)
violations = countMISViolations(N, pairs, sol, sel=sel_mis)
str_ket = '|' + ''.join([str(v) for v in sol]) + '>'
print('{0:3} {1:10} {2:20} {3:10} {4:10} {5:10}'.format(
str(i),
str(np.round(e, 3)),
str_ket,
str(sol.count(sel_mis)),
str(sol.count(1 - sel_mis)),
str(violations)))
# load the instance
data = None
with open('instance.pickle', 'rb') as f:
data = pickle.load(f)
print('Kanto League map problem successfully unpickled!')
G, V, E = data
print('{0:5} {1:20}'.format('#', 'Location'))
for i, val in enumerate(V):
print('{0:5} {1:20}'.format(i, val))
print()
pairs = makePairs(V, E)
print()
print(pairs)
print(len(pairs))
print()
N = len(V)
if len(sys.argv) < 2:
raise Exception('Dude u haz got tell which Hamiltonian operatyr pleashe')
H = None
sel_mis = 0
# minimum vertex cover Hamiltonian, |1> is selected vertex and |0> is not selected vertex
if sys.argv[1] == 'HMVC':
print('Minimum Vertex Cover, with constant energy shift term')
H = HMVC(N, pairs)
# minimum vertex cover Hamiltonian, |1> is selected vertex and |0> is not selected vertex, no constant shift version
if sys.argv[1] == 'HMVC_':
print('Minimum Vertex Cover, without constant energy shift term')
H = HMVC_(N, pairs)
# maximum independent set Hamiltonian, |1> is selected vertex and |0> is not selected vertex
if sys.argv[1] == 'HMIS':
print('Maximum Independent Set, with constant energy shift term')
H = HMIS(N, pairs)
sel_mis = 1
# maximum independent set Hamiltonian, |1> is selected vertex and |0> is not selected vertex, no constant shift version
if sys.argv[1] == 'HMIS_':
print('Maximum Independent Set, without constant energy shift term')
H = HMIS_(N, pairs)
sel_mis = 1
if H is None:
raise Exception('Broh, options to use are HMVC, HMVC_, HMIS, HMIS_')
diagonalizeHamiltonian(H, N, 5, pairs, sel_mis=sel_mis)
print()
# operators
from qm import Sz, x
# utilities
from utils import osum
# minimum vertex cover Hamiltonian, |1> is selected vertex and |0> is not selected vertex
def HMVC(N, pairs):
Hq = osum([(1. - x(N, i))*(1. - x(N, j)) + (1. - x(N, j))*(1. - x(N, i)) for i, j in pairs])
Hl = 0.99*osum([x(N, i) for i in range(N)])
return Hq + Hl
# minimum vertex cover Hamiltonian, |1> is selected vertex and |0> is not selected vertex, no constant shift version
def HMVC_(N, pairs):
Hq = osum([Sz(N, i) + Sz(N, j) + Sz(N, j)*Sz(N, i) for i, j in pairs])
Hl = - 0.5*osum([Sz(N, i) for i in range(N)])
return Hq + Hl
# maximum independent set Hamiltonian, |1> is selected vertex and |0> is not selected vertex
def HMIS(N, pairs):
Hq = osum([x(N, i)*x(N, j) + x(N, j)*x(N, i) for i, j in pairs])
Hl = - 0.99*osum([x(N, i) for i in range(N)])
return Hq + Hl
# maximum independent set Hamiltonian, |1> is selected vertex and |0> is not selected vertex, no constant shift version
def HMIS_(N, pairs):
Hq = 2.*osum([-Sz(N, i) - Sz(N, j) + Sz(N, i)*Sz(N, j) for i, j in pairs])
Hl = 0.5*osum([Sz(N, i) for i in range(N)])
return Hq + Hl
import pickle
# Pokemon Kanto region map instance to try our quantum adiabatic solver
f = frozenset
V = f({
'Indigo Plateau',
'Pallet Town',
'Viridian City',
'Pewter City',
'Cinnabar Island',
'Cerulean City',
'Saffron City',
'Celadon City',
'Lavender Town',
'Vermillion City',
'Fuschia City',
'Special Region'})
E = f({
f({'Indigo Plateau', 'Viridian City'}),
f({'Pallet Town', 'Viridian City'}),
f({'Viridian City', 'Pewter City'}),
f({'Pewter City', 'Cerulean City'}),
f({'Cerulean City', 'Saffron City'}),
f({'Saffron City', 'Celadon City'}),
f({'Celadon City', 'Fuschia City'}),
f({'Fuschia City', 'Special Region'}),
f({'Fuschia City', 'Cinnabar Island'}),
f({'Vermillion City', 'Special Region'}),
f({'Cinnabar Island', 'Pallet Town'}),
f({'Saffron City', 'Lavender Town'}),
f({'Saffron City', 'Vermillion City'}),
f({'Lavender Town', 'Cerulean City'}),
f({'Lavender Town', 'Special Region'})
})
G = (V, E)
# from now on we only need the list form and we can abandon the set description
V = list(V)
E = list(E)
data = G, V, E
with open('instance.pickle', 'wb') as f:
pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)
print('Kanto League map problem successfully pickled!')
import numpy as np
# data visualization
from plotting import plotEnergy, plotFidelities, plotLegend, plotFidelitiesAcc
tmax = 25.
T = tmax*np.pi
res = 300
times = np.linspace(0., T, res)
energies = np.load('energies.npy')
plotEnergy(times, energies, 'Energy measurement', res, tmax, filename='energy.png')
runs = 20
fidelities = np.load('fidelities.npy')
titles = ['Solution #1', 'Solution #2', 'Solution #3', 'All solutions']
plotFidelities(fidelities, titles, res, runs, 25., filename='fidelities.png')
plotLegend('Legend', res, runs, 25., filename='legend.png')
plotFidelitiesAcc(fidelities, titles, res, runs, 25., filename='fidelities_acc.png')
import numpy as np
import itertools
import matplotlib.pyplot as plt
from matplotlib.figure import figaspect
# tmax is in the units of pi
def plotEnergy(times, measurements, suptitle, res, tmax, filename=None):
w, h = figaspect(.5)
fig, ax = plt.subplots(1, 1, constrained_layout=True, figsize=(w, h))
fig.suptitle(suptitle)
xticks = np.linspace(0., tmax*np.pi, 11)
xlabels = ['$'+str(np.round(v/np.pi, 2))+'\\pi$' for v in xticks]
ax.set_ylabel('$E$')
ax.set_xlabel('$t$')
ax.grid()
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels)
ax.plot(times, measurements, linestyle='solid', color='tab:blue', label='$\\left<\\psi(t)|H(t)|\\psi(t)\\right>$')
ax.legend()
if filename is not None:
fig.savefig(filename)
else:
plt.show()
plt.close()
# tmax is in the units of pi
def plotFidelities(fidelities, suptitles, res, runs, tmax, filename=None):
cmap = plt.get_cmap('plasma')
w, h = figaspect(1.)
fig, axs = plt.subplots(2, 2, constrained_layout=True, figsize=(2*w, h))
xticks = np.linspace(0., res, 11)
xlabels = ['$'+str(np.round(v/res, 2))+'$' for v in xticks]
coordinates = list(itertools.product([0, 1], repeat=2))
for j, cors in enumerate(coordinates):
row, col = cors
ax = axs[row, col]
suptitle = suptitles[j]
ax.set_ylabel('$f$')
ax.set_xlabel('$\\lambda$')
ax.grid()
ax.set_title(suptitle)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels)
ax.set_ylim(bottom=0., top=1.)
for run in range(runs):
color = cmap(float(run)/float(runs))
ax.plot(fidelities[j+1, run, :], linestyle='solid', c=color)
if filename is not None:
fig.savefig(filename)
else:
plt.show()
plt.close()
# tmax is in the units of pi
def plotFidelitiesAcc(fidelities, suptitles, res, runs, tmax, filename=None):
cmap = plt.get_cmap('plasma')
w, h = figaspect(1.)
fig, ax = plt.subplots(1, 1, constrained_layout=True, figsize=(w, h))
xticks = np.linspace(0., res, 11)
xlabels = ['$'+str(np.round(v/res, 2))+'$' for v in xticks]
coordinates = list(itertools.product([0, 1], repeat=2))
suptitle = suptitles[-1]
ax.set_ylabel('$f$')
ax.set_xlabel('$\\lambda$')
ax.grid()
ax.set_title(suptitle)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels)
ax.set_ylim(bottom=0., top=1.)
for run in range(runs):
color = cmap(float(run)/float(runs))
ax.plot(fidelities[-1, run, :], linestyle='solid', c=color)
if filename is not None:
fig.savefig(filename)
else:
plt.show()
plt.close()
def plotLegend(suptitle, res, runs, tmax, filename=None):
cmap = plt.get_cmap('plasma')
w, h = figaspect(.125)
fig, ax = plt.subplots(1, 1, constrained_layout=True, figsize=(w, h))
ax.set_xlabel('$\\tau$')
xticks = np.linspace(1., tmax, runs)
xlabels = ['$'+str(np.round(v, 2))+'\\pi$' for v in xticks]
ax.grid()
ax.set_title(suptitle)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels)
ax.set_yticks([])
ax.set_xlim([xticks[0]-1., xticks[-1]+1.])
for run in range(runs):
color = cmap(float(run)/float(runs))
ax.axvline(linewidth=16, x=xticks[run], c=color)
if filename is not None:
fig.savefig(filename)
else:
plt.show()
plt.close()
import numpy as np
# qutip tools
from qutip import tensor
# operators
from qutip import sigmax, sigmaz, qeye
# Identity operator
def Is(i): return [qeye(2) for j in range(0, i)]
# Pauli-X acting on ith degree of freedom
def Sx(N, i): return tensor(Is(i) + [sigmax()] + Is(N - i - 1))
# Pauli-Z acting on ith degree of freedom
def Sz(N, i): return tensor(Is(i) + [sigmaz()] + Is(N - i - 1))
# a quantum binary degree of freedom
def x(N, i):
return 0.5*(1.-Sz(N, i))
import numpy as np
import pickle
import itertools
# utilities
from qutip import basis, tensor, mesolve, Options, qobj_list_evaluate, expect
# operators
from qm import Sx
# Hamiltonians
from hamiltonians import HMIS_
# solver for Minimum Vertex Cover we developed previously
from solvers import solveAll
# utilities
from utils import osum, makePairs
# load the instance
data = None
with open('instance.pickle', 'rb') as f:
data = pickle.load(f)
print('Kanto League map problem successfully unpickled!')
G, V, E = data
pairs = makePairs(V, E)
N = len(V)
# quantum adiabatic computation section
tmax = 25.
T = tmax*np.pi
res = 300
Hi = osum([Sx(N, i) for i in range(N)])
Hf = HMIS_(N, pairs)
H = [[Hi, (lambda t, args: 1. - t/T)], [Hf, (lambda t, args: t/T)]]
def measureEnergy(t, psi):
# evaluate time-dependent part
H_t = qobj_list_evaluate(H, t, {})
return expect(H_t, psi)
minus = (basis(2, 0) - basis(2, 1)).unit()
psi0 = tensor([minus for i in range(N)])
times = np.linspace(0., T, res)
opts = Options(store_final_state=True)
result = mesolve(H, psi0, times, e_ops=measureEnergy, options=opts)
psif = result.final_state
measurements = result.expect
np.save('energies.npy', measurements)
def extractFidelities(N, Hf, ket, top):
# produce all possible binary configurations of length N
configurations = list(itertools.product([0, 1], repeat=N))
fidelities = np.zeros(2**N)
energies = np.zeros(2**N)
for i, configuration in enumerate(configurations):
refket = tensor([basis(2, value) for value in configuration])
fidelity = np.abs(refket.overlap(ket))**2.
fidelities[i] = fidelity
energy = expect(Hf, refket)
energies[i] = energy
# get the indices highest fidelities
best = (-fidelities).argsort()[:top]
# return those configurations and corresponding fidelities and energies
top_configurations = [configurations[c] for c in best]
top_fidelities = [fidelities[c] for c in best]
top_energies = [energies[c] for c in best]
return top_configurations, top_fidelities, top_energies
configurations, fidelities, energies = extractFidelities(N, Hf, psif, 10)
print('Fidelities')
print('{0:3} {1:20} {2:10} {3:10}'.format('#', '|psi>', 'Fidelity', 'Energy'))
for i, (configuration, fidelity, energy) in enumerate(zip(configurations, fidelities, energies)):
str_ket = '|' + ''.join([str(v) for v in configuration]) + '>'
print('{0:3} {1:20} {2:10} {3:10}'.format(str(i), str_ket, str(np.round(fidelity, 5)), str(np.round(energy, 5))))
print()
print('Top 3 fidelities add up to')
print(np.round(np.sum(fidelities[0:3]), 5))
lst, min = solveAll([], 0, 9999, V, G)
sol_kets = []
for sol in lst:
# print(sol)
sll = [1 for i in range(N)]
for loc in sol:
idx = V.index(loc)
sll[idx] = 0
ket = tensor([basis(2, v) for v in sll])
sol_kets.append(ket)
def measureEnergyFidelities(t, psi):
# evaluate time-dependent part
H_t = qobj_list_evaluate(H, t, {})
# get current energy of the system
energy = expect(H_t, psi)
# get fidelity of all the solutions
fidelities = []
for ket in sol_kets:
fidelity = np.abs(psi.overlap(ket))**2.
fidelities.append(fidelity)
return tuple([energy] + fidelities + [np.sum(fidelities)])
runs = 20
fidelities = np.zeros((5, runs, res))
for run, tmax in enumerate(np.linspace(1., 25., runs)):
T = tmax*np.pi
times = np.linspace(0., T, res)
result = mesolve(H, psi0, times, e_ops=measureEnergyFidelities, options=opts)
measurements = result.expect
for i, measurement in enumerate(measurements):
for j in range(5):
# first one is energy, next 3 are individual fidelities and last is summed fidelities
fidelities[j, run, i] = measurement[j]
np.save('fidelities.npy', fidelities)
def isVC(sol, G):
V, E = G
cover = set(sol)
# print(E, type(E))
for edge in E:
if len(list(edge.intersection(cover))) == 0:
return False
return True
def solveAll(sol, idx, minT, C, G):
lst = []
min = int(minT)
if isVC(sol, G):
size = len(sol)
if size < min:
lst = [sol]
min = size
if idx < len(C):
# try to include the next element and see what happens...
lstL, minL = solveAll(sol + [C[idx]], idx + 1, min, C, G)
# try to skip the next element and see what happens...
lstR, minR = solveAll(sol, idx + 1, min, C, G)
# whatever turned out to be better...
if minL < min:
lst, min = lstL, minL
elif minL == min:
lst += lstL
if minR < min:
lst, min = lstR, minR
elif minR == min:
lst += lstR
return lst, min
import numpy as np
from qutip import expect
# operators
from qm import Sz
# convenient sum notation that works with generic type objects
def osum(lst): return np.sum(np.array(lst, dtype=object))
# a function that counts how many Maximum Independent Set constraints have been violated
def countMISViolations(N, pairs, sol, sel=1):
violations = 0
for i, j in pairs:
if sol[i] == sel and sol[j] == sel:
violations += 1
return violations
# a function that given a ket extracts spin configuration under a list form
def extractSpinConfiguration(N, ket):
sol = []
for i in range(N):
val = expect(Sz(N, i), ket)
if val > 0:
sol += [0]
else:
sol += [1]
return sol
# simple utility function that makes list of pairs of indices
def makePairs(V, E):
pairs = []
print('{0:20} {1:20} {2:5} {3:5}'.format(
'Location A', 'Location B', '#A', '#B'))
for edge in E:
dest_i, dest_j = tuple(edge)
i, j = V.index(dest_i), V.index(dest_j)
print('{0:20} {1:20} {2:5} {3:5}'.format(dest_i, dest_j, i, j))
pairs += [(i, j)]
return pairs
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.