levelsets_simulations.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# General stuff\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "from IPython.core.debugger import Tracer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgoAAAFkCAYAAABB1xPiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XmclXP/x/HXt70JoywJhSJK7mrGEi1EUm5Zs0y5beUO\nZYmULWSpKPuW+84SMWRJ2akQpWiGkhZ+dpUoSbSo5vv74zPdWuZMc2au61xneT8fj+sx5sx1ru/n\nXI05n/NdPl/nvUdERESkJJWiDkBERESSlxIFERERiUmJgoiIiMSkREFERERiUqIgIiIiMSlREBER\nkZiUKIiIiEhMShREREQkJiUKIiIiEpMSBREREYkp9ETBOberc+5J59wS59xK59xM51xO2O2KiIhI\nxVUJ8+LOue2BKcBE4BhgCbAPsCzMdkVERCQYLsxNoZxzQ4FDvfeHh9aIiIiIhCbsoYcuwAzn3Bjn\n3GLnXKFzrmfIbYqIiEhAwu5RWAV44A7geeBg4B6gl/f+yRLO3wEbovgWWB1aYCIiIumnBrAn8Kb3\nfmlQFw07UVgDfOS9b7vRY/cAB3rvW5dwfjfgqdACEhERSX/dvfdPB3WxUCczAouAuZs9Nhc4Ocb5\n3wKMHj2aJk2ahBiWbKxv377cddddUYeRUXTPE0/3PPF0zxNr7ty5nHnmmVD8XhqUsBOFKcC+mz22\nL/BdjPNXAzRp0oScHK2gTJTs7Gzd7wTTPU883fPE0z2PTKBD92FPZrwLaOWcu9o516h4aKEncH/I\n7YqIiEgAQk0UvPczgJOAPOAz4FrgUu/9M2G2KyIiIsEIe+gB7/1rwGthtyMiIiLB014PQl5eXtQh\nZBzd88TTPU883fP0EOryyHgV7wFRUFBQoAkwIiIicSgsLCQ3Nxcg13tfGNR11aMgIiIiMSlREBER\nkZiUKIiIiEhMShREREQkJiUKIiIiEpMSBREREYlJiYKIiIjEpERBREREYlKiICIiIjEpURAREZGY\nlCiIiIhITEoUREREJCYlCiIiIhKTEgURERGJSYmCiIiIxKREQURERGJSoiAiIiIxKVEQERGRmJQo\niIiISExKFEQkaaxbB/Pmwfffg/dRRyMioERBRJLEuHHQqBE0aQJ77AGHHQaffhp1VCISaqLgnLvB\nOVe02TEnzDZFJPU89hicdBIccABMnAjPPw+rVkHr1va9iESnSgLamA0cBbji79cloE0RSREffQS9\nekGPHvDww1Cp+ONL585w8smWQEybBk2bRhunSKZKxNDDOu/9L977n4uPXxPQpoikgHXroGdP+Mc/\n4IEH/k4SALKyrGehfn0480xYuza6OEUyWSIShX2ccwucc18550Y75+onoE0RSQEjR8Ls2daTUK3a\nlj/fZht44gmYNQsGD058fCISfqIwDTgHOAa4ANgLmOycqxVyuyKS5Natg9tug9NPh9zc2Ofl5sKA\nATB0qK2GEJHECjVR8N6/6b1/wXs/23v/NnAsUBs4Lcx2RST5jRkD334LV1219XOvugq22w6uuy70\nsERkM4mYzPg/3vvlzrkvgL1LO69v375kZ2dv8lheXh55eXlhhiciCTR8OHTqBM2bb/3cbbeFQYPg\nwguhf39o1iz8+ESSWX5+Pvn5+Zs8tnz58lDacj6BVU2cc9sA3wPXe+/vL+HnOUBBQUEBOTk5CYtL\nRBJr5kxo0QLGj4cuXcr2nL/+gr33hnbtYPTocOMTSUWFhYXk2jherve+MKjrhl1HYZhzrp1zbg/n\n3GHAWGAtkL+Vp4pIGhs1CnbayXoUyqpaNbjySnjmGfjmm/BiE5FNhT2ZcXfgaWAe8AzwC9DKe780\n5HZFJEmtXQtPPQXdu0PVqvE9t0cPqFMHhg0LJzYR2VLYkxnzvPe7e+9reu8beO+7ee/1WUAkg02a\nBD//DGedFf9zs7Lg0kvh0Ufhl1+Cj01EtqS9HkQkoV56Cfbay+YolMcFF4Bz8MgjwcYlIiVToiAi\nCVNUZBMYTzjB3uzLY4cdIC8PHnrIajGISLiUKIhIwsyYAQsXwoknVuw6ffpY8aVXXgkmLhGJTYmC\niCTMuHE2GbF164pdJycHDj0U7t9ikbWIBE2JgogkzFtvwTHHQJUASr316WNbUM/RxvUioVKiICIJ\n8euvUFAAHToEc72uXaFuXfUqiIRNiYKIJMS774L3wSUK1arBv/9tu0uGVLlWRFCiICIJMnGilWBu\n0CC4a/bqBWvWWLIgIuFQoiAiCTFhQnC9CRvsthucfLINPxQVBXttETFKFEQkdD/+CF98AUcdFfy1\n+/Sxa0+cGPy1RUSJgogkwMSJVmCpffvgr92mDfzjH5rUKBIWJQoiErrJk+GAA6yqYtCcs16Fl1+G\nb78N/voimU6JgoiEburUihdZKk23bpCdbWWdRSRYShREJFRLl8K8eeEmCrVqwXnnwciRsGpVeO2I\nZCIlCiISqg8/tK+HHRZuOxdeCMuWwdNPh9uOSKZRoiAioZoyBXbZBfbcM9x29t4bjjsOhg2D9evD\nbUskkyhREJFQbZifUN5tpeNxzTUwfz689FL4bYlkCiUKIhKav/6Cjz4Kf9hhg1atbAnm4MFWLlpE\nKk6JgoiE5tNPYfXqcCcybu7aa6GwEN54I3FtiqSzADZ7FREp2fTptnlTixaJa/PII6FtWxgwADp2\nhMqVg7v2rFnw3nu2CVXDhtC5M9SuHdz1RZKRehREJDQzZkDz5lC9euLadA7uuAM++wweeyyYa77z\nDhx8sL2Wfv3g3nuhe3eoXx9uv12TJyW9KVEQkdDMmAEHHpj4dg86yIowDRwIK1aU/zqrV9sOlUce\nCVWq2CTJFSvg559h4UI4/3y46io480xYuza4+EWSiRIFEQnFH3/A3LnRJAoAQ4bA77/DlVeW7/nf\nfWf7SIwaBQ8/bMs8TzjBhlIA6tWDu+6C556D55+Hyy4LLnaRZKJEQURC8ckntvIgqkShQQMbgnj4\nYXj99fie++abkJNjVSWnToV//zv28s5TToEHH7RjzJiKxy2SbBKWKDjnrnLOFTnn7kxUmyISnRkz\noEYNaNo0uhh69bIJh926WRnprVm/HgYNsucccggUFFjCsDXnn28Jw8UXW3Ihkk4Skig45w4C/g3M\nTER7IhK9GTOgZUsb24+Kc1bSebfdoFOn0pOFBQvg6KMtUbjxRnjlFahTp+xt3XcfrFljRZ9E0kno\niYJzbhtgNNAT+C3s9kQkOXz8sU0qjNr228Nrr9nGUa1awZNPwrp1f/986VIr0LTffpZITJwI118P\nleL861ivnk2efOQR+OqrYF+DSJQS0aPwAPCy935SAtoSkSTw22/w5ZfRzU/YXIMGtjlVp05w1lnW\nw9ChAxx6KOy6q/UinHsuzJljlR3L68ILYaed4Oabg4tdJGqhJgrOuTOAFsDVYbYjIsmlsNC+Jkui\nALDddvDMMzbv4LzzrFDSPvvA8OHw449WG2H77SvWRlYWXH219Vp8910wcYtELbTRQ+fc7sDdQAfv\nfVwrjPv27Ut2dvYmj+Xl5ZGXlxdghCISlhkzYJttoHHjqCPZUk5O2SYoltd559nQxYgRtkRTJAz5\n+fnk5+dv8tjy5ctDacv5kHZOcc6dALwIrAc2LCyqDPjix6r7zRp3zuUABQUFBeSE+X+yiITq9NNh\n0SKYPDnqSKLRt6/1KvzwA9SsGXU0kikKCwvJzc0FyPXeFwZ13TCHHiYAB2BDD82LjxnYxMbmmycJ\nIpI+Zs60FQ+Z6qKLbJKk6ipIOggtUfDe/+m9n7PxAfwJLPXezw2rXRGJ1sqVNpHxH/+IOpLo7LOP\nTYp84omoIxGpuERXZlQvgkiamzMHiooyO1EAW13xzjs2/CCSyhKaKHjvj/TeX57INkUksWbOtBoE\n++8fdSTROuUUq0z51FNRRyJSMdrrQUQCNWuWdb1nZUUdSbS23RZOPFGJgqQ+JQoiEqhZszTssMGp\np8Ls2TZnQyRVKVEQkcB4b0MPShTMMcfY8sixY6OORKT8lCiISGAWLIBly5QobJCVZcmCEgVJZUoU\nRCQws2bZ1+bNo40jmZx0EkybBgsXRh2JSPkoURCRwMyaZXsqNGgQdSTJ47jjoHJlGDcu6khEykeJ\ngogEZsP8BOe2fm6mqFMHjjhCiYKkLiUKIhKYWbM07FCSY4+F996DVauijkQkfkoURCQQq1fD/Pma\nyFiSTp3s/rz3XtSRiMRPiYKIBGLOHFi/XolCSZo0gd13hzffjDoSkfgpURCRQMyaZXMTmjWLOpLk\n45z1KrzxRtSRiMRPiYKIBGLWLGjUCLbZJupIklOnTjBvHnz7bdSRiMRHiYKIBEKlm0t31FG2TFLD\nD5JqlCiISIWpdPPWbb89HHIITJgQdSQi8VGiICIV9tNPsGSJEoWtad8e3n0XioqijkSk7JQoiEiF\nqXRz2bRvbwnV559HHYlI2SlREJEKmznTJjHuuWfUkSS3ww6DatXgnXeijkSk7JQoiEiFzZoFBxwA\nlfQXpVQ1a0KrVkoUJLXof2sRqTCteCi79u2tQqPmKUiqUKIgIhXy118wd67mJ5RV+/awbJkN14ik\nAiUKIlIhc+fCunXqUSirVq2gRg0NP0jqUKIgIhWyYcXDAQdEG0eqqF7dJjUqUZBUoURBRCpk1izY\nay/YbruoI0kdRxwBkyfbJloiyS7URME5d4FzbqZzbnnxMdU51ynMNkUksVSRMX5t28Lvv6uegqSG\nsHsUfgAGADlALjAJGOecaxJyuyKSIFrxEL+DD4aqVeH996OORGTrQk0UvPeveu/f8N5/5b3/P+/9\ndcAfQKsw2xWRxFi82A4lCvHJyoLcXPjgg6gjEdm6hM1RcM5Vcs6dAWQBHyaqXREJz2ef2VctjYxf\nmzbWo+B91JGIlC70RME518w5twJYAzwInOS9nxd2uyISvpkz7dNxw4ZRR5J62raFBQvgu++ijkSk\ndInoUZgHNAcOBh4CnnDO7ZeAdkUkZLNmQbNmULly1JGknsMOs68afpBkVyXsBrz364Cvi7/9xDl3\nMHApcGGs5/Tt25fs7OxNHsvLyyMvLy+0OEUkfrNmwYEHRh1FatpxR2jSxBKFM8+MOhpJNfn5+eTn\n52/y2PLly0NpK/REoQSVgOqlnXDXXXeRk5OToHBEpDzWroU5c6BHj6gjSV1t26pHQcqnpA/PhYWF\n5ObmBt5W2HUUBjvn2jrn9iieqzAEOBwYHWa7IhK++fNtnweteCi/Nm2slsLSpVFHIhJb2HMUdgZG\nYfMUJmC1FDp67yeF3K6IhEylmyuuTRv7OnVqtHGIlCbUoQfvfc8wry8i0Zk1Cxo0gNq1o44kde25\nJ+y2mw0/dOkSdTQiJdNeDyJSLirdXHHOWa+C5ilIMlOiICLlotLNwWjTBj7+GFatijoSkZIpURCR\nuC1ZAgsXKlEIQtu2toLk44+jjkSkZEoURCRuGyYyqnRzxTVrZlt0a/hBkpUSBRGJ28yZULMm7LNP\n1JGkvsqVrUqjEgVJVkoURCRuM2eqdHOQWre2JZLr10cdiciWlCiISNy04iFYbdrA8uVWfEkk2ShR\nEJG4bCjdrPkJwTn4YKhSRcMPkpyUKIhIXObNs9LNShSCk5UFublKFCQ5KVEQkbjMnGlfNfQQrDZt\nYMqUqKMQ2ZISBRGJy6xZsMcesP32UUeSXtq0ge+/t0MkmShREJG4zJypYYcwtG5tX9WrIMlGiYKI\nxEWJQjh22gn23VfzFCT5KFEQkTJbvNgOJQrh0AZRkoyUKIhImW2YyKhEIRxt2sBnn8Fvv0Udicjf\nlCiISJnNnAnbbAMNG0YdSXpq0wa8hw8/jDoSkb8pURCRMps5Ew44ACrpL0coGjWCunU1/CDJRf+7\ni0iZqXRzuJzTPAVJPkoURKRM1qyxqoyanxCuNm3go4/sfoskAyUKIlImc+fCunVKFMLWpg2sXg2F\nhVFHImKUKIhImWxY8XDAAdHGke5atIBatTT8IMlDiYKIlElBATRuDNtuG3Uk6a1KFTjkECUKkjyU\nKIhImRQW2g6HEr4NG0QVFUUdiYgSBREpg/Xr4ZNPICcn6kgyQ5s2sHQpzJ8fdSQiIScKzrmrnXMf\nOed+d84tds6Ndc41DrNNEQneF1/AypXqUUiUVq2sVoWGHyQZhN2j0Ba4DzgE6ABUBd5yztUMuV0R\nCVBBgX1t2TLaODLFttvapMb33486EhGoEubFvffHbvy9c+4c4GcgF1CuLJIiCgutauD220cdSeZo\n1w5efNFKOjsXdTSSyRI9R2F7wAO/JrhdEamAggLNT0i0I4+E77+Hr7+OOhLJdAlLFJxzDrgb+MB7\nPydR7YpIxRQV2URGzU9IrHbtbJ7CO+9EHUl8iopg2jR47DEYNcrqb3gfdVRSEaEOPWzmQaAp0Hpr\nJ/bt25fs7OxNHsvLyyMvLy+k0EQklv/7P1ixQj0KiZadbcnZpEnQs2fU0Wyd95Yc3Hrrlr0gzZrB\nHXdAx47RxJaO8vPzyc/P3+Sx5cuXh9KW8wlI9Zxz9wNdgLbe++9LOS8HKCgoKCBHf5VEkkJ+PnTr\nBkuWwA47RB1NZrnqKvtUvnBhcs9TWL4czjgD3ngDTjsNeve2olHr19uEzCFD4L337PG777aiUhK8\nwsJCcq3rL9d7H1gR8NCHHoqThBOA9qUlCSKSnAoLYY89lCREoX17+Okn24wrWS1ebMMk06bB66/D\ns8/a99WrQ1YWHHOMDZ88+CCMGAEnnwxr10YdtcQj7DoKDwLdgW7An865usVHjTDbFZHgFBRofkJU\n2rSxT9/JOk/h99+hc2f45RerJNmpU8nnOQcXXgivvGK9Duedp6qTqSTsHoULgO2Ad4GFGx2nhdyu\niATAe5VujlKtWtaFP2lS1JFsqagI8vLgq6/gzTehadOtP6dTJ3jySRg9GoYPDz9GCUbYdRRUIlok\nhX39tY0/a8pQdI480rrti4psFUSyuOMOeO01O+LZUfT0020lxDXXwKGHQtu24cUowUiiXzsRSTbT\np9vXgw6KNo5M1qGD7ftQGNjUtIqbNs3e6AcMsKGHeN10kyUJ555rpcEluSlREJGYpk2DffbRRMYo\nHXoobLedTRRMBmvWwDnn2HDUzTeX7xpVqsAjj8CPP8KgQYGGJyFQoiAiMU2fbmPkEp2qVa1XIVkS\nhaFDbUjqkUcstvJq3BhuuMGGMJKpt0S2pERBREq0erVVZGzVKupIpHNnS9p+jbj4/bx5MHgw9O8P\n++9f8ev16wf77QeXXabqjclMiYKIlOiTT2y9uxKF6HXqZJMZ33oruhi8hwsugAYN4Nprg7lm1arW\no/D++/DSS8FcU4KnREFESjR9OtSoAf/4R9SRyO6728qCKIcfxo2z6or33Qc1awZ33WOOsaN/f/jr\nr+CuK8FRoiAiJZo2zSasVWQcWoLTqZMVK4qiUNHatbbC4eijYxdVqojhw60ew8iRwV9bKk6JgoiU\naNo0DTskk86d4eefbUgo0UaOhC+/hNtvD+f6zZpZ8aYhQ2xVhSQXJQoisoWffoLvvlOikExat7Zl\nkuPHJ7bdFSvgxhvhX/+CFi3Ca2fgQFiwAB59NLw2pHyUKIjIFjYUWlKikDyqVYMuXeCFFxLb7rBh\nVp3zllvCbWe//axXYfBg9SokGyUKIrKFqVNh111tEp0kj1NOgc8/h/nzE9PewoW2KuGyy6B+/fDb\nGzjQ2nzkkfDbkrJToiAiW3j/fdXgT0bHHGNbNyeqV+HGG22Fw9VXJ6a9/fazvSCGDYN16xLTpmyd\nEgUR2cTKlTBjBrRrF3UksrmsLDj22MQkCp9/bp/sBw6E7Ozw29ugf3/49lt4/vnEtSmlU6IgIpuY\nPt2Ww6lHITl17Wolj7/5Jtx2rroK9twTLrww3HY216IFdOxoKyxUrTE5KFEQkU28/z7Urh1MiV4J\n3j//CbVqwejR4bXx7rvwyiu2XLFatfDaiaV/f1sGOnFi4tuWLSlREJFNTJ4MbdpAJf11SErbbAOn\nngqPPx5O8aWiIrjiCjj4YGsnCkceCTk54dVtkPjoT4GI/M/atfDhhxp2SHbnnGM7OH7wQfDXfuop\nG9q4805wLvjrl4Vz1qvw9tvRFJiSTSlREJH/KSy0yYyayJjc2raFhg3hsceCve7KlXDNNbYMs3Xr\nYK8dr1NOgb32Uq9CMlCiICL/8847Nv6dkxN1JFKaSpXg7LPhuefg99+Du+5dd8HixXDbbcFds7yq\nVIHLL7fX+O23UUeT2ZQoiMj/TJgARxyhjaBSwXnnWQXDoHoVfvoJhg6FPn2gUaNgrllR555rSzPv\nvjvqSDKbEgURAWDVKhvz7tAh6kikLHbfHU47De65B9avr/j1+vWD6tXhuusqfq2g1KoFvXvbplTL\nlkUdTeZSoiAiAEyZYp9QlSikjr59rZ7C2LEVu86kSTaJcdgwqFMnmNiC0qePVWl86KGoI8lcShRE\nBLBhh7p1VT8hlRx4IBx1FAwaVP5ehdWrrahS27a2miLZ7Lyzzce4916LVRIv1ETBOdfWOTfeObfA\nOVfknDs+zPZEpPwmTLDehKiWxEn53HorzJ4N+fnle/4119iW4iNGJO+//RVXwM8/h1tkSmILu0eh\nFvApcBGgYpwiSWrpUlsaqWGH1HPIIXDiiTa34I8/4nvuW2/ZSoehQ6Fp03DiC0LjxnDCCbaTZRhF\npqR0oSYK3vs3vPfXe+/HAUmaq4rIW29ZXf2jj446EimP4cPtE/e115b9OQsWWJf+0UfDJZeEF1tQ\n+vWDefPg1VejjiTzaI6CiPDyy9CyJey2W9SRSHk0agS33AL33WdDSFvz559w/PG2DPaJJ1KjXHfr\n1nDooTbhUhIrBX49RCRM69bB669Dly5RRyIVcemltuti164wZ07s837/HTp3hi++sARxl10SF2NF\nXXmlbVo2fXrUkWQWJQoiGW7KFPjtNzjuuKgjkYqoXBnGjIEGDaB9e5g6dctzZs2yDb9mzbLhpubN\nEx9nRRx/POyzj3oVEq1K1AGUpG/fvmRnZ2/yWF5eHnl5eRFFJJK+XnnFPlXm5kYdiVTUdtvZ1swn\nnmgJwemnQ6dO9rOJE+GZZ2Dffa2wVrNm0cZaHpUrW1nniy6Cr75KngqSUcjPzyd/s6Uuy5cvD6Ut\n531iFiM454qAE73340s5JwcoKCgoIEfF5kUSYr/9bA39f/8bdSQSlHXr4OGH4YEHYO5ce2zffaFn\nT7j4YqvAmKpWrYI99rAtsB94IOpokkthYSG5lvHneu8Lg7pu2HUUajnnmjvnWhQ/1LD4+/phtisi\nZTNvHsyfr/kJ6aZKFSt9PGeOLZn880/7t95QpjmV1axp1RofewyWLIk6mswQ9hyFA4FPgAKsjsId\nQCEwKOR2RaQMnn3Wuqs7dow6EglLrVqQlRV1FMG66CL7qh6FxAi7jsJ73vtK3vvKmx3nhdmuiGyd\n9zZmfeKJUKNG1NGIlN2OO9rumffdBytWRB1N+tOqB5EM9dln1h19+ulRRyISv/79LUm4//6oI0l/\nShREMtSzz0Lt2irbLKmpQQObnDlsmNWGkPAoURDJQN5bonDyyVCtWtTRiJTP1VfDypW2s6SER4mC\nSAaaPNnWoZ95ZtSRiJTf7rtDr162WdRvv0UdTfpSoiCSgUaOtAp3hx8edSQiFXPVVbB6te2CKeFQ\noiCSYZYtg+eft/Fdpz1dJcXVq2d1Fe64AxYtijqa9KREQSTDPPWUVe47++yoIxEJxrXX2hLfgQOj\njiQ9KVEQySDew4gRVomxbt2ooxEJxvbbw403wqOPwsyZUUeTfpQoiGSQN96Azz+3LYlF0kmvXtC4\nMVxxhSXEEhwlCiIZZNgwOOggaNcu6khEglW1Ktx5p+2SOWZM1NGkFyUKIhmioADeeccq2mkSo6Sj\nY4+FU06BSy6xSbsSDCUKIhliyBBo2BBOOinqSETCc++9tlxywICoI0kfShREMsD06fDCC3D99VC5\nctTRiIRn111h6FD473+tB00qTomCSJrz3j5dNWumSoySGXr1giOOgH/9C5YujTqa1KdEQSTNvfYa\nvPeefcpSb4JkgkqVYPRoWLXKCotpFUTFKFEQSWN//gkXXwzt29tEL5FMsdtuVlfhpZfgoYeijia1\nKVEQSWPXX29lbf/zH610kMxzwgmWKF96Kbz7btTRpC4lCiJpato0uPtuGDQI9t476mhEonHHHbb5\n2SmnwP/9X9TRpCYlCiJpaMkSOO00K650+eVRRyMSnapV4bnnYKed4Oij4Ycfoo4o9ShREEkz69dD\nt242keu556BKlagjEolW7drw9ts2qfGoo2DhwqgjSi1KFETSiPdwwQVWxvaZZ6B+/agjEkkO9evD\npEmWQB92GMyfH3VEqUOJgkia8B769YORI+Gxx+yTk4j8rWFDmDoVatWyZGHChKgjSg1KFETSwNq1\ntl78zjvh/vvhrLOijkgkOdWvDx98ALm50LGjbU+9bl3UUSU3JQoiKe6nn6BTJ3jySRg1Cnr3jjoi\nkeRWuza8/rolCTffbJN+P/446qiSV+iJgnOut3PuG+fcKufcNOfcQWG3KZIJvLf9Gw44AD7/HN56\nSz0JImVVubLVGZk+3b4/+GBbKTRnTrRxJaNQEwXn3OnAHcANQEtgJvCmc27HMNsVSXeffGJzELp2\nhdat4bPPrLa9iMTnwAOtN+GRRyxpaNbMqpi+9JIN6Un4PQp9gYe990947+cBFwArgfNCblck7Xhv\n1eW6dIGcHKu4+OqrMHasrREXkfKpUgXOOw++/NJ2nVyyxLZj33ln6N4dnnoKvvkmc/eMCG2FtXOu\nKpALDN7wmPfeO+cmAIeG1a5Iuvn2W3j6advkZu5c+8Tz6KO2E2TVqlFHJ5I+qlWDHj3smDnThvZe\nftn+/wOoWxdatoTGjWGffWwVxS67WEKx005QvXq08YclzFIsOwKVgcWbPb4Y2Le0Jw4ePJfddrMs\nr1o1+2O44esOO9SgRYum1K4NdepAdrbtFLaxOXPmsHr16pjXr1evHvXq1Yv581WrVjF37txSX1yT\nJk2oWbP6ZsDWAAAaP0lEQVRmzJ8vWrSIRYsWxfx5jRo1aNq0aalt6HWYTHwdv/4Kzz9vycH770NW\nFpx8MtxwwyL23nsRztlwQ7K/jlhS7d8jFr2Ov6Xj6zjxRDuWLbP/3z77DH78sR5vvVWPhx7acmhi\n222hVq1VVKkylxo1oGZNqFHD3rsqV7b3qkqVoHbtJlSrVpPKldnk2LAfy8qVi1i1KvbrqFKlBnXq\n/P06Njzvl19KfWnl570P5QDqAUXAIZs9fhvwYYzn5AC+9KOptw4gO6pW9b5RI++POsr7Hj28Hz7c\n+z32aFrqNW644QZfmtmzZ28lBvzs2bNLvcYNN9xQ6vObNm1a6vO9975pU72OTHsdw4bN9l262O91\npUred+zo/ZNPer9iRWq9jnT599Dr0OuI9TrWrvX++++9nzHD+9de8/7xx70fNsz7iy7a+uto1Wq2\nb9/e+3btvG/d2vtDD/X+4IP/PnbdtfTXUaNGU7/XXk/77Owumxy1arXbcE6OD/D93PmQBl2Khx5W\nAqd478dv9PjjQLb3/qQSnpMDFIwePZomTZr87/GiIlvnumYN/PVXDXbcsSnLlsGvv8LixfDdd9Y9\n+803NmN15co5wGr22QcOPdQme7VsaRkbpGZmWxK9jr+l8uv44QcYMwbGj4c//mhCq1Y16dbNZmDX\nrZs6r2NjqfzvsTG9DqPX8bdkfh2FhYXk5uYC5HrvC0sNIg6hJQoAzrlpwHTv/aXF3zvge+Be7/2w\nEs7PAQoKCgrIyckpV5vr18MXX9gs1gkTbMnY4sWw++62dKxXL2jQoAIvSiQg8+fDddfZOGjt2vDv\nf1vRpEaNoo5MRFJRWIlC2Kse7gTOd86d5ZzbDxgBZAGPh9Vg5crQpIklBU88YZt/TJtmM8UffND+\nCPfsCV9/HVYEIqX77Te48ELYf39bjjViBPz4IwwZoiRBRJJPqImC934M0A+4CfgE+AdwjPc+rCkX\nW6hUCQ45xJKEH36AoUPhlVcsmRg4EFauTFQkIvDaa5YgPP003H679X79+9826UlEJBmFXpnRe/+g\n935P731N7/2h3vsZYbcZyzbbwBVXWG/CgAH2h7plSyteIxKm9evh2mvhn/+E5s2tkuLll9uMaBGR\nZJaRez1kZcFNN9k62Vq1oFUrePjhqKOSdLVqlRVvGTrUjldftTkzIiKpICMThQ322w8+/BDOPx8u\nuMB6G9avjzoqSScrVlgvwoQJNuQ1YMDfa55FRFJBmAWXUkL16rYtb+PG0Lcv/PwzPP7430spRcpr\nzRo4/ngoLLTVN23aRB2RiEj8Mj5R2OCSS2zNevfu9onvsceULEj5FRXZypsPP7TeBCUJIpKqlChs\n5PTT7Wu3blYe+u67o41HUtcNN8Bzz1kZZiUJIpLKlChs5vTTYelS6N3bNvy45JKoI5JU8+abcOut\ncMsttj+DiEgqU6JQgosusiWUl10GTZtChw5RRySpYuFC29WxUye46qqooxERqbiMXvVQmttuswSh\ne3f74y+yNd5bifAqVawq6Oa7moqIpCL9KYuhcmXb4rdKFcjL07JJ2bpnn7UlkA89BDvuGHU0IiLB\nUKJQip13hvx8eP99uPfeqKORZLZ0KVx8MZx6qu1fLyKSLpQobEW7djah8ZprrC6/SEkGDbK6CUoo\nRSTdKFEog8GDYbfd4LzzNAQhW5o3zzYdu/Za2GWXqKMREQmWEoUyyMqyAkxTpmhPCNlS//62d8Ol\nl0YdiYhI8JQolFHbttajMHCgjUeLgFVefPllGDJEO0GKSHpSohCHwYNh3Tq4/vqoI5FkMWgQ7L//\n31U9RUTSjRKFONStCzfeCCNG2BbVktmmT7cqjAMHqmaCiKQv/XmLU58+ttPk5ZdHHYlE7aaboEkT\n6No16khERMKjRCFOVavC0KEwaZLtCiiZ6fPP4bXX4OqrtcuoiKQ3JQrlcPzx0KqV1VbwPupoJAr3\n3gv16mlugoikPyUK5eCczXL/+GMYOzbqaCTRli6FJ5+0zcOqVYs6GhGRcClRKKcjjoCOHa3Ijoow\nZZb//heKimwDKBGRdKdEoQJuvtmq8r3wQtSRSKKsXQv332+7iu60U9TRiIiET4lCBRx8sPUq3HKL\nfcJMJ97b0r8ePeCww6B9e6tAOH9+1JFF6+WXYcEC2/9DRCQThJYoOOeucc5Ncc796Zz7Nax2onbd\ndfDZZ7a9cLr48ks48kjo1MkqD+69t316fvxxaNoU+vWzT9aZ6JFHLEFs3jzqSEREEiPMHoWqwBjg\noRDbiFzbtrbD5C23pMcKiEmT7I3wxx9t+d/nn8MTT8CYMfDDDzaJ8557bOXHH39EHW1iLVgAb7xh\nvSwiIpkitETBez/Ie38P8FlYbSSL666zFRBvvx11JBUzcSIceywccoi9ns6dbYXHBtWr2/DD66/D\nBx/Y0sB166KLN9Eef9z2czjjjKgjERFJHM1RCECHDvYp/NZbo46k/D79FE44wVZzjBsH228f+9wO\nHeD5520OQ//+CQsxUkVF8OijcOqpsN12UUcjIpI4ShQC4BxcdRVMngwffRR1NPFbtgxOPtlKU7/w\ngvUcbM0xx8Dw4XDXXZlRofK99+DrrzXsICKZJ65EwTk3xDlXVMqx3jnXOKxgk9nxx9ukv+HDo44k\nPt7bm99vv1mSUKtW2Z97ySU26fHcc+H338OLMRmMGmX/vm3aRB2JiEhiVYnz/OHAY1s55+tyxvI/\nffv2JTs7e5PH8vLyyMvLq+ilQ1O5sm0U1aePffJs2DDqiMomP9+qSz7/POy1V3zPrVTJuuObNLFJ\njkOGhBNj1FatghdftH/fjedsiIhEJT8/n/z8/E0eW758eShtOR/yVH3n3NnAXd77OmU4NwcoKCgo\nICcnJ9S4wrByJTRoAN262V4AyW7RIth/fxtG2Oz3LS433mhJwty5qZMgxeOFF2yHyHnzYN99o45G\nRKRkhYWF5ObmAuR67wuDum6YdRTqO+eaA3sAlZ1zzYuPODq3U0tWFvTubWvtf02ByhFXXGG7Yd53\nX8Wu07+/1Vm4/vpg4ko2+fnQsqWSBBHJTGFOZrwJKARuALYp/u9CIDfENiPXu7ft/TBiRNSRlO79\n9+0N8LbbYMcdK3atrCzbbjk/34o1pZPff7diWkk86iUiEqow6yic672vXMIxOaw2k8HOO8PZZ9vQ\nw5o1UUdTsvXr4eKLbUnnWWcFc80ePey1Dx0azPWSxUsv2b+jaieISKbS8sgQXH45/PwzPPVU1JGU\nbORImDnTkplKAf0G1KgBV15pVRy//TaYayaD/Hyrvlm/ftSRiIhEQ4lCCPbd15ZLDh+efJtFrVgB\nAwfCv/5lFRiD1KuXFSOq6JyHZLF8udWIOO20qCMREYmOEoWQ9OtnqwBeey3qSDZ1zz32BnjzzcFf\nu1Yt6NnTlkz++Wfw10+011+3EtXHHx91JCIi0VGiEJLWre0T+7BhUUfyt19/tXguuAD22COcNi64\nwBKRp58O5/qJNH68rXZo0CDqSEREoqNEISTO2Zh9MpV1vv12m8h4zTXhtbHXXnDccfDAA6m9m+ba\ntdYbpN4EEcl0ShRCdOKJyVPWedEim7x42WVQt264bfXubZMlP/ww3HbCNHmy9YyccELUkYiIREuJ\nQog2lHV+4QUr6xylW2+1zZ769Qu/raOPtqGNxx8Pv62wjBtnKx1atIg6EhGRaClRCNk550CdOnDn\nndHF8M038J//wIABpW8fHZRKlWxVxbPP2j4JqcZ7m59w/PHa20FERIlCyGrWtI2iHn0UliyJJoZB\ngyxZufjixLV59tlW1fCllxLXZlBmzYLvvtP8BBERUKKQEL1729cHH0x823PmwJNPWu2EeLaQrqi9\n97aVH6NGJa7NoIwfD9tuC0ccEXUkIiLRU6KQADvuCOeeC/ffn/iu+IEDbaz9/PMT2y5Yr8Lbb8OC\nBYlvuyLGjYPOnaFatagjERGJnhKFBLn8cht6SOQn7I8+ghdfhJtuiuZN77TTrN3RoxPfdnn9+CMU\nFGi1g4jIBkoUEqRRIzj1VNut8a+/EtPmNdfA/vtD9+6JaW9z2dk2zv/MM9G0Xx4vv2yrVTp3jjoS\nEZHkoEQhgQYOtA2Tnngi/LYmTrTjllvsjS8qp58On34KX3wRXQzxGD8eDj8cateOOhIRkeSgRCGB\nmjWzXoVbbw23V8F7uPpqKyEddRd65842ifK556KNoyxWrIBJk7TaQURkY0oUEuz668PvVRg7Fj7+\nGAYPjr4OQM2a9sb77LPRxlEWb75pCZwSBRGRvylRSLANvQq33AJr1gR//b/+st6Eo4+GI48M/vrl\ncdpp8NlntptmMhs/Hg44wParEBERo0QhAjfdZLPr77sv+Gvfey989RXccUfw1y6vTp2sLkEyDz+s\nWwevvhr9UI2ISLJRohCB/faDXr2sVyHIao0//WRJyEUX2SfjZFGjhr0BJ/Pww5Qptg23hh1ERDal\nRCEiN95okw4HDQrumldfbXULgrxmUE47zapEfv551JGUbNw42HVXyM2NOhIRkeSiRCEiO+1kyyUf\nfBBmzKj49d5913ZrvPXW5Fza17EjbLcdjBkTdSRb2rAJVJcutqGViIj8TX8WI3TZZdC8OfTsCWvX\nlv86f/xhJaIPPzyaUs1lUb26deu/8ELUkWxpzhyb16H5CSIiW1KiEKEqVWDkSJg9G26/vfzXufJK\n+OUX26EymT8Rd+1qQw/Jtvph/Hir9dC+fdSRiIgknyR+W8kMOTkwYADccAO8/378z3/mGRgxAoYP\nh4YNg48vSB07wjbbJF+vwrhxtjKjRo2oIxERST6hJArOuT2ccyOdc18751Y65750zt3onKsaRnup\nbtAgOOwwK3f8ww9lf94nn0CPHraXQ69e4cUXlJo14bjj4Pnno47kbz/9BNOna7WDiEgsYfUo7Ac4\n4HygKdAXuAC4NaT2UlqVKrZ0sFo1OOaYsi2ZnDfPPqHvvz/85z/RV2Asq65dYeZM+PLLqCMxL79s\nwzX//GfUkYiIJKdQEgXv/Zve+x7e+4ne+2+9968Aw4GTw2gvHdSrB2+9ZUlC69alv5G+9x60aQO7\n7AKvvw5ZWYmLs6I6d7Z4k2X4Yfx4u5c77BB1JCIiySmRcxS2B35NYHspp3FjmDrVegdatoShQ2H5\n8r9/vmAB9OkDRx0FLVpYwpBqb3BZWXDssckx/PDnnzBhglY7iIiUJiGJgnNub6APMCIR7aWyvfe2\nMfPzz4frrrN6C/vvb0lE/frw5JOWQLz+OtSpE3W05dO1KxQUwDffRBvH22/D6tWanyAiUhrnvS/7\nyc4NAQaUcooHmnjvv9joObsB7wKTvPelTrlzzuUABe3atSM7O3uTn+Xl5ZGXl1fmWNPBggU2I3/e\nPJvH0KKFffrd7NaknBUrYOed4eaboV+/6OI491z46KPkrRYpIhJLfn4++fn5mzy2fPlyJk+eDJDr\nvS8Mqq14E4UdgK11dn/tvV9XfP6uwDvAVO/9uWW4fg5QUFBQQE5OTpnjktRz0kmwaBFMmxZN++vX\n2xyPnj1hyJBoYhARCVJhYSG5Voc+0EShSjwne++XAkvLcm5xT8Ik4GPgvPhDk3TWtSuceSZ8/z00\naJD49j/80CaOathBRKR0YdVR2BUbbvgO6A/s7Jyr65yrG0Z7knqOO86Wg774YjTtv/SS9Sgcckg0\n7YuIpIqwJjMeDTQEjgJ+ABYCi4q/ipCdbXUgolj94D2MHWvzPZK55LWISDIIq47CKO995c2OSt77\nymG0J6mpa1eYMsUmbSbS7Nnw9dc2T0JEREqnz1MSmeOPt9UcY8cmtt2XXrItr7UJlIjI1ilRkMjU\nrg0dOiR++GHsWCv6VK1aYtsVEUlFShQkUl27wuTJsHhxYtr77jvbTEvDDiIiZaNEQSK1YUJhooYf\nxo2znoROnRLTnohIqlOiIJHacUebK5Co4YexY22vjO22S0x7IiKpTomCRK5rV3j3Xfjll3Db+fln\nG+bQsIOISNkpUZDInXii1TYYNy7cdp5/3oY5TtZm5yIiZaZEQSJXty60axf+8EN+vhV5SrWtuUVE\noqREQZJC164wcSL8+ms41//hB/jgAzjjjHCuLyKSrpQoSFI46STb0XH8+HCuP2YMVK9uqyxERKTs\nlChIUth1V2jdGp59Npzr5+fDP/+p1Q4iIvFSoiBJo3t3eOstWBjw1mFffgkFBZCXF+x1RUQygRIF\nSRp5eVYMadSoYK87apT1JBx7bLDXFRHJBEoUJGlkZ9ukxkcfteWSQVi3Dh57zHorsrKCuaaISCZR\noiBJpUcP+L//sxUKQXjjDRvK6NkzmOuJiGQaJQqSVNq1g4YN4eGHg7neyJHQsiXk5ARzPRGRTKNE\nQZJKpUrQu7etfqjopMaFC+GVV9SbICJSEUoUJOn06AE1asADD1TsOvffb/MSuncPJi4RkUykREGS\nTna2JQsjRsDKleW7xp9/2vN79rTriYhI+ShRkKR0ySWwbBk8/nj5nv/oo7B8OVx6aaBhiYhkHCUK\nkpQaNrR9GQYPhtWr43vuqlUwdKjVZdhjj3DiExHJFEoUJGkNGgQ//QQPPRTf8x58EBYvhhtvDCUs\nEZGMokRByM/PjzqEEu2zD5x3HtxyC/zyS9mes2QJDBliz9t773Djq4hkvefpTPc88XTP00NoiYJz\nbpxz7jvn3Crn3ELn3BPOuXphtSfll8z/M99yi329/PKynX/55VBUBDffHF5MQUjme56udM8TT/c8\nPYTZozAJOBVoDJwMNAKeC7E9SUM77wzDh8Po0fDii6Wf+8or8OSTcMcdULduYuITEUl3VcK6sPf+\nno2+/cE5NxQY65yr7L1fH1a7kn7OOQdeew3OOgsaNYLmzbc8Z84cOPNM6NLFzhcRkWAkZI6Cc64O\n0B2YoiRB4uWcLZNs3Bg6dNhyH4ipU+HII6FBA+tRcC6SMEVE0lJoPQoAxb0IfYAs4EPguK08pQbA\n3LlzwwxLNrN8+XIKCwujDmOrhg+Hfv2gbVvbE2LffeHLL+G99+CAA+DOO+Grr6KOsmxS5Z6nE93z\nxNM9T6yN3jtrBHld5+PYz9c5NwQYUMopHmjivf+i+Pw6QB1gD+AG4HfvfcxkwTnXDXiqzAGJiIjI\n5rp7758O6mLxJgo7ADts5bSvvffrSnjubsAPwKHe++mlXP8Y4FsgzjI7IiIiGa0GsCfwpvd+aVAX\njStRqFBDzjXAEoAjvPeTE9KoiIiIVEgoiYJz7mDgIOADYBmwN3ATsBPQzHu/NvBGRUREJHBhrXpY\nidVOmADMA/4LfIr1JihJEBERSREJG3oQERGR1KO9HkRERCQmJQoiIiISU8ITBedcb+fcN8WbRU1z\nzh20lfOPcM4VOOdWO+e+cM6dnahY00U899w5d5Jz7i3n3M/OueXOuanOuY6JjDcdxPt7vtHzWjvn\n1jrnVKUmTuX421LNOXerc+7b4r8vXzvnzklQuGmhHPe8u3PuU+fcn8WbBT5SXG9HysA519Y5N945\nt8A5V+ScO74Mz6nwe2hCEwXn3OnAHVjxpZbATOBN59yOMc7fE3gFmAg0B+4BRjrnjk5EvOkg3nsO\ntAPeAjoDOcA7wMvOuRJ2WJCSlOOeb3heNjAKmwQscSjnPX8OaA+ci21elwfMDznUtFGOv+etsd/v\n/wJNga7AwcB/EhJweqiFLQy4CCtwWKrA3kO99wk7gGnAPRt974Afgf4xzr8NmLXZY/nAa4mMO5WP\neO95jGvMBq6L+rWkylHee178uz0I+8NbGPXrSKWjHH9bOgG/AttHHXuqHuW451cAX272WB/g+6hf\nSyoeQBFw/FbOCeQ9NGE9Cs65qkAultkA4C3qCcChMZ7Wii0/Xb1ZyvmykXLe882v4YBtsT+qshXl\nvefOuXOBvbBEQeJQznveBZgBDHDO/eicm++cG+acC7RGfroq5z3/EKjvnOtcfI26wKnAq+FGm9EC\neQ9N5NDDjkBlYPFmjy8GdonxnF1inL+dc656sOGlpfLc881diXV3jQkwrnQW9z13zu0DDMbqsxeF\nG15aKs/veUOgLbA/cCJwKdYV/kBIMaabuO+5934qcCbwrHPuL2ARVpCvT4hxZrpA3kO16kFiKt6k\nayBwqvd+SdTxpCPnXCVsI7QbvPcb9r7URtnhq4R13Xbz3s/w3r8BXA6crQ8h4XDONcXGyG/E5j8d\ng/WiPRxhWFIGoW4zvZklwHqg7maP1wV+ivGcn2Kc/7v3fk2w4aWl8txzAJxzZ2CTjLp6798JJ7y0\nFO893xY4EGjhnNvwabYSNurzF9DRe/9uSLGmi/L8ni8CFnjv/9josblYkrY7kCIblkemPPf8KmCK\n9/7O4u9nO+cuAt53zl3rvd/8k69UXCDvoQnrUfBWurkAOGrDY8Xj30cBU2M87cONzy/Wsfhx2Ypy\n3nOcc3nAI8AZxZ+0pIzKcc9/B5oBLbBZyc2BEVjp8+ZAiTutyt/K+Xs+BdjVOZe10WP7Yr0MP4YU\natoo5z3PAjbfWbgIm72vXrRwBPMemuBZmqdh+0CcBeyHdTktBXYq/vkQYNRG5+8JrMBmbu6LLQn5\nC+gQ9YzTVDnKcc+7Fd/jC7DMc8OxXdSvJVWOeO95Cc/XqoeQ7zk27+Y74FmgCbYseD4wIurXkipH\nOe752cCa4r8tewGtgY+AqVG/llQ5in9vm2MfLIqAy4q/rx/jngfyHhrFC70I2256FZbVHLjRzx4D\nJm12fjssc10FfAn8K+p/rFQ74rnnWN2E9SUcj0b9OlLpiPf3fLPnKlFIwD3Haie8CfxRnDTcDlSP\n+nWk0lGOe94b+Kz4nv+I1VWoF/XrSJUDOLw4QSjx73NY76HaFEpERERi0qoHERERiUmJgoiIiMSk\nREFERERiUqIgIiIiMSlREBERkZiUKIiIiEhMShREREQkJiUKIiIiEpMSBREREYlJiYKIiIjEpERB\nREREYvp/Hm3Z115INRkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f32855a2a58>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create an example dataset\n",
    "\n",
    "pm=[0.15,0.3,0.55,0.65,0.8]\n",
    "v=0.04\n",
    "ph=[-0.2,-0.3,0.5,0.4,0.2]\n",
    "x=np.linspace(0,1,1000)\n",
    "y=np.zeros(x.shape)\n",
    "for m,h in zip(pm,ph):\n",
    "    y+=norm.pdf(x,loc=m,scale=v)*h\n",
    "\n",
    "plt.plot(x,y,'b')\n",
    "plt.plot([0,1],[0,0],'k--')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Build neighbor structure\n",
    "\n",
    "neigh=np.ndarray(len(x),dtype=object)\n",
    "for nb in np.arange(0,len(x)):\n",
    "    if nb>1 and nb<len(x)-1:\n",
    "        neigh[nb]=[nb-1,nb+1]\n",
    "    elif nb==1:\n",
    "        neigh[nb]=nb+1\n",
    "    else:\n",
    "        neigh[nb]=nb-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data=y\n",
    "ind=np.argsort(data)\n",
    "clust=np.zeros(data.shape[0],dtype=int)\n",
    "nclust=0;\n",
    "start=0\n",
    "stop=len(data)-1\n",
    "\n",
    "while start-stop<1:\n",
    "\n",
    "    # Retrieve cluster assignment of neighbors\n",
    "    start_neigh_clust=clust[neigh[ind[start]]]\n",
    "    stop_neigh_clust=clust[neigh[ind[stop]]]\n",
    "\n",
    "    # How many clusters is it touching?        \n",
    "    start_neigh_unique=np.unique(start_neigh_clust)\n",
    "    start_neigh_valid=start_neigh_unique[start_neigh_unique>0] # Remove border and unassigned vertices\n",
    "    stop_neigh_unique=np.unique(stop_neigh_clust)\n",
    "    stop_neigh_valid=stop_neigh_unique[stop_neigh_unique>0] # Remove border and unassigned vertices\n",
    "\n",
    "    # Is this a new cluster?\n",
    "    if len(start_neigh_valid)==0:\n",
    "        nclust+=1\n",
    "        clust[ind[start]]=nclust\n",
    "        start+=1\n",
    "    if len(stop_neigh_valid)==0:\n",
    "        nclust+=1\n",
    "        clust[ind[stop]]=nclust\n",
    "        stop-=1\n",
    "\n",
    "    # No, is it a border vertice (does it touch more than one cluster or only a border)?\n",
    "    if len(start_neigh_valid) > 1:\n",
    "        clust[ind[start]]=-1\n",
    "        start+=1\n",
    "    if len(stop_neigh_valid) > 1:\n",
    "        clust[ind[stop]]=-1\n",
    "        stop-=1\n",
    "\n",
    "    # No, then assign it to the unique cluster it touches\n",
    "    if len(start_neigh_valid)==1 and len(stop_neigh_valid)==1:\n",
    "        start_dist=np.abs(data[ind[start]]-np.mean(data[clust==start_neigh_valid[0]]))\n",
    "        stop_dist=np.abs(data[ind[stop]]-np.mean(data[clust==stop_neigh_valid[0]]))\n",
    "        if start_dist>stop_dist:\n",
    "            clust[ind[start]]=start_neigh_valid[0]\n",
    "            start+=1\n",
    "        elif stop_dist>start_dist:\n",
    "            clust[ind[stop]]=stop_neigh_valid[0]\n",
    "            stop-=1\n",
    "        else:\n",
    "            clust[ind[start]]=start_neigh_valid[0]\n",
    "            clust[ind[stop]]=stop_neigh_valid[0]\n",
    "            start+=1\n",
    "            stop-=1\n",
    "    else:\n",
    "        if len(start_neigh_valid)==1:\n",
    "            clust[ind[start]]=start_neigh_valid[0]\n",
    "            start+=1\n",
    "        if len(stop_neigh_valid)==1:\n",
    "            clust[ind[stop]]=stop_neigh_valid[0]\n",
    "            stop-=1\n",
    "\n",
    "if np.sum(clust==0)>0:\n",
    "    raise ValueError(\"Unassigned vertice detected.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgoAAAFkCAYAAABB1xPiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XmcjeX/x/HXdc4ZQppCSdqzxK+ioYRQqfiW7KpBKdpX\nrdKCKFqVdtU36StTKe3KUnYJM5ZRQ0q0IVtD1jnnXL8/zlTIMMu5z3XmnPfz8bgfzZy5z3W9z904\n5zP3fd3XZay1iIiIiOyNz3UAERERiV8qFERERKRAKhRERESkQCoUREREpEAqFERERKRAKhRERESk\nQCoUREREpEAqFERERKRAKhRERESkQCoUREREpECeFwrGmCOMMf8zxqwzxmw1xiw0xqR53a+IiIiU\nXMDLxo0xBwMzgS+AVsA6oCaw0ct+RUREJDqMl4tCGWMeARpba1t41omIiIh4xutLDxcB84wx7xhj\n1hhjsowxV3ncp4iIiESJ12cUtgEWeBJ4FzgdGAZca6393172r0zkEsUKYLtnwURERBLPAcCxwHhr\n7fpoNep1obADmGOtbbbLY8OAhtbapnvZvyvwpmeBREREEl83a+3oaDXm6WBGYBWQs8djOUDHAvZf\nATBq1Cjq1KnjYSzZ1W233cZTTz3lOkZS0TGPPR3z2NMxj62cnBy6d+8O+Z+l0eJ1oTATqL3HY7WB\nlQXsvx2gTp06pKXpDspYSU1N1fGOMR3z2NMxjz0dc2eieune68GMTwFnGGP6GmNOyL+0cBXwnMf9\nioiISBR4WihYa+cBHYB0IBu4D7jVWvuWl/2KiIhIdHh96QFr7ThgnNf9iIiISPRprQchPT3ddYSk\no2Meezrmsadjnhg8vT2yqPLXgMjMzMzUABgREZEiyMrKokGDBgANrLVZ0WpXZxRERESkQCoURERE\npEAqFERERKRAKhRERESkQCoUREREpEAqFERERKRAKhRERESkQCoUREREpEAqFERERKRAKhRERESk\nQCoUREREpEAqFERERKRAKhRERESkQCoUREREpEAqFERERKRAKhRERESkQCoUREREpEAqFERERKRA\nKhRERESkQCoURCRu5OXlkZ2dzQ8//IC11nUcEUGFgojEiffee48jjzyWU045hRo1alC/fkMyMzNd\nxxJJep4WCsaY/saY8B7bt172KSKlz6uvvkrnzp35/fdGwDTgQ775Bpo0OZMJEya4jieS1GJxRmEx\nUBU4PH87MwZ9ikgp8fXXX3PttdcB1wPvAc2AtoRCM8nLa0m7dh1ZvHix25AiSSwWhULQWrvWWvt7\n/rYhBn2KSCmQl5fH5Zf3xJg0YBhgdvnpAVj7Nnl5x9Op0yVs377dUUqR5BaLQqGmMeZXY8wPxphR\nxpijYtCniJQCL7/8MsuWLSEUegVI2cseFQiFRrNs2TIGDx4c63gigveFwmzgCqAVcB1wHDDNGFPB\n435FJM7t2LGDBx98GGu7A/X2sedJWNuHIUMeZfny5bGKJyL5PC0UrLXjrbXvWWsXW2snAhcAhwAX\ne9mviMS/UaNGsXbtauDeQuzdFziUvn3v8ziViOzJxPpeZWPMHGCitfZf/+JN5EJlZvPmzUlNTd3t\nZ+np6aSnp8copYh4yVrLySefyrffHom1nxTyWa8A17Bw4UJOOeUUL+OJxL2MjAwyMjJ2eyw3N5dp\n06YBNLDWZkWrr5gWCsaYA4GfgH7W2uf28vM0IDMzM5O0tLSY5RKR2JozZw6NGjUCxgH/KeSz8ggE\natO2bQPee2+Mh+lESqesrCwaNGgAUS4UvJ5H4XFjTHNjzDHGmCbA+0AekLGfp4pIAhsx4nUCgerA\n+UV4VgrB4D28//57LF261KtoIrIHrwczHgmMBpYAbwFrgTOstes97ldE4tTOnTsZPfptgsFugL+I\nz+6Bz1eVxx9/wotoIrIXXg9mTLfWHmmtLWetPdpa29Va+6OXfYpIfPviiy/YtGkDUJwxR2UJhW7h\njTf+x++//x7taCKyF1rrQURiasyYdwkEarLvWyL35RpCIR/Dhw+PZiwRKYAKBRGJmVAoxAcffEww\n2IHdZ2EsisqEw5fx7LMvsnPnzmjGE5G9UKEgIjEze/ZsNm5cC7QrYUs3s3btKsaOHRuNWCKyDyoU\nRCRmxo0bRyBQBWhUwpZOwudrwTPPPB+NWCKyDyoURCRmPvroM4LBVhT9bod/C4dv5KuvZpCdnV3y\nYCJSIBUKIhITa9asYfHi+UDrKLXYnkCgGs8//0KU2hORvVGhICIxMWnSpPyvzotSiykEg1fzxhuj\n2LRpU5TaFJE9qVAQkZiYNOkLAoGTgapRbPVqtm/fxptvvhnFNkVkVyoURCQmJkz4kmCwZZRbPRK4\niGeffZFYL3AnkixUKIiI51asWMFvv60Ezo5629ZeT05ONjNmzIh62yKiQkFEYmDKlClEJlhq5kHr\n5xII1NSgRhGPqFAQEc9NnTqVQOAU4BAPWvcRDN7Au+++y6pVqzxoXyS5qVAQEc9NnjyDYNCLswl/\nuQJry2r9BxEPqFAQEU+tXr2alSu/B870sJeDCYd78OyzL7Jjxw4P+xFJPioURMRTM2fOzP/Ky0IB\n4BY2bPid0aNHe9yPSHJRoSAinpo1axaBwNFAdY97qo0xbXj00Sd1q6RIFKlQEBFPTZ/+FcFgk5j0\nZe3dLF36DZ9++mlM+hNJBioURMQzO3bsYP78TKBxjHo8E5+vCQ8++LDOKohEiQoFEfHMggULCAZ3\nErtCwRAOP8C8ebOZMGFCjPoUSWwB1wFEJHF9/fXX+HxlCYfrxbDXVvj9TejT5z7OO+88fL7o/D1k\nrSUrK4sZM2bwxx9/cPzxx3PBBRdQuXLlqLQvEq9UKIiIZ77++muMORUoE8NeDaHQEBYubMFbb71F\n165dS9zihAkTuPPOe8jOno/PdwA+38EEg2soW/YA+vd/gLvvvhu/3x+F7CLxR5ceRMQzM2fOIRQ6\nzUHPzTGmPbfffjebN28udiubN2/msst60KpVK7755kDgM8LhTQSDq4BV7NhxI/feex+XXJJOXl5e\n1NKLxBMVCiLiiY0bN+ZPtHS6k/6tfYq1azdy7733Fev53377LfXrNyQj431gBOHwVKA1kJK/R1Xg\nceB9xo79gBtvvCkquUXijQoFEfHEvHnz8r9ycUYB4FjC4Yd57rlnizyw8aOPPqJhw0asXJlCKDQP\nuILIolZ70w5rX+KVV15m1KhRJcwsEn9iVigYY+4xxoSNMUNj1aeIuDNv3jz8/oOAmg5T3ILPdz6X\nXNKV5cuX73fvcDjMoEGDaN++Pdu3n0coNBuoVYh+egLduOGGm1mzZk1JQ4vElZgUCsaY04BrgIWx\n6E9E3Js7dx7WNsDtiUsf4fBoNm+uxFlnncuPP/5Y4J6rVq2iVasL6NevP9YOwNp3gQOL0NfTbN3q\nK/alDpF45fm/YGPMgcAo4CrgD6/7E5H4MHv2PMLhBq5jAJUJhSaxapWfU089jdGjRxMKhf7+aW5u\nLo8//ji1atVh8uQFwOdAP4r+9liFUKg/I0aMICcnJ4r5RdyKRan/PPCxtfbLGPQlInFg3bp1rFr1\nExAPhQLA0QSDs9m06Sy6devG4YcfSbNmzTn11NM49NCq9OlzL3/+2ZVQ6Fvg/BL0cy1+f3UGDnwo\nWsFFnPO0UDDGXArUB/p62Y+IxJesrKz8r9Kc5thd5fzLCXNYt+50ZsyYzoIF88jLuwprVwAvAJVK\n2EdZgsE+vP32W/u8zCFSmng24ZIx5kjgaeBca22RbjC+7bbbSE1N3e2x9PR00tPTo5hQRLySlZWF\n31+RUKiG6yh7cRrwMPBR/ve3EN2VLa/E5+vHM888y1NPaey2eCMjI4OMjIzdHsvNzfWkL+PVwinG\nmHbAWCDEP/cV+QGb/1hZu0fnxpg0IDMzM5O0tHj6S0REiqJLl4sZO3Y14fA011EKsBg4Of/rpRTu\nzoaiuIcDD3yJNWt+o3z58lFuW2TvsrKyaNCgAUADa23W/vYvLC8vPUwi8i+xPlAvf5tHZGBjvT2L\nBBFJHHPnziccTuZi/xr+/DOXMWPGuA4iUmKeFQrW2i3W2m933YAtwHprrYYEiySoTZs25c/ImMyF\nwvH4fOcyfPh/XQcRKbFY3+CsswgiCW7hwr+mS6nvNIdr4fAVfPXVdA1qlFIvpoWCtfYca+3tsexT\nRGJr/vz5+HxlgTquozjWHp+vAqNHj3YdRKREtNaDiETV/PkL8Pn+j38WT0pWFQiH2/O//2Xsf1eR\nOKZCQUSiat68BQSDp7qOEScuZunSb/jmm29cBxEpNhUKIhI1eXl5LFnyDck+PuEfrfD7KzJ27FjX\nQUSKTYWCiERNTk4OweBOIndDC5QlFGrDO++oUJDSS4WCiETNP3c8nOI0R3xpz+LFC1i5cqXrICLF\nokJBRKJm4cKFpKQcC6Tub9ck0hpjUvj4449dBxEpFhUKIhI18+cvJC9P4xN2dxDGtOCDD1QoSOmk\nQkFEombBgkVofMK/hcNtmDp1Cn/++afrKCJFpkJBRKJi9erVbNjwOxqfsDcXEAzuZPLkya6DiBSZ\nCgURiYp/BjLqjMK/1SQQOIHPPvvMdRCRIlOhICJRsWjRIvz+A4HjXEeJS8FgKz75ZLzrGCJFpkJB\nRKJi0aJFRFaW19vK3rXi55+X88MPP7gOIlIk+hctIlGRmbmIUEjjEwp2Fsb4mThxousgIkWiQkFE\nSmznzp0sW5ZD5IyC7N1B+HyNmDhxkusgIkWiQkFESmzp0qUEg3nojod9C4XOZdKkyYTDYddRRApN\nhYKIlFh2dnb+VzqjsG8t2bRpAwsWLHAdRKTQVCiISIktWrSIlJSjgINdR4lzjfD5yvHFF1+4DiJS\naCoURKTEFixYRDCoyw77VxZoypdfTnEdRKTQVCiISIktWJCNtbrsUBjh8NlMmzadYDDoOopIoahQ\nEJES2bhxI2vW/IIGMhbW2Wzdupn58+e7DiJSKCoURKRENJCxqBri85VnypQproOIFIoKBREpkezs\nbIxJAWq7jlJKpABNmDx5qusgIoWiQkFESiSyxkMdIh+AUhjhcAumTZtOKBRyHUVkvzwtFIwx1xlj\nFhpjcvO3WcaY1l72KSKxNX9+NsGgLjsUTXO2bNmUvz6GSHzz+ozCz0AfIA1oAHwJfGiMqeNxvyIS\nA+FwmMWLs9H4hKI6HWPKMG3aNNdBRPbL00LBWvuptfZza+0P1trvrbX3A38CZ3jZr4jExsqVK9m2\n7U90x0NRHYDPdzrTpk13HURkv2I2RsEY4zPGXAqUB76KVb8i4h3d8VB8oVAzpk6dgbXWdRSRffK8\nUDDGnGSM2QzsAF4AOlhrl3jdr4h4Lzs7G7//EKC66yilUDPWr1/D999/7zqIyD7F4ozCEqAecDrw\nIvCGMebEGPQrIh5btOivGRmN6yilUBPAMGPGDNdBRPYp4HUH1togsDz/2/nGmNOBW4HrC3rObbfd\nRmpq6m6Ppaenk56e7llOESm6rKxFhMPnuI5RSqUSCJzCtGnTufLKK12HkVImIyODjIyM3R7Lzc31\npC8T6+tjxpgvgJXW2p57+VkakJmZmUlaWlpMc4lI0ezYsYPy5SsQDj8PXOs6ThEt5p9xFUuBWo5y\n3Myxx47nxx+/c9S/JJKsrCwaNGgA0MBamxWtdr2eR2GwMaaZMeaY/LEKQ4AWwCgv+xUR7+Xk5BAO\nh9BAxpJoyooVy1izZo3rICIF8nqMwmHASCLjFCYRmUvhfGvtlx73KyIe+2eyIBUKxXcmALNmzXKc\nQ6RgXs+jcJW19nhrbTlr7eHWWhUJIgkiOzublJRjgYquo5RiRxIIHK0BjRLXtNaDiBTLwoXZ5OXp\nbEJJBYNNmTp1pusYIgVSoSAixTJ//iI0I2M0NGXhwiy2bdvmOojIXqlQEJEiW79+PevWrULjE6Kh\nKcFgHnPnznUdRGSvVCiISJH9M5CxntMcieFk/P6KzJypyw8Sn1QoiEiRLVq0CJ+vLFDDdZQE4Mfa\nM5g+XYWCxCcVCiJSZJFCoS4xmNw1KYTDTZk5cxbhcNh1FJF/UaEgIkWWmbmIYFCXHaKnKZs2bWTp\n0qWug4j8iwoFESmSUChETs5iND4hmhoBPo1TkLikQkFEimTZsmXs3Lkd3fEQTRUJBE5RoSBxSYWC\niBSJ7njwRjDYlClTVChI/FGhICJFsnDhQgKBakAV11ESTGSBqN9//911EJHdqFAQkSKZP38hoZDO\nJkRfU0ALREn8UaEgIkUyf/4irFWhEH1HEwgcpXEKEndUKIhIoW3YsIHVq39G4xO8oQWiJB6pUBCR\nQtNARq81Zf78eVogSuKKCgURKbQFCxbkT91cy3WUBHUmwWAe8+bNcx1E5G8qFESk0BYuXIjPdzKa\nutkrkQWiZsyY4TqIyN9UKIhIoc2du0BTN3vKTzjcmGnTVChI/FChICKFsnPnTpYu/QaNT/CWtWcy\nY8ZMLRAlcUOFgogUSk5ODsFgHnCq6ygJrhl//pnL4sWLXQcRAVQoiEghLViwIP8rnVHwViOMSWH6\n9Omug4gAKhREpJDmz59PSkpNoKLrKAmuHD5fQ6ZOneY6iAigQkFECmnevPnk5dV3HSMphELNmTx5\nOtZa11FEVCiIyP6Fw+H8Sw8anxAbzVi3bhU//PCD6yAi3hYKxpi+xpg5xphNxpg1xpj3jTGaqUWk\nlPnxxx/ZsmUTkOY6SpJoChimTdPlB3HP6zMKzYBngUbAuUAKMMEYU87jfkUkirKysvK/0hmF2DiY\nQKAeU6dOdR1ExNvp1ay1F+z6vTHmCuB3oAGgGUVESomsrCwCgeoEg4e5jpI0gsGzmDTpfdcxRGI+\nRuFgwAIbYtyviJTAvHlZBIM6mxBbLfjtt5WsWLHCdRBJcjGbsN0YY4CngRnW2m9j1a+IlIy1lrlz\nM4GbXUdJMs0Bw9SpUzn22GNdhym0YDDI119/TU5ODn6/n3r16nHqqacS+QiQ0iiWK7u8ANQlMkpn\nn2677TZSU1N3eyw9PZ309HSPoolIQX766Sdyc9cTuWIosVOJQKAekydPpkePHq7D7Fc4HOaVV15h\nwICHWL36l91+VrNmHR5/fAjt2rVzlC7xZGRkkJGRsdtjubm53nRmrfV8A54DVgJH72e/NMBmZmZa\nEYkP7733ngUs/GbBJsiWnf+asLA0DvIUtPW21aodbcPhsOtfg31at26dPfvsc/OPZ3cLsy3stLDN\nwkRrzPkWsFdddbXdsWOH67gJKzMz86/f6zQbxc9wz8coGGOeA9oBZ1trf/K6PxGJrnnz5hEIHAFU\ncx0lCZ3DqlU/xfV8CqtWraJx42ZMm7YAmAT8j8iNbinAAcC5WPs58F9ee20kbdt2YMeOHQ4TS1F5\nPY/CC0A3oCuwxRhTNX87wMt+RSR65s7NJBjUZQc3mgM+vvzyS9dB9io3N5eWLVuxfHkuodBMoGUB\nexqgJ+HwJ0yc+AXdu1+u1TFLEa/PKFwHHARMAX7bZbvY435FJAqstcyZMw84zXWUJJWK3386kyZ9\n4TrIv4RCITp3voTvvvuZUGgiUJi59M4jHB7Nu+++w+DBg72OKFHiaaFgrfVZa/172d7wsl8RiY7l\ny5ezadMGoKHrKEkrFGrJxIlfxt1f4IMGPcQXX0wkFHqHyDj1wuoI9Kdfv36aUKqU0FoPIlKguXPn\n5n+lMwrunMcff6xj4cKFroP8bdasWQwcOBBr+wHnFaOFBzCmGd269WDz5s3RjidRpkJBRAo0Z84c\nUlKOB6q4jpLEGuPzVWDChAmugwCwdetWunXrgTGNgPuK2YqfcPh1Vq1aR9++90YznnhAhYKIFGjW\nrDnk5elsgltlsPZsxo0b7zoIAP37D+Cnn34mHB5ByabiOY5weBAvvPA88+bNi1Y88YAKBRHZq7y8\nPBYsyCJyq5u4ZG0rZs2a4fw0fXZ2NkOHDiUc7gfUjkKLN+PzncQNN9z811w6EodUKIjIXmVnZ7Nj\nxzbgdNdRhNYEg3lOb5O01nLddTdiTE3gzii1GiAUGsbcubN55513otSmRJsKBRHZq6+//hpjAkQm\nTBW3ahAI1ODTT8c5S/D2228za9Z0QqFngDJRbPlsjGnLXXf1ZefOnVFsV6JFhYKI7NXs2bPx++sB\n5VxHESAYvJCPPvrUySn6bdu2cccdfTCmLcW7y2HfrB3ML7+s5OWXX45621JyKhREZK+mTfuKYLCx\n6xjytzasWfMrCxYsiHnPTz/9NKtW/Ya1T3jUw/9hbXcefPBhtm3b5lEfUlwqFETkX9atW8eKFcsA\nFQrxozl+/0F8/PHHMe117dq1PPzwI1h7A1DTw576sX79Wl566SUP+5DiUKEgIv/y1Vdf5X+131Xh\nJWbKEApdwLvvfhDTXgcNeojt2wEe8LinE7D2cgYPfkxnFeKMCgUR+ZeZM2fmrxh5tOsospv2ZGfP\nZ8WKFTHpbfny5bz44ouEQn2IzaRb97J+/e+89tprMehLCkuFgoj8y9SpMwkGmxJZ9U/ixwX4fGUZ\nO3ZsTHq7774HsLYK0Dsm/UENrL2UwYMfIy8vL0Z9yv6oUBCR3ezYsYPMzLnoskM8qoi1rRk92vs5\nBxYsWMBbb40mFOoPlPe8v3/05bfffmL06NEx7FP2RYWCiOxm7ty55OXtAJq5jiJ7Ye3FZGZ+7fnl\nh7vv7ovfXxPo6Wk//3YSxlzI4MGPxd2KmclKhYKI7Gb69On4/RWBU1xHkb1qi89XjoyMDM96mDJl\nChMnfk4oNBhI8ayfgljbh++++5bPPvss5n3Lv6lQEJHdTJkyjXC4KSVb8Ee8cyDhcHtGjPifJ5Mv\nWWu5/fa78ftPAzpFvf3CORO/vxGPPPK4o/5lVyoURORvwWCQGTNmYm1z11Fkny5n2bIc5s6dG/WW\n33nnHebPn0so9BjuBrMaQqG7mDFjqievUYpGhYKI/C0rK4utWzcDZ7mOIvt0HoHAkbz66n+j2uqO\nHTu4666+GNMG978D7QkETuDJJ4c6ziEqFETkb1OmTMHnqwA0cB1F9slPMNiTUaNGR3Xp6WeffZZf\nfvkJax+LWpvF5ycY7M2YMWNYuXKl6zBJTYWCiPxt4sQvsbYZ0V0dULxxNdu3b+ONN96ISmtr1qxh\nwIBBWHsdUCcqbZbclRhzEE8/Pcx1kKSmQkFEANi5cyczZkzH2rNdR5FCORJrO/Dkk8Oichthnz73\nsH17AHiw5NGipgKh0LW8/PKrbNq0yXWYpKVCQUSAyLLS27dvBVq6jiKFdgc//riMDz4o2foP06dP\nZ+TI1/Nvh6wcnWhRcxPbtm3jv/+N7ngMKTwVCiICwMSJEwkEKgOnuo4ihXYGPt9Z9O8/qNhnFXbs\n2EHPntfg958BXB3VdNFRHWvTefLJYQSDQddhkpKnhYIxppkx5iNjzK/GmLAxpq2X/YlI8Y0bN4Fg\nsCX6+6F0CYcfZPHiBYwZM6ZYz+/ffwA//PADodArxO//+9v59deVMVvjQnbn9W9FBWABcAMQ/ZlB\nRCQqNmzYwPz5c4FWrqNIkTXHmDbceec9bN26tUjPnDp1Ko899hjWPgic5E28qKiPz3cOjz76pCeT\nTMm+eVooWGs/t9b2s9Z+iJahE4lbEydOzH8DPt91FCkGa4fy22+/MWBA4Qcirl69mi5d0jGmOXC3\nd+GiJBy+g6ysOcyYMcN1lKQTr+eZRCSGPv10HIHAycCRrqNIsdQkHB7AE088wZQpU/a795YtW7jw\nwnZs2ADh8GjA73nCkmtNIFCXxx57wnWQpKNCQSTJhUIhPvnkM4LBC11HkRK5G2Na0KFDF7777rsC\n9/rzzz9p06YdCxZ8Qyj0MVAtdhFLxEcweAeffvoxS5cudR0mqahQkP3buhWysiL/ldiJ0XGfO3cu\nGzeuBZKrUChH5P6OcmxzHSVK/ITDY9i8+VCaNm3B7NmzIw/v8nu0dOlSGjduxrRpXxMOj6P0zcDZ\nDZ+vKk888aTrIPFn61bIyfGk6bhcHu62224jNTV1t8fS09NJT093lCjJLVkCDRpAZiakpblOkzxi\ndNw/+ugj/P7KhEKNPesjHp0IZAFpLGc+9VzHiZLKhEKT2bixPU2aNKVr166k167Fhf36MeCii3j4\ns/FYewzh8ExK5zLiZQmFbuX11/szaNBADj/8cNeBnMnIyNh9qfHcXHKnTfOmM2ttTDYgDLTdzz5p\ngM3MzLQSRzIzrYXIfyV2YnTca9asa6GHBZtEW7Y9Nf+bUxkbB3mive20MMwGArV3eZ1YeMTCljjI\nV5Jto/X7K9o+ffp4+u+i1MnMtJlgiWxp1kbv89vreRQqGGPqGWPq5z90fP73R3nZr4gUztKlS1m2\n7Fugg+soElUpwC0Eg0t2eewRoA9Q3k2kqDmYUOh6nnvuRf744w/XYZKC12MUGgLzgUwiVc6TRM72\nxdNk4iJJ67333stfLVK3RUpp0ptt23bw/PPPuw6SFLyeR2GqtdZnrfXvsfX0sl8RKZyMjDGEwxcS\nGdonUlpUIxzuxeOPP8WWLVtch0l4uutBJEktXbqUxYsXABe7jiJSDHezaVMuL730kusgCU+FgkiS\nysjIwO+vCFzgOopIMRyDtVcwePBjOqvgMRUKIknIWsvIkW8SCnVGlx2k9LqXjRs38MILL7gOktBU\nKIgkoa+++ooVK74HuruOIlICx2FtTwYPfpTNmze7DpOwVCiIJKHXXhtBIHAMcJbrKCIl9ACbNv3J\n0KFDXQdJWCoURJLM5s2bGT36LYLBK9BbgJR+RxIO38Sjjz7B77//7jpMQtK7hEiSGT16NNu3bwV0\nl7Ikir7s3Okv0jLbUngqFESSiLWWYcOeBy4CjnYdRyRKKhMK3c/w4cP59ttvXYdJOCoURJLIlClT\nyMnJxtqbXUcRibKb8fmO5eabe2OtdR0moahQEEkijz76BH7/KcA5rqOIRFlZgsGn+PLLiYwdO9Z1\nmISiQkEkSSxatIjx48cRCt0JGNdxRDxwEca04cYbb9XtklGkQkEkSQwc+BCBwLHApa6jiHjG2udY\nu3Yj99wHkXz3AAAfa0lEQVTT13WUhKFCQSQJLFy4kPfeG0MweB+RJYhFEtUxhMODeeGF55k6darr\nMAlBhYJIErj77r4EAjWBHq6jiMTAzfh8zejWrQebNm1yHabUU6EgkuA+//xzJkz4jGBwCDqbIMnB\nRzg8ktWrN3LttdfpLogSUqEgksC2b9/OjTfeis/XAujoOo5IDB1HKDSct97KYPjw4a7DlGoqFEQS\n2EMPPcyPP/5IOPwCutNBks+lwA3cfPMtzJw503WYUkuFgkiCmjNnDkOGDMHa+4C6ruOIOPIU4XBj\nLrqoPT/88IPrMKWSCgWRBJSbm0uXLukYkwbc6zqOiENlCIffY9OmQ2jZshWrVq1yHajUUaEgkmBC\noRBdu3bn11/XEwq9hQYwilQhFJrAL79sp3nzc1QsFJEKBZEEYq3l1lt789ln4/KLhONdRxKJE8cS\nCk1mxYrNNG7cTJchikCFgkiCsNZy77338fzzz2HtC0Br15FE4kxNgsEZ/PKLn4YNz9CETIWkQkEk\nAQSDQW666WYeeWQI8CRwretIInHqWEKhWWzadArnnNOShx56iGAw6DpUXFOhIFLKrV27ltatL+TF\nF18CXgZudx1JJM5VJhweTzjcl379+nP66U3IyspyHSpueV4oGGNuNMb8aIzZZoyZbYw5zes+RZLF\nhx9+SN26pzBlShbWjgeudh1JpJQIAIOwdiaLFm2jYcOGXHbZ5Sxbtsx1sLjjaaFgjLmEyHnQ/sCp\nwEJgvDGmipf9iiS6xYsX85//XEj79u1Zv74BodAioKXrWCKl0BmEQvOx9nneemsStWvXpm3b9owb\nN06XJPJ5fUbhNmC4tfYNa+0S4DpgK9DT435FEo61lhkzZtChQydOPvlkJk1aAryLtR8D1VzHEynF\nAsD1BIM/YO1wPvvsRy688EKqVDmcK664goyMDH799VfXIZ0JeNWwMSYFaAAM/usxa601xkwCGnvV\nr0iiWbt2LaOeeoqXX36NJUsWEwjUAl4hGLwcKOM6nkgCKQdcTTB4FTCf3NwxvPnmJ4wcORKAww6r\nTsOG9albtw41atTg2GOPpXr16lStWpVKlSrh9/udpveKZ4UCUAXwA2v2eHwNUHtfT2zQoMFu3xtj\nKFeuHOXLl6dixYqkpqZyyCGHUKVKFSpXrozPpzGZXkoDegGPPPIIPx96qOs4SWHt2rWUy8lhJNCq\ndWvm//2TEwgGzwHm529SPIt3+foRYJKrIDH0JPCT6xClTjDYHKgFzOT3339l3LhfGTfuU9ex/mXK\n0KFU9KhtLwuFqLHWsnXrVrZu3cq6desK9Ry//2h8vsoeJ0t8weAP1Leb6AWMGfM52Sk1XEdKWNZu\nIRj87u/vT93lZ0cddSypqYfkf/d1THMlonA4TKMDm8CcWRx/3EbyKnzlOpJn6oZOhpxsjj/u4IR+\nnd47In+LyMvLY82a3/jjjw3uIu3ilmHDONijtr0sFNYBIaDqHo9XBVbv64mZmZmkpaXt9lgwGGTL\nli1s3ryZ3NxcNm7cyNq1a1mzZg0///wzK1euJCdnGUuW5LB162ZCoZ8wpgLB4PnARUALSkldFGe+\nAU4CwPI6eXlaqjj6lgDP4vf/D4AmTZrRvXs6l9SsCeedR1ZmJuzx70GiICsLGjTg3XffSuzjmyyv\nM9llZZHVoAEN9r9nkXn2yWmtzTPGZBIZiv0RgDHG5H//TFHbCwQCpKamkpqaypFHHlngfuFwmOXL\nlzN79mymTp3Kp5+OZdWqYQQCVfOv6V4PHFes1yQSXdkY0x9r36dy5apcf/2t9OzZk+OOy//91H3d\nIhIHvP4Teyjwen7BMIfIXRDlgde96tDn81GjRg1q1KhB9+7dsdaSmZnJqFGjGDHiVTZvfhJIx9oB\ngE6jiwvrgb7Aqxx11HH06/cq3bt3p2zZsq6DiYj8i6ejAK217wB3AgOJjLw6BWhlrV3rZb+7MsbQ\nsGFDnn76aVat+oVnnhnGYYdNwZi6QB8id2uKxMrH+P11OfDAdxg27Gm+/34JvXr1UpEgInHL89sF\nrLUvWGuPtdaWs9Y2ttbO87rPgpQvX56bbrqJH39cxoAB95OSMoxAoB6Rkx0iXgoCdwFtadXqdL77\nLodbbrmFlBQtAS0i8S0p7yssV64c/fr1Y/HiRdSrVwljmgLPA9Z1NElIW/D52uHzPcVTTz3FJ598\nRLVqmiBJREqHpCwU/lKrVi2++moGN998A3ATcDORGzVEomUTfv/5lC07jc8+G0fv3r2JjOkVESkd\nkrpQAEhJSWHYsGEMHz4cY17EmO5EThOLlNRW/P42lCv3DZMnT+L88893HUhEpMg0sUC+a665hipV\nqnDxxZcQDoO1o4hMLClSHCGM6UZKSiYTJ35Bo0aNXAcSESmWpD+jsKuOHTvy9ttvAe8QuQyhMQtS\nXA8AHzFmzNucccYZrsOIiBSbCoU9dOrUiZdfHg68SGRudJGi+ggYwqOPPkKbNm1chxERKREVCntx\n1VVXcc899wB3A+Ndx5FS5Vf8/iu56KJ23Hnnna7DiIiUmAqFAjz00EO0atUav78r8LPrOFIqWHy+\nXlSufAAjRvxXdzeISEJQoVAAv9/P6NGjOOywCvnFgm6blP0ZSTg8nhEjXqFyZa1cKiKJQYXCPlSq\nVIm3336TcHgW8JjrOBLX1uP330nXrt244IILXIcREYkaFQr70axZM/r0uRtj+gPZruNI3LqfcuWC\nDB2qAbAiklhUKBTCgAEDqFmzJn7/VegShPzbYox5mYED+1O1alXXYUREokqFQiGULVuWESNeIRye\nC7zgOo7EGZ+vD0cffRw33nij6ygiIlGnQqGQmjRpwlVXXY3ffz+wxnUciRvTCYfH8dhjgylTpozr\nMCIiUadCoQiGDBnMgQcGgL6uo0ic8Pv783//V4/OnTu7jiIi4gkVCkVQuXJlBg8eBLwOZDlOI+5N\nJxSazKBB/fH59E9JRBKT3t2K6JprrqFWrTr4/XegtSCSm883mDp1TqZdu3auo4iIeEaFQhEFAgGG\nDn2MUGgK8JnrOOLMQsLhz7nvvj46myAiCU3vcMVwwQUX0LRpc/z+vkDYdRxxYihHHHE0F198sesg\nIiKeUqFQDMYYHntsCKHQIiJLUktyWYUxGdx++y2kpKS4DiMi4ikVCsXUpEkTWrX6D4HAADQJU7J5\nibJly9CrVy/XQUREPKdCoQQeemggweBS4C3XUSRm8ggEXuaKKy7j4IMPdh1GRMRzKhRKoGHDhvzn\nPxcSCDxE4p1VsMA44LpdHnsOyHETJ258SDC4mhtuuMF1EBGRmPCsUDDG3GuMmWmM2WKM2eBVP64N\nGNCPYHAJ8J7rKFG0FL+/CXAhJ574B+efdx4AB6cuBP4P6A3scJjPHZ9vOI0aNeHkk092HUVEJCa8\nPKOQQmSk34se9uHc6aefTsuW5+P3DyYx5lWYgM/XkGOO2cCXX35JTk42jzzyCADjx4/jiSceJxB4\nEb+/NbDZbdSYW044PInrr7/GdRARkZjxrFCw1j5orR1GEqzNfP/9fQmFFhI5VV+afY7PdxGtWrVg\n4cJMzj777N1+mpKSwh133MHkyV9wwAFZ+HwdgTw3UZ14nfLlK9KlSxfXQUREYkZjFKKgRYsWnHba\nGfj9j7iOUgLz8fk60br1+XzwwVgOPPDAAvc888wz+eSTD/H5pgJ3xi6iU2ECgdfp1u1Sypcv7zqM\niEjMqFCIAmMM9913D6HQDGCW6zjFsJ5AoD316tVlzJi3C7UK4llnncVTTw0FngE+9Tyhe5MJBn/m\nyiuvdB1ERCSmilQoGGOGGGPC+9hCxphaXoWNZxdddBE1apyIMY+5jlJEFmN6UaHCn3zwwXtF+mv5\nxhtvzJ9L4hpgk3cR48IbHHdcTc444wzXQUREYipQxP2fAEbsZ5/lxczyt9tuu43U1NTdHktPTyc9\nPb2kTXvG5/Nxzz13ctVVVwNLgdquIxXSCKz9kJEjP+Doo48u0jONMbz88kvUrl2XYHAAMNSThO5t\nxe8fyxVX3IUxxnUYEREyMjLIyMj454HcXHK96sxa6+kG9AA2FHLfNMBmZmba0mj79u22SpXDLVxj\nwZaC7Vfr96fayy/vse8XtmWLtZmZkf/uxeDBg60xAQvL4uA1ebGNtoD9/vvvo/9LU4LjLiWULMc3\nWV5nstuyxWaOGmWJ3H6XZqP4Oe7lPApHGWPqAccAfmNMvfytgld9ula2bFl6974Jn+8NYK3rOPtl\nTG9SU8vmjzXYh/LlIS0t8t+96N27N1WrHo4x93uQ0j1jMmjYsBEnnHBCbDvez3GXEkqW45ssrzPZ\nlS8Pdep40rSXgxkHAllAf+DA/K+zgAYe9unc9ddfT5kyPuAF11H2YxLWjmHYsCepVKlSiVoqV64c\nAwf2w9q3gcXRiRc3NgKf0717/F72EhHxkpfzKFxprfXvZZvmVZ/xoFKlSvTqdSWBwPPANtdxChAk\nELiVJk2a0a1bt6i02KNHD6pXPwZjBkelvfjxPhDUctIikrR0e6QHevfuTSi0DhjlOkoBhhMK5fDc\nc8OiNjivTJky9O17F/A2URjPGjd8vnc488wWVKtWzXUUEREnVCh4oEaNGrRp05ZA4Ckg7DrOHjbh\n9w/g8st7cOqpp0a15SuvvJKDDjoEGBbVdt1Zj7WTuPRSzcQoIslLhYJH7rrrDoLBHGC86yh7eBy/\n/08eemhQ1FsuX748N910HX7/ayTGvAofAmE6derkOoiIiDMqFDxy5plncuqpp+HzxdPcAqvx+YZy\n2223cuSRR3rSw/XXX4+124A3PGk/lny+d2natDlVq1Z1HUVExBkVCh4xxnD33bcTDk8CFrqOk+9h\nKlQoQ58+fTzroXr16rRv355A4EVK92qaf2DtJC6+WGcTRCS5qVDwUKdOnahW7SjiY8bCFRgznHvu\nuYtDDjnE056uv/46gsFvgZme9uOtT7E2jw4dOrgOIiLilAoFD6WkpHDnnb0xJgP41XGaQRxyyCHc\ncsstnvd0zjnncNRRxwH/9bwvrxgzlrS00z27RCMiUlqoUPDYVVddRYUK5YmssujKMowZyQMP9N3n\n8tHR4vP5uOaanvj9Y4DNnvcXfdsw5nO6dNHZBBERFQoeO+igg7jhhmvx+1/C1Z0AxjzIYYcdznXX\nXRezPnv06EE4vBV4N2Z9Rs8kwuGttG/f3nUQERHnVCjEwC233IIx24DhDnr/BmtH07//fRxwwAEx\n6/Woo46iRYtz8PlGxqzP6PmAE06ozYknnug6iIiIcyoUYqB69epcdtll+RMw7Yhp38b058gjj6FX\nr14x7RfgiisuIxyeCvwU876LL0Qg8DGdO+tsgogIqFCImT597iYUWg3E8i/sLKx9j4ED+1GmTJkY\n9hvRsWNHypYtB4yOed/FN5tgcC1t27Z1HUREJC6oUIiR2rVr06lTZwKBR4FgTPr0+R7g+ONrcdll\nl8Wkvz1VrFiRdu3aEghkOOm/eD6iUqXDaNSokesgIiJxQYVCDN1//30Eg8uJzV/YMwiHxzF48EAC\ngUAM+tu7rl3TCQYXAd86y1AUgcBHtG/fBr/f7zqKiEhcUKEQQ/Xq1ctfLGoQ3p5VsPj993DSSfXp\n0sXtgkatW7fmwANTgbec5iicZQSDS7joootcBxERiRsqFGLswQf7Ewx+j7dLUH9EKDSTJ554BJ/P\n7f/ismXL0qlTewKBt4n/KZ0/JiWlLOedd57rICIicUOFQoylpaXRvn1HAoEH8eYOiDwCgT6cffa5\nnH/++R60X3SXXHIJweB3QLbrKPvk939My5YtqVChgusoIiJxQ4WCA4MHP0Qo9BPwogetv0Qo9B1D\nhz6OMcaD9ouuZcuWHHTQIcDbrqPsw0bC4em0a6fLDiIiu1Kh4ECdOnXo2bMnfv8gYGMUW16H39+f\nnj17Ub9+/Si2WzJlypShY8d2BALvEr+XH8ZjbYg2bdq4DiIiEldUKDgyaNBAypTZCQyIYqt9KVcu\nzODBD0exzejo0qVL/uWHb1xHKcDH/N//1dMiUCIie1Ch4Ei1atUYMOABjHkeWBiFFmcAr/Loo4M5\n7LDDotBedEWu/R9EfK79EMTv/4wOHXTZQURkTyoUHOrduze1ap2I3381JbtdciuBQC9OP71xTBd+\nKoqyZcvSvv1FBAJjXUfZi68IhTbqtkgRkb1QoeBQmTJlGDHiFcLhTODxErR0Nz7fT4wc+Zrz2yH3\npWPHjgSD2cAy11H28AmVKh1Gw4YNXQcREYk78fupkiQaN25Mnz53Y0w/YHYxWngPeJ6hQ5+I+9UO\nW7dunb/2Q3ydVQgEPqFt2wviusgSEXHFk3dGY8wxxphXjTHLjTFbjTHLjDEDjDEpXvRX2g0cOJDT\nTjuNQKAz8FsRnrkQn68HnTp14YYbbvAqXtSUL1+eCy74D35/PBUKywkGv9VlBxGRAnj1J9SJgAGu\nBuoCtwHXAfE3HD8OpKSk8P7771KlisHv/w+wvhDP+g6/vzUnn1ybkSNHxM2cCfvTsWMHQqE5wC+u\no+T7lECgjGZjFBEpgCeFgrV2vLW2l7X2C2vtCmvtJ8ATQEcv+ksERxxxBF98MZ6DDvoNv78F8OM+\n9p6F39+M448/hAkTPitVMwm2adMGvz8AfOg6CgA+38e0aNGCihUruo4iIhKXYnlR9mBgQwz7K3Xq\n1q3LrFnTqV59G35/feApYPMue6wC7sSY5px2Wk1mzZoWl7dC7svBBx/M2Wefg88XD5cfNmHtFM3G\nKCKyDzEpFIwxNYCbgJdi0V9pduKJJ7JoURa9eqXj892Fz1eFlJR6pKTUAapTrtxLDBr0INOmTaZK\nlSqu4xZL584dsXYqhbvE4qXxWJun2RhFRPbBWFv4KXWNMUOAPvvYxQJ1rLXf7fKc6sAU4Etr7bX7\naT8NyGzevDmpqam7/Sw9PZ309PRCZ00EP//8Mx988AFLliwhJSWF+vXr0759ew4++GDX0Upk9erV\nHHHEEVj7GnCFwyQ9OPHELHJy4nuxKhGRPWVkZJCRkbHbY7m5uUybNg2ggbU2K1p9FbVQqAxU3s9u\ny621wfz9jwAmA7OstVcWov00IDMzM5O0tLRC55LSp3HjM/n660pY+5GjBCH8/qrcddfVDBkyxFEG\nEZHoycrKokGDBhDlQiFQlJ2ttesp5Pni/DMJXwJzgZ5FjyaJrHPnDsyZcx/W/gkc6CDBV4RC62nb\ntq2DvkVESg+v5lE4gsjlhpXA3cBhxpiqxpiqXvQnpU/Hjh0Jh3cAnzlK8BGVKh1Go0aNHPUvIlI6\neDWY8TzgeKAl8DORWYRWUbTZhCSBHXfccZx0Un1czdIYCHxIhw4XaTZGEZH98GoehZHWWv8em89a\n6/eiPymdLr64I37/J8D2GPe8hGDwO9q1axfjfkVESh/9OSXOdOrUiVDoT2BSjHv+gAMOKM+5554b\n435FREofFQriTN26dalR40Tg3Zj26/d/QOvWrShXrlxM+xURKY1UKIhTl17aGb//Q2BnjHr8lVDo\nazp27BCj/kRESjcVCuJU586dCYX+IHInbSx8gN8f0GyMIiKFpEJBnDrllFM47riawJiY9OfzvcfZ\nZ5/DIYccEpP+RERKOxUK4pQxhq5dL8bvH4v3lx/WEg5PpUuXTh73IyKSOFQoiHMXX3xx/uWHCR73\n9D7GQIcOGp8gIlJYKhTEuZNPPpmaNesAb3naj8/3Di1anM2hhx7qaT8iIolEhYI4Z4zh8su74vN9\nAGz1qJc1WDuZrl0v8ah9EZHEpEJB4kJ6ejrh8BbAq9Ukx+Dz+ejYsaNH7YuIJCYVChIXTjjhBE47\n7Qx8vv950r7f/yatWrWmcuX9rZIuIiK7UqEgcePKKy/H2vHA6ii3/AOh0Gy6d+8a5XZFRBKfCgWJ\nG5dccgl+vx94M8otv0H58hW1CJSISDGoUJC4UalSJTp0aE8g8Bpgo9RqmEBgJOnpF1O+fPkotSki\nkjxUKEhcufrqqwgGvwW+ilKLXxIMrqRnz55Rak9EJLmoUJC40rJlS44++njgpai0Z8xwatWqS+PG\njaPSnohIslGhIHHF5/Nx443X4vO9A6wtYWurgA+46abrMMZEIZ2ISPJRoSBxp1evXgQCBhhewpZe\nomzZslx++eXRiCUikpRUKEjcqVy5MldccTmBwHPA9mK2sg2//wWuuupKUlNToxlPRCSpqFCQuHTH\nHXcQCv0OjCxmCyMIhzfQu3fvaMYSEUk6KhQkLtWqVYvOnbsQCAyh6MtP7yAQeIRLLrmEE044wYt4\nIiJJQ4WCxK0BA/oTCv0MvFzEZ75MKPQr/fo94EUsEZGkokJByMjIcB1hr+rWrUuPHj3w+wcAGwr5\nrA0EAgPp0aMHderU8TBdycTrMU9kOuaxp2OeGDwrFIwxHxpjVhpjthljfjPGvGGMqeZVf1J88fyP\neciQwRxwQB7Qp5DPuJcyZXYwePDDXsYqsXg+5olKxzz2dMwTg5dnFL4EugC1gI7ACcAYD/uTBHT4\n4Yfz5JOPAa8Cn+xn73HAcB5//BGqVVNNKiISDZ4VCtbaYdbaOdban621s4FHgDOMMX6v+pTEdM01\n19CmTVv8/u7AtwXstQS/vzutW1/AddddF8t4IiIJLSZjFIwxlYBuwExrbSgWfUriMMYwatQb1K59\nDIHAOcCsPfaYg99/NiecUI2MjDfx+TT0RkQkWgJeNm6MeQS4CShPZJWfNvt5ygEAOTk5XsaSPeTm\n5pKVleU6xn4988xQeve+ncWLmwJnA7WB74EvqFWrLs888zTLly93G7KQSssxTyQ65rGnYx5bu3x2\nHhDNdo21hV/O1xgzhH2PKrNAHWvtd/n7VwIqAccA/YFN1toCiwVjTFfgzUIHEhERkT11s9aOjlZj\nRS0UKgOV97PbcmttcC/PrQ78DDS21n69j/ZbASso/ty9IiIiyegA4FhgvLV2fbQaLVKhUKKOjDma\nSAFwlrV2Wkw6FRERkRLxpFAwxpwOnAbMADYCNYCBwKHASdbavKh3KiIiIlHn1fDwrUTmTpgELAFe\nARYQOZugIkFERKSUiNmlBxERESl9dMO5iIiIFEiFgoiIiBQo5oWCMeZGY8yP+YtFzTbGnLaf/c8y\nxmQaY7YbY74zxvSIVdZEUZRjbozpYIyZYIz53RiTa4yZZYw5P5Z5E0FRf893eV5TY0yeMUaz1BRR\nMd5byhhjHjbGrMh/f1lujLkiRnETQjGOeTdjzAJjzJb8xQL/mz/fjhSCMaaZMeYjY8yvxpiwMaZt\nIZ5T4s/QmBYKxphLgCeJTL50KrAQGG+MqVLA/scSWQnoC6AeMAx41RhzXizyJoKiHnOgOTAB+A+Q\nBkwGPjbG1ItB3IRQjGP+1/NSgZFEBgFLERTzmI8hMsXnlUQWr0sHlnocNWEU4/28KZHf71eAukBn\n4HTg5ZgETgwViNwYcAORCQ73KWqfodbamG3AbGDYLt8b4Bfg7gL2fxRYtMdjGcC4WOYuzVtRj3kB\nbSwG7nf9WkrLVtxjnv+7/SCRN94s16+jNG3FeG9pDWwADnadvbRuxTjmdwDL9njsJuAn16+lNG5A\nGGi7n32i8hkaszMKxpgUoAGRygYAG0k9CWhcwNPO4N9/XY3fx/6yi2Ie8z3bMEBFIm+qsh/FPebG\nmCuB44gUClIExTzmFwHzgD7GmF+MMUuNMY8bY6I6R36iKuYx/wo4yhjzn/w2qgJdgE+9TZvUovIZ\nGstLD1UAP7Bmj8fXAIcX8JzDC9j/IGNM2ejGS0jFOeZ7uovI6a53opgrkRX5mBtjagKDiczPHvY2\nXkIqzu/58UAz4P+A9sCtRE6FP+9RxkRT5GNurZ0FdAfeNsbsBFYRmZDvJg9zJruofIbqrgcpUP4i\nXQ8AXay161znSUTGGB+RhdD6W2t/+Othh5GShY/Iqduu1tp51trPgduBHvojxBvGmLpErpEPIDL+\nqRWRs2jDHcaSQvB0mek9rANCQNU9Hq8KrC7gOasL2H+TtXZHdOMlpOIccwCMMZcSGWTU2Vo72Zt4\nCamox7wi0BCob4z5669ZH5GrPjuB8621UzzKmiiK83u+CvjVWvvnLo/lECnSjgR+2Ouz5C/FOeb3\nADOttUPzv19sjLkBmG6Muc9au+dfvlJyUfkMjdkZBRuZujkTaPnXY/nXv1sCswp42le77p/v/PzH\nZT+KecwxxqQD/wUuzf9LSwqpGMd8E3ASUJ/IqOR6wEtEpj6vB+x1pVX5RzF/z2cCRxhjyu/yWG0i\nZxl+8ShqwijmMS8P7LmycJjI6H2dRfNGdD5DYzxK82Ii60BcDpxI5JTTeuDQ/J8PAUbusv+xwGYi\nIzdrE7klZCdwrusRp6VlK8Yx75p/jK8jUnn+tR3k+rWUlq2ox3wvz9ddDx4fcyLjblYCbwN1iNwW\nvBR4yfVrKS1bMY55D2BH/nvLcUBTYA4wy/VrKS1b/u9tPSJ/WISB3vnfH1XAMY/KZ6iLF3oDkeWm\ntxGpahru8rMRwJd77N+cSOW6DVgGXOb6f1Zp24pyzInMmxDay/aa69dRmrai/p7v8VwVCjE45kTm\nThgP/JlfNDwGlHX9OkrTVoxjfiOQnX/MfyEyr0I116+jtGxAi/wCYa/vz159hmpRKBERESmQ7noQ\nERGRAqlQEBERkQKpUBAREZECqVAQERGRAqlQEBERkQKpUBAREZECqVAQERGRAqlQEBERkQKpUBAR\nEZECqVAQERGRAqlQEBERkQL9P2dgVEWbPQZyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3210154a20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "clid=np.unique(clust)\n",
    "clid=clid[clid!=-1]\n",
    "for cl in clid:\n",
    "    xi=np.where(clust==cl)\n",
    "    xc=np.hstack((x[xi],x[xi][::-1]))\n",
    "    yc=np.hstack((y[xi],np.zeros(len(xi[0]))))\n",
    "    plt.fill(xc,yc,'b')\n",
    "for xi in x[np.where(clust==-1)]:\n",
    "    plt.plot([xi,xi],[-1,1],'r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "How many steps? 1000\n",
      "None\n",
      "> \u001b[0;32m<ipython-input-34-870a3e550981>\u001b[0m(18)\u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m     16 \u001b[0;31m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m     17 \u001b[0;31m        \u001b[0;31m# Retrieve cluster assignment of neighbors\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m---> 18 \u001b[0;31m        \u001b[0mstart_neigh_clust\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mclust\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mneigh\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mind\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mstart\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m     19 \u001b[0;31m        \u001b[0mstop_neigh_clust\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mclust\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mneigh\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mind\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mstop\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m     20 \u001b[0;31m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\n",
      "ipdb> q\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgoAAAFkCAYAAABB1xPiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xd4HOW9/v/3o2pLLrLljotc5W4sGRdcMKGacAIhJxBR\nEuDLLwRIQkwoaZwUEiABBwOHxMkJJ4ABEXwcEiAEDDa4V627LRe59yIhF8lq+/z+GBls47VVZna2\n3K/r2ovVenaeewdp97MzTzHWWkRERETOJsHvACIiIhK5VCiIiIhISCoUREREJCQVCiIiIhKSCgUR\nEREJSYWCiIiIhKRCQUREREJSoSAiIiIhqVAQERGRkFQoiIiISEieFwrGmE7GmKnGmEPGmDJjzEpj\nTI7X7YqIiEjjJXm5c2NMBjAfmAlcBRwCegMlXrYrIiIi7jBeLgpljHkSGGWtvcSzRkRERMQzXl96\n+A9gmTHmTWPMfmNMwBhzl8dtioiIiEu8PqNQDlhgEvB/wHDgWeBua+3Us2yfiXOJYhtwwrNgIiIi\nsacJkAV8YK097NZOvS4UKoAl1tqxpzz2LDDMWjv6LNvfDLzmWSAREZHYd4u19nW3duZpZ0ZgL7D+\njMfWAzeE2H4bwKuvvkq/fv08jCWnmjhxIs8884zfMeKKjnn46ZiHn455eK1fv55bb70Vaj9L3eJ1\noTAfyD7jsWxge4jtTwD069ePnByNoAyXli1b6niHmY55+OmYh5+OuW9cvXTvdWfGZ4CRxpgfG2N6\n1l5auAv4b4/bFRERERd4WihYa5cBXwXygNXAT4H7rbVveNmuiIiIuMPrSw9Ya98D3vO6HREREXGf\n1noQ8vLy/I4Qd3TMw0/HPPx0zGODp8Mj66t2DYiCgoICdYARERGph0AgQG5uLkCutTbg1n51RkFE\nRERCUqEgIiIiIalQEBERkZBUKIiIiEhIKhREREQkJBUKIiIiEpIKBREREQlJhYKIiIiEpEJBRERE\nQlKhICIiIiGpUBAREZGQVCiIiIhISCoUREREJCQVCiIiIhKSCgUREREJSYWCiIiIhKRCQUREREJS\noSAiIiIhqVAQERGRkFQoiEjEKC+v4u23NzJ37naCQet3HBFBhYKIRIhHHllAs2b7ue66Powb143m\nzVfz+uvr/I4lEvc8LRSMMT83xgTPuOkvX0ROc8cds/nd70YSDO4DVgFLKStL5ZZbuvHb3y7zO55I\nXEsKQxtrgMsAU/tzdRjaFJEo8b//u4aXXhoFLABG8/lbxQlgLT/6UR/69t3Eddf19i2jSDwLR6FQ\nba09GIZ2RCTKlJdXcc89yUARMILPiwSAJkBfYC833VTFwYMnaN68iR8xReJaOPoo9DbG7DbGFBlj\nXjXGdAlDmyISBe64Yz6Vlb2BRCD5LFukA1BR0Ydrr10QzmgiUsvrQmERcDtwFfAdoDswxxiT7nG7\nIhLhjh+vYNq0HsASoM85tuwFzGfOnFHMm7c9POFE5DOeFgrW2g+stdOttWustR8C1wCtgBu9bFdE\nIt/EiQsIBrsCHeqw9UXAMW69VYWCSLiFo4/CZ6y1pcaYjThfEUKaOHEiLVu2PO2xvLw88vLyvIwn\nImFirWXq1A7AMmBYHZ6RBgTYvn0c06dv4Gtfy/Y2oEiEy8/PJz8//7THSktLPWnLWBu+SU2MMc2A\nHcB/WWv/+yz/ngMUFBQUkJOTE7ZcIhJer766mttuGwQUALl1fFYVcJAuXTazY8c478KJRKlAIEBu\nbi5ArrU24NZ+vZ5H4SljzDhjTDdjzMXAWzh/7fnneaqIxLCnnjoAHAIurMezkoGt7Nx5MTNnbvEm\nmIh8gdedGTsDrwOFwBvAQWCktfawx+2KSIQ6fryS1asHAetxRjvUxzDgU777XfVVEAkXT/soWGvV\nqUBETvPMMyux9iKgpAHPTgUWUVg4inXrDtK/f1uX04nImbTWg4iE1auvHgN2cu4hkecyCLDcf/8q\n90KJSEgqFEQkbKqqati4sS+wjdNnYayP1kCAWbOyKSurdC2biJydCgURCZuXXlqDtR1xplNpjPYE\ng5159NHFbsQSkXNQoSAiYfPSSyVAMdCvkXvqBazhxRebNj6UiJyTCgURCZvly9viDIKq72iHszlG\naekw/vGPjS7sS0RCUaEgImFRWHiA8vJ+NLxvwplygQP89Kd7XdqfiJyNCgURCYspUwpx3nLOOYN7\nPSQD61i3LpedO4+4tE8ROZMKBREJiw8+qAa2A27OfdAHSOWHP1zp4j5F5FQqFEQkLDZt6gbscnmv\nnYDF/POfXaipCd+6NSLxRIWCiHhu4cKd1NT0BFI82HsKlZVZPPvsGg/2LSIqFETEcy++uA0IAr09\n2PswoJBJkyo82LeIqFAQEc/NmWOBTUCGB3tPAHaxZ08OS5Yc8GD/IvFNhYKIeG7r1gsALz/EhwGl\nTJyoORVE3KZCQUQ8tW7dfqqre+IMZ/RKBhBg4cL+fPqpLkGIuEmFgoh4aurUzbX3sjxuqTvWZjBx\n4lKP2xGJLyoURMRTM2eewLns0MHjlnoAS3n99U5UV2uopIhbVCiIiKcKC1sDW8PUWiqVlT149NGC\nMLUnEvtUKIiIZ44ereDo0b5AZZhaHAIU8NxzzQgGdVZBxA0qFETEM9OmbQKaAplhatEANZSV9eXx\nxwNhalMktqlQEBHPvPfeIaAC9xaCqouLgBX85jctqKoKurZXay1vvLGOr33tE770pVncc89cioqK\nXdu/SKRSoSAingkEEoAivJm6ORQDBDlxojf33uvOCIjf/W4p6elrycvrz9//PoqPPx7MlClj6dUr\nhQkTZlFZWeNKOyKRSIWCiHhm166OwGEfWs4BZvO//5vFnj1HG7yXgweP0rv3xzzyyEWUlwMEcN42\n2wAHgQDvvz+enj0XUVZW5UZwkYijQkFEPLFzZwlVVT2BJJ8S9CAYbMaVVzZsCeqPPtpE587b2bx5\nFLAAGIBTgJycOKotMA5Yxq5dw7noovluhBaJOCoURMQT06ZtxnmL6ehTgi7AItauHcNvflO/jo2/\n+MVCrrgik8rKDGA/cDHOJY2zGQ4sZt268Xzve3MalVgkEoWtUDDG/MgYEzTG/D5cbYqIf2bOLAXK\nga4+prgUWMKjj/bgww93nnfrqqogX/rSbH75yxE4cz9kAN3q0M4YYAkvvNCX9eu1MJXElrAUCsaY\ni4BvAw07BygiUWfVqhRgC/6euEwAemHtASZMMHz00a6QWwYC+2nffhUffzwWmAMMBZrVo62eWJvK\nV7+6vnGRRSKM53/BxphmwKvAXcCnXrcnIpFh797OQCQMH2wNpFNTU84VVzTnnnsWUFHx+SiFrVuP\nMGHCHHJz0ygp6QAUAOOp/9tjJrCKDRtG88EHm8+7tUi0CEep/wLwjrV2VhjaEpEIUFR0mJqaHni7\nYmR9XIDzQb6GKVMupmnTYjIylpGevooePVJ5//0RwDIgFWcehoYaDhTzne/scSGzSGTwtDuyMeYb\nwIU4i8WLSJyYPn0rzgez1wtB1UdrYDSwFmvXU1pqcfpQbMDpY3CpC22kAhvYtm008+btZMyYLi7s\nU8RfnhUKxpjOwGTgcmttvQYYT5w4kZYtW572WF5eHnl5eS4mFBGvfPzxUeAo/nZkDGUAYIFBtT9v\nwN2RGbnAMSZO3MjSpSoUxBv5+fnk5+ef9lhpaaknbRlrvVk4xRhzHfB3oIbPxxUl4vyF1gCp9ozG\njTE5QEFBQQE5OTme5BIR73XpsoBdu9JxFmmKRGs4vVDo4/L+5wH9OHSoKZmZaS7vW+TsAoEAubm5\nALnWWtcWO/Gyj8JHOH+JF+K8WwzBuQj4KjDkzCJBRGLHvn3tie++y12BTH70Iy1MJdHPs0LBWnvc\nWrvu1BtwHDhsrdX4IZEYtW/fEaqrs4icjox+6Aos5803m/odRKTRwj3AWWcRRGLcO+8U4VxlbOt3\nFJ+d4MiRocyfv8PvICKNEtZCwVr7JWvtA+FsU0TCa9aswzjdkLJ8TuK3wcAJfvlLzakg0U1rPYiI\nq1autMBO4vvSA0A6sILZszv7HUSkUVQoiIirduzIwFmCWSCJyso+vPfeJr+DiDSYCgURcU1lZRXH\nj/fEufQgzmCvMiZP3u53EJEGU6EgIq6ZPXszzgyIzf2OEiFSgbXMm9fO7yAiDaZCQURc88EHJ9c4\n0IyEnwtSXj6YJUtCr1wpEslUKIiIaxYvPgEcAVr4HSWC9AOqmTRJox8kOqlQEBHXbNyYBuib8+la\nAKv46KNmfgcRaRAVCiLimsOHOwLH/I4RgY5RXDyAgwd1bCT6qFAQEVds27aPmppeeLx6fZTqBjTl\nhRfW+B1EpN5UKIiIK959dwtOkaAe/l/UDdjLW28d9zuISL2pUBARV8yd+ykQBDr5HSVCbWHdum5+\nhxCpNxUKIuKKVasMsBu9rYSSQnV1L2bP1uRLEl30Fy0irtixoxWauvlcsoFq/vznbX4HEakXFQoi\n0mgVFZWUlfUEqv2OEsFaAOuYPTvF7yAi9aJCQUQabf78IqAtoLkCzq2EPXt6U10d9DuISJ2pUBCR\nRpsxY1/tPS2pfG6tsLYN06cX+h1EpM5UKIhIoy1eXI4z0ZKmbj63bKCC11/f63cQkTpToSAijbZh\nQxNgp98xokAqsIFFi5r6HUSkzlQoiEijHTzYEWcxKDm/Ug4c6EdlpTp+SnRQoSAijbJ/fwnV1b3Q\n20ldtQZa8eabG/wOIlIn+ssWkUb517+2AMlApt9RokQfoJy//e2A30FE6kSFgog0yuzZxbX3uvia\nI3okA4UsXpzudxCROlGhICKNsmKFBXbhfABK3ZRy8GA2lZU1fgcROS9PCwVjzHeMMSuNMaW1twXG\nmKu9bFNEwmvbtgxg33m3k1O1Alry1lsb/Q4icl5en1HYCTwC5AC5wCzgn8aYfh63KyJhEAwGOXKk\nG1Dhd5Qokw1UMm2aCiyJfJ4WCtbaf1lr37fWFllrN1trf4YzK8tIL9sVkfBYtmwn0B5I8ztKlGkC\nbGTxYl2ukcgXtj4KxpgEY8w3cN5RFoarXRHxzvvv766919HXHNGpmD17ehAMWr+DiJyT54WCMWag\nMeYozrnJPwBftdZqonORGLBwYRlwHOesgtRPOsFgJ2bP3uZ3EJFzCscZhUJgCDAc+CPwijGmbxja\nFRGPrV+fDGwHjN9RolBvAPLzd/icQ+TckrxuwFpbDWyp/XG5MWY4cD9wT6jnTJw4kZYtW572WF5e\nHnl5eZ7lFJH627u3HVB83u3kbFoARcyerSWnpf7y8/PJz88/7bHS0lJP2jLWhvf6mDFmJrDdWnvn\nWf4tBygoKCggJycnrLlEpH6OHaugeXOApcAYn9PU1xpgUO39DTizJfphHsnJ7ais9Kt9iSWBQIDc\n3FyAXGttwK39ej2PwuPGmLHGmG61fRWeAC4BXvWyXRHx3vvvb8VZDTHD7yhRLImqqj6sX6/pnCVy\ned1HoR3wMk4/hY9w5lK40lo7y+N2RcRjM2cerL3Xzdcc0c05dlOnFvmcQyQ0r+dRuMta28Na29Ra\n28FaqyJBJEYEAlXAbqC531GiWEdgLx9+WO53EJGQtNaDiDRIUVFznEJBGmcb69e38TuESEgqFESk\nQUpKOgNlfseIAVUcP96X4mKdVZDIpEJBROpt48ZigsGOQIrfUWJAOyCF117TAlESmVQoiEi9vfPO\n9tp7HXzNERt6A6W8+26J30FEzkqFgojU27x5R3EuO3T1O0oMSAQ2snx5M7+DiJyVCgURqbc1axJx\nJlz1fHLXOHGUgwd7U1WlWRol8qhQEJF627WrDZq62U0tgZb861/b/A4i8gUqFESkXk6cqOHEiSyg\nxu8oMaQPUMW0aXv9DiLyBSoURKReZszYgTN1cwu/o8SQ5sB6Fi5M9DuIyBeoUBCRevnww/2197L8\njBGDDrJjRxe/Q4h8gQoFEamXZcsqgZ1Apt9RYkwyNTUXEAgc8juIyGlUKIhIvWze3AzY5XeMGNQD\ngNde2+pzDpHTqVAQkXopLtbUzd7oDGxm1qwqv4OInEaFgojUWWFhCcGgM+WweGEHGza09TuEyGlU\nKIhInb3zzo7ae5q62RtByst7snevFoiSyKFCQUTqbO7cI8BRoLvfUWJUJyCBV18t8juIyGdUKIhI\nna1dmwRsRVM3eyUb2M9775X6HUTkMyoURKTOdu/W1M3echaIWrlSk1lJ5FChICJ1cvRoJRUVWYD1\nO0qMK6OkpBfl5VogSiKDCgURqZP33tsKJAOt/Y4S49oATfm//9vidxARQIWCiNTRRx8drL3Xw9cc\nsa8/cIy///2A30FEABUKIlJHS5cGgS04CxiJd5oChSxZkux3EBFAhYKI1NGWLRnAHr9jxIlP2bu3\nJzU16g8i/lOhICLnVVUV5OjRnkCF31HiRHOsbc2HH+44/6YiHvO0UDDG/NgYs8QYc8QYs98Y85Yx\npo+XbYqI+2bN2gWkAxq2Fx7ZQA35+Tv9DiLi+RmFscDzwAjgcpwu0zOMMU09bldEXPTuu3tr72X5\nGSOOZAAbmDfP+B1ExNvp1ay115z6szHmduAAkAvM87JtEXHP4sVVwDZUKITTQbZv74m1YFQviI/C\n3UchA2e2Fk3tJhJFNm1qAezyO0acSaWmpjNz5ui4i7/CNmG7McYAk4F51tp14WpXRBqnpsby6afd\ngcV+R4kzTj+Fl1/ewSWXdPY7TJ1VVlbz17+uZdGiT0lMTODyy9tw4419SUjQaZFoFc6VXf6AM5PI\n6PNtOHHiRFq2bHnaY3l5eeTl5XkUTURCmTNnD3AB0MzvKHGmFbCGTz6Jjg/Ympogd9wxh9df705N\nzZDPHn/xRfjWtzby4x8f4he/uNjHhLElPz+f/Pz80x4rLfVmMTFjrffjdI0x/w38BzDWWhtyvI8x\nJgcoKCgoICcnx/NcInJ+EycuYfLk4cA+oIPfcVyyBhhUe38DEKmDsWaSkDCIqqq2Ef2NfOvWYoYN\n20xx8XBgEc7omN5ADbAOZ7GrIfTv/zFLl44mLS3Fx7SxKxAIkJubC5BrrQ24tV/P+yjUFgnXAZee\nq0gQkci0cOEJYCexUyREk3SCwXbMmBG5b52rVu2nb98DFBf3AZYDI3FOHicDTYAcYDCwgHXrxpKV\ntZxjxzQfRzTxeh6FPwC3ADcDx40x7WtvTbxsV0Tcs3FjC5xCQcKvH1DFK69E5vHfubOU4cOLqazM\nAD4FhobY0gAXA6s4eHAoAwYsprpaq2NGC6/PKHwH5xzUJzhzv5683ehxuyLigmDQUlLSHTjhd5Q4\n1RJYx5w54exOVjdVVTUMHbqBiooLgOPUbehsDrCcHTvG8eUvf+xpPnGPp4WCtTbBWpt4ltsrXrYr\nIu745JMdOB9WWgjKP8Xs2dM74r6BX3PNLA4fHg5sBnrW45kjgLnMmHEJzz+/wptw4iqt9SAiIb31\n1slFoLS0tH9aYW0mf/vbRr+DfObPf17FRx+Nx5k3ryEdzy8G1jFxYgb79x91NZu4T4WCiIQ0b141\nsBXI9DtKHOsPlPLyywf9DgJASUkZ992XjjNaZGQD95IItKKmpgMTJixyL5x4QoWCiIS0aVMmmpHR\nbynAWhYvbu13EACuvXYe1dVdcBYJa0zfiS7AMpYvv5RXX9UcfJFMhYKInNXRo1UcP94LqPI7inCC\nI0f6sn27v6fp3357IwsWXAIsBbq7sMeRwDbuuaeCYND7OX2kYVQoiMhZTZ++BefbbGR8k41v3YFE\nnn9+g28JgkHLbbcdx5l46yKX9poEHOfYsaE89NB8l/YpblOhICJn9a9/HQLKcdYcEH91Bwp5+23/\nzu489NA8jhwZChzGKSDdMgQI8NxznTl2rNLF/YpbVCiIyFktW5aI02Gtqd9RBIDdFBX1oqYm/Kfo\njxwp59lnO+HMvOjF9PotqK7uyp13LvBg39JYKhRE5Kx27uyM8+1RIkMGwWBb3ngj/JcfbrttNjU1\n3fHuMlQvYBHTp/fm00/LPWpDGkqFgoh8wdq1h6mp6Yy7p5ilcQYBpfzP/+wPa6tbtx7i7bdzgGVA\nNw9b6kww2JG77tJZhUijQkFEvuD114tq77nRs13ckQKsZtGiC8La6k03LcOZib+Xxy11BRbz1lt9\ndVYhwqhQEJEvmDmzDNgLhPdDSc6voqIXM2eGZ26LhQu3s3TpWCBAeEa/dCIY7MDddy8JQ1tSVyoU\nROQL1q7NBLbgrPonkWMIcJxJk7aFpbVbby0CKoALw9Kec2ljMdOn96C8XPN3RAoVCiJympKSCo4d\nywY0VC3yNAdWMHt2G89bmj69kC1bLgHWA2met/e5dtTUdOGBB9RXIVKoUBCR00yduhHnenh7v6PI\nWVnKyvoyd663lx++851SYDcw3NN2vqgXEOCvf21LTU1krZgZr1QoiMhp3n67BCgB+vgdRc7qQuAY\nv/rVFs9amDx5BYcOjQD2AMmetRNaMhUV/XnySfVViAQqFETkNMuXN8eZaKkxC/6Id5oBBcyeneXJ\n5Es1NZaf/CQVWAOMcH3/dTMQKOTppxN9al9OpUJBRD5TXl5NcXEfoMzvKHJOTaiq6spf/lLo+p6/\n//3FlJf3w+mj4ldnVgOU8umnF/Hqq2t9yiAnqVAQkc/k52/AWT64rd9R5JyGAdv5/e9LXN1rSUkF\nf/pTV2A+3kzVXB/DgB389KelPucQFQoi8plp0w4Bx9BCUJEuEShi48bB7N7t3tLTN9+8hJqaNkRG\nR9ZEYCs7dgxn3rzwzBshZ6dCQUQ+s2RJC2Atmro5GvQDUrn//hWu7G316oO8//6FwDy8n4WxroYB\nR/je97zruCnnp0JBRAAoLa2kuLgfzhkFiXwdgUX8859ZVFU1fhjh9ddvxumXEK7JleoiHVjOihU5\nbN9+xO8wcUuFgogA8OKLhUAToIPfUaTOMqiu7sJPftK4YYSTJ69iy5ZRwGrCM1VzffQHkrj33lV+\nB4lbKhREBIDp0z8F9uOc0pboMAgo4PnnMxt8VqG0tIKHH24OLAfGuRnOJR2BBbz/fjbHjlX7HSYu\neVooGGPGGmPeNsbsNsYEjTFf8bI9EWm4FSvaAhvR94dok0BFRW8mTmzYWYXLL19KVVVHnPkZIvX/\nfUeCwUx++MPlfgeJS17/VqQDK4B7AfdnBhERV2zeXEJZWTb6M41GQ4GFTJnSmQMH6jf/xdNPr2LZ\nsjE4HRh7exHOJf2ABbz8chtPJpmSc/O0ULDWvm+t/S9r7T/RMnQiEeu55wpx3g4ipbe71E9Hampa\nM2FCoM7PWL78AI880h5YClzqWTL3JFJR0Z2nnlrnd5C4E6nnmUQkjN55B5xVAjv5nEQaJgtYQCBw\nMZMmrTzv1gcOHGf06GKCwWqgK86cBZFuBLCcp57yO0f8UaEgEucqKmrYvr0fzgJAEr0uBZbx4INZ\nvPfetpBb7dlzjD59iigvvwAoJTImV6qLBKCY4uIBvPnmdr/DxBUVCnJ+ZWUQCDj/lfAJ03F/8cVC\nrM0AMj1tJ9I0xbm635Ryv6O4JBHn0tE+rr02nT//eY3z8Cm/R2+/vY0ePfZQWtoV2Iwz9DCajAGK\neOSRQ34HiTxlZbB+vSe7jsjl4SZOnEjLli1PeywvL4+8vDyfEsW5wkLIzYWCAsjxe/73OBKm4/7K\nK4dxhkUO8qyNSNQXCAA5bGE5Q/yO45LWQDXW7uXuuwfy+OMLuSHnAEPe+ju/7/wtVu0aA2zF+f89\n1N+oDZIKbGLbtitZtOggI0fG75ok+fn55Ofnf/5AaSmlc+Z40lZEFgrPPPMMOfpAEgmLQKAzUAhc\n4ncUcUU7oBUwh+3bO/PM9lHAdaTuKsRZ7GkEkOZnwEYaCRzm3ns3EwjEb6HwhS/PgQCB3FxyPWjL\n63kU0o0xQ4wxJ+cE7VH7cxcv2xWRuvn3v7dSVdUDZ0ZGiR3JwHigF4NJ5zhp9OefOP0YorlIAMgA\nlrN8+RCKirSyZDh43UdhGM50XwU4A7Qn4Zzt+6XH7YpIHTz99E6cDm2D/Y4iHkmkjLSY6YdxkvP7\n+u1vn3+EhzSep5cerLWzUYdJkYhkLcyf3xlnTjRddpBo0gGYzccfD2TfvuN06JDud6CYpg9xkTj1\nzjvbqKjoQYR2VRI5j55Y24w776z7JFPSMCoUROLUU0/tAEoAdRyWaNQZZ7Gogezde9zvMDFNhYJI\nHAoGLYsWdcNZVrip33FEGqgn1jbljjt0VsFLKhRE4tAf/7iW6upuOCsGikSrrsB8ZswYxM6dR/0O\nE7NUKIjEoWeeOQpsAS4836YiEa4f1qZy2206q+AVFQoicWbHjqMUFV0IFKG3AIl+nYDFzJ6dy+rV\nmtrZC3qXEIkzDz+8EmekQz+/o4i4ZDBQwY03bvA7SExSoSASR4JByz/+0QFYhtNrXCQWtAZWUFg4\nir//vcjvMDFHhYJIHPnd71ZRUdELSPE7iojLnAWv7rqrlGDQ+h0mpqhQEIkjTz0FsAbNnSCxJxU4\nRElJDg8+uMTvMDFFhYJInHjjjU0UFw8BDgPG7zgiHhgBzOfZZ3uwY4eGS7pFhYJInHjwwcPAduBi\nv6OIeKgbwWAKEyas8jtIzFChIBIHXn11E7t3jwS24SxBLBKrOgMFrFs3msmTtbqkG1QoiMQ4a+H7\n3y8HNqCzCRIfxgMBfvjDdmzbpksQjaVCQSTGPfZYASUlg4FD6GyCxIcEoA3BYCpjxhRqFEQjqVAQ\niWGHD5/gscfaAQF0NkHiS1dgLbt3X8Qtt8zzO0xUU6EgEsO+/OUCqqvbAM3RSAeJP2OBT3jjjZE8\n//xqv8NELRUKIjHqxRfXsnjxSGAe0NvvOCI+uRhYyf33X8DMmdv9DhOVVCiIxKAtW45w992tgHU4\nHbtE4lUK0B1rP+WqqxJYtmy/34GijgoFkRhTWVnDiBFF1NSkAC1QB0aRTCCZmhq4+OJyli9XsVAf\nKhREYkgwaBk6dDGHDg0GtgDd/I4kEiG6ANVUVaVw0UVVfPjhDr8DRQ0VCiIxIhi0jBw5j3XrLgZm\nA8P9jiQfxlTTAAAdOklEQVQSYboDNdTUlHPVVc2YNEkTMtWFCgWRGFBeXs2AAQtZunQsMBP4kt+R\nRCJUFyATa7fw4IOD+NKX5nDiRLXfoSKaCgWRKLdmzSEuuGA1hYXDgE+Ay3xOJBLpWgMXAh/z8cej\nad26iKlTC/0OFbE8LxSMMfcZY7YaY8qNMYuMMRd53aZIPLAWHnlkCYMHG0pKLgBWoxEOInWVhFNU\nr6O8PMg3v9mXbt0W8+9/awjlmTwtFIwxNwGTgJ8DQ4GVwAfGmDZetisS66ZN20ybNiv53e+GY+1G\nwAK5fscSiUKDcOYZ+ZgdOy7gmmu60L59AY8+WkB5uS5JgPdnFCYCf7LWvmKtLQS+A5QBd3rcrkjM\nCQYtzz+/mg4dlnHjjb0oLk4H5gMjgfY+pxOJZknApTjDKOdx4EAqv/51Lmlpx+jRYyH33beAhQv3\nYeN0yYgkr3ZsjEnG+Yrz+MnHrLXWGPMRMMqrdkViTcHqo/zXo/P46KOOVFYOwlkFcjbOn1Evf8OJ\nxJSmwDicM3TrgT1s3dqJP/yhH3/4AyQkHCAzcwddupTTu7dhwICmZGe3JDu7Fb17Z5CWluhvfI94\nVigAbYBE4MyZLfYD2ed6Ym7uH3A6m1QCFUAVKSmW1FRIT6+iRYsaWrc2tG+fTIcOTUlM1Bz2XsoB\n/h/w5JNPsrNtW7/jxIXdu8vYvzSLZGYz9/ZxwHHg78AvgXSck4Fv+hkxyq055f6TwEd+BQmjSYDm\nDmiY40ANwWBfDh4cyMGDfQgEeuLM+niqI8Cx2u3LTrlVADV1uJ16yuLM0xfn/rfrLh5NMo8Av23Y\nSzwHLwuFRlgDZOC8GTq3yso8KivzOHoU9u07ddtKnF/+bbW39cByEhIOkpioGekaq7q6iAvtEf4f\nMG3a+6xO1jdYrwSDNdTU9AVuAa4GErmSj+jGrVRfsJxWrVJ9Thg7gsEgI5pdDEsW0KN7CVXpC/2O\n5Jn+NYNg/Wp6dM+I6dcZHmuB6QBUVgbZt6+GI0dSgHa1t7Y4hXxa7e3k/VSc783JQJPa+6feEmr/\ne9KZX37NWe5/UHtzzFj2FgkcbuwLPCtjPbroUnvpoQz4mrX27VMefwloaa396lmekwMUFBQUkJOT\n89njwaAzTvzQoePs2VPG7t3ObdeuMnburGT7dsuePakcOtScsrJuOP9zAAqB3ThFxxAiti6KaGsZ\nykACQA7TWc4NfgeKQUU4RW4ukEGzZmuYMKGEX33d0vfGS6CgAE75exCXBAKQmxv7xzdeXme8CwQI\n5Oae7NKca60NuLVrzz45rbVVxpgCnPEnbwMYY0ztz8/VZ18JCZCenkR6eku6dWt5zm2rqoLMnLmd\nf/xjFzNnJrBlyyCCwXY4BUMhTu/Wrg14RSJu2wAcBkZiTCtGjVrFb37TnfHjBzr/HHDt71xEpMG8\n/or9e+Cl2oJhCc4oiDTgJa8aTE5O4Oqru3H11c4c9zU1lqlT1/P73x9m9epcnLMNnwBZtTeRcCsG\nVgFjSExszs03L+C55y4iI2Oc38FERL7A0+GR1to3gQeBXwHLgcHAVdbag162e6rERMPtt/dj1aox\n7N+fwvXXLyYhoT/QEWeq27JwRREBFuN0bBrK9dcvpKSkLa+8MoaMDPU/EJHI5PnMjNbaP1hrs6y1\nTa21o6y1y7xuM5R27dJ4660xHDjQnEsuWQyMAfZxeg9oES9U4xSmI2jbdjfLlp3grbfG0ry5OtyK\nSGSLy7UeMjOb8skn43jvvT00bVqJ02/hE7445ETEDceBAmA81147j337csnN1QRJIhId4rJQOGnC\nhO4cONCTwYMX48yRPwtnLKuIW47gdKIdxK9/vZJ33hlDQoLm/RCR6BHXhQJAs2bJrFw5jptumo9T\nLCzAOU0s0lhlwGagF1OmbOWnP9XQNBGJPnFfKJz0xhujeeCBZTjT4i5AZxakcWqAANCfKVN2cffd\nA/wOJCLSICoUTjFp0ggmTiwARuPMpa8+C9JQs4GL+fGP16hIEJGopkLhDL///QhuvnkR8CWcPgsi\n9bUQGM/VV8/j8ceH+R1GRKRRVCicxWuvjWb48Dk4y44u9TuORJW9QG/atVvOu++O9TuMiEijqVAI\nYe7c0WRmrgC64Uz/LHI+FthOQkINc+d216qmIhITVCiEkJKSyPz53UlICAIHUedGOb85wEh+8pNd\n9OnT2u8wIiKuUKFwDtnZrZg8eT/OzNMf+x1HIloxMICuXRfz2GO5focREXGNCoXz+N73hjBs2EKc\nkRAb/I4jEWslkMK77/bwO4iIiKtUKNTBjBnDSE4+CJxAlyDkizYCY/nyl1cyaFBbv8OIiLhKhUId\ntGqVyuTJpcAQYK7fcSTilJCYuI/XXx/hdxAREdepUKije+8dRJ8+83D6KxzwO45EjBXACO67byct\nWqT4HUZExHUqFOrhH//oDyShZanlc5Wkpm5m0iSdTRCR2KRCoR769WvN1762CrgEWO93HPHdCmA4\n991XTFKS/pREJDbp3a2epk4dSUrKNqAcrQUR7ypJSSniySc1TbOIxC4VCvXUtGkSDz9cAuQAS/yO\nI74pBIZz5537SU7Wn5GIxC69wzXAL3+ZS/Pma4F0IOh3HPHFXhIS9vP00xf5HURExFMqFBogIcHw\n+OMWGAgs8DuOhN0+YARXXrmJ9PRkv8OIiHhKhUIDffe7A8nMLAA6okmY4s1qIIHnnx/kdxAREc+p\nUGiESZPSgJ7orEI8qQIG0LdvAb16tfQ7jIiI51QoNMK3vtWPNm0CQDti76yCxemsuZWNzOdSZrGb\nSqDI51x+Wwp04rHH2vkdREQkLDwrFIwxPzHGzDfGHDfGFHvVjt9+/esmQDaw0O8oLtoCBIDhpKT0\nI6vjcdpykBJzGZAFzAIq/Azoo0TS09fyn//Z2+8gIiJh4eUZhWTgTeCPHrbhu7vv7k9Gxiogg9iY\nV2EZ0IqkpEyeeGIlJ070ZM27mbzJTRTP3sqXvzwfGIsz2dBRf6OG3XZgGF//eqnfQUREwsazQsFa\n+0tr7bM4Pb9i2k9/enIERLTPq7AEGEibNlvYurUNP/rREIz5/F+bpSfx7rvjeOaZQqA/sA7nmn28\n2ASc4IknLvQ7iIhI2KiPggseeGAw6elrgVS/ozSC8+Hfps1atmwZROfOzUJu+YMfDOK3v90C5ALz\nwhXQZ0GgN717r6RDhzS/w4iIhI0KBRckJBi+973jwIXAKr/jNEAxkEbTprtYt64fzZuffxXEhx8e\nwg03LAQuBRZ5HTACrAC6MXFiC7+DiIiEVb0KBWPME8aY4DluNcaYPl6FjWS/+tUwkpO3E33X7S2w\nAWOaM2NGc9q2rfu35TffHEPr1iuArsARrwJGiOMkJW3j7rsH+B1ERCSskuq5/dPAX8+zzZYGZvnM\nxIkTadny9DHqeXl55OXlNXbXnklOTuDGG3fx2msjcQ5BD78j1dFsYDwPPbSEMWOG1+uZiYmGd95p\ny+jRGTijPi7zImAEKAOGMGbMchISsvwOIyJCfn4++fn5nz9QWopX3ayNtd721DfGfAt4xlrbug7b\n5gAFBQUF5OTkeJrLCyUlFWRmHsfaFcCX/I5TB/uAVLKyNrB168jQm5WVQWEh9O0LaV8843DZZXOZ\nNWs4sBdn+GSsmQuMZcaMHVxxRdfwNXue4y6NFC/HN15eZ7wrKyPw1lvk3norQK61NuDWrr2cR6GL\nMWYI0A1INMYMqb2le9Wm31q1SmX8+LXACOCQ33HqYDPGVDNjRva5N0tLg5yckG8yb745jISEYmCn\n+xEjQgJpaRvCWyTAeY+7NFK8HN94eZ3xLi0N+vXzZNdedmb8Fc6sPT8HmtXeD+B0lY9Zf/zjQCAR\nWOl3lPNYBozhrrs20bt3q0btKTOzKTfdtAW4GNjoRrgIUgLkcvnlB/wOIiLiCy/nUbjDWpt4ltsc\nr9qMBNnZrejffxkwCCj3O04I1UAzmjVbxx//OMqVPU6ZchEJCYeB3a7sL3KsBprw6KNx2UdXRETD\nI73wzDNdgDZE7gRM84C+vPBCIomJ5rxb10WLFilcd91GYAyww5V9RoamtGixkmHD2vsdRETEFyoU\nPHDlld1o334p0Blnop5IcgQYQI8e8/nmN8/TN6Ge/vSnHIw5SuwsHFUMXMhVV2nKZhGJXyoUPPLo\no01wlqB2reOpS5YALXjtNfeHb7Ztm8bIkWtwJp6KhXkVVgLJPPKINx2ERESigQoFj9xzz2DS0lYD\n55/lMHz2AxeRm7uMkSM7etLCpEl9gJY4nSWjnSE9fSO5uW39DiIi4hsVCh5JSDB8+9ulwGAiZyTA\nKiCR117zbnbBUaM60KFDAdCJ6F5N81PgQsaOLfE7iIiIr1QoeOiJJ0aQmLgV50PHbzuB0Ywfv5Ls\n7AxPW/r+9xOBvkTnuhcnLQMyuP/+LL+DiIj4SoWCh5o0Seb664twpo7Y53OajRhTydSpQzxv6aGH\nLiQxcRfOHATRqoLk5H1cdZVGO4hIfFOh4LEXXhiO84HZ6CUwGmErMJZrrllzzuWj3ZKUlMCYMVuB\noUTfIlngzH8xkNzc3Rh3Ro+KiEQtFQoea9++BaNHL8cZCeDXh+ZOEhI+5eWXh4Wtxcce64XTqbEg\nbG26ZwnQjbvuUidGEREVCmHgTOsM/kzrvBkYzde/XkhmZpOwtTp2bEdatFgDROP88sUYc4zbbgvz\n2g4iIhFIhUIYDBrUkb59F+F08KsIc+uHSEzcy1/+co7VIT1yww1HcPpnRNO0zjVAd7Kzi0iJpJGt\nIiI+UaEQJs89lwW0JrzTOhcCI7n99i00axb+T73HHhsEVBI5w0PrIgBcyE03NfU7iIhIRFChECZX\nXNGDzp0XAVk4izKFwwmSkzfywgsXh6m903Xu3JwuXVYA3kzu5I1dQDXf/W5Pv4OIiEQEFQphNHly\nJtCF8JxVWA1cyPe/f4DU1KQwtHd2N998ck6Fzb5lqA9jMunYcRNt2iT6HUVEJCKoUAijr30tm44d\n5+EUC16eVbBAIk2bruK3v/XnbMJJP/rREOAYsM3XHHWzGWtHMGFCld9BREQihgqFMJs8uTVOobDY\nw1aWAv35+c8rSEz0939xRkYq3buvxnnNkT6l83oglQce6OV3EBGRiKFCIcxuvLE/nTotALrhzQiI\nKqAtrVoV8PDD4Zs34VxuuSUFyCbSOzUak0qzZlsZMCAah3SKiHhDhYIP/vzn9kAHvDmrsBjowp/+\nlI6JkGkFH3poEM5kUzv9jnIOJViby9ixxX4HERGJKCoUfPDlL/ekX7+5wBDcXTDqMDCE/v1n8/Wv\n93Vxv43TokUKPXqsxjmLEqmXHwJAJvfe29nvICIiEUWFgk/efLN/7T03V1gsBCp5663BLu7THbfe\nmgL0JnJHP1SQlHSAa67RIlAiIqdSoeCTgQPbc801AWA07ly7Xw2M5hvfWE2fPpG3RsEPfzgIKCUy\nZ2msBvoyePA2EvQXISJyGr0t+mjatNGkpm6o/akxwyXLgFY0bx5g6tRxLiRzX4sWqWRlrQYi8dT+\nSqAHt93W0u8gIiIRR4WCj9LSUpgypRroAcxvxJ4KgAymT29JUlLk/i+96aZEoBeRN6fCAeA43/52\nb7+DiIhEnMj9VIkTt98+mFGj5uJcgljbgD0sBsZy003LuOKKyJ52+KGHhuCMftjhd5QzdKJbt3Wk\npenPQUTkTJ68Mxpjuhlj/mKM2WKMKTPGbDLG/MIYk+xFe9Fu5swxNG++BsgE9tXjmRuBgXTpMpfX\nX7/Em3AuysxMo3PnlThDQyPFdmAIX/lKZAwlFRGJNF59heoLGOD/A/oDE4HvAL/xqL2o1rRpMrNn\ndyAhoQo4DtRlLP9WIJOmTbdSUJBDQkJ0fNDdcEMN0If6FUReKgJqeOCBbL+DiIhEJE8KBWvtB9ba\n/2etnWmt3WatfRd4GrjBi/ZiwdChHXjnnQqMaYEzOuBckxOtAVqQknKYQKA9bdumhyekCx55ZAhO\nMbTV7yi10sjIWEtWVnO/g4iIRKRwXpTNoG5flePWNdf04t//PkJSUhBoBczDuaZ/0n6cTo99adFi\nK2vWZNK3b+QNhTyXTp0yaNeuAGjtdxTgCDCE8eNL/Q4iIhKxwlIoGGN6Ad8FpoSjvWh21VU9KSpq\nw8CBi4GRQCrO6fFtQFtgIFdf/Ql79w6hd+9MH5M23Fe+Uomz9kOJz0lWA035/vezfM4hIhK5jLV1\nn1LXGPME8Mg5NrFAP2vtZzMIGWMuAD4BZllr7z7P/nOAgnHjxtGy5elj2vPy8sjLy6tz1liwdOke\nnnpqI4WFkJQEubmJ/Oxng+jWLcPvaI2yefN+evdugXMJ5SIfk8wlJaUjFRVaLVJEokt+fj75+fmn\nPVZaWsqcOXMAcq21Abfaqm+hkInTNf9ctlhrq2u37wR8DCyw1t5Rh/3nAAUFBQXk5OTUOZdEn9at\n51BS0h7nzIIfaoASRo5cw8KF433KICLinkAgQG5uLrhcKCTVZ2Nr7WGclYfOq/ZMwixgKXBn/aNJ\nLJsw4Sivvz4Op2OjH50x1wKDueOO6Lx8IyISLl7No9AJ53LDduBhoJ0xpr0xRivuCAAPPzwAqAA2\n+ZSgBGP2c+edA3xqX0QkOtTrjEI9XIEzL3EPPh/nZ3D6MCR61KZEkSFDsmjefDZHj/pVO3YnO3s9\nSUmqXUVEzsWreRRettYmnnFLsNaqSJDPXH75pzhzc50Ic8tbga584xtNwtyuiEj00eT24puHHsoG\nqoDNYW55J1DKD34wJMztiohEHxUK4ptRo/qSlrYQ766AhdKeCy5YRcuWTcPcrohI9FGhIL669NLD\nOGs/VIapxX1ANtdfX/dhwSIi8UyFgvjqwQd7A0HCN/phE1DBQw8NClN7IiLRTYWC+OqSSwbQpMli\nnEEx4dCMzMyVdOvWKkztiYhENxUK4itjDJdeeojwXH44BAxmwoRwj7IQEYleKhTEdw891BvnjMJ6\nj1sqBIL87Gf9PW5HRCR2qFAQ340f34+mTZfh/eWHlrRuHSA7u43H7YiIxA4VCuI7YwyXXVYM9APK\nPGrlIDCI667zav8iIrFJhYJEhB//uB/O7N6rPWphA3CC//ovTbIkIlIfKhQkIlx8cRbNmq3Cu6VA\nMmnfPkBWVmuP9i8iEptUKEjEuPbaUmAocMDlPe8A+vGNbwRd3q+ISOxToSAR44knBuMMkdzg8p63\nAqX87Gc5Lu9XRCT2qVCQiJGV1YquXZcCXXFWJHdDEMimb98AbdqkubRPEZH4oUJBIsr3vpcGdAPW\nuLTHVUAHHn5YfRNERBpChYJElB/8IIekpO3AEZf2WElq6jpuv32wS/sTEYkvKhQkoiQlJXDVVUXA\ncOBwI/e2H8jhK185gDHhWktCRCS2qFCQiPP88xcCJ2j8lM4bgDKefnpo40OJiMQpFQoScbp3b83A\ngUuA/jgFQ0OUA4MYMmQZXbu2dC+ciEicUaEgEWny5B5ABrC0gXtYBrTgued6uhdKRCQOqVCQiHTZ\nZd3p0mUB0Iv6Lz9dAfSmW7eFjBvXzf1wIiJxRIWCRKw//akD0BZYXM9nLgXaMGVKB/dDiYjEGRUK\nQn5+vt8RzmrChF706TMfGACU1PFZJUA/+vSZz9VX9/IuXCNF6jGPZTrm4adjHhs8KxSMMf80xmw3\nxpQbY/YYY14xxnT0qj1puEj+Y54+vV/tvcI6PmMDkMC0aX09SuSOSD7msUrHPPx0zGODl2cUZgFf\nB/oANwA9gWketicxaODAdnzzm6uAUTgdFM9lGTCSW29dzuDB7b0PJyISBzwrFKy1z1prl1hrd1pr\nFwFPAiONMV6tIywx6q9/vYSOHecDfYGiEFttAXrTtu1iXnppfNiyiYjEurD0UTDGtAZuAeZba2vC\n0abEjoQEw7JlA2nSZCfOkMnVZ2yxDmhGSso+li7NJjFRXW9ERNyS5OXOjTFPAt8F0oCFwLXneUoT\ngPXrGzsjn9RHaWkpgUDA7xjn9dprJdx22wrKynoDr+AMm0wB+tGkySqmTm3L4cPHOdzYmZ/DIFqO\neSzRMQ8/HfPwOuWzs4mb+zXW1n05X2PME8Aj59jEAv2stRtrt28NtMZZDvDnwBFrbchiwRhzM/Ba\nnQOJiIjImW6x1r7u1s7qWyhkApnn2WyLtbb6LM+9ANgJjLLWnnVgfO3+rwK20fC5e0VEROJREyAL\n+MBa69q51XoVCo1qyJiuOAXAeGvtnLA0KiIiIo3iSaFgjBkOXATMw5kBpxfwK5xp9gZaa6tcb1RE\nRERc51X38DKcuRM+wpkp53+AFThnE1QkiIiIRImwXXoQERGR6KMB5yIiIhKSCgUREREJKeyFgjHm\nPmPM1trFohYZYy46z/bjjTEFxpgTxpiNxphvhStrrKjPMTfGfNUYM8MYc8AYU2qMWWCMuTKceWNB\nfX/PT3neaGNMlTFGs9TUUwPeW1KMMb8xxmyrfX/ZYoy5PUxxY0IDjvktxpgVxpjjtYsFvlg7347U\ngTFmrDHmbWPMbmNM0BjzlTo8p9GfoWEtFIwxNwGTcCZfGgqsBD4wxrQJsX0W8C4wExgCPAv8xRhz\nRTjyxoL6HnNgHDADmADkAB8D7xhjhoQhbkxowDE/+byWwMs4nYClHhp4zKcBlwJ34Cxel4ez/KjU\nQQPez0fj/H7/D9Af+E9gOPDnsASODek4AwPuxZng8Jxc+wy11obtBiwCnj3lZwPsAh4Osf1vgVVn\nPJYPvBfO3NF8q+8xD7GPNcDP/H4t0XJr6DGv/d3+Jc4bb8Dv1xFNtwa8t1wNFAMZfmeP1lsDjvkP\ngU1nPPZdYIffryUab0AQ+Mp5tnHlMzRsZxSMMclALk5lA4B1Un+Es4bw2Yzki9+uPjjH9nKKBh7z\nM/dhgOY4b6pyHg095saYO4DuOIWC1EMDj/l/4KxL/ogxZpcxZoMx5iljjKtz5MeqBh7zhUAXY8yE\n2n20B74O/MvbtHHNlc/QcF56aAMkAvvPeHw/0CHEczqE2L6FMSbV3XgxqSHH/EwP4ZzuetPFXLGs\n3sfcGNMbeBxnfvagt/FiUkN+z3sAY4EBwPXA/Tinwl/wKGOsqfcxt9YuAG4F/maMqQT24kzI910P\nc8Y7Vz5DNepBQqpdpOtR4OvW2kN+54lFxpgEnIXQfm6tLTr5sI+R4kUCzqnbm621y6y17wMPAN/S\nlxBvGGP641wj/wVO/6ercM6i/cnHWFIHni4zfYZDQA3Q/ozH2wP7QjxnX4jtj1hrK9yNF5MacswB\nMMZ8A6eT0X9aaz/2Jl5Mqu8xbw4MAy40xpz8NpuAc9WnErjSWvuJR1ljRUN+z/cCu621x055bD1O\nkdYZKDrrs+SkhhzzHwHzrbW/r/15jTHmXmCuMean1tozv/lK47nyGRq2MwrWmbq5ALjs5GO1178v\nAxaEeNrCU7evdWXt43IeDTzmGGPygBeBb9R+05I6asAxPwIMBC7E6ZU8BJiCM/X5EOCsK63K5xr4\nez4f6GSMSTvlsWycswy7PIoaMxp4zNOAM1cWDuL03tdZNG+48xka5l6aN+KsA/FNoC/OKafDQNva\nf38CePmU7bOAozg9N7NxhoRUApf73eM0Wm4NOOY31x7j7+BUnidvLfx+LdFyq+8xP8vzNerB42OO\n0+9mO/A3oB/OsOANwBS/X0u03BpwzL8FVNS+t3QHRgNLgAV+v5ZoudX+3g7B+WIRBH5Q+3OXEMfc\nlc9QP17ovTjLTZfjVDXDTvm3vwKzzth+HE7lWg5sAm7z+39WtN3qc8xx5k2oOcvtf/1+HdF0q+/v\n+RnPVaEQhmOOM3fCB8Cx2qLhd0Cq368jmm4NOOb3Aatrj/kunHkVOvr9OqLlBlxSWyCc9f3Zq89Q\nLQolIiIiIWnUg4iIiISkQkFERERCUqEgIiIiIalQEBERkZBUKIiIiEhIKhREREQkJBUKIiIiEpIK\nBREREQlJhYKIiIiEpEJBREREQlKhICIiIiH9/1HB8GtyL0iSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f32104b8ba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data=y\n",
    "ind=np.argsort(data)\n",
    "clust=np.zeros(data.shape[0],dtype=int)\n",
    "nclust=0;\n",
    "start=0\n",
    "stop=len(data)-1\n",
    "\n",
    "while start-stop<1:\n",
    "\n",
    "    nsteps = int(input('How many steps? '))\n",
    "    \n",
    "    for ns in np.arange(0,nsteps):\n",
    "    \n",
    "        if ind[start]==0:\n",
    "            Tracer()()\n",
    "    \n",
    "        # Retrieve cluster assignment of neighbors\n",
    "        start_neigh_clust=clust[neigh[ind[start]]]\n",
    "        stop_neigh_clust=clust[neigh[ind[stop]]]\n",
    "\n",
    "        # How many clusters is it touching?        \n",
    "        start_neigh_unique=np.unique(start_neigh_clust)\n",
    "        start_neigh_valid=start_neigh_unique[start_neigh_unique>0] # Remove border and unassigned vertices\n",
    "        stop_neigh_unique=np.unique(stop_neigh_clust)\n",
    "        stop_neigh_valid=stop_neigh_unique[stop_neigh_unique>0] # Remove border and unassigned vertices\n",
    "\n",
    "        # Is this a new cluster?\n",
    "        if len(start_neigh_valid)==0:\n",
    "            nclust+=1\n",
    "            clust[ind[start]]=nclust\n",
    "            start+=1\n",
    "        if len(stop_neigh_valid)==0:\n",
    "            nclust+=1\n",
    "            clust[ind[stop]]=nclust\n",
    "            stop-=1\n",
    "\n",
    "        # No, is it a border vertice (does it touch more than one cluster or only a border)?\n",
    "        if len(start_neigh_valid) > 1:\n",
    "            clust[ind[start]]=-1\n",
    "            start+=1\n",
    "        if len(stop_neigh_valid) > 1:\n",
    "            clust[ind[stop]]=-1\n",
    "            stop-=1\n",
    "\n",
    "        # No, then assign it to the unique cluster it touches\n",
    "        if len(start_neigh_valid)==1 and len(stop_neigh_valid)==1:\n",
    "            start_dist=np.abs(data[ind[start]]-np.mean(data[clust==start_neigh_valid[0]]))\n",
    "            stop_dist=np.abs(data[ind[stop]]-np.mean(data[clust==stop_neigh_valid[0]]))\n",
    "            if start_dist>stop_dist:\n",
    "                clust[ind[start]]=start_neigh_valid[0]\n",
    "                start+=1\n",
    "            elif stop_dist>start_dist:\n",
    "                clust[ind[stop]]=stop_neigh_valid[0]\n",
    "                stop-=1\n",
    "            else:\n",
    "                clust[ind[start]]=start_neigh_valid[0]\n",
    "                clust[ind[stop]]=stop_neigh_valid[0]\n",
    "                start+=1\n",
    "                stop-=1\n",
    "        else:\n",
    "            if len(start_neigh_valid)==1:\n",
    "                clust[ind[start]]=start_neigh_valid[0]\n",
    "                start+=1\n",
    "            if len(stop_neigh_valid)==1:\n",
    "                clust[ind[stop]]=stop_neigh_valid[0]\n",
    "                stop-=1\n",
    "        \n",
    "        if start-stop>1:\n",
    "            break\n",
    "                \n",
    "    clid=np.unique(clust)\n",
    "    clid=clid[np.logical_and(clid!=-1,clid!=0)]\n",
    "    for cl in clid:\n",
    "        xi=np.where(clust==cl)\n",
    "        xc=np.hstack((x[xi],x[xi][::-1]))\n",
    "        yc=np.hstack((y[xi],np.zeros(len(xi[0]))))\n",
    "        plt.fill(xc,yc,'b')\n",
    "    for xi in x[np.where(clust==-1)]:\n",
    "        plt.plot([xi,xi],[-1,1],'r')\n",
    "    plt.plot(x,y)\n",
    "    plt.show()        \n",
    "        \n",
    "if np.sum(clust==0)>0:\n",
    "    raise ValueError(\"Unassigned vertice detected.\")\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}