03-FitCurve.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3ca35305",
   "metadata": {},
   "source": [
    "# Fit a curved line\n",
    "\n",
    "This notebook will help you plot the **strength-duration curve** for your earthworm.\n",
    "\n",
    "**Note**: If you're new to Jupyter Notebooks and/or Python, please be sure to work through the Introduction notebook before this one."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5200e250",
   "metadata": {},
   "source": [
    "## Step 1. Define values to plot\n",
    "With matplotlib imported, we can now use the scatter function by calling `plt.scatter()`. However, we need to define what to plot first. One straightforward way to think about this is to define an x variable and a y variable. Below, there are **lists** of values (defined in brackets `[ ]`) assigned to `x` and `y`. Replace these with your own values, depending on what you'd like to plot on the x and y axis of your strength duration curve.\n",
    "\n",
    "**Note**: Remember that in a scatterplot, each dot has both an x and a y value. Therefore, these lists should be the same length. The coordinate for each point will be the values at the same **index** in `x` and `y`. For example, the coordinate for the very first point will be `x[0],y[0]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ed7b96f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add your data points here\n",
    "x = [1,2,3,4,5,6]\n",
    "y = [1,0.8,0.5,0.45,0.4,0.4]\n",
    "\n",
    "%whos"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2915a9d8",
   "metadata": {},
   "source": [
    "## Step 2. Fit a hyperbolic curve to hypothetical data\n",
    "\n",
    "The strength-duration curve should follow a **hyperbolic** function:\n",
    "\n",
    "$$\n",
    "a = r + \\frac{r*c}{t}\n",
    "$$\n",
    "\n",
    "where **a** is the stimulus amplitude (in amps or volts), **r** is the rheobase, **c** is the chronaxie, and **t** is time (the width of the stimulus).\n",
    "\n",
    "We'll first define that function below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "22b1924c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the hyperbolic curve function\n",
    "def func(t, r, c):\n",
    "    return r+((r*c)/t) #hyperbolic function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d43e7ddb",
   "metadata": {},
   "source": [
    "To demonstrate what this function looks like, we'll use `np.linspace()` to generate a list of 100 numbers from 0 to 1, with 0.01 spacing. We'll also define a sample rheobase and chronaxie. Then, from those points, we'll fit a hyperbolic curve using `scipy.optimize.curve_fit()`, which uses non-linear least squares to fit a function to data. [[curve_fit documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0633d11c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvAAAAILCAYAAACO86igAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAABYlAAAWJQFJUiTwAABRA0lEQVR4nO3deZhcVZ3/8fc3a2chIawJOxIgiEBMEJVNFgeV3W0WF4QZRZxR0RmdcRyX4OiMOjo/cRdEUUEZN0QQRAQEggxLAFlGIiBhSUgCZN86nfT5/XFvdyqV6u7qqttdXdXv1/Pc59662znVt5L+1Olzz42UEpIkSZKaw4hGV0CSJElS9QzwkiRJUhMxwEuSJElNxAAvSZIkNREDvCRJktREDPCSJElSEzHAS5IkSU3EAC9JkiQ1EQO8JEmS1EQM8JIkSVITMcBLkiRJTcQAL0mSJDWRUY2uwFASEU8Ak4AFDa6KJEmSWts+wKqU0r79PdAAv7VJ48aN2+Gggw7aodEVkSRJUuv64x//yPr162s61gC/tQUHHXTQDvPmzWt0PSRJktTCZs+ezb333ruglmPtAy9JkiQ1EQO8JEmS1EQM8JIkSVITMcBLkiRJTcQAL0mSJDURA7wkSZLURAzwkiRJUhMxwEuSJElNxAAvSZIkNREDvCRJktREDPCSJElSEzHAS5IkSU3EAC9JkiQ1EQO8JEmS1EQM8JIkSVITMcAPEZs3b2bt2rWklBpdFUmSJA1hBvghYIcddmDUqFFMnDiR1atXN7o6kiRJGsJGNboCgpEjR3Yvb9iwgUmTJjWwNsWbM6fysiRJkvrPAD8EtLW1dS9v2LChgTUZGBdcsGXZAC9JklQfu9AMAa0e4CVJklQcA/wQYICXJElStQzwQ4ABXpIkSdUywA8BBnhJkiRVywA/BBjgJUmSVC0D/BBggJckSVK1DPBDgAFekiRJ1TLADwEGeEmSJFXLBzkNAa0e4KdNa3QNJEmSWocBfgho9QC/aFGjayBJktQ67EIzBLR6gJckSVJxDPBDQGmAb29vb2BNJEmSNNQZ4IcAW+AlSZJULfvADwGtHuCvvnrL8mmnNa4ekiRJrcAAPwS0eoA//fQtyyk1rh6SJEmtwC40Q0CrB3hJkiQVxwA/BBjgJUmSVC0D/BBggJckSVK1DPBDwNixY7uXDfCSJEnqjQF+CLAFXpIkSdUywA8BBnhJkiRVywA/BBjgJUmSVC0D/BBggJckSVK1DPBDgAFekiRJ1fJJrENAqwf4WbMaXQNJkqTWYYAfAlo9wM+b1+gaSJIktQ670AwBrR7gJUmSVBwD/BBQHuBTSg2sjSRJkoYyA/wQMGrUKEaOHAlAZ2cnmzZtanCNJEmSNFTZB36IaGtrY+3atUDWCj969OgG16g4F120ZfnccxtXD0mSpFZggB8iygP8dttt1+AaFefd796ybICXJEmqj11ohghvZJUkSVI1mi7AR8SCiEg9TIsbXb9aGeAlSZJUjWbtQrMS+FKF9WsGuR6FMcBLkiSpGs0a4FeklOY0uhJFMsBLkiSpGk3XhaZVGeAlSZJUjWZtgR8bEW8D9gLWAg8At6aUNje2WrUzwEuSJKkazRrgpwI/KFv3RESck1K6pa+DI2JeD5tm1F2zGhngJUmSVI1m7ELzXeBEshA/ATgE+BawD3BdRBzWuKrVzgAvSZKkajRdC3xK6YKyVQ8B50XEGuCfgDnA6/s4x+xK6/OW+VkFVLPfSgN8e3t7I6ogSZKkJtB0Ab4X3yQL8Mc2uiK1aOUW+FNPbXQNJEmSWkcrBfjn8vmEhtaiRq0c4K++utE1kCRJah3N2Ae+J6/I539uaC1q1MoBXpIkScVpqgAfEQdFxDYt7BGxD/DV/OVlg1qpghjgJUmSVI1m60LzV8A/RcStwJPAamA/4BSgDbgW+ELjqlc7A7wkSZKq0WwB/mbgQOClwFFk/d1XAHPJxoX/QUopNax2dWjlAD9nTuVlSZIk9V9TBfj8IU19PqipGbVygL+gZOBPA7wkSVJ9mqoPfCsbO3Zs93KrBXhJkiQVxwA/RLRyC7wkSZKKY4AfIgzwkiRJqoYBfogwwEuSJKkaBvghwgAvSZKkahjghwgDvCRJkqphgB8iDPCSJEmqhgF+iDDAS5IkqRoG+CHCAC9JkqRqNNWTWFtZKwf4d72r0TWQJElqHQb4IaKVA/xFFzW6BpIkSa3DLjRDRCsHeEmSJBXHAD9EGOAlSZJUDQP8EDFq1ChGjMgux+bNm9m0aVODayRJkqShyD7wQ0RE0NbWxrp164CsFX7ixIkNrlUxzj13y7L94SVJkupjgB9CWjXAX3zxlmUDvCRJUn3sQjOE2A9ekiRJfTHADyEGeEmSJPXFAD+EGOAlSZLUl8L6wEfEi4Fjgb2AnYD1wFLgfuDWlNLqospqVQZ4SZIk9aWuAB8RewDnAn8LTOtaXbZbAjZHxG+BbwDXpJRSPeW2KgO8JEmS+lJTgI+IHYA5wLuB0cAC4IfA3cBiYBkwDtgRmAG8EjgOeA0wPyL+KaV0XX1Vbz2lAb69vb2BNZEkSdJQVWsL/GPAWODbwPdSSnf1dUBETAL+mqzF/pqI+GBK6cs1lt+SbIGXJElSX2oN8D8A/iOltKTaA1JKq4CLgIsi4kygrfcjhh8DvCRJkvpSU4BPKZ1fT6EppV/Uc3yrMsBLkiSpLz6JdQhp1QD/yU82ugaSJEmto5AAHxEjgbEppXVl608AzgDWARellJ4oorxW1aoBfs6cRtdAkiSpdRT1IKcvAMsiYnLXioj4a+AG4H3AvwB3RcSeBZXXklo1wEuSJKk4RQX4Y4GbU0orS9Z9ElgBnAX8M7A98I8FldeSDPCSJEnqS1EBfk+yoSUBiIgXAQcCX0kpXZZS+gJwHfDagsprSQZ4SZIk9aWom1gnAatKXh9F9gTWX5esexg4vqDyWtLYsWO7l1spwJ922pblq69uXD0kSZJaQVEB/llg35LXrwbWA/NK1k0ENhVUXktq1Rb4a65pdA0kSZJaR1EB/n+B0yPiVGAD8CbgxpRSR8k++wILCyqvJbVqgJckSVJxiuoD/x/5ua4CrgfGAJ/p2hgRbcAxwJ0FldeSDPCSJEnqSyEt8CmlByPi5cA78lX/k1K6u2SXlwI3AT8qorxWZYCXJElSXwp7EmtK6UHgQz1suwN4fVFltSoDvCRJkvpSSBeaiLgpIs7qY5+3RcRNRZTXqgzwkiRJ6ktRfeCPA/bpY5+9gVcVVF5LMsBLkiSpL0UF+GqMw2Eke2WAlyRJUl8K6wNP9uCmbUREAHsBJwNPF1heyzHAS5IkqS81t8BHRGdEbI6IzfmqOV2vSyeyVvc/AzOBK+qvcusywEuSJKkv9bTA38qWVvdjgaeABRX22wy8ANwIfLuO8lpeqwb4b32r0TWQJElqHTUH+JTScV3LEdEJfDel9KkiKjVctWqAP/fcRtdAkiSpdRTVB35fYEVB5xq2WjXAS5IkqThFPYn1ySLOM9wZ4CVJktSXmgJ8RHyixvJSSunfazy25Y0ePZqIIKXEpk2b2LRpE6NGFTlQkCRJkppdrelwTo3HJcAA34OIoK2tjfXr1wPQ3t7eEgF+9uwty/PmNa4ekiRJraDWdHh8obVQt9IAv2HDBiZMmNDgGtXv3nsbXQNJkqTWUVOATyndUnRFlLEfvCRJknpT84OcNDBKA3x7e3sDayJJkqShqNAO1hGxF3AW8FJge2AlcC/wA0eqqY4t8JIkSepNYQE+It4FfBkYA0TJpjOBj0XE+Skln8nZBwO8JEmSelNIF5qIOBH4JtAOfAY4ATgon38a2AB8Ld9PvTDAS5IkqTdFtcB/GFgNzE4pPV6yfj7wu4j4HjAv3+/GgspsSQZ4SZIk9aaom1iPAH5cFt675et/ku+nXhjgJUmS1JuiAvw44Pk+9nku30+9MMBLkiSpN0UF+CfJ+rv35njgqYLKa1kGeEmSJPWmqD7wVwL/HBFfBz6aUlrRtSEiJgH/TtZ95vMFldeyWjHA//KXja6BJElS6ygqwP8ncDpwHvDWiPgD8CwwFTgMmAQ8ku+nXrRigD/ttEbXQJIkqXUU0oUmpbQKOBK4GBgJHA28GTiG7EvCxcBR+X7qxdixY7uXWyXAS5IkqTg1t8BHxH3At4DLU0qrU0orgXdHxHuBA4HJZE9inZ9S6iiktsNAK7bAS5IkqTj1tMAfBnwNWBQR346IVwCklDpSSg+llG7P54b3fjDAS5IkqTf1BPijgO/ny38L3B4R90fEP0TE5PqrNjy1YoDfbbctkyRJkupTc4BPKd2RUjoH2A34e+A+4FDgy2St8pdGxFHFVHP4aMUA/+yzWyZJkiTVp+6bWPP+799MKR0OvBT4JrAROAu4NSIeiojzI2JKvWVVEhFvi4iUT+8ciDIGUysGeEmSJBWnqAc5AZBS+kNK6R+AacDZwO+BFwP/DSyMiB8UWV5E7Al8FVhT5HkbyQAvSZKk3hQa4LuklDaklL6fUjoGOAiYC7QBbymqjIgI4LvAC2St/i3BAC9JkqTeFPUgp21ExA5k3Wj+jqwVHrJhJYvyfuAE4Lh83hIM8JIkSepN4S3wEXFCRPwIWAh8ETgYuJNspJrdCyrjIOCzwIUppVuLOOdQYYCXJElSbwppgY+IqcA5ZK3t+wIBLAcuAi5KKT1cRDl5WaOAHwBPAR+t8Rzzetg0o9Z6FcUAL0mSpN7U8yTWAE4B3gmcDIwkC+5zyYL7T1NKA5FAP0E22s3RKaX1A3D+hjLAS5IkqTf1tMA/RTYGfJDdSPoDstb2R4qoWCUR8XKyVvcvppTuqPU8KaXZPZx/HjCr1vMWwQAvSZKk3tQT4HcHbiFrbf9ZSmljMVWqLO86833gT8DHB7KsRjLAS5IkqTf1BPgDU0qPFlaTvk0EDsiXN2Q9eLZxcURcTHZz6wcGq2JFasUAf889ja6BJElS66g5wA9yeAdoBy7pYdsssn7xc4H5QM3daxqtFQP87IodliRJklSLmgJ8RHwV+PeU0pIaj3890JZS+lG1x+Q3rL6zh/PNIQvw30spfbuWOg0VrRjgJUmSVJxax4F/C/DniPhGfmNpnyJickS8OyLuBX4K7Fhj2S3NAC9JkqTe1NqFZjrwKeBc4NyIeBq4HbgHeJZsDPg2spA+A3gF8DJgLPBH4NSU0nX1Vb01jRkzpnu5o6ODzs5ORowo/HlbkiRJalI1BfiU0jLgvRHxOeA84Gzgb/Iple0ewGbgRuDrwDUppc5aK9xDfeYAc4o8Z6NEBG1tbd2t7+3t7YwbN67BtapP6f3GqfzTIUmSpH6p60msKaWngX8D/i0iDgaOBvYia3lfDywFHgBuSymtqrOuw0ZpgN+wYUPTB3hJkiQVp64AXyql9DDwcFHnG87sBy9JkqSe2Ll6CDLAS5IkqScG+CHIAC9JkqSeGOCHIAO8JEmSemKAH4IM8JIkSeqJAX4IMsBLkiSpJwb4IcgAL0mSpJ4Y4IcgA7wkSZJ6MiABPiKmRMSeA3Hu4cAAL0mSpJ4UFuAjYmJEfDEiFgPPA0+UbHt5RFwbEbOKKq+VjR07tnu5FQL8woVbJkmSJNWnkCexRsRkYC5wMHA/WYA/qGSXB4FjgL8B7i2izFbWai3wu+3W6BpIkiS1jqJa4P+NLLyfnVKaBfykdGNKaR1wC3BiQeW1tFYL8JIkSSpOUQH+DcD1KaXv97LPk8DuBZXX0gzwkiRJ6kkhXWiAPYCf9bHPGmByQeW1tFYL8IsWbVm2O40kSVJ9igrwq4Fd+thnX7K+8epDqwX43Uv+7pJS4+ohSZLUCorqQnM3cGpEbFdpY0RMA04mu9FVfWi1AC9JkqTiFBXgLwR2BK6NiNLRZ8hf/wRoA75cUHktzQAvSZKknhTShSaldH1EXAB8EngI6ACIiOeBKUAA/5JS+n0R5bU6A7wkSZJ6UtiDnFJKF5ANE/lLYDmwGUjAtcCrU0r/VVRZrc4AL0mSpJ4UdRMrACmlm4GbizzncGSAlyRJUk8Ka4FXcQzwkiRJ6okBfggywEuSJKknNXWhiYhOsv7t/ZVSSoV222lFBnhJkiT1pNYwfSvbBvgpwKFAJ/A0sBiYCuxJ1tL/ANnNreqDAV6SJEk9qSnAp5SOK32dP6jp98DPgQ+nlJ4o2bYv8AXgpcBra67pMGKAlyRJUk+K6gP/OWB5SulNpeEdIH/9JmBlvp/60GoBPqUtkyRJkupTVIB/DXB9TxtTSinfbgt8FUoDfHt7ewNrIkmSpKGmqAC/HTC5j30m5/upD63WAi9JkqTiFBXg/wj8VUTsWWljROwN/BXwfwWV19IM8JIkSepJUUM6/hfwQ+C+iPgy2Sg1S4BdgVcB7yNrgf+vgspraWPHju1ebm9vJ6VERDSwRvWZN2/L8uzZjauHJElSKygkwKeUrshHovks8MmyzQF0AB9KKf1PEeW1uohg7Nix3f3f29vbt2qVbzaHH75l2RtZJUmS6lPYQ5VSSv8vIn4OvI1syMjJZCPP3AtcnlJ6sqiyhoO2trbuAL9hw4amDvCSJEkqTqFPRc1D+meKPOdw1dbWxsqVKwH7wUuSJGmLom5iVcG8kVWSJEmVFNICHxHHVrtvSunWIspsdQZ4SZIkVVJUF5rfAdXenjiyoDJbmgFekiRJlRQV4D9F5QC/PfAy4EjgarIbWlUFA7wkSZIqKWoYyTm9bY+Is4GvAP9WRHnDgQFekiRJlQzKTawppUuB/wX+YzDKawWlD3MywEuSJKnLYI5Ccz9Q9c2uw50t8JIkSaqk0HHg+7DnIJfX1FopwE+b1ugaSJIktY4BD9QRMRI4B3gTMHegy2sVrRTgFy1qdA0kSZJaR1HjwP+5l/Pvms83Ah8torzhoJUCvCRJkopTVAv8CCoPI9kBPAjcBXwlpfTHgspreQZ4SZIkVVLUMJL7FHEebWGAlyRJUiXeVDpEtVKAv/rqLcunnda4ekiSJLWCIvvAfyml9OVe9vkH4J9SSi8qosxW10oB/vTTtyynSh2tJEmSVLWixoHfB9i+j322B/YuqLyW10oBXpIkScUZzAc5bUc2Eo2qYICXJElSJTV3oYmIvcpWbV9hHcBIYC/gjUBPw02qjAFekiRJldTTB34BWw8deX4+9SSAf6yjvGHFAC9JkqRK6gnw3ycL8AGcBTwA3F9hv83AC8CNKaXf1FHesGKAlyRJUiU1B/iU0tldyxFxFnBlSulTRVRKBnhJkiRVVtSDnAbzZthhwQAvSZKkSgzeQ1RpgG9vb29gTSRJkjSU1NQCHxHfIev//tGU0pL8dTVSSunvailzuLEFXpIkSZXU2oXmbLIA/zlgSf66GgkwwFehlQL8rFmNroEkSVLrqDXA75vPF5a9VkFaKcDPm9foGkiSJLWOmgJ8SunJ3l6rfq0U4CVJklQcb2IdogzwkiRJqqTWm1j3qrXAlNJTtR47nBjgJUmSVEmtfeAXkN2Q2l+pjjKHlbFjx3Yvb9iwgZQSEdHAGtXuoou2LJ97buPqIUmS1ApqDdPfp7YAryqNGDGCMWPGsHHjRgA2bty4VahvJu9+95ZlA7wkSVJ9ar2J9eyC61G1iPgccDhwALATsB54EvgF8NWU0guNqlvR2traugP8hg0bmjbAS5IkqTjNeBPrB4EJwA3AhcDlwCZgDvBAROzZuKoVy37wkiRJKld4f/Q8QL8UmAysBO5LKT1dYBGTUkrbpNmI+AzwUeBfgb8vsLyGMcBLkiSpXGEt8BGxf0TcQHaD65XApfl8QUTcEBEHFFFOpfCe+3E+37+IcoYCA7wkSZLKFdICHxHTgd8DOwKPA3OBxcBU4GjgRGBuRByZUnqsiDIrOC2fPzBA5x90BnhJkiSVK6oLzX+Shffzga+llDq7NkTECOB9wP8D/gP4yyIKjIgPARPJuuocTvZF4QHgs1UcO6+HTTOKqFtRyoeSlCRJkooK8CcC16aUvlK+IQ/zF0bEa4BXF1QewIeAXUte/xo4O6X0XIFlNNTEiRO7l1esWNG4ikiSJGnIKCrAjwHu72Of+4BjCiqPlNJUgIjYFTiSrOX9vog4NaV0bx/Hzq60Pm+Zn1VUHes1derU7uUlS5Y0sCaSJEkaKoq6ifUPwPQ+9pnOAPRPTyktSSldCZxE1o3n+0WX0SilAX7x4sUNrIkkSZKGiqJa4P8DuDIiXpdSuq58Y0ScArweOLOg8raRUnoyIv4PmBkRO6WUnh+osgZLqwT4U09tdA0kSZJaR1EBfkfgOuCaiLgRuBVYQtZH/VXACcDVwE4RcVbpgSmlIlvMd8vnmws8Z8O0Sheaq69udA0kSZJaR1EB/lIgAUF2o2qlm1VPZ8tQj+T7JvrR5SUfS35JSmll2foRwL8DuwC/Tykt70/lh6pWaYGXJElScYoK8OcUdJ6+nAz8Z0TMBZ4AXmBLK/+LyMaef9cg1WXAGeAlSZJUrpAAn1L6XhHnqcJvyW6GPRp4KbA9sBb4E/AD4MsppWWDVJcBZ4CXJElSuaJa4AdFSukh4L2Nrsdg2XHHHRkxYgSdnZ0sW7aM9vb2rR7u1CzmzKm8LEmSpP5rqgA/3IwcOZJddtmlu/V96dKl7Lnnng2uVf9dcMGWZQO8JElSfYoaB56IOCIifhYRj0dEe0RsrjBtKqq84cJuNJIkSSpVSAt8RLwJuILsC8EC4C7AsF4AA7wkSZJKFdWFZg7ZzaSnpJTmFnROYYCXJEnS1orqQjMd+JHhvXit8jAnSZIkFaOoAL8Y6CjoXCphC7wkSZJKFRXgfwL8RUSMKeh8yhngJUmSVKqoAP9JYAXw44jYu6BzCgO8JEmStlbUk1jXRcS5wM3AnyNiBbCy8q5pvyLKHC523XXX7mUDvCRJkgppgY+Io4HfA1OAzcA6ICpMhY07P1zYAi9JkqRSRQ0j+TlgNHAW8MOUUmdB5x32Jk+ezNixY2lvb2ft2rWsWbOGiRMnNrpa/fKudzW6BpIkSa2jqAB/GNkwkpcVdD7lIoKpU6fy5JNPAlkr/PTp0xtcq/656KJG10CSJKl1FNWlZQ2wrKBzqYxjwUuSJKlLUQH+WuBVBZ1LZewHL0mSpC5FBfiPAJMi4msRMaGgcypngJckSVKXovrAXwGsBs4DzoqIP9HzMJInFlTmsNHsAf7cc7cs2x9ekiSpPkUF+ONKlicAL+1hv1RQecNKs48Ff/HFW5YN8JIkSfUp6kFOju8+gJq9BV6SJEnFMXg3AQO8JEmSugxKgI+I7SLiPRExbzDKazUGeEmSJHUpqg98RRHxCuBdwF+S9Y23D3wNSvvAL1myhJQSEdHAGkmSJKlRCm+Bj4jJEfHeiPgDcDtwTr7pUuDoossbDsaPH8+kSZMA6OjoYPny5Q2ukSRJkhqlsBb4iDgKOBd4E9AGdDUR/xr4q5TS6qLKGo6mTp3KqlWrgKwbzQ477NDgGkmSJKkR6mqBj4gpEfGBiHgIuBV4O9n47/8NHJrv9ozhvX72g5ckSRLU0QIfEZcBbyBrbd8I/Jysm8yvU0qb830KqKLAAC9JkqRMPV1o3gJ0Ap8DPpdSWlFIjVRRsz/MSZIkScWoJ8CvAiYBHwJmRsSlwC9SSu1FVExba+YW+E9+stE1kCRJah31BPhpwF+T3bj6GuAkYGVEXAF8N6V0dwH1U66ZA/ycOY2ugSRJUuuo+SbWlNL6lNJ3U0qvBA4BvkY2zvt5wP9GxMP5azvCF6CZA7wkSZKKU8g48Cmlh1NK7wd2A94B/B44iCy8nx0RV0fEGRExsojyhqPSAL9kyZIG1kSSJEmNVOiDnFJK7SmlH6SUjgFeDFwIrABOIRul5qkiyxtObIGXJEkSDMCTWLuklB5JKX0Q2B14G3AbWb951WDnnXfuHpbzueeeY9OmTQ2uUfVOO23LJEmSpPoU9iTWnqSUNgI/BH4YEfsPdHmtavTo0ey0004899xzpJR47rnnmDatOb4PXXNNo2sgSZLUOgasBb6SlNKjg1leq7EbjSRJkgY1wKs+PsxJkiRJBvgmYgu8JEmSDPBNxAAvSZIkA3wTMcBLkiTJAN9EfJiTJEmSBjzAR8SUiJgw0OUMB7bAS5IkqZAAHxEnRsTnI2JKybpdIuIW4HlgWUT8dxFlDWcGeEmSJBXVAv8+4A0ppeUl674AHAM8DrwAnB8Rf1lQecOSAV6SJElFPYn1MOCWrhcRMQ54E3BDSuk1EbEd8CBwHvDjgsocdqZMmcLo0aPp6Ohg5cqVrF+/nnHjxjW6Wn361rcaXQNJkqTWUVSA3wVYVPL65UAbcClASml1RFwDvL6g8oalESNGsMsuu7Bw4UIgu5F1n332aWylqnDuuY2ugSRJUusoqgtNO1DaFHwMkIBbS9atAnYoqLxhy240kiRJw1tRAf4J4ISS128EHk0pLSxZtyfZDa2qgwFekiRpeCsqwH8POCQi7oyI24BDgB+W7XMoML+g8oYtA7wkSdLwVlQf+G8ArwD+CgjgauBzXRsj4iVkof4TBZU3bDXjw5xmz96yPG9e4+ohSZLUCgoJ8CmlDuAtEXFe9jKtLttlMfBSYEER5Q1nzdgCf++9ja6BJElS6yiqBR6AlNKqHtY/j/3fC9GMAV6SJEnFKaoPvAaJAV6SJGl4K6QFPiL+XOWuKaW0XxFlDle77rpr97IBXpIkafgpqgvNCLJx38ttD0zOlxcBHQWVN2yVt8CnlIiIBtZIkiRJg6mom1j36WlbREwHvgxMAF5TRHnD2cSJExk/fjzr1q1jw4YNrFq1ismTJ/d9oCRJklrCgPeBTyk9BrwB2B345ECX1+oiwn7wkiRJw9ig3MSaUtoA3AD8zWCU1+qacSx4SZIkFWMwR6HZBEztcy/1qTTAP/vssw2siSRJkgbboAT4iNgJeD3w9GCU1+r23nvv7uX58+c3sCaSJEkabEUNI/mJXs6/J3AG2Wg0/1pEecPdYYcd1r18//33N64iVfrlLxtdA0mSpNZR1DCSc/rYvgr4dErp8wWVN6w1W4A/7bRG10CSJKl1FBXgj+9hfSewHHgkpbSpoLKGvRe/+MWMHj2ajo4OnnjiCVauXOlQkpIkScNEUePA31LEeVSdMWPG8OIXv5g//OEPADzwwAMcc8wxDa6VJEmSBsNgjkKjAs2cObN7uRm60UiSJKkYNbXAR8RetRaYUnqq1mO1RTP1g99tty3LixY1rh6SJEmtoNYuNAuAVMNxqY4yiYgdyYajPAU4hOzprhuBB4HvAt9NKXXWev5mUtoC39WVZqhyqHpJkqTi1Bqmv09tAb5ebwa+ATwL3Aw8BewKvAH4NvC6iHhzSqkRdRtUpS3wDz30EB0dHYwePbqBNZIkSdJgqCnAp5TOLrge1foTcDrwq9KW9oj4KHAX8EayMP+zxlRv8Oywww7stddePPXUU7S3tzN//nxe8pKXNLpakiRJGmBNdRNrSummlNLV5d1kUkqLgW/mL48b9Io1SDP1g5ckSVIxmirA96Ejnw+b8eYdiUaSJGn4KWQc+Ij4TpW7ppTS3xVRZln5o4Cz8pe/rmL/eT1smlFYpQZBM93IKkmSpGIU9STWs/vYnoDI54UHeOCzwEuAa1NK1w/A+Yek8hb4lBIR0bgKSZIkacAVFeD37WH99sDLgI8Dvwc+UlB53SLi/cA/AY8Ab6/mmJTS7B7ONQ+YVVztBtY+++zDpEmTWLVqFc8//zyLFi1i9913b3S1JEmSNIAK6QOfUnqyh+kPKaVvA0cDrwVeXUR5XSLivcCFwP8Bx6eUlhV5/qFuxIgRHHrood2v7QcvSZLU+gblJtaU0tPA1cD5RZ0zIj4AfAV4iCy8Ly7q3M3EG1klSZKGl6K60FRjCbB/ESeKiH8h6/d+P/AXKaXnizhvM2qGG1nvuafRNZAkSWodgxLgI2IkcAKwsoBzfRz4FDAPOGm4dZsp1wwt8LMr3nEgSZKkWhQ1jOSxvZx/T+AcYCbw7TrLeQdZeN8M3Aa8v8KoKwtSSpfWU04zOfjggxk5ciSbN2/mscceY/Xq1Wy33XaNrpYkSZIGSFEt8L8jGyKyJwHcCny4znK6RrsZCXygh31uAS6ts5ym0dbWxowZM3j44YdJKfHggw9y5JFHNrpakiRJGiBFBfhPUTnAdwLLgbtSSnfVW0hKaQ4wp97ztJqZM2fy8MMPA1k3GgO8JElS6yokwOfBWg0yc+ZMLr/8cmBo3sha2ssp9fZ3GkmSJPVpUIaR1MBqhhtZJUmSVAwDfAs47LDDupcfeOABNm3a1MDaSJIkaSAVFuAj4oiI+FlEPB4R7RGxucJkshwAO++8M7vtthsAGzZs4NFHH21wjSRJkjRQihpG8k3AFWRfCBYAdwGG9UE0c+ZMFi1aBGTdaA466KAG10iSJEkDoagW+DnAWuDYlNKLUkrHpJSOrzQVVJ7KNMMTWSVJklS/ogL8dOBHKaW5BZ1P/VTaD94bWSVJklpXUQF+MdBR0LlUA0eikSRJGh6KCvA/Af4iIsYUdD7103777ceECRMAWLJkCYsXL25wjSRJkjQQigrwnwRWAD+OiL0LOqf6YeTIkRx66KHdr+0HL0mS1JoKCfAppXXAucAxwJ8j4oWI+HOF6fEiylNlpf3g582b18CaSJIkaaAUEuAj4mjg98AUYDOwDogKkw+OGkAvf/nLu5evu+66BtZkawsXbpkkSZJUn0LGgQc+B4wGzgJ+mFLqLOi86odTTjmFESNG0NnZye9//3uWLl3KLrvs0uhqkT9jSpIkSQUoqkX8MLJhJC8zvDfOzjvvzJFHHglAZ2cn11xzTYNrJEmSpKIVFeDXAMsKOpfqcOaZZ3YvX3XVVY2riCRJkgZEUQH+WuBVBZ1LdTjjjDO6l3/zm9+wdu3aBtYms2jRlkmSJEn1KSrAfwSYFBFfi4gJBZ1TNZg+fToHH3wwABs2bOA3v/lNg2sEu+++ZZIkSVJ9igrwVwCrgfOAxRExLyJuqjDdWFB56oXdaCRJklpXUaPQHFeyPAF4aQ/7pYLKUy/OOOMMPvOZzwBw9dVXs2nTJkaNKupSS5IkqZGKepDTiCqnkUWUp97Nnj2b3fP+KsuWLeP2229vcI0kSZJUFB+s1IJGjBix1c2sv/jFLxpXGUmSJBXKAN+iygN8SvZekiRJagU1dYyOiGPzxbtSShtKXvcppXRrLWWqf4477jgmTZrEqlWrWLBgAQ8++CCHHnpoo6slSZKkOtV6Z+PvyG5IPQj4U8nratgPfhCMGTOGk08+mSuuuALIWuEN8JIkSc2v1gD/KbLA/nzZaw0hZ555ZneAv+qqq/jEJz7R4BpJkiSpXjUF+JTSnN5ea2h43etex+jRo+no6ODee+/lqaeeYq+99mp0tSRJklSHQm5ijYi9ImJSH/tsFxGmx0E0adIkTjjhhO7Xv/zlLxtYG0mSJBWhqFFongDO72Of9+f7aRCVPpW1UcNJprRlkiRJUn2KCvCRTxpiTj/99O7l3/3udyxfvryBtZEkSVK9BnMc+KnA2kEsT8Buu+3GEUccAcDmzZu59tprG1wjSZIk1aPWUWiIiLPKVs2ssA6yYSP3At4GPFhreardmWeeyV133QXAZZddxlvf+tYG10iSJEm1ilqf0BkRnVQ3dGRX15p1wBtSSr+pqcBBEBHzZs2aNWvevHmNrkqhHn/8cfbff//up7H+4Q9/GNQx4Ut/nLNnD1qxkiRJQ9bs2bO59957700p9Tsd1dwCD5yTzwP4DvAL4KoK+20GXgDuSCmtqKM81Wi//fbjjW98Iz/96U8B+NznPsfll18+aOUffviWZW9klSRJqk/NLfBbnSTiZuC7KaXv11+lxmnVFniAe+65h5e97GUAjBw5kkcffZR99913UMqOktubDfCSJEn1tcAXchNrSun4Zg/vre7www/nxBNPBLKbWb/4xS82uEaSJEmqxYCNQhMRp0fElyLiwoh440CVo+p95CMf6V6+5JJLWLp0aQNrI0mSpFrUHOAj4rSIuDUiXlVh23eBK8ke3vQ+4McR8bPaq6kinHjiicyaNQuADRs28JWvfKXBNZIkSVJ/1dMCfzowC7izdGVEnAq8g2zUmU8D/wL8GTgzIv6mjvJUp4jYqhX+a1/7GqtXr25gjSRJktRf9QT4I4DbUkobytb/LdnwkueklD6RUvov4BhgA+AA5A32hje8genTpwOwfPlyLr744gbXSJIkSf1RT4CfCjxcYf2xwAqgu8tMSmkx8CvgpXWUpwKMHDmSD3/4w92v//u//5uNGzc2sEaSJEnqj3oC/BRgq+QXEXsBOwBz07bjUz4B7FhHeSrIWWedxa677grAwoULB3VMeEmSJNWnngC/GtijbF3XOJb39XBMeXcbNUBbWxsf/OAHu19//vOfp7Ozs4E1kiRJUrXqCfAPAqdExMSSda8n6/8+t8L++wLP1lGeCnTeeecxadIkAB555BF++ctfDlhZ06ZtmSRJklSfegL85WTdaG6JiPdHxFfJblJdDNxcumNEBHA08H91lKcCTZ48mfe85z3drz/xiU/Q0dExIGUtWrRlkiRJUn3qCfCXANeT3Zj6/4C/BzYB56eUNpfteyLZTa+/raM8FewDH/gA48aNA+DBBx/kc5/7XINrJEmSpL7UHOBTSp3AKcDbgW+Sjfn+8pTSTyvsvhNwITBw/TTUb1OnTuXTn/509+tPfepTPPxwpYGFJEmSNFTU0wJPSqkzpXR5Sukf8jHf7+9hvytSSh9MKS2spzwV7/zzz+flL385AB0dHfzd3/0dmzeX/wFFkiRJQ0VdAV7Nb+TIkVxyySWMHj0agDvvvJMLL7yw0DKuvnrLJEmSpPoY4MXBBx/Mxz/+8e7XH/vYx3jssccKO//pp2+ZJEmSVB8DvAD4yEc+wqGHHgrA+vXreec73+nY8JIkSUOQAV4AjB49mu985zuMHDkSgFtuuYWLLrqowbWSJElSOQO8us2ePZsPf/jD3a//+Z//maeffrqBNZIkSVI5A7y28slPfpIDDzwQgNWrV3POOeewadOmBtdKkiRJXQzw2kpbWxuXXHIJ2cNz4cYbb+S9730vKaUG10ySJElggFcFRx11FB/72Me6X3/rW9/i85//fANrJEmSpC4GeFV0wQUX8Na3vrX79Uc+8hGuuOKKBtZIkiRJYIBXDyKCSy65hOOOO6573Tve8Q5uu+22xlVKkiRJBnj1bOzYsfz85z/noIMOAmDjxo2cccYZzJ8/v8E1kyRJGr4M8OrVlClTuPbaa9l1110BWL58OSeffDJLly6t+hyzZm2ZJEmSVB8DvPq0zz77cM011zB+/HgA/vznP3PKKafw/PPPV3X8vHlbJkmSJNXHAK+qHH744VxxxRWMGJF9ZO655x6OOuoonnjiiQbXTJIkaXgxwKtqp512Gt/85je7x4j/05/+xJFHHsl9993X4JpJkiQNHwZ49cu73vUufvKTnzB27FgAFi9ezLHHHssNN9zQ4JpJkiQNDwZ49dsb3/hGbrjhBrbffnsA1qxZw8knn8xll11Wcf+LLtoySZIkqT5NF+Aj4k0R8ZWIuC0iVkVEiojKyVED5phjjmHu3LnsscceAGzatIm3v/3tfPaznyWltNW+7373lkmSJEn1aboAD3wMeC8wE1jY2KoMbwcffDB33HEHL3nJS7rX/eu//iunnnpqv4aZlCRJUvWaMcB/EDgAmAS8p8F1Gfb22GMPbrvtNl71qld1r7v22ms55JBDuO666xpYM0mSpNbUdAE+pXRzSunRVN5PQw2z/fbbc/311/OP//iP3euWLl3KySefzPnnn9/AmkmSJLWepgvwGprGjh3LF7/4Ra6//nqmTp3avf7LX/5yA2slSZLUeoZlgI+IeZUmYEaj69bsTjrpJB544AFOO+20itvXrFkzyDWSJElqLcMywGtg7bzzzlx11VV8/etfp62tbattM2bM4Ic//OE2I9VIkiSpOsMywKeUZleagEcaXbdWERG85z3vYd68eVutX7hwIW9961s59thjfYKrJElSDYZlgNfgefGLX1xx/dy5c5k9ezbvfve7eeaZZwa5VpIkSc3LAK9B9eEPf5jRo0cDkFLioosuYr/99uO8887jiSeeaHDtJEmShj4DvAbcqadumT7/+c/z0EMP8brXva57+8aNG/nWt77F/vvvzznnnMOjjz7awNpKkiQNbQZ4Dbirr94yARxwwAFce+21XH/99Rx55JHd+23evJlLL72UGTNm8Ja3vIW77767QTWWJEkaupouwEfEmRFxaURcCnwkX/3KrnUR8YUGVk/9cNJJJzF37lxuuukmjj/++O71nZ2d/OhHP+KII47g8MMP55JLLmHdunUNrKkkSdLQ0XQBHpgJvCOfXpOve1HJujc1plqqRURw/PHHc9NNNzF37lxe85rXbLV93rx5vPOd72T33XfnAx/4APPnz29QTSVJkoaGpgvwKaU5KaXoZdqn0XVUbY466ih+/etfc/fdd/OOd7yDsWPHdm9bsWIFF154ITNmzODII4/kq1/9KkuXLm1gbSVJkhqj6QK8ms+cOVumahx++OFceumlLFy4kC984QtMnz59q+133HEH73vf+9htt9143etex2WXXcbq1auLrrYkSdKQFD4Rc4uImDdr1qxZ5Q8fUn0itizX8nHr7Ozkt7/9Ld/4xje45ppr2LRp0zb7jBs3jte+9rWcccYZnHLKKey000511FiSJGlgzZ49m3vvvffe/GGi/WILvIa8ESNGcNJJJ3HllVeyaNEivv71r3PUUUdttc/69eu58sorOfvss9l111059thj+cIXvuCQlJIkqeXYAl/CFviBUW8LfE8WLFjAj370Iy6//HIefvjhHvebPn06r371q/mLv/gLjj/+eKZMmVJcJSRJkmpQTwu8Ab6EAX5gDFSAL/XII4/wy1/+kquuuoo77riDnj7XI0aM4PDDD+fVr341J5xwAq94xSuYMGHCwFRKkiSpBwb4ghjgB8ZgBPhSS5cu5Ve/+hVXXXUVN9xwQ69jyI8aNYpZs2Zx9NFHc8wxx3DUUUex8847D3wlJUnSsGaAL4gBfmAMdoAv1d7ezh133MENN9zADTfcwD333NNj63yXAw44gCOOOKJ7Ouyww2hraxukGkuSpOHAAF8QA/zAaGSAL7ds2TJuvvlmfvvb33Lbbbf12ne+y+jRoznssMM4/PDDeelLX8rMmTN5yUtewvjx4wehxpIkqRUZ4AtigB8YQynAl1u2bBm33347c+fO5bbbbuOee+6ho6Ojz+NGjBjBgQceyMyZMznssMM45JBDOPjgg9lrr72I0jcsSZJUgQG+IAb4gTGUA3y59evXc//993P33Xdz1113cdddd/VrKMqJEydy8MEHd08zZsxgxowZ7L333owcOXIAay5JkppJPQF+1EBUSGpW48aN45WvfCWvfOUru9ctX76cu+++m/vuu4/777+f+++/n/nz51fsS79mzRruvPNO7rzzzq3Wjx07lunTp3PggQdy4IEHsv/++zN9+nT2228/pk2bZqu9JEmqmgFeA+5d72p0DeozZcoUTjrpJE466aTudWvXruWhhx7ivvvu44EHHuDhhx/moYceYtmyZRXP0d7ezsMPP1yxz/348eN50YtexPTp09l33327p3322Yd9993XYS4lSdJW7EJTwi40qkdKiSVLlnQH9Ycffpj58+czf/58Fi9eXPN5d955Z/bee2/22muv7mnPPffsnu+yyy52z5EkqcnYhUYaAiKCqVOnMnXqVE488cSttq1cubI7zM+fP5/HH3+cxx9/nMcee4zly5f3et7nnnuO5557jnvuuafi9lGjRjFt2jR233139thjD3bffXd23313pk2bxrRp09htt92YNm0akydPtquOJEktwAAvDYLJkyd3jytfbtmyZd2B/oknnuCJJ55gwYIFPPHEEzz55JN9joqzadMmnn76aZ5++ule92tra2PatGlMnTqVXXfddZv5Lrvs0j1NmjTJsC9J0hBlgJcabIcddmCHHXbgZS972TbbNm/ezLPPPstTTz1VcXrmmWd44YUXqipnw4YN3V8Q+jJmzBh22WUXdt55Z3beeWd22mkndtppp62Wd9xxx60mH3YlSdLgMMBrwJ177pbliy5qXD2a0ciRI9ljjz3YY489OPLIIyvus379ehYtWsTChQt55plnWLhwIQsXLuTZZ5/tnhYtWsS6deuqLnfjxo0888wzPPPMM1UfM378+O4wv8MOOzBlypRt5lOmTGH77bffZm4ffkmSqmeA14C7+OItywb44o0bN4799tuP/fbbr8d9UkqsXr2aZ599liVLlrBkyRIWL1681fy5555j6dKlLF26lLVr1/a7HuvWrWPdunV9duWpZOLEiWy//fZMnjyZ7bffvnt58uTJTJo0aZvlSZMmbTONHTu23+VKktSMDPDSMBAR3UH3wAMP7HP/devWdQf6559/nueee47nn39+q+UXXnhhq2nTpk0112/NmjWsWbOmXy3+5UaPHs12223X6zRx4kQmTpy41XLXNGHChK2WJ0yYwKhR/hcpSRp6/O0kaRvjx49n7733Zu+9965q/64W/hdeeIFly5Z1T8uXL99qecWKFSxfvnyr5ZUrVxZS546Oju6yijJmzJjuMF86jR8/fpvl8ePH9zqNGzdum/m4ceMYPXp0YfWVJA0PBnhJdStt4d933337dWxnZyerV69mxYoVrFixgpUrV3Yvr1q1ipUrV3bPu6bVq1ezevVqVq1a1b2tnr8A9GTjxo1s3Lixz6E+6zFy5MjuMD9+/Hja2tq6X1dabmtr657KX5dPY8eO3Wa5dD569GhHG5KkJmSAl9RQI0aM6O7jXm2Lf7mUEu3t7d3Bvnxas2ZN97xrefXq1axdu5Y1a9Z0z7umtWvXsnbtWjo7Owt+t9vavHlzd7mN0BXmq5nGjBnT63LpvGsqf93XNHr06O5lb26WpMoM8JKaXkR0tzTvvPPOhZwzpcSGDRu6w3xXyO+6WXft2rXd87Vr17J+/frubaVT6fry5fXr17N58+ZC6lurDRs2sGHDhobWoScRsVWorzTvaV35ttJp1KhR26zrbXvX677mXcul60uXR44c6V88JBXCAC9JFUREd/eVnXbaacDK6ejo6A7zlaYNGzZsM++auvZpb2/fan3XVLq+fLm9vX1Auh0VqesvK+3t7Y2uSmG6Qn1P0+jRoxk5cmSf+3VN5fv29LqaeX+3lS9Xel3ttpEjRzJixAi/4EhVMsBLUgN1tdxOmjRp0MvevHlzd0AuDfbl08aNG3tdrjQvn7qO6ejooKOjY5ttXeu65u3t7aSUBv1nMtA2bdo05L84NdKIESN6Dfj92Va+rvR113IR60aMGFFxXW/79HRMX/v3dZ6epr62V7NPRPglawgxwEvSMDVy5MjuUXKGok2bNm0V+EuDf9f68m2l63ta19HRsdW5e5q6wnb5/uXrKs1L99+8eTMdHR2Dck9Fs+vs7KSzs5OOjo5GV0W96OvLQF9T15eBnl5Xu670dekXjFqO22uvvbjgggsa/aOtmgFekjQkdXUDGTduXKOrUojOzk42b97cHe67gn1X4C9d3/W6dLl8W+lUvr7r3Js3b+7eVum8pdtLX5cf09t5ept6OmfXz6J08gtO8+j6otVKDj30UAO8VOqTn2x0DSSp8bpa/EaPHt0yX0qKlFLa6ktOabAvD/p9fXHo7Ziu5b7Wdb3ubb9K80rnSilttb2n/ft6Xb6+vIxK28rXda0v/XlXWl9p/1bs1tZlxIgRja5CvxjgNeDmzGl0DSRJQ11EdPczHzNmTKOrox6UB/y+Qn9nZ+dW2yvt3/UFp5r1pV88ur5UlG8rX650rvJ9p0yZ0ugfbb8Y4CVJklSV0i9aapzm+ntBq5szByKqm849d9vjzz23+uMrNYufdlr1x1900bbHz55d/fFXX73t8bvtVv3x8+Zte3y1x0bAokVbH7toUf+OLzdvXvXH7rbbtsdffXX1x8+eve3xF11U/fGnnbbt8X72/Oz52fOz52fPz95w/ew1YVcBA7wkSZLURAzwkiRJUhOJVr6juL8iYt6sWbNmzav0pyJJkiSpILNnz+bee++9N6VUoY9Y72yBlyRJkpqIAV6SJElqIgZ4SZIkqYkY4CVJkqQmYoCXJEmSmogBXpIkSWoiBnhJkiSpiRjgJUmSpCZigJckSZKaiAFekiRJaiIGeEmSJKmJGOAlSZKkJmKAlyRJkpqIAV6SJElqIgZ4SZIkqYkY4CVJkqQmEimlRtdhyIiIF8aNG7fDQQcd1OiqSJIkqYX98Y9/ZP369ctSSjv291gDfImIeAKYBCwYoCJm5PNHBuj8Gjq81sOH13r48FoPH17r4aOR13ofYFVKad/+HmiAH0QRMQ8gpTS70XXRwPJaDx9e6+HDaz18eK2Hj2a91vaBlyRJkpqIAV6SJElqIgZ4SZIkqYkY4CVJkqQmYoCXJEmSmoij0EiSJElNxBZ4SZIkqYkY4CVJkqQmYoCXJEmSmogBXpIkSWoiBnhJkiSpiRjgJUmSpCZigJckSZKaiAG+ThGxR0R8JyIWRUR7RCyIiC9FxJR+nmeH/LgF+XkW5efdY6Dqrv6p91pHxISIeGtE/DAiHomItRGxOiLuiYh/iogxA/0eVJ2i/l2XnfPYiNgcESkiPl1kfVW7Iq91RMzK/30/k59rSUTcEhFnDUTd1T8F/r4+OiKuyo/fEBFPRcS1EfHagaq7qhcRb4qIr0TEbRGxKv8/97Iaz1X474Ki+CCnOkTEfsDvgV2Aq4BHgCOA44H5wFEppReqOM+O+XkOAG4C7gZmAGcAS4FXppT+PBDvQdUp4lrn/7lfBywDbgYeA6YApwNT8/OfmFLaMEBvQ1Uo6t912Tm3Ax4AdgImAp9JKX2syHqr/4q81hHxXuBCYDnwK2AhsAPwEuCZlNJfF/4GVLUCf1+/B/g6sBa4EngG2AN4AzAe+FhK6TMD8R5UnYi4HzgMWEN2fWYAl6eU3tbP8xT+u6BQKSWnGifgeiAB7ytb/9/5+m9WeZ5v5ft/sWz9+/P1v270ex3uUxHXGpgJvBUYU7Z+O2Befp5/avR7He5TUf+uy479DtkXt4/m5/h0o9+nU6H/h58EdObn267C9tGNfq/DfSro//DRwApgPXBg2baDgA3AOmBso9/vcJ7IAvb+QADH5df3skZ8ZgZysgW+Rvk3s8eABcB+KaXOkm3bAc+SfXh2SSmt7eU8E8la2TuBaSml1SXbRgB/BvbOy7AVvgGKutZ9lPEW4HLgmpTSaXVXWjUZiGsdEWcAvwDeDowCvost8A1X5LWOiD8A04G9UiNb5FRRgb+vdwUWAw+klA6rsP0B4BBgJz8HQ0NEHEf2F+9+tcAPxu/9etkHvnbH5/PflF5YgDyE307257RX9HGeVwDjgNtLw3t+nq4WndLyNPiKuta96cjnm+o4h+pX6LWOiF2Ai4FfpJRq6oOpAVPItY6IlwCHAr8BlkXE8RHxofy+lhPzhhg1VlH/rpcCzwEHRMT+pRsi4gCyVt/7De8tYTB+79fF/1hqd2A+/1MP2x/N5wcM0nk0cAbjGv1tPv91HedQ/Yq+1heT/T97Xj2V0oAo6lq/LJ8vBX5Hdh/TfwFfAH4L3B8R02uvpgpQyLVOWZeFfyD7Nz0vIr4XEf8ZEd8n6wb5MPDmAuqrxhvy2WxUowpuAZPz+coetnet336QzqOBM6DXKL/57bXA/WR9pdU4hV3riPhbshuU/yqltKT+qqlgRV3rXfL535HduHoKMBfYFfgE8DbgVxFxSEppY821VT0K+3edUvpJRCwCfgSUji60hKx7nF1dW8OQz2a2wEsNFBFvAL5E1q/yjSmljt6PUDOIiH3IrutPUko/bmxtNMC6fo+OBP46pXRtSmlVSulRsoB3D1kr3RsbVUEVJyLeRvaXldvIblwdn89vBL4KXNG42mk4McDXruvb1+QetnetXzFI59HAGZBrFBFnkv1nvxQ4zpuUh4SirvV3yEaq+PsC6qSBUdS17tq+OKV0R+mGvMvFVfnLI/pZPxWnkGud93P/DllXmbenlB5JKa1PKT1CdpP6PODN+Y2Tam5DPpsZ4Gs3P5/31P+p6waXnvpPFX0eDZzCr1FEvBn4CdmfXV+VUprfxyEaHEVd61lkXSueyx8ikiIikf2JHeDf8nW/qKu2qkfR/4ev6GH78nw+rrpqaQAUda1PIhtK8pYKNzZ2ArfmL2fXUkkNKUM+m9kHvnY35/OTImJEhSGGjiIbD/Z/+zjP/5K11B0VEdtVGEbypLLyNPiKutZdx7wV+B5Zf9njbXkfUoq61t8n+9N6uf2BY8nud5gH3FdvhVWzIv8PXwvsExETKgwp95J8/kQBdVZtirrWY/P5zj1s71rvvQ7Nr9Df+wPBFvgapZQeJxs2bB+yu9JLXQBMAH5Q+p95RMyIiBll51kD/CDff07Zed6bn/96Q17jFHWt8/XvIAt3TwHHel2HlgL/Xb8/pfTO8oktLfC/ytd9bcDejHpV4LVeB1wCtAGfjogo2f8Q4Gyy4WF/Wvy7UDUK/D/8tnz+pog4tHRDRMwE3kT2gJ+bCqu8BlREjM6v9X6l62v5zAw2H+RUhwqP2f0j8HKy8UP/BBxZOh5s/id0UkpRdp4d8/McQPYP/y6ym2LOIOsffWT+YVKDFHGtI+J4spufRpD1o3y6QlErUkpfGph3oWoU9e+6h3OfjQ9yGjIK/D98EnAL2dOW7yQbI3pX4A1kXWc+kFK6cIDfjnpR4LX+DnAOWSv7lcCTZCHvTGAM8KWU0gcH9t2oN/n9ZWfmL6cCryEbHajrC9jzKaUP5fvuQ/bXsSdTSvuUnadfn5lBV9QjXYfrBOxJ9gv5WbJ/0E+SjT4xpcK+ify+pgrbdgAuzI/fmJ/vO8AejX6PTsVca7KWuNTHtKDR79OpuH/XFfbt+gx8utHv0anYaw1MBD5D9ou9naxP/G+Akxr9Hp2Ku9ZkT988m2zM/+Vkf11ZRjYKzV83+j06Jch6M1T1e5bsy1ePv3v785kZ7MkWeEmSJKmJ2AdekiRJaiIGeEmSJKmJGOAlSZKkJmKAlyRJkpqIAV6SJElqIgZ4SZIkqYkY4CVJkqQmYoCXJEmSmogBXpIkSWoiBnhJkiSpiRjgJUmSpCZigJckSZKaiAFe0pAVEXMiIkXEcY2uS6mIuDSv1z6NrkuRIuLs/H2d3Y9j+v2ziIjj8mPm1FDN3s77vYhYGhETijxvq4mIBRGxoMZjd4+I9RHx6YKrJakfDPCSChMRIyPiXRFxS0Qsi4iOPFA9EBHfjojTy/bvd2DUFhFxdf7ze10P2+fn27/Xw/YL8u2fKLhe++TnvbTI8/ZR5suAtwOfTSmtLdvWFhEfiog7I2JlRGyMiGcjYl5EfDUiXjVY9Wx2KaWFwDeBf4yIPRtdH2m4GtXoCkhqDRExErgGeC2wAvgV8AwwBjgYeAswA/hlP077VeAK4Kki69pCbgROBU4ArivdEBF7AAcACTi+h+NPzOe/zedXAv8LPFt4TQfeZ4BVwDdKV0bEROAWYBawGPhZPp8IHAacC2yf76Pq/BfwPuDjZD8/SYPMAC+pKH9DFt7/ALwqpbSydGNEjAde3p8TppSeB54vrIat56Z8fkKFbV3rfgq8OSL2Tyk92rUx72ZyBLAauAsgv2Yry0801EXEAcCrgW+nlNaXbf4AWXj/DXBaSmlj2bFTgIMGo56tIqW0KCJuAN4SER8u/7cuaeDZhUZSUY7M55dW+oWeUlqXUrq563VE/A74bv7yu3mXi65pn3yfin3g83W/i4hdI+I7EbEkItZGxO8j4ph8nwkR8V8R8WREtEfEwxHx5vJ69dbPvj9dQfrq112p33FEjImI90fEvRGxPCLW5ftdFRGv7qtM4EHgOWBmHkRLnQCsAz5X8rrUMcBo4NaU0qa8Pj12aYqIV0fEbfnPeVlE/CIiZlTYbw7wRP7yHWXXtdJ5Z0bEryJiRf7+b4mII8v368PfAgH8T4VtXef6Rnl4B0gpLU8p/b5CvUZFxN9HxP9GxKq8bvdFxHsjouLvzog4IiL+JyIW5p+5ZyPiNxHxlxX2/cuIuDXv0rM+Ih6MiH+NiLEV9l2QT12f6afy8z8WEf8SEVHhmMjr+nBEbMjr9NWImNxD3fv7WbwCmAD8daXzSRpYtsBLKsoL+fyAKve/lKyrzRnAVcD9JdtWVHH89sDtZC3IPwJ2IAsT10fEK4Fv5euuIQuqfwP8T0Q8nVL63yrrONAuJavXQ8D3gfXAbsDRZH/N+G2PRwIppRQRNwN/SdZN5uclm08A5gL3kv0V40Syn0mXru4zN/ZVyYh4E1k43pjPn83reAfwQNnuvyO7NueT/TXmFyXb7i/b93Dgn/PzfBvYC3gjcGNEzEwpze+rbrlXA5vJuv+U6+/nkogYDVwNvAaYD/wQ2ED2M/4K2V+S3l52zLvIuu9sJusm9iiwC9l7/HvgxyX7/gfwr2TX5YfAGuB1wH8Ar4mIkyp82RgNXE/2+bgO2AScCXwWaAMuKNv/S8D7ya7VRUAH2b+1l5N1ays//6X077N4ez7/C7b+XEkaDCklJycnp7on4KVkoaAT+AHwBmDvPo45m6yP9tk9bJ+Tbz+ubH3Kp28CI0rWvz1fv4wsgLWVbDsm33ZlNWXk2/bJt11atv7SfP0+JeuOy9fN6eG9LAAWlLyenP+s7gFGVth/xyp/7ufm5X61ZN3++bqP5K9/AiwFomSfefk+h/Z2Pcj6ir9AFgAPLyv7/5Vci9KfRcWfW4Wf1TbXHnh3vv7rVb7/CWRh9sEetp+an68d+DpwCjCtj3N2fSa+UnptgJHAJfm2M0rWvzj/+SwDDq5wvj1Kll+ZH/8UMLVk/aj8M5uAj1b47CTgWmBcyfpdyL7srgBGl6w/Mt//MWCHkvVtZF+WUhGfRWA5sLSa6+Tk5FTsZBcaSYVIKd0HvA1Yks9/BiyIiBci4sqIOK3gItcBH04pdZas+yFZmJsCnJ9S2lBSv9vIgtDMgutRq0TW7aOdLDxtvTGlF7Y5orKuFvTSLjJdy1195G8GdgYOge5+3zPJut882Mf5zyD7S8YPU0r3lG2bQ3195m9PKV1atu47ZNfwiCrPsTtZsK54421K6RqyvwasB95D9heZRXn3lssj4tjS/fPuMe8ju9H1gymlzSXn2gz8E9m1e2vJYe8hC+D/nlJ6uEIdnil5+bf5/NMppcUl+2zKz90JvLOH9/r+VNLHP6W0lOyvV5OBA0v2OyeffyaltKxk/w1kLf/bVJHaPouLgZ0joq2H7ZIGiF1oJBUmpfTjiLiSrKvB0WSt8keT/an/zIj4PlmLayqguD+llFaXlb85IpYAE1JKf65wzEL6eSPtQEkprYqIq4HTgPsj4mfAbcCdKaV1/TjP4xHxFHBQRExLKT1LFuBXkbWyQxbgydc/QNYCPgK4qYprMSufbzNKS0ppZUTcD9Q6DGP5FwJSSh35NSzv09+THfP58p52SCl9OSK+Tdbd40iyz+WRZCMjvSUi/j2l1DWU5gFkX1geBT5WoXs5ZF8GSm98fUU+v67CvuW6fp43lW9IKf0pIp4B9o2IyWnre0lWppQeq3C+p/N56c+rx2tG1q1qc+mKOj6LXV8OdiIbcUrSIDHASypUSqmDbMSP30D38JJvJGtZPYtsqMJfFFBUTy2/m/rYNpT+3/sr4F/IgmRXH+YNEfFT4EMppSVVnudGslbXEyLih2RfoG7taj1OKf0xD8UnkvWNLh8+sjddNz32VJfFPayvxooe1m8ia1WvRleLdK+twHkQvSqfiIgxwLuAC4GPR8TPU0r3s+ULwf7AJ3s55cSS5e3z+cIq6tv18+xpqM5nye4F2J6tP8creth/Uz4v/Xn1eM1SSpsiotLITrV8Fsfl8/KRfyQNMLvQSBpQKaXNKaUfk/WXhspDHjZSV5eBSsF++4LOU/FcKaX1KaU5KaUDyELb28haSN9GNvxjtUqHkzyErLvMzWX7/A44Nv9C1XUN+ryBlS0hctcetk+tvpoDYmk+37HXvcqklDamlL5GdgM0bPmZdL3fK1NK0cu0b8npVuTz3asouuv8Pf3cppXtV4ser1lEjCJrMd9KjZ/FHcm+QCzrYbukAWKAlzRYurq7lPZJ6PpTfrWtrQOhq+tFpadKHl7EeSJiOltaRStKKT2dUrqcbOSTx4CjI6LaUFoa4Mv7v3e5GZhE1k3iIOCJlNIT9O3efL5NN5l8SMKZFY4ZzOv6LFlf/gP72rEH5Z/LR8gC+Svy0Wiq0TX6TcUn4pa5L58fV74h/5zsQXZtVlRZdiU9XjOyLm29XpdqPouRPSBrd+CBgrrESeoHA7ykQkTE30TEX1QaIzsippJ1VwC4tWRT181xew10/XpxVz4/J2+dBCCyx8R/ovIhFT1C1u/8jIjYpeQ844Avl+8cETtHxCEVzjOBrHvGJrYd6q+ilNKivPx9yG6SXEY2hGOprhb5T+XzalrfIetyspysr3j5F5o5VP5ispzsxsgBv655eLwV2CkPwFuJiPMi4hXbHgmRjWPf9WyAW/PzbSIbfWYa8OX8+pUfNy0iXlyy6htk1+vjZeu79t+j5OV38vnHImLnkn1GAl8g+718SQ9vt1qX5vN/i4gdSspoA/6zQv1q+Sy+jOyLQPlfeiQNgqHUF1RSc3s52WgfiyNiLlse5rMv2dB948jCYOmf4+8gG03mA3kLX1d/6q+kQXq6Y0rpzoi4FTgWuCsibiLrenAa2bjblVrmK52nIyIuJHu8/H35zbyjyG6cXJRPpXbP93uQ7MbSp8layE8l617x5fKbdPtwIzCDrAvNz8tbRfMbJBfl27v2r+Z9rYmIc8nGf78tIkrHgX8JWfA9tsIxdwLHRMTlwJ/Ix0dPKZWPG1+En5HdZ9HVYlzqtcA3InuI1u1kP+exZH3cX0M2vvqXU0p3lxzz78BhwHnAaflnYiHZsI37A0cB/wb8H0BK6f8i4u/JhjW9LyKuIrsJdkeyoLuK7L4EUkq/j4jPk41//1Dex3wtWev9S8i6rfxXPT+MlNLtEfEVstF0usroGgd+Odv2v6/ls3hSPv9ZPXWVVKNGjV/p5OTUWhNZ0P0HsptU55OFlo1kYeFasr60Iyoc91qyIL+GsjHF6X0c+N/1UI8FlIxxXbbtd+SNtmXrtwcuJutP3U72MJtz6cc48Pn6AD4CPJ6/96eAzwPjy+uVl/kJsq4uC/Nyn83r+DeUjNle5c//9SU/v/f2sM9l+fZOYJcK28+mh3H5yb6IzCX7wrWc7MvYjF5+FtPJxjV/IS+v+7z0c8z8Kt77GLIbNu+ssO0AsuEZryML92vzn/VTZA++OrWHcwbZcwVuJPuLxsb8Os0FPgrsWeGYV5IF2qX5/ouAXwNvqrDvX+fnWk32kKiHyb4UtPXn50HP/0YCeC/wx/z9LgK+RvYXk7o+i2R/JXgauL/o/0ecnJyqmyIlu65JkppbRPwr2ZNMZ6XsmQQaIPkzHX4JvD2ldFmj6yMNRwZ4SVLTy/t3zye7qbLoh4YpF9nA+PPIukQdkQwRUkN4E6skqeml7CmjbwfuiYgJja5PC5tK1vr+LsO71Di2wEuSJElNxBZ4SZIkqYkY4CVJkqQmYoCXJEmSmogBXpIkSWoiBnhJkiSpiRjgJUmSpCZigJckSZKaiAFekiRJaiIGeEmSJKmJGOAlSZKkJmKAlyRJkpqIAV6SJElqIgZ4SZIkqYn8f/Nfi5wIFRF9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 261,
       "width": 376
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Import numpy & curve fitting package\n",
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "# Import & configure plotting packages\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline \n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "# Define pulse\n",
    "pulse_width = np.linspace(0.01,1,100)\n",
    "rheobase = .5 # in V\n",
    "chronaxie = 0.1 # in ms\n",
    "\n",
    "amplitude = rheobase+((rheobase*chronaxie)/pulse_width)\n",
    "\n",
    "# Use our function & fit the curve!\n",
    "sample_popt, sample_pcov = curve_fit(func,pulse_width,amplitude)\n",
    "\n",
    "# Plot a regression line, using the slope & intercept\n",
    "plt.plot(pulse_width, func(pulse_width, *sample_popt), 'black')\n",
    "\n",
    "# Plot chronaxie & rheobase lines\n",
    "plt.vlines(chronaxie,amplitude.min(),amplitude.max(),color='blue',linestyle='--')\n",
    "plt.hlines(rheobase,pulse_width.min(),pulse_width.max(),color='red',linestyle='--')\n",
    "\n",
    "# Label the axes & show plot\n",
    "plt.ylabel('Stimulus Amplitude (Volts)')\n",
    "plt.xlabel('Stimulus Width (Seconds)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42a2181d",
   "metadata": {},
   "source": [
    "## Step 3. Fit a hyperbolic curve to *your* data\n",
    "\n",
    "Above, you can see a perfectly fit curve, with the <font color='red'><b>rheobase<b/></font> and  <font color='blue'><b>chronaxie</b></font>  labeled.\n",
    "    \n",
    "Now, let's do the same, but with your data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed40b786",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transform x & y data into arrays\n",
    "x = np.array(x)\n",
    "y = np.array(y)\n",
    "\n",
    "# Fit the curve with your data!\n",
    "popt, pcov = curve_fit(func,x,y)\n",
    "\n",
    "# Generate additional x-values to make a smoooth curve\n",
    "xnew = np.linspace(x.min(),x.max(), 100)\n",
    "\n",
    "# Plot a regression line, using the slope & intercept\n",
    "plt.plot(xnew, func(xnew, *popt), 'black')\n",
    "plt.scatter(x,y)\n",
    "\n",
    "# Optional -- can you figure out how to add rheobase & chronaxie?\n",
    "\n",
    "# Add your labels here!\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcd3e0e3",
   "metadata": {},
   "source": [
    "It is also important to know how well our curve fits. To do so, we'll calculate the residuals and then the sum of squares."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04b4ffad",
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Calculate the distance (on y axis) from data to curve\n",
    "residuals = y - func(x, *popt)\n",
    "\n",
    "# Calculate the residual sum of squares\n",
    "ss_res = np.sum(residuals**2)\n",
    "\n",
    "# Calculate the total sum of squares\n",
    "ss_tot = np.sum((y-np.mean(y))**2)\n",
    "\n",
    "# Calculate r-squared\n",
    "r_squared = 1 - (ss_res / ss_tot)\n",
    "\n",
    "print('R^2 = ',r_squared)\n",
    "print('Rheobase: ' , popt[0])\n",
    "print('Chronaxie:  ', popt[1])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}