{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 223,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x125da2160>]"
      ]
     },
     "execution_count": 223,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU1fn48c+TfRJCMlkIkIQESATBBSGgrFVxQaVorVWsVdyKVmv129/Xalur1m9trd/W2vrthkulFhU33HdEQVQkYd/DFkgISQhJCNln5vz+yAQjJjCZ9U7yvF+vvGbmzp25DzczDyfnPuccMcaglFIq/ESEOgCllFLe0QSulFJhShO4UkqFKU3gSikVpjSBK6VUmIoK5sHS0tJMbm5uMA+plFJhr6io6IAxJv3o7UFN4Lm5uRQWFgbzkEopFfZEpKSr7dqFopRSYUoTuFJKhSlN4EopFaY0gSulVJjSBK6UUmFKE7hSSoUpTeBKKRWmNIEHwYayOpZuqwp1GEqpXkYTeIDVNrZy3dMrufvldaEORSnVywR1JGZf9MCbm6iqbyFCoM3pIjpS/89USvmHZpMA+mhLBa+sKmN4egIuA/vrmkMdklKqF9EEHiCHmtv4xSsbGJGRyD0XjQKgtKYpxFEppXoT7UIJkAff3EzV4RbmXTOO/nHRAJTVagJXSvmPJvAAWLqtioWFe/nRmcM5JSuZFocTgDJtgSul/Ei7UPzscIuDn7+ynrwB/bh9ej4AsVGRDEiMpay2McTRKaV6E22B+9nv3t5MeV0TL/1oEnHRkUe2Z9pt2oWilPIrbYH70WfbD7BgxR5umDKUsUPsX3suM9mmXShKKb/SBO4nDS0O7nplHUPTEvh/5434xvOZdhv7aptxuUwIolNK9UaawP3kf9/bSmlNEw9fdsrXuk46ZCXbaHW6OHC4JQTRKaV6I03gfvDlroM8/dlu5kzMZXxuSpf7ZNptAOzVbhSllJ9oAvdRU6uTn720luwUGz+b8c2ukw6ZyfGA1oIrpfznuAlcREaIyJpOP4dE5A4RSRGRD0Sk2H1rP9579UZ/fH8ru6sb+f13TyE+pvuino4WuF7IVEr5y3ETuDFmqzFmjDFmDDAOaAQWAXcDi40x+cBi9+M+paikhieX7+Kq04cwaXjaMfftFxtFki1aa8GVUn7T0y6U6cAOY0wJcDEw3719PnCJPwOzuua29q6TwUk2fn7hiR69RksJlVL+1NMEPht4zn0/wxhT7r6/H8jo6gUiMldECkWksKqq9yxq8OiHxeyoauB3l55Mv1jPxkPpYB6llD95nMBFJAaYBbx49HPGGAN0WeBsjJlnjCkwxhSkp6d7HaiVrN1by7ylO7iiIJtpJ3j+b+pogbefLqWU8k1PWuAXAKuMMRXuxxUiMgjAfVvp7+CsqMXh5M6X1jIgMY5fzvSs66RDlt1GQ6uTuqa2AEWnlOpLepLAr+Sr7hOA14E57vtzgNf8FZSVPbZ4O9sqDvO7S08+Mk2sp7LclSg6L7hSyh88SuAikgCcC7zSafNDwLkiUgyc437cq20oq+Pvn+zg0rGZnDVyQI9fr7XgSil/8ujqmzGmAUg9als17VUpfUKrw8V/v7iWlIQY7p05yqv30FpwpZQ/6XSyHvrbx9vZsr+eeVePIzk+xqv3sMdHY4uO1Ba4UsovdCi9BzaXH+L/PtrOrFMHc97ogV6/j4i0lxJqC1wp5QeawI/D6TLc/fI6kmzR3D9rtM/vl5msteBKKf/QBH4cz365h7Wldfxq5ihSErzrOuks026jtEaH0yulfKcJ/Biq6lt4+N0tTBqeysVjBvvlPTOTbdQ0ttHY6vDL+yml+i5N4Mfw27c309zm5IGLT0JE/PKeWVqJopTyE03g3fhsxwEWrS7jpmnDyRvQz2/vm5nsHsyj/eBKKR9pAu9Cq8PFr17dQHaKjR+fnefX99ZacKWUv2gdeBceX7aTHVUN/Ova8V2ub+mLAYlxREWIVqIopXymLfCj7Klu5C+Li5kxeqBXw+WPJzJCGJQcpy1wpZTPNIF3Yozhvtc3EBUh3DfLu+HyntBacKWUP2gC7+S9jftZsrWK/zr3BAYl2QJ2nMzkeG2BK6V8pgncraHFwa/f2MTIgYlcOyk3oMfKstuoqG+m1eEK6HGUUr2bJnC3Rz/cRnldMw9+5ySiIgN7WjLtNoyB/XXNAT2OUqp30wRO+2RVTy3fzezx2YzLSQn48bKO1ILrkHqllPf6fAJ3uQz3vLqBJFs0d80YGZRjai24Usof+nwCf7FoL0UlNdx9wUjsfpisyhODkmyI6Mo8Sinf9OkEfrChld+9s4XxuXYuG5sVtOPGREUwIDFW18ZUSvmkTyfwh97ZzOFmB7+55GQiIvwzWZWnMpN1YQellG/6bAIv3H2QFwpLuWHKUEYMTAz68TPt8dqFopTySZ9M4G1OF79ctIHBSXH8ZHp+SGLITLZRXteEy2VCcnylVPjzKIGLSLKIvCQiW0Rks4hMFJEUEflARIrdt/ZAB+sv/1q+i60V9dw3azQJsaGZzyvTbqPNaaisbwnJ8ZVS4c/TFvifgXeNMSOBU4HNwN3AYmNMPrDY/djyymqbePTDYqaPHMB5ozJCFkdHLXiZ1oIrpbx03AQuIknANOBJAGNMqzGmFrgYmO/ebT5wSaCC9Kdfv74RlzHcP2u031bZ8UZHLbhWoiilvOVJC3woUAX8S0RWi8gTIpIAZBhjyt377Ae6bM6KyFwRKRSRwqqqKv9E7aXFmyt4f1MFt52dT3ZKfEhjyTzSAtcErpTyjicJPAoYC/zdGHMa0MBR3SXGGAN0eTXOGDPPGFNgjClIT0/3NV6vNbU6ue/1jeQN6McPpw4LWRwdEmKjSI6P1lJCpZTXPEngpUCpMWaF+/FLtCf0ChEZBOC+rQxMiP7x2EfFlNY08ZtLTiImyhrFNzovuFLKF8fNZMaY/cBeERnh3jQd2AS8Dsxxb5sDvBaQCP1ge2U9jy/byaVjMzljWGqowzkiy66DeZRS3vO0hu42YIGIxAA7getoT/4viMgNQAlweWBC9E1Tq5OfvrCW+JgofnHhiaEO52syk+NZVnwAY0xIL6gqpcKTRwncGLMGKOjiqen+Dce/XC7DHQtXs76sjnlXF5DWLzbUIX1Npt1GY6uT2sa2oE2kpZTqPazRGRwgD727hfc2VvCri0ZxbghrvrujlShKKV/02gS+YEUJ85bu5JqJOVw3OTfU4XQp60gtuA7mUUr1XK9M4J9sq+Le1zZy1oh07p05yrL9yx0tcB3Mo5TyRq9L4Fv2H+LWBas4ISORx74/NuDrW/oiOT6a+JhI7UJRSnnFutnNC5WHmrn+XytJiI3kqWsL6Beiiao8JSI6L7hSymvWznA90Njq4MZ/F1Lb1MYLN01kUJIt1CF5JNOug3mUUt7pFS1wp8tw+/Nr2FBWx19mn8ZJmUmhDsljOhpTKeWtXpHAf/f2Zj7YVMGvZo7iHAuWCx5Lpt1GbWMbDS2OUIeilAozYZ/An/mihCc+3cW1k3K5bvLQUIfTY1oLrpTyVlgn8CVbK7nvtQ2cPXIAv5o5KtTheKWjFlwvZCqleipsE/jm8kP8eMEqRg7sz2NXnkZkkFeV95fM5PZ5yUu1Ba6U6qGwTOAVh5q5/umV9IuL4slrC0K2rqU/DEiMJTpStAWulOqxsMt8ja0Obpi/krowKxfsTkSEMFgrUZRSXgirFnhHueCmfYd47MrwKhc8lvbBPDofilKqZ8Iqgf/WXS5478xRTD8xvMoFj0VrwZVS3gibBP7M57t50l0ueG0YlgseS6bdRsWhFloczlCHopQKI2GRwJdsqeS+1zdyzonhWy54LB214OW1zSGORCkVTiyfwI0xPLV8FycO6s+fZ4dvueCxZNp1MI9SqucsX4UiIjx+TQH1zY6wLhc8lix3LbiWEiqleiIsMmJcdCRx0ZGhDiNgBibFIaKDeZRSPWP5LpS+ICYqgozEOG2BK6V6xKMWuIjsBuoBJ+AwxhSISAqwEMgFdgOXG2NqAhNm79c+L7jWgiulPNeTFvhZxpgxxpgC9+O7gcXGmHxgsfux8pLWgiulesqXLpSLgfnu+/OBS3wPp+/KtNsor23G6TKhDkUpFSY8TeAGeF9EikRkrntbhjGm3H1/P9Dl0EgRmSsihSJSWFVV5WO4vVdmsg2Hy1BZr7XgSinPeJrApxhjxgIXALeKyLTOTxpjDO1J/huMMfOMMQXGmIL09HTfou3FMnVecKVUD3mUwI0xZe7bSmARMAGoEJFBAO7bykAF2Rdk6co8SqkeOm4CF5EEEUnsuA+cB2wAXgfmuHebA7wWqCD7go4WeKm2wJVSHvKkjDADWCQiHfs/a4x5V0RWAi+IyA1ACXB54MLs/eJjokhJiNEWuFLKY8dN4MaYncCpXWyvBqYHIqi+KjPZpi1wpZTHdCSmhejCDkqpntAEbiHtozGbaC/qUUqpY9MEbiGZyTaa21wcbGgNdShKqTCgCdxCdF5wpVRPaAK3kI6VeXQwj1LKE5rALSRLW+BKqR7QBG4hSbZoEmIitZRQKeURTeAWIiJHKlGUUup4NIFbTHstuCZwpdTxaQK3GG2BK6U8pQncYjKT46lrauNwiyPUoSilLE4TuMXovOBKKU9pAreYI7XgvWiB4zanS/+iUCoANIFbTHYvmxf8sx0HOPeRT5jx6FKd40UpP/NkPnAVRGn9YomJjAj7LpTaxlZ++/ZmXigsxRYdSVObkx1VDeQN6Bfq0JTqNbQFbjEREcLg5DhKw7QSxRjD62v3cc4jn/DyqjJu/tZwXrx5IgBFJQdDHJ1SvYu2wC0o0x6eteBltU3cs2g9S7ZWcUpWEvOvn8DowUkYY7DHR1O4u4Yrxg8JdZhK9RqawC0oM9nGkq1VoQ7DY06XYf5nu/nD+1sB+NXMUVw7KZfICAHaR5iOy7FTVFITyjCV6nU0gVtQZnI8VfUtNLc5iYuODHU4x7Rp3yF+/so61pbWceaIdH5zyUlk2eO/sd+4nBQ+3FzJwYZWUhJiQhCpUr2PJnAL6qgFL69rZmhaQoij6Vpzm5M/Ly5m3tKdJNui+fPsMcw6dTDuxa+/oSDXDkBRSQ3njsoIZqhK9VqawC2o87zgVkzgy7cf4BeL1lNS3cj3xmXxy4tOJDn+2K3qkzOTiI4UCksOagJXyk88TuAiEgkUAmXGmJkiMhR4HkgFioCrjTG6FpgffDUvuLUG89Q0tPLg25t5qaiU3NR4nr3xdCblpXn02rjoSE7KTKJot/aDK+UvPSkjvB3Y3Onx74E/GWPygBrgBn8G1pcNTIojQqwznN4Yw2tryjjnkU94dXUZt5w5nHfvmOZx8u5QkGNnXVkdLQ5ngCJVqm/xKIGLSBZwEfCE+7EAZwMvuXeZD1wSiAD7oujICDL6W6MWvL65jVsWrOL259eQZbfxxm1T+NmMkV5dXB2Xk0Krw8WGsroARKpU3+NpF8qjwM+ARPfjVKDWGNMxwUUpkNnVC0VkLjAXYMgQrQH2lBXmBS+uqOemZ4ooOdjIzy8YyY1Thx0pDfTGuJz2C5mFu2sYl5PirzCV6rOO2wIXkZlApTGmyJsDGGPmGWMKjDEF6enp3rxFnxTqecHfXLePi/+6nEPNDp698XRu+tZwn5I3QHpiLLmp8VoPrpSfeNICnwzMEpELgTigP/BnIFlEotyt8CygLHBh9j2ZyTbeWleOw+kiKjJ4Mx60OV089M4Wnvx0F+Ny7PztqrFk9I/z2/uPy0nh462VGGO6LTlUSnnmuJnBGPNzY0yWMSYXmA18ZIy5ClgCXObebQ7wWsCi7IMy7TYcLkNFfUvQjllZ38xVj6/gyU93ce2kXJ774Rl+Td7Q3o1S3dDK7mprVdgoFY58adrdBfxURLbT3if+pH9CUsCR0YzB6gcv3H2QmX/5lHVltTx6xRjunzWamCj/t/w7BvQU7taJrZTyVY8G8hhjPgY+dt/fCUzwf0gKjl7YIXAX/IwxPP3Zbh58azNZdhvzr5/AiYP6B+x4een96B8XRVFJDd8ryA7YcZTqC3QkpkV1Ho0ZKI2tDn7+ynpeW7OPc04cwB8vH0OSLTpgx4P26XLH5dgp1AuZSvlME7hF2WIiSU2ICVglyq4DDdz8TBHbKuu58/wR/Ohbw4nwscrEUwW5KSzZupXaxtbjDsFXSnVPF3SwsEy7LSBLq72/cT+zHvuUyvpm5l83gVvPygta8oav6sFX7dFWuFK+0ARuYZnJ/q0Fd7oM//veFuY+U8TQ9ATeuG0K004Ifm3+qVnJREUIhTovilI+0S4UC2tf2ME/NdMHG1r5yXOr+XT7AWaPz+b+WaNDNte4LSaS0YP7az+4Uj7SBG5hmXYbzW0uqhtaSesX6/X7rC+t4+b/FFF1uIWHLj2Z2RNCP6XBuJwUFqwoodXhCki5olJ9gX5zLMwflSgfbKrg8n9+DsBLN0+0RPKG9nrwFoeLjft0YiulvKUJ3MIyj8wL7l0C//fnu7npmULyM/rx6q2TOSUr2Y/R+aYg56sVepRS3tEEbmFZyd6NxnS5DA++tYl7X9vI2SMH8PzcM0hP9L4LJhAG9I8jO8WmCVwpH2gfuIX1t0XRLzaqRy3w5jYnP31hDW+v38+ciTnc++3RPs8iGCgFOSl8uv2ATmyllJe0BW5hIkJmso3SGs8mfjrY0MpVT6zg7fX7ueeiE7l/lnWTN8DYHDtV9S3sPRj6hSuUCkfaArc4Twfz7D7QwHVPr6Sstom/XTWWC08eFITofNPRD15YcpAhqfEhjkap8KMtcIvzZDBPUUkNl/79M2obW3nuh6eHRfIGOCEjkcTYKK0HV8pLmsAtLstuo77ZwaHmti6ff3dDOd9//Av6x0Xxyi2Tw2qpssgI4bQce69fqb6uqY0dVYdDHYbqhTSBW9yRUsKjulGMMTyxbCc/WrCK0YP78/KPJjE0LSEUIfqkIMfOtsp66pq6/g8qnLlchoUr93DWHz7mgkeXUXmoOdQhqV5GE7jFdTWYx+ky/PqNTfzmrc3MGD2QZ394Bqk+jNQMpYIcO8bA6l42sdWavbV852/Luevl9WTbbbQ6Xby0qjTUYaleRhO4xR09mKep1cnN/yni6c92c+OUofz1+2NDNqeJP5yanUxkhPSaevADh1u466V1XPLX5ZTXNfPoFWN49dbJnD40hYUr9+JymVCHqHoRrUKxuLSEWGKiIiirbeLA4RZumF/I+tJafj1rNHMm5YY6PJ8lxEZx4qDEsJ+Z0OF08cwXJTzywTaaWp3cNG0Yt03Pp19s+1fsyglDuGPhGr7YVc2k4Wkhjlb1FprALS4ior0WfOXug7yzoZyq+hb+eXUB547KCHVoflOQ0946bXO6iI4Mvz8Kv9hZzf2vb2TL/nqm5qdx37dHkzeg39f2mXHSQPq/FsXClXs1gSu/Cb9vSx+UmWxj9Z5amlqdPD93Yq9K3tC+wENTm5PN5YdCHUqPlNc1cdtzq5k97wvqmx384wfj+Pf1E76RvAHioiP5zmmZvLNhP7WNrSGIVvVGmsDDwOlDUxg5MJFFt0xmTLZ1JqTyl69Wqg+PbpQWh5O/fbyd6X/8hPc27uf26fl8+NNvMeOkgcecEuCK8UNodbhYtLosiNGq3uy4XSgiEgcsBWLd+79kjLlPRIYCzwOpQBFwtTFGmxYBcNv0fG6bnh/qMAJmUJKNzGQbRXtquJ6hoQ7nmJZsreSBNzax60AD543K4FczR5Gd4tko0lGD+3NqVhLPf7mXayfl6vwvymeetMBbgLONMacCY4AZInIG8HvgT8aYPKAGuCFwYarebqx7QI8x1qzS2FPdyI3zC7nuXysRYP71E5h3TYHHybvDFeOHsLWinrWlOg+68t1xE7hp1zGMLNr9Y4CzgZfc2+cDlwQkQtUnFOTY2X+o2a9rgPpq78FGnv9yD7c9t5pz/vQJn+04wN0XjOTdO6bxLS/XEv32qYOwRUfy/Jd7/Byt6os8qkIRkUjau0nygL8CO4BaY4zDvUspkNnNa+cCcwGGDLHGajDKesZ1WuAhyx6aia0OHG7hsx3VfLb9AMt3HDgyS2J6YizfHZvJ7dNPYGBSnE/HSIyLZuYpg3h97T7umTnqSJmhUt7w6NNjjHECY0QkGVgEjPT0AMaYecA8gIKCAmv+faxCbuTARBJiIincXcPFY7psC/jd4RYHX+6qZvn2apZvP8CW/fUAJMZFccawVG6YPJTJeWnkDejn1/7q2ROG8GJRKW+t28cV47VRo7zXo//+jTG1IrIEmAgki0iUuxWeBeildeW1qMgIThtiD+jMhK0OF6v31LB8R3vCXru3FofLEBMVQUGOnTvPH8HkvDROGtyfqADWo48dkkz+gH489+XeXp/AN+07xH9WlLCnupEn5hSE9ahhK/KkCiUdaHMnbxtwLu0XMJcAl9FeiTIHeC2Qgareb1yOncc+Kqa+uY3EuGi/ve+60lr+8P42Vu46SFObkwiBk7OSmTttGJPz0hiXYw9qYhERrhifzW/e2syW/YcYObB/0I4dDM1tTt7ZUM5/vthDUUkNEQIuA1/uOsg0L68dqK550gIfBMx394NHAC8YY94UkU3A8yLyG2A18GQA41R9QEGuHZdpnwhqar5/vuiNrQ5uWbCK5jYXlxdkMTkvjdOHpZJk899/EN64dGwWD7+7lYUr93Lft0eHNBZ/Kalu4NkVe3ihcC81jW0MTUvgnotO5KJTBvGthz9mWXGVJnA/O24CN8asA07rYvtOYEIgglJ905jsZCKkfUCPvxL4I+9vo7SmiRdumsiEodaZKz0lIYbzRmewaHUZd80YGbZdCw6ni4+2VPKfFXtYuq2KyAjh3BMz+MEZOUwankqEe0m/cTl2lhUfCHG0vY9eAleWkRgXzYiB/f02M+G60lqeWr6LKycMsVTy7jB7/BDeXFfOexv3B+3Crb9U1jez8Mu9PPflHvbVNZPRP5Y7zsln9vghXVbqTD0hjYff3UplfTMDEn2r5FFf0QSuLKUgx84rq0pxOF0+XUhsc7q4++X1pPWL5e4LPC6aCqpJw1PJTrGxcOXesEjgxhg+31nNgi/28N7G/Thchqn5adz77dGcc+KAY/6+pual8zBbWb79AN85LSuIUfdumsCVpRTk2nnmixK27K/npMwkr9/nyU93san8EH+/amzI+7u7ExEhXFGQzR/e30ZJdQM5qdZcUam5zcmzK/awYEUJO6oaSLJFc93kXL5/eo7Hq0CNHtwfe3w0y7ZpAvcnncxKWUrnAT3eKqlu4E8fbOPcURnMOGmgv0ILiMvGZRMhsHDl3lCH0iWXy/DjZ1fzwJubSIyL5g/fO5UVv5jOLy8a1aMl/CIihMl5aXy6/YBlp0sIR5rAlaVkJtsY2D/O63pwYwy/WLSe6MgI/ufikyw/YdTApDjOGjGAF4vau42s5pEPtvHh5grunTmKV2+dzGXjsry+4DotP53K+ha2VegCz/6iCVxZiogwLsfOKi8T+Muryli+vZq7Zozwedh7sFwxPpuq+haWbK0KdShf88baffzfku3MHp/NdZNzfX6/KfntC1ksK7bWvzOcaQJXljMux05ZbRPldT2b2OrA4RZ+89YmxuXYuer0nABF539njxzAgMRYFq60zgRX60vruPOltYzPtfOAn/6SGZxsY3h6gpYT+pEmcGU53i7w8D9vbqKhxcFDl558pP44HERFRnDZuCw+2lLJ/rrmUIdDZX0zc58pJCU+hr//YBwxUf5LE1Pz01mxq5rmNqff3rMv0wSuLOfEQf2xRUf26ELmkq2VvLZmH7ecmUd+RmIAowuMywuycRl4qSi0FzNbHE5ufqaI2sY2Hp9TQFq/WL++/9T8NJrbXF53kamv0wSuLCc6MoIx2ckUlhz0aP+GFgf3LNpA3oB+3HLW8ABHFxi5aQlMHJbKwsK9uFyhqdIwxvDLRRtYtaeWP3zvVEYP9r6MszunD0slKkJYqt0ofqEJXFlSQa6dzeX1NLQ4jrvvIx9so6y2id9dejKxUeE5JB1g9oRs9h5s4rMd1SE5/lPLd/NSUSk/mZ7PRacMCsgx+sVGMTbHzqfb9UKmP2gCV5Y0NseO02VYu7f2mPut3VvLv5bv4qrThzA+13rD5Xvi/NEDSbJF83wILmYu3VbFg29t4vzRGdwR4PVXp+alsaHsENWHWwJ6nL5AE7iypLFD7IhwzHrwNqeLu15eR3piLHdZdLh8T8RFR/Kd0zJ5f2MFBxuCtz74zqrD/PjZVZyQkcgjl48J+AXgqe4ZCZeH6C+N3kQTuLKkJFs0JwxIPGYCf3zZTrbsr+fXs06ivx/nDw+l2ROyaXW6WLQ6OOujHGpu48Z/FxIVGcHj1xSQEIQl3k7OTCLJFs2ybdqN4itN4MqyxuXaWV1Sg7OLi3q7DzTw5w+LmTF6oOWHy/fEyIH9GZOdzPNf7gn4kHOny/CT51azp7qRv101luyU4KxFGhkhTBqeqsPq/UATuLKsghw79S0OtlXUf217x3D5mKgIfn1x71gMobPZ47MprjzMqj3H7v/31cPvbuHjrVU8cPFJnDEsNaDHOtrU/HTK65rZUaXD6n2hCVxZVkFO+0XJo7tRXiwq5bMd1dx9wUgy+ofHcPmemHnqYOJjIgM6MvOVVaX8c+lOrpmYw/dPD/66nFOPDKvXckJfaAJXlpWdYiOtXyxFu7+qB6+qb+HBtzYzITeFK3vpgsD9YqOYdepg3lhbTn1zm9/ff/WeGu5+ZT0Th6Xyq5mj/P7+nshOiSc3NZ5PNYH7RBO4siwRoSDHTtGer1rgD7y5iaZWJ78Ns+HyPXXF+Gya2py8sbbcr++7v66Zm54pIqN/LH+7aizRPiya4asp+Wl8vrOaVof1ZmEMF5rAlaUV5NrZe7CJykPNLNlSyRtr93HrWXnkDegX6tACakx2MiMyEv3ajdLc5mTuM4U0tDh44prx2BNi/Pbe3pian05jq5PVe3RYvbc0gStL61jg4ZNtVdzz6gbyB/TjR2eG53D5nhARrhifzdrSOjbtO+Tz+xljuOvldawvq+PR2acxYuNbkLgAAA0gSURBVGDo54uZODyVyAjRfnAfHDeBi0i2iCwRkU0islFEbndvTxGRD0Sk2H1rD3y4qq8ZPTiJ2KgIfv3GJvbVNfHQd0/26+x4Vnbp2ExioiL80gr/xyc7eW3NPv77vBGcOyrDD9H5rn9cNGOyk1m2XRO4tzz5JjiA/2eMGQWcAdwqIqOAu4HFxph8YLH7sVJ+FRMVwanZyRxucXD1GTmMywnv4fI9kRwfw4zRA1m0uqxH06/WNrayak8NLxeV8sf3t3LLgiIefm8LM08ZxC0W++tlSl4a60prqW0M3sjT3uS4w66MMeVAuft+vYhsBjKBi4Ez3bvNBz4G7gpIlKpPO3/0QKoPt3Dn+SNCHUrQzR6fzetr9/Huhv1cctpXK9c3tDjYXd3ArgMN7KpqYFfH/QMN1DZ+VbkSGSFk2W18Z0wmD37nZMstMTfthDT+vLiYz3ZUc+HJgZlAqzfr0bhZEckFTgNWABnu5A6wH7DG32Wq17lhylCun5xrueQTDGcMS2VISjx/+3g7K3ZVH0nSFYe+PhHUwP5xDE1L4MKTBzEsLYHc1ASGpieQbY+3dJfTqVnJJMZGsay4ShO4FzxO4CLSD3gZuMMYc6jzl8kYY0SkyzGxIjIXmAswZEjvrNtVgdcXkze0r+Y+Z1Iu//PmJg4cbiU3NZ4peekMS3cn6bQEctPiiY8J/BwmgRAVGcHE4aks3dY+rL6v/p695dFvXUSiaU/eC4wxr7g3V4jIIGNMuYgMAiq7eq0xZh4wD6CgoEAnPlCqh66fnMvs8dlBmWgqFKbmp/H+pgp2VzcyNC0h1OGEFU+qUAR4EthsjHmk01OvA3Pc9+cAr/k/PKWUiPTa5A3t9eAAn+pq9T3mSefYZOBq4GwRWeP+uRB4CDhXRIqBc9yPlVKqR3JS48my23SZNS94UoXyKdBdx9R0/4ajlOprRISp+Wm8ubYch9NFVAiH94cbPVNKqZCbmp9OfYuDtaWBnUK3t9EErpQKuUnDUxGBpdt6ZzdKIGaVBE3gSikLSI6P4ZSsZD7tZcPqiyvqufPFtUx4cDGlNY1+f//ee2lbKRVWpual8fdPdnCouS3s1zhdufsg//xkBx9uriQuOoLZ44cQE4C+fU3gSilLmJqfxv8t2c7nO6o5f3T4rXPqchk+3FzBP5fupKikhpSEGP7rnBO4emIOKQGaulcTuFLKEk4bYic+JpJlxVVhlcBbHE5eXV3GP5fuZGdVA9kpNh64eDTfG5eNLSYyoMfWBK6UsoSYqAgmDksNm2XWDjW3seCLPTy1fBdV9S2MHtyfx648jQtOGhi0UkhN4Eopy5iSn8biLZXsPdhIdkp8qMPp0v66Zv61fBcLVuzhcIuDqflp/OnyMUzOSw36XC6awJVSltExrH5Z8QG+f7q1Jr8rrqhn3tKdvLqmDKfLMPOUwcydNoyTMpNCFpMmcKWUZQxPT2BQUhzLiqv8msAPNrTS0OKgxeGixeGk1eGixeE66tb5tfudnyuuPMxHW9orSr4/YQg3Th1mib8QNIErpSyjY1j9uxv243QZIiN865IwxvDbtzfz+LJdXsYDMZER2ONjuOOcfK6ZmBuwihJvaAJXSlnKlPx0XigsZV1pLacN8X6pXafLcM+rG3juyz1cNi6LCUNTiI2KIDYq0n0bQYz7cYz7cWx0BDGREcRGRxITGUF0pFh6jnJN4EopS5mSl4YIfFp8wOsE7nC6+O8X1/Lqmn3ccuZw7jx/hKUTsbd0KL1SylJSEmIYPbi/16vVtzic3PrsKl5ds487zx/Bz2aM7JXJGzSBK6UsaEpeOqtKajjc4ujR65panfzw30W8t7GC+749ilvPygtQhNagCVwpZTnT8tNwuAwrdlZ7/Jr65jbm/OtLlhVX8fvvnsx1k4cGMEJr0ASulLKccbl24qIjWObhqMzaxlZ+8MQKikpqePSKMVwx3lo15IGiFzGVUpYTGxXJ6UNTWebBOplV9S1c/eQKdlY18PerxnJeGM2j4ittgSulLGlqfho7qhrYV9vU7T77apu44p+fU1LdyJPXFvSp5A2awJVSFvXVavVdd6OUVDfwvX98TlV9C/++YcKR/fsSTeBKKUs6IaMfAxJjWdpFN0pxRT3f+8fnNLQ6ePaHZzA+NyUEEYaeJnCllCWJCFPy01i+/QAulzmyfUNZHVfM+wKXgYVzJ3JyVugmkwq14yZwEXlKRCpFZEOnbSki8oGIFLtvvR/vqpRS3Zian0ZNYxsb9x0CoKikhisf/4K4qAhevHkiIwYmhjjC0PKkBf40MOOobXcDi40x+cBi92OllPKryXlpACzbXsVn2w9w9ZMrSE2I4YWbJzI0LSHE0YXecRO4MWYpcPCozRcD89335wOX+DkupZRiQGIcIwcm8tyXe7j26ZVk2W28cNNEsuyhn8rVCrztA88wxpS77+8HMrrbUUTmikihiBRWVR2/plMppTqbdkI6ew82cUJGP56fO5EB/eNCHZJl+DyQxxhjRMQc4/l5wDyAgoKCbvdTSqmuzJmUiwjcelYe/eOiQx2OpXjbAq8QkUEA7ttK/4WklFJfyUy28fMLTtTk3QVvE/jrwBz3/TnAa/4JRymllKc8KSN8DvgcGCEipSJyA/AQcK6IFAPnuB8rpZQKouP2gRtjruzmqel+jkUppVQP6EhMpZQKU5rAlVIqTGkCV0qpMKUJXCmlwpQmcKWUClNiTPAGR4pIFVDi5cvTAM8WyAsNjc83Gp9vND7fWD2+HGPMN1asCGoC94WIFBpjCkIdR3c0Pt9ofL7R+Hxj9fi6o10oSikVpjSBK6VUmAqnBD4v1AEch8bnG43PNxqfb6weX5fCpg9cKaXU14VTC1wppVQnmsCVUipMWS6Bi8gMEdkqIttF5BuLJYtIrIgsdD+/QkRygxhbtogsEZFNIrJRRG7vYp8zRaRORNa4f+4NVnzu4+8WkfXuYxd28byIyF/c52+diIwNYmwjOp2XNSJySETuOGqfoJ4/EXlKRCpFZEOnbSki8oGIFLtv7d28do57n2IRmdPVPgGK739FZIv797dIRJK7ee0xPwsBjO9+ESnr9Du8sJvXHvO7HsD4FnaKbbeIrOnmtQE/fz4zxljmB4gEdgDDgBhgLTDqqH1uAf7hvj8bWBjE+AYBY933E4FtXcR3JvBmCM/hbiDtGM9fCLwDCHAGsCKEv+v9tA9QCNn5A6YBY4ENnbY9DNztvn838PsuXpcC7HTf2t337UGK7zwgyn3/913F58lnIYDx3Q/8twe//2N+1wMV31HP/xG4N1Tnz9cfq7XAJwDbjTE7jTGtwPPAxUftczEw333/JWC6iEgwgjPGlBtjVrnv1wObgcxgHNuPLgb+bdp9ASR3LI8XZNOBHcYYb0fm+oUxZilw8KjNnT9j84FLunjp+cAHxpiDxpga4ANgRjDiM8a8b4xxuB9+AWT5+7ie6ub8ecKT77rPjhWfO29cDjzn7+MGi9USeCawt9PjUr6ZII/s4/4Q1wGpQYmuE3fXzWnAii6enigia0XkHREZHdTAwADvi0iRiMzt4nlPznEwzKb7L04ozx9AhjGm3H1/P5DRxT5WOY/X0/4XVVeO91kIpB+7u3ie6qYLygrnbypQYYwp7ub5UJ4/j1gtgYcFEekHvAzcYYw5dNTTq2jvFjgVeAx4NcjhTTHGjAUuAG4VkWlBPv5xiUgMMAt4sYunQ33+vsa0/y1tyVpbEfkl4AAWdLNLqD4LfweGA2OActq7KazoSo7d+rb8d8lqCbwMyO70OMu9rct9RCQKSAKqgxJd+zGjaU/eC4wxrxz9vDHmkDHmsPv+20C0iKQFKz5jTJn7thJYRPufqp15co4D7QJglTGm4ugnQn3+3Co6upXct5Vd7BPS8ygi1wIzgavc/8l8gwefhYAwxlQYY5zGGBfweDfHDfX5iwIuBRZ2t0+ozl9PWC2BrwTyRWSou5U2G3j9qH1eBzqu+F8GfNTdB9jf3H1mTwKbjTGPdLPPwI4+eRGZQPs5Dsp/MCKSICKJHfdpv9i14ajdXgeucVejnAHUdeouCJZuWz6hPH+ddP6MzQFe62Kf94DzRMTu7iI4z70t4ERkBvAzYJYxprGbfTz5LAQqvs7XVL7TzXE9+a4H0jnAFmNMaVdPhvL89Uior6Ie/UN7lcQ22q9Q/9K97QHaP6wAcbT/6b0d+BIYFsTYptD+5/Q6YI3750LgZuBm9z4/BjbSflX9C2BSEOMb5j7uWncMHeevc3wC/NV9ftcDBUH+/SbQnpCTOm0L2fmj/T+ScqCN9n7YG2i/prIYKAY+BFLc+xYAT3R67fXuz+F24Logxred9v7jjs9gR1XWYODtY30WghTfM+7P1jrak/Kgo+NzP/7Gdz0Y8bm3P93xmeu0b9DPn68/OpReKaXClNW6UJRSSnlIE7hSSoUpTeBKKRWmNIErpVSY0gSulFJhShO4UkqFKU3gSikVpv4/nrz0wwyVofoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Fit a model to some data\n",
    "# Model is:\n",
    "#    prediction = M0 * exp(-R2*TE)*(1-exp(-R1*TR))\n",
    "#    where M0,R1,R2 are unknown parameters and TE/TR are experimental parameters\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "\n",
    "TEs = np.array([10,40,50,60,80]) # TE values in ms\n",
    "TRs = np.array([.8,1,1.5,2])     # TR in seconds (I know this is bad)\n",
    "\n",
    "# All combinations of TEs/TRs\n",
    "comb    = np.array([(x,y) for x in TEs for y in TRs])\n",
    "TEs,TRs = comb[:,0],comb[:,1]\n",
    "\n",
    "# function for our model\n",
    "def forward(p):\n",
    "    M0,R1,R2 = p\n",
    "    return M0*np.exp(-R2*TEs)*(1-np.exp(-R1*TRs))\n",
    "\n",
    "# simulate data using model \n",
    "true_p    = [100,1/.8,1/50]\n",
    "data      = forward(true_p)\n",
    "snr       = 50\n",
    "noise_std = 100/snr\n",
    "noise     = np.random.randn(data.size)*noise_std\n",
    "data      = data + noise\n",
    "\n",
    "plt.plot(data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now for the fitting\n",
    "# we need a cost function:\n",
    "\n",
    " \n",
    "# always a good idea to calculate gradient\n",
    "def forward_deriv(p):\n",
    "    M0,R1,R2 = p\n",
    "    E1,E2    = np.exp(-R1*TRs),np.exp(-R2*TEs)\n",
    "    dE1      = -TRs*E1\n",
    "    dE2      = -TEs*E2\n",
    "    \n",
    "    # f = M0*E2*(1-E1)\n",
    "    dfdM0 = E2*(1-E1)\n",
    "    dfdR1 = M0*E2*(-dE1)\n",
    "    dfdR2 = M0*dE2*(1-E1)\n",
    "    return np.array([dfdM0,dfdR1,dfdR2])\n",
    "\n",
    "def forward_deriv2(p):\n",
    "    M0,R1,R2 = p\n",
    "    E1,E2    = np.exp(-R1*TRs),np.exp(-R2*TEs)\n",
    "    dE1      = -TRs*E1\n",
    "    dE2      = -TEs*E2\n",
    "    ddE1     = (TRs**2)*E1\n",
    "    ddE2     = (TEs**2)*E2\n",
    "    \n",
    "    #dfdM0 = E2*(1-E1)\n",
    "    dfdM0dM0 = np.zeros(E1.shape)\n",
    "    dfdM0dR1 = E2*(-dE1)\n",
    "    dfdM0dR2 = dE2*(1-E1)\n",
    "\n",
    "    #dfdR1 = M0*E2*(-dE1)\n",
    "    dfdR1dM0 = E2*(-dE1)\n",
    "    dfdR1dR1 = M0*E2*(-ddE1)\n",
    "    dfdR1dR2 = M0*(dE2)*(-dE1)\n",
    " \n",
    "    #dfdR2 = M0*dE2*(1-E1)\n",
    "    dfdR2dM0 = dE2*(1-E1)\n",
    "    dfdR2dR1 = M0*dE2*(-dE1)\n",
    "    dfdR2dR2 = M0*ddE2*(1-E1)\n",
    "\n",
    "    return np.array([[dfdM0dM0,dfdM0dR1,dfdM0dR2],\n",
    "                     [dfdR1dM0,dfdR1dR1,dfdR1dR2],\n",
    "                     [dfdR2dM0,dfdR2dR1,dfdR2dR2]])\n",
    "\n",
    "\n",
    "# cost function is mean square error divided by 2\n",
    "def cf(p):\n",
    "    pred = forward(p)\n",
    "    return np.mean((pred-data)**2)/2.0\n",
    "\n",
    "def cf_grad(p):\n",
    "    pred  = forward(p)\n",
    "    deriv = forward_deriv(p)\n",
    "    return np.mean( deriv * (pred-data)[None,:],axis=1)\n",
    "\n",
    "def cf_hess(p):\n",
    "    pred   = forward(p)\n",
    "    deriv  = forward_deriv(p)\n",
    "    deriv2 = forward_deriv2(p)\n",
    "    \n",
    "    H = np.zeros((len(p),len(p)))\n",
    "    for i in range(H.shape[0]):\n",
    "        for j in range(H.shape[1]):\n",
    "            H[i,j] = np.mean(deriv2[i,j]*(pred-data) + deriv[i]*deriv[j])\n",
    "    return H\n",
    "    \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 226,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fitted = [9.90012687e+01 1.28692542e+00 1.97590584e-02]\n",
      "true   = [100, 1.25, 0.02]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xUVfr48c+TEEIoSSgBSYJShFioGgtixRJEvxAVECtr/VpXdxUX/f4su19XcVndta0uVlRUFGlfW0QQFQsaDEWU0ISV0AIYaiDt/P44NxBCQiYzd+beSZ7365XXzNy5c+/DTfJwcu5zzhFjDEoppaJPjNcBKKWUCo4mcKWUilKawJVSKkppAldKqSilCVwppaJUk0ierF27dqZz586RPKVSSkW9+fPnbzbGpFTfHtEE3rlzZ3JzcyN5SqWUinoisqam7dqFopRSUUoTuFJKRSlN4EopFaU0gSulVJTSBK6UUlEqolUojdG0vALG5eSzrqiY1OQERmdlkN0vzeuwlFINgCbwMJqWV8C9UxZTXFoOQEFRMfdOWQygSVwpFTLtQgmjcTn5TvI2xFABQHFpOeNy8r0NTCnVIGgCD6N1RcXEU8KEuMeYGPfIAduVUipU2oUSRp2S4rh/9985I3YRO00zwABCanKC16EppRoAbYGHS0U5b7Z7mXNjf2BBRTdayh4S2U1CXCyjszK8jk4p1QBoAg+Higr4v9+TXvARPx57N5PjLwKgT+IuHr24l97AVEq5QhO424yBnHsh7w0440/0HH4/D189CIDXh6Vq8lZKuUYTuNtmPwzznoeTb4Uz77XbktLt47a13sWllGpwNIG76csn4Mu/w/G/g6y/gojd3rI9xDTRBK6UcpUmcLfMGw+z/gy9hsMFT+xP3gAxsZCYCtsLvItPKdXgaAJ3Q95E+Gg0ZFwA2c/ZhF1dYrq2wJVSrtIEHqolU2HGbdBtIAx/BWLjat4vSRO4UspdmsBDsSwH3rseOp0El06EJvG175uUBtvX2RJDpZRygSbwYK36HCZdBR16wuWToGnzQ++flA4VpbBrU2TiU0o1eJrAg/Hrd/DWZdC2G1w1FZol1f2ZRC0lVEq5q84ELiIZIrKgytd2EblTRNqIyEwRWe48to5EwJ5bvxDeGAatOsBV06B5m8A+p7XgSimX1ZnAjTH5xpi+xpi+wPHAbmAqMAaYZYzpDsxyXjdshfnw+kUQ3wqunm6TeKCSnBGYmsCVUi6pbxfK2cBKY8waYCgwwdk+Ach2MzDf2foLvDYUJBZGzYDkw+v3+WbJ0LSl1oIrpVxT3wQ+EnjLed7BGLPeeb4BqLE5KiI3ikiuiOQWFhYGGabHthXY5F22B66eZvu+60sEEtNg26/ux6eUapQCTuAi0hQYArxb/T1jjMFOdn0QY8x4Y0ymMSYzJSUl6EA9Ne0m2L0VrpwCHY4N/jhJ6fY/A6WUckF9WuDnAz8YYzY6rzeKSEcA57Fh1setngu/fAFn3Qdpx4V2rKQ07QNXSrmmPgn8MvZ3nwDMAEY5z0cB090KylfmjIWWHSDzmtCPlZhu68DL9oZ+LKVUoxdQAheRFsC5wJQqm8cC54rIcuAc53XDsnourP4SBtwJcS4sg1ZZSqg3MpVSLghoTUxjzC6gbbVtW7BVKQ2Xm61vqFJKWABturpzTKVUo6WLGtfGaX0v7jmGmx7/hnVFxaQmJzA6KyP4VXWSOtlH7QdXSrlAE3ht5oxlT3w7rlp4LEWlxQAUFBVz75TFAMEl8cRU+7hdE7hSKnQ6F0pNnNb3v8uHUFR64NzexaXljMvJD+64cQnQvJ2WEiqlXKEJvCZzxkKL9jy38/Qa315XVBz8sbWUUCnlEk3g1a3+ylaenHonbZNrnmUwNTmEipSkTlqFopRyhSbw6j63rW+Ov4bRWRkkxB3YhZIQF8vorIzgj5+oLXCllDv0JmZVq7+yoy6zHoGmzcnuZxdpGJeT704VCtha8L3bYc+2wOYRV0qpWmgCr6pK67tSdr+00BJ2dVVrwTWBK6VCoF0olSpb36feWffyaKGorAXXfnClVIg0gVeqofUdFomVLXCdVlYpFRpN4LC/9T3gjvC2vgFaHWYXhdBacKVUiDSBg9P6ToHMa8N/rphYOyJTK1GUUiHSBL7ma6f1Hea+76qS0rUPXCkVMk3gcyLY+q6kS6sppVzQuBP4mq/hl88j2/oGpwW+DioqIndOpVSD07gTuBetb7AJvLwEdkXpIs9KKV9ovAncq9Y37C8l1GlllVIhaLwJ3KvWN+xfWk0rUZRSIWicCXxf6zsCdd812ZfAtRJFKRW8xpnAvWx9AyS0hrjm2gJXSoWk8SXwNd9UaX238CYGEdsPrn3gSqkQBJTARSRZRCaLyFIR+VlE+otIGxGZKSLLncfW4Q7WFZEcdXkoSenaAldKhSTQFviTwMfGmKOAPsDPwBhgljGmOzDLee1va76BVXO8bX1XSkrTPnClVEjqTOAikgScDrwEYIwpMcYUAUOBCc5uE4DscAXpGr+0vsFOK7tzI5SVeB2JUipKBdIC7wIUAq+ISJ6IvCgiLYAOxpj1zj4bgA41fVhEbhSRXBHJLSz0cOBKZev7lN973/oGpxbcwI51XkeilIpSgSTwJsBxwHPGmH7ALqp1lxhjDGBq+rAxZrwxJtMYk5mSkhJqvMH7fCw0bwcnXOddDFVpKaFSKkSBJPC1wFpjzDzn9WRsQt8oIh0BnMdN4QkxNNPyCrj5kWdg1RyeKbmAaUuKvA7J0sE8SqkQ1bkmpjFmg4j8KiIZxph84GzgJ+drFDDWeZwe1kiDMC2vgIemzOc1eZHNksizO8+AKYsB3F3nMhg6nF4pFaJAFzW+HZgoIk2BVcA12Nb7OyJyHbAGGBGeEIM37uOl3M94esf8wg0lf6SYZlBazricfO8TeNPmkNBGW+BKqaAFlMCNMQuAzBreOtvdcNyVtXMql8R9yROlw5hZsT/8dUXFHkZVRVK69oErpYLWcEdirprDfXETySnP5OnyAyscU5MTPAqqGh3Mo5QKQcNM4L+thnd/x+5WXbmP2zBV/pkJcbGMzsrwLraqktK1D1wpFbSGl8BLdsHbV4CpIPF373D/xSeSlpyAAGnJCTx6cS/v+78rJabBnm2wd4fXkSilolCgNzGjgzEw7RbY9BNc8S607UZ2Wx9UnNSmai14+6O8jUUpFXUaVgt87hPw0zQ45yE48hyvo6mb1oIrpULQcBL4shyY9b/Qc5gdLh8NKhO49oMrpYLQMBL45uXw3vVwWC8Y8rSdbzsatDwMJEZb4EqpoER/At+zDd66DGLjYOREb5ZIC1ZsE2jVUWvBlVJBie6bmBUVMOW/YesqGDUDkg/3OqL6S0qHbb96HYVSKgpFdwt8ziOw7CMYNBY6n+p1NMFJTIPt2gJXStVf9Cbwn6bDF+Og35Vw4g1eRxO8yuH0psbZeJVSqlbRmcA3LoGpN0P6CXDBE9Fz07ImSelQvhd2bfY6EqVUlIm+BL57K7x9OcS3ghGvQ5N4ryMKzb5acO0HV0rVT3Ql8PIymHwtbF8Hl74BiR29jih0++YF135wpVT9+L4KZVpeAeNy8llXVMxfW0zi8vLPYMgz0OkEr0NzR1In+6i14EqpevJ1Ap+WV8C9UxZTXFrO0Ji5XF4+nYkVWbSQgWTX/fHo0LwNNGmmCVwpVW++7kIZl5NPcWk5PWUVj8W9wLcVR/NgyRWMy8n3OjT3iOi84EqpoPg6gduVcwwPxr3GZpK4peQOymjinxV13KK14EqpIPi6CyU1OYGComJuLvkDbWQ7W0nct71BSeoEK2d7HYVSKsr4ugU+OiuDhLhYNpPEMmNv9vlqRR23JKXBjvVQXup1JEqpKOLrBJ7dL41HL+7l3xV13JKUDhibxJVSKkABdaGIyGpgB1AOlBljMkWkDTAJ6AysBkYYY35zO8DsfmkNL2FXV1kLvm1tdE7IpZTyRH1a4GcZY/oaYzKd12OAWcaY7sAs57UKxr5acL2RqZQKXChdKEOBCc7zCdBwSrMjLqmyBa7D6ZVSgQs0gRvgExGZLyI3Ots6GGMqO203AB1q+qCI3CgiuSKSW1hYGGK4DVTTFpDQWksJlVL1EmgZ4anGmAIRaQ/MFJGlVd80xhgRqXE+VGPMeGA8QGZmps6ZWptEHcyjlKqfgFrgxpgC53ETMBU4EdgoIh0BnMdN4QqyUaicF1wppQJUZwIXkRYi0qryOXAe8CMwAxjl7DYKmB6uIBuFpDTtA1dK1UsgXSgdgKliF01oArxpjPlYRL4H3hGR64A1wIjwhdkIJKbBniLYuxPiW3odjVIqCtSZwI0xq4A+NWzfApwdjqAapcpSwu0FkNLARpoqpcLC1yMxG5WkKoN5lFIqAJrA/WLf0mqawJVSgdEE7hetOgKiteBKqYBpAveL2DibxLUFrpQKkCZwP0lK0wSulAqYJnA/0aXVlFL1oAncTyqXVjM644BSqm6awP0kqROU7YHdW7yORCkVBTSB+4nWgiul6kETuJ9oLbhSqh40gftJopPAtRZcKRWAQOcDV5HQoh3ExjeoFvi0vALG5eSzrqiY1OQERmdlNPw1TpWKEE3gfiLSoGrBp+UVcO+UxRSXlgNQUFTMvVMWA2gSV8oF2oXiN0npDaYLZVxO/r7kXam4tJxxOfkeRaRUw6IJ3G8a0NJq64qK67VdKVU/msD9JikddqyH8jKvIwlZanJCvbYrpepHE7jfJKWBqbBJPMqNzsogIS72gG0JcbGMztIFK5RygyZwv2lApYTZ/dJ49OJepCUnIEBacgKPXtxLb2Aq5RKtQvGbBjaYJ7tfmiZspcJEW+B+09CG0/+2Gt65Gt64RCfpUspl2gL3m/hW0Cwp+hN4yS748gn4+mko32u3FS6F9kd7G5dSDUjALXARiRWRPBF533ndRUTmicgKEZkkIk3DF2Yjk+ifWvBpeQUMGDubLmM+YMDY2UzLqyMuY2DRu/B0Jnz5dzhmKFw/27637OPwB6xUI1KfLpQ7gJ+rvH4M+Icx5kjgN+A6NwNr1JLSYduvXkexbyRlQVExhv0jKWtN4usWwMuDYMr10LI9XPsJXPICpB8PHftAviZwpdwUUAIXkXTgAuBF57UAA4HJzi4TgOxwBNgoJaXBNu9b4AGPpNxZCDNuh/FnwtaVMORpuOEzOPyk/fv0GARrv4NdOte5Um4JtAX+T+AeoMJ53RYoMsZUjjZZC9RYaiAiN4pIrojkFhYWhhRso5GUDsVboWS3p2HUOZKyrAS+eRaePh4WvAn9b4Xb58NxV0NMtR+tHoNsffuKmWGOWqnGo84ELiIXApuMMfODOYExZrwxJtMYk5mSkhLMIRofn9SCH3Ik5fJP4blTIOc+6HQC3PwNZP3V3oCtSce+0LKD9oMr5aJAWuADgCEishp4G9t18iSQLCKVVSzpgPd/8zcU+2rBve0Hr2kk5VFxm3gv6UmYeAmYcrhsElwxGVJ6HPpgMTHQIwtWzLItd6VUyOpM4MaYe40x6caYzsBIYLYx5grgM2CYs9soYHrYomxs9tWCe/t/YtWRlC0p5uEW7/Jhk7s57LdcOOfPcMu3kDHIToMbiB6DYO92+M834Q1cqUYilDrwPwFvi8jDQB7wkjshKVqlAuKLWvDsfmlkt/wZpt8HOzdC3yvg7AehVYf6H6zrmXbBimUfQ9cz3A5VqUanXgncGDMHmOM8XwWc6H5IiiZNbX/xdu8TOGu+hrcvh3bdYeRbtiQwWE1bQJfTIf8jyHok8Ja7UqpGOpTer5J8MC/4xiXw1khIPhyunhFa8q6UMQh++wU2Lw/9WEo1cprA/crrWvCi/9j5S+Kaw1VToEVbd47bPcs+ajWKUiHTBO5XSZ1sC9yLCaB2b4XXL7Z16Fe+Z1vgbknuBB16aQJXygU6mZVfJaZBWTEU/wbN2wR9mHqvCl+yCyYOty3wq6ZCh2ODPnetemTB3H/Y/yhC+Lcp1dhpC9yvXJgXvN5zmZSXwru/g3U/wLCXoPOAoM99SBnn2xryFbPCc3ylGglN4H7lwrzg9VoV3hiY8XtY/glc8Dgc/V9Bn7dOqcdBixTtRlEqRJrA/Sqpk30MYTh9vVaFn/VnWPgmnHkvZF4b9DkDEhNjb2aumGlb/UqpoGgC96vm7SC2aUjD6QNeFf7b52yf9PHXwBl/Cvp89dIjC/Zsg1/nReZ8SjVAmsD9KiYGElNDKiUMaFX4xZPh4zFw1IW26yRSg2u6nWX/g8r/KDLnU6oB0gTuZ5WlhEGqc1X4lZ/B1JvgiAFwyUsQE3vI47kqvhV0PhWW5UTunEo1MFpG6GeJabDmq5AOUeuq8OsWwKQroV0PGPkmxDUL6TxB6XE+fDQatqyEtt0if36lopy2wP0sKR22r4OK8rr3rY+tq2DiMEhobQfqJCS7e/xA9Wj4ozKn/bCW3//1n9x1358CW1NUqXrQFrifJaXZeukdG/aXFYZq5yZ4/SL7n8KVUyCxozvHDUbrI6D9MbYfvP+t3sURDhUVfPvxG3SZ9yRPyQpoCmdv68a9U+xc6IccTKVUgLQF7meVpYRuTWq1Z7ud32TnJrji3boXYYiEHll2fvDiIq8jcUd5GSx6B547hZO/u51ks52HS6+g3AgXxc6tvQ5fqSBoAvezRKeV5sa0smV7bZ/3xiUw4jVIzwz9mG7ocT5UlMHKKB+VWbYXcl+BZ46HKTcAcEfJrQwseZwXyy/gy4reZMd+hVBRa32+UvWlCdzPXBhOD0BFha02+eVzGPosdD839Njckp4JzdtGbzXK3p3w9TPwZB94/077bxn5Jtz8NbmJ51COreyZUn4a6bKZk2KW1lqfr1R9aR+4nzVLhPjE0KaVLS6CD/4IS6bYZdD6XuZefG6IiYXu59kbmeVlEBslP5LFv8G88TDvOfu8y+lw0fPQ5Yx9tfSjszK4d8piikvL+aTieHaYBIY3+YrYrMs9Dl41FFHy29J4bWvagYW5eYz64oPAZhOsauVsmH6bvQk68P/BgDvCG2ywemTBwrdg7fdwRH+vozm0HRvh22fh+5egZCdkDIZT/widTjho18rvk50NEr5o0p8h8h1xPXUGRuUOTeA+Ni2vgDbbW9DabDpgNkGoo4qhZBfMfBC+f8HWeV8/E9JcWE0nXLqdDTFNYNlH/k3gv62Br5+CH16HilI49mI49Q9wWM9DfuyAOvxVLeG1IZD/IfS8JAJBq4ZO+8B9bFxOPr+WtyVVtuzbVmcVw6/fw/On2eR98i3w31/4O3mD7So6YoA/+8GNgc/HwVP9YP4E6DMSbsu10+3WkbwP0vk0e2N64aTwxKoaHU3gPrauqJh1pi1tZQfxlByw/SBlJfDpn+Hl8+wMf6Peh0GPQlyU3DDLOB8Kl8LWX7yOZL+KcvjgLvjsYTj2IrhjIQx5KvhRozEx0Gs4rPjUlnIqFaI6E7iINBOR70RkoYgsEZE/O9u7iMg8EVkhIpNEpGn4w21cUpMTWGfsWpRVW+EHVTFs+BFeGAhzn4C+l8PNX0GX0yIZauj2jcr0SSu8bC9MvgZyX+L12IvoknsRA57LD30kZZ+RdnDWj++5E6dq1AJpge8FBhpj+gB9gUEicjLwGPAPY8yRwG/AdeELs3EanZXBltgUADo6CfyA2QQryuHLJ2D8mbBzI1z2ti0TbJboUcQhaNMV2mXYfnCvVQ54+mk6Yyuu5v5dwzFI3SsaBaL90dCxDyx82714VaNVZwI31k7nZZzzZYCBwGRn+wQgOywRNmLZ/dK4Mssua5Ymmw+cTXDLSnh5kF2I4ajBcMu3thsimvXIgtVf2QTqlR0b4dXB8J9v+HPcHTxfMuiAt10ZSdl7JKxfAJuWhnYc1egF1AcuIrEisgDYBMwEVgJFxpgyZ5e1QI1lESJyo4jkikhuYWGhGzE3Kued3A+Acee246sxA8numwrfvQDPnwqb8+HiF2H4BGjR1uNIXZBxvq3wWDnbm/NvXWXvIWxZCZdN4tUdJ9W4W8gjKXsNA4mFRQ2/FT4tr4ABY2fTZcwHOplXGASUwI0x5caYvkA6cCJwVKAnMMaMN8ZkGmMyU1JSggyzEWsSDy072JV5thXYiag+vBsO729b3b2HR24RhnBLPxGaJbvaDx5wAlm/EF46z7b+R/0fdD8n8BWN6qtle+g2EBa9a0fJNlD1XlRb1Vu9qlCMMUXAZ0B/IFlEKuvI0wH9roRLYhqsmgP/6g+/fgcX/sNOA5uY6nVk7optYkdlLs9xZQrdgBPIqs/hlQugSTO4NmffPDEBrWgUrD4j7Rw3a+aGfiyfqtei2ioogVShpIhIsvM8ATgX+BmbyIc5u40CpocryEYvKd22wNsfDTfPtYsON5RWd3U9smD3FiiYH/KhAkogS6baudGTO8F1nxwwQ2OdKxqF4qgLoGmrBl0TXq9FtVVQAhmJ2RGYICKx2IT/jjHmfRH5CXhbRB4G8oCXwhhn43bGPbZ/uPelkV32zAtHnmP7h/M/gk4nhnSoOhPIdy/Ah6Oh00lw+dt2gYtqal3RKFRxCXDMUPhpOgweB02bu38Oj6UmJ1BQw/dAJ/NyTyBVKIuMMf2MMb2NMT2NMX9xtq8yxpxojDnSGDPcGLM3/OE2Uof1svXdDT15g10d6IhTXOkHr7UPO6kZfPaIvZfQYxBcNbXG5B12fS6Fkh12aH0DVNkFJVTQit2Ai11QCtCRmMqPegyCTUug6D8hHaamPuwWccIbHd+Gzx+DflfCpW941/o94lRITPd1TXgoVSTZXWHKsXP5KuGPzIu/lT5Jxe51QSlAE7jyox5O7XWIrfDqfdhdkmKZmf4yXVa/Y2cQHPKMt9PXxsRA7xG2bNKHQ+uDqiIpL4Wf34eJI+CfPTl66dOkHt6d5rKX6Wds0OTtMk3gyn/aHQltj7T94CHK7pfGV2MG8stDp/LZYU+Tun4mDBoL5zzojxvBlUPrF0+ue98Iq1cVyZaV8OlD8I9jYdIVsGERnHaXnT/mmg8h9ThY1HBv2HpFp5NV/tRjEHw33q54E98ytGNtK4A3L7WTZV3ykh1I4xcpGdCxr50Pvf8tXkdzgDpvApfugZ9nwA+vweov7c3nHoPguKvtzeiqf930HgEfj7GjT9sHPIxE1UFb4MqfegyC8hJY9VnwxyjZBXPGwjOZdpTl5ZP8lbwr9RlpW6ybfvY6kgPUdhP4tMSN8OE98HiGXf9z269w9gPwhyVw2ZuQMejgrqmel9gEv/idCETeeGgCV/50+MkQn2SXWquvigpY8CY8fTzMedSuAXrzXDjybPfjdENPZ2i9z25mVr0J3IJiRsbOZnr8A7xW8geY/4ptZV89A27Ps90liR1rP1jL9tDtrAY/+jTStAtF+VNsHHQ/B5Z9Yn/hYwJsa/zyJXzyP3ZofOpxMOwV/67yU6llik2Gi9+Fsx8M/N8aZtn90sBUsPbDcVxT+jYtZC/bWx0JA8baMQnN67k0XK8RMPVG+HWe/78nUUITuPKvHoPsvNnr8iC9jlWFtqyEmQ/A0vdtad7FL9iWrU+SYZ36XAqTc2xfctczvI7G2rWF7J//CGWf2BkvT7uLxLTjg7/5e9QFENfc3szUBO4KTeDKv448ByTGzhFeWwLfvRU+/5tdQq5JMxh4P/S/NXpWIqqUMRjiE203ih8S+JqvYfJ1sHszDP47nHB96FU78S3hqAvt9AXn/w2a6BowoYqS5olqlJq3gU4n19wPXlYC3/zLrlX53b+h7xVw+w9w+t3Rl7zBGVo/xFZ1lOz2Lo6KCvji7/DqhXYmzOtmwok3uFdy2XsE7CmCFTPdOV4jpwlc+VuPLNiwGLatta+NsQNF/nUS5NwLqf3gprl2rcpWHbyNNVS9R0LJTlj6gTfn31kIEy+B2f9r52n57y8gta+75+h6FjRvpzXhLtEErvytcpWhZTm2L/zVC+xAkZg4uGKyncekw7HexuiWIwZAUidvFnr45Uu7SMjqr+DCf8Kwl8OzNF9sE1vKmf8xFBe5f/xGRhO48rd2PaB1Z5j9sF37s3ApXPA43Py1LQ/0w2hKt1QdWr9jQ2TOWVEOcx6D14ZAfCu4YTZkXhPe69prBJTvtd1FKiSawJW/idhBICU7YcAd8Ps8e0PNyzlMwqn3SDAVkRlav2MjvJ4Ncx6BXsPhxjlwWM/wnzftOGjTDRbpoJ5QaQJX/nfmfXDPKjj3L9Asyetowiulh+3XD3c3ysrP4PkB8Ov3dlKvi/4d+pQFgRKxdeSr5+6/t6GCoglc+V9sE/vnfWPRe6S9cbvxJ/ePXV5mu6Nevwiat4UbP4Pjrop8V1Tv4YDx5SRe0UQTuFJ+UzlviNut8O3rbF/3F+Ns2eUNs+0yfV5o0xXST7CjT1XQNIEr5TctU+wN2kXvuLK4MwDLP7VVJuvybHdJ9rPQtIU7xw5W70th44+w4Udv44himsCV8qPel8KO9fDLF6EdZ/0ieO96W9/d8jC48XM7+6EfHHsRxDTRGQpDoAlcKT/KON8OrQ9mwIsxsGqO7ef+92l2YYwBd8ANs+xNUr9o0Q66nW37wXWGwqBoAle+FsqajFFt36r1M+y85oEoL7OTf40/A14bCht+ZMnRfyBLnqfLrJMY8Pg3/rt+vUfA9gJY85XXkUSlOhO4iHQSkc9E5CcRWSIidzjb24jITBFZ7jx6sKy3asiCWpOxIekzEkp32akDDqVkN3z3Ajx9HEy+1ib8/3qKGWflMOzHk8nfFuPf65cxGJq21KH1QQqkBV4G3GWMOQY4GbhVRI4BxgCzjDHdgVnOa6VcU681GRuiw09hd0Iq86b/q+a/QHZvtaMo/9kTPrzbLppw6US49Xs4fhSPfbra/9evaXM4+r/sXxqle7yOJurUOZzNGLMeWO883yEiPwNpwFDgTGe3CcAc4E9hiVI1SnWuydjATVu4ng27TuIGmUY7fqOgCO6dspjmuws4b9tkyHsdSndDj/NtH/fhJx9Qzx0116/3CLsm6PIc222kAlav8cgi0hnoB8wDOjjJHWADEOVTwSm/SU1OoKCGZFPbWo0NzbicfOJLB3BT/FSGxn7NNxXHcsGnxaEAAAtwSURBVCPvM3DmtxATaytVTrm91kWCo+b6dTkDWnawZZOawOsl4JuYItISeA+40xizvep7xhgDmFo+d6OI5IpIbmFhYUjBqsal6pqMlRLiYhmdleFRRJG1rqiYVSaVBRVduavJu3wQfx8DY/J4qWww3LnI1nIfYoX3qLl+MbF2LpZlObZbSAUsoAQuInHY5D3RGDPF2bxRRDo673cENtX0WWPMeGNMpjEmMyUlxY2YVSOR3S+NRy/uRVpyAgKkJSfw6MW97FqNjUBlS/nFsgtYb9owtnQkp+x9mtdaXgeJqXV+PqquX6/hUFEKP033OpKoIrbxfIgdRATbx73VGHNnle3jgC3GmLEiMgZoY4y551DHyszMNLm5uS6ErVTDV1mFU/VGZEJcrH+TcCiMgWdPsvOzXPuR19G4o6LCjjT95XM7IGvYy0HP6SMi840xmdW3B9IHPgC4ClgsIgucbfcBY4F3ROQ6YA0wIqjIlFI1qkzS43LyWVdUTGpyAqOzMhpe8gZnhsIRdjWgov9A8uFeRxScrb/YhL1qjk3au7fY7e0yYFvBIbu8glFnC9xN2gJXStXqtzXwZG84+wE47S6vownMrs37E/aqz6Fojd3eqiN0PdN+dTk9oC6vQwmlBa6UUuHX+gg4vD8snASn/tGfqy3t3Qn/+WZ/wt642G6PT4Iup0H/22zSbtcdRJiWV8C4fy1lXVFeWP6C0gSulPKP3iPg/T/AhkXQsY/X0Vjb1kLeRJu0135vb7bGNrV19wPvtws1d+xz0CpR1e9hVI6EBVxL4prAlVL+cUw2fHiPrQn3OoFvK4C5T8APr0F5qY2n/63Q9QzodLIdRXoIhxpJrAlcKdXwNG8D3c+zMxSe+xdbIx5p29fBl0/ADxPs+qT9rrRdOq2PqNdhIjESVhO4Uspfeo+A/A9sFUe3syJ33u3rbYt7/qs2cfe9wt5MrWfirhSJkbA6naxSyl96DLJzoUdqubXt6+GjP8GTfSD3ZTsL5O3zYchTQSdviMxIWG2BK6X8Ja4ZHDMElkyHwX+vs685aDs2wNx/wvxXbB9338tti7tNF1cOH4k6fk3gSin/6X0p5L0Byz6yizy7accG+OpJ29ouL4W+l8Fpd9eYuKflFYSUgLP7pYV14JUmcKWU/xxxKiSmwaJ3Q0rgVRNwz6Q9/CP9c45cM8km7j6Xwel3QZuutX423GWAodIErpTyn5gYm7i//Rfs2gIt2tb7EDYBLyK17Ffua/IZV+75lLgVZaw5fAhHZD8Abbsd8vORKAMMlSZwpZQ/9b4Uvn4KlkyBE28I/HMV5bD2e3a9/xwfyLd0jd9AuRGmVpzG02XZlBV24as6kjdEx4IYmsCVUv50WE9of6wd1FNXAi8ttkPbl74Pyz6GXYUMN7F8a47h5dLzmVl+PBtpA4AEmICjYUEMTeBKKf/qPRw+fcjO8lf9JuPurXYRiPwPYMUsu7xcfCJ0PxcyBnPh+/Es23ZwpXSgCXh0VkaN0/n6aUEMTeBKKf/q5STwxe/CGffAb6th6YeQ/yGs+RpMObRKtSWAGYOh82nQpCkAt5TVPJ96oAk4Gqbz1elklVL+9uqFsHkZtGi/f/a/9sfYhH3UYOjYz970rEGoZYB+odPJKqWi0/G/gyk3QJtucN5fbdKupfSvunDXYXtNE7hSyt96DbOr1cfGeR2J7+hcKEop/9PkXSNN4EopFaU0gSulVJTSBK6UUlFKE7hSSkWpOhO4iLwsIptE5Mcq29qIyEwRWe48tg5vmEoppaoLpAX+KjCo2rYxwCxjTHdglvNaKaVcNS2vgAFjZ9NlzAcMGDubaXkFXofkK3UmcGPMF8DWapuHAhOc5xOAbJfjUko1cpXzcRcUFWPYPx+3JvH9gu0D72CMWe883wB0qG1HEblRRHJFJLewsDDI0ymlGptDzcetrJBvYho7mUqtE6oYY8YbYzKNMZkpKSmhnk4p1UhEw3zcXgs2gW8UkY4AzuMm90JSSqnap33103zcXgs2gc8ARjnPRwHT3QlHKaWs0VkZJMTFHrDNb/Nxe63OyaxE5C3gTKCdiKwFHgTGAu+IyHXAGmBEOINUSjU+0TAft9d0PnCllPK52uYD15GYSikVpTSBK6VUlNIErpRSUUoTuFJKRSlN4EopFaUiWoUiIoXYssNgtAM2uxiO2zS+0Gh8odH4QuP3+I4wxhw0lD2iCTwUIpJbUxmNX2h8odH4QqPxhcbv8dVGu1CUUipKaQJXSqkoFU0JfLzXAdRB4wuNxhcajS80fo+vRlHTB66UUupA0dQCV0opVYUmcKWUilK+S+AiMkhE8kVkhYgctFiyiMSLyCTn/Xki0jmCsXUSkc9E5CcRWSIid9Swz5kisk1EFjhfD0QqPuf8q0VksXPug6Z+FOsp5/otEpHjIhhbRpXrskBEtovIndX2iej1E5GXRWSTiPxYZVsbEZkpIsudx9a1fHaUs89yERlV0z5him+ciCx1vn9TRSS5ls8e8mchjPE9JCIFVb6Hg2v57CF/18MY36Qqsa0WkQW1fDbs1y9kxhjffAGxwEqgK9AUWAgcU22fW4DnnecjgUkRjK8jcJzzvBWwrIb4zgTe9/AargbaHeL9wcBHgAAnA/M8/F5vwA5Q8Oz6AacDxwE/Vtn2N2CM83wM8FgNn2sDrHIeWzvPW0covvOAJs7zx2qKL5CfhTDG9xBwdwDf/0P+rocrvmrvPw484NX1C/XLby3wE4EVxphVxpgS4G1gaLV9hgITnOeTgbNFRCIRnDFmvTHmB+f5DuBnINpmlx8KvGasb4HkyuXxIuxsYKUxJtiRua4wxnwBbK22uerP2AQgu4aPZgEzjTFbjTG/ATOBQZGIzxjziTGmzHn5LZDu9nkDVcv1C0Qgv+shO1R8Tt4YAbzl9nkjxW8JPA34tcrrtRycIPft4/wQbwPaRiS6Kpyum37AvBre7i8iC0XkIxE5NqKB2QWmPxGR+SJyYw3vB3KNI2Ektf/ieHn9ADoYY9Y7zzcAHWrYxy/X8VrsX1Q1qetnIZxuc7p4Xq6lC8oP1+80YKMxZnkt73t5/QLitwQeFUSkJfAecKcxZnu1t3/Adgv0AZ4GpkU4vFONMccB5wO3isjpET5/nUSkKTAEeLeGt72+fgcw9m9pX9baisj/AGXAxFp28epn4TmgG9AXWI/tpvCjyzh069v3v0t+S+AFQKcqr9OdbTXuIyJNgCRgS0Sis+eMwybvicaYKdXfN8ZsN8bsdJ5/CMSJSLtIxWeMKXAeNwFTsX+qVhXINQ6384EfjDEbq7/h9fVzbKzsVnIeN9Wwj6fXUUR+B1wIXOH8J3OQAH4WwsIYs9EYU26MqQBeqOW8Xl+/JsDFwKTa9vHq+tWH3xL490B3EenitNJGAjOq7TMDqLzjPwyYXdsPsNucPrOXgJ+NMU/Uss9hlX3yInIi9hpH5D8YEWkhIq0qn2Nvdv1YbbcZwNVONcrJwLYq3QWRUmvLx8vrV0XVn7FRwPQa9skBzhOR1k4XwXnOtrATkUHAPcAQY8zuWvYJ5GchXPFVvadyUS3nDeR3PZzOAZYaY9bW9KaX169evL6LWv0LWyWxDHuH+n+cbX/B/rACNMP+6b0C+A7oGsHYTsX+Ob0IWOB8DQZuAm5y9rkNWIK9q/4tcEoE4+vqnHehE0Pl9asanwDPOtd3MZAZ4e9vC2xCTqqyzbPrh/2PZD1Qiu2HvQ57T2UWsBz4FGjj7JsJvFjls9c6P4crgGsiGN8KbP9x5c9gZVVWKvDhoX4WIhTf687P1iJsUu5YPT7n9UG/65GIz9n+auXPXJV9I379Qv3SofRKKRWl/NaFopRSKkCawJVSKkppAldKqSilCVwppaKUJnCllIpSmsCVUipKaQJXSqko9f8BXjmgSwZu2twAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# get ready to minimize\n",
    "p0 = [200,1/1,1/70] # some random guess\n",
    "method = 'trust-ncg'\n",
    "\n",
    "kw_args = {'x0':p0,'method':method,'jac':cf_grad,'hess':cf_hess}\n",
    "\n",
    "result = minimize(cf,**kw_args)\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(data,'o')\n",
    "plt.plot(forward(result.x))\n",
    "print('fitted = {}'.format(result.x))\n",
    "print('true   = {}'.format(true_p))\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}