Commit

lineage/Analyze.py

@@ -191,7 +191,7 @@ def Results(tHMMobj: tHMM, LL: float) -> dict[str, Any]:
     results_dict["total_number_of_lineages"] = len(tHMMobj.X)
     results_dict["LL"] = LL
     results_dict["total_number_of_cells"] = sum(
-        [len(lineage.output_lineage) for lineage in tHMMobj.X]
+        [len(lineage) for lineage in tHMMobj.X]
     )
 
     true_states_by_lineage = [

lineage/CellVar.py

@@ -38,23 +38,6 @@ def __init__(self, parent: Optional["CellVar"], state: Optional[int] = None):
         self.right = None
         self.obs = None
 
-    def divide(self, T: np.ndarray, rng=None):
-        """
-        Member function that performs division of a cell.
-        Equivalent to adding another timestep in a Markov process.
-        :param T: The array containing the likelihood of a cell switching states.
-        """
-        rng = np.random.default_rng(rng)
-        # Checking that the inputs are of the right shape
-        assert T.shape[0] == T.shape[1]
-
-        # roll a loaded die according to the row in the transtion matrix
-        left_state, right_state = rng.choice(T.shape[0], size=2, p=T[self.state, :])
-        self.left = CellVar(state=left_state, parent=self)
-        self.right = CellVar(state=right_state, parent=self)
-
-        return self.left, self.right
-
     def isLeafBecauseTerminal(self) -> bool:
         """
         Returns true when a cell is a leaf with no children.
@@ -79,12 +62,6 @@ def isLeaf(self) -> bool:
         # otherwise, it itself is observed and at least one of its daughters is observed
         return False
 
-    def isRootParent(self) -> bool:
-        """
-        Returns true if this cell is the first cell in a lineage.
-        """
-        return self.parent is None
-
 
 @dataclass(init=True, repr=True, eq=True, order=True)
 class Time:

lineage/HMM/E_step.py

@@ -138,10 +138,7 @@ def get_beta(
     )  # MSD of the respective lineage
     ELMSD = EL * MSD
 
-    cIDXs = np.arange(first_leaf)
-    cIDXs = np.flip(cIDXs)
-
-    for pii in cIDXs:
+    for pii in range(first_leaf - 1, -1, -1):
         ch_ii = np.array([pii * 2 + 1, pii * 2 + 2])
         ratt = (beta[ch_ii, :] / MSD_array[ch_ii, :]) @ T.T
         fac1 = np.prod(ratt, axis=0) * ELMSD[pii, :]

lineage/LineageTree.py

@@ -70,14 +70,19 @@ def rand_init(
         )  # roll the dice and yield the state for the first cell
         first_cell = CellVar(parent=None, state=first_state)  # create first cell
         full_lineage = [first_cell]  # instantiate lineage with first cell
+        pIDX = 0
 
-        for cell in full_lineage:  # letting the first cell proliferate
-            if cell.isLeaf():  # if the cell has no daughters...
-                # make daughters by dividing and assigning states
-                full_lineage.extend(cell.divide(T, rng=rng))
+        # fill in the cells
+        while len(full_lineage) < desired_num_cells:  # letting the first cell proliferate
+            parent = full_lineage[pIDX]
+            states = rng.choice(T.shape[0], size=2, p=T[parent.state, :])
 
-            if len(full_lineage) >= desired_num_cells:
-                break
+            full_lineage.append(CellVar(parent=parent, state=states[0]))
+            full_lineage.append(CellVar(parent=parent, state=states[1]))
+            parent.left = full_lineage[-2]
+            parent.right = full_lineage[-1]
+
+            pIDX += 1
 
         # Assign observations
         for i_state in range(pi.size):
@@ -140,7 +145,7 @@ def get_Emission_Likelihoods(X: list[LineageTree], E: list) -> list:
     EL = []
     ii = 0
     for lineageObj in X:  # for each lineage in our Population
-        nl = len(lineageObj.output_lineage)  # getting the lineage length
+        nl = len(lineageObj)  # getting the lineage length
         EL.append(ELstack[ii : (ii + nl), :])  # append the EL_array for each lineage
 
         ii += nl

lineage/figures/figure18.py

@@ -49,7 +49,7 @@
         tmp_lineage = LineageTree.rand_init(
             pi, T, E, desired_num_cells, censor_condition=3, desired_experiment_time=96
         )
-        while len(tmp_lineage.output_lineage) < 3:
+        while len(tmp_lineage) < 3:
             tmp_lineage = LineageTree.rand_init(
                 pi,
                 T,
@@ -59,7 +59,7 @@
                 desired_experiment_time=96,
             )
         population.append(tmp_lineage)
-        nn += len(tmp_lineage.output_lineage)
+        nn += len(tmp_lineage)
 
     # Adding populations into a holder for analysing
     list_of_populations.append(population)

lineage/states/StateDistributionGaPhs.py

@@ -87,8 +87,9 @@ def assign_times(self, full_lineage: list):
         This is used in the creation of LineageTrees
         """
         # traversing the cells by generation
-        for cell in full_lineage:
-            if cell.isRootParent():
+        for ii, cell in enumerate(full_lineage):
+            # handle root separately
+            if ii == 0:
                 cell.time = Time(0, cell.obs[2] + cell.obs[3])
                 cell.time.transition_time = 0 + cell.obs[2]
             else:

lineage/states/StateDistributionGamma.py

@@ -135,8 +135,9 @@ def assign_times(self, full_lineage: list[CellVar]):
         This is used in the creation of LineageTrees.
         """
         # traversing the cells by generation
-        for cell in full_lineage:
-            if cell.isRootParent():
+        for ii, cell in enumerate(full_lineage):
+            # if root
+            if ii == 0:
                 cell.time = Time(0, cell.obs[1])
             else:
                 cell.time = Time(

lineage/states/stateCommon.py

@@ -14,10 +14,9 @@ def basic_censor(cells: list):
     """
     Censors a cell if the cell's parent is censored.
     """
-    for cell in cells:
-        if not cell.isRootParent():
-            if not cell.parent.observed:
-                cell.observed = False
+    for cell in cells[1:]:
+        if not cell.parent.observed:
+            cell.observed = False
 
 
 def bern_estimator(bern_obs: np.ndarray, gammas: np.ndarray):
@@ -114,6 +113,7 @@ def gamma_estimator(
             method="SLSQP",
             constraints=linc,
         )
+        print(res.nit)
 
     assert res.success
     return np.exp(res.x)

lineage/tests/test_StateDistribution.py

@@ -114,10 +114,10 @@ def test_censor(self):
         as expected.
         """
         for lin in self.population:
-            for cell in lin.output_lineage:
-                if not cell.isRootParent:
-                    if not cell.parent.observed:
-                        self.assertFalse(cell.observed)
+            # Skip root parent
+            for cell in lin.output_lineage[1:]:
+                if not cell.parent.observed:
+                    self.assertFalse(cell.observed)
 
 
 @pytest.mark.parametrize("dist", [StateDistribution, StateDistPhase])