debug simulation

dynamo/simulation/Gillespie.py

@@ -60,9 +60,10 @@ def Gillespie(
     """
 
     gene_num, species_num = C0.shape[0:2]
-    adata_no_splicing, P = None, None
+    adata_no_splicing, P, layers_no_splicing = None, None, None
 
     if method == "basic":
+        gene_num = 2
         if t_eval is None:
             steps = (t_span[1] - t_span[0]) // dt  # // int; %% remainder
             t_eval = np.linspace(t_span[0], t_span[1], steps)
@@ -75,10 +76,10 @@ def Gillespie(
             si=si,
             be=be,
             ga=ga,
-            C0=np.zeros((5, 1)),
+            C0=[np.zeros((1, 5))] * 2,
             t_span=[0, 50],
             n_traj=1,
-            t_eval=None,
+            t_eval=t_eval,
             report=verbose,
         )  # unfinished, no need to interpolate now.
         uu, ul, su, sl = [np.transpose(trajs_C[:, :, i + 1, :].reshape((gene_num, -1))) for i in range(4)]
@@ -108,20 +109,22 @@ def Gillespie(
             }
         )
         obs.set_index("Cell_name", inplace=True)
+        true_beta, true_gamma, delta, eta = be, ga, None, None
     elif method == "simulate_2bifurgenes":
         gene_num = 2
+        param_dict = {
+            "a": [20, 20],
+            "b": [20, 20],
+            "S": [20, 20],
+            "K": [20, 20],
+            "m": [3, 3],
+            "n": [3, 3],
+            "beta": [1, 1],
+            "gamma": [1, 1],
+        }
         _, trajs_C = simulate_2bifurgenes(
-            a1=20,
-            b1=20,
-            a2=20,
-            b2=20,
-            K=20,
-            n=3,
-            be1=1,
-            ga1=1,
-            be2=1,
-            ga2=1,
             C0=np.zeros(4),
+            param_dict=param_dict,
             t_span=t_span,
             n_traj=n_traj,
             report=verbose,
@@ -150,6 +153,7 @@ def Gillespie(
             }
         )
         obs.set_index("Cell_name", inplace=True)
+        true_beta, true_gamma, delta, eta = param_dict["beta"], param_dict["gamma"], None, None
     elif method == "differentiation":
         gene_num = 2
 
@@ -413,6 +417,7 @@ def Gillespie(
             }
         )
         obs.set_index("cell_name", inplace=True)
+        true_beta, true_gamma = [beta, beta], [gamma, gamma]
     elif method == "oscillation":
         gene_num = 2
 
@@ -658,6 +663,7 @@ def Gillespie(
             }
         )
         obs.set_index("cell_name", inplace=True)
+        true_beta, true_gamma = [beta, beta], [gamma, gamma]
     else:
         raise Exception("method not implemented!")
     # anadata: observation x variable (cells x genes)
@@ -668,8 +674,8 @@ def Gillespie(
     var = pd.DataFrame(
         {
             "gene_short_name": ["gene_%d" % (i) for i in range(gene_num)],
-            "true_beta": [beta, beta],
-            "true_gamma": [gamma, gamma],
+            "true_beta": true_beta,
+            "true_gamma": true_gamma,
             "true_eta": [eta, eta],
             "true_delta": [delta, delta],
         }
@@ -682,17 +688,21 @@ def Gillespie(
         var.copy(),
         layers=layers.copy(),
     )
-    adata_no_splicing = anndata.AnnData(
-        scipy.sparse.csc_matrix(E.astype(int)).copy(),
-        obs.copy(),
-        var.copy(),
-        layers=layers_no_splicing.copy(),
-    )
     if P is not None:
         adata.obsm["protein"] = P
-        adata_no_splicing.obsm["protein"] = P
     # remove cells that has no expression
     adata = adata[np.array(adata.X.sum(1)).flatten() > 0, :]
-    adata_no_splicing = adata_no_splicing[np.array(adata_no_splicing.X.sum(1)).flatten() > 0, :]
+
+    if layers_no_splicing is not None:
+        adata_no_splicing = anndata.AnnData(
+            scipy.sparse.csc_matrix(E.astype(int)).copy(),
+            obs.copy(),
+            var.copy(),
+            layers=layers_no_splicing.copy(),
+        )
+        if P is not None:
+            adata_no_splicing.obsm["protein"] = P
+        # remove cells that has no expression
+        adata_no_splicing = adata_no_splicing[np.array(adata_no_splicing.X.sum(1)).flatten() > 0, :]
 
     return adata, adata_no_splicing

dynamo/simulation/ODE.py

@@ -80,7 +80,7 @@ def hill_act_grad(x: float, A: float, K: float, n: float, g: float = 0) -> float
 
 
 def toggle(
-    ab: Union[np.ndarray, Tuple[float, float]], t: Optional[float] = None, beta: float = 5, gamma: float = 1, n: int = 2
+    ab: Union[np.ndarray, Tuple[float, float]], beta: float = 5, gamma: float = 1, n: int = 2
 ) -> np.ndarray:
     """Calculates the right-hand side (RHS) of the differential equations for the toggle switch system.
 
@@ -105,7 +105,7 @@ def toggle(
     return res
 
 
-def Ying_model(x: np.ndarray, t: Optional[float]=None):
+def Ying_model(x: np.ndarray):
     """network used in the potential landscape paper from Ying, et. al:
     https://www.nature.com/articles/s41598-017-15889-2.
     This is also the mixture of Gaussian model.
@@ -124,7 +124,7 @@ def Ying_model(x: np.ndarray, t: Optional[float]=None):
     return ret
 
 
-def jacobian_Ying_model(x, t=None):
+def jacobian_Ying_model(x):
     """network used in the potential landscape paper from Ying, et. al:
     https://www.nature.com/articles/s41598-017-15889-2.
     This is also the mixture of Gaussian model.
@@ -333,14 +333,12 @@ def ode_neurongenesis(
 
 def neurongenesis(
     x,
-    t=None,
     mature_mu=0,
     n=4,
     k=1,
     a=4,
     eta=0.25,
     eta_m=0.125,
-    eta_b=0.1,
     a_s=2.2,
     a_e=6,
     mx=10,
@@ -386,14 +384,10 @@ def neurongenesis(
     return dx
 
 
-def hsc():
-    pass
-
-
 def state_space_sampler(ode, dim, seed_num=19491001, clip=True, min_val=0, max_val=4, N=10000):
     """Sample N points from the dim dimension gene expression space while restricting the values to be between min_val and max_val. Velocity vector at the sampled points will be calculated according to ode function."""
 
-    seed(seed)
+    seed(seed_num)
     X = np.array([[uniform(min_val, max_val) for _ in range(dim)] for _ in range(N)])
     Y = np.clip(X + ode(X), a_min=min_val, a_max=None) if clip else X + ode(X)