{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Python implementation of the SW1Pers algorithm described in [0]\n",
    "\n",
    "References\n",
    " - [0] SUPPLEMENTARY INFORMATION SW1PerS: Sliding Windows and 1-Persistence Scoring; DiscoveringPeriodicity in Gene Expression Time Series Data\n",
    "\n",
    "Conventions:\n",
    " - point clouds/matrices are stored in ndarrays, where the last index indexes in a pomt coordinate tuple\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "import random\n",
    "import math\n",
    "\n",
    "import numpy as np\n",
    "from numpy import linalg\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "import scipy.interpolate\n",
    "\n",
    "import dionysus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moving_average(g, k = 2):\n",
    "    \"\"\"calculates the moving average as used in [0].\n",
    "    f: 1 dimensional numpy array, input feature array\n",
    "    k: maximum size\n",
    "    \n",
    "    in [0], n is called S and i is called\n",
    "    \n",
    "    TODO: the paper states the method in general can be used for non uniform sampled values. how can this filter work then??\n",
    "    TODO: implement in in O(n) and not O(n*k) or in numpy (convolve?)\n",
    "    \"\"\"\n",
    "    n = len(g)\n",
    "    res = np.empty(n)\n",
    "    for i in range(n):\n",
    "        l = min(i, n-i-1, k)\n",
    "        res[i] = np.sum(g[i-l : i+l+1]) / (2*l+1)\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert all(moving_average(np.array([1,2,3,4,5]), k=0) == np.array([1,2,3,4,5]))\n",
    "assert all(moving_average(np.array([1,2,3,4,5]), k=2) == np.array([1, 2, 3, 4, 5]))\n",
    "assert all(moving_average(np.array([1,2,1,2,1]), k=1) == np.array([1,4/3, 5/3, 4/3, 1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_center_normalize(xs):\n",
    "    \"\"\"mean centering from [0] along axis 1. step 3 from the SW1Pers pipeline\n",
    "    xs: two dimensional numpy array.\n",
    "     axis 0: different points\n",
    "     axis 1: a point\n",
    "    \n",
    "    returned points are on a xs.shape[0] dimensional unit square.\n",
    "    \n",
    "    TODO: does this mean centering destroy everything, if there is a constant component in the signal?\n",
    "        ah ok, as this is a sliding window, mean(x) is a sliding average so this makes sense\n",
    "    \"\"\"\n",
    "    \n",
    "    xs = xs - np.mean(xs, axis=1)[:,None]\n",
    "    xs = xs / linalg.norm(xs, axis=1)[:,None]\n",
    "    return xs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_pointcloud(t, g, dim, window_size, num_points):\n",
    "    \"\"\"step 2 from the SW1Pers pipeline\n",
    "    t: times, must be ascending. if the data equidistant, you can set it to np.arange(len(g))\n",
    "       t and window_size are scaled, so that t[0], t[-1] correspond to 0, 2*PI. then window_size is called w in [0]\n",
    "    g: time series, data\n",
    "    dim: dimension, called M+1 in [0]\n",
    "    window_size: w, scaled with t\n",
    "    num_points: the number of equidistant samples from SWf used to generate the point cloud\n",
    "    \n",
    "    TODO: what are the boundary parameters for splining?\n",
    "        \"\"\"\n",
    "    # rescale t and window_size to 2π\n",
    "    scale_t = math.tau / (t[-1]-t[0])\n",
    "    t -= t[0]\n",
    "    #print(t, scale_t)\n",
    "    t *= scale_t\n",
    "    window_size *= scale_t\n",
    "    window = np.linspace(0, window_size, dim)\n",
    "    starts = np.linspace(0, math.tau - window_size, num_points)\n",
    "    f = scipy.interpolate.CubicSpline(t, g)\n",
    "    return f(window[None, :] + starts[:, None])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 1. 2.] 3.141592653589793\n"
     ]
    }
   ],
   "source": [
    "cloud = make_pointcloud(np.array([0.,1.,2.]), np.array([0.,1.,0.]), dim=3, window_size = 0.1, num_points = 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<mpl_toolkits.mplot3d.art3d.Line3D at 0x7f1ea88f0af0>]"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAADzCAYAAABNPxCMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABzEElEQVR4nO19d3gc1dX+O9vU26pbXZYs26qWZYNNS6jBgEsAGz5IICSElgRD+BECCaEEQnoIJPCRBiahumDANh0+AhhssNV779pdrcr2en9/SHc8O9rVzuzOyrK97/PosXd36u6ce8895z3vYQghCCOMME4dyI73BYQRRhgLi7DRhxHGKYaw0YcRximGsNGHEcYphrDRhxHGKYaw0YcRxikGhZ/Pw/m8MMIIPZiFPFl4pg8jjFMMYaMPI4xTDGGjDyOMUwxhow8jjFMMYaMPI4xTDGGjDyOMUwxhow8jjFMMYaMPI4xTDGGjDyOMUwxhow8jjFMMYaMPI4xTDGGjDyOMUwxhow8jjFMMYaMPI4xTDGGjDyOMUwz+6unDCAHcbjesVitkMhkUCgXkcjkYZkFLqsM4hRE2+gUEIQQulwsOhwMOhwNut5s1doVCwf6FB4EwQgnGT7OLsHKORCCEYHx8HNHR0ZDJZHA4HB6fEULmDAJKpRIKhQIymSw8CJzcCCvnnGxwu92w2WxoaGjw+jnDMB6uvkwmg8vlgsViwfT0NPr6+jAxMQGbzQaXy4VwV6IwgkHYvQ8hCCFwOp1wOp2sYXM/8zV7Mwzj8ZlGo4FCoWDfYxjGYzkQ9gTCEIOw0YcIbrfbY91O/wKZpem+crkcwLHBhC4RwoNAGGIQNnqJwQ3WAXNnbSnAPyYhhA0O0s9pTIAuF8KDQBgUYaOXEHx3nm9oMpnMY6afz8Xnwt82XC+AHpc/CLhcLsTGxoYHgTDCgTyp4Ha7YbfbfRo8BTV0sUYnZllABwH6xzAMjh49CpPJBIPBgOnpaZjNZtjtdrjd7nBg8BRDeKYPElx3nh+s4yOYNX0woIOMQjHzcxNCYLfbYbPZ2M+USiW7HAjFkiSMxYOw0QcBajzcYN18CNTo6bmkAj8oCAB2ux12ux0A2PQhNyYQxsmDsNEHCOrOi3HXgzH6UIFed3gQOHUQ/vVEggbJ6urq4HA4RAXFgnHvF2qw4KYHuQZut9vx5ZdfYmpqCtPT07BYLGxKMowTC+GZXgS47rzZbBb9wFPjJYRAq9UiMjISsbGxi3r9zPUELBYLm4Gw2Wyw2WwAZjwBShkOewKLH2GjFwg+GYaffhMChmFgt9vR2toKhUIBl8sFk8mEmJgYJCUlISkpCVFRUYt6EOAbNB3EuIOAXC5nlwJcJmEYiwNho/cDX7n3QFxuh8OBxsZGlJSUIDExkd3fZDJhYmICnZ2dsFgsiIuLYweByMjIRRkLoPBGFKKlwzTeQQeBcAXh4kDY6OeBNyothUwmE+zeE0LQ09OD6elplJWVIS0tjQ2SMQyD2NhYxMbGIicnB4QQGAwGTExMoLW11SOYFhsbC5VKJf2NSoj5BgEAmJqaQlRUFOLj48ODwHFC2Oi9gE+l9bZGFTr72u12NDQ0IDY2FqmpqYiIiJh3e4ZhEB8fj/j4eOTl5cHtdqOlpQUWiwWNjY1wuVxITExEUlISEhMT2dz7YgV/EJiYmIDL5fIYvMKewMJicT8xxwE0Ou9yueZNxQlZ0+v1erS0tKC4uBhpaWloamoS7abLZDJERkYiMTERycnJcLlcmJycxMTEBHp7e8EwDBITE6FWqxEfH+9Bx12MIISwmQH62u12w2KxeAQNw4NA6BA2eg7E5N4ZhvHp3hNC0N3dDZ1Oh9WrVyMyMpLdhxp9oOt0uVyO5ORkJCcnA5iJE0xOTkKr1aKzsxMKhYKNB8TFxS26SDq/3oB+z/Q6w4NA6BE2esxf9+4LvoyWimXEx8djzZo1khjdfAOEUqlEamoqUlNT2fNPTExgeHgYBoMBERERSEpKYsU3jrfB+LsGIYMAlygUHgTE45Q3eqpqA4grg/U204+Pj6O1tRXLli1jjZC/T6ij8BEREcjIyEBGRgYAwGKxYGJiAna7HYcOHTru6UGxA4+3QcDlcuHQoUOorq4GEJYWE4tT2uhpnryurg5r1qwR9bBw1/SEEHR1dUGv13u483wcD0ZeVFQUoqKiMDQ0hJqaGpjN5nnTg6GG2+0Oyvvh/kZyuZwdBJxOJ/t5WFBkfpySRs9156nxin0wqCHabDbU19cjMTERNTU1IamykwoMwyAmJgYxMTHIzs72mh5MSEhgMwOhSA9KscTgHsNbijCsKjQ/Tjmj5+fe5XJ5QPxxmUyGqakptLe3o6SkBCkpKX73Od5Gz4e39OD09DQmJiYwODg4Jz0oBaQw+vm8BaGqQqfyIHDKGP18MlZiDZEQAr1eD7vdjpqaGsFu8WIvuJHJZEhMTERiYiIKCgrmpAfNZjO6urqQlJSEhISEgNKDUhm9mNiLL1Uhu90Og8GAzMxMj7qBk30QOCWMfr7cu9gf2Gq1or6+HgzDoLCwUNQ6eLHN9P7ATw8eOnQI8fHx0Ol06OrqCig9KJV7H2hcgDsIOJ1OjI+PIzk5eY6gCPUETkZBkZPe6AOpe/cFrVaL9vZ2LF++HAaDIaCCG7qP0WgEwzCIjo4Wtd/xBMMwgtKDSUlJPqsHgzFYimCDgRQul8uDKAR4qgoBxyoITyZVoZPW6MXIWPmD2+1GZ2cnpqenUVNTg4iICBiNxoCM3u12o6uri9Wy79BZ0GZQ4vKqTCzNSoVSqQz4OhcavtKD/f39MBqNiI6Ohlqt9kgPinHNfUEqo/d2nFNBVeikNHqxMlbzgbrzycnJWL16tUfUWGwA0O12o6enB5HxanQiE7u/HEXzqB2AHV8M9eH/rR5GlAIsrZaumxfLTO8PND24ZMkSEELY9GBXVxfMZjPi4uJgtVpht9uDygxIMXDQ4/jLtgAnn6rQSWf0oXDnV6xYAbVa7fEZbT0lFJOTk3i/oR9fTUTg86ER2JxulKTF4J4LlyI1VoWfvNaKZzuj8cQVK2A2TnusmxmGQUJCwqJg1AmFr/RgY2MjOjo64HA4EB8fD7VaLTo9KMUSAZhx78Ucx9cgYDAY0Nvbi5KSEnYQGBsbQ2Zm5oJwH8TipDH6QKi0vuB2u9HR0QGDwYA1a9Z4fSCFzr5j01Y893+teLN5HGNmghilDVuqMvHNqgyUpEax12lzunHv6224f38XfrNlOZsCtNls6OzshF6vh0ajQXR0NJKSkqBWqxe94AYXND0YERGBiooKMAwzb3pwvupBKd37YAqUuNkfOhBRL/PHP/4xHn30USxfvjzo65QaJ4XR09z7wYMHcfrppwdlCBaLBfX19UhNTfVw5/mYr8rO6XLj/zrG8cqXg/i4cxxuAtTkJuC6/CicXRiH4oI8OJ1OD0/hsvJ0jJvs+P37PUiJVeInFywFwzCIiIhgZ8SMjIw5jLr4+Hh2EFjstfbAsei9v/QgrR70lh4M5Zo+ENCAINezNJvNiI2NDfrYocAJbfT83Hsw7i/Vrevo6MDKlSuRlJQ07/be1vS942bsPDKE12pHoDXakaACtlYk4/qzl6EgJQb9/f3zHvO607KhMdjx/KEhpMVF4IZ1OXPOyXWZ3W43jEYj9Ho9mpqa4HQ6F32tva/fyFf1IF3myOVydnAT65b7gtvtluQ7okbPhclkChu91PAlYxUopZYKVfhy573tQwiBxe7C281j2HlkGIf7JiGXMTgtJwZXFwP/87VKJCUmzNmHew98zsBd5xdCZ7Tjjx/0ICVGhY0V6T6XEjKZjGXU5efnz5ktZTIZm0ILNhAote6+P/CrB+12O5senJiYYJl186UH/cHlcvkVNRF6HP4gRLUPFyNOSKP3JWNFKbVi1mlmsxkmkwmpqalYsWKFoIeHEIJ2rRW7anX47ysDMNpcyFNH4Y5zC1EaY0S8kqC0tGbOLCIk4i9jGPzyshLozQ78Yl871DFK5At8Lvmzpd1ux+TkJEZHR2E2m1FXV8em0GJiYk6YeAAAqFQqpKenIz09HWNjYzAYDFAoFB7pQbGxjlDGBqTyIkKBxXlVPuAv905164Qa/djYGDo7OxETE4OcnBy/D8qk2YE36kew8+gwWkeNUMkZXFyWjiuqs1CaqkJDQwOWpC3xeSyhwT+VQoY/XbES33m+DnfuasavLszAijT/JJ45x1GpkJaWhrS0NBgMBixbtgx6vR69vb2s+0kHAX9R5sWUOXC73YiIiMCSJUvmTQ/6qx6UyuidTuccgs9i+a684YQxeiG5d5pG80dwcbvdaGtrg8Viwdq1a9HQ0OBzBna7Cb7oncDOI0N4p0ULu9ON0iVxuOtr2ahSu7CmshQajQa1tc0oKytDQkKC1+MA4ph1sREKPHVVOa59rha/eH8Uf9iQhSVLBO3qE1FRUcjKykJWVhYIIWw8gFthR1No/O9wMT3IfGP1lh7k3xsNeCYlJbHLN29r8UCvh3+cxfR98XFCGL3Q3LsQhVqz2Yz6+npkZGRg+fLlrMfA3290yordtcPYdXQEgxMWxEcqsLV6Ca6ozsKKzDh2fdnW1gaj0SgoFiCWZJMSq8LTV5Xh2n8dwX3vjODFnCykxEoToWcYBnFxcYiLi2Mr7KampqDX69mAI5cktJjgL0/v7d5oenBoaIhND5rNZkliFUImmsWERW30YnPv/spkR0dH0d3djdLSUo8Hmc2vOt34qF2HV48M4ZPZVNvpBUnYfu5SXLAiFZHKY6O5w+GARqNBXl4eqqurBY3qgTDr8pOj8dAFS3D3W0O49aVG/OtbFYiJkP5n4wb9gGPRc41Gg46ODiiVSthsNhgMhuPelUfsetlbenBqagrj4+Nob2+HXC73mR4UAr7HECzjMNRYtEYfCJXWF0vO5XKhra0NNpsNa9asmTMqDxlceO2DXhxoHYfe5EB6fAS+f1Y+rli1BDnquWtpKosVFxeHwsJCwffENXoxRrM8NRL3nZOOBz8YxfZdzfjrtjIo5aGle/Kj50ajEQ0NDWzgjMpu0cDZQiLYtbhcLodarUZ0dDSWLVsGuVzuMz0opHqQb/Qmk0lQIdXxwqI0er7yiVAD8eamm0wm1NfXY8mSJR7ReZPNiQNNM6m2owNTkMsYnFuSgiurs3BmUTLkMu8VYt3d3RgfH0dpaSmGhoZE3RfX6A0GAywWC5KSkgQ9wGuyo/Hgpcvwszfa8bM32vCrTcshW8DZNiIiApGRkSgtLQUhhO3K097eDpvNFjClNhBITc7xlR4cGRlBW1ub3+pBvtEbjcZFm6MHFpnRE0JgsVgwNDQkKJrOB9+9HxkZQU9Pj4c7Xz80hVe+HMK+xjGY7S4UpETj+qp4XFmTh6KcdJ/H5jatoFpzgVbZ9fX1YXh4GLGxsejq6mIfKrVa7TWVRgeLTRUZ0BrtePzDXqTGRuCu84V7GcGCL1HF7crDV9xxu90eJCF+ZFvKawkGvgYPbnoQ8F49SAeB6OjoOXn6xUzMARaR0dPcu91uh0ajQW5uruhjUPfe5XKhtbUVDocDa9eunSmAmLbhN++0482GMUQpZWyqrTonAe3t7UiM8r2Om5ycRFNTE9u0gp4rkCo7jUaDpKQkrF69mn3o6ENFU2lxcXFQq9VeqbXfXZcDrcGO574YRGqsCtedni36ewoU8wVQuWtmp9OJyclJ6PV6dHd3syQatVotidFLTZ/1B1/Vg93d3TCbzXA6nYiNjYVCoUBkZCTMZrMgYs5bb72F22+/He3t7Z0A/k4IeYz7OcMwXwOwF0DP7Fu7CSEPzX72DQCPA5B723c+HHej51NplUqlqOo1LmQyGSvplJWVhZycHDjdBP/6rA9//rAbTjfBbecU4Ib1eYiNVHjs582ACSHo7+/HyMgIVq1a5bFOExuUM5lMaG9vR2RkJMrKyuByudhz8h8qg8EAvV7PtrFSqVSIiIhgH9K7L1gKncmO373fjZRYFS4pSwvo+xIDMfeqUCiQkpLiUTREvQBa20D5AdHR0aJnbamMHhCvnOQtPfjll1+ycaP6+nrs3bsXSqUSGo2GnST4cLlcuO222/Duu+9i6dKlKwEcZhjmdUJIM2/T/xJCLuVdgxzAXwBcAGBwnn294rgavS8Zq0CEKoGZtdTk5CRWrVqF+Ph4HOqdwEP7WtGhMeGc4mT8bEMJcr0E5rwVzzidTjQ2NkKlUmHt2rVexRaEGgKNgOfn58NoNM67La1G41Jru7u7MT09jSNHjkChUECtVuOnX8/ChMmBn73RhqRoJdYXzl8rECyCcamp2EZaWhpMJhOKiopYL4ASaeggIIQWK6XRBwuaVcrLy4NcLkdJSQl0Oh3ef/99XHXVVdi4cSO2b98+Z79Dhw6hqKgIhYWFIITYGYZ5CcAmAEIMdy2ATkJI9+w1iNn3+Bm9r9x7IH3fXS4XWlpaYDKZkJeXBysTgYd2NeKN+lFkJUbir1dX4tySFJ8PLZ8eazAY0NDQgPz8fCzxwYgR4t4TQtDZ2YmpqSmsWbMGJpMJBoNB1L3J5XLExcUhIiICubm5sNls0Ov1GBsexLcKLRiblGH7q43429WlqMxV+z9ggJBqLS6TyRAdHY3o6GiPOnu9Xo/m5mY4nU4PkpC31NxiI75wB6HIyEhkZWXhoosuwr333utzHxq34mAQwGleNl3HMEwdgGEAdxFCmgBkARgQsK9XLLjRS1n3DhxLJWVnZyM2Lh57mvR4vrYHdqcbt55TgO+fmY8o1fzrNq4BDw0Noa+vDxUVFfMGY/zN9Ha7HfX19UhISGBLdE0mU0D3yD1XREQEMjMzkZmZiZWEIK9wAt97uQW3vtyEe9eqsGyJmg0ySd3MMlhD82as3jwbShKiJbY0HhAfH8/+Votlpgfm3peQslofzw7/zSMA8gghRoZhNgB4DUAxAG8/hOBReUGNXkoZK+CYgZaXl6Nn2o2f7K5Ht96GM4uS8fMNJchPFpYrlclkrDvvcrnY4J+/fXwZ/dTUFBobGz0Cf4D0ApcMw6AwU42/X7sK395Ri780MXhiaQImJyfZKjsaEJSiyi4URs8HzaFTpSKHw4GJiQmMjo6ivb0dERERsFqtsFgsi1ZExGg0spF/X8jOzsbAwIDHW5iZzVkQQqY5/9/PMMxfGYZJwczMnjPfvvNhwYyeBuukkLFyOp1obp5ZvlRV1+DpT/rx90/7oI6S475zUvGtr1eIOr7D4cDAwAAKCwsFpwp9VcwNDg5iYGBgTuCP7hMK3fvClGg8ubUMN/6nHve9M4h/XluJoiI57HY79Ho9BgcHYTab0djY6CFUudAIZIZWKpVs0RBwTORkaGgInZ2diI2NZT0BMdJUodQcFJKyW7NmDTo6OtDT04PCwkIVgKsA/A93G4ZhMgCMEUIIwzBrAcgAjAOYBFDMMEwBgCFv+86HkBt9oO68r1mBrrdzc3OhZ+Jw5d+/QqfWhMtXLcGNNWo4zNOiDH50dBQDAwNIT08XlSbkGyKNK7jdbqxdu9arax0oI08IqrLj8btvrsDtrzbhzl3NeGJrKVQqFatWazQakZ+fD71ezxJq6No5KSnJr2ezUDO9P0RFRUGlUmHFihVQKBQwGo1z2nL5KhriQqolgrd7EuLeKxQKPPnkk7jooosAoAXAPwkhTQzD3Dx73KcBXAHgFoZhnAAsAK4iMw+Qk2GYHwB4GzMpu3/OrvUFIaRG76vu3R+ocXC3J4RgaGgIAwMDWFFahue/0uLp/7YhNVaFv11bhbOLU2Z05IzCIv9utxvt7e0wm80oKipi1U2FgnttFosFdXV1WLLEd1kt977EQuh+5xQn4/4NxfjFvg78Yl87HrmsxCuhJjc316PApq+vDwzDsG61N+rpYjF64JjBcgtruPdEiTQAFkRyi39PQgU0NmzYgA0bNgDAUvrerLHT/z8J4Elv+xJC9gPYH8g1h8To+bl3se68XC73YDk5nU40NTVBJpMhbWkZbny5BQ1D09hSlYn7Li5B3GzOXShhxmq1oq6uDmlpaSgpKYFGo4HVag3gTgGdToe2tjaUlpb67fcm9ZreG75ZlQmt0Y4n/68PqbERuOPcAq/b8QtsuMo0BoMBUVFRrNssFY881Ew6/j05nU5MTEx4KAvTgU2pVIZMQOOUY+TR3PvRo0dRWVkZ0I9MjV6pVHq48wfHGDzyzGFEKOT487ZyXLTSM1giRJaaGilX1joQdh0hMx1re3p62AYY/rAQRg8A3z8jF1qDHf88OIDUWBWuXZvldx8u9ZTLOuvs7ITVakVUVBTbAy7QMlIpU21CjqNQKOZ05NHr9RgYGMD09DRcLheGhoaCKho60fTxAImNnpt7N5vNAR9HLpfD6XRiYGAAg4ODKChZid98MID9jWM4vSAJv/lmKdLj5wZt5jNeQmZ6yE9MTMwxUrFG73A40NDQAEIIVq9eLXjGCLV7z93+pxcVYdzkwG/e7UJKrBLJIs/HF+AcGRnB0NAQ6uvrAWBOGk0IpGpSESi46U6DwYCenh5W7txqtXqQhIQWDfky+ri4uFDcgiSQzOhpOg6YMaJA9OooGIZBa2srIiMjEZuzAt96vgnDU1bccd5S3HhmvtcKOMB3PT03Z15TUzPnwRNj9NTzKCwshNVqFd0sgRrv6Ogo2zwxFMq1chmDxzYvx/dfqMe9r7dhe5USawI8lkwmQ1xcHOLj47F8+XKvaTTqNs9Hq/UnfrGQIIRApVIhJyeHLRqiJKGhoaF5i4a4OKVnev66nbroYo1+enqmu0tOTg6OTEXjl7u/QnKMCs9/ZzVW5ybOu683935iYgLNzc1YtmwZ6+Z520+I0Q8PD6O3t5cl7nR3dwu+L+BYmq+1tRVmsxkZGRlsTp2bn+aXbwbqIUQoZPjzlaW4/vk6PFFrxtoqI5ZnBPYwcl1zfhqNW4BisVh8zpiLiUnHfzZlMhkSEhKQkJDAFg3RQGd3d7fPGntvz7jFYjl16um5xkNddKFuEiEEAwMDGBoaQmpqKp4+asS+5n6cWZSM336zFOoY/8fhzvSEEPT19WFsbAzV1dXzrtn8GT01VLvdLoi44wsOhwNTU1NISkpCZWUlHA6HR1EKlaoyGo0elXbBICFKiaeuKsfWZw7hlpca8e/rq5CVGFirJV8GS2m1WVlZXmdMbnBtsRi9v+i9QqGYoyzMDXRGRkZ6xIS4IIQImuwEVNldA+Ansy+NAG4hhNTNftYLwADABcBJCKkRcNsz9yZ0Q7GgM70QOBwONDU1QalUYu3ateju7kZ+oh0//Fohbj2nADIf7jwfdKZ3OBxobGxEREQE1qxZ49elnG8mpZH+9PR0wRLZ3jA1NYWGhgZERkZi6dKlcwYZ7nqTX2lnt9uhUCgwMTGBhIQE0S5yRnwE7lytwm+PuHDziw3YcV0VkqLFBeOEehreZkwaQR8fHwchhDWY4ynDLTZlxw900nLokZERWK1W2Gw2JCUlCfbKBFbZ9QA4hxAywTDMxQCegSfH/uuEEJ3gm5hFyIxeoVDA6XT63W5qagpNTU0oKChAZmYmgJkB48oKNftaKBiGgdPpxOHDhz2O5w++Znoqi+WtgaUYDA8Ps3Th9vZ2v9vz+eharRYjIyNstR41GjHptKxYGZ7cuhw3vtCAH7zciL9dU4FoPzUJXATqmnMj6OPj49BqtVAoFHO0A4RW2EmFYJRwGYZhvRv6zCUmJkKv1+Omm27C4OAgtm/fjosuuojm4edASJUdIeQzzi6fY4ZuGzQkNXpva3pfoO736OgoqqqqPB5eMV4CF0NDQ7BYLFi/fr2oQArf6Akh6O3thUajwerVqwPuPEqltqk2H32PQqgRyeVyREdHo6ioCMDMGlqv17PpNKHMulU5Cfj15uW4c1cz7trdgsevXClYa08qco5KpfLQq+dX2FEFXl/BM6lSnlIKcSgUCnaQPnDgAM466yxs3LgRbW1tPvcTUWVH8V0ABzivCYB3GIYhAP6XEPKM0GsO6Uzvy3C57re3WnWxRu9yudDc3AxCCGJiYkRHTrlG73Q6WTdcyNLAlzHY7XbU1dUhOTmZldoW4vl4A//43NJUPrOOW2QTFxc3Z9/zSlJw3zeK8PCBTjy0vwMPXbpswVxs/nflrcKOr7jDD24uVG96oaAiJxR2ux3R0dE477zzcN555/ncT2CVHQCAYZivY8boz+S8fQYhZJhhmDQA7zIM00oI+VjINYd0Te/tIafSU0uXLkVGRobXfWUyGcvm8wcqfJmdnY3s7GwcPHhQ9LVSozcajaivr5+3jp6/nzejp1V2/IxBMOQcX/t5Y9bRIhsqV61Wqz08jK3VS6Az2vHUf/uRGqfCj77mnbXHP3+oabj8tlz84GZsbKxkGvxStZ3ip6Wphp4/CKmyAwCGYSoA/B3AxYSQcfo+IWR49l8NwzB7MCOssfBGP597T13msbExrxVoHhelUAgi91Ad+7KyMsTHxwd83TKZDDabDfX19SgvLxdMrKAzD3fGGBoaQn9/v6RVdmLALbIhnE4vVqsVhw8fZtNON52RA43Bjr99OoDU2AhcXRNk+xwBEDtw8IObRqMRWq0WZrMZhw8f9iu2MR+k6m4TqBKuwCq7XAC7AXyLENLOeT8GgIwQYpj9/4UAHhJ6zSF17ymf3W63o7GxEVFRUV7deT780WnpWtlqtXrVsRcDWnhjt9uxbt06UcfiGjF//e7tIeQbvVAjCIbJRwtSNBoNqqur2Uh6Z2cnNi5RYUgXjV+93Ql1tAIXrfSttSfFTB+Ma07vRalUwmg0orS01ENsw9+yxtu1hEJcUygxR2CV3f0AkgH8dfZ+aGouHcCe2fcUAF4ghLwl9JpD7t5TcgxfUMLfvr6MntZTp6WlsWtlPoQ+oHR2p9xrsYMHXRbYbDbU1dUhJSXF5zUBc9fmC52uksvlHoKVFosF96h1uHtfP+55rRXG8VGcVZKBpKSkkPSyk4KRRwcOvtgGf1njrxlHqIxeqBIuIKjK7nsAvsffj8xo41UGes0hc+9lMhn0ej30er1fcgwfvoxeq9Wivb0dK1euZNewfAjtXEsHo5KSEqSkpGBsbEzw9VEwDIOpqSl0dHSwxxGzbyi2FYOoqCgszcvBP7+TgW8/V4vffT6NpJgIJMyuNalRBbN04iKUAwd/WcNvxsHPcITSvV+sfekpQjLT2+12dHV1wel0Yt26daJHVG/xAK7A5HwsP3+cf8KRtRY7GPFhs9nQ0dHhN0YhBUIZC0iIUuLpq8tx7bO1eOS/E/j39VVIiZZDr9djeHgYra2tkMvlUKlUsFqtAacwpZrphZCt+M04+NoBLpcLCQkJiI2NDeqaTrTuNkAIjJ7OoDk5OZiYmAjoC+XSaakLnpiYyApMzof52lW7XC62Ln/NmjUBj/RcWm51dXVABk8VhYQsKRZiGZCZEImnri7H9TtqcdOLDdjx7SoPBlp/fz8mJyfR0tIiKJ/uDVLFBcQ+U94yHPX19dDpdBgYGPAgO4nV3eN3txGimnO8IanRm0wmdHR0YPXq1SCEQKcTzRAEcMxw+S640H29setoai8nJwfZ2YETm+j6PTU1FYmJiQENapQLYDKZ2ABUcnIy4uPj561QCzWWpcXg8StLcfOLDfjBK0342zXliFLKwTAMIiMjkZiYiLy8vDn5dKVSyRrNfNTa4x0MpFCpVFAqlVi2bBlUKhUsFosH2Yn25fMW2+CD772ccu59bGws1qxZA4Zh2CYWgUAmk7EdYQKJB/CNntJXy8rKgsrzUo4BHYRoTb0YuN1uHD58GDk5OVi5ciUb7KRudExMDJKTk6FWqxeUlkqxJi8Rj21ejh/vasHde1rxxytWQiHzlC/j59OtVisbRTeZTKzRUIUa7r0H695LVZ7Lldzik52mp6dZsQ0ArJfgq+6BOwiZTCZkZfkXLTmekNy95z4YgTDQqECF2+0WxIjjg5vuExML8AdvKrc0Ty8Uer0eZrMZa9euRUJCAux2u0eZKg1AjY+Po7m5GS6XC0lJSYiMjAy4608guGB5Ku69yIFH3u7ELw904BcbiuedpSMjIz2otdRoBgcHARwT3JBilpa6Yy0f3L58wDEJbjpx8LUD+IO+mOj98ULIUnaBdKqhxTeBCFRwz0sVfPjNJvzB24NN1++0GSZ3/Somf06Dh9HR0UhMTPS6HzcAlZeXx3oBo6OjmJiYQENDA7sUCDSYJhRX1SyBxmjD3z4dQFqcCluKhXkdDMN4VNlRoxkZGYFOp0NUVBTsdnvAElULTcP1JsHNbclls9kwNjbGagcIVc2hZbUulwtdXV33eCmrZTDToHIDADOA6wkhR2Y/C7h5JRDClJ0YEEIwODiIwcFBVFZWIiYmRrRABYVcLofBYBDNDfCW6qPr97S0NOTl5QWkuON2u9HS0gKXy4Wamhp88cUXgu+FVqhFRkZCLpcjLy8Per2eHYQSExNZ5Z1QKNL88Jx8lq6rdKbi4mXiJaC4RtPZ2YmoqCiWyGS329mAoNCOPMezeSUwk+bMyspCVlYWXC4XDh8+DLPZjKGhIezcuZNth/a1r33N58DMLavNzs5GRETE1V7Kai/GTDebYswU4jwF4DQmyOaVQIjcezEzPG1cIZPJfOrFi4HRaMTY2BhqampERdX5ngldvy9fvpxdu/Lh7179DRpCQfejunU5OTlsoJOy6yIiIth1tlSNLBiGwf0blmHc5MCTn2uRGCnHlUEWd0ZFRUGtVrP3QFNpPT09Xgts+FhMklvAzKBWUFDA/t1yyy344osvsGPHDrz88ssoLi6esw+3rHYW3hpQbgKwg8w8YJ8zDJPIMEwmgHwE0bwSOM5da2kfumAj6sCxZhNWqxVFRUWi02jc9Tn1OvwFEecz+unpaTQ0NIgm7QgFn11Hy20pGYXrBQQzkCpkDH67ZQWu+9eXeOzjMSzNSUd1TmDBUL5rzmfV0QKbvr4+j1p7tVrNxmMWUx87frqORvsffvhhtgzaGwSW1XprUpnl433BzSuBEBu9t4IUipGREfT09Mxb4CI0xcNtNhEoSYb2s+vo6IDL5RKUx/fl3tN7q6qqkiSoI8R74kag6Qw6Pj7OptTsdjtMJlNA/eCjVXI8cF4G7jowjB+80oQd365EUar4+/I3S3tTDxofH0djYyMru+VyuRZNHtwbq08IOUdgWa2vJpVBNa8EQuze85tWAMcCYzabbV69Obqvv+opfrOJ/v7+gCLdhBDU1dUhMzNTsCvurYCmo6MDRqMxKC29YMGfQa1WK7766itWuFJMOyuKhEg5HrlwCX58YAQ3v9iAf1+/Chnx4lKKYvL03Fp7ruxWX18ftFotdDqdaPUgqRFoowuBZbW+mlSqfLwvGCF9KrlNK4AZF7S+vh4ZGRl+9eZ8yVlTEELQ3d0NvV7voWMvpOEFHxMTE5icnMTy5ctFLTO4cQCn04m6ujrExcVh1apVi0YAEphJqalUKpSXl8+hpNIBIjk52S+xJjNehaevLsP1z9fhlpca8Oy3KpEQJbxIKRhyDg1qmkwmREVFIS4uDnq9Hh0dHQH15ZMC3mZ62hhkPnDLamdz+t4aUL4O4Aeza/bTAEwRQkYYhtEiiOaVQIiNnquTR/Oc8xXLcDFfpR3N5cfExMxpNiFGgAMAq8BLGXFiQJcvJpMJdXV1onT5xJ5HqgeVT0m12WwYHx/3INYkJyfPYaNRgy1Jj8XjV5Ti5pca8KNXm/C/V5cjUrmwNFwqH8ZXDxofH/dQD0pOTvYaEAxVsQ2Fv5gDt6x29hl/xUtZ7X7MpOs6MZOy+87sZ0E1rwQWYKZ3Op1ob2/H9PS0KIKMrxmb22zCm/KOPw+Bwu12o7m5mSUBtba2il4WMAyD6elpNjYhdNBYCEqtUERERLDEGipfPT4+jv7+fg/j4V7z2vxEPLpxOe7e04J79rbi999c6bMBCRehIud4G8j4ijt0KRARERGyslpCiODfllNWCwCPzO7PLaslAG7zti8JonklEEJGHkVzczPS0tIEE2QovM30/GYT3iDEvaey1hkZGcjNzQXDMKLJRIQQTExMsOt3oZRZOmuLLa1diIGCK19dWFjI1qj39/djYmKCLUZRq9X4xspUjJvseOydLjz6did+9o0iv/ckRbpNyHfnTXFnfHwcTU1NcLlciIuLg9PpDNr4fc30i2lp5w0hm+n1ej00Gg1yc3PnTV/4AtfoxTSb8EeYoUU8fFlrMZRaWq1nt9uRm5sriiNPDXh6ehqdnZ1ISEjwu54+XuDWqPf09IBhGJjNZpZeuz5FjWuqU/GfIyNIi1PhpjPz5j3e8aiy46oH5efnw+l0sgzHL7/8kqXVUn6D2Ao7/kx/IiAkXWu7u7uh0+mQlZUVcMqKGj23rbSQZhO+3Hsy20FneHjYa/5dbJvrzMxMqNVq0UFDhmEwMjKC/v5+FBcXw2w2o7e3F2az2WM9fbwi//MhOjoaaWlpLL1Wr9fjslwzukdkePL/+qB0WnDt+kKfS7jjVVrLhUKhQEJCAgwGA1asWAGLxYLx8XGPCjuhvwHf6K1W66JuZ0Uh+ZPV398Pu92OmpoaDA4OBlxpJ5fLMTU1hfb2dlHNJry593T9TgjxmX8XYvSUpUevZ3h4WFTQkMx2RhkdHcWaNWvgdruRkJCAzMxMj2BUb28vFAoFUlJSkJycvGDuvT9wDVapVLL19k+WuHDri/V4/DMNiGUalakza2waHKVGuhiMnn+MqKgoVknZm5w4vQ9vunsul8vDyxOqhHu8IbnR03prYMZwxRgFBSEEk5OTMJvNWLNmjajiEr7xelu/C9mPD28qt2KMkdbQA0B5eTkUCgXb5ZeenxuMslqt7AxksVjgcrkwPj4eNMMuUMxnsCqFHH/eVoEb/l2Hp+vNePqqUsRHz7jRbW1tiI6ORnJyMpxOZ8hltIWAzx2h8CcnHhMTwy4FIiIi5hznRFDNAUIcyJPL5awirlA4nU40NjbC5XIhNzdXdDUZ1733tX73Bl9GT9VyqfIu1+UTavRWqxW1tbXIzs72ePDne3gjIyPZwg6LxYKmpia2ukulUknOsw8W0So5/rKtDN9+rha372rBjm9XYfnyVBBCYDabMT4+DrPZjNraWnb2DKQvn1QzvZCB05vuHi17pqlohULBFguJLavV6/XYtm0b3nvvvQ4AvQC2EkImuNswDJMDYAeADABuAM8QQh6f/ewBADcC0M5ufu9sZH9ehNzoxbj33GYThJCAvARKp+3v78fw8LDgtlTeovcOhwN1dXVISkpCSUlJQFV2dElA+QljY2MghMDtdsPpdLIzl0wm8/kwy2QyKJVKtniDrkOpdHeoq+0AYTNscoxqRmvvudpZ1l4V0uIi2EIhjUaDiooKTE9Pe/Tlo6IhQgYwqd17oeCXPbtcLjQ0NMBgMOCrr77CxMQE3nnnHbhcLsHeyGOPPYbzzjsP7777bjHDMPcAuAfHutRSOAH8mBByhGGYOABfMQzzLqeq7o+EkN+JuZcFI+f4A21cQbn4o6Ojor0EACxnm7alEuoK86P3dABaunQp0tPTfe4z30xPG1fyhTOdTidkMhkUCgXcbjcIIXC5XOx3JZfL2YGAnocL7jqUSlfRajtqRFLX3At9kHOSovDUtjJ859/1uPmlRjz7rUrERyrYYyiVSrZQiMY46ABGO7/Op70nhdH7cu/FQC6XQ6lUIj8/HzExMRgZGcEHH3yAI0eOoLKyEnfffTeuvfbaeY+xd+9efPTRR/TlcwA+As/oCSEjAEZm/29gGKYFM0U3gqvq+FgQGu58oO4zXb9TFphQkg0X1I2Wy+UoKysTtfbjMvk0Gg06Ozv9drvxNdNTxR6DweCxJKCGMzIygoyMDERFRbEPNtf4uX316Dl8DS5c6SquEdGae+pKL2QgcGVmHP54+Urc9nIjbn+1CU9fXY4IxdwBjCtVRUttudp7KpWKXUPTQVNKVl+w4EbvMzMzccYZZyAlJQX3338/jEaj3/3HxsZYBucsxXZe8QeGYfIBrALAFWX4AcMw3wbwJWY8gglv+3IRUvd+viaWwLF68+Tk5Dl8dbEcer1ej5aWFqxYsQKtra2iHwx6vq6uLpbP74896G2mpwG76Ohoj3uiBl1cXAyNRsOuC5OTk5GSksKub7mDAB0ATCYTgJnlBtV186XVxjUip9OJyclJaDQamM1mNDQ0sAOEWP09sca2vjAJv7ysBPfsbcVP97bit1tW+N2Hr73HF6ykEmPBDmBS9bHz1d1GoVCwclvnn38+RkdH5+z7yCOPiDoXwzCxAHYB2E4ImZ59+ykAD2Omyu5hAL8HcIO/Yy0IDdcb/CndCo0HUHnm0dHRoNpKAzMlscnJyXP4/L7AN3qLxYLa2lrk5uZ6iCNSgyeEICoqCvn5+SxRZHx8HENDQ2hpaUFsbCybplOpVJDJZNDpdOju7sbKlSs9PAvqovoaAACwab+UlBRMTU2hoKCAZaa53W5BKrzcexA7kF5SlgadyY7fvdeNX7/bhfOTxBkrV6XG7XazA9jRo0ehUCjYAUJsuXAoW1rxA3nvvfeez/3T09MxMjKCzMxMzApkaLxtxzCMEjMG/x9CyG76PiFkjLPN3wC8KeS6F9y9p0bqr9mEEKOnLaoZhglIRJPCYrGgs7MT0dHRWLHC/4xEwTVCfsCOgs7YdHsuFAqFh7a8wWCATqdDbW0tGIaBQqGAzWbDqlWrPGZmrhfAPT59AH15AXz9PW4zi9jYWDagFoyAKB/XnZYNrcGO574YhL1YgbVrAzsOrQOgsRqqwEvLhYMh1QQK/kBoMpl8qix5w8aNG/Hcc8/hnnvuAYDrAOzlbzOrlfcPAC2EkD/wPsucXfMDwBYAjULOu6BGL6bZhD+j5wpn5ObmBnyNdFmQm5srOnBIZ3oasOMOYjRCT4tM/M1E/PrxlpYWGI1GREVF4ciRI0hISEBKSgrUajUUCgVr2PQcXOMX6gVwVXgpP72xsRGEENYLoKSUYNbSd55XAK3Rhl1NWlTVjWJzpfcW5WLAVeClstXcKjs6gHmrsguVzp63mX4+3HPPPdi6dSt++tOfdgDoB3Dl7DGXYEbwcgOAMwB8C0ADwzC1s7vS1NxvGIapwox73wvgJiHnDemanv+FiGk2MZ/RU0MVWqbrC7SsdvXq1TCbzYLaY/MxNTUFh8MxJ2AnxuC5oDyFuLg4lnZMmWLU1VcqlUhNTUVKSgqio6M9BgB6bu4A4C/LwOenU4otJaXExcXBZrMF3seAYfDLy0rQOzKOB/a1IzlGhbOKhDEsBR2fJ1vNl93i6/BLZfT871SoEi5FcnIy3n//fWBG/JJ73GHMlNWCEPIJvKvlgBDyLXFXPIOQzPT8hyyQZhO+lgZ9fX0YGxsLav1OC3icTifrcVitVlHZAqfTiba2NhBC5gTsAjV4q9XKDozcunwuU6y4uBgWi4VVDKJprpSUFJYkQtNJ1PjHxsagUqnY7IQ/L4BLsaXLjra2NnR0dLDdbHzVqvuCUi7DbVUqPNmswI93N+Mf11SgPEuaxph88KvsqBdAC4VcLhfi4+MlyQRwIXamP14IqXtPCIHVakV/f7/oZhP8dBh/aeBvpPb1g9rtdtTW1iI1NRX5+fnsNkILboBjAbuMjAxMTk7OidDTc4t5oKanp1leP52xfCEqKgo5OTkeqri0oy+lvKakpCAiIgJDQ0OYmJhAeXk5O5Byg4H+iEF02REbG0vlmj1q1ePj49kBx99aOkrB4K/byvCt52px26zWXn5yaLnqDHNMhx841sdOo9Ggr68v4FiGt+dLaG/6442QGT39cgGguro6oE411Fug6/esrCy+iqjPfb3lYqlC7bJly5Camup1H3+gWYeVK1ciKioKExMzaVEy25CSHksMNBoNuru7UVlZKbpgg6uKS6miOp0OjY2NMBqNiIyMRElJCZRKJRiGmZMSpAMV/T+fGERBH3IuNZWKbuh0OlZ6y19EPSV2hrX3bZa1twopsdIFDv1BpVIhIiICRUVFiIyM9IhliMloeFsimM1mUe798UJIjH56ehr19fUoLi5Gd3d3UHlVSjQRs3731q6aMv4qKyu9jsZCjJ4W3dCAnc1m84iki53daSZDp9Nh9erVgjrYzgcaoY+KisLk5CSysrIQFxfHRujj4uLYlKBSqZwTC/BGDPJHD+bOolR6iyvASSPq3N8iTx2Fv2wrww3/ntHa+9e3KhEb4f9RlIpgxO1jx49lCO0r6C0DcEq796Ojoyz1tL+/P6AUCSEENpsNnZ2dotfvXAOm7Dgq1+XLsOYzekKIB2uQ68Y6nU6PaLlQ0C4vbrcbq1atkiyabLfbWa9oyZIlAMAWjExPT0On06G/vx8Mw7AeAu3R7o0YRL0AborQ17XypbdoqXBPTw8rw202mxEdHY2yJXH4w+Ur8cNXmrB9ZzP+uq0MKsX830Go8usUQvoKqtVqJCQk+JS/PmVn+pKSEtaAKEFHzHqJrt+pfl0gSwPKZa+vr0dMTAyqq6vnNUpfcln0GLGxsaiqqvJYv8tkMkRFReHQoUNsSi05OdnvAEePmZSU5BFXCBYmkwkNDQ0oLi6eky/mrm2XLl0Ku90OnU6Hnp4emEwmJCQkIDU1FWq1GnK53MML0Gg0sFqtUCqV7ADgLxjIL1O1WCw4cuQIy65LTEzEiuRkPLChCD97swM/e6MNj21eDtk834VUgTchg4c3XgO3kaVSqYTT6YTVamUnJJvNJorpKKTKbvZaegEYALgAOAkhNbPvqwG8jJmuNz735yNk0XsKsZV2NEiWk5MDo9EY0I8sk8lYhdq8vDx2xvN3zfyZnl4L/xjUFQaA0tJSEELYlFpPTw9UKhU7i/LJRxaLha0k9FXIEwgmJyfR2tqK0tJSQbONSqXymJVp0U5XV5fH9dPmk6tXr2YHM1/EoPkMia6lKyoq2PONj48j0zaBq5ZH4qVmLRIiZbj3G8t8/uZSptrEPldUgjs1daZcWKPRYHBwkK1x2L9/P5RKJRwOh2CvVmCVHcXXCSE63nv3AHifEPKYgP2P3YugqwsCYirt6PqdNq4YGBgI6AdyOBxoaWlBZWWl4BQh372nATt6LRTeAnYMw3jMamazGTqdDi0tLXA4HGw0HQDLLxB6XUIwNjaGvr4+VFVVBZTGpGw3qjlgsVig1Wpx9OhR2O12ZGZmwmAwsKW7fGIQdwlAj8f3AriimPzzlZWZYWHa8dKRMbiN47iqKtVrqfDxbl7J3VepVCI+Ph7FxcVwuVwYGBjA7t27sW7dOpSXl2PHjh1+jyOkys4PNgH4mtj9Q270QmZ6Qgh6e3uh0Wg8Gld465Dj7zh9fX0wGAyiDYtr9PyAHYXQgF10dDRyc3ORm5vL8us7OzsxNTWFlJQUWK1WxMTESFL00dfXh/HxcVRXV0umqxcZGQmz2czyAiYnJzE2Nsaq4FAvICIiYk6JsDdiEF06+frOoqOj8dCWSljRildatCjKlmE10bFdbmlGQEqjDxbcQLFcLseWLVvwxz/+EUePHoVWq/Wz9wxEVNkRAO8wDEMA/C8h5JnZ99MpDVdIlR7Fgsz08xm9y+VCY2MjlErlnPU7v0POfHC73WhqagLDMMjIyBBtAJRQ1NbWNidgFwzhRi6Xw2KxQCaT4ayzzmJnUaqDx2XWiQENLjqdTlRVVUlmDG63G42NjYiJiUFhYSEYhvFwa2lKsKGhAW63m/ViqBYenxhEl0I2m439v7dYgIxh8OjGEkyY7fjNR8P4y7ZSrFtTzKruNDc3sxV2k5OTHtp7xwP8QB53QOKmgyWqsjuDEDI8a9TvMgzTSgj5OMBLP75rerPZjLq6Op/UXKHxAJvNxpJlcnNz0dnZKboW3+l0ssbJD9gFavBu90xvenpMmUwGlUqFhIQEFBUVwWq1QqvVssw6fpmtL9CBMjY2FsuW+V4Di4XD4UB9fT3S0tK88iG4wS0uXXdgYICl63pLCVqtVrS2tiInJ2deYpBKIcOfrijFd56vw/adzfjXtZUoXRKHmJgY5ObmYnJyEj09PRgdHfUgItGqxIUE3+gtFovX4jEpquxmabkghGgYhtkDYC2AjwGM0aKb+fbnY0Hce29ren7jSV/7+jPeqakpNDY2evSRF1uLTwcfriQVEBzDjhpQSkqKT0HOyMhID2adXq/HyMgIW/WWmprKGhAFTcktWbLEo3w3WFBtg7y8PMEBRj5dl5sSlMlkbDqwo6MDy5YtY38fb8QgYOZ3i1HJ8NRVZbj2uVrc+nIjnr+uCrnqGWOiegElJSUeKTUxxBqpcv385qpGo1F0jl5glV0MANmsak4MgAsBPDT78euz+z3ma39vWBCjt9ls7Gu6ftdqtR7rd1/7zme8vuSoxFBquQG75uZjCkTBMOxoo87CwkKkpQlaZkEul3u40QaDAVqt1q8BSQGa6lu2bJlgqXE++ClBm82G4eFh1NfXQ6VSQaPRsIbJTwnyiUEJEQz+unUlrv93A25+qQE7vl2FlFiVhwsttFSYP2hKlfbjy18HQsEVWGWXDmDP7DUrALxACHlr9hCPAXiFYZjvcvf3h5C79wqFglV+oRVkKpUKNTU1fo3J14w9H1kGEC61NTg4iMHBwTnkn0AZdsBM6qylpQWlpaWiG2JScMtsqQENDAygvr4eERER0Ol0YBhGEiHMqakpNDc3o6ysTFJiic1mw+joKNauXYvo6Og5KUEay6AuMZ8YlKeW4fHLl+Oml5tx60sN+Pv/lLFLAm/wVSpMqeBqtRopKSmIjIwMmYCGWKMXWGXXDaDS2/6EkHEA54k6KRYweu9v/T7fvlxQtzk+Pt5j7c2FP/eeBuwsFsucuv5gDH50dJTVxpdSlJIy29atWweVSsWSRNra2hATE8NG08Wua6mYZlVVlaRS2nq9Hu3t7R7H9ZYSbGlpgd1uZw2SnxJcnZ+M321ZgdtfbcKP97Ti3jNn9vcVDKTwVSo8MDCA6elpuFwuaDSaOZ15xYCfVTKZTCdEowtggaL3JpMJR48eFVVaC8w1eqpQ66tjLQVX5JIP2keeP2gQMtNxVKvVQq1Wi5oNCCHo6enB1NSUpKkzYKZjkFarRXV1NfuAcgtsjEYjdDod6urq2M9SU1P99sYbGRnB4OAgqqurJQ2CaTQa9PT0zFH74SIqKopNadJYBjclSGMZEREROHd5Gh68zIWfvd6KR21G/HFblYeEuL/6AMAz9mA2m9Ha2gqj0YiBgQEwDMMuA8T0E+TP9CdKowsgxO49IQQjIyOYmprCGWecIVqMkWv0tHS0oqLCrxvqy72n3kZ+fr5HvTpdUy5fvhyjo6Po7OxETEwM64LONxu43TMtsxQKhU/PIxAQQtDR0QG73e6Tm8+d0QoKClhqbXd3N0wm05w6e4re3l5MTEygurpa0m45w8PDGBoa8hig/IEfyzCZTNBqtR4pwdVJCnyzSIHdnU787QsN7rqgyGt9AOBfK4AQgoiICBQWFgKYCYzSVmJi+gnyC7pOlLJaIIQzPV2/y+VyJCQkiDZ44JjR9/T0QKfTCa7J9+beU7UdvrfBjdAnJCQgMTGRnUEpK00mk7EPJteFo+XDaWlpQUl28UFrD6Kjo1FaWip4IOFTaycmJqDT6dDR0YGoqCikpKRgenoabrcblZWVkua5abVgMAMJNzhHm2R2d3djYGAAF+cpYXJH4Z8HB5AcrcB3zsgXpBjE9wL4M7RKpWIFN/iyW/OVCkuxpj9eCInR2+12HD58GLm5uUhLS0NtbW3AxxocHERSUpJghVpgbvReSMCOe2zuDFpYWAibzTYnnx4XF4eenh4UFRXNqc0PBrSrTkZGhuDYhzdQnTiqeW80GtHU1ASHw4GIiAj09PQgNTXVa2NGMSCEsJ6FlCQhYCbmYDAYcMYZZ0ChUGDFyik0PdeAZz/rRXmEjl3mxMTEzEsM4moF0CWBN3iT3fJVKixECXexIiRGr1KpUFFRgZiYGHbtJRZWqxU9PT2Ijo7GypUrRe1LjZ4G7GgfOvojiSXcREREeHSU6evrQ2trK5RKJTQaDQghgqrr/IGKhSxdulTSgcTlcqGjo4MVEXU4HKzwhdFonFNhJxT0+3W73SgvL5dUempoaAijo6Ooqqpi3WwLE4FhgxM3npGH8vIsNhtAKcO09n2+lKDFYgHDMB7xAF+Yr1TYbDZjZGSEZVOaTCZRBVS0wq63txednZ3vwnsfuxLMVNFRFAK4nxDyJybAPnZACNf0dNTzVbI6H2juPDs726Ozq1DQEZ2qyFZWVkrCsANm+NI6nY6NpE9NTUGr1aK7uxsRERFsHEBs9J6mzoJJ9XkDlQfLzc1lg59KpdLDpaX30NXVJfgeaCwjIiLCa5+/YNDf34/x8XFUVVV5DEIvfzkMBgy2rl6CiIiIOZr43HvgVznK5XJWIKOsrAwARMUC+KXCX3zxBeRyOTo7O/H666/j8OHDcLvdHqW284FW2N1zzz1gGOZ9eKmQI4S0AagCAIZh5ACGAOzhbPJHIrKPHRDCNT1fHFMoqCteXV0Nq9WKkZER/zvxYLfbodVqUVpa6hHlD4ZhRwhBV1cXTCaTR5kpdQeLi2d44lqtFk1NTXC5XGwk3Z+AJB00pE6dUZKQt/p6Cv7DTCsEm5qa2A48qampHiw3l8vloQcgJXp6ejA9PT0n5mB1uLDz6AjOW56CzARPo+JX7XGrHO12O0vTHR4eRnV1NRtf8kYMEqIVQM9JB52ioiJs374d9fX1WL9+PXbs2MEOLL4QQIXdeQC6CCF98x5YAEKeshMKqiRjs9lYV9zhcIjm0NOAXVxcnGQGTwNrkZGRqKiomLdaLC8vD3l5eawLTUUqkpKSkJqaiqSkJI+HaWBgAGNjY6Ii3kJAhTbFeg7eKgQptz4+Ph5JSUkYHh5GZmampDRgOqharVaUl5fPMbgDTRpMWhz4nzX+z8m9B5fLhd7eXlbnoK2tbQ6vwZtikNAmIvR8kZGR+NGPfoR169YJmuzE9rEDcBWAF3nvie5jBywSo6d88uTkZCxfvtxDoVYMh57q2FdWVqKjo4N9nzuSiw000WvLzMwUFVjju9BcxVpKqJmenobdbg9IOHQ+jI+Po6OjIyChTS74HXjo7CmTyTA2NgaXy4XU1NSgvRPKsHS5XF6zFYQQ/OfwIIpSY7AmL1HUsSkTcP369VAqlR68BhqLoQFNb1oBvpqI8MGN3tPPJexjpwKwEcBPOW8H1McOWCD3nqrSeHuwDQYDGhoaUFRUNIenLrTKju8lcKO3wazfjUYjGhsb53WPhYAfSTcYDGwkPTo6GgMDA3PSgYGCsgKlJt1YrVZ0dXWhtLQUycnJrPY+14WmrDqxy6aWlhbI5XK2uQcf9UPTaB4x4v4N4ioKaaBv1apV7HfB5TVQb6y/v5/1ZGiVIL+LEPdZomXC3GCgNyVcKSrsZnExgCOE07uOBNjHDligmd6XGMbY2Bi6urpQUVHhNccpxOgpLTchIYH1EuggE4zB09myrKxM0vyr0+lER0cHsrOzkZOT4zUdmJqaioSEBNHXzM2VS8kKpN2JuMIkfO398fFxtkKQX2LrCzQYGBkZiaVLl/q83xcODyE2Qo7LKoRHx7VaLcsM9DX4cb0xruQZN0dPB2OaEqSqTJTcQ4OB/f39ogZZIRV2HFwNnmvPBNjHDlhgo6cPAF2/TU5OzqtQ68/ozWYzamtrvdJybTYbpqenRXVhoRgcHMTIyMi8VNJAQPXxCgoKWK+Gnw7kdrGNj49nKanzpdKo4q/VapU8V05jA+Xl5T4HP7lc7lHswi2x5eryc/PYbrcbDQ0NbO8+X9AZ7XirWYNtq7MQoxL2uHINXmichBYw0Ry91WpliU0WiwVJSUlITExEX18fCgsL2ZSq2+3Gs88+C5lMJmpyoBV2//jHPwDgAnivsAPDMNGzn/P71AXUxw4AGD9Bh4CLj6k0NAA29xwbG8v2b4+KivKb6iGE4ODBg1i/fv2cz6ieni+GnUajwejoKMxmM9RqNVJTU/1WpVHjsVgsKC0tlZSiSo1HqIwXnXm0Wi3Gx8d9ptLobKlUKiUV1ACOFc5UVlYGvG6nxqPT6WC1Wtl8en9/P6s1MB+e/m8v/vxhD/bdehoKUvwvfzQaDXp7e0UZvD+43W7odDq0trayrMHExES43W4cPHgQr7zyCl5//fVglmfS/WgCsKAzPeW+8/u3+4KvB3hgYADDw8NzGHbcGngagHK73dDr9RgdHUVbWxvi4uKQlpY2Z/akajQxMTGSE01oNZuYwBp35vGVDlSr1eju7oZarZY8dSakcEYIIiMjPTwZnU7H6hYoFAoolUqfyjdOtxsvfzmMMwqTBBn82NgYW+UoZSbE7Xajv78fJSUlSE9Ph8lkQktLC370ox+hv78fN9xwAzo7O1FRUSHZOUOJBTP6iYkJDA0NzauU4w80YGe321FTUyOIYUcFKGhVGpdMExkZyeagW1tbPRpESIWhoSGMjIwEHVjjpwNHR0dRW1sLmUyG6OhojI+Pz0kHBopACmeEgBpPcXExMjIy2PoGXxWCH7TpMGaw4f4Ny/wemxp8VVWVpNfscrlQW1uL7OxslnEXExODgYEBxMfHo6mpCYcOHcLo6OgJY/Qhc+9pswlCCL766itYrVbU1NSIZqp99tlnWL9+PctJT0pKYgUbgeAYdiaTCYODgxgaGkJUVBSWLFkiWRSdS+YpKyuTdKlA6bpFRUVQq9VsOnBiYkJwdaAv0GBgZWWlpNdMmYH5+fle1YRohaBWq2VptQ9/Og2N0YW3fng65DLfv+vo6CgGBwc9KLtSgBr8kiVLPKoy33zzTfzpT3/Cvn37gmqVzsHJ497T9abdbkdBQUHAwhK0cQU/YBdsSs5isWBiYgKnnXYaFAoFG0WnKSg+E00ouOvs+cg8gcBgMKCxsdEjNsAvrPFXHegNoSycocKlS5cuZfX/+eBXCH7VOYKjQyO4sliBxoZ6D9ltLkZGRjA0NBQSg6f8DK7Bv/322/jDH/4gpcEvOEI205vNZnz11VfsqK5QKAKqGvv4448hl8tRXl7uwSwLhmEHHGPCVVRUzHG7KRNNq9WyTR7S0tIEuc80hZiamippuS1wLLBWXl4uqKKLpgO1Wu286UBu4YyvXHmgsFqtqK2tFa2/9+C+NrxWN4oPbl8HJTnmBRBC2GXA9PQ0W5QjpVdCDT49Pd0j9vTBBx/gwQcfxP79+yUtiMICz/QhM3qNRgOHw4GUlBQMDQ3B4XCIDjb19/ejvb0d69ev95ipghGtpOwvh8OBlStX+t2fFnNoNBpMTEywKrUpKSlzZhar1cqKdEjZsgo41sWmsrIyoMAaTQdqtVpMT0+z6cCkpCS0tbWx7ZulNHjK/V++fLmoOI7B6sTX/vgZvlGaikc2rvD4jEuoMZlMyMjICKhC0Bfcbjfq6uqQmprqMUl9/PHHuO+++7Bv3755VZsCxMnh3qekpLCSVXK5HFarVfC+breb7REWHx8vSdMJ4JiwR1xcnOD0FreYg6tS29fXB6VSybrPDocDjY2NWLFiRcCBSl8YGBiARqMJKirNz6VPTU1Bo9GgqakJKpUKeXl5sNlskmn7GY1GNDQ0BFQ1+FrdKCwOF65ZM9czpE00lUolzjrrLBgMBpZ5F0yVIzDz3FHZcq7Bf/bZZ7j33nvxxhtvhMLgFxwLmrITAm7AbsWKFTh69KgklFqr1Yr6+nrk5OR4rNHEgK9Sy+35ZjabkZ2dDYVCIZnMMg0Gms1mSdtZ01xzV1cXiouLoVarA6oO9AUad5iP0OMLbkLw4peDqMqOx8rMubJog4OD0Gg0bKCRW11Hu+/QCkEaBxASl6FkIbVa7dHo4/Dhw7jrrrvwxhtvSFpgdDyxIEYvtIklDdgtXbqUdY9pbXwwBk+JMVLPwlFRUVAoFFAoFDj99NMxNTWFrq4uWCwWqNVqpKWlBUSnBY55OzKZTHLeAI2kcxtbiK0O9IWpqSm2eWggWZCD3RPoHbfgN1vmsvQGBgbmzSzExMQgJiaGvQ9u9x26nFGr1XOWZbSVV2Jiokcc5ujRo7j99tuxZ88erx1/TlSEtOCGQshMTxl2/IAdVbYN1OA1Gg26u7uDrjjjg0a7DQYDqwsXExODJUuWsAqvlE4rVpnG5XKhoaEBCQkJkvavB46l+3wVEfmrDpwvHTgxMYG2tragGHwvHB5EcowKF670DJT19/dDr9cL1vbjd9/h8jO4uvuRkZFoampCXFwc8vLy2P0bGhpwyy23YNeuXfPShE9EhCyQRwhhVW/MZjOrg+4N/f39GBkZQVVVlUeQihCCwcFB9PX1ISEhAWlpaYLlqQkhbM65oqJCcpIJrQwTQiWmqi56vR5RUVFsHMDbNdHljdT16sCxTjZiA2sAPNKBOp1uTjqQrqv5v6EYDE5YcNETn+Pms/Lww68Xsu/39fVhYmICFRUVkixx6LJMp9NhenoaMTExKC4uZr2y5uZm3HDDDXjllVewfPnyoM8nACdH9J5r9DabDY2NjVi9erXHNtSFdTqdc7juXNFKAGwEXa/XIzY2FmlpaUhJSfE6c1LmHk1BSZlzdjqdqK+vh1qtRl5enugyUirxrNVqPeSfo6Ki2Og/t6BDKggpnBEDbjrQZDKx33VKSkrAnsnv3u3Ec58P4r3b1yE9fmbg6O3txdTUlFdhjWBACEFzczPbUFSr1eK///0vDhw4gO7ubvznP//B6aefLtn5/ODkMHoAbA87p9OJr776Cqeddhr7mcPhQG1tLZKTk1FQUCCYYUcj6BqNBjqdDpGRkUhLS2NnTlrQk5iYKLlrTIOBXL25YI/HzaPbbDYsW7aM1lhLcMUzoPn9iooKybuw0Pr9nJwc6PV6j3Sgt/WzL1gdLnz9T5/h9IIk/PGKGampnp4eGAwGlJWVSW7wLS0tUKlUHiW97e3tuPXWW1FWVob6+nps374dV111lWTnnQcnR8qOC1+davjCGUIi9NwIelFREUwmEzQaDY4ePQqGYWC1WlFQUBCUfLQ3UEGNkpISyZhYtGttXFwcmpubkZ+fz9ZzJyUlIS0tLeh+dXQdK3WZMHBMsZbW73Pr0sWKhe5v0mDK4sT/1Mwsabq7u2E0GkNi8FTJmGvwfX19uO666/CPf/wDNTU17LaB4q233sLtt98Ol8uF733ve7RunsVHH32ETZs2oaCgAHV1dbUAdhNCHvJ6MIkRUqOn6jlc46UtqvmdagJl2MXExKCgoABqtRqNjY1IT0/H2NgYK0+clpYWtB65WCacGNBqturqatYoaACNtnqKi4tj6+rFUE25raukjGkAvhVrhVQH8tOBhBD859AgitNiUJOXyKYpQ2HwbW1tkMlkHkSkwcFBXH311fjf//1f1uDpvQQCl8uF2267De+++y6ys7OxZs0abNy4cY6U+1lnnYU333wTmFW8XSgsqEZeX18fRkdH57So5hp8ID/y2NgYent7UV1dzUaNaQFHR0cHS0FNS0sT3dxhZGQEAwMDIZkpBwcH2ZmSa5R8ea3p6WlWGILOnKmpqfNejxQdZ3zBl2KtNwgRC+0zydAyasQvNixjxTHLysokXeJQJiYAD2LWyMgItm3bhj//+c+SreEPHTqEoqIiVl3nqquuwt69e0X3bwgVFsToaZOBqakprFmzZk6bIW9dZoSAEIK+vj7o9fo5hsMt4KBcetrcgebQ59NzI4Sgt7cXk5OTHpLXUoA2vJyensaqVavmPTbDHOv7TpczWq2WbcFMBwDqgYSycIaShSwWS0CBNV/pwCc+HkO0gkEeNLBYFCEx+M7OTrjdbg/h1bGxMVx55ZX4/e9/j7PPPluy8w0NDXnk9bOzs/HFF1/M2e7gwYOorKxEfX39AQB3EUKaJLuIeRBy995ms6Gurg4ymczDXQuWYUfTZjKZzO/DzVV1paIaVM/NWyqQZhUYhpG85xu3uCWQFBQloOTn57P6/u3t7aw3YzaboVAoJCf0cBVrpTBK6s24I2LxlWYUFxVGgnHZYbE48NVXXwmuDhRy3V1dXXA4HB7FRDqdDldeeSV+9atf4dxzzw3qHN7OyQf/+6qurkZfXx9d5jwB4DXw+tSHCiE1eqPRiNraWhQVFaG7u5v9MoI1eFrJRuWWxOzPF9WgqcCOjg5Wmnp0dJRt5CCl4VB1ntjYWA9NgEChUqnYhgs0v2+328EwDFpbW9kIerCDFg1+yWQyyavwXv1qGE43wTeWRmPt2nJ2opBKLLS7uxs2mw0rV65k99Xr9bjyyivx4IMP4qKLLpLsXiiys7MxMDDAvh4cHJwjzsKrGN3PMMxfGYZJIYToJL8gHkKasmtqakJ6ejri4uLw5ZdfsiSZYEpiaeVWYWGhVzGGQEEIgV6vR1NTE2QyGWJiYthUoBRS0nSgSktLk5zSSRl8NE3JbfOk1+uDEtYQqlgbCOxOF8794yfIS1Dg3zeu93psX9WBQtKB3d3dMJvNHlr6k5OTuPzyy3H33Xdjy5Ytkt0LF06nE8uWLcP777+PrKwsrFmzBi+88AJKS0vZbUZHR5Genk5t4DQAOwHkkWBSBgIR0pm+pKTEo0uIw+Fge9sFYvCTk5NoaWmRvN8bMMNWo5LXarWaXTvTpUlqairS0tICqt6ipBuuCq5UcDqdbO03TVPyKwONRiOb1lQoFKzr7O9eKCed6sRLCUII/v1hPfQWN365yTer0Vt1oJB0YG9vL6taRI89PT2NrVu34o477giZwQMzy8knn3wSF110EVwuF2644QaUlpbi6aefBgDcfPPN2LlzJ5566ik6cP0ZwFULYfBAiGd6viJuRkYGkpOTA5otKAmkoqJCsvJPCsoZ96VxT0k0Go2G7eoiNBVIqa9S5vcpvBXO+AOloGq1WjaFRu+F+7vQfnWUeSgl6HLh/o/0mHLIcOAH88th+QLtWce9l9TUVIyPj2N6etojhmQ0GrF161bceOONuOaaayS9Hwlw8jDynE4nWyGn1+vR19cHu92OlJQUpKenz3nQvF7AbKSbUjGllEQCxA8mDoeDHQCsViv7oHkr36SdaMvKyuZ0PwkWVJEmmO47NIWm0WhYqXCa1uR7D1KBsuEGjQTb39LgrvOX4ob1wSsMcdtvWywWZGRksBkam82Gbdu24dprr8V3vvMdCe5CcpxcjDwasEtKSoJarYbT6WRbClssFnam8ZY/p+tJhUKBqqoqyaPR3HSf0MFEqVSyqUAq6UzLN7ksOr1ej87OTsk70QLHOs4EWyrMTaFxKwPHxsYQHx8PlUoFl8slWbqSy3f/TONChEKGb64KTNuAD0rBjoyMRE1NDbsMuPrqq6HVanHmmWdi48aNkpzrREdIZ/q77roLmZmZ2LRpk1c+OTWasbExmEwmlkCTkJDgEfiSWmuOupdSFuTQvLNGo2G13JYtW4a0tLSQNM2QqnCGC24v+8jISLbRBq1v4HZ6FQtq8BEREUjNysW5fzqIDWXpePgyaarYBgcHodVqPVKsNpsN11xzDaqqqhAdHY3PP/8cb7zxhqSTh0Q4edz73t5e7Nq1C6+99hoA4NJLL8XmzZuRnZ3tdQDQ6/UYGxvD1NQUHA4H8vLyQpI247ZTCpX3QHn0er0e0dHRrNEEszwJZeHMfIq1tL5Bp9OBYRjROXRCCJqamhAVFYXCwkI8/8UgHnunEztvrPGqjiMW1DvhimvY7XZcf/31OPvss3HHHXdI8jv749NTHD58GKeffjpefvllXHHFFUIOffIYPXsQQjA8PIxdu3Zhz549sFqtuPTSS9mCA+4Potfr0dbWhuzsbBgMBkxNTSEhIQHp6elBN3OgbadD0dSCklecTqeH98CNnut0OlbAIS0tTdSsSSPWwdSr+4IYxVqaQ6fCp/7ozW63G01NTYiOjsbSpUvhJgQb/vIFUmJU+Pd3qoO+9uHhYVaLgRq80+nEDTfcgJqaGvzkJz+RxOBdLheWLVvmwad/8cUX51BrXS4XLrjgAkRGRuKGG244dY3e44CEQKPRYM+ePdi1axempqawYcMGbNq0CZ9++iny8/NxxhlnsA82X402Pj6ebUslZgCgUfRg2057A01tRUVF+VWUNZvN7BKAYRiWCzDfup8WzkjdvYVeTyCKtcAxqXCNRgOj0ThHWstbyu/jjnHc/GI9fvfNldhQFpxi8MjICIaHhz0M3uVy4aabbkJJSQnuv/9+yTy5gwcP4oEHHsDbb78NAPjVr34FAPjpT3/qsd2f/vQnKJVKHD58GJdeeumiNPoFLbgBZuiI6enpuPnmm3HzzTdDp9Phtddew9atW8EwDLZs2YKMjAyW+cXPOU9NTWFsbAydnZ1+xTQoaH4/FFF0KqohpBkjMFOAkp+fj/z8fNhsNmg0GrS0tLBCjvz0WSgLZ2hAMFDeA5/ezJXWio2NhcViYfUSAMDucuN373UhMyEC568ITiRkbGyMbXLBNfgf/vCHKCgokNTgAWF8+qGhIezZswcffPABDh8+LNm5pcaCGz0fKSkpMJvN2LJlC+655x68+eabePjhh9Hf348LLrgAmzdvZjnq3LJNWn1GS1OjoqKQnp4+Z91MK/BWrVoleX6f1hUEKqoRERHB9ninKSea1UhOTobdbofT6ZS8cAYITrHWG7iVgS6XC0ePHoVcLodOp8PU1BTS0tLweqcNnVoT/npVOVTywO+H26iS/tZutxt33HEHUlJS8PDDD0serBPCp9++fTt+/etfSz44S40Fd++9wVtaaHp6Gvv27cOuXbvQ0dGBc889F5s3b8bq1avnGIC3dXN6ejqsVismJydRXl4eMrdYbOcWIaDeg9lshlwul0xQg4LyByoqKiTXB6Da8UlJSSypx2w2o7ZrGLfsHUBVmhwPXJATsM6BRqNBX1/fHIO/++67IZfL8fjjj0s+QALC3PuCggJ2cNDpdIiOjsYzzzyDzZs3+zv8yb2mDwRmsxn79+/H7t270dDQgHPOOQebN2/Gaaed5nVUNRqNaG5uhtlsZmMAYgNn84GmzUJBB6bcBNpxhhDCpgInJydZ7jm/1bZQSKFYO9+119XVITk52WOpQwjBDc/XomnEgN3fWwWZzcCSm8QU01BNAW7TD7fbjZ///OewWCz461//GhKDB4Tx6bm4/vrrw2v6YBAdHY0rrrgCV1xxBaxWK959913s2LED27dvxxlnnIEtW7Zg/fr1UCgUcDgc6OrqQnJyMtasWQOr1YqxsTGWQ08DZ4G6+uPj4+jo6JBcUhuYWzgDzLiQXEENLvc8KiqKjWkI8WTGx8dZwpDUSx1K201JSZlTUPR6/Ri+6J3E/RuWITs5DkAcS24aHx8XJBVOxTe4Bk8IwcMPP4zJyUn8/e9/D5nBA8L49CcKToiZ3hfsdjs++OAD7Ny5E59//jmqqqpQX1+Pxx9/3EOEk8JqtUKj0UCj0YAQwnoAQmc8qqJTVVUlmddA4a1wZj5QZV26pFEoFOyA5i2lp9Fo0NvbG5Jrpw0f09LS5lz7hNmOS/5yCPnJUfj3d6oh8zGbu91uttUWn9swPT3NymvTayeE4LHHHkNvby+effbZRb+O9oOwex8I2tracNlll2HlypXo7OxEdXU1Nm3ahHPPPderEXDzzU6n028RTV9fH8bHx1FRUSE5/z+Qwhk+LBaLBxuQO6CNjo6yg5XUsQ1fHV4p7t3bgjcbxrDr+zUoThMWMOQOaKOjo6zYaUZGBqKiokAIwR/+8Ac0NjbiP//5j+S/x3FA2OgDwTvvvIPU1FSsWrUKLpcLn3zyCXbt2oUPP/wQZWVl2LRpE84//3yvLjktohkbG5tTEASA1dkrLS2V3IWUonCGD+6AZjabAQBlZWUBt9jyBZfLhdraWmRmZnolO33eM4Ebnq/FjWfk4o7zloo+/sTEBNrb27F8+XI2U/PEE0/AYrHAbDbjwIEDknstxwlho5cSbrcbX3zxBXbu3Il3330XxcXF2LJlCy688EKvqSpaEKTRaGCxWOB2uxEXFxcSg5eqcMYXaLfbjIwM6HQ6SXrsUdDliC+Dtzld2PT0YRAC7L15DSKV4tzvyclJtLa2egiSEkLwxBNPsO2iOzo68PHHH0seTD0OCBt9qOB2u3H06FG8+uqreOutt5CXl4dNmzbh4osvRkJCgse2NG2mVCpZd5NbEBTsjEkzAKEgDAHeO8PwVWgSExORlpYmmt7sdDpRW1uLrKwsnx2Af/1OJ577fAB/v7YS6wvFpTRpE0xuwJEQgmeffRZvvvkm9uzZg8jISNjt9qBmen9c+r179+LnP/85ZDIZFAoF/vSnP+HMM88M+HzzIGz0CwFKEX311Vexf/9+pKenY9OmTbjkkktgs9lw9OhRrFq1in2ouQVB3DLapKQk0QMATZuFonCGquFSmShfxsynN8fFxbH05vmCYk6nE0ePHkVOTo5PQtKh3gl8Z0ctttUswf0bSkRdP7frLTfA+vzzz2Pnzp14/fXXJUk1CuHSG41Glh1ZX1+PrVu3orW1Nehze0E4ZbcQkMlkqKioQEVFBR566CG0tLRg586duPTSS6HRaHDNNddgzZo1rLQXt+8cpZyOjo6ira1NVEFQKDvOEELQ0dEBp9PpV7GWT2+ma+b5UoG0FVlubq7PgKPB6sRP97YgVx2Fu84vEnX909PTXg3+5Zdfxosvvoh9+/ZJxi0Qok3PXf6ZTKbFWJIbEE5Zo+eCYRisXLkS119/PXbv3o0nn3wSbW1tuPrqqxEZGYnLLrsMmzZtYoUMuZRT7ozZ3t4+b0FQKDvOUI0AhmFEK9Z609bna+olJSWhpaVlXoMHgEff6oBm2o5/f2cVolXC1/EGg4FlCXINe/fu3fjXv/6Fffv2ScoeFKpNv2fPHvz0pz+FRqPBvn37JDv/8UTY6DnIysrCm2++yeaa77nnHvT19WHXrl349re/DZlMxmoCZGVliSoIGhoaClnhDJ/FF8yMxDAMYmNjWZlui8WC0dFRHDp0CJGRkWzk3Nuy5EDTGPbWj+KWs/JQmZ3g5ejeQfsE8pc7b775Jp566ins27dP8riHEC49AGzZsgVbtmzBxx9/jJ///Od47733JL2O44Gw0XMgl8s9yCUMwyA/Px8//vGPceedd7KaADfddBNsNhurCUCFPvgFQQaDAWNjY+wMTGm1UoLGJqiRSg25XA6tVovy8nK2pXNbWxub2kxLS0NsbCyGp6z4xZttqMyKx81n5ws+vtFoRENDw5w6gLfffht/+MMfsG/fvpBkNoRo03Nx9tlno6urCzqdbo7IyImGUzaQFwyoJsDu3buxe/duTE1N4ZJLLsGmTZtQXFzs0Zixra0NLpcLubm5cwqCaHvtQBFKxVrgGGmooKAAqamepbBOp5MV1Zw2mvD7o24MGpzY/f01yFELC07SlCW/0u/999/Hww8/jH379s05r1QQwqXv7Oxktf6PHDmCyy67DIODg6FY24ej9ycaqCbA7t27odFocPHFF2PDhg144YUXcN11181ZY9M1s1arZemzYguC5qO+SgG73Y6jR496lc/i47fvduBfBwexfW08yhPsggKbZrMZdXV1cwz+448/xs9+9jO8+eabAZUri8H+/fuxfft2lkt/3333eXDpf/3rX2PHjh1QKpWIiorCb3/723DKzhv0ej22bduG3t5e5Ofn45VXXpmj997W1oZt27axr7u7u/HQQw9h+/bteOCBB/C3v/2NHeEfffRRbNiwQexlHDdMTExg165d+MUvfoGsrCx8/etfx+bNm302fLRYLBgbG4NWqxXcVENInjwYUL28oqIivyzBj9p1uPWlBmxdvQQPXFIyJxXoTeiEGjyfo/DZZ5/h7rvvxptvvim5nNkix4lt9HfffTfUajXuuecePPbYY5iYmMCvf/1rn9u7XC5kZWXhiy++QF5eHh544AHExsbirrvuEnvqRYOf/OQnKC4uxtatWz00Ac477zxs3rwZ1dXVXgcAIQVBDocDR48eDYqnPx+owRcXF/vVCRicsOCKv32JJYmRePGGakQoPAOUNK5BlzWRkZFITEzEyMjInLLkQ4cO4Y477sDrr78ueduvEwAnttGXlJTgo48+QmZmJkZGRvC1r30NbW1tPrd/55138OCDD+LTTz8FgJPC6N1u9xyjppoAO3fuRHNzM6sJsHbtWq/RfLvdzg4AtCAoKSkJbW1tKCwsDMlaV5RAptOFa/51BAN6K169sQa5av/5c71ej4aGBqhUKqhUKrZVlU6nw2233Ya9e/eyJcWnGBbU6CUvQB4bG2NdzszMTGg0mnm3f+mll3D11Vd7vPfkk0+ioqICN9xwAyYmJqS+xJDD2yxONQFeeuklHDp0CBdccAGeffZZrFu3DnfeeSc+/vhjOJ1OdnuVSoXs7GxUV1dj1apVkMlkOHLkCJxOJ6anp2E0GiXNBFCDLykp8WvwhBA8tL8dzSNGPLZ5hSCDt1qtaG9vR2VlJdatW4fS0lKYTCZcffXV2Lx5My688EI4HA6pbieMeRCQ0Z9//vkoKyub87d3715Rx7Hb7Xj99ddx5ZVXsu/dcsst6OrqYqu3fvzjHwdyiYsakZGR2LhxI55//nkcOXIEmzZtwquvvor169fjRz/6ET744AMPA3A4HBgZGUF1dTVOO+00REdHo6urC1988QU6OzsxPT0d1ABADX758uWC+u29+OUQ9tTO5OO/XuI/fUWXDFzF3cjISLjdbrhcLrzxxhtYuXIlduzYEfA9ULz11lsoKSlBUVERHnvssTmf/+c//2GZmOvXr0ddXV3Q5zzRcFzd+7179+Ivf/kL3nnnHY/3aTCws7MTWq0WAwMDXh/G/Px8xMXFQS6XQ6FQ4Msvv/TYf75g4mKEw+HA//3f/2Hnzp345JNPsHr1aqxevRoff/wxnnzyyTnVZLRDEJWgDqQgyGKxoK6uTpQE9jstGrzTrMVvvrnSpygGBTX4ZcuWefwGbW1tuO666/DCCy+grKxM0Hn9QQif/rPPPsOKFSuQlJSEAwcO4IEHHvDKxFtgnNhr+v/3//4fkpOT2UCeXq/Hb37zG6/bXnXVVbjooos8mgqOjIzgj3/8I9RqNSIiIvDss8/iG9/4htdgYH5+Pr788ss5KSWxwcTFCJfLheeeew733HMPCgoKsHTpUlYTwBv/PJCCIGrwK1asmFNlKAVo2o8fFOzq6sI111yDHTt2oKqqSrLzCdWmp5iYmEBZWRmGhoYku4YAsbCkfkLIfH+iodPpyLnnnkuKiorIueeeS8bHxwkhhAwNDZGLL76Y3c5kMhG1Wk0mJyc99r/22muJSqUiy5cvJ5dddhk5evQoWbZsmddz5eXlEa1WO+f9ZcuWkeHhYUIIIcPDwz73X+y4+eabSXt7O3G5XOTTTz8ld9xxBykvLyeXX345+fe//000Gg0xmUxz/gwGA+nv7yeHDx8m7733Hjl8+DDp7+8nBoOB3Uar1ZL33nuPDA8Pez1GsH8TExPkww8/JP39/R7vNzc3k4qKCnL48GHJv69XX32VfPe732Vf79ixg9x2220+t//tb3/rsf1xhD87lPRvUZJzEhMTMTk5yb5OSkryGtArKChgZ7KbbroJ3//+90XtfyLC7XbjyJEjePXVV/H2228jPz+f1QTwJiZBOGq6tENQQkIC+vv7UVZWFhIBCppWLCws9PDCBgcHsXXrVjz99NM4/fTTJT8v/U7+/ve/A5gpxz106BCeeOKJOdt++OGHuPXWW/HJJ59I3vEoAJwapbXnn38+RkdH57z/yCOPCD7Gp59+iiVLlkCj0eCCCy7A8uXLcfbZZ0t5mYsOMpkMNTU1qKmpwa9+9Ss0NDRg586duOSSS5CRkcFqAtD1M8MwHgVBtBZAqVSir69PUIcgMaDltwUFBR4GPzIygm3btuGJJ54IicEDwvn09fX1+N73vocDBw4sBoNfePhxBY4LAnHPf/GLX5Df/va3Ae9/osPtdpOmpiby4IMPkjVr1pALL7yQ/OUvfyF9fX2sa63RaMh7771HRkdHidFoJCMjI6S2tpa8//775LPPPiOdnZ1kamoqYJd+amqKfPTRR6Snp8fj/e7ubrJq1Sry/vvvh/Q7cDgcpKCggHR3dxObzUYqKipIY2OjxzZ9fX1k6dKl5NNPPw3ptYjEgrr3oRMKDwIbN27Ec889BwB47rnnsGnTpjnbmEwmGAwG9v/vvPMOGwWm++v1epx99tmsJ+DNxR8YGMDXv/51rFixAqWlpXj88cfZzx544AFkZWWhqqoKVVVV2L9/fyhuVxJQTYD7778fX3zxBZ588klMTU1h27ZtuPTSS/HII4/gyiuvRGlpKdthNj4+HsXFxTjttNNQWFgIs9mMr776CkePHsXw8LCovDmlBufm5iItLY19X6fT4corr8Rjjz2Gc889NxS3zoKrTb9ixQps3bqV1aannPqHHnoI4+PjuPXWW1FVVYWampqQXtNixKJc04+Pj2Pr1q3o7+9Hbm4uXn31VajVagwPD+N73/se9u/fj+7ubmzZsgXAzAP3P//zP7jvvvs89j9y5AjUajUOHz6MZ555xmsUf2RkhM2BGwwGrF69Gq+99hpWrlx5UrADCSF4++238d3vfhdFRTNKNlQTYMmSJV4j+2ILgmjvuuzsbI8iGb1ej8svvxz3338/LrnkktDc4MmBEztlt5gglhIMAJs2bcIPfvADXHDBBSeF0QPAww8/jK1bt2LZsmUYGhrCrl27sGfPHtjtdlYVKC8vz+sA4K8giMpgL1myxKP4Z3JyEpdffjl+8pOfCOnldqojbPRSQWwUv7e3F2effTYaGxsRHx+PBx54AM8++yzi4+NRU1OD3//+9ycEyUcIyGxQj2oCGAwGVhPAl/oOvyAoJSUFOp0OWVlZHgGz6elpXHHFFbj99ts92JZh+ETY6MVgvizAddddJ9jojUYjzjnnHNx333345je/CWCmjiAlJQUMw+DnP/85RkZG8M9//jMk93G8QTUBdu3aBa1Wiw0bNmDjxo0+9fYsFguOHj0KAKyOXkREBOLj47F161bceOONuOaaaxb6Nk5UnNjknMUEoVF8u91OLrzwQvL73//e57F6enpIaWlpSK5zsUGv15Nnn32WXHbZZWTVqlXkJz/5CTl48CBL7pmeniaffPIJaWtrIyaTiUxOTpKOjg6yefNmsmTJErJp0ybS2NhI3G738b6VEwXh6L1UEJIFIITgu9/9LlasWIE777yTff+tt97C0qVL2cKNPXv2eHDECSH40Y9+hKKiIlRUVODIkSMe+85X9LHYkZSUhOuuuw6vv/46PvroI1RUVLCqMffeey8uvvhiGI1GVrFHqVRCrVbDZDLh7rvvxtatW/HAAw9Ar9cHfS3+vsvW1lasW7cOERER+N3vfhf0+U4J+BkVTmgIoQT/97//JQBIeXk5qaysJJWVleT1118nhYWFZNOmTWTlypUkMjKSfO1rX2O9BkII2bdvH/nGN75B3G43OXjwIFm7di0hhBCn00kKCwtJV1cXmytuampa+JsPAfR6PVm7di1Zt24dKS8vJz/84Q/J+++/T3Q6HdmwYQN54oknJJ3dhXyXY2Nj5NChQ+Tee+9leRonIBZ0pj+p1XCTk5Px/vvvz3l/yZIlbM79zDPPBOHFNQ4ePIiioiK89tprAI4VbnCj03v37sW3v/1tMAyD008/HZOTkxgZGUFvb6/fJgonKtrb27Ft2zbceeedsFqteOedd/Cvf/0LBw4cwC233ILbbrtNUtFIIQ0paDrxZNGkXwic1EYfKIQ0QvC2zdDQkOAmCiciTjvtNJx22mkAjmkCbNy4EVNTU4iPj5dcJfZk/i6PJ8JG7wX8mR+Y2wjB1zZC9j3ZEIqyXEB4Q4owxCFs9F4gpHDD1zZ2u11UE4UwfENsQ4owhOGkjt4HijVr1qCjowM9PT2w2+146aWXsHHjRo9tNm7ciB07doAQgs8//xwJCQnIzMycs+/f/vY3PPPMMwHJN+Xn56O8vPyU5YgL+R3CCAB+In2nLPbt20eKi4tJYWEh+eUvf0kIIeSpp54iTz31FCFkpqrt1ltvJYWFhaSsrMxDFILuW1BQQJKSkuaNPn/66adEr9cTQgjZv38/mwUgxLdIyKkEf7/DyMgIycrKInFxcSQhIYFkZWWRqamp43nJgWBBo/cnPCNvMSNY+SZfcmBhnHQ4sSWwwzgGXxF+X/jHP/6Biy++mH3NMAwuvPBCrF69Gs8880xIrzWMUwfhQF4I4c2L8hV9/vDDD/GPf/wDn3zyCfveqagMFEboEZ7pQwix8k179+71kG+i26alpWHLli04dOhQ6C86jJMffhb9YQSBYOSbjEYjmZ6eZv+/fPlykpWVRZYuXUp+9atfzTnXhx9+SOLj41kq8YMPPsh+duDAAbJs2TKf+4Zx3LGggbyw0YcY/qLP3/3ud0liYiJrrKtXryaEENLV1UUqKipIRUUFWbFihd8swIcffkguueSSOedfTLUA/gYft9tNfvjDH5KlS5eS8vJy8tVXXx2HqzwuCBt9GJ747LPPyIUXXsi+fvTRR8mjjz7qsY0voxey70JAyODjq4jpFMCCGn14TX8CQGgW4ODBg6isrMTFF1+MpqYmUfuGGtziGZVKxRbPcOGriCkMibHQo0z4T/wfgCsB/J3z+lsAnuBtEw8gdvb/GwB0CN13ge7hCi/X8SRvmzcBnMl5/T6AmuP9/Z9sf+GZ/sTAIIAczutsAMPcDQgh04QQ4+z/9wNQMgyTImTfBYK3XCU/pylkmzCCRNjoTwwcBlDMMEwBwzAqAFcBeJ27AcMwGcwsCYBhmLWY+W3Hvex7I4DvMwzTyTDMPfwTMQzz/xiGqZ39a2QYxsUwjHr2s16GYRpmP/tS5D0IGXwWywB1cuN4uxrhP2F/mHHZ2wF0Abhv9r2bAdw8+/8fAGgCUAfgcwDrfew7DqAQgGp225XznPMyAB9wXvcCSAnw+hUAugEUcM5dytvmEgAHMDPjnw7g0PH+3k/GP3/c+zBOIjAMsw7AA4SQi2Zf/xQACCG/8rH9CwA+JIT8bfZ1L2bW2LoAz78BwJ8AyAH8kxDyCMMwN89ew9OznsqTAL4BwAzgO4QQsR5FGH4QNvpTCAzDXAHgG4SQ782+/haA0wghP/CybTRm3O0iQoh+9r0eABOYWWf/LyEkXBBwAiLMvT+1ICZQdhmAT6nBz+IMQsgwwzBpAN5lGKaVEPKx5FcZRkgRDuSdWhATKLsKwIvcNwghw7P/agDsAbA2BNcYRogRNvpTC36zAADAMEwCgHMA7OW8F8MwTBz9P4ALATQuyFWHISnC7v0pBEKIk2GYHwB4G8eCaU3cYNrsplsAvEMIMXF2TwewZzYrqADwAiHkrYW7+jCkQjiQF0YYpxjC7n0YYZxiCBt9GGGcYggbfRhhnGIIG30YYZxiCBt9GGGcYggbfRhhnGIIG30YYZxiCBt9GGGcYvj/yM7t2j618qgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.plot(cloud[:,0], cloud[:,1], cloud[:,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<mpl_toolkits.mplot3d.art3d.Line3D at 0x7f1ea942a640>]"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAADzCAYAAABNPxCMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABzEElEQVR4nO19d3gc1dX+O9vU26pbXZYs26qWZYNNS6jBgEsAGz5IICSElgRD+BECCaEEQnoIJPCRBiahumDANh0+AhhssNV779pdrcr2en9/SHc8O9rVzuzOyrK97/PosXd36u6ce8895z3vYQghCCOMME4dyI73BYQRRhgLi7DRhxHGKYaw0YcRximGsNGHEcYphrDRhxHGKYaw0YcRxikGhZ/Pw/m8MMIIPZiFPFl4pg8jjFMMYaMPI4xTDGGjDyOMUwxhow8jjFMMYaMPI4xTDGGjDyOMUwxhow8jjFMMYaMPI4xTDGGjDyOMUwxhow8jjFMMYaMPI4xTDGGjDyOMUwxhow8jjFMMYaMPI4xTDGGjDyOMUwz+6unDCAHcbjesVitkMhkUCgXkcjkYZkFLqsM4hRE2+gUEIQQulwsOhwMOhwNut5s1doVCwf6FB4EwQgnGT7OLsHKORCCEYHx8HNHR0ZDJZHA4HB6fEULmDAJKpRIKhQIymSw8CJzcCCvnnGxwu92w2WxoaGjw+jnDMB6uvkwmg8vlgsViwfT0NPr6+jAxMQGbzQaXy4VwV6IwgkHYvQ8hCCFwOp1wOp2sYXM/8zV7Mwzj8ZlGo4FCoWDfYxjGYzkQ9gTCEIOw0YcIbrfbY91O/wKZpem+crkcwLHBhC4RwoNAGGIQNnqJwQ3WAXNnbSnAPyYhhA0O0s9pTIAuF8KDQBgUYaOXEHx3nm9oMpnMY6afz8Xnwt82XC+AHpc/CLhcLsTGxoYHgTDCgTyp4Ha7YbfbfRo8BTV0sUYnZllABwH6xzAMjh49CpPJBIPBgOnpaZjNZtjtdrjd7nBg8BRDeKYPElx3nh+s4yOYNX0woIOMQjHzcxNCYLfbYbPZ2M+USiW7HAjFkiSMxYOw0QcBajzcYN18CNTo6bmkAj8oCAB2ux12ux0A2PQhNyYQxsmDsNEHCOrOi3HXgzH6UIFed3gQOHUQ/vVEggbJ6urq4HA4RAXFgnHvF2qw4KYHuQZut9vx5ZdfYmpqCtPT07BYLGxKMowTC+GZXgS47rzZbBb9wFPjJYRAq9UiMjISsbGxi3r9zPUELBYLm4Gw2Wyw2WwAZjwBShkOewKLH2GjFwg+GYaffhMChmFgt9vR2toKhUIBl8sFk8mEmJgYJCUlISkpCVFRUYt6EOAbNB3EuIOAXC5nlwJcJmEYiwNho/cDX7n3QFxuh8OBxsZGlJSUIDExkd3fZDJhYmICnZ2dsFgsiIuLYweByMjIRRkLoPBGFKKlwzTeQQeBcAXh4kDY6OeBNyothUwmE+zeE0LQ09OD6elplJWVIS0tjQ2SMQyD2NhYxMbGIicnB4QQGAwGTExMoLW11SOYFhsbC5VKJf2NSoj5BgEAmJqaQlRUFOLj48ODwHFC2Oi9gE+l9bZGFTr72u12NDQ0IDY2FqmpqYiIiJh3e4ZhEB8fj/j4eOTl5cHtdqOlpQUWiwWNjY1wuVxITExEUlISEhMT2dz7YgV/EJiYmIDL5fIYvMKewMJicT8xxwE0Ou9yueZNxQlZ0+v1erS0tKC4uBhpaWloamoS7abLZDJERkYiMTERycnJcLlcmJycxMTEBHp7e8EwDBITE6FWqxEfH+9Bx12MIISwmQH62u12w2KxeAQNw4NA6BA2eg7E5N4ZhvHp3hNC0N3dDZ1Oh9WrVyMyMpLdhxp9oOt0uVyO5ORkJCcnA5iJE0xOTkKr1aKzsxMKhYKNB8TFxS26SDq/3oB+z/Q6w4NA6BE2esxf9+4LvoyWimXEx8djzZo1khjdfAOEUqlEamoqUlNT2fNPTExgeHgYBoMBERERSEpKYsU3jrfB+LsGIYMAlygUHgTE45Q3eqpqA4grg/U204+Pj6O1tRXLli1jjZC/T6ij8BEREcjIyEBGRgYAwGKxYGJiAna7HYcOHTru6UGxA4+3QcDlcuHQoUOorq4GEJYWE4tT2uhpnryurg5r1qwR9bBw1/SEEHR1dUGv13u483wcD0ZeVFQUoqKiMDQ0hJqaGpjN5nnTg6GG2+0Oyvvh/kZyuZwdBJxOJ/t5WFBkfpySRs9156nxin0wqCHabDbU19cjMTERNTU1IamykwoMwyAmJgYxMTHIzs72mh5MSEhgMwOhSA9KscTgHsNbijCsKjQ/Tjmj5+fe5XJ5QPxxmUyGqakptLe3o6SkBCkpKX73Od5Gz4e39OD09DQmJiYwODg4Jz0oBaQw+vm8BaGqQqfyIHDKGP18MlZiDZEQAr1eD7vdjpqaGsFu8WIvuJHJZEhMTERiYiIKCgrmpAfNZjO6urqQlJSEhISEgNKDUhm9mNiLL1Uhu90Og8GAzMxMj7qBk30QOCWMfr7cu9gf2Gq1or6+HgzDoLCwUNQ6eLHN9P7ATw8eOnQI8fHx0Ol06OrqCig9KJV7H2hcgDsIOJ1OjI+PIzk5eY6gCPUETkZBkZPe6AOpe/cFrVaL9vZ2LF++HAaDIaCCG7qP0WgEwzCIjo4Wtd/xBMMwgtKDSUlJPqsHgzFYimCDgRQul8uDKAR4qgoBxyoITyZVoZPW6MXIWPmD2+1GZ2cnpqenUVNTg4iICBiNxoCM3u12o6uri9Wy79BZ0GZQ4vKqTCzNSoVSqQz4OhcavtKD/f39MBqNiI6Ohlqt9kgPinHNfUEqo/d2nFNBVeikNHqxMlbzgbrzycnJWL16tUfUWGwA0O12o6enB5HxanQiE7u/HEXzqB2AHV8M9eH/rR5GlAIsrZaumxfLTO8PND24ZMkSEELY9GBXVxfMZjPi4uJgtVpht9uDygxIMXDQ4/jLtgAnn6rQSWf0oXDnV6xYAbVa7fEZbT0lFJOTk3i/oR9fTUTg86ER2JxulKTF4J4LlyI1VoWfvNaKZzuj8cQVK2A2TnusmxmGQUJCwqJg1AmFr/RgY2MjOjo64HA4EB8fD7VaLTo9KMUSAZhx78Ucx9cgYDAY0Nvbi5KSEnYQGBsbQ2Zm5oJwH8TipDH6QKi0vuB2u9HR0QGDwYA1a9Z4fSCFzr5j01Y893+teLN5HGNmghilDVuqMvHNqgyUpEax12lzunHv6224f38XfrNlOZsCtNls6OzshF6vh0ajQXR0NJKSkqBWqxe94AYXND0YERGBiooKMAwzb3pwvupBKd37YAqUuNkfOhBRL/PHP/4xHn30USxfvjzo65QaJ4XR09z7wYMHcfrppwdlCBaLBfX19UhNTfVw5/mYr8rO6XLj/zrG8cqXg/i4cxxuAtTkJuC6/CicXRiH4oI8OJ1OD0/hsvJ0jJvs+P37PUiJVeInFywFwzCIiIhgZ8SMjIw5jLr4+Hh2EFjstfbAsei9v/QgrR70lh4M5Zo+ENCAINezNJvNiI2NDfrYocAJbfT83Hsw7i/Vrevo6MDKlSuRlJQ07/be1vS942bsPDKE12pHoDXakaACtlYk4/qzl6EgJQb9/f3zHvO607KhMdjx/KEhpMVF4IZ1OXPOyXWZ3W43jEYj9Ho9mpqa4HQ6F32tva/fyFf1IF3myOVydnAT65b7gtvtluQ7okbPhclkChu91PAlYxUopZYKVfhy573tQwiBxe7C281j2HlkGIf7JiGXMTgtJwZXFwP/87VKJCUmzNmHew98zsBd5xdCZ7Tjjx/0ICVGhY0V6T6XEjKZjGXU5efnz5ktZTIZm0ILNhAote6+P/CrB+12O5senJiYYJl186UH/cHlcvkVNRF6HP4gRLUPFyNOSKP3JWNFKbVi1mlmsxkmkwmpqalYsWKFoIeHEIJ2rRW7anX47ysDMNpcyFNH4Y5zC1EaY0S8kqC0tGbOLCIk4i9jGPzyshLozQ78Yl871DFK5At8Lvmzpd1ux+TkJEZHR2E2m1FXV8em0GJiYk6YeAAAqFQqpKenIz09HWNjYzAYDFAoFB7pQbGxjlDGBqTyIkKBxXlVPuAv905164Qa/djYGDo7OxETE4OcnBy/D8qk2YE36kew8+gwWkeNUMkZXFyWjiuqs1CaqkJDQwOWpC3xeSyhwT+VQoY/XbES33m+DnfuasavLszAijT/JJ45x1GpkJaWhrS0NBgMBixbtgx6vR69vb2s+0kHAX9R5sWUOXC73YiIiMCSJUvmTQ/6qx6UyuidTuccgs9i+a684YQxeiG5d5pG80dwcbvdaGtrg8Viwdq1a9HQ0OBzBna7Cb7oncDOI0N4p0ULu9ON0iVxuOtr2ahSu7CmshQajQa1tc0oKytDQkKC1+MA4ph1sREKPHVVOa59rha/eH8Uf9iQhSVLBO3qE1FRUcjKykJWVhYIIWw8gFthR1No/O9wMT3IfGP1lh7k3xsNeCYlJbHLN29r8UCvh3+cxfR98XFCGL3Q3LsQhVqz2Yz6+npkZGRg+fLlrMfA3290yordtcPYdXQEgxMWxEcqsLV6Ca6ozsKKzDh2fdnW1gaj0SgoFiCWZJMSq8LTV5Xh2n8dwX3vjODFnCykxEoToWcYBnFxcYiLi2Mr7KampqDX69mAI5cktJjgL0/v7d5oenBoaIhND5rNZkliFUImmsWERW30YnPv/spkR0dH0d3djdLSUo8Hmc2vOt34qF2HV48M4ZPZVNvpBUnYfu5SXLAiFZHKY6O5w+GARqNBXl4eqqurBY3qgTDr8pOj8dAFS3D3W0O49aVG/OtbFYiJkP5n4wb9gGPRc41Gg46ODiiVSthsNhgMhuPelUfsetlbenBqagrj4+Nob2+HXC73mR4UAr7HECzjMNRYtEYfCJXWF0vO5XKhra0NNpsNa9asmTMqDxlceO2DXhxoHYfe5EB6fAS+f1Y+rli1BDnquWtpKosVFxeHwsJCwffENXoxRrM8NRL3nZOOBz8YxfZdzfjrtjIo5aGle/Kj50ajEQ0NDWzgjMpu0cDZQiLYtbhcLodarUZ0dDSWLVsGuVzuMz0opHqQb/Qmk0lQIdXxwqI0er7yiVAD8eamm0wm1NfXY8mSJR7ReZPNiQNNM6m2owNTkMsYnFuSgiurs3BmUTLkMu8VYt3d3RgfH0dpaSmGhoZE3RfX6A0GAywWC5KSkgQ9wGuyo/Hgpcvwszfa8bM32vCrTcshW8DZNiIiApGRkSgtLQUhhO3K097eDpvNFjClNhBITc7xlR4cGRlBW1ub3+pBvtEbjcZFm6MHFpnRE0JgsVgwNDQkKJrOB9+9HxkZQU9Pj4c7Xz80hVe+HMK+xjGY7S4UpETj+qp4XFmTh6KcdJ/H5jatoFpzgVbZ9fX1YXh4GLGxsejq6mIfKrVa7TWVRgeLTRUZ0BrtePzDXqTGRuCu84V7GcGCL1HF7crDV9xxu90eJCF+ZFvKawkGvgYPbnoQ8F49SAeB6OjoOXn6xUzMARaR0dPcu91uh0ajQW5uruhjUPfe5XKhtbUVDocDa9eunSmAmLbhN++0482GMUQpZWyqrTonAe3t7UiM8r2Om5ycRFNTE9u0gp4rkCo7jUaDpKQkrF69mn3o6ENFU2lxcXFQq9VeqbXfXZcDrcGO574YRGqsCtedni36ewoU8wVQuWtmp9OJyclJ6PV6dHd3syQatVotidFLTZ/1B1/Vg93d3TCbzXA6nYiNjYVCoUBkZCTMZrMgYs5bb72F22+/He3t7Z0A/k4IeYz7OcMwXwOwF0DP7Fu7CSEPzX72DQCPA5B723c+HHej51NplUqlqOo1LmQyGSvplJWVhZycHDjdBP/6rA9//rAbTjfBbecU4Ib1eYiNVHjs582ACSHo7+/HyMgIVq1a5bFOExuUM5lMaG9vR2RkJMrKyuByudhz8h8qg8EAvV7PtrFSqVSIiIhgH9K7L1gKncmO373fjZRYFS4pSwvo+xIDMfeqUCiQkpLiUTREvQBa20D5AdHR0aJnbamMHhCvnOQtPfjll1+ycaP6+nrs3bsXSqUSGo2GnST4cLlcuO222/Duu+9i6dKlKwEcZhjmdUJIM2/T/xJCLuVdgxzAXwBcAGBwnn294rgavS8Zq0CEKoGZtdTk5CRWrVqF+Ph4HOqdwEP7WtGhMeGc4mT8bEMJcr0E5rwVzzidTjQ2NkKlUmHt2rVexRaEGgKNgOfn58NoNM67La1G41Jru7u7MT09jSNHjkChUECtVuOnX8/ChMmBn73RhqRoJdYXzl8rECyCcamp2EZaWhpMJhOKiopYL4ASaeggIIQWK6XRBwuaVcrLy4NcLkdJSQl0Oh3ef/99XHXVVdi4cSO2b98+Z79Dhw6hqKgIhYWFIITYGYZ5CcAmAEIMdy2ATkJI9+w1iNn3+Bm9r9x7IH3fXS4XWlpaYDKZkJeXBysTgYd2NeKN+lFkJUbir1dX4tySFJ8PLZ8eazAY0NDQgPz8fCzxwYgR4t4TQtDZ2YmpqSmsWbMGJpMJBoNB1L3J5XLExcUhIiICubm5sNls0Ov1GBsexLcKLRiblGH7q43429WlqMxV+z9ggJBqLS6TyRAdHY3o6GiPOnu9Xo/m5mY4nU4PkpC31NxiI75wB6HIyEhkZWXhoosuwr333utzHxq34mAQwGleNl3HMEwdgGEAdxFCmgBkARgQsK9XLLjRS1n3DhxLJWVnZyM2Lh57mvR4vrYHdqcbt55TgO+fmY8o1fzrNq4BDw0Noa+vDxUVFfMGY/zN9Ha7HfX19UhISGBLdE0mU0D3yD1XREQEMjMzkZmZiZWEIK9wAt97uQW3vtyEe9eqsGyJmg0ySd3MMlhD82as3jwbShKiJbY0HhAfH8/+Votlpgfm3peQslofzw7/zSMA8gghRoZhNgB4DUAxAG8/hOBReUGNXkoZK+CYgZaXl6Nn2o2f7K5Ht96GM4uS8fMNJchPFpYrlclkrDvvcrnY4J+/fXwZ/dTUFBobGz0Cf4D0ApcMw6AwU42/X7sK395Ri780MXhiaQImJyfZKjsaEJSiyi4URs8HzaFTpSKHw4GJiQmMjo6ivb0dERERsFqtsFgsi1ZExGg0spF/X8jOzsbAwIDHW5iZzVkQQqY5/9/PMMxfGYZJwczMnjPfvvNhwYyeBuukkLFyOp1obp5ZvlRV1+DpT/rx90/7oI6S475zUvGtr1eIOr7D4cDAwAAKCwsFpwp9VcwNDg5iYGBgTuCP7hMK3fvClGg8ubUMN/6nHve9M4h/XluJoiI57HY79Ho9BgcHYTab0djY6CFUudAIZIZWKpVs0RBwTORkaGgInZ2diI2NZT0BMdJUodQcFJKyW7NmDTo6OtDT04PCwkIVgKsA/A93G4ZhMgCMEUIIwzBrAcgAjAOYBFDMMEwBgCFv+86HkBt9oO68r1mBrrdzc3OhZ+Jw5d+/QqfWhMtXLcGNNWo4zNOiDH50dBQDAwNIT08XlSbkGyKNK7jdbqxdu9arax0oI08IqrLj8btvrsDtrzbhzl3NeGJrKVQqFatWazQakZ+fD71ezxJq6No5KSnJr2ezUDO9P0RFRUGlUmHFihVQKBQwGo1z2nL5KhriQqolgrd7EuLeKxQKPPnkk7jooosAoAXAPwkhTQzD3Dx73KcBXAHgFoZhnAAsAK4iMw+Qk2GYHwB4GzMpu3/OrvUFIaRG76vu3R+ocXC3J4RgaGgIAwMDWFFahue/0uLp/7YhNVaFv11bhbOLU2Z05IzCIv9utxvt7e0wm80oKipi1U2FgnttFosFdXV1WLLEd1kt977EQuh+5xQn4/4NxfjFvg78Yl87HrmsxCuhJjc316PApq+vDwzDsG61N+rpYjF64JjBcgtruPdEiTQAFkRyi39PQgU0NmzYgA0bNgDAUvrerLHT/z8J4Elv+xJC9gPYH8g1h8To+bl3se68XC73YDk5nU40NTVBJpMhbWkZbny5BQ1D09hSlYn7Li5B3GzOXShhxmq1oq6uDmlpaSgpKYFGo4HVag3gTgGdToe2tjaUlpb67fcm9ZreG75ZlQmt0Y4n/68PqbERuOPcAq/b8QtsuMo0BoMBUVFRrNssFY881Ew6/j05nU5MTEx4KAvTgU2pVIZMQOOUY+TR3PvRo0dRWVkZ0I9MjV6pVHq48wfHGDzyzGFEKOT487ZyXLTSM1giRJaaGilX1joQdh0hMx1re3p62AYY/rAQRg8A3z8jF1qDHf88OIDUWBWuXZvldx8u9ZTLOuvs7ITVakVUVBTbAy7QMlIpU21CjqNQKOZ05NHr9RgYGMD09DRcLheGhoaCKho60fTxAImNnpt7N5vNAR9HLpfD6XRiYGAAg4ODKChZid98MID9jWM4vSAJv/lmKdLj5wZt5jNeQmZ6yE9MTMwxUrFG73A40NDQAEIIVq9eLXjGCLV7z93+pxcVYdzkwG/e7UJKrBLJIs/HF+AcGRnB0NAQ6uvrAWBOGk0IpGpSESi46U6DwYCenh5W7txqtXqQhIQWDfky+ri4uFDcgiSQzOhpOg6YMaJA9OooGIZBa2srIiMjEZuzAt96vgnDU1bccd5S3HhmvtcKOMB3PT03Z15TUzPnwRNj9NTzKCwshNVqFd0sgRrv6Ogo2zwxFMq1chmDxzYvx/dfqMe9r7dhe5USawI8lkwmQ1xcHOLj47F8+XKvaTTqNs9Hq/UnfrGQIIRApVIhJyeHLRqiJKGhoaF5i4a4OKVnev66nbroYo1+enqmu0tOTg6OTEXjl7u/QnKMCs9/ZzVW5ybOu683935iYgLNzc1YtmwZ6+Z520+I0Q8PD6O3t5cl7nR3dwu+L+BYmq+1tRVmsxkZGRlsTp2bn+aXbwbqIUQoZPjzlaW4/vk6PFFrxtoqI5ZnBPYwcl1zfhqNW4BisVh8zpiLiUnHfzZlMhkSEhKQkJDAFg3RQGd3d7fPGntvz7jFYjl16um5xkNddKFuEiEEAwMDGBoaQmpqKp4+asS+5n6cWZSM336zFOoY/8fhzvSEEPT19WFsbAzV1dXzrtn8GT01VLvdLoi44wsOhwNTU1NISkpCZWUlHA6HR1EKlaoyGo0elXbBICFKiaeuKsfWZw7hlpca8e/rq5CVGFirJV8GS2m1WVlZXmdMbnBtsRi9v+i9QqGYoyzMDXRGRkZ6xIS4IIQImuwEVNldA+Ansy+NAG4hhNTNftYLwADABcBJCKkRcNsz9yZ0Q7GgM70QOBwONDU1QalUYu3ateju7kZ+oh0//Fohbj2nADIf7jwfdKZ3OBxobGxEREQE1qxZ49elnG8mpZH+9PR0wRLZ3jA1NYWGhgZERkZi6dKlcwYZ7nqTX2lnt9uhUCgwMTGBhIQE0S5yRnwE7lytwm+PuHDziw3YcV0VkqLFBeOEehreZkwaQR8fHwchhDWY4ynDLTZlxw900nLokZERWK1W2Gw2JCUlCfbKBFbZ9QA4hxAywTDMxQCegSfH/uuEEJ3gm5hFyIxeoVDA6XT63W5qagpNTU0oKChAZmYmgJkB48oKNftaKBiGgdPpxOHDhz2O5w++Znoqi+WtgaUYDA8Ps3Th9vZ2v9vz+eharRYjIyNstR41GjHptKxYGZ7cuhw3vtCAH7zciL9dU4FoPzUJXATqmnMj6OPj49BqtVAoFHO0A4RW2EmFYJRwGYZhvRv6zCUmJkKv1+Omm27C4OAgtm/fjosuuojm4edASJUdIeQzzi6fY4ZuGzQkNXpva3pfoO736OgoqqqqPB5eMV4CF0NDQ7BYLFi/fr2oQArf6Akh6O3thUajwerVqwPuPEqltqk2H32PQqgRyeVyREdHo6ioCMDMGlqv17PpNKHMulU5Cfj15uW4c1cz7trdgsevXClYa08qco5KpfLQq+dX2FEFXl/BM6lSnlIKcSgUCnaQPnDgAM466yxs3LgRbW1tPvcTUWVH8V0ABzivCYB3GIYhAP6XEPKM0GsO6Uzvy3C57re3WnWxRu9yudDc3AxCCGJiYkRHTrlG73Q6WTdcyNLAlzHY7XbU1dUhOTmZldoW4vl4A//43NJUPrOOW2QTFxc3Z9/zSlJw3zeK8PCBTjy0vwMPXbpswVxs/nflrcKOr7jDD24uVG96oaAiJxR2ux3R0dE477zzcN555/ncT2CVHQCAYZivY8boz+S8fQYhZJhhmDQA7zIM00oI+VjINYd0Te/tIafSU0uXLkVGRobXfWUyGcvm8wcqfJmdnY3s7GwcPHhQ9LVSozcajaivr5+3jp6/nzejp1V2/IxBMOQcX/t5Y9bRIhsqV61Wqz08jK3VS6Az2vHUf/uRGqfCj77mnbXHP3+oabj8tlz84GZsbKxkGvxStZ3ip6Wphp4/CKmyAwCGYSoA/B3AxYSQcfo+IWR49l8NwzB7MCOssfBGP597T13msbExrxVoHhelUAgi91Ad+7KyMsTHxwd83TKZDDabDfX19SgvLxdMrKAzD3fGGBoaQn9/v6RVdmLALbIhnE4vVqsVhw8fZtNON52RA43Bjr99OoDU2AhcXRNk+xwBEDtw8IObRqMRWq0WZrMZhw8f9iu2MR+k6m4TqBKuwCq7XAC7AXyLENLOeT8GgIwQYpj9/4UAHhJ6zSF17ymf3W63o7GxEVFRUV7deT780WnpWtlqtXrVsRcDWnhjt9uxbt06UcfiGjF//e7tIeQbvVAjCIbJRwtSNBoNqqur2Uh6Z2cnNi5RYUgXjV+93Ql1tAIXrfSttSfFTB+Ma07vRalUwmg0orS01ENsw9+yxtu1hEJcUygxR2CV3f0AkgH8dfZ+aGouHcCe2fcUAF4ghLwl9JpD7t5TcgxfUMLfvr6MntZTp6WlsWtlPoQ+oHR2p9xrsYMHXRbYbDbU1dUhJSXF5zUBc9fmC52uksvlHoKVFosF96h1uHtfP+55rRXG8VGcVZKBpKSkkPSyk4KRRwcOvtgGf1njrxlHqIxeqBIuIKjK7nsAvsffj8xo41UGes0hc+9lMhn0ej30er1fcgwfvoxeq9Wivb0dK1euZNewfAjtXEsHo5KSEqSkpGBsbEzw9VEwDIOpqSl0dHSwxxGzbyi2FYOoqCgszcvBP7+TgW8/V4vffT6NpJgIJMyuNalRBbN04iKUAwd/WcNvxsHPcITSvV+sfekpQjLT2+12dHV1wel0Yt26daJHVG/xAK7A5HwsP3+cf8KRtRY7GPFhs9nQ0dHhN0YhBUIZC0iIUuLpq8tx7bO1eOS/E/j39VVIiZZDr9djeHgYra2tkMvlUKlUsFqtAacwpZrphZCt+M04+NoBLpcLCQkJiI2NDeqaTrTuNkAIjJ7OoDk5OZiYmAjoC+XSaakLnpiYyApMzof52lW7XC62Ln/NmjUBj/RcWm51dXVABk8VhYQsKRZiGZCZEImnri7H9TtqcdOLDdjx7SoPBlp/fz8mJyfR0tIiKJ/uDVLFBcQ+U94yHPX19dDpdBgYGPAgO4nV3eN3txGimnO8IanRm0wmdHR0YPXq1SCEQKcTzRAEcMxw+S640H29setoai8nJwfZ2YETm+j6PTU1FYmJiQENapQLYDKZ2ABUcnIy4uPj561QCzWWpcXg8StLcfOLDfjBK0342zXliFLKwTAMIiMjkZiYiLy8vDn5dKVSyRrNfNTa4x0MpFCpVFAqlVi2bBlUKhUsFosH2Yn25fMW2+CD772ccu59bGws1qxZA4Zh2CYWgUAmk7EdYQKJB/CNntJXy8rKgsrzUo4BHYRoTb0YuN1uHD58GDk5OVi5ciUb7KRudExMDJKTk6FWqxeUlkqxJi8Rj21ejh/vasHde1rxxytWQiHzlC/j59OtVisbRTeZTKzRUIUa7r0H695LVZ7Lldzik52mp6dZsQ0ArJfgq+6BOwiZTCZkZfkXLTmekNy95z4YgTDQqECF2+0WxIjjg5vuExML8AdvKrc0Ty8Uer0eZrMZa9euRUJCAux2u0eZKg1AjY+Po7m5GS6XC0lJSYiMjAy4608guGB5Ku69yIFH3u7ELw904BcbiuedpSMjIz2otdRoBgcHARwT3JBilpa6Yy0f3L58wDEJbjpx8LUD+IO+mOj98ULIUnaBdKqhxTeBCFRwz0sVfPjNJvzB24NN1++0GSZ3/Somf06Dh9HR0UhMTPS6HzcAlZeXx3oBo6OjmJiYQENDA7sUCDSYJhRX1SyBxmjD3z4dQFqcCluKhXkdDMN4VNlRoxkZGYFOp0NUVBTsdnvAElULTcP1JsHNbclls9kwNjbGagcIVc2hZbUulwtdXV33eCmrZTDToHIDADOA6wkhR2Y/C7h5JRDClJ0YEEIwODiIwcFBVFZWIiYmRrRABYVcLofBYBDNDfCW6qPr97S0NOTl5QWkuON2u9HS0gKXy4Wamhp88cUXgu+FVqhFRkZCLpcjLy8Per2eHYQSExNZ5Z1QKNL88Jx8lq6rdKbi4mXiJaC4RtPZ2YmoqCiWyGS329mAoNCOPMezeSUwk+bMyspCVlYWXC4XDh8+DLPZjKGhIezcuZNth/a1r33N58DMLavNzs5GRETE1V7Kai/GTDebYswU4jwF4DQmyOaVQIjcezEzPG1cIZPJfOrFi4HRaMTY2BhqampERdX5ngldvy9fvpxdu/Lh7179DRpCQfejunU5OTlsoJOy6yIiIth1tlSNLBiGwf0blmHc5MCTn2uRGCnHlUEWd0ZFRUGtVrP3QFNpPT09Xgts+FhMklvAzKBWUFDA/t1yyy344osvsGPHDrz88ssoLi6esw+3rHYW3hpQbgKwg8w8YJ8zDJPIMEwmgHwE0bwSOM5da2kfumAj6sCxZhNWqxVFRUWi02jc9Tn1OvwFEecz+unpaTQ0NIgm7QgFn11Hy20pGYXrBQQzkCpkDH67ZQWu+9eXeOzjMSzNSUd1TmDBUL5rzmfV0QKbvr4+j1p7tVrNxmMWUx87frqORvsffvhhtgzaGwSW1XprUpnl433BzSuBEBu9t4IUipGREfT09Mxb4CI0xcNtNhEoSYb2s+vo6IDL5RKUx/fl3tN7q6qqkiSoI8R74kag6Qw6Pj7OptTsdjtMJlNA/eCjVXI8cF4G7jowjB+80oQd365EUar4+/I3S3tTDxofH0djYyMru+VyuRZNHtwbq08IOUdgWa2vJpVBNa8EQuze85tWAMcCYzabbV69Obqvv+opfrOJ/v7+gCLdhBDU1dUhMzNTsCvurYCmo6MDRqMxKC29YMGfQa1WK7766itWuFJMOyuKhEg5HrlwCX58YAQ3v9iAf1+/Chnx4lKKYvL03Fp7ruxWX18ftFotdDqdaPUgqRFoowuBZbW+mlSqfLwvGCF9KrlNK4AZF7S+vh4ZGRl+9eZ8yVlTEELQ3d0NvV7voWMvpOEFHxMTE5icnMTy5ctFLTO4cQCn04m6ujrExcVh1apVi0YAEphJqalUKpSXl8+hpNIBIjk52S+xJjNehaevLsP1z9fhlpca8Oy3KpEQJbxIKRhyDg1qmkwmREVFIS4uDnq9Hh0dHQH15ZMC3mZ62hhkPnDLamdz+t4aUL4O4Aeza/bTAEwRQkYYhtEiiOaVQIiNnquTR/Oc8xXLcDFfpR3N5cfExMxpNiFGgAMAq8BLGXFiQJcvJpMJdXV1onT5xJ5HqgeVT0m12WwYHx/3INYkJyfPYaNRgy1Jj8XjV5Ti5pca8KNXm/C/V5cjUrmwNFwqH8ZXDxofH/dQD0pOTvYaEAxVsQ2Fv5gDt6x29hl/xUtZ7X7MpOs6MZOy+87sZ0E1rwQWYKZ3Op1ob2/H9PS0KIKMrxmb22zCm/KOPw+Bwu12o7m5mSUBtba2il4WMAyD6elpNjYhdNBYCEqtUERERLDEGipfPT4+jv7+fg/j4V7z2vxEPLpxOe7e04J79rbi999c6bMBCRehIud4G8j4ijt0KRARERGyslpCiODfllNWCwCPzO7PLaslAG7zti8JonklEEJGHkVzczPS0tIEE2QovM30/GYT3iDEvaey1hkZGcjNzQXDMKLJRIQQTExMsOt3oZRZOmuLLa1diIGCK19dWFjI1qj39/djYmKCLUZRq9X4xspUjJvseOydLjz6did+9o0iv/ckRbpNyHfnTXFnfHwcTU1NcLlciIuLg9PpDNr4fc30i2lp5w0hm+n1ej00Gg1yc3PnTV/4AtfoxTSb8EeYoUU8fFlrMZRaWq1nt9uRm5sriiNPDXh6ehqdnZ1ISEjwu54+XuDWqPf09IBhGJjNZpZeuz5FjWuqU/GfIyNIi1PhpjPz5j3e8aiy46oH5efnw+l0sgzHL7/8kqXVUn6D2Ao7/kx/IiAkXWu7u7uh0+mQlZUVcMqKGj23rbSQZhO+3Hsy20FneHjYa/5dbJvrzMxMqNVq0UFDhmEwMjKC/v5+FBcXw2w2o7e3F2az2WM9fbwi//MhOjoaaWlpLL1Wr9fjslwzukdkePL/+qB0WnDt+kKfS7jjVVrLhUKhQEJCAgwGA1asWAGLxYLx8XGPCjuhvwHf6K1W66JuZ0Uh+ZPV398Pu92OmpoaDA4OBlxpJ5fLMTU1hfb2dlHNJry593T9TgjxmX8XYvSUpUevZ3h4WFTQkMx2RhkdHcWaNWvgdruRkJCAzMxMj2BUb28vFAoFUlJSkJycvGDuvT9wDVapVLL19k+WuHDri/V4/DMNiGUalakza2waHKVGuhiMnn+MqKgoVknZm5w4vQ9vunsul8vDyxOqhHu8IbnR03prYMZwxRgFBSEEk5OTMJvNWLNmjajiEr7xelu/C9mPD28qt2KMkdbQA0B5eTkUCgXb5ZeenxuMslqt7AxksVjgcrkwPj4eNMMuUMxnsCqFHH/eVoEb/l2Hp+vNePqqUsRHz7jRbW1tiI6ORnJyMpxOZ8hltIWAzx2h8CcnHhMTwy4FIiIi5hznRFDNAUIcyJPL5awirlA4nU40NjbC5XIhNzdXdDUZ1733tX73Bl9GT9VyqfIu1+UTavRWqxW1tbXIzs72ePDne3gjIyPZwg6LxYKmpia2ukulUknOsw8W0So5/rKtDN9+rha372rBjm9XYfnyVBBCYDabMT4+DrPZjNraWnb2DKQvn1QzvZCB05vuHi17pqlohULBFguJLavV6/XYtm0b3nvvvQ4AvQC2EkImuNswDJMDYAeADABuAM8QQh6f/ewBADcC0M5ufu9sZH9ehNzoxbj33GYThJCAvARKp+3v78fw8LDgtlTeovcOhwN1dXVISkpCSUlJQFV2dElA+QljY2MghMDtdsPpdLIzl0wm8/kwy2QyKJVKtniDrkOpdHeoq+0AYTNscoxqRmvvudpZ1l4V0uIi2EIhjUaDiooKTE9Pe/Tlo6IhQgYwqd17oeCXPbtcLjQ0NMBgMOCrr77CxMQE3nnnHbhcLsHeyGOPPYbzzjsP7777bjHDMPcAuAfHutRSOAH8mBByhGGYOABfMQzzLqeq7o+EkN+JuZcFI+f4A21cQbn4o6Ojor0EACxnm7alEuoK86P3dABaunQp0tPTfe4z30xPG1fyhTOdTidkMhkUCgXcbjcIIXC5XOx3JZfL2YGAnocL7jqUSlfRajtqRFLX3At9kHOSovDUtjJ859/1uPmlRjz7rUrERyrYYyiVSrZQiMY46ABGO7/Op70nhdH7cu/FQC6XQ6lUIj8/HzExMRgZGcEHH3yAI0eOoLKyEnfffTeuvfbaeY+xd+9efPTRR/TlcwA+As/oCSEjAEZm/29gGKYFM0U3gqvq+FgQGu58oO4zXb9TFphQkg0X1I2Wy+UoKysTtfbjMvk0Gg06Ozv9drvxNdNTxR6DweCxJKCGMzIygoyMDERFRbEPNtf4uX316Dl8DS5c6SquEdGae+pKL2QgcGVmHP54+Urc9nIjbn+1CU9fXY4IxdwBjCtVRUttudp7KpWKXUPTQVNKVl+w4EbvMzMzccYZZyAlJQX3338/jEaj3/3HxsZYBucsxXZe8QeGYfIBrALAFWX4AcMw3wbwJWY8gglv+3IRUvd+viaWwLF68+Tk5Dl8dbEcer1ej5aWFqxYsQKtra2iHwx6vq6uLpbP74896G2mpwG76Ohoj3uiBl1cXAyNRsOuC5OTk5GSksKub7mDAB0ATCYTgJnlBtV186XVxjUip9OJyclJaDQamM1mNDQ0sAOEWP09sca2vjAJv7ysBPfsbcVP97bit1tW+N2Hr73HF6ykEmPBDmBS9bHz1d1GoVCwclvnn38+RkdH5+z7yCOPiDoXwzCxAHYB2E4ImZ59+ykAD2Omyu5hAL8HcIO/Yy0IDdcb/CndCo0HUHnm0dHRoNpKAzMlscnJyXP4/L7AN3qLxYLa2lrk5uZ6iCNSgyeEICoqCvn5+SxRZHx8HENDQ2hpaUFsbCybplOpVJDJZNDpdOju7sbKlSs9PAvqovoaAACwab+UlBRMTU2hoKCAZaa53W5BKrzcexA7kF5SlgadyY7fvdeNX7/bhfOTxBkrV6XG7XazA9jRo0ehUCjYAUJsuXAoW1rxA3nvvfeez/3T09MxMjKCzMxMzApkaLxtxzCMEjMG/x9CyG76PiFkjLPN3wC8KeS6F9y9p0bqr9mEEKOnLaoZhglIRJPCYrGgs7MT0dHRWLHC/4xEwTVCfsCOgs7YdHsuFAqFh7a8wWCATqdDbW0tGIaBQqGAzWbDqlWrPGZmrhfAPT59AH15AXz9PW4zi9jYWDagFoyAKB/XnZYNrcGO574YhL1YgbVrAzsOrQOgsRqqwEvLhYMh1QQK/kBoMpl8qix5w8aNG/Hcc8/hnnvuAYDrAOzlbzOrlfcPAC2EkD/wPsucXfMDwBYAjULOu6BGL6bZhD+j5wpn5ObmBnyNdFmQm5srOnBIZ3oasOMOYjRCT4tM/M1E/PrxlpYWGI1GREVF4ciRI0hISEBKSgrUajUUCgVr2PQcXOMX6gVwVXgpP72xsRGEENYLoKSUYNbSd55XAK3Rhl1NWlTVjWJzpfcW5WLAVeClstXcKjs6gHmrsguVzp63mX4+3HPPPdi6dSt++tOfdgDoB3Dl7DGXYEbwcgOAMwB8C0ADwzC1s7vS1NxvGIapwox73wvgJiHnDemanv+FiGk2MZ/RU0MVWqbrC7SsdvXq1TCbzYLaY/MxNTUFh8MxJ2AnxuC5oDyFuLg4lnZMmWLU1VcqlUhNTUVKSgqio6M9BgB6bu4A4C/LwOenU4otJaXExcXBZrMF3seAYfDLy0rQOzKOB/a1IzlGhbOKhDEsBR2fJ1vNl93i6/BLZfT871SoEi5FcnIy3n//fWBG/JJ73GHMlNWCEPIJvKvlgBDyLXFXPIOQzPT8hyyQZhO+lgZ9fX0YGxsLav1OC3icTifrcVitVlHZAqfTiba2NhBC5gTsAjV4q9XKDozcunwuU6y4uBgWi4VVDKJprpSUFJYkQtNJ1PjHxsagUqnY7IQ/L4BLsaXLjra2NnR0dLDdbHzVqvuCUi7DbVUqPNmswI93N+Mf11SgPEuaxph88KvsqBdAC4VcLhfi4+MlyQRwIXamP14IqXtPCIHVakV/f7/oZhP8dBh/aeBvpPb1g9rtdtTW1iI1NRX5+fnsNkILboBjAbuMjAxMTk7OidDTc4t5oKanp1leP52xfCEqKgo5OTkeqri0oy+lvKakpCAiIgJDQ0OYmJhAeXk5O5Byg4H+iEF02REbG0vlmj1q1ePj49kBx99aOkrB4K/byvCt52px26zWXn5yaLnqDHNMhx841sdOo9Ggr68v4FiGt+dLaG/6442QGT39cgGguro6oE411Fug6/esrCy+iqjPfb3lYqlC7bJly5Camup1H3+gWYeVK1ciKioKExMzaVEy25CSHksMNBoNuru7UVlZKbpgg6uKS6miOp0OjY2NMBqNiIyMRElJCZRKJRiGmZMSpAMV/T+fGERBH3IuNZWKbuh0OlZ6y19EPSV2hrX3bZa1twopsdIFDv1BpVIhIiICRUVFiIyM9IhliMloeFsimM1mUe798UJIjH56ehr19fUoLi5Gd3d3UHlVSjQRs3731q6aMv4qKyu9jsZCjJ4W3dCAnc1m84iki53daSZDp9Nh9erVgjrYzgcaoY+KisLk5CSysrIQFxfHRujj4uLYlKBSqZwTC/BGDPJHD+bOolR6iyvASSPq3N8iTx2Fv2wrww3/ntHa+9e3KhEb4f9RlIpgxO1jx49lCO0r6C0DcEq796Ojoyz1tL+/P6AUCSEENpsNnZ2dotfvXAOm7Dgq1+XLsOYzekKIB2uQ68Y6nU6PaLlQ0C4vbrcbq1atkiyabLfbWa9oyZIlAMAWjExPT0On06G/vx8Mw7AeAu3R7o0YRL0AborQ17XypbdoqXBPTw8rw202mxEdHY2yJXH4w+Ur8cNXmrB9ZzP+uq0MKsX830Go8usUQvoKqtVqJCQk+JS/PmVn+pKSEtaAKEFHzHqJrt+pfl0gSwPKZa+vr0dMTAyqq6vnNUpfcln0GLGxsaiqqvJYv8tkMkRFReHQoUNsSi05OdnvAEePmZSU5BFXCBYmkwkNDQ0oLi6eky/mrm2XLl0Ku90OnU6Hnp4emEwmJCQkIDU1FWq1GnK53MML0Gg0sFqtUCqV7ADgLxjIL1O1WCw4cuQIy65LTEzEiuRkPLChCD97swM/e6MNj21eDtk834VUgTchg4c3XgO3kaVSqYTT6YTVamUnJJvNJorpKKTKbvZaegEYALgAOAkhNbPvqwG8jJmuNz735yNk0XsKsZV2NEiWk5MDo9EY0I8sk8lYhdq8vDx2xvN3zfyZnl4L/xjUFQaA0tJSEELYlFpPTw9UKhU7i/LJRxaLha0k9FXIEwgmJyfR2tqK0tJSQbONSqXymJVp0U5XV5fH9dPmk6tXr2YHM1/EoPkMia6lKyoq2PONj48j0zaBq5ZH4qVmLRIiZbj3G8t8/uZSptrEPldUgjs1daZcWKPRYHBwkK1x2L9/P5RKJRwOh2CvVmCVHcXXCSE63nv3AHifEPKYgP2P3YugqwsCYirt6PqdNq4YGBgI6AdyOBxoaWlBZWWl4BQh372nATt6LRTeAnYMw3jMamazGTqdDi0tLXA4HGw0HQDLLxB6XUIwNjaGvr4+VFVVBZTGpGw3qjlgsVig1Wpx9OhR2O12ZGZmwmAwsKW7fGIQdwlAj8f3AriimPzzlZWZYWHa8dKRMbiN47iqKtVrqfDxbl7J3VepVCI+Ph7FxcVwuVwYGBjA7t27sW7dOpSXl2PHjh1+jyOkys4PNgH4mtj9Q270QmZ6Qgh6e3uh0Wg8Gld465Dj7zh9fX0wGAyiDYtr9PyAHYXQgF10dDRyc3ORm5vL8us7OzsxNTWFlJQUWK1WxMTESFL00dfXh/HxcVRXV0umqxcZGQmz2czyAiYnJzE2Nsaq4FAvICIiYk6JsDdiEF06+frOoqOj8dCWSljRildatCjKlmE10bFdbmlGQEqjDxbcQLFcLseWLVvwxz/+EUePHoVWq/Wz9wxEVNkRAO8wDEMA/C8h5JnZ99MpDVdIlR7Fgsz08xm9y+VCY2MjlErlnPU7v0POfHC73WhqagLDMMjIyBBtAJRQ1NbWNidgFwzhRi6Xw2KxQCaT4ayzzmJnUaqDx2XWiQENLjqdTlRVVUlmDG63G42NjYiJiUFhYSEYhvFwa2lKsKGhAW63m/ViqBYenxhEl0I2m439v7dYgIxh8OjGEkyY7fjNR8P4y7ZSrFtTzKruNDc3sxV2k5OTHtp7xwP8QB53QOKmgyWqsjuDEDI8a9TvMgzTSgj5OMBLP75rerPZjLq6Op/UXKHxAJvNxpJlcnNz0dnZKboW3+l0ssbJD9gFavBu90xvenpMmUwGlUqFhIQEFBUVwWq1QqvVssw6fpmtL9CBMjY2FsuW+V4Di4XD4UB9fT3S0tK88iG4wS0uXXdgYICl63pLCVqtVrS2tiInJ2deYpBKIcOfrijFd56vw/adzfjXtZUoXRKHmJgY5ObmYnJyEj09PRgdHfUgItGqxIUE3+gtFovX4jEpquxmabkghGgYhtkDYC2AjwGM0aKb+fbnY0Hce29ren7jSV/7+jPeqakpNDY2evSRF1uLTwcfriQVEBzDjhpQSkqKT0HOyMhID2adXq/HyMgIW/WWmprKGhAFTcktWbLEo3w3WFBtg7y8PMEBRj5dl5sSlMlkbDqwo6MDy5YtY38fb8QgYOZ3i1HJ8NRVZbj2uVrc+nIjnr+uCrnqGWOiegElJSUeKTUxxBqpcv385qpGo1F0jl5glV0MANmsak4MgAsBPDT78euz+z3ma39vWBCjt9ls7Gu6ftdqtR7rd1/7zme8vuSoxFBquQG75uZjCkTBMOxoo87CwkKkpQlaZkEul3u40QaDAVqt1q8BSQGa6lu2bJlgqXE++ClBm82G4eFh1NfXQ6VSQaPRsIbJTwnyiUEJEQz+unUlrv93A25+qQE7vl2FlFiVhwsttFSYP2hKlfbjy18HQsEVWGWXDmDP7DUrALxACHlr9hCPAXiFYZjvcvf3h5C79wqFglV+oRVkKpUKNTU1fo3J14w9H1kGEC61NTg4iMHBwTnkn0AZdsBM6qylpQWlpaWiG2JScMtsqQENDAygvr4eERER0Ol0YBhGEiHMqakpNDc3o6ysTFJiic1mw+joKNauXYvo6Og5KUEay6AuMZ8YlKeW4fHLl+Oml5tx60sN+Pv/lLFLAm/wVSpMqeBqtRopKSmIjIwMmYCGWKMXWGXXDaDS2/6EkHEA54k6KRYweu9v/T7fvlxQtzk+Pt5j7c2FP/eeBuwsFsucuv5gDH50dJTVxpdSlJIy29atWweVSsWSRNra2hATE8NG08Wua6mYZlVVlaRS2nq9Hu3t7R7H9ZYSbGlpgd1uZw2SnxJcnZ+M321ZgdtfbcKP97Ti3jNn9vcVDKTwVSo8MDCA6elpuFwuaDSaOZ15xYCfVTKZTCdEowtggaL3JpMJR48eFVVaC8w1eqpQ66tjLQVX5JIP2keeP2gQMtNxVKvVQq1Wi5oNCCHo6enB1NSUpKkzYKZjkFarRXV1NfuAcgtsjEYjdDod6urq2M9SU1P99sYbGRnB4OAgqqurJQ2CaTQa9PT0zFH74SIqKopNadJYBjclSGMZEREROHd5Gh68zIWfvd6KR21G/HFblYeEuL/6AMAz9mA2m9Ha2gqj0YiBgQEwDMMuA8T0E+TP9CdKowsgxO49IQQjIyOYmprCGWecIVqMkWv0tHS0oqLCrxvqy72n3kZ+fr5HvTpdUy5fvhyjo6Po7OxETEwM64LONxu43TMtsxQKhU/PIxAQQtDR0QG73e6Tm8+d0QoKClhqbXd3N0wm05w6e4re3l5MTEygurpa0m45w8PDGBoa8hig/IEfyzCZTNBqtR4pwdVJCnyzSIHdnU787QsN7rqgyGt9AOBfK4AQgoiICBQWFgKYCYzSVmJi+gnyC7pOlLJaIIQzPV2/y+VyJCQkiDZ44JjR9/T0QKfTCa7J9+beU7UdvrfBjdAnJCQgMTGRnUEpK00mk7EPJteFo+XDaWlpQUl28UFrD6Kjo1FaWip4IOFTaycmJqDT6dDR0YGoqCikpKRgenoabrcblZWVkua5abVgMAMJNzhHm2R2d3djYGAAF+cpYXJH4Z8HB5AcrcB3zsgXpBjE9wL4M7RKpWIFN/iyW/OVCkuxpj9eCInR2+12HD58GLm5uUhLS0NtbW3AxxocHERSUpJghVpgbvReSMCOe2zuDFpYWAibzTYnnx4XF4eenh4UFRXNqc0PBrSrTkZGhuDYhzdQnTiqeW80GtHU1ASHw4GIiAj09PQgNTXVa2NGMSCEsJ6FlCQhYCbmYDAYcMYZZ0ChUGDFyik0PdeAZz/rRXmEjl3mxMTEzEsM4moF0CWBN3iT3fJVKixECXexIiRGr1KpUFFRgZiYGHbtJRZWqxU9PT2Ijo7GypUrRe1LjZ4G7GgfOvojiSXcREREeHSU6evrQ2trK5RKJTQaDQghgqrr/IGKhSxdulTSgcTlcqGjo4MVEXU4HKzwhdFonFNhJxT0+3W73SgvL5dUempoaAijo6Ooqqpi3WwLE4FhgxM3npGH8vIsNhtAKcO09n2+lKDFYgHDMB7xAF+Yr1TYbDZjZGSEZVOaTCZRBVS0wq63txednZ3vwnsfuxLMVNFRFAK4nxDyJybAPnZACNf0dNTzVbI6H2juPDs726Ozq1DQEZ2qyFZWVkrCsANm+NI6nY6NpE9NTUGr1aK7uxsRERFsHEBs9J6mzoJJ9XkDlQfLzc1lg59KpdLDpaX30NXVJfgeaCwjIiLCa5+/YNDf34/x8XFUVVV5DEIvfzkMBgy2rl6CiIiIOZr43HvgVznK5XJWIKOsrAwARMUC+KXCX3zxBeRyOTo7O/H666/j8OHDcLvdHqW284FW2N1zzz1gGOZ9eKmQI4S0AagCAIZh5ACGAOzhbPJHIrKPHRDCNT1fHFMoqCteXV0Nq9WKkZER/zvxYLfbodVqUVpa6hHlD4ZhRwhBV1cXTCaTR5kpdQeLi2d44lqtFk1NTXC5XGwk3Z+AJB00pE6dUZKQt/p6Cv7DTCsEm5qa2A48qampHiw3l8vloQcgJXp6ejA9PT0n5mB1uLDz6AjOW56CzARPo+JX7XGrHO12O0vTHR4eRnV1NRtf8kYMEqIVQM9JB52ioiJs374d9fX1WL9+PXbs2MEOLL4QQIXdeQC6CCF98x5YAEKeshMKqiRjs9lYV9zhcIjm0NOAXVxcnGQGTwNrkZGRqKiomLdaLC8vD3l5eawLTUUqkpKSkJqaiqSkJI+HaWBgAGNjY6Ii3kJAhTbFeg7eKgQptz4+Ph5JSUkYHh5GZmampDRgOqharVaUl5fPMbgDTRpMWhz4nzX+z8m9B5fLhd7eXlbnoK2tbQ6vwZtikNAmIvR8kZGR+NGPfoR169YJmuzE9rEDcBWAF3nvie5jBywSo6d88uTkZCxfvtxDoVYMh57q2FdWVqKjo4N9nzuSiw000WvLzMwUFVjju9BcxVpKqJmenobdbg9IOHQ+jI+Po6OjIyChTS74HXjo7CmTyTA2NgaXy4XU1NSgvRPKsHS5XF6zFYQQ/OfwIIpSY7AmL1HUsSkTcP369VAqlR68BhqLoQFNb1oBvpqI8MGN3tPPJexjpwKwEcBPOW8H1McOWCD3nqrSeHuwDQYDGhoaUFRUNIenLrTKju8lcKO3wazfjUYjGhsb53WPhYAfSTcYDGwkPTo6GgMDA3PSgYGCsgKlJt1YrVZ0dXWhtLQUycnJrPY+14WmrDqxy6aWlhbI5XK2uQcf9UPTaB4x4v4N4ioKaaBv1apV7HfB5TVQb6y/v5/1ZGiVIL+LEPdZomXC3GCgNyVcKSrsZnExgCOE07uOBNjHDligmd6XGMbY2Bi6urpQUVHhNccpxOgpLTchIYH1EuggE4zB09myrKxM0vyr0+lER0cHsrOzkZOT4zUdmJqaioSEBNHXzM2VS8kKpN2JuMIkfO398fFxtkKQX2LrCzQYGBkZiaVLl/q83xcODyE2Qo7LKoRHx7VaLcsM9DX4cb0xruQZN0dPB2OaEqSqTJTcQ4OB/f39ogZZIRV2HFwNnmvPBNjHDlhgo6cPAF2/TU5OzqtQ68/ozWYzamtrvdJybTYbpqenRXVhoRgcHMTIyMi8VNJAQPXxCgoKWK+Gnw7kdrGNj49nKanzpdKo4q/VapU8V05jA+Xl5T4HP7lc7lHswi2x5eryc/PYbrcbDQ0NbO8+X9AZ7XirWYNtq7MQoxL2uHINXmichBYw0Ry91WpliU0WiwVJSUlITExEX18fCgsL2ZSq2+3Gs88+C5lMJmpyoBV2//jHPwDgAnivsAPDMNGzn/P71AXUxw4AGD9Bh4CLj6k0NAA29xwbG8v2b4+KivKb6iGE4ODBg1i/fv2cz6ieni+GnUajwejoKMxmM9RqNVJTU/1WpVHjsVgsKC0tlZSiSo1HqIwXnXm0Wi3Gx8d9ptLobKlUKiUV1ACOFc5UVlYGvG6nxqPT6WC1Wtl8en9/P6s1MB+e/m8v/vxhD/bdehoKUvwvfzQaDXp7e0UZvD+43W7odDq0trayrMHExES43W4cPHgQr7zyCl5//fVglmfS/WgCsKAzPeW+8/u3+4KvB3hgYADDw8NzGHbcGngagHK73dDr9RgdHUVbWxvi4uKQlpY2Z/akajQxMTGSE01oNZuYwBp35vGVDlSr1eju7oZarZY8dSakcEYIIiMjPTwZnU7H6hYoFAoolUqfyjdOtxsvfzmMMwqTBBn82NgYW+UoZSbE7Xajv78fJSUlSE9Ph8lkQktLC370ox+hv78fN9xwAzo7O1FRUSHZOUOJBTP6iYkJDA0NzauU4w80YGe321FTUyOIYUcFKGhVGpdMExkZyeagW1tbPRpESIWhoSGMjIwEHVjjpwNHR0dRW1sLmUyG6OhojI+Pz0kHBopACmeEgBpPcXExMjIy2PoGXxWCH7TpMGaw4f4Ny/wemxp8VVWVpNfscrlQW1uL7OxslnEXExODgYEBxMfHo6mpCYcOHcLo6OgJY/Qhc+9pswlCCL766itYrVbU1NSIZqp99tlnWL9+PctJT0pKYgUbgeAYdiaTCYODgxgaGkJUVBSWLFkiWRSdS+YpKyuTdKlA6bpFRUVQq9VsOnBiYkJwdaAv0GBgZWWlpNdMmYH5+fle1YRohaBWq2VptQ9/Og2N0YW3fng65DLfv+vo6CgGBwc9KLtSgBr8kiVLPKoy33zzTfzpT3/Cvn37gmqVzsHJ497T9abdbkdBQUHAwhK0cQU/YBdsSs5isWBiYgKnnXYaFAoFG0WnKSg+E00ouOvs+cg8gcBgMKCxsdEjNsAvrPFXHegNoSycocKlS5cuZfX/+eBXCH7VOYKjQyO4sliBxoZ6D9ltLkZGRjA0NBQSg6f8DK7Bv/322/jDH/4gpcEvOEI205vNZnz11VfsqK5QKAKqGvv4448hl8tRXl7uwSwLhmEHHGPCVVRUzHG7KRNNq9WyTR7S0tIEuc80hZiamippuS1wLLBWXl4uqKKLpgO1Wu286UBu4YyvXHmgsFqtqK2tFa2/9+C+NrxWN4oPbl8HJTnmBRBC2GXA9PQ0W5QjpVdCDT49Pd0j9vTBBx/gwQcfxP79+yUtiMICz/QhM3qNRgOHw4GUlBQMDQ3B4XCIDjb19/ejvb0d69ev95ipghGtpOwvh8OBlStX+t2fFnNoNBpMTEywKrUpKSlzZhar1cqKdEjZsgo41sWmsrIyoMAaTQdqtVpMT0+z6cCkpCS0tbWx7ZulNHjK/V++fLmoOI7B6sTX/vgZvlGaikc2rvD4jEuoMZlMyMjICKhC0Bfcbjfq6uqQmprqMUl9/PHHuO+++7Bv3755VZsCxMnh3qekpLCSVXK5HFarVfC+breb7REWHx8vSdMJ4JiwR1xcnOD0FreYg6tS29fXB6VSybrPDocDjY2NWLFiRcCBSl8YGBiARqMJKirNz6VPTU1Bo9GgqakJKpUKeXl5sNlskmn7GY1GNDQ0BFQ1+FrdKCwOF65ZM9czpE00lUolzjrrLBgMBpZ5F0yVIzDz3FHZcq7Bf/bZZ7j33nvxxhtvhMLgFxwLmrITAm7AbsWKFTh69KgklFqr1Yr6+nrk5OR4rNHEgK9Sy+35ZjabkZ2dDYVCIZnMMg0Gms1mSdtZ01xzV1cXiouLoVarA6oO9AUad5iP0OMLbkLw4peDqMqOx8rMubJog4OD0Gg0bKCRW11Hu+/QCkEaBxASl6FkIbVa7dHo4/Dhw7jrrrvwxhtvSFpgdDyxIEYvtIklDdgtXbqUdY9pbXwwBk+JMVLPwlFRUVAoFFAoFDj99NMxNTWFrq4uWCwWqNVqpKWlBUSnBY55OzKZTHLeAI2kcxtbiK0O9IWpqSm2eWggWZCD3RPoHbfgN1vmsvQGBgbmzSzExMQgJiaGvQ9u9x26nFGr1XOWZbSVV2Jiokcc5ujRo7j99tuxZ88erx1/TlSEtOCGQshMTxl2/IAdVbYN1OA1Gg26u7uDrjjjg0a7DQYDqwsXExODJUuWsAqvlE4rVpnG5XKhoaEBCQkJkvavB46l+3wVEfmrDpwvHTgxMYG2tragGHwvHB5EcowKF670DJT19/dDr9cL1vbjd9/h8jO4uvuRkZFoampCXFwc8vLy2P0bGhpwyy23YNeuXfPShE9EhCyQRwhhVW/MZjOrg+4N/f39GBkZQVVVlUeQihCCwcFB9PX1ISEhAWlpaYLlqQkhbM65oqJCcpIJrQwTQiWmqi56vR5RUVFsHMDbNdHljdT16sCxTjZiA2sAPNKBOp1uTjqQrqv5v6EYDE5YcNETn+Pms/Lww68Xsu/39fVhYmICFRUVkixx6LJMp9NhenoaMTExKC4uZr2y5uZm3HDDDXjllVewfPnyoM8nACdH9J5r9DabDY2NjVi9erXHNtSFdTqdc7juXNFKAGwEXa/XIzY2FmlpaUhJSfE6c1LmHk1BSZlzdjqdqK+vh1qtRl5enugyUirxrNVqPeSfo6Ki2Og/t6BDKggpnBEDbjrQZDKx33VKSkrAnsnv3u3Ec58P4r3b1yE9fmbg6O3txdTUlFdhjWBACEFzczPbUFSr1eK///0vDhw4gO7ubvznP//B6aefLtn5/ODkMHoAbA87p9OJr776Cqeddhr7mcPhQG1tLZKTk1FQUCCYYUcj6BqNBjqdDpGRkUhLS2NnTlrQk5iYKLlrTIOBXL25YI/HzaPbbDYsW7aM1lhLcMUzoPn9iooKybuw0Pr9nJwc6PV6j3Sgt/WzL1gdLnz9T5/h9IIk/PGKGampnp4eGAwGlJWVSW7wLS0tUKlUHiW97e3tuPXWW1FWVob6+nps374dV111lWTnnQcnR8qOC1+davjCGUIi9NwIelFREUwmEzQaDY4ePQqGYWC1WlFQUBCUfLQ3UEGNkpISyZhYtGttXFwcmpubkZ+fz9ZzJyUlIS0tLeh+dXQdK3WZMHBMsZbW73Pr0sWKhe5v0mDK4sT/1Mwsabq7u2E0GkNi8FTJmGvwfX19uO666/CPf/wDNTU17LaB4q233sLtt98Ol8uF733ve7RunsVHH32ETZs2oaCgAHV1dbUAdhNCHvJ6MIkRUqOn6jlc46UtqvmdagJl2MXExKCgoABqtRqNjY1IT0/H2NgYK0+clpYWtB65WCacGNBqturqatYoaACNtnqKi4tj6+rFUE25raukjGkAvhVrhVQH8tOBhBD859AgitNiUJOXyKYpQ2HwbW1tkMlkHkSkwcFBXH311fjf//1f1uDpvQQCl8uF2267De+++y6ys7OxZs0abNy4cY6U+1lnnYU333wTmFW8XSgsqEZeX18fRkdH57So5hp8ID/y2NgYent7UV1dzUaNaQFHR0cHS0FNS0sT3dxhZGQEAwMDIZkpBwcH2ZmSa5R8ea3p6WlWGILOnKmpqfNejxQdZ3zBl2KtNwgRC+0zydAyasQvNixjxTHLysokXeJQJiYAD2LWyMgItm3bhj//+c+SreEPHTqEoqIiVl3nqquuwt69e0X3bwgVFsToaZOBqakprFmzZk6bIW9dZoSAEIK+vj7o9fo5hsMt4KBcetrcgebQ59NzI4Sgt7cXk5OTHpLXUoA2vJyensaqVavmPTbDHOv7TpczWq2WbcFMBwDqgYSycIaShSwWS0CBNV/pwCc+HkO0gkEeNLBYFCEx+M7OTrjdbg/h1bGxMVx55ZX4/e9/j7PPPluy8w0NDXnk9bOzs/HFF1/M2e7gwYOorKxEfX39AQB3EUKaJLuIeRBy995ms6Gurg4ymczDXQuWYUfTZjKZzO/DzVV1paIaVM/NWyqQZhUYhpG85xu3uCWQFBQloOTn57P6/u3t7aw3YzaboVAoJCf0cBVrpTBK6s24I2LxlWYUFxVGgnHZYbE48NVXXwmuDhRy3V1dXXA4HB7FRDqdDldeeSV+9atf4dxzzw3qHN7OyQf/+6qurkZfXx9d5jwB4DXw+tSHCiE1eqPRiNraWhQVFaG7u5v9MoI1eFrJRuWWxOzPF9WgqcCOjg5Wmnp0dJRt5CCl4VB1ntjYWA9NgEChUqnYhgs0v2+328EwDFpbW9kIerCDFg1+yWQyyavwXv1qGE43wTeWRmPt2nJ2opBKLLS7uxs2mw0rV65k99Xr9bjyyivx4IMP4qKLLpLsXiiys7MxMDDAvh4cHJwjzsKrGN3PMMxfGYZJIYToJL8gHkKasmtqakJ6ejri4uLw5ZdfsiSZYEpiaeVWYWGhVzGGQEEIgV6vR1NTE2QyGWJiYthUoBRS0nSgSktLk5zSSRl8NE3JbfOk1+uDEtYQqlgbCOxOF8794yfIS1Dg3zeu93psX9WBQtKB3d3dMJvNHlr6k5OTuPzyy3H33Xdjy5Ytkt0LF06nE8uWLcP777+PrKwsrFmzBi+88AJKS0vZbUZHR5Genk5t4DQAOwHkkWBSBgIR0pm+pKTEo0uIw+Fge9sFYvCTk5NoaWmRvN8bMMNWo5LXarWaXTvTpUlqairS0tICqt6ipBuuCq5UcDqdbO03TVPyKwONRiOb1lQoFKzr7O9eKCed6sRLCUII/v1hPfQWN365yTer0Vt1oJB0YG9vL6taRI89PT2NrVu34o477giZwQMzy8knn3wSF110EVwuF2644QaUlpbi6aefBgDcfPPN2LlzJ5566ik6cP0ZwFULYfBAiGd6viJuRkYGkpOTA5otKAmkoqJCsvJPCsoZ96VxT0k0Go2G7eoiNBVIqa9S5vcpvBXO+AOloGq1WjaFRu+F+7vQfnWUeSgl6HLh/o/0mHLIcOAH88th+QLtWce9l9TUVIyPj2N6etojhmQ0GrF161bceOONuOaaayS9Hwlw8jDynE4nWyGn1+vR19cHu92OlJQUpKenz3nQvF7AbKSbUjGllEQCxA8mDoeDHQCsViv7oHkr36SdaMvKyuZ0PwkWVJEmmO47NIWm0WhYqXCa1uR7D1KBsuEGjQTb39LgrvOX4ob1wSsMcdtvWywWZGRksBkam82Gbdu24dprr8V3vvMdCe5CcpxcjDwasEtKSoJarYbT6WRbClssFnam8ZY/p+tJhUKBqqoqyaPR3HSf0MFEqVSyqUAq6UzLN7ksOr1ej87OTsk70QLHOs4EWyrMTaFxKwPHxsYQHx8PlUoFl8slWbqSy3f/TONChEKGb64KTNuAD0rBjoyMRE1NDbsMuPrqq6HVanHmmWdi48aNkpzrREdIZ/q77roLmZmZ2LRpk1c+OTWasbExmEwmlkCTkJDgEfiSWmuOupdSFuTQvLNGo2G13JYtW4a0tLSQNM2QqnCGC24v+8jISLbRBq1v4HZ6FQtq8BEREUjNysW5fzqIDWXpePgyaarYBgcHodVqPVKsNpsN11xzDaqqqhAdHY3PP/8cb7zxhqSTh0Q4edz73t5e7Nq1C6+99hoA4NJLL8XmzZuRnZ3tdQDQ6/UYGxvD1NQUHA4H8vLyQpI247ZTCpX3QHn0er0e0dHRrNEEszwJZeHMfIq1tL5Bp9OBYRjROXRCCJqamhAVFYXCwkI8/8UgHnunEztvrPGqjiMW1DvhimvY7XZcf/31OPvss3HHHXdI8jv749NTHD58GKeffjpefvllXHHFFUIOffIYPXsQQjA8PIxdu3Zhz549sFqtuPTSS9mCA+4Potfr0dbWhuzsbBgMBkxNTSEhIQHp6elBN3OgbadD0dSCklecTqeH98CNnut0OlbAIS0tTdSsSSPWwdSr+4IYxVqaQ6fCp/7ozW63G01NTYiOjsbSpUvhJgQb/vIFUmJU+Pd3qoO+9uHhYVaLgRq80+nEDTfcgJqaGvzkJz+RxOBdLheWLVvmwad/8cUX51BrXS4XLrjgAkRGRuKGG244dY3e44CEQKPRYM+ePdi1axempqawYcMGbNq0CZ9++iny8/NxxhlnsA82X402Pj6ebUslZgCgUfRg2057A01tRUVF+VWUNZvN7BKAYRiWCzDfup8WzkjdvYVeTyCKtcAxqXCNRgOj0ThHWstbyu/jjnHc/GI9fvfNldhQFpxi8MjICIaHhz0M3uVy4aabbkJJSQnuv/9+yTy5gwcP4oEHHsDbb78NAPjVr34FAPjpT3/qsd2f/vQnKJVKHD58GJdeeumiNPoFLbgBZuiI6enpuPnmm3HzzTdDp9Phtddew9atW8EwDLZs2YKMjAyW+cXPOU9NTWFsbAydnZ1+xTQoaH4/FFF0KqohpBkjMFOAkp+fj/z8fNhsNmg0GrS0tLBCjvz0WSgLZ2hAMFDeA5/ezJXWio2NhcViYfUSAMDucuN373UhMyEC568ITiRkbGyMbXLBNfgf/vCHKCgokNTgAWF8+qGhIezZswcffPABDh8+LNm5pcaCGz0fKSkpMJvN2LJlC+655x68+eabePjhh9Hf348LLrgAmzdvZjnq3LJNWn1GS1OjoqKQnp4+Z91MK/BWrVoleX6f1hUEKqoRERHB9ninKSea1UhOTobdbofT6ZS8cAYITrHWG7iVgS6XC0ePHoVcLodOp8PU1BTS0tLweqcNnVoT/npVOVTywO+H26iS/tZutxt33HEHUlJS8PDDD0serBPCp9++fTt+/etfSz44S40Fd++9wVtaaHp6Gvv27cOuXbvQ0dGBc889F5s3b8bq1avnGIC3dXN6ejqsVismJydRXl4eMrdYbOcWIaDeg9lshlwul0xQg4LyByoqKiTXB6Da8UlJSSypx2w2o7ZrGLfsHUBVmhwPXJATsM6BRqNBX1/fHIO/++67IZfL8fjjj0s+QALC3PuCggJ2cNDpdIiOjsYzzzyDzZs3+zv8yb2mDwRmsxn79+/H7t270dDQgHPOOQebN2/Gaaed5nVUNRqNaG5uhtlsZmMAYgNn84GmzUJBB6bcBNpxhhDCpgInJydZ7jm/1bZQSKFYO9+119XVITk52WOpQwjBDc/XomnEgN3fWwWZzcCSm8QU01BNAW7TD7fbjZ///OewWCz461//GhKDB4Tx6bm4/vrrw2v6YBAdHY0rrrgCV1xxBaxWK959913s2LED27dvxxlnnIEtW7Zg/fr1UCgUcDgc6OrqQnJyMtasWQOr1YqxsTGWQ08DZ4G6+uPj4+jo6JBcUhuYWzgDzLiQXEENLvc8KiqKjWkI8WTGx8dZwpDUSx1K201JSZlTUPR6/Ri+6J3E/RuWITs5DkAcS24aHx8XJBVOxTe4Bk8IwcMPP4zJyUn8/e9/D5nBA8L49CcKToiZ3hfsdjs++OAD7Ny5E59//jmqqqpQX1+Pxx9/3EOEk8JqtUKj0UCj0YAQwnoAQmc8qqJTVVUlmddA4a1wZj5QZV26pFEoFOyA5i2lp9Fo0NvbG5Jrpw0f09LS5lz7hNmOS/5yCPnJUfj3d6oh8zGbu91uttUWn9swPT3NymvTayeE4LHHHkNvby+effbZRb+O9oOwex8I2tracNlll2HlypXo7OxEdXU1Nm3ahHPPPderEXDzzU6n028RTV9fH8bHx1FRUSE5/z+Qwhk+LBaLBxuQO6CNjo6yg5XUsQ1fHV4p7t3bgjcbxrDr+zUoThMWMOQOaKOjo6zYaUZGBqKiokAIwR/+8Ac0NjbiP//5j+S/x3FA2OgDwTvvvIPU1FSsWrUKLpcLn3zyCXbt2oUPP/wQZWVl2LRpE84//3yvLjktohkbG5tTEASA1dkrLS2V3IWUonCGD+6AZjabAQBlZWUBt9jyBZfLhdraWmRmZnolO33eM4Ebnq/FjWfk4o7zloo+/sTEBNrb27F8+XI2U/PEE0/AYrHAbDbjwIEDknstxwlho5cSbrcbX3zxBXbu3Il3330XxcXF2LJlCy688EKvqSpaEKTRaGCxWOB2uxEXFxcSg5eqcMYXaLfbjIwM6HQ6SXrsUdDliC+Dtzld2PT0YRAC7L15DSKV4tzvyclJtLa2egiSEkLwxBNPsO2iOzo68PHHH0seTD0OCBt9qOB2u3H06FG8+uqreOutt5CXl4dNmzbh4osvRkJCgse2NG2mVCpZd5NbEBTsjEkzAKEgDAHeO8PwVWgSExORlpYmmt7sdDpRW1uLrKwsnx2Af/1OJ577fAB/v7YS6wvFpTRpE0xuwJEQgmeffRZvvvkm9uzZg8jISNjt9qBmen9c+r179+LnP/85ZDIZFAoF/vSnP+HMM88M+HzzIGz0CwFKEX311Vexf/9+pKenY9OmTbjkkktgs9lw9OhRrFq1in2ouQVB3DLapKQk0QMATZuFonCGquFSmShfxsynN8fFxbH05vmCYk6nE0ePHkVOTo5PQtKh3gl8Z0ctttUswf0bSkRdP7frLTfA+vzzz2Pnzp14/fXXJUk1CuHSG41Glh1ZX1+PrVu3orW1Nehze0E4ZbcQkMlkqKioQEVFBR566CG0tLRg586duPTSS6HRaHDNNddgzZo1rLQXt+8cpZyOjo6ira1NVEFQKDvOEELQ0dEBp9PpV7GWT2+ma+b5UoG0FVlubq7PgKPB6sRP97YgVx2Fu84vEnX909PTXg3+5Zdfxosvvoh9+/ZJxi0Qok3PXf6ZTKbFWJIbEE5Zo+eCYRisXLkS119/PXbv3o0nn3wSbW1tuPrqqxEZGYnLLrsMmzZtYoUMuZRT7ozZ3t4+b0FQKDvOUI0AhmFEK9Z609bna+olJSWhpaVlXoMHgEff6oBm2o5/f2cVolXC1/EGg4FlCXINe/fu3fjXv/6Fffv2ScoeFKpNv2fPHvz0pz+FRqPBvn37JDv/8UTY6DnIysrCm2++yeaa77nnHvT19WHXrl349re/DZlMxmoCZGVliSoIGhoaClnhDJ/FF8yMxDAMYmNjWZlui8WC0dFRHDp0CJGRkWzk3Nuy5EDTGPbWj+KWs/JQmZ3g5ejeQfsE8pc7b775Jp566ins27dP8riHEC49AGzZsgVbtmzBxx9/jJ///Od47733JL2O44Gw0XMgl8s9yCUMwyA/Px8//vGPceedd7KaADfddBNsNhurCUCFPvgFQQaDAWNjY+wMTGm1UoLGJqiRSg25XA6tVovy8nK2pXNbWxub2kxLS0NsbCyGp6z4xZttqMyKx81n5ws+vtFoRENDw5w6gLfffht/+MMfsG/fvpBkNoRo03Nx9tlno6urCzqdbo7IyImGUzaQFwyoJsDu3buxe/duTE1N4ZJLLsGmTZtQXFzs0Zixra0NLpcLubm5cwqCaHvtQBFKxVrgGGmooKAAqamepbBOp5MV1Zw2mvD7o24MGpzY/f01yFELC07SlCW/0u/999/Hww8/jH379s05r1QQwqXv7Oxktf6PHDmCyy67DIODg6FY24ej9ycaqCbA7t27odFocPHFF2PDhg144YUXcN11181ZY9M1s1arZemzYguC5qO+SgG73Y6jR496lc/i47fvduBfBwexfW08yhPsggKbZrMZdXV1cwz+448/xs9+9jO8+eabAZUri8H+/fuxfft2lkt/3333eXDpf/3rX2PHjh1QKpWIiorCb3/723DKzhv0ej22bduG3t5e5Ofn45VXXpmj997W1oZt27axr7u7u/HQQw9h+/bteOCBB/C3v/2NHeEfffRRbNiwQexlHDdMTExg165d+MUvfoGsrCx8/etfx+bNm302fLRYLBgbG4NWqxXcVENInjwYUL28oqIivyzBj9p1uPWlBmxdvQQPXFIyJxXoTeiEGjyfo/DZZ5/h7rvvxptvvim5nNkix4lt9HfffTfUajXuuecePPbYY5iYmMCvf/1rn9u7XC5kZWXhiy++QF5eHh544AHExsbirrvuEnvqRYOf/OQnKC4uxtatWz00Ac477zxs3rwZ1dXVXgcAIQVBDocDR48eDYqnPx+owRcXF/vVCRicsOCKv32JJYmRePGGakQoPAOUNK5BlzWRkZFITEzEyMjInLLkQ4cO4Y477sDrr78ueduvEwAnttGXlJTgo48+QmZmJkZGRvC1r30NbW1tPrd/55138OCDD+LTTz8FgJPC6N1u9xyjppoAO3fuRHNzM6sJsHbtWq/RfLvdzg4AtCAoKSkJbW1tKCwsDMlaV5RAptOFa/51BAN6K169sQa5av/5c71ej4aGBqhUKqhUKrZVlU6nw2233Ya9e/eyJcWnGBbU6CUvQB4bG2NdzszMTGg0mnm3f+mll3D11Vd7vPfkk0+ioqICN9xwAyYmJqS+xJDD2yxONQFeeuklHDp0CBdccAGeffZZrFu3DnfeeSc+/vhjOJ1OdnuVSoXs7GxUV1dj1apVkMlkOHLkCJxOJ6anp2E0GiXNBFCDLykp8WvwhBA8tL8dzSNGPLZ5hSCDt1qtaG9vR2VlJdatW4fS0lKYTCZcffXV2Lx5My688EI4HA6pbieMeRCQ0Z9//vkoKyub87d3715Rx7Hb7Xj99ddx5ZVXsu/dcsst6OrqYqu3fvzjHwdyiYsakZGR2LhxI55//nkcOXIEmzZtwquvvor169fjRz/6ET744AMPA3A4HBgZGUF1dTVOO+00REdHo6urC1988QU6OzsxPT0d1ABADX758uWC+u29+OUQ9tTO5OO/XuI/fUWXDFzF3cjISLjdbrhcLrzxxhtYuXIlduzYEfA9ULz11lsoKSlBUVERHnvssTmf/+c//2GZmOvXr0ddXV3Q5zzRcFzd+7179+Ivf/kL3nnnHY/3aTCws7MTWq0WAwMDXh/G/Px8xMXFQS6XQ6FQ4Msvv/TYf75g4mKEw+HA//3f/2Hnzp345JNPsHr1aqxevRoff/wxnnzyyTnVZLRDEJWgDqQgyGKxoK6uTpQE9jstGrzTrMVvvrnSpygGBTX4ZcuWefwGbW1tuO666/DCCy+grKxM0Hn9QQif/rPPPsOKFSuQlJSEAwcO4IEHHvDKxFtgnNhr+v/3//4fkpOT2UCeXq/Hb37zG6/bXnXVVbjooos8mgqOjIzgj3/8I9RqNSIiIvDss8/iG9/4htdgYH5+Pr788ss5KSWxwcTFCJfLheeeew733HMPCgoKsHTpUlYTwBv/PJCCIGrwK1asmFNlKAVo2o8fFOzq6sI111yDHTt2oKqqSrLzCdWmp5iYmEBZWRmGhoYku4YAsbCkfkLIfH+iodPpyLnnnkuKiorIueeeS8bHxwkhhAwNDZGLL76Y3c5kMhG1Wk0mJyc99r/22muJSqUiy5cvJ5dddhk5evQoWbZsmddz5eXlEa1WO+f9ZcuWkeHhYUIIIcPDwz73X+y4+eabSXt7O3G5XOTTTz8ld9xxBykvLyeXX345+fe//000Gg0xmUxz/gwGA+nv7yeHDx8m7733Hjl8+DDp7+8nBoOB3Uar1ZL33nuPDA8Pez1GsH8TExPkww8/JP39/R7vNzc3k4qKCnL48GHJv69XX32VfPe732Vf79ixg9x2220+t//tb3/rsf1xhD87lPRvUZJzEhMTMTk5yb5OSkryGtArKChgZ7KbbroJ3//+90XtfyLC7XbjyJEjePXVV/H2228jPz+f1QTwJiZBOGq6tENQQkIC+vv7UVZWFhIBCppWLCws9PDCBgcHsXXrVjz99NM4/fTTJT8v/U7+/ve/A5gpxz106BCeeOKJOdt++OGHuPXWW/HJJ59I3vEoAJwapbXnn38+RkdH57z/yCOPCD7Gp59+iiVLlkCj0eCCCy7A8uXLcfbZZ0t5mYsOMpkMNTU1qKmpwa9+9Ss0NDRg586duOSSS5CRkcFqAtD1M8MwHgVBtBZAqVSir69PUIcgMaDltwUFBR4GPzIygm3btuGJJ54IicEDwvn09fX1+N73vocDBw4sBoNfePhxBY4LAnHPf/GLX5Df/va3Ae9/osPtdpOmpiby4IMPkjVr1pALL7yQ/OUvfyF9fX2sa63RaMh7771HRkdHidFoJCMjI6S2tpa8//775LPPPiOdnZ1kamoqYJd+amqKfPTRR6Snp8fj/e7ubrJq1Sry/vvvh/Q7cDgcpKCggHR3dxObzUYqKipIY2OjxzZ9fX1k6dKl5NNPPw3ptYjEgrr3oRMKDwIbN27Ec889BwB47rnnsGnTpjnbmEwmGAwG9v/vvPMOGwWm++v1epx99tmsJ+DNxR8YGMDXv/51rFixAqWlpXj88cfZzx544AFkZWWhqqoKVVVV2L9/fyhuVxJQTYD7778fX3zxBZ588klMTU1h27ZtuPTSS/HII4/gyiuvRGlpKdthNj4+HsXFxTjttNNQWFgIs9mMr776CkePHsXw8LCovDmlBufm5iItLY19X6fT4corr8Rjjz2Gc889NxS3zoKrTb9ixQps3bqV1aannPqHHnoI4+PjuPXWW1FVVYWampqQXtNixKJc04+Pj2Pr1q3o7+9Hbm4uXn31VajVagwPD+N73/se9u/fj+7ubmzZsgXAzAP3P//zP7jvvvs89j9y5AjUajUOHz6MZ555xmsUf2RkhM2BGwwGrF69Gq+99hpWrlx5UrADCSF4++238d3vfhdFRTNKNlQTYMmSJV4j+2ILgmjvuuzsbI8iGb1ej8svvxz3338/LrnkktDc4MmBEztlt5gglhIMAJs2bcIPfvADXHDBBSeF0QPAww8/jK1bt2LZsmUYGhrCrl27sGfPHtjtdlYVKC8vz+sA4K8giMpgL1myxKP4Z3JyEpdffjl+8pOfCOnldqojbPRSQWwUv7e3F2effTYaGxsRHx+PBx54AM8++yzi4+NRU1OD3//+9ycEyUcIyGxQj2oCGAwGVhPAl/oOvyAoJSUFOp0OWVlZHgGz6elpXHHFFbj99ts92JZh+ETY6MVgvizAddddJ9jojUYjzjnnHNx333345je/CWCmjiAlJQUMw+DnP/85RkZG8M9//jMk93G8QTUBdu3aBa1Wiw0bNmDjxo0+9fYsFguOHj0KAKyOXkREBOLj47F161bceOONuOaaaxb6Nk5UnNjknMUEoVF8u91OLrzwQvL73//e57F6enpIaWlpSK5zsUGv15Nnn32WXHbZZWTVqlXkJz/5CTl48CBL7pmeniaffPIJaWtrIyaTiUxOTpKOjg6yefNmsmTJErJp0ybS2NhI3G738b6VEwXh6L1UEJIFIITgu9/9LlasWIE777yTff+tt97C0qVL2cKNPXv2eHDECSH40Y9+hKKiIlRUVODIkSMe+85X9LHYkZSUhOuuuw6vv/46PvroI1RUVLCqMffeey8uvvhiGI1GVrFHqVRCrVbDZDLh7rvvxtatW/HAAw9Ar9cHfS3+vsvW1lasW7cOERER+N3vfhf0+U4J+BkVTmgIoQT/97//JQBIeXk5qaysJJWVleT1118nhYWFZNOmTWTlypUkMjKSfO1rX2O9BkII2bdvH/nGN75B3G43OXjwIFm7di0hhBCn00kKCwtJV1cXmytuampa+JsPAfR6PVm7di1Zt24dKS8vJz/84Q/J+++/T3Q6HdmwYQN54oknJJ3dhXyXY2Nj5NChQ+Tee+9leRonIBZ0pj+p1XCTk5Px/vvvz3l/yZIlbM79zDPPBOHFNQ4ePIiioiK89tprAI4VbnCj03v37sW3v/1tMAyD008/HZOTkxgZGUFvb6/fJgonKtrb27Ft2zbceeedsFqteOedd/Cvf/0LBw4cwC233ILbbrtNUtFIIQ0paDrxZNGkXwic1EYfKIQ0QvC2zdDQkOAmCiciTjvtNJx22mkAjmkCbNy4EVNTU4iPj5dcJfZk/i6PJ8JG7wX8mR+Y2wjB1zZC9j3ZEIqyXEB4Q4owxCFs9F4gpHDD1zZ2u11UE4UwfENsQ4owhOGkjt4HijVr1qCjowM9PT2w2+146aWXsHHjRo9tNm7ciB07doAQgs8//xwJCQnIzMycs+/f/vY3PPPMMwHJN+Xn56O8vPyU5YgL+R3CCAB+In2nLPbt20eKi4tJYWEh+eUvf0kIIeSpp54iTz31FCFkpqrt1ltvJYWFhaSsrMxDFILuW1BQQJKSkuaNPn/66adEr9cTQgjZv38/mwUgxLdIyKkEf7/DyMgIycrKInFxcSQhIYFkZWWRqamp43nJgWBBo/cnPCNvMSNY+SZfcmBhnHQ4sSWwwzgGXxF+X/jHP/6Biy++mH3NMAwuvPBCrF69Gs8880xIrzWMUwfhQF4I4c2L8hV9/vDDD/GPf/wDn3zyCfveqagMFEboEZ7pQwix8k179+71kG+i26alpWHLli04dOhQ6C86jJMffhb9YQSBYOSbjEYjmZ6eZv+/fPlykpWVRZYuXUp+9atfzTnXhx9+SOLj41kq8YMPPsh+duDAAbJs2TKf+4Zx3LGggbyw0YcY/qLP3/3ud0liYiJrrKtXryaEENLV1UUqKipIRUUFWbFihd8swIcffkguueSSOedfTLUA/gYft9tNfvjDH5KlS5eS8vJy8tVXXx2HqzwuCBt9GJ747LPPyIUXXsi+fvTRR8mjjz7qsY0voxey70JAyODjq4jpFMCCGn14TX8CQGgW4ODBg6isrMTFF1+MpqYmUfuGGtziGZVKxRbPcOGriCkMibHQo0z4T/wfgCsB/J3z+lsAnuBtEw8gdvb/GwB0CN13ge7hCi/X8SRvmzcBnMl5/T6AmuP9/Z9sf+GZ/sTAIIAczutsAMPcDQgh04QQ4+z/9wNQMgyTImTfBYK3XCU/pylkmzCCRNjoTwwcBlDMMEwBwzAqAFcBeJ27AcMwGcwsCYBhmLWY+W3Hvex7I4DvMwzTyTDMPfwTMQzz/xiGqZ39a2QYxsUwjHr2s16GYRpmP/tS5D0IGXwWywB1cuN4uxrhP2F/mHHZ2wF0Abhv9r2bAdw8+/8fAGgCUAfgcwDrfew7DqAQgGp225XznPMyAB9wXvcCSAnw+hUAugEUcM5dytvmEgAHMDPjnw7g0PH+3k/GP3/c+zBOIjAMsw7AA4SQi2Zf/xQACCG/8rH9CwA+JIT8bfZ1L2bW2LoAz78BwJ8AyAH8kxDyCMMwN89ew9OznsqTAL4BwAzgO4QQsR5FGH4QNvpTCAzDXAHgG4SQ782+/haA0wghP/CybTRm3O0iQoh+9r0eABOYWWf/LyEkXBBwAiLMvT+1ICZQdhmAT6nBz+IMQsgwwzBpAN5lGKaVEPKx5FcZRkgRDuSdWhATKLsKwIvcNwghw7P/agDsAbA2BNcYRogRNvpTC36zAADAMEwCgHMA7OW8F8MwTBz9P4ALATQuyFWHISnC7v0pBEKIk2GYHwB4G8eCaU3cYNrsplsAvEMIMXF2TwewZzYrqADwAiHkrYW7+jCkQjiQF0YYpxjC7n0YYZxiCBt9GGGcYggbfRhhnGIIG30YYZxiCBt9GGGcYggbfRhhnGIIG30YYZxiCBt9GGGcYvj/yM7t2j618qgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "cloud = mean_center_normalize(cloud)\n",
    "ax.plot(cloud[:,0], cloud[:,1], cloud[:,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_edges(pointcloud, dmax):\n",
    "    \"\"\"\n",
    "    pointcloud: array of points (one point each along axis 1)\n",
    "    build list [(len, i0, j0,), (len, i1, j1), ...]) of indices belonging to points specifying edges, sorted by length.\n",
    "    the length is capped to dmax.\n",
    "    \"\"\" \n",
    "    # TODO: evaluate and use clever algorithm (space partitioning?) for not calculating all n*n distance squares\n",
    "    n = len(pointcloud)\n",
    "    dmat = linalg.norm(pointcloud[:,None,:] - pointcloud[None,:,:], axis = 2)\n",
    "    edges = []\n",
    "    for j in range(n):\n",
    "        for i in range(j):\n",
    "            if dmat[i,j] <= dmax:\n",
    "                edges.append((dmat[i,j], i,j))\n",
    "    edges.sort()\n",
    "    return edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "es = make_edges(cloud, dmax=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "es"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "# currently unused, homologies and rips complex are computed by the library dionysus\n",
    "\n",
    "class UnionFind:\n",
    "    \"\"\"Union-Find Disjoint Set structure with n entries in range(0, n)\"\"\"\n",
    "    def __init__(self, n, randomseed = \"thisisnotcryptographicallyimportantbutresultsshouldbereproducible\"):\n",
    "        self.n = n\n",
    "        self.r = random.Random(randomseed)\n",
    "        self.parent = list(range(n))\n",
    "        self.ncomponents = n\n",
    "    def find(self, i):\n",
    "        \"\"\"finds a representant of entry i\"\"\"\n",
    "        p = i\n",
    "        while p != self.parent[p]:\n",
    "            p = self.parent[p]\n",
    "        # path compression. non recursive as python has a very limited stack <3\n",
    "        while i != self.parent[i]:\n",
    "            nexti = self.parent[i]\n",
    "            self.parent[i] = p\n",
    "            i=nexti\n",
    "        return p\n",
    "    def unite(self, i, j):\n",
    "        \"\"\"unites the components of i and j.\n",
    "        return if an new connection was made\"\"\"\n",
    "        if random.randrange(2):\n",
    "            i,j = j,i\n",
    "        ii = self.find(i)\n",
    "        jj = self.find(j)\n",
    "        self.parent[ii] = jj\n",
    "        newconnection = ii != jj\n",
    "        if newconnection:\n",
    "            self.ncomponents -= 1\n",
    "        return newconnection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test\n",
    "uf = UnionFind(4)\n",
    "uf.unite(0, 1)\n",
    "uf.unite(3, 2)\n",
    "uf.unite(0, 3)\n",
    "assert uf.ncomponents == 1\n",
    "assert uf.find(0) == uf.find(1) == uf.find(2) == uf.find(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "# currently unused, homologies and rips complex are computed by the library dionysus\n",
    "\n",
    "def uf_label_edges(npoints, edges):\n",
    "    uf = UnionFind(npoints)\n",
    "    # if no components are connected, the edge produces a cycle.\n",
    "    # this is the default, as there can be much more edges than points.\n",
    "    labels = [1] * len(edges)\n",
    "    for (ei,(_d, i, j)) in enumerate(edges):\n",
    "        if uf.ncomponents == 1:\n",
    "            break # all done, additional edges keep adding cycles\n",
    "        negative = uf.unite(i,j) # reduction of components\n",
    "        if negative:\n",
    "            labels[ei] = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pcscore(pointcloud, dmax, n=2, m=2):\n",
    "    \"\"\"sw1pers-score for the pointcloud,\n",
    "    pointcloud should be an output of mean_center_normalize, as angles are used as distances.\n",
    "    \"\"\"\n",
    "    # disanzen als winkel\n",
    "    cdm = scipy.spatial.distance.pdist(pointcloud, 'cosine') # 1-cos metric\n",
    "    cdm = np.arccos(1-cdm) # angles\n",
    "    # TODO: laesst sich irgendwie scipy.spatial.KDTree.sparse_distance_matrix nutzen?\n",
    "    # aber wahrscheinlich egal, weil die homologieberechnung die komplexitaet dominiert\n",
    "    \n",
    "    fil = dionysus.fill_rips(cdm, k = 2, r = dmax)\n",
    "    hom = dionysus.homology_persistence(fil)\n",
    "    # TODO: schauen ob k=2 spezialisiert ist, mit algorithmus aus [0] vergleichen\n",
    "    \n",
    "    dgm = dionysus.init_diagrams(hom, fil)[1]\n",
    "    p = max(dgm, key=lambda p: p.death - p.birth)\n",
    "    # TODO: warum sind die exponenten so unsymmetrisch?\n",
    "    return (p.death**n-p.birth**m) / (3**(n/2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.428097320481896"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = np.linspace(0, 5 * math.tau)\n",
    "fs = np.sin(t)\n",
    "cloud = make_pointcloud(t, fs, dim=3, window_size=2, num_points=100)\n",
    "cloud = mean_center_normalize(cloud)\n",
    "\n",
    "pcscore(cloud, dmax=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uf.unite(0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uf.ncomponents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on built-in function fill_rips in module dionysus._dionysus:\n",
      "\n",
      "fill_rips(...) method of builtins.PyCapsule instance\n",
      "    fill_rips(data: array, k: int, r: float) -> dionysus._dionysus.Filtration\n",
      "    \n",
      "    returns (sorted) filtration filled with the k-skeleton of the clique complex built on the points at distance at most r from each other\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(dionysus.fill_rips)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "fil = dionysus.fill_rips(cloud, k=2, r=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = dionysus.homology_persistence(fil)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Reduced matrix with 2651 columns"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = dionysus.init_diagrams(h, fil)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD6CAYAAACvZ4z8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUM0lEQVR4nO3da4xc9XnH8e/DgoNDiCmxSfCtOAmB0nDVFlCCGiGVW/rChLYKaZRKTVLLUlGUF6AQRYqipm2qkKpqJFLXQlTqC2RFFaaWYjBppYoXhNaLcCCmNnEcUtabBEMCIcFcDE9fzCwMuzO7Z+7nnPl+pJVnzjn/mWeP1z+G/znPfyMzkSTV1wnjLkCSNFwGvSTVnEEvSTVn0EtSzRn0klRzBr0k1VyhoI+IayPiYEQciohb2+zfHBGPRsS+iJiJiCuKjpUkDVcsdx99REwBTwBXAbPAXuATmfl4yzHvAH6TmRkRFwDfzsxzi4xtZ/Xq1XnWWWf1/l1J0oR5+OGHn8nMNe32nVhg/KXAocw8DBARO4DNwBthnZm/bjn+FCCLjm3nrLPOYmZmpkBpkiSAiPhJp31Fpm7WAU+1PJ9tblv4Jh+LiAPAd4BPdzNWkjQ8RYI+2mxbNN+TmTsz81zgeuCr3YwFiIgtzfn9maNHjxYoS5JURJGgnwU2tDxfD8x1OjgzHwDeFxGruxmbmdszczozp9esaTvNJEnqQZGg3wucHRGbImIFcCOwq/WAiHh/RETz8SXACuDZImMlScO17MXYzDweETcBe4Ap4M7M3B8RW5v7twF/BPxZRLwKHAM+no3bedqOHdL3IklqY9nbK8dheno6vetGqod7HjnCbXsOMvfcMdaetpJbrjmH6y/2noxBi4iHM3O63b4it1dKUk/ueeQIX7z7MY69+hoAR547xhfvfgzAsB8hl0CQNDS37Tn4RsjPO/bqa9y252DPr5mZfOu/DvHY7PP9ljcxDHpJQzP33LGuti8nM/mH//ghX7/vIPfsO9JPaRPFoJc0NGtPW9nV9qXMh/w3//OHfHx6A1/66O/0W97EMOglDc0t15zDypOm3rJt5UlT3HLNOV29zsKQ/9oN53PCCe36MdWOF2MlDc38Bdd+7rox5Ptn0EsaqusvXtfzHTaG/GA4dSOplAz5wTHoJZWOIT9YBr2kUjHkB885ekk9G/TyBob8cBj0knoy6OUNDPnhcepGUk8GubyBIT9cBr2kngxqeQNDfvgMekk9GcTyBob8aBj0knrS7/IGhvzoeDFWUk/6Wd7AkB8tg15Sz3pZ3sCQHz2nbiSNjCE/Hga9pJEw5MfHqRtJbQ2y69WQHy+DXtIig+x6NeTHz6kbSYsMquvVkC8Hg17SIoPoejXky8Ogl7RIv12vhny5GPSSFumn69WQLx8vxkpapNeuV0O+nAx6SW112/VqyJeXUzeS+mbIl5tBL6kvhnz5OXUjTRg7XiePQS9NEDteJ5NTN9IEseN1Mhn00gSx43UyFQr6iLg2Ig5GxKGIuLXN/k9GxKPNrwcj4sKWfU9GxGMRsS8iZgZZvKTu2PE6mZYN+oiYAm4HrgPOAz4REectOOzHwEcy8wLgq8D2BfuvzMyLMnN6ADVL6pEdr5OpyMXYS4FDmXkYICJ2AJuBx+cPyMwHW45/CFg/yCIlDYYdr5OpSNCvA55qeT4LXLbE8Z8B7m15nsD9EZHAP2fmwk/7kkbIjtfJUyTo2/2NZtsDI66kEfRXtGz+cGbORcQZwHcj4kBmPtBm7BZgC8DGjRsLlCVp2Az5eihyMXYW2NDyfD0wt/CgiLgAuAPYnJnPzm/PzLnmn08DO2lMBS2Smdszczozp9esWVP8O5A0FIZ8fRQJ+r3A2RGxKSJWADcCu1oPiIiNwN3ApzLziZbtp0TEqfOPgauBHwyqeEnDYcjXy7JTN5l5PCJuAvYAU8Cdmbk/IrY2928Dvgy8C/hWRAAcb95h825gZ3PbicBdmXnfUL4TSQNZ3sCQr5/IbDvdPlbT09M5M+Mt91I3Fi5vAI1bJ792w/mFw96Qr66IeLjTLex2xko10e/yBoZ8fRn0Uk30s7yBIV9vBr1UE70ub2DI159BL9VEL8sbGPKTwfXopZrodnkDQ35yGPRSjRRd3sCQnyxO3UgTxpCfPAa9NEEM+cnk1I1UAXa8qh8GvVRyg/iF3ob8ZHPqRio5O17VL4NeKjk7XtUvg14qOTte1S+DXio5O17VLy/GSiVnx6v6ZdBLFWDHq/rh1I1UE4a8OjHopRow5LUUg16qOENey3GOXhqzfpY3MORVhEEvjVE/yxsY8irKqRtpjHpd3sCQVzcMemmMelnewJBXtwx6aYy6Xd7AkFcvDHppjLpZ3sCQV6+8GCuNUdHlDQx59cOgl8ZsueUNDHn1y6kbqcQMeQ2CQS+VlCGvQXHqRhqiXrteDXkNkkEvDUmvXa+GvAbNqRtpSHrpejXkNQwGvTQk3Xa9GvIaFoNeGpJuul4NeQ1ToaCPiGsj4mBEHIqIW9vs/2REPNr8ejAiLiw6Vqqrol2vhryGbdmLsRExBdwOXAXMAnsjYldmPt5y2I+Bj2TmLyPiOmA7cFnBsVItFel6NeQ1CkXuurkUOJSZhwEiYgewGXgjrDPzwZbjHwLWFx0r1dlSXa+GvEalyNTNOuCpluezzW2dfAa4t9uxEbElImYiYubo0aMFypKqy5DXKBUJ+nY/fdn2wIgraQT9F7odm5nbM3M6M6fXrFlToCypmgx5jVqRqZtZYEPL8/XA3MKDIuIC4A7gusx8tpuxUtXY8aoqKfKJfi9wdkRsiogVwI3ArtYDImIjcDfwqcx8opuxUtXMd7weee4YyZsdr/c8cmTJcYa8xmXZoM/M48BNwB7gf4FvZ+b+iNgaEVubh30ZeBfwrYjYFxEzS40dwvchjYwdr6qaQmvdZOZuYPeCbdtaHn8W+GzRsVKV2fGqqrEzVuqSHa+qGoNe6pIdr6oalymWumTHq6rGoJd6YMerqsSpG2mADHmVkUEvDYghr7Iy6KUBMORVZs7RSwt0u7yBIa+yM+ilFt3+Qm9DXlXg1I3UopvlDQx5VYVBL7UouryBIa8qMeilFkWWNzDkVTUGvdRiueUNDHlVkRdjpRZLLW9gyKuqDHppgXbLGxjyqjKnbqRlGPKqOoNeWoIhrzpw6kYTpZuuV0NedWHQa2J00/VqyKtOnLrRxCja9WrIq24Mek2MIl2vhrzqyKDXxFiu69WQV10Z9JoYS3W9GvKqMy/GamJ06nrdfNFaQ161ZtBroizsevWTvCaBUzeaWIa8JoVBr4lkyGuSOHWjWrDjVerMoFfl2fEqLc2pG1WeHa/S0gx6VZ4dr9LSDHpVnh2v0tIMelWeHa/S0rwYq8qz41VaWqGgj4hrgX8EpoA7MvPvFuw/F/gX4BLgS5n5jZZ9TwIvAK8BxzNzejClS2+y41XqbNmgj4gp4HbgKmAW2BsRuzLz8ZbDfgF8Dri+w8tcmZnP9FmrVIghL71VkTn6S4FDmXk4M18BdgCbWw/IzKczcy/w6hBqlAoz5KXFigT9OuCpluezzW1FJXB/RDwcEVs6HRQRWyJiJiJmjh492sXLSw2GvNRekTn6dv9Ssov3+HBmzkXEGcB3I+JAZj6w6AUztwPbAaanp7t5fdVckeUNDHmpsyKf6GeBDS3P1wNzRd8gM+eafz4N7KQxFSQVMr+8wZHnjpG8ubzBPY8ceeMYQ15aWpGg3wucHRGbImIFcCOwq8iLR8QpEXHq/GPgauAHvRarybPc8gaGvLS8ZaduMvN4RNwE7KFxe+Wdmbk/IrY292+LiPcAM8A7gdcj4vPAecBqYGdEzL/XXZl531C+E9XSUssbGPJSMYXuo8/M3cDuBdu2tTz+GY0pnYV+BVzYT4GabGtPW8mRNmF/5qqTDXmpIJdAUKm1W97g5BNP4Pz1qwx5qSCXQFCpLVze4MxVJ3P++lXs2f9zQ14qyKBX6c0vb+CcvNQbp25UCYa81DuDXqVnyEv9cepGY2PHqzQaBr3Gosgv9DbkpcFw6kZjYcerNDoGvcbCjldpdAx6jUWnX+htx6s0eAa9xsKOV2l0vBirsbDjVRodg15jY8erNBpO3WisDHlp+Ax6jY0hL42GQa+xMOSl0XGOXkOx1PIGhrw0Wga9Bm6p5Q02X7TWkJdGzKDXwHVa3uDr9x3g8DO/MeSlEXOOXgPXcXmD518y5KUxMOg1cJ2WNwAMeWkMDHoNXLvlDQAu33S6IS+NgXP0Grj5u2u+ft8B5p5/CWiE/F1/cbkhL42Bn+g1FJsvWssfT28AGtM1hrw0Pga9Bs775KVyMeg1UIa8VD7O0atrnbpeDXmpnAx6daVT12tm8uNnXzTkpRIy6NWVTl2vX/73/bzw8nFDXioh5+jVlU5dr4a8VF4GvbrSqev17SumDHmppAx6daVd1+tUBH99/QcNeamkDHp15fqL1/G3H/sgp76tcXnn7SumuO1PLuCGS9aPuTJJnRQK+oi4NiIORsShiLi1zf5zI+J7EfFyRNzczVhVy/zdNfNz8j/4yjWGvFRyywZ9REwBtwPXAecBn4iI8xYc9gvgc8A3ehirivA+eamainyivxQ4lJmHM/MVYAewufWAzHw6M/cCr3Y7VtVgyEvVVeQ++nXAUy3PZ4HLCr5+P2M1Yna8SvVUJOjb/YvOgq9feGxEbAG2AGzcuLHgy2tQ7HiV6qtI0M8CG1qerwfmCr5+4bGZuR3YDjA9PV30PyQaEDtepfoqMke/Fzg7IjZFxArgRmBXwdfvZ6xGyI5Xqb6W/USfmccj4iZgDzAF3JmZ+yNia3P/toh4DzADvBN4PSI+D5yXmb9qN3ZI34v6sPa0lRxpE/Z2vErVV2hRs8zcDexesG1by+Of0ZiWKTRW5XPLNee8ZY4e7HiV6sLOWAF2vEp15jLFAhZ3vDpdI9WHn+jlffJSzRn0E86Ql+rPoJ9ghrw0GZyjnxALlze4+eoP2PEqTQiDfgK0W97gln97lOOvpyEvTQCnbiZAu+UNjr+eNkNJE8KgnwCdljd48ZXXDHlpAhj0E6DTL/Re12G7pHox6CfAzVd/gBMXfHJfedIUt1xzzpgqkjRKXoytufmO1/k5+RdfeY11Lb9URFL9GfQ15n3yksCpm9oy5CXNM+hryJCX1Mqpmxpo7Xo9c9XJnL9+FXv2/9yQlwQY9JW3sOt17vmXmHv+JS7fdLohLwlw6qby2nW9AvzfL1405CUBBn3ldep6/enzL424EkllZdBX3JmrTm67vVM3rKTJY9BXWGZy/vpVi7bb9SqplUFfUfO3UO7Z/3Mu33Q6a1edTNBYv+ZrN5xv16ukN3jXTQV5n7ykbviJvmIMeUndMugrxJCX1AuDviIMeUm9co6+xOaXNjjy3DFOfduJvPDycUNeUtcM+pJauLTBCy8fZyqCy957uiEvqStO3ZRUu6UNXsvk7+9/YkwVSaoqg76kOi1t0Gm7JHVi0JdUpyUMXNpAUrcM+pK65ZpzWHnS1Fu2ubSBpF54Mbak5pcwmP+FImv9hd6SemTQl9j1F68z2CX1rdDUTURcGxEHI+JQRNzaZn9ExDeb+x+NiEta9j0ZEY9FxL6ImBlk8ZKk5S37iT4ipoDbgauAWWBvROzKzMdbDrsOOLv5dRnwT80/512Zmc8MrGpJUmFFPtFfChzKzMOZ+QqwA9i84JjNwL9mw0PAaRFx5oBrlST1oEjQrwOeank+29xW9JgE7o+IhyNiS6+FSpJ6U+RibLt+++zimA9n5lxEnAF8NyIOZOYDi96k8R+BLQAbN24sUJYkqYgiQT8LbGh5vh6YK3pMZs7/+XRE7KQxFbQo6DNzO7AdICKORsRPCn4PnawGqnBdoCp1grUOS1VqrUqdMJm1/nanHUWCfi9wdkRsAo4ANwJ/uuCYXcBNEbGDxkXY5zPzpxFxCnBCZr7QfHw18FfLvWFmrilQ15IiYiYzp/t9nWGrSp1grcNSlVqrUidY60LLBn1mHo+Im4A9wBRwZ2buj4itzf3bgN3AR4FDwIvAnzeHvxvYGRHz73VXZt438O9CktRRoYapzNxNI8xbt21reZzAX7YZdxi4sM8aJUl9qPNaN9vHXUBBVakTrHVYqlJrVeoEa32LaHwYlyTVVZ0/0UuSqGDQ97nuzpJjS1brSNcIKlDruRHxvYh4OSJu7mZsieos2zn9ZPPv/dGIeDAiLiw6tmS1lu28bm7WuS8iZiLiiqJjS1TnYM9pZlbmi8ZdPz8C3gusAL4PnLfgmI8C99Jo4roc+O+iY8tSa3Pfk8DqEp3XM4DfA/4GuLmbsWWos6Tn9EPAbzUfX1fyn9W2tZb0vL6DN6elLwAOlPRntW2dwzinVftE38+6O0XGlqXWUVu21sx8OjP3Aq92O7YkdY5akVofzMxfNp8+RKPRsNDYEtU6akVq/XU20xI4hTe79Mv2s9qpzoGrWtD3s+5OkbGDVKU1gvo5N6M8r/2+V5nP6Wdo/N9dL2P71U+tUMLzGhEfi4gDwHeAT3cztgR1woDPadV+8Ug/6+4UGTtII1kjaED6OTejPK/9vlcpz2lEXEkjPOfnaMv4s9o4cHGtUMLzmpk7aTRr/j7wVeAPio4dkH7qhAGf06p9ou9n3Z0iYwdpYGsEAfNrBA1LP+dmlOe1r/cq4zmNiAuAO4DNmflsN2MHqJ9aS3leW2p7AHhfRKzudmyf+qlz8Od0GBcihvVF4/9ADgObePMCx+8uOOYPeesFzv8pOrZEtZ4CnNry+EHg2nHW2nLsV3jrxdiRndc+6yzdOQU20lg25EO9fp8lqLWM5/X9vHmR8xIaa3RF2X5Wl6hz4Od0KH8Zw/yicafKEzSuaH+puW0rsLX5OGj8RqwfAY8B00uNLWOtNK7Uf7/5tb8ktb6HxqeUXwHPNR+/c9Tntdc6S3pO7wB+Cexrfs2U+Ge1ba0lPa9faNayD/gecMU4zmuvdQ7jnNoZK0k1V7U5eklSlwx6Sao5g16Sas6gl6SaM+glqeYMekmqOYNekmrOoJekmvt/DnM7GboYNMcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dionysus.plot.plot_diagram(ds[1], show = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Diagram with 50 points, Diagram with 20 points, Diagram with 1816 points]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(ds[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "'dionysus._dionysus.DiagramPoint' object is not iterable",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-72-1ee975404b16>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m: 'dionysus._dionysus.DiagramPoint' object is not iterable"
     ]
    }
   ],
   "source": [
    "a = list(ds[1])[0]\n",
    "list(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "npoints = 200\n",
    "R = 10\n",
    "r = 2\n",
    "def toruspoint(angles):\n",
    "    theta, phi = angles\n",
    "    d = R + r*np.cos(theta)\n",
    "    h = r*np.sin(theta)\n",
    "    return np.array([d*np.cos(phi), h, d*np.sin(phi)])\n",
    "\n",
    "angles = np.random.uniform(low = 0, high = math.tau, size = (2, npoints))\n",
    "points = toruspoint(angles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "fil = dionysus.fill_rips(points.transpose(), k = 3, r = 10)\n",
    "hom = dionysus.homology_persistence(fil)\n",
    "dgms = dionysus.init_diagrams(hom, fil)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWC0lEQVR4nO3df2zc9X3H8ecrJmlc2tQMPErs0lRdCTCgSWdoNaZB02oBEiAwpLGWoUGnKFLVZX80sEyjBfUP1uWfqKJtFEUVINAQWrOoZKUZapWylhLmKCEhDZnQoClJt5gSt0pwQ+K898edW+e4s793/t59f9zrIVn13X2wP9/Gefqbz31/KCIwM7Pim5X1BMzMLB0OuplZSTjoZmYl4aCbmZWEg25mVhJnZfWNzzvvvFiwYEFW397MrJB27tz5RkT013sts6AvWLCA4eHhrL69mVkhSfpZo9e85GJmVhIOuplZSTjoZmYl4aCbmZWEg25mVhKZHeViVmvLrkOs23aAw6NjzO/rZc3ShaxYPJD1tMwKw0G3XNiy6xBrN+9l7OQ4AIdGx1i7eS+Ao26WkJdcLBfWbTvw25hPGDs5zrptBzKakWXt6PG3eeCpffym5ufCGnPQLRcOj4419byV29Hjb/PZTTt4fMdB9v/i11lPpzAcdMuF+X29TT1v5TUR81dGjrHpziEWX3hO1lMqDAfdcmHN0oX0zu4547ne2T2sWbowoxlZFmpj/qcX1b1kiTWQOOiSeiTtkrR1ijFXShqXdFs607NusWLxAA/eejkDfb0IGOjr5cFbL/cbol3EMZ+5Zo5yWQ3sB+bVe1FSD/BVYFsK87IutGLxgAPepRzzdCTaQ5c0CCwDNk0x7AvAt4EjKczLzLqEY56epEsu64F7gNP1XpQ0ANwCbJjqi0haKWlY0vDIyEgz8zSzEnLM0zVt0CUtB45ExM4phq0H7o2IKQ8YjYiNETEUEUP9/f6DM+tmjnn6kqyhXw3cJOkGYC4wT9JjEXHHpDFDwBOSAM4DbpB0KiK2pD1hMys+x7w9pg16RKwF1gJIuhb4Yk3MiYgPTXwu6WFgq2NuZvU45u3T8nHoklZJWpXmZMys3Bzz9mrq4lwRsR3YXv287hugEfHXM52UmZWPY95+PlPUzNrOMe8MB93M2sox7xwH3czaxjHvLAfdzNrCMe88B93MUueYZ8NBN7NUOebZcdDNLDWOebYcdDNLhWOePQfdzGbMMc8HB93MZsQxzw8H3cxa5pjni4NuZi1xzPPHQTezpjnm+eSgm1lTHPP8ctDNLDHHPN8cdDNLxDHPPwfdzKblmBeDg25mU3LMiyNx0CX1SNolaWud126WtEfSbknDkv4k3WmaWRYc82Jp5p6iq4H9wLw6r30f+E5EhKQrgCeBi1OYn5llxDEvnkR76JIGgWXApnqvR8SxiIjqw7OBqDfOzIrBMS+mpEsu64F7gNONBki6RdLLwL8DdzcYs7K6JDM8MjLS7FzNrAMc8+KaNuiSlgNHImLnVOMi4t8i4mJgBfCVBmM2RsRQRAz19/uHxCxvHPNiS7KHfjVwk6TXgCeAJZIeazQ4Ip4FPizpvHSmaGad4JgX37RBj4i1ETEYEQuA24EfRMQdk8dI+gNJqn7+MWAO8Ms2zNfM2sAxL4dmjnI5g6RVABGxAfhz4E5JJ4Ex4C8mvUlqZjnmmJeHsuru0NBQDA8PZ/K9zazCMS8eSTsjYqjeaz5T1KxLOebl46CbdSHHvJwcdLMu45iXV8tvippZPmzZdYh12w5weHSM+X29rFm6kBWLB+qOdczLzUE3K7Atuw6xdvNexk6OA3BodIy1m/cCvCPqjnn5ecnFrMDWbTvw25hPGDs5zrptB854zjHvDg66WYEdHh2b9nnHvHs46GYFNr+vd8rnHfPu4qCbFdiapQvpnd1zxnO9s3tYs3ShY96F/KaoWYFNvPFZe5TLNRf1O+ZdyEE3K7gViwfOOKLFe+bdy0suZiXimHc3B92sJBxzc9DNSsAxN3DQzQrPMbcJDrpZgTnmNpmDblZQjrnVctDNCsgxt3oSB11Sj6RdkrbWee2zkvZUP56T9NF0p2lmExxza6SZE4tWA/uBeXVeexW4JiKOSroe2Ah8PIX5mdkkjrlNJdEeuqRBYBmwqd7rEfFcRBytPnweGExnemY2wTG36SRdclkP3AOcTjD2c8DT9V6QtFLSsKThkZGRhN/azBxzS2LaoEtaDhyJiJ0Jxn6SStDvrfd6RGyMiKGIGOrv9w+kWRKOuSWVZA39auAmSTcAc4F5kh6LiDsmD5J0BZUlmesj4pfpT9Ws+zjm1oxp99AjYm1EDEbEAuB24Ad1Yn4hsBn4q4j477bM1KzLOObWrJYvnytpFUBEbAC+BJwLfEMSwKmIGEplhmZdyDG3VigiMvnGQ0NDMTw8nMn3Nsszx9ymImlnox1mnylqliOOuc2Eg26WE465zZSDbpYDjrmlwUE3y5hjbmnxTaKtlLbsOsS6bQc4PDrG/L5e1ixdeMaNlPPCMbc0OehWOlt2HWLt5r2MnRwH4NDoGGs37wXIVdQdc0ubl1ysdNZtO/DbmE8YOznOum0HMprROznm1g4OupXO4dGxpp7vNMfc2sVBt9KZ39fb1POd5JhbOznoVjprli6kd3bPGc/1zu5hzdKFGc2owjG3dvObolY6E2985ukoF8fcOsFBt1JasXggN0e0OObWKV5yMWsjx9w6yUE3axPH3DrNQTdrA8fcsuCgm6XMMbesOOhmKXLMLUsOullKHHPLWuKgS+qRtEvS1jqvXSzpJ5JOSPpiulM0yz/H3PKgmePQVwP7gXl1XnsT+FtgRQpzMisUx9zyItEeuqRBYBmwqd7rEXEkIv4LOJni3MxyzzG3PEm65LIeuAc4PZNvJmmlpGFJwyMjIzP5UmaZc8wtb6YNuqTlwJGI2DnTbxYRGyNiKCKG+vv9w2/F5ZhbHiXZQ78auEnSa8ATwBJJj7V1VmY55phbXk37pmhErAXWAki6FvhiRNzR3mlZt/E9QM1mruWrLUpaBRARGyS9HximcgTMaUl/B1waEb9OZZZWar4HqFk6mgp6RGwHtlc/3zDp+f8FBtOcmHWPqe4BmpegO+ZWBD5T1DLne4CapcNBt8z5HqBm6XDQLXO+B6hZOnwLOsuc7wFqlg4H3XLB9wA1mzkvuZhN4phbkTnoZlWOuRWdg26GY27l4KBb13PMrSwcdOtqjrmViYNuXcsxt7Jx0K0rOeZWRg66dR3H3MrKQbeu4phbmTno1jUccys7n/pvuZfG3Ywcc+sGDrrlWhp3M3LMrVt4ycVybaq7GSXhmFs3SRx0ST2SdknaWuc1SfqapFck7ZH0sXSnad1qJnczcsyt2zSzh74a2N/gteuBj1Q/VgLfnOG8zIDW72bkmFs3ShR0SYPAMmBTgyE3A49GxfNAn6QLUpqjdbFW7mbkmFu3Svqm6HrgHuC9DV4fAH4+6fHr1ed+0fLMzGj+bkaOuXWzaYMuaTlwJCJ2Srq20bA6z0Wdr7WSypIMF154YfJZWldLejcjx9y6XZIll6uBmyS9BjwBLJH0WM2Y14EPTHo8CByu/UIRsTEihiJiqL/ff9ksPY65WYKgR8TaiBiMiAXA7cAPIuKOmmHfAe6sHu3yCeBXEeHlFusIx9ysouUTiyStAoiIDcB3gRuAV4C3gLtSmZ1ZVaOzRR1zs99pKugRsR3YXv18w6TnA/h8mhMzm9DobNHjJ07x+I6DjrlZlU/9L5k0rnuSN43OFr3/qX1IcszNqhz0Eknjuid51Ois0JPjwaN3X+mYm1X5Wi4lMtPrnuRVo7NCzz17jmNuNomDXiIzue5JntU7W3ROzyzuW35pRjMyyycvueRAWuve8/t6OVQn3tNd9yTvViwe4PiJU9z/1D5Ojgfnnj2H+5ZfWuhlJLN2cNAzlua695qlC8/4WjD9dU+mm1se3mA9evxtHt9xEEleMzebgpdcMpbmuveKxQM8eOvlDPT1ImCgr5cHb728pQhP/KI5NDpG8LtfNFt2HWr6a82EjzM3S8576BlLe9076XVPpjPVL5pO7aU75mbN8R56xlq93ne7Zf0Gq2Nu1jwHPWOtXO+7E7L8ReOYm7XGQc9YmuveacrqF41jbtY6r6HnQFrr3mlq9sYSaXDMzWbGQbeGOvmLxjE3mzkvuVjmHHOzdDjolinH3Cw9DrplxjE3S5eDbplwzM3S56BbxznmZu3hoFtHOeZm7TPtYYuS5gLPAu+qjv/XiPhyzZhzgG8BHwZ+A9wdES+lP12bTl6ukFiPY27WXkmOQz8BLImIY5JmAz+S9HREPD9pzD8AuyPiFkkXA18HPtWG+doUproUL3T2JKFajrlZ+00b9IgI4Fj14ezqR9QMuxR4sDr+ZUkLJJ0fEf+X5mRtag1vpvydfZw4dTqze4065madkWgNXVKPpN3AEeCZiNhRM+RF4Nbq2KuADwKDdb7OSknDkoZHRkZmNHF7p0ZXQhwdO5nZvUYdc7POSRT0iBiPiEVUIn2VpMtqhvwTcE41+l8AdgGn6nydjRExFBFD/f3+i522Zq+E2O5L4TrmZp3V1LVcImJU0nbgOuClSc//GrgLQJKAV6sf1ib13vxsdAu6ubNncfStk+/4Gu28FK5jbtZ50+6hS+qX1Ff9vBf4NPByzZg+SXOqD/8GeLYaeWuDRreHA+peivfLN/5hRy+F65ibZSPJHvoFwCOSeqj8AngyIrZKWgUQERuAS4BHJY0DPwU+164J29S3h/vx3y9p+EZnJ45ycczNspPkKJc9wOI6z2+Y9PlPgI+kOzVrpJXbw3XiUriOuVm2fKZoAeXxPqSOuVn2HPQCytt9SB1zs3zwHYsKKIvbwzXimJvlh4NeUHm4D6ljbpYvXnKxljjmZvnjPfQCycuVFB1zs3xy0Atiqisp+qqJZgZecimMqU4m6hTH3CzfHPSCaOVkojQ55mb556AXRJYnEznmZsXgoBdEVicTOeZmxeE3RQsii5OJHHOzYnHQCyStk4mSHP7omJsVj4PeZZIc/uiYmxWT19C7zHSHPzrmZsXloHeZqQ5/dMzNis1B7zKNDnM8f95cx9ys4Bz0LlPv8Me5Z82iZ5Ycc7OCS3KT6LmSXpD0oqR9kh6oM+Z9kp6aNOau9kzXZmrF4oEzbiT9/nlzOfc972Lk2AnH3KzgkhzlcgJYEhHHJM0GfiTp6Yh4ftKYzwM/jYgbJfUDByQ9HhFvt2PSNjMThz96zdysXJLcJDqAY9WHs6sfUTsMeK8kAe8B3gROpThPS5ljblY+idbQJfVI2g0cAZ6JiB01Qx4CLgEOA3uB1RFxus7XWSlpWNLwyMjIzGZuLXPMzcopUdAjYjwiFgGDwFWSLqsZshTYDcwHFgEPSZpX5+tsjIihiBjq73dEsuCYm5VXU0e5RMQosB24rualu4DNUfEK8CpwcRoTtPQ45mblluQol35JfdXPe4FPAy/XDDsIfKo65nxgIfA/qc7UZsQxNyu/JEe5XAA8IqmHyi+AJyNiq6RVABGxAfgK8LCkvYCAeyPijXZN2uprdNEtx9ysOyQ5ymUPsLjO8xsmfX4Y+LN0p2bNaHTRreMnTvH4joOOuVkX8NUWS6LRRbfuf2ofkhxzsy7goJdEo4tunRwPHr37SsfcrAv4Wi4l0eiiW+eePccxN+sSDnpJ1Lvo1pyeWdy3/NKMZmRmneagl8SKxQP847JLmN0joLJn/s+3XdHWe46aWb54Db0kjh5/m8d3HESS18zNupT30EvAx5mbGTjoheeYm9kEB73AHHMzm8xBLyjH3MxqOegF5JibWT0OesE45mbWiINeII65mU3FQS8Ix9zMpuOgF4BjbmZJOOg555ibWVIOeo455mbWDAc9pxxzM2tWkptEz5X0gqQXJe2T9ECdMWsk7a5+vCRpXNLvtWfK5eeYm1krkuyhnwCWRMRHgUXAdZI+MXlARKyLiEURsQhYC/wwIt5Me7LdwDE3s1YluUl0AMeqD2dXP2KK/+QvgX+Z+dS6j2NuZjORaA1dUo+k3cAR4JmI2NFg3LuB64BvN3h9paRhScMjIyMtTrmcHHMzm6lEQY+I8epyyiBwlaTLGgy9Efhxo+WWiNgYEUMRMdTf72BN+NVbJx1zM5uxpo5yiYhRYDuVvfB6bsfLLU3rndPDh/rPdszNbEamXUOX1A+cjIhRSb3Ap4Gv1hn3PuAa4I7UZ1lyc86axdc/87Gsp2FmBZfknqIXAI9I6qGyR/9kRGyVtAogIjZUx90C/EdEHG/PVM3MbCpJjnLZAyyu8/yGmscPAw+nNTEzM2uOzxQ1MysJB93MrCQcdDOzknDQzcxKwkE3MysJB93MrCRUufZWBt9YGgF+lsk3b815wBtZT6INvF3FUtbtgvJuW9rb9cGIqHtKeWZBLxpJwxExlPU80ubtKpaybheUd9s6uV1ecjEzKwkH3cysJBz05DZmPYE28XYVS1m3C8q7bR3bLq+hm5mVhPfQzcxKwkE3MysJB71K0lxJL0h6UdI+SQ80GHetpN3VMT/s9DxbkWTbJL1P0lOTxtyVxVxbUb3n7S5JW+u8Jklfk/SKpD2SCnMnkWm267PV7dkj6TlJH81ijq2YarsmjblS0rik2zo5t5mYbrs60Y4kN7joFieAJRFxTNJs4EeSno6I5ycGSOoDvgFcFxEHJf1+RnNt1rTbBnwe+GlE3Fi9S9UBSY9HxNuZzLg5q4H9wLw6r10PfKT68XHgm9X/LYKptutV4JqIOCrpeipvvJVhu6jeTOerwLZOTioFDberU+3wHnpVVByrPpxd/ah9x/gzwOaIOFj9b450cIotS7htAbxXkoD3AG8Cpzo3y9ZIGgSWAZsaDLkZeLT6/8HzQJ+kCzo2wRZNt10R8VxEHK0+fJ7KDdxzL8GfF8AXgG8Dhfj7BYm2qyPtcNAnqf6TaTeVH6RnImJHzZCLgHMkbZe0U9KdHZ9kixJs20PAJcBhYC+wOiJOd3aWLVkP3AM0musA8PNJj1+vPpd365l6uyb7HPB0W2eTnvVMsV2SBqjcznJDvddzbD1T/3l1pB0O+iQRMR4Ri6js7Vwl6bKaIWcBf0TlN/FS4D5JF3V2lq1JsG1Lgd3AfGAR8JCkuv8kzgtJy4EjEbFzqmF1nsv1sboJt2ti7CepBP3etk9shhJu13rg3ogY78ysZi7hdnWkHQ56HRExCmwHrqt56XXgexFxPCLeAJ4FCvNmFEy5bXdR+SdhRMQrVNZoL+7s7Jp2NXCTpNeAJ4Alkh6rGfM68IFJjwep/Cskz5JsF5KuoPJP/Jsj4pednWJLkmzXEPBEdcxtwDckrejkJFuQ9Oew/e2ICH9UTq7qB/qqn/cC/wksrxlzCfB9Kr9t3w28BFyW9dxT2rZvAvdXPz8fOAScl/Xcm9jGa4GtdZ5fRmU5QsAngBeynmtK23Uh8Arwx1nPMc3tqhnzMHBb1nNN6c+rI+3wUS6/cwHwSPUd9lnAkxGxVdIqgIjYEBH7JX0P2ENlrWxTRLyU3ZQTm3bbgK8AD0vaSyV+90ZlT6Jwarbru8ANVOL3FpV/iRRSzXZ9CTiXyh4swKko6JUKa7arNLJoh0/9NzMrCa+hm5mVhINuZlYSDrqZWUk46GZmJeGgm5mVhINuZlYSDrqZWUn8P0Y975eSzxxCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dionysus.plot.plot_diagram(dgms[2], show = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on built-in function uniform:\n",
      "\n",
      "uniform(...) method of numpy.random.mtrand.RandomState instance\n",
      "    uniform(low=0.0, high=1.0, size=None)\n",
      "    \n",
      "    Draw samples from a uniform distribution.\n",
      "    \n",
      "    Samples are uniformly distributed over the half-open interval\n",
      "    ``[low, high)`` (includes low, but excludes high).  In other words,\n",
      "    any value within the given interval is equally likely to be drawn\n",
      "    by `uniform`.\n",
      "    \n",
      "    .. note::\n",
      "        New code should use the ``uniform`` method of a ``default_rng()``\n",
      "        instance instead; please see the :ref:`random-quick-start`.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    low : float or array_like of floats, optional\n",
      "        Lower boundary of the output interval.  All values generated will be\n",
      "        greater than or equal to low.  The default value is 0.\n",
      "    high : float or array_like of floats\n",
      "        Upper boundary of the output interval.  All values generated will be\n",
      "        less than or equal to high.  The default value is 1.0.\n",
      "    size : int or tuple of ints, optional\n",
      "        Output shape.  If the given shape is, e.g., ``(m, n, k)``, then\n",
      "        ``m * n * k`` samples are drawn.  If size is ``None`` (default),\n",
      "        a single value is returned if ``low`` and ``high`` are both scalars.\n",
      "        Otherwise, ``np.broadcast(low, high).size`` samples are drawn.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    out : ndarray or scalar\n",
      "        Drawn samples from the parameterized uniform distribution.\n",
      "    \n",
      "    See Also\n",
      "    --------\n",
      "    randint : Discrete uniform distribution, yielding integers.\n",
      "    random_integers : Discrete uniform distribution over the closed\n",
      "                      interval ``[low, high]``.\n",
      "    random_sample : Floats uniformly distributed over ``[0, 1)``.\n",
      "    random : Alias for `random_sample`.\n",
      "    rand : Convenience function that accepts dimensions as input, e.g.,\n",
      "           ``rand(2,2)`` would generate a 2-by-2 array of floats,\n",
      "           uniformly distributed over ``[0, 1)``.\n",
      "    Generator.uniform: which should be used for new code.\n",
      "    \n",
      "    Notes\n",
      "    -----\n",
      "    The probability density function of the uniform distribution is\n",
      "    \n",
      "    .. math:: p(x) = \\frac{1}{b - a}\n",
      "    \n",
      "    anywhere within the interval ``[a, b)``, and zero elsewhere.\n",
      "    \n",
      "    When ``high`` == ``low``, values of ``low`` will be returned.\n",
      "    If ``high`` < ``low``, the results are officially undefined\n",
      "    and may eventually raise an error, i.e. do not rely on this\n",
      "    function to behave when passed arguments satisfying that\n",
      "    inequality condition. The ``high`` limit may be included in the\n",
      "    returned array of floats due to floating-point rounding in the\n",
      "    equation ``low + (high-low) * random_sample()``. For example:\n",
      "    \n",
      "    >>> x = np.float32(5*0.99999999)\n",
      "    >>> x\n",
      "    5.0\n",
      "    \n",
      "    \n",
      "    Examples\n",
      "    --------\n",
      "    Draw samples from the distribution:\n",
      "    \n",
      "    >>> s = np.random.uniform(-1,0,1000)\n",
      "    \n",
      "    All values are within the given interval:\n",
      "    \n",
      "    >>> np.all(s >= -1)\n",
      "    True\n",
      "    >>> np.all(s < 0)\n",
      "    True\n",
      "    \n",
      "    Display the histogram of the samples, along with the\n",
      "    probability density function:\n",
      "    \n",
      "    >>> import matplotlib.pyplot as plt\n",
      "    >>> count, bins, ignored = plt.hist(s, 15, density=True)\n",
      "    >>> plt.plot(bins, np.ones_like(bins), linewidth=2, color='r')\n",
      "    >>> plt.show()\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(np.random.uniform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "# summary method, full pipeline\n",
    "\n",
    "def score(t, g, dim, window_size, cloud_points, timeseries_average_k = 0, score_n=2, score_m=2, dmax = 0.3):\n",
    "    \"\"\"SW1Pers-score of time line (t,g).\n",
    "    t, g, dim, window_size, cloud_points: parameters for make_pointcloud\n",
    "    timeseries_average_k: k of moving average. with k=0, averaging does nothing.\n",
    "    score_n, score_m: n,m exponents given to pcscore\n",
    "    dmax: maximum edge length in rips complex. TODO lookup default value\n",
    "     \n",
    "    TODO: take parameter for expected periods, calculate other defaults from it \n",
    "    TODO: t,g as one parameter that can have shape (n,) (equidistant) or shape (n,2) (explicit times)\n",
    "    \"\"\"\n",
    "    g = moving_average(g, timeseries_average_k)\n",
    "    cloud = make_pointcloud(t,g,dim,window_sizecloud_points)\n",
    "    cloud = mean_center_normalize(cloud)\n",
    "    return pcscore(cloud, dmax=dmax, n=score_n, m=score_m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: create the test data from the paper to compare results"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}