{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PEtab import and yaml2sbml\n",
    "\n",
    "[PEtab](https://petab.readthedocs.io/en/stable/) is a format for specifying parameter estimation problems in systems biology. This notebook illustrates how the PEtab format can be used together with the ODE simulation toolbox [AMICI](https://amici.readthedocs.io/en/latest/) to define ODE based parameter estimation problems for pyABC. Then, in pyABC we can perform exact sampling based on the algorithms introduced in [this preprint](https://www.biorxiv.org/content/10.1101/2020.01.30.927004v1.abstract).\n",
    "\n",
    "To use this functionality, you need to have (at least) PEtab and AMICI installed. Further, this notebook uses [yaml2sbml](https://yaml2sbml.readthedocs.io/en/latest/). You can install all via:\n",
    "\n",
    "    pip install pyabc[petab,amici,yaml2sbml]\n",
    "\n",
    "AMICI may require some [external dependencies](https://github.com/ICB-DCM/AMICI/blob/master/INSTALL.md)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# install if not done yet\n",
    "!pip install pyabc --quiet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pyabc\n",
    "import pyabc.petab\n",
    "import amici.petab_import\n",
    "\n",
    "%matplotlib inline\n",
    "pyabc.settings.set_figure_params('pyabc')  # for beautified plots\n",
    "\n",
    "# folders\n",
    "dir_in = 'models/'\n",
    "dir_out = 'out/'\n",
    "os.makedirs(dir_out, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate PEtab problem from YAML\n",
    "\n",
    "In this section, we use the tool [yaml2sbml](https://yaml2sbml.readthedocs.io/) to generate a PEtab problem based on a simple human-editable YAML file stored under `dir_in`, combining it with \"measurement\" data as generated in a later section. `yaml2sbml` is a simple way of manually generating models. The PEtab import below works independent of it with any valid PEtab model. We use the common conversion reaction toy model, inferring one parameter and the noise standard variation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "YAML file is valid ✅\n"
     ]
    }
   ],
   "source": [
    "import shutil\n",
    "import yaml2sbml\n",
    "import petab\n",
    "\n",
    "#  check yaml file\n",
    "model_name = 'cr'\n",
    "yaml_file = dir_in + model_name + '.yml'\n",
    "yaml2sbml.validate_yaml(yaml_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The format allows to compactly define ODEs, parameters, observables and conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "odes:\n",
      "  - stateId: x1\n",
      "    rightHandSide: (- theta1 * x1 + theta2 * x2)\n",
      "    initialValue: 1\n",
      "\n",
      "  - stateId: x2\n",
      "    rightHandSide: (theta1 * x1 - theta2 * x2)\n",
      "    initialValue: 0\n",
      "\n",
      "parameters:\n",
      "  - parameterId: theta1\n",
      "    parameterName: $\\theta_1$\n",
      "    nominalValue: 0.08\n",
      "    parameterScale: lin\n",
      "    lowerBound: 0.05\n",
      "    upperBound: 0.12\n",
      "    estimate: 1\n",
      "\n",
      "  - parameterId: theta2\n",
      "    parameterName: $\\theta_2$\n",
      "    nominalValue: 0.12\n",
      "    parameterScale: lin\n",
      "    lowerBound: 0.05\n",
      "    upperBound: 0.2\n",
      "    estimate: 0\n",
      "    \n",
      "  - parameterId: sigma\n",
      "    parameterName: $\\sigma$\n",
      "    nominalValue: 0.02\n",
      "    parameterScale: log10\n",
      "    lowerBound: 0.002\n",
      "    upperBound: 1\n",
      "    estimate: 1\n",
      "\n",
      "observables:\n",
      "  - observableId: obs_x2\n",
      "    observableFormula: x2\n",
      "    observableTransformation: lin\n",
      "    noiseFormula: noiseParameter1_obs_x2\n",
      "    noiseDistribution: normal\n",
      "\n",
      "conditions:\n",
      "  - conditionId: condition1\n"
     ]
    }
   ],
   "source": [
    "with open(yaml_file, 'r') as f:\n",
    "    print(f.read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We combine the YAML model with \"measurement\" data (which can be generated as in a later section) to create a PEtab problem. This generates a stand-alone PEtab folder under `dir_out`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32mChecking SBML model...\n",
      "\u001b[0m\u001b[32mChecking measurement table...\n",
      "\u001b[0m\u001b[32mChecking condition table...\n",
      "\u001b[0m\u001b[32mChecking observable table...\n",
      "\u001b[0m\u001b[32mChecking parameter table...\n",
      "\u001b[0m\u001b[32mOK\n",
      "\u001b[0m\u001b[0m"
     ]
    }
   ],
   "source": [
    "# convert to petab\n",
    "petab_dir = dir_out + model_name + '_petab/'\n",
    "measurement_file = model_name + '_measurement_table.tsv'\n",
    "yaml2sbml.yaml2petab(\n",
    "    yaml_file, output_dir=petab_dir, sbml_name=model_name,\n",
    "    petab_yaml_name='cr_petab.yml',\n",
    "    measurement_table_name=measurement_file)\n",
    "\n",
    "# copy measurement table over\n",
    "_ = shutil.copyfile(\n",
    "    dir_in + measurement_file, petab_dir + measurement_file)\n",
    "\n",
    "petab_yaml_file = petab_dir + 'cr_petab.yml'\n",
    "# check petab files\n",
    "!petablint -v -y $petab_yaml_file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import PEtab problem to AMICI and pyABC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We read the PEtab problem, create an AMICI model and then import the full problem in pyABC:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read the petab problem from yaml\n",
    "petab_problem = petab.Problem.from_yaml(petab_yaml_file)\n",
    "\n",
    "# compile the petab problem to an AMICI ODE model\n",
    "amici_dir = dir_out + model_name + '_amici'\n",
    "if amici_dir not in sys.path:\n",
    "    sys.path.insert(0, os.path.abspath(amici_dir))\n",
    "model = amici.petab_import.import_petab_problem(\n",
    "    petab_problem, model_output_dir=amici_dir, verbose=False)\n",
    "\n",
    "# the solver to numerically solve the ODE\n",
    "solver = model.getSolver()\n",
    "\n",
    "# import everything to pyABC\n",
    "importer = pyabc.petab.AmiciPetabImporter(petab_problem, model, solver)\n",
    "\n",
    "# extract what we need from the importer\n",
    "prior = importer.create_prior()\n",
    "model = importer.create_model()\n",
    "kernel = importer.create_kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once everything has been compiled and imported, we can simply call the model. By default, this only returns the log likelihood value. If also simulated data are to be returned (and stored in the pyABC datastore), pass `return_simulations=True` to the importer. `return_rdatas` returns the full AMICI data objects including states, observables, and debugging information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'llh': 22.37843729780134}\n"
     ]
    }
   ],
   "source": [
    "print(model(importer.get_nominal_parameters()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can inspect the prior to see what parameters we infer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<Distribution\n",
      "    theta1=<RV name=uniform, args=(), kwargs={'loc': 0.05, 'scale': 0.06999999999999999}>,\n",
      "    sigma=<RV name=uniform, args=(), kwargs={'loc': -2.6989700043360187, 'scale': 2.6989700043360187}>>\n"
     ]
    }
   ],
   "source": [
    "print(prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run analysis with exact inference\n",
    "\n",
    "Now we can run an analysis using pyABC's exact sequential sampler under the assumption of measurement noise. For details on the method check the [noise assessment notebook](noise.ipynb). If instead standard distance-based ABC is to be used, the distance function must currently be manually defined from the model output. Here we use a population size of 100, usually a far large size would be preferable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:Sampler:Parallelizing the sampling on 4 cores.\n",
      "INFO:History:Start <ABCSMC id=1, start_time=2021-02-16 13:09:39.040226>\n",
      "INFO:ABC:Calibration sample before t=0.\n",
      "INFO:ABC:t: 0, eps: 17.221164423753635.\n",
      "INFO:ABC:Acceptance rate: 100 / 390 = 2.5641e-01, ESS=9.9976e+01.\n",
      "INFO:ABC:t: 1, eps: 8.610582211876817.\n",
      "INFO:ABC:Acceptance rate: 100 / 373 = 2.6810e-01, ESS=9.0797e+01.\n",
      "INFO:ABC:t: 2, eps: 4.305291105938409.\n",
      "INFO:ABC:Acceptance rate: 100 / 474 = 2.1097e-01, ESS=8.0792e+01.\n",
      "INFO:ABC:t: 3, eps: 2.1526455529692043.\n",
      "INFO:ABC:Acceptance rate: 100 / 668 = 1.4970e-01, ESS=9.2427e+01.\n",
      "INFO:ABC:t: 4, eps: 1.0763227764846022.\n",
      "INFO:ABC:Acceptance rate: 100 / 506 = 1.9763e-01, ESS=9.4616e+01.\n",
      "INFO:ABC:t: 5, eps: 1.0.\n",
      "INFO:ABC:Acceptance rate: 100 / 208 = 4.8077e-01, ESS=6.3028e+01.\n",
      "INFO:pyabc.util:Stopping: minimum epsilon.\n",
      "INFO:History:Done <ABCSMC id=1, duration=0:01:26.428840, end_time=2021-02-16 13:11:05.469066>\n"
     ]
    }
   ],
   "source": [
    "sampler = pyabc.MulticoreEvalParallelSampler()\n",
    "\n",
    "temperature = pyabc.Temperature()\n",
    "acceptor = pyabc.StochasticAcceptor()\n",
    "\n",
    "abc = pyabc.ABCSMC(model, prior, kernel,\n",
    "                   eps=temperature,\n",
    "                   acceptor=acceptor,\n",
    "                   sampler=sampler,\n",
    "                   population_size=100)\n",
    "# AMICI knows the data, thus we don't pass them here\n",
    "abc.new(pyabc.create_sqlite_db_id(), {})\n",
    "h = abc.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's it! Now we can visualize our results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[<AxesSubplot:ylabel='log10($\\\\sigma$)'>, <AxesSubplot:>],\n",
       "       [<AxesSubplot:xlabel='log10($\\\\sigma$)', ylabel='$\\\\theta_1$'>,\n",
       "        <AxesSubplot:xlabel='$\\\\theta_1$'>]], dtype=object)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVUAAAFgCAYAAAALu+owAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABG8UlEQVR4nO3de5xU5ZXw+9/q6hvdTTfd3K8CKioiEm0FE8cbalBjiBliNIyXxHnJfCa+553byeV43txmPG9yciaZ5J3cmGgUQ0TGyGCUUZFoHI2iaBBBRJGAcocG+n6r7nX+2E81RXVV36r2rtrV6/v57E9X7b1r16qmWP3svZ9nPaKqGGOMyYyCbAdgjDH5xJKqMcZkkCVVY4zJIEuqxhiTQZZUjTEmgyypGmNMBllSNSYJEVkkIjtEZKeIfDXJ9hIRecRt3ygi07MQpslBllSNSSAiEeDHwHXAbOBWEZmdsNtdwHFVPQP4AfDdYKM0ucqSqjG9XQzsVNVdqtoBrAIWJ+yzGHjQPX4UWCgiEmCMJkflbVJdtGiRArbkxxK0ycCHcc/3unVJ91HVKFAPjE48kIgsE5FNIrLp3HPPzfbv0ZYAvpN5m1SPHj2a7RCMQVWXq2qtqtaOGDEi2+GYAORtUg2L5557jueeey7bYZhT7QOmxj2f4tYl3UdECoEqoC6Q6ExOs6SaZQ0NDTQ0NGQ7DHOq14AzRWSGiBQDtwCPJ+zzOHCHe7wE+J1adSIDFGY7gOFu8eLF2Q7BJFDVqIjcDTwNRID7VXWbiHwb2KSqjwP3AQ+JyE7gGF7iNcaSaq441txBcWEBFSX2T5ILVHUdsC5h3dfjHrcBnwk6LpP77H9wlj377LMA/PO2EmZPrOSfbz4/yxEZY9JhSTXLWltb6Yh2s/1AOxG7wm1M6Nl/4yy78cYbGXvuRwH4oK4ly9EYY9JlSTUHvPnhCQAa2qLUt3RmNxhjTFosqWbZM888w+4tL/c8/+CYtVaNCTNLqlnW2dlJXWMr50ysBCypGhN2llSzbN7HruL5lil8Yu5EwJKqMWFnSTXLYtdTLz1jDDXlxXx43JKqMWFmXaqy7PWXnuOS4gbOnjiSqTVlfGgtVWNCzVqqWXaksZ2a8mJKCiNMrR5hp//GhJwl1Szq6laebZzImLMvBmBaTRn7jrcS7erOcmTGmKGypJpFfzraTEtHF+dNGQV4STXarRyob8tuYMaYIbOkmkXHWzpYULSHEzs2Al5SBey6qjEhZkk1ixrbOunSAkpLigGYGkuq1gPAmNCypJpFjW1RXotO5ZLLrgRgYlUphQViN6uMCTFLqlnU2BYFYKSroVoYKWDSqBF8cKw1m2EZY9JgSTWLGtuifLRoNy8990zPumk1ZdZSzSIRqRGR9SLynvtZnWK/LhHZ7JbEqVbMMGZJNYsa2zrppJCR5WU962wAQNZ9FdigqmcCG9zzZFpVdZ5bPhlceCbXWVLNoqb2KO8VzeDqq6/uWTelegTHmjto6+zKYmTD2mLgQff4QeBT2QvFhFFokqqIlIrIqyLypohsE5FvZTumdDW2RRlZeupI4dHlXk+AuuaObIRkYLyqHnCPDwLjU+xXKiKbROQVEflUqoOJyDK336YjR45kOlaTg0KTVIF24CpVPR+YBywSkQXZDSk9jW2dnN/1HmvXru1ZV+2S6nFLqr4RkWdFZGuS5ZSpbd2U06mmnT5NVWuBzwH/IiKnJ9tJVZeraq2q1o4dOzazH8TkpNAUVHFf8Cb3tMgtoZ5nvbEtSklRGZWVlT3rrKXqP1W9OtU2ETkkIhNV9YCITAQOpzjGPvdzl4g8D3wEeN+PeE24hKmliohERGQz3hd9vapuTNgeqlOtxrYo7WPP5sorr+xZZy3VrHscuMM9vgNYm7iDiFSLSIl7PAb4GPB2YBGanBaqpKqqXao6D5gCXCwicxK2h+pUq7G9k5GlRaess5Zq1n0HuEZE3gOuds8RkVoR+YXb5xxgk4i8CTwHfEdVLakaIESn//FU9YSIPAcsArZmO56hamqLUnnoDR57bBef/vSnAagsLSJSINZSzRJVrQMWJlm/CfhL9/gPwHkBh2ZCIjQtVREZKyKj3OMRwDXAO1kNKg2qSmNblOLyKkaPHt2zvqBAqC4rspaqMSEVppbqROBBEYng/TFYrapPZDmmIWvr7CbardScPo/LLz/1xnF1WbG1VI0JqdAkVVXdgneHNS80tncC9OqnClBTXswxS6rGhFJoTv/zTayYSt3WF3j00UdP2VZTXsyxFkuqxoRRaFqq+SaWVEeNHsuEsRWnbLOWqjHhZUk1SxrbvNP/2R+5mIum15yyraa8mBMtHXR1K5ECyUZ4xpghstP/LGlyLdWKkuTXVLsV6ls7gw7LGJMmS6pZEjv9f3XDE6xevfqUbTVuAIBdAjAmfOz0P0sa3On/1KlTKS2KnLLNkqox4WVJNUua2r2W6uV/dmmv66bVZZZUjQkrO/3Pksa2KOXFkaQ3okZXWFI1JqwsqWZJY5tXTOXhhx/m4YcfPmVbrKV63PqqGhM6dvqfJU3tUSpKC5kxY0avbaVFEcqLI9Q1WVI1JmwsqWZJbCqVBQuST15QXV5sLVVjQshO/7OkoS3aq5ZqvNHlxVapypgQsqSaJU1tnYwsLWTlypWsXLmy1/bqcqtUZUwYWVLNksa2KCNLCpk1axazZs3qtd3G/2eHiHzGzdbbLSK1fey3SER2iMhOEflqkDGa3Bb4NVURKQfaVHVYT2wfu6Z60UVzk26vKbOkmiVbgU8DP0+1g6vp+2O8Qul7gddE5HGbUsVAAC1VESkQkc+JyJMichivWv8BEXlbRL4nImf4HUOu6ezqprWzq89rqjUVxbR2dtHaMaz/9gROVber6o5+drsY2Kmqu1S1A1gFLO7nNWaYCOL0/zngdOBrwARVnaqq44BLgVeA74rIXwQQR85obj9ZTGXFihWsWLGi1z41sVFV1gMgF00GPox7vtet6yVsM/ya9AVx+n+1qvYqt6Sqx4DfAL8RkdRNtjwUK6YysrSQmeeem3SfnvH/TR1MHjUisNiGAxF5FpiQZNM9qtprSup0qOpyYDlAbW2tZvLYJjf5nlRjCVVEvqmq3+xrn+EiVkxlZGkRF865MOk+PUnVWqoZp6pXp3mIfcDUuOdT3DpjAr1R9XU3C2oN8AawSlWPB/j+OaMprqWayslKVe2BxGQG5TXgTBGZgZdMbwE+l92QTK4IskuVAm3A03h/5f8gIucH+P45I/70/4EHHuCBBx7otU8sqR5vHlaN+KwTkZtEZC9wCfCkiDzt1k8SkXUAqhoF7sb7Lm/Hm9l3W7ZiNrklyJbqO6r6Dff4URF5APgZcFWAMeSEkzOpFjFv3ryk+1SWFlEgVlQlaKq6BliTZP1+4Pq45+uAdQGGZkIiyKR6VEQuVNXXAVT1XREZG+D754z4qVRmpEiqBQVCtfVVNSZ0gkyq/wewSkReB94C5gJ/CvD9c0ZD3Ol/V5fXDzUSifTaz4qqGBM+gV1TVdU3gXlArHjoc8CtQb1/Lmlsi1IcKaC0KMJDDz3EQw89lHS/6rIia6kaEzK+t1RFRFRVAVS1HXjSLUn36eM4U4EVwHi8m17LVfWH/kTtr0ZXTAXgggsuSLlfdVkxe+paggrLGJMBgYyoEpH/LiLT4leKSLGIXCUiDwJ3DOA4UeDvVXU2sAD4kojM9iFe38UKVAPMnTuXuXNTjP8vL7Z+qsaETBDXVBcBXwAeFpGZwHFgBF5Cfwb4F1X9Y38HUdUDwAH3uFFEtuMNDQxdEYtYMRWAzk6vJ0BRUe9BZbHyf6qKSO+5rIwxuSeIEVVtwE+An7jhqGOAVlU9MdRjish04CPAxoT1y4BlANOmTev9whzR2NbJyBIvicZqqd5555299qspKybarTS1913Q2hiTOwK5+y8iZ+NV8YkVndgnImtV9Z0hHKsCr2bA36hqQ/y2sIyzbmyLMq2mDIDa2pQlO6mOGwBgSdWYcAii9N9X8EqjCfCqWwSve9Wgivu6lu5vgJWq+limYw1KY9xUKnPmzGHOnDlJ96sp9/ax66rGhEcQLdW7gHMTi6aIyPeBbcB3BnIQ8S4q3gdsV9XvZzzKADXE3f1va2sDoLS0tNd+o2JTVVu3KmNCI4i7/93ApCTrJ7ptA/Ux4DbgKhHZ7Jbr+3tRrlGNXSP1kuqqVatYtWpV0n17aqpaUjUmNIJoqf4NsEFE3uNkYd9pwBnAfx/oQVT1RbzLBqHW3NGF6skKVfPnz0+5b881VTv9NyY0fG+pqupTwCzgW3hVfZ4Gvgmcpar/6ff755rGuFqqAOeccw7nnHNO0n0rSwuJFIi1VI1vVq5cyfTp0ykoKGD69OlJZ/Y1gxPI3X9V7cabOuUUIvJ5Vf1lEDHkisaEWqotLd6IqbKysl77inhFVaylavywcuVKli1b1vMd3LNnD8uWLQNg6dKl2Qwt1LI9RfW3svz+gYu1VCtKvKS6evVqVq9enXL/mnIb/2/8cc899/Qk1JiWlhbuueeeLEWUH4LoUrUlxfIW3jj+YeVkS9U7/b/kkku45JJLUu7vtVStULXJvA8++GDA6+0ywcAFcfo/Hvg43vDUeAL8IYD3zymxpFrpTv/POuusPvevKS9m5+Em3+MyHhH5DN41/3OAi1V1U4r9dgONQBcQVdXUozhy1LRp09izZ0/S9fHsMsHgBHH6/wRQoap7EpbdwPMBvH9OSWypNjU10dSUOmmOsmuqQdsKfBp4YQD7Xqmq88KQUONbmmPGjGHMmDFJEyp4STO+NWqXCQYniLH/d/WxbdhNlnby7r/3q3/00UeB5GP/wbumerylk+5upaAg9D3Kcp6qbgfyqoBNYkuzrq6u39fEt0YHc5nABFv53+C1VAsEyoq9Sv+XXnppn/tXlxXT1a00tkWpKrPx/zlEgWdERIGfu7oTveRCkZ9kLc2BaGlp4S/+4i9Sbs/lokXZFFhSFZG/S7K6HnhdVTcHFUe2NbZ1UlFS2NMSOuOMM/rcv2eq6pYOS6oZIiLPAhOSbLpHVdcO8DCXquo+ERkHrBeRd1S11yWDXCjy40eLsqysjHvvvTfjx80HQbZUa93yW/f8E8AW4K9E5N9V9f8NMJasaUwo41dfXw9AVVVV0v1jo6qONXcwY0y5/wEOA6p6dQaOsc/9PCwia4CLGdh12MCluiE1VKeddhr33nuv3aRKIch+qlOAC1T171X174ELgXHAZcCdAcaRVfEFqgHWrFnDmjW9ZkTuUWNFVXKOiJSLyMjYY+BavBtcOenee+9NOrhkqHbv3m0JtQ9BtlTHAe1xzzuB8araKiLtKV6Td+LnpwK47LLL+ty/xsb/B0pEbgL+NzAWeFJENqvqx0VkEvALVb0er5vgGncJpxD4tRuOnZNiCfCOO+7omb3X+CfIpLoS2Cgia/H6qH4C+LX7Sx+6KVGGqrEtyvjKk2X+Zs6c2ef+VlQlWKq6Buh16qCq+4Hr3eNdwPkBh5aWpUuX8tJLL/HTn/4026HkvcCSqqr+o4j8J14JP4C/iutYPWzOJRrbopwx7uSv/fhxb0xEdXV10v3LiyMURYRjzTaqyqRn3bp1aR8jn7qa+SXoLlWdeDVU1T0edhJP/9eu9W42p+qn2lNUxa6pmjSl2wvgvPPOY+HChXzrW9+iqqqKhQsXct5552UouvwR2I0qEfkfeJcAxuBdX/2ViAy4nmo+OFmg+uTd/yuuuIIrrriiz9fVlBdTZ0nVpCmdfqXnnXceN954I6NGjQK8Xiu//e1veeuttzIUXf4I8u7/XcB8Vf2Gqn4dWAD8twDfP+vao910dmlPhSqA6dOnM3369D5fN6aihCNNw+ZenvFJOr0AFi5cSHFx8SnrOjs72bBhQyZCyytBJlXBKz4R00UeVPIfjAY3RLUy7vT/6NGjHD16tM/Xja8s5XBDm6+xmfy3dOlSli9fTnn54Ps7p+pHHetnbU4KMqn+Eu/u/zdF5FvARuD+AN8/6xKLqQA88cQTPPHEE32+bkJVCYcb2+nqztlZt01ILF26tM8CPqmkSp5DOVa+CyypuhlQPw/UAUeBO1T1B0G9fy5IrPoP3mnVwoUL+3zdhMpSurqVOrsEYDLktNNOG9T+GzZsoKPj1Ov6HR0dPP3005kMKy/4fvdfRBrx7vb3rIrbpqpa6XcMuaIpSUt16tSp/b5unOvXerChreexMem49957+yyWkih2Q2rhwoVUVVVRX1/Phg0baGho8CvE0Aqi9N9Iv98jLBKnUgE4fPgwAOPGjUv5ugkukR5qsJaqyYylS5dy++230909sFniZ8+ezfvvv3/K3f6ysjKWL09anGtYy/YcVcNKstP/devW9dspe0LVyZaqMZnyxS9+cUD7/epXv2Lbtm0sX76c0047DRHhtNNOY/ny5VYDIAmrpxqgk3f/T57+X3PNNf2+bkxFCZEC4VC9JVWTOT/5yU8A+hy6+qtf/aoncS5dutSS6ACEpqUqIveLyGERydlqQP2JtVQr4lqqkydPZvLkyX2+LlIgjK0osZaqybif/OQnFBQkTwMFBQWWRIcgNEkVeABYlO0g0tHYFqW8OEIkblqUgwcPcvDgwX5fO76yhEOWVI0PUl1XHej1VnOq0CRVV1X9WLbjSEdTe+cprVSAp556iqee6r9q3PjKUkuqxhepulcNttuV8YQmqQ6EiCwTkU0isunIkSPZDqcXr0D1qVOiLFq0iEWL+m+AT6gq5aBdUzU+SDZ81aZLGbq8SqqqulxVa1W1duzYsdkOp5fEqv8AEyZMYMKEZNMlnWp8ZSkNbVFaO6zIsMms2PBVu7OfGXmVVHOdV/bv1Jbqvn372LdvX7+vnVBp3aqCICLfE5F3RGSLiKwRkVEp9lskIjtEZKeIfDXgMDNu6dKl7N69m+7ubpsuJU2WVAOUrKW6fv161q9f3+9rY7MF2CUA360H5qjqXOBd4GuJO4hIBPgxcB0wG7hVRGYHGqXJWaHppyoiDwNXAGNEZC/wDVW9L7tRDU5dcwfVCdNMX3/99QN67YSqEgAON1pS9ZOqPhP39BVgSZLdLgZ2umlVEJFVwGKG0bRAJrXQJFVVvTXbMaTjREsH9a2dnFZzatm1voanxrOWalZ8AXgkyfrJwIdxz/cC85MdQESWAcsgvSLRJjxCk1TDbnddCwCnjT71LuuHH3r/N/srrDKytIjy4ohdU80AEXkWSHZ38B5VXev2uQeI4s1WMWSquhxYDlBbW2u1G4cBS6oB2VPXDMD0Mae2VGOV01PNURVvfJX1Vc0EVb26r+0icifebL8LVTVZItwHxP8VnOLWGWNJNSi7j3ot1Wk1p7ZUP/GJTwz4GONHllqlKp+JyCLgy8DlqtqSYrfXgDNFZAZeMr0F+FxAIZocZ3f/A7KnrpmJVaWUFkVOWT9mzBjGjBkzoGPYAIBA/CswElgvIptF5GcAIjJJRNYBqGoUuBt4GtgOrFbVbdkK2OQWa6kGZHddc6/rqQC7d+8G6HfyP3BzVTW20d2tFBQMq+m9AqOqZ6RYvx+4Pu75OqDvmo1mWLKWakD21LUwfXTvCdeef/55nn/++QEdY9KoUjq7lEPWrcqYnGUt1QA0tHVS19zBaUmS6uLFiwd8nHMmejPPvL2/gYlVIzIWnzEmc6ylGoAPXHeq6UlO/6urq6murh7Qcc6ZWIkIbN1n8wIZk6ssqQZgd4ruVAC7du1i165dAzpORUkhM8aUs3W/zbVuTK6y0/8A7EnR8R/ghRdeAGDmzJkDOtacSVVs2h3qsrLG5DVLqgHYfbSZcSNLKCvu/eu+6aabBnWsOZMrefzN/dQ1tTO6oiRTIRpjMsRO/wOQ6s4/QFVVFVVVVQM+1pxJ3r7b9tt1VWNykSXVAKTqowqwc+dOdu7cOeBjneuSql1XNSY32em/z1o6ohxubE96kwrgxRdfBOCMM5L2Oe+lqqyIqTUj2GY9AIzJSZZUfRYb85+qpbpkSbJynX07b3KVtVSNyVF2+u+zzR+eAODsCZVJt1dUVFBRUTGoY547qYo9dS3Ut3amG54xJsMsqfrslV11jB1Zwuljk5/+79ixgx07dgzqmHMme9dV37abVcbkHEuqPlJVXt5VxyUzRyOSvADKyy+/zMsvvzyo486Z5LV63/jgeNoxGmMyy66p+uj9I80caWznktNHp9zn5ptvHvRxR1eUcM7ESl549whfunJgN7iMMcGwlqqPXt5VB8AlM1Mn1bKyMsrKkt/E6svls8by+p7jNLbZdVVjcoklVR+98n4dE6tKU975B9i+fTvbt28f9LEvnzWWaLfyh/fr0gnRGJNhllR9oqq80s/1VICNGzeycePGQR//wtOqqSgp5PfvHkknTJNARL4nIu+IyBYRWSMio1Lst1tE3nKzA2wKOEyTw+yaqk/ePdREXXMHC/q4ngpwyy23DOn4xYUFfPT00fx+xxFUtc/EbQZlPfA1VY2KyHeBrwFfSbHvlap6NLjQTBhYS9UnL7/v/V/r63oqQGlpKaWlpUN6j8vPGsu+E628f6R5SK83vanqM24OKoBX8GZKNWbALKn64HhzBz9/YRezxlcwtabvm1Bbt25l69atQ3qfy84cC2CXAPzzBeA/U2xT4BkReV1ElqU6gIgsE5FNIrLpyBH7dxoOQpVURWSRiOwQkZ0i8tVsx5OMqvIP//4mdU0dfP/mef3uv2nTJjZtGtoluak1ZZw+tpxVr35gs6wOgog8KyJbkyyL4/a5B4gCK1Mc5lJVvQC4DviSiFyWbCdVXa6qtapaO3bs2Ix/FpN7QnNNVUQiwI+Ba4C9wGsi8riqvp1s/25VWjqiyTZllKr3s7Orm/0n2njyrf1seOcw37xxds/Ip74sXbo0rfe/54ZzuPvXf+QT//u/+N5nzufcSZWMKIoQyaPZVpPVoU2Hql7d13YRuRP4BLBQNfYv3OsY+9zPwyKyBrgYeCGjgZpQCk1SxfvS7lTVXQAisgpYDCRNqtv2NzD7608HGN5Ji86dwB0fnT6gfYuKitJ6r6vOHs/jd3+MLz70Op//5WtpHStX7f7ODYG9l4gsAr4MXK6qLSn2KQcKVLXRPb4W+HZgQZqcFqakOhn4MO75XmB+/A7u2tYygNGTZ/C1684OJDARiBQUMKGylMnVIzhvctWA78Zv2bIFgLlz5w75/c8YN5K1d1/Khu2HaGyL0trRRXfyBpbp378CJcB692/4iqr+lYhMAn6hqtcD44E1bnsh8GtVfSpbAZvcEqak2i9VXQ4sB6itrdUvXn56liPq3xtvvAGkl1TBmxRw8bzJmQhpWFPVpON+VXU/cL17vAs4P8i4THiEKanuA6bGPZ/i1oXabbfdlu0QjDEZFKa7/68BZ4rIDBEpBm4BHs9yTGmLRCJEIpFsh2GMyZDQtFTdCJe7gaeBCHC/qm7Lclhp27x5MwDz5s3LahzGmMwITVIFUNV1wLpsx5FJllSNyS+Sohte6InIEWBPFt56DBCW8eBhibVUVedkO4h0iUgjMLhpHoIXlu9EtuM8qqqLkm0IVUt1MFQ1K8NXRGSTqtZm470HKyyx5lEVqB25/vsO03ciV+MM040qY4zJeZZUjTEmgyypZt7ybAcwCGGJNSxx9icMnyMMMUIOx5m3N6qMMSYbrKVqjDEZZEnVGGMyyJKqDwY6eVy2ichnRGSbiHSLSE52T8nlwuT9xSYiJSLyiNu+UUSmu/XTRaTVTRq4WUR+FveaC92EgjtF5EeS5uRjacS4NC6+ze47Ms9te94dM7ZtnM8xXiYib4hIVESWJGy7Q0Tec8sdcesz+nscFFW1JcMLXn3NQvf4u8B3sx1TijjPAc4Cngdqsx1PkvgiwPvATKAYeBOYne24Bhob8NfAz9zjW4BH3OPpwNYUx30VWAAI3lQu12UjxoR9zgPej3uese/LAGOcDswFVgBL4tbXALvcz2r3uDrTv8fBLtZS9YGGZPI4Vd2uqrk8wqenMLmqdgCxwuS5YCCxLQYedI8fBRb21WISkYlApaq+ol5mWAF8KgdivNW91g/9xqiqu1V1C9Cd8NqPA+tV9ZiqHsebCXeRD7/HQbGk6r++Jo8zfUtWmDxXisYOJLaefdwf2XogNr3uDBH5o4j8XkT+LG7/vf0cM8gYYz4LPJyw7pfu1P9/pnlqnc6/carXZvr3OCh5O0zVbyLyLDAhyaZ7VHWt26e/yeN8N5A4TeAOANNUtU5ELgT+Q0TOzXZQyYjIfKBFVeOn/F2qqvtEZCTwG+A2vNagwZLqkGkGJo8LQn9x5rhcLkw+kNhi++wVkUKgCqhz34d2AFV9XUTeB2a5/eMvFaX7eYccY9z2W0hoperJSQ8bReTXeKfwQ02q6fwb7wOuSHjt82T+9zgodvrvg7jJ4z6pKSaPMwOSy4XJBxLb40DsjvQS4HeqqiIyVrzZgRGRmcCZwC5VPQA0iMgCd0p9O5DO2cSQY3SxFQA3E3c9VUQKRWSMe1yE13DYytCl82/8NHCtiFSLSDXeDeKnffg9Dk5Qd8SG0wLsxLvWs9ktP8t2TCnivAnvelM7cMh9IbMeV0KM1wPv4t0hvifb8fQXG96sqp90j0uBf3ffh1eBmW79nwPb3HfjDeDGuGPW4iWp9/EmIZRsxOi2XYE38WH88cqB14Et7jP8EIj4HONF7nvajNeK3hb32i+42HcCn/fr9ziYxYapGmNMBtnpvzHGZJAlVWOMySBLqsYYk0GWVI0xJoMsqRpjTAZZUjXGmAyypGqMMRlkSdUYYzLIkqoxxmSQJVVjjMkgS6rGGJNBllSNMSaDAk2qQ53gS0TmicjL4k1St0VEPhtk3MYYM1CBValy9SPfBa7BK+P1GnCrqr4dt890oBL4B+BxVX3UrZ8FqKq+JyKT8EqPnaOqJ1K9X2FklpYVf8GnT2MSNbR+NemUGldfW6Z1dV1JX7P5jY6nVXWRr4HlkEWLFulTTz2V7TBMZqScQibIyv89E3wBiEhsgq+epKqqu922Uyb4UtV34x7vF5HDwFjgRKo3s9rQuaGurpvnX5qWdNuoETvHBBxOVh09ejTbIZgABJlUk03SNX+wBxGRi/Gmsn0/ybZlwDIAYdSQgjQZpiBdduk+Xzz33HMAXHnllVmOJHeFao4qN/XsQ8Adqpo4XS2quhxYDhApmGLVt3OAANKdzmSbJpc0NDRkO4ScF2RSTWsSNxGpBJ7Em27hlQzHZvyiINFsB2EyZfHixdkOIecFeV425Am+3P5rgBWxm1cmJBQKupIvJne0dXax/0RrtsPIC4ElVVWNAnfjzYC4HVitqttE5Nsi8kkAEblIRPYCnwF+LiLb3MtvBi4D7hSRzW6ZF1TsJg0K0qVJF5M7HvzDbq75/u9p7ej7r92zzz7Ls88+G1BU4RToNVVVXQesS1j39bjHr3HqfN2x9b8CfuV7gCbjBJCoJdBcd6C+jeaOLt7ce4IFM0en3K+11Vqz/QnVjSoTQgpip/o5r6ndu/C9afexPpPqjTfeGFRIoWV9XYy/3I2qZEsuEZH7ReSwiGyNW/dNEdkXd8np+rhtX3MjA3eIyMezE3XmNLV5/yCv7T6e5UjCz5Kq8V1Irqk+ACQb3fUDVZ3nlnUAIjIb70brue41P3EjBkMr1lJ9Y89xurpT/9s888wzPPPMM0GFFUqWVI2/3Ol/siWXqOoLwLEB7r4YWKWq7ar6J2An3ojB0GpsjyLi/XznYOq+qJ2dnXR2dgYYWfhYUjX+UiCaYgmHu10Rn/tFpNqtSzY6cHKyF4vIMhHZJCKbjhw54nesQ9bY1sn5U0YBsKmPSwA33HADN9xwQ0BRhZMlVeMrAaRLki4h8FPgdGAecAD458EeQFWXq2qtqtaOHTs2w+FlTlNblLMnjGRiVSmv7R5og90kY0nV+EuB7hTLAIjIbhF5y90o2uTW1YjIehF5z/2sdutFRH7kbiBtEZEL0gpd9ZCqdrkh0f/GyVP8tEYH5qKm9igjSwupnV7Da7uPkap63VNPPYVV2uqbJVXjuwy0VK90N4pq3fOvAhtU9Uxgg3sOcB1wpluW4bU0hx63V2si5iYg1jPgceAWESkRkRnu/V5N572yqatbaenooqKkiIumV3OooZ29x60/6lBZP1XjLwUyf6q/GLjCPX4QeB74ilu/Qr1m1isiMkpEJqrqgf4OKCIPu2OOcaP6vgFc4UbuKbAb+CKAGwm4Gq9sZRT4kqrm2K23gYt1p6ooLeSCad5l4zf3nmBqTVmvfRctGjblb4fMkqrxlwKp082Y2Cm9s9xVGks8wjMiosDP3fbxcYnyIDDePU51A6nfpKqqtyZZfV8f+98L3NvfccOgsd27mz+ypJCJVaUA1DV1ZDOkULOkanwmkLqe6tG4U/pULlXVfSIyDlgvIu/Eb1RVdQnXDFGsj+rI0kKqRhQhAnXNyZPqk08+CWA9APpgSdX4S0GiQz/9V9V97udhEVmDd7PoUOy03l33POx2z7sbSEGIP/0vjBRQNaKI4ymSalFRUZChhZLdqDL+il1TTbb0Q0TKRWRk7DFwLd7NoseBO9xudwBr3ePHgdtdL4AFQP1ArqcOd42xpFritbFqyoo51pI8qV577bVce+21gcUWRtZSNf4b+o2q8cAaEQHvu/prVX1KRF4DVovIXcAevNKQ4FVAux5vhFML8Pl0wh4uGuNO/wGqy4tTtlRN/yypGn+poEOco8pNEnl+kvV1wMIk6xX40pDebBiLnf6PLPVO7avLitl7PPnEmb/97W8Bq1bVF0uqxn/hGD01bDW5u/89p//lRby1L3lLdcSIEYHFFVaWVI2/FFBLqrmsqc0rplJW7BXa8k7/O1FV3KWXHldffXU2QgwVS6rGZ312qTI5oKEtSkVJYU8CHV1eTEdXN80dXT2tVzNw9m03/krj7r8JRlN7lJFxybO6rBgg6c2qtWvXsnbt2l7rzUn2Z8j4bqg3qkwwmtqiPTepAGrKvaR6rLmj11DVysrKQGMLI0uqxl9qp/+5rqk9SkVpXEs1llST9FW98sorA4srrCypGl+pgtqpfk5rbOtklDvlB6/zPyQ//Tf9syaE8V9XQfLF5ITGVC3VJEn1scce47HHHgsstjCylqrxVxqd/00wmtqiVMYl1crSQiIFkjSpjh6devpq47Gkanyn3ZZUc1lTe/SUrlMiQnVZMceTXFO9/PLLgwwtlCypGn/5U6TaZEh81f94NeVFSVuqpn+WVI2/7PQ/p8WX/YtXXeaNqkr06KOPArBkyRL/gwspS6rGd2rDVHNWfNX/eDXlxbx3uKnX/hMmTAgkrjALtAkhIotEZIeb7fKrSbZfJiJviEhURJYkbLvDzZ75nojckfhak8Ps7n/Oakoo+xdTk6L836WXXsqll14aSGxhFdg3W0QiwI/xZrycDdwqIrMTdvsAuBP4dcJra/AmYpuPV/n9G7FpiU1ui/VTTbYMhIhEROSPIvKEez5DRDa6P8yPiEixW1/inu9026f796nyR6rT/5py70ZVd7fNVDNYQTYXLgZ2quouVe0AVuHNftlDVXer6hZ6zwr/cWC9qh5T1ePAesCmdQwFQbsLki4D9D+A7XHPvwv8QFXPAI4Dd7n1dwHH3fofuP1MPxKr/sdUlxXTrdDQdup11dWrV7N69erA4gujIJNqqpkuM/ZaEVkmIptEZJNq85ADNRmk3tj/ZEt/RGQKcAPwC/dcgKuAR90uDwKfco8Xu+e47QslsW6d6SWx6n9MTYoBAFOmTGHKlCnBBBdSeXWjyk1fvBwgUjDFzltyQXp3//8F+DIw0j0fDZxQ1ah7Hv/HtecPr6pGRaTe7X90qG8+HCRW/Y+JjapK7Kv60Y9+NJjAQizIlmo6M13aLJkhpUB3tyRdgDGxMwu3LIu9TkQ+ARxW1dezFftwkFj1PyY2/v9Ykm5Vpm9BtlRfA84UkRl4CfEW4HMDfO3TwP8Td3PqWuBrmQ/RZJz2WfrvqKrWptj2MeCTInI9UApUAj8ERolIoWutxv9xjf3h3SsihUAVUJehT5G3GhOq/sdUl3st18QeAA8//DAAt956azABhlBgLVX3n+BuvAS5HVitqttE5Nsi8kkAEblIRPYCnwF+LiLb3GuPAf+Il5hfA77t1pmcN7QbVar6NVWdoqrT8f4A/05VlwLPAbHudonTU8e62i1x+9sloH40JlT9j6lJUf5vxowZzJgxI7D4wijQa6qqug5vGuH4dV+Pe/waXusj2WvvB+73NUCTeZkv/fcVYJWI/BPwR+A+t/4+4CER2Qkcw0vEph9N7VEqE66nAowoilBSWNDrRtWCBQuCCi208upGlck9SvojqlT1eeB593gXXve8xH3a8M5wzCA0tUWTzkMlItSUF9v4/yGwpGp8ZmP/c1li1f943vj/U5PqypUrAVi6dKnvsYWVJVXjL7XSf7kssep/vOryol5dqmbNmhVEWKFm33bju+4uSbrkGhG5X0QOi8jWuHU1IrLe1ZxYH+uBIp4fuWGxW0TkguxFPnSJVf/jjRpRzInWU7tUXXTRRVx00UVBhBZallSNr1TTHqYapAfoPfz5q8AGVT0T2OCeg1fD4ky3LAN+GlCMGZVY9T/eqLIiTrRYP9XByslvtskvYUmqqvoCXs+BePHDXxOHxa5Qzyt4/WcnBhJoBjWmuFEF3jXVEwlFVVasWMGKFSuCCi+U7Jqq8ZdCd7hvVI1X1QPu8UFgvHucqh7Fgbh1uFFiywCmTZvmb6SDFO3qprWzd9X/mFFlRXSrd4mgaoS3z7nnnhtkiKFkSdX4TNDu3Lt+OhSqqiIyqAEF8fUoamtrc2owQnN7F9C7mEpM7AbWiZaOnqR64YUXBhNciIW6CWFyn7qWarIlJA7FTuvdz8NufejrUcTK+qXuUuUlUruuOjih+Wab8ArLNdUU4oe/Jg6Lvd31AlgA1MddJgiFnqr/Ka6pjnJJNb5b1QMPPMADDzzge2xhZqf/xl/aU5Eq54nIw8AVeNWz9uLNNvEdYLWI3AXsAW52u68Drgd2Ai3A5wMPOE0np1JJdU01dvp/sqU6b9483+MKO0uqxleZGKYaFFVNVXppYZJ9FfiSvxH5q7Hf0/+T11RjLKn2z5Kq8Z0NU81NqaZSiYn1Xz0e11Lt6vJubkUikaSvMZZUjd9U6A7P9dNhJXb6n6rzf2GkgMrSwlNaqg899BAAd955p+/xhZUlVeO7EN2UGlYaU8ykGq+6/NShqhdcEMrRuIGypGp85XWpCsc11eGmqS1KpEAYUZT6VH7UiKJTTv/nzp0bRGihZk0I4zPv9D/Z0u8rRUpF5FUReVNEtonIt9z6GSKy0RUzeUREit36Evd8p9s+3d/PFm5N7cmr/scbVVZMfdzpf2dnJ52d1m+1L5ZUjb80rX6q7cBVqno+MA9Y5PqEfhf4gaqeARwH7nL73wUcd+t/4PYzKTS0daa8SRVTXXZqS3XlypU9NVVNcpZUja+82VSH1lJ1xUqa3NMityhwFfCoW59Y5CRW/ORRYKH01Qwb5praoimHqMaMKis+pfN/bW0ttbWp5mo0kIGkKiJfyUQgJk9p8lqq7jpryimqY0QkIiKb8YaHrgfeB064iSThZCETiCty4rbXA6N9/Xwh1tQ+kKRaRGNblGhXNwBz5sxhzpw5QYQXWoO+USUiq+Of4p2W2WmWSSrWUk2hrymqvderdgHzRGQUsAY4O6MBDmONbVHGVCSv+h8zyhVSqW/tZHRFCW1tbQCUlpb6Hl9YDaWl2qCqN7vlM8CzmQ7K5BdVSboM7hh6Am966kvwapfGGgTxhUx6ipy47VVAXQY+Ql7y5qdKPkQ1ptpNVR27rrpq1SpWrVrle2xhNpSkem/C83syEYjJU5rW3f+xroWKiIwArgG24yXXJW63xCInseInS4DfueGkJonGAV5TBahv9a6rzp8/n/nz5/seW5j1e/rvuqV8CTgdryr6ZhH5raruAVDVxErpxpwijWGqE4EHRSSC1wBYrapPiMjbwCoR+Sfgj8B9bv/7gIdEZCfed/WW9CLPb41tnSkrVMXETv+PN3st1XPOOcf3uMJuINdU1wI/Ap4C7se7TPZ/isgTwN+paruP8ZmQU2XIVapUdQvwkSTrdwEXJ1nfBnxmSG82zHREu2mPdg+gS5UrquJGVbW0tABQVlbmb4AhNpAmRERV71PVDcAxVf1veK3W3biK5sakNvTTf+Of5p6yf/20VMtjhaq90//Vq1ezevXqvl4y7A2kpfqsiNytqv+K10qNdVf5noi862t0IVOsXqKI4LXMOvC6oXQNbgaOvOK1VC2B5pqT4/77vlE1sqSQSIH09FW95JJLfI8t7AaSVP8O+JqIbAImub6ELXh3Ye3OqulXiKZOGTYa210t1X5O/0WEUSNOTlV91lln+R5b2PWbVFW1G7hXRH4AXI3XL7Ua2Ird+Qcg4roHTewuB2BqdwUAx8S73PxBQSMAHdLFSPWuURW6Ky8Nbp9W6Qou4CBZ6b+c1NQ2sNN/8AYAxJJqU5M3wK2iosK/4EJuwJ3/VbUFr8vK40N9MxFZBPwQiAC/UNXvJGwvAVYAF+K1gj+rqrtFpAj4BXCBi3mFqv6vocYRlInntnDWlY2UVnXT1CC89HwhB7f13dk63/TT+d9kSeOgkurJoaqPPuqNDrZ6qqkFVvrPdYv5MV5fw73AayLyuKq+HbdbT0EMEbkFb6TWZ/Hu6Jao6nkiUga8LSIPq+ruoOIfrNHnNjH9hnoKXQ4dWaUsvK6TP2gHf3p7OCVWa6nmoliB6v5O/8ErqrL/hDeS6tJLL/U1rnwQZD3Vi4GdrjsMIrIKrwBGfFJdDHzTPX4U+FdXEEOBcjdKZgTQATQEFHe/YjematQbuvexghKiV+2FhNxZVAwLruigfMsERql3g+C9iHdpYGfkBAAd0h1M0AGJTVFtcktje/8FqmNGlRXz9n7vv9sZZ5zha1z5IMhve0+xCye+EEavfRIKYjwKNAMHgA+A/y/ZoAMRWRYrzqHanPlPMBhV0aSri6vy9NppH1STLyZ7YpP+VfZz9x9OLVRdX19PfX29r7GFXViaEBcDXcAkYAbw9yIyM3EnVV2uqrWqWitSHnSMp6pP3gLoqB9mE6YpdHUXJF1M9jS1RSksEEoK+/93qC4vprWzi7bOLtasWcOaNWsCiDC8gjz97yl24cQXwkjcZ29CQYzPAU+paidwWEReAmqBXb5HPQBdXvddmsX7ax4pUPT3Y9AbDtId9xuWTqHkuXHUFpRwwDVky9w/QbF6yTbvTv/tmmpO8oqp9F31P2bsyBIADje0c9lll/kdWugF+W1/DTjTTYVRjDcuO7EnQaqCGB/gFSZGvCboAuCdQKIeosK3K5nxZhHFLQIKRc0FjH5hFGyrynZowXLXVJMtJnsGUkwlZlLVCAD217cyc+ZMZs7sdZJo4gTWUlXVqIjcDTyN16XqflXdJiLfBjap6uOkLojxY+CXIrINr4brL9248Jw2dl8hY/cV0tzsjZNubBzB0SzHFDRrqeamxrYoFSX9X08FmDjKuwF7oL6V48ePA1BdXe1bbGEX6GyqqroOWJew7utxj5MWxHBTauRkoYzK7mKqXIf+Cd1e8pw10xtodvnt6wEomuBd2K9/cxrP/eZKAN7YWQnAPvE6U3fka+d/sOunOaipvf8KVTE9LdUTbaxd65VPtn6qqdkU1cZXNvY/NzW2RZlQObDq/SOKI4wqK+JAfSt/ecUV/gaWByypDlGseMo4LePCzlEAXDDaG3XykY++CUDBJ3YD0DTVa4WWHaqkvNwrnVajNe6n1wrodl2MjtIK5FMRlp75qAb/SpGpeCPsxuP1VV6uqj8UkRrgEWA6XrW0m1X1uOvT/EPgerz6FHeq6htpf4Q8FLtRNVATq0Zw4EQb06ef52NU+cGaEMZXml6Xqijw96o6G+/m5JdEZDbwVWCDqp4JbHDPAa4DznTLMuCnmf48+aKpLTqg0VQxk6pK2V/fxtGjRzl6dLjdGRgcS6rGd2lMUX0g1tJU1Ua8qVQmc+pU1IlTVK9wU1u/gjeX1cQMf5y84N39H9iNKvBuVh2ob+WJJ57giSee8DGy8LPT/yGK9U1tJ0rsFlNFmVdxqrDQW1N41Ot7Gqk52eG/yG2bVuAlleqod8PqgwLvMsDrhV4roE7afIw+SEL3ICf5S3oUb1qfjwAbgfGqesBtOoh3eQBSj9o7QBpEZDfQiDcAJeoNLkl+CSKd9wlKe7SLjq7uAXepAu/0/0RLJ5dedgUlRcNsAMsgWUvV+Mob+y9JF2BMbFixW5YlO4aIVAC/Af5GVRtOPb4qrni6z65U1XlxU2qnugSR83oKVA/m9N91q5KRY5g6dWo/ew9v1lI1vlL67FJ1NC5JJeXKPv4GWKmqj7nVh0RkoqoecKf3h936gYzay5TFwBXu8YPA88BXfHqvjBpMLdWYia5b1Xu79zFSqxk3bpwvseUDS6pDFLs7f6ighZeKvP/TB/Z7p/J7H7wGgEX7xwAw67K3AKg/UE1Li/cXv9WNRq1zw1I7XWOrMN9OHtKY+M/dzb8P2K6q34/bFBt59x16T1F9t6uANh+oj7tMkA4FnhERBX6uqstJfQki8TMsw7tpxrRp0zIQSvqG1FJ1SfWtV57n4LYS66faB0uqxleK0DXEpAp8DLgNeEtENrt1/xdeMl0tIncBe4Cb3bZ1eN2pduJ1qfr8UN84waWquk9ExgHrReSUIdKqqi7h9uIS8HKA2tranOgnd7jRu14fG9M/EOOrvH0LppzPNRfb6X9fLKmmqUO6ORDx+p7WFXhf1t3d3l/1zevOB+CKV+YAMGFsPU3NXku13f332l7oXSI8UuD1T22UjmACD4oOfUSVqr4IpMrIC5Psr8CXhvRmfcexz/08LCJr8KqmpboEkfP2n/C+a5NHjRjwa0oKI4ypKOZIdxmTJydW7DTx8uxc0+Qapc8bVTlPRMpFZGTsMXAt3vxs8cV/4i9B5Lx9J9ooighjKgbeUgXvuurRI4c5ePCgT5HlB0uqxl/qnf4nW0JiPPCiiLwJvAo8qapP4V2CuEZE3sObEPM7fRwjp+w/0crEqhEUFAzu32BiVSmVR7fy1FNP+RRZfrDT/wyK1ULdF/FmHWh09VX31HszT844VsPp4v3KmxNqsJ5ws6rmz/BUj3f3PzQJtBc3/c/5SdbXkeQSRBgcqG9lYtXAxv3HmzRqBOvfn8aXF33Uh6jyh7VUje9C3lLNO/tPtA3qemrMxKpS9rWXUDZqtA9R5Q9LqsZX3th/S6q5ItrVzcGGNiYNJamOGsEYaWbru7szH1gesdN/HzUUeHfyd4g3evFQQQvvu54BsT9n9e5uf76d9sezSf5yx+HGdrq6dUhJdVJVKbVFH/Lqi/V8dO6ZPkSXHyypGl+F/Zpqvol1p4oNOx2MiaNG8ErnacyfZQm1L5ZUAxBrhdZJW88NqYjrfplvE/31otBlLdWcsb/e60s9lJbquJElNDCC492DT8jDiV1TNb5ShC5NvpjgxVqqQ7n7XxQp4Mzydg4e8KucQn6wpGp8pXgt1WSLCd7+E61UlhYOqpZqvPPkQ2T/1gxHlV/s9D9gsUsBXYFUq8sNlkBzx/4TrUM69Y9pnzSPXUebMxhR/rGWqvGV4lV2TraY4O0bYh/VmLFjxvKnxgLUunSkZEnV+M6Sau44UJ9eS3W01jMyeoKG1mgGo8ovdvpvfBW7pmqyr7k9yomWTiYOoTtVTNuHb/GRwkYONLRSVTa067L5zlqqxld2+p87DtQPvuRfoosuv5YXO2dw4ES+zKGWeZZUje+6Uyz9EZH7ReSwiGyNW1cjIutF5D33s9qtFxH5kYjsFJEtInKBH58lzPadGHof1Zgzp46nSUvY7xK06c2SqvFVmi3VB4BFCetSTbh3HXCmW5YBP00r8Dx0cjTV0JNq0+F9TI40cLDeWqqpWFI1PlO6Uiz9vlL1BeBYwurFeBPt4X5+Km79CvW8AoxyFfmNc+BEKwUC4wcxjUqil156kdqSg+y30/+UAk2qIrJIRHa4U7ReU/qKSImIPOK2b3Rzvce2zRWRl0Vkm4i8JSI2Vi4E+mmpDmiK6gSpJtybDHwYt99et8447x9pZkp1GYWRof+3v+mmmzhQPZeDDXb6n0pgd/9FJAL8GLgG7wv/mog8rqpvx+12F3BcVc8QkVuA7wKfFZFC4FfAbar6poiMBjqDit2kJ2UFLu1/iuq+9DXhnult+8EGzpk4Mq1jVFVVMbammu0HGjIUVf4JsqV6MbBTVXepagewCu+ULV78qd2jwEI3TfG1wBZVfRO8quuqajeQQ8BrqQ7t9D+FQ7HT+oQJ9/YB8dN8TnHrDNDW2cXuo82cNaEyrePs3LmTCZzgQH2bDQBIIcikOpDTs559VDUK1AOjgVmAisjTIvKGiHw5gHhNBviQVFNNuPc4cLvrBbAAqI+7TDDsvXeoiW6Fcyak11J98cUXKTz6Lq2dXdS32sliMmG5UVUIXAosdT9vEpFe8wOJyLLY9TlVG5+cK7ok+dIfEXkYeBk4S0T2ishdpJ5wbx2wC9gJ/Bvw1z58lNDaftA7XT97Ynot1SVLljD7o9cA2M2qFIIcUTWQ07PYPnvdddQqoA6vVfuCqh4FEJF1wAV4XWp6qOpyYDlApGCKnZvkAAWiQ2yVquqtKTb1+oOq3rnol4b0RsPAOwcaGVEUYVpNWVrHqaioYNp4r4V6sKGV2ZPSS9L5KMiW6mvAmSIyQ0SKgVvwTtnixZ/aLQF+5/6zPA2cJyJlLtleDryNyXmaRpcqkznvHGxg1oSRRAY5LXWiHTt20F7ntYWspZpcYC1VVY2KyN14CTIC3K+q20Tk28AmVX0cuA94SER24vVPvMW99riIfB8vMSuwTlWfDCp2k558nn8rDFSVdw42cu3s8f3v3I+XX34ZBSIF42wAQAqBFlRR1XV4177i13097nEb8JkUr/0VXrcqEyIKdFurNKuONLVzrLmDs9K8SQVw8803A/BvP3zZhqqmYFWqjO/sVD+73jnQCMDZaXanAigr867JTq4ewW4rVp1UWO7+m5Cya6rZ907szn8GWqrbt29n+/btLJg5ms0fnuBYc0fax8w3llSNrxSISnfSxQTjnQONTKgspbq8OO1jbdy4kY0bN3Lt7Al0K/zuncP9v2iYsdN/4ztrlWbX2wcaODvN4akxt9xyCwAlJSVMrCrlmW0HWXLhlIwcO19YS9X4KtZPNdli/Ldp9zHeOdjIgpmjM3K80tJSSktLERGumT2eF947QmuHjRiPZ0nV+Mq7ptqddMkH/VVeyyZV5Z+e3M64kSXcfslpGTnm1q1b2brVqxl+7ewJtHV281/vHcnIsfOFJVXjKxWIiiZdwi6u8tp1wGzgVhGZnd2oTnryrQNs/vAE/3DtWZQVZ+ZK36ZNm9i0aRMA82fWMLK0kPVvH8rIsfOFXVM1vsuXVmkSPZXXAEQkVnkt6Wi/LlUa2vwpQhJfMKqxrZMP6lr47lPvcPaEkfx5Bq95Ll26tOdxUaSAq84ex/rth1jzx718ZGo11WXF4AZtSXqDt3JaZWnqSQ8tqRpfKUo0f5Nqsspr81Pt/Pb+BuZ+8xnfg4opLBAe+PzFaQ9NjVdUdGoyuf2S0/jdO4f520fezNh7hMHu79yQcpslVeO74Tyiys1msAxgzOQZ/N83nOPnewFQVhzhtJoyTh9XwfjKzE6QsWXLFgDmzp0LwIWn1bD569ey42AjW/aeoNndtBrOtVYtqRpfef1U8/Y/WL+V1+Irp9XW1upf/tnM4KLzwRtvvAGcTKoAkQJh9qRKq1jlWFI1vsrz0/+eymt4yfQW4HPZDclft912W7ZDyHl299/4SoFOupMuA5HjXZaiQKzy2nZgtapuy25U/opEIkQikWyHkdOspWp8pSidMrTO4QOcLDKrklVey2ebN28GYN68eVmNI5dZUjW+UqBj6Kf/g+qyZPxnSbV/kq936UTkCLDH57cZAxz1+T3CEsM7qroocYOIPOW2J1MKxFc6Xu5u7MReuwRYpKp/6Z7fBsxX1bszFnmARKQR2JHtOPqRC9+ngch2nEeTfd8hj1uqqjrW7/fwJhgc+rz1eRZD0i9YqvXD1I5s/1v1Jxe+TwORy3HajSqTywYyWaQxOcWSqsllA5ks0pickren/wFZ3v8uvsvbGFJNFunHewUkF/6t+hOGGCGH48zbG1XGGJMNdvpvjDEZZEnVGGMyyJJqGkTkeyLyjohsEZE1IjIqS3F8RkS2iUi3iATazSSXh5H6rb/PLiIlIvKI275RRKa79dNFpFVENrvlZ3GvuVBE3nKv+ZFIelVJ04hxaVx8m913a57b9rw7ZmzbOJ9jvExE3hCRqOu7HL/tDhF5zy13xK3P6O9xUFTVliEuwLVAoXv8XeC7WYrjHOAs4HmgNsD3jQDvAzOBYuBNYHa2/11y5bMDfw38zD2+BXjEPZ4ObE1x3FeBBXilnv8TuC4bMSbscx7wftzzjH3PBhjjdGAusAJYEre+Btjlfla7x9WZ/j0OdrGWahpU9Rn1imoAvILXjzIbcWxX1WyM1OkZRqqqHUBsGOlwMJDPvhh40D1+FFjYV4tJRCYClar6inqZYQXwqRyI8Vb3Wj/0G6Oq7lbVLdBrvPPHgfWqekxVjwPrgUU+/B4HxZJq5nwB7y/icJKs8v3kLMUStIF89p593B/feiA2rekMEfmjiPxeRP4sbv+9/RwzyBhjPgs8nLDul+7U/3+meWqdznco1Wsz/XscFOun2g8ReRaYkGTTPaq61u1zDxAFVmYzDhMaB4BpqlonIhcC/yEi52Y7qGREZD7Qoqpb41YvVdV9IjIS+A1wG15r0GBJtV+qenVf20XkTuATwEJ3qpGVOLJkOA8jHchnj+2zV0QKgSqgzn1P2gFU9XUReR+Y5faPv4SU7u9zyDHGbb+FhFaqqu5zPxtF5Nd4p/BDTarpfIf2AVckvPZ5Mv97HBQ7/U+DiCwCvgx8UlVbsh1PFgznYaQD+eyPA7E70kuA36mqishYVysWEZkJnAnsUtUDQIOILHCn1LcD6ZyFDDlGF1sBcDNx11NFpFBExrjHRXgNiq0MXTrfoaeBa0WkWkSq8W4cP+3D73Fwgrojlo8LsBPvms5mt/wsS3HchHfdqB045L5YQb339cC7eHdw78n2v0nAv/denx34Nt4fWfBKG/67+568Csx06/8c2Oa+M28AN8YdsxYvSb0P/Ctu1GPQMbptVwCvJByvHHgd2OI+ww+BiM8xXuS+3814rehtca/9got9J/B5v36Pg1lsmKoxxmSQnf4bY0wGWVI1xpgMsqRqjDEZZEnVGGMyyJKqMcZkkCVVY4zJIEuqxpi0iEhERH7oyk++5QY0DFuWVDNARJrSeO39InJYRLYmrE9aY1JERrgiHJF+jlssIi+4oYfG+OlreCPCzgV+hFdOcNiypJp9DwCL4le4hPlj4DpgNnCriMx2m78APKaqXX0dVL0yahvwKgwZ4wsRKQduUtUfulV/As7IYkhZZ0k1g0Tk70Rkq1v+Jm79/3StzhdF5GER+YfYNlV9ATiWcKi+akwuJW4cs4ic71qkb7vq7Coi33ab/8Ptb4xfrgamxmYBAO6n9/d5WLFTwwxxJdw+D8zHqza+UUR+j/c7/nPgfKAIb6z36/0cLlmdyPmu4MRMVd3t3rMUeAS4XVVfFZF/xBvL/Q33uq1446aN8cs84Ouq+jMAEfkFsMVdV70HqFLVJX28Pu9YSzVzLgXWqGqzqjYBjwF/BnwMWKuqbaraCPw2jfcYA5yIe3418IaqvuqebwFq1BV0cJcIOlzdS2P8UA20gFfBCq9S1G/dWdZdWY0sSyyp5qZUNSZb8VqiMXOAt+KeX4DXEo5XArT5EKMx4FWXWuAe/y3wpKr+KYvxZJ0l1cz5L+BTIlIWu3jv1r0E3CgipSJSgVd/sj9Ja0yqNw9PxJ32g1cGbS6AiMwCPs2ptS9HA0dVtTMzH9GYXh4GLhCRnXjfxb/LcjxZZ9dUM0RV3xCRB/BqUgL8QlX/CCAij+Odmh/Ca1nWx14nIg/j1a0cIyJ7gW+o6n0icjdeEd4IcL+qbnMveQbvUsOzeF/oT7ruWEeBW1U1vmr7lcCTPnxcYwBwf+gXJK53f9DvBT4iIl9T1f8VeHBZYvVUAyAiFaraJCJlwAvAMlVNPE0f6LEuAP5WVW8bwL6PAV9V1XeH8l7GmMGzlmowlrt+pqXAg0NNqNDTIn5ORCJ99VV1lw3+wxKqMcGylqoxxmSQ3agyxpgMsqRqjDEZZEnVGGMyyJKqMcZkkCVVY4zJIEuqxhiTQZZUjTEmg/5/gB2uD68H2kwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyabc.visualization.plot_kde_matrix_highlevel(\n",
    "    h, limits=importer.get_bounds(),\n",
    "    refval=importer.get_nominal_parameters(), refval_color='grey',\n",
    "    names=importer.get_parameter_names(),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate data\n",
    "\n",
    "This section needs only be run if one wants to generate new synthetic data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Comment this line to run the code\n",
    "if False:\n",
    "    import sys\n",
    "    import pandas as pd\n",
    "    import amici\n",
    "    import amici.petab_import\n",
    "    import sys\n",
    "    import importlib\n",
    "\n",
    "    # check yaml file\n",
    "    model_name = 'cr_base'\n",
    "    yaml_file = dir_in + model_name + '.yml'\n",
    "    yaml2sbml.validate_yaml(yaml_file)\n",
    "\n",
    "    # convert to sbml\n",
    "    sbml_file = dir_out + model_name + '.xml'\n",
    "    yaml2sbml.yaml2sbml(yaml_file, sbml_file)\n",
    "\n",
    "    # convert to amici\n",
    "    amici_dir = dir_out + model_name + '_amici/'\n",
    "    sbml_importer = amici.SbmlImporter(sbml_file)\n",
    "    sbml_importer.sbml2amici(model_name, amici_dir)\n",
    "\n",
    "    # import model\n",
    "    if amici_dir not in sys.path:\n",
    "        sys.path.insert(0, os.path.abspath(amici_dir))\n",
    "    model_module = importlib.import_module(model_name)\n",
    "    model = model_module.getModel()\n",
    "    solver = model.getSolver()\n",
    "\n",
    "    # measurement times\n",
    "    n_time = 10\n",
    "    meas_times = np.linspace(0, 10, n_time)\n",
    "    model.setTimepoints(meas_times)\n",
    "\n",
    "    # simulate with nominal parameters\n",
    "    rdata = amici.runAmiciSimulation(model, solver)\n",
    "\n",
    "    # create noisy data\n",
    "    np.random.seed(2)\n",
    "    sigma = 0.02\n",
    "    obs_x2 = rdata['x'][:, 1] + sigma * np.random.randn(n_time)\n",
    "    obs_x2\n",
    "\n",
    "    # to measurement dataframe\n",
    "    df = pd.DataFrame(\n",
    "        {'observableId': 'obs_x2',\n",
    "         'simulationConditionId': 'condition1',\n",
    "         'measurement': obs_x2,\n",
    "         'time': meas_times,\n",
    "         'noiseParameters': 'sigma'\n",
    "         })\n",
    "\n",
    "    # store data\n",
    "    df.to_csv(dir_in + model_name[:-5] + '_measurement_table.tsv',\n",
    "              sep='\\t', index=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}