Merge pull request #105 from ucl-cssb/CRM

CRM update

Author: Fontanapink

examples/CRM/examples-bayes-CRM.ipynb

@@ -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:

examples/CRM/examples-sim-CRM.ipynb

@@ -3,6 +3,7 @@
 
 import random
 import numpy as np
+import pandas as pd
 import matplotlib.pyplot as plt
 
 cols = ["red", "green", "blue", "royalblue", "orange", "black"]
@@ -45,8 +46,7 @@ def plot_CRM(observed_species, observed_resources, timepoints, csv_file=None):
         # Use a different color index for resources (continuing from where
         # species left off)
         color_idx = observed_species.shape[1] + resource_idx
-        # Use modulo to cycle through colors if we have more entities than
-        # colors
+
         color_idx = color_idx % len(cols)
 
         label = f'Resource {resource_idx + 1}'
@@ -109,6 +109,85 @@ 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):
+    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)
+
+    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
+    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))
+
+    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]))
+
+    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}')
+
+        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.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.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)
+    plt.show()
+
+
 def plot_gMLV(yobs, sobs, timepoints):
     # fig, axs = plt.subplots(1, 2, layout='constrained')
     fig, axs = plt.subplots(1, 2)