diff --git a/notebooks/18_visualise_param_scan.ipynb b/notebooks/18_visualise_param_scan.ipynb index 8902c94e..ae1edc41 100644 --- a/notebooks/18_visualise_param_scan.ipynb +++ b/notebooks/18_visualise_param_scan.ipynb @@ -4739,9 +4739,13 @@ "# np.power(df['Sensitivity'], -2) + 100 * np.power(df['Precision'], -2) + 0.01 * np.power(df['Fold change'], -2)\n", "# )\n", "\n", - "df['Adaptation'] = df['Sensitivity'] * np.where(df['Precision'] > 1e2, 1e2, df['Precision']) \n", - "df['Adaptation'] = np.where(df['Adaptation'] == np.inf, 0, df['Adaptation'])\n", - "df['Adaptation'] = np.where(df['Adaptation'] == -np.inf, 0, df['Adaptation'])\n", + "def calc_adaptation(s, p):\n", + " adaptation = s * np.where(p > 1e2, 1e2, p)\n", + " adaptation = np.where(adaptation == np.inf, 0, adaptation)\n", + " adaptation = np.where(adaptation == -np.inf, 0, adaptation)\n", + " return adaptation\n", + "\n", + "df['Adaptation'] = calc_adaptation(df['Sensitivity'], df['Precision'])\n", "\n", "fig = plt.figure(figsize=(8, 7))\n", "ax = plt.subplot(1,1,1)\n", diff --git a/notebooks/22_adaptation_autograd.ipynb b/notebooks/22_adaptation_autograd.ipynb new file mode 100644 index 00000000..2f670974 --- /dev/null +++ b/notebooks/22_adaptation_autograd.ipynb @@ -0,0 +1,343 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Optimising for adaptation as a circuit function through automatic gradient calculation\n", + "\n", + "In order to find the derivative of adaptation, we need to be able to differentiate the dynamic simulation and the analytics calculated from it that are used to calculate adaptability. We need to be able to find the derivative of the sensitivity and precision with respect to the circuit topology. \n", + "\n", + "1. Set up a simple test circuit simulation environment\n", + "2. Re-write the dynamic simulation to keep track of the max / min of all species\n", + "3. Try to differentiate that" + ] + }, + { + "cell_type": "code", + "execution_count": 85, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The autoreload extension is already loaded. To reload it, use:\n", + " %reload_ext autoreload\n" + ] + } + ], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Imports " + ] + }, + { + "cell_type": "code", + "execution_count": 86, + "metadata": {}, + "outputs": [], + "source": [ + "from jax import jacrev\n", + "import numpy as np\n", + "import jax\n", + "import jax.numpy as jnp\n", + "import pandas as pd\n", + "import seaborn as sns\n", + "import matplotlib.pyplot as plt\n", + "import diffrax as dfx\n", + "\n", + "from functools import partial\n", + "import os\n", + "import sys\n", + "\n", + "jax.config.update('jax_platform_name', 'cpu')\n", + "\n", + "if __package__ is None:\n", + "\n", + " module_path = os.path.abspath(os.path.join('..'))\n", + " sys.path.append(module_path)\n", + "\n", + " __package__ = os.path.basename(module_path)" + ] + }, + { + "cell_type": "code", + "execution_count": 87, + "metadata": {}, + "outputs": [], + "source": [ + "from synbio_morpher.srv.parameter_prediction.simulator import make_piecewise_stepcontrol\n", + "from synbio_morpher.utils.misc.type_handling import flatten_listlike\n", + "from synbio_morpher.utils.results.analytics.timeseries import calculate_adaptation, compute_sensitivity, compute_precision" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Test environment for example circuits" + ] + }, + { + "cell_type": "code", + "execution_count": 88, + "metadata": {}, + "outputs": [], + "source": [ + "def make_species_bound(species_unbound):\n", + " return sorted(set(flatten_listlike([['-'.join(sorted([x, y])) for x in species_unbound] for y in species_unbound])))\n", + "\n", + "\n", + "# RNA circuit settings\n", + "species_unbound = ['RNA_0', 'RNA_1', 'RNA_2']\n", + "species_bound = make_species_bound(species_unbound)\n", + "species = species_unbound + species_bound\n", + "species_signal = ['RNA_0']\n", + "species_output = ['RNA_1']\n", + "idxs_signal = [species.index(s) for s in species_signal]\n", + "idxs_output = [species.index(s) for s in species_output]\n", + "signal_onehot = np.array([1 if s in idxs_signal else 0 for s in np.arange(len(species))])\n", + "\n", + "# Initial parameters\n", + "n_circuits = 3\n", + "k = 0.00150958097\n", + "N0 = 200\n", + "y00 = np.array([[N0, N0, N0, 0, 0, 0, 0, 0, 0]]).astype(np.float32)\n", + "y00 = np.repeat(y00, repeats=n_circuits, axis=0)\n", + "\n", + "# Simulation parameters\n", + "signal_target = 2\n", + "t0 = 0\n", + "t1 = 100\n", + "ts = np.linspace(t0, t1, 500)\n", + "dt0 = 0.0005555558569638981\n", + "dt1_factor = 5\n", + "dt1 = dt0 * dt1_factor\n", + "max_steps = 16**4 * 10\n", + "\n", + "# Reactions\n", + "rates = np.array([[1e-4, 1e-4, 1e1],\n", + " [1e-4, 1e-6, 1e-4],\n", + " [1e1, 1e-4, 1e-4]])\n", + "rates = np.random.randint(-6, 2, size=(n_circuits, len(species_unbound), len(species_unbound)))\n", + "rates = np.exp(rates)\n", + "\n", + "inputs = np.array([\n", + " [2, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 1, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 0, 1, 0, 0, 0, 0, 0, 0],\n", + " [0, 2, 0, 0, 0, 0, 0, 0, 0],\n", + " [0, 1, 1, 0, 0, 0, 0, 0, 0],\n", + " [0, 0, 2, 0, 0, 0, 0, 0, 0],\n", + "])\n", + "outputs = np.array([\n", + " [0, 0, 0, 1, 0, 0, 0, 0, 0],\n", + " [0, 0, 0, 0, 1, 0, 0, 0, 0],\n", + " [0, 0, 0, 0, 0, 1, 0, 0, 0],\n", + " [0, 0, 0, 0, 0, 0, 1, 0, 0],\n", + " [0, 0, 0, 0, 0, 0, 0, 1, 0],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1],\n", + "])\n", + "\n", + "# Rates\n", + "reverse_rates = np.array(list(map(lambda r: r[np.triu_indices(len(species_unbound))], rates)))\n", + "forward_rates = np.ones_like(reverse_rates) * k" + ] + }, + { + "cell_type": "code", + "execution_count": 89, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 0.98, 'Jacobian of system')" + ] + }, + "execution_count": 89, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8gAAAF2CAYAAACyDbEuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1lklEQVR4nO3deXxTdb7/8Xfa0rTSBQplKZQWcGEpqLeIQhlQQZQBUWYGRUEWRRgsFoY7KOBVcYGq4yCIyOI4wNWyCCoqIFKRRTbZREQdWQSsbBWFLiCptOf3hz9yqU2gaZOcnuT1fDzyR05P8v2c0MebvHuy2AzDMAQAAAAAQJALMXsAAAAAAACqAgoyAAAAAACiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSKMgAAAAAAEiiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSKMgAAItbs2aNbDabFi9efMl9Bw4cqOTkZN8PVQmFhYUaPHiw6tWrJ5vNppEjR5o9EgAAQYOCDACosDlz5shms2nbtm1mjxIwJk6cqDlz5mjYsGF64403dN9995k9kiRp+fLlGj9+vNljAADgU2FmDwAAgL+89tprKikpMXuMi/rkk090ww036MknnzR7lFKWL1+uadOmUZIBAAGNM8gAgKBRrVo12e12s8e4qNzcXNWoUcPsMQAACEoUZACA1+zatUsDBw5UkyZNFBERoXr16un+++/XTz/9VGbfw4cP64EHHlBCQoLsdrsaN26sYcOGqaioyLnPd999p969eysuLk6XXXaZbrjhBi1btszl2sXFxRo3bpzq1aun6tWrq2fPnsrJySm1j6v3IL/44otq3769atWqpcjISKWmprp8P7PNZtPw4cO1ZMkSpaSkyG63q2XLllqxYkW5Hpvc3Fw98MADqlu3riIiInT11Vdr7ty5zp+ffy/1gQMHtGzZMtlsNtlsNh08eNDtfWZnZ6tDhw6qUaOGoqKidNVVV2ncuHGSfnsvc/Xq1TVixIgyt/vhhx8UGhqqzMxMSdKvv/6qp556SldccYUiIiJUq1YtdejQQdnZ2c7Hbdq0ac7H4fzlvJKSEk2ePFktW7ZURESE6tatq6FDh+rkyZOl1k1OTlaPHj20Zs0atWnTRpGRkWrVqpXWrFkjSXrnnXfUqlUrRUREKDU1VZ9//nm5HlsAALyFl1gDALwmOztb3333nQYNGqR69erpq6++0qxZs/TVV19p8+bNzlJ15MgRtW3bVqdOndKQIUPUrFkzHT58WIsXL9aZM2cUHh6u48ePq3379jpz5owyMjJUq1YtzZ07Vz179tTixYvVq1evUmtPmDBBNptNjz76qHJzczV58mR16dJFO3fuVGRkpNuZp0yZop49e6pv374qKirSggUL1Lt3by1dulTdu3cvte/69ev1zjvv6KGHHlJ0dLRefvll/fnPf9b333+vWrVquV3jl19+0Y033qh9+/Zp+PDhaty4sRYtWqSBAwfq1KlTGjFihJo3b6433nhDf/vb39SwYUP993//tyQpPj7e5X1+9dVX6tGjh1q3bq2nn35adrtd+/bt04YNGyRJUVFR6tWrlxYuXKhJkyYpNDTUedv58+fLMAz17dtXkjR+/HhlZmZq8ODBatu2rfLz87Vt2zbt2LFDt9xyi4YOHaojR44oOztbb7zxRplZhg4dqjlz5mjQoEHKyMjQgQMH9Morr+jzzz/Xhg0bVK1aNee++/bt07333quhQ4eqX79+evHFF3X77bdrxowZGjdunB566CFJUmZmpu666y59++23Cgnh7/kAAD8xAACooNmzZxuSjK1btxqGYRhnzpwps8/8+fMNSca6deuc2/r372+EhIQ4b3ehkpISwzAMY+TIkYYk49NPP3X+rKCgwGjcuLGRnJxsFBcXG4ZhGKtXrzYkGQ0aNDDy8/Od+7711luGJGPKlCnObQMGDDCSkpJKrff7mYuKioyUlBTj5ptvLrVdkhEeHm7s27fPue2LL74wJBlTp051/QD9f5MnTzYkGW+++Wapddq1a2dERUWVmjspKcno3r37Re/PMAzjpZdeMiQZP/74o9t9PvroI0OS8eGHH5ba3rp1a6NTp07O61dfffUl10xPTzdcPW349NNPDUlGVlZWqe0rVqwosz0pKcmQZGzcuLHMjJGRkcahQ4ec22fOnGlIMlavXn3RuQAA8Cb+JAsA8JoLz9SePXtWJ06c0A033CBJ2rFjh6TfXo67ZMkS3X777WrTpk2Z+zh/lnn58uVq27atOnTo4PxZVFSUhgwZooMHD+rrr78udbv+/fsrOjraef0vf/mL6tevr+XLl5d75pMnTyovL09/+MMfnPNeqEuXLmratKnzeuvWrRUTE6PvvvvuomssX75c9erV0z333OPcVq1aNWVkZKiwsFBr16696O1dOf8+5ffee8/tB4916dJFCQkJysrKcm7bvXu3du3apX79+pW6r6+++kp79+71eI5FixYpNjZWt9xyi06cOOG8pKamKioqSqtXry61f4sWLdSuXTvn9euvv16SdPPNN6tRo0Zltl/qsQUAwJsoyAAAr/n55581YsQI1a1bV5GRkYqPj1fjxo0lSXl5eZKkH3/8Ufn5+UpJSbnofR06dEhXXXVVme3Nmzd3/vxCV1xxRanrNptNl19++UXfwytJS5cu1Q033KCIiAjFxcUpPj5e06dPd857oQsL3Hk1a9Ys815bV8dyxRVXlHmpsLtjKY+7775baWlpGjx4sOrWras+ffrorbfeKlWWQ0JC1LdvXy1ZskRnzpyRJGVlZSkiIkK9e/d27vf000/r1KlTuvLKK9WqVSuNHj1au3btKtcce/fuVV5enurUqaP4+PhSl8LCQuXm5pba//ePYWxsrCQpMTHR5fZLPbYAAHgT70EGAHjNXXfdpY0bN2r06NG65pprFBUVpZKSEt12221V8uuVPv30U/Xs2VMdO3bUq6++qvr166tatWqaPXu25s2bV2b/C9/HeyHDMHw9ahmRkZFat26dVq9erWXLlmnFihVauHChbr75Zq1cudI5a//+/fWPf/xDS5Ys0T333KN58+apR48ezgIqSR07dtT+/fv13nvvaeXKlfrXv/6ll156STNmzNDgwYMvOkdJSYnq1KlT6iz1hX7/Hmp3j2FVemwBAMGLggwA8IqTJ09q1apVeuqpp/TEE084t//+Zbvx8fGKiYnR7t27L3p/SUlJ+vbbb8ts/89//uP8+YV+v45hGNq3b59at27tdo23335bERER+uijj0p9/dPs2bMvOpunkpKStGvXLpWUlJQ6i+zuWMorJCREnTt3VufOnTVp0iRNnDhRjz32mFavXq0uXbpIklJSUnTttdcqKytLDRs21Pfff6+pU6eWua+4uDgNGjRIgwYNUmFhoTp27Kjx48c7C/KFn1p9oaZNm+rjjz9WWlraRT8MDQAAK+Al1gAArzh/BvD3Z/wmT55c6npISIjuvPNOffDBB9q2bVuZ+zl/+z/+8Y/asmWLNm3a5PzZ6dOnNWvWLCUnJ6tFixalbve///u/KigocF5fvHixjh49qm7dul10ZpvNpuLiYue2gwcPasmSJRc/WA/98Y9/1LFjx7Rw4ULntnPnzmnq1KmKiopSp06dPL7Pn3/+ucy2a665RpLkcDhKbb/vvvu0cuVKTZ48WbVq1SrzmPz+a7iioqJ0+eWXl7qf6tWrS5JOnTpVat+77rpLxcXFeuaZZ8rMc+7cuTL7AwBQlXEGGQDgFTExMerYsaNeeOEF/frrr2rQoIFWrlypAwcOlNl34sSJWrlypTp16qQhQ4aoefPmOnr0qBYtWqT169erRo0aGjNmjObPn69u3bopIyNDcXFxmjt3rg4cOKC33367zPt54+Li1KFDBw0aNEjHjx/X5MmTdfnll+vBBx90O3P37t01adIk3Xbbbbr33nuVm5uradOm6fLLLy/3e3DLY8iQIZo5c6YGDhyo7du3Kzk5WYsXL9aGDRs0efLkUh8uVl5PP/201q1bp+7duyspKUm5ubl69dVX1bBhw1IfbCZJ9957rx555BG9++67GjZsWKmvXZJ+++CsG2+8UampqYqLi9O2bdu0ePFiDR8+3LlPamqqJCkjI0O33nqrQkND1adPH3Xq1ElDhw5VZmamdu7cqa5du6patWrau3evFi1apClTpugvf/lLBR41AAD8j4IMAKiw82d7z589njdvnh5++GFNmzZNhmGoa9eu+vDDD5WQkFDqdg0aNNBnn32mxx9/XFlZWcrPz1eDBg3UrVs3XXbZZZKkunXrauPGjXr00Uc1depUnT17Vq1bt9YHH3xQ5vuJJWncuHHatWuXMjMzVVBQoM6dO+vVV1913p8rN998s15//XU999xzGjlypBo3bqznn39eBw8e9GpBjoyM1Jo1azRmzBjNnTtX+fn5uuqqqzR79mwNHDiwQvfZs2dPHTx4UP/+97914sQJ1a5dW506ddJTTz1V6v3F0m+PZdeuXbV8+XLdd999Ze4rIyND77//vlauXCmHw6GkpCQ9++yzGj16tHOfP/3pT3r44Ye1YMECvfnmmzIMQ3369JEkzZgxQ6mpqZo5c6bGjRunsLAwJScnq1+/fkpLS6vQ8QEAYAabwadfAAAq6OWXX9aIESO0b9++Ul9/hKqnV69e+vLLL7Vv3z6zRwEAoMriPcgAgArbunWrqlevXuEPmYJ/HD16VMuWLXN59hgAAPwfXmINAPDY22+/rTVr1igrK0uDBw9WWBj/nVRFBw4c0IYNG/Svf/1L1apV09ChQ80eCQCAKo1nNAAAj/39739XQUGBHnjgAb300ktmjwM31q5dq0GDBqlRo0aaO3eu6tWrZ/ZIAABUabwHGQAAAAAA8R5kAAAAAAAkUZABAAAAAJBEQQYAAAAAQBIFGQAAAAAASRRkAAAAAAAkUZABAAAAAJBEQQYAAAAAQBIFGQAAAAAASRRkAAAAAAAkUZABAAAAAJBEQQYAAAAAQBIFGQAAAAAASRRkAAAAAAAkUZABAAAAAJBEQQYAAAAAQBIFGQAAAAAASRRkVFJycrIGDhzo93UPHjwom82mOXPm+H1tACgP8hEAXCMfUZVRkOHS/v37NXToUDVp0kQRERGKiYlRWlqapkyZol9++cXs8Vxavny5xo8f79FtvvnmG912222KiopSXFyc7rvvPv3444++GRBAQAiGfNyyZYseeughpaamqlq1arLZbL4bDkDACPR8LCkp0Zw5c9SzZ08lJiaqevXqSklJ0bPPPquzZ8/6dlD4jc0wDMPsIVC1LFu2TL1795bdblf//v2VkpKioqIirV+/Xm+//bYGDhyoWbNmSZIcDodCQkJUrVo1v85oGIYcDoeqVaum0NBQSdLw4cM1bdo0lfdX+ocfftC1116r2NhYZWRkqLCwUC+++KIaNWqkLVu2KDw83JeHAMCCgiUfx48fr4kTJ6p169YqKCjQnj17yn1bAMEpGPKxsLBQ0dHRuuGGG9SjRw/VqVNHmzZt0ty5c9WxY0d98skn/EExAISZPQCqlgMHDqhPnz5KSkrSJ598ovr16zt/lp6ern379mnZsmXObXa7/ZL3efr0aVWvXt2rc9psNkVERFTqPiZOnKjTp09r+/btatSokSSpbdu2uuWWWzRnzhwNGTLEG6MCCBDBlI/Dhg3To48+qsjISA0fPlx79uzx0nQAAlGw5GN4eLg2bNig9u3bO7c9+OCDSk5O1pNPPqlVq1apS5cu3hgVZjKAC/z1r381JBkbNmwo1/5JSUnGgAEDnNdnz55tSDLWrFljDBs2zIiPjzdq1Kjh/Pny5cuNjh07GlFRUUZ0dLTRpk0bIysry+39ndepUyejU6dOzusHDhwwJBmzZ882DMMwBgwYYEgqc7mYOnXqGL179y6z/corrzQ6d+5cruMHEDyCKR8vlJ6e7tH+AIJPsObjebt27TIkGS+//LLHt0XVwxlklPLBBx+oSZMmpf4yVhEPPfSQ4uPj9cQTT+j06dOSpDlz5uj+++9Xy5YtNXbsWNWoUUOff/65VqxYoXvvvbdS6w0dOlRHjhxRdna23njjjUvuf/jwYeXm5qpNmzZlfta2bVstX768UvMACDzBko8A4Klgz8djx45JkmrXrl2peVA1UJDhlJ+fr8OHD+uOO+6o9H3FxcVp1apVzvd35OXlKSMjQ23bttWaNWtKvbzF8ML72tq1a6crr7xS2dnZ6tev3yX3P3r0qCSVegnQefXr19fPP/8sh8NRrpcAAQh8wZSPAOAJ8lF64YUXFBMTo27dulV6JpiPggyn/Px8SVJ0dHSl7+vBBx90hpskZWdnq6CgQGPGjCnz3g8zPszg/CcpuirA5+f75ZdfKMgAJAVXPgKAJ4I9HydOnKiPP/5Yr776qmrUqGH2OPACvuYJTjExMZKkgoKCSt9X48aNS13fv3+/JCklJaXS9+0NkZGRkn77FMXfO/8x/ef3AYBgykcA8EQw5+PChQv1P//zP3rggQc0bNgws8eBl1CQ4RQTE6OEhATt3r270vdV0XLp7q+BxcXFlRmnjPMvrT7/UusLHT16VHFxcZw9BuAUTPkIAJ4I1nzMzs5W//791b17d82YMcNn68D/KMgopUePHtq/f782bdrk1ftt2rSpJF0yPGvWrKlTp06V2X7o0KFLruHJS20aNGig+Ph4bdu2rczPtmzZomuuuabc9wUgOARLPgKAp4ItHz/77DP16tVLbdq00VtvvaWwMN61GkgoyCjlkUceUfXq1TV48GAdP368zM/379+vKVOmeHy/Xbt2VXR0tDIzM50vYT7vwg9ZaNq0qTZv3qyioiLntqVLlyonJ+eSa5z/rjxXAenKn//85zL3vWrVKu3Zs0e9e/cu130ACB7BlI8A4IlgysdvvvlG3bt3V3JyspYuXcpb8gIQf+5AKU2bNtW8efN09913q3nz5urfv79SUlJUVFSkjRs3atGiRRo4cKDH9xsTE6OXXnpJgwcP1nXXXad7771XNWvW1BdffKEzZ85o7ty5kqTBgwdr8eLFuu2223TXXXdp//79evPNN51/QbyY1NRUSVJGRoZuvfVWhYaGqk+fPm73HzdunBYtWqSbbrpJI0aMUGFhof7xj3+oVatWGjRokMfHCCCwBVM+Hjp0yPmVJ+dfafPss89KkpKSknTfffd5fJwAAlew5GNBQYFuvfVWnTx5UqNHj9ayZcvKPA7t2rXz+DhRxZj5Jcyouvbs2WM8+OCDRnJyshEeHm5ER0cbaWlpxtSpU42zZ88693P3Re9bt251eb/vv/++0b59eyMyMtKIiYkx2rZta8yfP7/UPv/85z+NBg0aGHa73UhLSzO2bdt2yS96NwzDOHfunPHwww8b8fHxhs1mK9cXve/evdvo2rWrcdlllxk1atQw+vbtaxw7dqx8DxKAoBQM+bh69WpDksvLhWsBwIUCPR/P397d5cJjgnXZDMMLXyIGAAAAAIDF8R5kAAAAAABEQQYAAAAAQBIFGQAAAAAASRRkAAAAAAAkUZABAAAAAJBEQQYAAAAAQJIU5u8FS0pKdOTIEUVHR8tms/l7eQAWZxiGCgoKlJCQoJCQwPobH/kIoDLIRwBwr7wZ6feCfOTIESUmJvp7WQABJicnRw0bNjR7DK8iHwF4A/kIAO5dKiP9XpCjo6MlSZ2aDFNYiN2va59Mre3X9aqC/CTz/oJsP2nOuiXh5qwrSWFnDPMWDxLFRWf1ddYzziwJJOSjf/16mXlnoc6ZtHbkTT+asq4k/bos+H7H/I189I1gzMdVz/zbtLU7Zj5gyrpm5bLE80d/KW9G+r0gn39ZTFiIXWGh/g240PAIv65XFYTazSvIoSYVVZt/f61KCT1HwPlLIL7Ejnz0r5Jw836HDLs5a4dWNy8gS4Lwd8ws5KN3BWM+xkSb+fzRnMfbrFyWeP7ob5fKyMB6gwoAAAAAABVEQQYAAAAAQBRkAAAAAAAkUZABAAAAAJBEQQYAAAAAQBIFGQAAAAAASRRkAAAAAAAkUZABAAAAAJBUwYI8bdo0JScnKyIiQtdff722bNni7bkAwJLIRwBwjXwEYAUeF+SFCxdq1KhRevLJJ7Vjxw5dffXVuvXWW5Wbm+uL+QDAMshHAHCNfARgFR4X5EmTJunBBx/UoEGD1KJFC82YMUOXXXaZ/v3vf/tiPgCwDPIRAFwjHwFYhUcFuaioSNu3b1eXLl3+7w5CQtSlSxdt2rTJ5W0cDofy8/NLXQAg0JCPAOAa+QjASjwqyCdOnFBxcbHq1q1banvdunV17Ngxl7fJzMxUbGys85KYmFjxaQGgiiIfAcA18hGAlfj8U6zHjh2rvLw85yUnJ8fXSwKAJZCPAOAa+QjALGGe7Fy7dm2Fhobq+PHjpbYfP35c9erVc3kbu90uu91e8QkBwALIRwBwjXwEYCUenUEODw9XamqqVq1a5dxWUlKiVatWqV27dl4fDgCsgnwEANfIRwBW4tEZZEkaNWqUBgwYoDZt2qht27aaPHmyTp8+rUGDBvliPgCwDPIRAFwjHwFYhccF+e6779aPP/6oJ554QseOHdM111yjFStWlPngBQAINuQjALhGPgKwCo8LsiQNHz5cw4cP9/YsAGB55CMAuEY+ArACn3+KNQAAAAAAVkBBBgAAAABAFGQAAAAAACRRkAEAAAAAkERBBgAAAABAEgUZAAAAAABJFGQAAAAAACRV8HuQveFkam2Fhkf4dc24LT/6db0L/dw23pR1I3MNU9aVpJIwmynrhp0275gBbyAf4Su/LjHvsS6OMOf/hFCHif8n8N+R15GP/pE6fpgp65opGJ8/VrvTvN9tM/8/uhTOIAMAAAAAIAoyAAAAAACSKMgAAAAAAEiiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSKMgAAAAAAEiiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSKMgAAAAAAEiqQEFet26dbr/9diUkJMhms2nJkiU+GAsArId8BADXyEcAVuFxQT59+rSuvvpqTZs2zRfzAIBlkY8A4Br5CMAqwjy9Qbdu3dStWzdfzAIAlkY+AoBr5CMAq/C4IHvK4XDI4XA4r+fn5/t6SQCwBPIRAFwjHwGYxecf0pWZmanY2FjnJTEx0ddLAoAlkI8A4Br5CMAsPi/IY8eOVV5envOSk5Pj6yUBwBLIRwBwjXwEYBafv8TabrfLbrf7ehkAsBzyEQBcIx8BmIXvQQYAAAAAQBU4g1xYWKh9+/Y5rx84cEA7d+5UXFycGjVq5NXhAMBKyEcAcI18BGAVHhfkbdu26aabbnJeHzVqlCRpwIABmjNnjtcGAwCrIR8BwDXyEYBVeFyQb7zxRhmG4YtZAMDSyEcAcI18BGAVvAcZAAAAAABRkAEAAAAAkERBBgAAAABAEgUZAAAAAABJFGQAAAAAACRRkAEAAAAAkERBBgAAAABAEgUZAAAAAABJUpjZA/jTz23jTVs7bsuPpqx7/MY6pqwrSSHnDFPWPVfdZsq6khR22pxjNpMRGtjrBQvy0b/Myopqd5rzWEuSlpj0O2befwkIEORjcOD5o3+Z8XyuvGtyBhkAAAAAAFGQAQAAAACQREEGAAAAAEASBRkAAAAAAEkUZAAAAAAAJFGQAQAAAACQREEGAAAAAEASBRkAAAAAAEkUZAAAAAAAJFGQAQAAAACQREEGAAAAAECShwU5MzNT1113naKjo1WnTh3deeed+vbbb301GwBYBvkIAO6RkQCswqOCvHbtWqWnp2vz5s3Kzs7Wr7/+qq5du+r06dO+mg8ALIF8BAD3yEgAVhHmyc4rVqwodX3OnDmqU6eOtm/fro4dO3p1MACwEvIRANwjIwFYhUcF+ffy8vIkSXFxcW73cTgccjgczuv5+fmVWRIALIF8BAD3LpWR5CMAs1T4Q7pKSko0cuRIpaWlKSUlxe1+mZmZio2NdV4SExMruiQAWAL5CADulScjyUcAZqlwQU5PT9fu3bu1YMGCi+43duxY5eXlOS85OTkVXRIALIF8BAD3ypOR5CMAs1ToJdbDhw/X0qVLtW7dOjVs2PCi+9rtdtnt9goNBwBWQz4CgHvlzUjyEYBZPCrIhmHo4Ycf1rvvvqs1a9aocePGvpoLACyFfAQA98hIAFbhUUFOT0/XvHnz9N577yk6OlrHjh2TJMXGxioyMtInAwKAFZCPAOAeGQnAKjx6D/L06dOVl5enG2+8UfXr13deFi5c6Kv5AMASyEcAcI+MBGAVHr/EGgBQFvkIAO6RkQCsosKfYg0AAAAAQCChIAMAAAAAIAoyAAAAAACSKMgAAAAAAEiiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSPPweZG8qaBiiULt/+3nET+Z9B9/xG+uYsm54oXnHXBRlM2XdsNN816I/2YoDez0zkI/wlV+XxJs9gv/xX0JAIR/9o/asTaasK0knhrQzZd1gfP5o5v8JZrSE8j6H5AwyAAAAAACiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSKMgAAAAAAEiiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSKMgAAAAAAEiiIAMAAAAAIImCDAAAAACAJA8L8vTp09W6dWvFxMQoJiZG7dq104cffuir2QDAMshHAHCPjARgFR4V5IYNG+q5557T9u3btW3bNt18882644479NVXX/lqPgCwBPIRANwjIwFYRZgnO99+++2lrk+YMEHTp0/X5s2b1bJlS68OBgBWQj4CgHtkJACr8KggX6i4uFiLFi3S6dOn1a5dO7f7ORwOORwO5/X8/PyKLgkAlkA+AoB75clI8hGAWTz+kK4vv/xSUVFRstvt+utf/6p3331XLVq0cLt/ZmamYmNjnZfExMRKDQwAVRX5CADueZKR5CMAs3hckK+66irt3LlTn332mYYNG6YBAwbo66+/drv/2LFjlZeX57zk5ORUamAAqKrIRwBwz5OMJB8BmMXjl1iHh4fr8ssvlySlpqZq69atmjJlimbOnOlyf7vdLrvdXrkpAcACyEcAcM+TjCQfAZil0t+DXFJSUuo9IgCA35CPAOAeGQmgKvLoDPLYsWPVrVs3NWrUSAUFBZo3b57WrFmjjz76yFfzAYAlkI8A4B4ZCcAqPCrIubm56t+/v44eParY2Fi1bt1aH330kW655RZfzQcAlkA+AoB7ZCQAq/CoIL/++uu+mgMALI18BAD3yEgAVlHp9yADAAAAABAIKMgAAAAAAIiCDAAAAACAJAoyAAAAAACSKMgAAAAAAEiiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSpDDTFj4jhRb7d81zkTb/LniBsF8MU9YtijLvmEPOmbNuiWm/1YB3kI+BL7ww+P5PALyBfPSPE0PambIuUBVwBhkAAAAAAFGQAQAAAACQREEGAAAAAEASBRkAAAAAAEkUZAAAAAAAJFGQAQAAAACQREEGAAAAAEASBRkAAAAAAEkUZAAAAAAAJFGQAQAAAACQVMmC/Nxzz8lms2nkyJFeGgcAAgP5CACukY8AqrIKF+StW7dq5syZat26tTfnAQDLIx8BwDXyEUBVV6GCXFhYqL59++q1115TzZo1vT0TAFgW+QgArpGPAKygQgU5PT1d3bt3V5cuXbw9DwBYGvkIAK6RjwCsIMzTGyxYsEA7duzQ1q1by7W/w+GQw+FwXs/Pz/d0SQCwBPIRAFwjHwFYhUdnkHNycjRixAhlZWUpIiKiXLfJzMxUbGys85KYmFihQQGgKiMfAcA18hGAldgMwzDKu/OSJUvUq1cvhYaGOrcVFxfLZrMpJCREDoej1M8k138BTExMVPNhExVqL19IBoKwX8r9MAeMkHPmrFvi8esiYCXFRWf15ezHlJeXp5iYGLPHcSIfKy4Y8zG80JxjLoqymbIu/IN8DDzBmI+Ar5Q3Iz2qEp07d9aXX35ZatugQYPUrFkzPfroo2XCTZLsdrvsdrsnywCA5ZCPAOAa+QjASjwqyNHR0UpJSSm1rXr16qpVq1aZ7QAQTMhHAHCNfARgJRX+HmQAAAAAAAJJpd+tuWbNGi+MAQCBh3wEANfIRwBVFWeQAQAAAAAQBRkAAAAAAEkUZAAAAAAAJFGQAQAAAACQREEGAAAAAEASBRkAAAAAAEkUZAAAAAAAJFGQAQAAAACQJIWZPYA/hf1imLb2uUibKeuaecwlJv12hRead8xFUeb8OwOVRT76l1lZUXvWJlPWlaQTQ9qZtjZQGeRjcOD5I87jDDIAAAAAAKIgAwAAAAAgiYIMAAAAAIAkCjIAAAAAAJIoyAAAAAAASKIgAwAAAAAgiYIMAAAAAIAkCjIAAAAAAJIoyAAAAAAASKIgAwAAAAAgiYIMAAAAAIAkDwvy+PHjZbPZSl2aNWvmq9kAwDLIRwBwj4wEYBVhnt6gZcuW+vjjj//vDsI8vgsACEjkIwC4R0YCsAKPkyksLEz16tXzxSwAYGnkIwC4R0YCsAKP34O8d+9eJSQkqEmTJurbt6++//77i+7vcDiUn59f6gIAgYh8BAD3PMlI8hGAWTwqyNdff73mzJmjFStWaPr06Tpw4ID+8Ic/qKCgwO1tMjMzFRsb67wkJiZWemgAqGrIRwBwz9OMJB8BmMVmGIZR0RufOnVKSUlJmjRpkh544AGX+zgcDjkcDuf1/Px8JSYmqvmwiQq1R1R06QoJ+6XCh1pp5yJtpqxr5jGbJbzQvGMuijLn3zmYFBed1ZezH1NeXp5iYmLMHsct8rH8yEf/qT1rk2lrnxjSzrS1g4VV8lG6dEaSj78hH/2H54+Br7wZWalPR6hRo4auvPJK7du3z+0+drtddru9MssAgOWQjwDg3qUyknwEYJZKfQ9yYWGh9u/fr/r163trHgAICOQjALhHRgKoqjwqyH//+9+1du1aHTx4UBs3blSvXr0UGhqqe+65x1fzAYAlkI8A4B4ZCcAqPHqJ9Q8//KB77rlHP/30k+Lj49WhQwdt3rxZ8fHxvpoPACyBfAQA98hIAFbhUUFesGCBr+YAAEsjHwHAPTISgFVU6j3IAAAAAAAECgoyAAAAAACiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSKMgAAAAAAEiiIAMAAAAAIImCDAAAAACAJCnMtIXPGgotMfy65tlaNr+ud6GIn/x7rMGsKMq8f2fAG8hH+MqJIe1MWzvknDnrlpj2TAe+QD7CV3j+iPM4gwwAAAAAgCjIAAAAAABIoiADAAAAACCJggwAAAAAgCQKMgAAAAAAkijIAAAAAABIoiADAAAAACCJggwAAAAAgCQKMgAAAAAAkijIAAAAAABIqkBBPnz4sPr166datWopMjJSrVq10rZt23wxGwBYCvkIAK6RjwCsIsyTnU+ePKm0tDTddNNN+vDDDxUfH6+9e/eqZs2avpoPACyBfAQA18hHAFbiUUF+/vnnlZiYqNmzZzu3NW7c2OtDAYDVkI8A4Br5CMBKPHqJ9fvvv682bdqod+/eqlOnjq699lq99tprF72Nw+FQfn5+qQsABBryEQBcIx8BWIlHBfm7777T9OnTdcUVV+ijjz7SsGHDlJGRoblz57q9TWZmpmJjY52XxMTESg8NAFUN+QgArpGPAKzEZhiGUd6dw8PD1aZNG23cuNG5LSMjQ1u3btWmTZtc3sbhcMjhcDiv5+fnKzExUa0GTVBoeEQlRvfc2Vo2v653oYifyv0wA7iI4qKz+nL2Y8rLy1NMTIzZ4ziRjxVHPgaHkHPmrFvi0ZvJrI189A3yEQgM5c1Ij84g169fXy1atCi1rXnz5vr+++/d3sZutysmJqbUBQACDfkIAK6RjwCsxKOCnJaWpm+//bbUtj179igpKcmrQwGA1ZCPAOAa+QjASjwqyH/729+0efNmTZw4Ufv27dO8efM0a9Yspaen+2o+ALAE8hEAXCMfAViJRwX5uuuu07vvvqv58+crJSVFzzzzjCZPnqy+ffv6aj4AsATyEQBcIx8BWInHH13Ro0cP9ejRwxezAIClkY8A4Br5CMAqPDqDDAAAAABAoKIgAwAAAAAgCjIAAAAAAJIoyAAAAAAASKIgAwAAAAAgiYIMAAAAAIAkCjIAAAAAAJIoyAAAAAAASJLCzB7An8LOmLf2uUibKeuG/WKYsi78y2biP7Nhzq82vIx8hK+VmPSMw1ZizrqSZHAaIiCQjwhUtWdtMm3tE0Pambb2pRDdAAAAAACIggwAAAAAgCQKMgAAAAAAkijIAAAAAABIoiADAAAAACCJggwAAAAAgCQKMgAAAAAAkijIAAAAAABIoiADAAAAACCJggwAAAAAgCQKMgAAAAAAkjwsyMnJybLZbGUu6enpvpoPACyBfAQA98hIAFYR5snOW7duVXFxsfP67t27dcstt6h3795eHwwArIR8BAD3yEgAVuFRQY6Pjy91/bnnnlPTpk3VqVMnrw4FAFZDPgKAe2QkAKvwqCBfqKioSG+++aZGjRolm83mdj+HwyGHw+G8np+fX9ElAcASyEcAcK88GUk+AjBLhT+ka8mSJTp16pQGDhx40f0yMzMVGxvrvCQmJlZ0SQCwBPIRANwrT0aSjwDMUuGC/Prrr6tbt25KSEi46H5jx45VXl6e85KTk1PRJQHAEshHAHCvPBlJPgIwS4VeYn3o0CF9/PHHeueddy65r91ul91ur8gyAGA55CMAuFfejCQfAZilQmeQZ8+erTp16qh79+7engcALI18BAD3yEgAVZ3HBbmkpESzZ8/WgAEDFBZW4c/4AoCAQz4CgHtkJAAr8Lggf/zxx/r+++91//33+2IeALAs8hEA3CMjAViBx3++69q1qwzD8MUsAGBp5CMAuEdGArCCCn+KNQAAAAAAgYSCDAAAAACAKMgAAAAAAEiiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSKMgAAAAAAEiqwPcgW1nYL+Z9997ZWjZT1g37xZRl4WeGOb9ekiRbSWCvFyzIRwQqg1MBqCTyEYHqxJB2pq1txvO58q7JfxsAAAAAAIiCDAAAAACAJAoyAAAAAACSKMgAAAAAAEiiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSKMgAAAAAAEiiIAMAAAAAIImCDAAAAACAJAoyAAAAAACSPCzIxcXFevzxx9W4cWNFRkaqadOmeuaZZ2QYhq/mAwBLIB8BwDXyEYCVhHmy8/PPP6/p06dr7ty5atmypbZt26ZBgwYpNjZWGRkZvpoRAKo88hEAXCMfAViJRwV548aNuuOOO9S9e3dJUnJysubPn68tW7b4ZDgAsAryEQBcIx8BWIlHL7Fu3769Vq1apT179kiSvvjiC61fv17dunVzexuHw6H8/PxSFwAINOQjALhGPgKwEo/OII8ZM0b5+flq1qyZQkNDVVxcrAkTJqhv375ub5OZmamnnnqq0oMCQFVGPgKAa+QjACvx6AzyW2+9paysLM2bN087duzQ3Llz9eKLL2ru3LlubzN27Fjl5eU5Lzk5OZUeGgCqGvIRAFwjHwFYiUdnkEePHq0xY8aoT58+kqRWrVrp0KFDyszM1IABA1zexm63y263V35SAKjCyEcAcI18BGAlHp1BPnPmjEJCSt8kNDRUJSUlXh0KAKyGfAQA18hHAFbi0Rnk22+/XRMmTFCjRo3UsmVLff7555o0aZLuv/9+X80HAJZAPgKAa+QjACvxqCBPnTpVjz/+uB566CHl5uYqISFBQ4cO1RNPPOGr+QDAEshHAHCNfARgJR4V5OjoaE2ePFmTJ0/20TgAYE3kIwC4Rj4CsBKP3oMMAAAAAECgoiADAAAAACAKMgAAAAAAkijIAAAAAABIoiADAAAAACCJggwAAAAAgCQKMgAAAAAAkjz8HmRvMAxDklRcdNbfS5uq2GEzZ90iw5R1ETxsJf5d73x2nM+SQEI++nld8hEBhnwMPOQjApW/nz9K5c9Im+HnFP3hhx+UmJjozyUBBKCcnBw1bNjQ7DG8inwE4A3kIwC4d6mM9HtBLikp0ZEjRxQdHS2bzbO/iuXn5ysxMVE5OTmKiYnx0YRVS7Adc7Adr8Qxe3rMhmGooKBACQkJCgkJrHeJkI+e4ZgD/5iD7Xgl8tEd8tEzHDPHHIgqe7zlzUi/v8Q6JCSk0n/VjImJCYpfggsF2zEH2/FKHLMnYmNjfTCN+cjHiuGYA1+wHa9EPv4e+VgxHHNwCLZjrszxlicjA+vPiwAAAAAAVBAFGQAAAAAAWawg2+12Pfnkk7Lb7WaP4jfBdszBdrwSxwzvCMbHlGMOfMF2vFJwHrOvBeNjyjEHh2A7Zn8dr98/pAsAAAAAgKrIUmeQAQAAAADwFQoyAAAAAACiIAMAAAAAIImCDAAAAACAJAsV5GnTpik5OVkRERG6/vrrtWXLFrNH8pnMzExdd911io6OVp06dXTnnXfq22+/NXssv3ruuedks9k0cuRIs0fxqcOHD6tfv36qVauWIiMj1apVK23bts3ssXyiuLhYjz/+uBo3bqzIyEg1bdpUzzzzjPicQO8IlowkH8nHQEVG+k6w5KNERpKPgcnf+WiJgrxw4UKNGjVKTz75pHbs2KGrr75at956q3Jzc80ezSfWrl2r9PR0bd68WdnZ2fr111/VtWtXnT592uzR/GLr1q2aOXOmWrdubfYoPnXy5EmlpaWpWrVq+vDDD/X111/rn//8p2rWrGn2aD7x/PPPa/r06XrllVf0zTff6Pnnn9cLL7ygqVOnmj2a5QVTRpKP5GOgIiN9I5jyUQrujCQfyUevMSygbdu2Rnp6uvN6cXGxkZCQYGRmZpo4lf/k5uYakoy1a9eaPYrPFRQUGFdccYWRnZ1tdOrUyRgxYoTZI/nMo48+anTo0MHsMfyme/fuxv33319q25/+9Cejb9++Jk0UOII5I8nHwBRs+WgYZKSvBHM+GkbwZCT5GNj8nY9V/gxyUVGRtm/fri5duji3hYSEqEuXLtq0aZOJk/lPXl6eJCkuLs7kSXwvPT1d3bt3L/XvHajef/99tWnTRr1791adOnV07bXX6rXXXjN7LJ9p3769Vq1apT179kiSvvjiC61fv17dunUzeTJrC/aMJB8DU7Dlo0RG+kKw56MUPBlJPpKP3hTmk3v1ohMnTqi4uFh169Yttb1u3br6z3/+Y9JU/lNSUqKRI0cqLS1NKSkpZo/jUwsWLNCOHTu0detWs0fxi++++07Tp0/XqFGjNG7cOG3dulUZGRkKDw/XgAEDzB7P68aMGaP8/Hw1a9ZMoaGhKi4u1oQJE9S3b1+zR7O0YM5I8jFwBVs+SmSkLwRzPkrBk5HkI/nobVW+IAe79PR07d69W+vXrzd7FJ/KycnRiBEjlJ2drYiICLPH8YuSkhK1adNGEydOlCRde+212r17t2bMmBGQAffWW28pKytL8+bNU8uWLbVz506NHDlSCQkJAXm88D3yMXAFWz5KZCS8LxgyknwkH31xzFW+INeuXVuhoaE6fvx4qe3Hjx9XvXr1TJrKP4YPH66lS5dq3bp1atiwodnj+NT27duVm5ur//qv/3JuKy4u1rp16/TKK6/I4XAoNDTUxAm9r379+mrRokWpbc2bN9fbb79t0kS+NXr0aI0ZM0Z9+vSRJLVq1UqHDh1SZmZmwAa6PwRrRpKP5GOgISO9L1jzUQqejCQff0M+eleVfw9yeHi4UlNTtWrVKue2kpISrVq1Su3atTNxMt8xDEPDhw/Xu+++q08++USNGzc2eySf69y5s7788kvt3LnTeWnTpo369u2rnTt3Bly4SVJaWlqZr17Ys2ePkpKSTJrIt86cOaOQkNKRExoaqpKSEpMmCgzBlpHkI/kYqMhI7wu2fJSCLyPJx9+Qj17mk4/+8rIFCxYYdrvdmDNnjvH1118bQ4YMMWrUqGEcO3bM7NF8YtiwYUZsbKyxZs0a4+jRo87LmTNnzB7NrwL9Uwi3bNlihIWFGRMmTDD27t1rZGVlGZdddpnx5ptvmj2aTwwYMMBo0KCBsXTpUuPAgQPGO++8Y9SuXdt45JFHzB7N8oIpI8nH35CPgYeM9I1gykfDICMNg3wMRP7OR0sUZMMwjKlTpxqNGjUywsPDjbZt2xqbN282eySfkeTyMnv2bLNH86tADzjDMIwPPvjASElJMex2u9GsWTNj1qxZZo/kM/n5+caIESOMRo0aGREREUaTJk2Mxx57zHA4HGaPFhCCJSPJx9+Qj4GHjPSdYMlHwyAjDYN8DET+zkebYRiGb85NAwAAAABgHVX+PcgAAAAAAPgDBRkAAAAAAFGQAQAAAACQREEGAAAAAEASBRkAAAAAAEkUZAAAAAAAJFGQAQAAAACQREEGAAAAAEASBRkAAAAAAEkUZAAAAAAAJFGQAQAAAACQREEGAAAAAECS9P8A08cbT6RnZxwAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "def one_step_de_sim_expanded(spec_conc, inputs, outputs, forward_rates, reverse_rates):\n", + " concentration_factors_in = jnp.prod(\n", + " jnp.power(spec_conc, (inputs)), axis=1)\n", + " concentration_factors_out = jnp.prod(\n", + " jnp.power(spec_conc, (outputs)), axis=1)\n", + " forward_delta = concentration_factors_in * forward_rates\n", + " reverse_delta = concentration_factors_out * reverse_rates\n", + " return (forward_delta - reverse_delta) @ (outputs - inputs)\n", + "\n", + "\n", + "# bb = partial(bioreaction_sim_dfx_expanded,\n", + "# t0=t0, t1=t1, dt0=dt0,\n", + "# signal=None, signal_onehot=signal_onehot,\n", + "# forward_rates=forward_rates,\n", + "# inputs=inputs,\n", + "# outputs=outputs,\n", + "# solver=dfx.Tsit5(),\n", + "# saveat=dfx.SaveAt(\n", + "# ts=ts),\n", + "# max_steps=max_steps,\n", + "# stepsize_controller=make_piecewise_stepcontrol(\n", + "# t0=t0, t1=t1, dt0=dt0, dt1=dt1)\n", + "# )\n", + "\n", + "Jbb = jax.vmap(jacrev(partial(one_step_de_sim_expanded,\n", + " forward_rates=forward_rates[0],\n", + " inputs=inputs,\n", + " outputs=outputs)))\n", + "\n", + "\n", + "sol_jac = Jbb(y00, reverse_rates=reverse_rates)\n", + "\n", + "fig = plt.figure(figsize=(4*n_circuits, 4))\n", + "for idx_circuit in range(n_circuits):\n", + " ax = fig.add_subplot(1, n_circuits, idx_circuit+1)\n", + " ax.set_title(f'Circuit {idx_circuit}')\n", + " plt.imshow(sol_jac[idx_circuit])\n", + " plt.clim([sol_jac.min(), sol_jac.max()])\n", + "plt.suptitle('Jacobian of system')" + ] + }, + { + "cell_type": "code", + "execution_count": 104, + "metadata": {}, + "outputs": [], + "source": [ + "def one_step_de_sim_expanded(t, spec_conc, args, inputs, outputs, forward_rates, reverse_rates):\n", + " concentration_factors_in = jnp.prod(\n", + " jnp.power(spec_conc, (inputs)), axis=1)\n", + " concentration_factors_out = jnp.prod(\n", + " jnp.power(spec_conc, (outputs)), axis=1)\n", + " forward_delta = concentration_factors_in * forward_rates\n", + " reverse_delta = concentration_factors_out * reverse_rates\n", + " return (forward_delta - reverse_delta) @ (outputs - inputs)\n", + "\n", + "\n", + "\n", + "def wrap(y0,\n", + " reverse_rates,\n", + " solver=dfx.Tsit5(),\n", + " saveat=dfx.SaveAt(\n", + " ts=ts),\n", + " max_steps=max_steps,\n", + " stepsize_controller=make_piecewise_stepcontrol(\n", + " t0=t0, t1=t1, dt0=dt0, dt1=dt1)):\n", + " term = dfx.ODETerm(\n", + " jax.jacfwd(\n", + " partial(one_step_de_sim_expanded,\n", + " forward_rates=forward_rates[0],\n", + " inputs=inputs,\n", + " outputs=outputs,\n", + " reverse_rates=reverse_rates)\n", + " )\n", + " # partial(bioreaction_sim_expanded,\n", + " )\n", + " return dfx.diffeqsolve(term, solver,\n", + " t0=t0, t1=t1, dt0=None,\n", + " y0=y0.squeeze(),\n", + " saveat=saveat, max_steps=max_steps,\n", + " stepsize_controller=stepsize_controller)\n", + " \n", + "\n", + "y01 = y00.copy()\n", + "y01[:, np.array(idxs_signal)] = y00[:, np.array(idxs_signal)] * signal_target\n", + "# sol_signal = wrap(y01[0], reverse_rates[0])\n", + "sol_signal = jax.vmap(wrap)(y01, reverse_rates)" + ] + }, + { + "cell_type": "code", + "execution_count": 105, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 0.98, 'Jacobian of system')" + ] + }, + "execution_count": 105, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9oAAAGHCAYAAABPmCpHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABELElEQVR4nO3de3gTZd7/8U8KbSiUtBZoQ6UtCCogVF0QiPiAj1QO1jMeUOS0CIpFQHYRcRUPq9RFH0VdFV13gVXwgAdcEcQuR5HKSZGDioIoKKRFsUkBaaG9f3/sjyyxLTbttEmT9+u65rqamXtmvncu/Ox+k5mMzRhjBAAAAAAALBEV7AIAAAAAAAgnNNoAAAAAAFiIRhsAAAAAAAvRaAMAAAAAYCEabQAAAAAALESjDQAAAACAhWi0AQAAAACwEI02AAAAAAAWotEGAAAAAMBCNNoAAEhasWKFbDab3njjjd8cO3z4cLVu3br2i6qBgwcP6uabb5bT6ZTNZtOECROCXRIAABGDRhsAEHSzZ8+WzWbThg0bgl1K2Jg2bZpmz56tMWPG6KWXXtKQIUOCXZIkadGiRbr//vuDXQYAALWqYbALAACgvvnb3/6msrKyYJdxUsuWLVOPHj103333BbsUP4sWLdIzzzxDsw0ACGt8ow0AQICio6Nlt9uDXcZJFRQUKCEhIdhlAAAQkWi0AQAhZ/PmzRo+fLhOO+00NWrUSE6nU7///e/1008/lRv7ww8/aOTIkUpJSZHdblebNm00ZswYlZSU+MZ88803uvbaa5WYmKjGjRurR48eeu+99yo8d2lpqe6++245nU41adJEl19+ufbs2eM3pqJ7tB977DGdf/75atasmWJjY9WlS5cK7/e22WwaO3asFixYoE6dOslut+uss87S+++/X6X3pqCgQCNHjlRycrIaNWqks88+W3PmzPFtP36v+a5du/Tee+/JZrPJZrPp22+/rfSYubm5uuCCC5SQkKC4uDideeaZuvvuuyX9517vJk2aaPz48eX2+/7779WgQQPl5ORIko4ePaoHHnhAp59+uho1aqRmzZrpggsuUG5uru99e+aZZ3zvw/HluLKyMs2YMUNnnXWWGjVqpOTkZN1yyy36+eef/c7bunVrXXrppVqxYoW6du2q2NhYde7cWStWrJAkvfXWW+rcubMaNWqkLl266NNPP63SewsAgFW4dBwAEHJyc3P1zTffaMSIEXI6ndq2bZteeOEFbdu2TR9//LGvOdu7d6+6deumwsJCjR49Wu3bt9cPP/ygN954Q4cPH1ZMTIzy8/N1/vnn6/Dhwxo3bpyaNWumOXPm6PLLL9cbb7yhq666yu/cDz/8sGw2myZPnqyCggLNmDFDmZmZ2rRpk2JjYyut+cknn9Tll1+uwYMHq6SkRK+++qquvfZaLVy4UFlZWX5jV69erbfeeku33XabmjZtqqeeekoDBw7U7t271axZs0rP8csvv+jCCy/Ujh07NHbsWLVp00bz58/X8OHDVVhYqPHjx6tDhw566aWXdMcdd6hVq1b6wx/+IElq0aJFhcfctm2bLr30UmVkZOjBBx+U3W7Xjh079NFHH0mS4uLidNVVV+m1117T448/rgYNGvj2feWVV2SM0eDBgyVJ999/v3JycnTzzTerW7du8nq92rBhgz755BNdfPHFuuWWW7R3717l5ubqpZdeKlfLLbfcotmzZ2vEiBEaN26cdu3apb/+9a/69NNP9dFHHyk6Oto3dseOHbrxxht1yy236KabbtJjjz2myy67TDNnztTdd9+t2267TZKUk5Oj6667Ttu3b1dUFN8vAADqiAEAIMhmzZplJJn169cbY4w5fPhwuTGvvPKKkWRWrVrlWzd06FATFRXl2+9EZWVlxhhjJkyYYCSZDz/80LetqKjItGnTxrRu3dqUlpYaY4xZvny5kWROPfVU4/V6fWNff/11I8k8+eSTvnXDhg0z6enpfuf7dc0lJSWmU6dO5qKLLvJbL8nExMSYHTt2+NZ99tlnRpJ5+umnK36D/r8ZM2YYSebll1/2O4/L5TJxcXF+daenp5usrKyTHs8YY5544gkjyezfv7/SMUuWLDGSzOLFi/3WZ2RkmN69e/ten3322b95zuzsbFPR//348MMPjSQzd+5cv/Xvv/9+ufXp6elGklmzZk25GmNjY813333nW//8888bSWb58uUnrQsAACvx0S4AIOSc+M3xkSNH9OOPP6pHjx6SpE8++UTSfy4zXrBggS677DJ17dq13DGOf+u9aNEidevWTRdccIFvW1xcnEaPHq1vv/1Wn3/+ud9+Q4cOVdOmTX2vr7nmGrVs2VKLFi2qcs0///yzPB6P/ud//sdX74kyMzPVtm1b3+uMjAw5HA598803Jz3HokWL5HQ6dcMNN/jWRUdHa9y4cTp48KBWrlx50v0rcvw+7nfeeafSH3jLzMxUSkqK5s6d61u3detWbd68WTfddJPfsbZt26avv/464Drmz5+v+Ph4XXzxxfrxxx99S5cuXRQXF6fly5f7je/YsaNcLpfvdffu3SVJF110kdLS0sqt/633FgAAK9FoAwBCzoEDBzR+/HglJycrNjZWLVq0UJs2bSRJHo9HkrR//355vV516tTppMf67rvvdOaZZ5Zb36FDB9/2E51++ul+r202m9q1a3fSe5wlaeHCherRo4caNWqkxMREtWjRQs8995yv3hOd2Aged8opp5S7F7miuZx++unlLoGubC5Vcf3116tnz566+eablZycrEGDBun111/3a7qjoqI0ePBgLViwQIcPH5YkzZ07V40aNdK1117rG/fggw+qsLBQZ5xxhjp37qxJkyZp8+bNVarj66+/lsfjUVJSklq0aOG3HDx4UAUFBX7jf/0exsfHS5JSU1MrXP9b7y0AAFbiHm0AQMi57rrrtGbNGk2aNEnnnHOO4uLiVFZWpv79+4fkY7U+/PBDXX755erVq5eeffZZtWzZUtHR0Zo1a5bmzZtXbvyJ9zmfyBhT26WWExsbq1WrVmn58uV677339P777+u1117TRRddpA8++MBX69ChQ/Xoo49qwYIFuuGGGzRv3jxdeumlvkZWknr16qWdO3fqnXfe0QcffKAXX3xRTzzxhGbOnKmbb775pHWUlZUpKSnJ71vzE/36HvPK3sNQem8BAJGLRhsAEFJ+/vlnLV26VA888ICmTp3qW//ry5FbtGghh8OhrVu3nvR46enp2r59e7n1X375pW/7iX59HmOMduzYoYyMjErP8eabb6pRo0ZasmSJ32O/Zs2addLaApWenq7NmzerrKzM71vtyuZSVVFRUerTp4/69Omjxx9/XNOmTdOf/vQnLV++XJmZmZKkTp066dxzz9XcuXPVqlUr7d69W08//XS5YyUmJmrEiBEaMWKEDh48qF69eun+++/3Ndon/sr4idq2bat///vf6tmz50l/dA4AgPqAS8cBACHl+DeSv/4GcsaMGX6vo6KidOWVV+rdd9/Vhg0byh3n+P6XXHKJ1q1bp7y8PN+2Q4cO6YUXXlDr1q3VsWNHv/3++c9/qqioyPf6jTfe0L59+zRgwICT1myz2VRaWupb9+2332rBggUnn2yALrnkErndbr322mu+dceOHdPTTz+tuLg49e7dO+BjHjhwoNy6c845R5JUXFzst37IkCH64IMPNGPGDDVr1qzce/Lrx6/FxcWpXbt2fsdp0qSJJKmwsNBv7HXXXafS0lL9+c9/LlfPsWPHyo0HACCU8Y02ACCkOBwO9erVS9OnT9fRo0d16qmn6oMPPtCuXbvKjZ02bZo++OAD9e7dW6NHj1aHDh20b98+zZ8/X6tXr1ZCQoLuuusuvfLKKxowYIDGjRunxMREzZkzR7t27dKbb75Z7n7nxMREXXDBBRoxYoTy8/M1Y8YMtWvXTqNGjaq05qysLD3++OPq37+/brzxRhUUFOiZZ55Ru3btqnyPclWMHj1azz//vIYPH66NGzeqdevWeuONN/TRRx9pxowZfj/iVlUPPvigVq1apaysLKWnp6ugoEDPPvusWrVq5fcDcpJ044036s4779Tbb7+tMWPG+D1uS/rPD5RdeOGF6tKlixITE7Vhwwa98cYbGjt2rG9Mly5dJEnjxo1Tv3791KBBAw0aNEi9e/fWLbfcopycHG3atEl9+/ZVdHS0vv76a82fP19PPvmkrrnmmmq8awAA1D0abQBA0B3/9vn4t9nz5s3T7bffrmeeeUbGGPXt21eLFy9WSkqK336nnnqq1q5dq3vvvVdz586V1+vVqaeeqgEDBqhx48aSpOTkZK1Zs0aTJ0/W008/rSNHjigjI0PvvvtuuedbS9Ldd9+tzZs3KycnR0VFRerTp4+effZZ3/EqctFFF+nvf/+7HnnkEU2YMEFt2rTRX/7yF3377beWNtqxsbFasWKF7rrrLs2ZM0der1dnnnmmZs2apeHDh1frmJdffrm+/fZb/eMf/9CPP/6o5s2bq3fv3nrggQf87r+W/vNe9u3bV4sWLdKQIUPKHWvcuHH617/+pQ8++EDFxcVKT0/XQw89pEmTJvnGXH311br99tv16quv6uWXX5YxRoMGDZIkzZw5U126dNHzzz+vu+++Ww0bNlTr1q110003qWfPntWaHwAAwWAz/DoIACDInnrqKY0fP147duzwe+wVQs9VV12lLVu2aMeOHcEuBQCAkMU92gCAoFu/fr2aNGlS7R/zQt3Yt2+f3nvvvQq/zQYAAP/FpeMAgKB58803tWLFCs2dO1c333yzGjbkf5ZC0a5du/TRRx/pxRdfVHR0tG655ZZglwQAQEjj/9EAAILmj3/8o4qKijRy5Eg98cQTwS4HlVi5cqVGjBihtLQ0zZkzR06nM9glAQAQ0rh0HEHXunXrav+IT018++23stlsmj17dp2fG8B/7Nq1Sz/++KP+9re/+R77hP8KlXwcPny4jDH67rvv+OVvACEjVDISqAiNNmrNzp07dcstt+i0005To0aN5HA41LNnTz355JP65Zdfgl1ehRYtWqT7778/oH2++OIL9e/fX3FxcUpMTNSQIUO0f//+2ikQQFiIhHxct26dbrvtNnXp0kXR0dGy2Wy1VxyAsBLuGVlWVqbZs2fr8ssvV2pqqpo0aaJOnTrpoYce0pEjR2q3UNQZfnUcteK9997TtddeK7vdrqFDh6pTp04qKSnR6tWr9eabb2r48OF64YUXJEnFxcWKiooq9zzW2maMUXFxsaKjo32PFBo7dqzvcUJV8f333+vcc89VfHy8xo0bp4MHD+qxxx5TWlqa1q1bp5iYmNqcAoB6KFLy8f7779e0adOUkZGhoqIiffXVV1XeF0DkioSMPHjwoJo2baoePXro0ksvVVJSkvLy8jRnzhz16tVLy5Yt48PJMMA92rDcrl27NGjQIKWnp2vZsmVq2bKlb1t2drZ27Nih9957z7fObrf/5jEPHTpk+WWlNptNjRo1qtExpk2bpkOHDmnjxo1KS0uTJHXr1k0XX3yxZs+erdGjR1tRKoAwEUn5OGbMGE2ePFmxsbEaO3asvvrqK4uqAxCuIiUjY2Ji9NFHH+n888/3rRs1apRat26t++67T0uXLlVmZqYVpSKYDGCxW2+91UgyH330UZXGp6enm2HDhvlez5o1y0gyK1asMGPGjDEtWrQwCQkJvu2LFi0yvXr1MnFxcaZp06ama9euZu7cuZUe77jevXub3r17+17v2rXLSDKzZs0yxhgzbNgwI6nccjJJSUnm2muvLbf+jDPOMH369KnS/AFEjkjKxxNlZ2cHNB5AZIrUjDxu8+bNRpJ56qmnAt4XoYdvtGG5d999V6eddprfp3TVcdttt6lFixaaOnWqDh06JEmaPXu2fv/73+uss87SlClTlJCQoE8//VTvv/++brzxxhqd75ZbbtHevXuVm5url1566TfH//DDDyooKFDXrl3LbevWrZsWLVpUo3oAhJ9IyUcAqI5Iz0i32y1Jat68eY3qQWig0YalvF6vfvjhB11xxRU1PlZiYqKWLl3qu/fF4/Fo3Lhx6tatm1asWOF3yY6x4L4/l8ulM844Q7m5ubrpppt+c/y+ffskye+ypuNatmypAwcOqLi4uEqXNQEIf5GUjwAQKDJSmj59uhwOhwYMGFDjmhB8NNqwlNfrlSQ1bdq0xscaNWqULyAlKTc3V0VFRbrrrrvK3RcTjB+MOP6rlxU10sfr++WXX2i0AUiKrHwEgEBFekZOmzZN//73v/Xss88qISEh2OXAAjzeC5ZyOBySpKKiohofq02bNn6vd+7cKUnq1KlTjY9thdjYWEn/+cXLXzv+aIbjYwAgkvIRAAIVyRn52muv6Z577tHIkSM1ZsyYYJcDi9Bow1IOh0MpKSnaunVrjY9V3Sa1sk8mS0tLa1JOOccvGT9+CfmJ9u3bp8TERL7NBuATSfkIAIGK1IzMzc3V0KFDlZWVpZkzZ9baeVD3aLRhuUsvvVQ7d+5UXl6epcdt27atJP1mAJ9yyikqLCwst/677777zXMEcvnQqaeeqhYtWmjDhg3ltq1bt07nnHNOlY8FIDJESj4CQHVEWkauXbtWV111lbp27arXX39dDRtyV284odGG5e688041adJEN998s/Lz88tt37lzp5588smAj9u3b181bdpUOTk5vkuzjzvxhyzatm2rjz/+WCUlJb51Cxcu1J49e37zHMefs1hRyFZk4MCB5Y69dOlSffXVV7r22murdAwAkSOS8hEAAhVJGfnFF18oKytLrVu31sKFC7ndMAzxsQks17ZtW82bN0/XX3+9OnTooKFDh6pTp04qKSnRmjVrNH/+fA0fPjzg4zocDj3xxBO6+eabdd555+nGG2/UKaecos8++0yHDx/WnDlzJEk333yz3njjDfXv31/XXXeddu7cqZdfftn3aebJdOnSRZI0btw49evXTw0aNNCgQYMqHX/33Xdr/vz5+t///V+NHz9eBw8e1KOPPqrOnTtrxIgRAc8RQHiLpHz87rvvfI+5OX7lz0MPPSRJSk9P15AhQwKeJ4DwFikZWVRUpH79+unnn3/WpEmT9N5775V7H1wuV8DzRIgJ5kO8Ed6++uorM2rUKNO6dWsTExNjmjZtanr27Gmefvppc+TIEd+49PR0M2zYMN/rWbNmGUlm/fr1FR73X//6lzn//PNNbGyscTgcplu3buaVV17xG/N///d/5tRTTzV2u9307NnTbNiwwfTu3dv07t3bN2bXrl1Gkpk1a5Zv3bFjx8ztt99uWrRoYWw2m6nKfyJbt241ffv2NY0bNzYJCQlm8ODBxu12V+1NAhCRIiEfly9fbiRVuJx4LgD4tXDPyOP7V7acOCfUXzZjLHh4HAAAAAAAkMQ92gAAAAAAWIpGGwAAAAAAC9FoAwAAAABgIRptAAAAAAAsRKMNAAAAAICF6uVztMvKyrR37141bdpUNpst2OUAqIeMMSoqKlJKSoqiosLnM0fyEUBNhWs+SmQkgJoJJB/rZaO9d+9epaamBrsMAGFgz549atWqVbDLsAz5CMAq4ZaPEhkJwBpVycd62Wg3bdpU0n8m6HA4glwNgPrI6/UqNTXVlyfhgnwEUFPhmo8SGQmgZgLJx3rZaB+/1MfhcBCSAGok3C4dJB8BWCXc8lEiIwFYoyr5GF433gAAAAAAEGQ02gAAAAAAWIhGGwAAAAAAC9FoAwAAAABgIRptAAAAAAAsRKMNAAAAAICFaLQBAAAAALBQjRrtRx55RDabTRMmTPCtO3LkiLKzs9WsWTPFxcVp4MCBys/P99tv9+7dysrKUuPGjZWUlKRJkybp2LFjNSkFAEIK+QgAFSMfAUSCajfa69ev1/PPP6+MjAy/9XfccYfeffddzZ8/XytXrtTevXt19dVX+7aXlpYqKytLJSUlWrNmjebMmaPZs2dr6tSp1Z8FAIQQ8hEAKkY+AogUNmOMCXSngwcP6ne/+52effZZPfTQQzrnnHM0Y8YMeTwetWjRQvPmzdM111wjSfryyy/VoUMH5eXlqUePHlq8eLEuvfRS7d27V8nJyZKkmTNnavLkydq/f79iYmJ+8/xer1fx8fHyeDxyOBy/Od4Yo1+OlgY6TQD1TGx0A9lstiqNDTRHqop8BBCKyMfA50Y+ApGjqhkZSIY0rE4h2dnZysrKUmZmph566CHf+o0bN+ro0aPKzMz0rWvfvr3S0tJ8QZmXl6fOnTv7QlKS+vXrpzFjxmjbtm0699xzy52vuLhYxcXFfhMMxC9HS9Vx6pKA9gFQ/3z+YD81jqlWrFmGfAQQiiIxH6WaZST5CESO2sjIgI/26quv6pNPPtH69evLbXO73YqJiVFCQoLf+uTkZLndbt+YE0Py+Pbj2yqSk5OjBx54INBSAaBOkY8AULFg5KNERgIInoAa7T179mj8+PHKzc1Vo0aNaqumcqZMmaKJEyf6Xnu9XqWmplZ5/9joBvr8wX61URqAEBIb3SBo5yYfAYSySMxHqWYZST4CkaM2MjKgRnvjxo0qKCjQ7373O9+60tJSrVq1Sn/961+1ZMkSlZSUqLCw0O9Tyfz8fDmdTkmS0+nUunXr/I57/Fclj4/5NbvdLrvdHkipfmw2W9AvlwIQ3shHAKhYsPJRqllGko8AaiKgXx3v06ePtmzZok2bNvmWrl27avDgwb6/o6OjtXTpUt8+27dv1+7du+VyuSRJLpdLW7ZsUUFBgW9Mbm6uHA6HOnbsaNG0AKBukY8AUDHyEUAkCuhjuqZNm6pTp05+65o0aaJmzZr51o8cOVITJ05UYmKiHA6Hbr/9drlcLvXo0UOS1LdvX3Xs2FFDhgzR9OnT5Xa7dc899yg7O7tG38oAQDCRjwBQMfIRQCSy/HqYJ554QlFRURo4cKCKi4vVr18/Pfvss77tDRo00MKFCzVmzBi5XC41adJEw4YN04MPPmh1KQAQUshHAKgY+Qgg3FTrOdrBVlvPdwQQOcI1R8J1XgDqTjjnSDjPDUDtCyRDArpHGwAAAAAAnByNNgAAAAAAFqLRBgAAAADAQjTaAAAAAABYiEYbAAAAAAAL0WgDAAAAAGAhGm0AAAAAACxEow0AAAAAgIVotAEAAAAAsBCNNgAAAAAAFqLRBgAAAADAQjTaAAAAAABYiEYbAAAAAAAL0WgDAAAAAGAhGm0AAAAAACxEow0AAAAAgIVotAEAAAAAsBCNNgAAAAAAFqLRBgAAAADAQjTaAAAAAABYiEYbAAAAAAAL0WgDAAAAAGAhGm0AAAAAACxEow0AAAAAgIUCarSfe+45ZWRkyOFwyOFwyOVyafHixb7tF154oWw2m99y6623+h1j9+7dysrKUuPGjZWUlKRJkybp2LFj1swGAIKEfASAypGRACJNw0AGt2rVSo888ohOP/10GWM0Z84cXXHFFfr000911llnSZJGjRqlBx980LdP48aNfX+XlpYqKytLTqdTa9as0b59+zR06FBFR0dr2rRpFk0JAOoe+QgAlSMjAUQcU0OnnHKKefHFF40xxvTu3duMHz++0rGLFi0yUVFRxu12+9Y999xzxuFwmOLi4iqf0+PxGEnG4/FUu24Aka0ucoR8BFAf1VWOkJEA6ptAMqTa92iXlpbq1Vdf1aFDh+RyuXzr586dq+bNm6tTp06aMmWKDh8+7NuWl5enzp07Kzk52beuX79+8nq92rZtW6XnKi4ultfr9VsAIFSRjwBQOTISQCQI6NJxSdqyZYtcLpeOHDmiuLg4vf322+rYsaMk6cYbb1R6erpSUlK0efNmTZ48Wdu3b9dbb70lSXK73X4BKcn32u12V3rOnJwcPfDAA4GWCgB1inwEgMqRkQAiScCN9plnnqlNmzbJ4/HojTfe0LBhw7Ry5Up17NhRo0eP9o3r3LmzWrZsqT59+mjnzp1q27ZttYucMmWKJk6c6Hvt9XqVmppa7eMBQG0gHwGgcmQkgEgS8KXjMTExateunbp06aKcnBydffbZevLJJysc2717d0nSjh07JElOp1P5+fl+Y46/djqdlZ7Tbrf7fqXy+AIAoYZ8BIDKkZEAIkmNn6NdVlam4uLiCrdt2rRJktSyZUtJksvl0pYtW1RQUOAbk5ubK4fD4bt0CADCBfkIAJUjIwGEs4AuHZ8yZYoGDBigtLQ0FRUVad68eVqxYoWWLFminTt3at68ebrkkkvUrFkzbd68WXfccYd69eqljIwMSVLfvn3VsWNHDRkyRNOnT5fb7dY999yj7Oxs2e32WpkgANQF8hEAKkdGAog0ATXaBQUFGjp0qPbt26f4+HhlZGRoyZIluvjii7Vnzx79+9//1owZM3To0CGlpqZq4MCBuueee3z7N2jQQAsXLtSYMWPkcrnUpEkTDRs2zO+ZiQBQH5GPAFA5MhJApLEZY0ywiwiU1+tVfHy8PB4P99oAqJZwzZFwnReAuhPOORLOcwNQ+wLJkBrfow0AAAAAAP6LRhsAAAAAAAvRaAMAAAAAYCEabQAAAAAALESjDQAAAACAhWi0AQAAAACwEI02AAAAAAAWotEGAAAAAMBCNNoAAAAAAFiIRhsAAAAAAAvRaAMAAAAAYCEabQAAAAAALESjDQAAAACAhWi0AQAAAACwEI02AAAAAAAWotEGAAAAAMBCNNoAAAAAAFiIRhsAAAAAAAvRaAMAAAAAYCEabQAAAAAALESjDQAAAACAhWi0AQAAAACwEI02AAAAAAAWotEGAAAAAMBCATXazz33nDIyMuRwOORwOORyubR48WLf9iNHjig7O1vNmjVTXFycBg4cqPz8fL9j7N69W1lZWWrcuLGSkpI0adIkHTt2zJrZAECQkI8AUDkyEkCkCajRbtWqlR555BFt3LhRGzZs0EUXXaQrrrhC27ZtkyTdcccdevfddzV//nytXLlSe/fu1dVXX+3bv7S0VFlZWSopKdGaNWs0Z84czZ49W1OnTrV2VgBQx8hHAKgcGQkg4pgaOuWUU8yLL75oCgsLTXR0tJk/f75v2xdffGEkmby8PGOMMYsWLTJRUVHG7Xb7xjz33HPG4XCY4uLiKp/T4/EYScbj8dS0fAARqi5yhHwEUB/VVY6QkQDqm0AypNr3aJeWlurVV1/VoUOH5HK5tHHjRh09elSZmZm+Me3bt1daWpry8vIkSXl5eercubOSk5N9Y/r16yev1+v7RLMixcXF8nq9fgsAhCryEQAqR0YCiAQBN9pbtmxRXFyc7Ha7br31Vr399tvq2LGj3G63YmJilJCQ4Dc+OTlZbrdbkuR2u/0C8vj249sqk5OTo/j4eN+SmpoaaNkAUOvIRwCoHBkJIJIE3GifeeaZ2rRpk9auXasxY8Zo2LBh+vzzz2ujNp8pU6bI4/H4lj179tTq+QCgOshHAKgcGQkgkjQMdIeYmBi1a9dOktSlSxetX79eTz75pK6//nqVlJSosLDQ7xPJ/Px8OZ1OSZLT6dS6dev8jnf8FyWPj6mI3W6X3W4PtFQAqFPkIwBUjowEEElq/BztsrIyFRcXq0uXLoqOjtbSpUt927Zv367du3fL5XJJklwul7Zs2aKCggLfmNzcXDkcDnXs2LGmpQBASCEfAaByZCSAcBbQN9pTpkzRgAEDlJaWpqKiIs2bN08rVqzQkiVLFB8fr5EjR2rixIlKTEyUw+HQ7bffLpfLpR49ekiS+vbtq44dO2rIkCGaPn263G637rnnHmVnZ/NpI4B6jXwEgMqRkQAiTUCNdkFBgYYOHap9+/YpPj5eGRkZWrJkiS6++GJJ0hNPPKGoqCgNHDhQxcXF6tevn5599lnf/g0aNNDChQs1ZswYuVwuNWnSRMOGDdODDz5o7awAoI6RjwBQOTISQKSxGWNMsIsIlNfrVXx8vDwejxwOR7DLAVAPhWuOhOu8ANSdcM6RcJ4bgNoXSIbU+B5tAAAAAADwXzTaAAAAAABYiEYbAAAAAAAL0WgDAAAAAGAhGm0AAAAAACxEow0AAAAAgIVotAEAAAAAsBCNNgAAAAAAFqLRBgAAAADAQjTaAAAAAABYiEYbAAAAAAAL0WgDAAAAAGAhGm0AAAAAACxEow0AAAAAgIVotAEAAAAAsBCNNgAAAAAAFqLRBgAAAADAQjTaAAAAAABYiEYbAAAAAAAL0WgDAAAAAGAhGm0AAAAAACxEow0AAAAAgIVotAEAAAAAsBCNNgAAAAAAFgqo0c7JydF5552npk2bKikpSVdeeaW2b9/uN+bCCy+UzWbzW2699Va/Mbt371ZWVpYaN26spKQkTZo0SceOHav5bAAgSMhHAKgcGQkg0jQMZPDKlSuVnZ2t8847T8eOHdPdd9+tvn376vPPP1eTJk1840aNGqUHH3zQ97px48a+v0tLS5WVlSWn06k1a9Zo3759Gjp0qKKjozVt2jQLpgQAdY98BIDKkZEAIo3NGGOqu/P+/fuVlJSklStXqlevXpL+82nkOeecoxkzZlS4z+LFi3XppZdq7969Sk5OliTNnDlTkydP1v79+xUTE/Ob5/V6vYqPj5fH45HD4ahu+QAiWG3nCPkIoL6qixwhIwHUR4FkSI3u0fZ4PJKkxMREv/Vz585V8+bN1alTJ02ZMkWHDx/2bcvLy1Pnzp19ASlJ/fr1k9fr1bZt2yo8T3Fxsbxer98CAKGMfASAypGRAMJdQJeOn6isrEwTJkxQz5491alTJ9/6G2+8Uenp6UpJSdHmzZs1efJkbd++XW+99ZYkye12+wWkJN9rt9td4blycnL0wAMPVLdUAKhT5CMAVI6MBBAJqt1oZ2dna+vWrVq9erXf+tGjR/v+7ty5s1q2bKk+ffpo586datu2bbXONWXKFE2cONH32uv1KjU1tXqFA0AtIx8BoHJkJIBIUK1Lx8eOHauFCxdq+fLlatWq1UnHdu/eXZK0Y8cOSZLT6VR+fr7fmOOvnU5nhcew2+1yOBx+CwCEIvIRACpHRgKIFAE12sYYjR07Vm+//baWLVumNm3a/OY+mzZtkiS1bNlSkuRyubRlyxYVFBT4xuTm5srhcKhjx46BlAMAIYN8BIDKkZEAIk1Al45nZ2dr3rx5euedd9S0aVPf/TDx8fGKjY3Vzp07NW/ePF1yySVq1qyZNm/erDvuuEO9evVSRkaGJKlv377q2LGjhgwZounTp8vtduuee+5Rdna27Ha79TMEgDpAPgJA5chIABHHBEBShcusWbOMMcbs3r3b9OrVyyQmJhq73W7atWtnJk2aZDwej99xvv32WzNgwAATGxtrmjdvbv7whz+Yo0ePVrkOj8djJJU7LgBUldU5Qj4CCBe1kSNkJIBwEEiG1Og52sHCMxAB1FS45ki4zgtA3QnnHAnnuQGofXX2HG0AAAAAAOCPRhsAAAAAAAvRaAMAAAAAYCEabQAAAAAALESjDQAAAACAhWi0AQAAAACwEI02AAAAAAAWotEGAAAAAMBCNNoAAAAAAFiIRhsAAAAAAAvRaAMAAAAAYCEabQAAAAAALESjDQAAAACAhWi0AQAAAACwEI02AAAAAAAWotEGAAAAAMBCNNoAAAAAAFiIRhsAAAAAAAvRaAMAAAAAYCEabQAAAAAALESjDQAAAACAhWi0AQAAAACwEI02AAAAAAAWotEGAAAAAMBCATXaOTk5Ou+889S0aVMlJSXpyiuv1Pbt2/3GHDlyRNnZ2WrWrJni4uI0cOBA5efn+43ZvXu3srKy1LhxYyUlJWnSpEk6duxYzWcDAEFCPgJA5chIAJEmoEZ75cqVys7O1scff6zc3FwdPXpUffv21aFDh3xj7rjjDr377ruaP3++Vq5cqb179+rqq6/2bS8tLVVWVpZKSkq0Zs0azZkzR7Nnz9bUqVOtmxUA1DHyEQAqR0YCiDQ2Y4yp7s779+9XUlKSVq5cqV69esnj8ahFixaaN2+errnmGknSl19+qQ4dOigvL089evTQ4sWLdemll2rv3r1KTk6WJM2cOVOTJ0/W/v37FRMT85vn9Xq9io+Pl8fjkcPh+M3xpqxMvxw5UN1pAqgnYhslyhZVtc8PA82RQJGPAEJJKOWjVD8yknwEIkdVMzKQDGlYk4I8Ho8kKTExUZK0ceNGHT16VJmZmb4x7du3V1pami8k8/Ly1LlzZ19ASlK/fv00ZswYbdu2Teeee2658xQXF6u4uNhvgoH45cgBdZ//vwHtA6D+WXvtcjVu3DzYZUgiHwGEllDKR6l+ZCT5CESO2sjIav8YWllZmSZMmKCePXuqU6dOkiS3262YmBglJCT4jU1OTpbb7faNOTEgj28/vq0iOTk5io+P9y2pqanVLRsAah35CACVIyMBRIJqf6OdnZ2trVu3avXq1VbWU6EpU6Zo4sSJvtderzegoIxtlKi11y6vjdIAhJDYRonBLkES+Qgg9IRKPkr1JyPJRyBy1EZGVqvRHjt2rBYuXKhVq1apVatWvvVOp1MlJSUqLCz0+0QyPz9fTqfTN2bdunV+xzv+i5LHx/ya3W6X3W6vTqmSJFtUVEhdLgUgfJGPAFC5+pSR5COAmgjo0nFjjMaOHau3335by5YtU5s2bfy2d+nSRdHR0Vq6dKlv3fbt27V79265XC5Jksvl0pYtW1RQUOAbk5ubK4fDoY4dO9ZkLgAQNOQjAFSOjAQQaQL6Rjs7O1vz5s3TO++8o6ZNm/ruh4mPj1dsbKzi4+M1cuRITZw4UYmJiXI4HLr99tvlcrnUo0cPSVLfvn3VsWNHDRkyRNOnT5fb7dY999yj7OzsGn0rAwDBRD4CQOXISAARxwRAUoXLrFmzfGN++eUXc9ttt5lTTjnFNG7c2Fx11VVm3759fsf59ttvzYABA0xsbKxp3ry5+cMf/mCOHj1a5To8Ho+RZDweTyDlA4CP1TlCPgIIF7WRI2QkgHAQSIbU6DnawVIXz3cEEN7CNUfCdV4A6k4450g4zw1A7QskQ6r9eC8AAAAAAFAejTYAAAAAABai0QYAAAAAwEI02gAAAAAAWIhGGwAAAAAAC9FoAwAAAABgIRptAAAAAAAsRKMNAAAAAICFaLQBAAAAALAQjTYAAAAAABai0QYAAAAAwEI02gAAAAAAWIhGGwAAAAAAC9FoAwAAAABgIRptAAAAAAAsRKMNAAAAAICFaLQBAAAAALAQjTYAAAAAABai0QYAAAAAwEI02gAAAAAAWIhGGwAAAAAAC9FoAwAAAABgIRptAAAAAAAsRKMNAAAAAICFAm60V61apcsuu0wpKSmy2WxasGCB3/bhw4fLZrP5Lf379/cbc+DAAQ0ePFgOh0MJCQkaOXKkDh48WKOJAECwkY8AUDHyEUCkCbjRPnTokM4++2w988wzlY7p37+/9u3b51teeeUVv+2DBw/Wtm3blJubq4ULF2rVqlUaPXp04NUDQAghHwGgYuQjgEjTMNAdBgwYoAEDBpx0jN1ul9PprHDbF198offff1/r169X165dJUlPP/20LrnkEj322GNKSUkJtCQACAnkIwBUjHwEEGlq5R7tFStWKCkpSWeeeabGjBmjn376ybctLy9PCQkJvpCUpMzMTEVFRWnt2rUVHq+4uFher9dvAYD6iHwEgIpZnY8SGQkgeCxvtPv3769//vOfWrp0qf7yl79o5cqVGjBggEpLSyVJbrdbSUlJfvs0bNhQiYmJcrvdFR4zJydH8fHxviU1NdXqsgGg1pGPAFCx2shHiYwEEDwBXzr+WwYNGuT7u3PnzsrIyFDbtm21YsUK9enTp1rHnDJliiZOnOh77fV6CUoA9Q75CAAVq418lMhIAMFT64/3Ou2009S8eXPt2LFDkuR0OlVQUOA35tixYzpw4ECl9+XY7XY5HA6/BQDqO/IRACpmRT5KZCSA4Kn1Rvv777/XTz/9pJYtW0qSXC6XCgsLtXHjRt+YZcuWqaysTN27d6/tcgAgZJCPAFAx8hFAfRfwpeMHDx70fbooSbt27dKmTZuUmJioxMREPfDAAxo4cKCcTqd27typO++8U+3atVO/fv0kSR06dFD//v01atQozZw5U0ePHtXYsWM1aNAgfjESQL1GPgJAxchHABHHBGj58uVGUrll2LBh5vDhw6Zv376mRYsWJjo62qSnp5tRo0YZt9vtd4yffvrJ3HDDDSYuLs44HA4zYsQIU1RUVOUaPB6PkWQ8Hk+g5QOAMaZ2coR8BBAOwjUfa2tuACJHIBliM8aYOu3sLeD1ehUfHy+Px8O9NgCqJVxzJFznBaDuhHOOhPPcANS+QDKk1u/RBgAAAAAgktBoAwAAAABgIRptAAAAAAAsRKMNAAAAAICFaLQBAAAAALAQjTYAAAAAABai0QYAAAAAwEI02gAAAAAAWIhGGwAAAAAAC9FoAwAAAABgIRptAAAAAAAsRKMNAAAAAICFaLQBAAAAALAQjTYAAAAAABai0QYAAAAAwEI02gAAAAAAWIhGGwAAAAAAC9FoAwAAAABgIRptAAAAAAAsRKMNAAAAAICFaLQBAAAAALAQjTYAAAAAABai0QYAAAAAwEI02gAAAAAAWCjgRnvVqlW67LLLlJKSIpvNpgULFvhtN8Zo6tSpatmypWJjY5WZmamvv/7ab8yBAwc0ePBgORwOJSQkaOTIkTp48GCNJgIAwUY+AkDFyEcAkSbgRvvQoUM6++yz9cwzz1S4ffr06Xrqqac0c+ZMrV27Vk2aNFG/fv105MgR35jBgwdr27Ztys3N1cKFC7Vq1SqNHj26+rMAgBBAPgJAxchHABHH1IAk8/bbb/tel5WVGafTaR599FHfusLCQmO3280rr7xijDHm888/N5LM+vXrfWMWL15sbDab+eGHH6p0Xo/HYyQZj8dTk/IBRLDazhHyEUB9Fa75aAwZCaBmAskQS+/R3rVrl9xutzIzM33r4uPj1b17d+Xl5UmS8vLylJCQoK5du/rGZGZmKioqSmvXrq3wuMXFxfJ6vX4LANQn5CMAVKy28lEiIwEEj6WNttvtliQlJyf7rU9OTvZtc7vdSkpK8tvesGFDJSYm+sb8Wk5OjuLj431LamqqlWUDQK0jHwGgYrWVjxIZCSB46sWvjk+ZMkUej8e37NmzJ9glAUBIIB8BoHJkJIBgsbTRdjqdkqT8/Hy/9fn5+b5tTqdTBQUFftuPHTumAwcO+Mb8mt1ul8Ph8FsAoD4hHwGgYrWVjxIZCSB4LG2027RpI6fTqaVLl/rWeb1erV27Vi6XS5LkcrlUWFiojRs3+sYsW7ZMZWVl6t69u5XlAEDIIB8BoGLkI4Bw1DDQHQ4ePKgdO3b4Xu/atUubNm1SYmKi0tLSNGHCBD300EM6/fTT1aZNG917771KSUnRlVdeKUnq0KGD+vfvr1GjRmnmzJk6evSoxo4dq0GDBiklJcWyiQFAXSMfAaBi5COAiBPoT5ovX77cSCq3DBs2zBjzn0c03HvvvSY5OdnY7XbTp08fs337dr9j/PTTT+aGG24wcXFxxuFwmBEjRpiioqIq18CjGQDUVG3kCPkIIByEaz4aQ0YCqJlAMsRmjDF12tlbwOv1Kj4+Xh6Ph3ttAFRLuOZIuM4LQN0J5xwJ57kBqH2BZEi9+NVxAAAAAADqCxptAAAAAAAsRKMNAAAAAICFaLQBAAAAALAQjTYAAAAAABai0QYAAAAAwEI02gAAAAAAWIhGGwAAAAAAC9FoAwAAAABgIRptAAAAAAAsRKMNAAAAAICFaLQBAAAAALAQjTYAAAAAABai0QYAAAAAwEI02gAAAAAAWIhGGwAAAAAAC9FoAwAAAABgIRptAAAAAAAsRKMNAAAAAICFaLQBAAAAALAQjTYAAAAAABai0QYAAAAAwEI02gAAAAAAWIhGGwAAAAAAC1neaN9///2y2Wx+S/v27X3bjxw5ouzsbDVr1kxxcXEaOHCg8vPzrS4DAEIO+QgAFSMfAYSbWvlG+6yzztK+fft8y+rVq33b7rjjDr377ruaP3++Vq5cqb179+rqq6+ujTIAIOSQjwBQMfIRQDhpWCsHbdhQTqez3HqPx6O///3vmjdvni666CJJ0qxZs9ShQwd9/PHH6tGjR22UAwAhg3wEgIqRjwDCSa18o/31118rJSVFp512mgYPHqzdu3dLkjZu3KijR48qMzPTN7Z9+/ZKS0tTXl5epccrLi6W1+v1WwCgPiIfAaBiVuejREYCCB7LG+3u3btr9uzZev/99/Xcc89p165d+p//+R8VFRXJ7XYrJiZGCQkJfvskJyfL7XZXesycnBzFx8f7ltTUVKvLBoBaRz4CQMVqIx8lMhJA8Fh+6fiAAQN8f2dkZKh79+5KT0/X66+/rtjY2Godc8qUKZo4caLvtdfrJSgB1DvkIwBUrDbyUSIjAQRPrT/eKyEhQWeccYZ27Nghp9OpkpISFRYW+o3Jz8+v8J6c4+x2uxwOh98CAPUd+QgAFbMiHyUyEkDw1HqjffDgQe3cuVMtW7ZUly5dFB0draVLl/q2b9++Xbt375bL5artUgAgpJCPAFAx8hFAfWf5peN//OMfddlllyk9PV179+7VfffdpwYNGuiGG25QfHy8Ro4cqYkTJyoxMVEOh0O33367XC4XvxgJIOyRjwBQMfIRQLixvNH+/vvvdcMNN+inn35SixYtdMEFF+jjjz9WixYtJElPPPGEoqKiNHDgQBUXF6tfv3569tlnrS4DAEIO+QgAFSMfAYQbmzHGBLuIQHm9XsXHx8vj8XCvDYBqCdccCdd5Aag74Zwj4Tw3ALUvkAyp9Xu0AQAAAACIJDTaAAAAAABYiEYbAAAAAAAL0WgDAAAAAGAhGm0AAAAAACxEow0AAAAAgIVotAEAAAAAsBCNNgAAAAAAFqLRBgAAAADAQjTaAAAAAABYiEYbAAAAAAAL0WgDAAAAAGAhGm0AAAAAACxEow0AAAAAgIVotAEAAAAAsBCNNgAAAAAAFqLRBgAAAADAQjTaAAAAAABYiEYbAAAAAAAL0WgDAAAAAGAhGm0AAAAAACxEow0AAAAAgIVotAEAAAAAsBCNNgAAAAAAFgpqo/3MM8+odevWatSokbp3765169YFsxwACBnkIwBUjHwEUB8ErdF+7bXXNHHiRN1333365JNPdPbZZ6tfv34qKCgIVkkAEBLIRwCoGPkIoL6wGWNMME7cvXt3nXfeefrrX/8qSSorK1Nqaqpuv/123XXXXSfd1+v1Kj4+Xh6PRw6H4zfPVVJSIs/u7yypG0Doik9LV0xMTJXGBpojdYl8BGA18vE/Apkb+QhEjqpmZCAZ0tCq4gJRUlKijRs3asqUKb51UVFRyszMVF5eXrnxxcXFKi4u9r32er0Bnc+z+ztt3t2/+gUDqBcy9L5atDs92GXUCPkIoDZEYj5KNctI8hGIHLWRkUG5dPzHH39UaWmpkpOT/dYnJyfL7XaXG5+Tk6P4+HjfkpqaWlelAkCdIh8BoGKB5qNERgIInqB8ox2oKVOmaOLEib7XXq83oKCMT0tXht6vjdIAhJD4tPRgl1DnyEcAVRGJ+SjVLCPJRyBy1EZGBqXRbt68uRo0aKD8/Hy/9fn5+XI6neXG2+122e32ap8vJiam3l8uBSAykI8AULFA81GqWUaSjwBqIiiXjsfExKhLly5aunSpb11ZWZmWLl0ql8sVjJIAICSQjwBQMfIRQH0StEvHJ06cqGHDhqlr167q1q2bZsyYoUOHDmnEiBHBKgkAQgL5CAAVIx8B1BdBa7Svv/567d+/X1OnTpXb7dY555yj999/v9wPXABApCEfAaBi5COA+iJoz9GuiVB+viOA+iFccyRc5wWg7oRzjoTz3ADUvkAyJCj3aAMAAAAAEK5otAEAAAAAsBCNNgAAAAAAFqLRBgAAAADAQjTaAAAAAABYiEYbAAAAAAALBe052jVx/IlkXq83yJUAqK+O50c9fMLhSZGPAGoqXPNRIiMB1Ewg+VgvG+2ioiJJUmpqapArAVDfFRUVKT4+PthlWIZ8BGCVcMtHiYwEYI2q5KPN1MOPK8vKyrR37141bdpUNputSvt4vV6lpqZqz549v/lw8VAVDnOQmEcoCYc5SNWbhzFGRUVFSklJUVRU+NxFQz7W3zlIzCPUhMM8yEd/gWZkOPwbkMJjHuEwB4l5hJpA5xFIPtbLb7SjoqLUqlWrau3rcDjq9T8GKTzmIDGPUBIOc5ACn0e4fVMjkY/hMAeJeYSacJgH+fgf1c3IcPg3IIXHPMJhDhLzCDWBzKOq+RheH1MCAAAAABBkNNoAAAAAAFgoYhptu92u++67T3a7PdilVFs4zEFiHqEkHOYghc88giUc3r9wmIPEPEJNOMwjHOYQTOHy/oXDPMJhDhLzCDW1OY96+WNoAAAAAACEqoj5RhsAAAAAgLpAow0AAAAAgIVotAEAAAAAsBCNNgAAAAAAFqLRBgAAAADAQhHRaD/zzDNq3bq1GjVqpO7du2vdunXBLumkcnJydN5556lp06ZKSkrSlVdeqe3bt/uNOXLkiLKzs9WsWTPFxcVp4MCBys/PD1LFv+2RRx6RzWbThAkTfOvqyxx++OEH3XTTTWrWrJliY2PVuXNnbdiwwbfdGKOpU6eqZcuWio2NVWZmpr7++usgVlxeaWmp7r33XrVp00axsbFq27at/vznP+vEhw6E4jxWrVqlyy67TCkpKbLZbFqwYIHf9qrUfODAAQ0ePFgOh0MJCQkaOXKkDh48WIezCG3kY/CRj8FFPpKPJ1OfMpJ8DC3kY/CETD6aMPfqq6+amJgY849//MNs27bNjBo1yiQkJJj8/Pxgl1apfv36mVmzZpmtW7eaTZs2mUsuucSkpaWZgwcP+sbceuutJjU11SxdutRs2LDB9OjRw5x//vlBrLpy69atM61btzYZGRlm/PjxvvX1YQ4HDhww6enpZvjw4Wbt2rXmm2++MUuWLDE7duzwjXnkkUdMfHy8WbBggfnss8/M5Zdfbtq0aWN++eWXIFbu7+GHHzbNmjUzCxcuNLt27TLz5883cXFx5sknn/SNCcV5LFq0yPzpT38yb731lpFk3n77bb/tVam5f//+5uyzzzYff/yx+fDDD027du3MDTfcUMczCU3kY/CRj8FHPpKPlalvGUk+hg7ykXw0xpiwb7S7detmsrOzfa9LS0tNSkqKycnJCWJVgSkoKDCSzMqVK40xxhQWFpro6Ggzf/5835gvvvjCSDJ5eXnBKrNCRUVF5vTTTze5ubmmd+/evqCsL3OYPHmyueCCCyrdXlZWZpxOp3n00Ud96woLC43dbjevvPJKXZRYJVlZWeb3v/+937qrr77aDB482BhTP+bx66CsSs2ff/65kWTWr1/vG7N48WJjs9nMDz/8UGe1hyryMbjIx9BAPpKPlanvGUk+Bg/5GDrzCGY+hvWl4yUlJdq4caMyMzN966KiopSZmam8vLwgVhYYj8cjSUpMTJQkbdy4UUePHvWbV/v27ZWWlhZy88rOzlZWVpZfrVL9mcO//vUvde3aVddee62SkpJ07rnn6m9/+5tv+65du+R2u/3mER8fr+7du4fUPM4//3wtXbpUX331lSTps88+0+rVqzVgwABJ9WceJ6pKzXl5eUpISFDXrl19YzIzMxUVFaW1a9fWec2hhHwMPvIxNJCP5GNFwiEjycfgIR9Dax4nqst8bGhd2aHnxx9/VGlpqZKTk/3WJycn68svvwxSVYEpKyvThAkT1LNnT3Xq1EmS5Ha7FRMTo4SEBL+xycnJcrvdQaiyYq+++qo++eQTrV+/vty2+jKHb775Rs8995wmTpyou+++W+vXr9e4ceMUExOjYcOG+Wqt6N9YKM3jrrvuktfrVfv27dWgQQOVlpbq4Ycf1uDBgyWp3szjRFWp2e12KykpyW97w4YNlZiYGLLzqivkY3CRj6EzD/Lxv8jH/6rvGUk+Bhf5GFrzOFFd5mNYN9rhIDs7W1u3btXq1auDXUpA9uzZo/Hjxys3N1eNGjUKdjnVVlZWpq5du2ratGmSpHPPPVdbt27VzJkzNWzYsCBXV3Wvv/665s6dq3nz5umss87Spk2bNGHCBKWkpNSreQAnIh+Di3wEQhf5GFzkI6Qw/9Xx5s2bq0GDBuV+iTA/P19OpzNIVVXd2LFjtXDhQi1fvlytWrXyrXc6nSopKVFhYaHf+FCa18aNG1VQUKDf/e53atiwoRo2bKiVK1fqqaeeUsOGDZWcnBzyc5Ckli1bqmPHjn7rOnTooN27d0uSr9ZQ/zc2adIk3XXXXRo0aJA6d+6sIUOG6I477lBOTo6k+jOPE1WlZqfTqYKCAr/tx44d04EDB0J2XnWFfAwe8jG05kE+/hf5+F/1OSPJx+AjH0NrHieqy3wM60Y7JiZGXbp00dKlS33rysrKtHTpUrlcriBWdnLGGI0dO1Zvv/22li1bpjZt2vht79Kli6Kjo/3mtX37du3evTtk5tWnTx9t2bJFmzZt8i1du3bV4MGDfX+H+hwkqWfPnuUejfHVV18pPT1dktSmTRs5nU6/eXi9Xq1duzak5nH48GFFRfn/596gQQOVlZVJqj/zOFFVana5XCosLNTGjRt9Y5YtW6aysjJ17969zmsOJeRj8JCPoZUr5CP5WJH6mJHkY2jMQSIfQ20eJ6rTfKzhD7mFvFdffdXY7XYze/Zs8/nnn5vRo0ebhIQE43a7g11apcaMGWPi4+PNihUrzL59+3zL4cOHfWNuvfVWk5aWZpYtW2Y2bNhgXC6XcblcQaz6t534q5HG1I85rFu3zjRs2NA8/PDD5uuvvzZz5841jRs3Ni+//LJvzCOPPGISEhLMO++8YzZv3myuuOKKoD/W4NeGDRtmTj31VN/jGd566y3TvHlzc+edd/rGhOI8ioqKzKeffmo+/fRTI8k8/vjj5tNPPzXfffddlWvu37+/Offcc83atWvN6tWrzemnn87ja/4/8jF0kI/BQz6Sj5WpbxlJPoYO8pF8NCYCHu9ljDFPP/20SUtLMzExMaZbt27m448/DnZJJyWpwmXWrFm+Mb/88ou57bbbzCmnnGIaN25srrrqKrNv377gFV0Fvw7K+jKHd99913Tq1MnY7XbTvn1788ILL/htLysrM/fee69JTk42drvd9OnTx2zfvj1I1VbM6/Wa8ePHm7S0NNOoUSNz2mmnmT/96U+muLjYNyYU57F8+fIK/1sYNmxYlWv+6aefzA033GDi4uKMw+EwI0aMMEVFRUGYTWgiH0MD+Rg85CP5eDL1KSPJx9BCPgZPqOSjzRhjqv79NwAAAAAAOJmwvkcbAAAAAIC6RqMNAAAAAICFaLQBAAAAALAQjTYAAAAAABai0QYAAAAAwEI02gAAAAAAWIhGGwAAAAAAC9FoAwAAAABgIRptAAAAAAAsRKMNAAAAAICFaLQBAAAAALDQ/wMRi6A9MsIVcAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig = plt.figure(figsize=(4*n_circuits, 4))\n", + "for idx_circuit in range(n_circuits):\n", + " ax = fig.add_subplot(1, n_circuits, idx_circuit+1)\n", + " ax.set_title(f'Circuit {idx_circuit}')\n", + " plt.plot(sol_signal.ts[idx_circuit], sol_signal.ys[idx_circuit], label=species)\n", + "plt.suptitle('Jacobian of system')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def compute_adaptability_full(ts):\n", + " \"\"\" ts: time series with dimensions [t, species] \"\"\"\n", + " \n", + " x0 = ts[0]\n", + " x1 = ts[1]\n", + " \n", + " p = compute_precision(x0, x1)" + ] + } + ], + "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.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/synbio_morpher/scripts/vis_6_scatter/run_vis_6_scatter.py b/synbio_morpher/scripts/vis_6_scatter/run_vis_6_scatter.py index 677b0882..75c88794 100644 --- a/synbio_morpher/scripts/vis_6_scatter/run_vis_6_scatter.py +++ b/synbio_morpher/scripts/vis_6_scatter/run_vis_6_scatter.py @@ -17,7 +17,7 @@ from synbio_morpher.utils.misc.string_handling import prettify_keys_for_label from synbio_morpher.utils.misc.scripts_io import get_search_dir from synbio_morpher.utils.results.analytics.naming import get_true_names_analytics -from synbio_morpher.utils.results.analytics.timeseries import calculate_robustness +from synbio_morpher.utils.results.analytics.timeseries import calculate_adaptation from synbio_morpher.utils.results.experiments import Experiment, Protocol from synbio_morpher.utils.results.result_writer import ResultWriter from synbio_morpher.utils.results.visualisation import visualise_data @@ -60,10 +60,10 @@ def get_selection(m): cols_x = [c for c in cols if 'sensitivity' in c] cols_y = [c for c in cols if 'precision' in c] - data['robustness'] = calculate_robustness( + data['adaptation'] = calculate_adaptation( data[cols_x].to_numpy().squeeze(), data[cols_y].to_numpy().squeeze()) - hue = 'robustness' if ( - ~data['robustness'].isna()).sum() > 0 else 'overshoot' + hue = 'adaptation' if ( + ~data['adaptation'].isna()).sum() > 0 else 'overshoot' for m in list(data['mutation_num'].unique()) + ['all']: data_selected = data @@ -84,7 +84,7 @@ def get_selection(m): cols_x=cols_x, cols_y=cols_y, plot_type='scatter_plot', - out_name=f'robustness_m{m}{extra_naming}{text_log}', + out_name=f'adaptation_m{m}{extra_naming}{text_log}', hue=hue, use_sns=True, log_axis=log_opt, diff --git a/synbio_morpher/utils/results/analytics/timeseries.py b/synbio_morpher/utils/results/analytics/timeseries.py index 81be76f5..8437ae8a 100644 --- a/synbio_morpher/utils/results/analytics/timeseries.py +++ b/synbio_morpher/utils/results/analytics/timeseries.py @@ -90,8 +90,9 @@ def compute_sensitivity_simple(starting_states, peaks, signal_factor): numer, signal_factor)) # type: ignore -def calculate_robustness(s, p): - """ s = sensitivity, p = precision """ +def calculate_adaptation(s, p): + """ Adaptation = robustness to noise + s = sensitivity, p = precision """ return np.log(log_distance(s=s, p=p) * np.log(sp_prod( s=s, p=p, sp_factor=(p / s).max(), s_weight=(np.log(p) / s))))