fit_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Short example of Model fitting\n",
    "\n",
    "Here we fit a MRI-style model to some simulated data.\n",
    "\n",
    "The model is: $\\textrm{Signal} = M_0\\exp\\left[-R_2\\textrm{TE}\\right]\\left(1-\\exp\\left[-R_1\\textrm{TI}\\right]\\right)$ \n",
    "\n",
    "The parameters that we will be fitting are $(M_0,R_1,R_2)$. \n",
    "\n",
    "Basic imports:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import minimize\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section:\n",
    "\n",
    "- defining a numpy array\n",
    "- double list comprehension\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now we define our forward model\n",
    "\n",
    "In this section:\n",
    "\n",
    "- inline function definition\n",
    "- random number generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x121b13978>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXyU9bX48c/JvpBMSDJZIOxZQBAQIrILglZcilqt0KrQ2lJfta12udbWe/3dW9tbl7Zqba+W1gWXuqFWcUGpIrLIEpawCCRhCUsSkhASspB1vr8/ZoIxJGSS2ZPzfr3yyswzz8xzXk8mJ0/OfL/fI8YYlFJKBZ4gXweglFKqZzSBK6VUgNIErpRSAUoTuFJKBShN4EopFaBCvHmwxMREM3ToUG8eUimlAt7WrVvLjTHW9tu9msCHDh1KTk6ONw+plFIBT0QKO9quJRSllApQmsCVUipAaQJXSqkApQlcKaUClCZwpZQKUJrAlVIqQGkCV0qpAKUJ3At2HqtkTV6Zr8NQSvUymsA9bPfxKhYu3civ39zl61CUUr2MV2di9jVHTtax+Nkt1Da2UN9cT4vNEBwkvg5LKdVL6BW4h5ysaWDRs5tparGxeOpQWmyGsuoGX4ellOpFNIF7QF1jM99dlkNR5RmeWZzNjIxEAIqrzvg4MqVUb6IJ3M2aWmzc+dI2dh2r5ImFFzFxSDwplggASqrqfRydUqo30Rq4GxljuO+tXazeX8bvrh/DFaNTAEi1RAJQrAlcKeVGegXuRn9alcdrOcf4yWXpfPuSIWe3948KJTwkSEsoSim30gTuJi9uLOSJTwq4OXsQP7088yuPiQiplgi9AldKuZUmcDdYubuE+9/ezZyRSfzu+jGInDtUMMUSoTVwpZRbaQJ30ZbDFfzkle2MTYvjiW9dREhwx6c01RKpV+BKKbfSBO6C/BPV3P7cFtLiInlm8cVEhXX+mXCKJYITp+ux2YwXI1RK9WaawHuouOoMi57ZTHhoMMu+O4n46LDz7j/AEkGzzVBeo5N5lFLuoQm8B6rONLH4mS2crm/m2cUXMyg+qsvnpOhQQqWUm2kC76b6pha+/3wOB8tr+NutExkz0OLU81Idk3k0gSul3EUn8nRDi83ws9d2sPlQBY8vGM+09ESnn/vlbEwdC66Uco8ur8BFJEtEdrT5Oi0id4tIvIisEpF8x/f+3gjYV4wxPPDuF7y/q4T7rhrF/PEDu/X8hOgwwoKD9ApcKeU2XSZwY8x+Y8x4Y8x4YCJQB7wF3At8bIzJAD523O+1nlpzkOc2HOb26cP4/szh3X6+iJCik3mUUm7U3Rr4HOCAMaYQmA8sc2xfBlznzsD8yTu5RTy0ch/XjhvAfVeN6vHr6GQepZQ7dTeBLwBedtxONsYUO26XAMkdPUFElohIjojklJUFXlsxYwyP/zuPMQNj+cNNYwlyoSFDqiWC4tNaA1dKuYfTCVxEwoCvA6+3f8wYY4AOZ6gYY5YaY7KNMdlWq7XHgfrK3uJqDpTVsuDiwYSHBLv0WimWCE5UNehkHqWUW3TnCnwesM0Yc8Jx/4SIpAI4vpe6Ozh/sGJnEcFBwrwxKS6/1gBLJI0tNk7WNrohMqVUX9edBL6QL8snAO8Aixy3FwFvuysof2GMYUVuEdPTE0noF+7y62ljB6WUOzmVwEUkGrgceLPN5geBy0UkH5jruN+r7DhaybFTZ7h23AC3vN6Xk3m0Dq6Ucp1TE3mMMbVAQrttJ7GPSum1VuQWExYcxBWjO/x8ttvOXoGf1itwpZTrdCp9J1pshnd3FjEry0psRKhbXjMxOpzQYKGoUhO4Usp1msA7seVwBaXVDW4rnwAEBQnJsRE6nV4p5RaawDuxIreIyNBg5oxKcuvrams1pZS7aALvQFOLjQ92lzD3guTzNmnoiRRLpNbAlVJuoQm8A+sLyqmobeTasaluf+0Bjitw+9wnpZTqOU3gHViRW0xMRAiXZrl/5miKJYLGZhsVOplHKeUiTeDt1De18NGeEr42OsXlqfMd0cYOSil30QTezpq8Mqobmt06+qSt1tZqOhtTKeUqTeDtrMgtIj46jKkjErreuQfOXoHrB5lKKRdpAm+jrrGZj/eWMm9MCqHBnjk1if3CCQkSiit1LLhSyjWawNv4995SzjS1eKx8AhB8djKPXoErpVyjCbyNFblFJMeGc/HQeI8eR1urKaXcQRO4Q9WZJtbsL+PqCwcQ7ELXHWekWCJ0Mo9SymWawB0+2lNCY4uNa8e5f/JOe/bJPGd0Mo9SyiWawB1W7CwmrX8k4wfFefxYKZZI6ptsVNY1efxYSqneSxM4cLKmgfUF5Vw7bgAini2fgE7mUUq5hyZw4IPdJbTYDNeO9dzok7a+bOygQwmVUj2nCRz76JMR1mhGpcZ45Xh6Ba6Ucoc+n8BLqurZfLjCa+UTgKSYCIKDRMeCK6Vc0ucT+Hu7ijEGrvFS+QTsk3mSYsK1tZpSyiV9PoGvyC3igtRY0pP6efW49rHgWgNXSvVcn07gRyvq2HG00qNT5zujrdWUUq5yKoGLSJyILBeRfSKyV0SmiEi8iKwSkXzH9/6eDtbdVuwsAuAaD3Te6UqqJZIS7cyjlHKBs1fgjwMrjTEjgXHAXuBe4GNjTAbwseN+QFmRW8xFg+MYFB/l9WOnWiKoa2zh9Jlmrx9bKdU7dJnARcQCzASeBjDGNBpjKoH5wDLHbsuA6zwVpCcUlFazt/i018Z+t5dydl1wrYMrpXrGmSvwYUAZ8KyIbBeRf4hINJBsjCl27FMCJHf0ZBFZIiI5IpJTVlbmnqjdYEVuMSJwtQ/KJ6BjwZVSrnMmgYcAE4AnjTEXAbW0K5cYeyG3w2KuMWapMSbbGJNttbq/SXBPGGNYsbOIS4bFkxwb4ZMYUrW1mlLKRc4k8GPAMWPMJsf95dgT+gkRSQVwfC/1TIju90XxaQ6W1fpk9Ekra0w4QaJX4EqpnusygRtjSoCjIpLl2DQH+AJ4B1jk2LYIeNsjEXrAitxiQoKEeWN8Uz4BCA0OwhoTrq3VlFI9FuLkfj8GXhKRMOAg8B3syf81EbkdKAS+6ZkQ3csYw4rcIqZnJBIfHebTWFIskdrYQSnVY04lcGPMDiC7g4fmuDccz9t+tJLjlWf42eWZvg6F1NgICspqfB2GUipA9bmZmCtyiwgLCeLy0R0OmvGq1DhtbqyU6rk+lcBbbIb3dhYzO8tKbESor8Mh1RJBTUMzp+u1M49Sqvv6VALfdOgkpdUNPh190laKDiVUSrmgTyXwFbnFRIUFc9nIJF+HAuhkHqWUa/pMAm9qsfHB7mLmjkomKszZwTee1ZrAS6p0KKFSqvv6TAJfV1BOZV2T35RPwN6ZR3Qyj1Kqh/pMAl+RW0RMRAgzMxN9HcpZYSFBJPYLp1g78yileqBPJPD6phY+2nOCK0enEB4S7OtwviLVEkGxTuZRSvVAn0jgn+4vo6ah2a/KJ61SLRFaA1dK9Yh/fJrnIVV1Tby+9SjPrDtEfHQYU0ck+Dqkc6RaItlw4KSvw1BKBaBemcC/KDrN858f5l87jlPfZCN7SH9+fkUWIcH+9w9HiiWC6vpmahqa6RfeK38cSikP6TUZo7HZxso9JTy/4TA5haeICA3iuvEDuXXKEEYPsPg6vE61HUqYnhTj42iUUoEk4BN4SVU9/9x8hJc3H6GsuoEhCVH859WjuGniICxRvp8u35WU2C8n82gCV0p1R0AmcGMMmw9V8PznhazcU4LNGGZnJXHrlCFcmmElKEh8HaLTBsTZp9PrWHClVHcFVAKvbWjmre3HeeHzQvafqMYSGcrt04dxyyVDGJzg/c7y7pAUGw7oeihKqe4LiAR+sKyGFzYWsjznGNUNzYweEMvD3xjLteMGEBnmX+O6uys8JJjEfmEU61BCpVQ3BUQC//Vbu9haeIqrLkzltilDmTA4DpHAKZN0JcUSoSUUpVS3BUQCf2D+GOKiwrDGhPs6FI9ItURytKLO12EopQKM/w2M7kBGckyvTd7gmE6vV+BKqW4KiATe26VYIqg600RdY7OvQ1FKBRBN4H5AGzsopXrCqQQuIodFZJeI7BCRHMe2eBFZJSL5ju/9PRtq75USq63VlFLd150r8NnGmPHGmGzH/XuBj40xGcDHjvuqBwbE6RW4Uqr7XCmhzAeWOW4vA65zPZy+KTlWW6sppbrP2QRugI9EZKuILHFsSzbGFDtulwDJHT1RRJaISI6I5JSVlbkYbu8UERpMfHQYRXoFrpTqBmfHgU83xhwXkSRglYjsa/ugMcaIiOnoicaYpcBSgOzs7A73UfZFrbQGrpTqDqeuwI0xxx3fS4G3gEnACRFJBXB8L/VUkH3BgDgdC66U6p4uE7iIRItITOtt4ApgN/AOsMix2yLgbU8F2RekaGs1pVQ3OVNCSQbecqw9EgL80xizUkS2AK+JyO1AIfBNz4XZ+6VaIjlV18SZxpaAX6BLKeUdXSZwY8xBYFwH208CczwRVF/U2tih5HQ9wxKjfRyNUioQ6ExMP5F6diy4llGUUs7RBO4nUi06G1Mp1T2awP1E296YSinlDE3gfiIyLJi4qFAtoSilnKYJ3I/oZB6lVHdoAvcjA+IitYSilHKaJnA/Yp/MowlcKeUcTeB+JDU2gpO1jdQ3tfg6FKVUANAE7kdSHJ15TpzWq3ClVNc0gfuRAXH2seBaB1dKOUMTuB9pvQLXOrhSyhmawP2ITuZRSnWHJnA/Eh0eQmxEiE7mUUo5RRO4n0m16FhwpZRzNIH7mdQ4HQuulHKOJnA/k2rR1mpKKedoAvczKbGRlNc00NCsk3mUUuenCdzPpDqGEpaebvBxJO7z7s4i/rq6wNdhKNXrONMTU3nRl5156hkUH+XjaFxjjOGpNQd5aOU+AG6amEaSY6ikUsp1egXuZ1qvwAN9KKHNZvjNu1/w0Mp9TB4eD8CneWU+jkqp3kUTuJ9J6QWt1Rqbbdz96g6eXX+Y704bxkvfm0xybDif7i/1dWhK9SpOJ3ARCRaR7SLyruP+MBHZJCIFIvKqiIR5Lsy+o194CDHhIQE7EqWmoZnvPreFd3KLuHfeSP7rmlEEBwmzs5JYm1dOU4vN1yEq1Wt05wr8LmBvm/sPAY8aY9KBU8Dt7gysL0uNiwjIEkp5TQMLl27k84Mn+cNN47jj0hGICACzspKobmhma+EpH0epVO/hVAIXkTTgauAfjvsCXAYsd+yyDLjOEwH2RSmWyIAroRw5WceNT24gv7Sav982kRsnpn3l8WnpCYQGC6u1jKKU2zh7Bf4YcA/Q+v9vAlBpjGl23D8GDOzoiSKyRERyRCSnrEw/xHJGamxgTebZfbyKG57cQOWZJv75/clcNjL5nH1iIkK5eGg8n+7T94BS7tJlAheRa4BSY8zWnhzAGLPUGJNtjMm2Wq09eYk+J8USQVlNA43N/l8v3lBQzoKlGwkLFpbfMYUJg/t3uu/srCT2n6jmeGXglYeU8kfOXIFPA74uIoeBV7CXTh4H4kSkdRx5GnDcIxH2QamWCIyB0mr/vgp/d2cRi5/dwoC4CN744VTSk2LOu//skfY/4DoaRSn36DKBG2N+ZYxJM8YMBRYAnxhjvg2sBm507LYIeNtjUfYxqXH+P5Rw2YbD/Pjl7YxNs/D6D6aS6hj+eD4jrP1I6x/Jai2jKOUWrowD/yXwMxEpwF4Tf9o9IakvJ/P4XwI3xvCHD/fz/97Zw5yRybz4vUuwRIU69VwR+3DC9QXlfWqtl3X55Ty0ch/GGF+HonqZbiVwY8ynxphrHLcPGmMmGWPSjTE3GWN6z+IdPuavrdWaW2zc+8Yu/rK6gAUXD+KpWyYQERrcrdeYPdLKmaYWNh+q8FCU/mXnsUq+/3wOT356gPzSGl+Ho3oZnYnph2LCQ4gOC6bIj8aCn2ls4Y4Xt/FqzlF+fFk6v7/hQkKCu//2mTI8kbCQoD5RRjl2qo7bl+UQG2n/qOgzXUpAuZkmcD8kIqTG+c9Y8IraRm59ehMf7zvBb+aP5udXZJ2doNNdkWHBTBme0Os/yDxd38Ttz+VQ39TCi7dfQnpSP9ZoAldupgncT/lDY4fdx6u4942dTHvwE3Yeq+IvCydw25ShLr/u7CwrB8trOVxe63qQfqipxcadL23jQFkNT90ykYzkGGZmWNl8qIL6pr5T+1eepwncT6XE+qa1Wn1TC29uO8b1/7eea55Yx792HOfr4wbwzo+ncfXYVLccY1ZWEtA7hxMaY7j/7d2szS/nf6+/kGnpiQDMzEykodnGxoMnfRyh6k10PXA/lWqJoLS6nuYWW49qzd115GQdL20u5LUtRzlV18TwxGjuv+YCvjEhzelRJs4amhjN8MRoVu8vY/G0YW59bV/722cHeXnzUe6cPYJvXjzo7PbJwxMIDwnis7zys3/AlHKVJnA/lWKJxGagtLqBAXFdj7HuiRab4dP9pbywsZA1eWUEiXDFBcncMnkIU0ck9LjO7YxZWUm8uKmQM40tRIZ1bySLv3pvZzEPfrCPa8cN4OeXZ33lsYjQYCYNi+ezfK2DK/fRBO6n2nbmcXcCL69p4LWco7y08QjHK8+QFBPOTy7LYOGkwWeHMHra7JFWnll/iM8Plne4dkqg2Vp4ip++toPsIf155MaxBAWd+8fv0kwrv31vL8crzzDQQ3+UVd+iCdxPpbp5LLgxhq2Fp3hhYyHv7yqmqcUwZXgC9109issvSCbUC2WatiYNiycyNJjV+8oCPoEfOVnHkudzSLVEsPS27E7Hxs/MtMJ7e/ksr4yFkwZ7OUrVG2kC91OpsfYrNHesC77hQDm/WfEF+0qqiQkP4duXDOGWyYO7XLvEk8JDgpmWnsjq/aUYYzxarvGkqromvvPcZppthmcXX0x8dOd9TTKS+pFqidAErtxGE7ifio0MITI02OWhhEcr6vjBC1uJjw7j9zdcyPzxA4gK848f++yRVv699wQHymp8+sekpxqbbfzgxRyOVpzhhdsnMdza77z7iwgzM6y8v7vYax9Oq95N30F+yj6Zx7WhhI3NNn70z20AvHj7JSycNNhvkjd8OZwwEGdlGmO4982dbDxYwcM3juWS4QlOPW9mppXq+mZ2HK30cISqL9AE7sfsk3l6XkL5/Qd7yT1WxSM3jmNQfJQbI3OPgXGRZCXHBGSXnic+KeDNbcf56dxMrruow14mHZqenkiQ6LR65R6awP1YSmzPp9Ov3F3Cs+sP851pQ7lyTIqbI3OfWSOtbDlcQXV9k69Dcdq/th/nT6vyuGHCQH4yJ71bz7VEhTJ+UBxr8ss9FJ3qSzSB+7FUSwQnqhtosXVvGdKjFXXcszyXsWkWfjVvlIeic4/ZWUk0tRjWFwTGDMXNhyq4Z/lOJg+P58Ebxvbow9eZmVZ2HqvkVG2jByJUfYkmcD+WGhdBi81QVu38Sr2tdW8D/GXhBMJC/PtHPHFIf2LCQwJiWv3BshqWvJBDWnwkf7slu8fndmamFWNgbYFehSvX+Pdvdx/3ZWMH5+vgD36wz1H3HsvgBP+re7cXGhzEjMwvhxP6q4raRr773BaCRXhu8SSXlhcYlxaHJTJU6+DKZZrA/VhKbPdaq324p4Rn1h9i8dShXDnGPQtPecOsrCROnG5gb3G1r0PpUH1TC0uez6Goqp6lt2W7/IcxOEiYnpHI2vwyv/6jpfyfJnA/1p3Wakcr6viP1x1176tGejo0t5qVaW927K+jUf740X5yCk/x6DfHM3FIf7e85qUZVk6cbmD/Cf/8o6UCgyZwPxYXFUp4SFCXJZTGZhs/enn72bp3eEhgLQ6VFBvBmIGxflkHP155hmUbCrlpYprbltMFmJFpX2Z2zX4to6ie0wTux0SEAXGRXV6BP7RyH7lHKwOm7t2R2VlJbC08RVWdfw0nfGxVHgB3X57p1tdNtdjHwOvqhMoVmsD9XFeNHT7aU8LT6wKv7t3erKwkbAa/SmgFpdW8se0Yt0we4pHVA2dmJrLl0CnqGpvd/tqqb+gygYtIhIhsFpFcEdkjIv/j2D5MRDaJSIGIvCoina/io3rsfK3VjlbU8YvXc7lwYODVvdsbPyiOuKhQv6qD/+HDPCJDg7lz9giPvP7MTCuNLTY2HazwyOur3s+ZK/AG4DJjzDhgPHCliEwGHgIeNcakA6eA2z0XZt+VYongxOn6cybznK17G/jrtwKv7t1ecJBwaaaVNfvLsHVz4pIn5B6tZOWeEr43YzgJ/cI9coyLh8YTERqkzY5Vj3WZwI1djeNuqOPLAJcByx3blwHXeSTCPi41LpJmm+FkzVcn87TWvR8O4Lp3e7OzkjhZ28iu41W+DoWHP9xHfHQY35vhuZZvEaHBTB6eoOPBVY85VQMXkWAR2QGUAquAA0ClMaa1eHcMcH5FH+W01NhzhxK21r0XTRnCvAsDt+7d3sxMKyK+H064Lr+c9QUnuXN2OjER7u0H2t7MDCsHy2s5WlHn0eOo3smpBG6MaTHGjAfSgEmA0wVXEVkiIjkiklNWplca3ZXSbix4a917zMBYfn21f69z0l3x0WGMHxTHah8OrTPG8MiH+xhgieDbl3i+6cJMxxh4LaOonujWKBRjTCWwGpgCxIlI6+LSacDxTp6z1BiTbYzJtlqtLgXbF33ZWu0Mjc02ftyL6t4dmZ2VxM5jlZTXOL/+izt9uKeE3GNV3H15Zqet0dxphDWagXGRWkZRPeLMKBSriMQ5bkcClwN7sSfyGx27LQLe9lSQfVl8dBhhIUEUV9Xz8Mp97DhayUM3jmVIQrSvQ/OI2VlJGOOb9bKbW2w88uF+RlijuaEba3y7QkSYmWllw4GTNLXYvHJM1Xs4cwWeCqwWkZ3AFmCVMeZd4JfAz0SkAEgAnvZcmH2XiJBqieC9XcX8w1H3vqoX1b3bGz0glsR+4T4po7y5/TgHymr5j69lebXd2aWZidQ0NLP9iHbpUd3TZX8tY8xO4KIOth/EXg9XHpYSG8GmQxW9su7dXlCQMCvLyqovTni1b2R9UwuPrcpjXJqFr432bgOMqemJBAcJa/JKmTQs3qvHVoFNZ2IGgKEJ0cSEh/Taund7s7OSqDrT5NW+kS9tOkJRVT33XDmyR00aXBEbEcqEwXF8lqfrg6vu0QQeAH599ShW/nRmr617tzc9w35F6q3hhDUNzfx1dQHT0hOYlp7olWO2NzPDyu6iqnPG+yt1PprAA4AlMtQja3H4K0tkKBOH9Pdat/p/rD1IRW0j93zNd8sRtHbpWaddelQ3aAJXfml2VhJfFJ/ucVNnZ52saeDvnx1k3pgUxg2K8+ixzmfMQAv9o0J1eVnVLZrAlV+aPbJ1gotnyyj/9+kBzjS18PMr3LtcbHcFBwkzMqx8ll/uF2vBqMCgCVz5pazkGFItER4toxyvPMMLnxdy48Q00pNiPHYcZ83MtFJe08DektO+DkUFCE3gyi+JCLOyklhXUE5js2cmuDz+b3uzhrvm+vbqu9XMDEeXHp2VqZykCVz5rdlZVmoamskpdP962QWl1Szfeoxbp3imWUNPJMVGMCo1ttdMqz/T2MKuY1W8sfUYr2w+og2cPaDLiTxK+cq09ERCg4VP95cxdYR7h/f98SN7s4YfzvJMs4aempmZyDPrDlHb0Ex0eGD8etY1NnOgtJa8E9XklVZTcKKGvNJqjp06Q9uc3dRi49YpQ30WZ28UGO8Q1SdFh4dwybAEVu8r5ddXuW8Gau7RSj7YXcLdczM81qyhpy7NsPK3NQf5/MBJ5l6Q7OtwvsKZRB0aLAxP7Me4tDhunDCIzOR+ZCT344F39/Lb9/YyZUSCX3ze0FtoAld+bVaWld++t5ejFXUMindP44pHPtzvaNYw3C2v504Th/YnMjSYNXllfpHA6xqbeXrtIV7feoyjp+rOJuqw4CCGW6MZlxbHTRMHkZHUj4zkGIYkRBHawfIHj9w0lisfW8tdr+zgrR9OIyxEq7fuoAlc+bXZI5P47Xt7+TSvjFsnD3H59dYXlLOuoJz/uuYC+vlhiSI8JJipIxJ83ty5ucXG8q3H+NOqPEqrG7g008qNE9PITO5HelIMQxOiurVOTVJMBA99Yyzffz6HP67az6/m9e41fbzF/97BSrUxPDGawfFRfLSnhAUXD+rw6s5ZxhgeXrmPgXGRXmnW0FMzM618vK+UwpO1Xl8+wRjDJ/tKefCDfeSX1jBxSH+evGUCE4e4vsjW5Rcks3DSYJZ+dpBZmUlMGZHghoj7Nv0/Rvk1EWHehSmszS/not+s4vvP5/DCxkIKT9Z2+7VamzXcNTfDK80aeqq1S4+3R6PkHq1kwdKN3L4shxab4W+3TmT5HVPckrxb/dc1oxiaEM3PX9tBVV2T2163r9IrcOX3fnFFFhcNimNNXjmf5ZWx6osTAAxJiGJGRiIzM6xMGZFw3v6Vrc0a0pP6ea1ZQ08NTYhiUHwka/LKvDJqo/BkLY98uJ93dxaT2C+MB64b4/J/O52JCgvhsZvH840nN/Cfb+/mzwvGe331x95EE7jye6HBQVw5JpUrx6RijOFQeS1r8+3J/M1tx3lx4xFCgoQJg/vbE3qmlTEDLQQHfZkYWps1PHXLBK82a+gJEeHSTCtvbTtOY7PNYx/4VdQ28sQn+by4sZCQoCB+MieDJTOHe/yzgXGD4rh7bgZ/+CiPy0Zauf6iNI8erzcTbw6uz87ONjk5OV47nur9GpttbC08xdr8Mj7LL2P3cfs09LioUKan26/OLxkez8KlG7HGhPOvO6cFxBXfR3tKWPLCVl7+/mS314rrm1p4Zv0hnlx9gNrGZm6+eDA/nZtBUmyEW49zPi02w4Kln7OvuJr375rhthFGvZWIbDXGZJ+zXRO46k1O1jSwrqCcz/LKWZtfRmn1l+trv/S9S3y23nd3Vdc32Wv+M4fzyyvds8xti83w5jb7yJLiqnrmjkrml1dmkZHsm3HZRyvqmPf4Wi5IjeXlJZO/8h+T+qrOEriWUFSvktAvnPnjBzJ//ECMMew/Uc3avHJsxgRM8gaIibCvib5mf5nLCdwYw5q8Mh78YB/7SqoZNyiOx24ezyXDfTsKZMEi16AAAA2sSURBVFB8FL+ZP5qfvZbLU2sOcOfsdJ/GE4g0gateS0QYmRLLyJRYX4fSIzMzrTzy4X7KqhuwxnR/xmhDcwvv7yrmufWHyT1WxZCEKP76rQlcdWGK35SRrr9oIJ/sK+XRVXnMyEhkbJrv1mQPRP79aY5SfdiljuGEa7s5qaf0dD1/WpXHtAdX89NXc6lpaOaB68aw6qeXcvXYVL9J3mD/I/u76y4kKSacu1/ZQV1js69DCih6Ba6Un7ogNZaE6DDW5JVxw4SuR2psP3KK5zYc5r2dxbQYw2VZSSyaOpTp6YkE+XF92RIVyh++OY5v/2MTv31vL/97/YW+DilgdJnARWQQ8DyQDBhgqTHmcRGJB14FhgKHgW8aY055LlSl+pagIGFmppU1eWXYbKbDJNy+TBITHsJtU4Zy25QhDE0MnCbYU0cksmTGcP722UFmZyVxuR+sAxMInLkCbwZ+bozZJiIxwFYRWQUsBj42xjwoIvcC9wK/9FyoSvU9MzMTeWv7cfYUnebCNMvZ7aWn63lp0xFe2nSE8poGhlujeWD+aK6fkOaXa7w442dXZLI2v5xfvrGTcYNmkBTjvWGNgarLn7QxphgodtyuFpG9wEBgPjDLsdsy4FM0gSvlVjMyHNPq88u4MM1ytkzy/q5imm2G2VlJLA6AMokzwkOC+fPC8Vz953Xcs3wnzy6+2K/q9f6oW3+qRWQocBGwCUh2JHeAEuwllo6eswRYAjB4sP8uIKSUP0rsF87oAbEs33qMj744Qe7RSmLCQ7h1cuCVSZyRnhTDfVeP4v639/DCxkJu0wYQ5+X0KBQR6Qe8AdxtjPlK11Vjnw3U4YwgY8xSY0y2MSbbarW6FKxSfdGckUkcKq+lur6JB+aP5vNfz+H+ay/odcm71a2ThzAry8rv3ttL/olqX4fj15yaiSkiocC7wIfGmD85tu0HZhljikUkFfjUGJN1vtfRmZhKdd+ZxhbyS6sZM8AS8GUSZ5VW13PlY2tJiY3grTunEh7iv6tHekNnMzG7vAIXexHqaWBva/J2eAdY5Li9CHjbHYEqpb4qMiyYsWlxfSZ5g70BxMPfGMsXxaf506o8X4fjt5wpoUwDbgUuE5Edjq+rgAeBy0UkH5jruK+UUm4x94JkvnWJvQHE5wdO+iyOQ+W13LM8lxc3FnK4vBZvrh/VFV3MSinlt+oam7nmz+s409TCB3fNIC4qzKvHb2y2cd1f17O35PTZfqBp/SOZnp7INMdXfLTnY9LFrJRSAScqLITHFtgbQNz7xi6evGWCV4cWPvrvPL4oPs3fb8tmhDWadQXlrMsv571dxbyy5SgAowfEnk3ok4bFe7Xbk16BK6X83t8/O8jv3t/Lb68bwy1uaG7tjE0HT7Lg7xtZcPEgfn/D2K881txiY9fxKtbl25tkbztyiqYWQ1hIENlD+jMtPZHp6YnnNBbpKV0PXCkVsGw2w+LntrDp4Ene+dF0slI8u4b56fom5j22ltBg4b2fzCC6i9mtdY3NbD5UcTah7yuxD3+0RIYydUQC09ITuWZsao9LQJrAlVIBray6gXmPryU+OpS375xOZJjnShU/e3UHb+cW8fodU5gwuH+3n19W3cCGA/Zyy/qCcoqq6vn0F7N6PHZfa+BKqYBmjQnn0ZvHcevTm3ngvS88tmrhuzuLeHP7ce6ak9Gj5A32WNs2Fjl8so4hCe5vG6frgSulAsaMDCs/uHQ4/9x0hA92FXf9hG4qrjrDfW/tZvygOH50mXs6BIkIwxKjPfLhqyZwpVRA+cUVWYwbFMcv39jJsVN1bntdm83wi9dzaWy28ejN4wkN9v/06P8RKqVUG6HBQTyx4CJsBu5+ZQfNLTa3vO6zGw6zvuAk9197AcMCZJ0ZTeBKqYAzOCGK310/hpzCUzz+cb7Lr7e/pJqHVu5j7qhkFlw8yA0ReocmcKVUQJo/fiA3TUzjL6sLXJpq39Dcwl2vbCc2IoQHv3FhQK1BrglcKRWw/mf+aIYlRnP3q9upqG3s0Wv88aM89pVU8/CNY0nsF+7mCD1LE7hSKmBFhYXwxMKLOFXbxD3Lc7u90NSGA+X8fe1Bvn3JYC4bGXh9ODWBK6UC2ugBFn511Uj+vbeUZRsOO/28qromfv5aLkMTornv6lGeC9CDNIErpQLe4qlDmTMyif99fx97iqqces5/vb2b0uoGHrt5PFFhgTmnURO4UirgiQiP3DSO/tGh/Pjl7dQ1Np93/7d3HOed3CLumpPBuEFxXorS/TSBK6V6hfjoMB69eTyHymv573f2dLrf8coz/Oe/djNhcBw/nDXCixG6nyZwpVSvMXVEIj+anc5rOcd4J7fonMdtNsPPX9uBzWZ49ObxhATAbMvzCezolVKqnbvmZDBxSH9+/eYujpz86lT7f6w7yMaDFfy/a0czJCEwZluejyZwpVSvEhIcxOMLxiMCP35lO02OqfZfFJ3mDx/m8bXRydyUnebjKN1DE7hSqtdJ6x/FQ98YS+7RSv74UR71TS3c/ep2LFGh/P6GsQE12/J8AnPsjFJKdeGqC1NZOGkwT605wJ6iKvJO1PDcdy72ShNib+nyClxEnhGRUhHZ3WZbvIisEpF8x/eerXqulFIedP81F5CR1I+1+eXcNmUIs7KSfB2SWzlTQnkOuLLdtnuBj40xGcDHjvtKKeVXIsOC+dutE7nj0hH8al5gzrY8ny4TuDHmM6Ci3eb5wDLH7WXAdW6OSyml3GK4tR/3zhvp0R6avtLTDzGTjTGt/YxKgE5XgRGRJSKSIyI5ZWVlPTycUkqp9lwehWLsy391ugSYMWapMSbbGJNttVpdPZxSSimHnibwEyKSCuD4Xuq+kJRSSjmjpwn8HWCR4/Yi4G33hKOUUspZzgwjfBn4HMgSkWMicjvwIHC5iOQDcx33lVJKeVGXE3mMMQs7eWiOm2NRSinVDTqVXimlApQmcKWUClDS3SagLh1MpAwo7OHTE4FyN4bjbhqfazQ+12h8rvH3+IYYY84Zh+3VBO4KEckxxmT7Oo7OaHyu0fhco/G5xt/j64yWUJRSKkBpAldKqQAVSAl8qa8D6ILG5xqNzzUan2v8Pb4OBUwNXCml1FcF0hW4UkqpNjSBK6VUgPK7BC4iV4rIfhEpEJFzOv2ISLiIvOp4fJOIDPVibINEZLWIfCEie0Tkrg72mSUiVSKyw/F1v7ficxz/sIjschw7p4PHRUT+7Dh/O0Vkghdjy2pzXnaIyGkRubvdPl49f660DBSRRY598kVkUUf7eCi+R0Rkn+Pn95aIxHXy3PO+FzwY33+LyPE2P8OrOnnueX/XPRjfq21iOywiOzp5rsfPn8uMMX7zBQQDB4DhQBiQC1zQbp8fAk85bi8AXvVifKnABMftGCCvg/hmAe/68BweBhLP8/hVwAeAAJOBTT78WZdgn6Dgs/MHzAQmALvbbHsYuNdx+17goQ6eFw8cdHzv77jd30vxXQGEOG4/1FF8zrwXPBjffwO/cOLnf97fdU/F1+7xPwL3++r8ufrlb1fgk4ACY8xBY0wj8Ar29m1ttW3nthyYIyLijeCMMcXGmG2O29XAXmCgN47tRvOB543dRiCudW13L5sDHDDG9HRmrluYnrcM/BqwyhhTYYw5Bazi3N6xHonPGPORMabZcXcjkObu4zqrk/PnDGd+1112vvgceeObwMvuPq63+FsCHwgcbXP/GOcmyLP7ON7EVUCCV6Jrw1G6uQjY1MHDU0QkV0Q+EJHRXg3M3h3pIxHZKiJLOnjcmXPsDQvo/BfHl+cPnGsZ6C/n8bvY/6PqSFfvBU/6kaPE80wnJSh/OH8zgBPGmPxOHvfl+XOKvyXwgCAi/YA3gLuNMafbPbwNe1lgHPAE8C8vhzfdGDMBmAfcKSIzvXz8LolIGPB14PUOHvb1+fsKY/9f2i/H2orIfUAz8FInu/jqvfAkMAIYDxRjL1P4o4Wc/+rb73+X/C2BHwcGtbmf5tjW4T4iEgJYgJNeic5+zFDsyfslY8yb7R83xpw2xtQ4br8PhIpIorfiM8Ycd3wvBd7C/q9qW86cY0+bB2wzxpxo/4Cvz5+DMy0DfXoeRWQxcA3wbccfmXM48V7wCGPMCWNMizHGBvy9k+P6+vyFADcAr3a2j6/OX3f4WwLfAmSIyDDHVdoC7O3b2mrbzu1G4JPO3sDu5qiZPQ3sNcb8qZN9Ulpr8iIyCfs59sofGBGJFpGY1tvYP+za3W63d4DbHKNRJgNVbcoF3tLplY8vz18bzrQM/BC4QkT6O0oEVzi2eZyIXAncA3zdGFPXyT7OvBc8FV/bz1Su7+S4zvyue9JcYJ8x5lhHD/ry/HWLrz9Fbf+FfZREHvZPqO9zbPsN9jcrQAT2f70LgM3AcC/GNh37v9M7gR2Or6uAO4A7HPv8CNiD/VP1jcBUL8Y33HHcXEcMreevbXwC/NVxfncB2V7++UZjT8iWNtt8dv6w/yEpBpqw12Fvx/6ZysdAPvBvIN6xbzbwjzbP/a7jfVgAfMeL8RVgrx+3vgdbR2UNAN4/33vBS/G94Hhv7cSelFPbx+e4f87vujfic2x/rvU912Zfr58/V790Kr1SSgUofyuhKKWUcpImcKWUClCawJVSKkBpAldKqQClCVwppQKUJnCllApQmsCVUipA/X+7M5sdyYOTvgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\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]   # M0, R1=1/T1,R2=1/T2\n",
    "data      = forward(true_p)\n",
    "snr       = 50\n",
    "noise_std = true_p[0]/snr\n",
    "noise     = np.random.randn(data.size)*noise_std\n",
    "data      = data + noise\n",
    "\n",
    "# quick plot of the data\n",
    "plt.figure()\n",
    "plt.plot(data)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have the data and our forward model we are almost ready to begin fitting.\n",
    "\n",
    "We need a cost function to optimise. We will use mean squared error. \n",
    "\n",
    "In this section:\n",
    "\n",
    "- '**' operation\n",
    "- np.mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can set up our optimisation.\n",
    "\n",
    "In this section:\n",
    "\n",
    "- scipy minimize\n",
    "- dictionary\n",
    "- keyword arguments\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ready to minimize\n",
    "p0 = [200,1/1,1/70]  # some random initial guess\n",
    "method = 'powell' # pick a method. scipy has loads!\n",
    "\n",
    "\n",
    "kw_args = {'x0':p0,'method':method}\n",
    "result  = minimize(cf,**kw_args)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the data with the model prediction.\n",
    "\n",
    "In this section\n",
    "\n",
    "- printing\n",
    "- text formatting\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fitted = [1.02562035e+02 1.16194491e+00 1.98071179e-02]\n",
      "true   = [100, 1.25, 0.02]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xUZdbA8d9JJZRJ6JCgUoKgBAWNFTsoWAlBsYtlF8vquu6Ka1tXXXflXfZd674iVlRWQQVEFJFiL2goShdEERJKBEINkPK8fzwzEEJCJjN35t6ZnO/nw2dm7ty593AJJ8+c+xQxxqCUUir2JLgdgFJKqdBoAldKqRilCVwppWKUJnCllIpRmsCVUipGJUXzZK1atTIdO3aM5imVUirmzZkz51djTOvq26OawDt27EhBQUE0T6mUUjFPRFbVtF1LKEopFaM0gSulVIzSBK6UUjFKE7hSSsUoTeBKKRWjotoLpSGaNK+QkdOWUVRSSmZGGsP7dyOvd5bbYSml4oAm8AiaNK+QeyYsoLSsAoDCklLumbAAQJO4UipsWkKJoJHTllFaVkETSmlNCQClZRWMnLbM5ciUUvFAW+ARVFRSSmN28VbKQwCcu2fE3u1KKRUuTeARlJWeyr07H+eIhF/YZZIBAwiZGWluh6aUigNaQomg5zp9xHmJ37Ck8lAaSRkZbCctOZHh/bu5HZpSKg5oAo+UJe9yxNKn+eWQgbyWeikAvXw7eDS/p97AVEo5QhN4JKxfBBNuhKxcDr1mNH8fOgCAlwdnafJWSjlGE7jTdm6C1y+H1GZw6WuQ3Ah8mfa9rYXuxqaUiit6E9NJFeXw5lDYtg6uex987e32pm1BEjWBK6UcpQncSR/eBz99CnmjoEPuvu0JidCsPWwtci82pVTc0RKKU+a+CrNHwYm/g16XH/i+L1Nb4EopR2kCd8Ivs2HKHdD5TDj74Zr38WXCFk3gSinnaAIP15ZCGHcVZBwCl7wEibVUpdI72BKKMdGNTykVtzSBh6OsFMZdaR8vex3Smte+ry8TykuhdHP04lNKxTVN4KEyBib/Hormw+DnoE33g++vXQmVUg7TBB6qL5+CBePhrPug27l17+/rYB+1J4pSyiGawEOxfAbM+CscmQen3hncZ7QFrpRyWJ39wEWkGzCuyqbOwAPAK/7tHYGfgSHGmLgr8FZfUeehk1Po98X10LYH5P0fiAR3oKZtQRK0Ba6UckydLXBjzDJjTC9jTC/gWGAnMBG4G5hpjOkKzPS/jiuBFXUKS0oxwNaSjXSe+Vt2m0S47L+Q0iT4gyUm2cE82pVQKeWQ+pZQ+gI/GmNWAQOBMf7tY4A8JwPzgsCKOgAJVPJE8tMcwnru4I+QcWj9D6iDeZRSDqpvAr8MeN3/vK0xZq3/+TqgbU0fEJFhIlIgIgXFxcUhhumOqivn3Jk0nrMS5/Ng+VCmbu0c2gF9mVpCUUo5JugELiIpwEXAm9XfM8YY7HIzBzDGjDbG5Bpjclu3bh1yoG4IrJxzfsLX3JI0mdfK+zK2ol/oK+r4OtgWuA7mUUo5oD4t8HOBucaY9f7X60WkPYD/cYPTwblteP9uNEs23Js8lvmVnXmofGh4K+r4MqFsJ+wqcTZQpVSDVJ8Efjn7yicAk4Gh/udDgXecCsor8npn8dJxa8iSjTxVnk+bjGbhraiztyuhllGUUuELajpZEWkCnA3cWGXzCGC8iNwArAKGOB+ey4wht+g1aNWNF265DxLC7Daf7h/Ms6XQdkNUSqkwBJXAjTE7gJbVtm3E9kqJXys/hnUL4KKnwk/eoIN5lFKO0pGYB/Plk3YAzlGXOnO8pu10MI9SyjGawGuzbiH8OAtOuBGSUp05ZmKSTeLaAldKOUATeG2+fAqSm0Du9c4eVwfzKKUcogm8JlvWwMK34JhrDj7Hdyh0MI9SyiGawGsye5QdbHPizc4fO72D7YWig3mUUmHSBF7dri1Q8DL0yIPmhzl/fF8mlO2w51FKqTBoAq9uzsuwZxuc/PvIHF8H8yilHKIJvKryPfD1KOh0GmT2isw5fP5RnHojUykVJk3gVS18G7YVRa71DZrAlVKO0QQeYIztOtjmSMjuF7nzNGsHiJZQlFJh0wQesGImbFgEJ98W/DJpoUhMtqM7tQWulAqTJvCAL5+0S57lXBz5c6Vn6dJqSqmwaQIHKJoPP30CJ9wESSmRP58O5lFKOUATONjad0pTOPba6JzPl6UJXCkVNk3gJb/Aook2eadlROecvizb11wH8yilwqAJ/Otn7E3LSAybr40O5lFKOaBhJ/DSzTBnDPTI37daTjRoX3CllAMadgIveMnOS9InggN3apLuT+DaE0UpFYaGm8DLd9tZBzufCe16RvfcTXUwj1IqfA03gX8/HravtwN3oi0pBZq20RKKUiosDTOBV1baroNte0KXs9yJwZelCVwpFZaGmcBXTIdfl0V+2PzB6GAepVSYgkrgIpIhIm+JyFIRWSIiJ4lICxGZLiLL/Y8Orz0WQV88aVvAOfnuxaCDeZRSYQq2Bf4E8IExpjtwNLAEuBuYaYzpCsz0v/a+wjmw6nPb7zsx2b040rNg91bYtdW9GJRSMa3OBC4i6cBpwAsAxpg9xpgSYCAwxr/bGCAvUkE66sunINUHxwx1N469fcG1Fa6UCk0wLfBOQDHwkojME5HnRaQJ0NYYs9a/zzqgbU0fFpFhIlIgIgXFxcXORB2qTT/B4nfssPlGPndj2TsaU29kKqVCE0wCTwKOAZ4xxvQGdlCtXGKMMUCNy6wbY0YbY3KNMbmtW7cON97wfP0MSGJ0h83XRlvgSqkwBZPA1wBrjDGz/a/fwib09SLSHsD/uCEyITpk5yaY9yr0vGRf69dNzdrbR22BK6VCVGcCN8asA1aLSDf/pr7AYmAyECgkDwXeiUiETvn2BSjbCSff6nYkVlIKNNHBPEqp0CUFud9twFgRSQFWAtdhk/94EbkBWAUMiUyIDijbBd88a9e6bNvD7Wj2SdeuhEqp0AWVwI0x84HcGt7q62w4EfL9G7Cj2J1h8wfjy4KNP7odhVIqRsX/SMwta2DGQ5B5DHQ63e1o9qejMZVSYYjvBF6+B968FirKIH+0e8Pma+PLgt1bYPc2tyNRSsWg+E7g0/8Ca76FgU9Bq65uR3Mg7UqolApD/CbwhRPsfN8n3Aw9BrkdTc10MI9SKgzxmcCLf4DJt0GH4+Hsh92Opnbp2gJXSoUu/hL4nh0w/hpISoVLXrb9rb0qMJhHl1ZTSoUg2H7gscEYmHIHFC+Fqyfsa+F6VVIqNGmtJRSlVEjiK4HPeQm+Hwdn3Lt3pZ1J8woZOW0ZRSWlZGakMbx/N/J6eyixa1dCpVSI4ieBF86FqX+GLn3htOGATd73TFhAaVmF3aWklHsmLADwThL3dYDNP7kdhVIqBsVHDXznJhg/1M4tkv8cJNi/1shpy/Ym74DSsgpGTlvmRpQ182VqCUUpFZLYT+CVlTDxJti2FoaMgSYt975VVFJa40dq2+4KXybs2gK7t7sdiVIqxsR+Av/iMVg+Dfr/AzrsP11LZkZajR+pbbsr0jvYR62DK6XqKbYT+E+fwqxHIGcwHP/bA94e3r8bacmJ+21LS05keP9uB+zrGh3Mo5QKkedvYtbai2TrWnjremjZFS58ssZ5TgI3Kj3fCwW0Ba6UqjdPJ/DaepFIZRkDv7sJ9uyEoVMgtWmtx8jrneWthF1dM03gSqnQeDqB19aLZMf7D0DFV5D/PLTp7lJ0DkluBI1bwdY1bkeilIoxnq6B19RbpH/Ct1xR8Q4c9xs46hIXoooAHcyjlAqBpxN49d4ih8k6RiaPYrFk214n8SK9gyZwpVS9eTqBV+1Fksoenkl+gkoSWd3vGTuPSLzwZdqVg5RSqh48XQOv2ovk99tHc2TCKr48cRT9+xzvcmQO82XCrhI7k2JKE7ejUUrFCE+3wMEm8S/OKeTSpI/htOGcPOByt0Nyni8wmGetu3EopWJKUAlcRH4WkQUiMl9ECvzbWojIdBFZ7n9sHpEIjYHlH9oFic+4JyKncN3evuBaRlFKBa8+LfAzjTG9jDGB8ep3AzONMV2Bmf7XzhOBS8bAZWMhIbHu/WORDuZRSoUgnBLKQGCM//kYIC/8cGqRkACpzSJ2eNftXdxYh9MrpYIXbAI3wIciMkdEhvm3tTXGBIq264C2NX1QRIaJSIGIFBQXF4cZbpxKbgSNW+rSakqpegm2F8opxphCEWkDTBeRpVXfNMYYETE1fdAYMxoYDZCbm1vjPgodzKOUqregWuDGmEL/4wZgInA8sF5E2gP4HzdEKsgGwaeDeZRS9VNnAheRJiLSLPAcOAdYCEwGhvp3Gwq8E6kgGwRfpvZCUUrVSzAllLbARLHTtSYB/zXGfCAi3wLjReQGYBUwJHJhNgC+TCjdbGdYTGnsdjRKqRhQZwI3xqwEjq5h+0agbySCapACK/NsWwstu7gbi1IqJnh+JGaDoSvzKKXqSRO4VwT6gmtXQqVUkDSBe4W2wJVS9aQJ3CuS0yCthXYlVEoFTRO4l/iytAWulAqaJnAvSdcErpQKniZwL9Hh9EqpetAE7iW+TNi5EcoOXMxZKaWq0wTuJXtX5tFWuFKqbprAvUQXdlBK1YMmcC/RhR2UUvWgCdxLdDCPUqoeNIF7SUpjSGuuJRSlVFA0gXuNL0vnQ1FKBUUTuNf4MrWEopQKSrBrYqpo8WVB4Vy3o3DMpHmFjJy2jKKSUjIz0hjevxt5vbPcDkupuKAJ3Gt8WbDzVyjbZVerj2GT5hVyz4QFlJZVAFBYUso9ExYAaBJXygFaQvGaQE+UbbF/I3PktGWUllWQQhkZbAOgtKyCkdOWuRyZUvFBE7jXpAf6gsd+Ai8qKeU4Wcr0lOF8kHo3iVTs3a6UCp8mcK/xxUkCLyvl0SavMy7lb7SQbbSTzRwtPwKQmZHmcnBKxQdN4F4TKKFsWeNuHOFYUwCjTuWyind5w5zNObv/SYURTk/8nrTkRIb37+Z2hErFhaATuIgkisg8EZnif91JRGaLyAoRGSciKZELswFJaQKNMmKzBV6+G2Y8BC+cbWdUvHoSjQc9TkJGB74zXeiXvIBH83vqDUylHFKfXii3A0sAn//1/wCPGWPeEJFRwA3AMw7H1zD5smIvga/9HibeBBsWQe+roP8/oFE6efh7nHx8CXw8gh6Hx3bPGqW8JKgWuIh0AM4Hnve/FuAs4C3/LmOAvEgE2CD5MmFrjJRQKsrgk3/Cc2fa7o9XjIeB/4FG6fvvl90PMLDyI1fCVCoeBVtCeRy4C6j0v24JlBhjyv2v1wA1fi8WkWEiUiAiBcXFxWEF22Ckx0gLfMNSeL4ffPR36DEIbvkaDu9f876Zve08LytmRjdGpeJYnQlcRC4ANhhj5oRyAmPMaGNMrjEmt3Xr1qEcouHxZcGOYltT9qLKCvjiCXj2NNiyGoa8AoOfh8Ytav9MQiJ0OQtWzIDKytr3U0oFLZgaeB/gIhE5D2iErYE/AWSISJK/Fd4B0Ak8nFJ1YYcWndyNpbqNP8Kkm2H1bDjiQjj/MWga5C/m7H6w8G1YvxDaHxXZOD1CpxJQkVRnC9wYc48xpoMxpiNwGTDLGHMl8BFwsX+3ocA7EYuyofFiX/DKSpg9Gp7pA8VLIf85GPJq8MkbbAscbCu8AQhMJVBYUoph31QCk+ZpW0c5I5x+4H8G/igiK7A18RecCUl5LoFvWQOvXARTh0PHU2yt+6ghIFK/4zRrB217Npg6eGAqgap0KgHlpHpNZmWM+Rj42P98JXC88yGpfSUUD/REKS2BVwbCtnVw0VPQ++r6J+6qsvvCV0/Drq3QyFf3/jGstikDdCoB5RQdielFqU1tNzy3W+AV5ax/6UrKNv7MpdvuoM+HHZg0P8yYsvtBZTn8/JkzMXpYbVMG6FQCyimawL3KA4N5Voy9g7YbPuf+suuYbY5wpoZ7yAmQ0rRB1MGH9+9GWnLiftt0KgHlJE3gXuXLdHc+lLmvkr3yFV4q78+4ijP3bg67hpuUAp1OtwncGAcC9a683lk8mt+TrIw0BMjKSNOpBJSjdEEHr/Jl2eHpblj1FUy5g88qevJI+VUHvB12DTe7Lyx7DzaugFZdwzuWx+X1ztKErSJGW+Be5cuCHRugfE90z1vyC4y7Cpofxt8b30UFiQfsEnYNN7uvfWwAZRSlIkkTuFftXZlnbfTOuXs7vH65nd/k8je4acCxkanhNu8ILbM1gSsVJi2heNXelXkKoflhIR8m6JGAlZUw8UbYsBiufBNadSWvlX0rIiMJs/vBnJfttLPJ2itDqVBoAvcqBwbz1GtR4Y//AUunQP9H/TMHsne/iNRws/vB7FGw6st9JRWlVL1oCcWr9g7mCb3LXtAjARe+DZ+OtPN4n3hzyOerl8P6QGJqgxmVqVQkaAL3qtRmkJoOW0JP4EGNBCycC5NugUNPgvP/Hd4oy/pIaQwd+2gdXKkwaAL3Ml9mWC3wOkcCbl0Lb1wBTVrbiamSUkM+V0iy+8Gvy2zPl3hW8ovtmqmUwzSBe5kvM6wa+EFHApaVwrgr7Zwkl79ev1kFndIl0J0wTssopSUw/QF4KhdeOtf1kbUq/mgC97L0rLBa4LWOBOyVCZN/D4VzIP9ZaNfTuZjro3U38HWIvzJKRRnMfhae7G0Xvuh6NmBg0SS3I1NxRnuheJkvC7b7B/MkpYR0iBp7kXz+GCwYD2febxdlcIuI7YGycIJNeonJ7sXiBGNg6Xu21b3pR+h4KpzzCGT2glGn2pvFJ93idpQqjmgL3Mt8mYCB7eucO+ayqTDjIeiRD6fd6dxxQ5XdD/ZsgzXfuh1JeArnwEvn2bJUQiJcPg6GvmuTN0BOPhQWwOafXQ1TxRdN4F4W6EoYRk+U/axfDG//BtofbVeOj1aPk4PpfDpIomfLKJPmFdJnxCw63f0efUbMOnAmxpJf7DV97iz49Qfbk+fmr6DbgP2vb498+7hoYvSCV3FPE7iX+TrYxzDq4Hvt2AivXwopTexNy5TG4R/TCY3S7RSzHkzgB10SbdeWfTcol7wLp/4Jfj8PjrsBEmuoTDY/DLJybRlFKYdoAveyqosbh6N8D4y/Brath8v+u++4XpHdF9Z+Z+v9HlLTQKiyst2sfO8xeKKXvUGZkw+3zYG+D9S9wlDOYFi3AH5dHsGoVUOiCdzLGvkgpVl4LfCdm+DNa2HV53ZJtA65joXnmMBQ+h9nuRtHNfsPhDKck/AtH6bcxR/Ln4O2PWDYJzBoFKR3CO6APfIAsTdtlXKAJnCvC6cr4c+fw6hTYPk0GDACjr7U2dic0u5oaNzKc2WUwICnzlLEuJS/MTrlMSpJYHjyvfvfoAyWLxMOO9mWUeJ8MQsVHZrAvS6UwTwVZTDzb/DyBZDUCG6YHr05TkKRkGBb4T/OsrMiesTw/t04KrmQ8SkP01XWcH/ZdeSZf9HnvKtCvwGck29Hn25Y7GywqkGqM4GLSCMR+UZEvhORRSLykH97JxGZLSIrRGSciITWUVkdnC+zfr1QNv0ELw6Az/4Fva+EGz+FrGMiF59TsvvBzo2wdr7bkeyVl7mFN9MeBUnkkj0P8lGzi3gkv1d4szMeMRAkQW9mKkcEM5BnN3CWMWa7iCQDn4vIVOCPwGPGmDdEZBRwA/BMBGNtmHwdYPv64Aa6fDcO3vuTTRAXv2Rbe7Giy1mA2GH1XviFs34RjLmQ1NRUUodNYWarbGeO27S1XRN04QQ46y/e6MqpYladLXBjbfe/TPb/McBZwFv+7WOAvIhE2NAFBvNsO8hgnl1b4e3fwsRh0C4Hbv48tpI3QJNWtqbshTr4ugW2/JSYCte+B04l74CcfNj8ExTNc/a4qsEJqgYuIokiMh/YAEwHfgRKjDHl/l3WALpyayTUtbDD6m/tjcqFb8OZ98HQKZBxaPTic1J2P1jzDZRudi+GdQtgzEV2laBrp0DLLs6fo/sFkJAMi7Q3igpPUAncGFNhjOkFdACOB7oHewIRGSYiBSJSUFxcHGKYDdjepdXW7L+9sgI+GQkv9gcMXDcVTr+r5kEksaJLXzCVsPITd86/9nsYcyEkN45c8gZo3MI/B8xET920VbGnXr1QjDElwEfASUCGiASyRQegxjttxpjRxphcY0xu69YuTFka62oazLNljU00Hz1i+xbf9DkceoI78Tmpw3F2EQs3yihF8+01TWlqk3eLzpE9X498+0s51ueAUa4KphdKaxHJ8D9PA84GlmAT+cX+3YYC70QqyAYt1WeTSiCBL5oEz5xsRy7mjYLBL9jh6PEgMcnOjbJiZnT7SRfNg1cG2mt97RRo0Sny5+x2ru3iqb1RVBiCaYG3Bz4Ske+Bb4HpxpgpwJ+BP4rICqAl8ELkwmy4Js0v4qeyDD764gsm/20wvDkUWnSx3QN7XR5/vRiy+8G2IiheGp3zFc7dP3k37xid8zby2XnCF0+y5TClQlBnwdQY8z3Qu4btK7H1cBUhgcmUnqU5ZyZ+R2W58KzJo92xDzGwZUe3w4uMwLD6FTOgzRGRPVfhHHh1kP0GM3SKnXAqmnIG24mwVn0BnU6L7rmjpXw3fD/e/kLu91Bs36PxIB2J6WGByZTmmq6srmzNlWX38uieIfxz+kq3Q4uc9A7Q+ojI18HXzIFXBkGjDNtVMNrJG6Brf0huEp9llNLN8Nn/wuM9YfKt8NXT8POnbkcVdzSBe1hgMqXHyy/m1D1P8FVlj/22x63svrDqS9izIzLHX1MAr+ZB4+Y2ebvV7TKlsa2FL55sB2rFg5JfYOrd8O8eMPNhO+nXFW/aSdl0Ei/HaQL3sDpXlY9X2X2hYo+djMtpq7+1ZZPGLfzJ+xDnz1EfOYOhdJN7XSedUjQf3rreTrP77XN2qb6bvoCrJ8Lh50D382HJZDu1sXKMJnAPO+iq8vHs0JMhKc35Msrqb/zJuyVc+37w08BGUnZf23UyFgf1GAPLZ9jul6NPhx8+tJOm3f6df7HsnH375uTbRTBWfuRevHFI7yh4WGDSpJHTllFUUkpmRhrD+3cLbzKlWJDcCDqd6mwC/+VreG0wNG1re5t4ZVGLpFQ44gJ7M/OCx+xrryvfAwvfgi+fsrMqNmsPZz8Mx15be5fWzmfa+w0L34bD+0c13HimCdzjalxVviHI7gfLP4RNK8MfVPPDNPv13mvJO6BHPswfa39hdT/f7Whqt2sLFLwEs0fBtrXQ5kjIewZyLoakmicjnTSvcG8D5MkmuQxY9C7JF5baqQpU2LSEorwpu599XDEz9GNsWgn/vQz+OwTSD7E1b68lb7CDl9JaePcm385NMO0+e2Nyxl+hVVe48m24+UvodcVBk3fVNUXfKD2O5IqdzP5wXHTjj2PaAlfe1KKzHVSzYiYc/9v6fXbPDvjs3/YrfkIS9HsQTrzFu+WJxGQ48iL4/k3Ys9M7C06DnVb39cvsnPQ9BsHJtwW9ElH1NUW/rjySX42P7XPGwfnXRijghkVb4MqbROzkVj99ageDBMMY24p9+ni7oMWRF8FtBXDKHd5N3gE5g6Fsh13+zisWT6Z8dF9+LdnGoF1/pc+KK5m0Pvj5jKp3d60gkfcrTuDkigLYvb2WT6n60ASuvCu7n01qv3xd977rF9veEG9dB2nN7eyMg5/3ZsmkJof1sTV6LwzqqayEjx6F8VezqDyL83Y9wjyTTWFJKfdMWMCkecGtEFVTd9cpFSeSJnvghw+cjrpB0gSuvKvTqVRKEmPHvkinu9+jz4hZByaP0hKY+mc7J/q6BXDev+DGT+ziwbEkIRGOzIPl0+0CHW7ZvR3GXw2fjGBqwhkM2X0/G2i+9+3SsgpGTlsW1KFq6ga7MOlIShu18W69P8ZoAleeNWnxVr6p7MYxZXMwsH8LsLIS5r4CTx0Ls5+FY66B2+baenlCYp3H9qScwVC+C5ZNdef8m36CF86GZe9D/0e5Zedv2c2BNyiDHQmc1zuLR/N7kpWRhgBZGWn8I/9o0npdDCum21++Kix6E1N51shpy7ig/CjuSX6dtmxiPS0oLatgytR3yft2HBTNhUNOgKveDvrGmqd1OM6ugbrwbTj60uiee+UndqZLY+z17HIWmZ/MorCGZF2fkcA1doNdMxi+/j/7i6LXFeFG3qBpC1x5VlFJKZ9UHg3AaYnf04ot/DPpWZ7f82fYWgiDnoXrp8VH8gZISICcQfDjLNt1LxqMsd9gXh0ETdrAb2f5F5iO4EjgrGPt/DNeqPfHOE3gyrMyM9JYag5hvcngxsQpzEr9I4MSP+e/iQPh1gI4+rL4mw+9Rz5UlsHSKZE/V/lumHwbTL0Lup4Dv5mx3zJyNZVAHs3vGf7AMhH791z5MezYGN6xGjgtoSjPGt6/G/dMWMBHFb24LOljPq3oyQiuY9gFA+yCCEGoOhIwJqYiyOwNzTvZm3zHXBO582xbb29Wrp4Np95pF8ROOLA9F7GRwDmD4YvH7QRXudc5f/wGQhO48qxA4nj2g+t4c+sZrPMdxfAB3YNOKIGRgIHBJIGboFWP7TkiNrl9/m/YXgxNI7CObOFcGHeVnbP74pfsRFPR1q4ntMy2k3hpAg+ZJnDlaeG0AKuPBIR93eA8m8DBJtTP/gVL3oHjfhPWoap/A3mix3Jyv3sAmrS29w/aH+VQ0PUU+EX16Uj7baBZW3fiiHFaA1dxq7bubp5fEKPNkdC6e9h9pavORSJUcvX2F8mdcxe/pufAsI/dS94BPfLBVMJiXQ89VJrAVdyK2QUxAjf5Vn0JW4tCPkzgG0hztvJ88r+4KeldXivvy+Dtd0GTVg4GHKI23aFND+2NEgZN4CpuxfSCGDn5gIFFk0I+RFlJEfcmjeWL1Ns5NWEB95Vdz/3lN/DLlnLn4gxXziBY/TVsWeN2JDFJE7iKWxHrBhcNrbraG32htE5LVsN7d/JZoz9wQ+L7fFB5HAP2jGBshZ2i11PfQHr4b6AumuhuHDGqzpuYInII8ArQFjDAaGPMEyLSAhgHdAR+BoYYYzZHLlSl6i+mF8TIGQwzHoTNP9updeuy8Ufbe+W7NwBh7f+uwSMAAAwBSURBVGF5DFt5Cj+U7evJ4rlvIC27QPte9hfVybe5HU3MCaYFXg78yRhzJHAi8DsRORK4G5hpjOkKzPS/Vko5JdjW6YYl8PZv4OlcWPAW5N4At8+n43XPc0v+2d7/BpIzGIrm2QU4VL2IMaZ+HxB5B3ja/+cMY8xaEWkPfGyMOeiv9tzcXFNQUBBysEo1NJueOJXiku0MKH3kwIFIRfP93Q3fheQmcNwNcNKtsdclr2Q1PJ4DZ/0FTrvT7Wg8SUTmGGNyq2+vVz9wEekI9AZmA22NMWv9b63Dllhq+swwYBjAoYceWp/TKdWgTZpXyJKNvbgn4RU6SRErSzK5Z8ICWmyax2lrX7Yz+qWmw2l32dXgG7dwO+TQZBxiJyVbNFETeD0FfRNTRJoCbwN/MMbsN2Gxsc34GpvyxpjRxphcY0xu69YRGFWmVJwaOW0Zk/YcT6URLkj4mpMSFvECD3HaZ1dC4RzbYr1jAZx1X+wm74Ae+bB+IRQHN9e4soJK4CKSjE3eY40xgdEF6/2lE/yPGyITolINU1FJKetpwbemG7clTeT1lL+TnVDEI2VXwh0LbWu1UbrbYTqjRx4gnlvoYdK8QvqMmFX7giIHY4wtc330Dxh9Buza4nh8wfRCEeAFYIkx5t9V3poMDAVG+B91OJVSDsrMSKOwpJQXy8/Fl7STseV9ebPidFplpHN/ShO3w3NWs3bQ8RTbG+WMu/fOMunmZGQhzaVTvht+/swuyrFsqp32WBJsiWj7Bsd/4QZTA+8DXA0sEJH5/m33YhP3eBG5AVgFDHE0MqUauMBsjNPKjmPanuMAD3YDdFJOPky5w5ZS2vV0fTKyoOfSKd1sl8Jb+h6smAl7tkFyYzuv+pn3weH9Izbytc4Eboz5HKht0uW+zoajlAoIJImYmg43HEcMhPfutK3wdj1dn4zsoHPpbP4Zlr5vVxVa9SWYCrsgRk4+dDsPOp8OyZEfMKWzESrlYTE9EKm+mrSEzmfYOnjfv7o+GVmghAUgVNJTfuLsxDmcmzwPnlhld2rdHfrcDt3Ph8xjDphTPdIlIE3gSinvyBkM79wChXP3S6BVRWsqgOH9u/HAhLkMqZzKb5Lep51spsIIm1vkQu6NcPiA/VYwqi4aJSCdC0Up5R3dz4fEFFg0wd3JyCoryUv8kq+a/Zn7k8eyojKLh5Nv54PzvqDVbTPgpN8dNHnDwWvoTtEWuFLKO9IyILsfLJxA3h1/A1y4B/DTp/DhX2DtfJq0OwoG/4dTupzJKfU8TDRKQJrAlVLe0iPf3hxcPZu83idF7x7A+sUw46+w/ENIPwQGjYael9S4VmgwolEC0hKKUspbup0LSWnRW+hhaxG8cyuM6mMXeT77b3BrARx9acjJG6IzH722wJVS3pLaFA4/BxZPggEjIDFCaWrXVvjiCfjqP7Yb4Im3wKl/cmxagmh0A9UErpTynpzBdq3MVZ/broVOqiiDOS/DxyNg56+QczH0/Utwc67XU6S7gWoCV0p5T9dzIKWpLaN0PsOZYxpjp96d8SBs+hEOOwXOeRiyjnXm+C7QGrhSynuS0+yIxiXvQvme8I/3y2x4sT+MvxoSk+GK8XDtlJhO3qAtcKWUV+Xkw4LxsPJjWxOvr1+X2zLMksmw9jto2g4ufBJ6XRm5unqUxcffQikVf7qcZWfvWzQhuARuDGxYDIsn28RdvMRu73Ac9P8HHHstxNksjprAlVLelJQKR1wIi96BC3ZBcqMD9zEG1s63SXvJZNi4AhA47GQ495/Q/QJIj9+5ZDSBK6W8q0c+zHsNVsyAIy6w2yorobBgX3mk5BeQROh0qh3i3v0CaNrG3bijRBO4Usq7Op0OjVvaWnhac3/Sfhe2FUFCMnQ5064J2v382F9WLgSawJVS3pWYBEcOhIIXbfJOamTnSjniQbtQQlqG2xG6ShO4UsrbTroVTKVtjXc9x47UVIAmcKWU17XsAhc+4XYUnqQDeZRSKkZpC1wpFbfcXNU+GjSBK6Xiktur2kdDnSUUEXlRRDaIyMIq21qIyHQRWe5/bB7ZMJVSqn6isaSZ24Kpgb8MDKi27W5gpjGmKzDT/1oppTzD7VXto6HOBG6M+RTYVG3zQGCM//kYIM/huJRSKiy1LV0WrVXtoyHUXihtjTFr/c/XAW1r21FEholIgYgUFBcXh3g6pZSqH1dXtY+SsLsRGmMMYA7y/mhjTK4xJrd169bhnk4ppYKS1zuLR/N7kpWRhgBZGWk8mt8zbm5gQui9UNaLSHtjzFoRaQ9scDIopZRyQqSXNHNbqC3wycBQ//OhwDvOhKOUUipYwXQjfB34CugmImtE5AZgBHC2iCwH+vlfK6WUiqI6SyjGmMtreauvw7EopZSqB50LRSmlYpQmcKWUilFiewFG6WQixcCqED/eCvjVwXCcpvGFR+MLj8YXHq/Hd5gx5oB+2FFN4OEQkQJjTK7bcdRG4wuPxhcejS88Xo+vNlpCUUqpGKUJXCmlYlQsJfDRbgdQB40vPBpfeDS+8Hg9vhrFTA1cKaXU/mKpBa6UUqoKTeBKKRWjPJfARWSAiCwTkRUicsBKPyKSKiLj/O/PFpGOUYztEBH5SEQWi8giEbm9hn3OEJEtIjLf/+eBaMXnP//PIrLAf+6CGt4XEXnSf/2+F5FjohhbtyrXZb6IbBWRP1TbJ6rXL5wlA0VkqH+f5SIytKZ9IhTfSBFZ6v/3mygiGbV89qA/CxGM70ERKazyb3heLZ896P/1CMY3rkpsP4vI/Fo+G/HrFzZjjGf+AInAj0BnIAX4Djiy2j63AKP8zy8DxkUxvvbAMf7nzYAfaojvDGCKi9fwZ6DVQd4/D5gKCHAiMNvFf+t12AEKrl0/4DTgGGBhlW3/BO72P78b+J8aPtcCWOl/bO5/3jxK8Z0DJPmf/09N8QXzsxDB+B4E7gzi3/+g/9cjFV+19/8XeMCt6xfuH6+1wI8HVhhjVhpj9gBvYJdvq6rqcm5vAX1FRKIRnDFmrTFmrv/5NmAJEGuTDQ8EXjHW10CGf073aOsL/GiMCXVkriNM6EsG9gemG2M2GWM2A9M5cO3YiMRnjPnQGFPuf/k10MHp8warlusXjGD+r4ftYPH588YQ4HWnzxstXkvgWcDqKq/XcGCC3LuP/4d4C9AyKtFV4S/d9AZm1/D2SSLynYhMFZEeUQ3Mro70oYjMEZFhNbwfzDWOhsuo/T+Om9cPglsy0CvX8XrsN6qa1PWzEEm3+ks8L9ZSgvLC9TsVWG+MWV7L+25ev6B4LYHHBBFpCrwN/MEYs7Xa23OxZYGjgaeASVEO7xRjzDHAucDvROS0KJ+/TiKSAlwEvFnD225fv/0Y+13ak31tReQ+oBwYW8subv0sPAN0AXoBa7FlCi+6nIO3vj3/f8lrCbwQOKTK6w7+bTXuIyJJQDqwMSrR2XMmY5P3WGPMhOrvG2O2GmO2+5+/DySLSKtoxWeMKfQ/bgAmYr+qVhXMNY60c4G5xpj11d9w+/r5rQ+UlaT2JQNdvY4ici1wAXCl/5fMAYL4WYgIY8x6Y0yFMaYSeK6W87p9/ZKAfGBcbfu4df3qw2sJ/Fugq4h08rfSLsMu31ZV1eXcLgZm1fYD7DR/zewFYIkx5t+17NMuUJMXkeOx1zgqv2BEpImINAs8x97sWlhtt8nANf7eKCcCW6qUC6Kl1paPm9evimCWDJwGnCMizf0lgnP82yJORAYAdwEXGWN21rJPMD8LkYqv6j2VQbWcN5j/65HUD1hqjFlT05tuXr96cfsuavU/2F4SP2DvUN/n3/Yw9ocVoBH2q/cK4BugcxRjOwX7dfp7YL7/z3nATcBN/n1uBRZh76p/DZwcxfg6+8/7nT+GwPWrGp8A//Ff3wVAbpT/fZtgE3J6lW2uXT/sL5K1QBm2DnsD9p7KTGA5MANo4d83F3i+ymev9/8crgCui2J8K7D148DPYKBXVibw/sF+FqIU36v+n63vsUm5ffX4/K8P+L8ejfj8218O/MxV2Tfq1y/cPzqUXimlYpTXSihKKaWCpAlcKaVilCZwpZSKUZrAlVIqRmkCV0qpGKUJXCmlYpQmcKWUilH/D+HkpOq5JKufAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optional: use gradients and hessians to help with the optimisation\n",
    "\n",
    "In this example the forward model is simple enough that calculating the derivatives of the cost function is relatively easy to do analytically. Below is an example of how you could define these and use them in the fitting\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# gradient of the forward model\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",
    "# hessian of the forward model\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",
    "    dfdM0dM0 = np.zeros(E1.shape)\n",
    "    dfdM0dR1 = E2*(-dE1)\n",
    "    dfdM0dR2 = dE2*(1-E1)\n",
    "\n",
    "    dfdR1dM0 = E2*(-dE1)\n",
    "    dfdR1dR1 = M0*E2*(-ddE1)\n",
    "    dfdR1dR2 = M0*(dE2)*(-dE1)\n",
    " \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": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fitted = [1.02576306e+02 1.16153921e+00 1.98044006e-02]\n",
      "true   = [100, 1.25, 0.02]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xUZfb48c9Jg1CSUCKQoFKCKAQEjKigiFKthGBfFZW17K5tdXHB3Z9lXdfC17XuyqKyYgVBmqIigqwVNBQpIgooQmhRSGgBUp7fH88EQkjIZObO3DuT8369eM3MnTtzD0M4ufPc5zxHjDEopZSKPDFuB6CUUiowmsCVUipCaQJXSqkIpQlcKaUilCZwpZSKUHHhPFjz5s1NmzZtwnlIpZSKeIsWLfrFGJNaeXtYE3ibNm3Izc0N5yGVUiriicj6qrbrEIpSSkUoTeBKKRWhNIErpVSE0gSulFIRShO4UkpFqLDOQqmLpi/JY8zs1WwqKCItJZGRgzqS3T3d7bCUUlFAE3gITV+Sx+ipyykqLgUgr6CI0VOXA2gSV0oFTYdQQmjM7NUUFZfSmL20YDsARcWljJm92uXIlFLRQM/AQ2hTQRENKWJKwgOUIZx34LGD25VSKliawEMoPbke9+19go4xG9ln4gEDCGkpiW6HppSKAjqEEkIvHv8RA2MXsbLseOpLMU3YRWJ8LCMHdXQ7NKVUFNAEHiorp3Pi98+z/rgc3qh3GQAnJ+3lkZwuegFTKeUITeChsGU5TP8dtO7J8deO5eHhgwF4eViaJm+llGM0gTttz6/w5lVQPwUufxXi6kGSL2nvzHM3NqVUVNGLmE4qLYbJw2H3VrjhfWjc0m5vdAxILOzc5G58SqmoogncSR+Mhp8+haHjIP2UQ9tjYqFxK03gSilH6RCKUxa9DF+/AL1ug5MvP/L5pDQo3Bj2sJRS0UsTuBPWfwmz/gQZ/aH/g1Xvk5yuZ+BKKUdpAg9W4UZ46xpIOQ6GvWiHS6qS5EvgxoQ3PqVU1NIEHowDe2HiVVCyH66cCIlNqt83KQ1KiqBoR/jiU0pFNU3ggTIGZt4Gm5fZM+/UE46+v04lVEo5TBN4oD5/ClZMgX73wQmDat7/YALXcXCllDM0gQfi+w/howchcxic+Uf/XpOUZm91JopSyiE1zgMXkY7ApAqb2gH3Aa/4trcBfgIuM8ZE3QBv5Y46D/aKp//nI6BlF7j4ORDx740at9RiHqWUo2o8AzfGrDbGdDPGdANOAfYC04BRwFxjTAdgru9xVCnvqJNXUIQBdhX8Qvu5N7KPeLjiDUho4P+bxcTaJK4JXCnlkNoOofQD1hpj1gNDgAm+7ROAbCcD84LyjjoAMZTxTPxzpLONP5q7IOXY2r9hUhrs1CEUpZQzapvArwDe9N1vYYzZ7Lu/BWhR1QtE5CYRyRWR3Pz8/ADDdEfFzjn3xE2ib+w33F9yHR/sbBvYGyZpMY9Syjl+J3ARSQAuBiZXfs4YY7DtZo5gjBlnjMkyxmSlpqYGHKgbyjvnXBjzJbfEvcMrJQN4s7Rf4B11tJhHKeWg2pyBnwcsNsZs9T3eKiKtAHy325wOzm0jB3WkcbxhdPwbLC1rz99Krgmuo05SGhTv1WIepZQjapPAr+TQ8AnATGC47/5wYIZTQXlFdvd0/nvqBtLlV54pyaFFSuPgOuok61xwpZRz/FpOVkQaAgOAmytsfhR4S0RGAOuBy5wPz2XGkJX3OqSeyPjf3QsxQU6br1jM0zIz+PiUUnWaXwncGLMHaFZp26/YWSnRa93HsHW5ne8dbPKGQ8U8Wk6vlHKAVmIezRfPQqMW0NWhLxeNWoLEaAJXSjlCE3h1tiyHtfPgtJttX0snxMbZJK5j4EopB2gCr84Xz0J8Q8i6wdn3TUrTM3CllCM0gVelcCOseBtOGX70Nb4DkZwOhZrAlVLB0wRelQXP22Kb03/n/HtrMY9SyiGawCvbVwiLJkDnobZNmtOS0qB4jz2OUkoFQRN4ZYtehgO7bHf5UNDOPEoph2gCr6jkACwYC237QFq30BxDO/MopRyiCbyiFVNg1ybodUfojqHFPEoph2gCL2eMnTp4TCfICGGBaeOWgOhMFKVU0DSBl1szF7Z9a8e+/W2TFojYeO3Mo5RyhCbwcl88DY1bQeYloT+WFvMopRygCRxg01L48RM47RaISwj98TSBK6UcoAkc7Nh3QmPIuj48x0tqbcfAtZhHKRUETeAFP8PKabZsvn5yeI5ZXsyzf2d4jqeUikqawBc8by9ahqJsvjoHpxLqhUylVODqdgIv2mHL5jOHQXLr8B23/Fg6lVApFYS6ncBz/2uHMkJVNl8dLeZRSjmg7ibwkv2wcCy0OwdadgnvsRu3AkSHUJRSQam7CXzZW7B7K/S+PfzHjo23rdp2bgz/sZVSUaNuJvCyMjt1sEUXewbuhqQ0PQNXSgWlbibwNXPgl9WhL5s/Gk3gSqkg+ZXARSRFRKaIyHciskpEzhCRpiIyR0R+8N063HsshD5/xi7rmpnjXgzJrXUWilIqKP6egT8NfGCMORE4GVgFjALmGmM6AHN9j70vbxGs/8zO+46Ndy+OpDTbOGKfFvMopQJTYwIXkWSgD/ASgDHmgDGmABgCTPDtNgHIDlWQjvriWaiXBD2GuxuHNnZQSgXJnzPwtkA+8F8RWSIiL4pIQ6CFMWazb58tQIuqXiwiN4lIrojk5ufnOxN1oHb8BN/OgFOug/pJ7sZyMIHrTBSlVGD8SeBxQA/geWNMd2APlYZLjDEGqHJlJmPMOGNMljEmKzU1Ndh4g/Plv0Fiw1s2Xx0tp1dKBcmfBL4R2GiMWeh7PAWb0LeKSCsA3+220ITokL3bYcmr0OXSQ8nTTVrMo5QKUo0J3BizBdggIh19m/oB3wIzgfKB5OHAjJBE6JTcl6B4L/S61e1IrLgEaHSMltMrpQIW5+d+twGvi0gCsA64Hpv83xKREcB64LLQhOiA4n2wcBxk9IcWnd2O5pCkNJ1KqJQKmF8J3BizFMiq4qkQdv910LKJsGcb9HKhbP5oktLh17VuR6GUilDRX4lZuBHm/g3SekDbPm5Hc7ikdB0DV0oFLLoTeMkBmHy9XXkwZ5x7ZfPVSUqD/YWwf5fbkSilIlB0J/A598HGr+DiZ6F5B7ejOZIW8yilghC9CXzlNFj4vO007+aaJ0ejjR2UUkGIzgT+yw8w41ZofSoMeMjtaKqX7DsD15koSqkARF8CP7AHJl0DcfXg0pftfGuvatzK3uoQilIqAP7OA48MxsC7d0H+d3D12+FtVByIuHrQMFWHUJRSAYmuBL7oZTvnu+9oyLBT1KcvyWPM7NVsKigiLSWRkYM6kt093d04K0pK1wSulApI9CTwTUvg/XugfT/ocw9gk/foqcspKi4FIK+giNFTlwN4J4knpcOOH92OQikVgaJjDLxoB7x1LTQ8BnJegBj71xoze/XB5H1w1+JSxsxe7UaUVUtK0zNwpVRAIj+Bl5XBtFtg52Z70bJhs4NPbSooqvIl1W13RXI67CuE/bvdjkQpFWEiP4F//hR8/wEMehiOPfWwp9JSEqt8SXXbXaHFPEqpAEV2Av/xU5j3EHTOgZ43HfH0yEEdSYyPPWxbYnwsIwd1PGJf12gxj1IqQJ6/iFntLJJdW2DKDdAsAy5+psp1TsovVHp+FgroGbhSqtY8ncCrm0UiZSUM+eYWOLAbhs+Eeo2rfY/s7uneStiVHSzm0TNwpVTteDqBVzeLZPd790HpF3bGyTEnuRSdQ+LrQ4PmmsCVUrXm6THwqmaLDIz5mt+UToesEdDVu02AaiVZ1wVXStWepxN45dkix8lW/i/+P6ySDBj8iEtRhUBSui5opZSqNU8n8IqzSOpxgOfjn6IMYUP/5+06ItFCi3mUUgHw9Bh4xVkkt+8eR+eY9Xx52vMM7N3T5cgclpQO+wrsSooJDd2ORikVITx9Bg42iX8+aBOXx82Hs/7EGedd5XZIztOphEqpAPiVwEXkJxFZLiJLRSTXt62piMwRkR98t01CEqExttKybR84596QHMJ1WsyjlApAbc7AzzHGdDPGZPkejwLmGmM6AHN9j50nApe+Ape/DjGxNe8fiZL1DFwpVXvBDKEMASb47k8AsoMPpxoxMVA/KWRv77rGvjNwnYmilKoFfxO4AT4UkUUiUr7oSAtjzGbf/S1Ai6peKCI3iUiuiOTm5+cHGW6Uiq8PDZrpEIpSqlb8nYVypjEmT0SOAeaIyHcVnzTGGBExVb3QGDMOGAeQlZVV5T4KX2ceHUJRSvnPrzNwY0ye73YbMA3oCWwVkVYAvtttoQqyTtDWakqpWqoxgYtIQxFpXH4fGAisAGYCw327DQdmhCrIOkGLeZRSteTPEEoLYJrY5VrjgDeMMR+IyNfAWyIyAlgPRMnCJC5JSrOt4Q7shYQGbkejlIoANSZwY8w64OQqtv8K9AtFUHVScmt7u2szNGvvbixKqYjg+UrMOqO8mKdwo7txKKUihiZwr9ByeqVULWkC9wotp1dK1ZImcK+IT4TEpprAlVJ+0wTuJVrMo5SqBU3gXpKsxTxKKf9pAveSpDRd0Eop5TdN4F6SlAZF26H4yGbOSilVmSZwL0nyFfPoOLhSyg+awL1EpxIqpWpBE7iXaDGPUqoWNIF7iZ6BK6VqQRO4lyQ0gMQmOhNFKeUXTeBeo8U8Sik/aQL3Gu3Mo5Tyk789MVW4JKVB3iK3o3DM9CV5jJm9mk0FRaSlJDJyUEeyu6e7HZZSUUETuNckpcPeX6B4n+1WH8GmL8lj9NTlFBWXApBXUMToqcsBNIkr5QAdQvGaZF9i2xX54+BjZq+mqLiU+uynGYUAFBWXMmb2apcjUyo6aAL3moNTCSM/gW8qKKJXzAo+qjeSD+qNIpbSg9uVUsHTBO415cU8kT6V8MAexjR8jTcS/kESe0mVQrrJGgDSUhJdDk6p6KAJ3GuioZjn5wUw9kwuKX2PV8rOY8D+xyk1Qp/YZSTGxzJyUEe3I1QqKvidwEUkVkSWiMi7vsdtRWShiKwRkUkikhC6MOuQhIZQPyUyh1CK98GH/w/GD4ayEhj+LklDnyAuJZ0lpgMD4pfzSE4XvYCplENqMwvlDmAVkOR7/BjwpDFmooiMBUYAzzscX90UiXPB8xbD9N9B/ndwyvUw8CGo15hsfDNO/ncpfPwPOp1Qz+1IlYoafp2Bi0hr4ALgRd9jAc4Fpvh2mQBkhyLAOikpLXISeMkBmPcwvNgf9u2Eq9+Gi56Ceo0P3699P8DA2o9dCVOpaOTvEMpTwD1Ame9xM6DAGFPie7wRqPJ7sYjcJCK5IpKbn58fVLB1RnKElNNvXQkvngufPA5dL4PffwkZ/aveN62bbdq85qPwxqhUFKsxgYvIhcA2Y0xA5YHGmHHGmCxjTFZqamogb1H3JKXDnnwo2e92JFUrLYFPn4D/nA27tsAVb8DQsZCYUv1rYmKh/bmwdi6UlVW/n1LKb/6MgfcGLhaR84H62DHwp4EUEYnznYW3BiLkO38EqDgXvGlbd2OpLP97mH6LLffvlA0X/BMaNvPvtRn9YcUU2LocWp0c2jg9QpcSUKFU4xm4MWa0Maa1MaYNcAUwzxjzG+Bj4BLfbsOBGSGLsq7xYmOHsjL48l/wn7Ng+zq4ZDxcNsH/5A32DBzqzDBK+VICeQVFGA4tJTB9iZ7rKGcEMw/8z8BdIrIGOyb+kjMhKc8l8IKfYcKFMPteaNcXfr8AMofV/n0at4CWXWHNXKcj9KTypQQq0qUElJNqtZiVMWY+MN93fx3Q0/mQ1KEhlI3uxgFQtANeyYbd22DIv6HbVSAS+Ptl9IMvnoV9hVA/2bk4Pai6JQN0KQHlFK3E9KJ6jWxyc/sMvLSEbeOvpHj7ei7ZdRe9Z7di+tIgY8rob4t8fvzEmRg9rLolA3QpAeUUTeBe5YHOPGtfu51j8r/k3uIbyDUdnRnDbd0TEhrXiXHwkYM6khgfe9g2XUpAOUkTuFclpUGhi0Moi16m/Y+v82LJeUwu7Xtwc9BjuHEJ0O5sWDMPjAk+Tg/L7p7OIzldSE9JRID0lERdSkA5Shs6eFVSOmxe5s6xf/ocZt3N/0q78kjJVUc8HfQYbkY/+O5d+OUHSD0huPfyuOzu6ZqwVcjoGbhXJaXDnm22VD2cdvwEb10DTdryjwb3UErsEbsEPYbbvp+9rQPDKEqFkiZwryqfiRLOzjz7d8GbV9qLjFdO5HeDe4RmDLfJ8dCsgyZwpYKkQyhelVxhLniTNgG/jd+VgGVlMPVmyF8NV0+B5hlkN7dPhaSSMKM/LPovFBdBvM7KUCoQmsC9yoFinlo1Ff7477B6Fgx+7FDFJCEcw83oDwuft+PtHapZAEspdVQ6hOJV5UMoQcxE8bsScNlkuzhVj2vhtJsDPl6ttOkNcfXt4lZKqYBoAveqeo2hXnDFPH5VAuYtgpm3wnG94PwngquyrI34RDi+t46DKxUETeBeFmRjhxorAXduhjevgkbHwOWv2jna4ZTRH375HnasD+9xw237j/Djp25HoaKQJnAvCzKBH7USsLgIJl5lZ55c8SY0bB5stLWX4ZtOGK3DKHu3wwf3wnOnwoSL3C3MUlFJE7iXBdmZp9pKwG5pMONW2LQEhr0ALTOdi7k2mp8AycdG3+qEJfvt0rvPdIcF/4aTLgQMrJzudmQqyugsFC9LSrerAJYcCHh4o8pZJJ8+YRsrnPv/4MQLHAg0QCL2LHz521BaDLHx7sXiBGPg2xnw0f22IKrdOTDw7/YX5PZ1sOJt6HWr21GqKKJn4F6WlAYY2L3Fuff87j2Y+xBkXgJn3e3c+wYqoz8c2AUbvnI7kuBs+ApeGgiTh0N8A9vc+drph77dZA6DTYvteLhSDtEE7mXlc8ELHergsnUlTL3RNhge8lz4ZpwcTds+EBPn2dko05fk0fvRebQdNYvej847ciXG7T/C5OvgpQFQsB4uegZu+ezI5s6dh9rblVPDEreqGzSBe9nBYh4HEvieX+DNKyChkW1C7JXqx/rJcOxpnkzgR22JVrQDZv8F/tUTVn8AZ/8ZblsMpwy3DZwrSznOLqW7YlrY/x4qemkC97KKzY2DUXIA3roWdm21ybv8fb2i/bmwZZmNz0OqKoQqKd7P+ln/B093sxcqu14Gty+Gc+61jTiOJnOYbeic/30Io1Z1iSZwL6ufZJsfBHMGvne7Td7rP4ch/4LWpzgXn1PKhxvWznM3jkoOL4QyDI75ijkJI7mjZLwdhrrlU/uZ+vsLsdMQQHQYRTlGE7jXJacHnsB//ASe91U7nvc4dL3U2dic0rIrNEz13DBKecFTe8ljcsKDjE14igPEcXf8X+Ga6dCyS+3eMKkVtDnTzkaJ8mYWKjw0gXtdUlrth1BKi+GjB2HCxZDQAH77UfjWOAlETIxdI3ztPCgrrXn/MBk5qCNd4/N4K+FvtJPNjC4eQY4Zw1nnB9HYufNQW326dYWzwao6qcYELiL1ReQrEflGRFaKyIO+7W1FZKGIrBGRSSIS5jrsOiIprXazULavg/GD4LN/Qver4ab/2a/7XpfRH4q2w+albkdyUHarHUxOfJhSSeCSAw/wSeML+XtOt+BWZ+w0BCQWVugwigqeP4U8+4FzjTG7RSQe+ExE3gfuAp40xkwUkbHACOD5EMZaNyW1ht1b/St0+WYizLrbzoK49OVDU9ciQftzALFVmekeGKffshwmXEy9+g055uZ3+LhZe2fet2Fz2xN0xdvQ7z5vTOVUEavGM3Bj7fY9jPf9McC5wBTf9glAdkgirOvKi3l2HaWYZ18hvH0jTLvZjiff8nlkJW+wiS2tmzfGwTcvs2uXxCfCde+CU8m7XOYwO2d802Jn31fVOX6NgYtIrIgsBbYBc4C1QIExpsS3y0ZAO7eGQk1zwTd8BWPPsmd05/zVJpyUY8MXn5My+sPGr+0ca7ds/gZeuRjiG9rPsmk7549x4gUQE6/DKCpofiVwY0ypMaYb0BroCZzo7wFE5CYRyRWR3Pz8/ADDrMOSq0ngZaXwvzEwfjBg4IYP4OyRVReRRIqM/mDKYN18d46/aanvwm+j0CVvgMQm9u+6cpptZadUgGo1C8UYUwB8DJwBpIhI+Rh6a6DKU0RjzDhjTJYxJis1NTWoYOukqop5CjbAyxfaNmiZObZ0+9ie7sTnpPQs28TCjWGUTUvglSFQL8mXvNuG9niZOfaX8oaFoT2Oimr+zEJJFZEU3/1EYACwCpvIL/HtNhyYEaog67R6SfaMsHwmysrpMLa3rVwcOg6GvWjL0aNBbBy07wtr5oV3nnTe4sOTdxBNpP3W8TzbUk6LelQQ/DkDbwV8LCLLgK+BOcaYd4E/A3eJyBqgGfBS6MKsu6Yv3cSPxSnM/+IL3nlomF3trlmGrQI8+XK3w3NeRn/YtQm2rQrP8fIWwSvZ9pfg9bOgyfHhOW69xnDCIPsL2UNz31VkqXEaoTFmGdC9iu3rsOPhKkTKF1P6D03oG/sNZSXCWDOUVj0eYEjTNm6HFxrtfV161nwELTqF9lgbF8GrQyExxXfx97jQHq+yzjl2/fCfPrNTC6NRcRF88ybkr4aBD9tvWcoxWonpYeWLKS02Hfi5LJWriv/Cowcu5fE569wOLXSS0yH1pNCPg2/MhVezoUETuG5W+JM3QIeBdnhsxdvhP3ao7fkV5j8GT2bCu3+EhWPhx/luRxV1NIF7WPliSk+VXEKfA0+zoKzTYdujVkY/+PlL2L+75n0DseFre+bdoJkvebs07TKhgR0LXzXTFmpFg+0/wqw/wZOdYf4/IL0HXD3VXl/QpXQdpwncw2rsKh+tMvpD6QE7tOC0DV8dnryTWzt/jNrIHGbnvbs1ddIpGxfZVS+f7QGLXrZ/r98vgN9Mtr+QT7wAVr1j+4Uqx2gC97CjdpWPZsedYduSOd2t/ueF8GoONEqF6987NMfeTe3PtVMnI3EYpazMNrP47/nw4rmwdj70uh3uXA7Z/4JjTjq0b+Yw2F/ouSWDI51eUfCw8kWTxsxezaaCItJSEhk5qGNwiylFgvj60OYsZ8fBf14Arw2DRi3sBUuvNLWIqwcnXWSHUYr32b+715Xsh2WT4Ivn4JfVdr2eQf+AHtfa2TVVadfXFjCteNsOGylHaAL3uCq7ytcFGf3gh9nw69rg1yJZ/T68/Vto3BKGv2vX5faSzKGw9DX7C+ukC92OpnpFOyB3PCz8j11grUUXyHnBrrtTzUJr05fkHTwBeaphFud/O4v4i/ba8X8VNB1CUd7kRJeeX9fCG5fbXqApx3kzeQO0PduOyXu1qGfPr/DBaPhnZ5j7N2jR2Ta0uOVT21LuKMm7Yk/RSUWnEl+6l6/mTAxv/FFMz8CVNzVtZysi13wEPW+s3WsP7IFPn4AvnoXYBBjwEJx2C8R5dMn62Hi7Tvg3E23sCQ3djuiQzctg4lWwa7Mdx+51m9+diCr3FF1Q1ol8k8zuRW/BBTeEKuI6Rc/AlTeJ2LPwHz/xf+aCMXaM9blTbQLvPBRuzYXet3s3eZfrnAPFe+H72W5HcsjKaZS8OJBthXu5uOgBen9/BdM3N/X75ZWnu5YRw3ulPelVmgv7dzkdbZ2kCVx5V0Z/m9R+XlDzvltX2jW8p9wADZrCDbMhZ5w3h0yqcnwvaNTSG7NRyspg3t9h8nUsLzmWC/Y9xDLTjryCIkZPXc70Jf51iKpquus7pWdQX4rt7BUVNE3gyrvanEWZxPH66+NpO2oWvR+dd2TyKCqA9+6xa6JvXQEX/NO2kTvudHdiDlRMLHTOhh/mwL6d7sWxbydM+g18MoZ3Yvtx+f6/kE/KwaeLiksZM3u1X29V1TTYb+NOoqh+C2/8oooCmsCVZ03/tpCvyjrSo3gRBg4/Aywrg0UTbOHI1y/AKdfBbYvh1BGRuyZ65jAo3Q+r33Pn+NvXwUsD7DDOeY9z+54bOMCRFyj9rQTO7p7OIzldSE9JRID0lET+kXMyid0vtdc23GzcESX0IqbyrDGzV3NhSVdGx79JC7azlaYUFZcy6/13yP56km1JduzpcP4YaNXV7XCD1/pUSD7Wduo5+YrwHnvtxzD5Onvt4Zqp0K4vaR/PI6+KZF2bSuAqp8FuzIEvn4PvZtnG2ypgegauPGtTQRH/KzsZgD6xy2hOIY/H/YcXDvzZNrjIecF2IoqG5A02eXYeaitQ924PzzGNgS//Da/lQONWcOPHtuiGEFYCp/eAlOO1pZwDNIErz0pLSeQ7cyxbTQq3xL7DvHp3kR37Ga/FDoXbcu0c5Gjr6p6ZA2Uldt2QUCvZDzP+ALNHQ8fz4bdzDutEVNUQyCM5XYIvLBOxf89182HPL8G9Vx2nQyjKs0YO6sjoqcv5uLQbV8TNZ37pyTzGddx84aDqS7YrqVgJGBFLEbTqZufAr5wKpwwP3XF2bYFJV9sm0mf/Gc4eBTFHns+FrBI4cxh89qRdQiBL54QHShO48qzyxDH2gxt4a2dftiZ1ZeTgE/1OKOWVgOXFJOUXQSu+t+eI2OT26ROwexs0Osb5Y+QtgolXw74CuOwVW0QUbi0yoVkHO4yiCTxgmsCVpwVzBli5EhAOTYPzbAIHW9TzyRjbrae2VaiVVP4G8nSn1WR9cz80bgEj5kDLTIeCrqXyX1T/e8x+G2jc0p04IpyOgauoVd10N883xGjRyXYlCvIiX8W1SIQyrt39IlmLR5Hf5GS4cb57ybtcZg5gbF9QFRBN4CpqRXRDjMwc25Wo0L+qx6qUfwNpyk7Gx4/h5rhZTCgZwLBdI6FhMweDDVBqRzuU4tVFvCKAJnAVtSK6IUZn39npt4GfnZYUbOKvca/yWb076BWzglHFv+X+kuvZUOih9m2dh8KGhVCwwe1IIpImcBW1QjYNLhyaZ0DLroENo+xYD+/exaf17+C62Nm8X9aT8w48ysTScwGPfQPJzLG3K7VfZiBqvIgpIscCrwAtAAOMM8Y8LSJNgUlAG+An4DJjjNbGKk+J6IYYmcPgo/thx092ad2a/LLGTs1bNhEQ8trkcNPaM/mhuPnBXTz3DT2KM7oAAAv0SURBVKRpO0jrbtdG6X2729FEHH/OwEuAu40xnYDTgT+ISCdgFDDXGNMBmOt7rJRySueh9rams9OtK+0qjP86FVZMgVNvhDu+oe11L/CHnP7e/waSOQw2L7UNOFStiDGmdi8QmQE85/vT1xizWURaAfONMUf91Z6VlWVyc3MDDlapumb70334pWAng4oePrIQKW+xnS/+3buQ0AhO/S2c8YfQzB0PpcKN8GRnOPev0Gek29F4kogsMsZkVd5eq3ngItIG6A4sBFoYYzb7ntqCHWKp6jU3ATcBHHfccbU5nFJ12vQleaz69WRGx7xCW9nEuoI0Rk9dTtNfF9Nny8t2Rb/6ybaK8rSb7TrokSi5tV2UbMU0TeC15PdFTBFpBLwN3GmMOWzBYmNP46s8lTfGjDPGZBljslJTU4MKVqm6ZMzs1Uw/0JMyI1wU8yW9YlYwngfp89nVsGkJ9Lsf7lwB54yO3ORdLjMHtq2Ebd+5HUlE8esMXETiscn7dWNM+WXxrSLSqsIQyrZQBalUXbSpoAhDU74yJ3Jb3DTipIwtpgl/K76G++58zFu9M4PVKRs+GGXnhB9zr9vRHBTUWjplZfYX7epZ8MOHtql2YkrNr6sFf2ahCPASsMoY888KT80EhgOP+m5nOBqZUnVcWkoieQVFjC8ZTFLcXl4v6ceU0j40T0nmvmhK3mBL+4/vbWej9B19cJVJNxcjC2gtneJ9to/r6lm2bdzuLSCxtmXenl/Cn8CB3sA1wHIRWerbdi82cb8lIiOA9cBljkamVB1Xvhrjh8Wn8uGBUwEPTgN0UuYwePdO2LIcWnV1fTEyv9fS2fMr/DDbdlJaMw+K99iLyhn9oOMF0GFAyIa4akzgxpjPgOoWXe7nbDhKqXLlSSKilsMNxkkXw6y77Vl4q66uL0Z21LV0fl1rE/Z378GGBWDKbEOMky+3SbvtWRBXL+Qx6mqESnlYRBci1VbDZtD+HDsO3v8B1xcjKx/CAhDK6C5r6B+7mPPiF8OzG+1OLTLhrLttQ4y07kc0GAn1EJAmcKWUd3TOgRm/h7xFhyXQisK1FMDIQR15YOoiriibxYi490mVQkpMDDuanQpZt0HHwUetkA3HEJCuhaKU8o4TL4DYBFgx1d3FyMrKyJZP+KLRSEbFT2RlWRseiP8js8//nNRbP4TTb6lxeYOjDQE5Rc/AlVLekZgCGf1h5TSyB/4dcOEawNp58OF9sHU5DVp1g8teoG/bPvSt5duEYwhIE7hSylsyh9kLhBsWkN29V/iuAWxZDnPuh7VzIeU4GPaSHdKpoleoP8IxBKRDKEopbzlhMMQl2tko4VCYB9N+B2PPsv1CBz4Mt+ZCl0sCTt4QnvXo9QxcKeUt9RrBCYNsT9DBj0FsiNLUvkL47ClY8G87DbDXrXZGSWITR94+HNNANYErpbwnM8d2I/rpUzu10EklByB3vG2oXLQdulxmV0JscryzxyH000A1gSulvKfDQFvNuHKqcwncGHtWP/dB2L4O2vaBAQ9BWjdn3t8FOgaulPKe+ERbHPPtTHvGHKz1X8JLA2DycIirD7+ZAtfOjOjkDXoGrpTyqsxhsPwtWDcfThhY+9fnr7a/AFbNsDNMGreCi5+DbldBTGzNr48AmsCVUt7U/lzbsGLF2/4lcGNse7lvZ8CqmZDvW1u8dU8Y/Cj0GA4JDUIbc5hpAldKeVNcApx0EaycYZdpja9/5D7G2DW3V820iXv7OpAYuzRt1gg46UJISgt/7GGiCVwp5V2dc2DJa7Bmjk3mYBslbPzal7RnQuHPds3tdmdDr9vhxAuhUd3o/qUJXCnlXW3PhgbNYPlkOz/72xmw6h3YtdmumdL+XOg7CjqeF/lt5QKgCVwp5V2xcdBpiJ23/e0MW6HZoT+cNMQW+9RPcjtCV2kCV0p5W6/b7G3bs213m2hrJxcETeBKKW9r2g4ufNLtKDxJC3mUUipC6Rm4UipqudnVPhw0gSulopLbXe3DocYhFBEZLyLbRGRFhW1NRWSOiPzgu3Vm/UWllHJIOFqauc2fMfCXgcGVto0C5hpjOgBzfY+VUsoz3O5qHw41JnBjzCfA9kqbhwATfPcnANkOx6WUUkGprnVZuLrah0Ogs1BaGGM2++5vAVpUt6OI3CQiuSKSm5+fH+DhlFKqdlztah8mQU8jNMYYwBzl+XHGmCxjTFZqat1Yn0Ap5b7s7uk8ktOF9JREBEhPSeSRnC5RcwETAp+FslVEWhljNotIK2Cbk0EppZQTQt3SzG2BnoHPBIb77g8HZjgTjlJKKX/5M43wTeBLoKOIbBSREcCjwAAR+QHo73uslFIqjGocQjHGXFnNU/0cjkUppVQt6FooSikVoTSBK6VUhBI7CzBMBxPJB9YH+PLmwC8OhuM0jS84Gl9wNL7geD2+440xR8zDDmsCD4aI5BpjstyOozoaX3A0vuBofMHxenzV0SEUpZSKUJrAlVIqQkVSAh/ndgA10PiCo/EFR+MLjtfjq1LEjIErpZQ6XCSdgSullKpAE7hSSkUozyVwERksIqtFZI2IHNHpR0Tqicgk3/MLRaRNGGM7VkQ+FpFvRWSliNxRxT59RaRQRJb6/twXrvh8x/9JRJb7jp1bxfMiIs/4Pr9lItIjjLF1rPC5LBWRnSJyZ6V9wvr5BdMyUESG+/b5QUSGV7VPiOIbIyLf+f79polISjWvPerPQgjje0BE8ir8G55fzWuP+n89hPFNqhDbTyKytJrXhvzzC5oxxjN/gFhgLdAOSAC+ATpV2uf3wFjf/SuASWGMrxXQw3e/MfB9FfH1Bd518TP8CWh+lOfPB94HBDgdWOjiv/UWbIGCa58f0AfoAayosO1xYJTv/ijgsSpe1xRY57tt4rvfJEzxDQTifPcfqyo+f34WQhjfA8Cf/Pj3P+r/9VDFV+n5J4D73Pr8gv3jtTPwnsAaY8w6Y8wBYCK2fVtFFdu5TQH6iYiEIzhjzGZjzGLf/V3AKiDSFhseArxirAVAim9N93DrB6w1xgRamesIE3jLwEHAHGPMdmPMDmAOR/aODUl8xpgPjTElvocLgNZOH9df1Xx+/vDn/3rQjhafL29cBrzp9HHDxWsJPB3YUOHxRo5MkAf38f0QFwLNwhJdBb6hm+7AwiqePkNEvhGR90Wkc1gDs92RPhSRRSJyUxXP+/MZh8MVVP8fx83PD/xrGeiVz/EG7DeqqtT0sxBKt/qGeMZXMwTlhc/vLGCrMeaHap538/Pzi9cSeEQQkUbA28CdxpidlZ5ejB0WOBl4Fpge5vDONMb0AM4D/iAifcJ8/BqJSAJwMTC5iqfd/vwOY+x3aU/OtRWRvwAlwOvV7OLWz8LzQHugG7AZO0zhRVdy9LNvz/9f8loCzwOOrfC4tW9blfuISByQDPwalujsMeOxyft1Y8zUys8bY3YaY3b77r8HxItI83DFZ4zJ891uA6Zhv6pW5M9nHGrnAYuNMVsrP+H25+eztXxYSapvGejq5ygi1wEXAr/x/ZI5gh8/CyFhjNlqjCk1xpQBL1RzXLc/vzggB5hU3T5ufX614bUE/jXQQUTa+s7SrsC2b6uoYju3S4B51f0AO803ZvYSsMoY889q9mlZPiYvIj2xn3FYfsGISEMRaVx+H3uxa0Wl3WYC1/pmo5wOFFYYLgiXas983Pz8KvCnZeBsYKCINPENEQz0bQs5ERkM3ANcbIzZW80+/vwshCq+itdUhlZzXH/+r4dSf+A7Y8zGqp508/OrFbevolb+g50l8T32CvVffNv+hv1hBaiP/eq9BvgKaBfG2M7Efp1eBiz1/TkfuAW4xbfPrcBK7FX1BUCvMMbXznfcb3wxlH9+FeMT4F++z3c5kBXmf9+G2IScXGGba58f9hfJZqAYOw47AntNZS7wA/AR0NS3bxbwYoXX3uD7OVwDXB/G+NZgx4/LfwbLZ2WlAe8d7WchTPG96vvZWoZNyq0qx+d7fMT/9XDE59v+cvnPXIV9w/75BftHS+mVUipCeW0IRSmllJ80gSulVITSBK6UUhFKE7hSSkUoTeBKKRWhNIErpVSE0gSulFIR6v8DbG2r/hb2J+4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "# 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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The end."
   ]
  },
  {
   "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
}