docs/examples/ex_tikhonov.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": "# Example: Matrix Equations using Tikhonov Regularization\nThis is an example of the use of matrix expressions in ``symfit`` models. This is illustrated by performing an inverse Laplace transform using Tikhonov regularization, but this could be adapted to other problems involving matrix quantities.\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%% md\n"
        }
      }
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "pycharm": {
          "is_executing": false
        }
      },
      "outputs": [],
      "source": "from symfit import (\n\tvariables, parameters, Model, Fit, exp, laplace_transform, symbols, \n\tMatrixSymbol, sqrt, Inverse, CallableModel\n) \nimport numpy as np\nimport matplotlib.pyplot as plt\n"
    },
    {
      "cell_type": "markdown",
      "source": "Say $f(t) \u003d t * exp(- t)$, and $F(s)$ is the Laplace transform of $f(t)$. Let us first evaluate this transform using ``sympy``.  \n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%% md\n"
        }
      }
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "pycharm": {
          "is_executing": false,
          "metadata": false,
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "text": [
            "[F(s; ) \u003d (s + 1)**(-2)]\n"
          ],
          "output_type": "stream"
        }
      ],
      "source": "t, f, s, F \u003d variables(\u0027t, f, s, F\u0027)\nmodel \u003d Model({f: t * exp(- t)})\nlaplace_model \u003d Model(\n\t{F: laplace_transform(model[f], t, s, noconds\u003dTrue)}\n)\nprint(laplace_model)\n"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "pycharm": {
          "is_executing": false,
          "metadata": false,
          "name": "#%% md\n"
        }
      },
      "source": "Suppose we are confronted with a dataset $F(s)$, but we need to know $f(t)$. This means an inverse Laplace transform has to be performed. However, numerically this operation is ill-defined. In order to solve this, Tikhonov regularization can be performed.\n\nTo demonstrate this, we first generate mock data corresponding to $F(s)$ and will then try to find (our secretly known) $f(t)$."
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "pycharm": {
          "is_executing": false,
          "metadata": false,
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "data": {
            "text/plain": "\u003cmatplotlib.legend.Legend at 0x2b51a9f8e80\u003e"
          },
          "metadata": {},
          "output_type": "execute_result",
          "execution_count": 6
        },
        {
          "data": {
            "text/plain": "\u003cFigure size 432x288 with 1 Axes\u003e",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEMCAYAAAAxoErWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFTxJREFUeJzt3X2QXuV53/HvDyFFgEkgaB1cFpAMCkWmVLJ3HFMGg4PMAG0RdVODYzsyxcgzDk7aqu2QpuN4iBmHxCUtU+xITSk4xrwYM60mocWF4IIxeLTYioRENGhkbDbQImQDqpGsl179Y1ea1eo50kraPc++fD8zGvac597zXM8i+O11Xu47VYUkSZ0c0+0CJEkTlyEhSWpkSEiSGhkSkqRGhoQkqZEhIUlqZEhIkhoZEpKkRoaEJKmRISFJanRstws4WnPmzKm5c+d2uwxJmjSeffbZ16qqZzRjJ31IzJ07l/7+/m6XIUmTRpIfjnasp5skSY0MCUlSI0NCktRo0l+TkDQ17Nq1i4GBAXbs2NHtUqaM2bNn09vby8yZM4/4GIaEpAlhYGCAE088kblz55Kk2+VMelXF1q1bGRgYYN68eUd8nNZONyW5M8mrSZ5reD1Jbk+yKcnaJO9uqzZJ3bdjxw5OOeUUA2KMJOGUU0456s6szWsSdwGXH+T1K4D5Q3+WAV9uoSZJE4gBMbbG4ufZWkhU1RPAjw8yZAnwlRr0DHBSkncc6ribt/x0rEqUJI0wke5uOg14adj2wNA+SeromhVPc82Kp7tdxpQ2kUKiU19UHQcmy5L0J+nftWvXOJclabp78cUXOe6441i4cOF++5cvX86CBQu44YYbWLhwIbNmzeK11147ovdYsWIF73jHO1i4cCELFy7k4x//OADbt2/n4osvZs+ePR2/b+fOnbz//e9n9+7dR/S+hzKR7m4aAE4ftt0LvNxpYFWtBFYC/OKZ53YMEkk6Gvfffz+PPfYYZ599Nh/+8Ic566yzWLNmzb7XN2/ezFNPPcWGDRv27TuaeeTWrl3L5z//ea6//vr99t9555186EMfYsaMGR2/b9asWVx66aXcf//9fPSjHz3i928ykTqJVcBvDN3l9D7gjap65VDf9M6eE8a/MknTyu23387GjRu57rrreOihhw54fePGjVx88cX88Ic/ZNGiRfz0p0d/bXTdunUHdCoA99xzD0uWLNm3fffdd/Oe97yH888/n4suugiAq6++mnvuueeoa+iktU4iyb3AJcCcJAPA7wEzAarqT4CHgSuBTcBbwHVt1SZJe7311lvcfPPNvPDCC+zevZuvfOUrB4w555xzWLp0KXPnzuWTn/zkQY930UUXsW3btgP2f/GLX2Tx4sX7ttevX891113HMcccw5w5c3j00UfZuXMnmzdv3tehbNu2jVtvvZU1a9Ywa9YsXn/9dQDOO+88Vq9efRSfullrIVFVHznE6wX8ZkvlSFJHjz/+OGeeeSYnn3wyAD09Pbz44osHjFu3bt1+v+E3efLJJw855qWXXuLUU09l7dq1++1/7bXXOOmkk/Ztz5gxg+3bt7N8+XKWLl1KX1/fvv2zZs1i27ZtnHjiiYd8v8MxkU43SVLXPfHEE1x22WWHHLd+/Xre9a53HXLcRRddtO9i9PA/jz766L4xa9eu7Xis4447br+H4Y4//niee+45LrzwQpYtW8aXvvSlfa/97Gc/Y/bs2Yes53BNpAvXktR13/nOd1i+fDlbtmyhp6eHN998k7feemu/Mdu2bWPmzJkcf/zxhzzeaDqJdevWdQyJk08+mT179rBjxw5mz57NCy+8wPz587n22mvZsGHDvgDZunUrPT09RzVHUxM7CUka8q1vfYuZM2dyyy23cNZZZ3HFFVfw1a9+9YDf0J977jnOO++8MXvfdevWsWDBgo6vXXbZZXz7298G4JZbbuGcc87h3e9+Nz/4wQ/49Kc/DQyeIrvyyivHrJ7h7CQkacgll1zCJZdccsD+kdckLrjgAr7+9a+P2fse7M6kG2+8kdtuu43Fixdz1113dRzzta99jS984QtjVs9whoSkSev+T13QyvvMmDGDN954g4ULF+73rMRe27dv54ILLmDXrl0cc8zYnqBZtGgRH/jAB9izZ0/HZyV27tzJ1VdfzTnnnDOm77tXBm8qmrz6+vrKNa6lye/555/n3HPP7XYZU06nn2uSZ6uqbzTf7zUJSVIjQ0KS1MiQkDRhTPbT3xPNWPw8DQlJE8Ls2bPZunWrQTFG9i5ferQP2Hl3k6QJobe3l4GBAbZs2dLtUqaM2bNn09vbe1THMCQkTQgzZ85k3rx53S5DI3i6SZLUyJCQJDUyJCRJjQwJSVIjQ0KS1MiQkCQ1MiQkSY0MCUlSI0NCktTIkJAkNTIkJEmNDAlJUiNDQpLUyJCQJDUyJCRJjaZMSFyz4mmuWfF0t8uQpCllyoSEJGnsGRKSpEathkSSy5NsTLIpyU0dXj8jyeNJvp9kbZIr26xPkrS/1ta4TjIDuAP4IDAArE6yqqo2DBv2b4EHqurLSRYADwNzR3P8+z91wRhXLElqs5N4L7CpqjZX1U7gPmDJiDEF/PzQ178AvNxifZKkEVrrJIDTgJeGbQ8AvzJizOeAbyb5DHACsLid0iRJnbTZSaTDvhqx/RHgrqrqBa4E/izJATUmWZakP0n/li1bxqFUSRK0GxIDwOnDtns58HTS9cADAFX1NDAbmDPyQFW1sqr6qqqvp6dnnMqVJLUZEquB+UnmJZkFXAusGjHmR8ClAEnOZTAkbBUkqUtaC4mq2g3cCDwCPM/gXUzrk9yc5KqhYcuBG5L8FXAv8ImqGnlKSpLUkjYvXFNVDzN4W+vwfZ8d9vUG4MI2a5IkNfOJa0lSI0NCktTIkJAkNTIkJEmNDAlJUiNDQpLUyJCQJDUyJCRJjQwJSVIjQ0KS1MiQkCQ1MiQkSY0MCUlSI0NCktTIkJAkNTIkJEmNDAlJUiNDQpLUaEqGxDUrnuaaFU93uwxJmvSmZEhIksbGsd0uYDzc/6kLul2CJE0JdhKSpEaGhCSpkSEhSWpkSEiSGhkSkqRGhoQkqZEhIUlqZEhIkhoZEpKkRoaEJKlRayGR5PIkG5NsSnJTw5gPJ9mQZH2Sr7VVmySps1bmbkoyA7gD+CAwAKxOsqqqNgwbMx/4HeDCqvpJkre3UZskqVlbncR7gU1VtbmqdgL3AUtGjLkBuKOqfgJQVa+OxRs7bbgkHbm2QuI04KVh2wND+4b7ZeCXkzyV5Jkkl7dUmySpQVtThafDvhqxfSwwH7gE6AWeTHJeVb1+wMGSZcAygDPOOGNsK5Uk7dNWSAwApw/b7gVe7jDmmaraBfwgyUYGQ2P1yINV1UpgJUBfX9/IsNmPa0tI0pFr63TTamB+knlJZgHXAqtGjPmvwAcAksxh8PTT5pbqkyR10EpIVNVu4EbgEeB54IGqWp/k5iRXDQ17BNiaZAPwOPCvqmprG/VJkjpL1UHP1kx4fX191d/f3+0yJGnSSPJsVfWNZqxPXEuSGhkSkqRGhoQkqZEhIUlqdNghkeSEobmYJElT3CFDIskxSX49yV8keRX4a+CVoZla/2hoYj5J0hQ0mk7iceAsBmdoPbWqTq+qtwMXAc8Af5DkY+NYoySpS0YzLcfioaky9lNVPwa+AXwjycwxr0yS1HWH7CQ6BcSRjJEkTT6HPcFfkpuBGcAaYE1VvTDmVUmSJoTDDomq+mySXwIWAf84yVlVdcPYlyZJ6rZRh0SS/wX8w6p6E/hHwGzgtqGV5iRJU9DhPCdxUlW9meQ9DC41ejLwn8anLEnSRHA4p5t2JTkW+A3g1qp6IMmkmn5171rXLkQkSaNzOCFxO/BXDJ5mumlo39vGvCJJ0oQx6pCoqq8keQjYU1Xbk5wNPD1+pY09OwhJOjyHDIkkqaGViarq/+7dX1WbgOtGjpEkTR2jmpYjyWeSnDF8Z5JZSX41yd3A0vEpT5LUTaM53XQ58E+Be5O8E/gJg9clZgDfBP64qtaMX4mSpG4ZTUisqKqlwJeG5miaA2yvqtfHtzRJUreN5nTT+cO+/ouqesWAkKTpYTQhMfyCdM94FSJJmnhGc7rp1CSfYPAZiYxvOZKkiWQ0IfE5oI/B2117k6wD1g/92VBV3xi/8iRJ3XTIkKiqlcO3k/QyeJ3i7wBXM7jwkCRpCjqSqcIHgAHg4bEvR5I0kRzOLLBTzjUrnt436Z8k6UCH3UlMJc7lJEkHN607CUnSwRkSkqRGhoQkqVGrIZHk8iQbk2xKctNBxv1akkrS12Z9kqT9tRYSSWYAdwBXAAuAjyRZ0GHcicBvAd9tqzZJUmdtdhLvBTZV1eaq2gncByzpMO73gT8EdrRYmySpgzZD4jTgpWHbA0P79kmyCDi9qv78YAdKsixJf5L+LVu2jH2lkiSg3ZDoNDngvhlmkxwD/DGw/FAHqqqVVdVXVX09PWMzMa0P1knSgdoMiQHg9GHbvcDLw7ZPBM4DvpXkReB9wCovXktS97T5xPVqYH6SecDfANcCv773xap6g8FV7wBI8i3gX1ZVfxvF+fS1JB2otU6iqnYDNwKPAM8DD1TV+iQ3J7mqrTokSaPX6txNVfUwI2aPrarPNoy9pI2aJEnNfOJaktTIkJAkNTIkJEmNDAlJUiNDQpLUyJDowKevJWmQISFJajSt17hu4tPXkjTITkKS1MiQkCQ1MiQkSY0MCUlSI0NCktTIkDgEn5mQNJ15C+wheDuspOnMTkKS1MiQkCQ1MiQkSY0MCUlSI0PiMHink6TpxrubDoN3OkmabuwkJEmNDAlJUiND4ih4jULSVOc1iaPgNQpJU52dhCSpkSEhSWpkSIwRr09Imoq8JjFGvD4haSqyk5AkNWotJJJcnmRjkk1Jburw+r9IsiHJ2iSPJTmzrdokSZ21EhJJZgB3AFcAC4CPJFkwYtj3gb6qOh94EPjDNmqTJDVrq5N4L7CpqjZX1U7gPmDJ8AFV9XhVvTW0+QzQ21JtY86L2JKmirYuXJ8GvDRsewD4lYOMvx747+Na0TjyIrakqaKtkEiHfdVxYPIxoA+4uPFgyTJgGcAZZ5wxFvVJkjpo63TTAHD6sO1e4OWRg5IsBn4XuKqqftZ0sKpaWVV9VdXX09Mz5sVKkga1FRKrgflJ5iWZBVwLrBo+IMkiYAWDAfFqS3WNO69PSJrMWjndVFW7k9wIPALMAO6sqvVJbgb6q2oV8EfA24CvJwH4UVVd1UZ948nrE5Ims9aeuK6qh4GHR+z77LCvF7dViyRpdHziWpLUyJCQJDUyJCRJjQyJLvLOJ0kTnVOFd5F3Pkma6OwkJgi7CkkTkZ3EBGFXIWkispOYgOwqJE0UdhITkF2FpInCTmKCs6uQ1E12EhOcXYWkbrKTkCQ1MiQkSY0MCUlSI0NCktTIkJiCvCNK0ljx7qYpyDuiJI0VO4lpxA5D0uGyk5hGhncYe8PCrkPSwRgS09TIcDA0JHViSAiwy5DUmSGhAxgOkvbywrWOiBfBpenBTkJHxG5Dmh7sJNQauw9p8rGTUGvsPqTJx05Ck4adiNQ+OwlNGqPtRIbfwuvtvNLRSVV1u4aj0tfXV/39/d0uQ5ImjSTPVlXfaMZ6ukmS1MiQkCQ1ajUkklyeZGOSTUlu6vD6zyW5f+j17yaZ22Z9kpp548D01NqF6yQzgDuADwIDwOokq6pqw7Bh1wM/qaqzk1wL3Apc01aNkppN9ov/3sRwZNq8u+m9wKaq2gyQ5D5gCTA8JJYAnxv6+kHgPyZJTfar65K6znA4Mm2ebjoNeGnY9sDQvo5jqmo38AZwSivVSZIO0GZIpMO+kR3CaMaQZFmS/iT9W7ZsGZPiJEkHajMkBoDTh233Ai83jUlyLPALwI9HHqiqVlZVX1X19fT0jFO5kqQ2Q2I1MD/JvCSzgGuBVSPGrAKWDn39a8Bfej1CkrqntQvXVbU7yY3AI8AM4M6qWp/kZqC/qlYB/xn4sySbGOwgrm2rPknSgVqdu6mqHgYeHrHvs8O+3gH8kzZrkiQ184lrSVIjQ0KS1MiQkCQ1mvRThSfZBmzsdh1dMgd4rdtFdJGf38/v5z8yZ1bVqJ4fmAqLDm0c7bzoU02S/un62cHP7+f387fx+T3dJElqZEhIkhpNhZBY2e0Cumg6f3bw8/v5p7dWPv+kv3AtSRo/U6GTkCSNk0kbEodaCnUqS3J6kseTPJ9kfZLf7nZN3ZBkRpLvJ/nzbtfStiQnJXkwyV8P/T2YNivqJPnnQ3/vn0tyb5LZ3a5pPCW5M8mrSZ4btu8Xk/zPJC8M/fPk8Xr/SRkSw5ZCvQJYAHwkyYLuVtWq3cDyqjoXeB/wm9Ps8+/128Dz3S6iS/4D8D+q6m8Df5dp8nNIchrwW0BfVZ3H4GShU30i0LuAy0fsuwl4rKrmA48NbY+LSRkSDFsKtap2AnuXQp0WquqVqvre0NfbGPwfxMhV/qa0JL3A3wf+tNu1tC3JzwPvZ3DWZKpqZ1W93t2qWnUscNzQmjPHc+C6NFNKVT3BgevqLAHuHvr6buDq8Xr/yRoSo1kKdVpIMhdYBHy3u5W07t8D/xr4f90upAveCWwB/svQ6bY/TXJCt4tqQ1X9DfBF4EfAK8AbVfXN7lbVFb9UVa/A4C+NwNvH640ma0iMapnTqS7J24BvAP+sqt7sdj1tSfIPgFer6tlu19IlxwLvBr5cVYuAnzKOpxsmkqFz70uAecDfAk5I8rHuVjW1TdaQGM1SqFNakpkMBsQ9VfVQt+tp2YXAVUleZPBU468m+Wp3S2rVADBQVXu7xwcZDI3pYDHwg6raUlW7gIeAv9flmrrh/yR5B8DQP18drzearCExmqVQp6wkYfB89PNVdVu362lbVf1OVfVW1VwG/93/ZVVNm98mq+p/Ay8lOWdo16XAhi6W1KYfAe9LcvzQfweXMk0u2o8wfKnnpcB/G683mpQT/DUthdrlstp0IfBxYF2SNUP7/s3Qyn+aHj4D3DP0S9Jm4Lou19OKqvpukgeB7zF4l9/3meJPXie5F7gEmJNkAPg94A+AB5Jcz2BwjtuKnj5xLUlqNFlPN0mSWmBISJIaGRKSpEaGhCSpkSEhSWpkSEiSGhkSkqRGhoQ0xpIsTfJskrVJnux2PdLR8GE6aQwlOZHBGXkXVtXOJCdNs2m8NcXYSUhjaw9wHPDvkvTtDYgkv9/dsqQjY0hIY6iq3gLOA54CVib5dJJTmaTzpEn+xZXGUJL5VfUCcN/QkrKzGVwUas3Bv1OamOwkpLH1u0k2JvkegwvjfAlYiCGhScpOQhpDVfWJkfuSnA280H410tHz7iZJUiNPN0mSGhkSkqRGhoQkqZEhIUlqZEhIkhoZEpKkRoaEJKmRISFJamRISJIa/X/OI+Fb3I2mmgAAAABJRU5ErkJggg\u003d\u003d\n"
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": "epsilon \u003d 0.01  # 1 percent noise\ns_data \u003d np.linspace(0, 10, 101)\nF_data \u003d laplace_model(s\u003ds_data).F\nF_sigma \u003d epsilon * F_data\nnp.random.seed(2)\nF_data \u003d np.random.normal(F_data, F_sigma)\n\nplt.errorbar(s_data, F_data, yerr\u003dF_sigma, fmt\u003d\u0027none\u0027, label\u003dr\u0027$\\mathcal{L}[f] \u003d F(s)$\u0027)\nplt.xlabel(r\u0027$s_i$\u0027)\nplt.ylabel(r\u0027$F(s_i)$\u0027)\nplt.xlim(0, None)\nplt.legend()\n"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "pycharm": {
          "is_executing": false,
          "metadata": false,
          "name": "#%% md\n"
        }
      },
      "source": "We will now invert this data, using the procedure outlined in \\cite{}.\n"
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "pycharm": {
          "is_executing": false,
          "metadata": false,
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "text": [
            "[W(I, M; a) \u003d (I + a**(-2)*M)**(-1),\n c(Fs, W; ) \u003d -W*Fs,\n d(c; ) \u003d (c.T*c)**(1/2)]\n"
          ],
          "output_type": "stream"
        }
      ],
      "source": "N_s \u003d symbols(\u0027N_s\u0027, integer\u003dTrue)  # Number of s_i points\n\nM \u003d MatrixSymbol(\u0027M\u0027, N_s, N_s)\nW \u003d MatrixSymbol(\u0027W\u0027, N_s, N_s)\nFs \u003d MatrixSymbol(\u0027Fs\u0027, N_s, 1)\nc \u003d MatrixSymbol(\u0027c\u0027, N_s, 1)\nd \u003d MatrixSymbol(\u0027d\u0027, 1, 1)\nI \u003d MatrixSymbol(\u0027I\u0027, N_s, N_s)\na, \u003d parameters(\u0027a\u0027)\n\nmodel_dict \u003d {\n    W: Inverse(I + M / a**2),\n    c: - W * Fs,\n    d: sqrt(c.T * c),\n}\ntikhonov_model \u003d CallableModel(model_dict)\nprint(tikhonov_model)"
    },
    {
      "cell_type": "markdown",
      "source": "A ``CallableModel`` is needed because derivatives of matrix expressions sometimes cause problems.\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%% md\n"
        }
      }
    },
    {
      "cell_type": "markdown",
      "source": "Build required matrices, ignore s\u003d0 because it causes a singularity.\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%% md\n"
        }
      }
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "outputs": [
        {
          "name": "stdout",
          "text": [
            "d [[0.01965714]]\n"
          ],
          "output_type": "stream"
        }
      ],
      "source": "I_mat \u003d np.eye(len(s_data[1:]))\ns_i, s_j \u003d np.meshgrid(s_data[1:], s_data[1:])\nM_mat \u003d 1 / (s_i + s_j)\ndelta \u003d np.atleast_2d(np.linalg.norm(F_sigma))\nprint(\u0027d\u0027, delta)\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%%\n",
          "is_executing": false
        }
      }
    },
    {
      "cell_type": "markdown",
      "source": "Perform the fit\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%% md\n"
        }
      }
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "outputs": [],
      "source": "model_data \u003d {\n\tI.name: I_mat,\n\tM.name: M_mat,\n\tFs.name: F_data[1:],\n} \nall_data \u003d dict(**model_data, **{d.name: delta})\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%%\n",
          "is_executing": false
        }
      }
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "pycharm": {
          "is_executing": false,
          "metadata": false,
          "name": "#%%\n"
        }
      },
      "outputs": [
        {
          "name": "stdout",
          "text": [
            "\nParameter Value        Standard Deviation\na         5.449374e-02 None\nStatus message         Optimization terminated successfully.\nNumber of iterations   14\nObjective              \u003csymfit.core.objectives.LeastSquares object at 0x000002B51BEDA390\u003e\nMinimizer              \u003csymfit.core.minimizers.BFGS object at 0x000002B51AB3A518\u003e\n\nGoodness of fit qualifiers:\nchi_squared            3.3130494560395173e-19\nobjective_value        1.6565247280197586e-19\nr_squared              -inf\n"
          ],
          "output_type": "stream"
        }
      ],
      "source": "fit \u003d Fit(tikhonov_model, **all_data)\nfit_result \u003d fit.execute()\nprint(fit_result)\n"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "pycharm": {
          "is_executing": false,
          "metadata": false,
          "name": "#%% md\n"
        }
      },
      "source": "Check the quality of the reconstruction\n"
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "outputs": [
        {
          "name": "stdout",
          "text": [
            "(100,) (100,)\n"
          ],
          "output_type": "stream"
        },
        {
          "data": {
            "text/plain": "\u003cmatplotlib.legend.Legend at 0x2b51bf2a6d8\u003e"
          },
          "metadata": {},
          "output_type": "execute_result",
          "execution_count": 23
        },
        {
          "data": {
            "text/plain": "\u003cFigure size 432x288 with 1 Axes\u003e",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAENCAYAAAD0eSVZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHv1JREFUeJzt3Xt4VPW97/H3N3cI4ZogmHBTMYqIgDnqLtZWpRZsudTTVnzUWnuhbkuvPj3bnna3e7vrtrbdvWwPWmm31VrUqrWWdlN121qVVixBIjelUkCIgITIHUJu3/PHmoQhJJkhmcyay+f1PPOstX7rNzPfafGTNb9Za/3M3RERkcySE3YBIiKSeAp3EZEMpHAXEclACncRkQykcBcRyUAKdxGRDKRwFxHJQAp3EZEMpHAXEclAeWG9cWlpqY8dOzastxcRSUsrV67c7e5lsfqFFu5jx46luro6rLcXEUlLZvZmPP00LCMikoEU7iIiGUjhLiKSgUIbcxcRiaWpqYna2loaGhrCLiXpioqKqKioID8/v0fPV7iLSMqqra2lpKSEsWPHYmZhl5M07k59fT21tbWMGzeuR6+hYRkRSVkNDQ0MGzYsq4IdwMwYNmxYr76xxAx3M7vPzHaZ2dou9puZ/aeZbTSz1WY2tcfViIh0kG3B3qa3nzueI/f7gRnd7J8JjI885gP39KoiERHptZjh7u4vAO9002UO8HMPLAcGm9nIWK+7qe5Q/FWKiMhJScSYezmwLWq7NtJ2AjObb2bVZlbd1NSUgLcWEel79957LyNHjmTy5MlMnjyZ66+/PuZzjhw5wnve8x5aWlo63d/Y2Mgll1xCc3NzossFEhPunQ0MeWcd3X2Ru1e5e1VPT+8REUm21atX861vfYuamhpqamp48MEHYz7nvvvu46qrriI3N7fT/QUFBVx++eX88pe/THS5QGLCvRYYFbVdAWxPwOuKiKSENWvWMHny5JN6zuLFi5kzZ0779gMPPMD555/PpEmTePe73w3A3LlzWbx4cUJrbZOI89yXAAvM7BHgQmCfu++I9aTTyooT8NYiIn1v3bp13HjjjeTk5FBaWsqzzz7bbf/GxkY2bdpE251vDxw4wJ133klNTQ0FBQXs3bsXgIkTJ7JixYo+qTlmuJvZw8B7gVIzqwW+CeQDuPuPgaXAlcBG4DBwY59UKiLZ7fe3ws41iX3NEefCzG9322Xbtm2MGDGC1atXd7rf3U84bXH37t0MHjy4fTs3N5cjR45wyy23cMMNN1BVVdXeXlBQwIEDBygpKenlhzlezHB392ti7HfgswmrSEQkhaxevZpzzjnnuLadO3dy1VVXMXv2bK699lruueceDh8+TGNjI3fffTf9+vU77gKk/v37s3btWn77298yf/58PvWpT3HzzTcDcPToUYqKihJet24/ICLpIcYRdl9Zs2bNCeG+atUq5s2bx+c//3kWLlzIkSNHGDx4MJs2bQJgyJAhtLS00NDQQFFREW+88Qbjx49n3rx5rF+/vj346+vrKSsr6/H9Y7qjcBcR6caaNWuYO3fucW01NTXtbatWrWLhwoUUFhYe1+eKK65g2bJlTJ8+ndtvv52XXnqJ4uJizjnnHH7yk58A8Nxzz3HllVf2Sd0KdxGRbnR2Nssbb7xBZWUlAHPmzOHjH/84o0aN4rLLLmPGjOCC/gULFvD973+f6dOnc//993f62g899BB33HFHn9StcBcROUn33Xdf+/qsWbOYNWvWCX2mTJnCpZdeSktLS6fnujc2NjJ37tz2PxKJpnAXEekjn/jEJ7rcV1BQwMc+9rE+e2/d8ldEJAMp3EVEMpDCXUQkAyncRUQykMJdRCQDKdxFRDKQwl1EJAMp3EUko1x970tcfe9LYZcROoW7iEgMnU2zF/Y0erEo3EVEYuhsmr2wp9GLReEuIhJDZ9PshT2NXiy6t4yISAwdp9lbunRp6NPoxaJwFxHpRmfT7G3fvj30afRi0bCMiEg3Optmr6tp9KZNm8b8+fO5++672/f11TR6sSjcRUS60dk0e9HT6EEweUdxcTHz5s3jgx/8YFKm0YtF4S4i0o01a9YwYcKEE9rbptEDuP3226msrGTq1Kls3ry5ffLrvpxGLxaNuYuIdKOrs13CnkYvFoW7iGSUX37mH5LyPmFPoxeLwl1EpIfCnEYvFo25i4hkIIW7iEgGUriLSEpz97BLCEVvP3fo4a7bc4pIV4qKiqivr8+6gHd36uvre3XxU+g/qF545EXe1fA8+O/ALOxyRCSFVFRUUFtbS11dXdilJF1RUREVFRU9fn7o4V7WsouLGpZBw17oNyTsckQkheTn5zNu3Liwy0hLcQ3LmNkMM9tgZhvN7NZO9o82s+fMbJWZrTazuC/J2p1bFqzsq433KSIiEkPMcDezXGAhMBOYAFxjZh2vxf068Ki7TwHmAXcTpy99+PJgReEuIpIw8Ry5XwBsdPdN7t4IPALM6dDHgYGR9UHA9rgrGBQZU1K4i4gkTDxj7uXAtqjtWuDCDn3+BXjGzD4HFAPT466geDjk5CvcRUQSKJ4j985OYel4XtI1wP3uXgFcCTxoZie8tpnNN7NqM6tu//U7JwcGlSvcRUQSKJ5wrwVGRW1XcOKwyyeBRwHc/SWgCCjt+ELuvsjdq9y9qqys7NiOQaMU7iIiCRRPuK8AxpvZODMrIPjBdEmHPluBywHM7GyCcI//xNRBFQp3EZEEihnu7t4MLACeBl4jOCtmnZndZmazI91uAT5tZq8CDwMf95O5pGxQBRzYDi3NJ/0BRETkRHFdxOTuS4GlHdq+EbW+HpjW4yoGVYC3woEdMHhU7P4iItKt0O8tA+h0SBGRBEuRcI8crSvcRUQSIjXCfWB5sNy3rft+IiISl9QI98IBwU3DdOQuIpIQqRHuEIy7738r7CpERDJCCoW7LmQSEUmUFAr3Co25i4gkSGqFe8M+aNgfdiUiImkvtcIdNO4uIpIAqRPuA3Uhk4hIoqROuLdfpapxdxGR3kqdcC8ZAZarI3cRkQRInXDPyQ2uVFW4i4j0WuqEO+i+7iIiCZKC4a4xdxGR3kq9cN+/A1pbwq5ERCStpV64tzbBwV1hVyIiktZSK9wHjw6We7aEWoaISLpLqXBf8MzBYGX338ItREQkzaVUuO/OHc5RCqFuQ9iliIiktZQK90dumkbhiErYrXAXEemNlAp3AMoqoU7DMiIivZF64V5aCfu2QuOhsCsREUlbqRfuZWcGS/2oKiLSYykY7mcFSw3NiIj0WOqF+9DTICdPP6qKiPRC6oV7bn4Q8DodUkSkx1Iv3AFKz1S4i4j0QmqGe1klvLMJmhvDrkREJC2laLifBd4SBLyIiJy01Az30rbTITU0IyLSE3GFu5nNMLMNZrbRzG7tos9HzWy9ma0zs4d6VVXp+GCpcXcRkR7Ji9XBzHKBhcD7gFpghZktcff1UX3GA18Fprn7HjMb3quqCoph0GiFu4hID8Vz5H4BsNHdN7l7I/AIMKdDn08DC919D4C79362jbJKNr++iqvvfanXLyUikm3iCfdyIHpi09pIW7QzgTPN7M9mttzMZvS6srJKypu3Ya4p90RETlbMYRnAOmnzTl5nPPBeoAJ40cwmuvve417IbD4wH2D06NHdv2vpmRTQSFmLptwTETlZ8Ry51wKjorYrgO2d9PmNuze5+2ZgA0HYH8fdF7l7lbtXlZWVdf+uZZUA3PW+4jhKFBGRaPGE+wpgvJmNM7MCYB6wpEOfJ4FLAcyslGCYpncnqbfdQOzttb16GRGRbBQz3N29GVgAPA28Bjzq7uvM7DYzmx3p9jRQb2brgeeAr7h7fa8q6zcYhoyD7TW9ehkRkWwUz5g77r4UWNqh7RtR6w58OfJInFOnQG11Ql9SRCQbpOYVqm1OnRzMynSod18CRESyTYqH+5RguWNVuHWIiKSZ1A73kecFS427i4iclNQO96JBMPR02K4jdxGRk5Ha4Q7BuPuOV8OuQkQkraRBuE+Bfdvg0O6wKxERSRupH+4jJwdLjbuLiMQtDcK97UdVjbuLiMQr9cO9aCAMOwN26MhdRCReqR/uEIy768hdRCRu6RHuIyfD/rfgoG7/KyISj/QI97YrVfWjqohIXNIj3EdOAkxDMyIicUqPcC8sYVveaF79y1NhVyIikhbSI9yB9QXnUtm0Hlqawi5FRCTlpU24v/8DH6bIG3QrAhGROKRNuDNmWrDc8mK4dYiIpIH0CfcBw6H0TNjy57ArERFJeekT7gBjL4aty6GlOexKRERSWnqF+5hp0HgAdq4OuxIRkZSWXuE+9uJguWVZuHWIiKS49Ar3khHBzExvatxdRKQ76RXuEBy9v/kStLaEXYmISMpKz3A/ug92rgm7EhGRlJV+4d52vruGZkREupR+4T6onJ25I1nxpyVhVyIikrLSL9yBNYVTOLexBpoawi5FRCQlpWW4v2/ux4P7zOiUSBGRTqVluDPuEsjvD3/7fdiViIikpPQM9/x+cNqlsOEpcA+7GhGRlJOe4Q5QOQP218Lba8OuREQk5cQV7mY2w8w2mNlGM7u1m34fNjM3s6rEldiF8e8Plhs0O5OISEcxw93McoGFwExgAnCNmU3opF8J8Hng5UQX2amSU6D8fI27i4h0Ip4j9wuAje6+yd0bgUeAOZ30+zfgO0Dyzk+snAlvrYQDbyftLUVE0kE84V4ObIvaro20tTOzKcAod/9ddy9kZvPNrNrMquvq6k662BOcOTNYvvF0719LRCSDxBPu1klb+ykqZpYD/AC4JdYLufsid69y96qysrL4q+zKKedQlzucFU8v7v1riYhkkHjCvRYYFbVdAWyP2i4BJgJ/MrMtwEXAkqT8qGrGysKLmHR0FRw90OdvJyKSLuIJ9xXAeDMbZ2YFwDyg/cYu7r7P3Uvdfay7jwWWA7PdvbpPKu5gxrwFFHIUXvttMt5ORCQtxAx3d28GFgBPA68Bj7r7OjO7zcxm93WBMY26AIaMg1cfCbsSEZGUkRdPJ3dfCizt0PaNLvq+t/dlnQQzmHQ1PH8n7HsLBpXHfo6ISIZL3ytUo036KOCw5rGwKxERSQmZEe7DToeKC2D1L3WvGRERMiXcAc67Gnat1/R7IiJkUrifcxXk5AdH7yIiWS5zwr3/UDjz/cG4e0tz2NWIiIQqc8Id+O7bU+Hg2/A33SlSRLJbRoX7K4UXUpc7HF7+cdiliIiEKqPC/eGbLqbssgWw5UXYqUk8RCR7ZVS4AzD1Y8H8qi/fE3YlIiKhybxw7zcEzpsHqx+DQ7vDrkZEJBSZF+4AF94ELUdh5c/CrkREJBSZGe5llXD6ZfDXn0JzY9jViIgkXWaGO8BFN8PBnSy869thVyIiknSZG+5nTGdz3un874MPQUtT2NWIiCRV5oa7GeM+8u+MaNkBNQ+FXY2ISFJlbrhDcDuC8ip4/jvQfDTsakREkiazw90MLvsa7K+FV34edjUiIkmT2eEOcNqlMPpd8ML3oOlI2NWIiCRF5oe7GVz2dTi4E/66KOxqRESSIvPDHWDsNFYWXsCRZ++A/TvCrkZEpM9lR7gD53/mXvrltMIzXw+7FBGRPpc14c7Q0+DiL8Hax2HzC2FXIyLSp7In3AEu/iIMHgNLv6ILm0Qko2VXuOf3g5l3Qt3r/OKH/xR2NSIifSa7wh2gciZUfoDrjiyGXa+HXY2ISJ/IvnAHmPVDKCiGJz6lu0aKSEbKznAfMBxm3wU718Cf/j3sakREEi47wx3grA/AlOtpXfZDvnnXT8KuRkQkobI33AFm3EHOkDH8a/OP4FB92NWIiCRMdod7YQl8+GdwcBc8dgO0NIddkYhIQsQV7mY2w8w2mNlGM7u1k/1fNrP1ZrbazP5gZmMSX2ofKZ8Ks34EW16E//nnsKsREUmImOFuZrnAQmAmMAG4xswmdOi2Cqhy90nA48B3El1on5p8DVz4j7D8bqh5OOxqRER6LZ4j9wuAje6+yd0bgUeAOdEd3P05dz8c2VwOVCS2zCS44t9YW3AejU9+DrYsC7saEZFeiSfcy4FtUdu1kbaufBL4fW+KCkVuPhO/+CQFpafBw9fAjlfDrkhEpMfiCXfrpM077Wh2HVAFfLeL/fPNrNrMquvq6uKvMln6D4Xrfw1Fg+DBq2D3xrArEhHpkXjCvRYYFbVdAWzv2MnMpgNfA2a7e6cTlrr7InevcveqsrKyntTb9waVBwEP8OCHYO+27vuLiKSgeMJ9BTDezMaZWQEwD1gS3cHMpgD3EgT7rsSXmWSl4+G6X3Fofz11d10G9X8PuyIRkZMSM9zdvRlYADwNvAY86u7rzOw2M5sd6fZdYADwmJnVmNmSLl4ufZw6meL5v6essAV+NhPeXh92RSIicTP3TofP+1xVVZVXV1eH8t4nZdfr8PM50HIUrv0VVJwfdkUiksXMbKW7V8Xql91XqMZj+FnwiaegcCDcfyWsfSLsikREYlK4x2PoOPj0H2HkZHj8RvjTnRDSNx4RkXgo3ONVXAo3LIHzrgluE/zYDdCwL+yqREQ6pXA/GXmFMPceeN9t8Nrv4N5LYPuqsKsSETmBwv1kmcG0L8CNS4NJtv/rClh+D7S2hl2ZiEg7hXtPjb4IblpGdd4UeOpWeGAWvLMp7KpERACFe+/0H0rVrU/DnIWwczXcMw2W/xhaW8KuTESynMK9t8xgynVw83IYMw2e+idY9B7Y+nLYlYlIFlO4J8qgcrj2MfjI/XD4HbjvCvj1TbDvrbArE5EspHBPJDM450Pw2b/y6+KraXr1cbhrKjzzz0Hgi4gkicK9LxQO4ENfWUT+F1bChLnwl7vgPyfD89+BI3vDrk5EsoDCvS8NGQNX3Qs3LQvG45+7HX4wEZ79VziYgvezF5GMoXBPhhET4ZqHg5AfP53WZT+g6Xtnw5OfhZ1rwq5ORDJQXtgFZJUR58JH7ifn0jfIefnHUPMQ1PwCxlwMVTfCWR+E/KKwqxSRDKBb/obp8Dvwys9h5c9gzxboNxTOmxc8RkwKfqAVEYkS7y1/Fe6poLUVNj8fhPzrS6G1CcrOhkkfhYlXwZCxYVcoIilC4Z6uDr8D637Nhmd+SmVTZPankZNhwhw4e1YwBaCIZC2FeybY8yas/w2sfxLeWhm0DT0dKmfC+CuC+9vkFYZbo4gklcI90+zdBn97Cjb8nua/P08ezZBfDGMvhtMvg3GXwPCzNU4vkuHiDXedLZMuBo+CCz4NF3yavKMHYMsy2PgsbPwDvPF00Ke4LDiffsy7gqP6UyZCTm64dYtIKBTu6aiwJBiaqZwZbO95E7a8CJtfgC1/DoZxAApKoHwqVPyv4FE+FQYMD69uEUkahXsmGDImeEy5LtjeuxW2LoetL0FtNSz7AXjkNsQDy+HUKTDyvOC8+xHnBm0azhHJKAr3TDR4dPCY9NFgu/EQbK+BHTXBtIDbV8Hr/w0Ev7cctAEMGDUJTpkQjNuXnQWllcG8sQp9kbSkcM8GBcUwdlrwaHP0ALy9HnauZsDb62DXeg5XP0R/P3ysT78hMGx8cPrlsDNg2Okw9DQYMg4KByT/c4hI3BTu2aqwBEZfGDwi+rvD/u1Q9zrs/ltkuTH44bZm8fHPLx4eXFw1ZGwwJDRo1LFvDAPLdRsFkZAp3OUYs2DSkUHlcMblx+9r2A97NkP934O5YvdsCR5bl8Pax8E7TBBeXBaE/KAKGHhq8CgZGfU4BQoHathHpI8o3CU+RQODH2FHnnfivpam4Ih/75uwrzaYfWrftmD9nU2w+UU4uu/E5+X1C0J+wCnBH4MBw4NvBMWlwXZxKfQvDZb9hui0TpGToHCX3svNP3bGTleOHoQDO+HAjshjJxx8+9ij/u/w5l/gSBczVllOEPD9hkL/YdB/aGS77TE4WBYNDtaLItuFAyFX/8wl++hfvSRH4QAoPANKz+i+X0tTcH+dQ3XB43A9HNoNh3cH7Yfrg8febbDj1aCt+Uj3r1kwIAj5ooHHLwtLjq0XDAi2CwcE1wcUDgjaCoqD9oJiyCvSMJKkDYW7pJbc/GCopuSU+J/T1AANe+HInmAaw4a90LAvWD+6P1hv2Bv8btCwL/jjsGdLcMZQw/7YfxzaWE4Q+Pn9g7Av6B/cAqKgf9CW3//49fx+HZZFwTKvKGiLXuYVBfvzijT8JAmhcJf0l18E+SOgZETPnt/SDI2RoG88BI0Hg+BvPBhsHz0Y7G88fGx/0+HI9sHgeQd2BvuajkQeh078kTleOfmRwC/sYlkAuYWR7cLIegHkRh7HtRUGfzDzCiP78yNtbesFUev5x9Zz2tbzIuv5kJOnby5pJK5wN7MZwI+AXOCn7v7tDvsLgZ8D5wP1wNXuviWxpYr0kdy8Y2P3ieIOLY1RYX84WDY3RNYbgm8MTQ1BW3NDZP/RY+0tRyPbDZFlZL3xMLTsgebGYLulMdjXvjyauM/RUU5eEPo5+ScGf9vyuPX84JtI9L7j+uSe2H5CW2TbcqO2czts5wXfrE7onxNZdtF23LKz9pzI63bo296Wun/sYoa7meUCC4H3AbXACjNb4u7ro7p9Etjj7meY2TzgTuDqvihYJC2YHTuy7jc4ue/tDq3NxwK/PfSbIttHg28rXa23tvWLtLc2RbVH9rU2B+sd97U2R+1rhtaWyB+kg5H2SN/2fs3BrTHa97U9mnr+zSfZjvtDkHss9I/7I5BzbP8Zl8OsH/Z5WfEcuV8AbHT3TQBm9ggwB4gO9znAv0TWHwf+n5mZh3U/YZFsZnZsmCWduQcB39J07A9Aa9Qyus1bu9jfybLTfa0ntntrMEtaZ33bHsf1benQ3tYvur01uL1HEsQT7uXAtqjtWuDCrvq4e7OZ7QOGAbsTUaSIZKG2o1/9wNwjOXH06WxQqeMReTx9MLP5ZlZtZtV1dXXx1CciIj0QT7jXAqOitiuA7V31MbM8YBBwwtUo7r7I3avcvaqsrKxnFYuISEzxhPsKYLyZjTOzAmAesKRDnyXADZH1DwN/1Hi7iEh4Yo65R8bQFwBPE5wKeZ+7rzOz24Bqd18C/BfwoJltJDhin9eXRYuISPfiOs/d3ZcCSzu0fSNqvQH4SGJLExGRnopnWEZERNKMwl1EJAMp3EVEMpCFdVKLmR0ANoTy5qmhlOy+yCubP382f3bQ5+/t5x/j7jHPJQ/zrpAb3L0qxPcPlZlV6/Nn5+fP5s8O+vzJ+vwalhERyUAKdxGRDBRmuC8K8b1TgT5/9srmzw76/En5/KH9oCoiIn1HwzIiIhkolHA3sxlmtsHMNprZrWHUEAYzG2Vmz5nZa2a2zsy+EHZNYTCzXDNbZWa/C7uWZDOzwWb2uJm9Hvl38A9h15RMZvalyL/9tWb2sJkVhV1TXzKz+8xsl5mtjWobamb/Y2ZvRJYJnN/xmKSHe9S0fTOBCcA1ZjYh2XWEpBm4xd3PBi4CPptFnz3aF4DXwi4iJD8CnnL3s4DzyKL/HcysHPg8UOXuEwluRJjpNxm8H5jRoe1W4A/uPh74Q2Q74cI4cm+fts/dG4G2afsynrvvcPdXIusHCP7DLg+3quQyswrgA8BPw64l2cxsIHAJwV1UcfdGd98bblVJlwf0i8z70J8T54bIKO7+AifObTEHeCCy/gAwty/eO4xw72zavqwKOAAzGwtMAV4Ot5Kk+yHwf4A0mf04oU4D6oCfRYalfmpmxWEXlSzu/hbwPWArsAPY5+7PhFtVKE5x9x0QHPABw/viTcII97im5MtkZjYA+BXwRXffH3Y9yWJmHwR2ufvKsGsJSR4wFbjH3acAh+ijr+SpKDK2PAcYB5wKFJvZdeFWlbnCCPd4pu3LWGaWTxDsi939ibDrSbJpwGwz20IwHHeZmf0i3JKSqhaodfe2b2uPE4R9tpgObHb3OndvAp4A3hVyTWF428xGAkSWu/riTcII93im7ctIZmYE462vufv3w64n2dz9q+5e4e5jCf5//6O7Z82Rm7vvBLaZWWWk6XJgfYglJdtW4CIz6x/5b+FysugH5SjR05LeAPymL94k6TcO62ravmTXEZJpwPXAGjOribT938hMV5IdPgcsjhzYbAJuDLmepHH3l83sceAVgjPHVpHhV6ua2cPAe4FSM6sFvgl8G3jUzD5J8AevT2ax0xWqIiIZSFeoiohkIIW7iEgGUriLiGQghbuISAZSuIuIZCCFu4hIBlK4i4hkIIW7ZCUz+4yZ7TCzmsjjQTPrZ2bPR25L3dlzCszshcgdDUVSmsJdstUk4OvuPjnyuB74BPCEu7d09oTILar/AFydxDpFekThLtnqXKCmQ9u1RN3nw8xuMLOVZrbazF6MND8Z6SeS0nT7AclKZlYPvEVwX/ndwJXAVncfEdlfQnCv/cnu3mhmg919b2TIZqe7l4VVu0g8NHYoWcfMRhEE9KSotlOB6FmRWoB+wH+Y2QPuXg3g7i1m1mhmJZHZtERSkoZlJBtNAjreifQI0D5Zs7sfBiYCfwYWmdnNUX0LgYa+LlKkNxTuko3OpUO4u/seINfMigDMbLy7H3L3R4DfEQl+MxsGtE02IZKyFO6Sjc6l80kyngEujqx/zcw2mNkrBNPC3R1pvxTQ/fcl5ekHVZEIM5sCfDlyWmRXfZ4AvuruG5JXmcjJ05G7SIS7rwKe6+4iJuBJBbukAx25i4hkIB25i4hkIIW7iEgGUriLiGQghbuISAZSuIuIZCCFu4hIBlK4i4hkoP8PD3f6obVq7DcAAAAASUVORK5CYII\u003d\n"
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": "ans \u003d tikhonov_model(**model_data, **fit_result.params)\nF_re \u003d - M_mat.dot(ans.c) / fit_result.value(a)**2\nprint(ans.c.shape, F_re.shape)\n\nplt.errorbar(s_data, F_data, yerr\u003dF_sigma, label\u003dr\u0027$F(s)$\u0027, fmt\u003d\u0027none\u0027)\nplt.plot(s_data[1:], F_re, label\u003dr\u0027$F_{re}(s)$\u0027)\nplt.xlabel(r\u0027$x$\u0027)\nplt.xlabel(r\u0027$F(s)$\u0027)\nplt.xlim(0, None)\nplt.legend()",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%%\n",
          "is_executing": false
        }
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "pycharm": {
          "is_executing": false,
          "metadata": false,
          "name": "#%% md\n"
        }
      },
      "source": "Reconstruct $f(t)$ and compare with the known original.\n"
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "outputs": [
        {
          "data": {
            "text/plain": "\u003cmatplotlib.legend.Legend at 0x2b51c00fda0\u003e"
          },
          "metadata": {},
          "output_type": "execute_result",
          "execution_count": 25
        },
        {
          "data": {
            "text/plain": "\u003cFigure size 432x288 with 1 Axes\u003e",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAENCAYAAADjW7WQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VNX5+PHPM1v2PWFLgIRdFgGNoIAIiihWgVqtuNSlLlVB22prbf2qrbb9WWvVqtSt4lKllOJGW9QqbgUECassRiBsAQJZyZ5JMuf3xx0whECGkOQmM8/79ZrXzL333JlnXJ65Ofec54gxBqWUUqHBYXcASiml2o8mfaWUCiGa9JVSKoRo0ldKqRCiSV8ppUKIJn2llAohmvSVUiqEaNJXSqkQoklfKaVCiMvuABpLTk426enpLTo3OzsbgIEDB7ZiREop1fGtWrWqwBiT0ly7Dpf009PTycrKatG5EyZMAODTTz9tvYCUUqoTEJGdgbTT7h2llAohmvSVUiqEaNJXSqkQ0uH69JVS6kTV1taSm5tLdXW13aG0ufDwcNLS0nC73S06X5O+UqrTy83NJSYmhvT0dETE7nDajDGGwsJCcnNzycjIaNF7aPeOUqrTq66uJikpKagTPoCIkJSUdFJ/0WjSV0oFhWBP+Iec7PfU7p2W2v455GdDWiZ0HQZO/UeplOr4NFOdKG8lfHg/rPzrt/vcUTB4Klz8BLgj7ItNKaWaoUn/RORtgAU3QME3cOZMGHUT7FltXfWvfg0O5sKVf4ewGLsjVUqpJmmffqDqamDelVBdCj94By78PST2gWGXwdSn4LvPw85l8Np0qCyyO1qllE3uvvtuBg8ezM0338w555xDfX09YI0w+sc//gGA1+tl/Pjx1NXVtXt8mvQDteoVKNkF330W+k4EoLq2no827efu+eu4eV1fPhnxJ8y+9fC371o/EkqpkJKTk8PSpUvZtGkTI0aM4NJLL8XpdAKwePFiVq9eDYDH4+G88847/CPQnrR7JxA1ZfDZo5AxHvpMxBjDMx9v5fnPcyivqSM23EVcpJsbNnVhsmMmL+x7goPv/5a4ix+2O3KlQs5v/rWRTXtLW/U9B/eI5cFLhhy3TXZ2NpMmTaKuro6RI0cC8PbbbwOwZMkS7rrrLuLj4/nggw94++23mT59Or/85S+5+uqrWzXW5mjSD8TyZ6GyAM77NYjw7Kdb+dOH3zB5cFeuPrM3Z/VJwuNysL2ggsWbT+GtxWuZlvU0qxPGcdrYC+yOXinVDgYOHMh1111Heno61157Lb169eJQmfhx48Zxxhln8NhjjzF06FAA6uvrWblyZbvHqUm/ORWFsPQpGHQxpJ3OGyt28uj72Uwb0YMnvj8Ch+PbMbMZyVHcdHYfdvV7joIXxhH/wY95sXoBN5831MYvoFRoae6KvC199dVXTJs2jYKCAuLj4484lp2dfcRaH06nE4/HQ1lZGTEx7Tf4Q/v0m7PkcaitgHPv5z/r9/F/72zg3EFdeOzy4Uck/IZ6de9G/JUv0sexD9cnDzF3xa52DlopZYeNGzcyZMgQIiIijpg1W1hYSFxc3FH1cmpqaggPD2/XGDXpH4+vHtb9HQZPoyKuH/e/u4HhafH85erTcDuP/48ubMBEfJk3cZ3rv7y+cBHLtha0U9BKKTuUlZXhdruJjIwkISGB+vr6w4l/+/bt9OjR44j2hYWFpKSktLhwWksFlPRF5EIRyRaRrSJybxPHbxWRr0RkrYgsEZHB/v3pIlLl379WRJ5r7S/QpnJXQmUhnHIJr36xg6IKLw9eMphwtzOg0x3n3oeExfJgxD+59fVV5OSXt228SinbbNiw4XB/PcDkyZNZsmQJAIMGDaKgoIChQ4eybNkyAD755BMuuuiido+z2aQvIk5gNjAFGAxceSipNzDXGDPMGDMCeBR4vMGxbcaYEf7Hra0VeLvIXgQOF2VpE3jh8xwmDkxhZK+EwM+PTETOvovRdVmMdmziplezqPS2/7hcpVTbO+uss/jnP/95eHvWrFm8+uqrAERHR/Pll1+yYcMGxowZA8DcuXO55ZZb2j3OQK70RwFbjTE5xhgvMA+Y1rCBMabh+KgowLReiDbKfh96j+GVVUWUVNby0/MHnPh7jP4RxKbyeMJb5BSU8+j72a0fp1Kqwxk5ciQTJ048PDmrIa/Xy/Tp04+4sdteAkn6qcDuBtu5/n1HEJGZIrIN60r/zgaHMkRkjYh8JiJnn1S07akoBwqyqcqYzIv/y2HSKV05NS2++fMac0fAxPuIKVzPo4O28cqyHXy5XWfsKhUKfvjDHx6enNWQx+Ph2muvtSGiwJJ+U0NUjrqSN8bMNsb0BX4B/J9/9z6glzFmJHAXMFdEYo/6AJFbRCRLRLLy8/MDj74tZb8PwLyDQyitruMnk/q3/L2Gz4AuQ7is5GXSEzzcs2AdVd6jf/2VUqqtBZL0c4GeDbbTgL3HaT8PmA5gjKkxxhT6X68CtgFH9ZEYY14wxmQaYzJTUlICjb1tffMeJnkgz673cd6gLgxNjWv5ezmccO59OEq28+Jpu9hRWMlj/9VuHqVU+wsk6a8E+otIhoh4gBnAwoYNRKThZfB3gC3+/Sn+G8GISB+gP5DTGoG3qeqDsHMZed0mcqCshqkjejR/TnMGTIEug+n/zQtcO7onc5ZuZ+Pegyf/vkopdQKaTfrGmDpgFvABsBmYb4zZKCIPichUf7NZIrJRRNZideNc598/HlgvIuuABcCtxpiO36G99SPw1fF+7Qg8LgfnDupy8u/pcMDZd0P+1/wiYxsJkR4e+tcmjAmOe95Kqc4hoDIMxphFwKJG+x5o8PrHxzjvTeDNkwnQFtnvYyKT+Ov2JMb3TyQmvJUmTwz5LnzyO6JWPMFPJ73G/e9u5ION+7lwaLfWeX+llGqGzshtyo7/UdxtLHtKa7loWCsmZIcTxv0U9q3jqqQtDOgaze8XbaamTm/qKqXahyb9xsryoGwfq+r64HYK553StXXf/9QZEJuKc8nj3H/xYHYVVfLK0h2t+xlKKdscWkTljjvuOG67qqqqw4ustOcCK5r0G9u7FoC381IY1y+ZuIhWrovh8sCYO2DXMs4O38F5g7rw9MdbKarwtu7nKKXaXcNFVJ5++unjtp0zZ87hRVbac4EVTfqN7VuLQfi0tDtThnVvm88Y+QMIi4MvnuGXFw2i0lvH859va5vPUkq1i+zsbM455xx27tzJyJEjqaioOG77N954g2nTph1eYGXBggWMGDGC7du3M336dN544402iVOTfmN711IY0RuvI4LJg1u5a+eQsGjIvB42L6Sfu4jpI1J5ddkODpRVN3uqUqpjOrSIysMPP8yaNWuIioo6Zluv10tOTg7p6emHF1h59913Wbt2LRkZGQwdOrTNFljRRVQa27uGtXWDOLNPEvGRnrb7nFE/gi9mw4rnufO8+3h33V7+8sk2fj3VvgUglAoK790LeV+17nt2GwZTHmm22aFFVBoyxiByZGGDxoustOcCK3ql31BZHpTnsbSyJ2f1TWrbz4pLhSGXwurXSI+u4/LT05i7Yhd7Sqra9nOVUm3m0CIqeXl5jBkzhkceeYTc3Fx+9atf8ZOf/ITbb78d4IhFVtp7gRW90m/IfxP3K18GP+99AiWUW+qsmfDVfFj9GnecdxNvrd7DMx9v4f9demrbf7ZSwSqAK/K20HARlc8++4wZM2Zw5513Mnv2bKqqqoiPjycnxypI0HCRlfZeYEWv9BvauwaDsMWRwfCeLaioeaJ6jID0s2H5c6TGuLlyVE/mZ+Wyq7Cy7T9bKdWqGi6isnbtWs4//3wA1qxZwyOPPMKvf/1rXnvttcPtDy2y0t4LrGjSb2jfWnJdPemT2jXg1bFO2pm3Q2kufP0vbpvQD6cIz+lIHqU6nYaLqGzZsuVwH/20adO4/vrrueeee3j//fcPtz+0yEp7L7Ci3TsNmL1rWeXtT2Z7dO0cMuACSEiHFc/Tbch3uSwzjQVZudx5bn+6xbXvgslKqdYxZ86cw68vueQSLrnkkqPaNFxkpWHN/bZeYEWv9A8p3YeU57G2PoPTeye23+c6nNZInl1fwN613HZOX+qN4cX/dfxipEqpk9PUIittvcCKJv1D9n17E/f09rzSBxh5NXiiYcXz9EyMZNrwHryxYieF5TXtG4dSKuhp0j9k71p8OKhIOIWUmLD2/ezwOBhxFWxYAOUHuH1iX2rqfMxZur1941BKBT1N+n5m7xpySGVIeissmNISo26Bei9kvUy/LjFMGdqN15bt5GBVrT3xKKWCkiZ9v/p961lf35vM9Hbu2jkkuT/0mwRZL0Gdl9sn9KOspo65K3bZE49SnUyoLEh0st9Tkz6AtwJX+T62+Xq078idxkbfBuX7YdO7DE2N4+z+ycxZul3r7SvVjPDwcAoLC4M+8RtjKCwsPKmZugEN2RSRC4E/A07gr8aYRxodvxWYCdQD5cAtxphN/mO/BG70H7vTGPNBi6NtK0XWSJkDnjT6pkTbF0ffcyGxL3z5PJx6OT8a35drXlrBO2v2cMUZveyLS6kOLi0tjdzcXPLz8+0Opc2Fh4eTlpbW4vObTfr+hc1nA+cDucBKEVl4KKn7zTXGPOdvPxV4HLhQRAZjLaQ+BOgBfCQiA4wxHevStdCaDBXZbQAOhzTTuA05HDD6R/DePbBnFWP7ncaQHrE8/3kOl5/e097YlOrA3G43GRkZdofRKQTSvTMK2GqMyTHGeIF5wBFl5IwxpQ02o4BDf2NNA+YZY2qMMduBrf7361Bq87cAkNxrkM2RAMOv9A/ffAER4Ufn9CUnv4KPNu+3OzKlVBAIJOmnArsbbOf69x1BRGaKyDbgUeDOEzz3FhHJEpEsO/48q9j7NftMIhmpbVQ//0SEx8KIq2HDm1B+gIuGdiMtIYLnP9fJWkqpkxdI0m+qT+GouyXGmNnGmL7AL4D/O8FzXzDGZBpjMlNSUgIIqXXVF2xjh68bA7u2bt3qFht1C/hqYdUruJwObj67D6t2FpO1o8juyJRSnVwgST8X6NlgOw3Ye5z284DpLTzXFhGl29lJd9KTj73STbtK7mcN31xpDd+8PDONuAi3lmZQSp20QJL+SqC/iGSIiAfrxuzChg1EpH+Dze8AW/yvFwIzRCRMRDKA/sCXJx92K6oqJrKuhLKoXridHWgE6+hboTwPNi8k0uPimjN78d9N+9lecPx1N5VS6niazXLGmDpgFvABsBmYb4zZKCIP+UfqAMwSkY0isha4C7jOf+5GYD6wCXgfmNnxRu5YV8++pH42B9JI3/Os4ZsrngPgurPScTsczFmipRmUUi0X0Dh9Y8wiYFGjfQ80eP3j45z7O+B3LQ2wrVXvzyYciOrWNmVMW6zh8M3cVXRJO51pI3rwz1W7uev8ASREteH6vUqpoNWB+jPsUbJ7M/VG6JLewZI+WEXYPDGHr/ZvOrsP1bU+3lix0+bAlFKdVcgnfe+BLeSaFAZ0b+OF0FsiLAZGXgMb34ayPAZ2i2H8gBReWbaT6tqO1UumlOocQj7pu0ty2Cnd6ZkYaXcoTRt1M/jqIOtlAG4+O4OC8hoWru1wg6CUUp1AaCd9Y4iv2kVJRC+cHbXEQVJf6D/ZX32zhnH9khnULYY5S7cHfXEppVTrC+2kX36ACFNFXXwfuyM5vjNvhYp82PAWIsKN4zL4Oq+MpVsL7Y5MKdXJhHTSL9/7NQBhXfo309JmfSZCyimwfDYYw9QRPUiODuOvS3SyllLqxIR00s/fuRGAxN6DbY6kGSJw5m2Q9xXsWEKYy8m1Z/Xm0+x8tuwvszs6pVQnEtJJv3LfN9QYF70yOuBwzcZO/T5EJsHyZwG4enQvwlwOXUdXKXVCQjrpS9E2culKj4QOUnPneNwRkPlDyF4EhdtIig7j0tPSeGv1HgrLa+yOTinVSYR00o+p2EFBWE9EOujIncbOuAkcLljxPAA3jkunps7HG7qOrlIqQKGb9I0hpS4Pb0wnWoYwphsM/R6seR2qSujXJYYJA1N47YsdOllLKRWQkE365SUHCMeLJPRsvnFHctbtUFsBq18F4KZxfSgo97JwnU7WUko1L2STfv7urQCEJ/W2OZIT1H04ZIyH5c9BnZex/ZKsyVpLdLKWUqp5IZv0S/KsUS+xXTvhYspj7oSyvbBRJ2sppU5MyCb9qgKrUmWXtL42R9IC/SZByiBY9rRO1lJKnZCQTfqmJJca4yYuubvdoZw4ERhzB+zfADmfEOZycp1O1lJKBSCgpC8iF4pItohsFZF7mzh+l4hsEpH1IrJYRHo3OFYvImv9j4WNz7WLq3wv+c4UxNFJf/eGXQ7RXa2rfeDqM3sT5nLwkq6spZQ6jmYznog4gdnAFGAwcKWINK5bsAbINMacCiwAHm1wrMoYM8L/mEoHEV2TR5mnq91htJwrzFpZa9vHkPcViVEevnd6Gm+t2UN+mU7WUko1LZDL3FHAVmNMjjHGC8wDpjVsYIz5xBhT6d9cDqS1bpityxhDUt0BvFE97A7l5GT+EDzRsPTPANw4LgNvnY+/LdeVtZRSTQsk6acCuxts5/r3HcuNwHsNtsNFJEtElovI9BbE2OryD5aTQjHEd+jfpuZFJMDp18OGN6FoO31Topl0ShdeX64raymlmhZI0m+qRkGTA8JF5BogE/hjg929jDGZwFXAkyJy1HAZEbnF/8OQlZ+fH0BIJ2f/nu04xRCW2Ilm4x7LWbOs0gz+vv2bzu5DUYWXN1fn2hyYUqojCiTp5wINp62mAUdN/xSRScB9wFRjzOFOZWPMXv9zDvApMLLxucaYF4wxmcaYzJSUlBP6Ai1Rss+62RnTGcfoNxbbHYbPsEozlB9gdEYiw1LjeGnJdnw+nayllDpSIEl/JdBfRDJExAPMAI4YhSMiI4HnsRL+gQb7E0QkzP86GRgLbGqt4FuqIt8qUJac2sFXzArU2J9AvReW/wUR4aazM8jJr+Djrw80f65SKqQ0m/SNMXXALOADYDMw3xizUUQeEpFDo3H+CEQD/2w0NPMUIEtE1gGfAI8YY2xP+r5iK+mHJQVB9w5Y6+gOngYrX4Lqg1w0rDup8RG88D+drKWUOpIrkEbGmEXAokb7HmjwetIxzlsGDDuZANuCq3wvZRJDjKcT1NEP1LifwqZ34MsXcY//GTeMTee3/9nM2t0ljOgZb3d0SqkOopPOTDo5kdV5HOzMY/Sb0mME9J8MX8yGmnJmjOpFTLiLFz/Xq32l1LdCLul763wk1R2gurOP0W/K+Hugqgiy5hAd5uKaM3vz3oZ97CqsbP5cpVRICLmkv7ekih5SALHHm2rQSfU8A/pMhGVPgbeSG8ak43SIFmJTSh0Wckl/z/4DxEklnmC5idvYOfdART6sfpUuseFMH5HK/KzdFFV47Y5MKdUBhFzSLwqmMfpN6T0Geo+DJU9CbTW3jO9Dda2P177YYXdkSqkOIOSSfmX+DgDigjXpg3W1X54Hq1+jf9cYJp3ShVeX7aDSW2d3ZEopm4Vc0q8rtsoIOTp73Z3jyRgPvcfC/x4DbyW3TehLcWUt81fubv5cpVRQC7mk7yzbSz0OiO5mdyhtRwQm3gfl+yHrJU7vncgZ6Qm8+L/t1Nb77I5OKWWjkEv6UVX7KHOngDOgeWmdV/pYayTPkiegppxbz+nLnpIq/r3+qLJJSqkQElJJv8pbT3J9PpURQXyV39C5/weVhfDl80wc2IWBXWN49tNtWohNqRAWUkn/QFk13aWQupggHKPflLRM6H8BLH0Kh7eUH53Th2/2l/NJthZiUypUhVTSzyuporsU4QjGiVnHcu59UF0Cy57mkuE9SI2PYPYnWzFGr/aVCkUhlfSLig4QJrWEJYRQ0u8+HIZcCl/Mxl2Zz63n9GH1rhK+yCm0OzKllA1CKulXFOwBICo5hJI+WH379V747A9cntmTlJgwZn+y1e6olFI2CKmkX1NsJf2IxCAstnY8SX2ttXRXvUJ46Q5uPjuDpVsLWbOr2O7IlFLtLKSSvq80DwCJ6W5zJDYYfw+4wuHjh7l6dG/iI916ta9UCAoo6YvIhSKSLSJbReTeJo7fJSKbRGS9iCwWkd4Njl0nIlv8j+taM/gT5ajYb72IDrJa+oGI6QpjZsHGt4kqWMcNYzL4aPMBNu8rtTsypVQ7ajbpi4gTmA1MAQYDV4rI4EbN1gCZxphTgQXAo/5zE4EHgdHAKOBBEUlovfBPTFjVAaockRAWbVcI9hpzB0Qmwwf/x/Vn9SY6zMUzH+vVvlKhJJAr/VHAVmNMjjHGC8wDpjVsYIz5xBhzaKWO5cChwjYXAB8aY4qMMcXAh8CFrRP6iTHGEOktoNKdbMfHdwxhMdYQzl3LiNv5HteN6c2iDfv4Zn+Z3ZEppdpJIEk/FWhYqSvXv+9YbgTea+G5baa0uo4kSvBGpNjx8R3HyGuhy2D48AFuOjOVSLeTpxZvsTsqpVQ7CSTpSxP7mpzZIyLXAJnAH0/kXBG5RUSyRCQrPz8/gJBO3IHSarpQTH0o9uc35HTBBb+D4h0kbJjDtWPS+c9X+9iiV/tKhYRAkn4u0LPBdhpwVNUuEZkE3AdMNcbUnMi5xpgXjDGZxpjMlJS2uRLPK62mi5TgjAmRujvH0/dcqzzD549xy2kxRLidPK19+0qFhECS/kqgv4hkiIgHmAEsbNhAREYCz2Ml/IaFXT4AJotIgv8G7mT/vnZXWFRElNQQHmpj9I9l8m+htpKE5X/gB2f15l/r97L1QLndUSml2lizSd8YUwfMwkrWm4H5xpiNIvKQiEz1N/sjEA38U0TWishC/7lFwMNYPxwrgYf8+9pdeUEuAFFJQbx4yolIGQCjfgSrX+O2/qWEu7RvX6lQEFBReWPMImBRo30PNHg96TjnzgHmtDTA1lLrn43riQ/BiVnHMuFe2LCA+E9+yfVjnuG5z7czc2I/BnaLsTsypVQbCZkZufX+2bhBvWLWiQqPhfMfhj2ruCNhOdEeF098+I3dUSml2lDIJP3Ds3FjQnz0TmOnfh96nUXk57/l9rOSeH9jHl/lHrQ7KqVUGwmZpB9WnU+teCA83u5QOhYRuOiPUFXMjTV/Iz7SzeMfZtsdlVKqjYRE0vf5DNHeAio8SVaSU0fqNgxG34pnzas8MLycT7LzWbXTlvvtSqk2FhJJv7DCSzIleMO72B1KxzXxPohNZXruo3SNcvDo+9m6upZSQSgkkv5+/8QsX6jPxj2esGj4zmM48jfzbJ9lrNhexKfftM3saKWUfUIo6RfjjNWRO8c1cAqccgkjt7/AmfEH+cN7X+Pz6dW+UsEkJJJ+QfFB4qQytNbGbakpjyION0/HvMrXeaW8u26P3REppVpRSCT9ikJrNm5kkib9ZsX2gMkPkZK/nLsTl/Gn/35DTV293VEppVpJSCR9b4lV480Vq7NxA3L6DZAxntu8r+Ar3s3ry3fZHZFSqpWERNKvP+ifjasTswIjAlOfxiXwXNyrPL34Gw5W1todlVKqFYRE0ndU+At/agmGwCWkw/m/4dSaVVzg/ZCnP9ZibEoFg5BI+mHVB6jHCZFJdofSuWTeCOln85uwN1j8xZfsLKywOyKl1EkK+qRfV+8jpraQSk8SOIL+67YuhwOm/wWPy8kfXc/yh0Ub7Y5IKXWSgj4LFlV66SIl1ISH+Nq4LRXfC8d3HiNTvqbn13NYuUPLMyjVmQV/0q+wkn5dpJZgaLFTr6B+0FR+5v4nr739b+p1wpZSnVbwJ/1yLylSoiN3ToYIzql/pj48gTuLH+HN5VqFU6nOKqCkLyIXiki2iGwVkXubOD5eRFaLSJ2IXNboWL1/CcXDyyi2p6KyCpKlVMfon6zIRMIuf5G+jr24PvyVDuFUqpNqNumLiBOYDUwBBgNXisjgRs12AdcDc5t4iypjzAj/Y2oTx9tUZYk1XDMsXodrnizpO5GiEbdzqVnMB/Nn2x2OUqoFArnSHwVsNcbkGGO8wDxgWsMGxpgdxpj1gK8NYjwpNSXWxKzIBE36rSH5kt+wK2ooU7Y/wtbsr+wORyl1ggJJ+qnA7gbbuf59gQoXkSwRWS4i05tqICK3+Ntk5ee3bjnf+jLrSt+pffqtw+km/ppXMeKAf96Aqa22OyKl1AkIJOk3tdTUiQzf6GWMyQSuAp4Ukb5HvZkxLxhjMo0xmSkprTy0ssL/IxKlQzZbS2z3fqw77Xf0q9tCzut32B2OUuoEBJL0c4GeDbbTgL2BfoAxZq//OQf4FBh5AvGdNGdlgfUiKrk9Pzbojb34et6JvIy+O+dT8eUbdoejlApQIEl/JdBfRDJExAPMAAIahSMiCSIS5n+dDIwFNrU02Jbw1BRSixvC49rzY4OewyEMuuaPrPCdgvu9n8J+na2rVGfQbNI3xtQBs4APgM3AfGPMRhF5SESmAojIGSKSC1wOPC8ihzLAKUCWiKwDPgEeMca0a9KPrC2iwp2gC6K3gUE9Ellx+mOU+CKofv0qqCq2OySlVDNcgTQyxiwCFjXa90CD1yuxun0an7cMGHaSMbaYz2eIqSumOjLRrhCC3k1TzuSujffwdNkD+BbciOPqf4LDaXdYSqljCOoZuSVVtSRKKbXhWl2zrUR6XMz43uXcX3s9jm2LYfFv7A5JKXUcQZ30iypqSJaDmEgdudOWJgzsgnf4D3i9fhIs/TOs/6fdISmljiGok35hWQ1JlCIxWmytrT1w8WCe8dzEV66hmHdnwu6VdoeklGpCUCf90pJCwqQOd6wm/bYWH+nhwekjuLZ8FqXuFJh3JRTvtDsspVQjQZ30K/0lGMLjtdhae5gyrDtnDRvAFWU/pb7WC3O/D9UH7Q5LKdVAUCd970GrBENkgpZgaC8PTxtKQURv7vPcgyncCvOvgzqv3WEppfyCOun7yvYD4InTYmvtJSk6jEcuPZV5BRm8l34v5HwCC+8AowuvKNURBHXSp+JQCQYdvdOeJg3uyhWZPZm1eTB7Rv4U1s+Dj35td1hKKYI86buq/MXWInWcfnv7v4tPoUd8BFd+PR7viOtg6ZOw/Dm7w1Iq5AV10vfUFFJ78sWbAAAb9ElEQVTmiAWn2+5QQk5MuJsnrxhBbkkVv6i6FjPwInj/F7B+vt2hKRXSgjrpR9UWUelOsDuMkJWZnshPJg3g7XX7eafvbyH9bHj7Vsh+z+7QlApZQZv0jTHE1JdQ7dGuHTvNnNiP0RmJ3PfvLWw//6/Qfbg1omf753aHplRICtqkX1pdRyKl1EVoHX07OR3CkzNG4HE5mPXmFqqvmA+JGTB3Buxabnd4SoWcoE36RRVeUuQgJlKTvt26x0Xwp8uHs3FvKb9ZvA+ufRdiusHrl2m5BqXaWdAm/eLSUmKlEofW3ekQzjulKzMn9uXvX+5mfnYtXP9vazWz1y+FPavtDk+pkBG0Sb+s0CrB4InT2bgdxV3nD2RM3yTuf2cDG8uj4Lp/QUQ8/G067Flld3hKhYSAkr6IXCgi2SKyVUTubeL4eBFZLSJ1InJZo2PXicgW/+O61gq8OdWH6+7obNyOwukQnrpyJPGRbm57fTUlnq5w/X8gPB5emw67VtgdolJBr9mkLyJOYDYwBRgMXCkigxs12wVcD8xtdG4i8CAwGhgFPCgi7TKG8lDdnZhELbbWkSRHh/HsNaeTd7CamXNXUxuTBje8Z82a/tt3YcdSu0NUKqgFcqU/CthqjMkxxniBecC0hg2MMTuMMesBX6NzLwA+NMYUGWOKgQ+BC1sh7maZcivph+mVfodzWq8Efn/pMJZuLeR3/9kMcalwwyLr+fXvwZaP7A5RqaAVSNJPBXY32M717wvEyZx7UqTSX4JB6+50SJednsZN4zJ4ZdkO/v7lLms0z/WLILkf/H0GbHjT7hCVCkqBJH1pYl+gJRMDOldEbhGRLBHJys/PD/Ctj89VVUANYeCJapX3U63v3imDGD8ghfvf2cDSrQUQnWL18adlwoIbIetlu0NUKugEkvRzgZ4NttOAvQG+f0DnGmNeMMZkGmMyU1Ja58o8vKaIUmcCSFO/O6ojcDkdPHPVSPqmRHPr31aRnVcG4XFwzVvQbxL8+yfw6R+0LLNSrSiQpL8S6C8iGSLiAWYACwN8/w+AySKS4L+BO9m/r81F1hVr3Z1OIDbczcs3nEGEx8kNL3/J/tJq8ETClX+H4VfCp7+Hf/0Y6uvsDlWpoNBs0jfG1AGzsJL1ZmC+MWajiDwkIlMBROQMEckFLgeeF5GN/nOLgIexfjhWAg/597W52PpiasK17k5n0CM+gjnXn8HBqlpueHklZdW1VmXU6c/C2XfD6lfhH1eDt8LuUJXq9AIap2+MWWSMGWCM6WuM+Z1/3wPGmIX+1yuNMWnGmChjTJIxZkiDc+cYY/r5H+3SSVvlrSeRg9Rp0u80hqbGMfvq0/hmfxk3vZpFdW291TV33gPwnT/Blv/Cy1OgdJ/doSrVqQXljNySymqSKMVE6sidzmTCwC786fvD+XJHEbPmrqa23j8C+Iyb4Mp5ULgNXjwX9q23N1ClOrGgTPplJQW4xIdEabG1zmbaiFQemjaUjzYf4J4F6/H5/DdxB1wAP3zfuvqfcyF8/R97A1WqkwrKpF9VYk3McsbolX5n9IMze/PzCwby9po93PtWg8TfbRjctBi6DIJ5V8Fnj+rIHqVOkMvuANpC9UFrrL9Hk36nNXNiP2rqfDy1eAsOEX7/3WE4HAKx3a2x/P/6MXzyO9i/Aab9BcKi7Q5ZqU4hKJN+bVkBABFxmvQ7s59O6o/PZ3jmk604HMJvpw21Er87Ar77PHQdCh89CPnfwBWvW7N5lVLHFZTdO74KK+lHJmhZ5c5MRLh78gBun9CXuSt2cc+b66k7dHNXBMbeaU3kKt8PL07Ufn6lAhCUSZ/KQgCiE3QBlc5ORPj5BQO56/wBLFiVy53z1uCta1DXr+9E+NHnkNTX6uf/7/1QX2tfwEp1cEGZ9B1VRdTgxqH9vEFBRLjzvP7cf/FgFn2Vx82vZVHpbTBDN74n3PA+ZP4Qlj0Fr3wHDubaF7BSHVhQJn13TTEHJVbr7gSZG8dl8IfvDeN/W/K58sUVFJTXfHvQHQ4XPwHfewn2b4Tnzobs9+wLVqkOKiiTvsdbTLkj1u4wVBu44oxePHfN6WTnlfK9Z5exo6BRaYZhl8Etn1m1+f8+Axb9HGqr7AlWqQ4oKJN+RO1BKl1xdoeh2sjkId2Ye/OZlFXXcemzy1i5o1E5p+R+1nj+s2bBly9Ys3j3b7QnWKU6mKBM+lH1B6nWCptB7bReCbx52xjiItxc9eJy5q/cfWQDVxhc8Du4+k2oKIAXJsDSp8BXb0u8SnUUQZn0Y0wptWGa9INdRnIU79w+ltEZSdzz5noe/vemb4d0HtJ/Etz+BfSfDB/eD69eAsU7bIlXqY4g6JK+QwyxpoK6cE36oSAu0s0rN5zB9WPSeWnJdn7w0pfkl9Uc2Sgq2Zq8Nf1Zq1jbX8bAly+Cr/GSzkoFv6BL+lEeK/GbCC2rHCpcTge/njqExy4fzupdxXznqf8d3c8vAiOusq76e42GRT+D16ZCUY49QStlk6BL+rFhTgAcUZr0Q81lp6fxzsyxRHqcXPnCcv7y6dZvi7UdEt/TmsU79WnYt8666l/6Z12ZS4WMgJK+iFwoItkislVE7m3ieJiI/MN/fIWIpPv3p4tIlYis9T+ea93wjxYbbo3Nd8VoWeVQdEr3WBbeMY4LhnTj0fezuealFeQdrD6ykQicdi3cvhz6ngsfPmCVcdi7xp6glWpHzSZ9EXECs4EpwGDgShEZ3KjZjUCxMaYf8ATwhwbHthljRvgft7ZS3McUG2Zd2WmFzdAVG+7mmatG8uj3TmXNrhIu/PPn/Gd9EytuxaXCjDfg+69Z9XtemAj/+RlUH2z/oJVqJ4Fc6Y8CthpjcowxXmAeMK1Rm2nAq/7XC4DzROyZDhvrsZJ+RLzW3QllIsL3z+jJv+8cR6/ESGbOXc3MuaspqvA2bgiDp8GslTDqZsh6CZ7OhHX/0Fr9KigFkvRTgYaDoHP9+5ps419I/SBwqFM9Q0TWiMhnInL2ScbbrDi3NSIjKl6v9BX0TYnmrdvG8PMLBvLfjXlMfuIz/rVuL6ZxQg+Pg4v+CDd/DHFp8PYt1gpd+9bZE7hSbSSQpN/UFXvjS6BjtdkH9DLGjATuAuaKyFH1EUTkFhHJEpGs/Pz8AEI6tjhPHRUmjNgYLcOgLC6ng5kT+/GvO8bRPS6CO/6+hutfXsnuosqjG/cYac3mnfo0FG6F58+xFmwpP9D+gSvVBgJJ+rlAzwbbacDeY7URERcQBxQZY2qMMYUAxphVwDZgQOMPMMa8YIzJNMZkpqSc3BV6nKuWYmKICXef1Puo4DOoWyzvzBzLg5cMJmtHEec/8RlPLd5CdW2jWboOh3Wj945VMPpWWPM6PHUaLHkCaqubfnOlOolAkv5KoL+IZIiIB5gBLGzUZiFwnf/1ZcDHxhgjIin+G8GISB+gP9CmA6PjXV4OSixOh1bYVEdzOoQbxmbw0d3ncO6gLjz+4TdMevwz3vtq39FdPhHxMOURa5RP+jj46NfwTCasm6cTu1Sn1WzS9/fRzwI+ADYD840xG0XkIRGZ6m/2EpAkIluxunEODescD6wXkXVYN3hvNcY0mjXTuuKdNZQ7tNiaOr7ucRH85erTmXvzaKI8Lm57YzXff/4LVu0sPrpxcn+4ah5c+y5EJsLbP4Lnx8OWj/Rmr+p05KirG5tlZmaarKysFp07YcIEXhu3g22Rw5n4q3dbOTIVrOrqfcxbuZsnP9pCQXkNFwzpys8mD6R/15ijG/t8sPEtWPwQlOyEXmPgvPuh95j2D1ypBkRklTEms7l2QTcjN8FRRbU73u4wVCficjq45szefPbzCdx9/gCWbi1k8pOfc8ff17Blf9mRjR0Oq2b/rCy46DGrjMPLU+Bvl8LulfZ8AaVOQFAlfZf4iHFUa4VN1SJRYS7uOK8/n98zkdvO6cvizfuZ/OTnzHxjNV/lNpqw5fJY4/rvXAPnPwT71sJLk+D178HuL+35AkoFIKiSfqzbGoVRH55ocySqM0uM8nDPhYNY8otzue2cvnz+TT6XPLOEa/66gs++yT/yhq8nEsb+GH683kr+e9fAS+fDKxdDzmfa5686nKBK+nFuq2iWidRia+rkHUr+S395LvdOGcQ3+8u4bs6XTHr8M/72xQ4qahoUaQuLtpL/T76CC34PBVusKp5/nQSbFuriLarDCKqkf6gEg2iFTdWKYsPd3HpOX5b84lyeuGI4UWEu7n93I6N/v5j73v6KDXsadP14ouCsmfDjdfCdx6GyAOb/AJ45A7LmgLeJCWFKtSOX3QG0ptgwf4XNaK2wqVqfx+XguyPTmD4ildW7SnhjxU4WrMrljRW7GJoay/dOS2PaiFQSozzgDoczboTTr4fNC2HJk/Dvn8LihyHzBjjjZojtbvdXUiEoqIZs/vG60fw842v+O+VzJo8e3sqRKXW0g5W1vLUmlwWrctm4txSXQ5gwMIVLhvdg0ildiQrzX1cZAzuXwfK/wNf/AYfTKvQ26hboOdoq/KbUSQh0yGZwXel7rFmSEbF6pa/aR1ykmxvGZnDD2Ay+zivlzVW5LFy3l482HyDc7eC8QV2ZPKQrEwd1ITZ9LKSPhcJtsPIlq7zDhjeh2zDIvNEaChrWxNwApVpRUF3pz5t5Ohcl72H7TV8zvKeO1Vf28PkMK3cU8a/1e/lg437yy2pwO4Wz+iZz7sAUzh3UlV5JkVBTDuv/YfX1798AnmgYdrlV96fHSL36Vyck0Cv9oEr67/1kOINjKzF3riU9OaqVI1PqxPl8hjW7i3l/Qx6LNx8gp6ACgL4pUZzdP4XxA5IZnZ5IVP5aq5b/xrehrhq6DoPTfgBDLwMdmKACEJJJ/4ufDSYiMpzev1hOQpSnlSNT6uRtL6hg8eb9fL6lgBU5hdTU+XA5hOE94zmzTyJjUt2cXraY8PWvW7X8HW4YcAEMvxL6nw+uMLu/guqgQjLpb7q3H3nuVM75zadaZVN1eNW19WTtKGbZtgK+yClkfe5B6n0Gh1hloL/TtZDJ3o/J2LcIV1U+hMdbN39P/b5V88cRVCOu1UkKyRu58a5qNvsiNeGrTiHc7WRc/2TG9bcGHlTU1LFqZzFZO4tZtbOI2ZvC+aP3ApxM4oLwzVzlWM6otfPxrH6VuqhuOIZOxzH0e5CaqT8AKmDBlfQd1RTXRtgdhlItEhXmYvyAFMYPsBYSqvcZthwoY+2uEtblpvOHPWPZdTCf8SaLi+uXM2H5S4SteI4SVwq53c6ldsDFJA+eQGpiDA698FHHEDxJ31tJpKOWg3XhdkeiVKtwOoRB3WIZ1C2WGaN6AeCt8/HN/nP5Oq+MJ3P3Eb3zQwYXf8qZu98mIvcfFC+OZqEZyVfRYynsNo7uXVJIT4qkV2IUvZIi6RYbrn8Jh7jgSfo1Zeyvj6VQk74KYh6Xg6GpcQxNjYPT04AzgF9RXFxM7rpFuLa8xwV5nzK98n/U5rj4cusgPq4fwQu+4WwzPXA7HXSPiyA1PoLUhAh6xEfQPS6c7nHhdIsLp2tMOPGRbkSHiwatoLqR2//OV3FXFbHpxZ+2clRKdSL1dbB7OXzzAeabD5CCbADKI3qwJWY0q1wj+Nw7iOxSFwfKao4qBOpxOegSE0ZKTBgp0WEkx4SRHB1GcrSHxCgPiZEeEqM9JER6iI90E+Zy2vAlVWOteiNXRC4E/gw4gb8aYx5pdDwMeA04HSgErjDG7PAf+yVwI1AP3GmM+eAEvscJ8bnCcdTrwtUqxDld1pq+6eOQyQ9D8Q7YupjobR8zMucjRnrf5iYEeoyg/vSzKe4yml3Rp7Kvys3+0mr2l1aTX1bDgbIadhZWsmpnMUWV3mNWiY70OEmI9BAb4SYuwkVchJvYcDexEW5iwl3EhLuJCXMRE+4iKsxFdLiL6DAXkR6n/9mFx6U3ottLs0nfv7D5bOB8IBdYKSILjTGbGjS7ESg2xvQTkRnAH4ArRGQw1kLqQ4AewEciMsAY0yZ1Zn2uMBx1mvSVOkJCulX87Ywbob4W9qyGnE8h51OcK54l2fcUyeKE7sOtZR/7jLHqAUV9W86k3mcoqvBSXOmlsNx7+HVJpZfiylpKKms5WFXLwSovOwoqKa2upbSqlgpvYP+ruxxChMdJpMdJhNtJhMdFhNtBhMdJuMtJuNtJmNtBuNvaDnM7CHM5CHM58bis14efndZrt9N6eFxy+LXbKbgcDtwuB26H4HI6cDoEt1OsZ4cj6G+CB3KlPwrYaozJARCRecA0oGHSnwb82v96AfCMWJ2C04B5xpgaYLt/4fRRwBfH+rDs7GwmTJhwgl8DfOLEjL6L/bu3t+h8pUJRmGMwQ+PKGRFfzqnF2ZySuwbPF88AsLsyjI2lUWwqjWJTaSTbKyKoN4ElxDAgBUhG8Dk9+FxhGIcHnzMMn8uDcbjxOT0Ypwefw4NxuvE53FQ43ZQ73PgcLozTg3G4MA6XtX344caI0ypa1xaMD4wPMT4wBuHQtvl2P8Y61sRrMca/r8G2f9+3x/zbGDBY+zC4q4pJyF3SNt/LL5CknwrsbrCdC4w+VhtjTJ2IHASS/PuXNzo3tfEHiMgtwC0AYWEtm3Hoc1k3cE1NBehkXKUCUuNzsKo4llXFsQB4HD4GxlQyJLaCoXEVjE4s5cJuRQBU1wtbyiPJLrMe35RFsLsyHB/H/iEQDM76Gpz1Na0euwHrR0Cc/h8DJ4jTv/3ts7XP0WDb/1oc/tcOfxs53BYR/zlWG5AG7QX4tg0iGMRq1/A1gjnUxnFov4D/+dttDu9rD4Ek/aYiady7d6w2gZyLMeYF4AWwbuR++umnAYR1JJ/PcM75FyLhPj79+MTPV0o1wRjrnsCeVYTnZjFs7xqG5a2H2nzruCsCug62KoV2HQpdh0CXwRChBQ9b7u4WnRXoiKtAkn4u0LPBdhqw9xhtckXEBcQBRQGe2yocDmmTqwmlQpoIJGZYj2GXWft89ZCfDXnrYd9663njO7DqlW/Pi+kOKYOsR3J/SBkIyQMgKkWrh9oskKS/EugvIhnAHqwbs1c1arMQuA6rr/4y4GNjjBGRhcBcEXkc60Zuf+DL1gpeKWUDh9O6uu86GIbPsPYZA2X7IG8D5G+GA/7H6lehtsESkWGxkNQXEvtaPyQJ/h+U+N7WD4WWk2hzzSZ9fx/9LOADrCGbc4wxG0XkISDLGLMQeAn4m/9GbRHWDwP+dvOxbvrWATPbauSOUspGIhDbw3oMmPztfp8PSvdAQba1eEzhVmvR+NyVsPEt/81OP6cH4tIgrqf/kQZxqf73TYWYblbROf1L4aQE1eSsQ6N2WnJPQCnVzuq8cHA3FG+H4p1QshNKdsHBXOtRlsdRtwBd4RDd1f/oYj2iuljDS6OSITLJekQkQkSCtVZxiAjJKptKqU7E5bG6epL6Nn28zgvleXBwj/XXQvl+qwupLA/KD1h/OexcBlXFNDE+xP8ZEVbyD4/79hEWA+Gx1kplYTH+52hwR4Inynp2R1g/MO5Iaw0DV7j17PRYj07cDaVJXynVMbk8EN/LehxPfR1UFUFFPlQWQWWh9agqhuoS//NB61GeB4VboLoUasqgpYM/HC5rgRun2//aZd3rEKf1g+AfBvrtEM2GPxLfjtPHP/Yf/5wAug2DK+e2LKYAadJXSnVuTte3XT0nqr4WvOXWesW1ldbDWwm1VVBX5X+usZawrKuBeq/1qKsBX631g+OrtUY0+eq+TeC+ejD1DZK7jyNGsDcYr2/9OPh/KBIzWukfyrFp0ldKhS6n2+r+iUiwO5J203k7ppRSSp0wTfpKKRVCNOkrpVQI0aSvlFIhRJO+UkqFEE36SikVQjTpK6VUCNGkr5RSIaTDFVwTkXxg50m8RTJQ0ErhdBah9p1D7fuCfudQcTLfubcxJqW5Rh0u6Z8sEckKpNJcMAm17xxq3xf0O4eK9vjO2r2jlFIhRJO+UkqFkGBM+i/YHYANQu07h9r3Bf3OoaLNv3PQ9ekrpZQ6tmC80ldKKXUMQZP0ReRCEckWka0icq/d8bQ1EekpIp+IyGYR2SgiP7Y7pvYiIk4RWSMi/7Y7lvYgIvEiskBEvvb/+z7L7pjamoj81P/f9QYR+buIBN1ityIyR0QOiMiGBvsSReRDEdnif271Qv9BkfRFxAnMBqYAg4ErRWSwvVG1uTrgbmPMKcCZwMwQ+M6H/BjYbHcQ7ejPwPvGmEHAcIL8u4tIKnAnkGmMGQo4gRn2RtUmXgEubLTvXmCxMaY/sNi/3aqCIukDo4CtxpgcY4wXmAdMszmmNmWM2WeMWe1/XYaVCFLtjartiUga8B3gr3bH0h5EJBYYD7wEYIzxGmNK7I2qXbiACBFxAZHAXpvjaXXGmM+Boka7pwGv+l+/Ckxv7c8NlqSfCuxusJ1LCCTAQ0QkHRgJrLA3knbxJHAP4LM7kHbSB8gHXvZ3af1VRKLsDqotGWP2AI8Bu4B9wEFjzH/tjarddDXG7APrwg5owcK/xxcsSV+a2BcSw5JEJBp4E/iJMabU7njakohcDBwwxqyyO5Z25AJOA541xowEKmiDP/k7En8/9jQgA+gBRInINfZGFTyCJennAj0bbKcRhH8ONiYibqyE/4Yx5i2742kHY4GpIrIDqwvvXBF53d6Q2lwukGuMOfRX3AKsH4FgNgnYbozJN8bUAm8BY2yOqb3sF5HuAP7nA639AcGS9FcC/UUkQ0Q8WDd9FtocU5sSEcHq591sjHnc7njagzHml8aYNGNMOta/44+NMUF9BWiMyQN2i8hA/67zgE02htQedgFnikik/7/z8wjym9cNLASu87++Dni3tT/A1dpvaAdjTJ2IzAI+wLrTP8cYs9HmsNraWOAHwFcista/71fGmEU2xqTaxh3AG/4LmhzgBpvjaVPGmBUisgBYjTVKbQ1BODtXRP4OTACSRSQXeBB4BJgvIjdi/fhd3uqfqzNylVIqdARL945SSqkAaNJXSqkQoklfKaVCiCZ9pZQKIZr0lVIqhGjSV0qpEKJJXymlQogmfaWOQUT+JCKbRORFEfnMX8IbEUkTkSv8rz0i8rm/GqRSHZ4mfaWaICJ9gLHGmMHAWuAtY0y9//B5+Ovf+Et5LwausCVQpU6QzshVqhF/nZuPsMqU5Pl3f9cYs0NExmHVQykByoDvArHA/zPGXGRHvEqdCE36SjVBRH4L7ABeA3YZY7o1OPY+8DNjzAb/thPIM8ak2BGrUidCu3eUatowYB2QjHVV39BAIPvQhr/bxysiMe0XnlItozeflGraEGAjEAYcXpRbRJKwVnKqbdQ+DKhuv/CUahm90leqEf8Ve60xptIYUww4ReRQ4s+g0QI9/h+C/CZ+CJTqcDTpK3W0ocCGBtv/Bcb5X3+NVf98g4gcWs1pIqDrGKhOQW/kKtUMERkJ3GWM+cExjr8F/NIYk93UcaU6Er3SV6oZxpg1wCeHJmc15F/N6h1N+Kqz0Ct9pZQKIXqlr5RSIUSTvlJKhRBN+kopFUI06SulVAjRpK+UUiFEk75SSoUQTfpKKRVC/j8knUuV7qxh0gAAAABJRU5ErkJggg\u003d\u003d\n"
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": "t_data \u003d np.linspace(0, 10, 101)\nf_data \u003d model(t\u003dt_data).f\nf_re_func \u003d lambda x: - np.exp(- x * s_data[1:]).dot(ans.c) / fit_result.value(a)**2\nf_re \u003d [f_re_func(t_i) for t_i in t_data]\n\nplt.axhline(0, color\u003d\u0027black\u0027)\nplt.axvline(0, color\u003d\u0027black\u0027)\nplt.plot(t_data, f_data, label\u003dr\u0027$f(t)$\u0027)\nplt.plot(t_data, f_re, label\u003dr\u0027$f_{re}(t)$\u0027)\nplt.xlabel(r\u0027$t$\u0027)\nplt.xlabel(r\u0027$f(t)$\u0027)\nplt.legend()",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%%\n",
          "is_executing": false
        }
      }
    },
    {
      "cell_type": "markdown",
      "source": "Not bad, for an ill-defined problem.\n",
      "metadata": {
        "pycharm": {
          "metadata": false,
          "name": "#%% md\n"
        }
      }
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.1"
    },
    "stem_cell": {
      "cell_type": "raw",
      "source": "",
      "metadata": {
        "pycharm": {
          "metadata": false
        }
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 1
}