From 63adb45c02f79f045df35539164385c0e8290af0 Mon Sep 17 00:00:00 2001 From: Mike Arpaia Date: Thu, 24 Sep 2026 16:16:39 -0600 Subject: [PATCH] Diagnose out-of-plane motion and document tutorial dimensionality --- docs/architecture/0014-native-controllers.md | 8 +- docs/development/planar-mechanics-followup.md | 40 ++ docs/tutorials/README.md | 2 + docs/tutorials/biophysics-and-growth.md | 2 + docs/tutorials/discrete-state-and-contacts.md | 2 + docs/tutorials/flow-solvers.md | 2 + docs/tutorials/intracellular-dynamics.md | 2 + docs/tutorials/microfluidics.md | 2 + docs/tutorials/planarity.md | 78 ++++ docs/tutorials/signaling.md | 2 + docs/tutorials/simbol.md | 2 + python/src/microsimulator/division.py | 9 +- python/tests/test_planarity.py | 109 ++++++ scripts/diagnose_planarity.py | 350 ++++++++++++++++++ 14 files changed, 607 insertions(+), 3 deletions(-) create mode 100644 docs/development/planar-mechanics-followup.md create mode 100644 docs/tutorials/planarity.md create mode 100644 python/tests/test_planarity.py create mode 100644 scripts/diagnose_planarity.py diff --git a/docs/architecture/0014-native-controllers.md b/docs/architecture/0014-native-controllers.md index 732dd4c..be68f52 100644 --- a/docs/architecture/0014-native-controllers.md +++ b/docs/architecture/0014-native-controllers.md @@ -28,8 +28,12 @@ A controller-backed model resumes through `resume(context, checkpoint)`. The ope 1. compute and validate host regulation; 2. apply cell attributes, species, and fixed-state updates; 3. apply division requests and division callbacks; -4. execute `Simulation.step(dt)` for growth and typed rate plans; and -5. execute exactly `MechanicsConfig.passes` contact/relaxation passes when mechanics is configured. +4. apply requested cell removals; +5. execute `Simulation.step(dt)` for growth and typed rate plans; +6. apply flow drift when enabled in the mechanics configuration; and +7. execute exactly `MechanicsConfig.passes` contact/relaxation passes when mechanics is configured and cells remain. + +`UniformLengthDivision.jitter_z` controls only the random perturbation added to daughter orientation: `None` disables jitter, `False` adds XY-only jitter, and `True` adds XYZ jitter. Native normalization can change an inherited nonzero Z component even when the added Z perturbation is zero. Division places daughters along the parent's three-dimensional axis; subsequent contact relaxation and flow drift remain three-dimensional. See the [planarity diagnostics and tutorial audit](../tutorials/planarity.md) for reproducible examples and the limits of finite-height confinement. The standard payload records a stable model ID and version, completed-step counter, model JSON state, random stream, and every mechanics parameter. `NativeController.from_checkpoint` validates and restores that payload while the checkpoint's native state retains rate plans, signal grids, geometry, and lineage. Exact mechanics passes are a new explicit native-controller contract. The legacy adapter separately preserves the intent of `max_substeps` as a bounded new-contact frontier: it performs at most `max_substeps - 1` solves and stops when rediscovery produces no contact identity not seen earlier in the biological step. This distinction is checkpointed and supported by recorded colony trajectories; native models never inherit the legacy heuristic silently. diff --git a/docs/development/planar-mechanics-followup.md b/docs/development/planar-mechanics-followup.md new file mode 100644 index 0000000..8589de8 --- /dev/null +++ b/docs/development/planar-mechanics-followup.md @@ -0,0 +1,40 @@ +# Proposed follow-up: explicit planar cell mechanics + +Status: specification only; not implemented. The [planarity investigation](../tutorials/planarity.md) found expected three-dimensional responses, including Z-directed degenerate contacts and inherited tilt. `jitter_z=False` must retain its current orientation-perturbation meaning. A strict planar tutorial needs a separately selected native mechanical mode. + +## Public contract and state ownership + +Add a native mechanical dimensionality setting with `spatial_3d` as the backward-compatible default and `planar_xy` with an explicit finite `plane_z`. This is simulation state shared by Python, CPU, Metal, and CUDA, not a viewer option or controller-only callback. Preserve the existing three-coordinate geometry API and existing capsule length, radius, and volume conventions. + +For `planar_xy`, all mobile cell centers satisfy `position.z == plane_z` and all directions satisfy `direction.z == 0`, within a documented float32 representation tolerance after every public geometry mutation and completed native stage. The plane's stored value is its canonical float32 representation. Reject non-finite or unrepresentable heights. Initialization and explicit geometry edits reject appreciably off-plane centers, directions with nonzero Z beyond the declared input tolerance, and directions with no finite nonzero XY component; canonicalize accepted roundoff to exact stored plane Z and normalized XY direction. Do not silently flatten tilted checkpoint geometry or switch modes on an occupied simulation. + +## Translation, rotation, growth, and division + +Solve only the two in-plane translation components and rotation about Z. Assemble/project mechanical Jacobians, forces, increments, and residuals in those degrees of freedom; projecting a completed 3D solve afterwards is insufficient because its contact response was computed in a different space. Fixed cells obey the same planar geometry validation and remain stationary. + +Growth changes length without introducing Z. Division uses the planar parent axis to place daughters, keeps both centers on the configured plane, and preserves the existing volume and lineage contracts. A model requesting an XYZ geometry edit or incompatible division jitter must receive a clear validation error, so a mistakenly configured tutorial does not appear to work through silent filtering. XY-only or absent division jitter works unchanged. + +## Contacts and external boundaries + +Compute planar closest-segment contacts and normals in XY. Crossing rods and coincident parallel rods must receive deterministic in-plane separating directions instead of the 3D cross-product/fallback directions. Specify canonical axis signs, cell-ID ordering, tie breaks, and normal sign conventions centrally; test cell order reversal and both equivalent representations of a rod axis. The same overlapping input and solver parameters must produce matching contact identities and tolerance-equivalent corrections on all available backends. + +The first implementation should support Z-extruded lateral planes, axis-aligned boxes, and Z-aligned cylinders whose planar cross-sections are well-defined and whose finite-height caps fully clear a capsule on the chosen plane. Validate this compatibility when adding constraints and when adding or resizing cells. Reject tilted planes, spheres, or cap intersections until their planar mechanical semantics are explicitly implemented. A floor/ceiling touching a planar capsule must not create an unsatisfiable out-of-plane force. This mode models discs/capsules constrained to a plane within a compatible device; it does not replace finite-height 3D confinement. + +## Flow and chemical fields + +Project sampled drift velocity onto XY before applying mobile-cell translation. Deliberately ignore the normal component as a kinematic constraint and document that this is not a resolved reaction-force or momentum-conservation model. Preserve in-plane interpolation, obstacle handling, boundary behavior, and fixed-cell exclusions. Keep signal grids and chemistry independently three-dimensional: a 3D concentration field may be sampled at `plane_z`; a shallow grid does not automatically opt the cells into planar mode. + +## Checkpoint, resume, and diagnostics + +Version the native checkpoint schema to persist dimensionality and canonical plane height. Old checkpoints migrate explicitly to `spatial_3d`; missing or invalid fields in the new schema fail validation. Validate planar geometry and compatible constraints before exposing a restored simulation, without silently repairing incompatible data. Preserve mode and plane across CPU/Metal/CUDA restoration, controller restart, and clone/export paths. Reject a resume request whose requested mechanics mode conflicts with saved state. + +Expose the mode and plane in scene/analysis metadata so diagnostics can distinguish an intended invariant from a visual appearance. Coordinate that schema change with the existing metadata owner; do not infer dimensionality from channel labels, camera view, signal-grid depth, or jitter configuration. + +## Acceptance and validation + +- Shared native fixtures cover separated cells, overlapping parallel rods, crossing rods, order-reversed pairs, arbitrary in-plane orientations, division, long growth runs, and mixed mobile/fixed cells. Check both center Z and direction Z after each relevant stage, including zero-duration controller steps. +- Test all public construction and geometry-edit paths: accepted roundoff canonicalizes consistently; inherited tilt and invalid heights fail clearly; default 3D behavior and its existing conformance fixtures remain unchanged. +- Verify planar contacts separate overlap in XY with finite residuals and deterministic identities, including degenerate ties. Check force/rotation consistency and convergence rather than only final projection onto the plane. +- Exercise each supported boundary, each rejected incompatible boundary, cell growth approaching a cap, and a prescribed flow with nonzero Z velocity. In-plane drift is preserved and out-of-plane drift is suppressed only in planar mode. +- Round-trip checkpoints with each mode, migrate old data, reject malformed/inconsistent planar state, and compare resumed trajectories against uninterrupted execution. Run the same fixture contract on CPU, Metal, and CUDA; report unavailable accelerator hardware instead of substituting CPU. +- Add an opt-in tutorial using the new mode and a paired finite-height 3D example. Documentation explains the distinct mechanical assumptions and retains the current definition of `jitter_z`. diff --git a/docs/tutorials/README.md b/docs/tutorials/README.md index e55816a..1399f97 100644 --- a/docs/tutorials/README.md +++ b/docs/tutorials/README.md @@ -28,6 +28,8 @@ These lessons develop the biological rules used within devices and in standalone The [analysis tutorial](analysis.md) covers checkpoints, contact graphs, and quantitative output. Continue with [analysis recipes](../analysis/recipes.md) for reproducible Parquet/Zarr datasets and Polars queries. +The [dimensionality audit and planarity diagnostic](planarity.md) explain XY-only division jitter, finite-height confinement, and reproducible causes of out-of-plane cell motion. + ## Working with the examples Teaching models are under [`examples/tutorials`](../../examples/tutorials). Scenario parameters are JSON values passed with `--parameter`; every command in the tutorials can be run from the repository root. diff --git a/docs/tutorials/biophysics-and-growth.md b/docs/tutorials/biophysics-and-growth.md index 0b72a64..acac890 100644 --- a/docs/tutorials/biophysics-and-growth.md +++ b/docs/tutorials/biophysics-and-growth.md @@ -2,6 +2,8 @@ This tutorial introduces cell geometry, growth, division, lineage, cell types, mechanical constraints, and competition. Its five runnable scenarios are defined in `examples/tutorials/biophysics.py`. +The `basics`, `two_types`, and `competition` scenarios add XY-only division jitter but retain unrestricted 3D mechanics. `short_cells` adds XYZ jitter. `box` also adds XYZ jitter and has a floor and four lateral walls, with no ceiling. None guarantees a planar colony; see [division jitter and out-of-plane motion](planarity.md) for the full contract and reproducible diagnostics. + ## 1. A founder that grows and divides Run the basic model: diff --git a/docs/tutorials/discrete-state-and-contacts.md b/docs/tutorials/discrete-state-and-contacts.md index c7b8e5e..722b2f0 100644 --- a/docs/tutorials/discrete-state-and-contacts.md +++ b/docs/tutorials/discrete-state-and-contacts.md @@ -2,6 +2,8 @@ This tutorial uses plasmid segregation and conjugation to show how discrete biological state, stochastic events, and contact-dependent behavior fit into a MicroSimulator model. +Both models start with centers at Z=0 and axes in XY, and division adds no orientation jitter. Neither has mechanical Z confinement: daughters inherit the parent axis and contact relaxation remains three-dimensional. See [division jitter and out-of-plane motion](planarity.md). + New founders preserve their requested length unless it exceeds the single sampled division target. This also handles the rare short Gaussian target in the conjugation model without rejection sampling. Checkpoint restoration keeps stored lengths and targets. See [founder initialization](biophysics-and-growth.md#length-and-volume). ## 1. Incompatible plasmid segregation diff --git a/docs/tutorials/flow-solvers.md b/docs/tutorials/flow-solvers.md index 27ec792..0f04742 100644 --- a/docs/tutorials/flow-solvers.md +++ b/docs/tutorials/flow-solvers.md @@ -2,6 +2,8 @@ The [pillar-channel model](../../examples/tutorials/pillar_channel.py) combines cylindrical walls, a depth-integrated flow calculation, attached founder lineages, and released daughters: +Cells retain three-dimensional mechanics inside the channel walls at Z=±3. XY-only division jitter and depth-integrated flow do not impose a planar cell constraint; fixed founders remain attached while released daughters can move and tilt within the finite-height chamber. See the [dimensionality audit](planarity.md). + ```console uv run microsimulator view --model examples/tutorials/pillar_channel.py --seed 7 --dt 0.01 --backend metal --open ``` diff --git a/docs/tutorials/intracellular-dynamics.md b/docs/tutorials/intracellular-dynamics.md index 1bcb155..619f2d8 100644 --- a/docs/tutorials/intracellular-dynamics.md +++ b/docs/tutorials/intracellular-dynamics.md @@ -2,6 +2,8 @@ This tutorial introduces intracellular concentrations, growth dilution, typed rate equations, gene-circuit feedback, and quantitative time-course analysis. The runnable scenarios are collected in `examples/tutorials/gene_expression.py`. +All five scenarios use XY-only division jitter and start with centers at Z=0, but add no mechanical walls. Their cells retain three-dimensional translations and rotations. See [division jitter and out-of-plane motion](planarity.md) before treating a planar-looking trajectory as a strict 2D model. + All five gene-expression scenarios request a founder centerline length of 3.5 and cap it at the one sampled division target. Initial concentrations are unchanged; the smaller biomass can change total initial amount. See [founder initialization and volume conventions](biophysics-and-growth.md#length-and-volume). ## The native species contract diff --git a/docs/tutorials/microfluidics.md b/docs/tutorials/microfluidics.md index 36d68b9..bea083d 100644 --- a/docs/tutorials/microfluidics.md +++ b/docs/tutorials/microfluidics.md @@ -2,6 +2,8 @@ This tutorial connects device geometry, flowing media, and cell biology in runnable MicroSimulator models. The [modeling guide](../microfluidics.md) introduces the workflow and the choice of flow solver. Four examples cover the range: +These models use XY-only division jitter and finite-height 3D confinement. A thin cavity can encourage a monolayer, but it does not force a common center Z or eliminate tilt; walls are soft constraints whose residual depends on relaxation tolerance and passes. The [dimensionality audit](planarity.md) lists each device and reproduces these distinctions. + The microfluidic-trap, Danino, biopixel, and pillar tutorial founders request centerline length 3.5, capped at their single sampled target in [3.2, 3.8]. Attachment, position, radius, and concentrations are preserved. This affects new construction only; saved geometry is restored unchanged. See [founder initialization and volume conventions](biophysics-and-growth.md#length-and-volume). | Model | Device | Demonstrates | diff --git a/docs/tutorials/planarity.md b/docs/tutorials/planarity.md new file mode 100644 index 0000000..e1366ca --- /dev/null +++ b/docs/tutorials/planarity.md @@ -0,0 +1,78 @@ +# Division jitter, confinement, and out-of-plane motion + +MicroSimulator's native cell mechanics is three-dimensional. An initially planar colony and XY-only division jitter can remain planar in a particular run, but neither imposes a planar constraint. A shallow flow calculation or a one-layer signal grid does not change the mechanical degrees of freedom. + +`UniformLengthDivision(jitter_z=False)` draws the usual three random perturbations and replaces the new Z perturbation with zero. It adds the resulting vector to each daughter's inherited direction, then the engine normalizes that direction. An inherited nonzero Z component therefore survives, and its normalized value can change when X or Y changes. Daughter centers are placed along the parent's three-dimensional axis before jitter is applied. `jitter_z=None` disables orientation jitter entirely; `jitter_z=True` includes the new Z perturbation. All three settings leave subsequent contact relaxation and flow drift three-dimensional. + +## Tutorial audit + +The following contracts describe the checked-in models after the founder-initialization correction in issue #22. “XY-only” refers to the added division perturbation, not a mechanical constraint. No model below implements strict 2D mechanics. + +| Model or scenario | Initial plane and division jitter | Mechanical confinement and dimensional contract | +| --- | --- | --- | +| `biophysics.py`: `basics`, `two_types`, `competition` | Centers Z=0, directions in XY; XY-only jitter | No walls. All translations and rotations remain 3D. | +| `biophysics.py`: `short_cells` | Centers Z=0; XYZ jitter | No walls. Orientation leaves XY at division by design. | +| `biophysics.py`: `box` | Centers Z=0.5; XYZ jitter | Floor at Z=0 and four lateral walls; no ceiling. This is an open 3D box. | +| `gene_expression.py`: all five scenarios | Centers Z=0, directions in XY; XY-only jitter | No mechanical walls; unrestricted 3D. | +| `signaling.py`: `single_gene`, `communication` | Centers Z=0, directions in XY; XY-only jitter | Only lateral Y walls at ±16; no Z confinement. | +| `signaling.py`: `mutualism` | Centers Z=0, directions in XY; XY-only jitter | No mechanical walls; unrestricted 3D. | +| `simbol_circuits.py`: all six circuits | Centers Z=0, directions in XY; XY-only jitter | No mechanical walls; unrestricted 3D. | +| `plasmid_segregation.py`, `conjugation.py` | Centers Z=0, directions in XY; no orientation jitter | Division inherits the parent axis; contact mechanics remains unrestricted 3D. | +| `pillar_channel.py` | Centers Z=0, directions in XY; XY-only jitter | Box walls at Z=±3 and cylindrical pillars. Attached cells are fixed; released daughters move in a finite-height 3D chamber. | +| `danino_clock.py`, `examples/microfluidic_trap.py` | Centers Z=0, directions in XY; XY-only jitter | Trap walls at Z=±3. These permit nonzero Z and tilt within a finite-height 3D device. | +| `biopixel_trap.py` | Centers Z=0.825, directions in XY; XY-only jitter | Trap floor Z=0 and roof Z=1.65; the channel is taller. The thin trap encourages a monolayer but does not force identical center Z or zero tilt. | +| `examples/culture_dish.py` | Centers Z=0, directions in XY; XY-only jitter | Inside cylinder with caps at Z=±1. Finite-height confinement permits nonzero Z and tilt. | + +“Monolayer” in the device tutorials describes a physical modeling intention supported by a thin cavity. Walls are soft numerical constraints, and residual overlap depends on the contact solver tolerance and number of relaxation passes. Their presence does not guarantee exact wall bounds after one pass or eliminate out-of-plane degrees of freedom. + +## Reproduce and locate a departure + +Run the diagnostic from the repository root after installing the development environment: + +```console +uv run python scripts/diagnose_planarity.py --backend cpu --seed 17 --dt 0.02 --output results/planarity-fixtures.json +uv run python scripts/diagnose_planarity.py --backend cpu --seed 17 --dt 0.02 --model examples/tutorials/biophysics.py --scenario basics --steps 1000 --max-cells 128 --output results/planarity-basics.json +uv run python scripts/diagnose_planarity.py --backend metal --seed 17 --dt 0.02 --model examples/tutorials/biopixel_trap.py --plane-z 0.825 --steps 1000 --max-cells 128 --output results/planarity-biopixel.json +``` + +Use `--plane-z 0.5` for the box scenario. Other model parameters accept JSON, for example `--parameter circuit='"bba_0001"'`. The script rejects unavailable backends instead of falling back. A division batch can exceed `--max-cells`; that limit stops the next biological step rather than truncating the model's division requests. + +The JSON records the backend, seed, timestep, model path and supplied parameters, model/controller state, source provenance, mechanics settings, initial and final geometry, and constraints. Every stage records maximum absolute center Z, displacement from the selected reference plane, absolute direction Z, and the change in each Z component for surviving cell IDs. The first event exceeding `1e-6` includes before/after geometry. This is a diagnostic detection tolerance, not a global engine guarantee. For an initially tilted cell, the first event is initialization; later stages still report their changes. + +The probe wraps native Python entry points for the selected simulation in a dedicated process and restores them in `finally`. It observes division, geometry edits, removal, growth/chemistry, flow drift, and contact/constraint relaxation without duplicating the controller or drawing random numbers. Contact and wall relaxation use the same native entry point: the isolated contact-only and wall-only fixtures distinguish their causes. Arbitrary native work performed internally by a custom extension is outside these Python stage boundaries. Use this as a headless diagnostic, not inside a multithreaded application. + +## Findings and classification + +The seven fixtures in `scripts/diagnose_planarity.py` and their shared-backend assertions in `python/tests/test_planarity.py` establish the following causes independently of any long colony trajectory. The measurements below were reproduced on CPU and Apple Metal with seed 17; NVIDIA CUDA hardware was unavailable for this investigation. The same tests select CUDA when its runtime is available. + +| Fixture | Observation | Classification | +| --- | --- | --- | +| Separated planar rods | Two nonoverlapping X-oriented rods retain center Z=0 and direction Z=0 through relaxation. | Planar control case; no departure. | +| Crossing rods | Two rods at the same center, directions X and Y, centerline length 4 and radius 0.5, produce a Z-directed normal. One relaxation moves their centers to approximately ±0.4. | Expected 3D contact response. Crossing rods at zero centerline separation are an overlapping initialization if intended as a nonoverlapping planar colony. | +| Coincident parallel rods | Two identical X-oriented rods at one center use the deterministic degenerate-contact fallback; the normal points in Z and relaxation separates centers to approximately ±0.4. | Expected 3D degeneracy handling, also an overlapping initialization. | +| Planar division | An intentionally oversized planar parent divides with XY-only jitter without creating a Z component. | Planar division control; no departure. The fixture intentionally bypasses tutorial founder capping. | +| Inherited tilt | A parent initialized with direction `(1, 0, 0.2)` has normalized Z≈0.196116. Division places daughters at Z≈±0.245145. Both geometry-edit calls request exactly zero added Z, although native normalization changes the daughters' direction Z. | Expected inherited geometry and normalization; an initially tilted setup cannot test preservation of a planar state. | +| Finite-height walls | A horizontal radius-0.5 rod at Z=0.8 intersects the ceiling at Z=1. One relaxation moves its center to approximately 0.6. Repeated default solves stop with approximately 0.003704 penetration; tightening the residual tolerance to `1e-8` reduces penetration below `1e-6`. The center remains near Z=0.5. | Expected soft-wall convergence within a finite-height 3D space. Default residual tolerance is `0.005`; physical confinement neither implies exact wall projection nor the plane Z=0. | +| Vertical flow | Prescribed Z velocity 0.2 over dt=0.02 moves a center from Z=1 to approximately 1.004, with no division. | Expected 3D advection. XY-only jitter does not filter flow. | + +For zero closest-point separation, the CPU, Metal, and CUDA contact implementations first try the cross product of rod axes, then transverse center separation, then a deterministic perpendicular fallback. Two crossing XY axes have a Z-directed cross product; coincident X axes use a fallback that can also point in Z. Choosing these directions is consistent with the existing 3D contact contract. Suppressing them globally would change three-dimensional mechanics. + +Representative tutorial runs used seed 17, dt=0.02, a maximum of 1000 steps, and a stop threshold of 128 cells. CPU and Metal reached the same step/count endpoints below; agreement here is not a claim of bitwise trajectory equivalence. Geometry was checked after each instrumented stage, including both center displacement and direction Z. + +| Model/scenario | Reference Z | Completed steps / final cells | First departure | +| --- | ---: | ---: | --- | +| Biophysics: basics | 0 | 312 / 128 | None above diagnostic tolerance | +| Biophysics: two_types | 0 | 149 / 128 | None above diagnostic tolerance | +| Biophysics: competition | 0 | 351 / 130 | None above diagnostic tolerance | +| Biophysics: short_cells | 0 | 325 / 128 | Division geometry edit at time≈0.02: direction Z≈0.000531725 | +| Biophysics: box | 0.5 | 327 / 128 | Division geometry edit at time≈0.02: direction Z≈0.000531725 | +| Gene expression: constitutive | 0 | 327 / 128 | None above diagnostic tolerance | +| Signaling: single_gene | 0 | 179 / 128 | None above diagnostic tolerance | +| Pillar channel | 0 | 394 / 131 | None above diagnostic tolerance | +| Biopixel trap | 0.825 | 478 / 128 | None above diagnostic tolerance | + +The two departures in this matrix are the tutorials' intentional XYZ jitter. The source audit covers additional scenarios without claiming that they all received these trajectory runs. The unchanged plane in the other runs provides reproducibility evidence for those particular initial states and durations, not proof of planar invariance. + +No solver defect was demonstrated. Issue #22 separately corrected oversized tutorial founders, which could previously divide immediately; the diagnostics here use that corrected initialization. Luiza's older-version observation cannot be assigned to one particular mechanism without its original model, geometry, seed, and trajectory. The fixtures nevertheless show several reproducible paths to Z motion even with `jitter_z=False`. + +If a lesson requires every center to remain at a specified Z and every axis to remain in XY, it requires an explicit planar mechanics feature. The [proposed planar-mechanics contract](../development/planar-mechanics-followup.md) defines that follow-up without changing the meaning of division jitter. diff --git a/docs/tutorials/signaling.md b/docs/tutorials/signaling.md index e617be3..f8e42ef 100644 --- a/docs/tutorials/signaling.md +++ b/docs/tutorials/signaling.md @@ -2,6 +2,8 @@ This tutorial introduces extracellular grids, diffusion, cell-grid exchange, sender-receiver communication, and two-strain mutualism. Run the scenarios in `examples/tutorials/signaling.py` with a small time step such as `0.01`. +All scenarios use XY-only division jitter with three-dimensional mechanics. `single_gene` and `communication` have lateral Y walls only; `mutualism` has no mechanical walls. Signal-grid depth does not constrain cell Z or tilt. See [division jitter and out-of-plane motion](planarity.md). + Each new founder requests centerline length 3.5 and is capped at its sampled division threshold. Its radius, position, cell type, and initial concentrations are preserved. See [founder initialization](biophysics-and-growth.md#length-and-volume). ## Grid geometry and units diff --git a/docs/tutorials/simbol.md b/docs/tutorials/simbol.md index 9b3908d..6bffc75 100644 --- a/docs/tutorials/simbol.md +++ b/docs/tutorials/simbol.md @@ -2,6 +2,8 @@ SimBOL connects an SBOL 3 design to simulator-specific code through a summarized JSON representation. This tutorial presents typed MicroSimulator versions of six BioBrick circuit examples and a spatial quorum-sensing clock. +The six circuit models start in XY and use XY-only division jitter without mechanical walls. The Danino clock uses a finite-height trap. Both retain three-dimensional mechanics; see [division jitter, confinement, and out-of-plane motion](planarity.md). + These are explicit example models, not a general SBOL-to-rate-plan import path. The [source reference](../compatibility/tutorial-source-provenance.md#simbol-source-workflow) describes how they relate to the SimBOL notebook, generated Python, and JSON fixtures. ## Run the six circuits diff --git a/python/src/microsimulator/division.py b/python/src/microsimulator/division.py index a8c756e..6105614 100644 --- a/python/src/microsimulator/division.py +++ b/python/src/microsimulator/division.py @@ -40,7 +40,14 @@ def capped_founder_length(requested: float, target: float) -> float: @dataclass(frozen=True, slots=True) class UniformLengthDivision: - """Divide above per-cell thresholds sampled from one uniform distribution.""" + """Divide above per-cell thresholds sampled from one uniform distribution. + + ``jitter_z=None`` disables orientation jitter; ``False`` adds XY-only + perturbations; ``True`` adds XYZ perturbations. XY-only jitter leaves the + inherited Z component unchanged before native direction normalization, + which can change its normalized value. Neither option constrains daughter + positions, contact mechanics, or flow drift to a plane. + """ minimum: float maximum: float diff --git a/python/tests/test_planarity.py b/python/tests/test_planarity.py new file mode 100644 index 0000000..1a99098 --- /dev/null +++ b/python/tests/test_planarity.py @@ -0,0 +1,109 @@ +"""Shared native-backend regressions for the documented three-dimensional contract.""" + +from __future__ import annotations + +import math +import runpy +from itertools import pairwise +from pathlib import Path + +import pytest +from microsimulator import ( + BackendKind, + ModelContext, + NativeController, + Simulation, + backend_available, + build_model, +) +from microsimulator.scene import capture_scene, dumps_scene + +ROOT = Path(__file__).resolve().parents[2] +DIAGNOSTIC = runpy.run_path(str(ROOT / "scripts" / "diagnose_planarity.py")) + + +@pytest.mark.parametrize("backend", list(BackendKind)) +def test_three_dimensional_diagnostic_fixtures(backend: BackendKind) -> None: + if not backend_available(backend): + pytest.skip(f"{backend} runtime unavailable; no fallback") + results = DIAGNOSTIC["fixtures"](backend, 17, 0.02) + for name in ("separated_planar", "planar_division"): + assert results[name]["first_out_of_plane"] is None + assert all( + stage["max_center_displacement_from_plane"] < 1e-6 for stage in results[name]["stages"] + ) + assert all(stage["max_direction_z"] < 1e-6 for stage in results[name]["stages"]) + for name in ("crossing", "coincident_parallel"): + result = results[name] + assert any(abs(normal[2]) > 0.99 for normal in result["contact_normals"]) + assert result["first_out_of_plane"]["stage"] == "contact_and_constraint_relaxation" + negative, positive = sorted(cell["center"][2] for cell in result["final_geometry"]) + assert math.isclose(negative, -0.4, abs_tol=1e-6) + assert math.isclose(positive, 0.4, abs_tol=1e-6) + inherited = results["inherited_tilt"] + assert inherited["first_out_of_plane"]["stage"] == "initialization" + assert any( + stage["max_center_displacement_from_plane"] > 0.2 + for stage in inherited["stages"] + if stage["stage"] == "division" + ) + jitter = [stage for stage in inherited["stages"] if stage["stage"] == "geometry_edit"] + assert len(jitter) == 2 + assert all(stage["requested_direction_delta_z"] == 0 for stage in jitter) + assert any(stage["max_stage_direction_change_z"] > 1e-6 for stage in jitter) + confined = results["finite_height_constraints"] + default = confined["default_wall_penetration_by_pass"] + tight = confined["tight_wall_penetration_by_pass"] + assert math.isclose(default[0], 0.3, abs_tol=1e-6) + assert 1e-6 < default[-1] < 0.005 + assert tight[-1] < 1e-6 + assert all(after <= before + 1e-7 for before, after in pairwise(tight)) + final = confined["final_geometry"][0] + axial_z_extent = abs(final["direction"][2]) * final["length"] / 2 + final["radius"] + assert final["center"][2] + axial_z_extent <= 1 + 1e-6 + assert final["center"][2] > 0.49 # walls permit nonzero Z; they do not impose z=0 + flow = results["vertical_flow"] + assert flow["first_out_of_plane"]["stage"] == "flow_drift" + assert math.isclose(flow["final_geometry"][0]["center"][2], 1.004, abs_tol=1e-6) + + +def test_diagnostics_preserve_rng_and_restore_native_methods_on_failure() -> None: + model_path = ROOT / "examples/tutorials/biophysics.py" + traced, _ = build_model( + model_path, ModelContext(BackendKind.CPU, 0, seed=17, parameters={"scenario": "two_types"}) + ) + normal, _ = build_model( + model_path, ModelContext(BackendKind.CPU, 0, seed=17, parameters={"scenario": "two_types"}) + ) + assert isinstance(traced, NativeController) + assert isinstance(normal, NativeController) + trace = DIAGNOSTIC["Trace"](traced.simulation) + originals = { + name: getattr(Simulation, name) + for name in ( + "divide", + "divide_equal", + "remove_cell", + "step", + "apply_flow_drift", + "relax_cell_mechanics", + "set_cell_geometry", + ) + } + with trace.instrument(): + for _ in range(5): + traced.step(0.02) + for _ in range(5): + normal.step(0.02) + assert dumps_scene(capture_scene(traced.simulation)) == dumps_scene( + capture_scene(normal.simulation) + ) + assert traced.controller_state() == normal.controller_state() + assert traced.simulation.cell_count > len(trace.initial["cells"]) + assert any(stage["stage"] == "geometry_edit" for stage in trace.events) + assert all(getattr(Simulation, name) is method for name, method in originals.items()) + with pytest.raises(RuntimeError, match="probe"), trace.instrument(): + raise RuntimeError("probe") + assert all(getattr(Simulation, name) is method for name, method in originals.items()) + assert trace.initial["cells"] + assert "constraints" in trace.report() diff --git a/scripts/diagnose_planarity.py b/scripts/diagnose_planarity.py new file mode 100644 index 0000000..7aad6b9 --- /dev/null +++ b/scripts/diagnose_planarity.py @@ -0,0 +1,350 @@ +"""Record stage-local Z motion without changing the 3D solver or its random stream. + +Run in a dedicated diagnostic process: native Python entry points are temporarily +instrumented and restored on exit. No controller implementation is duplicated. +""" + +from __future__ import annotations + +import argparse +import json +import math +import random +from contextlib import contextmanager +from pathlib import Path + +from microsimulator import ( + BackendKind, + CellInit, + GridBoundary, + GridBoundaryKind, + GridShape, + MechanicsConfig, + ModelContext, + NativeController, + PlaneConstraintInit, + SignalGridSpec, + SignalGridVelocityField, + Simulation, + StepPlan, + UniformLengthDivision, + Vec3, + backend_available, + build_model, +) +from microsimulator.scene import capture_scene, dumps_scene + + +def cells(simulation): + return [ + dict( + id=c.id, + center=[c.position.x, c.position.y, c.position.z], + direction=[c.direction.x, c.direction.y, c.direction.z], + length=c.length, + radius=c.radius, + ) + for c in simulation.cells() + ] + + +class Trace: + def __init__(self, simulation, plane_z=0.0, tolerance=1e-6): + self.simulation = simulation + self.plane_z = plane_z + self.tolerance = tolerance + self.events = [] + self.first_event = None + self.initial = json.loads(dumps_scene(capture_scene(simulation)))["frame"] + self.record("initialization", [], cells(simulation)) + + def record(self, stage, before, after, requested_z=None): + entry = dict( + stage=stage, + time=self.simulation.time, + cells=len(after), + max_center_z=max((abs(c["center"][2]) for c in after), default=0.0), + max_center_displacement_from_plane=max( + (abs(c["center"][2] - self.plane_z) for c in after), default=0.0 + ), + max_direction_z=max((abs(c["direction"][2]) for c in after), default=0.0), + ) + previous = {c["id"]: c for c in before} + entry["max_stage_center_change_z"] = max( + ( + abs(c["center"][2] - previous[c["id"]]["center"][2]) + for c in after + if c["id"] in previous + ), + default=0.0, + ) + entry["max_stage_direction_change_z"] = max( + ( + abs(c["direction"][2] - previous[c["id"]]["direction"][2]) + for c in after + if c["id"] in previous + ), + default=0.0, + ) + if requested_z is not None: + entry["requested_direction_delta_z"] = requested_z + self.events.append(entry) + if ( + self.first_event is None + and max(entry["max_center_displacement_from_plane"], entry["max_direction_z"]) + > self.tolerance + ): + self.first_event = dict(**entry, before=before, after=after) + + @contextmanager + def instrument(self): + names = { + "divide": "division", + "divide_equal": "division", + "remove_cell": "removal", + "step": "growth_and_chemistry", + "apply_flow_drift": "flow_drift", + "relax_cell_mechanics": "contact_and_constraint_relaxation", + "set_cell_geometry": "geometry_edit", + } + originals = {} + + def wrapped(name, original): + def call(simulation, *args, **kwargs): + if simulation is not self.simulation: + return original(simulation, *args, **kwargs) + before = cells(simulation) + requested_z = None + if name == "set_cell_geometry": + cell_id = args[0] if args else kwargs.get("cell_id") + direction = args[2] if len(args) >= 3 else kwargs.get("direction") + if cell_id is not None and direction is not None: + previous = simulation.cell(cell_id) + requested_z = direction.z - previous.direction.z + result = original(simulation, *args, **kwargs) + self.record(names[name], before, cells(simulation), requested_z) + return result + + return call + + try: + for name in names: + originals[name] = getattr(Simulation, name) + setattr(Simulation, name, wrapped(name, originals[name])) + yield + finally: + for name, original in originals.items(): + setattr(Simulation, name, original) + + def report(self): + return dict( + plane_z=self.plane_z, + tolerance=self.tolerance, + initial_geometry=self.initial["cells"], + constraints=self.initial["constraints"], + backend=self.initial["backend"], + first_out_of_plane=self.first_event, + final_geometry=cells(self.simulation), + stages=self.events, + ) + + +def add(simulation, position=(0, 0, 0), direction=(1, 0, 0), length=4.0, radius=0.5): + cell = CellInit() + cell.position, cell.direction = Vec3(*position), Vec3(*direction) + cell.length, cell.radius, cell.growth_rate = length, radius, 0.0 + return simulation.add_cell(cell) + + +def fixtures(backend, seed, dt): + results = {} + for name in ("separated_planar", "crossing", "coincident_parallel"): + simulation = Simulation(backend) + add(simulation) + add( + simulation, + position=(0, 3, 0) if name == "separated_planar" else (0, 0, 0), + direction=(0, 1, 0) if name == "crossing" else (1, 0, 0), + ) + normals = [ + [c.normal.x, c.normal.y, c.normal.z] for c in simulation.find_cell_contacts().contacts + ] + trace = Trace(simulation) + with trace.instrument(): + simulation.relax_cell_mechanics(*MechanicsConfig().native_parameters()) + results[name] = dict( + classification="expected 3D contact behavior", + mechanics=MechanicsConfig().to_json(), + contact_normals=normals, + **trace.report(), + ) + for name, direction in (("planar_division", (1, 0, 0)), ("inherited_tilt", (1, 0, 0.2))): + simulation = Simulation(backend) + founder = add(simulation, direction=direction) + rng = random.Random(seed) + policy = UniformLengthDivision(3.0, 3.0, jitter_z=False) + state = {} + policy.initialize(state, rng, (founder,)) # intentionally oversized diagnostic founder + controller = NativeController( + simulation, + model_id="planarity-diagnostic", + model_version=1, + rng=rng, + state=state, + regulate=lambda step, policy=policy: StepPlan(divisions=policy.requests(step)), + on_division=policy.on_division, + ) + trace = Trace(simulation) + with trace.instrument(): + controller.step(0.0) + results[name] = dict( + classification="expected inherited geometry and normalized XY jitter", + controller_step_dt=0.0, + division_target=3.0, + jitter_z=False, + mechanics=MechanicsConfig().to_json(), + **trace.report(), + ) + simulation = Simulation(backend) + add(simulation, position=(0, 0, 0.8)) + for height, normal in ((-1.0, 1.0), (1.0, -1.0)): + plane = PlaneConstraintInit() + plane.point, plane.inward_normal = Vec3(0, 0, height), Vec3(0, 0, normal) + simulation.add_plane_constraint(plane) + trace = Trace(simulation, plane_z=0.8) + + def relax_to_tolerance(config): + residuals = [] + for _ in range(20): + residuals.append( + max( + ( + max(0.0, -contact.signed_separation) + for contact in simulation.find_external_contacts().contacts + ), + default=0.0, + ) + ) + if residuals[-1] <= 1e-6 or (len(residuals) > 1 and residuals[-1] == residuals[-2]): + break + simulation.relax_cell_mechanics(*config.native_parameters()) + return residuals + + default_config = MechanicsConfig() + tight_config = MechanicsConfig(residual_rms_tolerance=1e-8) + with trace.instrument(): + default_residuals = relax_to_tolerance(default_config) + default_geometry = cells(simulation) + tight_residuals = relax_to_tolerance(tight_config) + results["finite_height_constraints"] = dict( + classification="expected 3D wall relaxation; finite height is not strict 2D", + default_mechanics=default_config.to_json(), + tight_mechanics=tight_config.to_json(), + default_wall_penetration_by_pass=default_residuals, + tight_wall_penetration_by_pass=tight_residuals, + default_final_geometry=default_geometry, + **trace.report(), + ) + simulation = Simulation(backend) + shape = GridShape() + shape.x, shape.y, shape.z = 1, 1, 3 + spec = SignalGridSpec() + spec.shape, spec.signal_count, spec.diffusion, spec.advection = shape, 1, [0.0], [Vec3()] + boundary = GridBoundary() + boundary.kind, boundary.values = GridBoundaryKind.FIXED, [0.0] + spec.z_lower, spec.z_upper = boundary, boundary + field = SignalGridVelocityField() + field.x_faces, field.y_faces, field.z_faces = [0.0] * 6, [0.0] * 6, [0.2] * 4 + spec.velocity_field = field + simulation.configure_signal_grid(spec) + add(simulation, position=(0, 0, 1), length=1.0) + trace = Trace(simulation, plane_z=1.0) + with trace.instrument(): + simulation.apply_flow_drift(dt) + results["vertical_flow"] = dict( + classification="expected prescribed 3D advection", + prescribed_velocity=[0.0, 0.0, 0.2], + drift_dt=dt, + **trace.report(), + ) + return results + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--backend", choices=("cpu", "metal", "cuda"), default="cpu") + parser.add_argument("--seed", type=int, default=17) + parser.add_argument("--dt", type=float, default=0.02) + parser.add_argument("--model", type=Path) + parser.add_argument("--scenario", help="Convenience alias for a JSON string scenario parameter") + parser.add_argument("--parameter", action="append", default=[], metavar="NAME=JSON") + parser.add_argument("--steps", type=int, default=1000) + parser.add_argument("--max-cells", type=int, default=128) + parser.add_argument("--plane-z", type=float, default=0.0) + parser.add_argument("--output", type=Path) + args = parser.parse_args() + if not math.isfinite(args.dt) or args.dt <= 0 or args.steps < 1 or args.max_cells < 1: + parser.error("dt, steps and max-cells must be positive and finite") + if not math.isfinite(args.plane_z): + parser.error("plane-z must be finite") + if not args.model and (args.scenario is not None or args.parameter): + parser.error("scenario and parameters require a model") + backend = getattr(BackendKind, args.backend.upper()) + if not backend_available(backend): + parser.error(f"{args.backend} backend unavailable; no fallback performed") + result = dict( + diagnostic_version=1, + seed=args.seed, + dt=args.dt, + backend=args.backend, + fixtures=fixtures(backend, args.seed, args.dt), + ) + if args.model: + parameters = {} + for parameter in args.parameter: + try: + key, value = parameter.split("=", 1) + if not key or key in parameters: + raise ValueError("parameter names must be nonempty and unique") + parameters[key] = json.loads(value) + except (ValueError, json.JSONDecodeError) as error: + parser.error(f"invalid parameter {parameter!r}: {error}") + if args.scenario is not None: + if "scenario" in parameters: + parser.error("provide scenario once, using --scenario or --parameter") + parameters["scenario"] = args.scenario + model, provenance = build_model( + args.model, + ModelContext(backend, 0, seed=args.seed, parameters=parameters), + ) + trace = Trace(model.simulation, args.plane_z) + completed = 0 + with trace.instrument(): + for _ in range(args.steps): + if model.simulation.cell_count >= args.max_cells: + break + model.step(args.dt) + completed += 1 + result["tutorial"] = dict( + model=str(args.model), + parameters=parameters, + model_state=model.controller_state().get("model"), + controller_kind=model.controller_state().get("kind"), + mechanics=model.controller_state().get("mechanics"), + completed_steps=completed, + stop_reason="max_cells" if model.simulation.cell_count >= args.max_cells else "steps", + provenance=provenance, + requested_steps=args.steps, + max_cells=args.max_cells, + **trace.report(), + ) + encoded = json.dumps(result, indent=2, allow_nan=False) + "\n" + if args.output: + args.output.parent.mkdir(parents=True, exist_ok=True) + args.output.write_text(encoded) + else: + print(encoded, end="") + + +if __name__ == "__main__": + main()