inclusion probs

gistfile1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from itertools import chain, combinations, permutations\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "def powerset(iterable):\n",
    "    \"powerset([1,2,3]) --> [{} {1} {2} {3} {1,2} {1,3} {2,3} {1,2,3}]\"\n",
    "    s = list(iterable)\n",
    "    l = list(chain.from_iterable(combinations(s, r) for r in range(len(s)+1)))\n",
    "    return [frozenset(x) for x in l]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "256"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N=8\n",
    "\n",
    "p = 10* np.random.dirichlet(np.ones(N), 1)[0] # 1 distributions with N elements, unnormalized\n",
    "\n",
    "#enumerate the power set of N elements, with indexes, no order\n",
    "U = set(range(N))\n",
    "power_U = powerset(U)\n",
    "len(power_U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "264ba94aedb44153aead983e92a32066",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/256 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def Z(subset, p=p):\n",
    "    return np.sum([np.exp(p[i]) for i in subset])\n",
    "\n",
    "def p_permutation(seq, U=U, p=p):\n",
    "    prob = 1.0 if len(seq) > 0 else 0.0\n",
    "    remaining_choices = U.copy()\n",
    "    for e in seq:\n",
    "        prob *= np.exp(p[e]) / Z(remaining_choices)\n",
    "        remaining_choices -= {e}\n",
    "    return prob\n",
    "\n",
    "sampling_design = { e: 0.0 for e in power_U }\n",
    "for subset in tqdm(power_U):\n",
    "    subset_permutations = permutations(subset)\n",
    "    sampling_design[subset] = sum(p_permutation(perm) for perm in subset_permutations)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.000000000000007"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check that the sum of the sampling design is N\n",
    "\n",
    "sum(sampling_design.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999998"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now we decide for a sample size\n",
    "k = 3\n",
    "for subset, prob in sampling_design.items():\n",
    "    if len(subset) != k:\n",
    "        sampling_design[subset] = 0.0\n",
    "sum(sampling_design.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "frozenset({0, 1, 2}) 0.00915007531645386\n",
      "frozenset({0, 1, 3}) 0.0006042388924569105\n",
      "frozenset({0, 1, 4}) 0.0006911671590201153\n",
      "frozenset({0, 1, 5}) 0.0020004510397720254\n",
      "frozenset({0, 1, 6}) 0.002046551387552016\n",
      "frozenset({0, 1, 7}) 0.003562312702199033\n",
      "frozenset({0, 2, 3}) 0.006122983517974756\n",
      "frozenset({0, 2, 4}) 0.007007812568070215\n",
      "frozenset({0, 2, 5}) 0.020505587726829533\n",
      "frozenset({0, 2, 6}) 0.020987926938527125\n",
      "frozenset({0, 2, 7}) 0.037207963043798684\n",
      "frozenset({0, 3, 4}) 0.00046344608495887897\n",
      "frozenset({0, 3, 5}) 0.0013412044996397867\n",
      "frozenset({0, 3, 6}) 0.0013721044521920545\n",
      "frozenset({0, 3, 7}) 0.002387730964937475\n",
      "frozenset({0, 4, 5}) 0.001534197785464144\n",
      "frozenset({0, 4, 6}) 0.0015695466350435388\n",
      "frozenset({0, 4, 7}) 0.0027315133260899837\n",
      "frozenset({0, 5, 6}) 0.004545867621371258\n",
      "frozenset({0, 5, 7}) 0.007922693688366761\n",
      "frozenset({0, 6, 7}) 0.008106043270404786\n",
      "frozenset({1, 2, 3}) 0.010379221379189615\n",
      "frozenset({1, 2, 4}) 0.011878081227389424\n",
      "frozenset({1, 2, 5}) 0.034713009720573576\n",
      "frozenset({1, 2, 6}) 0.0355280473321789\n",
      "frozenset({1, 2, 7}) 0.06290262578201915\n",
      "frozenset({1, 3, 4}) 0.0007837179077472724\n",
      "frozenset({1, 3, 5}) 0.00226834651082399\n",
      "frozenset({1, 3, 6}) 0.002320622452453476\n",
      "frozenset({1, 3, 7}) 0.004039543483677852\n",
      "frozenset({1, 4, 5}) 0.002594719919284697\n",
      "frozenset({1, 4, 6}) 0.0026545203445850767\n",
      "frozenset({1, 4, 7}) 0.004621016805466166\n",
      "frozenset({1, 5, 6}) 0.007687375341343069\n",
      "frozenset({1, 5, 7}) 0.013398244865166543\n",
      "frozenset({1, 6, 7}) 0.013708167097873892\n",
      "frozenset({2, 3, 4}) 0.00794950908560696\n",
      "frozenset({2, 3, 5}) 0.023255386389736604\n",
      "frozenset({2, 3, 6}) 0.02380221118476094\n",
      "frozenset({2, 3, 7}) 0.04218635734886808\n",
      "frozenset({2, 4, 5}) 0.02660710081524312\n",
      "frozenset({2, 4, 6}) 0.02723246584100966\n",
      "frozenset({2, 4, 7}) 0.04825103172577559\n",
      "frozenset({2, 5, 6}) 0.07929283520988797\n",
      "frozenset({2, 5, 7}) 0.1398750132485948\n",
      "frozenset({2, 6, 7}) 0.14312073207320758\n",
      "frozenset({3, 4, 5}) 0.0017396636457933927\n",
      "frozenset({3, 4, 6}) 0.0017797484714791268\n",
      "frozenset({3, 4, 7}) 0.0030974870884876438\n",
      "frozenset({3, 5, 6}) 0.005154540312203458\n",
      "frozenset({3, 5, 7}) 0.008983525004797457\n",
      "frozenset({3, 6, 7}) 0.009191404822916726\n",
      "frozenset({4, 5, 6}) 0.005896057917664715\n",
      "frozenset({4, 5, 7}) 0.010275920793631833\n",
      "frozenset({4, 6, 7}) 0.010513679944056052\n",
      "frozenset({5, 6, 7}) 0.030456650285382546\n"
     ]
    }
   ],
   "source": [
    "#print non-zero possible samples\n",
    "for subset, prob in sampling_design.items():\n",
    "    if prob > 0:\n",
    "        print(subset, prob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Compute the inclusion probabilities of each element\n",
    "inclusion_probabilities = { e: 0.0 for e in U }\n",
    "for subset, prob in sampling_design.items():\n",
    "    for e in subset:\n",
    "        inclusion_probabilities[e] += prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: 0.14186141862112295,\n",
       " 1: 0.2275320566672267,\n",
       " 2: 0.8179559774756961,\n",
       " 3: 0.15922299350070246,\n",
       " 4: 0.17987240509186758,\n",
       " 5: 0.4300483923415712,\n",
       " 6: 0.436967098936094,\n",
       " 7: 0.6065396573657187}"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inclusion_probabilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#check that the sum of the inclusion probabilities is k\n",
    "\n",
    "sum(inclusion_probabilities.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.03288434, 0.08287621, 0.26631456, 0.0449561 , 0.05779543,\n",
       "       0.15518749, 0.15716   , 0.20282588])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#print normalized p\n",
    "p/sum(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "frozenset({1, 2, 6})\n",
      "frozenset({1, 3, 7})\n",
      "frozenset({2, 5, 7})\n",
      "frozenset({2, 5, 6})\n",
      "frozenset({1, 5, 7})\n",
      "frozenset({2, 5, 6})\n",
      "frozenset({4, 5, 7})\n",
      "frozenset({1, 2, 7})\n",
      "frozenset({2, 5, 7})\n",
      "frozenset({1, 2, 7})\n"
     ]
    }
   ],
   "source": [
    "#sample 10 times from the sampling design\n",
    "for i in range(10):\n",
    "    sample = np.random.choice(list(sampling_design.keys()), p=list(sampling_design.values()))\n",
    "    print(sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fdf63f612b0>"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACW+ElEQVR4nOzdeXxcdb34/9eZPZNM9j1Nm+77vlGWQqG1gKIi/uSqV4SLKAgCoiJchCqLIAoiyvIVRdCL4nW5ioAgFAqylO57m+5N2+zbLJl9zvn9MUmaPbMlM0nez8cjj7Yz55z5JM3MvOfzeb/fH0XTNA0hhBBCiCTRJXsAQgghhBjbJBgRQgghRFJJMCKEEEKIpJJgRAghhBBJJcGIEEIIIZJKghEhhBBCJJUEI0IIIYRIKglGhBBCCJFUhmQPIBKqqlJdXY3NZkNRlGQPRwghhBAR0DQNp9NJaWkpOl3/8x8jIhiprq6mvLw82cMQQgghRAxOnjzJuHHj+r1/RAQjNpsN2r+ZzMzMZA9HCCGEEBFwOByUl5d3vo/3Z0QEIx1LM5mZmRKMCCGEECPMYCkWksAqhBBCiKSSYEQIIYQQSSXBiBBCCCGSakTkjAghRKQ0TSMYDBIKhZI9FCFGPb1ej8FgiLvthgQjQohRw+/3U1NTg9vtTvZQhBgzrFYrJSUlmEymmK8hwYgQYlRQVZVjx46h1+spLS3FZDJJk0QhhpCmafj9fhoaGjh27BhTp04dsLHZQCQYEUKMCn6/H1VVKS8vx2q1Jns4QowJaWlpGI1GTpw4gd/vx2KxxHQdSWAVQowqsX4yE0LEJhHPOXnWCiGEECKpJBgRQggxqIqKCh577LGEXOv48eMoisKOHTsScr1kuOCCC7j11lvjusaGDRtQFIXW1tZ+j3nuuefIzs7u/Pf3v/99FixY0Pnvq6++mk9/+tNxP06ySTAihBBJdPXVV6MoCoqiYDKZmDJlCvfeey/BYLDzmGeeeYb58+eTkZFBdnY2Cxcu5MEHH+z3mqn+Zl9eXk5NTQ1z5sxJ9lBS3pVXXsnBgwf7vf9nP/sZzz33XOe/+wqSzj77bGpqasjKyhrSscZDEliFECLJLr74Yn7zm9/g8/l49dVXufHGGzEajdx55508++yz3HrrrTz++OOcf/75+Hw+du3axZ49e5I97Jjp9XqKi4uTPYw++f3+uEpUEy0tLY20tLR+748kwDCZTCn78+4gMyNCCJFkZrOZ4uJiJkyYwA033MDq1at56aWXAHjppZf43Oc+x7XXXsuUKVOYPXs2n//853nggQcivn7HNP369etZsmQJVquVs88+m8rKym7H/eMf/2Dp0qVYLBby8/O5/PLL+7xeXzMvra2tKIrChg0bAGhpaeGLX/wiBQUFpKWlMXXqVH7zm9/0e/4777zDsmXLMJvNlJSUcMcdd3SbHbrgggu4+eabuf3228nNzaW4uJjvf//7A37fHUsYP/jBDygoKCAzM5Prr78ev9/f7bo33XQTt956K/n5+axduzai8QAEg0FuuukmsrKyyM/P5+6770bTtM77f/e737FkyRJsNhvFxcV84QtfoL6+vtc433//febNm4fFYuGss87qFmj2XKbp73vs+Ps777zDz372s87ZtuPHj/e5TPPee+9x3nnnkZaWRnl5OTfffDNtbW2d9z/55JNMnToVi8VCUVERn/3sZwf8WcdLghEx5gVCarKHIIaApmm4/cGkfHV9Q4pFWlpa5xtmcXExGzdu5MSJE3H/TO666y4eeeQRtmzZgsFg4L/+678673vllVe4/PLLufTSS9m+fTvr169n2bJlMT/W3Xffzb59+/jnP//J/v37eeqpp8jPz+/z2NOnT3PppZeydOlSdu7cyVNPPcWvf/1r7r///m7HPf/886Snp/PRRx/x8MMPc++99/LGG28MOI7169ezf/9+NmzYwB/+8Af++te/8oMf/KDXdU0mE++//z5PP/10VOMxGAxs2rSJn/3sZzz66KP86le/6rw/EAhw3333sXPnTv72t79x/Phxrr766l5j/M53vsMjjzzC5s2bKSgo4LLLLiMQCET0c+7qZz/7GStWrOC6666jpqaGmpoaysvLex135MgRLr74Yq644gp27drFH//4R9577z1uuukmALZs2cLNN9/MvffeS2VlJa+99horV66MejzRkGUaMaY5vAG8/hCFmbHVxovU5QmEmHXP60l57H33rsVqiv7lVdM01q9fz+uvv843vvENANatW8dnPvMZKioqmDZtGitWrODSSy/ls5/9bNQllQ888ADnn38+AHfccQcf//jH8Xq9WCwWHnjgAf7jP/6j2xv1/Pnzo/4eOlRVVbFw4UKWLFkC7Qmw/XnyyScpLy/nF7/4BYqiMGPGDKqrq/nud7/LPffc0/l9zps3j3Xr1gEwdepUfvGLX7B+/XrWrFnT77VNJhPPPvssVquV2bNnc++99/Kd73yH++67r/O6U6dO5eGHH+4856677opoPOXl5fz0pz9FURSmT5/O7t27+elPf8p1110H0C3YmzRpEo8//jhLly7F5XKRkZHRed+6des6v4fnn3+ecePG8X//93987nOfi+pnnpWVhclkwmq1Drgs8+CDD/LFL36xM7dk6tSpncuATz31FFVVVaSnp/OJT3wCm83GhAkTWLhwYVRjiZbMjIgxrcnlx+6J/hOIEIn08ssvk5GRgcVi4ZJLLuHKK6/sXIIoKSnhww8/ZPfu3dxyyy0Eg0G+/OUvc/HFF6Oq0c3qzZs3r/PvJSUlAJ3LBjt27OCiiy5K2Pd0ww038OKLL7JgwQJuv/12Pvjgg36P3b9/PytWrOjWMfecc87B5XJx6tSpPsff8T30tezR1fz587s1wVuxYgUul4uTJ0923rZ48eKYxnPWWWd1O2bFihUcOnSoc1+krVu3ctlllzF+/HhsNltnIFhVVdXt8VasWNH599zcXKZPn87+/fsH/L7isXPnTp577jkyMjI6v9auXdvZxXjNmjVMmDCBSZMm8aUvfYkXXnhhyLdYkJkRMaY1uXzJHoIYImlGPfvuXZu0x47GqlWreOqppzCZTJSWlmIw9H5pnjNnDnPmzOHrX/86119/Peeddx7vvPMOq1ativhxjEZj59873kQ7ApqBkiR76pgZ6Loc1XNZ4ZJLLuHEiRO8+uqrvPHGG1x00UXceOON/OQnP4n4cQYaf8f3EG1A1pf09PS4r9FTW1sba9euZe3atbzwwgsUFBRQVVXF2rVru+WsJIPL5eJrX/saN998c6/7xo8fj8lkYtu2bWzYsIF//etf3HPPPXz/+99n8+bNA+avxEOCETFmBUIqdk8ARQFV1dDpZB+T0URRlJiWSpIhPT2dKVOmRHz8rFmzoP0NL1HmzZvH+vXrueaaawY9tqCgAICamprO6fu+yogLCgr48pe/zJe//GXOO+88vvOd7/QZjMycOZO//OUvaJrWGSS9//772Gw2xo0bF9f3tXPnTjweT2ewtXHjRjIyMvrMpYh2PB999FG38zZu3MjUqVPR6/UcOHCApqYmHnrooc7H2rJlS5+Pt3HjRsaPHw/tib8HDx5k5syZMX2/JpNp0B2rFy1axL59+wb8nTMYDKxevZrVq1ezbt06srOzeeutt/jMZz4T07gGI8s0YsxqbvOjaaCq4PQGIzhDiOF3ww03cN999/H+++9z4sQJNm7cyFVXXUVBQUG36f14rVu3jj/84Q+sW7eO/fv3s3v3bn70ox/1eWxaWhpnnXUWDz30EPv37+edd97he9/7Xrdj7rnnHv7+979z+PBh9u7dy8svv9zvG+zXv/51Tp48yTe+8Q0OHDjA3//+d9atW8dtt90Wd6txv9/Ptddey759+3j11VdZt24dN91004DXjXQ8VVVV3HbbbVRWVvKHP/yBn//859xyyy3QZYbh5z//OUePHuWll17ivvvu6/Px7r33XtavX8+ePXu4+uqryc/PH7SRWX8qKir46KOPOH78OI2NjX3OHH33u9/lgw8+4KabbmLHjh0cOnSIv//9750JrC+//DKPP/44O3bs4MSJE/z2t79FVVWmT58e05giIcGIGLMauyzRSN6ISFWrV69m48aN/H//3//HtGnTuOKKK7BYLKxfv568vLyEPc4FF1zAn/70J1566SUWLFjAhRdeyKZNm/o9/tlnnyUYDLJ48WJuvfXWXpUmJpOJO++8k3nz5rFy5Ur0ej0vvvhin9cqKyvj1VdfZdOmTcyfP5/rr7+ea6+9tleAE4uLLrqIqVOnsnLlSq688ko++clPDloSHOl4rrrqKjweD8uWLePGG2/klltu4atf/Sq0zwo999xz/OlPf2LWrFk89NBD/S5RPfTQQ9xyyy0sXryY2tpa/vGPf8Tc6+Tb3/42er2eWbNmdS4N9TRv3jzeeecdDh48yHnnncfChQu55557KC0tBSA7O5u//vWvXHjhhcycOZOnn36aP/zhD8yePTumMUVC0eKtQRsGDoeDrKws7HY7mZmZyR6OGCXePdiAPxj+1FCUaWHuuNTtTigG5/V6OXbsGBMnTox551Axulx99dW0trbyt7/9LdlDGdUGeu5F+v4tMyNiTHJ4A52BCECrJ7kJZUIIMZZJMCLGpCZX9+DDF1DxBgZO+hJCCDE0og5G3n33XS677DJKS0tRFCWi6a8NGzawaNEizGYzU6ZM6bapjxDJ0FdJr+SNCDG6PPfcc7JEM0JEHYy0tbUxf/58nnjiiYiOP3bsGB//+MdZtWoVO3bs4NZbb+UrX/kKr7+enM6IQnSU9PYkwYgQQiRH1EX4l1xyCZdccknExz/99NNMnDiRRx55BNrrt9977z1++tOfdm5IJMRwanKFS3p7kmBECCGSY8hzRj788ENWr17d7ba1a9fy4Ycf9nuOz+fD4XB0+xIiURr76brq9AZQ1ZQvLhNCiFFnyIOR2tpaioqKut1WVFSEw+HA4/H0ec6DDz5IVlZW59dAnfKEiFZzW9+VM6oarrIRQggxvFKymubOO+/Ebrd3fnXd0EiIePQs6e1JlmqEEGL4DfnGDcXFxdTV1XW7ra6ujszMzH43ZjKbzZjN5qEemhiDGp0Db4wnwYgQQgy/IZ8ZWbFiBevXr+922xtvvJHQPRWEiFRTP0s0HSQYEaK7ri0cjh8/jqIofW6KN9Q2bNiAoii0trYm5HoXXHABt956a0KuJeIXdTDicrnYsWNH5y/jsWPH2LFjR2f/+zvvvJOrrrqq8/jrr7+eo0ePcvvtt3PgwAGefPJJ/vd//5dvfvObifw+hBhUIKTiGCTYkOZnYrg1NDRwww03MH78eMxmM8XFxaxdu5b3338/2UPrpby8nJqaGubMmZPsocTtr3/9a78b14nhF/UyzZYtW1i1alXnv2+77TYAvvzlL/Pcc89RU1PTbWOeiRMn8sorr/DNb36Tn/3sZ4wbN45f/epXUtYrhl1/Jb09tboDFGfph2NIQnDFFVfg9/t5/vnnmTRpEnV1daxfv56mpqZkD60XvV5PcXFxsoeRELm5uckegugi6pmRCy64AE3Ten11dFV97rnn2LBhQ69ztm/fjs/n48iRI1x99dWJ+w6EiFB/Jb09yVKNGC6tra38+9//5kc/+hGrVq1iwoQJLFu2jDvvvJNPfvKTncc9+uijzJ07l/T0dMrLy/n617+Oy+XqvP+5554jOzubl19+menTp2O1WvnsZz+L2+3m+eefp6KigpycHG6++WZCoTMzfxUVFdx33318/vOfJz09nbKysgEbWvZcpulYOlm/fj1LlizBarVy9tlnU1lZ2e28+++/n8LCQmw2G1/5yle44447WLBgwYA/m1dffZVp06aRlpbGqlWrOH78eK9j3nvvPc477zzS0tIoLy/n5ptvpq2trfP+J598kqlTp2KxWCgqKuKzn/1s5309l2lqamr4+Mc/TlpaGhMnTuT3v/89FRUVPPbYY53HKIrCr371Ky6//HKsVitTp07lpZdeGvD7EJFJyWoaIYZCfyW9PUkwMkpoGvjbkvMV4WboGRkZZGRk8Le//Q2fr/9gWafT8fjjj7N3716ef/553nrrLW6//fZux7jdbh5//HFefPFFXnvtNTZs2MDll1/Oq6++yquvvsrvfvc7/t//+3/8+c9/7nbej3/8Y+bPn8/27du54447uOWWW3jjjTei+lHfddddPPLII2zZsgWDwcB//dd/dd73wgsv8MADD/CjH/2IrVu3Mn78eJ566qkBr3fy5Ek+85nPcNlll7Fjx47OAKarI0eOcPHFF3PFFVewa9cu/vjHP/Lee+9x0003Qfss/s0338y9995LZWUlr732GitXruz3Ma+66iqqq6vZsGEDf/nLX/jlL39JfX19r+N+8IMf8LnPfY5du3Zx6aWX8sUvfpHm5uaofl6ityGvphEiFdg9A5f0duXyBQipGnqdMuTjEkMo4IYflibnsf+7Gkzpgx5mMBh47rnnuO6663j66adZtGgR559/Pv/xH//BvHnzOo/r+gm+oqKC+++/n+uvv54nn3yy8/ZAIMBTTz3F5MmTAfjsZz/L7373O+rq6sjIyGDWrFmsWrWKt99+myuvvLLzvHPOOafzjX7atGm8//77/PSnP2XNmjURf7sPPPAA559/PgB33HEHH//4x/F6vVgsFn7+859z7bXXcs011wBwzz338K9//avbzE5PHd9HR+fu6dOns3v3bn70ox91HvPggw/yxS9+sfNnM3XqVB5//HHOP/98nnrqKaqqqkhPT+cTn/gENpuNCRMmsHDhwj4f78CBA7z55pts3ryZJUuWAPCrX/2KqVOn9jr26quv5vOf/zwAP/zhD3n88cfZtGkTF198ccQ/L9GbzIyIMaGvjfH6o6rhbqxCDIcrrriC6upqXnrpJS6++OLOjUW7bij65ptvctFFF1FWVobNZuNLX/oSTU1NuN3uzmOsVmtnIEJ7c8mKigoyMjK63dbz037PysYVK1awf//+qL6HroFTSUkJQOfjVFZWsmzZsm7H9/x3T/v372f58uUDjnPnzp0899xznbNLGRkZrF27FlVVOXbsGGvWrGHChAlMmjSJL33pS7zwwgvdfl5dVVZWYjAYWLRoUedtU6ZMIScnZ8DvNT09nczMzD5nUER0ZGZEjAmDlfT21OoOkG01Ddl4xDAwWsMzFMl67ChYLBbWrFnDmjVruPvuu/nKV77CunXruPrqqzl+/Dif+MQnuOGGG3jggQfIzc3lvffe49prr8Xv92O1hh/LaDR2u6aiKH3epqqRzRBGo+vjKEp4RnEoHqcrl8vF1772NW6++eZe940fPx6TycS2bdvYsGED//rXv7jnnnv4/ve/z+bNm8nOzo75cYfrZzrWSDAiRj1/cPCS3p4kb2QUUJSIlkpS0axZszp7e2zduhVVVXnkkUfQ6cKT2f/7v/+bsMfauHFjr3/PnDkzYdefPn06mzdv7tbyYfPmzQOeM3PmzF6JoT3HuWjRIvbt28eUKVP6vY7BYGD16tWsXr2adevWkZ2dzVtvvcVnPvOZXmMMBoNs376dxYsXA3D48GFaWlqi+l5F7GSZRox6zW2RlfR2JcGIGA5NTU1ceOGF/M///A+7du3i2LFj/OlPf+Lhhx/mU5/6FLQvFwQCAX7+859z9OhRfve73/H0008nbAzvv/8+Dz/8MAcPHuSJJ57gT3/6E7fcckvCrv+Nb3yDX//61zz//PMcOnSI+++/n127dnXOoPTl+uuv59ChQ3znO9+hsrKS3//+992WrQC++93v8sEHH3DTTTexY8cODh06xN///vfOBNaXX36Zxx9/nB07dnDixAl++9vfoqoq06dP7/V4M2bMYPXq1Xz1q19l06ZNbN++na9+9aukpaUNOE6ROBKMiFEv0pLervxBFY9fmp+JoZWRkcHy5cv56U9/ysqVK5kzZw5333031113Hb/4xS8AmD9/Po8++ig/+tGPmDNnDi+88AIPPvhgwsbwrW99iy1btrBw4ULuv/9+Hn300YT2gfriF7/InXfeybe//W0WLVrEsWPHuPrqq7FYLP2eM378eP7yl7/wt7/9jfnz5/P000/zwx/+sNsx8+bN45133uHgwYOcd955LFy4kHvuuYfS0nDScnZ2Nn/961+58MILmTlzJk8//TR/+MMfmD17dp+P+dvf/paioiJWrlzJ5ZdfznXXXYfNZhtwnCJxFE2L9jPj8HM4HGRlZWG328nMzEz2cMQI8+7BhograbqaU5ZFcZa8EI0UXq+XY8eOMXHiRHkDiVBFRQW33nrrsLdFX7NmDcXFxfzud78b1seNxqlTpygvL+9MHhb9G+i5F+n7t+SMiFEtmpLenlo9fglGhIiT2+3m6aefZu3atej1ev7whz/w5ptvRt3LZKi99dZbuFwu5s6dS01NDbfffjsVFRUD9iYRiSPBiBjVoinp7cnulrwRIeKlKAqvvvoqDzzwAF6vl+nTp/OXv/yF1atXJ3to3QQCAf77v/+bo0ePYrPZOPvss3nhhRd6Vc+IoSHBiBjVoi3p7crlC0rzMzGq9dViPdHS0tJ48803h/xx4rV27VrZMy2JJIFVjFqxlPR2pWnEdb4QQojISDAiRq1YSnp7apVgRAghhpwEI2LUiqWktyfpNzLyjIACQSFGlUQ85yQYEaNWPPkiHSQYGTk6Eg37239ECDE0Op5z8ST7SgKrGJXsngCBGEt6uwoEVdz+IFaTPFVSnV6vJzs7u3PTMqvVKt0zhRhCmqbhdrupr68nOzsbvV4f87XkFVaMSvGU9PZk9wQkGBkhiouLocuOsUKIoZednd353IuVvMKKUSkRSzQdWt0BSrLSEnY9MXQURaGkpITCwkICAVliE2KoGY3GuGZEOkgwIkadeEt6e5K8kZFHr9cn5AVSCDE8JIFVjDqJKOntqs0XJBiKP/9ECCFE3yQYEaNOIkp6u9I0cHiDCb2mEEKIMyQYEaOKpmkJzRfpIEs1QojR6mSzG5cvuR+4JBgRo4rDG0xISW9Pre7EBzhCCJFsRxtcVNY6k94sUBJYxaiSyJLermRmRAgx2hysc1LVlBpNAmVmRIwqQ7FEAxAMabQleRpTCCESQdM09lbbUyYQQYIRMZokuqS3J5kdEUKMdKqqseuUnZpWb7KH0o0EI2LUSHRJb0+tbglGhBAjVzCksv1kKw3OoVnOjofkjIhRI9ElvT3JzIgQYqQKhFR2nGzFnqIfqiQYEaPCUJX0duX2h5ufGfQyoSiEGDm8gRDbq1pTOu9NXlXFqDBUJb1daZrMjgghRhaPP8TWEy0pHYggwYgYLYZ6iaaDBCNCiJHC5Quy5UQzHn8o2UMZlCzTiFGhyTU8TclaJRgRQowAdk+A7VUtBEPJbWYWKQlGxIjnD6o4vcMTJAxl6bAQQiRCc5ufnadaCY2QQARZphGjwVCX9HYVDGlJ38NBCCH6U+/0suNky4gKRJBgRIwGw5Uv0kHyRoQQqai61cPuU3bUoc3lHxISjIgRbThKentK1Tp9IcTYdbLZzb5qx7DNEiea5IyIEc3hGfqS3p5aPbKDrxAidRxpcHGsoS3Zw4iLzIyIEa2xbfjbGrt9IQKhETgPKoQYdSprnSM+EEGCETHSDVdJb0+SNyKESCZN09hz2s7J5tTZeTceEoyIEWuod+kdiAQjQohk6dh5t9aeWjvvxkNyRsSI1ZSEJZoOEowIIZIhGFLZecpOyzAn7g81CUbEiJWsJRragxFN01AUJWljEEKMLYGQyvaq1lHZfFGWacSIlIyS3q5C0vxMCDGMvIEQW463jMpABAlGxEiVjJLenmSpRggxHEbKzrvxkGBEjEjJKOntSYIRIcRQG0k778ZDckbEiJTMfJEO0olVCDGU7O4A20+OnJ134yEzI2LESWZJb1dufwh/kpeKhBCjU3Obn21jJBBBghExEiWypHdbVQtvHaiP+XxZqhFCJNpI3Xk3HrJMI0acRC3RuP1BfvnuUYKqxqT8dCry06O+ht0ToMBmTsh4hBCiutXD/pqRu+FdrGRmRIwoiSzp3X6ylaAafsbvrXHEdA2ZGRFCJEpV08jeeTceEoyIESWRJb2bjzd3/n1fdWzBiKO9+ZkQQsTjSIOLg3XOZA8jaSQYESNKokp6Xb4g+6vPPPEPN7jwBaIvnQup0vxMCBGf0bLzbjwkGBEjSqLyRXZUtRLSNMblpJGXbiKkahysd8V0rVYp8RVCxGC07bwbDwlGxIjhC4YSVtLbsUSztCKXWSWZAOyttsd0LckbEUJEazTuvBsPCUbEiNGcoMRVpzfA/tpwjsiSCTnMKg0HI/tiTGJNhZ4nQoiRIxhS2X6ylQZn8jtJpwoJRsSIkaglmu1VragajM+1MpHTnG06hgJUt3ppdUf/GNL8TAgRKX9QZVtVKy1J3OgzFUkwIkYETdNodCXmU0THEs3yCTYWb/gS573/Jc7KCc+KxDo70uqRFxYhxMC8gfCGdzKb2psEI2JEcHiCCWmL7PAEONBePndJxhHM3kZ0WpBPp+8FYG8cJb5CCNEftz846nfejYcEI2JESFRJ77aqFjQNKvKsTG15p/P25ep2gPbOh9EHPZLEKoToj9MbYMvxllG/8248YgpGnnjiCSoqKrBYLCxfvpxNmzYNePxjjz3G9OnTSUtLo7y8nG9+85t4vZJBLCKXqHyRLSdaAFg6IYeC02923j7OvpV0QwiHN8ipVk/U13V4gtL8TAjRi90dYOuJFskrG0TUwcgf//hHbrvtNtatW8e2bduYP38+a9eupb6+783Gfv/733PHHXewbt069u/fz69//Wv++Mc/8t///d+JGL8YAxJV0mv3BKhsX6JZk12DxVNH0JCOz5KPIeTh07knIMZurCFVwynTr0KILppcPrZVjZ2dd+MRdTDy6KOPct1113HNNdcwa9Ysnn76aaxWK88++2yfx3/wwQecc845fOELX6CiooKPfexjfP7znx90NkWIDokq6d16IrxEMyk/vXOJpql4JU3F5wHwMdMeiCOJ1S7Nz4QQ7eodXnaeaiWkSiASiaiCEb/fz9atW1m9evWZC+h0rF69mg8//LDPc84++2y2bt3aGXwcPXqUV199lUsvvbTfx/H5fDgcjm5fYuxK3BJNuIpmScWZJZqGsos6g5F53i0AHKxzEghFP6UqeSNCCNp33t192o4qKzMRM0RzcGNjI6FQiKKiom63FxUVceDAgT7P+cIXvkBjYyPnnnsumqYRDAa5/vrrB1ymefDBB/nBD34QzdDEKJWokt5Wt59DdeF27xfkO8nYcwhVMdBYcj6KpqKhkOM6zFSLg0PeTA7Xu5jZ3pk1UhKMCCGqmtxjesO7WA15Nc2GDRv44Q9/yJNPPsm2bdv461//yiuvvMJ9993X7zl33nkndru98+vkyZNDPUyRouyeQELWW7ecaEEDJhekM6X5XQBaCpcRNGURMOfgyJ0LwGezD0KMJb4efwhfULLlhRirxvrOu/GIamYkPz8fvV5PXV1dt9vr6uooLi7u85y7776bL33pS3zlK18BYO7cubS1tfHVr36Vu+66C52udzxkNpsxm83RfSdiVGpM1BLN8fYqmopcCqvfAKCh9MxyY1PxeWQ17+JcZSewJPa8EU+AQps+IWMWQowclbVO2fAuDlHNjJhMJhYvXsz69es7b1NVlfXr17NixYo+z3G73b0CDr0+/GItpZBiME0JWKJpbvNzuMGFApxdrJLVGO4p0lB20ZnHKV4JwFTXZvSEONnsxumNftlFkliFGFtk593EiHqZ5rbbbuOZZ57h+eefZ//+/dxwww20tbVxzTXXAHDVVVdx5513dh5/2WWX8dRTT/Hiiy9y7Ngx3njjDe6++24uu+yyzqBEiL74giGc3vjLZTsSV6cUZjCl5T0UNOw5c/BZSzqPceTOJWDMxBRwcFHmKTRgf030062SNyLE2CE77yZOVMs0AFdeeSUNDQ3cc8891NbWsmDBAl577bXOpNaqqqpuMyHf+973UBSF733ve5w+fZqCggIuu+wyHnjggcR+J2LUSVRJb9clmjNVNGu6HaPpDDQXnUPRqX/yCet+/uWYwL4aB8sm5kb1WA5vAFXV0OmUhIxdCJGagiGVnafssuFdgkQdjADcdNNN3HTTTX3et2HDhu4PYDCwbt061q1bF9sIxZjV6Iz/Sd7k8nG0sQ0FWFZmInfP+wA0lK3ufWzxuRSd+idLg9uAi9lXHW4NryiRBxaqCk5fkKw0Y9xjF0KkJn9QZcfJVtmTKoFkbxqRkjRNoykB+9F0tH+fVmRjkv0j9Kofd8YE2jKn9Dq2o99IsWsf+ToXzW4/dY7oxyB5I0KMXrLz7tCQYESkpESV9G4+Hs4XWVqRQ8HpcBVNfdlq6GO2w2ctxpU1DQWVz2QfBmBvtT2msQshRh/ZeXfoSDAiUlIiSnobnD6ON7lRFFg8LoP86vASYkPpRf2e01R0LgCrjbsgxtbwEowIMfrIzrtDS4IRkZISUdLbUUUzo8jGeNcOjAEHfnMu9ryF/T9uSXipZrZnC6BRWeckGGVPZ28ghDcgL1hCjBay8+7Qk2BEpJxElfRu7trorKOKpvRC0PVfUt6av4SQ3kK6v5FF5mq8AZVjDW1RP7asJwsxOsjOu8NDghGRchKxMV6dw0tVsxudAgvLsyioDjfq61nS25OqN9NSuByAyzPD+y3FslTTKsGIECOe7Lw7fCQYESknEcFIRxXNzOJMSr2HsLhrCBqsNBf23Sm42+O3V9Wco+2AGPepkbwRIUY22Xl3eEkwIlJKwkp626tollTkdDY6ayo+D9VgGfTcjmBkQttOrHg51tSG2x/dspGzvfmZEGLkqWpyt/cZSvZIxg4JRkRKSURJb63dy8kWD3pFYeH4nC5dV3s3OuuLO6MCd/o49FqQizMOoWlwoDa61vCqSkLyXoQQw+twvey8mwwSjIiUkoiS3s3tVTQzS20UBKqx2StRFT2NJRdEdgFFobl9duTjafsA2CdLNUKMepW1To43Rp+wLuInwYhIKQkp6e2ooplwZi+aloJlBE1ZkY+jPRhZFNgGMSexyp4VQowEsvNu8sW0N40QQyERJb3VrR5Ot3rQ6xQWjs+m4L1wMNJY1n+js740F56FqhjI8Z5kolLLMWcxDU4fBTZzxNeQmREhhpaqagRVjZCqEVTV9j+1M3+GNEKaRkhVCaoawZDW4xi127EieSQYESkjEVU0He3fZ5dmkqXayW4Kz2zUl0aWL9IhZMzAnr+InIZNXJ5ZyaP2YvbVODjfVhDxNXwBFW8ghMXYf18TIcYSTesdLPQZRPQXYITO3K5qmlS6jCISjIiUEW8womlaZ0nv0gm5FFS/jaKpOHJm40svjfp6jcXnkdOwiVWG3TzK+eFgZFrkwQjtsyMSjIiRStN6ziT0CBJCGqrWd7DQ+zxVggfRLwlGREpIRElvdauXGrsXg05hQXk2BZvaG51FOSvSobn4PNj9CDO82zESZH+NA1XV0Ol6b7LXH7snQFHm4OXEQiTTwTonTm+g1zKGNPsSw0USWEVKSERJb0cVzZyyLNJ1PnLr3oOOXXpj4Myegc+SjzHk4WzjIdz+EMebo8u0b3VL3ohIbW5/kKomNy1tAZzeIG5/CH9QlUBEDCsJRkRKiLekV9O0znyRpRNyyKt9D33Ihzu9nLasabFdVNF17uL7qYz9EEOJr8snzc9Eaqtu9SR7CEJIMCJSQ7wlvadaPNQ5fBj1CvPLs7s3OlMiX1bpNa72Et8V2naIocRXVcHhldkRkZpUVaO61ZvsYQghwYhIvkSU9HYs0cwtyyJNr1FQ8zZE0XW1P81F56ChUOI9QiEtHGlowxsIRXUNKfEVqarR5cMflKxSkXwSjIikS0QVzeb2RmdLJuSS3bgFo9+O35xDa96iuK4dsOTiyJ0DwKXWfYRULepW0ZI3IlLVKVmiESlCghGRdPEGI1XNbhqcPkx6HfPGZZ1Zoim9CHTxl9U2FYWXai4274EYlmpkZkSkIo8/RHMCevsIkQgSjIikSkRJb8esyLxxWVgMuu7BSAI0lYSDkfn+7ehQo05i9QdVPP7olnaEGGqnZVZEpBAJRkRSxVvSG250Fs4XWVKRg611P2nuakL6NJqLzknIGB258wkYbaSFHMxXjlJt99Liju4TpcyOiFSiaRo1dglGROqQYEQkVbwlvceb3DS6/JgNOuaWnVmiaSo+F9WQmGZjms5Ac9HZAFzWUeIrSzViBGtw+fAFJHFVpA4JRkRSxVvSu6W9t8j8cdmYDXoKqruU9CZQR4nvBfqdEEO/kdYoZ1KEGEqnW2RWRHTnirOiMV4SjIik8QbiK+nVNI3N7XvRLKnIweI6ia31AKqip7HkggSO9EwwUuE9QCYu9tU4ULXIl5dcvqB0tBQpwRsI0dwmwbEIC6oqL26u4spfbkzqDK4EIyJp4n1BPNrYRnNbeIlmTmkWBdXhvWhaC5YSMOckaJRhPmsJrswp6FC5wLgPpzcY1adLTQOHLNWIFHC61UMUcbQYxeyeAI++cZA399dT1ezmrQN1SRuLBCMiaeIt6d3SXkWzoDwbk0FHYWcVTWKXaDo0FYdbw3/Cuhckb0SMQJqmSft3AcCRBhf3vbyPg3UuLEYdj105n8sXjkvaeCQYEUkRb0mv2qWKZmlFLkZfM9mNWwBoKEtMSW9PTcUrAVge2gFoUeeNSDAikq3R5ZfE1TFO0zTeOdjAw69X0uoJUJxl4a5LZ3LRzKKkjsuQ1EcXY1a8Jb1HGly0uAOkGfXMLs0kv+pvKJqKI3sW3vSyhI61Q2vBUkJ6C1nBBqYppzhYryMQUjHqI4vpWyUYEUkmvUXGtkBI5YWPqnjvcCMAi8Zn81/nTMRijL85ZLxkZkQkRbwlvV2XaIx6HQWn34AhnBUBUPVmWgqWQns31kBI41CdK+LzA0EVtz+5Geti7PIGQnFXr4mRq8nl40evHeC9w40oClyxqIwbzp+cEoEIEoyIZGmM40UxvEQTDkaWVuSgC3rIq3sfhqCkt6eOpZo1JmkNL0aWaklcHbP21zi475X9HG9yk27Sc+tFU7lkTglKHDuaJ5os04hh5w2E4qppP1zvwu4JL9HMKskkr3Y9+pAXT/o4XFkzEjrWnjqSWGcG9pCGN6ZgpCQrbWgGJ0Q/womr3mQPQwwzTdN4fW8df9l+Ck2D8blWvn7BZPIzzMkeWi8SjIhh1xRnSe/m9kZnC8dnY9DrKDgdLultKF0NQxzpu22T8FjLSHOfZrluPxuaF+L0BrBZjBGdLzv4imRoavPjDcj+SGOJNxDiNx8cZ2v7LPLZk/P4z+UTMBlSc0EkNUclRrV41q1VVet8ci2ryEVRg+RXvwVA/RAv0QCgKJ0N0D6etg+iXKpp8wUJhqSaQQwv6bg6ttTavfzw1f1sPdGCXqfwxeXjuebsipQNRJBgRAw3TdPianZ2sN6JwxvEatIzo8RGVuM2TP5W/KZs7PmLEjrW/nQs1azU7YAoW8NrGjiS3HZZjC2+YCiuHC0xsmyvauGBV/dTbfeSlWbkOx+bzqrphSmVH9IXWaYRwyrekt7N7VU0i8bnYNCdaXTWWHohmm54fp2bC1egKgaKAqcpV+rYV2NE07SIn+x2T4DcdNOQj1MIgOpWrySujgGqqvHSzmpe3l0DwNTCDK4/fzJZaZEtISebzIyIYRVPSW9I1dhWdaaKBk0bso3xBhyHyYY9fyEAq/S7aXEHqHVEnhwoFTViOEnH1dHP5Qvy+FuHOgORi2YU8q2PTRsxgQgSjIjhFs90cWWtE6c3SIbZwIziTDLsB0hrO0VIb6Gp6JyEjnMwTUXhpZqLLeHW8HujWKqRYEQMlyaXD49fEldHs6pmN/e/so891Q5Meh3XnjuRzy8bj0E3st7eR9ZoxYgWb0lvR/v3ReOz0esUCtqXaJqKz0U1DG+5bFNJuN/IInUXRoJRJbFK8zMxXKTj6uj24dEmHvrnARpdfvIzTNxxyQxWTMpL9rBiIjkjYtjEU9IbVFW2VbVC+140wJBvjDcQZ/ZM/OZcLL5mFusOsr12DkFVjfjTSKs7gNUkTz8xdCRxdfQKqir/u+UUbx2oB2BOaSZfOW8SGeaR+5oiMyNi2MRT0ltZ68TlC2KzGJhWZMPSdgpb6340RUdj6QUJHWdEFF1nVc1q4x58QZWjDW0Rny5LNWKo1bR6UaWKfNSxewI88q+DnYHIJ+aWcPOFU0d0IIIEI2K4xFvS21FFs3h8TvsSTbi3SGv+EgLm3ISNMxodreEvMu6GKPuNSDAihpokro4+h+td3PvyPg7Vu0gz6rnxgsl8emEZOl1ql+1GYmSHUmLEaHXHXtIbDKmdVTRLKnIAOqtohqXRWT86kmYnBo9QQCv7qtP59ILIdgzuaH5miHDHXyGi0dzmxy2Jq6OGpmlsqGzgxS0nCakapVkWvn7BFIqzLMkeWsLIK6EYFk1tsS/R7K914vaHyEozMq3QhtHXQnbDZgAaSodul97BBCx5OHLmAHCebhfHmtpo80WWmKppMjsiho50XB09/EGV33xwnBc2VRFSNRZPyOG/L505qgIRJBgRwyWe/iIde9EsHp+DTqeQX7MBnRbCmT0Db0Z5AkcZvY68kbWWvWgaHKh1RnyuBCNiKPiDKg0u2RRvNGh0+XjotQN8cKQJRYHPLhrH9SsnYTHqkz20hJNgRAy5eEp6AyGV7e1VNJ1LNEmsoumpY5+as9mJDlXyRkTS1dg9krg6CuyttnPfy/uoanaTYTZw2+ppXDynOOXbusdKckbEkIunpHdfjQNPIER2mpEphRnogh7yav8NCcoX0esUTAZdzI2h7HkLCBozsAUczFGOsa8m8n4nEoyIoSC9RUY2TdP4555a/m/HaTQNKvKs3HD+ZPIyzMke2pCSmREx5OIp6e1YollSkYNOUcit+wB9yIvHWoYre2bcYyvKtMS19qrpjDQXrgBglX4XDU4fDc7Ivt9gSIs4x0SISLS0+XH7JHF1pPL4Qzz1zhH+uj0ciJw7JZ/vXjxj1AciSDAihlo8Jb2BkMqOk+1LNBN6NDoruwgSMF05LjeN4sz4EsE6lmrWmPdAlCW+rTI7IhJIZkVGrhq7hwf+uZ9tVa3odQpfOmsCX14xAeMYqbgbG9+lSJp4Snr3nLbjDajkWI1MKkgHNUR+Tbi/SCI2xstMM5JpMZJuNpBhiX3FsiMYmRmqJJM29kWzT41bghGRGIGQSr1TEldHom1VLTzw6n5q7V5yrEZuXzud86cVjNr8kL5IzogYUvGU9HY0OltSkYtOUchu3IbJ10LAlEVr/pK4xzYu50x+R1GmBZfXFdN1vOlltNkmke48ytm6vbxbm4mqahE1IpK8EZEo0nF15FFVjb/tOM2re2oBmFaUwddWTh5Ru+0misyMiCEVa0mvP6iy81T7XjQTelTRlKxC08UXRxv0CkVdlmfiX6pp78Zq2IXbH+J4U2St4dt8QQIheQcR8ZMlmpHF5Q3y2PpDnYHImplF3LZm2pgMRJBgRAyleEp6d5+24wuq5KWbmJifDpp2JhhJwBJNaXYa+i4zF2kmPZlxvAh09BtZZdgFaFHljThkdkTEqdXtl2ToEeREUxv3vbKPfTUOTAYd1503kSuXlke80eZoNHa/czHk4inp3XLiTBWNoiik2w9ibTtJSG/ufOOPR1l27xLcoszYM9ZbCpYR0pvJVxuZopyWJFYxrE5Jx9UR4/0jjTz02gGa2vwU2MzceckMlk/MS/awkk6CETFkYi3p9QVC7DxlB2Bpjyqa5qJzUQ3WuMaVk24ivY8dLoviWKpRDRZaCpYBcL5uJ0ca2vAGIiuxlLwREY9ASI24nFwkTzCk8sJHJ/jN+8cJhDTmlmXxvUtnUp4T3+vZaCHBiBgSqqrFPDOy+7Qdf1ClIMPMhLzwEzWRG+OV5/TdmMxi1JNtjX2pprkoPGOz2rSHkKpxsC6y1vAOTwBNi63iSIhau5eQKr8/qazV7efH/6rk7coGAC6bV8I3LpzS54eisUqCETEk7J4AoRhLejefOLNDr6IomNuqyWzZi6boaCxZFde4zEYdBbb+l2PimR1pLAmX+C7W9mHBx94IS3yDIY022WFVxEiWaFLboTon972ynyMNbaQZ9Xzjwil8akEZujFUthsJCcvEkIi1pNcbCLG75xJN+6xIa94iApbcuMZVmp02YO1+YaaZg3VOYpmocNsm47WWYHHXcJZuP/tqsiI+1+4JkCGfkkSU7O6AJK6mKE3TeOtAPf+75RQhTaMsO42vXzA5rg88o1lMMyNPPPEEFRUVWCwWli9fzqZNmwY8vrW1lRtvvJGSkhLMZjPTpk3j1VdfjXXMYgSItaR31yk7/pBKoc1MeW54OSVRVTSK0nfialdmg55sqynmB2hsb4C2UreLGrs34u6zre7Yk33F2HWq1Z3sIYg++IIhfv3+Mf6w+SQhTWNZRS7/fcmMlA1E8mrexfzabSSzUU3Uwcgf//hHbrvtNtatW8e2bduYP38+a9eupb6+vs/j/X4/a9as4fjx4/z5z3+msrKSZ555hrKyskSMX6SgeEp6N/eoojH47WQ3bIYEBCP5GeaItt6OZ6+ajm6sq427AdgfYVWNJLGKaAVDKvUOSVxNNQ1OHw/98wAbjzajU+BzS8Zx3XkTMUfw2pMMRVUvM/+96zHt+C1sez5p44h6XvjRRx/luuuu45prrgHg6aef5pVXXuHZZ5/ljjvu6HX8s88+S3NzMx988AFGYzg5sKKiIhFjFykq1sRVj7/LEk1FeDkmv2YDOi2EM2s6nozxcY1rXD+Jqz0V2sxU6mL7kNBceDaqome8dppxSgP7anI5Z0r+oOe5fSECIXXM7EMh4lcjiaspZ89pO7/891Hc/hA2i4GvrZzEjOLMZA+rX+MO/Y7p2+9HQSMw+7MYF/5n0sYS1Suf3+9n69atrF595hOqTqdj9erVfPjhh32e89JLL7FixQpuvPFGioqKmDNnDj/84Q8JhfpP2PP5fDgcjm5fYuSItaR356lWgqpGcaaFcdmJXaKxmvQR73xp1OvIiXGpJmSyYc9bAO1LNftqHKgRJqDI7IiIhnRcTR2qpvHyrmp+tv4Qbn+Iifnp3P3xWakbiGgak/b8jBnb70NBo2rKl/B+4knQJ6/7a1TBSGNjI6FQiKKiom63FxUVUVtb2+c5R48e5c9//jOhUIhXX32Vu+++m0ceeYT777+/38d58MEHycrK6vwqLy+PZpgiieIp6d1yvHsVjS7oJa/235CAYKQswlmRDvEs1TS3L9Ws0u/C6Q1GXO3QKpvmiQjZPYGYl0JFYrn9QZ7ccIS/7ahGA1ZOzef2tdPJTY8x92yoaSrTt/2ASfueAODI7Fs4uPB7oCR3VnbIH11VVQoLC/nlL3/J4sWLufLKK7nrrrt4+umn+z3nzjvvxG63d36dPHlyqIcpEiTWkl63P8ie6u5LNLn1H2IIuvFaS3Bmz4p5TDodlGRFF4wUZJiJtTNzRxLrOfo9GAhGvIuvzIyISJ2Wct6UUN3q4YFX97PjZCsGncJVKyZw1YqKlF1uVUJ+5mz8FuVHfo+Gwv5F3+fY7BvD2f1JFlXOSH5+Pnq9nrq6um6319XVUVxc3Oc5JSUlGI1G9PozyTszZ86ktrYWv9+PydQ7ejSbzZjNsbfmFskTa0nvjpPhJZrSLEtnxUvB6TcAqC9dHdeTpdBmwWSI7sXBoNeRl26OqbOlM2c2fnMOVl8Li5RD7KvJ5eI5fT8/unJ4w83PxtK24SJ6wZBKndOb7GGMeVtONPOb94/jC6rkWI3ccMFkJuVnJHtY/dIF3cz74Bvk1/4bVWdk77KHqRv/8WQPq1NUr9Amk4nFixezfv36zttUVWX9+vWsWLGiz3POOeccDh8+jNolG/DgwYOUlJT0GYiIkS3Wkt4zSzTtfUTUEAXVb0EClmhibbcc81KNoqOpvRvrSv0uDtY58QcHz4YNhTRc0jNCDKLW4Y25oaCIX0jV+PPWUzz9zlF8QZUZxTbu/vislA5EDL5WFr1zNfm1/yakT2PHuf8vpQIRYlmmue2223jmmWd4/vnn2b9/PzfccANtbW2d1TVXXXUVd955Z+fxN9xwA83Nzdxyyy0cPHiQV155hR/+8IfceOONif1ORNLFWtLb5guyt70EdsmEHACymrZj8jUTMGXRWrAk5jFlWAxkxdjiPT/D3G1n32h0lPheaNhFUNU4VB9Za3hZqhGDkSWa5HF6Azz25kFe2xvOkVw7q4hvrp4W147fQ83srmXJ218ku2kHflM2Wy94nuYEbDaaaFGX9l555ZU0NDRwzz33UFtby4IFC3jttdc6k1qrqqrQdVlsLy8v5/XXX+eb3/wm8+bNo6ysjFtuuYXvfve7if1ORNLFmri6/WQrITXcobC0fYmmY2O8xpLz0XSxP9EjLefti16nkJ9hps4R/ZR4x5N9FsfIx86+GgezSwfvyNrqDjAuJ6bhijHA4Q3glMTVpDje2MaT7xyhuc2P2aDj6rMrOvPbUpXVeYyF7/wXae7TeNOK2L7yN7RlTUn2sPoUU//pm266iZtuuqnP+zZs2NDrthUrVrBx48ZYHkqMII0x7hy65Xi40dnSivZ3YU3r3BivoWxNzOPR6xWK4+x4WJQVWzDit+TjyJ5FZus+ztXtZmt1MSwe/DyHzIyIAcisSHK8d6iR//noBEFVo8hm5usXTIm6Qm+42Vr2svDdazH5mmmzTWT7ymfxpqdus9HUTPkVI46qajTH0NLc5Q2yvya8hNGRL5LuOIzVVUVIZ+rMvYhFSZYFQ5xZ7fnpZvT6+JZqztfv5GSLJ6JAw+0PRZRfIsaekKpRG0NgLGIXCKn8buMJnvvwOEFVY/64LO76+MyUD0Ry6j9i8dv/icnXjCNnNltW/T6lAxEkGBGJEmtJ77aTLYQ0jfKctM5ZjI5GZ81F5xAypsc8pnExJq52pdMpFETYLK2njmDkAv1uFFRpDS/iIomrw6vF7efHr1fyzsEGFOBTC0q5cdUUrKbU3tCy4PQbLHj3WgzBNpoLlrP1gt8RsOQle1iDkmBEJERjjF1XO6pouq69dpT0xlNFk201JmwX3Firaux5Cwga0snBwWzleGeS7qDnSTAi+iBLNMOnstbJvS/v42hjG1aTnm9cOIXL5pWiS/Gy+5Jjf2HeB99Ar/qpL1vDjpW/ImRM3SqfriQYEQkRS0mv0xvgQG17FU17vojZXUNWyx40FBpKL4x5PImYFemQazVhjLJPCYCmN9FcFC55P1+3i/01DrQIWsPbPbKDr+jO6Q1IPtEw0DSNN/bV8cgblTi9QcblpPG9j89k3rjsZA9tUBMO/IrZm+9E0VROT/wsu1f8DFU/cvp1STAi4uYNhGiLoT/GtqpWVA0m5FkptHUs0YR72LTmL4p5atFo0FFoS9yTMK6lmqKOvJFdtLgD1NgHX/N3eIIRBS1i7JB9aIaeLxDiV+8d449bTqJqsHxiLndeMqPztSllaRpTdj7M1F0PA3B8xnXsX/IAmi61l5N6GlmjFSkp1pLeze1VNB29RYAuVTSxL9GUZVvQxdgfpD/FWRaqY3hDaGov8V2sO4gNN/tqHJ3ly/0JqRpOX5BMS+r2LhDDJ6Rq1EYQxIrY1Tu9PPH2EU63etAp8Lkl5Vw0ozDluyErapAZW++h7NifATg477tUzbg22cOKiQQjIm6xlPTaPQEq69qraCaE80UMfjs59ZsAaCiNLRhRFCjLTtwSTYccqxGTQRd1pYs3o5w220TSncc4W7eXfdUlrJ5ZNOh5dndAghEBQJ3DS3CMJq5qmkZQ1fAHVfwhtfufQZVA+999oe7/7u/4nrcFQuFre4MhNA0yLQauP38y04psyf7WB6UL+Ziz8ZsUnn4TTdGxb8kD1Ey8ItnDipkEIyIusZb0bjvRgqbBxPx0CtqXVPJr3kGnBXFlTsVjmxDTeHLTTaSZ9BEcGR1FUSjKtHCy2R31uU3F55HuPMb5up38oG45wZA6aMmx3RNA9qoWtG/GlmpUTev2xt/xpu4Lhbr9u8+goEdwEGj/09fj3x3HDFcYNqUgg6+dP4kca+pvU6IPuJj/3g3kNnxESGdiz4qfxtWTKRVIMCLi0hpjSe/mE30s0bTni8SzRJPIxNWeijLNMQcj4w/9lgsMu/hvb4ijjW2DfvKSihoB4PIFaXUP7e9Cjd3D+4eb8AVD4YCg6wxDfwFEEmZqdAqYDDpMel34z65/1+swtv9pbr/PqO99TF/nmQzhc7LSjCm/LANg9Dax8N9fIbNlL0FDOjvPfZqWwuXJHlbcJBgRcWmKoaS31e3nUJ0LugQjupCPvNp3AaiPMRixGPXkZwzdp5psqwmLUY83EIrqvJaCZYR0JkrVRiYr1eyrLh00GPH4Q/iCIcyGxM/yiJFjqMt5VVXj6XeOxpUga9ApEQUJ/QUCPe8zGhTMen17QKF03mfQSb2Fpe0UC9/5L9Jdx/GZ89ix8tc4c2Yle1gJIcGIiEssJb3bqlrRgEn56eS1V6nk1H+IIdiGN60IZ86cmMZSlpM25J9sijLNnGiKbnZENaTRWrCUvLr3OV+3i7drpvHphYN3Q7R7AhTaJBgZq1RVo8Y+tMHIxmNNnG71YDXpuWhG4SBBwpmZh66zDolOFhd9S7cfYuG712Dx1OOxlrHt/N/gsVUke1gJI8GIiFmsJb2bO/ei6drorEsVTQwBhU4HpdlDX4JXmGmJOhihvaomHIzs5DdNl9DmC5I+SFM2hyeQ+mWFYsjUOYc2cTUQUvn7jmoALplTzCVzSobssUR8shq3s+C9r2L023FlTmX7yl/jsxYne1gJJfNeImaxdF1tcfs5VB9eolnckS+ihuLOFynIsAzLkkZWmhFrDAmyTcUrAViuP4BJ83Og1jnoOUOdKyBS21Av0bxzsIGmNj/ZaUYunFE4pI8lYpdb+28WvXM1Rr+d1rwFbFn1wqgLRJBgRMSjKYYlmq0nwu3fpxRkkJsezu/Iat6J2ddEwGijpWBZTGMZN4wbVxXGsBNwW+YUvGnFWPCzXLeffRG0hnd6g6jq2CzpHOvahjhx1RsI8cruGgAum18quUkpqqjqZRa8dz36kIfG4vPYdv5zBM2p3w02FhKMiJjEWtJ7ZommaxVNeImmsWQVmi763hrpZgM56cNXjleUGUM3VkXpbIC2UreLfdWDByMdzc/E2DPUHVff2FeH0xukyGbmnCmpv4naWDTu8AvM2fgtdGqA2vGfYOc5T6Eahq5aMNkkGBExiaWkt8nl40hDG0rXJRpNozDOjfGGc1YEwGYxYjXHvlRzvm4XDS4f9c5IWsPLUs1YE05cHbqOq05vgNf31QLw6YVlUqWSajSNiXt/zoxtP0BB4+SU/2TP8p+g6VO//0k85LdQxCSWkt6tVeElmqlFGWS3NxZKdxzB6jpBSGfqnDmIhl6nxLyrbjyKY1iqaS5agabomKo7TSmNEc2OSL+Rsafe6SMQZaffaLy6uxZvQGV8rvXMhwKRGjSV6dvvY/LenwNwZPY3qFx4Nyij/6169H+HYkjEUtK7+Xg4GFk6oWsVTXhWpKVoRUxbXRdlWjAO0s10KBTFEIwETVnYcxcAsFK/K6K8EUliHXtOt0ZfrRWpJpePtyvrAbhiURm6EdDka6xQQn7mfPRtyg//DxoKBxbew7HZ34ipunAkkmBERC2Wkt5Gl49jjW0oCizqtjFeuIqmPsa9aMblDu8STYd0s4EMS/SV8U3F7bv46nZyoNZJaJAEVW8gFHWTNTFyuf1BWtqGLgB9aWc1QVVjRrGNWSWZQ/Y4Ijq6oJv5799AcdXLqIqBPWc9wqmp/5nsYQ0rCUZE1GIp6d3SPisyvchGVlo4SdXsriWreRcaCo2lF0Z9zcw0Y1I3k4tlqaZjKepc3R78fh/Hm9oGPUfyRsaOoSznrW718MHRJgA+s7BsRLQ+HwsMvlYWvXMN+bX/JqS3sPPcp6kb/4lkD2vYSTAiohZLSW/HXjRdG53lV78FgD1vIf60gqivOdyJqz3FslTjyJmD35SNTfGwQDkc0VKN5I2MDaqqUT2Eiav/t/00mgYLx2czqSD6JVGReCZPHUs2/CfZTdsJmLLYdv7zNJWsTPawkkKCERGVWEp6651eTjS5w0s048/UyBd27boaJYNeiSkYSKQ0k57MtChnZnR6movOAeB8fWQlvhKMjA0NrqFLXD3S4GL7yVYUBS5fMPhWBGLopTmPs/Stz5NhP4g3rZAtq36PPX9hsoeVNBKMiKjEUtLbsUQzo9iGrX1ZxeB3kFO/EWLcGK80Ow19CuyJEdNSTUk4b2SlbhdHG9oGzQlxeAPS/GwMGKreIpqm8ZdtpwA4e1IepdnJnVEUYGvZx5K3Pk9a2yncGRPYcuGLtGVNTfawkkqCERGVWEp6Oxuddamiyat9F50WxJU5JabNnpK9RNOhMIYGaE1F4byRubpjZGl2KusGbg2vquFurGL08vhDNMew/BmJvdUODta5MOgUPjm/dEgeQ0Quu34Tizf8J2ZfE47sWWy58A9408cle1hJJ8GIiEq0Jb21Di8nWzzoFFg0vmvX1dgbneWkm7CaUmOPR4tRT7Y1uqUaf1ohzuwZ6NA4V7ebvbJUM+YNVTmvqmn8dftpAFbNKOzcJVskR/7p9Sx8978wBFy0FCxj6wW/w2/JT/awUoIEIyJisZT0bmmfFZlZktlZCquE/OTXvAsxBiPlKTIr0iGW3JXOEl/9TkliHeNUVaO6dWgSV7ccb6Gq2Y3FqOPSOaNvc7WRpOTYX5n3wU3oVT/1pRexfeWvCZlsyR5WypBgREQsppLeE70bneXWf4gh2IY3rRBHzpyormc26iiwpdanu8JMc9R9ibq2hq+zu2luG3jGqdUzNFP4IvkaXT78Q5C4GlRV/rYjPCuydnZxZ76WGH7jK3/N7M13oNNCVFdcwe6zf46qT63XsWSTYERELNqS3hq7h1MtHvQ6hYVdqmgKTocbnTWUro66zXFpdlrK9UcwG/RRb9TXmreQoCGdfMXBLKVq0NkRX0CV5mej1KkhSlx9/3AT9U4fNouBNTOLhuQxxCA0jSm7fsy0nT8C4Pj0r7Bv6Q/RdKmxzJxKJBgREYmlpLejimZWSSbp5vYnn6Z2dl2NdolGUaAsRSsBol2q0fQmWgrPgvZurFLiOzZ5/CFaBpkVi4UvGOKlndUAfGJuCRZj9Bs7ivgoapCZW+6i4sAzABya9x0Oz799zLR3j5YEIyIisZT0djQ6W1JxJnE1s3kXZm8DQWMGLQXLorpefoY5ZV9UC21mot38tKMba8c+Nao28M9XgpHR53Srh0H+22Py1oF67J4A+RkmVk6LvqGgiI8u5GPuh7dQduzPaIqOfUse4MSM65I9rJQmwYiISLQlvadbPVS3ejHoFBaWd210Fq6iaSy5IOotsVOlnLcvRr2O3PTo1oA7klgXKwfB5+BU88DT9RKMjC6aplFjT/wSTZsvyD/31ALwqfllSdlIcizTB1ws+Pd1FJ5+A1VnZNeKx6me9P8le1gpT35LRUSiLentqKKZXZrZrQy3oL3rarSNzqwmfcqXJRZF2XPEkzGetowKjEqIs3V72VtjH/B4pzQ/G1UaXD58gcQnrv5zTy1uf4iy7DSWT8yN4AyRKEZvM4s3XEVu/UaChnS2r/w1DeM+luxhDcpk0CU9aJVgRAwq2pJeTdPY3F5Fs6TLXjRWxxHSncdQdcbOapJIlaXwrEiHgozol2qaO5ZqdLsGTWJV1XA3VjE6DMWmeC1uP+sP1AFw+cIydCnQpXisMLdVs+Ttz5PZsge/OZetq/6nMy8slZVkW1gxOS/pS+ASjIhBRVvSe6rVQ609vESzYFzXKprwrEhz4QpCxsg36tLpGBEtrA16HflRzt40dvQb0e3iUJ1z0BJPWaoZHbyB0KDl3LF4eVcNgZDGlIIM5o/LSvj1Rd/S7YdZ+taVpDuP4bGWsuXCP+DMmZ3sYQ0ozaRn4fhsZpdmJX1WBAlGRCSiLentqKKZW5ZFmulMtB3rxniFNktKPFkiEW1VTUvhclSdkXJdA+VaDYfqB24NL8HI6DAUiat1Di//PtQAwBWLylKuBH60ymzayeK3v4DFU4crcwpbLnwRt21isofVL0WB8XlWzpqUl1JL3yPjFV4kTbQlvZqmdeaLdK2iMXnqyGreiYZCQ+lFUY2hPMca1fHJlJ9hRq+P/E1ANVhpyV8CEZb4SjAy8mmaRs0QdFz9247TqBrMK8tiapF09hwOubXvseidL2Pyt2LPnc/WVS/gs6Zup9sMi4ElFblMK7KlxEajXUkwIgYUbUnvyWYPdU4fRr3C/K5LNNVvA2DPm48/LfJSwwyLgawo935JJr1OoSDKTxtnurEO3hpemp+NfI0uf8L/D080tbH5eAsKcPmisoReW/St8OSrLHjvaxiCbpqKzmXb+c8RMOdEcObw0+lgUkE6yypyyUpLzddTCUbEgKLNF+noLTKvLLtbQlTnxnil0S3RpHI5b3+i3cm3o9/IWbr91LfYB539aHXL7MhIdnoIOq52bIa3bGLuiJpJHKnKDv+euR9+E50aoLb8Unac+zQhY3qyh9WnbKuR5RPzmFSQkdIJzRKMiAFFE4yEl2g6qmjOfELQB1zk1m+EKPNF9HqFkqyRF4zkp5sxRLFU05Y1DW9aIWmKn6W6SvYPMjsiSzUjlzcQirpnz2AO1DrYW+1Aryh8eoHMigwpTWPivieZue37KGicnPwF9ix/JOqeScNBr1eYXmxjSUXumQ7YKUyCEdEvbyCE2xf5dPKJZjcNLh8mg455ZWcy+fNr3kGnBmizTcKdOSni65VkWVJuXTMSOp0S3WZ+ikJzUUdVzeBLNRKMjFzVCU5c1TSNv24Lz4qsnJafcptIjiqayrQdDzB5z2MAHJ11I5WL1oEu9bpC52WYWDEpj/LckTNLJsGI6FfUSzTHO5ZosjB3W6KJrdHZuBE83RxtVU1jSTgYWanbxb5qB9oA71guX4CQND8bcTRNozrBiavbT7ZytLENk0HHJ+aVJvTa4gxFDTD7o9sZf+i3AFQu/B5H59yScvvMGA065pRlsXB8TtL7hkRLghHRr2i6rnZdolnapdGZEvKTX7MBgIayNRFfL9tqJGMETC32Jy/dhNEQ+dOrufBsNHRM150izVNLjb3/Ny1VDXdjFSNLU1tiE1dDqsb/teeKrJlZlLKJiSOdLuhh/vtfp6TqJVTFwJ7lP+Hk1KuSPaxeirMsrJiUR3FWdB+EUoUEI6JPqqpFtZvosaY2mtr8mA065pRldt6e0/ARhmAbPkshjty5EV9vJM+KACiKQmEUU+ZBczb2vHkAnKffxd5BSnwliXXkSXTH1Q+PNlFj95Ju0rN2dlFCry3CDH47i969hvyadwjpLew890lqJ3wy2cPqxmLUs2B8NnPKsjBF8QEo1YzckYsh1eL2R7UUsLl9VmT+uGzMhr4anV0ESmS/biaDLqo38lQV7VJNU7HkjYxWvmAo6mXPgQRCKi/tqAbg0rkl3fZ/Eolh8tSz+O3/JLtxGwFjJtvO/w1NJRcke1jdjMtN46xJuVF3fk5FEoyIPjVFMSuiahpbO5doutTZayoFp9dDlPkipdmWlC5Bi1SO1YjZGPlTrKPfyHm6PRypayUY6r81vAQjI0t1qzehiatvV9bT7PaTYzWyanph4i4sAEhznmDJW/+BzV6Jz1LIllUvYM9fnOxhdbKa9SypyGFGcSaGEdKdejCj47sQCRfNp7hjjW00uzuWaM5U0WQ278bsrSdoSKelYHlE11IUKMse2Us0HcJLNZHPjjhy5uA3ZZOpuJkROsSRhrZ+j/UHVTx+aX42UlQnsLeIxx/i1d21AHxyfumInppPRRkt+1jy9uextp3CnTGezRf+gbbs6ckeFrQ3L6vIT+esiXlkW1OvnDge8lssevH4oyvp7aiiWVCe3W0PmY4qmsaS8yOuw89NN3Xbz2akK4qmAZpOT3PR2QCcr5elmtGiyeVLaOD4r321uHxBijMtnD05P2HXFZDdsJklG/4Ts7cRZ/YMtlz4It6M8mQPC4DMNCNLK3KZUpjazctiJcGI6CWaWRG1nyoagILq6DfGG+mJqz1lW01Rldg1ddnFd2+1fcBjWz2J3/VVJF4iO646PAH+ta8OgMsXlo3IPjypKr/6LRa++18YAi5aCpay5YIX8FuSH+zpdQrTimwsrcjBZhm9FVOS9SR6iSZf5Ei9i1ZPgDSjntmlZ6porI6jZDiOoOqMNBafH9G10kx68jNG19Qj7bMjJ5rcER3b0Rp+nnIUR1MtLt+0fkuc7VJRk/ISnbj6yu4afEGVijwri8ZnR3CG6EbT0Kl+dEEP+pAXXciLPuQhq3E707ffh04L0VB6IbvPegzVkPwS2Zx0E7NKMkfVbHF/JBgR3URb0ttRRdNriaY6nLjaXHgWIVNkO4iWZqeNym3Pi7IsEQcj/rQinFnTsdkrOVe3hwO1C1gyIbfPY12+ICFVk0/HKaym1Yvafx5yVBpdPjYcbADgikXjRtVzRVED6IJe9O3BgS7kC//Zfpsu5EEf8vUIItr/7PO2jut0v00X8qLQfyZxdcXl7F/yAJouuW+NBn14NqQ0e+RthxErCUZEN9GU9KqqxtaqPqpouuSLRLoxnk4XrqIZjTItRqwmPe4I8waais/DZq9kpX4XL1U7+g1GNC08bZ+TPvpmk0aLRCau/n1HNSFVY2aJjZklmRGcESdN7fYmrg95er2x60OeLn8/c9uZIKJ70NDx7+63+dBpwaH/fnpQdUZUvZmQPo2QIY2aCZ/i2KwbI25BMFQKM81ML7Z1a5EwFkgwIrqJZonmUL0LuyeA1aRnVpcXR5OnnqymHQA0lF0Y0bUKbZZR/eQrzLRwvLH/6piumorPo6LyV6zU7eKhajuapvX7KbhVgpGU1dzmjzgAHcypFjcbjzYB8JmF4xJyzQ4Fp99k4r4nMQSc3Wcg1OHPSdJQCBnSOoMEVW8hpLegGtr/1Ke1/2kh1Hmb5cyxho6/m9uv0/U8M6qh43wzmi618i/MRh3Ti21RVeCNJhKMiG6iWd/eciJcRbOwPLtbrXtB9VsoaNhz5+NPi6wzZNkon44szoo8GGnNX0xQn0ZhqJUC92HqndP7baAmFTWpK5EdV/9v+2k0YPGEHCbmJ26r+ozWA8z58NZBA4+QztTljbxLENDzNkP7m363oOFMcHEmIOgacJjbr5MWDhBG0fJTpEqz05halNFtqXuskWBEdIqmpDekamw50U8VTZQb46WbDaP+032G2UC62UCbb/DpaE1voqXwLApq3m7vxnqeBCMjjD+o0uBKzKZ4h+qd7DxlR6fA5QvKEnJNAF3Qzdz2QKSxeCXHZt3QOYMQ0qd1mY2wJH3pYrSymvTMLMkc9a9/kZDfMNEpmlmRg3VOnN4g6SY9M0rOJKjqAy5y6z+EKEp6x+WM7lmRDtH0HOla4jtQv5FAUMXtH/71djGwGrsnIYmrmqbx123hzfDOmZyf0E3QZmy7j3TnUbxpRexd/jD2/MU4c2bhzpyML72UgDkX1WCVQGQIKApU5FtZPilPApF28lsmOkWTL9IxK7JofA4G3Zlfo7zaf6NTA7TZJuLOnDzodfQ6hZIRustktKJ5I+kIRhbrKqmqqR8wqVhmR1JPonqL7D5t51C9C6Ne4bL5pQm5JkDxiZcoPf4XNEXHnuWPEDD3nSQtEi/DYmDpxFymFNqkEq4LCUYERFnSG1I1trYHI0virKIpyrSMmr0VBmM1GbBZIlsZ9dgm4E4fj0kJsTC0m+NN/eebyA6+qaWlzR9VB+P+qJrGX7eHZ0UunFFIboI+Qac5jzNj6zoAjs66kdbCZQm5rhiYTgdTCjNYPjGXzFHcvCxWY+NdQAwqmpLeA7UOXL4gGWYDM4rPVNEoIT/5NRsgmiWa3LGxRNMhmp18OxqgrdTtYl91/0s1MjOSWhI1K7L5WDOnWjykGfVcMqckIddUQn7mbvwmhmAbLQXLODbz6wm5rhhYTrqRsyblUZGfPqr6wySSBCMCol2iaW90tnhCTrdpxpyGzRgDTnyWfOx58we9TpbVOOY+IUS1VFMS3sX3At0O9g3QGr7NFxxwh18xfAIhlXpn/ImrwZDK33ZUA3DxnOJ+u/BGa+quh8ls2YvfnMOe5T8B3egtp08Fer3CjBIbiyfkYjVJvchAYgpGnnjiCSoqKrBYLCxfvpxNmzZFdN6LL76Ioih8+tOfjuVhxRCKNHk1qKpsa290tmRCzyWaNwBoKL0ooqS30V7O2xeLUU+WNbIArKVgOSHFwHhdA2rTEbyBvqf+NQ0cXkliTQWJ6rj670ONNLh8ZFoMrJ5RmIihkX96PeMP/RaAvUsfwmctTsh1Rd8KbGZWTMobdfttDZWog5E//vGP3Hbbbaxbt45t27Yxf/581q5dS319/YDnHT9+nG9/+9ucd9558YxXDIFoSnoP1Dhp84ewWQxML+rS5l1TO1vAR7JEY9ArFEexZDGaFEXY1ChkTMdesASAc5SdHKh19nusLNWkhkQs0XgDIf6xKzwrctm8UsxRbLTYH7O7hlmb7wTgxLRraCpdFfc1Rd9MBh1zx2Uxvzw7qk0yx7qog5FHH32U6667jmuuuYZZs2bx9NNPY7VaefbZZ/s9JxQK8cUvfpEf/OAHTJo0Kd4xiwSLpqR38/Fwo7MlE3K6bWOd2bIHi6eOoCGd5sIVg16nNDttVG6DHYnCTHPEfZ26lfgOkDfS6pYdfJOt1e2PqI/MYN7cX4fDG6Qgw8x5U+PfNVZRg8zZeBsmfyv23LkcnvutuK8p+laSbWHF5LyocsNEWFTBiN/vZ+vWraxefeaTr06nY/Xq1Xz44Yf9nnfvvfdSWFjItddeG9Hj+Hw+HA5Hty8xdCLNFwmGVLafbAXotV9KRxVNY8lKNP3gWf9jpbdIXyxGPdkRLtV0BCMrdPs4VN3U73EyM5J8pxLQcdXlDfL63joAPr2gNCGVZhP3/oKcxq0EjRnsOeunET0/RXTSTHoWjs9mdmnWmO6iGo+ofmqNjY2EQiGKirq3+C4qKqK2trbPc9577z1+/etf88wzz0T8OA8++CBZWVmdX+Xl5dEMU0QhmpLefTUO3P4QWWlGphZmdLuv4HTHEs2aQa+Tm2Ea88lckX5ycmVNx2spwKr4KG/bSXM//1fBkJaQT+UiNoGQSoMz8hnG/vxzTw2eQIhxOWksnRh/74+cug+ZuP8pAPYvvg9Pxvi4rynOUBQYn2flrEl55GVE3tRQ9DakIZzT6eRLX/oSzzzzDPn5kU833nnnndjt9s6vkydPDuUwx7SmtshLejd3qaLpusSS5jxOhuMQqmKgqXjloNcZNwYTV3sqtFkiW6pRFJo7l2p2Solviqq1eyN+HvWnuc3P+gPh3LvPLCxDF2cJqNHbxJyPvo2CxqlJn6Nu/Mfjup7oLt1sYMmEXKYVSfOyRIjq42l+fj56vZ66urput9fV1VFc3Dsz+8iRIxw/fpzLLrus8za1PdXcYDBQWVnJ5Mm9u3SazWbMZokyh5o/qHKgNrIlsEBIZUf7Es3SXlU04VmRlsLlBE0Db21uNuoosMn/rcmgIyfdRLNr8FmppuJzKT3+V1bqdrGuxsG5/eQRtLoDlEqglxSJSFz9x85qgqrG1MIM5pZlxXcxTWX2ptsxextwZU7l4IK74h6fCNPpoCIvnYq89DGb9zYUopoZMZlMLF68mPXr13fepqoq69evZ8WK3kmLM2bMYPfu3ezYsaPz65Of/CSrVq1ix44dsvySRJqmsfu0HV8gsjrEvdUOPIEQOVYjk3ss0RRWt5f0RlBFU5adJk1/2kW6VNNcdA4aCjN1J2msOY6q9f0JXGZGksPuDuCKs7S6xu7hvSONAFyxaFzcz5EJlc+SX/tvQnoLu1c8hmqQIDURsqxGlk3MY1JBhgQiCRb1wv1tt93Gl7/8ZZYsWcKyZct47LHHaGtr45prrgHgqquuoqysjAcffBCLxcKcOXO6nZ+dnQ3Q63YxvI40tEWcK0KXKprFE3K6TR+bvI1kNW6Hjv4iA1AU5JN7F4U2M5U6Bu1LETDnYM+dR3bzThYHt3GyeTkT8npvI+/2h5ufjZX2+qniVKs77mv8bUc1mgYLxmUzpUewH63Mph1M3v0oAJUL7qIta2rc4xvr9HqFKQUZjMuRD1NDJepg5Morr6ShoYF77rmH2tpaFixYwGuvvdaZ1FpVVYVOJy+GqazB6eN4Y/97nfTkD3ZZoqnonlSXX/0WChr23LmDNlEqsJml7r4Lo15HbrqZxggSH5uLzyW7eSfn63bxdo2jz2BE08KzI5JIN3yCIZV6R3yJq8ca29h6ogUFuHxhWVzXMvgdzN14GzotSG35pVRP+lxc1xOQl2FiZkmmvHYNsZhKGm666SZuuummPu/bsGHDgOc+99xzsTykSBCPP8TeAVqL92VPtR1fUCXXamJifvc3wWg2xhuLHVcHU5xpiSgYaSpeyaR9T3CubjdPnW7pd68SCUaGV00CElf/uv0UAGdNyqMsnpJ3TWPmlu+R1nYKd3o5+xffR8QNbUQvRoOOaUUZlGTJ69ZwGNv1lWNMSNXYeaqVYCi6F8+OvWiWVHRfotEHXOTWfQARlPRaTXp5k+xDfoYJvU4Z9A3NkTsXn8FGdtBJWuMufMHpmA29P6m1St7IsIo3cXVftYP9NU70OoVPLSiN61plR/9I0anXUBUDe1b8lJDJFsFZoi/FWRamFdkwGWSWf7jIT3oMOVDriDrRzhcMsfNUe6Oziu5VNHm176FX/bRlVNCW2bsqqivZn6FvBr2OvIzBm1BpOgMtxecAcK6yg0N1rj6Pc0gwMmzsnvgSVzVN4y/tsyIXTCsgP45gPb21kmk7HgDg8Lxv4cidF/O1xjKzUcf88mzmlGVJIDLM5Kc9Rpxu9VDTGv1uortPh5do8jNMTOyRp1BQ3b5EU7Z6wOlgnS7cJln0LdI9eprbe7icr9vFvpq+S7KDIQ2XND8bFqfj7Li6taqFE01uzAYdH5/b97JbJHRBN3M33oo+5KOx5Hyqpl0T17jGqnG5aayYlCetB5JElmnGAIc3QGWE/UR66lyimZDbLYtcUQPkV4fzgwYr6S3KtEiL5AHkZZjR6xVCgyyfNRWfC8A85Qinqk8BfZfG2z2BhG05L/oWDKnUOaMP7juEVI3/234agI/NKiIzLbLtAfoyffv9ZDiO4E0rZO/SH0W0Y3ayKQroFKXzz/AXKO1/6nRd/62gQGdjMZ2ioNPR45z28wa6Zse1upzb9Xh5jUouecUa5QIhld2n7DFta+4LhNh1Kpzs2nOJJrthM8aAA585D3vu/AGvMy5blmgGotcpFGSYqbUP/ObmsxZjt00ly3mISY4t2D1LyOrjTczuDkiy8BCrdXgHDR4H8v6RRuocPjLMBj42a+AqtIEUnfgHZcf+jKbo2Lv8EQKW2FrIZ1uNmAy6yN7MuwQLER/f7TikPFb0IsHIKLfntB2PPxTTuTtP2fGHVAoyzEzI7R5QFHZsjFd6Iej6L3mzWQxkRbgp3FhWlGkZNBgBaC05jyznIVa2L9WsmJTX+xiP7OA71OJZovEHVf6xsxqAS+cWk2aKrWQ0zXmCmVvvAeDYzK/TUrg8pusU2MzML8+O6VwhEkXmpUaxow0umiJoN96fLSfCjc6WVOR0/ySjaZ0lvfWDLNGMy5VZkUjkpZsw6Af/tNhU0r5PjX4n+073XaLt9oUIhGKYChMRcXgDOONIXH27sp4Wd4DcdBOrphfGdA0l5GfuxlsxBNtoKVjKsVlfj+k6VrOe2aUDb+EgxHCQYGSUanL5OBZFY7OevIEQu9vf7Ho2OrO17MXiqSVosNJSdHa/19DrlYiTM8c6nU6JKHGuNX8JAZ2FIqWVQM0eNGkNP+zimRVx+4O8srsGgE/NL405T2Hqrp+Q2bIXvymbPcsfQdNFP8lt0CssKM+Wjr0iJchv4SjkDYTYU+2gn/epiOw81UogpFFkM1PeoxFTx6xIU/FKVH3/b6ClWWmym2UUIgncVL25czp+UWAr1f0s7UgwMjRCqkatI/bE1df31uH2hyjNsvS5xBaJ/Oq3GH8o3Dxy37IfDdr5uD+zS7OwmmSlXqQGCUZGGVXV2HXKTiAY3zT95s5GZ7m9ks0Ku5b0DiCubpJjUG66CWMEvQ1aOpZqdLvYV913lZQEI0MjnsRVuyfAG/vDO55fvrAspo3WzO5aZm26A4AT066msXRVTGOZVJAuJawipUgwMspU1jnjbnzl8YfY07lE072KJs15ggz7QVTFQGPJ+f1eIyfdKOWlUVIUhcII3iAa2/uNLNFVcqS6ts9j7J5Av0s4InbVcXRcfXlXNf6gyqT8dBbEkDCqqEHmbLwNk78VR84cDs/9dkzjKLCZmVQQ32Z8QiSaBCOjSI3dE3cjJoAdp1oJqhrFWZZeJaIdjc5aCpcRNGX1e40yKeeNSSRLNZ6MCTjTyjArQXIbNvWZrBqS5mcJ5/QGsLtjC/TrnV7ePdgIwGcWlcVU2jpx3xPkNG4haEhn94rH0PSDd+7tKd1skIRVkZIkGBklnN4AB2qcCbnW5uPhKpqlE3J6vWgWnF4Pg2yMZzLoIvqEL3rLthoxGwd5WioKraXh2ZFztB0caei7Nbws1SRWPPvQ/H1HNSFNY3ZpJjOKow8Gcuo3MnHfkwDsX3IfnozxUV/DoFeYX54lCasiJclv5SgQbG9sFu/uobRn++9tz0NY0qOKxuhtIrtpGwANZRf1e43S7LSY1sNFeKmmKILZkebicN7IygFaw0swkjghVYuoD0xfqprdfHQsHOB/ZmFZ1Ocbvc3M2fhtFDROT/wsdeM/EfU1FAXmlEnCqkhdEoyMAnurHbhjbGzW0/aTrYRUjdLsvpZo3kbRVBw5c/BZ+95LQ1FgnCSuxqXIFkEwUngWIfRM1NXReqqyz2NiXVIQvdU5vFHvdt2ho+370oocJvTY32lQmsrsTbdj9tbjypxM5cK7YxrDpIKMuDbiE2KoSTAywh1vbKPB6UvY9c4s0fRuK10QQRVNXoYZizG2jpIiLMtqHLQrZ8iYQVPeIgAmOzb1uXus2x/CH2dVlQiLNXH1YJ2T3aft6BT49ILoZ0XGH/wN+bXvEtKb2b3iZ6iG6AP9okwLE/OjDIKEGGYSjIxgLW3+fvMFYuHyBdlfHc476bkXjT7QRm7tezBI11WZFUmMoszBP8Xa2/NGVup2cqCfjRBlqSZ+Ll+Q1hhmmTRN4y/bTgFw3tSCiJbfusps2smUXY8AcHDBXbRlTYt6DBkWA7MkYVWMABKMjFAdHVITWb25o6qVkKYxLieNkqzuQUVu3XvoVT/ujPG0ZU7t8/w0k5689Ogz/EVvhRG8cTW1542crdtLZXVTn8dIMBK/WGdFdp6yc6ShDZNex2Xz+l7W7I/e72TOxm+i04LUjbuE05OujPrxDXqF+eOypfGgGBEkGBmBNE1jz2l7wqfgO5ZolkzI6XVfx8Z4DaWrw4khfSjLTpPdOBMk02LEOshSjSt7Bm3GXNIVH6bqzX32FZFgJD6qqsUUjKiq1pkrctHMQrKtUQTpmsasLd/D2nYKT/o49i25v9/nXH8UBeaWZcW8CZ8Qw02CkRHoUL0rpmnjgTi9Afa3T/X33ItGUQPk12yAAZZodLpwFY1InKKsQWZHFB0t7bMjC/1bqe8jd8jhleZn8ahzxpa4+tGxZk63erCa9Fw8O7p27WVH/0jRqX+iKgZ2n/VTQiZb1I8/uSCDPElYFSOIBCMjTJ3DS1WTO+HX3V7ViqrB+Fxrr7Xt7IYtGP12/OZc7HkL+zy/0GbBFEErcxG5SHIMWsrCXXBX9tMaXpqfxSeWJoKBkMrfd4ZnRS6ZU0x6FJ2I0+0HmbbjAQAOz/0Wjrz5UT9+cZaFCklYFSOMvHuMIG2+YL89JeK1+UT/SzQF1R2Nzi4EXd/TvpK4mngZZsOgb2TNRWejojBLd4KaU8f6PCbRs2hjRVuMiavvHGyg0eUnK83IhTMKIz5PF/Qw98Nb0Yd8NBavpGr6NVE/dobFwMwSSVgVI48EIyNEqH0DvFg36RqIwxPgQG24iqbnEg2a1rlLb0PZmj7PTzcbolsTFxErHmSpJmDOpdE2C4Cixg/6bHwneSOxiaXjqjcQ4pXdNQB8cn4pZkPkORvTt99PhuMwPkshe5c9DEp0L89Gg04SVsWIJcHICLG/xkHbEE23b6tqQdOgIs/aaydPW+t+0tzVhPRpNBeu6PN8mRUZOpGU+DrLwiW+K7QdHGts63V/vBsnjkWqqlETQ8fVN/bV4fQGKbSZOWdKXsTnFVW9TNmxP6GhsOesnxCw9O7zMxBJWBUjnQQjI8DJZnfMragjseVECwBL+mp0dvoNaC8jVQ29P6XrdQolgyVaiphZTQZslkGWatr7jZyn283+6pZe90vzs+jVO30EovyZOb0BXt8X3kX50wvKMOgie3lNc1Uxc0u4s+qxWV+npfCsqMc7pTCDXCmrFyOYBCMpzu4OcKg+MRvg9Xl9T4DKur4bnQGdSzT9VdEUZ1lk460hNthSjSN3Ph5dBjmKi+Cp7X0e0+rxD9HoRqdYlmhe3V2LN6AyPtfa53OpL0rIz5wPb8UQbKMlfwnHZt0Y9eMWZ1mibzMvRIqRd5EU5g+q7DrdijqEH2q3nggv0UzMT++1d0WaqwqbvRJV0dNYckGf58sSzdAbrKpG0xloKAgvoU1xbMTTxz5FslQTObc/SEtbdMFbk8vH25X10L4Zni7CviBTdv+ErJY9+E3Z7DnrUTRddBvZ2SRhVYwSEoykKE3T2H3aji8wtNPrW9qraJb2OSsSrqJpLVhK0Jzd6/4sqxGbxTik4xNgMerJtg78c3aVd5T47uyc6epKklgjF0s570s7qwmqGtOLbMyOsP16fvXbTDj4HAD7lj2EzxpdPxKjQcf8cklYFaODBCMp6khDW9SfzqLV4vZzqC68t83i8X2V9A68MZ7MigyfwWZHOlrDL1AOc+zkqV73OzxBaX4WAVXVqI4yP6u61cMHR8Pt+K9YVBZRF2Kzu5ZZm74LQNXUL9NYemFUj9mRsCqbUorRQoKRFFTv9HK8j6qIRHL5gjzx9mE0YHJBeq9ujUZfM9mNW8PjKe0djBgNuoi2uheJUZhpHrAjuM9aQmPaRPSKRmbt+73uD6kaTml+NqgGV/SJq/+34zSaBgvHZzOpIGPwE9QQcz76FiZ/K46c2Rya952oxzm10CYJq2JUkWAkxbj9wT47aSaS0xvgkX9VcrzJTYbZwH8un9DrmPzqDSiaiiNnNr700l73l2ZZ0Mn08LAxG/SD9nLp2MV3gXcrTa7ereHt0vxsUNEmrh5tcLG9qhVFgcsXlEV0zqR9T5DTsJmgIZ3dZ/0UTR9dUFGcZWF8njWqc4RIdRKMpJCOxmax7IURqVa3nx+/XsnJFg+ZFgPf+dh0ynN7v7B1lPQ29DErAlAmSzTDbrCeI/aO1vD6Xeyrtve+X/JGBuTxh2h2Rb40qmkaf9kWbvt+9qS8iPZmyqn/iIn7nwRg/+J78dgqohpjZpqRWZKwKkYhCUZSyIFaBy7v0E2lN7eFA5Fqu5ccq5Hb187oM6jQBT3k1YWn+hvKLup1f26GCaspuqx/Eb9Cm2XApZrWgqX4FTMlSjPOk7t73S/ByMBOt0a359PeageVdU4MOoVPzu89e9iT0dvM7I++haKpVFdcQd2Ey6J6PJNBx7xxWTIjKUYlCUZSxOlWDzWtQ9fYrMHp4+HXD1Dn9JGXbuL2tTP67V+RV/ce+pAXd/o4XFnTe90viavJYTLoyBkgT0DVm6nJWQJASeMHqD0SVj3+EL5g77JfEZ7lqI7i+adqGn/dHp4VWTW9cPAdcjWN2Zu/i8VTjytzMgcW3R3V+HQ6SVgVo5sEIynA4Q1QWTt0eSK1Di8Pv36ARpefQpuZ29dO79X2vasze9GspudHcbNRR4FsTZ40xYNU1bjbS3zPUrdT1dz7k77MjvStwemLqkvtluMtVDW7sRh1XDp38JLc8Qd/Q37NO4R0Jnaf9RiqIbqcj6mFtgEDUSFGOglGkiwQUtl9yj5kjc1Ot3r48euVtLgDFGdZuH3t9AE/xSlqkPzqt6GfjfHKstMiKl0UQ6PAZmagLuMtJeEk1mW6Axw+Vdfrfkli7dupKBJXg6rK33aEZ0XWzioetNdOZvMupux+BICDC+6iLbv3bONASrItfeZ1CTGaSDCSZHtO2/vsmJkIVc1ufvx6JXZPgHE5adz+semDVmRkN27F5G/Fb87Bnrew232KQkRJemLoGPU68tL7Dybdtom0GIsxK0GMpz7sdb/MjHTn9gc51hhdT5/3DzdR7/RhsxhYM6towGP1fidzPvwmOjVA3biLOT35P6IaX2aakZnFkrAqRj8JRpLoaIOLpiiy96O6dqOLn/yrEpcvSEWelW9/bDqZaYN3S+1YomksubBXa+oCm1nWrFPAgA3QFIWGonADtCmOjb1yRBzeAKo6tpufdQQgHx1t4oPDTRypdxFpPzhfMMRLO6sB+PjckoGfD5rGzK13Y207iSd9HPuX3N9r2XMgkrAqxhIpiUiSJpevz+3eE+FQnZOfvXUIb0BlckE6t1w0NbLqF00bcGO8cTkyVZwKCmxm9DqFUD9BhXv8BXDqT5yr7OIfdS7mlGV13qeq4PQFyYogMB1N3P4gdQ4fdQ5vXBVrbx2ox+4JkJ9h4vxpBQMeW3rsTxSffBVVMbD7rEcJmiKf4dDpYN44SVgVY4cEI0ngDYTYU+2I+NNYNPbXOPj524fxB1WmF9n4xoVTIn5By7AfIM19mpA+jeaic7rdZzXppeNjitDrFPIzzNQ5+q7+aClaQRA9k3U11FUdYE7Z8m73292BMRGMJCoA6dDmC/LPPbUAfGp+GcYBdqtOtx9i+vb7ATgy95s48hZE9VjTimyDLqkKMZpIMDLM1PbGZtG2nI7EntN2nthwmEBIY3ZJJl9fNRmzIfJPVh2zIk3F56Iaui8FyKxIainK7D8YCRkzqLbNY7xzO7k17wE9gpFRnDfS5gtS5/BS7/QlvGfPa3trcftDlGWnsXxibr/H6YIe5n54K/qQl8bi8zgx/dqoHqc0O02eb2LMkWBkmFXWOYdkO/ftVS08/e5RQqrG/HFZXH/+5AE/ufWl8HTfG+PpdQol2bIPTSrJzzCj1yuE+unW6yxbCQe2M9+3hVa3v9un7NEWjHQEIHUOH21DtP9Oq9vP+v31AFy+sGzAPI5pOx4gw3EIn6WAvcseBiXy52GW1ciMYltCxizESCIJrMOoxu6JaXvywWw61szT74QDkcUTcrghhkDE0nYKW+t+VEVPY8kF3e4rzDRHfT0xtHQ6ZcB+L67y8P/h2bq9VJ5u7nafNxDCGxjZzc9cviBHG1x8eKSJD480cbShbcgCEYB/7KrBHwrnYM0fl9XvcUVVrzDu6P+iobBn+Y8JWPIifgyzUcfcMklYFWOTzIwME6c3wIEaZ8Kv+/6RRp774DiaBmdNyuWasyeij+HFrOD0egBa85cQMOd0u0+mjFNTUaaF2n62u3dmz8ShzyEz1EKo6iOY+qlu9zs8gRGXHOnqWIIZwhmQvtQ5vPz7UAMAVywa12+fnTRXFTO3hjurHpt5Ay1FZ0f8GDodzCvLHnH/J0IkigQjw6CjsVl/1Q+xeudgA7/beAKA86bk86WzJsT8qaqgnyUam8UwJpIdR6K8dBMGvdL3xoqKjuq8s8is/ydlTe+jaZ/s9iba6glQOEg311Tg9Aaod4aTUN2+5Mzm/G3HaVQt3I59WlHfSyhKyM+cjd/EEHDRkr+YY7NviuoxphdnkmWV55kYuyQYGQb7qh24E9zY7M39dby4+SQAF04v5D+WlaOLsTOq0ddCTuNm6GNjvHHS+TFl6XQKhTYL1f10D/VNuBDq/8lZ6na2tHq7bYqYynkjTm+AOoePemfyApAOJ5ra2Hy8BYDPLCzr97gpux8lq3k3AVMWe856tFePnoGU5aRRJs0ExRgnwcgQO97YRoPTl9Br/nNPTefW5WtnF/HZAaaOI5Ff8zaKpuLMnok3fVzn7Qa9MuheKCK5ijLN/QYjrSXnAjBbd4K/VB2lLGd2533O9uZnqZKf0BmAOLwJD9zj0bEZ3vKJuf22ZM+rfpsJB58FYO/Sh/BZSyK+frbVyPR+ZluEGEskGBlCLW1+jjS4EnY9TdP4x66azg6Ql80r4ZPzS+PeK6a/RmclWWkx5Z+I4ZObbsJk0PW5yVvAksdJy3TKvZVknHoX5p8JRlQVnN5gUpcGUjUA6XCg1sHeagd6ReFTC0r7PMbsrmX25jsAqJp6FY09ZhYHYjbqmCsdVoUAqaYZOt5AiN2n7QlrbKZpGn/ZdrozEPnMwjI+taAs7kBEF/SQV/se9JEvMi5Hpo5TnaIoFGb2X1XTXNLeGt65iUCoe8CSjKUahzfA4XonHxxu5KOjzRxvbEvJQETTNP7aPvu4clo+hbY+ZgjVELM/+g4mXwuO7Fkcmnd7xNcPd1jNjqoPkBCjmQQjQ0DTNPactke1Jflg13tx80le2xvu/njlknIunRv5VPBAcus+QB/y4kkfhytrRuftOelG0s0ycTYSFPX1RtnOP2EVAOcqOzlab+92X6tnaPZF6qlrALLpaDPHG90pGYB0teNkK0cb2zAZdHxiXt+zIhP3P0Vuw0cEDensXvEYmj7yjqkzijMlMVyILuTdZggcqnfRmqCt2lVN4382nuDdQ40AfHH5eFZNL0zItena6Kx0dbdNvKScd+TISTdhNurwBXoHv478BbgVK7m4cB3bCiVnZr+GcmbE7gnQ4Aw3IhuqXamHiqpqnbkia2YW9Rk0ZNdvYtK+XwBwYPH38dgqIr7+uNw02f1aiB4kGEmwOoeXqiZ3Qq4VUjWe++A4Hx5tQlHg6hUVnDMlPyHXBlDUIPnV4f4iXatoTAbdgA21ROopyrT0+Xun6YycyFrGzNYN5NW9B5wJRnwBFW8glLDeFnZPgPr2VuwjLQDp6oOjTdTYvaSb9KydXdTrfqOvmTkffQtFU6mu+Ay1Ez7V53X6kpMuCatC9EWCkQRq8wXZV+NIyLWCqsqv/n2MLSda0CnwlXMnsWyA/TCiZWk7TcX+pzH5W/GbsmnNX9x5X2l2miTVjTBFtr6DEQB3+fnQuoH5vi2c9gbJsJx52tvjbH5mdweob58BGeldXWnvCfTSjnBe1qVzS3rvdq1pzNp0JxZPHW22iVQuvDvia1uMeuaUZcWd5yXEaCTBSIKE2jfA62+vkGgEQir/792j7DjZil6n8LWVk1g0PieCMweX0XqACQd+RdHJV9Bp4TePU1O+2NkXQVEkcXUkyrIaSTPp+5yR8Iy/AHb/gIXKYf596jTzpkzovM/uCVAUZfm23R2gzhnuhDoaApCu3q6sp9ntJ8dq7HM5tPzQ8xTUvE1IZ2L3ip8RMqZHdF2dDuaOy5KEVSH6IcFIguyvcSSkRbU/qPLkhsPsqXZg0Cl8/YLJzBuXHd9FNY3shs1UHPgl+bXvdt7cVHQ2J2ZcR3PhmbbVeRlmaUk9QhVlmjne2Ht2xJteRo1xPCWBKvTH34UpX+q8L9LcptEcgHTw+EO8ujucJH7Z/FJMhu75/bbm3Uzd9WMADi74b1zZM/q8Tl9mlkjCqhADkWAkAU42u/vdIyQavkCIn799mAO1Tkx6HTetmsKs0szYL6ipFJx+k4oDz5DVvDN8k6KjbtxaTky/DmfunF6nyKzIyFWUaekzGAGoyT+Hkpoqxrd8gKb9Z+dSgcvXf/OzVre/sxV7X8mxo82/9tXi8gUpyjRzzuTuuVn6gIu5G7+JTg1QN24tpyd/PuLrludaKcmS55UQA5FgJE52d4BD9fFvgOfxh/jZ+kMcbnBhNui45aKp/e6DMRgl5KfkxN+ZUPkr0p3HAAjpTNRMvIIT06/FkzG+z/PSTHry0iMvTxSpxWYxYjXr+2yhHph4IdT8gWXqDt52eCluf3NU1XDpbbY1/P/e6vZ3tmIfCwFIB4cnwL/21QFw+cKy7s3+NI2ZW+7G6qrCYy1j/5L7u1WeDSQn3ci0ooyhGrYQo4YEI3HwBUPsOt2KGudrtssX5LE3D3K8yU2aUc+tq6cyuSD6FzB9wMW4Iy9Sfug5LJ56AALGTE5N+QInp16F3zJwJU5Zdpok141wRZkWjjW09bq9rXg5PkyUKU00HttF8YLlnfdVt3rHZADS1Su7a/AFVSbkWVncIz+r9NifKT75CqqiZ89ZjxI0ZUV0TYtRz9yybHlOCRGBmJqePfHEE1RUVGCxWFi+fDmbNm3q99hnnnmG8847j5ycHHJycli9evWAx48U4cZmjrhfvJ3eAI/8q5LjTW4yzAa+/bFpUQciJk8Dk3f9hHNfPp+pux7G4qnHm1bEwfnf5b1PbODI3NsGDUR0OqT3wSjQ315CqsHC0fQFAGSefqfbfdWtHk42u8dsINLo8rHhYAMAVyzsvs9Tuv0w07ffB8CROd/Enr8womvqdQrzyrN65Z0IIfoW9TPlj3/8I7fddhvr1q1j27ZtzJ8/n7Vr11JfX9/n8Rs2bODzn/88b7/9Nh9++CHl5eV87GMf4/Tp04kYf9IcaXDR0hZfB8tWt58fv17JyRYPmRYD3/nYdCbkRZadD5DmPM6MLXdzziurmHjglxgDTtpsk9i79EHev3Q9VdOvJWSMLLAptFnkhXMUSDcbupXudmVvbw0/zbWJkJqgfQpGOLc/yJ+3niKkaswssXXL0dIFvczdeCv6kJemonM5MeMrEV93ZkkmmRZJWBUiUoqmRbd7yvLly1m6dCm/+EW4+6CqqpSXl/ONb3yDO+64Y9DzQ6EQOTk5/OIXv+Cqq66K6DEdDgdZWVnY7XYyM+NI6EyQeqeXXSftERzZv+Y2P4/8q5I6p48cq5FvrZlOcVZkJZa25t1UHHiGwlOvoxD+72vNW8iJGdfRUHohKNEHFUsqcjrzBsTIdqyxjSP1vTdoTGs9zDn/uhSfZuR3KzdQUVKQlPElU5svyKF6F5W1TirrnJxsdtPxAnjXpTOZmH/mw8CMLfcw7uiL+Cz5fPSxlwadXewwPs8ac76XEKNNpO/fUeWM+P1+tm7dyp133tl5m06nY/Xq1Xz44YcRXcPtdhMIBMjNTVwDr+Hk9gfZVx1fY7MGp49H3qik0eUnL93Etz82nQLbIB1PNY3cuvepOPBLcus3nrlWyQWcmHEdrflLIk6q6yndbJBAZBQpzrT0GYx4sibToCukQK3Hf/Q9KLk8KeMbTi5vkEP14cCjstbJqRYPPT99FdnMXDSzqFsgUnjyn4w7+iIaCnuX/zjiQCQn3cTUQklYFSJaUQUjjY2NhEIhioq6t0guKiriwIEDEV3ju9/9LqWlpaxevbrfY3w+Hz6fr/PfDkdiuprGq6OxWTCOxma1Di+P/KuSFneAQpuZb62ZRt4ArdcVNUjhqdeZcOAZMlv3AaAqBurGf5zj079CW/b0mMfSQcp5R5c0k57MNCOOnnvPKAoncs6ioOkliurfA0ZfMOL0BjhY56KyzsnBunDw0VNxpoXpxTamFWUwvcjWKxC3uE4ya8tdAByf+TWai86J6LHTTHrmSodVIWIyrNU0Dz30EC+++CIbNmzAYul/SeLBBx/kBz/4wXAOLSIHah24vLE3Njvd6uGRf1Xi8AYpzrLw7TXT+p2R0AW9lB7/C+Mrn8XadhKAkD6N05M+x4lp1+BL73sn0Wjp9QolES4PiZGjKNPcOxgBvBNWQdNLLPBtYa8/2Lvd+Qjj8AQ4WNc+81HnpLq1d7+fkiwL04tsTCuyMb3YNmDzMSXkZ+7Gb2IIuGjNW8TR2TdHNA69TmHeOElYFSJWUb0S5efno9frqaur63Z7XV0dxcXFA577k5/8hIceeog333yTefPmDXjsnXfeyW233db5b4fDQXl5eTRDTbhTLW5q+nihi1RVs5tH3ziIyxdkXE4at62eRmYfL4oGv53ywy9Qfui3mHzNAPhN2ZycehWnpnyRgDkxbeEBDHqFGcWZGPTyAjraFGVaOFTXe6kmMH4lwW06JutqePN4JZOnzU7K+GJl9wSorHV2BiA1fTQbLM0OBx/Ti2xMLRo4+Ohp8p6fktW8i4Api91nPdq5TcJgZpVmYpOEVSFiFlUwYjKZWLx4MevXr+fTn/40tCewrl+/nptuuqnf8x5++GEeeOABXn/9dZYsWTLo45jNZszm1Nk11uENf/qK1dFGF4+9eQi3P0RFnpVbV08jw9z9R2921zL+4G8oO/pHDMFwF01P+jhOTLuG6omfRTUkdimlwGZmerFNWr+PUhajnmyrsVe795DJxlHzLKb59mCp2gApHoy0uv3tSy7hpNNaR+/gY1xOWnjWoyi89BJrUJBX8w4Vlb8GYN/SH0Y8+zghzxr1/j5CiO6inqO97bbb+PKXv8ySJUtYtmwZjz32GG1tbVxzzTUAXHXVVZSVlfHggw8C8KMf/Yh77rmH3//+91RUVFBbG977ISMjg4yM1E/0CoRUdp+yx9zY7FCdk5+9dQhvQGVyQTq3XDS129R4uv0wEyp/RXHVP9Cp4TcOZ9Z0js/4KvXll0T8ySxSZqOO6cU2Cm3y4jnaFWVa+tx7pq7wXKad3MP4lo2EuDEpY+tPc5s/POvRPvtR5/R1u19pDz6mF7fPfBTa+i1ljobJU8fsTbcDcHLKf9JQtiai83IzTEyRhFUh4hb1s/jKK6+koaGBe+65h9raWhYsWMBrr73WmdRaVVWFTndm2v+pp57C7/fz2c9+ttt11q1bx/e///1EfA9DJtzYzN7nTqiR2F/j4OdvH8YfVJleZOMbF07pnInIatxGxYFnKKhe33l8c8FyTsy4jqbi82KujBlIWU4aUwozMMqyzJhQmGnmYJ2TnsX76qQL4eTTLFF38U+Hi9zM5L2ZNrl8nQmnlXVOGvoIPspzreFll2IbUwszSDcnOM9FDTHno+9g8rXgzJ7Jofnfjeg0qySsCpEwUfcZSYZk9Rk52uDiaB+ttSOx57SdJzYcJhDSmF2SyddXTcasV8iveYcJB35JTuNWADQUGsrWcHzGdTjy5if4OwizmvXMKsmU8t0xaOuJlt7N+TSVJX9eRrbm4NdTn2DCwshmARKh0eXrLLM9WOek0dV9bIoCE3KtncmmUwszhibJVlMx+lqwuGsornqJCQefI2iwsmnN/+G2TRz0dL1OYenE3F7LrUKI7oakz8hY0uTycawxtkBke1ULT797lJCqMX9cFjecN57yUy8x4cCvyHAcAkDVGamZ8GlOTL8Wd+akBI8+TKeD8bnpTMpP73NXVjH6FWdZegcjio5DGctY6nyT7Op3YYiCEU3TaHT5uwUfTT3GolNgQl56Z5ntlEQEH5qG0d+CxV2L2V2DxdP+p7v2zN89tZ3Loh0OLPp+RIEIwOzSTAlEhEggeTb1wRsIsafa0Wt6OxKbjjXzq/eOompwTrmF75W8x8TXn8PirgEgaEjn1OQvUDXtKvxpRYNeL1ZZViMzim2S4T/GFWSYqdTRK+fJUbYSDrzJDNcmTqtaQoJVTdOod/o6K10O1rpodncPPvSKwoT2DqUzisPBR1RJ1JqGwW/H4qnB7K5tDzDa/+6pxeKuweypRR/yDX4pFPyWArzWIurKP05txacjGkJFfjqFkrAqREJJMNKD2t7YLBCMPmP1/SONPPfBcXI0B3flv8snW1/B2BBuG++zFFA19cucmvx5QqahaxWt1ytMzs+gPFd24BVgMujIsZpo6rkcMuVCOHAPs5VjbK89RXFp9KXzmqZR5/SdKbWtddLao7eJXlGoyLd2JpxOLhgg+NA0DAEn5vagwuI+E2R0zGZY3LXoQ70bmfXFZ8nHl1aM11qMN60En7UYr7UEb1oxPmsxPkshmj66pcu8DBOTCyLfP0oIERkJRnqorHP22SxqMO8cbODtjzbzff0r/IfxXcyu8Cczd8YEjk//CrUVn0bVD225cl6GiZklmVKuK7opzrL0CkZUayFH9JOZHDqCcvQtKP3yoNfRNI1ah7c9+Agnndp7Bh86hUn56Z1NxiYXpGNu/33UB1xY3Mcwe2rag43azsDD3B5oGIKRLY36zTntAUZJe7BR3P73EnxpRXjTiqMONAZjNemZIwmrQgwJCUa6qLF7ON1H++jB7N32Pksqn2GdaSMGRQUN7DlzODHjq9SXrQHd0AYHRoOO6UW2iDfaE2NLQYYZXR9LNSdzz2ZywxFKGt5HpXcwomkaNXZv56ZyB+ucOHp0IDboFCYVpDO3wMDCLDdTLa1k+I+Hg43TNZgPn1lKMQR6N2Hri9+UHZ7F6JjNSAvPaHS9bagD+570eoV55dlSiSbEEJFgpJ3TG+BATRSNzTSNnIaPSN/8C1a3bYL2eKOx6FxOzPwqLQXLh6Q8t6fiLAvTimzShlr0y6DXkZdu7lU265+4Chp+x3z/Njb6AxiNBmpavZ1ltgfrnDi9QSz4KFGamaM0UWZoZna6i6kWO2X6ZnKDDaS5azEejmz/qIAx88zshbWPmY204oQ3+EsESVgVYmjJs6tLY7OQGkHGqhqisPpNJuz/JVktuwEIaQo7MlcRWP4NXLnD09EyzaRnRrFtwE32hOhQnGXpFYzoypfh2pRGnuLgwGtP0eiBnGA9M5VmLlSaKFWaKTE3k6P0CNJ97V89BI0ZeNPCwYWvYzaj/c+O20LGkZdvUZGfLk0ChRhiEowA+6oduAdpbKaE/JQc/z8mHHyWdOcxALyakf8NXUDV9P/irMWLh2WsigLjc61MKshAL+W6IkL5GWb0OqV7wK03ccCykCXeD/iW9xfhDmP9FF8FDdYugcWZmY2uAUfIOPo6kebbzNJhVYhhMOaDkRNNbb0+MXal9zsZd+QPjD/0PGZvAwBtOhu/9l/E88G1rF4yhzWzhq5Et6sMi4FZpZlkSrmuiJJep5CfYaaux94uznnX0Lz1IEF9Gn5rCWSV4bf2rDwpIWi0DcuyYyqxmvXMLh2+JotCjGVjPhjp2YSpg8lTz/iDzzPu6B86E++8acX8xfwpHqhdjhsLX1w+nlXTC4d8jDodTMrPYEKeVTL5RcyKsnoHI6GK89lW8UHSxpSq9HqF+eMkYVWI4TLmg5GerM5jTDjwK0pO/K2zQ6MrcwrHpn2FH52ew3vHHCjA1WdXcO6U/CEfT066iZkltqFpiS3GlPx0M3q9QiiU8jtAJN3s0szE74EjhOiXPNvaZTbtDG9cd/oNFMIv1q35izg+46vUFq3kV++dYMuJFnQKXHvuRJZPzBvS8Rj0ClOLbJRlp15lgRiZdDqFggwztXZvBEePXZMKJGFViOE2toMRTSPz1AYqtj9FbsNHnTc3lF7I8RnXYc9fTCCk8v/ePcqOk63odQpfWzmJReNzhnRYRZkWphVnYDZI8zKRWMVZFglGBlBgMzOpQBJWhRhuYzcYCQXg12uYUr0dAFUxUDvhMk5M/wptWVMB8AdVntxwmD3VDgw6ha9fMJl547KHbEhmo47pxTb5VCaGTK7VhEGvEJSlml4kYVWI5Bm7wYjeCHlTCNVXcmrS56iadg0+a0nn3d5AiF+8fZgDtU5Meh03rZrCrCF8oRqXm8aUggwMkjAnhpBOp1Bos1DdGn2n4dHMoFdYUJ4tzz8hkmTsBiMAa+5lz7y7aAhYu93s9gd5fP1hDje4MBt03HLRVKYVDc3mdlaznlklmWRbE7uPhhD9Kc6SYKQrRYHZpVmSJC5EEo3tZ19mKaHWFgicKe91+YI89uZBjje5STPquXX1VCYPwRqyTgcT8tKZmJeekO3bhYhUjtWIyaDDH8PO1COZooQ7F2eYDVhNBjLMBtLNetJNBnkOCpFkYzsY6cHpDfDoGwc52eIhw2zgm6unMiEv8e2rs6xGZpbIXhciORRFoTDTzKnm0Tk7otNBmrFLsGE2kG42YDXqJegQIkXJu2G7VrefR984SLXdS6bFwLfWTKcsJ7FltXq9wpSCDMblpEnzMpFUxZmWER+M6HVK50xHepdZDqtJL88vIUYYCUaA5jY/j/yrkjqnjxyrkW+tmU5xVmIrWvJtZmYU27AYpVxXJF+21YTZqMMXSP2lGr1OwWoKz3BkmA1YzeEAJM0oQYcQo8WYD0bqHF4efv0AjS4/eekmvv2x6RTYErcTrskQLtctypRyXZFaijItVDW5kz2MTnq90p7P0WW2w2QgzSQBvBCj3ZgORo42uLjr//bQ3Oan0GbmW2umkZeRuECkJNvCtCKb7G8hUlKyghFDZ9DRPa9DZg2FGLvGbDDS5gvy+Wc20tzmpzjLwrfXTEtYea3VpGdGSSa56VKuK1JXVpqRNJMejz80JNc3GnRkmPVdgo7wrIcEHUKInsZsMJJuDiepPrnhMDdfOJXMNGPc11QUmJBnZWJ+BnrJ2hcjQFGmheONbXFdw2TQncnn6LLEYjLIjKAQIjJjNhgB+NzSciryrTg8wbivZbMYmFmaSaYl/qBGiOFSlGmOOBgxG88EHeF8jvDyiixDCiHiNaaDESDu9s96ncKkgnTG51ols1+MODaLEatZj9t3ZqnGYtST3l6xYjUbyDCF8zqkVboQYqiM+WAkHrkZJmYWZ0q2vxjRphXZCITUztwOWWIUQgw3CUZiYNArTCuyUZqd2KZoQiRDfgIryIQQIhYSjESpOMvC1KIMzAaZDRFCCCESQYKRCFmMeqYX2xLaEE0IIYQQEowMSlFgXI6VyQXpksAnhBBCDAEJRgaQbjYwqySTLKuU6wohhBBDRYKRPuh0UJGXTkVeumw5LoQQQgwxCUZ6yLYamVmSSbpZfjRCCCHEcJB33HZ6vcLUwgzG5ViTPRQhhBBiTJFgBCiwmZlebJMNvIQQQogkGPPByLQiGxmyJCOEEEIkzZivVZVARAghhEiuMR+MCCGEECK5JBgRQgghRFJJMCKEEEKIpJJgRAghhBBJJcGIEEIIIZJKghEhhBBCJJUEI0IIIYT4/9u7u5Cm4z2O45/p2pSappbm8qkHy8oc5lLMIkgpRKRuTMJgZV0UkzQJoqt15bwJekAsKywIsQi0BzAzMyNIfGKgBZYlJVlKUDp3kbH9zsUhOdLpnA7U/1vn/3nBQCds799P2b7+95+K4jBCREREojiMEBERkSgOI0RERCSKwwgRERGJ4jBCREREojiMEBERkSgOI0RERCTKKB3wM5RSAICpqSnpFCIiIvpJ3563vz2P/8hfMYx4vV4AQHx8vHQKERER/Y+8Xi/Cw8N/+HWD+m/jyh8gEAhgbGwMFosFBoPhl93u1NQU4uPjMTo6irCwsF92u38Tve+B3tcP7gHXr/P1g3vwW9evlILX64XVakVQ0I/PDPkrjowEBQUhLi7ut91+WFiYLn8A/5Xe90Dv6wf3gOvX+frBPfht6/9PR0S+4QmsREREJIrDCBEREYnS9TBiNpvhcrlgNpulU8TofQ/0vn5wD7h+na8f3IM/Yv1/xQmsRERE9P9L10dGiIiISB6HESIiIhLFYYSIiIhEcRghIiIiUboeRmpqapCUlISQkBBkZWWhu7tbOkkzjx8/RmFhIaxWKwwGA5qbm6WTNOV2u7Fx40ZYLBZER0dj165dGBoaks7STG1tLdLS0mb/yFF2djZaWlqks8RUV1fDYDCgoqJCOkUzJ0+ehMFgmHNJSUmRztLUu3fvsHfvXkRFRSE0NBTr169Hb2+vdJZmkpKSvvsZMBgMcDqdmrfodhi5fv06Kisr4XK50N/fD5vNhh07dmBiYkI6TRM+nw82mw01NTXSKSI6OzvhdDrR1dWFtrY2fP36Fdu3b4fP55NO00RcXByqq6vR19eH3t5ebNu2DTt37sSzZ8+k0zTX09ODCxcuIC0tTTpFc+vWrcP79+9nL0+ePJFO0synT5+Qk5ODefPmoaWlBc+fP8epU6cQEREhnaaZnp6eOd//trY2AEBRUZH2MUqnMjMzldPpnP3c7/crq9Wq3G63aJcEAKqpqUk6Q9TExIQCoDo7O6VTxERERKhLly5JZ2jK6/Wq5ORk1dbWprZu3arKy8ulkzTjcrmUzWaTzhBz/PhxtXnzZumMP0p5eblasWKFCgQCmt+3Lo+MzMzMoK+vD3l5ebPXBQUFIS8vD0+fPhVtIxmTk5MAgMjISOkUzfn9fjQ2NsLn8yE7O1s6R1NOpxMFBQVzHgv05OXLl7BarVi+fDlKSkrw9u1b6STN3L59G3a7HUVFRYiOjkZ6ejouXrwonSVmZmYG165dQ2lp6S/9h7Q/S5fDyMePH+H3+xETEzPn+piYGHz48EGsi2QEAgFUVFQgJycHqamp0jmaGRgYwIIFC2A2m3Ho0CE0NTVh7dq10lmaaWxsRH9/P9xut3SKiKysLFy5cgX37t1DbW0tRkZGsGXLFni9Xuk0Tbx+/Rq1tbVITk5Ga2srDh8+jCNHjuDq1avSaSKam5vx+fNn7Nu3T+T+/4r/2kv0OzmdTgwODurq9XIAWL16NTweDyYnJ3Hz5k04HA50dnbqYiAZHR1FeXk52traEBISIp0jIj8/f/bjtLQ0ZGVlITExETdu3MCBAwdE27QQCARgt9tRVVUFAEhPT8fg4CDOnz8Ph8Mhnae5y5cvIz8/H1arVeT+dXlkZNGiRQgODsb4+Pic68fHx7FkyRKxLtJeWVkZ7t69i46ODsTFxUnnaMpkMmHlypXIyMiA2+2GzWbDmTNnpLM00dfXh4mJCWzYsAFGoxFGoxGdnZ04e/YsjEYj/H6/dKLmFi5ciFWrVmF4eFg6RROxsbHfDd5r1qzR1UtV37x58wYPHjzAwYMHxRp0OYyYTCZkZGSgvb199rpAIID29nbdvWauV0oplJWVoampCQ8fPsSyZcukk8QFAgF8+fJFOkMTubm5GBgYgMfjmb3Y7XaUlJTA4/EgODhYOlFz09PTePXqFWJjY6VTNJGTk/Pd2/lfvHiBxMREsSYp9fX1iI6ORkFBgViDbl+mqayshMPhgN1uR2ZmJk6fPg2fz4f9+/dLp2lienp6zm9AIyMj8Hg8iIyMREJCgmibFpxOJxoaGnDr1i1YLJbZc4XCw8MRGhoqnffbnThxAvn5+UhISIDX60VDQwMePXqE1tZW6TRNWCyW784Pmj9/PqKionRz3tCxY8dQWFiIxMREjI2NweVyITg4GHv27JFO08TRo0exadMmVFVVYffu3eju7kZdXR3q6uqk0zQVCARQX18Ph8MBo1FwJND8/Tt/kHPnzqmEhARlMplUZmam6urqkk7STEdHhwLw3cXhcEinaeLfrR2Aqq+vl07TRGlpqUpMTFQmk0ktXrxY5ebmqvv370tnidLbW3uLi4tVbGysMplMaunSpaq4uFgNDw9LZ2nqzp07KjU1VZnNZpWSkqLq6uqkkzTX2tqqAKihoSHRDoP65wMzERERkQhdnjNCREREfw4OI0RERCSKwwgRERGJ4jBCREREojiMEBERkSgOI0RERCSKwwgRERGJ4jBCREREojiMEBERkSgOI0RERCSKwwgRERGJ4jBCREREov4B+JJRJp0eB2EAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#use priority sampling to compute inclusion probabilities\n",
    "def ps_incl_probs():\n",
    "    uniforms = [np.random.uniform() for i in range(N)]\n",
    "    k_is = [uniforms[i] / p[i] for i in range(N)]\n",
    "    S = sorted(range(N), key=lambda i: k_is[i])[:k+1]\n",
    "    τ = k_is[S[k]]\n",
    "    S = S[:k]\n",
    "    pi_i = [min(1, p[i] * τ) for i in range(N)]\n",
    "    return pi_i\n",
    "\n",
    "# compute 100 times and compare against inclusion_probabilities\n",
    "reps = 10000\n",
    "pi_i = np.zeros((N,reps))\n",
    "for i in range(reps):\n",
    "    pi_i[:,i] = ps_incl_probs()\n",
    "#plot both inclusion probabilities and variance of ps\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(np.mean(pi_i, axis=1), label='PS Inclusion probabilities')\n",
    "plt.plot(np.asarray(list(inclusion_probabilities.values())), label='Sampling design')\n",
    "plt.fill_between(range(N), np.mean(pi_i, axis=1) - np.std(pi_i, axis=1), np.mean(pi_i, axis=1) + np.std(pi_i, axis=1), alpha=0.3)\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.4352892128123305"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#compute mean sum\n",
    "np.mean(np.sum(pi_i, axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.0"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(np.asarray(list(inclusion_probabilities.values())))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}