<jupyter_start><jupyter_text>Implementation of a Quantum Approximation Optimization Algorithm (QAOA).This example is based on the exampled provided by Rigetti [Rigetti](https://github.com/rigetticomputing/grove/blob/master/examples/IsingSolver.ipynb).This code finds the global minima of an Ising model with external fields of the form$$f(x)= \Sigma_i h_i x_i + \Sigma_{i,j} J_{i,j} x_i x_j.$$Two adjacent sites $i,j$ have an interaction equal to $J_{i,j}$. There is also an external magnetic field $h_i$ that affects each individual spin. The discrete variables take the values $x_i \in \{+1,-1\}$.The reference hamiltonian will be $$H_{b}=\sum_{i=0}^{N-1}\sigma^X_i$$with N the number of qubits for the problem (in this small case, 4). This state has the ground state define for a Walsh-Hadamard state$$|\Psi(0)\rangle = |+>_{N-1}\otimes|+>_{N-2}\otimes \dots \otimes |+>_{0} = \frac{1}{\sqrt{2^N}}\sum_{i=0}^{2^N-1}|i\rangle$$ So, the time evolution will be applied starting from this ground state.You will find the minima of the following Ising model$$f(x)=x_0+x_1-x_2+x_3-2 x_0 x_1 +3 x_2 x_3.$$Which corresponds to $x_{min}=[-1, -1, 1, -1]$ in numerical order, with a minimum value of $f(x_{min})=-9$.  Remember that, as Variational Quantum Eigensolver (VQE), this is a hybrid algorithm. Part of the code is executed in the CPU (the optimisation) and part in the QPU (the calculus of the expected values of the Hamiltonian).   1. Import the needed packagesImport:1. ProjecQ Simulator2. Operations to be used. Because this is mainly a time evolution, the most important is TimeEvolution3. The optimization function from Scipy.<jupyter_code>import projectq
from projectq.backends import Simulator
from projectq.ops import All, Measure, QubitOperator, TimeEvolution,H
from scipy.optimize import minimize<jupyter_output><empty_output><jupyter_text>2. Define the functions for the optimization This functions will calculate the expectation value for a hamiltonian *H* after applying the Time Evolution of the Hamiltonian composed by the reference and cost hamiltonians, for a selected number of times. However, the time evolution of each hamiltonian is executed for a different *time* $\theta_i$ which are the optimisation parameters.<jupyter_code>def Expectation_H(theta,nqubits, steps,base_ham,cost_ham,eval_ham):
    """
    Args:
        theta (float): array of variational parameters for ansatz wavefunction
        nqubits: number of qubits to use for this Hamiltonian
        steps: number of times that the time evolution is repeated
        base_ham: the base hamiltonian
        cost_ham: the cost hamiltonian
        hamiltonian (QubitOperator): Hamiltonian to evaluate
    Returns:
        energy of the wavefunction for parameters
    """
    # Create a ProjectQ compiler with a simulator as a backend


    eng = projectq.MainEngine(backend=Simulator(gate_fusion=True, rnd_seed=1000))
    wavefunction = eng.allocate_qureg(nqubits)

    # Initialize to Walsh-Hadamard state
    All(H) | wavefunction
    
    #Get the parameters from the optimizer
    alfa=theta[:steps]
    gamma=theta[steps:]
    
    #Apply the "time evolution"  a number of times (steps)
    for i in range(steps):
        TimeEvolution(gamma[i], cost_ham) | wavefunction
        TimeEvolution(alfa[i], base_ham) | wavefunction
        
    # flush all gates
    eng.flush()
    # Calculate the energy.
    # The simulator can directly return expectation values, while on a
    # real quantum devices one would have to measure each term of the
    # Hamiltonian.
    energy = eng.backend.get_expectation_value(eval_ham, wavefunction)

    # Measure in order to return to return to a classical state
    # (as otherwise the simulator will give an error)
    All(Measure) | wavefunction
    
    
    del eng
    
    return energy<jupyter_output><empty_output><jupyter_text>Helper function to compose the real hamiltonians from their terms.<jupyter_code>def compose_ham(Base,hamiltonian):
    
    H_o=None
    for i in hamiltonian:
        if (H_o is None):
            H_o=i
        else:
            H_o=H_o+i
    H_b=None
    for i in Base:
        if (H_b is None):
            H_b=i
        else:
            H_b=H_b+i
    return H_b,H_o<jupyter_output><empty_output><jupyter_text>This function returns the most probable state (which is is solution)<jupyter_code>def State_H(theta,nqubits, steps,base_ham,cost_ham):
    
    """
    Args:
        theta (float): variational parameter for ansatz wavefunction
        hamiltonian (QubitOperator): Hamiltonian of the system
    Returns:
        energy of the wavefunction for parameter theta
    """
    import numpy as np
    H_b,H_o=compose_ham(base_ham,cost_ham)
    # Create a ProjectQ compiler with a simulator as a backend
    from projectq.backends import Simulator
    eng = projectq.MainEngine(backend=Simulator(gate_fusion=True, rnd_seed=1000))
    wavefunction = eng.allocate_qureg(nqubits)

    # Initialize to Walsh-Hadamard state
    All(H) | wavefunction
    #print("Theta:",theta)
    alfa=theta[:steps]
    gamma=theta[steps:]
    #print(steps)
    for i in range(steps):
        TimeEvolution(gamma[i], H_o) | wavefunction
        TimeEvolution(alfa[i], H_b) | wavefunction
        
    # flush all gates
    eng.flush()
    maxp=0.0
    maxstate=None
    for i in range(2**nqubits):
        bits=np.binary_repr(i,width=len(wavefunction))
        statep=eng.backend.get_probability(bits[-1::-1],wavefunction)
        if (maxp < statep):
            maxstate=bits[-1::-1]
            maxp=statep
            
    All(Measure) | wavefunction
    eng.flush()
    
    del eng
    
    return maxstate,maxp<jupyter_output><empty_output><jupyter_text>Function to calculate the expectation values of each term of the hamiltonian. This step can be executed in parallel<jupyter_code>def variational_quantum_eigensolver(theta,nqubits,steps,Base,hamiltonian):
    #print("Theta:",theta)
    vqe=0.
    
    H_b,H_o=compose_ham(Base,hamiltonian) 
    
    for i in hamiltonian:
        vqe+=Expectation_H(theta,nqubits,steps,H_b,H_o,i)
    print("VQE:",vqe)    
    return vqe<jupyter_output><empty_output><jupyter_text>4. OptimizeThis is the main part. Starting from a defined Ising Hamiltonian, find the result using an optimizer  The input for the code in the default mode corresponds simply to the parameters $h_i$ and $J_{i,j}$, that we specify as a list in numerical order and a dictionary. The code returns the bitstring of the minima, the minimum value, and the QAOA quantum circuit used to obtain that result.<jupyter_code>J = {(0, 1): -2, (2, 3): 3}
h = [1, 1, -1, 1]
num_steps=10
import numpy as np
#
# if the number os steps is 0, select them as twice the number of qubits
if num_steps == 0:
    num_steps = 2 * len(h)

nqubits = len(h)
hamiltonian_o=[]
hamiltonian_b=[]
for i, j in J.keys():
    hamiltonian_o.append(QubitOperator("Z%d Z%d"%(i,j),J[(i, j)]))

for i in range(nqubits):
    hamiltonian_o.append(QubitOperator("Z%d"%i,h[i]))

for i in range(nqubits):
    hamiltonian_b.append(QubitOperator("X%d"%i,-1.0))

betas = np.random.uniform(0, np.pi, num_steps)[::-1]
gammas = np.random.uniform(0, 2*np.pi, num_steps)
theta_0=np.zeros(2*num_steps)
theta_0[0:num_steps]=betas
theta_0[num_steps:]=gammas

minimum = minimize(variational_quantum_eigensolver,theta_0,args=(nqubits,num_steps,hamiltonian_b,hamiltonian_o),
                  method='Nelder-Mead',options= {'disp': True,'ftol': 1.0e-2,'xtol': 1.0e-2,'maxiter':20})<jupyter_output>VQE: -3.1267108885017905
VQE: -3.0746127151692515
VQE: -2.4703078260129576
VQE: -2.992830279033867
VQE: -3.328404601259681
VQE: -3.0947854574467586
VQE: -2.9046451732570784
VQE: -3.671034058669843
VQE: -2.8959402609395615
VQE: -3.021877167297663
VQE: -3.2702033489490994
VQE: -1.2537536430481186
VQE: -4.732157941905806
VQE: -1.9466821811644195
VQE: -2.878178215938844
VQE: -3.0799305618827497
VQE: -2.910286964127516
VQE: -1.8677697906524502
VQE: -3.0074331998096366
VQE: -3.960728630494119
VQE: -3.7770706378717254
VQE: -4.514548657387155
VQE: -2.6776213636174955
VQE: -4.514024405976839
VQE: -4.918323312752894
VQE: -5.817911089057661
VQE: -2.2842284936619692
VQE: -3.5376680087503636
VQE: -3.668135887828318
VQE: -4.902031488366625
VQE: -4.441174337021387
VQE: -3.2189336009442426
VQE: -4.988852793021318
VQE: -4.572076042091683
VQE: -5.406066757494104
VQE: -5.430014975251307
VQE: -5.361876257918651
VQE: -5.939176668760072
VQE: -6.47308100271566
VQE: -6.157242214265662
VQE: -3.6792837753790777[...]<jupyter_text>And calculate now the most probable state<jupyter_code>maxstate,maxp=State_H(minimum.x,nqubits,num_steps,hamiltonian_b,hamiltonian_o)<jupyter_output><empty_output><jupyter_text>Ok. This is the end. Show the results<jupyter_code>print([(-1 if int(i)==1 else 1) for i in maxstate], " with probability %.2f"%maxp)<jupyter_output>[-1, -1, 1, -1]  with probability 0.51
{'fraction_non_alphanumeric': 0.10265646105357946, 'fraction_numerical': 0.08667266996848266, 'mean_word_length': 4.054038680318544, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 34, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '531976', 'license': 'Apache-2.0', 'n_tokens_mistral': 3614, 'n_tokens_neox': 3155, 'n_words': 864, 'path': 'Notebooks/QAOA.ipynb', 'repo_name': 'gomeztato/QuantumCourse', 'size': '8884'}