Instantly share code, notes, and snippets.

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

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 == '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 == '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 == '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 == '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-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 +=  else: sol +=  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