Simulations
InterSIM bridge utilities and semi-synthetic trajectory generators for multi-omics simulation studies.
R dependency
The bridge is optional and requires Rscript plus the R InterSIM package:
install.packages(
"InterSIM",
repos = c("https://cran.r-universe.dev", "https://cloud.r-project.org")
)
Check availability before running a simulation:
from motco.simulations import check_intersim_available
availability = check_intersim_available()
if not availability.available:
print(availability.message)
Example
from motco.simulations import InterSIMParams, run_intersim
result = run_intersim(
InterSIMParams(
seed=1203,
n_sample=100,
cluster_sample_prop=(0.3, 0.3, 0.4),
delta_methyl=1.0,
delta_expr=1.0,
delta_protein=1.0,
p_dmp=0.1,
)
)
methylation = result.methylation
expression = result.expression
proteomics = result.proteomics
clusters = result.clusters
Semi-synthetic trajectory generation
The semi-synthetic trajectory generator samples directly from the numpy-native
InterSIM reimplementation (no R at runtime) and returns a MOTCO-ready dataset with:
- aligned methylation, gene expression, and proteomics matrices
- sample metadata containing
sample_id, group, and stage
- truth metadata recording trajectory mode, per-stage/group differential
indicators, per-omic δ, and the generator seed
Group A is a random baseline trajectory; group B is a deterministic transform of
group A's per-stage differential indicators selected by trajectory_mode.
from motco.simulations import (
SemiSyntheticTrajectoryParams,
generate_semisynthetic_trajectory,
)
dataset = generate_semisynthetic_trajectory(
SemiSyntheticTrajectoryParams(
seed=99,
trajectory_mode="magnitude",
n_samples=120,
n_stages=3,
group_effect_size=0.2,
group_ratio=0.5,
),
)
sample_metadata = dataset.metadata
truth = dataset.truth
Supported trajectory modes (all governed by the unified group_effect_size knob,
where 0 is the null for every mode):
| Mode |
Injected group-specific pattern |
none |
Identical to baseline; useful for Type I error scenarios |
translation |
Constant observed-space location offset |
magnitude |
Scales δ (size); see the magnitude_kind table below |
orientation |
One global per-omic feature permutation (rotation) |
shape |
Permutes interior stages only (bend); requires at least three stages |
magnitude_kind selects how magnitude scales group B (the motco simulate --magnitude-kind
flag has the same default):
magnitude_kind |
Construction |
'joint' (default) |
Scales every omic's δ (methylation, expression, proteomics) by the same 1 + e. Size-pure in the joint standardized space, not only within each block. This is the construction the Phase 5 magnitude re-measurement validated. |
'all' |
Scales the methylation δ alone at every stage. Size-pure within the methylation block but rotates the concatenated trajectory. |
'extremes' |
Scales only the endpoint stages' methylation indicators. |
'extremes' and the shape mode's shape_kind='magnitude' remain methylation-only scalings; the
joint default does not extend to them.
Explicit group × stage sizes. group_stage_sizes replaces proportional sizing
(n_samples × stage_sample_prop × group_ratio) with an exact table: one tuple of per-stage
counts per group, rows ordered as group_labels and columns as stages. Every group-stage cell gets
exactly the stated number of samples, and truth metadata records the realized table under
truth["group_stage_sizes"]. Use it when the imbalance varies by stage, as in a real cohort:
SemiSyntheticTrajectoryParams(
seed=1,
n_stages=3,
group_stage_sizes=((11, 10, 28), (9, 10, 12)), # 80 samples, the SEA-AD MTG design
)
The table must have two rows of n_stages integers, each at least 1. Setting it together with a
non-default n_samples, stage_sample_prop or group_ratio is an error that names the
conflicting settings. When the table is absent (the default None), generation, truth metadata,
cell ids and parameter signatures are byte-identical to the generator before the field existed:
hash payloads omit the field while it is None.
Migration. Before default-joint-magnitude-construction the default was 'all'. Code that
relied on the default and wants the old construction must pass magnitude_kind="all" (or
--magnitude-kind all). The committed historical study configs under
examples/trajectory_power_study/ already pin "all", so their cells, signatures and datasets are
unchanged.
Evaluation harness
The evaluation harness runs one SemiSyntheticTrajectoryDataset through MOTCO integration and trajectory testing. It is the per-replicate layer used before larger Type I error or power grids.
from motco.simulations import (
SimulationEvaluationParams,
evaluate_semisynthetic_trajectory,
)
evaluation = evaluate_semisynthetic_trajectory(
dataset,
SimulationEvaluationParams(
integration_method="concat",
permutations=0,
),
)
observed_delta = evaluation.pair_statistics["delta"]
truth = evaluation.truth_metadata
Supported integration methods:
| Method |
Behavior |
concat |
Column-binds methylation, expression, and proteomics matrices after deterministic per-feature standardization by default |
snf |
Builds per-omic affinity matrices, fuses them with SNF, and uses spectral embedding as the trajectory outcome matrix |
Block selection. integration_params["layers"] names the omic blocks the evaluation measures,
for every integration method: e.g. {"layers": ["methylation", "expression"]}. All three blocks are
still generated; only the selected ones enter pooled preprocessing, integration, every
realized-geometry checkpoint (per-block scopes for the selected blocks only, joint scopes over them)
and orientation attribution. The selection is normalized to canonical order (methylation,
expression, proteomics) and must be non-empty, duplicate-free and known; integration metadata
records it as integration_params.layers. When the key is absent, all three blocks are measured and
the evaluation result is byte-identical to the harness before the key existed (no layers metadata
entry is added).
Orientation estimator. orientation_estimator selects how each group's orientation is measured
for the observed angle, every RRPP draw and every realized-geometry checkpoint angle. The default
is "net_displacement" (the unit vector from the first to the last stage mean), adopted after the
net-displacement power comparison; it is part of the
cell id and parameter signature and is recorded in runtime_metadata["orientation_estimator"]. Pass
orientation_estimator="pc1" to opt out to the principal-axis estimator: "pc1" normalizes to
None, which is left out of every hash payload and record, so a PC1 evaluation is byte-identical to
the harness before the selection existed. Every committed study config and historical diagnostic pins
"pc1". Only angle depends on the selection. The stats primitives (estimate_difference, RRPP)
keep "pc1" as their default.
Set permutations=0 for observed statistics only. When permutations > 0, the harness runs RRPP and computes upper-tail empirical p-values with plus-one correction:
p = (1 + count(null >= observed)) / (1 + n_permutations)
The result includes observed delta, angle, and shape matrices, scalar two-group pair statistics, optional p-values, latent matrix metadata, generator truth metadata, runtime metadata, group/stage levels, and the trajectory contrast. Shape pair statistics and p-values are reported as unavailable for datasets with fewer than three stages.
Realized-geometry decomposition
Generated semi-synthetic datasets also expose exact group-stage population means in integration units: methylation M-values and native expression/proteomics units. Evaluation uses one pooled fitted preprocessing artifact for both diagnostics and PLS, then reports:
| Checkpoint |
Scopes |
Interpretation |
population_native |
each omic |
Exact construction and biological-cascade geometry without sampling noise |
population_standardized |
each omic and joint |
Exact means under the replicate's pooled fitted feature scaling |
observed_standardized |
each omic and joint |
Geometry actually presented to integration, including finite sampling |
pls_latent |
joint |
Geometry retained in the fitted production PLS representation |
Every applicable scope records both group path lengths, pairwise delta, angle, and shape, and explicit availability flags. Joint native-space geometry is intentionally omitted because the unstandardized omic blocks use incompatible units. Raw distances across feature and latent checkpoints are not scale-equivalent; use the decomposition to locate changes in behavior rather than subtracting distances across spaces.
See the Phase 2 findings for the revised orientation and shape construction contracts and the matched-seed characterization command.
Grid orchestration
The grid orchestration layer enumerates parameter cells, runs local replicates through the evaluation harness, persists one JSONL row per replicate, resumes completed work, and summarizes rejection rates for Type I error or power studies.
from pathlib import Path
from motco.simulations import (
SemiSyntheticTrajectoryParams,
SimulationEvaluationParams,
SimulationRunConfig,
enumerate_type_i_grid,
run_simulation_grid,
summarize_rejection_rates,
)
grid = enumerate_type_i_grid(
baseline_generator_params=SemiSyntheticTrajectoryParams(seed=2, n_samples=60),
evaluation_params=SimulationEvaluationParams(integration_method="concat", permutations=99),
axes={
"generator.n_samples": [60, 120],
"generator.group_ratio": [0.5, 0.7],
},
n_replicates=3,
base_seed=2026,
)
records = run_simulation_grid(
grid,
config=SimulationRunConfig(output_path=Path("simulation-results.jsonl")),
)
summaries = summarize_rejection_rates(records, alpha=0.05)
Each SimulationCell stores a stable cell_id, phase, SemiSyntheticTrajectoryParams, SimulationEvaluationParams, replicate count, base seed, and metadata such as the varied axis. Axis names use explicit namespaces: generator.<field> or evaluation.<field>.
Initial persistence is JSON Lines. Each row records cell and replicate IDs, deterministic seeds, a parameter signature, status, p-values, pair statistics, realized-geometry diagnostics, truth metadata, runtime metadata, cell metadata, and optional error details. The diagnostic schema version participates in the parameter signature. With resume=True, completed rows with matching parameter signatures are skipped. A matching cell/replicate with a different parameter signature raises unless overwrite=True; legacy rows remain readable but have no realized-geometry payload.
summarize_rejection_rates groups completed replicate rows by cell and statistic, then reports available replicate count, rejection count, rejection rate, Monte Carlo standard error, and unavailable replicate count. Missing statistics, such as shape p-values for two-stage datasets, remain unavailable rather than being counted as non-significant.
API
motco.simulations.InterSIMParams
dataclass
Parameters for the InterSIM R package.
Parameters left as None are omitted from the R call, allowing InterSIM
defaults to apply. The covariance parameters currently support InterSIM's
native "indep" string option or None.
Source code in src/motco/simulations/intersim.py
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71 | @dataclass(frozen=True)
class InterSIMParams:
"""Parameters for the InterSIM R package.
Parameters left as ``None`` are omitted from the R call, allowing InterSIM
defaults to apply. The covariance parameters currently support InterSIM's
native ``"indep"`` string option or ``None``.
"""
seed: int
n_sample: int | None = None
cluster_sample_prop: tuple[float, ...] | None = None
delta_methyl: float | None = None
delta_expr: float | None = None
delta_protein: float | None = None
p_dmp: float | None = None
p_deg: float | None = None
p_dep: float | None = None
sigma_methyl: str | None = None
sigma_expr: str | None = None
sigma_protein: str | None = None
cor_methyl_expr: float | None = None
cor_expr_protein: float | None = None
|
motco.simulations.InterSIMResult
dataclass
Normalized InterSIM simulation output.
Source code in src/motco/simulations/intersim.py
74
75
76
77
78
79
80
81
82
83 | @dataclass(frozen=True)
class InterSIMResult:
"""Normalized InterSIM simulation output."""
methylation: pd.DataFrame
expression: pd.DataFrame
proteomics: pd.DataFrame
sample_ids: pd.Index
clusters: pd.Series
metadata: dict[str, Any] = field(default_factory=dict)
|
motco.simulations.InterSIMAvailability
dataclass
Availability result for the external R InterSIM dependency.
Source code in src/motco/simulations/intersim.py
| @dataclass(frozen=True)
class InterSIMAvailability:
"""Availability result for the external R InterSIM dependency."""
available: bool
message: str
rscript_path: str | None = None
|
motco.simulations.SemiSyntheticTrajectoryParams
dataclass
Parameters for semi-synthetic trajectory generation.
group_ratio is the proportion assigned to the first label in
group_labels within every stage. group_effect_size is the unified
effect knob described in the module docstring (0 is null for all modes).
p_dmp is the per-stage probability that a methylation feature is
differential (InterSIM's p.DMP); expression/protein indicators are
derived from it via the cross-omic maps. baseline_continuity (ρ) is the
baseline stage-program continuity described in the module docstring: 0 draws
each stage's program independently, higher values make a CpG's differential
status persist across adjacent stages with correlation ρ while holding every
stage's marginal at p_dmp. delta_* are the per-omic
mean-shift sizes (InterSIM's delta.*). shape_kind selects the
single-interior-stage perturbation used by shape. magnitude_kind
selects whether magnitude scales every omic's effect size together
(joint, the default), group B's methylation effect alone at all
stages (all, a uniform methylation δ scale), or the methylation effect
only at the extreme stages (extremes, first and last) leaving interior
stages at the baseline effect.
group_stage_sizes is the explicit alternative to n_samples /
stage_sample_prop / group_ratio: one tuple of per-stage counts per
group, rows ordered as group_labels and columns as stages, realized
exactly. It may not be combined with a non-default proportional setting.
None (the default) keeps proportional sizing, byte-identical to the
generator before the table existed.
surgery_censoring is the policy for pool-limited surgeries described in
the module docstring: "error" (default) refuses to realize a partial
surgery, "clamp" realizes the largest surgery the pool allows. Use
"clamp" only to reproduce a historical run or for exploratory work where
a partial surgery is acceptable — under it the requested effect no longer
labels the realized construction.
Source code in src/motco/simulations/semisynthetic.py
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153 | @dataclass(frozen=True)
class SemiSyntheticTrajectoryParams:
"""Parameters for semi-synthetic trajectory generation.
``group_ratio`` is the proportion assigned to the first label in
``group_labels`` within every stage. ``group_effect_size`` is the unified
effect knob described in the module docstring (``0`` is null for all modes).
``p_dmp`` is the per-stage probability that a methylation feature is
differential (InterSIM's ``p.DMP``); expression/protein indicators are
derived from it via the cross-omic maps. ``baseline_continuity`` (ρ) is the
baseline stage-program continuity described in the module docstring: 0 draws
each stage's program independently, higher values make a CpG's differential
status persist across adjacent stages with correlation ρ while holding every
stage's marginal at ``p_dmp``. ``delta_*`` are the per-omic
mean-shift sizes (InterSIM's ``delta.*``). ``shape_kind`` selects the
single-interior-stage perturbation used by ``shape``. ``magnitude_kind``
selects whether ``magnitude`` scales every omic's effect size together
(``joint``, the default), group B's methylation effect alone at *all*
stages (``all``, a uniform methylation δ scale), or the methylation effect
only at the *extreme* stages (``extremes``, first and last) leaving interior
stages at the baseline effect.
``group_stage_sizes`` is the explicit alternative to ``n_samples`` /
``stage_sample_prop`` / ``group_ratio``: one tuple of per-stage counts per
group, rows ordered as ``group_labels`` and columns as stages, realized
exactly. It may not be combined with a non-default proportional setting.
``None`` (the default) keeps proportional sizing, byte-identical to the
generator before the table existed.
``surgery_censoring`` is the policy for pool-limited surgeries described in
the module docstring: ``"error"`` (default) refuses to realize a partial
surgery, ``"clamp"`` realizes the largest surgery the pool allows. Use
``"clamp"`` only to reproduce a historical run or for exploratory work where
a partial surgery is acceptable — under it the requested effect no longer
labels the realized construction.
"""
seed: int
trajectory_mode: TrajectoryMode = "none"
n_samples: int = 600
n_stages: int = 3
group_effect_size: float = 0.0
group_ratio: float = 0.5
group_labels: tuple[str, str] = ("A", "B")
p_dmp: float = 0.2
baseline_continuity: float = 0.0
delta_methyl: float = 2.0
delta_expr: float = 2.0
delta_protein: float = 2.0
shape_kind: ShapeKind = "relocate"
magnitude_kind: MagnitudeKind = "joint"
surgery_censoring: SurgeryCensoring = "error"
stage_sample_prop: tuple[float, ...] | None = None
group_stage_sizes: tuple[tuple[int, ...], tuple[int, ...]] | None = None
|
motco.simulations.SemiSyntheticTrajectoryDataset
dataclass
MOTCO-ready semi-synthetic trajectory dataset.
Source code in src/motco/simulations/semisynthetic.py
171
172
173
174
175
176
177
178
179
180 | @dataclass(frozen=True)
class SemiSyntheticTrajectoryDataset:
"""MOTCO-ready semi-synthetic trajectory dataset."""
methylation: pd.DataFrame
expression: pd.DataFrame
proteomics: pd.DataFrame
metadata: pd.DataFrame
truth: dict[str, Any] = field(default_factory=dict)
population_trajectories: PopulationTrajectories | None = None
|
motco.simulations.SimulationEvaluationParams
dataclass
Parameters for evaluating one semi-synthetic trajectory dataset.
Source code in src/motco/simulations/evaluation.py
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145 | @dataclass(frozen=True)
class SimulationEvaluationParams:
"""Parameters for evaluating one semi-synthetic trajectory dataset."""
integration_method: IntegrationMethod = "concat"
integration_params: Mapping[str, Any] = field(default_factory=dict)
permutations: int = 0
n_jobs: int | None = 1
seed: int | None = None
progress: bool = False
include_null_distributions: bool = False
group_col: str = "group"
stage_col: str = "stage"
attribution: AttributionDiagnosticSettings = field(default_factory=AttributionDiagnosticSettings)
# Orientation estimator for the observed ``angle``, every RRPP draw and every
# realized-geometry checkpoint. The default is net displacement; it is hashed
# into cell ids and signatures and recorded in ``runtime_metadata``. ``None``
# is the hash-omitted spelling of PC1 (``grid._SIGNATURE_OMIT_WHEN_NONE``) and
# an explicit ``"pc1"`` normalizes to it, so PC1 cell ids, signatures and
# records equal those produced before the field existed.
orientation_estimator: str | None = "net_displacement"
def __post_init__(self) -> None:
estimator = self.orientation_estimator
if estimator is None:
return
if estimator not in EVALUATION_ORIENTATION_ESTIMATORS:
raise SimulationEvaluationError(
f"Unknown orientation_estimator {estimator!r}; allowed: "
f"{list(EVALUATION_ORIENTATION_ESTIMATORS)}."
)
if estimator == "pc1":
object.__setattr__(self, "orientation_estimator", None)
@property
def resolved_orientation_estimator(self) -> str:
"""The estimator name the statistics are computed with."""
return self.orientation_estimator or "pc1"
|
resolved_orientation_estimator
property
The estimator name the statistics are computed with.
motco.simulations.SimulationEvaluationResult
dataclass
Result from evaluating one semi-synthetic trajectory dataset.
Source code in src/motco/simulations/evaluation.py
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227 | @dataclass(frozen=True)
class SimulationEvaluationResult:
"""Result from evaluating one semi-synthetic trajectory dataset."""
observed_deltas: np.ndarray
observed_angles: np.ndarray
observed_shapes: np.ndarray
pair_statistics: dict[str, float]
p_values: dict[str, float]
latent_matrix_metadata: dict[str, Any]
truth_metadata: dict[str, Any]
runtime_metadata: dict[str, Any]
evaluation_params: SimulationEvaluationParams
group_levels: list[str]
stage_levels: list[str]
contrast: list[list[int]]
null_distributions: dict[str, list[float]] | None = None
realized_geometry: dict[str, Any] = field(default_factory=dict)
attribution_diagnostics: dict[str, Any] = field(default_factory=dict)
# Compact per-statistic description of the permutation null. Produced on
# every RRPP run — unlike ``null_distributions``, which stays opt-in — so a
# record can locate its own observed statistic against its own critical
# value without retaining the full draw vectors.
null_summary: dict[str, dict[str, float]] = field(default_factory=dict)
# Eigenspectrum of the centered stage-mean configuration in the *evaluated
# latent space*, pooled and per contrast group, plus — whenever permutations
# ran — a summary of the pooled relative eigengap over the permutation draws.
# A recorded covariate: the geometry audit
# (``docs/reports/geometry-audit-2026-09-01.md``) found the relative eigengap
# predicts how wide a replicate's own ``angle`` permutation null will be, and
# therefore whether its orientation is resolvable at all. It qualifies
# reporting; it enters no statistic and no decision rule.
config_spectrum: dict[str, Any] = field(default_factory=dict)
|
motco.simulations.SimulationTrajectoryDesign
dataclass
Trajectory design objects derived from sample metadata.
Source code in src/motco/simulations/evaluation.py
183
184
185
186
187
188
189
190
191
192 | @dataclass(frozen=True)
class SimulationTrajectoryDesign:
"""Trajectory design objects derived from sample metadata."""
model_full: np.ndarray
model_reduced: np.ndarray
ls_means: np.ndarray
contrast: list[list[int]]
group_levels: list[str]
stage_levels: list[str]
|
motco.simulations.LatentIntegrationResult
dataclass
Integrated latent/outcome matrix and metadata.
artifacts carries the ephemeral fitted objects for the PLS path and is
never serialized; only metadata reaches the persisted record.
Source code in src/motco/simulations/evaluation.py
170
171
172
173
174
175
176
177
178
179
180 | @dataclass(frozen=True)
class LatentIntegrationResult:
"""Integrated latent/outcome matrix and metadata.
``artifacts`` carries the ephemeral fitted objects for the PLS path and is
never serialized; only ``metadata`` reaches the persisted record.
"""
matrix: pd.DataFrame
metadata: dict[str, Any]
artifacts: PLSIntegrationArtifacts | None = None
|
motco.simulations.SimulationCell
dataclass
One simulation parameter cell with one or more replicates.
Source code in src/motco/simulations/grid.py
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61 | @dataclass(frozen=True)
class SimulationCell:
"""One simulation parameter cell with one or more replicates."""
cell_id: str
phase: str
generator_params: SemiSyntheticTrajectoryParams
evaluation_params: SimulationEvaluationParams
n_replicates: int = 1
base_seed: int = 0
metadata: Mapping[str, Any] = field(default_factory=dict)
def __post_init__(self) -> None:
if self.n_replicates < 1:
raise SimulationGridError("n_replicates must be at least 1.")
if not self.cell_id:
raise SimulationGridError("cell_id must be a non-empty string.")
if not self.phase:
raise SimulationGridError("phase must be a non-empty string.")
|
motco.simulations.SimulationGrid
dataclass
Collection of simulation cells.
Source code in src/motco/simulations/grid.py
64
65
66
67
68
69
70
71
72
73
74
75 | @dataclass(frozen=True)
class SimulationGrid:
"""Collection of simulation cells."""
cells: tuple[SimulationCell, ...]
metadata: Mapping[str, Any] = field(default_factory=dict)
def __post_init__(self) -> None:
ids = [cell.cell_id for cell in self.cells]
duplicates = sorted({cell_id for cell_id in ids if ids.count(cell_id) > 1})
if duplicates:
raise SimulationGridError(f"cell_id values must be unique; duplicates: {duplicates}.")
|
motco.simulations.SimulationReplicateResult
dataclass
One persisted row for a simulation cell replicate.
Source code in src/motco/simulations/grid.py
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121 | @dataclass(frozen=True)
class SimulationReplicateResult:
"""One persisted row for a simulation cell replicate."""
cell_id: str
phase: str
replicate_index: int
replicate_seed: int
generator_seed: int
evaluation_seed: int | None
parameter_signature: str
status: Literal["completed", "failed"]
p_values: dict[str, float | None] = field(default_factory=dict)
pair_statistics: dict[str, float | None] = field(default_factory=dict)
realized_geometry: dict[str, Any] = field(default_factory=dict)
truth_metadata: dict[str, Any] = field(default_factory=dict)
runtime_metadata: dict[str, Any] = field(default_factory=dict)
cell_metadata: dict[str, Any] = field(default_factory=dict)
# Additive Phase 4 fields. They default to empty so records written before
# they existed still load, and so a failed diagnostic is distinguishable
# from one that was never requested.
integration_metadata: dict[str, Any] = field(default_factory=dict)
attribution_status: str = "not_requested"
attribution_diagnostics: dict[str, Any] = field(default_factory=dict)
# Per-statistic description of this replicate's own permutation null. Empty
# for records written before the field existed and for evaluations that ran
# no permutations; ``runtime_metadata["permutations"]`` separates the two.
# Deliberately absent from ``parameter_signature`` — it summarizes draws that
# already happened and changes no generation, integration, or permutation
# behavior, so pre-change shards stay resumable.
null_summary: dict[str, dict[str, float]] = field(default_factory=dict)
# Latent stage-mean configuration spectrum (pooled and per group) and, when
# permutations ran, the pooled eigengap over those draws. Empty for records
# written before the field existed; a *recorded* degenerate spectrum is a
# populated block whose ``relative_eigengap`` is ``None``, which is how the
# two are told apart. Unlike ``null_summary`` this **is** in
# ``parameter_signature``: a resumed pre-change shard would otherwise produce
# a merged set where only some records carry the covariate, and a missing
# covariate is indistinguishable from a degenerate one at analysis time.
config_spectrum: dict[str, Any] = field(default_factory=dict)
diagnostic_error_type: str | None = None
diagnostic_error_message: str | None = None
error_type: str | None = None
error_message: str | None = None
|
motco.simulations.SimulationRunConfig
dataclass
Runtime options for local grid execution.
Source code in src/motco/simulations/grid.py
178
179
180
181
182
183
184
185
186 | @dataclass(frozen=True)
class SimulationRunConfig:
"""Runtime options for local grid execution."""
output_path: Path | None = None
persistence_format: PersistenceFormat = "jsonl"
resume: bool = True
overwrite: bool = False
error_policy: ErrorPolicy = "raise"
|
motco.simulations.SimulationSummaryResult
dataclass
Rejection-rate summary for one cell and statistic.
Source code in src/motco/simulations/grid.py
124
125
126
127
128
129
130
131
132
133
134
135
136
137 | @dataclass(frozen=True)
class SimulationSummaryResult:
"""Rejection-rate summary for one cell and statistic."""
cell_id: str
phase: str
statistic: str
alpha: float
completed_replicates: int
available_replicates: int
rejected_replicates: int
rejection_rate: float | None
monte_carlo_se: float | None
unavailable_replicates: int
|
motco.simulations.check_intersim_available(rscript='Rscript')
Check whether Rscript and the R InterSIM package are available.
Source code in src/motco/simulations/intersim.py
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131 | def check_intersim_available(rscript: str = "Rscript") -> InterSIMAvailability:
"""Check whether Rscript and the R InterSIM package are available."""
rscript_path = which(rscript)
if rscript_path is None:
return InterSIMAvailability(
available=False,
message=f"Rscript dependency is missing or not on PATH: {rscript}",
)
cmd = [
rscript_path,
"-e",
'quit(status = if (requireNamespace("InterSIM", quietly = TRUE)) 0 else 1)',
]
proc = subprocess.run(cmd, capture_output=True, text=True, check=False)
if proc.returncode != 0:
details = _format_process_output(proc)
return InterSIMAvailability(
available=False,
message="R package InterSIM is not installed or cannot be loaded." + details,
rscript_path=rscript_path,
)
return InterSIMAvailability(
available=True,
message="Rscript and R package InterSIM are available.",
rscript_path=rscript_path,
)
|
motco.simulations.run_intersim(params, *, rscript='Rscript', check_dependency=True)
Invoke R InterSIM and return aligned omics matrices plus cluster metadata.
Source code in src/motco/simulations/intersim.py
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163 | def run_intersim(
params: InterSIMParams,
*,
rscript: str = "Rscript",
check_dependency: bool = True,
) -> InterSIMResult:
"""Invoke R InterSIM and return aligned omics matrices plus cluster metadata."""
if check_dependency:
availability = check_intersim_available(rscript)
if not availability.available:
raise InterSIMDependencyError(availability.message)
rscript_cmd = availability.rscript_path or rscript
else:
rscript_cmd = which(rscript) or rscript
with tempfile.TemporaryDirectory(prefix="motco-intersim-") as tmp:
output_dir = Path(tmp)
proc = subprocess.run(
_build_rscript_command(params, output_dir=output_dir, rscript=rscript_cmd),
capture_output=True,
text=True,
check=False,
)
if proc.returncode != 0:
raise InterSIMRuntimeError(
f"InterSIM R process failed with exit code {proc.returncode}."
f"{_format_process_output(proc)}"
)
return _load_result(output_dir, params=params)
|
motco.simulations.generate_semisynthetic_trajectory(params, *, reference=None)
Generate a semi-synthetic trajectory dataset using the numpy generator.
Source code in src/motco/simulations/semisynthetic.py
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237 | def generate_semisynthetic_trajectory(
params: SemiSyntheticTrajectoryParams,
*,
reference: IntersimReference | None = None,
) -> SemiSyntheticTrajectoryDataset:
"""Generate a semi-synthetic trajectory dataset using the numpy generator."""
_validate_params(params)
ref = reference if reference is not None else load_reference()
rng = np.random.default_rng(params.seed)
if params.group_stage_sizes is not None:
group_a_sizes, group_b_sizes = (list(row) for row in params.group_stage_sizes)
else:
stage_sizes = _stage_sizes(params)
group_a_sizes, group_b_sizes = _group_stage_sizes(stage_sizes, params)
methyl_a = _baseline_methyl(rng, ref, params)
methyl_b, deltas_b, transform_meta = _transform_group_b(rng, ref, params, methyl_a)
deltas_a = (params.delta_methyl, params.delta_expr, params.delta_protein)
indicators_a = _derive_group(methyl_a, ref)
indicators_b = _derive_group(methyl_b, ref)
gen_a = _generate_group(rng, ref, indicators_a, deltas_a, group_a_sizes)
gen_b = _generate_group(rng, ref, indicators_b, deltas_b, group_b_sizes)
population_trajectories = _build_population_trajectories(
ref, params, indicators_a, indicators_b, deltas_a, deltas_b
)
methylation, expression, proteomics, metadata = _assemble(
ref, params, gen_a, gen_b, group_a_sizes, group_b_sizes
)
truth = _build_truth(params, indicators_a, indicators_b, deltas_a, deltas_b, transform_meta)
if params.group_stage_sizes is not None:
truth["group_stage_sizes"] = {
params.group_labels[0]: list(group_a_sizes),
params.group_labels[1]: list(group_b_sizes),
}
return SemiSyntheticTrajectoryDataset(
methylation=methylation,
expression=expression,
proteomics=proteomics,
metadata=metadata,
truth=truth,
population_trajectories=population_trajectories,
)
|
motco.simulations.integrate_semisynthetic_dataset(dataset, params, *, preprocessor=None, observed_blocks=None)
Construct the molecular latent space from aligned omics layers.
snf and pls build genuine latent spaces (the measurement substrate
for trajectory geometry); concat is a standardized-feature-concatenation
baseline, not a constructed latent space. See the module docstring for the
latent-space architecture.
Source code in src/motco/simulations/evaluation.py
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468 | def integrate_semisynthetic_dataset(
dataset: SemiSyntheticTrajectoryDataset,
params: SimulationEvaluationParams,
*,
preprocessor: FittedOmicsPreprocessor | None = None,
observed_blocks: Mapping[OmicsLayer, pd.DataFrame] | None = None,
) -> LatentIntegrationResult:
"""Construct the molecular latent space from aligned omics layers.
``snf`` and ``pls`` build genuine latent spaces (the measurement substrate
for trajectory geometry); ``concat`` is a standardized-feature-concatenation
baseline, not a constructed latent space. See the module docstring for the
latent-space architecture.
"""
_validate_dataset(dataset, params.group_col, params.stage_col)
layers = _measured_layers(params.integration_params)
if preprocessor is None:
preprocessor = fit_omics_preprocessor(dataset, layers)
else:
try:
preprocessor.require_layers(layers)
except ValueError as exc:
raise SimulationEvaluationError(str(exc)) from exc
if observed_blocks is None:
observed_blocks = preprocessor.transform_dataset(dataset)
method = params.integration_method
if method == "concat":
return _concat_integration(dataset, params.integration_params, observed_blocks)
if method == "snf":
return _snf_integration(dataset, params.integration_params)
if method == "pls":
return _pls_integration(
dataset,
params.integration_params,
stage_col=params.stage_col,
observed_blocks=observed_blocks,
preprocessor=preprocessor,
)
if method == "pca":
return _pca_integration(dataset, params.integration_params, observed_blocks)
if method == "pls_within_group":
return _pls_within_group_integration(
dataset,
params.integration_params,
group_col=params.group_col,
stage_col=params.stage_col,
observed_blocks=observed_blocks,
)
raise SimulationEvaluationError(f"Unsupported integration_method: {method!r}.")
|
motco.simulations.build_simulation_trajectory_design(metadata, *, group_col='group', stage_col='stage')
Build full/reduced model matrices, LS means, and two-group contrast.
Source code in src/motco/simulations/evaluation.py
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505 | def build_simulation_trajectory_design(
metadata: pd.DataFrame,
*,
group_col: str = "group",
stage_col: str = "stage",
) -> SimulationTrajectoryDesign:
"""Build full/reduced model matrices, LS means, and two-group contrast."""
_validate_metadata(metadata, group_col, stage_col)
design_frame = metadata[[group_col, stage_col]].reset_index(drop=True).copy()
design_frame[group_col] = design_frame[group_col].astype(str)
design_frame[stage_col] = design_frame[stage_col].astype(str)
group_levels = _sort_levels(pd.unique(design_frame[group_col]).tolist())
stage_levels = _sort_levels(pd.unique(design_frame[stage_col]).tolist())
if len(group_levels) != 2:
raise SimulationEvaluationError(f"Expected exactly two groups, found {len(group_levels)}: {group_levels}.")
if len(stage_levels) < 2:
raise SimulationEvaluationError(
f"Expected at least two stages, found {len(stage_levels)}: {stage_levels}."
)
_validate_group_stage_combinations(design_frame, group_col, stage_col, group_levels, stage_levels)
model_full = get_model_matrix(design_frame, group_col=group_col, level_col=stage_col, full=True)
model_reduced = get_model_matrix(design_frame, group_col=group_col, level_col=stage_col, full=False)
ls_means = build_ls_means(group_levels, stage_levels, full=True)
n_stages = len(stage_levels)
contrast = [[group_index * n_stages + stage_index for stage_index in range(n_stages)] for group_index in range(2)]
return SimulationTrajectoryDesign(
model_full=model_full,
model_reduced=model_reduced,
ls_means=ls_means,
contrast=contrast,
group_levels=group_levels,
stage_levels=stage_levels,
)
|
motco.simulations.evaluate_semisynthetic_trajectory(dataset, params=None, *, latent=None)
Evaluate one semi-synthetic trajectory dataset through MOTCO routines.
latent optionally supplies an integration already built for this
dataset with the same integration_method and integration_params
(from :func:integrate_semisynthetic_dataset), so evaluations that differ
only in measurement settings — the orientation estimator, say — share one
space fit. Omitted (the default), the integration is built here.
Source code in src/motco/simulations/evaluation.py
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371 | def evaluate_semisynthetic_trajectory(
dataset: SemiSyntheticTrajectoryDataset,
params: SimulationEvaluationParams | None = None,
*,
latent: LatentIntegrationResult | None = None,
) -> SimulationEvaluationResult:
"""Evaluate one semi-synthetic trajectory dataset through MOTCO routines.
``latent`` optionally supplies an integration already built for this
dataset with the same ``integration_method`` and ``integration_params``
(from :func:`integrate_semisynthetic_dataset`), so evaluations that differ
only in measurement settings — the orientation estimator, say — share one
space fit. Omitted (the default), the integration is built here.
"""
params = params or SimulationEvaluationParams()
_validate_evaluation_params(params)
start = perf_counter()
preprocessor = fit_omics_preprocessor(dataset, _measured_layers(params.integration_params))
observed_blocks = preprocessor.transform_dataset(dataset)
if latent is None:
latent = integrate_semisynthetic_dataset(
dataset,
params,
preprocessor=preprocessor,
observed_blocks=observed_blocks,
)
elif latent.metadata.get("integration_method") != params.integration_method:
raise SimulationEvaluationError(
f"Supplied latent integration was built with {latent.metadata.get('integration_method')!r}, "
f"not integration_method={params.integration_method!r}."
)
design = build_simulation_trajectory_design(
dataset.metadata,
group_col=params.group_col,
stage_col=params.stage_col,
)
observed_deltas, observed_angles, observed_shapes, observed_spectra = estimate_difference(
Y=latent.matrix,
model_matrix=design.model_full,
LS_means=design.ls_means,
contrast=design.contrast,
return_spectra=True,
orientation_estimator=params.resolved_orientation_estimator,
)
config_spectrum = _config_spectrum_block(observed_spectra, design.group_levels)
shape_available = len(design.stage_levels) >= 3
pair_statistics = _extract_pair_statistics(
observed_deltas,
observed_angles,
observed_shapes,
shape_available=shape_available,
)
p_values: dict[str, float] = {}
null_distributions: dict[str, list[float]] | None = None
null_summary: dict[str, dict[str, float]] = {}
if params.permutations > 0:
dist_delta, dist_angle, dist_shape, dist_eigengap = RRPP(
latent.matrix,
design.model_full,
design.model_reduced,
design.ls_means,
design.contrast,
permutations=params.permutations,
n_jobs=params.n_jobs,
progress=params.progress,
seed=params.seed,
return_eigengaps=True,
orientation_estimator=params.resolved_orientation_estimator,
)
null_values = _extract_null_distributions(
dist_delta,
dist_angle,
dist_shape,
shape_available=shape_available,
)
p_values = {
statistic: _empirical_p_value(null_values[statistic], observed)
for statistic, observed in pair_statistics.items()
if statistic in null_values and np.isfinite(observed)
}
null_summary = _summarize_null_distributions(null_values)
# Only the summary survives: the per-permutation eigengaps are discarded
# here, exactly as the per-permutation statistics are.
config_spectrum["permutation_pooled_eigengap"] = _summarize_permutation_eigengaps(
dist_eigengap
)
if params.include_null_distributions:
null_distributions = null_values
runtime_seconds = perf_counter() - start
runtime_metadata: dict[str, Any] = {
"runtime_seconds": runtime_seconds,
"permutations": params.permutations,
"n_jobs": params.n_jobs,
"seed": params.seed,
"progress": params.progress,
"shape_available": shape_available,
}
if params.orientation_estimator is not None:
runtime_metadata["orientation_estimator"] = params.orientation_estimator
latent_checkpoint = _LATENT_CHECKPOINTS.get(params.integration_method)
realized_geometry = calculate_realized_geometry(
dataset,
preprocessor=preprocessor,
observed_blocks=observed_blocks,
latent_matrix=latent.matrix if latent_checkpoint is not None else None,
group_col=params.group_col,
stage_col=params.stage_col,
latent_checkpoint=latent_checkpoint or "pls_latent",
orientation_estimator=params.resolved_orientation_estimator,
).to_dict()
# The fitted estimator and standardized joint matrix live only inside this
# call: attribution consumes them here, and only the compact JSON-safe
# record survives into the result.
attribution_diagnostics = _attribution_diagnostics(dataset, params, latent)
return SimulationEvaluationResult(
observed_deltas=observed_deltas,
observed_angles=observed_angles,
observed_shapes=observed_shapes,
pair_statistics=pair_statistics,
p_values=p_values,
latent_matrix_metadata=latent.metadata,
# Drop the raw indicator arrays — they are not JSON-serializable for the
# study's JSONL persistence; the per-stage/group counts and transform
# notes (which are) are retained for downstream characterization.
truth_metadata={k: v for k, v in dataset.truth.items() if k != "indicators"},
runtime_metadata=runtime_metadata,
evaluation_params=params,
group_levels=design.group_levels,
stage_levels=design.stage_levels,
contrast=design.contrast,
null_distributions=null_distributions,
realized_geometry=realized_geometry,
attribution_diagnostics=attribution_diagnostics,
null_summary=null_summary,
config_spectrum=config_spectrum,
)
|
motco.simulations.enumerate_type_i_grid(*, baseline_generator_params, evaluation_params=None, axes=None, n_replicates=1, base_seed=0)
Enumerate a null Type I grid with baseline plus one-factor-at-a-time axes.
Source code in src/motco/simulations/grid.py
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270 | def enumerate_type_i_grid(
*,
baseline_generator_params: SemiSyntheticTrajectoryParams,
evaluation_params: SimulationEvaluationParams | None = None,
axes: Mapping[str, Sequence[Any]] | None = None,
n_replicates: int = 1,
base_seed: int = 0,
) -> SimulationGrid:
"""Enumerate a null Type I grid with baseline plus one-factor-at-a-time axes."""
null_generator = replace(baseline_generator_params, trajectory_mode="none", group_effect_size=0.0)
cells = [
make_simulation_cell(
phase="type_i_baseline",
generator_params=null_generator,
evaluation_params=evaluation_params,
n_replicates=n_replicates,
base_seed=base_seed,
metadata={"varied_axis": None, "varied_value": None},
)
]
for axis, values in (axes or {}).items():
baseline_value = _get_axis_value(null_generator, evaluation_params, axis)
for value in values:
if _to_jsonable(value) == _to_jsonable(baseline_value):
continue
generator_params, eval_params = _apply_axis_value(
null_generator,
evaluation_params or SimulationEvaluationParams(),
axis,
value,
)
cells.append(
make_simulation_cell(
phase="type_i_ofat",
generator_params=generator_params,
evaluation_params=eval_params,
n_replicates=n_replicates,
base_seed=base_seed,
metadata={"varied_axis": axis, "varied_value": value},
)
)
return SimulationGrid(cells=tuple(cells), metadata={"grid_type": "type_i"})
|
motco.simulations.enumerate_power_grid(*, baseline_generator_params, evaluation_params=None, trajectory_modes, effect_sizes, axes=None, n_replicates=1, base_seed=0, evaluation_resolver=None, cell_metadata=None)
Enumerate power cells from trajectory modes, effect sizes, and optional axes.
evaluation_resolver may specialize the evaluation parameters per cell
(for example, to enable attribution diagnostics only on selected cells).
Because it runs before cell construction, the resolved settings enter the
cell id and parameter signature rather than being applied after the fact.
cell_metadata is merged into every emitted cell's metadata.
Source code in src/motco/simulations/grid.py
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360 | def enumerate_power_grid(
*,
baseline_generator_params: SemiSyntheticTrajectoryParams,
evaluation_params: SimulationEvaluationParams | None = None,
trajectory_modes: Sequence[str],
effect_sizes: Sequence[float],
axes: Mapping[str, Sequence[Any]] | None = None,
n_replicates: int = 1,
base_seed: int = 0,
evaluation_resolver: CellEvaluationResolver | None = None,
cell_metadata: Mapping[str, Any] | None = None,
) -> SimulationGrid:
"""Enumerate power cells from trajectory modes, effect sizes, and optional axes.
``evaluation_resolver`` may specialize the evaluation parameters per cell
(for example, to enable attribution diagnostics only on selected cells).
Because it runs *before* cell construction, the resolved settings enter the
cell id and parameter signature rather than being applied after the fact.
``cell_metadata`` is merged into every emitted cell's metadata.
"""
if not trajectory_modes:
raise SimulationGridError("trajectory_modes must contain at least one mode.")
if not effect_sizes:
raise SimulationGridError("effect_sizes must contain at least one value.")
cells: list[SimulationCell] = []
eval_params = evaluation_params or SimulationEvaluationParams()
extra_metadata = dict(cell_metadata or {})
def resolve(phase: str, mode: str, effect_size: float) -> SimulationEvaluationParams:
if evaluation_resolver is None:
return eval_params
return evaluation_resolver(phase, mode, float(effect_size))
for mode, effect_size in itertools.product(trajectory_modes, effect_sizes):
generator = replace(
baseline_generator_params,
trajectory_mode=mode, # type: ignore[arg-type]
group_effect_size=float(effect_size),
)
primary_eval = resolve("power_primary", mode, float(effect_size))
cells.append(
make_simulation_cell(
phase="power_primary",
generator_params=generator,
evaluation_params=primary_eval,
n_replicates=n_replicates,
base_seed=base_seed,
metadata={
"trajectory_mode": mode,
"effect_size": float(effect_size),
"varied_axis": None,
**extra_metadata,
},
)
)
for axis, values in (axes or {}).items():
ofat_eval = resolve("power_ofat", mode, float(effect_size))
baseline_value = _get_axis_value(generator, ofat_eval, axis)
for value in values:
if _to_jsonable(value) == _to_jsonable(baseline_value):
continue
generator_params, axis_eval_params = _apply_axis_value(
generator,
ofat_eval,
axis,
value,
)
cells.append(
make_simulation_cell(
phase="power_ofat",
generator_params=generator_params,
evaluation_params=axis_eval_params,
n_replicates=n_replicates,
base_seed=base_seed,
metadata={
"trajectory_mode": mode,
"effect_size": float(effect_size),
"varied_axis": axis,
"varied_value": value,
},
)
)
return SimulationGrid(cells=tuple(cells), metadata={"grid_type": "power"})
|
motco.simulations.run_simulation_replicate(cell, replicate_index, *, evaluator=None, error_policy='raise')
Run one replicate, using a fake evaluator in tests or the default harness in production.
Source code in src/motco/simulations/grid.py
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503 | def run_simulation_replicate(
cell: SimulationCell,
replicate_index: int,
*,
evaluator: Evaluator | None = None,
error_policy: ErrorPolicy = "raise",
) -> SimulationReplicateResult:
"""Run one replicate, using a fake evaluator in tests or the default harness in production."""
_validate_error_policy(error_policy)
replicate_seed = derive_replicate_seed(cell, replicate_index)
generator_params = replace(cell.generator_params, seed=replicate_seed)
evaluation_seed = cell.evaluation_params.seed if cell.evaluation_params.seed is not None else replicate_seed
evaluation_params = replace(cell.evaluation_params, seed=evaluation_seed)
evaluator = evaluator or _default_evaluator
start = perf_counter()
try:
result = evaluator(generator_params, evaluation_params)
except Exception as exc:
if error_policy == "raise":
raise
return SimulationReplicateResult(
cell_id=cell.cell_id,
phase=cell.phase,
replicate_index=replicate_index,
replicate_seed=replicate_seed,
generator_seed=generator_params.seed,
evaluation_seed=evaluation_params.seed,
parameter_signature=parameter_signature(cell),
status="failed",
runtime_metadata={"runtime_seconds": perf_counter() - start},
cell_metadata=dict(cell.metadata),
error_type=type(exc).__name__,
error_message=str(exc),
)
attribution = dict(result.attribution_diagnostics or {})
attribution_status = str(attribution.get("status", "not_requested"))
return SimulationReplicateResult(
cell_id=cell.cell_id,
phase=cell.phase,
replicate_index=replicate_index,
replicate_seed=replicate_seed,
generator_seed=generator_params.seed,
evaluation_seed=evaluation_params.seed,
parameter_signature=parameter_signature(cell),
status="completed",
p_values=_finite_mapping(result.p_values),
pair_statistics=_finite_mapping(result.pair_statistics),
realized_geometry=result.realized_geometry,
truth_metadata=result.truth_metadata,
runtime_metadata=result.runtime_metadata,
cell_metadata=dict(cell.metadata),
integration_metadata=dict(result.latent_matrix_metadata),
attribution_status=attribution_status,
attribution_diagnostics=attribution,
null_summary={
statistic: dict(entry) for statistic, entry in (result.null_summary or {}).items()
},
config_spectrum=dict(result.config_spectrum or {}),
diagnostic_error_type=attribution.get("error_type"),
diagnostic_error_message=attribution.get("reason") if attribution_status == "failed" else None,
)
|
motco.simulations.run_simulation_grid(grid, *, config=None, evaluator=None)
Run a grid locally with optional JSONL persistence and resume support.
Source code in src/motco/simulations/grid.py
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546 | def run_simulation_grid(
grid: SimulationGrid,
*,
config: SimulationRunConfig | None = None,
evaluator: Evaluator | None = None,
) -> list[SimulationReplicateResult]:
"""Run a grid locally with optional JSONL persistence and resume support."""
config = config or SimulationRunConfig()
if config.persistence_format != "jsonl":
raise SimulationGridError("Only JSONL persistence is currently supported.")
completed = (
_completed_index(config.output_path)
if config.output_path and config.resume and not config.overwrite
else {}
)
records: list[SimulationReplicateResult] = []
for cell in grid.cells:
signature = parameter_signature(cell)
for replicate_index in range(cell.n_replicates):
key = (cell.cell_id, replicate_index)
existing_signature = completed.get(key)
if existing_signature is not None:
if existing_signature != signature:
if not config.overwrite:
raise SimulationGridError(
f"Existing result for {cell.cell_id} replicate {replicate_index} has a different "
"parameter signature. Set overwrite=True to rerun."
)
else:
continue
record = run_simulation_replicate(
cell,
replicate_index,
evaluator=evaluator,
error_policy=config.error_policy,
)
records.append(record)
if config.output_path is not None:
append_replicate_results(config.output_path, [record])
return records
|
motco.simulations.read_replicate_results(path)
Read replicate records from a JSONL file.
Source code in src/motco/simulations/grid.py
558
559
560
561
562
563
564
565
566
567
568
569 | def read_replicate_results(path: Path) -> list[SimulationReplicateResult]:
"""Read replicate records from a JSONL file."""
if not path.exists():
return []
records = []
with path.open("r", encoding="utf-8") as handle:
for line in handle:
if not line.strip():
continue
records.append(_replicate_result_from_dict(json.loads(line)))
return records
|
motco.simulations.append_replicate_results(path, records)
Append replicate records to a JSONL file.
Source code in src/motco/simulations/grid.py
549
550
551
552
553
554
555 | def append_replicate_results(path: Path, records: Iterable[SimulationReplicateResult]) -> None:
"""Append replicate records to a JSONL file."""
path.parent.mkdir(parents=True, exist_ok=True)
with path.open("a", encoding="utf-8") as handle:
for record in records:
handle.write(json.dumps(_to_jsonable(record), sort_keys=True, allow_nan=False) + "\n")
|
motco.simulations.summarize_rejection_rates(records, *, alpha=0.05, statistics=('delta', 'angle', 'shape'))
Summarize p-value rejection rates by cell and statistic.
Source code in src/motco/simulations/grid.py
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610 | def summarize_rejection_rates(
records: Sequence[SimulationReplicateResult],
*,
alpha: float = 0.05,
statistics: Sequence[str] = ("delta", "angle", "shape"),
) -> list[SimulationSummaryResult]:
"""Summarize p-value rejection rates by cell and statistic."""
if not (0 < alpha < 1):
raise SimulationGridError("alpha must be between 0 and 1.")
groups: dict[tuple[str, str], list[SimulationReplicateResult]] = {}
for record in records:
if record.status != "completed":
continue
groups.setdefault((record.cell_id, record.phase), []).append(record)
summaries: list[SimulationSummaryResult] = []
for (cell_id, phase), group_records in sorted(groups.items()):
for statistic in statistics:
p_values = [record.p_values.get(statistic) for record in group_records]
available = [float(p) for p in p_values if p is not None and math.isfinite(float(p))]
rejected = sum(p < alpha for p in available)
rate = rejected / len(available) if available else None
se = math.sqrt(rate * (1.0 - rate) / len(available)) if rate is not None else None
summaries.append(
SimulationSummaryResult(
cell_id=cell_id,
phase=phase,
statistic=statistic,
alpha=alpha,
completed_replicates=len(group_records),
available_replicates=len(available),
rejected_replicates=rejected,
rejection_rate=rate,
monte_carlo_se=se,
unavailable_replicates=len(group_records) - len(available),
)
)
return summaries
|
Phase 4 diagnostics and study gate
Phase 4 adds two optional, additive layers on top of the existing harness: a
bounded orientation-attribution diagnostic computed from the same fitted PLS
estimator that produced the trajectory scores, and a predeclared gate that turns
the study's operating characteristics into an explicit proceed / hold /
indeterminate decision for the paper-grade study.
Both are off by default. SimulationEvaluationParams.attribution is a disabled
AttributionDiagnosticSettings, and a study configuration without an
attribution or matched_seeds block behaves exactly as before. Enabling
attribution requires integration_method="pls"; concat and snf produce no
fitted estimator to condition on and are rejected with a descriptive error.
Every gate threshold lives in the study configuration's acceptance.gate block,
never in the summary or report code. See
examples/trajectory_power_study/README.md for the config fields, the matched-
seed and shared-anchor semantics, the gate rules, and the full output inventory.
motco.simulations.AttributionDiagnosticSettings
dataclass
Optional orientation-attribution diagnostic settings for one evaluation.
Attribution is disabled by default: the zero-argument instance carries
enabled=False and no bootstrap work, so evaluations that do not ask for
diagnostics behave exactly as before. When enabled the settings are frozen
into the cell signature, so a diagnostic change cannot silently resume into
an existing shard.
Source code in src/motco/simulations/evaluation.py
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104 | @dataclass(frozen=True)
class AttributionDiagnosticSettings:
"""Optional orientation-attribution diagnostic settings for one evaluation.
Attribution is **disabled by default**: the zero-argument instance carries
``enabled=False`` and no bootstrap work, so evaluations that do not ask for
diagnostics behave exactly as before. When enabled the settings are frozen
into the cell signature, so a diagnostic change cannot silently resume into
an existing shard.
"""
enabled: bool = False
bootstrap_replicates: int = 0
bootstrap_seed: int | None = 0
top_k: int = 20
zero_tolerance: float = 1e-12
|
motco.simulations.attribution_diagnostics.compute_attribution_diagnostics(dataset=None, *, model, features, original_scales, settings, group_col, stage_col, selected_components, feature_order_signature, methylation_units='mvalue', layers=OMIC_LAYERS, metadata=None, groups=None, stages=None, truth=None, original_basis=None)
Run frozen-model attribution and reduce it to a compact JSON-safe record.
See :func:run_attribution_diagnostics for the arguments; this returns only
the compact record.
Source code in src/motco/simulations/attribution_diagnostics.py
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215 | def compute_attribution_diagnostics(
dataset: SemiSyntheticTrajectoryDataset | None = None,
*,
model: Any,
features: pd.DataFrame,
original_scales: np.ndarray | None,
settings: Any,
group_col: str,
stage_col: str,
selected_components: int | None,
feature_order_signature: str | None,
methylation_units: str = "mvalue",
layers: Sequence[OmicsLayer] = OMIC_LAYERS,
metadata: pd.DataFrame | None = None,
groups: Sequence[str] | None = None,
stages: Sequence[str] | None = None,
truth: TruthDrivers | None = None,
original_basis: Mapping[str, str] | None = None,
) -> dict[str, Any]:
"""Run frozen-model attribution and reduce it to a compact JSON-safe record.
See :func:`run_attribution_diagnostics` for the arguments; this returns only
the compact record.
"""
record, _ = run_attribution_diagnostics(
dataset,
model=model,
features=features,
original_scales=original_scales,
settings=settings,
group_col=group_col,
stage_col=stage_col,
selected_components=selected_components,
feature_order_signature=feature_order_signature,
methylation_units=methylation_units,
layers=layers,
metadata=metadata,
groups=groups,
stages=stages,
truth=truth,
original_basis=original_basis,
)
return record
|
motco.simulations.attribution_diagnostics.derive_truth_driver_features(dataset, feature_names, stage_levels, layers=OMIC_LAYERS)
Aligned truth drivers from changed group-stage differential patterns.
A feature is a truth driver for the transition s -> s+1 when the two
groups' population mean changes over that transition differ, i.e. when
delta_A * (v_A[s+1] - v_A[s]) != delta_B * (v_B[s+1] - v_B[s])
for that feature's per-omic effect size delta and per-stage differential
indicator v. Because expression and protein indicators are re-derived
from methylation through the CpG→gene→protein incidence maps, propagated
cascade drivers satisfy this definition on their own layers — a real
downstream driver is therefore never counted as a false positive. The
definition also covers effect-size-only constructions such as
magnitude_kind='all', where the indicators are identical between groups
but the per-omic delta is not.
The global driver set is the union over transitions. layers names the
measured blocks, in canonical order, that make up the joint feature order.
Source code in src/motco/simulations/attribution_diagnostics.py
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169 | def derive_truth_driver_features(
dataset: SemiSyntheticTrajectoryDataset,
feature_names: Sequence[str],
stage_levels: Sequence[str],
layers: Sequence[OmicsLayer] = OMIC_LAYERS,
) -> TruthDrivers:
"""Aligned truth drivers from changed group-stage differential patterns.
A feature is a truth driver for the transition ``s -> s+1`` when the two
groups' *population* mean changes over that transition differ, i.e. when
``delta_A * (v_A[s+1] - v_A[s]) != delta_B * (v_B[s+1] - v_B[s])``
for that feature's per-omic effect size ``delta`` and per-stage differential
indicator ``v``. Because expression and protein indicators are re-derived
from methylation through the CpG→gene→protein incidence maps, propagated
cascade drivers satisfy this definition on their own layers — a real
downstream driver is therefore never counted as a false positive. The
definition also covers effect-size-only constructions such as
``magnitude_kind='all'``, where the indicators are identical between groups
but the per-omic delta is not.
The global driver set is the union over transitions. ``layers`` names the
measured blocks, in canonical order, that make up the joint feature order.
"""
truth = dataset.truth
indicators = truth.get("indicators")
group_labels = truth.get("group_labels")
deltas = truth.get("deltas")
names = tuple(str(name) for name in feature_names)
n_transitions = max(len(stage_levels) - 1, 0)
transition_ids = [
f"{stage_levels[index]}->{stage_levels[index + 1]}" for index in range(n_transitions)
]
if not indicators or not group_labels or not deltas:
empty = np.zeros(len(names), dtype=bool)
return TruthDrivers(
feature_names=names,
global_mask=empty,
per_transition={transition_id: empty.copy() for transition_id in transition_ids},
available=False,
definition="unavailable: dataset truth carries no differential indicators",
)
label_a, label_b = (str(label) for label in group_labels)
per_transition: dict[str, np.ndarray] = {}
blocks: list[list[np.ndarray]] = [[] for _ in transition_ids]
for layer in layers:
layer_index = OMIC_LAYERS.index(layer)
ind_a = np.asarray(indicators[label_a][layer], dtype=float)
ind_b = np.asarray(indicators[label_b][layer], dtype=float)
delta_a = float(deltas[label_a][layer_index])
delta_b = float(deltas[label_b][layer_index])
for index in range(n_transitions):
step_a = delta_a * (ind_a[:, index + 1] - ind_a[:, index])
step_b = delta_b * (ind_b[:, index + 1] - ind_b[:, index])
blocks[index].append(np.abs(step_a - step_b) > _TRUTH_TOLERANCE)
for index, transition_id in enumerate(transition_ids):
mask = np.concatenate(blocks[index]) if blocks[index] else np.zeros(0, dtype=bool)
if mask.shape[0] != len(names):
raise AttributionDiagnosticError(
f"Truth indicators ({mask.shape[0]}) do not align to the joint feature order ({len(names)})."
)
per_transition[transition_id] = mask
global_mask = (
np.logical_or.reduce(list(per_transition.values()))
if per_transition
else np.zeros(len(names), dtype=bool)
)
return TruthDrivers(
feature_names=names,
global_mask=global_mask,
per_transition=per_transition,
available=True,
definition=(
"features whose group-stage differential mean change differs between groups, "
"including CpG->gene->protein propagated effects"
),
)
|
motco.simulations.study.AttributionSelector
dataclass
Which enumerated cells receive orientation-attribution diagnostics.
The selector is resolved into per-cell evaluation parameters during
enumeration, so eligibility and the effective bootstrap settings are part of
the cell's parameter signature rather than something inferred after a run.
Selection is never conditioned on an observed p-value.
Source code in src/motco/simulations/study/config.py
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156 | @dataclass(frozen=True)
class AttributionSelector:
"""Which enumerated cells receive orientation-attribution diagnostics.
The selector is resolved into per-cell evaluation parameters during
enumeration, so eligibility and the effective bootstrap settings are part of
the cell's parameter signature rather than something inferred after a run.
Selection is never conditioned on an observed p-value.
"""
enabled: bool = False
trajectory_modes: tuple[str, ...] = ("orientation",)
phases: tuple[str, ...] = ("power_primary",)
nonzero_effects_only: bool = True
effect_sizes: tuple[float, ...] | None = None
bootstrap_replicates: int = 100
bootstrap_seed: int | None = 0
top_k: int = 20
zero_tolerance: float = 1e-12
def __post_init__(self) -> None:
if not self.enabled:
return
unknown = sorted(set(self.trajectory_modes) - _TRAJECTORY_MODES)
if unknown:
raise StudyConfigError(f"attribution.trajectory_modes has unknown mode(s): {unknown}.")
if not self.trajectory_modes:
raise StudyConfigError("attribution.trajectory_modes must contain at least one mode.")
if not self.phases:
raise StudyConfigError("attribution.phases must contain at least one phase.")
if self.effect_sizes is not None:
if not self.effect_sizes:
raise StudyConfigError("attribution.effect_sizes must be null or a non-empty sequence.")
for value in self.effect_sizes:
if float(value) < 0:
raise StudyConfigError(f"attribution.effect_sizes must be non-negative; got {value}.")
if self.bootstrap_replicates < 0:
raise StudyConfigError("attribution.bootstrap_replicates must be non-negative.")
if self.top_k < 1:
raise StudyConfigError("attribution.top_k must be at least 1.")
if not (self.zero_tolerance > 0):
raise StudyConfigError("attribution.zero_tolerance must be positive.")
def settings(self) -> AttributionDiagnosticSettings:
"""Effective evaluation-level settings implied by this selector."""
return AttributionDiagnosticSettings(
enabled=True,
bootstrap_replicates=int(self.bootstrap_replicates),
bootstrap_seed=None if self.bootstrap_seed is None else int(self.bootstrap_seed),
top_k=int(self.top_k),
zero_tolerance=float(self.zero_tolerance),
)
def selects(self, *, phase: str, trajectory_mode: str | None, effect_size: float | None) -> bool:
"""Whether a cell with this identity is eligible for attribution."""
if not self.enabled:
return False
if phase not in self.phases:
return False
if trajectory_mode is None or trajectory_mode not in self.trajectory_modes:
return False
if effect_size is None:
return False
value = float(effect_size)
if self.nonzero_effects_only and value == 0.0:
return False
if self.effect_sizes is not None:
return any(abs(value - float(candidate)) < 1e-12 for candidate in self.effect_sizes)
return True
|
settings()
Effective evaluation-level settings implied by this selector.
Source code in src/motco/simulations/study/config.py
129
130
131
132
133
134
135
136
137
138 | def settings(self) -> AttributionDiagnosticSettings:
"""Effective evaluation-level settings implied by this selector."""
return AttributionDiagnosticSettings(
enabled=True,
bootstrap_replicates=int(self.bootstrap_replicates),
bootstrap_seed=None if self.bootstrap_seed is None else int(self.bootstrap_seed),
top_k=int(self.top_k),
zero_tolerance=float(self.zero_tolerance),
)
|
selects(*, phase, trajectory_mode, effect_size)
Whether a cell with this identity is eligible for attribution.
Source code in src/motco/simulations/study/config.py
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156 | def selects(self, *, phase: str, trajectory_mode: str | None, effect_size: float | None) -> bool:
"""Whether a cell with this identity is eligible for attribution."""
if not self.enabled:
return False
if phase not in self.phases:
return False
if trajectory_mode is None or trajectory_mode not in self.trajectory_modes:
return False
if effect_size is None:
return False
value = float(effect_size)
if self.nonzero_effects_only and value == 0.0:
return False
if self.effect_sizes is not None:
return any(abs(value - float(candidate)) < 1e-12 for candidate in self.effect_sizes)
return True
|
motco.simulations.study.MatchedSeedPolicy
dataclass
Opt-in matched generator seeds shared across primary cells.
When enabled, every primary power cell draws its generator seed from a
shared (seed family, replicate index) key instead of its own cell id, so
requested-effect comparisons are paired at the generated-reference level.
Negative-control and OFAT cells keep their own seed families so they stay
independent. The policy version enters the parameter signature.
Source code in src/motco/simulations/study/config.py
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181 | @dataclass(frozen=True)
class MatchedSeedPolicy:
"""Opt-in matched generator seeds shared across primary cells.
When enabled, every primary power cell draws its generator seed from a
shared ``(seed family, replicate index)`` key instead of its own cell id, so
requested-effect comparisons are paired at the generated-reference level.
Negative-control and OFAT cells keep their own seed families so they stay
independent. The policy version enters the parameter signature.
"""
enabled: bool = False
version: int = 1
primary_family: str = "primary"
shared_zero_effect_anchor: bool = True
def __post_init__(self) -> None:
if not self.enabled:
return
if self.version < 1:
raise StudyConfigError("matched_seeds.version must be at least 1.")
if not str(self.primary_family).strip():
raise StudyConfigError("matched_seeds.primary_family must be a non-empty string.")
|
motco.simulations.study.Phase4GateConfig
dataclass
Every parameter the Phase 4 gate consumes.
Summary and report code reads all thresholds from here and hard-codes none
of its own, so a gate can be re-specified by editing the committed config.
Source code in src/motco/simulations/study/config.py
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257 | @dataclass(frozen=True)
class Phase4GateConfig:
"""Every parameter the Phase 4 gate consumes.
Summary and report code reads all thresholds from here and hard-codes none
of its own, so a gate can be re-specified by editing the committed config.
"""
enabled: bool = False
alpha: float = 0.05
control_se_tolerance: float = 2.0
monotonicity_se_tolerance: float = 2.0
min_power_at_top: float = 0.80
confirmation_se_threshold: float = 1.0
max_marginal_exceedances: int = 1
control_modes: tuple[str, ...] = ("none", "translation")
rules: tuple[GateRule, ...] = ()
require_complete_records: bool = True
def __post_init__(self) -> None:
if not self.enabled:
return
if not (0 < self.alpha < 1):
raise StudyConfigError("acceptance.gate.alpha must be between 0 and 1.")
for name in ("control_se_tolerance", "monotonicity_se_tolerance", "confirmation_se_threshold"):
value = float(getattr(self, name))
if value < 0:
raise StudyConfigError(f"acceptance.gate.{name} must be non-negative.")
if not (0 <= self.min_power_at_top <= 1):
raise StudyConfigError("acceptance.gate.min_power_at_top must be between 0 and 1.")
if self.max_marginal_exceedances < 0:
raise StudyConfigError("acceptance.gate.max_marginal_exceedances must be non-negative.")
unknown = sorted(set(self.control_modes) - _TRAJECTORY_MODES)
if unknown:
raise StudyConfigError(f"acceptance.gate.control_modes has unknown mode(s): {unknown}.")
if not self.control_modes:
raise StudyConfigError("acceptance.gate.control_modes must contain at least one mode.")
if not self.rules:
raise StudyConfigError("acceptance.gate.rules must declare every mandatory and descriptive pair.")
if not any(rule.role == "mandatory_power" for rule in self.rules):
raise StudyConfigError("acceptance.gate.rules must declare at least one 'mandatory_power' rule.")
seen: set[tuple[str, str]] = set()
for rule in self.rules:
key = (rule.trajectory_mode, rule.statistic)
if key in seen:
raise StudyConfigError(f"acceptance.gate.rules has a duplicate entry for {key}.")
seen.add(key)
def power_floor(self, rule: GateRule) -> float:
"""Minimum power at the top effect for ``rule``."""
return float(self.min_power_at_top if rule.min_power_at_top is None else rule.min_power_at_top)
|
power_floor(rule)
Minimum power at the top effect for rule.
Source code in src/motco/simulations/study/config.py
| def power_floor(self, rule: GateRule) -> float:
"""Minimum power at the top effect for ``rule``."""
return float(self.min_power_at_top if rule.min_power_at_top is None else rule.min_power_at_top)
|
motco.simulations.study.GateRule
dataclass
One predeclared mode/statistic gate role.
Source code in src/motco/simulations/study/config.py
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203 | @dataclass(frozen=True)
class GateRule:
"""One predeclared mode/statistic gate role."""
trajectory_mode: str
statistic: str
role: str
min_power_at_top: float | None = None
def __post_init__(self) -> None:
if self.trajectory_mode not in _TRAJECTORY_MODES:
raise StudyConfigError(f"gate rule trajectory_mode {self.trajectory_mode!r} is unknown.")
if self.statistic not in _STATISTICS:
raise StudyConfigError(f"gate rule statistic {self.statistic!r} is unknown.")
if self.role not in _GATE_ROLES:
raise StudyConfigError(f"gate rule role {self.role!r} is unknown; expected one of {sorted(_GATE_ROLES)}.")
@property
def name(self) -> str:
return f"{self.role}[{self.trajectory_mode},{self.statistic}]"
|
motco.simulations.study.enumerate_study(config)
Build the combined Type I + power grid for a study.
Source code in src/motco/simulations/study/enumerate.py
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87 | def enumerate_study(config: StudyConfig) -> SimulationGrid:
"""Build the combined Type I + power grid for a study."""
resolver = _attribution_resolver(config)
type_i = enumerate_type_i_grid(
baseline_generator_params=config.generator,
evaluation_params=config.evaluation,
axes=config.axes,
n_replicates=config.n_replicates,
base_seed=config.base_seed,
)
cells = list(type_i.cells)
cells.extend(_negative_control_cells(config))
cells.extend(_power_cells(config, resolver))
cells.extend(_design_point_cells(config, resolver))
if config.matched_seeds.enabled:
cells = _assign_seed_families(cells, config)
_require_unique_ids(cells)
_require_surgery_headroom(cells)
grid_metadata: dict[str, Any] = {
"grid_type": "study",
"study_config_metadata": dict(config.metadata),
"matched_seed_policy": {
"enabled": config.matched_seeds.enabled,
"version": config.matched_seeds.version,
"primary_family": config.matched_seeds.primary_family,
"shared_zero_effect_anchor": config.matched_seeds.shared_zero_effect_anchor,
},
}
if config.design_grid.enabled:
grid_metadata["design_grid"] = {
"axes": {axis: list(values) for axis, values in config.design_grid.axes.items()},
"n_points": len(config.design_points()),
"baseline_point": _to_jsonable(config.baseline_design_point()),
}
if config.design_grid.explicit:
grid_metadata["design_grid"]["points"] = _to_jsonable(config.design_points())
return SimulationGrid(cells=tuple(cells), metadata=grid_metadata)
|
motco.simulations.study.summarize_realized_geometry(records)
Long-form geometry by mode, effect, checkpoint, scope, and statistic.
Path lengths are summarized alongside delta / angle / shape.
Unavailable values (a null path with no direction, or shape with fewer
than three stages) are counted rather than coerced to zero, and each row
carries the measurement space its numbers live in so feature-space and
latent-space magnitudes are never silently compared.
Source code in src/motco/simulations/study/phase4.py
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256 | def summarize_realized_geometry(records: Sequence[SimulationReplicateResult]) -> pd.DataFrame:
"""Long-form geometry by mode, effect, checkpoint, scope, and statistic.
Path lengths are summarized alongside ``delta`` / ``angle`` / ``shape``.
Unavailable values (a null path with no direction, or ``shape`` with fewer
than three stages) are counted rather than coerced to zero, and each row
carries the measurement space its numbers live in so feature-space and
latent-space magnitudes are never silently compared.
"""
buckets: dict[tuple, list[float | None]] = {}
identity: dict[tuple, dict[str, Any]] = {}
for record in records:
if record.status != "completed":
continue
geometry = record.realized_geometry or {}
meta = dict(record.cell_metadata)
mode = str(meta.get("trajectory_mode") or geometry.get("requested", {}).get("trajectory_mode") or "none")
effect_size = meta.get("effect_size")
effect_size = float(effect_size) if effect_size is not None else None
for checkpoint, scopes in (geometry.get("checkpoints") or {}).items():
for scope, values in scopes.items():
for statistic, value in _geometry_statistics(values):
key = (mode, effect_size, record.cell_id, checkpoint, scope, statistic)
buckets.setdefault(key, []).append(value)
identity.setdefault(
key,
{
"trajectory_mode": mode,
"effect_size": effect_size,
"cell_id": record.cell_id,
"phase": record.phase,
"checkpoint": checkpoint,
"scope": scope,
"measurement_space": MEASUREMENT_SPACES.get(checkpoint, "unknown"),
"statistic": statistic,
},
)
rows: list[dict[str, Any]] = []
for key, values in buckets.items():
available = [value for value in values if value is not None and math.isfinite(value)]
row = dict(identity[key])
row.update(
{
"n_replicates": len(values),
"n_available": len(available),
"n_unavailable": len(values) - len(available),
"mean": _mean(available),
"sd": _sd(available),
"median": _median(available),
}
)
rows.append(row)
frame = pd.DataFrame(
rows,
columns=[
"trajectory_mode",
"effect_size",
"cell_id",
"phase",
"checkpoint",
"scope",
"measurement_space",
"statistic",
"n_replicates",
"n_available",
"n_unavailable",
"mean",
"sd",
"median",
],
)
if frame.empty:
return frame
frame["checkpoint_order"] = frame["checkpoint"].map(
{name: index for index, name in enumerate(CHECKPOINT_ORDER)}
)
return frame.sort_values(
["trajectory_mode", "effect_size", "checkpoint_order", "scope", "statistic"]
).drop(columns=["checkpoint_order"]).reset_index(drop=True)
|
Angle-pivotality diagnostic
motco.simulations.pivotality.association_table(records, *, statistics=STATISTICS, null_targets=NULL_TARGETS)
Per cell and statistic, how strongly the null tracks the observed value.
A pivotal statistic shows no association: the null mean, spread, and
critical value are properties of the design, not of the replicate's realized
signal. A correlation near +1 with a slope near 1 on q95 is the
non-pivotal case that costs power.
Source code in src/motco/simulations/pivotality.py
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261 | def association_table(
records: Iterable[SimulationReplicateResult],
*,
statistics: Sequence[str] = STATISTICS,
null_targets: Sequence[str] = NULL_TARGETS,
) -> list[AssociationRow]:
"""Per cell and statistic, how strongly the null tracks the observed value.
A pivotal statistic shows no association: the null mean, spread, and
critical value are properties of the design, not of the replicate's realized
signal. A correlation near +1 with a slope near 1 on ``q95`` is the
non-pivotal case that costs power.
"""
rows: list[AssociationRow] = []
for cell, group in group_by_cell(records).items():
for statistic in statistics:
for target in null_targets:
x, y = _paired(group, statistic, target)
r, low, high = _correlation_with_interval(x, y)
rows.append(
AssociationRow(
cell=cell,
statistic=statistic,
null_target=target,
n_replicates=int(x.size),
correlation=r,
correlation_ci_low=low,
correlation_ci_high=high,
slope=_slope(x, y),
mean_observed=float(x.mean()) if x.size else None,
mean_null_target=float(y.mean()) if y.size else None,
)
)
return rows
|
motco.simulations.pivotality.rejection_split_table(records, *, alpha=0.05, statistics=STATISTICS)
Split the observed statistic — and its critical value — by rejection outcome.
inverted flags the Phase 4 pattern: non-rejecting replicates carrying a
larger mean observed statistic than rejecting ones. Reporting the mean
critical value beside it answers the follow-up the spec requires — whether
replicates with a larger observed statistic also carry a proportionally
larger bar.
Source code in src/motco/simulations/pivotality.py
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439 | def rejection_split_table(
records: Iterable[SimulationReplicateResult],
*,
alpha: float = 0.05,
statistics: Sequence[str] = STATISTICS,
) -> list[RejectionSplitRow]:
"""Split the observed statistic — and its critical value — by rejection outcome.
``inverted`` flags the Phase 4 pattern: non-rejecting replicates carrying a
*larger* mean observed statistic than rejecting ones. Reporting the mean
critical value beside it answers the follow-up the spec requires — whether
replicates with a larger observed statistic also carry a proportionally
larger bar.
"""
if not (0 < alpha < 1):
raise PivotalityError("alpha must be between 0 and 1.")
rows: list[RejectionSplitRow] = []
for cell, group in group_by_cell(records).items():
for statistic in statistics:
rejecting: list[tuple[float, float | None]] = []
non_rejecting: list[tuple[float, float | None]] = []
for record in group:
p_value = record.p_values.get(statistic)
observed = _observed(record, statistic)
if p_value is None or observed is None or not math.isfinite(float(p_value)):
continue
critical = _null_feature(record, statistic, "q95")
target = rejecting if float(p_value) < alpha else non_rejecting
target.append((observed, critical))
mean_obs_rej = _mean_first(rejecting)
mean_obs_non = _mean_first(non_rejecting)
rows.append(
RejectionSplitRow(
cell=cell,
statistic=statistic,
alpha=alpha,
n_rejecting=len(rejecting),
n_non_rejecting=len(non_rejecting),
mean_observed_rejecting=mean_obs_rej,
mean_observed_non_rejecting=mean_obs_non,
mean_critical_rejecting=_mean_second(rejecting),
mean_critical_non_rejecting=_mean_second(non_rejecting),
inverted=(
None
if mean_obs_rej is None or mean_obs_non is None
else mean_obs_non > mean_obs_rej
),
)
)
return rows
|
motco.simulations.pivotality.standardized_counterfactual_table(records, *, alpha=0.05, statistics=STATISTICS, null_control_modes=('none', 'translation'))
Rejection rate under a z calibrated on the null-control cells.
Per replicate, z = (observed − null_mean) / null_sd. The threshold is the
1 − alpha quantile of z pooled over the null-control replicates for
that statistic, and every cell — controls included — is then scored against
it. Reporting the controls' own rates is what makes a power gain bought with
an inflated Type I rate visible.
This borrows a reference distribution that does not exist outside a
simulation. It is a diagnostic, never a deployable test; see the module
docstring.
Source code in src/motco/simulations/pivotality.py
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528 | def standardized_counterfactual_table(
records: Iterable[SimulationReplicateResult],
*,
alpha: float = 0.05,
statistics: Sequence[str] = STATISTICS,
null_control_modes: Sequence[str] = ("none", "translation"),
) -> list[StandardizedRow]:
"""Rejection rate under a ``z`` calibrated on the null-control cells.
Per replicate, ``z = (observed − null_mean) / null_sd``. The threshold is the
``1 − alpha`` quantile of ``z`` pooled over the *null-control* replicates for
that statistic, and every cell — controls included — is then scored against
it. Reporting the controls' own rates is what makes a power gain bought with
an inflated Type I rate visible.
This borrows a reference distribution that does not exist outside a
simulation. It is a diagnostic, never a deployable test; see the module
docstring.
"""
if not (0 < alpha < 1):
raise PivotalityError("alpha must be between 0 and 1.")
groups = group_by_cell(records)
control_modes = set(null_control_modes)
thresholds: dict[str, float | None] = {}
for statistic in statistics:
pooled = [
z
for cell, group in groups.items()
if cell.trajectory_mode in control_modes
for record in group
if (z := replicate_z(record, statistic)) is not None
]
thresholds[statistic] = float(np.quantile(pooled, 1.0 - alpha)) if pooled else None
rows: list[StandardizedRow] = []
for cell, group in groups.items():
for statistic in statistics:
zs = [z for record in group if (z := replicate_z(record, statistic)) is not None]
p_values = [
float(p)
for record in group
if (p := record.p_values.get(statistic)) is not None and math.isfinite(float(p))
]
threshold = thresholds[statistic]
rows.append(
StandardizedRow(
cell=cell,
statistic=statistic,
alpha=alpha,
n_replicates=len(zs),
as_specified_rate=(
float(np.mean([p < alpha for p in p_values])) if p_values else None
),
standardized_rate=(
float(np.mean([z >= threshold for z in zs]))
if zs and threshold is not None
else None
),
z_threshold=threshold,
mean_z=float(np.mean(zs)) if zs else None,
is_null_control=cell.trajectory_mode in control_modes,
)
)
return rows
|
motco.simulations.pivotality.replicate_z(record, statistic)
(observed − null_mean) / null_sd for one replicate, or None.
Source code in src/motco/simulations/pivotality.py
451
452
453
454
455
456
457
458
459
460 | def replicate_z(record: SimulationReplicateResult, statistic: str) -> float | None:
"""``(observed − null_mean) / null_sd`` for one replicate, or ``None``."""
observed = _observed(record, statistic)
mean = _null_feature(record, statistic, "mean")
sd = _null_feature(record, statistic, "sd")
if observed is None or mean is None or sd is None or sd <= 0.0:
return None
z = (observed - mean) / sd
return z if math.isfinite(z) else None
|
motco.simulations.pivotality.write_pivotality_tables(records, out_dir, *, alpha=0.05, statistics=STATISTICS)
Write the association, rejection-split, and counterfactual tables as CSV.
Source code in src/motco/simulations/pivotality.py
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590 | def write_pivotality_tables(
records: Sequence[SimulationReplicateResult],
out_dir: Path,
*,
alpha: float = 0.05,
statistics: Sequence[str] = STATISTICS,
) -> dict[str, Path]:
"""Write the association, rejection-split, and counterfactual tables as CSV."""
out_dir = Path(out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
tables: Mapping[str, list[dict[str, object]]] = {
"pivotality_association": rows_to_records(
association_table(records, statistics=statistics)
),
"pivotality_spectrum": rows_to_records(
spectrum_association_table(records, statistics=statistics)
),
"pivotality_rejection_split": rows_to_records(
rejection_split_table(records, alpha=alpha, statistics=statistics)
),
"pivotality_standardized": rows_to_records(
standardized_counterfactual_table(records, alpha=alpha, statistics=statistics)
),
}
paths: dict[str, Path] = {}
for name, table in tables.items():
path = out_dir / f"{name}.csv"
if not table:
path.write_text("", encoding="utf-8")
paths[name] = path
continue
with path.open("w", encoding="utf-8", newline="") as handle:
writer = csv.DictWriter(handle, fieldnames=list(table[0]))
writer.writeheader()
writer.writerows(table)
paths[name] = path
return paths
|
motco.simulations.study.summarize_pls_selection(records)
Selected component counts, effective CV settings, and AUROC per cell.
Source code in src/motco/simulations/study/phase4.py
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343 | def summarize_pls_selection(records: Sequence[SimulationReplicateResult]) -> pd.DataFrame:
"""Selected component counts, effective CV settings, and AUROC per cell."""
buckets: dict[str, dict[str, Any]] = {}
for record in records:
if record.status != "completed":
continue
meta = dict(record.cell_metadata)
entry = buckets.setdefault(
record.cell_id,
{
"cell_id": record.cell_id,
"phase": record.phase,
"trajectory_mode": str(meta.get("trajectory_mode") or "none"),
"effect_size": meta.get("effect_size"),
"completed_replicates": 0,
"missing_integration_metadata": 0,
"selected_lv": [],
"cv_mean_auroc": [],
"cv_settings": set(),
},
)
entry["completed_replicates"] += 1
integration = record.integration_metadata or {}
if not integration or integration.get("integration_method") != "pls":
entry["missing_integration_metadata"] += 1
continue
selected = integration.get("selected_lv")
if selected is None:
selected = (integration.get("integration_params") or {}).get("selected_lv")
if selected is None:
entry["missing_integration_metadata"] += 1
else:
entry["selected_lv"].append(int(selected))
auroc = integration.get("cv_mean_auroc")
if auroc is not None and math.isfinite(float(auroc)):
entry["cv_mean_auroc"].append(float(auroc))
params = integration.get("integration_params") or {}
entry["cv_settings"].add(
tuple(
(key, params.get(key))
for key in ("cv1_splits", "cv2_splits", "n_repeats", "max_components", "random_state")
)
)
rows: list[dict[str, Any]] = []
for entry in buckets.values():
selected = entry.pop("selected_lv")
auroc = entry.pop("cv_mean_auroc")
settings = entry.pop("cv_settings")
row = dict(entry)
row.update(
{
"n_selected_lv": len(selected),
"selected_lv_mean": _mean(selected),
"selected_lv_median": _median(selected),
"selected_lv_mode": _mode(selected),
"selected_lv_min": min(selected) if selected else None,
"selected_lv_max": max(selected) if selected else None,
"cv_mean_auroc_mean": _mean(auroc),
"cv_settings": _format_settings(settings),
"cv_settings_consistent": len(settings) <= 1,
}
)
rows.append(row)
frame = pd.DataFrame(rows)
if frame.empty:
return frame
return frame.sort_values(["trajectory_mode", "effect_size", "cell_id"]).reset_index(drop=True)
|
motco.simulations.study.summarize_attribution(records)
Attribution availability, retention, stability, and truth recovery.
Cross-replicate agreement is computed from the persisted top-k identifiers:
top_k_jaccard is the mean pairwise overlap of the selected sets across
replicates, and sign_agreement is the mean, over features selected in at
least two replicates, of the fraction agreeing with that feature's modal
sign. Both are None when fewer than two replicates carry the transition.
Source code in src/motco/simulations/study/phase4.py
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465 | def summarize_attribution(records: Sequence[SimulationReplicateResult]) -> pd.DataFrame:
"""Attribution availability, retention, stability, and truth recovery.
Cross-replicate agreement is computed from the persisted top-k identifiers:
``top_k_jaccard`` is the mean pairwise overlap of the selected sets across
replicates, and ``sign_agreement`` is the mean, over features selected in at
least two replicates, of the fraction agreeing with that feature's modal
sign. Both are ``None`` when fewer than two replicates carry the transition.
"""
groups: dict[tuple, dict[str, Any]] = {}
for record in records:
if record.status != "completed":
continue
meta = dict(record.cell_metadata)
mode = str(meta.get("trajectory_mode") or "none")
effect_size = meta.get("effect_size")
effect_size = float(effect_size) if effect_size is not None else None
diagnostics = record.attribution_diagnostics or {}
status = record.attribution_status or diagnostics.get("status") or "not_requested"
if status == "not_requested":
continue
per_transition = {
(entry["transition_id"], component)
for entry in diagnostics.get("transitions", [])
for component in _COMPONENTS
}
if not per_transition:
key: tuple[Any, ...] = (mode, effect_size, record.cell_id, None, None)
entry_group = _attribution_group(groups, key, mode, effect_size, record)
entry_group["eligible"] += 1
entry_group["failed"] += 1 if status == "failed" else 0
continue
retention = {entry["transition_id"]: entry.get("retention", {}) for entry in diagnostics.get("transitions", [])}
recovery = {
(entry["transition_id"], entry["component"]): entry
for entry in diagnostics.get("truth_recovery", [])
}
stability = {
(entry["transition_id"], entry["component"]): entry for entry in diagnostics.get("stability", [])
}
selected: dict[tuple[str, str], list[tuple[str, int]]] = {}
for entry in diagnostics.get("top_features", []):
selected.setdefault((entry["transition_id"], entry["component"]), []).append(
(str(entry["feature"]), int(entry.get("sign", 0)))
)
for transition_id, component in sorted(per_transition):
key = (mode, effect_size, record.cell_id, transition_id, component)
entry_group = _attribution_group(groups, key, mode, effect_size, record)
entry_group["eligible"] += 1
entry_group["computed"] += 1
transition_retention = retention.get(transition_id, {})
if component == "observed" and transition_retention.get("available"):
_append(entry_group["retention_cosine"], transition_retention.get("cosine"))
_append(entry_group["retention_norm_ratio"], transition_retention.get("norm_ratio"))
_append(entry_group["residual_fraction"], transition_retention.get("residual_fraction"))
recovery_entry = recovery.get((transition_id, component))
if recovery_entry and recovery_entry.get("available"):
_append(entry_group["precision"], recovery_entry.get("precision"))
_append(entry_group["recall"], recovery_entry.get("recall"))
_append(entry_group["selected_count"], recovery_entry.get("selected_count"))
stability_entry = stability.get((transition_id, component))
if stability_entry:
_append(entry_group["sign_stability"], stability_entry.get("mean_sign_stability"))
_append(
entry_group["top_k_selection_frequency"],
stability_entry.get("mean_top_k_selection_frequency"),
)
features = selected.get((transition_id, component))
if features:
entry_group["selections"].append(features)
rows: list[dict[str, Any]] = []
for entry in groups.values():
selections = entry.pop("selections")
row = {
"trajectory_mode": entry["trajectory_mode"],
"effect_size": entry["effect_size"],
"cell_id": entry["cell_id"],
"transition_id": entry["transition_id"],
"component": entry["component"],
"eligible_replicates": entry["eligible"],
"computed_replicates": entry["computed"],
"failed_replicates": entry["failed"],
"availability_rate": (entry["computed"] / entry["eligible"]) if entry["eligible"] else None,
"retention_cosine_mean": _mean(entry["retention_cosine"]),
"retention_norm_ratio_mean": _mean(entry["retention_norm_ratio"]),
"residual_fraction_mean": _mean(entry["residual_fraction"]),
"precision_mean": _mean(entry["precision"]),
"recall_mean": _mean(entry["recall"]),
"selected_count_mean": _mean(entry["selected_count"]),
"bootstrap_sign_stability_mean": _mean(entry["sign_stability"]),
"bootstrap_top_k_frequency_mean": _mean(entry["top_k_selection_frequency"]),
"top_k_jaccard": _mean_pairwise_jaccard(selections),
"sign_agreement": _sign_agreement(selections),
}
rows.append(row)
frame = pd.DataFrame(rows)
if frame.empty:
return frame
return frame.sort_values(
["trajectory_mode", "effect_size", "cell_id", "transition_id", "component"],
na_position="first",
).reset_index(drop=True)
|
motco.simulations.study.localize_off_diagonal(records, *, scope='joint', materiality_threshold=None, rule='absolute', dust_atol=MATERIALITY_DUST_ATOL)
First checkpoint at which each off-diagonal response becomes material.
Every checkpoint is judged on its own scale-free quantity, so a
standardized-feature distance is never weighed against a PLS-latent one:
delta is divided by that checkpoint's mean group path length, giving a
dimensionless relative size difference;angle is divided by 180, the estimator reporting orientation in
degrees on [0, 180], giving the fraction of the maximum possible
orientation difference; andshape is already a Procrustes distance on centroid-scaled configurations.
Two rule values decide when that excess counts as material:
"absolute" (the default, and the rule every already-reported run was
judged under): material when the excess is at least
materiality_threshold, defaulting to
:data:DEFAULT_MATERIALITY_THRESHOLD. One absolute cut cannot serve all
three statistics — shape responses run ~1e-2 while the cut is 5e-2, so
a shape response the RRPP test rejects at high rate still reports
not_material. Kept as the default only so committed reports keep
reproducing."null_dispersion": material when the excess is at least
materiality_threshold (defaulting to
:data:DEFAULT_MATERIALITY_DISPERSION_UNITS) times the anchor's dispersion
at the same checkpoint, which puts the three statistics in common units.
Where that dispersion is degenerate — zero or dust, as it is at every
population checkpoint for an analytically null construction — the rule
falls back to "excess exceeds dust_atol" rather than dividing by dust.
materiality_basis records which path decided each row.
Comparing an excess rather than a ratio matters throughout: an exactly-null
construction leaves the anchor's angle at roughly 1e-8, and any ratio against
that is meaningless.
The emitted classification labels where a response first appears —
construction, sampling/preprocessing, or projection. It is descriptive, not
causal, and gates nothing.
Source code in src/motco/simulations/study/phase4.py
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773 | def localize_off_diagonal(
records: Sequence[SimulationReplicateResult],
*,
scope: str = "joint",
materiality_threshold: float | None = None,
rule: str = "absolute",
dust_atol: float = MATERIALITY_DUST_ATOL,
) -> pd.DataFrame:
"""First checkpoint at which each off-diagonal response becomes material.
Every checkpoint is judged on its **own** scale-free quantity, so a
standardized-feature distance is never weighed against a PLS-latent one:
- ``delta`` is divided by that checkpoint's mean group path length, giving a
dimensionless relative size difference;
- ``angle`` is divided by 180, the estimator reporting orientation in
**degrees** on ``[0, 180]``, giving the fraction of the maximum possible
orientation difference; and
- ``shape`` is already a Procrustes distance on centroid-scaled configurations.
Two ``rule`` values decide when that excess counts as material:
- ``"absolute"`` (the default, and the rule every already-reported run was
judged under): material when the excess is at least
``materiality_threshold``, defaulting to
:data:`DEFAULT_MATERIALITY_THRESHOLD`. One absolute cut cannot serve all
three statistics — ``shape`` responses run ~1e-2 while the cut is 5e-2, so
a ``shape`` response the RRPP test rejects at high rate still reports
``not_material``. Kept as the default only so committed reports keep
reproducing.
- ``"null_dispersion"``: material when the excess is at least
``materiality_threshold`` (defaulting to
:data:`DEFAULT_MATERIALITY_DISPERSION_UNITS`) times the anchor's dispersion
at the same checkpoint, which puts the three statistics in common units.
Where that dispersion is degenerate — zero or dust, as it is at every
population checkpoint for an analytically null construction — the rule
falls back to "excess exceeds ``dust_atol``" rather than dividing by dust.
``materiality_basis`` records which path decided each row.
Comparing an excess rather than a ratio matters throughout: an exactly-null
construction leaves the anchor's angle at roughly 1e-8, and any ratio against
that is meaningless.
The emitted classification labels where a response first appears —
construction, sampling/preprocessing, or projection. It is descriptive, not
causal, and gates nothing.
"""
geometry = summarize_realized_geometry(records)
columns = [
"trajectory_mode",
"effect_size",
"statistic",
"first_material_checkpoint",
"classification",
"materiality_threshold",
"normalized_value",
"normalized_null",
"normalized_excess",
"measurement_space",
]
# The legacy rule keeps the legacy schema exactly. Adding columns under it
# would change every already-committed `phase4_localization.csv`, and a
# frozen report is asserted to regenerate byte-identically; the extra
# columns describe the new rule's reasoning, so they travel with it.
if rule != "absolute":
columns = [
*columns[:5],
"materiality_rule",
"materiality_threshold",
"materiality_basis",
"normalized_value",
"normalized_null",
"normalized_excess",
"null_dispersion",
"excess_in_dispersion_units",
"measurement_space",
]
if rule not in MATERIALITY_RULES:
raise Phase4SummaryError(
f"Unknown materiality rule {rule!r}; expected one of {MATERIALITY_RULES}."
)
if materiality_threshold is None:
materiality_threshold = (
DEFAULT_MATERIALITY_THRESHOLD
if rule == "absolute"
else DEFAULT_MATERIALITY_DISPERSION_UNITS
)
if geometry.empty:
return pd.DataFrame(columns=columns)
anchor_cells = {
record.cell_id
for record in records
if record.cell_metadata.get("zero_effect_anchor")
or (record.phase.startswith("type_i_") and not record.cell_metadata.get("trajectory_mode"))
}
scoped = geometry[geometry["scope"] == scope]
normalized = _normalized_geometry(scoped)
null_reference = {
(checkpoint, statistic): value
for (cell_id, checkpoint, statistic), value in normalized.items()
if cell_id in anchor_cells
}
dispersion = _normalized_geometry(scoped, column="sd")
null_dispersion = {
(checkpoint, statistic): value
for (cell_id, checkpoint, statistic), value in dispersion.items()
if cell_id in anchor_cells
}
rows: list[dict[str, Any]] = []
grouped: dict[tuple[str, Any, str], dict[str, float]] = {}
identity: dict[tuple[str, Any, str], dict[str, Any]] = {}
for row in scoped[scoped["statistic"].isin(_STATISTICS)].itertuples(index=False):
if row.cell_id in anchor_cells:
continue
key = (str(row.trajectory_mode), row.effect_size, str(row.statistic))
value = normalized.get((str(row.cell_id), str(row.checkpoint), str(row.statistic)))
if value is not None:
grouped.setdefault(key, {})[str(row.checkpoint)] = value
identity.setdefault(key, {"cell_id": str(row.cell_id)})
for key in sorted(grouped, key=lambda k: (k[0], k[2], k[1] if k[1] is not None else -1.0)):
mode, effect_size, statistic = key
if DIAGONAL_STATISTIC.get(mode) == statistic:
continue
by_checkpoint = grouped[key]
first_checkpoint: str | None = None
observed: float | None = None
null_value: float | None = None
excess: float | None = None
spread: float | None = None
excess_units: float | None = None
basis: str | None = None
for checkpoint in CHECKPOINT_ORDER:
value = by_checkpoint.get(checkpoint)
reference = null_reference.get((checkpoint, statistic))
if value is None:
continue
baseline = 0.0 if reference is None else reference
candidate_excess = value - baseline
candidate_spread = null_dispersion.get((checkpoint, statistic))
candidate_units: float | None = None
if rule == "absolute":
candidate_basis = "absolute_excess"
material = candidate_excess >= materiality_threshold
elif candidate_spread is not None and candidate_spread > dust_atol:
# Usable null dispersion: judge in its units.
candidate_basis = "null_dispersion"
candidate_units = candidate_excess / candidate_spread
material = candidate_units >= materiality_threshold
else:
# Degenerate null (analytically zero at population checkpoints):
# any excess above dust is material. Never divide by dust.
candidate_basis = "degenerate_null_dust_floor"
material = candidate_excess > dust_atol
if material:
first_checkpoint = checkpoint
observed = value
null_value = reference
excess = candidate_excess
spread = candidate_spread
excess_units = candidate_units
basis = candidate_basis
break
localization_row: dict[str, Any] = {
"trajectory_mode": mode,
"effect_size": effect_size,
"statistic": statistic,
"first_material_checkpoint": first_checkpoint,
"classification": (
CHECKPOINT_CLASSIFICATION.get(first_checkpoint, "unclassified")
if first_checkpoint is not None
else "not_material"
),
"materiality_threshold": materiality_threshold,
"normalized_value": observed,
"normalized_null": null_value,
"normalized_excess": excess,
"measurement_space": (
MEASUREMENT_SPACES.get(first_checkpoint) if first_checkpoint is not None else None
),
}
if rule != "absolute":
localization_row["materiality_rule"] = rule
localization_row["materiality_basis"] = basis
localization_row["null_dispersion"] = spread
localization_row["excess_in_dispersion_units"] = excess_units
rows.append(localization_row)
frame = pd.DataFrame(rows, columns=columns)
if frame.empty:
return frame
return frame.sort_values(["trajectory_mode", "statistic", "effect_size"]).reset_index(drop=True)
|
motco.simulations.study.type_i_inflation_bound(alpha, se_tolerance, n_available)
One-sided inflation bound alpha + k * sqrt(alpha (1 - alpha) / n).
Source code in src/motco/simulations/study/phase4.py
| def type_i_inflation_bound(alpha: float, se_tolerance: float, n_available: int) -> float:
"""One-sided inflation bound ``alpha + k * sqrt(alpha (1 - alpha) / n)``."""
if n_available <= 0:
raise Phase4SummaryError("n_available must be positive to compute an inflation bound.")
return float(alpha + se_tolerance * math.sqrt(alpha * (1.0 - alpha) / n_available))
|
motco.simulations.study.evaluate_phase4_gate(gate, summaries, records, *, expected_units=None)
Aggregate the predeclared gates into proceed/hold/indeterminate.
Gate multiplicity is predeclared, not discovered after the run: each control
cell contributes one one-sided test per statistic, and at a true rate of
exactly alpha each exceeds its bound with probability around 0.023. One
marginal exceedance therefore must not decide the phase. When exactly one
control statistic exceeds its bound by less than
confirmation_se_threshold Monte Carlo standard errors and no other
mandatory gate fails, the decision is indeterminate with a named
confirmation re-run. Two or more exceedances, or any exceedance of at least
one standard error, is hold.
proceed is never emitted when mandatory records or diagnostics are
incomplete.
Source code in src/motco/simulations/study/phase4.py
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353 | def evaluate_phase4_gate(
gate: Phase4GateConfig,
summaries: Sequence[SimulationSummaryResult],
records: Sequence[SimulationReplicateResult],
*,
expected_units: int | None = None,
) -> Phase4GateDecision:
"""Aggregate the predeclared gates into ``proceed``/``hold``/``indeterminate``.
Gate multiplicity is predeclared, not discovered after the run: each control
cell contributes one one-sided test per statistic, and at a true rate of
exactly alpha each exceeds its bound with probability around 0.023. One
marginal exceedance therefore must not decide the phase. When exactly one
control statistic exceeds its bound by less than
``confirmation_se_threshold`` Monte Carlo standard errors and no other
mandatory gate fails, the decision is ``indeterminate`` with a named
confirmation re-run. Two or more exceedances, or any exceedance of at least
one standard error, is ``hold``.
``proceed`` is never emitted when mandatory records or diagnostics are
incomplete.
"""
if not gate.enabled:
raise Phase4SummaryError(
"Phase 4 gate evaluation requires acceptance.gate.enabled in the study configuration."
)
operating = build_operating_frame(summaries, records)
observations: list[GateObservation] = list(evaluate_type_i_inflation(gate, summaries, records))
for rule in gate.rules:
if rule.role == "mandatory_power":
observations.append(evaluate_power_rule(gate, rule, operating))
elif rule.role == "mandatory_control":
observations.append(evaluate_control_rule(gate, rule, operating))
else:
observations.append(evaluate_descriptive_rule(rule, operating))
completeness = evaluate_completeness(records, expected_units=expected_units)
if gate.require_complete_records:
observations.append(completeness)
control_failures = [
observation
for observation in observations
if observation.kind in {"type_i_inflation", "control"} and observation.met is False
]
marginal = [
observation
for observation in control_failures
if _is_marginal(observation, gate.confirmation_se_threshold)
]
material = [observation for observation in control_failures if observation not in marginal]
other_failures = [
observation
for observation in observations
if observation.kind in {"power", "completeness"} and observation.met is False
]
missing_evidence = [
observation
for observation in observations
if observation.kind in {"power", "type_i_inflation", "control", "completeness"}
and observation.met is None
]
confirmation_runs: tuple[dict[str, Any], ...] = ()
if material or other_failures:
decision = "hold"
rationale = "Mandatory gate(s) failed: " + "; ".join(
f"{observation.rule} ({observation.detail})" for observation in material + other_failures
)
if len(marginal) > gate.max_marginal_exceedances:
rationale += (
f"; plus {len(marginal)} marginal control exceedance(s) beyond the "
f"{gate.max_marginal_exceedances} allowed"
)
elif len(marginal) > gate.max_marginal_exceedances:
decision = "hold"
rationale = (
f"{len(marginal)} control statistics exceed their inflation bounds, above the "
f"{gate.max_marginal_exceedances} marginal exceedance(s) the confirmation rule allows: "
+ "; ".join(f"{observation.rule} ({observation.detail})" for observation in marginal)
)
elif marginal:
observation = marginal[0]
decision = "indeterminate"
confirmation_runs = (
{
"cell_id": observation.cell_id,
"statistic": observation.statistic,
"trajectory_mode": observation.trajectory_mode,
"effect_size": observation.effect_size,
"reason": "single marginal control exceedance awaiting confirmation",
"action": (
f"re-run cell {observation.cell_id} with independent seeds and re-evaluate "
f"the {observation.statistic} control before any Phase 5 decision"
),
},
)
rationale = (
f"Exactly one control statistic ({observation.rule}) exceeds its bound by less than "
f"{gate.confirmation_se_threshold:g} Monte Carlo SE and no other mandatory gate fails; "
"a confirmation re-run is required before a Phase 5 decision."
)
elif missing_evidence:
decision = "indeterminate"
rationale = "Incomplete evidence: " + "; ".join(
f"{observation.rule} ({observation.detail})" for observation in missing_evidence
)
else:
decision = "proceed"
rationale = "All mandatory gates met with complete eligible diagnostics."
return Phase4GateDecision(
decision=decision,
rationale=rationale,
observations=tuple(observations),
confirmation_runs=confirmation_runs,
)
|
motco.simulations.study.Phase4GateDecision
dataclass
Aggregated Phase 4 gate outcome for the paper-grade study.
Source code in src/motco/simulations/study/phase4.py
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234 | @dataclass(frozen=True)
class Phase4GateDecision:
"""Aggregated Phase 4 gate outcome for the paper-grade study."""
decision: str
rationale: str
observations: tuple[GateObservation, ...]
confirmation_runs: tuple[dict[str, Any], ...] = ()
def to_frame(self) -> pd.DataFrame:
"""Structured rows for CSV/JSON output."""
rows = [
{
"rule": observation.rule,
"kind": observation.kind,
"trajectory_mode": observation.trajectory_mode,
"statistic": observation.statistic,
"cell_id": observation.cell_id,
"effect_size": observation.effect_size,
"met": observation.met,
"detail": observation.detail,
**{f"obs_{key}": value for key, value in observation.observations.items()},
}
for observation in self.observations
]
return pd.DataFrame(rows)
|
to_frame()
Structured rows for CSV/JSON output.
Source code in src/motco/simulations/study/phase4.py
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234 | def to_frame(self) -> pd.DataFrame:
"""Structured rows for CSV/JSON output."""
rows = [
{
"rule": observation.rule,
"kind": observation.kind,
"trajectory_mode": observation.trajectory_mode,
"statistic": observation.statistic,
"cell_id": observation.cell_id,
"effect_size": observation.effect_size,
"met": observation.met,
"detail": observation.detail,
**{f"obs_{key}": value for key, value in observation.observations.items()},
}
for observation in self.observations
]
return pd.DataFrame(rows)
|
motco.simulations.study.build_phase4_frames(gate, summaries, records, *, expected_units=None)
Build every Phase 4 structured output plus the gate decision.
Source code in src/motco/simulations/study/report.py
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114 | def build_phase4_frames(
gate: Phase4GateConfig,
summaries: Sequence[SimulationSummaryResult],
records: Sequence[SimulationReplicateResult],
*,
expected_units: int | None = None,
) -> Phase4ReportFrames:
"""Build every Phase 4 structured output plus the gate decision."""
decision = evaluate_phase4_gate(gate, summaries, records, expected_units=expected_units)
return Phase4ReportFrames(
operating=build_operating_frame(summaries, records),
geometry=summarize_realized_geometry(records),
pls_selection=summarize_pls_selection(records),
attribution=summarize_attribution(records),
localization=localize_off_diagonal(records),
gate=decision.to_frame(),
decision=decision,
)
|
motco.simulations.study.write_phase4_report(frames, out_dir)
Write the Phase 4 tables as CSV and the gate decision as JSON.
Source code in src/motco/simulations/study/report.py
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166 | def write_phase4_report(frames: Phase4ReportFrames, out_dir: Path) -> dict[str, Path]:
"""Write the Phase 4 tables as CSV and the gate decision as JSON."""
out_dir = Path(out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
paths = {
"phase4_operating": out_dir / "phase4_operating.csv",
"phase4_geometry": out_dir / "phase4_geometry.csv",
"phase4_pls_selection": out_dir / "phase4_pls_selection.csv",
"phase4_attribution": out_dir / "phase4_attribution.csv",
"phase4_localization": out_dir / "phase4_localization.csv",
"phase4_gate": out_dir / "phase4_gate.csv",
}
frames.operating.to_csv(paths["phase4_operating"], index=False)
frames.geometry.to_csv(paths["phase4_geometry"], index=False)
frames.pls_selection.to_csv(paths["phase4_pls_selection"], index=False)
frames.attribution.to_csv(paths["phase4_attribution"], index=False)
frames.localization.to_csv(paths["phase4_localization"], index=False)
frames.gate.to_csv(paths["phase4_gate"], index=False)
decision_path = out_dir / "phase4_gate_decision.json"
decision_path.write_text(
json.dumps(
{
"decision": frames.decision.decision,
"rationale": frames.decision.rationale,
"confirmation_runs": list(frames.decision.confirmation_runs),
"observations": [
{
"rule": observation.rule,
"kind": observation.kind,
"trajectory_mode": observation.trajectory_mode,
"statistic": observation.statistic,
"cell_id": observation.cell_id,
"effect_size": observation.effect_size,
"met": observation.met,
"detail": observation.detail,
"observations": observation.observations,
}
for observation in frames.decision.observations
],
},
indent=2,
default=_json_default,
)
+ "\n",
encoding="utf-8",
)
paths["phase4_gate_decision"] = decision_path
return paths
|
Render every Phase 4 figure alongside the existing study outputs.
contract is forwarded to :func:render_attribution_stability; without it
every figure is byte-for-byte the pre-contract output.
Source code in src/motco/simulations/study/report.py
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195 | def render_phase4_figures(
frames: Phase4ReportFrames, out_dir: Path, *, contract: ReportContract | None = None
) -> dict[str, Path]:
"""Render every Phase 4 figure alongside the existing study outputs.
``contract`` is forwarded to :func:`render_attribution_stability`; without it
every figure is byte-for-byte the pre-contract output.
"""
out_dir = Path(out_dir)
return {
"phase4_geometry_plot": render_geometry_checkpoints(
frames.geometry, out_dir / "phase4_geometry.png"
),
"phase4_selected_components_plot": render_selected_components(
frames.pls_selection, out_dir / "phase4_selected_components.png"
),
"phase4_attribution_stability_plot": render_attribution_stability(
frames.attribution, out_dir / "phase4_attribution_stability.png", contract=contract
),
}
|