Fix growth, division, and physics invariants

SyMBac/physics/colony.py

@@ -73,6 +73,17 @@ def handle_cell_overlaps(self, newly_born_cells_map: dict[SimCell, SimCell]) ->
                         # If the daughter's head is overlapping with the mother
                         if info.shape in mother_shapes:
 
+                            mother_can_trim = (
+                                len(mother.physics_representation.segments)
+                                > mother.config.MIN_LENGTH_AFTER_DIVISION
+                            )
+                            daughter_can_trim = (
+                                len(daughter.physics_representation.segments)
+                                > daughter.config.MIN_LENGTH_AFTER_DIVISION
+                            )
+                            if not mother_can_trim or not daughter_can_trim:
+                                break
+
                             mother_removed_segment = mother.physics_representation.remove_tail_segment()
                             daughter_removed_segment = daughter.physics_representation.remove_head_segment()
                             if mother_removed_segment is None or daughter_removed_segment is None:

SyMBac/physics/config.py

@@ -103,6 +103,23 @@ class CellConfig:
 
     def __post_init__(self):
 
+        if isinstance(self.GRANULARITY, bool) or not isinstance(self.GRANULARITY, int) or self.GRANULARITY <= 0:
+            raise ValueError("GRANULARITY must be an integer greater than 0.")
+        if not np.isfinite(self.SEPTUM_DURATION) or self.SEPTUM_DURATION <= 0:
+            raise ValueError("SEPTUM_DURATION must be finite and greater than 0.")
+        if (
+            isinstance(self.MIN_LENGTH_AFTER_DIVISION, bool)
+            or not isinstance(self.MIN_LENGTH_AFTER_DIVISION, int)
+            or self.MIN_LENGTH_AFTER_DIVISION < 2
+        ):
+            raise ValueError("MIN_LENGTH_AFTER_DIVISION must be an integer of at least 2.")
+        if (
+            isinstance(self.SEED_CELL_SEGMENTS, bool)
+            or not isinstance(self.SEED_CELL_SEGMENTS, int)
+            or self.SEED_CELL_SEGMENTS < 2
+        ):
+            raise ValueError("SEED_CELL_SEGMENTS must be an integer of at least 2.")
+
         if self.MAX_LENGTH_STD < 0:
             raise ValueError("MAX_LENGTH_STD must be non-negative.")
         if self.WIDTH_STD < 0:
@@ -162,13 +179,17 @@ class PhysicsConfig:
     COLLISION_SLOP: Optional[float] = None # Amount of overlap between shapes that is allowed. To improve stability, set this as high as you can without noticeable overlapping. It defaults to 0.1.
 
     def __post_init__(self):
+        if isinstance(self.ITERATIONS, bool) or not isinstance(self.ITERATIONS, int) or self.ITERATIONS <= 0:
+            raise ValueError("ITERATIONS must be an integer greater than 0.")
+        if not np.isfinite(self.DAMPING) or not 0.0 <= self.DAMPING <= 1.0:
+            raise ValueError("DAMPING must be finite and between 0 and 1.")
         if self.THREADS > 2:
             raise ValueError("THREADS cannot be greater than 2.")
         if not self.THREADED and self.THREADS != 1:
             raise ValueError("If THREADED is False, THREADS must be 1.")
         if self.THREADED and self.THREADS != 2:
             raise ValueError("If THREADED is True, THREADS must be 2.")
-        if self.DT <= 0:
-            raise ValueError("DT must be greater than 0.")
+        if not np.isfinite(self.DT) or self.DT <= 0:
+            raise ValueError("DT must be finite and greater than 0.")
         if self.COLLISION_SLOP is not None and self.COLLISION_SLOP < 0:
             raise ValueError("COLLISION_SLOP must be non-negative when provided.")

SyMBac/physics/growth_manager.py

