qaoa-code-sample.py

#! /usr/bin/env python3

from qiskit import QuantumCircuit, transpile
from qiskit.providers.aer import *
from qiskit.providers.aer.extensions import *
import qiskit
import random
import numpy as np
import sys
def E3LIN2_tuples(n, D, max_iterations=10000):
    counts = {} # count the number of terms variable n appears in
    for k in range(n):
        counts[k] = 0

    vars_remaining = set(range(n))

    terms = set()
    count = 0
    while len(vars_remaining) >= 3:
        count += 1
        if count > max_iterations:
            return None
        
        x1, x2, x3 = random.sample(list(vars_remaining),3)
        x1,x2,x3 = sorted((x1,x2,x3))
        
        if not (x1, x2, x3) in terms:        
            terms.add((x1,x2,x3))
            counts[x1] += 1
            counts[x2] += 1
            counts[x3] += 1
            
            if counts[x1] == 4:
                vars_remaining.remove(x1)
            if counts[x2]== 4:
                vars_remaining.remove(x2)
            if counts[x3] == 4:
                vars_remaining.remove(x3)

            
    return terms, counts, vars_remaining

def E3LIN2(n, D):
    #randomly generate E3LIN2 instances until we find one that meets the constraints given in the paper
    terms = None
    while terms == None:
        ret = E3LIN2_tuples(n,D)
        if type(ret) == tuple:        
            terms, counts, vars_remaining = ret
            if len(vars_remaining) != 1:
                terms = None

    polynomial = {}
    for term in terms:
        polynomial[term] = random.choice([-1, 1])

    return polynomial

def C(polynomial, z):
    #see page 12 of paper
    val = 0.

    for term in polynomial:
        u,v,w = term
        zu, zv, zw = int(z[u]),int(z[v]),int(z[w])        
        val += polynomial[(u,v,w)]*(1-2*zu)*(1-2*zv)*(1-2*zw)/2.        
    return val

def QAOA_circuit_U(num_qubits, gamma, polynomial):
    circ = QuantumCircuit(num_qubits, num_qubits)

    for q in range(num_qubits):
        circ.h(q)
        
    for term in polynomial:
        circ.cx(control_qubit=term[1],target_qubit=term[0])
        circ.cx(control_qubit=term[2],target_qubit=term[0])
        if polynomial[term] < 0:
            circ.x(term[0])
        circ.p(gamma,term[0])
        if polynomial[term] < 0:
            circ.x(term[0])
        circ.cx(control_qubit=term[1],target_qubit=term[0])
        circ.cx(control_qubit=term[2],target_qubit=term[0])

    for q in range(num_qubits):
        circ.h(q)
        circ.s(q)
        circ.h(q)
                    
    return circ


def QAOA(n, polynomial):
    time_s = 0.
    num_shots = 40000
    mixing_time = 10000
    threads = 10
    error = 0.05
    gammas = np.array([0.1])*np.pi
    expectation_vals = np.zeros_like(gammas)
    
    for i, gamma in enumerate(gammas):
        circ = QAOA_circuit_U(n, gamma, polynomial)
        circ.measure(list(range(n)), list(range(n)))        
        extended_stabilizer_simulator = AerSimulator(method='extended_stabilizer',
                                                     max_parallel_threads=threads,
                                                     extended_stabilizer_sampling_method="metropolis",
                                                     extended_stabilizer_metropolis_mixing_time=mixing_time)
        
        t_circ = transpile(circ, extended_stabilizer_simulator)        
        result = extended_stabilizer_simulator.run(t_circ, shots=num_shots,extended_stabilizer_approximation_error=error).result()
        counts = result.get_counts(0)
        for bitstring in counts:
            reversed_string = bitstring[::-1] # qiskit gives us the bitstrings "backwards" with string[0] being the outcome of the measurement on qubit[-1]
            expectation_vals[i] += C(polynomial, reversed_string) * (counts[bitstring]/float(num_shots))

        print(i, gamma, expectation_vals[i])    

        
if __name__ == "__main__":
    print("python version:", sys.version_info)
    print("qiskit version:", qiskit.__version__)
    print("aer version:", qiskit.providers.aer.__version__)
    n=50
    D=4
    #polynomial = E3LIN2(n, D)
    #for reporducability we use a fixed polynomial, uncomment the code above and comment out the fixed polynomial if you want to use a randomly generated one
    polynomial = {(17, 42, 46): -1, (33, 35, 40): 1, (2, 8, 15): -1, (13, 36, 48): 1, (10, 25, 35): -1, (11, 28, 49): -1, (13, 14, 48): 1, (5, 6, 26): -1, (7, 26, 31): -1, (21, 34, 37): 1, (8, 24, 45): -1, (19, 26, 31): 1, (5, 27, 33): 1, (3, 23, 36): -1, (0, 18, 28): 1, (22, 27, 29): 1, (31, 38, 45): -1, (3, 15, 37): -1, (10, 13, 39): -1, (4, 30, 47): -1, (11, 38, 45): 1, (5, 42, 43): 1, (2, 6, 29): 1, (20, 37, 42): 1, (3, 12, 32): -1, (19, 31, 39): -1, (12, 15, 41): -1, (9, 33, 41): -1, (18, 33, 37): 1, (34, 45, 48): 1, (1, 24, 28): 1, (13, 17, 49): -1, (4, 8, 40): 1, (1, 7, 30): 1, (25, 38, 39): -1, (27, 36, 44): -1, (10, 20, 43): 1, (18, 21, 47): 1, (22, 24, 41): 1, (1, 24, 49): 1, (15, 19, 43): -1, (23, 30, 48): 1, (3, 10, 25): 1, (20, 40, 49): -1, (7, 22, 38): -1, (16, 20, 26): 1, (21, 28, 32): -1, (1, 17, 18): -1, (21, 44, 46): 1, (6, 11, 35): -1, (5, 7, 23): 1, (12, 22, 27): -1, (32, 41, 44): -1, (2, 14, 29): 1, (0, 14, 32): -1, (2, 34, 47): -1, (39, 43, 46): -1, (9, 34, 35): 1, (4, 16, 19): 1, (6, 9, 40): 1, (4, 16, 25): 1, (11, 12, 23): 1, (17, 36, 42): 1, (8, 14, 29): 1, (16, 44, 47): -1, (9, 30, 46): 1}
    
    QAOA(n, polynomial)