{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python implementation of the Kneedle algorithm\n",
    "Finding a “Kneedle” in a Haystack: Detecting Knee Points in System Behavior  \n",
    "Ville Satopa, Jeannie Albrecht, David Irwin, and Barath Raghavan  \n",
    "https://www1.icsi.berkeley.edu/~barath/papers/kneedle-simplex11.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finding the knee from figure 2 from the paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def figure2():\n",
    "    x = np.linspace(0.0, 1, 10)\n",
    "    with np.errstate(divide='ignore'):\n",
    "        return x,np.true_divide(-1, x + 0.1) + 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 0:  Raw input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x,y = figure2()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "invert = False\n",
    "if invert:\n",
    "    x = [max(x) - x_ for x_ in x]\n",
    "    y = [max(y) - y_ for y_ in y]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not np.array_equal(np.array(x), np.sort(x)):\n",
    "    raise ValueError('x needs to be sorted')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 1: Fit a spline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.interpolate import interp1d,UnivariateSpline,InterpolatedUnivariateSpline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = len(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Ds = the finite set of x- and y-values that define a smooth curve, \n",
    "# one that has been fit to a smoothing spline.\n",
    "uspline = interp1d(x,y)\n",
    "\n",
    "Ds_x = np.linspace(np.min(x), np.max(x), N)\n",
    "Ds_y = uspline(Ds_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHPlJREFUeJzt3Xl0VNeB5/HvVWnfd4R2oQ2w2QUIL4ljO/Fu7HTmJOnu\nrJ3DeMaxM6fTJ+skOZPuTKfTPd1Z3EmG6XZnPJ0Zp0+8YBOIg+O2nZhdrGaRECCQhIR2tC9VdecP\nKQ7BAglUqlf16vc5h3Mo1dN9v4tKP56eXr1rrLWIiIh7RDkdQEREAkvFLiLiMip2ERGXUbGLiLiM\nil1ExGVU7CIiLqNiFxFxGRW7iIjLqNhFRFwm2omdZmdn29LSUid2LSISturq6rqstTkzbedIsZeW\nlrJ//34ndi0iEraMMedms51OxYiIuIyKXUTEZVTsIiIuo2IXEXEZFbuIiMuo2EVEXEbFLiLiMgG7\njt0Y4wH2A63W2gcDNa6ISLiwfj9Dg5cY7Oti6FI3owM9jA/2MDHUg3+4FzvSR/Gdn6Fg0U3zmiOQ\nb1D6HHACSA3gmCIiQeX3+Ri41MNgXxfD/V2MTZWzd6h3spxH+4ga6yd6/BIxE/3EewdI9A2QZIdI\ntkMkGx/JVxvbGo6euyU8it0YUwg8AHwL+PNAjCkiMhdjo8P0dbUx0N3GSF8HY4Pd+IZ6scN9MNpH\n1Nglosf7ifVOlbN/kGQ7SLIdJs1Y0q4y7oT1MGCSGDLJjHiSGY1OZSChEG9sKv74dExCOlEJ6UQn\nZRCbnEl8ShaJadkkpWWRkprBCo9n3uceqCP27wJfAFKutoExZhOwCaC4uDhAuxWRSDE8eIm+rnYG\ne9oZ6bvIeH8HvsFO7FAXnpEeYsd6SJjoJdnXR5q/n2QzwgJgwTRjjdoYBkwyw1HJjHhSGIrNpi9m\nEb64tMlyjk8jKjGDmKQM4lKySEjNJCE1i5T0bBKTUsmMiiIz2P8A12HOxW6MeRDosNbWGWPuuNp2\n1trNwGaAmpoaO9f9ikj4sn4/gwN99He3MdDTzmhfBxP9HXgHuzBDnXhGJ4s60dtHsrePNNtPohkj\ncZqxxm00fSaVAU86w9HptCUW0RyfiU3MIioph5jUXBLSc0lMzSYpLZvkjGziE5KID/qsgycQR+y3\nAg8bY+4H4oFUY8y/Wmv/NABji0gYGRroo7utif6O84x0t+Dta8UMdRI92k3seC+JE32k+PpIt/2k\nmIlpf8QftTH0mXQGPGmMxGTQl1TGufgsbGIWnuQcYlNziU/PJTkjj9TshaSkZpAbFUVusCcbwuZc\n7NbaLwNfBpg6Yv8LlbqIu/h9Pno6L9B38RyDnc2M9bTg779A9GAb8SMXSZnoJNPfTSrDJF3xucM2\njr6oNAY96QzFZtMbV8XZ+ExIyiYqOYe4tBwS0vNIzsgjLTuPxKRU8qKiyHNkpu7gyG17RSR0jI4M\n0d12nksd5xjubsbb2wL9bcQOt5M41knaRCdZtods4yP7ss/zWUO3yaAvOpvehBIuJqzDn7KQ6PQC\nErIKSc0tJmthKUkp6dOeQpH5E9Bit9a+DrweyDFF5MZYv5/+3k562idPjYz1tOC7dAHPYBtxIxdJ\nHu8k09dFBgMUAAWXfe6wjaM7Kov+mBxaU1fSlLQQk7qQ2MwiknOKSF9QQmZuAbkxsToFEoJ0xC4S\nxnxeLxebG+luPsFwWwO25wzx/U1kjLWwwHeRNDPxrsv2ukmj15PNQOwCuhJX4E/Ow5NeQEJmISm5\nxWTklZKalklRlN6YHq5U7CIhzjsxzsXm03SfP87IxVPY7tMkDJwjY6yFPF87+cZH/tS2wzaOds9C\nuhPKaEu+HVLzickoICm7mLQFxWTllZAVF0+WozOS+aZiFwkB3olx2s810NN88vflPXiOzNEW8vwX\nKTC+d06VDNs42qLz6UoopzXtLjxZ5STnV5NTsoTsvGIW6Ug74qnYRYJkYnyM9vNT5d3egOk5Q/zA\nObLGWljg76DQ+Cic2nbIxtMenU9nUiUtqe/Hk1NB8sIqcouXkJVXRLnKW65BxS4SQN6JcdqaTtDT\nXP9OeScMniNzrIU8fwdFxk/R1LaDNoH26HwuJlXTnHYPnuwKUvKryClZSlZugcpbbpiKXeQG+bxe\nzjccpKthD/7Wg6T1Had4/DRFZuyd8h6wCbRHF3AxeQnNafcTnV1OytRpk8ycfCpU3jIPVOwis+Dz\nemk+dYjO+t1/UOJlZowyJs97n4ut4MiCR/DkLyelYDG5JUvIyF5IpcpbgkzFLnKF35f4HvytB94p\n8VIzRimXl/hGogtXk1O1nsKK5SyJ1reThAa9EiWiXV7ivtaDpPcdu2qJewpWkVtdqxKXkKdXp0SM\nK0s8re8YJe8q8fLfl3jVegorV6jEJezoFSuuNOsSz30YT+Fqlbi4il7F4gpd7ec5u2crvpYD1y7x\nglXkVK+nqHKlSlxcS69sCUvW7+fs8X1c3Pc8Wa2vUeVtIJvJEj+vEpcIp1e7hI3xsVHq92xn+MjL\nlHS9ySI6WQQ0RFexq+QxclY/TNlN61msEpcIp+8ACWl9Xe2ceut5PKe2UzWwj2VmhBEbS33SGs6X\nP86iDR+kKr/E6ZgiIUXFLiHnfMMhLux9gdRzr1I9foy1xtJFOsez7iZ26QNUb3iQlUlXXTddJOKp\n2MVx3olxGvb/mv5DWyjoeINie4Fi4LSnjL1FnyJr9SNUrLiNbI/H6agiYUHFLo4YuNRDw1svYE9u\np6J/F0sZZNx6OJmwkgtln6C49lHKS6opdzqoSBhSsUvQXGiq5/yu50hq+hXVo0dYY3z0ksKptFvx\nLL6Pyls2sjwt0+mYImFPxS7zxu/zcergG/Qc3EJe2+uU+ZvIB85FFXJg4UdIW7WRqjV3sVZXsYgE\nlL6jJKBGhgao3/ky48e3sqj3Larpw2ujqI+7id0lf07B+g9SUrEMXcciMn9U7DJnXRfOcfqtnxN3\n5lcsHq5jpZlgwCbQkFpLU+W9VN76KDdlLXA6pkjEULHLDbtw9iQtW/4bq3t/yXrj54JZwKHcR0he\n/hBV6+5hTVy80xFFIpKKXa5be3Mj5178Jqu7tpJFFPsXfIgFd2yidPEa8rWohIjjVOwya10XznH6\nhW+yquNFMrEcyNlI2aNfp7agzOloInIZFbvMqKejlYbn/pKV7T9nDT4OZN5P0SNfZ31JtdPRRGQa\nKna5qkvdFzn+3LdY0fosaxnnQPoHyN/4DdYtusnpaCJyDSp2eZf+vm6OPffX3Hz+X1nPKAdT30f2\ng99gbfVKp6OJyCyo2OUdQwN9HH3uOyxp+gkbGOJA8u1kPPAN1ixd63Q0EbkOKnZhZGiAwy/8HdWN\nT1NLP4cSakm+9+usXnGr09FE5Aao2CPY6MgQh178LhX1m6mljyPxNXTc/VVW1tzpdDQRmQMVewQa\nHxvl4JYfUHb8h9TSw7HYFXTetZnl6+9xOpqIBICKPYJMjI9x8OUfUXT0KdbTycmYpXS893vcfNvD\nTkcTkQBSsUcAn9fLgV9sZuGh77HOttMQXUXnbX/Dsvc8itE7RUVcZ87FbowpAp4BFgAW2Gyt/d5c\nx5W58/t8HHzlJ2Tv/3vW+ls47VnEoQ0/ZsWdH1ahi7hYII7YvcDnrbUHjDEpQJ0xZoe19ngAxpYb\nYP1+Du74KRl7/o41/iaaooo5UPs9Vr7/Y0RpeTkR15tzsVtr24C2qb8PGGNOAAWAij3IrN/Pkdf/\njaS3vsNq32maTT77a/6WVfd+mlItZiESMQL63W6MKQVWAXsCOa5cm/X7efu3W4h989us8J7kglnA\nvpXfYtUDmyiKiXU6nogEWcCK3RiTDDwH/Bdrbf80z28CNgEUFxcHarcR79jObZjX/zvLxo/STjZ7\nb/4Gqx5+nPzYOKejiYhDAlLsxpgYJkv9p9ba56fbxlq7GdgMUFNTYwOx30g2PHiJhh9+mJXDu+gk\ngz1LvszKjU+SF5/odDQRcVggrooxwD8DJ6y1fz/3SDKTsdFhTv/gEZaNHmR3+ZOs/NAXWZ+Y7HQs\nEQkRgThivxX4GHDUGHNo6mNfsdZuC8DYcoWJ8TGOf/+PWDV2gL0r/4raR59wOpKIhJhAXBXzW8AE\nIIvMwOf1cvipP6ZmeCd7lnyZ9Sp1EZmG3qUSJqzfT90/foKa/lfZtehJ1n/4S05HEpEQpWIPA9bv\nZ8+P/yPrereyq/DTbPj4XzodSURCmIo9DOx++vPUdvwbu3M/TO2n/4fTcUQkxKnYQ9yuZ77Ghpan\n2ZvxIOsf+7Hu8SIiM1JLhLA9P/s2G858n/2pd7Pm8f+tUheRWVFThKi9L/yA9Sf+moOJt7Dis/8X\nj+71IiKzpGIPQXXb/oU1h77Gkfg1LH3yOWJ0ewARuQ4q9hBz+LVnWb7n8zTELqXyiS3E6RYBInKd\nVOwh5O3fvsTiNz5LU8wiCh5/mYSkFKcjiUgYUrGHiJN7d7Box2e44Mkn57GtpKZnOR1JRMKUij0E\nNB7+LfnbPk53VBYpm7aSnp3ndCQRCWMqdoedO1FH1gsfYZgkYj71Etl5ule9iMyNit1BrWeOkfiz\nP8KHB++fvkhecaXTkUTEBVTsDmlvbiTqmY1E42Xww89RWHGz05FExCVU7A7oaj/PxNMPkWwH6X70\nWUqX1DgdSURcRMUeZJe6LzKw+SGy/N20PvB/qFhxm9ORRMRlVOxBNHCph4s/epB8Xytn7v5fLF73\nfqcjiYgLqdiDZGRogOanHqJs4jQnbv8BN9++0elIIuJSKvYgGBsd5tQPNlI9fowj677Dyrs/6nQk\nEXExFfs8+93i08tH66hb+U3WPPAZpyOJiMup2OfR7xafXjW8k93VX2Tdo086HUlEIoCKfZ5Yv5+6\nH35ycvHpss9S+9GvOB1JRCKEin0eTC4+/Rjrel5mV8En2fCJbzkdSUQiiIp9Hux5+i+o7fgZu3P+\nA7V/9g9OxxGRCKNiD7Bdz3yN2pZ/Zm/GA6x77H9qnVIRCTq1TgD9bvHpupQ7WfP4M0R5PE5HEpEI\npGIPkH0vPvXO4tPLn3hWi0+LiGNU7AFQt+1fWH3wv3I0bhVLnvi5Fp8WEUep2Ofo94tPL6H8iS3E\nJyQ5HUlEIpyKfQ7eWXw6uoyCx7eSmJzmdCQRERX7jTq571UW7fgMbZ6F5PynX2jxaREJGSr2G9B4\n+C3yf/ExeqIySf6MFp8WkdCiYr8B9qUnGSEBzye3kJ1f4nQcEZE/oGK/Th2tZ6n0NXJm0Z+wsKTa\n6TgiIu+iYr9OZ3c+D0BezcMOJxERmV5Ait0Yc68xpt4Y02iM+VIgxgxVcWd30EYOpYvXOB1FRGRa\ncy52Y4wH+EfgPmAp8FFjzNK5jhuKRocHqRo6wPns23UPGBEJWYFop3VAo7X2jLV2HHgWcOWCnvV7\ntpNoxoi/6X6no4iIXFUgir0AaL7sccvUx/6AMWaTMWa/MWZ/Z2dnAHYbfKPHtjFs46iuVbGLSOgK\n2vkEa+1ma22NtbYmJycnWLsNGOv3U9z1G+qT1ui2ASIS0gJR7K1A0WWPC6c+5ipNJ+tYSCfjZXc7\nHUVE5JoCUez7gEpjTJkxJhb4CPBSAMYNKe37J6dUdssHHU4iInJtc75puLXWa4z5LPAK4AGettYe\nm3OyEJPe/GsaPeVUFJQ5HUVE5JoCshqEtXYbsC0QY4Wivq52qsaPs7foU1Q4HUZEZAa6GHsWGndt\nwWMsmSsfcjqKiMiMVOyz0fAKPaRSueq9TicREZmRin0G3olxqgZ2czrtFi1OLSJhQcU+g4a610hl\nCM/ie52OIiIyKyr2GVw6/DIT1kPlLa68S4KIuJCKfQYLL75BffwyUtIynY4iIjIrKvZruHD2JKX+\nZgaL73I6iojIrKnYr6F5zwsAFKx/1OEkIiKzp2K/hoSmV2k2+RRVLHM6iojIrKnYr2JooI/FI4do\nzX2P01FERK6Liv0qGnZtJdZ4Sb5Z914XkfCiYr+KiZO/ZNAmULXuHqejiIhcFxX7NKzfT1nPb2lI\nWUtsXLzTcURErouKfRqnj+4kh1685R9wOoqIyHVTsU+j88DL+K1h0YZHnI4iInLdVOzTyGp9jVMx\nVWTnFc28sYhIiFGxX6GrvZkqbwM9Be9zOoqIyA1RsV/hzK4XAchZ/bDDSUREboyK/QrRjb+ig0zK\nl21wOoqIyA1RsV9mfGyUqsF9NGXeionSP42IhCe112Ua9r5CshkhdonebSoi4UvFfpnBo79gzMZQ\nteEBp6OIiNwwFftlCjvfpD5hJYnJaU5HERG5YSr2Kc2nDlNo2xgpu9vpKCIic6Jin9K6d/IyxyIt\nqiEiYU7FPiX5/Gs0RRWTX1rtdBQRkTlRsQP9fd1Ujx6lbcEdTkcREZkzFTtwaudLxBgf6SsfdDqK\niMicqdgBf/0vuUQSlat1fxgRCX8RX+w+r5fySzs5lbqB6JhYp+OIiMxZxBf7qUNvkEk/VGpRDRFx\nh4gv9t5DW/FZQ+UtWlRDRNwh4os9t+116mNvIi1rgdNRREQCIqKL/WLLacp9Z7hUdKfTUUREAmZO\nxW6M+VtjzEljzBFjzAvGmPRABQuGpqlFNRau3ehwEhGRwJnrEfsO4GZr7XKgAfjy3CMFT9zZHVww\nuZRUr3Y6iohIwMyp2K21v7LWeqce7gYK5x4pOEaHB6keqqM56zYtqiEirhLIRvs0sD2A482r+j3b\nSTDjJNyke6+LiLtEz7SBMeZVIG+ap75qrd0ytc1XAS/w02uMswnYBFBcXHxDYQNp9Ng2hm0cVbX3\nOR1FRCSgZix2a+01b1BujPkk8CBwl7XWXmOczcBmgJqamqtuFwzW76e46zfUJ61hVUKSk1FERAJu\nrlfF3At8AXjYWjscmEjzr+lkHQvpZGLR+52OIiIScHM9x/4UkALsMMYcMsb8OACZ5l37vsnLHEs3\naFENEXGfGU/FXIu1tiJQQYIpveXfafSUU1FQ5nQUEZGAi7jr/Pq62qkaP07XwjucjiIiMi8irtgb\nd76Ix1gyVz/sdBQRkXkRccXOqV/RTRoVK253OomIyLyIqGL3ToxTNbCbM+m3EOXxOB1HRGReRFSx\nN+z/NakM4Vl8r9NRRETmTUQV+6UjWxm3Hio36Py6iLhXRBX7wotv0BC/jJS0TKejiIjMm4gp9gtn\nT1Lqb2aw5Jp3SBARCXsRU+zNe14AoGCd1jYVEXeLmGJPaHqVZpNPUcUyp6OIiMyriCj2oYE+Fo8c\nojX3PU5HERGZdxFR7A27thJrvCQv06IaIuJ+EVHsEye2M2ATqFr7AaejiIjMO9cXu/X7Ket9i1Mp\na4mNi3c6jojIvHN9sZ8+upMcevFV3ON0FBGRoHB9sXfWvYTfGhZt0GWOIhIZXF/sWRf+nVMxVWQt\nKHQ6iohIULi62Lvam6nyNtBT8D6no4iIBI2ri/3Mrsl3m+au2ehwEhGR4HF1sUc37qCDTBbdXOt0\nFBGRoHFtsY+PjVI9uJezmbdholw7TRGRd3Ft49XveYUkM0rc0vudjiIiElSuLfaht3/BmI2hqlbF\nLiKRxbXFXtj5JicTVpKYnOZ0FBGRoHJlsTefOkyhbWO07P1ORxERCTpXFnvr3hcBKK7Vu01FJPK4\nsthTzv2as1ElLCypdjqKiEjQua7Y+/u6qRp7m/a89zodRUTEEa4r9lM7XyLG+Ehf8aDTUUREHOG6\nYvfXb6ePZCpX6/4wIhKZXFXsPq+X8ku7aEytJTom1uk4IiKOcFWxnzr0Bpn0Q5UW1RCRyOWqYu89\n+DJeG0XlBt3NUUQil6uKPbf9DRpil5KWtcDpKCIijnFNsV9sOU257wz9RXc6HUVExFEBKXZjzOeN\nMdYYkx2I8W5E09SiGgvX6t2mIhLZ5lzsxpgi4APA+bnHuXFxZ1/lgsmluHqVkzFERBwXiCP2fwC+\nANgAjHVDRocHqR6qoznrdi2qISIRb04taIzZCLRaaw/PYttNxpj9xpj9nZ2dc9ntu9Tv2U6CGSfh\n5gcCOq6ISDiKnmkDY8yrQN40T30V+AqTp2FmZK3dDGwGqKmpCejR/eixbQzbOKrW3xvIYUVEwtKM\nxW6tvXu6jxtjlgFlwGFjDEAhcMAYs85a2x7QlNfK5/dT0vUb6pPWsCohKVi7FREJWTd8KsZae9Ra\nm2utLbXWlgItwOpgljpA08k68uhkonxWPziIiLhe2P+msX3f5KIaZRsedTiJiEhomPFUzGxNHbUH\nXXrLazR6yqnId2T3IiIhJ6yP2Pu62qkaP0Fnvm7RKyLyO2Fd7I07X8RjLFmrHnI6iohIyAjrYjen\nXqGbNCpW3O50FBGRkBG2xe6dGKdyYA+n028lyuNxOo6ISMgI22Jv2P9rUhkierEW1RARuVzYFnv/\n4ZcZtx4qNzzsdBQRkZAStsWe1/Em9fHLSUnLdDqKiEhICctiv3D2JKX+ZoZK7nI6iohIyAnLYm/e\nM7moRuE6vdtURORKYVnsCU2v0mzyKay42ekoIiIhJ+yKfWigj8Ujh2jNfa/TUUREQlLYFXvDrq3E\nGi/Jy7WohojIdMKu2CdObGfAJlC9VrfpFRGZTlgVu/X7Ket9i1Mp64iJjXM6johISAqrYj99dCc5\n9OKr0LtNRUSuJqyKvbPuJfzWUH7LI05HEREJWWFV7NHpBezPuI/M3AKno4iIhKyAraAUDGs/+Dng\nc07HEBEJaWF1xC4iIjNTsYuIuIyKXUTEZVTsIiIuo2IXEXEZFbuIiMuo2EVEXEbFLiLiMsZaG/yd\nGtMJnLvBT88GugIYJxxozpFBc44Mc5lzibU2Z6aNHCn2uTDG7LfW1jidI5g058igOUeGYMxZp2JE\nRFxGxS4i4jLhWOybnQ7gAM05MmjOkWHe5xx259hFROTawvGIXUREriFki90Yc68xpt4Y02iM+dI0\nzxtjzPennj9ijFntRM5AmsWc/2RqrkeNMTuNMSucyBlIM835su3WGmO8xpgPBTNfoM1mvsaYO4wx\nh4wxx4wxbwQ7Y6DN4nWdZox52RhzeGrOn3IiZyAZY542xnQYY96+yvPz21/W2pD7A3iA08AiIBY4\nDCy9Ypv7ge2AAWqBPU7nDsKcbwEypv5+XyTM+bLtXgO2AR9yOvc8f43TgeNA8dTjXKdzB2HOXwH+\nZurvOUAPEOt09jnO+z3AauDtqzw/r/0Vqkfs64BGa+0Za+048Cyw8YptNgLP2Em7gXRjzMJgBw2g\nGedsrd1pre2dergbKAxyxkCbzdcZ4AngOaAjmOHmwWzm+8fA89ba8wDW2kiYswVSjDEGSGay2L3B\njRlY1to3mZzH1cxrf4VqsRcAzZc9bpn62PVuE06udz5/xuT/+OFsxjkbYwqAR4EfBTHXfJnN17gK\nyDDGvG6MqTPGfDxo6ebHbOb8FLAEuAAcBT5nrfUHJ55j5rW/wmrNU5lkjHkfk8V+m9NZguC7wBet\ntf7JAzrXiwbWAHcBCcAuY8xua22Ds7Hm1T3AIeBOoBzYYYz5jbW239lY4StUi70VKLrsceHUx653\nm3Ayq/kYY5YD/wTcZ63tDlK2+TKbOdcAz06VejZwvzHGa619MTgRA2o2820Buq21Q8CQMeZNYAUQ\nrsU+mzl/Cvi2nTz53GiMOQssBvYGJ6Ij5rW/QvVUzD6g0hhTZoyJBT4CvHTFNi8BH5/67XItcMla\n2xbsoAE045yNMcXA88DHXHIEN+OcrbVl1tpSa20p8HPgP4dpqcPsXtdbgNuMMdHGmERgPXAiyDkD\naTZzPs/kTygYYxYA1cCZoKYMvnntr5A8YrfWeo0xnwVeYfK36k9ba48ZYx6bev7HTF4hcT/QCAwz\n+b9+2JrlnL8OZAE/nDqC9dowvoHSLOfsGrOZr7X2hDHml8ARwA/8k7V22kvmwsEsv8Z/CfzEGHOU\nyatEvmitDes7Phpj/h9wB5BtjGkBvgHEQHD6S+88FRFxmVA9FSMiIjdIxS4i4jIqdhERl1Gxi4i4\njIpdRMRlVOwiIi6jYhcRcRkVu4iIy/x/d2MzG2lbu54AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5bf7f9a6d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,y);\n",
    "plt.plot(Ds_x,Ds_y);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 2: Normalize the spline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def normalize(a):\n",
    "    \"\"\"return the normalized input array\"\"\"\n",
    "    return (a - min(a)) / (max(a) - min(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# x and y normalized to unit square\n",
    "xsn = normalize(Ds_x)\n",
    "ysn = normalize(Ds_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 3: Calculate the difference curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# the difference curve\n",
    "xd = xsn\n",
    "yd = ysn - xsn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VFX+//HXJ40USCAJRWrovYgBVEC6gmXtDddeFvtW\ndXe/us3d1d+u69pW7F1QkbWtDQgIilSlJ4TQE0oKkEAg/fz+ODcyxJQBJnMzM5/n4zEPZubezHzu\nJLxzcu6554gxBqWUUsElzO0ClFJK+Z6Gu1JKBSENd6WUCkIa7kopFYQ03JVSKghpuCulVBDScA9h\nIvJHEXnTud9ZRA6JSLiP32ObiEz08Wu+KiIPO/dHi8hGX76+L9X3GYtIWxFZKCIHReQxsV4Rkf0i\nsszdylWg03BvRE6w5YpInMdzt4jIAhfLqpUxZocxprkxptLtWo6HMWaRMaZ3Y7y2E7aPiUiBc5t1\nMq9Xy2d8G5APxBtjfgWMAiYBHY0xw0+uehXqNNwbXzhw78m+iBM0+v3yr7OBnwKDgfbAcz5+/S7A\nBnP0SsIuwDZjTPHxvpCIRPi0MpcF2/G4QcOi8f0D+LWItKxto4icKSLLRaTQ+fdMj20LROSvIvIN\ncBjo5jz3sIgsdv7E/1hEkkTkLREpcl4jxeM1nhCRnc62lSIyuo46UkTEiEiEiJzhvHb1rUREtjn7\nhYnIAyKy2WnNvisiiR6vc62IbHe2/b6+D0ZEzhWRDU63RI6I/Np5fqyIZIvI70Qk3/kL6Jo6XmOs\niGR7PN4mIr8WkTXOZ/qOiER7bD9fRFaJyAHnMxxUT4nlwBFgjzGm1Bgzp77jcV6/q4h85RzTHCDZ\nY5vnZ/wqcD1wn/MZ/wx4Eaj+7P/UUL3Osd4vImuAYud124vI+yKSJyJbReQej/3/6Hy/XnfqWy8i\nqR7bO4nIbOdrC0TkaY9tN4lIutguoy9EpEs9n8Eop9YDzs/eDc7zC0TkFo/9bhCRrz0eGxG5U0Q2\nAZtE5FkR+WeN1/5QRH7p3K/zWBVgjNFbI92AbcBEYDbwsPPcLcAC534isB+4FogArnYeJznbFwA7\ngP7O9kjnuSygO5AAbAAynfeJAF4HXvGo4adAkrPtV8AeINrZ9kfgTed+CmCAiBrHEAl8BfzdeXwv\nsAToCDTDtmZnONv6AYeAs5xt/wIqgIl1fD67gdHO/VbAUOf+WOfr/uW8zhigGOjtbH/V4/McC2TX\n+MyXYVvaiUA6MM3ZdiqQC4zA/kV1vbN/szrqaw8UOe8X5uX3/FuPus8CDtb1GXseh/P4BuBrj8f1\n1uvcXwV0AmKwjbWVwENAFNAN2AKc4/H9LgHOdV7v78ASZ1s4sBp4HIgDooFRzrYLsT9zfbE/R/8H\nLK7j+Ls4x3w19mcnCRji8fN8Sz3Ha4A5zvctxvn8dgLi8TNyxPm+1HusejPacveTh4C7RaR1jefP\nAzYZY94wxlQYY2YAGcAFHvu8aoxZ72wvd557xRiz2RhTCHwGbDbGzDXGVADvYUMBAGPMm8aYAufr\nH8OGzvH0UT+J/c9a3QqfBvzeGJNtjCnFBsZlYv+Mvgz4xBiz0Nn2IFBVz2uXA/1EJN4Ys98Y812N\n7Q8a22L+CvgfcIW3NRtjdhlj9gEfA0Oc528DnjPGLDXGVBpjXgNKgdNrvoCIRAJfAHdgQ+VFcbrF\nRORrEbmglq/pDAzzqHuh8/4nypt6nzTG7DTGHHHeu7Ux5s/GmDJjzBbgBeAqj/2/NsZ8amy//xvY\nLieA4djQ/I0xptgYU2KMqW5VT8P+ck93fsb+Bgypo/U+FZhrjJlhjCl3fvZWHccx/90Ys885nkXY\nwK/+a/My4FtjzC4vjzWkabj7gTFmHfAJ8ECNTe2B7TWe2w508Hi8s5aX3Otx/0gtj5tXP3C6KNKd\nLooD2NZ+Ml5wugrGAlONMdUh3QX4r/Mn9wFsy7gSaOsczw/1Gtt3XFDPW1yKbUVud7oyzvDYtt8c\n2/e83Xl9b+zxuH+Yo59HF+BX1bU79Xeq43XHA1HGmDeBK4Gu2ICPB/oAX9fyNe3rqPtEeVPvzhr7\nt6+x/++w35tqNT+baOcXcydguxPetdXxhMdr7gOEY39Oq3UCNh/fYR7D8+fHADOxfwWA/cXxlkdN\nDR1rSNOTFv7zB+A74DGP53Zhf0g9dQY+93h8wtN2iu1fvw+YAKw3xlSJyH7sf0xvvvYv2D/Nizw2\n7QRuMsZ8U8vX7Mb+6V79OBb7Z3mtjDHLgQudVvJdwLvYcABoJSJxHkHZGVjXUN0N2An81RjzVy/2\nre4GwxhTIiI/AeYDy4GZxpj9tXzN7jrqPtHvoTf1er72TmCrMabnCb5XZxGJqCXgq+t4q5avq+11\n6hrpUwzEejxuV8s+NT+rGcCXIvIItnvqYo/3OdFjDQnacvcTY0wW8A7gedLnU6CXiEx1ToZdie23\n/sRHb9sC23edB0SIyENAfENfJCKdsEF7nTEms8bm6cBfq/8kF5HWInKhs20WcL5zQi0K+DN1/IyJ\nSJSIXCMiCU53UxE/7sL5k7PfaOB8bJfTyXgBmCYiI8SKE5HzRKRFLft+jW3V/llEqvuz5wO9sC3e\nHzHGbAdWeNQ9imO72BqzXrDnGg6KPckaIyLhIjJARIZ58V7LsL+cHnHeJ1pERjrbpgO/FZH+ACKS\nICKX1/E6bwETReQK52c6SUSqu8VWAZeISKyI9ABubqgoY8z32OGiLwJfGGMO+OBYQ4KGu3/9GXuy\nCgBjTAE2tH6F7b64DzjfGJPvo/f7AvtXQCa2e6CE2rt5apqA/fN2lhwdMbPe2fYE8BG2NXUQe3J1\nhHM864E7gbexQbEfyK754h6uBbaJSBG2X9dzRMwe5+t3YQNjmjEmw4va62SMWQHcCjztvHYW9qRe\nbfsWYodCnu7UsBn7V8hw4EYRubWOt5mK/Tz2Yf9ae90f9Tr7V2J/noYAWzkaiglevFcl9hdRD+xJ\n/GxsdxTGmP8CjwIzne/VOmBKHa+zA9vV9ivsZ7CKo/36jwNl2G7E1zjaxdKQt7EDBt72xbGGiuqz\n0Eo1GSIyFjvCpKPbtSgVqLTlrpRSQUjDXSmlgpB2yyilVBDSlrtSSgUh18a5Jycnm5SUFLfeXiml\nAtLKlSvzjTE1r3b/EdfCPSUlhRUrVrj19kopFZBExKurnrVbRimlgpCGu1JKBSENd6WUCkIa7kop\nFYQ03JVSKgg1GO4i8rLYRZ5rnW7Vma3uSRHJEru02VDfl6mUUup4eNNyfxWYXM/2KUBP53Yb8OzJ\nl6WUUupkNDjO3RizUDwWXK7FhcDrzqopS0SkpYicYozZ7aMalVIqoJSUV1J4pJyiI+UU1rgVHalg\naJeWjO7Z4HVIJ8UXFzF14Ng5wrOd534U7iJyG7Z1T+fOnX3w1kop5XvGGIrLag/oorpCu6Tih/tl\nFfUtHQy3j+0eEOHuNWPM88DzAKmpqTpjmVKqUVVUVrH/cDn7issoKC7lwOHyulvUJRU/PF90pJyK\nqrojSgTioyOJj4kgISaShJhI2iVEkxATSXxMJPHRkT88X32Lr/43OoKI8MYfy+KLcM/h6LqXAB2d\n55RSyqfKKqrYf7iMgkNlPwT20ftl7Csu9bhfRuGRcuqa+DYiTI4GrvNv58RYEmIi6g/nmEhaNIsg\nLKzBpYhd5Ytw/wi4S0RmYpcXK9T+dqWUN0rKK20YH7JBva/YI6gP/TiwD5bUXLvbChNoFRtFYlwU\nSc2j6NsunsS4o4+r77eKjfohrGOjwhFp2gF9MhoMdxGZAYwFkkUkG7suZPWq8NOxizyfi13f8TBw\nY2MVq5Rq+soqqthbVMLeohL2FJWwp7CEvIOlP7SmfwjsQ2UUl1XW+hoRYUKruCiSnFAe2LHlD/cT\nPZ63wd2MljGRTb4l7W/ejJa5uoHtBrsoslIqiBljKCqpsKFdeDS49xSVsLf636IS8g+V/ehro8LD\njmlJpyTFkhTX7JhW9Q+BHdeM+JiIoG5V+4NrU/4qpZqOisoq8g+VOYF9xAnt0h8F+ZHyH7e0W8VG\n0i4hhnbxzRjUMYG28dG0i4+mbUI0pyTY+wkxkRrWfqbhrlSQKy6tOKZ1XX1/d+HRrpO8g6XUHBwS\nGS60aRFNu4Ro+p0Sz7jebTglwYZ2OyfA28Q3Izoy3J0DU/XScFcqwBljyDtUyvaCw2zNL2Z7QTHb\n8g+zraCYHfsO13oSskV0hA3ohGh6tW1Bu4ToH1rc7RLsLTE2SvuxA5iGu1IBwBhD3sFSJ7xtcG9z\nQnx7QfExJyYjwoROibF0SYrltC6taOd0j3iGd2yU/tcPdvodVqqJMMawt6iUbQW29b3VCe7qQPfs\n744IEzo7AT6iWyIpSXF0SYqla3Ic7VvGEOmHi2RU06bhrpQfVVUZ9h4s+aHbZFtBMdud+zUDPDLc\ntsC7JsVxZvdkUpJjSUmKIyUpjvYto/1ylaMKXBruSvmYMYY9RSVszT/abVLd+t6+r5iS8qPzjkSF\nh9E5KZaUpFhG9kgmJTmOlKRYJ8BjCNc+b3WCNNyVOgnGGHIOHGFdTiFrcwpZm1PEupxC9hUfHesd\nFRFGl8RYuiTFcVavZLokxdE12XajnJKgAa4ah4a7Ul5qKMgjwoRebVswqW9bBnSIp3vr5nRJjuOU\n+GgddaL8TsNdqVoYY9hVWMLa7AN1BnnP6iDvmMDADgn0addCx3yrJkPDXYW8o0FeyNqcAz8K8nCn\nRT6xbxsGdmypQa4Cgoa7CimeQb4up5A1OYV1B3mHBAZ0SKDvKfEa5CrgaLiroFUzyNc6N88g79mm\nORP6tGFQRw1yFVw03FXQKC6t4NvNBazaafvJ1+UUUlBLkA90+sg1yFUw03BXAW134RHmpecyN30v\nizcXUFZR9UOQj3eCfECHBPppkKsQo+GuAooxhvW7ipizYS/zMvayLqcIgM6Jsfx0RBcm9m3D0C6t\nNMhVyNNwV01eSXkl324uYG76Xual57KnqAQRGNq5FfdN7s2kvm3p0aa5zheulAcNd9Uk5R8qJS0j\nl7kb9vJ1Vj6HyyqJjQpndM9kftW3F+P6tCG5eTO3y1SqydJwV02CMYZNuYeYm76XuRv28v3OAxgD\n7eKjuWRoByb0bcsZ3ZK0u0UpL2m4K9eUV1axfOs+5jjdLTv2HQZgYIcE7p3Qk4l929K/fbx2tyh1\nAjTclV8VHi5nQWYuc9NzWbAxl4MlFURFhDGyexI/G9ONCX3a0i4h2u0ylQp4Gu6q0W3LL7bdLel7\nWb5tP5VVhuTmUUwZ0I4JfdsyumeyrgyklI/p/yjlc5VVhu937GeuM/48K/cQAL3btuBnZ3VjYr+2\nDOnYUmdKVKoRabgrnygurWDRpjzmbMhl/sZc9hWXEREmjOiWyDUjOjOxb1s6Jca6XaZSIUPDXZ2U\ngyXlvPLNNl5YtIWDJRUkxEQyrndrJvRty5jerYmPjnS7RKVCkoa7OiGHyyp4bfF2nlu4mQOHy5nU\nry03jezKsJRWuranUk2Ahrs6LiXllby5ZDvTv9pM/qEyxvZuzS8n9WJQx5Zul6aU8qDhrrxSWlHJ\nO8t38nRaFrkHSxnVI5lfTOrFaV1auV2aUqoWGu6qXuWVVcxamc1T8zaxq7CE4SmJPHn1qZzeLcnt\n0pRS9dBwV7WqqKzig1W7eHLeJnbsO8yQTi159LJBjOqRrFeMKhUANNzVMSqrDJ+s2cUTczexJb+Y\nAR3iefmGVMb1bqOhrlQA8SrcRWQy8AQQDrxojHmkxvYE4E2gs/Oa/zTGvOLjWlUjqqoyfLF+D4/P\nzSRz7yH6tGvBc9eextn92mqoKxWAGgx3EQkHngEmAdnAchH5yBizwWO3O4ENxpgLRKQ1sFFE3jLG\nlDVK1cpnjDHMTc/l8TmZbNhdRPfWcTw99VTOHXCKXkGqVADzpuU+HMgyxmwBEJGZwIWAZ7gboIXY\nJl5zYB9Q4eNalQ8ZY/gqM4/H52SyOruQlKRYHr9yMD8Z3IFwDXWlAp434d4B2OnxOBsYUWOfp4GP\ngF1AC+BKY0xVzRcSkduA2wA6d+58IvUqH1iclc+/5mSyYvt+OrSM4f9dOohLhnbQi4+UCiK+OqF6\nDrAKGA90B+aIyCJjTJHnTsaY54HnAVJTU42P3lt5afm2fTz25UaWbNlHu/hoHr5oAFekdiIqQkNd\nqWDjTbjnAJ08Hnd0nvN0I/CIMcYAWSKyFegDLPNJleqkrNp5gMe+3MiiTfkkN2/GHy7ox9XDO+uq\nRkoFMW/CfTnQU0S6YkP9KmBqjX12ABOARSLSFugNbPFloer4rcsp5PE5mczLyCUxLorfnduHa09P\nISZKQ12pYNdguBtjKkTkLuAL7FDIl40x60VkmrN9OvAX4FURWQsIcL8xJr8R61b12LjnII/PyeTz\n9XtIiInkN+f05vozU2jeTC9rUCpUePW/3RjzKfBpjeeme9zfBZzt29LU8dqcd4h/z93EJ2t20Twq\ngnsn9OTm0V112l2lQpA25YLA9oJinpi3iQ++zyE6Mpzbx3TntrO60TI2yu3SlFIu0XAPYMYYHp+7\niWfmZxERJtw8qivTxnQnqXkzt0tTSrlMwz2A/XvuJp6ct4mLT+3Ab6f0oU18tNslKaWaCA33APX8\nws08MW8TV6R25JFLBulUAUqpY+jVKwHojSXb+dunGZw/6BT+rsGulKqFhnuAeX9lNg9+sI6Jfdvw\n+JVDdB4YpVStNNwDyGdrd/ObWasZ2SOJp6cOJVLnglFK1UHTIUDMz8jlnpnfM7RzK164LlWnDlBK\n1UvDPQB8u7mAaW+upHe7Frx84zBio/Q8uFKqfhruTdx3O/Zz82vL6ZwYy+s3jdCrTZVSXtFwb8LW\n7yrkhpeX0bpFM966ZQSJcXrFqVLKOxruTVRW7iGue2kZzZtF8NYtI/QCJaXUcdFwb4J2FBzmmheX\nICK8devpdGwV63ZJSqkAo+HexOwuPMI1Ly2htKKKN28ZTtfkOLdLUkoFIA33JiT/UCnXvLiU/cXl\nvH7TcPq0i3e7JKVUgNJwbyIOHC7jpy8uZdeBI7xy4zAGdWzpdklKqQCm4d4EHCqt4PpXlrMlr5gX\nrktlWEqi2yUppQKcXg3jsiNlldz86nLW5RTy7DVDGd2ztdslKaWCgLbcXVRaUcm0N1eybNs+/nXF\nYM7u387tkpRSQULD3SUVlVXcO2MVX2Xm8cglA7lwSAe3S1JKBRENdxdUVRnum7WGz9fv4aHz+3Hl\nsM5ul6SUCjIa7n5mjOHBD9cx+/scfn12L24a1dXtkpRSQUjD3Y+MMfzt03TeWrqD28d2585xPdwu\nSSkVpDTc/eiJeZt4YdFWrj+jC/ed0xsRXUVJKdU4NNz95PmFm/n33E1cdlpH/nBBfw12pVSj0nD3\ngzedBa3PG3QKj16qC1orpRqfhnsjm/1dNg9+uI4Jfdrw+BW6oLVSyj803BvRZ2t38+v3VnNGtySe\nuWYoURH6cSul/EPTppHM32gXtB7SqaUuaK2U8jsN90bw7eYCpr2xkl5tW/DKjcOJa6ZT+Cil/Mur\ncBeRySKyUUSyROSBOvYZKyKrRGS9iHzl2zIDx/c79nOLs6D1GzePICFGF7RWSvlfg01KEQkHngEm\nAdnAchH5yBizwWOflsB/gMnGmB0i0qaxCm7KNuwq4vqXl5Hcohlv6oLWSikXedNyHw5kGWO2GGPK\ngJnAhTX2mQrMNsbsADDG5Pq2zKYvK/cQ1760lLhmEbx58wja6oLWSikXeRPuHYCdHo+znec89QJa\nicgCEVkpItfV9kIicpuIrBCRFXl5eSdWcRN0zILWt4ygU6IuaK2UcpevTqhGAKcB5wHnAA+KSK+a\nOxljnjfGpBpjUlu3Do5FKfYUlnDNS0soKbcLWndr3dztkpRSyquVmHKATh6POzrPecoGCowxxUCx\niCwEBgOZPqmyibILWi9hf3E5b90yQhe0Vko1Gd603JcDPUWkq4hEAVcBH9XY50NglIhEiEgsMAJI\n922pTUvh4XKufWkZOQeO8NL1qQzupAtaK6WajgZb7saYChG5C/gCCAdeNsasF5Fpzvbpxph0Efkc\nWANUAS8aY9Y1ZuFue3r+JjbtPchLNwxjRLckt8tRSqljeHV1jTHmU+DTGs9Nr/H4H8A/fFda0zY3\nPZeRPZIZ0ys4zh0opYKLXqF6ArbkHWJrfjHj+4TkcH6lVADQcD8BaRl2GL+Gu1KqqdJwPwFpGbn0\nattcx7MrpZosDffjVFRSzrKt+xinrXalVBOm4X6cFmXmU1FlmNCnrdulKKVUnTTcj1NaRi4JMZEM\n7azj2pVSTZdONH4cqqoMCzbmMqZXayLCXf69WFkBm76E716DwhxITIFWXSGxKyR2s/cTOkKYLhKi\nVCjScD8Oq7MPUFBcxoS+Lva3F2bDd6/Dd2/AwV3QvC20GwS5GZD5BVSWHd03LBJadrZhn9jVCX/n\nfssuEKkzVyoVrDTcj0NaRi5hgv8vXKqsgKw5sPJV21o3BnpMgHP/H/SaDOHOgiBVlVC0C/ZvhX1b\nYN9W5/5W2LEEyg56vKhAfAcn9FNq/ALoCtEJ/j1GpZRPabgfh3npuZzWpRUtY/20CEdhtm2hf/8G\nFOXYVvqoX8LQa20g1xQWDi072VvXs47dZgwcLvAIfI/wz/wcimtMwRyb9ONunur7ca1BpNEOWyl1\n8jTcvbSnsIQNu4u4f3Kfxn2jqkrYNAdWvnK0ld59PEx+BHpPOdpKP14iEJdsb52G/Xh76UHYv80G\n/r4tR1v8O5fCuvfBVB3dNzKu9hZ/+1MhRk80K9UUaLh7qfqq1Ebrby/MsS30796Aomynlf4LGHpd\n7a10X2vWAtoNtLeaKsrgwI6jgV8d/vmb7C+iylK7X3gU9JgEAy+FXlMgSi/yUsotGu5eSsvYS4eW\nMfRs48PFOKoqIWuu7UvP/Ny2jruPh8l/P7lWuq9FREFyD3urqarKntgtyILML2H9bNj4P9u67z0F\nBl4G3SfY11BK+Y2GuxdKyiv5JquAy1M7Ir7oay7aZVvo371uW+lxbWDkz20rPbHryb++P4WF2SGX\nCR2h21g4+y+wfTGsmwUbPrT/RidA35/YoE8ZrcMzlfIDDXcvfLulgCPllSc3UVhVJWTNs33p1a30\nbuNg8t+g97lNp5V+ssLCoetoezv3n7B5vg349f+13U7N20L/i2HApdBxmJ6YVaqRaLh7IS09l5jI\ncE4/kUU5inY7femvQ+FOO9Jk5L0w9PrAa6Ufr/BI6HW2vZUdtieI182CFa/A0ul2DP6AS+2t7QAN\neqV8SMO9AcYY0jLswhzRkV52J1RVwuY025e+8TMwlbaVfvbDtpUeiv3PUbHQ/yJ7KymEjP/ZUTjf\nPAlfPw7JvW23zYBLIam729UqFfA03BuQufcQOQeOcNf4Wk4m1lS0G75/02ml77Ct9DPvhtOut0MG\nlRWdAEOm2ltxPmz4ANbNhvl/tbf2p9qQ738JJHRwu1qlApKGewPmZewFYFzvOvrbq6qcVvorHq30\nsXD2n6H3eaHZSj8ecckw7BZ7K8yxo23WzoIv/w++fBC6nAkDLoF+F9l9lVJeEWOMK2+cmppqVqxY\n4cp7H4/Lpy/mcFkl/7tn9LEbSg/B0mdtK/3ADohNhlN/ake8aLfCySvYbLtt1s6C/I0g4dB9HAy4\nDPqcB9HxbleolCtEZKUxJrWh/bTlXo/9xWWs3L6fO8fV0iUz94+w/AXoOgYm/gn6nK+tdF9K6g5j\n7oOzfgN719mgX/c+fDANwpvZk7QDLoNe50BkjNvVKtXkaLjXY+GmPKpMLWulVpTC2vdg4OVw6Yvu\nFBcqRI5eOTvhD5C93Lbm1/8X0j+GqOa2JT/gMtuyD5YhpUqdJA33esxLzyUpLorBHWvMl5L5BZQc\ngMFXuVNYqBKBTsPtbfLfYdsiG/TpH8GadyCmFfS7EFJvglMGu12tUq7ScK9DRWUVCzbmMqlfO8LC\naoy/Xj0TmreDrmNdqU1hL5bqNtbezvsXbJ5ng37Nu3YIaspoO1KpxyR7Fa1SIUbDvQ7f7ThAUUnF\njycKKy6wF+OM+BmE68fXJERE2Xlsek+xY+hXvmYvknr7CkjuBWfcCYOu0sVJVEjRJk0d5mXsJSJM\nGN2zxvC79bOhqhwGX+1OYap+0Qkw8h64dzVc8qI92frxvfB4f1jwiB1Xr1QI0HCvQ1p6LsO7JtIi\nusYJutUz7aXy7Qa4U5jyTngkDLocbvsKrv8EOqbCgr/bkP/4XsjLdLtCpRqVhnstdu47zKbcQz8e\nJZO/CXJW6InUQCJiJzGb+g7cudx+71bNgGeGwdtXwtZFdkEUpYKMhnstji7M0fbYDatngoTZIZAq\n8LTuBRc8Ab9YD2MesMMqXzsfnh8Da96DynK3K1TKZzTcazEvI5euyXF0TY47+mRVlR2J0W0ctGjn\nXnHq5DVvDeN+a0P+gieg/AjMvgWeGGInMispdLtCpU6ahnsNxaUVLNlc8OMumR2L7WRgeiI1eETG\nwGk3wB1LYeq7dgrmOQ/Cv/rD57+z00ooFaC8CncRmSwiG0UkS0QeqGe/YSJSISKX+a5E//omK5+y\nyiom1Az31TOPXg2pgktYmJ3G4IZP7AnY3pPtUMonhsB7N0LOSrcrVOq4NRjuIhIOPANMAfoBV4tI\nvzr2exT40tdF+lNaRi7Nm0WQmpJ49MnyI7D+A3v1oy76HNzaD7FTSvx8DZxxh13j9oXx8PIUOwd9\nVZXbFSrlFW9a7sOBLGPMFmNMGTATuLCW/e4G3gdyfVifX1UvzHFWr2SiIjw+moz/QdlBHSUTShI6\n2sVVfrEezvkbFGbDzKnwdCosf9GuLKVUE+ZNuHcAdno8znae+4GIdAAuBp6t74VE5DYRWSEiK/Ly\n8o631ka3flcRuQdLGd+nxiiZNe9AfEfoMsqdwpR7ouPtFa73fA+XvWwvkvrfr+x4+bSH4VDAtmVU\nkPPVCdV/A/cbY+r9m9UY87wxJtUYk9q6dWsfvbXvpGXkIgJje3vUdijXLmw96AqdoySUhUfY1aFu\nTYMbP4M/QwdGAAAUBklEQVTOZ8DCf9qQ//BOyE13u0KljuHN5Cg5QCePxx2d5zylAjPFLnCcDJwr\nIhXGmA98UqWfzMvIZXDHliQ3b3b0ybWz7OpK2iWjwF4U1eVMe8vPgiX/gVVv2+UVe0yEM+6yk5np\nYt/KZd40RZcDPUWkq4hEAVcBH3nuYIzpaoxJMcakALOAOwIt2PMOlrJ654FaRsnMsGt6tu7tTmGq\n6UruAef/y/bLj/s/2L0G3rgIpo+2V8HqRVHKRQ2GuzGmArgL+AJIB941xqwXkWkiMq2xC/SXBRtt\n3+k4z3DfuwH2rNGx7ap+cUkw5jfw87Xwk6ftxHIfTIOnhzkhX+F2hSoEeTVnrTHmU+DTGs9Nr2Pf\nG06+LP9Ly8ilbXwz+rf3WJtzzUwIc/palWpIZDQMvdaupZv5Ocz/mw35hf+AsQ/Yn6OwcLerVCFC\nzxACZRVVLNqUz/g+bZDqvtKqSjvdQI9JEJdc/wso5UnEzi3/s4Vw5Vv2StjZt8J/TrfrwOpYeeUH\nGu7A8m37OFRacewQyK1fwcHdeiJVnTgR6Hs+/GwRXP6anXRu1k3w7Jmw4UMNedWoNNyxa6VGRYQx\nskfS0SdXvwPNEqDXZPcKU8EhLAz6XwS3L4ZLX4KqCnj3OnjuLEj/RKccVo1Cwx1Iy9jLGd2SiI1y\nTkGUHrKLLg+4WJdmU74TFg4DL4M7l8IlL0D5YXjnGjvl8MbPNeSVT4V8uG/JO8S2gsPHrpWa8Yn9\njzdIu2RUIwgLtxfF3bkMLnrWTjE840o7h82muRryyidCPtyrF+YY19sj3FfPgJZdoPPpLlWlQkJ4\nBAyZCnetgJ88Zdd3fetSeOls2DxfQ16dlJAP93npufRq25xOic5sj4U5sOUrO7ZdrzJU/hAeCUOv\ng7tXwvmPQ1GOvRjqlXPtMoBKnYCQDveiknKWb9t37CiZte8Bxv7ZrJQ/RURB6k12krJz/wn7t9pl\nAF89H7Yvdrs6FWBCOtwXZeZTUWWO9rcbYxfl6DQCkrq7W5wKXRHNYPitcM8qmPwo5G2EV6bA6xfC\nzmVuV6cCREiHe1pGLgkxkZzaqaV9Ys8ayEvXse2qaYiMhtOnwb2r4ey/wp518NIkePNSyNbVoVT9\nQjbcK6sMCzbmMrZ3ayLCnY9h9UwIj4L+F7tbnFKeomLhzLvs6lAT/wQ538GL4+HtK2HXKrerU01U\nyIb76uwDFBSXHV0Iu7LC9rf3mgwxrdwtTqnaRMXBqJ/bkJ/wEOxYYsfIz5hqZ6RUykPIhvv8jFzC\nBMb0chbm2JwGxXnaJaOavmYtYPSvbMiP+z1s+xqeGw3vXGtnMlWKEA73eem5pHZJpGVslH1i9QyI\nSbQThSkVCKITYMx9NuTH3G/Hxj97Jrx3oz0Jq0JaSIb77sIjbNhdxPjqUTIlhbDxU3tpeESUu8Up\ndbxiWsK439mQH/1LyPwCnhkB799qV4tSISkkw31+hl2c+4f+9g0fQkWJTjegAltsou2L//laGHmv\nnUbjmeHw8b1wcI/b1Sk/C8lwT8vYS8dWMfRs09w+sXomJPWEDkPdLUwpX4hLgkl/gnvX2PHy378F\nT54KaQ9DSZHb1Sk/CblwLymv5OusfCZUL8yxfzts/wYGX6nTDajg0rw1THkU7lpmR4Et/IcN+aXP\nQUWZ29WpRhZy4f7tlgJKyquOrpW65l3776Ar3StKqcaU2A0ufwVunQ9t+sJn99numnXv6+RkQSzk\nwj0tPZeYyHBO75bkTDcwA1JGQ8vObpemVOPqMBSu/xiumQWRsXZVqBfGwdaFblemGkFIhbsxhrSM\nXEb1TCY6MhxyVsK+zdpqV6FDBHpOgmmL7Fzyh/LgtQvgrcth73q3q1M+FFLhnrn3EDkHjhwdJbN6\nBkREQ78L3S1MKX8LC7dzyd+9Aib9GXYuhWdHwgd3QGG229UpHwipcJ+XsRdwFuaoKLN9jn3Oh+h4\nlytTyiWRMXbY5D2r7Pw1a2fBU6fBnIfgyH63q1MnIaTCfX5GLv3bx9MuIRo2fWl/eHW6AaXsGPmz\nH7Yt+f4XwzdPwhNDYPFTUF7idnXqBIRMuO8vLmPl9v1M8OySiWsD3ca5W5hSTUnLznDxdNsn3+E0\n+PL/4OlUey1IVZXb1anjEDLh/lVmHlUGxvdtC4f32Uu0B15u17FUSh2r3UC4djZc96Ft1f/3Z/Dc\nWZA1z+3KlJdCJtzTMnJJbh7FoA4JsH42VJVrl4xSDek2Fm5dAJe+BKVF8OYldkUonUe+yQuJcK+o\nrHIW5mhDWJjA6negTX/bOlFK1S8szE6qd9dymPyInTv++TEw62bYv83t6lQdQiLcV27fT1FJhe1v\nL9gM2ct0ugGljldEMzj9drh3lZ1PPuN/8FQqfP5bKC5wuzpVQ0iEe9rGXCLDhVE9k+2JIQmDgVe4\nXZZSgSk6wc4+ec93MORqWDodnhwCix6DssNuV6ccXoW7iEwWkY0ikiUiD9Sy/RoRWSMia0VksYgM\n9n2pJy4tPZfhXRNpERUOa2ZC1zEQf4rbZSkV2OLbw0+egtu/hS4jYd6f7Rj57163y1YqVzUY7iIS\nDjwDTAH6AVeLSL8au20FxhhjBgJ/AZ73daEnaue+w2zKPcT4Pm1h5xI4sAMGX+12WUoFjzZ9YOpM\nuPEzSOgAH90N00fCxs90YjIXedNyHw5kGWO2GGPKgJnAMdfrG2MWG2OqL2dbAnT0bZknLi0jF3AW\n5lg9EyLjoO/5LlelVBDqcibcPAeueB0qy2HGVfDKuZC9wu3KQpI34d4B2OnxONt5ri43A5/VtkFE\nbhORFSKyIi8vz/sqT8K8jFy6JcfRNSEM1n8A/X5iV5FXSvmeiJ2r6c6lcN5jUJAFL06Ad6+DfVvd\nri6k+PSEqoiMw4b7/bVtN8Y8b4xJNcaktm7d2pdvXavi0gqWbC6wrfaNn0FpoY5tV8ofwiNh2C1w\nz/cw9rewaY6dQ37OQ3bNYtXovAn3HKCTx+OOznPHEJFBwIvAhcaYJjEu6pusfMoqq2y4r3kH4jvY\nuduVUv7RrDmMfQDuXgkDLoNvnoAnh8KKl/WkayPzJtyXAz1FpKuIRAFXAR957iAinYHZwLXGmEzf\nl3li0jJyadEsgtTWlbblMPByO9WpUsq/4tvDxc/CbQsguRd88gt4bjRsTnO7sqDVYLgbYyqAu4Av\ngHTgXWPMehGZJiLTnN0eApKA/4jIKhFx/QxK9cIcZ/VqTVT6bDCV2iWjlNvanwo3fmpPupYVwxsX\nw9tXQv4mtysLOl7NmmWM+RT4tMZz0z3u3wLc4tvSTs76XUXkHiy1a6WumAmnDLbrRyql3FV90rXn\nOfYCqIX/hP+cbvvox9xvJypTJy1or1Cdl56LCExI2ge7V+nYdqWamshoGPVze9L11Gth2fPw5Kmw\nZLodSqlOStCGe9rGXIZ0akmrrNkg4fZkjlKq6WneGi74N0z7GtoPgc/vh/+cARs/14ugTkJQhnve\nwVJW7zzAhF5JsOZd6DHR/gAppZqutv3h2g/g6nfs4xlXwhsX6cLdJygow33+RntV6vkJm6EoR0+k\nKhUoRKD3ZLjjW5j8qJ03fvoo+PjncMg/Fz4Gi+AM94xc2sVH0yX7Y2gWD72nuF2SUup4hEfC6dNs\nf/zw2+D7N+CpofD1v6Gi1O3qAkLQhXtZRRULM/M4u1cLZMNH0P8iu8K7UirwxCbClEfhjiV27pq5\nf4Cnh9mpRLQ/vl5BF+7Ltu6juKySK+JWQ3mxjpJRKhgk94Sp79g++ag4eO96OynZru/drqzJCrpw\nT8vIJSoijL65/7MruXc63e2SlFK+0n0c/GwRnP845GfC8+Pgv7dD0W63K2tygircjTHMy9jLuV0M\n4du+gkFX2fUflVLBIzwCUm+yK0GNvAfWzbL98Qse1ZWgPARV8m3JL2Z7wWGujVsOpkpHySgVzKIT\nYNKf4c5l0HMSLPgbPJ1qhz9XVbldneuCKtznOwtzDCz4DDoOg6TuLleklGp0iV3tXDU3fApxyTD7\nVnhpIuxY6nZlrgqqcJ+XnsuU5HyiCtK11a5UqEkZCbcugIuehaJd8PLZ8N6NdmnNEBQ04V5UUs7y\nbfu4scUSCIuE/pe4XZJSyt/CwmDIVDt//Jj77SI9T6XaxbtLD7pdnV8FTbgvyszHVFVw6v4vodc5\nOrOcUqEsKg7G/Q7uXmFnoFz0GDx1Gnz3Rsj0xwdNuM/L2Ms5MelEluTr2HallJXQES59AW6ZBy27\nwEd32TVddy53u7JGFxThXlll+GpjHje3WAoxraDn2W6XpJRqSjqmws1fwiUvwMHd9oTrf6fBwT1u\nV9ZogiLcV2cfoLT4AIOLv4EBl0JElNslKaWaGhEYdAXctQJG/RLWvW+7ar5+PCjnqwmKcE9Lz+Xc\niOVEVJbYC5eUUqouzZrDxD/AnUuh61kw9492Jaggmz8+OMI9I5frYr+FxO72zy+llGpIYje4egb8\n9H0Ii7Dzx791edCs5xrw4b678AiFuzczoGyNPZEq4nZJSqlA0mMi3L4Yzvkb7FxqW/Ff/B5Kityu\n7KQEfLinZeRyYfg39sGgK9wtRikVmMIj4Yw77fj4wVfDt8/Y/vjv3wzYoZMBH+7z0/dyZdQ3mC5n\nQqsubpejlApkzdvAhU/DrWnQKgU+vNOOrMle4XZlxy2gw72kvJIDm5fRxeQgeiJVKeUrHYbCTV/A\nxc9BYY4dG//f2wNq6GRAh/u3mws433xFVViUXXFJKaV8JSzMzlF19woY9QtnauHT4JsnoKLM7eoa\nFNDh/tWGHH4S/i2mz3l2+k+llPK1Zi1g4h/tUn8po2HOQ/aka+aXbldWr4ANd2MMR9K/IFEOEj5E\npxtQSjWypO4wdSZc8z5IGLx9uTN0MsvtymoVsOGeufcQY0rmURKVCN3Hu12OUipU9HSGTp79MGz/\n1rbiv3ywyQ2dDNhw/3rtJiaEfUflgMvsMCallPKXiCg482671N/gK2Hxk3YVqFVvN5mhkwEb7hVr\nZtNMKohLvcbtUpRSoap5G7jwGTt0smVn+OB2Z+jkSrcrC8xw319cxmmFX5Af2w1OGex2OUqpUNfh\nNLjpS7hoOhRmw4vj4YM74OBe10oKyHBf/v1KUsMyKet3hU43oJRqGsLCYMjV9irXkT+3C3U/dRp8\n86QrQye9CncRmSwiG0UkS0QeqGW7iMiTzvY1IjLU96UeVfH9DKoQ2o26tjHfRimljl+zFjDpT3bW\nyZSRMOdBePYM2DTHr2U0GO4iEg48A0wB+gFXi0i/GrtNAXo6t9uAZ31c5w8qKioZWPA5WXGnEday\nY2O9jVJKnZyk7jD1Hbhmln381mXw1hVQsNkvb+9Ny304kGWM2WKMKQNmAhfW2OdC4HVjLQFaisgp\nPq4VgI0r5tKJvRzpe1ljvLxSSvlWz0lw+7cw6S+wfTE8M8JOTNbIvAn3DsBOj8fZznPHuw8icpuI\nrBCRFXl5ecdbKwBhYWGsiR5Gt7N0LhmlVICIiIKR99j++EFX2EnJGvstG/0dPBhjngeeB0hNTT2h\nJU/6Dp8Ewyf5tC6llPKLFm3hov/45a28abnnAJ08Hnd0njvefZRSSvmJN+G+HOgpIl1FJAq4Cvio\nxj4fAdc5o2ZOBwqNMbt9XKtSSikvNdgtY4ypEJG7gC+AcOBlY8x6EZnmbJ8OfAqcC2QBh4EbG69k\npZRSDfGqz90Y8yk2wD2fm+5x3wB3+rY0pZRSJyogr1BVSilVPw13pZQKQhruSikVhDTclVIqCIk9\nF+rCG4vkAdtP8MuTgXwflhMI9JhDgx5zaDiZY+5ijGnd0E6uhfvJEJEVxphUt+vwJz3m0KDHHBr8\ncczaLaOUUkFIw10ppYJQoIb7824X4AI95tCgxxwaGv2YA7LPXSmlVP0CteWulFKqHhruSikVhJp0\nuDe1hbn9wYtjvsY51rUislhEBrtRpy81dMwe+w0TkQoRCfg1Fr05ZhEZKyKrRGS9iHzl7xp9zYuf\n7QQR+VhEVjvHHNCzy4rIyyKSKyLr6tjeuPlljGmSN+z0wpuBbkAUsBroV2Ofc4HPAAFOB5a6Xbcf\njvlMoJVzf0ooHLPHfmnY2Ukvc7tuP3yfWwIbgM7O4zZu1+2HY/4d8KhzvzWwD4hyu/aTOOazgKHA\nujq2N2p+NeWWe5NamNtPGjxmY8xiY8x+5+ES7KpXgcyb7zPA3cD7QK4/i2sk3hzzVGC2MWYHgDEm\n0I/bm2M2QAsREaA5Ntwr/Fum7xhjFmKPoS6Nml9NOdx9tjB3ADne47kZ+5s/kDV4zCLSAbgYeNaP\ndTUmb77PvYBWIrJARFaKyHV+q65xeHPMTwN9gV3AWuBeY0yVf8pzRaPml18XyFa+IyLjsOE+yu1a\n/ODfwP3GmCrbqAsJEcBpwAQgBvhWRJYYYzLdLatRnQOsAsYD3YE5IrLIGFPkblmBqSmHeyguzO3V\n8YjIIOBFYIoxpsBPtTUWb445FZjpBHsycK6IVBhjPvBPiT7nzTFnAwXGmGKgWEQWAoOBQA13b475\nRuARYzuks0RkK9AHWOafEv2uUfOrKXfLhOLC3A0es4h0BmYD1wZJK67BYzbGdDXGpBhjUoBZwB0B\nHOzg3c/2h8AoEYkQkVhgBJDu5zp9yZtj3oH9SwURaQv0Brb4tUr/atT8arItdxOCC3N7ecwPAUnA\nf5yWbIUJ4Bn1vDzmoOLNMRtj0kXkc2ANUAW8aIypdUhdIPDy+/wX4FURWYsdQXK/MSZgpwIWkRnA\nWCBZRLKBPwCR4J/80ukHlFIqCDXlbhmllFInSMNdKaWCkIa7UkoFIQ13pZQKQhruSikVhDTclVIq\nCGm4K6VUEPr/YW0T1SM0TXwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5bf7f53cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(\"Normalized spline & difference curve\");\n",
    "plt.plot(xsn,ysn);\n",
    "plt.plot(xd,yd);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 4: Identify local maxima  \n",
    "of the difference curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.signal import argrelextrema"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# local maxima for knees\n",
    "xmx_idx = argrelextrema(yd, np.greater)[0]\n",
    "xmx = xd[xmx_idx]\n",
    "ymx = yd[xmx_idx]\n",
    "#Dmx = np.stack((xmx,ymx))\n",
    "\n",
    "# local minima\n",
    "xmn_idx = argrelextrema(yd, np.less)[0]\n",
    "xmn = xd[xmn_idx]\n",
    "ymn = yd[xmn_idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VVW6x/Hvm96AJCR0QmgiSO8KKrahjYpt7Ag6Ktbx\nzp3iOI4zdyxXp6kzqAyDXUf0WhhExjJWAiJdEBASaug5KZBG6rp/rA0cQkIOcJJ9yvt5nvNwzt47\nZ787J/mxs/baa4kxBqWUUqElwu0ClFJK+Z+Gu1JKhSANd6WUCkEa7kopFYI03JVSKgRpuCulVAjS\ncA8wIrJVRC5swvcfIyI7mur9T6COGSLym2bYT4aIlIhIpJ/ez4hID+f5UccgIneIyF5nf61FZJSI\nZDuvJ/lj/0r5KsrtAlR4MsZMa6b9bAeSmui9Dx+DiEQDfwFGGmO+dZb9HphujHm6Kfav1PHombtS\n/tEWiAPWei3rUue1z0QkpE68Qu14goGGewATkVgReUpEdjmPp0Qk1mv9pSKySkQOiMgmERnnLJ8q\nIutFpFhENovI7SewTyMidzrNCcUi8rCIdBeRRc5+3hKRGGfbFBGZJyJ5IlLoPO/krEsVkR0icrHz\nOklEckRksvP6JRF5xHk+xtn2FyKyT0R2i8gkEZkgIhtFpEBEHvCqcbiIfC0iRc620w/VVM/xZDrH\nFOW8/sI5poXO8X0sImnH+X783NnHLhG5uc66l0TkERE5DdjgLC4Skc9EZBPQDXjfaZaJFZFWIvK8\n8347na+NdN5rilPTkyKSD/zOWX6z81kWishHItKlzmc1zfmsikTkGRERr/W3ev0crBORwc7yDiLy\njvO5bRGRe49z/PEi8mcR2SYi+0Uky1l2TPOeeDUpisjvRORtEXlNRA4AD4hIuYikem0/SEQ8Yv/q\nOe6xqpNgjNFHAD2ArcCFzvPfA4uBNkA6sAh42Fk3HNgPXIT9T7ojcLqzbiLQHRDgXKAMGOysGwPs\nOM7+DfAvoCVwBlABfIoNqlbAOuAmZ9vWwBVAAtAC+D9gjtd7/QDY49T/D+Btr3UvAY941VQNPARE\nA7cCecA/nfc9AygHujrbDwFGYpsVM4H1wH0NHE+mc0xRzusvgE3AaUC88/rxBr52HLAX6AskOvUY\noEc9x3DUfup+ls7r94C/O+/VBlgC3O6sm+J8D+5xjiseuBTIAXo7yx4EFtX5rOYByUCG8z0b56y7\nCtgJDHN+Dnpg/5KIAJY73+sY53PdDIxt4HvwjPM96ghEAmcBsdTzc8TRP7u/A6qASc4+44HPgFu9\ntv8jMMN5ftxj1cdJZInbBeijzgdy9C/IJmCC17qxwFbn+d+BJ318zznAT5znx/xS1tnWAKO8Xi8H\nfun1+s/AUw187UCgsM6yvwFrnKBp7bXcOxjHYMM70nndwqljRJ06JjWw3/uA9xpYl8mx4f6g1/o7\ngQ8b+NoX8Ap+7H8IJxXu2GabCiDea/21wOfO8ynA9jr7/zdwi9frCOx/1F28PqvRXuvfAu53nn90\n6DOv854j6tnPr4AX69k2wvlcBtSz7pifI44N96/qrP8x8JnzXIBc4BxfjlUfJ/7QZpnA1gHY5vV6\nm7MMoDM2/I8hIuNFZLHTnFEETAAabHqox16v5+X1vE5y9pMgIn93/mQ/AHwFJMvRPVNmYs98XzLG\n5B9nn/nGmBqvfdRXx6H9nuY0Ae1x9vvYCR7fHq/nZTR8wbUDNoAO2dbAdr7ogv2rZLfThFKE/Q+6\njdc2ufV8zdNe2xdgQ7Gj1zYNHUtDPx9dgA6H3tN53wew//nUlYa9jlDvz5kP6h7PO8CZItIeOAeo\nBRZ41dXYsaoToOEe2HZhf+gPyXCWgf3F6V73C8S2yb8D/Aloa4xJBuZjf1H87b+BXtgz7JbYX1gO\n7csJ+ZnAK8Cd4nQh9IPngO+Bns5+H6Bpjm83NiQPyTiF98rFnrmnGWOSnUdLY8wZXtvUHaI1F9ts\nk+z1iDfGLPJxf8f8fDjLt9R5zxbGmAn1bOsBDjbwPqXY5jjg8GedXmebo47HGFMIfAxcDVwHzDbO\naTqndqyqHhruge0N4EERSXcu+j0EvOasex6YKiIXiEiEiHQUkdOx7aix2PbXahEZj237bgotsGfU\nRc6Fst/WWf8A9hf8Zmz76ivin/7mLYADQIlzzHf44T3r8xYwRUT6iEgCxx6fz4wxu7HB9mcRael8\nZt1F5NzjfNkM4FcicgaAc0H2Kh93OQv4mYgMEauHc4FyCVAsIr90LoxGikhfERlWT8212KapvzgX\nYSNF5EznBGIjECciE50Log9if+4a809gMnCl89wfx6rqoeEe2B4BlgGrse3WK5xlGGOWAFOBJ7EX\nVr/Etk8WA/dig6kQe4Y0t4nqewp7ocyDvfD74aEVIjIE+Ckw2WlueQIb9Pf7Yb8/wx5XMfZC7Zt+\neM9jGGP+jT3Gz7AX+z47xbecjP3Pdx32s3kbaH+c/b+H/b7NdpqfvgPG+7IjY8z/AY9iA7QYe90l\n1fksfoi9PrIF+9nNwl4sr8/PsD97S7FNJU8AEcaY/djrFbOw11NKAV9ujpsL9AT2GOd+gFM9VlU/\nOfJXkVJKqVChZ+5KKRWCNNyVUioEabgrpVQI0nBXSqkQ5NpgPmlpaSYzM9Ot3SulVFBavny5xxhT\n956CY7gW7pmZmSxbtsyt3SulVFASEZ/ulNZmGaWUCkEa7kopFYI03JVSKgRpuCulVAjScFdKqRDU\naLiLyAtipz77roH1IiJ/FTuF2upDU3kppZRyjy9n7i9hpxtryHjsKG89gduwY20rpZRyUaP93I0x\nX4lI5nE2uRR4xRl0f7GIJItIe2f8aqWUClrGGMoqaygqr6KwtJL95VUUllVSVFZFUVklnVISmDQo\nMCeL8sdNTB05ejqtHc6yY8JdRG7Dnt2TkXEqk9oopdSJOVhVQ1HZ0eFc5IT1/qOWV1FUXklhWRX7\ny6qorKlt8D0v7N02pMPdZ8aYmdhp1xg6dKgOJK+UOimlFdV4SirwlFSQX1J5VCAfDm4nsA+dbR+s\najikY6IiSEmIJiUhhlbx0XRLSyI5IZrkhBiSE6JJSYimVXwMKc6ylIRoWsZHExftj4nFmoY/wn0n\nR88z2clZppRSPjHGUFpZg6e4grySCjzFNrjzSiptiB9+XYGnuJLyqpp63yc6Umwgx9ug7pyaQH8n\nkFs5y2xoR5McH0NKov03PiZwQ/pk+SPc5wJ3i8hsYASwX9vblVLGGIorqp1grjx8pm0DvJI8J7AP\nPeo7sxaB1IQY0pJiSW8Ry5CMBNKSYklrEWv/TbLrkp2z7oSYSESaYq704NNouIvIG8AYIE1EdmAn\nCY4GMMbMAOYDE7BzTJZh5/VUSoWwkopqdheVs7OonF1FB9m9v/xwWOeVVB4+066oPjawIwRSE48E\ndte0xMMhfSi405NiSWsRQ2pCDFGRejvOyfClt8y1jaw3wF1+q0gp5arqmlr2FVewyyu8dxWVe70u\n58DB6qO+xga2Deu0pBi6pyU6Z9dHQvxQeKcmxhAZoWfXTc21IX+VUs3PGMOBg9WHw9oG9pHw3r3/\nIHsOHKSm9uj+Dq3io+mQHE+nlHiGd02lQ3I8HZLj6ZgcR4fkeNKTYvUMO8BouCsVQiqra9l74ODh\nM2zv8N69356Fl1QcfdYdHSm0bxVPh+Q4RnRLpWNy/OHXHZPjaZ8cT1KsRkWw0U9MqSBjjGHPgYNs\n2lfKprySw4/NeaXsOXAQU6eTcevEGNonx5HZOpGzuqfR0TnrPhTeaUmxRGgzScjRcFcqQFVU17DV\nU2bDe9+hEC9lc14JpZVHugK2iI2iW5skzuzWms6pCUeFd/tW8SHZzU81TsNdKZcVlFYeE+Cb8krI\nLSjDu+m7Y3I83dITuWpoZ7q3SaJ7eiI90pNIbxGr3f/UMTTclWoGNbWG3IKyI80oXk0qhWVVh7eL\niYqgW1oifTu24tKBHemenkj39CS6pSeSEKO/rsp3+tOilB+VV9aQs6/kqLbwTftK2eIpPWqMktaJ\nMXRPT2Jc33Z0T0+ie5skeqQn0SE5XrsJKr/QcFfqJBlj2JpfxsrthazcXsSK7YV8v6f4cDfCyAgh\nIzWB7umJjOmV7oR4It3SkkhJjHG5ehXqNNyV8lHxwSpW79jPyu2FrNhexMrthYebVBJjIhmYkcwd\n53bnjA4t6dEmiYzWCcRG6cVM5Q4Nd6XqUVtr2OwpORziK7cXsWFv8eFuhj3aJHFh77YM7pLCoIxk\nerZpoc0pKqBouCsF7C+rYtWOIlZsK2RlbhGrthcevsW+ZVwUgzJSGNe3HYMyUhjYOZlW8dEuV6zU\n8Wm4q7BTU2vI3lfMim1FThNLIZvySgE7RsppbVswsX8HBmUkMzgjhW5piXqTjwo6Gu4q5BWUVh51\n0fPb3KLDNwGlJsYwqHMylw/uxKDOyfTvnKy32quQoD/FKuRsyy/lq415h9vLt+aXAbb3Sp/2Lbli\nSKfDZ+UZqQl6A5AKSRruKiRs9ZTywZrdzF+zm7W7DgCQ3iKWwRnJXDs8g0EZKfTr2EpvxVdhQ8Nd\nBa36An1QRjIPTuzND/q0o3NqvJ6Vq7Cl4a6CyvECfXy/9nRMjne5QqUCg4a7CniHAv2D1btZt1sD\nXSlfaLirgLTFU8p8DXSlTpqGuwoYGuhK+Y+Gu3LV8QJ9Qr/2dNBAV+qkaLirZnco0Oet3s16J9AH\na6Ar5Vca7qpZbM4rsWfoa/ZooCvVDDTcVZPZ4inlg9W7NNCVcoGGu/K73IIynvjwe+at3g1ooCvl\nBg135TfFB6t45vNNvLBwCxEC95zfg2uHZ2igK+UCDXd1yqpranlzWS5/+Xgj+aWVXD64Iz8f24v2\nrTTUlXKLhrs6JV9uzOPRD9axcW8Jw7um8uLE3vTvlOx2WUqFPQ13dVI27i3m0Q/W8+XGPLq0TmDG\nDUMYe0ZbHahLqQCh4a5OiKekgic/2cgbS7aTGBvFgxN7M/nMTGKiItwuTSnlRcNd+eRgVQ0vLtzK\nM5/nUF5Vw+QzM7n3gp6kJsa4XZpSqh4+hbuIjAOeBiKBWcaYx+usbwW8BmQ47/knY8yLfq5VucAY\nw7zVu3n839+zs6icC3u34f7xvenRJsnt0pRSx9FouItIJPAMcBGwA1gqInONMeu8NrsLWGeMuVhE\n0oENIvK6MaaySapWzWLF9kIembeOFduL6N2+JX+4sj+jeqS5XZZSyge+nLkPB3KMMZsBRGQ2cCng\nHe4GaCH2aloSUABU+7lW1Ux2FJbxhw83MPfbXaS3iOUPV/TniiGdiIzQi6VKBQtfwr0jkOv1egcw\nos4204G5wC6gBXC1Maa27huJyG3AbQAZGRknU69qQsUHq3jui03MyrI3Id17fg9uP7c7ibF6aUap\nYOOv39qxwCrgfKA78ImILDDGHPDeyBgzE5gJMHToUOOnfatTVF1Ty1vLdvCXTzbgKank8kEd+dnY\nXnpnqVJBzJdw3wl09nrdyVnmbSrwuDHGADkisgU4HVjilypVk/lqYx6PfrCeDXuLGZaZwgtThulN\nSEqFAF/CfSnQU0S6YkP9GuC6OttsBy4AFohIW6AXsNmfhSr/yt5bzKPz1/PFhjwyUhN47vrBjOvb\nTm9CUipENBruxphqEbkb+AjbFfIFY8xaEZnmrJ8BPAy8JCJrAAF+aYzxNGHd6iTll1Tw5H828saS\nXBJiIvn1hN5MPqsLsVGRbpemlPIjn9rcjTHzgfl1ls3wer4L+IF/S1P+dLCqhpcWbeWZz3Ioq6rh\nhhEZ/OTC0/QmJKVClHaDCHHGGOav2cPjH64nt6CcC05vw68m6E1ISoU6DfcQVlJRzY9fXsrizQWc\n3q4Fr90ygtE99SYkpcKBhnuIqq6p5e5/rmDp1kIemdSXa4dn6E1ISoURDfcQZIzht3PX8sWGPB67\nrB/XjdAbxpQKNzpOawj6x4LNvP7Ndqad212DXakwpeEeYj5YvZvH5n/PxP7t+cXYXm6Xo5RyiYZ7\nCFm+rZD/emsVQ7qk8OerBhChbexKhS0N9xCxLb+UW19ZRodWcfxj8lDiovWmJKXCmYZ7CCgsrWTq\ni0sxxvDi1OF6Y5JSSnvLBLuK6hpuf3U5OwrLef3WEXRNS3S7JKVUANBwD2K1tYaf/99qlmwt4K/X\nDmJYZqrbJSmlAoQ2ywSxv3yykbnf7uLnY3txyYAObpejlAogGu5B6q2luUz/PIdrhnXmzjHd3S5H\nKRVgNNyDUFa2hwfeW8PZPdN4eFJfHYNdKXUMDfcgs2FPMXe8tpwebZJ45vrBREfqR6iUOpYmQxDZ\ne+AgU19cQnxMJC9MGUbLuGi3S1JKBSjtLRMkSiuqueXlpRSVV/HW7Wfq5NVKqePSM/cgUFNruPeN\nlazbdYDp1w2ib8dWbpeklApweuYe4Iwx/P79tXz6/T4evvQMzj+9rdslKaWCgJ65B7gXFm7l5a+3\ncevZXbnxzEy3y1FKBQkN9wD24Xd7eOSDdYzv245fje/tdjlKqSCi4R6gVuUWcd+bKxnQKZknrx6o\nw/cqpU6IhnsAyi0o48cvLyW9RSyzbtLhe5VSJ04vqAaY/WVVTHlxCVU1htlThpOWFOt2SUqpIKRn\n7gGksrqW219bxvaCMv5+4xB6tElyuySlVJDSM/cAYYzh/ndWs3hzAU9dPZCR3Vq7XZJSKojpmXuA\nePrTbN5duZOfXnQakwZ1dLscpVSQ03APAG8v38FT/8nmyiGduOf8Hm6Xo5QKARruLlu0ycOv3l3N\nWd1b89hl/XT4XqWUX2i4uyh7bzG3v7qczNaJPHfDEGKi9ONQSvmHT2kiIuNEZIOI5IjI/Q1sM0ZE\nVonIWhH50r9lhp684gqmvrSUuOhIXpw6jFbxOnyvUsp/Gu0tIyKRwDPARcAOYKmIzDXGrPPaJhl4\nFhhnjNkuIm2aquBQUF5Zw49fXkp+SSVv3j6STikJbpeklAoxvpy5DwdyjDGbjTGVwGzg0jrbXAe8\na4zZDmCM2effMkNHTa3hJ7NXsnrnfv567SD6d0p2uySlVAjyJdw7Arler3c4y7ydBqSIyBcislxE\nJtf3RiJym4gsE5FleXl5J1dxkHv0g/V8vG4vD/2wDxf10eF7lVJNw19X8KKAIcBEYCzwGxE5re5G\nxpiZxpihxpih6enpftp18Hhp4RZeWLiFqaMymTqqq9vlKKVCmC93qO4EOnu97uQs87YDyDfGlAKl\nIvIVMADY6JcqQ8B/1u3l9/PWcVGftjw4sY/b5SilQpwvZ+5LgZ4i0lVEYoBrgLl1tvkXMFpEokQk\nARgBrPdvqcFrzY793PPGSvp2bMXT1wwkUofvVUo1sUbP3I0x1SJyN/AREAm8YIxZKyLTnPUzjDHr\nReRDYDVQC8wyxnzXlIUHk4fmfkdKQjSzbhpKQowO56OUano+JY0xZj4wv86yGXVe/xH4o/9KCw1F\nZZWsyi3iJxf0pE2LOLfLUUqFCb0lsol9vSkfY+Dsnmlul6KUCiMa7k1sQY6HpNgo7c+ulGpWGu5N\nLCvbw8hurYmO1G+1Uqr5aOI0oe35ZWwvKNMmGaVUs9Nwb0ILcuxduKM13JVSzUzDvQktzPHQvlUc\n3dIS3S5FKRVmgrLT9ZgxY9wuoVEGIXfoXSQU5HDeeQ83+/5jImq5qG0BE9rlU2OE3PI4tpfFklsW\nR25ZLLsOxlJj9GYqpdzwxRdfNPk+gjLcg0FlYltqo+KJ27+1WfebGlPFpA55XNIhn+SYajaVxFFS\nHcmZrfczsX314e1qDOwqjyW3LPZw8O9wnhdURgEa/EoFs6AM9+b4X+9UPfN5Dn/8aAMfvvw0aUmx\nTb/D3ath8bOw5m2orYZe42HknXTPHA2Hpu4rL4L8HPBkE5mfQ+f8bDp7cqBgE1QfPPJesS2hdQ9I\n62n/PfQ8tTvE6NjzSgWDoAz3YJCV7aFP+5ZNG+y1NbDxIxvqWxdAdCIMnQojpkHr7sduH58MnYba\nx1HvUwv7c23wO+FPfjZsWwSr3zx621ad7Xu37nkk/NN6QstOEKGXcJQKFBruTaC8sobl2wqZMiqz\naXZQUQKrXofFz0HhFhu4Fz0MgyfbAD9RERGQ0sU+elxw9LrKMntm78k+OvxXvwkVB45sFxVnz+zT\nehwJ/vRe0G6Ahr5SLtBwbwLfbMmnsqaW0T383AWyKBeW/B2WvwIV+6HTMLjgIeh9CUQ20UcZkwDt\n+tmHN2OgZJ8T+NlHwn/Pd7B+Hpgau11yFxh4PQy8FpIzmqZGpdQxNNybQFa2h5ioCIZ3TfXPG+Yu\nhcXPwDpnpOU+l8DIu6DzMP+8/8kQgRZt7SNz1NHraqqgcCvsXA6r/glfPAZf/C90GwODboDTJ0J0\nvAtFKxU+NNybQFaOh2GZKcRFR578m9RUw/p/wdfPws5lENsKzrwLht8GyZ0b/3o3RUbbZpm0njDg\nGijcBt++AStfh3dugbhW0PdKG/QdBh254KuU8hsNdz/bV3yQ7/cU84txvU7uDcqLYMXL8M1MOLAD\nUrvB+D/CwOsgNsm/xTaXlC4w5n445xf2wu/K1+w1g2XPQ5s+NuT7Xw2JeievUv6i4e5ni3LyATi7\nxwnOEZu/yV4gXfVPqCqFzLNh4p+g59jQuSAZEQHdzrWP8j/C2ndt0H/0AHzyEJw2DgbdCD0ubLpr\nCEqFCf0N8rMF2R5SEqI5o0PLxjc2xp7Jfv0sbPzQNmf0vRJG3gHt+zd9sW6KT4ahN9vHvvU25L+d\nDd/Pg6S2tjln4A2Qfsw860opH2i4+5ExhqycPM7qkUbE8eZJra6A796x/dP3rIGENDj3FzD0FnuB\nMty06Q1jH4ULfwfZH9ugXzQdFj4NnYbbZpszLoM4H/7DVEoBGu5+tSmvhL0HKhruAlnqgaXPw9JZ\nULoP0nvDJX+Dfj+CaJ2Cj8ho25Pm9IlQvNf2pV/5Grx/L3x4P/S51AZ9l1F6EVapRmi4+9GCbA/A\nseFeuBW++hOsfgtqKqDHRXDmndDtPA2phrRoC6PuhbPugR3LYNVrsOYd2+smJdM22Qy8Flp1crtS\npQKShrsfZWV7yGydQOfUOuOvvHmjvcln4HW2PT39JHvShCMR25+/8zAY+7+w/n1Y+Sp8/gh8/ih0\nP8+ezfeaqH/9KOVFw91PqmpqWbw5n8sGdzx6hScH9qyGcU/AyGnuFBcqYhJgwNX2UbDlSN/5t2+G\nuGTod5UN+vYD9C8iFfZCpI+d+1blFlFaWcPoul0g171n/+1zSfMXFcpSu8J5D8B9q+HG92z3yRWv\nwMxzYcZoe59AVbnbVSrlGg13P1mQ7SFC4MzurY9esXYOdB4JLTu4U1ioi4iE7ufDlc/DzzbAhD9B\nRBT8++fwVH/b66ayzO0qlWp2Gu5+kpWdR/9OybSKjz6y0JMDe7+DMya5V1g4iU+B4bfC7V/ClPm2\ni+XHv4an+9tulRUlbleoVLPRcPeDAwer+HbHfs6uOxH2oSaZ3tok0+wyR8FNc+Hmj+yIlp88ZEN+\nwV+gotjt6pRqchrufvD1pnxqas2xXSDX/gs6j4BWHev/QtX0MkbaNvlbPoEOg+HT/4Gn+sFXf4SD\nBxr/eqWClIa7HyzM8ZAQE8mgjJQjC/M3wd410EebZAJC5+Fww9vw48/sXa+fPQJP9YUvnrCDtSkV\nYjTc/SAr28OIrqnERHl9O9ce6iVzqTtFqfp1GgLXvwW3fWHvdP3iMXvh9fPHoLzQ7eqU8hsN91O0\ns6iczZ5SRves2wVyjj1D1CaZwNRhEFz7Btz+FXQ9G758Ap7sB58+DGUFblen1CnTcD9FWdl5AEdf\nTM3fZAcEO+Myl6pSPms/AK55HaYthB7nw4I/2Tb5//wOSvPdrk6pk+ZTuIvIOBHZICI5InL/cbYb\nJiLVInKl/0oMbFk5+bRpEUvPNl4TaaybY//VJpng0a4v/OgVuONr6PkDyHrKhvzHv4GSPLerU+qE\nNRruIhIJPAOMB/oA14pInwa2ewL42N9FBqraWsPCHA+je6Qh3re7r31Pm2SCVds+cNWLcNc3cPoE\n+Hq6DfmPfm1HqlQqSPhy5j4cyDHGbDbGVAKzgfpOSe8B3gH2+bG+gLZu9wEKSisZXW+TjPaSCWrp\nveCKWXDXEvsX2OJnbT/5f98PB3a7XZ1SjfIl3DsCuV6vdzjLDhORjsBlwHPHeyMRuU1ElonIsry8\n4P9TNyunniF+tUkmtKT1hMv/Dncvg75XwJKZ8PQAmP9z2L/T7eqUapC/Lqg+BfzSGFN7vI2MMTON\nMUONMUPT009wjtEAtDDHQ6+2LWjT0muo2bVzoNMwHWc81LTuDpOehXuWQ/8fwbIX4K8DYd5PoSi3\n8a9Xqpn5Eu47gc5erzs5y7wNBWaLyFbgSuBZEQnpdomDVTUs2VLAKO+z9oLNdnhfvXEpdKV2hUun\nwz0r7Pj8K16Bvw6C9++Dwm1uV6fUYb6E+1Kgp4h0FZEY4BpgrvcGxpiuxphMY0wm8DZwpzFmjt+r\nDSDLthZSUV17dBfItdokEzZSusDFT8O9K2HwZFj1OvxtMMy9R0NeBYRGw90YUw3cDXwErAfeMsas\nFZFpIhK2s08syMkjOlIY0S31yMJ1c6DjUEju3PAXqtCS3Bl++Be4dxUMvRm+fROmD4WPH9RhDZSr\nfJqJyRgzH5hfZ9mMBradcuplBb6sbA+DM1JIiHG+hQVbYPe38INH3C1MuaNVR5jwRxh1n53+b9F0\nO7n3uffDsFvs5N9KNSO9Q/UkFJRWsnbXgaObZLSXjAIb8pOetcMatOsHH/4SnhkB6+eBMW5Xp8KI\nhvtJWOh0gTzqYuraQ00yGS5VpQJK+/4weS5c95adGerN6+GlH8LOFW5XpsKEhvtJyMr20DIuiv6d\nku2Cgi2we5XeuKSOJgKnjYU7FsHEP0Pe9/CP8+Dd22D/DrerUyFOw/0EGWPIyvFwVvc0IiOcIQe0\nSUYdT2QUDPux7Vkz+r/sX3l/GwL/+R+dMEQ1GQ33E7Q1v4ydReVHDzmwdg50HKJNMur44lrChb+D\ne5bZqRfQlognAAAQzUlEQVSz/mK7Ty59Hmqq3a5OhRgN9xN0aIjfw0MOHGqS0RuXlK+SM+CKf8Ct\nn0HrnvDBT2HGKNj4sV50VX6j4X6CFmR76JQST5fWCXbBun/Zf7VJRp2ojkNg6ny4+jWoqYR/XgWv\nTrIDzyl1ijTcT0B1TS1fb8rn7J5eQ/yum2MnXk7p4m5xKjiJQO+L4c5vYNzj9l6JGWfDnLt09El1\nSjTcT8C3O/ZTXFHN6B7OoGeFW2HXSu0lo05dVAyMvMNedD3zLlj9pm2P//x/obLU7epUENJwPwEL\nczyIwFndW9sFh5tkNNyVn8SnwNhH4e6ldkaoLx+Hvw6GFa9CbY3b1akgouF+ArKyPfTt0IqUxBi7\nYK02yagmktoVfvQy3PyxHb9m7t3w93Ng02duV6aChIa7j0oqqlmxvfBIF8jCrbBrhTbJqKaVMQJu\n+QSufBEqiuHVy+C1K2HfercrUwFOw91H32zOp7rWcPahLpDaS0Y1FxHoe7ltqrnoYchdAs+dZceQ\nLwmbWS3VCdJw91FWjoe46AgGd0mxC9bOgQ6DICXT1bpUGImKhVH32ouuw26Fla/aiUK++hNUlbtd\nnQowGu4+ysr2MCwzlbjoSDsZw64VeiFVuSOxNUz4g+0+2fVc+OxhO5zBt7Oh9rgzXaowouHugz37\nD5K9r+TIEL+HmmS0vV25Ka0HXPtPmPIBJKbDe7fD8xfZ7rkq7Gm4+yDLGeL3cP/2dXOg/UBtklGB\nIXM03Po5THoOirbDzPNse3xZgduVKRdpuPtgYY6HtKQYTm/Xwv7y7FyuZ+0qsERE2Am771kGI6bZ\nibv/NhiWvaj948OUhnsjvIf4jYgQvXFJBba4VjD+cZi2ANJ7w7z7YNYFsGO525WpZqbh3ogNe4vJ\nK6440r99rdMkk9rV3cKUOp62Z9hByS7/BxzYZQN+7j1Qmu92ZaqZaLg3Iivbtref3TPNaZJZpk0y\nKjiIQP8fwd3L7Hg1K193xo+fpU01YUDDvRELsj10T0+kfat4bZJRwSmupR2v5o6FdtLuD/7bTveX\nu8TtylQT0nA/jorqGpZsKeDsnk4vmbVzoP0AbZJRwalNb7jpfbjieXtn6/MX2aGFS/Lcrkw1AQ33\n41ixrYjyqhpG9UiDolzbJKNn7SqYiUC/K+1QBmfdC6tnw/Qh8M1MneovxGi4H0dWTh6REcLIbql6\n45IKLbEt4AcPwx1f22E0/v1zmDkGti92uzLlJxrux5GV7WFQ52RaxEXbG5fa9YfUbm6XpZT/pJ8G\nN86Bq16G8gJ4YSy8N00HJAsBGu4N2F9Wxeqd+20XyP07YMdSPWtXoUnE/mzfvRRG/xTWvG3Hqln8\nnDbVBDEN9wYs2uTBGBjdI017yajwEJMIF/4W7lwMnYbBh/fbCUK2LnS7MnUSNNwbsCDHQ1JsFAM6\nJ9teMu36QevubpelVNNL6wE3vANXvwYVB+ClCfDOrVC8x+3K1AnQcG9AVraHkd1aE12yC3YsgTMu\nc7skpZqPCPS+GO5aAuf83F5z+ttQWDQdaqrcrk75wKdwF5FxIrJBRHJE5P561l8vIqtFZI2ILBKR\nAf4vtflszy9je0GZvStVm2RUOItJgPMftE01Xc6Ej38NM0bDlgVuV6Ya0Wi4i0gk8AwwHugDXCsi\nfepstgU41xjTD3gYmOnvQpvT4SF+e6Zpk4xSYH/+r3sLrnkDqsrg5R/C2zfbcWtUQPLlzH04kGOM\n2WyMqQRmA0dNHGqMWWSMKXReLgY6+bfM5pWVk0f7VnF0iy60TTJ61q6Ubao5fYJtqjn3flg/zzbV\nZD0F1ZVuV6fq8CXcOwK5Xq93OMsacgvw7/pWiMhtIrJMRJbl5QXmLc81tYaFOfmM7pGGrH/fLtT2\ndqWOiI6H834Fd30DXc+B//wWZozSppoA49cLqiJyHjbcf1nfemPMTGPMUGPM0PT0dH/u2m++27mf\n/eVVtklm3Rxoq00yStUrtStcN9s211RX2KaaOXfqsMIBwpdw3wl09nrdyVl2FBHpD8wCLjXGBO2n\ne6i9/ey2lZD7DZxxaSNfoVSYO22sveA6+qew+k2YPhRW/ROMcbuysOZLuC8FeopIVxGJAa4B5npv\nICIZwLvAjcaYjf4vs/lkZXvo3b4lqds+tAv6aJOMUo2KSbA3QN2+ANJ6wpw74OWLIS+o4yCoNRru\nxphq4G7gI2A98JYxZq2ITBORac5mDwGtgWdFZJWILGuyiptQeWUNy7cV2i6Qa+dA2772hg6llG/a\n9oGpH8LFT8Oe1bYt/vPHoOqg25WFnShfNjLGzAfm11k2w+v5j4Ef+7e05vfNlnwqa2o5v0M1LFls\n+/cqpU5MRAQMmQK9JsBHD8CXT9jxan74JHQ71+3qwobeoeplYY6HmKgIhpR+ZRdok4xSJy+pDVwx\nC254F0wNvHKJHXGy1ON2ZWFBw93LgmwPQ7ukEP39XG2SUcpfelxgL7ie/TN7Bj99KKx4VS+4NjEN\nd0decQXf7ylmbEYt5C7WG5eU8qfoeLjgNzBtAaSfDnPvhpcmQt4GtysLWRrujoVOF8iLzDd2gY7d\nrpT/tekNU+bDJX+DvWvhuVHw2aN6wbUJaLg7FmR7SEmIpv3Oj6DNGbY7l1LK/yIiYPBkuHsZ9L0c\nvvoDPHcmbPrc7cpCioY7YIxhYY6H8ZkguYv1rF2p5pCUDpfPtNP8Abw6yY4bXxKYQ5MEGw13YFNe\nCXsOHOTyuOWA0fZ2pZpT9/PsRN3n/ALWvmcvuC5/GWpr3a4sqGm4Y5tkAPoWfQ5t+thJg5VSzSc6\nDs7/Ndyx0P4Ovn+vnQFq33q3KwtaGu7YIQeGpBwkbpcO76uUq9J7wZQP4JLpkPe9nRjk099DVbnb\nlQWdsA/3qppaFm/OZ2rqasBoe7tSbouIgME32guu/a6CBX+GZ0dCzqduVxZUwj7cV+UWUVpZw1kV\nWU6TTC+3S1JKASSmwWUzYPJckEh47XJ4+xYo3ut2ZUEh7MN9QbaHtlJIime5NskoFYi6nQt3LHJm\nf5oLzwyDZS/qBddGhH24Z2XncUvr7xBtklEqcEXH2dmfpi20E+jMuw9eHAd717ldWcAK63A/cLCK\nb3fsZ3zEYkjvrU0ySgW69NNgyjyY9Bx4suHvZ8N/fqcXXOsR1uG+eFM+qbUFdDqwSs/alQoWIjDw\nOnvBtf/VkPUkPHeWzuFaR1iHe1aOh0tiltsmGW1vVyq4JLaGSc/C5H+BqbVzuM69B8qL3K4sIIR3\nuGd7uCp+mR2lrs3pbpejlDoZ3cbYO1zPugdWvgbPjID177tdlevCNtx3FpVT7NlJr4o1etauVLCL\nSYAfPAK3fgaJ6fDmDfZxYLfblbkmbMN9YbaHcZFLtJeMUqGkwyC47XO44Lew8WN7Fr/8pbCcGCRs\nw31BjodJMUsx6afbMaaVUqEhMhrO/qntG9+uH7z/E3j5Ysjf5HZlzSosw7221vB9djaDzDpEm2SU\nCk1pPeCm9+Hip2H3atujJutJqKlyu7JmEZbhvm73AUZULCJCm2SUCm0RETBkCtz1DfS40PaJ/8d5\nsGuV25U1ubAM94U5HiZGfEN16mnaJKNUOGjZHq55HX70KpTsg3+cDx//BirL3K6syYRluK/ekM3w\nyO+J6neZ26UopZpTn0vsWfyg62HRX21TzeYv3a6qSYRduB+sqiE99yMiqdUukEqFo/gUO0H3TU5f\n+FcugX/dBeWF7tblZ2EX7su2FvIDFlPasps2ySgVzrqeA3d+DaPug1VvwPThsHZOyHSbDLtwX7F+\nIyMi1hPT73I7RoVSKnxFx8NF/2NvfmrRDv7vJph9PRzY5XZlpyzswj1yw/tEiiG6/+Vul6KUChQd\nBsKtn8NFv4dNn9qbn5a9ENRjxodVuBeUVjKo+AsK4jPtrEtKKXVIZBSM+om9+an9AJj3X3YwMk+O\n25WdlLAK96VrNzBC1lPR6xJtklFK1a91d3ux9ZLpsPc726NmwZ+D7uansAr30lXvESmGNiOudrsU\npVQgE7GTdN+1BHqNg09/DzPHwM4VblfmM5/CXUTGicgGEckRkfvrWS8i8ldn/WoRGez/Uk+NMYaM\nPR+zJ7ozke3OcLscpVQwaNEOfvQKXP06lHpg1gXw0a+hstTtyhrVaLiLSCTwDDAe6ANcKyJ1G6zH\nAz2dx23Ac36u85Rtz93OoJrv2JcxXptklFInpvcP7c1Pg2+Cr6fDs2fCps/druq4fDlzHw7kGGM2\nG2MqgdnApXW2uRR4xViLgWQRae/nWk/J7sVvESmG1sN+5HYpSqlgFJ8MFz8FUz6AiCh4dRLM+6nb\nVTXIl3DvCOR6vd7hLDvRbRCR20RkmYgsy8vLO9FaT0lMizRWJJ1Lh9OGNOt+lVIhJnO07VFz9n9D\nale3q2lQVHPuzBgzE5gJMHTo0Ga9DWzw+Kkwfmpz7lIpFaqi4+CCh9yu4rh8OXPfCXT2et3JWXai\n2yillGomvoT7UqCniHQVkRjgGmBunW3mApOdXjMjgf3GmPCdvFAppVzWaLOMMaZaRO4GPgIigReM\nMWtFZJqzfgYwH5gA5ABlgLZ/KKWUi3xqczfGzMcGuPeyGV7PDXCXf0tTSil1ssLqDlWllAoXGu5K\nKRWCNNyVUioEabgrpVQIEuPSlFIikgdsa+bdpgGeZt5ncwrl49NjC16hfHxuHFsXY0x6Yxu5Fu5u\nEJFlxpihbtfRVEL5+PTYglcoH18gH5s2yyilVAjScFdKqRAUbuE+0+0CmlgoH58eW/AK5eML2GML\nqzZ3pZQKF+F25q6UUmFBw10ppUJQSIZ7KEzo3RAfju1655jWiMgiERngRp0nq7Hj89pumIhUi8iV\nzVnfqfDl2ERkjIisEpG1IvJlc9d4snz4uWwlIu+LyLfOsQXNyLEi8oKI7BOR7xpYH5h5YowJqQd2\nWOJNQDcgBvgW6FNnmwnAvwEBRgLfuF23H4/tLCDFeT4+WI7N1+Pz2u4z7EilV7pdtx8/u2RgHZDh\nvG7jdt1+PLYHgCec5+lAARDjdu0+Ht85wGDguwbWB2SehOKZe0hM6N2ARo/NGLPIGFPovFyMnRUr\nWPjy2QHcA7wD7GvO4k6RL8d2HfCuMWY7gDEmWI7Pl2MzQAsRESAJG+7VzVvmyTHGfIWttyEBmSeh\nGO5+m9A7AJ1o3bdgzyiCRaPHJyIdgcuA55qxLn/w5bM7DUgRkS9EZLmITG626k6NL8c2HegN7ALW\nAD8xxtQ2T3lNLiDzpFknyFbNR0TOw4b7aLdr8bOngF8aY2rtSWBIiQKGABcA8cDXIrLYGLPR3bL8\nYiywCjgf6A58IiILjDEH3C0rdIViuIfyhN4+1S0i/YFZwHhjTH4z1eYPvhzfUGC2E+xpwAQRqTbG\nzGmeEk+aL8e2A8g3xpQCpSLyFTAACPRw9+XYpgKPG9tInSMiW4DTgSXNU2KTCsg8CcVmmVCe0LvR\nYxORDOBd4MYgPONr9PiMMV2NMZnGmEzgbeDOIAh28O3n8l/AaBGJEpEEYASwvpnrPBm+HNt27F8k\niEhboBewuVmrbDoBmSchd+ZuQnhCbx+P7SGgNfCsc3ZbbQJ01Lq6fDy+oOTLsRlj1ovIh8BqoBaY\nZYypt/tdIPHxc3sYeElE1mB7lfzSGBMUwwCLyBvAGCBNRHYAvwWiIbDzRIcfUEqpEBSKzTJKKRX2\nNNyVUioEabgrpVQI0nBXSqkQpOGulFIhSMNdKaVCkIa7UkqFoP8HMytuQYkfca8AAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5bc931b390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title(\"local maxima in difference curve\");\n",
    "plt.plot(xsn,ysn);\n",
    "plt.plot(xd,yd);\n",
    "plt.hlines(ymx, plt.xlim()[0], plt.xlim()[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 5: Calculate thresholds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Sensitivity parameter S\n",
    "# smaller values detect knees quicker\n",
    "S = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def threshold(ymx_i):\n",
    "    return ymx_i - (S * np.diff(xsn).mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.42528736])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Tmx = threshold(ymx)\n",
    "Tmx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 6: knee finding algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If any difference value (xdj, ydj), where j > i, drops below the threshold y = T|mxi\n",
    "for (x|mxi, y|mxi) before the\n",
    "next local maximum in the difference curve is reached,\n",
    "Kneedle declares a knee at the x-value of the corresponding\n",
    "local maximum x = x|xi.  \n",
    "**If the difference values reach\n",
    "a local minimum and starts to increase before y = T|mxi\n",
    "is reached, we reset the threshold value to 0 and wait for\n",
    "another local maximum to be reached.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mxmx_iter = np.arange(xmx_idx[0], len(xsn))\n",
    "xmx_idx_iter = np.append(xmx_idx, len(xsn))\n",
    "\n",
    "knee_, norm_knee_, knee_x = 0.0, 0.0, None                                                                                                                                                                 \n",
    "for mxmx_i in range(len(xmx_idx_iter)):                                                                                                                                                                    \n",
    "\n",
    "    # stopping criteria for exhasuting array\n",
    "    if mxmx_i == len(xmx_idx_iter) - 1:                                                                                                                                                                    \n",
    "         break\n",
    "        \n",
    "    # indices between maxima/minima\n",
    "    idxs = ( mxmx_iter > xmx_idx_iter[mxmx_i] ) * ( mxmx_iter < xmx_idx_iter[mxmx_i+1] )\n",
    "    between_local_mx = mxmx_iter[np.where(idxs)]\n",
    "\n",
    "    for j in between_local_mx:\n",
    "        if j in xmn_idx:\n",
    "\n",
    "            # reached a minima, x indices are unique\n",
    "            # only need to check if j is a min\n",
    "            if yd[j + 1] > yd[j]:\n",
    "                Tmx[mxmx_i] = 0\n",
    "            elif yd[j + 1] <= yd[j]:\n",
    "                print('If this is a minima, how would you ever get here: {}'.format(j,knee))\n",
    "       \n",
    "        if yd[j] < Tmx[mxmx_i]:\n",
    "            if not knee_x:\n",
    "                knee_x = j\n",
    "            # declare a knee\n",
    "            knee = x[xmx_idx[mxmx_i]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HXJ/tGAiQhwbCEfTGAQAREFERQUFtEqftS\n17qgtv3qT/vtt5u23/ZrbW2tuFBq3aHuWosbq6jsCAgk7LuYCQkkkD2Z8/vjDBIikAmZyc3c+Twf\nj3mQmbnkfIblzeXczz1HjDEopZRylwinC1BKKRV4Gu5KKeVCGu5KKeVCGu5KKeVCGu5KKeVCGu5K\nKeVCGu5KKeVCGu5KKeVCGu5KKeVCUU4NnJaWZrKzs50aXimlQtLKlSv3G2PSGzvOsXDPzs5mxYoV\nTg2vlFIhSUR2+nOcTssopZQLabgrpZQLabgrpZQLabgrpZQLabgrpZQLabgrpZQLabgrpZQLabgr\npZQLabi73WOP2YdSKqxouLvd++/bh1IqrGi4K6WUC2m4K6WUC2m4K6WUC2m4K6WUCzUa7iLynIh4\nRGTdCd4XEXlCRLaIyFoRGRL4MtUpW7DAPpRSYcWfM/fngQkneX8i0Mv3uB14uvllKaWUao5GN+sw\nxnwqItknOWQS8KIxxgBLRKStiHQ0xuwLUI2qOY70uN9/v7N1KBVijDFU1ngpqaihpKKG0soaSspr\nvn1eUlFDl/YJXD60k9OlHlcgdmLKAnbXe77H99p3wl1Ebsee3dOlS5cADK0adaTHXcNdhSFjDGXV\ndTaM6wVz6ZGwrvjuo7SihpKKWkoraqiu8570+1/QP8PV4e43Y8x0YDpAbm6uacmxlVKhraK6jqKy\nKorLqikqq6b4cDXFZdUUl1fXC+WjP9qz7VrqvCeOGhFIjosmJd4+kuOj6JgS5/va91q9948eF01y\nXBRRka23JyUQ4b4X6FzveSffa0opdVxHzqiLD1cfE9hFh6spLquy4e17FPlCvKKm7rjfKypCjgnd\nlIQYuqQmkhIfdfT1uGOD+ciPbWKjiIiQFv70LSMQ4f4eMFVEZgHDgRKdb1cqvBhjKK2s9YXx8cK5\nwWtl1VTXHn/KIzYqgrSkWNonxtA+MYae6Un266QYUhNjaJ9o30v1vdYmNgoRdwZ0czQa7iIyExgD\npInIHuBXQDSAMeYZYDZwEbAFKAduClaxSqmWV1lTR0FpJd+UVPJNaaXv6yoKSivZV1JBQWkVnkOV\n1NQdf/ojMSaS9kk2lDOS4+jXMdkX0vaR6nvvyGsJMZEa1gHgT7fM1Y28b4C7A1aRCiztcVcnYIzh\nYHkN35T6Qrukkn0lvvD2hXlBaSUHymu+83MTYiLJTIkjMzmO4d3a0yE5jrSkI2Ede0x4x0VHOvDp\nVIteUFVKtYyaOi+eQ1X2bPuYM+5jv65qMDUiAqmJsWSmxNKpXTy52e3ITI4jIznu2zDPSInTqZAQ\noOHudtrn7kper+Hrkgp2FpWzo6iMHfvL2FFUzr6SCr4pqaKorArTYJYkJiqCTF9ID+rUlgtP94V2\nchyZKbFkpsTToU0s0a24A0T5T8Pd7bTPPWTVeQ1fH7QBvr2ojJ37y2yQF5Wzq7j8mAuSsVERdGmf\nQFa7eHJOSzl6pu07285MjqNtQrSebYcRDXelHHQkwLfvL2OnL7h3+EJ8d3HFMTfRxEZFkJ2aSI/0\nRM7v24HstES6piaQnZpIZnKca1v61KnRcFcqyGrrvOw9WMGOonJ2FpX5grzcF+Dlx3SZxEdH0jU1\ngV4d2jCufwbdUhPpmppIt7REOrSJ1QBXftNwVypASsprWPd1CZsLDn0b5DuKytldXE5tvbskE2Ii\n6ZqaSJ+MNlzQP5NuaQnHBLhOnahA0HBX6hQcCfKv9pbw1R77467i8m/fT4iJJDs1kX4d2zAxJ5Ps\nVDuF0i0tkXQNcNUCNNzdTvvcm+2YIN9bwrq9JewsOhrkndrFMyArhSvP7MyArBT6ZrbRAFeO03BX\nqp6SihrW7y1hbSNBfkWuDfIBWSm0S4xxsGKljk/D3e20z/2EjgT5V74w1yBXbqLh7nba5w4cG+RH\nzsh3aJArF9NwV65TUV3Hl7sPfHuhs2GQZ7WNZ2CnFH7gC/KcrBTaa5Arl9FwV67gKa1kbr6HORsK\n+GzL/m/XTMlqa8/INchVuNFwVyHJGEPevkPMyStgbl4Ba/aUAHZ65ephXRjdJ51BndpqkKuwpeGu\nQkZVbR1LthUzN6+AuXke9h6sQATO6NyWBy7sw7h+GfTOSNIWRKXQcHe/EO9zLy6rZn6+hzl5BXy6\nqZCy6jrioyMZ1SuN+87vxXl9O5DeJtbpMpVqdTTcVatijGFrYdm30y0rdx7AayAjOZZJg7MY168D\nI3uk6QYQSjVCw93tQqDPvbbOy/IdB5ibV8CcvIJvO1v6d0xm6thejOvXgZzTUnTRLKWaQMPd7Vpp\nn3tpZQ0LNxYyJ6+ABRsLKamoISYygrN6pHLLqG6M7ZdBVtt4p8tUKmRpuKsWs6uo3E635BewdFsx\ntV5D+8QYxvXLYHz/DozqlU5SrP6RVCoQ9G+SChqv17B6z0HmbLDTLZsKDgPQs0MSt57TnXH9OjC4\nSzsidbpFqYDTcFcBV15dy/Nf7OC5z3aw/3AVkRHCmdnt+J+L+zGuXwbZaYlOl6iU62m4q4CprKnj\n5SU7eXrBVorKqhndO53LhmQxpncHUhKinS5PqbCi4e52LdDnXlVbx6xlu5k2fwueQ1WM6pnGT8b3\nZmjXdkEfWyl1fBru6pTV1Hl5fcUenpy3ma9LKhmW3Z4nrh7MiO6pTpemVNjTcHe7IPS519Z5efvL\nvTwxbzO7iysY3KUtj04ZxNk9U/XWf6VaCQ13twtgn3ud1/D+2q/565zNbNtfRk5WMg//MIcxfdI1\n1JVqZfwKdxGZAPwViARmGGP+0OD9FOBloIvvez5mjPlngGtVDvF6DR+u/4bHP9nEZs9h+ma24dnr\nh3JB/wwNdaVaqUbDXUQigWnAeGAPsFxE3jPGbKh32N3ABmPM90QkHdgoIq8YY6qDUrVqEcYY5uR5\n+PMnm8jbV0qP9ESevGYwF+V01KUAlGrl/DlzHwZsMcZsAxCRWcAkoH64G6CN2NO4JKAYqA1wraqF\nGGNYuKmQxz/ZxJo9JWSnJvD4lYP4/qAsveFIqRDhT7hnAbvrPd8DDG9wzJPAe8DXQBvgSmOMt+E3\nEpHbgdsBunTpcir1qiD7Yst+/vTJJlbuPEBW23gevXwglw3JIioywunSlFJNEKgLqhcCq4GxQA/g\nExFZZIwprX+QMWY6MB0gNzfXBGhsdTJ+9rkv31HMnz7eyJJtxWQmx/HbS3O4IrczMVEa6kqFIn/C\nfS/Qud7zTr7X6rsJ+IMxxgBbRGQ70BdYFpAqVdCs3n2QP328kUWb95OWFMuvvtefq4d10fXSlQpx\n/oT7cqCXiHTDhvpVwDUNjtkFnA8sEpEMoA+wLZCFqlN0gj73dXtLePyTTczN99A+MYb/vqgv14/I\nJj5GQ10pN2g03I0xtSIyFfgI2wr5nDFmvYjc4Xv/GeAR4HkR+QoQ4EFjzP4g1q381aDPPf+bUh7/\nZBMfrS8gJT6aBy7sw40js3WpXaVcxq+/0caY2cDsBq89U+/rr4ELAluaCqSthYf5y5zNvL/2a5Ji\norjv/F7cck43kuN0QS+l3EhP11xuZ2xb/trpLN7580LioiO5c3QPbj+3O20TYpwuTSkVRBruLrZs\nezE3DrwRrwi3jOrGHaN7kJoU63RZSqkWoOHuUmv3HOTm55dzWvUhXs57jY6Pvu90SUqpFqTh7kL5\n35Ryw3PLaJcYzSs/vYLMlBucLkkp1cL0DhWX2VZ4mOtmLCMuKpJXbx1BZkqc0yUppRyg4e4iu4vL\nuXbGUowxvHzrcDq3T7B97kd63ZVSYUPD3SU8pZVc94+llFXV8tItw+nZIcm+8f77R3vdlVJhQ8Pd\nBYrLqrl2xlL2H6rihZuH0f+0ZKdLUko5TC+ohriSihqu/8dSdhWX88LNwxjcRTelVkrpmXtIK6uq\n5ebnl7Op4BDPXj9UN6ZWSn1Lz9xDVGVNHbe9uILVuw8y7ZrBjOnTwemSlFKtiIZ7CKqp83L3K6tY\nvK2IP18xiAk5HU98sJ/ruSul3EWnZUJMndfw43+tZm6+h99emsPkwZ2cLkkp1QppuIcQr9fw4Jtr\n+c/affz8on5cO7xr4z9J+9yVCksa7iHCGMNv/r2eN1bu4cfjenHbud39+4na565UWNJwDwHGGP7v\nw428sHgnt5/bnfvO7+V0SUqpVk7DPQRMm7+FZxZu5drhXfjZxL6IiNMlKaVaOQ33Vu65z7bz2Meb\nuGxwFo9MytFgV0r5RcO9FZu1bBcPv7+BiTmZPDplIBERGuxKKf9on3sr9e7qvfzs7a8Y0yedv141\nmKjIU/x3WPvclQpLeubeCn28/ht++toahndrzzPXDSUmSn+blFJNo6nRyizaXMjUV79kQFYKM248\nk7joyOZ9Q+1zVyosabi3Isu2F3Pbiyvo0SGJF24aRlJsAGbNtM9dqbCk4d5KfLuhddt4XrplGCkJ\n0U6XpJQKYRrurUD9Da1fvXUEaUmxTpeklApxGu4O0w2tlVLB4Fe4i8gEEdkoIltE5KETHDNGRFaL\nyHoRWRjYMt1pz4Fyrmu4obVSSgVAo1fsRCQSmAaMB/YAy0XkPWPMhnrHtAWeAiYYY3aJiO4c0QhP\naSXXzljK4apaZt4+4uiG1oGmfe5KhSV/ztyHAVuMMduMMdXALGBSg2OuAd4yxuwCMMZ4Alumu9Tf\n0Pr5m4dx+mkpTpeklHIZf8I9C9hd7/ke32v19QbaicgCEVkpIjcEqkC3qb+h9Ywbz2RIsDe01j53\npcJSoC6oRgFDgYuBC4FfiEjvhgeJyO0iskJEVhQWFgZo6NBRXn10Q+tnrh/KWT1aYENr7XNXKiz5\nE+57gc71nnfyvVbfHuAjY0yZMWY/8CkwqOE3MsZMN8bkGmNy09PTT7XmkHRkQ+svdx3giasGc55u\naK2UCiJ/wn050EtEuolIDHAV8F6DY94FRolIlIgkAMOBvMCWGrqObGj9+ZYiHvvBICYOOMmG1kop\nFQCNdssYY2pFZCrwERAJPGeMWS8id/jef8YYkyciHwJrAS8wwxizLpiFh5KH/73h2w2tLxuiG1or\npYLPr8VLjDGzgdkNXnumwfM/An8MXGnuUF3r5e0v93L5kE5cN8KPDa2VUioAdD33IFu+o5jDVbVM\nzMl0pgDtc1cqLOnyA0E2N89DTFQEI3u2QGeMUkr5aLgH2bz8Akb2SCUhxqH/JGmfu1JhScM9iLYV\nHmZHUTnn93Ww7VH73JUKSxruQTQv367CcJ6T4a6UCkt6QTWI5uZ56JPRhk7tHFztMa4KjIDXCxH6\nb7lS4ULDPUhKK2tYvqOY287t3vKDV5fD+rdh5fMwwrd45+8yoV1XaNcN2neD9t2Pft22C0TpBiFK\nuYmGe5As2rSfWq9hbEtOyRRsgJX/hDX/gqoSSO0FW06Duki4aTIc2A7FO2Dn51B9uN5PFEjpDO2z\njx/+sW1a7jMopQJCwz1I5uYX0DYhmsGd2wZ3oJqKo2fpu5dCZAz0nwRDb4KuI+Ee+e7PMQbKCqF4\nuy/wtx39Ov8/UL7/2OMT048GfTtf8B/5OjEN5DhjKKUcpeEeBHVew8KNhYzpnU5UZJDmuT15sOKf\nsHYWVJZAak+44Hcw6GpIbKSnXgSSOthHl+Hffb+y1Bf6vuA/8vXOL2Dta4A5emxMUr3QbxD+yVkQ\nERnQj62U8o+GexCs2XOQorLqwHfJ1FTA+nd8Z+lL7Fl6v+9D7k3Q9ezjn0Ef6XG//37/x4lLho6D\n7OM7NVTCwV3fPeP35MHGD8Bbc/TYqHjoNR4GTIFeF0B0fJM+rlLq1Gm4B8G8PA+REcLo3gFa1tiT\n75tLn1nvLP23MOiaxs/Sj/S4NyXcTyY6DtJ720dD3joo3Xs08Pethbz37COmDfS7BHKmQPfREBkd\nmHqUUsel4R4E8/I9DO3ajrYJMaf+TWoqYMO79ix912KIiIb+37dz6dmjWuc8d0Sk7bxp2wUYbV+b\n+Cjs+BTWvQkb/m3/gUpItdcFcqZAl7O0RVOpINBwD7B9JRVs2FfKQxP7nto38OTbQF8zEyoPQvse\nMP4ROOMae/Ey1ERGQY+x9nHxn2HLHPjqDVg9E1Y8Z+flT58MOZfDaYNb5z9aSoUgDfcAm59vtw9s\n0pIDNZW+s/R/Hj1L7/c9O5eefY57Ai8qFvpebB9Vh2HThzbolz4Li5+0F2Jzptg5+vQ+TlerVEjT\ncA+wefkFdG4fT88OSY0fXLjx6Fl6xQEbbuMftnPpSS7fhjA2yYb4gClQXgz579ugX/QYfPooZOTY\ns/mcy+3NV0qpJtFwD6DKmjo+27KfK3M7Iyc6266ptBcYVz5vbyaKiLYXGof6ztIDPf8cCuu5J7SH\nITfYx6EC27e/7g2Y+xv76DTM/iPQ/1Jok+F0tUqFBA33AFq8rYjKGi9j+x0ngAo3+c7SX7Vn6e26\nwbjfwBnXuv8svSnaZMCIO+zjwA5Y95a9GPvB/4MPH7L/AA6YYqet4ts5Xa1SrZaGewDNy/OQEBPJ\n8G7tj33jo5/bOeWIKOh7iW8u/dyW6RI5lT731qJdNpzzU/vw5NuQX/cGvHcPvP9T6DnOBn2fiRCT\n6HS1SrUqGu4BYoxhXr6Hs3umERdd767MigOwbLq92ejiP9m7QltSoPvcndKhL4z9OZz33/D1l76g\nfws2fQDRCTbgcy63ga+LoCml4R4omwoOs/dgBfeM7XnsG+vfhrpqOOe/Wj7Y3UgEsobYx/hHbHfR\nujfsnbvr3oS4FDtlM/BKd3UaKdVEGu4BMje/ADjOxhxrZkF6v+Pfyq+aJyICss+2j4mPwrYFtuNm\n/bvw5cu24+asu217ZVQzbihTKgTprYEBMi/PQ05WMhnJcUdfLNpqV2ocdKWeQQZbZLRdx+ayZ+GB\nzTBpml0O4Z074S8DYNGfbMulUmFCwz0ADpRVs2rXAcb2bdAls/Y1QGDAFY7UFbai42HwdXDXYrju\nTejQD+Y+DI+fDrMfsAueKeVyOi0TAAs3FeI1HLsxhzH25qTuoyEly7niQqHPPVhE7AXWnuPgm3Ww\neJpdJnn5DHuX7Fn3HH/JY6VcQM/cA2Buvoe0pBgGZqUcfXHXEji4EwZe5Vxh6qjMHJj8NPz4Kzj7\nx7B9ETx3AcwYZy/GeuucrlCpgArJM/cxz4/5zmuX9L6E+0fe3+LvGyPs3jyVPlmVRETIt+//dP8u\nxksEl616iorVzzpW35Uf7Aag7L47HRm/1b2f3JExe+YT36ELEw4nMeWbNWS9fiMH49vSdvTPYPB1\njJl5SeutX993xfsLfrjgO8cEml9n7iIyQUQ2isgWEXnoJMedKSK1IjIlcCW2blUVWXi9cXTNOLon\naYzXy9iygyxKSKHC4Z2IzlpTxFlrihytoTWqiIjk7eR0rs/qzy/Su1EWmwQfPgiP9+f24r2k1VY7\nXaJSzSLGmJMfIBIJbALGA3uA5cDVxpgNxznuE6ASeM4Y88bJvm9ubq5ZsWJFM0pvHX4/O4/nPt/O\nql+Mp02cbwOK9e/A6zfC9W/bpW6dNGaM/TGc5979tXsZfPE3u4iZRNqbokZOhcwBTlem1LdEZKUx\nJrex4/w5cx8GbDHGbDPGVAOzgEnHOe4e4E3A06RKQ9y8fA/Du6UeDXawve1tOkK30c4Vppqu8zC4\n8iW4ZxXk3gx5/4ZnRsGLk2DzHHuRXKkQ4U+4ZwG76z3f43vtWyKSBUwGng5caa3frqJyNnsOH3vj\nUtl+2PIJDPiBbg4dqtp3g4sehZ+uh3G/tkszv3I5PHUWrHoJaqucrlCpRgWqW+YvwIPGGO/JDhKR\n20VkhYisKCwsDNDQzpnnuyv1mI051r0J3loYdLVDVamAiW8Ho34C962Fyc/ahd/emwqP58DCP0KZ\nXstQrZc/3TJ7gc71nnfyvVZfLjDLt4Z5GnCRiNQaY96pf5AxZjowHeyc+6kW3VrM21hI9/REstPq\nrUi4ZiZkDoSM/s4VVp/OtTdfVAwMusquV7N9IXzxJMz/rb3r9YyrYcTdkNaz8e+jVAvy58x9OdBL\nRLqJSAxwFfBe/QOMMd2MMdnGmGzgDeCuhsHuNmVVtSzZWsTYPvXO2gs32hULB2lvuyuJQPcxcN0b\ncNcSu9zwly/Dk7kw82rY8bnOy6tWo9FwN8bUAlOBj4A84DVjzHoRuUNE7gh2ga3VZ1v2U13nZWy/\neuG+Zpavy6IVdYI+9tjRNd1V4HToB5OehJ+sh3MfsDetPX8R/P08301RJ52hVCroGm2FDJZQb4V8\n6M21/GftPlb9cjzRkRH2L/NfBti/9NedtAu0ZWkrZMuoLrdTcounQfFWuyLlmJ/ZZQ500TgVQIFs\nhVQNHNmY49ze6TbYAXZ+BqV7dEomXMUkwJm3wNTlcNkMqKmAf10Lz54LGz/Q6RrV4jTcT8H6r0vx\nHKo6dqGwNbMgNtmeqanwFREJA38Ady+DS5+BqlKYeZWdrtn8iYa8ajEa7qdgbp4HERjTx7exdXU5\nbHgX+n/fLjerVGSU7aSZugK+/ySUF8ErU+Af42HrPA15FXQa7qdgXn4BZ3RuS2qSb6/O/P9A9WHt\nbVffFRkNQ66HqSvhkr9A6T54aTL8cyJs/9Tp6pSLabg3UeGhKtbsKTn2xqU1MyGlM3QZ6VxhJ7Jg\ngV5MbQ2iYiD3Jrh3FVz0GBzYCS98D56/xLZQKhVgGu5NNH+jXTrn2yUHDn0D2+bbG1wi9JdTNSIq\nFobdBvd+afd93b/JtlC+8H3YtdTp6pSLaBo10bw8D5nJcfTvmGxf+Op1MN7W2yWjfe6tU3QcDP8R\n3LcGLvxf8Gywm4e8dBnsCd0WYdV6aLg3QXWtl0WbCxnbrwNypHd5zSzIGgppvZwt7kTef98+VOsU\nHQ9n3W1DfvzD9g7nGefDK1fYr5U6RRruTbBsezFl1XVHlxz45isoWKcXUlXzxSTC2ffBj9fC+b+E\n3Uth+hi7rMG+tU5Xp0KQhnsTzM0vIDYqgrN7ptkX1syCiGg4/TJnC1PuEdsGzvkvu9fref8DOz+H\nZ8+Bf10HBeudrk6FEA13Px25K3Vkj1TiYyKhrtbOt/e6ABJTnS5PuU1cMox+wC43PPoh2LYQnh4J\nr/8QPPlOV6dCgIa7n7btL2NnUfnRu1K3L4DDBa33Qqpyh/i2cN7P7Jz8Offbu1yfGgFv3gr7Nztd\nnWrFNNz9NC+vQQvkmlkQ1xZ6X+hgVX7QPnd3SGgP5//CnsmffZ+9cW7aMHjrR1C01enqVCuk4e6n\nefke+ma2oVO7BKg6BHnvQ85ltm9ZqZaSmArjf2NDfsRddtmLJ8+Ed+6G4u1OV6daEQ13P5RU1LB8\nR/HRs/YN70FtRWh0yWifuzslpcOFv7PTNcN/ZK//PJkL7/8UDofVHvXqBDTc/bBocyG1XnN0yYE1\nM6F9d+h0prOF+UP73N2tTQZM+D3ctxqG3AirXoC/ngHzf2//h6nCloa7H+ble2ibEM3gLu3g4G7Y\n8RkMvEo3YVCtR/JpcMmf7VLDvcbBwj/AE4Nh2d+hrsbp6pQDNNwbUec1LNhYyJje6URGCHz1GmBg\n4BVOl6bUd6X2gCtehFvnQlpvmH0/TBtut/7TZYbDioZ7I1bvPkhxWTVj+2XYvxxrZkGXs6B9N6dL\nU+rEOuXCD/8DV/8LImPg9RthxjhdgTKMaLg3Yn6+h8gIYXSvdLvWx/5N2tuuQoMI9JkAd35uNwwp\n/dquQPnqleDJc7o6FWQa7o2Ym+9haNd2pCRE27P2yFjof6nTZflP+9xVRKTdMOSelXbdmp1f2Ltd\n351qA1+5kob7SewrqSBvX6ntkqmrgXVvQJ+J9q5BpUJNTIJdt+be1TD8Dnuy8sQQmPMbqCxxujoV\nYBruJzEv3/YLn9+vA2yZY/fBDIXe9vq0z101lJhq2yfvWQH9LoHP/mzbJxc/BbVVTlenAkTD/STm\n5Xno3D6eHulJtrc9IQ16nu90WU2jfe7qRNplw+Uz4PaF0HEgfPQzeyPU2tfB63W6OtVMGu4nUFlT\nx+db93N+3wyk8iBs/BAGTLEbHivlJqedATe8C9e9BbEp8Nat8PcxsHW+05WpZtBwP4HFW4uorPHa\nVSDXvwN1VXafVKXcquf58KNPYfJ0KD8AL10KL03WzUJClIb7CczNLyAhJpLh3dvbC09pfeC0wU6X\npVRwRUTAoCth6nK44LewdxU8e65dffLgLqerU03gV7iLyAQR2SgiW0TkoeO8f62IrBWRr0TkCxEZ\nFPhSW44xhvn5hYzqmUZs6U7YvcT2tutyAypcRMfByHvsmjVn3wvr34a/DYWPfg7lxU5Xp/zQaLiL\nSCQwDZgI9AeuFpH+DQ7bDow2xgwAHgGmB7rQlrSx4BB7D1bYKZm1rwESussNaJ+7ao74dnbj7ntX\nwYAfwOJp8MQZ8NlfoKbC6erUSfhz5j4M2GKM2WaMqQZmAZPqH2CM+cIYc8D3dAnQKbBltqy5Rzbm\n6JNup2S6nQMpIf2RlGqelE5w6VP2btfOw2HOr+yZ/JevgLfO6erUcfgT7lnA7nrP9/heO5FbgA+a\nU5TT5ud7GJCVQkbJWjiwPfR62+vTPncVSBmnw7Wvw43/hqQMePcueOYc2PSxLkzWygT0gqqInIcN\n9wdP8P7tIrJCRFYUFhYGcuiAKS6rZtWuA3ZjjjUzISoe+n3P6bJOnfa5q2Dodi7cNg+m/BNqyuHV\nH9jOmoINTlemfPwJ971A53rPO/leO4aIDARmAJOMMUXH+0bGmOnGmFxjTG56evqp1Bt0Czd58BoY\n1ysF1r9lgz22jdNlKdX6iNitJu9eBhP+YBfWe+Zs+PeP4XDrPHkLJ/6E+3Kgl4h0E5EY4CrgvfoH\niEgX4C3gemPMpsCX2XLm5ReSlhRLTtliu97GIO1tV+qkomJgxJ1w75dw5m2w6kX42xB70VWXM3BM\no+FujKkTTutmAAAMcUlEQVQFpgIfAXnAa8aY9SJyh4jc4Tvsl0Aq8JSIrBaRFUGrOIhq6rws3Ojh\nvD7pRKz9FyRlQrcxTpelVGhIaA8XPQp3LbF7Hsz5FUwbZjfx1vn4Fhflz0HGmNnA7AavPVPv61uB\nWwNbWstbufMApZW1TOgeDf/52J6NRPr1S6SUOiK9N1z7GmydZ/viX7sBup5tN/TWGwFbjN6hWs/8\nfA/RkcKoyoXgrbX7pIY67XNXTukxFn60CC7+MxRuhOnnwdt3Quk+pysLCxru9czN9zC8WyqxG16H\njAGQmeN0SUqFtsgoOPMWexPUyHvsngh/GwIL/g+qy52uztU03H12FZWzxXOYyZ3LYO9K91xI1T53\n1RrEpcAFj8DdS6HnOFjwv77lhV/T5YWDRMPdZ15+AQDjauaDRNhbrd1A+9xVa9K+O1z5EvxwNiSm\nwVu3wT/Gwe5lTlfmOhruPnPzPfRIiydl89t2rrBNptMlKeVe2WfDbQtg0lNQshf+MR5ev0lXngwg\nDXegrKqWpduK+WHW11Cy2x0XUpVq7SIiYPC1duPuc/8fbJwNf8uFuQ9D1SGnqwt5Gu7AZ1v2U13n\n5cLaBRCTBH0vdrokpcJHbBKM/bkN+f6TYNGf7KJkq17URcmaQcMdu1dqelwd6bs/hP6X2l3ilVIt\nK6UTXP53uHUutO0K790D00fD9k+driwkhX24e72G+Rs93Jm5Eak+5J4umSO0z12Fmk65cMvHcPk/\noOIgvPA9mHUtFG11urKQEvbhvv7rUjyHqrjIuxCSO0HXUU6XpJQSsRvST10OY38B2xbAtOH2jteK\ng05XFxLCPtzn5heQLgfJKPzc7rYU4bJfEu1zV6EsOh7Ovd/Oxw+60rcT1GBY9neoq3W6ulbNZUnW\ndPPzPdyVugoxXrtPqtton7tygzaZMGka/Gih3TBk9v12eeHNc5yurNUK63D3HKpkzZ4SLjGfwmlD\nIL2P0yUppU6m4yC7C9SVr9jlhF+5HF6eAvs3O11ZqxPW4b4gv5A+sov0sk3uPGtXyo1EoN8ldpOQ\nC34Lu5fCUyPsfHxlidPVtRphHe7z8j1cn7AEExEFOZc7XY5SqimiYuxiZPessvscL57m27T7ZV2v\nhjAO96raOj7fXMD3ZRHSc7xd50IpFXqS0mHSk3ZP13bd4N27Ycb5sHu505U5KmzDfdn2Ys6oXUNy\nbZG7p2S0z12Fi6whcPNHMPlZKP3aLkj29h1w6BunK3NE2Ib7vHwPU6I+w8SlQO8JTpejlAqEiAh7\nsnbPChj1E1j3pp2qCcP9XMMy3I0xfLFhJxMiVyCnT4boOKdLCh7tc1fhKLYNjPu13c81+xy7n+tT\nI2DTR05X1mLCMty3FpaRU7KQWFNpL8S4mfa5q3CW2gOumQXXvgkSCa9eETatk2EZ7vPzPUyOXERt\nSlfoPNzpcpRSwdZrHNz5hW2d3LUEnjoLPv4fqCx1urKgCctw/3L9OkZGbiDqjKttz6xSyv2OtE7e\nu8ouZfDF33ytk6+4snUy7MK9pKKGrntnE4Fx3wqQSqnGJXWwSxncNg/adYV377KdNXtWOl1ZQIVd\nuC/a5GFyxKcc7jDU7ueolApPWUPh5o/h0megZA/MGAvv3AWHCpyuLCDCLtw3ffkZvSP2knDmdU6X\n0jK0z12pE4uIgDOutqtOnn0frH3NTtV8/gTUVjtdXbOEVbjXeQ2ZO9+lVqKJyJnsdDlKqdYitg2M\nfxjuXgpdR8Inv4Cnz4JNHztd2SkLq3BfvbOQC7yL8HQ8D+LbOV1Oy9A+d6X8l9oDrn0NrnndPn/1\nB/DKFSG5C5Rf4S4iE0Rko4hsEZGHjvO+iMgTvvfXisiQwJfafDuWvkealJIy4nqnS2k52ueuVNP1\nvgDuXAzjH4GdX9hdoD75JVQdcroyvzUa7iISCUwDJgL9gatFpH+DwyYCvXyP24GnA1xnQKRtfYfS\niBQS++tyA0qpRkTFwNn32vn4gVfA53+18/GrXw2J1kl/ztyHAVuMMduMMdXALGBSg2MmAS8aawnQ\nVkQ6BrjWZtn3zTeMqF7CrtMm2N80pZTyR5sMuPQpuHUupHSCd+6Ef4yHva27ddKfcM8Cdtd7vsf3\nWlOPcdTORa8SKzWkDA+jKRmlVOB0yoVb5sClT8PBXfD3sTD7AaerOqGolhxMRG7HTtvQpUuXlhya\nqIQUvkwcxRmnn92i4yqlXCQiAs64BvpeAp/+ERJSna7ohMQYc/IDRM4Cfm2MudD3/GcAxpjf1zvm\nWWCBMWam7/lGYIwxZt+Jvm9ubq5ZsWJF8z+BUkqFERFZaYzJbew4f6ZllgO9RKSbiMQAVwHvNTjm\nPeAGX9fMCKDkZMGulFIquBqdljHG1IrIVOAjIBJ4zhizXkTu8L3/DDAbuAjYApQDNwWvZNUkR3rc\n77/f2TqUUi2q0WmZYNFpmRYyZoz9UZcgUMoVAjkto5RSKsRouCullAtpuCullAtpuCullAu16E1M\nygF6IVWpsKRn7kop5UIa7kop5UIa7kop5UIa7kop5UIa7kop5UIa7kop5UIa7kop5UIa7kop5UIa\n7kop5UKOLfkrIoXAzhYeNg3Yr+OF3FhuH8/Nn62lx3PzZzuiqzEmvbGDHAt3J4jICn/WQdbxWtdY\nbh/PzZ+tpcdz82drKp2WUUopF9JwV0opFwq3cJ+u44XkWG4fz82fraXHc/Nna5KwmnNXSqlwEW5n\n7kopFRZcGe4iMkFENorIFhF56Djvi4g84Xt/rYgMCfJ4fUVksYhUicj9QR7rWt9n+kpEvhCRQUEe\nb5JvvNUiskJERgVzvHrHnSkitSIyJVhjicgYESnxfbbVIvLLUx3Ln/HqjblaRNaLyMJgjiciD9T7\nbOtEpE5E2gdprBQR+beIrPF9tptOZZwmjNdORN72/dlcJiI5zRjrORHxiMi6E7wf0DwJGGOMqx5A\nJLAV6A7EAGuA/g2OuQj4ABBgBLA0yON1AM4EfgfcH+SxRgLtfF9PbIHPlsTR6b2BQH4wx6t33Dxg\nNjAliJ9tDPB+C/65bAtsALoc+XMT7F/Lesd/D5gXxM/238D/+b5OB4qBmCCO90fgV76v+wJzm/Fr\neS4wBFh3gvcDlieBfLjxzH0YsMUYs80YUw3MAiY1OGYS8KKxlgBtRaRjsMYzxniMMcuBmlMcoylj\nfWGMOeB7ugToFOTxDhvfn3AgEWjORRx/fu8A7gHeBDwtMFag+DPeNcBbxphdYP/cBHm8+q4GZgZx\nLAO0ERHBnhAUA7VBHK8/9gQAY0w+kC0iGacymDHmU1+9JxLIPAkYN4Z7FrC73vM9vteaekwgxwuU\npo51C/aMIqjjichkEckH/gPcHMzxRCQLmAw83Yxx/BrLZ6Tvv9ofiMjpQR6vN9BORBaIyEoRuSHI\n4wEgIgnABOw/mMEa60mgH/A18BVwnzHGG8Tx1gCXAYjIMKArzTvRaW49Lc6N4a4AETkPG+4PBnss\nY8zbxpi+wKXAI0Ee7i/Ag80IhqZYhZ0iGQj8DXgnyONFAUOBi4ELgV+ISO8gjwl2SuZzY8zJzk6b\n60JgNXAacAbwpIgkB3G8P2DPoFdj/6f3JVAXxPFanSinCwiCvUDnes87+V5r6jGBHC9Q/BpLRAYC\nM4CJxpiiYI93hDHmUxHpLiJpxphTWW/Dn/FygVn2f/ekAReJSK0xpqnB2+hYxpjSel/PFpGngvzZ\n9gBFxpgyoExEPgUGAZuCNN4RV3HqUzL+jnUT8AffFN4WEdmOnQtfFozxfL93N4G94AlsB7adwlgB\nqccRTk/6B/qB/QdrG9CNoxdbTm9wzMUcewFkWTDHq3fsr2neBVV/PlsXYAswsoV+LXty9ILqEOwf\nagn2r6Xv+Oc59Quq/ny2zHqfbRiwK5ifDTttMdd3bAKwDsgJ5q8lkIKdT04M8p+Tp4Ff+77O8P05\nSQvieG3xXbAFbsPOiTfn70I2J76gGrA8CeTD8QKC8qHs1etN2CvqP/e9dgdwh+9rAab53v8KyA3y\neJnYs7JS4KDv6+QgjTUDOID9L/BqYEWQP9uDwHrfWIuBUcEcr8Gxz3OK4e7nZ5vq+2xrsBenm/UP\npj+fDXgA2zGzDvhxC4z3Q2BWc8bx89fyNOBj39+3dcB1QR7vLN/7G4G38HWQneJYM4F92IaIPdjp\nzqDlSaAeeoeqUkq5kF5QVUopF9JwV0opF9JwV0opF9JwV0opF9JwV0opF9JwV0opF9JwV0opF9Jw\nV0opF/r/uvS9Z3wI/0wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5bc8f9d160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xticks(np.arange(0,1.1,0.1))\n",
    "plt.plot(xsn,ysn);\n",
    "plt.plot(xd,yd);\n",
    "plt.hlines(Tmx[0], plt.xlim()[0], plt.xlim()[1], colors='g', linestyles='dashed');\n",
    "plt.vlines(xmx, plt.ylim()[0], plt.ylim()[1], colors='r', linestyles='dashed');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vertical, red dashed line represents the x value of the knee point. The horizontal greeb dashed line represents the threshold value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.22222222222222221"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "knee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.55555555555555558"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# normalized x value where the knee was determined\n",
    "xsn[knee_x]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py36",
   "language": "python",
   "name": "py36"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}