@@ -21,7 +21,11 @@ def grow(cell: SimCell, dt: float):
             cell: The `SimCell` object to be grown.
             dt: The time step for the current simulation frame.
         """
-        if not cell.is_dividing and cell.length >= cell.max_length: #LENGTH_FIX # TODO: is this actually necessary?
+        enough_segments_to_divide = (
+            cell.physics_representation.num_segments
+            >= 2 * cell.config.MIN_LENGTH_AFTER_DIVISION
+        )
+        if not cell.is_dividing and cell.length >= cell.max_length and enough_segments_to_divide:
             return
 
         if cell.physics_representation.num_segments < 2:
@@ -34,28 +38,21 @@ def grow(cell: SimCell, dt: float):
         cell.physics_representation.growth_accumulator_head += half_growth
         cell.physics_representation.growth_accumulator_tail += half_growth
 
-        # Stretch the head joint by adjusting the anchor on the first segment.
-        # This pushes the head segment outwards.
+        while cell.physics_representation.growth_accumulator_head >= cell.config.GROWTH_THRESHOLD:
+            cell.physics_representation.add_head_segment()
+            cell.physics_representation.growth_accumulator_head -= cell.config.GROWTH_THRESHOLD
+
+        while cell.physics_representation.growth_accumulator_tail >= cell.config.GROWTH_THRESHOLD:
+            cell.physics_representation.add_tail_segment()
+            cell.physics_representation.growth_accumulator_tail -= cell.config.GROWTH_THRESHOLD
+
         first_pivot_joint = cell.physics_representation.pivot_joints[0]
         first_pivot_joint.anchor_a = (
             cell.config.JOINT_DISTANCE / 2 + cell.physics_representation.growth_accumulator_head,
-            0
+            0,
         )
-
-        # Stretch the tail joint by adjusting the anchor on the last segment.
-        # This pushes the tail segment outwards.
         last_pivot_joint = cell.physics_representation.pivot_joints[-1]
         last_pivot_joint.anchor_b = (
             -cell.config.JOINT_DISTANCE / 2 - cell.physics_representation.growth_accumulator_tail,
             0,
         )
-
-        # If the head has grown enough, insert a new segment.
-        if cell.physics_representation.growth_accumulator_head >= cell.config.GROWTH_THRESHOLD:
-            cell.physics_representation.add_head_segment()
-            cell.physics_representation.growth_accumulator_head = 0.0
-
-        # If the tail has grown enough, insert a new segment.
-        if cell.physics_representation.growth_accumulator_tail >= cell.config.GROWTH_THRESHOLD:
-            cell.physics_representation.add_tail_segment()
-            cell.physics_representation.growth_accumulator_tail = 0.0

SyMBac/physics/segments.py

@@ -106,10 +106,14 @@ def radius(self) -> float:
     @radius.setter
     def radius(self, new_radius: float) -> None:
         if self.use_unsafe_radius_set:
-            self.shape.unsafe_set_radius(new_radius)  # More efficient but is it okay?
+            self.shape.unsafe_set_radius(new_radius)
         else:
             friction, filter = self.shape.friction, self.shape.filter
             self.space.remove(self.shape)
             self.shape = pymunk.Circle(self.body, new_radius)
             self.shape.friction, self.shape.filter = friction, filter
             self.space.add(self.shape)
+
+        self.body.moment = pymunk.moment_for_circle(self.body.mass, 0, new_radius)
+        if self.shape.space is not None:
+            self.shape.space.reindex_shape(self.shape)

SyMBac/physics/simulator.py

@@ -108,6 +108,7 @@ def __init__(
         self.space: pymunk.Space = space
         self.space.threads = physics_config.THREADS
         self.space.iterations = physics_config.ITERATIONS
+        self._baseline_iterations = physics_config.ITERATIONS
         self.space.gravity = physics_config.GRAVITY
         self.space.damping = physics_config.DAMPING
         self.dt: float = physics_config.DT
@@ -322,7 +323,10 @@ def step(self) -> None:
             hook(self)
 
         if self.adaptive_iterations:
-            self.space.iterations = max(10, int(self.space.iterations * 0.9))
+            self.space.iterations = max(
+                self._baseline_iterations,
+                int(self.space.iterations * 0.9),
+            )
             self.max_joint_impulse *= 0.9
 
         self.frame_count += 1

SyMBac/simulation.py

@@ -24,6 +24,20 @@ def _atomic_pickle_dump(payload, output_path):
             os.remove(temp_path)
 
 
+def _handle_lysis(sim, lysis_p):
+    if lysis_p <= 0:
+        return
+
+    candidates = list(sim.colony.cells)
+    np.random.shuffle(candidates)
+    lysis_threshold = norm.ppf(lysis_p)
+    for cell in candidates:
+        if len(sim.colony.cells) <= 1:
+            break
+        if norm.rvs() <= lysis_threshold:
+            sim.colony.delete_cell(cell)
+
+
 class Simulation:
 
     """
