qubo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10のメジャーな組合せ最適化問題におけるQUBO行列生成コードの実装\n",
    "\n",
    "- 参考：[A Tutorial on Formulating and Using QUBO Models](http://meta-analytics.net/references/QUBO%20Tutorial%20%20Version%201-4.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import annotations\n",
    "from itertools import combinations\n",
    "from dwave_qbsolv import QBSolv\n",
    "from typing import Set, Tuple, Container, Dict, Sequence\n",
    "from dataclasses import dataclass\n",
    "import networkx as nx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- qubo行列は`np.ndarray`で表すことにする\n",
    "- しかしdwave_qbsolvに渡すには`dict`にする必要があるので，helper functionを用意しておく"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_dict(q: np.ndarray) -> Dict[Tuple[int, int], float]:\n",
    "    \"\"\"Returns q matrix in the form of dict from np.ndarray\n",
    "    \n",
    "    The dict can be directly passed to dwave_qbsolv.QBSolv\n",
    "    \"\"\"\n",
    "    assert len(q.shape) == 2 and  q.shape[0] == q.shape[1]\n",
    "    return dict(zip(zip(*np.nonzero(q)), q[q!=0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- test code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_to_dict():\n",
    "    q = np.array([[1,2],\n",
    "                  [2,4]])\n",
    "    assert to_dict(q) == {(0,0):1, (0,1):2, (1,0):2, (1,1):4}\n",
    "test_to_dict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `np.ndarray`のQUBO行列から解を得るfunctionも用意しておく"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qbsolv_from_ndarray(q: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:\n",
    "    \"\"\"Returns lowest energy samples from dwave qbsolv.(for test util)\"\"\"\n",
    "    res = QBSolv().sample_qubo(to_dict(q)).record\n",
    "    i = res.energy == res.energy.min()\n",
    "    return res.sample[i], res.energy[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  The Number Partitioning Problem\n",
    "\n",
    "- [Partition problem \\- Wikipedia](https://en.wikipedia.org/wiki/Partition_problem)\n",
    "- 数の集合を2つに分割したとき，それぞれの和が一致するような分け方はあるか？\n",
    "- 集合の要素にqubitを$x_{i}$を割り当てる\n",
    "    - 分割後の集合kに入る場合は$x_{i}=k$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubo_number_partition(number_set: Container[float]) -> np.ndarray:\n",
    "    \"\"\"Return Q matrix of number partitioning problem\n",
    "    \n",
    "    Math:\n",
    "        q_{ii} = s_{i}(s_{i}-c) \\\\\n",
    "        q_{ij} = q_{ji} = s_{i}s_{j} \\\\\n",
    "        c = \\sum_i{s_{i}}\n",
    "    \"\"\"\n",
    "    number_set = np.array(list(number_set))\n",
    "    c = number_set.sum()\n",
    "    q = np.outer(number_set, number_set)\n",
    "    np.fill_diagonal(q, number_set*(number_set-c))\n",
    "    return q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qubo_number_partition():\n",
    "    s = np.array([3,1,1,2,2,1])\n",
    "    sample, _ = qbsolv_from_ndarray(qubo_number_partition(s))\n",
    "    assert s[sample[0] == 1].sum() == 5\n",
    "    assert s[sample[0] == 0].sum() == 5\n",
    "test_qubo_number_partition()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Max-Cut Problem\n",
    "\n",
    "- https://en.wikipedia.org/wiki/Maximum_cut\n",
    "- 先程と同じく，各ノードに対してqubitを割り当て，2つの集合に分けることを考える\n",
    "- エッジ(i,j)に対して，$(x_{i}-x_{j})^2=1$のとき，$(i,j)$がカットに含まれることを利用する\n",
    "- 後々グラフを扱う問題がいくつか出てくるので，Graph Classを用意する\n",
    "    - データ形式のチェックなどはこれに任せる．（とりあえずは最低限のチェックのみ実装する）\n",
    "    - networkxとのやり取り機能をつけておく"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "@dataclass\n",
    "class Graph:\n",
    "    edges: Collection[Collection[int]] # edge is represented as (i, j)\n",
    "    n_nodes: int\n",
    "    init_value: Mapping[int, int] = None # init_value[i] indicates init value of node i\n",
    "    nodes: Sequence[float] = None # value of nodes. len(nodes) must be n_nodes\n",
    "     \n",
    "    def __post_init__(self):\n",
    "        if self.nodes is not None: assert len(self.nodes) == self.n_nodes \n",
    "        if self.edges is None: return\n",
    "        if len(self.edges) == 0: \n",
    "            self.edges = None\n",
    "            return\n",
    "        self.edges = np.array(list(self.edges))\n",
    "        assert self.edges.max() < self.n_nodes\n",
    "        assert self.edges.shape[1] == 2\n",
    "        if self.init_value is not None:\n",
    "            assert all(np.array(list(self.init_value)) < self.n_nodes)\n",
    "        \n",
    "    @classmethod\n",
    "    def from_networkx(cls, graph: nx.Graph) -> Graph:\n",
    "        return cls(edges = graph.edges, n_nodes = graph.number_of_nodes())\n",
    "    \n",
    "    def to_networkx_graph(self) -> nx.Graph:\n",
    "        nxg = nx.Graph()\n",
    "        nxg.add_nodes_from(range(self.n_nodes))\n",
    "        nxg.add_edges_from(self.edges)\n",
    "        return nxg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_graph():\n",
    "    edges = {(1,2), (3,4)}\n",
    "    graph = Graph(edges, 5)\n",
    "    assert graph.edges.shape == (2,2)\n",
    "    nxg = graph.to_networkx_graph()\n",
    "    _graph = Graph.from_networkx(nxg)\n",
    "    assert np.array_equal(graph.edges, _graph.edges)\n",
    "    assert graph.n_nodes == _graph.n_nodes\n",
    "test_graph()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- max-cut 実装"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubo_max_cut(g: Graph) -> np.ndarray:\n",
    "    n_nodes = g.n_nodes\n",
    "    q = np.zeros((n_nodes, n_nodes))\n",
    "    i, j = g.edges.T\n",
    "    np.add.at(q, (i,i), -1)\n",
    "    np.add.at(q, (j,j), -1)\n",
    "    q[i,j] += 1\n",
    "    q[j,i] += 1\n",
    "    return q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- test\n",
    "    - test用のhelper functionを用意しておく"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_contain(x: np.ndarray, y: np.ndarray) -> bool:\n",
    "    \"\"\"If y \\in x\"\"\"\n",
    "    assert len(x.shape) -1 == len(y.shape)\n",
    "    return any(np.array_equal(xi, y) for xi in x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qubo_max_cut():\n",
    "    edges = [(0,1), (0,2), (1,3), (2,3), (2, 4), (3, 4)]\n",
    "    g = Graph(edges=edges, n_nodes=5)\n",
    "    q = qubo_max_cut(g)\n",
    "    assert q.tolist() == [[-2,1,1,0,0],\n",
    "                          [1,-2,0,1,0],\n",
    "                          [1,0,-3,1,1],\n",
    "                          [0,1,1,-3,1],\n",
    "                          [0,0,1,1,-2]]\n",
    "    sample, _ = qbsolv_from_ndarray(q)\n",
    "    assert is_contain(sample, np.array([0,1,1,0,0]))\n",
    "test_qubo_max_cut()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# known penalties\n",
    "\n",
    "- 上の2つのモデルは，制約式がなかったため非常に簡単にQUBOモデルにできた\n",
    "- しかし制約がある場合はどうしたらよいか？一般的に下の表が役に立つ （[A Tutorial on Formulating and Using QUBO Models](http://meta-analytics.net/references/QUBO%20Tutorial%20%20Version%201-4.pdf)\n",
    "![](img/known_penalties.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Minimum Vertex Cover Problem (MVC)\n",
    "\n",
    "- [Vertex cover \\- Wikipedia](https://en.wikipedia.org/wiki/Vertex_cover)\n",
    "- グラフのすべての枝について，端点の少なくとも一方が属するような頂点の部分集合を求める問題"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubo_mvc(g: Graph, penalty: float = 10.):\n",
    "    \"\"\"The Minimum Vertex Cover Problem (MVC)\n",
    "    \n",
    "    Math:\n",
    "        y = \\sum_{j \\in V}x_{j} + P\\left(\\sum_{(i,j) \\in E}\\left(1-x_{i}-x_{j}+x_{i}x_{j}\\right)\\right)\n",
    "    \"\"\"\n",
    "    q = np.diagflat(np.ones(g.n_nodes))\n",
    "    i, j = g.edges.T\n",
    "    np.add.at(q, (i,i), -penalty)\n",
    "    np.add.at(q, (j,j), -penalty)\n",
    "    q[i,j] += penalty/2.\n",
    "    q[j,i] += penalty/2.\n",
    "    return q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qubo_mvc():\n",
    "    edges = [(0,1), (0,2), (1,3), (2,3), (2, 4), (3, 4)]\n",
    "    g = Graph(edges=edges, n_nodes=5)\n",
    "    q = qubo_mvc(g, penalty=8.)\n",
    "    sample,_ = qbsolv_from_ndarray(q)\n",
    "    assert is_contain(sample, np.array([0,1,1,0,1]))\n",
    "test_qubo_mvc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Weighted Minimum Vertex Cover Problem (W-MVC)\n",
    "\n",
    "- MVCの重み付きバージョン"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubo_wmvc(g: Graph, penalty: float = 10.):\n",
    "    \"\"\"The Minimum Vertex Cover Problem (MVC)\n",
    "    \n",
    "    Math:\n",
    "        y = \\sum_{j \\in V}w_{j}x_{j} + P\\left(\\sum_{(i,j) \\in E}\\left(1-x_{i}-x_{j}+x_{i}x_{j}\\right)\\right)\n",
    "    \"\"\"\n",
    "    q = np.diagflat(np.ones(g.n_nodes)*g.nodes)\n",
    "    i, j = g.edges.T\n",
    "    np.add.at(q, (i,i), - penalty*g.nodes[i])\n",
    "    np.add.at(q, (j,j), - penalty*g.nodes[j])\n",
    "    q[i,j] += penalty/2.\n",
    "    q[j,i] += penalty/2.\n",
    "    return q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- test\n",
    "    - めんどいので重みは全部1にする"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qubo_wmvc():\n",
    "    edges = [(0,1), (0,2), (1,3), (2,3), (2, 4), (3, 4)]\n",
    "    g = Graph(edges=edges, n_nodes=5, nodes=np.ones(5))\n",
    "    q = qubo_wmvc(g, penalty=8.)\n",
    "    sample, _ = qbsolv_from_ndarray(q)\n",
    "    assert is_contain(sample, np.array([0,1,1,0,1]))\n",
    "test_qubo_wmvc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remarks about the scalar penalty P\n",
    "\n",
    "> Generally, there is a ‘Goldilocks region’ of considerable size that contains penalty values that work well. A little preliminary thought about the model can yield a ballpark estimate of the original objective function value. Taking P to be some percentage (75% to 150%) of this estimate is often a good place to start. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Set Packing Problem\n",
    "- [Set packing \\- Wikipedia](https://en.wikipedia.org/wiki/Set_packing)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubo_set_pack(a: np.ndarray, weight: np.ndarray, penalty=8.) -> np.ndarray:\n",
    "    \"\"\"Set Packing Problem\n",
    "    \n",
    "    Math:\n",
    "        min\\ \\sum_{j=1}^n w_{j}x_{j} \\\\\n",
    "        st \\\\\n",
    "        \\sum_{j=1}^n a_{ij}x_{j} \\le 1\n",
    "    \"\"\"\n",
    "    assert len(a.shape) == 2\n",
    "    assert a.shape[-1] == len(weight)\n",
    "    q = -np.diagflat(weight)\n",
    "    c = np.einsum(\"ij,ik->ijk\",a,a).astype(np.float) / 2. # constraint\n",
    "    c *= (1 - np.eye(a.shape[-1]))[None,...] # zeros diag\n",
    "    q += penalty*c.sum(0)\n",
    "    return q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qubo_set_pack():\n",
    "    a = np.array([[1,0,1,1],\n",
    "                  [1,1,0,0]])\n",
    "    w = np.ones(4)\n",
    "    q = qubo_set_pack(a, w, penalty=6.)\n",
    "    assert q.tolist() == [[-1,3,3,3],\n",
    "                         [3,-1,0,0],\n",
    "                         [3,0,-1,3],\n",
    "                         [3,0,3,-1]]\n",
    "    assert is_contain(qbsolv_from_ndarray(q)[0], np.array([0,1,1,0]))\n",
    "test_qubo_set_pack()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Max 2-Sat Problem\n",
    "- [Boolean satisfiability problem \\- Wikipedia](https://en.wikipedia.org/wiki/Boolean_satisfiability_problem)\n",
    "- 面白いことに，Clauseの数によらないモデルとなる\n",
    "- Clauseのデータ型を定義しておく\n",
    "    - リテラルが否定なら sign = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "@dataclass\n",
    "class Clauses:\n",
    "    literals: np.ndarray[int]\n",
    "    signs: np.ndarray[bool]\n",
    "    def __post_init__(self):\n",
    "        assert len(self.literals) == len(self.signs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 実装"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubo_max2sat(c: Clauses):\n",
    "    n = c.literals.max()+1\n",
    "    q = np.zeros((n,n))\n",
    "    i,j=c.literals.T\n",
    "    si,sj=c.signs.T\n",
    "    np.add.at(q, (i[sj], i[sj]), ((-1)**si)[sj])\n",
    "    np.add.at(q, (j[si], j[si]), ((-1)**sj)[si])\n",
    "\n",
    "    offdiag=np.zeros_like(q)\n",
    "    np.add.at(offdiag, (i,j), (-1)**(si^sj) / 2.)\n",
    "    return offdiag + offdiag.T + q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qubo_max2sat():\n",
    "    l = np.array([[0,1],[0,1],[0,1],[0,1],[0,2],\n",
    "                  [0,2],[1,2],[1,3],[1,2],[1,2],\n",
    "                  [2,3],[2,3]])\n",
    "    s = np.array([[True,True],[True,False],[False,True],[False,False],\n",
    "                  [False,True],[False,False],[True,False],[True,True],\n",
    "                  [False,True],[False,False],[True,True],[False,False]])\n",
    "    q = qubo_max2sat(Clauses(l,s))\n",
    "    assert q.tolist() == [[1,0,0,0],\n",
    "                          [0,0,-0.5,0.5],\n",
    "                          [0,-0.5,0,1],\n",
    "                          [0,0.5,1,-2]]\n",
    "    sample, _ = qbsolv_from_ndarray(q)\n",
    "    assert is_contain(sample, np.array([0,0,0,1]))\n",
    "test_qubo_max2sat()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Set Partitioning Problem (SPP)\n",
    "- [Partition of a set \\- Wikipedia](https://en.wikipedia.org/wiki/Partition_of_a_set)\n",
    "- この問題ではqubitは集合の要素ではなく，部分集合を表す\n",
    "- SPPはクラスタリングともみなすことができ，機械学習との関連が資料の中で紹介されている"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubo_spp(cost: np.ndarray, set_flag: np.ndarray, p=10.):\n",
    "    \"\"\"The Set Partitioning Problem (SPP)\n",
    "    \n",
    "    Args:\n",
    "        set_flag: 2-rank array. (i,j) = 1 if element i in set j.\n",
    "    \"\"\"\n",
    "    assert cost.shape[0] == set_flag.shape[1]\n",
    "    b = np.ones(set_flag.shape[0])\n",
    "    return np.diagflat(cost - p*2*b.dot(set_flag))  + p*np.einsum(\"ij,ik->jk\",set_flag,set_flag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qubo_spp():\n",
    "    c = np.array([3,2,1,1,3,2])\n",
    "    a = np.array([[1,0,1,0,0,1],\n",
    "                  [0,1,1,0,1,1],\n",
    "                  [0,0,1,1,1,0],\n",
    "                  [1,1,0,1,0,1]])\n",
    "\n",
    "    q=qubo_spp(c, a)\n",
    "    assert q.tolist() == [[-17,10,10,10,0,20],\n",
    "                        [10,-18,10,10,10,20],\n",
    "                        [10,10,-29,10,20,20],\n",
    "                        [10,10,10,-19,10,10],\n",
    "                        [0,10,20,10,-17,10],\n",
    "                        [20,20,20,10,10,-28]]\n",
    "    sample, _ = qbsolv_from_ndarray(q)\n",
    "    assert is_contain(sample, np.array([1,0,0,0,1,0]))\n",
    "test_qubo_spp()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Graph Coloring\n",
    "- [Graph coloring \\- Wikipedia](https://en.wikipedia.org/wiki/Graph_coloring)\n",
    "- 応用範囲が広い\n",
    "    - スケジューリング，基板デザイン，レジスタ割付\n",
    "    - 以前[数独の実装](https://qiita.com/tamurahey/items/105f1dc9ee9a3bc01f15)をした\n",
    "- ノードiが色jで塗られるとき，qubitを1とする\n",
    "    - つまり見通しを良くするために，qubitは行列で，qubo行列は4階のndarrayで表現する"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubo_graph_coloring(g: Graph, n_color, penalty=10.) -> np.ndarray:\n",
    "    \"\"\"Graph Coloring\"\"\"\n",
    "    n = g.n_nodes\n",
    "    q=np.zeros((n,n_color,n,n_color))\n",
    "    \n",
    "    # each nodes are colored by one color\n",
    "    i=range(n)\n",
    "    q[i, :, i, :] = penalty*(np.ones((n_color,n_color)) - 2*np.eye(n_color))\n",
    "    \n",
    "    # adjacent nodes are not colored with same color\n",
    "    i, j = g.edges.T[...,None]\n",
    "    k=range(n_color)\n",
    "    q[i,k,j,k]+=penalty/2.\n",
    "    q[j,k,i,k]+=penalty/2.\n",
    "    \n",
    "    return q.reshape(n*n_color,n*n_color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qubo_graph_coloring():\n",
    "    edges = [[0,1],[0,4],[1,2],[1,3],[1,4],[2,3],[3,4]]\n",
    "    n_nodes, n_color=5, 3\n",
    "    g = Graph(edges=edges, n_nodes=n_nodes)\n",
    "    ans=np.zeros(n_nodes*n_color)\n",
    "    ans[[1,3,8,10,14]]=1\n",
    "    q=qubo_graph_coloring(g, n_color, penalty=4.)\n",
    "    assert is_contain(qbsolv_from_ndarray(q)[0],ans)\n",
    "test_qubo_graph_coloring()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> This approach to coloring problems has proven to be very effective for a wide variety of coloring instances with hundreds of nodes.\n",
    "\n",
    "らしい.n_color色に対してn_color qubitを用意しているが，log_2(n_color)ですむような構成法はないか？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# General 0/1 Programming\n",
    "\n",
    "> Many important problems in industry and government can be modeled as 0/1 linear programs with a mixture of constraint types.\n",
    "\n",
    "- 不等式制約をslack変数の導入によってQUBOに入れている点が面白い"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubo_general01(cost, a, b, sign, penalty=10.):\n",
    "    \"\"\"General 0/1 Programming\n",
    "    \n",
    "    Args:\n",
    "        a, b: constraint (ax=b)\n",
    "        sign: 0 is equality, 1 is >=, -1 is <=\n",
    "    \"\"\"\n",
    "    slack = - sign * b\n",
    "    slack[sign==0] = 1.\n",
    "    slack[sign == 1] += a[sign==1,:].sum()\n",
    "    slack = np.ceil(np.log2(slack)).astype(np.int)\n",
    "\n",
    "    ni = a.shape[0]\n",
    "    for i, ns in enumerate(slack):\n",
    "        if ns:\n",
    "            t = np.zeros((ni,ns))\n",
    "            t[i] = -sign[i]*2**np.arange(ns)\n",
    "            a=np.concatenate([a,t],-1)\n",
    "            \n",
    "    return - np.diagflat(np.concatenate([cost, np.zeros(sum(slack))])+2*penalty*b.dot(a)) + penalty*np.einsum(\"ij,ik->jk\",a,a) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qubo_general01():\n",
    "    c = np.array([6,4,8,5,5])\n",
    "    a = np.array([[2,2,4,3,2],\n",
    "                  [1,2,2,1,2],\n",
    "                  [3,3,2,4,4]])\n",
    "    b = np.array([7,4,5])\n",
    "    s = np.array([-1,0,1])\n",
    "    q=qubo_general01(c,a,b,s)\n",
    "    sample, _=qbsolv_from_ndarray(q)\n",
    "    assert is_contain(sample[:,:5], np.array([1,0,0,1,1]))\n",
    "test_qubo_general01()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quadratic Assignment Problem (QAP)\n",
    "\n",
    "- [Quadratic assignment problem \\- Wikipedia](https://en.wikipedia.org/wiki/Quadratic_assignment_problem)k\n",
    "- 工場の配置など"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubo_qap(flow:np.ndarray, distance: np.ndarray, penalty=10.):\n",
    "    \"\"\"Quadratic Assignment Problem (QAP)\"\"\"\n",
    "    n = len(flow)\n",
    "    q = np.einsum(\"ij,kl->ikjl\",flow, distance).astype(np.float)\n",
    "\n",
    "    i = range(len(q))\n",
    "    q[i,:,i,:] += penalty\n",
    "    q[:,i,:,i] += penalty\n",
    "    q[i,i,i,i] -= 4*penalty\n",
    "    return q.reshape(n**2,n**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qubo_qap():\n",
    "    flow = [[0,5,2],\n",
    "            [5,0,3],\n",
    "            [2,3,0]]\n",
    "    distance = [[0,8,15],\n",
    "                [8,0,13],\n",
    "                [15,13,0]]\n",
    "    q = qubo_qap(flow, distance, penalty=200.)\n",
    "    sample, _ = qbsolv_from_ndarray(q)\n",
    "    assert is_contain(sample, np.array([1,0,0,0,1,0,0,0,1]))\n",
    "test_qubo_qap()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quadratic Knapsack Problem (QKP)\n",
    "\n",
    "- 2次のナップサック問題\n",
    "\n",
    "> widespread application in such areas as project selection and capital budgeting.\n",
    "\n",
    "> In such settings, a set of attractive potential projects is identified and the goal is to identify a subset of maximum value (or profit) that satisfies the budget limitations.\n",
    "> The quadratic version of this problem arises when there is an interaction between pairs of projects affecting the value obtained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubo_qkp(value:np.ndarray, a:np.ndarray, b:float, penalty=10.) -> np.ndarray:\n",
    "    \"\"\"Quadratic Knapsack Problem (QKP)\"\"\"\n",
    "    n = len(value)\n",
    "    nslack = np.ceil(np.log2(b))\n",
    "    slack = 2 ** (np.arange(nslack))\n",
    "    a = np.concatenate([a,slack])\n",
    "    q = penalty*(np.outer(a,a) - 2*b*np.diag(a))\n",
    "    q[:n,:n] -= value\n",
    "    return q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_qubo_qkp():\n",
    "    value = np.array([[2,4,3,5],\n",
    "                      [4,5,1,3],\n",
    "                      [3,1,2,2],\n",
    "                      [5,3,2,4]])\n",
    "    a = np.array([8,6,5,3])\n",
    "    b = 16\n",
    "    q = qubo_qkp(value, a, b, penalty=10.)\n",
    "    sample, _ = qbsolv_from_ndarray(q)\n",
    "    n = len(a)\n",
    "    assert is_contain(sample[:,:n], np.array([1,0,1,1]))\n",
    "test_qubo_qkp()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The QUBO approach to QKP has proven to be successful on problems with several hundred variables and as many as five knapsack constraints."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "qubogen",
   "language": "python",
   "name": "qubogen"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}