{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# D-Wave Examples\n",
    "\n",
    "This notebook shows how to formulate and solve QUBO problems using D-Wave's dimod software package. The problems and formulations are either from the paper *Ising formulations of many NP problems* by Andrew Lucas (2014), or from *A Tutorial on Formulating and Using QUBO Models*, by Fred Glover, Gary Kochenberger, and Yu Du (2019).\n",
    "\n",
    "## 1. Common functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dimod\n",
    "import time\n",
    "import numpy as np\n",
    "from dwave.system import DWaveSampler, AutoEmbeddingComposite\n",
    "\n",
    "def solve_qubo(Q,\n",
    "               sampler=\"CPU\", # CPU or QPU\n",
    "               k=10,\n",
    "               chain_strength=None):\n",
    "    \"\"\"\n",
    "    Given an upper triangular matrix Q of size NxN, solves the quadratic unconstrained binary\n",
    "    optimization (QUBO) problem given by\n",
    "    \n",
    "        minimize sum(x[i] * Q[i,j] * x[j]\n",
    "                     for i in range(N),\n",
    "                     for j in range(i+1, N))\n",
    "    \n",
    "    Uses dimod.SimulatedAnnealingSampler, which solves the problem k times through simulated\n",
    "    annealing (on a regular CPU). This method returns the best solution found.\n",
    "    \"\"\"\n",
    "    assert isinstance(Q, np.ndarray)\n",
    "    assert sampler in [\"CPU\", \"QPU\"]\n",
    "    n = Q.shape[0]\n",
    "    nz = len(Q[Q!=0])\n",
    "    print(\"Solving QUBO problem (%d vars, %d nz) on %s...\" % (n, nz, sampler))\n",
    "    \n",
    "    start = time.time()\n",
    "    if sampler == \"CPU\":\n",
    "        sampler = dimod.SimulatedAnnealingSampler()\n",
    "        response = sampler.sample_qubo(Q, num_reads=k)\n",
    "    else:\n",
    "        if chain_strength is None:\n",
    "            chain_strength = int(10 * np.max(np.abs(Q)))\n",
    "        sampler = AutoEmbeddingComposite(DWaveSampler(solver=dict(qpu=True)))\n",
    "        response = sampler.sample_qubo(Q, num_reads=k, chain_strength=chain_strength)\n",
    "    elapsed = time.time() - start\n",
    "    \n",
    "    print(\"Solved in %.2f seconds\" % elapsed)\n",
    "    solution = min(response.data([\"sample\", \"energy\"]), key=lambda s: s.energy)\n",
    "    return solution, response"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Number partitioning problem\n",
    "\n",
    "Given a set $S=\\left\\{s_1,\\ldots,s_n\\right\\}$ of $n$ positive integers, find a partition $(S_0, S_1)$ of $S$ minimizing the difference between $\\sum_{s \\in S_0} s$ and $\\sum_{s \\in S_1} s$.\n",
    "In this formulation, the upper triangular matrix $Q$ is given by\n",
    "$$\n",
    "    Q_{i,j} = \\begin{cases}\n",
    "        s_i (s_j - c) & \\text{if } i=j \\\\\n",
    "        2 s_i s_j & \\text{if } i>j,\n",
    "    \\end{cases}\n",
    "$$\n",
    "where $c = \\sum_{s \\in S} s$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Create random instance and Q matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import randint\n",
    "\n",
    "# Pick random integer numbers between 1 and 100\n",
    "np.random.seed(42)\n",
    "S = randint(low=1, high=11).rvs(30)\n",
    "\n",
    "# Generate Q matrix\n",
    "c, n = sum(S), len(S)\n",
    "Q = np.zeros((n, n), dtype=int)\n",
    "for i in range(n):\n",
    "    Q[i,i] = S[i] * (S[i] - c)\n",
    "    for j in range(i + 1, n):\n",
    "        Q[i,j] = 2 * S[i] * S[j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Solve on the CPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving QUBO problem (30 vars, 465 nz) on CPU...\n",
      "Solved in 18.87 seconds\n",
      "sum(S0)       86\n",
      "sum(S1)       86\n"
     ]
    }
   ],
   "source": [
    "# Solve problem\n",
    "solution, response = solve_qubo(Q)\n",
    "\n",
    "# Display solution\n",
    "S0 = [S[i] for (i, xi) in solution.sample.items() if xi >  0.5]\n",
    "S1 = [S[i] for (i, xi) in solution.sample.items() if xi <= 0.5]\n",
    "print(\"sum(S0) %8d\" % sum(S0))\n",
    "print(\"sum(S1) %8d\" % sum(S1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Solve on the QPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving QUBO problem (30 vars, 465 nz) on QPU...\n",
      "Solved in 17.30 seconds\n",
      "sum(S0)       86\n",
      "sum(S1)       86\n"
     ]
    }
   ],
   "source": [
    "# Solve problem\n",
    "solution, qresponse = solve_qubo(Q, sampler=\"QPU\")\n",
    "\n",
    "# Display solution\n",
    "S0 = [S[i] for (i, xi) in solution.sample.items() if xi >  0.5]\n",
    "S1 = [S[i] for (i, xi) in solution.sample.items() if xi <= 0.5]\n",
    "print(\"sum(S0) %8d\" % sum(S0))\n",
    "print(\"sum(S1) %8d\" % sum(S1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Inspect QPU solution and embedding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'http://127.0.0.1:18000/?problemId=72672b81-b943-4dc6-ade7-24b0adce7c5d'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import dwave.inspector\n",
    "dwave.inspector.show(qresponse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}