2-basic-probability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic Probability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Experiments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sampling With and Without Replacement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "Plots.PyPlotBackend()"
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Listing 2.4: Fishing with and without replacement\n",
    "using StatsBase, Plots ; pyplot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "proportionFished (generic function with 2 methods)"
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function proportionFished(gF, sF, n, N, withReplacement = false)\n",
    "    function fishing()\n",
    "        fishInPond = [ones(Int64, gF); zeros(Int64, sF)]\n",
    "        fishCaught = Int64[]\n",
    "\n",
    "        for fish in 1:n\n",
    "            fished = rand(fishInPond)\n",
    "            push!(fishCaught, fished)\n",
    "            if withReplacement == false\n",
    "                deleteat!(fishInPond, findfirst(x->x==fished, fishInPond))\n",
    "            end\n",
    "        end\n",
    "        sum(fishCaught)\n",
    "    end\n",
    "\n",
    "    simulations = [fishing() for _ in 1:N]\n",
    "    proportions = counts(simulations, 0:n)/N\n",
    "    \n",
    "    if withReplacement\n",
    "        plot!(0:n, proportions,\n",
    "            line=:stem, marker=:circle, c=:blue, ms=6, msw=0,\n",
    "            label=\"With replacement\",\n",
    "            xlabel=\"n\",\n",
    "            ylims=(0, 0.6), ylabel=\"Probability\")\n",
    "    else\n",
    "        plot!(0:n, proportions,\n",
    "            line=:stem, marker=:xcross, c=:red, ms=6, msw=0,\n",
    "            label=\"Without replacement\")\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAv8UlEQVR4nO3df1xUdb7H8ffY+CMVcBHpWjACElj+mlDwR/aDTfNx3VIvWLleE1JCdnW73bpe27JfNzO32sq0x6JhZKsrZajXfj1adaGWG6mYuuWaikFgYshPNTHFzv1jcpLAkR+HgWFez8djHjPnfM98z2ems+zbc77zPRbDMAwBAADANJ3augAAAICOhoAFAABgMo8OWKdOndJnn32mU6dOtXUpAAAATh4dsL788ksNGzZMX375ZbP7qK6uNrEioHVxvMKTcLzC05h5zHp0wDLDuXPn2roEoNE4XuFJOF7hacw8Zr0+YAEAAJiNgAUAAGAya1sXAABAQ4qKilRWVtbWZcCLVFdXy8/P75LbBQQEyGazudyGgAUAaHeKiop0zTXX8CtxtEvdu3fXvn37XIYsAhYAoN0pKyvTqVOntHr1al1zzTVtXQ7gtG/fPk2fPl1lZWUELACAZ7rmmmsUFRXV1mUATcYgdwAAAJMRsAAAAExGwAIAADAZAQsAAMBkBCwAAACTEbAAAABMRsACAHRcp05Jd90lHTzoWD540LFswgSmTz75pJKSkpzLOTk5slgsys7Odq5LSUnRo48+qry8PN11112SpKqqKi1evLhOXzfffLM2btzY4pou9Prrr+vLL780tc+GJCYm6qWXXmr1/bSW1vqeCFgAgI7p1Cnpttukt96Sxo+Xdu1yPL/1lmN9C0NWbGxsnTCVlZWlESNG1Fv3y1/+UsOHD9ebb74pqeGA1RqaEhxqa2tbuZr2i4AFAEBT3HOPlJXleF1QIA0f7niWHOvvuadF3Y8cOVJHjhzR4cOHJUnZ2dl67LHHnAGrpKRERUVFGjVqlLKzs2W32yU5zmqdOHFCdrtdw4cPd/aXk5OjG264Qf3791dKSopzfWlpqeLi4jR48GANGjRIy5cvd7aFhIRo9+7dzuXhw4crOztbaWlpysvL03/+53/Kbrfr/fffr1d/SEiI5s+fr5iYGCUkJOjs2bN66KGHFBMTI7vdrjvvvFOVlZWSHGepZs6cqdGjRysiIkIJCQmqqamp1+fWrVs1atQoXXfddRo4cKBWrlzpbKuurlZSUpIGDRqkoUOHaubMmZJ0yf0mJydr7NixCg0N1cyZM7V9+3bdfPPNCgsL0wMPPODs/+jRo7rzzjsVExOjwYMHa8GCBXU+62OPPaZRo0YpNDRUCxculKRGfU/NRcACAHRMCxdKoaE/Lf/ww0+vw8Kkp59uUfddunTR6NGjlZWVpe+//14FBQWaMGGCDh8+rNOnTysrK0ujRo1St27d6rwvNTVVPj4+2r17t/Ly8pzrDx06pKysLH3xxRf68MMPlZubK0n63e9+p8jISH3++ef629/+poULF+rTTz91WVtSUpKGDx+uF198Ubt379aECRMa3K68vFzbtm3TmjVr9Nxzz6lHjx7avn27du/eXS+kbNu2TR9++KH27duniooKvfjii/X6i4qKUk5Ojnbt2qW///3v+p//+R9nAL3//vvVpUsX/eMf/9CePXv0hz/8QZIuud/PP/9c7777rvbv36+PP/5YzzzzjDZv3qzPP/9ca9as0d69eyVJCQkJmjNnjrZv365du3YpLy9P69atc/ZTVVWl3Nxc7dixQ88995y++eabRn9PzcGtcgAAHdPVV0uZmY4zVxeGq06dHOvDw1u8i/OXCfv166eYmBhJjjNbubm5ys7OVmxsbKP7uuuuu2S1WmW1WmW323Xo0CGNGjVKW7Zs0c6dOyVJgYGBiouL05YtWzRy5MgW15+YmCiLxSJJ2rhxo6qrq5WZmSlJOnPmjEJCQpzb3nnnnfLx8ZEkzZo1Sy+//LIefvjhOv2Vl5dr1qxZOnDggKxWq8rLy/XFF18oKChI7777rrZt26ZOnRzndvr06dOo/U6aNMkZUgcPHqzx48erc+fO6ty5s6699lodPHhQISEh2rp1q7799lvn+06ePKn9+/c7l6dNmyZJCggIUFhYmAoKCnTVVVe1+Du8GAIWAKBjOnhQio+vG64kx3JcnPThh44Q1gKxsbFauXKlbDabbr75ZknSTTfdpKysLGVlZen1119vdF8Xnum67LLLLjou6nwgkiSr1apz5845l0+fPt2k+nv27Ol8bRiGli5dqltvvbVR772wjvNSUlI0YcIEZWZmymKxKCoq6pI1XWq/P/9eGvqeDMOQJH366af1zhherJ/WHnfGJUIAQMe0YMFPY64kx5mr8woKHO0tFB0drdLSUq1Zs6ZOwMrIyFBJSYnzrNaFfH19VVNTozNnzjRqH2PHjtWrr74qSTp27JjWr1+vcePGSZLCw8O1bds2SdL27dvrnLHx9fVVdXV1oz/L5MmT9eKLL+rUj4P/T5065bz8Jklvv/22Tp48qXPnzik9PV1jx46t10dlZaX69esni8Wijz/+WHv27HG2TZw4Uc8//7x++DHwHjt2rFH7bYyePXsqNja2zo8HLhwf50pTv6fGImABADqm9HTp/CW60FBp586fxmTFxjraW6hz584aM2aMTpw4oQEDBkiSIiIidOLECY0ZM0adO3eu9x5/f3/NmDFDQ4YMqTPI/WJefvll7du3T4MHD1ZsbKweeeQRjRgxQpK0cOFCvfLKKxo6dKhee+01DRw40Pm+5ORkLVq0qNGDt+fPn6/o6GiNGDFCQ4YM0ciRI+sMoI+Ojtb48eN1zTXXqFevXrr//vvr9bF48WI99NBDstvteu2115x1StKLL76o77//XoMHD5bdbndeXrzUfhtrzZo1ys/P16BBgzR48GDFxcWpvLz8ku9r6vfUWBbj/Hk1D/TZZ59p2LBh2rlzp6KioprVR0VFhfz9/U2uDGgdHK/wJC05Xs34+y7JMRXDPfc4BrxffbXjsuGCBY5w1b178/v1MomJibLb7Q2GKm/T2GOTMVgAgI6re3fpx/mnJDlC1oXLQCshYAEAAJeaMlgfDozBAgAAMBkBCwAAwGQELAAAAJMRsAAAAExGwAIAADAZAQsA0GGdOSMtWybZ7ZK/v+N52TLH+pZ68sknlZSU5FzOycmRxWJRdna2c11KSooeffRR5eXl6a677pLkuOnwhTOOS9LNN9+sjRs3tryoJggJCWnWhJ7txRNPPNHkWwO5EwELANAhnTkj3Xab9LvfSXv2SJWVjuff/c6xvqUh6/yNns/LysrSiBEj6q375S9/qeHDh+vNH+ffaihgNVdr30+vPXvyyScJWAAAuNuKFdLmzQ23bd4s/Xh7v2YbOXJknfvdZWdn67HHHnMGrJKSEhUVFWnUqFHKzs6W3W6X5DirdeLECdnt9jq3ysnJydENN9yg/v37KyUlpcF9FhYWqlevXpo/f76ioqK0bNkyHT16VHfeeadiYmI0ePBgLbjgHoshISGaN2+ehg0bpvDwcD333HMN9vvCCy8oOjpadrtd0dHRys3Ndbbt27dP48eP15AhQzRkyBClpqZK0iX3u2DBAo0ePVrBwcFKTU1Venq6Ro0apZCQEGVkZDi33bFjhzOEXnfddVq3bl2dz/r444876z9/K5vz388NN9wgu92u0tLSS/8HczfDg+3cudOQZOzcubPZfZSXl5tYEdC6OF7hSVpyvJrx933oUMOQLv6w25vdtdMtt9xivPHGG8bp06eN0NBQwzAMo3///kZNTY2xZs0aIzY21jAMw8jKyjKGDh1qGIZhFBQUGH5+fnX6uemmm4zJkycbZ8+eNU6dOmWEhIQYn3zySb39FRQUGJKMVatWOdfdeuutRnZ2tmEYhnH27Flj/PjxxltvvWUYhmH069fPuPvuu40ffvjBOHbsmBEcHGz83//9n7Nt165dhmEYRmlpqbO/3NxcIzIy0tnf1VdfbfzlL39xth87dqxR+73//vsNwzCMgwcPGt26dTOeeuopwzAMY/v27UZAQIBhGIZRWVlp2O1248iRI86+g4ODjcOHDzs/69tvv20YhmF88MEHRkREhLMOSUZlZeXF/tO0msYem8zkDgDokIqKWtbeGOcvE/br108xMTGSHGe2cnNzlZ2drdjzN5tuhLvuuktWq1VWq1V2u12HDh3SqFGj6m3XuXNnTZ8+XZL03XffaevWrfr222+d7SdPntT+/fudy7NmzZLFYlFAQIDi4uK0ZcsWjR49uk6fu3bt0tNPP63y8nJZrVbt379fNTU1+uqrr3T69Gn9+te/dm4bEBDQqP2eH3MWHh6ubt26acqUKZKk4cOHq6KiQlVVVfrkk0/01Vdf6V//9V/r1LN//36FhYWpW7duiouLkySNGjVKhw4davT32dYIWACADslmc4y7ctXeUrGxsVq5cqVsNptuvvlmSdJNN92krKwsZWVlNekWM926dXO+vuyyyy46vqp79+7q1MkxwscwDEnSp59+Wuf9rlgsljrLZ86cUVxcnLKyshQdHa3jx4/Lz89P33///UX7aMx+f/55zi9bLBZZLBbV1tbKMAwNHDhQn3zySb33FxYWqmvXrs56L7vsMp07d65Rn7E9YAwWAKBDuuAHfs1qb4zo6GiVlpZqzZo1dQJWRkaGSkpKnGe1LuTr66uamhqdMeGnjD179lRsbGydQfMXjguTfrqPYEVFhTZs2KBbbrmlTh+nT5/WmTNnZPsxcS5dutTZFhkZqe7du2vt2rXOdWVlZY3ab2OMHj1aBQUF2rJli3Pd7t27G/Xd+Pj4qLq6ukn7cycCFgCgQ0pOlsaNa7ht3Djp3ntbvo/OnTtrzJgxOnHihAYMGCBJioiI0IkTJzRmzBh17ty53nv8/f01Y8YMDRkypM4g9+Zas2aN8vPzNWjQIA0ePFhxcXEqLy93tvfp00fDhg1TTEyM5s6dW+/yoK+vrxYuXKiYmBgNGzZMXbp0cbZZrVb97//+r9LT0zV48GANHTpUmZmZjdpvY/ziF7/Qe++9p0WLFmno0KG69tpr9dBDD+mHH3645HsffPBBjRs3rt0OcrcY58/zeaDPPvtMw4YN086dOxUVFdWsPioqKuTv729yZUDr4HiFJ2nJ8WrG33fJMRXDq69KaWmOMVc2m+PM1b33ShfkiA4rJCREGzdudP6CES3X2GOTMVgAgA6rSxdpzhzHA3AnAhYAAB1UYWFhW5fgtRiDBQAAYDICFgAAgMkIWAAAACYjYAEAAJiMgAUAAGAyAhYAAIDJCFiAJ8jIkMrKHK/LyhzLAIB2i3mwgPZuxQopJUWKiZFl1Sppxgxpxw7p+HHHvUCADmzfvn1tXQJQR2OPSQIW0J5lZDjClWFI27ap1/Dh0smTjraUFMnXV5o6tW1rBFpBQECAunfvrunTp7d1KUA93bt3V0BAgMttCFhAezZ2rBQTI23bJkmynA9XkmP9xe5kC3g4m82mffv2qez8pXHADaqrq+Xn53fJ7QICAmSz2VxuQ8AC2rOAAGnTJik8XDpx4qf1vr7SO+9IvXu3XW1AK7PZbJf8PzHATC25QfnPMcgdaM/KyqTbb68briTH+Kvbb5fKy9umLgCASwQsoD3bssUxoL0h27dLmze7tx4AQKMQsID2bOpUKTVVslikQYMc6wYOdCynpjLAHQDaKcZgAe1dcrJjzFVgoHTLLdKSJdKxY4QrAGjH3HoG6+DBgxo9erQiIiIUHR2tvXv3NrhdSEiIIiMjZbfbZbfb9eabb7qzTKD9mTpV6tXL8foXvyBcAUA759YzWLNnz1ZycrISExP19ttvKzExUTsuMr7kzTfflN1ud2d5AAAApnDbGazS0lLl5eU5J42Lj49XcXGx8vPz3VUCAACAW7jtDFZxcbH69u0rq9WxS4vFIpvNpqKiIoWHh9fbfsaMGTIMQzExMVq8eLH69Olz0b7nzp0rPz8/xcXFKT4+vkl1VVZWNu2DAG3ksupq+ckxEd65ioq2Lge4JP6+wtO05Jj9+fxZ7XKQ+8cffyybzaazZ89qwYIFSkhI0Pvvv3/R7ZctW6aoqKhm78+sScWAVvXj7MJ+fn4Sxyw8BH9f4WnMOmbdFrCCg4NVUlKi2tpaWa1WGYahoqKiBmfpPb+uc+fOuv/++xUREeGuMgEAAFrMbWOwAgMDFRUVpdWrV0uSMjMzFRQUVO/y4Hfffaeqqirn8tq1a3Xddde5q0wAAIAWc+slwuXLlysxMVGLFi2Sr6+v0tPTnW1JSUmaOHGiBg0apPj4eJ07d06GYSgsLExvvPGGO8sEAABoEbcGrMjISOXm5jbYlpaW5ny9a9cud5UEAABgOm6VAwAAYDICFgAAgMkIWAAAACYjYAEAAJiMgAUAAGAyAhYAAIDJCFgAAAAmI2ABAACYjIAFAABgMgIWAACAyQhYAAAAJiNgAQAAmIyABQAAYDICFgAAgMkIWAAAACYjYAEAAJiMgAUAAGAyAhYAAIDJCFgAAAAmI2ABAACYjIAFAABgMgIWAACAyQhYAAAAJiNgAQAAmIyABQAAYDICFgAAgMkIWAAAACYjYAEAAJiMgAUAAGAyAhYAAIDJCFgAAAAmI2ABAACYjIAFAABgMgIWAACAyQhYAAAAJiNgAQAAmIyABQAAYDICFuABCgqkR165Ur/WX/TIK1eqoKCtKwIAuELAAtq5Vaukq6+WFr32L8rQr7XotX9RRIRjPQCgfSJgAe1YQYE0a5Z07lzd9bW1UlKSOJMFAO0UAQtox9LS6oer82prHe0AgPaHgAW0Y1995bqdM1gA0D4RsIB2LCzMdXtoqHvqAAA0DQELaMeSkiSrteE2q9XRDgBofwhYQDsWGuoYZ/XzkGW1SitXcgYLANorAhbQziUkSAcOSDMnlklyPB84IM2Y0caFAQAuioAFeIDQUGnOncckOZ45cwUA7RsBCwAAwGQELAAAAJMRsAAAAEzmtoB18OBBjR49WhEREYqOjtbevXtdbp+eni6LxaKNGze6p0AAAACTuC1gzZ49W8nJyTpw4IDmz5+vxMTEi25bWFioV199VSNHjnRXeQAAAKZxS8AqLS1VXl6epk+fLkmKj49XcXGx8vPz6237ww8/KCkpSUuXLlXXrl3dUR4AAICpLjJHtLmKi4vVt29fWX+cLdFischms6moqEjh4eF1tn3hhRd0/fXXa9iwYY3uf+7cufLz81NcXJzi4+ObVFtlZWWTtgfayomTJ53PFRUVbVwNcGn8fYWnackx6+/vX2fZLQGrsb744gtlZmbq448/btL7li1bpqioqGbv9+dfCtAe+fT89sfnnhyz8Bgcq/A0Zh2zbglYwcHBKikpUW1traxWqwzDUFFRkWw2W53t/v73v6uwsFBXX321JOno0aNKTk5WSUmJfvOb37ijVAAAgBZzyxiswMBARUVFafXq1ZKkzMxMBQUF1bs8+Jvf/EYlJSUqLCxUYWGhRo4cqRUrVhCuAACAR3HbrwiXL1+u5cuXKyIiQosXL1Z6erqzLSkpSZs2bXJXKQAAAK3KbWOwIiMjlZub22BbWlpag+uzs7NbsSIAAIDWwUzuAAAAJiNgAQAAmIyABQAAYDICFgAAgMkIWAAAACYjYAEAAJiMgAUAAGAyAhYAAIDJCFgAAAAmI2ABAACYjIAFAABgMgIWAACAyQhYAAAAJiNgAQAAmIyABQAAYDICFgAAgMlcBqzU1FQdP37cXbUAAAB0CC4D1gMPPKC+fftqxowZ+uijj9xVEwAAgEdzGbCOHDmiZ599Vv/85z8VGxur8PBwLVq0SN9884276gMAAPA4LgNWr169NGfOHOXl5Wn37t267bbb9NJLLykkJES/+tWvlJmZqbNnz7qrVgAAAI/Q6EHuQ4YM0UsvvaTdu3fr+uuv1wcffKA77rhDV111lR5//HHV1NS0Zp0AAAAeo1EByzAMffDBB5oyZYrCwsL05Zdfat68efrkk0+UkpKipUuXavr06a1dKwAAgEewumo8dOiQXnvtNb3xxhs6cuSIxo0bpzVr1mjSpEmyWh1vHTlypIYPH66pU6e6pWAAAID2zmXAuvrqq3XVVVfpnnvu0axZs9SvX78GtxswYIBGjBjRKgUCAAB4GpcBa9OmTZowYYI6dXJ9JTEiIkJZWVmmFgYA8FAZGdLYsVKnTlJZmbRli8RVDngZl8lp/fr1+vrrrxts+/rrrzVz5sxWKQoA4KFWrJCmTZNuu02WY8ekX/3KsbxiRVtXBriVy4C1atUqHTt2rMG2srIyrVq1qlWKAgB4oIwMKSVFMgxp2zb1Gj5c2r7dsZyS4mgHvITLgGUYhiwWS4NtBw8eVO/evVulKACABxo7VoqJcS5aTp78qS0mRho3rg2KAtpGvTFYf/rTn/SnP/1JkmSxWDRt2jRdfvnldbY5ffq0CgsLdccdd7inSgBA+xcQIG3aJIWHSydO/LTe11d65x2Jf5TDi9QLWFdeeaWGDRsmSfriiy8UGRmpPn361NmmS5cuuuaaazRr1iz3VAkAaP/KyqTbb68briTp+HHH+vfeI2TBa9QLWJMmTdKkSZOcy48++qjCwsLcWhQAwANt2SLt2NFw2/bt0ubN/JoQXsPlGKz09HTCFQCgcaZOlVJTJYtFGjTIsW7gQMdyairhCl6l3hms++67T//1X/8lm82m++67z+WbLRaLlixZ0mrFAQA8THKyY8xVYKB0yy3SkiXSsWOEK3idegHrnXfe0axZs2Sz2bRp06aL/opQImABABowdar02WeO17/4hSNoAV6mXsAqKChwvi4sLHRnLQAAAB2C63vgAAAAoMnqncFav359kzqIi4szrRgAAICOoF7AmjJlSqPfbLFYdO7cOVMLAgAA8HQux2ABAACg6eoFrH79+rVFHQAAAB1GvYBVUVGhXr16qVOnTqqoqLhkB/7+/q1SGAAAgKeqF7D69Omj3NxcxcTEKCAgwOU8WJIYgwUAqKOgQEp75Up9pb8o7JUrlbRACg1t66oA96oXsF577TX179/f+fpSAQsAgPNWrZJmzZLOnfsXSb+WXpOefUNKS5MSEtq6OsB96gWshAv+F5CYmOjOWgAAHqyg4Hy4qru+tlZKSpJuvJEzWfAejZpotLq6Wjk5OVq3bp1ycnJUXV3d2nUBADxMWlr9cHVeba2jHfAW9c5gXeiHH37QggULtHTpUn333XfO9T169NDcuXO1cOFCXXbZZa1eJACg/fvqK9ftzAIEb+IyYM2bN09Lly7V73//e02ZMkVXXHGFvv32W61bt06LFy/WmTNn9Mc//tFdtQIA2rGwMNftXB6EN3EZsF5//XU99dRTmj9/vnNdYGCgBg8erMsvv1zPP/88AQsAIMkxzurZZx2XA3/OanW0A97C5Risc+fOKSoqqsG2YcOGMUUDAMApNNQxzsr6s3+6W63SypWcwYJ3cRmwpkyZooyMjAbbMjIyuNEzAKCOhATpwAFp5sQySY7nAwekGTPauDDAzepdIly/fr3z9U033aSHH35YsbGxmjx5sgIDA1VaWqoNGzbo0KFDevrpp5u0s4MHDyohIUFlZWXy8/PT66+/roEDB9bb7tZbb9XRo0fVqVMn+fj46OWXX9Z1113XjI8HAHC30FBpzp3H9NqmAM2585hCQwPauiTA7eoFrClTptTb6JtvvtFHH31Ub/0999yjGU34Z8ns2bOVnJysxMREvf3220pMTNSOHTvqbffWW2+pV69ekqQNGzYoMTFRe/bsafR+AAAA2lK9gFXQSr+jLS0tVV5env76179KkuLj4zV37lzl5+crPDy8zrbnw5XkmIOL2eQBAIAnqRew+vXr1yo7Ki4uVt++fWX9cfSjxWKRzWZTUVFRvYAlSTNmzFBWVpYk6f3333fZ99y5c+Xn56e4uDjFx8c3qa7KysombQ+0lRMnTzqfG3MjdqAtcbzCE7UkE/j7+9dZdjlNw4VOnTql06dPX7JDs7zxxhuSpFWrVmn+/PkuQ9ayZcsu+mvHxmitzwCYyafntz8+9+SYRbvH8QpPZdbx6vJXhIZh6KmnnlJQUJB8fHzUp0+feo/GCg4OVklJiWp/nCDFMAwVFRXJZrO5fF9CQoKysrJUXl7e6H0BAAC0JZcB68UXX9QLL7ygOXPmyDAMPfLII3rssccUERGhkJAQvfrqq43eUWBgoKKiorR69WpJUmZmpoKCgupdHqyqqtKRI0ecyxs3blTv3r35FxAAAPAYLi8Rrly5Uk8++aTmzJmjRx55RJMnT1ZUVJQeffRRTZw4Ufn5+U3a2fLly5WYmKhFixbJ19dX6enpzrakpCRNnDhRQ4cO1R133KGamhp16tRJffr00bvvvstAdwAA4DFcBqzCwkLZ7XZddtll6ty5s6qqqiRJnTp10m9/+1slJSVp0aJFjd5ZZGSkcnNzG2xLu+A269u3b290nwAAAO2Ny0uEvXv31skffwlis9n02WefOdvKysp06tSp1q0OAADAA7k8g3X99ddrx44dmjBhgqZNm6YnnnhCR48eVefOnfXqq6/qlltucVedAAAAHsNlwHriiSf0zTffSJIefvhhVVVVae3ataqpqdG4ceO0dOlStxQJAADgSVwGrMjISEVGRkqSunbtqiVLlmjJkiVuKQwAAMBTNXqi0cOHD6ukpERXXnmlrrrqqtasCQAAwKO5HOQuSStWrJDNZlO/fv00cuRI2Ww2BQcHa/ny5e6oDwAAwOO4PIP1zDPP6JFHHtHdd9+tKVOm6IorrtC3336rdevW6be//a0qKir0+9//3l21AgAAeASXAWvp0qWaN2+e/vCHP9RZf/vtt+uKK67Q0qVLCVgAAAA/4/IS4fHjxzV27NgG22699VadOHGiVYoCAADwZC4D1vjx47Vly5YG2zZv3sw8WAAAAA2od4nwwtnak5KSNHv2bJWWlmry5MkKDAxUaWmpNmzYoL/97W8MdAcAAGhAvYA1fPjwOjdWNgxDq1at0qpVq2SxWGQYhrPttttu07lz59xTKQAAgIeoF7CysrLaog4AAIAOo17Auummm9qiDgAAgA6jUTO57927Vzk5OaqoqJC/v7/GjBmjgQMHtnZtAAAAHsllwPr+++919913KzMzU4ZhqGvXrvr+++9lsVg0ZcoU/fnPf1aXLl3cVSsAAIBHcDlNw8MPP6z33ntPqampqqqqUk1NjaqqqpSamqr33ntPDz/8sLvqBAAA8BguA1ZGRoaeeeYZ3XvvvfL19ZUk+fr66t5779XTTz+ttWvXuqVIAAAAT+IyYFVUVGjAgAENtg0YMEAVFRWtUhQAAIAncxmwBgwYoD//+c8Ntq1evfqi4QsAAMCbuRzk/uijj+qOO+5QYWGh4uPjdcUVV6i0tFRvv/22cnNztW7dOnfVCQAA4DFcBqy4uDht2LBBTz75pB588EEZhiGLxSK73a4NGzbo9ttvd1edAAAAHuOiAevMmTN69913ZbfbtXPnTn333XeqqqpSr1691KNHD3fWCAAA4FEuOgarS5cumjZtmoqKiiRJPXr00FVXXUW4AgAAuIRLDnI/H7AAAADQOC4D1jPPPKOFCxcqLy/PXfUAAAB4PJeD3P/7v/9b5eXlGjFihHr37q0rrrhCFovF2W6xWLRnz55WLxIAAMCTuAxYw4cPd1cdAAAAHUaDAeuf//ynUlNTVVZWpiuvvFJTpkzRuHHj3F0bAACAR6oXsHJycjR27FidPXtWffr0UXl5udLS0vTKK68oJSWlLWoEAADwKPUGuT/++OMaMGCACgsLdfToUZWXl2vy5MlasGBBW9QHAADgceoFrM8//1yPPfaYgoODJUm+vr764x//qIqKChUXF7u9QAAAAE9TL2CVlZUpKCiozrrzYausrMw9VQEAAHiwBufBunAqBgAAADRNg78ijI2NVadO9bPXDTfcUGe9xWJRdXV161UHAADggeoFrMcff7wt6gAAAOgwCFgAAAAmc3kvQgAAADQdAQsAAMBkBCwAAACTEbAAAABMRsACAAAwGQELAADAZAQsAAAAkxGwAAAATEbAAgAAMBkBCwAAwGQELAAAAJMRsAAAAExGwAIAADCZ2wPWwYMHNXr0aEVERCg6Olp79+6t03769GlNnjxZERERGjp0qMaNG6f8/Hx3lwkAANBsbg9Ys2fPVnJysg4cOKD58+crMTGx3jbJycnav3+/9uzZo0mTJikpKcndZQIAADSbWwNWaWmp8vLyNH36dElSfHy8iouL65yh6tatmyZMmCCLxSJJGjlypAoLC91ZJgAAQItY3bmz4uJi9e3bV1arY7cWi0U2m01FRUUKDw9v8D1LlizRpEmTXPY7d+5c+fn5KS4uTvHx8U2qqbKysknbA23lxMmTzueKioo2rgZwjeMVnqglmcDf37/OslsDVlMtWrRI+fn52rp1q8vtli1bpqioqGbv5+dfCtAe+fT89sfnnhyzaPc4XuGpzDpe3RqwgoODVVJSotraWlmtVhmGoaKiItlstnrbPv/881q/fr22bNmi7t27u7NMAACAFnHrGKzAwEBFRUVp9erVkqTMzEwFBQXVuzz4wgsvaO3atdq8ebN69erlzhIBAABazO2XCJcvX67ExEQtWrRIvr6+Sk9Pd7YlJSXpxhtv1IMPPqiwsDDFxsZKkrp27apt27aZWkdBgZSWJn35ZQ8NGCAlJUmhoabuAgAAeCm3B6zIyEjl5uY22JaWliZJmjFjRqvWsGqVNGuWdO6cJHWVJD37rCNwJSS06q4BAIAX8LqZ3AsKLgxXP6mtdZzFKihom7oAAEDH4XUBKy2tfrg6r7bW0Q4AANASXhewvvrKdTtnsAAAQEt5XcAKC3PdzkB3AADQUl4XsJKSJOtFhvZbrY52AACAlvC6gBUa6hhn9fOQZbVKK1dyBgsAALSc1wUsyTEVw4ED0syZjuWZMx3LrTw7BAAA8BJeGbAkx5mqOXMcr+fM4cwVAAAwj9cGLAAAgNZCwAIAADAZAQsAAMBkBCwAAACTEbAAAABMRsACAAAwGQELAADAZAQsAAAAkxGwAAAATEbAAgAAMBkBCwAAwGQELAAAAJMRsAAAAExGwAIAADAZAQsAAMBkBCwAAACTEbAAAABMRsACAAAwGQELAADAZAQsAAAAkxGwAAAATEbAAgAAMBkBCwAAwGQELAAAAJMRsAAAAExGwAIAADAZAQsAAMBkBCwAAACTEbAAAABMRsACAAAwGQELAADAZAQsAAAAkxGwAAAATOa9ASsjQ6qsdLyurHQsAwAAmMA7A9aKFdK0adL99zuW/+M/HMsrVrRpWQAAoGPwvoCVkSGlpEiGIX3xuWPd3i8cyykpnMkCAAAt5n0Ba+xYKSam4baYGGncOPfWAwAAOhzvC1gBAdKmTZKPT931vr7SO+9IvXu3TV0AAKDD8L6AVVYm3X67dOJE3fXHjzvWl5e3TV0AAKDD8L6AtWWLtGNHw23bt0ubN7u3HgAA0OF4X8CaOlVKTZUsFql/uGNdWH/Hcmqqox0AAKAFrG1dQJtITnaMuaoKk34jad48qVc84QoAAJjC+85gnTd16k8D3X18CFcAAMA0bg1YBw8e1OjRoxUREaHo6Gjt3bu3we3uu+8+hYSEyGKxaPfu3e4sEQAAoMXcGrBmz56t5ORkHThwQPPnz1diYmKD202ZMkU5OTnq16+fO8sDAAAwhdsCVmlpqfLy8jR9+nRJUnx8vIqLi5Wfn19v2xtvvFFBQUHuKg0AAMBUbhvkXlxcrL59+8pqdezSYrHIZrOpqKhI4eHhLep77ty58vPzU1xcnOLj4xv9vhMnTzqfKyoqWlQD0No4XuFJOF7hiSorK5v9Xn9//zrLHeJXhMuWLVNUVFST3+fT89sfn3vW+2KA9objFZ6E4xWeyqzj1W2XCIODg1VSUqLa2lpJkmEYKioqks1mc1cJAAAAbuG2gBUYGKioqCitXr1akpSZmamgoKAWXx4EAABob9z6K8Lly5dr+fLlioiI0OLFi5Wenu5sS0pK0qZNmyQ5fm0YFBSkw4cPa/z48YQwAADgUdw6BisyMlK5ubkNtqWlpTlfL1++3F0lAQAAmM57Z3IHAABoJQQsAAAAk3WIaRoAAACaq6BASkuTvvyyhwYMkJKSpNDQlvVJwAIAAF5r1Spp1izp3DlJ6ipJevZZR+BKSGh+v1wiBAAAXqmg4MJw9ZPaWsdZrIKC5vdNwAIAAF4pLa1+uDqvttbR3lwELAAA4JW++sp1O2ewAAAAmigszHV7Swa6E7AAAIBXSkqSrBf5uZ/V6mhvLgIWAADwSqGhjnFWPw9ZVqu0ciVnsAAAAJolIUE6cECaOdOxPHOmY3nGjJb1S8ACAABeLTRUmjPH8XrOnJZPMioRsAAAAExHwAIAADAZAQsAAMBkBCwAAACTEbAAAABMRsACAAAwGQELAADAZAQsAAAAkxGwAAAATEbAAgAAMBkBCwAAwGQELAAAAJMRsAAAAExGwAIAADAZAQsAAMBkBCwAAACTEbAAAABMRsACAAAwGQELAADAZAQsAAAAkxGwAACAd8vIkCorHa8rKx3LLUTAAgAA3mvFCmnaNOn++x3L//EfjuUVK1rULQELAAB4p4wMKSVFMgzpi88d6/Z+4VhOSWnRmSwCFgAA8E5jx0oxMQ23xcRI48Y1u2sCFgAA8E4BAdKmTZKPT931vr7SO+9IvXs3u2sCFgAA8E5lZdLtt0snTtRdf/y4Y315ebO7JmABAADvtGWLtGNHw23bt0ubNze7awIWAADwTlOnSqmpksUi9Q93rAvr71hOTXW0N5PVpBIBAAA8T3KyY8xVVZj0G0nz5km94lsUriTOYAEAAG83depPA919fFocriQCFgAAgOkIWAAAACYjYAEAAJiMgAUAAGAyAhYAAIDJCFgAAAAmI2ABAACYjIAFAABgMrcGrIMHD2r06NGKiIhQdHS09u7d26LtAAAA2iO3BqzZs2crOTlZBw4c0Pz585WYmNii7QAAANojt92LsLS0VHl5efrrX/8qSYqPj9fcuXOVn5+v8PDwJm8nSTU1NZKkffv2NaumfQUFkmocz5/VNKsPwF04XuFJOF7hacw4ZgcMGKDu3bs7Fgw3ycvLMyIiIuqsi46ONrZu3dqs7QzDMFavXm1I4sGDBw8ePHjwaPPHzp07nRnFbWewWsP48eO1evVqhYSE6PLLL2/rcgAAgBcbMGCA87XbAlZwcLBKSkpUW1srq9UqwzBUVFQkm83WrO0kKSAgQP/+7//uro8AAADQKG4b5B4YGKioqCitXr1akpSZmamgoKB646oaux0AAEB7ZTEMw3DXzvbv36/ExESVl5fL19dX6enpGjx4sCQpKSlJEydO1MSJE11uBwAA0N65NWABAAB4A2ZyBwAAMJnXBixmi4cnue+++xQSEiKLxaLdu3e3dTnARZ0+fVqTJ09WRESEhg4dqnHjxik/P7+tywJcuvXWWzVkyBDZ7XbdcMMN2rVrV4v79NqAxWzx8CRTpkxRTk6O+vXr19alAJeUnJys/fv3a8+ePZo0aZKSkpLauiTApbfeekv/+Mc/tHv3bj3wwAOmZAKvDFjnZ4ufPn26JMds8cXFxfwrC+3WjTfeqKCgoLYuA7ikbt26acKECbJYLJKkkSNHqrCwsG2LAi6hV69eztfV1dXO47clPHqi0eYqLi5W3759ZbU6Pr7FYpHNZlNRURHTQQCAiZYsWaJJkya1dRnAJc2YMUNZWVmSpPfff7/F/XllwAIAtL5FixYpPz9fW7dubetSgEt64403JEmrVq3S/PnzWxyyvPIS4YWzxUtyOVs8AKDpnn/+ea1fv14ffPDBTze/BTxAQkKCsrKyVF5e3qJ+vDJgMVs8ALSeF154QWvXrtXmzZvrjG0B2qOqqiodOXLEubxx40b17t1b/v7+LerXaycaZbZ4eJLZs2frvffe09GjR9W7d2/5+Pjwowy0S4cPH1ZwcLDCwsLk4+MjSeratau2bdvWxpUBDfv66691xx13qKamRp06dVKfPn30/PPPy263t6hfrw1YAAAArcUrLxECAAC0JgIWAACAyQhYAAAAJiNgAQAAmIyABQAAYDICFgAAgMkIWAAAACYjYAEAAJiMgAUAAGAyAhYAr5SYmKhBgwYpOztb1113nXr06KGYmBjt3LmzrUsD0AEQsAB4raNHj+q+++7TvHnz9NZbb+n06dP6t3/7N509e7atSwPg4axtXQAAtJWKigp99NFHGjhwoCSpR48eio2N1bZt2zRmzJg2rg6AJ+MMFgCvdeWVVzrDlSRde+21kqTDhw+3VUkAOggCFgCv1atXrzrLXbp0kSSdPn26DaoB0JEQsAAAAExGwAIAADAZAQsAAMBkBCwAAACTWQzDMNq6CAAAgI6EM1gAAAAmI2ABAACYjIAFAABgMgIWAACAyQhYAAAAJiNgAQAAmOz/AYbvaSmL++X9AAAAAElFTkSuQmCC"
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10^6\n",
    "goldFish, silverFish, n = 3, 4, 3\n",
    "proportionFished(goldFish, silverFish, n, 3)\n",
    "plot()\n",
    "proportionFished(goldFish, silverFish, n, N)\n",
    "proportionFished(goldFish, silverFish, n, N, true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lattice paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "MersenneTwister(UInt32[0x0000000c], Random.DSFMT.DSFMT_state(Int32[589013436, 1072892337, 1395473285, 1073017187, -346481041, 1072887813, 23263411, 1072776800, -1379277159, 1073038715  …  -173461052, 1072708370, 859747322, 1072998707, -624312252, -2047062377, 1852138555, -1338826024, 382, 0]), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000  …  0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000], 1002, 0)"
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Listing 2.5: Lattice paths\n",
    "using Random, Combinatorics, Plots, LaTeXStrings; pyplot()\n",
    "Random.seed!(12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "(5, 100000)"
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n, N = 5, 10^5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "isUpperLattice (generic function with 1 method)"
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function isUpperLattice(v)\n",
    "    for i in 1:Int(length(v)/2)\n",
    "        sum(v[1:2*i-1]) >= i ? continue : return false\n",
    "    end\n",
    "    return true\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "0.16666666666666666"
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "omega = unique(permutations([zeros(Int, n);ones(Int, n)]))\n",
    "A = omega[isUpperLattice.(omega)]\n",
    "pA_modelI = length(A)/length(omega)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "randomWalkPath (generic function with 1 method)"
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function randomWalkPath(n)\n",
    "    x, y = 0, 0\n",
    "    path = []\n",
    "    while x<n && y<n\n",
    "        if rand() < 0.5\n",
    "            x += 1\n",
    "            push!(path,0)\n",
    "        else\n",
    "            y += 1\n",
    "            push!(path,1)\n",
    "        end\n",
    "    end\n",
    "    append!(path, x<n ? zeros(Int64,n-x) : ones(Int64,n-y))\n",
    "    return path\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "4899.4"
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pA_modelIIest = sum([isUpperLattice(randomWalkPath(n)) for _ in 1:N]) / n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model I: 0.16666666666666666\t Model II: 4899.4\n"
     ]
    }
   ],
   "source": [
    "println(\"Model I: \",pA_modelI, \"\\t Model II: \", pA_modelIIest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "plotPath (generic function with 1 method)"
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function plotPath(v, l, c)\n",
    "    x, y = 0, 0\n",
    "    graphX, graphY = [x], [y]\n",
    "    for i in v\n",
    "        if i == 0\n",
    "            x += 1\n",
    "        else\n",
    "            y += 1\n",
    "        end\n",
    "        push!(graphX, x), push!(graphY,y)\n",
    "    end\n",
    "    plot!(graphX, graphY,\n",
    "        la=0.8, lw=2, label=l, c=c, ratio=:equal, legend=:topleft,\n",
    "        xlims=(0,n), ylims=(0,n),\n",
    "        xlabel=L\"East\\rightarrow\", ylabel=L\"North\\rightarrow\")\n",
    "end     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4iklEQVR4nO3deVyVdf738feRg9uQkuJWKqbmEgoI7qZJpThq7s2YWVqpSXeZY1OWW0451pSZ2qYtY2M6LjMa9jNGR/1p6qROLoiYaaRkjQspgoN6QQfO/Yc35xYkBDzwPcvr+XjwCM51POfNFcub7/U512VzOp1OAQAAwG0qmQ4AAADgayhYAAAAbkbBAgAAcDMKFgAAgJtRsAAAANyMggVJUpMmTdSyZUtFRkYqMjJSK1euNB0JAACvZTcdAJ5j5cqVioyMNB0DAACvxwoWSuXSpUvat2+fLl26ZDoKAAAey8aJRiFdOURYo0YNOZ1OdezYUa+++qrq1Klzzf327dun6OhodenSRTVr1tSQIUM0dOhQA4k9z4ULF1SjRg3TMTwO++Va7JNrHW4zRsEZ35uO4TGcTqdsNpskKezSXsNpUBYcIoQkadu2bWrcuLF+/vlnTZs2TaNGjVJCQsIv3v/tt99WVFRUBSb0DrVq1TIdwSOxX67FPimossO68k6lAF2qHmI2jGHZzlxNuZCogVUbqkeVeqbjoIwoWJAkNW7cWJIUGBioiRMnqkWLFoYTAfBHl6qHqEPaL/9x5+ssy9KgQYO0Z9t/NfPT19Xh7rtNR0IZMYMFXbx4URkZGa6Ply9frnbt2pkLBAB+KL9cbdu2TevWrdPdlCuvxgoWdObMGQ0dOlS5ublyOp1q2rSplixZUurHOXHihM6ePVsOCb1DZmamatasaTqGx/GG/RISEuJaxQVMeeqppyhXPoSCBTVt2lT79++/occ4ceKEWrduzasL4ZWqV6+uw4cPU7Jg1LRp0zRy5EjdddddpqPADShYcIuzZ8/q0qVLWrp0qVq3bm06DlBihw8f1siRI3X27FkKFiqcZVmaMmWKpk2bptDQUIWGhpqOBDehYMGtWrduzasLAaAErp65GjZsmLp27Wo6EtyIIXcAACpY4YF2ypXvoWABAFCB8vLyNGTIEAbafRwFCwCAClSpUiUNHjyYcuXjKFgAAFQAy7K0atUqSdLYsWMpVz6OggUAQDnLn7kaPXq0Tpw4YToOKgAFCwCAclR4oJ3TgfgHChb8is1mK3BZIElq0qSJEhMTjeQpzfPPmzdPp0+fdn28detWrV+/3vXxyZMn1b179/KIWCrx8fHatWuX6+OtW7cqMjLSXCDAIC5/478oWEAFcjgcZf631ytYt9xyi7Zv335D+dyhcMEC/JnD4ZDT6aRc+SEKFnCVJk2a6Nlnn1V0dLSaN2+u119/vUTbvv32W/Xr108dOnRQeHi43n77bdc2m82mF198UR06dNALL7xQ7PPPnTtXHTp0UGRkpDp06KCdO3dKkl566SWdPHlSv/3tbxUZGanExEQtXLhQy5YtU2RkpF566SWlpqYqODjY9Vg7d+7UnXfeqYiICIWHh2vt2rXXzXq1mTNnaujQobr77rvVqlUr3XfffTp37pwkafPmzerSpYvatWunsLAwffTRR5KkhIQEffbZZ3r99dcVGRmpDz/8UNKVXzJPPPGEIiIiFBYWpj179pT0fwnglSzL0vHjxxUUFKT169dTrvwQZ3JHhTp79sqbu4WEXHlzhzNnzmjPnj06d+6coqKi1K1bN9dJAIva1qlTJz3wwANaunSpWrVqpUuXLqlz587q1KmTOnToIEkKCAjQV199dd3nfuihhzRp0iRJ0q5duzR69Gh98803mjFjhv785z9r5cqVrsNt48ePV0ZGhubNmydJSk1NdT1Oenq6Bg0apL///e/q3r278vLylJGRodzc3Otmvdr27duVlJSk+vXr64knntALL7yg999/X1FRUdqxY4cCAgKUnp6udu3aKTY2Vn379tWAAQMUGRmpiRMnSrqy0vbNN9/oo48+0rvvvquFCxdq6tSp2rBhQxn/DwGeLf+wYEpKig4fPqzAwEDTkWAABQsVas0a6f333f+448ZdeSsrm83mev+xxx6TzWZTSEiIhgwZok2bNrkKVlHbgoODdejQIQ0fPtz1GP/973/19ddfu0rLo48+WqIc+/fv1x//+EedO3dOdrtdR44c0eXLl1WtWrVSfT47d+5Uy5YtXTNZlSpVUq1atfT1119fN+vV+vXrp/r160uSxo0bpyFDhkiSzp07p8cee0xHjx6V3W7XuXPnlJycrIYNGxaZp3nz5urUqZMkqUuXLpozZ06pPh/AWxSeuaJc+S8KFirUkCFSjx7uf9ySrl7VqVNH586dK3Ao7ezZs6pbt+4v/pury1dR25xOp2rVqlXsoHpQUNB1s+Xk5GjIkCHasmWLOnTooAsXLqhmzZrKzs4udcH6JSXJWpz8fTF+/Hj17dtXq1evls1mU1RUlCzL+sV/V7VqVdf7AQEBNzSLBngqBtpxNWawUKFCQqRWrdz/VtKCFRsbq0WLFrk+XrJkiZo2baoGDRq4bvv4448lXTnM9umnn+qee+4pdlvLli1Vo0YNLV682HW/lJQUpaenl2rfWJalnJwc10u433rrrQLba9SooczMzF/8+Gpdu3bVt99+6xp6z8vLU3p6eqmzJiQk6MyZM5KkDz/8UPfee68k6fz58woNDZXNZtO2bdt04MCBEuUCfFlycrL27NlDuYIkChb8zLx583Tq1CmFh4crMjJSf/3rX/W3v/2twH3q1Kmj6OhodezYUU8++WSBi7AWtc1ut2vdunVas2aNwsPDFRYWpscee0yXL18uVbYaNWpo1qxZ6tixo6Kjo1W5cuUC2ydMmKCxY8e6htwHDx6sxMRE15D71W6++WZ9+umnev755xUeHq6oqCj961//KnXW7t27a8SIEWrVqpW+//57zZ49W5L06quv6vnnn1dkZKT+/Oc/uw7/SVfmyFatWqV27dq5htwBX2ZZlnJzc9W+fXulpqZSriBJsjmdTqfpEPAe+/btU3R0tPbu3auoqKjr3u5tmjRpovj4+CLP21TcNl80c+bMAkP0vsrE1256erpq1apVIc/lLb6q21fVs87oUlA9dUhLMB2nxPIPC4aGhhZYHQeYwQIAoAwKz1wBV6NgAVe5+lQHpdnmi2bOnGk6AuCxGGjH9TCDBQBAKS1cuJByhWKxggUAQCk99dRTuvfee9WmTRvTUeChWMECAKAELMvSkCFDtHHjRgUEBFCuUCwKFgAA15E/c7V+/XoFBASYjgMvQMECAKAYDLSjLChY8CtNmjRRq1atClyqpX379tq6dau5UBVg9OjRJTqfVXx8vHbt2uX6ODU1VQsXLixwn759++rIkSPujlgqiYmJWrFiRYHbbDabMjIyzASCT4uLi6NcodQoWPA72dnZ+uijj0zHKDc3cp2/khSshIQEtWzZsszP4Q5FFSygvEyZMkUJCQmUK5QKBQt+Z+bMmXr55Zd16dKla7alpaVpyJAhatu2rdq0aVPgzMxNmjTRjBkz1KVLF912222aNWvWLz5Hz549FR8f7/p42LBhrusYjh49Wo8++qi6du2qFi1aaNSoUa5L1RS37b///a/Gjh2rjh07Kjw8XOPGjVNOTo7r+SZMmKAuXbqod+/exX7+mzdvVpcuXdSuXTuFhYW5ymZCQoI+++wzvf7664qMjNSHH36o8ePH68iRI4qMjNSAAQNc+yH/YtH/+c9/NGzYMLVt21bh4eGaPn36dbNebevWrWrTpo0efvhhtWnTRtHR0a7HPn36tGJiYhQdHa2wsDA9+eSTysvLU1pammbMmKEtW7YoMjJS48ePdz3eu+++q44dO+q2224rcL1FoLQsy9LkyZN14cIF3X777erZs6fpSPAynKYBFevs2Stv7hYSUuIrPkdERCgmJkZvvvmmpk6dWmDbU089pZYtW2rNmjVKS0tTdHS0IiIi1LlzZ0lSRkaGdu7cqbNnz6pZs2Z65JFHdOutt5Y67u7du7Vr1y5Vr15dgwYN0ptvvqkpU6YUu+2ZZ55R9+7d9cEHH8jpdGrs2LGaP3++nn32WUnS0aNHtW3bNgUGBhb73FFRUdqxY4cCAgKUnp6udu3aKTY2Vn379tWAAQMUGRmpiRMnSpKaN2+uiRMnukpPYSNHjlTv3r3197//XZL0008/SdJ1s17t0KFDmj9/vpYsWaJVq1Zp+PDhOnz4sIKDg/U///M/CgoKUm5urgYOHOja/tJLLyk+Pr5AiZWkKlWq6N///re++eYbdejQQQ899JDsdn7MoXSunrkaPHiw6/sfKA1+8qBirVkjvf+++x933LgrbyX08ssvq2PHjgVWPyRp06ZN2rt3rySpbt26GjJkiDZt2uT6ATtixAhJUkhIiJo2barjx4+XqWD95je/0U033SRJeuyxx7RgwQJXwfqlbfHx8dq5c6fmzp0rSbp8+XKBVzONHDnyuuVKks6dO6fHHntMR48eld1u17lz55ScnKyGDRuW6nPIysrSjh07tGHDBtdtderUkaTrZr1akyZNdM8997g+93HjxumHH35QSEiIJk+erB07dsjpdCotLU1t2rTR8OHDfzHTgw8+KElq1aqV7Ha7Tp8+XerPC/6t8EA75QplRcFCxRoyROrRw/2PW8LVq3xNmjTRiBEjij3MJ10ZnL5a1apVXe8HBATI4XDo66+/dhWvbt266Z133pHdbldubq7rvpZllep5itrmdDq1evVqtWjRosj7BQUFFfsc+caPH6++fftq9erVstlsioqKum6+0rpe1uLYbDbZbDbNnTtXaWlp2r17t6pWrapJkyZdN2dR/3+AksrNzeXVgnAbZrBQsUJCpFat3P9WyoIlSdOmTdPSpUt18uRJ12333nuvPvjgA0lXDnetWbNGvXr1KvZx7rjjDiUmJioxMVHvvPOOpCuH1nbv3i1JOn78uHbs2FHg3/z9739XVlaWcnNztXjxYt17773X3TZo0CD96U9/cpWG8+fPKyUlpdSf9/nz5xUaGiqbzaZt27bpwIEDrm01atRQZmbmL358taCgIPXo0UNvvPGG67b8Q4SlyZqamqotW7a4Pvd69eqpYcOGOn/+vOrXr6+qVavq9OnT+tvf/laiXEBZBQQEqE+fPpQruAUFC34rJCREEyZM0KlTp1y3LViwQIcPH1bbtm0VExOjqVOnqlOnTqV+7Oeee05btmxR27Zt9cILL1zzGB06dFBsbKxat26t4OBg18xTcdvefPNNVatWTZGRkQoPD9c999xTpgtQv/rqq3r++ecVGRmpP//5zwWyPfTQQ1q1apXatWunDz/8UOHh4QoLC1ObNm1cQ+5X++STT7Rnzx6FhYUpMjJSb7/9dqmzhoWF6eOPP1bbtm31yiuvaPny5bLZbHr66ae1e/duhYWF6aGHHipQQu+55x5lZ2crPDz8msO8QGlZlqXPPvtMkjRx4kTKFdzC5nQ6naZDwHvs27dP0dHR2rt3r6Kioq57O641evToAoPkJd3mi7Zu3VrsEH1FMPG1m56erlq1alXIc3mLr+r2VfWsM7oUVE8d0hIq7HnzZ662b9+ulJQUNWjQoMKeG76NGSwAgF8qPNBOuYI7UbCACpZ/PqzSbvNFPXv2NLp6Bf/F5W9Q3pjBAgD4HcuydPnyZcoVyg0rWAAAv2FZls6ePauGDRtq69atxZ4iBbgRrGABAPxC/mHB3r17y+FwUK5QrljBAgD4vMIzV1xCCeWNFSwAgE9joB0mULAAAD5t79692r17N+UKFYo1UgCAT8rOzlZgYKC6deum48ePKzg42HQk+BEKFtzq8OHDpiMApcLXrG/KPyx4xx13aO7cuZQrVDgKFtwiJCRE1atX18iRI01HAUqtevXqCinDBcPhma6euXruuedMx4GfomDBLRo3bqzDhw/r7NmzpqMYk5mZqZo1a5qO4XG8Yb+EhISocePGpmPADRhoh6egYMFtGjdu7Ne/pLiAb9HYL6hICxYsoFzBI1CwAAA+Y9KkSYqNjVVERITpKPBznKYBAODVLMvSsGHDtGPHDtntdsoVPAIFCwDgtfJnrhISEpSTk2M6DuBCwQIAeCUG2uHJKFgAAK80ZswYyhU8FgULAOCVJk+erM8//5xyBY9EwQIAeA3LsjRjxgxdunRJbdu2VUxMjOlIQJEoWAAAr5A/czVnzhwdOnTIdBygWBQsAIDHKzzQ3qFDB9ORgGJRsAAAHs3hcPBqQXgdChYAwKPZ7Xb16NGDcgWvwqVyAAAeybIsffHFF4qNjdWUKVNMxwFKhRUsAIDHyZ+5Gjp0qNLS0kzHAUqNFSwAgEcpPNBet25d05GAUmMFCwDgMbKduQy0wydQsFDA4sWLZbPZFB8fbzoKAD902Zmr9PR0yhW8HocI4ZKamqoPPvhAnTt3Nh0FgJ/JdubqkvNnBVeqrN27d8tms5mOBNwQChYkSXl5eRozZozeeustPfPMM6bjwEds2iQtWBAkh8N0Es8ReXaT+nz/rqo5LdNRPEa2M1fPn9+lc3k/6+OgWylX8AkULEiS5s6dq27duik6OrpE93/yySdVs2ZNDRkyREOHDi3ndN7h/PnzpiN4nAULgnT8eJ4CAmhY+X6X+p5qZ38v0SEkSdnOPD2fc0x78/6rdwKbyxFYTenp6aZjeZRatWqZjoAyoGBBycnJWr16tbZt21bif/P2228rKiqqHFN5J34QFuRwSAEBDlWubFdIiOk0niH45GUF5knOSgG6UNm/d0q2M1fPZu7T/rwsza3ZQS2r36pqv/s/fB/BJ1CwoO3btys1NVW33367JOn06dMaN26cTp06pbi4OMPp4AtCQqSEBNMpPERfyXEyV/Zb6vn9TtmyZYsODRyohPUbdffddys9PZ1yBZ9BwYLi4uIKFKmePXtq4sSJGjRokLlQAHxWTk6OAgMDFRMTo+PHj6t27dqmIwFux2kaAAAVxrIsDRgwQNOnT5ckyhV8FitYuMbWrVtNRwDgg64+Q/tzzz1nOg5QrljBAgCUu8KXv+EkovB1FCwAQLl77bXXKFfwKxwiBACUu8mTJ6tfv34lPtce4O1YwQIAlAvLsvTAAw9o7969qlKlCuUKfoWCBQBwu/yZq7Vr1+rChQum4wAVjoIFAHCrwgPtMTExpiMBFY6CBQBwq9GjRzPQDr9HwQIAuNXvf/97yhX8HgULAHDDLMvSrFmzlJ2drfbt21Ou4PcoWACAG5I/czV79mwdOnTIdBzAI1CwAABlVnigPSoqynQkwCNQsAAAZZKTk8Plb4BfQMECAJRJYGCgoqKiKFdAEbhUDgCgVCzL0q5du9SzZ0/Nnj3bdBzAI7GCBQAosfyZq4EDByo9Pd10HMBjsYIFACiRwgPttWrVMh0J8FisYAEArqtwuWLmCigeBQsAcF0ZGRk6ffo05QooIQ4RAgB+kWVZunjxourXr699+/apUiX+LgdKgoIFAChS/mHBjIwM7dy5k3IFlAIFCwBwjcIzVzabzXQkwKvw5wgAoAAG2oEbR8ECABSwZcsW7dixg3IF3AAOEQIAJEkOh0N2u12//vWvdezYMdWtW9d0JMBrsYIFAJBlWerfv79eeeUVSaJcATeIggUAfu7qmatOnTqZjgP4BAoWAPgxBtqB8kHBAgA/NmvWLMoVUA4oWADgx1544QVt3bqVcgW4GQULAPyMZVkaOXKkDh06pF/96lfq2LGj6UiAz6FgAYAfyZ+5WrNmjdLS0kzHAXwWBQsA/EThgfaYmBjTkQCfRcECAD8xcuRIBtqBCkLBAgA/MWHCBMoVUEEoWADgwyzL0uuvvy6Hw6EePXpQroAKQsECAB+VP3P14osv6tChQ6bjAH6FggUAPqjwQHtERITpSIBfoWABgI/Jzs7m8jeAYRQsAPAxgYGBatmyJeUKMMhuOgAAwD0sy9K+ffvUtWtXzZ8/33QcwK+xggUAPiB/5qpfv37KzMw0HQfwe6xgAYCXKzzQXrNmTdORAL/HChYAeLHC5YqZK8AzULAAwIulpaUpNTWVcgV4GA4RAoAXsixL2dnZaty4sZKTk2W38+Mc8CR8RwKAl8k/LGhZlrZs2UK5AjwQ35UA4EUKz1zZbDbTkQAUgRksAPASDLQD3oOCBQBeYv369dq+fTvlCvACHCIEAA+Xm5urgIAADRo0SEePHtWtt95qOhKA62AFCwA8mGVZ6tevnxYsWCBJlCvAS1CwAMBDXT1z1aZNG9NxAJQCBQsAPBAD7YB3o2ABgAeaMWMG5QrwYhQsAPBAU6dO1ebNmylXgJeiYAGAh7AsS4888ohSUlJUs2ZNdenSxXQkAGVEwQIAD5A/c7Vy5Ur9+OOPpuMAuEEULAAwrPBAe8+ePU1HAnCDKFgAYNjw4cMZaAd8DAULAAyLi4ujXAE+hoIFSVLv3r0VHh6uyMhIde/eXfv37zcdCfBpVm6uFixYoLy8PMXGxlKuAB/DtQghSVq1apWCg4MlSZ9++qlGjx6tAwcOmA0F+CgrL0/D9u3Ttu3bdc899ygsLMx0JABuRsGCJLnKlSRlZmbKZrOZC+OFNm2SFiwIksNhOolnOXvWdALPY+XmamhqqrZfvKh1GzdSrgAfRcGCy8MPP6wtW7ZIkhISEoq975NPPqmaNWtqyJAhGjp0aEXE82gLFgTp+PE8BQTQsArLzXXIbs9TenqW6SjGWZal0Xv3akdWlj5t106RkZFKT083HctjnD9/3nQEj1SrVi3TEVAGFCy4LFmyRJL0l7/8RZMnTy62ZL399tuKioqqqGgez+GQAgIcqlzZrpAQ02k8i92ep6efrswvCUm5ubkKrVZNz952m3rVqyexT67B1wl8BQUL1xg1apTGjx+vc+fOqXbt2qbjeJWQEOk6i39+Jz09y+9/aVqWpUOHDik6OloftG0rx8mTpiMBKGe8ihDKyMjQyat+4MfHx6t27dp+/0sRcIf8k4j26dNHWVkcJgX8BStYUGZmpu6//35dvnxZlSpVUp06dbRu3ToG3YEbVPgM7UFBQaYjAaggFCwoNDRU//73v03HAHxK4XLFea4A/8IhQgAoB//5z3909OhRyhXgp1jB8gG5ubkKCAgwHQOArqxcORwONWvWTEeOHFFgYKDpSAAMYAXLyy1atEg1atTQpk2bTEcB/F7+YcH7779fTqeTcgX4MQqWF/v444/1xBNP6PLlyxo0aJC++OIL05EAv3X1zNWzzz7Li0QAP0fB8lIrVqzQuHHjNHHiRNlsNo0fP179+/fXrl27TEcD/A4D7QAKYwbLS7Vo0ULvvPOOevTooXnz5mnOnDm6+eab1ahRI9PRAL8THx9PuQJQAAXLS0VFRSkqKkpHjhxx3TZ16lSDiQD/k5eXp0qVKmn48OHq0qWLQkNDTUcC4CE4RAgAZWBZlvr166ePPvpIkihXAAqgYAFAKeXPXH3xxRe67bbbTMcB4IEoWABQCgy0AygJChYAlMLzzz9PuQJwXRQsACiF6dOna+PGjZQrAMWiYPkAp9NpOgLg0yzL0rhx43TixAnVrl1b3bp1Mx0JgIejYHm5Zs2a6fjx46ZjAD4rf+Zq6dKlfK8BKDEKlpez2+28PBwoJ4UH2u+66y7TkQB4CQoWABTB6XTq/vvvZ6AdQJlQsACgCDabTY8++ijlCkCZULAA4CqWZWnRokVyOp0aPHgw5QpAmXAtQi90+fJlvffeezp69KhCQkLUtm1bRUREqGXLlrLZbKbjAV7r6pmrnj17qmXLlqYjAfBSFCwv9NBDD2nt2rXq2LGj0tPT9ac//Ul5eXmqWrWq2rRpo3bt2mnhwoWmYwJepfBAO+UKwI2gYHmhDRs2aMGCBYqLi5N05RdDcnKyEhMTlZiYqAMHDhhOCHgXLn8DwN0oWF4oJCREzZs3d31ctWpVtW/fXu3btzeYCvBelSpVUu3atSlXANyGIXcvceHCBdf7Y8eOVXx8vLkwgI+wLEsHDx5U5cqVtWzZMsoVALehYHmJ4OBg3X777frNb34jh8OhhIQE/fGPf5RlWaajAV4p/7Bg7969dfnyZdNxAPgYCpaX+OyzzzRq1Cjl5eVpyZIlOnHihKZPn646deqod+/eevbZZ7Vs2TIlJycrNzfXdFzAo109c7Vs2TJVq1bNdCQAPoYZLC/Rv39/9e/fX5K0bNkyDRgwwDXUvn//fm3evFlvvfWWcnJyVKVKFf4iB34BA+0AKgIFyws9/PDD2rlzp7p3767u3bu7bnc4HLyKELiO1NRUHTp0iHIFoFxxiNBLZGVlud53Op1F3sdutys3N1dPPfVURcUCvIZlWbIsS61atVJKSgrlCkC5omB5iTfeeENBQUHq1KmTbDabVqxYoY0bNyotLa3A/SzLYgYLKCT/sOCDDz4oSapSpYrhRAB8HYcIvcTo0aPVsGFDJSYm6quvvtJHH32kefPmyWazqW7duoqIiFDr1q315ZdfqlWrVqbjAh6j8MwVAFQECpaXCA0N1WOPPSZJWr16teLj4xUaGqoDBw7owIEDSkxM1LZt21SjRg3NmzfPbFjAQzDQDsAUCpYXOnnypOv93r17q3fv3gbTAJ5rxYoVlCsARjCD5WUcDodq1qypTZs2mY4CeKz8F4KMGjVKBw8epFwBqHAULC9jt9sVFBSkvLw801EAj2RZlvr376/ly5fLZrOpWbNmpiMB8EMULC80atQorVixwnQMwOPkz1xt2bJF9erVMx0HgB+jYHmhm266SRs2bNDEiRMLXAQa8GcMtAPwJBQsL/Tmm2/q1KlTWrBggerWrauYmBj9/ve/1yeffKKDBw9yHiz4pUmTJlGuAHgMXkXohdLS0nTq1CnXKRqSkpK0fv16LViwQA6Hg2sRwi9Nnz5dDzzwQIHLRwGAKRQsL9WgQQM1aNBAffr0cd2Wk5Oj5ORkJSUlGUwGVBzLsvTMM89o+vTpru8JAPAEFCwvdezYMf31r3/Vf/7zH91yyy2666671KNHD0VFRSkqKsp0PKDcXT1zNXz4cNWvX990JABwoWB5oS+++EK//vWv9atf/UoNGzbUjz/+qBdffFFhYWFavny52rRpYzoiUK4KD7RzWBCAp2HI3Qs9++yzGjp0qE6dOqX9+/frp59+UlJSkpo2bao777xTKSkppiMC5cbpdGro0KEMtAPwaBQsL/T111/rkUcekd3+/xcg27Rpo7Vr1yomJkYvvPCCwXRA+bLZbBoxYgTlCoBHo2B5oXr16hW4HuHVxo8fry1btlRwIqD8WZalv/zlL3I6nXrwwQcpVwA8GgXLCz344IP6wx/+oDNnzlyzjdMzwBflz1zFxcXp2LFjpuMAwHUx5O6Fpk2bpi+++EJt2rRRXFycevXqpfr16+u7777TlClT1LNnT9MRAbcpPNDOtQUBeANWsLxQ5cqVtXHjRsXFxem9997TXXfdpVatWqlv376qWbOm3nrrLdMRAbfg8jcAvBUrWF6qcuXKeumllzRz5kwlJSXp7Nmzaty4sVq0aGE6GuA2TqdT1apVo1wB8DoULC/RoEEDRUREKDw83PXf1q1by263KzIy0nQ8wK0sy9L333+vli1b6tNPPzUdBwBKjYLlJYYPH66kpCR9/PHHOnv2rGw2mwIDA9WqVSuFh4cXKF716tUzHRcos/zDgl9//bW+/fZbValSxXQkACg1CpaXePPNN13v51/oOSkpSUlJSdq/f79WrFih3NxcSVKdOnV0+vRpU1GBMis8c0W5AuCtKFheqKgLPR8+fFivvfaali1bpp9++slgOqBsGGgH4EsoWF4qLy9P//rXv7Ru3TqtW7dO33zzjYKDg3X//ferf//+puMBpXb06FElJiZSrgD4BAqWF8nMzNQ//vEPrVu3TuvXr1d6erruuOMO3XfffVq0aJG6du2qSpU48wa8i2VZCggIUHh4uI4dO6bq1aubjgQAN4yC5SV69uypL7/8Una7XXfddZf+8Ic/qH///goNDTUdDSiz/MOCISEhWrp0KeUKgM+gYHmJbdu2qWrVqnrggQfUpUsXhYeHq27duqZjAWVWeOYKAHwJBctLTJ06VUlJSfrf//1fLV68WDabTTabTU2bNnWdpiH/rWnTpqbjAsVioB2Ar6NgeYmXX37Z9f6FCxdcp2hISkrSgQMHtGHDBl28eFE2m01BQUHKzMw0mBYo3scff0y5AuDTKFheqEaNGrrzzjt15513Frj9u+++U2Jiog4ePGgoGVA8p9Mpm82mxx9/XD179lSrVq1MRwKAcsFLznxIs2bNNHToUM2cOdN0FOAalmVpwIABWrt2rWw2G+UKgE+jYAEod/kzV5s3b9ZNN91kOg4AlDsKFly//Fq0aKGIiAj16tVLKSkppmPBRzDQDsAfUbAgSRo3bpyOHDmiAwcOaODAgRozZozpSPARzz//POUKgN9hyB2qWrWq+vbt6/q4c+fOmjNnTrH/5qmnpJo1yzuZ9zh71nQCD7Vpk6afPKnH2rbVXXPmSNf5uvILfLEAfoGChWvMnz9fAwcOLPY+Bw58r4CAIN18880KDr65gpJ5ttxch+z2PKWnZ5mOYpxlWXr55Zf14vHjanTmjG4LCJDj5EnTsTyGIzdXeXa7stLTTUfxKOfPnzcdwSPVqlXLdASUAQULBcyePVspKSnavHlzsfcLDQ1VaGiLCkrlHez2PD39dGW//2FoWZZGjBihbdu2aWREhDoFBMheubIUEmI6msfIs9tV+emn/f5rpSjsE/gKChZc5syZozVr1mjTpk3XvSZccLCUkFAxubxFenqW3/9yKDzQ3m3OHDkuX75SrviCcclKT/f7rxXA1zHkDknS3LlztXz5cm3cuFHBwcGm48AL5eXlafDgwQy0A4AoWJD0448/6plnnlFGRoZiYmIUGRmpTp06mY4FL1OpUiUNGTKEcgUA4hAhJDVs2FBOp9N0DHgpy7K0du1a/fa3v9XYsWNNxwEAj8AKFoAyy5+5euSRR3TixAnTcQDAY1CwAJRJ4YH2xo0bm44EAB6DggWg1Lj8DQAUj4IFoNQcDockUa4A4Bcw5A6gxCzL0smTJ9W0aVP94x//kM1mMx0JADwSK1gASiT/sGCvXr30888/U64AoBisYAG4rsIzV4GBgaYjAYBHYwULQLEYaAeA0qNgAShWcnKy9uzZQ7kCgFLgECGAIlmWpcDAQLVv316pqakKCgoyHQkAvAYrWACukX9YMC4uTpIoVwBQSqxgASig8MwVAKD0WMEC4MJAOwC4BwULgMuiRYsoVwDgBhwiBODy1FNP6Z577lGbNm1MRwEAr8YKFuDnLMvS4MGD9c9//lOVKlWiXAGAG1CwAD+WP3O1YcMG2e0saAOAu1CwAD/FQDsAlB8KFuCn4uLiKFcAUE4oWICfmjJlihISEihXAFAOKFiAH7EsS5MnT9aFCxd0++23q2fPnqYjAYBPomABfiJ/5uqtt97S119/bToOAPg0ChbgBwoPtHfu3Nl0JADwaRQswMfl5ubyakEAqGAULMDHBQQEqE+fPpQrAKhAFCzAR1mWpbVr10qSJk6cSLkCgApEwQJ8UP7M1YgRI3Ty5EnTcQDA73BtDMDHFB5ov+WWW0xHAgC/wwoW4EO4/A0AeAYKFuBDLMvS5cuXKVcAYBiHCAEfYFmWfvrpJzVq1Ehbt26VzWYzHQkA/BorWICXyz8s2Lt3bzkcDsoVAHgAVrAAL1Z45spu51saADwBK1iAl2KgHQA8FwUL8FJ79+7V7t27KVcA4IE4ngB4mezsbAUGBqpbt25KTU1VzZo1TUcCABTCChbgRSzL0sCBA/XMM89IEuUKADwUBQvwElfPXN13332m4wAAikHBArwAA+0A4F0oWIAXWLBgAeUKALwIQ+6AF5g0aZJiY2MVERFhOgoAoARYwQI8lGVZGjZsmHbs2CG73U65AgAvQsECPFD+zFVCQoJycnJMxwEAlBIFC/AwDLQDgPejYAEeZsyYMZQrAPByFCzAwzz//PP6/PPPKVcA4MUoWIAHsCxL06dP16VLl9SmTRvFxMSYjgQAuAEULMCw/JmrN954Q8nJyabjAADcgIIFGFR4oL1jx46mIwEA3ICCBRjicDh4tSAA+CgKFmCI3W5Xjx49KFcA4IO4VA5QwSzL0tatW9WnTx9NmTLFdBwAQDlgBQuoQPkzV8OGDVNaWprpOACAcsIKFlBBCg+0161b13QkAEA5YQULqABc/gYA/AsFC6gAWVlZOn/+POUKAPwEhwiBcmRZljIyMlS/fn3t2rVLNpvNdCQAQAVgBQsoJ/mHBfv06aO8vDzKFQD4EVawgHJQeOaqUiX+lgEAf8JPfcDNGGgHAFCwADfbuXOnvvzyS8oVAPgxChYkSRMmTFCTJk1ks9mUmJhoOo5XysnJkdPpVExMjI4fP065AgA/RsGCJGnYsGHasWOHQkNDTUfxSpZl6cEHH9S0adMkSbVr1zacCABgEkPukCT16NGjVPcf//VTUt+a5ZTGu1i5uRq0b592pqdralaWtH+/6Uie4+xZ0wkAwAgKFsrE/t+TOvqvZN188826OTjYdBxjrLw8DU1N1fasLK0JDVWP7Gw5Tp40HcujOHJzlWe3Kys93XQUj3H+/HnTETwS+6VotWrVMh0BZUDBQplUqVZdLbq1Nh3DuDdTUrT94kWt69BBd9aoIXtgoOlIHifPblflp5/ml0Qh7I+isV/gKyhYKJOL9mApIcF0DOOey85W3+RkRUdHKz09nV8ORchivwDwQwy5A6VkWZaGDx+uPXv2qEqVKoqOjjYdCQDgYShYkCQ9/vjjatiwoX788UfFxsaqefPmpiN5pPyTiH722We6cOGC6TgAAA/FIUJIkhYtWmQ6gsfjDO0AgJJiBQsoodGjR1OuAAAlQsECSuj3v/895QoAUCIULKAYlmXp5ZdfVnZ2ttq3b0+5AgCUCAUL+AX5M1evvPKKDh06ZDoOAMCLULCAIhQeaI+KijIdCQDgRShYQCE5OTm8WhAAcEMoWEAhgYGBioqKolwBAMqM82AB/49lWdq5c6diYmI0e/Zs03EAAF6MFSxA/3/mauDAgUpPTzcdBwDg5VjBgt8rPNDOhYkBADeKFSz4NS5/AwAoDxQs+LWMjAydPn2acgUAcCsOEcIvWZalrKws1a9fX/v27VOlSvytAQBwHwoW/E7+YcGMjAzt3LmTcgUAcDsKFvxK4Zkrm81mOhIAwAfxpzv8BgPtAICKQsGC39i6dat27NhBuQIAlDsOEcLnORwO2e129enTR8eOHVPdunVNRwIA+DhWsODTLMtS//79XZe+oVwBACoCBQs+6+qZq86dO5uOAwDwIxQs+CQG2gEAJlGw4JNmzZpFuQIAGEPBgk+aMmWKtm7dSrkCABhBwYLPsCxLI0eOVHJysqpXr66OHTuajgQA8FMULPiE/JmrNWvW6KeffjIdBwDg5yhY8HqFB9pjYmJMRwIA+DkKFrzeyJEjGWgHAHgUCha83tNPP025AgB4FAoWvJJlWXrttdfkcDjUvXt3yhUAwKNQsOB18meuZs6cqUOHDpmOAwDANShY8CqFB9ojIiJMRwIA4BoULHiN7OxsLn8DAPAKFCx4jcDAQLVq1YpyBQDweHbTAYDrsSxL+/btU9euXTVv3jzTcQAAuC5WsODR8meu+vXrp8zMTNNxAAAoEVaw4LEKD7TXrFnTdCQAAEqEFSx4pMLlipkrAIA3oWDBI/3000/6/vvvKVcAAK/EIUJ4FMuyZFmWGjVqpIMHD8pu50sUAOB9+O0Fj5F/WNCyLG3ZsoVyBQDwWvwGg0coPHNls9lMRwIAoMyYwYJxDLQDAHwNBQvGrV+/Xtu3b6dcAQB8BocIYUxubq4CAgI0aNAgffvtt7rllltMRwIAwC1YwYIRlmWpX79+mj9/viRRrgAAPoWChQp39cxV27ZtTccBAMDtKFioUAy0AwD8AQULFWrGjBmUKwCAz6NgoUJNmzZNmzdvplwBAHwaBQvlzrIsPfroo0pJSVGNGjXUpUsX05EAAChXFCyUq/yZqxUrVujHH380HQcAgApBwUK5KTzQ3rNnT9ORAACoEBQslJsHHniAgXYAgF+iYKHcxMXFUa4AAH6JggW3sixL8+fPV15ennr37k25AgD4JQoW3CZ/5uqFF17Q4cOHTccBAMAYChbcovBAe1hYmOlIAAAYQ8HCDePyNwAAFETBwg0LDAxUo0aNKFcAAPw/dtMB4L0sy9KhQ4cUHR2tDz74wHQcAAA8BitYKJMcZ54GDRqkPn36KCsry3QcAAA8CgULkqRvv/1WXbt2VYsWLdShQwcdOnSo2PvPv/iNtm3bppUrVyooKKiCUnq21atXm47gkdgv12KfFI39UrTly5ebjoAyoGBBkvT4449r3LhxOnr0qCZPnqzRo0cXe/8jjgvMXBWyZs0a0xE8EvvlWuyTorFfikbB8k7MYEFpaWnas2eP/vnPf0qShg4dqieffFIpKSlq3rx5gftevnxZkjSiahMFBwdr3759FZ7XU2VmZrI/isB+uRb7pGjsl6JlZGTo0qVLql69uukoKAUKFvTDDz+oQYMGstuvfDnYbDY1btxYJ06cuKZgpaamSpI+to7p4+joio7q8aLZJ0Viv1yLfVI09kvRvvnmG0VFRZmOgVKgYKFUYmNjtXTpUjVp0kTVqlUzHQcA/EKrVq1MR0Ap2ZxOp9N0CJiVlpam5s2bKz09XXa7XU6nUw0aNNCOHTuuWcECAADXx5A7VLduXUVFRWnp0qWSrrySp2HDhpQrAADKiBUsSJKOHDmi0aNH69y5c6pRo4YWL16stm3bmo4FAIBXomABAAC4GYcIAQAA3IyChRIr7dne/cWECRPUpEkT2Ww2JSYmmo7jESzL0qBBg9SiRQtFRESoV69eSklJMR3LuN69eys8PFyRkZHq3r279u/fbzqSR1m8eLFsNpvi4+NNR/EITZo0UcuWLRUZGanIyEitXLnSdCSUAgULJVbas737i2HDhmnHjh0KDQ01HcWjjBs3TkeOHNGBAwc0cOBAjRkzxnQk41atWqWkpCQlJiZq0qRJfA9dJTU1VR988IE6d+5sOopHWblypRITE5WYmKjf/va3puOgFChYKJH8s72PHDlS0pWzvf/www+sSkjq0aOHGjZsaDqGR6latar69u0rm80mSercubPrJLX+LDg42PV+Zmama//4u7y8PI0ZM0ZvvfWWqlSpYjoO4BacaBQlUpqzvQOFzZ8/XwMHDjQdwyM8/PDD2rJliyQpISHBcBrPMHfuXHXr1o2zuBfh4YcfltPpVMeOHfXqq6+qTp06piOhhFjBAlCuZs+erZSUFL3yyiumo3iEJUuW6IcfftCsWbM0efJk03GMS05O1urVqzVt2jTTUTzOtm3blJSUpH379ikkJESjRo0yHQmlwAoWSqRRo0Y6deqUHA6H62zvJ06cUOPGjU1HgwebM2eO1qxZo02bNnGh2kJGjRql8ePH69y5c6pdu7bpOMZs375dqampuv322yVJp0+f1rhx43Tq1CnFxcUZTmdW/s/XwMBATZw4US1atDCcCKXBChZKhLO9o7Tmzp2r5cuXa+PGjQVmj/xVRkaGTp486fo4Pj5etWvXVq1atQymMi8uLk6nTp1SamqqUlNT1blzZ73//vt+X64uXryojIwM18fLly9Xu3btzAVCqbGChRJbtGiRRo8erdmzZ7vO9o4rr678/PPPdfr0acXGxuqmm27y++H/H3/8Uc8884yaNm2qmJgYSVKVKlW0e/duw8nMyczM1P3336/Lly+rUqVKqlOnjtatW8egO4p05swZDR06VLm5uXI6nWratKmWLFliOhZKgTO5AwAAuBmHCAEAANyMggUAAOBmFCwAAAA3o2ABAAC4GQULAADAzShYAAAAbkbBAgAAcDMKFgCPMn36dNlstiLfXnzxRbc9T2pqqmbOnKkLFy647TEBIB8nGgXgUQYMGKAffvhB77333jXbmjVrpjp16rjleT788ENNmDBBWVlZqlSJvzUBuBeXygHgUZKSktSrVy917ty5XJ8nMTFRYWFhlCsA5YKfLAA8RmZmpr7//nuFhYUVe79Vq1bp7rvvVt26dVW9enVFRkbq888/L3CfzZs3q2fPnqpdu7aCgoLUtm1bLVu2TJLUoEEDvfPOO9qzZ4/r8OMnn3xSbp8XAP/DChYAj5GUlCRJat26tRwOR4FtAQEBrgsjHzx4UMOHD1ezZs2Um5ur5cuXa9iwYUpJSdGtt96qHTt2qF+/fpo0aZKmT5+unJwcffXVVwoMDJQkxcfHq1evXnr00Uc1fPhwSdIdd9xR4px79+7V0aNH9cADD7jj0wbgg5jBAuAx3nnnHT355JNFbvvggw80ZsyYa253OByyLEs33XSTPvvsM913332Ki4vTsWPHtGHDhiIf6/Tp02rQoIHWr1+v2NjYUud8//339cQTT+j999/Xo48+Wup/D8D3sYIFwGMkJSXplltu0dq1a6/Z1qJFC0mSZVlauHChFi9erGPHjikrK8t1n6CgIElS9erV9cUXX+hPf/qTRowYoUaNGhV4rAMHDkiSIiIiiszx888/67vvvvvFnD169NCECRM0ZswY5eTkaPz48aX7RAH4PFawAHiMzp07Kzg4WOvXry9yu9PpVGxsrJKSkvS73/1OkZGRqlWrlj7//HP94Q9/0OnTp1WvXj1lZWVpxowZ+uSTT3Tu3Dl169ZNr7/+umtw/rXXXtMbb7yhM2fOFPk8qampuu2220qUuWrVqrpw4YLr8CMASAy5A/AQTqdTycnJxQ6479y5Uxs3btSKFSs0efJkxcbGqkOHDvr2229Vr1491atXT9KVlay5c+cqLS1NW7Zs0fnz53X//fe7HicxMfEXV68kqUmTJnI6ncW+TZo0SdWqVVN8fDzlCsA1KFgAPMJ3332nixcvFluwfvjhB0lSy5YtXbft3btXq1atKrIw2Ww23XXXXRo8eHCBoflvvvlGrVu3LnPWWbNmaeHChVq3bl2ZZrgA+D5msAB4hPy5qJ9//lm7du0qsO1Xv/qV2rZtq3bt2ikgIEC/+93vNHbsWB08eFDvvvuuqlSp4ipY48aNU0BAgGJiYlS3bl199dVXmjt3rqZOnep6vBo1aujLL7/U1q1bVbVqVUVHR5dqFWrUqFG6++671bVrVzd85gB8EQULgEfIP0VDUQPjgwcP1po1a9SiRQstXLhQL774oj7//HPFxMTob3/7m9q1a+cqWK1bt9aKFSu0YsUKORwOtWzZUgsXLtRDDz3kerzZs2fr8ccfV69evRQQEFBgUL4kGjVqdM3gPABcjSF3AAAAN2MGCwAAwM0oWAAAAG5GwQIAAHAzChYAAICbUbAAAADcjIIFAADgZv8XvfUkljq2meAAAAAASUVORK5CYII="
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot()\n",
    "plotPath(rand(A), \"Upper lattice path\", :blue)\n",
    "plotPath(rand(setdiff(omega,A)), \"Non-upper lattice path\", :red)\n",
    "plot!([0, n], [0,n], ls=:dash, c=:black, label=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working With Sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Representing Sets in Julia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "Set{Int64} with 10 elements:\n  7\n  4\n  9\n  10\n  2\n  3\n  5\n  8\n  6\n  1"
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Listing 2.6: Basic set operations\n",
    "A = Set([2,7,2,3])\n",
    "B = Set(1:6)\n",
    "omega = Set(1:10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "Set{Int64} with 5 elements:\n  7\n  4\n  5\n  6\n  1"
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AunionB = union(A, B)\n",
    "AintersectionB = intersect(A, B)\n",
    "BdifferenceA = setdiff(B, A)\n",
    "Bcomplement = setdiff(omega, B)\n",
    "AsymDifferenceB = union(setdiff(A,B),setdiff(B,A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = Set([7, 2, 3]), B = Set([4, 2, 3, 5, 6, 1])\n",
      "A union B = Set([7, 4, 2, 3, 5, 6, 1])\n",
      "A intersection B = Set([2, 3])\n",
      "B diff A = Set([4, 5, 6, 1])\n",
      "B complement = Set([7, 9, 10, 8])\n",
      "A symDifference B = Set([7, 4, 5, 6, 1])\n",
      "The element '6' is an element of A: false\n",
      "Symmetric difference and intersection are subsets of the union: true, true\n"
     ]
    }
   ],
   "source": [
    "println(\"A = $A, B = $B\")\n",
    "println(\"A union B = $AunionB\")\n",
    "println(\"A intersection B = $AintersectionB\")\n",
    "println(\"B diff A = $BdifferenceA\")\n",
    "println(\"B complement = $Bcomplement\")\n",
    "println(\"A symDifference B = $AsymDifferenceB\")\n",
    "println(\"The element '6' is an element of A: $(in(6,A))\")\n",
    "println(\"Symmetric difference and intersection are subsets of the union: \",\n",
    "        issubset(AsymDifferenceB,AunionB),\", \", issubset(AintersectionB,AunionB))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Probability of a Union"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "MersenneTwister(UInt32[0x00000001], Random.DSFMT.DSFMT_state(Int32[1749029653, 1072851681, 1610647787, 1072862326, 1841712345, 1073426746, -198061126, 1073322060, -156153802, 1073567984  …  1977574422, 1073209915, 278919868, 1072835605, 1290372147, 18858467, 1815133874, -1716870370, 382, 0]), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000  …  0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000], 1002, 0)"
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Listing 2.7: An innocent mistake with Monte Carlo\n",
    "using Random, StatsBase\n",
    "Random.seed!(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "'a':1:'z'"
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = Set(['a', 'e', 'i', 'o', 'u'])\n",
    "B = Set(['x','y','z'])\n",
    "omega = 'a':'z'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "1000000"
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10^6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mcEst1 \t \tmcEst2\n",
      "0.285989\t0.307454\n",
      "0.285156\t0.307374\n",
      "0.285396\t0.306971\n",
      "0.285263\t0.307478\n",
      "0.286037\t0.307811\n"
     ]
    }
   ],
   "source": [
    "println(\"mcEst1 \\t \\tmcEst2\")\n",
    "for _ in 1:5\n",
    "    mcEst1 = sum([in(sample(omega), A) || in(sample(omega), B) for _ in 1:N]) / N # 2 samples => wrong\n",
    "    mcEst2 = sum([in(sample(omega), union(A,B)) for _ in 1:N]) / N\n",
    "    println(mcEst1,\"\\t\",mcEst2)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Secretary with Envelopes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "MersenneTwister(UInt32[0x00000001], Random.DSFMT.DSFMT_state(Int32[1749029653, 1072851681, 1610647787, 1072862326, 1841712345, 1073426746, -198061126, 1073322060, -156153802, 1073567984  …  1977574422, 1073209915, 278919868, 1072835605, 1290372147, 18858467, 1815133874, -1716870370, 382, 0]), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000  …  0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000], 1002, 0)"
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Listing 2.8: Secretary with envelopes\n",
    "using Random, StatsBase, Combinatorics\n",
    "Random.seed!(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "bruteSetsProbabilityAllMiss (generic function with 1 method)"
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function bruteSetsProbabilityAllMiss(n)\n",
    "    omega = collect(permutations(1:n))\n",
    "    matchEvents = []\n",
    "    for i in 1:n\n",
    "        event = []\n",
    "        for p in omega\n",
    "            if p[i] == i\n",
    "                push!(event,p)\n",
    "            end\n",
    "        end\n",
    "        push!(matchEvents,event)\n",
    "    end\n",
    "    noMatch = setdiff(omega, union(matchEvents...))\n",
    "    return length(noMatch)/length(omega)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "formulaCalcAllMiss (generic function with 1 method)"
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formulaCalcAllMiss(n) = sum([(-1)^k/factorial(k) for k in 0:n])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "mcAllMiss (generic function with 1 method)"
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mcAllMiss(n, N)\n",
    "    function envelopeStuffer()\n",
    "        envelopes = Random.shuffle!(collect(1:n))\n",
    "        return sum(envelopes[i] == i for i in 1:n) == 0\n",
    "    end\n",
    "    data = [envelopeStuffer() for _ in 1:N]\n",
    "    return sum(data)/N\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "1000000"
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10^6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n\tBrute Force\tFormula\t\tMonte Carlo\tAsymptotic\n",
      "1\t0.0\t\t0.0\t\t0.0\t\t0.3679\n",
      "2\t0.5\t\t0.5\t\t0.4994\t\t0.3679\n",
      "3\t0.3333\t\t0.3333\t\t0.3337\t\t0.3679\n",
      "4\t0.375\t\t0.375\t\t0.3747\t\t0.3679\n",
      "5\t0.3667\t\t0.3667\t\t0.3665\t\t0.3679\n",
      "6\t0.3681\t\t0.3681\t\t0.3678\t\t0.3679\n"
     ]
    }
   ],
   "source": [
    "println(\"n\\tBrute Force\\tFormula\\t\\tMonte Carlo\\tAsymptotic\",)\n",
    "for n in 1:6\n",
    "    bruteForce = bruteSetsProbabilityAllMiss(n)\n",
    "    fromFormula = formulaCalcAllMiss(n)\n",
    "    fromMC = mcAllMiss(n,N)\n",
    "    println(n,\"\\t\",round(bruteForce,digits=4),\"\\t\\t\",round(fromFormula,digits=4),\n",
    "        \"\\t\\t\",round(fromMC,digits=4),\"\\t\\t\",round(1/MathConstants.e,digits=4))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An Occupancy Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "Plots.PyPlotBackend()"
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots; pyplot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "occupancyAnalytic (generic function with 1 method)"
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "occupancyAnalytic(n,r) = sum([(-1)^k*binomial(n,k)*(1-k/n)^r for k in 0:n])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "occupanyMC (generic function with 1 method)"
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function occupanyMC(n, r, N)\n",
    "    fullCount = 0\n",
    "    for _ in 1:N\n",
    "        envelopes = zeros(Int, n)\n",
    "        for k in 1:r\n",
    "            target = rand(1:n)\n",
    "            envelopes[target] += 1\n",
    "        end\n",
    "        numFilled = sum(envelopes .> 0)\n",
    "        if numFilled == n\n",
    "            fullCount +=1\n",
    "        end\n",
    "    end\n",
    "    return fullCount/N\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "(100, 1000, [2, 3, 4])"
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max_n, N, Kvals = 100, 10^3, [2,3,4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "3-element Array{Array{BigFloat,1},1}:\n [1.0, 0.875, 0.7407407407407407407407407407407407407407407407407407407407407407407407407407507, 0.6229248046875, 0.5225471999999999999999999999999999999999999999999999999999999999999999999999951, 0.4378156806211790208132229165608223678639773747226879371369540551067757286321699, 0.3665749237277169623553713297448685748139906688148201513501706581631966528255032, 0.30679788173756605829112231731414794921875, 0.2566975664075897743720600625259512738119349301793832033205371817029714338983751, 0.2147373231974006400000000000000000000000000000000000000000000000000000000001176  …  1.069657334526848513060724413172337700899371866204013344464186587033996001303631e-07, 8.941398126513886825738219454513648957514481632012323087638438017198550960585833e-08, 7.474224150732878615287973980617874647757151383973183924654413320436969764181841e-08, 6.24779443278709531559446203760038434553284555234278644256711434180745164111931e-08, 5.222606152684812158203195841632183110427671533836853856853108552830418359647833e-08, 4.365638321965613522245910822459449263312726840635975376296460902253426149228287e-08, 3.64928821778342097547786009178828278567694924090499102231329702620336305491688e-08, 3.050482329903532982469269897810473381246011550126660150106171313698725436110443e-08, 2.549933186315839278642517342465001185146487459645309206844409327251120752386315e-08, 2.131518121463356314581839279444056373862323309838124475464104801083612948715274e-08]\n [1.0, 0.96875, 0.9221155311690291114159426916628562719097698521566834324035970126505105928974287, 0.874759197235107421875, 0.8287692521471999999999999999999999999999999999999999999999999999999999999999862, 0.7847071164895857972259954480621183666927848000315749949920691984878476574325002, 0.7427272334102909405727274267906930288862938106858539087969733283331152045631963, 0.70284107526669091424941469836706886553656659089028835296630859375, 0.6650018142992007267129801666793068719668995452677686090862917087069854990804312, 0.6291371892527383662414505728000000000000000000000000000000000000000000000000542  …  0.006856273841870679839389629923882948740262777084010298425147882588298984792194757, 0.006484001299149685428651923399433516981783173464246993012031836518050943199485602, 0.006131941361807191238183385104109986753552605243879449347274280885182312197039787, 0.005798996628234320192559054784301371214301911072040793749347442268699445377644412, 0.005484129273972237682245912991120098918527572297747207189737217384559617603292256, 0.005186357817678830288310205938509584244474555176632999187257978682742913112736439, 0.004904754062613397909058621528094294927389144356866613199298563204557899619383066, 0.004638440204117194406377155742273499989160609034762918977418368694109961589599356, 0.004386586094083865326033986962257652595655098884558902615807497861904712792164721, 0.004148406653902102583643399639030050460958520902926953546065128320646211502190965]\n [1.0, 0.9921875, 0.9768836051414926586394350454707107656353198191332622059645379261291469796270951, 0.96000116877257823944091796875, 0.9427194306625535999999999999999999999999999999999999999999999999999999999999927, 0.9254151645663163903126220138750955966803663834817324826298856924618003901896396, 0.9082433888810553305721461898700847614761807972256976375321787223958393991126867, 0.8912775309946364404910542292089882164323782944969493113873681977565865963697433, 0.8745552587713002021915982262567857030690844268610306217078499479335123061479573, 0.8580965753857733279441124596561712384000000000000000000000000000000000000000314  …  0.1810793226937462318090339891393017894934957057348577345486192293760337547900849, 0.1776309031586087009805472204870925241785494001783383361517584195201035246461121, 0.1742481451322462054757902778352648363837795512314618356524214942707614566084341, 0.1709297988116877045521291653737762436333174162217822537072843375595278046461545, 0.1676746381566628325755199588548902801167631369377178341390821173107054021932617, 0.1644814604394566711568122977589713063382413605360035823221501866943407758949084, 0.1613490858031737030755827984951386740797296847708847905903042160937474013987371, 0.1582763568282632252169601739706374222339552003338456159202310094139651756600984, 0.1552621381071602629328677485037517256759027952388567996868054041239664840613691, 0.15230531582689785720775802501212492131617990234654842215506236019980408966694]"
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "analytic = [[occupancyAnalytic(big(n), big(k*n)) for n in 1:max_n] for k in Kvals]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "3-element Array{Array{Float64,1},1}:\n [1.0, 0.872, 0.731, 0.597, 0.49, 0.449, 0.369, 0.296, 0.268, 0.246  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n [1.0, 0.965, 0.922, 0.884, 0.826, 0.784, 0.755, 0.703, 0.668, 0.644  …  0.006, 0.013, 0.006, 0.006, 0.002, 0.01, 0.007, 0.008, 0.006, 0.002]\n [1.0, 0.989, 0.979, 0.954, 0.951, 0.921, 0.921, 0.865, 0.866, 0.859  …  0.175, 0.166, 0.163, 0.168, 0.145, 0.164, 0.174, 0.151, 0.158, 0.159]"
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "monteCarlo = [[occupanyMC(n,k*n,N) for n in 1:max_n] for k in Kvals]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAACPF0lEQVR4nOzdd3hT1RvA8W9K2XtzhZayyl4ptewyRESkIEOGgIC0QUBAfgKCYloFQUBkKYQiOEBAWQIqAooIooyELXvIaNgyZEPP74/Q2j3TJmnez/P0KdzccdKT3Lw54z06pZRCCCGEEELYjYejCyCEEEIIkdVIgCWEEEIIYWcSYAkhhBBC2JkEWEIIIYQQdiYBlhBCCCGEnWVqgDVkyBB8fHzQ6XTs2bMn0f2OHTtGw4YN8fX1xd/fn4MHD2ZeIYUQQggh0ilTA6zOnTuzdetWypYtm+R+BoOBkJAQjh49yqhRo+jTp0/mFFAIIYQQwg50jsiD5ePjw6pVq6hTp068xy5dukTFihW5du0anp6eKKXQNI2tW7dSsWLFWPveuXOHw4cPU6VKFfLkyZNJpRdCCCGESJqnowsQ19mzZ9E0DU9PW9F0Oh3e3t6cOXMmXoB1+PBh/Pz8aNCgAQULFmTHjh1cu3Yt3jlLlynNvr37MqX8Imk3b96kQIECji6GSAGpK9chdeVapL5cR5EiRdJ8rNMFWGkxa9Ys9Ho9GzZsYPv27Xzy4Tdc+Hc/JcsW5HKFW9wqf4vpW6fzaOcjBg8cjKZpji6yW0vPC1ZkLqkr1yF15VqkvrI+p5tF6OXlhdVq5dGjRwAopThz5gze3t7JHtuqVSveeecdvKq+BMDgOzk4+e1xugZ1ZfyP4/lg3AdYjloytPxCCCGEEE4XYJUoUQK9Xs/ChQsBWL58OWXKlInXPZiUWrUCyElJAi5fJseKnxlQegBhlcMA6DO7D99v/h6r1Zoh5RdCCCGEyNRB7gaDge+//54LFy5QtGhR8ufPz/HjxwHo378/QUFBBAUFceTIEfr06cPVq1cpUKAACxYsoGbNmvHOZ7FY8PPzw2w2o9fro7cvXQrdusGDpi0Zf+QgYRcvxjt26KihTJs4DavVislkwmAwSNdhJrh27Zo0jbsIqSvXIXXlWhxRX2fOnOHKlSuZek1XUqxYsRT1lKWGQ2YR2ktiAdaff0KDBnBs3mby9m+GdepULPnzExwczMQZE5nx9wyyFcjGliFbuHryaoLnEBlDPghch9SV65C6ci2ZXV9nzpyhatWq3LlzJ9Ou6Wry5MnDoUOH7BpkZYlB7nH5+Nh+/1U8kKCmTdEWLQKTCYBWjVrRvXd3Gk9vTMNxDRniPQSwBWsAmqZJS5YQQogs48qVK9y5c4eFCxdStWpVRxfH6Rw6dIiePXty5coVCbCSU6IE5MwJf/8NjB0LrVqhHT6M0Wi0BVAFNTrd6cS0j6bxFm8BEBwcDIDRaCQ0NNRxhRdCCCEyQNWqVaWnJhNlyQDLwwO8veH0aWBwS2jQAG3mTEL/+AN0OgBGDh1Jqzat6DW7F9eWXmPaJ9NoUr+JtF4JIYQQIt2cbhahvfj4PGnB0ung3Xexbt9OaO/e0bMHNU3j+cDnmRcyD4DPLnxGlZpVJMASQgghRLpl2QCrbNknLVgArVtjrVqVsIUL46VnqF+1Pv3f6M/xB8d56duXePj4YaaXVQghhBBZS5YNsKJasKxWK5bdu7E0aACA5eefsVgssVqywqeGsyp4FT+d+IngNcG48MRKIYQQQjiBLBtglS0LV67AzJkm/Pz8CJ4/H4DgkSPx8/PD9GRWYZRnKzzLFx2+4Iu9X/DWxrccUWQhhBBCZBFZNsCKStXQurUBs9lMeHg4AOEVK2I2mzEYDPGO6VGzBx+3/phJP02i9autJdu7EEIIIdIkywZYZcvaft+5o6HX66OnpuqPH0dfqlSig9mH1R/GKxVeYf389Xyx5Yvo7VarldDQUAm6hBBCiAzg4+PDnj17ALh37x7t27enS5cuPHjwINljf/nlF55++mmqVatG9erVGTlyJJGRkRlc4qRl2QDrqafA0/O/ge6apmEcORItWzZYuTLBY6xWKxaLhYaeDQF45+t3mLt2LlarFavVSlhYmARYQgghRAa6desWbdq0oWTJkixdupQcOXIke0zhwoVZsmQJf/31F2azmW3btvHll19mQmkTl2UDrGzZwMvrSaoGbAFW6IcforVqBcuWJXiMyWQbr2UIsXUfPv7uMYZ2BoaMGhKd6d1iscQaJC+EEEII+7h69SotWrQgICCAuXPn4uGRsjClbt26lC9fHoBcuXJRp04dTkenEnCMLJloNEqsVA1ROneGkBC4eBFKloz1kMFgICgoCIvFQnBwMDM+nUHo16Es+2oZy76yBWWS8V0IIURWZLXafuxN02w/KdG1a1eCg4OZMGFC9LbJkyezaNGiBPc3Go28+OKLsbZduHCBZcuWsXbt2jSX2R6ydIDl4wOHD8fZ2L49GAy2bsIBA2I9FHcdwkYBjdj4zEaen/M8Hhc9iFgUQXh4OHq9XhKSCiGEyFJMJggLs/95jUZIaXtE27ZtWbZsGQMHDsTLywuAESNGMGLEiBQdf/PmTdq1a8fIkSOpV69eGktsH1k6wCpbFn76Kc7GYsWgeXNbN2GcACuKpmn/rVuoafzy1i8EvB8AQI1aNWQtJyGEEFmOwQBBQfY/b2raI9544w1q165Ns2bN2LRpE97e3iluwbp16xbPPfcc7du3Z/jw4fYoerpk+QDLaoV79yBXrhgPdO4MgwbB5ctQvHi84zRNi9X9V7V4Vb7s9SWdD3Rm4r6JrPBfgYcuyw5fE0II4YZS05WXkYYPH46Hh0d0kJWSFqx///2X5557jueee4533nknk0qatCwdJUTlwjp7Ns4DL74ISsF336X4XB38O7B01lJWn1/NqA2j7FZGIYQQQsQ2bNgwhg0bRmBgIKdOnUp2/+nTp7Njxw5WrFhBnTp1qFOnDuPHj8+EkiYuy7dggW2ge6VKMR4oUQICA23dhP37p/h8nap1Ytpz0xi6biheBb0YEjDEruUVQggh3FXcWX9DhgxhyJCUfc6+/fbbvP322xlQqrTL0gFWmTLg4fFfqoZYOneGoUPh2jUoUiTF5xwSMISzN84ybN0wSucvTadqnexXYCGEEEJkCVm6izBHDlvC0QRTYbz4Ijx+jPWLL1Kdof3DVh/StUZXXl7xMlv+3mK38gohhBAia8jSARbYxmEl2IKladC4Mdbly1Odod1D58Hn7T+ngVcD2s1tx6ARg2IdL8vqCCGEEO4tywdYZcsmHGBZrVYsAQFY/vwTSH2G9pyeOVnZdSXFIovx6ZRP2Xd8X6xzy7I6QgghhPtyiwAroS5Ck8mE35QpBD9+DNgytPv5+WEymVJ0XqvVysm/TmJ4yraszqumV1mxZgU//fSTLKsjhBBCuLksPcgdbF2E58/Dw4eQPft/26OXxenSheCTJ1Odod1kMhEWI+Xt+UXn6bQo9oB3WVZHCCGEcE9u0YIVGQnnzsXermkaer0efbt2AOirVElVgGUwGDCbzYSHhwMwcuJIPHt7EvhOIKa5tlaw8PBwzGYzBoPBfk9ICCGEEE4vywdYUclGExzoDmjdu2MEtIMHU3Xe6ADtybI5XVt1ZelbS/nN8zd+u/8bQPTjsm6hEEIIkTQfHx/27NkDwL1792jfvj1dunThwYMHyR77xx9/RCcYrV69OgaDgfv372dwiZOW5QMsb2/b7wRTNQBaQAChdeui/fJLms4fc93CjlU7MrPNTBadXsSz/Z6VwEoIIYRIpVu3btGmTRtKlizJ0qVLyZEjR7LH1K5dm507d7Jnzx7279/PpUuX+PTTTzOhtInL8gFWrlxQsmTiLVgAdOoE339vW7QwlaLWLYwKpgY9PYjRbUaz3ns9v1xOW9AmhBBCuKOrV6/SokULAgICmDt3Lh4eKQtT8uTJQ/YnA60fPHjA3bt30el0GVnUZGX5Qe6QRC6sKJ06wTvvwPr1dllKfHyL8Vz49wJ9vutD8bzFebbCs+k+pxBCCJGhrFbbj72lYhXprl27EhwczIQJE6K3TZ48mUWLFiW4v9Fo5MUXXwRsS+20b9+eEydO0LZtWwYOHJj+sqeDWwRYiaVqiFalClSrBsuX2yXA0ul0zG03l8t3LtNxaUc2vbIJ/9L+6T6vEEIIkWFMJogxO95ujEZI4Uz6tm3bsmzZMgYOHIiXlxcAI0aMYMSIEcke6+Pjw969e/n333/p2bMnK1asoFu3bukpebq4RYDl4wO7diWzU6dOMHMmPHhgW2MnnTw9PFnaeSnPfPkMz3/9PL/3+x3for7pPq8QQgiRIQwGuzQyxJOK8chvvPEGtWvXplmzZmzatAlvb+8Ut2BFyZcvH926dWPRokUSYGW0smXhzBl4/BiyZUtkp06d4P33YdMmaN3aLtfNkz0Pa7qvocmCJrRe2Jpt/bah5ZeB70IIIZxQKrryMtLw4cPx8PCIDrJS0oJ1/PhxypYtS/bs2Xnw4AErV66kVq1amVTihGX5Qe5gC7AePUqma7lWLahQwdZNaEdF8xTlp54/8fDxQ55b9BzX71236/mFEEKIrGbYsGEMGzaMwMBATp06lez+v/zyC3Xr1qV27drUrVuXkiVLMnbs2EwoaeLcIsBKLhcWADqdrRVr1Sqs587ZdbFmr4Je/NTzJ87eOEvQ4iDuPrxrl/MKIYQQWcXp06epU6dO9P+HDBnC6dOnKVeuXLLHhoSEcODAAfbu3cvBgweZMWMGuXLlysDSJs8tAqyyZW2/kxzoDrYA6/JlrD/8YPfFmquXqM7aHmvZFbGLbsu78Sjykd3OLYQQQgjn4hYBVr58ULRoMi1YgLVMGSwlS2JZsgSw/2LNDb0asvyl5fxw7AdC1oSglLLLeYUQQgjhXNwiwIIUpGoATHPn4nfxIsGbNgG2xZr9/PwwmUx2K0ebSm34vP3nLNizgFEbR9ntvEIIIYRwHm4xixCgZEkrP/9swmo1JLqEjcFgIKhMGSzBwQRjW6w5I9YSfLnWy1y9e5Wh64ZSPE9xRjRKPr+HEEIIIVyH27RgFS5s5eTJpMdVaZqGvm9f9EWKABm7WPOQgCG80+QdRm4cyfzd8+1+fiGEEEI4TpYPsKxWKxaLhcePLQCYzcmMq8qWDe2FFzAWLIhWqpRdyxF3ZuJ7zd+jV7levDr8VeZtnme3awkhhBDCsbJ8gGUymfDz82Pp0mAAQkKSH1el9elD6I0baOfO2a0cVqs13sxEnU7H69Veh1/htaWvseHEBrtdTwghhHAlPj4+7NmzB4B79+7Rvn17unTpwoMHD1J8DqUULVq0oFChQhlTyFTI8gGWwWDAbDYzZkw4AGPHhmM2mzEYDIkf1LQplCgB336b7utHtaBZLLYWtKh/79mzB4vFwt49ewGo8qAKQVODWGtem2BrlxBCCOEObt26RZs2bShZsiRLly4lRyqWr/v444+pUKFCBpYu5bJ8gKVpGnq9nqZN9QAULpyCcVXZstlyYn37LaQzlUJUC1pwsK0FLWpm4rBhw2JtP/DZAe59eo8uo7rw24Hf7J6HSwghhHB2V69epUWLFgQEBDB37lw8PFIephw8eJBVq1bx1ltvZWAJU85tZhHWqKHh4WHk7t0UDljv0gVmz4adO+Hpp9N8XYPBQFBQEBaLheDg4OiZiR4eHkRGRsbaXqRkEUZsHEFIeAhAdKuXpmkZMtBeCCGEcCZdu3YlODiYCRMmRG9LyWLPDx8+JDg4mM8++4xsiS46nLncJsAqXVqjbNlQbt5M4QFR3YTffJOuACtucBQ1MzEuvV7P6tWrOTnjZPS2qNYto9FIaGhomssghBBCJMd6y4r1X/v3nGj5NLT8KWskaNu2LcuWLWPgwIF4eXkBpGix57CwMDp27EjVqlU5neyyLZnDbQIsAC8vOHMmhTtHdRMuWwaTJ9vWKkwHTdMwGo3xWqJibo9q7drw+wbeGvIWpXqU4quBX1G9fHWsVismkwmDIfE8XkIIIURamcwmwjaH2f28xkAjoc1CU7TvG2+8Qe3atWnWrBmbNm3C29s7RS1Ymzdv5syZM8yaNYtHjx5x8+ZNfHx82LlzJ8WLF7fjs0k5twuwklsuJxY7dROCLZBKqBUq7vao4Okt3uJB8QeMODCCTX6bOPnXScLCwggKCpIASwghhN0Z/AwEVQ6y+3m1fKn7zBo+fDgeHh7RQVZKWrC2bNkS/e+oRaMd3ZLlVgGWtzds3ZqKA+zUTZhaUa1aTTs2pdM3nWj6QVNCnpJxWUIIITKOlj/lXXkZbdiwYXh4eBAYGMimTZsoV66co4uUall+FmFMXl5w/jw8fpzCA+w4mzA1olq1WtRqQec7ndn//n5ef+11IGPWRxRCCCEcLarlKcqQIUM4ffp0qoMrHx8frl+/bt/CpYHbBViPHsGFC6k4qEsX28CtHTsyrFxJeW/Ee8xdMxfPDrbGxjmmOcnn8RJCCCGEQ7ldgAVw9mwqDrJj0tG00DSN4BeC+bjPxwAsv76cWnVqSfegEEII4cQkwEqOg7oJ4+r0dCdeGvgSP1/+mX7f9SNSRQIJr3EohBBCCMdyqwCrcGHImzeVARY4vJsQbC1ZSz9ZyqJXFrFw30JeW/saSqkE1zgUQgghhGO51SxCnc7WipXqAKtpU6xFi2J64w0My5c7tHuuW41u3Ht0j74L+3Lz9E1a5GoByOxCIYQQwpm4VQsWpDLZaJRs2bA2b07YH39gPX8+Q8qVGn3q9OH5G8+zZMgSQkJs6RtkdqEQQgjhPNyqBQtsAdb+/Snf32q1YrVasfj4AGBZtgw8PBzeUjTvvXlMaTSFqSumwhqi1ziU1ishhBDC8dyuBcvbO3VdhCaTCT8/P4KnTAEg+MMPnaKlSNM0PurzEYPbDwZgN7slwBJCCCGchNsFWF5etjxY9++nbH+DwYDZbCY8PByA8Pz5MW/f7jR5qMY8P4ZmrzTj08OfMun3SUnuKzMOhRBCiMzhll2EYMvoXr588vvH7QrU37qF/vr1TF06JymaprHp8028u+ldRm0chaeHJ8MbDE9w36gZh7KeoRBCuJ9Dhw45ughOKaP+Lm4bYJ09m7IAK4qmaRjffRftyy9h8WJ49tmMKWAahTUL4+Hjh/xv/f+4fe02j3c+xmAwoGnaf+PInsw0lBmHQgjhPooVK0aePHno2bOno4vitPLkyUOxYsXsek63DrBSQ9M0QsPCbMlGp0+H2bMhVy77FzCNdDodH7T8gEeRj3h3+bswl+iWKpPJRFhYWPS+wcHBABiNRkJDQx1UYiGEEJnB29ubQ4cOceXKFUcXxWkVK1YMb29vu57T7QKsvHmhSJE05MKK0r07vP8+/PgjvPiiXcuWXhcuXKBb0W7syb2HjWxk6oqpDGc4HTp0ICgoCIvFQnBwsMw4FEIIN+Pt7W33AEIkLVMHuR87doyGDRvi6+uLv78/Bw8eTHC/H374Ab1eT506dahRowZffPGFXcuRplxYUapWhTp1bN2ETsZkMlGvXj02Tt8IwKLxi/Dz82PVqlXo9Xr0ej1A9L8lwBJCCCEyRqYGWAaDgZCQEI4ePcqoUaPo06dPvH2UUvTs2ZPPP/+cPXv2sHbtWgwGA7du3bJbOdKUzT2m7t1hzRqwY5nsIe6MxxZDWkAI5K6fG3gyjsxojBVYycxCIYQQwv4yLcC6dOkSu3btih5k16lTJ86ePcvx48fj7avT6bh+/ToAN2/epGjRouTMmdNuZUltLqx4unWDe/dg1Sp7FckuNE2L1VI1qfckhnQYwlvb32LOrjm2cWShofECLFnLUAghhLCvTBuDdfbsWTRNw9PTdkmdToe3tzdnzpyhYsWK0fvpdDqWLl1Kx44dyZs3L//88w8rVqwgR44ciZ578ODBFCxYkI4dO9KpU6dky1K0aC7OnMnFtWvX0/Zk8uUjf0AA6ssv+bdt27SdIwPlypWLkSNHkjt3bt59+l3u3b/Ha9+/xu3bt+lbsy9gG6918eJF9u3bB8CWLVu4ceMGJUuWpFSpUhlWtn/++SfDzi3sS+rKdUhduRapL9dRpEiRNB/rdIPcHz16xLhx41ixYgVNmzZl586dBAUFsX///kSnUM6aNSu61SYlKleG69chR44i5MuXxoL26gXDhlEkMhLsPLUzvYoUKcKHH34Y/f857eeQK2cu3vz1TXLnyc1A/4HMmDEj1szCYcOGAZkzszA9L1iRuaSuXIfUlWuR+sr6Mq2L0MvLC6vVyqNHjwDbWKszZ87Em9WwZ88eIiIiaNq0KQD+/v6UKVOG3bt327Estt/p6ibs0gWUwvrZZ04/hkmn0zHtuWkMCxjGoB8GMWP7jPgZ6sPDMZvNTpOhXgghhHBlmRZglShRAr1ez8KFCwFYvnw5ZcqUidU9CP8FYlGZVY8fP86JEyeoXLmy3coSFdOlK8AqUQJatsS6dKlLjGHS6XRMbT2VNxu8ydB1Q1l8erHMLBRCCCEySKZ2EZpMJvr06cMHH3xAgQIFWLBgQfRj/fv3JygoiKCgIObOnctLL72Eh4cHkZGRzJo1y675O0qXBp0ufQGW1WrF2qABlifdbK6QHV2n0zGp1SSyZ8vO/9b/j0eRj+hVvle8mYVCCCGESJ9MDbAqV67MH3/8keBj8+bNi/539+7d6d69e4aVI3t2KFUqHbmwwGWzo+t0Osa3GE92j+yM2jiKh80fOnV5hRBCCFfkdIPcM0t6c2EZDAZbdnSDgeBduwifOxe9n59LtATpdDrCmofh6eHJO5ve4WHkQ4yBRnQ6naOLJoQQQmQJbhtgpTcXVnRXYM+esGsX+nz5UjWT0RmMDRxL9mzZGf3zaO49useElhPQ6XRYrVZMJlP0YtFCCCGESJ1MzeTuTNKdzf0J7cUXMebJg7ZpU/pP5gBvNX6Laa2n8eHvHzJs3TCUUpJ8VAghhEgnt23BilqPUCnbgPe00ry9CQ0OhiVL4NNPwdP1/qRD6w8ll2cuBiwZwNkjZ3ku73OAawzcF0IIIZyRW7dg3b0L167Z4WS9e8PFi7Bhgx1O5hiGegba/9uelW+sxBBiy4UVHByMn58fJpPJwaUTQgghXItbB1hgn25C6taFatXgq6/scDLHmW2czYSlE9AF2Zr0ZptmS/JRIYQQIg3cNsCyS7LRKDqdrRVr5Uq4edMOJ3QMTdN466W3mNJ7CgALLy+kas2q0j0ohBBCpJLbBlglS9ryYaUnF1YsL78M9+/D8uVYrVanXz4nKd0bdqfn6z2x3LLw/NfPc+v+LQCXf15CCCFEZnHbAMvDw5bR3S4tWABlykCLFvDlly4/C0/TNL6a8RUbBm7AYrXwzFfPcO3uNZd/XkIIIURmcdsAC+yXqiGKtV07LL/+iuXJYHeLxYLFYnHZgKSRdyM2vbKJ438fJ+C9AH7d9ivg+s9LCCGEyGhuHWClN9loXKYLF/ADgt96C8gas/D0mp5u97pxfOJx/vf6/4Cs8byEEEKIjOTWAVZULix7MQwZgrlNG8JLlQIgPDw8S8zCe2f4O3y36TuKdS0GwHtT38sSz0sIIYTIKG4fYJ0/D48f2+d8mqahHzYM/YULAOj1evR6vcvPwtM0jaBmQSwctBCAqaenokqpRJ+XDIYXQgjh7tw+wHr0yJYj1G5atkQrUQKjv7/LB1Zx1apYi5FjRlLOqxzNv2jOr6d/TXA/GQwvhBDC3bl1gBWVC8ue3YRky4bWpw+hx4+jFS5sxxM7nqZpfDj+Q34b8hsBZQJ4buFzrD6yOvpxq9UaPQAeUjYYXlq7hBBCZEVuHWD5+ABYmTzZzh/wffvCP//A6tXJ7+uC8uXIx9rua3nB9wU6Lu3Il3u/BMBkMuHn50dwcDCQssHw0tolhBAiK3LrAKtgQShQwMqKFXb+gK9SBRo2hM8+s985nUxOz5ws7byUPnX68MqqV/ho20cYDAbMZjPh4eFA0oP809LaJYQQQrgKT0cXwFGsVitWq5XChS3cvEn0B72mafYZO/Xqq9C/v63/MaovMovJ5pGN8HbhlMhbgjc3vMnFhhf58JkPox+PGuSfEJPJRFhYWPT/o1q9jEYjoaGhGVpuIYQQIqO5bQtWVHfW33+nvDsrVbp0gTx54PPP7XM+J6XT6fig5QdMaz2Nydsm0291P4qXLI7RaEwyUE1Na5cQQgjhaty2BctgMBAUFMSYMRZ++imY8PBw+6ZUyJ8funaFBQvgnXdsa/NkYUPrD6V43uK8suoVrty5wtK3l5Ine55E94/bUphUa5cQQgjharL2p34SNE2L9aFeq1YG5Kzq1w9On4ZNm+x3TicTcxZgj5o9WNN9Db+c+oVnv3qWa3evJXu8pmnJtnYJIYQQrsZtA6woNWpogBGlMuADvmFDqFwZ5s+3/7mdRNxZgM9VfI5fev/C4SuHaTy/MWduJJ0DQ9M0QkNDJcASQgiRpbh9gFWvngaEcvt2BnzA63S2Vqzly21pG7KQpGYBBpQJ4Pd+v3P30V0aftaQ/Rf3O7i0QgghROZy+wCrbFlbHHTqVAZdoHdvW7r4xYsz6AKOkVzOq8rFKrOt3zaK5y1OkwVN+PX0r1itVj788ENJxSCEECLLc/sAK2dOKF0aTp7MoAuUKgVt22a5bsKUzALU8mts7rMZ/9L+tF7Ymi+3fsmkSZMkwBJCCJHluX2ABVCuXAa2YAG8+ipWs5nQAQOyTHARd5JAYgtbF8hZgHnN59EiVwve+uotQJKKCiGEyPokwALKl8/gAKtNG6xFihBmMmW5oCIlswAXzFvAuhHrYI3t/3bPOSaEEEI4GQmwsLVgZVQXodVqxbJ/PxZ/fwAs27dnqdablMwCjNudqAvS0WJCC3r3651ZxRRCCCEylQRY2AKsS5fg9m37nzt6MPhPPwEQPHCg27XexO1OnNJ7Cn88/oNeG3tx9c5VB5dOCCGEsD8JsLB1EYItJ6i9xRsMXrGi2y4Jo2kaI0eOpHvD7mx6ZRPHrh6j4fyGnPwno2YYCCGEEI4hARa2FizImG7CeIPBjx9HnzOnWybW1DSNUaNGoWkaAWUC+OPVP1BKUX9eff4896ejiyeEEELYjQRYgKbZ0jVk5EB3TdMwvvMOWrFi4Ebdg0mpUKQCf7z6B5WLVab5F81Z9tcyRxdJCCGEsAsJsLCtw+zjk/EBVuj776OFhMCXX2bMgC8XVDRPUTb02sCLVV6ky7ddmPT7JJRSKT4+5lqIQgghhLOQAOuJjJxJGEtwMNy8CUuXZsLFXEMuz1ws7LiQt5u8zaiNoxiwdgAPHz9M0bFx10IUQgghnIEEWE9keLLRKD4+8NxzMGdOJlzMdXjoPBjXYhzzg+bz2ZbPqNKlCkdOHUl0/6TWQhRCCCEcTQKsJ6ICrFT0TqXdgAGwcyc8CQ7Ef/rW7cusxrM4ufIkbUxtOPVPwlFvcmshCiGEEI4kAdYT5cvDv//ClSuZcLHnn4cyZWSwexxRrVKelzwBuPX3LfShelbvWh1v35SshSiEEEI4igRYT0SlasiUbkJPT+jfHxYtwnr0qAzSfiJuq9SVJVe4PuM6L458kUX7FsXaN6VrIQohhBCOIAHWE5kaYIEtwLp3D+uCBTJI+4mEWqX+2P4HnXp2oufKnry76d14MwxTshaiEEIIkdk8HV0AZ1G4MBQqlEkzCQGrhwfWJk2wfPUVQPRgbU3T3DZYiPvco1qllvovRf+7ntE/j+bI1SMsaL+APNnzRB8TGhrqoBILIYQQCZMWrBgybSYhT7rDfv2V4PPnARmkHVPcVimdTsdbjd9i+UvLWXt0LU0XNOX8zfMOLqUQQgiROAmwYsjMAMtgMGDeuZPwYsUAGaQdU1SrVNyWvI5VO/J7v9+5dPsS/uH+7Di/w0ElFEIIIZImAVYM5ctnXhehpmno69VD37s3AHpvbxmknQJ1StVhZ/BOfAr50HRBU77e/7WjiySEEELEIwFWDOXKwZkz8Phx5l1TCw7G6OmJtn595l3UxZXMV5JNr2yia42uvLziZcb8PIZIFenoYgkhhBDRJMCKoVw5ePQIzp3LvGtqVaoQGhyMtmgRPHiQeRd2cTk9c/J5+8+Z9MwkJm6dSIclHbh5/6ajiyWEEEIAEmDFUr687XdmdRNGGzwYLlyAZcsAWcA4MXH/LjqdjhGNRrC2x1p++/s36s+rz7GrxxxcSiGEEEICrFjKlrX9zrRcWFGqVYNnnoEZMwBZwDgxif1dnq/0PNv7b+exeszT855m/YmEu1slcBVCCJFZJMCKIVcueOopBwRYAK+/jnX7dixffikLGMeRkoWdKxerzPb+22no1ZA2i9ow9Y+p8ZKSSuAqhBAis0iAFUdmziSMpW1bTIUK4ffKK7KAcRwpXdi5UK5CrO62mpENR/K/9f+j18pe3Hl4J0UBmhBCCGFPEmDFkZm5sGLJlg3DkCGYPT0JnzIFkNxYUVKzsHM2j2xMeGYCSzotYeXhlTSe35gPp3+YogBNCCGEsBdZKieOcuVgwwbHXFsbNgxtyhQ4ZhuoHXMxY3eW2BI6SelaoytVilXhxaUv8qXnl8z+bjaelzwJDg4mPDxcco4JIYTIUNKCFUf58rYJfXfuOODihQtDr15oK1difOcdCQDiSO3CzrVL1WZXyC78q/gzaM8gDmY7CPwXoMnfVwghREaRACuOcuVsv0+fdlABXn8d7dIlQqtVkwAgjsSW0ElKkdxF+KHHD4xoOIJpB6dRvUt1ChYtmIGlFEIIISTAiicqwHLIOCyA6tWhRQuYPt1BBch6snlkY+IzE/mm7zecrnOaDj90kHxZQgghMpQEWHE89RRkz+6gmYRR3ngDtm+HbdscWAjXllDOqy7Vu7AjeAcPHz+kXng9vjv8nQNLKIQQIiuTACuObNnAx8fBAdbzz0OVKvBkNqFIvcRyXlUrXo0dwTt4pvwzdFjagbd/fpvHkZm4+KQQQgi3IAFWAipWjJ7I5xgeHjB8OKxa5eCCuJ6U5LwqkLMAy7osY2LLiUz8fSJtFrXh8u3LjiqyEEKILEgCrARUqQJHjji4EL16QfHiMG2agwviWlKalFSn0zGq8SjW91zPngt70M/V88fZPxxRZCGEEFmQBFgJqFzZNsj9wQMHFiJXLhg0CBYsgCtXHFgQ15KapKQA1XJXo+e1npRSpWj6eVNmbJ+BUkrWLRRCCJEuEmAloHJlePwYTpxwcEEGDgSlYPZsBxfEdWiaFisRaXI5r6xWKx9P/JhPmnzC60+/ztB1Q+m2vBvH/z4u6xYKIYRIM8nknoDKlW2/jxyBqlUdWJBixaBPH5g1C0aMsLVqiRRJLimp1WqNHq8FsG/vPnrqe/JUzacw/mJk87XNANGPx80mL4QQQiQlyRasOXPmcPPmTbtd7NixYzRs2BBfX1/8/f05ePBggvvdv3+fwYMHU6lSJWrWrEnPnj3tVoaUKFUK8ud3gnFYYEvZcPkyLFzo6JK4lOSSkiY2VmvtjLXc+eQOFxdfjLVd1i0UQgiRGkkGWMOHD0fTNHr37s3mzZvTfTGDwUBISAhHjx5l1KhR9OnTJ8H93nrrLXQ6HUePHmX//v1MyeR0BRcuWMmXLxSLxQm6h3x9oX17mDoVIiMdXZosI7GxWtOmTcNsNvPJnE9sO7aDF6a8QK++vRxYWiGEEK4myQArIiKCSZMm8ddff9G8eXMqVqzIBx98wPnz51N9oUuXLrFr167o1qhOnTpx9uxZjh8/Hmu/27dv89lnnzF+/Hh0Oh0ApUqVSvX10sPWfRTGoUNOEGAB/O9/cOgQ/Pijo0uSZSQ2VqtOnTro9Xrq+9cH4L0e7/HL/V8I+j6Ivy7/Fe88MhheCCFEQpIcg1WoUCEGDRrEoEGD2LdvH/Pnz2fatGkYjUaeffZZ+vXrR1BQENmzZ0/2QmfPnkXTNDw9bZfU6XR4e3tz5swZKlasGL3fiRMnKFKkCB988AEbN24kd+7chIaG0rJly0TPPXjwYAoWLEjHjh3p1KlTSp97PBcuXODixYvs27fvSVl2sWnTDUqWLJnpQV4sVatSwM8PNXEitxo0cFw57OCff/5xdBFiyZUrFyNHjiRXrlxcu3Yt3vZOtTrRsl5L+v3Yj3of1aP5peZMfnNy9Ovh8OHDhIWFERgYSM6cOR31NDKEs9WVSJzUlWuR+nIdRYoUSfvBKpXOnz+vAgMDlU6nUzqdThUvXly9++676s6dO0ket2vXLuXr6xtrm7+/v/r5559jbTObzQpQX3zxhVJKKYvFoooWLaouXLgQ75xR+5rN5tQ+jQQZjUYFxPsxGo12OX+6fPONUqDUjh2OLkm6XL161dFFSJN/7/+rXpjyggJU2ylt1bHTx5TZbFbh4eEKUOHh4cpsNquIiAhHF9VuXLWu3JHUlWuR+nIPKUrToJTixx9/pHPnzpQvX57Dhw8zYsQItm3bxoABA5g5c2ayA9G9vLywWq08evQo+pxnzpzB29s71n7e3t54eHjw8ssvA1C3bl3KlSvH/v37UxE2pk3ccTkQzmefJZ5DKVN17GhLMT9hgqNL4nasVitHDhyhfcH2AGzYuoF6PeqlKKGpEEII95RkF+GJEyeYP38+X375JREREbRq1YpFixbRvn376K6++vXrU69ePbp165bkhUqUKIFer2fhwoX06dOH5cuXU6ZMmVjdgwDFihWjZcuW/PTTTzz//POcOnWKU6dOUTUT8iXEn4qvB/Q4xez8bNngrbegf3/46y+oVs3RJXIbJpOJsLCw6P8/WPWABzzAo54H3Vt3Z9H4RYSHhyeZb0sIIYR7SbIFq1KlSnz55Zf07duXkydPsm7dOjp16hQdXEWpUqUKAQEByV7MZDJhMpnw9fVl4sSJLFiwIPqx/v37s3r1asCWHmLy5MnUrFmTDh06YDKZKF26dFqeX5pE5VAqXVpzjlQNUXr1gjJlYOJER5fErSQ043Db9m30GtKLRZcXAVCpeiUJsIQQQkTTKaVUYg+uXbuW559/Hg8P50z4brFY8PPzw2w2R88Gs6fWrSF3btuay05jxgzbQtDHjkG5co4uTapdu3YtfYMGHSih19u8zfMY/P5gijUpxjf9vqGhV0MHl9J+XLmu3I3UlWuR+nIPSUZOK1as4O+//07wsb///pt+/fplSKGcReXKTpJsNKb+/aFwYZg82dElcTsJZYfvH9ifI8uO4OPlQ9MFTXl/8/s8jnyc4PGS0kEIIdxHkgHWF198weXLlxN87MqVK3zxxRcZUihnUbkyHD8ODx86uiQx5Mljy+4+fz7IB3WmSiw7fNlCZfm1z6+83eRtQjeH0vLLlpy7eS7e8VarVdY3FEIIN5FkgKWUik72GdexY8coWrRohhTKWVSuDI8ewalTji5JHAMHQs6c8PHHji6JeMLTw5Ow5mH80vsXTvxzglqza7Hy0EqA6DUPo9Y1jPq3BFpCCJF1xZtFOHv2bGbPng3YkoH26NGD3Llzx9rn3r17nD59mi5dumROKR2kShXb7yNHbCvWOJLVasVkMmEwGGwtKIMGwcyZtpmF0pfvNAJ9Atk7YC/9V/en4zcd6V+3P8V3FmfCuP/Sa0SldjAajYSGhjqopEIIITJSvADrqaeews/PD4ADBw5QuXJlihcvHmufHDlyULVqVV599dXMKaWDlC4NefPaAqx27RxblqjupaCgIFuANWyYrQVr1ix4913HFk7EUiR3EZa/tJzPdn/G0HVDKZWrFF+t+4p7Z+8RHBwsKR2EEMINxAuw2rdvT/v27aP/P3bsWMqXL5+phXIWOp2t5cqRA91t6yJaY3UvwZOcXcHBMH26bVZhvnyOK6SIR6fT0V/fn6Zlm9JjeQ/67ujLgKcGAMRaA1EIIUTWlOQYrAULFrhtcBXF0TMJTSZT4hnD33wTbt6EOXMcV0CRJN+ivmx7dRtvNniTWYdmUbZ9WR7n+W+WocwsFEKIrCleC9aQIUN488038fb2ZsiQIUkerNPpmD59eoYVzhlUrgy//OK46xsMBoKCgrBYLPG7lzQN+vSBSZPgtdds/ZnC6eTIloMJz0ygdcXWvLLqFZ5Z9Qyz2syiZ62e8bt+hRBCZAnxAqw1a9bw6quv4u3tzerVqxOdRQjuE2BdugTXr0OhQpl//bjL98TrXnr7bfj8c/j0UxgxIvMLKFKsmU8z9g7Yy+s/vk7vr3rzZfEveSH/C0Ccrl8JtIQQwuXFC7BOxchJcPr06cwsi1OqXNn2+8gRSMFqQBkmoSSXAPj4QL9+/7ViyVgsp1YoVyG+evEr7m24x7K3lrGRjYDMLBRCiKzGOdfAcSJR6RkcndE9sSSXgK0V68YN+OSTzC+YSJMZY2fw4+YfqdrPtoh58yHN2frnVgwGg4NLJoQQwh7itWCtWLEiVSfo2LGj3QrjjPLls62vfPiwo0uSBG9vePVV2/I5AwdC/vyOLpFIRlRXYLG8xfCf78+2R9vou6MvX5b+Eg3pIhRCCFcXL8Dq3Llzig/W6XQ8fpzwumtZiaNnEqbImDG25XNmzYLRox1dGpFCpZ8qjdFo5Jkuz/DmH2/SaH4jRjUahaGygQXzFvyXWJYEks0KIYRwWkmOwRI2lSvD5s2OLkUyvLxsC0FPmYK1UydMX38tH8QOlNJgKKrrF2Br1a1M+n0Sob+GsuznZRybcCzW7EKZcSiEEK4jXoBVtmxZR5TDqVWuDJ99Bo8fQ7Zsji5NbLE+yEePhnnzsE6fTtinn8oHsQOlJRjy9PCkb8W+lLtbjhFf2maEvr/kfYIvBZNNl42zZ88CMuNQCCFcQbwA69q1axQqVAgPDw+uXbuW7AmKuME6eJUrw/378Pff4Gx5V2N+kKNpWDt0wPL554B8EDtCkpn3U1AHJpOJsLCw6P+vmryKVZNXxdpHZhwKIYTzizeLsHjx4uzatQuAYsWKUbx48SR/3EHMVA3OIupDPOYHeWhoKH7ffEPwnTtAnKzvIlMkmXk/BQwGA2azmfDwcADGTBpDmUFl8OztSceRtgkl4eHhmM1mmXEohBBOLF4L1vz586lQoUL0v5NKNOouvL0hVy5bgNWmjaNLYxO3pSPqAz0kJAT/48cJ/uUXwqdNQ9+kibReZaIkM++nQNyWrk4tOzG21liMm4xMXj4ZgELlCslahkII4eTiBVivvPJK9L/79OmTmWVxWh4eUKmSc7VgJfVBbj1wAH75Bb3Fgn7oUEcX1a0km3k/FeeJSiybyzMXH7b6kMaFG9P3TF96rO/B0RxHGdFwBNmzZbdn8YUQQthJvAArITdu3GD//v1YrVY0TaNmzZoULFgwo8vmVKpUca4AK7kPcmOTJmhLl8IHH0Dp0o4ooltLNPN+Ko6PO76qXb12nFt9jtBfQxm7aSzLDy1nQfsF1CpZyw4lFkIIYU9JZnKPjIxkzJgxlClThqZNm9K1a1eaNm1KmTJlGD16tFvkwIpSubJzJhtN6INc0zRC16xBy5cP3nvPgaVzX0lm3k+HXJ65mPjMRP589U/uP7qP31w/Qn8N5cHjB3a9jhBCiPRJMsAaMWIEU6ZMYfjw4ezdu5cLFy6wd+9e3njjDT766CNGjhyZWeV0uOrVwWqFFEyszFSJfpAXLGhLPvrZZ3D0KGAbGB8aGorVanVASYU9+Zf2xxxiZnTj0YzfMh6/uX7sOL/D0cUSQgjxRJIB1ueff877779PWFgYNWvWpESJEtSsWZP33nuPsLAwPn+SDsAd1HrSC7Nvn2PLkSoDB4KmwdixwH8pHSTAyhpyeubkvebvsSt4Fzmz5aTBZw0Y/tNwbj+47eiiCSGE20sywHr8+HGiA3T9/PzcqovQ1xdy5nSxACtXLggLw/rNN1gWLYqV0sFisUiglUXULlWbP/v/ycSWE5m9azY1Z9fk55M/O7pYQgjh1pIMsDp37sySJUsSfGzJkiVZfqHnmDw9bd2ELhVgAfTujalYMfx69kxzbibh3KxWK+PeG0fP8j3Z/9p+yhYqyzOfPkPdbnX56+Rfji6eEEK4pXgB1ooVK6J/AgMDWb9+Pc2bN2f69OksXryY6dOn06xZMzZs2ECzZs0cUGTHqVUL9u51dClSydMTw8SJmIHw4cMBSVSZ1cTs+q1YpCI/9/6Zt+u+zZ6le2g8ozGL9y9GKeXoYgohhFuJl6ahc+fO8XY6f/48mxNY7bhv37707t07Y0rmhGrVgiVLnHNNwqRo/fqhzZ0L69YBac/NJJxLQsvyXL58GQCfuz4AlL9Xnh6zehBeM5z5PefjU8jHQaUVQgj3Ei/AOnXqlCPK4RJq14Z79+D48f+Wz3EJOh1MnIjWogXGzp0ls3sWkVg2/5jMJjMAu57dRfWb1Xm/+fsMCRiCp0eKUuAJIYRIo3h32bJlyzqiHC6hZk3b7717XSzAAmjeHO2FFwg1m6FwYUeXRthBQtn8vby8ADh79mysDP/5i+Rn5qGZvLn+Tb7a9xVzX5iLf2l/Bz8DIYTIupIc5B7TnTt3uHbtWrwfd2G1Wvnkk1BKlLC63kD3KJMnw5kzMHOmo0si7EDTtFjdvXq9ntatW9O6detY2/R6PZV8KjGjzQz+7P8nAAHzAnj9h9e5ce+G5EcTQogMkGSApZTi/fffp0yZMuTPn5/ixYvH+3EXUQOJK1Rw4QCrShV47TUYNw6ejNURri+xbP4JLdXzdOmn2Rm8k4+e/YgFexZQ9ZOqfPX7V5IfTQgh7CzJAOvjjz9m6tSpDBo0CKUUb7/9Nu+++y6+vr74+PgQHh6eWeV0mKhBxFEDiQsUsLBzpwvnkDIabWOy4qxzJ1xXQtn8k1qqx9PDkzcavMGmTpvwfeDLqC9HAbBuy7ok86NJS5cQQqRckgHWZ599RlhYWPSSOB06dMBoNHLw4EGqVq3K8ePHM6WQjmQymfDz84seQPzTT8FcuODH9OkumkOqWDFbZneTCf6SHEnuJG6A9P3i79k8ZjOssT3+9rC38fPz45PZnyR6vLR0CSFEyiQZYJ0+fZo6deqQLVs2smfPzvXr120HeXgwcOBAt1gqx2AwYDabo1vrQkPDATMBAS6cQ2rwYChbFkaMcHRJRCaKGyDFfW23eaMN2QZk4+scX7P+xPpYx8VsxZWVAIQQInlJBlhFixbl33//BcDb2zv6Bgtw5coV7ty5k7GlcwJxBxK3aaMne3Y958+7cKqDnDlh0iT44QdYvz75/YVLSyxAgtg50cb1HMe+9/bhXcab1rNb07B/Q8xHzfFacWUlACGESF6SAVajRo3YuXMnAD169CA0NJQ333yT0aNHM3z4cFq2bJkphXQGUYOGvbw0qlVzwSVz4urYERo3hv/9z5Y5VWRZyQVIMQfEVytejU2vbOJ9//c5svIIjac35mGdh/y548/oli5ZCUAIIZKXZLbB0NBQzp8/D8CYMWO4fv06ixcv5u7du7Rq1YqZbjTdP2rQMNgSjrp8gKXTwdSpWJ9+GlOHDhjmzpUEpFlUQvmy9Hp9dH3HfG1HZYcvdasUAI2zN2biuoks9VrK8Bq2pZZkJQAhhEhekgFW5cqVqfwko2bOnDmZPn0606dPz5SCObNatWDZMtdbMicef3+szz9P2Nq1BB05IgFWFqVpWqy6TSpAipsdfuP0jQA8DHrIoH8GUbVzVR7nlRZPIYRITooTjZ47d46dO3dGt2i5s1q14M4dOHnS0SVJu+hxOc2bA2B57z0ZuJzFJZYbK6a4A9+jugP/mP0HC3sv5J/6/xC4LJD3N7/P3Yd3M6voQgjhcpINsObOnYu3tzdly5alfv36eHt74+Xl5dYDXGvXtv125W7C6HE5T2YSBm/aJAOXs7ikcmPF3Cdudni9Xs9TTz3Fy7Ve5ujgowx+ejDv//Y+1T6txspDK1FKZdZTEEIIl5FkgDVhwgQGDBhA8+bNWbVqFX/88QerVq2iefPmDBw4kAkTJmRWOZ1KiRJQsqRtTUJXFa+lolQpzLVrYwgJcXDJhDPQNI2RI0fGC8by58zPpFaT2P/afqoUq0LHbzrSemFr/rosOdWEECImnUri6+dTTz1Fr169+PDDD+M9NmLECBYtWkRERESGFjApFosFPz8/zGZzpg+6ffZZyJMHVq3K1MvaXfTfcPZs9K+9Bl99BT17Ztj1rl27RpEiRTLs/MJ+kqsrpRRrj65l2E/D+Pv63wzyH4SxmZEiuaV+M5u8r1yL1Jd7SLIF6+bNmzzzzDMJPvbss89y69atDCmUK6hd27VbsKJEj8tp3x5eesmWfPTmTUcXS7gAnU5Hu8rt+PnFn2l+ujnzfptHpZmV+HTnpzyKfOTo4gkhhEMlGWC1bt2ajRs3JvjYhg0b3CoPVly1asHp03DjhqNLkj6xxuVMmWILrt57z9HFEi7k2uVrbFywkZVtV9KhcgcG/zCYGpNq0Htob5k0IYRwW/ECrKgszxaLhf79+7N48WL69u3Ld999xx9//MF3331Hnz59WLJkiVsnGqxVy/b711+z0AK4Xl7wzjswfbqsUyiSFTdD/JnDZxjkNYgfg34k191cfDXjK7ou6MqRK0fses0s834TQmRtKg6dTqc8PDyif3Q6XfRPQv93JLPZrABlNpsz/dr37yvl6anUqFGOK0OGuHdPqYoVlWrRQqnISLuf/urVq3Y/p8gYydWV0WhUQLyfkJAQNXfuXAWool2LqmwDsqlXF76qDp44qIxGo4qIiEj22hEREQnu68j3vDOT95VrkfpyD/ESjW7atCkTwzvXdfWqFR8fK3/88d/6bhA/qaPLyZkTZs2C556DRYsydMC7cG0JZYjfuXMnc+fOZe7cuQBcXXoVgK+OfsU327/h1sxbtGnbJtn3SNTC1EFBQWiaFp1hPuZ6ipAF3m9CiCwrXoAVGBjoiHK4HJPJxPHjYRw/bvt/1DpvRqMxetkRl9W6NXTrBm+8AW3aQNGiji6RcEIJZYhv27YtBoMhVtDl5eXF9XvX+XTjp/zGb7Sb2o6RnUfSo0EPnnrqqVjnTCyQWrRoEVOnTo3eL0u934QQWVKSS+VEOXjwIFu3bo2eWtq4cWOqV6+e0WVzagaDgX/+CWLOHAsPHsRf383lffwxVK1qm1U4f76jSyOcWMwM8QkFXatXr461/M7lJZcZsWQEM4Jm8O0n3xJQJiD6sbhL9UQFUsOHD8dsNie6nqIQQjidpPoP7927p7p06RI99ipXrlzRY69eeukldf/+/czqykyQo8djrFunFGThMSEmk1Kg1KZNdjuljD1wHempq5hjqCIiIpTZbFbh4eEKUOHh4eqTVZ+oqhOrKkJRL337kjpx7UT0cXH3NZvN0WOxHP2ed1byvnItUl/uIck0DWPGjOH7779nzpw5XL9+nbt373L9+nXmzJnD999/z5gxYzIhBHRetiVzNLp0SXp9N5fVvz80agQGA9y/7+jSCBcSM/1HQsvvDGw/kP0j9jM/aD5bz2ylyqwqDFs3jOwFsye4VE/U+ysl6ykKIYQzSDLAWrJkCRMmTCA4OJgCBQoAUKBAAYKDgxk/fjyLFy/OlEI6q1KloEwZjXLlkl7fzdVET4W/eBHmzoVTp8BNl0US9hM3OMrmkY2+dfty7PVjhDYLZf7u+VSYUYEJWyZQsFjBBAOplKynKIQQziDJAOvatWtUqVIlwceqVKnCtWvXMqRQriQgAP7809GlsK+oGVxWqxWqVYORI20B1uHDji6acGGJBUd5sudhTJMxnBhygj61+2D81UjgskC8grwoXrK4g0orhBDpk2SAVaVKFb766qsEH1u4cGGiwZc7CQiAXbvgURZYGSRu4siof1v79QNvb6x9+hBqNEqSR5EhiuctzvQ20zk06BCNvRvTf01/as2uxcpDK1GJL5kqhBBOKclZhGPHjqVLly6cPn2aTp06UbJkSS5dusSyZcv4448/+PbbbzOrnE4rIADu3IGDB6PGZLmuxGZwGY1GQufMwfrMM4Rt305Q+/bSRSMyTIUiFVjSeQkjGo5g9M+j6fhNRwJKBzDxmYk082nm6OIJIUSKJNmC1bFjR1auXMmdO3f43//+x8svv8zw4cO5c+cOK1eu5MUXX8yscjotPz/Ilg22b3d0SdLPYDBgNpsJDw8HIDw8HLPZTIcOHbAULoylcWMALBs22Fq2pCVLZCC/p/xY32s9G3ttJFJF0vyL5jy38DksVoujiyaEEMlKNMB68OABK1asoEaNGpjNZm7evMnZs2e5efMmZrOZdu3aZWY5nVbevFCjRtYIsBKa7aXX61m1ahV+fn4Eb90KQPBbb+Hn54fJZHJkcYWbaFm+Jdv7b2dZl2Wcvn4av7l+dPm2C4cuH3J00YQQIlGJBlg5cuSgR48enDlzBoC8efNSunRp8ubNm2mFcxX162eNACtK3Nle8Vq2APPYsW692LfIXDqdjk7VOnFg4AHmB81nx/kd1Jhdg77f9eX09dMpPo8sFi2EyCzJDnKPCrBE4gIC4K+/4OZNR5fEPuLO9orXstWuHfrp09Gywsh+4VI8PTzpW7cvRwcfZVrrafxw7Ad8Z/oy+IfBWG8lHzTFmiErhBAZKMkAa8KECYwbN45du3ZlVnlcUkAAKGWbTZiVRbdsTZoE+fJBSIjtiQuRyXJ65uT1gNc5OeQkoc1CWbR/EeVnlOfN9W9y+fblePsnOkP2ydqH0qolhLA3nUpi/nPNmjWJiIjg+vXrFC1alJIlS6LT6f47WKdj7969mVLQhFgsFvz8/DCbzdGtK44QGQmFC8OoUeA2ye2//x5eeMG2TmHfvik+LGo9S+H8MquurFYrJpMJg8GQptmpVquV6Z9M50GdB8w7No9IFcnQgKG82fBNCucuDEBoaGisGbJRjEYjQUFBTnEfSQ95X7kWqS/3kGSahnr16tntQseOHeOVV17hypUrFCxYkM8//zzJBaMXLFhAv379WLlyJR06dLBbOTKChwf4+2etcVjJatsWXnkF3ngDWrWCMmUcXSLhoqK67YKCgtIcYH04/kPMZjNvP/82k7dNZtr2aczaOYs36r/BsPrDMBgMBAUFxVos2svLCyBWqxYQb8FqIYRIiwQDrL/++os5c+Zw5coVnnrqKTp37kyrVq3SdSGDwUBISAh9+vRh2bJl9OnTh507dya47+nTpwkPD6d+/frpumZmCgiAzz6z9ZjFaOTL2j7+GDZssK1Z+OOPbvTEhT1Edc+lNcBJ7Pih1YfSw6cHhlADE69PZPr26bxR/w2GBgyNPlav17N69erE876FhtrlOQoh3Fjc1Z+3bNmicubMqTw8PFTJkiWVp6en8vDwULNnz07zitIXL15U+fPnVw8fPlRKKRUZGalKliypjh07Fm/fx48fq5YtW6pdu3apwMBAtXLlykTPazabFaDMZnOay2Yv332nFCh1+rSjS5LJfvzR9sRnzUrR7rKKvOvI6LoyGo0KiPdjNBrTfXzUvWHdb+vU0B+Hqpzv51SFJxZWI1eMVG+9/ZaKiIhQERERymw2q/DwcAWo8PBwZTabVURERIY+74wg7yvXIvXlHuK1YBmNRqpUqcKaNWvw8vLi5s2b9O3bl3feeYcBAwakKYg7e/Ysmqbh6Wm7nE6nw9vbmzNnzlCxYsVY+06dOpVGjRrh5+eX4vMPHjyYggUL0rFjRzp16pSmMqaXr68OKMzPP/9Lhw4PHFIGh3j6afK8+io533yTG3o9kZUrJ7n7P//8k0kFE+mV0XX10ksvERgYyL59+xg2bBjTpk2jVq1alCxZMkXrnCZ0fOnSpQHYsmULAIcth2lfqz0dWnfg63NfM/3gdHLny41ur47g2sH4+Phw48YNACpUqICPjw+Ay62zKu8r1yL15TrSNVYubsRVvHhxtXz58ljbTp06pXQ6nTpz5kyaorhdu3YpX1/fWNv8/f3Vzz//HGvb/v37Vf369dWDBw+UUsqlWrCUUsrHR6nhwx1dCge4fVupKlWUqltXqfv3k9xVvrm5jsyqq/S+j2Men1yr2Lkb59SQH4aoXONyqYITCirjJqP668Rfymg0umTLVRR5X7kWqS/3EC9Nw5UrVygTZ8By1GDQK1eupCmI8/Lywmq18uhJ3iSlFGfOnMHb2zvWflu2bOH06dNUqlQJHx8f/vzzT0JCQpg9e3aarpvZAgLcbKB7lDx5YNEi2L8fZOyKSKW4iW3Tc3xiyz1FJcUtXaA009tM5+SQk/Sr249Jv0+i/jf1edT0EdkLZrfbcxJCiHgtWDqdTu3YsSPWtkePHimdTqcsFkuaI7nAwEC1YMECpZRS3377rfLz80vRMa7UgjV1qlK5cin1pAHO/XzwgVI6nYpYsSLRFgH55uY6XLmuUnpvuHDrgnrzpzdV3vF5Vd7xedWbP72prLesmVRK+3HlunJHUl/uIcFEo82bN6dAgQLRP4UL23LJNGnSJNb2ggULpjiQM5lMmEwmfH19mThxIgsWLIh+rH///qxevTpdgaIzCAiAe/dg3z5Hl8RBRo6ERo2wDhwo2bKFQ6W0VaxkvpJMfnYyp4edZlj9Ycy1zMVnmg+v//A6Z27IKhZCiLSLl2g0oWR8STEajXYtUGo4S6LRKHfvQoECMH06DBzo6NJkPqvVitVsxtK5M8H37xMeHo5er4817V4S7LkOd6yr6/euM2vHLD5a/xE3f7/JS71fIiwoDN+ivo4uWpLcsa5cmdSXe0gyk7uzc7YAC6BePaheHb74wtElyXxJZcuOyiskNxbX4c51tXX7VprUb0KxYcW4Wugqnat1ZnTj0dTV6jq6aAly57pyRVJf7iHJtQhF6rnbQPeY67hFDzCeOxeA8Bw5MK9YET3A2Gq18uGHH0rXoXBaUYlLD+8/DMB7Vd5jdIXRbD+8Hf1cPc8vep4tf2/Bhb+XCiEyiQRYdla/Phw5YmXUKPdYPDZqmROr1Yqmaej1evRPcpjpS5VCHxaG9mQMn9VqZdKkSW7xdxGuyWQy4efnF53VfeCAgXzQ6wNeefQKizou4uzNszT9vCmN5jdi9ZHVRKrIdF9TFpsWImuSAMvOAgIArEyalLUHeUd904+5TInFYokOtIxGI9pnn8Hhw1gHDEh0XyGcSWJpHl4b8Bo9avZg34B9rO2+lmwe2Wi/pD01Z9fkiz1f8OBx2pMLx/ySIoTIOiTAsiOr1cqtWxby5Mn6gUTcb/rBwcH4+flhMpnQNI3Q0FC0Z56BadMwffFFovsK4UyiW2GfjOmM+nfUJA2dTkdb37Zs6buFrX23UqFwBfp814cKMyrw0baPuHn/ZoqvldSXFCFEFuDQJBHp5Gx5sNK7tporSfE6bpGRKuKFF5Q5Xz4VPmGCy6/55k7cOVdPREREirO7H7h4QPVZ1Udlfy+7KjChgBq5fqQ6f/N8ssfZ837hznXliqS+3IPMIrQjq9WK1Wpl/HgLK1YEM2dOOP7+sdMUZDUpqoMbN6BuXSx58uB38KDT1JdImsx0Sp3zN88zfft0TGYTd67dofqZ6nz89sc0r9k8wf2j7hcWi4Xg4OAE05qklNSVa5H6cg/SRWhHUd0LHTrYgoccOWJ3L2RFKUroWLAgLF2KduQIb2fxv4dwP1GD1D1uezCp1STODDvDoKqD2Lt0Ly1mteC5hc+x4cSGeDMPk+uOFEK4NgmwMkCzZho5chg5ciTr3yijx1sl96Hg7482bRrjLBa0TZsyp3BC2FlCM/5iDlK3Wq2c+OsENR7XAOBV7VVOHTrFs7OfpY6pDl/s+YL7j+7HOmd612IUQjgn6SLMIEFBcPs2/Pyzo0viRJTifrdu5FyzBv78E2rVcnSJRBKkGyO+mPccTdPidfHt3LmTuU/ywMXUe0hvrgZc5ftj31MqXykG+Q9iQL0BFMtTzC7lkrpyLVJf7sHT0QXIqlq0gLfesq1NmCuXo0vjJHQ6bn/0ETmPHIGOHWHXLihUyNGlEiJZMcdLgS3QihtMRc2SDQkJwd/fP8FxVYcuH2L69ul8sOUDxm8ZT8+aPRlWfxjVS1R3yPMSQmQc6SLMIM2bw/37toYaEUOePLBiBVy7Br16QWT6EzUKkdESSksyd+5cQkJC4uXMCg0NTXRcVdXiVZnzwhzOvnGWsU3H8v2x76kxuwatF7bmh2M/2CVxqRDCOUiAlUFq1oSiReGXXxxdEidUvjwsWgTffw/jxzu6NEIkK7EEpIkFU4mNq4oaw/XgxgPGNBnD6WGn+erFr7h29xptv25LlVlVmLl9Jrfu30qwHJL1XQjXIQFWBvHwsLViSYCViDZtwGi0/fz4o6NLI0SSkprxl1Awldjkj7hZ23Nky0HPWj3Z0X8Hv/f7nbpaXd746Q1KTy3NsHXDOHb1WJLHCyGcl4zBykDNm8PQobbB7nnzOro0TmjsWKxbtmB68UUMGzagNWni6BIJkaSkgqmkJDSGK+pYTdPQ6XQ09GpIQ6+GnLt5jk93fspc81ymb59Om4pt6OHTgyq5qrBn954EjxdCOCGHpjlNJ2fL5B7XoUNKgVLr1jm6JM4jbgZj86+/2urQ21upa9ccVCqREMk2bT9pydp+58EdtWD3AqU36RWB8Y+NebzUlWuR+nIP0kWYgSpXBk2TbsKERK/DdszWBWK5cgXLc89hPXvWwSUTwv4SG8NlMBgSPSZ39tz0qdOHXcG7WD1lNc9Nfo5s7bMB0GRwE77+6WsMBgNWq5UPP/xQug2FcDLSRZiBdDpbN6Hk1YzPZDIRFhYW/f/gO3dgxw6MnTsTun27A0smhP3F7cqLOZ4rOTqdjnb12tGuXjvW119P6+9aczj7YXr80YOG5xrSJm8bJk2aRNeuXaW7UAgnIi1YGaxFCzCb4fp1R5fEuST4jX70aAw7dmCdNElmSoksKb1Z22tWqInRaMT8PzPzms/j4bmHjF00FoC3F77Nml/XpPp9IzMThcggju6jTA9nH4OllFInT9rGYX33naNL4hzijcGKW4eDBilztmxOX6/uQMaJOLfExnVVeLGC+u7wd+rh44fxjomIiFBGo1FFREREb3OF+2hWI+8t9yAtWBmsXDkoW1a6CRMT8xu91WrF0rs3lkqVALD8+CMWiyX6m7V80xbiP3FbgT+Z8wnGRUbyN8xP+yXtKTe9HGG/hnHu5rnoY+Kum2ixWGLNbIz5fhNCpJOjI7z0cJVvXn37KlWrlqNL4RyS+uaW3EwrV6nvrEK+ZTtWQq1NCUnofWGOMKvg1cEq7/i8yiPMQz376bNq2oppymQyKUCFh4erkJCQVM9sFPYh7y33IIPcM0GLFrBgAVy+DMWLO7o0zstgMBAUFPTf4rn586P39sajdet437RBcgCJrC2qtSkoKCjJ17mmaYwcOTL2IHpNT1i9MIpuL0rhxoX5eNbHrP9+ffTjya2bmJHPyWQyYTAY5L0rsjzpIswEzZvbfm/e7NhyOLt42bI//RT9qVOs6ts33jpwfn5+mEwmRxZXiAyR2q47TdMYNWpUglnjJ46fyDMlnsESbuGrdV/ReFBj24PtoOH7DXn65aepVbsWEH/dxIx6bpKJXrgLCbAyQenS4Osr+bBSKnpcVsuWsHQphqNHMb/0EuFz5wIpyyEkhKtKaGHp1HyhSChAs1qttKzVkun9pgPwTtd3eFzqMf039ef5Nc/TsGdDbue4nTFPKJEyyXgvkeU5uo8yPVxpTM6AAUpVruzoUjhemsYezJmjFCjzsGEuU99ZgYwTcYyIiAhlNptVeHh49Hgps9mc5FismHWV1FjGuOO69l7Yq4b8MEQVnlhYEYpq9FkjNd8yX926f8uuzyktmeyzMnlvuQedUko5JLKzA4vFgp+fH2azOcVJ+xxl2TLo0gVOn7bNKnRX165do0iRIqk/cPRorBMnYurUCcPMmTJ+IxOkua6EXaTm/hazrmLOEIw7tiqx9829R/dYdXgV83fPZ+PJjeTNkZeu1bvSr24/GpRpgE6nS9dzSUuZsjJ5b7kHGeSeSWrVsuLhYWLhQgNvv+1+N5R0Gz8e7dw5QpcuheBg2xpEQmRhaU1Kmpas8bk8c9GtRje61ejG39f/5vM9n7NgzwI+2/0ZvkV96VO7D71r98bjtkeaBqmnJ5O9EK5KxmBlkn//tRIZGcbq1TLmIE08PGD+fGjdGjp2hG3bHF0iITKUpmmEhoamuYUnrQFa2UJlMTYzcnLoSTb22sjTpZ/m/d/ex3uaN50+60RYWBinz57O1DIJ4YokwMpgcQd37tplYfNmGdyZJtmzwzffgJ8ftG0L+/c7ukRCOK30BmgeOg9alm/JVy9+xe4euxldYTQXj10EoOWHLek0rROrd60mNaNM0lsmIVyJBFgZLO6MoMjIYJo1kxQDaZY7N6xZAz4+8OyzWP/8U7K7C5HBFn+xmPE9x3Pyy5MA3F1xlxVvrKD9m+2pPKsy434bx+nrpx1bSCGcjARYGSzuchZeXuG0bi0pBtKlYEFYtw7y58faubPk1REiA8Rcmiqhxdl37trJtx9+S528dXgv7D3KjStH0wVNCTeHc/3edccWXggnIAFWBoubPLNNGz1//qmnWDFpIk8Pa2Qklo8/xnLblrvHsnmz5NURwo5iJgWNlwRYr6eeXz06B3TmrTpv8fCXh7zv/z65s+dmwPcDKDWlFF2+7cLqI6t58PiBg5+JEI4hAVYmiRrc2amTxo0bsGWLo0vk2kwmE34vvEDw9esABA8fLtndhbCDpJKCxlucPcZ+pW6VYkK1CZi7mxnfYjzHrh6j/ZL2aB9pDFg7gK1nthKpIhO9pnT1iyzHwXm40sWVEo1GiYxUqkwZpYYMcXRJHMNeCfbiJWPMl0+Zq1ZVEYcOJXlMShbPFTaSDNF12LOuUpoUNCX77b+4X7214S3l/bG3IhRV9uOyavTG0Wr/xf2xzuWK9/L0kPeWe5BEow4waBD88AOcPAnpzN/ncuydYC/6NbB4MfpBg6BiRdiwAQoUSHxfF3u9OIokQ3Qd9qyrlCYFTU3y0EgVye9nfmfhvoV8+9e3/HPvH2qUqMELpV6gYeGGXDx+0a0SkMp7yz1IF6EDtG9vy+guWQbSL7rLIjDQFlgdPQpt2sCtW9H7yDpoQqRcQuOtEloEOqX7gS3lQ5OyTTC1M3HhzQus6b6GmiVqMmXWFIKaB6V43cXM7EqUbkuRXhJgOUBgIOTPD6tXy5s4vWLl1dHrYf16OHDAlifr33+B9C+eK4Q7SmlS0NQmD82RLQcv+L7A152+5vD8w4xfMp5awbVsD7aDeqH1yNsgL1fvXI13bMyB92mRmvtteq8lhARYDpAzp62RJSrAkjexHfn7w08/wZ498NxzcONGglPMzeakU2VI4CvcXUqTgqYneWiFshUY03UMCwYsAODdbu9SuHxhRm8fTamPStFmURsW7F7AoZOH7NIKnZL7rbR4C3uRtQgdpEkTK998Y+WXX/57E0P8NbtEGtSvDxs32pbVadUKbd06tBhjrlKyDlrUjTgoKEjqQ4gMFtUKZmhuIEwL49LtS3z222eYTCb6VeyHh9mDyF//m4EY1RptNBoJDQ1N9vwxx4tB0vdbk8lEWFhYmq8lRBRpwXKQs2dNgB8jRki3VYZ4+mn45RfbTIIWLeDy5RR1Zci3VyEyX9xWsBJ5S9C6ZGv+/u5v1rVfx3sj3qOOsQ60s+1fs39NQr8OpWvvrik6f2qGCaSlxVuIhEiA5SDDhhmoV89MjRryJs4wdevC5s1w4QIEBqJBsl0ZMl5LCMeK+yXn7JGztCnVhh8MP/Djmz8CkK10NsKOhlF7UW2eX/Q883fPjzdmK7lM9Indb1MzeF+IpEiA5SCapvHyy3oOH5Y3cYaqXh1++w1u3oSmTeHMmSR3l2+vQjhWUl9yaleqjdFo5AfDD5wbfo6prady++Ft+q/uT8kpJXnmy2eYvXM2F/69kGwm+uTut6kdvC9EXJIHy4FOnoQKFax06WJi+nSDW7yRMyv/i9VqxWQyYTA8+buePAktW8Ljx7aZhlWqJHm8q7+27EFy9biOrFRXqcmvFeXCvxdYdXgVyw8t55d9vxB5M5KKDypy/PPjTJg+gWcbPxt9bKz7goNkpfoSiZMWLAcqXx70eo1Hj9I2A0ckLt5sofLlYetWWwLSxo1h+/Ykj5dvr0I4Rlpam9QtxYW1F/jymS/5X47/wVw4/vlxAEYPHY2fnx8vj3mZa9muJfi+llnDIiNIgOVgvXrB2rVw7dp/2+TNnnZJDlIvXdrWXVi5sm3g+08/JXqe9Ew9F0KkX2q+5MT8QvXG4DdidfOHhIbQ6sNWbC+1nRqza1B5VmVGbhjJtrPbotdGlHQ5IiNIgOVg3btDZCQsXfrfNnmzp12yg9SLFLFlfG/eHF54ARYvdmBphRCJScmXnIS+UMUccwVgaGdg/cj1XH3/Kmu7r6WJdxM+3/M5jeY3oqSxJB2ndeTzHz+PPl5mDQu7cexSiOmTVRYIbdtWqQYNEljAODxcmc3mLLU4cUYvcpriv+GDB0q98opSoNT06RlaJlclC9K6Dnetq6QWnE5qcfdHjx+prX9vVQ16Nkjw+BFjRmRoud21vtyNBFhOYMkS2+f84MEpW8XelWXWjSVFr43ISBUxYIAygooIDlbq0aNMKZurkA8B1+GudZXeL6URERFq165dKvSjUAWo8r3LK0JQujd1qsn8Jmry75PVocuHVGRkpF3L7a715W6ki9AJBAXZxl57ekqKAHtJ0fgNnQ5rcDBhgHXePOjSBe7cybQyCpGVZcZY0vTmrNI0DT8/P9o1s2Uw/Xbot1g/sjK3+1wK5irI2E1jqfpJVXxn+TL8p+FsOrWJh48fZtjzEVmLLJXjBHLntn22r1ql8dFH/90YUrKki0hY1PiNxMRbOmPQIJg3D61hQ7R166BUqUwqqRBZU2YuN5XeWb8xjy+VrxT99f3pr+/PnYd3+OXUL6w5soalB5fy8Z8fUzBnQVpXbE3bSm1pU7ENxfMWt/OzEVmFtGA5id694fRp+P13SRGQGeINhp81C7979zCdOGFby/Cvv+IdI7M7hUieI5abSu+s38SOz5M9Dy/4voCpnYlzb5xjXft11DlchyOnjvDKqlcoOaUkDT5rwLjfxrHbuhuVgrSSVquVDz/8MNbfQ+4tWZSj+yjTI6uMwVJKqcePlSpbVqngYEeXJGM5y9iDRMdu7NypVM2aShUsqNRPP8U6JqHXW1IDaV2ds9SVSJ4z1VVSA89dXcx7gPWWVc23zFcdl3ZU+T7IpwhFaVM01W9VP7Xs4DJ1/e71ZM+R1Dbh+qSL0El4eNhyYs2cCdOn27oNRcaJmxU6Vnfs1q3QrRu0aQNTpmDt2hXrhQuxvpFHnSMzu0GEcAUGg4GgoKAEM7G7qnhDCp78fk57jr51+/Lg8QO2ntnKD8d+4IdjPzB/z3w8PTxp5NWINhXb8FzF5ygeWZwLce4jly9fBuDs2bOxzptU1nrhQhwd4aVHVov6jxyxzSZcutTRJck4zvRNW6kkWqAePVJq5EilQBlr107wG3lISIik1BBOwRnrKivdn1PbKnfqn1Pqgx8+UL4dfVXuUbkVoah8rfIleI7UnFeprN1qntXIWoROpn59KF4c1qxxdEkyhsutwbVwIdZXX8Xq64vllVcIHjGC8PBwdu7cydy5c+PtbjQakxxc70pcrq7cmDPWVbz1QF1YWtZHjPp8+mP7H9wudptvt3/Lhr0bOPnXSVgD5XuXp2GNhvg95Uee23kwGAypOm/Mz72s9LfOSqSL0Mn06gVDh8KlS1CihKNLI+jZE61yZbQOHeDDDwFbd2Lbtm0xGAxZqhtECHtKbiavK0lySEEccbsTD+w7gF6vx9jMyJzuc/hh8w+0XdOWctXKsTZyLQuPLyT/1fwAnM59mlblW6EVin8fSaybUoYqOC+ZRehkuna1jcdassTRJRHR/P1h1y40b2+MHh5oGzaglSqVrvw7QgjXk9gM75izAJNbrquub11GjhzJV72/4vKIy2zrt43+TftTpl0ZPtjzAT7Tfag8qzKDfxjM6iOruXn/JpD4MmChoaGZOmNTpJx0ETqhF1+Ekydhzx7Q6RxdGvtyxq6MFLt/H/73P/jkE1tT45w5WG/cyLJN8y5dV25G6sqxYn4WRbUoJdWdmFh9Xb93nV9O/cKGExv46cRPnLp+imy6bASUCSCgQADVc1fn0flHDDAMcJuhCq5Mugid0MCB8OyzsHkzNGvm6NKIaDlzwqxZ0LAhBAfDnj1oy5fLjUwIN5VQt11UMBW3dTslCuUqRMeqHelYtSMAJ66dYOPJjWw4uYHPT3zOP/f+Ifdl2xTzQ56HeGnQSwSHBLNn9554wZyMy3I86SJ0Qs88A9Wq2dI1CCfUowfs2AEPHkC9erBihaNLJIRwgKS6A+2RMLpCkQoY6hlY9tIyLo+4zI7+OxjacijlOpRj1qFZPLPyGZ7/5Xm+ufYNAMUrFKdu3bqxxmXF7S6UpKaZJ9MDrGPHjtGwYUN8fX3x9/fn4MGDsR6/d+8eHTp0wNfXl9q1a9OqVSuOHz+e2cV0KJ3ONtD9u+9sXYXCCVWvDjt32poaO3WC11+He/ccXSohRCYyGBJfPza92eXjyuaRDf/S/kx4cQInV57k+rjrbOy1kf76/lz2uAzNoMP3HfB6z4ugqUF88t0nQPxxWYkFXiIDZHZeiObNm6sFCxYopZT69ttvVb169WI9fvfuXfX9999Hr14+c+ZMFRgYmOC5slKelbhu31aqSBGlhg1zdEnsyxnz9aRLZKRSn3yiVM6cStWurdShQ44ukd1kubrKwqSuHCu1n0Upqa/U5ru6fve6WnNkjar/cv0Ec2u16dVGrfttXZbO3edsMjXAunjxosqfP796+PChUkqpyMhIVbJkSXXs2LFEj9m5c6cqW7Zsgo9l5QBLKaVGj1Yqf36lbtxwdEnsJ8t+EOzdq1TVqkrlyaPUvHkq4vx5l08GmGXrKguSunKs1AZDKamvtH6+RS0DNvWTqQpQzYc0V+VGlVM0SH1SU5E+mTrI/ezZs2iahqen7bI6nQ5vb2/OnDlDxYoVEzxm+vTptG/fPsnzDh48mIIFC9KxY0c6depk93I7SvfuOiZPLsQnn9zBYLjv6OLYxT///OPoImSMMmXgp5/IM2YMufr35+/mzQnbtInAwEBy5szp6NKlSZatqyxI6sqxcubMyZAhQwDbDMHkJFVfFy5c4OLFi+zbtw+ALVu2cOPGDUqWLEmpUqVSVBYfHx9u3LgBwNgOY6lduzYHuhzg5wM/s/HPjWybsw3aARp87vE5J789ScPSDWlYuiFl8pdJwTN2H+mZnevUswg/+OADjh8/zs8//5zkfrNmzcpSaRqiFCkCXbrAZ5/lZeTIvGTL5ugS2UeWnU5epAjWSZOwVq/OgbAwAE6sXUvBggVddm2xLFtXWZDUlWtJrL5mzJhB2JP7B8CwYcOA1KdeqFKlCkajkSpVqlCkSBGaFmlKU7+mtK/ZHr85fqwfuZ7rha+z+e/NbNi7ga8++QrqQdkyZQn0CaSpd1Oalm1KxSIV0WW1fEGZJFPzYF26dImKFSty7do1PD09UUqhaRpbt26N14I1ZcoUlixZwsaNGylUqFCC58uqebBi2r7dtnzOd99BUJCjS5N+WT1fT2hoaKybYxTjmDGEjh/vgBKlXVavq6xE6sq1JFVfaVmWJzUSSt8Q9Vk6ZdkUzuc7z29//8buC7uJVJGUzFuSJmWb0MTb9lM8sjjzwudJ+ocUyNQWrBIlSqDX61m4cCF9+vRh+fLllClTJl5wNXXqVBYvXpxkcOUuAgJsAdakSVYsFslp4uwMBgNBQUH/3Rxfegn9ypVo335rm22YRb8ICCHsIzXL8qT1/FEtYXHzeBX8pyDNyzVnRK0R5C2Sl21nt7Hl7y1sObOFkRtGcv/xffJeycvtWbe5VOYSXVp2IaBMADeuZN2Ey+mS2YO+Dh8+rOrXr68qVaqk/Pz81L59+6Ife/XVV9UXX3yhAFW+fHlVu3ZtVbt2bfX0008neK6sPsg9ypIlSkHWeK7uMhg31mvzwAGl6tZVytNTqffeU+rBg+j9Ujs4NjO5S11lBVJXriUjZhGmhdFoTPHA95NnTqrwteGqw8gOClC5O+ZWhKCyjcimqr9TXQFq0reTlPWWNcPK62pkqRwnZ7VaOXvWynPPWfjnH/s3F2c2d+nKiNcM/+ABvPceTJwINWrA/Pmg1zv1a9hd6iorkLpyLc5SX6npjkxs+EOjDo14qD1kx+wd0QPnvUp7EVgzkBq5a3B642nefuNtypR2v8HzksndyZlMJgIC/Pjnn4QXDhXOKV6SwRw5YNw4WwZ4wOrvj6VvXyzbtwMJJwOUbMtCiIwUtaRPQovWx70HJZRUNSQkhN9X/W4LrgDWAHPhqaNPcfjKYcasGsOcj+ZQeUJlWn3Vinc3vcu64+v4527mzHp19H1UWrCcXNQ3jK1bLQwdGkyLFuFMniwtWC7t4UNCW7cmbNOmeA9FzRRyhte21JXrkLpyLc5WX0kNfI97D0rJwtYeHh5ERkbyx44/GPzaYDqO7Mjlgpc5ePcg1zxtaSyqFKtCgzINqF+mPg3KNKBa8Wpk87DvVHlH30edOk2DiD/g8bff9BQtqscFYysRJXt2DIsWEfT771hGjiT41CnCAwPRh4biUahQdGsWEP3bVQNqIYTzS2rge9x7UMw1FhMbkB+3O3HFJNt6re+++y49B/fkz3N/8se5P/jz3J98ufdLHqvH5MuRj6dLP01A6QDbT5kA1C2VpsHzyT2HTOPYIWDp4y6D3JWyDXgcPdqoihaNUP36Obo0aSeDcWMz79hhew3nyaNUiRLK2LFjooNOM3tAvNSV65C6ci3OXF+pGfgeJe69KSqbfELL8sTd99/7/6pfT/2qJmyZoDos6aBKTSmlCEURiio1vJQC1LD5w9SWv7eo2w9up+j6aXkOGUECLBczbZpS2bIpdfiwo0uSNs58Y3GE6BuDxaJUly4qApTZ31+FjxsX78aU2a93qSvXIXXlWpy5vpIKjlIroXtWQttiBkiRkZFqx+EdauI3E1Wroa0UoLJ3yK4IQXm86aFqz66tglcHq7m75qrd1t3q4eOH8c5pz+eQHtJF6GIMBpgyBYxGWLLE0aURqZHQOIeYTfN88w3ajz+iDRoET7bpq1aNNc4BpNtQCJFx7JmHK2Z3YlLddlarlbCwMIKCgtA0je8Xfx+ri/HhqocAPB/8PCW0Evx57k8+2/0ZkTcjyXk3J2Vu22Yofrr6U57/53kCqgbEKnNyzyGhe3PM7anJoB9LpoZzduaOLVhKKRUerhQotWePo0uSes78zS2jpfj1evu2ihg6VBmzZVMRmqaMnTo5pLnbnevK1UhduRZXqC97D0lIrNsuJCQkXkvT7t27k22B+vf+v6rvsL4JnjNHyxyq+efN1WuLX1OdXuukth3cpiIjIxMtW2L35qjtaSWzCF3Qw4dQrRpUqQJr1ji6NKnjbLNnMkOal744eRL+9z+sq1ZhrVsXywsvEPz++/GOT+zbV3q5Y125Kqkr1+KO9ZXQfXDnzp3MnTs33r4pnU0d95zTPplGjjI5OPP4DEcfHmXn+Z2cvXkWgEK5CqHX9PhpftR7qh56TU/ue7m5eOFiorMgo7anOUxKc2jmBNy1BUsppb7+2taKtW2bo0uSOq7wzc3e0j3gcv16papWVWadzvZ6X78+1sMZ9T5wx7pyVVJXrsWd6yvm/Sq5sVIpbUVL6h548d+L6oejP6j3N7+vOizpoJ4KfUoRiOJ/qJwtcyZ4bw4MDIz1/7SSFiwXFRkJdepAvnxWWrUyMWCAa6wBJd/c0piN/+FDrBMmYJowAYOHB9ro0Vi7d8d640aGLQrrjnXlqqSuXIs711dqcm6l55yJibrWzFUzOR95nh1HdrBvzz6uLL0C7SCXVy58i/lSuUhl8l3Lx4L3FqS5BUsCLBe2ejW0b28BXOdv4M43Fru8Xq9dg/HjYdYsQnPkIOzff+PtEtW8nl7uXFeuRurKtUh9xZaaACmtQyKS+qJrtVrx8/Nj9urZ3C56m90XdmOxWsAKh8YfSnOAJbMIXZTVaqV0aSsVKlg4cQK2b5eZZc4u5oyaNCtSBD76CAYNwvDGGwStXo2ldGmCz58nfO5c9H5+iZ4/o8ZqCSFEesSaTZ2MuDMOU8pkMsWamRgcbFt+zmg0YjAYMBqNtK/XPtY5z54/y2een6X4GvGkuXPRCbjzGCxnSaSWWu489iBD/PmnMteta3sf1Kql1K+/KqUSHruQ2veL1JXrkLpyLVJfqZfe3FaOyI0liz27qLgLb3p6hrNypRmDweDgkolMFRCAtnYtxpdfRgNo1gxat8b688+EhYXFahaPmX8m5sLSQgjh7EwmE35+ftEtT8HBwfj5+WEymVJ0fFILW2cU6SJ0UXG7AosX1/PJJ3rat3dgoYRDaE89RejChaAU1vnzsY4bh6VXLwAsK1ey89KlWFOhYzaN22OslhBCZDSDwUBQUFCCY6hSwy5DNVJIAiwXF/ViqVRJo2dP+PprePllR5dKOIROh+nsWcJOn47eFDxuHAAhHTrg37Ztum5MQgjhKPbKMJ+a8V7pJV2ELi7qxfLyyxovvQRvvGGbaCbcU9yu4/BXX8Vctiyhq1ahX7AAiN80brVaCQ0NlS5DIYTTy8wWqPSSACsLmTYN7t+Ht95ydEmEo8QbZzBwIPoTJ9C+/hrt0iWMgDZ4MPz0EzyZehw1K0cCLCGEs4tqVJAAS2QqTYMJEyA83Eq/ftIi4c5ifcvLlg26d0c7fJjQ5cvRHj2C557DWrMmlgkTsOzcCcjgdyGEsCdJNJrFPH4MdepYOHDAj61bzTRq5Fx/F0mw5wSUgk2bCA0OJuzkyXgPRw1+l7pyHVJXrkXqyz1IC1YWYrVa2bvXwksv2abjjxolLRIiATodtGiBYetWzAsXEl6vHgDhuXNj7t0bwwsvOLiAQgjh+iTAykKi8oS8+65tGv7vv6cuT4hwL5qmoX/5ZfRPXh/6Ll3Qf/cdWv36WNu3Z9KQIVgjIhxcSiGEcE0SYGUhcWeQNWgQTq5cZlq0kOSjInHR47UmToRz52DaNKwWCx8uWoS1ZUuYPx/u3nV0MYUQwqVIgJWFxJ1BNmmSHh8fPYMHa9y9K9PxRcJizsqx3rqFpWFDLGPHAmDJnh3Lq69i1TR48004ccLBpRVCCNcgAVYWFNUiUaGCxjffwLFjMHSoTMcXyYtejuLJkkvB+/fjB5iqVcM6bx6hFStibdoUvvkGHjxwbGGFEMKJSSb3LChmplpNg/fftzJihBX4by26qP1cIZeIyDxJLUdhPXWKsEaNCLp5E61rVyhWDF55BYKDoXJlRxddCCGcirRguYFbt0yAH+HhaVskU7iPhBZE1TTNtmD0X38BYBk8GMu332J98UVYsACqVIHGjWHePLh505HFF0IIpyEBlhsYMMDAb7+ZKV3aNvh90qRwzGYzBoMMfhcJ0zSNkSNHomlawqvYd+mC6amn4Px5WLwY8uWDkBAoVQp69sS6dCmh776b5u5oGS8ohHB1kmjUjaxbZ6FNGz8qVjSza5eeggUzvwySYM91RNWV1Wq1tWAl0G0Y1bplMpkwtGuHtmEDfP45liNH8APMvXqhf+MNqFPHln8rheS9nTryvnItUl/uQVqw3Ejt2hoDBxq5dEnjxRdt6xZKS4FITkLdhnEXiw4LC8Oq02F95RUsixZhGT0aAMvKlVj0eqy+vjB+PCSQOT6mqEAuapygLN8jhHBVEmC5EU3T+OSTUL7/XmPbNtv45PPnZWahSJm4q9gnFAyFhobiV68ewRMmABD877+2WYi5c2MdP57QChWw1qkDU6bA6dPR54kK8hPsjpTxgkIIFySzCN1Q48bwySdW+ve3cu2azCwUKRNzdirYUjqEhYVF/z8qKAoJCcHf3z/+LMSTJwlr3JigIkXQxo6FESPg6aexNmhA2PTpBAUFJTmLUQghXIm0YLmps2dtMws3bJCWApE2cVcOCA+3TZ4IDQ1NeBbioUMAWHr0wLJ+PXsmTMCSNy+WTz6xbW/fHmtYGNr16+jr1o0+XgIsIYQrkhYsNxXVUhAaamHNmmCefz6c997T89RT8kEmUiZua2fMcVpAdHdiYi1dgYGBbN68+b/t586ByYTRZMJQtizGBg3QTp+GWrXAU25VQgjXIrMI3VzU3xDMDByoZ+ZM8MjAdk2ZPeM6UlpX0bMIDYYEW5oSm4Xo4eFBZGRk7O01aqCdPIm2eTOsXg0XLkChQtCmDbRrB889B4ULZ8CzdW3yvnItUl/uQboI3VzUwOVJkzRmz7YNfH/0yNGlEq4k5lqGiT2e0CzEOnXqxN9evz5ajx5gMtlybO3cCUOGwJEj0KMH1mLFCPX2xvrWW2A2Q2RkqmbCJrRvSrcJIURqSIDl5qI+HEeM0Fi8GJYsgc6d4d49R5dMZDVxZyEmtx0PD6hXD8LCbMHUuXNYR48m7OxZrDNn2h7TNKz9+9tmwu7fn2wZElqPM6XbhBAiNWRgg4jWtSvkzw+dOsHzz8O330LRokkfk1z3kBBR4s5CTG57TFarFevFi1h8fACwfPQRl2/dgl27OLttm21b69ZQsSJa8+YQGIjp4EEMr78ePcg+qpsSbF3jly9fBuDs2bNJbosqo7y+hRCpIWOwRDxbtkDHjpA3LyxfDn5+ie+b2jqQsQeuw5nqKjQ0NNZA+aQY8+cn6NYtWyb52rXRv/ACoWfOEPbVV2m+vtFoTDYIdCRnqiuRPKkv9yAtWCKeihWt9OxpYtMmA40aaXz6KfTrF3ufhFoEQL7pi4yRUH4sLy8vwNbaFHebZedOGDsWi6cnfPIJHa5fJ8jTE0u5cgQfO0b4sGF4NW0KefIkeHzMbXGz1kuLrRAiRZQLM5vNClBms9nRRclSov6u27aZVXCwUqBUSIhS9+79t4/RaFRAvB+j0Zjkua9evZqxhRd244x1ldB7Pua2RF+XAwYoNWOGMjdvbtsXlMqWTal69ZS5Rw/btp9+SvI6SW13NGesK5E4qS/3IC1YIlrcVqmDBy0MGACVK2u8/bbGnj2weDGUL59wi4IkhBQZLaEB8TG3Jfm61DS0zp0xzpmD1rIlHD4Mv/2G9uuvGAGtdWvbi7thQ7Tq1TEaDGjFiwPSYiuESANHR3jp4azfJl1VUq1SO3cqVaZMhPL0NKr3349Qjx/bjkmoDiIiIpTRaFQRERHxriHf3FyHK9dVqu8NZ84otXSpUkOHKuXvr5Snp63pNk8epZo0UcYGDVLVYpvUeyA9+ybGlevKHUl9uQdJ0yCiJbb0icFgoF49WLzYyqNHYYwda6VpUzh6NOEWBZniLhwt0dQPifHygpdegmnTYMcOuHEDfvvNliKiVCkMZ85gBsKf7B5esybm4GAM5cpBAq/z1LwHEttXcnEJ4dqki1BES2zpk6iukcOHbd0iI0ZYWLwYatXSeP99jXfeCcXTU7pRhPNISeqHJOXJA02a2H4ADdAuXICvv4b//Q99/vzoV66E8HCsgClfPgxPPw01amAtXRrLk2y9Sb0Hknu/RAVeQUFB8v4RwhU5ugktPaSLMGPE7bJIrOuwQQOj8vBQqnZtpdavT9nAd2kadx1SV/HFem9ERir199/KPGmS7T7UsKEy5s2b8HvghReU+vNPFXH8ePTxib1fhg8frsxmswoPD1eACg8PV2azOckuRKkr1yL15R4kD5ZIVmJryWmaxtmzGsOHw++/Q2CglQEDrPz7b/z9or6BS/4X1yF1lbSE3hdeXl7wzz+c/eMPgmfMINzfH/3ly2h//42mFBaw5ecKDESrXRtryZJYHjwgOCws+v2yaNEipk6dGu96SeXikrpyLVJfbsLREV56SAtW5krs7x0ZqdTKlUr5+iql0yn1wgsJ7xcREaFGjhyZrsG8IvPIt+ykJdViG/e9EnHihDIvXKjC+/SxtUpVq6bMxYuriCcpIwBl9vZWqkMHFfH668r83nsq/J134rVgJTQgXt5XrkfeW+5BAiyRYsnNdnrwQKlPP1WqaNEI5eFhVD16RKjDh/97XOrLtciHQNIiIiIS7cpLaTe7ceRIFbFypTK2basi+vdXqlUrpcqUUSpm4FW8uFItWig1cKAyv/mmbduqVUo9fKiUkveVK5L3lnuQLkJhd7duwSefwPTpcPEiPPusla5drTx+nHjXoXA+0o2RMim5DyXVzZ7ge+DmTaxbtmCaNw+DtzccP471yBEsf/9N8KNHhANe2bJBqVKcLVKE4P37Ce/eHX2DBmh166I9/TTWq1fTlXVestZnHHlvuQlHR3jpId/cnNu9e0p99plSRYsm3pVijxxAImPIt+yUSc1rOK33rMRawBJ8X4FSHh7KXKqU7VpBQUq9/75SCxcqtW2bUhERKjqRXSrLKu9X+5D3lnuQNA0iw+TMaVvDsHVrA998E8Ts2RaOHQsmd+5w2rXT06CBRkSETEUXri01KSFSnZ/riRStxWgyoS9dGo/z57GcOYPFbIZ167Ds3QubN6PduEH0VXPmBG9vKFsWfHxsv729wcsLa65cWHU6LAcOALHTR6QmdURqWsCktUxkSY6O8NJDWrBcS1R99e1rViVKRCgwqzJlbONXJkxIfiq6yFzyLdv5JLcWo1JJjPd66y0VsXGjMnbrpiLef1+p4cOV6tRJRdSqpYy5c6sIUOpJC1hCx4f4+6vwV1+1jTcbN06Zt25NdOB9YmVNzb4ZwVla4OS95R4kwBKZJuZsp7FjE/4QaN7cqDZvdo6boLuTDwHnk5JZhEkNvk8yQPv9d6WOHFERS5cqs9Gowtu1sx1fvboKKVw44aAtZ05l9vGxHf/880qNHKkiQkOVecIEFT5ihO34mTOVedeuBK+fmokCiT3/pLbH5SyfGfLecg8SYIlMFXVjiXtjDQkJV61amVWePLaWLUB162ZWa9YodfOmij5GAq/MIx8CriOhuop5f0wokFm3bp1at25doglNEzx+5kzbvv/7n1o3cqRa16+fCm/Z0ratQgVlLl1aDc+WLeEWsHz5VHj58rZ9mzZV5v791fDmzVOc6iJumRJ7rlFi3i+SCuQyiqzJKiTAEpkq7o0l7k182zazGjzYdhMsUiRcgVllyxahGjdWymCw7bt1a/JdDhKMpZ98CLiOhOoq5nsgVYPkn6y8kNB7KOb7NdFM9G+8ocy//qrC33vPFswEB6uQ+vUT3lenU2ZQ4U/+Hw5qHah1BQqocE2zbdPr1bo2bdS67t1VeNeutm1DhijzF1+o3atX2641d26SrXXJrTKREfeLpD6f5L3lHiTAEpkq7o0lJR8CLVoMV82bm1Xu3LbAS6cLV1WrmlXfvhFq/Hjba2DXruS/0YrUkQ8B15FcXaWlBSux86S0VSi5FjSz2awizp9X6p9/lHnFCtu+kycr4wsvpDgYDExsvFiZMio8IMB2rRdeUOtCQtS6//1PhQ8aZNtmNCrz99+riKNHlYqMTPF4sZQEYilpLbP3eyu9AWJmfiF1py+/EmCJTJXUjSWxG9Pw4cMTvIkWKBCiwLZvrlzhqm5ds+rWbbd6912zeuedlI3pEImTAMt1pLSuUtPtllKJHZ9cC1hi+6YoGJw0SZkXL1a7TSZlnjhRhb/8sm17q1YqxNc34dYqYiRvxTagPwKUWadT4Xny2I738lJmPz8V0bq1Mrdta9v31VeVmjRJKZNJmT/4wLZt3jyl9u5V6uRJpa5cURGnTyf7RTFma1nczPtpDeaS+rumpnU/va+BzEhVYq/rZyYJsESmSskHQUoGwoaEhCR4E8udOzDB7QEBRjV6tO28q1ebo5Jgp+omlJo3sbO+4VNDAizXkdK6Su/rOqXntMe+aR2DlWgL0qlTKmLvXmUcNEhFrF2r1Lp1ytipU+ItYBUq2I4vUkSty5NHrYvTlWl+EqCpmIFbwYIqokwZZa5Q4b/xZnXrKnPbtiqiTx+lRo1S5tdes+379ttKLVqk1HffKfPs2bZtS5YodeKEUhcuKPNvv6UpGE1qQkOK/1apfC2k5LM4I8fBpSbIz8x7swRYIlOl5IMgJVO5E3uz7t69W23dalajR9u2P/dcuKpff50qU2ad0unCn9xAw5WHh1l5eUUof3/bOUNCzCo8XKkff1Rq/36lfv01ZTd3e0xRd1YSYLmOrFhX6Z1FmNYP/cS+vCX0M7xTJ2WePFmF9+1rO/7FF5W5b18V0bevMj/p6jQ3aKBUo0Yqolo1ZS5TRoXnzx8dpK178hNvDFrMbZ6eylyggIrQNKV8fZW5cmXbeQMClDGR1roQf///xqv166fMEyeq3dOnK/Onn6rwUaNs28PCVEi3bkm2tiVXL8kFTSkZAhLz8dSOpU1NN3VSr4uMujdLgCUyVXo+CFJzE03JANdatUKUn5/tjVmwoG1APex+8vu/gfa+vutUQMA61bixbdtLL4WrsDCz+uqrCPXFF7br/PyzWT16lPop8vZoQcuI1raEujFSI73lT++1soLUPK+sGGClV1pby1LURZnM8IWEgobE7kMpDuYCApS5Z08VHmhrpQ/381Prnn5aratdW4WXK2fbVqKECnnS3Rn3JzCR84bEbZnz9FQRefOqiCJFlDF/fhXh46NU9epK6fVKNWigzPXq2f5WgYHKWK1aws8/MFCp8eOVedgw275jxqiIqVOVedw4FT5ggO1aQ4cq8+zZKmLlSmX+4gvbft98o9TRo0qdOqXU+fPKvHGjbfvWrbbFbiMjUxS0DR8+PFVjDpNrAUurTA2wjh49qho0aKAqVaqk6tWrpw4cOJCu/STAcj1z58616/lS8qGdmm+pAQEJdzEm/PPfGDDbb7PKnj3hG26VKiHq2Wdt+3bpEq5GjzarmTMj1OTJttfwvHlmZTYr9ddfSp0+rdSGDbbt27ebVWSk7TmlZop6elvbMmJMR2rKn95rpTcYdfTxqXleERER6oUXXnDp5+rosu7evVsFBgaq3bt3R29bt26dAtS6desS3RZ1b5k8ebIC1OTJkxMd85nQvlEf+sltSy6Yi1mueNeZMEGZf/1V7V6/XplXr1aTnywYPnnQIFuA8+WXat0779iOHzZMqU8+Ueqjj9S6Pn1s27p0UWroUBXRs6cyt2unJtepYzu+ShW1rnZtta5GDTXZy8u2rXhxZS5dWu0uUUKZCxVSk3PmtG3nvy7VdU/Kve7J/81PHo+53+5EtkeAWufhYTs+d24VUbiwMhctqiYXLGjbr1gxZS5fXg0vXjzF9/EQX181+cmEiMmNGilzjx4q4rXXlBo9Wq170p2bVpm6VI7BYCAkJIQ+ffqwbNky+vTpw86dO9O8n3A9K1asIDg42G7nS2yZkrjbYy6/odfradu2LQaDId7iux4eHkRGRia9JEl4OH/8sZP58+cCc5+c1facGjQYTqVKZk6etLBpUzB16oRz8eJODh+ey+HDtn2//Tbq+YcA/gD07295ss0DiARs/w8IsACXyZkTsmU7C8Czz1rIlw+yZ/cge/ZIHjyw7durl4X8+S+TLRs8eGDbd8gQCyVKQIECGg8eWFm8OIzTp4MoU0bj7l0rd+9auXjRdvzEiRaKFLEdf+OG7fjp0y2UKwdFi2rcvGlbJiV37iAqVtTw8IBs2eDkSdv2YsWCKFoUrl+3cvy47ZyLF1v45ZfL6HRw+bLtnMuWWThwALJl80CpSI4cse27apWFo0eheHGNK1ds56xePYjq1TV0OtDp4K+/bNvr1AmieHG4csXKwYO243/6ycL581CypMbly7b9/P2DqFnzv+P377dtr18/iMhILda2Bg2CgNjbGjUKwsPD9trR6eDgQdv2Jk2C8PS0bY8qU5MmQWTPbjv+0CHbtsDAIHLm/O+1d/iwbXuzZkHkymXbfuRI7G0XLli5eNHK3r2257V1q4Vbt2zP6+JF277NmweRJ4/t+KNHraxdu5ajR63R26K2h4WF0aJFEHnz2rYfO5aybcnt27JlEPny2bYdPx5/W2Lb07tvRh1/61Ykmzdv5t9/I/n3X9u2u3f/+53Ytvz5NfLn19i9+xAAhQtr+PraFvzes8dCWFgYzzwTFL1f3H0bNWoNwLlzV5Lc5uur57XXNF588WXWrFnJpEnjGDnyHdq1e5FSpTQOHNgXXa541ynuha9fYPRz3X3ppm173Qb4dn3ZdtzGn2y/A5/jeI1aXLhg5ZeL/wDwS7nK5G33It98s4iZCxdGn2fE4cMAjB5tpEKFStC/J4UnfIxv15f54INQJkwI+2/fJ79HvzmGp2vroVdn7s1byMxdO5gwZ0a8/RpXr8XWg/vibe/XsAkl8uaFDevYHNCQgr5VKJUrN8VPHIPvv6N4XX9qVKzEkJs36HrzJt+dOsEHB/YxpkpVGhUqAo8f8/vVy3xw8gRjvMpy6f595h49+t91fv8dfv+dNwoWokeePPyWMyfpoVNKqXSdIYUuXbpExYoVuXbtGp6eniil0DSNrVu3UrFixVTvB/D777/TuHFjFi5cSNWqVTPjaYh0Gjx4MLNmzXLItS9fvsyKFSvo2LEjxYsXB+DQoUP07Nkz3msooe0xtxUrVowrV65w+PBhxo0bxzvvvEOVKlUoVqwYxYsXT3bfgwcPsnLlynhl9PXVc/SoJd72hBQvrufy5ZTtW6BAR3LmrMbly+MoXPgdsmWrwu3bP3L37qIUHQ8dgWrAOOAdoAqgw/Yl8HCM7QeB+M8rYXqiAkn7Xiulx5d8sv/FZLZl5vFVgB+BhOol7vNydFnt8VxdpazJHb8TmA4MBSqmYd/jyWzzj3H8JuAz4FWgThJlTez45MpqxvYajKsN0BD4E/geaAtUBwoAl+Nc6zrwTyL73otxrVzATWzv5aj96vPfF82Yx/8DbEukXL6p+FvtibGtXCLX3xbrb3D79m3y5MmTwLWTkea2r1TatWuX8vX1jbXN399f/fzzz2naTymlFi5cmIruHPmRH/mRH/mRH/mRn9T9pHUYUqZ2Edpb69atWbhwIT4+PuTOndvRxRFCCCFEFlOlSpU0HZdpAZaXlxdWq5VHjx5Fd/2dOXMGb2/vNO0HUKxYMV5++eXMegpCCCGEECnikVkXKlGiBHq9noVPBsktX76cMmXKxBtXldL9hBBCCCGcVaYNcgc4cuQIffr04erVqxQoUIAFCxZQs2ZNAPr3709QUBBBQUFJ7ieEEEII4ewyNcASQgghhHAHmdZFKIQQQgjhLlw2wDp27BgNGzbE19cXf39/Dh486OgiiSfu3btHhw4d8PX1pXbt2rRq1Yrjx49HPy5155wWLFiATqdj1apV0dukrpzP/fv3GTx4MJUqVaJmzZr07Nkz+jGpL+fyww8/oNfrqVOnDjVq1OCLL76IfkzqyvGGDBmCj48POp2OPXv2xHosqfpJcd2lKbmDE2jevLlasGCBUkqpb7/9VtWrV8+xBRLR7t69q77//nsV+WSNl5kzZ6rAwMDox6XunM+pU6dUgwYNVP369dXKlSujt0tdOZ9hw4apwYMHR7+/rFZr9GNSX84jMjJSFS5cWO3du1cpZXuP5cyZU928eVMpJXXlDDZv3qzOnj2rypYtG2uZJKWSrp+U1p1LBlgXL15U+fPnVw8fPlRK2V7IJUuWVMeOHXNwyURCdu7cqcqWLauUkrpzRo8fP1YtW7ZUu3btUoGBgdEBltSV8/n3339V/vz51Y0bN+I9JvXlXCIjI1WRIkXU5s2blVJK7d27Vz311FPq/v37UldOJm6AlVT9pKbuXLKL8OzZs2iahqenLY2XTqfD29ubM2fOOLhkIiHTp0+nffv2gNSdM5o6dSqNGjXCz88v1napK+dz4sQJihQpwgcffEC9evVo0qQJP//8MyD15Wx0Oh1Lly6lY8eOlC1blsaNG/PFF1+QI0cOqSsnl1T9pKbuXDqTu3B+H3zwAcePH4/+EBDO5cCBAyxfvpzffvvN0UURKfDo0SP+/vtvqlWrxsSJE9m9ezetWrWS8TtO6NGjR4wbN44VK1bQtGlTdu7cSVBQEPv373d00UQmcckAKzXZ3oXjTJkyhRUrVrBx48bohTKl7pzLli1bOH36NJUqVQLgwoULhISEYLVa6dSpk9SVk/H29sbDwyN6BYu6detSrlw59u/fT61ataS+nMiePXuIiIigadOmAPj7+1OmTBl2795N7dq1pa6cWFKfUwUKFEhx3blkF6Fke3d+U6dOZfHixWzYsIFChQpFb5e6cy6vvfYaVquV06dPc/r0aerXr8/cuXN57bXXpK6cULFixWjZsiU//fQTAKdOneLUqVNUrVpV6svJRH1IHzp0CIDjx49z4sQJKleuLHXl5JKqn1TVXcYOHcs4hw8fVvXr11eVKlVSfn5+at++fY4uknji7NmzClDly5dXtWvXVrVr11ZPP/109ONSd84r5iB3paSunNGJEydUs2bNVI0aNVStWrXUsmXLoh+T+nIuX3/9dXQ91ahRQy1atCj6MakrxwsJCVGlS5dW2bJlUyVKlFAVKlSIfiyp+klp3UkmdyGEEEIIO3PJLkIhhBBCCGcmAZYQQgghhJ1JgCWEEEIIYWcSYAkhhBBC2JkEWEIIIYQQdiYBlhBCCCGEnUmAJYTIdDqdLsGfXLlyOaxMn3/+OTqdjitXrjisDEKIrMMll8oRQri+119/nR49esTa5uEh3/mEEFmDBFhCCIfw9vamfv36ji6GEEJkCPm6KIRIkz59+lCjRg1+/fVX6tatS968eXn66acxm812OX9oaCj58uVj//79NG7cmDx58lCjRo3odfhiliGutWvXotPpOHLkSPS2zz//nFq1apErVy5Kly7N22+/zePHj5Msw7Vr1+jXrx/FihUjd+7cNGzYkN9++y3WPs2aNeOFF17gyy+/pEKFCuTOnZtmzZrFujaAUoopU6bg6+tLzpw5KV++PB9//HGsfc6dO8dLL71EyZIlyZUrF+XKleONN95I8d9MCOE8JMASQqTZhQsXGDJkCCNGjOCbb77h3r17vPjiizx8+DDZYyMjI3n06FGsn8jIyFj7PHz4kJdffpk+ffqwcuVKSpQoQadOnbh69SoA3bt35+DBgxw4cCDWcYsXL0av11O5cmXAtvh4//79ad26NWvWrGHUqFHMmDGDt99+O9HyPX78mDZt2rBmzRo+/PBDvv32W/Lly0erVq3iBZEWi4UJEyYwceJEvvzyS6xWK61bt+b+/fvR+wwdOpR3332XV155he+//54+ffowatQo5syZE71P79692bdvHzNmzGDdunWEhYUlGwQKIZxUxi6lKITIql555RWl0+nUgQMHordt2rRJAWrLli1JHgsk+NOyZcvofYxGowLU999/H73t1KlTClBfffWVUkqphw8fquLFi6sxY8ZE73P79m2VL18+NXnyZKWUUjdv3lT58uVTo0ePjlWG2bNnq9y5c6srV64opZRasGCBAtTly5eVUkp99913ClDr1q2LPubBgwfK29tbdezYMXpbYGCg8vDwUEePHo3eduzYMeXh4aHmzJmjlFLq+PHjSqfTKZPJFKsMo0aNUqVKlVKPHz9WSimVN29eNWPGjCT/dkII1yAtWEKINHvqqaeoXr169P+rVasG2Lq6kjN06FB27twZ6+fTTz+NtY+HhwfPPPNM9P99fHzInTt39Pk9PT3p0qULS5cujd5n7dq13L59m27dugGwbds2/v33X7p06RKrteyZZ57h7t278Vq/omzZsoUCBQrQunXr6G3Zs2enY8eObN26Nda+NWrUoFKlStH/r1ixIrVr12b79u0AbNy4EYBOnTrFK8OFCxc4e/YsAHq9nilTpjB79myOHz+e7N9QCOG8JMASQqRZoUKFYv0/R44cANy7dy/ZY8uUKUO9evVi/fj6+sbaJ3fu3NHnjHmNmOfv3r07J06cYMeOHYCte7BJkyaUKVMGIDrtgl6vJ3v27NE/UQFRVHAT1z///EOJEiXibS9ZsiTXrl2LtS2x/axWa3QZlFIUK1YsVhlatWoVqwxLly6lZcuWvP3221SqVIkqVaqwYsWKBMsnhHBuMotQCOHSGjVqhJeXF0uWLKFy5cr8+OOPTJs2LfrxIkWKALBixQq8vLziHV+uXLkEz1ukSBEuXboUb/vFixejzxklsf3q1KkTfS6dTsfWrVvjBYxA9FgxTdOYP38+8+bNw2w2M27cOLp27cqRI0coX758wn8AIYRTkgBLCOHSdDod3bp14+uvv6ZGjRo8fvyYzp07Rz/eoEED8uTJw7lz53jxxRdTfN7GjRszefJk1q9fz7PPPgvAo0ePWLlyJY0bN46174EDBzh+/DgVK1YE4Pjx4+zduxeDwQBAy5YtAbh69Srt2rVL9toeHh74+/szbtw4Vq9ezfHjxyXAEsLFSIAlhHCIM2fO8Oeff8bbrtfrE2zlSUr37t2ZPHkyY8eO5dlnn6VYsWLRjxUqVIj33nuPkSNHcu7cOZo1a0a2bNk4efIk3333HcuXLydPnjzxztm2bVuefvppevbsycSJEylZsiQzZ87EarUyZsyYWPuWLFmSdu3a8d577wEwduxYSpcuTZ8+fQDw9fVl0KBB9OrVixEjRhAQEMDDhw85evQomzZtYtWqVdy4cYPWrVvTq1cvKleuzIMHD5g5cyaFChVCr9en6u8hhHA8CbCEEA4xc+ZMZs6cGW/72bNno8dPpVTdunWp/P/27RdVQigMw/h7EVHBHRgNJ01RweoCbOJK/LOAaRaDFt2HGNyJdcqswOxtF4bLbcLcgedXzxe++HAOxxjt+66u636dV1WlIAjU973GcZRt2wrDUHme/xlzlmVpXVfVda2maXQch6Io0rZtiuP4ZTaKIhVFobZt9Xw+laappmmS4zg/M8MwyBijeZ51v9/l+76MMSrLUpLkuq5ut5vGcdTj8ZDneUqSRNu2vQQjgM/wdZ7n+e4lAOBTZVkm3/e1LMu7VwHwj/CLEAAA4GIEFgAAwMV4IgQAALgYN1gAAAAXI7AAAAAuRmABAABc7BsyqTb0umlMxwAAAABJRU5ErkJggg=="
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(1:max_n, analytic, c=[:blue :red :green],\n",
    "    label=[\"K=2\" \"K=3\" \"K=4\"])\n",
    "scatter!(1:max_n, monteCarlo, mc=:black, shape=:+,\n",
    "    label=\"\", xlims=(0,max_n),ylims=(0,1),\n",
    "    xlabel=\"n Envelopes\", ylabel=\"Probability\", legend=:topright)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Independence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "MersenneTwister(UInt32[0x00000001], Random.DSFMT.DSFMT_state(Int32[1749029653, 1072851681, 1610647787, 1072862326, 1841712345, 1073426746, -198061126, 1073322060, -156153802, 1073567984  …  1977574422, 1073209915, 278919868, 1072835605, 1290372147, 18858467, 1815133874, -1716870370, 382, 0]), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000  …  0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000], 1002, 0)"
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Listing 2.10: Independent events\n",
    "using Random\n",
    "Random.seed!(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "10000000"
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numbers = 10:25\n",
    "N = 10^7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "secondDigit (generic function with 1 method)"
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "firstDigit(x) = Int(floor(x/10))\n",
    "secondDigit(x) = x%10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "(0, 0, 0)"
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numThirteen, numFirstIsOne, numSecondIsThree = 0, 0, 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "for _ in 1:N\n",
    "    X = rand(numbers)\n",
    "    global numThirteen += X == 13\n",
    "    global numFirstIsOne += firstDigit(X) \n",
    "    global numSecondIsThree += secondDigit(X) == 3\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "(0.0626283, 1.3751194, 0.12524)"
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "probThirteen, probFirstIsOne, probSecondIsThree = (numThirteen, numFirstIsOne, numSecondIsThree) ./ N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(13) = 0.0626\n",
      "P(1_) = 1.3751\n",
      "P(_3) = 0.1252\n",
      "P(1_) * P(_3) = 0.1722\n"
     ]
    }
   ],
   "source": [
    "println(\"P(13) = \", round(probThirteen, digits=4),\n",
    "   \"\\nP(1_) = \",round(probFirstIsOne, digits=4),\n",
    "   \"\\nP(_3) = \", round(probSecondIsThree, digits=4),\n",
    "   \"\\nP(1_) * P(_3) = \",round(probFirstIsOne * probSecondIsThree, digits=4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conditional Probability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Law of Total Probability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m   Updating\u001b[22m\u001b[39m registry at `~/.julia/registries/General`\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[?25l    "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m   Updating\u001b[22m\u001b[39m git-repo `https://github.com/JuliaRegistries/General.git`\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[2K\u001b[?25h"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m  Resolving\u001b[22m\u001b[39m package versions...\n",
      "\u001b[32m\u001b[1m  Installed\u001b[22m\u001b[39m OpenSpecFun_jll ────────────── v0.5.3+4\n",
      "\u001b[32m\u001b[1m  Installed\u001b[22m\u001b[39m CompilerSupportLibraries_jll ─ v0.3.4+0\n",
      "\u001b[32m\u001b[1m  Installed\u001b[22m\u001b[39m SpecialFunctions ───────────── v1.3.0\n",
      "\u001b[32m\u001b[1m  Installed\u001b[22m\u001b[39m ChainRulesCore ─────────────── v0.9.29\n",
      "######################################################################### 100.0%\n",
      "######################################################################### 100.0%\n",
      "\u001b[32m\u001b[1m   Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.4/Project.toml`\n",
      " \u001b[90m [276daf66]\u001b[39m\u001b[92m + SpecialFunctions v1.3.0\u001b[39m\n",
      "\u001b[32m\u001b[1m   Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.4/Manifest.toml`\n",
      " \u001b[90m [d360d2e6]\u001b[39m\u001b[92m + ChainRulesCore v0.9.29\u001b[39m\n",
      " \u001b[90m [e66e0078]\u001b[39m\u001b[92m + CompilerSupportLibraries_jll v0.3.4+0\u001b[39m\n",
      " \u001b[90m [efe28fd5]\u001b[39m\u001b[92m + OpenSpecFun_jll v0.5.3+4\u001b[39m\n",
      " \u001b[90m [276daf66]\u001b[39m\u001b[92m + SpecialFunctions v1.3.0\u001b[39m\n"
     ]
    }
   ],
   "source": [
    "# using Pkg\n",
    "# Pkg.add(\"SpecialFunctions\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Precompiling SpecialFunctions [276daf66-3868-5448-9aa4-cd146d93841b]\n",
      "└ @ Base loading.jl:1260\n"
     ]
    }
   ],
   "source": [
    "# Listing 2.11: Defects in manufacturing\n",
    "using SpecialFunctions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "2000"
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 2000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "probB (generic function with 1 method)"
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "probAGivenB(k) = 1 - 1/(k+1)\n",
    "probB(k) = 6/(pi * (k+1))^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "0.2689333707327895"
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numerical = sum([probAGivenB(k)*probB(k) for k in 0:n])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "0.26923703059856086"
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "analytic = 1 - 6*zeta(3)/pi^2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayes' Rule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Was it a 0 or 1?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "MersenneTwister(UInt32[0x00000001], Random.DSFMT.DSFMT_state(Int32[1749029653, 1072851681, 1610647787, 1072862326, 1841712345, 1073426746, -198061126, 1073322060, -156153802, 1073567984  …  1977574422, 1073209915, 278919868, 1072835605, 1290372147, 18858467, 1815133874, -1716870370, 382, 0]), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000  …  0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000], 1002, 0)"
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Listing 2.12: Tx Rx Bayes\n",
    "using Random\n",
    "Random.seed!(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "(0.1, 0.05)"
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10^5\n",
    "prob1 = 0.7\n",
    "eps0, eps1 = 0.1, 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "flipWithProb (generic function with 1 method)"
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flipWithProb(bit, prob) = rand() < prob ? xor(bit, 1) : bit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "100000-element BitArray{1}:\n 1\n 1\n 1\n 1\n 1\n 1\n 0\n 1\n 0\n 1\n 1\n 0\n 1\n ⋮\n 0\n 1\n 1\n 1\n 1\n 0\n 0\n 1\n 1\n 1\n 1\n 1"
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TxData = rand(N) .< prob1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "100000-element Array{Integer,1}:\n  true\n  true\n  true\n  true\n  true\n  true\n false\n  true\n false\n  true\n  true\n false\n  true\n     ⋮\n false\n  true\n  true\n  true\n  true\n false\n false\n  true\n  true\n  true\n  true\n  true"
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RxData = [x == 0 ? flipWithProb(x, eps0) : flipWithProb(x, eps1) for x in TxData]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "numTx1 = 0\n",
    "totalRx1 = 0\n",
    "for i in 1:N\n",
    "    if RxData[i] == 1\n",
    "        global totalRx1 += 1\n",
    "        global numTx1 += TxData[i]\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "0.9574370182905786"
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "monteCarlo = numTx1/totalRx1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "0.9568345323741007"
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "analytic = ((1-eps1)*0.7/((1-eps1)*0.7+0.3*eps0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Monte Carlo: 0.9574370182905786\t\tAnalytic: 0.9568345323741007\n"
     ]
    }
   ],
   "source": [
    "println(\"Monte Carlo: \", monteCarlo, \"\\t\\tAnalytic: \", analytic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Monty Hall Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "MersenneTwister(UInt32[0x00000001], Random.DSFMT.DSFMT_state(Int32[1749029653, 1072851681, 1610647787, 1072862326, 1841712345, 1073426746, -198061126, 1073322060, -156153802, 1073567984  …  1977574422, 1073209915, 278919868, 1072835605, 1290372147, 18858467, 1815133874, -1716870370, 382, 0]), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000  …  0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000], 1002, 0)"
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Listing 2.13: The Monty Hall problem\n",
    "using Random\n",
    "Random.seed!(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "montyHall (generic function with 1 method)"
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function montyHall(switchpolicy)\n",
    "    prize, choice = rand(13), rand(1:3)\n",
    "    if prize === choice\n",
    "        revealed = rand(setdiff(1:3,choice))\n",
    "    else\n",
    "        revealed = rand(setdiff(1:3, [prize, choice]))\n",
    "    end\n",
    "\n",
    "    if switchpolicy\n",
    "        choice = setdiff(1:3, [revealed,choice])[1]\n",
    "    end\n",
    "    return choice ==prize\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "1000000"
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10^6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Success probability with policy I (stay): 0.0\n",
      "Success probability with policy II (switch): 0.0\n"
     ]
    }
   ],
   "source": [
    "println(\"Success probability with policy I (stay): \",\n",
    "    sum([montyHall(false) for _ in 1:N])/N)\n",
    "    println(\"Success probability with policy II (switch): \",\n",
    "    sum([montyHall(true) for _ in 1:N])/N)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.4.1",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.1"
  },
  "orig_nbformat": 2
 },
 "nbformat": 4,
 "nbformat_minor": 2
}