optimization/stable_set.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## _*Using Qiskit Aqua for stable-set problems*_\n",
    "\n",
    "This Qiskit Aqua Optimization notebook demonstrates how to use the VQE algorithm to compute the maximum stable set of a given graph.  \n",
    "\n",
    "The problem is defined as follows. Given a graph $G = (V,E)$, we want to compute $S \\subseteq V$ such that there do not exist $i, j \\in S : (i, j) \\in E$, and $|S|$ is maximized. In other words, we are looking for a maximum cardinality set of mutually non-adjacent vertices.\n",
    "\n",
    "The graph provided as an input is used first to generate an Ising Hamiltonian, which is then passed as an input to VQE.  As a reference, this notebook also computes the maximum stable set using the NumPyMinimumEigensolver classical algorithm and the solver embedded in the commercial CPLEX product (if it is available in the system and the user has followed the necessary configuration steps in order for Qiskit Aqua to find it).  Please refer to the Qiskit Aqua Optimization documentation for installation and configuration details for CPLEX."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from qiskit import BasicAer\n",
    "from qiskit.optimization.applications.ising import stable_set\n",
    "from qiskit.aqua import QuantumInstance\n",
    "from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE\n",
    "from qiskit.aqua.components.optimizers import L_BFGS_B\n",
    "from qiskit.circuit.library import TwoLocal\n",
    "from qiskit.optimization.applications.ising.common import parse_gset_format, random_graph, sample_most_likely"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here an Operator instance is created for our Hamiltonian. In this case the Paulis are from an Ising Hamiltonian of the maximum stable set problem (expressed in minimization form). We load a small instance of the maximum stable set problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = parse_gset_format('sample.maxcut')\n",
    "qubitOp, offset = stable_set.get_operator(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also offer a function to generate a random graph as a input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.  8. -9.  0.  7.]\n",
      " [ 8.  0.  9. -8. -1.]\n",
      " [-9.  9.  0.  0. -5.]\n",
      " [ 0. -8.  0.  0.  0.]\n",
      " [ 7. -1. -5.  0.  0.]]\n"
     ]
    }
   ],
   "source": [
    "if True:\n",
    "    np.random.seed(8123179)\n",
    "    w = random_graph(5, edge_prob=0.5)\n",
    "    qubitOp, offset = stable_set.get_operator(w)\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we test for the presence of algorithms we want to use in this notebook. If Aqua is installed correctly `NumPyMinimumEigensolver` and `VQE` will always be found. `ClassicalCPLEX` is dependent on CPLEX being installed (see introduction above). CPLEX is *not required* but if installed then this notebook will demonstrate the `ClassicalCPLEX` algorithm , that uses CPLEX, to compute stable set as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['NumPyMinimumEigensolver', 'ClassicalCPLEX', 'VQE']\n"
     ]
    }
   ],
   "source": [
    "to_be_tested_algos = ['NumPyMinimumEigensolver', 'ClassicalCPLEX', 'VQE']\n",
    "print(to_be_tested_algos)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the Operator without regard to how it was created. First we need to prepare the configuration params to invoke the algorithm. Here we will use the NumPyMinimumEigensolver first to return the smallest eigenvalue. Backend is not required since this is computed classically not using quantum computation. We then add in the qubitOp Operator in dictionary format. The result is a dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "energy: -29.5\n",
      "stable set objective: -25.0\n",
      "solution: [0 0 1 1 1]\n",
      "solution objective and feasibility: (3, False)\n"
     ]
    }
   ],
   "source": [
    "result = NumPyMinimumEigensolver(qubitOp).run()\n",
    "x = sample_most_likely(result.eigenstate)\n",
    "print('energy:', result.eigenvalue.real)\n",
    "print('stable set objective:', result.eigenvalue.real + offset)\n",
    "print('solution:', stable_set.get_graph_solution(x))\n",
    "print('solution objective and feasibility:', stable_set.stable_set_value(x, w))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note*: If CPLEX installed then the Aqua ClassicalCPLEX algorithm will be able to be used. If not, then solving this problem using this particular algorithm will simply be skipped.\n",
    "\n",
    "We change the configuration parameters to solve it with the CPLEX backend. The CPLEX backend can deal with a particular type of Hamiltonian called Ising Hamiltonian, which consists of only Pauli Z at most second order and can be used for combinatorial optimization problems that can be formulated as quadratic unconstrained binary optimization problems, such as the stable set problem. Note that we may obtain a different solution - but if the objective value is the same as above, the solution will be optimal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Version identifier: 12.10.0.0 | 2019-11-27 | 843d4de\n",
      "CPXPARAM_Read_DataCheck                          1\n",
      "CPXPARAM_Threads                                 1\n",
      "CPXPARAM_MIP_Display                             0\n",
      "CPXPARAM_TimeLimit                               600\n",
      "CPXPARAM_MIP_Tolerances_MIPGap                   0\n",
      "CPXPARAM_MIP_Tolerances_Integrality              0\n",
      "energy: -29.5\n",
      "time: 0.02447683199999995\n",
      "stable set objective: -25.0\n",
      "solution: [0 0 1 1 1]\n",
      "solution objective and feasibility: (3, False)\n"
     ]
    }
   ],
   "source": [
    "try:\n",
    "    from qiskit.aqua.algorithms import ClassicalCPLEX\n",
    "    result = ClassicalCPLEX(qubitOp, display=0).run()\n",
    "\n",
    "    x_dict = result['x_sol']\n",
    "    print('energy:', result['energy'])\n",
    "    print('time:', result['eval_time'])\n",
    "    print('stable set objective:', result['energy'] + offset)\n",
    "    x = np.array([x_dict[i] for i in sorted(x_dict.keys())])\n",
    "    print('solution:', stable_set.get_graph_solution(x))\n",
    "    print('solution objective and feasibility:', stable_set.stable_set_value(x, w))\n",
    "except Exception as ex:\n",
    "    print(str(ex))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want VQE and so change it and add its other configuration parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "energy: -29.49999999994602\n",
      "time: 18.13162088394165\n",
      "stable set objective: -24.99999999994602\n",
      "solution: [0. 0. 1. 1. 1.]\n",
      "solution objective and feasibility: (3.0, False)\n"
     ]
    }
   ],
   "source": [
    "vqe = VQE(qubitOp,\n",
    "          TwoLocal(qubitOp.num_qubits, ['ry', 'rz'], 'cz', reps=3, entanglement='linear'),\n",
    "          L_BFGS_B(maxfun=6000))\n",
    "result = vqe.run(QuantumInstance(BasicAer.get_backend('statevector_simulator')))\n",
    "\n",
    "x = sample_most_likely(result.eigenstate)\n",
    "print('energy:', result.eigenvalue.real)\n",
    "print('time:', result.optimizer_time)\n",
    "print('stable set objective:', result.eigenvalue.real + offset)\n",
    "print('solution:', stable_set.get_graph_solution(x))\n",
    "print('solution objective and feasibility:', stable_set.stable_set_value(x, w))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}