Merge pull request #90 from AsierO/master

QAOA for Ising models with external fields

examples/IsingSolver.ipynb

@@ -0,0 +1,151 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This Notebook is a short example of how to use the Ising solver implemented using the QAOA algorithm. We start by declaring the import of the ising function."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from grove.ising.ising_qaoa import ising\n",
+    "from mock import patch"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This code finds the global minima of an Ising model with external fields of the form\n",
+    "$$f(x)= \\Sigma_i h_i x_i + \\Sigma_{i,j} J_{i,j} x_i x_j.$$\n",
+    "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\\}$.\n",
+    "\n",
+    "In order to assert the correctness of the code we will find the minima of the following Ising model\n",
+    "$$f(x)=x_0+x_1-x_2+x_3-2 x_0 x_1 +3 x_2 x_3.$$\n",
+    "Which corresponds to $x_{min}=[-1, -1, 1, -1]$ in numerical order, with a minimum value of $f(x_{min})=-9$. \n",
+    "\n",
+    "This Ising code runs on quantum hardware, which means that we need to specify a connection to a QVM or QPU. Due to the absence of a real connection in this notebook, we will mock out the response to correspond to the expected value. In order to run this notebook on a QVM or QPU, replace cxn with a valid PyQuil connection object.\n",
+    "\n",
+    "\n",
+    "\n",
+    "  "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with patch(\"pyquil.api.SyncConnection\") as cxn:\n",
+    "        cxn.run_and_measure.return_value = [[1,1,0,1]]\n",
+    "        cxn.expectation.return_value = [-0.4893891813015294, 0.8876822987380573, -0.4893891813015292, -0.9333372094534063, -0.9859245403423198, 0.9333372094534065]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "J = {(0, 1): -2, (2, 3): 3}\n",
+    "h = [1, 1, -1, 1]\n",
+    "\n",
+    "solution, min_energy, circuit = ising(h, J, connection=cxn)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "It is also possible to specify the Trotterization order for the QAOA algorithm used to implement the Ising model. By default this value is equal to double the number of variables. It is also possible to change the verbosity of the function, which is True by default. There are more advanced parameters that can be specified and are not described here.  "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "solution_2, min_energy_2, circuit_2 = ising(h, J, num_steps=9, verbose=False, connection=cxn)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For large Ising problems, or those with many and close suboptimal minima, it is possible for the code to not return the global minima. Increasing the number of steps can solve this problem.\n",
+    "\n",
+    "Finally, we will check if the correct bitstring was found, corresponding to the global minima, in both runs."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert solution == [-1, -1, 1, -1], \"Found bitstring for first run does not correspond to global minima\"\n",
+    "print(\"Energy for first run solution\", min_energy)\n",
+    "assert solution_2 == [-1, -1, 1, -1], \"Found bitstring for second run does not correspond to global minima\"\n",
+    "print(\"Energy for second run solution\", min_energy_2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If the assertions succeeded, and the energy was equal to $-9$, we have found the correct solution for both runs. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    ""
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 2",
+   "language": "python",
+   "name": "python2"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2.0
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}

grove/ising/__init__.py

@@ -0,0 +1 @@
+# Ising code init

grove/ising/ising_qaoa.py

