{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Suppose $f(t,x) = x_0 e^{x_1 t}$\n",
      "\n",
      "with data\n",
      "\n",
      "    t: 0     1    2     3\n",
      "    y: 2.0   0.7  0.3   0.1\n",
      "\n",
      "Use Gauss Newton to fit the data."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy as sp\n",
      "import matplotlib.pyplot as pt\n",
      "import scipy.linalg as la"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/lib/python2.7/dist-packages/numpy/oldnumeric/__init__.py:11: ModuleDeprecationWarning: The oldnumeric module will be dropped in Numpy 1.9\n",
        "  warnings.warn(_msg, ModuleDeprecationWarning)\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def residual(x):\n",
      "    return y - x[0] * np.exp(x[1] * t)\n",
      "\n",
      "def jacobian(x):\n",
      "    return np.array([\n",
      "        ...,\n",
      "        ...\n",
      "        ]).T"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<!--\n",
      "def jacobian(x):\n",
      "    return np.array([\n",
      "        -np.exp(x[1] * t),\n",
      "        -x[0] * t * np.exp(x[1] * t)\n",
      "        ]).T\n",
      "-->\n",
      "(Edit this cell for solution)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# data\n",
      "t = np.array([0.0, 1.0, 2.0, 3.0])\n",
      "y = np.array([2.0, 0.7, 0.3, 0.1])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# initial guess\n",
      "x = np.array([1, 0])\n",
      "#x = np.array([0.4, 2])\n",
      "\n",
      "def plot_iterate(x):\n",
      "    pt.plot(t, y, 'ro', markersize=20, clip_on=False)\n",
      "    T = np.linspace(t.min(), t.max(), 100)\n",
      "    Y = x[0] * np.exp(x[1] * T)\n",
      "    pt.plot(T, Y, 'b-')\n",
      "\n",
      "plot_iterate(x)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEFCAYAAAAG45eHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFWpJREFUeJzt3X9MVGe+x/HPMMAfVGI1C3QXyW2j2DJKEW43c5P1B6bj\nmrItmqxNJFk01U2oSbdpQrLbbtKrtF5qE7FpY2NsbtIEmoKrfxR0lciYRTe4yG5lN9BpVmyWW6Qt\nQjZuaLObOR7O/UNhRWBmYA4y8/B+JSdhOM88833ytJ+ePsx5jsdxHEcAACOlLHQBAID5Q8gDgMEI\neQAwGCEPAAYj5AHAYIQ8ABgsdb46tixLn9TXq/fCBTk3byrFtjXm9cqTna3CQEDbKiuVlpY2Xx8P\nAJDkcft78uFwWPWHD6uvqUk7QiGV2La895y3JV31enXS59Pqigrtqq5Wenq6myUAAO5yNeRHR0f1\n62ef1YsdHVpj21Hb93q9Or5+vWpPn1ZmZqZbZQAA7nIt5EdHR/Xqli2quXJF35vF+4YlHfD7dait\njaAHAJe5EvLhcFjVW7Zo/6VLswr4ccOS3ti4UXVtbSzdAICLIn67ZmBgQJs3b9aaNWu0du1avffe\ne9O2e6a0VC2XLulpSd1zKCJLUlVHh+rr6ubwbgDATCKGfFpamt555x199tln6uzs1Pvvv6/PP/98\nUpuWlhb19/To/yR9IGnfHAtZa9u61tgoy7Lm2AMA4H4RQ/6RRx7RunXrJElLlixRQUGBvvrqq0lt\njh45or3ffSdJ8ku6JWlojsU8HwqpuaFhju8GANwv5puh+vv71d3dLb/fP/n3fX3acM+y/gpJN+ZY\nTIltqycYnOO7AQD3iynkv/32W+3YsUPvvvuulixZMvlkODylE88ci/FKcoaH5/huAMAUThThcNj5\n8Y9/7LzzzjvTnv/P73/faZQc5+7xuOR8c8/r8WOl5IiDg4ODY1bHypUro8V0RBGv5B3H0d69e+Xz\n+fTKK69M2yY/O1v1d3/ulPSwpJxp2n0R44heDwTkOE7SHfv371/wGhgb42N85h1ffPFFpJiOKuLe\nNR0dHfroo4/05JNPqri4WJJUW1urL7/8UpJUVVWl1T6fRv7yF62S9JCkD+MoxpbkycqKowcAwL0i\nhvz69es1NjYWsYPCQED/8Zvf6P0YtjGI5qrXq8JAIO5+AAB3xL3V8LbKSp3y+dyoRSd9Pm2rrHSl\nrwettLR0oUuYNyaPTWJ8yc708cXLlW0N/vett/Rfr7+utRGu5j26s+Y+k16vV51vvqmfv/ZavOUA\ngDE8Ho/iiekHtndNpJAflvTGpk2qO3+evWsA4B7xhrwrT4ZKT09X7Zkz2u/3a2SW7x3fhbL29GkC\nHgBc9sD2k5/uSp795AEgsoRYrrlXOBxWw5Ej+uvHH+v5e54MNR7yPBkKAGKXcCE/zrIsNTc0qCcY\nlDM8rDeDQb0eCMiTlcUzXgEgRgkb8lM+KM5CAWAxSog/vAIAEhMhDwAGI+QBwGCEPAAYjJAHAIMR\n8gBgMEIeAAxGyAOAwQh5ADAYIQ8ABiPkAcBghDwAGIyQBwCDEfIAYDBCHgAMRsgDgMEIeQAwGCEP\nAAYj5AHAYIQ8ABiMkAcAgxHyAGAwQh4ADEbIA4DBCHkAMBghDwAGI+QBwGCEPAAYjJAHAIMR8gBg\nMEIeAAxGyAOAwQh5ADAYIQ8ABiPkAcBghDwAGIyQBwCDEfIAYDBCHgAMFjXk9+zZo5ycHBUWFk57\nvr29XUuXLlVxcbGKi4t18OBB14sEAMxNarQGL7zwgn7xi19o165dM7bZtGmTWlpaXC0MABC/qFfy\nGzZs0LJlyyK2cRzHtYIAAO6Je03e4/Ho8uXLKioqUllZmUKhkBt1AQBcEHW5JpqSkhINDAwoIyND\n586d0/bt23Xt2rVp2x44cGDi59LSUpWWlsb78QBglPb2drW3t7vWn8eJYa2lv79fzz33nHp6eqJ2\n+Nhjj+nTTz/V8uXLJ3+Qx8OyDgDMUrzZGfdyzdDQ0EQBXV1dchxnSsADABZG1OWaiooKXbx4USMj\nI8rLy1NNTY0sy5IkVVVV6dSpUzp27JhSU1OVkZGhpqameS8aABCbmJZrXPkglmsAYNYWfLkGAJC4\nCHkAMBghDwAGI+QBwGCEPAAYjJAHAIMR8gBgMEIeAAxGyAOAwQh5ADAYIQ8ABiPkAcBghDwAGIyQ\nBwCDEfIAYDBCHgAMRsgDgMEIeQAwGCEPAAYj5AHAYIQ8ABiMkAcAgxHyAGAwQh4ADEbIA4DBCHkA\nMBghDwAGI+QBwGCEPAAYjJAHAIMR8gBgMEIeAAxGyAOAwQh5ADAYIQ8ABiPkAcBghDwAGIyQBwCD\nEfIAYDBCHgAMRsgDgMEIeQAwGCEPAAYj5AHAYIQ8ABiMkAcAgxHyAGCwqCG/Z88e5eTkqLCwcMY2\nL7/8svLz81VUVKTu7m5XCwQAzF3UkH/hhRfU2to64/mzZ8/q+vXr6uvr0wcffKB9+/a5WiAAYO6i\nhvyGDRu0bNmyGc+3tLRo9+7dkiS/369bt25paGjIvQoBAHMW95r84OCg8vLyJl6vWLFCN27ciLdb\nAIALUt3oxHGcSa89Hs+07TyeA/e8Kr17AAD+rf3u4Y64Qz43N1cDAwMTr2/cuKHc3Nxp2zrOgXg/\nDgAMV6p7L4A9npq4eot7uaa8vFz19fWSpM7OTj388MPKycmJt1sAgAuiXslXVFTo4sWLGhkZUV5e\nnmpqamRZliSpqqpKZWVlOnv2rFatWqWHHnpIH3744bwXDQCIjce5f0F9vj7I45mydg8AiCze7OSO\nVwAwGCEPAAYj5AHAYIQ8ABiMkAcAgxHyAGAwQh4ADEbIA4DBCHkAMBghDwAGI+QBwGCEPAAYjJAH\nAIMR8gBgMEIeAAxGyAOAwQh5ADAYIQ8ABiPkAcBgUR/kjcXHsix9Ul+v3gsX5Ny8qRTb1pjXK092\ntgoDAW2rrFRaWtpClwkgBjzIGxPC4bDqDx9WX1OTdoRCKrFtee85b0u66vXqpM+n1RUV2lVdrfT0\n9IUqF1gU4s1OQh6SpNHRUf362Wf1YkeH1th21Pa9Xq+Or1+v2tOnlZmZ+QAqBBYnQh5xGx0d1atb\ntqjmyhV9bxbvG5Z0wO/XobY2gh6YJ4Q84hIOh1W9ZYv2X7o0q4AfNyzpjY0bVdfWxtINMA/izU6+\nXbPI1R8+rBc7OuYU8JKUJamqo0P1dXVulgXAJYT8ImZZlvqammJag49krW3rWmOjLMtyqTIAbiHk\nF7HmhgbtCIVc6ev5UEjNDQ2u9AXAPYT8ItYTDKokzqv4cSW2rZ5g0JW+ALiHkF/EnJs3J30PPh5e\nSc7wsEu9AXALIb+Ipbh0FT/R3+3brvYHIH6E/CI25nXrOv5uf6nskgEkGkJ+EfNkZ8uta3lbkicr\ny6XeALiFkF/ECgMBXXXpav6q16vCQMCVvgC4h5BfxLZVVuqUz+dKXyd9Pm2rrHSlLwDuIeQXsbS0\nNOVXVKg3zqv5Xq9Xqysq2H4YSEDsXbPIubJ3zaZNqjt/nr1rgHnA3jWIS3p6umrPnNF+v18js3zv\n+C6UtadPE/BAguJKHpLYTx5IVGw1DNeEw2E1HDmiv378sZ7nyVBAQiDk4TrLstTc0KCeYFDO8LBS\nbt/WWGqqPFlZPOMVeMAIeQAwGH94BQDMiJAHAIMR8gBgMEIeAAxGyAOAwQh5ADBY1JBvbW3VE088\nofz8fL399ttTzre3t2vp0qUqLi5WcXGxDh48OC+FAgBmL+KjfGzb1ksvvaRgMKjc3Fz98Ic/VHl5\nuQoKCia127Rpk1paWua1UADA7EW8ku/q6tKqVav06KOPKi0tTTt37lRzc/OUdtzkBACJKWLIDw4O\nKi8vb+L1ihUrNDg4OKmNx+PR5cuXVVRUpLKyMoVCofmpFAAwaxGXazweT9QOSkpKNDAwoIyMDJ07\nd07bt2/XtWvXpm174MCBiZ9LS0tVWlo6q2IBwHTt7e1qb293rb+Ie9d0dnbqwIEDam1tlSS99dZb\nSklJ0a9+9asZO3zsscf06aefavny5ZM/iL1rAGDW5nXvmqeeekp9fX3q7+9XOBzWiRMnVF5ePqnN\n0NDQRAFdXV1yHGdKwAMAFkbE5ZrU1FQdPXpUW7dulW3b2rt3rwoKCnT8+HFJUlVVlU6dOqVjx44p\nNTVVGRkZampqeiCFAwCiY6thAEhgbDUMAJgRIQ8ABiPkAcBghDwAGIyQBwCDEfIAYDBCHgAMFvFm\nKADJx7IsfVJfr94LF+TcvKkU29aY1ytPdrYKAwFtq6xUWlraQpeJB4SboQBDhMNh1R8+rL6mJu0I\nhVRi2/Lec96WdNXr1UmfT6srKrSrulrp6ekLVS5iFG92EvKAAUZHR/XrZ5/Vix0dWmPbUdv3er06\nvn69ak+fVmZm5gOoEHNFyAOL3OjoqF7dskU1V67oe7N437CkA36/DrW1EfQJjJAHFrFwOKzqLVu0\n/9KlWQX8uGFJb2zcqLq2NpZuEhR71wCLWP3hw3qxo2NOAS9JWZKqOjpUX1fnZllIIIQ8kKQsy1Jf\nU1NMa/CRrLVtXWtslGVZLlWGRELIA0mquaFBO1x6pvLzoZCaGxpc6QuJhZAHklRPMKiSOK/ix5XY\ntnqCQVf6QmIh5IEk5dy8Oel78PHwSnKGh13qDYmEkAeSVIpLV/ET/d2+7Wp/SAyEPJCkxrxuXcff\n7S+VXU5MRMgDScqTnS23ruVtSZ6sLJd6QyIh5IEkVRgI6KpLV/NXvV4VBgKu9IXEQsgDSWpbZaVO\n+Xyu9HXS59O2ykpX+kJiIeSBJJWWlqb8igr1xnk13+v1anVFBdsPG4q9a4Ak5sreNZs2qe78efau\nSVDsXQMsYunp6ao9c0b7/X6NzPK947tQ1p4+TcAbjCt5wADsJ28uthoGIOnO0k3DkSP668cf63me\nDGUMQh7AJJZlqbmhQT3BoJzhYaXcvq2x1FR5srJ4xmsSIuQBwGD84RUAMCM2qwCABGJZlj6pr1fv\nhQtybt6UJP13ICBPdvaclttYrgGABBAOh1V/+LD6mpq0454/nHskOZr7H84JeQBYYJG+Ajse8vea\nzVdgCXkAWECjo6N6dcsW1Vy5Mu1dy9OFvPTvm9kOtbVFDHpCHgAWSCzbUswU8tLdbSk2blRdW9uM\nSzd8uwYAFkj94cN6saNjTvsOSVKWpKqODtXX1c3YhpAHgAVgWZb6mppi2oYikrW2rWuNjbIsa9rz\nhDwALIDmhgbtCIVc6ev5UEjNDQ3TniPkAWAB9ASDKnHpYewltq2eYHDac4Q8ACwA5+ZNufUodq8k\nZ3h42nOEPAAsgBSXruIn+rt9e/rfu/opAICYjLn0EPaJ/lKn36WGkAeABeDJzpZb1/K2JE9W1rTn\nCHkAWACFgYCuunQ1f9XrVWEgMO05Qh4AFsC2ykqd8vlc6eukz6dtlZXTniPkAWABpKWlKb+iQr1x\nXs33er1aXVEx4/bDUUO+tbVVTzzxhPLz8/X2229P2+bll19Wfn6+ioqK1N3dHVfBALBY7Kqu1vEf\n/Ugjc3z/sKTj69drV3X1jG0ihrxt23rppZfU2tqqUCikxsZGff7555PanD17VtevX1dfX58++OAD\n7du3b47lJrf29vaFLmHemDw2ifElu2QeX3p6umrPnNF+v3/WQT++C2Xt6dMR95WPGPJdXV1atWqV\nHn30UaWlpWnnzp1qbm6e1KalpUW7d++WJPn9ft26dUtDQ0OzLDf5JfM/aNGYPDaJ8SW7ZB9fZmam\nDrW1qWbjRn0W49JNr9erNzZtirrNsBQl5AcHB5WXlzfxesWKFRocHIza5saNGzEVCgC4E/R1bW3q\nPHhQvyws1B+93ilfr7Ql/dHr1S8LC9X55puqO38+asBLUZ7x6vF4Yirw/n3iY30fAOCO9PR07X31\nVVnV1WpuaNCZYPDOVgXB4J1nvGZlqTAQ0P/M8hmvciL4wx/+4GzdunXidW1trXPo0KFJbaqqqpzG\nxsaJ148//rjzzTffTOlr5cqVju7sfc/BwcHBEeOxcuXKSDEdVcQr+aeeekp9fX3q7+/XD37wA504\ncUKNjY2T2pSXl+vo0aPauXOnOjs79fDDDysnJ2dKX9evX4/0UQCAeRAx5FNTU3X06FFt3bpVtm1r\n7969Kigo0PHjxyVJVVVVKisr09mzZ7Vq1So99NBD+vDDDx9I4QCA6B7YM14BAA+e63e8mnzzVLSx\ntbe3a+nSpSouLlZxcbEOHjy4AFXOzZ49e5STk6PCwsIZ2yTrvEnRx5fMcydJAwMD2rx5s9asWaO1\na9fqvffem7Zdss5hLONL1jn817/+Jb/fr3Xr1snn8+m1116btt2c5y6uFf373L5921m5cqXzt7/9\nzQmHw05RUZETCoUmtfntb3/rPPPMM47jOE5nZ6fj9/vdLGHexDK23/3ud85zzz23QBXG59KlS87V\nq1edtWvXTns+WedtXLTxJfPcOY7jfP311053d7fjOI4zOjrqrF692ph/9xwntvEl8xx+9913juM4\njmVZjt/vd37/+99POh/P3Ll6JW/yzVOxjE3SlK+TJosNGzZo2bJlM55P1nkbF218UvLOnSQ98sgj\nWrdunSRpyZIlKigo0FdffTWpTTLPYSzjk5J3DjMyMiRJ4XBYtm1r+fLlk87HM3euhrzJN0/FMjaP\nx6PLly+rqKhIZWVlCrn0kN5EkKzzFiuT5q6/v1/d3d3y+/2Tfm/KHM40vmSew7GxMa1bt045OTna\nvHmzfPftThnP3EX8ds1smXzzVCw1lpSUaGBgQBkZGTp37py2b9+ua9euPYDqHoxknLdYmTJ33377\nrXbs2KF3331XS5YsmXI+2ecw0viSeQ5TUlL05z//Wf/4xz+0detWtbe3q7S0dFKbuc6dq1fyubm5\nGhgYmHg9MDCgFStWRGxz48YN5ebmulnGvIhlbJmZmRP/2/XMM8/Isiz9/e9/f6B1zpdknbdYmTB3\nlmXppz/9qX72s59p+/btU84n+xxGG58Jc7h06VL95Cc/0Z/+9KdJv49n7lwN+XtvngqHwzpx4oTK\ny8sntSkvL1d9fb0kRbx5KtHEMrahoaGJ/9p2dXXJcZwpa2vJKlnnLVbJPneO42jv3r3y+Xx65ZVX\npm2TzHMYy/iSdQ5HRkZ069YtSdI///lPtbW1qbi4eFKbeObO1eUak2+eimVsp06d0rFjx5SamqqM\njAw1NTUtcNWxq6io0MWLFzUyMqK8vDzV1NTIsixJyT1v46KNL5nnTpI6Ojr00Ucf6cknn5wIiNra\nWn355ZeSkn8OYxlfss7h119/rd27d2tsbExjY2OqrKzU008/7VpucjMUABiMx/8BgMEIeQAwGCEP\nAAYj5AHAYIQ8ABiMkAcAgxHyAGAwQh4ADPb/GbZdVInFTjAAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x4906b90>"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# evaluate this cell in-place many times (Ctrl-Enter)\n",
      "\n",
      "J = jacobian(x)\n",
      "r = residual(x)\n",
      "s = la.lstsq(...)[0]\n",
      "x = x + s\n",
      "print la.norm(r)\n",
      "\n",
      "plot_iterate(x)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.0446775329872\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEFCAYAAAAG45eHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4FOX5//H3sklaAwioEJRzOQdCCGhTlUPQIEo5WVFJ\nLVBADbbIj5qrin6rgGjUlmg91S/6bamBCihVAYVIggY0GFAONhCUQI2GgBC0WKDaLJv5/fEYJCTZ\nbLKbzB4+r+uaiyQzmb3nmvbO4z3P3I/DsiwLEREJSc3sDkBERBqPkryISAhTkhcRCWFK8iIiIUxJ\nXkQkhCnJi4iEsIjGOrHL5eL1zEx2b9yIdfQozdxuKpxOHO3aEZeczPjJk4mMjGysjxcREcDh73ny\n5eXlZC5aRNGKFUwsLGSQ243zrP1uYIfTySuxsfRKSWFKWhpRUVH+DEFERL7j1yR/4sQJ7hszhpl5\nefRzu+s8frfTyeIhQ0hfu5aWLVv6KwwREfmO35L8iRMnmDtyJAu2buWievxeGTA/MZFHs7OV6EVE\n/MwvSb68vJy0kSOZt3lzvRJ8pTLgwWHDyMjOVulGRMSPPM6uKSkpYcSIEfTr14/+/fvz1FNP1Xjc\ndUlJrNm8mauBnQ0Ioi2QmpdHZkZGA35bRERq4zHJR0ZG8sQTT7Bnzx7y8/N59tln2bt3b5Vj1qxZ\nQ3FBAZ8BzwN3NDCQ/m43+5Yvx+VyNfAMIiJyLo9Jvn379gwcOBCAFi1a0LdvXw4dOlTlmGcef5wZ\np04BkAgcB440MJgbCwtZvXRpA39bRETO5fXLUMXFxezcuZPExMSqPy8qYuhZZf2OwMEGBjPI7aYg\nJ6eBvy0iIufyKsmfPHmSiRMn8uSTT9KiRYuqO8vLq53E0cBgnIBVVtbA3xYRkWqsOpSXl1vXXHON\n9cQTT9S4f/DFF1vLwbK+23qD9cVZ31du53OhBWjTpk2btnps3bt3rytNe+RxJG9ZFjNmzCA2NpY5\nc+bUeEzPdu3I/O7rfKA1EFPDcf/mSw7Rvs4ruj85Gcuygm6bN2+e7THo2nR9ur7Q2w4cOOApTdfJ\nY++avLw8li1bxoABA0hISAAgPT2dzz//HIDU1FR6xcZy7KOP6AE0B5Z4ON+jzOVJav5jAablgaNt\n23pegoiI1MZjkh8yZAgVFRUeTxCXnEyXl1/mWS/aGCxlMnfzezpwqMb9O5xO4pKT6zyPiIh4x+dW\nw+MnT2ZVbKxXx05jCY8yt9b9r8TGMn7yZF9DskVSUpLdITSaUL420PUFu1C/Pl/5pa3B/z3yCD+5\n/376exjNO4AvaEdf9vIPBtCR0ir7dzud5C9cyK333utrOCIiIcPhcOBLmm6y3jUOzIPV3/J7TtGc\nP/HrM/vKgAeHDydjwwb1rhEROUtAJHmouwtlZZIv4yL68DHbGUxXPlMXShERDwImyYPnfvKVSR7g\nfh7kIB1Jc96mfvIiIh4EVJIHU7pZ+vjjfPLSS9x41spQlUneDWxqdgFjHEXce+cr3PPYNJVoRERq\nEXBJvpLL5WL10qUU5ORglZWxMCeH+5OTcbRtS1xyMh8fmkpBgZOVKxvj00VEQkPAJvlqH3ROoKdO\nQY8esH49fNfoUkREzuFrkvd5nnxDNW8Oc+fC/ffbFYGISOizbSQP8O230KsXLF8OV17ZFFGIiASX\noB3JA/zwh7BggRnRN82fGhGR8GJrkgeYMgW++grefNPuSEREQo/tSd7phPR0uPde8KLHmYiI1IPt\nSR5g3Dg4/3xYtszuSEREQoutD17P9t57cMst8MknplYvIiJB/uD1bEOGwIAB8Kc/2R2JiEjoCJiR\nPEBhIQwfDvv2QZs2TRGViEhgC9o3XmuTmgotW8KiRU0QlIhIgAu5JP/FF9CvH3z4IXTr1gSBiYgE\nsJCpyVdq3x7+3/+D++6zOxIRkeAXcCN5MM3LevWC116DH/+4kQMTEQlgITeSB9O87MEHIS1N7Q5E\nRHwRkEke4Je/hH//G1591e5IRESCV0CWayq9/TbcequZWqkXpEQkHIVkuabSVVeZF6SefNLuSERE\nglNAj+QBiorg8sth924z80ZEJJyE3Dz5mqSlmfr8Cy/4OSgRkQAXFkn++HHo3RuysiAhwc+BiYgE\nsJCuyVdq3dpMqZw9W1MqRUTqIyiSPJhZNqdOmfVgRUTEO0FRrqmUlwc33wwffwwtWvgpMBGRABYW\n5ZpKV14JI0bAww/bHYmISHAIqpE8wKFDZu78++9Dz55+CExEJICF1Uge4JJL4O67TadKPYQVEfEs\n6JI8wJw58OmnsHq13ZGIiAS2oCvXVHr7bZg2zfS1ad7cb6cVEQkoYVeuqXTVVXDFFXoIKyLiSdCO\n5OH7h7DvvQd9+vj11CIiASFsR/JgHsL+z//ArFl6CCsiUpOgTvIAd94Jx47BihV2RyIiEniCulxT\naetWuP562LMH2rRplI8QEbFFWHSh9MasWeByweLFjfYRIiJNTkn+O19/DbGx8PLLpv2BiEgoCOsH\nr2dr1Qr++EdITYXycrujEREJDCGT5AEmToQuXeAPf7A7EhGRwBAy5ZpKn30Ggwdr7ryIhAaVa87R\npQvMmwe33QYVFXZHIyJirzqT/PTp04mJiSEuLq7G/bm5ubRq1YqEhAQSEhJ46KGH/B5kff3qV+B2\na6aNiEid5Zp3332XFi1aMGXKFAoKCqrtz83N5fHHH2fNmjWeP6iJyjWVCgth+HDYuRM6dmyyjxUR\n8atGL9cMHTqUNnW8YdSUydtbsbHmbdiZM9XyQETCl881eYfDwZYtW4iPj2f06NEUFhb6Iy6/mDsX\nSkpg2TK7IxERsUeErycYNGgQJSUlREdHs379eiZMmMC+fftqPHb+/Plnvk5KSiIpKcnXj/coKgr+\n+lcYNQqSk+Hiixv140REfJabm0tubq7fzufVFMri4mLGjh1bY03+XN26dWP79u1ccMEFVT+oiWvy\nZ3vgAdi1y6wk5XDYEoKISIPYPoXyyJEjZwLYtm0blmVVS/B2+93voLgYXnrJ7khERJpWneWalJQU\nNm3axLFjx+jUqRMLFizA5XIBkJqayqpVq3juueeIiIggOjqaFQHY8zcqCpYsgdGjzYpSKtuISLgI\nuTdePXngAdi+Hd54Q2UbEQkOtpdrgsnvfgeHD8Of/2x3JCIiTSOsRvIAu3fDiBGwbRt062Z3NCIi\nnmkkX0/9+8Pdd8Mvf6neNiIS+sIuyQPcdZdJ8E88YXckIiKNK+zKNZX++U9ITISNG2HAALujERGp\nmco1DfSjH5nFRW65Bb791u5oREQaR9iO5ME0LrvpJujQwSwdKCISaLSQt4++/BIGDjTTKq+5xu5o\nRESqUrnGRxdeaN6GnTYNysrsjkZExL/CfiRf6Z57YM8eWLtWb8OKSODQSN5PFi6Eo0fh6aftjkRE\nxH80kj/LgQPwk5/Ahg2QkGB3NCIiGsn7Vffu8OSTMGkSnDxpdzQiIr7TSL4G06aZ6ZV//avdkYhI\nuNNIvhE884xpYKYkLyLBTiP5WuzZA0lJkJsL/frZHY2IhCuN5BtJv36waBHceCOcOmV3NCIiDaOR\nfB2mTQO3G158UfPnRaTpaSTfyJ59FnbtgueftzsSEZH600jeC0VFcOWV8OabcNlldkcjIuFEI/km\n0LMnLF5s6vPHjtkdjYiI9zSSr4e774aPPoJ168DptDsaEQkHGsk3ofR0KC+HefPsjkRExDtK8vUQ\nEQErV8LSpfDqq3ZHIyJSN5VrGmD7drj2Wr0oJSKNT+UaGwweDBkZcP31cPy43dGIiNROI3kfzJ5t\nple+8YYexIpI49BI3kYZGeZB7Ny5dkciIlIzJXkfREbCK6/Aa6+ZtgciIoFG5Ro/KCw0HSvXrDEr\nS4mI+IvKNQEgNhaWLIEbboDPP7c7GhGR7ynJ+8lPfwppaTBmDJw4YXc0IiKGyjV+ZFkwcyYcPAir\nV5uXp0REfKFyTQBxOMzSgeXlZlQvImI3JXk/q5xxs2EDPP203dGISLhTQaERtG5tOlUOGQKdO8P4\n8XZHJCLhSjX5RvThh3DddeaN2MREu6MRkWCkmnwAu/RSM7VywgQ4cMDuaEQkHCnJN7IxY0z/+Wuv\nhaNH7Y5GRMKNknwTmDkTJk0yc+lPnrQ7GhEJJ6rJNxHLgttvh5IS0/4gKsruiEQkGPiaO5Xkm9Dp\n0/Czn0GrVqahWTP9d5SI1EEPXoNIRASsWAGffgq/+Y0Z3YuINCYl+SYWHW2mVG7aBAsX2h2NiIQ6\nvQxlg9at4a23zMtSbdrAnXfaHZGIhColeZvExEB2Ngwdamr0U6bYHZGIhCIleRt17Wp63Fx1FZx3\nHtx4o90RiUioqbMmP336dGJiYoiLi6v1mNmzZ9OzZ0/i4+PZuXOnXwMMdX37QlYWzJplavUiIv5U\nZ5KfNm0aWVlZte5ft24d+/fvp6ioiOeff5477rjDrwGGg/h4WLsWpk+HnBy7oxGRUFJnkh86dCht\n2rSpdf+aNWuYOnUqAImJiRw/fpwjR474L8Iw8eMfw9//Dj//OeTm2h2NiIQKn6dQlpaW0qlTpzPf\nd+zYkYMHD/p62rA0dCi8/LKpzW/aZHc0IhIK/PLg9dy3sRwOR43HzZ8//8zXSUlJJCUl+ePjQ0pS\nknlhauJEePVVk/hFJHzk5uaS68f/nPeqrUFxcTFjx46loKCg2r6ZM2eSlJTEpEmTAOjTpw+bNm0i\nJiam6geprUG9ZGfDLbeYEo4SvUj4sr2twbhx48jMzAQgPz+f1q1bV0vwUn8jR8Lf/mZ63ahGLyIN\nVedIPiUlhU2bNnHs2DFiYmJYsGABLpcLgNTUVABmzZpFVlYWzZs3Z8mSJQwaNKj6B2kk3yBvvw03\n32xKOFdfbXc0ItLU1IUyDGzeDDfcAMuWwahRdkcjIk3J9nKNNL5hw2D1apg8GV5/3e5oRCSYqK1B\nkLjiCli/3qwudeqUeSgrIlIXJfkgMngwbNwI11xjEv3tt9sdkYgEOiX5INOvn3lRauRI+Ne/4J57\n7I5IRAKZHrwGqdJSM6IfPRp+/3uo5f0zEQlyml0Txr780tToY2Ph+efN8oIiElo0uyaMXXih6VpZ\nWmpemvrPf+yOSEQCjZJ8kGvRwrQpbtUKkpPN6F5EpJKSfAiIioIXXzQ9bq68EoqL7Y5IRAKFqrgh\nolkzeOwx6NDBJPo1a8yUSxEJb3rwGoJeew1SU+Evf4ExY+yORkR8oQevUs3115s6/e23w5/+ZHc0\nImInjeRD2D//aaZYjhoFGRngdNodkYjUl+bJi0fHj5vlBKOiYPlyOP98uyMSkfpQuUY8at0a1q2D\nzp3NA9lPP7U7IhFpSkryYSAy0tTmb7vNdLPUIuEi4UNJPkw4HDB7NixdCjfdBM89Z3dEItIUVJMP\nQ/v3w/jxMGQIPPUU/OAHdkckIrVRTV7qrUcPeP99OHoUkpJM7xsRCU0ayYexigp49FF49lmzUPjQ\noebnLpeL1zMz2b1xI9bRozRzu6lwOnG0a0dccjLjJ08mMjLS3uBFwoSmUIrPsrJg6lT47W9P0+q/\nv2f/yhVMLCxkkNvN2VPr3cAOp5NXYmPplZLClLQ0oqKi7ApbJCwoyYtfFBSc5OorSkk4WcArTOd8\nTng8frfTyeIhQ0hfu5aWLVs2UZQi4Uc1efHZiRMn+N/bktlxMp6ufMllfEAB/T3+Tn+3mwc2bWLu\nyJGcOOH5D4KI2EdJPsyVl5dz35gxLNi6lY78l8XM5H94mKt4m/9jBp7GD22B+Vu3ct+YMZSXlzdV\nyCJSD0ryYS5z0SJm5uVx0Vk/m8JSNjGcPzKHX7CME7So9ffbAql5eWRmZDR6rCJSf0ryYczlclG0\nYgX93O5q+2LZyzZ+TDT/YTDb2UFCrefp73azb/lyXC5XY4YrIg2gJB/GVi9dysTCwlr3R/MNL3A7\nD/IA15JFBndRgaPGY28sLGT10qWNFaqINJCSfBgryMlhUA2j+HNNYiVbSWQVE7mO9XxBTLVjBrnd\nFOTkNEaYIuIDJfkwZh09irct5rtRzGaG8WO2MZBdvM74KvudgFVW5vcYRcQ3SvJhrJkXo/izRXKa\nhTzAq/yMNDK4lRc4SfPvz3f6tL9DFBEfKcmHsYoGLhV1Be+zi4FYOIjnI95liDlfhNaFFwk0+n9l\nGHO0a4cbvC7ZnK0lJ/kzt7KGsdzMSibxEuddUODvEEXERxrJh7G45GR2+Ljw6zjW8g8G8A9HF158\n949s3eqn4ETEL5Tkw9j4yZNZFRvr83ku4ksG9V/I7xe1ZPx4+O1v4Ztv/BCgiPhMST6MRUZG0jMl\nhd0+juZ3O530Sknh5z+PoKAASkogPh7efddPgYpIg6kLZZgrLy8nbeRI5m3eXKW1gbfKgAeHDydj\nw4YqbYdfew1mzYJx40zP+lat/BaySFhRF0rxSVRUFOlvvMG8xESO1fN3y4D5iYmkr11bra/89dfD\nnj1mYZJ+/UzSF5Gmp5G8AKbd8H1jxjAzL6/GXjbnqk8/+U2b4PbboU8fePpp6NzZX1GLhD6N5MUv\nWrZsSUZ2NvkPPcTdcXF84HRybqp3Ax84ndwdF0f+woVkbNjg1YIhw4fDP/4BgwfDoEGQkQHqZSbS\nNDSSl2pcLherly6lICcHq6yMZqdPUxERgaNtW5/XeC0qgl//Gg4fNmvLDhvm5+BFQoyW/5OgY1nw\n97/DXXeZUf4f/gDt29sdlUhgUrlGgo7DARMnQmEhdOgA/fvDokWgxaVE/E9JXmzTooWZXrllC7zz\nDsTFwbp1dkclElpUrpGAsW4d/OY30K2beTjbr5/dEYnYT+UaCRmjR0NBAVx7LSQlwa9+BWpRL+Ib\nJXkJKFFRMGcOfPwxREZC376Qng7/+Y/dkYkEJyV5CUgXXghPPgn5+bBzJ/TuDUuWQD3XOREJe3Um\n+aysLPr06UPPnj157LHHqu3Pzc2lVatWJCQkkJCQwEMPPdQogUp46tEDXnnFbEuWwIAB8PrrZhqm\niNTN44NXt9tN7969ycnJoUOHDlx22WUsX76cvn37njkmNzeXxx9/nDVr1nj+ID14FR9ZFqxfD/fe\nC9HR8PDDcNVVdkcl0rga9cHrtm3b6NGjB127diUyMpJJkyaxevXqascpeUtTcDjMw9mdO+HOO2Hm\nTJPk8/LsjkwkcHlM8qWlpXTq1OnM9x07dqS0tLTKMQ6Hgy1bthAfH8/o0aMpLCxsnEhFvtOsGfz8\n5+Zlql/8Am65xczI0apUItV5XOPV4XDUeYJBgwZRUlJCdHQ069evZ8KECezbt6/GY+fPn3/m66Sk\nJJKSkuoVrMjZIiJg+nST6JcsgRtvNG/PzpsHiYl2RyfSMLm5ueTm5vrtfB5r8vn5+cyfP5+srCwA\nHnnkEZo1a8Y999xT6wm7devG9u3bueCCC6p+kGry0sj++1+T7NPTTVvj+++HoUPtjkrEN41ak7/0\n0kspKiqiuLiY8vJyVq5cybhx46occ+TIkTMBbNu2DcuyqiV4kabwgx+YOv3+/XDzzTBtmmmA9tZb\nmo0j4ctjuSYiIoJnnnmGUaNG4Xa7mTFjBn379mXx4sUApKamsmrVKp577jkiIiKIjo5mxYoVTRK4\nSG2iomDGDJg6FVasgLQ08wdg7lz42c/AxyVtRYKKetdIyKuogDffhEceMW0S7roLfvlLOO88uyMT\nqZv6yYt4ybLgvfdMW+P8fNMb51e/grZt7Y5MpHZqUCbiJYfDPIhdvdqsO3vwIPTqBbfdZhYdFwlF\nSvISlvr0gRdegH37zMLiyclwzTWmrFNRYXd0Iv6jco0IZvrlyy+bpmjHj8OsWaZu37q13ZFJuFO5\nRsQPfvADmDwZPvgAMjPN27PduplSzq5ddkcn0nBK8iJncTjgiitg+XLYuxe6doWxY+Hyy+HFF+Gb\nb+yOUKR+VK4RqcPp02ZpwsWLzQj/llvg1lvNmrQijU1TKEWaUHEx/OUvZuvY0ST7m26C88+3O7Lv\nuVwuXs/MZPfGjVhHj9LM7abC6cTRrh1xycmMnzyZyMhIu8MULynJi9jA7TbtEl54Ad55B8aNM83S\nhg0zXTLtUF5eTuaiRRStWMHEwkIGud2c/XKvG9jhdPJKbCy9UlKYkpZGVFSUPcGK15TkRWx29Cj8\n7W9mdH/ypHmAO2WKWdWqqZw4cYL7xoxhZl4e/bxYI3G308niIUNIX7uWli1bNkGE0lBK8iIBwrLM\ngiaZmebBbffupg3yTTfBRRc13ueeOHGCuSNHsmDrVurzMWXA/MREHs3OVqIPYEryIgHI5TLlnL/9\nzTy0HTYMUlJMWadFC/99Tnl5OWkjRzJv8+Z6JfhKZcCDw4aRkZ2t0k2A0jx5kQAUGQljxpgR/cGD\nZkGTZcvMw9pJk8xi5N9+6/vnZC5axMy8vAYleIC2QGpeHpkZGb4HIwFJSV6kkbVsaWr069aZXvdJ\nSfDUU3DxxaZ+v2ZNwxK+y+WiaMUKr2rwnvR3u9m3fDkul8un80hgUpIXaUIXXWQWNnn7bfOy1U9+\nAhkZ0L69Wbf273+HU6e8O9fqpUuZ6Kc1lW8sLGT10qV+OZcEFiV5EZu0bw+//rXpiPnJJ2YVq//9\nX7jkEpgwwbxh+9VXtf9+QU4Og3wcxVca5HZTkJPjl3NJYFGSFwkAMTGQmgrZ2eaFqxtuMC2Ru3Uz\n5Z0nnoADB6r+jnX0KP5a5MoJWGVlfjqbBBKPy/+JSNNr08bU6idPNr1yNm40Cf+xx+CCC8wD3TFj\nAJd/Z6s1O33ar+eTwKAkLxLAzjvv+6ReUQHbt8Mbb8CcObC34DU+YT2jWce1ZBHDUZ8+qyJC6SAU\nqVwjEiSaNYPLLoMFC2DHDrjjp3MZyVusZSy9+YRBbOc+HmYTwyinfr1p3IBD6yCGJCV5kSB1xbhL\niXdmsoobKaMtTzGbZlSQRgZtKWMsa3iaWXxMb+oq7OxwOolLTm6SuKVp6Y1XkSDlcrn43eDBPFZQ\nUG3fMS5kI1ezgWvIIZnTRJBMDlezkRG8QycOVjn+7rg4Ht6+Xd0pA5CvuVNFOJEgFRkZSc+UFHYX\nFtL/nKmUF/ElN/MyN/MyFnCA7uSQzFrGkkYGrTnOCN5hBO8Q0+w9eqWkKMGHKI3kRYJYQ3rXVOBg\nN/15hxG8RRJvR1xFxy7nk5TkYNgw02enSxezSpbYTw3KRMKcr10o09/K5rPPWrJpE7z7rtkiImDI\nELjySvNvXBw4/TUpX+pFSV5E/NpP3rJMj528PHjvPbMdOgSJiWb928svN1+3adNYVyNnU5IXEcCU\nbpY+/jifvPQSN/p5ZahjxyA/H7ZsMdv27aajZmLi91tcnOm+Kf6lJC8iVbhcLlYvXUpBTg5WWRnN\nTp+mIiICR9u2flvj9fRp2LMH3n8ftm0zC5wXF0N8vJnLX7n17GnfcoihQkleRALCv/9tRvgffggf\nfGC2L7+EhAQYPNhsCQnQu7fq+/WhJC8iAevLL83budu3m3937IAvvjClnYEDTdIfOBD694foaLuj\nDUxK8iISVL7+Gnbtgo8+Mmvi7toFH38MnTubcs+AAWaLizNTOcOt3ONyuXg9M5PdGzdiHT3Kwo0b\nuf/qq3G0a9egcpuSvIjYzuUyif6jj6Cg4Pvt+HHo188k/P79zdexsWZVrVCbx19eXk7mokUUrVjB\nxLMenDsAi4Y/OFeSF5GA9a9/mQe8u3ebpF9YaL53uUyy79u36hasI39PU2Ark/zZPE2BPZeSvIgE\nnbIyk/D37jVbYaFZXevYMTOjp08f6NXLPOTt3dv8rHVru6OuWV0vs9WU5OH7l9kezc72mOiV5EUk\nZJw8aco++/aZpH/2v82bm2RfufXoYbbu3aFVK3vi9aYtRW1JHkyif3DYMDKys2st3SjJi0jIsywz\nq+eTT8zbvEVFZjtwwHwfHW2Sfffu8KMfmX+7dTNfX3JJ4035/L/0dC5/4AGPbyl7SvJgSjf5Cxdy\n67331vz7SvIiEs4sC44cMcn+n/8024ED8Omn5uuvvoJOnUzS79rVbF26fP/vxRc37I+Ap1bRZ6sr\nyYPnVtFK8iIiHnzzDXz2mXmjt7jYJP+zv//Xv6BDBzMFtEsX8wehc2fzb+VWUzlo1V/+Qpfbb+ey\nOnoNeZPkP3A6+ez555k4fXr131eSFxFpuG+/hZIS+Pxzk/xLSqpvDofp9VO5degAO7L+zG0fvk5n\nSulAKW0po1kN6dybJO8GHkxJYcFLL1Xbp0VDRER88MMffv8wtyaWZVo+lJRAaSkcPGi2w4fa8wIz\nKaUDh7iE47QmhiNczGEu4dCZf+GhOmNwAlZZWY37lORFRBqRw2HKNa1amRe6Klm5i5h/KPfM9+VE\n8gXtKaUDh7n4zOatZqdP1/hzJXkRERtUnPO0NgoXnSmhMyVVfv6wt+eLqDmdB+G7YSIiwc/Rrh11\nL+/iHTfgaNu2xn1K8iIiNohLTmaHnybg73A6iUtOrnGfkryIiA3GT57MqthYv5zrldhYxk+eXOM+\nJXkRERtERkbSMyWF3T6O5nc7nfRKSam1/XCdST4rK4s+ffrQs2dPHnvssRqPmT17Nj179iQ+Pp6d\nO3f6FLCISLiYkpbG4iuv5FgDf78MWDxkCFPS0mo9xmOSd7vdzJo1i6ysLAoLC1m+fDl79+6tcsy6\ndevYv38/RUVFPP/889xxxx0NDDe45ebm2h1CownlawNdX7AL5uuLiooi/Y03mJeYWO9EX9mFMn3t\nWo995T0m+W3bttGjRw+6du1KZGQkkyZNYvXq1VWOWbNmDVOnTgUgMTGR48ePc+TIkXqGG/yC+X9o\ndQnlawNdX7AL9utr2bIlj2Zns2DYMPZ4WbrZ7XTy4PDhdbYZhjqSfGlpKZ06dTrzfceOHSktLa3z\nmIMHD3oVqIiImESfkZ1N/kMPcXdcHB84ndWmV7oxPWrujosjf+FCMjZsqDPBQx0vQzm8XF/r3J40\n3v6eiIgYUVFRzJg7F1daGquXLuWNnBzTqiAnhweSk3G0bUtccjIP13ONVywP3n//fWvUqFFnvk9P\nT7ceffTKFO23AAAEMklEQVTRKsekpqZay5cvP/N97969rS+++KLaubp3725h+uxo06ZNmzYvt+7d\nu3tK03XyOJK/9NJLKSoqori4mEsuuYSVK1eyfPnyKseMGzeOZ555hkmTJpGfn0/r1q2JiYmpdq79\n+/d7+igREWkEHpN8REQEzzzzDKNGjcLtdjNjxgz69u3L4sWLAUhNTWX06NGsW7eOHj160Lx5c5Ys\nWdIkgYuISN2arJ+8iIg0Pb+/8RrKL0/VdW25ubm0atWKhIQEEhISeOihuvtAB4rp06cTExNDXFxc\nrccE632Duq8vmO8dQElJCSNGjKBfv37079+fp556qsbjgvUeenN9wXoPv/32WxITExk4cCCxsbHc\nW8tarQ2+dz5V9M9x+vRpq3v37tann35qlZeXW/Hx8VZhYWGVY958803ruuuusyzLsvLz863ExER/\nhtBovLm2d955xxo7dqxNEfpm8+bN1o4dO6z+/fvXuD9Y71uluq4vmO+dZVnW4cOHrZ07d1qWZVkn\nTpywevXqFTL/37Ms764vmO/hqVOnLMuyLJfLZSUmJlrvvvtulf2+3Du/juRD+eUpb64NCNolDocO\nHUqbNm1q3R+s961SXdcHwXvvANq3b8/AgQMBaNGiBX379uXQoUNVjgnme+jN9UHw3sPo6GgAysvL\ncbvdXHDBBVX2+3Lv/JrkQ/nlKW+uzeFwsGXLFuLj4xk9ejSFhYVNHWajCdb75q1QunfFxcXs3LmT\nxMTEKj8PlXtY2/UF8z2sqKhg4MCBxMTEMGLECGLP6U7py73z68pQofzylDcxDho0iJKSEqKjo1m/\nfj0TJkxg3759TRBd0wjG++atULl3J0+eZOLEiTz55JO0aNGi2v5gv4eeri+Y72GzZs3YtWsXX3/9\nNaNGjSI3N5ekpKQqxzT03vl1JN+hQwdKSr5fuqqkpISOHTt6PObgwYN06NDBn2E0Cm+urWXLlmf+\ns+u6667D5XLx1VdfNWmcjSVY75u3QuHeuVwubrjhBn7xi18wYcKEavuD/R7WdX2hcA9btWrFT3/6\nUz788MMqP/fl3vk1yZ/98lR5eTkrV65k3LhxVY4ZN24cmZmZAB5fngo03lzbkSNHzvy13bZtG5Zl\nVautBatgvW/eCvZ7Z1kWM2bMIDY2ljlz5tR4TDDfQ2+uL1jv4bFjxzh+/DgA33zzDdnZ2SQkJFQ5\nxpd759dyTSi/POXNta1atYrnnnuOiIgIoqOjWbFihc1Rey8lJYVNmzZx7NgxOnXqxIIFC3C5XEBw\n37dKdV1fMN87gLy8PJYtW8aAAQPOJIj09HQ+//xzIPjvoTfXF6z38PDhw0ydOpWKigoqKiqYPHky\nV199td/ypl6GEhEJYVr+T0QkhCnJi4iEMCV5EZEQpiQvIhLClORFREKYkryISAhTkhcRCWFK8iIi\nIez/A9NMN23LWZ54AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x4932750>"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<!--\n",
      "s = la.lstsq(J, -r)[0]\n",
      "-->\n",
      "(Edit this cell for solution.)"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}