{{ message }}

Instantly share code, notes, and snippets.

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`.

This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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()   This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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!') This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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')
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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()
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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))
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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