CellGraph

dycpm/analysis.py

@@ -58,7 +58,31 @@ def neighbor_index(state, same, other):
             values.append(count / total)
     return np.mean(values)
 
-def cluster_number(state, type_id):
-    """Calculate number of clusters of a given cell type"""
-    for cell_id in cells.keys():
-        pass
+class CellGraph(nx.Graph):
+    """Graph of cell connectivity"""
+
+    def __init__(self, state, type_id):
+        """Construct graph of a given cell type from a state"""
+        super().__init__()
+        lattice = state.lattice()
+
+        visited_epi = []
+        for cell_id, cell in state.cells.items():
+            if cell.type_id == type_id:
+                self.add_node(cell_id)
+                visited_epi.append(cell_id)
+
+                visited_cached = visited_epi.copy()
+                for vid in cell.voxel_ids:
+                    ref = state.voxel_coordinates(vid)
+                    for i in range(-1, 2):
+                        for j in range(-1, 2):
+                            for k in range(-1, 2):
+                                xyz = np.array([i, j, k], dtype=int) + ref
+                                if np.any(xyz < 0) or np.any(xyz >= state.side): continue
+
+                                cid = lattice[tuple(xyz)]
+                                if cid in visited_cached: continue
+                                visited_cached.append(cid)
+
+                                if state.cell_type(cid) == type_id: self.add_edge(cell_id, cid)

tests/test_dycpm.py

@@ -3,7 +3,7 @@ import dycpm
 import dycpm.analysis
 import pytest
 import json
-import numpy as np
+import numpy as np, networkx as nx
 
 def test_cell():
     """Test `dycpm.Cell`"""
@@ -213,3 +213,50 @@ def test_analysis_neighborIndex():
     assert dycpm.analysis.neighbor_index(state, [1], [0, 2]) == 2 / 3
     assert dycpm.analysis.neighbor_index(state, [2], [1]) == 1 / 3
     assert dycpm.analysis.neighbor_index(state, [2], [0, 1]) == 2 / 3
+
+def test_analysis_cellGraph():
+    """ Test `dycpm.analysis.cellGraph()`"""
+    state = dycpm.State(
+        time=0, side=10, action_probabilities=3 * [0.3], death_probability=1e-4,
+        growth_rate=2 * [1.6], volume_elasticity=2 * [1], surface_tension=5 * [10],
+        division_distribution=dycpm.DIVISION_VOLUME_DISTRIBUTION,
+        cells={
+            1: {
+                'type_id': 1, 'parent_id': 0, 'age': 0,
+                'preferred_volume': 1,
+                'division_volume': dycpm.sample_division_volume(),
+                "voxel_ids": [4 * 10**2 + 4 * 10 + 4]
+            },
+            2: {
+                'type_id': 1, 'parent_id': 0, 'age': 0,
+                'preferred_volume': 1,
+                'division_volume': dycpm.sample_division_volume(),
+                "voxel_ids": [4 * 10**2 + 4 * 10 + 5]
+            },
+            3: {
+                'type_id': 2, 'parent_id': 0, 'age': 0,
+                'preferred_volume': 1,
+                'division_volume': dycpm.sample_division_volume(),
+                "voxel_ids": [4 * 10**2 + 5 * 10 + 4]
+            },
+            4: {
+                'type_id': 2, 'parent_id': 0, 'age': 0,
+                'preferred_volume': 1,
+                'division_volume': dycpm.sample_division_volume(),
+                "voxel_ids": [0 * 10**2 + 0 * 10 + 0]            
+            },
+            5: {
+                'type_id': 2, 'parent_id': 0, 'age': 0,
+                'preferred_volume': 1,
+                'division_volume': dycpm.sample_division_volume(),
+                "voxel_ids": [1 * 10**2 + 0 * 10 + 0]            
+            }
+        }
+    )
+
+    epi = dycpm.analysis.CellGraph(state, dycpm.EPI)
+    pre = dycpm.analysis.CellGraph(state, dycpm.PRE)
+    assert len(epi.nodes) == 2
+    assert len(pre.nodes) == 3
+    assert nx.number_connected_components(epi) == 1
+    assert nx.number_connected_components(pre) == 2