Merge branch 'master' into CRM

Author: ChaniaClare

examples/gLV/examples-Rutter-Dekker.ipynb

@@ -90,7 +90,7 @@
     }
    ],
    "source": [
-    "# read in pickled simulated parameters, mu, M, epsilon\n",
+    "# read in pickled simulated parameters, mu, M, epsilon, created in examples-sim-gLV.ipynb\n",
     "num_species = 5\n",
     "with open(\"params-s5.pkl\", \"rb\") as f:\n",
     "    params = pickle.load(f)\n",
@@ -661,20 +661,6 @@
     "#inference.plot_posterior(idata)\n",
     "\n",
     "\n",
-    "\n",
-    "\n",
-    "#nX = num_species\n",
-    "#n_obs = times.shape[0] - 1\n",
-    "#noise_stddev = 0.1\n",
-    "\n",
-    "# Params for shrinkage on M_ij (non diagonal elements)\n",
-    "#DA = nX*nX - nX\n",
-    "#DA0 = 3     # expected number of non zero entries in M_ij\n",
-    "#N = n_obs - 2\n",
-    "\n",
-    "#inference = infergLVbayes(X, F, mu_prior, M_prior, DA=DA, DA0=DA0, N=N, noise_stddev=noise_stddev)\n",
-    "#idata = inference.run_bayes_gLV_shrinkage()\n",
-    "\n",
     "# print summary\n",
     "summary = az.summary(idata, var_names=[\"mu_hat\", \"M_ii_hat\", \"M_ij_hat\", \"M_hat\", \"sigma\"])\n",
     "print(summary[[\"mean\", \"sd\", \"r_hat\"]])\n",
@@ -705,7 +691,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "id": "c6d6c2df",
    "metadata": {
     "ExecuteTime": {
@@ -951,20 +937,6 @@
     "#inference.plot_posterior_pert(idata)\n",
     "\n",
     "\n",
-    "\n",
-    "\n",
-    "#nX = num_species\n",
-    "#n_obs = times.shape[0] - 1\n",
-    "#noise_stddev = 0.1\n",
-    "\n",
-    "# Params for shrinkage on M_ij (non diagonal elements)\n",
-    "#DA = nX*nX - nX\n",
-    "#DA0 = 3     # expected number of non zero entries in M_ij\n",
-    "#N = n_obs - 2\n",
-    "\n",
-    "#inference = infergLVbayes(X, F, mu_prior, M_prior, DA=DA, DA0=DA0, N=N, noise_stddev=noise_stddev, epsilon=epsilon)\n",
-    "#idata = inference.run_bayes_gLV_shrinkage_pert()\n",
-    "\n",
     "# print summary\n",
     "summary = az.summary(idata, var_names=[\"mu_hat\", \"M_ii_hat\", \"M_ij_hat\", \"M_hat\", \"epsilon_hat\", \"sigma\"])\n",
     "print(summary[[\"mean\", \"sd\", \"r_hat\"]])\n",

examples/gLV/examples-bayes-gLV.ipynb

@@ -459,8 +459,10 @@ def run_inference(self) -> None:
         with bayes_model:
             # Priors for unknown model parameters
 
+
             sigma = pm.HalfNormal('sigma', sigma = 0.1, shape=(1,))  # Same sigma for all responses
 
+
             # Conditionally define parameters based on whether priors are
             # provided
 

mimic/model_infer/infer_CRM_bayes.py

@@ -279,7 +279,7 @@ def _run_inference_large(self, **kwargs) -> None:
             c2 = pm.InverseGamma("c2", 2, 8)
             tau = pm.HalfCauchy("tau", beta=tau0)
             lam = pm.HalfCauchy("lam", beta=1, shape=(ndim, ndim))
-            A = pm.Normal('A', mu=A_prior_mu, sigma=tau * lam *
+            A = pm.Normal('A', mu=A_prior_mu, sigma=tau * lam * \
                           at.sqrt(c2 / (c2 + tau**2 * lam**2)), shape=(ndim, ndim))
 
             # If noise covariance is provided, use it as a prior
@@ -438,14 +438,14 @@ def _run_inference_large_xs(self, **kwargs) -> None:
             c2_A = pm.InverseGamma("c2_A", 2, 1)
             tau_A = pm.HalfCauchy("tau_A", beta=tau0_A)
             lam_A = pm.HalfCauchy("lam_A", beta=1, shape=(nX, nX))
-            Ah = pm.Normal('Ah', mu=A_prior_mu, sigma=tau_A * lam_A *
+            Ah = pm.Normal('Ah', mu=A_prior_mu, sigma=tau_A * lam_A * \
                            at.sqrt(c2_A / (c2_A + tau_A**2 * lam_A**2)), shape=(nX, nX))
 
             tau0_B = (DB0 / (DB - DB0)) * 0.1 / np.sqrt(N)
             c2_B = pm.InverseGamma("c2_B", 2, 1)
             tau_B = pm.HalfCauchy("tau_B", beta=tau0_B)
             lam_B = pm.HalfCauchy("lam_B", beta=1, shape=(nS, nX))
-            Bh = pm.Normal('Bh', mu=0, sigma=tau_B * lam_B *
+            Bh = pm.Normal('Bh', mu=0, sigma=tau_B * lam_B * \
                            at.sqrt(c2_B / (c2_B + tau_B**2 * lam_B**2)), shape=(nS, nX))
 
             if noise_cov_prior is not None:

mimic/model_infer/infer_VAR_bayes.py

@@ -10,12 +10,18 @@
 
 
 def plot_gLV(yobs, timepoints):
-    # fig, axs = plt.subplots(1, 2, layout='constrained')
+    # fig, axs = plt.subplots(1, 2, layout='constrained')  # Optional
+    # alternative
     fig, axs = plt.subplots(1, 1)
     for species_idx in range(yobs.shape[1]):
-        axs.plot(timepoints, yobs[:, species_idx], color=cols[species_idx])
+        label = f'Species {species_idx + 1}'  # Add a label for each species
+        axs.plot(timepoints, yobs[:, species_idx],
+                 color=cols[species_idx], label=label)
+
     axs.set_xlabel('time')
     axs.set_ylabel('[species]')
+    axs.legend()  # Ensure the legend is called on the correct axes
+    plt.show()
 
 
 def plot_CRM(observed_species, observed_resources, timepoints, csv_file=None):
@@ -102,49 +108,80 @@ def plot_CRM(observed_species, observed_resources, timepoints, csv_file=None):
 
     return fig, ax
 
-def plot_CRM_with_intervals(observed_species, observed_resources, species_lower, species_upper, 
-                           resource_lower, resource_upper, times, filename=None):
+
+def plot_CRM_with_intervals(
+        observed_species,
+        observed_resources,
+        species_lower,
+        species_upper,
+        resource_lower,
+        resource_upper,
+        times,
+        filename=None):
     fig, ax = plt.subplots(figsize=(12, 8))
-    
+
     # Plot median trajectories
     for i in range(observed_species.shape[1]):
-        ax.plot(times, observed_species[:, i], label=f'Species {i+1}', linewidth=2)
-    
+        ax.plot(times, observed_species[:, i],
+                label=f'Species {i+1}', linewidth=2)
+
     for i in range(observed_resources.shape[1]):
-        ax.plot(times, observed_resources[:, i], label=f'Resource {i+1}', linewidth=2, linestyle='--')
-    
-    # Add confidence ribbons 
+        ax.plot(times,
+                observed_resources[:,
+                                   i],
+                label=f'Resource {i+1}',
+                linewidth=2,
+                linestyle='--')
+
+    # Add confidence ribbons
     for i in range(observed_species.shape[1]):
-        ax.fill_between(times, species_lower[:, i], species_upper[:, i], 
-                       alpha=0.2, color=plt.cm.tab10(i))
-    
+        ax.fill_between(times, species_lower[:, i], species_upper[:, i],
+                        alpha=0.2, color=plt.cm.tab10(i))
+
     for i in range(observed_resources.shape[1]):
-        ax.fill_between(times, resource_lower[:, i], resource_upper[:, i], 
-                       alpha=0.2, color=plt.cm.tab10(i + observed_species.shape[1]))
-    
+        ax.fill_between(times,
+                        resource_lower[:,
+                                       i],
+                        resource_upper[:,
+                                       i],
+                        alpha=0.2,
+                        color=plt.cm.tab10(i + observed_species.shape[1]))
+
     if filename:
         true_data = pd.read_csv(filename)
         true_times = true_data['time'].values
-        
+
         for i in range(observed_species.shape[1]):
             col_name = f'species_{i+1}'
             if col_name in true_data.columns:
-                ax.scatter(true_times, true_data[col_name], 
-                          marker='o', s=30, color=plt.cm.tab10(i), label=f'True {col_name}')
-        
+                ax.scatter(
+                    true_times,
+                    true_data[col_name],
+                    marker='o',
+                    s=30,
+                    color=plt.cm.tab10(i),
+                    label=f'True {col_name}')
+
         for i in range(observed_resources.shape[1]):
             col_name = f'resource_{i+1}'
             if col_name in true_data.columns:
-                ax.scatter(true_times, true_data[col_name], 
-                          marker='s', s=30, color=plt.cm.tab10(i + observed_species.shape[1]), 
-                          label=f'True {col_name}')
-    
+                ax.scatter(
+                    true_times,
+                    true_data[col_name],
+                    marker='s',
+                    s=30,
+                    color=plt.cm.tab10(
+                        i + observed_species.shape[1]),
+                    label=f'True {col_name}')
+
     ax.set_xlabel('Time', fontsize=14)
     ax.set_ylabel('Concentration', fontsize=14)
-    ax.set_title('Consumer-Resource Model Dynamics with 95% Credible Intervals', fontsize=16)
+    ax.set_title(
+        'Consumer-Resource Model Dynamics with 95% Credible Intervals',
+        fontsize=16)
     ax.legend(loc='best', fontsize=12)
     ax.grid(True, alpha=0.3)
-    
+
     plt.tight_layout()
     if filename:
         plt.savefig(f"{filename.split('.')[0]}_with_intervals.png", dpi=300)