notebooks/BANDIT_utils.ipynb

{
    "cells": [
        {
            "attachments": {},
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "# バンディットの実験に便利なコード\n",
                "\n",
                "タイトルの通りです。よく使う関数をまとめます。\n",
                "\n",
                "## 確率的バンディット\n",
                "\n",
                "参考：\n",
                "* [Bandit Algorithms](https://tor-lattimore.com/downloads/book/book.pdf)\n",
                "\n",
                "確率的バンディット（の環境）は次で定義されます。\n",
                "\n",
                "* $\\mathcal{A}$：選択可能な行動の集合\n",
                "    * 有限な場合はアームの数＝$A$とします\n",
                "* $\\nu=\\left(P_a: a \\in \\mathcal{A}\\right)$：各行動に設定された報酬の分布\n",
                "\n",
                "また、$\\mathcal{A}$が有限であり、次のような性質を満たす分布$\\mathcal{M}_a$が存在する場合、そのバンディットの集合は「Unstructuredである」と言います。\n",
                "\n",
                "$$\n",
                "\\mathcal{E}=\\left\\{\\nu=\\left(P_a: a \\in \\mathcal{A}\\right): P_a \\in \\mathcal{M}_a \\text { for all } a \\in \\mathcal{A}\\right\\}=\\times_{a \\in \\mathcal{A}} \\mathcal{M}_a\n",
                "$$\n",
                "\n",
                "以下のUnstructuredな確率的バンディットはよく出てきます。\n",
                "\n",
                "|           Name          |                       Symbol                       |                                                                     Definition                                                                    |\n",
                "|:-----------------------:|:--------------------------------------------------:|:-------------------------------------------------------------------------------------------------------------------------------------------------:|\n",
                "|        Bernoulli        |            $\\mathcal{E}_{\\mathcal{B}}^A$           |                                   $\\left\\{\\left(\\mathcal{B}\\left(\\mu_i\\right)\\right)_i: \\mu \\in[0,1]^A\\right\\}$                                   |\n",
                "|         Uniform         |            $\\mathcal{E}_{\\mathcal{U}}^A$           |        $\\left\\{\\left(\\mathcal{U}\\left(a_i, b_i\\right)\\right)_i: a, b \\in \\mathbb{R}^A\\right.$ with $a_i \\leq b_i$ for all $\\left.i\\right\\}$       |\n",
                "|  Gaussian (known var.)  | $\\mathcal{E}_{\\mathcal{N}}^A\\left(\\sigma^2\\right)$ |                           $\\left\\{\\left(\\mathcal{N}\\left(\\mu_i, \\sigma^2\\right)\\right)_i: \\mu \\in \\mathbb{R}^A\\right\\}$                           |\n",
                "| Gaussian (unknown var.) |            $\\mathcal{E}_{\\mathcal{N}}^A$           | $\\left\\{\\left(\\mathcal{N}\\left(\\mu_i, \\sigma_i^2\\right)\\right)_i: \\mu \\in \\mathbb{R}^A \\; \\text{ and } \\; \\sigma^2 \\in[0, \\infty)^A\\right\\}$ |\n",
                "|     Finite variance     |  $\\mathcal{E}_{\\mathbb{V}}^A\\left(\\sigma^2\\right)$ |                             $\\left\\{\\left(P_i\\right)_i: \\mathbb{V}_{X \\sim P_i}[X] \\leq \\sigma^2\\; \\forall i\\right\\}$                             |\n",
                "\n",
                "基本のバンディットをいくつか実装しておきます。"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 1,
            "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"
                    ]
                },
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "平均報酬 [0.72542226 0.58597934 0.9922459  0.51303005 0.8881017 ]\n",
                        "アームを引いて出た報酬 1.0\n",
                        "平均報酬 [0.08185661 0.58915615 0.24576998 0.5304918  0.2302388 ]\n",
                        "アームを引いて出た報酬 0.16982543\n"
                    ]
                }
            ],
            "source": [
                "import jax.numpy as jnp\n",
                "from jax.random import PRNGKey\n",
                "import jax\n",
                "from typing import NamedTuple\n",
                "from typing_extensions import Self\n",
                "\n",
                "key = PRNGKey(0)\n",
                "\n",
                "A = 5  # アームの数\n",
                "A_set = jnp.arange(A)  # アームの集合\n",
                "\n",
                "\n",
                "# BernoulliBanditを実装します\n",
                "key, _ = jax.random.split(key)\n",
                "rew_mean = jax.random.uniform(key=key, shape=(A,))\n",
                "\n",
                "\n",
                "class BernoulliBandit(NamedTuple):\n",
                "    A_set: jnp.array  # 行動集合\n",
                "    rew_mean: jnp.array  # 各アームの報酬の期待値\n",
                "\n",
                "    @property\n",
                "    def A(self) -> int:  # アームの数\n",
                "        return len(self.A_set)\n",
                "   \n",
                "    @staticmethod\n",
                "    def pull(bandit: Self, key: PRNGKey, arm: int):\n",
                "        key, _ = jax.random.split(key)\n",
                "        rew = jax.random.bernoulli(key, bandit.rew_mean[arm])\n",
                "        return key, jnp.array(rew, dtype=jnp.float32)\n",
                "\n",
                "bernoulli_bandit = BernoulliBandit(A_set, rew_mean)\n",
                "print(\"平均報酬\", bernoulli_bandit.rew_mean)\n",
                "print(\"アームを引いて出た報酬\", bernoulli_bandit.pull(bernoulli_bandit, key, 2)[1])\n",
                "\n",
                "\n",
                "# GaussianBanditを実装します\n",
                "key_mean, key_var, key = jax.random.split(key, num=3)\n",
                "rew_mean = jax.random.uniform(key=key_mean, shape=(A,))\n",
                "rew_var = jax.random.uniform(key=key_var, shape=(A,))\n",
                "\n",
                "\n",
                "class GaussianBandit(NamedTuple):\n",
                "    A_set: jnp.array  # 行動集合\n",
                "    rew_mean: jnp.array  # 各アームの報酬の期待値\n",
                "    rew_var: jnp.array  # 各アームの報酬の期待値\n",
                "\n",
                "    @property\n",
                "    def A(self) -> int:  # アームの数\n",
                "        return len(self.A_set)\n",
                "\n",
                "    @staticmethod\n",
                "    def pull(bandit: Self, key: PRNGKey, arm: int):\n",
                "        key, _ = jax.random.split(key)\n",
                "        rew = jax.random.normal(key) * jnp.sqrt(bandit.rew_var[arm]) + bandit.rew_mean[arm]\n",
                "        return key, jnp.array(rew, dtype=jnp.float32)\n",
                "\n",
                "\n",
                "gaussian_bandit = GaussianBandit(A_set, rew_mean, rew_var)\n",
                "print(\"平均報酬\", gaussian_bandit.rew_mean)\n",
                "print(\"アームを引いて出た報酬\", gaussian_bandit.pull(gaussian_bandit, key, 2)[1])"
            ]
        },
        {
            "attachments": {},
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## 実装したアルゴリズム\n",
                "\n",
                "* ``pull_all``：各アームを順番に引くアルゴリズム\n",
                "* ``ETC``：[ETCアルゴリズム](BANDIT_basics.ipynb)\n",
                "* ``UCB``：[UCBアルゴリズム](BANDIT_basics.ipynb)\n",
                "    * 証明：[BANDIT_UCB_regret_proof.ipynb](BANDIT_UCB_regret_proof.ipynb)\n",
                "* ``EXP3``：[EXP3アルゴリズム](BANDIT_adversarial.ipynb)\n",
                "    * 証明：[BANDIT_adversarial.ipynb](BANDIT_adversarial.ipynb)\n",
                "    * 敵対的バンディットですが、同じ確率的バンディットにも実装できます。"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 2,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "text/plain": [
                            "<matplotlib.legend.Legend at 0x7f24c81083a0>"
                        ]
                    },
                    "execution_count": 2,
                    "metadata": {},
                    "output_type": "execute_result"
                },
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGdCAYAAAAMm0nCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfZklEQVR4nO3dd3hUVf7H8feUzKRXSIPQIfQOIfbCEpSfFV1FVFBsLKiIqy66ltVV7A1RdhVFXUBx7eCqCAiidAjdAFJCSUIJyaRnMnN+f7CMZkElkGRSPq/nmWeZe8/c+72HZebjPefeazHGGERERETqEau/CxARERGpKgUYERERqXcUYERERKTeUYARERGRekcBRkREROodBRgRERGpdxRgREREpN5RgBEREZF6x+7vAmqK1+tl3759hIWFYbFY/F2OiIiInABjDAUFBSQmJmK1/vp5lgYbYPbt20dSUpK/yxAREZGTsHv3bpo3b/6r6xtsgAkLCwOOdEB4eLifqxEREZET4XK5SEpK8v2O/5oGG2CODhuFh4crwIiIiNQzvzf9Q5N4RUREpN5RgBEREZF6RwFGRERE6p0GOwfmRBhjqKiowOPx+LuURsNms2G323Vpu4iInJJGG2DKy8vJysqiuLjY36U0OsHBwSQkJOBwOPxdioiI1FONMsB4vV527NiBzWYjMTERh8OhMwK1wBhDeXk5Bw4cYMeOHbRv3/43b1IkIiLyaxplgCkvL8fr9ZKUlERwcLC/y2lUgoKCCAgIYNeuXZSXlxMYGOjvkkREpB5q1P/5q//69w/1u4iInCr9koiIiEi9owAjIiIi9Y4CTCPTqlUrXnzxRd97i8XCJ598ckrb/Pbbb7FYLOTl5QEwbdo0IiMjT2mbIiIiv6VKAWbixIn069ePsLAwYmNjufTSS8nIyKjU5pxzzsFisVR63XbbbZXaZGZmMmTIEIKDg4mNjeWee+6hoqKiUptvv/2W3r1743Q6adeuHdOmTTu5I2xARo4c6etTh8NBu3btePTRR4/pu5qwZMkSbDYbQ4YMqfF9iYiI/J4qBZiFCxcyZswYli5dyty5c3G73QwaNIiioqJK7W6++WaysrJ8r6efftq3zuPxMGTIEMrLy/nhhx94++23mTZtGg899JCvzY4dOxgyZAjnnnsu6enpjBs3jptuuomvvvrqFA+3/hs8eDBZWVls3bqVu+++m0ceeYRnnnmmxvc7depUbr/9dhYtWsS+fftqfH8iIlJ3Ldt+iOumLqO4vOb/A/rXVOky6i+//LLS+2nTphEbG8uqVas466yzfMuDg4OJj48/7ja+/vprNm3axDfffENcXBw9e/bkscce47777uORRx7B4XAwZcoUWrduzXPPPQdAp06dWLx4MS+88AJpaWlVPcbfZYyhxO2fu/EGBdiqdA8ap9Pp69vRo0fz8ccf89lnnzFhwgTOOeccevbsWWmI6NJLLyUyMvKUzmAVFhby/vvvs3LlSrKzs5k2bRr333//SW9PRETqJ4/X8OqCbbzwzRa8BiYv2MY9aR39Ussp3QcmPz8fgOjo6ErLp0+fzr/+9S/i4+O56KKLePDBB333W1myZAndunUjLi7O1z4tLY3Ro0ezceNGevXqxZIlSxg4cGClbaalpTFu3LhfraWsrIyysjLfe5fLdcLHUeL20Pkh/5zd2fRoGsGOk/9rCAoK4tChQ9VY0bFmzZpFx44dSU5O5tprr2XcuHFMmDBBN/8TEWlEclyljJ+VzvfbjvzmXN67GX86p53f6jnpX06v18u4ceM4/fTT6dq1q2/5NddcQ8uWLUlMTGTdunXcd999ZGRk8NFHHwGQnZ1dKbwAvvfZ2dm/2cblclFSUkJQUNAx9UycOJG//e1vJ3s49Y4xhnnz5vHVV19x++231+i+pk6dyrXXXgscGcLKz89n4cKFnHPOOTW6XxERqRsWbz3Ine+t4VBROUEBNh69pAtX9k3ya00nHWDGjBnDhg0bWLx4caXlt9xyi+/P3bp1IyEhgfPPP5+ffvqJtm3bnnylv2PChAmMHz/e997lcpGUdGKdGxRgY9Oj1T80daL7rorZs2cTGhqK2+3G6/VyzTXX8Mgjj9RMcUBGRgbLly/n448/BsBut3PVVVcxdepUBRgRkQauwuPl5XlbmbRgG8ZAx/gwXrmmN+1iQ/1d2skFmLFjxzJ79mwWLVpE8+bNf7NtSkoKANu2baNt27bEx8ezfPnySm1ycnIAfHM74uPjfct+2SY8PPy4Z1/gyNwQp9N5MoeDxWI5pWGc2nTuuefy2muv4XA4SExMxG7/uW6r1YoxplJ7t9t9SvubOnUqFRUVJCYm+pYZY3A6nbzyyitERESc0vZFRKRuys4v5faZq1mx8zAAw/q34OGLOhNYxf/wrilVugrJGMPYsWP5+OOPmT9/Pq1bt/7dz6SnpwOQkJAAQGpqKuvXr2f//v2+NnPnziU8PJzOnTv72sybN6/SdubOnUtqampVym2QQkJCaNeuHS1atKgUXgCaNm1KVlaW773H42HDhg0nva+KigreeecdnnvuOdLT032vtWvXkpiYyMyZM0962yIiUnct3HKAC1/+jhU7DxPqtPPS1T2ZeHm3OhNeoIpnYMaMGcOMGTP49NNPCQsL881ZiYiIICgoiJ9++okZM2Zw4YUXEhMTw7p167jrrrs466yz6N69OwCDBg2ic+fOXHfddTz99NNkZ2fz17/+lTFjxvjOoNx222288sor3Hvvvdx4443Mnz+fWbNmMWfOnGo+/IblvPPOY/z48cyZM4e2bdvy/PPP+24udzJmz57N4cOHGTVq1DFnWoYOHcrUqVOPucePiIjUX26Pl2e/zuAfC7cD0DkhnFeH96ZVkxA/V3asKp2Bee2118jPz+ecc84hISHB93r//fcBcDgcfPPNNwwaNIiOHTty9913M3ToUD7//HPfNmw2G7Nnz8Zms5Gamsq1117L9ddfz6OPPupr07p1a+bMmcPcuXPp0aMHzz33HG+88UaNXELdkNx4442MGDGC66+/nrPPPps2bdpw7rnnnvT2pk6dysCBA487TDR06FBWrlzJunXrTqVkERGpI/bllXD1P5f6wst1A1ry0Z9Oq5PhBcBi/nfSRAPhcrmIiIggPz+f8PDwSutKS0vZsWMHrVu3JjAw0E8VNl7qfxGRuuXrjdnc++E68ordhAXaeeaK7gzumuCXWn7r9/uX6sfMVREREal2bo+Xp/7zI28s3gFAt2YRTL6mNy1igv1c2e9TgBEREWmEMg8Vc8d7a0jfnQfATWe05t7BHXHY68dznhVgREREGpkvN2Rxz7/XUVBaQXignWev7MGgLsd/BFBdpQAjIiLSSJS6PTz5nx+Z9sNOAHq3iOTlYb1oHlX3h4z+lwKMiIhII7DzYBFjZqxm474jzwq89aw2/DktmQBb/Rgy+l8KMCIiIg3cnHVZ3PfhOgrLKogOcfDclT04t2Osv8s6JQowIiIiDVSp28PfPt/IzOW7AejXKoqXh/UiIeL4j+WpTxRgREREGqDtBwoZM2MNm7NcWCww+uy2jP9DB+z1dMjofynAiIiINDAfrd7DXz/ZQHG5h5gQBy9e3ZMz2zf1d1nVqmHEsEZk5MiRWCyWY16BgYHHXf7L17fffgvAhx9+yDnnnENERAShoaF0796dRx99lNzcXP8enIiInJKScg/3/nst42etpbjcQ0rraL6488wGF15AZ2DqpcGDB/PWW29VWmaxWPjlUyHuvPNOXC5XpXbR0dE88MADPPXUU9x111088cQTJCYmsnXrVqZMmcK7777LnXfeWWvHISIi1Scju4DbZ65mS04hFguMO78DY89rh81q8XdpNUIBph5yOp3Ex//2DYeCgoIoKyur1G758uU88cQTvPjii5WCSqtWrfjDH/5wSk+uFhER/5m1cjcPfbqBUreXJqFOXr66J6e1a+LvsmqUAgyAMeAu9s++A4LBUjvpePr06YSGhvKnP/3puOsjIyNrpQ4REakehWUVPPzpRj5cvQeAM9s34YWretIk1OnnymqeAgwcCS9PJPpn3/fvA0fVHlU+e/ZsQkNDK2/m/vu5//77f/NzW7dupU2bNgQEBFS5TBERqVs2Z7kYM3012w8WYbXA+D904E/ntMPaQIeM/pcCTD107rnn8tprr1VaFh0d/buf++UcGRERqZ+MMby/YjcPf7aRsgovCRGBvHhVT1LaxPi7tFqlAANHhnHu3+e/fVdRSEgI7dq1q/LnOnTowOLFi3G73ToLIyJSDxWUuvnLR+uZsy4LgHOSm/L8H3sSHeLwc2W1T5dRw5E5KI4Q/7xqaf4LwDXXXENhYSGvvvrqcddrEq+ISN21cV8+F7/yPXPWZWG3WvjLBR15c0S/RhleQGdg6qWysjKys7MrLbPb7TRp8tszzlNSUrj33nu5++672bt3L5dddhmJiYls27aNKVOmcMYZZ+gyahGROsYYwztLdvH4nM2Ue7wkRgQy6Zre9GkZ5e/S/EoBph768ssvSUhIqLQsOTmZH3/88Xc/+9RTT9GnTx8mT57MlClT8Hq9tG3bliuuuIIRI0bUVMkiInISXKVu/vLhOr5Yf+Q/Wgd2iuWZK3oQ1UjPuvySxTTQmZ0ul4uIiAjy8/MJDw+vtK60tJQdO3bQunVrAgMD/VRh46X+FxH5fem787hj5hoyc4sJsFmYcEEnbji9FZZanHrgD7/1+/1LOgMjIiJShxhjePP7nTz5n824PYZmkUFMHt6bnkmR/i6tTlGAERERqSMOF5Vz34fr+HpTDgAXdotn4uXdiQjSlaP/SwFGRESkDlideZix01ezL7+UAJuFvw7pzPWpLRv8kNHJUoARERHxI6/XMHXxDp768kcqvIZWMcFMGtabbs0j/F1anaYAIyIi4ie5ReWMn5XOtxkHABjSPYEnL+9GWKCGjH6PAoyIiIgfrNyZy+0z15CVX4rTbuXB/+vM8JQWGjI6QQowIiIitcjrNbz67TZe+GYrHq+hTZMQJg/vTaeEX79kWI6lACMiIlJLDhWWMe79dL7behCAS3om8vhl3Qh16ue4qtRjIiIiteCHbQcZ9346+wvKCAyw8tglXbmyb5K/y6q3FGBERERqkMdrmDR/Ky/N24ox0LZpCK8O70NyfJi/S6vX9DTqeuScc85h3LhxxyyfNm0akZGRvvcul4sHHniAjh07EhgYSHx8PAMHDuSjjz7i6JMjzjnnHCwWi+8VFxfHlVdeya5du2rpaEREGr4cVynXv7mMF785El6u6pvE7NvPVHipBgowDUxeXh6nnXYa77zzDhMmTGD16tUsWrSIq666invvvZf8/Hxf25tvvpmsrCz27dvHp59+yu7du7n22mv9WL2ISMPx3dYDXPjSd3y/7RBBATaeu7IHT13RnSCHzd+lNQgaQmpg7r//fnbu3MmWLVtITEz0Le/QoQPDhg2r9PDE4OBg4uPjAUhISGDs2LHceuuttV6ziEhDUuHx8vK8rUxasA1joFNCOJOG9aJdbKi/S2tQFGA48uCskooSv+w7yB5Ubdf8e71e3nvvPYYPH14pvBwVGvrr/3hyc3OZNWsWKSkp1VKLiEhjlJ1fyu0zV7Ni52EAhvVvwcMXdSYwQGddqpsCDFBSUULKDP/8cC+7ZhnBAcHVsq2DBw9y+PBhOnbseELtX331Vd544w2MMRQXF9OhQwe++uqraqlFRKSxWbjlAHe9n05uUTmhTjuPX9aVS3o283dZDZbmwDQgRyfonqjhw4eTnp7O2rVrWbx4Me3atWPQoEEUFBTUUIUiIg2P2+Nl4hebGfHmcnKLyumSGM7s289QeKlhOgPDkWGcZdcs89u+T1R4eHilSbhH5eXlERERQdOmTYmMjOTHH388oe1FRETQrl07ANq1a8fUqVNJSEjg/fff56abbjrhukREGqt9eSXcPnMNq3YdGTK6bkBLHhjSSUNGtUABBrBYLNU2jFOTkpOT+frrr49Zvnr1ajp06IDVauXqq6/m3Xff5eGHHz5mHkxhYSGBgYHY7cf/a7fZjvyDKynxz3wgEZH65KuN2dz34Tryit2EBdp55oruDO6a4O+yGg0NIdUjo0ePZsuWLdxxxx2sW7eOjIwMnn/+eWbOnMndd98NwOOPP05SUhIpKSm88847bNq0ia1bt/Lmm2/Sq1cvCgsLfdsrLi4mOzub7Oxs1q5dy+jRowkMDGTQoEH+OkQRkTqvvMLLY7M3ceu7q8grdtOtWQRzbj9T4aWW6QxMPdKmTRsWLVrEAw88wMCBAykvL6djx4588MEHDB48GIDo6GiWLl3Kk08+yd///nd27dpFVFQU3bp145lnniEiIsK3vddff53XX38dgKioKLp3784XX3xBcnKyX45PRKSuyzxUzJ3vr2FNZh4AN53RmnsHd8Rh1/mA2mYxVZ35WU+4XC4iIiLIz88nPLzyEz5LS0vZsWMHrVu3rnRfFKkd6n8RqY++3JDFPf9eR0FpBeGBdp69sgeDusT7u6wG57d+v39JZ2BERER+Q6nbw5P/+ZFpP+wEoE/LKF66uifNo+r+3MmGTAFGRETkV+w8WMSYGavZuM8FwK1nteHPackE2DRk5G8KMCIiIscxZ10W9324jsKyCqJDHDx7ZXfO6xjn77LkvxRgREREfqHU7eFvn29k5vLdAPRrFcXLw3qREHHi9+2SmqcAIyIi8l/bDxQyZsYaNme5sFhg9NltGf+HDtg1ZFTnKMCIiIgAH63ew4OfbKCo3ENMiIMXr+7Jme2b+rss+RUKMCIi0qiVlHt4+LMNzFq5B4D+raOZNKwXceG6zUNdpgAjIiKN1pacAm6fsYaMnAIsFhh3fgfGntcOm9Xi79LkdyjAiIhIo/TByt08+OkGSt1emoQ6efnqnpzWrom/y5ITpAAjIiKNSmFZBQ99uoGPVu8F4Mz2TXj+jz1pGub0c2VSFZpWXc+MHDkSi8VyzGvw4MHs27ePqKgoXn755UqfWbZsGQEBAb4nWU+bNs33OavVSvPmzbnhhhvYv3+/7zMXX3wxLVq0IDAwkISEBK677jr27dtXq8cqIlLdNme5uHjSYj5avRerBf48qANv39Bf4aUe0hmYemjw4MG89dZblZY5nU6ioqKYNGkSt956KxdccAHt27enpKSEESNGcNNNN1V6ynR4eDgZGRl4vV7Wrl3LDTfcwL59+/jqq68AOPfcc7n//vtJSEhg7969/PnPf+aKK67ghx9+qNVjFRGpDsYY3l+xm4c/20hZhZeEiEBevKonKW1i/F2anCQFmHrI6XQSH3/8B4hde+21fPTRR4wcOZLvvvuOCRMm4Ha7eeaZZyq1s1gsvm0kJiZyxx138OCDD1JSUkJQUBB33XWXr23Lli35y1/+wqWXXorb7SYgIKDmDk5EpJoVlLqZ8NF6Zq/LAuCc5KY8/8eeRIc4/FyZnAoFGI4kc1NS4pd9W4KCsFiqd7b7lClT6Nq1K8OHD+eDDz5g/vz5hIaG/uZngoKC8Hq9VFRUHLMuNzeX6dOnc9pppym8iEi9snFfPmNnrGHHwSLsVgt3D0rm1rPaYNVVRvWeAgxgSkrI6N3HL/tOXr0KS3DVnmg6e/bsYwLJ/fffz/333w9AbGwsjz32GLfddhujR4/mrLPO+s3tbd26lSlTptC3b1/CwsJ8y++77z5eeeUViouLGTBgALNnz65SnSIi/mKM4Z0lu3h8zmbKPV4SIwKZdE1v+rSM8ndpUk2qNIl34sSJ9OvXj7CwMGJjY7n00kvJyMio1Ka0tJQxY8YQExNDaGgoQ4cOJScnp1KbzMxMhgwZQnBwMLGxsdxzzz3H/Jf/t99+S+/evXE6nbRr145p06ad3BE2QOeeey7p6emVXrfddptvvcfjYdq0aQQHB7N06dLjnlXJz88nNDSU4OBgkpOTiYuLY/r06ZXa3HPPPaxZs4avv/4am83G9ddfjzGmxo9PRORUuErdjJmxmoc/20i5x8vATnF8ceeZCi8NTJXOwCxcuJAxY8bQr18/KioquP/++xk0aBCbNm0iJCQEgLvuuos5c+bwwQcfEBERwdixY7n88sv5/vvvgSM/rkOGDCE+Pp4ffviBrKwsrr/+egICAnjiiScA2LFjB0OGDOG2225j+vTpzJs3j5tuuomEhATS0tKquQuODOMkr15V7ds90X1XVUhICO3atfvV9c8++yzbt29n5cqVnH322TzxxBM89NBDldqEhYWxevVqrFYrCQkJBB2njiZNmtCkSRM6dOhAp06dSEpKYunSpaSmpla5ZhGR2rAm8zB3vpdOZm4xATYLEy7oxA2nt6r2oXqpA8wp2L9/vwHMwoULjTHG5OXlmYCAAPPBBx/42mzevNkAZsmSJcYYY7744gtjtVpNdna2r81rr71mwsPDTVlZmTHGmHvvvdd06dKl0r6uuuoqk5aWdsK15efnG8Dk5+cfs66kpMRs2rTJlJSUnPjB1hEjRowwl1xyya+u37Bhg3E6nebDDz80xhjz73//2wQEBJi1a9f62rz11lsmIiKiSvvdtWuXAcyCBQtOourK6nP/i0jd5PV6zRvfbTft7p9jWt4325w2cZ5Zk3nY32XJSfit3+9fOqX7wOTn5wMQHR0NwKpVq3C73QwcONDXpmPHjrRo0YIlS5YAsGTJErp160ZcXJyvTVpaGi6Xi40bN/ra/HIbR9sc3cbxlJWV4XK5Kr0aqrKyMrKzsyu9Dh48SEVFBSNGjODyyy/n8ssvB2Do0KEMHTqUkSNHHnco6XiWLVvGK6+8Qnp6Ort27WL+/PkMGzaMtm3b6uyLiNQ5h4vKufXdVTw2exNuj+HCbvF8ceeZ9EyK9HdpUoNOOsB4vV7GjRvH6aefTteuXQHIzs7G4XAQGRlZqW1cXBzZ2dm+Nr8ML0fXH133W21cLhclv3K10MSJE4mIiPC9kpKSTvbQ6rwvv/yShISESq8zzjiDJ554gr179/LKK69Uaj958mSysrJ8Q3S/Jzg4mI8++ojzzz+f5ORkRo0aRffu3Vm4cCFOp272JCJ1x+rMwwx5+Tu+3pRDgM3C3y7uwuRrehMRpCsmG7qTvgppzJgxbNiwgcWLF1dnPSdtwoQJjB8/3vfe5XI1yBAzbdq035zQ/L9zXeDIGbKsrCzf+5EjRzJy5Mhf3Ua3bt2YP3/+qZQpIlKjvF7D1MU7eOrLH6nwGlo3CWHSsF50bRbh79KklpxUgBk7diyzZ89m0aJFNG/e3Lc8Pj6e8vJy8vLyKp2FycnJ8d00LT4+nuXLl1fa3tGrlH7Z5n+vXMrJySE8PPy4k03hyM3ddHZARKThyy0q5+5Z6SzIOADA/3VPYOLl3QgL1FmXxqRKQ0jGGMaOHcvHH3/M/Pnzad26daX1ffr0ISAggHnz5vmWZWRkkJmZ6Zs7kZqayvr16ys9d2fu3LmEh4fTuXNnX5tfbuNoG82/EBFp3FbuzGXIy9+xIOMATruVv1/alUnDeim8NEJVOgMzZswYZsyYwaeffkpYWJhvzkpERARBQUFEREQwatQoxo8fT3R0NOHh4dx+++2kpqYyYMAAAAYNGkTnzp257rrrePrpp8nOzuavf/0rY8aM8Z1Bue2223jllVe49957ufHGG5k/fz6zZs1izpw51Xz4IiJSH3i9hskLtvHCN1vwGmjTJITJw3vTKSHc36WJn1QpwLz22msAnHPOOZWWv/XWW745FS+88AJWq5WhQ4dSVlZGWloar776qq+tzWZj9uzZjB49mtTUVEJCQhgxYgSPPvqor03r1q2ZM2cOd911Fy+99BLNmzfnjTfeqJF7wIiISN12qLCMce+n893WgwBc1qsZj13alVCnbibfmFmMaZi3VnW5XERERJCfn094eOWEXlpayo4dO2jdujWBgYF+qrDxUv+LyIn6YdtB7nw/nQMFZQQGWHnskq5c2bfhXaAhP/ut3+9fatTxtYFmtzpP/S4iv8fjNUyav5WX5m3FGGgXG8rka3qTHB/2+x+WRqFRBpijT1QuLi7+1auapOYUFxcD6MnWInJc+12l3PleOku2HwLgqr5JPHJxF4IcNj9XJnVJowwwNpuNyMhI35VQwcHBek5GLTDGUFxczP79+4mMjMRm05eRiFS2eOtBxr2/hoOF5QQF2Pj7pV0Z2qf5739QGp1GGWDg53vO/PJybqkdkZGRvv4XEQGo8Hh58ZutTP52G8ZAx/gwXh3emzZNQ/1dmtRRjTbAWCwWEhISiI2Nxe12+7ucRiMgIEBnXkSkkuz8Um6fuZoVOw8DMKx/Eg9f1IXAAH1XyK9rtAHmKJvNph9UERE/WZCxn7tnrSW3qJxQp52Jl3fjoh6J/i5L6oFGH2BERKT2uT1env0qg38s2g5A54RwXh3em1ZNQvxcmdQXCjAiIlKr9uaVcMfMNazadWTIaERqSyZc2ElDRlIlCjAiIlJrvt6YzZ8/WIurtIIwp52nr+jOBd0S/F2W1EMKMCIiUuPcHi9Pf/kjr3+3A4AeSZG8fHVPWsZoyEhOjgKMiIjUqMxDxYyduZp1e/IBGHVGa+4b3BGH3ernyqQ+U4AREZEa8+WGbO7591oKSiuICArg6Su6k9ZF94GSU6cAIyIi1a7U7eGJLzbzzpJdAPRuEcmka3rTLFKPb5HqoQAjIiLVatehIv40fTUb97kAuPnM1tw7uCMBNg0ZSfVRgBERkWrz+dp93P/RegrKKogKDuCFq3pyTnKsv8uSBkgBRkRETlmp28OjszcxY1kmAH1bRjHpml4kRGjISGqGAoyIiJySnw4UMnbGGjZnubBYYOy57bjj/PYaMpIapQAjIiIn7aPVe3jg4w2UuD3EhDh48eqenNm+qb/LkkZAAUZERKqspNzDI59t5P2VuwE4rW0Mz/+xJ/ERgX6uTBoLBRgREamSLTkF/Gn6arbtL8RigTvPb8/t57XHZrX4uzRpRBRgRETkhP171R4e/OTIkFFsmJMXrurJ6e2a+LssaYQUYERE5HcVlVXwwMfr+SR9HwBntm/CC1f1pEmo08+VSWOlACMiIr/px2wXf5q+mu0HirBZLdx5fnvGnNtOQ0biVwowIiJyXMYYZizP5LHZmyh1e4kPD2TSNb3o1yra36WJKMCIiMixCkrd3P/xBj5fe2TI6OwOTXn+jz2I0ZCR1BEKMCIiUsn6Pfnc8d4adhw8MmR03+BkbjqjDVYNGUkdogAjIiLAkSGjd5fu4u+zN1Pu8ZIYcWTIqE9LDRlJ3aMAIyIi5Be7uf/j9cxZnwXAoM5xPH1FdyKDHX6uTOT4FGBERBq5tbvzGDNjNXsOl2C3WphwYSduPL0VFouGjKTuUoAREWmkjDFM+2EnT3yxGbfH0CI6mJeH9aJnUqS/SxP5XQowIiKNUF5xOX/+YC3fbN4PwAVd43nqiu6EBwb4uTKRE6MAIyLSyKzJPMzYGWvYm1eCw27lgQs7cX1qSw0ZSb2iACMi0kh4vYZ/LNrOs19n4PEaWsYEM/ma3nRtFuHv0kSqTAFGRKQROFxUzvhZ6SzIOADA/3VPYOLl3QjTkJHUUwowIiIN3LLth7jzvXSyXaU47Fb+dnEXru6XpCEjqdcUYEREGiiv1/Dawp947usMvAZaNwnhlWt60SVRQ0ZS/ynAiIg0QAcKyhg/K53vth4E4PLezfj7pV0JduhrXxoG/T9ZRKSBWfLTIe54bw0HCspw2q08dklXruzbXENG0qAowIiINBAer2Hygm28+M0WvAY6xIUy+ZretI8L83dpItVOAUZEpAHY7yrljvfWsHR7LgBX9mnO3y7poiEjabD0/2wRkXru+20HufO9dA4WlhHisPHYpV25vHdzf5clUqMUYERE6qkKj5cXvtnC5AU/AdAxPozJw3vTtmmonysTqXkKMCIi9VB2fil3zFzD8p1HhoyG9U/iof/rQpDD5ufKRGqHAoyISD2z4Mf9jJ+VzuFiN6FOOxMv78ZFPRL9XZZIrVKAERGpJ9weL899vYUpC48MGXVOCGfy8N60bhLi58pEap8CjIhIPbDncDF3zFzD6sw8AEaktuT+IZ1w2jVkJI2TAoyISB33zaYc7v5gLfklbsKcdp66ojsXdkvwd1kifqUAIyJSR5VVeHj6ywymLt4BQI+kSF4Z1ouk6GA/VybifwowIiJ10O7cYsbOWM3aPfkA3Hh6a/5yQUccdqufKxOpGxRgRETqmP+sz+LeD9dRUFpBRFAAz17Zgz90jvN3WSJ1igKMiEgdUVbh4fE5m3lnyS4AerWI5JVretMsMsjPlYnUPQowIiJ1wK5DRYydsYb1e48MGd16dhv+PCiZAJuGjESORwFGRMTPPlu7jwkfrqOo3ENUcADPX9WTc5Nj/V2WSJ2mACMi4ielbg9/n7OJfy3NBKB/q2hevLoniRoyEvldCjAiIn6wbX8hY6avJiOnAIsFxpzTjnED22PXkJHICVGAERGpZZ+s2cv9H6+nuNxDk1AHz/2xJ2d3aOrvskTqlSpH/UWLFnHRRReRmJiIxWLhk08+qbR+5MiRWCyWSq/BgwdXapObm8vw4cMJDw8nMjKSUaNGUVhYWKnNunXrOPPMMwkMDCQpKYmnn3666kcnIlKHFJdXcM8Haxn3fjrF5R5S28TwxZ1nKryInIQqn4EpKiqiR48e3HjjjVx++eXHbTN48GDeeust33un01lp/fDhw8nKymLu3Lm43W5uuOEGbrnlFmbMmAGAy+Vi0KBBDBw4kClTprB+/XpuvPFGIiMjueWWW6pasoiI323NKeBP01ezdX8hFgvcfm477hzYAZvV4u/SROqlKgeYCy64gAsuuOA32zidTuLj44+7bvPmzXz55ZesWLGCvn37AjBp0iQuvPBCnn32WRITE5k+fTrl5eW8+eabOBwOunTpQnp6Os8//7wCjIjUK8YYPli5h4c/20iJ20NsmJOXh/ViQJsYf5cmUq/VyGyxb7/9ltjYWJKTkxk9ejSHDh3yrVuyZAmRkZG+8AIwcOBArFYry5Yt87U566yzcDgcvjZpaWlkZGRw+PDh4+6zrKwMl8tV6SUi4k9FZRWMn7WWez9cR4nbwxntmjDnjjMVXkSqQbUHmMGDB/POO+8wb948nnrqKRYuXMgFF1yAx+MBIDs7m9jYyvc3sNvtREdHk52d7WsTF1f5ttlH3x9t878mTpxIRESE75WUlFTdhyYicsI27XNx8SuL+XjNXmxWC/ekJfP2jf1pGub8/Q+LyO+q9quQrr76at+fu3XrRvfu3Wnbti3ffvst559/fnXvzmfChAmMHz/e997lcinEiEitM8YwY3kmf/t8E+UVXuLCnUwa1pv+raP9XZpIg1Ljl1G3adOGJk2asG3bNs4//3zi4+PZv39/pTYVFRXk5ub65s3Ex8eTk5NTqc3R9782t8bpdB4zWVhEpDYVllVw/0fr+WztPgDO6xjLM1d0JyZU300i1a3G75i0Z88eDh06REJCAgCpqank5eWxatUqX5v58+fj9XpJSUnxtVm0aBFut9vXZu7cuSQnJxMVFVXTJYuIVNmGvfkMefk7Plu7D5vVwoQLOvLG9X0VXkRqSJUDTGFhIenp6aSnpwOwY8cO0tPTyczMpLCwkHvuuYelS5eyc+dO5s2bxyWXXEK7du1IS0sDoFOnTgwePJibb76Z5cuX8/333zN27FiuvvpqEhMTAbjmmmtwOByMGjWKjRs38v777/PSSy9VGiISEakLjDG8u2Qnl7/6A7sOFdMsMoj3bxnArWe3xapLpEVqjMUYY6rygW+//ZZzzz33mOUjRozgtdde49JLL2XNmjXk5eWRmJjIoEGDeOyxxypNys3NzWXs2LF8/vnnWK1Whg4dyssvv0xoaKivzbp16xgzZgwrVqygSZMm3H777dx3330nXKfL5SIiIoL8/HzCw8OrcogiIickv8TNff9ex5cbj1xcMLBTHM9d2YOI4AA/VyZSf53o73eVA0x9oQAjIjVp7e48xsxYzZ7DJQTYLNyb1pGbzmyNxaKzLiKn4kR/v/UsJBGRKjDGMHXxDp768kfcHkNSdBCTr+lN9+aR/i5NpFFRgBEROUH5xW7u/mAt32w+clVkWpc4nr6iBxFBGjISqW0KMCIiJ2DVrlzumJnO3rwSHDYrD/5fJ64d0FJDRiJ+ogAjIvIbvF7D699t5+mvMvB4DS2ig3l1eG+6Novwd2kijZoCjIjIrzhUWMa9/17HvB+P3Hzz4h6J/P2yroQHashIxN8UYEREjmP5jlxun7maHFcZDruVRy7qwrD+SRoyEqkjFGBERH7B4zX8Y9FPPPf1FjxeQ9umIUwa1pvOibodg0hdogAjIvJfBwrKuOv9dBZvOwjAZb2a8fdLuxLi1FelSF2jf5UiIsDS7Ye4Y+Ya9heUERRg428Xd+HKvs01ZCRSRynAiEij5vEaXp63lZfnb8UYaB8byqvDe9M+LszfpYnIb1CAEZFGa39BKePeS+eHnw4BcEWf5jx6SReCHfpqFKnr9K9URBqlRVsOMH5WOgcLywl22Hjism5c2quZv8sSkROkACMijUqFx8uL32xl8rfbMAaS48KYPLw37WJD/V2aiFSBAoyINBr78koY9346y3fkAjA8pQUP/l9nAgNsfq5MRKpKAUZEGoUFGfsZ/346h4vdhDhsPHF5Ny7pqSEjkfpKAUZEGjS3x8tzX29hysKfAOjaLJxXhvWmVZMQP1cmIqdCAUZEGqy9eSWMnbGaNZl5AFyf2pIHhnTCadeQkUh9pwAjIg3SN5tyuPuDteSXuAkLtPP00O5c0C3B32WJSDVRgBGRBqW8wsuT//mRN7/fAUD35hFMvqY3SdHBfq5MRKqTAoyINBi7c4sZO2M1a/fkAzDqjNbcOzhZQ0YiDZACjIg0CF+sz+K+D9dRUFpBRFAAz13Zg4Gd4/xdlojUEAUYEanXyio8PDFnM28v2QVArxaRTBrWi+ZRGjISacgUYESk3tpxsIjbZ65mw14XALed3ZY/D+qA3Wb1c2UiUtMUYESkXvp87T7+8uE6iso9RAYH8Pwfe3BeRw0ZiTQWCjAiUq+UlhQx790nOJC5jbuBuOhAzkluSsiO+bDD39WJNDI9robEnn7ZtQKMiNQbmVvXUf7eCIZ4tv/87VUMrPFnVSKNWPO+CjAiIr9l5Wev0WnVI4RYSjlMOIWdryYpJszfZYk0bk07+m3XCjAiUqeVFBWw4fWb6Zf3H7DAJkc3Ykf+i6TEVv4uTUT8SAFGROqsXZtX4f3gBvp5d+ExFpa3vJn+1z+BzR7g79JExM8UYESkzjFeLys+fpnu6/5OoMXNQSLJHvQKqadf5O/SRKSOUIARkTqluDCfja/fTP/8r8AC6529SbjhbbrGt/B3aSJShyjAiEidsX3DMgI+HEk/s+/IkFGbP5Fy7WNYbXqWkYhUpgAjIn53ZMjoJbqve5xAi5v9RHMgbTKpp13o79JEpI5SgBERvyrIz2XLGzfSv2ABWGBtUH9ajnqXLk3i/V2aiNRhCjAi4jc/rfsB58c30sdkUWGsrGg7lpThj2jISER+lwKMiNQ64/WyfNZT9Nz8HE6Lm2yaknfRP0jte76/SxORekIBRkRqVf7hg/z0xkhSir4DC6QHp9J61NvEx+hBjCJy4hRgRKTWbFn9LSGf30Jvk0O5sbE6eTwpV9+PxWr1d2kiUs8owIhIjTNeL8tm/p0+W14kwOJhnyWWwotfZ0Dvc/xdmojUUwowIlKj8g/lsH3qDQwo/h4ssDr0LNre+CaJ0U39XZqI1GMKMCJSY35cPpeoL26lF4coN3bWdL6H/lfeqyEjETllCjAiUu28Hg/LZz5Kn62TCLB42G1JpPTSf5LS80x/lyYiDYQCjIhUq9z9e9nz5ggGlK4AC6wKO4/km98kNDzK36WJSAOiACMi1Wbzsq+I+c9tdCeXUhPA2i730f+KuzVkJCLVTgFGRE6Z1+Nh2bsP0G/HP7BbvGRam+EZ+hYpXVL8XZqINFAKMCJySg7l7GHfm9eRWrYaLLAyfCCdbp5KSFikv0sTkQZMAUZETtqGxZ8R+82ddCOXEuNgfc+H6HfJGA0ZiUiNU4ARkSrzVFSw/O2/kJL5BlaLYZc1CXPlNPp36uvv0kSkkVCAEZEqObBvJzlvjyC1LB0ssDxqCF1HvUZwaIS/SxORRkQBRkRO2PqFH9FswZ10xUWxcbKpz9/of/Fof5clIo2QAoyI/K4Kdzkr3/4L/Xe/idVi+MnWGvsfp9E3uae/SxORRkoBRkR+U86en8h9+1oGuDeBBZbFXEKPm14jMCjE36WJSCOmACMiv2rdgn+TtPAuOuGiwASR0e9RUv7vFn+XJSKiACMix3KXl7HqzbsYkD0dgG22tgQNf5e+bbr4uTIRkSMUYESkkuzd28h7+1oGVGwGYFmTy+kx6hUNGYlInaIAIyI+a77+F61/uI+OFOIimG2pT5OSdp2/yxIROUaVb5e5aNEiLrroIhITE7FYLHzyySeV1htjeOihh0hISCAoKIiBAweydevWSm1yc3MZPnw44eHhREZGMmrUKAoLCyu1WbduHWeeeSaBgYEkJSXx9NNPV/3oROSElJeVsvS1W+n1wxgiKWSrrR0F18+nt8KLiNRRVQ4wRUVF9OjRg8mTJx93/dNPP83LL7/MlClTWLZsGSEhIaSlpVFaWuprM3z4cDZu3MjcuXOZPXs2ixYt4pZbfp4Y6HK5GDRoEC1btmTVqlU888wzPPLII/zzn/88iUMUkd+yd/tmdj5zJgNy3gNgadwwWt77Pc3adPJzZSIiv8GcAsB8/PHHvvder9fEx8ebZ555xrcsLy/POJ1OM3PmTGOMMZs2bTKAWbFiha/Nf/7zH2OxWMzevXuNMca8+uqrJioqypSVlfna3HfffSY5OfmEa8vPzzeAyc/PP9nDE2nwVn85zbgeijPm4XCT/3CCWf3Vu/4uSUQauRP9/a7WJ67t2LGD7OxsBg4c6FsWERFBSkoKS5YsAWDJkiVERkbSt+/Pz0wZOHAgVquVZcuW+dqcddZZOBwOX5u0tDQyMjI4fPjwcfddVlaGy+Wq9BKR4ystKWLZ5FH0WnIHYZYSfgzoTPGNC+k16Fp/lyYickKqNcBkZ2cDEBcXV2l5XFycb112djaxsbGV1tvtdqKjoyu1Od42frmP/zVx4kQiIiJ8r6SkpFM/IJEGaM+2Dex99gxSDvwbgKXxw2l7z7fEt2jv58pERE5cg3nm/YQJE8jPz/e9du/e7e+SROqcVXPeIPLdgbT1bOcw4aw963UG3PYqAQ6nv0sTEamSar2MOj4+HoCcnBwSEhJ8y3NycujZs6evzf79+yt9rqKigtzcXN/n4+PjycnJqdTm6Pujbf6X0+nE6dSXsMjxlBYXsu6N2+if+zlYYHNAF6JHvEuP5m39XZqIyEmp1jMwrVu3Jj4+nnnz5vmWuVwuli1bRmpqKgCpqank5eWxatUqX5v58+fj9XpJSUnxtVm0aBFut9vXZu7cuSQnJxMVFVWdJYs0eLu3rmXfs6fTP/dzvMbCkmY30P7eb4lTeBGReqzKAaawsJD09HTS09OBIxN309PTyczMxGKxMG7cOP7+97/z2WefsX79eq6//noSExO59NJLAejUqRODBw/m5ptvZvny5Xz//feMHTuWq6++msTERACuueYaHA4Ho0aNYuPGjbz//vu89NJLjB8/vtoOXKQxWPHJZGL+9QfaeHdyiAg2nv82qTe/iD3A8fsfFhGpy6p6edOCBQsMcMxrxIgRxpgjl1I/+OCDJi4uzjidTnP++eebjIyMSts4dOiQGTZsmAkNDTXh4eHmhhtuMAUFBZXarF271pxxxhnG6XSaZs2amSeffLJKdeoyamnMigtdZtkLVxvzcLgxD4ebDY+fafbv3eHvskREfteJ/n5bjDHGj/mpxrhcLiIiIsjPzyc8PNzf5YjUmh2bVmD59w208u7Gaywsa3kL/a9/AptdTw4RkbrvRH+/9Y0m0oAs/3gS3dIfJchSzkEiyf7DJFLPuNjfZYmIVDsFGJEGoNB1mB/fuJn+rrlggXWBfWl+4zt0jW3m79JERGqEAoxIPbdj4zLs/x5JX7MPj7GwvPVtpFz3OFabzd+liYjUGAUYkXrKeL0s//B5emx4kkCLmxxiyL3wNVJT0vxdmohIjVOAEamHCvJz2fLGjaQULAALrA3sR4tR79KpacLvf1hEpAFQgBGpZ7at/R7nJzfRx+zDbWysan87/Yc9pCEjEWlUFGBE6gnj9bJ81lP02vwsDksF2TQh7//+wYB+A3//wyIiDYwCjEg94Mo7xLY3biClcCFYYE3wabS+8S06Njn+s8FERBo6BRiROm7L6m8J++wmenPgyJBR8l2kXP0AFmuDeZi8iEiVKcCI1FHG62XZe0/QO+N5HBYP+yxxFPzfFAb0Pc/fpYmI+J0CjEgdlH8ohx1TRzCgeAlYYHXIWbS7eRqJkTH+Lk1EpE5QgBGpY35cOY/I2bfSkwOUmQDWdBxPylV/0ZCRiMgvKMCI1BFej4dl0x+h70+TCbB42GNJoPSyqQzocbq/SxMRqXMUYETqgMMHssh883pSS5YfGTIKPZv2N0+jeUS0v0sTEamTFGBE/GzTkv8Q89UYenCIUhPAuu4P0O+yOzVkJCLyGxRgRPzE6/Gw7N0H6L9jCjaLYbclEfcV0+jfNcXfpYmI1HkKMCJ+cDB7N/umjSS1dCVYYEVEGp1v+ichYZH+Lk1EpF5QgBGpZRu+/5z4uWPpTh4lxsH6ng/S/7I7/F2WiEi9ogAjUksq3OWsePev9N/1T2wWww5rSyxXvkn/Tn39XZqISL2jACNSCw7u20XOtOGklq8HCyyPvJBuN/+ToJAwf5cmIlIvKcCI1LD1iz4lcf7tdCGfIhPIpt4P0/+SP/m7LBGRek0BRqSGVLjLWfHWn0nd9zYA262tsF/9Nv069PRvYSIiDYACjEgN2L93BwenXUuqewMAy6IvpsdNrxEYHOrnykREGgYFGJFqtnb+e7RY9Gc6U0ChCSIj5XFSLhzl77JERBoUBRiRauIuL2PlW+NJzfoXAD/Z2uAY9g592nXzc2UiIg2PAoxINcjalUH+v0aS6t4EwLKmV9Bz1CScgcF+rkxEpGFSgBE5RelzZ9Dm+z+TQBEFJoitAyaScsEN/i5LRKRBU4AROUnlZaWsfvNOBuS8B8AWewdCr/0XvVsl+7kyEZGGTwFG5CTs2/EjhdOvY0DFFgCWxl5F71Ev43AG+rkyEZHGQQFGpIrWfPU2bZdMIJEi8ghl5xnPMmDgMH+XJSLSqCjAiJygstJi0t8YQ8rBjwD40d6JyOvfpWeL9n6uTESk8VGAETkBe7dvpGT6daR4fgJgScK19L3heQIcTj9XJiLSOCnAiPyOlXNep9Pyv9LMUsphwsg8+wVSz73S32WJiDRqCjAiv6K0pIi1b/yJlEOfgAU2BXQl+vq36ZHUzt+liYg0egowIsexKyMd7/sjSPHuxGssLGs+kn4jn8Ye4PB3aSIiggKMyDFWfjaFzqseIthSxiEi2HvuC6SeM9TfZYmIyC8owIj8V3FhPhveuI3+eV+ABTY6ehB347/oHt/C36WJiMj/UIARAXZtXoX5YAT9vbuPDBm1GEX/EU9hs+ufiIhIXaRvZ2nUjNfLyk8m0WXt4wRbyjhAFPvTXiX1tAv9XZqIiPwGBRhptIoK8tj8+k30c80FC6x39iLhhnfpEp/k79JEROR3KMBIo/TT+qXYPx5FX+8eKoyVFW1G03/4oxoyEhGpJ/RtLY2K8XpZ/uHz9NzwJE6Lm/1Ec/CC10gdMNjfpYmISBUowEijUZCfS8bUm0lxfQMWSA8aQIsb3qRzbDN/lyYiIlWkACONwra1iwn65Eb6mhwqjJWV7e8g5ZqHsVit/i5NREROggKMNGjG62X5v5+l18ancFgqyKIp+UNeY0D/P/i7NBEROQUKMNJg5R8+yPY3RpBStBgssCb4NNrc9A4J0U39XZqIiJwiBRhpkLauWUTIZzfRy+RQbuys7nAHKcMe1JCRiEgDoQAjDYrxelk28zF6b3kJh8XDPkscRZe8wYBeZ/m7NBERqUYKMNJg5OceYPsb1zOg+IcjQ0YhZ9DmprdJjGri79JERKSaKcBIg/Dj8rlEfDGaXhyg3NhZ0+U++l/xZw0ZiYg0UAowUq95PR6WT3+Evj+9gt3iZa8ljpLL3iSlxxn+Lk1ERGqQAozUW7n797L7rRsYULIMLLAy7HySb3qDZhHR/i5NRERqmAKM1Eubln5Jky9H04NcykwA6V0n0H/oXRoyEhFpJBRgpF7xVFSwfPpD9Nv+GnaLl13W5lRcNpWUbgP8XZqIiNQiBRipNw5m7ybrretILVtzZMgo/A90uvkNQsIi/V2aiIjUsmo/3/7II49gsVgqvTp27OhbX1paypgxY4iJiSE0NJShQ4eSk5NTaRuZmZkMGTKE4OBgYmNjueeee6ioqKjuUqUe2fj9HJhyBt3K1lBsnKzo/ih9xs1SeBERaaRq5AxMly5d+Oabb37eif3n3dx1113MmTOHDz74gIiICMaOHcvll1/O999/D4DH42HIkCHEx8fzww8/kJWVxfXXX09AQABPPPFETZQrdZinooLlb99HSuZUrBbDTmsSlivfpl+nPv4uTURE/KhGAozdbic+Pv6Y5fn5+UydOpUZM2Zw3nnnAfDWW2/RqVMnli5dyoABA/j666/ZtGkT33zzDXFxcfTs2ZPHHnuM++67j0ceeQSHw1ETJUsddDA7k5w3ryW1fC1YYHnkhXS9aQrBoRH+Lk1ERPysRi7Z2Lp1K4mJibRp04bhw4eTmZkJwKpVq3C73QwcONDXtmPHjrRo0YIlS5YAsGTJErp160ZcXJyvTVpaGi6Xi40bN9ZEuVIHrfv2QyxTzqBL+VqKjZOVvZ+i/7iZCi8iIgLUwBmYlJQUpk2bRnJyMllZWfztb3/jzDPPZMOGDWRnZ+NwOIiMjKz0mbi4OLKzswHIzs6uFF6Orj+67teUlZVRVlbme+9yuarpiKQ2VbjLWTHtHlL3TgNgh7Ul1qveoW9yT7/WJSIidUu1B5gLLrjA9+fu3buTkpJCy5YtmTVrFkFBQdW9O5+JEyfyt7/9rca2LzUve/c2ct8dSWr5egCWxVxKj5teJTAoxM+ViYhIXVPjd/2KjIykQ4cObNu2jfj4eMrLy8nLy6vUJicnxzdnJj4+/pirko6+P968mqMmTJhAfn6+77V79+7qPRCpUWvnzyJw6tl0Ll9PoQliVb/nSLn9bYUXERE5rhoPMIWFhfz0008kJCTQp08fAgICmDdvnm99RkYGmZmZpKamApCamsr69evZv3+/r83cuXMJDw+nc+fOv7ofp9NJeHh4pZfUfe7yMpb843Z6LLqZSArZamvH4evm0mfITf4uTURE6rBqH0L685//zEUXXUTLli3Zt28fDz/8MDabjWHDhhEREcGoUaMYP3480dHRhIeHc/vtt5OamsqAAUfupDpo0CA6d+7Mddddx9NPP012djZ//etfGTNmDE6ns7rLFT/K2pWB693rSK3IAGBZk6H0vOkVnIHBfq5MRETqumoPMHv27GHYsGEcOnSIpk2bcsYZZ7B06VKaNm0KwAsvvIDVamXo0KGUlZWRlpbGq6++6vu8zWZj9uzZjB49mtTUVEJCQhgxYgSPPvpodZcqfpT+zUxaL76bBIpwEcK2lMdJueAGf5clIiL1hMUYY/xdRE1wuVxERESQn5+v4aQ6pKy0mDVvjmPA/vcB2GLvQOi1/yKxVbKfKxMRkbrgRH+/9SwkqTX7dvxI0fRrGVCxFYClsX+k96hJOJyBfq5MRETqGwUYqRVrvnqbtksmkEgReYSy88znGHD+1f4uS0RE6ikFGKlRpSVFrJ06lpSDHwHwo70Tkde/S88W7f1cmYiI1GcKMFJj9m7fSMn060jx/ATAkoRr6XvD8wQ4dDWZiIicGgUYqRErZ/+TTisepJmllMOEkXn2C6See6W/yxIRkQZCAUaqVWlJEWvfGE3KoU/BApsDuhA94l16NG/r79JERKQBUYCRarMrI52KWSNJ8ezAaywsS7qBfiOewh7g8HdpIiLSwCjASLVY+dlrdF71MMGWMnIJZ+95L5F69uX+LktERBooBRg5JcWF+WyYOpr+h+eABTY4exI/4h26Jbb0d2kiItKAKcDISdu5eSWWD0bS37v7yJBRy5vpf/1EbHb930pERGqWfmnkpCz/eBLd0h8lyFLOfqLZ/4eXSD3jYn+XJSIijYQCjFRJoeswP75xE/1d34AF1gX2odkN79A1rrm/SxMRkVPkNV5KK0o5UHIAt8eNx3ioMBV4vJ4jf/ZWcLj0MNvytlHhrSCtVRrJ0f55lp0CjJyw7RuWEfDhSPqafVQYKyta3ULK9U9gtdn8XZqISKPg9rr5audX/JT3E1lFWVR4K/AaL8YYvMaLFy8Y8OL1LTeYn9vgpcJbQaG7kLKKMso8ZZR6Sin3lFPuKcdjPFWqp31UewUYqbuM18vyD5+nx4YnCbS42U80uRf+g9SUQf4uTUSkXirzlJFflo/b68btcVPuLWfb4W2s2b8GV7mLCm8FFd4K31mPvLI8DpQcIL8snzJPWY3XF2QPItAWiM1qw2axYbfasVls2Kw2gu3BJEcnE2wPpmW4/y7YUICR31SQn8uWN24kpWDBf4eM+pE06l06Nk3wd2kiIjXGGEOpp5SSihLyy/J9Zyd+GSoOlR5i/YH1HCg+QIU5cibk6FDL0ZfXeDlcepiSihLcXjcV3grcXjdF7iK8xntStcUExnBei/NICkvCaXNitVixYMFisfj+bLVYK723WCxY+XlZmCOMQFsgTrvzyP/anDhsDgKsAThsDoLtwVgslmru1eqlACO/atva73F+chN9zD7cxsaq9nfQf9iDGjISkVphjKHIXUS5t9w3xJFTnOObm+E1Xt/LYzxHhk28Xt/wyS/PYvxvuHCVu0jfn052UfaRZV7Pz9sx3iP785bX6PFZLVYcVgcBtgACrAHEBcfRJ64PCSEJ2K32Sq+QgBDig+MJd4STGJqIzarvYQUYOYbxeln2/pP0/vFZHBYP2TQl76J/MKDv+f4uTUTqKFe5i+yibNweN26vm3JPOfnl+RSWFx4ZJvnv2YcidxH7i/dTUF5AQXmBLzT4zlh4f36fVZRFQXmBvw+NkIAQAm2BlQKFzWLDYXPQvUl3WkW0+nmIxWLDarFit9qxWqzYLDbCHeGEOkKxW+0EWAOwW+2EOcKICYyp82c56jIFGKnElXeIba+PZEDRIrDAmuDTaDNqGvExcf4uTUSqydGzE0fPXhgMHuNh86HNzMucR7mnvFKwOBoqvMZLbmkuRe6iI2c3zJEzHOWecg6UHKjRmu0WOwG2AJoGNSXIHuQLB1aLtdLrl8t+OW/jf4OF0+akS0wX2kW2q7T+l+2iAqN8+5K6RwFGfDJWzidszq30NvspNzZWJ48n5er7sVj1j1ekrsoqzGJlzkoK3YXsL95PsbvYN8+izFPGnsI9FJUfGYY5OmE0ryyvylebnIgoZxROu/PIsIg1gFBHKBHOCAKsAb4zD4H2QGKDYwl3hBPmCMNusWO1Wn1nL34ZJKICo2gV0co3z0PklxRg5MiQ0XuP0zvjBRwWD/sssRRe/DoDep/j79JEGpT8snzyyvJ8wyxu75EwcfSMxtG5GEfPbpRUlBy5zLWilJ/yfiKzIJOSipJKE0mL3EXVWmO/+H6kxKdgs/4cJH4ZKsIcYb5Q8sshkfiQeKICo6q1FpHfogDTyOUdzGbHmzcwoPgHsMDq0LNod9M0EiNj/F2aSJ1izJFhlnLPkTMZJRUl7CnYw5bDW9hTuIftedvZU7jn50mjXu/PV6b8N5gUVxTXSG0dojrQMrwlMYExhDnCfJNCA6wBxIXEERMY43sfYAsg3BFOSEBI5eEX/jsEo8mhUk8owDRiPy6fS9QXt9KLQ5QbG2s63UP/P96nISNpEIwxlHnKcJW7yCnK8d2s63DZYXbm76TcW34kbPzijEe5p5zSilJKPaUcKD6Aq9zlu/qlxHPkzMepCgkI8Q2xBNgCCHOEEe4I//leG1YbdsuR/w20BRJoP3KJa5OgJnRt0pWQgJAjE0kt9kpDMiKNjQJMI+T1eFg+8zH6bH2ZAIuH3ZZEyi57g5Qep/u7NBFfYHB73eSX5VNUUcTugt0cKjlESUWJL1QcvflXbmkurjJXpXtslFSUkFOcc9L32TgRVouVxJBEEkITSI5KpmV4S9pGtsVpcx4zkfTo8EtkYCThjvAaq0mkMVGAaWRy9+9lz5vXM6B0JVhgVdh5JN/8JqHhGruW6uU1Xso8Zb5LaovcRazIXsGB4gPsLdz7811I/xs8jt5tNLMgs1rrsGAhLiSOIHsQTpuTYHswbSLbEGwP9p3tODqPw2Fz+G7uFemMJCYwBofNgcPmIMgeRJA9yHfmxG6x6xJYET9SgGlENi/7ipj/3EZ3cik1Aazt+hf6Dx2vISM5IR6vh+352/kx90cK3YWUVJRQWF7oO+vxywe9bc7dTE5RDhXm1IZcnDYnEc4IwgLCaBPZhtCA0CP35LAH+m4AFu4IJyYo5phJpXHBcUQ4I3y3QxeRhkUBphHwVFSw/N0H6Lfzn9gtXjKtzfAMfYuULin+Lk3qELfHzefbP+envJ/ILsrmUOkhit3FHCo5hKvcRamn9JS277A6aBfVjk7RnYgLiSM2KNZ3JuNo8AgKCKJdZDvf5bUKHiLyaxRgGrhDOXvY9+Z1pJatBgusDB9Ip5unEhIW6e/SpBqUecrYnrfddzv0o2dDCsoLKHQXUlpRSnFFsW9iapG7iNySXA6VHvIN21SYIxNZC92Fv7u/IHsQXWK6+G7wFRoQSoA14MhQjPXnINImog2dojsR7gzHYXVgt2q4RUSqlwJMA7Zh8WfEfXMH3ThMiXGwodfD9Lt0rL/LkuNwe9yUeEooqyjzBY6teVvZenhrpbki+4v3szN/p+/eIKc6RPO/whxhXNL2EuJD4okLiSPEHkJUYJQvsIQ7wrFb9bUhIv6nb6IGyFNRwfJ3JpCy63WsFsMuaxLmymn069TX36U1GuWeI1fHlFaUUu4tp6yijG1529hXtO/IxFbPkStlduTvILMgk4MlB096X06bkzYRbQi0B/ourw20BxIdGE2gPdA3+dRpcxIcEEx0YDTRgdE4bA7fBNajZ1BiAmMItAdWY0+IiNQMBZgG5sC+neS8PYLUsnSwwPKo/6PrqFcJDo3wd2kNmtvrZu7OuWw5vIV5mfPY6dp5UtuxYCE4IBinzUm4I5x+8f2IC46r9Jj7DlEdiAqM8gWSsIAwDc+ISKOjANOArF/4Ec0W3ElXXBQbJ5v6/I3+F4/2d1n1mtvrZrdrN1lFWZRWlFLiKeFw6WF25O8gvyyfck85WUVZZBzOOOazR8+EHL0MNyYwhi4xXX5+VowtgMSQRNpEtqFVeCuC7cGaKyIicoIUYBqACnc5K96+j5Tdb2G1GH6ytcb+x2n0Te7p79LqnExXpi+MZBVlsbdwLwXlBRS5i45caVNR6nsIXlZRFgXlBSe87WB7MGcnnU33Jt0Z3HowMYExCiMiIjVEAaaey9nzE7lvX0uqexNYYFnMJfS46TUCg0L8XZrfVHgr2OXaxfb87RwuPUxJRQnF7mLWHljL9/u+r/L2guxBNA9rToj9yP1HwhxhtApvRUxQjO9W792bdCchNEETXEVEaom+beuxtQs+oMXCu+hEAYUmiB/7P0bKkJv9XVatMMaQW5rL17u+ZnfBbg6XHmaXaxe5pblkFWX95i3k20S0ISQghJCAENpFtiPSGUlIQAihjlDfZFeH1UFkYCQJIQlEOCOwWnSzPxGRukQBph5yl5exauqdDMiZCcA2W1sCh71D33Zd/VxZzdhbuJdMV6bv5mrZRdksyFzA/pL9v/qZIHsQbSPaEhsc67sSJy44jrOTzqZzTOdarF5ERGqCAkw9k717G3lvX8uAis0ALGtyOT1GvdKghozcHjfZRdnMzZzLd3u+Y2XOyl9t2yq8FWc3P5swRxgtI1oSFxxH89DmRAdG6y6uIiINmAJMPbLm63/R5od7iacIF8FsS32alLTr/F1WlRljKHIXsS1vGz/m/khOcQ4HSw6yPX87O/J2UOA+duJsq/BWNAtrRkxgDE2DmtIivAVprdIICWg4wU1ERE6cAkw9UF5WyuqpdzBg//sAbLW3J/iad+ndppOfKzu+gvICsouyOVBygEMlh9jl2nXkJm6F+yiuKCa7KJsyT9nvbqdPXB/OaHYGF7a+kMTQxFqoXERE6gsFmDpu7/aNFM0YyYCKLQAsjRtG7xtfxOH0/91S9xfvZ/X+1ewr3Ed+WT6Zrky25W074Zu4hQaE0iu2F4mhicSHxJMQkkD7qPbEBccREhCiK3pERORX6ReiDlvz1du0++E+mllKcBHC9tOfZcAfrvFbPdlF2Xy/93uWZi1lf/F+1h1cR4X3+M/iiXBGEBMYQ2xwLBHOCLo16Uar8FaEOkJpEtSE2OBYguxBtXwEIiLSUCjA1EGlJUWsnXo7KQc/BAv8GNCZyOvfpWdSu1rZv9vjPnLlT0Em2/O2sz1/O6tyVpFZkHlM27YRbekY05FIZyRNgprQKboTnWI6ER0YXSu1iohI46QAU8fs2baBspnXkeLZDsCShGvpe8PzBDicNbZPYwzpB9JZsHsB6fvT2XhwI+Xe8mPa2Sw2OkZ3ZEDCADpEdaBDVAfaRdVOqBIREfklBZg6ZNUXU0le9gChlhIOE07m2c+Teu6V1bqPnKIcNh7ayI+5P7Ll8BayirLILsomtzS3UrtAWyAtw1vSLLQZHaI70CWmC33i+hDmCKvWekRERE6GAkwdUFpcyLrXb6X/4dlggU0BXWky8l/0aNb6lLZrjGF7/nYyXZmsyFnBkn1L2Ja37bhtHVYHaa3S6J/Qnx5Ne9AqvJWe4yMiInWWAoyf7d66FvfM6+nv3YnXWFjWfAT9Rj6DPcBxUts7WHKQ9P3pZBzO4D87/sMu165K6y1YaBfVjk7RnUiOSqZleEuaBjelZXhL3VNFRETqDQUYP1rxyWQ6r3mUEEsph4hg3/mTSD3rkiptw+P1sDl3M9/v/Z7Pt39+TGBx2py0DG9JtybdSE1MJSU+hcjAyGo8ChERkdqnAOMHJUUFrH/9FvrnfQEW2OjoTtzIf9EtseUJfX5f4T7WH1zPnO1zWJmzkoLyyneu7RDVgY7RHenepDtprdIUWEREpMFRgKllOzatgH+Por93Fx5jYXmrW+l/3ePY7L/9V+Eqd7Fw90LmZc5jXua8SuvCAsLoE9eHs5LO4vwW5+sSZhERafAUYGqJ8XpZ+ckkuq79O0GWcg4SSfagV0g9/aLK7YzBc+gQxTu3s33HajLX/8DBrO0UFB7G7jF09kB3j4UIWwjR9ghiHTGEWJ1QcRjMJ7j4BJefjlFERBqXpmPHEHLaaX7ZtwJMLSh0HWbz1Fvpl/8VWGBdYF8SRrxF2wN5HJ41i/LtOyhZu5bSnCy8+w9gqfAA4ATa//dVmQEK/vvaQ0ktHouIiMhRFYcP+23fCjA1bPuGZQR8OJJ+Zh8eY2FV5PW0KG9G3rARHMzOPqa9BfACh8KhMMyOaRZHZIv2xEclEREagyXAgSUgAEuAHYvdDnY7FnvAkT/brLV+fCIi0ngFdevmt30rwNQQ4/Wy4qMX6b7+CZy42XOwCa6sLoRtmsvRvFpht7C5OeyLhsymFvbE2YlJak/rtr25oP1FnB7TGZvV5tfjEBERqYsUYGqAK+8QW9+4kX4F31KU7eSnH1thzSnHwlYA1re08J++Fta3slDmsNC9aXcuaXsJf2j5B6ICo/xbvIiISD2gAFPNtq1djPODUbTaVkDGzlhMnh0r5XgssKibhY9Os3Ig2ka/uH6Mb3EeAxIG0Dqite56KyIiUgUKMNXEeL0se/9JOi5+ib3LIygsiACg3A5ze1n4KsVBu46p/KnlHzg36VydaRERETkFCjDVIP/wQbZOuQ7rms3sWR+FzVg4HAL/PsNK4Zk9SOs+lI9aDdKDEEVERKpJnQ4wkydP5plnniE7O5sePXowadIk+vfv7++yKslY9S3bp99GxHcQVRAEwPdd7ey76QJu7T+K5OhkP1coIiLS8NTZAPP+++8zfvx4pkyZQkpKCi+++CJpaWlkZGQQGxvr7/IwXi+fTb2Tgm++psd6K3YvHAy3sPvGP/DHGx8j3BHu7xJFREQaLIsxxvi7iONJSUmhX79+vPLKKwB4vV6SkpK4/fbb+ctf/vK7n3e5XERERJCfn094ePWGiaycXXzyzFC6LSwi5r+PIdrbO4kek96gaUyLat2XiIhIY3Kiv9918gxMeXk5q1atYsKECb5lVquVgQMHsmTJkuN+pqysjLKyMt97l6tmbqg//+UJ5M35hHP++9Dn3EgbcY89yvkDL9OVRCIiIrWkTt669eDBg3g8HuLi4iotj4uLI/s4d68FmDhxIhEREb5XUlJStddVUeHm8Jef0mkXeKyQc05XBnyzlHZ/uFzhRUREpBbVyQBzMiZMmEB+fr7vtXv37mrfh90eQPPLR7KybyBR09/inCkfYAsNrfb9iIiIyG+rk0NITZo0wWazkZOTU2l5Tk4O8fHxx/2M0+nE6XTWeG0pN91Lyk331vh+RERE5NfVyTMwDoeDPn36MG/ePN8yr9fLvHnzSE1N9WNlIiIiUhfUyTMwAOPHj2fEiBH07duX/v378+KLL1JUVMQNN9zg79JERETEz+psgLnqqqs4cOAADz30ENnZ2fTs2ZMvv/zymIm9IiIi0vjU2fvAnKqavA+MiIiI1IwT/f2uk3NgRERERH6LAoyIiIjUOwowIiIiUu8owIiIiEi9owAjIiIi9Y4CjIiIiNQ7CjAiIiJS7yjAiIiISL2jACMiIiL1Tp19lMCpOnqDYZfL5edKRERE5EQd/d3+vQcFNNgAU1BQAEBSUpKfKxEREZGqKigoICIi4lfXN9hnIXm9Xvbt20dYWBgWi6XatutyuUhKSmL37t16xlINU1/XDvVz7VA/1w71c+2oyX42xlBQUEBiYiJW66/PdGmwZ2CsVivNmzevse2Hh4frH0ctUV/XDvVz7VA/1w71c+2oqX7+rTMvR2kSr4iIiNQ7CjAiIiJS7yjAVJHT6eThhx/G6XT6u5QGT31dO9TPtUP9XDvUz7WjLvRzg53EKyIiIg2XzsCIiIhIvaMAIyIiIvWOAoyIiIjUOwowIiIiUu8owFTR5MmTadWqFYGBgaSkpLB8+XJ/l1RnTZw4kX79+hEWFkZsbCyXXnopGRkZldqUlpYyZswYYmJiCA0NZejQoeTk5FRqk5mZyZAhQwgODiY2NpZ77rmHioqKSm2+/fZbevfujdPppF27dkybNq2mD6/OevLJJ7FYLIwbN863TP1cPfbu3cu1115LTEwMQUFBdOvWjZUrV/rWG2N46KGHSEhIICgoiIEDB7J169ZK28jNzWX48OGEh4cTGRnJqFGjKCwsrNRm3bp1nHnmmQQGBpKUlMTTTz9dK8dXV3g8Hh588EFat25NUFAQbdu25bHHHqv0bBz1ddUtWrSIiy66iMTERCwWC5988kml9bXZpx988AEdO3YkMDCQbt268cUXX1T9gIycsPfee884HA7z5ptvmo0bN5qbb77ZREZGmpycHH+XVielpaWZt956y2zYsMGkp6ebCy+80LRo0cIUFhb62tx2220mKSnJzJs3z6xcudIMGDDAnHbaab71FRUVpmvXrmbgwIFmzZo15osvvjBNmjQxEyZM8LXZvn27CQ4ONuPHjzebNm0ykyZNMjabzXz55Ze1erx1wfLly02rVq1M9+7dzZ133ulbrn4+dbm5uaZly5Zm5MiRZtmyZWb79u3mq6++Mtu2bfO1efLJJ01ERIT55JNPzNq1a83FF19sWrdubUpKSnxtBg8ebHr06GGWLl1qvvvuO9OuXTszbNgw3/r8/HwTFxdnhg8fbjZs2GBmzpxpgoKCzD/+8Y9aPV5/evzxx01MTIyZPXu22bFjh/nggw9MaGioeemll3xt1NdV98UXX5gHHnjAfPTRRwYwH3/8caX1tdWn33//vbHZbObpp582mzZtMn/9619NQECAWb9+fZWORwGmCvr372/GjBnje+/xeExiYqKZOHGiH6uqP/bv328As3DhQmOMMXl5eSYgIMB88MEHvjabN282gFmyZIkx5sg/OKvVarKzs31tXnvtNRMeHm7KysqMMcbce++9pkuXLpX2ddVVV5m0tLSaPqQ6paCgwLRv397MnTvXnH322b4Ao36uHvfdd58544wzfnW91+s18fHx5plnnvEty8vLM06n08ycOdMYY8ymTZsMYFasWOFr85///MdYLBazd+9eY4wxr776qomKivL1+9F9JycnV/ch1VlDhgwxN954Y6Vll19+uRk+fLgxRn1dHf43wNRmn/7xj380Q4YMqVRPSkqKufXWW6t0DBpCOkHl5eWsWrWKgQMH+pZZrVYGDhzIkiVL/FhZ/ZGfnw9AdHQ0AKtWrcLtdlfq044dO9KiRQtfny5ZsoRu3boRFxfna5OWlobL5WLjxo2+Nr/cxtE2je3vZcyYMQwZMuSYvlA/V4/PPvuMvn37cuWVVxIbG0uvXr14/fXXfet37NhBdnZ2pT6KiIggJSWlUj9HRkbSt29fX5uBAwditVpZtmyZr81ZZ52Fw+HwtUlLSyMjI4PDhw/X9GHWCaeddhrz5s1jy5YtAKxdu5bFixdzwQUXAOrrmlCbfVpd3yUKMCfo4MGDeDyeSl/wAHFxcWRnZ/upqvrD6/Uybtw4Tj/9dLp27QpAdnY2DoeDyMjISm1/2afZ2dnH7fOj636rjcvloqSkpCYOp8557733WL16NRMnTjxmnfq5emzfvp3XXnuN9u3b89VXXzF69GjuuOMO3n77beDnfvqt74js7GxiY2Mrrbfb7URHR1fp76Kh+8tf/sLVV19Nx44dCQgIoFevXowbN47hw4cD6uuaUJt9+mttqtrnDfZp1FK3jBkzhg0bNrB48WJ/l9Lg7N69mzvvvJO5c+cSGBjo73IaLK/XS9++fXniiScA6NWrFxs2bGDKlCmMGDHCz9U1LLNmzWL69OnMmDGDLl26kJ6ezrhx40hMTFRfi4/OwJygJk2aYLPZjrlyIycnh/j4eD9VVT+MHTuW2bNns2DBApo3b+5bHh8fT3l5OXl5eZXa/7JP4+Pjj9vnR9f9Vpvw8HCCgoKq+3DqnFWrVrF//3569+6N3W7HbrezcOFCXn75Zex2O3FxcernapCQkEDnzp0rLevUqROZmZnAz/30W98R8fHx7N+/v9L6iooKcnNzq/R30dDdc889vrMw3bp147rrruOuu+7ynWFUX1e/2uzTX2tT1T5XgDlBDoeDPn36MG/ePN8yr9fLvHnzSE1N9WNldZcxhrFjx/Lxxx8zf/58WrduXWl9nz59CAgIqNSnGRkZZGZm+vo0NTWV9evXV/pHM3fuXMLDw30/JqmpqZW2cbRNY/l7Of/881m/fj3p6em+V9++fRk+fLjvz+rnU3f66acfcxuALVu20LJlSwBat25NfHx8pT5yuVwsW7asUj/n5eWxatUqX5v58+fj9XpJSUnxtVm0aBFut9vXZu7cuSQnJxMVFVVjx1eXFBcXY7VW/nmy2Wx4vV5AfV0TarNPq+27pEpTfhu59957zzidTjNt2jSzadMmc8stt5jIyMhKV27Iz0aPHm0iIiLMt99+a7Kysnyv4uJiX5vbbrvNtGjRwsyfP9+sXLnSpKammtTUVN/6o5f3Dho0yKSnp5svv/zSNG3a9LiX995zzz1m8+bNZvLkyY3q8t7j+eVVSMaon6vD8uXLjd1uN48//rjZunWrmT59ugkODjb/+te/fG2efPJJExkZaT799FOzbt06c8kllxz3MtRevXqZZcuWmcWLF5v27dtXugw1Ly/PxMXFmeuuu85s2LDBvPfeeyY4OLjBXtp7PCNGjDDNmjXzXUb90UcfmSZNmph7773X10Z9XXUFBQVmzZo1Zs2aNQYwzz//vFmzZo3ZtWuXMab2+vT77783drvdPPvss2bz5s3m4Ycf1mXUtWHSpEmmRYsWxuFwmP79+5ulS5f6u6Q6Czju66233vK1KSkpMX/6059MVFSUCQ4ONpdddpnJysqqtJ2dO3eaCy64wAQFBZkmTZqYu+++27jd7kptFixYYHr27GkcDodp06ZNpX00Rv8bYNTP1ePzzz83Xbt2NU6n03Ts2NH885//rLTe6/WaBx980MTFxRmn02nOP/98k5GRUanNoUOHzLBhw0xoaKgJDw83N9xwgykoKKjUZu3ateaMM84wTqfTNGvWzDz55JM1fmx1icvlMnfeeadp0aKFCQwMNG3atDEPPPBApUtz1ddVt2DBguN+J48YMcIYU7t9OmvWLNOhQwfjcDhMly5dzJw5c6p8PBZjfnFrQxEREZF6QHNgREREpN5RgBEREZF6RwFGRERE6h0FGBEREal3FGBERESk3lGAERERkXpHAUZERETqHQUYERERqXcUYERERKTeUYARERGRekcBRkREROodBRgRERGpd/4fpIVkmpE1SyMAAAAASUVORK5CYII=",
                        "text/plain": [
                            "<Figure size 640x480 with 1 Axes>"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                }
            ],
            "source": [
                "import jax.numpy as jnp\n",
                "from jax.random import PRNGKey\n",
                "from functools import partial\n",
                "import jax\n",
                "\n",
                "\n",
                "@partial(jax.jit, static_argnames=(\"num_round\"))\n",
                "def pull_all(bandit, key: PRNGKey, num_round: int):\n",
                "\n",
                "    best_rew_mean = jnp.max(bandit.rew_mean)\n",
                "    A = bandit.A\n",
                "\n",
                "    def body_fn(round, args):\n",
                "        key, rew_sum, num_arm_pulled, pseudo_rew_gaps = args\n",
                "\n",
                "        # アームを順番に選択します\n",
                "        arm = round % A\n",
                "        key, rew = bandit.pull(bandit, key, arm)\n",
                "\n",
                "        # 記録の更新\n",
                "        num_arm_pulled = num_arm_pulled.at[arm].add(1.0)\n",
                "        rew_sum = rew_sum.at[arm].add(rew)\n",
                "\n",
                "        # 結果の記録\n",
                "        rew_gap = best_rew_mean - bandit.rew_mean[arm]\n",
                "        pseudo_rew_gaps = pseudo_rew_gaps.at[round].set(rew_gap)\n",
                "        return key, rew_sum, num_arm_pulled, pseudo_rew_gaps\n",
                "\n",
                "    pseudo_rew_gaps = jnp.zeros(num_round)\n",
                "    rew_sum = jnp.zeros(A)\n",
                "    num_arm_pulled = jnp.ones(A)\n",
                "    args = key, rew_sum, num_arm_pulled, pseudo_rew_gaps\n",
                "    args = jax.lax.fori_loop(0, num_round, body_fn, args)\n",
                "    key, rew_sum, num_arm_pulled, pseudo_rew_gaps = args\n",
                "    return key, rew_sum, num_arm_pulled, pseudo_rew_gaps\n",
                "\n",
                "\n",
                "@partial(jax.jit, static_argnames=(\"m\", \"num_round\"))\n",
                "def ETC(bandit, key: PRNGKey, m: int, num_round: int):\n",
                "\n",
                "    best_rew_mean = jnp.max(bandit.rew_mean)\n",
                "    A = bandit.A\n",
                "\n",
                "    # 最初にmA回引きます\n",
                "    key, rew_sum, num_arm_pulled, pseudo_rew_gaps = pull_all(bandit, key, m * A)\n",
                "    greedy_arm = jnp.argmax(rew_sum / num_arm_pulled)\n",
                "\n",
                "    def body_fn(round, args):\n",
                "        key, pseudo_rew_gaps = args\n",
                "\n",
                "        # 貪欲に引きます\n",
                "        key, _ = bandit.pull(bandit, key, greedy_arm)\n",
                "\n",
                "        # 結果の記録\n",
                "        rew_gap = best_rew_mean - bandit.rew_mean[greedy_arm]\n",
                "        pseudo_rew_gaps = pseudo_rew_gaps.at[round].set(rew_gap)\n",
                "        return key, pseudo_rew_gaps\n",
                "\n",
                "    # 貪欲に引きます\n",
                "    pseudo_rew_gaps = jnp.hstack((pseudo_rew_gaps, jnp.zeros(num_round - m * A)))\n",
                "    args = key, pseudo_rew_gaps\n",
                "    args = jax.lax.fori_loop(m * A, num_round, body_fn, args)\n",
                "    key, pseudo_rew_gaps = args\n",
                "    return key, pseudo_rew_gaps\n",
                "\n",
                "\n",
                "@partial(jax.jit, static_argnames=(\"num_round\"))\n",
                "def UCB(bandit, key: PRNGKey, num_round: int):\n",
                "\n",
                "    best_rew_mean = jnp.max(bandit.rew_mean)\n",
                "    A = bandit.A\n",
                "\n",
                "    # 最初に全部1回ずつ引きます\n",
                "    key, rew_sum, num_arm_pulled, pseudo_rew_gaps = pull_all(bandit, key, A)\n",
                "\n",
                "    def body_fn(round, args):\n",
                "        key, rew_sum, num_arm_pulled, pseudo_rew_gaps = args\n",
                "        rew_mean = rew_sum / num_arm_pulled\n",
                "\n",
                "        # UCBスコアを計算してアームを引きます\n",
                "        delta = 1 / (num_round ** 2)\n",
                "        ucb = rew_mean + jnp.sqrt(2 * jnp.log(1 / delta) / num_arm_pulled)\n",
                "        arm = jnp.argmax(ucb)\n",
                "        key, rew = bandit.pull(bandit, key, arm)\n",
                "\n",
                "        # 記録の更新\n",
                "        num_arm_pulled = num_arm_pulled.at[arm].add(1.0)\n",
                "        rew_sum = rew_sum.at[arm].add(rew)\n",
                "\n",
                "        # 結果の記録\n",
                "        rew_gap = best_rew_mean - bandit.rew_mean[arm]\n",
                "        pseudo_rew_gaps = pseudo_rew_gaps.at[round].set(rew_gap)\n",
                "        return key, rew_sum, num_arm_pulled, pseudo_rew_gaps\n",
                "\n",
                "    # UCBを実行します\n",
                "    pseudo_rew_gaps = jnp.hstack((pseudo_rew_gaps, jnp.zeros(num_round - A)))\n",
                "    args = key, rew_sum, num_arm_pulled, pseudo_rew_gaps\n",
                "    args = jax.lax.fori_loop(A, num_round, body_fn, args)\n",
                "    key, rew_sum, num_arm_pulled, pseudo_rew_gaps = args\n",
                "    return key, rew_sum, num_arm_pulled, pseudo_rew_gaps\n",
                "\n",
                "\n",
                "@partial(jax.jit, static_argnames=(\"num_round\"))\n",
                "def EXP3(bandit, key: PRNGKey, num_round: int):\n",
                "\n",
                "    best_rew_mean = jnp.max(bandit.rew_mean)\n",
                "    A = bandit.A\n",
                "    eta = 2 * jnp.sqrt(jnp.log(A) / (num_round * A))\n",
                "\n",
                "    def body_fn(round, args):\n",
                "        key, rew_sum, pseudo_rew_gaps = args\n",
                "\n",
                "        key, sample_key = jax.random.split(key)\n",
                "        arm = jax.random.categorical(sample_key, eta * rew_sum)\n",
                "        key, rew = bandit.pull(bandit, key, arm)\n",
                "\n",
                "        # 記録の更新\n",
                "        prob = jax.nn.softmax(eta * rew_sum)[arm]\n",
                "        rew_t_estimate = 1 - (1 - rew) / prob\n",
                "        rew_sum = rew_sum.at[arm].add(rew_t_estimate)\n",
                "\n",
                "        # 結果の記録\n",
                "        rew_gap = best_rew_mean - bandit.rew_mean[arm]\n",
                "        pseudo_rew_gaps = pseudo_rew_gaps.at[round].set(rew_gap)\n",
                "        return key, rew_sum, pseudo_rew_gaps\n",
                "\n",
                "    # EXP3を実行します\n",
                "    rew_sum = jnp.zeros(A)\n",
                "    pseudo_rew_gaps = jnp.zeros(num_round)\n",
                "    args = key, rew_sum, pseudo_rew_gaps\n",
                "    args = jax.lax.fori_loop(A, num_round, body_fn, args)\n",
                "    key, rew_sum, pseudo_rew_gaps = args\n",
                "    return key, rew_sum, pseudo_rew_gaps\n",
                "\n",
                "\n",
                "# テスト\n",
                "\n",
                "import matplotlib.pyplot as plt\n",
                "key = PRNGKey(0)\n",
                "num_round = 10000\n",
                "\n",
                "key, *_, num_arm_pulled, pseudo_rew_gaps = pull_all(bernoulli_bandit, key, num_round)\n",
                "plt.plot(jnp.cumsum(pseudo_rew_gaps), label=\"Pull All\")\n",
                "\n",
                "key, *_, pseudo_rew_gaps = ETC(bernoulli_bandit, key, int(num_round / 10), num_round)\n",
                "plt.plot(jnp.cumsum(pseudo_rew_gaps), label=\"ETC\")\n",
                "\n",
                "key, *_, pseudo_rew_gaps = UCB(bernoulli_bandit, key, num_round)\n",
                "plt.plot(jnp.cumsum(pseudo_rew_gaps), label=\"UCB\")\n",
                "\n",
                "key, *_, pseudo_rew_gaps = EXP3(bernoulli_bandit, key, num_round)\n",
                "plt.plot(jnp.cumsum(pseudo_rew_gaps), label=\"EXP3\")\n",
                "\n",
                "plt.legend()"
            ]
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "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
}