@@ -0,0 +1,129 @@
+
+"""
+Finding the minimum energy for an Ising problem by QAOA.
+"""
+import pyquil.api as api
+from grove.pyqaoa.qaoa import QAOA
+from pyquil.paulis import PauliSum, PauliTerm
+from scipy.optimize import minimize
+import numpy as np
+
+CXN = api.SyncConnection()
+
+
+def energy_value(h, J, sol):
+    """
+    Obtain energy of an Ising solution for a given Ising problem (h,J).
+
+    :param h: External magnectic term of the Ising problem. List.
+    :param J: Interaction term of the Ising problem. Dictionary.
+    :param sol: Ising solution. List.
+    :return: Energy of the Ising string.
+    :rtype: Integer or float.
+
+    """
+    ener_ising = 0
+    for elm in J.keys():
+        if elm[0] == elm[1]:
+            raise TypeError("""Interaction term must connect two different variables""")
+        else:
+            ener_ising += J[elm] * int(sol[elm[0]]) * int(sol[elm[1]])
+    for i in range(len(h)):
+        ener_ising += h[i] * int(sol[i])
+    return ener_ising
+
+
+def print_fun(x):
+    print(x)
+
+
+def ising_trans(x):
+    # Transformation to Ising notation
+    if x == 1:
+        return -1
+    else:
+        return 1
+
+
+def ising(h, J, num_steps=0, verbose=True, rand_seed=None, connection=None, samples=None,
+          initial_beta=None, initial_gamma=None, minimizer_kwargs=None,
+          vqe_option=None):
+    """
+    Ising set up method
+
+    :param h: External magnectic term of the Ising problem. List.
+    :param J: Interaction term of the Ising problem. Dictionary.
+    :param num_steps: (Optional.Default=2 * len(h)) Trotterization order for the
+                  QAOA algorithm.
+    :param verbose: (Optional.Default=True) Verbosity of the code.
+    :param rand_seed: (Optional. Default=None) random seed when beta and
+                      gamma angles are not provided.
+    :param connection: (Optional) connection to the QVM. Default is None.
+    :param samples: (Optional. Default=None) VQE option. Number of samples
+                    (circuit preparation and measurement) to use in operator
+                    averaging.
+    :param initial_beta: (Optional. Default=None) Initial guess for beta
+                         parameters.
+    :param initial_gamma: (Optional. Default=None) Initial guess for gamma
+                          parameters.
+    :param minimizer_kwargs: (Optional. Default=None). Minimizer optional
+                             arguments.  If None set to
+                             {'method': 'Nelder-Mead',
+                             'options': {'ftol': 1.0e-2, 'xtol': 1.0e-2,
+                                        'disp': False}
+    :param vqe_option: (Optional. Default=None). VQE optional
+                             arguments.  If None set to
+                       vqe_option = {'disp': print_fun, 'return_all': True,
+                       'samples': samples}
+    :return: Most frequent Ising string, Energy of the Ising string, Circuit used to obtain result.
+    :rtype: List, Integer or float, 'pyquil.quil.Program'.
+
+    """
+    if num_steps == 0:
+        num_steps = 2 * len(h)
+
+    n_nodes = len(h)
+
+    cost_operators = []
+    driver_operators = []
+    for i, j in J.keys():
+        cost_operators.append(PauliSum([PauliTerm("Z", i, J[(i, j)]) * PauliTerm("Z", j)]))
+
+    for i in range(n_nodes):
+        cost_operators.append(PauliSum([PauliTerm("Z", i, h[i])]))
+
+    for i in range(n_nodes):
+        driver_operators.append(PauliSum([PauliTerm("X", i, -1.0)]))
+
+    if connection is None:
+        connection = CXN
+
+    if minimizer_kwargs is None:
+        minimizer_kwargs = {'method': 'Nelder-Mead',
+                            'options': {'ftol': 1.0e-2, 'xtol': 1.0e-2,
+                                        'disp': False}}
+    if vqe_option is None:
+        vqe_option = {'disp': print_fun, 'return_all': True,
+                      'samples': samples}
+
+    if not verbose:
+        vqe_option['disp'] = None
+
+    qaoa_inst = QAOA(connection, n_nodes, steps=num_steps, cost_ham=cost_operators,
+                     ref_hamiltonian=driver_operators, store_basis=True,
+                     rand_seed=rand_seed,
+                     init_betas=initial_beta,
+                     init_gammas=initial_gamma,
+                     minimizer=minimize,
+                     minimizer_kwargs=minimizer_kwargs,
+                     vqe_options=vqe_option)
+
+    betas, gammas = qaoa_inst.get_angles()
+    most_freq_string, sampling_results = qaoa_inst.get_string(
+        betas, gammas)
+    most_freq_string_ising = [ising_trans(it) for it in most_freq_string]
+    energy_ising = energy_value(h, J, most_freq_string_ising)
+    param_prog = qaoa_inst.get_parameterized_program()
+    circuit = param_prog(np.hstack((betas, gammas)))
+
+    return most_freq_string_ising, energy_ising, circuit

grove/tests/ising/test_ising.py

@@ -0,0 +1,28 @@
+from grove.ising.ising_qaoa import ising
+from grove.ising.ising_qaoa import energy_value
+import numpy as np
+from mock import patch
+
+
+def test_energy_value():
+    J = {(0, 1): 2.3}
+    h = [-2.4, 5.2]
+    sol = [1, -1]
+    ener_ising = energy_value(h, J, sol)
+
+    assert(np.isclose(ener_ising, -9.9))
+
+
+def test_ising_mock():
+    with patch("pyquil.api.SyncConnection") as cxn:
+        # Mock the response
+        cxn.run_and_measure.return_value = [[1, 1, 0, 1]]
+        cxn.expectation.return_value = [-0.4893891813015294, 0.8876822987380573, -0.4893891813015292, -0.9333372094534063, -0.9859245403423198, 0.9333372094534065]
+
+    J = {(0, 1): -2, (2, 3): 3}
+    h = [1, 1, -1, 1]
+    p = 1
+    most_freq_string_ising, energy_ising, circuit = ising(h, J, num_steps=p, vqe_option=None, connection=cxn)
+
+    assert most_freq_string_ising == [-1, -1, 1, -1]
+    assert energy_ising == -9