Untitled1-Copy1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8ccb822e-296f-426c-959c-11be781a6ce6",
   "metadata": {},
   "source": [
    "# Tetracycline Resistance Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f71b51fb-f313-4d53-8755-213680c3db0d",
   "metadata": {},
   "source": [
    "(c) 2021 Tom Röschinger. This work is licensed under a \n",
    "[Creative Commons Attribution License CC-BY 4.0](https://creativecommons.org/licenses/by/4.0/). \n",
    "All code contained herein is licensed under an \n",
    "[MIT license](https://opensource.org/licenses/MIT)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "id": "a1b8fd92-1cd7-4972-9c9e-5c1fe66ddebe",
   "metadata": {},
   "outputs": [],
   "source": [
    "using SymPy, Polynomials, ColorSchemes, Turing, Distributions, LinearAlgebra, resistance_evolution, ForwardDiff, CairoMakie, CSV, DataFrames\n",
    "\n",
    "resistance_evolution.plotting_style.default_makie!()\n",
    "CairoMakie.activate!()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd8ca2ab-5edf-41d8-b0a2-9f0438dff58d",
   "metadata": {},
   "source": [
    "In this notebook we go through everything bla bla"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c04d93f6-a8f4-4d5e-b462-c9d07ae7c42f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(a, λ, λ_0, κ_t, K_d, Δr, j, V_0, a_ex, K_M, P_out, P_in)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@syms a λ λ_0 κ_t K_d Δr j V_0 a_ex K_M P_out P_in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "febf9ac8-498c-44c9-828e-c4bcb2480907",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$- Δr \\left(- \\frac{λ}{λ_{0}} + 1\\right) + \\frac{a λ}{K_{d} κ_{t}}$"
      ],
      "text/plain": [
       "     ⎛  λ     ⎞    a⋅λ  \n",
       "- Δr⋅⎜- ── + 1⎟ + ──────\n",
       "     ⎝  λ₀    ⎠   K_d⋅κₜ"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq1 = (a * λ) / (κ_t * K_d) - Δr * (1 - λ/λ_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dbdb7624-c4af-494d-99f1-d8dc5ea7b303",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$K_{M} P_{in} a_{ex} + a^{2} \\left(- P_{out} - λ\\right) + a \\left(- K_{M} P_{out} - K_{M} λ + P_{in} a_{ex} - \\frac{V_{0} λ}{λ_{0}}\\right)$"
      ],
      "text/plain": [
       "               2                 ⎛                              V₀⋅λ⎞\n",
       "K_M⋅Pᵢₙ⋅aₑₓ + a ⋅(-Pₒᵤₜ - λ) + a⋅⎜-K_M⋅Pₒᵤₜ - K_M⋅λ + Pᵢₙ⋅aₑₓ - ────⎟\n",
       "                                 ⎝                               λ₀ ⎠"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq2 = a^2 * (-λ - P_out) + a * (-λ * K_M  + P_in * a_ex - P_out * K_M - V_0 * λ/λ_0 ) + P_in * a_ex * K_M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2aa1f909-54e0-4065-8d74-9bdc88197a5c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\frac{K_{d} Δr κ_{t} \\left(- λ + λ_{0}\\right)}{λ λ_{0}}$"
      ],
      "text/plain": [
       "K_d⋅Δr⋅κₜ⋅(-λ + λ₀)\n",
       "───────────────────\n",
       "        λ⋅λ₀       "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr1 = solve(eq1, a)[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3be1138e-c9ef-4d09-a1ee-aec4e4cf4ad6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$K_{M} P_{in} a_{ex} + \\frac{K_{d}^{2} Δr^{2} κ_{t}^{2} \\left(- P_{out} - λ\\right) \\left(- λ + λ_{0}\\right)^{2}}{λ^{2} λ_{0}^{2}} + \\frac{K_{d} Δr κ_{t} \\left(- λ + λ_{0}\\right) \\left(- K_{M} P_{out} - K_{M} λ + P_{in} a_{ex} - \\frac{V_{0} λ}{λ_{0}}\\right)}{λ λ_{0}}$"
      ],
      "text/plain": [
       "                                                                        ⎛     \n",
       "                 2   2   2                      2   K_d⋅Δr⋅κₜ⋅(-λ + λ₀)⋅⎜-K_M⋅\n",
       "              K_d ⋅Δr ⋅κₜ ⋅(-Pₒᵤₜ - λ)⋅(-λ + λ₀)                        ⎝     \n",
       "K_M⋅Pᵢₙ⋅aₑₓ + ─────────────────────────────────── + ──────────────────────────\n",
       "                              2   2                                           \n",
       "                             λ ⋅λ₀                                            \n",
       "\n",
       "                         V₀⋅λ⎞\n",
       "Pₒᵤₜ - K_M⋅λ + Pᵢₙ⋅aₑₓ - ────⎟\n",
       "                          λ₀ ⎠\n",
       "──────────────────────────────\n",
       "λ⋅λ₀                          \n",
       "                              "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr2 = subs(eq2, a=>expr1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c013e50e-0092-4bf1-910d-b54751b21aa1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$λ^{2} λ_{0}^{2} \\left(K_{M} P_{in} a_{ex} + \\frac{K_{d}^{2} Δr^{2} κ_{t}^{2} \\left(- P_{out} - λ\\right) \\left(- λ + λ_{0}\\right)^{2}}{λ^{2} λ_{0}^{2}} + \\frac{K_{d} Δr κ_{t} \\left(- λ + λ_{0}\\right) \\left(- K_{M} P_{out} - K_{M} λ + P_{in} a_{ex} - \\frac{V_{0} λ}{λ_{0}}\\right)}{λ λ_{0}}\\right)$"
      ],
      "text/plain": [
       "       ⎛                                                                      \n",
       "       ⎜                 2   2   2                      2   K_d⋅Δr⋅κₜ⋅(-λ + λ₀\n",
       " 2   2 ⎜              K_d ⋅Δr ⋅κₜ ⋅(-Pₒᵤₜ - λ)⋅(-λ + λ₀)                      \n",
       "λ ⋅λ₀ ⋅⎜K_M⋅Pᵢₙ⋅aₑₓ + ─────────────────────────────────── + ──────────────────\n",
       "       ⎜                              2   2                                   \n",
       "       ⎝                             λ ⋅λ₀                                    \n",
       "\n",
       "  ⎛                              V₀⋅λ⎞⎞\n",
       ")⋅⎜-K_M⋅Pₒᵤₜ - K_M⋅λ + Pᵢₙ⋅aₑₓ - ────⎟⎟\n",
       "  ⎝                               λ₀ ⎠⎟\n",
       "──────────────────────────────────────⎟\n",
       "        λ⋅λ₀                          ⎟\n",
       "                                      ⎠"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr3 = expr2 * λ^2 * λ_0^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7b591a88-cfab-4ebe-9bf3-b13e6e50be4c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$K_{M} K_{d} P_{out} Δr κ_{t} λ^{2} λ_{0} - K_{M} K_{d} P_{out} Δr κ_{t} λ λ_{0}^{2} + K_{M} K_{d} Δr κ_{t} λ^{3} λ_{0} - K_{M} K_{d} Δr κ_{t} λ^{2} λ_{0}^{2} + K_{M} P_{in} a_{ex} λ^{2} λ_{0}^{2} - K_{d}^{2} P_{out} Δr^{2} κ_{t}^{2} λ^{2} + 2 K_{d}^{2} P_{out} Δr^{2} κ_{t}^{2} λ λ_{0} - K_{d}^{2} P_{out} Δr^{2} κ_{t}^{2} λ_{0}^{2} - K_{d}^{2} Δr^{2} κ_{t}^{2} λ^{3} + 2 K_{d}^{2} Δr^{2} κ_{t}^{2} λ^{2} λ_{0} - K_{d}^{2} Δr^{2} κ_{t}^{2} λ λ_{0}^{2} - K_{d} P_{in} a_{ex} Δr κ_{t} λ^{2} λ_{0} + K_{d} P_{in} a_{ex} Δr κ_{t} λ λ_{0}^{2} + K_{d} V_{0} Δr κ_{t} λ^{3} - K_{d} V_{0} Δr κ_{t} λ^{2} λ_{0}$"
      ],
      "text/plain": [
       "                    2                             2                  3        \n",
       "K_M⋅K_d⋅Pₒᵤₜ⋅Δr⋅κₜ⋅λ ⋅λ₀ - K_M⋅K_d⋅Pₒᵤₜ⋅Δr⋅κₜ⋅λ⋅λ₀  + K_M⋅K_d⋅Δr⋅κₜ⋅λ ⋅λ₀ - K_\n",
       "\n",
       "             2   2                2   2      2        2   2  2        2       \n",
       "M⋅K_d⋅Δr⋅κₜ⋅λ ⋅λ₀  + K_M⋅Pᵢₙ⋅aₑₓ⋅λ ⋅λ₀  - K_d ⋅Pₒᵤₜ⋅Δr ⋅κₜ ⋅λ  + 2⋅K_d ⋅Pₒᵤₜ⋅Δ\n",
       "\n",
       " 2   2           2        2   2   2      2   2   2  3        2   2   2  2     \n",
       "r ⋅κₜ ⋅λ⋅λ₀ - K_d ⋅Pₒᵤₜ⋅Δr ⋅κₜ ⋅λ₀  - K_d ⋅Δr ⋅κₜ ⋅λ  + 2⋅K_d ⋅Δr ⋅κₜ ⋅λ ⋅λ₀ -\n",
       "\n",
       "    2   2   2     2                      2                            2       \n",
       " K_d ⋅Δr ⋅κₜ ⋅λ⋅λ₀  - K_d⋅Pᵢₙ⋅aₑₓ⋅Δr⋅κₜ⋅λ ⋅λ₀ + K_d⋅Pᵢₙ⋅aₑₓ⋅Δr⋅κₜ⋅λ⋅λ₀  + K_d⋅\n",
       "\n",
       "          3                 2   \n",
       "V₀⋅Δr⋅κₜ⋅λ  - K_d⋅V₀⋅Δr⋅κₜ⋅λ ⋅λ₀"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr4 = factor(expr3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6b245d82-2d91-4a63-9dfb-7ca833aadf29",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$- K_{d}^{2} P_{out} Δr^{2} κ_{t}^{2} λ_{0}^{2} + λ^{3} \\left(K_{M} K_{d} Δr κ_{t} λ_{0} - K_{d}^{2} Δr^{2} κ_{t}^{2} + K_{d} V_{0} Δr κ_{t}\\right) + λ^{2} \\left(K_{M} K_{d} P_{out} Δr κ_{t} λ_{0} - K_{M} K_{d} Δr κ_{t} λ_{0}^{2} + K_{M} P_{in} a_{ex} λ_{0}^{2} - K_{d}^{2} P_{out} Δr^{2} κ_{t}^{2} + 2 K_{d}^{2} Δr^{2} κ_{t}^{2} λ_{0} - K_{d} P_{in} a_{ex} Δr κ_{t} λ_{0} - K_{d} V_{0} Δr κ_{t} λ_{0}\\right) + λ \\left(- K_{M} K_{d} P_{out} Δr κ_{t} λ_{0}^{2} + 2 K_{d}^{2} P_{out} Δr^{2} κ_{t}^{2} λ_{0} - K_{d}^{2} Δr^{2} κ_{t}^{2} λ_{0}^{2} + K_{d} P_{in} a_{ex} Δr κ_{t} λ_{0}^{2}\\right)$"
      ],
      "text/plain": [
       "     2        2   2   2    3 ⎛                      2   2   2               ⎞ \n",
       "- K_d ⋅Pₒᵤₜ⋅Δr ⋅κₜ ⋅λ₀  + λ ⋅⎝K_M⋅K_d⋅Δr⋅κₜ⋅λ₀ - K_d ⋅Δr ⋅κₜ  + K_d⋅V₀⋅Δr⋅κₜ⎠ \n",
       "\n",
       "   2 ⎛                                        2                 2      2      \n",
       "+ λ ⋅⎝K_M⋅K_d⋅Pₒᵤₜ⋅Δr⋅κₜ⋅λ₀ - K_M⋅K_d⋅Δr⋅κₜ⋅λ₀  + K_M⋅Pᵢₙ⋅aₑₓ⋅λ₀  - K_d ⋅Pₒᵤₜ⋅\n",
       "\n",
       "  2   2        2   2   2                                            ⎞     ⎛   \n",
       "Δr ⋅κₜ  + 2⋅K_d ⋅Δr ⋅κₜ ⋅λ₀ - K_d⋅Pᵢₙ⋅aₑₓ⋅Δr⋅κₜ⋅λ₀ - K_d⋅V₀⋅Δr⋅κₜ⋅λ₀⎠ + λ⋅⎝- K\n",
       "\n",
       "                    2        2        2   2         2   2   2   2             \n",
       "_M⋅K_d⋅Pₒᵤₜ⋅Δr⋅κₜ⋅λ₀  + 2⋅K_d ⋅Pₒᵤₜ⋅Δr ⋅κₜ ⋅λ₀ - K_d ⋅Δr ⋅κₜ ⋅λ₀  + K_d⋅Pᵢₙ⋅aₑ\n",
       "\n",
       "          2⎞\n",
       "ₓ⋅Δr⋅κₜ⋅λ₀ ⎠"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr5 = collect(expr4, λ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "078ea630-8fb0-46ed-a939-689818dfbbe3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(expr5.coeff(λ, 4))\n",
    "expr5.coeff(λ, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "67d1c776-a0c5-4d91-8ed6-854deb76d14f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K_M*K_d*Δr*κ_t*λ_0 - K_d^2*Δr^2*κ_t^2 + K_d*V_0*Δr*κ_t\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$- K_{d} Δr κ_{t} \\left(- K_{M} λ_{0} + K_{d} Δr κ_{t} - V_{0}\\right)$"
      ],
      "text/plain": [
       "-K_d⋅Δr⋅κₜ⋅(-K_M⋅λ₀ + K_d⋅Δr⋅κₜ - V₀)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(expr5.coeff(λ, 3))\n",
    "factor(expr5.coeff(λ, 3), K_d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e148c2b9-f245-490d-8fab-41f64f201b48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K_M*K_d*P_out*Δr*κ_t*λ_0 - K_M*K_d*Δr*κ_t*λ_0^2 + K_M*P_in*a_ex*λ_0^2 - K_d^2*P_out*Δr^2*κ_t^2 + 2*K_d^2*Δr^2*κ_t^2*λ_0 - K_d*P_in*a_ex*Δr*κ_t*λ_0 - K_d*V_0*Δr*κ_t*λ_0\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$K_{M} P_{in} a_{ex} λ_{0}^{2} - K_{d}^{2} \\left(P_{out} Δr^{2} κ_{t}^{2} - 2 Δr^{2} κ_{t}^{2} λ_{0}\\right) - K_{d} \\left(- K_{M} P_{out} Δr κ_{t} λ_{0} + K_{M} Δr κ_{t} λ_{0}^{2} + P_{in} a_{ex} Δr κ_{t} λ_{0} + V_{0} Δr κ_{t} λ_{0}\\right)$"
      ],
      "text/plain": [
       "              2      2 ⎛       2   2       2   2   ⎞       ⎛                  \n",
       "K_M⋅Pᵢₙ⋅aₑₓ⋅λ₀  - K_d ⋅⎝Pₒᵤₜ⋅Δr ⋅κₜ  - 2⋅Δr ⋅κₜ ⋅λ₀⎠ - K_d⋅⎝-K_M⋅Pₒᵤₜ⋅Δr⋅κₜ⋅λ₀\n",
       "\n",
       "               2                                 ⎞\n",
       " + K_M⋅Δr⋅κₜ⋅λ₀  + Pᵢₙ⋅aₑₓ⋅Δr⋅κₜ⋅λ₀ + V₀⋅Δr⋅κₜ⋅λ₀⎠"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(expr5.coeff(λ, 2))\n",
    "factor(expr5.coeff(λ, 2), K_d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "fd88f485-4b40-44d3-badf-a2db6916f4b7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-K_M*K_d*P_out*Δr*κ_t*λ_0^2 + 2*K_d^2*P_out*Δr^2*κ_t^2*λ_0 - K_d^2*Δr^2*κ_t^2*λ_0^2 + K_d*P_in*a_ex*Δr*κ_t*λ_0^2\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$K_{d} Δr κ_{t} λ_{0} \\left(- K_{M} P_{out} λ_{0} + K_{d} \\left(2 P_{out} Δr κ_{t} - Δr κ_{t} λ_{0}\\right) + P_{in} a_{ex} λ_{0}\\right)$"
      ],
      "text/plain": [
       "K_d⋅Δr⋅κₜ⋅λ₀⋅(-K_M⋅Pₒᵤₜ⋅λ₀ + K_d⋅(2⋅Pₒᵤₜ⋅Δr⋅κₜ - Δr⋅κₜ⋅λ₀) + Pᵢₙ⋅aₑₓ⋅λ₀)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(expr5.coeff(λ, 1))\n",
    "factor(expr5.coeff(λ, 1), K_d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "048d2217-6ca1-47c2-b842-468accbca45d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-K_d^2*P_out*Δr^2*κ_t^2*λ_0^2\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$- K_{d}^{2} P_{out} Δr^{2} κ_{t}^{2} λ_{0}^{2}$"
      ],
      "text/plain": [
       "    2        2   2   2\n",
       "-K_d ⋅Pₒᵤₜ⋅Δr ⋅κₜ ⋅λ₀ "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(expr5.coeff(λ, 0))\n",
    "expr5.coeff(λ, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ee7694f2-3671-42b6-84b1-e9bd88711b21",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$- K_{d}^{2} P_{out} Δr^{2} κ_{t}^{2} λ_{0}^{2} + λ^{3} \\left(K_{M} K_{d} Δr κ_{t} λ_{0} - K_{d}^{2} Δr^{2} κ_{t}^{2} + K_{d} V_{0} Δr κ_{t}\\right) + λ^{2} \\left(K_{M} K_{d} P_{out} Δr κ_{t} λ_{0} - K_{M} K_{d} Δr κ_{t} λ_{0}^{2} + K_{M} P_{in} a_{ex} λ_{0}^{2} - K_{d}^{2} P_{out} Δr^{2} κ_{t}^{2} + 2 K_{d}^{2} Δr^{2} κ_{t}^{2} λ_{0} - K_{d} P_{in} a_{ex} Δr κ_{t} λ_{0} - K_{d} V_{0} Δr κ_{t} λ_{0}\\right) + λ \\left(- K_{M} K_{d} P_{out} Δr κ_{t} λ_{0}^{2} + 2 K_{d}^{2} P_{out} Δr^{2} κ_{t}^{2} λ_{0} - K_{d}^{2} Δr^{2} κ_{t}^{2} λ_{0}^{2} + K_{d} P_{in} a_{ex} Δr κ_{t} λ_{0}^{2}\\right)$"
      ],
      "text/plain": [
       "     2        2   2   2    3 ⎛                      2   2   2               ⎞ \n",
       "- K_d ⋅Pₒᵤₜ⋅Δr ⋅κₜ ⋅λ₀  + λ ⋅⎝K_M⋅K_d⋅Δr⋅κₜ⋅λ₀ - K_d ⋅Δr ⋅κₜ  + K_d⋅V₀⋅Δr⋅κₜ⎠ \n",
       "\n",
       "   2 ⎛                                        2                 2      2      \n",
       "+ λ ⋅⎝K_M⋅K_d⋅Pₒᵤₜ⋅Δr⋅κₜ⋅λ₀ - K_M⋅K_d⋅Δr⋅κₜ⋅λ₀  + K_M⋅Pᵢₙ⋅aₑₓ⋅λ₀  - K_d ⋅Pₒᵤₜ⋅\n",
       "\n",
       "  2   2        2   2   2                                            ⎞     ⎛   \n",
       "Δr ⋅κₜ  + 2⋅K_d ⋅Δr ⋅κₜ ⋅λ₀ - K_d⋅Pᵢₙ⋅aₑₓ⋅Δr⋅κₜ⋅λ₀ - K_d⋅V₀⋅Δr⋅κₜ⋅λ₀⎠ + λ⋅⎝- K\n",
       "\n",
       "                    2        2        2   2         2   2   2   2             \n",
       "_M⋅K_d⋅Pₒᵤₜ⋅Δr⋅κₜ⋅λ₀  + 2⋅K_d ⋅Pₒᵤₜ⋅Δr ⋅κₜ ⋅λ₀ - K_d ⋅Δr ⋅κₜ ⋅λ₀  + K_d⋅Pᵢₙ⋅aₑ\n",
       "\n",
       "          2⎞\n",
       "ₓ⋅Δr⋅κₜ⋅λ₀ ⎠"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "69dd2108-d567-4948-b668-e1af749388a4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "solve_polynomial (generic function with 1 method)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "function solve_polynomial(aex, λ0, κt, Kd, KM, j, V0, Δr)\n",
    "    c = zeros(4)\n",
    "    c1 = κt * Kd * Δr\n",
    "    c[4] = c1 * (KM * λ0 - c1 + V0)\n",
    "    c[3] = KM * j * aex * λ0^2 - c1^2 * (j - 2λ0) - c1 * λ0 * (-KM * j  + KM * λ0 + j * aex + V0)\n",
    "    c[2] = c1 * λ0 * (-KM * j * λ0 + c1 * (2j - λ0) + j * aex * λ0)\n",
    "    c[1] = -Kd^2*j*Δr^2*κt^2*λ0^2\n",
    "    pol = Polynomial(c)\n",
    "    return Polynomials.roots(pol)\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "f95af6c4-0aad-4b31-a725-7167166a5810",
   "metadata": {},
   "outputs": [],
   "source": [
    "_aex = 0  # µM\n",
    "_λ0 = 0.68  # h**-1\n",
    "_κt = 0.06  # µM**-1 h**-1\n",
    "_Kd = .1  # µM\n",
    "_KM = 100 # µM\n",
    "_j = 10# h**-1\n",
    "Pin = 10\n",
    "Pout = 3\n",
    "_V0 = 0\n",
    "_Δr = 46.5  # µM\n",
    "\n",
    "\n",
    "a_ex_range = [0, 0.1, 0.26, 0.4, 0.6, 1]\n",
    "sol_list = zeros(ComplexF64, 3, length(a_ex_range))\n",
    "for (i, x) in enumerate(a_ex_range)\n",
    "    # Pack parameters together\n",
    "    args = (x, _λ0, _κt, _Kd, _KM, _j, _V0, _Δr)\n",
    "    # Find roots\n",
    "    sol_list[:, i] = solve_polynomial(args...)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "id": "9b6fbeec-2494-4325-a50e-832c83aa72a0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{ComplexF64}:\n",
       "    -12.309686480501533 + 0.0im\n",
       " -0.0028014943665923427 + 0.0im\n",
       "      0.552410494838448 + 0.0im"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_list[:, 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "id": "72c8e5dd-dcbf-470f-aa71-ddc3a175a8af",
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "UndefVarError: scatter not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: scatter not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope",
      "   @ In[120]:1",
      " [2] eval",
      "   @ ./boot.jl:368 [inlined]",
      " [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
      "   @ Base ./loading.jl:1281"
     ]
    }
   ],
   "source": [
    "p = scatter(\n",
    "    ylim=[0, 1],\n",
    "    xscale=:log,\n",
    "    legendtitle = \"V0\",\n",
    "    ylabel=\"λ/λ0\",\n",
    "    xlabel=\"a_ex\",\n",
    "    title=\"Is this the right font?\"\n",
    ")\n",
    "\n",
    "V_0_list = [0, 1000, 5000, 10000]\n",
    "for (_V0, c) in zip(V_0_list, palette(:BuPu_6)[2:end])\n",
    "    y = []\n",
    "    x = Float64[]\n",
    "    for i in 1:length(a_ex_range)\n",
    "        args = (a_ex_range[i], _λ0, _κt, _Kd, _KM, _j, _V0, _Δr)\n",
    "        solutions = solve_polynomial(args...)\n",
    "        _y = [imag(x) == 0 ? real(x)/_λ0 : missing for x in solutions]\n",
    "        push!(y, _y...)\n",
    "        _x = ones(length(_y)) .* a_ex_range[i]\n",
    "        push!(x, _x...)\n",
    "    end\n",
    "    scatter!(p, x, y, color=c, label=\"$_V0\")\n",
    "end\n",
    "savefig(p, \"res_model.pdf\")\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38d60a49-3233-4633-aa1e-0b64df2668cc",
   "metadata": {},
   "source": [
    "## Writing Generative Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9735d3b0-0ce5-4125-8472-ba23259cdd5e",
   "metadata": {},
   "source": [
    "Now we are going to write down the model in a Bayesian way and try out various priors.\n",
    "\n",
    "\\begin{align}\n",
    "\\lambda_0 &= 0.68 h^{-1},\\\\[.5em]\n",
    "\\kappa_t &= 0.06 \\mu M^{-1} h^{-1},\\\\[.5em]\n",
    "K_d &= 0.1 \\mu M,\\\\[.5em]\n",
    "K_M &= 10 \\mu M,\\\\[.5em]\n",
    "\\Delta r &= 46.5 \\mu M,\\\\[.5em]\n",
    "\\log j &\\sim \\text{Norm}(2, 2),\\\\[.5em]\n",
    "j &= 10^{\\log j} \\times 1 s^{-1},\\\\[.5em]\n",
    "\\log P_\\mathrm{in} &\\sim \\text{Norm}(2, 2),\\\\[.5em]\n",
    "P_\\mathrm{in} &= 10^{\\log P_\\mathrm{in}} \\times 1 s^{-1},\\\\[.5em]\n",
    "\\log P_\\mathrm{out} &\\sim \\text{Norm}(2, 2),\\\\[.5em]\n",
    "P_\\mathrm{out} &= 10^{\\log P_\\mathrm{out}} \\times 1 s^{-1},\\\\[.5em]\n",
    "\\log V_0 &\\sim \\text{Norm}(2, 2),\\\\[.5em]\n",
    "V_0  &= 10^{\\log V_0 } \\times 1 \\mu M s^{-1}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "id": "57338bcc-1378-46c2-a110-29f1b42e6ba0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fit_growth (generic function with 4 methods)"
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@model function fit_growth(growth_rates, aex_arr, aex_ppc)\n",
    "    \n",
    "    y_ppc = zeros(Float64, length(aex_ppc))\n",
    "\n",
    "    # Priors\n",
    "    log_j ~ Normal(2, 2)\n",
    "    log_V0 ~ Normal(2, 2)\n",
    "    log_sigma ~ Normal(-3, 1)\n",
    "    \n",
    "    # Transform Parameters\n",
    "    sigma = float(10^log_sigma)\n",
    "    j = 10^log_j\n",
    "    V0 = 10^log_V0\n",
    "    \n",
    "    # Fixed Parameters\n",
    "    λ0 = 0.68  \n",
    "    κt = 0.06\n",
    "    Kd = .1\n",
    "    KM = 10\n",
    "    Δr = 46.5\n",
    "    \n",
    "    for (i, aex) in enumerate(aex_arr)\n",
    "        if typeof(V0) == ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}\n",
    "            x = [aex, λ0, κt, Kd, KM, ForwardDiff.value(j), ForwardDiff.value(V0), Δr]\n",
    "        else\n",
    "            x = [aex, λ0, κt, Kd, KM, j, V0, Δr]\n",
    "        end\n",
    "        solutions = solve_polynomial(x...)\n",
    "        \n",
    "        y = [imag(x) == 0 ? real(x)/λ0 : missing for x in solutions]\n",
    "        y = filter(x -> !ismissing(x) && (x > 0), y)\n",
    "        if length(y) == 3\n",
    "            deleteat!(y, 2)\n",
    "            opt_y = y[argmin([(λ - root)^2 for root in y])]\n",
    "        else\n",
    "            opt_y = y[1]\n",
    "       end\n",
    "        \n",
    "\n",
    "        growth_rates[i] ~ Normal(opt_y, sigma)\n",
    "    end\n",
    "    for (i, aex) in enumerate(aex_ppc)\n",
    "        if typeof(V0) == ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}\n",
    "            x = [aex, λ0, κt, Kd, KM, ForwardDiff.value(j), ForwardDiff.value(V0), Δr]\n",
    "        else\n",
    "            x = [aex, λ0, κt, Kd, KM, j, V0, Δr]\n",
    "        end\n",
    "        solutions = solve_polynomial(x...)\n",
    "        \n",
    "        y = [imag(x) == 0 ? real(x)/λ0 : missing for x in solutions]\n",
    "        y = filter(x -> !ismissing(x) && (x > 0), y)\n",
    "        if length(y) == 3\n",
    "            deleteat!(y, 2)\n",
    "            opt_y = y[argmin([(λ - root)^2 for root in y])]\n",
    "        else\n",
    "            opt_y = y[1]\n",
    "        end\n",
    "        if typeof(sigma) == ForwardDiff.Dual{ForwardDiff.Tag{Turing.TuringTag, Float64}, Float64, 3}\n",
    "            y_ppc[i] = rand(Normal(opt_y / λ0, ForwardDiff.value(sigma)), 1)[1]\n",
    "        else\n",
    "            y_ppc[i] = rand(Normal(opt_y, sigma), 1)[1]\n",
    "        end\n",
    "\n",
    "    end\n",
    "    return V0, j, sigma, y_ppc\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "id": "681f4bb3-6176-45e4-a180-896ca98c983b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-9.999999999999993, -0.002801494366592343, 0.6800000000000003]\n",
      "[-12.309686480501533, -0.0028014943665923427, 0.552410494838448]\n",
      "[-16.0798081795807, -0.0028014943665923436, 0.4228906168566823]\n",
      "[-19.419271557329584, -0.0028014943665923427, 0.3501676146772568]\n",
      "[-24.22436502617409, -0.0028014943665923427, 0.2807091121956212]\n",
      "[-33.893388953043114, -0.0028014943665923427, 0.20062909641230958]\n"
     ]
    }
   ],
   "source": [
    "y = Float64[]\n",
    "for i in 1:length(a_ex_range)\n",
    "    args = (a_ex_range[i], _λ0, _κt, _Kd, _KM, _j, 0, _Δr)\n",
    "    solutions = solve_polynomial(args...)\n",
    "    println(solutions)\n",
    "    _y = [imag(x) == 0 ? real(x) : missing for x in solutions]\n",
    "    push!(y, _y[3])\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "id": "c96c1b38",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.6800000000000003, 0.552410494838448, 0.4228906168566823, 0.3501676146772568, 0.2807091121956212, 0.20062909641230958]\n"
     ]
    }
   ],
   "source": [
    "println(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "id": "9a280afe-43ce-4379-8102-a4cdc0f6a08a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Found initial step size\n",
      "│   ϵ = 0.4\n",
      "└ @ Turing.Inference /Users/tomroeschinger/.julia/packages/Turing/S4Y4B/src/inference/hmc.jl:188\n",
      "\u001b[32mSampling: 100%|█████████████████████████████████████████| Time: 0:00:03\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1000×1 Matrix{Tuple{Float64, Float64, Float64, Vector{Float64}}}:\n",
       " (174.1290183867386, 31.874618172280446, 0.24942362560416556, [1.4709298906318506, 0.936139745349281, 0.10219344464106694, 0.34672555851941034, 0.42045604740067466, 0.17990436747677926])\n",
       " (55.85772837153071, 158.85534418392677, 0.22954920547333937, [0.8863023209069691, 1.01421182276424, 0.49191596870647136, 0.5067035374771646, 0.32074984820112734, 0.5774356906669439])\n",
       " (0.18242382651841124, 327.9587831574582, 0.22156181454578242, [0.9437184840165844, 0.8384417484111363, 0.9544609569628102, 0.4918368890842335, 0.3388227514245791, 0.5707538469314668])\n",
       " (0.17642235126784842, 36.102619074545125, 0.18930307944134567, [0.8720148194154137, 1.0230908560478258, 0.3931006102587782, 0.6685835616282074, 0.46309628255463975, 0.2889525898743723])\n",
       " (9.527664649315895, 71.6534521006404, 0.2153472834500859, [1.0007873051560066, 0.7368232124515822, 0.38917934614752436, 0.048500049243827115, 0.4713118214566264, 0.2938498765258395])\n",
       " (31.01176604660996, 90.66970503762195, 0.23089080709222795, [0.9296630778970681, 0.614958541817575, 0.6986569691612498, 0.7003036633874575, -0.10053370289167607, 0.16001132223187212])\n",
       " (6.725637118623289, 91.06375354314206, 0.2651711321420357, [0.5145287911770373, 0.6666418680263252, 0.27878736414557054, 0.4958238456749066, 0.03573106946979798, 0.1708992706663609])\n",
       " (128.46719305213503, 32426.384806621907, 0.17634021287937784, [0.9689843382303719, 0.5506862196961078, 0.6179119443330878, 0.8234031151750509, 0.5347077816291013, 0.3666930625040402])\n",
       " (0.02741040391317572, 55.24804834857196, 0.17133141781318256, [1.0501349264195827, 0.8010875127152678, 0.9208831865573821, 0.7543233713902062, 0.20299866285418558, 0.3992548254165417])\n",
       " (0.3110448517800365, 395.23347757261854, 0.22685265827701745, [0.6117242664655635, 0.8567750842175312, 0.6428988155059293, 0.7624902580725167, 0.526775455071108, 0.12023744808214057])\n",
       " (24.55630996333377, 19.044280170396654, 0.2031240472482805, [1.0746994739279583, 0.6339112330455101, 0.5074905086651572, 0.7436405111869471, 0.26696365192674065, 0.4502270828378483])\n",
       " (6.286005241765871, 32.384936402793706, 0.19863807271836592, [1.2553208767777075, 0.9215867328909604, 0.9658221512079115, 0.5124799113459934, 0.4693885194129933, 0.14726456388650003])\n",
       " (2.8871660633922507, 4.302197978994898, 0.18187072786443667, [0.904079462365157, 0.5780400926115608, 0.6134165808648734, 0.48627082625706186, 0.6647745546748742, 0.47216151804139384])\n",
       " ⋮\n",
       " (0.009939905845205838, 9.952301345449062, 0.16917996411773992, [1.0970926038722908, 0.8796805764439554, 0.644900539474085, 0.3744242143974218, 0.28582437743443745, 0.5030881625237015])\n",
       " (15.444471601670855, 264.336095122655, 0.13484810669273223, [0.8269295450890077, 0.8384523651708307, 0.8296196803905768, 0.7186833834158313, 0.33287488280739463, 0.24549710897553775])\n",
       " (2.5977685095740948, 41.861930618380754, 0.12501019131841046, [0.9097557659132102, 0.9121203940879243, 1.038528165708378, 0.40072769486497556, 0.4858557842269716, 0.2577625435250624])\n",
       " (9.357819541002764, 29.61800039009533, 0.14652360293213115, [0.8736946270298649, 0.9434627863774645, 0.5296601561180283, 0.6017180299791307, 0.12933900466265458, 0.24465346120955017])\n",
       " (1124.725019835313, 1525.1585577588166, 0.32802470867854466, [1.460352779666127, 0.8106139112340004, 0.808987443771909, 0.7947961894423752, 0.3853802161516232, 0.5794042825867578])\n",
       " (91.31541824573561, 22450.453601990583, 0.27404251065541124, [1.0860839786959908, 0.3907434130276563, 0.43635605629624596, 0.32784085658164874, -0.10759181196913875, 0.8440616570800642])\n",
       " (3390.042970131226, 1.660799094246383e7, 0.19955784451475475, [1.0563267665703606, 0.7025259969440805, 0.6220328866163505, 0.8587151853527574, 0.17670681584646436, 0.2436615354948573])\n",
       " (373.55821130348625, 3650.9511837041587, 0.24948140455796605, [1.0082576493992796, 0.34090671618863544, 0.5327004083393363, 0.7213165493518563, 0.2512400751559964, 0.04560712460153504])\n",
       " (469.6795894086615, 31630.44459887798, 0.26343194592176367, [0.8276996967028182, 0.7837776936533081, 0.4512472128376189, 0.6083373914695498, 0.007891944221569747, 0.548147639168734])\n",
       " (634.9317382078159, 31087.04142681108, 0.2952129456805483, [0.9673076683663587, 0.6512033560294714, 0.7241466641336552, 0.7581330303107501, 0.9579720995564633, 0.34840629895335545])\n",
       " (0.010186841931134307, 13692.316277451102, 0.14565392816061404, [0.8128529315259567, 0.6899876712443191, 0.48890978491592696, 0.6174233099032729, 0.619386859033698, 0.4712673881289403])\n",
       " (213985.2277450963, 5154.737503304414, 0.26079471593840803, [1.191205075527066, 1.5724158600475704, 1.0125150714733555, 1.0699064341633662, 0.5712524936032315, 0.3924336729189202])"
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chn = sample(fit_growth(y, a_ex_range, a_ex_range), NUTS(0.65), 1000)\n",
    "chains_params = Turing.MCMCChains.get_sections(chn, :parameters)\n",
    "gen = generated_quantities(fit_growth(y, a_ex_range, a_ex_range), chains_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "id": "b29fa9e4-83d7-44a6-a386-a0beadcf4c08",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Vector{Float64}:\n",
       " 0.8207717098019167\n",
       " 0.6041826605226885\n",
       " 0.5164043525282879\n",
       " 0.41925971636038145\n",
       " 0.30429665699064173"
      ]
     },
     "execution_count": 164,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_ppc = mean(hcat([x[end] for x in gen]...), dims=2)[2:end] |> vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "id": "46154f1a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVxU593//3NmYRa2QWBEARFEFlmk1ZJEGyUiasU1LtEQF9R417TpnebXJtrkURNz/6r33XxTv0lvk5tGQBONUdGIS8QtEqOJsTY6gAiCigubiMi+8/1jUuutJDo4zJnl9fwjDzgzc+bTw3HeveY6n3OJXV1dAgAAjkomdQEAAEiJIAQAODSCEADg0AhCAIBDIwgBAA6NIAQAODSCEADg0AhCAIBDIwgBAA6NIAQAODSCEADg0AhCAIBDIwgBAA6NIAQAODRpgvDo0aNHjx6V5K0BALibQpJ3PXr0aG19Q3DYkEffVT+91+3a+sbm5kffFayHVqNWOzlV366VuhCr0MfdTRAEjoadkclkPl59SiurpC7EKmjVanc3lzIzHQ0/H71Jz+erUQCAQyMIAQAOjSAEADg0ghAA4NAIQgCAQyMIAQAOjSAEADg0ghAA4NAIQgCAQyMIAQAOjSAEADg0ghAA4NAIQgCAQyMIAQAOzeaDsORW099OXe/o7JK6EACATZJmPUKzqKhvPVB06+vrTQcKb/zt22vrpgwe4q2VuigAgI2x1SDMzL/5/+2/2Nzeafy1uLp5QnrOvJi+r8UN0CptfpgLALAYm8yMwqqmf//8Xylo1CUI2/OqEtINxy7flqowAIDNsckg3GSobO/ovH97e2fH6IHuz20vWJZZVN3UbvnCAAA2xyaDsLCqsdvtrR1C0tC+GXPDC282jv7wzCZDpYULAwDYHJsMQqXsB8t2kovDfV2z5kcui+3/x0OX520vuHa7xZK1AQBsi00G4U/7u3S73UOjCPRQC4KgkMteeKz/4eTo9s7Op1IN606W0l8BAOiWTQbh/J/09XJW3r/9N0/4KmTinV8Heqg3zwpfPS7wg1NlEz/KPVveYMEaAQC2oSdB2Nra2u325ubmRyvmYfXRKD6dHR6h/1fXoEYpWz7Kf8kwn3ueKYrCzAivw8nR4V7aaZvyVh4paWjtsEyRAACbYHIf4dKlS3NzcxsbG9PT02NiYowbW1tbZ86cef369bq6unfffXfChAnmrvNeIV6az+dHFVU33exQqoX2YHeli0r+Q0/2dlauTRw0ObzPawcvHyy6tXpc4OiB7r1dIQDAJpg2IszOzq6srDxx4kRKSsorr7xyZ/vevXs1Gs3p06czMjJef/11cxfZPVEUBntqpkf6jBig+5EUvCM+yOOLRUMnh3kuyChYlll0s5H+CgCAiSPCY8eOGUd7sbGx+fn5d7YHBwdfunTp3LlzBw4cGDJkyD2vqqioKC8vv3tLQ0ODUqHQqFQ9LftfRFF0Uiq6uh5qVxqVsCoheGakz0v7CuPWn/1jfOC8mH6i+OAXwsKclAq5XGaWM8QOyOUyQRA4GnZGlImiKPJnNXJSKiQ8GqYFYXV1dXBwsPFnlUrV1NSk0WgEQQgKCnJxcVmyZMnFixdTU1PvedXOnTvff//9u7f4+/sPGzbMQ+f2CJV/TyaKzlqNVqt5+JeM1rmdCu237vil1z8v2HW++v0ZUaH67i9DhVREQRBE0UPZzSVRDkgURUEQOBp2RhQEURTN8jFoB6Q9GmJXlwl9BatWrerXr9/zzz8vCMKAAQOuXLli3P7WW2/J5fI//OEPdXV14eHhFy9edHJy+pH9vPHGG7X1DS//7vePUrpRP73X7dr6xh5dp1NS07zi4OVT1+qWxfZ78QlfpYyxobXQatRqJ6fq27VSF2IV+ri7CYLA0bAzMpnMx6tPaWWV1IVYBa1a7e7mUmamo+Hnozfp+abNEcbGxh47dkwQhMLCwsDAQEEQurq6Wlpa6urqFAqFIAhyubylpaW93Tam3wJ06k0zw1aPC9zwXUXixtwzZfRXAIDDMS0Ix40b19DQMHv27Llz565Zs0YQhOLi4oCAgJdeemnLli3jxo2LiYl5/fXXtVqbWQ7pTn/FEG/t1E15y7Mu1dNfAQCOxLQ5QplMlpGRUVlZ6eHhoVQqBUEIDg42Xghz+vTp69ev9+nTx4ZS8A4vZ+XaxEFTwj1fO3jp8Pqat8YOnDDYQ+qiAACW0JOGer1er7xv3l4URT8/P1tMwTvGBOmOLBr6dITXLzMvLMssqmpsk7oiAECvs8lbrPUejVK2YpT/3nmRV263jP7QsMlQacq1RAAA20MQdiNCr92VNOTlkb5vfXFl5pZzRdVNUlcEAOgtBGH3FDJx8TCfAwui1ArZhPScd45fa2P9CgCwRwThjxmgU22aFbY2cdDGMxUTN+b+o7Re6ooAAGZGED7YpFDP7CVDh/Vzmb75HP0VAGBnCMKH4q5SrBkfuOWZ8K+v1satN3x+oVrqigAA5kEQmuAJf9eshVHPRnu/sLs4eUdBWX336zICAGwIQWgatUL28ki/vc9F3GhoG5uas/50eScNFgBgywjCnhii136WFPHySN8/f3Vt5if5F2725JbfAABrQBD2kLG/4uDCKK1SNm5Dzuovr7Z2MDQEANtDED4Sf3fVx7PC3ksM2ppzY2ya4cQVFsoBABtDEJrBpFDPo0uiR/i7zd16fnnWpfoW+isAwGYQhOZh7K/YOifsm2t1o1MNewvprwAA20AQmtNjfm77F0QmRXu/uLsoeUdBaV2L1BUBAB6AIDSz7/sr5kfebGwfm5ZLfwUAWDmCsFeEe2s/SxryWpz/219dm/FJfmEV61cAgJUiCHuLTBSTovVHF0d7aRXjN+bSXwEA1okg7F19XZz+Ni0kZcrgHXlV8WmG4/RXAICVIQgtISFYd3hR1JggXdK28y/tLb7V1C51RQCA7xGEFuKmUrw5JmDrM+Fnyxvi0wzb86qkrggAIAgEoYXF+rnuXxj13FD9K/svJu8ouF7L+hUAIDGC0NJUcvHlkX6fL4i61dz+VOrZdSdLOzq5iAYAJEMQSiPUS7NzbsTKMQHvnSyd/HFeTkWD1BUBgIMiCCUjikJStD57UbS/u2rKx3mrv7zaQn8FAFgcQSgxvYvT/0wdnDI1ZGdeVXyq4djl21JXBACOhSC0Csb+ivhBunkZBS/tLa6mvwIALIUgtBauKsWbYwK2zQk3VNBfAQCWQxBal5/5uh5YEPn8MJ9Xsy4uyCi4Rn8FAPQygtDqKOSyFx7rfyg5urm9cwz9FQDQywhCKxXood4yO/xPCYHvf1s66aM8Q3m91BUBgH0iCK2XKAozI7yOLBoa6qWZuuncyiMljW2dUhcFAPaGILR23s7KtYmDUqaG7C+sTkg3fEl/BQCYFUFoGxKCdYeTo8cO8pifUbAss+hmI/0VAGAeBKHNcFHJ3xwTsH1OeOHNxrj1ZzYZKqWuCADsAUFoY4b7umbNj1wW23/locuzt5y7eKtZ6ooAwLYRhLbnTn+FXCaOT89Zd7K0nf4KAOgpgtBWDfRQb54Vvnpc4AenyhI/yj1bzvoVANATBKEN+76/Ijk63Es7bVPeyiMlDa0dUhcFADaGILR5Xs7KtYmDPpw+OOtCdUJ6Tjb9FQBgCoLQTsQHeXyxaOjkMM8F9FcAgCkIQvuhUcpWjPLfMy+ipKZ59IdnNxkqu7iGBgAehCC0N5F658znIn470nfVkZLZn+YXVzdJXREAWDWC0A4pZOLiYT4HFkYp5eKEDbnvHL/WRn8FAPwAgtBuBejUm2eF/WVi0IbvKhI35n5XxvoVANANgtDOTQr1PJwcPcRbO23TueVZl+rprwCA/40gtH/G/or0GSHZJbfj1hv2X7gldUUAYEUIQkfxVKDuSHL0jAivX2ZeSN5RUF7fKnVFAGAVCEIHYuyv2DsvsrKhPT41h/4KABAIQgcUodfuShry8kjft764MnPLuSL6KwA4NoLQERn7Kw4ujFIrZBPSc+ivAODICELH5e+u2jQrbG3ioI1nKiZuzP1HKf0VABwRQejoJoV6Zi8ZOqyfy/TN55ZnXapvob8CgGMhCCG4qxRrxgd++kzY11dr41IN+wqrpa4IACyHIMT3Hvd3y1oY9Wy096/2FCfvKCijvwKAYxC7pLiC/o033rjzX1ibs6W1S7edLahseHNCyIs/D5SJotQVAUAvUkj1xrX1DdfKKx99P/30Xrdr6xubmx99VzDylAnbZod8knPjtX3nN5268p/jAwd7aixcg1ajVjs5Vd+utfD7Wqc+7m6CIHA07IxMJvPx6lNaWSV1IVZBq1a7u7mUmelo+PnoTXo+X42iGzJRTIrWH1oY5ewkH7chd/WXV1s76K8AYJ8IQvwgP3fVRzND30sM2ppzY2ya4cQVRiQA7BBBiAeYFOqZvSR6hL/b3K3nX9pbXNNMfwUAu0IQ4sHcVIo14wO3zgk7W94Qn2bYU3BT6ooAwGwIQjysx/zcPl8QmRTt/Zs9xck7CkrrWqSuCADMgCCECdQK2csj/fbNj7zZ2D42LXf96fJOFrAAYOMIQpgszFv7WdKQ1+L83/7q2oxP8gurWL8CgA0jCNETxv6Ko4ujvbXK8RvprwBgwwhC9FxfF6eUaYNTpgzekVcVn2Y4Tn8FABtEEOJRJQTrDi+KGhOkS9p2/qW9xbea2qWuCABMQBDCDNxUijfHBGx9JvxsRX18mmF7HneNAmAzCEKYTayf68EFUUuG+byadXFhRuH1WtavAGADCEKYk0Iue+Gx/vvmR9W0tD2VenbdydKOTi6iAWDVCEKYX6iXZufciJVjAt47WTr547ycigapKwKAH0QQoleIopAUrc9eFD3AXTXl47yVR0oa2zqlLgoAukEQohfpXZw+mDo4ZWrI5wXV49Jzjl2+LXVFAHAvghC9zthfET9INy+j4KW9xdX0VwCwJgQhLMFVpXhzTMC2OeGGigb6KwBYFYIQlvMzX9cDCyKfH+bzatbFBRkF1+ivAGAFCEJYlLG/4lBydHN75xj6KwBYAYIQEgj0UG+ZHf6nhMD3vy1N/CjXUF4vdUUAHBdBCGmIojAzwuvIoqFhXtqpm86tPFLS0NohdVEAHBFBCCl5OyvXJg7627SQrAvV4zbkfEl/BQCLIwghvbGDdEcWDZ0U6jk/o2BZZtHNRvorAFgOQQiroFXKVozyz5gzpPBm4+gPz274R1kX19AAsAiCEFZkmK9L1vzIFx7rt2J/0S9S/37xVrPUFQGwfwQhrIuxv+L4suFymTg+Peed49faO7hJKYBeRBDCGgV6aPYsHLZ6XGD6dxWJH+edLWf9CgC9hSCElfq+vyI5OtxLO21THv0VAHoJQQir5uWsXJs4aP30kAMXqhPSc7LprwBgbgQhbMCYIN2RRUMnh3kuyChYlllU1dgmdUUA7AdBCNugUcpWjPLfMy+ipKY57kPDJkMl/RUAzIIghC2J1DtnPhfx25G+q46UzNqSX1zdJHVFAGweQQgbo5CJi4f5HFwY7aQQJ2zIfef4tTbWrwDwCAhC2KQBOtXmWWF/mRi04buKxI2535WxfgWAHiIIYcMmhXpmL4n5aT+XaZvOLc+6VE9/BQDTEYSwbTq1fM34wPQZIdklt+PWG/ZfuCV1RQBsDEEIe/BUoO5IcvSMCK9fZl5I3lFQXt8qdUUAbAZBCDth7K/YOy/yRkNbfGrO+tPlnTRYAHgIBCHsSoRe+1lSxMsjff/81bVZW/KL6K8A8CAEIezNP/srojQK2YT0HPorAPw4hdQFAL3C31318aywPQU3Xz90eV9h9X+ND/ppf5ce7OfzC9Vpp8vPVTYp5EJUX5d/+1m/nwe4mb1aABJiRAh7NinUM3vJ0OH9XadvPrc861J9i2n9Fe99U7r0swtfX6273dJ+s7H96KWaZ7flb8250UvVApAEQQg7565SrBkf+OkzYV9frR2dathXWP2QLyypaX7n+LV7NnZ1CSuPlNxuaTd3mQAkQxDCITzu75a1MCop2vvXu4uSdxSUPUR/xcHimvbuJhfrWzuOsRoUYEcIQjgKtUL28ki/vfMjqxofqr+iquEHF3u68cMPAbA5BCEcS7i3dldSxGtx/m9/dW3mJ/kXbv5gf4W3s/KHHtK7OPVOdQAkQBDC4chEMSla/8Xi6D4axbgNuau/vNra0c3QMGGQh0LezT8QF5WcC0cBe0IQwkH5uDh9OD3kvcSgrTk3xqYZTlypvecJA3Sq34/0vWejTBTfih/orqLvCLAf/HuGQ5sU6jlqoPufjl6du/X89HDPN+IH6tTyO4++8Fj/wZ6atO8qzlU2KGSySL3zL2N9HvdnOAjYFYIQjs5NpVgzPnB6hOfyrMtjUs/+YfSAmRFedx5NCPZICPaQsDwAvY2vRgFBEITH/Nz2L4x6bqj+lf0Xk3cUlNa1SF0RAAshCIHvqeTiyyP99s2PrG5qj0/LWX+6vIOblAIOgCAE/pcwb+3OZ4e8Hjfg7a+uzdySX1jF+hWAnetJELa2dn9Xjvb29i5WgIPtM/ZXHF0c7a1Vjt+Yu/rLqy3d9VcAsA8mB+HSpUvj4uJiYmLOnDlz9/Y1a9Y8+eSTI0aMyMrKMl95gGT6ujilTBv87sSgbbk3xqfn3N9fAcA+mBaE2dnZlZWVJ06cSElJeeWVV+5sz8/PP3To0IkTJ3bu3Llr1y5zFwlIZnKY5xeLox/zc5m79fzznxWW1DRLXREAMzMtCI8dOzZhwgRBEGJjY/Pz8+9s37dv35QpU3bv3v33v//9vffeM3ONgKTcVYr/HB+0d35ETXNHXGrOyiMldaw+AdgR0/oIq6urg4ODjT+rVKqmpiaNRiMIQllZ2e7du2fPnl1UVJSWlpaRkXH3q9atW7du3bq7twQEBAwfPty3r/ejFS8IgiCKooe7q4e766PvClZFFEWNWiV1Ff/i29d7fNTA3XkVv83M23W++vWxwb8eGSiXiRZ4a1EUBUGwqqMBsxBF0Swfg/ZBwqNhWhDqdLq6ujrjz62trcYUFARBq9U+88wzq1at6uzs9PX1ramp0el0d141c+bMUaNG3b2fDRs2NLW03qwxw6RLH51bQ2NTSyurAdgVlZPSSamoa7C6KzZH+GqOLf1p+umylfsL3z9+eVVCUPygPr39pq7OGkEQrPBo4FHIZKKHm6tZPgbtgMpJ6azVVJvpaHh5uJv0fNOCMDY2dvPmzc8//3xhYWFgYKAgCF1dXa2trT/5yU8yMzMFQWhoaGhra3N2dr77VXq9Xq/X373F2dm5tr6hucUMPctdXV1tbe1m2RWsh0wmymUyq/2zzh/qNTlEt/br60mf5o0Y4PbmmIDBnpreezutWiUIgtUeDfSMTCbr6uriz2okE0UJj4Zpc4Tjxo1raGiYPXv23Llz16xZIwhCcXFxQEDA5MmTb9y4MXLkyJiYmDVr1iiVP7h+DWAfPDSKN8cE7JkX0dbRNW5DLhOHgO0ybUQok8kyMjIqKys9PDyMaRccHFxeXi4Iwp49eyorK11cXLRaba9UClifSL3ztjnhB4tq3vyiZOe5qn9/wnfhT/paZuIQgLn0pKFer9d3O+bT6/WkIBxQQrDu6KKof3/C9/8cvz423fDFpRqpKwJgAm6xBpiBQi5bPMzn+PNDRw3UJe8ofHbb+Qs3ubYFsA0EIWA2xonDgwuj5aJonDisZeIQsHoEIWBmgz3VH80MTZky+HBxzciUM+tPl7ezigVgxQhCoFcYJw5fGuH3f45fT0jLOXKRiUPAShGEQG/518RhoPuinUwcAlaKIAR6FxOHgJUjCAFLYOIQsFoEIWA5TBwCVoggBCzq/onDwiomDgEpEYSABO5MHCpEcfzG3OVZl6qbmDgEpEEQApIZ7KneODN049Mh316vG/0hE4eANAhCQGJPDnQ/sCCSiUNAKgQhID0mDgEJEYSAtWDiEJAEQQhYFyYOAQsjCAFrdGfi8J3j14e/9/X+whtSVwTYLYIQsFLfTxwujUkI9nzm47NMHAK9hCAErJpOLf9zYuipF59Qypg4BHoFQQjYgBBv5w0zQjc+HXKqtG70h2fWnSxtY+IQMBOCELAZTw50z5ofuXz0gP85VZaQlnP44i2pKwLsAUEI2BKFXJYUrc9eEjM60H3JzgvPbjtfwMQh8GgIQsD26NTyN8cEHEyOUsrECRtylmddutnIxCHQQwQhYKuC+2g2zAjdOCP0VGnd6A/PMnEI9AxBCNg248Thijh/48ThoWJuVQqYhiAEbN7dE4fPf1bIxCFgEoIQsBN3Jg6dZDImDoGHRxACdiW4jyZ9RggTh8DDIwgBO3T3xOHYVAMTh8CPIAgB+3Rn4jAuSGecODx/o1HqogBrRBAC9sw4cXgoOcpJJvvFxlwmDoH7EYSA/Rv0z4nDvzNxCNyHIAQcxZMD3ffPj1wR55/y97KxqYY9BTelrgiwCgQh4ECME4dHF8dMCOnzmz3FTBwCAkEIOCCdWr5ilP/B5CiVXDROHFY1tkldFCAZghBwUIP6aNKeDjVOHMZ9aFh3srS1g4lDOCKCEHBoTw50378gyjhxmJDGxCEcEUEIODqFTEyK1mcvGWqcOJy7NZ+JQzgUghCAIAiCu0phnDhUK2RMHMKhEIQA/oWJQzggghDAve6eOBzLxCHsHUEIoBt3Jg5/8c+Jw3wmDmGnFFIXAMB6GScOn4ny/vOxaxM25E4P93z9qQFeWuU/Suu/Kqm9Wtsc4K4aHaiL6ussdaVAzxGEAB4gyEP9/pTgY5dvrzp6ZfSHhuA+6n+U1d959L++uvZvP+v32ugBElYIPAq+GgXwUJ4c6P75/MingtzvTkFBELq6hA++Lcs4VyVVYcAjIggBPCyFTPyhmcKPzlRYuBjAXAhCAA+rs6vrYnVztw9dqOp+O2D9CEIAD0smigqZ2O1Dcj5LYLM4eQGYYLivS7fba1s61p0s7eyi+x62hyAEYILfjvC7f1Cokou/eqzfu9+UPre94HptqySFAT1GEAIwQayfa8rUEB8Xpztb/NxV6TPCfv9z/4MLo9o7OsemGzYZKhkZwobQRwjANAnBuqeCYgqrmq7XtgzQqYP7qOUyURAEf3fVp88M2ZxTuepIyecF1f81IbC/q0rqYoEHY0QIwGQKmThEr00I9gj10sjv+qZUFIWkaP3BhdHN7V1j03I3GSolLBJ4SAQhADMboFNtnRP2Wpz/m0dK5m0vKKtn1hBWjSAEYH4yUUyK1h9cGNXY1hGfmsPQENZMsjlCjVrVR+f26PuRiaKLVqNWOz34qbAdCrlcLpOZ5QyxAyqlUhAEWzwafXRuB//N+92vSlYeLv7ict1fpw3p78as4fdEQRRF0Rb/rL1BIZPLpDsakgVhe3tHc7MZvjBRq1St7e2trSylbVecnJRKhcIsZ4gdkMvkgiDY7tF4Ibb/mED3X2cWDP+/J1bGBy34aT+pK7IKokxUq5xs989qXk5OSoXSbP/ktWq1Sc+XLAjb2tsbm81wTyZ3N5fW1jaz7ApWRBRkosif1UitchIEwaaPxgBX+Y654Smnypbvv7A3/8Z/jh/Y18XRv8WRyWTuLs42/Wc1L41aZa6j0UcwbWTJHCEAS1DIxBce679vflRlQ2t8ag6rVcB6EIQALCfUS5P5XMQLj/X7/ecXl2UW3Wxsl7oigCAEYFn/HBpGltQ0j00z7CuslroiODqCEIAEwry1mUlDnh/u8+vdRcsyi6qbGBpCMgQhAGko5LIXHuu/d37kpVtN8WmG/RduSV0RHBRBCEBK4d7aPc9FPD/MZ9nuomWZRbcYGsLiCEIAEjMODXfOHXK+qiE+zXCAoSEsiyAEYBVi+jnvXxA9K9J7aeaFZZlFNc0dUlcER0EQArAWKrm4YpT/jrlDzt9oHJN69mARQ0NYAkEIwLr8tL/L/oVRsyK9n991YVlm0e0WZg3RuwhCAFbHODTMmDMk/0bjmNScQ8U1UlcEe0YQArBSw3xd9i+InBnhteSzwuVZl+pbmTVEryAIAVgvtUK2YpT/9jnhJ67WJqQZviqplboi2CGCEIC1G+7remBh1JRwr+e2n1+edamBoSHMiiAEYAOMQ8Ntc8KPX60dm55z4gpDQ5gNQQjAZvzM1/XAgqgpYZ7PbmNoCLMhCAHYEo1StmKU/0czQo9erklIz/n6ap3UFcHmEYQAbM+TA90PLYweFeA+59P85VmXGts6pa4INowgBGCTXFTyNeMDN84M/eJSzbj0nJPXmDVEDxGEAGzY6IHuh5Kjfj7AbfaW88uzLjUxNITpCEIAts1VpVgzPnDDjJAjF2sSP8o9U9YgdUWwMQQhAHsQF6g7vCgq1td1+ua81V9ebe3okroi2AyCEICdMA4N108P2ZFX9YuNOYbyeqkrgm0gCAHYlTFBusOLon7W33XqpnMMDfEwCEIA9sZNpVgzPvDD6YMz8qombszJqWDWED+GIARgn+KDPA4viorUO0/5OG/1l1fbOhkaonsEIQC75a5SrE0c9N6kQZ/mVE7cmMvQEN0iCAHYuUmhnoeSo4M81FM2nWNoiPsRhADsn5dW+T9TB7+XGLTFUJm4MTevslHqimBFCEIAjmJSqOfh5OgAnWrSx3mrv7za3sFtaCAIBCEAh+LlrPzbtJD3EoM+MVRO/+TchZvNUlcE6RGEAByOcWiod3H6xQbDupOlHcwaOjaCEIAj8nZWrp8WsjZx0Pvflj79ybmi6iapK4JkCEIAjmtSqOeRRUM9tcoJ6TkMDR0WQQjAoXk7K1Onh6xNHLTuZNmMLeeKGRo6HoIQAIxDwygPtWLChtx1J0s7uxgaOhCCEAAEQRD0Lk6p00PfiA9475vSmZ/kX7rFBaWOgiAEgO+JopAUrT+YHK1SiOM35DA0dBAEIQD8L35uTptnha8cE/DuN6UzP8m/zNDQ3hGEAHAv49Dw0MIopVwctyFn/elyRoZ2jCAEgO75uau2zA5fOSbgz19dS9p+/nptq9QVoVcQhADwg76fNVwY1d7ROTbdsMlQydDQ/hCEAPAA/u6qT58Z8nrcgFVHSuZtP19a1yJ1RTAnghAAHsw4NDywML4HbDwAAA+XSURBVKq5vWtsWu4mQ6XUFcFsCEIAeFgBOvXWOWGvxfm/eaRk3vaCsnpmDe0BQQgAJpCJYlK0fu+8yJrmtnFphpRTZS0dTBvaNoIQAEzmJBd1amVtS+dbR68M/su3kz7Ky61skLoo9BBBCACmKa9vnb753NFLNcb7znR1CWfL66d8nJdbQRbaJIIQAEzz7telNxra7tnY1tE1e0v+mTKy0PYQhABgmmMlt7vd3tjWMXVT3vKsS/WtHRYuCY+CIAQA09S3dJ9zHV3C36aHZJfcjltv+PxCtYWrQo8RhABgGj93p263612cxg3SfbEo+tlo719lFiXvoL/CNhCEAGCaGUO8fmS7WiF7eaTf3vmRNxvb41Nz1p8uZy0nK0cQAoBpnovpOzGkzz0bY/1cXxrhe+fXcG/tZ0lDXovzf/urazM+yS+sarJsjTCBogevaW1tdXLq/puBxsZGtVotk5GvAOyWQiZ+MGXwnoLqw8W3Sm639HNVPhmgmxXppZCJdz/N2Ho/Nkj3x8Ml4zfmLh3u8/JIP5Vc/KHdQiomB+HSpUtzc3MbGxvT09NjYmLufuj27dsxMTEHDx4MDg42X4UAYHVEUZgc1mdy2L3jwvv1dXH6n6mDDxbVvH7o0t6C6tUJA58c6G6BCvHwTBu6ZWdnV1ZWnjhxIiUl5ZVXXrnn0V/96letrcwMA8C9EoJ1h5Kj4gfp5mUULMssutnYLnVF+BfTRoTHjh2bMGGCIAixsbH5+fl3P5SWlhYaGlpeXt7tq44dO3b3litXrvTr19/DzdX0gu8lE0VnrVrlpHz0XcF6KBRyuUxuljPEDjgplYIgcDRsnYcgvDvN49mf+v1mV/5TqWdXjR/8gpcnf1YjhUIuE0WpjoZpQVhdXX3na0+VStXU1KTRaARBKCws3Lp16549e7Kzs+9/VWtra319/d1bOjo6BFEQZeb5rlwURHPtClZCFEXRfGeIrRNFQRA4GnbiiYEeJ379eMrJq8v3FW7Pqfy/k0MHeztLXZT0REEUpDvJTQtCnU5XV1dn/Lm1tdWYgoIgLF++3N/ff/Xq1ZcvX163bt2rr77at2/fO6+Kj4+Pj4+/ez9vvPFGbX1DdU3toxUvCILQT+9V39jU2Nz86LuC9dBq1Gonp+rbZjhD7EAfdzdBEDga9mTuEI+f+0X/8cjVx//6zbLYfi8+4at07P+jo1Wr3d1czBIKgiBofdQmPd+0OcLY2Fjjl5yFhYWBgYGCIHR1dbW0tPzud7+bOHFiZGSks7NzSEiIWm1aEQDgaAJ06qylj/1lYtDGMxUTN+b+o7T+wa9B7zAtCMeNG9fQ0DB79uy5c+euWbNGEITi4uKAgIARI0ZMmzZt2rRp3t7eY8eOdXfnmigAeLBJoZ7ZS4YO6+cyffO55VmXfujmbehVpn01KpPJMjIyKisrPTw8lEqlIAjBwcF3XyBz6NAhMxcIAHbNXaVYMz5w2hDPFQcuj041rBjlPzOi+zvXoJf0pPNdr9cbUxAAYBaP+7t9viAyKdr7lf0Xk3cUXK+lFc1yuAUMAFgF401K982PvNXcPib17LqTpR2d3KTUEghCALAiYd7anXMj/v+EwPe/LZ30UZ6hnItoeh1BCADWRRSFmRFeh5OjQ700UzedW3mkpIGVfnsTQQgA1kjv4rQ2cdDfpoVkXahOSM/Jvnxb6orsFkEIANZr7CDdF4uGTg7zXJBRsCyzqKqxTeqK7BBBCABWTaOUrRjlv2dexJXbLaM/NGwyVLLQr3kRhABgAyL1zruShrw80vetL67M3HLuwk3uK2k2BCEA2AaFTFw8zOfgwiitUj5uQ87qL6+2djA2NAOCEABsib+76qOZoe8lBm3NuTE2zXD8Cndjf1QEIQDYnkmhntlLop8K0iVtO//S3uJbTaz023MEIQDYJDeV4s0xAVufCT9bUR+fZtieVyV1RbaKIAQAGxbr53pwQdSSYT6vZl1ckFFwjZuUmo4gBADbppDLXnis/6Hk6Jb2Tm5S2gMEIQDYg0AP9Sezw/+UEPjBqbKJH+WeKWuQuiKbQRACgJ0w3qT0SHJ0uJd26qa85VmX6rlJ6UMgCAHArng5K9cmDkp7OiS75PZT6w2fX6iWuiJrRxACgB0aE6Q7khw9N9r7hd3FyTsKyuq5iOYHEYQAYJ80StnLI/32PhdR1dgWn5qz/nR5J3cp7Q5BCAD2bIheuysp4rU4/7e/ujbzk/zCqiapK7I6BCEA2DmZKCZF648ujvbSKsdvzF395dUWblJ6F4IQABxCXxenlGmDU6YM3plXFZ9qOMZKv/9EEAKAA0kI1h1eFBU/SDcvo+ClvcXV3KSUIAQAR+OqUrw5JmD7nPCcyobRH55hpV+CEAAc0XBf16z5kcti+688dPmZT88VVzvuRTQEIQA4KONNSg8vilbIZRM25L5z/FqbQ96klCAEAIcWoFNvnhX2l4lBG89UJG7M/UdpvdQVWRpBCAAQJoV6Hl0c89N+LtM3n1uedam+xYFuUkoQAgAEQRB0avma8YGfPhP2zbW60amGvYWOcpNSghAA8C+P+7vtXxCZFO394u6i5B0FpXUtUlfU6whCAMD/olbIXh7pt29+ZHVT+1PrDXa/0i9BCADoRpi39rNnI/44JuC9k6WTP87LqbDblX4JQgBA90RRSIrWZy+KDvHUTPk4b+WRksa2TqmLMj+CEADwY/QuTmsTB6VMDcm6UJ2Qbsi2u5uUEoQAgAdLCNYdWTR0UqjngoyCZZlFVY1tUldkNgQhAOChaJWyFaP898yLuHK7Je5Dg93cpJQgBACYIFLvvCtpyG9H+r71xZVZW/KLbP8mpQQhAMA0Cpm4eJjPwYVRaoU4IT3nnePXWm15pV+CEADQE/7uqo9nha1NHPTRmcqxaYYTV2qlrqiHCEIAQM9NCvXMXhI9wt9t7tbzL+0tvmWDK/0ShACAR+KmUqwZH7htTvjZivr4NMP2vCqpKzINQQgAMINYP9eDC6KWDPN5NeviwozCa7WtUlf0sAhCAIB5GFf6PZQc3dzeMSb1rK3cpJQgBACYU6CH+pPZ4X9KCPzgVFniR7lny639JqUEIQDAzERRmBnhdTg5OsxLO21T3sojJQ2t1rvSL0EIAOgV3s7KtYmD1k8POXChOm69Yf+FW1JX1D2CEADQi8YE6Y4sGvp0hNcvMy8k7ygor7e6i2gIQgBA79IoZStG+e+bF1nV2BafmrP+dHmnNd2llCAEAFjCEL12V1LEH+L83/7q2sxP8i/ctJablBKEAAALkYliUrT+i8XRXlrluA25q7+82mIFNyklCAEAFuXj4pQybXDKlME786rGphmOSb3SL0EIAJBAQrDu8KKoMUG6eRkFyz47X9Ug2UU0BCEAQBquKsWbYwK2zQk/U1b3RdFNqcogCAEAUhrq47Juapins1OZRJ0VCkneVRAEhUKhVqkefT+iKCqVCnWXGXYF66FUKuRymVnOEDsgl8sEQeBo2BmZTBRFkT9r6t9LVx+9XNP8/eJNTwxwf3dyyEAPjSVrkCwInZQKF6360fcjCoLayUmpkD/6rmA95HK5TJSZ5QyxAwq5QhAEjobdEQWH/7OuP3Xt1f1Fd2/5+srtpzflnPrNE+5qpcXKkCwIG5uaq26Z4UqhfnqvuobGxubmR98VrIdWo1Y7OVXfttUFr82rj7ubIAgcDTsjk8l8vPqY5WPQRnV0dr1x4ML926/WNK89WvSrx/r3eM9+PnqTns8cIQBAAhdvtVQ1tHX70Klr9ZashCAEAEigpeMH16NobrfoUhUEIQBAAgHuKoW8+wwa5GnRi2UIQgCABFxViimhfe7fLhPFudHelqyEIAQASOM/xgb8zNf17i0KuWz1uIGRemdLliHZVaMAAAfnqlJsmxOeVXTLUNl8u6XT11lMDPUM9LB0SwlBCACQjFwmTgzpMzNa7e7mUlZZJUkNfDUKAHBoBCEAwKERhAAAh0YQAgAcGkEIAHBoBCEAwKERhAAAh0YQAgAcGkEIAHBo0txZ5vLlywWFF955+8+PvitXZ21zS2tbe/uj7wrWQ6lUKOTypuYWqQuxChq1ShAEjoadEUXRRaupa2iUuhCroFQo1Concx2NKZMS4+LiHv750gRhTEyMIAhuLma4rep33303aNAgNze3R98VrI3ShVsACoIgXL58WRCEgQMHSlwHeoFZPgbtQG1t7fn8cz/5yU8keXexq6tLkjc2lxEjRrz99tsjRoyQuhCgt7zxxht3/gvYpRMnTvzud787ceKEJO/OHCEAwKERhAAAh0YQAgAcGkEIAHBoNn9V3uLFi/39/aWuAuhFJl0IDtgif3//xYsXS/XuNn/VKAAAj4KvRgEADo0gBAA4NGsPws7Ozqampru3tLa2/sjz29ra7vmyNysra9euXb1SHGAOpp7k95/S169f37Jly5YtW2pqanqlRMBMfvxs7/bj2gKf4dYehJs3b/7v//7vO78uXbo0Li4uJibmzJkz9z/56tWrTz/9dEVFxZ0t27dvr6io0Gg0KSkpligXMJ1JJ3m3p/SpU6daWloef/xxrVZriYqBnvqRs73bc9syn+FWHYSTJ0/+5S9/eefX7OzsysrKEydOpKSkvPLKKzdu3Dj6TwUFBYIg+Pj43HOvtW+++WbChAkJCQnffvutpasHHoKpJ3m3p7Sfn19JScmqVauMzwGs04+f7d2e25b5DLfq9ondu3f/9a9/bW5uNv567NixCRMmCIIQGxubn58vl8vVarXxIaVSafyvRqO5ew96vb6iokIURS8vL8vWDjwUU0/ybk/p4cOHDx8+PDc3NzMzMyoqyuL/I4CH8uNn+9ixY+8/ty3zGW7VQXiP6urq4OBg488qlUqj0Tz++OP3PCc0NNSYhTqd7vz584sWLVq7dq0oii+++KKlywVM98CT/J5T2nie5+XlnTx5sqKi4ve//70ERQM9cs/Z/uyzz37wwQf3nNuW+Qy3pSDU6XR1dXXGn1tbW+8Z/Bn94he/MP5w56qB//iP/7BMecCje+BJ7uXldfcpbTzPfXx84uPjLVYkYBb3nO1+fn73n9uCRT7DrXqO8B6xsbHHjh0TBKGwsDAwMFDqcgDz4ySH47Ces92WgnDcuHENDQ2zZ8+eO3fumjVrpC4HMD9OcjgO6znbbe8Wa5WVlR4eHsYLBwC7xEkOx2ENZ7vtBSEAAGZkS1+NAgBgdgQhAMChEYQAAIdGEAIAHBpBCABwaAQhAMChEYQAAIdGEAIAHBpBCABwaAQhAMChEYQAAIdGEAIAHBpBCABwaAQhAMChEYQAAIdGEAIAHNr/A8CGwYT2xbeFAAAAAElFTkSuQmCC",
      "text/plain": [
       "Figure()"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig = Figure(resolution=(600, 400))\n",
    "ax = Axis(fig[1, 1], xscale=log10)\n",
    "scatter!(ax, a_ex_range[2:end], y[2:end] ./ _λ0)\n",
    "\n",
    "lines!(ax, a_ex_range[2:end], y_ppc)\n",
    "\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "id": "ff4b3603",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Vector{Float64}:\n",
       " 0.0\n",
       " 0.1\n",
       " 0.26\n",
       " 0.4\n",
       " 0.6\n",
       " 1.0"
      ]
     },
     "execution_count": 144,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a_ex_range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "id": "17e55124",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Vector{Float64}:\n",
       " 1.0000000000000002\n",
       " 0.8123683747624234\n",
       " 0.6218979659657092\n",
       " 0.5149523745253777\n",
       " 0.4128075179347371\n",
       " 0.2950427888416317"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y ./ _λ0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "id": "e4bec8c9",
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "ArgumentError: \"20220714_r2_allstrains_tc/exp_analysis_summary.csv\" is not a valid file or doesn't exist",
     "output_type": "error",
     "traceback": [
      "ArgumentError: \"20220714_r2_allstrains_tc/exp_analysis_summary.csv\" is not a valid file or doesn't exist",
      "",
      "Stacktrace:",
      " [1] CSV.Context(source::CSV.Arg, header::CSV.Arg, normalizenames::CSV.Arg, datarow::CSV.Arg, skipto::CSV.Arg, footerskip::CSV.Arg, transpose::CSV.Arg, comment::CSV.Arg, ignoreemptyrows::CSV.Arg, ignoreemptylines::CSV.Arg, select::CSV.Arg, drop::CSV.Arg, limit::CSV.Arg, buffer_in_memory::CSV.Arg, threaded::CSV.Arg, ntasks::CSV.Arg, tasks::CSV.Arg, rows_to_check::CSV.Arg, lines_to_check::CSV.Arg, missingstrings::CSV.Arg, missingstring::CSV.Arg, delim::CSV.Arg, ignorerepeated::CSV.Arg, quoted::CSV.Arg, quotechar::CSV.Arg, openquotechar::CSV.Arg, closequotechar::CSV.Arg, escapechar::CSV.Arg, dateformat::CSV.Arg, dateformats::CSV.Arg, decimal::CSV.Arg, truestrings::CSV.Arg, falsestrings::CSV.Arg, stripwhitespace::CSV.Arg, type::CSV.Arg, types::CSV.Arg, typemap::CSV.Arg, pool::CSV.Arg, downcast::CSV.Arg, lazystrings::CSV.Arg, stringtype::CSV.Arg, strict::CSV.Arg, silencewarnings::CSV.Arg, maxwarnings::CSV.Arg, debug::CSV.Arg, parsingdebug::CSV.Arg, validate::CSV.Arg, streaming::CSV.Arg)",
      "   @ CSV ~/.julia/packages/CSV/jFiCn/src/context.jl:236",
      " [2] #File#25",
      "   @ ~/.julia/packages/CSV/jFiCn/src/file.jl:221 [inlined]",
      " [3] File",
      "   @ ~/.julia/packages/CSV/jFiCn/src/file.jl:162 [inlined]",
      " [4] read(source::String, sink::Type; copycols::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})",
      "   @ CSV ~/.julia/packages/CSV/jFiCn/src/CSV.jl:91",
      " [5] read(source::String, sink::Type)",
      "   @ CSV ~/.julia/packages/CSV/jFiCn/src/CSV.jl:87",
      " [6] top-level scope",
      "   @ In[171]:12",
      " [7] eval",
      "   @ ./boot.jl:368 [inlined]",
      " [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
      "   @ Base ./loading.jl:1281"
     ]
    },
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mThe Kernel crashed while executing code in the the current cell or a previous cell. Please review the code in the cell(s) to identify a possible cause of the failure. Click <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info. View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "files = [\n",
    "    \"20220714_r2_allstrains_tc\",\n",
    "    \"20220715_r3_allstrains_tc\",\n",
    "    \"20220722_r1_allstrains_tc\",\n",
    "    \"20220723_r2_allstrains_tc\",\n",
    "    \"20220725_r3_allstrains_tc\"\n",
    "    ]\n",
    "\n",
    "df = DataFrame()\n",
    "\n",
    "for file in files\n",
    "    _df = CSV.read(file * \"/exp_analysis_summary.csv\", DataFrame)\n",
    "    insertcols!(_df, 5, :tc => map(x -> parse(Float64, split(x, \"_\")[1]), _df.pos_selection))\n",
    "    insertcols!(_df, 5, :rep => map(x -> parse(Int64, split(x, \"_\")[2]), _df.strain))\n",
    "    insertcols!(_df, 5, :date => parse(Int64, split(file, \"_\")[1]))\n",
    "    _df.strain = map(x -> split(x, \"_\")[1], _df.strain)\n",
    "    df = vcat(df, _df)\n",
    "end\n",
    "first(df, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a676643",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0-rc1",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}