AUC_sweep_with_Linfinity_balls.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from sklearn import metrics\n",
    "import scipy as sp\n",
    "from scipy.spatial import distance_matrix\n",
    "\n",
    "from diffusion_functions import *\n",
    "from semi_supervised_manifold_learning import *\n",
    "from animate_diffusion import animate_diffusion\n",
    "\n",
    "# commands for displaying animations in jupyter notebook\n",
    "from IPython.display import HTML\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Concentric rectangles in increasing dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def format_axes(ax, titlestring):\n",
    "    # figure formatting\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "    ax.spines[\"right\"].set_visible(False)\n",
    "    ax.spines[\"bottom\"].set_visible(False)\n",
    "    ax.spines[\"left\"].set_visible(False)\n",
    "    ax.axis('off')\n",
    "    if not titlestring==None:\n",
    "        ax.set_title(titlestring)\n",
    "    return\n",
    "\n",
    "def plot_label_comparison_colorful(ax, label_vector, data_matrix, titlestring=None):\n",
    "    threshold = 0 \n",
    "    label_estimates = label_vector#make_sweep_cut(label_vector, threshold)\n",
    "\n",
    "    im = ax.scatter(data_matrix[:,0],data_matrix[:,1],c=label_estimates)\n",
    "    plt.colorbar(im, ax=ax)\n",
    "    \n",
    "    format_axes(ax,titlestring)\n",
    "    return\n",
    "\n",
    "def plot_label_comparison_binary(ax, label_vector, data_matrix, titlestring=None, \n",
    "                                 objective_function = sweep_cut_classification_error):\n",
    "    \n",
    "    cut_val, threshold = find_min_sweepcut(label_vector,100,objective_function, orthogonality_constraint = 'auto')\n",
    "    label_estimates = make_sweep_cut(label_vector, threshold)\n",
    "    \n",
    "    classification_error = sweep_cut_classification_error(label_estimates)\n",
    "    orthogonality_error = np.abs(np.sum(label_estimates)/len(label_estimates))\n",
    "\n",
    "    im = ax.scatter(data_matrix[:,0],data_matrix[:,1],c=label_estimates)\n",
    "    #plt.colorbar(im, ax=ax)\n",
    "    \n",
    "    subtitle = f'\\n Threshold = {threshold:.3f}. Cut objective = {cut_val:.3f} \\n Class. error = {classification_error:.3f} \\n Orthog. error = {orthogonality_error:.3f}'\n",
    "    format_axes(ax,titlestring+subtitle)\n",
    "    return\n",
    "\n",
    "def plot_AUC_hist(AUC_vals, titlestring = None):\n",
    "    plt.rcParams.update({'font.size': 15})\n",
    "    # get bin parameters for consistent scaling of both datasets without displaying\n",
    "    _, first_bins, _ = plt.hist([[v[1] for v in AUC_vals], [v[0] for v in AUC_vals]])\n",
    "    plt.clf()\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize = (6, 5))\n",
    "    ax.hist([v[1] for v in AUC_vals], bins = first_bins, alpha=0.5, edgecolor = 'black', label = 'graph')\n",
    "    ax.hist([v[0] for v in AUC_vals], bins = first_bins, alpha=0.5, edgecolor = 'black', label='hypergraph')\n",
    "#     ax.hist([v[1] for v in AUC_vals], alpha=0.5, edgecolor = 'black', label = 'graph')\n",
    "#     ax.hist([v[0] for v in AUC_vals], alpha=0.5, edgecolor = 'black', label='hypergraph')\n",
    "    ax.legend()\n",
    "    \n",
    "    # figure formatting \n",
    "    ax.set_title(titlestring)\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "    ax.spines[\"right\"].set_visible(False)\n",
    "    ax.tick_params(axis='x', labelsize=15)\n",
    "    ax.tick_params(axis='y', labelsize=15)\n",
    "    plt.show()\n",
    "    return\n",
    "\n",
    "def compare_ROC_curves(hypergraph_x, graph_x):\n",
    "    n = len(hypergraph_x)\n",
    "    labels = np.hstack([np.full(shape=int(n/2),fill_value = -1),np.full(shape=int(n/2),fill_value = 1)])\n",
    "    \n",
    "    plt.plot(figsize=(15, 15))\n",
    "    for x, legend_val in [(graph_x, 'graph'),(hypergraph_x, 'hypergraph')]:\n",
    "        fpr, tpr, thresholds = metrics.roc_curve(labels, x, pos_label=1)\n",
    "        plt.plot(fpr, tpr, label = legend_val)\n",
    "    \n",
    "    graph_auc_score = metrics.roc_auc_score(labels, graph_x)\n",
    "    hypergraph_auc_score = metrics.roc_auc_score(labels, hypergraph_x)\n",
    "    \n",
    "    plt.legend()\n",
    "    plt.title('ROC curves \\n hyper AUC score = '+str(hypergraph_auc_score)+' \\n graph AUC score = '+str(graph_auc_score))\n",
    "    plt.show()\n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization in 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANsAAAFlCAYAAABiNJ3xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABn3UlEQVR4nO29e3Bc13kn+DtNEo2HJcqWxYbQIEBFImQ9QteMhUZmRiYAMwQgzzr2eGfKsZjEDvFIKpOaKStVm0S0uCnJnI2TCl3xxq4YIBkrNXTiyqyziWIbj9Bs0NmqsCmlYkmUI4p2BLIbAqhENinjRZF99o/Gd/Ddc899dd/bD7B/VSywu+/j3HPPd7739wkpJeqoo47oEav0AOqo41ZBndjqqKNMqBNbHXWUCXViq6OOMqFObHXUUSbUia2OOsqErZW46Xvf+165a9euSty6jjoixQsvvPAvUsq7TL9VhNh27dqF559/vhK3rqOOSCGEmHP6rS5G1lFHmVAntjrqKBPqxFZHHWVCndjqqKNMqBNbHXWUCXViq6OOMqFObHXUUSbUia2OOsqEOrFtAugJwPWE4OpEndhqHOl0GlNTU4rApJSYmppCOp2u7MDqsKFObDUMKSVWV1dx9uxZRXBTU1M4e/YsVldX6xyuylCR2Mg6woEQAoODgwCAs2fP4uzZswCAnp4eDA4OQghRyeHVoaHO2WocnOAIdUIzo9K6bZ3YahwkOnJwHa6OAqpBt60TWw2D62g9PT04fPgwenp6LDpcHdWj29Z1tiqHlNIiEvLPQgg0NjZadDQSKRsbG6tKlKQFTWPSP0cJmpdsNmvTbaWUmJ2dRV9fX+TjqBObA9wWebmQTqexurqqCIl25MbGRrU4+vr61NjoLz++GggunU7j4sWLaG9vV5vB5OQkcrkc7rvvvrIsdABIJpPI5XLqs5QSmUxGEV3Uc3XLEpsbMflZ5OUYH4k+QMHowUVGncNVw5idnmNlZQW5XE4tdFrkANDe3u5roYex+enHZzIZpFKpshmUap7YinkJbguzt7fX9yKPEm5m/YGBAZs4FtWYS13kQggMDQ0BKCxuGiM9i5+FXupGouu2fAzl5Pw1TWzFvAQvjgEAAwMDAKrDdxWPxy2fBwYGMD09jXg8DiGEev4o/G1hcUsiOOJmBD9jC8Lh6XjT58bGRqRSKZsxJJvNep4bFmqW2IK+BIKXI3h2dharq6sYGBiw7IC0uIsZZzEvMJ1OY2VlBXNz1voxR48exdLSEhKJBLZs2YL5+Xn1/PpCKoXQip1fOpf/ls/nbe4JoOCi8Bqj0/vSxT8ydDhtDlJKzM3NYXFxUUkH4+PjyOVymJqaQjwex9raWqRieM0Sm5/oCaeFTudyYqJr0QLTF/mFCxfQ29uLWMy/t6RYzkALnThBU1MTVlZWAABLS0toamrC4uIiUqkUhBCW5+fws5idUGx0iv7M+XweX/nKV3DlyhV1PulsnJD9EJyT+Efzms1mlV7INweap8XFRbS2tirpgD43NDRgbW0tctWhZGITQuwE8KcAEgAkgHEp5R+Wel2f93YkGi+9zOQIHhwcxMDAAObm5rCwsKB+a21txcLCAqanp30v3lI4AzdV53I5RWiElZUVtTtPTU1ZLGyJREItIr+L2QlO8+vG0fRnPnHiBK5cuYKWlhY89NBDasxtbW0QQtgkBtO8mBz3TvOaTCaNmwM/75lnnlG/DQwMIBaLKakgStUhDM52A8BvSCn/QQhxG4AXhBAzUspXQrg2AGdRTEqJ48ePW44lpyW3fum7HH85/Dc6tqury0Jso6OjmJ6etvmudJO7yQQPFPcChRAYHh7G008/bfydjCT6dYjjDQ4OGsccBE7RKbRA+XFcYpBS2rjtQw89pMZEVkAS3bj4yzdEuvb09LR6X0SsxBlNRMV1Qz7XbhtH0I2lGJRMbFLKNwC8sf7/t4UQ3weQBBAKsTlxKHpRRFBtbW3I5/Nqsrq7u3Hjxg1HvczJEQwAa2trljFMT0/bFhiNi8bBDRdra2tqwei6H7ckusG00DnGx8fR2dmp/ES6iMX9bcWAW/BaW1vVhkMidldXF/r7+23vh+7PQeMjIiCioetls1nlAiBCPH78OK5du4YHH3wQ8XgcPT09yOfzOHHiBO69916kUikbUQFwlFicfuOczencauJsCkKIXQD+DQC7AlEEvEQx/hLOnTunzmtra0MsFsPi4qLlejRxfcwRvD5uywvhuygtCACWF0PjIhGTRE/6nEwm8dprr+HmzZuWMRw9ehQf+MAH0N/f7/rcNI5kMokbN26oZ0mlUkrR37p1q6uFrZRFQiIeF6G5iN3Z2Yl8Pu/4fvTnITz84osYPHEC4vJlDO7cidZPfAJ/xaSQVCoFAOozvQv+jpPJpG28JNHQ5sPHQu/aJM3Mzc1hZGQEMzMz6n12dXXZdLgwCE6EFRcmhHgXgFkAR6SU3zD8PgZgDAA6Ojo+oBsgnMAXHkE3ggAwilu0UEznOcEPxyJiNS0sug9fHIlEQnEh+jw2NuZqbOEcg/5PJuy9e/eqMZEhhZ5tcnLS8rnURZLP5y0bDj2fLmno85BKpTA0NKTGAxQI7SPPPYeGd95Rx8nmZnxjaAgv79ljm8ObN2/aytRzo5BOOMlkEslkEkNDQzY9HYDNcDMxMYGFhQW1cV+4cEF9pvdehJvjBSnlI6bfQuFsQohtAP4fACdNhAYAUspxAOMA8Mgjj/imcD+ytJO4RRNn0sucFmFfX59tgXEOx40b+rj4+ABgfn4euVwOi4uLijMlEgncf//9nlZNzn37+/ttsYT0DOl02rL4aaGFFRsZi8UC6zpEaHy8yWQSHz9xAoIRGgCI5WXsO3XKRmw05zqGhoZc1YDe3l7b2Lh7gP4fi8UsojFBv25V6WyiMJrjAL4vpTxa+pCscJOlAaidk+94hEQioXQkOt6P9YsWGODuVnAichqfycDhxdE4+Lj0MdJnJ5E4rEXipcuYfudjsQRKj40Z77H92jVlRSWMj4/b1AB+b/1+9J3TPJl+c9tIoogUCoOz/QcAvwjgJSHEP65/96SU8lulXlgPs9E5VDweRy6XQ3d3tyK01tZWxGIxCCGQy+Us5npu/fLye7lxVN14wHU1bnI3iehTU1OK+4QFtwVWCrzmX+f4JgnCshl0dAAG9WH1rrtshMU/m3xzuhMagOO7pGfRN1k6h4OMJsWIkF4Iwxr5dwAiiWGy7YwahyJrZHt7u0WZ12VuIpAgjkuvHZ3GFY/H0dnZadHtOjs7LQ5WznVJfzERXNThQkHhNv+NjY2IxWJobGzEx5aWsOdXfkUZPe7+xCfwIybGqmc4cgQYGwOWl9U9rm/bhm89+igA2Ljbu971LuUyIORyOfXu/b5Lk0Wbsg5yuZxtoyQjEOneJjdHUfMZloEkCB555BEZpD+bm5/Ny3jCz/NS9vn1nXZ00/VNf2dnZ3Hx4kWlsANQUQ6UVsLHV61R+4DHJnDyJOTYGAQjINncDDE+Dhw4YL/YyZOQTz4JXL6Mq7ffjlP79tn0NQJ3OtN9Afh693ysTu+SrJpkvebXSiaTuPfeewOHcEVuIIkaTmKSl/HEtLv5iXn02tH165v+9vX1KccsN2rQZ05cwEaYWDabxfDwsG2n1p+nnFzPVUw9dMhCaEDB6IFDh8zEduAAxIEDyOfz+Pq6NZCju7sbjz32mNGgpetffpzQ/N3p+jd3kJPLgdDW1hZ6CFdNEJsT/DoiuV/Mb8yjX8OD267vtEh1/+HAwIAiqFwup4wqenB0NXI9XLoU7HtsRIXohAZAvQfT5qZfw68T2o0wefqPPo6wsz9qlti8lHedwxUT8+j1uVjRz2m3dYqKqIb8Okc4GD3Q0eF4ChmrdD9oa2urRdJwWtRB3j0dPzk5abnG5OSkzT1B4BE5uiRUiqW3Zgv+OIl6PT09xt0wFouhq6vL8t3o6Kjj8V7g3KmYIjK0AbiBdm56rrNnz+Lpp5+26Y4VxZEjQHOz9bvm5sL3DiBjFRkhqFDRwsKCJVbS6dmCvHuKnyVR8fDhw2pTO378OKSUlpw2OodiN3VfXymFlGqWswHBRD3AX8yjX7jpAn6IgEQpDuJqus9wcHAw8iDZokF62aFDBdGxo6NAaCZ9bR1+dWI3+H33+n35XwAqa4KiesiAtmPHDmzbtg25XC5wYIQTaprYAP+iHg9yFUKoSktA+VJQCLoYNDAwoNI+9GtQqFHUQbIl4cABV+IyoRhi0eHHvyiEsBic9E0xnU47xn8mk8lQK5fVPLG5gYt6yWTSFleXTCYjSUHxs8PSzm4KSyId0hQcHcYOWy3wQyxh3cdpU+zv70dvb6/iaKaUnVI2BI5NTWy6qEdOZtNEBkVQJV1HX581BtOUZaA70MPaYW81eG2KXvGfHKXM96YmNiC6pMAwiIAiMLyuEYbIdavCz6YIlEdM3/TEFmVSYBhE4Pca5RK5Nhu8NkWgfGL6pia2UkU9PwiDCOqEFD10cz13L5RLTN/UxFbXd+ogI1kmk1Hvn0K0KCigXGJ6TQQimxAkQr7aounriBbFBqCHAeESiFyTESRBe23VxbRbB6a1QalPHJUwMNUcsZUaJlXH5oXb2rhw4YLl2FLCropFzelspYZJ1bF54bQ2TEnFlQgKqDnOBlgnlVAntDoA89rYvXu3hdCoBGK5jWQ1SWxOvrO6CFkHRflzvPbaa7hx4wZOnDih1I39+/crK2S51k3NiZHl8J1FhbpVNFrQ2qCwvO7ubly+fBkLCwsqb66lpQVSSszMzKhaJrxURZSoOWITQigxgAiL8sKq2XdWzTVGCLW+GXC/KmV56FhaWlLFc3nyajKZjPx5a06MTKfTqlIxLVqSw6tl0eoI04J68iSwaxcQixX+njwZzhiDulOqFX19fRgcHFRZ2G4gQqOislFvLDXF2fTaHVVXIsABYVlQT560VoKbmyt8BgKnk1lQq/PqBr0MQnd3t6UfBEc5CA2owQgSvyXMqhFSSkuF5MOHDwca865d5nIfnZ3A66+XPrZanVcO/hzJZBLt7e22xis6qMVWGPUhN1UESa2a/cOwoBZRyMo3anVedXC9bXh42FL2oqWlxXJsa2srgEI5ivHxceTzeQDRidA1R2y1aPbXLahU4IbrcH7gVLDKpZBV4DFyVPu8OoH0NkoMJeL7wAc+gGQyie7ubvT09KCrqwupVAotLS1YXFzE9PR0pBFJNaez1aLZP6zsA0P1bq9CVr5Qq/PqBj5eHtXPjWhkYBscHHQsi3DLptjUcspMGGkcRRSy8oVanle/4POufx9VNr9tDLVmIAFq3x9UrbhV5zVM49CmMpAAmz9lximr2OnzyZMyFN/bZp9XE8LSp/2gpsTIWwF6ww1go+9Yb2+vre7/yZMSw8M3sbZWeJVh+d50mLgeULlmH2GhnCJ0TYqR1Qg/IpjXMbqPKJlMQohCnUtqVE+tg2lxtLau4sqVJtt4OjokXn89HO5kCjWjYN/h4WFb+Fku1xe6XmlCmNn6YYnQbmJknbOFAD9xj36OoV2VGinyZopzc3NYXFxEKpVCe3u7spxduXLYOKbLl4Hjx4+XHGBrii6hRoLARqUy2iSWlj6GL31JYnm5sFBNnDaMhR0k1lQ/lsok0LGm+0fBoevEViL8hDoBwTrRJJNJtZiBjR4Ara2tlrAvANi+/SquXr3DNq7bb7+KXC5XVICtbjU1tU6ifma6ufxXfmWPIjQCb9d2+vTpwA0GTePzO5/6sfF4HBcuXLB0Fo2ipa8JdWIrEX7jHsOIjRwdHQUATExMqO/27TuFb37zo0pnA4Bt265j375TRQXYmjiGqYZHY2MjLl68aPkun8+7Rrnk83m10IHiYzCDxJrqxxJaW1ttPcGj1jmrnthqwRztx0/j5xjShXK5nK1XG1AIrr106RIWFxfR2tqK0dHR9T4Bf4W/+7sP48qVRmzffhX79p3Cnj0vY2jIHnvpNp9SSly8eFFxVU4MeqjTP/3TP+HKlSuW786dO4ft2/+9kdPu3LnRALG5udlWtsDUAdYNQXxjpmMXFhZUM5NyxYBWNbEVI2tXAk6hTvwF+jmG46WXXrJ9R8G0iUQCo6OjqkY9MIXHHjtpET0B4Nvf/jYee+wxdX0/IhyJsJwYgEIeWE9PD+LxOF599VVLo3mOj3/8eXz96/ssouS2bdfxyCPP4ezZl20NEIHCwu/s7Az0Pv3OOX9OJ5QrQqZqia2aZG2vcfqtJe8VDiWEwMGDBzE+Pm5ZzDt27LBwkbGxMVuEOjemAAU9j4iTcrv8iHBDQ0PIZrOYn59X17/77rshpVSZzYuLi2hpacH27dstxwHAL/7iFuzbJ/Abv7GKK1fi2LkTeOSR57Bnz8sAgPvuuw8AbB1HKT8xrDknFwkXFROJBLq6uvDaa69Z7u+06YUtVVUtsVWTrO01Tj9+Gr++nFgshvvvv99CbLt27bIQG29LLIRAU1OTchU0NjZiZWVFEej8/DympqaQzWaxsLCA1tZWm57DF7qUEjdv3rQ8o5RScR9uLFlaWrLNx9raGj75yTzuvDONs2fP2jjZP/zDP2BZa3jv1mq5mDkHYNuoKdgYABYXF5FIJLB161b1LICVw0WRWV/1fjY9B4yjXLJ2GD60IMfooUPARs6V3mKKOJxcL1zDd/GOjg5LHlcymcTBgwctjRf37t1ryXyfnJw0lhMgHTEWiyGfz1uuwbkpUHgv+/fvx7Fjx4w9zJubmy0E19raiq6uLvT399vu6wYv/dM0j7Qp0Xhp7E1NTYqI3Din15qrWT9bNcjafnc4+o2bzPXF4OXLMb1kMphks1lMT09j//79AKD6PcfjcfT39ysuxztn6npVMpm0zSeJlnNzcxgZGcH58+cBFBbfysqKOi4Wi1msk/pzEMfL5XJobGzEli1b0NXVZSG2kZERRYD6Ig6qs3nNp8koQvoo79MnpTS+yyhqk1ZtbKS+8J566imV7Edwi13ziif0Owa/tUP81vBwG5dJPBoeHkZ3d7eKJDl27Bj279+PtbU1Ven39OnT6hqxWEy5CHTMzc2phhK8aXxzczMWFhbwuc99DktLSzZCo7GZ3gnFEZIYODw8rIxXeg/zmZkZdHV1+Wo8Xyry+bxtY7lx44bls5SFokBra2u298IJjlDq5l61nI0vPNLRSOfo6upSiw2wT0JY8rbfHc6vk9XPuHQrqxACjz32GKSUmJiYUEQBbIhlnCs4SQPkSuDGCHq2eDyOM2fOqGOJ0MhPR6LliRMncO+99zrqSlykdRPDuI5I1wiT0E6fPm0xpg0MDNiMTsCGGOlkHAm7r1/VEhtgXXic8GKxmNqJ9B0xSHSBHwTxoQHORBlkXCbxSAiB0dFRi67ExTFOaJlMBi0tLRYDxtzcHFKplIUohBDYv38/ZmZmfM8HbQpuxOJlwOCWVNPzlgLiqLQxDwwMYGpqShFaW1ubxYJqIjJgw3pMG04YybRVTWzAxmSYdnzTQ9P3UkpbeFExk+R3h/MiylL1AJOuBMBmMo/H44rj0TMTR9y6daslfCufzysdSrcaNjU1IZPJKKOHWzSK6Tu/7ytszM7OQkqpxFvanG677TY88MADAGAhtmw2q/7PJY/GxkakUikIITA7OxtKJkDV6mwmeBkYCLOzs8ZjZ2dnAfjX33RxyC3XyYkodZ2sGD2Aj0PXW6lQzenTpzE1NYXe3l5VWwMAzpw5g9HRUfT09AAoiE40Lh7RQVySdGNdZyNCI6unH/h9X2GBpIdMJmMb4/ve9z4AsOisqVQKuVwOU1NTyOfzSvKYnJxEb2+v0pNpLgYHBwOpITqqnrMFhS6uEYhggjjD/frQ/DhZnfQpP3qAEMLCsUicPnr0KBYXF3Hs2DHlC5ubm8Pu3btVlkBPT49l3DQuzlkbGhpw/fp1NQ5uNeTjJN25XOW6g4JLNbr7IhaL2Spp0wYSj8dVNE42m7VxdNqoS33eTUdsgDPnunnzZmBnuB9xyA9R+iVIJ5APqrOzU137oYceQiaTwbVr17C0tKQ4FBGJKTLDSdTlBpaZmRmbaHn27FnlUgAKZvR8Ph9KrcWwoY+BpBGTcYaiYmjN6BkXxN3CCJ7YlMTW1NRkU4QBgBzpQf0lfsQhL6L0yyXdoBsn6DpEPHpkBjmh3SyVxFn5s+lWYLo+L9fd2NhoifoII8IiLHA9DIDS4XTjDBlT6PnIR8nhx5HtF5uK2EjB3b9/P1599VXH44LE4QWBF1GWYjTQyyU4pb5wkBOcd2yhLPD29nZlROIW376+Pss4dU4IwBLJQp+roVw5ETxlTXArYiqVQm9vr+V4J6NVIpGwuAnCirLaNMTGdTUSd/SwIEKQOLywUYzRgD+blNKyiBKJhPEcSmN56aWXsLy8rCL2SQchwiAjAZVaIEJx89mR/gZUV/fXYqQH04bS2dlpITY9m6JYbBpi03cpwC5WffCDH8T169dL9peUG/RsuvJOOzD51LZu3YobN25g69at6tmXl5eV7kZxjfTMfIGZHPVcp0wkEujs7FT3n5ubw/33328ZZzXMZ1DpwbShkMVycHBQOfRLdWgDNWb694LJtM5x/fp1DAwMVKTFazHQxRcKoCVQ1Podd9yBtrY2vP/970dra6stLGlkZETpbtPT0xa/EYHrXsAGl0gmkxaxipfrroam8DqChOnpRqu9e/eqJFnu4qBsiltaZ9N1AzLrc7S2tmJkZAQzMzMWRViPYqg2mEK7Xtda1ZBuQjpTJpPBI488YkvOPHbsmMrqJv1F380nJyeVLkZl8/r6+tDb22vJJiBwN0RQy2pUSKfTuHjxItrb29VmQsWJTK4KkzFoaWkJiURCGYAopEvX94pBzRKbKYubIiVIP6H4uJmZGUt30monNCmt5Qkotu/KlStGPZR8RlJKW2skcgdQJAXX2YhgWlpalGhIFby4ZdGk13R1dSk3RDGW1bAhpcTKyoolqp/728gg5Ga04moIcfIw9dCaJDZTnCEPVCbO1dvbq3ZqclpWm+hoyskCnMsT3HbbbRZi48q7vomQf4ljaGgI6XQara2t2L17Nzo6OtSCJELTfUsAbJyQ58AB5QvHcgJtOEBB53LTR03n0l+vONhSUJPE5hZnyEVEncCqjdDcsgD4wiGQ7qTnupF18PLly5br5/N5mxmbxCrKFhgcHLTcgyeAkjTg1xlf6fklgtOjR/wSTBSR/hw1SWyA/11I5xpun8uJYrITtm7dqhY8sMH9stkspqamlBOfwtJIpOTcii9Ek45LIB2msbFRBeVykZGSLitNYBxSSlt7X8AfwZQa4eMHoRCbEOIEgP8NwBUp5cNhXNMLQXehKGpKlAJuzte5cz6fx4kTJ2z+HSq6Sudz7sezj4kYifjoOL4x6T2mdV2Q8r9I/y3GP1hOEKFxzkw6mx+CCSPCxwthcbavAvgjAH8a0vVcEXQXCjvHLUzosXhk5CATtJ7Aqcft6WITf/bh4WH1vb4x6frd8vKyReSk2pQUgVKNc8chxEbhI26NBArhW3rgeDEhd6UiFGKTUp4RQuwK41p+EHQXctPxKm000e9NpubGxkasrq6qyHTiTk1NTTbHMwfn7pyD80Bc3YwPbLhIKAscKFgct2zZEvrcFSPO+zmHXBXAxrw2NjZi586d6vtKSjRl09mEEGMAxgCgI4Qm0EF3oagtTUGhE4G++GkhEMEQwZkiPLxqUdLGFI/HLY5nqjRF5vJjx45ZxkDR8GHOXTHifJBzdILUy2e4ceWoVY2yOZyklONSykeklI/cddddoVwziB7hJ7mznOBEoI/h/Pnztkj9dDpteT4n7m6Kjunr68PAwADW1taQyWTwyiuvoLu7G8lkEplMBu+88w5aWlqUj1JPkjUVzylm7rg471VAqZRz+BzRnJw9exZPP/20YxR/Kffxi9DqRq6LkX/jx0BS7v5sblygUqIkzTvpYmTtI8MEr+XvNs4gIhl3/BPIsd3W1oadO3fadnWeLeA1d37GosdcAt4iaTHn6Ofz2qOHD9t7IIRxH8C9bmR1h1KEhCBcoBygsndAQQcj0/yZM2cwNjamCOCZZ57x3BC8uDnfTIUQtjJ3xM1GRkZsPrPBwUH09/f7mju/pfy4/kzwWszFnMPnwC9XLuU+fhCW6f/PAPQBeK8QIgvg/5RSHg/j2mEhakuTX3hZRokgeBWtYnLeqNa9lBJNTU3o7e21NDHkcMrvIw7H547rcJyjBemXFtRxXKyz2cs4pMfI1oRTW0r5yTCuEzWqwVfkZRkFYHM0+11YPOeNnNhAwX3AfVB6Ja2JiQmV1c0RpBq0H4tlMY7jUpzNunFoenpaRcXQZ94RqSac2nUEg5N1D/AfGuV0TcDqvAY2QrDa2tqQz+ctRYNIh6OEWrpWUN+kH4tlMY7jUp3NfX19KlKGczRTLZqondpV31hjM8JNEZ+dncXKyoqldFwQ87NuDOA4fPgw0um0pUcbrzYGwDGTQh9nKYYPnjNn+uz0XEF9c8WMr9T7uBlI6pytzNDFFWouSItg//79KsGTiKsUYwAHL+5D19toqAgbF6VMCo4wRD3SJ/n3XhtKqSqAX19hlKrGLWGNrBR0qUEXVwYGBvCDH/xA9Qt74403MDExgUwmg4sXLyr/FhWX9boXL25DRVkBqM/kQ9JBzm+TT8pPM5Mg1t4g/izT/BWLavCz1jlbRPAyLtBLbm9vRy6Xc2zvRPUwvMQZfcHPzs6q88gaCZj1D7q2affn2diTk5OO3MqvtdevMSXMaI5yGD/8YNMSW6myd6n39mtc8DJqUG1+PxV5+YLnBM3vpc+BXiJPT1FJJBLKgieEcK3HoTu3nZzdXiJdUOOMF8ph/PCDTUlslU6n8bt7E/L5vOO1fvSjH9nKzHnd2/R/02d9UUspbZW7FhcXMT09rX7r6elxrcfhZ+69/FlB588P/HLeKLHpdLZyxLj5gZ9oBCkLOVh63RCO5eVldHd3ezp93T57jZH0NJ4LNjY2pppz0G9+wqq85l4X6ZyalUQRzVFpP2vNcDa/YmEUu2Kx4/WKRpBSqoiO7u5uXLp0yaa7AfbcM45SubibSOeWK+d2LcB97oM0K+EIM5qjEqgJzuY37o4Qxa4YBH5273Q6jenpaUsnz5/85CfG6+m16/l9gnBx02enRV1spL+fue/r6zPGYepiphf3KyfCsIxWPWcrRlmu9K6oK+Q0bqCwe/Nn6unpwf79+zExMWHpFErYsWOHpaiPzuX8cnETB6RYSdIJ+dxSCfcgnTdpAfqZezeRrtwGDS+pKSwbQNUTW1CxsFrMvKSQcwcudx5TrzD+TEAhpKq9vR1AwSK5bds2pFIpS0iV/sK9nLVOG1Ymk0EymVTpPdS1M5VKWcra+em8mU6nsbKyAiGsZfCoHJ9pXH7mL2qDhh8XTViW0aonNiBYlnW1mHkJbi9qYGDA8kzERTiampqwd+9ex64xpurGPFIEcN+w4vG4Ejmpa2cikcD9999v4b6AcyVp6trJiZesl6lUCu3t7UXNfdQGDb+EFJYNoCaIzVRyzU0srAYzL78vYK5vaSojJ6VUC5oIz+k6vEQd+ePy+bz6PZ/PY35+XpXe1jcsGkMmk4EQAg0NDWhublYmf6CQIUBtqUwdW4kzNDQ0oLW11VKRuLW11VJVuVJwEhP9ElKQzd4NVW8gOX36NCYmJtRuQz2fvZTlSpt5+X11g4EedU4GgEwmg4mJCeV344YgIYRyLhOGhoYUF6fmj1QV+vLlyzh37hxyuRxWV1eNBg9KOSFx9rvf/a6ttPnCwgLW1tZU9Ag3vHDO8Nprr9niKOnciuDkSWDXLshYDKutrZAnTwKwG9f8umjCCPWqas4mZaFgCy8rzsuMm+oZul2rEhElphdFTQz5DkpdL0k3470Kenp6cPPmTVtBHpO4mMlkbGUPKLj57NmzSkSkUgdSStuurcNrx+fHcJjaDJcFJ08CY2PA8jIEgKYrV3BjeBhbAEzdeadFTATcDTph2gCqmtj0l0nZy7zMuB8i0kOSgI0o897e3sgWg9uLomege8diMUvNEQIRCzWU12uTAM5+MQCqsUg8Hletnnbt2oWBgQFIKZHNZnHixAlfz+MU66jrnvzeFWk8eegQoHHorWtr+PF//a84+5nPKEKjdCMvQgrLBlDVxAa4y8t+Q4NI1Mlms0gmk0rPoQzmpqYmm54X1tjdXpRubKB0Fz0QmGo5EqHxtBh64WTK10EiNwd/Tl6yvKGhAS+88IKxWytg1pOllI4lzINKH6Hh0iXj19uvXgWwEZZGBiIvQgrLBlD1xOYkLw8MDPiyJOnckZT3ZDKpWiT19PRYkiiD+E68EORFmZ6Vg5cu4NchQstkMnjxxYfxd3/3Ybz5ZiNuv/0q9u07hT17NsTKVCpl6V5KFYRJRKdOpbFYzNY7QN/xOeemUgv6387OzrIGgQMAOjoArRE9AFzdvh0ALNXM6Dm83k8YNoCqJjYvedlvX2cTd+QWM6c0+bDgdS2uO3ARU8+S1kUy/repqQmXL38Q3/zmXqytFV7r1at34LnnPgIA2LPnZeO9qUQ558INDQ344Q9/iO7ubsRiMVXOTo/4F0Iowtq9ezc6OzuVyNvQ0FCR6mVSSogjR5TORrgRj+PUvn3Gc8plTKtqYvMjhun6gkkhd+MYvFFgOeIn3URfvQsmLeSuri5bZV99jH19ffjUpyTW1qzfv/NOA06d2ofR0RZLVD+BG1n6+jbqdeRyOQvHo02AR/xLKZWY2tnZaRk316vLBTW3jz8OAUA++SRw6RLefve7kfnYx/ByZ6c69vz582hsbER/f3/ZxlfVxAa4i2GnT5+29XWemJhAV1eXmkRe7IW4lr7gCFESGnEvXgGLh0Hx0C6+yXBDEOCulF++bP7+2rXtlg2ppaUFTzzxhNHIwvVBvxKDn2Ojhs1B/fjjyvLY1NSElZUVpfNSIdwLFy6gt7e3bBtC1RMbYGbz+XxemcZpErnY1dvbizNnzuDixYsqBIkIrbm5GXfccQeuXr1qiUeMKn6Sc7PBwUE1Dp7SYgrWpU2GRzLoOgaHg6qCnTutn5944gmjkYUQNGInDIdvqXAj/IaGBuUHJCmGJIZyct6qd2o7IRaLoaury1I9mItdQgisrq4il8thfn7eks6yvLyMGzduYGlpCa2trdi7d29kUeV6ZD5g3zzcFjJlPOTzeYvYefr0advxR44Azc3W75qbJT7xiRct31EyKC1Q3SAUxIkblsNXv6bbZyc4Oag/9KEP2SpBj46OllWEBGqEszmhv78fvb29lurB3GLHDSg8QZNaMt1zzz2BxLRi4LTjcpiqXgFWQp2bm7P416gQD18wBw4U/h46VLB+79wpsX//abS0fNdmYMpmsxgeHnbUb718TyajThhB36VE2LtZrnX3RCX8fzVNbKaYSZpEirY3OVx5TcZSfSd+YBK1qJGFlFLpcKSr0aLSI0u4CORkVj9wYIPopASOH/8hcrkN7uAnH003SsXjcYvPzGTUCSPou5QIe7dNguavklkgQA0XaSXjiKm6Lw/Q1UttA4Ws6Mcee6xsk8wXAoHGxVNRcrmcivynsVGCJ+feAGwGFX4v/TMvPQ7Adg+nMXuFK/kdQxCY5sqvwcWJK5J6UQy3DAqx2Yq08phJqv5EJudEIqHCsGhHI9ACP3funApviprg3BarHt1B0S38XLdmGHqxU9MiMoVx+XluneMD/i2OpcxpKQYXN8t1NWSB1KSBhCYrmUxicXHR0lqps7NTmbC7uros542OjqKnpwfJZNLSLjfqsToVMNXH197ebjHS8D7aiUTCcuzExASWl5ctx5tKIoRhwHAyPEQxf6WO18lBXS7HtRtqkrMRTM3fKeYtn8/b0jsopUQIUdbJ1ndcAEalXUppy95OJpO4ceOGCivj4vKWLVuQSCQsxycSCYtuRQRbir7iRABhE1yYEfbViJomNn3iKeZNj3owvbRKjdWtfACNn2N4eBjpdBq7du1Si42sktls1laNa3FxURlOgEJYWjKZ9DRgnDy5YcXs6Ci4EQ4csAcF6D3OSiUA630FfuEXOvDoo/Acby2iJolN3wG5fC9EoWZGNpu1KPBSStdqvuUaNy8fwCNaUqmUjVMD7s0wpJS2+En9XCrmQzDpKyz9C0DBMT42Bnz/+6/g0Ucvqch4IjT6XOpcmu77hS88iAceeKDi+lUUqEliE2IjAFavJvzKK6/ggQceUDs6AFXcJpVKqdi+SiST6sYGIiy+kE2VrgAzNzZ1mQHgGZ2iP6ch/QvLy8DExD3Ytu0vbBxNz8UjBJ1Dp/seOiSU+8I03lpFTRIbD4DVTf1vv/02AHsqCV90lSxP7mZtC+Kz4v4vTnBUNly/thsc0r/w5puNNh2SN+4odQ6d7nvpkgSwOQiMo+atkUBhJ6cdN5VKoampyValihYFVYLysuJFBTdrW1+fe/FS/Trk/kilUjh8+LAqQWe6ths6Oszf79xpd4LTBuA0hysrK5Zz3O7tdN+77lp1LMBby6hJzgYUFuLw8LClyyYXtdzCdkw1G8sRqV6Mtc1pPMQJyTfHRbiWlhZ84AMf8EzLIRTSvySWlzd+L8RUfs8xJcfkezP5Cd243ZEjEsPDN1X+HQDE4zfw6KPfwupqS/mTTiNGzRKbVxycnj2sh+3oYVzlUMKDiope6O3ttYQ3cbFybW1NxYZ6XbugHwk8+aTE5ctQGd4tLYWEU6eNQReHyU8IOPettt5XANiCJ55YwZUrjdi+vXDf0dGWkt9HJXRyL9RkuJZXCNEbb7yB69evY2RkBDMzM7bQH918Td+Xy+oV5kLwCm8Kem0prT259WtRzZbe3l7bfSmzO+i86vc8fPhwSe+hwjq5Y7hWzepsTlEZ8Xgcd999NxYWFvAnf/InNmvlzZs3ceLECUWYlWjaEGY0gxD2epLcUhiU0HRpgc/H7OwshBDYu3evIrREIqF053Pnztnm26uUXTERI24pOHpKU7l1cjfUrBjpFQcHWM3rhOeffx5tbW1VU568VHhlq/vlom7SAhG0k8hKIXCZTMbWa84tlaUYHdaLaxUTy1ku1CyxAe5xcE61DAFg586d2L9/v+W8ctbLCEuM9MpWl1Li+vXrvsQpL31SL5dAoIJJACzGlKeeesoz0iSoDus3BcfNvVJJ1DSxOUFK51qG1BFmZmYGUm40d4+ijJ0JYegTtKgo2FpKact32717N65fv+66MAHrhuUmLfDPfBEvLCxgamrKtpApDhVwlxi87snhl2uVK5YzKDYdsemiST6ft4g2//zP/4x77rlHvSjqAuMnQbFUcB8f4M9ip0Mn1r6+Pqyurlr8a5StzsVpfWECMBK5mz5pWsSJRMLC0XQjiR+JIYgO68W1qjmYuWaJzUkU46KJlBLnzp1Twb3nz5/Hm2++iTfffFOd5xbWFCb0bi+cAPxWDnYSo3RfGNeT9IVJBBh0g/FaxITHHntM/V+v+hyG+OzFtcJ2r4SJmjT9+y07zkOK6Ds945lQqrnZDXyhmjLHAf/EbjL1A7B1CO3pKXQ0/e3ffhnj4524enW78mNRwdagG4w+7/l83hYI7ZS9HZb47JU1zjkc/z8fT5TSi5vpv+Y4m18lGbCXg3PS44BoZXpd19BB2eZ+7m3iVrzMAd3njTfewOjoafzP/7kX77zTAMBeITno85rm0622ByHIO/N6dr9ci/6vN1Up1ucWBleuOWJzUpJ5eyJ9QvmOyIOWAdg+R0lwThbSxcVFjI2lMTPTj8uXhSWfTIdJjOLjnZ2dRT6fx913341vfOMRRWiEd95pwF/+5cfwjW98HF/+8iqOHm3E44/798dxQ0aQhR+WOT6IQSUsIr9lemqb4GQVI13FNKF8YczOzqrfyBoJRCvTu3HWF198GH/zNx/E9euFe1M+GWAlOC+9iXxhpKdeu7bdYSxbAABXrjRhePgGvve98/jwh38U2BJbjCUxDHO8X4OKHyL3IriwCBaoUWIz7e660UHPueILg/4PbLyoKI0jXjrbqVP7cP36Nst3hbwu2PK6gvjCbr/93+Hq1Ttcx7a2thXj45340IfeKEo08rvwozDH+xHtZmdnbVEj+Xwes7OzKuTMjUOFyZVrLlyLL1xqlt7T02NbwGRW56kaumhjkvGjACeSrq4uFd5EcOJApnyvvj73NBy+OPbtO4Vt2657ju/ate0240KY0DlyGCFyVCmazqV78PctpcTKyorNWnvu3DksLy/7DuPic0ooZpOoOc6mm/YzmQy6u7vR3NxsaeJHzROi9p35BXFTKSVeffVVy2/vfvdP8NZbt9nOccr38vKFUVNEsjpSvzYgr0RIjrvuWgXQqM4PO2g3bHN8ENHO6drke/UbKB0GV645YgM2Fi7B1JVmcXGxcj2dYRZxgIIPjMZG3Livbxp/8zc/ZxElm5sLRpKg9yS/GzU5TKUk9uz5PaRSKZw+3YYvfvFhrK1tENy2bdfx6KPfwtRUS0n6iBeC6Hde8CvacSJ3Ct3zS2hhOMlrktiAjR3L1EeaUO4uJQS/PdjI57dnz8t43/veh5MnH7JVtwoCIYSlbS2wQfSNjY343d/9aQAvYXy8E9eubcfOncAnPvEKWlpextmzzo3qw4Jf/c7vtfwYXHp7e43tjwleHCpMrlyzxAaY2TvH2tpa2UVILxFncHDQaJn84Acv43Ofe9DVwOBlDKCG7NwFQpnpJBJ++MM/woc+9AbbCN6Pyck3LBuWrp9UI/yIdpzTl+LyCYsr1yyx6eyd9DcA6nMl4uG8RBwAnjUtTUTkJ2LGjx6jLxwAtjSkyclJxSV5l5xq0H1pHH5EuzBdPmFw5ZolNn0in332WSQSCXR2dqKxsRF79+7FpUuXsLCwUPYF4iXiRJlWAvjrGErX5r0E9E0rkUiozpxBDSdhRFw4wY9oZ3Lz8P+Xw+Wjo2aJDbBOHu/tTP4T3iWm3KKkm4gTRCwhP5GpLLmJiII4joUoNL7XRSqgUDRocXHRNVDACeUoS+A2h/r9AXOGA51XLoRCbEKIIQB/CGALgGNSyt8N47o+7w0All1dj+TXESXxBRFx9MVi4gZUQVkvS6776vi9ObwMAHzD4sTGe24HMZyEGXHhBZNoV877B0XJxCaE2ALgSwD2A8gCOCeE+Gsp5SulXjvgOIy7ehjFRIOOw4+Y6Gf3p3O5aMfvw1GqiVonUkr+5POpu1FMCzeIOBsFSNcsNo0pSoTB2VIALkopfwgAQog/B/BRAGUlNtOuTop+uXY5k55A4AtNTyL1Gpc+Ru43cjIG+DVRuxHp3Nyc5diJiQlLYqqbaBZWHGRQSLlRvJbDqVNrOREGsSUBXGafswB6HI6NBG4LhsK5ot5lnfQE3vWSxhq0UGw2m7U9r6mxRRBdkIMTKQBLa2HyCVLeGhGcW4Z5FHGQfiGEtTUyoZIBDoSyGUiEEGMAxgCgwykOqfhru+7qvb29tl02TLjpCby7qM7BvArF0qLl7X/5JkKma30u3D5z6LlewIYhoaurC52dnWpMvJgQOeNNm0OYERfFgDYzE2erRNN6jjCILQdgJ/vcvv6dBVLKcQDjQCFTO4T7WuC0qwPA8ePHLcdSAGtTU1MoepubnsLrjJi+18flFG4URkwhhx/nO40BKLSpGh0dtWS6mxZulGPmY3fSHbnOpnO2zaCznQOwWwhxDwpE9vMAHg/huoFhmkjiDIA9aiCVStn8LqXc20lP0b/nBJhMJjE8POy4+xcrGvodL+DPkGGKenESDaMaM+BtWOI6m85ZTTpblP5AHSUTm5TyhhDi1wFMoWD6PyGlPF/yyEIA32V1i15bW5vqdBOGddJJTzFxsOnpaVy+XFBzk8kkpJSW3d/0HG6fi4VfQ0YxomEUY/Zr1qeGI7SR0txms1nMzs6q91zuMuWh6GxSym8B+FYY1wob3DLIiW3nzoLkG4Z10suix53r1N4KKJR9a2xsVLoELYooX7hp3BycW/HolMbGRqRSqchEQz/wa9bv7e1VeWzcfZLL5dDe3q7WQ7n9cTUdQRIE+qIK0zrppqdwcQawEvzQ0JCtarCfFx6G6OPFreLxONbW1tTzUFQOcQaTaBi1SObXrC+EUAWQ3LL3y+0P3PTEZlpUk5OTtij3UifXTU+hndTkOP6XfxnAl7/ci9/5HWqZtITRUecXHmZFZZ7uwzcIIjQ/8Zhhjksfn/45iFnfJCLH43GbRdLNHxo2ao7Ygu6eOteh7zjC8gG5jYsInHORiYklfPObeaytNQHYKDX38Y9vc+RopYo+nCj6+vqQz+ctpdf1DcKv8cTPuIKOTydaKhPvx6xvEpGpLwKNUd90gWj9gTVFbMXunlxvK4cPyOTgzuVyluDhwcFBfOpTq5aum0Ch1NwTT6zg8cfNfbRLEX04UUgpld+OYi/5Lk9cxC1cSx9XNpu1jCuRSEBKiXQ6jb6+PgvRFLOZAPBl1ncTkXV9j89f1P7AmiG2Und1+q0cPiDTOHO5nFowhEJdEDuuXGk0JkKSKFVsKBQnikwmo3b17u5uCCFU8mljYyOklK6tqExIJpOW3LjFxUVVAmJlZUWV2PMK8wKccwH9mPXddOh4PG6s4FwOo0/NEFupuzrBTbeKcpymIrI7dvRhcdFOcDt2rDoGLQOwpfkHFX3a29stRHH58mUVVrawsIBUKoW5uTlFKHorKspxMz27CQsLC+q6ZLRw2iC9NhO/m6XpPZvcMDo3r+ts6yhlV9ev4/a5VJjGaSoi+0u/dDe+9KU9WuN44A/+YKOMgc4pub+w2Ix03ShARMQ5RjqdhhDC1orKVNdFF9v4c3P4aWLi5Y4IslnqOrRThjy/RpRujJoitkoGuAaBaZymIrKDg3vQ0PB9TEzcgzffbERHh8CRIxJ33jmNdNrcSZOg5+rpoUh8QXLi4jqaqVDS4GAhLWltbQ0jIyP43Oc+p37bvXu3UYTkYptOyCZ4EZqfXED9/l5wEy3L5SusGWKrdICrX3iNk0Av+tFHL2Hbtr+wHWsqeWAyVpCzFoAiTm5hBIDXXnsNO3fuxODgIOLxOHbs2IGXX37ZOH6KG81kMrYUm9deew19fX2Ks3GC7u3txbe//W2ll5EY6nQPp6iTKAkiahXCCzVDbNWwM/mBaZwm3xAtOD96qIlT8lQXIrbJyUk0NjYqE3dPTw9u3ryJ+fl5zM/PAyjk0l25cgUAVGHbRCKBxcVFtLS0qIyCRCJhs/hxUVhPygWA+fl5JJNJDA0N4atf/aplvN3d3bh06RKWl5ddN8hy6NRun6NEzRAbUPmdyS/4OHnKhxNHdtNDdU7Jc8t4qose+0kGGW4Q0LnrHXfcgYcffhivv/46AGD79u14+OGH0djYiPvvv9/WzZRzSzeLaz6fx/Xr1rLnly5dslSodtsgiyGIcgYUF4uaIjYg3J0pyhfENwQ3jgzYI0uI6+nnU8tcU6oLYA0F04lRJ7REImHheIlEAsPDw2rcTmUSSIR04sgDAwOWpFnijnqFat2lQSjmHZiaNHJHfbUQXs0RW1iIOuKbc7a+9UgN+swJzkm/y2azSiSj8ynDe21tzXIv0rOckM/nbd/p+tTY2Jil3IGbbgy413yhZNfBwUHLpqBbMsMKPeNcNh6PW8RonfDckM/nLePTP5eKW5LYwgh7cgMtIiIMEufoM3/xxAFoxyf97tq1a8jlcmpRk9mauIWJEGiB6613eRMJXdwECv3h7r57TVlEf+EXOvDoozByYtLVnJJfyWCTTCZt3JG6xRBhhdWN1GSt5WK0V9iYlBJf/epXcf36dVVj5ebNmzh27Bji8Tg+/elPe47DD24pYtOte+SjCjPimy8iIgwyjnBCIU5ELgEyPHD9zjS+eDyOjo4OCyFcvnzZQpR0r66uLvzgBz9ALpdTi09KiZdffll1/Ck0Ytxo6jE3B3zhCw/igQcesOnGwAYn1p+JnqOnp8fmVuBEThuK3kuulHfg5Nc0lW84ffq0JZshn8/j+PHjuHLlCm7cuIGJiQmMjIzg6NGjWF5eRmtra2gc7pYhNlO8ov5CwzC2mHZak9OY7uMVnqRnJ8zOztrGuHPnTsTjccRiMZt+17cek0h+uOnpabWIRkdHcffdaw6NGIWtEaM+Xv5MxEUo9tHUO4Cc5Dxo2KsOix+YrLUcnNB4MPLAwIBqLQYATU1NWFhYUL7F5uZmjIyMhCZK1lwzxGLAuQ2JOibfl5fu4xec4HToi8l0LOciHJOTk47PQU1E+tYj92mB0PX7+/stxhYSl5ziM02NGN3GS5yK8t446PPAwIDigFNTU8jn85iYmLAdG+Qd6DrmU089hdbWVts18/m8iqukMTzzzDNYXFxEIpFAKpXCysqK5bzPfOYz2LJli+VepaAmOVsxaTYmDgLAVrUKKJ3Due20pgBj0+IksauY0nxuFlvdOtfRIaD5rte/9/9sVP5BCOFpXHGKG+Xl8ejYoJEhxF25GM3z8gYGBtS9OTo7O7F//36bLvuFL3wBTzzxBLZs2RKKAa3miM3JguXVccUk1xOhhekg5zst12f4jgpscC89z21ychJnz260MDYZKfyU5nPbgPj3R44AY2MF0ZHg1IhRfzZOIHNzc9i9e7ctd9A0r/p7IC5b7DvgG4guRhM3op7jutgKFET1F154wTZHy8vLOHr0KJ544gnMzMyUbECrKWJzs2C1tra6dlwx7cgmYgxDZ6MXHo/H0dnZabFGUpcdIQopLdxMTucnk0ncd999lrwv3UjBcfz4ceUm0E3oeu6YvlBILzt0CJ6NGIWwlokj3xuPjhkdHbX46hobGz05Odfhin0H3M2iWx25MUQXWyl65ubNm+o7er6tW7dieXlZ6XClGtBqiticxEE9lEg3IQP+kkbDcny6vXDd9E1FaPQxOiVYOj0HhWwNDQ0pbplMJpWZ3q1k+IED/ruc9vf3q4xpv++AEGV8q58Sd1zE5Fy5paUFS0tL6lr0W0NDA7773e+q70vdjGuK2ACzOMgnzkmPKXdcpa4/mf4GNX2TNa+lpUU5qgcGBvD666/jrbfesiWExmIxJeJx/cXL5+Q1HyTyBX0HXtE0xb4Hvz47XcQk9082m7UQ2/T0NPbv34+ZmRnLfUrNMKk5YnMSQ5xMyDTRumEgLLGRjyvoonWKwnA6T0qJ1tZWzM3N4dy5c4obm6LrST8hEY/EPL7Qi43gCPoOOKJ4D343Lv3edO78/LxrCcKwOHBNERsXQ3QrIgXTEshCxkN13Kx0paDYRUuhRPq43TgbX1Tceka6B0G3uhL4GIuJ4DCJgsePH1cLlOP48eOWeEv+HG6fi4HfjUsnND8lCMPiwDVFbDQ5yWRSPTCJUZQisrS0pKx+eoFUfZLCCIItdtFyByuP5nfaPd1y2wCoiHqyaJoSQwErMRcTwaEvUKBQXTqXy9mK0eZyubIl97q5JLhTWn8f1JzExGnD5sA1RWxAYXL0xb11a+ExHnroIQAbURemiA1CWIHIxSxaKaXFwar7h/Ssa68aJECBs5EfSb9nIpHA2NiY0qmklBaXR9AIDt0A9NhjjyGXy2F+ft4Sp0mEWS5Cc3JJUJEiv24jQtgcuOaIzWtxA/YQJxNHCzMQOeii1Z+Bx/CZdmKnGiSJRAKdnZ3IZDJYXFzE5OSkMoqQ6+GFF17A4uKi4jAUKDw7O+sY7aGb6p04LF+8IyMjePrpp9Vx1EchakKje7i5JDo7O10bUHq5jcJCzREbwbS4Ae/cMPrLowlKDUR2EmH8EJwfHYMTJoG4FS38bDaL5uZmy3MAUD25KYNACKHcDabCsXwxurVI1iUMPZRpcnJSEVw54OSSMOldQdxGYW4WNUds6XQaKysrtkkgx65p8WSzWbS3t1sWDTmZOUohtKC+oyAEaiJMyj+je9BxdG36P695T5yHxjk7O+ua1OpVMNW0CegR/ib9JyqYXBK666EYl0VYqClik1Kq7iQAbMo4AFuIE4lNTuIDRzHKfDG+o6AE6hV14WVxc1qAXmZ4P7qo6drpdNqie0YpmnF4bWCluCzCQE0RmxCF7iS5XM5CQFwZ10OceA82k/gQhh9FX7Rf+5rAk08O4vJlYQx/CkKgpUZdeC3Ar31NaKFaG6k1boTKI3M4yHjDdSW3/tthwWue9ERS3afmND9houqJzfRyhoeHfSnjtCBMi4Z6RkspMTs7W7Ifhc45eZICewuf5+YKnwErwblxFf17anYf1OfjtQD/9V8HMTYmVBDyxlglDhwwcwIyp09NTalNj/vbqCZlKpUqi2hGMM0T6eUUhGyqevb6669bNl0KBAfC53BVTWxOhVyy2azlOG4E4edevHgR7e3tanFwrK2tIZ/P49y5cxYdpNTJPXTIGkEPUDKmv/hDk3l6bW3N4g7w6/Px4qCHDgnjWJ94YhWf/GRccQIqeaf7L5PJpMUYQ7X+c7kcDh48aLEKOzXlCAvUp4DX0yS9nERXfYM7c+YMtm7dagwEj8JlUbXEppu8eSEXYCP2jzqnAFYRZ2VlRb14vWYjAIsIGuYO5pR0qX/v5PPJZrNK/3SzjvkdrxsHdRrrlSuNmJ6espjTeaUsvcEjN8YABdcLL/QDWHXMsBHElcONSF6B4GGjaonNyeTd0tKCBx98UFnYenp6bDsR19UymYyF0PT6GGStCwsdHfBMxvRaHMlkMlQRzGRAkVI6Jo7u2LFqmXPy/+nl85yspnx+n3rqqaKSQoPAy/dqul8x55SKqi6LwCeEsLS0hHPnzqkUkoGBAVXuDbDqaSY/jz6JYZVCIBw5Uki+5NCTMem5KOP66aeftuhV+rj9vnz9OZyeK51OY2pqCkeOSONY/+APrC6R/fv3O1bS4vcg0Y2DrH09PT2RRpOY1orXvBVzTimoamKTUhpDk/jv/OWSKJZOpx3PpUV9+PBhtdh1gvO7aE04cAAYHwc6OwEhCn/Hx+36mtOLBsyOea/xEQHRb/l8Xs0FP4dz1TvvnMJXvpLHe97zNgCJ97znbfzxH+extGRNsDx69CjOni1kj3d3d1vqiExOTqr51muB0PwSwVXC7O/27oo5pxRUrRhJE0HWLQC24Np8Pu8oilGTB8BeL5Em02TVCyNmkidjOoU7ub3oYqI6uK5nKlQ6NTWFpqYm9PX1Kf9j4bpn8d/+24Yr5Ic/LIyltbUVIyMjOHbsGBYWFtDc3Iy7775bzWlraysuX76s0lNoHpPJpHpGmt9sNoszZ85ERmx+zP56CBzNd7EulWJQtcTGLWmmgqC0a5pbMRUcqy0tLbjjjjssHCSbzaKpqcnilNWV5rBiJt0CX3kiJ7+PW+0RwD2qg3Q9AgU5nzhxQpVf4GI2h17SnOqC8EaIzz//vPrdlB/X29urgg70oIL29vbIfGxuVle9dB7fPMudUFy1xAZYe2Hr5n4ppQq25QYQOn5tbQ1LS0t4+OGHAVh1OD3Cgv8/LKXZi3DdfGem2iN+ojoA53r//Fn5Dk4YHx+3fJ6amsLQ0JCxtwCHHg7Fw8PK5WMD3DuNOr0D3Zrq16VSLKqa2AjkQKWEUQqg7e7uttW956UBvEq8meAWNREEfgjXb76UaUPwE4TNoZe/A2ARzynxlH/HF6wT9HCosOavGOj3KKbqcpTjrHpi00UEqgjc3d2N+fl5VVpbrxEIhFddt9jwHa+FZzKdFzs+ruvxtlL8PA7auADYXCE0lng8bgu1cmoxrPc10MdXLoLTUUni11H1xAZYxUkSzXp6enDvvffi5s2bWFhYQIehqqhp1yWY9IdS4xB1hEm4XuMjXc9UqJTSbDh4lIVpfDTGdDpt2ehMcahcB3XqAVDM/IWBsN9BKagJYgM2dn2To5sWwNmz9tokwMau62VhdFO0gyrNYROun/HR4ueGJSGEWmzJZBLDw8O+i9pIaa3DyKMqdAMTbYZ0vqkHQKUILYiVMsoxiqh8Cm545JFHJLdsEfxwH/qeByIfPnzY0STuZvlzkt39jsMLUfWA88ulvdwZxHlM4wMQeOym92KSKsqJYp69lHcjhHhBSvmI8bdqITa/C5PvVgSTZYmO5dcynRP1IgiLcKMah+l7wNkH5WTgqeQcm6CPDdhYG3raj58N2C/ciK0qxEi//i0/ohmHbumrBgtZmPcMQshO43D63mTFI78d3YuXSCi3g9gNtHFzMZGknP7+/lB7wwVBVRCbX/9WsTpVNSnJxUInJN7Uj8AXf6nPZdqgyEFMjmrKZQPKX3HaCXzjnpubs5Q9oMz8/vX2WeXegKuC2AD/3Icr7Pw8L0Krll23GJjy+igci+qrkFk+lUqFphfqGxS3LALWzTDoe4kK5Bsk4w8548k6SgnDgHNxqKjGXDWByEGCQoOIZk7cMOoo9LDAd2qaDzLtJxIJVR6CCI2sstS/uth78g2KHN16TqC+MKMUmYOAol449Bw8/nxuQelhoio4W9Tcp1p23WLgJmLrOWZOkTNOcNL59A2KX5ujWkVxk/8QsGaLV0LsrQpiC9O/5XYPt8/VDJOIbSr1QPDjztAbuevWX9qgAHsFZs5Bne5XKfCNm2eXA8DExIQKsK7EBlw1YmRfX5/NGDI4OFiS3rFZIKXE8ePHLd+Nj48jk8mgra1NGSkIU1NTqqWUU67bhQsXLGITLVBd/KQ0J6o3QpnuFAhu2gx1Mayc7iUhBBoaGizV0z772c9airE6ZT5EvWFUBWcj1DL3iQpSFpJgKVdNj1G8efOm2sWFELh27ZqyxO3evRtra2vKcggUuBAP5wpi/SXQ9yarZ1SOfL9Ip9O4fv06du/eben6StXUKqmnVxWx6agWh3AlIYRAU1OTTXQDCuFX9957L4QQiniWlpbQ3NxsKc6TSqXQ3t7uqfP5sf46HQcE85dG8R75/XmoGP/Mw7PKjaoltkrvkNUErj9xYqPeZ7zGPQAsr9en0y1wus7nNzrfr8ThZMwhZ3jU79Gvv7ZSqBqdjcNk7nbSKWodQfQbJ9cIj4jQ4ZTrRr3gdNP35OSkZQym8biNmS94AteVon6PpvtXA6EBVcrZwtyhohJFw7huMfGgbqW1TZicnFTip34u5ziDg4OWOibARgwhH4+fRvFOzvBycBq/0UKVUFFK4mxCiP8ihDgvhMgLIYzBlyVcu+QdSrfC0YugilPFIozrBuHebo55nuBJ4Uj8byaTQTabNZ7b1dWldJh0Oo0bN26obqH5fF5xPxqP15ipsJDOMU3OcNN8lAp9U3JyVke1LrxQKmd7GcDHAXwlhLFYUGo8o19lvZhxhXFdJ+6tV2im67n5hSjBMx6PWyxw9NmrromUhRLni4uLNi7EuR+VrNNLLCSTSXU9Uz19vSI1tfeieiUhprd4+mv196cbUfL5fGRGlJKITUr5fSB8E30YESVRKcthXtfkrNZFG74InQwVnBDprykNxnSu6Zk4yAlMi5SXFiQkk0n1f30s5GbQrYO5XA5SFloEh7EJmu7Pn03/TM+qW2d1sTlMlM1AIoQYE0I8L4R4/s033/Q6NnA8o0lpj0pZDuu6Ju5drFGILybTX9O9dVD6DAcZNuiZ9fLtpnuY/HQk7lJ15La2Npw7d85WDTqMjduPNXX//v2W76jqc5TGG0/OJoT4WwCthp8OSSn/yu+NpJTjAMaBQvKox7Guu7UOr/qMHCSKrj+b3+HbxldqxLiJe1PVsKiNCab5mpycxPnz5y3H6fVDAPuckUhJx5j8dDxZE9gQ3ebn59Vx5bQYnj59GhcuXLB897nPfQ5AtMYbT2KTUv5s6Hd1QVD/mpsORSE6yWQSyWTS4hSWUqKpqckWBeElxoQVNG3i3qZqV2G/eKf5onu2tLTgoYceUuNIJBLIZrOq4f0rr7xiu14ikcDCwoLjOGOxmLGHOUe5gpp5ehJVfSZCAwocLqoxVJXpvxjjgxBCtTbSFfuGhgZLeFMqlUJPT48ycSeTycCOcz9KuF/o+gVdnyPsReiksySTSbS1tUEIoUripVIpNVfU8P7tt98GYA9GpmrLbpKH3k63ra0NIyMjZc0vjMVi6OrqAlDw/3FCA4CZmZnKcTY3CCH+E4D/G8BdAL4phPhHKaXZu+rveq7Kq4kDkSVNNy9zpVwIYVkYgDly3a+i7qWEB31mep5yJbmaDDPDw8OW301JotTwnsdmArBYFnXwDVRvp5vP5yGlLEt6C0d/fz/27t1rIbTPfvazmJmZiZToS7VG/iWAvwxpLADMC4GUaycOZDIvk8mazLj6NXlbJi8dyY8DNAyjSxRpRqaxA+5Zyvpc0fe0yQBWUZfCxpyey+n98Dr85dTZpJSYmZmxfDczM2NpC1x1nC0KmIwPvOOozoFI+TZxNnqRZILm4MVI3UoxlCtGM6hRyOtaQghb4RsyhOi9sL2iUaampmzFc/Tf3cZJoht/R1QbpBI1Svj60fsBRBmsXFXE5iZKuYX88Ha0hNbWVsTjcaTTaUv7JCrNffbsWSOH4gsnKse4jlIIWh8DJYUODAxYxDcqz85z03QOqpcb5887NzeHkZERJWqZSr67ZQPoVmFyA5Q7Ct9LgohyPFVFbG4ToXer4QuTdDZ9gXR0dECIjfQTXpqbrGzUX8ykI/G/UZnjSyFot0JAACzim97iieab/urRKKYIkJmZGcvGRtEqgLPoVU5d1C/C1LmDoGqKtHLoC0z30wDWBe/GGfbu3Ws8d2BgAGfOnHE8D9ioCAzAWOk3LPAFaXo+r3P0Rdzc3KzSbHT4GXvQ+ffi8LdSupSo9iKtOnSF3km0AQo7o9NONTs7q8QV3eDiVocCsPYw0zekycnJkhvf6/fVx2iKj9TPNXFdXZzW4ceVoP9GKTxO4/Oah0pxEj8oZ/R/VeazcTiJlhQCxF+gThQknk1MWHtEX7hwQdXoMFkV+T3Onj1rqVpF4UqllDwz1QTRx0i/nz592jFCnRMcYWRkREX9E5qbm/HZz3626HJtFM1vGp9fhG29DQPljv6vSs6mw7QzUigWt97pQbtc33j4xRex79QpbL96FVe3b8dLP/gB9nz+864ma9NuTr8Va0UzRZ1TL7WWlhY8+OCDys+Vz+dx6dIl1ajQyRLLQT2wOZaXl32btvWd/Tvf+Q5ee+01o4GJxlSstdTpczlQLuMXR00QG2AXLXUrmGmiyOT83ulpfOS559DwzjsAgDuuXsVtX/wixPvfv9FpXoNb/GMUWQOJRAJbtmzBuXPnkEqlVPTG4uIiWlpabKIiD5zVLYO6zkaWWRIHgxhciNB0AxNdM6j+51VCr1xweg8VjY2sRvidKCLKfadOKUIjbFlbAw4dMhKbHwuaHzjt4CauOTo6ipmZGczPz9ui6h966CFbr2yTZfBnf/ZnHXuuUdltN46m7/ScsHgp756eHuzfvx9btmxxfFbA21paDm7iBjcHfhSoSWIDvCeKE8zgtWvmi1y65HjtUqM53Cxw1PWFg4t5eliZDlPNfeI6DQ0NlgRUcup7jdttA9OrcMXjcUsMoYk7eRFvubiJE/xE0oSNmiU2rzQXTjDYudNMWIbWwIRSLGhe+gCl0ujcx2RwmJubw+LiopHDEngUBE/kJPhdOKYNzBRR4oc7BSHechMabYQ8QF0IgWw2W72xkZWCX0epIpj/8T+AsTGA6TGyuRniyBHLNb0sZkEXLWAWcymgV+c+5GTnWFxcRCKRUAHVpjhDWrg8uJquG0RMM21gZLxpbW21tF9qbm72xZ38EG+50msA60ZIfchpznp6epBMJm+d2Eg/CCLmCSE29LJDh4BLl7C6Ywde/aVfwp7HH4eAvfxAWGP0CujVueaJEycAwLaBbN26VYlrepkBJw4RVEwzbWDHjx9HLpdTRV+np6eRz+fR0tKCpaUly/lOMaU6tz569CiWlpbKEk3ixz/JK03rkTVho+r9bE7o6wvYG+DAAeD11yFv3kT6q1/F/9vSUlT5Ab9wEnPp+ibH8X333Wf0J95333020djEIeizbsDxs4hnZ2ctVbiklLhx4wYA4I477lA+unPnzmFpacnmy9P9hsQ9MpkMuru7C+I8gKWlJQunpmcMm5uQD433PJiamsLp06eN/kmufkSFmuRshGLEvHKYfIuNLPejJ/b19dl6QvPrewVXO413dXVVJdQCBcIlEXZ+ft4m3posngCU/3NgYEAlnp47dw7ARp8CHvAbRTSJnkOndx+VUuL69euB5igM1DSxFYuoTb6cS1AkfT6fRzKZVJ+dRFY/G0gsFrOJ0f/6r4P48pd78Tu/04jt269ibOxufP7ze3yJafoGxCNmdFGVwPtr07lEaHS/gwcPWs7libxez1gKKOiBxF8ag94DodyB0bcksXlZMku9NnGJtrY2y+Lr7u5Wn0v1K3EuePIkMDYmsLzcBAC4evUOfOlL2/H+9wt88pMbBhU/pn/dmKHPUyKRQGdnJzKZjKq81dTUZHEDAOZaI5RfGDXIv6oHYy8vL6sAgM7OTktmA2CtLVk3kISAqFM+aFdtbW1V4hOBi1NhEDadf+iQxdAKAFheFnjySYk779wwqHg5tZ0skVQwidwQACz1SUgfozHpQdVPPfWURYyOWlxzstoChdjRLVu2uCbYRhXRcssRWxgOazfw/DonhL3YHHzzuHQJlk3Fi9CcknYPHjyodnwiQCI6U9SOyXATdckB/XlM2ftAIXZULzxbrvjIW47YgGhTPtx2VULYynhHB6DV0gEAbN9+FYA10dZJZ3NL2uXGjNHRUUd3QyVLDujP09DQYIsR5S4M/TnLEdFSs6b/UhGVku62q1IpvWLSXPR78HOPHJFobrZea9u269i37xSAjba/bukjTq6U/v5+27Nx6O4MTrQU+Eym/XKVQCBrI+loTz31FHp6etRnPTWrGFdJMbglOVuU4DobJ7hEImHRA3Rxym/aSTqdxsWLF9He3q4WyXveM4lPfjKO5577d3jzzUbcfvtVfPzjz6Oz82UVh8gDlN04nNPnwFE7EUgNfqH7I3n/ung8jv7+fstmVa74yDqxhQyus+mLcteuXQDsO6eedmLqi0bXXllZsfTIpvi+nTuBZ5+9bvFzTU/fUMTAx1PMIgoctaOdW27oRD87OwspCxXM0uk0VlZWIKXE/Pw8crlcoAJGxaJObCEjqAGG6s7zIj1khNC5kBCFEuUUQ8ktfiQe0a4dhS+RHOp8PJXuU+0GrktevHjR0uiRfIm33XabKp/up4BRKagTWwTwK0pxLujWF03Hzp07jREdPGctCl+iKW0oyhZLYSKZTCKXy9lyBd9++228/fbbltjIW0pnq0TafNj39CNK6dYwDjJPm0DxfhyUvc0JLUxfYiXKCIQFkggA2IiNQO8hyqzxqiO2SpQ9q2SpNZMTGIDFPE2QslDRWHeWAwXOduzYMXR1daG/vz90X2K5zeRRwG2Meu/xTV+DpBK7Z6V3bFNlLb0vGhdHm5qa0NbWhh//+McWH9LWrVstomQUVsEo9MBygHN6AJbiTz95z3sw3dcH4ne3TA2SSuyeldyxuU9OT850KqjT29uL1dVVzM/PW9wLN27csOl5YVsFo4wpjRrZbBYA8NGlJfz0N79ZqEED4La33sJHnnsOAPDynj2RPkvVmZHK6WSs5D3pviTukY5GTmASB73O4ejq6orMMqjrgYcPHw7FQV8OCCFw3333IZVK4f1f/7oiNELDO+9g36mNAIConqWqOBtQmd2zkjt2MeIemeD1aA5eRzNsRB1TGjXUPF++bPx9+7VravMAboEaJFFH5FfLPXUEFfdI/Cz3mKshOqQUCCEcA0lFR0fkm0dVEVslds9a3LErOeZqiA4pCUeO2Io/obkZOHIk8s2jJrrY1KKfrRyoxTGXE47zc/KkKv6Ejo4CATpUxg4KUctdbEyfN8s9S0UtjrlccPWdHjgQGnEFQdVZI+uoo1Rw32mUFdSCoio5Wx11lIJqjXapc7Y6Ni0q4Tt1wy3B2eqGhFsLlK+m4/jx4xgeHq5ztqgQqLvkyZPArl1ALFb4e/JkOYdaRwig3LVMJqOaZlAnoFwuh8nJyYrpbJua2AIpyoXiiwWHp5SFv2NjdYKrQfBOPkR0wEZpCp5UyhE1EValny1M6BHfgIOivGuXuURVZyfw+uuRj7OO8ECpSDx3jXrWUexoVGlVbn62Tc3ZgABBxm7FFw0o965YRzCYfJCco1XCNbDpDSS+g4ydii8aGiZWMtm0DneYJBkAtrjRSrgGNjVn4xOfSqUsaSE2RfnIkUKMHMd6zJx+zWp0mN4K8CtNUO4avXMykND3QGXSqmqfs7nEuVHAbjKZVJM4OFjoPZbL5TA7O7vBibSGiU4xc9XqMN3sMEkTk5OTaGpqUu+QXDr33XcfkskkhoaGVP0RKSWam5vV+zGlKEWdVlXbxEYWRIrgJgsioIiEMpt1MYJ6kXGfm3z8cQiNuJx8dLVYHqBWYSpdQV1RU6mU4nBclNdrRsZiMfT29gLYKEVR7rZRtU1s5vYthe8ZdzNxIs7tAHOrXyfdjAqhckSxKwZxxm9mx73TO6TfADjWjTERKi9FwTugAtGmKNU2sfm0IJo4UXt7u2uRHwCOhYCo9kfYuyJfJEGMMH6PrRaCLGYTMb1D0r/dRHk3sZ+3ieK5bFHNS20Tm08LoskiKaX0fFlOLykej1ua6YWxK3KCATYIPZvNYnh42LgZ0MLwUx2sWiyoxW4iQKHcHIduIHHa6NzE/nLOS20Tm0vWLcGt7AFZqQhOu6LpJRVbHsC0qwN2Lkrf53I5PP300wA2elI/++yzaG1tVfekFlVOG0cp5frC5IZBxqEfK6W0tB/O5/O2gqtOojwZUzgmJycxODhY78/mGz4siEIUusrwxTcwMKAskhz6ywpSCMjPS3HbRU1cNJFIqKaDNJ5MJmOrK2lqUaXXmyzGghr2ru8l0unj4McSeN83oKB7c87Pz6Mxc2PK0NCQii7J5XI4ePAgpJSW8VC0SdiiZO372Q4cKIRT5fOFv5o1MZ1Oq64ufMFQJxinsmxhl27z8s8BUE0dCJ2dnZbPFFg7OjqqxvL0008rPZJDH2NQv5Iff2IxUTSmccTjcdWfm64zNTWF2dlZ45hjsZiqwUJR/AMDA6oP3OzsrC34HIDaXPlmOjs7a9w4Z2dnPZ8lKGqbs3nASWzJZDJIJpOeBXPCLKpD59/9ne+g8+d/Hrh6FT+zfTvuHhvDnnWxUa+MTMTFxSUnY4GXwSYIl+bjBczccHZ2tiiuZxoH7+Kji/lOY+bmfd4XWwihzid9lxqRcJGc+oSvrq7axNG6GFkEvBYMHcOP5ZPb19cHefIkxD33AJcuFcqdHTli88X5Hs/XvoY9X/oSxLqOecfVq3jXF78IuWcPJpaWLJWRT5w4gVwuhznNAJTNZpHP5zEzM2P53s2M7aa3AsENC4CzpdZtkfJx6BWgm5ubbaIcAM8xmzZUJ33XZI1OJBLGdxVFJNCmJTY3k7G+uPix/LM8eRLCw2keZPeTTz6pCI2wdW0NP/71X8fCZz6jFmAsFsPBgwcxPj6OxcVF9T1x5aNHj2J5edm2CPVe0XwjKYZLu3HDYnRA0p/JdUKN7U39xwcHB3HixAkkk0nXMTttqLpEYBJ5dQIkpFIpNDU1ha6zbUpi8zIZc/HJyQiwsLCAT/zmb6JJIw6xvFwgmgMHAhkMpJSAUzXeq1cBbLSJklIiFouhq6sLQghL0/VLly75csjqC4WLXfS7l87mxQ2LiaLp7+9Hb2+v4mh0vt5sniI8uru71XdOYzZtqG2nT+O/j49j+9WruLp9O07t24fUyAiGhoYcrdH8ehRtEiZKIjYhxO8D+AiA6wB+AOCXpZQ/DmFcRcPLZEyWJ6BgkHBzXDdeuWK+yeXLtsXoJjrR92s7dqCRWRcJV7dvBwC101ODwQ996EPo6+uzLUzukAXCrUzMubwbNwSK70VNPQ04cSwvL6OlpQVLS0sACjpoIpFQxgrazNzEU8LDL76IB557Dg3vvAOgIK5/5Lnn8P177wWGhpSoqVuj9fURtkWyVM42A+C3pZQ3hBCfB/DbAH6z9GEVD12sIOh6WmNjo6WxuWkxrx4+jCYDwV29/Xb8IdMFnF4K9cpWfa5///eRP3gQsRs31DE3t2zBqX37HBvNmxZmMa4HP2Z8/Zje3l5MTk6qxc7ntpSS7SbxNJFIoLOz0yL6bd26FZlMxrceSGP5yXvfqwiN0PDOO7j/T/8U4vd+T12L4mMpdI+ukUqlqq/Nr5SSh03/PYD/XNpwwoGXnsb/T43NOaSUOHPmDHqPHsX1X/5ly4uTzc04tW+f8bocvFc26Sjv/NM/4bF83uJvicViuO/ee/Hyus6iN5oPakU0wY8zGTAbPfTFXqql1ks85SDXjCmQgGDiwO/60Y+M9yZJRT/HckxjI3p7e0MnNCBcne0ggK+HeL2i4bVAuSFkZWXFZvrNZDLo7u7G8ZUVvPsjH1GN81Z37MD/9x//I15m/i/TwpfS2iublP8Pfvvb2Kq16BXvvIOuZ58Ffu3X1Hc6oRXLQdQ9PKyyfBPyOoYQVAfkY9GJgxtJnAwbbvqxbSxOYXyxWCFT5MAB2zn0/FEQGcGT2IQQfwug1fDTISnlX60fcwjADQCO1XGEEGMAxgCgw5D9HBaCLlCnyT1//jyWl5eR27MHL+/ZY4nmSCQSGBsbU7qUfl0nUZYMITp03ZATsBcH8RtO5ccq6+cY/Zpun52gL3QyBnV0dNiucf78eRUl46YfWz6bwvgA4OZNizW52PEXi5IL/gghPg3gVwDsk1IuexwOIPqCP8UEu9ICe/jFFzGQTuNdb72lrFgv79mjjm9pacETTzyhrIZu1sh8Po9nnnlGff7vX/gC7jAQ3I+3b8ff//mf28Q7J/GpmOBivgkRdK7l55ggCBJXqW+SAwMDyvVR1FhOngQ+9akCgemIsIiTiKrgjxBiCMD/AeDn/BJaOdDX12fkNCaC6O3tVaLKwy++iI889xxue+stCGxYsR5+8UV1PBGa5bq5nK3eZD6fx/j4uOVe333sMVzfts3y3c14HHNjYxbO9bGlJfR9+tMQW7YAu3ZBfO1rlnOIsPyWZ/ATehZ2eFqgep2wi5exWAxjxIXWEYjoDxwohPCZ4JSaFTFK1dn+CEAcwMz6JPy9lPJXSx5VCPAjItACICPA4IkTEAYr1r5TpxR3445joBAVomeLy7ExvPS972GxpQWJRAJbtmzB/Pw8/uF978N1pgNe3b4dc2Nj+Onf/V3L9XiUicmRTs/jV8fyI47S/4txfOvcC/AXYaKfS74t+q3ksgUBijiVA5u+bqQXLA7wLVsgDPMhAUx+61sW87B66Q71JlcTCaS/+lXleKawIWAjqlzvHup2PSfRR0ppufbhw4cBwLKYnRa3SaxzOsb0vVt8ZG9vr6tI6iYC83N1vTuwKGlKwRoft2xcQcRdL7iJkZsygiQILMq6w064umMHmpqa1M5r2ekdRJLGK1csPikOsogaOUeA+pUmq+vx48dViNP09LTy8Zn0OCdDiv7ZqdgOZU4AUM54zr2cjC1ergiak5KDwH2kYJUzefSW52wWGHZCub4TCqdYSBdOJP/5nz13aAC+r8c5mymChZvMyXqql3AIumO7WXdTqZSF4ICNKIxsNgshhO03k1vD9DvdOyyOE/TZip2vOmfzC8NOKAzJqBZ49GgOvEP7yD6ncZgcs0RwZMUrhdDoPm5RNpQbSNCJ3sn9Qn+9XBH6WMJEEL03DNSJTUfQFrAeooruU/J0/vqsX2m69tDQkM1BD5TurHXyv9FvHHT/trY27Ny5E4ODgypKh5I7ASgxlIMMIPp1w+ZoHEF9i6WgTmxhwINAA+/QAQhe95GZUGqZPdO1JycnLQYjPdyK90GjBM1UKqXiLXlwOOd82WwW7e3tZdGhnJ6t1PlyQp3YNgG47kE6Gv9bShS7m16TTCYtRVI5pqamVEViLqoRkSWTSQtRDQ5uROIHTUotFmGFw/nFLUFsUSvalQbX36jMHlkH6XOxUexuemc8HreUFaD8MN4TjRMc537Dw8Pq+vR3aGgIwEZGQdQ6VFE6dQnY9MRWLfUSowbX3+ivW7R8sdcGrHpnOp1WaSpELEIIZLNZle3sJqpx0PXLpUN5PVvY2NTE5ie1ZLNxOLe/HEG5vZPe2dfXZ4n8AOyl5IKIauXUobyeLWxsamIrt2m3VhBFPUinz0FEtXLrUOXGpiY2oLym3VpAJbi9X1Gt3DpUubHpI0jCThvZDKj2Oallg5ZbBEntV0R2QdhpI5sFnGMQqoXQgPLpUOXGpiY2J7GEIhk2y0sMCicjxK26+ZQLm15nK6dptxaw2Y0Q1YxNT2zA5hVLisFmN0JUM24JYqvDijq3rww2tc5WhzPq3L78qBNbHXWUCXViq6OOMqFObHXUUSbUia2OOsqEOrHVUUeZUCe2OuooE+rEVkcdZUKd2Oqoo0yoE1sddZQJFclnE0K8CcBQ9jcSvBfAv5TpXvUxVP8YgGjH0SmlvMv0Q0WIrZwQQjzvlMxXH8OtN4ZKjqMuRtZRR5lQJ7Y66igTbgViG/c+JHLUx1BANYwBqNA4Nr3OVkcd1YJbgbPVUUdVYNMRmxDivwghzgsh8kIIR4uTEGJICPGqEOKiEOK3Qh7De4QQM0KI19b/vtvhuJtCiH9c//fXId3b9bmEEHEhxNfXfz8rhNgVxn0DjuHTQog32bOPRDCGE0KIK0KIlx1+F0KIL66P8UUhxL8Neww2SCk31T8ADwC4H0AawCMOx2wB8AMAPwWgAcD3ADwY4hh+D8Bvrf//twB83uG4n4T87J7PBeDXAPzx+v9/HsDXKzCGTwP4o4jXwV4A/xbAyw6/fxjAtwEIAD8D4GzUa3PTcTYp5fellK96HJYCcFFK+UMp5XUAfw7goyEO46MAnl3//7MAPhbitd3g57n42P4XgH0i3JoIUc+tL0gpzwB4y+WQjwL4U1nA3wO4Qwhxd5Rj2nTE5hNJAJfZ5+z6d2EhIaV8Y/3/CwASDsc1CiGeF0L8vRDiYyHc189zqWOklDcAXAVwZwj3DjIGAPjf18W3/yWE2Bni/f0i6jVgQ01W1xJC/C2AVsNPh6SUf1XpMfAPUkophHAy+XZKKXNCiJ8C8B0hxEtSyh+EPdYqxHMA/kxKuSaE+BUUOO2HKjymyFGTxCal/NkSL5EDwHfT9vXvQhmDEGJRCHG3lPKNddHkisM1cut/fyiESAP4NyjoO8XCz3PRMVkhxFYA2wH8awn3DDwGKSW/3zEUdNxyo+Q1EBS3qhh5DsBuIcQ9QogGFAwFoVgD1/HXAD61/v9PAbBxWyHEu4UQ8fX/vxfAfwDwSon39fNcfGz/GcB35LrFICR4jkHTjX4OwPdDvL9f/DWAX1q3Sv4MgKtM9I8GUVtgyv0PwH9CQf5eA7AIYGr9+zYA39KsURdQ4CSHQh7DnQBOAXgNwN8CeM/6948AOLb+/38P4CUUrHUvARgO6d625wLwNICfW/9/I4C/AHARQAbAT0XwDrzG8H8BOL/+7KcBvC+CMfwZgDcAvLO+HoYB/CqAX13/XQD40voYX4KD5TrMf/UIkjrqKBNuVTGyjjrKjjqx1VFHmVAntjrqKBPqxFZHHWVCndjqqKNMqBNbHXWUCXViq6OOMqFObHXUUSb8//d8HzZpObF3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# generate new data\n",
    "inner_sidelengths=[1, 3, 1, 1, 1, 1, 1, 1]\n",
    "outer_sidelengths=[2, 4, 2, 2, 2, 2, 2, 2]\n",
    "_, data_matrix = generate_concentric_highdim_rectangles(inner_sidelengths = inner_sidelengths[:2],\n",
    "                                                        outer_sidelengths = outer_sidelengths[:2],\n",
    "                                                        verbose = False)\n",
    "n = data_matrix.shape[0]\n",
    "\n",
    "num_rand_seeds = int(0.05*n)\n",
    "x0 = np.full(shape=(n, 1), fill_value=0)\n",
    "random_seeds = np.random.choice(np.arange(n), size=num_rand_seeds)\n",
    "x0[random_seeds[random_seeds < n / 2]] = -1\n",
    "x0[random_seeds[random_seeds > n / 2]] = 1\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (6,6))\n",
    "# formatting\n",
    "unlabeled_idxs = (x0==0).reshape(n,)\n",
    "plt.scatter(data_matrix[unlabeled_idxs,0], data_matrix[unlabeled_idxs,1],marker=\"x\",c='grey')\n",
    "pos_idxs = (x0==1).reshape(n,)\n",
    "plt.scatter(data_matrix[pos_idxs,0], data_matrix[pos_idxs,1],marker=\"o\",c='red')\n",
    "neg_idxs = (x0==-1).reshape(n,)\n",
    "plt.scatter(data_matrix[neg_idxs,0], data_matrix[neg_idxs,1],marker=\"o\",c='blue')\n",
    "ax.set_aspect('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## AUC Histograms for varying dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n",
      "Using degree regularizer\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAFICAYAAACREMOwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvL0lEQVR4nO3deZgU1dn38e/NKqNhUVZZAtEoAeI67vERoqhgMGwKJj7KYlCfuCbBiCIOmmg0xhAN+kqiIoniAoobIC4gCioawQUEgwEV2QRBRIb9fv+omklPT/dMT9HL9PD7XFdfTZ86axfTd9c5VdXm7oiIiERRK9cdEBGR/KUgIiIikSmIiIhIZAoiIiISmYKIiIhEpiAiIiKRKYhIysysyMzczNqnmH+5mc3KbK/Sz8zah+MsqiitpjCz8Wamc/0lEgWRPGdmTcysOPyA+99c9ydXzOyIMMi1z3VfpHoIv8R4BY+/5bqPNUGdXHdA9tjPgfrAMmAI8I/cdqeMQ4FsfcM9ArgRmAUsz0D9nwINgJ0ZqDvXfgFckutOZMBVwH4J0n8JHA88m9Xe1FAKIvlvKDATeBoYY2bfc/f/5LhPALj7tlz3IV08uLXD1lz3IxPcfQewI9f9SDd3nxKfZmYNgL8Cq4Cp2e5TTaTprDxmZkcRfAN/CHiE4FvykCR5PZz7/rGZvWFmW8xshZn9NtzexMzuN7O14bbnzOzAJE3va2Z3mdnqcCrtLTM7NUGbCddEzKzQzJ4ys3Vmts3MlpjZ9WZWJy7frLCOA81sopltCPv2gpkdEpOvCHgwfDkzZrpifMXvIJjZj8xsTjiONWb2VxJ8e61sncTMzjWzBWE9S81scJinnZlNMrOvzOwbM/unmX0nQf2tzOxeM/vMzLab2UozG2dmzePylaxLHWpmt4T7cJuZvWdmPRPUe4GZzTOzjWb2rZn9x8weNrNmMXkSromY2WHhflpvZlvNbJGZXWNmtePyjQ/71Cgcw9ow/xwzO66yfZBl/YFGwEPuXhOPKrNORyL5bSiwGZjs7t+a2XPAhWY2yt13J8h/JNALGAdMAM4F/mBmW4ELCaaBioCDgSvCPKclqGcCsAu4DfgOcDEw3cx6uPtLFXXYzM4CngSWAn8CvgJOAG4iCIjnxBXZF5gNvAlcB3QArgSeNrMu7r4rrK8VMAy4BfgoLPtJJX05DngJ+CYcy0ZgYDi+qvgJwXTQPeF4hgIPmNn2sD+vhH0/hiDIbwUuiulHO+ANoB5wf9jvg4FLgW5mVujuX8e1+RDB0cMdYbmrgClmdoi7Lw/r/d8w32vAKKAYaAv0BJoDXyYbkJkVAq+GbYwFVhP837kNOJxgGjXeC2GdNwEHAL8CnjezDu7+TbK2wvbqEny4p8Td16WaN85QginW+yOWl3jurkcePoB9gA3A+Ji0nxL8gfRIkN+B3cBxMWn1CA7rdwN3xeW/MyxzaExaUZj2FlAvJr0NQTD7KK6O5cCsuD6vJggKdeLyXh3W3TUmbVaYdk1c3uFh+hkxaYPiy6fwHs4FtgOHxL0n88K6imLS21eQ9i3w3Zj0ZgSBYjfwq7g2nwzb3C8m7WlgLdAmLm8hwdFlbJsl++A5wGLSjwnTb41ra1P8e53gfRhPOGMXkzYnbPuwmDQDHg/bOTW+PHBPXB3nhOkXp7AvuoZ5U3pE/Js5ONwnMzP997k3PTSdlb/6Ao0JvmmWmErwTTDhlBbwhru/VfLC3bcTfGAacFdc3tfC5+8nqOfPYdmSelYADwMdzewHFfS5O9CCYOqpsZk1LXnw3/np0+PK7E7Qt1cq6FtKwmmiE4Cn3f3jkvRwXH+uYnVT3P3TmDq+BJYQ9H1sXN7XgLoEAQgza0RwJPMMsDXuPVlOcMQW/54A/MXDT8awzbcJAnnse/I1UACcZWaW6mDC9+ZE4Bl3fz+mDQd+H77sk6Bo/PtWlf30HsH/j1QfUQwl+L+uo5A00nRW/hpKEDBWmNnBMekzgHPMrKmXP+RPtOC+IXxeliT9gARlPkqQtih8/l6S7QAlAeaBJNshCDKxVrp7/IL2+gr6lqrvhc+LE2xblCCtIsne11Ve/uSC+Pf1UIK1yaHhI9X6E6Wtp+x7cgvwP8AUYL2ZvQpMAx7ziqeXOoTPCxNs+4ggOH4vwbYyfXL39WHsqnQ/ufsGgqnFjAjXcS4kmLKclKl29kYKInnIzDoA3Qi+VX2cJNv5wJi4tF3J6vRgbSFhc1XtXwVK6hoOLEiSZ2Xc66R9Jr192xPJ+phK30ue/0nZo8pYxVWou/Q9cfd/m1kn4NTwcQrwN2C0mf2Pu1e4ZlRVe/J/yMzqAftXoa3VqeYN9SRYNxub4EuJ7AEFkfw0mOAP8xcE36zi/Y5gSmtMhtr/AcH0Q6xO4XNFpxf/O3z+1itZgI+gqtejlBx5dUywrVOCtExZStD3ehl4TwiPhKaGD8IzuJ4nWPT+ZZJiJe9N5wTbOhIcOaX7NPITCU5VT1VVv0CUnMjw9yqWk0ooiOQZM6tFsIj8gbsn/IMws85AkZkdE86Vp9vVZvZkybqImbUBfgYscfdkU1kQnL2zFrjWzB5z96/i+t2AYBG4wjN5ktgcPqf0bdbd15jZm8BPwzOaPg77UI9gkT8rwimfqUBfMzve3d+M3R6uZTQN11mqJMmU5rvhc9L3yd3XmtlcoFd4BtyHMX0ZEWZ7qqr9qUTJmkjamVlLgiORd919QSba2JspiOSf0wlO06xocXAywVk8Q4FMBJE6wGtmNpHgFN9LCK7mvqKiQh6chnwBwRz9EjN7gOCbeGOCb7h9CRZsZ0Xo09sEc/XXm1kTgjOmlsWeSJDAr8K25pjZWP57im+2/y4uBV4HZpvZBGA+wbf97xGccTeBYH9W1Qwz20iwmP85wfs8iODIp7I7G1xJcIrva+F7s5rgBIAzgEfc/eUI/Ukqw2siFxLsUx2FZICCSP4pWXx9MlkGd//QzD4GBprZ1e6eaE59T1xAEDiuJfhgeh8Y5O4vVlbQ3V8ws2PCsucTnA67geDaiDvDuqrM3T8zsyHAb4F7Cc6AeojgdORkZd4ws+7AH8L+fE2w6Hov8EGUfkTh7p+b2dEEff8pwfuyleCD/1mC02qjuJfgWqCLCY481hMEqMvdvcKpI3d/x8xOBEYD/0dwvc5/wj7+KWJ/cmUIwbrSI7nuSE1kMWcJioiIVImuExERkcgUREREJDIFERERiUxBREREIlMQERGRyPL+FN8zzzzTp0+fnutuiIjkk7TdMijvj0TWrYv6swIiIrKn8j6IiIhI7iiIiIhIZAoiIiISmYKIiIhEpiAiIiKRKYiIiEhkeX+diMjeYNOmTaxdu5YdO3bkuitSzdWtW5fmzZvTsGHDrLSnICJSzW3atIk1a9bQunVrGjRoQPADgyLluTvFxcV88cUXAFkJJJrOEqnm1q5dS+vWrSkoKFAAkQqZGQUFBbRu3Zq1a9dmpU0FEZFqbseOHTRo0CDX3ZA80qBBg6xNfSqIiOQBHYFIVWTz/4vWRCQrbrzldlaszs19ztq0bMro667JSdsiNV1KQcTMDgaGAycAnYHX3L1rzPauwMwkxWe4+xkV1D0euDDBph+4++JU+ifV34rV6+jS66KctP3hs3/PSbtSPbVv357+/ftzxx135LorNUKqRyKdgZ7Am0DdBNvfJQgwsdoBjwHTUqh/MTA4Lm15in0T2Svl6uhOR3YSK9Ug8qy7Pw1gZpOAprEb3X0TQYApZWYnA7uBx1Oo/1t3f7PybCJSIldHd9k8sisuLtZJBdVcSgvr7r47Qt3nAa+6+8oIZUWkBvrrX/9K27Zt2Xfffenduzcvv/wyZsasWbOAYEH4zjvv5KqrrqJZs2b88Ic/BOD555+ne/fupRfRHX/88cyYMaNM3UVFRTRt2pQ5c+Zw1FFHsc8++3DEEUfw+uuvJ+zLn//8Z9q0aUOTJk0YOHAgGzduzOTQa6yMnJ1lZocARwITUyzSycw2mdk2M3vdzE7JRL9EJHeeeuopLr/8cs4++2yeeuopDjvsMIYOHVou3x//+EdWrVrFP/7xD+666y4Ali1bRq9evfjHP/7B5MmTOfHEE+nRowdz5swpU3bLli2cf/75XHLJJTzxxBM0btyYHj16sHr16jL5Hn/8cV5++WXGjRvHbbfdxnPPPcd1112XucHXYJk6O2sgsAOYnELe+cBbwCKgGfBr4EUz+5G7z8tQ/0Qky2655RZ69uzJ2LFjATj99NNZt24d9957b5l8rVq14rHHHiuTdtlll5X+e/fu3XTr1o2FCxdy//33c9JJJ5VuKy4u5ve//z0/+9nPAOjWrRvt2rVjzJgx/OEPfyjNV7duXaZMmUKdOsFH4KJFi3j00Ue555570jvovUCmrhMZSHBW1leVZXT3v7j7ve7+qrtPAk4FvgCSfi0ws2Fm9o6ZvfPll1+mr9cikhE7d+5k/vz5nH322WXS418D9OzZs1zaihUruPDCC2ndujV16tShbt26zJgxg48//rhc3j59+pT+e7/99qN79+7Mm1f2+2i3bt1KAwhAp06ddG+yiNIeRMzscOAHpD6VVYa7bwGmAkdVkGecuxe6e2GzZs2idVREsmbdunXs2rWL+L/XRH+/LVq0KPN69+7dnH322cydO5ebbrqJmTNn8vbbb9OjRw+2bt1aJu9+++1XbiG+efPmrFq1qkxa48aNy7yuV68e7s62bduqOrS9XiamswYCxcDTe1CHhw8RqQGaNm1K7dq1iZ85SDSTEH+19dKlS5k/fz7Tpk3jzDPPLE0vLi4uV3bz5s3lzuhau3YtrVq12tMhSBKZmM4aSHBK8OYohc2sAXAW8K+09kpEcqZOnToceeSRPP102e+WzzzzTKVlS4JF/fr1S9M+/fTTcovqJZ566qnSf2/evJkXX3yRY489Nkq3JQWpXrFeQHCxIUBroKGZ9Q9fTw2noDCz44H2wNVJ6rkAeAA4yN0/NbNGwHPAP4GlBNefXA0cCJwTZUAiUj2NGDGCfv36cdlll3H22WczZ84cnn/+eQBq1Ur+fbZjx460adOGX//619x8881888033HjjjbRu3bpc3gYNGnD99dezefNmDjzwQO644w62b9/OlVdembFx7e1SPRJpDjwRPo4HOsW8bh6TbyDwNcmvUq8F1AZKjle3AV8CIwnWQcYBG4FT3P2dVAchItVf3759ueuuu5gyZQq9e/fm7bffLr31SEW/e1G/fn2efPJJ6tSpQ//+/bnhhhsYMWIEp5xS/kqAgoICJkyYwD333EO/fv3YsGEDU6dO1XRWBqV0JOLuy/nvB39F+a4Crqpg+3hgfMzrrUDfVPogImW1adk0J/cFa9OyaeWZkrj88su5/PLLS1//7ne/Y5999uHQQw8Fgh9VSuSYY44pd4bVoEGDEuY9+eSTWbBgQdI+LF++vFzaoEGDktYnFdNdfEXyVL7dv+rLL7/k1ltvpVu3bhQUFPDaa69x2223MXToUN3aJI8piIhIVtSrV4/FixczYcIEvv76a1q1asWVV17JzTffnOuuyR5QEBGRrGjUqBFTp07NWP1FRUUUFRVlrH5JTL9sKCIikSmIiIhIZAoiIiISmYKIiIhEpiAiIiKRKYiIiEhkCiIiknGDBg2isLAw193IKyU/91vdKYiIiEhkuthQJE+NuWUUG9d8nvV2G7doy1XX3ZT1djMt/ndIqmLXrl3s2rWLevXqpblX1Z+CiEie2rjmc4p6H5L1doumlP9J2lS9+OKL/PrXv+aTTz7hyCOP5L777qNz586ce+65rF27llmzZpVtq6iIe++9lxUrVvDFF1/QoUMHHn74YaZNm8aUKVNo0KABv/zlL7nxxhvLlPvwww/57W9/y+zZswE488wzufvuu2nZsiUAs2bNolu3bkyfPp2xY8fyyiuvMGDAAO6//37ef/99LrnkEt59910OOugg/vjHPzJq1Ci6dOnC+PHjgWB67sMPP2TkyJFcf/31fPzxx7zyyiscfPDBXH/99cyaNYtVq1bRtm1bzj33XEaNGlUaYJYvX57yOADmz5/PpZdeyvvvv8+hhx7KXXfdxcknnxx5H6SbprNEJCs+++wzhg8fzvXXX8/EiRNZu3YtAwYMwN0ZOnQos2fPZtmyZaX53Z2HHnqI888/n7p165amDx8+nIKCAiZNmsQvfvELRo8ezdixY0u3L126lJNOOomtW7fyz3/+k/Hjx7Nw4UJ69epV7i7BQ4cO5fDDD+eZZ55h6NChbNmyhTPOOIPi4mImTpzIyJEjufrqq/nss8/KjWf58uVcc801jBgxgmnTptGhQwfWrVvH/vvvz5133sn06dMZPnw4Dz74YJk7F6c6DoAtW7Zw4YUXcvHFFzN58mTq169P37592bJlS+T9kG46EhGRrPjqq6+YM2cO3//+94Hgt9P79OnDkiVL6N69O23atGH8+PGMHj0agJkzZ7J8+XIGDx5cpp7OnTtz3333AXDGGWewdu1abrnlFi699FJq1arF6NGjadmyJdOmTSv99n/YYYfRsWNHpk6dyllnnVVa1znnnFPmBpBjx45l/fr1vPPOO6U/enXQQQdx3HHHlRvP+vXreemllzjiiCNK09q0aVP6GykAJ510Evvuuy9Dhgzh7rvvLjPdVdk4IJhiGzNmDD/+8Y8BaNWqFUceeSSzZ88u81PBuaQjERHJivbt25cGEIBOnToBsGLFCmrVqsXgwYOZMGFC6dHC+PHjKSwspEuXLmXq6dOnT5nXffv2ZeXKlaxYsQKAl156iT59+lCrVi127tzJzp076dChA+3bt+edd8r+1l1sQAF4++23Ofroo8v8auKxxx5LixYtyo2ndevWZQIIBEdPY8aMoVOnTjRo0IC6devy85//nG3btpU7mqlsHBDc+bhr164J37PqQkFERLKicePGZV6XfCvfunUrAIMHD+bTTz9l5syZfPPNN0yePJkhQ4aUq6d58+YJX69atQqAdevWcdttt1G3bt0yj//85z98/nnZExHig8Pq1atp1qxZuTYTpSUKLGPGjOE3v/kNffr04emnn2bevHmlU1Ql40x1HADf+c53yvx0cPx7Vh1oOktEqoX27dtz2mmnMX78eJYtW8bu3bs577zzyuVbu3ZtwtclP4G7//7706dPHy666KJyZeOvuzAr+4OtLVu2ZMmSJeXKffnll+XS4ssCPPHEE/Tv35/f//73pWmLFi0qly+VceQLBRERqTaGDh3KkCFDWLhwIb179y539ALw1FNPcemll5a+fvLJJ2nVqhVt2rQB4NRTT2XhwoUcffTRCT/oK3LMMcfwyCOP8MUXX5ROac2bN481a9akVL64uJj69euXSXv44YcT5q1sHPlC01kiUm307t2bffbZh3fffbfcgnqJhQsXcvHFFzNjxgxGjhzJ/fffz4gRI0qnfYqKivjggw8466yzmDRpErNmzeLhhx9m0KBB5U4hjjd48GAOOOAAfvKTnzBlyhQmTpzI+eefT7NmzcpMKyXTvXt3HnvsMe655x5eeOEFLrjgApYuXRppHPkiv3orIjVa/fr16dGjB23btuW0005LmOf2229n06ZN9OvXj/vuu48bbriByy67rHT7IYccwptvvklBQQHDhg2jR48e3HjjjdSvX5+DDz64wvYLCgqYPn06DRo0YMCAARQVFXH77bfTuHFjGjZsWGn/R40axXnnncfIkSM577zzqFevHnfddVekceQLTWeJ5KnGLdru0YV/e9JuVZVcpBerffv25a7b2LlzJzNnzmTIkCFJv5E3atSIiRMnVthex44dmTRpUtLtXbt2Ldd2icMPP5y5c+eWvl62bBnLli3j8MMPL01LNB6A/fbbjwcffLBceqK2KhtHsp/7TdbvXFEQEclTNenWI9u3b+e9997jkUceYf369Vx88cU568utt97KgQceyHe/+10+++wzbr31Vpo1a0a/fv1y1qfqLKUgYmYHA8OBE4DOwGvu3jUuz3Lgu3FF17h7yxTq/ynwO+D7wH+A0e7+WCp9E5H8t3LlSo499liaN2/Offfdl9PFZTNj9OjRrFy5kvr163PyySdzxx13pDSdtTdK9UikM9ATeBOoW0G+R4C7Y15vr6xiM/sRMBm4B7gibGeimW1w9xkp9k9E8liiqa0oedLh2muv5dprr81Y/dkaR7akGkSedfenAcxsEpDsJver3P3NKvbhBmC2u18Rvp5pZp2BUYCCiIhINZbS2VnuvjsTjZtZfaAb8HjcpkeBE8ysUSbaFRGR9Ej3Kb5DzWy7mX1tZpPMLH6NJN5BBNNji+PSPwr7lv37XItUQzVp+kMyL5v/X9J5dtbTBGsmK4AfADcCr5nZD9396yRlmoTPG+PSN8RtF9lr1a1bl+LiYgoKCnLdFckTxcXFZW6fn0lpCyLufmXMy9fMbC6wABgMjElXOwBmNgwYBtCuXbt0Vl3j3XjL7axYvS7r7S744EO69Mp6szVC8+bNS2/D0aBBgyrfykP2Hu5OcXExX3zxRcIbRGZCxq4TcfcPzWwJcFQF2UqOOOLXPprEbY+vexwwDqCwsFDH+VWwYvU6uvQqf2O6TJs775Kst1lTlJxaunLlSnbs2JHj3kh1V7duXVq0aJG1U5IzfbGhh49kPgF2AB2BV2PSOwK7gexfjitSDTVs2FDXKUi1lLF7Z5lZF4Jg8K9kedx9GzATOCdu0wDgjQrWUkREpBpI9Yr1AoKLAAFaAw3NrH/4eirBabrnA88BKwmCx0jgM2B8TD0XAA8AB7n7p2HyzcAsMxsDTAnb6QlUj99+FBGRpFKdzmoOPBGXVvK6A/B5mGcM0BhYD0wHrnP3TTFlagG1gdKVQXd/PQxIvwMuBZYBP9PV6iIi1V9KQcTdlxPzwZ/EqSnUM56YI5OY9CkERyEiIpJH9HsiIiISmYKIiIhEpiAiIiKRKYiIiEhkCiIiIhKZgoiIiESmICIiIpEpiIiISGQKIiIiEpmCiIiIRKYgIiIikSmIiIhIZAoiIiISmYKIiIhEpiAiIiKRKYiIiEhkCiIiIhKZgoiIiESmICIiIpEpiIiISGQKIiIiEpmCiIiIRFYn1x0QEdlTY24ZxcY1n+ek7cYt2nLVdTflpO3qIKUgYmYHA8OBE4DOwGvu3jVmeyvgV8DpwEHABuAVYIS7r6yk7iLgxgSberj79FT6JyJ7t41rPqeo9yE5abtoysc5abe6SPVIpDPQE3gTqJtg+9FAH+DvwFtAC6AImGtmXdx9cyX1fw2cGZf2UYp9ExGRHEk1iDzr7k8DmNkkoGnc9teBju6+syTBzN4FlgD9gIcqqX+nu7+ZYl9ERKSaSCmIuPvuSrZvTJD2sZltAQ6M1jUREanuMnZ2lpkdBhQAqUwYNjazdWa2w8zmm1nfTPVLRETSJyNBxMxqAX8B/g08U0n2pcA1wDkEU18rgckKJCIi1V+mTvG9leBMrlPcfUdFGd39n7GvzexZYC4wCngyURkzGwYMA2jXrl06+isiIhGk/UjEzP6P4HTgC939raqWd3cnCB6HmVntJHnGuXuhuxc2a9ZszzosIiKRpTWImFk/4G7gGnd/bA+q8vAhIiLVWNqCiJl1BR4G7nb3O/agHiNYG3nP3Xelp3ciIpIJqV6xXkBwsSFAa6ChmfUPX08FvgtMARYDj5nZ8THFv3T3T8J6TgFeBk5191fDtFeByWHZfYFfAMcBvSOPSkREsiLVhfXmwBNxaSWvOxB86DcCDidYFI/1EDAo/LcBtcPnEkuBq4BWwG7gXeAsd5+WYt9ERCRHUr3YcDllP/jjjQ8fldUzK74edx+aSh9ERKT60a3gRUQkMgURERGJTEFEREQiUxAREZHIFERERCQyBREREYlMQURERCJTEBERkcgUREREJDIFERERiUxBREREIlMQERGRyBREREQkMgURERGJTEFEREQiUxAREZHIFERERCQyBREREYlMQURERCJTEBERkcgUREREJDIFERERiUxBREREIkspiJjZwWZ2n5m9b2a7zGxWgjxmZteZ2edmVmxms83siBTr/6mZfWBmW81skZkNqNowREQkF1I9EukM9ASWAB8nyXMtcANwG9AL2Ay8ZGYtK6rYzH4ETAZmAj2A54GJZnZ6in0TEZEcSTWIPOvubd39HGBh/EYz24cgiNzq7n9195eAcwAHLquk7huA2e5+hbvPdPfhwHRgVMqjEBGRnEgpiLj77kqynAg0BB6PKfMt8CzB0UVCZlYf6BZbLvQocIKZNUqlfyIikhvpWljvCOwC/h2X/lG4LZmDgLrA4gTlagGHpKl/IiKSAekKIk2Aze6+Ky59A1BgZvUqKAewMUG52O0iIlIN1cl1B6Iws2HAMIB27drluDci1cuYW0axcc3nOWm7cYu2XHXdTTlpW3IjXUFkA7CfmdWOOxppAmxx9+0VlAOIX/toEre9DHcfB4wDKCws9GhdFqmZNq75nKLeuZkJLpqS7ORNqanSNZ21GKgNHByX3pHy6x2xPgF2UH7dpCOwm+SnE4uISDWQriAyF9hEcFovAGZWQHC9yLRkhdx9G8H1IefEbRoAvOHuX6epfyIikgEpTWeFAaFn+LI10NDM+oevp7r7FjP7A3CDmW0gOPr4FUGQujumnguAB4CD3P3TMPlmYJaZjQGmhO30BM7cg3GJiEgWpLom0hx4Ii6t5HUHYDnwB4KgMQI4AHgH6O7ua2LK1CKY9rKSBHd/PQxIvwMuBZYBP3P3GVUaiYiIZF1KQcTdlxPzwZ8kjwO/Dx/J8owHxidIn0JwFCIiInlEd/EVEZHIFERERCQyBREREYlMQURERCJTEBERkcgUREREJDIFERERiSwv7+IrUhXzFyxg6BXXZL3dNi2bMvq67LebS+8tmE/RlYOz3u6iDxZAjm46ubdTEJEar3j7Trr0uijr7X747N+z3mau+fYtObmDcO95c7LepgQ0nSUiIpEpiIiISGQKIiIiEpmCiIiIRKYgIiIikSmIiIhIZAoiIiISmYKIiIhEpiAiIiKRKYiIiEhkCiIiIhKZgoiIiESmICIiIpEpiIiISGRpCyJmNsvMPMnjhCRl2ifJ/2i6+iUiIpmTzt8T+T+gYVzaTcCRwNuVlP0NEPuDAOvS2C8REcmQtAURd18U+9rM6gGFwGPuvrOS4kvc/c109UVERLIjk2siZwJNgIkZbENERHIok0FkILACeC2FvA+a2S4zW2Vmd5pZgwz2S0RE0iQjv7FuZgXA2cB97u4VZN0GjAVmAJuArsBvgYOAn2aibyIikj4ZCSJAL2BfKpnKcvdVwGUxSbPMbA1wj5kd7u7vJSpnZsOAYQDt2rVLT49FRKTKMjWdNRBY6u7vRCg7KXw+OlkGdx/n7oXuXtisWbNIHRQRkT2X9iBiZo2AHkRfUPe4ZxERqaYycSTSB6hP9CDSP3z+V3q6IyIimZKJNZGBwHvu/lH8BjNbCrzq7kPD10XAdwguNNwE/A8wHHjS3d/PQN9ERCSN0nokYmZNgVOBZLctqQPUjnm9GDgFeBCYCvwM+GP4LCIi1Vxaj0TcfR1Qt4Lt7eNeP0rygCMiItVcpk7xFdnrLXtvDkVXDs56u4s+WAC9D8l6u7J3UhARyZDa27+hKAcf5r3nzak8k0ia6PdEREQkMgURERGJTEFEREQiUxAREZHIFERERCQyBREREYlMQURERCJTEBERkcgUREREJDIFERERiUxBREREIlMQERGRyBREREQkMt3FV0RkD7y3YH5ObvnfuEVbrrrupqy3G09BRERkD/j2LTm55X/RlI+z3mYims4SEZHIFERERCQyBREREYlMQURERCJTEBERkcgUREREJLK0BREzG2RmnuBxSSXlGpnZg2a2wcy+NrOHzeyAdPVLREQyJxPXifwYKI55/Z9K8j8OHAJcBOwGbgOmACdnoG8iIpJGmQgib7v75lQymtkJwOnAKe4+O0z7AnjLzE5z95cy0D8REUmTXK+J9ADWlAQQAHefBywLt4mISDWWiSDyiZntNLMlZnZxJXk7AosTpH8UbhMRkWosndNZq4AbgHlAbWAg8P/MrMDd/5ykTBNgY4L0DcD30tg3ERHJgLQFEXd/AXghJmmame0DjDSzv7j77nS1ZWbDgGEA7dq1S1e1IiJSRZleE5kE7A+0T7J9A9AoQXqTcFtC7j7O3QvdvbBZs2Z73EkREYkm00HE457jLSbx2keytRIREalGMh1E+gPrgE+TbJ8GtDSzH5UkmFkhwXrItAz3TURE9lDa1kTMbDLBovr7BAvrA8LHFSXrIWa2FHjV3YcCuPsbZjYDmGBmv+G/Fxu+rmtERESqv3SenbUEGAK0BQxYBFzg7v+Ia692XLkBwJ+BBwiOjJ4Drkhjv0REJEPSeXbWdcB1leRpnyBtIzA4fIiISB7J9RXrIiKSxxREREQkMgURERGJTEFEREQiUxAREZHIFERERCQyBREREYksE79sKJW48ZbbWbF6XU7aXvDBh3TplZOmRaQGUhDJgRWr19Gl10U5aXvuvEty0q6I1EyazhIRkcgUREREJDIFERERiUxBREREIlMQERGRyBREREQkMgURERGJTEFEREQiUxAREZHIFERERCQyBREREYlMQURERCJTEBERkch0F1+p8XatW8bcv43IervffLUm622KZFvagoiZnQP8L3A00AhYAtzh7hMrKecJkt9y9+PT1TfZuzWwbfzqjA5Zb/dX82dlvU2RbEvnkcivgGXA1cA6oCfwiJk1dfe7Kyn7J2BSzOtv0tgvERHJkHQGkV7uHvtzfa+Y2YEEwaWyILLc3d9MY19ERCQL0rawHhdASswHDkxXGyIiUr1k+uysE4CPU8hXZGY7zWydmT1gZvtnuF8iIpIGGTs7y8xOBXoDQyrJ+hDwLPAlUAjcABxuZse6+65M9U9ERPZcRoKImbUHHgGedvfxFeV190ExL2eb2UfAVKAXMCVJ/cOAYQDt2rXb4/6KiEg0aZ/OCqeipgGfAj+PUMV0YDNwVLIM7j7O3QvdvbBZs2bROioiInssrUHEzAqA54B6wE/cfUtV63D3kutGEl0/IiIi1Ug6LzasAzwBfB840d3XRqznTGA/4F/p6puIiGRGOtdE7iG4wPBK4AAzOyBm23x332ZmLwO4+6lQurZRCLxEcIHiUcBIYB7wfBr7JiIiGZDOIHJ6+PyXBNs6AMuB2nHpnwAXAv2AhsBqYAJwg87MEhGp/tIWRNy9fQp5usa9fhl4OV19EBGR7Nqr7+J74y23s2J1ogvtM2vBBx/SpVfWm82pXN1JF8C3bs5Ju9u2bePp6S9mvd31GzZkvU3Ze+3VQWTF6nV06XVR1tudO++SrLeZa7m6ky7A/72Rm5lRd6f5IUnPVM+Ynbvey3qbsvfSj1KJiEhkCiIiIhKZgoiIiESmICIiIpEpiIiISGQKIiIiEpmCiIiIRKYgIiIikSmIiIhIZAoiIiISmYKIiIhEpiAiIiKRKYiIiEhkCiIiIhKZgoiIiESmICIiIpEpiIiISGQKIiIiEpmCiIiIRKYgIiIikSmIiIhIZGkNImbWycxeNrMtZrbSzG4ys9oplGtkZg+a2QYz+9rMHjazA9LZNxERSb866arIzJoALwGLgJ8CBwF/IghUIysp/jhwCHARsBu4DZgCnJyu/omISPqlLYgAlwANgL7uvgl40cwaAkVmdnuYVo6ZnQCcDpzi7rPDtC+At8zsNHd/KY19FBGRNErndFYP4IW4YPEoQWA5pZJya0oCCIC7zwOWhdtERKSaSmcQ6Qgsjk1w98+ALeG2lMuFPqqknIiI5Fg6g0gTYGOC9A3htnSXExGRHDN3T09FZjuA4e4+Ji59BTDB3a9LUu5F4Ft37x2X/k/ge+5+YoIyw4Bh4ctDgSURu90UWBexbHVU08YDNW9MNW08UPPGVNPGA+XHtM7dz0xHxelcWN8ANEqQ3iTcVlG5ZlUp5+7jgHFV7WA8M3vH3Qv3tJ7qoqaNB2remGraeKDmjammjQcyO6Z0TmctJm4Nw8zaAgUkXvNIWi6UbK1ERESqiXQGkWnAGWb2nZi0AUAx8Gol5Vqa2Y9KEsysEPheuE1ERKqpdAaR/wdsA540s9PCdYsi4M7Y037NbKmZ3V/y2t3fAGYAE8ysr5n1Bh4GXs/CNSJ7PCVWzdS08UDNG1NNGw/UvDHVtPFABseUtoV1CG57AvwVOIHgjKu/A0Xuvismz3JglrsPiklrDPwZ6EMQ2J4DrnD3mra4JSJSo6Q1iIiIyN4l7+/iG/Wmj2HZvmb2tpkVm9l6M5tuZvvG5fmpmX1gZlvNbJGZDcjMSMq0mbExmdl4M/MEj4xd2BllPGZWlKSfbmYj4vLmxT5KdUz5so/CcoVmNsPMvgofL5nZcQny5cU+CstVOqY820edw/FsMbN1Znavme2XIF+0feTuefsgOA14JcGNH7sT3L/rW+B3KZS9CNgK3AR0JZhKuxtoFJPnR8BO4C6gG/BHghtEnp7HYxpPcDeA4+Me+1Sn8QBtEvTxD4ADR+TjPqrCmPJlH7UlmLZ+BTgrfMwCNgHfzdN9lOqY8mUfNQLWAK8BvYBBwFpgSly+yPsoIzswWw9gBMG1JA1j0q4huNVKwwrKNQW+AX5RSf0vAK/EpU0lWPTP1zGNB96p7vsoSV3PAx/l6z6qwpjyYh+FH2S7KPslpUmYdmk+7qMqjClf9tEIggDYOCatF8EXl8J07KN8n86KetPHc8Pnh5JlMLP6BBH58bhNjwInmFmiCyvTIWNjypGo4ynDgt+X6Q5MjEnLt31URqIx5UjU8dQl+Pb6bUza5jDNIC/3UaVjypGo4zmCINhtjEl7kSCInAV7vo/yPYhEvenjcQS3ShlqZivMbIeZvWVmsbdYOYjgP1T8BY8fEbxvh+xp55PI5JhKdDKzTWa2zcxeN7OUP/giiDqeeP0I9kfsB26+7aN4icZUIh/20eQwz5/MrLmZNSc4y3ID8ESYJ9/2USpjKpEP+2gfYHtc2k6CqaofhK/3aB/lexCJevPGlgT33BoJ/Jbg8O5bYLqZtYipmwT1b4jbnm6ZHBPAfODX4fafA7UJfvvl2D3ueWLpusHmQOBdd/93XN0kqL+67qN4icYEebKP3H0lwTfYfgTz7muAvsAZ7v5lTN0kqL9a7qMUxwR5so+ApcDhZlY3Ju1ogv7uH1M3CepPaR+l895Z+cSA/YBz3H06gJnNBT4FLgNuyGHfokppTO7+lzKFzKYCC4HrgN5Z7G/KzKwVwSH7b3Pdl3SpaEz5so/CMTwB/IvgpA6AXwLPm9mJ4TflvJLqmPJlHwF/A64E7jazIuAA4B6CNZ7d6Wgg349E9uSmj05w1gUA4Vzjv4BOMXlIUH+TuO3plskxlePuWwgW0I6K0NdURB1PrHMJguRjCeomQf3VdR/FSjamcqrxPhpOMA3S392nh19e+hF8QP0mpm4S1F9d91EqYyqnuu4jd19McMfz84BVwPvAPGABsDqmbhLUn9I+yvcgEvWmjx8R/AHHL5QZ/43OnwA74usPX+8GPo7W5UplckzJePjIhKjjiTWQ4CyRz+PS820fxUo2pmSq4z7qCCx09x0lCe6+neAb+UFhUr7to1TGlEx13Ee4+wNAC+Aw4ECCmYmDgTfDLHu0j/I9iES96eNz4XO3koTwDISjgfcA3H0bMBM4J67sAOANd/96z7qeVMbGlIiZNSA4S+NfUTtciajjAcDM2hOcf19u8TkP9xFQ8ZiS5K+u++hToIuZ1StJCM/06QIsh7zcR5WOKZFqvI8AcPet7v6Bu68Bzif47H883LZn+yhb5zln6NzpJgSHaC8CpxEctm0m7gIcgsWl++PSpoRlLyTY+a8CXwJNYvKUXIAzhuDivdvJzkVSGRkTweHqa8DFwKnhf5I3CW6cWVjdxhOmX0vwLalpkvrzah9VNqZ82kcEX1B2EFzrchbwE4IPux3A4fm4j1IZU57to4bAbeFYziC4uHUHMChdf0cZ2YHZfBDM979CEJFXATcDtePyLAfGx6XtB9wLrA/LvgT8MEH9vYEPw/8gi4GB+TomgtP9ngQ+D8fzNTAdOL46jidMXwBMr6T+vNlHlY0p3/YRwYfobOCr8PEq0DWf91FlY8qnfQTsS3CX9K/Ccm8DvZPUH2kf6QaMIiISWb6viYiISA4piIiISGQKIiIiEpmCiIiIRKYgIiIikSmIiIhIZAoiIiISmYKIiIhEpiAiIiKR/X/DQ9RWA3Lh7QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pts_per_community = 300\n",
    "n = 2*pts_per_community\n",
    "labels = np.hstack([np.full(shape=int(n/2),fill_value = -1),np.full(shape=int(n/2),fill_value = 1)])\n",
    "\n",
    "\n",
    "node_weight_method = 'gaussian_to_central_neighbor'\n",
    "order = 2\n",
    "inner_sidelengths=[1, 3, 1, 1, 1, 1, 1, 1]\n",
    "outer_sidelengths=[2, 4, 2, 2, 2, 2, 2, 2]\n",
    "\n",
    "# experimental procedure parameters\n",
    "k = 5\n",
    "num_trials = 50\n",
    "PPR_iterations = 50\n",
    "\n",
    "for ambient_dim in [7]:\n",
    "    dimension_dependent_data_generation = lambda verbose: generate_concentric_highdim_rectangles(verbose = False,\n",
    "                                                        inner_sidelengths = inner_sidelengths[:ambient_dim],\n",
    "                                                        outer_sidelengths = outer_sidelengths[:ambient_dim])\n",
    "    AUC_vals = []\n",
    "    for trial in range(num_trials):\n",
    "        graph_x, hypergraph_x, _ = compare_estimated_labels('PPR',\n",
    "                                    generate_data = dimension_dependent_data_generation,\n",
    "                                    k = 5,\n",
    "                                    num_iterations = PPR_iterations,\n",
    "                                    diffusion_step_size=None,\n",
    "                                    titlestring=None,\n",
    "                                    node_weight_method=node_weight_method,\n",
    "                                    order = order\n",
    "                                    )\n",
    "\n",
    "        graph_auc_score = metrics.roc_auc_score(labels, graph_x)\n",
    "        hypergraph_auc_score = metrics.roc_auc_score(labels, hypergraph_x)\n",
    "            \n",
    "        AUC_vals.append((hypergraph_auc_score, graph_auc_score))\n",
    "    plot_AUC_hist(AUC_vals,titlestring = f'Ambient dimension = {ambient_dim}') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sorry for the delay. Here's a comprehensive look at how the classification performance varries with different node weights (unweighted, Gaussian L2, Gaussian L1), different distributions (square vs asymmetric boxes), and dimension. All histograms are over 50 independent trials.\n",
    "\n",
    "We compare to the graph PPR vector for the corresponding Gaussian-weighted graph Laplacian.\n",
    "\n",
    "Scales are not consistent across all the histograms, so I advise zooming in to see. I think there are a couple of takeaways:\n",
    "1. Gaussian kernels do seem to increase separation between graph and hypergraph in higher dimensions, compared to unweighted.\n",
    "2. Both graph and hypergraph perform better on the square boxes than the rectangular.\n",
    "3. In high dimensions (dim 7), the Gaussian kernels seem to improve performance of both methods, compared to unweighted. It's difficult to say whether the Gaussian L1 or Gaussian L2 has a stronger effect--both appear comparable. Gaussian L1 seems to really shine in both data distributions, but it's difficult for me to come up with an intuitive reason. It could just be a statistical fluke.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}