Analysis module added

dycpm/__init__.py

@@ -167,6 +167,11 @@ class State:
         rem = voxel_id % factor
         return [rem % self.side, rem // self.side, quo]
 
+    def cell_type(self, cell_id):
+        """Return cell type"""
+        if cell_id == 0: return 0
+        return self.cells[cell_id].type_id
+
     def cell_positions(self):
         """Return cell positions"""
         return np.asarray([
@@ -185,3 +190,10 @@ class State:
     def volume(self):
         """Return total volume of all cells"""
         return np.sum([len(cell.voxel_ids) for cell in self.cells.values()])
+
+    def lattice(self):
+        """Return values of the CPM lattice"""
+        values = np.zeros((self.side, ) * 3)
+        for cell_id, cell in self.cells.items():
+            values.flat[cell.voxel_ids] = cell_id
+        return values

dycpm/analysis.py

+"""Tools for analysis of CPM states"""
+import json
+import numpy as np
+from . import State
+
+def neighbor_index(state, same, other):
+    """Calculate the neighbour index for a given cell value
+
+    Arguments:
+        state (dycpm.State): State of the CPM
+        same List[int]: List of type IDs considered as same value
+        other List[int]: List of type IDs considered as the other
+    """
+    lattice = state.lattice()
+
+    values = [] # sample of fractions
+    for cell_id, cell in state.cells.items():
+        if cell.type_id not in same: continue
+        visited = [cell_id]
+        count = 0
+        total = 0
+
+        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: continue
+                        visited.append(cid)
+
+                        type_id = state.cell_type(cid)
+                        if type_id in other:
+                            count += 1
+                        elif type_id not in same:
+                            continue
+                        total += 1
+        if count == 0 and total == 0:
+            values.append(0)
+        else:
+            values.append(count / total)
+    return np.mean(values)

tests/test_dycpm.py

 """Test `dycpm` module"""
 import dycpm
+import dycpm.analysis
 import pytest
 import json
 import numpy as np
@@ -104,3 +105,111 @@ def test_simulation():
 
     assert state.time == 2
     assert state.side == 100
+
+def test_lattice():
+    """Test `State.lattice()"""
+
+    side = 100
+    volume = round(np.mean(dycpm.DIVISION_VOLUME_DISTRIBUTION) / 2)
+    partition = round(np.cbrt(volume));
+    margin = round(np.floor((side - partition) / 2));
+
+    state = dycpm.State(
+        time=0, side=side, 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': volume,
+                'division_volume': dycpm.sample_division_volume(),
+                "voxel_ids": [round(
+                    margin + (n - 1) % partition**2 % partition
+                    + side * (margin + (n - 1) % partition**2 // partition)
+                    + side**2 * (margin + (n - 1) // partition**2)
+                ) for n in range(volume)]
+            }
+        }
+    )
+
+    lattice = state.lattice()
+    assert lattice.shape == (side, ) * 3
+    assert lattice[49, 49, 49] == 1
+    assert lattice[0, 0, 0] == 0
+    assert np.count_nonzero(lattice.ravel() == 1) == volume
+
+def test_analysis_neighborIndex():
+    """ Test `dycpm.analysis.neighbor_index()`"""
+    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]
+            }
+        }
+    )
+    assert dycpm.analysis.neighbor_index(state, [1], [2]) == 0.5
+    assert dycpm.analysis.neighbor_index(state, [1], [0, 2]) == 2 / 3
+    assert dycpm.analysis.neighbor_index(state, [2], [1]) == 1
+    assert dycpm.analysis.neighbor_index(state, [2], [0, 1]) == 1
+
+
+    """ Test `dycpm.analysis`"""
+    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]            
+            }
+        }
+    )
+    assert dycpm.analysis.neighbor_index(state, [1], [2]) == 0.5
+    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