elapsed time lag for the same dataset.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# version : Ripser 1.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nml: 58 sec. gap: -.-- sec.\n",
      "fps: 22.9 sec. gap: 35.2 sec.\n",
      "sfl: 3.92 sec. gap: 54.1 sec.\n",
      "cps: 60.5 sec. gap: -2.45 sec.\n",
      "fps) bottleneck distance: 0\n",
      "sfl) bottleneck distance: 0\n",
      "cps) bottleneck distance: 0\n"
     ]
    }
   ],
   "source": [
    "#example 1) spiral\n",
    "import time\n",
    "import numpy as np\n",
    "from sklearn.metrics.pairwise import pairwise_distances\n",
    "\n",
    "import tadasets\n",
    "from ripser import ripser\n",
    "from persim import bottleneck\n",
    "\n",
    "def dpoint2pointcloud(X, i, metric):\n",
    "\n",
    "    ds = pairwise_distances(X, X[i, :][None, :], metric=metric).flatten()\n",
    "    ds[i] = 0\n",
    "    return ds\n",
    "\n",
    "def get_nearby_perm(X, n_perm=None, distance_matrix=False, metric=\"euclidean\"):\n",
    "    \n",
    "    if not n_perm:\n",
    "        n_perm = X.shape[0]\n",
    "    \n",
    "    idx_perm = np.zeros(n_perm, dtype=np.int64)\n",
    "    if distance_matrix:\n",
    "        dpoint2all = lambda i: X[i, :]\n",
    "    else:\n",
    "        dpoint2all = lambda i: dpoint2pointcloud(X, i, metric)\n",
    "    ds = np.zeros((n_perm, n_perm))\n",
    "    ds = dpoint2all(0)\n",
    "    dperm2all = [ds]\n",
    "    \n",
    "    for i in range(1, n_perm):\n",
    "        dsn0 = np.nonzero(ds)\n",
    "        dsn0 = dsn0[0]\n",
    "        idx0 = np.argmin(ds[dsn0])\n",
    "        idx = dsn0[idx0]\n",
    "        idx_perm[i] = idx\n",
    "        \n",
    "        dperm2all.append(dpoint2all(idx))\n",
    "        ds = np.where(ds == 0, 0, dperm2all[-1])\n",
    "        \n",
    "    dperm2all = np.array(dperm2all)\n",
    "    return (idx_perm, dperm2all)\n",
    "\n",
    "N = 2000\n",
    "t = np.linspace(0, 2*np.pi, N+1)[0:N]\n",
    "X = np.zeros((N, 2))\n",
    "X[:, 0] = np.cos(t) * t\n",
    "X[:, 1] = np.sin(t) * t\n",
    "\n",
    "# 1. compute with original dataset\n",
    "tic = time.time()\n",
    "rep = ripser(X)\n",
    "toc = time.time()\n",
    "normal = toc-tic\n",
    "dgms_full = rep['dgms']\n",
    "print(\"nml: %.3g sec. gap: -.-- sec.\"%(normal))\n",
    "\n",
    "# 2. use furthest point sampling\n",
    "tic = time.time() \n",
    "res = ripser(X, n_perm=N)\n",
    "toc = time.time()\n",
    "fps = toc-tic\n",
    "dgms_sub = res['dgms']\n",
    "print(\"fps: %.3g sec. gap: %.3g sec.\"%(fps, normal-fps))\n",
    "# the required time is shorten\n",
    "# can the elapsed time be short if distance matrix is from greedy permutation?\n",
    "\n",
    "# 3. use 'dm' from shuffle point\n",
    "pwd = pairwise_distances(X)\n",
    "a = np.random.permutation(N)\n",
    "dm0 = pwd[:, a]\n",
    "dm = dm0[a, :]\n",
    "tic = time.time() \n",
    "res = ripser(dm, distance_matrix=True)\n",
    "toc = time.time()\n",
    "sfl = toc-tic\n",
    "dgms_sfl = res['dgms']\n",
    "print(\"sfl: %.3g sec. gap: %.3g sec.\"%(sfl, normal-sfl))\n",
    "# random sample take a time shorter than original dataset too (look remark 1 : moving shuffle zone)\n",
    "# then, why the original dataset need so long time?\n",
    "# (the most shortly time is elapsed)\n",
    "\n",
    "# 4. use 'dperm2all' from closest point sampling\n",
    "idx_perm, dperm2all = get_nearby_perm(X)\n",
    "cpsdm = dperm2all[:, idx_perm]\n",
    "tic = time.time() \n",
    "res = ripser(cpsdm, distance_matrix=True)\n",
    "toc = time.time()\n",
    "cps = toc-tic\n",
    "dgms_cps = res['dgms']\n",
    "print(\"cps: %.3g sec. gap: %.3g sec.\"%(cps, normal-cps))\n",
    "# the time lag is disappeared\n",
    "# do it need a long time when points of the dataset that neighboring numbers is nearby?\n",
    "\n",
    "bdist = bottleneck(dgms_full[0], dgms_sub[0])\n",
    "print(\"fps) bottleneck distance: %.3g\"%(bdist))\n",
    "bdist = bottleneck(dgms_full[0], dgms_sfl[0])\n",
    "print(\"sfl) bottleneck distance: %.3g\"%(bdist))\n",
    "bdist = bottleneck(dgms_full[0], dgms_cps[0])\n",
    "print(\"cps) bottleneck distance: %.3g\"%(bdist))\n",
    "# check) has no bottleneck distance\n",
    "# let's check another example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nml: 27.1 sec. gap: -.-- sec.\n",
      "fps: 5.62 sec. gap: 21.5 sec.\n",
      "sfl: 5.36 sec. gap: 21.8 sec.\n",
      "cps: 10.7 sec. gap: 16.4 sec.\n",
      "fps) H0 bottleneck distance: 0\n",
      "fps) H1 bottleneck distance: 0\n",
      "sfl) H0 bottleneck distance: 0\n",
      "sfl) H1 bottleneck distance: 0\n",
      "cps) H0 bottleneck distance: 0\n",
      "cps) H1 bottleneck distance: 0\n"
     ]
    }
   ],
   "source": [
    "#example 2) infty_sign \n",
    "import time\n",
    "import numpy as np\n",
    "from sklearn.metrics.pairwise import pairwise_distances\n",
    "\n",
    "import tadasets\n",
    "from ripser import ripser\n",
    "from persim import bottleneck\n",
    "\n",
    "def dpoint2pointcloud(X, i, metric):\n",
    "\n",
    "    ds = pairwise_distances(X, X[i, :][None, :], metric=metric).flatten()\n",
    "    ds[i] = 0\n",
    "    return ds\n",
    "\n",
    "def get_nearby_perm(X, n_perm=None, distance_matrix=False, metric=\"euclidean\"):\n",
    "    \n",
    "    if not n_perm:\n",
    "        n_perm = X.shape[0]\n",
    "    \n",
    "    idx_perm = np.zeros(n_perm, dtype=np.int64)\n",
    "    \n",
    "    if distance_matrix:\n",
    "        dpoint2all = lambda i: X[i, :]\n",
    "    else:\n",
    "        dpoint2all = lambda i: dpoint2pointcloud(X, i, metric)\n",
    "    ds = np.zeros((n_perm, n_perm))\n",
    "    ds = dpoint2all(0)\n",
    "    dperm2all = [ds]\n",
    "    \n",
    "    for i in range(1, n_perm):\n",
    "        dsn0 = np.nonzero(ds)\n",
    "        dsn0 = dsn0[0]\n",
    "        idx0 = np.argmin(ds[dsn0])\n",
    "        idx = dsn0[idx0]\n",
    "        idx_perm[i] = idx\n",
    "       \n",
    "        dperm2all.append(dpoint2all(idx))\n",
    "        ds = np.where(ds == 0, 0, dperm2all[-1])\n",
    "        \n",
    "    dperm2all = np.array(dperm2all)\n",
    "    return (idx_perm, dperm2all)\n",
    "\n",
    "N = 2000\n",
    "X = tadasets.infty_sign(n=N, noise=0.1)\n",
    "\n",
    "# 1. original datasets\n",
    "tic = time.time()\n",
    "rep = ripser(X)\n",
    "toc = time.time()\n",
    "normal = toc-tic\n",
    "dgms_full = rep['dgms']\n",
    "print(\"nml: %.3g sec. gap: -.-- sec.\"%(normal))\n",
    "\n",
    "# 2. use furthest point sampling\n",
    "tic = time.time() \n",
    "res = ripser(X, n_perm=N)\n",
    "toc = time.time()\n",
    "fps = toc-tic\n",
    "dgms_sub = res['dgms']\n",
    "print(\"fps: %.3g sec. gap: %.3g sec.\"%(fps, normal-fps))\n",
    "# the required time is shorten too\n",
    "\n",
    "# 3. use 'dm' from shuffle point\n",
    "pwd = pairwise_distances(X)\n",
    "a = np.random.permutation(N)\n",
    "dm0 = pwd[:, a]\n",
    "dm = dm0[a, :]\n",
    "tic = time.time() \n",
    "res = ripser(dm, distance_matrix=True)\n",
    "toc = time.time()\n",
    "sfl = toc-tic\n",
    "dgms_sfl = res['dgms']\n",
    "print(\"sfl: %.3g sec. gap: %.3g sec.\"%(sfl, normal-sfl))\n",
    "# the same result is happened as two unconnected circle\n",
    "# as I had expected, if nearby, the time is long?\n",
    "\n",
    "# 4. use 'dperm2all' from closest point sampling\n",
    "idx_perm, dperm2all = get_nearby_perm(X)\n",
    "cpsdm = dperm2all[:, idx_perm]\n",
    "tic = time.time() \n",
    "res = ripser(cpsdm, distance_matrix=True)\n",
    "toc = time.time()\n",
    "cps = toc-tic\n",
    "dgms_cps = res['dgms']\n",
    "print(\"cps: %.3g sec. gap: %.3g sec.\"%(cps, normal-cps))\n",
    "# the elapsed time is longer than shuffle and furthest point sampling\n",
    "# but it is half time of computing original dataset.\n",
    "# probably, there is other reasons...\n",
    "\n",
    "for dim, (I1, I2) in enumerate(zip(dgms_full, dgms_sub)):\n",
    "    bdist = bottleneck(I1, I2)\n",
    "    print(\"fps) H%i bottleneck distance: %.3g\"%(dim, bdist))\n",
    "for dim, (I1, I2) in enumerate(zip(dgms_full, dgms_sfl)):\n",
    "    bdist = bottleneck(I1, I2)\n",
    "    print(\"sfl) H%i bottleneck distance: %.3g\"%(dim, bdist))\n",
    "for dim, (I1, I2) in enumerate(zip(dgms_full, dgms_cps)):\n",
    "    bdist = bottleneck(I1, I2)\n",
    "    print(\"cps) H%i bottleneck distance: %.3g\"%(dim, bdist))\n",
    "# check) has no bottleneck distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nml: 28.9 sec. gap: -.-- sec.\n",
      "sfl: 27.2 sec. gap: 1.67 sec. (shuffle : 500 points)\n",
      "sfl: 27.3 sec. gap: 1.59 sec. (shuffle : 1000 points)\n",
      "sfl: 18.7 sec. gap: 10.2 sec. (shuffle : 1500 points)\n",
      "sfl: 5.41 sec. gap: 23.4 sec. (shuffle : 2000 points)\n"
     ]
    }
   ],
   "source": [
    "# remark 1 : moving shuffle zone (use : infty_sign, without check bottleneck distance)\n",
    "import time\n",
    "import numpy as np\n",
    "from sklearn.metrics.pairwise import pairwise_distances\n",
    "\n",
    "import tadasets\n",
    "from ripser import ripser\n",
    "from persim import bottleneck\n",
    "\n",
    "def dpoint2pointcloud(X, i, metric):\n",
    "\n",
    "    ds = pairwise_distances(X, X[i, :][None, :], metric=metric).flatten()\n",
    "    ds[i] = 0\n",
    "    return ds\n",
    "\n",
    "def get_nearby_perm(X, n_perm=None, distance_matrix=False, metric=\"euclidean\"):\n",
    "    \n",
    "    if not n_perm:\n",
    "        n_perm = X.shape[0]\n",
    "    \n",
    "    idx_perm = np.zeros(n_perm, dtype=np.int64)\n",
    "    \n",
    "    if distance_matrix:\n",
    "        dpoint2all = lambda i: X[i, :]\n",
    "    else:\n",
    "        dpoint2all = lambda i: dpoint2pointcloud(X, i, metric)\n",
    "    ds = np.zeros((n_perm, n_perm))\n",
    "    ds = dpoint2all(0)\n",
    "    dperm2all = [ds]\n",
    "    \n",
    "    for i in range(1, n_perm):\n",
    "        dsn0 = np.nonzero(ds)\n",
    "        dsn0 = dsn0[0]\n",
    "        idx0 = np.argmin(ds[dsn0])\n",
    "        idx = dsn0[idx0]\n",
    "        idx_perm[i] = idx\n",
    "       \n",
    "        dperm2all.append(dpoint2all(idx))\n",
    "        ds = np.where(ds == 0, 0, dperm2all[-1])\n",
    "        \n",
    "    dperm2all = np.array(dperm2all)\n",
    "    return (idx_perm, dperm2all)\n",
    "\n",
    "N = 2000\n",
    "X = tadasets.infty_sign(n=N, noise=0.1)\n",
    "\n",
    "# 1. original datasets\n",
    "tic = time.time()\n",
    "rep = ripser(X)\n",
    "toc = time.time()\n",
    "normal = toc-tic\n",
    "dgms_full = rep['dgms']\n",
    "print(\"nml: %.3g sec. gap: -.-- sec.\"%(normal))\n",
    "\n",
    "# 2. use 'dm' from shuffle point (moving shuffle range)\n",
    "for i in range(500, N+1, 500):\n",
    "    pwd = pairwise_distances(X)\n",
    "    a = np.random.permutation(i)\n",
    "    dm0 = np.zeros((N, N))\n",
    "    dm0 = pwd[0:i, a]\n",
    "    dm[0:i, 0:i] = dm0[a, 0:i]\n",
    "    dm[0:i, i:N] = pwd[a, i:N]\n",
    "    dm[i:N, 0:i] = pwd[i:N, a]\n",
    "    dm[i:N, i:N] = pwd[i:N, i:N]\n",
    "    tic = time.time() \n",
    "    res = ripser(dm, distance_matrix=True)\n",
    "    toc = time.time()\n",
    "    sfl = toc-tic\n",
    "    dgms_sfl = res['dgms']\n",
    "    print(\"sfl: %.3g sec. gap: %.3g sec. (shuffle : %i points)\"%(sfl, normal-sfl, i))\n",
    "# the bigger shuffle range spread, the shorter time elapsed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "got perm = \n",
      "[ 0 29 28 27 26 25 24 23  7  6  5  2  1  4  3 21 22  8  9 11 10 12 13 14\n",
      " 15 16 18 19 20 17]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAHgCAYAAABej+9AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeVxVZeLH8c8BlcXrjqZiirtSiSQaamaalToq6aghLWPZlJOVWmn9Wsxf02RTk2tm06It/lK7ZIoKbZg5uVRYtrgkFppa5oKIioDA8/sDYUBRUbn33Hv5vl8vXnHPfe69XxqmL+ec5zzHMsYgIiIi3snP7gAiIiJy4VTkIiIiXkxFLiIi4sVU5CIiIl5MRS4iIuLFVOQiIiJerIqdH25Z1lxgALDPGHN5Gc9fCywF0k5uWmyMefps7xkSEmLCwsIqOKmIiIjrHTp0iF9++QXggDGmfnleY2uRA28CLwFvn2XMf4wxA8r7hmFhYaSkpFxsLhEREbfKzMykZcuWdO/enTVr1uws7+tsLXJjzGrLssLszCAiIuIJatasSXJyMs2bN6dmzZrlfp03nCPvalnWd5ZlJVmWdZndYURERCqS0+nkn//8JwAdOnSgRo0a5/V6Ty/yb4BmxpgIYBawpKxBlmXdbVlWimVZKfv373drQBERkQvldDoZMWIEy5YtIzc394Lew6OL3BiTaYw5evL7RKCqZVkhZYx71RgTZYyJql+/XHMDREREbFVU4tHR0SQlJVGtWrULeh+PLnLLshpalmWd/L4LhXkP2ptKRETk4pxa4ud7OL0kuy8/WwBcC4RYlrUbeAqoCmCMeQUYCvzNsqw84DgQa3S7NhER8XKHDx+mW7durFix4qJKHMDytV6MiooyuvxMREQ80cGDB6lXrx4A+fn5+Pv7lznOsqwNxpio8rynRx9aFxER8RVOp5PmzZuzfv16gDOW+PlSkYuIiLhY0TnxiIgILrusYq+kVpGLiIi4UFGJd+3alcTExIs+J34qFbmIiIiLfPnlly4tcVCRi4iIuEznzp159tlnXVbioCIXERGpcAkJCfz666/4+fkxceJEl5U4qMhFREQqlNPpZMiQITz22GNu+TwVuYiISAUpObFtzpw5bvlMFbmIiEgFcPXs9DNRkYuIiFyk/Px8/vWvf7m9xMHmtdZFRES8nTEGf39/kpKSqFq1qltLHLRHLiIicsGcTicxMTFkZ2dTt25dt5c4qMhFREQuSNE58UOHDpGXl2dbDhW5iIjIeTp1YpvD4bAti4pcRETkPCxevNiW2elnoiIXERE5Dy1atOBPf/qTR5Q4qMhFRETKZfPmzQB07NiRpUuXekSJg4pcRETknOLj4+nQoQPz5s2zO8ppVOQiIiJnER8fT2xsLF27dmXo0KF2xzmNilxEROQMSpa4p5wTP5WKXEREpAy7du3illtu8egSBy3RKiIiUqZLL72UxYsXc80113hsiYP2yEVEREpxOp18+OGHAPzpT3/y6BIHFbmIiEixohXbXnzxRYwxdscpFxW5iIgIpZddXbx4MZZl2R2pXFTkIiJS6Z26drqnH04vSUUuIiKV3ueff+6VJQ6atS4iIpVYTk4OAQEBzJw5k+zsbIKDg+2OdN60Ry4iIpWS0+kkPDycnTt34ufn55UlDipyERGphIrOiTdu3Ji6devaHeeiqMhFRKRS8eaJbWVRkYuISKXx8ccf+1SJg4pcREQqkauuuoq//e1vPlPioCIXEZFKYOXKlWRlZVGrVi1mzZrlMyUOKnIREfFxTqeTG264gcmTJ9sdxSVU5CIi4rNKTmx78skn7Y7jEipyERHxSb42O/1MVOQiIuJzsrKyGDt2rM+XOGiJVhER8UHBwcF89tlnNG7c2KdLHLRHLiIiPsTpdPLEE09gjKFt27Y+X+KgIhcRER9RdE581apV5OTk2B3HbVTkIiLi9YpKPDo6mqSkJAIDA+2O5DYqchER8WqnlnhlOJxekopcRES8Wn5+PldffXWlLHFQkYuIiJf6448/AIiNjWXlypWVssRBRS4iIl7I6XTSvHlzPvvsMwD8/CpvnVXen1xERLxS0TnxK6+8kqioKLvj2E5FLiIiXqOyT2wri4pcRES8wjfffKMSL4OKXEREvEJkZCRTp05ViZ9CRS4iIh5tyZIlpKamYlkWDzzwgEr8FCpyERHxWE6nk6FDh/LEE0/YHcVjqchFRMQjlbyf+Ouvv253HI+lIhcREY9TssR9/X7iF0tFLiIiHsUYw8svv6wSL6cqdgcQEREpYozBsiwSEhIAVOLloD1yERHxCE6nk+uvv55jx45Ro0YNlXg5qchFRFwsPj6ebt26Ua9ePQIDA2nbti3PPPMMubm5dkfzGEXnxHNycigoKLA7jlfRoXURERc7ePAgvXr1YsKECdSuXZuvvvqKyZMns3fvXl566SW749lOE9sujmWMsTtDhYqKijIpKSl2xxAROaM9e/bQvHlzTpw4QWZmZqUurg8++IBhw4apxE9hWdYGY0y57gijQ+siIm42YcIEAgMD7Y7hEdq3b8/gwYNV4hdBRS4i4ib5+fl88sknLF++HH9/fwAsy7I5lT02btyIMYZ27drhdDpV4hdBRS4i4ibBwcHccMMNHDlyhPDwcLvj2MbpdBIVFcWcOXPsjuITVOQiIm4yfvx4mjZtyvPPP8+3335rdxxblJzYdtttt9kdxyeoyEVE3ODgwYO8+uqrzJkzhwkTJnDLLbcA8Msvv9iczH00O901VOQiIm7w+OOPc9VVV9G/f38AmjVrBsDOnTvtjOU2e/fu5fbbb1eJu4CuIxcRcYE/MrP5IvUAjWsHUTN7L3PnzmX16tVkZGQAsGnTJgDq1q3L8ePHCQoKsjOuyzVs2JDly5fTpUsXlXgFs/U6csuy5gIDgH3GmMvLeN4CZgD9gSxgpDHmm7O9p64jFxG77T6UxU2z13DgaOHKbfnvT2D39i1nHD9q1CifvU1nfHw8AEOHDrU5iXfxpuvI3wT6nuX5fkDrk193A5riKCIeb9l3vxeXOEBO/faEhYURGBhI9erVadWqFX369AEgMTGRCRMm2BXVpeLj44mNjWX27Nn42uJjnsTWIjfGrAbSzzIkBnjbFFoP1LYsq5F70omIXJhaQVVLPXb0vJMP133P8ePHOXr0KKmpqcWT3Xr06EHbtm3tiOlSRSXetWtXEhISKu318u5g9x75uYQCu0o83n1yWymWZd1tWVaKZVkp+/fvd1s4EZGyDItqQvVq/sWPCwx88O0eGxO5V8kS18Q21/P0Ii/rT7jTjs8YY141xkQZY6Lq16/vhlgiImdW1d+P9o1qltpWJ7j0XvrIkSMxxuBwONwZzS1SUlJU4m7k6UW+G7i0xOMmwG82ZRERKbdnh1xB07rBWBZc164Bf+kWZncklzt+/DgAU6ZM4ZNPPlGJu4mnF3kCcLtVKBo4bIz53e5QIiLn0uaSGqye2IvUZ/rxxsjOBFb1P/eLvFh8fDytW7cmNTUVy7J0Uxg3svU6csuyFgDXAiGWZe0GngKqAhhjXgESKbz0bDuFl5/dYU9SEZELU8Xf0/eXLl7Jc+INGza0O06lY2uRG2NGnON5A4xxUxwRETlPmthmP9//U9GHbN++nXvuuYeIiAj8/f259tprSz2/atUqLMsq8+vGG2+0J7SI+Kzk5GSVuAfQEq1eZNOmTSQmJhIdHU1ubu5pz1955ZWsW7eu1LZff/2Vm2++mX79+rkrpohUEtHR0YwbN46nnnpKJW4j7ZF7kYEDB7Jr1y6cTieXXXbZac/XrFmT6OjoUl87d+7Ez8+P4cOH25C4YjidTgYNGkRoaCgOh4NOnTqxYMGC08a99tprtG7dmsDAQDp16kRycrINaUV83yeffEJmZibVq1fnX//6l0rcZipyL+Lnd/7/cy1YsICePXvSuHFjFyRyj6lTp+JwOJg2bRoJCQn06tWLuLg4Zs2aVTxm4cKFjB49mttvv52kpCQuu+wyBgwYwI8//mhjchHf43Q66devH08++aTdUeQkW2+a4gqV5aYpQ4cO5cCBA6xateqMY1JTU2nTpg2vvvoqf/3rX90XroIdOHCAkJCQUtvi4uJYt24daWlpALRt25bu3bszd+5cAAoKCoiIiCAiIoL58+e7PbOIL9L9xN3Hm26aIi60YMECqlatyp///Ge7o1yUU0scIDIykn379gHwyy+/sG3btlKnD/z8/Bg2bBhJSUluyyniy1TinktF7sMWLlzIDTfcQN26de2OUuHWrl1LeHg4AFu3bgWgXbt2pca0b9+e9PR0tP6+yMXJzs5m4sSJKnEPpSL3Ud999x1btmxhxIizXqpvi3NdRgeQkZHBnXfeSd26dXE4HPTr14/t27cDhZe8LF26lDFjCpcYOHToEAC1a9cu9R516tQp9byIXJjAwEBWrlypEvdQKnIftXDhQoKCgoiJibE7ymmKLqNr06YNbdq0KXPMzTffzEcffcSMGTN49913OXjwINdddx0//PADcXFxxMTEMHLkyFKvOfU2iUXzP3T7RJEL43Q6GTduHMYYmjdvrhL3ULqO3IMVFBg+2rSXY7n53HDZJdQMrHruF520aNEiBg4c6JF3Vho4cGDxHxhFk/ZKWrduHR9//DHJycn07t0bgKuuuoqwsDD69OlD06ZNS01gK9rzzsjIoFatWsXbMzIygNP31EXk3EqeE8/OziYoKMjuSHIGKnIP9sDCb1n+feE9Ylp97uDdOyJZnfwxAHv27CEzM5P4+HgA+vfvT3BwMADr168nLS2NqVOn2hP8HM51Gd3GjRupUqUKPXv2LN5Wo0YNqlSpQmZmJitWrKB69erFzxWdG9+6dSvNmjUr3r5161bq1q2Lbm0rcn5OndimEvdsKnIPdehYbnGJA2zfd5SPUn7i9mHDSo0bdvJxWloaYWFhQOFh9Vq1anntam7Z2dlYfv78+ZX1hNYJ4ol+bbnr1uHk5ORQs2ZNGjRoUGp8ixYtaNOmDU6ns3gp2oKCguLrXUWk/OLj4zU73cuoyD1UUDV/Aqv6kX2ioHhbeJuWlOe6/+nTpzN9+nRXxnOpI9XqcSI3h6++2Ui1+mF8/MrTbP9PIsHBwWRkZLB+/frisZGRkQQEBDB58mRuvfVWwsLC6N69O2+99Rapqam8++67Nv4kIt4nMDCQnj17smTJEpW4l9BkNw8VWNWf6Td3pEZgFar6W9zfuxWdmtlzGVl5Zpn//vvv3HHHHcXLqEZGRvJ///d/F/R5jpadqFLrEtI/fIkTB3ez87s1AGRlZZGfn0/Xrl2Lv37/vfCoxYgRI3jllVd488036du3L99//z3Lly/n8ssvv+CfW6Qy2bNnDwADBgzg008/VYl7Ee2Re7C+lzei7+WNyC8w+PvZN/P6XDdrKSgoYNCgQRw8eJDnn3+ehg0bEh8fz6233kpwcDCDBw8+r8+7pl0jGsRMZF/CC/z2+mgArr76alq3bs3KlSvZsWNHma/761//6tUr2InYxel0cuutt7J06VL69u2rKz28jIrcC9hZ4nDuWebbtm0jJSWFhIQEBg4cCMB1113Hl19+ycKFCxk8eDC70rPIKzA0D6l+2vufKiqsLgsev40VA6+jytG93NK1OeFt2zBgwACio6Mr/gcUqcRKTmzr3r273XHkAqjI5ZzONcv8xIkTAKUu/YLCy76MMbz48U/MWlm4mMvQTk3417CIc37mNW3qc02b+kDh2NTUVD799FOWLVt2AT+BiJRFy676BhW5XLTLL7+cq666ikmTJvHaa69xySWXsHjxYtasWcP/OZfy8MkSB4jfsJshHUL47YfC+6af6TK6v//977Rr146QkBB++OEH/v73vxMbG8v1119vy88o4mt+/PFHlbiPUJHLRbMsi6SkJGJiYopXaqtatSrz5s2j2zU9Ye3KUuMP7t9ffNlckVMvozt48CDjxo3jwIEDXHrppTz88MM89NBD7vmBRCqByy67jFdeeYWbb75ZJe7lVORy0QoKCrjttts4ePAgixYtokGDBiQmJjJq1CiWLKnHHd3DmLdmBwD9r2jIgO4R57yMztsvoRPxVIsXL6ZVq1Z06NCBu+66y+44UgFU5HLRli9fzooVK9i2bRutW7cG4Nprr2XXrl1MnDiR77//nps7X8qJPMPloTU1I1bEJkXnxAcOHMgHH3xgdxypILqOXC7a1q1bCQ4OLi7xIpGRkfz8888AtGtYkyua1FKJi9ikqMSjo6N5++237Y4jFUhFLqcxxrDwq1/532WbWLn1j3OOb9asGVlZWfz000+ltm/YsKF42VgRsU/JEk9KStI5cR+jQ+tympnJ25n26TYA5q3Zwcyh7cnd8Q1Q9izz/v3707RpU2666SYmTZpE/fr1WbFiBe+99x6zZ8+27ecQkcI/zN98880KK3Gn08k777zDhg0bOHz4MG3btuXhhx9mxIgRFZRYzpdVnrW7vUlUVJRJSUmxO4ZXG/TSF3y/+3Dx435hfrzyt/5lji2aZb59+3b+53/+hzVr1pCZmUnLli259957ufvuu3U4XcQmBQUF+Pn5cfz4cfLy8ipkT7xr1640b96cm266iZCQEBITE3nxxReZOXMm999/fwWkFgDLsjYYY6LKM1Z75HKalvUdpYr8yvC255xl3qpVK5xOp6ujiUg5OZ1Opk2bRmJiIrVr166w9122bBkhISHFj3v37s1vv/3G1KlTVeQ20TlyOc3kgZfxpw6NaNXAwR3dw7jz6uZ2RxKR81B0TtzPzw9/f/8Kfe+SJV4kMjKSffv2VejnSPlpj1xOUyu4KrPjrrQ7hohcADsmtq1du5bw8HCXf46UTUUuIuIjli5d6vYST05OZunSpcydO9flnyVl06F1EREfERERwfDhw91W4jt27CAuLo6YmBhGjhzp8s+TsqnIRUS83Ndff01BQQFhYWG8++67binx9PR0+vXrR9OmTZk/f77LP0/OTEUuIuLFnE4nXbt2Zdq0aS79nOwT+fyRmU1BgSErK4sBAwaQm5vLihUrqF69uks/W85O58hFRLxUyfuJ33333S77nA070xn1VgoZWSfo0NhBTtI/SU1NZc2aNTRo0MBlnyvloyIXEfFCJUvc1fcTn5ywmYysEwB8NncKR7/7iBkzZpCens769euLx0VGRhIQEOCyHFI2FbmIiJc5cOAAd955p1tKHOD4ifz/fp/2LQBjx449bVzRSo/iXjpHLiLiZUJCQvjwww/dUuIA9/Vqhb9f4VLLV078P3alH8MYc9qXStwe2iMXEfESTqeTY8eOMXLkSLp37+62z70pMpTLQ2vya3oWHS+tQ93q1dz22XJu2iMXEfECRefE582bR0FBgds/v1WDGvRud4lK3AOpyEVEPFzJiW3Lly/Hz0//6Zb/0m+DiIgHc+fsdPFOKnIREQ+2bds2lbiclXWu+0x7m6ioKJOSkmJ3DBGRi3L06FEcDgcAubm5VKumc9OViWVZG4wxUeUZqz1yEREP43Q6adGiBT/++COASlzOSkUuIuJBis6Jt23blmbNmtkdR7yAilxExENoYptcCBW5iIgH+Pzzz1XickFU5CIiHiA6OppHHnlEJS7nTUUuImKjDz/8kIMHDxIQEMA//vEPlbicNxW5iBs4nU4GDRpEaGgoDoeDTp06sWDBglJjFi1axJAhQ2jUqBGWZfHmm2/aE1bcJj4+ngEDBvD444/bHUW8mIpcxA2mTp2Kw+Fg2rRpJCQk0KtXL+Li4pg1a1bxmPj4eHbs2MGAAQNsTCruEh8fT2xsLF27duWFF16wO454MS0II+IGBw4cICQkpNS2uLg41q1bR1paGgAFBQX4+flx9OhRatSowbx58xg5cqQNacXVSpa4zolLWbQgjIiHObXEASIjI9m3b1/xY90Io3LIzc3liSeeUIlLhdH9yEVssnbtWsLDw+2OIW5WrVo1kpOTqVmzpkpcKoR2AURskJyczNKlSxkzZozdUcRNnE4nd911F/n5+YSGhqrEpcKoyEXcbMeOHcTFxRETE6Nz4JVE0YptP/30E9nZ2XbHER+jIhdxo/T0dPr160fTpk2ZP3++3XGkAr355ptYlnXa11//+tdSK7ZVr17d7qjiY3SOXMRF8gsMezOzqRtcjaBq/mRlZTFgwAByc3NZsWKF/oPuo1auXElQUFDx908++STdunXTxDZxGRW5iAtkZp/gtte/5Lvdh6kZWIVXbonkH+PuIDU1lTVr1tCgQQO7I4qLdO7cufg+4llZWfTp04f4+HiVuLiMilzEBd5eu4Pvdh8GIDM7j9i/3MWudYnMmDGD9PR01q9fXzw2MjKSgIAANm/ezObNm4vPoaakpOBwOKhfvz49e/a05eeQC7Njxw7CwsLo3bs3vXr1wrIsuyOJD1ORi7hATl5Bqcf7t34FwNixY08bm5aWRlhYGO+99x7/+7//W7x99uzZzJ49m549e7Jq1SqX5pWK06RJEw4fPkzjxo2ZNGkS99xzj92RxMdpspuIC4zo0pSGNQMB8PezWLhyA8aYMr/CwsIAmDx5cpnPq8S9Q6NGjbj55ps5cuQI4eHhXHPNNYwePZpp06bZHU18nPbIRVygce0gPhzXg+92H6ZJnSBa1nfYHUlcLDMzk/j4+FIT2woKCnjmmWcYO3asVu4Tl9FvloiL1A6uRs829VXilcC2bdtKXWJWNLFt6NChpKens2PHDnsDik9TkYuIXKQ2bdrw1ltvnfESM012E1fSoXURkQuQm1fAvc+8zPGA2vxlUB9uueWW08a8//77hISE0KxZMxsSSmWhIhcRuQA3PzqVJVMfIahFJ9Zm1qHpN3Po37sHHTp0ID8/n0WLFrFo0SJmzpyp8+PiUrYWuWVZfYEZgD/wujHmuVOeHwm8AOw5ueklY8zrbg0pInIKp9PJkmmPEBDajpCBEzAGrNqNmTt3Lrt27cIYQ3h4OG+//Ta33Xab3XHFx9lW5JZl+QOzgeuB3cDXlmUlGGM2nzJ0kTHmPrcHFBEpQ9ENUBq26kDV/o/hFxAMwPj/eYqYjqE2p5PKyM7jPV2A7caYX4wxucBCIMbGPCIiZ2WMwel00rVrV75a/SlDo1tzZdPaPNqvnUpcbGPnofVQYFeJx7uBq8oY92fLsq4BtgHjjTG7yhgjIuJSeXl5VKlShfnz55Obm4vD4WDqzfXsjiVi6x55WddjmFMeLwPCjDEdgE+Bt8p8I8u627KsFMuyUvbv31/BMUWksnM6nXTu3Jn9+/dTrVq14puiiHgCO4t8N3BpicdNgN9KDjDGHDTG5Jx8+BrQqaw3Msa8aoyJMsZE1a9f3yVhRaRyKjon7nA4CAwMtDuOyGnsLPKvgdaWZTW3LKsaEAsklBxgWVajEg8HAVvcmE9EKrmiEj91xTYRT2JbkRtj8oD7gI8oLOj3jDGbLMt62rKsQSeHPWBZ1ibLsr4DHgBG2pPWu23fvp177rmHiIgI/P39ufbaa8sc98MPPzBgwABq1apFjRo16NKlCxs2bHBvWBE3ysvL47nnnqN169YEBATQpEkTxo8fD8Dy5ctV4uIVbL2O3BiTCCSesm1Sie//B/gfd+fyNZs2bSIxMZHo6Ghyc3PLHLNx40Z69OhBTEwMixYtAuDrr7/m+PHj7owq4lZ33HEHycnJPPXUU7Rr145du3axeXPhFbCdOnXitttuY+bMmSpxN9u+fTsvvPAC69ev58cff6RHjx66C+BZaGW3SmDgwIHExBRe2Td06FAOHDhw2pjRo0czcOBA5s+fX7ytb9++bsso4m4ffvghCxcu5LvvviM8PLx4+9q1a8nLy6NRo0bMmzfPxoSVV3l2PuS/tG5gJXCu5SE3b97Ml19+yf333++mRCL2mzt3Lr179y5V4k6nk2uuuYbnnnvuLK8UVxs4cCC7du3C6XRy2WWX2R3H46nIhS+//BKAQ4cOERERQZUqVWjZsiVvvPGGzclEXOfLL7+kTZs23HfffdSsWZOAgACGDx9OZGQkY8eOtTtepaa16c+P/m0Je/fuBeD222/nlltu4ZNPPqFv377cddddJCYmnuPVIt5p7969vPnmm2zcuJExY8Zw4sQJAgICMMboOnHxKipyoaCgAIC77rqLiRMn0qtXL2bPnk2vXr2YMmWKzelEXMMYgzGGt99+m1deeYVu3bqxZMkSNmzYwMqVK+2OJ1JuKnKhbt26APTq1avU9t69exfP4BXxNXXq1OGKK66gRYsWfPLJJyQlJXHDDTdQrVo1/d6LV9GsdR+Vl1/AviM5hDgCqFbl7H+vtW/fvsztxhidqxKfFRISUnxaKSoqCig8OqXfe/E2+m31QXsyjtNn6ud0e24l177wGT/vP3rW8d26daNOnTokJyeX2p6cnExERIQro4q4VUGB4dXVP/Oncc+zZcsWfv311+IyB1i9ejUnTpzQ7714Fe2R+6CZn6ay42AWAL8dzmZKwnf0rfU7AHv27CEzM5P4+HgA+vfvT3BwMJMmTWLixInUrl2bzp078/7777N69Wo+//xz234OkYo2/dNtTHn5TQ4kPE+1Rm1w5B1m8ODBPPbYYxw5coRHHnmEPn36cPXVV9sdtVI5npvPs4lb2Px7Jt1bhTDuutb4+ZV1Xy0pi4rcB2Xn5Zd6nJG+n2F3DSu1bdiwwsdpaWmEhYUxbtw4CgoKmDVrFpMnT6Zt27bEx8fTo0cPt+UWcbX3nE4OJDxPQGg7GgydzHVhAfzx0SvExsZSrVo1YmJimDZtmt0xK51/JG5m/vpfAdiw8xCBVh4NMgrnKZxp50P+S0Xug0Zd3ZyVW/ZxJCePoKr+PDQkmvgJp94h9nQPPvggDz74oBsSitijNseKS9wvIJjuV7bnrrG6xNJum37LLPX4m607eOP+s+98yH+pyH1Qhya1+eTBnmz5PZPWlzhoUkd/vUrllpGRQe3atVny8jM8uzyWLX8co2vLEO7s3tzuaAJ0bxnCt79mFD/u17UDr5tz73xIIRW5j2pYK5CGtXTvZBGn08ndd99NcnIyV155JU8P1kQ2TzP++jbUCqrKlt8z6dYqhD93amJ3JK+iIhcRn1Xyfg6kn8AAACAASURBVOKtW7e2O46cgb+fxV+vaWF3DK+ly89EKpDT6WTQoEGEhobicDjo1KkTCxYsKH4+MzOTp556ii5dulCrVi0aNmzI4MGD2bZtm42pfVPJEtf9xMWXqchFKtDUqVNxOBxMmzaNhIQEevXqRVxcHLNmzQLg119/5bXXXuPGG28kPj6ef//73/z+++9cddVV7Nq1y+b0vmPNmjUqcak0LONjEwqioqJMSkqK3TGkkjpw4AAhISGltsXFxbFu3TrS0tI4duwYfn5+BAUFFT+fnp5O06ZNmTBhAk899ZS7I/ukvLw8pkyZwrhx41Ti4pUsy9pgjIkqz1jtkYtUoFNLHCAyMpJ9+/YBUL169VIlDoVr3Tdr1qx4jFy4FStW8Pvvv1OlShWefPJJlbhUCipyERdbu3Yt4eHhZ3x+//79bN++/axj5NycTicxMTE8/vjjdkcRcSvNWhdxoeTkZJYuXcrcuXPPOOahhx7C4XAQGxvrxmS+peTEthkzZtgdR8StVOQiLrJjxw7i4uKIiYlh5MiRZY6ZM2cO8+fP5/3336devXruDegjNDtdKjtNdhNxgfT0dLp3747D4WDVqlVUr179tDEJCQkMGTKEKVOmMGHCBBtSer+8vDyioqKoUaOGSlx8yvlMdtMeuUgFOJ5beKOaoGr+ZGVlMWDAAHJzc1mxYkWZJb527VpiY2MZPXq0SvwiVKlShY8//pigoCCVuFRaKnKRi/TK5z/z/IdbMcCDfVry4dSHSE1NZc2aNTRo0OC08Zs2bWLAgAH07duXmTNnuj+wD3A6nbz//vu88847Zf47FqlMVOQiF2FXehb//HArRWeoHntoHEe/+4gZM2aQnp7O+vXri8dGRkZy+PBh+vbti8Ph4IEHHuCrr74qfr5mzZqauV4OJc+J5+TkULVqVbsjidhKRS5yEY7m5FFymsnxtG8BGDt27Glj09LS2LFjB7t37wagV69epZ7v2bMnq1atcllWX3DqxDaHw2F3JBHbqchFLkLbS2rQu10DVm4tXMzl1unLmPuXzvj5WWWODwsLw9cmmLrL+++/r9npImVQkYtcBD8/i9duj2L1tv0YDNe0rn/GEpeLExoayo033sjChQtV4iIl6PIzEfFoqampugWpVDpaa11EfEJ8fDzh4eG8++67dkcR8VgqchHxSPHx8cTGxhIdHc3AgQPtjiPisVTkIuJxikpcE9tEzk1FLiIeJS0tTbPTRc6DZq2LiEdp3rw5ixYt4vrrr1eJi5SD9shFxCPEx8cXL4gzZMgQlbhIOanIRcR2TqeT2NhYnnvuOS2YI3KeVOQiYquSy646nU4sq/SCOtdeey2WZZX5tW7dOptSi3gOnSMXEducunZ6WYfTX375ZTIzM0ttmzRpEt9++y2dO3d2V1QRj6UiFxHbJCUlnXN2+ql3hMvNzSUlJYWbb76ZKlX0nzAR/b9ARNzuxIkTVK1alddee43s7GyqV69e7td++OGHHDp0iBEjRrgwoYj30DlyEXErp9NJhw4d+O233/D39z+vEgdYuHAhoaGh9OjRw0UJRbyLilxE3KbonHhISMgFXV6WlZXFsmXLuPnmm0+bFCdSWanIRcQtyjOx7VyWLVvG0aNHdVhdpAQVuYi4XFJSUoUsu7pw4UJatWpFVFS57u4oUimoyEXE5Tp37syoUaMuqsQPHz5c/AeBiPxXuWetW5bVDQgr+RpjzNsuyCQiPmL16tVcddVVhISE8O9//7tcrzmSfYJpn6TyW8ZxYjo2pt8VjQD44IMPyMnJUZGLnKJce+SWZb0D/Au4Guh88kvHtkTkjJxOJ7179+bvf//7eb1u/KLvmLsmjQ837eXed79h3c8HgcLD6hEREbRv394VcUW8Vnn3yKOAcKNFkEWkHEpObHvkkUfO67UpO9OLvzcGvvn1EK1rGZKTk8/7jwKRyqC8Rf4j0BD43YVZRMQHXOzs9I6X1mbVT/uLH0c0qU1ISAgnTpyo6KgiPsE62062ZVnLAAPUADoCXwE5Rc8bYwa5OuD5ioqKMikpKXbHEKmUjhw5QsuWLWnbtu0FT2w7nHWC5z7cwm8Z2cR0bMyQK5u4IKmIZ7Msa4MxplynsM+1R/6vCsgjIpVEjRo1SE5OJiws7IJnp9cKrsqUIR0qOJmI7zrrZDdjzOfGmM+B/kXfl9zmnogi4umcTidTpkwB4IorrrjgEheR81fe68ivL2Nbv4oMIiLeqeic+IoVK8jNzbU7jkilc9ZD65Zl/Q24F2hhWdb3JZ6qAax1ZTAR8XxFJR4dHU1SUhLVqlWzO5JIpXOuc+TvAknAFODREtuPGGPSy36JiFQGp5a4DqeL2OOsRW6MOQwcBkYAWJbVAAgEHJZlOYwxv7o+ooh4oiNHjtCtWzdWrFihEhexUXlXdhtoWVYqkAZ8DuygcE9dRCqZAwcOAHDnnXfy2WefqcRFbFbeyW7PANHANmNMc+A6YI3LUomIR3I6nTRv3py1awunyPj7+9ucSETKW+QnjDEHAT/LsvyMMZ9RuECMiFQSRefEIyIiuOKKK+yOIyInlbfIMyzLcgD/Af7PsqwZQJ7rYvmO7du3c8899xAREYG/vz/XXnvtaWPCwsKwLKvUV8OGDd0fVuQMNLFNxHOVt8hjgCxgHPAh8DMw0FWhfMmmTZtITEykTZs2tGnT5ozj4uLiWLduXfFXYmKiG1OKnNlXX32lEpdimzdv5rrrriM4OJjGjRszadIk8vPz7Y5VqZXrpinGmGOWZTUDWhtj3rIsKxjQybFyGDhwIDExMQAMHTq0eKLQqRo1akR0dLQ7o4mUS1RUFFOmTGH06NEq8Uru0KFD9OnTh/DwcJYuXcrPP//MQw89REFBAc8884zd8SqtchW5ZVl/Be4G6gItgVDgFQonvclZ+PmV96CHiGdJSEggIiKCZs2aMWHCBLvjiAd45ZVXOH78OIsXL6ZmzZpcf/31ZGZmMnnyZCZOnEjNmjXtjlgplbdlxgDdgUwAY0wq0MBVoSqjuXPnUq1aNWrVqsXQoUPZuXOn3ZGkEnM6nQwZMoTHH3/c7ijiQZKSkrjxxhtLFXZsbCzHjx/n888/tzFZ5VbeIs8xxhQvomxZVhUKb28qFSAmJobZs2eTnJzMCy+8wLp16+jRoweHDx+2O5pUQiXvJz5nzhy744gH2bp1K+3atSu1rWnTpgQHB7N161abUkm5Dq0Dn1uW9RgQZFnW9RSuv77MdbEqlxkzZhR/36NHD7p160bHjh2ZN28e48aNszGZVDYlS/xC7ycuvuvQoUPUrl37tO116tTh0KFDNiQSKP8e+aPAfuAH4B4gEXjCVaEqu8svv5y2bdvyzTff2B1FKpGCggKmTp2qEpezsizrtG3GmDK3i3uUd9Z6gWVZS4Alxpj9FfXhlmX1BWZQOAP+dWPMc6c8HwC8DXQCDgI3G2N2VNTnezr9H0PcxRiDn58fiYmJVKlSRSUuZapTpw4ZGRmnbT98+HCZe+riHmfdI7cKTbYs6wCwFfjJsqz9lmVNutgPtizLH5hN4X3Nw4ERlmWFnzJsFHDIGNMKmAb882I/19Ve/88vRPzvx0Q/m8xnW/dd0Hv8+OOP/PTTT3Tq1KmC04mczul0MmjQILKzs6lTp45KXM6oXbt2p50L37VrF8eOHTvt3Lm4z7n2yMdROFu9szEmDcCyrBbAHMuyxhtjpl3EZ3cBthtjfjn5vgspXHhmc4kxMcDkk9/HAy9ZlmUZYzxyot2Pew7zzIotABw+foIx737DFw91Z+UnHwGwZ88eMjMziY+PB6B///589tlnzJ8/nwEDBtC4cWO2bt3KM888Q9OmTRk5cqRdP4pUEiXPieflabFGKc0Yw9PLN7Psu99pUieIq665jjdensGRI0eK/+BbtGgRQUFB9OzZ0+a0lde5ivx24HpjTPEqJsaYXyzLuhX4mMK95AsVCuwq8Xg3cNWZxhhj8izLOgzUA0qtqmJZ1t0UXudO06ZNLyLSxdl/JKfU46zcfH7Z9RvDhg0rtb3ocVpaGpdeein79u1j3LhxZGRkUK9ePfr27cuzzz6razLFpU6d2OZwOOyOJB7m/W/2MG/NDgAOHM3hREgXAgICGDJkCI888gi//PILkydP5sEHH9R/r2x0riKvWrLEixhj9luWVfUiP7usE8Cn7mmXZwzGmFeBVwGioqJs21uPCqtDs3rB7DyYBUDPNvWJurwt5zqAkJyc7I54IsUWL16s2elyTnsOHS/1eH9OVZKTk7nvvvsYOHAgtWvXZvz48UyePNmegAKcu8hzL/C58tgNXFricRPgtzOM2X3y2vVaQPpFfq7L1Aisygf3didh4x6CqvkzOLKJ3ZFEytSyZUsGDBjAO++8oxKXM7rhskuY8/l2sk8UADAwojHh4eGsXLnS5mRSknW2vUXLsvKBY2U9BQQaYy54r/xkMW+jcJnXPcDXQJwxZlOJMWOAK4wxoy3LigWGGGOGn+19o6KiTEpKyoXGEvFpmzZtIjw8XFdESLlt+T2TTzf/QWidIAZHhup3x00sy9pgjIkqz9izzlo3xvgbY2qW8VXjYkr85HvnAfcBHwFbgPeMMZssy3rasqxBJ4e9AdSzLGs78CCF17OLyAVwOp1EREQwd+5cu6OIF2nfqCb3X9eaIVc2UYl7qPKu7OYSxphECheXKbltUonvs4Fhp75ORM5PyYltw4ef9aCWiHgZ3ZpLxMdp2VUR36YiF/Fhu3fv5tZbb1WJi/gwWw+ti4hrNWnShCVLlnD11VerxEV8lIpcxAfFx8cTFBTEn/70J/r162d3HBFxIR1aF/Ex8fHxxMbGMm3atHMuRiQi3k9FLuJDikq8a9eufPDBB7pcSKQSUJGL+IiSJa6JbSKVh4pcxEesXr3alhIvug1qaGgoDoeDTp06sWDBglJjRo8eTbt27XA4HNSpU4drrrmGTz/91G0ZRXyZJruJeLmcnBwCAgKYPn062dnZBAcHu/Xzp06dSvPmzZk2bRohISEkJiYSFxfHgQMHuP/++wE4fvw49913H23btiU3N5c33niDfv368Z///Ifo6Gi35hXxNWdda90baa11qUzi4+OZOHEiK1euJCwszJYMBw4cICQkpNS2uLg41q1bR1paWpmvyc/Pp3nz5tx0003MnDnTHTFFvEqFrbUuIp6r6Jx4aGgo9erVsy3HqSUOEBkZyb59+874Gn9/f2rXrk1u7sXeRFFEdGhdxAt5+sS2tWvXEh4eXmqbMYb8/HwOHz7M22+/TWpqKq+//rpNCUV8h4pcxMt8/PHHHl3iycnJLF269LS7rC1atIgRI0YAUL16dRYtWkSXLl3siCjiU3RoXcTLXHXVVdx7770eWeI7duwgLi6OmJgYRo4cWeq5G2+8ka+//pqkpCQGDx5MbGwsq1atsiWniC/RZDcRL7Fy5Uqio6PdPiu9vNLT0+nevTsOh4NVq1ZRvXr1s46/7rrrOHHiBKtXr3ZTQhHvocluIj7G6XRyww038NRTT9kdBYCcvHx+2H2Y3zKOA5CVlcWAAQPIzc1lxYoV5yxxKJwQ98svv7g6qojP0zlyEQ9X8n7ikyZNsjsOR3PyGPHqen7Yc5gqfhbPxoTz1tNjSE1NZc2aNTRo0OCc72GMYd26dTRv3twNiUV8m4pcxIOVLHFPOSf+wTe7+WHPYQDyCgz3jrmXgxuSmDFjBunp6axfv754bGRkJF999RUvvvgiQ4YMoWnTphw8eJC33nqL9evXs2zZMrt+DBGfoSIX8VBZWVmMGzfOo0ocgFNuxHLk528AGDt27GlD09LSuPTSS6lSpQqPPfYY+/fvp379+nTs2JEvvviCrl27uiWyiC9TkYt4qODgYD777DMaNWrkOSUO/PnKUOI37Oa7XRlU9bdYtHIDN0WGnvU18fHxbkonUvmoyEU8jNPpZOPGjTzzzDO0adPG7jinCa5WhfjRXUn94yj1HNW4pGag3ZFEKjUVuYgHKXlOPCcnh8BAzyzJqv5+hDeuaXcMEUGXn4l4jFMntnlqiYuIZ1GRi9ho4cKFXHnllQQGBjJ8+HDq1q3LG2+84VHnxEXEs6nIRWySkJDAiBEj6NatGxMmTKBdu3YEBQURGxtLQUGB3fFExEtoiVYRm8TGxrJ582a+//57oHCRlGXLlhETE8PmzZtp3769zQlFxC5aolXEC6SlpbFp0yZWrlwJgGVZ1K5dGygsdRFf43Q6GTRoEKGhoTgcDjp16sSCBQtKjcnJyeHBBx+kYcOGBAcH06NHD7RzdnYqchEbOJ1OUlJSMMaQmppKZmYm27Zt44knnqBXr16n3ctbxBdMnToVh8PBtGnTSEhIoFevXsTFxTFr1qziMQ888ABvvPEGkydP5v3338fhcNCnTx927txpY3IPZ4zxqa9OnToZEU/23nvvGX9/f3P11Veb1157zQQEBBjAAKZbt27m0KFDdkcUcYn9+/eftm3EiBEmLCzMGGPMrl27jL+/v3njjTeKn8/OzjaNGzc2Y8aMcVtOTwCkmHL2nvbIRdxo48aNxZeYPfroo4wfP56xY8fy2WefsXDhQtLT0xk8eDD5+fl2RxWpcCEhIadti4yMZN++fQD88MMP5Ofn06dPn+LnAwIC6NGjBytWrHBbTm+jBWFE3CgiIoIZM2Zw++2307NnTwYNGsQ///nP4uc7duxIu3btWLp0KUOGDLExqYh7rF27tvhUUnZ2NgDVqlUrNSYgIICdO3eSlZVFcHCw2zN6Ou2Ri7jBkiVL+Omnn7AsizFjxlCjRg22bt1Kx44dS41r27YtQUFB/PzzzzYlFXGf5ORkli5dypgxYwBo1aoVAF9//XXxGGMMX3/9NcYYMjIybMnp6bRHLuJiRSu2DR48mKETX2T2Zz8TXM2fSxo34Ztvvik1dsuWLRw/fpywsDB7woq4yY4dO4iLiyMmJoaRI0cCcMUVV9C9e3cefvhhGjduTNOmTZk6dSrbtm0DwN/f38bEnkvXkYu4UFGJR0dHM33eIobN3UjR/+VOfLec3z/6N+PHj6dfv3788ccfPP300+Tm5vLjjz9SvXp1e8OLuEh6ejrdu3fH4XCwatWqUr/r27dvZ9iwYWzcuBGAyy67jBtvvJFZs2Zx7Ngxqlataldst9J15CIeoGSJJyUlcehEFUr+3Vylw5947sUZfPLJJ8TExDBx4kQ6duxIcnKySlx8VlZWFgMGDCA3N5cVK1ac9rveqlUrvv32W37++We2bt3KDz/8QE5ODldeeWWlKfHzpSL3Ydu3b+eee+4hIiICf39/rr322rOOHzduHJZl8fDDD7snoA8zxjBnzpziEq9RowadmtXhkpoBxWOualGPCePu4/vvv+fYsWPs2bOHRYsW0aJFCxuTi1SslVv/4J53Unj8gx/YfziLYcOGkZqaSlJSEg0aNDjj61q0aEHbtm05ePAg7733HqNGjXJjau+ic+Q+bNOmTSQmJhIdHU1ubu5Zx27evJm5c+dSs6ZuTXmxjDFYlkVCQgLGmOIboNStXo0P7u3Oeym7CK7mz23RYViWZXNaEdf5fncGf317A/kFhYeiFkx9grQvEpkxYwbp6emsX7++eGxkZCQBAQHMnDmTevXqERoaSmpqKlOmTOGKK65QkZ+FityHDRw4kJiYGACGDh3KgQMHzjj2gQceYOzYsbzzzjvuiueTnE4nc+bMISEhAYfDcdrzjWsHMa5PGxuSibjfxl0ZxSUOsOuHwuIeO3bsaWPT0tIICwvj+PHjPPHEE/z22280aNCAuLg4nnrqKfz8dAD5TPRvxoeV9xc/Pj6eLVu28Oijj7o4kW8rOieem5urtdJFgCtCa+FX4qBTzJQPzrg6WdGVGo888ghpaWnk5OSwa9cu/vnPf+ra8XPQHnkld/z4cR566CGee+45TbC6CKdObNP9xEUgsmkdXr6lE++l7CLEUY2JfdvZHcknqcgruSlTptCoUSNuvfVWu6N4rSVLlqjERc6g7+UN6Xt5Q7tj+DQVeSWWlpbGv/71L1auXKlJVxehffv2DB48mLlz56rERcTtdI68Env00Ufp168f7dq1IyMjg4yMDAoKCsjJySEjI0Pnec9h48aNGGNo27YtTqdTJS4itlCRV2I//fQTixcvpk6dOsVfu3bt4qWXXqJOnTrs2bPH7ogey+l0EhUVxcsvv2x3FBGp5HRo3YccOJrDix//RPqxXGK7NKVX2zMvtgDw+uuvc/To0VLbYmNj6dmzJ3/729+oX7++K+N6raKJbV27duX222+3O46IVHIqch9y11spbNxVeHeg5C37WDTqStK+/QKAPXv2kJmZSXx8PAD9+/cnKur0ZXwDAwO59NJLz7kKXGVVssQTExN1OF1EbKci9xH5Baa4xAHyCgxrfviZscOGlRo37OTjosUXpPz27t3LX/7yF5W4iHgUFbmP8Pez6NCkFt/vPgxAFT+L3p0vP+8Jazt27HBBOt/QsGFDli9fTufOnVXiIuIxVOQ+5PW/RPHChz9xKCuXEV2acnloLbsj+QSn04kxhuHDh9O7d2+744iIlKIi9yENagTywrAIu2P4lKJz4j169GDYsGG63l5EPI4uPxM5g5IT2xISElTiIuKRVOQiZdDsdBHxFipykTJs3LhRJS4iXkHnyEVKyMrKIjg4mGeeeYacnBwCAwPtjiQiclbaIxc5yel00rp1a7Zt24ZlWSpxEfEKKnIR/ntOvEWLFjRq1MjuOCIi5aYil0pPE9tExJupyKVSW7lypUpcRLyailwqtejoaMaPH68SFxGvpSKXSunjjz/m8OHDBAcH88ILL6jERcRrqcil0omPj6d///48+eSTdkcREbloKnKpVOLj44mNjaVr16784x//qND3XrJkCR06dCAgIIDmzZszderUCn1/EZGyqMil0ihZ4hV9TnzNmjUMGTKELl26sGzZMu68804eeeQRpk+fXmGfISJSFut871ft6aKiokxKSordMcTD5OTk0L59e0JDQ10yse3GG2/k+PHjrF69unjbgw8+yJtvvsnevXupVq1ahX6eiPg2y7I2GGOiyjNWe+RSKQQEBLBy5UqXzU7fuHEjffr0KbXthhtu4NChQ6xbt67CP09EpIiKXHya0+nkgQcewBhDWFiYy2anZ2dnn7bXHRAQAMCWLVtc8pkiImBTkVuWVdeyrE8sy0o9+c86ZxiXb1nWxpNfCe7OKd6taMW2b7/9luPHj7v0s1q1asXXX39dattXX30FQHp6uks/W0QqN7v2yB8Fko0xrYHkk4/LctwY0/Hk1yD3xRNvd+qyq8HBwS79vNGjR7N06VJee+01Dh06xEcffcSLL74IgL+/v0s/W0QqN7uKPAZ46+T3bwE32ZRDfFB8fLzbl1298847GT16NH/729+oW7cuQ4YMYdKkSQBccsklLv98Eam87CryS4wxvwOc/GeDM4wLtCwrxbKs9ZZlqeylXIKCgrj22mtdXuLZJ/LZsPMQvx7Mwt/fn5deeon9+/fz/fff88cffxAdHQ1Q/E8REVdw2eVnlmV9CjQs46nHgbeMMbVLjD1kjDntPLllWY2NMb9ZltUCWAlcZ4z5uYxxdwN3AzRt2rTTzp07K+rHEC+ye/dumjRpAoAxBsuyXPZZR3PyGP7KOjb/nomfBf8YfAUjujQtNebOO+/kp59+Ys2aNS7LISK+6XwuP6viqhDGmD5nes6yrD8sy2pkjPndsqxGwL4zvMdvJ//5i2VZq4BI4LQiN8a8CrwKhdeRV0B88TJOp5Nbb72VJUuW0K9fP5eWOMCSb/ew+fdMAAoMPPnqB+xZfYKOHTuSmZnJggUL+Oijj/jiiy9cmkNExK5D6wnAX05+/xdg6akDLMuqY1lWwMnvQ4DuwGa3JRSvUTSxrUuXLlx99dVu+Ux/v9J/KFStWoVFixZx0003MXLkSLKyslizZg0dOnRwSx4Rqbxctkd+Ds8B71mWNQr4FRgGYFlWFDDaGHMX0B74t2VZBRT+wfGcMUZFLqWcOjvdXXcxGxwZyvsbdpOy8xBV/S2eu+cmYuaMcctni4iUpCVaxWtt2rSJiIgIt5d4kbz8AtIOHKNO9WqEOALc+tki4ts84hy5iKuFh4fz73//m+HDh9tyP/Eq/n60vkT3MRcRe2mJVvE6ixcv5rvvvsOyLEaNGmVLiYuIeAoVuXgVp9PJ8OHDeeqpp+yOIiLiEVTk4jVKTmx755137I4jIuIRVOTiFeyanS4i4ulU5OLxjDG8/fbbKnERkTJo1rp4tIKCAvz8/HA6neTl5eFwOOyOJCLiUbRHLh7L6XTSvXt3Dh06RGBgoEpcRKQMKnLxSEXnxKtUqUKVKjpwJCJyJipy8Tia2CYiUn4qcvEoCQkJKnERkfOgIhePEhERQWxsrEpcRKScVOTiEb766isKCgpo1qwZ8+fPV4mL+BCn08mgQYMIDQ3F4XDQqVMnFixYUGpMZmYm48aNIywsjODgYNq3b8/06dPxtRt7uYKKXGzndDrp1q0bL774ot1RRMQFpk6disPhYNq0aSQkJNCrVy/i4uKYNWtW8ZiRI0cyf/58HnvsMZYvX87QoUN58MEHmT59uo3JvYNuYyq2KprYFh0dTVJSkvbERXzQgQMHCAkJKbUtLi6OdevWkZaWRlZWFjVq1GD69Oncf//9xWOGDBnCnj17+PLLL90d2XbncxtT7ZGLbVTicrHKc8i2pOnTp2NZFkOHDnVjSjm1xAEiIyPZt28fAHl5eRQUFFCrVq1SY2rXrq1D6+WgC3TFFgcOHODOO+9UictFmTp1Ks2bN2fatGmEhISQmJhIXFwcBw4cKLVnB7Bv3z6efvpp6tevb1NaKWnt2rWEh4cDULNmTYYPH87zzz9Phw4daNmyJZ9//jnvvfces2fPtjmp59OhdbHN2rVrueKKK1TicsHOavf5MwAAH/5JREFUdci2pFGjRpGbm8uuXbsICQkhPj7enVGlhOTkZK6//nrmzp3LyJEjAcjJyeGWW27h/fffB8CyLKZMmcIjjzxiY1L7nM+hde2Ri1s5nU6OHj3KHXfcQbdu3eyOI17uTIdsly5dWmrb119/zXvvvcfWrVu55ZZb3BVPyrBjxw7i4uKIiYkpLnGA8ePH8+WXXzJv3jxatGjBF198weTJkwkJCWHUqFH2BfYCKnJxm6Jz4t27d+f222/H39/f7kjig0oesoXCu+fdd999TJw4kdDQUBuTSXp6Ov369aNp06bMnz+/ePv333/PnDlz+Pjjj7n++usBuOaaazhy5AgPP/wwd9xxB35+mtJ1Jvo3I25RctnV5cuXq8TFJZKTk1m6dCljxowp3jZv3jz27t3Lww8/bGOyyicrN49vfj3EH5nZhY+zshgwYAC5ubmsWLGC6tWrF4/dunUrAB07diz1HpGRkWRkZHDw4EH3BfdC2iMXl9Pa6eIOZR2yPXz4MI899hgzZ84kKCjI3oCVyL7MbIa+so5f07MIqOLHrJs7MOPRu0lNTWXNmjU0aNCg1PhmzZoB8M0333DjjTcWb9+wYQPVq1cv8xSK/JeKXFxu+/btKnFxqTMdsn322We59P/bu/Poqqq7/+PvTQiEJMxBZJIEngRkFEgJiNaIlAKKoDIVV4WfpVJBhGX1EYpPfwoq0FYsWm3V/hRwQLlBRUxii1FEESzDw2hRUCIEIpBACFMISfbvjwxNIMMNSe65w+e11l3knHvuud971k0+nLP32btDB4YOHUpWVhZQeKvTxYsXycrKonHjxro6VAde3/QDB0+cA+BCXgG/uu83pG1KYsmSJZw4cYJNmzaVbNunTx9iY2OJjY3l3nvvZd68eURFRfHFF1/w5z//mZkzZ2KMceqj+AT1Wpc6c/r06ZLgzs3NpUGDBg5XJP7i26OnOZ+bT892TcnJOc+QIUM4evQoGzduLHO2N3r06Ms6vpX2+eefc8MNN3ii5ICyeO23PJeyr2T52MtTOH/yx3K3PXDgAJGRkfz444889thjrF27luPHj9OxY0fuuecefvvb3wbk347q9FpXkEudcLlcTJ8+nZSUFHr27Ol0OeJHFn20l7+u+w6A+OgWHFv1JJs3/4sNGzYQExNTZtvdu3eTkZFRZt2sWbNo2rQpTzzxBH369LlsEBKpucwzFxj30ka+O36W0AZBvPTLftwYrfv3q0O3n4mjSreJR0ZGOl2O+JGTZ3NLQhxg1fNPcGbHPyq8ZNujR4/L9tGsWTMiIiKIj4/3RMkBqWV4QxIfvJHUzLO0bhxC87DAO6P2JAW51Cp1bJO6VM8YjIHiC4nnD/wvADNnzrxs2+JLtuKMkOAgul7dxOkyAoJuP5Nas379eoW41KmmocE8PLQLxX2f7v/bRxQUFGCtvexRUYivW7dOo7qJX1EbudSaCxcuMH/+fB599FGFuNSpI1nnybmYT6dW4U6XIlInNPuZeFRycjIZGRk0bNiQJ598UiEuda5ts0YKcZEiCnKpEZfLxciRI5k7d67TpYiIBCQFuVyx0h3b/vSnPzldjohIQFKQyxWpqHf6/v37mTp1Kr179yYoKKjcW3xefPFFbr31Vlq2bIkxhnXr1nm2eBERP6Igl2q7ePEiv//978vtnb5nzx6SkpKIiYm5bHCOYsuXL+fEiRNlxlQWEZEro/vIpdqCg4NJSUmhcePGl3VsGzlyJKNGjQJgzJgxl42qBYXTTNarV4/du3ezYsUKj9QsIuKvdEYubnO5XNx7773k5+fTtm3bcnunuzNnsOYVFhGpPfqLKm4pbhPft28fOTk5TpcjIiJFFORSpUs7toWFhTldkoiIFFGQS6USEhI07KqIiBdTkPuQqm7tys3NZdy4cXTq1IlGjRrRqlUrhg8fztatW6/4PSMiIhgyZIhCXETESynIfUhVt3bl5+djjGHOnDkkJibyyiuvcO7cOQYPHsz3339frfc6cOAAAPHx8SQnJyvERUS8lILch4wcOZJDhw7hcrno3r37Zc83atSId955h1//+tcMHjyY0aNHk5iYSE5ODu+//77b75OQkEBMTAzvvvsuAKZ4qqlyfPV9Jrf/5QuGL/mctV8frf6HEhGRGtF95D7kSm7bCgsLIyQkhNzcXLe2T0hIYMKECQwcOJCf/exnlW57OuciU5Zv4XROHgDT39pG4v0/YcfGdQAcPnyY7OzskikjR4wYQWhoKFu2bCE1NZVDhw4B8Nlnn5GRkUFkZCSxsW5N9iMiIkUU5H7IWkt+fj4ZGRk888wzBAUF8Ytf/KLK15UOcXfaxI+dvlAS4gC5eQXs/u4QY8eOLbNd8fKBAweIjIzkL3/5C8uWLSt5/vHHHwdg0qRJLF261M1PKSIioEvrfmnRokUEBwfTpk0bli1bRlJSEh07dqz0Nfv27at27/QOzUPp0vo/27VtGsIt/XtgrS33ERkZCcDSpUvLfV4hLoHE5XJx++23065dO8LDw+nXr1+ZkQ5TU1MxxpT76NKli4OVi7fRGbkfmjx5MkOGDCE9PZ0XX3yR2267jfXr19OtW7cKXxMdHc2yZcsYOXKk2x3bGtSvx4r7BvDahgPk5hcwaWAkTUKCa+tjiPi1xYsXExUVxbPPPktERARJSUlMnDiRjIwMZsyYQZs2bdi4cWOZ15w/f56hQ4cyfPhwh6oWb2SstU7XUKtiY2Ptli1bnC6jzhWPY17VzGF5eXl0796duLg4li9fftnzq1atol27dgwYMKCOKhWR8mRkZBAREVFm3cSJE9m4cWPJXSOXcrlcjBs3jk2bNhEXF+eJMsUhxpit1lq3Og3p0rqfq1+/Pj179iz39jOXy8X48eOZP3++A5WJBLZLQxygT58+HDt2rMLXrFixgqioKIW4lKEg92I7DmVxx4sbGPbn9azefviK9pGTk8O2bduIiooqs770sKtvv/12bZQrIjX05ZdfVtgElp2dTXJyslsdVyWwqI3cS+XmFXDv0s1kni28beyhlTuIalqfb7d+DpR/a9fq1atJTk5m2LBhtG3btqSNPD09nYceeqhk35eOna7BXkScl5KSwurVq3n11VfLff79998nJyeHCRMmeLgy8XYKci+VdT63JMQB8gssO/Yf5J5Kbu3q0qULb7zxBg899BAnT56kTZs2xMXFsWXLljIDyKxatUohLuJFUlNTmThxIqNGjWLy5MnlbrNixQq6d+9Oz549PVuceD0FuZdqFd6Q6zo0Y/uhLABahjVg2IBeVNY5MTIyksTExAqfz8vLo379+rz++utcuHCB8PDwWq9bRKrnxIkTDB8+nGuuuYY33nij3G0yMzP5+OOPS8ZcEClNbeReyhjD8l/1Z9aQaO77aSdW3X89rRo3vOL9uVwuYmNjOXbsGMHBwQpxES9w7tw5brvtNnJzc0lMTKxwiuCEhATy8vJ0WV3KpTNyL9YkJJhZQy6fHKW6SreJN2rUqBYqE5Er8c7mg+xMO0X/qBbc2qM1Y8eOZd++fWzYsIGrrrqqwtetWLGC/v3707lzZw9WK75CQe7n1LFNxDu8sv57nkr6NwBvfnWQ5795k5SkJJYsWcKJEyfYtGlTybZ9+vShYcPCK3BHjhzh888/55lnnnGkbvF+urTuxxITExXiIlQ9HCoUTtlb3nCoOTk5tVLDp9+UvT/8X1+sA2DmzJkMHDiwzCM9Pb1ku5UrVwIwbty4WqlD/I+C3I/17duXe+65RyEuAW/x4sWEh4fz7LPP8sEHH3DzzTczceJEnn/++TLb3XzzzWzcuLHMo/jMuKZiWpf9HZy3Yl2V8xIAzJo1i/z8fNq2bVsrdYj/0RCtfmjDhg3ExcVRv75aTkTAveFQ4+PjiYiIKBmbobadz83n8Q/2sPPwKeKiWvC7EdfSoL7OpaR8GqI1gLlcLm666SYWLFjgdCkiXuNKhkOtbY0aBLFoTC+SZ97I47d3V4hLrdE3yY+U7tg2a9Ysp8upsf379zN16lR69+5NUFAQ8fHxl21jreXpp5+mQ4cONGrUiJ/+9Kds377d88WKzylvONR//vOfhIaGEhoays9//nN27tzpUHUi7lOQ+wl/7J2+Z88ekpKSiImJISam/NvwFi5cyPz583n00UdZs2YN4eHhDBkyhB9//NHD1YovKR4Odfr06SXrbrrpJpYsWcI//vEPXn75ZQ4ePMiNN95Iamqqc4WKuKOizha++ujXr58NNCdPnrTNmjWzN9xwg83Ozna6nFqTn59f8vNdd91lb7rppjLPnz9/3jZp0sQ+8cQTJevOnDljIyIi7Ny5cz1VpviYAwcO2KuuusqOHj260u3S09Nt06ZN7cyZMz1S18qVK+3IkSNt27ZtbVhYmO3bt6996623LtsuNTXVTpgwwTZv3tw2atTI9urVyyYnJ3ukRvEcYIt1M/fUG8oPNGvWjJSUFKKjo/3iTLxYvXqVXzD68ssvyc7OLnNbTlhYGCNHjiQ5OZknn3yyrksUH+POcKjFrr76agYNGsS2bds8UtvixYuJiori2WefJSIigqSkJCZOnEhGRgYzZswA4NChQwwcOJDevXvz2muvERYWxvbt2zl//rxHahTvpCD3YS6Xi/T0dB588EH69u3rdDket3fvXoKCgoiOji6z/tprr+Wdd95xqCrxVu4Oh3opY0wdV1ZozZo1ZTrlDR48mCNHjrB48eKSIH/kkUfo3LkziYmJJf/RHTJkiEfqE++lNnIfVdwm7nK5yMvLc7ocR5w8eZLw8HCCgoLKrG/evDnnzp0jNze3gleKvysosMxb8zXXL0hh/EsbST2eXTIcanJycqXDoRY7evQoGzZsoF+/fh6ouOqe9adOneLdd99l2rRpVV6tksCiM3IfVBziAwYMICkpKaDvFy/vbMkWjY3gqTMp8T6urYd4dUPh/eFHTuVwy5338P0XFQ+H+s033zBnzhzGjh1Lx44dOXjwIAsWLKBevXqO3gFSumf9tm3buHjxIsYYBg0axL/+9S9at27N9OnTmT17tr7vASxwE8BHlQ7x5ORkv2oTr67mzZtz+vRp8vPzy5yVZ2VlERoaSnBwsIPViZMOnjhXZvnw7sLgnjlz5mXbHjhwgJYtW2KtZc6cOWRmZtK4cWPi4+N5//33ueaaazxS86WKe9a/+uqrACV3YkydOpX777+fJ598kk8//ZTHHnuMpk2bMm3aNEfqFOc5EuTGmLHA48C1QH9rbblDsRljhgFLgCDg79bahR4r0kv9+OOPfnWLWU107dqV/Px89u/fT5cuXUrW7927l65duzpYmTjtZ92u5pX1B8jNLwDgf974hLm3dqv0NUlJSZ4ozS2pqalMnDiRUaNGMXnyZAAKCgo/y/Dhw1m4sPBP4c0330xaWhoLFixQkAcwpxpadgN3Ausr2sAYEwS8AAwHugG/MMZU/pvox7KysgCYMWMGn376qV+G+Omci/zho73MeXcX2w9lVbn99ddfT5MmTXC5XCXrzp07x5o1axg+fHhdlipe7roOzXD9ZiAPDv4v/jS2N78bca3TJbmtop71LVq0AArDu7TBgweTlpZGdna2R+sU7+HIGbm19t9QZRtmf2C/tfb7om3fBkYBX9d5gV7G5XJx3333sXbtWmJjY/22Tfy+5VvZ+H0mAO/9bxoJU/qxd3Ph//UOHz5MdnZ2yTjYI0aMIDQ0lNmzZzN//nyaN29O165dWbx4MQUFBSW9fCVw9e7QjN4dmjldRrVU1rP+2mvL/89IcZ8QdYALXN6cCO2AQ6WW04C48jY0xtwH3Ac41p5VV0q3iZe+fOxv8vILSkIcIOdiAeu272fG2LFlthtbtHzgwAEiIyOZPXs2BQUFLFiwgMzMTGJjY1m7di2tW7f2aP0i1ZWbV8DvV+/mi/0ZdL26CQtGd2PyxMKe9Rs2bLisZ31kZCTdu3cnJSWFqVOnlqxPSUmhc+fOhIeHe/ojiJeosyA3xnwMXF3OU3Ottavd2UU568qdqs1a+zLwMhTOfuZ2kV4ukDq21Q+qR+dWYXx3/CwAxsCgPteWnG1UxBjD3LlzmTt3rifKFKk1r3z+PW9vLjxXSTt5nq+WP83uTyruWd+wYUPmz5/PXXfdxSOPPMLQoUNZt24dr7/+OsuXL3fqY4gXqLMgt9bWdJSCNKBDqeX2wJEa7tNnfPnllwET4sVenfwTnljzNSfP5fLLAR3pc01zp0sSqTM/ZJ4ts/zd9o1AxT3rIyMjueOOO1i+fDlPPfUUS5Ys4ZprruGFF17g7rvv9kjN4p0cnY/cGLMOeLi8XuvGmPrAt8AtwGFgMzDRWrunsn36y3zkeXl5LFy4kJkzZwZEiIsEmo+/PsqvX99C8Z/gOcO7MvWmzs4WJV6jOvOROxLkxpg7gOeBVkAWsN1a+3NjTFsKbzMbUbTdCODPFN5+9qq19qmq9u3rQf7hhx/Sr18/2rRp43QpIlLHvtiXwYbvMuh6dWNGXdfO6XLEi3h9kNclXw7y4jbxX/7yl7z22mtOlyMiIg6pTpDrfgUvUXo+8eeee87pckRExEcoyL1A6RAP5BHb9u/fz9SpU+nduzdBQUHEx8eXeT49PZ1HHnmE3r17Ex4eTocOHZg0aRJHjgRMH0gRkct4833kASEvL4+nn3464EMcYM+ePSQlJTFgwIByZy7bunUr7733HlOmTCEuLo6jR4/y+OOPc/3117N7927dRysiAUlt5A6y1mKM4fjx44SEhAR0iEPhWNLFo1ONGTOGjIwM1q1bV/J8VlYW4eHhZUa2+/bbb+nSpQtLly5l0qRJni5ZRKROqI3cB7hcLsaPH09ubi6tWrUK+BCHqoeYbNas2WXD08bExBAaGloyZ7OISKBRkDuguE08PT293EvI4r6dO3dy7ty5kjmbRUQCjYLcwy7t2KZ23StXUFDAzJkziY6OZujQoU6XIyLiCHV286BVq1apd3otmjNnDhs3buSzzz4jODjY6XJERByhM3IP6tChA8OGDVOI14IXX3yRP/7xjyxbtoy4uHInxRMRCQgKcg/49ttvAejfvz8ffvihQryUrHO5JO9KZ+sPJ91+zapVq5gxYwZ/+MMfGD9+fB1WJyLi/RTkdczlctGtWzfefPNNp0vxOsdPX+DW577g/je3cddfv+SFT/dX+Zp169Zx991388ADD/Dwww97oEoREe+mNvI6VLpj2+233+50OV5nzY4jHM46X7L814/30DpzOwCHDx8mOzubhIQEAEaMGMEPP/zA6NGj6dq1K+PHjy8zX3OrVq3o3FkzR4lI4FGQ1xENu1q18JCyX7+GeWcZO7bsvMpjx44FCudj/uqrrzh16hQ7duxg0KBBZbabNGkSS5curdN6RUS8kUZ2qwOpqanExMQQFxenEK/ExfwCpr+5jX9+fZTGIfX56939uCE6wumyREQcV52R3XRGXgciIyNZuXIlt9xyi0K8EsFB9Xj5nliycy4SGhxE/SB12RARqS795axFCQkJfPLJJwCMHj1aIe6mJiHBCnERkSukv561JCEhgQkTJrBo0SL8rblCRES8l4K8FhSH+MCBA0lISMAY43RJIiISIBTkNVQ6xNWxTUREPE1BXkPJyckKcRERcYx6rV+h3NxcGjRowMsvv0xOTg5hYWFOlyQiIgFIQX4FEhISeOyxx0hJSaFdu3YKcRERh1y8eJG0tDRycnKcLuWKhISE0L59+xrN4Kggr6bSbeJNmjRxuhwRkYCWlpZG48aNiYyM9LmOxtZaMjMzSUtLIyoq6or3ozbyalDHNhER75KTk0PLli19LsQBjDG0bNmyxlcTFORu+uijjxTiIiJeyBdDvFht1K4gd9NPfvITpkyZohAXERGvoiCvwmeffcaFCxdo2bIlf/vb3xTiIiLiVRTklXC5XNxyyy3MmzfP6VJERMSLvfTSS0ybNq3Muu7du7N37946f28FeQVKzyc+e/Zsp8sREREvtnPnTvr06VOynJOTw8GDB4mOjq7z91aQl6N0iKtNXEREqrJr1y769u1bZjkmJoagoKA6f2/dR36JM2fOMH36dIW4iIgfeuurg7y+6QdahAUzb1QPOrcKr5X97tmzhzvvvLOkF/qZM2e47bbbamXfVVGQXyI8PJxPPvmEjh07KsRFRPzI5tQT/O69XSXLv162hU8ejq/xfg8dOkSrVq3KtIc/8MADdOrUibNnzzJt2jQaNGhAfHw8d999d43f71K6tF7E5XLx1FNPAdCjRw+FuIiIn/nu2Jkyywcyz5KXX1Dj/e7cuZPu3buXWff111/Ts2dP3n33XcaMGcMrr7zCBx98UOP3Ko+CnP+0iX/00Ufk5uY6XY6IiNSBuE4tCW3wnzbrG6NbUT+o5jG4a9cuunXrVmbdnj176NWrF2lpaXTo0AGgztrLA/7S+qUd2xo0aOB0SSIiUgeiIsJYOXUgq7al0SK0AVNu7FQr+921a1eZ9vATJ05graV169a0b9+etLQ0rrvuOgoKan72X56ADnL1ThcRCSw92jWlR7umtbrPN998s8xyixYtOHbsGAB33nknDzzwAImJiYwcObJW37dYQAf52bNnGTRoEB9++KFCXEREal1YWBivvfZanb5HQLaRZ2RkADB58mQ+/fRThbiIiPisgAtyl8tFVFQUGzZsAKBevYA7BCIi4kcCKsWK28Svu+46evXq5XQ5IiIiNRYwQa6ObSIi4o8CIsg3b96sEBcR8VPWWqdLuGK1UXtABHm/fv1YtGiRQlxExM+EhISQmZnpk2FurSUzM5OQkJAa7cf44oevTGxsrN2yZQsAq1evpnfv3kRGRjpblIiI1ImLFy+SlpZGTk6O06VckZCQENq3b09wcHCZ9caYrdbaWHf24bf3kRe3iY8bN4633nrL6XJERKQOBAcHExUV5XQZjvLLS+vFIT5gwABeeuklp8sRERGpM353ab1z5872hx9+YMCAASQnJ6tNXEREfE51Lq373Rn50aNHFeIiIhIw/O6M3BhzHPjB6TrqQASQ4XQRPkLHqnp0vKpHx8t9OlbVU/p4dbTWtnLnRX4X5P7KGLPF3cssgU7Hqnp0vKpHx8t9OlbVc6XHy+8urYuIiAQSBbmIiIgPU5D7jpedLsCH6FhVj45X9eh4uU/Hqnqu6HipjVxERMSH6YxcRETEhynIvZQxZqwxZo8xpsAYU2EvRmPMMGPMN8aY/caY2Z6s0VsYY1oYY9YaY/YV/du8gu3yjTHbix4feLpOp1X1XTHGNDTGvFP0/FfGmEjPV+kd3DhWk40xx0t9n6Y4Uac3MMa8aow5ZozZXcHzxhjzXNGx3GmM6evpGr2JG8cr3hhzqtR36/dV7VNB7r12A3cC6yvawBgTBLwADAe6Ab8wxnTzTHleZTaQYq2NBlKKlstz3lp7XdHjds+V5zw3vyu/Ak5aa/8LeBZY5NkqvUM1fq/eKfV9+rtHi/QuS4FhlTw/HIguetwH/NUDNXmzpVR+vAA+L/XdmlfVDhXkXspa+29r7TdVbNYf2G+t/d5amwu8DYyq++q8zihgWdHPy4DRDtbirdz5rpQ+jgnALcYY48EavYV+r6rBWrseOFHJJqOA5bbQJqCZMaaNZ6rzPm4cr2pTkPu2dsChUstpResCTWtrbTpA0b9XVbBdiDFmizFmkzEm0MLene9KyTbW2jzgFNDSI9V5F3d/r+4qulScYIzp4JnSfJL+TlXfQGPMDmNMsjGme1Ub++00pr7AGPMxcHU5T8211q52ZxflrPPL2xAqO1bV2M011tojxphOwCfGmF3W2u9qp0Kv5853JWC+T1Vw5zisAVZYay8YY35D4ZWMwXVemW/S96p6tlE4POsZY8wI4H0KmyUqpCB3kLV2SA13kQaUPhNoDxyp4T69UmXHyhhz1BjTxlqbXnTJ7lgF+zhS9O/3xph1QB8gUILcne9K8TZpxpj6QFNq+RKgj6jyWFlrM0stvkKA9idwU8D8naoN1trsUj8nGWNeNMZEWGsrHLNel9Z922Yg2hgTZYxpAEwAAq43NoWfeVLRz5OAy65mGGOaG2MaFv0cAQwCvvZYhc5z57tS+jiOAT6xgTnQRJXH6pI23tuBf3uwPl/zAXBPUe/1AcCp4qYwuZwx5urivinGmP4U5nRmZa/RGbmXMsbcATwPtAISjTHbrbU/N8a0Bf5urR1hrc0zxjwA/AMIAl611u5xsGynLARWGmN+BRwExgIU3bb3G2vtFOBa4CVjTAGFvxgLrbUBE+QVfVeMMfOALdbaD4D/B7xujNlP4Zn4BOcqdo6bx+pBY8ztQB6Fx2qyYwU7zBizAogHIowxacD/BYIBrLV/A5KAEcB+4Bzwf5yp1Du4cbzGAPcbY/KA88CEqv5DrZHdREREfJgurYuIiPgwBbmIiIgPU5CLiIj4MAW5iIiID1OQi4iI+DAFuYgAZWaH22GM2WaMub5ofVtjTEIFr4k0xkwstTzZGPMXT9UsIgpyEfmP4tnhegNzgAVQOCKetXbMpRsXjf4WCUy89DkR8RwNCCMi5WkCnITCs27gQ2ttD2PMZOBWIAQIA0KBa40x2ykcb/wk0NYY8xHQGXjPWvvfHq9eJIAoyEWkWKOiQA4B2lDxJCADgV7W2hPGmHjgYWvtbVB4aR24jsJx7C8A3xhjnrfWHqpgXyJSQ7q0LiLFii+tdwWGAcsrmI98rbW2sslUUqy1p6y1ORSOZ9+xLooVkUIKchG5jLV2IxBB4Vj/lzpbxcsvlPo5H135E6lTCnIRuYwxpiuFE4ZUOusScBpoXPcViUhF9D9lESlW3EYOYIBJ1tr88q+ul9gJ5BljdgBLKeogJyKeo9nPREREfJgurYuIiPgwBbmIiIgPU5CLiIj4MAW5iIiID1OQi4iI+DAFuYiIiA9TkIuIiPgwBbmIiIgP+/8ajmVr8OuoWwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# remark 2 : what closest point sampling (use : infty_sign)\n",
    "# closest point sampling is the sampling that points list from a near point in order\n",
    "# please see below diagram and number of points\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from sklearn.metrics.pairwise import pairwise_distances\n",
    "\n",
    "import tadasets\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from persim import plot_diagrams\n",
    "\n",
    "def dpoint2pointcloud(X, i, metric):\n",
    "\n",
    "    ds = pairwise_distances(X, X[i, :][None, :], metric=metric).flatten()\n",
    "    ds[i] = 0\n",
    "    return ds\n",
    "\n",
    "def get_nearby_perm(X, n_perm=None, distance_matrix=False, metric=\"euclidean\"):\n",
    "    \n",
    "    if not n_perm:\n",
    "        n_perm = X.shape[0]\n",
    "    \n",
    "    idx_perm = np.zeros(n_perm, dtype=np.int64)\n",
    "    \n",
    "    if distance_matrix:\n",
    "        dpoint2all = lambda i: X[i, :]\n",
    "    else:\n",
    "        dpoint2all = lambda i: dpoint2pointcloud(X, i, metric)\n",
    "    ds = np.zeros((n_perm, n_perm))\n",
    "    ds = dpoint2all(0)\n",
    "    dperm2all = [ds]\n",
    "    \n",
    "    for i in range(1, n_perm):\n",
    "        dsn0 = np.nonzero(ds)\n",
    "        dsn0 = dsn0[0]\n",
    "        idx0 = np.argmin(ds[dsn0])\n",
    "        idx = dsn0[idx0]\n",
    "        idx_perm[i] = idx\n",
    "       \n",
    "        dperm2all.append(dpoint2all(idx))\n",
    "        ds = np.where(ds == 0, 0, dperm2all[-1])\n",
    "        \n",
    "    dperm2all = np.array(dperm2all)\n",
    "    return (idx_perm, dperm2all)\n",
    "\n",
    "N = 30\n",
    "X = tadasets.infty_sign(n=N, noise=0.1)\n",
    "\n",
    "originaldm = pairwise_distances(X)\n",
    "idx_perm, dperm2all = get_nearby_perm(X)\n",
    "cpsdm = dperm2all[:, idx_perm]\n",
    "print(\"got perm = \\n{0}\".format(idx_perm))\n",
    "# got perm(idx_perm) is the closest point sampling\n",
    "\n",
    "df = pd.DataFrame(X)\n",
    "\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "plot_diagrams(X)\n",
    "for k, v in df.iterrows():\n",
    "    ax.annotate(k, xy=(v[0], v[1]), size=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}