diff --git a/src/maple/core/calibration/__init__.py b/src/maple/core/calibration/__init__.py index 7c8d77d..f3a5beb 100644 --- a/src/maple/core/calibration/__init__.py +++ b/src/maple/core/calibration/__init__.py @@ -27,7 +27,10 @@ # Enums from maple.core.calibration.enums import ( + AssayModality, Compartment, + QuantityKind, + REQUIRES_REFERENCE, ExtractionMethod, Indication, IndicationMatch, @@ -66,6 +69,11 @@ SupportType, ) +from maple.core.calibration.readout import ( + Readout, + ReadoutReference, +) + # Experimental context from maple.core.calibration.experimental_context import ( ExperimentalContext, @@ -88,6 +96,31 @@ SourceRelevanceAssessment, SubmodelInput, UncertaintyType, + DistributionShape, + ExperimentalUnitType, + ObservedDistribution, + POPULATION_SPREAD_SOURCES, + ReportedStatistic, + QuantileConvention, + SpreadSource, + StatKind, +) + +# Cohort registry +from maple.core.calibration.cohort import ( + Cohort, + CohortRegistry, + EligibilityInterval, + load_cohorts, +) +from maple.core.calibration.registry_audit import ( + RegistryProblem, + check_registries, + covariance_blocks, + find_registry_problems, + resolve_n, + warn_merged_blocks, + warn_unused_cohorts, ) # Validators @@ -184,6 +217,30 @@ "SourceType", "ExtractionMethod", "enum_field_description", + "AssayModality", + "QuantityKind", + "REQUIRES_REFERENCE", + # Reported distribution + "DistributionShape", + "ExperimentalUnitType", + "ObservedDistribution", + "POPULATION_SPREAD_SOURCES", + "ReportedStatistic", + "QuantileConvention", + "SpreadSource", + "StatKind", + # Cohort registry + "Cohort", + "CohortRegistry", + "EligibilityInterval", + "load_cohorts", + "RegistryProblem", + "check_registries", + "covariance_blocks", + "find_registry_problems", + "resolve_n", + "warn_merged_blocks", + "warn_unused_cohorts", # Scenario "Intervention", "Scenario", @@ -193,6 +250,8 @@ "ConstantSourceType", "PopulationAggregation", "Observable", + "Readout", + "ReadoutReference", "ObservableConstant", "SupportType", "Submodel", diff --git a/src/maple/core/calibration/calibration_target_models.py b/src/maple/core/calibration/calibration_target_models.py index be27893..098173d 100644 --- a/src/maple/core/calibration/calibration_target_models.py +++ b/src/maple/core/calibration/calibration_target_models.py @@ -150,7 +150,7 @@ class CalibrationTargetEstimates(BaseModel): median: List[float] = Field( description=( - "COMPUTED median value from distribution_code, as a length-1 list " "(e.g. ``[42.0]``)." + "COMPUTED median value from distribution_code, as a length-1 list (e.g. ``[42.0]``)." ) ) ci95: List[List[float]] = Field( @@ -162,10 +162,12 @@ class CalibrationTargetEstimates(BaseModel): units: str = Field(description="Units of the observable (Pint-parseable)") # Sample size metadata - sample_size: int = Field( + sample_size: Optional[int] = Field( + default=None, description=( - "Sample size (n) for the measurement. CRITICAL for uncertainty quantification " - "and pooling across studies.\n\n" + "Sample size (n) for the measurement. Mechanistic targets only: a literature " + "target names a cohort, which owns the patient count, and declaring it here " + "too lets the two disagree.\n\n" "WHERE TO LOOK:\n" "- 'n = X' or 'N = X' in methods/results sections\n" "- Sample sizes in figure legends (e.g., 'n=5 per group')\n" @@ -176,17 +178,29 @@ class CalibrationTargetEstimates(BaseModel): "- Check figure error bars - if SEM is reported, can back-calculate n from SD/SEM\n" "- Use conservative estimate based on study type\n" "- Document uncertainty in sample_size_rationale" - ) + ), ) - sample_size_rationale: str = Field( + sample_size_rationale: Optional[str] = Field( + default=None, description=( - "Explanation of how sample_size was determined.\n\n" + "Explanation of how sample_size was determined. Required with sample_size.\n\n" "Examples:\n" "- 'n=5 per group stated in Methods section 2.3'\n" "- 'n inferred from figure error bars using SEM formula: n = (SD/SEM)²'\n" "- 'n estimated as typical for this study type; not explicitly reported'\n" "- 'n=3 replicates per condition, standard for in vitro T cell assays'" - ) + ), + ) + + n_evaluable: Optional[int] = Field( + default=None, + ge=1, + description=( + "Patients with an evaluable value for THIS readout, when fewer than the cohort's " + "n_c (a panel that failed on some patients, a fold change needing both timepoints). " + "Omit when every patient in the cohort contributes; consumers then read the " + "cohort's n_c, so the two cannot drift apart." + ), ) inputs: List[EstimateInput] = Field( @@ -250,130 +264,102 @@ class CalibrationTargetEstimates(BaseModel): " 'ci95_lower': np.array(lowers) * means.units,\n" " 'ci95_upper': np.array(uppers) * means.units,\n" " }\n\n" - "POPULATION SAMPLE (for hierarchical inference), gated by population_spread:\n" - "If population_spread='across_patient', you MUST ALSO return a 'samples' key\n" - "holding the across-patient population draw (a Pint Quantity array, one value per\n" - "patient-equivalent; 1-D for a scalar observable, or 2-D (n_patients, k) for a\n" - "joint/compositional/trajectory observable). Its empirical spread (IQR/SD) is the\n" - "population-variability (omega) signal, read directly with no parametric refit.\n" - "It MUST be the POPULATION sample (dispersion = genuine patient-to-patient\n" - "variability), not the sampling distribution of an estimator (e.g. a pooled-mean /\n" - "SEM draw). median_obs/ci95 are UNAFFECTED — they remain the center + measurement\n" - "uncertainty that non-hierarchical inference reads (and may legitimately differ\n" - "from the sample spread, e.g. a meta-analytic anchor CI vs a wider population).\n" - "If population_spread='center_only' (the default), do NOT return 'samples'." + "POPULATION SAMPLE (optional):\n" + "A reported width is a ROW of the observation vector, carried by " + "observed_distribution.statistics, not an error bar — population inference " + "builds its own error bars by resampling the predicted cohort. So a " + "population spread_source does not oblige this code to return anything " + "extra.\n" + "You MAY return a 'samples' key holding the across-patient population draw " + "(a Pint Quantity array, one value per patient-equivalent; 1-D for a scalar " + "observable, or 2-D (n_patients, k) for a joint/compositional/trajectory " + "one). When present it MUST be the POPULATION sample (dispersion = genuine " + "patient-to-patient variability), not the sampling distribution of an " + "estimator (e.g. a pooled-mean / SEM draw), and its median must match the " + "reported median. For a center/technical spread_source, do NOT return it." ) ) - population_spread: Literal["across_patient", "center_only"] = Field( - default="center_only", - description=( - "What this target's reported width MEANS, and whether it feeds the population-" - "variability (omega) signal in hierarchical inference. This is an OPT-IN " - "contract: a target contributes to omega only if it explicitly declares " - "'across_patient' AND returns a population `samples` array.\n\n" - " - 'center_only' (default): the width is NOT genuine across-patient spread — " - "e.g. a pooled-mean / SEM confidence interval (shrinks with n), a transcript-vs-" - "protein method bound, or an assumed CV with no dispersion data. The target " - "still constrains the population CENTER via `median`, but is EXCLUDED from omega " - "conditioning. distribution_code must NOT return a `samples` key (that would " - "contradict the declaration).\n" - " - 'across_patient': the width is real patient-to-patient variability, usable " - "as the omega signal. distribution_code MUST return a `samples` key holding the " - "across-patient population draw; the validator rejects the target otherwise. " - "The conservative default keeps a fabricated or sample-size-shrinking width from " - "silently polluting inferred diversity — you have to assert the spread is real." - ), - ) observed_distribution: Optional[ObservedDistribution] = Field( default=None, description=( - "Optional quantile-anchor representation of the reported distribution " - "(median plus whatever scale anchors the source gives: IQR edges, quartiles, " - "deciles, or a dense empirical quantile function). This is the general, " - "shape-preserving data layer shared with SubmodelTarget; the framework can " - "derive median/IQR/scale from it on demand. ADDITIVE and OPTIONAL — when " - "omitted, the target behaves exactly as before (median/ci95 from " - "distribution_code, omega gated by population_spread). Its `spread_source` " - "carries finer provenance than the two-valued population_spread; when both " - "are present they must agree on whether the width is genuine population spread " - "(see resolved_spread_source)." + "The reported distribution: every statistic the source actually printed, as a " + "flat list. The data layer shared with " + "SubmodelTarget, and the sole declaration of whether the reported width is " + "genuine population spread (`spread_source`). Required for " + "epistemic_basis='literature'; null for 'mechanistic', which asserts a constraint " + "rather than measuring anyone." ), ) @property def resolved_spread_source(self) -> SpreadSource: - """Unified spread provenance for downstream inference. - - Prefers the richer ``observed_distribution.spread_source`` when present; - otherwise maps the legacy two-valued ``population_spread`` literal onto the - shared enum. Lets consumers read one field regardless of which layer the - target was authored against. - """ + """Spread provenance for downstream inference.""" if self.observed_distribution is not None: return self.observed_distribution.spread_source - return ( - SpreadSource.ACROSS_PATIENT - if self.population_spread == "across_patient" - else SpreadSource.CENTER_ONLY - ) + return SpreadSource.CENTER_ONLY + + @property + def feeds_population_spread(self) -> bool: + """Whether this target's width feeds the population-variability signal.""" + od = self.observed_distribution + return od is not None and od.feeds_population_spread @model_validator(mode="after") - def validate_observed_distribution_consistency(self) -> "CalibrationTargetEstimates": - """When both spread declarations are present, they must not contradict. + def validate_units_are_carried_by_the_cohort(self) -> "CalibrationTargetEstimates": + """The cohort owns the unit accounting, so the distribution must not restate it. - ``observed_distribution`` is additive, so the legacy ``population_spread`` - contract still governs the ``samples`` gate; this only forbids the two from - disagreeing about whether the reported width is genuine population spread. + ``n_biological`` / ``experimental_unit_type`` / ``unit_group`` exist for submodel + targets, which have no cohort registry. A calibration target names a cohort, and + two places for one count is how the two drift apart. """ od = self.observed_distribution if od is None: return self - feeds = od.feeds_population_spread - declared = self.population_spread == "across_patient" - if feeds != declared: + set_here = [ + name + for name in ("n_biological", "n_technical", "experimental_unit_type", "unit_group") + if getattr(od, name) is not None + ] + if od.n_biological_is_floor: + set_here.append("n_biological_is_floor") + if set_here: raise ValueError( - f"observed_distribution.spread_source='{od.spread_source.value}' " - f"({'feeds' if feeds else 'does not feed'} population spread) contradicts " - f"population_spread='{self.population_spread}'. Align them: a population " - f"spread source (across_patient / biological_experimental) requires " - f"population_spread='across_patient'; a center/technical source requires " - f"'center_only'." + f"observed_distribution sets {set_here}, which a calibration target's cohort " + "already carries (cohort_id, n_c, n_c_is_floor). Set n_evaluable on " + "empirical_data if this readout has fewer evaluable patients than the cohort. " + "Those fields are for submodel targets, which have no cohort." ) return self @model_validator(mode="after") def validate_bounded_observable_uses_logit_normal(self) -> "CalibrationTargetEstimates": - """A bounded observable's population spread (``moments`` form) must use - ``shape: logit_normal``, not normal/lognormal. - - Parity with the SubmodelTarget validator of the same name: for a - fraction / proportion / probability / percent observable, ``normal`` puts - mass outside the bound and ``lognormal`` is unbounded above; ``logit_normal`` - keeps expanded quartiles in (0, 1). Only applies to the ``moments`` form — - the ``quantiles`` form carries the empirical shape directly and is exempt. + """A bounded observable that declares a shape must use ``logit_normal``. + + For a fraction / proportion / probability / percent observable, ``normal`` + puts mass outside the bound and ``lognormal`` is unbounded above; only + ``logit_normal`` keeps derived quartiles in (0, 1). A distribution with no + declared shape is exempt: it reports its quantiles directly and nothing is + expanded. """ od = self.observed_distribution - if od is None or od.moments is None: - return self - if od.moments.shape == DistributionShape.LOGIT_NORMAL: + if od is None or od.shape is None or od.shape == DistributionShape.LOGIT_NORMAL: return self BOUNDED_UNITS = {"percent", "%", "fraction", "proportion", "probability"} if (self.units or "").strip().lower() not in BOUNDED_UNITS: return self raise ValueError( f"Observable units='{self.units}' are a bounded fraction/percentage, but " - f"observed_distribution.moments uses shape='{od.moments.shape.value}'. " - "Bounded observables must use shape='logit_normal', which expands quartiles " - "in logit space so they never escape (0, 1); normal puts mass outside the " - "bound and lognormal is unbounded above. logit_normal requires center in " - "(0, 1) with center_type='median' — express a percent as a fraction (12% -> 0.12)." + f"observed_distribution declares shape='{od.shape.value}'. Bounded observables " + "must use shape='logit_normal', which expands quartiles in logit space so they " + "never escape (0, 1); normal puts mass outside the bound and lognormal is " + "unbounded above. Express a percent as a fraction (12% -> 0.12)." ) @field_validator("sample_size") @classmethod - def validate_sample_size_positive(cls, v: int) -> int: + def validate_sample_size_positive(cls, v: Optional[int]) -> Optional[int]: """Validate sample size is at least 1.""" - if v < 1: + if v is not None and v < 1: raise ValueError(f"sample_size must be at least 1, got {v}") return v @@ -692,6 +678,17 @@ class CalibrationTarget(BaseModel): ), ) + # --- Cohort (required for literature targets) --- + cohort_id: Optional[str] = Field( + default=None, + description=( + "Registered cohort whose patients this target's statistics are computed over, from " + "the project's cohort registry. Targets sharing a cohort form one covariance block " + "in population inference. Required when epistemic_basis='literature'; null when " + "'mechanistic', which asserts a constraint rather than measuring patients." + ), + ) + # --- Sources (LLM-generated, required for literature targets) --- primary_data_source: Optional[Source] = Field( default=None, @@ -1270,6 +1267,68 @@ def validate_source_relevance_warnings(self) -> "CalibrationTarget": return self + @model_validator(mode="after") + def validate_literature_target_is_placed(self) -> "CalibrationTarget": + """A literature target must say whose patients it measured, what it measured, and + how it was reported. + + All three are meaningless for ``epistemic_basis='mechanistic'``, which asserts a + constraint rather than measuring anyone, so all three are skipped there. + """ + if self.epistemic_basis != "literature": + return self + + missing = [] + if not self.cohort_id: + missing.append( + "cohort_id: name the registered cohort these patients belong to. Targets " + "sharing a cohort form one covariance block, so an unplaced target is either " + "silently independent of its siblings or silently merged with strangers." + ) + if self.observable.readout is None: + missing.append( + "observable.readout: say what was measured. Inference builds the " + "measurement-discrepancy design from its attributes, and its species " + "composition is the row's identity." + ) + if self.empirical_data.observed_distribution is None: + missing.append( + "empirical_data.observed_distribution: list the statistics the source " + "actually printed. Without them the reported shape is whatever " + "distribution_code fitted, and an SD cannot be told from an SE." + ) + if missing: + raise ValueError( + "epistemic_basis='literature' target is missing:\n- " + "\n- ".join(missing) + ) + return self + + @model_validator(mode="after") + def validate_n_is_carried_by_the_cohort(self) -> "CalibrationTarget": + """Whoever owns the patients owns the count. + + A literature target names a cohort, whose ``n_c`` is the count; a mechanistic + one measures nobody and states its own. Consumers resolve a literature n as + ``n_evaluable`` else the cohort's ``n_c``. + """ + has_n = self.empirical_data.sample_size is not None + if self.epistemic_basis == "literature": + if has_n: + raise ValueError( + f"sample_size={self.empirical_data.sample_size} is set on a literature " + "target, but its cohort already carries n_c. Two places for one count is " + "how they drift apart. Drop it, and set empirical_data.n_evaluable if " + "fewer patients contributed to this readout than the cohort holds." + ) + elif not has_n: + raise ValueError( + "epistemic_basis='mechanistic' requires sample_size and " + "sample_size_rationale: there is no cohort to read the count from." + ) + if has_n and not self.empirical_data.sample_size_rationale: + raise ValueError("sample_size requires sample_size_rationale.") + return self + @model_validator(mode="after") def validate_context_mismatch_justified(self) -> "CalibrationTarget": """A declared context mismatch must carry its justification. @@ -1552,25 +1611,21 @@ def validate_derivation_code(self) -> "CalibrationTarget": f"reported CI95[{i}] upper ({ci_rep[1]:.4g}) within 10% tolerance" ) - # --- Population sample, gated by population_spread (hierarchical omega) --- - # 'across_patient' MUST return a 'samples' key (the population draw feeding - # omega); 'center_only' MUST NOT (a population sample would contradict the - # declaration that the width is center / measurement uncertainty). - population_spread = self.empirical_data.population_spread + # --- Population sample, when the target supplies one --- + # A reported width is a ROW of the observation vector, not an error bar, + # so a population spread_source does not oblige the target to hand over a + # population array as well. A center / technical one still must not: a + # population sample would contradict the declaration. + feeds = self.empirical_data.feeds_population_spread + declared = self.empirical_data.resolved_spread_source.value has_samples = "samples" in result - if population_spread == "across_patient" and not has_samples: - raise ReturnStructureError( - "population_spread='across_patient' requires distribution_code to " - "return a 'samples' key — the across-patient population draw used as " - "the omega signal in hierarchical inference. Return one, or set " - "population_spread='center_only' if the width is not a real spread." - ) - if population_spread == "center_only" and has_samples: + if not feeds and has_samples: raise ReturnStructureError( - "population_spread='center_only' must NOT return a 'samples' key: a " + f"observed_distribution.spread_source='{declared}' is not a population " + "spread, so distribution_code must NOT return a 'samples' key: a " "population sample contradicts the declaration that the width is " - "center / measurement uncertainty. Remove 'samples', or set " - "population_spread='across_patient' if the spread is genuine." + "center / measurement uncertainty. Remove 'samples', or declare a " + "population spread_source if the spread is genuine." ) if has_samples: samples = result["samples"] @@ -1603,7 +1658,7 @@ def validate_derivation_code(self) -> "CalibrationTarget": if np.std(finite) == 0: raise ReturnStructureError( "samples has zero variance (all identical) — not a usable population " - "spread; set population_spread='center_only' if no real spread exists" + "spread; declare a center/technical spread_source if no real spread exists" ) # Tie the declared sample to the declared center: for a scalar (1-D) # target its median must match the reported median within MC tolerance. @@ -1617,9 +1672,6 @@ def validate_derivation_code(self) -> "CalibrationTarget": f"median ({median_reported[0]:.4g}) within 10% — the declared " f"'samples' array must be the population draw the median/CI summarize" ) - self._check_center_channel_is_not_population( - finite, ci95_reported[0], self.empirical_data.sample_size - ) except CalibrationTargetValidationError: # Re-raise all our custom validation errors @@ -1662,76 +1714,6 @@ def validate_derivation_code(self) -> "CalibrationTarget": return self - @staticmethod - def _check_center_channel_is_not_population( - finite: "np.ndarray", ci95_pair: list[float], sample_size: int - ) -> None: - """The two channels must not both carry the population spread. - - Calibration-target analogue of - ``SubmodelTarget.validate_center_channel_sem_scale``. A target that declares - ``population_spread='across_patient'`` uses two channels: ``median`` + - ``ci95`` pin the CENTER (so the interval must shrink with n — a bootstrap or - SEM-scale interval on the median), and ``samples`` carries the POPULATION - spread that hierarchical inference reads as omega. Returning the population's - own 2.5th / 97.5th percentiles as ``ci95`` encodes the spread TWICE: omega - gets it from ``samples``, and the flat likelihood reads it as measurement - noise, so the target is weighted as though a single simulated patient were - allowed to land anywhere in the cohort. - - This is the ``notes/calibration`` position — biological variability is the - theta term, not the noise — made enforceable. - - The trap is ``population.summarize()``, whose ``ci95_lower`` /``ci95_upper`` - ARE ``np.percentile(samples, 2.5 / 97.5)``. It is the obvious helper to call - and it silently produces a population-scale center channel. Pair - ``population.empirical_population()`` with ``population.bootstrap_median()`` - instead, which bootstraps the median and so shrinks with n. - - Detection compares the reported interval against the sample's own 95% range. - Agreement within 15% on both edges means the center channel is the population - range. ``sample_size=1`` (a single subject) is exempt: there the two - coincide legitimately. - """ - if sample_size is not None and sample_size <= 1: - return - lo_rep, hi_rep = float(ci95_pair[0]), float(ci95_pair[1]) - lo_pop = float(np.percentile(finite, 2.5)) - hi_pop = float(np.percentile(finite, 97.5)) - if lo_rep == 0 or hi_rep == 0: - return - - def _agrees(a: float, b: float) -> bool: - scale = max(abs(a), abs(b)) - return scale > 0 and abs(a - b) <= 0.15 * scale - - if not (_agrees(lo_rep, lo_pop) and _agrees(hi_rep, hi_pop)): - return - - raise ScaleMismatchError( - f"population_spread='across_patient' but ci95 = [{lo_rep:.4g}, {hi_rep:.4g}] " - f"is the POPULATION range of 'samples' " - f"([{lo_pop:.4g}, {hi_pop:.4g}] at the 2.5th/97.5th percentiles), not the " - f"uncertainty on the center.\n\n" - "The spread is then encoded twice — once in 'samples' (the omega signal " - "hierarchical inference reads) and once in ci95 (which flat inference reads " - "as measurement noise). The target is weighted as if one simulated patient " - "could land anywhere in the cohort, so it constrains far less than its n " - f"(n={sample_size}) justifies.\n\n" - "FIX — keep 'samples' exactly as it is (that channel is correct) and give " - "the center its own interval that shrinks with n. Pass the study's real " - "sample size to summarize(); it subsample-bootstraps the median and leaves " - "'samples' untouched:\n\n" - " return pop.summarize(samples, n=int(inputs['sample_size'].magnitude))\n\n" - "For a target built from per-patient values, pop.bootstrap_median(values, " - "rng=rng) gives the same center interval directly.\n\n" - "DO NOT clear this error by setting population_spread='center_only' unless " - "the reported width was never real across-patient spread (a pooled-mean / " - "SEM interval, or an assumed CV). That switch DELETES the population " - "channel — the target stops contributing to omega. Here the spread looks " - "genuine, so the center channel is what needs fixing, not the spread one." - ) - @model_validator(mode="after") def validate_source_refs(self) -> "CalibrationTarget": """Validator: Check all source_refs in empirical_data.inputs point to defined sources.""" diff --git a/src/maple/core/calibration/cohort.py b/src/maple/core/calibration/cohort.py new file mode 100644 index 0000000..0b858a8 --- /dev/null +++ b/src/maple/core/calibration/cohort.py @@ -0,0 +1,319 @@ +"""Cohorts: one study's patients, measured once. + +A calibration target names one row of the observation vector. A cohort owns an +n, an eligibility rule and a level, but it is not the unit of independence: the +block is, and population inference forms one covariance block per block, whose +off-diagonal comes from resampling whole patients. Cohorts that share people +declare a ``PatientBlock`` saying how many; everything else is its own block. +Cross-target checks live in ``registry_audit``. +""" + +from __future__ import annotations + +from pathlib import Path +from typing import Dict, FrozenSet, List, Optional + +import yaml +from pydantic import BaseModel, ConfigDict, Field, model_validator + + +class EligibilityInterval(BaseModel): + """An inclusion criterion the study applied, as an interval on a measurement. + + Stated in the reported units of the named target, since that is what the study + screened on. + """ + + model_config = ConfigDict(extra="forbid") + + target_id: str = Field( + description="Calibration target whose observable the criterion is applied to. " + "Eligibility is stated in that measurement's reported units." + ) + lo: Optional[float] = Field( + default=None, description="Inclusive lower bound in ``units``. Omit if one-sided." + ) + hi: Optional[float] = Field( + default=None, description="Inclusive upper bound in ``units``. Omit if one-sided." + ) + units: str = Field(description="Pint-parseable units the bounds are stated in.") + rationale: str = Field(description="Where the source states the criterion.") + + @model_validator(mode="after") + def _bounds_are_usable(self) -> "EligibilityInterval": + if self.lo is None and self.hi is None: + raise ValueError( + f"EligibilityInterval on '{self.target_id}' has neither lo nor hi; it selects " + "everyone. Drop it instead." + ) + if self.lo is not None and self.hi is not None and self.lo >= self.hi: + raise ValueError( + f"EligibilityInterval on '{self.target_id}' has lo={self.lo} >= hi={self.hi}, " + "which selects nobody." + ) + return self + + +class Cohort(BaseModel): + """One study's patients. Every literature-derived target names exactly one.""" + + model_config = ConfigDict(extra="forbid") + + cohort_id: str = Field( + description="Stable identifier referenced by ``CalibrationTarget.cohort_id``. Name it " + "for the patients, not for a measurement or scenario." + ) + description: str = Field( + description="Who these patients are and what measurement occasion the statistics come from." + ) + scenarios: List[str] = Field( + description="QSP scenario(s) this cohort's targets are evaluated under. More than one " + "when the cohort reports a contrast between scenarios in the same patients, such as a " + "paired pre/post fold change." + ) + n_c: int = Field( + ge=1, + description="Patients the reported statistics are computed over. This is the resampling " + "n, not the number enrolled or screened, and never a sum across studies.", + ) + n_c_is_floor: bool = Field( + default=False, + description="True when ``n_c`` is a lower bound rather than an exact count.", + ) + source_tag: str = Field( + description="The single source reporting this cohort. A target drawing on several " + "sources is a pooled estimate, not a cohort, and must be split." + ) + eligibility: List[EligibilityInterval] = Field( + default_factory=list, + description="Inclusion criteria the study applied. Empty when it reports its full sample.", + ) + notes: Optional[str] = Field(default=None, description="Anything the fields above miss.") + + @model_validator(mode="after") + def _well_formed(self) -> "Cohort": + if not self.scenarios: + raise ValueError(f"Cohort '{self.cohort_id}' declares no scenarios.") + if len(set(self.scenarios)) != len(self.scenarios): + raise ValueError(f"Cohort '{self.cohort_id}' repeats a scenario.") + return self + + +class Stratum(BaseModel): + """Patients labelled by every cohort they belong to, and how many there are. + + One line of a partition: a patient falls in exactly one stratum, and cohort + membership follows from which strata name it. + """ + + model_config = ConfigDict(extra="forbid") + + cohorts: List[str] = Field( + min_length=1, + description="Every ``cohort_id`` these patients belong to. A single entry means " + "patients only that cohort measured.", + ) + n: int = Field(ge=1, description="How many patients. Omit the stratum rather than write 0.") + + @model_validator(mode="after") + def _cohorts_distinct(self) -> "Stratum": + if len(set(self.cohorts)) != len(self.cohorts): + raise ValueError(f"Stratum {sorted(self.cohorts)} repeats a cohort_id.") + return self + + @property + def key(self) -> FrozenSet[str]: + return frozenset(self.cohorts) + + +class PatientBlock(BaseModel): + """Cohorts that share people. Declared only where they do. + + Resampling draws a block's patients once and evaluates each cohort's rows on + its own members, which is what carries the across-cohort correlation into the + covariance. Cohorts sharing nobody need no entry. + + ``strata`` counts the sharing out. It is the Venn regions of the member + cohorts, which is a complete description of the overlap: identities never + reach the resample, so the counts are all of it. Counting patients also keeps + the description realisable, since a partition of people always describes some + people, and settles three-way memberships that pairwise overlaps leave open. + + Omit ``strata`` when the overlap is known but uncounted, which is the usual + state for two papers reporting one trial. The block still says the cohorts are + not independent; it just cannot say by how much, so a consumer building a + covariance has to fall back and report that it did. + """ + + model_config = ConfigDict(extra="forbid") + + block_id: str = Field(description="Stable identifier. Name it for the study or trial.") + description: str = Field(description="Who these patients are and how the cohorts divide them.") + cohorts: List[str] = Field( + min_length=2, + description="The ``cohort_id``s this block spans. A block records sharing between " + "cohorts, so one cohort is not a block.", + ) + strata: Optional[List[Stratum]] = Field( + default=None, + description="The partition of the block's patients. Omit when the overlap is known " + "but nobody has counted it.", + ) + notes: Optional[str] = Field( + default=None, description="Where the source states the split, or why it does not." + ) + + @model_validator(mode="after") + def _well_formed(self) -> "PatientBlock": + if len(set(self.cohorts)) != len(self.cohorts): + raise ValueError(f"Block '{self.block_id}' repeats a cohort_id.") + if self.strata is None: + return self + if not self.strata: + raise ValueError( + f"Block '{self.block_id}' has an empty strata list. Omit the key to declare an " + "uncounted overlap; an empty partition describes no patients." + ) + keys = [s.key for s in self.strata] + if len(set(keys)) != len(keys): + raise ValueError( + f"Block '{self.block_id}' repeats a stratum. Two lines naming the same cohorts " + "describe one group of patients; sum them." + ) + declared = set(self.cohorts) + named = set().union(*keys) + if named - declared: + raise ValueError( + f"Block '{self.block_id}' has strata naming {sorted(named - declared)}, which " + "are not in its cohorts." + ) + if declared - named: + raise ValueError( + f"Block '{self.block_id}' spans {sorted(declared - named)} but no stratum places " + "them. Every member needs patients, or it does not belong to the block." + ) + return self + + @property + def is_quantified(self) -> bool: + """Whether the overlap is counted, and so whether a joint resample is defined.""" + return self.strata is not None + + @property + def members(self) -> FrozenSet[str]: + return frozenset(self.cohorts) + + @property + def n_patients(self) -> Optional[int]: + """Patients in the block, which is the number a joint resample draws.""" + return None if self.strata is None else sum(s.n for s in self.strata) + + def size_of(self, cohort_id: str) -> Optional[int]: + """Block patients belonging to ``cohort_id``. Equals that cohort's ``n_c``.""" + if self.strata is None: + return None + return sum(s.n for s in self.strata if cohort_id in s.key) + + def overlap(self, a: str, b: str) -> Optional[int]: + """Patients in both cohorts. Zero is a fact, not a missing value.""" + if self.strata is None: + return None + return sum(s.n for s in self.strata if a in s.key and b in s.key) + + +class CohortRegistry(BaseModel): + """The declared cohorts of a project, and the blocks over them. + + Counted blocks partition: a cohort belongs to at most one, since its patients + can be divided up only once. Uncounted blocks overlay, so a cohort may sit in + any number of them. That is what lets a study whose internal sharing is + reported sit inside a wider overlap nobody has counted. + """ + + model_config = ConfigDict(extra="forbid") + + cohorts: List[Cohort] = Field(default_factory=list) + blocks: List[PatientBlock] = Field( + default_factory=list, + description="Groups of cohorts that share patients. Cohorts sharing nobody are their " + "own block and need no entry.", + ) + + @model_validator(mode="after") + def _ids_unique_and_blocks_resolve(self) -> "CohortRegistry": + counts: Dict[str, int] = {} + for c in self.cohorts: + counts[c.cohort_id] = counts.get(c.cohort_id, 0) + 1 + dupes = sorted(k for k, v in counts.items() if v > 1) + if dupes: + raise ValueError(f"Duplicate cohort_id(s) in registry: {dupes}") + + block_ids = [b.block_id for b in self.blocks] + if len(set(block_ids)) != len(block_ids): + raise ValueError( + f"Duplicate block_id(s) in registry: " + f"{sorted({b for b in block_ids if block_ids.count(b) > 1})}" + ) + + by_id = {c.cohort_id: c for c in self.cohorts} + claimed: Dict[str, str] = {} + for b in self.blocks: + missing = sorted(set(b.cohorts) - set(by_id)) + if missing: + raise ValueError( + f"Block '{b.block_id}' names {missing}, which are not in the registry." + ) + if not b.is_quantified: + continue + for cid in b.cohorts: + # One partition per cohort: two would each claim all its patients. + if cid in claimed: + raise ValueError( + f"Cohort '{cid}' is counted by blocks '{claimed[cid]}' and " + f"'{b.block_id}'. Its patients can be divided up once; leave one " + "block's strata off to declare an uncounted overlap instead." + ) + claimed[cid] = b.block_id + size = b.size_of(cid) + if size != by_id[cid].n_c: + raise ValueError( + f"Block '{b.block_id}' places {size} patients in cohort '{cid}', which " + f"declares n_c={by_id[cid].n_c}. The strata are that cohort's patients, " + "so they have to add up to it." + ) + return self + + def get(self, cohort_id: str) -> Optional[Cohort]: + return self.as_dict().get(cohort_id) + + def as_dict(self) -> Dict[str, Cohort]: + return {c.cohort_id: c for c in self.cohorts} + + def counted_block_for(self, cohort_id: str) -> Optional[PatientBlock]: + """The block whose strata divide this cohort's patients, if one does.""" + for b in self.blocks: + if b.is_quantified and cohort_id in b.members: + return b + return None + + def uncounted_blocks_for(self, cohort_id: str) -> List[PatientBlock]: + """Declared overlaps involving this cohort that nobody has counted.""" + return [b for b in self.blocks if not b.is_quantified and cohort_id in b.members] + + @property + def uncounted_blocks(self) -> List[PatientBlock]: + """Every declared overlap without strata. A consumer building a covariance + cannot honour these and has to report that it drew their cohorts apart.""" + return [b for b in self.blocks if not b.is_quantified] + + +def load_cohorts(path: Path | str) -> CohortRegistry: + """Load a cohort registry YAML: a mapping with ``cohorts:`` and optional ``blocks:``, + or a bare list of cohorts.""" + path = Path(path) + if not path.exists(): + raise FileNotFoundError(f"Cohort registry not found: {path}") + raw = yaml.safe_load(path.read_text()) or {} + if isinstance(raw, list): + raw = {"cohorts": raw} + return CohortRegistry.model_validate(raw) diff --git a/src/maple/core/calibration/cross_scenario_loader.py b/src/maple/core/calibration/cross_scenario_loader.py index 6aaa405..1ad0bda 100644 --- a/src/maple/core/calibration/cross_scenario_loader.py +++ b/src/maple/core/calibration/cross_scenario_loader.py @@ -153,6 +153,8 @@ def load_cross_scenario_targets(yaml_dir: Path | str) -> pd.DataFrame: ci95_lower = float("nan") ci95_upper = float("nan") + # A literature contrast has no single n: each arm is its own cohort, and + # they are resampled independently. Only a mechanistic target states one. sample_size = empirical.get("sample_size", float("nan")) if isinstance(sample_size, list): # Vector sample_size doesn't apply to cross-scenario scalar diff --git a/src/maple/core/calibration/cross_scenario_target.py b/src/maple/core/calibration/cross_scenario_target.py index 2001e80..715f4ab 100644 --- a/src/maple/core/calibration/cross_scenario_target.py +++ b/src/maple/core/calibration/cross_scenario_target.py @@ -47,6 +47,8 @@ from maple.core.calibration.calibration_target_models import ( CalibrationTargetEstimates, ) +from maple.core.calibration.enums import QuantityKind +from maple.core.calibration.readout import Readout from maple.core.calibration.shared_models import SecondarySource, Source @@ -124,6 +126,27 @@ def compute_test_statistic(time, species_dict) -> float "required_species (series / param / template / missing)." ), ) + readout: Readout = Field( + description="What this arm measures. Its declared composition is what lets an arm " + "be compared against a standalone target on the same cohort." + ) + cohort_id: Optional[str] = Field( + default=None, + description=( + "Registered cohort whose patients this arm measures. Required when the " + "target is epistemic_basis='literature'. Two roles may not name one cohort: " + "the same patients under both conditions is a paired contrast, which belongs " + "in a single CalibrationTarget with an Observable.reference." + ), + ) + n_evaluable: Optional[int] = Field( + default=None, + ge=1, + description=( + "Patients with an evaluable value for this arm, when fewer than the cohort's " + "n_c. Omit when every patient in the cohort contributes." + ), + ) class CrossScenarioObservable(BaseModel): @@ -163,6 +186,13 @@ class CrossScenarioObservable(BaseModel): "(event times), 'cell/mm**2' (density differences)." ) ) + quantity_kind: QuantityKind = Field( + description=( + "What kind of quantity the reduction produces. Each arm declares its own " + "readout; this is the composed output. No `reference` is needed: the roles " + "and their cohorts are the reference." + ) + ) inputs: List[CrossScenarioInput] = Field( min_length=2, description=( @@ -172,6 +202,17 @@ class CrossScenarioObservable(BaseModel): ), ) + @model_validator(mode="after") + def validate_arms_share_an_assay(self) -> "CrossScenarioObservable": + """The arms are one measurement made in different people.""" + modalities = sorted({i.readout.assay_modality.value for i in self.inputs}) + if len(modalities) > 1: + raise ValueError( + f"Arms declare different assay_modality values {modalities}. A contrast " + "between assays measures the assays, not the treatment." + ) + return self + @model_validator(mode="after") def validate_unique_roles(self) -> "CrossScenarioObservable": roles = [inp.role for inp in self.inputs] @@ -183,8 +224,7 @@ def validate_unique_roles(self) -> "CrossScenarioObservable": dupes.append(r) seen.add(r) raise ValueError( - f"CrossScenarioObservable roles must be unique; " - f"duplicates: {sorted(set(dupes))}" + f"CrossScenarioObservable roles must be unique; duplicates: {sorted(set(dupes))}" ) return self @@ -269,7 +309,7 @@ class CrossScenarioCalibrationTarget(BaseModel): primary_data_source: Optional[Source] = Field( default=None, description=( - "Required when epistemic_basis='literature'. May be null for " "mechanistic targets." + "Required when epistemic_basis='literature'. May be null for mechanistic targets." ), ) secondary_data_sources: List[SecondarySource] = Field( @@ -295,6 +335,40 @@ def validate_epistemic_basis_consistency(self) -> "CrossScenarioCalibrationTarge ) return self + @model_validator(mode="after") + def validate_literature_arms_are_placed(self) -> "CrossScenarioCalibrationTarget": + """A measured contrast must say whose patients each arm is. + + The arms are different people, so the contrast is not a row of any one + cohort: it is a derived row over the cohorts it draws on, and inference + resamples each of them independently to get its variance. That is only + possible once each arm names its own. + """ + if self.epistemic_basis != "literature": + return self + unplaced = [i.role for i in self.observable.inputs if not i.cohort_id] + if unplaced: + raise ValueError( + f"epistemic_basis='literature' but these roles name no cohort: {unplaced}. " + "Each arm measures a different set of patients and must name the registered " + "cohort they belong to." + ) + cohorts = [i.cohort_id for i in self.observable.inputs] + repeated = sorted({c for c in cohorts if cohorts.count(c) > 1}) + if repeated: + raise ValueError( + f"Roles share cohort(s) {repeated}. One cohort under both conditions is the " + "same patients measured twice, so the contrast is paired and belongs in a " + "single CalibrationTarget whose observable declares a reference." + ) + if self.empirical_data.sample_size is not None: + raise ValueError( + f"sample_size={self.empirical_data.sample_size} is set on a literature " + "contrast, but its arms are different people with their own cohorts. A sum " + "across arms is not a resampling n. Declare n_evaluable per role instead." + ) + return self + @model_validator(mode="after") def validate_units_match_empirical(self) -> "CrossScenarioCalibrationTarget": if self.observable.units != self.empirical_data.units: diff --git a/src/maple/core/calibration/denominator_audit.py b/src/maple/core/calibration/denominator_audit.py new file mode 100644 index 0000000..230fae3 --- /dev/null +++ b/src/maple/core/calibration/denominator_audit.py @@ -0,0 +1,381 @@ +"""Cross-target denominator checks. + +Every other calibration validator is single-target: pydantic hands a model +validator one object, so it can only ask questions about that object. Some +defects are not properties of a target at all — they are properties of a *pair*. +This module holds the checks that need the whole loaded set, and is called from +:func:`maple.core.calibration.test_stats_loader.load_calibration_targets`. + +The motivating case (pdac, July 2026). Two targets computed the identical model +expression, ``Treg / (Treg + Th + Th_exh)``, and asserted different values for +it — 0.50 from a panel that counted Tregs against *polarised* CD4 only, 0.34 +from one that counted against *all* CD4. Both declared the denominator audit +correctly. Both wrote the contradiction out in English, in adjacent fields: + + treg_fraction_cd4 "...(Treg + Th1 + Th2 + Th17), excluding Th0 bystanders" + treg_fraction_hiraoka "CD4+ tumor-infiltrating T lymphocytes..." + +One says *excluding*, the other says *includes*. Nothing compared them, so both +entered the likelihood at full weight and became the top two misfit drivers — +57% of joint influence off 15% of the observables. The extraction agent did its +job; the pipeline had no place to notice. + +Note that the CIs of that pair *overlap*, so a value-agreement check would not +have caught it. The signal is the denominator prose, not the numbers. +""" + +from __future__ import annotations + +import ast +import re +import warnings +from dataclasses import dataclass, field +from typing import Any, Dict, List, Optional, Set, Tuple + + +# --------------------------------------------------------------------------- # +# Reading the observable's numerator and denominator out of its code # +# --------------------------------------------------------------------------- # +def _species_reached(node: ast.AST, env: Dict[str, Set[str]]) -> Set[str]: + """Species-dict keys reachable from ``node``, resolving local variables.""" + found: Set[str] = set() + for sub in ast.walk(node): + if ( + isinstance(sub, ast.Subscript) + and isinstance(sub.value, ast.Name) + and sub.value.id == "species_dict" + and isinstance(sub.slice, ast.Constant) + ): + found.add(sub.slice.value) + elif isinstance(sub, ast.Name) and sub.id in env: + found |= env[sub.id] + return found + + +def numerator_and_denominator(code: str) -> Tuple[Set[str], Set[str]]: + """Species on each side of the division(s) in an ``observable.code`` body. + + Returns ``(numerator, denominator)``. Both are empty when the code does not + divide, or fails to parse (syntax is reported by the schema validators). + + Real observables bind intermediates first — ``total_t = treg + cd8 + th``, + then ``treg / total_t`` — so a single forward pass records what each local + name reaches before the divisions are read. Names are resolved + transitively, which covers the chained case. + + An observable that divides more than once (a ratio of two fractions) unions + both sides; that is coarse but conservative, and no live target does it. + """ + try: + tree = ast.parse(code) + except SyntaxError: + return set(), set() + + env: Dict[str, Set[str]] = {} + for node in ast.walk(tree): + if ( + isinstance(node, ast.Assign) + and len(node.targets) == 1 + and isinstance(node.targets[0], ast.Name) + ): + env[node.targets[0].id] = _species_reached(node.value, env) + + numerator: Set[str] = set() + denominator: Set[str] = set() + for node in ast.walk(tree): + if isinstance(node, ast.BinOp) and isinstance(node.op, ast.Div): + numerator |= _species_reached(node.left, env) + denominator |= _species_reached(node.right, env) + return numerator, denominator + + +# --------------------------------------------------------------------------- # +# A — cross-target mapping collisions # +# --------------------------------------------------------------------------- # +@dataclass(frozen=True) +class MappingCollision: + """Two or more targets reducing to one model quantity.""" + + numerator: Tuple[str, ...] + denominator: Tuple[str, ...] + readout_time: Optional[float] + members: Tuple[str, ...] + experimental_denominators: Tuple[str, ...] + justified_by: Tuple[str, ...] = field(default=()) + + @property + def is_justified(self) -> bool: + return bool(self.justified_by) + + +def _mapping_key(target: Dict[str, Any]) -> Optional[tuple]: + observable = target.get("observable") or {} + readout = observable.get("readout") or {} + denominator = readout.get("denominator_species") or [] + if not denominator: + return None + return ( + tuple(sorted(readout.get("numerator_species") or [])), + tuple(sorted(denominator)), + observable.get("readout_time"), + ) + + +def find_mapping_collisions(targets: Dict[str, Dict[str, Any]]) -> List[MappingCollision]: + """Group ``{target_id: parsed_yaml}`` by model quantity; return the groups of >1. + + The key is (numerator species, denominator species, readout time) as declared + on the readout; ``check_code_matches_readout`` separately holds the code to that + declaration. Including the numerator matters: nine density targets divide by + ``V_T`` and are not in conflict, because they count different cells. + + Callers pass one scenario at a time. Across scenarios the same expression is + expected (a baseline and a day-21 arm), and is not a collision. + """ + groups: Dict[tuple, List[str]] = {} + for target_id, data in targets.items(): + key = _mapping_key(data) + if key is None: + continue + groups.setdefault(key, []).append(target_id) + + collisions = [] + for key, members in sorted(groups.items()): + if len(members) < 2: + continue + members = sorted(members) + observables = [(targets[m].get("observable") or {}) for m in members] + collisions.append( + MappingCollision( + numerator=key[0], + denominator=key[1], + readout_time=key[2], + members=tuple(members), + experimental_denominators=tuple( + (o.get("readout") or {}).get("experimental_denominator") or "" + for o in observables + ), + justified_by=tuple( + m + for m, o in zip(members, observables) + if (o.get("duplicate_mapping_justification") or "").strip() + ), + ) + ) + return collisions + + +def check_mapping_collisions(targets: Dict[str, Dict[str, Any]]) -> None: + """Raise when two targets compute one model quantity without saying why. + + No similarity threshold. Prose is too weak a discriminator to tune on: in + the pdac corpus two targets phrasing the SAME denominator differently + ("mm^2 of tumor tissue section" vs "mm^2 of intratumoral tumor section area, + pooled across the three cohorts") score 0.38 on token overlap, while the + genuinely contradictory Treg pair scores 0.06. Any cut between those is + fitted to one example. So a collision is always surfaced, and the author + declares which kind it is via ``duplicate_mapping_justification``. + + That is the right default anyway: two targets reducing to one quantity is a + claim worth stating explicitly, whether or not it is sound. + """ + unjustified = [c for c in find_mapping_collisions(targets) if not c.is_justified] + if not unjustified: + return + + blocks = [] + for c in unjustified: + lines = [ + f" {' + '.join(c.numerator) or ''} / {' + '.join(c.denominator)}" + + (f" at t={c.readout_time}" if c.readout_time is not None else ""), + ] + for member, exp in zip(c.members, c.experimental_denominators): + lines.append(f" - {member}") + lines.append(f" experimental_denominator: {exp or ''}") + blocks.append("\n".join(lines)) + + raise ValueError( + "Calibration targets compute the same model quantity without declaring why:\n\n" + + "\n\n".join(blocks) + + "\n\nEach group reduces to ONE model expression, so its members make one " + "claim several times over and each enters the likelihood at full weight. " + "Compare the experimental_denominator lines above:\n\n" + " - If they describe DIFFERENT experimental quantities, one of the model " + "mappings is wrong. Re-derive it against a denominator that matches what " + "the paper actually measured. This is the real defect the check exists " + "for — the pdac Treg pair mapped 'Tregs among polarised CD4' and 'Tregs " + "among all CD4' onto the same expression, and the disagreement showed up " + "as the top two misfit drivers instead of as a extraction bug.\n" + " - If they describe the SAME quantity in different cohorts, set " + "observable.duplicate_mapping_justification on at least one of them, " + "naming the cohorts. Then check the values are mutually compatible; a " + "disagreement beyond the CIs is a cross-study conflict to adjudicate, " + "not something to average away." + ) + + +# --------------------------------------------------------------------------- # +# C — a declared denominator bias should be heard # +# --------------------------------------------------------------------------- # +#: A declaration opening with "None" / "N/A" asserts there is no bias — the field +#: docstring offers exactly that ("None (denominator fully captured by model +#: species)"), so it is an answer, not an omission. +_NO_BIAS = re.compile(r"^\s*(none|n/?a)\b", re.I) + +#: A usable declaration says which way the bias runs, or how big it is. Prose is +#: matched loosely on purpose: this decides whether to *mention* a target, never +#: whether to reject one. +_DIRECTION_OR_MAGNITUDE = re.compile( + r"""( + \d+\s*[-–—]\s*\d+\s*(x|%|-?fold) # a range: "50-70%", "2-3x" + | \d+(\.\d+)?\s*(x|%|-?fold)\b # a multiplier: "3x", "10-fold" + | \b(higher|lower)\b + | \b(over|under)-?(predict|estimat) # "overpredict", "under-estimate" + | \binflat # "inflates the denominator" + | \bbias(es|ed)?\s+(up|down) + )""", + re.I | re.X, +) + + +@dataclass(frozen=True) +class DeclaredDenominatorBias: + target_id: str + experimental_denominator: str + unmodeled_components: str + states_direction_or_magnitude: bool + + +def collect_declared_biases( + targets: Dict[str, Dict[str, Any]], +) -> List[DeclaredDenominatorBias]: + """Targets whose own ``unmodeled_denominator_components`` admits a bias. + + Declarations that open with "None" are skipped — they assert the denominator + is fully captured, which is an answer rather than an omission. + """ + out = [] + for target_id in sorted(targets): + observable = targets[target_id].get("observable") or {} + unmodeled = (observable.get("unmodeled_denominator_components") or "").strip() + if not unmodeled or _NO_BIAS.match(unmodeled): + continue + out.append( + DeclaredDenominatorBias( + target_id=target_id, + experimental_denominator=( + (observable.get("readout") or {}).get("experimental_denominator") or "" + ).strip(), + unmodeled_components=unmodeled, + states_direction_or_magnitude=bool(_DIRECTION_OR_MAGNITUDE.search(unmodeled)), + ) + ) + return out + + +def warn_declared_biases(targets: Dict[str, Dict[str, Any]]) -> List[DeclaredDenominatorBias]: + """Surface declared denominator biases once per load; return them for diagnostics. + + ``unmodeled_denominator_components`` asks the author to document "expected + direction and magnitude of systematic bias", and authors fill it in. It then + had no consequence anywhere: the pdac Hiraoka target declared that its + experimental denominator includes bystander CD4+ cells the model does not + represent, and still entered the likelihood at full weight, unweighted and + unflagged. A field only a reader can act on is not a control. + + **One** warning, not one per target. In the pdac corpus 32 of 52 targets + declare a bias; a per-target warning would be 32 lines on every load and + would be tuned out within a week. + + Deliberately not an error, on two counts. Whether a declared bias should + downweight a target or gate it out of the joint is an inference-side policy + decision, not a schema one — hence the returned list, so the caller can act. + And measured against the live corpus, only 2 of 32 declarations state a + direction or magnitude, so requiring one would block 30 targets. That is a + real quality gap and a fair thing to enforce later, but it is a corpus + migration, not a validator flip. + """ + biases = collect_declared_biases(targets) + if not biases: + return biases + + vague = [b for b in biases if not b.states_direction_or_magnitude] + message = ( + f"{len(biases)} calibration target(s) declare an unmodeled denominator " + "component. Each enters the likelihood at full weight regardless — decide " + "whether they should be downweighted or reconciled first." + ) + if vague: + listing = "\n".join(f" - {b.target_id}" for b in vague) + message += ( + f"\n\n{len(vague)} of them name the missing components but never state the " + "expected DIRECTION or MAGNITUDE of the bias, which is what " + "unmodeled_denominator_components asks for. Without it the declaration " + "cannot be acted on, only read:\n" + f"{listing}" + ) + warnings.warn(message, UserWarning) + return biases + + +@dataclass(frozen=True) +class CodeReadoutMismatch: + """A target whose code divides by something other than its declared readout.""" + + target_id: str + declared: Tuple[str, ...] + in_code: Tuple[str, ...] + + +def find_code_readout_mismatches( + targets: Dict[str, Dict[str, Any]], +) -> List[CodeReadoutMismatch]: + """Targets whose ``observable.code`` disagrees with ``readout.denominator_species``. + + The declaration carries the row's identity, so the code has to compute what it + says. Only checked where the code divides: an observable that reaches its + denominator some other way (an area built from a constant) is not comparable + this way. + + A sum the model defines as an aggregate is named by that aggregate on both + sides, so this is a plain comparison. Listing the members instead is what + goes stale when a pool joins the rule, and is what this then reports. + """ + out = [] + for tid, data in sorted(targets.items()): + observable = data.get("observable") or {} + declared = set((observable.get("readout") or {}).get("denominator_species") or []) + _, in_code = numerator_and_denominator(observable.get("code") or "") + if not declared or not in_code or declared == in_code: + continue + out.append(CodeReadoutMismatch(tid, tuple(sorted(declared)), tuple(sorted(in_code)))) + return out + + +def warn_code_readout_mismatches( + targets: Dict[str, Dict[str, Any]], +) -> List[CodeReadoutMismatch]: + """Warn on each mismatch. Returns them.""" + found = find_code_readout_mismatches(targets) + for m in found: + warnings.warn( + f"{m.target_id}: readout.denominator_species={list(m.declared)} but the code " + f"divides by {list(m.in_code)}. Name the model aggregate that defines the sum, " + "on both sides.", + UserWarning, + ) + return found + + +__all__ = [ + "numerator_and_denominator", + "CodeReadoutMismatch", + "find_code_readout_mismatches", + "warn_code_readout_mismatches", + "MappingCollision", + "find_mapping_collisions", + "check_mapping_collisions", + "DeclaredDenominatorBias", + "collect_declared_biases", + "warn_declared_biases", +] diff --git a/src/maple/core/calibration/enums.py b/src/maple/core/calibration/enums.py index 41db654..e2834c2 100644 --- a/src/maple/core/calibration/enums.py +++ b/src/maple/core/calibration/enums.py @@ -200,6 +200,77 @@ class ExtractionMethod(str, Enum): """Other extraction method (specify in extraction_notes).""" +class QuantityKind(str, Enum): + """What kind of quantity an observable is. + + Read by inference to build the measurement-discrepancy design matrix. Two + observables of the same kind and modality share a row of it, and so share a + location and scale correction. + """ + + DENSITY = "density" + """Cells or mass per unit area or volume of tissue.""" + + FRACTION = "fraction" + """A part over a whole containing it, bounded in (0, 1).""" + + RATIO = "ratio" + """One quantity over another not containing it. Unbounded above.""" + + FOLDCHANGE = "foldchange" + """A quantity relative to its own value at a reference. Needs ``reference``.""" + + CONCENTRATION = "concentration" + """Amount of a soluble species per unit volume or per unit protein.""" + + INTENSITY = "intensity" + """A normalised or relative signal level with no absolute scale.""" + + TIME = "time" + """A duration.""" + + +#: Kinds that are meaningless without ``Observable.reference``. +REQUIRES_REFERENCE = frozenset({QuantityKind.FOLDCHANGE}) + + +class AssayModality(str, Enum): + """How an observable was measured.""" + + MIHC = "mihc" + """Multiplex immunohistochemistry with digital image analysis.""" + + IHC = "ihc" + """Single-plex chromogenic immunohistochemistry.""" + + MIF = "mif" + """Multiplex or double immunofluorescence.""" + + FLOW_CYTOMETRY = "flow_cytometry" + """Flow cytometry of digested tissue or whole blood.""" + + SCRNASEQ = "scrnaseq" + """Single-cell RNA sequencing.""" + + MULTIPLEX_IMMUNOASSAY = "multiplex_immunoassay" + """Bead-based multiplex immunoassay on tissue homogenate or serum.""" + + ELISA = "elisa" + """Single-analyte enzyme-linked immunosorbent assay.""" + + MULTIPHOTON_MICROSCOPY = "multiphoton_microscopy" + """Label-free multiphoton microscopy (SHG / TPEF).""" + + DIGITAL_HISTOPATHOLOGY = "digital_histopathology" + """Digital image analysis of stained whole sections or cores.""" + + IMAGING = "imaging" + """Clinical imaging: CT, MRI, ultrasound.""" + + CLINICAL = "clinical" + """Recorded clinical observation rather than an assay.""" + + # ============================================================================= # SOURCE RELEVANCE ENUMS # ============================================================================= diff --git a/src/maple/core/calibration/observable.py b/src/maple/core/calibration/observable.py index 5a1a1fe..14dae3c 100644 --- a/src/maple/core/calibration/observable.py +++ b/src/maple/core/calibration/observable.py @@ -10,12 +10,15 @@ from pydantic import BaseModel, ConfigDict, Field, model_validator -from maple.core.calibration.enums import ExtractionMethod, SourceType +from maple.core.calibration.enums import ( + ExtractionMethod, + SourceType, +) +from maple.core.calibration.readout import Readout # Import SubmodelInput directly (not under TYPE_CHECKING) so Pydantic can resolve it from maple.core.calibration.shared_models import SubmodelInput - # Support types for measurement output constraints SupportType = Literal["positive", "non_negative", "unit_interval", "positive_unbounded", "real"] @@ -384,6 +387,14 @@ class Observable(BaseModel): description="Pint-parseable units of the observable output (must match empirical_data.units)" ) + readout: Optional[Readout] = Field( + default=None, + description="What the experiment measured: quantity kind, assay, denominator, and " + "reference. The rest of this model is how to compute it. Required for a literature " + "target; null for a mechanistic one, which asserts a constraint rather than running " + "an assay.", + ) + readout_time: Optional[float] = Field( default=None, description=( @@ -473,30 +484,6 @@ class Observable(BaseModel): ), ) - experimental_denominator: Optional[str] = Field( - default=None, - description=( - "What the experimental measurement divides by. Required when the " - "observable is a density or fraction.\n\n" - "Examples:\n" - "- 'mm^2 of tumor tissue (whole section including stroma)'\n" - "- 'all cells in ROI (all nucleated cells)'\n" - "- 'CD3+ T cells (pan-T-cell marker)'\n" - "- 'CD45+ leukocytes'" - ), - ) - - model_denominator_species: Optional[List[str]] = Field( - default=None, - description=( - "Which model species compose the denominator in the observable code. " - "Required when experimental_denominator is set.\n" - "Format: 'compartment.species' (e.g., ['V_T.CD8', 'V_T.Th', 'V_T.Treg']).\n" - "For area-based denominators, list the species used to compute area " - "(e.g., ['V_T.C1'] when tumor area = C1 * area_per_cell)." - ), - ) - unmodeled_denominator_components: Optional[str] = Field( default=None, description=( @@ -511,14 +498,41 @@ class Observable(BaseModel): ), ) + duplicate_mapping_justification: Optional[str] = Field( + default=None, + description=( + "Why it is legitimate for this target to compute the SAME model " + "expression (same numerator species over the same denominator species, " + "at the same readout time) as another target in the same scenario.\n\n" + "Required only when such a collision exists — " + "``check_mapping_collisions`` raises otherwise. Two targets that reduce " + "to one model quantity are making one claim twice, and there are only " + "two ways that is sound:\n\n" + " - REPLICATE: different cohorts measuring the same thing. State the " + "cohorts and note that the two values must be mutually compatible; if " + "they disagree beyond their CIs that is a real cross-study conflict to " + "adjudicate, not something to average away.\n" + " - DELIBERATE POOLING: the targets are known duplicates kept on " + "purpose (e.g. a sensitivity check). Say so.\n\n" + "If instead the two targets measure genuinely DIFFERENT experimental " + "quantities, this field is the wrong fix — one of the model mappings is " + "wrong. That was the Treg case: one source counted Treg over " + "polarised CD4 only, the other over all CD4, and both mapped to " + "``Treg / (Treg + Th + Th_exh)``." + ), + ) + @model_validator(mode="after") def validate_denominator_fields(self) -> "Observable": """Validate denominator audit fields for density/fraction observables.""" - if self.experimental_denominator and not self.model_denominator_species: + if self.readout is None: + return self + + if self.readout.experimental_denominator and not self.readout.denominator_species: raise ValueError( - f"Observable has experimental_denominator='{self.experimental_denominator}' " - f"but model_denominator_species is not set. Specify which model species " - f"compose the denominator to complete the denominator audit." + f"readout.experimental_denominator='{self.readout.experimental_denominator}' " + "but readout.denominator_species is empty. Name the model species that " + "compose the denominator." ) # Density observables (units like cell/mm**2) must have denominator audit @@ -527,10 +541,10 @@ def validate_denominator_fields(self) -> "Observable": and "/" in self.units and self.support in ("positive", "non_negative") ) - if is_density and not self.experimental_denominator: + if is_density and not self.readout.experimental_denominator: raise ValueError( f"Observable with units='{self.units}' and support='{self.support}' " - f"is a density but experimental_denominator is not set. " + f"is a density but readout.experimental_denominator is not set. " f"Document what the experiment divides by (e.g., 'mm^2 of tumor " f"tissue including stroma') to enable denominator audit." ) @@ -561,9 +575,7 @@ def validate_reduction_choice(self) -> "Observable": ) if has_time and not self.readout_time_unit: - raise ValueError( - "readout_time_unit is required when readout_time is set " "(e.g. 'day')." - ) + raise ValueError("readout_time_unit is required when readout_time is set (e.g. 'day').") if not has_time and self.readout_time_unit: raise ValueError( "readout_time_unit is only meaningful with readout_time; " @@ -623,8 +635,7 @@ class SubmodelStateVariable(BaseModel): figure_id: Optional[str] = Field( None, description=( - "Figure identifier (e.g., 'Figure 2A', 'Fig. 3B'). " - "Required when source_type='figure'." + "Figure identifier (e.g., 'Figure 2A', 'Fig. 3B'). Required when source_type='figure'." ), ) diff --git a/src/maple/core/calibration/population.py b/src/maple/core/calibration/population.py index 567b733..62cb366 100644 --- a/src/maple/core/calibration/population.py +++ b/src/maple/core/calibration/population.py @@ -310,7 +310,7 @@ def summarize(samples, *, n=None, rng=None, n_boot=2_000): median / ci95 reduce over the patient axis (axis 0). **Pass ``n`` (the study's real sample size) whenever the target declares - ``population_spread='across_patient'``.** The two channels mean different things: + a population ``spread_source``.** The two channels mean different things: ``samples`` carries the POPULATION spread that hierarchical inference reads as omega, while ``median`` / ``ci95`` pin the CENTER, so the interval must shrink with n. Without ``n`` this returns the population's own 2.5th / 97.5th diff --git a/src/maple/core/calibration/readout.py b/src/maple/core/calibration/readout.py new file mode 100644 index 0000000..5f77ba7 --- /dev/null +++ b/src/maple/core/calibration/readout.py @@ -0,0 +1,93 @@ +"""The measurement a target reports, separate from how the model computes it.""" + +from __future__ import annotations + +from typing import List, Literal, Optional + +from pydantic import BaseModel, ConfigDict, Field, model_validator + +from maple.core.calibration.enums import REQUIRES_REFERENCE, AssayModality, QuantityKind + + +class ReadoutReference(BaseModel): + """What a relative readout is measured against.""" + + model_config = ConfigDict(extra="forbid") + + kind: Literal["timepoint", "scenario"] = Field( + description="Whether the reference is another time in this trajectory or another scenario." + ) + timepoint: Optional[float] = Field(default=None, description="Reference time.") + timepoint_unit: Optional[str] = Field( + default=None, description="Pint-parseable unit for ``timepoint``." + ) + scenario: Optional[str] = Field(default=None, description="Reference scenario name.") + + @model_validator(mode="after") + def _matches_kind(self) -> "ReadoutReference": + if self.kind == "timepoint": + if self.timepoint is None or not self.timepoint_unit: + raise ValueError("reference kind='timepoint' needs timepoint and timepoint_unit.") + if self.scenario is not None: + raise ValueError("reference kind='timepoint' must not set scenario.") + else: + if not self.scenario: + raise ValueError("reference kind='scenario' needs scenario.") + if self.timepoint is not None: + raise ValueError("reference kind='scenario' must not set timepoint.") + return self + + +class Readout(BaseModel): + """What was measured. Inference builds the measurement-discrepancy design from + these attributes, so two readouts agreeing on them share one correction.""" + + model_config = ConfigDict(extra="forbid") + + quantity_kind: QuantityKind = Field(description="What kind of quantity was reported.") + assay_modality: AssayModality = Field(description="How it was measured.") + numerator_species: List[str] = Field( + min_length=1, description="Model species composing the numerator." + ) + denominator_species: List[str] = Field( + default_factory=list, + description="Model species composing the denominator. Empty for an absolute quantity.", + ) + experimental_denominator: Optional[str] = Field( + default=None, + description="What the experiment divided by, in the paper's own words. Required for a " + "density or a fraction; ``denominator_species`` names the model side.\n\n" + "Examples:\n" + "- 'mm^2 of tumor tissue (whole section including stroma)'\n" + "- 'all cells in ROI (all nucleated cells)'\n" + "- 'CD3+ T cells (pan-T-cell marker)'", + ) + reference: Optional[ReadoutReference] = Field( + default=None, + description="What a relative quantity is measured against. Required when quantity_kind " + "is 'foldchange', forbidden otherwise.", + ) + + @model_validator(mode="after") + def _composition_is_usable(self) -> "Readout": + if set(self.numerator_species) == set(self.denominator_species): + raise ValueError( + "numerator_species and denominator_species are identical, which is constant at 1." + ) + return self + + @model_validator(mode="after") + def _reference_matches_quantity_kind(self) -> "Readout": + needs = self.quantity_kind in REQUIRES_REFERENCE + if needs and self.reference is None: + raise ValueError( + f"quantity_kind='{self.quantity_kind.value}' is a relative quantity but no " + "reference is declared. Say what it is measured against: a timepoint in this " + "trajectory, or another scenario." + ) + if not needs and self.reference is not None: + raise ValueError( + f"quantity_kind='{self.quantity_kind.value}' is an absolute quantity but a " + "reference is declared. Only relative quantities take one." + ) + return self diff --git a/src/maple/core/calibration/registry_audit.py b/src/maple/core/calibration/registry_audit.py new file mode 100644 index 0000000..e449b0c --- /dev/null +++ b/src/maple/core/calibration/registry_audit.py @@ -0,0 +1,467 @@ +"""Cross-target checks against the cohort registry. + +Pydantic validators see one target at a time. These need the whole loaded set: +whether a cohort_id resolves, whether a target pools several sources, whether two +targets give one cohort the same quantity twice, and whether a cohort's rows are +deterministic functions of each other. +""" + +from __future__ import annotations + +import warnings +from dataclasses import dataclass +from typing import Any, Dict, FrozenSet, List, Optional, Tuple + +from maple.core.calibration.cohort import CohortRegistry + + +@dataclass(frozen=True) +class RegistryProblem: + """One cross-target defect.""" + + kind: str + target_ids: Tuple[str, ...] + detail: str + + +def _observable(target: Dict[str, Any]) -> Dict[str, Any]: + return target.get("observable") or {} + + +def _estimates(target: Dict[str, Any]) -> Dict[str, Any]: + return target.get("empirical_data") or {} + + +def _source_refs(target: Dict[str, Any]) -> List[str]: + return sorted( + {i.get("source_ref") for i in _estimates(target).get("inputs") or [] if i.get("source_ref")} + ) + + +def _is_literature(target: Dict[str, Any]) -> bool: + return (target.get("epistemic_basis") or "literature") == "literature" + + +def _arms(target: Dict[str, Any]) -> List[Dict[str, Any]]: + """Per-role inputs of a cross-scenario target; empty for a scalar target.""" + return _observable(target).get("inputs") or [] + + +def _readout(target: Dict[str, Any]) -> Dict[str, Any]: + return _observable(target).get("readout") or {} + + +def _composition(readout: Dict[str, Any]) -> Optional[Tuple[tuple, tuple]]: + """Declared (numerator, denominator) species, or None when nothing is declared.""" + num = tuple(sorted(readout.get("numerator_species") or [])) + den = tuple(sorted(readout.get("denominator_species") or [])) + return (num, den) if num else None + + +def _cohort_ids(target: Dict[str, Any]) -> List[str]: + """Cohorts this target draws on: one for a scalar target, one per arm for a contrast.""" + cid = target.get("cohort_id") + if cid: + return [cid] + return [a["cohort_id"] for a in _arms(target) if a.get("cohort_id")] + + +def find_registry_problems( + targets: Dict[str, Dict[str, Any]], cohorts: CohortRegistry +) -> List[RegistryProblem]: + """Every resolvable defect in ``{target_id: parsed_yaml}`` against the registry.""" + problems: List[RegistryProblem] = [] + by_cohort = cohorts.as_dict() + + for tid, data in sorted(targets.items()): + cid = data.get("cohort_id") + + if cid and cid not in by_cohort: + problems.append( + RegistryProblem( + "unknown_cohort", (tid,), f"cohort_id '{cid}' is not in the cohort registry." + ) + ) + + if not _is_literature(data): + continue + + refs = _source_refs(data) + if len(refs) > 1 and cid: + problems.append( + RegistryProblem( + "pooled_target", + (tid,), + f"inputs cite {len(refs)} sources ({refs}) but the target names cohort " + f"'{cid}'. A pooled estimate is not a cohort: its patients were never " + "measured together, so there is no resampling distribution over them. " + "Split into one target per source.", + ) + ) + + cohort = by_cohort.get(cid) if cid else None + if cohort is None: + continue + + n_eval = _estimates(data).get("n_evaluable") + if n_eval is not None and n_eval > cohort.n_c: + problems.append( + RegistryProblem( + "n_evaluable_exceeds_cohort", + (tid,), + f"n_evaluable={n_eval} exceeds cohort '{cid}' n_c={cohort.n_c}.", + ) + ) + src = (data.get("primary_data_source") or {}).get("source_tag") + if src and src != cohort.source_tag: + problems.append( + RegistryProblem( + "source_disagrees_with_cohort", + (tid,), + f"primary_data_source '{src}' differs from cohort '{cid}' source_tag " + f"'{cohort.source_tag}'.", + ) + ) + + problems.extend(_cross_scenario_arms(targets, cohorts)) + problems.extend(_duplicate_rows(targets)) + problems.extend(_singular_blocks(targets)) + return problems + + +def _cross_scenario_arms( + targets: Dict[str, Dict[str, Any]], cohorts: CohortRegistry +) -> List[RegistryProblem]: + """Each arm of a contrast against the cohort it names.""" + by_cohort = cohorts.as_dict() + problems: List[RegistryProblem] = [] + + for tid, data in sorted(targets.items()): + arms = _arms(data) + if not arms: + continue + for arm in arms: + cid = arm.get("cohort_id") + if not cid: + continue + cohort = by_cohort.get(cid) + if cohort is None: + problems.append( + RegistryProblem( + "unknown_cohort", + (tid,), + f"role '{arm.get('role')}' names cohort '{cid}', which is not in the " + "cohort registry.", + ) + ) + continue + if arm.get("scenario") and arm["scenario"] not in cohort.scenarios: + problems.append( + RegistryProblem( + "scenario_not_in_cohort", + (tid,), + f"role '{arm.get('role')}' runs scenario '{arm['scenario']}' but " + f"cohort '{cid}' declares {cohort.scenarios}. The arm's patients were " + "not measured under that condition.", + ) + ) + n_eval = arm.get("n_evaluable") + if n_eval is not None and n_eval > cohort.n_c: + problems.append( + RegistryProblem( + "n_evaluable_exceeds_cohort", + (tid,), + f"role '{arm.get('role')}' has n_evaluable={n_eval}, above cohort " + f"'{cid}' n_c={cohort.n_c}.", + ) + ) + + named = [c for c in _cohort_ids(data) if c in by_cohort] + overlapping = sorted( + { + tuple(sorted((a, b))) + for block in cohorts.blocks + for a in named + for b in named + if a < b and {a, b} <= block.members and block.overlap(a, b) != 0 + } + ) + if overlapping: + problems.append( + RegistryProblem( + "paired_contrast_as_cross_scenario", + (tid,), + f"arms name cohorts that declare shared patients: {overlapping}. The " + "contrast is then within-patient, and a paired contrast belongs in a " + "single CalibrationTarget whose observable declares a reference. Arms of " + "a cross-scenario target must be disjoint sets of people, since its " + "variance comes from resampling each arm independently.", + ) + ) + + problems.extend(_redundant_arms(targets)) + return problems + + +def _arm_composition(arm: Dict[str, Any]) -> Optional[Tuple[tuple, tuple]]: + """An arm's declared composition, in the same shape as a scalar target's.""" + return _composition(arm.get("readout") or {}) + + +def _redundant_arms(targets: Dict[str, Dict[str, Any]]) -> List[RegistryProblem]: + """An arm duplicating a standalone target on its cohort double-counts that number. + + A cross-scenario term earns a likelihood only when its per-arm constituents are + deliberately kept out of the fit. Once a constituent is also a target, the fit + conditions on it and on the contrast, which is the same belief twice. + """ + scalar: Dict[Tuple[str, tuple], List[str]] = {} + for tid, data in targets.items(): + cid = data.get("cohort_id") + if not cid or _arms(data): + continue + key = _composition(_readout(data)) + if key: + scalar.setdefault((cid, key), []).append(tid) + + problems = [] + for tid, data in sorted(targets.items()): + for arm in _arms(data): + cid = arm.get("cohort_id") + key = _arm_composition(arm) + if not cid or not key: + continue + for other in sorted(scalar.get((cid, key), [])): + problems.append( + RegistryProblem( + "redundant_cross_scenario_arm", + tuple(sorted((tid, other))), + f"role '{arm.get('role')}' computes what target '{other}' already " + f"reports for cohort '{cid}'. Conditioning on the constituent and on " + "the contrast counts one measurement twice; drop one.", + ) + ) + return problems + + +def _duplicate_rows(targets: Dict[str, Dict[str, Any]]) -> List[RegistryProblem]: + """Two targets computing one model quantity for one cohort are one row reported twice. + + Keyed on the readout's declared composition, so it covers absolute quantities + as well as ratios. + """ + seen: Dict[Tuple[str, tuple, tuple, Any], List[str]] = {} + for tid, data in targets.items(): + cid = data.get("cohort_id") + comp = _composition(_readout(data)) + if not cid or comp is None: + continue + key = (cid, comp[0], comp[1], _observable(data).get("readout_time")) + seen.setdefault(key, []).append(tid) + return [ + RegistryProblem( + "duplicate_row", + tuple(sorted(members)), + f"cohort '{key[0]}' has {len(members)} targets computing {list(key[1])} / " + f"{list(key[2])} at t={key[3]}. One cohort reports a quantity once; a second " + "target is either a duplicate or belongs to another cohort.", + ) + for key, members in sorted(seen.items()) + if len(members) > 1 + ] + + +#: Relative agreement at which two reported fractions count as complementary. +_SUM_TOL = 1e-3 + + +def _reported_center(target: Dict[str, Any]) -> Optional[float]: + """The location the source printed: a reported median, else a mean, else ``median``.""" + od = _estimates(target).get("observed_distribution") or {} + stats = od.get("statistics") or [] + for want, p in (("quantile", 0.5), ("mean", None), ("geometric_mean", None)): + for s in stats: + if s.get("stat") == want and s.get("p") == p: + return float(s["value"]) + med = _estimates(target).get("median") + return float(med[0]) if med else None + + +def _singular_blocks(targets: Dict[str, Dict[str, Any]]) -> List[RegistryProblem]: + """Rows of one cohort that are deterministic functions of each other. + + Such a block is singular, and the ridge that keeps it invertible turns into a + huge direction in its inverse rather than an error. Two cheap detectors: a set + of fractions whose numerators partition their shared denominator, and a pair + whose printed centers sum to one. + """ + fractions: Dict[str, List[Tuple[str, frozenset, frozenset, Optional[float]]]] = {} + for tid, data in targets.items(): + cid = data.get("cohort_id") + readout = _readout(data) + comp = _composition(readout) + if not cid or comp is None or readout.get("quantity_kind") != "fraction": + continue + fractions.setdefault(cid, []).append( + (tid, frozenset(comp[0]), frozenset(comp[1]), _reported_center(data)) + ) + + problems: List[RegistryProblem] = [] + for cid, rows in sorted(fractions.items()): + problems.extend(_partitioning_rows(cid, rows)) + problems.extend(_centers_summing_to_one(cid, rows)) + return problems + + +def _partitioning_rows( + cid: str, rows: List[Tuple[str, frozenset, frozenset, Any]] +) -> List[RegistryProblem]: + """Fractions over one denominator whose numerators exactly partition it sum to 1.""" + by_den: Dict[frozenset, List[Tuple[str, frozenset]]] = {} + for tid, num, den, _ in rows: + by_den.setdefault(den, []).append((tid, num)) + + out: List[RegistryProblem] = [] + for den, members in sorted(by_den.items(), key=lambda kv: sorted(t for t, _ in kv[1])): + if len(members) < 2: + continue + union: set = set() + disjoint = True + for _, num in members: + if union & num: + disjoint = False + break + union |= num + if not disjoint or union != set(den): + continue + out.append( + RegistryProblem( + "singular_block", + tuple(sorted(tid for tid, _ in members)), + f"cohort '{cid}' has {len(members)} fractions over {sorted(den)} whose " + f"numerators partition it, so they sum to 1 by construction and the block " + "is singular. Drop one, or widen the denominator so they do not exhaust it.", + ) + ) + return out + + +def _centers_summing_to_one( + cid: str, rows: List[Tuple[str, frozenset, frozenset, Any]] +) -> List[RegistryProblem]: + """Two printed fractions summing to one are one number reported twice.""" + out: List[RegistryProblem] = [] + centers = sorted((tid, c) for tid, _, _, c in rows if c is not None) + for i, (tid_a, a) in enumerate(centers): + for tid_b, b in centers[i + 1 :]: + total = 1.0 if max(a, b) <= 1.0 else 100.0 + if abs(a + b - total) <= _SUM_TOL * total: + out.append( + RegistryProblem( + "singular_block", + tuple(sorted((tid_a, tid_b))), + f"cohort '{cid}' reports {a} and {b}, which sum to {total}. As data " + "these are one number reported twice and the block is singular, even " + "where the model observables are not exactly complementary.", + ) + ) + return out + + +def check_registries(targets: Dict[str, Dict[str, Any]], cohorts: CohortRegistry) -> None: + """Raise on any registry defect; warn on unused cohorts and on merged blocks.""" + problems = find_registry_problems(targets, cohorts) + warn_unused_cohorts(targets, cohorts) + warn_merged_blocks(targets, cohorts) + if not problems: + return + lines = [f"{len(problems)} registry problem(s):"] + for p in problems: + lines.append(f"\n[{p.kind}] {', '.join(p.target_ids)}\n {p.detail}") + raise ValueError("\n".join(lines)) + + +def warn_unused_cohorts(targets: Dict[str, Dict[str, Any]], cohorts: CohortRegistry) -> List[str]: + """Warn on cohorts no target refers to. Returns the unused ids.""" + used = {d.get("cohort_id") for d in targets.values() if d.get("cohort_id")} + unused = sorted(c.cohort_id for c in cohorts.cohorts if c.cohort_id not in used) + if unused: + warnings.warn( + f"Cohorts no target uses: {unused}. Stale entries drift out of step with the " + "corpus; remove them or add the target.", + UserWarning, + ) + return unused + + +def covariance_blocks( + targets: Dict[str, Dict[str, Any]], cohorts: CohortRegistry +) -> List[FrozenSet[str]]: + """Cohorts that must share one covariance block, as connected components. + + Blocks are assumed independent, so two cohorts must sit in one whenever a row + depends on both. Two ways that happens, and they need opposite draws. + + A target drawing on both joins disjoint sets of people: resample each cohort + independently, which puts the zeros between arms and the covariance around the + derived row. Cohorts that share patients are resampled once from the block's + patient set, each cohort's rows evaluated on its own members; drawing those + independently would restore the independence the block exists to deny. + + The return keeps the partition and not which relation joined a pair, nor + whether the sharing was counted. A caller building V reads the registry's + blocks for both: a counted block supplies the strata a joint draw needs, an + uncounted one only says that drawing its cohorts apart is wrong. + """ + parent = {c.cohort_id: c.cohort_id for c in cohorts.cohorts} + + def find(x: str) -> str: + while parent[x] != x: + parent[x] = parent[parent[x]] + x = parent[x] + return x + + def union(a: str, b: str) -> None: + ra, rb = find(a), find(b) + if ra != rb: + parent[max(ra, rb)] = min(ra, rb) + + # Counted and uncounted alike: both say the cohorts are not independent, and + # they differ in whether a joint resample is defined, not in the partition. + for block in cohorts.blocks: + members = [c for c in block.cohorts if c in parent] + for other in members[1:]: + union(members[0], other) + + for data in targets.values(): + named = [c for c in _cohort_ids(data) if c in parent] + for other in named[1:]: + union(named[0], other) + + blocks: Dict[str, set] = {} + for cid in parent: + blocks.setdefault(find(cid), set()).add(cid) + return sorted((frozenset(b) for b in blocks.values()), key=lambda b: sorted(b)) + + +def warn_merged_blocks( + targets: Dict[str, Dict[str, Any]], cohorts: CohortRegistry +) -> List[FrozenSet[str]]: + """Warn on blocks spanning several cohorts. Returns them.""" + merged = [b for b in covariance_blocks(targets, cohorts) if len(b) > 1] + for block in merged: + warnings.warn( + f"Cohorts {sorted(block)} form one covariance block: a row depends on more than " + "one of them, so they are not independent. Inference must draw them together.", + UserWarning, + ) + return merged + + +def resolve_n(target: Dict[str, Any], cohorts: CohortRegistry) -> Optional[int]: + """Patients behind this target's statistics: ``n_evaluable``, else the cohort's ``n_c``.""" + n_eval = _estimates(target).get("n_evaluable") + if n_eval is not None: + return int(n_eval) + cohort = cohorts.get(target.get("cohort_id") or "") + return cohort.n_c if cohort else None diff --git a/src/maple/core/calibration/shared_models.py b/src/maple/core/calibration/shared_models.py index f729df7..575e0c3 100644 --- a/src/maple/core/calibration/shared_models.py +++ b/src/maple/core/calibration/shared_models.py @@ -8,7 +8,7 @@ import math from enum import Enum -from typing import List, Literal, Optional, Union +from typing import List, Optional, Union from pydantic import BaseModel, ConfigDict, Field, field_validator, model_validator @@ -93,10 +93,8 @@ class SpreadSource(str, Enum): """Provenance of a reported spread — what kind of variability it measures. This drives whether a spread feeds the population-spread hyperparameter - (omega) in hierarchical inference or only the center's error budget. - Subsumes the calibration side's earlier two-valued ``population_spread`` - (``across_patient`` / ``center_only``) and adds the finer provenance the - submodel (in-vitro / ex-vivo) targets need. + (omega) in hierarchical inference or only the center's error budget. It is + the sole declaration of that, on both calibration and submodel targets. """ ACROSS_PATIENT = "across_patient" @@ -147,174 +145,118 @@ class ExperimentalUnitType(str, Enum): CLONAL = "clonal" # clones / passages of one line — treat as technical for spread -class QuantileAnchor(BaseModel): - """One (probability, value) point on an observed distribution's quantile function.""" +class StatKind(str, Enum): + """A statistic a source can print. What it is, not what a model does with it.""" - model_config = ConfigDict(extra="forbid") - - p: float = Field( - description="Probability level in the open interval (0, 1). " - "0.5 is the median; 0.25/0.75 are the IQR edges." - ) - value: float = Field(description="Observed value at this quantile (in ``units``).") - - @field_validator("p") - @classmethod - def _p_in_open_unit_interval(cls, v: float) -> float: - if not (0.0 < v < 1.0): - raise ValueError(f"quantile probability p must be in (0, 1), got {v}") - return v - - -# Standard normal quantiles used to expand a scalar scale into quartile anchors. -_Z_Q = 0.6744897501960817 # Phi^-1(0.75): the 0.25/0.75 quantile of N(0,1) -_Z_95 = 1.959963984540054 # Phi^-1(0.975): half-width of a 95% normal interval - - -class DistributionShape(str, Enum): - """Shape used to expand a reported center+scale into quantile anchors.""" - - NORMAL = "normal" - LOGNORMAL = "lognormal" - LOGIT_NORMAL = "logit_normal" # bounded to (0, 1): fractions, probabilities + QUANTILE = "quantile" # needs `p`; the median is p=0.5 + MEAN = "mean" + GEOMETRIC_MEAN = "geometric_mean" + SD = "sd" # dispersion of the sample + CV = "cv" # sd / mean, dimensionless + IQR = "iqr" # width, q75 - q25 + RANGE = "range" # width, max - min + SE = "se" # dispersion of the estimator, not the sample + CI95_LO = "ci95_lo" + CI95_HI = "ci95_hi" + MIN = "min" + MAX = "max" -class ScaleType(str, Enum): - """What kind of scale a reported dispersion value is. +class QuantileConvention(str, Enum): + """Which order statistic a package returns for a quantile, by Hyndman-Fan type. - Determines how it expands to quartiles and whether it needs ``n_biological``. + At small n the choice moves the answer: asking for a lower quartile of 9 points + gives the population's 0.30 quantile under ``type7`` and its 0.25 under + ``type6``, some 0.4 sampling standard errors apart. Papers do not state it, so + this is usually unknown and left unset. """ - SD = "sd" # population standard deviation (linear units) - SEM = "sem" # standard error of the mean: SD = SEM * sqrt(n_biological) - CV = "cv" # coefficient of variation, SD/mean (dimensionless) - IQR = "iqr" # full interquartile range, q75 - q25 (linear units) - CI95_HALFWIDTH = "ci95_halfwidth" # half-width of a 95% interval (linear units) + TYPE2 = "type2" # averaged inverted cdf; SAS + TYPE4 = "type4" # interpolated inverted cdf + TYPE6 = "type6" # SPSS, Minitab, Excel PERCENTILE.EXC + TYPE7 = "type7" # R, numpy, pandas, Excel PERCENTILE + TYPE8 = "type8" # median-unbiased; Hyndman and Fan's recommendation -class MomentSpread(BaseModel): - """A reported distribution given as center + scale + shape (mean +/- SD, etc.). +#: Statistics that locate the distribution. +LOCATION_STATS: frozenset = frozenset({StatKind.MEAN, StatKind.GEOMETRIC_MEAN}) +#: Statistics that describe a width, in the value's own units unless noted. +WIDTH_STATS: frozenset = frozenset({StatKind.SD, StatKind.CV, StatKind.IQR, StatKind.RANGE}) +#: Widths of an estimator rather than of the sample. These carry a population +#: width too, but only alongside the n they were computed over: sd = se * sqrt(n). +#: Kept apart from WIDTH_STATS so nothing reads one as the other by accident. +SAMPLING_WIDTH_STATS: frozenset = frozenset({StatKind.SE}) - This is the form most papers actually report. It is an alternative to explicit - ``quantiles`` on :class:`ObservedDistribution`: the framework expands it to - quartile anchors once, centrally, so extractors never hand-convert mean +/- SD - into q25/q50/q75 (an error-prone step that also duplicates snippet-validated - inputs). The imposed ``shape`` is recorded, so nothing is silent. - """ + +class ReportedStatistic(BaseModel): + """One number the source printed.""" model_config = ConfigDict(extra="forbid") - center: float = Field(description="Reported central value (see ``center_type``).") - center_type: Literal["mean", "median"] = Field( - default="mean", description="Whether ``center`` is the arithmetic mean or the median." + stat: StatKind = Field(description="Which statistic this is.") + value: float = Field( + description="The printed value, in the target's units (dimensionless for cv)." ) - scale: float = Field( - description="Reported dispersion value, in the same units as ``center`` " - "(dimensionless for ``cv``). Interpreted per ``scale_type``.", - ge=0.0, - ) - scale_type: ScaleType = Field(description="What kind of scale ``scale`` is (see ScaleType).") - shape: DistributionShape = Field( - description="Shape used to expand center+scale into quartiles. Records the " - "imposed shape explicitly." + p: Optional[float] = Field( + default=None, + description="Probability level in (0, 1). Required for stat='quantile', forbidden " + "otherwise. The median is p=0.5.", ) - def to_quartiles(self, n_biological: Optional[int] = None) -> tuple: - """Expand to (q25, q50, q75). ``n_biological`` is required for scale_type='sem'.""" - return _expand_moments(self, n_biological) + @model_validator(mode="after") + def _p_matches_stat(self) -> "ReportedStatistic": + if self.stat == StatKind.QUANTILE: + if self.p is None: + raise ValueError("stat='quantile' needs a probability level p.") + if not (0.0 < self.p < 1.0): + raise ValueError(f"quantile p must be in (0, 1), got {self.p}.") + elif self.p is not None: + raise ValueError(f"stat='{self.stat.value}' must not set p; only quantiles have one.") + if self.stat in (WIDTH_STATS | SAMPLING_WIDTH_STATS) and self.value < 0: + raise ValueError(f"stat='{self.stat.value}' is a width and cannot be negative.") + return self -def _expand_moments(m: "MomentSpread", n_biological: Optional[int]) -> tuple: - """Expand a center+scale+shape spec into (q25, q50, q75). +# Standard normal quantiles used to expand a scalar scale into quartiles. +_Z_Q = 0.6744897501960817 # Phi^-1(0.75) +_Z_95 = 1.959963984540054 # Phi^-1(0.975) - Handles the common, unambiguous cases; raises with a clear pointer for combos - that a single center+scale cannot determine (e.g. lognormal + bare IQR). - """ - center = m.center - def _linear_sd() -> float: - if m.scale_type == ScaleType.SD: - return m.scale - if m.scale_type == ScaleType.SEM: - if n_biological is None: - raise ValueError( - "scale_type='sem' needs n_biological to recover the population SD " - "(SD = SEM * sqrt(n))." - ) - return m.scale * math.sqrt(n_biological) - if m.scale_type == ScaleType.CI95_HALFWIDTH: - return m.scale / _Z_95 - if m.scale_type == ScaleType.CV: - return m.scale * abs(center) - raise ValueError(f"scale_type '{m.scale_type}' has no linear-SD form") # pragma: no cover - - if m.shape == DistributionShape.NORMAL: - # mean and median coincide. - if m.scale_type == ScaleType.IQR: - half = m.scale / 2.0 - return (center - half, center, center + half) - sd = _linear_sd() +class DistributionShape(str, Enum): + """Shape used to expand a center + scale into quantiles, when asked.""" + + NORMAL = "normal" + LOGNORMAL = "lognormal" + LOGIT_NORMAL = "logit_normal" # bounded to (0, 1): fractions, probabilities + + +def _quartiles_from_center_scale( + center: float, center_is_mean: bool, sd: float, shape: DistributionShape +) -> tuple: + """Expand a center and a normal-equivalent SD into (q25, q50, q75).""" + if shape == DistributionShape.NORMAL: return (center - _Z_Q * sd, center, center + _Z_Q * sd) - if m.shape == DistributionShape.LOGIT_NORMAL: - # Bounded to (0, 1): expand in logit space so quartiles never escape the - # bounds. logit(X) ~ Normal(mu_l, sigma_l); the reported scale sets sigma_l - # via the delta method at the median (d logit/dx = 1/(x(1-x))). + if shape == DistributionShape.LOGIT_NORMAL: if not (0.0 < center < 1.0): raise ValueError( - "shape='logit_normal' needs center in the open interval (0, 1) " - f"(it is a bounded fraction/probability), got {center}." - ) - if m.center_type != "median": - raise ValueError( - "shape='logit_normal' needs center_type='median' (the logit map is " - "applied at the median). Provide the median, or use explicit quantiles." + f"shape='logit_normal' needs a center in (0, 1), got {center}. Express a " + "percent as a fraction." ) - median = center - mu_l = math.log(median / (1.0 - median)) - if m.scale_type == ScaleType.IQR: - # linear IQR -> linear SD (normal-equivalent) -> logit-space sigma - sd_linear = m.scale / (2.0 * _Z_Q) - else: - sd_linear = _linear_sd() # sd/sem/ci95 -> linear SD; cv -> cv*|center| - sigma_l = sd_linear / (median * (1.0 - median)) + mu_l = math.log(center / (1.0 - center)) + sigma_l = sd / (center * (1.0 - center)) def _expit(z: float) -> float: return 1.0 / (1.0 + math.exp(-z)) - return (_expit(mu_l - _Z_Q * sigma_l), median, _expit(mu_l + _Z_Q * sigma_l)) + return (_expit(mu_l - _Z_Q * sigma_l), center, _expit(mu_l + _Z_Q * sigma_l)) # lognormal - if m.scale_type == ScaleType.IQR: - # median (IQR) is the common clinical form. Given the median and the IQR, - # solve for sigma_ln: IQR = median * 2 * sinh(Z_Q * sigma_ln). - if m.center_type != "median": - raise ValueError( - "shape='lognormal' with scale_type='iqr' needs center_type='median' " - "(median + IQR is well-determined; mean + IQR is not). Provide the median, " - "or use explicit quantiles." - ) - median = center - sigma_ln = math.asinh(m.scale / (2.0 * median)) / _Z_Q - return ( - median * math.exp(-_Z_Q * sigma_ln), - median, - median * math.exp(_Z_Q * sigma_ln), - ) - if m.scale_type == ScaleType.CV: - sigma_ln = math.sqrt(math.log(1.0 + m.scale**2)) - median = center if m.center_type == "median" else center / math.sqrt(1.0 + m.scale**2) - else: # sd / sem / ci95_halfwidth -> a linear SD, which needs the mean to form CV - if m.center_type != "mean": - raise ValueError( - "shape='lognormal' with a linear scale (sd/sem/ci95_halfwidth) needs " - "center_type='mean' to form CV=SD/mean. Use scale_type='cv' with a median " - "center, or provide explicit quantiles." - ) - sd = _linear_sd() - cv = sd / abs(center) - sigma_ln = math.sqrt(math.log(1.0 + cv**2)) - median = center / math.sqrt(1.0 + cv**2) + if center <= 0: + raise ValueError(f"shape='lognormal' needs a positive center, got {center}.") + cv = sd / abs(center) + sigma_ln = math.sqrt(math.log(1.0 + cv**2)) + median = center / math.sqrt(1.0 + cv**2) if center_is_mean else center return ( median * math.exp(-_Z_Q * sigma_ln), median, @@ -323,212 +265,220 @@ def _expit(z: float) -> float: class ObservedDistribution(BaseModel): - """General representation of a reported distribution. - - SD, SEM, CV, IQR, CI95, quartiles, deciles, and full samples are all partial - specifications of one distribution. This object is the unifying data layer, - authored in whichever form the paper reports and reduced to quartiles centrally: - - - ``moments``: center + scale + shape (mean +/- SD, median +/- IQR, CV, CI). The - dominant form in the literature; the framework expands it to quartiles so the - extractor never hand-converts (and never restates snippet-validated inputs). - - ``quantiles``: explicit quantile anchors, for sources that give - quartiles/percentiles/samples directly. - - Provide EXACTLY ONE. Derivations (``median``/``iqr``/``quantile``) work off - whichever form is present. This is orthogonal to ``spread_source``: the form - carries the *shape*, the ``spread_source`` tag carries the *provenance* (whether - that shape is genuine population spread). Both route the quantity in inference. + """What a source printed about a distribution, as a flat list of statistics. + + One entry is one number the paper reported. ``spread_source`` says whether the + width among them is genuine population variability. Nothing here says what a + model should do with any entry; that is the consumer's decision. + + ``quantile`` / ``median`` / ``iqr`` derive values on demand, using explicit + quantile entries where available and otherwise expanding a center and scale + through ``shape``. """ model_config = ConfigDict(extra="forbid") - quantiles: Optional[List[QuantileAnchor]] = Field( - default=None, - description="Quantile anchors: the median (p=0.5) plus whatever scale anchors the " - "source reports (IQR edges at minimum for a spread). Use this when the paper gives " - "quantiles/percentiles/samples directly. Provide EITHER quantiles OR moments.", + statistics: List[ReportedStatistic] = Field( + description="Every statistic the source printed for this quantity. Record what was " + "actually on the page: median plus quartiles, mean plus SD, an SE, a range. Do not " + "convert between them." + ) + spread_source: SpreadSource = Field( + description="Provenance of the spread these statistics describe (see SpreadSource)." ) - moments: Optional["MomentSpread"] = Field( + shape: Optional[DistributionShape] = Field( default=None, - description="Center + scale + shape form (mean +/- SD, median +/- IQR, CV, CI). Use " - "this when the paper reports moments rather than quantiles — the framework expands it " - "to quartiles centrally, so you never hand-convert. Provide EITHER quantiles OR moments.", + description="Shape assumed when deriving quantiles from a center and a scale. Only " + "needed when the source printed no quartiles and a consumer asks for them; recorded " + "here so the assumption is explicit rather than applied silently.", ) - spread_source: SpreadSource = Field( - description="Provenance of the spread these anchors describe (see SpreadSource). " - "Determines whether the spread feeds population omega or only the center." + quantile_convention: Optional[QuantileConvention] = Field( + default=None, + description="Which order statistic the source's software returned for a quantile " + "(see QuantileConvention). Set it only where the paper or its methods say so, which " + "is rare. Unset means unrecorded, and a consumer that needs one has to pick a " + "default and report what the choice costs.", ) n_biological: Optional[int] = Field( default=None, ge=1, - description="Number of BIOLOGICAL units (donors/animals/patients) the summary " - "is computed over. REQUIRED when spread_source declares a population spread " - "(across_patient / biological_experimental): it licenses the SD<->SEM round-trip " - "and sets per-target finite-sample noise. Distinct from technical replicates.", + description="Biological units (donors/animals/patients) the summary is computed over. " + "For submodel targets, where there is no cohort registry to carry it. Calibration " + "targets name a cohort instead and must leave this unset.", ) n_biological_is_floor: bool = Field( default=False, - description="True when ``n_biological`` is a LOWER BOUND, not an exact count — the " - "source reports the unit count as 'n>=8', 'at least 8 donors', 'n=8-12 across " - "conditions', etc. Consumers that weight panels by precision (finite-sample noise " - "~ 1/sqrt(n), inverse-variance moment weighting) must treat a floor conservatively: " - "an exact-looking n from a floor over-states precision and over-weights the panel. " - "Leave False only when the source gives an exact per-summary n.", + description="True when ``n_biological`` is a lower bound rather than an exact count.", ) n_technical: Optional[int] = Field( - default=None, - ge=1, - description="Number of technical replicates, when reported separately. Does not " - "license the SD = SEM*sqrt(n) recovery.", + default=None, ge=1, description="Technical replicates, when reported separately." ) experimental_unit_type: Optional[ExperimentalUnitType] = Field( default=None, - description="What one replicate is (biological/technical/clonal). Gates whether " - "an SEM can be converted to a population SD. REQUIRED when spread_source declares " - "a population spread (and must be 'biological' — technical/clonal spreads are not " - "population variability). Omit only for center_only / technical sources.", - ) - shape_assumption: Optional[str] = Field( - default=None, - description="Distributional shape imposed when anchors were expanded from a SCALAR " - "scale with no shape information (e.g. 'lognormal', 'normal'). Records the imposed " - "shape so it is explicit rather than silent. Omit when anchors come directly from " - "reported quantiles/samples.", + description="What one replicate is. Gates whether an SE can be read as a population " + "SD. Submodel side; calibration targets carry it on the cohort.", ) unit_group: Optional[str] = Field( default=None, - description="Name of the shared biological-unit group this observation belongs to. " - "DEFAULT (omitted) = this observable is its own group — the correct choice for the " - "vast majority of single-observable targets. Set a shared string ONLY across " - "observables measured on the SAME biological units (e.g. the same donor panel across " - "the doses of a dose-response, or one cohort followed over a time course). Grouped " - "observables are moment-matched JOINTLY for the population spread (omega), so a " - "donor's variation is treated as one shared random effect rather than K independent " - "measurements (which would spuriously shrink the spread by ~sqrt(K)). Do NOT share a " - "group across observables from DIFFERENT populations (different mouse lines, " - "genotypes, or proxy species) — those are distinct spreads and must stay separate. " - "A string (not a bool) so a target can carry two independent groups if it has two " - "unit sets. Members of one group must agree on n_biological / spread_source / " - "experimental_unit_type (they describe the shared batch).", - ) - - @field_validator("quantiles") - @classmethod - def _non_empty(cls, v: Optional[List[QuantileAnchor]]) -> Optional[List[QuantileAnchor]]: - if v is not None and len(v) == 0: - raise ValueError("observed_distribution.quantiles must have at least one anchor") - return v - - @model_validator(mode="after") - def _exactly_one_form(self) -> "ObservedDistribution": - """Exactly one of quantiles / moments must be given.""" - has_q = self.quantiles is not None - has_m = self.moments is not None - if has_q == has_m: - raise ValueError( - "observed_distribution must specify EXACTLY ONE of 'quantiles' or 'moments' " - f"(got quantiles={'set' if has_q else 'unset'}, " - f"moments={'set' if has_m else 'unset'})." - ) - return self + description="Shared biological-unit group, for observables measured on the SAME units " + "(one donor panel across the doses of a dose-response). Grouped observables share one " + "random effect rather than counting as independent measurements. Submodel side; " + "calibration targets name a cohort instead.", + ) @model_validator(mode="after") - def _validate_moments_derivable(self) -> "ObservedDistribution": - """Fail fast on a moments spec that a single center+scale cannot determine.""" - if self.moments is not None: - _expand_moments(self.moments, self.n_biological) # raises with a clear pointer - return self + def _statistics_well_formed(self) -> "ObservedDistribution": + if not self.statistics: + raise ValueError("observed_distribution.statistics must have at least one entry.") + + keys = [(s.stat, s.p) for s in self.statistics] + dupes = sorted( + { + f"{k[0].value}{'' if k[1] is None else f'@p={k[1]}'}" + for k in keys + if keys.count(k) > 1 + } + ) + if dupes: + raise ValueError(f"observed_distribution repeats statistic(s): {dupes}") - @model_validator(mode="after") - def _validate_quantile_function(self) -> "ObservedDistribution": - # Only applies to the explicit-quantiles form. - if self.quantiles is None: - return self - # Unique, sorted probability levels with non-decreasing values (a valid, - # non-crossing quantile function). - ps = [q.p for q in self.quantiles] - if len(set(ps)) != len(ps): - raise ValueError( - f"observed_distribution has duplicate probability levels: {sorted(ps)}" - ) - ordered = sorted(self.quantiles, key=lambda q: q.p) - prev = None - for q in ordered: - if prev is not None and q.value < prev.value: + quants = sorted( + ((s.p, s.value) for s in self.statistics if s.stat == StatKind.QUANTILE), + key=lambda t: t[0], + ) + for (p_lo, v_lo), (p_hi, v_hi) in zip(quants[:-1], quants[1:]): + if v_hi < v_lo: raise ValueError( - "observed_distribution quantiles must be non-decreasing in value with p: " - f"value {q.value} at p={q.p} is below value {prev.value} at p={prev.p}" + "quantiles must be non-decreasing in p: " + f"value {v_hi} at p={p_hi} is below value {v_lo} at p={p_lo}." ) - prev = q - # Keep anchors stored in probability order. - object.__setattr__(self, "quantiles", ordered) - # A quantiles-form spread that feeds omega needs at least two anchors (a scale, - # not just a center). The moments form always carries a scale. - if self.spread_source in POPULATION_SPREAD_SOURCES and len(self.quantiles) < 2: + if self.spread_source in POPULATION_SPREAD_SOURCES and not self._has_width(): raise ValueError( - f"spread_source='{self.spread_source.value}' declares a population spread " - "but only a single quantile anchor is given (no scale). Provide IQR edges " - "(p=0.25, 0.75) or set spread_source='center_only'." + f"spread_source='{self.spread_source.value}' declares a population spread but " + "no width statistic was reported. Provide quartiles, an sd/iqr/cv/range, or " + "declare a center-only spread_source." ) - return self - @model_validator(mode="after") - def _require_biological_provenance_for_spread(self) -> "ObservedDistribution": - """A population-spread claim must state its biological unit count and unit type. - - These are what separate a genuine cross-donor/patient spread from an SEM-scale - width over technical replicates — the exact conflation this schema exists to - prevent. Only enforced when spread_source feeds the population-spread magnitude; - center_only / technical / translation / assumed sources are exempt. Applies to - both the quantiles and moments forms. - """ - if self.spread_source not in POPULATION_SPREAD_SOURCES: - return self - # SD<->SEM recovery is only licensed for biological units. - if self.experimental_unit_type in ( - ExperimentalUnitType.TECHNICAL, - ExperimentalUnitType.CLONAL, + if self.quantile_convention is not None and not any( + s.stat == StatKind.QUANTILE for s in self.statistics ): raise ValueError( - f"experimental_unit_type='{self.experimental_unit_type.value}' cannot support " - f"spread_source='{self.spread_source.value}': a spread over " - "technical/clonal replicates is not population variability. Use " - "'center_only' (or 'technical'), or provide a biological n." + f"quantile_convention='{self.quantile_convention.value}' was set but no " + "quantile was reported; it describes how a quantile was computed." ) - if self.n_biological is None: + return self + + def _has_width(self) -> bool: + """Whether any entry carries a width, explicit or recoverable from one. + + An SE counts: it is the sample width divided by sqrt(n), so it determines + one given the n it was computed over. Recovering it needs that n, which + the derivation accessors take as an argument rather than assume. + """ + kinds = {s.stat for s in self.statistics} + if kinds & (WIDTH_STATS | SAMPLING_WIDTH_STATS): + return True + if len({s.p for s in self.statistics if s.stat == StatKind.QUANTILE}) >= 2: + return True + # A min/max pair is the sample range stated as its two endpoints, which is + # what a source prints when it gives an observed span rather than a width. + return {StatKind.MIN, StatKind.MAX} <= kinds + + # ---- Accessors -------------------------------------------------------- + + def get(self, stat: StatKind, p: Optional[float] = None) -> Optional[float]: + """The reported value for ``stat`` (at ``p`` for a quantile), or None.""" + for s in self.statistics: + if s.stat == stat and s.p == p: + return s.value + return None + + def center(self) -> Optional[float]: + """The reported center: the median if given, else the mean.""" + med = self.get(StatKind.QUANTILE, 0.5) + return med if med is not None else self.get(StatKind.MEAN) + + # ---- Derivations ------------------------------------------------------ + + def _anchor_pairs(self, n: Optional[int] = None) -> List[tuple]: + """Effective (p, value) quantile anchors, p-sorted. + + Explicit quantiles are used as reported. With fewer than two of them, a + center and a scale are expanded through ``shape``. ``n`` is the count the + statistics were computed over, needed only to widen an SE into a sample SD. + """ + quants = sorted( + ((s.p, s.value) for s in self.statistics if s.stat == StatKind.QUANTILE), + key=lambda t: t[0], + ) + if len(quants) >= 2: + return quants + + sd = self._normal_equivalent_sd(n) + center = self.center() + if sd is None or center is None: + missing = "no center" if center is None else "no width" + if center is not None and n is None and self.get(StatKind.SE) is not None: + missing = ( + "only an SE, which is a width per sqrt(n); pass n (the cohort's n_c) " + "to widen it into a sample SD" + ) raise ValueError( - f"spread_source='{self.spread_source.value}' declares a population spread " - "but n_biological is not set. A population-spread claim needs a biological " - "unit count (donors/animals/patients) to license the SD<->SEM round-trip " - "and set per-target finite-sample noise. Provide n_biological, or use " - "spread_source='center_only'/'technical'." + "observed_distribution cannot produce quantiles: it reports " + f"{[s.stat.value for s in self.statistics]}, which gives {missing}." ) - if self.experimental_unit_type is None: + if self.shape is None: raise ValueError( - f"spread_source='{self.spread_source.value}' declares a population spread " - "but experimental_unit_type is not set. State that the n counts biological " - "units (technical/clonal spreads are not population variability)." + "observed_distribution needs `shape` to expand a center and scale into " + "quantiles, since the source reported fewer than two quantiles." ) - return self - - # ---- Derivations (median / quantile / IQR / scale) -------------------- + center_is_mean = self.get(StatKind.QUANTILE, 0.5) is None + q25, q50, q75 = _quartiles_from_center_scale(center, center_is_mean, sd, self.shape) + return [(0.25, q25), (0.5, q50), (0.75, q75)] - def _anchor_pairs(self) -> List[tuple]: - """Effective (p, value) anchors from whichever form is present, p-sorted. + def _normal_equivalent_sd(self, n: Optional[int] = None) -> Optional[float]: + """A linear SD implied by whichever width statistic was reported. - The quantiles form returns its anchors; the moments form is expanded to - (q25, q50, q75) once, centrally. + An SE is used last and only with ``n``: it is the sample SD over sqrt(n), + so without n it is not a sample width at all. At n=10 the two differ by a + factor of 3.2, which is why this never guesses. """ - if self.quantiles is not None: - return [(q.p, q.value) for q in self.quantiles] # sorted by validator - q25, q50, q75 = _expand_moments(self.moments, self.n_biological) - return [(0.25, q25), (0.5, q50), (0.75, q75)] - - def quantile(self, p: float) -> float: - """Linearly interpolated value at probability level ``p`` (clamped to the anchor range).""" - anchors = self._anchor_pairs() + sd = self.get(StatKind.SD) + if sd is not None: + return sd + iqr = self.get(StatKind.IQR) + if iqr is not None: + return iqr / (2.0 * _Z_Q) + cv = self.get(StatKind.CV) + if cv is not None: + center = self.center() + return cv * abs(center) if center is not None else None + lo, hi = self.get(StatKind.CI95_LO), self.get(StatKind.CI95_HI) + if lo is not None and hi is not None: + return (hi - lo) / (2.0 * _Z_95) + se = self.get(StatKind.SE) + if se is not None and n is not None: + if n < 1: + raise ValueError(f"widening an SE needs a positive n, got {n}.") + return se * math.sqrt(n) + return None + + def population_sd(self, n: Optional[int] = None) -> Optional[float]: + """The sample SD the statistics imply, or None if none is recoverable. + + ``n`` is the count the statistics were computed over: a calibration + target's cohort ``n_c``, or a submodel target's ``n_biological``. It is + read only when the reported width is an SE. + """ + return self._normal_equivalent_sd(n) + + def quantile(self, p: float, n: Optional[int] = None) -> float: + """Value at ``p``: reported if the source printed it, else interpolated.""" + reported = self.get(StatKind.QUANTILE, p) + if reported is not None: + return reported + anchors = self._anchor_pairs(n) if p <= anchors[0][0]: return anchors[0][1] if p >= anchors[-1][0]: @@ -537,26 +487,56 @@ def quantile(self, p: float) -> float: if p_lo <= p <= p_hi: if p_hi == p_lo: return v_lo - frac = (p - p_lo) / (p_hi - p_lo) - return v_lo + frac * (v_hi - v_lo) + return v_lo + (p - p_lo) / (p_hi - p_lo) * (v_hi - v_lo) return anchors[-1][1] # pragma: no cover - def median(self) -> float: - """Value at p=0.5 (interpolated if not an explicit anchor).""" - return self.quantile(0.5) + def median(self, n: Optional[int] = None) -> float: + """Value at p=0.5, reported if given and interpolated otherwise.""" + return self.quantile(0.5, n) - def iqr(self) -> Optional[float]: - """Interquartile range (q75 - q25), or None if the anchors do not span it.""" - ps = [p for p, _ in self._anchor_pairs()] + def iqr(self, n: Optional[int] = None) -> Optional[float]: + """Interquartile range, reported if given, else from the quantile anchors.""" + reported = self.get(StatKind.IQR) + if reported is not None: + return reported + ps = [p for p, _ in self._anchor_pairs(n)] if min(ps) > 0.25 or max(ps) < 0.75: return None - return self.quantile(0.75) - self.quantile(0.25) + return self.quantile(0.75, n) - self.quantile(0.25, n) @property def feeds_population_spread(self) -> bool: - """Whether these anchors contribute a base population-spread magnitude (omega).""" + """Whether these statistics contribute a base population-spread magnitude (omega).""" return self.spread_source in POPULATION_SPREAD_SOURCES + def require_unit_provenance(self, label: str = "observed_distribution") -> None: + """Raise unless a population-spread claim states its unit count and unit type. + + Called by owners that carry unit accounting on the distribution itself + (submodel error models). Calibration targets name a cohort instead. + """ + if self.spread_source not in POPULATION_SPREAD_SOURCES: + return + if self.experimental_unit_type in ( + ExperimentalUnitType.TECHNICAL, + ExperimentalUnitType.CLONAL, + ): + raise ValueError( + f"{label}: experimental_unit_type='{self.experimental_unit_type.value}' cannot " + f"support spread_source='{self.spread_source.value}'. A spread over " + "technical/clonal replicates is not population variability." + ) + if self.n_biological is None: + raise ValueError( + f"{label}: spread_source='{self.spread_source.value}' declares a population " + "spread but n_biological is not set." + ) + if self.experimental_unit_type is None: + raise ValueError( + f"{label}: spread_source='{self.spread_source.value}' declares a population " + "spread but experimental_unit_type is not set." + ) + class TableExcerpt(BaseModel): """ diff --git a/src/maple/core/calibration/submodel_target.py b/src/maple/core/calibration/submodel_target.py index 20d03c1..01b6803 100644 --- a/src/maple/core/calibration/submodel_target.py +++ b/src/maple/core/calibration/submodel_target.py @@ -30,7 +30,6 @@ TableExcerpt, ) - # ============================================================================= # ENUMS # ============================================================================= @@ -961,6 +960,19 @@ def independent_variable(self) -> Optional[IndependentVariable]: """Backwards compatibility: access forward_model.independent_variable.""" return self.forward_model.independent_variable + @model_validator(mode="after") + def _spread_source_states_its_units(self) -> "Calibration": + """A population ``spread_source`` must state its unit count and unit type. + + Submodel targets carry unit accounting on the distribution itself, having no + cohort registry to hold it. + """ + for entry in self.error_model: + od = getattr(entry, "observed_distribution", None) + if od is not None: + od.require_unit_provenance(f"Error model '{entry.name}'") + return self + @model_validator(mode="after") def _unit_groups_consistent(self) -> "Calibration": """Validate ``observed_distribution.unit_group`` across error-model entries. @@ -3041,32 +3053,29 @@ def validate_center_channel_sem_scale(self) -> "SubmodelTarget": @model_validator(mode="after") def validate_bounded_observable_uses_logit_normal(self) -> "SubmodelTarget": - """A bounded observable's population spread (``moments`` form) must use - ``shape: logit_normal``, not normal/lognormal. + """A bounded observable that declares a ``shape`` must use ``logit_normal``. For a fraction / proportion / probability / percent observable, ``normal`` puts mass outside the bound and ``lognormal`` is unbounded above (a near-1 fraction's upper quartile escapes past 1). ``logit_normal`` expands the - quartiles in logit space so they stay in (0, 1). Only applies to the - ``moments`` form — the ``quantiles`` form carries the empirical shape (and - skew) directly and is exempt. + quartiles in logit space so they stay in (0, 1). A distribution with no + declared shape is exempt: it reports its quantiles directly and nothing is + expanded. """ BOUNDED_UNITS = {"percent", "%", "fraction", "proportion", "probability"} for entry in self.calibration.error_model: od = entry.observed_distribution - if od is None or od.moments is None: - continue - if od.moments.shape == DistributionShape.LOGIT_NORMAL: + if od is None or od.shape is None or od.shape == DistributionShape.LOGIT_NORMAL: continue if (entry.units or "").strip().lower() not in BOUNDED_UNITS: continue raise ValueError( f"Error model '{entry.name}' is a bounded observable (units='{entry.units}') " - f"but its observed_distribution.moments uses shape='{od.moments.shape.value}'. " - "Bounded fractions/percentages must use shape='logit_normal', which expands " - "quartiles in logit space so they never escape (0, 1); normal puts mass outside " - "the bound and lognormal is unbounded above. logit_normal requires center in " - "(0, 1) with center_type='median' — express a percent as a fraction (12% -> 0.12)." + f"but its observed_distribution declares shape='{od.shape.value}'. Bounded " + "fractions/percentages must use shape='logit_normal', which expands quartiles " + "in logit space so they never escape (0, 1); normal puts mass outside the " + "bound and lognormal is unbounded above. Express a percent as a fraction " + "(12% -> 0.12)." ) return self @@ -3261,19 +3270,19 @@ def validate_no_invisible_characters(self) -> "SubmodelTarget": # Characters that are invisible or cause issues INVISIBLE_CHARS = { # Zero-width characters - "\u200B", # Zero-width space - "\u200C", # Zero-width non-joiner - "\u200D", # Zero-width joiner - "\uFEFF", # Byte order mark / zero-width no-break space + "\u200b", # Zero-width space + "\u200c", # Zero-width non-joiner + "\u200d", # Zero-width joiner + "\ufeff", # Byte order mark / zero-width no-break space # Soft hyphen - "\u00AD", # Soft hyphen (invisible in most contexts) + "\u00ad", # Soft hyphen (invisible in most contexts) # Other problematic invisibles "\u2060", # Word joiner "\u2061", # Function application "\u2062", # Invisible times "\u2063", # Invisible separator "\u2064", # Invisible plus - "\u180E", # Mongolian vowel separator + "\u180e", # Mongolian vowel separator } def check_invisible(value: str, location: str) -> None: diff --git a/src/maple/core/calibration/test_stats_loader.py b/src/maple/core/calibration/test_stats_loader.py index 0d344c3..b511cf0 100644 --- a/src/maple/core/calibration/test_stats_loader.py +++ b/src/maple/core/calibration/test_stats_loader.py @@ -38,6 +38,14 @@ import pandas as pd import yaml +from maple.core.calibration.cohort import CohortRegistry, load_cohorts +from maple.core.calibration.denominator_audit import ( + check_mapping_collisions, + warn_code_readout_mismatches, + warn_declared_biases, +) +from maple.core.calibration.registry_audit import resolve_n + # Column order used by both loaders (calibration-only columns are the # canonical layout; the prediction loader appends ``is_prediction_only`` # for callers that need to distinguish rows after concat). @@ -104,7 +112,37 @@ def _gather_yaml_files(dirs: List[Path]) -> List[Path]: return [seen[name] for name in sorted(seen)] -def load_calibration_targets(yaml_dir: Path | str | List) -> pd.DataFrame: +def _as_registry(cohorts: Path | str | CohortRegistry | None) -> CohortRegistry | None: + if cohorts is None or isinstance(cohorts, CohortRegistry): + return cohorts + return load_cohorts(cohorts) + + +def _resolve_sample_size( + target_id: str, data: dict, registry: CohortRegistry | None +) -> float | int: + """The patients behind this target: its own ``sample_size``, else the cohort's.""" + declared = (data.get("empirical_data") or {}).get("sample_size") + if declared is not None: + return declared + if registry is None: + raise ValueError( + f"'{target_id}' declares no sample_size, so its n lives on the cohort, but no " + "cohort registry was passed. Call load_calibration_targets(..., cohorts=...) " + "with the path to cohorts.yaml." + ) + resolved = resolve_n(data, registry) + if resolved is None: + raise ValueError( + f"'{target_id}' declares no sample_size and its cohort_id " + f"'{data.get('cohort_id')}' does not resolve in the registry." + ) + return resolved + + +def load_calibration_targets( + yaml_dir: Path | str | List, cohorts: Path | str | CohortRegistry | None = None +) -> pd.DataFrame: """ Load calibration target YAMLs from one or more directories into a test-statistics DataFrame. @@ -124,6 +162,11 @@ def load_calibration_targets(yaml_dir: Path | str | List) -> pd.DataFrame: ``["calibration_targets/clinical_progression", "calibration_targets/mechanistic/clinical_progression"]``. + cohorts: The cohort registry, as a path to ``cohorts.yaml`` or a loaded + ``CohortRegistry``. Required whenever a literature target is loaded: + those carry no ``sample_size`` of their own, and n resolves to + ``n_evaluable`` else the named cohort's ``n_c``. + Returns: DataFrame with columns: test_statistic_id, required_species, model_output_code, @@ -145,10 +188,21 @@ def load_calibration_targets(yaml_dir: Path | str | List) -> pd.DataFrame: joined = ", ".join(str(d) for d in dirs) raise ValueError(f"No YAML files found in {joined}") - rows: List[dict] = [] + parsed: dict[str, dict] = {} for yaml_file in yaml_files: with open(yaml_file, "r") as f: - data = yaml.safe_load(f) + parsed[yaml_file.name] = yaml.safe_load(f) + + # Cross-target checks. These cannot live on the pydantic model, which only + # ever sees one target: a mapping collision is a property of a pair. + check_mapping_collisions(parsed) + warn_declared_biases(parsed) + warn_code_readout_mismatches(parsed) + + registry = _as_registry(cohorts) + rows: List[dict] = [] + for yaml_file in yaml_files: + data = parsed[yaml_file.name] target_id = data["calibration_target_id"] observable = data["observable"] @@ -180,7 +234,7 @@ def load_calibration_targets(yaml_dir: Path | str | List) -> pd.DataFrame: ci95_upper = float("nan") units = observable.get("units", "") - sample_size = empirical.get("sample_size", float("nan")) + sample_size = _resolve_sample_size(target_id, data, registry) # Generate wrapper code. The time-series reduction is declared on the # observable: readout_time (interpolate at that time) XOR diff --git a/src/maple/core/calibration/validators.py b/src/maple/core/calibration/validators.py index 73a84c1..45de465 100644 --- a/src/maple/core/calibration/validators.py +++ b/src/maple/core/calibration/validators.py @@ -15,7 +15,6 @@ import numpy as np import requests - # ============================================================================= # MODULE-LEVEL CACHE FOR PAPER TEXTS # ============================================================================= diff --git a/src/maple/prompts/calibration_target_prompt.md b/src/maple/prompts/calibration_target_prompt.md index 1a9a11f..4a807dd 100644 --- a/src/maple/prompts/calibration_target_prompt.md +++ b/src/maple/prompts/calibration_target_prompt.md @@ -25,10 +25,11 @@ You produce a single `CalibrationTarget` YAML with these top-level fields. **Req | `study_interpretation` | str | One-paragraph narrative: what is being measured, how the experimental context maps to the model, key methodological points | | `key_assumptions` | List[str] (≥1) | Biological + statistical assumptions made in extraction (e.g., cell-type equivalence, normal-distribution assumption) | | `key_study_limitations` | List[str] | Issues that bias estimates or limit generalizability (e.g., small cohort, figure-extracted values) | +| `cohort_id` | str | The registered cohort in `calibration_targets/cohorts.yaml` whose patients this measures. Required when `epistemic_basis: literature`. See "Cohort" below | | `observable` | dict | How to compute the observable from QSP species. See "Observable" section below | | `experimental_context` | dict | The source paper's context: `{species, system, indication, treatment, stage?, mouse_subspecifier?, cell_lines?, culture_conditions?, tissue_source?, assay_type?}` — describes WHERE the data came from, not the model target | | `scenario` | dict (optional) | Interventions + measurement timing. Omit if untreated baseline and no perturbations | -| `empirical_data` | dict | Computed `median`, `ci95`, `units`, `sample_size`, `inputs[]`, `assumptions[]`, `distribution_code`. See "Empirical Data" section | +| `empirical_data` | dict | Computed `median`, `ci95`, `units`, `inputs[]`, `assumptions[]`, `distribution_code`. See "Empirical Data" section | | `primary_data_source` | dict | Single paper with verified DOI. See "Source Requirements". Required when `epistemic_basis: literature` (default) | | `secondary_data_sources` | List[dict] | Reference values / conversion-factor sources. Empty list OK | | `epistemic_basis` | `"literature"` (default) or `"mechanistic"` | Use `"mechanistic"` only for biological-invariant priors with no primary measurement (live in `calibration_targets/mechanistic/`); requires deliberately wide CIs and rationale in `key_assumptions`. Otherwise leave at default | @@ -179,43 +180,49 @@ observable: rationale: "ORR per RECIST 1.1 (≥30% decrease in longest diameter)" ``` -### Population Spread (across-patient variability) +### Cohort -For hierarchical / virtual-patient inference, a target declares whether its reported width is genuine patient-to-patient spread (usable as the population-spread / omega signal) or just uncertainty on the mean. Default is `center_only` (excluded from omega); opt in explicitly. Two interoperable ways to declare it: +`cohort_id` names one entry of `calibration_targets/cohorts.yaml`: one study's patients, measured once. It is required for every `epistemic_basis: literature` target and must be omitted for `mechanistic` ones, which assert a constraint rather than measuring anyone. -1. **`population_spread` + a `samples` array.** Set `empirical_data.population_spread: across_patient` and have `distribution_code` ALSO return a `samples` key — the across-patient population draw (one value per patient-equivalent); its empirical spread is the omega signal. Keep the default `center_only` (and do NOT return `samples`) when the width is a pooled-mean / SEM CI that shrinks with n. `median_obs` / `ci95` are unaffected either way. +Population inference forms one covariance block per cohort, so the cohort decides which reported statistics are treated as correlated. Consequences when authoring: -2. **`observed_distribution`** (general representation, shared with submodel targets). Author it in whichever form the paper reports — **prefer `moments`** (mean +/- SD, median +/- IQR, CV, CI); the framework expands it to quartiles, so do not hand-convert: +- **One target, one source.** A target whose `inputs` cite several papers is a meta-analysis, not a cohort: its patients were never measured together, so there is no resampling distribution over them. Split it into one target per source. A validator rejects the pooled form. +- **Do not invent a cohort per readout.** Ten measurements on the same 9 patients are ten targets naming one cohort. +- **Per-readout evaluability goes on the target.** When fewer patients contributed to this statistic than the cohort holds (paired fold changes evaluable in 6 of 9), set `empirical_data.n_evaluable`. The cohort's `n_c` stays the patient count. +- If the cohort you need is not in the registry, add it there rather than inventing the id. + +### Reported Distribution + +`empirical_data.observed_distribution` records what the source actually printed. It is required for `epistemic_basis: literature`, and it is the sole declaration of whether the reported width is genuine patient-to-patient variability. + +**Record the statistics on the page. Do not convert between them.** A median with quartiles is three entries; a mean with an SD is two. Converting an SE to an SD, or a mean+SD to quartiles, destroys the distinction the consumer needs. ```yaml empirical_data: - population_spread: across_patient observed_distribution: - moments: - center: 17 - center_type: median - scale: 21 # full IQR here - scale_type: iqr # sd | sem | cv | iqr | ci95_halfwidth - shape: lognormal # lognormal | normal | logit_normal (for [0,1] fractions) - spread_source: across_patient # or center_only (SEM/CI on the mean; default) - n_biological: 40 - experimental_unit_type: biological + statistics: + - {stat: quantile, p: 0.5, value: 17} + - {stat: quantile, p: 0.25, value: 9} + - {stat: quantile, p: 0.75, value: 30} + spread_source: across_patient ``` -Use the `quantiles` form when the paper gives quartiles/percentiles/samples directly: - ```yaml observed_distribution: - quantiles: - - {p: 0.25, value: 9} - - {p: 0.5, value: 17} - - {p: 0.75, value: 30} + statistics: + - {stat: mean, value: 227.7} + - {stat: se, value: 15.0} # an SE, because that is what the paper printed spread_source: across_patient - n_biological: 40 - experimental_unit_type: biological + shape: lognormal ``` -Provide EXACTLY ONE of `moments` / `quantiles`. A population spread (`spread_source: across_patient`) REQUIRES `n_biological` + `experimental_unit_type: biological`. When both `observed_distribution` and `population_spread` are present they must agree that the width is (or is not) genuine spread — a validator enforces this. For a `[0,1]`-bounded observable (a fraction / proportion / response rate), use `shape: logit_normal` (with `center_type: median`) so expanded quartiles stay in-bounds — **a validator rejects a `percent`/`fraction`/`proportion`/`probability` observable whose `moments` shape is `normal` or `lognormal`.** If the source reports the unit count as a LOWER BOUND ("n≥40", "at least 40 patients"), set `n_biological_is_floor: true` (with the floor value for `n_biological`) so precision-weighting consumers don't over-weight the panel. If a target ever carries several observables measured on the SAME patients/units and you want them treated as one shared spread, tag them with a common `unit_group` string (rare for full-model calibration targets, which are usually one observable each — omit it otherwise). +`stat` is one of: `quantile` (needs `p`; the median is `p: 0.5`) · `mean` · `geometric_mean` · `sd` · `cv` · `iqr` · `range` · `se` · `ci95_lo` · `ci95_hi` · `min` · `max`. + +`spread_source` says what the width measures: `across_patient` (inter-individual spread in the target population) · `biological_experimental` · `translation` · `technical` · `assumed` · `center_only` (an SE or CI on the mean, which shrinks with n and is not spread). A `spread_source` that claims population spread requires at least one width statistic. + +`shape` is needed only when the source printed no quartiles and a consumer has to expand a center and a scale into them. Recording it here makes the assumption explicit rather than silent. Use `logit_normal` for a `[0,1]`-bounded observable so expanded quartiles cannot escape the bound; `lognormal` on a near-1 fraction pushes the upper quartile past 1. + +**`n_biological`, `n_technical`, `experimental_unit_type`, `unit_group` are submodel-target fields and must not appear on a calibration target.** The cohort owns the unit accounting; a validator rejects them here. ### Source Relevance Assessment @@ -296,10 +303,12 @@ When a paper reports `mean ± value`, you must determine whether the dispersion **Encoding:** -| Identified as | Input name prefix | `dispersion_type` | In `distribution_code` | +| Identified as | Input name prefix | `dispersion_type` | `observed_distribution.statistics` | |---|---|---|---| -| SD | `sd_*` | `sd` | use directly | -| SEM | `sem_*` | `se` | convert: `sd = sem * np.sqrt(n)` | +| SD | `sd_*` | `sd` | `{stat: sd, value: ...}` | +| SEM | `sem_*` | `se` | `{stat: se, value: ...}` | + +Record whichever was printed. `distribution_code` may still scale an SEM to an SD to fit the center's CI, but `observed_distribution` must keep the raw statistic: which one the paper printed is data, what a consumer does with it is that consumer's decision, and they differ by √n (a factor of 3.2 at n=10). Always populate `dispersion_type_rationale` with your evidence (label, √n test arithmetic, or CV check). @@ -339,6 +348,8 @@ When the calibration target requests a fold-change (pre-to-post treatment change 3. **If no paired pre/post data exists** for the requested comparison, clearly state this limitation rather than substituting a cross-arm ratio. +4. **Declare the reference.** Set `observable.readout.quantity_kind: foldchange` and `observable.readout.reference` to whatever the change is measured against — a timepoint of the same scenario, or another scenario. Both are required for a fold change and the reference must not be left implicit in the code. + ### Scaling Factor Red Flag **CRITICAL:** If you find yourself needing a dimensionless scaling factor > 10× in `observable.constants` to reconcile model output with literature data, this almost certainly indicates an extraction error. @@ -462,7 +473,25 @@ The top-level `observable:` block describes a single experimental observable. (T - `positive_unbounded`: Output must be > 0, no upper bound (fold-changes, ratios) - `real`: Any real value (log-ratios, change scores) -7. **observable.experimental_denominator** / **observable.model_denominator_species** — describe what the experiment normalizes by, and which model species form the matching model-side denominator. **Conditionally required:** when the observable is a density or per-mass concentration (units like `cell/mm**2`, `pg/mg`, `cell/g`, etc., with `support: positive`), validation REQUIRES `experimental_denominator` to be set. Omitting it triggers a `value_error: "Observable with units='pg/mg' and support='positive' is a density but experimental_denominator is not set"` failure. Optional only for unitless ratios (`support: unit_interval`) or absolute counts. +7. **observable.readout** (required for `epistemic_basis: literature`; omit for `mechanistic`) — what the EXPERIMENT measured, as opposed to how the model computes it. The rest of `observable` is the computation; this block is the measurement, and inference reads it to decide which rows share a measurement correction and which rows are the same quantity. A mechanistic target ran no assay, so it has no readout. + +```yaml +readout: + quantity_kind: fraction + assay_modality: mihc + numerator_species: ['V_T.CD8'] + denominator_species: ['V_T.CD8', 'V_T.Th', 'V_T.Treg'] + experimental_denominator: "CD3+ T cells (all T cell subsets)" + reference: # fold changes only; omit otherwise +``` + + - **`quantity_kind`**: `density` (cells or mass per area/volume of tissue) · `fraction` (a part over a whole containing it, in (0,1)) · `ratio` (one quantity over another not containing it) · `foldchange` (a quantity over its own value at a reference) · `concentration` · `intensity` · `time`. A percentage is a `fraction`, not its own kind; the compartment it is a percentage OF is named by the denominator fields. + - **`assay_modality`**: `mihc` · `ihc` · `mif` · `flow_cytometry` · `scrnaseq` · `multiplex_immunoassay` · `elisa` · `multiphoton_microscopy` · `digital_histopathology` · `imaging` · `clinical`. Describe the measurement, not the cell type: two readouts agreeing on kind and modality share one correction whatever they count. + - **`numerator_species`** (required) and **`denominator_species`** (empty for an absolute quantity) — which model species compose the row. This is the row's identity, so it must not depend on how you spelled the arithmetic in `code`: declare the components even when the code divides by a single pre-aggregated species, and expand that aggregate to the species the model currently sums into it. A cross-target audit holds `code` to this declaration and reports any disagreement. + - **`experimental_denominator`** — what the experiment divided by, in the paper's own words. **Conditionally required:** an observable that is a density or per-mass concentration (units like `cell/mm**2`, `pg/mg`, with `support: positive`) is rejected without it. Setting it requires `denominator_species` too. + - **`reference`** (required for `quantity_kind: foldchange`, forbidden otherwise) — what the fold change is measured against. Do not rely on a `species[0]` or `series[0]` convention inside the code. + - An earlier time in the same scenario: `{kind: timepoint, timepoint: 0.0, timepoint_unit: day}` + - Another scenario: `{kind: scenario, scenario: baseline_no_treatment}` 8. **observable.readout_time / observable.reduce_observable** (**required — set EXACTLY ONE**) — how the observable time-series from `observable.code` is reduced to the single scalar compared against `empirical_data`. There is NO implicit default; a target that sets neither is rejected. - **`readout_time`** (float) + **`readout_time_unit`** (str, e.g. `day`) — the common case: the series is linearly interpolated to this simulation time. For a treatment-arm biopsy at day 21, use `readout_time: 21.0`, `readout_time_unit: day`. For a **baseline / diagnosis-snapshot** measurement (treatment-naive resection, the reference state), set `readout_time: 0.0` explicitly (the trajectory's `t=0` is diagnosis). @@ -483,6 +512,12 @@ observable: return (cd8 / tumor_area).to('cell/mm**2') units: cell/mm**2 species: ['V_T.CD8', 'V_T.C1'] + readout: + quantity_kind: density + assay_modality: ihc + numerator_species: ['V_T.CD8'] + denominator_species: ['V_T.C1'] + experimental_denominator: "mm^2 of tumor tissue section (cancer cells + stroma)" constants: - name: area_per_cancer_cell value: 2.27e-4 @@ -497,8 +532,6 @@ observable: source_type: reference_db reference_db_name: pdac_stromal_fraction support: positive - experimental_denominator: "mm^2 of tumor tissue section (cancer cells + stroma)" - model_denominator_species: ['V_T.C1'] ``` **Example observable (dimensionless ratio — preferred when available):** @@ -512,10 +545,14 @@ observable: return (treg / total_t).to('dimensionless') units: dimensionless species: ['V_T.Treg', 'V_T.CD8', 'V_T.Th', 'V_T.CD8_exh', 'V_T.Th_exh'] + readout: + quantity_kind: fraction + assay_modality: mihc + numerator_species: ['V_T.Treg'] + denominator_species: ['V_T.Treg', 'V_T.CD8', 'V_T.Th', 'V_T.CD8_exh', 'V_T.Th_exh'] + experimental_denominator: "CD3+ T cells (all T cell subsets)" constants: [] support: unit_interval - experimental_denominator: "CD3+ T cells (all T cell subsets)" - model_denominator_species: ['V_T.Treg', 'V_T.CD8', 'V_T.Th', 'V_T.CD8_exh', 'V_T.Th_exh'] ``` **Example observable (compartment bridge via auxiliary parameter — serum→tumor TGFβ):** @@ -534,6 +571,10 @@ observable: return predicted_serum.to('nanomolar') units: nanomolar species: ['V_T.TGFb'] + readout: + quantity_kind: concentration + assay_modality: multiplex_immunoassay + numerator_species: ['V_T.TGFb'] constants: - name: active_fraction value: 0.10 @@ -684,7 +725,11 @@ When `observable.code` includes conversion factors (cells → volume, IHC score ### Sample Size -`empirical_data.sample_size` (int) and `empirical_data.sample_size_rationale` (str) are required. Look for `n =`, `N =`, figure legends, patient counts, replicate counts. If unreported, back-calculate from SD/SEM (if both given), or use a conservative type-based estimate and note the uncertainty in `sample_size_rationale`. +**A literature target must NOT set `sample_size`.** Its cohort carries the patient count as `n_c`, and stating it twice is how the two drift apart. Consumers resolve n as `n_evaluable` if set, else the cohort's `n_c`. A validator rejects `sample_size` on an `epistemic_basis: literature` target. + +`empirical_data.sample_size` (int) and `empirical_data.sample_size_rationale` (str) are **required for `epistemic_basis: mechanistic` only**, which names no cohort because it measures nobody. + +`empirical_data.n_evaluable` (int, optional) is the patient count behind THIS statistic when it is smaller than the cohort's `n_c` — paired fold changes evaluable in 6 of 9 arm patients, a stain scorable on 7 of 9 slides. Set it whenever the source states a per-readout denominator; leave it null when the whole cohort contributed. It may never exceed the cohort's `n_c`. **Reference example with all `empirical_data` fields:** ```yaml @@ -692,8 +737,7 @@ empirical_data: median: [149.94] # List[float] — length-1 (scalar target) ci95: [[100.79, 199.35]] # List[List[float]] — a single [lo, hi] pair units: cell / mm**2 - sample_size: 42 - sample_size_rationale: "n=42 patients with resected PDAC tumors, stated in Table 1" + n_evaluable: 38 # optional; omit when the whole cohort contributed inputs: [...] # List[EstimateInput] — paper-reported values used by distribution_code assumptions: [...] # Optional List[ModelingAssumption] — values NOT from the paper but needed for computation (e.g., n_mc_samples=10000, assumed_cv when not reported); each requires a rationale distribution_code: | diff --git a/src/maple/prompts/submodel_target_prompt.md b/src/maple/prompts/submodel_target_prompt.md index 2d3cd9d..e86a60c 100644 --- a/src/maple/prompts/submodel_target_prompt.md +++ b/src/maple/prompts/submodel_target_prompt.md @@ -331,34 +331,32 @@ def derive_observation(inputs, sample_size, rng, n_bootstrap): `observation_code` informs the parameter's CENTER. To ALSO inform how much patients *differ* — the population-spread (omega) signal used by virtual-patient inference — add an `observed_distribution` block to the error_model entry whenever the source reports a genuine spread. It is OPTIONAL and additive: omit it and the target behaves exactly as before. -Record the reported distribution in whichever form the paper gives, plus a provenance tag. **Prefer the `moments` form** — most papers report mean +/- SD (or median +/- IQR, CV, CI), and the framework expands it to quartiles for you. Do NOT hand-convert mean +/- SD into quartiles. +**Record the statistics the paper printed, as a flat list. Do NOT convert between them.** A mean with an SD is two entries; a median with quartiles is three. Converting an SE to an SD, or a mean+SD to quartiles, throws away the distinction the consumer needs. ```yaml observed_distribution: - moments: - center: 0.178 - center_type: mean # or median - scale: 0.09 - scale_type: sd # sd | sem | cv | iqr | ci95_halfwidth - shape: lognormal # lognormal | normal (how to expand center+scale) + statistics: + - {stat: mean, value: 0.178} + - {stat: sd, value: 0.09} spread_source: biological_experimental + shape: lognormal # only needed if a consumer must expand these into quartiles n_biological: 45 experimental_unit_type: biological ``` -Use the `quantiles` form only when the source reports quartiles/percentiles/samples directly: - ```yaml observed_distribution: - quantiles: - - {p: 0.25, value: 0.12} - - {p: 0.5, value: 0.178} - - {p: 0.75, value: 0.25} + statistics: + - {stat: quantile, p: 0.5, value: 0.178} + - {stat: quantile, p: 0.25, value: 0.12} + - {stat: quantile, p: 0.75, value: 0.25} spread_source: biological_experimental n_biological: 45 experimental_unit_type: biological ``` +`stat` is one of: `quantile` (needs `p`; the median is `p: 0.5`) · `mean` · `geometric_mean` · `sd` · `cv` · `iqr` · `range` · `se` · `ci95_lo` · `ci95_hi` · `min` · `max`. + **Choose `spread_source` by what the reported width actually measures:** | Reported width | `spread_source` | Feeds population spread? | @@ -370,14 +368,14 @@ observed_distribution: | No spread reported | omit the field | No — a wide default is used | **Rules (validated):** -- Provide EXACTLY ONE of `moments` or `quantiles`. +- A `spread_source` that claims population spread requires at least one WIDTH statistic (quartiles, or an `sd` / `iqr` / `cv` / `range`). An entry carrying only a center must declare `center_only`. - A population spread (`biological_experimental` / `across_patient`) REQUIRES `n_biological` and `experimental_unit_type: biological`. A spread over technical/clonal replicates is not population variability — use `technical` or `center_only`. -- `scale_type: sem` recovers the population SD as `SEM * sqrt(n_biological)` — so give `n_biological`. +- Record an SE as `{stat: se, ...}`; a consumer recovers the population SD as `SE * sqrt(n_biological)`, so give `n_biological`. Never do that conversion here. - If the source reports the unit count as a LOWER BOUND ("n≥8", "at least 8 donors", "n=8–12 across conditions"), set `n_biological_is_floor: true` and use the floor value for `n_biological`. An exact-looking n from a floor over-states precision and over-weights the panel in the finite-sample / inverse-variance weighting downstream. - Keep values in paper units; do unit conversion in code, not in the anchors. - Most in-vitro submodel data is a donor/animal spread that is a LOWER BOUND on PDAC patient spread — grade that transfer in `source_relevance.heterogeneity_transfer` (see Source Relevance below). - `observation_code` still just pins the CENTER (keep it SEM-scale); `observed_distribution` is the single source of population-spread truth. Do not double-encode spread in both. -- For a `[0,1]`-bounded observable (a fraction, proportion, or probability — e.g. a polarization fraction, a positive-cell %), use `shape: logit_normal` instead of `lognormal` in the `moments` form. It expands the quartiles in logit space so they can never escape `(0, 1)`; `lognormal` on a near-1 fraction would push the upper quartile past 1. `logit_normal` requires `center_type: median`. **This includes `percent` observables** — express the value as a fraction in `(0, 1)` first (12% → 0.12), since `logit_normal` needs the center inside the unit interval. A validator rejects a `percent`/`fraction`/`proportion`/`probability` observable whose `moments` shape is `normal` or `lognormal`. (The `quantiles` form is exempt — explicit anchors carry the empirical shape and skew directly, and are often the better choice for a visibly-skewed donor distribution.) +- `shape` is needed ONLY when the source printed no quartiles and a consumer has to expand a center and a scale into them. Recording it makes the assumption explicit rather than silent; omit it when you gave quartiles directly. For a `[0,1]`-bounded observable (a fraction, proportion, or probability — a polarization fraction, a positive-cell %), use `shape: logit_normal`: it expands in logit space so the quartiles can never escape `(0, 1)`, where `lognormal` on a near-1 fraction would push the upper quartile past 1. **This includes `percent` observables** — express the value as a fraction in `(0, 1)` first (12% → 0.12). A validator rejects a `percent`/`fraction`/`proportion`/`probability` observable that declares `shape: normal` or `lognormal`. **`unit_group` — only for multi-observable, SAME-unit targets.** When a single target has SEVERAL `error_model` entries measured on the SAME biological units (the same donor panel across the doses of a dose-response; one cohort followed over a time course), tag those entries with a shared `unit_group` string. This tells the hierarchical layer they share ONE biological random effect, so it moment-matches them jointly instead of treating each dose/timepoint as an independent measurement (which would spuriously shrink the population spread by ~sqrt(number-of-points)). @@ -386,14 +384,16 @@ error_model: - name: IFNg_at_0p5_kPa observed_distribution: unit_group: saitakis_donors # same 13 donors measured at every stiffness - moments: {center: 370.5, center_type: mean, scale: 316.4, scale_type: sd, shape: lognormal} + statistics: [{stat: mean, value: 370.5}, {stat: sd, value: 316.4}] + shape: lognormal spread_source: biological_experimental n_biological: 13 experimental_unit_type: biological - name: IFNg_at_2p0_kPa observed_distribution: unit_group: saitakis_donors # <- same tag: shared donor panel - moments: {center: 210.0, center_type: mean, scale: 180.0, scale_type: sd, shape: lognormal} + statistics: [{stat: mean, value: 210.0}, {stat: sd, value: 180.0}] + shape: lognormal spread_source: biological_experimental n_biological: 13 experimental_unit_type: biological diff --git a/tests/unit/core/test_calibration_target_validators.py b/tests/unit/core/test_calibration_target_validators.py index db5b0ad..277b870 100644 --- a/tests/unit/core/test_calibration_target_validators.py +++ b/tests/unit/core/test_calibration_target_validators.py @@ -41,7 +41,6 @@ from maple.core.calibration import CalibrationTarget, Observable from maple.core.calibration.calibration_target_models import CalibrationTargetEstimates - DEFAULT_CLINICAL_SOURCE_RELEVANCE = { "indication_match": "exact", "indication_match_justification": "Human PDAC resection specimens with quantitative IHC, directly matching the model indication.", @@ -103,7 +102,14 @@ def model_structure(): def golden_calibration_target_data(): """Complete valid CalibrationTarget data that passes all validators.""" return { + "cohort_id": "smith2020_resected", "observable": { + "readout": { + "quantity_kind": "ratio", + "assay_modality": "ihc", + "numerator_species": ["V_T.CD8"], + "denominator_species": ["V_T.C1"], + }, "code": ( "def compute_observable(time, species_dict, constants):\n" " cd8 = species_dict['V_T.CD8']\n" @@ -146,8 +152,14 @@ def golden_calibration_target_data(): "median": [1.0], "ci95": [[0.3737, 2.7]], "units": "dimensionless", - "sample_size": 42, - "sample_size_rationale": "n=42 patients in resected PDAC cohort, Table 1", + "observed_distribution": { + "statistics": [ + {"stat": "mean", "value": 1.0}, + {"stat": "sd", "value": 0.5}, + ], + "shape": "lognormal", + "spread_source": "center_only", + }, "inputs": [ { "name": "cd8_ratio_mean", @@ -257,6 +269,76 @@ def test_golden_yaml_passes_all_validators( assert target.observable.reduce_observable is None +class TestReadoutIsForMeasurements: + """A readout describes an assay, so only a literature target has one.""" + + def test_literature_target_requires_a_readout( + self, model_structure, golden_calibration_target_data, mock_crossref_success + ): + data = copy.deepcopy(golden_calibration_target_data) + data["observable"].pop("readout") + with pytest.raises(ValidationError, match="observable.readout"): + CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) + + def test_mechanistic_target_needs_none( + self, model_structure, golden_calibration_target_data, mock_crossref_success + ): + data = copy.deepcopy(golden_calibration_target_data) + data["epistemic_basis"] = "mechanistic" + data.pop("cohort_id") + data["observable"].pop("readout") + data["empirical_data"]["observed_distribution"] = None + data["empirical_data"]["sample_size"] = 1 + data["empirical_data"]["sample_size_rationale"] = "Asserted, not measured." + target = CalibrationTarget.model_validate( + data, context={"model_structure": model_structure} + ) + assert target.observable.readout is None + + +class TestSampleSizeOwnership: + """Whoever owns the patients owns the count.""" + + def test_literature_target_rejects_sample_size( + self, model_structure, golden_calibration_target_data, mock_crossref_success + ): + data = copy.deepcopy(golden_calibration_target_data) + data["empirical_data"]["sample_size"] = 42 + data["empirical_data"]["sample_size_rationale"] = "Table 1" + with pytest.raises(ValidationError, match="its cohort already carries n_c"): + CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) + + def test_literature_target_may_declare_n_evaluable( + self, model_structure, golden_calibration_target_data, mock_crossref_success + ): + data = copy.deepcopy(golden_calibration_target_data) + data["empirical_data"]["n_evaluable"] = 6 + target = CalibrationTarget.model_validate( + data, context={"model_structure": model_structure} + ) + assert target.empirical_data.n_evaluable == 6 + assert target.empirical_data.sample_size is None + + def test_mechanistic_target_requires_sample_size( + self, model_structure, golden_calibration_target_data, mock_crossref_success + ): + data = copy.deepcopy(golden_calibration_target_data) + data["epistemic_basis"] = "mechanistic" + data.pop("cohort_id") + with pytest.raises(ValidationError, match="requires sample_size"): + CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) + + def test_sample_size_requires_its_rationale( + self, model_structure, golden_calibration_target_data, mock_crossref_success + ): + data = copy.deepcopy(golden_calibration_target_data) + data["epistemic_basis"] = "mechanistic" + data.pop("cohort_id") + data["empirical_data"]["sample_size"] = 42 + with pytest.raises(ValidationError, match="requires sample_size_rationale"): + CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) + + class TestCalibrationTargetContextOptional: """Model-structure-dependent validators must DEFER when no validation context is supplied, and still ENFORCE when it is. @@ -1530,6 +1612,11 @@ class TestObservableDenominatorAudit: def _make_observable(self, **overrides): """Helper to create Observable with sensible defaults.""" base = { + "readout": { + "quantity_kind": "density", + "assay_modality": "ihc", + "numerator_species": ["V_T.CD8"], + }, "code": ( "def compute_observable(time, species_dict, constants):\n" " return species_dict['V_T.CD8']" @@ -1542,14 +1629,16 @@ def _make_observable(self, **overrides): "readout_time": 0.0, "readout_time_unit": "day", } + readout = dict(base["readout"], **overrides.pop("readout", {})) base.update(overrides) + base["readout"] = readout return Observable(**base) def test_observable_without_denominator_fields_passes(self): """Observable without denominator fields passes for non-density units.""" obs = self._make_observable() - assert obs.experimental_denominator is None - assert obs.model_denominator_species is None + assert obs.readout.experimental_denominator is None + assert obs.readout.denominator_species == [] def test_density_observable_without_experimental_denominator_fails(self): """Density observable (cell/mm**2) must declare experimental_denominator.""" @@ -1564,17 +1653,19 @@ def test_density_observable_with_denominator_audit_passes(self): obs = self._make_observable( units="cell / millimeter**2", support="positive", - experimental_denominator="mm^2 of tumor tissue (whole section including stroma)", - model_denominator_species=["V_T.C1"], + readout={ + "experimental_denominator": "mm^2 of tumor tissue (whole section)", + "denominator_species": ["V_T.C1"], + }, ) - assert obs.experimental_denominator is not None - assert obs.model_denominator_species == ["V_T.C1"] + assert obs.readout.experimental_denominator is not None + assert obs.readout.denominator_species == ["V_T.C1"] def test_experimental_denominator_without_model_species_fails(self): """Setting experimental_denominator without model_denominator_species fails.""" - with pytest.raises(ValidationError, match="model_denominator_species"): + with pytest.raises(ValidationError, match="denominator_species"): self._make_observable( - experimental_denominator="CD3+ T cells", + readout={"experimental_denominator": "CD3+ T cells"}, ) def test_fraction_with_full_denominator_audit_passes(self): @@ -1582,8 +1673,10 @@ def test_fraction_with_full_denominator_audit_passes(self): obs = self._make_observable( units="dimensionless", support="unit_interval", - experimental_denominator="all cells in ROI (all nucleated cells)", - model_denominator_species=["V_T.CD8", "V_T.Th", "V_T.Treg", "V_T.Mac_M1"], + readout={ + "experimental_denominator": "all cells in ROI (all nucleated cells)", + "denominator_species": ["V_T.CD8", "V_T.Th", "V_T.Treg", "V_T.Mac_M1"], + }, unmodeled_denominator_components=( "B cells (50-70% of LA cells) not modeled; model prediction " "will be ~2-3x higher than experimental value." @@ -1597,11 +1690,11 @@ def test_non_density_units_with_slash_no_cell_passes(self): units="nanomolarity", support="positive", ) - assert obs.experimental_denominator is None + assert obs.readout.experimental_denominator is None class TestCalibrationTargetPopulationSample: - """The optional declared 'samples' population draw + population_spread gate.""" + """The declared 'samples' population draw + the spread_source gate.""" # A lognormal population draw whose median matches the golden reported median (1.0). # @@ -1640,26 +1733,28 @@ def test_defaults_center_only( target = CalibrationTarget.model_validate( golden_calibration_target_data, context={"model_structure": model_structure} ) - assert target.empirical_data.population_spread == "center_only" + assert target.empirical_data.feeds_population_spread is False def test_across_patient_with_samples_validates( self, model_structure, golden_calibration_target_data, mock_crossref_success ): - data = self._with_code( - golden_calibration_target_data, self._GOOD_CODE, population_spread="across_patient" - ) + data = self._with_code(golden_calibration_target_data, self._GOOD_CODE) + _spread_source(data, "across_patient") target = CalibrationTarget.model_validate( data, context={"model_structure": model_structure} ) - assert target.empirical_data.population_spread == "across_patient" + assert target.empirical_data.feeds_population_spread is True - def test_across_patient_without_samples_rejected( + def test_across_patient_without_samples_accepted( self, model_structure, golden_calibration_target_data, mock_crossref_success ): - # Declaring across_patient but returning no samples is a hard error. - data = self._with_code(golden_calibration_target_data, population_spread="across_patient") - with pytest.raises(ValidationError, match="requires distribution_code to"): - CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) + # A reported width is a row of the observation vector, not an error bar, so + # declaring across_patient does not oblige the target to hand over an array. + data = _spread_source(self._with_code(golden_calibration_target_data), "across_patient") + target = CalibrationTarget.model_validate( + data, context={"model_structure": model_structure} + ) + assert target.empirical_data.feeds_population_spread is True def test_center_only_with_samples_rejected( self, model_structure, golden_calibration_target_data, mock_crossref_success @@ -1684,9 +1779,8 @@ def test_samples_median_mismatch_rejected( " out['samples'] = samples * 5.0\n" " return out", ) - data = self._with_code( - golden_calibration_target_data, code, population_spread="across_patient" - ) + data = self._with_code(golden_calibration_target_data, code) + _spread_source(data, "across_patient") with pytest.raises(ValidationError, match=r"median\(samples\)"): CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) @@ -1700,9 +1794,8 @@ def test_degenerate_samples_rejected( " out['samples'] = np.ones(10000) * mean.magnitude * mean.units\n" " return out", ) - data = self._with_code( - golden_calibration_target_data, code, population_spread="across_patient" - ) + data = self._with_code(golden_calibration_target_data, code) + _spread_source(data, "across_patient") with pytest.raises(ValidationError, match="zero variance"): CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) @@ -1713,6 +1806,12 @@ def test_degenerate_samples_rejected( # ============================================================================ +def _spread_source(data: dict, source: str) -> dict: + """Set the target's spread provenance. Unit accounting lives on the cohort.""" + data["empirical_data"]["observed_distribution"]["spread_source"] = source + return data + + def _bounded_cal_estimates(shape: str) -> dict: return { "median": [0.5], @@ -1736,22 +1835,19 @@ def _bounded_cal_estimates(shape: str) -> dict: " v = inputs['resp_fraction']\n" " return {'median_obs': v, 'ci95_lower': v * 0.6, 'ci95_upper': v * 1.4}" ), - "population_spread": "center_only", "observed_distribution": { - "moments": { - "center": 0.5, - "center_type": "median", - "scale": 0.1, - "scale_type": "sd", - "shape": shape, - }, + "statistics": [ + {"stat": "quantile", "p": 0.5, "value": 0.5}, + {"stat": "sd", "value": 0.1}, + ], + "shape": shape, "spread_source": "center_only", }, } class TestCalBoundedObservableLogitNormal: - """CalibrationTargetEstimates: bounded moments-form observable must use logit_normal.""" + """CalibrationTargetEstimates: a bounded observable declaring a shape must use logit_normal.""" def test_percent_with_normal_shape_raises(self): with pytest.raises(ValidationError, match="logit_normal"): @@ -1791,7 +1887,6 @@ def _estimates_with_input( f" v = inputs['{name}']\n" " return {'median_obs': v, 'ci95_lower': v * 0.6, 'ci95_upper': v * 1.4}" ), - "population_spread": "center_only", } @@ -1976,103 +2071,6 @@ def test_justified_mismatch_passes( CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) -class TestCalCenterChannelNotPopulation: - """Ported from SubmodelTarget.validate_center_channel_sem_scale. - - An across_patient target has two channels: `samples` is the population - spread (omega), median/ci95 pin the centre and must shrink with n. Returning - the population percentiles as ci95 encodes the spread twice. Measured - against the live corpus: 26 targets do exactly this, all via - population.summarize() without an n. - """ - - # Population draw whose median matches the golden reported median (1.0). - _POP = ( - "def derive_distribution(inputs, ureg):\n" - " import numpy as np, math\n" - " rng = np.random.default_rng(42)\n" - " mean = inputs['cd8_ratio_mean']\n" - " sigma_log = inputs['cd8_ratio_sigma_log']\n" - " mu_log = math.log(mean.magnitude)\n" - " samples = rng.lognormal(mu_log, sigma_log.magnitude, 10000) * mean.units\n" - "{body}" - ) - - _POPULATION_CI = _POP.format( - body=( - " ci95 = np.percentile(samples, [2.5, 97.5])\n" - " return {'median_obs': np.median(samples), 'ci95_lower': ci95[0],\n" - " 'ci95_upper': ci95[1], 'samples': samples}" - ) - ) - - _SEM_CI = _POP.format( - body=( - " from maple.core.calibration import population as pop\n" - " return pop.summarize(samples, n=42, rng=rng)" - ) - ) - - def _data(self, golden, code, **ed): - data = copy.deepcopy(golden) - data["empirical_data"]["distribution_code"] = code - data["empirical_data"]["population_spread"] = "across_patient" - data["empirical_data"].update(ed) - return data - - def test_population_range_as_ci95_rejected( - self, model_structure, golden_calibration_target_data, mock_crossref_success - ): - data = self._data(golden_calibration_target_data, self._POPULATION_CI) - with pytest.raises(ValidationError, match="is the POPULATION range"): - CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) - - def test_error_names_the_channel_that_needs_fixing( - self, model_structure, golden_calibration_target_data, mock_crossref_success - ): - """Clearing the error with population_spread='center_only' would delete - the omega channel. The message must steer to summarize(n=...) instead.""" - data = self._data(golden_calibration_target_data, self._POPULATION_CI) - with pytest.raises(ValidationError) as exc: - CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) - msg = str(exc.value) - assert "summarize(samples, n=" in msg - assert "DELETES the population channel" in msg - - def test_sem_scale_ci95_passes( - self, model_structure, golden_calibration_target_data, mock_crossref_success - ): - # summarize(n=42) bootstraps the median, so ci95 is far narrower than the - # population range; median and samples are unchanged. - data = self._data(golden_calibration_target_data, self._SEM_CI, ci95=[[0.888, 1.124]]) - target = CalibrationTarget.model_validate( - data, context={"model_structure": model_structure} - ) - lo, hi = target.empirical_data.ci95[0] - assert 0.8 < lo < 1.0 and 1.0 < hi < 1.3 - - def test_single_subject_exempt( - self, model_structure, golden_calibration_target_data, mock_crossref_success - ): - """With n=1 the centre's uncertainty IS the population spread; no double - encoding is possible, so the check must not fire.""" - data = self._data( - golden_calibration_target_data, - self._POPULATION_CI, - sample_size=1, - sample_size_rationale="single reported subject", - ) - CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) - - def test_center_only_targets_unaffected( - self, model_structure, golden_calibration_target_data, mock_crossref_success - ): - """A center_only target has no second channel, so a wide ci95 is fine.""" - CalibrationTarget.model_validate( - golden_calibration_target_data, context={"model_structure": model_structure} - ) - - class TestCalNoHardcodedValuesInDistributionCode: """Ported from SubmodelTarget.validate_no_hardcoded_values_in_observation_code, scoped to non-integer floats. @@ -2145,7 +2143,7 @@ def test_integers_are_exempt( corpus is a count, index or bound, so integers are not candidates.""" data = self._with_code( golden_calibration_target_data, - " n_mc = 200000\n" " vals = [inputs['cd8_ratio_mean'] for _ in range(1, 11)]\n", + " n_mc = 200000\n vals = [inputs['cd8_ratio_mean'] for _ in range(1, 11)]\n", ) CalibrationTarget.model_validate(data, context={"model_structure": model_structure}) @@ -2202,6 +2200,10 @@ def test_mechanistic_targets_exempt( """A mechanistic target encodes reasoning, not measured text.""" data = copy.deepcopy(golden_calibration_target_data) data["epistemic_basis"] = "mechanistic" + # No cohort to read n from once it stops being a measurement. + data.pop("cohort_id") + data["empirical_data"]["sample_size"] = 42 + data["empirical_data"]["sample_size_rationale"] = "assumed cohort scale" with patch( "maple.core.calibration.snippet_validator.load_paper_texts", return_value={"smith_2020": ("Nothing matching.", "pdf")}, diff --git a/tests/unit/core/test_cohort_registry.py b/tests/unit/core/test_cohort_registry.py new file mode 100644 index 0000000..a0b29a8 --- /dev/null +++ b/tests/unit/core/test_cohort_registry.py @@ -0,0 +1,794 @@ +"""Cohort registry, its cross-target checks, and the Observable measurement attributes.""" + +import warnings + +import pytest +import yaml +from pydantic import ValidationError + +from maple.core.calibration.cohort import ( + Cohort, + CohortRegistry, + PatientBlock, + Stratum, + load_cohorts, +) +from maple.core.calibration.enums import AssayModality, QuantityKind +from maple.core.calibration.observable import Observable +from maple.core.calibration.registry_audit import ( + check_registries, + covariance_blocks, + find_registry_problems, + resolve_n, + warn_merged_blocks, + warn_unused_cohorts, +) + +_CODE = ( + "def compute_observable(time, species_dict, constants):\n" + " return species_dict['V_T.CD8'] / species_dict['V_T.nucleated']\n" +) + + +def _cohort(**over): + base = dict( + cohort_id="li2022_arm_a", + description="Arm A patients, paired pre/post biopsy.", + scenarios=["baseline_no_treatment"], + n_c=9, + source_tag="Li2022_CancerCell_PDAC_AntiPD1", + ) + base.update(over) + return Cohort(**base) + + +def _registry(**over): + return CohortRegistry(cohorts=[_cohort(**over)]) + + +def _block(**over): + """Two cohorts of 9 and 10 sharing 6 patients, which is Li 2022 in miniature.""" + base = dict( + block_id="li2022", + description="One trial, two reported occasions.", + cohorts=["a", "b"], + strata=[ + Stratum(cohorts=["a", "b"], n=6), + Stratum(cohorts=["a"], n=3), + Stratum(cohorts=["b"], n=4), + ], + ) + base.update(over) + return PatientBlock(**base) + + +def _ab_registry(**over): + return CohortRegistry( + cohorts=[_cohort(cohort_id="a", n_c=9), _cohort(cohort_id="b", n_c=10)], + blocks=[_block(**over)], + ) + + +def _observable(**over): + readout = dict( + quantity_kind=QuantityKind.FRACTION, + assay_modality=AssayModality.MIHC, + numerator_species=["V_T.CD8"], + denominator_species=["V_T.nucleated"], + ) + readout.update(over.pop("readout", {})) + for key in ("quantity_kind", "assay_modality", "reference"): + if key in over: + readout[key] = over.pop(key) + base = dict( + readout=readout, + code=_CODE, + units="dimensionless", + species=["V_T.CD8", "V_T.nucleated"], + support="unit_interval", + readout_time=0.0, + readout_time_unit="day", + ) + base.update(over) + return Observable(**base) + + +def _target(**over): + base = dict( + cohort_id="li2022_arm_a", + epistemic_basis="literature", + observable={ + "code": _CODE, + "readout_time": 0.0, + "readout": { + "quantity_kind": "fraction", + "assay_modality": "mihc", + "numerator_species": ["V_T.CD8"], + "denominator_species": ["V_T.nucleated"], + }, + }, + empirical_data={"inputs": [{"source_ref": "Li2022_CancerCell_PDAC_AntiPD1"}]}, + primary_data_source={"source_tag": "Li2022_CancerCell_PDAC_AntiPD1"}, + ) + base.update(over) + return base + + +def _fraction(numerator, denominator, center=None, **over): + """A target whose code divides ``numerator`` by the sum of ``denominator``.""" + den = " + ".join(f"species_dict['{s}']" for s in denominator) + t = _target( + observable={ + "code": ( + "def compute_observable(time, species_dict, constants):\n" + f" return species_dict['{numerator}'] / ({den})\n" + ), + "readout_time": 0.0, + "readout": { + "quantity_kind": "fraction", + "assay_modality": "mihc", + "numerator_species": [numerator], + "denominator_species": list(denominator), + }, + }, + **over, + ) + if center is not None: + t["empirical_data"]["observed_distribution"] = { + "statistics": [{"stat": "quantile", "p": 0.5, "value": center}], + "spread_source": "across_patient", + } + return t + + +# --------------------------------------------------------------------------- # +# Cohort # +# --------------------------------------------------------------------------- # +class TestCohort: + def test_minimal_cohort_validates(self): + assert _cohort().n_c == 9 + + def test_no_scenarios_rejected(self): + with pytest.raises(ValueError, match="declares no scenarios"): + _cohort(scenarios=[]) + + def test_repeated_scenario_rejected(self): + with pytest.raises(ValueError, match="repeats a scenario"): + _cohort(scenarios=["a", "a"]) + + def test_overlap_is_not_a_cohort_field(self): + """Sharing is a property of a set of cohorts, so it lives on the block.""" + with pytest.raises(ValidationError): + _cohort(shares_patients_with=["arm_b"]) + + def test_zero_n_rejected(self): + with pytest.raises(ValueError): + _cohort(n_c=0) + + def test_source_tag_is_singular(self): + """A cohort is one study's patients, so pooling is unrepresentable.""" + with pytest.raises(ValueError): + _cohort(source_tag=["a", "b"]) + + def test_eligibility_needs_a_bound(self): + with pytest.raises(ValueError, match="neither lo nor hi"): + _cohort( + eligibility=[ + {"target_id": "cd8_fraction", "units": "dimensionless", "rationale": "x"} + ] + ) + + def test_eligibility_inverted_bounds_rejected(self): + with pytest.raises(ValueError, match="selects nobody"): + _cohort( + eligibility=[ + { + "target_id": "cd8_fraction", + "lo": 0.5, + "hi": 0.1, + "units": "dimensionless", + "rationale": "x", + } + ] + ) + + +class TestCohortRegistry: + def test_duplicate_ids_rejected(self): + with pytest.raises(ValueError, match="Duplicate cohort_id"): + CohortRegistry(cohorts=[_cohort(), _cohort()]) + + def test_block_members_must_resolve(self): + with pytest.raises(ValueError, match="not in the registry"): + CohortRegistry( + cohorts=[_cohort(cohort_id="a", n_c=9)], + blocks=[_block(cohorts=["a", "ghost"], strata=None)], + ) + + def test_counted_block_accepted(self): + assert len(_ab_registry().blocks) == 1 + + def test_load_from_yaml(self, tmp_path): + p = tmp_path / "cohorts.yaml" + p.write_text(yaml.safe_dump({"cohorts": [_cohort().model_dump()]})) + assert load_cohorts(p).get("li2022_arm_a").n_c == 9 + + def test_load_bare_list(self, tmp_path): + p = tmp_path / "cohorts.yaml" + p.write_text(yaml.safe_dump([_cohort().model_dump()])) + assert len(load_cohorts(p).cohorts) == 1 + + def test_missing_file_raises(self, tmp_path): + with pytest.raises(FileNotFoundError): + load_cohorts(tmp_path / "nope.yaml") + + +# --------------------------------------------------------------------------- # +# Readout # +# --------------------------------------------------------------------------- # +class TestReadoutComposition: + """The declared composition is the row's identity; the code is audited against it.""" + + def test_numerator_is_required(self): + with pytest.raises(ValidationError): + _observable(readout={"numerator_species": []}) + + def test_denominator_defaults_to_empty(self): + obs = _observable(readout={"denominator_species": []}) + assert obs.readout.denominator_species == [] + + def test_identical_numerator_and_denominator_rejected(self): + with pytest.raises(ValidationError, match="constant at 1"): + _observable( + readout={ + "numerator_species": ["V_T.CD8"], + "denominator_species": ["V_T.CD8"], + } + ) + + def test_experimental_denominator_needs_model_species(self): + with pytest.raises(ValidationError, match="denominator_species"): + _observable( + readout={ + "denominator_species": [], + "experimental_denominator": "all nucleated cells", + } + ) + + +class TestObservableMeasurementAttributes: + """The attributes inference reads to build the measurement-discrepancy design.""" + + def test_minimal_observable_validates(self): + assert _observable().readout.quantity_kind == QuantityKind.FRACTION + + def test_quantity_kind_is_required(self): + with pytest.raises(ValidationError): + _observable(quantity_kind=None) + + def test_assay_modality_is_required(self): + with pytest.raises(ValidationError): + _observable(assay_modality=None) + + def test_unregistered_quantity_kind_rejected(self): + with pytest.raises(ValidationError): + _observable(quantity_kind="luminosity") + + def test_unregistered_modality_rejected(self): + with pytest.raises(ValidationError): + _observable(assay_modality="vibes") + + def test_foldchange_requires_a_reference(self): + with pytest.raises(ValidationError, match="no reference is declared"): + _observable(quantity_kind=QuantityKind.FOLDCHANGE) + + def test_foldchange_with_a_timepoint_reference(self): + obs = _observable( + quantity_kind=QuantityKind.FOLDCHANGE, + reference={"kind": "timepoint", "timepoint": 0.0, "timepoint_unit": "day"}, + ) + assert obs.readout.reference.timepoint == 0.0 + + def test_foldchange_with_a_scenario_reference(self): + obs = _observable( + quantity_kind=QuantityKind.FOLDCHANGE, + reference={"kind": "scenario", "scenario": "baseline_no_treatment"}, + ) + assert obs.readout.reference.scenario == "baseline_no_treatment" + + def test_absolute_quantity_may_not_declare_a_reference(self): + with pytest.raises(ValidationError, match="absolute quantity"): + _observable(reference={"kind": "timepoint", "timepoint": 0.0, "timepoint_unit": "day"}) + + def test_timepoint_reference_needs_units(self): + with pytest.raises(ValidationError, match="needs timepoint and timepoint_unit"): + _observable( + quantity_kind=QuantityKind.FOLDCHANGE, + reference={"kind": "timepoint", "timepoint": 0.0}, + ) + + def test_scenario_reference_rejects_a_timepoint(self): + with pytest.raises(ValidationError, match="must not set timepoint"): + _observable( + quantity_kind=QuantityKind.FOLDCHANGE, + reference={"kind": "scenario", "scenario": "baseline", "timepoint": 1.0}, + ) + + def test_same_kind_and_modality_share_a_design_row(self): + """No identifier needed: the attributes are the row.""" + a = _observable() + b = _observable(readout={"numerator_species": ["V_T.Treg"]}) + assert (a.readout.quantity_kind, a.readout.assay_modality) == ( + b.readout.quantity_kind, + b.readout.assay_modality, + ) + + +# --------------------------------------------------------------------------- # +# Cross-target checks # +# --------------------------------------------------------------------------- # +class TestRegistryAudit: + def test_clean_corpus_passes(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + check_registries({"t1": _target()}, _registry()) + + def test_unknown_cohort_reported(self): + problems = find_registry_problems({"t1": _target(cohort_id="ghost")}, _registry()) + assert [p.kind for p in problems] == ["unknown_cohort"] + + def test_pooled_target_rejected(self): + """Several sources behind one target is a meta-analysis, not a cohort.""" + t = _target( + empirical_data={ + "inputs": [ + {"source_ref": "Golesworthy2022"}, + {"source_ref": "Liu2015"}, + {"source_ref": "Jansen2021"}, + ] + } + ) + pooled = [ + p for p in find_registry_problems({"t1": t}, _registry()) if p.kind == "pooled_target" + ] + assert len(pooled) == 1 + assert "Split into one target per source" in pooled[0].detail + + def test_mechanistic_target_exempt_from_pooling_and_cohort(self): + t = _target(cohort_id=None, epistemic_basis="mechanistic") + assert find_registry_problems({"t1": t}, _registry()) == [] + + def test_n_evaluable_may_not_exceed_cohort(self): + t = _target() + t["empirical_data"]["n_evaluable"] = 99 + kinds = {p.kind for p in find_registry_problems({"t1": t}, _registry())} + assert "n_evaluable_exceeds_cohort" in kinds + + def test_source_disagreeing_with_cohort_reported(self): + t = _target(primary_data_source={"source_tag": "Li2022_CancerCell"}) + kinds = {p.kind for p in find_registry_problems({"t1": t}, _registry())} + assert "source_disagrees_with_cohort" in kinds + + def test_duplicate_row_reported(self): + """One cohort reports a quantity once.""" + dup = [ + p + for p in find_registry_problems({"t1": _target(), "t2": _target()}, _registry()) + if p.kind == "duplicate_row" + ] + assert len(dup) == 1 + assert dup[0].target_ids == ("t1", "t2") + + def test_duplicate_row_keys_on_the_declared_composition(self): + """Different numerators over one denominator are different rows.""" + t2 = _target() + t2["observable"] = dict(t2["observable"]) + t2["observable"]["readout"] = dict( + t2["observable"]["readout"], numerator_species=["V_T.CD4"] + ) + problems = find_registry_problems({"t1": _target(), "t2": t2}, _registry()) + assert [p.kind for p in problems] == [] + + def test_same_expression_in_different_cohorts_is_fine(self): + cohorts = CohortRegistry(cohorts=[_cohort(), _cohort(cohort_id="other")]) + problems = find_registry_problems( + {"t1": _target(), "t2": _target(cohort_id="other")}, cohorts + ) + assert [p.kind for p in problems] == [] + + def test_check_raises_with_every_problem(self): + with pytest.raises(ValueError, match="registry problem"): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + check_registries({"t1": _target(cohort_id="ghost")}, _registry()) + + def test_unused_cohorts_warn(self): + cohorts = CohortRegistry(cohorts=[_cohort(), _cohort(cohort_id="unused")]) + with pytest.warns(UserWarning, match="no target uses"): + unused = warn_unused_cohorts({"t1": _target()}, cohorts) + assert unused == ["unused"] + + +class TestSingularBlocks: + """Rows of one cohort that are deterministic functions of each other.""" + + _CAF = ["V_T.iCAF", "V_T.myCAF"] + + def _kinds(self, targets, cohorts=None): + return [p.kind for p in find_registry_problems(targets, cohorts or _registry())] + + def test_numerators_partitioning_their_denominator_flagged(self): + targets = { + "icaf": _fraction("V_T.iCAF", self._CAF), + "mycaf": _fraction("V_T.myCAF", self._CAF), + } + problems = find_registry_problems(targets, _registry()) + assert [p.kind for p in problems] == ["singular_block"] + assert problems[0].target_ids == ("icaf", "mycaf") + assert "partition" in problems[0].detail + + def test_numerators_not_exhausting_the_denominator_are_fine(self): + targets = { + "cd8": _fraction("V_T.CD8", ["V_T.nucleated"]), + "cd4": _fraction("V_T.CD4", ["V_T.nucleated"]), + } + assert self._kinds(targets) == [] + + def test_centers_summing_to_one_flagged(self): + """The model observables need not be complementary for the data to be.""" + den = self._CAF + ["V_T.apCAF"] + targets = { + "icaf": _fraction("V_T.iCAF", den, center=0.284), + "mycaf": _fraction("V_T.myCAF", den, center=0.716), + } + problems = find_registry_problems(targets, _registry()) + assert [p.kind for p in problems] == ["singular_block"] + assert "one number reported twice" in problems[0].detail + + def test_centers_summing_to_one_hundred_flagged(self): + den = self._CAF + ["V_T.apCAF"] + targets = { + "icaf": _fraction("V_T.iCAF", den, center=28.4), + "mycaf": _fraction("V_T.myCAF", den, center=71.6), + } + assert self._kinds(targets) == ["singular_block"] + + def test_centers_falling_short_of_one_are_fine(self): + den = self._CAF + ["V_T.apCAF"] + targets = { + "icaf": _fraction("V_T.iCAF", den, center=0.28), + "mycaf": _fraction("V_T.myCAF", den, center=0.61), + } + assert self._kinds(targets) == [] + + def test_centers_in_different_cohorts_are_fine(self): + cohorts = CohortRegistry(cohorts=[_cohort(), _cohort(cohort_id="other")]) + den = self._CAF + ["V_T.apCAF"] + targets = { + "icaf": _fraction("V_T.iCAF", den, center=0.284), + "mycaf": _fraction("V_T.myCAF", den, center=0.716, cohort_id="other"), + } + assert self._kinds(targets, cohorts) == [] + + def test_center_falls_back_to_the_computed_median(self): + den = self._CAF + ["V_T.apCAF"] + a = _fraction("V_T.iCAF", den) + b = _fraction("V_T.myCAF", den) + a["empirical_data"]["median"] = [0.284] + b["empirical_data"]["median"] = [0.716] + assert self._kinds({"icaf": a, "mycaf": b}) == ["singular_block"] + + def test_a_non_fraction_pair_is_not_checked(self): + """Two densities summing to 1 in their own units is a coincidence, not a constraint.""" + code = "def compute_observable(time, species_dict, constants):\n return species_dict['V_T.CD8']\n" + targets = {} + for tid, center in (("a", 0.284), ("b", 0.716)): + t = _target(observable={"code": code, "readout_time": 0.0}) + t["empirical_data"]["median"] = [center] + targets[tid] = t + assert self._kinds(targets) == [] + + +def _arm(role, scenario, cohort_id, numerator=("V_T.CD8_TLA",), **over): + base = dict( + role=role, + scenario=scenario, + cohort_id=cohort_id, + required_species=list(numerator), + observable_code="def compute_test_statistic(t, s):\n return s['V_T.CD8_TLA']\n", + readout={ + "quantity_kind": "density", + "assay_modality": "mihc", + "numerator_species": list(numerator), + }, + ) + base.update(over) + return base + + +def _contrast(**over): + """A two-arm cross-scenario target, as parsed YAML.""" + base = dict( + epistemic_basis="literature", + observable={ + "inputs": [ + _arm("nivo", "gvax_nivo_neoadjuvant", "arm_b"), + _arm("urelumab", "gvax_nivo_urelumab_neoadjuvant", "arm_c"), + ] + }, + empirical_data={"inputs": []}, + ) + base.update(over) + return base + + +def _two_arms(**over): + b = _cohort(cohort_id="arm_b", scenarios=["gvax_nivo_neoadjuvant"], n_c=10) + c = _cohort(cohort_id="arm_c", scenarios=["gvax_nivo_urelumab_neoadjuvant"], n_c=8, **over) + return CohortRegistry(cohorts=[b, c]) + + +class TestCrossScenarioArms: + """A contrast over disjoint arms is a derived row over several cohorts.""" + + def _kinds(self, targets, cohorts): + return [p.kind for p in find_registry_problems(targets, cohorts)] + + def test_placed_arms_pass(self): + assert self._kinds({"inv": _contrast()}, _two_arms()) == [] + + def test_unknown_cohort_on_an_arm_reported(self): + t = _contrast() + t["observable"]["inputs"][1]["cohort_id"] = "ghost" + assert self._kinds({"inv": t}, _two_arms()) == ["unknown_cohort"] + + def test_arm_scenario_must_belong_to_its_cohort(self): + t = _contrast() + t["observable"]["inputs"][1]["scenario"] = "baseline_no_treatment" + assert self._kinds({"inv": t}, _two_arms()) == ["scenario_not_in_cohort"] + + def test_arm_n_evaluable_may_not_exceed_its_cohort(self): + t = _contrast() + t["observable"]["inputs"][1]["n_evaluable"] = 99 + assert self._kinds({"inv": t}, _two_arms()) == ["n_evaluable_exceeds_cohort"] + + def test_arms_sharing_patients_are_a_paired_contrast(self): + b = _cohort(cohort_id="arm_b", scenarios=["gvax_nivo_neoadjuvant"], n_c=10) + c = _cohort(cohort_id="arm_c", scenarios=["gvax_nivo_urelumab_neoadjuvant"], n_c=8) + reg = CohortRegistry( + cohorts=[b, c], + blocks=[ + PatientBlock( + block_id="trial", + description="One trial, two arms with crossover.", + cohorts=["arm_b", "arm_c"], + strata=[ + Stratum(cohorts=["arm_b", "arm_c"], n=4), + Stratum(cohorts=["arm_b"], n=6), + Stratum(cohorts=["arm_c"], n=4), + ], + ) + ], + ) + problems = find_registry_problems({"inv": _contrast()}, reg) + assert [p.kind for p in problems] == ["paired_contrast_as_cross_scenario"] + assert "resampling each arm independently" in problems[0].detail + + def test_disjoint_arms_of_one_block_are_not_a_paired_contrast(self): + """Two arms of one trial share a block and no patients; zero overlap is a fact.""" + b = _cohort(cohort_id="arm_b", scenarios=["gvax_nivo_neoadjuvant"], n_c=10) + c = _cohort(cohort_id="arm_c", scenarios=["gvax_nivo_urelumab_neoadjuvant"], n_c=8) + d = _cohort(cohort_id="baseline", n_c=18) + reg = CohortRegistry( + cohorts=[b, c, d], + blocks=[ + PatientBlock( + block_id="trial", + description="Both arms biopsied at baseline, neither crosses over.", + cohorts=["arm_b", "arm_c", "baseline"], + strata=[ + Stratum(cohorts=["arm_b", "baseline"], n=10), + Stratum(cohorts=["arm_c", "baseline"], n=8), + ], + ) + ], + ) + kinds = [p.kind for p in find_registry_problems({"inv": _contrast()}, reg)] + assert "paired_contrast_as_cross_scenario" not in kinds + + def test_arm_duplicating_a_standalone_target_reported(self): + """Conditioning on the constituent and the contrast counts one number twice.""" + standalone = _target(cohort_id="arm_b") + standalone["observable"] = dict( + standalone["observable"], + readout={ + "quantity_kind": "density", + "assay_modality": "mihc", + "numerator_species": ["V_T.CD8_TLA"], + }, + ) + kinds = self._kinds({"inv": _contrast(), "cd8_tla": standalone}, _two_arms()) + assert "redundant_cross_scenario_arm" in kinds + + def test_mechanistic_contrast_needs_no_cohorts(self): + t = _contrast(epistemic_basis="mechanistic") + for arm in t["observable"]["inputs"]: + arm.pop("cohort_id") + assert self._kinds({"inv": t}, _two_arms()) == [] + + +class TestCovarianceBlocks: + def test_independent_cohorts_are_their_own_blocks(self): + assert covariance_blocks({}, _two_arms()) == [frozenset({"arm_b"}), frozenset({"arm_c"})] + + def test_a_contrast_merges_the_cohorts_it_draws_on(self): + assert covariance_blocks({"inv": _contrast()}, _two_arms()) == [ + frozenset({"arm_b", "arm_c"}) + ] + + def test_shared_patients_merge_without_any_target(self): + reg = CohortRegistry( + cohorts=[ + _cohort(cohort_id="a", n_c=9), + _cohort(cohort_id="b", n_c=10), + _cohort(cohort_id="c"), + ], + blocks=[_block()], + ) + assert covariance_blocks({}, reg) == [frozenset({"a", "b"}), frozenset({"c"})] + + def test_an_uncounted_block_merges_the_same_way(self): + """It cannot say how much they overlap, only that they do.""" + reg = CohortRegistry( + cohorts=[_cohort(cohort_id="a", n_c=9), _cohort(cohort_id="b", n_c=10)], + blocks=[_block(strata=None)], + ) + assert covariance_blocks({}, reg) == [frozenset({"a", "b"})] + + def test_both_edge_kinds_compose_into_one_component(self): + b = _cohort(cohort_id="arm_b", scenarios=["gvax_nivo_neoadjuvant"], n_c=10) + c = _cohort(cohort_id="arm_c", scenarios=["gvax_nivo_urelumab_neoadjuvant"], n_c=8) + d = _cohort(cohort_id="arm_d", n_c=9) + reg = CohortRegistry( + cohorts=[b, c, d], + blocks=[_block(block_id="shared", cohorts=["arm_c", "arm_d"], strata=None)], + ) + blocks = covariance_blocks({"inv": _contrast()}, reg) + assert blocks == [frozenset({"arm_b", "arm_c", "arm_d"})] + + def test_a_scalar_target_merges_nothing(self): + cohorts = CohortRegistry(cohorts=[_cohort(), _cohort(cohort_id="other")]) + blocks = covariance_blocks({"t1": _target()}, cohorts) + assert blocks == [frozenset({"li2022_arm_a"}), frozenset({"other"})] + + def test_merged_blocks_warn(self): + with pytest.warns(UserWarning, match="one covariance block"): + merged = warn_merged_blocks({"inv": _contrast()}, _two_arms()) + assert merged == [frozenset({"arm_b", "arm_c"})] + + def test_independent_blocks_do_not_warn(self): + with warnings.catch_warnings(): + warnings.simplefilter("error") + assert warn_merged_blocks({}, _two_arms()) == [] + + +class TestResolveN: + def test_falls_back_to_cohort_n(self): + assert resolve_n(_target(), _registry()) == 9 + + def test_n_evaluable_wins(self): + t = _target() + t["empirical_data"]["n_evaluable"] = 6 + assert resolve_n(t, _registry()) == 6 + + def test_unknown_cohort_gives_none(self): + assert resolve_n(_target(cohort_id="ghost"), _registry()) is None + + +class TestPatientBlock: + def test_strata_count_the_overlap(self): + b = _block() + assert b.overlap("a", "b") == 6 + assert b.size_of("a") == 9 + assert b.size_of("b") == 10 + assert b.n_patients == 13 + + def test_zero_overlap_is_a_fact_not_a_gap(self): + """Cohorts in one block that share nobody, as two arms of one trial do.""" + b = _block( + cohorts=["a", "b", "c"], + strata=[ + Stratum(cohorts=["a", "c"], n=9), + Stratum(cohorts=["b", "c"], n=10), + ], + ) + assert b.overlap("a", "b") == 0 + assert b.n_patients == 19 + + def test_a_block_of_one_cohort_is_rejected(self): + with pytest.raises(ValidationError): + _block(cohorts=["a"], strata=[Stratum(cohorts=["a"], n=9)]) + + def test_repeated_stratum_rejected(self): + with pytest.raises(ValueError, match="repeats a stratum"): + _block( + strata=[ + Stratum(cohorts=["a", "b"], n=6), + Stratum(cohorts=["b", "a"], n=1), + ] + ) + + def test_stratum_naming_a_non_member_rejected(self): + with pytest.raises(ValueError, match="not in its cohorts"): + _block(strata=[Stratum(cohorts=["a", "b"], n=6), Stratum(cohorts=["z"], n=1)]) + + def test_member_with_no_patients_rejected(self): + with pytest.raises(ValueError, match="no stratum places them"): + _block(cohorts=["a", "b", "c"]) + + def test_empty_strata_list_rejected(self): + """Omitting the key declares an uncounted overlap; an empty list says nothing.""" + with pytest.raises(ValueError, match="empty strata list"): + _block(strata=[]) + + def test_zero_count_stratum_rejected(self): + with pytest.raises(ValidationError): + Stratum(cohorts=["a"], n=0) + + +class TestBlocksInTheRegistry: + def test_strata_must_add_up_to_each_n_c(self): + with pytest.raises(ValueError, match="declares n_c=10"): + CohortRegistry( + cohorts=[_cohort(cohort_id="a", n_c=9), _cohort(cohort_id="b", n_c=10)], + blocks=[ + _block(strata=[Stratum(cohorts=["a", "b"], n=6), Stratum(cohorts=["a"], n=3)]) + ], + ) + + def test_a_cohort_is_counted_by_at_most_one_block(self): + with pytest.raises(ValueError, match="Its patients can be divided up once"): + CohortRegistry( + cohorts=[_cohort(cohort_id="a", n_c=9), _cohort(cohort_id="b", n_c=10)], + blocks=[_block(), _block(block_id="other")], + ) + + def test_an_uncounted_block_may_overlay_a_counted_one(self): + """Li's internal split is reported; its overlap with a second paper is not.""" + reg = CohortRegistry( + cohorts=[ + _cohort(cohort_id="a", n_c=9), + _cohort(cohort_id="b", n_c=10), + _cohort(cohort_id="other_paper", n_c=7), + ], + blocks=[ + _block(), + _block( + block_id="one_trial", + cohorts=["a", "b", "other_paper"], + strata=None, + ), + ], + ) + assert reg.counted_block_for("a").block_id == "li2022" + assert [b.block_id for b in reg.uncounted_blocks] == ["one_trial"] + assert [b.block_id for b in reg.uncounted_blocks_for("other_paper")] == ["one_trial"] + assert reg.counted_block_for("other_paper") is None + + def test_an_uncounted_block_reports_no_counts(self): + b = _block(strata=None) + assert not b.is_quantified + assert b.n_patients is None + assert b.overlap("a", "b") is None + assert b.size_of("a") is None + + def test_duplicate_block_ids_rejected(self): + with pytest.raises(ValueError, match="Duplicate block_id"): + CohortRegistry( + cohorts=[_cohort(cohort_id="a", n_c=9), _cohort(cohort_id="b", n_c=10)], + blocks=[_block(strata=None), _block(strata=None)], + ) + + def test_cohorts_without_a_block_are_untouched(self): + assert _registry().counted_block_for("li2022_arm_a") is None + assert _registry().uncounted_blocks == [] diff --git a/tests/unit/core/test_cross_scenario_target.py b/tests/unit/core/test_cross_scenario_target.py index b26e8b3..f7f6dbe 100644 --- a/tests/unit/core/test_cross_scenario_target.py +++ b/tests/unit/core/test_cross_scenario_target.py @@ -19,7 +19,6 @@ CrossScenarioObservable, ) - # ============================================================================ # Fixtures # ============================================================================ @@ -34,7 +33,7 @@ def _scalar_empirical(): "units": "dimensionless", "sample_size": 1, "sample_size_rationale": ( - "Mechanistic prior; sample_size=1 denotes a single soft-prior " "assertion." + "Mechanistic prior; sample_size=1 denotes a single soft-prior assertion." ), "inputs": [], "assumptions": [], @@ -76,20 +75,31 @@ def invariance_target(): return { "cross_scenario_target_id": "cd8_intla_number_invariance_nivo_vs_urelumab", "observable": { - "code": ("def compute(inputs):\n" " return inputs['urelumab'] / inputs['nivo']\n"), + "code": ("def compute(inputs):\n return inputs['urelumab'] / inputs['nivo']\n"), "units": "dimensionless", + "quantity_kind": "ratio", "inputs": [ { "role": "nivo", "scenario": "gvax_nivo_neoadjuvant_zheng2022", "observable_code": code, "required_species": ["V_T.CD8_TLA", "V_T.CD8_TLA_act"], + "readout": { + "quantity_kind": "density", + "assay_modality": "mihc", + "numerator_species": ["V_T.CD8_TLA", "V_T.CD8_TLA_act"], + }, }, { "role": "urelumab", "scenario": "gvax_nivo_urelumab_neoadjuvant_heumann2023", "observable_code": code, "required_species": ["V_T.CD8_TLA", "V_T.CD8_TLA_act"], + "readout": { + "quantity_kind": "density", + "assay_modality": "mihc", + "numerator_species": ["V_T.CD8_TLA", "V_T.CD8_TLA_act"], + }, }, ], }, @@ -191,6 +201,11 @@ def _one_input(): "scenario": "s", "observable_code": "def compute_test_statistic(time, species_dict): return 0.0", "required_species": ["V_T.C1"], + "readout": { + "quantity_kind": "density", + "assay_modality": "mihc", + "numerator_species": ["V_T.C1"], + }, } @@ -200,6 +215,7 @@ def test_observable_requires_at_least_two_inputs(): CrossScenarioObservable( code="def compute(inputs): return inputs['only']", units="dimensionless", + quantity_kind="ratio", inputs=[_one_input()], ) @@ -213,6 +229,7 @@ def test_observable_rejects_duplicate_roles(): CrossScenarioObservable( code="def compute(inputs): return inputs['a']", units="dimensionless", + quantity_kind="ratio", inputs=[a, b], ) diff --git a/tests/unit/core/test_denominator_audit.py b/tests/unit/core/test_denominator_audit.py new file mode 100644 index 0000000..ba9b2bd --- /dev/null +++ b/tests/unit/core/test_denominator_audit.py @@ -0,0 +1,96 @@ +"""Cross-target denominator checks: mapping collisions, and code against declaration.""" + +import warnings + +from maple.core.calibration.denominator_audit import ( + find_code_readout_mismatches, + find_mapping_collisions, + numerator_and_denominator, + warn_code_readout_mismatches, +) + +_FRACTION = ( + "def compute_observable(time, species_dict, constants):\n" + " return species_dict['V_T.CD8'] / species_dict['V_T.nucleated']\n" +) + + +def _target(numerator=("V_T.CD8",), denominator=("V_T.nucleated",), code=_FRACTION, **over): + observable = { + "code": code, + "readout_time": 0.0, + "readout": { + "quantity_kind": "fraction", + "assay_modality": "mihc", + "numerator_species": list(numerator), + "denominator_species": list(denominator), + }, + } + observable.update(over) + return {"observable": observable} + + +class TestMappingCollisions: + def test_same_declared_composition_collides(self): + collisions = find_mapping_collisions({"a": _target(), "b": _target()}) + assert len(collisions) == 1 + assert collisions[0].members == ("a", "b") + + def test_different_numerator_does_not_collide(self): + assert find_mapping_collisions({"a": _target(), "b": _target(numerator=["V_T.CD4"])}) == [] + + def test_absolute_quantity_is_skipped(self): + """No denominator, no shared model quantity to collide over.""" + assert ( + find_mapping_collisions({"a": _target(denominator=[]), "b": _target(denominator=[])}) + == [] + ) + + def test_collision_keyed_on_declaration_not_code(self): + """Two spellings of one row still collide.""" + other_spelling = ( + "def compute_observable(time, species_dict, constants):\n" + " total = species_dict['V_T.nucleated']\n" + " return species_dict['V_T.CD8'] / total\n" + ) + collisions = find_mapping_collisions({"a": _target(), "b": _target(code=other_spelling)}) + assert len(collisions) == 1 + + +class TestCodeReadoutAgreement: + def test_agreement_is_silent(self): + assert find_code_readout_mismatches({"a": _target()}) == [] + + def test_disagreement_reported(self): + aggregated = ( + "def compute_observable(time, species_dict, constants):\n" + " return species_dict['V_T.CD8'] / species_dict['CD8_total_T']\n" + ) + found = find_code_readout_mismatches({"a": _target(code=aggregated)}) + assert len(found) == 1 + assert found[0].declared == ("V_T.nucleated",) + assert found[0].in_code == ("CD8_total_T",) + + def test_non_dividing_code_is_not_compared(self): + """An observable can reach its denominator without a division.""" + code = ( + "def compute_observable(time, species_dict, constants):\n" + " return species_dict['V_T.CD8'] * constants['area']\n" + ) + assert find_code_readout_mismatches({"a": _target(code=code)}) == [] + + def test_undeclared_denominator_is_not_compared(self): + assert find_code_readout_mismatches({"a": _target(denominator=[])}) == [] + + def test_mismatch_warns(self): + aggregated = _FRACTION.replace("V_T.nucleated", "CD8_total_T") + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + found = warn_code_readout_mismatches({"a": _target(code=aggregated)}) + assert len(found) == 1 + assert "divides by" in str(caught[0].message) + + +def test_parser_still_reads_a_division(): + num, den = numerator_and_denominator(_FRACTION) + assert (num, den) == ({"V_T.CD8"}, {"V_T.nucleated"}) diff --git a/tests/unit/core/test_observed_distribution.py b/tests/unit/core/test_observed_distribution.py index f80e6d8..ac0d5d7 100644 --- a/tests/unit/core/test_observed_distribution.py +++ b/tests/unit/core/test_observed_distribution.py @@ -1,34 +1,29 @@ #!/usr/bin/env python3 -""" -Tests for the shared quantile-anchor variability layer. - -Covers: -- QuantileAnchor / ObservedDistribution validators (probability range, non-crossing - quantile function, population-spread requires a scale, technical units cannot be a - population spread) -- ObservedDistribution derivations (median, quantile interpolation, IQR) -- SpreadSource / POPULATION_SPREAD_SOURCES routing semantics -- Wiring into CalibrationTargetEstimates (resolved_spread_source fallback + consistency) - and ErrorModel (additive, backwards-compatible) +"""Tests for the reported-statistics variability layer. + +Covers ReportedStatistic / ObservedDistribution validation, the on-demand +quantile derivations, unit provenance, SpreadSource routing, and the wiring into +CalibrationTargetEstimates and ErrorModel. """ import pytest from pydantic import ValidationError from maple.core.calibration.shared_models import ( + QuantileConvention, + DistributionShape, ExperimentalUnitType, - MomentSpread, ObservedDistribution, POPULATION_SPREAD_SOURCES, - QuantileAnchor, + ReportedStatistic, SourceRelevanceAssessment, SpreadSource, + StatKind, ) from maple.core.calibration.enums import HeterogeneityTransfer from maple.core.calibration.calibration_target_models import CalibrationTargetEstimates from maple.core.calibration.submodel_target import ErrorModel - _SOURCE_RELEVANCE = dict( indication_match="exact", indication_match_justification="exact PDAC match", @@ -40,31 +35,42 @@ tme_compatibility_notes="recapitulates target biology", ) - -# --------------------------------------------------------------------------- -# Fixtures / helpers -# --------------------------------------------------------------------------- - _DEFAULT_OBS_CODE = ( "def derive_observation(inputs, sample_size, rng, n_bootstrap):\n" " return rng.normal(0.0, 1.0, n_bootstrap)" ) +_Z_Q = 0.6744897501960817 +_Z_95 = 1.959963984540054 + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + + +def _q(p, value): + return ReportedStatistic(stat=StatKind.QUANTILE, p=p, value=value) + + +def _s(stat, value): + return ReportedStatistic(stat=stat, value=value) + def _median_iqr(q25, q50, q75, **kwargs): - # Inject the biological provenance a population spread now requires, so tests - # that aren't specifically about that requirement stay focused. + """A median with quartiles. Injects unit provenance for population sources.""" if kwargs.get("spread_source") in POPULATION_SPREAD_SOURCES: kwargs.setdefault("n_biological", 42) kwargs.setdefault("experimental_unit_type", ExperimentalUnitType.BIOLOGICAL) - return ObservedDistribution( - quantiles=[ - QuantileAnchor(p=0.25, value=q25), - QuantileAnchor(p=0.5, value=q50), - QuantileAnchor(p=0.75, value=q75), - ], - **kwargs, - ) + return ObservedDistribution(statistics=[_q(0.25, q25), _q(0.5, q50), _q(0.75, q75)], **kwargs) + + +def _center_scale(center, scale, shape, stat=StatKind.SD, center_stat="mean", **kwargs): + stats = [ + _q(0.5, center) if center_stat == "median" else _s(StatKind.MEAN, center), + _s(stat, scale), + ] + return ObservedDistribution(statistics=stats, shape=shape, **kwargs) def _cal_estimates(**overrides): @@ -82,262 +88,321 @@ def _cal_estimates(**overrides): # --------------------------------------------------------------------------- -# QuantileAnchor +# ReportedStatistic # --------------------------------------------------------------------------- +def test_quantile_needs_p(): + with pytest.raises(ValidationError, match="needs a probability level"): + ReportedStatistic(stat=StatKind.QUANTILE, value=1.0) + + @pytest.mark.parametrize("bad_p", [0.0, 1.0, -0.1, 1.5]) -def test_quantile_anchor_p_must_be_open_unit_interval(bad_p): - with pytest.raises(ValidationError): - QuantileAnchor(p=bad_p, value=1.0) +def test_quantile_p_must_be_in_open_unit_interval(bad_p): + with pytest.raises(ValidationError, match="p must be in"): + _q(bad_p, 1.0) + + +def test_non_quantile_must_not_set_p(): + with pytest.raises(ValidationError, match="must not set p"): + ReportedStatistic(stat=StatKind.SD, value=1.0, p=0.5) + + +def test_negative_width_rejected(): + with pytest.raises(ValidationError, match="cannot be negative"): + _s(StatKind.IQR, -1.0) -def test_quantile_anchor_valid(): - a = QuantileAnchor(p=0.5, value=3.0) - assert a.p == 0.5 and a.value == 3.0 +def test_negative_location_is_fine(): + """A mean may legitimately be negative; only widths may not.""" + assert _s(StatKind.MEAN, -3.0).value == -3.0 # --------------------------------------------------------------------------- -# ObservedDistribution validators +# ObservedDistribution validation # --------------------------------------------------------------------------- -def test_empty_quantiles_rejected(): - with pytest.raises(ValidationError): - ObservedDistribution(quantiles=[], spread_source=SpreadSource.CENTER_ONLY) +def test_empty_statistics_rejected(): + with pytest.raises(ValidationError, match="at least one entry"): + ObservedDistribution(statistics=[], spread_source=SpreadSource.CENTER_ONLY) -def test_duplicate_probability_levels_rejected(): - with pytest.raises(ValidationError, match="duplicate probability"): +def test_duplicate_statistic_rejected(): + with pytest.raises(ValidationError, match="repeats statistic"): ObservedDistribution( - quantiles=[QuantileAnchor(p=0.5, value=1.0), QuantileAnchor(p=0.5, value=2.0)], + statistics=[_s(StatKind.SD, 1.0), _s(StatKind.SD, 2.0)], spread_source=SpreadSource.CENTER_ONLY, ) +def test_duplicate_quantile_level_rejected(): + with pytest.raises(ValidationError, match="repeats statistic"): + ObservedDistribution( + statistics=[_q(0.5, 1.0), _q(0.5, 2.0)], spread_source=SpreadSource.CENTER_ONLY + ) + + +def test_same_stat_at_different_p_is_fine(): + d = ObservedDistribution( + statistics=[_q(0.25, 1.0), _q(0.75, 3.0)], spread_source=SpreadSource.CENTER_ONLY + ) + assert len(d.statistics) == 2 + + def test_crossing_quantile_function_rejected(): - # value must be non-decreasing in p with pytest.raises(ValidationError, match="non-decreasing"): ObservedDistribution( - quantiles=[QuantileAnchor(p=0.25, value=20.0), QuantileAnchor(p=0.75, value=10.0)], - spread_source=SpreadSource.ACROSS_PATIENT, + statistics=[_q(0.25, 20.0), _q(0.75, 10.0)], spread_source=SpreadSource.CENTER_ONLY ) -def test_population_spread_requires_a_scale(): - # single anchor cannot be a population spread (no width) - with pytest.raises(ValidationError, match="only a single quantile anchor"): +def test_population_spread_requires_a_width(): + with pytest.raises(ValidationError, match="no width statistic"): ObservedDistribution( - quantiles=[QuantileAnchor(p=0.5, value=1.0)], + statistics=[_q(0.5, 1.0)], spread_source=SpreadSource.ACROSS_PATIENT, + n_biological=10, + experimental_unit_type=ExperimentalUnitType.BIOLOGICAL, ) -def test_center_only_single_anchor_allowed(): +def test_scalar_width_satisfies_the_population_spread_requirement(): + """A median plus an SD is a width, even with no quartiles.""" d = ObservedDistribution( - quantiles=[QuantileAnchor(p=0.5, value=3.0)], spread_source=SpreadSource.CENTER_ONLY + statistics=[_q(0.5, 10.0), _s(StatKind.SD, 2.0)], + spread_source=SpreadSource.ACROSS_PATIENT, + n_biological=10, + experimental_unit_type=ExperimentalUnitType.BIOLOGICAL, ) - assert d.feeds_population_spread is False - assert d.iqr() is None + assert d.feeds_population_spread is True + + +def test_center_only_single_statistic_allowed(): + d = ObservedDistribution(statistics=[_q(0.5, 3.0)], spread_source=SpreadSource.CENTER_ONLY) + assert d.median() == 3.0 + + +def test_median_sd_and_se_together_are_representable(): + """The combination a single center+scale could not express.""" + d = ObservedDistribution( + statistics=[_q(0.5, 10.0), _s(StatKind.SD, 2.0), _s(StatKind.SE, 0.6)], + spread_source=SpreadSource.CENTER_ONLY, + ) + assert d.get(StatKind.SD) == 2.0 + assert d.get(StatKind.SE) == 0.6 + assert d.get(StatKind.QUANTILE, 0.5) == 10.0 + + +# --------------------------------------------------------------------------- +# Unit provenance +# --------------------------------------------------------------------------- @pytest.mark.parametrize("unit", [ExperimentalUnitType.TECHNICAL, ExperimentalUnitType.CLONAL]) def test_technical_unit_cannot_be_population_spread(unit): - with pytest.raises(ValidationError, match="not population variability"): - _median_iqr( - 10.0, - 15.0, - 20.0, - spread_source=SpreadSource.BIOLOGICAL_EXPERIMENTAL, - experimental_unit_type=unit, - ) + d = _median_iqr( + 10.0, + 15.0, + 25.0, + spread_source=SpreadSource.BIOLOGICAL_EXPERIMENTAL, + experimental_unit_type=unit, + ) + with pytest.raises(ValueError, match="not population variability"): + d.require_unit_provenance() -def test_anchors_stored_in_probability_order(): +def test_population_spread_requires_n_biological(): d = ObservedDistribution( - quantiles=[ - QuantileAnchor(p=0.75, value=20.0), - QuantileAnchor(p=0.25, value=10.0), - QuantileAnchor(p=0.5, value=15.0), - ], + statistics=[_q(0.25, 10.0), _q(0.75, 20.0)], spread_source=SpreadSource.ACROSS_PATIENT, - n_biological=30, experimental_unit_type=ExperimentalUnitType.BIOLOGICAL, ) - assert [q.p for q in d.quantiles] == [0.25, 0.5, 0.75] - - -def test_population_spread_requires_n_biological(): - with pytest.raises(ValidationError, match="n_biological is not set"): - ObservedDistribution( - quantiles=[QuantileAnchor(p=0.25, value=10.0), QuantileAnchor(p=0.75, value=20.0)], - spread_source=SpreadSource.BIOLOGICAL_EXPERIMENTAL, - experimental_unit_type=ExperimentalUnitType.BIOLOGICAL, - ) + with pytest.raises(ValueError, match="n_biological is not set"): + d.require_unit_provenance() def test_population_spread_requires_experimental_unit_type(): - with pytest.raises(ValidationError, match="experimental_unit_type is not set"): - ObservedDistribution( - quantiles=[QuantileAnchor(p=0.25, value=10.0), QuantileAnchor(p=0.75, value=20.0)], - spread_source=SpreadSource.ACROSS_PATIENT, - n_biological=30, - ) - - -def test_center_only_exempt_from_biological_provenance(): d = ObservedDistribution( - quantiles=[QuantileAnchor(p=0.5, value=3.0)], spread_source=SpreadSource.CENTER_ONLY + statistics=[_q(0.25, 10.0), _q(0.75, 20.0)], + spread_source=SpreadSource.ACROSS_PATIENT, + n_biological=8, ) - assert d.n_biological is None and d.experimental_unit_type is None + with pytest.raises(ValueError, match="experimental_unit_type is not set"): + d.require_unit_provenance() + + +def test_center_only_exempt_from_unit_provenance(): + d = ObservedDistribution(statistics=[_q(0.5, 3.0)], spread_source=SpreadSource.CENTER_ONLY) + d.require_unit_provenance() # does not raise # --------------------------------------------------------------------------- -# ObservedDistribution derivations +# Derivations # --------------------------------------------------------------------------- -def test_median_and_iqr(): - d = _median_iqr(10.0, 15.0, 25.0, spread_source=SpreadSource.ACROSS_PATIENT, n_biological=42) +def test_median_and_iqr_from_quantiles(): + d = _median_iqr(10.0, 15.0, 25.0, spread_source=SpreadSource.ACROSS_PATIENT) assert d.median() == 15.0 assert d.iqr() == 15.0 +def test_reported_iqr_wins_over_derivation(): + d = ObservedDistribution( + statistics=[_q(0.5, 10.0), _s(StatKind.IQR, 4.0)], + spread_source=SpreadSource.CENTER_ONLY, + shape=DistributionShape.NORMAL, + ) + assert d.iqr() == 4.0 + + def test_quantile_interpolation_and_clamping(): d = _median_iqr(10.0, 15.0, 25.0, spread_source=SpreadSource.ACROSS_PATIENT) - # midway between p=0.25 (10) and p=0.5 (15) at p=0.375 -> 12.5 assert d.quantile(0.375) == pytest.approx(12.5) - # clamp below/above the anchor range - assert d.quantile(0.01) == 10.0 - assert d.quantile(0.99) == 25.0 + assert d.quantile(0.01) == 10.0 # clamped to the lowest anchor + assert d.quantile(0.99) == 25.0 # clamped to the highest def test_iqr_none_when_range_not_spanned(): d = ObservedDistribution( - quantiles=[QuantileAnchor(p=0.4, value=1.0), QuantileAnchor(p=0.6, value=2.0)], - spread_source=SpreadSource.ACROSS_PATIENT, - n_biological=30, - experimental_unit_type=ExperimentalUnitType.BIOLOGICAL, + statistics=[_q(0.4, 1.0), _q(0.6, 2.0)], spread_source=SpreadSource.CENTER_ONLY ) assert d.iqr() is None +def test_center_prefers_median_over_mean(): + d = ObservedDistribution( + statistics=[_q(0.5, 3.0), _s(StatKind.MEAN, 5.0)], + spread_source=SpreadSource.CENTER_ONLY, + ) + assert d.center() == 3.0 + + # --------------------------------------------------------------------------- -# Moments form (mean +/- SD etc.) — native authoring, central expansion +# Expansion from a center and a scale # --------------------------------------------------------------------------- -def test_moments_lognormal_mean_sd_matches_hand_expansion(): - # k_C1_growth (Ahn 2017): VDT 132.3 +/- 132.1 days, lognormal, across 100 patients. - d = ObservedDistribution( - moments=MomentSpread(center=132.3, scale=132.1, scale_type="sd", shape="lognormal"), - spread_source=SpreadSource.ACROSS_PATIENT, - n_biological=100, - experimental_unit_type=ExperimentalUnitType.BIOLOGICAL, - ) - assert d.quantile(0.25) == pytest.approx(53.4, abs=0.2) - assert d.median() == pytest.approx(93.6, abs=0.2) - assert d.quantile(0.75) == pytest.approx(164.1, abs=0.2) +def test_normal_expansion_is_symmetric(): + d = _center_scale(10.0, 2.0, DistributionShape.NORMAL, spread_source=SpreadSource.CENTER_ONLY) + assert d.median() == pytest.approx(10.0) + assert d.quantile(0.75) - d.median() == pytest.approx(d.median() - d.quantile(0.25)) + assert d.iqr() == pytest.approx(2 * _Z_Q * 2.0) -def test_moments_normal_mean_sd(): - d = ObservedDistribution( - moments=MomentSpread(center=10.0, scale=2.0, scale_type="sd", shape="normal"), +def test_lognormal_expansion_is_log_symmetric(): + d = _center_scale( + 132.3, 132.1, DistributionShape.LOGNORMAL, spread_source=SpreadSource.CENTER_ONLY + ) + # q25 * q75 == median^2 for a lognormal. + assert d.quantile(0.25) * d.quantile(0.75) == pytest.approx(d.median() ** 2) + # A mean-anchored lognormal has its median below the mean. + assert d.median() < 132.3 + + +def test_lognormal_median_anchored_keeps_the_median(): + d = _center_scale( + 50.0, + 20.0, + DistributionShape.LOGNORMAL, + center_stat="median", spread_source=SpreadSource.CENTER_ONLY, ) - assert d.median() == 10.0 - # IQR = 2 * Phi^-1(0.75) * sd = 2 * 0.674489 * 2 - assert d.iqr() == pytest.approx(2.0 * 0.6744897501960817 * 2.0) + assert d.median() == pytest.approx(50.0) -def test_moments_sem_recovers_population_sd_with_n(): - # SD = SEM * sqrt(n); normal -> IQR = 2 * Z_Q * SD - d = ObservedDistribution( - moments=MomentSpread(center=1.97, scale=0.4, scale_type="sem", shape="normal"), - spread_source=SpreadSource.BIOLOGICAL_EXPERIMENTAL, - n_biological=10, - experimental_unit_type=ExperimentalUnitType.BIOLOGICAL, +def test_logit_normal_stays_in_the_unit_interval(): + d = _center_scale( + 0.95, + 0.10, + DistributionShape.LOGIT_NORMAL, + center_stat="median", + spread_source=SpreadSource.CENTER_ONLY, ) - expected = 2.0 * 0.6744897501960817 * 0.4 * (10**0.5) - assert d.iqr() == pytest.approx(expected) + assert 0.0 < d.quantile(0.25) < d.quantile(0.75) < 1.0 -def test_moments_cv_lognormal(): - d = ObservedDistribution( - moments=MomentSpread(center=100.0, scale=0.5, scale_type="cv", shape="lognormal"), +def test_logit_normal_center_must_be_in_the_unit_interval(): + d = _center_scale( + 1.5, + 0.1, + DistributionShape.LOGIT_NORMAL, + center_stat="median", spread_source=SpreadSource.CENTER_ONLY, ) - # median = mean / sqrt(1 + cv^2) - assert d.median() == pytest.approx(100.0 / (1.25**0.5)) + with pytest.raises(ValueError, match="center in"): + d.quantile(0.25) -def test_moments_iqr_normal_direct(): +def test_iqr_scale_expands_through_the_normal_equivalent_sd(): d = ObservedDistribution( - moments=MomentSpread(center=15.0, scale=10.0, scale_type="iqr", shape="normal"), + statistics=[_q(0.5, 15.0), _s(StatKind.IQR, 10.0)], + shape=DistributionShape.NORMAL, spread_source=SpreadSource.CENTER_ONLY, ) - assert d.iqr() == pytest.approx(10.0) - assert d.median() == 15.0 + assert d.quantile(0.75) - d.quantile(0.25) == pytest.approx(10.0) -def test_exactly_one_form_required(): - with pytest.raises(ValidationError, match="EXACTLY ONE"): - ObservedDistribution(spread_source=SpreadSource.CENTER_ONLY) - with pytest.raises(ValidationError, match="EXACTLY ONE"): - ObservedDistribution( - quantiles=[QuantileAnchor(p=0.5, value=1.0)], - moments=MomentSpread(center=1.0, scale=1.0, scale_type="sd", shape="normal"), - spread_source=SpreadSource.CENTER_ONLY, - ) +def test_cv_scale_uses_the_center(): + d = ObservedDistribution( + statistics=[_s(StatKind.MEAN, 100.0), _s(StatKind.CV, 0.5)], + shape=DistributionShape.NORMAL, + spread_source=SpreadSource.CENTER_ONLY, + ) + assert d.iqr() == pytest.approx(2 * _Z_Q * 50.0) -def test_moments_lognormal_median_iqr(): - # median (IQR) clinical form: recover q25/q75 that reproduce the median and IQR. +def test_ci95_bounds_give_a_scale(): d = ObservedDistribution( - moments=MomentSpread( - center=17.0, center_type="median", scale=21.0, scale_type="iqr", shape="lognormal" - ), + statistics=[ + _s(StatKind.MEAN, 10.0), + _s(StatKind.CI95_LO, 6.0), + _s(StatKind.CI95_HI, 14.0), + ], + shape=DistributionShape.NORMAL, spread_source=SpreadSource.CENTER_ONLY, ) - assert d.median() == pytest.approx(17.0) - assert d.iqr() == pytest.approx(21.0) - assert d.quantile(0.25) < 17.0 < d.quantile(0.75) + assert d.iqr() == pytest.approx(2 * _Z_Q * (8.0 / (2 * _Z_95))) -def test_moments_lognormal_mean_iqr_underdetermined_rejected(): - with pytest.raises(ValidationError, match="needs center_type='median'"): - ObservedDistribution( - moments=MomentSpread(center=5.0, scale=1.0, scale_type="iqr", shape="lognormal"), - spread_source=SpreadSource.CENTER_ONLY, - ) +def test_expansion_needs_a_shape(): + d = ObservedDistribution( + statistics=[_q(0.5, 10.0), _s(StatKind.SD, 2.0)], + spread_source=SpreadSource.CENTER_ONLY, + ) + with pytest.raises(ValueError, match="needs `shape`"): + d.quantile(0.25) -def test_moments_sem_without_n_rejected(): - with pytest.raises(ValidationError, match="needs n_biological"): - ObservedDistribution( - moments=MomentSpread(center=5.0, scale=1.0, scale_type="sem", shape="normal"), - spread_source=SpreadSource.CENTER_ONLY, - ) +def test_expansion_needs_a_width(): + d = ObservedDistribution( + statistics=[_q(0.5, 10.0)], + shape=DistributionShape.NORMAL, + spread_source=SpreadSource.CENTER_ONLY, + ) + with pytest.raises(ValueError, match="no width"): + d.quantile(0.25) -def test_moments_population_spread_requires_biological_provenance(): - # moments form is still subject to the biological-provenance rule - with pytest.raises(ValidationError, match="experimental_unit_type is not set"): - ObservedDistribution( - moments=MomentSpread(center=1.0, scale=0.2, scale_type="sd", shape="normal"), - spread_source=SpreadSource.ACROSS_PATIENT, - n_biological=20, - ) +def test_expansion_needs_a_center(): + d = ObservedDistribution( + statistics=[_s(StatKind.SD, 2.0)], + shape=DistributionShape.NORMAL, + spread_source=SpreadSource.CENTER_ONLY, + ) + with pytest.raises(ValueError, match="no center"): + d.quantile(0.25) -def test_moments_technical_unit_cannot_be_population_spread(): - with pytest.raises(ValidationError, match="not population variability"): - ObservedDistribution( - moments=MomentSpread(center=1.0, scale=0.2, scale_type="sd", shape="normal"), - spread_source=SpreadSource.BIOLOGICAL_EXPERIMENTAL, - n_biological=5, - experimental_unit_type=ExperimentalUnitType.TECHNICAL, - ) +def test_reported_quantiles_win_over_expansion(): + """Explicit quartiles are used as printed, not re-derived from the SD.""" + d = ObservedDistribution( + statistics=[_q(0.25, 1.0), _q(0.5, 2.0), _q(0.75, 9.0), _s(StatKind.SD, 0.1)], + shape=DistributionShape.NORMAL, + spread_source=SpreadSource.CENTER_ONLY, + ) + assert d.iqr() == pytest.approx(8.0) # --------------------------------------------------------------------------- @@ -346,11 +411,11 @@ def test_moments_technical_unit_cannot_be_population_spread(): def test_population_spread_sources_membership(): - assert SpreadSource.ACROSS_PATIENT in POPULATION_SPREAD_SOURCES - assert SpreadSource.BIOLOGICAL_EXPERIMENTAL in POPULATION_SPREAD_SOURCES + for s in (SpreadSource.ACROSS_PATIENT, SpreadSource.BIOLOGICAL_EXPERIMENTAL): + assert s in POPULATION_SPREAD_SOURCES for s in ( - SpreadSource.TECHNICAL, SpreadSource.CENTER_ONLY, + SpreadSource.TECHNICAL, SpreadSource.TRANSLATION, SpreadSource.ASSUMED, ): @@ -362,38 +427,53 @@ def test_population_spread_sources_membership(): # --------------------------------------------------------------------------- -def test_cal_resolved_spread_source_legacy_fallback(): - # No observed_distribution -> maps from legacy population_spread - assert _cal_estimates().resolved_spread_source == SpreadSource.CENTER_ONLY - assert ( - _cal_estimates(population_spread="across_patient").resolved_spread_source - == SpreadSource.ACROSS_PATIENT - ) +def test_cal_no_observed_distribution_is_center_only(): + e = _cal_estimates() + assert e.resolved_spread_source == SpreadSource.CENTER_ONLY + assert e.feeds_population_spread is False -def test_cal_resolved_spread_source_prefers_observed_distribution(): - od = _median_iqr(10.0, 15.0, 25.0, spread_source=SpreadSource.ACROSS_PATIENT, n_biological=42) - e = _cal_estimates(population_spread="across_patient", observed_distribution=od) +def test_cal_spread_source_comes_from_observed_distribution(): + # The cohort carries the unit accounting on the calibration side, so none is set here. + od = ObservedDistribution( + statistics=[_q(0.25, 10.0), _q(0.5, 15.0), _q(0.75, 25.0)], + spread_source=SpreadSource.ACROSS_PATIENT, + ) + e = _cal_estimates(observed_distribution=od) assert e.resolved_spread_source == SpreadSource.ACROSS_PATIENT + assert e.feeds_population_spread is True assert e.observed_distribution.median() == 15.0 -def test_cal_observed_distribution_contradiction_rejected(): - od = _median_iqr(10.0, 15.0, 25.0, spread_source=SpreadSource.ACROSS_PATIENT) - with pytest.raises(ValidationError, match="contradicts"): - _cal_estimates(population_spread="center_only", observed_distribution=od) - - -def test_cal_observed_distribution_center_only_consistent(): +def test_cal_center_only_observed_distribution_does_not_feed_omega(): + od = ObservedDistribution(statistics=[_q(0.5, 15.0)], spread_source=SpreadSource.CENTER_ONLY) + e = _cal_estimates(observed_distribution=od) + assert e.resolved_spread_source == SpreadSource.CENTER_ONLY + assert e.feeds_population_spread is False + + +@pytest.mark.parametrize( + "field,value", + [ + ("n_biological", 42), + ("n_technical", 3), + ("experimental_unit_type", ExperimentalUnitType.BIOLOGICAL), + ("unit_group", "donors"), + ("n_biological_is_floor", True), + ], +) +def test_cal_rejects_unit_accounting_the_cohort_owns(field, value): od = ObservedDistribution( - quantiles=[QuantileAnchor(p=0.5, value=15.0)], spread_source=SpreadSource.CENTER_ONLY + statistics=[_q(0.25, 10.0), _q(0.5, 15.0), _q(0.75, 25.0)], + spread_source=SpreadSource.ACROSS_PATIENT, + **{field: value}, ) - e = _cal_estimates(population_spread="center_only", observed_distribution=od) - assert e.resolved_spread_source == SpreadSource.CENTER_ONLY + with pytest.raises(ValidationError, match="cohort"): + _cal_estimates(observed_distribution=od) # --------------------------------------------------------------------------- -# ErrorModel (submodel) wiring — additive / backwards compatible +# ErrorModel (submodel) wiring # --------------------------------------------------------------------------- @@ -404,7 +484,7 @@ def test_submodel_observed_distribution_defaults_none(): def test_submodel_observed_distribution_set(): od = ObservedDistribution( - quantiles=[QuantileAnchor(p=0.25, value=1.0), QuantileAnchor(p=0.75, value=3.0)], + statistics=[_q(0.25, 1.0), _q(0.75, 3.0)], spread_source=SpreadSource.BIOLOGICAL_EXPERIMENTAL, n_biological=6, experimental_unit_type=ExperimentalUnitType.BIOLOGICAL, @@ -419,6 +499,26 @@ def test_submodel_observed_distribution_set(): assert em.observed_distribution.feeds_population_spread is True +def test_submodel_population_spread_must_state_its_units(): + """The requirement moved off ObservedDistribution onto the owner that carries n.""" + from maple.core.calibration.submodel_target import Calibration + + od = ObservedDistribution( + statistics=[_q(0.25, 1.0), _q(0.75, 3.0)], + spread_source=SpreadSource.BIOLOGICAL_EXPERIMENTAL, + ) + em = ErrorModel( + name="d1", + units="nM", + sample_size_input="n", + observation_code=_DEFAULT_OBS_CODE, + observed_distribution=od, + ) + cal = Calibration.model_construct(error_model=[em]) + with pytest.raises(ValueError, match="n_biological is not set"): + cal._spread_source_states_its_units() + + # --------------------------------------------------------------------------- # heterogeneity_transfer on SourceRelevanceAssessment # --------------------------------------------------------------------------- @@ -447,53 +547,40 @@ def test_heterogeneity_transfer_requires_justification(): # --------------------------------------------------------------------------- -# logit_normal shape (bounded [0, 1] fractions) +# logit_normal shape on a bounded observable # --------------------------------------------------------------------------- -def test_moments_logit_normal_stays_in_unit_interval(): - q25, q50, q75 = MomentSpread( - center=0.9, center_type="median", scale=0.15, scale_type="sd", shape="logit_normal" - ).to_quartiles() - assert 0.0 < q25 < q50 < q75 < 1.0 - assert q50 == pytest.approx(0.9) - - -def test_moments_logit_normal_near_one_does_not_escape(): - # lognormal would push the upper quartile past 1 here; logit_normal must not. - _, _, q75 = MomentSpread( - center=0.97, center_type="median", scale=0.05, scale_type="sd", shape="logit_normal" - ).to_quartiles() - assert q75 < 1.0 +def test_bounded_units_reject_a_normal_shape(): + od = ObservedDistribution( + statistics=[_q(0.5, 0.5), _s(StatKind.SD, 0.1)], + shape=DistributionShape.NORMAL, + spread_source=SpreadSource.CENTER_ONLY, + ) + with pytest.raises(ValidationError, match="logit_normal"): + _cal_estimates(units="percent", observed_distribution=od) -def test_moments_logit_normal_requires_median_center(): - with pytest.raises(ValidationError, match="center_type='median'"): - ObservedDistribution( - moments=MomentSpread( - center=0.5, center_type="mean", scale=0.1, scale_type="sd", shape="logit_normal" - ), - spread_source=SpreadSource.CENTER_ONLY, - ) +def test_bounded_units_accept_logit_normal(): + od = ObservedDistribution( + statistics=[_q(0.5, 0.5), _s(StatKind.SD, 0.1)], + shape=DistributionShape.LOGIT_NORMAL, + spread_source=SpreadSource.CENTER_ONLY, + ) + assert _cal_estimates(units="percent", observed_distribution=od) is not None -@pytest.mark.parametrize("bad_center", [1.5, 0.0, 1.0, -0.1]) -def test_moments_logit_normal_center_must_be_in_open_unit_interval(bad_center): - with pytest.raises(ValidationError, match=r"open interval \(0, 1\)"): - ObservedDistribution( - moments=MomentSpread( - center=bad_center, - center_type="median", - scale=0.1, - scale_type="sd", - shape="logit_normal", - ), - spread_source=SpreadSource.CENTER_ONLY, - ) +def test_no_declared_shape_is_exempt(): + """Reported quantiles are used as printed, so nothing is expanded.""" + od = ObservedDistribution( + statistics=[_q(0.25, 0.4), _q(0.5, 0.5), _q(0.75, 0.6)], + spread_source=SpreadSource.CENTER_ONLY, + ) + assert _cal_estimates(units="percent", observed_distribution=od) is not None # --------------------------------------------------------------------------- -# unit_group (shared-biological-unit panels) +# unit_group (shared-biological-unit panels, submodel side) # --------------------------------------------------------------------------- @@ -501,7 +588,8 @@ def _od( n_bio, group, unit=ExperimentalUnitType.BIOLOGICAL, src=SpreadSource.BIOLOGICAL_EXPERIMENTAL ): return ObservedDistribution( - moments=MomentSpread(center=100.0, scale=30.0, scale_type="sd", shape="normal"), + statistics=[_s(StatKind.MEAN, 100.0), _s(StatKind.SD, 30.0)], + shape=DistributionShape.NORMAL, spread_source=src, n_biological=n_bio, experimental_unit_type=unit, @@ -520,10 +608,10 @@ def _em(name, od): def _od_center(group): - # A center_only observed_distribution: does NOT feed omega, so it is exempt from - # the conditional-presence requirement even in a multi-entry target. + # Does NOT feed omega, so it is exempt from the conditional-presence rule. return ObservedDistribution( - moments=MomentSpread(center=100.0, scale=30.0, scale_type="sd", shape="normal"), + statistics=[_s(StatKind.MEAN, 100.0), _s(StatKind.SD, 30.0)], + shape=DistributionShape.NORMAL, spread_source=SpreadSource.CENTER_ONLY, unit_group=group, ) @@ -537,7 +625,8 @@ def _run_unit_group_check(ems): def test_unit_group_defaults_none(): od = ObservedDistribution( - moments=MomentSpread(center=1.0, scale=1.0, scale_type="sd", shape="normal"), + statistics=[_s(StatKind.MEAN, 1.0), _s(StatKind.SD, 1.0)], + shape=DistributionShape.NORMAL, spread_source=SpreadSource.CENTER_ONLY, ) assert od.unit_group is None @@ -548,19 +637,14 @@ def test_unit_group_consistent_members_ok(): def test_unit_group_unbalanced_n_allowed(): - # An unbalanced panel (per-point n differs — e.g. a donor missing at some doses) - # is still one population viewed several times; n_biological may differ. - _run_unit_group_check([_em("d1", _od(13, "donors")), _em("d2", _od(12, "donors"))]) + # A donor missing at some doses: one population viewed several times. + _run_unit_group_check([_em("d1", _od(13, "donors")), _em("d2", _od(11, "donors"))]) def test_unit_group_separate_groups_may_differ(): - # Distinct groups (e.g. two mouse lines) are independent — no consistency demand. _run_unit_group_check([_em("d1", _od(13, "lineA")), _em("d2", _od(20, "lineB"))]) -# --- conditional presence: required once 2+ population-spread observables exist --- - - def test_unit_group_required_when_multiple_population_entries(): with pytest.raises(ValueError, match="leave unit_group unset"): _run_unit_group_check([_em("d1", _od(13, None)), _em("d2", _od(13, None))]) @@ -571,17 +655,14 @@ def test_unit_group_multiple_population_tagged_same_ok(): def test_unit_group_multiple_population_tagged_different_ok(): - # Two separate populations, both tagged -> a conscious decision, allowed. _run_unit_group_check([_em("d1", _od(13, "lineA")), _em("d2", _od(20, "lineB"))]) def test_unit_group_single_population_entry_untagged_ok(): - # Only one population-spread observable -> no ambiguity, tag not required. _run_unit_group_check([_em("d1", _od(13, None)), _em("c1", _od_center(None))]) def test_unit_group_center_only_entries_exempt_from_presence(): - # Two entries, but neither feeds omega -> presence rule does not apply. _run_unit_group_check([_em("c1", _od_center(None)), _em("c2", _od_center(None))]) @@ -596,3 +677,86 @@ def test_unit_group_mismatched_spread_source_rejected(): ), ] ) + + +# --------------------------------------------------------------------------- +# A standard error is a width per sqrt(n) +# --------------------------------------------------------------------------- + + +def test_se_satisfies_the_population_width_requirement(): + """An SE determines a sample width given n, so it is not a center-only spread.""" + d = ObservedDistribution( + statistics=[_s(StatKind.MEAN, 9.84), _s(StatKind.SE, 1.16)], + spread_source=SpreadSource.ACROSS_PATIENT, + ) + assert d.feeds_population_spread + + +def test_se_widens_to_a_sample_sd_only_with_n(): + d = ObservedDistribution( + statistics=[_s(StatKind.MEAN, 9.84), _s(StatKind.SE, 1.16)], + spread_source=SpreadSource.ACROSS_PATIENT, + shape=DistributionShape.LOGNORMAL, + ) + assert d.population_sd() is None + assert d.population_sd(40) == pytest.approx(1.16 * 40**0.5) + + +def test_quantiles_from_an_se_need_n_and_say_so(): + d = ObservedDistribution( + statistics=[_s(StatKind.MEAN, 9.84), _s(StatKind.SE, 1.16)], + spread_source=SpreadSource.ACROSS_PATIENT, + shape=DistributionShape.LOGNORMAL, + ) + with pytest.raises(ValueError, match="pass n"): + d.median() + # With n the SE widens and the lognormal stays log-symmetric about its median. + assert d.quantile(0.25, 40) * d.quantile(0.75, 40) == pytest.approx(d.median(40) ** 2) + + +def test_a_reported_sd_is_preferred_over_an_se(): + """Both reported: the sample width is the printed one, not the widened SE.""" + d = ObservedDistribution( + statistics=[_s(StatKind.MEAN, 10.0), _s(StatKind.SD, 4.0), _s(StatKind.SE, 1.0)], + spread_source=SpreadSource.ACROSS_PATIENT, + ) + assert d.population_sd(100) == 4.0 + + +def test_negative_se_is_rejected(): + with pytest.raises(ValidationError, match="cannot be negative"): + ReportedStatistic(stat=StatKind.SE, value=-1.0) + + +def test_quantile_convention_defaults_to_unrecorded(): + """Papers do not state it, so the schema must not invent one.""" + d = _median_iqr(1.0, 2.0, 3.0, spread_source=SpreadSource.ACROSS_PATIENT) + assert d.quantile_convention is None + + +def test_quantile_convention_records_the_estimator(): + d = _median_iqr( + 1.0, + 2.0, + 3.0, + spread_source=SpreadSource.ACROSS_PATIENT, + quantile_convention=QuantileConvention.TYPE6, + ) + assert d.quantile_convention is QuantileConvention.TYPE6 + + +def test_quantile_convention_needs_a_quantile(): + with pytest.raises(ValidationError, match="no quantile was reported"): + ObservedDistribution( + statistics=[_s(StatKind.MEAN, 10.0), _s(StatKind.SD, 4.0)], + spread_source=SpreadSource.ACROSS_PATIENT, + quantile_convention=QuantileConvention.TYPE7, + ) + + +def test_quantile_convention_rejects_an_unknown_type(): + with pytest.raises(ValidationError): + _median_iqr( + 1.0, 2.0, 3.0, spread_source=SpreadSource.ACROSS_PATIENT, quantile_convention="type99" + ) diff --git a/tests/unit/core/test_submodel_target_validators.py b/tests/unit/core/test_submodel_target_validators.py index a6e21dd..4313b24 100644 --- a/tests/unit/core/test_submodel_target_validators.py +++ b/tests/unit/core/test_submodel_target_validators.py @@ -24,7 +24,6 @@ ) from maple.core.model_structure import ModelStructure, ModelParameter - # ============================================================================ # Default observation code — single source of truth for the current signature. # Tests that don't specifically test observation_code behavior should use this. @@ -1352,7 +1351,7 @@ def test_invisible_unicode_in_text_fails(self): input_value=10.0, ) # Add zero-width space in a text field - data["study_interpretation"] = "Test interpretation\u200Bwith invisible character" + data["study_interpretation"] = "Test interpretation\u200bwith invisible character" with pytest.raises(ValidationError) as exc_info: SubmodelTarget(**data) @@ -1985,10 +1984,10 @@ def derive_observation(inputs, sample_size, rng, n_bootstrap): def _population_observed_distribution() -> dict: return { - "quantiles": [ - {"p": 0.25, "value": 8.0}, - {"p": 0.5, "value": 10.0}, - {"p": 0.75, "value": 13.0}, + "statistics": [ + {"stat": "quantile", "p": 0.25, "value": 8.0}, + {"stat": "quantile", "p": 0.5, "value": 10.0}, + {"stat": "quantile", "p": 0.75, "value": 13.0}, ], "spread_source": "biological_experimental", "n_biological": 10, @@ -2071,12 +2070,12 @@ def test_well_behaved_positive_normal_no_warning(self): # ============================================================================ # Tests for validate_bounded_observable_uses_logit_normal (V-B) — bounded -# observable in the moments form must use shape='logit_normal'. +# observable declaring a shape must use shape='logit_normal'. # ============================================================================ def _percent_target(shape: str) -> dict: - """A percent-unit target with a moments-form population spread of the given shape.""" + """A percent-unit target whose population spread declares the given shape.""" data = make_algebraic_target( input_value=0.12, measurement_error_code=SEM_SCALE_OBSERVATION_CODE ) @@ -2086,13 +2085,11 @@ def _percent_target(shape: str) -> dict: em = data["calibration"]["error_model"][0] em["units"] = "percent" em["observed_distribution"] = { - "moments": { - "center": 0.12, - "center_type": "median", - "scale": 0.05, - "scale_type": "sd", - "shape": shape, - }, + "statistics": [ + {"stat": "quantile", "p": 0.5, "value": 0.12}, + {"stat": "sd", "value": 0.05}, + ], + "shape": shape, "spread_source": "biological_experimental", "n_biological": 10, "experimental_unit_type": "biological", @@ -2101,7 +2098,7 @@ def _percent_target(shape: str) -> dict: class TestBoundedObservableLogitNormal: - """V-B: bounded moments-form observable must use logit_normal.""" + """V-B: a bounded observable declaring a shape must use logit_normal.""" def test_percent_with_normal_shape_raises(self): with pytest.raises(ValidationError, match="logit_normal"): diff --git a/tests/unit/core/test_view_figure.py b/tests/unit/core/test_view_figure.py index 48cc078..def27c6 100644 --- a/tests/unit/core/test_view_figure.py +++ b/tests/unit/core/test_view_figure.py @@ -5,7 +5,6 @@ Tests HTML parsing, figure extraction, and fuzzy label matching. """ - from maple.core.tools.view_figure import ( _normalize_label, _label_matches, @@ -13,7 +12,6 @@ find_figure, ) - # ============================================================================ # Sample HTML fixtures # ============================================================================ diff --git a/tests/unit/test_test_stats_loader.py b/tests/unit/test_test_stats_loader.py index 7c62ed7..52bfdb0 100644 --- a/tests/unit/test_test_stats_loader.py +++ b/tests/unit/test_test_stats_loader.py @@ -192,6 +192,69 @@ def mixed_files_dir(temp_dir): # ============================================================================ +class TestSampleSizeFromTheCohort: + """A literature target carries no n; it resolves through the registry.""" + + def _placed_dir(self, temp_dir, **over): + import copy + + data = copy.deepcopy(BASELINE_YAML) + data["empirical_data"].pop("sample_size") + data["cohort_id"] = "smith2020" + data.update(over) + yaml_dir = temp_dir / "placed" + yaml_dir.mkdir() + with open(yaml_dir / "m1_m2_ratio.yaml", "w") as f: + yaml.dump(data, f) + return yaml_dir + + def _registry(self, temp_dir, n_c=113): + path = temp_dir / "cohorts.yaml" + with open(path, "w") as f: + yaml.dump( + { + "cohorts": [ + { + "cohort_id": "smith2020", + "description": "Resected PDAC.", + "scenarios": ["baseline_no_treatment"], + "n_c": n_c, + "source_tag": "Smith2020", + } + ] + }, + f, + ) + return path + + def test_resolves_from_the_cohort(self, temp_dir): + df = load_calibration_targets(self._placed_dir(temp_dir), cohorts=self._registry(temp_dir)) + assert df.iloc[0]["sample_size"] == 113 + + def test_n_evaluable_wins_over_the_cohort(self, temp_dir): + yaml_dir = self._placed_dir(temp_dir) + path = yaml_dir / "m1_m2_ratio.yaml" + data = yaml.safe_load(path.read_text()) + data["empirical_data"]["n_evaluable"] = 6 + with open(path, "w") as f: + yaml.dump(data, f) + df = load_calibration_targets(yaml_dir, cohorts=self._registry(temp_dir)) + assert df.iloc[0]["sample_size"] == 6 + + def test_missing_registry_raises(self, temp_dir): + with pytest.raises(ValueError, match="no cohort registry was passed"): + load_calibration_targets(self._placed_dir(temp_dir)) + + def test_unresolvable_cohort_raises(self, temp_dir): + yaml_dir = self._placed_dir(temp_dir, cohort_id="ghost") + with pytest.raises(ValueError, match="does not resolve"): + load_calibration_targets(yaml_dir, cohorts=self._registry(temp_dir)) + + def test_declared_sample_size_needs_no_registry(self, single_baseline_dir): + """A mechanistic target states its own n.""" + assert load_calibration_targets(single_baseline_dir).iloc[0]["sample_size"] == 30 + + class TestLoadCalibrationTargets: def test_load_single_baseline_yaml(self, single_baseline_dir): """Load one YAML, check DataFrame columns and values.""" @@ -499,7 +562,7 @@ def test_wrapper_no_constants_still_passes_empty_dict(self): """Wrapper with no constants still passes _constants={} to observable.""" code = _generate_wrapper_code( observable_code=( - "def compute_observable(time, species_dict, constants):\n" " return 42\n" + "def compute_observable(time, species_dict, constants):\n return 42\n" ), constants=[], readout_time=0.0, @@ -513,7 +576,7 @@ def test_wrapper_handles_scalar_result(self): """Wrapper that returns a plain scalar works correctly.""" code = _generate_wrapper_code( observable_code=( - "def compute_observable(time, species_dict, constants):\n" " return 42.0\n" + "def compute_observable(time, species_dict, constants):\n return 42.0\n" ), constants=[], readout_time=0.0,