examples/Analysis_results_showcase.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T13:39:02.358041Z",
     "start_time": "2018-01-04T13:38:58.059324Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Configuration read from /Users/jburgess/.threeML/threeML_config.yml\n"
     ]
    }
   ],
   "source": [
    "from threeML import *\n",
    "\n",
    "from threeML.analysis_results import *\n",
    "\n",
    "from threeML.io.progress_bar import progress_bar\n",
    "\n",
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T13:39:07.340091Z",
     "start_time": "2018-01-04T13:39:07.332791Z"
    }
   },
   "outputs": [],
   "source": [
    "# These are the same simulated dataset we use in the test of the XY plugin\n",
    "\n",
    "x = np.linspace(0, 10, 50)\n",
    "\n",
    "poiss_sig = [44, 43, 38, 25, 51, 37, 46, 47, 55, 36, 40, 32, 46, 37, 44, 42, 50, 48, 52, 47, 39, 55, 80, 93, 123, 135,\n",
    "             96, 74, 43, 49, 43, 51, 27, 32, 35, 42, 43, 49, 38, 43, 59, 54, 50, 40, 50, 57, 55, 47, 38, 64]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T13:39:14.172126Z",
     "start_time": "2018-01-04T13:39:13.881517Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using Poisson log-likelihood\n",
      "Best fit values:\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "WARNING RuntimeWarning: invalid value encountered in log\n",
      "\n",
      "\n",
      "WARNING RuntimeWarning: invalid value encountered in log\n",
      "\n",
      "\n",
      "WARNING RuntimeWarning: invalid value encountered in log\n",
      "\n",
      "\n",
      "WARNING RuntimeWarning: invalid value encountered in log\n",
      "\n",
      "\n",
      "WARNING RuntimeWarning: invalid value encountered in log\n",
      "\n",
      "\n",
      "WARNING RuntimeWarning: invalid value encountered in log\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>result</th>\n",
       "      <th>unit</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>parameter</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>(7.800 +/- 3.2) x 10^-1</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>(4.040 +/- 0.18) x 10</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>(7.200 +/- 0.6) x 10</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>4.990 +/- 0.030</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>(3.310 +/- 0.27) x 10^-1</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                        result  \\\n",
       "parameter                                                        \n",
       "fake.spectrum.main.composite.a_1       (7.800 +/- 3.2) x 10^-1   \n",
       "fake.spectrum.main.composite.b_1         (4.040 +/- 0.18) x 10   \n",
       "fake.spectrum.main.composite.F_2          (7.200 +/- 0.6) x 10   \n",
       "fake.spectrum.main.composite.mu_2              4.990 +/- 0.030   \n",
       "fake.spectrum.main.composite.sigma_2  (3.310 +/- 0.27) x 10^-1   \n",
       "\n",
       "                                                  unit  \n",
       "parameter                                               \n",
       "fake.spectrum.main.composite.a_1      1 / (cm2 keV2 s)  \n",
       "fake.spectrum.main.composite.b_1       1 / (cm2 keV s)  \n",
       "fake.spectrum.main.composite.F_2           1 / (cm2 s)  \n",
       "fake.spectrum.main.composite.mu_2                  keV  \n",
       "fake.spectrum.main.composite.sigma_2               keV  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Correlation matrix:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table id=\"table4902340112\">\n",
       "<tr><td>1.00</td><td>-0.83</td><td>-0.01</td><td>-0.05</td><td>-0.01</td></tr>\n",
       "<tr><td>-0.83</td><td>1.00</td><td>-0.17</td><td>0.05</td><td>-0.12</td></tr>\n",
       "<tr><td>-0.01</td><td>-0.17</td><td>1.00</td><td>-0.04</td><td>0.41</td></tr>\n",
       "<tr><td>-0.05</td><td>0.05</td><td>-0.04</td><td>1.00</td><td>-0.08</td></tr>\n",
       "<tr><td>-0.01</td><td>-0.12</td><td>0.41</td><td>-0.08</td><td>1.00</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       " 1.00 -0.83 -0.01 -0.05 -0.01\n",
       "-0.83  1.00 -0.17  0.05 -0.12\n",
       "-0.01 -0.17  1.00 -0.04  0.41\n",
       "-0.05  0.05 -0.04  1.00 -0.08\n",
       "-0.01 -0.12  0.41 -0.08  1.00"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of -log(likelihood) at the minimum:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>-log(likelihood)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>test</th>\n",
       "      <td>175.306988</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>175.306988</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       -log(likelihood)\n",
       "test         175.306988\n",
       "total        175.306988"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of statistical measures:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>statistical measures</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AIC</th>\n",
       "      <td>361.977613</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BIC</th>\n",
       "      <td>370.174092</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     statistical measures\n",
       "AIC            361.977613\n",
       "BIC            370.174092"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# A simple test of the likelihood minimization\n",
    "\n",
    "y = np.array(poiss_sig)\n",
    "\n",
    "xy = XYLike(\"test\", x, y, poisson_data=True)\n",
    "\n",
    "fitfun = Line() + Gaussian()\n",
    "\n",
    "fitfun.a_1.bounds = (-10, 10.0)\n",
    "fitfun.b_1.bounds = (-100, 100.0)\n",
    "fitfun.F_2 = 60.0\n",
    "fitfun.F_2.bounds = (1e-3, 200.0)\n",
    "fitfun.mu_2 = 5.0\n",
    "fitfun.mu_2.bounds = (0.0, 100.0)\n",
    "fitfun.sigma_2.bounds = (1e-3, 10.0)\n",
    "\n",
    "model = Model(PointSource('fake',0.0, 0.0, fitfun))\n",
    "\n",
    "data = DataList(xy)\n",
    "\n",
    "jl = JointLikelihood(model, DataList(xy))\n",
    "_ = jl.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T13:39:20.252892Z",
     "start_time": "2018-01-04T13:39:20.106735Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>result</th>\n",
       "      <th>unit</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>parameter</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>(7.800 +/- 3.2) x 10^-1</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>(4.04 -0.18 +0.19) x 10</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>(7.200 +/- 0.6) x 10</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>4.990 +/- 0.030</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>(3.31 -0.26 +0.28) x 10^-1</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                          result  \\\n",
       "parameter                                                          \n",
       "fake.spectrum.main.composite.a_1         (7.800 +/- 3.2) x 10^-1   \n",
       "fake.spectrum.main.composite.b_1         (4.04 -0.18 +0.19) x 10   \n",
       "fake.spectrum.main.composite.F_2            (7.200 +/- 0.6) x 10   \n",
       "fake.spectrum.main.composite.mu_2                4.990 +/- 0.030   \n",
       "fake.spectrum.main.composite.sigma_2  (3.31 -0.26 +0.28) x 10^-1   \n",
       "\n",
       "                                                  unit  \n",
       "parameter                                               \n",
       "fake.spectrum.main.composite.a_1      1 / (cm2 keV2 s)  \n",
       "fake.spectrum.main.composite.b_1       1 / (cm2 keV s)  \n",
       "fake.spectrum.main.composite.F_2           1 / (cm2 s)  \n",
       "fake.spectrum.main.composite.mu_2                  keV  \n",
       "fake.spectrum.main.composite.sigma_2               keV  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This functionality is still there, although it cannot be propagated\n",
    "# (error propagation assume Gaussian errors, i.e., symmetric errors)\n",
    "# In this case though errors are pretty symmetric, so we are likely in the case\n",
    "# where the MLE is actually normally distributed\n",
    "\n",
    "_ = jl.get_errors()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T13:39:23.744571Z",
     "start_time": "2018-01-04T13:39:23.739234Z"
    }
   },
   "outputs": [],
   "source": [
    "# This is the novelty. This access the new class of the results\n",
    "ar = jl.results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T13:39:26.451134Z",
     "start_time": "2018-01-04T13:39:26.191875Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: AstropyDeprecationWarning: \"clobber\" was deprecated in version 2.0 and will be removed in a future version. Use argument \"overwrite\" instead. [astropy.utils.decorators]\n",
      "WARNING:astropy:AstropyDeprecationWarning: \"clobber\" was deprecated in version 2.0 and will be removed in a future version. Use argument \"overwrite\" instead.\n"
     ]
    }
   ],
   "source": [
    "ar.write_to(\"test_mle.fits\", overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T13:39:37.104132Z",
     "start_time": "2018-01-04T13:39:36.972821Z"
    }
   },
   "outputs": [],
   "source": [
    "ar_reloaded = load_analysis_results(\"test_mle.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T13:39:38.243632Z",
     "start_time": "2018-01-04T13:39:38.227669Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>-log(likelihood)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>test</th>\n",
       "      <td>175.306988</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>175.306988</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       -log(likelihood)\n",
       "test         175.306988\n",
       "total        175.306988"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ar_reloaded.get_statistic_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-01-04T13:39:40.632560Z",
     "start_time": "2018-01-04T13:39:40.586380Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>value</th>\n",
       "      <th>negative_error</th>\n",
       "      <th>positive_error</th>\n",
       "      <th>error</th>\n",
       "      <th>unit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>0.781166</td>\n",
       "      <td>-0.307275</td>\n",
       "      <td>0.316438</td>\n",
       "      <td>0.311856</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>40.361719</td>\n",
       "      <td>-1.849585</td>\n",
       "      <td>1.827499</td>\n",
       "      <td>1.838542</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>71.561162</td>\n",
       "      <td>-5.899793</td>\n",
       "      <td>5.817131</td>\n",
       "      <td>5.858462</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>4.989653</td>\n",
       "      <td>-0.029269</td>\n",
       "      <td>0.030078</td>\n",
       "      <td>0.029673</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>0.330522</td>\n",
       "      <td>-0.026752</td>\n",
       "      <td>0.027020</td>\n",
       "      <td>0.026886</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                          value  negative_error  \\\n",
       "fake.spectrum.main.composite.a_1       0.781166       -0.307275   \n",
       "fake.spectrum.main.composite.b_1      40.361719       -1.849585   \n",
       "fake.spectrum.main.composite.F_2      71.561162       -5.899793   \n",
       "fake.spectrum.main.composite.mu_2      4.989653       -0.029269   \n",
       "fake.spectrum.main.composite.sigma_2   0.330522       -0.026752   \n",
       "\n",
       "                                      positive_error     error  \\\n",
       "fake.spectrum.main.composite.a_1            0.316438  0.311856   \n",
       "fake.spectrum.main.composite.b_1            1.827499  1.838542   \n",
       "fake.spectrum.main.composite.F_2            5.817131  5.858462   \n",
       "fake.spectrum.main.composite.mu_2           0.030078  0.029673   \n",
       "fake.spectrum.main.composite.sigma_2        0.027020  0.026886   \n",
       "\n",
       "                                                  unit  \n",
       "fake.spectrum.main.composite.a_1      1 / (cm2 keV2 s)  \n",
       "fake.spectrum.main.composite.b_1       1 / (cm2 keV s)  \n",
       "fake.spectrum.main.composite.F_2           1 / (cm2 s)  \n",
       "fake.spectrum.main.composite.mu_2                  keV  \n",
       "fake.spectrum.main.composite.sigma_2               keV  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# You can get a DataFrame with the results\n",
    "ar_reloaded.get_data_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from threeML.analysis_results import AnalysisResultsSet\n",
    "\n",
    "analysis_set = AnalysisResultsSet([ar, ar_reloaded])\n",
    "\n",
    "analysis_set.set_bins(\"testing\", [-1, 1], [3, 5], unit = 's')\n",
    "\n",
    "analysis_set.write_to(\"analysis_set_test.fits\", overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_set = load_analysis_results(\"analysis_set_test.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best fit values:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Value</th>\n",
       "      <th>Unit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>(7.8 +/- 3.2) x 10^-1</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>(4.04 +/- 0.18) x 10</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>(7.2 +/- 0.6) x 10</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>4.990 +/- 0.030</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>(3.31 +/- 0.27) x 10^-1</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                        Value  \\\n",
       "fake.spectrum.main.composite.a_1        (7.8 +/- 3.2) x 10^-1   \n",
       "fake.spectrum.main.composite.b_1         (4.04 +/- 0.18) x 10   \n",
       "fake.spectrum.main.composite.F_2           (7.2 +/- 0.6) x 10   \n",
       "fake.spectrum.main.composite.mu_2             4.990 +/- 0.030   \n",
       "fake.spectrum.main.composite.sigma_2  (3.31 +/- 0.27) x 10^-1   \n",
       "\n",
       "                                                  Unit  \n",
       "fake.spectrum.main.composite.a_1      1 / (cm2 keV2 s)  \n",
       "fake.spectrum.main.composite.b_1       1 / (cm2 keV s)  \n",
       "fake.spectrum.main.composite.F_2           1 / (cm2 s)  \n",
       "fake.spectrum.main.composite.mu_2                  keV  \n",
       "fake.spectrum.main.composite.sigma_2               keV  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Correlation matrix:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table id=\"table140619854942160\">\n",
       "<tr><td>1.00</td><td>-0.83</td><td>-0.01</td><td>-0.05</td><td>-0.01</td></tr>\n",
       "<tr><td>-0.83</td><td>1.00</td><td>-0.17</td><td>0.05</td><td>-0.12</td></tr>\n",
       "<tr><td>-0.01</td><td>-0.17</td><td>1.00</td><td>-0.04</td><td>0.41</td></tr>\n",
       "<tr><td>-0.05</td><td>0.05</td><td>-0.04</td><td>1.00</td><td>-0.08</td></tr>\n",
       "<tr><td>-0.01</td><td>-0.12</td><td>0.41</td><td>-0.08</td><td>1.00</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       " 1.00 -0.83 -0.01 -0.05 -0.01\n",
       "-0.83  1.00 -0.17  0.05 -0.12\n",
       "-0.01 -0.17  1.00 -0.04  0.41\n",
       "-0.05  0.05 -0.04  1.00 -0.08\n",
       "-0.01 -0.12  0.41 -0.08  1.00"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of -log(likelihood) at the minimum:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>-log(likelihood)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>test</th>\n",
       "      <td>175.307004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>175.307004</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       -log(likelihood)\n",
       "test         175.307004\n",
       "total        175.307004"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "analysis_set[0].display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best fit values:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Value</th>\n",
       "      <th>Unit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>(7.8 +/- 3.2) x 10^-1</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>(4.04 +/- 0.18) x 10</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>(7.2 +/- 0.6) x 10</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>4.990 +/- 0.030</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>(3.31 +/- 0.27) x 10^-1</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                        Value  \\\n",
       "fake.spectrum.main.composite.a_1        (7.8 +/- 3.2) x 10^-1   \n",
       "fake.spectrum.main.composite.b_1         (4.04 +/- 0.18) x 10   \n",
       "fake.spectrum.main.composite.F_2           (7.2 +/- 0.6) x 10   \n",
       "fake.spectrum.main.composite.mu_2             4.990 +/- 0.030   \n",
       "fake.spectrum.main.composite.sigma_2  (3.31 +/- 0.27) x 10^-1   \n",
       "\n",
       "                                                  Unit  \n",
       "fake.spectrum.main.composite.a_1      1 / (cm2 keV2 s)  \n",
       "fake.spectrum.main.composite.b_1       1 / (cm2 keV s)  \n",
       "fake.spectrum.main.composite.F_2           1 / (cm2 s)  \n",
       "fake.spectrum.main.composite.mu_2                  keV  \n",
       "fake.spectrum.main.composite.sigma_2               keV  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Correlation matrix:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table id=\"table140619854919888\">\n",
       "<tr><td>1.00</td><td>-0.83</td><td>-0.01</td><td>-0.05</td><td>-0.01</td></tr>\n",
       "<tr><td>-0.83</td><td>1.00</td><td>-0.17</td><td>0.05</td><td>-0.12</td></tr>\n",
       "<tr><td>-0.01</td><td>-0.17</td><td>1.00</td><td>-0.04</td><td>0.41</td></tr>\n",
       "<tr><td>-0.05</td><td>0.05</td><td>-0.04</td><td>1.00</td><td>-0.08</td></tr>\n",
       "<tr><td>-0.01</td><td>-0.12</td><td>0.41</td><td>-0.08</td><td>1.00</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       " 1.00 -0.83 -0.01 -0.05 -0.01\n",
       "-0.83  1.00 -0.17  0.05 -0.12\n",
       "-0.01 -0.17  1.00 -0.04  0.41\n",
       "-0.05  0.05 -0.04  1.00 -0.08\n",
       "-0.01 -0.12  0.41 -0.08  1.00"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of -log(likelihood) at the minimum:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>-log(likelihood)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>test</th>\n",
       "      <td>175.307004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>175.307004</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       -log(likelihood)\n",
       "test         175.307004\n",
       "total        175.307004"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# or display again the results. However, this time I choose to display\n",
    "# hpd intervals\n",
    "ar.display(\"hpd\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Propagating a+b, with a and b respectively:\n",
      "equal-tail: (7.7 -3.1 +3.2) x 10^-1, hpd: (7.7 -3.4 +2.9) x 10^-1\n",
      "equal-tail: (4.04 -0.19 +0.18) x 10, hpd: (4.04 +/- 0.18) x 10\n",
      "\n",
      "This is the result (with errors):\n",
      "equal-tail: (4.12 -0.16 +0.15) x 10, hpd: (4.12 -0.15 +0.16) x 10\n",
      "(39.58120551740343, 42.71993633743719)\n",
      "\n",
      "This is 50 * a/a:\n",
      "equal-tail: (5.0 +/- 0) x 10, hpd: (5.0 +/- 0) x 10\n",
      "\n",
      "This is arcsinh(a + 5*b) / np.log10(b) (why not?)\n",
      "equal-tail: 3.738 -0.017 +0.019, hpd: 3.738 +/- 0.018\n",
      "\n",
      "Ratio a/b:\n",
      "equal-tail: (5.2 -1.7 +4) x 10, hpd: (5.2 -2.4 +1.8) x 10\n",
      "\n",
      "Propagating using a custom function:\n",
      "equal-tail: (7.7 -1.5 +2.0) x 10, hpd: (7.7 -1.8 +1.7) x 10\n"
     ]
    }
   ],
   "source": [
    "# You can use the results for propagating errors non-linearly for analytical functions\n",
    "p1 = ar.get_variates(\"fake.spectrum.main.composite.a_1\")\n",
    "p2 = ar.get_variates(\"fake.spectrum.main.composite.b_1\")\n",
    "\n",
    "print(\"Propagating a+b, with a and b respectively:\")\n",
    "print(p1)\n",
    "print(p2)\n",
    "\n",
    "print(\"\\nThis is the result (with errors):\")\n",
    "res = p1 + p2\n",
    "print(res)\n",
    "\n",
    "print res.equal_tail_confidence_interval()\n",
    "\n",
    "# The propagation accounts for covariances, so for example this\n",
    "# has error of zero (of course) since there is perfect covariance\n",
    "print(\"\\nThis is 50 * a/a:\")\n",
    "print(50 * p1/p1)\n",
    "\n",
    "# You can use arbitrary (np) functions\n",
    "print(\"\\nThis is arcsinh(a + 5*b) / np.log10(b) (why not?)\")\n",
    "print(np.arcsinh(p1 + 5*p2) / np.log10(p2))\n",
    "\n",
    "# Errors can become asymmetric. For example, the ratio of two gaussians is\n",
    "# asymmetric notoriously:\n",
    "print(\"\\nRatio a/b:\")\n",
    "print(p2/p1)\n",
    "\n",
    "# You can always use it with arbitrary functions like:\n",
    "def my_function(x, a, b):\n",
    "    \n",
    "    return b*x**a\n",
    "\n",
    "print(\"\\nPropagating using a custom function:\")\n",
    "print(my_function(2.3, p1, p2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "equal-tail: (1.31 -0.10 +0.12) x 10^2, hpd: (1.31 -0.13 +0.09) x 10^2\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAFkCAYAAACjCwibAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xl8nFd59//PpXVGkiVLsmTZsePdjpOYEDsrIUmzQNiS\nsjxtYkjTAqW0JZTHLWWnvz7QBWjBlKWFAm2hEEMJS1gTCKGQEJIQOziQeI2XOF5ka5dG++j8/jhz\nS6PRaDSSR5pF3/frNS9p7rnn1pFs6Vz3Ode5jjnnEBEREZlKUbYbICIiIvlBQYOIiIikRUGDiIiI\npEVBg4iIiKRFQYOIiIikRUGDiIiIpEVBg4iIiKRFQYOIiIikRUGDiIiIpEVBg4iIiKRl2kGDmV1t\nZt8xs+NmNmJmtyQ5Z6OZ3WNmHWbWY2aPmtmyuNfLzezTZtZiZt1mdreZNZ7tNyMiIiKzZyYjDZXA\nr4G3ABM2rjCzNcCDwNPANcAm4INAf9xpHwdeDrwmds5S4BszaIuIiIjMETubDavMbAR4pXPuO3HH\ndgCDzrk/nOQ91cAZ4Dbn3LdixzYAe4ArnHOPzbhBIiIiMmsymtNgZoYfQThgZveaWbOZPWJmvxt3\n2hagBPhJcMA5tw94Frgyk+0RERGRzCnJ8PUagSrgncB7gXcALwW+aWa/45x7EGjCj0R0Jby3Ofba\nBGZWD9wEHGH8NIeIiIikFgJWAvc551rP5kKZDhqCkYtvO+c+Efv8STN7AfCn+FyHyRhJciRibgK+\nkpkmioiIzEuvA+46mwtkOmhoAYbx+Qnx9gBXxT4/BZSZWXXCaEMjfrQhmSMAX/7yl9m4cWPmWisp\nbdu2je3bt2e7GfOKfuZzTz/zuaef+dzas2cPt99+O8T60rOR0aDBOTdkZr8CNiS8tB44Gvt8Jz6w\nuAEIEiHXA+cCv5zk0v0AGzduZPPmzZlssqRQU1Ojn/cc08987ulnPvf0M8+as57en3bQYGaVwFr8\ndALAajO7CGhzzh0D/gn4qpk9CPwUn9PwCuBaAOdcl5l9AfiYmbUD3cAngF9o5YSIiEjumslIwyX4\nYMDFHh+NHf8i8Abn3LfN7E+B9wD/AuwDXu2cix9F2AZEgbuBcuBefN0HERERyVHTDhqccz9jiqWa\nzrn/Av4rxesDwFtjDxEREckD2ntCJrV169ZsN2He0c987ulnPvf0M89fZ1URcq6Y2WZg586dO5U8\nIyIiMg27du1iy5YtAFucc7vO5loaaRAREZG0KGgQERGRtChoEBERkbQoaBAREZG0KGgQERGRtCho\nEBERkbQoaBAREZG0KGgQERGRtChoEBERkbQoaBAREZG0KGgQERGRtChoEBERkbQoaBAREZG0KGgQ\nERGRtChoEBERkbQoaBAREZG0KGgQERGRtChoEBERkbQoaBAREZG0KGgQERGRtChoEBERkbQoaBAR\nEZG0KGgQERGRtChoEBERkbQoaBAREZG0KGgQERGRtChoEBERkbQoaBAREZG0KGgQERGRtChoEBER\nkbQoaBAREZG0TDtoMLOrzew7ZnbczEbM7JYU5342ds5fJByvNbOvmFmnmbWb2efNrHIm34CIiIjM\njZmMNFQCvwbeArjJTjKzVwKXAceTvHwXsBG4AXg5cA3w2Rm0RUSyqLsbBgez3QoRmSsl032Dc+5e\n4F4AM7Nk55jZOcAngJuAHyS8dl7s+Bbn3BOxY28Fvm9mb3fOnZpum0Rk7g0NwT/+I0QisGkTrF8P\na9bAkiVQpIlPkYI07aBhKrFA4kvAR5xze5LEFVcC7UHAEHM/ftTicuCeTLdJRDLv2DH/KCqC++6D\nH/4QFiyA886DbdugvDzbLRSRTMt40AC8Cxh0zn1qktebgNPxB5xzUTNri70mInngyBHo7YULLwQz\nGBmBtjb47W/h+HFYvTrbLRSRTMvoIKKZbQH+Anj9TN5OihwJEcktBw5AcbEPGMCPONTXQ1+fDxpE\npPBkeqThhUADcCxuWqIY+JiZ/V/n3GrgFNAY/yYzKwZqgeZUF9+2bRs1NTXjjm3dupWtW7dmpvUi\nkpahIXj6aUj4dcTMP44dy067ROa7HTt2sGPHjnHHOjs7M3b9TAcNXwJ+nHDsR7Hj/xl7/ktgoZld\nHJfXcAN+pOHRVBffvn07mzdvzmBzRWQmjh3zUxFNSSYUKyr8KISIzL1kN9K7du1iy5YtGbn+tIOG\nWD2FtfhOHmC1mV0EtDnnjgHtCecPAaeccwcAnHN7zew+4HNm9mdAGfBJYIdWTojkhyCfoaJi4msL\nFsCJE3455oIFc940EZlFM8lpuAR4AtiJz0H4KLAL+H+TnJ8sT+G1wF78qonvAT8H3jyDtohIFiTm\nM8SrroauLuU1iBSimdRp+BnTCDZieQyJxzqA26f7tUUk+ybLZwiUlflzjh/3yy9FpHCoBIuITEuQ\nz7BwYfLXg2TIZ5+d23aJyOxT0CAi05IqnyFQVQX79oHTImqRgqKgQUSmJVU+Q2DBAmhp8SMSIlI4\nFDSISNqmymcILFigZEiRQqSgQUTSNlU+Q6C01JeVPnFibtolInNDQYOIpC2dfIZAcbE/X0QKh4IG\nEUlbOvkMgQULYP9+iEZnv10iMjcUNIhIWgYH08tnCCxYAO3tcObM7LZLROaOggYRSctzz6WXzxCo\nqoKeHiVDihQSBQ0ikpbp5DOAn8YYGVHQIFJIFDSISFqefRaKitLLZwiUlcHBg7PXJhGZWwoaRCQt\nJ05AKDS99yxYAIcP+/oOIpL/FDSIyJRGRqC5efpBQ3U1dHbCyZOz0y4RmVsKGkRkSp2dPp8hHJ7e\n+yoq/PuU1yBSGBQ0iMiU2tuhr2/6QUNRkd+0SkGDSGFQ0CAiU2pvh/5+KC+f/ntDIZ/XICL5ryTb\nDRCR3Nfe7j8Wxd1mHDv2S37xiw+zcOFK6urWjj5qa1djNnZiOOxzGpyb3soLEck9ChpEZErJtrhe\ntuwKzjnnMp588r/51a/+lZERv0Tiyivfzotf/E+j54XDvshTd7dPjBSR/KWgQUSmdPKk37kynplx\n9dXv4eqr38PISJSurmP8+Md/zYED3x8XNIRC0NHhRysUNIjkN+U0iMiUTp5MnQRZVFTMwoUr2bDh\nlbS07CESOT36WjjskyiTjVaISH5R0CAiKfX1+ZGCdGo0rFx5LUVFJTQ3Pzl6LMiDCPIiRCR/aXpC\nRFIKVk7U1099bnX1Mt75zg7KyiqTXkdE8ptGGkQkpba26dVoSBYwlJSoKqRIIVDQICIptbdDNOo7\n/pkKll2KSH5T0CAiKSVOKwwPD0z7GuGwH7Ho789Qo0QkKxQ0iEhKZ86MFWXq62vjIx+pZ//+70/r\nGqGQn+JQXoNIflPQICIpnTgxls9w8OB9DA1FaGp6/rSuEQ77UQYtuxTJbwoaRGRS0SicPj0WNBw4\n8D2ami6muvqcaV2npASGhzXSIJLvFDSIyKQ6OvzW1qEQjIwMc+DAD1m//hVpvdc5RzQ6CIxNbyho\nEMlvChpEZFJBjYZw2G9Q1d/fnnbQ8IUvXMkDD7xv9LmZH7UQkfyloEFEJtXeDgMDfkvs/fu/R2Vl\nI0uXXpLWe2trV3H06M9Gn4fDPj9CRPKXggYRmVR7ux8hMPP5DOvWvXzctteprFhxLSdO7GRgoBvw\nUxynT/s8CRHJTwoaRGRSra3gHHR2PsuZM0+zbt3L037vihXX4lyUY8ceBvxIQ2+vz5MQkfykoEFE\nJnXyJJSVQXX1ct7+9mbWr08/aFi06DwqKxtHpyiCZZdKhhTJX9MOGszsajP7jpkdN7MRM7sl7rUS\nM/uwmT1pZj2xc75oZksSrlFrZl8xs04zazezz5vZxIL1IpI1zo1tiW1mVFY2UlKSxlaXMWbGihXX\ncPTozwGfFzEwoFoNIvlsJiMNlcCvgbcALuG1CuD5wP8DLgZeBWwA7kk47y5gI3AD8HLgGuCzM2iL\niMyS3l7o6kpvS+zJrFhxLcePP8bQUK+WXYoUgGlvQeOcuxe4F8As+DMw+loXcFP8MTO7E3jUzJY5\n554zs42xc7Y4556InfNW4Ptm9nbn3KmZfSsikknt7b70c2PjzK+xYsU1jIwMcfz4Y6xc+TuAz5MQ\nkfw0FzkNC/EjEkH60xVAexAwxNwfO+fyOWiPiKQh2GDqbEYaGhou4M1vfoJly64A/BSFll2K5K+z\n2Ox2amZWDnwIuMs51xM73ASMK/HinIuaWVvsNRHJAe3tMDKSektsF5ugHD/mOKaoqHjcPhWhEJw6\n5d832XtEJHfNWtBgZiXA1/EjCH+ezluYmCMxzrZt26ipqRl3bOvWrWzdunWmzRSRSUyVexCNwu7d\nfnUF+I+VlVBTA1VVyd8TDvs8iZ4eWLAgs+0VEdixYwc7duwYd6yzszNj15+VoCEuYFgOXB83ygBw\nCmhMOL8YqAWaU113+/btbN68OcOtFZFkTp+GohQTmG1tsHgx3Hqrr73w7LP+cfgwrFkDFRUT3xME\nDe3tChpEZkOyG+ldu3axZcuWjFw/40FDXMCwGrjOOZd4v/JLYKGZXRyX13ADfqTh0Uy3R0RmJtgS\n+6GHPsTQUC/XXfeBca+fOQNXXQU33jh2bGAA3vUun+yYLGgIhXxyZVsbnHvuLH8DIpJxM6nTUGlm\nF5lZMFG5OvZ8eWzE4BvAZuB2oNTMFscepQDOub3AfcDnzOxSM7sK+CSwQysnRHLD0JAPCsJhOHDg\n+7S1HRz3ejTq8x0SB/7Ky/2xyao+Fhf7fAYtuxTJTzNZPXEJ8ASwE5+D8FFgF742wzLg5tjHXwMn\ngJOxj1fGXeO1wF78qonvAT8H3jyj70BEMq6jw48IhEJw5sweFi3aOO71tjaorYXzz5/43o0bfXAw\nNDT59VXgSSQ/zaROw89IHWxMGYg45zrwIxEikoOCLbEXLDhDX18rDQ3jo4OWFnjBC6C+fuJ716+H\nurqxnIdEJSXQnDJ7SURylfaeEJEJ2tthcBA6Op4GoKFhbKQhGoXh4YlTE4HqarjwwvFFnB5//DM8\n8sjHAT96oVoNIvlJQYOITBDkJLS27sGsmLq6taOvtbX5EYYLLpj8/Zs2jeU9ABw//hi/+c1dgM+T\naGvzSZMikl8UNIjIBEHQcObMHurr11FcXDb6WkuLH0moq5v8/evX+3oNwfLw+voNtLbuwzlHOOzz\nJZQMKZJ/FDSIyATNzVBaCi0tT49LgoxG/WOqcilNTbBypQ8wABYt2sDAQBeRSPO4ZZcikl8UNIjI\nBM3NPvdgw4bf5cILbxs93trqRxhSTU2ALxG9ZYsPDpzzIw0ALS37KC31OREaaRDJP7O694SI5J/B\nQd+hh0Jw2WV3jnutpQWuucYvt5zK+vW+wFNvL9TVrcGsmNbWfaxceS0weS0HEcldGmkQkXE6O/1y\ny/Ly8cejUT9qcPHF6V1n1So/TdHSAsXFZdTWrqKlZR/gRyI0PSGSfxQ0iMg4HR3Jt8RubZ28oFMy\nJSV+iqKryz+vr19Pa6sPGsrK/G6XIpJfND0hIuN0dPgpisSRhulMTQTOO88HD4ODsH79zfT0+KpO\noZDPm9AW2SL5RUGDiIwTLJNM7MxHRmDt2onnp7J+va/p0NoKl1zyp6PHQyG/PXZfX/KNrUQkN2l6\nQkTGSZagGI36ICJZWehUKivhec+buFIiFPJTIEqGFMkvChpEZJzTp/2UQrxIxAcA0w0aANatG6sM\nGVDQIJKfFDSIyDg+QbGFo0cfJBodBHzQUFUFDQ3Tv15DAxQV+doMgZIS/1xBg0h+UdAgIqOGh33+\nQUvLT/iv/7qGwcEewOcfLF/uq0ROV0PDWL2GQJAvEeRPiEh+UNAgIqOCGg3d3XuorGwkHPYbTAwM\n+LoLM7FokR+liEQmvqZaDSL5RUGDiIwKajR0dj5NQ4MvyOCcf62paWbXLC2FJUvGjzSAr9XQ3HwW\njRWROaegQURGdXb6UYW2tj2jG1UNDvqOv7Fx5tddscIvrxwZidLefojh4X5CIZ8/EQQlIpL7FDSI\nyKiODnDO0d5+kLq6dcBYEuRMVk4EgveePv1bPvGJNZw8uYtQyFeL7O/PQMNFZE4oaBCRUX564hTD\nw/3U1a0BfBJkbS0sXDjz6zY0+OTHhQt9INLSso9QyI9qKBlSJH8oaBCRUWfOQH//IQAWLvSZj5GI\nT4I8m3LPwQqKoaEKamrOpbV1n2o1iOQhBQ0iMqq5GYaGngOgttYHDdGoX255NhYt8sWhIhGor99A\na+s+Skt9voSCBpH8oaBBRAAfHJw5A+eddyvveU+EsrKq0UqOZ5PPAH7zq6Ymv4Kivn4DLS37MPOj\nFwoaRPKHggYRAXxSYl+fL/FcWup3kert9dMKZxs0AKxc6a+/aNEG2toOMjISBSbuSyEiuUtBg4gA\nYzUaQqGxY8GeE2ez3DKweLFfXllXt5aRkSG6uo5RWur3uhCR/KCgQUSAsRoN8UFDTw8sXTr+2EwF\n+1bU1PhVGe3th0ZrNYhIfiiZ+hQRmQ98jQa/uVSgr2/m5aMTBSsoystXceed+1m4cCUtLf7rDgz4\nvAcRyW0aaRARYPKExCVLMnP9hgY/1dHXV0J9/TqKi0u17FIkzyhoEBHA724ZX4thaAiKizOTBAl+\niqOxcfweFEHQoAJPIvlBQYOIAD63IDEJ8mzLRydasWJ80FBW5oMTjTSI5AcFDSLCyIhfxbB//0f5\n/vffAvgkyJoaqK/P3NdZsmT8BlXByIaCBpH8oERIEaG7248AtLQ8hFkf4EcaNm0anxh5toIVFCMj\nY9d1TrUaRPKFRhpEZHS5ZXf3IWprVwN+2mDFisx+nYYGCIf9qoxAWZkvXy0iuU9Bg4jQ0QF9fY7O\nTh80OOdHADKZzwBjKygikbFj5eUq8CSSL6YdNJjZ1Wb2HTM7bmYjZnZLknM+YGYnzKzXzH5sZmsT\nXq81s6+YWaeZtZvZ582s8my+ERGZOV8NsoXBwR5qa1fT1+dHBDIdNFRU+M2renvh/vvfzW9+s4NQ\nCNra/OZVIpLbZjLSUAn8GngL4BJfNLN3AncCbwYuAyLAfWZWFnfaXcBG4Abg5cA1wGdn0BYRyYCO\nDujtPQxAbe3qjJaPTrRypQ8ajhx5gEOHfqRllyJ5ZNpBg3PuXufc3zjnvg1YklPeBnzQOfdd59xv\ngTuApcArAcxsI3AT8Ebn3OPOuYeBtwK3mVnTTL8REZm5tjaIRA4BsHDhKnp7/YhAVVXmv9aSJT4R\nsrZ2DW1tz6jAk0geyWhOg5mtApqAnwTHnHNdwKPAlbFDVwDtzrkn4t56P37U4vJMtkdE0nPqFPT3\nHyIcriMUqqG3F845Z3a+VrCCorZ2De3thygv91MTChpEcl+ml1w24Tv/xFzo5thrwTnj0p6cc1Ez\na4s7R0TmiHM+EbGp6RIWLfprAKLRzJWPTtTQ4ItIVVauprv7OMPDfUBY0xMieWCu6jQYSfIfpnvO\ntm3bqKmpGXds69atbN269exaJzKP9fT41Qxr176YRYtePHo8k0Wd4o2toPC7XXZ0HAbO10iDSAbs\n2LGDHTt2jDvWmcGIPNNBwyl857+Y8aMNjcATceeMS68ys2KglokjFONs376dzZs3Z6yxIhIst4QF\nC/zzYM+J2QoaKiuhrg7a233Q0Nb2DCUl52vZpUgGJLuR3rVrF1u2bMnI9TOa0+CcO4wPCm4IjplZ\nNT5X4eHYoV8CC83s4ri33oAPNh7NZHtEZGrt7T4RMRz2z3t7/dLI2QoazGD5cjBbQklJiPb2Q4RC\nPq9CRHLbtEcaYvUU1jK2cmK1mV0EtDnnjgEfB95nZgeBI8AHgeeAewCcc3vN7D7gc2b2Z0AZ8Elg\nh3NOfzZE5lh7u1/NUFzsnwc1GurqZu9rNjVBNFrExRe/kZqacwmF/C6bw8NQouL2IjlrJr+elwA/\nxecfOOCjseNfBN7gnPuImVXg6y4sBB4EXuqciy/d8lrgU/hVEyPA3filmiIyxxL3fejthdWrfXnn\n2RIEJC972adG29DT46dKFi2ava8rImdn2kGDc+5nTDGt4Zz7W+BvU7zeAdw+3a8tIpl35sz4Tan6\n+/30wWyqr/dfMxhZCIehpUVBg0iu094TIvPcyZN+CWRgZCTz5aMT1dWN37iqvNxvmKXdLkVym4IG\nkXlseNiPNPT27qaz89joRlWzlQQZSAwaLJYhpaBBJLcpaBCZxzo7fQ7Dz3/+Bn7+8w8yOOjv+mc7\naAiHobZ2/BbZ4MtZi0juUtAgMo8Fyy27uw9TW7ua3l7foc920AC+THV80FBW5qdKRCR3KWgQmcfa\n26G7u53+/vbRLbErK2Hhwtn/2kuW+OmRQCjkgwY3Ve1YEckaBQ0i81hHB0QiY1ti9/bC0qVjNRtm\nU12dDxAGByN0dh4jHIaurolTFiKSOxQ0iMxj8Vti19aupr8fli2bm68dLLv84Q//L1/72qsIhXzA\noGRIkdyloEFkHgu2xC4vryYUqgXGtq6ebXV1vlx1VdUa2tufIRz2+RUKGkRyl4IGkXns5EkfNNTW\nrsY5v+5xLpIgYWzZZUXFGvr7OxgaaiMaVdAgkssUNIjMU319QU7DIWpr19Df75MR5ypoCId9wmV5\n+WrA73bpnIIGkVymrWFE5qmODj8d8JrX/IBQqIdIZHZ3t0xk5pMuDxzwW2S3tz9DUdGlnDkzN19f\nRKZPIw0i81RQo6GysoRQaCG9vVBdDQsWzF0bli4Fs4WEw3W0tfm8hhMn5u7ri8j0KGgQmafa28dv\nRd3X5wsuBSWd50IwqlFb65MhQyE4fRqi0blrg4ikT0GDyDwV5A4EQcLgoA8a5lJdnf/6CxeOraDo\n7fXlrUUk9yinQWSeammZOKow19tSBysorrrqQyxYUEZJiW9Xe7t/TURyi4IGkXnq1Cm/ORX46QCz\nuUuCDARBQ0nJChYs8BUitUW2SO7S9ITIPDQy4oOGcNg/7+vzn8/13X1FBdTUaItskXyhoEFkHuru\nhkjE12UAn0cwl8stA2YTd7t0Tltki+QqBQ0i81B7O/T2Rvne917C4cM/pa8Pamt94DDXzjkHhobG\nnmuLbJHcpaBBZB5qb4e2tuc4cuQ+hof76e2F5cuz05bEKZGgVoO2yBbJPQoaROYhXz76GQDq6tYy\nPAxLlmSnLcGyy6A2Qyjkt8ju789Oe0RkcgoaROah9nbo6TmIWRELF64A5j6fIRCsoAiChHBYW2SL\n5CoFDSLzUHMz9PY+Q03NCkZGfH2EbAcNzz33W7797T/ErFtbZIvkKAUNIvPQyZPQ13eQuro19PX5\nBMhsFVOqrPR7XvT0RNi9+0t0dx9keFhBg0guUtAgMs8MDfmqi5HIM9TWrs1ajYZAsOyytHQdAG1t\nB2Ifs9MeEZmcggaReaajA/r6HF1dYyMNixaNVYfMBr/bZR3hcB2trQcoKkJbZIvkIAUNIvOMr9Ew\nyMUX/xnnnvvC0d0ts6m+3o841NWto63tgLbIFslR2ntCZJ5pb4dotJyXvOSfMIPf/hYWL85um4Kp\nkdratbS1HRi3RXZxcXbbJiJjNNIgMs8kbontXPZ3lAxWUFRXr6Ot7eDoskttkS2SWxQ0iMwz8QmG\nQ0NkdblloL7er+CorFxHJHIasy7VahDJQQoaROaZ48fHkh6zvXIiUFnpd7usqno+F164FbOItsgW\nyUEKGkTmkZERHzQEG1PlStBgBsuWQXn5+bzmNXdRU+NrWitoEMktChpE5pH2dr8tdnzQUF+f3eWW\ngcTdLkFBg0iuyXjQYGZFZvZBMztkZr1mdtDM3pfkvA+Y2YnYOT82s7WZbouIjNfSAr2944OGZcuy\n26ZAYl5Faam2yBbJNbMx0vAu4M3AnwPnAe8A3mFmdwYnmNk7gTtj510GRID7zKxsFtojIjG+EmQ7\nbW27GRkZJhqFpqZst8qrr4eiIhge9s/DYT+Voi2yRXLHbAQNVwL3OOfudc4965z7JvAjfHAQeBvw\nQefcd51zvwXuAJYCr5yF9ohITEsLnD79Y/7935/PwEAXkP18hkCwgqKvzz+vrPTTE1p2KZI7ZiNo\neBi4wczWAZjZRcBVwA9iz1cBTcBPgjc457qAR/EBh4jMkhMnIBI5SChUS2lpHcXFuRU0hMN++gR8\n0BCJqJy0SC6ZjaDhQ8DXgL1mNgjsBD7unPtq7PUmwAHNCe9rjr0mIrPAOXj2WRgYeGZ0z4lcWDkR\nKC+HhoaxkYbychgc9JUhRSQ3zEbQcCvwWuA24GLgD4G/NrM/mOJ9hg8mRGQW9Pb6wk6RyEHq6rK/\nu2Uyy5b5dvb2ttDX14pzChpEcsls7D3xEeAfnHNfjz1/ysxWAu8G/hs4hQ8QFjN+tKEReCLVhbdt\n20ZNTc24Y1u3bmXr1q0ZabhIIQtWTnR1PcOaNVfT2+uXOYZC2W7ZmKYmX0viX//1Ai699C0sXvw3\nHDuW7VaJ5I8dO3awY8eOccc6M5gYNBtBQwUTRwxGiI1qOOcOm9kp4AbgSQAzqwYuBz6d6sLbt29n\n8+bNGW+wyHxw5gxEIn309BwfHWnIleWWgWDZZbBx1erVfkrFubG9MkRkcslupHft2sWWLVsycv3Z\nCBq+C7zXzI4BTwGbgW3A5+PO+TjwPjM7CBwBPgg8B9wzC+0REfxIQ0/PIQBqa9fQ3Z07yy0D9fW+\nPkNt7TpaW/dSWQkdHX4VRS5No4jMV7OR03AncDd+1OBp/HTFvwF/E5zgnPsI8Engs/hVE2Hgpc65\nwVloj4jgcwMikWcAqKtbG/uYzRZNFKygqKpaR1vbAa2gEMkxGR9pcM5FgL+MPVKd97fA32b664tI\nckePwqpVN3PFFScIhZooKsq9oKG21i+1rKhYR19fG9FoG0NDdZw+DRs2ZLt1IqK9J0TmgaEhaG6G\nykpjwYIl9PcbFRW5FzSUlPgpk/LydQC0tx/UCgqRHKKgQWQeaG31w/y5trtlMsuXQ2mpnz5pbd1P\nWZlPhhSR7FPQIDIPJNuoqq7OBw65pqEBSkoWUF29jNOnn6Ky0gcNIyPZbpmIzMbqCRHJMS0tEI36\nlQngg4bQ+eBpAAAgAElEQVRzzslumyYTLLu8/fb7qak5h0gEurp8YapFi7LbNpH5TiMNIvNAS8v4\n58PDubfcMlBf70tIV1VtoKysSisoRHKIggaReeDYMShL2Hg+F/MZYOJul6WlPpFTyZAi2aegQaTA\nOeeDhiCfYXgYiorGpgFyTU0NLFgwtttlUAmyOXGLOxGZcwoaRApcZyd0d8NDD93O7t1fyumVE+CD\nhKVLx0YawE9XaAWFSPYpaBApcGfOQE/PMAcPfp3+/s6cDxrA74kxMDD2PFhBEY1mr00ioqBBpOC1\ntEBr6zNEo4M0Nl5AX5+vvJiLyy0DiaskKiuhp8fXmxCR7FHQIFLgWlqgq+tpABoaLsjp5ZaB+no/\nTRGMLGgFhUhuUNAgUuBOnoSenqcIh+uorGxkeBiWLMl2q1JbtMiPhJw8eYCvfe3V9PUdZ3hYKyhE\nsk1Bg0iBe/ZZiESeoqHhAsBwDhobs92q1ILdLgcHy9i791s0N+/GTCsoRLJNQYNIAevrC6YnfNAw\nOOhXIuR60BAkahYXn0t5eTXNzb+hvNzv1Cki2aOgQaSAtbT4lROdnftoaDifSMTnB+R60GDmV1D0\n9RmNjRdy+vSTVFbCc8/5OhMikh0KGkQKWEsLRCIDXHXVe1m58neIRGDhQv/IdUuX+gChsXETzc2/\nGV1BkVgSW0TmjoIGkQJ26hQUF1dy3XV/w+LFm+jpgVWrxqos5rKgYmVj4yZaWvYSCg1pBYVIlilo\nEClghw/7ktGBaNQP++eDpia/X0Zt7SZGRobo7NzHyIiSIUWySUGDSIGKRuHAAaiu9s+d8x9zPZ8h\n0NQEVVUQDm8CoLn5N5ipnLRINiloEClQzc3Q3u43fwK/kiIUgsWLs9uudC1YAA0NMDxcyxVX/CU1\nNeeyYAHs2QMjI9luncj8pKBBpEAdP+4TB6uq/PN8WTkRMIM1a3y7b7rpo5x77lVUV/tESBV5EskO\nBQ0iBer4cT8lUVzsn0ciPmCorMxuu6Zj2bKxaRXwow/d3X6rbxGZewoaRArUgQO+kFMgEoGVK7PW\nnBlZvNgncga1GYqLfRChvAaR7FDQIFKABgbgyBEYGHiCEyceB3xnm+sbVSUKkiF7esaOhcOwd2/2\n2iQynyloEClAJ05AVxc89dSH+NGP3k406nME8iWfIbBokV/9ER80VFf7kYbu7uy1S2S+UtAgUoCO\nH/erJdra/J4Tvb1QUZF/QUNxMaxe7adWAjU1PiB67rnstUtkvlLQIFKAjh2DaLSf1tb9NDZeMLpy\noqEh2y2bvhUrYGho7Hl5uX+uZEiRuaegQaQA7d8PfX07GRkZYtmyK4lEfD5DfGJkvmhq8h8jkXYO\nHrwX5xzFxfDMM9ltl8h8pKBBpMB0dwc5DQ9TWlrJ4sWb6OvLv5UTgaYmn/z4zDMP85WvvJSOjsNU\nV8O+fdrxUmSuKWgQKTDHj/s5/5aWh1m27HKKikqAsTv2fLN4sV9BsXDhCwDj6NGfU13tq12ePJnt\n1onMLwoaRArM8eMwOOg4ceJhli17AUNDPqEwX8pHJwqFfJGnoaFampou4ujRn40uw1Reg8jcUtAg\nUmCOHoWhoRaKi8tZvvzKvCsfncyaNX41yIoV13LkyP9i5peQHj2a7ZaJzC8KGkQKiHM+CbK+voFt\n255l7dqXEon42gb19dlu3cwtWeK/txUrrqWj4widnc9SWQlPPz2+zLSIzK5ZCRrMbKmZ/beZtZhZ\nr5ntNrPNCed8wMxOxF7/sZmtnY22iMwnra1w5szYdthmRiQCy5eP7UGRj5qaoKwMmpquBuDIkZ9R\nXQ2nTvncBhGZGxkPGsxsIfALYAC4CdgI/BXQHnfOO4E7gTcDlwER4D4zK8t0e0Tmk+PH/eqJYDts\ngMHB/F05EQjKSY+MLKKxcRNHjvwv1dXQ2am8BpG5NBsjDe8CnnXO/bFzbqdz7qhz7n7n3OG4c94G\nfNA5913n3G+BO4ClwCtnoT0i88bx4zAyAiV+wcTo0H0+5zOAHzlZtMgnP65adQMDA52UlkI0qqBB\nZC7NRtBwM/C4mf2PmTWb2S4z++PgRTNbBTQBPwmOOee6gEeBK2ehPSLzxuHD46chBgf9sH6+Bw1m\nPhmypwduuulj/P7v3w34723//iw3TmQemY2gYTXwZ8A+4MXAZ4BPmNntsdebAAc0J7yvOfaaiMzA\n8LDfDjvIZwDfyVZW5u9yy3jLlvmREzMbPVZd7StDDgxksWEi88hsBA1FwE7n3Pudc7udc/8OfA4f\nSKRi+GBCRGbg+HGfFBgfNEQiUFvrN3nKd01NfsQhvgpkTQ10dGjppchcKZmFa54E9iQc2wO8Ovb5\nKXyAsJjxow2NwBOpLrxt2zZqEv76bd26la1bt55Ne0UKwr590NU1wtKlfUAl4JMir77ad7b5LkiG\njETGgqBw2NdveOopWL8+u+0TyQU7duxgx44d4451dnZm7PqzETT8AtiQcGwDcBTAOXfYzE4BNwBP\nAphZNXA58OlUF96+fTubN29OdYrIvOQcPPEE9Pf/lg9/eDNvfOMvWbr0UpyDtQWymLmhwQcLPT1j\nQYOZH1l55BG4+eaxBFCR+SrZjfSuXbvYsmVLRq4/G9MT24ErzOzdZrbGzF4L/DHwqbhzPg68z8xu\nNrNNwJeA54B7ZqE9IgWvpcXnM/T3PwxAY+MFdHf7fIZ8X24ZKC6G1av96Em8xka/guLAgey0S2Q+\nyXjQ4Jx7HHgVsBX4DfBe4G3Oua/GnfMR4JPAZ/GrJsLAS51zg5luj8h8sHevn9tvbX2YJUsuprS0\ngs5Ov0xx6dJsty5zzj/frwiJrwJZUQH9/fDkk9lrl8h8MSsVIZ1zP3DOPc85V+Gcu8A59x9Jzvlb\n59zS2Dk3OecOzkZbROaD3/wGiorguef8JlXgd7q84IL8rgSZaN06X7gqGG3o6joOwMKF8NhjWkUh\nMtu094RInuvu9kFDKNRMe/szLF/+ApzzRZ4KJZ8hsGyZ34eirQ327Pkm27cvo6enmcZGv3pk375s\nt1CksCloEMlz+/f7PSd6e38JwPLlL6Cnx680KJR8hkBxMWze7AOlZcuuAODo0Z8RCsHQkKYoRGab\nggaRPPfUU35U4cSJX1BdvYyamuV0dPh8hnPOyXbrMm/9er9Korx8KXV16zhy5GeAr0fx2GN+CaaI\nzA4FDSJ5bHDQL7WsqYHDh+9n5crrgMLMZwisXQt1dX6KYtWqGzhw4Ps4N0JjIzQ3w57EKjEikjEK\nGkTy2MGDvqNctAjuuOMBrr/+73COgqrPkKiqygdEbW1w0UV30Nl5lEOH7qeszG9g9UTKEnEicjYU\nNIjksb17/XLDcBjC4Vpqas4d3W+i0PIZ4l1wgQ8QzjnnChoaLmDXrs8BUF8PO3dOrOUgIpmhoEEk\nT42MwOOP+yWI8YJ8hkKqz5BobOmlsXnzm9i79x4ikdMsWgRnzsDTT2e7hSKFSUGDSJ569ll47jl/\ndx2vq8sXQSrkkspLl/rll62tcNFFf0A4XMepU7spLfWv79qV3faJFCoFDSJ5at8+vw9D/EhDodZn\nSFRUBBdf7DevCofr+Mu/fI41a14E+LLSv/qVz/cQkcxS0CCSh6JRv7wwFBq/g+V8yGcIrF8PpaW+\nCmRR0diwyqJFfrTl7rt97QYRyRwFDSJ5aOdOX58hsQ5DIddnSLR2rZ+aaWsbf9zMb2y1cyc89FB2\n2iZSqBQ0iOSZwUH4wQ+Czw/h4nZv6u72KwsKOZ8hEA7D8543MWgAv4lVRQV8+9s+70FEMkNBg0ie\nefRRvzqgsbGFT31qA7t3fxGYP/kM8c4/f+z7TnTuuT5Z9J57xu+KKVKoRkb8Euzubh8snzwJR4/6\n34NMmQf3IyKFo6/PjzKUlsIzz3wd5xzr1r0M8PP48yWfIbBuna+G2d4+cRVJcTEsXw4//Slceils\n2pSdNookE436fJzEx9CQH00cHPSfDwz4z4OP/f3+70Bfn/88eASvj4z4a8c/zpzJXLsVNIjkkV/8\nwm9QtWEDPPjgl1m79iYqKxsBOHECLrnE32HPF4sX+w2sHnhgYtAAvtx0czN84xs+wAiF5r6Nkv+i\n0fEdd+IjsaOP/zgw4Dv4SAR6e/0j6PCjURge9o/Ejj4QJDrHj5YVF489ioomfl5ePv6cZCNxM6Wg\nQSRPdHfDD3/oRxMikUMcO/Ywr371XYD/AwRw7bXjV1MUOjO48UY/ZdPRAQsXQmvrfh5++J952cs+\nTXFxKatX+63D77kHXv1qRms5SOFxbuwOPlnnnqqzHxwc69CDTr23d+wuPujMgw4+8Y7ebPzvXlDO\n3cx35iUlvgMPPgafBx18/PHi4sz+Hmcyx0lBg0ieePBBOHIENm6E733vA4RCtWzYcAvgRxlWroSL\nLspqE7Ni/Xo/2vDQQz5o6O1t5Ykn/oPKykauv/7vKC+Hpib4+tf9qMMf/IHfEVOyz7mJHXf8Ixiy\nT/Z5b6+/ew86+aCjHxqaeNce/xgZmbyDj79rj38EHXpZWfLXi4rmT7CuoEEkD7S3w733+vn7Y8ce\nYPfuL3LzzZ+nrKySaNTXZ7j+en/XMt8Eow2/+pXP61i+/Equu+4DPPDA+1i16gZWrbqOhga/2uKn\nP/XJYa9/vQ82ZHqCpNOhIf8I7sLj59bjh/GDj0GnHj8839fnX082NJ84TJ+sg5+sYy8q8qNJodDk\nHbzMnIIGkRzX3w9f+xocOwYbNvTxuc+9mRUrruHii98A+Lvnpia47LIsNzSLzj8fnv98Hzicfz5c\nddU7OXTofr71rdv50z/dTUXFotHdMfftg3/+Z7jtNj+dU4jbhweCufj4Djx+zj1+3j3xbj4SGXsE\nd/HxnfzIyNh8fPy8fHAXHwzNB3f2yTrw4DHZHXxJyfy5g88XChpEclhfH3zhCz7Rb9Uq/4d048bX\n8Pzn/xFmhnPQ0gL/5//M7yH3oiI/2rBrlx91qaoq5lWv+m8+85mLuOeeN3DbbfdgZpSU+KDi2Wfh\ns5+F+++Hyy/3x1avzm4A4ZzveBMT7qaanw/u9IO7+OCRzlz8yMjYsHx8ol1wJ584zx7cqQd385WV\n/pz489TJFzYFDSI5KhKBz30Ofv5z36FVVQGEufHGD42e09EB1dVw5ZVZa2bOeN7z4MILfdLjeedB\ndfU5/O7v/idf/eot/OpXn+ayy+4EfKe2YoX/+R4/Dl/8oi8Ede65fvVJY6NfdVFX5wOxYMon6GwT\nO93gebLEumR38MFdfOJQfbL5+MSh++DuPfg+4pPtkiXZpbqTD+bv1cnnPudGGB7uZ2ioN+HRx7nn\nXpXyvbt2fYE9e76fsbYoaBDJQd3d8O//7pP71q71d3TJnDgBV101v2ozTKa4GF70InjySd8Jh8Ow\nYcPNXHbZW/nRj/6KFSuuYfHi542eX1npR2+c8534c8/5qQsz3+GGw/4RCo113kEnHdyhxz9PllEP\nE+fjYWLCXXwHn2w+vqRk7A5f8sfTT9892rlP7PB7ufjiN7Bs2eWTvn/fvu9w9923MjzcP+k573//\n0Li9VxK1tu6nvX3vWX0f8RQ0iOSYo0fhq1+FRx7xtQUqKpKf19fnO5NrrtHdYuDii/3qkn37fC0L\ngBe96CNUVS2hoeGCpO8x86M4fiRnbJogGAHo6hrLjg8+Bp14MH8fdOiJd/P6d8mOwcEe2tsPE40O\nMDzcH7tL74t97j9efPEbUna2jz76CZ555kcMDfUyPDy+0x8cjLB69Q38/u9/I2U7vvnN24lGBzAr\norS0YtyjpCTMeee9MuX7Gxou4MYbPzzhvfEPs9SR5Ite9GHWrr2VL31pS8rz0qWgQSQHOOeXU/7k\nJ/Dww37aYf16f6c7meeegzVrVOkwXmkp3HSTDxpaW33Bp5KSEFdf/e60r2Hmr1Na6qd+ZPqcc0Sj\ngwwNRWKddd/ox9LSChobL0z5/nvv3cbAQMe4Tnqs8+7j+uv/ngsu+L1J33/06IPcddfLUn6NCy/c\nSnn5gklfHxmJUlRUQlXVYkpKKigtDcc66krKyiqpr596+c3b395MaWmYoqJSbAYRZF3dGi6//C+m\n/b7ZpKBBJIsGB+HwYfjf//UjC11dsGSJL3+c6m9MV5ef/77xRhUrSnT55XDzzb4KZFkZLJi8X5AE\nvb2tHDv28IQ7a99x9zA4GOHFL/4nSksnGf4CfvCDt/L44/+Kc8nLEK5adQN33HF/ynacOrWLkZHh\n0btp33GHR+/QFy5ckfL9y5dfyRvf+AglJeUUF5dTWhqmpCRESYn/WFxcNmUnfuWV27jyym0pz5lK\nKFRzVu/PRQoaRObIyIjP7G9u9qMKBw7Anj2+BkNfHyxd6pPxgr9le/few6pV10+4GwoCjZtugquv\nnvvvI9cVFcHv/R50dsKPfjT1iA3AoUP3U1Nzblp3j7NlcDDCwEBnks46wtBQhPLyGtate2nKa3z9\n679HT8+pCUPxwfD8i1/8US699M8mff/p07/lq1+9ZfR5UVHpaMddVlZFWVklQ0N9KYOGjRtfTWPj\nBZSWVsbeG451+P5jODz1Mp8/+qOfTXlOKqHQwpS5AjJzChpEpiGY7w6WucUXtkkscjMw4IOEtjY/\nVN7e7o9FIv714mJfwXDJEp+3EH/j8/DD/8yPf/zX3HTTx7niireNHh8Z8UPvF10Er3tdYdcYOBul\npXDHHX5E5pFH/JLKVCMy9977Ns6ceZpwuI7KysVUVjZSVbWYiopGSksrWL/+FaxYMXmE1tFxlMce\n+yTR6BAjI0OMjAyPdtZB5/2a1+ygomLRpNf4yU/ew2OPfWLS15ctu2LKoCEcXjR6Nx5/Zx08li9/\nQcr3L1t2BX/1V6dGO/tUc/6TWbXqOlatum7a75P8oKBBClpQpjZx3XtiSdpk69+Dte6J1eyCZXGJ\nhW2GhycuhwM/RB48yst9oBAOJ59+aGs7yIMP/gO//vV/8sIXvmfCfOahQz7IeMMbNOw+lcpKeNOb\nfOD29NO+sNNkqw/e9KbHOXjwXlpb99HT00wk0kwkcpozZ/YwPNxHff26lEFDf38H+/d/l6KiUoqL\nSykqKonrrMOEQrW4Kfbn3rLlTaxZ82LKyipJTJgrK6tKeXcfeMUr/m3Kc1IpKSmnqmrxWV1DCpuC\nBskpk9Win2y9e/y5fX2+g+jpGVsHH+wkl2zt+/Dw+II28bvJBcVt4jPi4zecKS31HX9Q2CZ4bSbZ\n8s45jhz5Xx55ZDv793+Piop6Xv7yz3DJJW8ed15zs7/+HXfMr50sz0ZdHbz5zbB9Ozz1lF++mmyq\norQ0zMaNr5rx12lquog779x3Fi2FxsYLp0wQFElH4tLg4eHMXVtBg8xIUEI2KEWbatOZZI/4gjbx\nG84EgcBkdegT7+ZhrOOPL3CTuKNc4lax8XXqs7kszjnHF75wBcePP0Zj44XcfPPn2LTptZSWju/Z\nurvhzBnYuhUuvTRLjc1Ty5bBW98Kd93lK0bW1PhjWg4pMxXfISd+DD4PHvHPk32e+HGy1+L/vybW\n/ogf2Yw/J/j71tGRue9dQUOBCraIjR86j+/ck9WfT9bxx1esi98HPujYJytqM9V2scEdfLBlbOIa\n98QKdvlUpnZkJEpfXxuRyGkGB7tZtuyKSc81M1auvJ7rr/8HVq26fkJGt3N+g6X2drjuOrjlltz/\n/nPRypXw9rf7ctzf/ravGjlWZVNyzXQ716k63VTHgs42sdpmsjbFnxdfnyO+g052PP55aelYEa/S\n0rGbnODz+I/xj8SS3tN5vncv3HdfZv5tFDTMsWjUd7jJaszH37Un3sXHP49PxEu2ZWzQ+Sf+ssSX\nvp2sNG28+M498T9kcLy0dOJdfL5sF+ucY3CwJ0nxlr7Y+vJeli+/KuUc7/793+eRRz7G0FAvAwNd\nRCKn6e1tBfxfmHC4nne8oyVlO2688R+THh8e9issKiv9ds4vfamWV56NsjJ4yUt8qemvfx0efdT/\njBcu9NMYk+WZFKLEO+LEznOyTja+bHaqTjvZ3XcgnY4ZJnbGyTroxGPxI43xnW3wsaxs4vHg8+AG\nJv7vV7K9NxJHMeNHLac6N1ubcPX2Zu5aChpignmfZHfl8Y/J7tTjt4oNMueTJdsND09+Vx5fyz5x\nKCr+FxrGR6/xm8jEP0+MchP/Y8cfyxUjI9HRZWKpMs0B9u79Ni0t+2IdfSQWAERi68l7WLXqBl74\nwndN+v6BgS4+/OGFKb/G6173Q9aufcmkr5eUhKisbKSkpILy8gVUVDRQWdlIZWUDFRUNLFiwJPU3\nPImuLr+s8rzz4PbbfRKfZMayZfAXfwE7d/o8h927/R4UfX2+U6mo8GWcg0dp6cz/yCcbrk5195yq\nTPVkzyfrsJN1zomdd/zdcfxjss45eASJvfGdblnZWEccvBZ8TPY3KlkHG9+ZT7Zh1lQdea7frOS7\nvAoaOjrg1KmJHXWyYfb4z5Ml0cXfkQd3/cmi6cRH/Px5IPEXMf4XZLLPp7o7z+X/9M6NcPr0U+OW\nlPnOOjLaga9ffzN1dWsmvcb+/d/jgQfeN3p+cJcfjQ4CUFIS5r3vTR0eP/nklzly5Kdx68gXxNaS\nVxEO1xEO16d8f1lZJa95zVfjMtXHiscEGezl5alLAq5efQOrV9+Q8px0OeeDhRMn/P+DF7/Y5zAs\nTB3XyAwUF/utxC+7zP/uHz0KBw/65awnT/oaD21tY1NxMPmccqLE3SITO9xkHXKyDjyYikt1h5v4\nCIX835by8vF30okdfLK/O4mPxGnD+E5dW1bPX7MeNJjZu4G/Bz7unPvL2LFy4GPArUA5cB/w5865\n06mu9S//4ocRE4fZJ0sUSRzqmuyXIzECLi2dPBrOpQ7dl2odYHjY11cPPq+vX5fyfb/5zV2cPPnE\n6JB8fInXoaFeVqy4luuv/7sUX3eEz3zmeZO+XlpaQW3tmpRBQ2VlI+ee+8Kk9diDDnsqv//7d095\nTipFRSVceOGtZ3WNTIhG/cqIlhY/x755M1x7rU94VB2G2VdW5vf4WLfOTwGBH3Xo6PCPSGRs9DH4\nCMnvzhP/liTr2FNN9yWOAOZDDo/ML7MaNJjZpcCbgN0JL30ceCnwGqAL+DTwDSBlfTszv1VtLv5i\nOedwLpqyGIpzbnTXs/GFX8YKwFx00R0sWbJ50mvs2/ddvvOdN47emQdz5/H+5m+iKTcxOXbsYQ4d\n+vG4Km3BXXVVVRPV1ctTfq9FRSW88Y2/HFc4JujoS0pCU26gAnDOOZdxzjmXTXleoRke9vOLwdLQ\ngQH/f7ehAV79arjiCr+fRLb/P893wQ6XS2Y2uyRSsGYtaDCzKuDLwB8D7487Xg28AbjNOfez2LHX\nA3vM7DLn3GOTXbOiYuYbyAwPD9DX1zp6Vx7cXQd32oODES688NaUHd5DD32Igwd/mGSbUz+8vmHD\nLdx667cmfb+Z8c1vvo6RkSHA4uqhh0bro69e/aKUf6jq69dxxRX/d/SuPOj0S0rKR68zlZe97FNT\nnjOVVCsCClWy+eNkiaWJBZ8GB8euUVQ09v94wwafwb9kif+8durquiIiWTWbIw2fBr7rnHvAzN4f\nd/yS2Nf9SXDAObfPzJ4FrgQmDRp27/40hw5VxzrpYA7df7z66vewfv0rJm3MoUM/ZseOm1M2eMOG\nWygrq5z09bKyKqqrl4/rsIMdz0pLK6mrW5vy+gB//denKSkJk86GKcksWnQeV1/9nmm/L18ly/JO\ntYQqPpk0yEFJZ3nWZBndqZZZJSZ1BfPAoZC/S62s9FUbGxt9XkJNjX/U1fljWgkhIvlmVoIGM7sN\neD4+QEi0GBh0znUlHG8GmlJd98iRH1BRUTu6cUppaSXhcB3V1csIhVJniy1deilbt36P4uKyCRuo\nlJaGKStbMOU8+mWX3clll92Z8pypTNXOXDFZZx2/fDNZJ52YTJosa3yy/JPEPJTEjnqqZVfxnXhQ\nsjk+GztZdnfi8cR558mWV022jjoU8tfT9IKIFKKMBw1mtgyfs/Ai59zQdN5Ksgn6OLW166msHL/V\n6IYNr2TTpq1TXryqajHr1798Gs3JrlRLsiZbM51s6Way8yF5VcWg005c1jlZsZLEBK7y8vGddfyy\nq+DzIKs7VWc81RrnyZZfxXfkWnYlIvPRjh072LFjx7hjnZ2dGbu+TbWJyrQvaPa7wDeBKD4QACjG\nBwRR4CXA/cDC+NEGMzsCbHfO/UuSa24Gdm7dupP16ydPEsykye6007mbnuox9n2Nfa3Er52swliy\nu+rEyorBkqv4TZISO+74ZVeT3UHHf0zWWSfeZQcdtYiI5JZdu3axZcsWgC3OuV1nc63ZmJ64H9iU\ncOy/gD3Ah4DjwBBwA/AtADNbD5wL/DLVhaNRn1SW7O55ukPl8ZUQk3Xe8R33ZJ104rH4O+zgERQ9\nCdZOxw+ZJw6Jp7s0K9l7grLLursWEZHZkvGgwTkXAZ6OP2ZmEaDVObcn9vwLwMfMrB3oBj4B/CLV\nygnw9fcPHx5/5z1ZpbFkQ+Xx2xPHV3sL7sATy4om68AT78LjP9edtoiIFLK5qgiZOAeyDT9VcTe+\nuNO9wFumusif/Als2pS89GiyIXURERHJnDkJGpxz1yc8HwDeGnukbfVqv55dRERE5p4G1EVERCQt\nChpEREQkLQoaREREJC0KGkRERCQtChpEREQkLQoaREREJC0KGkRERCQtChpEREQkLQoaREREJC0K\nGkRERCQtChpEREQkLQoaREREJC0KGkRERCQtChpEREQkLQoaREREJC0KGkRERCQtChpEREQkLQoa\nREREJC0KGkRERCQtChpEREQkLQoaREREJC0KGkRERCQtChpEREQkLQoaREREJC0KGkRERCQtChpE\nREQkLQoaREREJC0KGkRERCQtChpEREQkLQoaREREJC0KGkRERCQtChpEREQkLQoaZFI7duzIdhPm\nHUoLXTUAAAY0SURBVP3M555+5nNPP/P8lfGgwczebWaPmVmXmTWb2bfMbH3COeVm9mkzazGzbjO7\n28waM90WOTv6xZ57+pnPPf3M555+5vlrNkYargY+CVwO3AiUAj8ys3DcOR8HXg68BrgGWAp8Yxba\nIiIiIhlSkukLOudeFv/czP4IOA1sAR4ys2rgDcBtzrmfxc55PbDHzC5zzj2W6TaJiIjI2ZuLnIaF\ngAPaYs+34IOVnwQnOOf2Ac8CV85Be0RERGQGMj7SEM/MDD8V8ZBz7unY4SZg0DnXlXB6c+y1ZEIA\ne/bsmZV2SnKdnZ3s2rUr282YV/Qzn3v6mc89/cznVlzfGTrba5lz7myvMfnFzf4NuAl4oXPuROzY\nVuA/nHPhhHMfA+53zr0nyXVeC3xl1hoqIiJS+F7nnLvrbC4wayMNZvYp4GXA1UHAEHMKKDOz6oTR\nhkb8aEMy9wGvA44A/bPQXBERkUIVAlbi+9KzMisjDbGA4XeBa51zhxJeqwbO4BMhvxU7th7YC1yh\nREgREZHclPGgwcz+FdgK3ALsj3up0znXH3fOS4HXA93AJ4AR59zVGW2MiIiIZMxsBA0j+NUSiV7v\nnPtS7Jxy4J/xwUU5cC/wFufc6Yw2RkRERDJmVhMhRUREpHBo7wkRERFJi4IGERERSUteBA1m9hYz\nO2xmfWb2iJldmu02Fap0NhyT2RP7+Y+Y2cey3ZZCZ2ZLzey/Yxvn9ZrZbjPbnO12FSozKzKzD5rZ\nodjP+6CZvS/b7SokZna1mX3HzI7H/o7ckuScD5jZidi/wY/NbO10vkbOBw1mdivwUeD/Ay4GdgP3\nmdmirDascKWz4ZjMglgw/Cb8/3GZRWa2EPgFMIAvQLcR+CugPZvtKnDvAt4M/DlwHvAO4B1mdmdW\nW1VYKoFfA28hyYIEM3sncCf+3+EyIILvT8vS/QI5nwhpZo8Ajzrn3hZ7bsAx4BPOuY9ktXHzQCw4\nOw1c45x7KNvtKVRmVgXsBP4MeD/whHPuL7PbqsJlZh8CrnTOXZvttswXZvZd4JRz7k1xx+4Gep1z\nd2SvZYUptpLxlc6578QdOwH8k3Nue+x5Nb6o4h865/4nnevm9EiDmZXiN7iK39zKAfejza3mSuKG\nYzI7Pg181zn3QLYbMk/cDDxuZv8Tm4bbZWZ/nO1GFbiHgRvMbB2AmV0EXAX8IKutmifMbBV+f6f4\n/rQLeJRp9KezumFVBiwCiplYXroZ2DD3zZlfJtlwTDLMzG4Dng9cku22zCOr8aM6HwX+Hj8d9wkz\n63fOfTmrLStcHwKqgb1mFsXftL7XOffV7DZr3mjC3wAm608n2yxyglwPGiZjJC8gJZn1r8D5+LsB\nmQVmtgwfmL3IOTeU7fbMI0XAY86598ee7zazC/CBhIKG2XEr8FrgNuBpfKD8L2Z2wjn331lt2fw2\nrf40p6cngBYgCixOOJ5qcyvJgLgNx37HOXcy2+0pYFuABmCnmQ2Z2RBwLfA2MxuMjfZI5p0E9iQc\n2wOcm4W2zBcfAf7ROfd159xTzrmvANuBd2e5XfPFKXyAcFb9aU4HDbE7r53ADcGx2B/RG/DzYzIL\n4jYcu84592y221Pg7gc24e+6Loo9Hsff7V7kcj1TOX/9golTnBuAo1loy3xRwcQ72hFyvB8qFM65\nw/jAIb4/rcZPzaXdn+bD9MTHgC+a2U7gMWAb/j/ff2WzUYUqYcOxiJkFUenohmOSOc65CH6odpSZ\nRYBW51zinbBkznbgF2b2buB/8H84/xi/5FVmx3eB95rZMeApYDP+7/nns9qqAmJmlcBa/IgCwOpY\nwmmbc+4Yfir0fWZ2EDgCfBB4Drgn7a+RDzcyZvbn+DW9i/FrUN/qnHs8u60qTOlsOCazy8weAH6t\nJZezy8xehk/OWwscBj7qnPuP7LaqcMU6tA8Cr8IPiZ8A7gI+6JwbzmbbCoWZXQv8lIl/w7/onHtD\n7Jy/Bf4EvzLuQfxmkQfT/hr5EDSIiIhI9mkuSURERNKioEFERETSoqBBRERE0qKgQURERNKioEFE\nRETSoqBBRERE0qKgQURERNKioEFERETSoqBBRERE0qKgQf7/dutAAAAAAECQv/UKAxRFALBIAwCw\nBNVNaCNZGnIbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe4f117ee50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This is an example of an error propagation to get the plot of the model with its errors\n",
    "# (which are propagated without assuming linearity on parameters)\n",
    "\n",
    "def go(fitfun, ar, model):\n",
    "    \n",
    "    # Gather the parameter variates\n",
    "    \n",
    "    arguments = {}\n",
    "\n",
    "    for par in fitfun.parameters.values():\n",
    "\n",
    "        if par.free:\n",
    "\n",
    "            this_name = par.name\n",
    "            \n",
    "            this_variate = ar.get_variates(par.path)\n",
    "            \n",
    "            # Do not use more than 1000 values (would make computation too slow for nothing)\n",
    "            \n",
    "            if len(this_variate) > 1000:\n",
    "                \n",
    "                this_variate = np.random.choice(this_variate, size=1000)\n",
    "            \n",
    "            arguments[this_name] = this_variate\n",
    "    \n",
    "    # Prepare the error propagator function\n",
    "    \n",
    "    pp = ar.propagate(ar.optimized_model.fake.spectrum.main.shape.evaluate_at, **arguments)\n",
    "    \n",
    "    # You can just use it as:\n",
    "    \n",
    "    print(pp(5.0))\n",
    "    \n",
    "    #Make the plot\n",
    "    \n",
    "    energies = np.linspace(0, 10, 100)\n",
    "\n",
    "    low_curve = np.zeros_like(energies)\n",
    "    middle_curve = np.zeros_like(energies)\n",
    "    hi_curve = np.zeros_like(energies)\n",
    "\n",
    "    free_parameters = model.free_parameters\n",
    "    \n",
    "    with progress_bar(len(energies), title=\"Propagating errors\") as p:\n",
    "        \n",
    "        with use_astromodels_memoization(False):\n",
    "        \n",
    "            for i, e in enumerate(energies):\n",
    "\n",
    "                this_flux = pp(e)\n",
    "                \n",
    "                low_bound, hi_bound = this_flux.equal_tail_confidence_interval()\n",
    "                \n",
    "                low_curve[i], middle_curve[i], hi_curve[i] = (low_bound, this_flux.median, hi_bound)\n",
    "\n",
    "                p.increase()\n",
    "\n",
    "    plt.plot(energies, middle_curve, '--', color='black')\n",
    "    plt.fill_between(energies, low_curve, hi_curve, alpha=0.5, color='blue')\n",
    "\n",
    "go(fitfun, ar, model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Mean acceptance fraction: 0.5499\n",
      "\n",
      "Maximum a posteriori probability (MAP) point:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Value</th>\n",
       "      <th>Unit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>(8.1 -3.5 +2.8) x 10^-1</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>(4.03 -0.17 +0.20) x 10</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>(7.1 -0.5 +0.7) x 10</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>4.987 -0.029 +0.031</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>(3.31 -0.27 +0.30) x 10^-1</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                           Value  \\\n",
       "fake.spectrum.main.composite.a_1         (8.1 -3.5 +2.8) x 10^-1   \n",
       "fake.spectrum.main.composite.b_1         (4.03 -0.17 +0.20) x 10   \n",
       "fake.spectrum.main.composite.F_2            (7.1 -0.5 +0.7) x 10   \n",
       "fake.spectrum.main.composite.mu_2            4.987 -0.029 +0.031   \n",
       "fake.spectrum.main.composite.sigma_2  (3.31 -0.27 +0.30) x 10^-1   \n",
       "\n",
       "                                                  Unit  \n",
       "fake.spectrum.main.composite.a_1      1 / (cm2 keV2 s)  \n",
       "fake.spectrum.main.composite.b_1       1 / (cm2 keV s)  \n",
       "fake.spectrum.main.composite.F_2           1 / (cm2 s)  \n",
       "fake.spectrum.main.composite.mu_2                  keV  \n",
       "fake.spectrum.main.composite.sigma_2               keV  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Values of -log(posterior) at the minimum:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>-log(posterior)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>test</th>\n",
       "      <td>-176.698279</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>total</th>\n",
       "      <td>-176.698279</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       -log(posterior)\n",
       "test       -176.698279\n",
       "total      -176.698279"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Exactly the same can be done for a Bayesian analysis\n",
    "# Let's run it first\n",
    "\n",
    "for parameter in ar.optimized_model:\n",
    "    \n",
    "    model[parameter.path].value = parameter.value\n",
    "\n",
    "model.fake.spectrum.main.composite.a_1.set_uninformative_prior(Uniform_prior)\n",
    "model.fake.spectrum.main.composite.b_1.set_uninformative_prior(Uniform_prior)\n",
    "model.fake.spectrum.main.composite.F_2.set_uninformative_prior(Log_uniform_prior)\n",
    "model.fake.spectrum.main.composite.mu_2.set_uninformative_prior(Uniform_prior)\n",
    "model.fake.spectrum.main.composite.sigma_2.set_uninformative_prior(Log_uniform_prior)\n",
    "\n",
    "bs = BayesianAnalysis(model, data)\n",
    "\n",
    "samples = bs.sample(20, 100, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "ar2 = bs.results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "ar2.write_to(\"test_bayes.fits\", overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ar2_reloaded = load_analysis_results(\"test_bayes.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(ar2_reloaded.samples, ar2.samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>value</th>\n",
       "      <th>negative_error</th>\n",
       "      <th>positive_error</th>\n",
       "      <th>error</th>\n",
       "      <th>unit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.a_1</th>\n",
       "      <td>0.810175</td>\n",
       "      <td>-0.352269</td>\n",
       "      <td>0.281939</td>\n",
       "      <td>0.317104</td>\n",
       "      <td>1 / (cm2 keV2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.b_1</th>\n",
       "      <td>40.266060</td>\n",
       "      <td>-1.714355</td>\n",
       "      <td>1.968567</td>\n",
       "      <td>1.841461</td>\n",
       "      <td>1 / (cm2 keV s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.F_2</th>\n",
       "      <td>70.624725</td>\n",
       "      <td>-4.970925</td>\n",
       "      <td>6.588271</td>\n",
       "      <td>5.779598</td>\n",
       "      <td>1 / (cm2 s)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.mu_2</th>\n",
       "      <td>4.987446</td>\n",
       "      <td>-0.029302</td>\n",
       "      <td>0.030792</td>\n",
       "      <td>0.030047</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fake.spectrum.main.composite.sigma_2</th>\n",
       "      <td>0.330793</td>\n",
       "      <td>-0.026754</td>\n",
       "      <td>0.029881</td>\n",
       "      <td>0.028317</td>\n",
       "      <td>keV</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                          value  negative_error  \\\n",
       "fake.spectrum.main.composite.a_1       0.810175       -0.352269   \n",
       "fake.spectrum.main.composite.b_1      40.266060       -1.714355   \n",
       "fake.spectrum.main.composite.F_2      70.624725       -4.970925   \n",
       "fake.spectrum.main.composite.mu_2      4.987446       -0.029302   \n",
       "fake.spectrum.main.composite.sigma_2   0.330793       -0.026754   \n",
       "\n",
       "                                      positive_error     error  \\\n",
       "fake.spectrum.main.composite.a_1            0.281939  0.317104   \n",
       "fake.spectrum.main.composite.b_1            1.968567  1.841461   \n",
       "fake.spectrum.main.composite.F_2            6.588271  5.779598   \n",
       "fake.spectrum.main.composite.mu_2           0.030792  0.030047   \n",
       "fake.spectrum.main.composite.sigma_2        0.029881  0.028317   \n",
       "\n",
       "                                                  unit  \n",
       "fake.spectrum.main.composite.a_1      1 / (cm2 keV2 s)  \n",
       "fake.spectrum.main.composite.b_1       1 / (cm2 keV s)  \n",
       "fake.spectrum.main.composite.F_2           1 / (cm2 s)  \n",
       "fake.spectrum.main.composite.mu_2                  keV  \n",
       "fake.spectrum.main.composite.sigma_2               keV  "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ar2.get_data_frame(\"equal tail\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "equal-tail: (7.8 +/- 3.2) x 10^-1, hpd: (7.8 -3.1 +3.3) x 10^-1\n",
      "equal-tail: (4.03 -0.18 +0.19) x 10, hpd: (4.03 -0.21 +0.15) x 10\n",
      "equal-tail: (4.11 -0.15 +0.16) x 10, hpd: (4.11 -0.17 +0.14) x 10\n"
     ]
    }
   ],
   "source": [
    "p1 = ar2.get_variates(\"fake.spectrum.main.composite.a_1\")\n",
    "p2 = ar2.get_variates(\"fake.spectrum.main.composite.b_1\")\n",
    "\n",
    "print(p1)\n",
    "print(p2)\n",
    "\n",
    "res = p1 + p2\n",
    "\n",
    "print(res)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "equal-tail: (1.30 -0.10 +0.11) x 10^2, hpd: (1.30 -0.11 +0.08) x 10^2\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAFkCAYAAACjCwibAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xmc5FV97//X6equqt632XeGWViGZaZlMwyISFAkrgky\niAvG5SbgNagx4aq/n4nXXH94FeISw5UsKjoGiYmoXEBQFFBmAsPqjMPsPXtPb9X7fn5/nPpWfau6\nuru6p6pr6ffz8ahHd3/rW9863TB1Pt9zPudzjLUWERERkamU5LoBIiIiUhgUNIiIiEhaFDSIiIhI\nWhQ0iIiISFoUNIiIiEhaFDSIiIhIWhQ0iIiISFoUNIiIiEhaFDSIiIhIWhQ0iIiISFqmHTQYYzYb\nYx40xhw1xowZY94yybn3RM/570nH640x3zPGRIwxHcaYe40xlTP5BURERGR2zGSkoRJ4AbgVmHDj\nCmPM24CLgaMpnv4+cDZwNfBm4Argnhm0RURERGZJ6XRfYK19GHgYwBhjUp1jjFkKfBW4Fngo6bmz\nosebrLXPR499FPiZMeaT1toT022TiIiIZF/GcxqigcR3gDuttbtSnHIZ0OEFDFGP4UYtLsl0e0RE\nRCQzpj3SkIa/BoastV+f4PlFQIv/gLV21BjTHn1uHGNMI2504iAwkLmmioiIFL0wsAp4xFrbdjoX\nymjQYIxpAv47sHEmL2fiHIlrge/NtF0iIiLCu3E5hTOW6ZGGy4H5wGFfukMA+Iox5i+stauBE8AC\n/4uMMQGgHjg5wXUPAtx3332cffbZGW6yTOT222/nrrvuynUz5pRC+5sPDkIwCKmzmwpDof3Ni4H+\n5rNr165d3HzzzRDtS09HpoOG7wA/Tzr2aPT4v0R//i1QZ4zZ6MtruBo30rBtgusOAJx99tls2rQp\nsy2WCdXW1urvPcsK6W/e1QVf+AIsWwavex2ccw4EArlu1fQV0t+8WOhvnjOnPb0/7aAhWk9hDa6T\nB1htjLkAaLfWHgY6ks4fBk5Ya/cAWGt/b4x5BPiWMebPgCDwNWCrVk6IFI5Dh+DgQdi3D377Wzj3\nXHj962HTJgiFct06EcmGmYw0vAb4JS7/wAJfjh7/NvCBFOenylO4Cfg6btXEGPAA8LEZtEVEcqS5\nGUZG4LzzoKcHXn4ZXnwRmprgU58q7CkLEUltJnUafsU0lmpG8xiSj3UCN0/3vUUkf+zdC2Vl7vuq\nKjjrLGhthT17oKMDGhpy2z4RyTztPSET2rJlS66bMOcUyt98cNAFB7W1icdratyow/HjuWnXTBTK\n37yY6G9euBQ0yIT0D3v2Fcrf/OhR6Ox0QYJfMOimLE4UUHZSofzNi4n+5oVLQYOITFtzM/T3Q0XF\n+OeMgWPHZr9NIpJ9ChpEZNoOHoSSktTJjhUVbupCRIqPggYRmZaxMdi50yU/plJVBSdPutwGESku\nChpEZFpaWtwqieQkSE9VlQsYCimvQUTSo6BBRKbl8GHo7obq6tTPh8MwMFBYKyhEJD0KGkRkWpqb\n3RTFRCWjjXEPBQ0ixUdBg4hMy65dUF4++TnhMOzfPzvtEZHZo6BBRNLW3e2mJybKZ/BUVbkRicHB\n2WmXiMwOBQ0ikrbDhyESGV/UKZmXDHlyos3uRaQgKWgQkbQ1N8Pw8NS7WFZUQF+fVlCIFBsFDSKS\ntn37oDSNbe5KSsBaJUOKFBsFDSKSluFh2L176qkJT1kZHDqU3TaJyOxS0CAiaTl2zG1SNVUSpKe6\n2q2gGB3NbrtEZPYoaBCRtDQ3uzyFysr0zq+qckmTp05lt10iMnsUNIhIWryphlSbVKVSWQm9vUqG\nFCkmChpEJC0HDqTeCnsipaWucqSSIUWKh4IGEZnS6KjbqCocnt7rSkpcbQcRKQ4KGkRkSp2dLp9h\nqvLRySor3TJNa7PTLhGZXQoaRGRKHR1u58rpjjRUV0Nbm3u9iBQ+BQ0iMqWODrePxHSDBq+ctJIh\nRYqDggYRmZI3UpDuyglPMOiKQiloECkOChpEZErt7TN/rTEKGkSKhYIGEZnS8eNu1GAy3d3Heeih\njzI01JtwPBhU0CBSLBQ0iMikvI2n/PkM27d/neee+1bCeS0tr/D88/fyzDN3JxwPh7VFtkixUNAg\nIpPq74eursTlltu3f53jx59LOO/MM69h1aqrOHz46YTj4bArJz0wMButFZFsUtAgIpPq6HCBgzfS\n0NV1lLa23ZxxxtXjzl269BKOHt2O9RVmCIfdyovOztlqsYhki4IGEZmUV6PBG2k4cOBxAM4446px\n5y5dejH9/W10dOyLHQuH3esVNIgUPgUNIjKpjg5XRrq01P184MDjLFp0IRUV88adu3TpxQAcObIt\ndqysDIaGFDSIFAMFDSIyKX81R2st+/c/zhlnvCHluRUVjTQ0rOHo0e2xY8a4h4IGkcKnoEFEJnXq\nVLyoU1vbbrq7j7J69fh8Bo/La9g27rhKSYsUvtJcN0BE8tuxY/F8hv37H6ekpIwVKzZPeP7GjX9K\nd/exhGOq1SBSHBQ0iMiEkrfEXrBgA6973ecIBisnfE2qBEnVahApDgoaRGRCkYjbEruqyv28atWV\nrFp15bSv46/VMN1Nr0QkfyinQUQmNNMtsZNp2aVIcVDQICITam+f2ZbYyUIhFzREIplpl4jkhoIG\nEZnQTLfETuZtka2RBpHCpqBBRCZ0Olti+3lBh4IGkcKmoEFEJpTOltjToVoNIoVNQYOIpJRqS+zp\nOHz4Nzz//L/Efi4r07JLkUI37aDBGLPZGPOgMeaoMWbMGPMW33Olxpj/zxjzkjGmJ3rOt40xi5Ou\nUW+M+Z4xJmKM6TDG3GuMmXjht4jMOv+W2Hv3PsKePQ9N6/V79z7CY499KrbjZTisAk8ihW4mIw2V\nwAvArYBNeq4CuBD4G2Aj8HZgPfDjpPO+D5wNXA28GbgCuGcGbRGRLPFvib19+1d59tl/nNbrly69\nmL6+Vjo7DwDuOp2dbjWGiBSmaRd3stY+DDwMYExiTrW1tgu41n/MGHMbsM0Ys8xae8QYc3b0nCZr\n7fPRcz4K/MwY80lrre5FRPKAf0vstrY9rFt3/bRe7+14efTodurrV8eChs5OWLgwGy0WkWybjZyG\nOtyIhJc3fSnQ4QUMUY9Fz7lkFtojImnwtsQuKRmhs/MADQ1rp/X6ysr51Nevjm2TrQJPIoUvq0GD\nMSYEfBH4vrW2J3p4EdDiP89aOwq0R58TkTzgrXTo7DzI2NgIjY3TCxogccfLYBCGhhQ0iBSyrO09\nYYwpBX6IG0H483RewvgciQS33347tbW1Cce2bNnCli1bZtpMEZmAtyV2W9segGmPNICboti160eM\njg4TCJRhjKpCimTT1q1b2bp1a8KxSAb/0WUlaPAFDMuB1/tGGQBOAAuSzg8A9cCkC7LuuusuNm3a\nlOHWikgq3pbYp07tIRAIUVu7fNw5HR3unImWZS5ZchGjo4OcOrWTRYsuiL1GRLIj1Y30jh07aGpq\nysj1Mx40+AKG1cBV1trkj4jfAnXGmI2+vIarcSMN2zLdHhGZPv+W2G1te2hoOBNjEmczu7vh6FE3\n7eAlTDY0QF0dBALunHnz1rNixeWMjQ0DqtUgUuimHTRE6ymswXXyAKuNMRfgchKOAf+OW3Z5PVBm\njPHypNuttcPW2t8bYx4BvmWM+TMgCHwN2KqVEyL5wb8ldklJgCVLXjPunNZWWLsW3vMe2LcPXnwR\nDh50gcR557mpjYqKedxyy5Ox14TDChpECtlMRhpeA/wSl39ggS9Hj38bV5/hj6LHX4ge93IVrgJ+\nHT12E/B13KqJMeAB4GMzaIuIZIG33HLePHjjG+9OeU53N7zmNXDOOe5x/fWwezf83d9BTw9UV49/\nTTjsrj046Ha+FJHCMpM6Db9i8lUXU67IsNZ2AjdP971FZHZ4HftEuQr9/a7TP+us+DFjYM0aN0UR\niUwcNEQi7rFgwfjnRSS/ae8JERlnqi2x29pcgaY1axKPl5a6QKKrK/XrVKtBpLApaBCRcSbq9D2d\nnbBpU+ophjVrXCKlTbGAWrUaRAqbggYRGae9HUom+HQYHnYjEOeck/r5FStcMDEwMP45b+RCQYNI\nYVLQICLjnDo1caJiW5tLkFy/PvXzy5e7ZZeT1ZNR0CBSmBQ0iEgCa91yyomChvZ2t6QyVaIjuHoN\nZ56ZGDRYO8bwcD/gajVoi2yRwqSgQUQS9PdDb68LGmxSYsLoKIyNuaBhMuvXu9wFz733XsKjj34S\ncMmQLS0TvFBE8pqCBhFJ0NXlllsGAoPceWcju3b9KPZcZ6ebevAvtUxlxQq3ksILHGpqltPRsRdw\nQUN7e2JQISKFQUGDiCTwgob+/v0MDHRQXt4Qe661FdatczkNk1m5Empr46swGhrW0N4eDxoGBrRx\nlUghUtAgIgm6utwKie7uxN0trXXBxIUXTly/wVNbC8uWxQOD+voz6ew8xOjocGxlhZIhRQqPggYR\nSeCNDrS376GsrILq6iWx49XVU09NeM45x+1fAW6kwdpRIpFDhEKq1SBSqBQ0iEiCeNCwl4aGNZjo\nsEJrq8tVWLYsveusXOlqPYyOuqDBu6Y3StHdnemWi0i2KWgQkQTt7W76ob19T2xqAtyoQVPTxEWf\nkq1Y4UYmurqgpmYZJSVlsbwGUE6DSCFS0CAiCbzCTv6gYWzMBRJLlqR/nQUL3CMScdtr19evpr19\nH+AFJdlovYhk00y2xhaRImWtCxoCgQEikcM0Nrqgob/fFW2aPz/9axkDGzbA3ujgwlvf+i9UVS0C\nXFBy6lSmWy8i2aagQURiBgehp8d16n/8xz9g6dJLADc1UVEx/e2sV61ygYi1sHz5ZbHj4bDLkbB2\n6pUYIpI/FDSISIxXo6GuLsySJTfEjvf2wtKlUFk5veutWAFVVS4Q8ZedDoXcsYEBN4IhIoVBOQ0i\nEuMFDcn7TvT2utUQ07VkCTQ0jE96DIXc+0y1BbeI5BcFDSIS09XlaigEg4nHrXUjDdNVWurqOiQH\nBwoaRAqTggYRienqcjkG/mWVY2Pu63SSIP2WLo1fwxMMuuBEQYNIYVHQICIxXV1uVMHPWzkx3SRI\nT2NjPBnS4yU/qlaDSGFR0CAiMalKO/f1uQTImY40NDbGpyOSaaRBpLAoaBCRmJaW1EmQ8+ZNf+WE\np7HRjVT098OhQ7/m6afvBNxoQ0fHaTZYRGaVggYRiWlthcHBfTz33LcYHu4H3EjDqlUzv2ZdnVt2\n2dcHx449xxNPfA5rLaGQC1JEpHAoaBARIJ6Y2Nb2a3760w/HNqoaG4PFi2d+3UDALb3s74eGhjMZ\nGemnp+d4rCpkcg6FiOQvBQ0iAsRrNPT3N1NZuZDS0nBs1cNMkyA9y5a5Qk7x3S73xQo8pcp1EJH8\npKBBRIB40NDX10xt7QpgZntOpOK9vr5+NWBob9+rWg0iBUhBg4gA8aChpyceNPT2ugTI0x1paGx0\nX0tKwtTULEsIGrTsUqRwKGgQEQC6u11+QVdXPGjo63OjBBUVp3dt/wqKhoY1dHTsVYEnkQKkoEFE\nAK+wkyUSSRxpmMmeE8n8QUN9/Zm0te2JVZ1U0CBSOLTLpYgArrDT0FArIyMDsaDBWrfy4XRVVbml\nl62tsHTpxQwOxqtIKWgQKRwKGkQE8JY/tlNbu4La2pWnveeEnzFuBcXhw9DU9CGamj4UO64CTyKF\nQ0GDiABuFGD+/PX8xV8cAtzUxOnsOZFs6VIYHk48pgJPIoVFOQ0iwsiIu+P3l5D2Vk5kYqQB4iso\n/EIhF6yISGFQ0CAidHe75Y/JQcOCBae/csLT2Oi23B4ZiR8LheLvLSL5T0GDiMRqNPiDhr6+zKyc\n8DQ2ugCkry9+TLUaRAqLggYRiQUNwWD8mLWnt+dEMv+yS4+qQooUFgUNIkJXl9uYqjSaGp2pPSf8\nQiGXH+EfaQgGFTSIFBIFDSIyrtPu73dTCZlKgvSsWBEfaRgbG6WkxC27VNAgUhgUNIjIuJyC3l4X\nNGRypAFg4UI3irFjxz/xxS/WYO1YtHR1Zt9HRLJDQYOI0NoKnZ1P8/Wvr6er6wh9fS5gKC/P7Pt4\nyy4rKuYxPNxHb28LxrhqlCKS/6YdNBhjNhtjHjTGHDXGjBlj3pLinL81xhwzxvQZY35ujFmT9Hy9\nMeZ7xpiIMabDGHOvMabydH4REZm5U6dgcHA/bW2vUl7ewMBAZpMgPY2NUFYG5eWuTHUk0kwwqAJP\nIoViJiMNlcALwK2ATX7SGPNXwG3AR4CLgV7gEWOMLy+b7wNnA1cDbwauAO6ZQVtE5DSNjkJ7OwwM\nNFNRMY+ysgpGRjKfzwDxZZdlZcsBiEQOEwq5oEVE8t+0y0hbax8GHgYwxpgUp3wM+Ly19ifRc94L\nnATeBtxvjDkbuBZostY+Hz3no8DPjDGftNaemNFvIiIz0tMDAwPQ3x/f3RKgtjbz71Vf76pM9vc3\nUlpaTiTSTGNj6joRIpJ/MprTYIw5A1gEPO4ds9Z2AduAy6KHLgU6vIAh6jHcqMUlmWyPiEzN67B7\ne13QYKPjh3V1mX+vQMBNe/T3G2prVxCJNKtWg0gByXQi5CJc538y6fjJ6HPeOQkzmNbaUaDdd46I\nzJJIxHXaPT3N1NSsYGjI1U/IxkgDwPLlbmSjtnYFXV0KGkQKyWytnjCkyH+YwTkikmGRCIyMWLq6\n3EjDwICbJshW0ODlSmikQaTwZHpr7BO4zn8hiaMNC4DnfeckrP42xgSAesaPUCS4/fbbqU36JNuy\nZQtbtmw5vVaLzGGRCAwPdzI01BMLGsLh7ExPQHzZZVPTrQwPd1ESvXVR0CBy+rZu3crWrVsTjkUy\nuLlLRoMGa+0BY8wJ3KqIlwCMMTW4XIVvRE/7LVBnjNnoy2u4GhdsbJvs+nfddRebNm3KZJNF5rz2\ndjAmwBvf+FWWLr2Iri5XoyFbSYmNjS4oqa/fmLCDpoIGkdOX6kZ6x44dNDU1ZeT60w4aovUU1uA6\neYDVxpgLgHZr7WHgbuAzxpi9wEHg88AR4McA1trfG2MeAb5ljPkzIAh8DdiqlRMis6+lBSora9i4\n8aOxnzNdCdLPv9ulP2hQgSeR/DeTkYbXAL/E5R9Y4MvR498GPmCtvdMYU4Gru1AHPAm8yVo75LvG\nTcDXcasmxoAHcEs1RWSWtbS4O3/P0JAr95wt1dVQVZU4shAKqcCTSCGYSZ2GXzFFAqW19nPA5yZ5\nvhO4ebrvLSKZNTjochqSpyLq67P3nsa4oMQfJKjAk0hh0N4TInOYt9zSG2nwajRka+WEZ/Fi976e\nUMiNPAwNTfwaEck9BQ0ic1gkQmy1BMDIiCvAlK2VE56GhniAAi5oGBgYv9umiOQXBQ0ic1gkQqyY\nE8RHHbI90pB8/XBYQYNIIVDQIDKHeZ20t4uMV9gp2yMNdXVQUgJHjjzPjh33Egq54EUrKETym4IG\nkTnMFXaKsHv3Txgc7GJgwG0oVZnljerr6tzowt69v+Thhz+GVwxWIw0i+U1Bg8gcduoUdHe/xA9+\n8Ba6uo4yOOjKPKfcvzaDvKAhHF7B8HAf/f3tgEYaRPKdggaROezkSRgaagagtnY5AwOwaBa2jauq\ncqMZoZDbijsSaSYQgNbW7L+3iMycggaROWpszHXSQ0PNlJc3EAxWYW18b4hsMsZVnSwriwcN4TCc\nUE1YkbymoEFkjuruhv5+6O93u1t6sr1ywrN4McACAoFgbLfL1lYYHZ2d9xeR6VPQIDJHeYWdenub\nqalZztiYO57tlROehgaAEmpqltPVdTi27LK7e3beX0SmT0GDyBzlFXbq6TlEbe3KWJGn2Rpp8IKT\n2toVsemJgQElQ4rkMwUNInNUJAKjo5ZI5BB1dStnrbCTx6vVMG/euQQCZaoKKVIAZrLLpYgUgUgE\nRka6GR0dorZ2BQMDUF4+u0FDOAyve93XYnUhxsYUNIjkMwUNInNUezsEgzV8+tP9WDvGoUNuuWUg\nMDvv7wUNXkEpcKsqND0hkr80PSEyR5086facMKaEkpJSBgfdltWzpaoKKipc0OAxxgUzIpKfFDSI\nzFEtLW6fCc/IiKudMFtKStz7+YOGUMgFMyKSnxQ0iMxB3tJGb0tsz2zlM3gWLkwMGsJhF8z4t80W\nkfyhoEFkDvKWW3ojDV4nPVs1Gjzz5iUGCKEQ9PRAX9/stkNE0qOgQWQO8go7eSMNQ0Muv2G2RxqS\ngxTVahDJbwoaROagSCQeKEB81CEXQYMx8dLR4bALZrTsUiQ/KWgQmYO8TtnbAturBjnb0xPessun\nnrqLf/7nyyktdQmZChpE8pOCBpE5KBKBvr5m/uVfrqCl5RUGB90og381xWzwgobR0QDHjj0LuA0w\nND0hkp8UNIjMQS0t0N+/n+bmJwkEggwMzO5yS091tatCWV6+itHRQXp63HpLBQ0i+UlBg8gcdPIk\nDA83A1BTs5yhodkt7OQpKYH58yEYXAVAZ+dBgkE4cWL22yIiU1PQIDLHjI5CWxsMDByiomI+ZWXl\nANTX56Y9ixdDaekqADo7DxAKuZEQEck/ChpE5pjubpf42N/fTF3dylidhNleOeGZNw8CgRrKyxvo\n7DxIOOymJ4aGctMeEZmYggaROcYr7NTbe4ja2hWMjLhNqmZ75YTHe9+6ujPo6DgQW3apvAaR/KOg\nQWSO8YKGnp5mamtXxoo85WqkwavVUFu7ikjkIKGQa5+WXYrkH22NLTLHRCIwOmqJRJqprV2RsxoN\nHm/Z5XnnfYSSkj5CITc1oaBBJP8oaBCZY1xnPMyll97O8uV/wMAA1NRAZWVu2uMFDfX119DQED+u\n6QmR/KOgQWSO6eiAkpIgV1/9BQD27oW1a+PVIWdbdXV8zwk/jTSI5B/lNIjMMSdOJFZ+HBhwyx5z\nJRBwtRr8QUMgoGWXIvlIQYPIHGKt64y93S09/mmBXFi0KDFoCIddASoRyS8KGkTmkN5e6OpypZv9\ncpUE6Zk/H8bG4j+Hw64Albf7pYjkBwUNInNIRwf098dHGnJdo8GT/P6hkGtnV1du2iMiqSloEJlD\nOjri22BD7rbETua9vzfa4BV4UjKkSH5R0CAyh3R0uI65NLpuKp+ChnAYTpzYz86dD8QKPGnZpUh+\nUdAgMoe4oGGIo0e3MzTUy8AAVFXlrkaDp77e5Vns2fMoDzxwI8aMYK1GGkTyjYIGkTmkpQX6+w9w\n772XcOzYswwMuC2xc1WjwVNTAxUVEA6fgbWjdHUdARQ0iOQbBQ0ic8ixYzAy0gwQKyGdyxoNnpIS\nWLAAyspWAdDZeRBjoLU1t+0SkUQKGkTmiJEROHUKhoYOAYaamqVA7ms0eBYvhtLSlQCxLbKPHs1x\no0QkQcaDBmNMiTHm88aY/caYPmPMXmPMZ1Kc97fGmGPRc35ujFmT6baISFxnp1vGODDQTHX1EgKB\nIJD7JEjPvHlQUhKmqmoxHR0HqKhw0ykjI7lumYh4sjHS8NfAR4A/B84CPgV8yhhzm3eCMeavgNui\n510M9AKPGGOCWWiPiADt7W5FQl+f290yX2o0eOrr3de6OrdFdnm5K0bV0ZHbdolIXDaChsuAH1tr\nH7bWNltrfwQ8igsOPB8DPm+t/Ym19hXgvcAS4G1ZaI+I4DrfwUHo7j5EXd3KvFlu6amvd0FMbe0Z\ndHQcoLzcjYy0teW6ZSLiyUbQ8BvgamPMWgBjzAXAHwAPRX8+A1gEPO69wFrbBWzDBRwikgUdHW6V\nRFdXMzU1K2KVIfMlaPBqNVRWrmJoqIdgEIaH3QiJiOSHbGyN/UWgBvi9MWYUF5h82lr7g+jziwAL\nJG9HczL6nIhkQVsbWDtGJHI4tnKirs4tdcwH9fUuaGhq+p+86U1fiB3XSINI/shG0PAu4CbgRmAn\ncCHw98aYY9ba707yOoMLJiZ0++23U1tbm3Bsy5YtbNmy5fRaLDIHHD8OoZDhL/+yBWMCNDfnR40G\nT0WFq9fQ3h5vkDHa7VJkOrZu3crWrVsTjkUyWPAkG0HDncDfWWt/GP35d8aYVcAdwHeBE7gAYSGJ\now0LgOcnu/Bdd93Fpk2bMt1ekaJnLZw4AeXlhnDYzUcMDuZHjQaPMa49/mWWFRVw5Eju2iRSaFLd\nSO/YsYOmpqaMXD8bOQ0VjB8xGPPey1p7ABc4XO09aYypAS7B5UOISIb19kJ39/gtsfOlRoNn0SIY\nGor/XF7ull36j4lI7mRjpOEnwKeNMYeB3wGbgNuBe33n3A18xhizFzgIfB44Avw4C+0RmfO8LbGr\nq93PNhrW50sSpCc5iCkvd/Ul2ttdQCEiuZWNoOE2XBDwDdyUwzHgm9FjAFhr7zTGVAD3AHXAk8Cb\nrLW6nxDJAi9o8EYa8q1Gg8er1WCtm66oqHC5GAoaRPJDxoMGa20v8PHoY7LzPgd8LtPvLyLjtbe7\nLbEDAfdzvtVo8NTVQVmZm44IhdwW3iMjWkEhki+094TIHJBcVTFfgwZvi+yBAfezMW7UQUGDSH5Q\n0CAyB7S0uLt2z8CAy2/IlxoNHq/A0549j/HNb57PyMggpaVadimSLxQ0iMwBx47B4OBufvjDG4hE\nDjMwkF81GjzBoEuGHB4O0NLyMpHIIcrL4fDhXLdMREBBg0jRGx6G1lbo63uFnTt/SFlZOUND+ZtY\nuGQJBINnANDRsZ/ycjc94U1ZiEjuKGgQKXLelti9vXsIhWopL2/E2vyr0eBZuBCCweWUloZpbd1N\nRYVrv/agEMk9BQ0iRa6jw92l9/TspbFxLa4ga/4lQXrq6sCYAA0Na2lr2015OfT1KRlSJB8oaBAp\nch0dbgljZ+ceGhrW5G2NBo9Xq6Gx8SxaW39PaSmMjipoEMkHChpEipy33LKjYy8NDWtjyy29zjnf\neLtd1tW5oAFcwqamJ0RyT0GDSJFra4ORkV66u4/R0LAmb2s0eOrqXK2G6ur19PQcZ3Cwi9JStwJE\nRHJLQYNIkTt2DIaG9gHERhqqq8dvXpUvamtd2+rrX8s113wJay3l5W63S5u8FZ6IzCoFDSJFzFq3\nd0NFRZhqNKBGAAAgAElEQVQLL/wAjY3rGBhwyy3zrUaDp6TEraAoKzuD1772k4TDtVRUxFeBiEju\nKGgQKWI9Pe6xaNE63vrWf6KiojGvazR4Fi9OrMugFRQi+UFBg0gR85Zb+qci8rlGg2fevMSpiPJy\n1WoQyQcKGkSKmLcldjjsfvY64nxNgvQkr+wIBLRxlUg+UNAgUsQ6Olxn622JPTTk9neYNy+37ZpK\nfb1r88hI4nEFDSK5paBBpIgld7J9fW6ov7ExN+1Jl7fbpT+voaxMyy5Fck1Bg0gRO3TIjSx4+vuh\nsrIwpifC4cTVElp2KZJ7ChpEitToKDQ3uyDB09fndpH0pivyVUUF1NS4kYaDB5/g+PHnqaiASMSt\nBhGR3FDQIFKk2tqgqwtGR5vp6joCuE542bIcNywNxrhaDf398Mgjt/Pss/+oFRQieUBBg0iRammB\n3l54/vnP84MfvC12fP78HDZqGpYscYmb8+adRVvb72NBg5IhRXJHQYNIkTp50k1RdHS43S3Hxtwd\nfL6vnPB47WxsXE9r625KSlw+Q0tLbtslMpcpaBApUidOuK/t7Ym7W+b7ygnP/PkuyGloOIve3pMM\nDHQSCLg8DRHJDQUNIkXq4EEoK+uju/soDQ1rCma5pWf+fJcQWVm5HoDW1t1UVcHevVpBIZIrChpE\nitDgIBw9CiMjbnfLxsa19Pe7HSSrqnLcuDR5QUMotA6A1tbfU1XlEiGVDCmSGwoaRIrQqVMuCXJw\ncA9AbKRh2bL83d0yWSjkVlAMDVVSW7uCtjY30tDd7XbuFJHZp6BBpAh5Kyd6e/cSCtVQUTGf4WFY\nujTXLZueVavcion588+lv7+DYNCVlvbyNURkdpXmugEiknknT7qv3soJEx1eKJSVE55Fi1z+wk03\n/Sz2O5SUwOHDOW6YyByloEGkCB075qYhrrvua/T3tzMy4jrbQkmC9Hg1Jaw1sWmVysp4MmShTLWI\nFAtNT4gUoQMHXOdaWhqmunoJ/f0uqbAQg4aKClf+2lNV5XI2urpy1y6RuUpBg0iR6elxnWrynhOF\ntNzSM3+++z38QUN1tZIhRXJFQYNIkfGSIJODhvnz3YqEQlJe7trd2xs/Fgy68tJKhhSZfQoaRIpM\nS4tbcVBeHj/W3w/Ll+euTadj5crEkQZj3OPo0dy1SWSuUtAgUmS8lRP+JMGxMbcSoRAtXjy+AmRF\nhUuGFJHZpaBBpMgcOQJlZfGfvQ630PIZPN4KirGx+LGqKpfT4J+2EJHsU9AgUkSsja+c8AwNuSCi\nkIOG8nI4cOA57r57JZ2dB5UMKZIjChpEikhHB3R2wsDAK9x33xuJRA7T11eYyy093gqK0tLlRCLN\nHDmyjXAYBgYUNIjMNgUNIkXEWznR0bGN/ft/Tnl5A/39rtOtr89162amstJVsrR2AXV1qzh6dFss\nX0NBg8jsUtAgUkRaWtx0xKlTzzFv3lkEg5X09cGSJRAI5Lp1M7dqFdENty7lyJFnADdlsW9fbtsl\nMtcoaBApIidOuFUTx4/vYPHiTYAbxl+2LMcNO02LF7tEyKVLL+H48R2Mjg5RVeX2oBgYyHXrROYO\nBQ0iReTQISgrG+HkyRdZvLgpdtxbgVCovPYvWXIJo6ODnDz5UiwZUkWeRGaPggaRIjEy4u68h4d3\nMTIywOLFm2LLFAs1CdIzfz6Ew1BXt5GSkjKOHNkW25NCeQ0isycrQYMxZokx5rvGmFZjTJ8x5kVj\nzKakc/7WGHMs+vzPjTFrstEWkbmitdXdeXd37wBg0aKNDAy4zrYYgobKShgcDLNo0YUcPfpMLBlS\nIw0isyfjW2MbY+qAp4HHgWuBVmAt0OE756+A24D3AQeA/wk8Yow521o7lOk2icwFx465zara25+j\nsXEdoVA1ra2FvdzSU1kJDQ0uQLjmmi8RDtcBbi+N/ftz3DiROSTjQQPw10CztfaDvmOHks75GPB5\na+1PAIwx7wVOAm8D7s9Cm0SK3v79Lllw3bo3smTJRsDtOdHQ4HaGLGTGuBUU+/bBOedcGTteVQUH\nD8LwcGIVTBHJjmxMT/wR8Kwx5n5jzEljzA5jTCyAMMacASzCjUQAYK3tArYBl2WhPSJFz1p45RU3\nqrB27XVs3HgL4Go2LF+euA9FoVqyBEZHE4/V1LhiVocP56ZNInNNNoKG1cCfAbuBPwT+EfiqMebm\n6POLAIsbWfA7GX1ORKaprc3tOVFXl3h8ZATOOCM3bco0bwWFf/OqykoXGO3Zk5s2icw12ZieKAG2\nW2s/G/35RWPMubhA4r5JXmdwwcSEbr/9dmpraxOObdmyhS1btpxGc0UK38GDEIkk7mQ5MgIlJYW7\nJXYybwXFwEB8229jXF7Diy/Ctdfmtn0i+WDr1q1s3bo14VgkEsnY9bMRNBwHdiUd2wW8I/r9CVyA\nsJDE0YYFwPOTXfiuu+5i06ZNk50iMid5+Qylvn/RXV0ul6HQCzt5vBUUvb3xoAFckueePW60pdAT\nPkVOV6ob6R07dtDU1DTBK6YnG9MTTwPrk46tJ5oMaa09gAscrvaeNMbUAJcAv8lCe0SKmrXw8suJ\nO1uCCxoWLCiejrS62iV1Jm+HXV/vNup69dXctEtkLslG0HAXcKkx5g5jzJnGmJuADwJf951zN/AZ\nY8wfGWPOA74DHAF+nIX2iBS1U6fccsvkfIaeHjjrrOJIggT3e5x7rqtF4RcIuOd27sxNu0TmkoxP\nT1hrnzXGvB34IvBZXB2Gj1lrf+A7505jTAVwD1AHPAm8STUaRKbv4EE3qrBkSfyYte6xcmXOmpUV\nq1e7r2Nj8Pzz3wKgqelD1NXBCy8QK2YlItmRlYqQ1tqHrLXnW2srrLXnWmv/OcU5n7PWLomec621\ndm822iJS7Lx8hlOnXmTnzn/HWkt/v5v3L5YkSM/Kla42Q08PNDc/yY4dLnBobHQVMbXrpUh2ae8J\nkQLm5TNUVcHLL3+PRx/9BMYYurpcDYOlS3PdwsxavNglRHZ2uh0vT5x4gZGRAcJhtyX47t25bqFI\ncVPQIFLAWlrchk11dXDs2LOx7bC7ulx9hmIbqg8EYMMG9/stW3YJY2PDHD/uFl1VVsKOHcQ26RKR\nzFPQIFLADhxwHWh5eQ+HDz/NypVXADA4COvW5bhxWbJ6tRthWbDgfAKBEEePbgPcFMWRI3D0aI4b\nKFLEFDSIFDBvs6aDBx9ndHSIdeuuZ3TUFXUqlvoMyVatctMxfX1Blix5DYcO/Rpw0zHd3Vp6KZJN\nChpECtTYGLz0kutAX331pzQ2rqehYU2sqFOxJUF6liyJ5zWsX/9W9u79vwwOdmOM27Tq5Zdz3UKR\n4qWgQaRAnTzpHrW1lj17fsa6ddcDbrpi3jz3KEaBAJxzjhtVOPfcGxgZGWD37gcBV/xp1y5XUltE\nMk9Bg0iBOnDAdZx9fc/T03OctWvfDMSLOpUU8b/uNWvcSEtt7UquuOKz1Ne7Ag4NDdDerikKkWzJ\nxt4TIjILvHyGkZEeVq26ihUrLsda15muWpXTpmWdv17DVVf9bex4aalLkty5Ey66KIcNFClSRXwv\nIlK8RkZcPkN1NaxceQXve98vCATKYhURizUJ0rN0qVst0dk5/rmGBti2zRV7EpHMUtAgUoBeeQWa\nm11CoJ9X1KnYg4bSUrcPRVfX+OcWLXK1K37969lvl0ixU9AgUoCeftqNNlRUJB7v6nJD98nHi9GZ\nZ7qpGGsTj5eUuFGIxx93+Q0ikjkKGkQKzNGjrvLhokXjnxsYgPXJG9MXqVWrXBXI5K2ywZWbPn4c\nnnxy1pslUtQUNIgUmG3boKPD3U37FXtRp2TLlk2c1xAIQH09/PznqZ8XkZlR0CBSQPr63N1zXR0Y\nk/hcd7dbUTBXgobSUlevIVVeA7giUMeOabRBJJO05FKkgOzY4fZXSLWvxKlT7vjChbPfrlxZswYe\nftjlNRgDjz12B/X1q2lq+hCBgAuufv5zuOIKqK3NdWtFTt/ICPT3xx8DA+4xNOQeg4PxrwMD7mZi\n167Mvb+CBpECMTbmVgQEAlBaann++X9l7do3UVW1iLEx9wFy6aXFXdQp2RlnuNGV7m63aqSjYx/7\n9j1CU9OHADfasHOnSxy97rocN1bmPGtdpz8wEO/Ukx/e8cFB92+6uzv+6Olxz42MwPCw++o9xsbG\njz6C+zzI5PJjBQ0iBWLPHnfHsGQJtLXt5sEHP8CWLT9l3bo309Hh5vDPPz/XrZxdK1a46pcvvuiC\nhg0bbuT++99Ja+tu5s1bT2mpG2F49FHYvNnVtRCZKa/D99/h+zt7f4c/MOCSdHt73bSi9/3QUGJn\n7z1GR+OdvjEuwLDWTcOVlbmv3vehkEsC9h83JnXQALB7d+b+BgoaRArEb37jPnxqauDXv36A0tIw\nZ5xxFQAtLXDxxalXVBSzkhI39bBjh7vzWrPmTVRWLuBXv/ob3vnO7wOuENTOnfDjH8ONN7oPWJlb\nRkbinXrywxvKT3709bm7e6+z9zr85Dt87y4/mbXeqGD8q/cIhdyyaH+nX1IycaefT/TPR6QAtLXB\nM8+4Yk6Dg90888zdXHDB+ygrq4jdpVx8cWF86GTahRe60ZcTJ2D58nKuueZ/85//+V42bvxTVq++\nmtJSFzj8+Mfuw33LFvdhLYVldHT88L1319/f7zp57/uuLrdpWVdXfEh/dDR+R+99nzys793he9/7\nO3qv8/fu8v3HC6XDzwQFDSIF4Le/dYmOGzbAU0/9PUND3Wze/GnAzVfOmwfnnZfjRuZIZaUbbfje\n99zKkfPPv5nnn7+Xhx66lf/2316ktDREY6P7UH/wQddR3HwzBIO5bvnc493xe8P7qb56AUBPT7zT\n7+6Oz+X7h/O97yHeaVvrOvGysvHD+uFw4p1/IDC3OvxMUNAgkudefdXdJdfXw9BQJ7/97ZdpavoI\ntbXLATcKcc01bqXAXHXxxfDTn7r6FQ0Nhuuu+wb33LOR3/72y2ze/D8AtydFSQn87Geuw3nve91d\no0zOS96bajjfu/v3nvfm8nt64nP6g4Pxjt4/zO/d7Xt3+d5cvr/DT+74/Q91+rNHQYNIHotE4Nvf\ndgWKzjkHnnjiK4yMDHD55XcA7gO6pASamnLc0BxbtsxNUzz1lAsOFizYwKWX3s7Ro9ux1mKivYpX\n3+Lhh11n9fa3F3ceiLWuc06VpZ8qkc/L1u/rSxzuHx5OHN73f/U6fP+dPrj/L/1z+d73ZWVQXp4Y\nDOhuv3AoaBDJU6Oj8P3vuyS+c86BgYEOnnnmbi666FaqqxcDLgFy0SL3/FxmDFx+uUsW9Xb6fP3r\nv0BJSWksYPDU1roS1I8+Cs89Bxs3wmWXuamffJiyGBsbfzfv3ZX7v/qX7Hkdvpew593d9/W5a/mH\n8r3XJy/R82pd+Dt5/zC+19knJ/dNlrUvxUdBg0ieeuwx+OUvXS0Cd0dWxzvecR/Lll0aO6ejA66+\n2n2Yz3XnneeCgePH3d8sEJg427G62p3f1ub+xk8+6ZZvXnaZ27diwQL3qKxMv0O0NnVxneS7eX9A\n4BXf8S/L6+9PfTc/Ouq+eo9U7UqVra/kvbmjv7+d/v4ORkb6GR7uiz2OHn01Y++hoEEkD+3aBT/8\noVteWVPjjhljWL/+LbFz+vpcZ3DhhTlqZJ4JBuHKK+Hee12nOlWRq5IStxpl/nzXeR8/Dvfd554r\nL3ed7Lx5brrDP7ceCLiOP1USnz8z3393Pzrqrps8b5+coe/v9P139d5Qf0lJvLNXh18YrLWMjAxQ\nUhIgEJh4KKuz8yC7dv2IoaGe2GN4uDeh83/3u/8vpaXhCa/xyCO38+KL3xl3vLS0KiO/CyhoEMk7\nJ0/Cd77jMsfPPnvy81asgLVrZ69t+e6ii+A//9OtKFmwIP3XhcNudAJcZ+4N9R89CgcPumPeNtxe\ndr43fO+fu/fu6isqxgcD6uTzm9e5+zvsoSHXaQeDlSxdevGkr/+3f3sHvb0tsdf5r2PtGNdf/39i\nlUpTaW/fxxNPfI5gsCr6qKSsrDL6tYJwuI6xsdFJ2/Da1/4lF1zwfsrKyikrq6CsrILS0nIOHdrD\nj3501Yz+LskUNIjkibExV4vhhz+EQ4dcnsJEHc3YmBvW/uM/VrEiv/nz3UqKhx5y38+kozbGdfoV\nFe4akt8GBiIcP76D4eE+37B84vD85s13EAxOfLf98MO3s337V7E2RZUmYMWKzdxyy68nbUd5eQOh\nUHW0o3cdf/z7yoRpxVRWr76aO+6YYPe1NC1YsCHl8WDw5Gld108fNyJ5oKPDBQtPPOHuXM89192d\nTuTIEZcAedFFs9bEgnHNNS7B8dAhl+OQbHR0mJde+i4XXvh+jJlDG3XMAu9u3eusgdjS4Ik88cTf\n0NNzPOnu3rtT7+Xyy+9g06Y/nfD1p079ju985/UJx4wpSbjTvvji2yYNGs46623Mm7fed5df5bvL\nryQUqpnyd3/LW+6d8pxioKBBJIeGhty+CfffD3v3wsqVU9db8MrbbtlS3MsFZ2rlSrjhBrjnHrdU\nNfnveeDAL3jwwQ/y4ovfYdOmD3HWWW8jGKzMTWNzaGRkkEikOdrBu2F4bzje67QvvPAWwuGJtwd9\n+ukv8V//9Q3fNfoBG3t+8eImPvzhZydtx+HDT9Pf304w6O7Kw+E6amqWUlbmOux589ZP+vpFizby\n0Y/upaysnNJSNywfCATHrZqZzKpVV7Jq1ZVpnz+XKWgQmWU9Pa5g086d7o74+PH46II31TA42E0k\ncmjccKO1sH+/G2F43etmv+2F4sorXRD20EPxjX08a9Zcy003/ZSnnvoi//EfN1NWVslZZ72N88+/\nmZUrr6CsrCJ3DY8aGxulpGSSoSbglVf+jZ6e477O3t2lj4y4n889912cd96WCV/f2rqLe+7ZOOHz\nZWUVrF375kmDhgULNnDeee+Ozbu7jr4i9qioaJzyd33Pex6d8pzJlJWV09Bw5mldQ9KnoEEkS4aH\n3Z3uqVMuMe/UKTh82AUMra0uL6G21t0Z+5dMjowMcv/976C1dTcf/egeSkvjZQuPH3d3zjfcoP0T\nJlNSAu96l5ui2L17fH7I2rXXsXbtdXR2HuTll7/PSy99l5df/h5guPTSv+Daa78y4bVHR4dob9+H\ntaNYO8bY2ChjY8OMjAwyOjrI6OgQK1ZcPumQ9u9+90Neeum74zLlve/r68/k1lt3Tvo7PvXU39He\nvi/WQScmzlVOeafd0LCW973vl7GOPt7xu0c6Uzdr176JtWvfNOV5UjwUNIikaWxsfH18/9a3PT0u\nSGhtdev/vXr5fX0ugDDGdfQ1NW7FQ6pOv6vrKD/96Uc4dOhJbr754YSAYXDQXfc974ln+svEampc\nqeg773SrIJYtG39OXd0qNm/+H1x++R20tLzMsWPPUle3atLrRiLN/MM/TF5N68Mffo7FizdN+LyX\ncFdVtSjFPHoVVVULp/z9PvKRF6Y1BJ8sGKxk1arXzfj1MjcpaJA5YXQ0saDORGV0va9eCV1/0R2v\n8/dX1RsaSnwfY9ySu2DQfa2pcXkHZWWTZ/IPDfXw9NNf4je/+RLBYBU33PDv4z7Q9+1zBYmuvTbz\nf59itW4d/MmfwD/9kyvJXTvBSLsxhoULz2fhwvOnvGZ19VI+8IGnMaYEYwIYU0IgUEYgEKK0NEQg\nEKSiYvJlFxs2vIsNG941k18poc0is01Bg+Qtf0efvB2uv7Ke/6t3Z+9V1uvri5+fXHAneYc8G8/f\nwpjU1fVKSxPr5peVTV1EaDJjY6O8+OJ3+MUvPk1/fzuXXvoXXH75HePmkU+dcu/7J3+i6o/T9YY3\nwIED8PjjbvRn6dLTq5lQVlbO8uWvzVwDRbLAqykyNhb/nMsEBQ2SMZN18gMDrlNP3tbWq6znDfN7\nd/T9/eNr5vu/9+95D4k/e0V3kmvnl5bm3w55bW27efDBD7Bhw41cffX/Sjk03t7uCjm94x3aY2Im\nAgH44Afd8ssHHnAJqGvX5sc+E1KcvM46ufx3qu+n+pr8Wef/zEv+/PO/P8QriHadXvmHBAoa5iDv\nf0hveN3/NdWw/fBw/OEvn5tqeD/V3fzoaOL/4H7e/9ypOnjvWCg0/m7/dO7us8lay/BwL729LfT2\ntmBMyaSV5ObNO3vS+e/WVhcwvPWtbpRBI9IzU1oKb3wjrF4N3/0uvPKKq6ZZX5/rlsl0+Dtjr0P1\n/zzVw3++dyfuPbzr+7/6O23/Z1jyRl/JbTQm3mEnlwBP/j4USnyEw/Gv4XC8dLn/tamul6pCqffY\nvdtt0JYJChrykLWJnbq/Y/ce3nC8//nkjt2/ta13vjcn772HF0B43yfvb58qkk31P2ny3XyhdPJ+\no6NDDA52MzTUzeBgF4ODXSxdesmkGx9t2/ZVfve7+xka6qa/v4O+vlOMjAzEnl+58kre//4nJny9\nMWbCgKGlxSU+vvOdLmCYrNiTpGfdOvjUp1whrV/8whXJqq93lR9DoalfX+z8neVknfNEd9LJnbH/\nZ//rvK/JQbA3TZjqcye5M07umFM9AgE3ouRNJXoPf3nv5I7Xu66XuOyd7x+9TH6P5NFL/3v4j/mv\n530/G/+uBwczdy0FDachVefu78D93/uH672H16F7Gfn+YX3/znb+Tt2/GY7/H1zynLy/U0/+n9x7\npPqHUEid/JEjzzAw0JlQLnZ4uD9aSrafM8/8Q1as+IMJX3/s2HPcf/87o8vcuhkdHRp3zsc/fpTq\n6iUTXiMUqqGh4UyCwWrC4ToqKxdQUTGfysoFVFbOp7Z25Yx+t+PH3ZDiDTfA299eGP89CkV1Ndxy\niys3/dJLsG2bq30xMuISJcvL3cObypoNk90xe58zyY9UQ9n+TjrVXXGq903+eaqOOdWxsrJ45+x9\nfvi31E7uSEOhxM441Qhj8p1z8h4fyddP9fBepxG6zCnaoMHbk36iO/HkvelTfe8l1iV36v47/VT/\nmCebk/IPdaXq2JO/DwbHP+c98u0fQnzDl27KyxsoKZn4f6/9+x/jwIFfMjTUHeu0vbKxw8N9NDau\n5Z3v3Drp+z3wwLuIRJoTjrmKcK4yXE3N0kmDhsrKBZx33k0Eg9WEQtXRZW/VhEI1scdUWfAXXvh+\nLrzw/ZOeMx19fdDc7D7wbroJrr8+//47FwNjYMMG93jHO9zw7UsvuUdXl0s87e93/3491k6++VRy\nZ58c1HvX8H5O7rBTBfZeB+3/dx8Mxoeu/UPa3lC2/w52sg54qu+9a6TagdN/3DtP/5/ODQUVNGzf\n7oZs/UP03t168jz74OD44fdUUblnooQSL7Eu+R9VcvTrde4T/ePMN6Ojwylqvbuv9fWraWxcN+Fr\n29v38dOffjhWiCY+pN+NtW4Xtttue5XGxom3Xzx69L94+eX7fB22+1pRMZ+ysgrq66eu8Pa+9z1B\nIFAWLU5TTiAQmtYytNra5Vx99d+lfX42dXe7wk+BgEt2vPZadyesD+LsC4fhggvcw1oXuHV0xB+9\nvfEbh95eV4/D/9nh8ZbZenfd/s52osTcVHfKqTpq/538VMt3RbIp60GDMeYO4AvA3dbaj0ePhYCv\nAO8CQsAjwJ9ba1smu9b997u97SFxGC1VRO09kiPuVPvSF5JIpJn29n2MjAxEH/0JnXYoVM1rX/vJ\nSa/xzW+eT0vLyxM+f8UV/w9XXfU3Ez5fWhqisnIBdXWrY0Vp/B1/MFhNdfXiSduwefMdbN58x+S/\n7BTq6wu3wtHYWLwYVCTiOq7XvMYtDzz/fO1cmSvGuLLTlZWpi0GJzHVZ/WgyxlwEfAh4Mempu4E3\nAe8EuoBvAP8ObJ7seitWpN61Lh9YaxkbG5k0aa6/v53t27/B0FC3bw4+vlHM4GAX73jHfcybd9aE\n13jppfv4xS8+nXTUxDrr+fPPnjJouOKKz0T3iE/cyS2+YczkKeU1NcumnDoQx9p4noq/kqS1bm59\n3jzYvNmNKpx9dn6OSomIeLIWNBhjqoD7gA8Cn/UdrwE+ANxorf1V9NgtwC5jzMXW2u3ZahO48q3J\nw/FDQ71YOzbp/DfAI498gpaWl2OJdq6z9zLtu7nssk/wh3/4pQlfPzzcx7PP/kNCR+3VfK+qWkhD\nw1oCgclTuJuaPsKGDTdSWhqOVqALR+vEpz9kcu65N6R97lznJah501teEupE9SNSJakGg/HkuqVL\n4cwzXRnolStdtUitihCRQpHNkYZvAD+x1v7CGPNZ3/HXRN/3ce+AtXa3MaYZuAyYMGhob/89xsQ3\ndYk/ujnvvHdTVzdxpvqzz97Do49+PLbHe7LKyoV88pMnJv2FAoEywuE6ysoWU1rqku3iSXPVLFo0\n8Y5x4O7QP/GJ45OeM5WKisa0do6bK5I7dX/HntzB+/Na/Dkt/ixzb9rLn4HunwJLTgSrqIg/Kitd\n9n11tfu5vByqqtySvoYGt9GUph1EpJBl5SPMGHMjcCEuQEi2EBiy1ibXqDoJLJrsuj/5ybvHHSst\nDRMMVrFixeZJg4alSy/m9a//AsFg9bjheHe3Xz3Vr8Ub3vDFKc8Rx79MzOuw0ym4kmrdt//OPbmo\nirWpa0X4a0aEw/GO3fvZn7TmXw6WanmY972/AIv/9YWWFyMiMlMZDxqMMctwOQvXWGuHp/NSwE52\nwvXX38fKlU0Ju8FNlkPgt3jxRhYvnnwkYK5Ite47ucjTRJ23//tkyXfoyRnjE9WK8DLNg8F4Z+wt\nK/N/9bLH/V/9x9Wpi4hkVzZGGpqA+cBzJj7RHgCuMMbcBrwRCBljapJGGxbgRhsm9MIL/4dXX03c\nyGfDhi2cd96WjDV+tnmdbKqOO7lgy2QduPc1Vbnmiao6TlT/wSvWUlkZ73j9d+ehkBt69y8v85aC\n+b/3Om9/x55qzbc6dRGRzNi6dStbtyYmqkcikYxd39jkwgSne0FjKoHkeYJ/BXYBXwSOAqdwiZD/\nEbnC79YAAAh8SURBVH3NOuD3wKWpEiGNMZuA5973vudYtWriPeozaaKOPJ078OTn4r9H/Nr+9/HX\nfZhsaah3Z+6/m04u8BIKTVwdLbmUaqq7dX+ZU3XmIiKFb8eOHTQ1NQE0WWt3nM61Mj7SYK3tBXb6\njxljeoE2a+2u6M//BHzFGNMBdANfBZ6eauWEt4ti8pz4ZPPlEz28zjw5Cc79DqmLOqUq7uS/I/cP\niyd36MnD6sk1yJMLuPjvzP1D/OrIRUQkV2Yrlzt5OON2YBR4AFfc6WHg1qkucvKkS6pLtZdCcnXG\niR5eh1xentgxp9rUJLmDnyhZTh25iIjMBbMSNFhrX5/08yDw0egjbR/+sCv1OtFuYv7NUtSRi4iI\nZFZBrRpfvRrWTrydgYiIiGSRitaKiIhIWhQ0iIiISFoUNIiIiEhaFDSIiIhIWhQ0iIiISFoUNIiI\niEhaFDSIiIhIWhQ0iIiISFoUNIiIiEhaFDSIiIhIWhQ0iIiISFoUNIiIiEhaFDSIiIhIWhQ0iIiI\nSFoUNIiIiEhaFDSIiIhIWhQ0iIiISFoUNIiIiEhaFDSIiIhIWhQ0iIiISFoUNIiIiEhaFDSIiIhI\nWhQ0iIiISFoUNIiIiEhaFDSIiIhIWhQ0iIiISFoUNIiIiEhaFDSIiIhIWhQ0iIiISFoUNIiIiEha\nFDSIiIhIWhQ0iIiISFoUNIiIiEhaFDSIiIhIWhQ0iIiISFoUNIiIiEhaFDSIiIhIWhQ0iIiISFoU\nNIiIiEhaFDTIhLZu3ZrrJsw5+pvPPv3NZ5/+5oUr40GDMeYOY8x2Y0yXMeakMeY/jDHrks4JGWO+\nYYxpNcZ0G2MeMMYsyHRb5PToH/bs09989ulvPvv0Ny9c2Rhp2Ax8DbgEeANQBjxqjCn3nXM38Gbg\nncAVwBLg37PQFhEREcmQ0kxf0Fp7nf9nY8z7gRagCXjKGFMDfAC40Vr7q+g5twC7jDEXW2u3Z7pN\nIiIicvpmI6ehDrBAe/TnJlyw8rh3grV2N9AMXDYL7REREZEZyPhIg58xxuCmIp6y1u6MHl4EDFlr\nu5JOPxl9LpUwwK5du7LSTkktEomwY8eOXDdjTtHffPbpbz779DefXb6+M3y61zLW2tO9xsQXN+ab\nwLXA5dbaY9FjW4B/ttaWJ527HXjMWvs/UlznJuB7WWuoiIhI8Xu3tfb7p3OBrI00GGO+DlwHbPYC\nhqgTQNAYU5M02rAAN9qQyiPAu4GDwEAWmisiIlKswsAqXF96WrIy0hANGN4KXGmt3Z/0XA1wCpcI\n+R/RY+uA3wOXKhFSREQkP2U8aDDG/AOwBXgL8KrvqYi1dsB3zpuAW4Bu4KvAmLV2c0YbIyIiIhmT\njaBhDLdaItkt1trvRM8JAf8bF1yEgIeBW621LRltjIiIiGRMVhMhRUREpHho7wkRERFJi4IGERER\nSUveBw3GmFuNMQeMMf3GmGeMMRfluk3FLJ0NxyR7on//MWPMV3LdlmJmjFlijPludNO8PmPMi8aY\nTbluV7EyxpQYYz5vjNkf/XvvNcZ8JtftKjbGmM3GmAeNMUejnyNvSXHO3xpjjkX/O/zcGLNmOu+R\n10GDMeZdwJeB/xfYCLwIPGKMmZfThhW3dDYckyyIBsQfwv1/LllijKkDngYGccXnzgY+AXTksl1F\n7q+BjwB/DpwFfAr4lDHmtpy2qvhUAi8At5JiQYIx5q+A23D/LS4GenF9ajDdN8jrREhjzDPANmvt\nx6I/G+Aw8FVr7Z05bdwcEQ3QWoArrLVP5bo9xcoYUwU8B/wZ8FngeWvtx3PbquJkjPkicJm19spc\nt2WuMMb8BDhhrf2Q79gDQJ+19r25a1nxiq5kfJu19kHfsWPAl6y1d0V/rsEVVXyftfb+dK6btyMN\nxpgy3OZW/o2tLPAY2thqNiVvOCbZ8Q3gJ9baX+S6IXPAHwHPGmPuj07B7TDGfDDXjSpyvwGuNsas\nBTDGXAD8AfBQTls1hxhjzsDt7+TvU7uAbUyjT83qhlWnaR4QYHxp6ZPA+tlvztwzwYZjkmHGmBuB\nC4HX5Lotc8Rq3IjOl4Ev4KbivmqMGbDW3pfTlhWvLwI1wO+NMaO4G9ZPW2t/kNtmzSmLcDeAqfrU\niTaLHCefg4aJGFIXj5LM+wfgHNwdgWSBMWYZLjC7xlo7nOv2zBElwHZr7WejP79ojDkXF0goaMiO\ndwE3ATcCO3FB8t8bY45Za7+b05bJtPrUvJ2eAFqBUWBh0vHJNraSDPFtOPY6a+3xXLeniDUB84Hn\njDHDxphh4ErgY8aYoehoj2TWcWBX0rFdwIoctGWuuBP4X9baH1prf2et/R5wF3BHjts1l5zABQin\n1afmbdAQvet6DrjaOxb9AL0aNz8mWeLbcOwqa21zrttT5B4DzsPdeV0QfTyLu+O9wOZzpnLheprx\nU5zrgUM5aMtcUcH4u9kx8rgPKjbW2gO4wMHfp9bgpufS7lPzfXriK8C3jTHPAduB23H/8/1rLhtV\nzJI2HOs1xnhRaWzDMckca20vbrg2xhjTC7RZa5PvhiUz7gKeNsbcAdyP+9D8IG65q2THT4BPG2MO\nA78DNuE+z+/NaauKjDGmEliDG1EAWB1NOm231h7GTYV+xhizFzgIfB44Avw47ffI9xsZY8yf49b0\nLsStP/2otfbZ3LaqeKWz4ZhklzHmF8ALWnKZPcaY63DJeWuAA8CXrbX/nNtWFa9oZ/Z54O244fBj\nwPeBz1trR3LZtmJijLkS+CXjP8O/ba39QPSczwEfxq2MexK3WeTetN8j34MGERERyQ+aTxIREZG0\nKGgQERGRtChoEBERkbQoaBAREZG0KGgQERGRtChoEBERkbQoaBAREZG0KGgQERGRtChoEBERkbQo\naBAREZG0KGgQERGRtPz/uODyThSSpbkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe1bfa05f50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "go(fitfun, ar2, model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "4px",
    "width": "254px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  },
  "widgets": {
   "state": {
    "c93e1f0200634f2ab013baac1d1e439a": {
     "views": [
      {
       "cell_index": 14
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}