qmut-tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantitative mutations using `qmut`\n",
    "The package qmut contains mainly a class called QuantitativeMutation to enable the simulation of quantitative mutations of whole genome metabolic reconstructions. \n",
    "\n",
    "## Installation\n",
    " 1. Check that you have `cobrapy` installed in your environment\n",
    " 2. Download and copy the `\\cobrapy_qmut` folder in your working directory.\n",
    "\n",
    "\n",
    "## Turotial\n",
    "In this tutorial we will show the main applications and methods of the QuantiativeMutation objects, briefly:\n",
    "\n",
    "\n",
    "1. Load a metabolic model and create its QuantitativeMutation object\n",
    "2. Using `Q.optimize()` and `Q.slim_optimize()` to find solutions and biomass production rate, respectively.\n",
    "3. Generating random media following the protocol of [Wang and Zhang](https://doi.org/10.1371/journal.pgen.1000329) and applying it to `Q`.\n",
    "4. Compute maximal bounds of `Q` across a list of media. And save/load them with `Q.load_bounds()` and `Q.save_bounds()`\n",
    "5. Check how gene reaction rules are quantitatively interpreted.\n",
    "6. Working with relative gene dosages and quantitative mutations.\n",
    "7. Compute individuals, and populations in a single or in multiple processors.\n",
    "\n",
    "\n",
    "### 0. First load some useful packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cobra\n",
    "import cobrapy_qmut\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Create qmut object and define standard medium and solver\n",
    "We create an instance of QuantitativeMutation using a path to a model file in json format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Academic license - for non-commercial use only - expires 2022-08-06\n",
      "Using license file /home/pyubero/gurobi.lic\n",
      "Loaded model iND750 from user file.\n",
      "Biomass reaction id is BIOMASS_SC4_bal\n",
      "The model has 750 genes and 1266 reactions.\n",
      "... number of exchange reactions:\t116\n",
      "... number of non-exchange reactions:\t1148\n"
     ]
    }
   ],
   "source": [
    "modelpath = './cobrapy_qmut/models/yeast_iND750.json'\n",
    "\n",
    "Q = cobrapy_qmut.QuantitativeMutation('./cobrapy_qmut/models/yeast_iND750.json', verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that `Q.model` can be accessed as a regular `cobra.Model`. For example, we can define the standard medium, and show the current solver used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solver being used:  <class 'optlang.glpk_interface.Model'>\n",
      "\n",
      "Standard medium composition:\n",
      "     EX_na1_e\t9999.0\n",
      "     EX_nh4_e\t9999.0\n",
      "      EX_o2_e\t2.0\n",
      "      EX_pi_e\t9999.0\n",
      "     EX_so4_e\t9999.0\n",
      "     EX_co2_e\t9999.0\n",
      "  EX_glc__D_e\t20.0\n",
      "     EX_h2o_e\t9999.0\n",
      "       EX_k_e\t9999.0\n"
     ]
    }
   ],
   "source": [
    "std_medium = {}\n",
    "std_medium.update( {  'EX_o2_e'  : 2 } )\n",
    "std_medium.update( {  'EX_co2_e' : 9999 } )\n",
    "std_medium.update( {  'EX_k_e'   : 9999 } )\n",
    "std_medium.update( {  'EX_na1_e' : 9999 } )\n",
    "std_medium.update( {  'EX_nh4_e' : 9999 } )\n",
    "std_medium.update( {  'EX_pi_e'  : 9999 } )\n",
    "std_medium.update( {  'EX_so4_e' : 9999 } )\n",
    "std_medium.update( {  'EX_h2o_e' : 9999 } )\n",
    "std_medium.update( {  'EX_glc__D_e' : 20 } )\n",
    "\n",
    "Q.model.medium = std_medium\n",
    "Q.model.solver = 'glpk'\n",
    "\n",
    "print('Solver being used: ', type(Q.model.solver))\n",
    "print('')\n",
    "print('Standard medium composition:')\n",
    "_ = [ print('%13s\\t%1.1f' % (key, value)) for key, value in Q.model.medium.items() ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Optimizing the model\n",
    "We can also compute the biomass production, or any other optimization problem we may have set in different ways.\n",
    "\n",
    "- Both original methods from the `cobra.Model` object work exactly as expected.\n",
    "    1. Through the `cobra.Model.optimize()` method\n",
    "    2. Through `cobra.Model.slim_optimize()` \n",
    "\n",
    "\n",
    "- Directly as `Q` methods\n",
    "    3. Through Q's own `Q.optimize()`, work exactly as `Q.model.optimize()` for now...\n",
    "    4. Through `Q.slim_optimize()`, which always returns the biomass production rate.\n",
    "\n",
    "<> Although both return the same object, `Q`'s own method applies the current gene dosages within a model context manager. By now, there is no difference between both methods but only because our Q instance is unbounded, just as Q.model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution 1 objective value:  0.513926\n",
      "Solution 2 objective value:  0.513926\n",
      "Q.model    objective value:  0.513926\n",
      "Q's biomass production rate: 0.513926\n",
      "\n",
      "Are the fluxes of both solutions equal?  True\n"
     ]
    }
   ],
   "source": [
    "solution1 = Q.optimize()\n",
    "solution2 = Q.model.optimize()\n",
    "mu1 = Q.slim_optimize()\n",
    "mu2 = Q.model.slim_optimize()\n",
    "\n",
    "print('Solution 1 objective value:  %1.6f' % solution1.objective_value )\n",
    "print('Solution 2 objective value:  %1.6f' % solution2.objective_value )\n",
    "print('Q.model    objective value:  %1.6f' % mu1 )\n",
    "print(\"Q's biomass production rate: %1.6f\" % mu2 )\n",
    "print('')\n",
    "print('Are the fluxes of both solutions equal? ', all(solution1.fluxes==solution2.fluxes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, `Q.slim_optimize()` will not always be faster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.55 ms ± 11.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "1.58 ms ± 4.64 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "13.5 ms ± 2.5 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "14.2 ms ± 133 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit Q.model.slim_optimize()\n",
    "%timeit Q.slim_optimize()\n",
    "\n",
    "%timeit Q.model.optimize()\n",
    "%timeit Q.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Generate random media \n",
    "The function `Q.random_medium()` generates a random medium dict following the protocol originally from [Wang and Zhang] but widely used by [others].\n",
    "\n",
    "First define a minimal composition of all media.\n",
    "\n",
    "However, generating a random media **DOES NOT** mean using it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Random medium composition:\n",
      "     EX_pap_e\t1.392\n",
      "  EX_rib__D_e\t13.664\n",
      "  EX_sbt__D_e\t1.116\n",
      "   EX_thmmp_e\t6.491\n",
      "  EX_thr__L_e\t5.063\n",
      "   EX_ttdca_e\t13.775\n",
      "     EX_cit_e\t1.088\n",
      "    EX_dann_e\t11.215\n",
      "    EX_dcyt_e\t14.384\n",
      "     EX_din_e\t16.331\n",
      "     EX_fum_e\t14.683\n",
      "     EX_gal_e\t10.480\n",
      "  EX_ile__L_e\t4.397\n",
      "       EX_k_e\t9999.000\n",
      "   EX_8aonn_e\t11.629\n",
      "      EX_o2_e\t2.000\n",
      "     EX_h2o_e\t9999.000\n",
      "     EX_co2_e\t9999.000\n",
      "     EX_nh4_e\t9999.000\n",
      "      EX_pi_e\t9999.000\n",
      "     EX_so4_e\t9999.000\n",
      "     EX_na1_e\t9999.000\n",
      "\n",
      "Is this random composition the same used in Q.model?  False\n"
     ]
    }
   ],
   "source": [
    "minimal_medium = {}        \n",
    "minimal_medium.update( {'EX_o2_e' :     2 }) #... oxygen\n",
    "minimal_medium.update( {'EX_h2o_e' : 9999 }) #... water\n",
    "minimal_medium.update( {'EX_co2_e' : 9999 }) #... carbon dioxyde\n",
    "minimal_medium.update( {'EX_nh4_e' : 9999 }) #... ammonia\n",
    "minimal_medium.update( {'EX_pi_e'  : 9999 }) #... phosphate\n",
    "minimal_medium.update( {'EX_so4_e' : 9999 }) #... sulphate\n",
    "minimal_medium.update( {'EX_k_e'   : 9999 }) #... potassium\n",
    "minimal_medium.update( {'EX_na1_e' : 9999 }) #... sodium\n",
    "\n",
    "rand_med = Q.random_medium(m = 0.1, minimal=minimal_medium, lower_lim=0, upper_lim = 20, method='linear')\n",
    "print('Random medium composition:')\n",
    "_ = [ print('%13s\\t%1.3f' % (key, value)) for key, value in rand_med.items() ]\n",
    "\n",
    "print('\\nIs this random composition the same used in Q.model? ', rand_med == Q.model.medium )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Using context managers to temporally modify `Q.model`\n",
    "You can use context managers as with `cobrapy` to temporally modify the model,e.g. the medium, and quickly test its viability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Biomass production in the random medium:   0.5357\n",
      "Biomass production in the standard medium: 0.5139\n"
     ]
    }
   ],
   "source": [
    "with Q.model:\n",
    "    Q.set_medium( rand_med )\n",
    "    print('Biomass production in the random medium:   %1.4f' % Q.slim_optimize() )\n",
    "print('Biomass production in the standard medium: %1.4f' % Q.slim_optimize() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will create a list of 100 random media and measure the growth rate in all of them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 100/100 [00:05<00:00, 16.97it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished computing.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "from tqdm import tqdm \n",
    "\n",
    "NMEDIA = 100\n",
    "\n",
    "mu           = np.zeros((NMEDIA,))\n",
    "random_media =[[],]*NMEDIA\n",
    "\n",
    "for jj in tqdm(range(NMEDIA)):\n",
    "    with Q.model:\n",
    "        new_medium       = Q.set_medium( Q.random_medium( m = 0.10, minimal = minimal_medium) )\n",
    "        mu[jj]           = Q.slim_optimize()\n",
    "        random_media[jj] = new_medium\n",
    "    \n",
    "print('Finished computing.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfFElEQVR4nO3de5gcZZ328e9NAgYNghIcIRwCGMEgyCGCosgED5uACqsoJ2FhxRg1u+riCrrKIiqLy/K+ctJs9MWAssRVERCieByQoxAEAiRAwAgBJBwDE4IQ+L1/1DNYaXqmO5mpmae77891zZWuqqerfk/XpO+pQz+tiMDMzCw364x0AWZmZvU4oMzMLEsOKDMzy5IDyszMsuSAMjOzLDmgzMwsSw4oa1mSQtLrRrqONSVplqQvj3Qd9Ug6UtKVI13HYEnqkXR0enyYpF+OdE225hxQtkYkHSzpOkkrJC1Ljz8pSRVv98U3nOFSVQBGxIyI+OpQr9fqi4jzIuI9I12HrTkHlDVN0jHAacApwGuBLmAG8DZgvX6eM2rYClwDkkaPdA2DpYL/D1vb8i+3NUXShsCJwCcj4scR8VQU/hgRh0XEX1O7OZK+LWmepBXAFElvSEdAT0i6TdL7U9ut07x10vR3JS0rbfMHkj4j6evAXsCZknolnVkq7V2S7pL0uKSz+juSk3SCpB+ndT4JHClpd0nXpBoelHSmpPVS+yvSU29O2zwozX+vpJvSc66WtFNpG8dKul/SU5LukPTOfmqZI+lr6XG3pKWSjklHpA9KOmqA/dAj6euSrgKeBraRdJSkhWm790j6eKn9gOuXtLGkiyU9KekPwLY129tT0vWSlqd/96yp5WvpdeiV9LO0vvPS+q6XNKGffkxIR6hHSbov7b8Zkt4s6Zb0+p5Z85x/TP18XNJlkrYqLXu3pEWpzjMBlZatdtpS0mlpm09Kmi9pr/5ebxthEeEf/zT8AaYCq4DRDdrNAZZTHFWtA2wALAa+SHGUtQ/wFLBdan8vsFt6fAdwD/CG0rJd0uMe4OiabQVwCbARsCXwMDC1n7pOAJ4DDkh1rQ/sBrwFGA1MABYCn6lZ/+tK07sCy4A9gFHAPwBLgJcB2wH3AZulthOAbQd4jb6WHnen1/VEYF1gX4rgeVU/z+1Jr8sOqe51gf0ogkXA3un5uzazfmAu8L/AK4A3AvcDV6ZlrwYeBw5P2zokTW9cqmVx2vaGwO3AncC7Uvtzge/1048J6fWdBYwB3gM8A1wIvAYYn17rvVP7A9K23pDW/SXg6rRsHPAkcGDq42dTn49Oy4/s61Oa/giwcVrPMcBfgDEj/X/MPy/98RGUNWsc8EhErOqbkf5yfkLSSknvKLW9KCKuiogXgJ2BscDJEfFsRPyWIlQOSW0vB/aW9No0/eM0vTXwSuDmBnWdHBFPRMS9wO/S9vpzTURcGBEvRMTKiJgfEddGxKqIWAL8N8UbfH8+Bvx3RFwXEc9HxDnAXylC7nmKoJokad2IWBIRdzeovc9zwIkR8VxEzAN6KQKvP3Mi4rZU93MRcWlE3B2Fy4FfUhxxDrh+FadfPwgcHxErIuJW4JzS8/YD7oqI76dtnQ8sAt5XavO9tO3lwM+BuyPi1+n35EfALg36/tWIeCYifgmsAM6PiGURcT/w+9LzPw78R0QsTOs+Cdg5HUXtC9wexZH9c8A3KUKnroj4QUQ8mvp0Kn/7A8My44CyZj0KjFPp2k1E7BkRG6Vl5d+l+0qPNwPuS2HV588UfyFDEVDdwDuAKyj+Kt87/fy+5nn1lN+InqYIw/6U60LS6yVdIukv6bTfSRRB3J+tgGNSKD8h6QlgC4qjpsXAZyiO1JZJmitpswa193m0HPxr0Y9pkq6V9Fiqad+afvS3/k0ojiLK6/tz6fFmNdN9y8eXph8qPV5ZZ3qgfqzJ87cCTiu97o9RHDGOT3W+2IeICGpeo7J0unNhOh34BMXR30D73UaIA8qadQ3F0cL+TbQtD5H/ALCFVr+YvyXFqSQoAmovipC6HLiS4vTg3mm63jrXVu06vk1xRDAxIl5JcRpyoLsR7wO+HhEblX5eno4siIj/iYi3U7yZBvCNIai5nhf7IellwE+A/wK60h8M8xr0o8/DFKfCtijN27L0+AGKvlCz/H6G333Ax2te+/Uj4mrgQUp9kCRW7xOlZXsBxwIfpjjNuRHFKelK70K1teOAsqZExBPAV4BvSTpQ0lhJ60jameL6RX+uozh183lJ60rqpjhFNDet9y6Kv5Q/AlwREU9S/BX9QVYPqIeAbYayTxTXx54EeiVtD3yiZnntNr8DzJC0hwqvkLSfpA0kbSdpnxQYz6Q+PT/E9dazHsUpqoeBVZKmUVzPaSgingcuAE6Q9HJJkyiuq/WZB7xe0qGSRqu4UWQSxSna4TYL+IKkHaC4aUfSh9KyS4EdJH0gHeH/M8VdpvVsQBHKDwOjJR1PcSrZMuSAsqZFxH8C/wJ8nuIC9kMU122OBa7u5znPAu8HpgGPAN8CjoiIRaVml1Ochrq3NC3gj6U2pwEHpju4Th+iLn0OOJTipo3vAD+sWX4CcE46rfThiLiB4jrUmRQ3CyymuAAPRUicTNHHv1Bc6P/iENXZr4h4iuIN+X9TTYcCF6/BKmZSnEb7C8XNG98rrftR4L0UNxI8SrHf3xsRjwxF7WsiIn5KcUQ6N52OvZXid4pUz4coXv9HgYnAVf2s6jKKa2V3UpyufIYBTgfayFJxutbMzCwvPoIyM7MsOaDMzCxLDigzM8tSZQEl6ew0tMqt/SyXpNMlLU5Dm+xaVS1mZtZ6qhwwcw7F3U7n9rN8GsXdNhMpho75dvp3QOPGjYsJEyYMqrAVK1bwilcMdGd0e+m0/kLn9dn9bX/t3Of58+c/EhGb1M6vLKAi4or+BopM9gfOTZ/6vlbSRpI2jYgHB1rvhAkTuOGGGwZVW09PD93d3YNaRyvptP5C5/XZ/W1/7dxnSbUjlgDVHkE1Mp7VP3+wNM17SUBJmg5MB+jq6qKnp2dQG+7t7R30OlpJp/UXOq/P7m/768Q+j2RA1RtapO6HsiJiNjAbYPLkyTHYvyLa+S+Rejqtv9B5fXZ/218n9nkk7+JbyurjZW1OMfaXmZnZiAbUxcAR6W6+twDLG11/MjOzzlHZKT5J51OMUD1O0lLg3ym+TIyImEUxEOW+FOOZPQ30+y2iZmbWeaq8i++QBssD+FRV2zczs9bmkSTMzCxLDigzM8uSA8rMzLI0kp+DGjEL7l/OkcddOtJlsOTk/Ua6BDOzbPkIyszMsuSAMjOzLDmgzMwsSw4oMzPLkgPKzMyy5IAyM7MsOaDMzCxLHfk5qFxMGKbPYh2z46oBP/flz2OZWY58BGVmZllyQJmZWZYcUGZmliUHlJmZZckBZWZmWXJAmZlZlhxQZmaWJQeUmZllyQFlZmZZckCZmVmWHFBmZpYlB5SZmWXJAWVmZllyQJmZWZYcUGZmliUHlJmZZckBZWZmWXJAmZlZlhxQZmaWJQeUmZllyQFlZmZZckCZmVmWKg0oSVMl3SFpsaTj6izfUNLPJN0s6TZJR1VZj5mZtY7KAkrSKOAsYBowCThE0qSaZp8Cbo+INwHdwKmS1quqJjMzax1VHkHtDiyOiHsi4llgLrB/TZsANpAkYCzwGLCqwprMzKxFKCKqWbF0IDA1Io5O04cDe0TEzFKbDYCLge2BDYCDIuLSOuuaDkwH6Orq2m3u3LmDqm3ZY8t5aOWgVtFSutZnwP7uOH7D4StmmPT29jJ27NiRLmPYuL/tr537PGXKlPkRMbl2/ugKt6k682rT8O+Am4B9gG2BX0n6fUQ8udqTImYDswEmT54c3d3dgyrsjPMu4tQFVXY9L8fsuGrA/i45rHv4ihkmPT09DPb3pJW4v+2vE/tc5Sm+pcAWpenNgQdq2hwFXBCFxcCfKI6mzMysw1UZUNcDEyVtnW58OJjidF7ZvcA7ASR1AdsB91RYk5mZtYjKznNFxCpJM4HLgFHA2RFxm6QZafks4KvAHEkLKE4JHhsRj1RVk5mZtY5KL8RExDxgXs28WaXHDwDvqbIGMzNrTR5JwszMsuSAMjOzLDmgzMwsSw4oMzPLkgPKzMyy5IAyM7MsOaDMzCxLDigzM8uSA8rMzLLkgDIzsyw5oMzMLEsOKDMzy5IDyszMsuSAMjOzLDmgzMwsSw4oMzPLkgPKzMyy5IAyM7MsOaDMzCxLDigzM8uSA8rMzLLkgDIzsyw5oMzMLEsOKDMzy5IDyszMsuSAMjOzLDmgzMwsSw4oMzPLkgPKzMyyNHqkC7CRN+G4S0e6BACWnLzfSJdgZhnxEZSZmWXJAWVmZllyQJmZWZYcUGZmliUHlJmZZanSgJI0VdIdkhZLOq6fNt2SbpJ0m6TLq6zHzMxaR2W3mUsaBZwFvBtYClwv6eKIuL3UZiPgW8DUiLhX0muqqsfMzFpLlUdQuwOLI+KeiHgWmAvsX9PmUOCCiLgXICKWVViPmZm1EEVENSuWDqQ4Mjo6TR8O7BERM0ttvgmsC+wAbACcFhHn1lnXdGA6QFdX125z584dVG3LHlvOQysHtYqW0rU+LdHfHcdvOGTr6u3tZezYsUO2vty5v+2vnfs8ZcqU+RExuXZ+lSNJqM682jQcDewGvBNYH7hG0rURcedqT4qYDcwGmDx5cnR3dw+qsDPOu4hTF3TOIBrH7LiqJfq75LDuIVtXT08Pg/09aSXub/vrxD5X+a61FNiiNL058ECdNo9ExApghaQrgDcBd2JmZh2tymtQ1wMTJW0taT3gYODimjYXAXtJGi3p5cAewMIKazIzsxZR2RFURKySNBO4DBgFnB0Rt0makZbPioiFkn4B3AK8AHw3Im6tqiYzM2sdlV6YiIh5wLyaebNqpk8BTqmyDjMzaz0eScLMzLLkgDIzsyw5oMzMLEsOKDMzy5IDyszMsjRgQEnKf/gBMzNrS42OoP7Q90DSGRXXYmZm9qJGAVUeT+9tVRZiZmZW1iigqhnq3MzMrIFG15i2l3QLxZHUtukxaToiYqdKqzMzs47VKKDeMCxVmJmZ1RgwoCLiz/DiV7NPTLPvjIjlFddlZmYdbsCASl+TMRs4APgTxam9rST9FJiRvsrdzMxsyDW6SeJLFF/JvkVE7BIROwNbUgTblyuuzczMOlijgPoA8LGIeKpvRnr8SeDvqyzMzMw6W6OAeiEinq6dGRG9+BZ0MzOrUKO7+ELSq1j9A7t9XqigHjMzM6BxQG0IzKd+QPkIyszMKtPoNvMJw1SHmZnZahqNZv4aSd+UdImkkyS9crgKMzOzztboJolzgRXAGcAGwOmVV2RmZkbja1CvjYh/S48vk3Rj1QWZmZlB44BSzV18o8rTEfFYlcWZmVnnWpO7+CL923cUFcA21ZVmZmadrKm7+CStAxwGbB0RJ0raEti0+vKsk0w47tIhW9cxO67iyLVc35KT9xuyOsxs7TW6SaLPWcBbgEPS9FPAmZVUZGZmRuNTfH32iIhdJf0RICIeTyOdm5mZVaLZI6jnJI0ijR4haRM81JGZmVWo2YA6Hfgp8BpJXweuBE6qrCozM+t4TZ3ii4jzJM0H3klxJ98BEbGw0srMzKyjNXsNiohYBCyqsBYzM7MXNXuKz8zMbFg5oMzMLEsOKDMzy1KlASVpqqQ7JC2WdNwA7d4s6XlJB1ZZj5mZtY7KAip9buosYBowCThE0qR+2n0DuKyqWszMrPVUeQS1O7A4Iu6JiGeBucD+ddr9E/ATYFmFtZiZWYupMqDGA/eVppemeS+SNB74e2BWhXWYmVkLavpzUGtBdeZFzfQ3gWMj4nmpXvO0Imk6MB2gq6uLnp6eQRXWtX4x2nWn6LT+wuD6PNjfr5HQ29vbknWvrU7rL3Rmn6sMqKXAFqXpzYEHatpMBuamcBoH7CtpVURcWG4UEbOB2QCTJ0+O7u7uQRV2xnkXceqCKruel2N2XNVR/YXB9XnJYd1DW8ww6OnpYbD/L1pJp/UXOrPPVb5rXQ9MlLQ1cD9wMHBouUFEbN33WNIc4JLacDIzs85UWUBFxCpJMynuzhsFnB0Rt0makZb7upOZmfWr0vM+ETEPmFczr24wRcSRVdZiZmatxSNJmJlZlhxQZmaWJQeUmZllyQFlZmZZ6qwPx5g1YcJxl450CQAsOXm/kS7BbET5CMrMzLLkgDIzsyw5oMzMLEsOKDMzy5IDyszMsuSAMjOzLDmgzMwsSw4oMzPLkgPKzMyy5IAyM7MsOaDMzCxLDigzM8uSA8rMzLLkgDIzsyw5oMzMLEsOKDMzy5IDyszMsuSAMjOzLDmgzMwsSw4oMzPLkgPKzMyy5IAyM7MsOaDMzCxLDigzM8uSA8rMzLLkgDIzsyw5oMzMLEsOKDMzy5IDyszMsuSAMjOzLFUaUJKmSrpD0mJJx9VZfpikW9LP1ZLeVGU9ZmbWOioLKEmjgLOAacAk4BBJk2qa/QnYOyJ2Ar4KzK6qHjMzay1VHkHtDiyOiHsi4llgLrB/uUFEXB0Rj6fJa4HNK6zHzMxaiCKimhVLBwJTI+LoNH04sEdEzOyn/eeA7fva1yybDkwH6Orq2m3u3LmDqm3ZY8t5aOWgVtFSutano/oL7dHnHcdv2HTb3t5exo4dW2E1eem0/kJ793nKlCnzI2Jy7fzRFW5TdebVTUNJU4CPAm+vtzwiZpNO/02ePDm6u7sHVdgZ513EqQuq7HpejtlxVUf1F9qjz0sO6266bU9PD4P9f9FKOq2/0Jl9rvJ/8FJgi9L05sADtY0k7QR8F5gWEY9WWI+ZmbWQKq9BXQ9MlLS1pPWAg4GLyw0kbQlcABweEXdWWIuZmbWYyo6gImKVpJnAZcAo4OyIuE3SjLR8FnA8sDHwLUkAq+qdhzQzs85T6Un6iJgHzKuZN6v0+GjgJTdFmJmZeSQJMzPLkgPKzMyy5IAyM7MsOaDMzCxLDigzM8uSA8rMzLLkgDIzsyw5oMzMLEsOKDMzy5IDyszMsuSAMjOzLDmgzMwsSw4oMzPLkgPKzMyy5IAyM7MsVfp9UGa29iYcd2nTbY/ZcRVHrkH7NbHk5P0qWa9ZIz6CMjOzLDmgzMwsSw4oMzPLkgPKzMyy5IAyM7MsOaDMzCxLvs3czKwFrcnHEKpU5ccQfARlZmZZckCZmVmWHFBmZpYlB5SZmWXJAWVmZllyQJmZWZYcUGZmliV/DsrMBtQJn7exPPkIyszMsuSAMjOzLDmgzMwsSw4oMzPLUqUBJWmqpDskLZZ0XJ3lknR6Wn6LpF2rrMfMzFpHZQElaRRwFjANmAQcImlSTbNpwMT0Mx34dlX1mJlZa6nyCGp3YHFE3BMRzwJzgf1r2uwPnBuFa4GNJG1aYU1mZtYiqvwc1HjgvtL0UmCPJtqMBx4sN5I0neIIC6BX0h2DrG0c8Mgg19Ey/rnD+gud1+dO6K++sdpk2/e3jiz7XLNf1tZW9WZWGVCqMy/Wog0RMRuYPRRFAUi6ISImD9X6ctdp/YXO67P72/46sc9VnuJbCmxRmt4ceGAt2piZWQeqMqCuByZK2lrSesDBwMU1bS4Gjkh3870FWB4RD9auyMzMOk9lp/giYpWkmcBlwCjg7Ii4TdKMtHwWMA/YF1gMPA0cVVU9NYbsdGGL6LT+Quf12f1tfx3XZ0W85JKPmZnZiPNIEmZmliUHlJmZZamtA6rThlpqor/dkpZLuin9HD8SdQ4VSWdLWibp1n6Wt9v+bdTfdtu/W0j6naSFkm6T9Ok6bdpmHzfZ37baxw1FRFv+UNyYcTewDbAecDMwqabNvsDPKT6P9RbgupGuu+L+dgOXjHStQ9jndwC7Arf2s7xt9m+T/W23/bspsGt6vAFwZ5v/H26mv221jxv9tPMRVKcNtdRMf9tKRFwBPDZAk3bav830t61ExIMRcWN6/BSwkGKkmbK22cdN9rejtHNA9TeM0pq2aRXN9uWtkm6W9HNJOwxPaSOmnfZvs9py/0qaAOwCXFezqC338QD9hTbdx/VUOdTRSBuyoZZaRDN9uRHYKiJ6Je0LXEgxkny7aqf924y23L+SxgI/AT4TEU/WLq7zlJbexw3625b7uD/tfATVaUMtNexLRDwZEb3p8TxgXUnjhq/EYddO+7ehdty/ktaleLM+LyIuqNOkrfZxo/624z4eSDsHVKcNtdSwv5JeK0np8e4U+//RYa90+LTT/m2o3fZv6sv/AxZGxP/pp1nb7ONm+ttu+7iRtj3FF3kPtTTkmuzvgcAnJK0CVgIHR7o1qBVJOp/irqZxkpYC/w6sC+23f6Gp/rbV/gXeBhwOLJB0U5r3RWBLaMt93Ex/220fD8hDHZmZWZba+RSfmZm1MAeUmZllyQFlZmZZckCZmVmWHFBmZpYlB5QNmqTn08jKN0u6UdKeaf5mkn480vVVLY0wfclaPvcASZNK0ydKetcQ1bWLpO8Oxbqa3N6E2pHWJc1Pn8vr7zknSLo//f7cJemCvtdD0v6SLiy1/YKkxaXp90m6uGb5YZLmSDqwzrY2kfSLQXbThpEDyobCyojYOSLeBHwB+A+AiHggIl7yRtEqJA3H5wQPAF4MqIg4PiJ+PUTr/iJwxhCta42l8eTuT4MXD+T/pt+ficAPgd9K2gS4Gnhrqd1bgSclvSZN7wlcVVr+HuCX/W0kIh4GHpT0tjXriY0UB5QNtVcCj8Pqf1FLGiPpe5IWSPqjpClp/pGSLpT0M0l/kjRT0r+kNtdKenVq9zFJ16ejtJ9Ienma/yFJt6b5V6R5O0j6Q/qr/BZJLxmrTFKvpFPTEd9v0hsiknoknSTpcuDTkt6Zalmg4vuYXpbaTZW0SNKVwAdK6z1B0udK07emN2okHZHquVnS99OR5vuBU1Kt25b/+h9g20skfSXVvkDS9nX6twGwU0TcXKrrHEm/TM//gKT/TM//hYohdvrWfZKkayTdIGlXSZdJulvpQ98qnJL6tkDSQf38LkwDfpGe81FJd6bX9zuSzqz3hIj4IUXIHJoCZbmk16XF4ymGAdozTe9JEWJIeiWwXnoOwDskXS3pnpqjqQuBw/qp1zLjgLKhsH56g10EfBf4ap02nwKIiB2BQ4BzJI1Jy94IHErxlSFfB56OiF2Aa4AjUpsLIuLN6ShtIfDRNP944O/S/PeneTOA0yJiZ2AyxXhttV4B3BgRuwKXU4zK0GejiNgbOAuYAxyU6h5N8Sn+McB3gPcBewGvbfQCqRh1+t+AfVKtn46IqymG6vnXdARxd6n9mHrbLq3ykVT7t4HP8VKTgdovNtwW2I/iKyp+APwurXtlmt/nvoh4K/D7VMOBFN+1dGJa/gFgZ+BNwLsoArbeV1xMBX4haTPgy2kd7wZeEqg1biy1uRrYU9J2wF3AtWl6NLATxRBfpDp+U1rHpsDbgfcCJ5fm30Cxz6wFOKBsKPSd4tue4k3pXEm1o0y/Hfg+QEQsAv4MvD4t+11EPNX3FzPwszR/ATAhPX6jpN9LWkDxF3Df1wxcBcyR9DGKIZ6gCLYvSjqWYuTnlXVqfoHidBIUb9ZvLy3rm78d8KeIuDNNn0PxpYHbp/l3pWFmfjDAa9NnH+DHEfFIeg0afa9Tf9vu0zeQ6Hz+9hqVbQo8XDPv5xHxHMXrOop0dMPqrzP8bQzHBRRfANi3b56RtBHFa3V+RDwfEQ9RBPybyxtScd1p84i4h+IPj8sj4rG0/R816Hv5d+cqiiOlPSn26x+APSi+iuKOiHgmtZtK8cWFfS6MiBci4nagqzR/GbBZg+1bJhxQNqQi4hpgHLBJzaJ6X4vQ56+lxy+Upl/gb+NFzgFmpr/4vwKMSdubAXyJYkTrmyRtHBH/Q3E0tRK4TNI+zZReeryiiZr7GyNsFav/v+o7StQAz6lnoG3D316j56k/pubK0rZXe05EvAA8VxrDrfw6l9dd3hfldo1qg+Io5cr0uJn2ZbtQHCVDOoJKP9ekL/IbQzEmYfn60+4U4dWnXHd5+2MoXhtrAQ4oG1LpesgoXjrC8hWkc/+SXk8xAOYda7DqDSgucK9L6RqCpG0j4rqIOB54BNhC0jbAPRFxOsXRwE511rcOxakrKE4vXlmnzSJgQukayOEURwuLgK0lbZvmH1J6zhKKr2VH0q7A1mn+b4APS9o4LXt1mv9U6luz227WQuB1DVutnSuAgySNStfu3sHq4QCrH9H8Adhb0qvSqbkP9rdiSR+kuNnh/DTrdoojnr2AP6Z5N1Gcxu27/rQDsCginm+i9tfz0lOflikHlA2FvmtQN1GcHvuHOm8W3wJGpVN0PwSOjIi/0rwvU3y76K8o3rz7nJIu1N9K8cZ5M3AQcGuqZ3vg3DrrWwHsIGk+xem3E2sbpNNHRwE/SnW/AMxK86cDl6abJP5cetpPgFenbX8CuDOt6zaK62uXS7oZ6Ps6hbnAv6abIfoCr99tN/NCpecvAjZMN0sMtZ8Ct1C81r8FPh8Rf6lp000K1Ii4HziJYv/9miJ0lpfafjb9/twFfITiOt3D6bmRnvdIOj0Ixam+bUgBRelmjCZMAS5tsq2NMI9mbh1JUm9EjB3pOqok6bPAUxExbJ+FStvdHPhOREwrzRubvgV2NEXAnR0RPx2i7f0KOKKZ74FScafn/hHx+FBs26rlgLKO1CEBNQb4UER8P4Na/oviTrsxFLeRf3q4v8conY58W0RcOJzbtbXngDIzsyz5GpSZmWXJAWVmZllyQJmZWZYcUGZmliUHlJmZZen/A7YOWRhPCazoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure( figsize=(6,4) )\n",
    "plt.hist(mu, density=True);\n",
    "plt.xlabel('Biomass production (mmol/gDW/h)')\n",
    "plt.ylabel('ePDF')\n",
    "plt.title('Growth rates in random media')\n",
    "plt.grid();\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Compute maximal bounds (lower and upper) in the standard medium and across random media\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000/1000 [00:23<00:00, 43.47it/s]\n",
      "100%|██████████| 1000/1000 [00:24<00:00, 40.94it/s]\n"
     ]
    }
   ],
   "source": [
    "NSAMPLES       = 1_000 # 1000 samples in 50 processores takes about 25s each\n",
    "MEDIA_RICHNESS = 0.10  # media richness scales with the number of components with which the minimal media is suplemented\n",
    "\n",
    "# Initialize empty outputs\n",
    "bounds_std = np.nan*np.zeros( (2,Q.N_NEX_RXNS) )\n",
    "bounds_rnd = np.nan*np.zeros( (2,Q.N_NEX_RXNS) )\n",
    "\n",
    "# Generate NSAMPLES random media\n",
    "random_media_list   = [ Q.random_medium( m = MEDIA_RICHNESS, minimal = minimal_medium) for _ in range(NSAMPLES) ]\n",
    "standard_media_list = [ std_medium for _ in range(NSAMPLES) ]\n",
    "\n",
    "# Compute bounds in standard medium\n",
    "Q.compute_bounds(media_list  = standard_media_list, \n",
    "                 pfba        = True, \n",
    "                 nprocessors = 50, \n",
    "                 chunksize   = 20 )\n",
    "\n",
    "\n",
    "# Extract bounds found in standard medium\n",
    "bounds_std[0,:] = [Q.upper_bounds[rxn_id] for rxn_id in Q.NEX_RXNS]\n",
    "bounds_std[1,:] = [Q.lower_bounds[rxn_id] for rxn_id in Q.NEX_RXNS]\n",
    "\n",
    "\n",
    "\n",
    "# Compute bounds in random medium\n",
    "Q.compute_bounds(media_list  = random_media_list, \n",
    "                 pfba        = True, \n",
    "                 nprocessors = 50, \n",
    "                 chunksize   = 20 )\n",
    "\n",
    "# Extract bounds found in random media\n",
    "bounds_rnd[0,:] = [Q.upper_bounds[rxn_id] for rxn_id in Q.NEX_RXNS]\n",
    "bounds_rnd[1,:] = [Q.lower_bounds[rxn_id] for rxn_id in Q.NEX_RXNS]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we compare the bounds obtained in standard vs random media, we find that the ones obtained in random med. tend to be larger than those obtained in standard med. The exceptions we find are \n",
    "- not great deviations from the diagonal, and \n",
    "- probably a problem due to using a small number of samples. \n",
    "\n",
    "For serious computations please use at least $10^4$ samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAGGCAYAAABmPbWyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAA9hAAAPYQGoP6dpAABItElEQVR4nO3deXxU1f3/8dcne8IOyqKyKQoqCAhUyShgK2Ct39al/bal7U+0da/7UqvtV9S2VqtWbdXaVktty7f9alu1VStoxQVwQ0AoghvIapA1QBKynd8f906YGSbkZjLJbO/n4zGPmbn3zJ3PZyaQT84591xzziEiIiIi+5eX6gBEREREMoGKJhEREZEAVDSJiIiIBKCiSURERCQAFU0iIiIiAahoEhEREQlARZOIiIhIACqaRERERAIoSHUAmcjMDDgI2JnqWERERKRVugAbXAKre6toSsxBwLpUByEiIiIJOQRY39oXqWhKzE6AtWvXUlpayuzZs5kyZQqFhYWpjqtD1NXVKecckYt5K2flnM1yMe/InKurq+nfvz8kOFKkoqkNunbtSmlpKWVlZXTt2jWnfgCVc27IxbyVs3LOZrmYdzJz1kRwERERkQBUNImIiIgEoKJJREREJADNaWonDQ0N1NXVpTqMdlFXV0dBQQE1NTU0NDSkLI7CwkLy8/NT9v4iIpJbVDQlmXOOTz75hO3bt6c6lHbjnKNv376sXbsWb8mq1OnevTt9+/ZNeRwiIpL9VDQlWbhg6t27N2VlZVn5y7yxsZFdu3bRuXNn8vJSM8LrnKOqqopNmzYB0K9fv5TEISIiuUNFUxI1NDQ0FUy9evVKdTjtprGxkdraWkpKSlJWNAGUlpYCsGnTJnr37q2hOhERaVeaCJ5E4TlMZWVlKY4kd4Q/62ydPyYiIulDRVM7yMYhuXSlz1pERDqKiiYRERGRAFQ0CQCTJk3iiiuu2Gf7E088od4cERHpMNW1DSxeu53q2tQtadMcTQSXtOeco6GhgYIC/biKiGSj6toGVlbsZECPMs56cB6rtlQxuFcZz1w+gdKi9DnJRz1NEtiMGTMYNWoUDz30EEcffTSdO3fmK1/5StSaVNOnT+f000/n5ptvpnfv3nTt2pULLriA2trapjbOOe644w4OPfRQSktLGTlyJI8//njT/rlz52JmPPfcc4wdO5bi4mJeeeWVjkxVREQ6SHVtA6fe+zKn3z+P0+57mVVbqgBYtaWKlRU7UxxdNP3pLq3ywQcf8Pjjj/PnP/+ZhoYGzjvvPC655BL+9Kc/NbV54YUXKCkp4cUXX2T16tWcc845HHDAAfz4xz8G4Ac/+AF/+9vfePDBBzn88MN5+eWX+eY3v8mBBx7IxIkTm45z3XXXceedd3LooYfSvXv3jk5VREQ6wDvrtjcVShsq99CvWwkbd9QwuFcZQ/t0SXF00VQ0pbFwd+XQPl3SpnuypqaG3/3ud3Tt2pWuXbvyi1/8gi984Qvcdddd9O3bF4CioiIeeeQRysrKOProo7nlllu49tprufXWW6murubuu+/m3//+N+PHjwfg0EMP5dVXX+Whhx6KKppuueUWJk+enJI8RUSk/VXXNnDdX99pej6wZyl/v/gE1myrSqvffWEZUzSZ2feBM4FhQDUwH/iec25lRBsDbgLOB3oArwOXOOf+E9GmGLgT+DpQCrwAXOycW9dBqQQS7q5Mt3HdAQMGcMghh1BZWQnA+PHjaWxsZOXKlU1F08iRI6PWqho/fjy7du1i7dq1bNq0iZqamn2KodraWkaPHh21bezYse2cjYiIpNLKip187PcyAdzx5ZH07FxEz85FKYyqeRlTNAETgfuBN/Hi/jEw28yOcs7t9ttcB1wFTAfeA34AzDGzoc658MDoPcB/AV8DtgB3Af80szHOubSZqr+yYuc+47qj+ndvt/fr2rUrO3bs2Gf79u3b6dq1a7OvC59ZF+QMOzOjsbERgKeffpqDDz44an9xcXHU806dOrV4TBERyVxD+3RhcK+ypg6CYw7pnuqQ9itjiibn3CmRz83sHGATMAZ42e9lugL4sXPub36bs4EKYBrwkJl1A74NfMs597zf5pvAWuBk4LmOyaZlsT9I7T2uO2zYMJ599tl9tr/55psMHTq06fmaNWvYsGEDnTt3BmDBggXk5eVxxBFHNLVZsmQJ1dXVTZc5ee211+jcuTOHHHIIPXr0oLi4mDVr1kQNxYmISO4pLcrnmcsnpN1UlOZkTNEURzf/fqt/PxjoC8wON3DO7TGzl4By4CG8Aqswps0GM1vmt4lbNPlDepHdIF3Au3RH+DT4uro6GhoacM7R2NjY1KOSqOIC45+XnsB7FTs5ok8XiguszcfcnwsvvJBf/vKXXHzxxZx33nmUlpby/PPP8/DDD/P73/+exsZGnHOUlJQwffp0brrpJhoaGrjsssv4yle+Qu/evZva1NbWcu6553LjjTfy8ccfc9NNN3HJJZcAXu/R1VdfzZVXXkl9fT0nnHAClZWVLFiwgE6dOnH22Wc35Rnkcwy/Z11dXbteey58mZZcu1xLLuatnHNDLuYMHZd3dW0D723axRG9O7dYCBUYHN23E9BIXV3yf89F5tzWvDOyaPJ7le4GXnXOLfM39/XvK2KaVwADI9rUOue2xWnTl+Z9H2+uVJTZs2c3zd2ZM2cOBQUF9O3bl127dkWdYt8Wg7vlUVezm7qapByuWT179uSZZ57hRz/6EVOnTqWmpoYhQ4Zw//33M3XqVCorK9mzZw+DBw/m85//PF/96lfZtm0bkydP5qc//WnTHKe6ujomTJjAgAEDmDhxIrW1tZx55plceeWVTW2uueYaunbtym233cbq1avp1q0bI0eObGpTVeUNS+7cubPFCwLX1tZSXV3Nyy+/TH19fft+SHjfcy7KxbyVc27IxZyhffOubYA73snn0xrjwBLHdcc0kA4dSHPmzGn6/ZIoc84lKZyOY2b3A18ATghP4DazcmAecJBzbmNE298A/Z1zp5jZNOB3zrnimOPNAT50zl3YzPvF62lat3nzZkpLS5kzZw6TJ0+moaGBtWvXMmjQIEpKSpKZclq4+eabefLJJ1m4cCE7d+6kS5cu+8xlOuecc9i+fTt///vfOySmmpoaVq9eTf/+/dv1M6+rq2v6ngsLC9vtfdJNLuatnJVzNuuIvJes28GXH3q96fnjFxzHyEO67ecV7Ssy5+rqag444ACAbs65ytYeK+N6mszsF8AXgQkxZ7x94t/3BTZGbO/N3t6nT4AiM+sR09vUG+9svLicc3uAPRExAFBYWNj0Q1dYWEheXh5mRl5eXos9JJkodtJ3ONfYNvG2t5fwZx75XbSnjnqfdJOLeSvn3JCLOUP75n30wT2i5uQefXAPCgtT39VUWFjY5hGJjCma/CG5XwBnAJOcc6timqzCK4omA4v81xThnXX3Pb/NQqDOb/N/fpt+wHC8M+9ERESkDTJtcndrZEzRhLfcwDTgS8BOMwvPQdrhnKt2zjkzuwe4wczeB94HbgCqgFkAzrkdZvYwcJeZbcGbRH4nsBR4vkOzyUAzZsxgxowZ+52cPXPmzI4LSERE0lJpUX67LpOTKplUNF3k38+N2X4OMNN/fAfegpUPsHdxyykRazQBXAnU4/U0hRe3nJ5OazSJiIhI+smYosk51+Lqic6b1T7DvzXXpga41L+JiIiIBJJ9s5VFRERE2oGKJhEREZEAVDSJiIiIBKCiSURERCQAFU0iIiIiAahoEgCmT5/O6aefnuowRERE0paKJkl7ybr4sYiISFuoaJIWvfTSS3zmM5+huLiYfv368f3vf7/p+j3/+Mc/6N69e9Mq4YsXL8bMuPbaa5tef8EFF/D1r3+96fn8+fOZMGECpaWl9O/fn8suu4zdu3c37R80aBA/+tGPmD59Ot26deO8887roExFRESap6JJ9mv9+vWceuqpjBs3jiVLlvDggw/yyCOPcOeddwIwYcIEdu7cyaJFiwCvwDrggAN46aWXmo4xd+5cJk6cCMDSpUuZOnUqZ555Ju+88w5/+ctfePXVV/nud78b9b4/+9nPGD58OAsXLuSHP/xhB2UrIiLSPBVN6ay2CtYt9O5T5IEHHqB///788pe/ZNiwYZx++unMmDGD+++/n8bGRrp168aoUaOYO3cu4BVIV155JUuWLGHnzp188sknvPfee0yaNAnwiqFp06ZxxRVXcPjhh1NeXs59993Ho48+Sk1NTdP7fvazn+Waa65hyJAhDBkyJAWZi4iIRFPRlK5qq+BXJ8BvP+vdp6hwevfddxk/fjxme69iU15ezq5du1i3bh0AkyZNYu7cuTjneOWVV/jSl77E8OHDefXVV3nxxRfp06cPw4YNA2DhwoXMnDmTzp07N92mTp1KY2Mjq1atanqPsWPHdmyiIiIiLciYa8/lnE3vwtYPvcdbP/SeHzKmw8NwzkUVTOFtQNP2SZMm8fDDD7NkyRLy8vI46qijmDhxIi+99BLbtm1rGpoDaGxs5IILLuCyyy7b570GDBjQ9LhTp07tkY6IiEjCVDSlq95HQs/DvIKp52He8xQ46qij+Otf/xpVPC1YsIAuXbpw8MEHA3vnNd1zzz1MnDgRM2PixIncdtttbNu2jcsvv7zpeMceeyz/+c9/NOQmIiIZR8Nz6aqoDC58Fb7zb+++qKzd33LHjh0sXrw46nb++eezdu1aLr30UlasWMGTTz7JjBkzuPjii8nL8358wvOa/vjHPzbNXZowYQJvv/121HwmgO9973ssWLCASy65hMWLF/P+++/z1FNPcemll7Z7fiIiIm2hnqZ0VlTWoUNyc+fOZfTo0VHbzj77bJ555hmuvfZaRo4cSc+ePTn33HO55pprotqddNJJvP32200FUo8ePTjqqKPYsGEDRx65t5fsmGOO4aWXXuLGG2/kxBNPxDnHYYcdxle/+tV2z09ERKQtVDQJADNnzmTmzJnN7n/jjTeaHjc2NlJZWRm1/84772xahiBs8eLFcY81btw4Zs+e3ex7rV69usV4RUREOpqG50REREQCUNEkIiIiEoCKJhEREZEAVDSJiIiIBKCiSURERCQAFU3tILxitrQ/fdYiItJRVDQlUWFhIQBVVam7wG6uCX/W4c9eRESkvWidpiTKz8+ne/fubNq0CYCysrJ9rtuWDRobG6mtraWmpqZpVfCO5pyjqqqKTZs20b17d/Lz81MSh4iI5A4VTUnWt29fgKbCKRs556iurqa0tDTlRWH37t2bPnMREZH2pKIpycyMfv360bt3b+rq6lIdTruoq6vj5ZdfZsKECSkdFissLFQPk4iIdBgVTe0kPz8/a3+h5+fnU19fT0lJieYSiYhIztBEcBEREZEAVDSJiIiIBKCiSURERCQAFU0iIiIiAahoEhEREQlARZOIiIhIACqaRERERAJQ0SQiIiISgIomERERkQBUNImIiIgEoKJJREREJAAVTSIiIiIBqGgSERERCUBFk4iIiEgAKppEREREAlDRJCIiIhKAiiYRERGRAFQ0iYiIiASgoklEREQkABVNIiIiIgGoaBIREREJQEWTiIiISAAFibzIzDoBE4EBQFHkPufcfUmIS0RERCSttLpoMrPRwDNAGdAJ2AocAFQBmwAVTSIiIpJ1Ehme+znwD6AnUA0cDwwEFgLXJC80ERERkfSRSNE0CrjLOdcANADFzrm1wHXAT5IYm4hIh6uubeDjnd69iEikRIqmOsD5jyvw5jUB7Ih4LCKScaprG/jiAwu4e1kBX3xggQonEYmSSNG0CBjrP34RuMXMvgHcAyxNUlwiIh1uZcVOVm+pAmD1lipWVuxMcUQikk4SKZpuADb6j38IbAEeBHoD5ycpLhGRDje0TxcG9SoDYFCvMob26ZLiiEQknbT67Dnn3FsRjz8FTk1qRCIiKVJalM9TF49n5t+fY/oZ4yktyk91SCKSRrS4pYhIhNKifAZ2QQWTiOwjUE+Tmb0NfM45t83MFrF3Ivg+nHPHJis4ERERkXQRdHjuSWCP//iJ9gmlZWY2AbgWGAP0A85wzj0Rsd+Am/DmVvUAXgcucc79J6JNMXAn8HWgFHgBuNg5t66D0hAREZEMFKhocs7dHO9xCnQClgC/A/4aZ/91wFXAdOA94AfAHDMb6pwLnwZzD/BfwNfwJrHfBfzTzMb4a0+JiIiI7COha8+linPuWeBZAK9TaS+/l+kK4MfOub/5287GW0tqGvCQmXUDvg18yzn3vN/mm8Ba4GTguQ5JRERERDJOoIngZrbNzLYGubV3wPsxGOgLzA5vcM7tAV4Cyv1NY4DCmDYbgGURbUREJINV1zaweO12LU4qSRe0p+mKiMe98Ia9ngMW+NvGA1OBW5MWWev19e8rYrZX4F0bL9ym1jm3LU6bvjTDnwdVHLGpC0BdXR0FBd5HWFdXl1jUGSicq3LOfrmYt3LObOFV3VdvqWJQrzKeujj+0hHZlHNr5GLekTm3NW9zrtkT4eK/wOyvwIvOuV/GbP8ucLJz7vQ2RRQ8DkfERHAzKwfmAQc55zZGtPsN0N85d4qZTQN+55wrjjnWHOBD59yFzbzXDLwJ5lFmzZpFWVlZkjISEZG2+ngn3L1sb3/AVcPrGag1SsVXVVXFtGnTALo55ypb+/pE5jRNBb4XZ/tzwE8TOF6yfOLf92XviuXgrVReEdGmyMx6xPQ29Qbm7+fYtwF3RzzvAqybMmUKpaWlzJkzh8mTJ1NYWNi2DDJEXV2dcs4RuZi3cs7snKtrG/jbxr09Tc0tUppNObdGLuYdmXN1dXWbjpVI0bQFOAP4Wcz20/19qbIKryiajHd9PMysCJjI3iJvId4FhycD/+e36QcMxzvzLi5/blR4yYWmSeiFhYVNP3SRj3OFcs4duZi3cs5MhYWFPHv5BFZW7GRony4tLlKaDTknIhfzLiwspL6+vk3HSKRougl42MwmsXdO0/HAKcB32hRNC8ysMzAkYtNgMxsFbHXOrTGze4AbzOx94H286+RVAbMAnHM7zOxh4C4z2wJsxVuzaSnwfHvGLiIi7W/rrlqeX1HBycP6aFV3SbpErj0308zeBS4DzgQMWA6EnHOvJzm+WGOBFyOeh4fMfo+3NtMdeAtWPsDexS2nRKzRBHAlUI/X0xRe3HK61mgSEclsW3fV8pmfPE99o6Mgz3jjhpPp2bko1WFJFklonSa/OPpGkmMJ8r5z8Yq05vY7YIZ/a65NDXCpfxMRkSzx/IoK6hu9k5vqGx3Pr6jgv8f2T3FUkk0SumCvmR1mZj8ys1lm1tvfdoqZHZ3c8ERERII5eVgfCvK8v6sL8oyTh/VJcUSSbVrd02RmE/FW5Z4HTMBbs2kTcAzenKYvJzNAERGRIHp2LuKNG05umtOkoTlJtkR6mn4K/MA5Nxmojdj+It4ilyIiIinRs3MR/z22vwomaReJFE0jgL/H2f4p3mrhIiIiIlknkaJpO9AvzvbRwPo2RSMiIiKSphIpmmYBt5tZX8ABeWYWwlvv6NFkBiciIiKSLhIpmm4E1uD1KnXGW6PpZbzLkPwoeaGJiIiIpI9EFresA75hZv+DNySXByxyzr2f7OBERERE0kVCi1sCOOc+BD5MYiwiIiIiaSuRdZoMby2mk4DexAzxOefOTE5oIiIiIukjkTlN9wJ/AAYDu4AdMTcREZG0Vl3bwMc7vXuRoBIZnvsmcKZz7plkByMiItLeqmsb+OIDC1i9pYC/bVzAs5dPoLQoP9VhSQZIpKdpB/BRsgMRERHpCCsrdrJ6SxUAq7dUsbJiZ4ojkkyRSNE0A7jJzEqTHIuIiEhSVNc2sHjt9rjDb0P7dGFQrzIABvUqY2ifLh0dnmSoRIbnHgO+Dmwys9VAXeRO59yxSYhLREQkIdW1DZx678us2lLF4F5lPBMz/FZalM9TF49n5t+fY/oZ4zU0J4ElUjTNBMYAfwQq8FYFFxERSQsrK3ayyh9+W+UPv43q3z2qTWlRPgO7oIJJWiWRoukLwFTn3KvJDkZERKSthvbpwuBeZU09TRp+k2RJpGhaC1QmOxAREZFkKC3K55nLJ7CyYidD+3RRb5IkTSITwa8G7jCzQUmORUREJClKi/IZ1b+7CiZJqkR6mv4IlAEfmlkV+04E75mMwERERETSSSJF0xXJDkJEREQk3bW6aHLO/b49AhERyQbVtQ2aSyOSpRLpaRIRkThaWh9IRDJbIhPBRUQkjnjrA4lI9lDRJCKSJOH1gQCtDySShTQ8JyKSJFofSCS7qWgSEUmi8PpAIpJ9Wl00mVkJcClwEtCbmCE+XbBXREREslEiPU2PAJOBx4E30AV7RUREJAckesHeU51z85IdjIhIqlXXNvDxTu++sLAw1eGISBpJ5Oy59YDOoxWRrFNd28Bp98/n7mUFnHb/fKprG+K2Wbx2e9x9IpLdEr1g7+1mNjDZwYiIpNI767azZms1AGu2VvPOuu1N+6prG3j9oy2ccu/LnH7/PE699+XAhZMKLZHskMjw3FtACfCRLtgrIrkgcqXvsPDilS2dKadVwkWyRyJF0/8CBwM3ABVoIriIZIljDunOgJ6lrNlazYCepRxzSHcgeqXvsKCLV8ZbJVxLEohkpkSKpnJgvHNuSbKDERFJpdKifP55STkz//4c088ob+oRCq/0vWpLFQN7lnLHl0dyzCHdA/UYRb5Wq4SLZLZEiqYVQGmyAxERSQelRfkM7EJUQdSWlb61SrhI9khkIvj1wF1mNsnMeplZ18hbsgMUEelIG7ZX89TH8Py7FVETt8MrfSdS9LTltSKSPhLpafqXf/9CzHbDm9+k/xVEJCOt31bNxLteAQp4YdYSDu6+guevOknFjogAiRVNJyU9ChGRNPCbVz6Ier5++x5N3BZpB+vXr2fevHlNt4ULF0btz8vLo6Eh/ZboaHXR5Jx7qT0CERFJtZqa6P+ki/NhQI+yqG3VtQ2an5QFtPJ7+9m5cyevv/56VFFUVVUVt+2RRx5JKBTiu9/9LqFQiCFDhmBmHRxxcIn0NGFm3YFvA0fiDcktBx5xzu1IXmgiIh2ra6foX557GuCsB+fx14tCrNlWxYAeZZz14DytuZThqmsb+OIDC1i9pYC/bVzAs/oeA6uvr2fp0qVRBdHatWvjtu3Tpw+hUIhbb72VUCjE6NGjKSoq6uCIk6vVRZOZjQWeA6rxLthrwFXAjWY2xTn3dnJDFBHpGJ9U7tln26otVZx238tsqNxDv24lbNxR07RdQ3eZaWXFTlb7a2et1vfYxDnHmjVrogqiJUviry5UUFBAKBTiW9/6FqFQiPHjx9OjR48OjrjjJdLT9HPgKeA851w9gJkVAL8F7gEmJC06EZEOdMGEw3hqycaobfl5xga/mNq4o4aDuhazoXKP1lzKYEP7dGFQrzJWb6liUA59jzt27ODVV19l1qxZ3HvvvcyfP589e/b9QwFgxIgRhEIhrrnmGkKhEIMGDUrrYbOOkkjRNJaIggnAOVdvZnfgXWJFRCQj7dpTv8+2hkYXVSiFh+o0pylzlRbl89TF4/1FTMdnxfdYV1fH4sWLo3qJNm7cGLdtz549Oemkk7j99tsJhUKMHDlS87oCSqRoqgQG4C1yGak/sLPNEYmIpEhl9b5FU0Ge8dhFIT7dtaepUOrZObPnZUj8RUzTlXOOjz76qKkYmj9/PsuWLYvbtqSkhFAoxHnnnUcoFOK4446jW7duTfvr6up45plnOPXUU1UoJSCRoukvwMNmdg0wH28i+AnAz/CuSycikpEqdtbss62+0fHprj2a8yLtZuvWrSxYsCCql6i50+1Hjx7NpEmTuOGGGwiFQgwYMKCDo81tiRRN1+AVSo9GvL4OeBBvtXARkYxUkr/vnI2BOTTnRZJvz549vP3228yfP7+pINq0aVPctgMGDCAUCnHvvfcSCoUYMWIE+fnp3xOWSxJZp6kWuNzMvg8chnf23AfOufiLMIiIZIif/mvlPtvuOOuYjBjCkY7nnOP999+P6iFasSJ25oqnc+fOUesRfeYzn6Fz584dHLG0VULrNAH4RdLSJMYiIpJSI/t35YUVm6O2dS5O+L9JyXCffvop8+fPj+olas7YsWOZOnUqt9xyC+Xl5Rx88MEdGKl0lED/G5jZ34Ie0Dl3ZuLhiIikzjsfb95n29d/8xqv33Byq3qbWlo1XKuKp15NTQ3Lly/nP//5D6+99hrz5s1j69atcdseeuihUWsSHXXUUeTlJXK9e8l0Qf+Eilzp24Az/G3hJQbGAN2BwMWViEi6+bR6322VNfW8s247xx3aK9AxqmsbOPXel5tdNTxyf79uJTx96Yk6Gy/JGhsbWbFiRdSw2QcffBC3bbdu3QiFQlx99dWEQiHGjRtHWVlZ3LYigYom59w54cdmdjvwf8CFzrkGf1s+8ADecgQiIjlrZcVOVvmrTcdbNTxy/8YdNZx238u8cM1J6nFqhYqKiqiC6PXXX2+27fHHH8+XvvSlplWr+/btC+jUe0lMIoP15wInhAsmAOdcg5ndjbcEwbXJCk5EJNX69yjhmEO6Bx5SG9qnC4N7lTX1NMWeeTe0T5eoy7FsqNyjy3hEqKqq4o033ogqiior4/89fvjhhxMKhfjOd75DKBRi6NChGjaTdpVI0VSAd6He2NNMjgT00yoiGWtE3y4s/SR6jd6SwnwWfbyVG59YFuhCvaVF+Txz+YRmC6zSonyevvTEpuvZ5dLlWBoaGli+fHlUQbRq1aq4bXv16kUoFOLGG28kFAoxZswYSkpKOjhikWiJFE2/Ax4xsyHAa/624/HWaPpdsgITEelo2/fU7rPt/U27mfbwG03Pg1yot7Qof7/7e3Yu4oVrTsq6yeDr16+PKogWLlzYbNtQKMRXvvKVpmGzAw88sAMjFUlMootbfgJcCfTzt20E7gDuSlJcIiIdrm/XEtZui38B06J8qG0gqmcol86C27VrF6+//npUUbR79+64bY888sioNYmGDBmii71KVkhkcctGvALpDjPr6m/TBHARyXilBc0XPjOnf4ayksKmAqmls+T2py2vbQ/19fV89NFHPPDAA02n369duzZu2z59+hAKhbj55psJhUIce+yxFBXp7D/JDW1atU3Fkohkk3fW79hn20HdShg7qBsDD+hMz05FTT1LLZ0ltz9teW0i1qxZE9VDtHjx4rjtCgoKotYjOv744+nZs2e7xSWSaVpdNJlZH+BO4HNAb7x1m5o457K7j1pEsla/rsVsr4m+ItSGHTU8taSGp5ZUcHD3YtZv9yZv//Wi0H7Pktufls6wa40dO3Y09Q6Fb3v2xB9iHDFiBKFQiKuuuopQKMTgwYMxM51+LxJQIj1NM4EBwK14c5lcMgMSEUmV44ccwLub1jS7f/12rxhZtaWKNduq9nuW3P60dIZdWF1dHYsXL466jMeGDRvitj3ooIMIhULcfvvthEIhRo4cqQJIJMkSKZpOAE50zi1OciwdyswuxltTqh/wH+AK59wrqY1KRFJpZ3VDy43YOxm8pbPk4glPHj+id2e61m3l8b/8o6kgWrZsWdzXlJSURK1HdNxxx9GtW7dWva+ItF0iRdNaYobkMo2ZfRW4B7gYmAdcADxrZkc555r/M1NEstrLK9a32Oa7nx3CJZOGNPUOrd9Wzf++uYavjxvAwT1Km9pt3bqVBQsWNBVE8+fPp76+Pu4xR40axcSJE7nhhhsIhUIMGDAgOQmJSFIlUjRdAfzUzC5wzq1Objgd5irgYefcb/3nV5jZVOAi4PupC0tEUmlTNSxgGgcWwKf1MJ5ZUfsL8oxzyweTTwOvvfYmz8yZy51/fIo9697l2qrtcY85YMAAQqEQ99xzD70OHcH3XtyO5XkF1xOXhLQSuEgGSaRo+gtQBnxoZlVAXeRO51xan2phZkV4Fxj+acyu2UB5M68pBoojNnUBb75BQYH3EdbV1cV5ZXYK56qcs19O5V1XxVs2jR4FjXywrZF5axspX30Gj63rQ/3WdU3Net229yVFJaXk9R1Gl2O/QPHBR/Ldr0zh+i+OavYtqmsb+NWyBazeUsWgXmUc2rMkLT7bnPqefbmYM+Rm3pE5tzVvc65187jN7Oz97XfO/b5NEbUzMzsIWA+EnHPzI7bfAJztnBsa5zUzgJtit8+aNUtXwxbJMJWVlaxYsYJ3332XFStWsGLFCpr7f3BMvzzeO+g0egw8ku+fcjj9eh8QtX9rDdy8KB9vxoLjptEN9GzhSh+1DbCxCvqVeQtmikjHqaqqYtq0aQDdElk2qdVFU6aLKJrKnXMLIrbfCHzLOTcszmvi9TSt27x5M6WlpcyZM4fJkyfnzJkqdXV1yjlHZFreNTU1vP3228yfP5/58+ezYMECtmzZErft4MGDGT9+POXHDOHE1Xdx9IF55OdFT9esr4ezDnySP3zn+GbPcNuwvZq/vLWOr449hIO6l8Ztk+4y7XtOhlzMGXIz78icq6urOeCAAyDBoqlNi1uaWSkQ9alnwIKXm4EGoG/M9t5ARbwXOOf2AE0Ln4QvB1BYWNj0Qxf5OFco59yRDnk3NjaycuXKqNPv33vvvbhtu3btSigU4sorryQUCjFu3Dg6deq0b8PaKnjsHKiOXxB9WtCDJy+ftN+4Bh5YyHWfP6q16aSldPieO1ou5gy5mXdhYWGzJ2MElcjilp2A24H/BnrFaZLWHc7OuVozWwhMBv4esWsy8GRqohKRioqKqILotddea7btcccdx2mnnUYoFKK8vJy+fWP/Bgqgtgp+dQJs/bDZJv16H9z644pI1kqkp+kO4CS80/UfBS4BDsY7bf/65IXWru4G/mBmbwELgPPxFuz8VUqjEslSVVVVvPnmm1Gn32/fvj1u28MPP5xQKMS5555LKBRi2LBh5OXlJT+oTe9GF0wDT4CPX2166gA75bZ9XyciOSuRoum/gP/nnJtrZo8ArzjnPjCzj4FvAH9KaoTtwDn3FzPrBfwP3uKWy4BTnXMfpzYykczT2NjI8uXLoy7j8dFHH8Vt27NnT0KhENdffz2hUIgxY8ZQWpqieUC9j4Seh3mFU5d+sGZB1G4DqIo/H0pEclMiRVNPYJX/uNJ/DvAq8GAyguoIzrkHgAdSHYdIutu6dSuPP/44r7/+OvPmzeOtt95qtm15eTlnnXUWoVCI8ePH07t37w6MtJWKyuDbs+GhCVDZzKKWy/4Kw8/o2LhEJG0lUjR9BAwCPgaW481tegOvB2p7sgITkfa3e/fupmIofNu1a1fctsOGDSMUCnHxxRcTCoU4/PDDm06KyFjbPm6+YAIo7dFxsYhI2kukaPodMBJ4CbgNeNrMLvWPdVUSYxORNmhoaGDp0qVRBdGaNfGvEtS7d29CoRAzZsygvLycY489luLiYurq6njmmWc49dRTg51pU1vlzRXqfaTXk5PuIofoug+EASF4J2IV8IGh1MUmImmn1UWTc+7nEY9fNLNhwFjgQ+fckmQGJyLNW7NmTdOk6nnz5rFo0aK47fLz8wmFQnzjG98gFApx/PHH06tXvBNf2yjybLSeh8GFr6Z/4VRU5sUZLvTm3Re9f9vqlIQlIumpVUWTmRXiXW7kAufcewD+BW51kVuRJKqsrOS1116L6iWqqamJ23b48OFRaxINHjw4NcNmkWejbf3Qe37ImI6Po7WKyvbGuTHm777Y5yKS01pVNDnn6sxsON7ZuCKSgPr6epYsWRJVEK1fH39eTb9+/QiFQtx2222EQiFGjRqVvgvSRQ519TzMe55pOvXd/3MRyWmJzGl6FPg2mbMmk0iHcc6xevXqqIJo6dKlcdsWFxdHrUd03HHH0b17944NOJlih7rSfWgunsM/i1v0iH8lObAjP5/qiEQkjSRSNBUB3zGzycBbwO7Inc45TQaXrLVt2zYWLlwYVRQ1d9XskSNHcuKJJzatSTRgwIDMP9usJZFDXZmmtgr++h3C35ABzLsHjpiSuphEJK0kUjQNB972Hx8Rs0/DdpKxamtrWbRoUVRBVFER93KE9O/fn1AoxN13300oFGLEiBEUFLTpUo6SapvehYaYeWP5mXkBXhFpH4mcPXdSewQi0p6cc3zwwQdRl/FYvnx53LadOnUiFApx0UUXNQ2blZSUtO7Ue8k8vY+Ekh5Qs23vtv5jUxePiKQd/WksWWHz5s0sWLAgqpfIufgdn2PGjOHkk0/mpptuory8nEMOOaTF4zc3BCdk3tpMzSkqg1N+Ck9csHfbgcNSF4+IpB0VTZL2ampq9plHtGVL/GuCDRo0iFAoxLRp0wiFQhx99NHk5+d3cMQ5JBPXZtqfuqro59Xb4rcTkZykoklSyjnHypUrowqi9957L27brl27Rq1HNG7cODp16tTBEUuUTF2bqTmHT8ZZPuYacHkF2FFfTHVEIpJGVDRJu6qoqGhasXrevHm89tprzbY97rjjOO200ygvL6e8vJx+/fp1YKSSkGxYmymstgoe/RLmGqjJ70r+RS9T2OmAVEclImlERZMkrLq6ep95RNu3b4/bdsiQIVFrEg0bNoy8vLyODViSLxvWZgqL6DUraaikftcmOGBwioMSkXSSUNFkZkcAk4DeQNRvPufcLW0PS1KtsbGR5cuXRxVEH330Udy2PXv2JBQKcf3111NeXs7YsWMpLdWp2jkjE9ZmCjJZvcdA6HowVK6nqrAnhT0GdmyMIpL2Wl00mdl5wIPAZuATotdmcoCKpgywcePGqILozTffbLZteXk5Z511FqFQiPHjx9OjRw+dfi+ZI8hk9doqeHgKVK7HWT5ldVtxMz8PF83L7N4zEUmqRHqafgDc6Jy7PdnBSHLs3r2b119/Paoo2rVrV9y2w4YNi1qT6PDDD29x1Wqdfi8ZJchk9Yg25hq8+20fZf7EdhFJqkSKph7AY8kORIJpaGhg2bJlUQXRxx9/HLftgQceSCgU4qabbiIUCnHsscdSXFzcwRGLpFiQyeoRbRzeJVSc5WMaohORCIkUTY8BU4BfJTkWAdauXRtVEC1atChuu/z8/Kj1iI4//nh69erVwdGKZIAgk9WLyuDbs2HeL7D59wB+j9OnK0Fn0ImIL5Gi6QPgVjM7HlgKRI3VOOfuS0ZguWjRokUce+yxABx99NGEQiGuuOIKQqEQhx56aPZf7FWkvbQ0WT08pyk8jBdWs6N94xKRjJJI0XQ+sAuY6N8iOUBFU4JGjx7d7KU/RCTJIs+oi5z3FKnyk46PS0TSViIX7NXCJSKS2WLPqPv27H3nNAF2xOQUByoi6USrC4pI7ok9o27bx968py/eT3gQ3AB2bUpRgCKSjgL1NJnZ3cAPnXO7/cfNcs5dlZTIRETaS7wz6orKYPiZuJd/hm1fjes+CMvky8KISNIFHZ4bDRRGPG6OJuSISPrb3xl14RMudOKFiMQIVDQ5506K91hEJGPFO6NuwyJs2yoA716LW4pIBM1pEhEBb3L4Exc1PXV5Bd716EREfCqaRETA61Xavnd1fWus9yaIi4j4VDSJiIA3t6nH3hVVXI/B8S+5IiI5K5HFLUVEsk9RGVw0n/o1b/La669x3BkXUhjvkisikrPU0yQiElZUhhtYzpbOw6BQBZOIRGt10WRmZ5vZFyKe32Fm281svplp1qSIZLa6Krrv/hDqqlIdiYikmUR6mm4AqgHMbDzwXeA6YDPw8+SFJiLSwWqrKPj1BCa+dzMFv57gnVEnIuJLpGjqD3zgPz4deNw592vg+8CJSYpLRKTjbViEbV8N4N1vWJTScEQkvSRSNO0CevmPpwDP+49rgNJkBCUiIiKSbhIpmuYAvzWz3wJHAE/7248GVicpLhGRjnfQaG+pAfwlBw7a31WjRCTXJFI0XQIsAA4EznLObfG3jwH+N1mBiYh0uKIy6s97iZeOuIn6816KviadiOS8Vq/T5Jzbjjf5O3b7TckISEQkpQrL2N7pMC05ICL7CFQ0mdkxQQ/onHsn8XBERERE0lPQnqbFgAPMv9+f/LYEJCIiIpKOgs5pGgwc6t+fBawCLgZG+7eLgQ/9fSIiIiJZJ1BPk3Ou6VLfZvYYcJlz7pmIJu+Y2VrgVuCJpEYoIiIikgYSOXtuBF5PU6xVwFFtC0dEREQkPSVSNL0L/MDMSsIbzKwY+IG/T0RERCTrtHrJAeBC4B/AWjNb4m8biTdB/LRkBSYiIiKSThJZp+kNMxsMfBMYhndG3V+AWc653UmOT0RERCQtJNLThHOuCvh1kmMRERERSVsJFU1mdgQwCehNzLwo59wtbQ9LREREJL20umgys/OAB4HNwCdEL3bpABVNIiIiknUS6Wn6AXCjc+72ZAcjIiIikq4SWXKgB/BYsgMRERERSWeJFE2PAVOSHYiIiIhIOktkeO4D4FYzOx5YCtRF7nTO3ZeMwERERETSSSJF0/nALmCif4vkABVNIiIiknUSWdxycHsEIiIiIpLOEpnTJCIiIpJzElmn6ZH97XfOnZt4OCIiIiLpKZE5TT1inhcCw4HuwL/bGpCIiIhIOkpkTtMZsdvMLA94APgoGUGJiIiIpJukzGlyzjUCPweuTMbx4jGzG81svplVmdn2ZtoMMLN/mNluM9tsZveZWVFMmxFm9pKZVZvZejP7HzOz9opbREREskNCF+xtxmFJPl6sIryFNRcA347daWb5wNPAp8AJQC/g94ABl/ptugJzgBeBccARwExgN3BXO8YuIiIiGS6RieB3x24C+gFfwCtS2oVz7ib//ac302QKcBTQ3zm3wW97NTDTzG50zlUC3wBKgOnOuT3AMjM7ArjKzO52zrlmji0iIiI5LpHhudExt2P87VcDVyQnrISMB5aFCybfc0AxMCaizUt+wRTZ5iBgUEcEKSIiIpkpkYngJ7VHIEnQF6iI3OCc22Zmtf6+cJvVMa+riNi3Kt6BzawYr/gK6wJQV1dHQYH3EdbV1cV5ZXYK56qcs18u5q2cc0Mu5gy5mXdkzm3N2xIdkTKzA4GheJdOec8592kCx5gB3NRCs3HOubciXjMduMc51z3mWL8GBjrnpsZsrwX+n3Puz2Y2G1jlnLsgYv/BwDpgvHPutdbEOWvWLMrKyloIX0RERNJBVVUV06ZNA+jmT9tplUTmNHUCfgH8P/YO7zWY2aPApc65qlYc7pfAn1toszrgsT4BjovcYGY98NaRqoho0zfmdb39+wqadxsQOZerC7BuypQplJaWMmfOHCZPnkxhYWHAUDNbXV2dcs4RuZi3clbO2SwX847Mubq6uk3HSuRst7vxLtT7X8A8f9sJeBfqvQu4KOiBnHObgc0JxBDPAuBGM+vnnNvob5sC7AEWRrT5iZkVOedqI9psYD/FmT8HqmkeVHiFgsLCwqYfusjHuUI5545czFs554ZczBlyM+/CwkLq6+vbdIxEJoKfBXzbOfesc67Svz0DnAd8uU3R7Ie/BtMoYACQb2aj/Ftnv8lsYDnwBzMbbWafA+4EfhPRBTcLr/iZaWbDzewM4AZAZ86JiIjIfiXS01RG/KGsTf6+9nILcHbE80X+/UnAXOdcg5l9AW9l8nlANV6RdE34Bc65HWY2GbgfeAvYhtdzFruMgoiIiEiURIqmBcDNZvb/nHM1AGZWijdRekEyg4vknJsOTG+hzRrgtBbaLAUmJC0wERERyQmJFE2XA/8C1pnZEryz50YBNcDU/bxOREREJGMlsk7TMjM7HPgmMAxvRfA/A39yzrVtWrqIiIhImkroWnF+cfSbJMciIiIikrYSWaepl3Nui/+4P95Zc6XAP5xzLyc5PhEREZG0EHjJATMbYWargU1mtsI//f9N4ErgfODfZnZ6ewQpIiIikmqtWafpDmAp3sKWc4F/As8A3YAewEPA9UmOT0RERCQttGZ4bhzwWefcO2a2GK936QHnXCOAmf0CiHvtNhEREZFM15qepp54127DObcL2A1sjdi/De+abCIiIiJZp7WXUYm91IguPSIiIiI5obVnz800s/CFa0uAX5nZbv95cfLCEhEREUkvrSmafh/z/I9x2jzahlhERERE0lbgosk5d057BiIiIiKSzlo7p0lEREQkJ6loEhEREQlARZOIiIhIACqaRERERAJQ0SQiIiISgIomERERkQBUNImIiIgEoKJJREREJAAVTSIiIiIBqGgSERERCUBFk4iIiEgAKppEREREAlDRJCIiIhKAiiYRERGRAFQ0iYiIiASgoklEREQkABVNIiIiIgGoaBIREREJQEWTiIiISAAqmkREREQCUNEkIiIiEoCKJhEREZEAVDSJiIiIBKCiSURERCQAFU0iIiIiAahoEhEREQlARZOIiIhIACqaRERERAJQ0SQiIiISgIomERERkQBUNImIiIgEoKJJREREJAAVTSIiIiIBqGgSERERCUBFk4iIiEgAKppEREREAlDRJCIiIhKAiiYRERGRAFQ0iYiIiASgoklEREQkABVNIiIiIgGoaBIREREJQEWTiIiISAAqmkREREQCUNEkIiIiEoCKJhEREZEAVDSJiIiIBKCiSURERCSAjCiazGyQmT1sZqvMrNrMPjSzm82sKKbdADP7h5ntNrPNZnZfnDYjzOwl/zjrzex/zMw6NqM2qK2CdQu9exEREekwBakOIKBheAXeBcAHwHDgN0An4BoAM8sHngY+BU4AegG/Bwy41G/TFZgDvAiMA44AZgK7gbs6KpmE1VbBr06ArR9Cz8PgwlehqCzYa3dvhpX/gqGnQKcD2jdOERGRLJQRRZNz7l/AvyI2fWRmQ4GL8IsmYApwFNDfObcBwMyuBmaa2Y3OuUrgG0AJMN05twdYZmZHAFeZ2d3OOddBKSVm07tewQTe/aZ34ZAxLb9u92a4axg01kFeIVy9QoWTiIhIK2VE0dSMbsDWiOfjgWXhgsn3HFAMjMHrXRoPvOQXTJFtbgMGAavivZGZFfvHCesCUFdXR0GB9xHW1dW1IZWAeg6hoMeh2LaPcD0Opb7nEAjwvrb8aQoa/XaNddQvfxo3alrCYYRz7ZCc00Qu5gy5mbdyzg25mDPkZt6RObc1b0v3zpV4zOww4G3gaufcb/1tvwYGOeemxLTdg9ez9L9mNhtY7Zw7P2L/QcB6oNw5t6CZ95sB3BS7fdasWZSVBRweS5L8xj10qV7HztJDaMgrbvkFQGF9JacsvZw8Gmgkn3+NuJe6gq7tHKmIiEh6qaqqYtq0aQDd/BGoVklpT1NzxUiMcc65tyJecxDeUN1j4YIpQrwK0GK2x7axZrZHug24O+J5F2DdlClTKC0tZc6cOUyePJnCwsL95ZFSDSdPpvG953BHTGVyWa82Hauuri4jck6mXMwZcjNv5aycs1ku5h2Zc3V1dZuOlerhuV8Cf26hzerwA79gehFYAJwf0+4T4LjIDWbWAygEKiLa9I15XW//voJm+MN5TUN64ZPtCgsLm37oIh+npW59YdzZST1k2ufcDnIxZ8jNvJVzbsjFnCE38y4sLKS+vr5Nx0hp0eSc2wxsDtLWzA7GK5gWAuc45xpjmiwAbjSzfs65jf62KXjFzsKINj8xsyLnXG1Emw1EFGciIiIisTJlnaaDgLnAWryz5Q40s75mFtlrNBtYDvzBzEab2eeAO4HfRIxbzsIromaa2XAzOwO4AUi/M+e0HpOIiEhaSfXwXFBTgCH+bV3MPgNwzjWY2ReAB4B5QDVekRRekgDn3A4zmwzcD7wFbMObq3Q36aQt6zGJiIhIu8iIosk5NxNvEcqW2q0BTmuhzVJgQlICS7baKm/tpfqaxNZjEhERkXaTEUVTTojsXeoxGHocCts+8nqaeh+Z6uhERERynoqmdBG52ve2VTD9GSgo8QomDc2JiIiknIqmdNH7SK9XKTyP6aDRKpZERETSiIqmdFFU5k343vSuepdERETSkIqmdFJUpgnfIiIiaSoj1mkSERERSTUVTSIiIiIBqGgSERERCUBFk4iIiEgAKppEREREAlDRJCIiIhKAiiYRERGRAFQ0iYiIiASgoklEREQkABVNIiIiIgHoMiptUFlZSV1dHVVVVVRWVlJYWJjqkDqEcs6NnCE381bOyjmb5WLekTlXV1e36VjmnEtSWLnDzA4G1qU6DhEREUnIIc659a19kYqmBJiZAQcBO4EueAXUIf7zXKCcc0cu5q2cc0Mu5gy5mXdszl2ADS6BAkjDcwnwP+j1AF79BMBO51xlyoLqQMo5N3KG3MxbOSvnbJaLecfJOeG8NRFcREREJAAVTSIiIiIBqGhquz3Azf59rlDOuSMX81bOuSEXc4bczDtpOWsiuIiIiEgA6mkSERERCUBFk4iIiEgAKppEREREAlDR1EZmVmxmi83MmdmomH0DzOwfZrbbzDab2X1mVpSiUJPCzJ4yszVmVmNmG83sD2Z2UEybrMnbzAaZ2cNmtsrMqs3sQzO7OTafbMoZwMxuNLP5ZlZlZtubaZNVOQOY2cX+d11jZgvN7MRUx5RMZjbB/842+P9nnR6z38xshr+/2szmmtnRKQq3zczs+2b2ppntNLNNZvaEmQ2NaZNVOQOY2UVm9o6ZVfq3BWb2+Yj9WZdzLP+7d2Z2T8S2Nuetoqnt7gA2xG40s3zgaaATcALwNeAs4K4OjS75XgT+GxiKl89hwOPhnVmY9zC8fycXAEcDVwIXAj8JN8jCnAGKgMeAB+PtzMaczeyrwD3Aj4HRwCvAs2Y2IJVxJVknYAnw3Wb2Xwdc5e8fB3wCzDGzLh0TXtJNBO4Hjgcm4y3oPNvMOkW0ybacwVv9+npgrH/7N/BkRIGQjTk3MbNxwPnAOzG72p63c063BG/A54F3gaMAB4yK2dcAHBSx7WtADdA11bEn8TP4ItAIFOZK3sC1wEe58F0D04HtcbZnXc7A68CDMdveBW5LdWztlK8DTo94bsBG4HsR24qB7cAFqY43STkf6Oc9IVdyjshrK/DtbM8Z6Ay8B5wMzAXuSeZ3rZ6mBJlZH+A3wLeAqjhNxgPLnHORvVDP4X1JY9o/wvZnZj2BbwDznXN1/uaszxvohvcfUFgu5Bwrq3L2hxXHALNjds0Gyjs+opQYDPQl4jNwzu0BXiJ7PoNu/n3432/W52xm+Wb2NbxexgVkf873A087556P2Z6UvFU0JcC8C9nMBH7lnHurmWZ9gYrIDc65bUCtvy9jmdntZrYb2AIMAL4UsTtr8wYws8OAS4FfRWzO6pybkW05HwDkE5OT/zwT80lEOM+s/Az8/7fvBl51zi3zN2dtzmY2wsx24S3o+CvgDOfccrI7568BxwLfj7M7KXmraIrgTxBzLdzG4v3S7Arc1sIh460cas1sT5lW5B32M7w5H1Pwhmge9f9DCkv7vBPIGfMmvP8LeMw599uYQ2Zlzi1I+5wTEBt7pueTiGz9DH4JHAN8Pc6+bMx5JTAKbz7Xg8DvzeyoiP1ZlbOZ9QfuBb7pnKvZT9M25V2QQGzZ7JfAn1tosxr4Ad4P4p7oWoG3zOxPzrmz8SaYHRe508x6AIXsW+mmWtC8AXDObQY2A++Z2bvAWrzPYwGZk3ercvYLphfxcjw/pl1W5tyCTMk5qM14fwDE/sXZm8zMJxGf+Pd98eZ+hGX8Z2Bmv8CbfznBObcuYlfW5uycqwU+8J++5U+Ovhy43d+WbTmPwcthYcTv5Xxggpl9F+/kJWhj3iqaIkQUA/tlZpfhFU5hB+HN5/gq3mRS8H653mhm/Zxz4S9oCl5X6cKkBZ0EQfNuRvins9i/z4i8W5OzmR2MVzAtBM5xzjXGNMm6nAPIiJyDcs7VmtlCvDOs/h6xazLwZGqi6nCr8IqIycAiaJrrNRH4XgrjSpjfA/4L4AxgknNuVUyTrMt5Pwzv/+lszfkFYETMtt8BK/AKxY9IRt6pnumeDTdgEPuePZcPLAWexxvK+hxej8wvUh1vG/L8DN6pmqOAgcBJeKdlfwAUZ2PeeAXx+/4/yIPx/krpC/TN8u96gP89/w+w0388CuicxTl/FW9O1rnAkcDPgV3AwFTHlsQcO0d8lw5vCY1RwAB///fwziY6AxgOzMJbUqVLqmNPMN8H/HwmRv7bBUoj2mRVzn5OPwFO9H83jcBbRqMBmJytOTfzOczFP3suWXmnPKlsuBGnaPK3DwD+iXd23Ra8v3iKUx1vG/Icgbfexxa8U8tX4Y2VH5yteeOdcu/i3bI1Zz+fmc3kPSlbc/ZzuhhviDLcYzYh1TElOb9JzXyvM/39BszAG76owTuzaHiq425DvnH/7QLTI9pkVc5+Tg9H/BxvwvvjZnI259zM5xBbNLU5b/MPJCIiIiL7obPnRERERAJQ0SQiIiISgIomERERkQBUNImIiIgEoKJJREREJAAVTSIiIiIBqGgSERERCUBFk4iIiEgAKppEpFXMbJCZOTMbleDr55rZPcmNKnXaM590+qzMbIaZLY54PtPMnkhdRCIdTxfsFUkzZjYTODti01bgTeA659w7KQkquc4E6trr4P7n1905d3p7vYcAcDl7L9gtkhPU0ySSnv4F9PNvnwPq8a7zlvGcc1udcztTHUe6ME/G/QHrnNvhnNue6jhEOpKKJpH0tMc594l/WwzcDvQ3swPDDcxshJn928yqzWyLmf3azDpH7N9naMfMnvB7YsLPV5vZDWb2iJntNLM1ZnZ+zGs+Y2aLzKzGzN4CRsfs72FmfzKzT/1Y3jezc5pLLDauIDHEOcaXzWxpRO7Pm1knM5uB10v3JX8I0ZnZJP81t5vZe2ZWZWYfmdmtZlYYccwZZrbYzL7lx7TDzP5sZl0i2nQys0fNbJeZbTSzq+PE9k0ze8vP5RMzm2VmvSP2T/Ljmup/nnuAE4McO857hWM+1//cdpnZg2aWb2bX+e+/ycxujHldN//nZZOZVfo/RyNj2lxvZhV+Hg8DJTH7o4bn/M/sipg2i/3vJPzcmdkFZvZP/3t418zGm9kQ/+dit5ktMLPDWspdJBVUNImkOb8Q+gbwAbDF31aG1xu1DRgHfAU4GfhlAm9xNRAuhh4AHjSzYf77dMLr4VoJjMG7QvidMa+/FTgK+DxwJHARsDlZMcQys37A/wKP+O83Cfgb3lDRncD/Ed1TN99/6U5guh/r5cB5wJUxhz8MOB04zb9NBK6P2P8z4CTgDGCK/95jYo5RBPwQGOkfazAwM04qdwDf93N4J+Cx4zkM77M/Bfg6cC7wNHCIH//3gB+Z2fHg9Wz5+/sCp/rv8Tbwgpn19Nv8N3AzcCMwFu+q8BcHiCWIHwKPAqOAFcAs4CHgNv+9ILGfY5H255zTTTfd0uiG9wu2Htjl3xywATg2os15eHOdOkVsOxVoAPr4z+cC98Qc+wlgZsTz1cAfIp4bUAFc6D8/H69QK4toc6Ef0yj/+VPAI63ILyqulmKI8/pj/fcfuJ/P74kAcVwLvBXxfAawG+gSse0O4DX/cWe8XqGvRuzvCVTFfs4x7zPOj7ez/3yS//xLEW0SPXa8mP8FrALyIratAK73H38W2AEUxxzrA+B8//F84MGY/a8Bi5v7nP3v8YqY1ywGZkQ8d8CtEc+P97edG7Hta0B1Kv7t6aZbSzf1NImkpxfx/hIfBRwHzAaeNbOB/v4jgSXOud0Rr5mH13s8tJXv1TS53DnngE+A8HBS+H2qItoviHn9g8DX/KGYO8ysvJXv31IMsZYALwBLzewxMzvPzHq09Ab+kN6r/pDVLrwesgExzVa76PlWGyPiOAyvF6kpf+fcVrxeuMj3GW1mT5rZx2a2E69IJM57vRXxONCxmxEbcwWw3DnXGLMtnMcYvCJtiz+ct8v/PAb7cYD3vcd+z7HPExV5MkOFf780ZluJmXVN0vuJJI2KJpH0tNs594F/ewP4NtAJr4cJvN4Y18xrw9sb2ffspkL2FXsmm2Pv/w0tnh3lnHsWGAjcAxyEN8wTO4TXkv3FEPt+DcBkvCGp5cClwEozG9zcwf2hqT8Dz+INu40GfoxXqASNo8XPwh/OnI3XQ/hNvF6mM/zdse8VWfC25Sy0eDHvL488vGJwVMxtKN4QYaIS+Xlz+9mm30+SdvRDKZIZHN4vpVL/+XJglP9LOizkt3nPf/4p3pweAMwsHxjeyvddDow0s9KIbcfvE5xznzrnZjrnvglcgTes126cZ55z7ia8AqiWvcVJLZAf85IQ8LFz7sfOubecc+/jFXqt8QHeL/em/P0eriMi2gwDDsAbCnvFObeC5nvMWnvsZHkbbz5TfURhHr6F56K9y77f8z7fe4zYn7eueL1XIllDRZNIeio2s77+7UjgF3hDKv/w9/8JqAF+b2bDzewkv80fnHPhIY9/A18wsy/4k6ofALq3Mo5ZeIXYw2Z2lJmdClwT2cDMbjGzL/lnQB2N15PzbqszDsjMjjPvbLuxZjYAb92nAyPeczVwjJkNNbMDzDtD7gNggJl9zcwOM7PL2FtkBeKc2wU8DPzMzD5nZsPx5vVEDoOtwSvaLjWzQ83si3gTn5Nx7GR5Hm+o7Qn/DL5BZlZuZj8ys/BE7HuBc/2z8o4ws5uBo1s47r+Bb5nZiX78v8ebYyeSNTJubRCRHHEK3hAKeGd9rQC+4pybC+CcqzKzqXi/3N7EmzD8V+CqiGM8gncG16N4E8t/jjdXKjDn3C4z+y/gV8AivJ6n7/nvFVaLd+bTIKAaeAVvMm97qQQm4PVodQU+Bq72hwkBfoM32fotvELzJOfck2b2c7yzsorxzh67FW8idWtc6x/zKbzv5S6gW3inc+5TM5sO/AS4DK9X5xq/fZuOnSzOOecXvz/G+xk5EG8O2cv4c4ycc3/xT/u/HW+pgb/izV2bup9D3wYcine25Q68YlE9TZJVzJtzKSIiIiL7o+E5ERERkQBUNImIiIgEoKJJREREJAAVTSIiIiIBqGgSERERCUBFk4iIiEgAKppEREREAlDRJCIiIhKAiiYRERGRAFQ0iYiIiASgoklEREQkABVNIiIiIgH8f0YCifhjei6dAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure( figsize=(6,4) , dpi=100)\n",
    "\n",
    "plt.plot( bounds_std[0,:], bounds_rnd[0,:],'.', ms=3)\n",
    "plt.plot( bounds_std[1,:], bounds_rnd[1,:],'.', ms=3)\n",
    "plt.plot( bounds_std[0,:], bounds_std[0,:],'k', lw=0.5)\n",
    "plt.plot( bounds_std[1,:], bounds_std[1,:],'k', lw=0.5)\n",
    "plt.grid()\n",
    "plt.xlabel('Bounds in standard medium')\n",
    "plt.ylabel('Bounds in random media')\n",
    "plt.legend(('Upper','Lower'))\n",
    "\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can save the bounds with an explicit description as a .txt file using `Q.save_bounds()`..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maxbounds exported as bounds_test.txt on 2022-02-16 18:27:21.771082.\n"
     ]
    }
   ],
   "source": [
    "Q.save_bounds('bounds_test.txt', description='This is only a toy file.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and then load them using `Q.load_bounds()`.\n",
    "\n",
    "To check that they were loaded correctly, compare the bounds of a reaction before and after loading the .txt file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Upper bound of ACLSm: 1.4712\n",
      "Random value assigned: 0.6593\n",
      "Maxbounds loaded from bounds_10.0.csv on 2022-02-16 18:27:21.782826.\n",
      "Value after loading bounds: 11.2277\n",
      "\n",
      "Some description of the loaded file:\n",
      "nsamples=10000.\n"
     ]
    }
   ],
   "source": [
    "reaction = 'ACLSm'\n",
    "print('Upper bound of %s: %1.4f' % (reaction, Q.upper_bounds[reaction]) )\n",
    "Q.upper_bounds[reaction] = np.random.rand()\n",
    "\n",
    "print('Random value assigned: %1.4f' % Q.upper_bounds[reaction])\n",
    "\n",
    "# Load bounds\n",
    "out = Q.load_bounds('bounds_10.0.csv')\n",
    "\n",
    "print('Value after loading bounds: %1.4f' % Q.upper_bounds[reaction])\n",
    "print('')\n",
    "print('Some description of the loaded file:')\n",
    "print(out.split(';')[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Observe the gene reaction rule translation in action\n",
    "All quantitative interpretation of gene reaction rules are stored as lambda functions in the dictionary `Q.GRR`.\n",
    "\n",
    "Here we show an example for a specific reaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The complete gene reaction rule of CYOR_u6m is:\n",
      "YHR001W_A and YJL166W and YGR183C and YBL045C and YFR033C and YEL024W and Q0105 and YDR529C and YOR065W and YPR191W\n",
      "\n",
      "The list of genes involved is:\n",
      "  Identifier\t Name\tDosage\n",
      "  ----------\t ----\t------\n",
      "     YBL045C\t COR1\t1.000\n",
      "       Q0105\t  COB\t1.000\n",
      "     YFR033C\t QCR6\t1.000\n",
      "     YGR183C\t QCR9\t1.000\n",
      "     YDR529C\t QCR7\t1.000\n",
      "     YJL166W\t QCR8\t1.000\n",
      "   YHR001W_A\tQCR10\t1.000\n",
      "     YOR065W\t CYT1\t1.000\n",
      "     YPR191W\t QCR2\t1.000\n",
      "     YEL024W\t RIP1\t1.000\n",
      "\n",
      "By introducing these gene dosages in the GRR,\n",
      " \bresults in a relative flux of the reaction: 1.000\n"
     ]
    }
   ],
   "source": [
    "reaction = 'CYOR_u6m'\n",
    "\n",
    "print('The complete gene reaction rule of %s is:' % reaction)\n",
    "print(Q.model.reactions.get_by_id(reaction)._gene_reaction_rule)\n",
    "\n",
    "print('\\nThe list of genes involved is:')\n",
    "print('  %10s\\t%5s\\t%s' % ('Identifier', 'Name', 'Dosage'))\n",
    "print('  %10s\\t%5s\\t%s' % ('----------', '----', '------'))\n",
    "_=[print('  %10s\\t%5s\\t%1.3f' % (gene.id, gene.name, Q.gene_dosage[gene.id])) for gene in Q.model.reactions.get_by_id(reaction).genes]\n",
    "\n",
    "print('\\nBy introducing these gene dosages in the GRR,\\n',\n",
    "      '\\bresults in a relative flux of the reaction: %1.3f' % Q.GRR[reaction](Q) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Working with relative gene dosages\n",
    "In light of the previous data, lets try to modify the gene dosage of one of the genes involved in the reaction. As the previous reaction is only composed of \"AND\" statements, the relative flux of the reaction should be the minimum of the dosages of those genes.\n",
    "\n",
    "When calling a GRR lambda function, we need as first argument the object itself, this is because in general, we will only access these functions from within the object itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dosage of YPR191W: 0.562\n",
      "\n",
      "By introducing these gene dosages in the GRR,\n",
      " \bresults in a relative flux of the reaction: 0.562\n"
     ]
    }
   ],
   "source": [
    "Q.set_dosage('YPR191W', np.random.rand() )\n",
    "\n",
    "print('Dosage of YPR191W: %1.3f' % Q.gene_dosage['YPR191W'])\n",
    "print('\\nBy introducing these gene dosages in the GRR,\\n',\n",
    "      '\\bresults in a relative flux of the reaction: %1.3f' % Q.GRR[reaction](Q) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To set up a random genotype, although you can call `Q.set_dosage()` and iterate through all genes in `Q.GENES` you can also use `Q.set_genotype()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYpElEQVR4nO3deZRnZX3n8feHBoIKgkqjrDYoLq0DKo1LNBFjojROgoonggQiUQlRHJ3oCOO+jEaHZDQK2hLCMOBCMgqISStJhqPoAIFGEWgRp9NsTas0grJohIbv/HFv66+L2sC+VU9Vv1/n1OEuT937/T1U/T713Hv7+aWqkCSpNVvMdgGSJI3HgJIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ1yYCSGpbk1Um+ObJ+Z5K9NtGx357klH55UZJKsuUmOvYefa0LNsXxtHkyoNSsJIcm+dckdyW5uV9+fZIMfN6NQqElVbVtVa2erE2SA5KsmcaxPlRVr90UdSW5Lsnvjhz7hr7WezfF8bV5MqDUpCRvAf4aOAF4DPBo4BjgucDWs1javLCpRkrSkAwoNSfJ9sD7gddX1Req6o7qfLuqDq+qX/TtfiPJXya5IcmPkixL8pB+3wFJ1iR5Sz/6+kGSo0bPkeT0JOuSXJ/knUm2SPJkYBnwnP4S1U+S7N8ff8uR7z8kyeX98nuTfCHJ3yW5I8m3kuw70naXJF/sz3Vtkv80yWt/VJJzk9ye5BLgcWP2V5LH98sHJfluf86bkrw1ycOArwC79PXf2Z9/Q42fSXI78Op+22fGlPAnSdb2/fWWkfOeluS/jaz/cpSW5AxgD+DL/fneNvaSYV/DuUluTbIqyetGjvXeJH/f//+4I8nKJEsm/gnR5sKAUoueA/wG8KUp2n0EeALwNODxwK7Au0f2PwbYvt/+GuCkJI/o932i37cX8HzgSOCoqrqabqR2UX+JaoequhT4MfB7I8f+I+CMkfWDgf8NPBL4HHBOkq2SbAF8GfhOX8cLgTcnefEEr+kk4N+BnYE/6b8m8rfAn1bVdsBTgfOr6i5gKbC2r3/bqlo7UuMXgB2Az05wzBcAewMvAo4fvWw3kao6ArgB+P3+fP99nGafB9YAuwCvAD6U5IUj+/8AOLOv7VzgxKnOq/nPgFKLdgRuqar1GzYkubAfzfw8yW/396FeB/znqrq1qu4APgQcOnKce4D3V9U9VbUcuBN4Yn/j/pXAf+1HZ9cBfwUcMUlN/4sulEjySODFdEG0wWX9aO8e4H8A2wDPBvYHFlbV+6vq7v7+0d+MqXPDa1wAHAK8u6ruqqqr+vNO5B5gcZKHV9VtVfWtSdpCF7rnVNV9VfXzCdq8rz/3lcD/BA6b4phTSrI78DzguKr696q6HDiFjfv7m1W1vL9ndQaw7/2PpM2N16HVoh8DOybZckNIVdVvAvSXlbYAFgIPBS4beWYiwOhTYz8eDTngZ8C2dAG4NXD9yL7r6UY4E/kMcHWSbYE/BL5RVT8Y2X/jhoWquq+vcxeg6C63/WSk7QLgG+OcYyHd7+SNI9uuH6fdBocA7wQ+nOQK4PiqumiS9jdOsm+8NtcD/2Ea3zOVXYANf0SMHnv0Mt4PR5Z/Bmwz+v9fmydHUGrRRcAv6C5JTeQW4OfAU/rLcDtU1fZVte00jn8L3ejjsSPb9gBu6pfvN8V/Vd3U1/Uyur/8zxjTZPcNC/1lvd2AtXRv+NeO1LhDVW1XVQeNU9c6YP3osfq6xlVVl1bVwcBOwDnA309U/xTbx30d/bk3XB68i+4Pgg0e8wCOvRZ4ZJLtxhz7pgnaS4ABpQZV1U+A9wGfTPKKJNv2DzA8DXhY3+Y+uktlH02yE0CSXSe5tzN6/Hvp3sw/mGS7JI8F/pxulATwI2C3JGOfFjwdeBvdqOLsMfv2S/Ly/qGAN9MF7MXAJcDtSY5L8pAkC5I8Ncn+E9R1FvDeJA9Nshj44/FeQ5KtkxyeZPv+suLtwIZHun8EPCrdwyYP1Lv6cz8FOAr4u3775cBBSR6Z5DH9axz1I7r7efdTVTcCFwJ/kWSbJPvQ3ROc6D6YBBhQalR/o/3P6QLhZro3wE8Dx9G92dEvrwIu7p9M+xfgidM8xRvpRgWrgW/S3U86td93PrAS+GGSW0a+52y6UdfZ/cMIo75Ed1/rNroR1sv7e1/3Ar9P9yDHtXSjt1PoHtAYz7F0lyF/CJxGdx9oIkcA1/Wv/Rj6e2RV9T26hxJW9/ftdpnkGGN9na5P/w/wl1X1T/32M+ge9LgO+Cd+FVwb/AXwzv58bx3nuIcBi+hGU2cD76mqf34AdWkzFD+wUJq+JP9G9+Tcv4xsey/w+Kr6o1krTJqHHEFJ05TkELp7LefPdi3S5sCn+KRpSPI1YDFwRH//S9LAvMQnSWqSl/gkSU2ac5f4dtxxx1q0aNFslyFJ2kQuu+yyW6pq4djtcy6gFi1axIoVK2a7DEnSJpJk3BlTvMQnSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlq0pybSUKSJrLo+H+ctXNf9+GXzNq55ytHUJKkJhlQkqQmeYlP895sXfbxko/063EEJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkv8OSpI2AadZ2vQcQUmSmuQISpLmuPk6enMEJUlqkgElSWqSl/ikgczXyy7STHEEJUlqkgElSWqSASVJapIBJUlq0mABleTUJDcnuWqC/Uny8SSrklyR5BlD1SJJmnuGfIrvNOBE4PQJ9i8F9u6/ngV8qv/v4Hy6ShrObP5+aX4ZbARVVRcAt07S5GDg9OpcDOyQZOeh6pEkzS2zeQ9qV+DGkfU1/bb7SXJ0khVJVqxbt25GipMkza7ZDKiMs63Ga1hVJ1fVkqpasnDhwoHLkiS1YDYDag2w+8j6bsDaWapFktSY2Qyoc4Ej+6f5ng38tKp+MIv1SJIaMthTfEk+DxwA7JhkDfAeYCuAqloGLAcOAlYBPwOOGqoWSdLcM1hAVdVhU+wv4A1DnV+SNLc5k4QkqUkGlCSpSX4e1GbEGTQ2H87moPnAEZQkqUkGlCSpSQaUJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUlOdTTDnIJGkqbHEZQkqUkGlCSpSQaUJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUkGlCSpSQaUJKlJBpQkqUkGlCSpSYMGVJIDk1yTZFWS48fZv32SLyf5TpKVSY4ash5J0twxWEAlWQCcBCwFFgOHJVk8ptkbgO9W1b7AAcBfJdl6qJokSXPHkCOoZwKrqmp1Vd0NnAkcPKZNAdslCbAtcCuwfsCaJElzxJABtStw48j6mn7bqBOBJwNrgSuBN1XVfWMPlOToJCuSrFi3bt1Q9UqSGjJkQGWcbTVm/cXA5cAuwNOAE5M8/H7fVHVyVS2pqiULFy7c1HVKkho0ZECtAXYfWd+NbqQ06ijgrOqsAq4FnjRgTZKkOWLIgLoU2DvJnv2DD4cC545pcwPwQoAkjwaeCKwesCZJ0hyx5VAHrqr1SY4FzgMWAKdW1cokx/T7lwEfAE5LciXdJcHjquqWoWqSJM0dgwUUQFUtB5aP2bZsZHkt8KIha5AkzU3OJCFJatKgIyhpg0XH/+NslyBpjnEEJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWrSpAGVZMuZKkSSpFFTjaAu2bCQ5BMD1yJJ0i9NFVAZWX7ukIVIkjRqqoCqGalCkqQxprrH9KQkV9CNpB7XL9OvV1XtM2h1kqTN1lQB9eQZqUKSpDEmvcRXVddX1fXAT4Gd+q+fjGyfVJIDk1yTZFWS4ydoc0CSy5OsTPL1B/MiJEnzz6QjqCRbAycDLwWupbu099gkZwPHVNXdk3zvAuAk4PeANcClSc6tqu+OtNkB+CRwYFXdkGSnX+/lSJLmi6kekngnsBWwe1U9vaqeBuxBF2zvmuJ7nwmsqqrVfZCdCRw8ps2rgLOq6gaAqrr5AdYvSZqnpgqolwOvq6o7Nmzol18PvGyK790VuHFkfU2/bdQTgEck+VqSy5IcOd6BkhydZEWSFevWrZvitJKk+WCqgLqvqn42dmNV3cnUj6BnnG1jv2dLYD/gJcCLgXclecI45zu5qpZU1ZKFCxdOcVpJ0nww1VN8leQRjB82903xvWuA3UfWdwPWjtPmlqq6C7gryQXAvsD3pzi2JGmemyqgtgcuY3qjobEuBfZOsidwE3Ao3T2nUV8CTuzn/NsaeBbw0amKliTNf5MGVFUterAHrqr1SY4FzgMWAKdW1cokx/T7l1XV1Um+ClxBNyI7paquerDnlCTNH1M9Zr4T8Hbg8XQh8uGqun26B6+q5cDyMduWjVk/AThhuseUJG0epnpI4nTgLuATwHbAxwevSJIkpr4H9Ziqeke/fF6Sbw1dkCRJMHVAZcxTfAtG16vq1iGLkyRtvh7IU3zV/3fDKKqAvYYrTZK0OZvWU3xJtgAOB/asqvcn2QPYefjyJEmbq6kektjgJODZwGH9+h3AiYNUJEkSU1/i2+BZVfWMJN8GqKrb+pnOJUkaxHRHUPf0H59RAEkWMvVUR5IkPWjTDaiPA2cDOyX5IPBN4EODVSVJ2uxN6xJfVX02yWXAC+me5HtpVV09aGWSpM3adO9BUVXfA743YC2SJP3SdC/xSZI0owwoSVKTDChJUpMMKElSkwwoSVKTDChJUpMMKElSkwwoSVKTDChJUpMMKElSkwwoSVKTDChJUpMMKElSkwwoSVKTDChJUpMMKElSkwwoSVKTBg2oJAcmuSbJqiTHT9Ju/yT3JnnFkPVIkuaOwQIqyQLgJGApsBg4LMniCdp9BDhvqFokSXPPkCOoZwKrqmp1Vd0NnAkcPE67NwJfBG4esBZJ0hwzZEDtCtw4sr6m3/ZLSXYFXgYsm+xASY5OsiLJinXr1m3yQiVJ7RkyoDLOthqz/jHguKq6d7IDVdXJVbWkqpYsXLhwU9UnSWrYlgMeew2w+8j6bsDaMW2WAGcmAdgROCjJ+qo6Z8C6JElzwJABdSmwd5I9gZuAQ4FXjTaoqj03LCc5DfgHw0mSBAMGVFWtT3Is3dN5C4BTq2plkmP6/ZPed5Ikbd6GHEFRVcuB5WO2jRtMVfXqIWuRJM0tziQhSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWrSoAGV5MAk1yRZleT4cfYfnuSK/uvCJPsOWY8kae4YLKCSLABOApYCi4HDkiwe0+xa4PlVtQ/wAeDkoeqRJM0tQ46gngmsqqrVVXU3cCZw8GiDqrqwqm7rVy8GdhuwHknSHDJkQO0K3DiyvqbfNpHXAF8Zb0eSo5OsSLJi3bp1m7BESVKrhgyojLOtxm2YvIAuoI4bb39VnVxVS6pqycKFCzdhiZKkVm054LHXALuPrO8GrB3bKMk+wCnA0qr68YD1SJLmkCFHUJcCeyfZM8nWwKHAuaMNkuwBnAUcUVXfH7AWSdIcM9gIqqrWJzkWOA9YAJxaVSuTHNPvXwa8G3gU8MkkAOuraslQNUmS5o4hL/FRVcuB5WO2LRtZfi3w2iFrkCTNTc4kIUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlqkgElSWqSASVJapIBJUlq0qABleTAJNckWZXk+HH2J8nH+/1XJHnGkPVIkuaOwQIqyQLgJGApsBg4LMniMc2WAnv3X0cDnxqqHknS3DLkCOqZwKqqWl1VdwNnAgePaXMwcHp1LgZ2SLLzgDVJkuaILQc89q7AjSPra4BnTaPNrsAPRhslOZpuhAVwZ5Jrfs3adgRu+TWPMZ/YHxuzPzZmf2zM/hiRj2yS/njseBuHDKiMs60eRBuq6mTg5E1RFECSFVW1ZFMdb66zPzZmf2zM/tiY/bGxIftjyEt8a4DdR9Z3A9Y+iDaSpM3QkAF1KbB3kj2TbA0cCpw7ps25wJH903zPBn5aVT8YeyBJ0uZnsEt8VbU+ybHAecAC4NSqWpnkmH7/MmA5cBCwCvgZcNRQ9YyxyS4XzhP2x8bsj43ZHxuzPzY2WH+k6n63fCRJmnXOJCFJapIBJUlq0rwOKKda2tg0+uPwvh+uSHJhkn1no86ZMlV/jLTbP8m9SV4xk/XNtOn0R5IDklyeZGWSr890jTNpGr8v2yf5cpLv9P0xU/fQZ1ySU5PcnOSqCfYP815aVfPyi+7BjH8D9gK2Br4DLB7T5iDgK3T/HuvZwL/Odt2z3B+/CTyiX166uffHSLvz6R7oecVs1z3LPx87AN8F9ujXd5rtume5P94OfKRfXgjcCmw927UP1B+/DTwDuGqC/YO8l87nEZRTLW1syv6oqgur6rZ+9WK6f5c2X03n5wPgjcAXgZtnsrhZMJ3+eBVwVlXdAFBV87lPptMfBWyXJMC2dAG1fmbLnBlVdQHd65vIIO+l8zmgJppG6YG2mS8e6Gt9Dd1fRPPVlP2RZFfgZcCyGaxrtkzn5+MJwCOSfC3JZUmOnLHqZt50+uNE4Ml0kwtcCbypqu6bmfKaM8h76ZBTHc22TTbV0jwx7dea5AV0AfW8QSuaXdPpj48Bx1XVvd0fyfPadPpjS2A/4IXAQ4CLklxcVd8furhZMJ3+eDFwOfA7wOOAf07yjaq6feDaWjTIe+l8DiinWtrYtF5rkn2AU4ClVfXjGaptNkynP5YAZ/bhtCNwUJL1VXXOjFQ4s6b7+3JLVd0F3JXkAmBfYD4G1HT64yjgw9XdhFmV5FrgScAlM1NiUwZ5L53Pl/icamljU/ZHkj2As4Aj5ulfxaOm7I+q2rOqFlXVIuALwOvnaTjB9H5fvgT8VpItkzyU7tMJrp7hOmfKdPrjBrrRJEkeDTwRWD2jVbZjkPfSeTuCqranWppx0+yPdwOPAj7ZjxrW1zydtXma/bHZmE5/VNXVSb4KXAHcB5xSVeM+djzXTfPn4wPAaUmupLvEdVxVzcuP4UjyeeAAYMcka4D3AFvBsO+lTnUkSWrSfL7EJ0mawwwoSVKTDChJUpMMKElSkwwoSVKTDChpAkkeneRzSVb3U/tclORlM1zDnTN5PqklBpQ0jn4C0HOAC6pqr6raj+4fa87nCXSlphhQ0vh+B7h79B/sVtX1VfUJgCQLkpyQ5NL+82/+tN9+QD+Z6heSfC/JZ/uwI8l+Sb7ej8bOG2+2537mgov6435gZHv6812V5Mokr+y375zkgv4zmq5K8lv99k8lWdF/TtH7Ro5zUF/XN/vP7/mHfvvD+s/8uTTJt5OMN7O7NKMMKGl8TwG+Ncn+19BN57I/sD/wuiR79vueDrwZWEz3eULPTbIV8Am6z5TaDzgV+OA4x/1r4FP9cX84sv3lwNPo5r77XeCEPuBeBZxXVRv2Xd63f0c/C8g+wPOT7JNkG+DTdPMsPo/uM4w2eAdwfn/eF/THf9gkr18a3Lyd6kjalJKcRDe7+939m/iLgH3yq0/Z3R7YG7gbuKSq1vTfdzmwCPgJ8FS6Ga+hmz5nvLnKngsc0i+fAXykX34e8Pmquhf4UbpPs92fbs64U/sAPKeqLu/b/2GSo+l+x3emC8stgNVVdW3f5vPA0f3yi4A/SPLWfn0bYA/m71x7mgMMKGl8K/lVUFBVb0iyI7Ci3xTgjVV13ug3JTkA+MXIpnvpfs8CrKyq50zj3OPNPzbu531U1QVJfht4CXBGkhOAbwBvBfavqtuSnEYXOJN9ZkiAQ6rqmmnUJ80IL/FJ4zsf2CbJn41se+jI8nnAn/UjF5I8YYpLYtcAC5M8p2+/VZKnjNPu/9I9jAFw+Mj2C4BX9ve+FtJ9BPclSR4L3FxVfwP8Ld3Hcj8cuAv4aT/L9tL+GN8D9kqyqF9/5ZjX88aR+2VPn+S1SDPCEZQ0jqqqJC8FPprkbcA6ujf94/omp9BduvtW/6a+DnjpJMe7u78c+PEk29P97n2MbqQ26k3A55K8ie6j5jc4G3gO8B26EdbbquqHSf4Y+C9J7gHuBI6sqmuTfLs/9mq60KOqfp7k9cBXk9zCxp9b9IG+niv613Md8B+n7ilpOM5mLm1GkmxbVXf2IXQS8P+q6qOzXZc0Hi/xSZuX1/UPbqyke7Dj07NbjjQxR1CSpCY5gpIkNcmAkiQ1yYCSJDXJgJIkNcmAkiQ16f8DmijNwXn3+gAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Slow method:\n",
    "#_ = [Q.set_dosage(gene, np.random.rand()) for gene in Q.GENES]\n",
    "\n",
    "# Fast method:\n",
    "Q.set_genotype( np.random.rand(Q.N_GENES) )\n",
    "\n",
    "# Recover the genotype\n",
    "g = [ value for key, value in Q.gene_dosage.items() ]\n",
    "\n",
    "# Plot a histogram of gene dosages in the current individual\n",
    "plt.figure( figsize=(6,4) )\n",
    "plt.hist(g, bins=11, density=True)\n",
    "plt.xlabel('Gene dosage')\n",
    "plt.ylabel('ePDF')\n",
    "plt.title('Genotype distribution')\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**WARNING:** Gene dosages can be modified any time (in and out of context managers), however their effect will not take place until a call to `Q.apply_dosage()`, which changes the model reaction bounds. We strongly advise to use **ALWAYS** this within a context manager."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Biomass production before calling apply_dosage() : 0.5139\n",
      "Biomass production after calling apply_dosage()  : 0.0257\n",
      "Biomass production outside the context manager   : 0.5139\n"
     ]
    }
   ],
   "source": [
    "print('Biomass production before calling apply_dosage() : %1.4f' % Q.slim_optimize())\n",
    "\n",
    "with Q.model:\n",
    "    Q.apply_dosage()\n",
    "    print('Biomass production after calling apply_dosage()  : %1.4f' % Q.slim_optimize())\n",
    "\n",
    "print('Biomass production outside the context manager   : %1.4f' % Q.slim_optimize())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To restore the genotype, you can call `Q.reset_dosage()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Biomass production before calling apply_dosage() : 0.5139\n",
      "Biomass production after calling reset_dosage()  : 0.5139\n",
      "Biomass production outside the context manager   : 0.5139\n"
     ]
    }
   ],
   "source": [
    "print('Biomass production before calling apply_dosage() : %1.4f' % Q.slim_optimize())\n",
    "\n",
    "\n",
    "Q.reset_dosage()\n",
    "with Q.model:\n",
    "    Q.apply_dosage()\n",
    "    print('Biomass production after calling reset_dosage()  : %1.4f' % Q.slim_optimize())\n",
    "\n",
    "print('Biomass production outside the context manager   : %1.4f' % Q.slim_optimize())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### You can use the functions in the `QuantitativeMutation` object to generate random genotypes, and set their clipping values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAEYCAYAAABBWFftAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtRklEQVR4nO3dffzm9Zz3/8fTTLEqihlpuzBxxro4K4xCLsLJr3IRv9NvFdI62cGWi10Wyy6W0+/HaQnL1jkyG5acTkSI6mSTq2gi1UhpRzSKmUgXspvJ6/fH5zM65uj4zveYmc/3uPh+H/fb7bjN5/N+vz/H8Tq+Ol4+r/fnKlWFJEmSJGn73WHcAUiSJEnSfGGBJUmSJEkdscCSJEmSpI5YYEmSJElSRyywJEmSJKkjFliSJEmS1BELLHUuyTlJXtizvjTJZUnuNOa4npbk4+OMQdLoJbkpyb2HGLcsSSVZ3NP2oiTvntMAB8fysiRvG/XnSvNJkpOS/F3P+kuS/KLNCXcfweffMckPktxzrj9ra7g/NPcssKZUkivbJLFTT9sLk5wzxrBm8lrgn6vq3wGS/GmSbya5eZh4kzw7yU+S/CbJZ5LcbYZx90hyapKrk1yf5BtJDt7UX1WnAw9Ksn9H30uaF9p88tt2p2PT633jjmtb9E/wAFTVzlW1dhvea0fgb4F39LQdmOSCNn9dkOTALWx/xySrktyQ5OdJ/qqvv9q8tulvfnJP90rguUnusbVxS/NF+xv5T31tb0ryL8NsX1Uvrqq3tNvtALwLeFKbE37ZfcS3swI4t6p+3saQJG9P8sv29T+SZNCGPRM+vXn57waNnWH7tyS5OMnGJG/q7XN/aO5ZYE23xcDLt+cN2h/7nP13kOSOwLFAbzL8FfBuYNbZ2SQPBP4ncAywO3Az8E8zDN8ZOB94KHA34EPAF5Ls3DPmVJqEJ2lzT213Oja9jh93QBPgSOCHVfUz+EPB9VmafLYbTY75bNs+yJuA/YB7AY8DXp3ksL4xB/T8zf9QGLYTUl8Entfh95EWst2BOwFrtnbD7dhXehHwkZ71FcDTgQOA/YGntGO2ZNeeHPGWrfjsK4BXA1+Yod/9oTlkgTXd3gG8Ksmu/R1JHpnk/PZIzvlJHtnTd06Styb5Bk3Bcu92luQvkvwoyY3tzMd9knyrnX39xKadiCS7Jfl8kg1JrmuX95ohxoOBX1fVuk0NVfV/quoTwNVDfMfnAJ+rqnOr6ibg74D/O8ku/QOram1VvauqrqmqW6tqJbAjcL+eYecATx7ic6UFL8mJST7Zs/72JF9udzYOTbIuyeuSXNseBXtOz9i7Jvlwmyd+kuRvN+2gJPmzJF9P8g9tDvlxksP7tv1gkmuS/CzJf0+yaLZtk7wVeDTwvvQcheudBU/y5CTfa/PaVf0zu30OB77as34ozcTWu6vqP6rqvUCAx8+w/fOAt1TVdVV1KfAB4M9m/8v/wTmYr6QZ9eShVyZZ3+aM5/f0n9Lmj/sCl7XNv07ylbZ/LveV9gHuA3y7J+RjgXdW1bp24uadbF1OGFpVfaiqvgjcOMOQczC/zBkLrOm2muYH8qrexjSn0H0BeC9wd5pD4l/I5ucbH0Mzc7EL8JO27TCaoz8Pp5n1WElT4OwNPAg4uh13B+CfaWZl9wF+C8x0OtF/5rakti0eCHx/00pV/RtwC3Df2TZMc+rOjjSzOJtcCixLcpftiElaKF4J7N8WNY8GXgAcW1XV9t8TWALsSbPjsDLJpgmNfwTuCtwbeCxNsfH8nvc+mCY3LAH+B/DB5A+nynwI2Aj8J+DBwJOAF862bVW9HvgacPwWjsL9po1lV5qdi5ckefoM378/fz0QuKjn+wNc1LZvJsluwB/Tk7/a5f6x56Y5ffDTSZb19V1KM9MtaWb3pMk1e9LkqPe3v78/qKrLue23t2tVPX4E+0r/GVhbVRt73m+zfRoG54R+P2mLyH9OsmSWsVvD/aE5ZIE1/d4AvDTJ0p62JwM/qqqPVNXGqjoV+CHw1J4xp1TVmrb/d23b26vqhqpaA1wCnNUeFbqe5lSVBwNU1S+r6lNVdXNV3Qi8lWYHapBdmXn2ZBg7A9f3tV1Pk+xm1CaMjwB/38a/yaZYdt2OmKT56DNJft3z+vOquhl4Ls2Ox78AL+09Gt36u/Zozldpdlb+tD3a9Czgb6rqxqq6kmam9pie7X5SVR+oqltpCqo9gN2T7E5z5OgVVfWbqloPnAAcNdu2w3zJqjqnqi6uqt9X1UU0p8kMm7+2Jh/t3NM/09jHAsuAP6E5ov/59Nxgo/3su874ZSQB/A54c1X9rqrOAG5i8zNXZjKn+0oM3v/pzyHXAzv3TC71uhZ4GM1k9kNpcsdHh/hew3J/aA4tnn2IJllVXZLk8zQ3kri0bf5jbptp2eQnNLM7m1w14O1+0bP82wHr9wRIcmeaHZ7DaK5DANglyaJ2h6fXdcxSDM3iJqB/duUubKFoS/JHwOeA86rq/+vr3hTLr7cjJmk+enpV/Z/+xqr6TpK1wD2AT/R1X1dVv+lZ/wlN/llCc/T4J319vTno5z2fcXO7f7EzzfWTOwDX9Oxz3IHNc9ZM284qzY1v3kYz07wjcEfgf88wvD9/bU0+uqmn/98Hja2qc9vFW5K8HLgBuD9wcdu+C7cv6KSF5FaafNBrB5qiapNf9h0lupnh8sGc7isxeP+nP4fcBbip76g4AO1lEas3fWaS42ny4l2q6oYBcW0t94fmkEew5oc3An/ObUnhapoZj177AD/rWb/dj3krvJJmdujgqroL8Ji2fdAMzEUMcTrfFqyh5xSZNLdaviNw+aDBaW6q8Rma7zrowtH7A1d2lJykeS/JcTS/uatpTofptVt67mRKk2euppl5/R2b56H+HDSTq4D/AJZU1a7t6y5VNdtpNJvMlts+BpwO7F1VdwVOYnDugtvnrzU0p0z2jt+fARfNV9V1wDVsforfAYPG9sXe+973Z/PTiaSF5qc0R3l77cvtC6NtMdf7ShfRXLfVezBjs30aZs8JvTbFMlO+2lruD80hC6x5oKquAP4X8LK26Qzgvmlub744ybOABwCf7+gjd6GZpfl1ew7zG7cw9jvArkn+MCOUZFGaZ2ItBu6Q5E5pbp+6qf/KJH/Wrn4UeGqSR7c7cm8GPt2emriZ9j0+2cb2vKr6/YB4HktzCF/SLNoLw/87zWmCx9DcBe/AvmF/n2TH9hqtpwD/uz2S/QngrUl2SXIv4K/Y/G6iA1XVNcBZwDuT3CXJHdqLyGc6ja/fL2iu+5rJLsCvqurfkxwEPHsLY89g89MHz6GZUX9Zmluwb7rG6yszbP9h4G/T3BjoT2gmwk6B5g6paW75vijNnU7fSbNjd2nP9uYrLXT/i+Y3tFebC/4LzSl8n5xlu2HM6b5Sezr1j4CDepo/DPxVkj2T/DHNhPUpg7ZPcnCS+7Xf++4014qds+myh/ba2Ctn+vwkO7T7WncAFrf7Wot6hphf5pAF1vzxZmAnaK6RotnReSXwS5pZ56dU1bUdfda7gT+imaU+D/jSTAOr6haa5PHcnuZjaIqgE2nu+PVbmrtrbboN8t3b96U9x/nFNIXWepqdo7/Y9EZpHiJ4Urv6SJrv/SSa4m/TcyMe3fPZR9Pc9l3S5j6XzZ+3chpNQfT2qvp+Vf0IeB3wkfZIMTSn6l1HMxP8UeDFVfXDtu+lNDeUWAt8nebI0aohY3kezel7P2jf/5M011kN4z3AM9PcYfC9A/r/AnhzkhtprmHtP+2x1+eAP2l3hDbls6e38f0a+G80p1beApDkOUl6Z6PfCPwbzWz7V4F3VNWmfLk7zc7jDTR/o2U0efp37XvdCTiC5hozaaF6M/BNmhxyHc1NbZ5TVZds7xuPYF8JbnvMTO/652hOA76E5rrVP+yTJFmT2+7Gem+a/asb27H/wW030IDmphrf2MJnf4Bm/+po4PXtcm8s7g/NoQw47VPqVHsDjq8BD66q384y9lHAcVV19JbGbWMcTwWOqao/7fq9pYUmyaHAv1TVTI9omBeSrAAeUFWvGPHnvpTmNMb+0zIlTYl2Mup7wBPao/NdvvdZwMvbR0Bs7bbuD80xCyxJ0lZbKAWWJElby1MEJUmSJKkjHsGSJEmSpI54BEuSJEmSOjLSBw0vWbKkli1bNsqPlDQCF1xwwbVVtXTccWyJ+UeanyY9/5h7pPlrpvwz0gJr2bJlrF69evaBkqZKki4e+jinzD/S/DTp+cfcI81fM+UfTxGUJEmSpI5YYEmSJElSRyywJEmSJKkjFliSJEmS1BELLEmSJEnqiAWWJEmSJHXEAkuSJEmSOmKBJUmSJEkdscCSJEmSpI4sHncAkibDCWdfvtn6Xz7xvmOKRNJC0p97wPwjzVcLZV/DI1iSJEmS1BGPYEmSJEmaCPPhKJcFliRJkqSRG3SK8HxggSVpYiVZBTwFWF9VD+ppfylwPLAR+EJVvXpMIUqSNDXm8prH+XDkqStegyVpkp0CHNbbkORxwJHA/lX1QOAfxhCXpCmVZO8k/5rk0iRrkrx8wJgkeW+SK5JclOQhPX2HJbms7XvtaKOXNA08giVpYlXVuUmW9TW/BHhbVf1HO2b9yAOTNM02Aq+squ8m2QW4IMnZVfWDnjGHA/u1r4OBE4GDkywC3g88EVgHnJ/k9L5tpS3ySM/WGeao26TdjdQCS9K0uS/w6CRvBf4deFVVnT9oYJIVwAqAffbZZ3QRSppYVXUNcE27fGOSS4E9gd4i6Ujgw1VVwHlJdk2yB7AMuKKq1gIk+Xg71gJLE2m+XuM06SywJE2bxcBuwMOBhwGfSHLvdkdoM1W1ElgJsHz58tv1S1rY2iPkDwa+3de1J3BVz/q6tm1Q+8ED3tfJHQ1t3EdfPKLWPQssSdNmHfDptqD6TpLfA0uADeMNS9I0SbIz8CngFVV1Q3/3gE1qC+2bNzi5s2DNVbEy7iJsW8zl0bNJPzJngSVp2nwGeDxwTpL7AjsC1441IklTJckONMXVR6vq0wOGrAP27lnfC7iaJt8MapfGrquiYxqLuUHGeWTOAkvSxEpyKnAosCTJOuCNwCpgVZJLgFuAYwedHihJgyQJ8EHg0qp61wzDTgeOb6+xOhi4vqquSbIB2C/JvsDPgKOAZ48ibqnfpB/FWcgssCRNrKo6eoau5440EEnzySHAMcDFSS5s214H7ANQVScBZwBHAFcANwPPb/s2JjkeOBNYBKyqqjUjjV5TxSJoss3VUa5ZC6wkewMfBu4J/B5YWVXvSfIm4M+57bqH11XVGZ1EJUmSNAeq6usMvpaqd0wBx83QdwZNASbNmUkvzCY9vkFGGfMwR7AGPi+i7TuhqnzIpyRJkiQxRIG1hedFSJIkSZJ6bNU1WH3PiziE5gLQ5wGraY5yXTdgG58FIUmSJE2paTwlcJzuMOzAAc+LOBG4D3AgzRGudw7arqpWVtXyqlq+dOnS7Y9YkiRJkibUUAXWoOdFVNUvqurWqvo98AHgoLkLU5IkSZIm36wF1kzPi0iyR8+wZwCXdB+eJEmSJE2PYa7Bmul5EUcnORAo4ErgRXMQnyRJkiRNjWHuIjjT8yJ8BoQkSZIk9Rj6JheSJEmSpC2zwJIkSZKkjlhgSZIkSVJHLLAkSZIkqSMWWJIkSZLUEQssSZIkSeqIBZYkSZIkdcQCS9LESrIqyfoklwzoe1WSSrJkHLFJkiQNYoElaZKdAhzW35hkb+CJwE9HHZCk6baliZu2/6+TXNi+Lklya5K7tX1XJrm47Vs92sglTQsLLEkTq6rOBX41oOsE4NVAjTYiSfPAKQyYuNmkqt5RVQdW1YHA3wBfrarePPS4tn/53IYpaVpZYEmaKkmeBvysqr4/xNgVSVYnWb1hw4YRRCdp0m1h4maQo4FT5zAcSfOQBZakqZHkzsDrgTcMM76qVlbV8qpavnTp0rkNTtK80uabw4BP9TQXcFaSC5Ks2MK2Tu5IC5gFlqRpch9gX+D7Sa4E9gK+m+SeY41K0nz0VOAbfacHHlJVDwEOB45L8phBGzq5Iy1si8cdgCQNq6ouBu6xab0tspZX1bVjC0rSfHUUfacHVtXV7b/rk5wGHAScO4bYJE0wj2BJmlhJTgW+BdwvybokLxh3TJLmvyR3BR4LfLanbacku2xaBp4EDLwToaSFzSNYkiZWVR09S/+yEYUiaZ5oJ24OBZYkWQe8EdgBoKpOaoc9Azirqn7Ts+nuwGlJoNl/+lhVfWlUcUuaHhZYkiRpwZht4qYdcwrN7dx729YCB8xNVJLmE08RlCRJkqSOWGBJkiRJUkcssCRJkiSpIxZYkiRJktQRCyxJkiRJ6ogFliRJkiR1xAJLkiRJkjpigSVJkiRJHbHAkiRJkqSOWGBJkiRJUkcssCRJkiSpIxZYkiRJktSRWQusJHsn+dcklyZZk+Tlbfvdkpyd5Eftv7vNfbiSJEmSNLmGOYK1EXhlVd0feDhwXJIHAK8FvlxV+wFfbtclSZIkacGatcCqqmuq6rvt8o3ApcCewJHAh9phHwKePkcxSpIkSdJU2KprsJIsAx4MfBvYvaqugaYIA+4xwzYrkqxOsnrDhg3bGa6khSTJqiTrk1zS0/aOJD9MclGS05LsOsYQJUmSNjN0gZVkZ+BTwCuq6oZht6uqlVW1vKqWL126dFtilLRwnQIc1td2NvCgqtofuBz4m1EHJWl6DZq46es/NMn1SS5sX2/o6TssyWVJrkjipRGSBhqqwEqyA01x9dGq+nTb/Iske7T9ewDr5yZESQtVVZ0L/Kqv7ayq2tiungfsNfLAJE2zU7j9xE2/r1XVge3rzQBJFgHvBw4HHgAc3V6TLkmbGeYuggE+CFxaVe/q6TodOLZdPhb4bPfhSdIW/TfgizN1eoqypH6DJm6GdBBwRVWtrapbgI/TXI8uSZsZ5gjWIcAxwON7DpcfAbwNeGKSHwFPbNclaSSSvJ7mLqcfnWmMpyhL2kaPSPL9JF9M8sC2bU/gqp4x69q223FyR1rYFs82oKq+DmSG7id0G44kzS7JscBTgCdUVY07HknzyneBe1XVTe2E8meA/Ri8LzQw/1TVSmAlwPLly81R0gKzVXcRlKRxS3IY8BrgaVV187jjkTS/VNUNVXVTu3wGsEOSJTRHrPbuGboXcPUYQpQ04SywJE2sJKcC3wLul2RdkhcA7wN2Ac5uT1k+aaxBSppXktyzvf6cJAfR7Cv9Ejgf2C/Jvkl2BI6iuR5dkjYz6ymCkjQuVXX0gOYPjjwQSfNGO3FzKLAkyTrgjcAOAFV1EvBM4CVJNgK/BY5qT0XemOR44ExgEbCqqtaM4StImnAWWJIkacGYYeKmt/99NEfKB/WdAZwxF3FJmj88RVCSJEmSOmKBJUmSJEkdscCSJEmSpI5YYEmSJElSRyywJEmSJKkjFliSJEmS1BFv0y5JkkbmhLMvH3cIkjSnPIIlSZIkSR2xwJIkSZKkjlhgSZIkSVJHLLAkSZIkqSMWWJIkSZLUEQssSZIkSeqIBZYkSZIkdcQCS5IkSZI6YoElSZIkSR2xwJI0sZKsSrI+ySU9bXdLcnaSH7X/7jbOGCVNl0F5pa//OUkual/fTHJAT9+VSS5OcmGS1aOLWtI0scCSNMlOAQ7ra3st8OWq2g/4crsuScM6hdvnlV4/Bh5bVfsDbwFW9vU/rqoOrKrlcxSfpClngSVpYlXVucCv+pqPBD7ULn8IePooY5I03WbIK73936yq69rV84C9RhKYpHnDAkvStNm9qq4BaP+9x0wDk6xIsjrJ6g0bNowsQEnzxguAL/asF3BWkguSrBhTTJIm3OJxByBJc6WqVtKe3rN8+fIacziSpkiSx9EUWI/qaT6kqq5Ocg/g7CQ/bI+I9W+7AlgBsM8++4wkXkmTwyNYkqbNL5LsAdD+u37M8UiaZ5LsD5wMHFlVv9zUXlVXt/+uB04DDhq0fVWtrKrlVbV86dKlowhZ0gSxwJI0bU4Hjm2XjwU+O8ZYJM0zSfYBPg0cU1WX97TvlGSXTcvAk4CBdyKUtLB5iqCkiZXkVOBQYEmSdcAbgbcBn0jyAuCnwP8zvgglTZsZ8soOAFV1EvAG4O7APyUB2NjeMXB34LS2bTHwsar60si/gKSJZ4ElaWJV1dEzdD1hpIFImje2kFc29b8QeOGA9rXAAbffQpI25ymCkiRJktSRWQusQU88T/KmJD9rn2R+YZIj5jZMSZIkSZp8wxzBOoXBTzw/oX2S+YFVdUa3YUmSJEnS9Jm1wJrtieeSJEmSpMb2XIN1fJKL2lMId5tpUJIVSVYnWb1hw4bt+DhJkiRJmmzbWmCdCNwHOBC4BnjnTAN92J4kSZKkhWKbCqyq+kVV3VpVvwc+wAxPMpckSZKkhWSbCqwke/SsPgOfZC5JkiRJsz9oeIYnnh+a5ECggCuBF81diJIkSZI0HWYtsGZ44vkH5yAWSZIkSZpq23MXQUmSJElSDwssSZIkSeqIBZYkSZIkdcQCS5IkSZI6YoElSZIkSR2xwJIkSZKkjlhgSZpKSf4yyZoklyQ5Ncmdxh2TpMmXZFWS9UkumaE/Sd6b5IokFyV5SE/fYUkua/teO7qoJU0TCyxJUyfJnsDLgOVV9SBgEXDUeKOSNCVOAQ7bQv/hwH7tawVwIkCSRcD72/4HAEcnecCcRippKllgSZpWi4E/SrIYuDNw9ZjjkTQFqupc4FdbGHIk8OFqnAfsmmQP4CDgiqpaW1W3AB9vx0rSZiywJE2dqvoZ8A/AT4FrgOur6qz+cUlWJFmdZPWGDRtGHaak6bQncFXP+rq2baZ2SdqMBZakqZNkN5qZ432BPwZ2SvLc/nFVtbKqllfV8qVLl446TEnTKQPaagvtt38DJ3ekBc0CS9I0+i/Aj6tqQ1X9Dvg08MgxxyRpflgH7N2zvhfNKcgztd+OkzvSwmaBJWka/RR4eJI7JwnwBODSMcckaX44HXheezfBh9OcgnwNcD6wX5J9k+xIc2Od08cZqKTJtHjcAUjS1qqqbyf5JPBdYCPwPWDleKOSNA2SnAocCixJsg54I7ADQFWdBJwBHAFcAdwMPL/t25jkeOBMmjuXrqqqNSP/ApImngWWpKlUVW+k2TGSpKFV1dGz9Bdw3Ax9Z9AUYJI0I08RlCRJkqSOeARLkiTNiRPOvnzcIUjSyHkES5IkSZI6YoElSZIkSR2xwJIkSZKkjlhgSZIkSVJHLLAkSZIkqSMWWJIkSZLUEQssSZIkSeqIBZYkSZIkdcQCS5IkSZI6YoElSZIkSR2xwJIkSZKkjlhgSZIkSVJHZi2wkqxKsj7JJT1td0tydpIftf/uNrdhSpIkSdLkG+YI1inAYX1trwW+XFX7AV9u1yVJkiRpQZu1wKqqc4Ff9TUfCXyoXf4Q8PRuw5KkLUuya5JPJvlhkkuTPGLcMUmafEkOS3JZkiuS3G6COMlfJ7mwfV2S5NYkd2v7rkxycdu3evTRS5oGi7dxu92r6hqAqromyT06jEmShvEe4EtV9cwkOwJ3HndAkiZbkkXA+4EnAuuA85OcXlU/2DSmqt4BvKMd/1TgL6uqd6L5cVV17QjDljRl5vwmF0lWJFmdZPWGDRvm+uMkLQBJ7gI8BvggQFXdUlW/HmtQkqbBQcAVVbW2qm4BPk5zVs5MjgZOHUlkkuaNbS2wfpFkD4D23/UzDayqlVW1vKqWL126dBs/TpI2c29gA/DPSb6X5OQkO/UPcoJHUp89gat61te1bbeT5M4016B/qqe5gLOSXJBkxZxFKWmqbWuBdTpwbLt8LPDZbsKRpKEsBh4CnFhVDwZ+w4Cb7TjBI6lPBrTVDGOfCnyj7/TAQ6rqIcDhwHFJHjPwQ5zckRa0YW7TfirwLeB+SdYleQHwNuCJSX5Ecx7z2+Y2TEnazDpgXVV9u13/JE3BJUlbsg7Yu2d9L+DqGcYeRd/pgVV1dfvveuA0mlMOb8fJHWlhm/UmF1V19AxdT+g4FkkaSlX9PMlVSe5XVZfR5KMfzLadpAXvfGC/JPsCP6Mpop7dPyjJXYHHAs/tadsJuENV3dguPwl480iiljRVtvUugpI0bi8FPtreQXAt8PwxxyNpwlXVxiTHA2cCi4BVVbUmyYvb/pPaoc8Azqqq3/RsvjtwWhJo9p8+VlVfGl30kqaFBZakqVRVFwLLxx2HpOlSVWcAZ/S1ndS3fgpwSl/bWuCAOQ5P0jww57dplyRJkqSFwgJLkiRJkjpigSVJkiRJHbHAkiRJkqSOWGBJkiRJUkcssCRJkiSpIxZYkiRJktQRn4MlSZI6ccLZl487BEkaO49gSZIkSVJHLLAkSZIkqSMWWJIkSZLUEQssSZIkSeqIBZYkSZIkdcQCS5IkSZI6YoElSZIkSR2xwJIkSZKkjlhgSZpaSRYl+V6Sz487FknTIclhSS5LckWS1w7oPzTJ9UkubF9vGHZbSQJYPO4AJGk7vBy4FLjLuAORNPmSLALeDzwRWAecn+T0qvpB39CvVdVTtnFbSQucR7AkTaUkewFPBk4edyySpsZBwBVVtbaqbgE+Dhw5gm0lLSAWWJKm1buBVwO/n2lAkhVJVidZvWHDhpEFJmli7Qlc1bO+rm3r94gk30/yxSQP3MptJS1wFliSpk6SpwDrq+qCLY2rqpVVtbyqli9dunRE0UmaYBnQVn3r3wXuVVUHAP8IfGYrtm0GOrkjLWgWWJKm0SHA05JcSXOazuOT/Mt4Q5I0BdYBe/es7wVc3Tugqm6oqpva5TOAHZIsGWbbnvdwckdawCywJE2dqvqbqtqrqpYBRwFfqarnjjksSZPvfGC/JPsm2ZEmf5zeOyDJPZOkXT6IZl/pl8NsK0ngXQQlSdICUVUbkxwPnAksAlZV1ZokL277TwKeCbwkyUbgt8BRVVXAwG3H8kUkTTQLLElTrarOAc4ZcxiSpkR72t8ZfW0n9Sy/D3jfsNtKUj9PEZQkSZKkjlhgSZIkSVJHLLAkSZIkqSPbdQ1We4vkG4FbgY1VtbyLoCRJ0mQ74ezLxx2CJE2kLm5y8biquraD95EkSZKkqeYpgpIkSZLUke0tsAo4K8kFSVYMGpBkRZLVSVZv2LBhOz9OkiRJkibX9hZYh1TVQ4DDgeOSPKZ/QFWtrKrlVbV86dKl2/lxkiRJkjS5tqvAqqqr23/XA6cBB3URlCRJkiRNo20usJLslGSXTcvAk4BLugpMkiRJkqbN9txFcHfgtCSb3udjVfWlTqKSJEmSpCm0zQVWVa0FDugwFkmSJEmaat6mXZIkSZI6YoElSZIkSR2xwJIkSZKkjlhgSZIkSVJHLLAkTZ0keyf51ySXJlmT5OXjjkmSJAm27zbtkjQuG4FXVtV32+fxXZDk7Kr6wbgDk+arE86+fNwhdCLJYcB7gEXAyVX1tr7+5wCvaVdvAl5SVd9v+64EbgRuBTZW1fJRxS1pelhgSZo6VXUNcE27fGOSS4E9AQssSTNKsgh4P/BEYB1wfpLT+yZnfgw8tqquS3I4sBI4uKf/cVV17ciCljR1PEVQ0lRLsgx4MPDtMYciafIdBFxRVWur6hbg48CRvQOq6ptVdV27eh6w14hjlDTlLLAkTa0kOwOfAl5RVTcM6F+RZHWS1Rs2bBh9gJImzZ7AVT3r69q2mbwA+GLPegFnJbkgyYqZNjL3SAubBZakqZRkB5ri6qNV9elBY6pqZVUtr6rlS5cuHW2AkiZRBrTVwIHJ42gKrNf0NB9SVQ8BDgeOS/KYQduae6SFzQJL0tRJEuCDwKVV9a5xxyNpaqwD9u5Z3wu4un9Qkv2Bk4Ejq+qXm9qr6ur23/XAaTSnHErSZiywJE2jQ4BjgMcnubB9HTHuoCRNvPOB/ZLsm2RH4Cjg9N4BSfYBPg0cU1WX97Tv1N61lCQ7AU8CLhlZ5JKmhncRlDR1qurrDD7VR5JmVFUbkxwPnElzm/ZVVbUmyYvb/pOANwB3B/6pOVj+h9ux7w6c1rYtBj5WVV8aw9eQNOEssCRJ0oJRVWcAZ/S1ndSz/ELghQO2WwscMOcBSpp6FliSJGkz8+WhwpI0Dl6DJUmSJEkdscCSJEmSpI5YYEmSJElSR7wGS5pyXishSZI0OTyCJUmSJEkd8QiWNEU8WiVpLphbJKk7FljShHAHR5IkafpZYEljYkElSZI0/1hgSSNgMSVJkrQwWGBJc8CCStKkMj9J0tyywJK2kzsrkiRJ2sQCS9pKFlSSJEmaiQWWtAUWU5KmnXlMkkbLAksLljsdkiRJ6poFlhYMCypJ8515TpLGb7sKrCSHAe8BFgEnV9XbOolK2gJ3IATmHwnMh9tittyRJG3/EcDNwJ9V1XeH2VaSYDsKrCSLgPcDTwTWAecnOb2qftBVcFp43FnQMMw/WojMj9tvyNxxOLBf+zoYOBE42LwjaVjbcwTrIOCKqloLkOTjwJGAiWae8f/UNYHMP5oa5tCJMkzuOBL4cFUVcF6SXZPsASwbYltJ2q4Ca0/gqp71dTQzPZtJsgJY0a7elOSyId9/CXDtdsQ3DsY8GsY8An+1dTHfay5jGcD8c3vGPBrGPAJzmH+GyR2Dxuw55LaAuWcKGPNoTF3MW5l7YIb8sz0FVga01e0aqlYCK7f6zZPVVbV8WwIbF2MeDWMejQmP2fzTx5hHw5hHYw5jHiZ3zDRmqLwD5p5JZ8yjsZBj3p4Cax2wd8/6XsDV2xeOJA3F/CNpWwyTO2Yas+MQ20oSd9iObc8H9kuyb5IdgaOA07sJS5K2yPwjaVsMkztOB56XxsOB66vqmiG3laRtP4JVVRuTHA+cSXO70lVVtaazyLbh0PoEMObRMObRmNiYzT8DGfNoGPNozEnMM+WOJC9u+08CzqC5RfsVNLdpf/6Wtu04RP+3Gg1jHo0FG3Oam+RIkiRJkrbX9pwiKEmSJEnqYYElSZIkSR0Ze4GV5LAklyW5IslrB/QnyXvb/ouSPGQccfbFNFvMz2ljvSjJN5McMI44+2LaYsw94x6W5NYkzxxlfDPEMmvMSQ5NcmGSNUm+OuoYB8Qz238bd03yuSTfb2N+/jji7IlnVZL1SS6ZoX/ifn9dMfeMhrlnNKYt97QxmX/MP3PK/DMa05Z/RpJ7qmpsL5qLRP8NuDfN7U+/Dzygb8wRwBdpnj/xcODbUxDzI4Hd2uXDpyHmnnFfobnA95mTHjOwK/ADYJ92/R5TEPPrgLe3y0uBXwE7jjHmxwAPAS6ZoX+ifn8j/t9qor67uWdyYjb3dBa3+cf8M9aYe8aZf+Y25onKP6PIPeM+gnUQcEVVra2qW4CPA0f2jTkS+HA1zgN2TbLHqAPtMWvMVfXNqrquXT2P5lkZ4zTM3xngpcCngPWjDG4Gw8T8bODTVfVTgKoad9zDxFzALkkC7EyTZDaONsyeYKrObWOYyaT9/rpi7hkNc89oTF3uAfOP+WfOmX9GY+ryzyhyz7gLrD2Bq3rW17VtWztmlLY2nhfQVMHjNGvMSfYEngGcNMK4tmSYv/N9gd2SnJPkgiTPG1l0gw0T8/uA+9M8nPJi4OVV9fvRhLdNJu331xVzz2iYe0ZjPuYemLzfYFfMP6Nh/hmN+Zh/tvv3t83PwepIBrT13zd+mDGjNHQ8SR5Hk2QeNacRzW6YmN8NvKaqbm0mGMZumJgXAw8FngD8EfCtJOdV1eVzHdwMhon5/wIuBB4P3Ac4O8nXquqGOY5tW03a768r5p7RMPeMxnzMPTB5v8GumH9Gw/wzGvMx/2z372/cBdY6YO+e9b1oqtutHTNKQ8WTZH/gZODwqvrliGKbyTAxLwc+3iaYJcARSTZW1WdGEuHtDfvfxrVV9RvgN0nOBQ4AxpVkhon5+cDbqjnJ94okPwb+BPjOaELcapP2++uKuWc0zD2jMR9zD0zeb7Ar5p/RMP+MxnzMP9v/+9vai7a6fNEUeGuBfbntwrgH9o15MptfaPadKYh5H5onwD9ynLFuTcx9409h/Bd6DvN3vj/w5XbsnYFLgAdNeMwnAm9ql3cHfgYsGfPfehkzX+g5Ub+/Ef9vNVHf3dwzOTGbezqN3fxj/hlbzH3jzT9zF/PE5Z+5zj1jPYJVVRuTHA+cSXMXklVVtSbJi9v+k2ju6nIEzY/2ZpoqeGyGjPkNwN2Bf2pnRTZW1fIJj3miDBNzVV2a5EvARcDvgZOrauAtNyclZuAtwClJLqb54b6mqq4dV8xJTgUOBZYkWQe8EdgBJvP31xVzz0TFPFHMPaNj/jH/TEDME8X8MxqjyD1pKzVJkiRJ0nYa910EJUmSJGnesMCSJEmSpI5YYEmSJElSRyywJEmSJKkjFliSJEmS1BELrHksye5JPpZkbZILknwryTNGHMNNo/w8SeNn7pE0LuYfTQILrHkqzUMoPgOcW1X3rqqHAkfRPI1akuaEuUfSuJh/NCkssOavxwO39D5Ir6p+UlX/CJBkUZJ3JDk/yUVJXtS2H5rknCSfTPLDJB9tExZJHprkq+2M0JlJ9uj/0CT7trNF5yd5S0972s+7JMnFSZ7Vtu+R5NwkF7Z9j27bT0yyOsmaJH/f8z5HtHF9Pcl7k3y+bd8pyar2c7+X5Mg5+atKmo25R9K4mH80GarK1zx8AS8DTthC/wrgb9vlOwKrgX1pnmx9Pc1szx2AbwGPonnC9TeBpe02z6J5Wnf/+54OPK9dPg64qV3+r8DZNE/53h34KbAH8Erg9e2YRcAu7fLdetrOAfYH7gRcBezb9p0KfL5d/n+B57bLuwKXAzuN+38HX74W2svcY+7x5WtcL/OP+WdSXovRgpDk/TTJ4paqehjwJGD/JM9sh9wV2A+4BfhOVa1rt7sQWAb8GngQcHY7qbMIuGbARx1Ck1AAPgK8vV1+FHBqVd0K/CLJV4GHAecDq5LsAHymqi5sx/9pkhXAYppk9ACapLe2qn7cjjmVJlnSfp+nJXlVu34nYB/g0qH/SJI6Z+6RNC7mH42LBdb8tYbbfuxU1XFJltDM1gAEeGlVndm7UZJDgf/oabqV5r+TAGuq6hFDfHYNaMvAgVXnJnkM8GTgI0neAXwNeBXwsKq6LskpNElj4Hv0vP9/rarLhohP0twx90gaF/OPJoLXYM1fXwHulOQlPW137lk+E3hJO3tCkvsm2WkL73cZsDTJI9rxOyR54IBx36C5oBTgOT3t5wLPas9/Xgo8BvhOknsB66vqA8AHgYcAdwF+A1yfZHfg8PY9fgjcO8mydv1Zfd/npT3nTD94C99F0twx90gaF/OPJoJHsOapqqokTwdOSPJqYAPND/c17ZCTaQ5/f7f9YW4Anr6F97ulPaT+3iR3pflv5900s0W9Xg58LMnLgU/1tJ8GPAL4Ps0sz6ur6udJjgX+OsnvgJtozmH+cZLvte+9liZxUVW/TfIXwJeSXAt8p+f939LGc1H7fa4EnjL7X0pSl8w95h5pXMw/5p9JkapBRzSlyZRk56q6qU0k7wd+VFUnjDsuSfObuUfSuJh/po+nCGra/Hl78ekamotT/+d4w5G0QJh7JI2L+WfKeARLkiRJkjriESxJkiRJ6ogFliRJkiR1xAJLkiRJkjpigSVJkiRJHbHAkiRJkqSO/P8aEe66CZKQzQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "g1 = np.array([Q.random_genotype( sigma=[1,0.2], method='normal' )      for _ in range(100)])\n",
    "g2 = np.array([Q.random_genotype( sigma=[0.05],  method='exponential' ) for _ in range(100)])\n",
    "g3 = np.array([Q.random_genotype( sigma=[0.5,1], method='uniform' )     for _ in range(100)])\n",
    "xx = np.linspace(0,1,51)\n",
    "\n",
    "\n",
    "plt.figure( figsize=(12,4) )\n",
    "plt.subplot(1,3,1)\n",
    "plt.hist( g1.flatten() ,bins=xx, alpha=0.5, density=True)\n",
    "plt.xlabel('Gene dosage')\n",
    "plt.title('Normal(1,0.2)')\n",
    "\n",
    "plt.subplot(1,3,2)\n",
    "plt.hist( g2.flatten() ,bins=xx, alpha=0.5, density=True)\n",
    "plt.xlabel('Gene dosage')\n",
    "plt.title('Exponential(0.05)')\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "plt.hist( g3.flatten() ,bins=xx, alpha=0.5, density=True)\n",
    "plt.xlabel('Gene dosage')\n",
    "plt.title('Uniform(0.5, 1)')\n",
    "\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Computing individuals and populations\n",
    "\n",
    "Imagine that we would like to compute the biomass production of many random genotypes growing in the standard medium. \n",
    "\n",
    "#### 7.0 Computing individuals by \"hand\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1951349383166596\n"
     ]
    }
   ],
   "source": [
    "# Generate a random genotype\n",
    "genotype = np.random.rand(Q.N_GENES)\n",
    "\n",
    "with Q.model:\n",
    "    # Set medium\n",
    "    Q.model.medium = std_medium\n",
    "\n",
    "    # Set genotype\n",
    "    Q.set_genotype( genotype )\n",
    "\n",
    "    # Apply dosage and compute solution\n",
    "    Q.apply_dosage()\n",
    "    print(Q.slim_optimize())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7.1 Computing individuals using `Q.compute_individual()`\n",
    "For robustness, we enter a context manager and set the bounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1951349383166597\n"
     ]
    }
   ],
   "source": [
    "with Q.model:\n",
    "    # We call these for robustness\n",
    "    Q.reset_dosage()\n",
    "    Q.apply_dosage()\n",
    "    \n",
    "    # Compute individual\n",
    "    _mu, _gen, _flx, _med, _wt = Q.compute_individual(genotype, std_medium)\n",
    "print(_mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7.2 Computing a population \"by hand\"\n",
    "and in standard medium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running 10 iterations have taken 1.9 seconds.\n"
     ]
    }
   ],
   "source": [
    "from datetime import datetime\n",
    "\n",
    "NSAMPLES  = 10\n",
    "genotypes = np.random.rand(NSAMPLES, Q.N_GENES)\n",
    "media_list  = [ std_medium for _ in range(NSAMPLES) ]\n",
    "\n",
    "\n",
    "tStart = datetime.now()\n",
    "with Q.model:\n",
    "    Q.reset_dosage()\n",
    "    Q.apply_dosage()\n",
    "    \n",
    "    for jj in range(10):\n",
    "        _mu, _gen, _flx, _med, _wt = Q.compute_individual( genotypes[jj,:] , media_list[jj] )\n",
    "    \n",
    "        \n",
    "tTotal = (datetime.now()-tStart).total_seconds()\n",
    "print('Running %d iterations have taken %1.1f seconds.' % (NSAMPLES, tTotal) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 7.2 Computing a population using `Q.compute_population()` in a parallel\n",
    "and in standard medium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing population of 10 samples...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10/10 [00:00<00:00, 16.65it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total infeasible from within 0\n",
      "Running 10 iterations have taken 0.9 seconds.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "NPROCESSORS = 10\n",
    "\n",
    "tStart = datetime.now()\n",
    "with Q.model:\n",
    "    Q.reset_dosage()\n",
    "    Q.apply_dosage()\n",
    "    \n",
    "    _mu, _gen, _flx, _med, _wt = Q.compute_population( genotypes , media_list, nprocessors = NPROCESSORS)\n",
    "    \n",
    "      \n",
    "tTotal = (datetime.now()-tStart).total_seconds()\n",
    "print('Running %d iterations have taken %1.1f seconds.' % (NSAMPLES, tTotal) )"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}