src/calibration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spike Calibration\n",
    "The accuracy of the double spike technique depends crucially on how well the double spike is calibrated. This notebook illustrates one method of calibration which uses measurements of mixtures of the double spike with a known standard material. See Figure 6 of [Klaver and Coath 2019, Geostand. Geoanal. Res.](https://doi.org/10.1111/ggr.12248) and Appendix F of the manuscript for further details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import doublespike as ds\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example we will consider Fe isotopes, with a $^{57}$Fe-$^{58}$Fe double spike."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "isodata_fe = ds.IsoData(\"Fe\")\n",
    "true_spike = [0.0025, 0.1060, 0.4413, 0.4502]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo simulated data\n",
    "We will use Monte Carlo simulation to produce some fictious data to use for the calibration. We first measure the pure double spike, taking an average over 100 integrations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.02382983 1.04017915 4.39233518 4.54365706]\n"
     ]
    }
   ],
   "source": [
    "n = 100\n",
    "spike_measurements = ds.monterun(isodata_fe, 1.0, true_spike, alpha=0.0, beta=0.8, n=n)\n",
    "spike_measurement = np.mean(spike_measurements, axis=0)\n",
    "print(spike_measurement)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will take a series of mixtures of spike and standard, again averaging over 100 integrations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.38272341 6.65847944 1.50674743 1.45205632]\n",
      " [0.28380651 5.0328508  2.33539185 2.34797347]\n",
      " [0.17637531 3.39969452 3.16088875 3.26305085]]\n"
     ]
    }
   ],
   "source": [
    "mixture_props = [0.3, 0.5, 0.7]\n",
    "mixture_betas = [2.0, 1.5, 1.7]\n",
    "mixture_measurements = [\n",
    "    ds.monterun(isodata_fe, prop, true_spike, alpha=0.0, beta=beta, n=n)\n",
    "    for prop, beta in zip(mixture_props, mixture_betas)\n",
    "]\n",
    "mixture_measurement = [np.mean(m, axis=0) for m in mixture_measurements]\n",
    "mixture_measurement = np.vstack(mixture_measurement)\n",
    "print(mixture_measurement)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing the calibration\n",
    "We now perform the calibration using both the pure spike measurement and the spike-standard mixture measurements as input. The calibration routine minimises the chi-squared misfit between data and model-predicted isotopic ratios, i.e. minimising\n",
    "\n",
    "$$ \\chi^2 = \\sum (\\mathbf{m}_\\text{obs} - \\mathbf{m}_\\text{model})^T \\; \\mathbf{W}\\; (\\mathbf{m}_\\text{obs} - \\mathbf{m}_\\text{model}) $$\n",
    "where the weight matrix $\\mathbf{W}$ is chosen as the inverse of the expected covariance matrix for `errormodel['measured']`. It may be possible to improve the calibration by adjusting the weighting to be more appropriate for the data of interest.\n",
    "\n",
    "The calibration routine returns the composition of the calibrated spike, estimates of the mixing proportions of spike and standard, the instrumental fractionation factors for both the spike measurements and the spike-standard mixture measurements, the chi-squared misfit, and the degrees of freedom for the chi-squared statistic. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'calibrated_spike': array([0.00250002, 0.10600085, 0.44130018, 0.45019895]), 'prop_mixture': [0.300000764129837, 0.5000000274237444, 0.7000021519071813], 'beta_mixture': array([2.00001085, 1.49999473, 1.69989461]), 'beta_spike': array([0.80041034]), 'misfit': 0.013707250567584737, 'df': 2}\n"
     ]
    }
   ],
   "source": [
    "cal = ds.spike_calibration(isodata_fe, spike_measurement, mixture_measurement)\n",
    "print(cal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since in this example we know the true composition of the spike, we can see how well the calibration routine has recovered the true composition of the spike:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True spike: [0.0025, 0.106, 0.4413, 0.4502]\n",
      "Calibrated spike: [0.00250002 0.10600085 0.44130018 0.45019895]\n",
      "Difference: [ 2.11636760e-08  8.49068513e-07  1.78117189e-07 -1.04834938e-06]\n"
     ]
    }
   ],
   "source": [
    "calibrated_spike = cal[\"calibrated_spike\"]\n",
    "print(\"True spike:\", true_spike)\n",
    "print(\"Calibrated spike:\", calibrated_spike)\n",
    "print(\"Difference:\", calibrated_spike - true_spike)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Validation of the calibration\n",
    "One way of checking the calibration is to run the double spike inversion on data produced from a series of mixtures of spike and standard in different proportions. Here we use Monte Carlo simulation to produce more fictious data for a series of such mixtures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 1000\n",
    "mixture_props = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])\n",
    "mixture_betas = 1.1 * np.ones_like(mixture_props)\n",
    "measured = [\n",
    "    ds.monterun(isodata_fe, prop, true_spike, alpha=0.0, beta=beta, n=n)\n",
    "    for prop, beta in zip(mixture_props, mixture_betas)\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we run this fictious data through the double spike inversion to estimate the natural fractionation factor $\\alpha$. The true value of $\\alpha$ should be zero as we are measuring mixtures of the standard and the spike."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.00026161827970858553, 9.780590101979544e-05, 8.761077848548851e-05, 7.224993386801511e-05, 0.00012707427804734892, 6.668255672932258e-05, -9.313806538432102e-05, 0.00023519598247875235, 0.0003340033239492854]\n"
     ]
    }
   ],
   "source": [
    "data = [ds.dsinversion(isodata_fe, m, calibrated_spike) for m in measured]\n",
    "alpha = [np.mean(d[\"alpha\"]) for d in data]\n",
    "prop = [np.mean(d[\"prop\"]) for d in data]\n",
    "alpha_err = [np.std(d[\"alpha\"]) for d in data]\n",
    "print(alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the recovered $\\alpha$ values as a function of the proportion of spike in the mixture. Here there is no systematic trend in mean $\\alpha$ with proportion of spike. One way of identifying a poorly-calibrated double spike is to see a systematic variation of $\\alpha$ with proportion of spike in the mixture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAEGCAYAAAC6i5gfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAm40lEQVR4nO3de5wcVZ338c+XRFhULgnEGBIwQQK7ATXAyOVZVDRIIq99THS5hF0kKHJR8LLeFh59AGGzgKjsslwUMUtgkQRBYEQwhgAqPFwyECCAxgwQJDHAbMJ9IZDwe/6oM6TSdGe6Z6arK5nv+/Xq11SdOnXqV109/euqrj5HEYGZmVmRNml1AGZmNvA4+ZiZWeGcfMzMrHBOPmZmVjgnHzMzK9zgVgewIdh2221j9OjRrQ7DzGyDcu+99/53RAyrtszJpw6jR4+mo6Oj1WGYmW1QJD1Ra5kvu5mZWeGcfMzMrHBOPmZmVjgnHzMzK5yTj5mZFc7Jx8zMCufkY2ZmhStF8pE0SdIiSZ2STqqyfDNJs9PyuyWNTuXbSLpV0kuSzq9YZ09JC9M650lSKh8qaa6kxenvkEJ20szM3tTy5CNpEHAB8AlgHHC4pHEV1Y4Gno2InYBzgbNT+avA/wW+UaXpi4BjgLHpMSmVnwTMi4ixwLw0b2ZmFQ778Z0c9uM7m9J2y5MPsBfQGRGPRcRrwCxgckWdycDMNH01MEGSIuLliLidLAm9SdIIYMuIuCuy0fIuA6ZUaWtmrtzMzApShuQzEngyN780lVWtExGrgeeBbXpoc2mNNodHxPI0/RQwvFoDko6V1CGpo6urq579MDOzOpUh+bRMOiuqOo54RFwcEW0R0TZsWNV+8czMrJfKkHyWAdvn5kelsqp1JA0GtgJW9NDmqBptPp0uy3Vfnnum15GbmVmvlCH5zAfGShojaVNgKtBeUacdmJamDwZuSWctVaXLai9I2ifd5XYkcH2Vtqblys3MrCAtH1IhIlZLOhGYAwwCZkTEw5JOBzoioh34KXC5pE5gJVmCAkDSEmBLYFNJU4ADI+IR4IvApcDmwE3pAXAWcJWko4EngEObvpNmZraOlicfgIi4EbixouyU3PSrwCE11h1do7wD2K1K+QpgQh/CNTOzPirDZTczMxtgnHzMzKxwTj5mZlY4Jx8zMyuck4+ZmRXOycfMzArn5GNmZoVz8jEzs8I5+ZiZWeGcfMzMrHBOPmZmVjgnHzMzK5yTj5mZFc7Jx8zMCufkY2ZmhXPyMTOzwpUi+UiaJGmRpE5JJ1VZvpmk2Wn53ZJG55adnMoXSZqYynaRdH/u8YKkr6Zlp0lallt2UFH7aWZmmZaPZCppEHAB8HFgKTBfUnsaCrvb0cCzEbGTpKnA2cBhksaRDam9K7AdcLOknSNiETA+1/4y4Npce+dGxPebvGtmZlZDGc589gI6I+KxiHgNmAVMrqgzGZiZpq8GJkhSKp8VEasi4nGgM7WXNwF4NCKeaNoemJlZQ8qQfEYCT+bml6ayqnUiYjXwPLBNnetOBa6sKDtR0oOSZkgaUi0oScdK6pDU0dXV1cj+mJlZD8qQfJpG0qbAJ4Gf54ovAt5LdlluOfCDautGxMUR0RYRbcOGDWt2qGZmA0oZks8yYPvc/KhUVrWOpMHAVsCKOtb9BHBfRDzdXRART0fEmoh4A/gJb71MZ2ZmTVaG5DMfGCtpTDpTmQq0V9RpB6al6YOBWyIiUvnUdDfcGGAscE9uvcOpuOQmaURu9lPAQ/22J2ZmVpeW3+0WEaslnQjMAQYBMyLiYUmnAx0R0Q78FLhcUiewkixBkepdBTwCrAZOiIg1AJLeQXYH3XEVm/yepPFAAEuqLDczsyZrefIBiIgbgRsryk7JTb8KHFJj3enA9CrlL5PdlFBZ/pm+xtuIw358JwCzj9u3yM2amZVaGS67mZnZAOPkY2ZmhXPyMTOzwjn5mJlZ4Zx8zMyscE4+ZmZWOCcfMzMrnJOPmVmLHfbjO9/8TeBA4eRjZmaFc/IxM7PCOfmYmVnhnHzMzKxwTj5mZlY4Jx8zMyuck4+ZmRXOycfMzApXiuQjaZKkRZI6JZ1UZflmkman5XdLGp1bdnIqXyRpYq58iaSFku6X1JErHypprqTF6e+Qpu+gmZmto+XJR9Ig4ALgE8A44HBJ4yqqHQ08GxE7AecCZ6d1x5ENqb0rMAm4MLXX7aMRMT4i2nJlJwHzImIsMC/Nm5lZgVqefIC9gM6IeCwiXgNmAZMr6kwGZqbpq4EJkpTKZ0XEqoh4HOhM7a1Pvq2ZwJS+74KZmTWiDMlnJPBkbn5pKqtaJyJWA88D2/SwbgC/kXSvpGNzdYZHxPI0/RQwvFpQko6V1CGpo6urq/G9MjOzmsqQfJplv4jYg+xy3gmSPlxZISKCLEm9RURcHBFtEdE2bNiwJodqZjawlCH5LAO2z82PSmVV60gaDGwFrFjfuhHR/fcZ4FrWXo57WtKI1NYI4Jl+3BczM6tDGZLPfGCspDGSNiW7gaC9ok47MC1NHwzcks5a2oGp6W64McBY4B5J75C0BYCkdwAHAg9VaWsacH2T9svMzGoY3OoAImK1pBOBOcAgYEZEPCzpdKAjItqBnwKXS+oEVpIlKFK9q4BHgNXACRGxRtJw4NrsngQGAz+LiF+nTZ4FXCXpaOAJ4NDCdrZEuscOmX3cvi2OxMwGopYnH4CIuBG4saLslNz0q8AhNdadDkyvKHsM+ECN+iuACX0M2czM+qAMl93MzGyAcfIxswFjIA5XXVZOPmZmVjgnHzMzK5yTj5VKWS+LlDWusvLzZT1x8jEzs8I5+ZhtwHyGYRsqJx8zMyuck4+ZmRXOycfMzArn5GNmZoWrO/lI2rOZgZiZ2cDRyJnPNElXS9qnu0DSD5sQk5mZbeQaST7PAG3ALyQtkvQ48O7mhGVmZhuzRoZUOALYJSJWSdoOOBNY0JywzMxsY9bImc+TwBiAiPhLREwDjmtKVGZmtlFrJPl8BbhG0uWSvibpB8DL/RGEpEnpUl6npJOqLN9M0uy0/G5Jo3PLTk7liyRNTGXbS7pV0iOSHpb0lVz90yQtk3R/ehzUH/tgZmb1q/uyW0Q8ImkP4ABgPPAUMLmvAUgaBFwAfBxYCsyX1B4Rj+SqHQ08GxE7SZoKnA0cJmkc2ZDauwLbATdL2plsSO2vR8R9krYA7pU0N9fmuRHx/b7GbmZmvdPQMNoRsQr4VXr0l72AzjT0NZJmkSW1fPKZDJyWpq8GzpekVD4rxfW4pE5gr4i4E1ieYn5R0h+AkRVtmplZi5ThR6Yjyb5P6rY0lVWtExGrgeeBbepZN12i2x24O1d8oqQHJc2QNKRaUJKOldQhqaOrq6vhnTIzs9rKkHyaRtI7gWuAr0bEC6n4IuC9ZJcOlwM/qLZuRFwcEW0R0TZs2LAiwjUzGzB6vOwm6WvrWx4Rff2h6TJg+9z8qFRWrc5SSYOBrYAV61tX0tvIEs8VEfGLXLxPd09L+glwQx/jNzOzBtVz5rNFerQBXyC7rDUSOB7Yox9imA+MlTRG0qZkNxC0V9RpB6al6YOBWyIiUvnUdDfcGGAscE/6PuinwB8qk6OkEbnZTwEP9cM+mJlZA3o884mI7wJI+h2wR0S8mOZPox9uPIiI1ZJOBOYAg4AZEfGwpNOBjohoJ0skl6cbClaSJShSvavIbiRYDZwQEWsk7Qd8Blgo6f60qf8TETcC35M0HghgCf6tkplZ4Rq522048Fpu/rVU1mcpKdxYUXZKbvpV4JAa604HpleU3Q6oRv3P9DVeMzPrm0aSz2Vkl7SuTfNTgEv7OyAzM9v4NfIj0+mSbgI+lIo+GxHu283MzBrW6I9M7wPua1IsZmY2QDQymJwkHSHplDS/g6S9mheamZltrBr5kemFwL7A4Wn+RbI+2czMzBrSyGW3vSNiD0kLACLi2fS7HDMzs4Y0cubzeuqBOgAkDQPeaEpUZma2UWsk+ZwHXAu8S9J04HbgX5sSlZmZbdTquuyWuqv5HXAvMIHsB5xTIuIPTYzNzMw2UnUln4gISTdGxPuAPzY5JjMz28g1ctntPkkfbFokZmY2YDR0txtwhKQlwMtkl94iIt7fjMDMzGzj1Ujymdi0KMzMbEBp5LLbn8n6dZsWEU+Q3XLdL71am5lZuVy3YBkL/vwcdz++kr896xauW1A5xmffuIeDJmr2wTOzDV8Z3yeuW7CMk3+xkNfWZD/lXPbcK5z8i4X9GlsjyWfviDgBeBWyHg4A93BQQxEHr7fK+GK3xvk4bvjK+j5xzpxFvPL6mnXKXnl9DefMWdRv2yhFDweSJklaJKlT0klVlm8maXZafrek0bllJ6fyRZIm9tRmGq777lQ+u1ldBBVx8HqjrC92KO+baRnj8nHcOJT1feIvz73SUHlvtLyHg5TQLgA+AYwDDpc0rqLa0cCzEbETcC5wdlp3HNmQ2rsCk4ALJQ3qoc2zgXNTW8+mtvtdEQevN8r6Yi/rm2lZ4/JxbFwZk2JZ3ye223rzhsp7o+7kExFXAN8CzgSWk/Vw8PN+iGEvoDMiHouI14BZwOSKOpOBmWn6amBC6nVhMjArIlZFxONAZ2qvaptpnY+lNkhtTumHfXiLIg5eb5T1xV7WN9OyxuXj2JiyJsWyvk98c+IubP62QeuUbf62QXxz4i79to1GxvP5GvBiRFwQEef3Y9c6I4Enc/NLU1nVOhGxGnge2GY969Yq3wZ4LrVRa1sASDpWUoekjq6uroZ3qoiD1xtlfbGX9c20rHH5ODamrEmxrO8TU3YfyZmffh+bDspSxMitN+fMT7+PKbtXfbvslUYuu20B/EbS7yWdKGmjvs06Ii6OiLaIaBs2bFjD6xdx8HqjrC/2sr6ZljUuH8fGlDUplvV9ArLYdt9ha/YeM5Q7TvpYv8fUyGW370bErsAJwAjgt5Ju7ocYlgHb5+ZHpbKqdSQNBrYCVqxn3VrlK4CtUxu1ttVvmn3wehtTGV/sZX0zLWtcPo6NKWtShHK+TxShkTOfbs8AT5G9kb+rH2KYD4xNd6FtSnYDQXtFnXZgWpo+GLglIiKVT013w40BxgL31GozrXNraoPU5vX9sA8blDK+2Mv6ZlrWuMDHsRFlTYoDWd3d60j6InAoMAz4OXBMRDzS1wAiYrWkE4E5wCBgRkQ8LOl0oCMi2oGfApdL6gRWkiUTUr2rgEeA1cAJEbEmxfuWNtMm/xmYJelfgAWpbSuBKbuP5Mp7/gzA7OP2bXE0a5U1rrIq4/PVnfy+dfWDvLbmDUZuvTnfnLhLy5PiQNZI327bA1+NiPv7O4iIuBG4saLslNz0q8AhNdadDkyvp81U/hjZ3XBmNoCUMSkOZHUnn4g4WdIH0hkFwO8j4oEmxWVmZhuxRm61/jJwBdn3PO8C/kvSl5oVmJmZbbwauez2ebL+3V4GkHQ2cCfwH80IzMzMNl6N3O0mIP8rrTWpzMzMrCGNnPn8J3C3pGvT/BRgRr9HZGZmG71Gbjj4oaTbgP1S0WcjYkFTojIzs41aI7/zmQl8JSLuS/NDJM2IiM81LTozM9soNfKdz/sj4rnumTSY3O79HpGZmW30Gkk+m0ga0j0jaSiNfWdkZmYGNJY8fgDcKal7DJ9DqNKzgJmZWU8aueHgMkkdZIOxAXy6P/p2MzOzgaehy2Yp2TjhmJlZnzTSvY4kHSHplDS/gyR30GlmZg1r5IaDC4F9gcPT/IvABf0ekZmZbfQauey2d0TsIWkBZLdap4HazMzMGtLImc/rkgYBASBpGPBGU6IyM7ONWiPJ5zzgWuBdkqYDtwP/2peNSxoqaa6kxenvkBr1pqU6iyVNy5XvKWmhpE5J50lSKj9H0h8lPSjpWklbp/LRkl6RdH96/Kgv8ZuZWe/UnXwi4grgW8CZwF+Ag4G+XnY7CZgXEWOBeWl+HenHrKcCe5ONQHpqLkldBBwDjE2PSal8LrBbRLwf+BNwcq7JRyNifHoc38f4zcysF3pMPpK2lHSypPOBHchuPNgE+CVwaB+3PxmYmaZnkvWUXWkiMDciVqYufeYCkySNALaMiLsiIoDLutePiN9ExOq0/l3AqD7GaWZm/aieM5/LgV2AhWQDyt1KdtYzJSIm93H7wyNieZp+Chhepc5I4Mnc/NJUNjJNV5ZX+hxwU25+jKQFkn4r6UO1ApN0rKQOSR1dXV117IqZmdWrnrvddoyI9wFIugRYDuwQEa/WswFJNwPvrrLo2/mZiAhJUU+b9ZL0bWA12fDfsDb2FZL2BK6TtGtEvFC5bkRcDFwM0NbW1q9xmZkNdPUkn9e7JyJijaSl9SaetM4BtZZJelrSiIhYni6jPVOl2jJg/9z8KOC2VD6qonxZru2jgL8DJqTLckTEKmBVmr5X0qPAzkBHvftjZmZ9V89ltw9IeiE9XgTe3z0t6S1nDA1qB7rvXpsGXF+lzhzgwDR+0BDgQGBOulz3gqR90l1uR3avL2kS2c0Rn4yI/+luSNKwdLs4knYku0nhsT7ug5mZNajHM5+IGNTE7Z8FXCXpaOAJ0g0MktqA4yPi8xGxUtIZwPy0zukRsTJNfxG4FNic7Hud7u92zgc2A+amu6/vSne2fRg4XdLrZL9ROj7XlpmZFaSl4/FExApgQpXyDrKbG7rnZwAzatTbrUr5TjW2dw1wTR9CNjOzftDIj0zNzMz6hZOPmZkVzsnHzMwK5+RjZmaFc/IxM7PCtfRuNzPrm9nH7dvqEMx6xWc+ZmZWOCcfMzMrnC+7WamU9TJSWeMy21A5+ZhZv3Oytp74spuZmRXOycfMzArn5GNmZoVz8jEzs8L5hoMByl8Im1krOfmY2YDhD13l0dLLbpKGSporaXH6O6RGvWmpzmJJ03Lle0paKKlT0nlpOG0knSZpmaT70+Og3Donp/qLJE1s/l6amVmlVn/ncxIwLyLGAvPS/DokDQVOBfYG9gJOzSWpi4BjgLHpMSm36rkRMT49bkxtjQOmArumuhdKauYw4WZmVkWrk89kYGaanglMqVJnIjA3IlZGxLPAXGCSpBHAlhFxV0QEcFmN9Su3NysiVkXE40AnWUIzM7MCtTr5DI+I5Wn6KWB4lTojgSdz80tT2cg0XVne7URJD0qakTtTqtXWW0g6VlKHpI6urq66d8jMzHrW9OQj6WZJD1V5TM7XS2cv0U+bvQh4LzAeWA78oNEGIuLiiGiLiLZhw4b1U1hmZgYF3O0WEQfUWibpaUkjImJ5uoz2TJVqy4D9c/OjgNtS+aiK8mVpm0/ntvET4IZcW9tXW8fMzIrT6stu7UD33WvTgOur1JkDHChpSLp8diAwJ12ue0HSPukutyO710+JrNungIdy25sqaTNJY8huUrinv3fKzMzWr9W/8zkLuErS0cATwKEAktqA4yPi8xGxUtIZwPy0zukRsTJNfxG4FNgcuCk9AL4naTzZZbwlwHEAEfGwpKuAR4DVwAkRsaape2hmZm/R0uQTESuACVXKO4DP5+ZnADNq1NutSvln1rPN6cD0XoZsZmb9oNWX3czMbABy8jEzs8K1+jsfM7MBbyD2OeczHzMzK5yTj5mZFc7Jx8zMCufkY2ZmhXPyMTOzwjn5mJlZ4Zx8zMyscP6dT5MNxPv3zcx64jMfMzMrnJOPmZkVzsnHzMwK5+RjZmaFc/IxM7PCtTT5SBoqaa6kxenvkBr1pqU6iyVNy5XvKWmhpE5J56XhtJE0W9L96bFE0v2pfLSkV3LLflTIjpqZ2TpafeZzEjAvIsYC89L8OiQNBU4F9gb2Ak7NJamLgGOAsekxCSAiDouI8RExHrgG+EWuyUe7l0XE8c3ZLTMzW59WJ5/JwMw0PROYUqXORGBuRKyMiGeBucAkSSOALSPirogI4LLK9dOZ0KHAlc0J38zMeqPVyWd4RCxP008Bw6vUGQk8mZtfmspGpunK8rwPAU9HxOJc2RhJCyT9VtKHagUm6VhJHZI6urq66twdMzOrR9N7OJB0M/DuKou+nZ+JiJAU/bz5w1n3rGc5sENErJC0J3CdpF0j4oXKFSPiYuBigLa2tv6Oy8xsQGt68omIA2otk/S0pBERsTxdRnumSrVlwP65+VHAbal8VEX5slzbg4FPA3vmYlkFrErT90p6FNgZ6Ghsr8zMrC9afdmtHei+e20acH2VOnOAAyUNSTcaHAjMSZfrXpC0T/pu58iK9Q8A/hgRb16akzRM0qA0vSPZTQqP9fdOmZnZ+rU6+ZwFfFzSYrJkcRaApDZJlwBExErgDGB+epyeygC+CFwCdAKPAjfl2p7KW280+DDwYLr1+mrg+FxbZmZWkJb2ah0RK4AJVco7gM/n5mcAM2rU261G20dVKbuG7NZrMzNroVaf+ZiZ2QDk5GNmZoVz8jEzs8I5+ZiZWeGcfMzMrHBOPmZmVjgnHzMzK5yTj5mZFc7Jx8zMCufkY2ZmhXPyMTOzwjn5mJlZ4Zx8zMyscE4+ZmZWOCcfMzMrXEvH8zEzs/Kafdy+TWu7pWc+koZKmitpcfo7pEa9aanOYknTcuXTJT0p6aWK+ptJmi2pU9Ldkkbnlp2cyhdJmti0nTMzs5pafdntJGBeRIwF5qX5dUgaCpwK7A3sBZyaS1K/TGWVjgaejYidgHOBs1Nb48iG194VmARcKGlQv+6RmZn1qNXJZzIwM03PBKZUqTMRmBsRKyPiWWAuWeIgIu6KiOU9tHs1MEGSUvmsiFgVEY8DnVRPXmZm1kStTj7Dc8njKWB4lTojgSdz80tT2fq8uU5ErAaeB7ZppC1Jx0rqkNTR1dXV036YmVkDmn7DgaSbgXdXWfTt/ExEhKRodjz1ioiLgYsB2traShOXmdnGoOnJJyIOqLVM0tOSRkTEckkjgGeqVFsG7J+bHwXc1sNmlwHbA0slDQa2AlbkyvNtLetpH8zMrH+1+rJbO9B999o04PoqdeYAB0oakm40ODCV1dvuwcAtERGpfGq6G24MMBa4p4/7YGZmDWp18jkL+LikxcABaR5JbZIuAYiIlcAZwPz0OD2VIel7kpYCb5e0VNJpqd2fAttI6gS+RrqLLiIeBq4CHgF+DZwQEWsK2VMzM3uTshMCW5+2trbo6OhodRhmZhsUSfdGRFvVZU4+PZPUBTzRhya2Bf67n8LpT46rMY6rMY6rMRtjXO+JiGHVFjj5FEBSR63s30qOqzGOqzGOqzEDLa5Wf+djZmYDkJOPmZkVzsmnGBe3OoAaHFdjHFdjHFdjBlRc/s7HzMwK5zMfMzMrnJOPmZkVzsmnn0ialAao65RUbVyiD0u6T9JqSQeXKK6vSXpE0oOS5kl6T0niOl7SQkn3S7o9jcXU8rhy9f5eUkgq5NbYOp6voyR1pefrfkmfL0Ncqc6h6TX2sKSfFRFXPbFJOjf3fP1J0nMliWsHSbdKWpD+Lw8qSVzvSe8RD0q6TdKoPm0wIvzo4wMYBDwK7AhsCjwAjKuoMxp4P3AZcHCJ4voo8PY0/QVgdkni2jI3/Ung12WIK9XbAvgdcBfQVoa4gKOA84t4XTUY11hgATAkzb+rLLFV1P8SMKMMcZF9wf+FND0OWFKSuH4OTEvTHwMu78s2febTP/YCOiPisYh4DZhFNnDdmyJiSUQ8CLxRsrhujYj/SbN3kfX0XYa4XsjNvgMo4s6YHuNKziAbHffVAmJqJK6i1RPXMcAFkQ0ESURU67m+VbHlHQ5cWZK4AtgyTW8F/KUkcY0DbknTt1ZZ3hAnn/7RmwHvitBoXEcDNzU1okxdcUk6QdKjwPeAL5chLkl7ANtHxK8KiKfuuJK/T5dErpa0fZXlrYhrZ2BnSXdIukvSpALiqjc2ILucBIxh7Rtrq+M6DTgidZp8I9lZWRniegD4dJr+FLCFpG16u0EnHwNA0hFAG3BOq2PpFhEXRMR7gX8GvtPqeCRtAvwQ+HqrY6nil8DoiHg/2VDzM3uoX5TBZJfe9ic7u/iJpK1bGVAVU4Grozw93B8OXBoRo4CDgMvTa6/VvgF8RNIC4CNkY6H1+jkrww5tDMo6SF1dcUk6gGxk2U9GxKqyxJUzC5jSzICSnuLaAtgNuE3SEmAfoL2Amw56fL4iYkXu2F0C7NnkmOqKi+wTdHtEvB4RjwN/IktGZYit21SKueQG9cV1NNnQL0TEncBfkXXu2dK4IuIvEfHpiNidNBJ1RDzX6y02+4usgfAg+3T3GNmpe/eXdbvWqHspxd1w0GNcwO5kXzSOLdPzlY8H+N9ARxniqqh/G8XccFDP8zUiN/0p4K6SxDUJmJmmtyW7tLNNGWJL9f4aWEL6wX0Z4iK79H1Umv4bsu98mhpfnXFtC2ySpqeTja3W+20W8YQPhAfZ6fGf0hv5t1PZ6WRnEwAfJPsU+DLZkN4PlySum4GngfvTo70kcf078HCK6db1JYEi46qoW0jyqfP5OjM9Xw+k5+uvSxKXyC5VPgIsBKYWEVe9x5Ls+5WzioqpzudsHHBHOpb3AweWJK6DgcWpziXAZn3ZnrvXMTOzwvk7HzMzK5yTj5mZFc7Jx8zMCufkY2ZmhXPyMTOzwjn5WN0kTcn3Li3p9PQD1WZu88rUZcw/9VDvpT5sY4mkt/yIT9Jpkr7R23br2O6NkraWNFrSQ71s45Pr6327jvV7vY+pJ+3zayzr9fGoY7tv7rOkS4vsJb63+nKMe2i3T8e/lQa3OgDrX5IGRRO6CZE0mKyXgRvIfrNBRJzS39up2Oa7gQ9GxE7N3E6rRMRBAH3pbiYi2oH2/oppQzAQ97mWDfm58JnPBiJ9cvqjpCsk/SF1Hvn2tGyJpLMl3QccIunwNBbOQ5LOzrXxUhrD5OE0LsewVD4+dfr4oKRrJQ1J5bdJ+jdJHWT9q30SOCeNf/Le/KdOSRPS+CMLJc2QtFkutu8qG8tooaS/rrJvfyXpP9PyBZI+mhb9BhiZtvehinXGSLozrfMvuXJJOift+0JJh6Xy/SXdkKt3vqSjck1+K9W/R9Jbkl3a319LulfS72vsx0e0dnyYBZK2SNv9naRfKRsr5Ufd/XRVO+OStGNa94N1bvPNs490PM6T9P8kPVbrjEDSt5WNX3M7sEuufH2vg7Y0va2yroW6bZ+WL5Z0ao3tfVPS/NTud6ssH5Ri7z5m/5Tb7r+n5/MhSXtV7nNFO2ekdgb1tM1Uv9rxeqey/43u1+vkVLf7/+/S9NxdIekAZR2mLs7Fdpqky9Nrc7GkY2rs7zm5+I6rUqfe7eWP//WSjkzTx0m6otp+l0aRv+z1o0+/Ph5N1tX636b5GcA30vQS4Ftpejvgz8AwsjPbW4ApaVkA/5imTyGN/wI8CHwkTZ8O/Fuavg24MBfDpeS6BuqeJ+t76klg51R+GfDVXGxfStNfBC6psm9fJ42lQtbdyZ9Tm6OBh2o8H+3AkWn6BOClNP33ZB1rDgKGp7ZGkHVseUNu/fNZ24XJEtb+ovvI7npkv37vfo7nkbr8AfYGbqkS0y9zx+ed6fnfn2zohR1TTHO7n8O03W2795MsESwAPtDANo/KHcdLycZc2YTsV/KdVervSdbTwNvJuu3vzO3j+l4HbWl6W9L4Mmnby4FtgM3TPnTX6z4eB5KNT6MU1w3Ah6vENDc3v3Vuuz9J0x/ufi1U2eeDyTrE/VHaTo/bXM/xGkwaSyrta2dqZzSwGnhfavNesv9BkQ0tcF3uNfNAej66uxPajtxrGTgW+E6a3gzoAMZU+X+vZ3v552J4ivdDZL0QDG31+9b6Hj7z2bA8GRF3pOn/AvbLLZud/n4QuC0iuiJiNXAF2T8uZGMJzc6vL2krsn/236bymbn6+XbXZxfg8Yj4U402fpH+3kv2T1VpvxQPEfFH4AmyrvjX529Z2xnk5RVtXRkRayLiaeC3ZM9JT67M/d03v0DSO4H/Bfxc0v3Aj8kSWqU7gB9K+jLZc7o6ld8T2Tgpa1L7+1VZdxhwPdmHgwca2Gal6yLijYh4hOzNqNKHgGsj4n8iGzOpPe1jT6+DWuZG1qnpK2THuXLfDkyPBcB9ZB8uKjsWfQzYUdJ/KBtyIT+W05UAEfE7YEtVv0T5f4GtIuL4yN6F69kmVD9eAv5V0oNkXU+NZO3z+HhELIyIN8i6MpqXtreQdV/X10fEKxHx32RdHe1V5Tk5Mh3Xu8mSd7X46t0e6Tl6muxD5a3A1yNiZZU2S8Pf+WxYKvtCys+/3A/tVdObdit197a8hv59zTXSN9Rq1r3M/Ffraauy3U2A5yJi/HqDiThL0q/I+si6Q9LEGu1Vi/t5srO0/ci+U6trm1XkeyVXg+vWkn/u1ve8VZsXcGZE/LhW4xHxrKQPABOB44FDgc/V2T7AfGBPSUPTG27VbUo6gWxwO4CDahyvfcg+COwZEa+nS4zd+5x/bt/Izb/Buq/rep6TL0XEnCr7klfv9vLeR9Z35HY9tN1yPvPZsOwgqftT+T8At1epcw/ZmBvbShpENjZI96fZTcguUby5fkQ8Dzyrtd+pfCZXv9KLZMMKVFoEjNba70rW10Y1vwf+EUDSzsAOqc31uYOsK3y61821dVi6rj6M7NP7PWRnU+MkbZY+PU+oaO+w3N878wvSGcLjkg5JMSq9Wa5D0nvTJ9Wzyd4Qu7+j2UvZd1SbpParHbfXyHqjPlLSP9S7zV74HTBF0uaStiDrMZweXgdLWDtEQ+X3SB+XNFTS5mQ3pNxRsXwO8Ll0JoekkZLela+g7HuvTSLiGrJxm/bILe7+zm4/4PkUZ6VfA2cBv0r7VHWbkY0PNT49/lLjeG0FPJMSz0eB91R/GtdrsrLvMbchu+w6v8pz8gVJb0vx7SzpHb3YzjrS90CfIOup/huSxvS1zWbymc+GZRFwgqQZZJ+OL6qsEBHLld16eSvZJ6xfRcT1afHLZG+E3wGeYe0b7jTgR8puYHgM+GyN7c8iGwzsy+TehCLiVUmfJbtENJjsn+1HDezXhcBFkhaSfco+KiJWSev94P4V4GeS/pnsclW3a8kumz1A9onzWxHxFICkq8i+l3ic7JJM3pB0qWUVWcKu9I8pxu8AbyN7Lh6oqPPV9IbVfZnkphTLfLLvmHYiOy7XVtuhiHhZ0t8Bc5XdqlzPNhsSEfdJmp3aeYZ13xhrvQ6+D1wl6VigcgTXe4BryMZ/+a+I6KjY3m8k/Q1wZzqeLwFHpG13Gwn8p9YOmHZybtmrygYvextrz4aq7dfPU+JpJzuT+VkP24Tqx2sL4JfptdgB/LHWNtfjQbLjvC1wRkp0o3PLLyG7bHafsgC76ON4Vcpu8PkJ8Nm0va8DMyR9LF2qKx33ar2BSC/eGyJitz608VJEvLP/orKeSNqf7Av9v2txKBscSbeRPXcdPdUtC0mnkd1s8f1Wx1J2vuxmZmaF85mPmZkVzmc+ZmZWOCcfMzMrnJOPmZkVzsnHzMwK5+RjZmaF+/+kFeTddtHD6gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(prop, alpha, yerr=alpha_err, fmt=\"o\")\n",
    "plt.xlabel(\"proportion of double spike in double spike-sample mix\")\n",
    "plt.ylabel(r\"Recovered $\\alpha$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compare the errors in $\\alpha$ obtained from the inversion with those expected from the linear error propagation. In this case we have an excellent match as the Monte-Carlo simulations are based on exactly the same error model as the linear error propagation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEaCAYAAABXZ4NKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+AElEQVR4nO3deXwddb3/8dc7e5qtbbrQpksKbYGWpUBZFESU3SsUroBFEFAUvai4XpWfV0Su3ituqBcUEVBEVkGhLMpSCghCaQu10FJK6EL3pmmbrc3++f0xk/aQnjQn25lzks/z8TiPzJnlO5+ZJOdz5jvf+X5lZjjnnHPJlhF1AM455wYnT0DOOeci4QnIOedcJDwBOeeci4QnIOecc5HwBOSccy4SnoCcc85FwhOQ6zFJ5ZIqJT0bvkbGLLtQUmVUMUg6SdJcSfMknZvsGCS9L+b9Ckk3JHP/4fwMSX+Q9A9JL0g6qL9icK4nsqIOwKW958zsvNgZkjKB84G1UcQgKR/4OnCmmTVFEQNQCZwUxvMH4KEk7x9gBpBrZh+Q9AHga8AV/RyHcwnzKyDXW8eH37D/R5LCeRcCfwbaIorhfcAu4BFJf5W0XwQxACApBzgG+EcE+18XhCABw4CtfblDSUslndSXZfZWKsbkOucJyPXGRmAycCIwCvj38OrnAuC+jitLWi1pl6S6mNfYvo4BGB3OOwv4HXBtBDG0OwWYa2ZtEex/K9AMLAf+D/h1X8ZgZtPN7Nlext6nUjEm1zlPQK7HzKzRzOot6FDwL8DhwMXA/e0fuHGcZWaFMa8N/RDDDuDFsPptLjA9ghjanU9wNRjF/k8DWszsQOBjwM/6K4ZkkeS3DQYQT0CuxyQVxbz9AFABTAMukfR3YIqkXyVY1lhJD4Y301dJuqoXMSwADg6rnmYAKyOIAUnZwNHAC1HsHxBQFU5vBUr6MobwSuqUDu+/IWmJpGpJ90nKk/QtSQ902PaX7X8bXe0zLPdbkpYA9ZKywvfrJdVKekvSyR1jknRw2ChjR1g1d3accveKt6tz5PqQmfnLXz16AWcCiwjub/wRyOqwfGGH96uBU+KUkxGWcw2QA+xPkDRO72kMwBeA54HngAMiiuFM4FcRnoMsgqrQ54CXgff3ZQwdywnfvwKMBYYDbwKfByYCO4GicL1MgmrD4xLZZ1juYmA8kA8cSNDAZWy4vLz9d9weE5BNkIj/X1juh4Fa4MCu4o36/2owvSIPwF+D5xX+w9cRVJHtAB4K5x8LvNth3auB3w+0GKLef1/G0EkCujjm/Y+Bm8PpF4BLwulTgXcS3WdY7qdj3k8GtrQnmngxEVwJbgIyYpbdA1ybSLz+Ss7L61Ndsp1jZk93mDcRGCtpR8y8TPqv5VjUMUS9//6MYVPM9E6CqwuAuwlaR/4R+ET4vjv73N2k38wqJH2FoHHJdElPAF+z997HGgustffei1wDlCUYr0sCT0AuFawFVpnZlEEcQ9T77+8Y/gz8TNI44FyCpvLd2ed7Rs40s7uBuyUVA78Frgc+GbPKBmC8pIyYJDQBWNG7w3B9yRshuFTwClAb3ljOl5Qp6RBJR0PwIKeChzkHcgxR77/LGHrDzCqBZ4HfEyScN3u6T0kHSvqwpFyggeCZr46tLucTXNF8U1K2gmeDzgLu7e2xuL7jCchFzsxagY8StFhbRdBi61b2tNoaD7w4kGOIev8JxtBbdxPcn2mvfuvpPnOBH4XrbiJ49unq2BUsaIJ/FkEDja0Ez0BdYmbL++ZQXF9QePPNuZSkoCeBfwGHmVnzYIwh6v071188ATnnnIuEV8E555yLhCcg55xzkfAE5JxzLhKegJxzzkXCH0QNjRgxwsrLyztdXrGljswMMWlEQfKCcs65frS5poEttY0cWtbzlvaLFi3aamYju15zb56AQuXl5SxcuLDT5RffOp+dTS385crjkxiVc871n+8+9AaPLtnAwmtO63EZktb0dFuvgktQQW4m9Y2tUYfhnHN9ZtvOJoYNyYls/56AElSYm01dY0vUYTjnXJ/ZsbOJYQWegFJeYW6mJyDn3ICyrb7Zr4DSQUFuFvWNLXjPEc65gWLHziaGDcmObP9JS0CSzgiHzq2Q9O04y3PDIXErJM2XVB7OL5U0T1KdpBtj1i+StDjmtVXSL8Jll4XD+7Yv+0xv4y/My6KlzWhs6djprnPOpR8zo6q+ieGF0V0BJaUVnKRM4CaCkRDXAQskzTGzZTGrXQ5sN7PJkmYTjO/xcYLu1r8LHBK+ADCzWoIedNv3sQj4S0x595nZF/vqGApzg1NV19hCXnZmXxXrnHORqGlooamljZGFuZHFkKwroGOACjNbGXaTfi8wq8M6s4A7wukHgJMlyczqzewFgkQUl6SpBF2y99fokRTkBAmo3u8DOecGgK11jQCMGAQJqIyYIXUJroI6Do27ex0zawGqgdIEy59NcMUTe4PmY5KWSHpA0vh4G0m6QtJCSQsrKyv3uYPCvCAB1TZ4AnLOpb+ttYMnAfW32cA9Me8fAcrN7DDgKfZcWb2Hmd1iZjPNbObIkft+kLe9Cs6vgJxzA0Fl+xVQ0cBvBbeeYETHduPCeXHXkZRFMCJiVVcFSzocyDKzRe3zzKzKzBrDt7cCR/U89EBBewJq8gTknEt/7VdAg+Ee0AJgiqRJ4eiOs4E5HdaZA1waTp8HPGOJtXm+kPde/SBpTMzbs4E36aX2KyCvgnPODQRb65rIzFCkzwElpRWcmbVI+iLwBJAJ3G5mSyVdByw0sznAbcCdkiqAbQRJCgBJq4FiIEfSOcBpMS3oLgA+0mGXV0k6G2gJy7qst8ewpwrOu+NxzqW/rXWNDC/IISNDkcWQtM5Izexx4PEO866JmW4Azu9k2/J9lLt/nHlXA1f3NNZ4CnKDptd+D8g5NxBU1jZG2gABBk4jhH7X3gy71hOQc24A2FrXyIgIH0IFT0AJy8gQRblZ1DY0Rx2Kc8712ta6pkgbIIAnoG4pzs+mZpdfATnn0puZUVnXyIgiT0Bpozg/mxq/AnLOpblU6IYHPAF1S3FeFtW7PAE559Lb1hR4CBU8AXVLUAXnCcg5l95SoRse8ATULSWegJxzA8DmMAGNKsqLNA5PQN1QnJdNjfeE4JxLc5urg8EF9iv2BJQ2SvKzqWtsoaXVB6VzzqWvTTUN5GVnUJyftL4I4vIE1A3tvyzvD845l8421TQwpiQfKbpueMATULcU5wVjp3tTbOdcOttc3cDo4mgbIIAnoG4pyQ8SkDfFds6ls001DZHf/wFPQN1SHCYg7w3BOZeu2tqMLTWNjC7xBJRW2u8BeRWccy5dbdvZRFNrm18BpRuvgnPOpbtNKdIEGzwBdcvuRgiegJxzaWpzTZCAvAouzQzJySQrQ34F5JxLW5tq/AooLUli6JActu9sijoU55zrkc3VDUgwMuKhGMATULcNL8hmW70nIOdcetpU08CIwlyyM6P/+I8+gjQzvCCH7fVeBeecS0/rd+yibGh+1GEASUxAks6Q9JakCknfjrM8V9J94fL5ksrD+aWS5kmqk3Rjh22eDctcHL5G7ausvjC8IIdtXgXnnEtT67bvYtywQZSAJGUCNwFnAtOACyVN67Da5cB2M5sM3ABcH85vAL4LfKOT4i8ysxnha0sXZfXasCE5XgXnnEtLbW3Gxh0NlA2mBAQcA1SY2UozawLuBWZ1WGcWcEc4/QBwsiSZWb2ZvUCQiBIVt6yeh7/H8IIcduxsorXN+qI455xLmsq6Rppa2xg3yKrgyoC1Me/XhfPirmNmLUA1UJpA2b8Pq9++G5NkelpWl4YX5NBm/iyQcy79rNu+C2DQXQH1l4vM7FDgA+Hrk93ZWNIVkhZKWlhZWZnQNsMLgjHUq7wazjmXZtZt3wnAuGFDIo4kkKwEtB4YH/N+XDgv7jqSsoASoGpfhZrZ+vBnLXA3QVVfwmWZ2S1mNtPMZo4cOTKhAxk2JEhA/iyQcy7drN8RXgENsiq4BcAUSZMk5QCzgTkd1pkDXBpOnwc8Y2ad3miRlCVpRDidDXwUeKMnZXVH+xWQN0RwzqWb9dt3MXRINgW50Y6E2i4pUZhZi6QvAk8AmcDtZrZU0nXAQjObA9wG3CmpAthGkKQAkLQaKAZyJJ0DnAasAZ4Ik08m8DTwu3CTTsvqrfYEtN0TkHMuzaTSM0CQpAQEYGaPA493mHdNzHQDcH4n25Z3UuxRnazfaVm91V4F5/eAnHPpZt32XRwwsiDqMHZL90YISZefk0l+dqZfATnn0oqZsX77LsqGpkYDBPAE1CPeG4JzLt1U1jayq7mV8hGegNJa0B+cJyDnXPpYtbUegImlXgWX1oYVeHc8zrn0sqYqeAZokieg9DZ8SDZb6zwBOefSx+qqerIyxNih0Q9E184TUA+MKs6jsq6RPnq0yDnn+t2aqp2MHz6ErBQYB6hd6kSSRkYV5dLU0kbNrpaoQ3HOuYSs2lpPeWnqNEAAT0A90j6U7Zba7nTQ7Zxz0TAz1lTVp1QDBPAE1COji4M61M01jRFH4pxzXdta10R9U6tfAQ0Eo/wKyDmXRtZUBU2wy0f4FVDaGxVeAW2p9Ssg51zqWx02wS73Krj0V5ibxZCcTLZ4FZxzLg2s3lpPZoZSZiC6dp6Aemh0cR6bvQrOOZcGVm2tZ/ywfLJTqAk2eALqsZFFuVT6FZBzLg2s2FzL5FFFUYexF09APTSqKNcbITjnUl5zaxurttYzdXRh1KHsxRNQD40uzmNLrfeG4JxLbWuq6mlpM6Z4Aho4RhXlsrOplbpG7w3BOZe6VmyuA2CKV8ENHKOK258F8vtAzrnU9fbmOiQ4YKRfAQ0Y+xUHzRk3Vft9IOdc6lqxpZbxw4aQn5MZdSh78QTUQ2VDgwS0fvuuiCNxzrnOVWyuS8kGCOAJqMf2K8lDgnU7PAE551JTS2sbK7fWpWQTbEhiApJ0hqS3JFVI+nac5bmS7guXz5dUHs4vlTRPUp2kG2PWHyLpMUnLJS2V9KOYZZdJqpS0OHx9pq+PJycrg9FFeX4F5JxLWaurdtLcakwZNYivgCRlAjcBZwLTgAslTeuw2uXAdjObDNwAXB/ObwC+C3wjTtE/NbODgCOA4yWdGbPsPjObEb5u7cPD2a1sWD4b/ArIOZeiVmyuBWDq6MF9BXQMUGFmK82sCbgXmNVhnVnAHeH0A8DJkmRm9Wb2AkEi2s3MdprZvHC6CXgVGNefB9HR2KH5rPcE5JxLUUs3VJOZoZR8BgiSl4DKgLUx79eF8+KuY2YtQDVQmkjhkoYCZwFzY2Z/TNISSQ9IGt/DuPepbGg+G6t30dbmD6M651LP0g01TBlVSF526rWAgwHQCEFSFnAP8CszWxnOfgQoN7PDgKfYc2XVcdsrJC2UtLCysrLb+y4blk9zq/mzQM65lLR0Qw3TxhZHHUankpWA1gOxVyHjwnlx1wmTSglQlUDZtwBvm9kv2meYWZWZtWeFW4Gj4m1oZreY2Uwzmzly5MhEjuM9xrU3xd6xs9vbOudcf9pS20BlbSPTx5ZEHUqnkpWAFgBTJE2SlAPMBuZ0WGcOcGk4fR7wjHXR0ZqkHxAkqq90mD8m5u3ZwJs9D71z7WNrrN/hD6M651LLsg01AExP4SugrO5uIKkAaDCz1kS3MbMWSV8EngAygdvNbKmk64CFZjYHuA24U1IFsI0gSbXvczVQDORIOgc4DagBvgMsB16VBHBj2OLtKklnAy1hWZd19zgTMdYfRnXOpailYQJK5Sq4LhOQpAyCZHARcDTQCORK2go8BvzWzCq6KsfMHgce7zDvmpjpBuD8TrYt7yy8Tta/Gri6q5h6qzA3i2FDsnl3m1fBOedSy7INNYwfnk9xXnbUoXQqkSq4ecABBB/o+5nZeDMbBZwAvAxcL+nifowxpU0sLWBNVX3UYTjn3Hss3VDN9DGpe/8HEquCO8XMmjvONLNtwIPAg5JSN8X2s0kjCpi/MpG2Es45lxzVO5tZXbWT82f2yxMofabLK6DY5CNppKS9movFS1CDRXlpARuqG2hoTviWmHPO9avF63YAcMT4oZHG0ZUuE5AC14b3fN4CVoT9rF3T1baDQfmIIQCsqfL7QM651PDau9uR4LB0T0DAV4HjgaPNbLiZDQOOJeh77av9Gl0aKC8tAGDVVr8P5JxLDa+9u4MDRxdRmNvths5JlUgC+iRwoZmtap8R9jhwMXBJfwWWLspHBAlotTdEcM6lADNj8dodzEjxqx9ILAFlm9nWjjPNrBIYtI0P2pXkZzO8IMdbwjnnUsKqrfVU72rmiAlDow6lS4kkoKYeLhs0ykuHeBWccy4lvPbuDgCOmDAs2kASkEgF4eGSamLexz78mdfH8aSl8hEFvPD2XheJzjmXdK++u53C3CwOGJmaQzDESqQZdqaZFce8imJeg74KDoLBnrbUNlK9c9C2RnfOpYhXVm3jqInDyMyI21FMSkmkGfbRkvaLeX+JpIcl/VLS8P4NLz1MDQd7WrGlNuJInHOD2da6Rt7eUsex+6fHR3Mi94B+S3ivR9KJwI+APxJ0BnpL/4WWPtqHu20f/tY556LwyqptABw7KaGxPCOXyD2gzLDbHYCPA7eYWXsXPIv7LbI0UjY0n4KcTFZs8gTknIvO/JVV5Gdncti41O4Drl0iV0CZ4QBxACcDz8QsS+2nnJJEEpNHF7Fic13UoTjnBrH5q7Yxs3wY2ZnpMdh1IlHeAzwn6WFgF/APAEmTgep+jC2tHDi6kLf9HpBzLiLb6ptYvqmWYyelx/0fSKwV3A+BrwN/AE6IGaU0A/hS/4WWXqaOLmJrXRNVdY1dr+ycc33slVVBr/zH7p8e938gwSo0M3s5zrwVfR9O+jpwv6AhwvJNtRw/OTfiaJxzg81zK7ZSmJuVFl3wtOtVRaGkT/VVIOlu+tjgpt/SDV4r6ZxLLjPj+RWVvP+A0rS5/wO9TEDA9/skigFgeEEOZUPzeX19TdcrO+dcH3qnsp71O3Zx4tS9hmtLaV1WwUla0tkiYHTfhpPepo8tZul6vwJyziXX8ysqAfjgQEtABEnmdGB7h/kC/tnnEaWxQ8pKeHLZZmobminK816KnHPJ8fzblew/ooDxw4dEHUq3JFIF9yhQaGZrOrxWA88muiNJZ0h6S1KFpG/HWZ4r6b5w+XxJ5eH8UknzJNVJurHDNkdJej3c5leSFM4fLukpSW+HP5PSLeyhZcF9oGUbvBrOOZccu5paeXllVdpVv0FiCegzZvZCvAVm9gkIhu3eVwGSMoGbgDOBacCFkqZ1WO1yYLuZTQZuAK4P5zcA3wW+Eafo3wCfBaaErzPC+d8G5prZFGBu+L7fTS8rBuANT0DOuST5x9uVNDS3ccrB6XdHJJEENE/SlyRNiJ0pKUfShyXdAVzaRRnHABVmttLMmoB7gVkd1pkF3BFOPwCcLElmVh8mwIYO+x8DFJvZy+GzSX8EzolT1h0x8/vVqKI89ivOY/HaHcnYnXNuMFpyP9xwCFw7FG44hC0v/omS/Oy06YA0ViL3gM4APg3cI2kSsINgHKBM4EngF2b2WhdllAFrY96vA47tbB0za5FUDZQCnQ20UxaWE1tmWTg92sw2htOb6KSxhKQrgCsAJkyYEG+VbjuqfBiLVm/rekXnnOuuJffDI1dB867gffVaPlZ9PYz9JtmZp0UbWw8k0hNCg5n92syOByYS9Ad3pJlNNLPPJpB8IhVeHVkny24xs5lmNnPkyL6pP505cRgbqhvYsGNXn5TnnHO7zb1uT/IJ5dPEx3bcHlFAvdOt54DMrNnMNprZjm7uZz0wPub9uHBe3HXCzk9LgKouyhzXSZmbwyq69qq6Ld2Mt8dmTgwugxeu6dho0Dnneql6XdzZeTs3xp2f6pL1yOwCYIqkSZJygNnAnA7rzGHPvaTzgGdi+p3bS1jFViPpuLARxCXAw3HKujRmfr87eEwRQ3IyvRrOOdf3SsbFna1O5qe6pCQgM2sBvgg8AbwJ3G9mSyVdJ+nscLXbgFJJFcDXiGm5Jmk18HPgMknrYlrQXQncClQA7wB/C+f/CDhV0tvAKeH7pMjKzGDG+KF+BeSc63snXwPZ+e+Z1ZKZF8xPQwmP5yPpKDNb1NMdmdnjwOMd5l0TM90AnN/JtuWdzF8IHBJnfhXBvapIzJw4jBvnVVDX2EJhrg+Z5JzrI4ddEPycex1WvY4NVkrpv/2QrPb5aaY7V0CXSnpA0nHtMyT9vB9iSnszy4fTZrDIr4Kcc33tsAto+tISjsi4nx8f/AB5R86OOqIe604C2gLMBP4S9miwCtivf8JKb0eXDycnM4MX3q6MOhTn3AD03IpKduxs5pwZZV2vnMK6Uz90MXCgmTVKGgv8L5DSTbCjkp+TyVETh/GPtzt7hMk553ruodfWU1qQwwlTRkQdSq905wpoLTAJwMw2mNmlwOf6JaoB4IQpI1i+qZbKWh8h1TnXd6p3NvP0m5v56GFj0mrsn3i6E/2XgQcl3Snpa5J+BtT3U1xp78QpwYOtL1b4VZBzru88+Oo6GlvauODo8V2vnOISTkBmtgw4kqAft3yCLm469ufmQtPHFjNsSDbP+30g51wfMTPufuVdDh8/dPcozOmsW22EzawReCx8uX3IyBDvnzyCf7y9lbY2IyNjnx2GO+dclxau2U7Fljp+/LHDog6lT6R3BWKKO+XgUVTWNvKvdTuiDsU5NwDcPf9dinKz+OjhY6IOpU94AupHHz5wNJkZ4sllm6MOxTmX5rbWNfLY6xs598gyhuQMjAfcPQH1o5Ih2Ry3/3CeXLop6lCcc2nuzpfW0NTSxqXvL486lD7Tna54coGPAeWx25nZdX0f1sBx2rT9+N6cpVRsqWPyqMKow3HOpaGG5lbufHkNpxw8igNGDpzPke5cAT1M0OqthaD5dfvL7cOp04Kx8J7wqyDnXA89+Oo6ttU38dkP7B91KH2qOxWJ48zsjH6LZIAaOzSfGeOH8si/NvCFD02OOhznXJppbTNue2EVh40r4ZhJ6Tfs9r505wron5IO7bdIBrBzjyhj+aZa3txYE3Uozrk089jrG1lZWc/nTjyAYOizgaM7CegEYFHYEekSSa9LWtJfgQ0kHz1sDJkZ4qHFHQeBdc65zrW2Gb98egUHji7izEMGXt/P3amCO7PfohjgSgtz+eDUkTz82ga+dfpB/lCqcy4hj/xrA+9U1vPri44ckJ8b3emKZ028V38GN5Ccc0QZm2oaePEd7xvOOde1ltY2fjX3bQ7ar4gzpg+8qx9IIAFJeiH8WSupJuZVK8lvaiTotGmjGTYkmz+97DnbOde1v7y2npVb6/nKKVMG5NUPJJCAzOyE8GeRmRXHvIrMrLj/QxwY8rIzueDo8Tz95hY2Vu+KOhznXAqrb2zhJ0+8xREThnL6AL36Ae8JIakuOmYibWbc88raqENxzqWwm597h8raRr770WkDruVbLE9ASTShdAgnTR3JPa+8S3NrW9ThOOdS0Podu7jl+ZWcffhYjpwwLOpw+lVCCUiBXo1+JOmMsAl3haRvx1meK+m+cPl8SeUxy64O578l6fRw3oGSFse8aiR9JVx2raT1Mcs+0pvY+9In3zeRytpGHluyMepQnHMp6H8ffxOAb515UMSR9L+EEpCZGfB4T3ciKRO4iaAp9zTgQknTOqx2ObDdzCYDNwDXh9tOA2YD04EzgF9LyjSzt8xshpnNAI4CdgJ/jSnvhvblZtbj2PvaSVNHMXV0Ib959h3a2izqcJxzKWTe8i08umQjV540mbKh+VGH0++6UwX3qqSje7ifY4AKM1tpZk0Eo6p2HE11FnBHOP0AcLKCys9ZwL1m1mhmq4CKsLxYJwPvpEOz8IwMceVJk3lrcy1zl2+JOhznXIqob2zhvx56gymjCvmPkw6IOpyk6E4COhZ4SdI7PegJoQyIvfO+LpwXdx0zawGqgdIEt50N3NNh3hfDOG+XFLciVdIVkhZKWlhZmbyhsz962BjGD8/nxnkVBBeXzrnB7mdPrmD9jl386GOHkpM1OG7Pd+coTwcOAD4MnAV8NPwZKUk5wNnAn2Nm/4Yg1hnARuBn8bY1s1vMbKaZzRw5cmR/h7pbVmYGn//gAfxr7Q6ef9sfTHVusHtl1Tb+8M9VXHzcBI6aOLA6HN2XbvWEAAwlSDpnAUO7UeW1HohtxDAunBd3HUlZQAlQlcC2ZwKvmtnuYUfNbLOZtZpZG/A79q6yi9x5R41j/PB8rv/bcr8X5NwgVtPQzFfvW8z44UO4+syDow4nqRJOQJK+DNwFjApff5L0pQQ3XwBMkTQpvGKZDczpsM4c4NJw+jzgmbDxwxxgdthKbhIwBXglZrsL6VD9Jil2wPRzgTcSjDNpcrMy+cZpB7JsYw0P/8s7KXVusLr24aVsqmngho/PoCB3YAy1najuHO3lwLFmVg8g6XrgJeD/utrQzFokfRF4AsgEbjezpZKuAxaa2RzgNuBOSRXANoIkRbje/cAygsHwvmBmrWEMBcCpwOc67PLHkmYABqyOszwlnHXYWH73j5X89IkVnHnIGPKyM6MOyTmXRA+9tp6/vLaeL588ZcA/8xOPEr0JLul14Ggzawjf5wELzGxAjBE0c+ZMW7hwYdL3+8+KrXzi1vncfPg7nLHpFqheByXj4ORr4LALkh6Pcy453txYw7m/fpHDyoZy12ePJTszPRseSFpkZjN7sm13roB+D8yX1P6szTkEVy2uF94/eQTfK1/Kict/AmoKZlavhUeuCqY9CTk34FTvaubzf1pEcV42N150RNomn95KuCcEglZmnyKoHtsGfMrMftF/oQ0el9TfwZD25NOueRfMvS6agJxz/aa1zfjafYvZsGMXv7n4SEYV5UUdUmQSugIyM5P0eFjd9mo/xzToZNZ20gihel1yA3HO9bsfPvYmc5dv4b9nTR9UTa7jSVZPCG5fSsZ1b75zLi3d9sIqbn9xFZ8+fhKffF951OFELlk9Ibh9OfkayH5vv0+WnR/Md84NCH97fSM/eGwZZ0zfj+/82+B63qczCVXBhfeArgBSvq+1tNTe0GDudVj1Ota3lfLm5K9wqjdAcG5AmLd8C1fd+xpHjB/KL2bPIHOAjnDaXd25B3TTQGlynZIOuyB4mfG/d7/G3xZv5M6jtnL85BFRR+ac64XnV1TyuT8t4sD9ivj9p47x5/1i+D2gFCOJ6887jANGFvLFu19l3fadUYfknOuhf1Zs5bN/XMgBIwv50+XHUpKfHXVIKaW794Be9ntA/a8wN4vffvIoWlqNz925iPrGlqhDcs5105NLN3HZHxZQXlrAXZ85lqFDcqIOKeV0tzfs/Umx3rAHqv1HFvKrC4/gzY01XHnXqz6Et3OpYMn9cMMhcO3Q4OeS++Oudv+CtXz+T4uYNqaYe684juEFnnzi6TIBSfom7O4N+xgzW9P+IkX7WBsoPnTQKH547qE8t6KSbz/4uo8d5FyUltwf9FBSvRawPT2WxCQhM+OmeRV888ElnDBlJHd95liGefLpVCJXQLNjpq/usOyMPozFxXHhMRP48slTePDVdXz/kWWehJyLytzrgh5KYsX0WNLQ3MpX7lvMT554i3NmjOXWS2YOut6tuyuRs6NOpuO9d/3gK6dMob6xhVtfWEWbGd8/ezpBy3jnXNJ01jNJ9To2Vu/iij8u4o0N1fzn6Qdy5UkH+P9oAhJJQNbJdLz3rh9I4jv/djAZGeKW51fS2mZcN+sQf5bAuWQqGRdWv71XQ8EYzvq/F9nV1MLvPjmTU6aNjiC49JRIAjpcUg3B1U5+OE34fvD2opdkkrj6zIPIkLj5uXeoqmviF7Nn+DMFziXLydcE93xiquGaM/L41o5zKBmexd2fPZapo4siDDD9dJmAzMw/4VKEJL595kGMLMrlB48t4xO/e5lbLz3aW9g4lwwdeiypzBjJDxrOI//IC3nk7GkMyfH7Pd3lZywNXX7CJMaW5PHl+xZz9o0vcPPFR3FIWUnUYTk34LUecj531BzNT554iywTP/z4oZx9+Niow0pbg3MUpAHgzEPHcP/n3kdrm/Hvv/kn9y/Yu27aOdd3KrbUcv7N/+S6R5dx3P7DeeKrJ3ry6SW/AkpjM8YP5dEvncCX713MNx9cwvxV27j27GkU5Xl3H871lbrGFm6aV8Ft/1jFkNxMfn7B4Zx7RJm3cusDnoDSXGlhLnd8+hh++fQKbpxXwcsrq/jZBYdz3P6lUYfmXFprazP++tp6rv/7crbUNvLvR5Zx9ZkHM7IoN+rQBgz5g42BmTNn2sKFC6MOo1cWrdnO1+9fzJptO7n0feV8/bSpfjXkXA+89E4V1/99OYvX7uDwcSV87+zpHDlhWNRhpSRJi8xsZk+2Tdo9IElnSHpLUoWkb8dZnivpvnD5fEnlMcuuDue/Jen0mPmrw05RF0taGDN/uKSnJL0d/hwUfzlHTRzGY1d9gIuPncgdL63m5J89xyP/2uC9JziXoFff3c5Ft77Mhb97mY3Vu/jJeYfx1yuP9+TTT5JyBSQpE1gBnAqsAxYAF5rZsph1rgQOM7PPS5oNnGtmH5c0DbgHOAYYCzwNTDWzVkmrgZlmtrXD/n4MbDOzH4XJbpiZfWtfMQ6EK6BYi9fu4L8eep031tdw/ORS/t9HDmb6WG8p51w8C1dv49fPvsMzy7dQWpDDlR+azEXHTvDn7BLQmyugZN0DOgaoMLOVAJLuBWYBy2LWmQVcG04/ANwYjsQ6C7jXzBqBVZIqwvJe2sf+ZgEnhdN3AM8C+0xAA82M8UN5+AsncNf8Nfz8qRX8269eYNaMsXz91AOZUDok6vCci1xbm/H0m5v57fMrWbRmO0OHZPOfpx/IZe8v9z7ckiRZZ7kMiG0nvI5gfKG465hZi6RqoDSc/3KHbcvCaQOelGTAb83slnD+aDPbGE5vAuL2jSHpCoKhxpkwYUIPDiu1ZWaIS95XzqwZZfz2uXe4/cVVPP76Rj5+9Hg+d+IBjB/uicgNPrUNzTz02np+/8/VrKysZ/zwfL5/9nTOnznOHyZNsnQ/2yeY2XpJo4CnJC03s+djVwiHE49bzxgmrFsgqILr/3CjUZKfzTfPOIhL31/OL+e+zX0L1nLPK2s567AxfP6kAzhov+KoQ3Su372xvpq75q/h4cUb2NnUyqFlJfzfhUdw5iH7kZXpj0RGIVkJaD0wPub9uHBevHXWScoCSoCqfW1rZu0/t0j6K0HV3PPAZkljzGyjpDHAlr4/pPQzujiP/zn3UK768BRue2Eld81/l4cWb+DEqSP55HETObn5OTKeuS7o9bdkXND3VXv3I86loaq6Rh5dspEHX13HknXV5GVnMOvwMi46bgKHjRsadXiDXrIaIWQRNEI4mSB5LAA+YWZLY9b5AnBoTCOEfzezCyRNB+5mTyOEucAUgo5QM8ysVlIB8BRwnZn9XdJPgKqYRgjDzeyb+4pxoDVCSMSOnU388aU13DV/DcfWzeX6nNvIp3HPCtn5cNavPAm5tLKrqZUnl23iodfW8/zbW2ltMw7ar4jZR4/n3CPHUZLvjyb0pd40Qkjac0CSPgL8AsgEbjezH0q6DlhoZnMk5QF3AkcA24DZMY0WvgN8GmgBvmJmf5O0P/DXsPgs4G4z+2G4filwPzABWANcYGbb9hXfYExA7Zpb22j52XTyd27Ya1lb8TgyvrY0zlbOJdGS+4OB3zq5Oq/e2cy8t7bw5LJNPPtWJTubWhlbksfZM8o454ixXs3cj9IiAaW6wZyAgGCM+zjDO7WZ+NrB8zjr8LEcP3mEN0t1ydc+FHbsaKTZ+Ww/+ac8aifwxNLNvLyyipY2Y1RRLqdOG81Zh4/lmPLhZPiYWf0uHZphu1TXyWBb1TmjeGb5Fh5avIH87Ew+OHUkp04bzYcPGuVj3bvk6GQo7PrHv8d3m37F/iML+OyJ+3PatNEcPm6oJ5004gnIBeIMtkV2PsPO+gELpp3Cyyu38dSyTTy9bAt/X7qJzAxx5IShHD95BCdMHsHh44eS7S2JXB/a2dTC4nd38L7qdcRLKWUZVTz9tROZPMoHgUtXXgUXGvRVcNBlPTuAmfH6+mqeXLqZ59+u5PX11ZhBQU4mx+1fyvsnj+Do8mEcPKbYE5Lrli21DSxavZ0Fq7ezaM023thQQ2ub8ULuVYzT1r03KBkPX30j+YG69/B7QH3AE1DP7NjZxEvvVPFCxVZerNjK6qqdAORlZ3D4uKEcNXEYR00cxhEThvnIrekqgS8m3bWltoGl62t4fX01r6+v5o311WysbgAgNyuDGeOHcnT5cI4qH8axtXMZ8sRX97o69xaaqcETUB/wBNQ3NlbvYtGa7Sxas51X12xn6YYaWtqCv7GyoflMG1vM9LHFTB9bwvSxxYwpydt7XJV++MDrtlSIIRV00gAg0Q//huZW3qms4+3NdazYXMuKzbW8vr6azTVBc38JJo0o4NCyEg4tK+HIicM4ZGwJOVkdrp7995GyPAH1AU9A/WNXUytL1u3gtbU7WLqhhqUbqlm1tZ72P7thQ7KZMrqIyaMKmTyykOPqn+GgBd8hoyXCb7u9/NAdUG44JG7jlNjqLzOjqr6JNVX1rKnaycrKelZsruXtLXWsqaon/P5BVoYoH1HAIWOLOSRMONPLSij0ftfSmiegPuAJKHnqG1tYvqkmSEjra6iorKNiSx3Vu5p5IecqxmXsXd9flzeGV//9H4wfPoSyofl7f0PuSwl86CZFKnzr76R5viGunDyXNVU7WVNVT31T6+5lmRmivHQIU0cXMWV0EVNHFzJ1dBHlpQX9+3tzkfBm2C6tFORmcdTE4Rw1cfjueWbG1romRvysKu42Q3Zt4pLbXwGCapv9ivMYP2wIo4pzGV2cx+jw56iivN3zevzNunpd9+b3h45XYdVrg/fQp0moobmVqvomquoaqaprorK2kY3VDWys3sXG6gZ+rBGMtsq9tttgpby1uZby0gKOmTSciaVDKC8tYELpEMYPG+KJxiXEE5BLCZKCoY47eR6prbiM+849jrXbd7F2207Wbt/J+u27WLqhhrlvbmFXc+te2xTkZDK8MIdhQ9pf2QwNp4cX7JkuzMuiMDeTgtwsCnKzKCoZh+JeAY3rj0OPr5NnX5h7XdyWibuaW6nZ1UJNQzM1u5qpaWimtqElnA5+Vu9qZmtdE1X1QbKpqmt8z5VLrBGFuYwpyeOREZ/hkq0/J8f2dNHUlpXPmLP+l2cOP6mvj9oNMp6AXGrp5HmkrFO/x7H7l+41hgcEH8B1jS1srmlkS00Dm2sb2FLTyKaaBrbXN7FtZzPbdzaxcmsdO+qbqW1s2WcIZ2ecxfXZt5Kvpt3zGsjltoxPsPD3r5CTlUFOVia5WRnBdGbG7ukMCYngJ+x+KHLPfBCi1YzWNqOl1Whpa6Olbc/71rY2ru3k2Ze26nWc/vPn2NXcSkNzK7uaWmloaaO1bd9V6blZGRTnZ1NakMOIwlwmTBhCaUEupYU5jCjMiZkOrh73XMGcAEv2f09VYIY3AHB9xBOQSy3tH2zduPchiaK8bIryspk8qrDLXTS3trEjTEo7djZT19hMXWMr9Y0t1De2UNc4hb+vH8VJ637D0OYtVGWO5N6iT/Fszkk01TfR2NxGU2sbTS1tNLa00tgSTDe1ttHTW6rZmSIrI4OsDJGZKT5PKWPY+17YtsyRTBldSF52JvnhKy87k6K8LIrzs4OfedkU52dTnJcVnpes3nWhdNgFnnBcv/BGCCFvhOD6gplhFty2bwun28L/sdj3WRkZZGaIrAzF7zrGW+K5NOGNEJxLEQqr2gAy41aiJagHV4LOpRtPQM6lKq/6cgOct5V0zjkXCU9AzjnnIuEJyDnnXCQ8ATnnnIuEJyDnnHOR8ATknHMuEp6AnHPORSJpCUjSGZLeklQh6dtxludKui9cPl9Secyyq8P5b0k6PZw3XtI8ScskLZX05Zj1r5W0XtLi8PWRpBykc865hCXlQVRJmcBNwKnAOmCBpDlmtixmtcuB7WY2WdJs4Hrg45KmAbOB6cBY4GlJU4EW4Otm9qqkImCRpKdiyrzBzH6ajONzzjnXfcm6AjoGqDCzlWbWBNwLzOqwzizgjnD6AeBkBWM1zwLuNbNGM1sFVADHmNlGM3sVwMxqgTeBsiQci3POuT6QrARUBsQOsLKOvZPF7nXMrAWoBkoT2TasrjsCmB8z+4uSlki6XdKweEFJukLSQkkLKyv3HnTLOedc/0n7RgiSCoEHga+YWU04+zfAAcAMYCPws3jbmtktZjbTzGaOHDkyGeE655wLJSsBrQfGx7wfF86Lu46kLKAEqNrXtpKyCZLPXWb2l/YVzGyzmbWaWRvwO4IqQOeccykkWQloATBF0iRJOQSNCuZ0WGcOcGk4fR7wjAWDFc0BZoet5CYBU4BXwvtDtwFvmtnPYwuSNCbm7bnAG31+RM4553olKa3gzKxF0heBJ4BM4HYzWyrpOmChmc0hSCZ3SqoAthEkKcL17geWEbR8+4KZtUo6Afgk8LqkxeGu/p+ZPQ78WNIMgnHBVgOfS8ZxOuecS5yPiBryEVGdc677ejMiato3QnDOOZeePAE555yLhCcg55xzkfAE5JxzLhKegJxzzkXCE5BzzrlIeAJyzjkXCU9AzjnnIuEJyDnnXCQ8ATnnnIuEJyDnnHOR8ATknHMuEp6AnHPORcITkHPOuUh4AnLOORcJT0DOOeci4QnIOedcJDwBOeeci4QnIOecc5HwBOSccy4SnoCcc85FImkJSNIZkt6SVCHp23GW50q6L1w+X1J5zLKrw/lvSTq9qzIlTQrLqAjLzOn3A3TOOdctSUlAkjKBm4AzgWnAhZKmdVjtcmC7mU0GbgCuD7edBswGpgNnAL+WlNlFmdcDN4RlbQ/Lds45l0KSdQV0DFBhZivNrAm4F5jVYZ1ZwB3h9APAyZIUzr/XzBrNbBVQEZYXt8xwmw+HZRCWeU7/HZpzzrmeyErSfsqAtTHv1wHHdraOmbVIqgZKw/kvd9i2LJyOV2YpsMPMWuKs/x6SrgCuCN82SnqjG8c0kI0AtkYdRIrwc7GHn4s9/FzscWBPN0xWAkpJZnYLcAuApIVmNjPikFKCn4s9/Fzs4ediDz8Xe0ha2NNtk1UFtx4YH/N+XDgv7jqSsoASoGof23Y2vwoYGpbR2b6cc85FLFkJaAEwJWydlkPQqGBOh3XmAJeG0+cBz5iZhfNnh63kJgFTgFc6KzPcZl5YBmGZD/fjsTnnnOuBpFTBhfd0vgg8AWQCt5vZUknXAQvNbA5wG3CnpApgG0FCIVzvfmAZ0AJ8wcxaAeKVGe7yW8C9kn4AvBaW3ZVb+uhwBwI/F3v4udjDz8Uefi726PG5UHDB4JxzziWX94TgnHMuEp6AnHPORWLQJaDedAk00CRwLr4maZmkJZLmSpoYRZzJ0NW5iFnvY5JM0oBtgpvIuZB0Qfi3sVTS3cmOMVkS+B+ZIGmepNfC/5OPRBFnf5N0u6QtnT0rqcCvwvO0RNKRCRVsZoPmRdBY4R1gfyAH+BcwrcM6VwI3h9OzgfuijjvCc/EhYEg4/R+D+VyE6xUBzxM8GD0z6rgj/LuYQtC4Z1j4flTUcUd4Lm4B/iOcngasjjrufjoXJwJHAm90svwjwN8AAccB8xMpd7BdAfWmS6CBpstzYWbzzGxn+PZlgmeqBqJE/i4A/pugn8GGZAaXZImci88CN5nZdgAz25LkGJMlkXNhQHE4XQJsSGJ8SWNmzxO0Tu7MLOCPFniZ4FnMMV2VO9gSULwugTp20/OeLoGA9i6BBppEzkWsywm+4QxEXZ6LsEphvJk9lszAIpDI38VUYKqkFyW9LOmMpEWXXImci2uBiyWtAx4HvpSc0FJOdz9PgEHeFY9LjKSLgZnAB6OOJQqSMoCfA5dFHEqqyCKohjuJ4Kr4eUmHmtmOKIOKyIXAH8zsZ5LeR/As4yFm1hZ1YOlgsF0B9aZLoIEmkXOBpFOA7wBnm1ljkmJLtq7ORRFwCPCspNUEddxzBmhDhET+LtYR9DrSbEEP9SsIEtJAk8i5uBy4H8DMXgLyCDoqHWwS+jzpaLAloN50CTTQdHkuJB0B/JYg+QzUen7o4lyYWbWZjTCzcjMrJ7gfdraZ9bgTxhSWyP/IQwRXP0gaQVAltzKJMSZLIufiXeBkAEkHEySgyqRGmRrmAJeEreGOA6rNbGNXGw2qKjjrRZdAA02C5+InQCHw57AdxrtmdnZkQfeTBM/FoJDguXgCOE3SMqAV+E8zG3C1BAmei68Dv5P0VYIGCZcNxC+sku4h+NIxIrzf9T0gG8DMbia4//URgvHadgKfSqjcAXiunHPOpYHBVgXnnHMuRXgCcs45FwlPQM455yLhCcg551wkPAE555yLhCcgt5ukcyRNi3l/Xfggan/u856w99yvdrFeXS/2sTp8XqXj/GslfaOn5Saw38clDZVU3lkvwgmUcfa+eudOYPseH6OkyyTd2MmyHv8+Etjv7mOW9AdJ5/XXvvpKb37HXZTbq99/qhtUzwENBJIyLRySvI/LzQLOAR4lGP4cM7umr/fTYZ/7AUeb2eT+3E9UzOwjAJKG9qKMOez98OOANhiPuTMD/Vz4FVCKCL9BLZd0l6Q3JT0gaUi4bLWk6yW9Cpwv6UJJr0t6Q9L1MWXUSbohHKNlrqSR4fwZYaeRSyT9VdKwcP6zkn4haSHwLeBs4CeSFks6IPbbp6STFYx58rqCsUFyY2L7vqRXw2UHxTm2PEm/D5e/JulD4aIngbJwfx/osM0kSS+F2/wgZr4k/SQ89tclfTycf5KkR2PWu1HSZTFFfjNc/xVJeyW88Hj/LmmRpH90chwfDGNdHB5HUbjf5yU9pmDcmJsV9B0X98pL0v7htkcnuM/dVyHh7+NXkv4paWVnVwaSviNphaQXgANj5u/r72BmOD1CQXdD7caHy9+W9L1O9vefkhaE5X4/zvLMMPb239lXY/b7y/B8viHpmI7H3KGc/w7Lyexqn+H68X5fhQr+N9r/XmeF67b///0hPHd3STpFQYerb8fEdq2kO8O/zbclfbaT4/1JTHyfi7NOovuL/f0/LOmScPpzku6Kd9xpJepxJvy1ezyNcoInqY8P398OfCOcXg18M5weS9D9x0iCK9hngHPCZQZcFE5fA9wYTi8BPhhOXwf8Ipx+Fvh1TAx/AM7r+J6ge5G1wNRw/h+Br8TE9qVw+krg1jjH9nWCp8gBDgrjzwuPubPxReYAl4TTXwDqwumPAU8RPJk+OixrDMFT2o/GbH8jwVPp7TF+J5y+pH09gp6M28/xXGBKOH0sQRdMHWN6JOb3Uxie/5MIhmfYP4zpqfZzGO53RPtxEiSD14DDu7HPy2J+j38A/kzwxXEawVABHdc/CngdGEIwTEBFzDHu6+9gZjg9gnBMm3DfGwl6g88Pj6F9vfbfx2kEY+IojOtR4MQ4MT0V835ozH5/F06f2P63EOeYzyPolePmcD9d7nMfv68soDjmWCvCcsqBFuDQsMxFBP+DIhhq4KGYv5l/hedjBMH/xVhi/paBK4D/CqdzgYXApDj/74nsL/ZcjA7j/QBB/3vDo/7c6u3Lr4BSy1ozezGc/hNwQsyy+8KfRwPPmlmlBcNF3EXwzwvQFrPen4ATJJUQ/MM/F86/I2b92HL35UBglZmt6KSMv4Q/FxH8Y3V0QhgPZrYcWEPQf9i+HA/cE07f2aGse8ys1cw2A88RnJOu3BPz832xCyQVAu8n6HJoMUH/d/HGMnkR+LmkqwjOaUs4/xULxoxpDcs/Ic62I4GHCb4g/Ksb++zoITNrM7NlBB9IHX0A+KuZ7TSzGsLqmwT+DjrzlJlVmdkugt9zx2M7LXy9BrxK8AWjY8ekK4H9Jf2fgqEbamKW3QO7x5spVvzqyu8CJWb2eQs+iRPZJ8T/fQn4H0lLgKcJhgxoP4+rzOx1C3qyXgrMDff3Ou/9u37YzHaZ2VZgHsG4QR3PySXh73U+QQKPF1+i+yM8R5sJvljOA75uZvsanyct+D2g1NKxX6TY9/V9UF48PSm3o/Zeslvp27+p7vQT1cJ7q5Tz9lFWx3IzgB1mNmOfwZj9SNJjBH1evSjp9E7Kixd3NcHV2gkE99gS2mccsT2S99VAibHnbl/nLd57Af9rZr/trHAz2y7pcOB04PPABcCnEywfgk5Bj5I0PPzQjbtPSV8gGCwP4COd/L6OI/gycJSZNYfVje3HHHtu22Let/Hev+tEzsmXzOyJOMcSK9H9xTqUoHf+sV2UnRb8Cii1TFAwpgjAJ4AX4qzzCvDBsK4+k2A8kvZvtRkE1RW7tzezamC79txj+WTM+h3VEgw90NFbQLn23DvZVxnx/AO4CEDSVGBCWOa+vMiejmAv6lDWx8N69pEE3+JfIbiqmiYpN/wWfXKH8j4e8/Ol2AXhlcIqSeeHMSr8wHwPSQeE31ivJ/hQbL9nc4yCe1YZYfnxfm9NwLkE34w/keg+e+B54BxJ+ZKKgLPCY9zX38Fqgmoy2PP30+5UScMl5RM0Unmxw/IngE+HV3RIKpM0KnYFBffBMszsQeC/CIZ2btd+D+8Egh6Uq+Mc09+BHwGPhccUd59mdpOZzQhfGzr5fZUAW8Lk8yFgYvzTuE+zFNzXLCWogl0Q55z8h6TsML6pkgp6sJ/3CO8LnQkcAXxD0qTelhk1vwJKLW8BX5B0O8G35N90XMHMNipoljmP4JvWY2b2cLi4nuDD8L+ALez50L0UuFlBo4aVdN5T7b0EPfteRcwHkZk1SPoUQXVRFsE/3M3dOK5fA7+R9DrBt+3LzKxR+x7p/MvA3ZK+RVB11e6vBFVo/yL45vlNM9sEIOl+gvsUqwiqZ2INC6tdGgmSdkcXhTH+F0Evv/eG+4j1lfBDq73K5G9hLAsI7jlNJvi9/DXeAZlZvaSPAk8paMacyD67xcxelXRfWM4W3vvh2NnfwU+B+yVdAXQc8fUV4EGC8V3+ZB2GoDCzJxUMQ/BS+PusAy4O992uDPh9mKABro5Z1iDpNYLj/zSdMLM/h8lnDsEVzd1d7BPi/76KgEfCv8WFwPLO9rkPSwh+zyOA/w6TXXnM8lsJqtBeVRBgJUHy7jEFjX5+B3wq3N/XgdslfTistktL3ht2igj/gB81s0N6UUadmRX2XVSuK5JOIrjJ/9GIQ0k7kp4lOHdpM66SpGsJGmD8NOpYBgKvgnPOORcJvwJyzjkXCb8Ccs45FwlPQM455yLhCcg551wkPAE555yLhCcg55xzkfj/NFQi1B1ZQowAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds.errorcurve(isodata_fe, calibrated_spike)\n",
    "plt.plot(prop, alpha_err, \"o\");"
   ]
  },
  {
   "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}