notebooks/RL_Markov_Decision_Process.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MDPについて\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MDPの生成について\n",
    "\n",
    "参考文献\n",
    "* [Average-reward model-free reinforcement learning: a systematic review and literature mapping](https://arxiv.org/abs/2010.08920)\n",
    "* [On the Generation of Markov Decision Processes](https://www.jstor.org/stable/2584329)\n",
    "* [Natural Actor Critic Algorithms](https://hal.inria.fr/hal-00840470/document)\n",
    "\n",
    "強化学習アルゴリズムの性能を比較したいとき、どんなMDPを使って評価するべきでしょうか？\n",
    "適当にMDPを生成しても正しい評価にはなりません。\n",
    "適当に生成したMDPが特定のアルゴリズムが得意なMDPになっていまう可能性は十分にあります。\n",
    "\n",
    "### Garnet\n",
    "\n",
    "MDPの適切な生成方法についての議論では[On the Generation of Markov Decision Processes](https://www.jstor.org/stable/2584329)がよく引用されます。\n",
    "評価用のMDPを作る際は以下の要素を考慮するべきだとこの論文では議論されています：\n",
    "1. 割引率 $\\gamma$: どんなアルゴリズムでも割引率が小さいほど良い収束性を示します。\n",
    "2. 最適性 $\\varepsilon$: どんなアルゴリズムでも$\\varepsilon$が小さいほど$\\varepsilon$-最適な方策を見つけるために必要なコストが大きくなります。\n",
    "3. 状態における選択可能な行動の数の平均と分散: \n",
    "4. 状態と行動における遷移可能な次状態の数: 一つの状態の情報を更新するコストは3と4によって決定されます。状態で選択可能な行動の数や次状態の数が多いほどコストが高いわけですね。\n",
    "5. 状態の数 $|S|$: 状態の数が多い & 一つの状態におけるコストが大きい とき、状態ごとに並列化して情報を処理するメリットが大きくなります。\n",
    "6. Mixingの早さ: どんな初期状態からスタートしても、別の$n$個の状態に行ってしまう確率が$n$にしたがって急激に高くなるとき、Mixingが早いと呼びます。早いMixingでは特定のアルゴリズム（Pre-Jacobiとか）の収束が早くなります。TODO [Dynamic Programming and Optimal Control](http://athenasc.com/dpbook.html)をよめ！\n",
    "7. 遷移行列が下三角行列に集中しているか？: 集中している場合の状態の並びをweighted lower triangular (WLT)と呼びます。WLTのとき、特定のアルゴリズム (Gauss-Seidelとか)の収束が早くなります。 TODO [Dynamic Programming and Optimal Control](http://athenasc.com/dpbook.html)をよめ！\n",
    "8. 非最適方策による遷移行列が最適方策のそれに近いか？: \n",
    "9. 非最適方策の価値が最適方策のそれに近いか？: 価値が似通っているとき、最適方策を見分けるのが難しくなります。つまり、最適方策に至るまでに必要なコストが高くなります。\n",
    "\n",
    "[Garnet](https://hal.inria.fr/hal-00840470/document) (Generic Average Reward Non-stationary Environment Testbed)はこの性質をいい感じにランダム化して生成できるようなMDPです（らしい。ほんまか？）。\n",
    "```\"Garnet problems do not correspond to any particular application, but are meant to be totally abstract or generic while remaining representative of the kind of MDPs that might be encountered in practice.\"```\n",
    "と良く言われています。\n",
    "\n",
    "Garnetの生成についてはいろいろパターンがありますが、共通して重要なパラメータが以下のBranching factorです：\n",
    "* $b$: それぞれの状態・行動について、遷移可能な次状態の数です。次状態は状態空間の中から重複なしでランダムに選択されます。遷移確率もランダムに選択されます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from functools import partial\n",
    "from typing import Optional, NamedTuple\n",
    "import jax.numpy as jnp\n",
    "\n",
    "def create_rew_P(S: int, A: int, b: int):\n",
    "    \"\"\" Garnet MDP (報酬関数, 遷移確率) を作成します\n",
    "\n",
    "    Args:\n",
    "        S (int): 状態空間のサイズ\n",
    "        A (int): 行動空間のサイズ\n",
    "        b (int): branching factor\n",
    "\n",
    "    Returns:\n",
    "        rew (float): 報酬関数\n",
    "        P (float): 遷移確率\n",
    "    \"\"\"\n",
    "    # ランダムな報酬関数を作ります\n",
    "    rew = np.random.rand(S, A)\n",
    "\n",
    "    # ランダムな遷移確率を作ります\n",
    "    P = np.random.rand(S * A, S)\n",
    "\n",
    "    for idx in range(S*A):\n",
    "        unreachable_states = np.random.choice(S, size=S-b, replace=False)\n",
    "        P[idx, unreachable_states] = 0\n",
    "    P = P / np.sum(P, axis=-1, keepdims=True)  # 正規化して確率にします\n",
    "    P = P.reshape(S, A, S)\n",
    "    np.testing.assert_almost_equal(P.sum(axis=-1), 1)\n",
    "    assert np.all(np.count_nonzero(P, axis=-1) == b)\n",
    "    return jnp.array(rew), jnp.array(P)\n",
    "\n",
    "\n",
    "class MDP(NamedTuple):\n",
    "    S_array: np.array  # 状態空間\n",
    "    A_array: np.array  # 行動空間\n",
    "    gamma: float  # 割引率\n",
    "    horizon: int  # ホライゾン\n",
    "    rew: float  # 報酬関数\n",
    "    P: float  # 遷移確率\n",
    "    optimal_Q: Optional[np.ndarray] = None  # 最適Q値\n",
    "    optimal_Q_idx: Optional[np.ndarray] = None  # 後でplot用に使います\n",
    "\n",
    "    @property\n",
    "    def S(self) -> int:  # 状態空間のサイズ\n",
    "        return len(self.S_array)\n",
    "\n",
    "    @property\n",
    "    def A(self) -> int:  # 行動空間のサイズ\n",
    "        return len(self.A_array)\n",
    "\n",
    "\n",
    "S = 50\n",
    "A = 2\n",
    "b = 2\n",
    "S_array = np.arange(S)\n",
    "A_array = np.arange(A)\n",
    "gamma = 0.99\n",
    "horizon = int(1 / (1 - gamma))\n",
    "\n",
    "# Garnet MDPを作ります\n",
    "t_rew, t_P = create_rew_P(S, A, b)\n",
    "garnet_MDP = MDP(S_array, A_array, gamma, horizon, t_rew, t_P)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GarnetMDPでValue Iterationを動かしてみましょう。\n",
    "\n",
    "実験の前に、いくつか便利な関数を導入します。\n",
    "\n",
    "1. ``compute_greedy_policy``: Q関数 ($S \\times A \\to \\mathcal{R}$) の貪欲方策を返します\n",
    "2. ``compute_optimal_Q``: MDPの最適Q関数 $q_* : S \\times A \\to \\mathcal{R}$ を返します。\n",
    "3. ``compute_policy_Q``: 方策 $\\pi$ のQ関数 $q_\\pi : S \\times A \\to \\mathcal{R}$ を返します。\n",
    "4. ``compute_optimality_gap``: 方策 $\\pi$ と最適解との差 $\\|q_* - q_{\\pi}\\|_\\infty$ を返します。この値が小さいほど、$\\pi$ はほぼ全ての状態行動対で最適な方策になっています。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "import jax\n",
    "\n",
    "\n",
    "# >>>> 1 >>>>\n",
    "@jax.jit\n",
    "def compute_greedy_policy(Q: np.ndarray):\n",
    "    greedy_policy = jnp.zeros_like(Q)\n",
    "    S, A = Q.shape\n",
    "    greedy_policy = greedy_policy.at[jnp.arange(S), Q.argmax(axis=1)].set(1)\n",
    "    assert greedy_policy.shape == (S, A)\n",
    "    return greedy_policy\n",
    "# <<<< 1 <<<<\n",
    "\n",
    "\n",
    "# >>>> 2 >>>>\n",
    "@partial(jax.jit, static_argnames=(\"S\", \"A\"))\n",
    "def _compute_optimal_Q(mdp: MDP, S: int, A: int):\n",
    "\n",
    "    def backup(optimal_Q):\n",
    "        greedy_policy = compute_greedy_policy(optimal_Q)\n",
    "        max_Q = (greedy_policy * optimal_Q).sum(axis=1)\n",
    "        next_v = mdp.P @ max_Q\n",
    "        assert next_v.shape == (S, A)\n",
    "        return mdp.rew + mdp.gamma * next_v\n",
    "    \n",
    "    optimal_Q = jnp.zeros((S, A))\n",
    "    body_fn = lambda i, Q: backup(Q)\n",
    "    return jax.lax.fori_loop(0, mdp.horizon + 1000, body_fn, optimal_Q)\n",
    "\n",
    "compute_optimal_Q = lambda mdp: _compute_optimal_Q(mdp, mdp.S, mdp.A)\n",
    "# <<<< 2 <<<<\n",
    "\n",
    "\n",
    "# >>>> 3 >>>>\n",
    "def compute_policy_Q(mdp: MDP, policy: np.ndarray):\n",
    "    S, A = policy.shape\n",
    "\n",
    "    def backup(policy_Q):\n",
    "        max_Q = (policy * policy_Q).sum(axis=1)\n",
    "        next_v = mdp.P @ max_Q\n",
    "        assert next_v.shape == (S, A)\n",
    "        return mdp.rew + mdp.gamma * next_v\n",
    "    \n",
    "    policy_Q = jnp.zeros((S, A))\n",
    "    body_fn = lambda i, Q: backup(Q)\n",
    "    return jax.lax.fori_loop(0, mdp.horizon + 1000, body_fn, policy_Q)\n",
    "# <<<< 3 <<<<\n",
    "\n",
    "# >>>> 4 >>>>\n",
    "@jax.jit\n",
    "def compute_optimality_gap(mdp: MDP, policy: np.ndarray):\n",
    "    policy_Q = compute_policy_Q(mdp, policy)\n",
    "    optimality_gap = jnp.abs(mdp.optimal_Q - policy_Q).max()\n",
    "    return optimality_gap\n",
    "# <<<< 4 <<<<\n",
    "\n",
    "garnet_MDP = garnet_MDP._replace(optimal_Q=compute_optimal_Q(garnet_MDP))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABMAUlEQVR4nO3dd3hUZd7/8c+Z9E4KSYAEEikJvcsCK+KqgApixYIKuM/KoyAg6gr7/DCCKwgqq4Jre3aBR7GsKIqs4iKKgqJUQURCR1oCAdIhZeb8/ggZCKFkwiQnmXm/risXc8qc+SZB/Xif871vwzRNUwAAAKj3bFYXAAAAAPcg2AEAAHgIgh0AAICHINgBAAB4CIIdAACAhyDYAQAAeAiCHQAAgIcg2AEAAHgIgh0AAICHINgBqHfmzp0rwzC0Z8+eOvHZffv2Vd++fWu9FgA4G8EOwCX75ZdfdM8996hJkyYKCAhQ48aNNXToUP3yyy+XdN2pU6fq448/dk+RtejgwYN66qmn9NNPP9XI9Tdt2qQRI0YoOTlZgYGBCg0NVadOnfTnP/9Zu3btqpHPBFA/GKwVC+BSfPTRR7rrrrsUFRWlP/7xj0pOTtaePXv0j3/8Q0ePHtV7772nm2++uVrXDg0N1W233aa5c+dW2G+321VSUqKAgAAZhuGG76Lq5s6dqxEjRmj37t1KSkqSJBUXF0uS/P39JUlr165V9+7dNWfOHA0fPtytn//mm2/qwQcfVExMjIYOHarU1FSVlpZq8+bN+vDDD3Xs2DGdOHFCPj4+bv1cAPWDr9UFAKi/du7cqXvvvVeXXXaZvv32WzVs2NB5bOzYsbriiit07733atOmTbrsssvc9rk+Pj51KriUB7qa9v333+vBBx9U7969tXjxYoWFhVU4/sILL+iZZ56plVoA1E3cigVQbc8995wKCwv1xhtvVAh1khQTE6PXX39dBQUFmjFjhnP/U089JcMwtHXrVg0ZMkTh4eGKjo7W2LFjdfLkSed5hmGooKBA8+bNk2EYMgzDOfp1rufckpKSNHDgQC1fvlzdunVTUFCQ2rdvr+XLl0sqG1ls3769AgMD1bVrV23YsKFCvZs2bdLw4cN12WWXKTAwUPHx8br//vt19OjRi/4cznzGbvny5erevbskacSIEc7a586dq7S0NPn5+enIkSOVrvHAAw+oQYMGFX4GZ5s8ebIMw9D8+fMrhTpJCgwM1NNPP10h9K5YsUK33367mjZtqoCAACUmJuqRRx7RiRMnKrx3+PDhCg0N1a5du9S/f3+FhISocePGmjJlirixA9QfBDsA1fbpp58qKSlJV1xxxTmP9+nTR0lJSfr3v/9d6diQIUN08uRJTZs2Tddff71efvllPfDAA87jb731lgICAnTFFVforbfe0ltvvaWRI0desJ4dO3bo7rvv1qBBgzRt2jQdP35cgwYN0vz58/XII4/onnvu0eTJk7Vz504NGTJEDofD+d6lS5dq165dGjFihGbNmqU777xT7733nq6//nqXgk3r1q01ZcoUSWVhrbz2Pn366N5771Vpaanef//9Cu8pLi7WggULdOuttyowMPCc1y0sLNRXX32lvn37KiEhocr1fPDBByosLNSDDz6oWbNmqX///po1a5buu+++Sufa7XYNGDBAcXFxmjFjhrp27aq0tDSlpaVV+fMAWMwEgGrIzs42JZmDBw++4Hk33nijKcnMzc01TdM009LSTEnmjTfeWOG8hx56yJRkbty40bkvJCTEHDZsWKVrzpkzx5Rk7t6927mvWbNmpiTz+++/d+774osvTElmUFCQuXfvXuf+119/3ZRkfv311859hYWFlT7n3XffNSWZ33777QU/+8orrzSvvPJK5/aaNWtMSeacOXMqXbNnz55mjx49Kuz76KOPKtVzto0bN5qSzHHjxlU6dvToUfPIkSPOr6Kiogt+X9OmTTMNw6jwMxk2bJgpyXz44Yed+xwOh3nDDTeY/v7+5pEjR85bG4C6gxE7ANWSl5cnSee8JXim8uO5ubkV9o8aNarC9sMPPyxJ+uyzz6pdU5s2bdSzZ0/ndo8ePSRJf/jDH9S0adNK+8/sIA0KCnK+PnnypLKysvS73/1OkrR+/fpq13S2++67Tz/++KN27tzp3Dd//nwlJibqyiuvPO/7yn9+oaGhlY5ddtllatiwofNr0aJFzmNnfl8FBQXKyspSr169ZJpmpdvRkjR69Gjna8MwNHr0aBUXF+vLL7907RsFYAmCHYBqKQ9s5QHvfM4XAFu2bFlhu3nz5rLZbJc0N92Z4U2SIiIiJEmJiYnn3H/8+HHnvmPHjmns2LGKi4tTUFCQGjZsqOTkZElSTk5OtWs62x133KGAgADNnz/fee3Fixdr6NChF+zwLf/55efnVzr2ySefaOnSpXr++ecrHfvtt980fPhwRUVFKTQ0VA0bNnQGyLO/L5vNVqnJpVWrVpJkyZyBAFxHVyyAaomIiFCjRo20adOmC563adMmNWnSROHh4Rc8zx3TlpyvU/Z8+80znp0bMmSIvv/+ez3++OPq1KmTQkND5XA4NGDAgArP4l2qyMhIDRw4UPPnz9eTTz6pBQsWqKioSPfcc88F39eiRQv5+vpq8+bNlY6VBzVf34r/Srfb7br22mt17NgxPfHEE0pNTVVISIgOHDig4cOHu/X7AlA3MGIHoNoGDhyo3bt3a+XKlec8vmLFCu3Zs0cDBw6sdGz79u0Vtnfs2CGHw+GcG05yT9iriuPHj2vZsmWaMGGCJk+erJtvvlnXXntttadouVjd9913n7Zt26Y1a9Zo/vz56ty5s9q2bXvB94SEhKhv37765ptvdODAgSrV8fPPP2vbtm164YUX9MQTT2jw4MG65ppr1Lhx43Oe73A4Kk1wvG3bNkmq8HsBUHcR7ABU2+OPP66goCCNHDmy0rQgx44d03//938rODhYjz/+eKX3vvLKKxW2Z82aJUm67rrrnPtCQkKUnZ3t/sLPUj6iZ57V/friiy9W63ohISGSdN7ar7vuOsXExGj69On65ptvLjpaV+7JJ5+U3W7XPffcc85bsmfXf67vyzRNvfTSS+f9jNmzZ1c4d/bs2fLz89PVV19dpRoBWItbsQCqrWXLlpo3b56GDh2q9u3bV1p5IisrS++++66aN29e6b27d+/WjTfeqAEDBmjVqlV6++23dffdd6tjx47Oc7p27aovv/xSM2fOVOPGjZWcnOxsfHCn8PBw9enTRzNmzFBJSYmaNGmi//znP9q9e3e1rte8eXM1aNBAr732msLCwhQSEqIePXo4n9nz8/PTnXfeqdmzZ8vHx0d33XVXla57xRVXaPbs2Xr44YfVsmVL58oTxcXF2rZtm+bPny9/f3/Fx8dLklJTU9W8eXM99thjOnDggMLDw/Xhhx9WeLbwTIGBgVqyZImGDRumHj166PPPP9e///1v/eUvf6k0TyGAOsrCjlwAHmLTpk3mXXfdZTZq1Mj08/Mz4+Pjzbvuusv8+eefK51bPt3Jli1bzNtuu80MCwszIyMjzdGjR5snTpyocO7WrVvNPn36mEFBQaYk59Qn55vu5IYbbqj0eZLMUaNGVdi3e/duU5L53HPPOfft37/fvPnmm80GDRqYERER5u23324ePHjQlGSmpaU5z6vKdCemaZqffPKJ2aZNG9PX1/ecU5+sXr3alGT269fvHD/RC9uwYYN53333mU2bNjX9/f3NkJAQs0OHDuajjz5q7tixo8K5W7ZsMa+55hozNDTUjImJMf/0pz85p045s6Zhw4aZISEh5s6dO81+/fqZwcHBZlxcnJmWlmba7XaXawRgDdaKBVCrnnrqKU2ePFlHjhxRTEyM1eVYZuPGjerUqZP+7//+T/fee6/V5Wj48OFasGDBOW/xAqg/eMYOACzw5ptvKjQ0VLfccovVpQDwIDxjBwC16NNPP9WWLVv0xhtvaPTo0c5GCwBwB4IdANSihx9+WJmZmbr++us1efJkq8sB4GF4xg4AAMBD8IwdAACAhyDYAQAAeIh6/Yydw+HQwYMHFRYWVmtLDwEAANQm0zSVl5enxo0by2a78JhcvQ52Bw8eVGJiotVlAAAA1Lh9+/YpISHhgufU62AXFhYmqewbDQ8Pt7gaAAAA98vNzVViYqIz91xIvQ525bdfw8PDCXYAAMCjVeWxM5onAAAAPATBDgAAwEMQ7AAAADwEwQ4AAMBDEOwAAAA8BMEOAADAQxDsAAAAPATBDgAAwEMQ7AAAADwEwQ4AAMBDEOwAAAA8BMEOAADAQxDsAAAAPATB7iJM09SerAIt2Zwhu8O0uhwAAIDzIthdhMOUBrz0rf777XXae7TA6nIAAADOi2B3ET42Q63iwiRJ6Rl5FlcDAABwfgS7Kkg5Fey2EuwAAEAdRrCrgpR4RuwAAEDdR7CrgtT4cEnStkyCHQAAqLsIdlXQKj5UkrTnaIFOltgtrgYAAODcCHZV0DA0QFEh/nKY0vbMfKvLAQAAOCeCXRUYhnFGA0WuxdUAAACcG8GuimigAAAAdR3BropSy4MdDRQAAKCOIthVUStG7AAAQB1HsKui8tUnDucV6XhBscXVAAAAVEawq6LQAF8lRgVJYgUKAABQNxHsXJASVzZRcTqdsQAAoA4i2LmABgoAAFCXEexcUD7lCbdiAQBAXUSwc0F5sNuWkSfTNC2uBgAAoCKCnQuSY0Lk52OooNiu/cdPWF0OAABABQQ7F/j52NS8Yagk5rMDAAB1D8HORTRQAACAuopg56KU+LIpT2igAAAAdQ3BzkXOETvmsgMAAHUMwc5F5WvG7jpSoOJSh8XVAAAAnEawc1HjiECFBfqq1GFqV1a+1eUAAAA4EexcZBiGUuLKb8fynB0AAKg7CHbVwAoUAACgLiLYVcPpBgqCHQAAqDsIdtVQPuUJwQ4AANQlBLtqKH/G7kD2CeWdLLG4GgAAgDIEu2qICPZTfHigJGkbK1AAAIA6gmBXTTRQAACAuoZgV000UAAAgLqGYFdNjNgBAIC6hmBXTa1ONVBsy8yTaZoWVwMAAECwq7YWsaHysRnKLizR4bwiq8sBAAAg2FVXoJ+PkqKDJXE7FgAA1A0Eu0uQ6pyoONfiSgAAAAh2l4QGCgAAUJcQ7C5BClOeAACAOoRgdwnKlxbbfjhfdgedsQAAwFoEu0vQNCpYQX4+Ki51aM/RAqvLAQAAXo5gdwlsNkOt4kIlcTsWAABYj2B3iWigAAAAdQXB7hKlMOUJAACoIwh2lyiVzlgAAFBHEOwuUfmasXuPFepEsd3iagAAgDcj2F2ihmEBig7xl2lK2w8zagcAAKxjabCz2+2aNGmSkpOTFRQUpObNm+vpp5+WadavOeFooAAAAHWBr5UfPn36dL366quaN2+e2rZtq7Vr12rEiBGKiIjQmDFjrCzNJSnxYfp+51GeswMAAJayNNh9//33Gjx4sG644QZJUlJSkt59912tXr3ayrJcRgMFAACoCyy9FdurVy8tW7ZM27ZtkyRt3LhRK1eu1HXXXXfO84uKipSbm1vhqy4ob6BIzyTYAQAA61g6YjdhwgTl5uYqNTVVPj4+stvteuaZZzR06NBznj9t2jRNnjy5lqu8uPJgdySvSMcKihUV4m9xRQAAwBtZOmL3r3/9S/Pnz9c777yj9evXa968eXr++ec1b968c54/ceJE5eTkOL/27dtXyxWfW0iAr5pGBUuStjJRMQAAsIilI3aPP/64JkyYoDvvvFOS1L59e+3du1fTpk3TsGHDKp0fEBCggICA2i6zSlLiw/TbsUKlZ+SpV/MYq8sBAABeyNIRu8LCQtlsFUvw8fGRw+GwqKLqo4ECAABYzdIRu0GDBumZZ55R06ZN1bZtW23YsEEzZ87U/fffb2VZ1cJcdgAAwGqWBrtZs2Zp0qRJeuihh3T48GE1btxYI0eO1JNPPmllWdWScqqBYntmnhwOUzabYXFFAADA2xhmfVvm4Qy5ubmKiIhQTk6OwsPDLa2lxO5Q2ye/ULHdoRV/vkqJp5opAAAALoUreYe1Yt3Ez8em5rGhkrgdCwAArEGwc6PTDRRMeQIAAGofwc6NaKAAAABWIti5UQpTngAAAAsR7NyovDN2d1aBikrtFlcDAAC8DcHOjRpFBCos0FelDlO7jhRYXQ4AAPAyBDs3MgyDFSgAAIBlCHZuRgMFAACwCsHOzVLiyyYOZMoTAABQ2wh2bsatWAAAYBWCnZu1ii0LdgdzTir3ZInF1QAAAG9CsHOziGA/NYoIlCRtY9QOAADUIoJdDaCBAgAAWIFgVwNYgQIAAFiBYFcDaKAAAABWINjVgFanlhZLz8yTaZoWVwMAALwFwa4GtIgNlY/NUM6JEmXmFlldDgAA8BIEuxoQ4Ouj5JgQSdJWJioGAAC1hGBXQ2igAAAAtY1gV0NS4wh2AACgdhHsaghz2QEAgNpGsKsh5cFux5F8ldodFlcDAAC8AcGuhiRGBivY30fFpQ7tOVpodTkAAMALEOxqiM1mqCXP2QEAgFpEsKtBpxsomPIEAADUPIJdDaKBAgAA1CaCXQ1yrhmbSbADAAA1j2BXg1qdCna/HStUYXGpxdUAAABPR7CrQTGhAYoJ9ZdpStsz860uBwAAeDiCXQ1jaTEAAFBbCHY1LCUuXBINFAAAoOYR7GrY6QYKpjwBAAA1i2BXw1pxKxYAANQSgl0NaxUXKsOQsvKLlZVfZHU5AADAgxHsaliwv6+aRgVLkrYxagcAAGoQwa4WpMSxAgUAAKh5BLtakMpzdgAAoBYQ7GpBSvypKU9YWgwAANQggl0tSIkPlSRtz8yTw2FaXA0AAPBUBLtakBQdIn9fmwqL7dp//ITV5QAAAA9FsKsFvj42tWhYNmq3NYOJigEAQM0g2NUSGigAAEBNI9jVkpRTwY4GCgAAUFMIdrUkhRE7AABQwwh2taQ82O3OKlBRqd3iagAAgCci2NWS+PBAhQf6yu4wtfNwgdXlAAAAD0SwqyWGYSj11ETF6Zl0xgIAAPcj2NUiZwMFz9kBAIAaQLCrRTRQAACAmkSwq0XMZQcAAGoSwa4WtYwrC3aHck4q50SJxdUAAABPQ7CrRRFBfmocEShJ2sZExQAAwM0IdrWMBgoAAFBTCHa1LKV8ypMMpjwBAADuRbCrZTRQAACAmkKwq2Wt4k7fijVN0+JqAACAJyHY1bLmsSHysRnKO1mqjNyTVpcDAAA8CMGulgX4+uiymBBJNFAAAAD3IthZgBUoAABATSDYWYAGCgAAUBMIdhYon/KEW7EAAMCdCHYWSDnVGbvzcL5K7Q6LqwEAAJ6CYGeBhMggBfv7qNju0J6jBVaXAwAAPATBzgI2m1FhPjsAAAB3sDzYHThwQPfcc4+io6MVFBSk9u3ba+3atVaXVeNooAAAAO7ma+WHHz9+XL1799ZVV12lzz//XA0bNtT27dsVGRlpZVm1onzKE0bsAACAu1Qr2Nntdi1cuFC//vqrJKl169a66aab5Ovr2uWmT5+uxMREzZkzx7kvOTm5OiXVO8xlBwAA3M3lW7G//PKLWrVqpWHDhmnhwoVauHChhg8frpYtW2rz5s0uXWvRokXq1q2bbr/9dsXGxqpz58568803z3t+UVGRcnNzK3zVV+Wdsb8dK1RhcanF1QAAAE/gcrD7r//6L7Vt21b79+/X+vXrtX79eu3bt08dOnTQAw884NK1du3apVdffVUtW7bUF198oQcffFBjxozRvHnzznn+tGnTFBER4fxKTEx0tfw6Izo0QDGhAZKkbZn5FlcDAAA8gWGapunKG4KCgrR27Vq1bdu2wv7Nmzere/fuOnHiRJWv5e/vr27duun777937hszZozWrFmjVatWVTq/qKhIRUVFzu3c3FwlJiYqJydH4eHhrnwbdcI9//ujVu7I0vRb2+uO7k2tLgcAANRBubm5ioiIqFLecXnErlWrVsrMzKy0//Dhw2rRooVL12rUqJHatGlTYV/r1q3122+/nfP8gIAAhYeHV/iqz2igAAAA7uRysJs2bZrGjBmjBQsWaP/+/dq/f78WLFigcePGafr06S49/9a7d2+lp6dX2Ldt2zY1a9bM1bLqJRooAACAO7ncFTtw4EBJ0pAhQ2QYhiSp/G7uoEGDnNuGYchut1/wWo888oh69eqlqVOnasiQIVq9erXeeOMNvfHGG66WVS+VN1AQ7AAAgDu4HOy+/vprt3149+7dtXDhQk2cOFFTpkxRcnKyXnzxRQ0dOtRtn1GXtYoLk2FIRwuKlZVf5GymAAAAqA6Xg92VV17p1gIGDhzoHAX0NkH+PmoWFaw9RwuVnpGnmBYEOwAAUH3VXnmisLBQv/32m4qLiyvs79ChwyUX5U1S4sO052ihtmbkqXeLGKvLAQAA9ZjLwe7IkSMaMWKEPv/883Mev9hzdagoJT5cX/ySqfSM+jvZMgAAqBtc7oodN26csrOz9eOPPyooKEhLlizRvHnz1LJlSy1atKgmavRoqXTGAgAAN3F5xO6rr77SJ598om7duslms6lZs2a69tprFR4ermnTpumGG26oiTo9VqtTnbHbMvPlcJiy2QyLKwIAAPWVyyN2BQUFio2NlSRFRkbqyJEjkqT27dtr/fr17q3OCyRFB8vf16YTJXbtO15odTkAAKAecznYpaSkOCcV7tixo15//XUdOHBAr732mho1auT2Aj2dr49NLWNDJbECBQAAuDQuB7uxY8fq0KFDkqS0tDR9/vnnatq0qV5++WVNnTrV7QV6A1agAAAA7uDyM3b33HOP83XXrl21d+9ebd26VU2bNlVMDNN1VAcNFAAAwB2qPY9dueDgYHXp0sUdtXitlPhwSdJWpjwBAACXoMrBLjs7W++++64efPBBSdLQoUN14sQJ53EfHx+9+eabatCggduL9HTla8buOVqokyV2Bfr5WFwRAACoj6r8jN2bb76plStXOrcXLVokm82miIgIRURE6Oeff9aLL75YEzV6vLjwAEUE+cnuMLXzSL7V5QAAgHqqysFuwYIFGjFiRIV9M2bM0Jw5czRnzhxNmzZNn3zyidsL9AaGYdBAAQAALlmVg92uXbuUkpLi3E5JSZG/v79zu2PHjtq+fbt7q/MiNFAAAIBLVeVgV1BQoJycHOf22rVrlZCQUOG4w+Fwb3VepHzEjrnsAABAdVU52F122WUXXFli7dq1Sk5OdktR3ogROwAAcKmqHOxuvvlm/b//9/+UmZlZ6VhGRobS0tJ08803u7U4b9LyVGdsRu5J5RSWWFwNAACoj6o83cmf//xnffjhh2rZsqXuvfdetWrVSpKUnp6ut99+W02aNNETTzxRY4V6uvBAPzVpEKQD2SeUnpmny5OjrC4JAADUM1UOdmFhYfruu+80ceJEvfvuu8rOzpYkNWjQQHfffbemTp2qsLCwmqrTK6TEh5UFu4xcgh0AAHCZSytPREZG6rXXXtOrr76qI0eOSJIaNmwowzBqpDhvkxIfpq+2HqaBAgAAVEu1lhQzDEOxsbHursXr0UABAAAuRZWbJ1DzWp1qoEjPzJNpmhZXAwAA6huCXR3SvGGofG2G8k6W6lDOSavLAQAA9QzBrg7x97XpsoYhkrgdCwAAXOdysNu1a1dN1IFTUuLDJbECBQAAcJ3Lwa5Fixa66qqr9Pbbb+vkSW4XutvpBopciysBAAD1jcvBbv369erQoYPGjx+v+Ph4jRw5UqtXr66J2rxSShxrxgIAgOpxOdh16tRJL730kg4ePKh//vOfOnTokH7/+9+rXbt2mjlzpnN+O1RPyqkRu11HClRid1hcDQAAqE+q3Tzh6+urW265RR988IGmT5+uHTt26LHHHlNiYqLuu+8+HTp0yJ11eo0mDYIU4u+jYrtDe7IKrC4HAADUI9UOdmvXrtVDDz2kRo0aaebMmXrssce0c+dOLV26VAcPHtTgwYPdWafXsNkMtYrndiwAAHCdyytPzJw5U3PmzFF6erquv/56/d///Z+uv/562WxlGTE5OVlz585VUlKSu2v1GqnxYdrwW7bSM/I0qKPV1QAAgPrC5WD36quv6v7779fw4cPVqFGjc54TGxurf/zjH5dcnLeigQIAAFSHy8Fu6dKlatq0qXOErpxpmtq3b5+aNm0qf39/DRs2zG1FepvyuezSM5nyBAAAVJ3Lz9g1b95cWVlZlfYfO3ZMycnJbinK25V3xu47dkIFRaUWVwMAAOoLl4Pd+Ranz8/PV2Bg4CUXBCkqxF8NwwIkSdsyuR0LAACqpsq3YsePHy9JMgxDTz75pIKDg53H7Ha7fvzxR3Xq1MntBXqr1PgwHckrUnpGnjo3jbS6HAAAUA9UOdht2LBBUtmI3c8//yx/f3/nMX9/f3Xs2FGPPfaY+yv0UilxYVqxPYsGCgAAUGVVDnZff/21JGnEiBF66aWXFB4eXmNF4fRzdukEOwAAUEUud8XOmTOnJurAWZzBLjNPpmnKMAyLKwIAAHVdlYLdLbfcorlz5yo8PFy33HLLBc/96KOP3FKYt2sZGybDkI4VFCsrv9jZTAEAAHA+VQp2ERERzhGjiIiIGi0IZYL8fZQUHaLdWQVKz8gj2AEAgIuqUrA78/Yrt2JrT0pcmHZnFWhrRq5+3zLG6nIAAEAd5/I8dqg9NFAAAABXVGnErnPnzlV+eH/9+vWXVBBOSz2jgQIAAOBiqhTsbrrpphouA+fS6lSw25aZJ7vDlI+NzlgAAHB+VQp2aWlpNV0HziEpOkQBvjadLHFo37FCJcWEWF0SAACow3jGrg7zsRlqGRcqSaxAAQAALsrlYGe32/X888/r8ssvV3x8vKKioip8wb1S4spW+KCBAgAAXIzLwW7y5MmaOXOm7rjjDuXk5Gj8+PG65ZZbZLPZ9NRTT9VAid7tdANFrsWVAACAus7lYDd//ny9+eabevTRR+Xr66u77rpL//u//6snn3xSP/zwQ03U6NXKpzzhViwAALgYl4NdRkaG2rdvL0kKDQ1VTk6OJGngwIH697//7d7q4Ax2e7IKdLLEbnE1AACgLnM52CUkJOjQoUOSpObNm+s///mPJGnNmjUKCGDZK3eLDQtQg2A/OUxpx+F8q8sBAAB1mMvB7uabb9ayZcskSQ8//LAmTZqkli1b6r777tP999/v9gK9nWEYSoljBQoAAHBxVZrH7kzPPvus8/Udd9yhpk2batWqVWrZsqUGDRrk1uJQJjU+TD/uPsYKFAAA4IJcDnZn69mzp3r27OmOWnAeKfFlU57QQAEAAC6kWsHu4MGDWrlypQ4fPiyHw1Hh2JgxY9xSGE5LiS+bpDg9gylPAADA+bkc7ObOnauRI0fK399f0dHRMozT65cahkGwqwGtTj1jl5lbpOzCYjUI9re4IgAAUBe5HOwmTZqkJ598UhMnTpTNxopktSEs0E9NGgTpQPYJpWfkqcdl0VaXBAAA6iCXk1lhYaHuvPNOQl0tO70CBc/ZAQCAc3M5nf3xj3/UBx98UBO14AJYgQIAAFyMy7dip02bpoEDB2rJkiVq3769/Pz8KhyfOXOm24rDaeXBjrnsAADA+VQr2H3xxRdKSUmRpErNE6gZ5cFuW0aeTNPkZw0AACpxOdi98MIL+uc//6nhw4fXQDk4n8tiQuVrM5RXVKqDOSfVpEGQ1SUBAIA6xuVn7AICAtS7d++aqAUX4O9rU/OGzGcHAADOz+VgN3bsWM2aNcvthTz77LMyDEPjxo1z+7U9BQ0UAADgQly+Fbt69Wp99dVXWrx4sdq2bVupeeKjjz5yuYg1a9bo9ddfV4cOHVx+rzdJiQ+TNtJAAQAAzs3lYNegQQPdcsstbisgPz9fQ4cO1Ztvvqm//vWvbruuJ0qlMxYAAFyAy8Fuzpw5bi1g1KhRuuGGG3TNNdcQ7C6ifGmxnUfyVWJ3yM+HSaIBAMBpLgc7d3rvvfe0fv16rVmzpkrnFxUVqaioyLmdm+tdTQQJkUEKDfBVflGpdmcVOIMeAACAVMVg16VLFy1btkyRkZHq3LnzBedQW79+fZU+eN++fRo7dqyWLl2qwMDAKr1n2rRpmjx5cpXO9USGYahVXKjW/5atrRl5BDsAAFBBlYLd4MGDFRAQ4Hztjslx161bp8OHD6tLly7OfXa7Xd9++61mz56toqIi+fj4VHjPxIkTNX78eOd2bm6uEhMTL7mW+iQlPlzrf8sum/KkY2OrywEAAHVIlYJdWlqa8/VTTz3llg+++uqr9fPPP1fYN2LECKWmpuqJJ56oFOqksjn0ygOmt6KBAgAAnI/Lz9hddtllWrNmjaKjoyvsz87OVpcuXbRr164qXScsLEzt2rWrsC8kJETR0dGV9uM05rIDAADn43Jb5Z49e2S32yvtLyoq0v79+91SFM4v5dRzdfuPn1B+UanF1QAAgLqkyiN2ixYtcr7+4osvFBER4dy22+1atmyZkpOTL6mY5cuXX9L7vUFkiL9iwwJ0OK9I2zLz1KVppNUlAQCAOqLKwe6mm26SVNaZOWzYsArH/Pz8lJSUpBdeeMGtxeHcUuLDdDivSOkZBDsAAHBalYOdw+GQJCUnJ2vNmjWKiYmpsaJwYanxYVqxPYsGCgAAUIHLzRO7d++uiTrggpT4cEnS1gzvmqAZAABcWLXWpFq2bJkGDhyo5s2bq3nz5ho4cKC+/PJLd9eG8yhvoEjPyJNpmhZXAwAA6gqXg93f//53DRgwQGFhYRo7dqzGjh2r8PBwXX/99XrllVdqokacpWVcqGyGdLywREfyiy7+BgAA4BVcvhU7depU/e1vf9Po0aOd+8aMGaPevXtr6tSpGjVqlFsLRGWBfj5Kig7RrqwCpWfkKTasakuyAQAAz+byiF12drYGDBhQaX+/fv2Uk5PjlqJwcSmsQAEAAM7icrC78cYbtXDhwkr7P/nkEw0cONAtReHiWIECAACczeVbsW3atNEzzzyj5cuXq2fPnpKkH374Qd99950effRRvfzyy85zx4wZ475KUQFrxgIAgLMZpottlVVdXcIwjCqvG1tdubm5ioiIUE5OjsLDw2v0s+qaXUfy9YcXvlGAr01bpgyQj82wuiQAAFADXMk7zGNXTzWLDlGgn00nSxz67VihkmNCrC4JAABYrFrz2ElSVlaWsrKy3FkLXOBjM9Qytvx2LBMVAwAAF4Nddna2Ro0apZiYGMXFxSkuLk4xMTEaPXq0srOza6hEnA8NFAAA4ExVvhV77Ngx9ezZUwcOHNDQoUPVunVrSdKWLVs0d+5cLVu2TN9//70iI1mUvrbQQAEAAM5U5WA3ZcoU+fv7a+fOnYqLi6t0rF+/fpoyZYr+9re/ub1InBtz2QEAgDNV+Vbsxx9/rOeff75SqJOk+Ph4zZgx45zz26HmlK8Zu+dogU6W2C2uBgAAWK3Kwe7QoUNq27bteY+3a9dOGRkZbikKVdMwLECRwX5ymNKOw/lWlwMAACxW5WAXExOjPXv2nPf47t27FRUV5Y6aUEWGYdBAAQAAnKoc7Pr376//+Z//UXFxcaVjRUVFmjRp0jnXkEXNSo0vm6iQKU8AAIBLzRPdunVTy5YtNWrUKKWmpso0Tf3666/6+9//rqKiIr311ls1WSvOgRE7AABQrsrBLiEhQatWrdJDDz2kiRMnqnwlMsMwdO2112r27NlKTEyssUJxbq3i6IwFAABlXFpSLDk5WZ9//rmOHz+u7du3S5JatGjBs3UWKh+xO5xXpOMFxYoM8be4IgAAYBWX14qVpMjISF1++eXurgXVEBrgq4TIIO0/fkLpmXn63WXRVpcEAAAsUu21YlF3sAIFAACQCHYegQYKAAAgEew8QgpTngAAABHsPEL50mLbMvOd3coAAMD7EOw8wGUNQ+TnYyi/qFQHsk9YXQ4AALAIwc4D+PnY1LxhqCQaKAAA8GYEOw9BAwUAACDYeYgUpjwBAMDrEew8BHPZAQAAgp2HKF8zdueRfBWXOiyuBgAAWIFg5yGaNAhSWICvSh2mdmcVWF0OAACwAMHOQxiGoVbOBgomKgYAwBsR7DwIDRQAAHg3gp0HoYECAADvRrDzIOVLizGXHQAA3olg50HKb8UeyD6hvJMlFlcDAABqG8HOgzQI9ldceIAkaVtmvsXVAACA2kaw8zAp8eGSeM4OAABvRLDzMKcbKJjyBAAAb0Ow8zA0UAAA4L0Idh7GOZddZp5M07S4GgAAUJsIdh6mRWyobIaUXViiI3lFVpcDAABqEcHOwwT6+SgpJkQSt2MBAPA2BDsPxAoUAAB4J4KdB0qJK5vyhBE7AAC8C8HOA51uoGDKEwAAvAnBzgOVB7vtmfmyO+iMBQDAWxDsPFDTqGAF+tlUVOrQ3qMFVpcDAABqCcHOA/nYDLWKo4ECAABvQ7DzUKxAAQCA9yHYeagUpjwBAMDrEOw8VGp82ZQn6ZkEOwAAvAXBzkO1ig+VJO05WqATxXaLqwEAALWBYOehGoYGKCrEX6Yp7Ticb3U5AACgFhDsPJRhGGc0UDBRMQAA3oBg58FooAAAwLsQ7DxYqnNpMYIdAADegGDnwVrFM5cdAADehGDnwcpXnziSV6RjBcUWVwMAAGoawc6DhQb4KjEqSBLP2QEA4A0Idh4uJe7URMV0xgIA4PEsDXbTpk1T9+7dFRYWptjYWN10001KT0+3siSPQwMFAADew9Jg980332jUqFH64YcftHTpUpWUlKhfv34qKCiwsiyPkkIDBQAAXsPXyg9fsmRJhe25c+cqNjZW69atU58+fSyqyrOUB7ttGXlyOEzZbIbFFQEAgJpSp56xy8nJkSRFRUVZXInnSI4JkZ+PoYJiuw5kn7C6HAAAUIPqTLBzOBwaN26cevfurXbt2p3znKKiIuXm5lb4woX5+djUvGGoJDpjAQDwdHUm2I0aNUqbN2/We++9d95zpk2bpoiICOdXYmJiLVZYf9FAAQCAd6gTwW706NFavHixvv76ayUkJJz3vIkTJyonJ8f5tW/fvlqssv5KiS+b8oQGCgAAPJulzROmaerhhx/WwoULtXz5ciUnJ1/w/ICAAAUEBNRSdZ7DOWLHXHYAAHg0S4PdqFGj9M477+iTTz5RWFiYMjIyJEkREREKCgqysjSPUr5m7K4jBSoudcjft04M1AIAADez9L/wr776qnJyctS3b181atTI+fX+++9bWZbHaRwRqLBAX5U6TO3Kyre6HAAAUEMsvxWLmmcYhlLiwrR273GlZ+Qp9dQzdwAAwLNwT85LsAIFAACej2DnJU43UBDsAADwVAQ7L1E+5QnBDgAAz0Ww8xIpcWUjdgeyTyj3ZInF1QAAgJpAsPMSEcF+ig8PlCRtZwUKAAA8EsHOi9BAAQCAZyPYeREaKAAA8GwEOy/CiB0AAJ6NYOdFWsWdHrFjcmgAADwPwc6LtIgNlY/NUM6JEmXmFlldDgAAcDOCnRcJ9PNRUnSwJCmdzlgAADwOwc7LpDonKs61uBIAAOBuBDsvQwMFAACei2DnZVKY8gQAAI9FsPMy5UuLbT+cr1K7w+JqAACAOxHsvEzTqGAF+fmouNShvccKrS4HAAC4EcHOy9hshlrFhUridiwAAJ6GYOeFaKAAAMAzEey8UApTngAA4JEIdl4olc5YAAA8EsHOC5WvGbv3WKEKi0strgYAALgLwc4LNQwLUHSIv0xT2nE43+pyAACAmxDsvBQNFAAAeB6CnZdiBQoAADwPwc5L0UABAIDnIdh5qfIGCm7FAgDgOXytLgDWKA92WflFOppfpOjQAIsrqnl7sgr04fr9WrolU5IUHeqvyGB/RYWc8WeIv6JOvY4K8VeDYD8F+vlYXDkAAFVDsPNSIQG+ahoVrN+OFSo9M0+9PDTY5Z4s0b83HdKH6/Zr7d7j1bpGiL9PWeA7IwBGVdj2qxAMI4P95WMz3PydAABwcQQ7L5YSH1YW7DLy1Kt5jNXluI3dYeq7HVlasG6/vvglQ0WlDkmSzZCuaNlQt3Rposhgfx0vLNbR/GIdLyzWsYIz/iwo0bHCYh0vKFapw1RBsV0FxSe0//iJKn2+YUgRQX6KCj4d9KJC/CqNBpZvR4b4KzzQV4ZBGAQAXBqCnRdLjQ/T0i2ZHtNAseNwnhasO6CPNxxQRu5J5/6WsaG6rWuCburcRHHhgVW+nmmayj1ZquMFxc6gdzoAluhYQZGOFZToePmxwmJlF5bINKXswhJlF5ZIWQVV+ixfm6EGpwJg1Fmjg2eOBkafEQiD/LlFDACoiGDnxTxhLrvswmJ9uvGgFqw/oI37sp37GwT7aXDHxrq1a4LaN4mo1miYYRiKCPJTRJCfkhRSpfeU2h3KOVFy1mhgyRmjgadD4tFT2wXFdpU6TGXlFykrv6jK9QX62ZwjfmffJo4M8VdEkJ98vGAUMDY8QN2aRTLiCQAi2Hm1lFMNFNsy8+RwmLLVk+fCSu0Ofbv9iBas268vtxxWsb3sVquPzdBVKQ11W9cEXZUaqwDf2h/R8vWxKTo0QNGhAWoRW7X3nCyxK7uwpMLt4GMFZ90eLg+Ip/YX2x06WeLQwZyTOphz8uIf4uF+3yJGT93YRi1iw6wuBQAsRbDzYkkxIfL3samw2K79x0+oaXSw1SVd0NaMXC1Yu18f/3SwwshW60bhuq1rggZ3aqyYetgEEujno/gIH8VHVO02sWmWPfdXHvIq3yYu+8o5UXZb2JOZkn7al62VO7I04MUVGtE7SWOubqmwQD+rSwMASxDsvJifj03NY0P166FcpWfm1clgdzS/SIs2HtSCdfv1y8Fc5/7oEH8N7tREt3ZtoraNIyyssPYZhqHQAF+FBvgqMaru/c5q229HCzVl8RZ9+Wum3lyxWx//dFB/uT5VN3Vqwu1ZAF6HYOflUuPDyoJdRq6ubRNndTmSpOJSh75OP6wF6/br662HVeooG3by8zF0dWqcbuuaoCtTGsrPh/m1ITWNDtb/Duumr9MPa8qnW7Q7q0CPvL9R83/4TZMHt/W64A/AuxHsvFxdaaAwTVO/HMzVgnX7tWjjQR0rKHYe65AQodu6JmhQh8aKDPG3sErUZVelxKpX82j9Y+VuzVq2Q2v3HtegWSs1tEczPdqvlRoE83cHgOcj2Hm5FIvXjD2cd1KfbCi71ZqeebqG2LAA3dy5iW7tmuBcJQO4mABfHz3Ut4Vu7txEUz/bqk83HtRbP+zV4k0H9Xj/VN3RPZHJowF4NIKdlyvvjN2VVaCiUnutdJKeLLFr2a+HtWDdPn27PUv2U7da/X1t6tem7Fbr71vEyJdbraimRhFBmnVXZ919eVM9tegXpWfm6S8Lf9a7q8tuz3ZpGml1iQBQIwh2Xq5RRKDCAn2Vd7JUu44UqHWj8Br5HNM09dO+bC1Yt1+fbjyo3JOlzmNdmjbQrV0TNLBDY0UE0c0I9+nZPFr/HvN7vfXDXs38zzb9fCBHt/z9e93WNUFPDEhVw7D610UNABdCsPNyhmEoNT5Ma/YcV3pGntuD3aGcE1q44YAWrNuvXUdOr8LQKCJQt3Rpolu6JKh5w1C3fiZwJl8fm0b0TtbADo01Y8lWfbBuf9lyc5szNO7aVrqvZzMacQB4DIIdlHIq2LmrgeJEsV3/2ZKhBev2a+WOLOdcaoF+Nl3XrpFu7ZKgns2jedYJtaphWICeu72j7urRVGmf/KKfD+To6cVb9P6a3/TUjW09ar1kAN6LYAelxJeN0qVn5F7kzPMzTVNr9x7Xh+v2a/GmQ8ovOn2r9fLkKN3WNUHXt2+k0AD+ysFaXZpG6uNRvfWvtfs0Y8lWbcvM191v/qgbOjTS/1zfWo0bBFldIgBUG/+VhVIvoTN2//FCfbT+gD5cv197jxY69ydGBemWzgm6tUtCnZz4GN7Nx2borsub6rp28Zq5dJve/mGv/r3pkL769bBG/6GF/uuKZEuWpAOAS0Wwg1qdWl/zYM5J5ZwouWgDQ0FRqT7fnKEP1+3Xql1HnftD/H10fftGuq1rgronRdWbtWfhvRoE+2vK4Ha6s3tTpS3arDV7juu5L9L1wdp9ShvUVlelVnHBXwCoIwh2UESwnxpFBOpQzkltz8xTt6SoSuc4HKZ+2H1UH647oM83H1JhsV2SZBhSr+bRurVLgga0i1ewP3+lUP+0aRyuf43sqU9+Oqipn/2qPUcLNWLuGl3TOlaTBrZRs+gQq0sEgCrhv8KQVNZAcSjnpLZmVAx2e7IK9NH6/fpw/QEdyD7h3J8cE6JbuzTRzV0S1IRnkuABDMPQTZ2b6Jo2cZq1bLv+sXK3vvz1sL7dnqWRfS7TQ31bKMif27MA6jaCHSSVBbvl6UeUnpGn3JMl+mzTIX24fr/W7DnuPCcs0FcDOzTWbV2bqEvTSBZYh0cKDfDVxOtb6/ZuiZr86S9asT1Ls77aoY/WH9D/3NBa17WL5+8+gDrLMM3yySjqn9zcXEVERCgnJ0fh4TUzsa63WLhhvx55f6PCA31VbHfoZIlDkmQzpCtaNtStXRPUr02cAv0YsYD3ME1TX/ySqacXb3GOWPduEa2nBrVVS5a6A1BLXMk7BDtIkn49lKvrXlrh3G4ZG6pbuybo5s5NFBceaGFlgPVOFNv16jc79do3O1Vc6pCvzdDwXkkae01LhQWyWgqAmkWwg8tM09SLX25X7skS3dSpiTokRHC7CTjLb0cL9fS/t2jplkxJUkxogCZel6qbOzehCxxAjSHYAUANWp5+WJM/3aLdWWXL5HVtFqnJN7ZVuyYRFlcGwBMR7ACghhWV2vXPlXs066vtKiy2yzCkuy9vqsf6pSgyxN/q8gB4EFfyDitfA0A1BPj66MG+zfXVo311Y8fGMk1p/o+/6aoXlmv+j3tld9Tb/2cGUI8xYgcAbvDDrqNK++QXpWeWLc3Xrkm4Jt/YTl2bRVpcGYD6jluxAGCBUrtDb/+wVy8s3aa8k6WSpFu7JOiJ61IUG0Z3OYDq4VYsAFjA18em4b2T9fVjfTWkW4Ik6cP1+3X189/of1fsUondYXGFADwdI3YAUEM2/HZcaYt+0ab9OZLK5oecfGNb9WoRY3FlAOoTbsUCQB3hcJj619p9mvFFuo4VFEuSbmjfSH+5oTXrLAOoEoIdANQxOYUlmrk0XW/9sFcOUwry89Goq5rrv664jKX6AFwQwQ4A6qgtB3P11KJftHrPMUlSs+hgpQ1qoz+kxllcGYC6imAHAHWYaZpatPGgpn72qzJziyRJf0iN1ZMD2ygpJsTi6gDUNQQ7AKgH8otKNeur7frnyt0qsZvy97HpgT6X6aGrmivY39fq8gDUEQQ7AKhHdhzO1+RPf9GK7VmSpMYRgfqfG9ro+vbxMgzD4uoAWI1gBwD1jGma+s+WTE35dIsOZJ+QJPVqHq3JN7ZVy7gwi6sDYKV6N0HxK6+8oqSkJAUGBqpHjx5avXq11SUBQK0yDEP928Zr2aNXauzVLRXga9P3O4/qupdW6OnFW3Qw+4TyTpaolEmOAVyA5SN277//vu677z699tpr6tGjh1588UV98MEHSk9PV2xs7AXfy4gdAE+171ihnl68Rf/ZklnpmL+PTYF+NgX5+yjY31eBfj4KOrUd5Od76k/bGcd8FORvU5C/b9nr8m3nuT4K9vcpO/fUto+NW8BAXVGvbsX26NFD3bt31+zZsyVJDodDiYmJevjhhzVhwoQLvpdgB8DTfbPtiKZ99qvSM/NUm/+29ve1OUPg2aHP+Wf567O3L7b/1D5bDYdH0zTlMCW7w5TDLP86tX1qn9005XCo7LXjrHNOvafsfDnPN01TdkfZOeapfWWvy/ZVOMf5+vQ5DtOUzTBksxmyGZKPzZBhGPIxDPnYdMZrQ8ap4z7GqXNs5z/HZpR9+Zy67unXhmw2yccwnJ97znMM8UxnHeVK3rG07aq4uFjr1q3TxIkTnftsNpuuueYarVq1qtL5RUVFKioqcm7n5ubWSp0AYJUrWzXUla0ayjRNFZU6dLLErsJiu06U2HWi2F5h2/n6jOMnSk4dKz7jfSUVzznzGuWKSx0qLnUo50RJjX1vAb62CoEv0NdHhqHTIak8jDnOFcDOCmPl55hyvq6/T5BbpzzsOUOn83XFAFgpMJafc8Z7bUZZ6PQG7z/QU0H+dWOicUuDXVZWlux2u+LiKk7MGRcXp61bt1Y6f9q0aZo8eXJtlQcAdYZhGAr0Kxs5axBcM59RHh7PDI4nKoXBUp0odpzaLj31p+PU/vJzHWccOytglpx+RrCo1KGiUoeyVXPh8WIM48yRrCqMap19jnHGqJpz5O3sUbizR950jmB6Ksg6RxJ1+vUZI4mnA+yp95w9sniuc069rkrQLa9LDlKxKxx16P8i6tVESRMnTtT48eOd27m5uUpMTLSwIgDwHGeGx5ricJg6WXqOUcViuwpPjRhWCEPGGYHpHEGqcvA6NWpk6IzXpwLZOW5XetOtR/O8o53nGB0tD4gVguMFblmfPap66ha0twjwrRO9qJIsDnYxMTHy8fFRZmbFh4MzMzMVHx9f6fyAgAAFBATUVnkAADez2QwF+/syAbMFDMOQr4/3BFlvZWnE9Pf3V9euXbVs2TLnPofDoWXLlqlnz54WVgYAAFD/WP6/TOPHj9ewYcPUrVs3XX755XrxxRdVUFCgESNGWF0aAABAvWJ5sLvjjjt05MgRPfnkk8rIyFCnTp20ZMmSSg0VAAAAuDDL57G7FMxjBwAAPF29W1IMAAAAl45gBwAA4CEIdgAAAB6CYAcAAOAhCHYAAAAegmAHAADgIQh2AAAAHoJgBwAA4CEIdgAAAB6CYAcAAOAhLF8r9lKUr4aWm5trcSUAAAA1ozznVGUV2Hod7PLy8iRJiYmJFlcCAABQs/Ly8hQREXHBcwyzKvGvjnI4HDp48KDCwsJkGEaNfU5ubq4SExO1b9++iy6+i7qP36fn4HfpWfh9eg5+l+5lmqby8vLUuHFj2WwXfoquXo/Y2Ww2JSQk1NrnhYeH8xfUg/D79Bz8Lj0Lv0/Pwe/SfS42UleO5gkAAAAPQbADAADwEAS7KggICFBaWpoCAgKsLgVuwO/Tc/C79Cz8Pj0Hv0vr1OvmCQAAAJzGiB0AAICHINgBAAB4CIIdAACAhyDYVcErr7yipKQkBQYGqkePHlq9erXVJcFF06ZNU/fu3RUWFqbY2FjddNNNSk9Pt7osuMmzzz4rwzA0btw4q0tBNRw4cED33HOPoqOjFRQUpPbt22vt2rVWl4VqsNvtmjRpkpKTkxUUFKTmzZvr6aefrtJSWHAPgt1FvP/++xo/frzS0tK0fv16dezYUf3799fhw4etLg0u+OabbzRq1Cj98MMPWrp0qUpKStSvXz8VFBRYXRou0Zo1a/T666+rQ4cOVpeCajh+/Lh69+4tPz8/ff7559qyZYteeOEFRUZGWl0aqmH69Ol69dVXNXv2bP3666+aPn26ZsyYoVmzZlldmtegK/YievTooe7du2v27NmSypYxS0xM1MMPP6wJEyZYXB2q68iRI4qNjdU333yjPn36WF0Oqik/P19dunTR3//+d/31r39Vp06d9OKLL1pdFlwwYcIEfffdd1qxYoXVpcANBg4cqLi4OP3jH/9w7rv11lsVFBSkt99+28LKvAcjdhdQXFysdevW6ZprrnHus9lsuuaaa7Rq1SoLK8OlysnJkSRFRUVZXAkuxahRo3TDDTdU+GcU9cuiRYvUrVs33X777YqNjVXnzp315ptvWl0WqqlXr15atmyZtm3bJknauHGjVq5cqeuuu87iyrxHvV4rtqZlZWXJbrcrLi6uwv64uDht3brVoqpwqRwOh8aNG6fevXurXbt2VpeDanrvvfe0fv16rVmzxupScAl27dqlV199VePHj9df/vIXrVmzRmPGjJG/v7+GDRtmdXlw0YQJE5Sbm6vU1FT5+PjIbrfrmWee0dChQ60uzWsQ7OB1Ro0apc2bN2vlypVWl4Jq2rdvn8aOHaulS5cqMDDQ6nJwCRwOh7p166apU6dKkjp37qzNmzfrtddeI9jVQ//61780f/58vfPOO2rbtq1++uknjRs3To0bN+b3WUsIdhcQExMjHx8fZWZmVtifmZmp+Ph4i6rCpRg9erQWL16sb7/9VgkJCVaXg2pat26dDh8+rC5dujj32e12ffvtt5o9e7aKiork4+NjYYWoqkaNGqlNmzYV9rVu3VoffvihRRXhUjz++OOaMGGC7rzzTklS+/bttXfvXk2bNo1gV0t4xu4C/P391bVrVy1btsy5z+FwaNmyZerZs6eFlcFVpmlq9OjRWrhwob766islJydbXRIuwdVXX62ff/5ZP/30k/OrW7duGjp0qH766SdCXT3Su3fvSlMPbdu2Tc2aNbOoIlyKwsJC2WwVo4WPj48cDodFFXkfRuwuYvz48Ro2bJi6deumyy+/XC+++KIKCgo0YsQIq0uDC0aNGqV33nlHn3zyicLCwpSRkSFJioiIUFBQkMXVwVVhYWGVno8MCQlRdHQ0z03WM4888oh69eqlqVOnasiQIVq9erXeeOMNvfHGG1aXhmoYNGiQnnnmGTVt2lRt27bVhg0bNHPmTN1///1Wl+Y1mO6kCmbPnq3nnntOGRkZ6tSpk15++WX16NHD6rLgAsMwzrl/zpw5Gj58eO0WgxrRt29fpjuppxYvXqyJEydq+/btSk5O1vjx4/WnP/3J6rJQDXl5eZo0aZIWLlyow4cPq3Hjxrrrrrv05JNPyt/f3+ryvALBDgAAwEPwjB0AAICHINgBAAB4CIIdAACAhyDYAQAAeAiCHQAAgIcg2AEAAHgIgh0AAICHINgBAAB4CIIdgDpvz549MgxDP/30k9WlOG3dulW/+93vFBgYqE6dOlldTrUYhqGPP/7Y6jIAuBHBDsBFDR8+XIZh6Nlnn62w/+OPPz7vcm2eLi0tTSEhIUpPT9eyZcvOec6RI0f04IMPqmnTpgoICFB8fLz69++v7777rparBeAtCHYAqiQwMFDTp0/X8ePHrS7FbYqLi6v93p07d+r3v/+9mjVrpujo6HOec+utt2rDhg2aN2+etm3bpkWLFqlv3746evRotT8XAC6EYAegSq655hrFx8dr2rRp5z3nqaeeqnRb8sUXX1RSUpJze/jw4brppps0depUxcXFqUGDBpoyZYpKS0v1+OOPKyoqSgkJCZozZ06l62/dulW9evVSYGCg2rVrp2+++abC8c2bN+u6665TaGio4uLidO+99yorK8t5vG/fvho9erTGjRunmJgY9e/f/5zfh8Ph0JQpU5SQkKCAgAB16tRJS5YscR43DEPr1q3TlClTZBiGnnrqqUrXyM7O1ooVKzR9+nRdddVVatasmS6//HJNnDhRN954o/O8mTNnqn379goJCVFiYqIeeugh5efnO4/PnTtXDRo00OLFi5WSkqLg4GDddtttKiws1Lx585SUlKTIyEiNGTNGdrvd+b6kpCQ9/fTTuuuuuxQSEqImTZrolVdeOef3W27fvn0aMmSIGjRooKioKA0ePFh79uxxHl++fLkuv/xyhYSEqEGDBurdu7f27t17wWsCqF0EOwBV4uPjo6lTp2rWrFnav3//JV3rq6++0sGDB/Xtt99q5syZSktL08CBAxUZGakff/xR//3f/62RI0dW+pzHH39cjz76qDZs2KCePXtq0KBBztGv7Oxs/eEPf1Dnzp21du1aLVmyRJmZmRoyZEiFa8ybN0/+/v767rvv9Nprr52zvpdeekkvvPCCnn/+eW3atEn9+/fXjTfeqO3bt0uSDh06pLZt2+rRRx/VoUOH9Nhjj1W6RmhoqEJDQ/Xxxx+rqKjovD8Lm82ml19+Wb/88ovmzZunr776Sn/+858rnFNYWKiXX35Z7733npYsWaLly5fr5ptv1meffabPPvtMb731ll5//XUtWLCgwvuee+45dezYURs2bNCECRM0duxYLV269Jx1lJSUqH///goLC9OKFSv03XffKTQ0VAMGDFBxcbFKS0t100036corr9SmTZu0atUqPfDAA157Kx6os0wAuIhhw4aZgwcPNk3TNH/3u9+Z999/v2maprlw4ULzzH+NpKWlmR07dqzw3r/97W9ms2bNKlyrWbNmpt1ud+5LSUkxr7jiCud2aWmpGRISYr777rumaZrm7t27TUnms88+6zynpKTETEhIMKdPn26apmk+/fTTZr9+/Sp89r59+0xJZnp6ummapnnllVeanTt3vuj327hxY/OZZ56psK979+7mQw895Nzu2LGjmZaWdsHrLFiwwIyMjDQDAwPNXr16mRMnTjQ3btx4wfd88MEHZnR0tHN7zpw5piRzx44dzn0jR440g4ODzby8POe+/v37myNHjnRuN2vWzBwwYECFa99xxx3mdddd59yWZC5cuNA0TdN86623zJSUFNPhcDiPFxUVmUFBQeYXX3xhHj161JRkLl++/IL1A7AWI3YAXDJ9+nTNmzdPv/76a7Wv0bZtW9lsp//1ExcXp/bt2zu3fXx8FB0drcOHD1d4X8+ePZ2vfX191a1bN2cdGzdu1Ndff+0cKQsNDVVqaqqksufhynXt2vWCteXm5urgwYPq3bt3hf29e/d2+Xu+9dZbdfDgQS1atEgDBgzQ8uXL1aVLF82dO9d5zpdffqmrr75aTZo0UVhYmO69914dPXpUhYWFznOCg4PVvHlz53ZcXJySkpIUGhpaYd+Ffl7l2+f7HjZu3KgdO3YoLCzM+fOLiorSyZMntXPnTkVFRWn48OHq37+/Bg0apJdeekmHDh1y6ecBoOYR7AC4pE+fPurfv78mTpxY6ZjNZpNpmhX2lZSUVDrPz8+vwrZhGOfc53A4qlxXfn6+Bg0apJ9++qnC1/bt29WnTx/neSEhIVW+pjsEBgbq2muv1aRJk/T9999r+PDhSktLk1Q2jcvAgQPVoUMHffjhh1q3bp3zObgzGztq4ud1tvz8fHXt2rXSz2/btm26++67JUlz5szRqlWr1KtXL73//vtq1aqVfvjhh2p/JgD3I9gBcNmzzz6rTz/9VKtWraqwv2HDhsrIyKgQ7tw599yZIaK0tFTr1q1T69atJUldunTRL7/8oqSkJLVo0aLClythLjw8XI0bN640Jcl3332nNm3aXPL30KZNGxUUFEiS1q1bJ4fDoRdeeEG/+93v1KpVKx08ePCSP6Pc2aHrhx9+cP68ztalSxdt375dsbGxlX5+ERERzvM6d+6siRMn6vvvv1e7du30zjvvuK1eAJeOYAfAZe3bt9fQoUP18ssvV9jft29fHTlyRDNmzNDOnTv1yiuv6PPPP3fb577yyitauHChtm7dqlGjRun48eO6//77JUmjRo3SsWPHdNddd2nNmjXauXOnvvjiC40YMaJCt2hVPP7445o+fbref/99paena8KECfrpp580duzYKl/j6NGj+sMf/qC3335bmzZt0u7du/XBBx9oxowZGjx4sCSpRYsWKikp0axZs7Rr1y699dZb523oqI7vvvtOM2bM0LZt2/TKK6/ogw8+OO/3MHToUMXExGjw4MFasWKFdu/ereXLl2vMmDHav3+/du/erYkTJ2rVqlXau3ev/vOf/2j79u3nDYoArEGwA1AtU6ZMqXTrr3Xr1vr73/+uV155RR07dtTq1avP2TFaXc8++6yeffZZdezYUStXrtSiRYsUExMjSc5RNrvdrn79+ql9+/YaN26cGjRoUOF5vqoYM2aMxo8fr0cffVTt27fXkiVLtGjRIrVs2bLK1wgNDVWPHj30t7/9TX369FG7du00adIk/elPf9Ls2bMlSR07dtTMmTM1ffp0tWvXTvPnz7/gdDKuevTRR7V27Vp17txZf/3rXzVz5szzTvESHBysb7/9Vk2bNtUtt9yi1q1b649//KNOnjyp8PBwBQcHa+vWrbr11lvVqlUrPfDAAxo1apRGjhzptnoBXDrDPPuBGABAvZeUlKRx48Zp3LhxVpcCoBYxYgcAAOAhCHYAAAAegluxAAAAHoIROwAAAA9BsAMAAPAQBDsAAAAPQbADAADwEAQ7AAAAD0GwAwAA8BAEOwAAAA9BsAMAAPAQBDsAAAAP8f8BG3ksJ6C6XF8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "@jax.jit\n",
    "def compute_targ_Q(mdp: MDP, Q: np.ndarray):\n",
    "    S, A = Q.shape\n",
    "    greedy_policy = compute_greedy_policy(Q)\n",
    "    max_Q = np.sum(greedy_policy * Q, axis=1).reshape(1, 1, S)\n",
    "    next_v = np.sum(mdp.P * max_Q, axis=-1)\n",
    "    assert next_v.shape == (S, A)\n",
    "    return mdp.rew + mdp.gamma * next_v\n",
    "\n",
    "\n",
    "def do_table_value_iteration(mdp: MDP, title):\n",
    "    optimality_gaps = []\n",
    "    \n",
    "    # >>>> Do value iteration >>>>\n",
    "    Q_table = np.zeros((mdp.S, mdp.A))\n",
    "    for _ in range(10):\n",
    "        Q_table = compute_targ_Q(mdp, Q_table)\n",
    "        greedy_policy = compute_greedy_policy(Q_table)\n",
    "        optimality_gaps.append(compute_optimality_gap(mdp, greedy_policy))\n",
    "    # <<<< Do value iteration <<<<\n",
    "\n",
    "    # >>>> Plot results >>>>\n",
    "    plt.plot(optimality_gaps)\n",
    "    plt.xlabel(\"Number of Samples\")\n",
    "    plt.ylabel(\"Optimality Gap\")\n",
    "    plt.title(\"Optimality Gap\")\n",
    "    # <<<< Plot results <<<<\n",
    "\n",
    "    plt.tight_layout()\n",
    "\n",
    "do_table_value_iteration(garnet_MDP, \"GarnetMDP\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('shumi-VTLwuKSy-py3.9')",
   "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.9.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "d6b7cac5e0d2ff733f340da4d53ae5ecfef7f7ad39623f5982b029a09306b36b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}