Code for https://mareknarozniak.com/2020/05/22/quantum-adiabatic-optimization/ tutorial article

README.md

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.

eigenenergies.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()

energy.png

fidelities.png

fidelities_acc.png

hamiltonians.py

# 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

instance.py

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!')

legend.png

plot.py

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')

plotting.py

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()

qm.py

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))

run.py

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)

solvers.py

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

utils.py

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

Hi @Mocha1208
I had the same error as you. I changed the osum function as follows: def osum(lst): return sum(lst), then the mesolve works. Not sure if this works for you as well, but you can try.