{ "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[petab,amici,yaml2sbml] --quiet" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os\n", "import sys\n", "\n", "import amici.petab_import\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "import pyabc\n", "import pyabc.petab\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", "\n", "import petab\n", "import yaml2sbml\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", "\n" ] } ], "source": [ "with open(yaml_file) 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[32mPEtab format check completed successfully.\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,\n", " output_dir=petab_dir,\n", " sbml_name=model_name,\n", " petab_yaml_name='cr_petab.yml',\n", " measurement_table_name=measurement_file,\n", ")\n", "\n", "# copy measurement table over\n", "_ = shutil.copyfile(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,\n", " model_output_dir=amici_dir,\n", " verbose=False,\n", " generate_sensitivity_code=False,\n", ")\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": [ ",\n", " sigma=>\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": [ "ABC.Sampler INFO: Parallelize sampling on 4 processes.\n", "ABC.History INFO: Start \n", "ABC INFO: Calibration sample t = -1.\n", "ABC.Population INFO: Recording also rejected particles: True\n", "ABC.Population INFO: Recording also rejected particles: True\n", "ABC INFO: t: 0, eps: 1.54589326e+01.\n", "ABC INFO: Accepted: 100 / 375 = 2.6667e-01, ESS: 9.9990e+01.\n", "ABC INFO: t: 1, eps: 7.72946632e+00.\n", "ABC INFO: Accepted: 100 / 439 = 2.2779e-01, ESS: 6.1330e+01.\n", "ABC INFO: t: 2, eps: 3.86473316e+00.\n", "ABC INFO: Accepted: 100 / 617 = 1.6207e-01, ESS: 9.1263e+01.\n", "ABC INFO: t: 3, eps: 1.93236658e+00.\n", "ABC INFO: Accepted: 100 / 568 = 1.7606e-01, ESS: 8.1187e+01.\n", "ABC INFO: t: 4, eps: 1.00000000e+00.\n", "ABC INFO: Accepted: 100 / 479 = 2.0877e-01, ESS: 9.0495e+01.\n", "ABC INFO: Stop: Minimum epsilon.\n", "ABC.History INFO: Done \n" ] } ], "source": [ "sampler = pyabc.MulticoreEvalParallelSampler()\n", "\n", "temperature = pyabc.Temperature()\n", "acceptor = pyabc.StochasticAcceptor()\n", "\n", "abc = pyabc.ABCSMC(\n", " model,\n", " prior,\n", " kernel,\n", " eps=temperature,\n", " acceptor=acceptor,\n", " sampler=sampler,\n", " population_size=100,\n", ")\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": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVUAAAFgCAYAAAALu+owAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABIeUlEQVR4nO3deZxddZnv+8+3dtWuMUOlKmQkEwkzMUIxqIhoAAOIURsVmkZUuqPH5hxtu6/K4bZ2ey6vo9dztO3bTukWGUwzdJQGITKDHBSQgJiBEAkhAUJCUhkqNY/P/WOtXdlV2TXutcd63q/XemXvNe2nKjtPfr/f+g0yM5xzzkWjJNcBOOdcMfGk6pxzEfKk6pxzEfKk6pxzEfKk6pxzEfKk6pxzEfKk6lwKklZI2ippm6SvpTheLunO8PizkhbkIEyXhzypOjeIpBjwA+Bi4GTgSkknDzrtWuCgmS0Gvgd8O7tRunzlSdW5o50FbDOz7WbWBdwBrBx0zkrglvD1WmC5JGUxRpenijaprlixwgDfimfLpjnAG0nv3wz3pTzHzHqAJqBu8I0krZK0XtL6U045Jde/Q9+y8H0s2qTa2NiY6xCcw8xWm1mDmTVUVlbmOhyXBUWbVAvF448/zuOPP57rMNxAu4Bjk97PDfelPEdSKTAF2J+V6Fxe86SaY4cPH+bw4cO5DsMN9BywRNJCSXHgCuDeQefcC1wTvr4ceMx8diIHlOY6gIlu5cqVuQ7BDWJmPZKuAx4EYsBNZrZZ0jeB9WZ2L/BT4DZJ24ADBInXOU+q+eJQWxexEjGpoizXoTjAzNYB6wbt+3rS6w7g49mOy+U/T6o59sgjjwCw+tVqZk+p4J+ueGeOI3LOpcOTao61t7cDsHM/dPb05Tga51y6PKnm2GWXXYaZ8eVnfk1piT83dK7Q+b/iPHC4o4fuXmNfSyf+ANm5wuYl1Rx76KGHaGrvBqCrp4/mzh4m+8Mq5wqWl1RzrLu7m9aOzv73+5o7hznbOZfvPKnm2KWXXsrsU9/d/77Rk6pzBa2gkqqkmKQ/SLov17FE6UBrV//rxpauYc50zuW7QmtT/SKwBZic60Ci8sADD7DzzUPAJAAaW7yk6lwhK5iSqqS5wKXAv+U6lqh1dPVRWRajRN6m6lyhK5ikCvwT8BVgyB7yyXNX7tu3L2uBpWPFihUcnHYS9ZPiTKsu95KqcwWuIJKqpA8Be83s+eHOS567cvr06VmKLn37W7uYVl1OfU3ck6pzBa5Q2lTfA3xY0iVABTBZ0s/N7C9yHFfa7r//fibt203ZzHfQ3dvHPn9Q5VxBK4iSqpldb2ZzzWwBwRRrjxVDQgUoKyujrRvqquNMryn3LlXOFbiCSKrF7MILL+R3HXOYVhOnflK5D1V1rsAVSvW/n5k9ATyR4zAi09LZQ1dvH3XVccCHqjpX6LykmmP3/eo+3l22g2nV5UyfVA54t6pckjRN0sOSXgn/rB3ivF5JL4bb4KVW3ATmSTXH+mJldFopddVx6muCpOrtqjn1NeBRM1sCPBq+T6XdzJaF24ezF57Ld55Uc2zGiWfyfM9c6mqSkqr3AMillcAt4etbgI/kLhRXiDyp5tj+1qBUOq063l/9976qOTXDzHaHr/cAM4Y4ryIcaPKMpI8MdbNCHJDi0lNwD6qKzWvP/4Zzy5qoq15BvLTEh6pmgaRHgJkpDt2Q/MbMTNJQXTHmm9kuSYuAxyRtNLNXB59kZquB1QANDQ3erWMC8KSaY10l5XSWVFAZjwH4UNUsMLMLhjom6W1Js8xst6RZwN4h7rEr/HO7pCeAdwJHJVU38Xj1P8cOTTmetyoX9b/3oao5dy9wTfj6GuCewSdIqpVUHr6uJxjx91LWInR5zZNqju1v7aKuJt7/fvqkch+qmlvfAi6U9ApwQfgeSQ2SEjOknQSsl/RH4HHgW2bmSdUBXv3PuSlv/4FjSkuAcwGYXlPO9n2tuQ1qAjOz/cDyFPvXA38Zvv4dcFqWQ3MFwkuqOXagN05p1ZE5t32oqnOFzZNqDpkZz3XMZNLCpf376mvi/UNVnXOFx5NqDrV19dLZ08e06iNtqrVVwetDrd25Css5lwZPqjl0oLWL95W9SvvW3/bvSyTVg23+sMq5QuRJNYcOd3RzwKqYUndklYLasNR6wJOqcwXJk2oOtXf1srFnFkuWNvTvq60Kpvw75EnVuYLkSTWHWrt6AagKR1MB/e2rB7xN1bmC5Ek1h9q7enh/fBt/eOLX/fsmV5RRIi+pOleoPKnmUGtnL/v6apg1e07/vpISMbUqzoFWT6rOFSJPqjnU1t3Lpp6ZnHnOOQP2T60q41CbV/+dK0QFk1QlVUj6vaQ/Stos6R9zHVO62ruCDv5V8YGjhad5SdW5glUwSRXoBD5gZu8AlgErJJ0z/CX5rbWzl+XxV7j3F/8xYP/Uqrj3U3WuQBVMUrVAS/i2LNwKeoB8e3cv+5jCokULB+yfVl3mSdW5AlUwSRVAUkzSiwQTBz9sZs/mOKS0tHX18GbZXM4Z1KZaWxXnYFu3T6riXAEqqKRqZr1mtgyYC5wl6dTk44W2HlBbZ++APqoJtdXBpCptYT9W51zhKKikmmBmhwgmB14xaP9qM2sws4bp06envDaftHX10tC9mTVr1gzYnxhV5U0AzhWegkmqkqZLmhq+rgQuBF7OaVBpau3qoaViOscff/yA/f2TqvioqqyT9PGwd0mfpIZhzlshaaukbZK+ls0YXX7LelKVVC3p6DrvyGYBj0vaADxH0KZ6X7TRZVd7Vy9tUxZw5plnDtifmFTFS6o5sQn4GPDkUCeE398fABcDJwNXSjo5O+G5fJfx5VQklQBXAFcBZxJ0jSqX1AjcD/zEzLaNdB8z20CwYmXRaOvqZWpY1U/m0//ljpltAZA03GlnAdvMbHt47h3ASnzxP0d2SqqPA8cB1wMzzexYMzuGYFGmZ4BvS/qLLMSRd9q6epix9/fceuutA/b3t6n6AIB8NQd4I+n9m+G+oxTaw1OXvmws/HeBmR3VOGhmB4BfAL+QdHRxbQJo6+qFY+ZyyilzB+yfUlmGBAd9qGpGSHoEmJni0A1mdtSS1Okws9XAaoCGhgbvIzcBZDypJhKqpH8ws38Y7pyJpr2rl/jMxZxxxikD9pfGSphc4QMAMsXMLkjzFruAY5Pezw33OZfVJaq/Hj61nwa8ANxhZgez+Pl5xcxo7eqhOp76r2BaddxLqvnrOWCJpIUEyfQK4M9zG5LLF9l8+m9AB/Agwf/yv5P0jix+fl7p7Omjz6Bzy2PcfPPNRx2fWlXmbao5IOmjkt4E3gXcL+nBcP9sSesAzKwHuI7gu7wFuMvMNucqZpdfsllSfdnMvhG+XivpZuDHwAeyGEPeSIyWmjxnMcuWHD1QYVpVnD2HO7Id1oRnZncDd6fY/xZwSdL7dcC6LIbmCkQ2S6qNks5IvDGzPwH5P+wpQ9rCaf9mLDiRZcuWHXV8alXcS6rOFaBsllT/G3CHpOeBjcBS4LUsfn5eaQ9LqhWlore3l1hs4HiIYKYqb1N1rtBkraRqZn8kmAf19nDX48CV2fr8fJNY9O/Vp9dx2223HXV8alWc9u5eOrp9UhXnCkk2RlTJwjnszKyTYBTV/UOdM1Ekqv/zlpzKkhk1Rx2fljRUddaUyqzG5pwbv6yMqJL0XyXNS94pKS7pA5JuAa7JQhx5JVH9X3ziySxduvSo44lRVb6sinOFJRttqiuAzwK3S1oEHAQqCRL6Q8A/mdkfshBHXklU/+Pqo7u7m7KygYPKpobj/30BQOcKSzZGVHUAPwR+GA5HrQfawzlRJ6zEon9PPvCfPFtawqc//ekBxxPVfy+pOldYsvKgStKJkr4K/G+CiVU+J+nEbHx2vkr0U132ztNpaDh62s6pPlG1y7A1a9awYMECSkpKWLBgwVGTpbvxycaDqq8SPOW/A/h9uHsuQfeqO8zsW5mOIR/1J9V3nEZ5aYolVXyiapdBa9asYdWqVbS1tQGwc+dOVq1aBcBVV12Vy9AKXjbaVK8FThk8aYqk7wKbgQmaVHuIlYi+7i46ekRFRcWA42WxEiZVlHpJ1WXEDTfc0J9QE9ra2rjhhhs8qaYpG9X/PmB2iv2zwmMTUmtnL1VlMe68807uuOOOlOcEq6p6UnXRe/3118e0341eNkqqXwIelfQKRyb2nQcsBv5rFj4/L7V39VJVHuPss88e8pxan6nKZci8efPYuXNnyv0uPdl4+v+ApOMJlqBIzI6+C3jOzCbscKG27l6q4qWcdNJJQ55TW1XG/hYvqbro3XjjjQPaVAGqqqq48cYbcxhVccjK038z6zOzZ8zsF+H2jJn1SvrMaO8h6VhJj0t6KVzt8ouZjDnT2jp7qCyL0dbWdlTbVsI0r/67iCSe9EuitLSUv/iLv6CyspK6ujokMX/+fFavXu3tqRHI5oQqqfwj8LNRntsD/K2ZvSBpEvC8pIfNrCAXW2vr6qW6PMZdd90FcFQ/VfCZqlw01qxZw2c/+1m6uoLvUm9vUEHcv38/VVVV3HbbbZ5MI5TxkqqkDUNsG4EZo72Pme02sxfC180EkwOnXGytELR191IZL+Vd73oX73rXu1KeM626jNauXjp7JmwriYvAF7/4xf6EOlhbWxvXXHON91WNUDaq/zOATwGXpdj2j+eGkhYQLFf9bDQhZl9bZw9VZTFOOOEETjjhhJTn+FDV7JP08bB5qU/S0aMyjpy3Q9JGSS9KWp/NGMdq//7h/5n19vZiZv19Vb/whS/4oIA0ZKP6fx9QY2YvDj4g6Ymx3kxSDcEqrF8ys8ODjq0CVkH+P8VsC5/+t7S0AFBTM/RMVQdau5gxueKo4y4jNgEfA34yinPfb2aNGY4nq9ra2vjRj37U/94HBYxdxkuqZnatmT01xLExLZYWzh3wC2CNmf0yxf1Wm1mDmTVMn57fiwq0dfVQFY+xdu1a1q5dm/IcH6qafWa2xcy25jqOqERRymxra+Nzn/ucl15HKdcPqkZNkoCfAlvM7Lu5jiddbV1Bl6pzzz13yHP651T1oar5yICHJBnwEzNbneqkXNeebrjhhkju09raSmtrK+Cl15FkLalK+nKK3U3A86maBlJ4D3A1sFFS4vz/Hi7AVlB6+4zOnj6q4jEWL1485HnTqo5MVO2iI+kRYGaKQzeY2T2jvM25ZrZL0jHAw5JeNrMnB58UJtvVAA0NDVmfiD1TI6R8SOvQsrnwXwPweYIn9nOAzxHMtfqvkr4y0sVm9pSZycyWmtmycCu4hApHZv2visdoamqiqakp5XlT+ydV8aQaJTO7wMxOTbGNNqFiZrvCP/cSrL56VqbiTUcmS8c7d+70JoEUsplU5wKnm9nfmtnfAmcAxwDnAZ/OYhw5l5j1vzJeyt13383ddx+1IjIA8dISaspLfahqnpFUHfaVRlI1cBHBA668c+ONN1JVVZWRe0ti586dA3oOeGLNblI9BuhMet8NzDCz9kH7i15i2r/qeIzzzjuP8847b8hzp1aVefU/iyR9VNKbwLuA+yU9GO6fLSlRM5oBPCXpjwTTWd5vZg/kJuLhXXXVVVxzTWZWKxq8rFyiSWCiy+aDqjXAs5LuAQR8CPj38H/6ghwVNV6tSdX/RYsWDXvutGofqppNZnY3QXV+8P63gEvC19uBd2Q5tHFZs2YNt9xyS9Y+z2e5ymJSNbP/IenXBA+cAD5vZolO0xOqtTtR/a+Kl3Lw4EEAamtrU57rQ1VdOlLNm5pJ+d4/PBuy3aWqm2AOVQtfT0it/Uk1xj33BM9GUo39B5hWVcaOxtZsheaKTJQlx9NOO43ly5czZcoUmpqaePTRR9m4cWP/cZ/lKpC1NtVwVqk1BAv/HQP8XNKEnE81sehfZTzG+eefz/nnnz/kuV5SdemIquR42mmncdlllzF16lQkMXXqVC677DJOO+00n+VqkGyWVK8FzjazVgBJ3waeBv6/LMaQF448qCplwewFw547rTpOc2cP3b19lMWy+VzRFYNU86aOx/Lly4nH4wP2xeNxPvjBD7Jhw4a07l1ssvmvVEDydEu94b4JJ7n639jYSGPj0MPHa32oqkvDVVddxerVKQd7jcmUKVNS7k81Z8VEl82k+jOCp///IOkfCWaYuimLn583kqv/9913H/fdd9+Q59ZW+0xVLj1RVMmHGqAyVLKdyLL59P+74axUiaf/14xyeGrRaUt6+r98+fJhz00MVT3g7aouDSUlJfT1jX+dzUcffZTLLrtsQBNAT0/PiN/fiSjjSVVSM8HT/v5dScfMzCZnOoZ8c6itm0nlpcRKxLHHHjvsuT5U1UXhc5/73IAp/cYq8ZQ/+en/8ccfz2mnnRZViEUjGwv/Tcr0ZxSafS2d1E8qB2Dv3r0AHHPMMSnP7Z+pyqv/Lg0//OEPAVi9enX/ciqjIYmysjK6urrYuHEjGzduRBKf//zn+S//5b9kKtyC5o+Tc6CxuZP6miBZrlu3jnXrhp4XxudUdVH54Q9/SE9PD/Pnzx/V+fF4nNtuu42bbrqJ+fPn93eduu222/qTtDtawcynWkwaWzo5fkZQgL/wwguHPbeiLEZVPObVfxeZ0QwIqKur4/vf/37/Qy7vfzp6nlRzoLGli3cfF1T/58wZee3C2qo4B7yk6iIyb948du7cOeTx+fPns2PHjuwFVGS8+p9lXT19NLV3U18TJNU9e/awZ8+eYa+prS7zLlUuMsNNB+hDTdPnSTXL9rcGsxzWTwraVB944AEeeGD4WeOm15Tz9uGOjMfmJobEgIBE22osFgPwoaYR8ep/ljU2B9X4REl1xYoVI14zt7aKP7xxKJNhuQnmqquu8uSZIZ5Us6yxJSyphkl15sxUSyUNNLe2kkNt3TR3dDOpoiyj8Tnn0lMw1X9JN0naKykvl60YrX1hUp0eJtVdu3axa9euYa+ZWxu0f715sD2zwTkkfUfSy5I2SLpb0tQhzlshaaukbZK+luUwXR4rmKQK3EywUGBB6y+phm2qDz/8MA8//PCw18ytrQQ8qWbJw8CpZrYU+BNw/eATJMWAHwAXAycDV0o6OatRurxVMNV/M3tS0oJcx5GuxuYuquMxquLBr/6SSy4Z8ZojSTV7M7hPVGb2UNLbZ4DLU5x2FrAtXFYFSXcAK5lgywK51AomqRaLxqQhqjD08NRk06rjVJbFvKSafZ8F7kyxfw7wRtL7N4GzU91A0ipgFfhSIxNFIVX/RyRplaT1ktbv27cv1+Gk1NjS2f+QCuCNN97gjTfeGOaKYPz13NpKL6lGRNIjkjal2FYmnXMD0EOwWsW4mdlqM2sws4bp06enG7orAEVVUjWz1cBqgIaGBhvh9JxobOlkYX11//tHH30UGHqNqoQgqXpJNQpmdsFwxyV9mmC13+U2eB3mwC4geXqxueE+54orqRaCxpYuzlwwrf/9hz70oVFdN7e2ihdeP5ShqFyCpBXAV4D3mdlQVYPngCWSFhIk0yuAP89SiC7PFUz1X9LtBGtanSDpTUnX5jqmserp7eNgW9eA6n99fT319fUjXju3tpKm9m4Od/hw1Qz7F2AS8LCkFyX9GEDSbEnrAMysB7gOeBDYAtxlZptzFbDLLwVTUjWzK3MdQ7oOtHZhxoAHVYmJKxYsWDDstcdOC/qq7jrYzuRZPgAgU8xs8RD73wIuSXq/Dhh6zkY3YRVMSbUYHOn4f2RJiieeeIInnnhixGu9r6pzhaFgSqrFYF/zwCGqACtXrhzq9AGOjKryHgDO5TNPqlnU2DJwMhWA2traUV1bW1VGVTzGGwe8pOpcPvPqfxYdGaJ6JKlu376d7du3j3it91V1rjB4STWLGps7qSgroToe69/35JNPArBo0aIRr59bW+Vtqs7lOU+qWZQYTSX1r9LNRz/60VFfP7e2kvU7DmQiNOdcRDypZlFjy8A+qgBTpkwZ9fVzays53NFDU1s3U6q8W5Vz+cjbVLNo8Lh/gG3btrFt27ZRXX/yrCABv/D6wchjc85Fw5NqlnT19PFaYyvz6wYuuPbUU0/x1FNPjeoeDQtqicdK+N2rjZkI0TkXAa/+Z8mmt5ro7OmjYf7ALlSXX55qus7UKspinD5/Kr/dtj/q8JxzEfGSapY8vyOosp+xYGBSrampoaamZtT3ec9x9by0+zAHW7sijc85Fw1PqlmyfucB5k2r4phJFQP2b926la1bt476Pu9eXAfA09u9tOpcPvKkmgVmxvodB2lYcPToqaeffpqnn3561PdaOncq1fEYv93m7arO5SNvU82CHfvb2N/aRcP8aUcd+8QnPjGme5XFSjh7UR1PvxqUVNdt3E2fGR9aOjuSWJ1z6fGkmgWJDvupSqpVVVVH7RvJu4+r47GX9/KFNc+zbuMeYiVi/rRqTps7+j6vzrnM8Op/Fjy/8yCTK0pZPP3oB1Jbtmxhy5YtY7rfu48LJrX+9aY9fP59x1FXHef/WvtHunr6IonXOTd+XlLNgud2HOCM+bWUlOioY88++ywAJ5100qjvd+LMSXzpgiWcuWAa71lczxnza/mrW9fzwye28aULjo8s7olI0neAy4Au4FXgM2Z2KMV5O4BmoBfoMbOGLIbp8pgn1Qw72NrFq/ta+djpc1Mev+KKK8Z8z5ISDUieF548g5XLZvMvj23jsnfM5rgUJWI3ag8D15tZj6RvA9cDXx3i3PebmT8xdAN49T/Dbn16JwDvXZJ6HaqKigoqKipSHhuLv//QyZSXlvC/Hhx99yx3NDN7KFyDCuAZgpVSnRs1T6oZtOtQOz/6zTYuPW0WS+dOTXnOpk2b2LRpU9qfVV9Tzl+dt4hfb9rDH3xugKh8Fvj1EMcMeEjS85JWDXUDSaskrZe0ft++fRkJ0uWXgkqqklZI2ippm6Sv5TqekfzPdVswg+svOXHIc9avX8/69esj+by/fO8i6qrjfOvXL5N6uXoHIOkRSZtSbCuTzrkB6AHWDHGbc83sdOBi4K8lnZfqJDNbbWYNZtYwffr0yH8Wl38Kpk1VUgz4AXAh8CbwnKR7zeylVOf3mdHW1ZPqUEaYBUWXPjPeONDGs9sPcN+G3XzpgiX960ulctVVV0UWQ015Kf9t+RK+ce9mvvXAy1xw0gyOnzGJWIkQoKOfkxWMqnh0X1Uzu2C445I+DXwIWG5D/O9kZrvCP/dKuhs4C3gysiBdwSqYpErwpd1mZtsBJN0BrARSJtXNbx3m5K8/mMXwjnb8jBo+d95xw55TVhbtvKhXnjWPR7a8zU9+s52f/GbkZVoKxY5vXZqVz5G0AvgK8D4zS7l2jaRqoMTMmsPXFwHfzEqALu8VUlKdA7yR9P5N4OzkE8K2rVUAdXMWcv3FQ1e7M0ECIWZOqeCkWZNZUFdFaWz4FpYNGzYAsHTp0khiiJeWcNu1Z3OwtYtnXzvAGwfaMAxvDRi1fwHKgYfDFRqeMbPPS5oN/JuZXQLMAO4Oj5cC/25mD+QqYJdfCimpjsjMVgOrARoaGuxz7xu+lJgPXnjhBSC6pJpQWx1nxakzI73nRGBmi4fY/xZwSfh6O/CObMblCkchJdVdwLFJ7+eG+wra1VdfnesQnHMRKqSn/88BSyQtlBQHrgDuzXFMaYvFYsRisZFPdM4VhIIpqYYjXK4DHgRiwE1mtjnHYaXtxRdfBGDZsmU5jcM5F42CSaoAZrYOWJfrOKLkSdW54qJi7SQuaR+wMwcfXQ8UynjwQoq1wsxOzXUQ6ZDUDOT7OOJC+U7kOs5GM1uR6kBBlVTHwsxyMnxF0vpCmbGo0GLNdQwR2Jrvv+9C+U7kc5yF9KDKOefynidV55yLkCfV6K3OdQBj4LFmVyH8DIUQI+RxnEX7oMo553LBS6rOORchT6rOORchT6oZIOk7kl6WtEHS3ZKm5jqmoUj6uKTNkvok5V0XlXyemHyk2CSVS7ozPP6spAXh/gWS2iW9GG4/TrrmDEkbw2v+WUpvFtw0YrwqKb4Xw+/HsvDYE+E9E8eOyXCM50l6QVKPpMsHHbtG0ivhdk3S/kh/j2NiZr5FvBHMr1kavv428O1cxzRMrCcBJwBPAA25jmdQbDGCFU0XAXHgj8DJuY5rtLEBXwB+HL6+ArgzfL0A2DTEfX8PnAOIYCmXi3MR46BzTgNeTXof2XdllDEuAJYCtwKXJ+2fBmwP/6wNX9dG/Xsc6+Yl1QywAlo8zsy2mFm+jvLpn5jczLqAxMTk+WA0sa0EbglfrwWWD1dikjQLmGxmz1iQGW4FPpIHMV4ZXpsJI8ZoZjvMbAPQN+jaDwIPm9kBMztIsBLuigz8HsfEk2rmDbd4nBteqonJ5+QolsFGE1v/OeF/sk1AXXhsoaQ/SPqNpPcmnf/mCPfMZowJnwRuH7TvZ2HV/+/TrFqn83c81LVR/x7HpGiHqWaapEeAVLNA32Bm94TnjLR4XFaMJlaXVbuBeWa2X9IZwH9KOiXXQaUi6WygzcySl/y9ysx2SZoE/AK4mqA06PCkOm4WweJx2TJSrHksnycmH01siXPelFQKTAH2h9+HTgAze17Sq8Dx4fnJTUXp/rzjjjHp+BUMKqXakUUPmyX9O0EVfrxJNZ2/413A+YOufYLof49j4tX/DEhaPO7DNsTicW5U8nli8tHEdi+QeCJ9OfCYmZmk6QpWB0bSImAJsN3MdgOHJZ0TVqk/BaRTkxh3jGFsJcAnSGpPlVQqqT58XUZQcNjE+KXzd/wgcJGkWkm1BA+IH8zA73FssvVEbCJtwDaCtp4Xw+3HuY5pmFg/StDm1Am8HX4pcx5XUnyXAH8ieEJ8Q67jGSk2glVVPxy+rgD+I/w+/B5YFO7/M2Bz+N14Abgs6Z4NBEnqVYJFCJWLGMNj5xMsfJh8v2rgeWBD+DN8H4hlOMYzw+9oK0EpenPStZ8NY98GfCZTv8exbD5M1TnnIuTVf+eci5AnVeeci5AnVeeci5AnVeeci5AnVeeci5AnVeeci5AnVeeci5AnVeeci5AnVeeci5AnVeeci5AnVeeci5AnVeeci1BWk+p4F/iStEzS0woWqNsg6ZPZjNs550Yra7NUhfNH/gm4kGAar+eAK83spaRzFgCTgb8D7jWzteH+4wEzs1ckzSaYeuwkMzs01OeVxo63qvhnM/TTuMEOt38t5ZIaF1xUZfv396a85sUXuh40sxUZDSyPrFixwh544IFch+GiMeQSMtmc+b9/gS8ASYkFvvqTqpntCI8NWODLzP6U9PotSXuB6cChoT7M54bOD/sb+/jNU/NTHptS/Up9lsPJqcbGxlyH4LIgm0k11SJdZ4/1JpLOIljK9tUUx1YBqwDE1HEF6aIlQL3edF8sHn/8cQDe//735ziS/FVQa1SFS8/eBlxjZoOXq8XMVgOrAWIlc3327XxgoKP+plyhOnz4cK5DyHvZTKppLeImaTJwP8FyC89EHJvLFAOlblJ1BWjlypW5DiHvZbNeNu4FvsLz7wZuTTy8coVDPZZyc7n31qF2Wjt7ch1GUclaUjWzHuA6ghUQtwB3mdlmSd+U9GEASWdKehP4OPATSZvDyz8BnAd8WtKL4bYsW7G7NIQl1VSby72P//hpPv/z5xltL6BHHnmERx55JMNRFbastqma2Tpg3aB9X096/RwD1+tO7P858POMB+giJ8NLpXmquaObXYfa2XWonV9v2sMlp80a8Zr29vYsRFbYCupBlStMXirNTzv3B90O46Ul/D/3vcT5J0ynKj58SrjsssuyEVpB874uLrMM1GspN5dbO/a3AvD3l57EW00d/Mtj23IcUXHwpOoyy0A9qTeXW4mS6p+dMZfL3jGbm3+3g76+4f+ze+ihh3jooYeyEV7B8qTqMs4fVOWnHY2tHDOpnKp4KWctnEZbVy97mzuHvaa7u5vu7u4sRViYvE3VZZQM1DvkMGmXQzv2t7KgrhqABXVV/ftmTqkY8ppLL700K7EVMi+puszrG2JzObVjfxsL6oNkOn9akFx3hu2sbvy8pOoyywCv6ued1s4e9jV3Mj8sqc6eWkFpifrbWYeSmGVrxYoJM7nYmHlJ1WVWWP1PtY2GpJikP0i6L3y/UNKz4Zy8d4aj7ZBUHr7fFh5fkLkfqvAlnvwnqv+lsRLm1laOmFTdyDypuszrUeptdL5IMAIv4dvA98xsMXAQuDbcfy1wMNz/vfA8N4RE8pwftqUGr6vZeWD46v+KFSu8lDoCT6ous0zQO8Q2AklzgUuBfwvfC/gAkJj/4RbgI+HrleF7wuPLw/NdCv0l1frq/n3z66rY2dg26iGrLjVPqi6zhq/+10tan7StGnT1PwFf4chjrTrgUDiPBARz8s4JX/fP1xsebwrPdynsaGylvqacmvIjj1Xm11XT3NnDwbahu0zdf//93H///dkIsWD5gyqXeUOXShvNrCHVAUkfAvaa2fOSzs9QZBPWjv1t/d2oEuZPO9Ktalp1POV1ZWVlGY+t0HlSdZllwPhm/n8P8GFJlwAVBGuXfR+YKqk0LI0mz8mbmK/3TUmlwBRgf5rRF62d+1s5d/H0AfsS3ate39/G6fNqU1530UUXZTy2QufVf5dh42tTNbPrzWyumS0gmHv3MTO7CngcSKy0ew1wT/j63vA94fHHbAyNg5JukrRX0qakfdMkPSzplfDP2nC/JP1z2NNgg6TTR/s5+aCtq4e3D3eysH5gSXVubRXSkfZWNz6eVF1mGVhvScptnL4KfFnSNoI205+G+38K1IX7vwwctQT6CG4GBj/W/hrwqJktAR5NuufFwJJwWwX8aIyflVNHnvxXD9hfURZj1uQKXh+mW9WvfvUrfvWrX2U0vkLn1X+XWcaonvQPewuzJ4AnwtfbCVbmHXxOB8Hk5uP9jCdT9G1dCZwfvr4ljOGr4f5bw5LwM5KmSpplZrvH+/nZtOdwBxB0+B9sXl3VsCXVysrKjMVVLDypuszrK9ieTTOSEuUeYEb4OtXKwHOAo5Jq8gq/8+bNy1ykY3C4PXi6P6Xy6IdOC+qqeWTL20Nee8EFF2QsrmLh1X+XWabgQVWqrYCEpdIxd+A0s9Vm1mBmDdOnTx/5gixIJNXJKZLqvLoqGlu6aPF1q8atsL7ZriBZr1JuBeDtcFn0xPLoe8P9aa0MnGtNiaRakbqkCkNPrHLPPfdwzz33pDzmAp5UXWYlulQVZkk1uUfB4J4Gnwp7AZwDNBVKeyoESbW8tISKsthRx+ZNO9KtKpXJkyczefLkjMZX6LxN1WWYCiKBSrqd4KFUfbii7zeAbwF3SboW2Emwqi8Ei1deAmwD2oDPZD3gNBxu70nZngowtzZ4ELXrUOoF/t7//vdnLK5i4UnVZZQZWAE8qDKzK4c4tDzFuQb8dWYjypym9u4hk+qUyjKq47Ehk6obmSdVl1lhP1WXP5rau1M+pAKQxOyplew6mDqp/vKXvwTgYx/7WMbiK3SeVF2GFUb1fyI53NHNjMlDL5kyp7ZyyJJqXZ3PUTMST6ous4xCedI/YTS1d3P8jElDHp8ztZI/vnEo5bH3ve99GYqqeHhSdRnn1f/8MlybKgQl1YNt3bR19VAV9xQxVv5td5llGnpzWdfXZ7R09jC5YuhkOWdq2AMgRbvq2rVrWbt27VH73RH+35DLOC+p5o/mjh7MUo+mSkgk1TcPtbNkUDPBzJkzMxpfMcjqt13SCklbwynTjppFSNJ5kl6Q1CPp8kHHrgmnYHtF0jWDr3X5KdGlKtXmsq9pmHH/CXPCvqpvpXhYde6553LuuedmJrgikbWSqqQY8APgQoIJKJ6TdK+ZvZR02uvAp4G/G3TtNILO2A0EY3SeD689mI3YXTrkJdU8crhj6HH/CcdMCparHqpblRteNr/tZwHbzGy7mXUBdxBModbPzHaY2QaOrEmU8EHgYTM7ECbShzl67kuXj6KfT9WlYTQl1ViJmDmlImW3qrvuuou77rorY/EVg2y2qaaaLu3sNK6dM/ik5GnWxNRxBemiZ32eQPPFaJIqBO2qqar/c+fOzUhcxaSoHlSZ2WpgNUCsZK6vs5sPrGBmpJoQhpv2L9mc2kqeefXoJb7e/e53ZySuYpLNIkQ606UV9FRrE5nh1f98MpaS6p7DHXT3Dm6JcyPJ5jf7OWCJpIWS4gSLud07ymsfBC6SVBsuvnZRuM/lOwuq/6k2l32HO7qJlYjq+NHT/iWbM7WSPoM9TR0D9t9+++3cfvvtmQyx4GWt+m9mPZKuI0iGMeAmM9ss6ZvAejO7V9KZwN1ALXCZpH80s1PM7ICk/0GQmAG+aWYHshW7S4c//c8nTe3dTK4oRRq+SSa5W9Wx046surpw4cKMxlcMstqmambrCOaiTN739aTXzxFU7VNdexNwU0YDdJEzg75x9kmVVAE8CZQTfFfXmtk3JC0k6D1SBzwPXG1mXZLKgVuBM4D9wCfNbEf6P0XxaBpmLtVk/aOqBj2sOuecczISVzHxIoTLvPEPU+0EPmBm7wCWASvCmfa/DXzPzBYDB4Frw/OvBQ6G+78XnueSHB5h3H/C7GGGqrrheVJ1mWWir7ck5TbipYGW8G1ZuBnwASAxAP0W4CPh65Xhe8LjyzVSPXeCGW4u1WQVZTHqa+JHlVTXrFnDmjVrMhVeUSiqLlUuPw3zUKpe0vqk96vDbnH9wpF4zwOLCUbkvQocMrPEcp/JfZb7+zOHbfhNBE0EjeONXdIJwJ1JuxYBXwemAn8F7Av3//eweSuvHW7v7m8vHcmMyRXsbe4csO/444/PRFhFxZOqy6zh51NtNLOGYS836wWWSZpK8BDzxGgDHJ6ZbSVoekgk+F1hHJ8haIL4X9mMJ12HO0ZX/QeorymnsWVgUj3zzDMzEVZR8eq/yyhD9PWVpNzGdB+zQ8DjwLuAqZISBYLkPsv9/ZnD41MIHlhFZTnwqpntjPCeWWNmI86lmqyuJs7+lq4MR1V8PKm6zEpj7L+k6WEJFUmVBJPxbCFIrolZzAYvHZ2Ywexy4LFwkb6oXAEkd9K8TtIGSTeF/adT/QyrJK2XtH7fvn2pTsma9u5eunuNyRVjK6km/wpvvfVWbr311kyFWBQ8qbqMS2Pqv1nA45I2EPRRftjM7gO+CnxZ0jaCNtOfhuf/FKgL938ZOGp6yfEKB6x8GPiPcNePgOMImgZ2A/871XVmttrMGsysYfr06VGFMy6H24Nm6NFX/+N09vTR0tnTv++UU07hlFNOyUh8xcLbVF1GJar/47o2mLHsnSn2byeY9Wzw/g7g4+P6sJFdDLxgZm+Hn/V24oCkfwXuy9DnRma0Q1QT6qrLAdjf0sWksHR7xhlnZCa4IuIlVZdZxTP135UkVf0lzUo69lFgU9YjGqOm/slURleWqp8UJNXBD6vc8Lyk6jJs/CXVfCGpmqA993NJu/9fScsI+s3uGHQsLx0ec0k1DkBj0sOqm2++GYBPf/rTkcZWTDypusyyYKhqITOzVoK22+R9V+conHEba/V/eoqS6rJlyyKPq9h4UnUZlZj6z+XeWJNqbVVQUk3uVuVJdWSeVF1mGQVf/S8WifWpJo2yS1W8tIQplWUDSqq9vb0AxGLDTx04kXlSdRlW+G2qxaKpvZtJ5aXESkY/HUJ9TZz9rUeS6m233QZ4m+pwPKm6jLLhh6m6LGru6GFSxdj+ydfVlA94UHX66adHHVbR8aTqMs5LqvmhpaNn1FX/hOk15WzZc7j//dKlS6MOq+h4UnUZ5tX/fNHS2UN1+djaQgeP/+/uDtply8rGlpwnEv+2u4wyX6MqbzR39lAzxpJqfU05Te3ddPUECwD6fKoj85KqyyyDPm9TzQutnT3MnTq6uVQT6mqCblUHWruYOaWChoZhZ2p0RFBSlfTVKAJxxSqaqf9c+lo6eqgpH+ODquqBAwBOPfVUTj311MhjKyZjLqlKuiv5LcEsPb4W0DDiNjCBdGnirKVu+IOqfBG0qY7tn/z0SYmhqkFS7egIlqyuqKiINrgiMp7q/2Ez+8vEG0k/ijAeV2wMbHSL/LkM6uszWjp7qBlrl6r+kmrwsOqOO+4AvJ/qcMaTVG8c9P6GKAJxxUo+TDUPtHUHI6EmjbGkmpipan9YUj377LOjDawIjfgblrQA+GuCCXkPAC9K+lViSQkzO5DRCAtYLCyh1VpQVUo0AxxU8AVtKenOTWBZZD5MNS+0dAQTTY+1pFodj1FeWsL+1qCketJJJ0UeW7EZzbf9HuBlgpUsLwTeATwp6QeSyjMZnCsOfX1KubnsaekM/gMfa5uqpGBZlXBV1ba2Ntra2iKPr5iM5jccM7OfAkg6YGZ/FS6q9jfAao6sCeSSxExMsqCRv7Yv+L+nhqCPYB/BXHjt4SrLvSrwufGG5Z3/80FzWFIda/UfgvH/jWFJ9a67gufU3qY6tNH8hh+RdJ2Z/QvBw1zCNde/I+lPGY3OFTwz6PM21Zxr7QzaVMda/Ydg/P+epuCp/7ve9a5I4ypGo/kNfxm4XtJ6YLakVUAbwVLBUS7/64qSl1TzQaL6P9Z+qhCUVDe/1QTACSecEGlcxWjEb7uZ9ZnZjcB5wCpgJnAGwZo8F2c2vMITtxLiVsLsvmoW9U1mUd9kjrFKjrFKKixGhcWopIxKyohTQrzIRwonHlSNp/O/pGMlPS7pJUmbJX0x3D9N0sOSXgn/rA33S9I/S9oWLh0dyZRKknZI2ijpxbBwMWQM+SpR/R9PUq2rKWd/S1fQLaulhZaWlqjDKyqj/hdtZm1mdq+ZfdPM/sbMfmRmh8byYZJWSNoafumPWj5YUrmkO8Pjz4Y9D5BUJumW8Iu9RdL1Y/ncXDnmlFbOuW43Df/3ayz9r69Td8rE/DKmMaKqB/hbMzsZOAf4a0knEyw9/aiZLQEe5chS1BcDS8JtFcEy0lF5v5ktM7PEOM2hYshLiWWmx1dSLaenz2hq72bt2rWsXbs26vCKStbG/kuKcaQHwZvAc5LuNbOXkk67FjhoZoslXUEwUuuTBMsOl5vZaZKqgJck3W5mO7IV/1gde3IHJ17aRiwePIQqn9rLwg8FrSV7tkygWdNN425TNbPdwO7wdbOkLcAcYCVwfnjaLcATwFfD/beamQHPSJoqaVZ4n6gNFUNeag2T6lif/sPAtarOPffcSOMqRtmse54FbDOz7WbWBdxB8MVMtpLgCwqwFlguSQQPyKrDXgeVQBdwmDxUZWVUWRmnv7+zP6EmxOLGwg8cYn5vNfN7q5nWV8G0vgpipv4+rcXGCEZUpdrGIqy1vBN4FpiRlCj3ADPC13OAN5IuezPcly4DHpL0fPhMgWFiGBz3KknrJa3ft29fBKGMT3NnD/HSEuKlY/8nXx9OqrKvuZPFixezePHiqMMrKtmcpSrVF37w8Iz+c8ysR1ITwSqWawkS7m6gCvibVIMOwi/8KgAxNeLwxyY+pSfl/pIh9hez3qH7pNYn2ihDq81s9eCTJNUAvwC+ZGaHg/9nA2ZmUsb7pJ1rZrskHQM8LOnl5IPDxRD+PKsBGhoactZ3rqWjZ1zdqSCYqBpgX0snTU3BA6spU6ZEFluxKZSp/84CeoHZQC3wfyQ9Ymbbk09K/gLHSubm5AtcFhb+e5tKKZ16dAItbxNnlwW/9uaeGgAOlATdVdrpzVKU2WM27NP/xqQ2ypQklREk1DVm9stw99uJar2kWcDecP8u4Niky+eG+9JiZrvCP/dKupvg+zhUDHlpPOP+E45U/7u4++67Ae+nOpxsVv9H84XvPyes6k8h6Lb158ADZtZtZnuB3wJ5PbFj86P1lPQOzOvqgdmbJt4gtDSe/gv4KbDFzL6bdOhejgw6uYZg1F9i/6fCXgDnAE3ptqdKqpY0KfEauIig58tQMeSl1s6xT/uXMKWyjLKY2NfcyXnnncd5550XcXTFJZtJ9TlgiaSFkuLAFQRfzGTJX9TLgcfChw6vAx+A/i/2OQRDZ/NW+8YpnLzVqOgwMKOqq5f5z1dQ90Y816FlV9j5P9U2Cu8BrgY+EHZnelHSJcC3gAslvQJcEL4HWAdsB7YB/wp8IYKfYAbwlKQ/Ar8H7jezB4aJIS81d4x92r+E/qGqLZ0sWrSIRYsWRRxdccla9T9sI70OeBCIATeZ2WZJ3wTWm9m9BKWS2yRtI5i85Yrw8h8AP5O0mWAO15+Z2YZsxT4aiclSaiwYijq9RJxTt4NzmmDGu4KBZ/sOL4bjYNLkoAP1679bAsDOWFB6bVfxjam2NDr/m9lTBH/fqSxPcb4RTP4TmbCJ6R0p9u9PFUO+aunsYebk8c+BWl9Tzr7mTg4ePAhAbW1ed8vNqay2qZrZOoLSRPK+rye97iDoPjX4upZU+10BMOj1EVU5l06bKgQ9APY2d3LPPUErh7epDq1QHlS5AmX4GlX5IJ02VQgeVr20+zDnf/j86IIqUp5U05ToXzrFgir8nL5qAI6b0cqkmUFVqeTEoPdXXWvQo2xuU3DOseuPA6DKgr+Gmr4ySsLaboeCngMFv/SKyUuqeaB5HOtTJQvaVLuYN28+JSX+n+RwPKm6jPI1qnKvq6ePzp6+tEuqvX3Ga7v2MKWyjPr6+ggjLC6eVNMUC0uWdX3BQ4CFvcGfs2a+TtW8YFhq2+LgnMrtwYOonq7g157ow17eFwxbraOCGX1VAOwraQfg9ZJmoJDnXPWSaq4lhqim16Ya1MQe/PU6quIxb1MdhidVl1EW9ChzOZTOZCoJiaQ6/7SzWDp3ahRhFS1Pqi7jfJLq3IoiqSZGVfVUTuPYY6OYTqF4eVJNU6L6X04s/DPQ11eCSoOHTL3VwVIUfS3hcL/dQXtUV2/iIVcwIKCMEmaFzQiJRwH7FDQDtKgwFwk0G3bsv8uClgiq/4nx/2/t3sPe2WUcc8wxkcRWjLwI4TIsaFNNtbnsaEljguqEyZWlxGMl7NvyLOvWrRv5ggnMS6ouo4Kn/15SzaVESXVSGiXVYKhqnOa6k7jwguOjCq0oeVJNUxdBFb8xrKa/WhJUkxbumMHBjcH8MXWzg6lfD785DYCO9qCKP6kimJXquJbgif/m0mZ2xcLqPkF1v7fg+6lCr3f+z6mWNCaoTlY/qZy3e+PMmeNtqsPxOpjLKPPqf85FUf2HoF217WAje/bsiSKsouUl1TQl+o8m+pX2Bat4U324nvn3vQeA47YF/7PvfWs6AK/umAlAZ3eQWGrDglwJ4rWSoFTbGc6tmigJFyp/UJV7zYmSajzNkmpNOdU7t/LAA3u9n+owvLjgMq6vTyk3lx2tnT1Ux2NpDy+dPqmc33Udy4UXfTCiyIqTl1RdRgUl1VxHMbG1dKQ3Q1VCfU2cxt5KyidPiyCq4uVJNSKJiU/2hs0AfyjdT8fLdQDMeymo7h8MzzkUTpYSC6v4veH7Q+qkWV0D7lcM+op0UcNC0ZLmDFUJ9ZPKqVcrW7fvpH6p9wAYilf/Xcb19aXeXHY0d/ZQU1GW9n2m15TTUPYG63/7mwiiKl6eVF1GGdBrSrkVAknHSnpc0kuSNkv6Yrj/HyTtGrTMS14K5lKNpX2f6ZPKeaZ7PvUnDl4E2SXz6n/EEr0B9pS00RpW67eGS6wkqvSJJ/uxQSuFtKm7qKr9EKymWigJdAg9wN+a2QvhAoDPS3o4PPY9M/tfOYxtVFo6eqivqUr7PrOnVnLIKmnsrYwgquLlSdVlXG8Bz1IVrsa6O3zdLGkLUFC934M21fSr/xVlMU6s6eLNN98AlqQfWJHy6r/LqKD6n3obiaSbJO2VtClp3zRJD0t6JfyzNtwvSf8saZukDZJOj/pnkbQAeCfwbLjruvCzbkrEkeKaVZLWS1q/b9++qEMaleaO7rSGqCZbWvIGemvTyCdOYJ5UM6RXxqGSTg6VdPJWSStvlbTSqHYa1d6/f39Jx4CtXb25DjsjxptUgZuBFYP2fQ141MyWAI+G7wEuJig+LQFWAT+KIvYESTXAL4Avmdnh8P7HAcsISrL/O9V1ZrbazBrMrGH69OlRhjQqXT19HO7oYWpV+iVVAJt3OuttQST3KlaeVF1GGdA7xDbitWZPEixVnmwlcEv4+hbgI0n7b7XAM8BUSbPSiz4gqYwgoa4xs1+Gsb1tZr1m1gf8K3BWFJ8VtbcPdwAwe0o07aDzZs1ge3OMju7iLABEwZNqFvTKBmwTSTrV/yHMCNs5AfYAM8LXc4A3ks57kwjaPiUJ+Cmwxcy+m7Q/OWF/FMjLOvHupiCpzpxSEcn9plkTM0qaeeNAWyT3K0b+oMpl3DBlmnpJ65Perzaz1aO9r5mZlPH/pd4DXA1slPRiuO+/A1dKWkbw/8YO4HMZjmNcdjcFg1FmRZRUD776Iu8sbWLn/jaWzJgUyT2LjSdVl1FG0CdpCI1m1jDGW74taZaZ7Q5Li3vD/buAY5POmxvuS4uZPQWk6hNWEDM174m4pHrJpZfx4x/8lrO8pDokr/67jLMhtnG6F7gmfH0NcE/S/k+FvQDOAZqSmgkmrN1NHUwqL2VSBCOqABbOOQbi1bzuSXVIXlJ1GWUYveNMoZJuB84naCZ4E/gG8C3gLknXAjuBT4SnrwMuAbYBbcBn0ou8OOxp6oislArw2muvccqkDk+qw/Ck6jJuvM+JzezKIQ4tT3GuAX89zo8qWrub2iNNqk8++STH9Tbz7P4pkd2z2GS1+i9phaStYQftr6U4Xi7pzvD4s2Fn68SxpZKeDsdfb5QU3TfFZYxxdO+HidgLIld2N3VE9pAK4KMf/SjVx7+bNw6209fnf4epZC2pSooBPyDopH0ywdPTkweddi1w0MwWA98Dvh1eWwr8HPi8mZ1CUCUszDWbJ5jEg6pUm8us7t4+9rV0MiuiPqoAU6ZMYf6serp6+ni7uSOy+xaTbJZUzwK2mdl2M+sC7iDosJ0suWP3WmB52E/wImCDmf0RwMz2m5n3Pi4QvWG76uDNZdbe5k7MoutOBbBt2zaqOxoBeH2/t6umks2kOprO2f3nmFkP0ATUAccDJulBSS9I+koW4nURCEZUeVLNhT1hH9Uo21Sfeuop3n7lRQB2+sOqlArlQVUpcC5wJsGT3UclPW9mjyafJGkVwbhvxNRsx+hSSLSpDnnQZcxbh4LqeZTV/8svv5zu3j6+8+2nfFTVELJZUh1N5+z+c8J21CnAfoJS7ZNm1mhmbQTdZ46ahSh58gqpOgM/ghu71KVUL6lmXtQd/wFqamqonTKZ2VMr2OHV/5SymVSfA5ZIWigpDlxB0GE7WXLH7suBx8KuMg8Cp0mqCpPt+4CXshS3S0PwoMpSbi6zdjd1UB2PMTmiaf8Atm7dytatWzljXi2PbXmbg61dkd27WGQtqYZtpNcRJMgtwF1mtlnSNyV9ODztp0CdpG3AlwmndTOzg8B3CRLzi8ALZnZ/tmJ36TGl3lxm7Tkc9FENnvVG4+mnn+bpp5/m8+cfR2tXLz/77WuR3btYZLVN1czWMWjMtJl9Pel1B/DxIa79OUG3KldAggdVxbVETKEI+qhGu/TJJz4RDGCrqqpixSkz+dnvdnDtexcxpTKaYbDFwMf+u4yyIar+Xv3PvN2Hoh2iCkEyraoK1ru67gOLae7o4Zbf7Yj0MwqdJ1WXcT6iKvt6evvY2xztaCqALVu2sGXLFgBOnTOF5Scew0+feo2mNh+Lk+BJ1WWU91PNjX0tnfRZtN2pAJ599lmeffbZ/vdfvuh4mju6+e7DWyP9nEJWKP1UXYEKnv57m2q2JWb8j7qkesUVVwx4f8rsKVx9znxue2YnH284llPn+EQrXlJ1Geb9VHPh8ZeDubsX1EfbX7uiooKKioGJ+ssXnUBtVZyv37PJJ1nBk6rLMBP0qC/lVgxGmnktF9440MZPntzOymWzWRhxUt20aRObNg1cjmtKZRnXX3ISL7x+iNX/Z3ukn1eIvPrvMirR+b8YJc28diHBqL/nJN1rZjkdmHLj/VuISXzt4hMjv/f69cGSYqeeeuqA/X92+hwee/ltvvXrl6mtKuOTZ86L/LMLhSdVl1HBzP/FUSpNoX/mNQBJiZnXUibVXjMOd0T7lNzCtWn6zNjX0sn6HQd5YPMe/u6i4yN/SAVw1VVXpdwvie99chktnc9z/S830tjSxZkLprGgvop4rARJRDgGIecmD7M8jSdVl3FF/KAq1cxrZw918ktvHWbpPzyU8aAW1Vfzl+9dlJF7l5UNnUzKS2P85C/O4LM3P8d3Hizu3gA7vnXpkMc8qbqMMqCvSKv/o5E8c1r9nIX835eelInPQEBdTZy5tZWcOHMyFWWxyD8HYMOGDQAsXbo05fHKeIx//6uzeaupgz+93cwbB9ro7TN6J9ADLE+qLqMM6Cnejv4jzrxmZquB1QANDQ2WqRJktrzwwgvA0EkVgiQ/Z2olc6ZG3/xQCDypuowyjO7irf73z7xGkEyvAP48tyFl1tVXX53rEPKed6lyGddDX8ptNPKxy1LCUDOv5TaqzIrFYsRimWlaKBZeUnUZ1YfRpfEtJ5avXZaSpZp5rZi9+OKLACxbtiynceQzT6ouowzoGn9H/zF1WXKZ50l1ZAom1i8+kvYBOzP8MfVAY4Y/oxBiAKgws1MH75T0AEGMKa8Bktc5Xh0+2Elcezmwwsz+Mnx/NXC2mV0XXdjZI6kZyPe+RvnyfRpJruNsNLMVqQ4UbUnVzKZn+jMkrTezhkx/Tr7HkIgj1f6hvngT1NZ8+LsaTr58n0aSz3H6gyqXz0azWKRzecWTqstno1ks0rm8UrTV/yxZPfIpGZcPMUAG4jCzHkmJLksx4KYC77KUL39XwymEGCGP4yzaB1XOOZcLXv13zrkIeVJ1zrkIeVJNg6TvSHpZ0gZJd0uamqM4Pi5ps6Q+SVntZpLPw0gzbaSfXVK5pDvD489KWhDuXyCpXdKL4fbjpGvOkLQxvOafpfRmIU0jxquS4nsx/G4tC489Ed4zceyYDMd4nqQXJPWEfZeTj10j6ZVwuyZpf6S/xzExM9/GuQEXAaXh628D385RHCcBJwBPAA1Z/NwY8CqwCIgDfwROzvXfS7787MAXgB+Hr68A7gxfLwA2DXHf3wPnAAJ+DVycixgHnXMa8GrS+8i+Z6OMcQGwFLgVuDxp/zRge/hnbfi6Nurf41g3L6mmwcwesmBSDYBnCPpR5iKOLWaWi5E6/cNIzawLSAwjnQhG87OvBG4JX68Flg9XYpI0C5hsZs9YkBluBT6SBzFeGV6bCSPGaGY7zGwDHDULzweBh83sgJkdBB4GVmTg9zgmnlSj81mC/xEnklQz38/JUSzZNpqfvf+c8D/fJqAuPLZQ0h8k/UbSe5POf3OEe2YzxoRPArcP2vezsOr/92lWrdP5Dg11bdS/xzHxfqojkPQIMDPFoRvM7J7wnBuAHmBNLuNwBWM3MM/M9ks6A/hPSafkOqhUJJ0NtJlZ8hKqV5nZLkmTgF8AVxOUBh2eVEdkZhcMd1zSp4EPAcvDqkZO4siRiTyMdDQ/e+KcNyWVAlOA/eH3pBPAzJ6X9CpwfHh+chNSur/PcceYdPwKBpVSzWxX+GezpH8nqMKPN6mm8x3aBZw/6NoniP73OCZe/U+DpBXAV4APm1lbruPJgYk8jHQ0P/u9QOKJ9OXAY2ZmkqaHc8UiaRGwBNhuZruBw5LOCavUnwLSqYWMO8YwthLgEyS1p0oqlVQfvi4jKFBsYvzS+Q49CFwkqVZSLcGD4wcz8Hscm2w9ESvGDdhG0KbzYrj9OEdxfJSg3agTeDv8YmXrsy8B/kTwBPeGXP+dZPn3ftTPDnyT4D9ZCKY2/I/we/J7YFG4/8+AzeF35gXgsqR7NhAkqVeBfyEc9ZjtGMNj5wPPDLpfNfA8sCH8Gb4PxDIc45nh97uVoBS9Oenaz4axbwM+k6nf41g2H6bqnHMR8uq/c85FyJOqc85FyJOqc85FyJOqc85FyJOqc85FyJOqc85FyJOqcy4tkmKSvh9OP7kxHNAwYXlSjYCkljSuvUnSXkmbBu1POcekpMpwEo7YCPeNS3oyHHroXCZdTzAi7BTgnwmmE5ywPKnm3s3AiuQdYcL8AXAxcDJwpaSTw8OfBX5pZr3D3dSCadQeJZhhyLmMkFQNfNTMvh/ueg1YnMOQcs6TaoQkfVnSpnD7UtL+vw9LnU9Jul3S3yWOmdmTwIFBtxpujsmrSBrHLOkdYYn0pXB2dpP0zfDwf4bnO5cpFwDHJlYBAG7i6O/zhOJVw4iEU7h9BjibYLbxZyX9huB3/GfAO4AygrHez49wu1TzRJ4dTjixyMx2hJ9ZAdwJfMrMfi/pfxCM5f5GeN0mgnHTzmXKMuDrZvZjAEn/BmwI21VvAKaY2eXDXF90vKQanXOBu82s1cxagF8C7wXeA9xjZh1m1gz8Ko3PqAcOJb2/AHjBzH4fvt8ATLNwQoewiaArnPfSuUyoBdogmMGKYKaoX4W1rGtzGlmOeFLNT0PNMdlOUBJNOBXYmPT+dIKScLJyoCMDMToHwexS54Sv/wa438xey2E8OedJNTr/B/iIpKpE432477fAZZIqJNUQzD85kpRzTFqwDk8srPZDMA3aUgBJxwMfY+Dcl3VAo5l1R/MjOneU24HTJW0j+C5+Ocfx5Jy3qUbEzF6QdDPBnJQA/2ZmfwCQdC9B1fxtgpJlU+I6SbcTzFtZL+lN4Btm9lNJ1xFMwhsDbjKzzeElDxE0NTxC8IX+cNgdqxG40sySZ21/P3B/Bn5c5wAI/6M/Z/D+8D/0G4F3SrrezP5n1oPLEZ9PNQsk1ZhZi6Qq4ElglZkNrqaP9l6nA39jZleP4txfAl8zsz+N57Occ2PnJdXsWB32M60AbhlvQoX+EvHjkmLD9VUNmw3+0xOqc9nlJVXnnIuQP6hyzrkIeVJ1zrkIeVJ1zrkIeVJ1zrkIeVJ1zrkIeVJ1zrkIeVJ1zrkI/f+ibQUSHnwBVgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pyabc.visualization.plot_kde_matrix_highlevel(\n", " h,\n", " limits=importer.get_bounds(),\n", " refval=importer.get_nominal_parameters(),\n", " 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": [ "# Change this line to run the code\n", "if False:\n", " import importlib\n", " import sys\n", "\n", " import amici\n", " import amici.petab_import\n", " import pandas as pd\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", " {\n", " 'observableId': 'obs_x2',\n", " 'simulationConditionId': 'condition1',\n", " 'measurement': obs_x2,\n", " 'time': meas_times,\n", " 'noiseParameters': 'sigma',\n", " }\n", " )\n", "\n", " # store data\n", " df.to_csv(\n", " dir_in + model_name[:-5] + '_measurement_table.tsv',\n", " sep='\\t',\n", " index=False,\n", " )" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.7" } }, "nbformat": 4, "nbformat_minor": 4 }