Problem 3.3, SP9.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SP 9, \"Mass Cash\": Problem 3\n",
    "\n",
    "The \"Mass Cash\" lottery randomly chooses $5$ of the numbers from $1, 2, ..., 35$ each day (without repeitions within the choice of $5$ numbers). Suppose that we want to know how long it will take until all numbers have been chosen.\n",
    "\n",
    "Let $a_{j}$ be the average number of _additional_ days needed if we are missing $j$ numbers (so $a_{0} = 0$ and $a_{35}$ is the average number of days needed to collect all $35$ numbers). Find a recursive formula for the $a_{j}$'s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = \"retina\"\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import dependencies\n",
    "import random\n",
    "import pandas as pd\n",
    "from functools import lru_cache\n",
    "from collections import OrderedDict\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from math import floor, ceil, factorial, comb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Solution \n",
    "\n",
    "The `generate_sim_data` function runs the simulation and produces a dataframe indexed on $j$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sim_mass_cash_iterate(balls: tuple, unused_balls: set, draws_per_day: int = 5) -> int:\n",
    "    \"\"\"\n",
    "    Simulation counting the number of days required to have appearance of each ball\n",
    "    at least once given each day's draw consists of 5 balls without replacement\n",
    "    \"\"\"\n",
    "    days = 0\n",
    "    while unused_balls:\n",
    "        draw = set(random.sample(balls, draws_per_day))\n",
    "        unused_balls -= draw\n",
    "        days += 1\n",
    "\n",
    "    return days    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_sim_data(num_balls: int = 35, draws_per_day: int = 5, N: int = 50_000) -> pd.DataFrame:\n",
    "    \"\"\"\n",
    "    Returns dataframe of simulation data parameterized by the arguments to the function\n",
    "    \"\"\"\n",
    "    # Mass Cash numbered balls\n",
    "    balls = tuple(x for x in range(1, num_balls + 1))\n",
    "\n",
    "    # Run simulation lots of times and average\n",
    "    trial_data = np.array(\n",
    "        [sim_mass_cash_iterate(balls, set(balls), draws_per_day=draws_per_day) \n",
    "         for _ in range(N)]\n",
    "    )\n",
    "    \n",
    "    # Generate data\n",
    "    data = OrderedDict()\n",
    "    unused_balls = list(balls)\n",
    "    while len(data) < num_balls:\n",
    "        data[len(unused_balls)] = np.array(\n",
    "            [sim_mass_cash_iterate(balls, set(unused_balls), draws_per_day=draws_per_day) \n",
    "             for _ in range(N)]\n",
    "        ).mean()\n",
    "        unused_balls.pop()\n",
    "\n",
    "    # Add base case\n",
    "    data[0] = 0.0\n",
    "    \n",
    "    # Gather data\n",
    "    df = pd.DataFrame.from_dict(data, orient='index', columns=['a_j'])\n",
    "    df['delta'] = df['a_j'] - df['a_j'].shift(-1)\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running the simulation yields the following estimate for $a_{35}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "27.71368"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Average number of days needed to collect all numbers\n",
    "df = generate_sim_data()\n",
    "df.loc[35, 'a_j']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's construct a plot for all $a_{j}$!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZoAAAPCCAYAAADxnW1SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAABYlAAAWJQFJUiTwAACMlElEQVR4nOzdd3hfd302/vtIsjzlvUc8Eo84sZM4ccIMIYu9AgEKFHjo01Joy2xpnxZaOn9dQAcUuoBSoJRNITSDmbAynGEnjp1hO/FesmVZtvb5/SFFiRMn8ZD11Xi9rkuXdD7nnO/3trGNfPuT9ynKsgwAAAAAAJyoqkoHAAAAAABgYFM0AwAAAABwUhTNAAAAAACcFEUzAAAAAAAnRdEMAAAAAMBJUTQDAAAAAHBSFM0AAAAAAJwURTMAAAAAACdF0QwAAAAAwElRNAMAAAAAcFIUzQAAAAAAnBRFMwAAAAAAJ0XRfJKKovhCURRfqHQOAAAAAIBKqal0gEFgyYoVK1YkeUOlgwAAAAAAnKTiRG6yoxkAAAAAgJOiaAYAAAAA4KQomgEAAAAAOCmKZgAAAAAAToqiGQAAAACAk6JoBgAAAADgpCiaAQAAAAA4KYpmAAAAAABOSk2lAwwlnZ2dqa+vT2NjY1paWlKWZaUjDQpFUWT48OGpq6vLxIkTU1Xl308AAAAAoC8pmvtIZ2dnNm/enEOHDlU6yqBTlmWam5vT3NycpqamzJkzR9kMAAAAAH1I0dxH6uvrc+jQodTU1GT69OkZPXq0MrSXdHZ2pqmpKTt27MihQ4dSX1+fyZMnVzoWAAAAAAwZms4+0tjYmCSZPn166urqlMy9qKqqKnV1dZk+fXqSR3+uAQAAAIC+oe3sIy0tLUmS0aNHVzjJ4PXIz+0jP9cAAAAAQN9QNPeRRx78ZyfzqVMURZJ4yCIAAAAA9DGtJ4PGI0UzAAAAANC3FM0AAAAAAJwURTMAAAAAACdF0QwAAAAAwElRNAMAAAAAcFIUzQAAAAAAnBRFMwAAAAAAJ0XRTMV88YtfzBvf+MYsWrQodXV1mTBhQlauXJnPfOYzlY4GAAAAAByHmkoHYGhqbGzMm9/85qxcuTLPfe5zM3Xq1OzatSvf+ta38ra3vS27du3K7/7u71Y6JgAAAABwDBTNVERRFNmyZUumT59+xPqf/dmfZdGiRfnMZz6jaAYAAACAAULR3E/M+71rKh3hmG36y5ec9GuMGTMmY8aMecL6jBkzMnPmzNTX1x+x/gd/8Ae55ZZbcsMNN5z0ewMAAAAAvUvRTEXs27cvH//4x3PNNddk/fr1OXDgQDo7O3vOn3feeUdc/zu/8zuprq7u65gAAAAAwDFQNNPnVq9enSuvvDI7d+7MhRdemNe//vWZNGlSampqsnHjxnzuc5/LOeecc8Q948ePr0xYAAAAAOBpKZr7id4YRzFQ/PIv/3L279+fH/7wh7nkkkuOOPeHf/iHSZILLrigZ23//v2ZMGFCbr311iPWAQAAAID+oarSARhaNm/enNWrV+eSSy55Qsm8f//+fOYzn0mSnH/++T3rd955Z2pqanL22Wf3ZVQAAAAA4BgpmulTI0aMSJJs2LAhbW1tPet79+7N6173umzZsiU1NTU599xze87deeedWbx4cc+9AAAAAED/YnQGfWrKlCm59NJL84Mf/CAXXXRRLr/88mzfvj3/+7//m+c///mpqqrK0qVLjyiV77zzziOKZwAAAACgf7GjmT73pS99Kb/8y7+czZs355Of/GQefPDBfPSjH83v//7vp7Oz8wlzmO+66y5FMwAAAAD0Y3Y00+emTJmSz33uc0c9V5blEcdtbW1Zu3atohkAAAAA+jE7munX1q5dm9bWVkUzAAAAAPRjimb6tTvvvDMzZ87M5MmTKx0FAAAAAHgSimb6tVtvvTUXXXRRpWMAAAAAAE9B0Uy/dPjw4axatSpf+9rX8qIXvajScQAAAACAp6Bopl/6yEc+kssvvzwve9nL8ta3vrXScQAAAACAp1BT6QBwNB/84AfzwQ9+sNIxAAAAABgiOjrLbNt/OA/tPZSW9o5cdua0SkcaUBTNAAAAAMCQ0N7Rma37D2fT3kN5aG9TNu3p+rxxb1M21x9KW0eZJJk7aZSi+TgpmgEAAACAQaOtozNb9h3Opr1N2bSnKQ/tPZRNe7s+b64/lPbO8mlfY8u+w2nr6MywapOHj5WiGQAAAAAYUFraO7K5/nDXruTu3ckbu0vlrfsPp+MYyuSjmVI3PPMmjcrcSaNzqLUj40Yqmo+VohkAAAAA6Hea2zqyuf5QT4H8yK7kTXubsm3/4Zxgl5xpY4dn7qTRmT9pdOZOHpV5k0Znbne5PGa4uvRE+Zlj0CjLE/zTBQAAAICKaG7r6B5x0T0z+TFfbz/QnBOte2aOG5G5k0Zn3uSuAnle99enTRyVUbUq0VPBz2ofKYoiZVmms7MzVVW23J8KjxTNRVFUOAkAAAAAj1WWZbY3NGfdjgO5d3tj1m4/kHXbD2TjnqYT2plcFMnMcSMzr3tH8iO7kudNHp3TJo7KiGHVvf+D4CkpmvvI8OHD09zcnKamptTV1VU6zqDU1NSUpOvnGgAAAIDKaG7ryP07D+be7Qe6CuXucrnhcNtxvU5VkcyeMKqrQO4ukudPHp25k0ZnzsSRGV6jTO5PFM19pK6uLs3NzdmxY0eSZPTo0SmKwu7bk1SWZcqyTFNTU8/PrSIfAAAA4NQryzI7D7Q8plBuzL3bD2TD7oPHvEu5KJK5E7t2IvfsSu7+PHvCqNTWmAwwUCia+8jEiRPT1NSUQ4cOZcuWLZWOM2iNGjUqEydOrHQMAAAAgEGlua0jD+w62D3yoqtQXrfjQPYdOvZdynUjanLm9LE5c0ZdzpwxNktmjM3iaXUZWWtn8mCgaO4jVVVVmTNnTurr69PY2JiWlhYPr+slRVFk+PDhqaury8SJE83ABgAAADhBZVlmV2PLEYXyvdsPZMOepnQc4zblokjmTRrdVShP7yqUz5xRl1njR/qv+wcxRXMfqqqqyuTJkzN58uRKRwEAAABgiGtp75ql/MjIi3u7x1/UN7Ue82vUDa/Jkkd2KHfvVl48vS6jatWOQ43/xQEAAABgkOrsLLPnYEu27j+cbfub83D9oazb0bVb+cHdB9N+rMOUk8ybNKq7TH50/MXsCXYp00XRDAAAAAADVFNLe7btP5yt+w9ne0Nzz9fbuovl7Q2H09ZxfONbR9dW94y7eGSn8pLpdRk9XJXIk/OrAwAAAAD6oY7OMrsaHymPm7vL48NHHDccPvaH8R3NaRNHHVEoL+3epVxVZZcyx0fRDAAAAAAV0Njclm37H78LuWsn8tb9h7PjQPMxP4DvqUwYNSwzx4/s+hg3ImdMq8vSGXVZPH1sxtilTC/xKwkAAAAAToE9B1uyYXfTUYvkbQ2H09jcftLvUVtdlRnjR2TmuK4iedb4EY+WyuNHZub4ER7MR5/wqwwAAAAATlJ9U2vWbG3I3VsbsnrL/qzZ0pBtDc0n/bqTx9R270R+tDh+bIk8efRwYy7oFxTNAAAAAHAcGg63dRfKDVmzdX9Wb2nIln2Hj/t1htdUZdb4JxbIj6zNGDciI4ZVn4IfAfQ+RTMAAAAAPInG5rbcs+1A1mxpyOqtDVmzZX827T10TPcOr6nK4ul1mT3hsTuSHymSR2Ti6NoUhd3IDA6KZgAAAABIcqi1vadUXtM9AmPDnqaUx/A8vtrqqpw5c2yWzxqXZbPHZfnscTljypjUVFed+uDQDyiaAQAAABhymts6snZ7907l7hEYD+w6mM5jKJVrqoosmVGXZbPGZ/nscVk2a1wWTatLbY1SmaFL0QwAAADAoNbS3pH1Oxq7CuXuERj37WxMxzG0ytVVRRZOHdNVKM8en+WzxmXx9Dqzk+FxFM0AAAAADBptHZ1Zv6Mxa7Z2jb9Ys6Uh63YcSFvH05fKRZEsnDqmZ6fy2bPGZemMsRlZq1SGp6NoBgAAAGBA6ewss7OxOQ/tPZSH6w9lc33X5417mrJuR2Na2zuP6XUWTBndPVO5q1heOmNsRg9Xl8GJ8DsHAAAAgH7nUGt7Hq4/lIcfUyY/1F0ob6k/nNaOYyuTHzF30qgsmzWue6by+Jw1a2zGjhh2itLD0KNoBgAAAKDPdXaW2X2wpWdXck+ZvLcpD9cfzp6DLSf82rMnjOwplJfPHpezZ47LuFFKZTiVFM0AAAAAnBKHWzuyed+ju5If+7G5/lBajnHExdFMHF2bORNHZe7EUTmt+2POxFFZPL0uE0fX9uKPAjgWimYAAAAATkhZdu1KfkKR3H28q/HEdyXXVBWZPWFkV5k86dEy+bSJozNn4sjUGXsB/YqiGQAAAIBjsq+pNbc9tC+3bqrPrZvqs257Yw63dZzw640fNewxBfJjPiaNyoxxI1NdVfRieuBUUjQDAAAA8ARlWWbLvsPdpXJXufzAroPH9RrVVUVmjR/ZUx6f9rgxF+NG2pUMg4WiGQAAAIB0dpZZv7Px0WJ5Y312HGh+2vvqRtQ8ZrTF6J4iee6kUZkxbkRqqqv6ID1QaYpmAAAAgCGoua0ja7Y25JaN9bltU31ue2hfGpvbn/Kemqoiy2aPy8p5E3PB3AlZMXdCJo8Z3keJgf5M0QwAAAAwBDQcbsvtD+3LLZu6iuW7NjektaPzKe8ZXVudFXMnZOW8iVk5b2LOnTM+I2ur+ygxMJAomgEAAAAGoe0Nh7t3K3fNV16/szFl+dT3TB4zPBfOf7RYXjK9zugL4JgomgEAAAAGuM7OMg/uPti9W3lfbtlYn637Dz/tfQsmj84F8x4tludOGpWiKPogMTDYKJoBAAAABpjW9s7cva0ht27senDfbQ/VZ/+htqe8p7qqyFkzx+aCuRNz4fwJOX/uxEypM18Z6B2KZgAAAIB+bl9Ta1Zvbchtm+pzy8b63Ll5f1ran3q+8shh1TnvtPE9u5XPO218Rg9XBQGnhj9dAAAAAPqRhkNtuXtbQ1ZvaciarfuzektDtux7+jEYE0fX5oK5E3Lh/Im5YN7EnDVzbIaZrwz0EUUzAAAAQIU0Nrfl7q0HcvfWhqze2pA1W/Zn095Dx3TvaRNHde9WnpCV8ydmweTR5isDFaNoBgAAAOgDTS3tWbv9QNdO5S37s3prQzbsbjqme2urq3LmzLE5d/a4rJzfNQpj2tgRpzgxwLFTNAMAAAD0ssOtHVm7/UBPobxmS0Me3H0wneXT3zususiS6WNz9qxxWT57XJbNGpdF0+pSW2MMBtB/KZoBAAAATkJzW0fW7WjMmu7RF6u3NOT+XQfTcQytcnVVkUXT6rJ81rgsm91VLC+eXpfhNdV9kByg9yiaAQAAAI5Ra3tn7tvZeMSD+tbvaEz7MZTKVUWycGpdlnXvUl42e1yWzhibEcOUysDAp2gGAAAAOIq2js7cv/NgT6G8ZmtD1m1vTGtH59PeWxTJgsmjs3z2+CzrHoGxdObYjKpVxQCDkz/dAAAAgCFv78GWrN/RmHt3NGb9jgNZt6Mx63c0pqX96UvlJJk/eXRPoXz2rHE5a+bY1I0YdopTA/QfimYAAABgyGhp78gDuw5m3fbGrOsulNftaMzuxpZjfo05E0dm+azxXTOVZ43LWbPGZdxIpTIwtCmaAQAAgEGnLMts3X8467uL5Hu3H8j6HY3ZsKfpmB7S94hZ40f2zFNeNqvrY8Lo2lOYHGBgUjQDAAAAA1pjc1tPobxux4Gerxub24/5NUYMq8qiaXVZMr0ui6ePzZnT67J4el0mjRl+CpMDDB6KZgAAAGBAaO/ozKa9TV2F8mNGX2zZd/i4XmfupFFZPK0uS2aMzZLpXeXy3EmjU11VnKLkAIOfohkAAADod3Y3tvTsTr63u1S+f9fBtB7jw/mSZNzIYT1F8pIZY7N4el0WT6vL6OHqEIDe5k9WAAAAoKK27j+cmzfszT3bDvSUy3sOth7z/TVVRc6YOqZn7MWSGXU5c/rYTBs7PEVhlzJAX1A0AwAAAH1q6/7D+cWDe/OLDXvzi417s7n+2EdfTB87Iktm1GXJ9O6xFzPqsmDymNTWVJ3CxAA8HUUzAAAAcEqdSLE8qrY6ix8ZezF9bM/X40fV9kFiAI6XohkAAADoVcdbLI8cVp0L5k3I+XMn5MzuB/TNmTAqVR7OBzBgKJoBAACAk3KixfIzFkzKMxZMyvLZ4zKs2ugLgIFM0QwAAAAcl5MtlpfNGmemMsAgo2gGAAAAntK2/Ye7SuUNe/OLDfV5uP7QU15/ZLE8MctmjVcsAwxyimYAAADgCIplAI6XohkAAACGOMUyACdL0QwAAABDRFmW2dvUmq37DueBXQdz80bFMgC9Q9EMAAAAg0RbR2d2NDRn6/7D2brvcLbuP5xt+w8fcdzS3vm0rzNiWFVWzpuoWAbgmCmaAQAAYIBoamk/ojR+fKG880BzOsvjf90Rw6pywdyJecaCrnJ5+WzFMgDHR9EMAAAA/UBZltlzsPUJO5Af+3XD4baTfp+xI2oyc/zIzJ4wKufOGadYBqBXKJoBAACgDzwy1mLLvieOs3jk+FjGWjyVokim1g3PrPEjM3P8yMyaMDKzuz/PHD8ys8aPTN2IYb30IwKARymaAQAAoJeVZZkt+w7n5o31uWXj3ty6aV827W1KeQJjLR6rtroqM8ePyKwJIx8tk7uL5FnjR2bGuJF2JgNQEYpmAAAAOEllWebB3U25pbtYvmVjfbY1NB/369SNqMms8SMz+7FF8oRHy+TJo4enqqo4BT8CADg5imYAAAA4Th2dZdbtONBdLHd97G1qfcp7Hhlr8dhdyLMfUybPHD8yY421AGCAUjQDAADA02jr6MzdWxt6SuVbN9XnQHP7U94zurY658+bmIvmT8yF8ydm2axxGTGsuo8SA0DfUjQDAADA4zS3deTOzft7iuVVD+3L4baOp7xn3MhhuXD+o8Xy0hljU1NtXjIAQ4OiGQAAgCHvYEt7bn9oX27unq981+aGtHZ0PuU9U+qGH1EsL5paZ34yAEOWohkAAIAhZ/+h1ty6aV/Pg/vu3nYgHZ3lU94ze8LIxxTLkzJv0qgUhWIZABJFMwAAAEPArsbmIx7ct25H49Pes2DK6Fw0f1Iumj8xK+dPzKzxI/sgKQAMTIpmAAAABpWyLLNl3+HcuunRYnnDnqanvKcokiXTx/aMwVg5b2Km1A3vo8QAMPApmgEAABiwyrLMjgPNWb2lIXdvbciarQ1Zs6Uhe5tan/K+6qoiy2aN6ymWL5g7MeNGDeuj1AAw+CiaAQAAGBDKsszOAy3dZfL+rs9bG7Ln4FOXyklSW1OVc+eMz0XzJ+ai+ZNy3mnjM3q4vxIDQG/x/6oAAAD0SzsPNGfNloas3vrobuXdjS3HdO+Y4TU577TxPQ/uWz57XEYMqz7FiQFg6BpQRXNRFJOSvCrJS5IsSzIrSWuSNUk+k+QzZVl2Pub6eUk2PsVL/ndZlq8/ZYEBAAA4JrsONPfsUF6zpevzrmMslUfXVufsWeOybNa4LJvd9XnepNGpqipOcWoA4BEDqmhOcnWSTybZnuSHSR5OMi3JVUn+LcmLiqK4uizL8nH33ZXkm0d5vbtPXVQAAACOZldjc9cO5S0HsmZr1wiMnQeOvVQ+a+ajhfKy2eMyX6kMABU30Irm+5K8PMk1j9u5/PtJbkny6nSVzl973H13lmX54b4KCQAAQJfdjS09Yy8eeWDfjgPNx3TvqNrqnDVzbJbNGp9ls7s+z588OtVKZQDodwZU0VyW5Q+eZH1HURSfSvLnSS7JE4tmAAAATrE9B7se1Hf3Y+Yqb284tlJ55LDuUvmRncqzxmXBlDFKZQAYIAZU0fw02ro/tx/l3MyiKN6eZFKSvUl+Xpbl6j5LBgAAMAht2384P7l/T268f3duf2hfth1jqTxiWFXX+IvHzFU+XakMAAPaoCiai6KoSfLm7sNrj3LJFd0fj73nR0neUpblw8f4Hque5NSSY4wJAAAwoDW1tOcXG/bmpvv35Kb7d+fB3U1Pe8/wmqru8Rfjsmz2+CybNS6nTxmdmuqqPkgMAPSVQVE0J/nLJGcn+W5Zltc9Zv1Qkj9N14MAN3SvLU/y4STPT/L9oijOLcvy6b87AgAAGGI6Osvcs60hN92/Jzfetzu3P7wvbR2Pf/b6o4bXVGXpI6Vy907lM6aMUSoDwBBQlOWTf5MwEBRF8a4kf59kXZJnl2VZfwz31CT5SZKLkrynLMu/P4n3X7VixYoVq1Y92YZnAACAgWPr/sP5yf27c+P9e/LTB/Zk/6G2J712eE1VLpw/Mc9dODnPOn1ylkyvUyoDwMB3QrOsBvSO5qIofjNdJfPaJJcdS8mcJGVZthdF8W/pKpov7n4NAACAIed4x2EsmV6XixdNyXMXTs7KeRMzYlh1HyUFAPqzAVs0F0XxniQfS3J3ukrmXcf5Eru7P4/uzVwAAAD9WUdnmbu3NuQnDxzbOIzJY4bn4oWT89xFk/PsMyZnat2IPkwLAAwUA7JoLorid9M1l/nOJFeUZbnnBF7mGd2fNzzlVQAAAAPciYzDuHjhlDxnYdc4jKI4of+CFgAYQgZc0VwUxYeS/EmSVUmufKpxGUVRrEhyZ1mWnY9bvyzJe7sPP3+qsgIAAFTCY8dh3Hj/7mx4mnEYZ84Ym+cunGwcBgBwwgZU0VwUxVvSVTJ3JLkpybuO8i/rm8qy/Gz31x9NsrAoip8l2dK9tjzJpd1ff6gsy5+d0tAAAACn2CPjMG7q3rV8x9OMw5hSNzzPPcM4DACg9wyoojnJ/O7P1Une8yTX/DjJZ7u//s8kr0qyMsmLkgxLsjPJl5N8vCzLm05VUAAAgFNp6/7Duem+3bnpgeMbh/HcRZOzeJpxGABA7yrK8sn/lZunVxTFqhUrVqxYtWpVpaMAAACDVFmW2Vx/OLduqs9tD9Xn5g312bDn6cdhXLxwcp67cEoumDfBOAwA4Fid0L9GD7QdzQAAAINeR2eZdTsO5NaN9bn1oX25bVN9dh5oecp7ptQNz3MXTs7FC6fk2WdMzpS64X2UFgBA0QwAAFBxzW0duXPz/ty2qT63btqX2x/al8aW9qe8Z3hNVS5aMKln1rJxGABAJSmaAQAA+tj+Q625bdO+3PpQfW7dWJ81Wxue8uF9SVI3vCYr5k7IynkTcsG8iTl3znjjMACAfkPRDAAAcAqVZZmt+7vmK9+6qWsMxn07Dz7tfVPrhmfl/Im5cN7EXDBvQpZMH5vqKjuWAYD+SdEMAADQizo7y6zf2dgzBuPWTfXZ3tD8tPedPmV0Lpw/MRfMnZgL50/M7AkjjcIAAAYMRTMAAMBJaG7ryJqtDV07ljfWZ9VD+3Kg+annK9dUFTl71risnDchK+dNzPlzJ2TSGA/vAwAGLkUzAADAcWg43JZVDz06BuOuzQ1p7eh8yntG11Z3z1fuGoNx7pzxGVXrr2MAwODhOxsAAICnsL3hcG7ZWN/18L5N9Vm/szHlUz+3L5PHDM+F8yfkgrkTs3LexJw5oy411VV9ExgAoAIUzQAAAN3KssymvYdyy8a9uXljfW7dVJ/N9Yef9r4Fk0fngu4xGCvnTczcSaPMVwYAhhRFMwAAMGQ98uC+WzbW55aN9bl5Y332HGx5ynuqq4qcNXNsd6k8IefPnZgpdeYrAwBDm6IZAAAYMto6Orse3NddLN+6qf5pH9w3YlhVVpz26G7l804bn9HD/VUKAOCxfHcEAAAMWodbO3LH5n09pfLtD+3P4baOp7ynbkRNLpw3MSvnT8yF8yfm7JnjUltjvjIAwFNRNAMAAIPGgea2rNq0LzdvrM8tG/dmzdaGtHU89ZP7Jo8Znou6S+WV8yZm8fS6VFeZrwwAcDwUzQAAwIC1u7Elt26q75mxfO+OAymfulfOnIkjs3LexO5yeVLmeXAfAMBJUzQDAAADxpZ9h3rGYNy8sT4bdjc97T0Lp47JyvldxfLKeRMzc/zIPkgKADC0KJoBAIB+qSzLPLi7qXu38t7cumlftu4//JT3VBXJWTPH9YzBWDlvQiaNGd5HiQEAhi5FMwAA0G9srj+U79+7M7/Y0LVreW9T61NeX1tdlXPmdBXLF86flBWnjU/diGF9lBYAgEcomgEAgIrp7CyzZmtDbli7Mzes3Zn1Oxuf8vpRtdU5f+6EXDiv6+F958wZnxHDqvsoLQAAT0bRDAAA9KmW9o78/MG9uWHtznzv3p3ZeaDlSa8dP2pYVs6b2FMsnzVzbGqqq/owLQAAx0LRDAAAnHINh9ryg/Vdu5Z/vH53mlo7jnpdbU1VnnPG5Dx/8ZRcOH9SFk4dk6qqoo/TAgBwvBTNAADAKbG5/lDPSIxbNtWno7M86nUTRg3LpUum5YqlU/PchVMyeri/pgAADDS+gwMAAHpFWR45b3ndjieftzx30qhccea0XLF0Ws6fO8E4DACAAU7RDAAAnLCW9o78YkN9bli7I99buys7DjQ/6bXnzhmfK5Z2lcsLp45JURiJAQAwWCiaAQCA49JwqC0/XL8rN9zbNW/5YEv7Ua+rranKs0+flCuWTs/lZ07N1LEj+jgpAAB9RdEMAAA8rS37HjNveWN92p9k3vL4UcNy6eKpuWLptFy8yLxlAIChwnd9AADAE5Rlmbu3HsgN93aVy/duP/Ck1542cVTPSIwLzFsGABiSFM0AAECSpLW9M7/YsDc3rN2Z7927M9sbnnze8jmzx3WXy9OzaJp5ywAAQ52iGQAAhrBdjc35+YNd5fKP1+9O45PNW66uyrPOmJQrlk7L5WdOyzTzlgEAeAxFMwAADBH1Ta1ZvWV/1mxpyOqtDVmzpSE7Djz5ruVxI4fl0iWPzlseY94yAABPwneKAAAwCDUcbsvdWxuyektD1mzdn9VbGrJl3+GnvW/OxJG54szpuXzp1KycNzHDzFsGAOAYKJoBAGCAO9jSnnu2NmTN1obctaUha7bsz6a9h47p3hHDqnL2zHF53qIpueKsaVk8rc68ZQAAjpuiGQAABpDDrR1Zu717p3L3CIwHdx9MWT79vbU1VVk6Y2yWzx6XZbPGZfns8Tl9yujU2LUMAMBJUjQDAEA/1dLekXXbG7vnKXeNv7h/18F0dD59q1xTVWTJjLosmzW+p1hePL3OKAwAAE4JRTMAAPQDbR2duW9nY88u5dVb9mf9jsa0dTx9qVxVJIum1XUVyrPHZ3l3qTxiWHUfJAcAAEUzAAD0uY7OMg/sOpjVW/ZnTfcD+9ZuP5DW9s6nvbcoktOnjMnyWeOybPa4LJ89LktnjMvIWqUyAACVo2gGAIA+sKuxOdfevSP/u2ZH7ty8P4fbOo7pvnmTRmX57EfHX5w1a1zGDPdtPAAA/YvvUAEA4BTZe7Al/3v3jlyzentu3rg3TzdaefaEkd2FclexfPbMcRk3aljfhAUAgJOgaAYAgF60r6k1192zI9es2Z6fPbj3SR/cN33siCybPS7ndM9VXjZrXCaOru3jtAAA0DsUzQAAcJIaDrXl+rU78p3V2/PTB/ak/SjlclEkF86bmJeeMzNXnDkt08eNqEBSAAA4NRTNAABwAhqb23LD2p25ZvX23Hj/7rR1HH3n8gVzJ+Sly2fkxctmZOpY5TIAAIOTohkAAI7RwZb2fP/enfnO6u358X2709reedTrzjttfF6ybEZesnxGZowb2ccpAQCg7ymaAQDgKRxqbc8P1u3KNau35wfrdqXlScrl5bPH9excnj1hVB+nBACAylI0AwDA4zS3deRH63fl26u35wf37srhto6jXnfWzLF5yfIZeemymTltknIZAIChS9EMAABJWto78uP1u3PNmu353tqdaWo9erm8ZHpdz1iMBVPG9HFKAADonxTNAAAMWa3tnfnJA7vzndXbc8M9O9PY0n7U686YOiYvXT4jL10+I2dMrevjlAAA0P8pmgEAGFLaOjrzswf35jt3bct19+zIgeajl8vzJ4/uLpdnZtG0MSmKoo+TAgDAwKFoBgBg0Gvv6MwvNtTnmjXbcu3dO7LvUNtRrztt4qi8dHnXWIylM8YqlwEA4BgpmgEAGJQ6OsvcsrE+31ndVS7vbWo96nWzxo/sKZeXzRqnXAYAgBOgaAYAYNDo7Cyz6uF9+c5d2/Ldu3dkd2PLUa+bMW5EXrysa+byuXPGK5cBAOAkKZoBABjQyrLMnZv35zurt+e7a7Zne0PzUa+bWje8p1xecdqEVFUplwEAoLcomgEAGHDKssw92w7k26u35ZrV27Nl3+GjXjd5TG1edHbXWIyV8yamWrkMAACnhKIZAIABoSzLrNvRmGtWb893Vm/Lpr2Hjnrd+FHD8qKzp+ely2fmovkTU1Nd1cdJAQBg6FE0AwDQrz2wqzHfvqurXH5wd9NRr6kbUZMXnDU9L10+I88+Y3KGKZcBAKBPKZoBAOh3Nu1pyndWb8t3Vm/Puh2NR71mzPCaXLF0Wl66fEaes3ByhtdU93FKAADgEYpmAAD6hc31h3LNmq6dy3dvPXDUa0YOq85lZ07NS5fPzCWLp2TEMOUyAAD0B4pmAAAqZtv+w/numu359urtuWvz/qNeM7ymKpcu6SqXn79kSkbV+hYWAAD6G9+lAwDQp3YdaM5312zPd1Zvz20P7TvqNbXVVbl40ZS87JwZuezMaRkz3LetAADQn/mOHQCAU27PwZb87907cs3qbbl5Y33K8onX1FQVec7CyXnp8pm5Yum0jBs5rO+DAgAAJ0TRDADAKbH/UGuuvXtHvrN6e36+YW86Op/YLldXFXnW6ZPy0uUz8oKzpmf8qNoKJAUAAE6WohkAgF5zoLkt19+zM99ZvS0/uX9P2o9SLhdFctH8iXnp8pl50dnTM2nM8AokBQAAepOiGQCAk7LzQHNuWLsz192zI7/YsDdtHUeZi5Fk5bwJPeXy1LEj+jglAABwKimaAQA4bg/uPpjr7tmR6+/ZmTs373/S686dMz4vXT4jL1k+IzPGjey7gAAAQJ9SNAMA8LQ6O8vctWV/rl+7M9ffsyMP7m560muXzRqXlyyfkZcsm5E5E0f1YUoAAKBSFM0AABxVa3tnfrFhb65fuyM3rN2ZnQdajnpddVWRi+ZPzAvOmp7Ll07LrPF2LgMAwFCjaAYAoMfBlvb8eP3uXL92R36wblcam9uPet3IYdV53qIpufKsabl0ydSMH1Xbx0kBAID+RNEMADDE7W5syffu7RqJ8dMH9qa1o/Oo100YNSyXnzktV541Pc85Y3JG1lb3cVIAAKC/UjQDAAxBm/Y05fq1XQ/zW/XwvpTl0a+bPWFkrlw6PVeeNS0XzJ2Qmuqqvg0KAAAMCIpmAIAhoCzL3L31QE+5vH5n45Nee+aMsbly6bS84KzpOXNGXYqi6MOkAADAQKRoBgAYpNo6OnPrxvpcv7ZrLMa2huajXldVJBfM63qY35VLp2XOxFF9nBQAABjoFM0AAIPIodb23Hjfnlx/z458f92uNBxuO+p1w2uq8tyFXQ/zu2zJ1EwaM7yPkwIAAIOJohkAYICrb2rtfpjfztx0/+60tB/9YX7jRg7LZUum5sqzpuXiRVMyqta3ggAAQO/wtwsAgAFo/6HWfGf19nz7rm25dVN9Op/kYX4zxo3ombe8cv7EDPMwPwAA4BRQNAMADBCt7Z354fpd+cbtW/ODdbvS2nH0ncuLpo3JlUun5wVnTc/Zs8Z6mB8AAHDKKZoBAPqxsixzx+b9+cbtW/Pt1duy/9ATZy4XRXL+aRNy5VnTcsXS6Zk/eXQFkgIAAEOZohkAoB/aXH8o37hja75xx9Zs3NN01GuWzx6Xq86blRcvn5GpdSP6OCEAAMCjFM0AAP1Ew+G2fHfN9nz99i25ddO+o14za/zIvPK8mXnVebNzxtQxfZwQAADg6BTNAAAV1NbRmR+v352v37El37t3V1rbnzh3eczwmrx42fS86rzZuWj+xFRVmbkMAAD0L4pmAIA+VpZl7trSkG/cviXfXr099U2tT7imuqrIxQsn51UrZueKM6dlZG11BZICAAAcG0UzAEAf2bLvUL55x9Z8/Y6t2bD76HOXz541Nq86b3Zefs7MTKkb3scJAQAAToyiGQDgFDrQ3Jb/XbM9X799a27eWH/Ua2aMG5FXnjcrV503Kwun1fVxQgAAgJOnaAYA6GVtHZ256f7d+drtW/O9tTvTcpS5y6Nrq/PCs2fk1Stm5aIFk1Jt7jIAADCAKZoBAHpBWZa5e+uBfO32Lfn2Xduy9yhzl6uK5LkLp+SqFbNyxdJpGVXrWzEAAGBw8LcbAICTsHX/4Xzzjq35xh1b88Cug0e9ZumMsblqxay8/JyZmTp2RB8nBAAAOPUUzQAAx6mxuS3/e/eOfOP2rfnFxr0pyydeM23s8Lzy3Fl51YpZWTJ9bN+HBAAA6EOKZgCAY9De0ZmbHtiTb9y+Ndev3ZHmtifOXR5VW50XnjU9r1oxK886fbK5ywAAwJChaAYAeBJlWeaebQfy9du35n/u2pY9B1uecE1VkTz7jMm5asWsXLl0ekYP9+0VAAAw9PibEADA42xvOJxv3rEt37hjS+7befS5y0um1+WqFbPyinNnZZq5ywAAwBCnaAYASHKwpT3X3r0j37hjS3724NHnLk+pG55XnDMzV62YnaUzzV0GAAB4hKIZABiy2js689MH9+Ybt2/JdffszOG2jidcM3JYdV5w1rS8asXsPPv0SamprqpAUgAAgP5N0QwADDlrtx3I12/fkm/dtS27G584d7kokmedPimvOm92Xnj29IwxdxkAAOAp+VsTADAk7Ghozrfu3Jpv3LE163Y0HvWaRdPG5KoVs/OKc2dmxriRfZwQAABg4FI0AwCDVlNLe667Z0e+fvvW/PTBPUeduzx5zPC84tyZedV5s3LWzLEpiqLvgwIAAAxwimYAYFDp6Czz0wf25Bt3bM21d+846tzlEcOqcuXS6blqxaw854zJ5i4DAACcJEUzADAo3Lv9QL5xx9Z8686t2Xng6HOXn7lgUl513qy88OzpqRsxrAIpAQAABidFMwAwYO080DV3+eu3P/nc5YVTx+RVK2bllefOyszx5i4DAACcCopmAGBAOdT6mLnLD+xJ51HnLtfm5efMylUrzF0GAADoC4pmAKDf6+gs8/MH9+brd2zJtXfvyKHWJ85dHl5TlSvPmp6rzpuV5y40dxkAAKAvKZoBgH5r/Y7GfP2OLfnWHduy40DzUa95xoKJuWrF7LzI3GUAAICKUTQDAP3Kodb2fGf19nzx5odz5+b9R73mjKlj8qrzZuWV583KLHOXAQAAKk7RDAD0C+t2HMgXb34437hjaxqb259wftLo2rzsnJm5asWsLJs1ztxlAACAfkTRDABUzOHWjlyzZnu+ePNDuf3h/U84X1tdlSuWTsurz5+V5y6ckmHmLgMAAPRLimYAoM/dt7MxX7z54Xz99i05cJTdy/Mnj84vXTgnr14xO5PGDK9AQgAAAI6HohkA6BPNbR357pqu2cu3PbTvCeeHVRd5wVnT84YLT8szT59kNAYAAMAAomgGAE6pB3Y15gs3P5yv3741DYfbnnB+7qRR+aULT8trzp+dyXYvAwAADEiKZgCg1zW3deTau3fkizc/nFs21T/hfE1V1+7lX7rwtDzr9EmpqrJ7GQAAYCBTNAMAveaBXQfzX7c8nK/dviX7Dz1x9/KciSPzSxeelqvPn5MpdXYvAwAADBaKZgDgpLS0P7p7+eaNT9y9XF1V5Iozp+UNF52W55wx2e5lAACAQUjRDACckA27u3Yvf3XVluw7yu7l2RMe2b08O1PHjqhAQgAAAPqKohkAOGat7Z257p6u3cs/37D3Ceerq4pctmRq3nDRabl44RS7lwEAAIYIRTMA8LQ27Wnq2b28t6n1CednjR+Z16+ck9eunJNpdi8DAAAMOYpmAOCoWts7c8PanfniLQ/lpw88cfdyVZFcumRa3njRabl40ZRU270MAAAwZCmaAYAjPLS3Kf91y+Z8ddXm7Dn4xN3LM8eNyOtWnpbXrpydGeNGViAhAAAA/Y2iGQBIWZb5/r278h8/35Sb7t/zhPNVRfL8xV2zly9ZPNXuZQAAAI6gaAaAIayzs8x19+zIP/zggdy7/cATzk8fOyKvXTknr185JzPH270MAADA0SmaAWAI6ugs89012/OPP7g/9+08eMS5okguWTQlb7hobp6/eEpqqqsqlBIAAICBQtEMAENIe0dnvrO6q2B+cHfTEedGDqvOGy86LW999rzMnjCqQgkBAAAYiBTNADAEtHV05lt3bssnfvhANu45smAeVVudNz9zXv7vc+dn8pjhFUoIAADAQKZoBoBBrLW9M1+/fUs+8aMHsrn+8BHn6obX5K3Pnpe3PXt+JoyurVBCAAAABgNFMwAMQi3tHfnKbVvyyR89mK37jyyYx46oydueMz//51nzM27UsAolBAAAYDBRNAPAINLc1pEv3fJwPvXjDdlxoPmIc+NHDcuvPndBfvmZczN2hIIZAACA3qNoBoBB4HBrR75w80P55xs3ZHdjyxHnJo2uza9evCBvesbcjBnu//oBAADoff62CQADWFNLez7/i4fyrzdtyJ6DrUecmzxmeH79eQvyhotOy6ha/5cPAADAqeNvnQAwADU2t+VzP38o/3bThuw71HbEuWljh+fXn3d6funC0zJiWHWFEgIAADCUKJoBYABpONyWz/50Uz79041pOHxkwTxz3Ii84/ln5OrzZyuYAQAA6FOKZgAYAPYfas2nf7Ixn/nppjS2tB9xbvaEkfmN55+RV6+YndqaqgolBAAAYChTNANAP1bf1Jp/u2lD/uNnm9LU2nHEuXmTRuU3nn9GXnnerAyrVjADAABQOYpmAOiHdje25F9v2pDP/+KhHHpcwbxgyuj81qVn5GXLZ6ZGwQwAAEA/oGgGgH5k54Hm/POPN+SLtzyU5rbOI84tnDomv3XZwrxk2YxUVxUVSggAAABPpGgGgH5g2/7D+dSPH8yXbt2c1vYjC+Yl0+vyrssW5oVnTU+VghkAAIB+SNEMABW0uf5QPvnjB/OV2zanraM84tzZs8bmXZcuzOVnTlMwAwAA0K8pmgGgArbsO5R//P4D+drtW9LeeWTBfM6c8Xn3ZWfk+YunpigUzAAAAPR/imYA6EP7mlrziR8+kM/9/KG0dhw5IuP8uRPy7ssW5rkLJyuYAQAAGFAUzQDQBw61tuczP92UT/3owTS2tB9x7qL5E/PuyxbmmadPUjADAAAwICmaAeAUauvozJdv25y/+9792d3YcsS5c+eMz+++cEmeefqkCqUDAACA3qFoBoBToCzLfHfNjvzt9euzcU/TEecWTBmdD7xgSV5w1jQ7mAEAABgUBlTRXBTFpCSvSvKSJMuSzErSmmRNks8k+UxZlp1Hue9ZST6Y5BlJRia5P8mnk/xjWZYdfZMegKHiZw/syV9duy53bWk4Yn3a2OF57+WL8przZ6emuqpC6QAAAKD3DaiiOcnVST6ZZHuSHyZ5OMm0JFcl+bckLyqK4uqyLMtHbiiK4hVJvpakOcl/J6lP8rIkH0vy7O7XBICTds+2hvzVtetz4327j1ivG1GTd15yRt76rHkZWVtdoXQAAABw6gy0ovm+JC9Pcs1jdy4XRfH7SW5J8up0lc5f614fm+Rfk3QkuaQsy9u61z+U5AdJXlMUxevLsvxSn/4oABhUHt57KB+5YX2+dee2I9Zra6ryf541L++45PSMH1VboXQAAABw6g2oorksyx88yfqOoig+leTPk1yS7qI5yWuSTEnyuUdK5u7rm4ui+GCS7yd5RxJFMwDHbc/Blnz8Bw/kCzc/lLaOnv+YJlVF8przZ+c9ly/KzPEjK5gQAAAA+saAKpqfRlv35/bHrF3a/fnao1x/Y5JDSZ5VFMXwsixbTmU4AAaPgy3t+bebNuRfb9yQptYjR/1fsXRaPvCCxVk4ra5C6QAAAKDvDYqiuSiKmiRv7j58bKm8uPvzfY+/pyzL9qIoNiY5K8mCJPc+zXusepJTS44vLQADVWt7Z/7rlofzD9+/P3ubWo84t3LehPzei5bk/LkTK5QOAAAAKmdQFM1J/jLJ2Um+W5bldY9ZH9f9ueFJ7ntkffwpygXAINDZWebbq7flI9ffl4frDx1xbtG0MfndFy7JpUumpiiKCiUEAACAyhrwRXNRFO9K8v4k65L88ql6n7Isz3+S91+VZMWpel8AKqcsy9x4/5789bXrcs+2A0ecmzluRN535eK86rxZqa5SMAMAADC0DeiiuSiK30zy90nWJrmsLMv6x13yyI7lcTm6R9b39346AAayuzbvz19duy4/e3DvEevjRw3Lbz7/jLzpGXMzYlh1hdIBAABA/zJgi+aiKN6T5GNJ7k5XybzrKJetT3JBkkVJjpix3D3XeX66Hh644ZSGBWDA2LD7YD5y/X25Zs32I9ZHDKvKrzxnft7+vNMzdsSwCqUDAACA/mlAFs1FUfxuuuYy35nkirIs9zzJpT9I8sYkL0zyX487d3GSUUluLMuy5RRFBWCA2HWgOX///fvzpVs3p6Oz7FmvriryupVz8u7LFmba2BEVTAgAAAD914Armoui+FCSP0nXDuUrjzIu47G+muSvkry+KIp/LMvytu7XGJHkz7qv+eSpzAtA/3aguS3/8uMN+fefbMzhto4jzr1k2Yy8/8pFWTBlTIXSAQAAwMAwoIrmoijekq6SuSPJTUneVRRPeADTprIsP5skZVkeKIriV9NVOP+oKIovJalP8vIki7vX/7tv0gPQnzS3deTzv3gon/jhA9l3qO2Ic89cMCm/+6IlOXfO+MqEAwAAgAFmQBXN6ZqpnCTVSd7zJNf8OMlnHzkoy/KbRVE8L8kfJHl1khFJHkjyviT/UJZlebQXAWBw6ugs8407tuZjN9yXrfsPH3HuzBlj83svWpKLF07OUf4hEwAAAHgSA6poLsvyw0k+fAL3/TTJi3s7DwADR1mW+cG6Xfnra9dn/c7GI87NmTgyv33l4rxs+cxUVSmYAQAA4HgNqKIZAE7Ezx7ck49ef19ue2jfEeuTRtfmty49I2+4aG5qa6oqlA4AAAAGPkUzAIPWrZvq89Hr78vPN+w9Yn10bXX+73MX5FcvXpAxw/1fIQAAAJwsf7sGYNC54+F9+egN9+Wm+/ccsT6susgbLjwtv3npwkypG16hdAAAADD4KJoBGDTu3tqQj95wX36wbtcR6zVVRa6+YHZ+4/lnZPaEURVKBwAAAIOXohmAAe/e7QfysRvuy/Vrdx6xXlUkrzpvdt592cKcNknBDAAAAKeKohmAAev+nY35u+/dn2vWbD9ivSiSl58zM+++bGEWTBlToXQAAAAwdCiaARhwNu5pyt9/7758665tKcsjz7142fS85/JFWTStrjLhAAAAYAhSNAMwYGyuP5S///79+cYdW9PReWTDfMXSaXnv5YuydObYCqUDAACAoUvRDEC/t3X/4Xz8Bw/kK7dtTvvjCubnL56S912xOMtmj6tQOgAAAEDRDEC/taOhOf/0owfypVs2p7Wj84hzz104Oe+5fFHOnzuhQukAAACARyiaAeh3dje25JM/ejCfv/mhtLYfWTBfNH9i3nfFoly0YFKF0gEAAACPp2gGoN+ob2rNP9/4YD73s4dyuK3jiHMrThuf91+5OM86fVKKoqhQQgAAAOBoFM0AVFzDobb8600b8pmfbkxT65EF8/LZ4/K+KxbleYumKJgBAACgn1I0A1AxB5rb8umfbMy/37QxjS3tR5xbOmNs3nfFolx25lQFMwAAAPRzimYA+lxTS3s++7NN+ZcbN6ThcNsR5xZNG5P3Xr4oLzhreqqqFMwAAAAwECiaAegzh1s78p+/2JRP/XhD6ptajzi3YMrovOfyRXnpshkKZgAAABhgFM0AnHLNbR354s0P559+9GD2HGw54tzcSaPy7ssW5uXnzExNdVWFEgIAAAAnQ9EMwCnT0t6RL9+6OR//4QPZeeDIgnnW+JF512Vn5KoVszNMwQwAAAADmqIZgF7X0Vnmq6s25x++/0C27j98xLnpY0fkNy89I6+9YE5qaxTMAAAAMBgomgHoVXdt3p8PfvPurNnacMT6lLrh+Y1LTs/rLzwtI4ZVVygdAAAAcCoomgHoFQ2H2vI316/LF25+OGX56PrE0bV5x/NOz5ueMTcjaxXMAAAAMBgpmgE4KWVZ5uu3b81ffPfe7G1q7VmvranKO553en7t4gUZPdz/3QAAAMBg5m/+AJyw+3Y25oPfvDu3bKw/Yv15i6bkT15xVuZOGl2hZAAAAEBfUjQDcNyaWtrzD9+/P//+k41p73x0TsaMcSPyRy9bmhecNT1FUVQwIQAAANCXFM0AHLOyLHPdPTvyx99em+0NzT3rNVVFfuU58/OuyxYakwEAAABDkDYAgGPy0N6m/NH/3JMfrd99xPqF8ybmT195dhZPr6tQMgAAAKDSFM0APKXmto78y40b8okfPpCW9s6e9Umja/P7Lz4zV62YZUwGAAAADHGKZgCe1I337c4f/c892binqWetKJI3XnRafufKJRk3algF0wEAAAD9haIZgCfY0dCcP71mba5Zvf2I9WWzxuXPXnl2zpkzvjLBAAAAgH5J0QxAj/aOznz2Z5vysRvuS1NrR8963YiafOAFi/OGi+amusqYDAAAAOBIimYAkiS3barPB795d9btaDxi/arzZuX/vfjMTKkbXqFkAAAAQH+naAYY4uqbWvOX/3tvvnzbliPWz5g6Jn/6irPzzNMnVSgZAAAAMFAomgGGqM7OMv992+b81bXrsv9QW8/6yGHVefflC/O2Z89PbU1VBRMCAAAAA4WiGWAIuntrQz74zbtz5+b9R6y/4Kxp+cOXnZVZ40dWJhgAAAAwICmaAYaQA81t+ej19+VzP9+UzvLR9TkTR+aPX35WLl0yrXLhAAAAgAFL0QwwBJRlmf+5a1v+7Jp7s7uxpWe9troqb3/egrzzkjMysra6ggkBAACAgUzRDDDIPbDrYP7wW3fnZw/uPWL9OWdMzh+/4qycPmVMhZIBAAAAg4WiGWCQOtzakY//8P78y40b0tbx6JyMqXXD86GXLs1Ll89IURQVTAgAAAAMFopmgEHoe2t35o/+555s3X+4Z62qSN76rPl57xULUzdiWAXTAQAAAIONohlgENlcfyh//O21+d69O49YX3Ha+PzpK8/OWTPHVSgZAAAAMJgpmgEGgdb2zvzbTzbkH75/f5rbOnvWx48alv/3oiW5+vw5qaoyJgMAAAA4NRTNAAPcvdsP5L3/fWfW7Wg8Yv11F8zJ775oSSaOrq1QMgAAAGCoUDQDDFAdnWX++cYH87Eb7jviYX9Lptflz191ds6fO7GC6QAAAIChRNEMMABt2tOU93/lrqx6aF/P2vCaqvzOCxbnrc+al5rqqgqmAwAAAIYaRTPAAFKWZb5w88P582vuzeG2jp71c+aMz0dfe05OnzKmgukAAACAoUrRDDBA7Ghozge+tjo33re7Z62mqsi7LluYd15yul3MAAAAQMUomgEGgG/duTUf+ubdOdDc3rO2cOqYfOx15+bsWeMqmAwAAABA0QzQr+1ras0Hv3V3rlm9vWetKJL/+5z5ef+VizNiWHUF0wEAAAB0UTQD9FM/XLcrH/ja6uxubOlZmz1hZP726nPyjAWTKpgMAAAA4EiKZoB+pqmlPX92zb35r1sePmL99Svn5IMvXZoxw/3RDQAAAPQv2gqAfuTWTfV5/5fvysP1h3rWJo8Znr969bJcdua0CiYDAAAAeHKKZoB+oLmtIx+74b78y00bUpaPrr/o7On581cty8TRtZULBwAAAPA0FM0AFXbPtoa877/vyvqdjT1rdSNq8qevODuvOHdmiqKoYDoAAACAp6doBqiQ9o7O/PONG/J337svbR2PbmN+zhmT8zdXL8+McSMrmA4AAADg2CmaASpg456mvO/Ld+aOh/f3rI0YVpXff/GZedNFc1NVZRczAAAAMHAomgH6UFmW+c9fPJS/+O69aW7r7Fk/77Tx+cjV52TBlDEVTAcAAABwYhTNAH1ke8PhfOCrq3PT/Xt61oZVF3nP5Yvy9osXpKa6qoLpAAAAAE6cohngFCvLMt+6c1s+9K2709jc3rO+eFpdPvq6c3LWzHEVTAcAAABw8hTNAKdQfVNrPvjNNfnumh09a0WR/NpzF+S9VyzKiGHVFUwHAAAA0DsUzQCnyPfv3Znf/dqa7DnY0rM2Z+LIfOTqc3Ph/IkVTAYAAADQuxTNAL2ssbktf/ade/Pft20+Yv2XLjwtf/CSMzNmuD96AQAAgMFF2wHQi27esDfv/8pd2bLvcM/alLrh+etXL8/zl0ytYDIAAACAU0fRDNALmts68pHr1+fffrIxZfno+kuWz8ifveLsTBhdW7lwAAAAAKeYohngJN29tSHv/e87c/+ugz1rY0fU5M9etSwvP2dmBZMBAAAA9A1FM8AJau/ozCd/9GD+/vv3p73z0W3MFy+akr9+9fJMHzeigukAAAAA+o6iGeAEPLj7YN735bty1+b9PWsjh1XnD15yZt540WkpiqJy4QAAAAD6mKIZ4Dh9846t+b2vr05zW2fP2orTxuejrz038yaPrmAyAAAAgMpQNAMco9b2zvzFd+/NZ3+2qWdtWHWR916xKG+/+PRUV9nFDAAAAAxNimaAY7DzQHPe+YXbs+qhfT1rp08ZnX/8pRVZOnNsBZMBAAAAVJ6iGeBp3Lxhb37ji3dkz8GWnrUXnT09f3P1ORkz3B+jAAAAABoSgCdRlmU+/dNN+Yvv3puOzjJJUlUkv/vCJfm1ixd44B8AAABAN0UzwFE0tbTn976+Jt++a1vP2qTRtfnHXzovzzpjcgWTAQAAAPQ/imaAx9mw+2B+/fOrct/Ogz1r58wZn0++cUVmjh9ZwWQAAAAA/ZOiGeAxrr9nR97/5bvS2NLes/aGi07LH71saYbXVFcwGQAAAED/pWgGSNLRWeajN6zPJ374YM9abU1V/uyVZ+e1F8ypYDIAAACA/k/RDAx59U2tefeX7shN9+/pWZs9YWQ+9abzc/ascRVMBgAAADAwKJqBIW3Nlob8+udXZev+wz1rFy+akr9/3bmZMLq2gskAAAAABg5FMzBk/fetD+dD37onre2dPWvvuvSMvPvyRamuKiqYDAAAAGBgUTQDQ05zW0f++Nv35L9u2dyzVjeiJh977bm5fOm0CiYDAAAAGJgUzcCQsnX/4bzj86uyektDz9qS6XX51JvOz7zJoyuYDAAAAGDgUjQDQ8ZP7t+T3/qv27PvUFvP2ivOnZn/76plGVXrj0MAAACAE6VZAQa9sizzyR8/mL+9bn06y661mqoiH3zJmXnLs+alKMxjBgAAADgZimZgUGtsbstvf+WuXHfPzp61KXXD809vXJGV8yZWMBkAAADA4KFoBgat+3c25u3/uSob9jT1rK2cNyGfeMOKTB07ooLJAAAAAAYXRTMwKH1n9bZ84Kurc6i1o2ft/zx7Xn7/xWdmWHVVBZMBAAAADD6KZmBQae/ozF/+77r820829qyNHFadv3z1srzi3FkVTAYAAAAweCmagUFjd2NLfvOLt+fmjfU9a/Mmjcqnfvn8LJk+toLJAAAAAAY3RTMwKKx6aF/e+YVV2XmgpWft8jOn5aOvOydjRwyrYDIAAACAwU/RDAxoZVnm8794KH/ynbVp6yiTJEWRvP+KRXnnJWekqqqocEIAAACAwU/RDAxYh1s78gffWJOv37G1Z238qGH5h9efl4sXTalgMgAAAIChRdEMDEgP7z2Ut39+Ve7dfqBn7exZY/PJN56fORNHVTAZAAAAwNCjaAYGnB+u25V3f+mOHGhu71m7+vzZ+dNXnp0Rw6ormAwAAABgaFI0AwNGZ2eZv//+/fmHH9yfsmscc2qrq/LHrzgrr185J0VhHjMAAABAJSiagQGh4VBb3vPfd+SH63f3rM0cNyL/9Kbzc+6c8ZULBgAAAICiGej/7t1+IG//z1V5uP5Qz9qzTp+Uf/yl8zJpzPAKJgMAAAAgUTQD/dwNa3fm3V+6I4daO3rWfv15p+e3r1yUmuqqCiYDAAAA4BGKZqBfKssy/3Ljhvzltet65jGPGV6Tv716eV549ozKhgMAAADgCIpmoN9pbe/MH3xjTb6yakvP2pyJI/Pvb1mZRdPqKpgMAAAAgKNRNAP9Sn1Ta37986tyy8b6nrWV8ybkU2863zxmAAAAgH5K0Qz0Gw/saszbPnvbEQ/9e/WK2fmLq87O8JrqCiYDAAAA4KkomoF+4cb7duc3vnh7GpvbkyRFkfzuC5fk7RcvSFEUFU4HAAAAwFNRNAMV97mfb8off3ttOjq7nvo3clh1/u715+YFZ02vcDIAAAAAjoWiGaiY9o7O/Ml31uZzP3+oZ23GuBH51zdfkLNnjatgMgAAAACOh6IZqIiGw235zS/enpvu39Ozds7scfnXN1+QqWNHVDAZAAAAAMdL0Qz0uYf2NuVtn701D+5u6ll7yfIZ+cjV52TEMA/9AwAAABhoFM1An7p5w978+udXZd+htp61d1+2MO+5fKGH/gEAAAAMUIpmoM98+bbN+YNvrElbR9dD/2prqvK3V5+Tl58zs8LJAAAAADgZimbglOvoLPNX167Lv9y4oWdt8pjh+Zc3n58Vp02oYDIAAAAAeoOiGTilmlra8+4v3Znv3buzZ23J9Lr8+1tXZtb4kRVMBgAAAEBvUTQDp8zW/YfzK5+9Net2NPasXX7m1Pzd68/LmOH++AEAAAAYLDQ9wClx+8P78mufW5U9B1t61t5+8YJ84IVLUl3loX8AAAAAg4miGeh1/3PXtvz2V+5Ka3tnkmRYdZE/f+WyvHblnAonAwAAAOBUUDQDvaYsy3zse/fnH75/f8/a+FHD8qk3nZ9nLJhUwWQAAAAAnEqKZqBXNLd15P1fuSvXrN7es3b6lNH59FtXZu6k0RVMBgAAAMCppmgGTtquA8351c/dlru2NPSsPXfh5Hz8DSsybuSwCiYDAAAAoC8omoGTcvfWhvzq527L9obmnrU3P3Nu/vClS1NTXVXBZAAAAAD0FUUzcMKuu2dH3vOlO3O4rSNJUl1V5I9etjRvfua8ygYDAAAAoE8pmoHjVpZlPvXjDfnr69alLLvW6obX5BNvXJGLF02pbDgAAAAA+pyiGTguLe0d+f2v352v3b6lZ+20iaPy6bdekDOm1lUwGQAAAACVomgGjtnegy359c+vyq2b9vWsXTh/Yj71pvMzcXRtBZMBAAAAUEmKZuCY3LezMb/yH7dmc/3hnrWrz5+dP3/VstTWeOgfAAAAwFCmaAae1o/W78pvffGONLa0J0mKIvm9Fy7Jr128IEVRVDgdAAAAAJWmaAaeVFmW+Y+fbcqffGdtOrsf+jeqtjp///rzcsXSaZUNBwAAAEC/oWgGjqqtozN//O178vlfPNyzNnPciPzrWy7IWTPHVTAZAAAAAP2Nohl4goZDbfmNL96enzywp2ftnDnj869vPj9T60ZUMBkAAAAA/ZGiGTjCnoMted0//zwP7m7qWXvZOTPzN69ZnhHDqiuYDAAAAID+StEM9Ghqac/bPnvrESXzey5fmHdfttBD/wAAAAB4UopmIEnXTOZ3fuH2rN7SkCSpKpKPve7cvOLcWRVOBgAAAEB/V1XpAEDllWWZ//f1Nfnxfbt71v7kFWcrmQEAAAA4JopmIB+5/r58ddWWnuPfuvSMvOkZcyuYCAAAAICBRNEMQ9x//uKhfPyHD/QcX33+7LzvikUVTAQAAADAQKNohiHsunt25I++dXfP8SWLp+QvrlrmwX8AAAAAHBdFMwxRt22qz7v+6450ll3H58wel39644oMq/bHAgAAAADHp1cbpaIoziuK4p1FUYx7zNrooij+oyiK/UVRbCuK4t29+Z7A8XtgV2N+5T9uS0t7Z5Jk3qRR+fe3rsyo2poKJwMAAABgIOrtrYu/m+QPyrJseMza/5fkl7vfa1KSjxZFceWJvkFRFK8piuIfi6K4qSiKA0VRlEVRfP5Jrp3Xff7JPr50ojlgoNp5oDlv+fStaTjcliSZPKY2//G2CzN5zPAKJwMAAABgoOrt7YsXJPnhIwdFUQxL8pYktyS5JMnEJHckeVeS60/wPT6Y5JwkB5NsSbLkGO65K8k3j7J+91HWYNA60NyWt3z6lmzdfzhJMqq2Op9+68rMnTS6wskAAAAAGMh6u2iemq7y9xEXJKlL8s9lWTYn2VYUxbeSvPAk3uO93e/xQJLn5THF9lO4syzLD5/Ee8KA19Lekbd/blXW7WhMktRUFfmnN67I8tnjKxsMAAAAgAGvt4vm8nGv+ZzutR8/Zm13kikn/AZl+dgd0yf6MjCkdHaWef+X78rPN+ztWfvLVy/PJYunVjAVAAAAAINFbxfNDyd5xmOOX5FkS1mWGx6zNjPJvl5+36czsyiKt6drRvTeJD8vy3J1H2eAivmL796b76ze3nP8Oy9YnNecP7uCiQAAAAAYTHq7aP5ykj8uiuKrSZqTPDPJ3z3umjOTPNjL7/t0ruj+6FEUxY+SvKUsy4eP5QWKolj1JKeOZUY0VMy/3bQh//aTjT3Hv/yMuXnnJadXMBEAAAAAg01VL7/ex5L8PMlVSd6Qrofw/ckjJ4uimJ9kZY4cpXEqHUryp0nOTzKh++ORuc6XJPl+URSegsag9T93bcufXXNvz/ELzpqWD7/8LGNnAAAAAOhVvbqjuSzLg0meXRTF2d1La8uy7HzsJekqoW/rzfd9ijy7kvzh45ZvLIriyiQ/SXJRkv+b5O+P4bXOP9p6907nFScZFXrdzx7Yk/d/+c6e4wvmTsjfv/68VFcpmQEAAADoXb29ozlJUpbl3d0fnY9b31SW5bfKstx6Kt73WJVl2Z7k37oPL65kFjgV1m47kF/7z1Vp6yiTJGdMHZN/e8sFGTGsusLJAAAAABiMerVoLorirqIo3lEURV1vvu4psrv7s9EZDCpb9h3KWz9zSw62tCdJpo0dnv9424UZP6q2wskAAAAAGKx6e0fz0iQfT7KtKIp/LYrigl5+/d70jO7PGyqaAnrRvqbWvOXTt2RXY0uSpG54Tf7jbRdm1viRFU4GAAAAwGDW20Xz7CQfStdu4V9JcnNRFLcVRfGrlXjoXlEUK4qieMKPsSiKy5K8t/vw832bCk6N5raO/N/P3ZYHdzclSWqrq/Ivb74gS6aPrXAyAAAAAAa73n4Y4M4kf5HkL4qieEGSX0vysiSfSvK3RVF8Icm/lGV554m+R1EUr0zyyu7D6d2fn1kUxWe7v95TluVvd3/90SQLi6L4WZIt3WvLk1za/fWHyrL82Ylmgf6io7PMu/7rjqx6aF/P2kdfd06eefqkCqYCAAAAYKjo1aL5scqyvC7JdUVRTEvX7uZfSfL2JG8viuLWdJXPXyrLsvk4X/rcJG953NqC7o8keSjJI0XzfyZ5VZKVSV6UZFiSnUm+nOTjZVnedJzvDf1OWZb5w2/dnevX7uxZ+9BLl+aly2dWMBUAAAAAQ0lRluWpf5OiKJK8PF3zm2d1L5dJ6pP8eVmWf3fKQ5wiRVGsWrFixYpVq1ZVOgpD1Md/cH/+9vr7eo5/7eIF+f0Xn1nBRAAAAAAMYMWJ3NTbM5qPUBTFrKIo/ihdu4y/nq5RF/+TrtEXf5qkI8lHiqL401OZAwarL9+2+YiS+RXnzszvvXBJBRMBAAAAMBT1etFcdHlxURTfSrIxyR+la2TFXyRZUJblK8uy/J+yLD+cZGGSVekaqwEchx+u35X/9/U1PcfPPmNS/uY156Sq6oT+0QkAAAAATlivzmguiuJD6SqN56Rri/WNSf4pydfLsmx//PVlWTYWRfHtJB/uzRww2N21eX/e+fnb09HZNfrmzBlj86k3nZ/amlP6HykAAAAAwFH19sMA/zjJgXSVy58sy3LtMdyzKsnnejkHDFqb9jTlbZ+9NYfbOpIks8aPzH/8n5WpGzGswskAAAAAGKp6u2j+9SRfKMuy6VhvKMvyu0m+28s5YFDa3diSN3/6luxtak2SjB81LJ/7lQszdeyICicDAAAAYCjr1aK5LMt/6c3XAx7V1NKeX/mPW/Nw/aEkyfCaqvz7W1bm9CljKpwMAAAAgKHOQFcYANo6OvPOL9ye1VsakiRVRfLxN6zI+XMnVDgZAAAAAPT+6IwURTE6yTuTvCDJrCTDj3JZWZbl6b393jAYlWWZ3/vamvz4vt09a3/6yrNzxdJpFUwFAAAAAI/q1aK5KIrxSX6SZGm6Hgo4NklDktokI7sv25akrTffFwazv71+fb52+5ae43ddekbeeNHcCiYCAAAAgCP19uiMD6arZP6VJI/8N/0fSzImybOS3J7kwSRn9vL7wqD0nz/flE/88MGe49ddMCfvvWJRBRMBAAAAwBP1dtH88iQ3lmX5mbIsy0cWyy6/SPLiJEuS/EEvvy8MOtfevSN/+D/39Bw/f/GU/Pmrzk5RFBVMBQAAAABP1NtF85wkqx5z3JnHzGguy3JXkv9N8vpefl8YVG7dVJ93femOPPLPNefMGZ9PvHFFaqo9vxMAAACA/qe3W6tD6SqXH9GQZPrjrtmZrocEAkdx/87G/N//uC2t7V2/leZNGpVPv+WCjKrt9Wd3AgAAAECv6O2ieXO6djU/Ym2Si4uieOz7PCfJjl5+XxgUdjQ05y2fviUNh7uelzl5TG0+97aLMmnM8Ke5EwAAAAAqp7eL5h8neV7x6BDZ/05yepLvFkXxG0VRfCXJM5J8t5ffFwa8A81teetnbsm2huYkyeja6nzmrRfmtEmjKpwMAAAAAJ5ab/+3+P+RpDbJ7HTtbv5UkkuTvDLJld3X/DTJB3v5fWFAa2nvyK997ras29GYJKmpKvLJN52fZbPHVTgZAAAAADy9Xi2ay7K8Pck7HnPcnuSqoijOT3JGkk1Jbi3LsvPorwBDT2dnmfd9+a78YkN9z9pfvXp5Ll40pYKpAAAAAODY9cnTxcqyXJVkVV+8Fww0f3v9+lyzenvP8QdeuDivPn92BRMBAAAAwPE5JUVzURRzk0xJUibZXZblw6fifWCg+97anfmnHz3Yc/yWZ87NO553egUTAQAAAMDx67WHARZFMbkoio8WRbE9yYYkNye5JcnGoii2FUXxN0VRTOyt94OBbnP9obzvy3f2HD9v0ZT84cvOyqPP0gQAAACAgaFXiuaiKBYmuS3Ju5NMS9KRZFeS3d1fT0/yviS3FUWxoDfeEway5raOvOMLq3KguT1JMnPciPzd685NdZWSGQAAAICB56SL5qIoqpJ8IclpSX6c5PIkY8qynFGW5fQkdUmuTHJjknlJPn+y7wkD3Z9+Z23u3nogSTKsusgn3rgiE0bXVjgVAAAAAJyY3tjRfGWSC5J8OcllZVn+oCzL1kdOlmXZUpbl95JcmuSrSS4qiuKKXnhfGJC+ecfWfOHmR8eW/8GLz8x5p02oYCIAAAAAODm9UTS/OklLkt8qy7J8sou6z/1mkrYkr+mF94UB5/6djfl/X1/Tc/yS5TPylmfNq1wgAAAAAOgFvVE0r0jy07Isdz/dhWVZ7kryk+57YEhpamnPr39+VQ63dSRJFkwZnb969XIP/wMAAABgwOuNonlOknuO4/p7kszthfeFAaMsy/ze19fkwd1NSZIRw6ryyTeenzHDayqcDAAAAABOXm8UzWOT7D+O6/en6wGBMGR8/hcP5dt3bes5/otXLcvi6X4bAAAAADA49EbRXJuk4ziu7+y+B4aEOzfvz598Z23P8S9deFquWjG7gokAAAAAoHf1RtGcJE/6EEAYyvY1teY3vnB72jq6foucNXNs/uhlSyucCgAAAAB6V28NiP1wURQf7qXXgkGhs7PM+758Z7buP5wkqRtRk0++8fyMGFZd4WQAAAAA0Lt6q2gujvN6O6AZ9D754wfzw/W7e44/cvU5OW3SqAomAgAAAIBT46SL5rIse2v8BgwaP3twTz5y/fqe47dfvCBXnjW9gokAAAAA4NRREkMv23mgOe/6rzvS2b1v/8J5E/PbL1hc2VAAAAAAcAopmqEXtXd05re+eEf2HGxNkkweU5t/fMN5GVbttxoAAAAAg5f2C3rR31y/Prdsqk+SVBXJP/zSeZk2dkSFUwEAAADAqaVohl5y/T078s8/3tBz/P4rF+dZp0+uYCIAAAAA6BuKZugFD+89lPd/5a6e4+cvnpJ3PO/0CiYCAAAAgL6jaIaT1NzWkXd8YVUam9uTJLPGj8xHX3tuqqqKCicDAAAAgL6haIaT9MffXpt7th1IkgyrLvKJN67IhNG1FU4FAAAAAH1H0Qwn4eu3b8l/3fJwz/GHXro0584ZX7lAAAAAAFABimY4Qet3NOYPvnF3z/HLzpmZX37G3AomAgAAAIDKUDTDCTjY0p53fGFVDrd1JElOnzI6/99Vy1IU5jIDAAAAMPQomuE4lWWZ3/3a6mzY3ZQkGTmsOp980/kZM7ymwskAAAAAoDIUzXCcPvfzh3LN6u09x39x1dlZNK2ugokAAAAAoLIUzXAc7nh4X/7smrU9x2+46LS86rzZFUwEAAAAAJWnaIZjtK+pNb/xhdvT1lEmSc6eNTZ/+NKlFU4FAAAAAJWnaIZj0NlZ5j3/fWe2NTQnScaOqMkn33h+RgyrrnAyAAAAAKg8RTMcg0/88IH8+L7dPccfee25mTNxVAUTAQAAAED/oWiGp/HTB/bkY9+7r+f41593eq5YOq2CiQAAAACgf1E0w1PY0dCcd3/pjnR2jWXOhfMn5revXFTZUAAAAADQzyia4Um0dXTmt/7r9uw52JokmTxmeD7+S+elptpvGwAAAAB4LI0ZPIm/uW59bt20L0lSVST/8EvnZurYERVOBQAAAAD9j6IZjuLau3fkX27c0HP8/isX51mnT65gIgAAAADovxTN8DgP7W3K73zlrp7jS5dMzTued3oFEwEAAABA/6ZohsdobuvIOz5/expb2pMks8aPzEdfe06qqooKJwMAAACA/kvRDI/x4f+5J2u3H0iS1FZX5ZNvWpHxo2ornAoAAAAA+jdFM3T76qot+dKtm3uOP/SypVk+e3zlAgEAAADAAKFohiTrdhzIB7+5puf45efMzJsuOq2CiQAAAABg4FA0M+Q1NrflHZ+/Pc1tnUmSM6aOyf931bIUhbnMAAAAAHAsFM0MaWVZ5ve+tiYb9zQlSUYOq84n37gio4fXVDgZAAAAAAwcimaGtM/+bFOuWbO95/gvX70sC6fVVTARAAAAAAw8imaGrFUP7cufX3Nvz/GbnnFaXnHurAomAgAAAICBSdHMkFTf1Jrf/OLtae8skyTLZ4/Lh166tMKpAAAAAGBgUjQz5HR0lnn3l+7I9obmJMm4kcPyiTesyPCa6gonAwAAAICBSdHMkPPxHzyQm+7f03P80deekzkTR1UwEQAAAAAMbIpmhpT1Oxrz99+/r+f4HZecnsvOnFbBRAAAAAAw8CmaGVL+5rp16R7LnJXzJuT9VyyqbCAAAAAAGAQUzQwZt2ysz/fu3dVz/OGXn5Waar8FAAAAAOBkadkYEsqyzF/+7709x688d2bOmjmugokAAAAAYPBQNDMkXL92Z25/eH+SZFh1kfdfubiygQAAAABgEFE0M+i1d3Tmb65b33P8pmfMzZyJoyqYCAAAAAAGF0Uzg97Xbt+SB3YdTJKMGV6T33z+GRVOBAAAAACDi6KZQe1wa0c+dsP9Pcdvv3hBJo0ZXsFEAAAAADD4KJoZ1D77s03ZcaA5STJ5zPD8ynPnVzgRAAAAAAw+imYGrf2HWvNPP3qg5/g9ly/MqNqaCiYCAAAAgMFJ0cyg9U8/ejCNze1JkvmTR+d1K+dUOBEAAAAADE6KZgalrfsP57M/29Rz/DsvWJxh1X65AwAAAMCpoHljUPrYDfeltb0zSXLOnPF50dnTK5wIAAAAAAYvRTODzrodB/K127f0HP/eC5ekKIoKJgIAAACAwU3RzKDzN9euT1l2fX3J4il55umTKhsIAAAAAAY5RTODys0b9ub763YlSYoi+cALllQ4EQAAAAAMfopmBo2yLPOX167rOX7VubOydObYCiYCAAAAgKFB0cygcd09O3PHw/uTJLXVVXnvFYsqGwgAAAAAhghFM4NCe0dn/ua6R3czv+kZczNn4qgKJgIAAACAoUPRzKDw1VVb8uDupiTJmOE1+c1Lz6hwIgAAAAAYOhTNDHiHWzvyse/d13P8689bkImjayuYCAAAAACGFkUzA95nfrYxOw+0JEmm1A3P254zv8KJAAAAAGBoUTQzoO1ras0nf/Rgz/F7Ll+YUbU1FUwEAAAAAEOPopkB7Z9+9EAam9uTJAsmj85rL5hT4UQAAAAAMPQomhmwtuw7lP/42UM9x7/zgsUZVu2XNAAAAAD0Na0cA9bHbrg/rR2dSZJz5ozPC8+eXuFEAAAAADA0KZoZkNbtOJCv37Gl5/j/vWhJiqKoYCIAAAAAGLoUzQxIf33t+pRl19fPXzwlz1gwqbKBAAAAAGAIUzQz4Pxiw978YN2uJElRJB944ZIKJwIAAACAoU3RzIBSlmX+8n/X9Ry/6rxZOXPG2AomAgAAAAAUzQwo192zM3du3p8kqa2uyvuuWFTZQAAAAACAopmBo72jM3993aO7md/8zLmZPWFUBRMBAAAAAImimQHkK6u2ZMPupiRJ3fCa/Mbzz6hwIgAAAAAgUTQzQBxu7cjHbriv5/jXLzk9E0bXVjARAAAAAPAIRTMDwqd/ujG7GluSJFPrhuf/PHteZQMBAAAAAD0UzfR7+5pa86kfPdhz/J7LF2VUbU0FEwEAAAAAj6Vopt/7xA8fSGNLe5JkweTRee0FsyucCAAAAAB4LEUz/dqWfYfyuZ8/1HP8gRcuTk21X7YAAAAA0J9o7OjXPnrDfWnt6EySnDtnfF5w1vQKJwIAAAAAHk/RTL917/YD+cYdW3uOf+9FS1IURQUTAQAAAABHo2im3/rra9elLLu+vnTJ1DxjwaTKBgIAAAAAjkrRTL/08wf35ofrdydJiqJrNjMAAAAA0D8pmul3yrLMX167ruf4qvNmZ8n0sRVMBAAAAAA8FUUz/c519+zIXZv3J0lqa6ryvisXVTYQAAAAAPCUFM30K+0dnfnra9f3HL/lmXMza/zICiYCAAAAAJ6Oopl+5cu3bcmGPU1JkroRNXnnJWdUOBEAAAAA8HQUzfQbh1rb83ffu6/n+B2XnJ4Jo2srmAgAAAAAOBaKZvqNz/x0U3Y1tiRJpo0dnv/zrPkVTgQAAAAAHAtFM/1CfVNrPvWjB3uO33P5ooysra5gIgAAAADgWCma6Rc+8cMH0tjSniRZMGV0rj5/doUTAQAAAADHStFMxW2uP5T//PlDPccfeMGS1FT7pQkAAAAAA4U2j4r72A33pbWjM0ly3mnj84KzplU4EQAAAABwPBTNVNTabQfyjTu39hz/3guXpCiKCiYCAAAAAI6XopmK+uvr1qUsu76+bMnUXLRgUmUDAQAAAADHTdFMxfz8wb350frdSZKiSD7wwiUVTgQAAAAAnAhFMxVRlmX+8tp1PcevXjE7i6fXVTARAAAAAHCiFM1UxLV378hdm/cnSWprqvLeKxZVNhAAAAAAcMIUzfS5to7O/M1163uO3/qseZk1fmQFEwEAAAAAJ0PRTJ/78m2bs2FPU5KkbkRN3nnJ6RVOBAAAAACcDEUzfepQa3v+7nv39xy/85IzMn5UbQUTAQAAAAAnS9FMn/r0TzZmd2NLkmT62BH5P8+eV9lAAAAAAMBJUzTTZ+qbWvOpH2/oOX7vFQszYlh1BRMBAAAAAL1B0Uyf+fgPHsjBlvYkyelTRufVK2ZXOBEAAAAA0BsUzfSJzfWH8p+/2NRz/IEXLklNtV9+AAAAADAYaProEx+94b60dZRJkhWnjc+VS6dVOBEAAAAA0FsUzZxy92xryDfv3Npz/HsvOjNFUVQwEQAAAADQmxTNnHJ/fe36lF2bmXP5mVNz4fyJlQ0EAAAAAPQqRTOn1M8e3JMf37c7SVJVJL/zgiUVTgQAAAAA9DZFM6fUp3+ysefrV6+YncXT6yqYBgAAAAA4FRTNnDK7Gpvzw/W7e47f+fwzKpgGAAAAADhVBlzRXBTFa4qi+MeiKG4qiuJAURRlURSff5p7nlUUxXeLoqgviuJwURSri6J4T1EU1X2Veyj6xu1b09HZNZz5wnkTM3/y6AonAgAAAABOhZpKBzgBH0xyTpKDSbYkecqhv0VRvCLJ15I0J/nvJPVJXpbkY0meneTqUxl2qCrLMl9ZtaXn+OoLZlcwDQAAAABwKg24Hc1J3ptkUZKxSd7xVBcWRTE2yb8m6UhySVmWv1KW5e8kOTfJz5O8piiK15/auEPTnZv354FdB5Mko2qr8+JlMyqcCAAAAAA4VQZc0VyW5Q/Lsry/LMvyGC5/TZIpSb5UluVtj3mN5nTtjE6epqzmxHz5tkd3M79k2YyMHj4QN88DAAAAAMdiwBXNx+nS7s/XHuXcjUkOJXlWURTD+y7S4He4tSPfuWtbz/HVF8ypYBoAAAAA4FQb7NtMF3d/vu/xJ8qybC+KYmOSs5IsSHLvU71QURSrnuTUU86IHoquu2dHGlvakyTzJo3KynkTKpwIAAAAADiVBvuO5nHdnxue5Pwj6+NPfZSh4yurNvd8ffUFc1IURQXTAAAAAACn2mDf0dxryrI8/2jr3TudV/RxnH5rc/2h/PSBvUmSqiK5asWsCicCAAAAAE61wb6j+ZEdy+Oe5Pwj6/tPfZSh4Wu3P/oQwOcsnJIZ40ZWMA0AAAAA0BcGe9G8vvvzosefKIqiJsn8JO1JNvRlqMGqs7PMV1c9WjS/9oLZFUwDAAAAAPSVwV40/6D78wuPcu7iJKOS/Kwsy5a+izR4/WLj3mzZdzhJMm7ksFx+5rQKJwIAAAAA+sJgL5q/mmRPktcXRXHBI4tFUYxI8mfdh5+sRLDB6Cu3Pbqb+RXnzsyIYdUVTAMAAAAA9JUB9zDAoihemeSV3YfTuz8/syiKz3Z/vacsy99OkrIsDxRF8avpKpx/VBTFl5LUJ3l5ksXd6//dN8kHtwPNbfnfu7f3HF99/pwKpgEAAAAA+tKAK5qTnJvkLY9bW9D9kSQPJfntR06UZfnNoiiel+QPkrw6yYgkDyR5X5J/KMuyPNWBh4JrVm9Pc1tnkmTJ9LqcPWtshRMBAAAAAH1lwBXNZVl+OMmHj/OenyZ58anIQ5ev3La55+urL5iToigqmAYAAAAA6EuDfUYzfeCBXY25/eH9SZKaqiKvPHdmZQMBAAAAAH1K0cxJ+8qqRx8CeNmZUzNpzPAKpgEAAAAA+pqimZPS3tGZr9++tefYQwABAAAAYOhRNHNSbrx/d3Y3tiRJptQNzyWLp1Q4EQAAAADQ1xTNnJQv3/ro2IyrzpuVmmq/pAAAAABgqNEKcsLqm1rz/XU7e46vvmB2BdMAAAAAAJWiaOaEffOOrWnrKJMk584ZnzOm1lU4EQAAAABQCYpmTkhZlvnybZt7jl97gYcAAgAAAMBQpWjmhNyz7UDW7WhMkowYVpWXnjOjwokAAAAAgEpRNHNCvvKY3cwvPGt6xo4YVsE0AAAAAEAlKZo5bi3tHfnWXdt6jq82NgMAAAAAhjRFM8fte2t3Zf+htiTJrPEj88wFkyqcCAAAAACoJEUzx+0rqx4dm/Ga82enqqqoYBoAAAAAoNIUzRyXHQ3NufG+3T3Hrzl/dgXTAAAAAAD9gaKZ4/K127eks+z6+pkLJmXOxFGVDQQAAAAAVJyimWNWlmW+umpLz/FrV9rNDAAAAAAomjkOqx7al417mpIkdcNr8sKzZlQ4EQAAAADQHyiaOWZfvu3RhwC+9JwZGVlbXcE0AAAAAEB/oWjmmBxqbc81q7f3HL/m/DkVTAMAAAAA9CeKZo7Jd9fsSFNrR5Lk9Cmjs+K08ZUNBAAAAAD0G4pmjslXHjM24+oL5qQoigqmAQAAAAD6E0UzT+uhvU25eWN9kqS6qshV582qcCIAAAAAoD9RNPO0vrpqS8/Xz1s0JVPHjqhgGgAAAACgv1E085Q6Ost87TFF82svmF3BNAAAAABAf6Ro5in97ME92dbQnCSZOLo2ly6ZVuFEAAAAAEB/o2jmKX35tkd3M7/i3JmprfFLBgAAAAA4ktaQJ9VwqC3X3bOj5/jq8+dUMA0AAAAA0F8pmnlS/7N6W1rbO5MkZ88am6Uzx1Y4Ef9/e/ceZOd51wf8+2h1v1iyZVmKL0HY8a4g5GYbMAklcTIYCgUcsNN2yqUMzDSFTgoDM+2U0sS0TOkMA4FAA9MO0BpmqO0Q0kIImSYmCZjL2I7jmkQr2/HdlmTJlqyVrMvuPv1jj1ZHslbWanfPe95zPp8ZzXve9+yufrGeeaL97qPfDwAAAAD6kaCZOd1939Ozr51mBgAAAADmImjmrMZ3H8qXnjmYJFk5sizf99bLG64IAAAAAOhXgmbO6q6u08zf/vVbs2ntygarAQAAAAD6maCZVzkxNZ0/fvDZ2fvbbriywWoAAAAAgH4naOZV7tm5N/smjidJtl20Ov/g2i0NVwQAAAAA9DNBM69y533PzL7+/uuuyMiy0mA1AAAAAEC/EzRzmhcOHcs943tn72+9XtsMAAAAAODcBM2c5o+/+GympmuS5Bu3X5yrt6xvuCIAAAAAoN8JmplVa82d9z09e3/b9Vc1WA0AAAAA0BaCZmZ96ZmDeWTvRJJkzYqRfNebX9dwRQAAAABAGwiamXVX12nm73rT67J+1fIGqwEAAAAA2kLQTJLk6Imp/O8vPTd7/74bDAEEAAAAAM6PoJkkyZ///e4cOjqZJPmazWvzTV97ScMVAQAAAABtIWgmSXLXfc/Mvr71uitTSmmwGgAAAACgTQTN5NkDr+SvHtuXJCkl+YHrtc0AAAAAAM6foJl87P5nUuvM6299w6W5fNOaZgsCAAAAAFpF0Dzkpqdr7r7/VNuM2264qsFqAAAAAIA2EjQPub99/MU89eKRJMlFq5fn5q/f2nBFAAAAAEDbCJqH3F33Pz37+nvfenlWrxhpsBoAAAAAoI0EzUNs4thk/uz/7Z69v+16bTMAAAAAgPkTNA+xP33oubxyYipJMrZ1Q9585caGKwIAAAAA2kjQPMTuvK97COCVKaU0WA0AAAAA0FaC5iH12AsTuf/Jl5Iky5eV3PK2KxquCAAAAABoK0HzkLr7/lOnmW/acVkuXb+qwWoAAAAAgDYTNA+hyanp/NEDp4Lm991gCCAAAAAAcOEEzUPoC4/sy56XjyVJLl2/Mu8a29JwRQAAAABAmwmah9Bd9z89+/q9b7siK0YsAwAAAADgwkkYh8xLh4/n/3557+z9bdpmAAAAAAALJGgeMp948Nkcn5pOkrzlqk0Z3bqh4YoAAAAAgLYTNA+Zu+4/NQTwtuuvbLASAAAAAGBQCJqHyN8/dzB//9zLSZJVy5fle95yecMVAQAAAACDQNA8RO6679Rp5u9447ZsXLOiwWoAAAAAgEEhaB4Sxyen84kHn529f58hgAAAAADAIhE0D4nPfGVPXjpyIklyxaY1efs1mxuuCAAAAAAYFILmIXHnfU/Pvv6B667IsmWlwWoAAAAAgEEiaB4Ce14+ms/temH2/tbrtc0AAAAAABaPoHkI/NEDz2a6zry+8epL8vrNa5stCAAAAAAYKILmAVdrzV33n2qbcZvTzAAAAADAIhM0D7gHnnopX33hcJJk/arl+Ydv2tZwRQAAAADAoBE0D7i77ntm9vV3v+l1WbtyeYPVAAAAAACDSNA8wI4cn8yfPPT87P37vvHKBqsBAAAAAAaVoHmAferh3Zk4NpkkuXrLulz3+osbrggAAAAAGESC5gF2532nhgDeev2VKaU0WA0AAAAAMKgEzQPqqf1H8jdffTFJsqwkP3CdthkAAAAAwNIQNA+oux84NQTwnaNbsvWi1Q1WAwAAAAAMMkHzAJqervnY/aeC5ttuuKrBagAAAACAQSdoHkD3PrY/zx54JUmyae2KvOfrLmu4IgAAAABgkAmaB9Bd958aAnjLW6/IquUjDVYDAAAAAAw6QfOAOfjKiXzq4d2z97fdYAggAAAAALC0BM0D5k8eei7HJqeTJF//uovyxss3NlwRAAAAADDoBM0D5s77uocAOs0MAAAAACw9QfMAOXDkeL709IEkyYqRku976xXNFgQAAAAADAVB8wAZ331o9vXo1g25ZN3KBqsBAAAAAIaFoHmA7NpzKmge27ahwUoAAAAAgGEiaB4gO7tONI9tFTQDAAAAAL0haB4g3SeaR51oBgAAAAB6RNA8IGqtp/Vo3iFoBgAAAAB6RNA8IHa/fDQvH51MkmxYvTzbLlrdcEUAAAAAwLAQNA+IM08zl1IarAYAAAAAGCaC5gHRHTSPGgQIAAAAAPSQoHlAjHcNAhzTnxkAAAAA6CFB84DY1R00O9EMAAAAAPSQoHkATE3XPLJnYvbeiWYAAAAAoJcEzQPgyf2Hc2xyOkmy9aJV2bR2ZcMVAQAAAADDRNA8ALrbZhgECAAAAAD0mqB5AOzcrT8zAAAAANAcQfMAOG0QoP7MAAAAAECPCZoHwGknmgXNAAAAAECPCZpb7uiJqTyx73CSpJTk2ssEzQAAAABAbwmaW+6xFyYyXWdef80la7Nm5UizBQEAAAAAQ0fQ3HLjXW0zRg0CBAAAAAAaIGhuufGuQYA79GcGAAAAABogaG650040C5oBAAAAgAYImltuV1fQPKZ1BgAAAADQAEFzi7189ESeO3g0SbJyZFm2X7qu4YoAAAAAgGEkaG6x7tPMV29ZlxUj/jgBAAAAgN6TTLaYQYAAAAAAQD8QNLeYQYAAAAAAQD8QNLfYuEGAAAAAAEAfEDS3VK01u7paZ4w50QwAAAAANETQ3FIvHDqWl46cSJKsX7U8V2xa03BFAAAAAMCwEjS3VPcgwNGt61NKabAaAAAAAGCYCZpb6rT+zNpmAAAAAAANEjS3VHfQPGoQIAAAAADQIEFzSxkECAAAAAD0C0FzC01P1+zaMzF7P+ZEMwAAAADQIEFzCz390pG8cmIqSXLp+lXZvH5VwxUBAAAAAMNM0NxCpw8CXN9gJQAAAAAAguZWMggQAAAAAOgnguYWGu8aBLjDIEAAAAAAoGGC5hZyohkAAAAA6CeC5pY5NjmVx/cdnr0XNAMAAAAATRM0t8zj+w5ncromSa66ZE3WrVrecEUAAAAAwLATNLdMd9uMMaeZAQAAAIA+IGhumdOCZoMAAQAAAIA+IGhuGYMAAQAAAIB+I2humfE9TjQDAAAAAP1F0NwiE8cm88xLryRJli8rufrS9Q1XBAAAAAAgaG6VXV2nma/esi4rl/vjAwAAAACaJ6lskV2nDQK8qMFKAAAAAABOETS3yM7uoHmrthkAAAAAQH8QNLdId+uM0a0GAQIAAAAA/UHQ3CLdQfMOrTMAAAAAgD4haG6JfRPHsm/ieJJk7cqRXHnxmoYrAgAAAACYIWhuie5BgNdu3ZBly0qD1QAAAAAAnCJobgmDAAEAAACAfjUUQXMp5YlSSp3j1+6m6zsfBgECAAAAAP1qedMF9NDBJB8+y/OJHtdxQcYNAgQAAAAA+tQwBc0Haq0farqICzE9XU/r0Ty6TesMAAAAAKB/DEXrjLZ79sArOXx8KklyybqV2bJ+VcMVAQAAAACcMkwnmleVUn4wyeuTHE7yUJLP11qnmi3rtZ3en3l9SikNVgMAAAAAcLphCpq3JbnjjGePl1J+tNb6udf65FLK/XO8tWPBlb2GnV1tM8YMAgQAAAAA+sywtM743STvyUzYvC7Jm5L8dpLtSf6slPKW5kp7bd0nmscMAgQAAAAA+sxQnGiutd5+xqOHk7y/lDKR5GeSfCjJe1/ja1x/tuedk87XLUKZcxrvPtFsECAAAAAA0GeG5UTzXH6rc/22Rqs4hxNT03nshYnZ+1GtMwAAAACAPjPsQfMLneu6Rqs4hyf2Hc6JqZokuWLTmmxYvaLhigAAAAAATjfsQfONnetXG63iHLoHAY5u1TYDAAAAAOg/Ax80l1K+rpTyqhPLpZTtSX6jc/v7PS1qHgwCBAAAAAD63TAMA/zHSX6mlPL5JE8mOZTkmiTfnWR1kk8m+eXmyju3nQYBAgAAAAB9bhiC5nuSjCV5W5J3ZKYf84Ekf5nkjiR31FprY9W9htNONG91ohkAAAAA6D8DHzTXWj+X5HNN13EhjhyfzFMvHkmSjCwruXpL384sBAAAAACG2MD3aG6zR/ZM5ORZ6+2b12b1ipFmCwIAAAAAOAtBcx8b72qbscMgQAAAAACgTwma+9h41yDA0a0bGqwEAAAAAGBuguY+dtogwG3rG6wEAAAAAGBuguY+1n2ieUzrDAAAAACgTwma+9RLh49n76FjSZLVK5bl9ZesbbgiAAAAAICzEzT3qe5BgNdetiEjy0qD1QAAAAAAzE3Q3KcMAgQAAAAA2kLQ3KfGDQIEAAAAAFpC0NyndhkECAAAAAC0hKC5D9VaTz/RrHUGAAAAANDHBM196PmDR3Po6GSSZOOaFdl60aqGKwIAAAAAmJuguQ+deZq5lNJgNQAAAAAA5yZo7kPjXf2ZRw0CBAAAAAD6nKC5DxkECAAAAAC0iaC5D+3cbRAgAAAAANAeguY+Mzk1nUdfmJi9FzQDAAAAAP1O0NxnnnzxSI5PTidJtl20OhvXrmi4IgAAAACAcxM095nTBwE6zQwAAAAA9D9Bc5/pDpp3CJoBAAAAgBYQNPeZ0040688MAAAAALSAoLnP7NrjRDMAAAAA0C6C5j5y9MRUnth/OElSSvKGy9Y3XBEAAAAAwGsTNPeRR/dOZLrOvN6+eV1WrxhptiAAAAAAgPMgaO4j3f2Zx/RnBgAAAABaQtDcR8a7+jOP6s8MAAAAALSEoLmPONEMAAAAALSRoLmP7Oo60TzmRDMAAAAA0BKC5j5x8MiJPH/waJJk5fJl2b55bcMVAQAAAACcH0Fzn9i199Rp5jdsWZ/lI/5oAAAAAIB2kGb2iZ27tc0AAAAAANpJ0NwndnUFzaMGAQIAAAAALSJo7hPjXYMAdzjRDAAAAAC0iKC5D9RaM959olnQDAAAAAC0iKC5D+w9dCwHXzmRJNmwanku37i64YoAAAAAAM6foLkPnHmauZTSYDUAAAAAAPMjaO4D4wYBAgAAAAAtJmjuAwYBAgAAAABtJmjuA040AwAAAABtJmhu2NR0zSN7TwXNY040AwAAAAAtI2hu2NMvHsnRE9NJki0bVuWSdSsbrggAAAAAYH4EzQ3b2dU2Y0zbDAAAAACghQTNDdu1R9sMAAAAAKDdBM0NG3eiGQAAAABoOUFzw8adaAYAAAAAWk7Q3KBjk1N5fN/h2ftrt65vsBoAAAAAgAsjaG7QY3sPZ2q6Jklef8narF25vOGKAAAAAADmT9DcIIMAAQAAAIBBIGhu0E6DAAEAAACAASBoblD3ieZRJ5oBAAAAgJYSNDdovOtE8w5BMwAAAADQUoLmhhw6eiLPHnglSbJipORrL13XcEUAAAAAABdG0NyQXXsmZl9fs2V9Voz4owAAAAAA2km62ZDuthmjBgECAAAAAC0maG5I9yDAMf2ZAQAAAIAWEzQ3pPtE85gTzQAAAABAiwmaG1BrzbgTzQAAAADAgBA0N2DfxPG8ePh4kmTdypFcsWlNwxUBAAAAAFw4QXMDuvszX7t1Q5YtKw1WAwAAAACwMILmBuzUnxkAAAAAGCCC5gbs2q0/MwAAAAAwOATNDdhpECAAAAAAMEAEzT02PV3ziKAZAAAAABggguYee/bAKzlyfCpJsnndyly6flXDFQEAAAAALIyguce6BwGOGgQIAAAAAAwAQXOP7dI2AwAAAAAYMILmHus+0SxoBgAAAAAGgaC5x3YJmgEAAACAASNo7qHjk9N57IWJ2ftrL1vfYDUAAAAAAItD0NxDj+87nMnpmiS5YtOabFi9ouGKAAAAAAAWTtDcQ+NdgwB3aJsBAAAAAAwIQXMPje9+efb1qKAZAAAAABgQguYeGt99qj/z2FZBMwAAAAAwGATNPbSrq3XGmBPNAAAAAMCAEDT3yOFjk3nqxSNJkpFlJVdvWddwRQAAAAAAi0PQ3COP7D3VNuPqS9dl1fKRBqsBAAAAAFg8guYeMQgQAAAAABhUguYeMQgQAAAAABhUguYeMQgQAAAAABhUguYe2bm7K2h2ohkAAAAAGCCC5h7YP3Es+yaOJUlWr1iW11+ytuGKAAAAAAAWj6C5B3btOdWfeXTrhixbVhqsBgAAAABgcQmae2B898uzr0e1zQAAAAAABoyguQfGu0407zAIEAAAAAAYMILmHnCiGQAAAAAYZILmJVZrPa1HsxPNAAAAAMCgETQvsecOHs3Esckkyaa1K7Jlw6qGKwIAAAAAWFyC5iV2ZtuMUkqD1QAAAAAALD5B8xIb361tBgAAAAAw2ATNS8wgQAAAAABg0Amal9i4QYAAAAAAwIATNC+hyanpPLb3VNB8rRPNAAAAAMAAEjQvoSf2H87xqekkyes2rs7GNSsarggAAAAAYPEJmpdQ9yDAMW0zAAAAAIABJWheQt2DAMe0zQAAAAAABpSgeQmN7zk0+3pU0AwAAAAADChB8xLatUfrDAAAAABg8Amal8grx6fyxP7DSZJlJXnDZesbrggAAAAAYGkImpfIo3snUuvM6+2XrsvqFSPNFgQAAAAAsEQEzUtkp0GAAAAAAMCQEDQvkV0GAQIAAAAAQ0LQvETGuwYB7jAIEAAAAAAYYILmJTLe1TpjVNAMAAAAAAwwQfMSOHDkePa8fCxJsnL5smzfvK7higAAAAAAlo6geQns6mqbce1l6zOyrDRYDQAAAADA0hI0L4HuthljBgECAAAAAANO0LwExvccmn09pj8zAAAAADDgBM1LYHz3qaDZIEAAAAAAYNAJmhdZrfW0oHmHoBkAAAAAGHCC5kW25+VjefnoZJJkw+rl2XbR6oYrAgAAAABYWoLmRbbzjEGApZQGqwEAAAAAWHqC5kW2yyBAAAAAAGDICJoX2c7dgmYAAAAAYLgImhfZaSeatwqaAQAAAIDBJ2heRFPTNY/smZi9HxU0AwAAAABDQNC8iJ7cfzjHJqeTJJdtWJWL161suCIAAAAAgKUnaF5EBgECAAAAAMNI0LyIThsEqG0GAAAAADAkBM2LqPtE86gTzQAAAADAkBA0L6LxrhPNOwTNAAAAAMCQEDQvkqMnpvLE/iNJklKSay8TNAMAAAAAw0HQvEgee2EiU9M1SfI1l6zNmpUjDVcEAAAAANAbguZF0t02Y9QgQAAAAABgiAiaF8l41yDAMf2ZAQAAAIAhImheJLt2C5oBAAAAgOEkaF4k3a0zxrTOAAAAAACGiKB5EUzVmucOHk2SrBxZlu2Xrmu4IgAAAACA3hE0L4JjJ6ZnX1+9ZV1WjPjPCgAAAAAMD4noIjh6Ymr2tf7MAAAAAMCwETQvgmOCZgAAAABgiAmaF8HRyVOtMwwCBAAAAACGjaB5EWidAQAAAAAMM0HzIpiarkmSdStHcsWmNQ1XAwAAAADQW4LmRTS6bUNKKU2XAQAAAADQU4LmRbRD2wwAAAAAYAgNTdBcSrmylPI7pZTnSinHSilPlFI+XEq5eLF+j1GDAAEAAACAIbS86QJ6oZRyTZJ7k1yW5BNJdib5piT/Osl3llLeUWvdv9DfxyBAAAAAAGAYDcuJ5v+amZD5A7XWW2qt/7bW+u4kv5pkLMkvLsZvMuZEMwAAAAAwhAY+aO6cZr45yRNJfvOMtz+Y5HCSHyqlrFvI73Pp+pXZvH7VQr4EAAAAAEArDXzQnOSmzvXTtdbp7jdqrYeS/FWStUluXMhvom0GAAAAADCshqFH81jnumuO9x/JzInn0SSfmeuLlFLun+Ott5zY/3Tu+aUfy/UfXX3hVQIAAAAANOyBBx74g1rrP5vv5w1D0Lyxcz04x/snn2+6wK+/rE4en3r+sS9/6fkL/ALQh3Z0rjsbrQIWl3XNoLK2GUTWNYPIumYQWdcMIuv6Ag1D0Lwoaq3Xn+35yZPOc70PbWRdM4isawaVtc0gsq4ZRNY1g8i6ZhBZ1xduGHo0nzyxvHGO908+P7D0pQAAAAAADJ5hCJrHO9fROd6/tnOdq4czAAAAAADnMAxB8z2d682llNP+95ZSNiR5R5IjSf6m14UBAAAAAAyCgQ+aa62PJfl0ku1JfvKMt29Psi7JHbXWwz0uDQAAAABgIAzLMMCfSHJvkl8vpbwnyVeSfHOSmzLTMuPnGqwNAAAAAKDVSq216Rp6opRyVZJfSPKdSTYneT7Jx5PcXmt9qcnaAAAAAADabGiCZgAAAAAAlsbA92gGAAAAAGBpCZoBAAAAAFgQQTMAAAAAAAsiaAYAAAAAYEEEzQAAAAAALIigGQAAAACABRE0X6BSypWllN8ppTxXSjlWSnmilPLhUsrFTdcGF6Kzhuscv3Y3XR+cSynl1lLKR0opXyilvNxZt7//Gp/z9lLKJ0spL5ZSXimlPFRK+alSykiv6oZzmc+6LqVsP8ceXkspf9jr+uFsSimbSyk/Xkr5eCnl0c7+e7CU8pellB8rpZz1+xN7Nv1svuvank1blFL+SynlM6WUpzvr+sVSyhdLKR8spWye43Ps1/S1+axr+/X8lVpr0zW0TinlmiT3JrksySeS7EzyTUluSjKe5B211v3NVQjzV0p5IsmmJB8+y9sTtdZf7mU9MB+llAeTvCXJRJJnkuxI8ge11h+c4+O/L8nHkhxN8r+SvJjke5KMJbm71npbD8qGc5rPui6lbE/yeJIvJfnjs3y5h2utdy9VrXC+SinvT/LRJM8nuSfJU0m2Jvn+JBszszffVru+SbFn0+/mu67t2bRFKeV4kgeSfDnJ3iTrktyY5IYkzyW5sdb6dNfH26/pe/NZ1/br+RM0X4BSyp8nuTnJB2qtH+l6/itJfjrJb9da399UfXAhOkFzaq3bm60E5q+UclNmgrhHk7wzM9/kzRXIXdT5uI2Z+cHgfZ3nq5N8Nsm3JPmntVY/naZR81zX2zPzl+D/UWv95z0sE+allPLuzHxD96e11umu59uS/F2Sq5LcWmv9WOe5PZu+dwHrenvs2bRAKWV1rfXoWZ7/YpJ/l+Sjtdaf6DyzX9MK81zX22O/nhetM+apc5r55iRPJPnNM97+YJLDSX6olLKux6UBDK1a6z211kfq+f309NYkW5L84cm/AHe+xtEk/75z+y+XoEyYl3mua2iFWutna63/pzuM6zzfneS3Orfv6nrLnk3fu4B1Da1wtjCu487O9dquZ/ZrWmGe65p5Wt50AS10U+f66bP8ReJQKeWvMhNE35jkM70uDhZoVSnlB5O8PjM/NHkoyedrrVPNlgWL6t2d66fO8t7nkxxJ8vZSyqpa67HelQWL4vJSyr9IsjnJ/iR/XWt9qOGa4Hyd6Fwnu57Zs2m7s63rk+zZtNX3dK7d69V+TdudbV2fZL8+T4Lm+RvrXHfN8f4jmQmaRyNopn22JbnjjGePl1J+tNb6uSYKgiUw5z5ea50spTye5I1Jrk7ylV4WBovg2zu/ZpVS/iLJj9Ran2qkIjgPpZTlSX64c9sdUtizaa1zrOuT7Nm0QinlZ5Osz0xbjBuSfGtmwrhf6vow+zWtcp7r+iT79XnSOmP+NnauB+d4/+TzTUtfCiyq303ynsyEzeuSvCnJbyfZnuTPSilvaa40WFT2cQbRkST/Mcn1SS7u/DrZ1/ldST6jrRd97peSfEOST9Za/7zruT2bNptrXduzaZufzUyr0J/KTBj3qSQ311pf6PoY+zVtcz7r2n49T4JmIElSa729019uT631SK314c5Qy19JsibJh5qtEIC51Fr31lr/Q631gVrrgc6vz2fmX1n9bZI3JPnxZquEsyulfCDJzyTZmeSHGi4HFsW51rU9m7aptW6rtZbMHEr6/sycSv5iKeW6ZiuDC3c+69p+PX+C5vk7+VO4jXO8f/L5gaUvBXri5ACTb2u0Clg89nGGRq11Msl/79zax+k7pZR/leTXknw5yU211hfP+BB7Nq1zHuv6rOzZ9LvOoaSPZyZk25zkf3a9bb+mlV5jXc/1OfbrOQia52+8cx2d4/2T0ynn6uEMbXPyn4345yAMijn38U4vxa/NzMCer/ayKFhC9nH6Uinlp5J8JMnDmQnjdp/lw+zZtMp5rutzsWfT92qtT2bmBylvLKVc2nlsv6bV5ljX52K/PgtB8/zd07neXEo57b9fKWVDkndkpofL3/S6MFgiN3au/kLAoPhs5/qdZ3nv25KsTXKvadgMEPs4faeU8m+S/GqSBzMTxu2d40Pt2bTGPNb1udizaYvLO9epztV+zSA4c12fi/36LATN81RrfSzJpzMzIO0nz3j79sz8JOOOWuvhHpcGF6yU8nVna2BfStme5Dc6t7/f06Jg6dydZF+Sf1JKueHkw1LK6iT/qXP70SYKgwtVSrnuzB+Ad56/J8lPd27t4/SFUsrPZ2ZI2v1J3lNr3XeOD7dn0wrzWdf2bNqglDJaSnlVG4xSyrJSyi8muSwzwfFLnbfs1/S9+a5r+/X8lVpr0zW0TinlmiT3ZmYBfiLJV5J8c5KbMtMy4+211v3NVQjzU0r5UGaGlXw+yZNJDiW5Jsl3J1md5JNJ3ltrPd5UjXAupZRbktzSud2W5Dsy85PlL3Se7au1/uwZH393kqNJ/jDJi0m+N8lY5/n7qv+DpGHzWdellL/ITPuue5M803n/zUne3Xn987XWk9/kQWNKKT+S5Pcyc1LoIznV07PbE7XW3+v6nFtiz6aPzXdd27Npg04bmP+c5C+TPJ5kf5KtSd6ZmaFpuzPzQ5Uvd33OLbFf08fmu67t1/MnaL5ApZSrkvxCZv5ZyOYkzyf5eJLbu36iB61QSnlnkvcneVtmwox1mRnS8GCSOzJzSt9mQd/q/LDkg+f4kCdrrdvP+Jx3JPm5JN+SmR+oPJrkd5L8eq31fP6pFCyp+azrUsqPJXlvkm9IcmmSFUn2JPnrJL9Ra/3CXF8Eeuk81nWSfK7W+q4zPs+eTd+a77q2Z9MGpZRvyMz3iN+a5Mokm5Iczszhuj/NzP77qkGX9mv62XzXtf16/gTNAAAAAAAsiB7NAAAAAAAsiKAZAAAAAIAFETQDAAAAALAggmYAAAAAABZE0AwAAAAAwIIImgEAAAAAWBBBMwAAAAAACyJoBgAAAABgQQTNAAAAAAAsiKAZAAAAAIAFETQDAAAAALAggmYAAAAAABZE0AwAAAAAwIIImgEAoEVKKT9cSqmllA80XQsAAJwkaAYAgHa5oXO9r9EqAACgS6m1Nl0DAABwnkopVyTZkOTRWutk0/UAAEAiaAYAAAAAYIG0zgAAgJYopbyt05/5vzVdCwAAdBM0AwBAe1zXuT7QaBUAAHAGQTMAALSHoBkAgL4kaAYAgPa4LslUkoeaLgQAALoJmgEAoAVKKSNJ3pzkK7XWV5quBwAAugmaAQCgHXYkWRttMwAA6EOCZgAAaAf9mQEA6FuCZgAAaAdBMwAAfUvQDAAA7XBdkprkwYbrAACAVym11qZrAAAAzqGUUpIcSLK71jrWcDkAAPAqTjQDAED/uybJRUm+2HQhAABwNoJmAADof9/cuf5do1UAAMAcBM0AAND/butcP9toFQAAMIflTRcAAAC8WillWZLbk7w1yT9K8oVa64NN1gQAAHMxDBAAAPpQKWVHkoeTPJ/kE0l+vtb6UrNVAQDA2QmaAQAAAABYED2aAQAAAABYEEEzAAAAAAALImgGAAAAAGBBBM0AAAAAACyIoBkAAAAAgAURNAMAAAAAsCCCZgAAAAAAFkTQDAAAAADAggiaAQAAAABYEEEzAAAAAAALImgGAAAAAGBBBM0AAAAAACyIoBkAAAAAgAURNAMAAAAAsCCCZgAAAAAAFuT/Az68eb35lyrAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 481,
       "width": 717
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot data\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.plot(df['a_j'], label='$a_j$')\n",
    "\n",
    "# Cosmetic adjustments\n",
    "ax = plt.gca()\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.spines['top'].set_visible(False)\n",
    "plt.xlim(left=0)\n",
    "plt.ylim(bottom=0)\n",
    "plt.xlabel('$j$')\n",
    "plt.ylabel('Days')\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytical Solution "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=’Observation’></a>\n",
    "#### Observation\n",
    "\n",
    "Notice that $a_{1} = E(T_{1})$ where $T_{1} \\sim \\text{FS}(p)$ and $p$ is the probability the last remaining ball is drawn in the daily lotto drawing. Here, I'm following the convention that a random variable $X \\sim \\text{FS}(p)$ counts the number of Bernoulli trials (where each trial \"succeeds\" with probability $p$) until the the first \"success\", _inclusive_ of the \"success\".\n",
    "\n",
    "This means that $p = \\frac{1}{7} \\implies E(T_{1}) = 7$ by symmetry because all arrangements of 35 balls are equally likely and the missing one has $5(\\frac{1}{35}) = \\frac{1}{7}$ probability of being among the first 5 drawn. This is confirmed by the simulation above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zooming out for a second, we notice that the structure of the problem is similar to the [\"coupon collector\" problem](https://en.wikipedia.org/wiki/Coupon_collector%27s_problem) from Chapter 4. \n",
    "\n",
    "The \"coupon collector\" problem asked how long we had to wait (in terms of trials, where each trial produces a coupon) until all coupon types (or varieties) were collected. The wrinkle in the \"Mass Cash\" problem is that each trial (day) is a set of $5$ balls whereas in the normal \"coupon collector\" problem each trial consists of one draw."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Strategy: Simple, Non-Trivial Case\n",
    "\n",
    "Hmmmm...I'm still a little stuck connecting the dots between the \"coupon collecting\" problem and my intuitive understanding of iterated expectation, so let's try to solve a simpler (but non-trivial) problem.\n",
    "\n",
    "Rather than $35$ balls and $5$ drawn daily without replacement, let's try to solve the same problem with $4$ balls and $2$ drawn daily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To check analytical result with numerical approximation\n",
    "df_simple = generate_sim_data(4, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to back into each of the $a_{j}$ by considering conditional expectations. We know that $a_{0} = 0$. Using the same reasoning above (see **[Observation](#Observation)**) we get $a_{1} = 2$.\n",
    "\n",
    "Now things get interesting. We notice that after Day 1, we _always_ get to a point where we need to calculate $a_{2}$ since the $2$ drawn balls, tautologically, are \"never-before-seen\" or new numbers.\n",
    "\n",
    "There are $3$ cases to condition on:\n",
    "\n",
    "1. $0$ new numbers are drawn (i.e., the exact same numbers are drawn on Day 2 as were drawn before \\[on Day 1\\])\n",
    "2. $1$ new number is drawn (i.e., the other number was already observed before \\[on Day 1\\])\n",
    "3. $2$ new numbers are drawn (i.e., we have observed all balls over Day 1 and Day 2)\n",
    "\n",
    "After we observe Day 2's result, we can use iterated expectation to solve for $a_{2}$ in terms of $a_{1}$ and $a_{0}$:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "a_{2} &= \\frac{4}{4!}(a_{2} + 1) + \\frac{16}{4!}(a_{1} + 1) + \\frac{4}{4!}(a_{0} + 1)\\\\\n",
    "a_{2} &= 1 + \\frac{1}{6}a_{2} + \\frac{2}{3}a_{1}\\\\\n",
    "a_{2} &= \\frac{14}{5} = 2.8\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.79608"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check a_2 result\n",
    "df_simple.loc[2, 'a_j']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We proceed similarly for $a_{3}$, even if we know in reality we will never observe a situation when we have exactly $3$ unobserved balls (from the fact that Day 1's result always leaves us with $2$ remaining, unseen balls).\n",
    "\n",
    "Looking at the pattern that's emerging it's clear that we're conditioning on the cases where there are $j \\geq 0$ balls remaining _after_ the last draw.\n",
    "\n",
    "But notice that the case that we need to see $3$ more balls again after the latest draw is impossible since we're drawing $2$ balls per day and $j=3$ implies we've only seen one ball. Even if we drew this ball again, we are guaranteed to see at least $1$ new ball. Furthermore, there's no way we can be done (i.e., having seen all the balls after the latest draw) because we are only drawing $2$ balls at a time and we're missing $3$ balls -- this means that the probability of being done ($j = 0$) is $0$.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "a_{3} &= 0(a_{3} + 1) + \\frac{12}{4!}(a_{2} + 1) + \\frac{12}{4!}(a_{1} + 1) + 0(a_{0} + 1)\\\\\n",
    "a_{3} &= \\frac{1}{2}\\left(\\frac{19}{5}\\right) + \\frac{1}{2}(3)\\\\\n",
    "a_{3} &= \\frac{17}{5} = 3.4\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.40352"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check a_3 result\n",
    "df_simple.loc[3, 'a_j']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating $a_{4}$ in the same way yields:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "a_{4} &= 0(a_{4} + 1) + 0(a_{3} + 1) + 1(a_{2} + 1) + 0(a_{1} + 1) + 0(a_{0} + 1)\\\\\n",
    "a_{4} &= \\frac{14}{5} + 1 = 3.8\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "This is like being in Day 0 and conditioning on all the events that could occur after Day 1's draw. As noted above, we will always observe $2$ new balls and always need to calculate the average number of days to get the remaining $2$ balls (and add it to the one day that transpired, i.e., Day 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.8044"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check a_4 result\n",
    "df_simple.loc[4, 'a_j']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General Analytical Solution\n",
    "\n",
    "With this example under our belt, we can transcribe the pattern we see into a recursive expression. Basically, **we are conditioning on the outcome of a draw knowing what $j$ is going into that draw**. Over time the number of \"remains-to-be-seen\" balls decreases (with enough days passing by!), so we know the outcomes ($k$ describes how many balls are left _after_ we see the draw) are bounded: either the draw gives new numbers ($k < j$) or the draw doesn't yield any new numbers ($k = j$).\n",
    "\n",
    "For \"Mass Cash\", we can formalize these ideas:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "a_{j} &= \\sum_{k=0}^{j} \\binom{5}{k}\\binom{30}{j-k} \\frac{j!(35-j)!}{35!} (a_{j-k} + 1)\\\\\n",
    "&= 1 + \\sum_{k=0}^{j} \\frac{\\binom{5}{k}\\binom{30}{j-k}}{\\binom{35}{j}} (a_{j-k})\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Notice that there is a Hypergeometric PMF embedded in the sum. Letting $N$ count the number of never-before-seen balls from the draw, $N \\sim \\text{HGeom}(5, 30, j)$. We can also express $a_{j-k} + 1$ as $E(T_{j} | N)$, where $T_{j}$ counts the number of days until the remaining $j$ balls are drawn at least once. Hence, the equivalence (by [\"Adam's Law\"](https://youtu.be/gjBvCiRt8QA?t=636)):\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "a_{j} =& \\sum_{k=0}^{j} \\binom{5}{k}\\binom{30}{j-k} \\frac{j!(35-j)!}{35!} (a_{j-k} + 1)\\\\\n",
    "\\Updownarrow&\\\\\n",
    "E(T_{j}) =& E(E(T_{j} | N))\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check our solution, let's compare $a_{j}$ with our numerical estimates. \n",
    "\n",
    "We can write the general formula as follows for the lotto scenario with $w$ balls and $v$ draws from these balls without replacement:\n",
    "\n",
    "$$\n",
    "a_{j} = \\sum_{k=0}^{j} \\binom{v}{k}\\binom{w-v}{j-k} \\frac{j!(w-j)!}{w!} (a_{j-k} + 1)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "@lru_cache\n",
    "def analytical_solution_mass_cash(w: int = 35, v: int = 5, j: int = 35) -> float:\n",
    "    \"\"\"\n",
    "    Calculates the analytical solution for expected number of trials until\n",
    "    all w unique items in a set are drawn where each trial is a random draw without\n",
    "    replacement of v items.\n",
    "    \n",
    "    Parameters:\n",
    "        w: number of unique items\n",
    "        v: number of draws without replacements\n",
    "        j: remaining number of items not yet drawn\n",
    "        \n",
    "    Returns:\n",
    "        a_j: expected number of additional days needed to collect all unique items\n",
    "    \"\"\"\n",
    "\n",
    "    # Base case\n",
    "    if j == 0:\n",
    "        return 0\n",
    "    \n",
    "    # Recursive case\n",
    "    ex_a_j_terms = 0\n",
    "    for k in range(1, j + 1):\n",
    "        p = comb(v, k) * comb(w-v, j-k) * factorial(j) * factorial(w-j) / factorial(w)\n",
    "        ex_a_j_terms += p * (analytical_solution_mass_cash(w, v, j-k) + 1)\n",
    "    \n",
    "    # Note: k = 0 case is separated because it would case stack overflow if \n",
    "    # included in the for-loop above\n",
    "    k = 0\n",
    "    p_j = comb(v, k) * comb(w-v, j-k) * factorial(j) * factorial(w-j) / factorial(w)\n",
    "    a_j = (ex_a_j_terms + p_j)  / (1 - p_j)\n",
    "        \n",
    "    return a_j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "    Analytical solution: 27.718\n",
      "    Numerical approximation: 27.714\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Check and compare solutions\n",
    "w = 35\n",
    "v = 5\n",
    "j = 35\n",
    "print(f\"\"\"\n",
    "    Analytical solution: {analytical_solution_mass_cash(w, v, j):0.3f}\n",
    "    Numerical approximation: {df.loc[j, 'a_j']:0.3f}\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Feel free to change `df` to `df_simple` and use `w=4`, `v=2`, `j=4` to verify the results from the simplified problem above.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE**: The simulation could be re-written to be more efficient!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}