@@ -207,6 +221,8 @@ def __init__(
 
         if isinstance(substeps, bool) or not isinstance(substeps, int) or substeps <= 0:
             raise ValueError("substeps must be an integer greater than 0.")
+        if isinstance(lysis_p, bool) or not 0.0 <= lysis_p <= 1.0:
+            raise ValueError("lysis_p must be in [0.0, 1.0].")
         if cell_config_overrides is not None and not isinstance(cell_config_overrides, dict):
             raise TypeError("cell_config_overrides must be a dict when provided.")
         if physics_config_overrides is not None and not isinstance(physics_config_overrides, dict):
@@ -459,13 +475,7 @@ def remove_out_of_bounds(sim):
                     sim.colony.delete_cell(cell)
 
         def handle_lysis(sim):
-            if self.lysis_p <= 0:
-                return
-            for cell in sim.colony.cells[:]:
-                if len(sim.colony.cells) <= 1:
-                    break
-                if norm.rvs() <= norm.ppf(self.lysis_p):
-                    sim.colony.delete_cell(cell)
+            _handle_lysis(sim, self.lysis_p)
 
         def apply_rigid_body_brownian_jitter(sim):
             if (

tests/test_colony_overlap.py

@@ -1,3 +1,7 @@
+from types import SimpleNamespace
+
+import pytest
+
 from SyMBac.physics.colony import Colony
 
 
@@ -12,25 +16,29 @@ def __init__(self, shape):
 
 
 class _PhysicsRepresentation:
-    def __init__(self, shapes, tail_return=None, head_return=None):
+    def __init__(self, shapes, minimum_segments):
         self.segments = [_Segment(shape) for shape in shapes]
-        self._tail_return = tail_return
-        self._head_return = head_return
+        self.minimum_segments = minimum_segments
         self.tail_calls = 0
         self.head_calls = 0
 
     def remove_tail_segment(self):
         self.tail_calls += 1
-        return self._tail_return
+        if len(self.segments) <= self.minimum_segments:
+            return None
+        return self.segments.pop()
 
     def remove_head_segment(self):
         self.head_calls += 1
-        return self._head_return
+        if len(self.segments) <= self.minimum_segments:
+            return None
+        return self.segments.pop(0)
 
 
 class _Cell:
-    def __init__(self, physics_representation):
+    def __init__(self, physics_representation, minimum_segments):
         self.physics_representation = physics_representation
+        self.config = SimpleNamespace(MIN_LENGTH_AFTER_DIVISION=minimum_segments)
 
 
 class _Space:
@@ -41,13 +49,28 @@ def shape_query(self, _shape):
         return [_QueryInfo(self._overlap_shape)]
 
 
-def test_handle_cell_overlaps_stops_cleanly_when_removal_returns_none():
+@pytest.mark.parametrize("mother_segments,daughter_segments", [(2, 3), (3, 2)])
+def test_handle_cell_overlaps_validates_both_cells_before_mutating(
+    mother_segments,
+    daughter_segments,
+):
+    minimum_segments = 2
     overlap_shape = object()
-    mother_pr = _PhysicsRepresentation(shapes=[overlap_shape], tail_return=None)
-    daughter_pr = _PhysicsRepresentation(shapes=[object()], head_return=None)
+    mother_pr = _PhysicsRepresentation(
+        shapes=[overlap_shape, *[object() for _ in range(mother_segments - 1)]],
+        minimum_segments=minimum_segments,
+    )
+    daughter_pr = _PhysicsRepresentation(
+        shapes=[object() for _ in range(daughter_segments)],
+        minimum_segments=minimum_segments,
+    )
+    mother = _Cell(mother_pr, minimum_segments)
+    daughter = _Cell(daughter_pr, minimum_segments)
 
     colony = Colony(space=_Space(overlap_shape), cells=[])
-    colony.handle_cell_overlaps({_Cell(daughter_pr): _Cell(mother_pr)})
+    colony.handle_cell_overlaps({daughter: mother})
 
-    assert mother_pr.tail_calls == 1
-    assert daughter_pr.head_calls == 1
+    assert mother_pr.tail_calls == 0
+    assert daughter_pr.head_calls == 0
+    assert len(mother_pr.segments) == mother_segments
+    assert len(daughter_pr.segments) == daughter_segments