{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "sys.path.insert(1, \"../../\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 二次割当問題\n",
    "\n",
    "二次割当問題は、以下のような問題です。\n",
    "\n",
    "```{card} 二次割当問題\n",
    "$N$ を正の整数とする。$N$ 個の工場建設候補地に $N$ 個の工場を建設することを考える。それぞれの工場は、どの建設候補地に建設されても良い。どの 2 つの工場にも、それらを行き来するトラックが走っており、その輸送量はあらかじめ分かっている。輸送量 x 走行距離の和を最小化するにはどうすればよいか。\n",
    "```\n",
    "\n",
    "応用として、集会の座席表を仲が良い人同士が近くなるように決定することなどが考えられます。\n",
    "\n",
    "## 定式化\n",
    "\n",
    "$N$ 個の工場建設候補地を 土地 $0$, 土地 $1$, ..., 土地 $N-1$ と表記し、$N$ 個の工場を工場 $0$, 工場 $1$, ..., 工場 $N-1$ と表記します。また土地 $i$ と土地 $j$ の距離を $D_{i, j}$, 工場 $k$ と工場 $l$ の間の輸送量を $F_{k, l}$ とします。\n",
    "\n",
    "### 変数\n",
    "\n",
    "$N \\times N$ 個のバイナリ変数 $q$ を用意し、$q_{i, k}$ は工場 $k$ を土地 $i$ に建設するかどうかを表すことにします。\n",
    "\n",
    "たとえば、$q$ が以下のような値をとるとき、土地 $0$ には 工場 $3$ が建設されます。\n",
    "\n",
    "```{csv-table} バイナリ変数テーブル\n",
    ":header-rows: 1\n",
    ":stub-columns: 1\n",
    "\n",
    "\"\", \"工場 0\", \"工場 1\", \"工場 2\", \"工場 3\", \"工場 4\"\n",
    "\"土地 0\", 0, 0, 0, 1, 0\n",
    "\"土地 1\", 0, 1, 0, 0, 0\n",
    "\"土地 2\", 0, 0, 0, 0, 1\n",
    "\"土地 3\", 1, 0, 0, 0, 0\n",
    "\"土地 4\", 0, 0, 1, 0, 0\n",
    "```\n",
    "\n",
    "### 制約条件\n",
    "\n",
    "バイナリ変数テーブルのそれぞれの行と列には 1 となる値がちょうど 1 個である必要があります。したがって、各行各列に one-hot 制約をかけます。逆に、これらがみたされていれば、どの土地にどの工場を建設するかが 1 通りに定まります。\n",
    "\n",
    "### 目的関数\n",
    "\n",
    "目的関数は、工場と工場の間の輸送量 x 距離の総和です。これは $q$ を用いて数式で表すと\n",
    "\n",
    "$$\n",
    "\\sum_{q_{i, k} = 1, q_{j, l} = 1} D_{i, j} \\ F_{k, l} = \\sum_{i, j, k, l} q_{i, k} \\ q_{j, l} \\ D_{i, j} \\ F_{k, l}\n",
    "$$\n",
    "\n",
    "となります。\n",
    "\n",
    "### 数式\n",
    "\n",
    "以上の定式化は、$N\\times N$ 個のバイナリ変数 $q$ を用いて\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\text{minimize} \\quad   &\\sum_{i, j, k, l} q_{i, k} \\ q_{j, l} \\ D_{i, j} \\ F_{k, l} \\\\\n",
    "\\text{subject to} \\quad &\\sum_k q_{i, k} = 1 \\quad \\text{for} \\quad i \\in \\{0, 1, \\ldots, N - 1\\}, \\\\ \n",
    "                        &\\sum_i q_{i, k} = 1 \\quad \\text{for} \\quad k \\in \\{0, 1, \\ldots, N - 1\\}, \\\\\n",
    "                        &q_{i, k} \\in \\{0, 1\\} \\quad \\text{for} \\quad i, k \\in \\{0, 1, \\ldots, N - 1\\}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "と書くことができます。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 問題の作成\n",
    "\n",
    "Amplify SDK による定式化を行う前に、問題を作成しておきます。簡単のため工場の数 $N$ を 10 とします。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=3, linewidth=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "土地間の距離を表す行列 $D$ を作成します。土地はユークリッド平面上にランダムに生成します。距離行列は 2 次元の {py:class}`numpy.ndarray` として作成します。名前は `distance` とします。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng()\n",
    "\n",
    "x = rng.integers(0, 100, size=(N,))\n",
    "y = rng.integers(0, 100, size=(N,))\n",
    "\n",
    "distance = (\n",
    "    (x[:, np.newaxis] - x[np.newaxis, :]) ** 2\n",
    "    + (y[:, np.newaxis] - y[np.newaxis, :]) ** 2\n",
    ") ** 0.5\n",
    "\n",
    "print(distance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "工場間の輸送量を表す行列 $F$ を作成します。2 次元の対称行列をランダムに作成し、名前は `flow` とします。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flow = np.zeros((N, N), dtype=int)\n",
    "for i in range(N):\n",
    "    for j in range(i + 1, N):\n",
    "        flow[i, j] = flow[j, i] = rng.integers(0, 100)\n",
    "\n",
    "print(flow)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Amplify SDK による定式化\n",
    "\n",
    "Amplify SDK による定式化を行います。定式化において、どの 2 つのバイナリ変数からなる 2 次項も目的関数に現れるので、{py:class}`~amplify.Matrix` クラスを使うと効率的な定式化を行うことができます。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 変数の作成\n",
    "\n",
    "{py:class}`~amplify.Matrix` クラスを使用した定式化を行うためには、{py:class}`~amplify.VariableGenerator` の {py:meth}`~amplify.VariableGenerator.matrix` メソッドを使用して変数を発行します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amplify import VariableGenerator\n",
    "\n",
    "gen = VariableGenerator()\n",
    "matrix = gen.matrix(\"Binary\", N, N)  # coefficient matrix\n",
    "q = matrix.variable_array  # variables\n",
    "\n",
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 目的関数の作成\n",
    "\n",
    "上で作成した `matrix` は {py:class}`~amplify.Matrix` クラスのインスタンスであり、\n",
    "\n",
    "* {py:attr}`~amplify.Matrix.quadratic`\n",
    "* {py:attr}`~amplify.Matrix.linear`\n",
    "* {py:attr}`~amplify.Matrix.constant`\n",
    "\n",
    "の 3 種類のプロパティを持ちます。\n",
    "\n",
    "{py:attr}`~amplify.Matrix.quadratic` は 2 次の項の係数を表す {py:class}`numpy.ndarray` で、その {py:attr}`~numpy.ndarray.shape` は今回は `(N, N, N, N)` となっています。`quadratic[i, k, j, l]`{l=python} は `q[i, k] * q[j, l]`{l=python} の係数に対応します。つまり {py:attr}`~amplify.Matrix.quadratic` には `quadratic[i, k, j, l] = distance[i, j] * flow[k, l]`{l=python} となるように 4 次元 NumPy 配列を設定する必要があります。\n",
    "\n",
    "{py:attr}`~amplify.Matrix.linear` と {py:attr}`~amplify.Matrix.constant` はそれぞれ線形項の係数と定数項を表しますが、今回使用する目的関数には 2 次の項しか含まれないため、設定しません。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "np.einsum(\"ij,kl->ikjl\", distance, flow, out=matrix.quadratic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 制約条件の作成\n",
    "\n",
    "[](#変数の作成) で作成した変数配列 `q` の各行各列に one-hot 制約をかけます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amplify import one_hot\n",
    "\n",
    "constraints = one_hot(q, axis=1) + one_hot(q, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(qap-model)=\n",
    "\n",
    "### 組合せ最適化モデルの作成\n",
    "\n",
    "目的関数と制約条件を組み合わせて、モデルを作成します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "penalty_weight = np.max(distance) * np.max(flow) * (N - 1)\n",
    "model = matrix + penalty_weight * constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "制約条件に `penalty_weight` をかけているのは、制約条件に重みをつけるためです。今回使用するソルバーである Amplify AE においては、制約条件に適切な重みを指定しないと、ソルバーが制約条件をみたそうとするよりも目的関数を小さくする方向に動いてしまい、良い解を発見することができなくなってしまいます。詳しくは [](penalty.md) を参照してください。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ソルバークライアントの作成\n",
    "\n",
    "Amplify AE を用いて組合せ最適化を実行するために、ソルバークライアントを作成します。Amplify AE に対応するソルバークライアントクラスは {py:class}`~amplify.FixstarsClient` クラスです。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amplify import FixstarsClient\n",
    "\n",
    "client = FixstarsClient()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Amplify AE の実行に必要な API トークンを設定します。\n",
    "\n",
    "```{tip}\n",
    "[ユーザ登録](https://amplify.fixstars.com/register)を行うと、評価・検証目的に使える API トークンを無料で入手できます。\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "skip-execution"
    ]
   },
   "outputs": [],
   "source": [
    "client.token = \"xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "client.token = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ソルバーのタイムアウト値を設定します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "\n",
    "client.parameters.timeout = datetime.timedelta(seconds=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 求解の実行\n",
    "\n",
    "作成した組合せ最適化モデルとソルバークライアントを使用してソルバーの実行を行い、二次計画問題の解を求めます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amplify import solve\n",
    "\n",
    "result = solve(model, client)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "目的関数の値は以下のように表示できます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result.best.objective"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最適解における変数の値は、NumPy の多次元配列の形式で以下のようにして取得できます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_values = q.evaluate(result.best.values)\n",
    "\n",
    "print(q_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 結果の確認\n",
    "\n",
    "matplotlib を用いて結果を可視化します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x, y)\n",
    "factory_indices = (q_values @ np.arange(N)).astype(int)\n",
    "\n",
    "for i, j in itertools.combinations(range(N), 2):\n",
    "    plt.plot(\n",
    "        [x[i], x[j]],\n",
    "        [y[i], y[j]],\n",
    "        c=\"b\",\n",
    "        alpha=flow[factory_indices[i], factory_indices[j]] / 100,\n",
    "    )"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "amplify-lrEeNlW4",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}