running MC

notebooks/23_Monte_Carlo_adaptability_2.ipynb

@@ -7,7 +7,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [
     {
@@ -33,7 +33,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
@@ -42,7 +42,7 @@
        "[cuda(id=0)]"
       ]
      },
-     "execution_count": 13,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -77,7 +77,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -98,7 +98,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -139,7 +139,7 @@
     "total_steps = 10\n",
     "total_samples = 100\n",
     "choose_max = 20\n",
-    "mutation_scale = 0.1\n",
+    "mutation_scale = 1\n",
     "n_circuits_display = 30\n",
     "N0 = 200\n",
     "y00 = np.array([[N0] * len(species_unbound) + [0] * len(species_bound)]).astype(np.float32)\n",
@@ -182,7 +182,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -207,7 +207,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -226,7 +226,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -239,14 +239,14 @@
     "    idxs_next = jnp.argsort(adaptability[:, idxs_output], axis=0)[-choose_max:]\n",
     "    # summary_data = {}\n",
     "    # summary_data['adaptability'] = adaptability\n",
-    "    return params[idxs_next], adaptability\n",
+    "    return params[idxs_next], adaptability, sensitivity, precision\n",
     "\n",
     "\n",
     "def mutate(parents: jnp.ndarray, n_samples_per_parent, mutation_scale):\n",
     "    # Generate mutated samples from each parent\n",
     "    mutated = jax.tree_util.tree_map(\n",
-    "        lambda x: x + mutation_scale * np.random.randn(n_samples_per_parent, *x.shape), parents)\n",
-    "    return mutated\n",
+    "        lambda x: x + x * mutation_scale * np.random.randn(n_samples_per_parent, *x.shape), parents)\n",
+    "    return mutated.reshape(mutated.shape[2], mutated.shape[0] * mutated.shape[1], *mutated.shape[3:])\n",
     "\n",
     "\n",
     "def simulate(y00, reverse_rates, sim_func, t0, t1, tmax, threshold):\n",
@@ -279,19 +279,21 @@
      "text": [
       "\n",
       "\n",
-      "Starting step 1 out of 10\n",
+      "Starting iteration 1 out of 10\n",
       "\n",
-      "\n",
-      "Steady states:  300  iterations.  163  left to steady out.  0:00:20.077289\n",
-      "Steady states:  600  iterations.  81  left to steady out.  0:00:39.409788\n"
+      "\n"
      ]
     }
    ],
    "source": [
     "gen_y = reverse_rates\n",
     "n_samples_per_parent = total_samples//choose_max\n",
-    "params_all = np.zeros((total_steps, total_samples, *reverse_rates.shape))\n",
-    "adaptability_all = np.zeros((total_steps, total_samples))\n",
+    "\n",
+    "# [total_steps, n_outputs, n_samples, n_reactions]\n",
+    "params_all = np.zeros((total_steps, len(species_output), *reverse_rates.shape))\n",
+    "adaptability_all = np.zeros((total_steps, total_samples, len(species)))\n",
+    "sensitivity_all = np.zeros((total_steps, total_samples, len(species)))\n",
+    "precision_all = np.zeros((total_steps, total_samples, len(species)))\n",
     "\n",
     "sim_func = jax.vmap(partial(bioreaction_sim_dfx_expanded,\n",
     "                            t0=t0, t1=t1, dt0=dt0,\n",
@@ -306,8 +308,9 @@
     "\n",
     "    print(f'\\n\\nStarting iteration {step+1} out of {total_steps}\\n\\n')\n",
     "\n",
-    "    ys0, ts0, ys1, ts1 = simulate(y00, gen_y, sim_func, t0, t1, tmax, threshold_steady_state)\n",
-    "    next_starting, adaptability = choose_next(params=gen_y, sol=(ys0, ys1), idxs_signal=idxs_signal, idxs_output=idxs_output,\n",
+    "    ys0, ts0, ys1, ts1 = simulate(\n",
+    "        y00, gen_y, sim_func, t0, t1, tmax, threshold_steady_state)\n",
+    "    next_starting, adaptability, sensitivity, precision = choose_next(params=gen_y, sol=(ys0, ys1), idxs_signal=idxs_signal, idxs_output=idxs_output,\n",
     "                                              use_sensitivity_func1=use_sensitivity_func1, choose_max=choose_max, n_samples_per_parent=n_samples_per_parent)\n",
     "    gen_z = mutate(next_starting, n_samples_per_parent, mutation_scale)\n",
     "    params_all[step] = gen_z\n",