diff --git a/PUBLIC_SNAPSHOT_MANIFEST.sha256 b/PUBLIC_SNAPSHOT_MANIFEST.sha256 index 7ac5ba8..999ba34 100644 --- a/PUBLIC_SNAPSHOT_MANIFEST.sha256 +++ b/PUBLIC_SNAPSHOT_MANIFEST.sha256 @@ -10,12 +10,12 @@ b245147442a93a898ee7deb552f89c3744ef79c9864147ce16739bdbe3d73705 CONTRIBUTING.m c8c9b410f264c2b10136c11cadf98cc88f0226388b6150df7932231d8b20820b README.md 08ca2198b1ccadb2fa14a57241a0926786ccddb0494b7b2de17671f5c465ea1a SECURITY.md ea53afdda1f95592d41033b587609b227ddbbd3255537064173ab76f0406c748 docs/calibration.md -4a7701edadf5680a2a6e355593ea93b73e4294ddc150454b3a60216fc9cd96ab docs/data_contracts.md -fa8c6654114a41d73cbcf9f766443f513a578b9650205d58e750f94a5512b26d docs/file_format.md +24b834ffc13d1adddb6a2b2696542ac4a28e1b1abecf556156ce26d16ff23660 docs/data_contracts.md +2dac1b577e09b6335d720edd6c81364f14310246dbb341b0e0f6fa64bf85489b docs/file_format.md 732ba1632cacf42e1ca2429d508adbbfd0c22faeea49ad7ac61250387aa91cbb docs/manuscript_reproduction.md 1fb721d4e828a6ab1ab40857cac70379cafbed9c5fface23874c65f64380b07d docs/public_release.md 2dd72b343c93b7ea1605eea17429ffb7e0949abd42c399b6bc26bfc507651467 docs/replicated_studies.md -92dfe0eb82344212371be94ea7537ba61afb527a1a207fddea179da6ad5e37f3 docs/scientific_contracts.md +40cfe4a0c3244f8784afad01f13bb0cde4d5ba82fba9eb71d14e1d7ad4c6f2a1 docs/scientific_contracts.md f5e13d205b9aec0382592897dba524ac5f56be5068d45bc881c0a225cb4a2f7e docs/support.md ca5265d3a1d4e383c75223321dd816761dfbdb6c45270a99e7a8a936d3caf496 docs/templates/paqxos_raw_intensity_report_template.txt 927740f6071333ec7c9717a83a1b73e9fd250d91dc13bccff42bcb022dc94570 docs/testing.md @@ -33,11 +33,11 @@ dddc8ff632e28f0149cd173bc5a3f6183dd5152af7880b96799d9842cb7de5f5 scripts/build_ afa1dc545653f63ae2c263e20381ce0f58815677289a0f0010af10faaec4e518 src/diffractomorph_pipeline/assay/timecourse.py e24eeb18138d29575020003b11e0f52de63b385c3974890de4f0f7825af93612 src/diffractomorph_pipeline/band_routing.py c041ec714aad40a3498fe4b4f7750fcb2daf476ec1adf6e13da47677ee0a4041 src/diffractomorph_pipeline/batch.py -96e749974bab245c4b4d7513104be90e44f2fa26a15d0226052beb8a964313c6 src/diffractomorph_pipeline/cli.py +1e78213f7cf2dbfc59252378e61e246aa288b924f51129d50a90a4f37cfd5803 src/diffractomorph_pipeline/cli.py 15938e8d4ad0a82b1059955b345d27e2594ea4f064839659a72a8911cad22cfb src/diffractomorph_pipeline/config.py 6fce3b759452e08345aa6ce96b3e81e8e572e7622a1afa405bbb3a620e49c392 src/diffractomorph_pipeline/data/examples/synthetic_minimal/README.md 956d3543d4c20da9dfbd0c7f4bdaa62d3cc575c33ecdf2238c2d9f03873f78c9 src/diffractomorph_pipeline/data/examples/synthetic_minimal/compound_x_run.csv -02562e83cd24e3d778e3df00d7a17b1212ceb2c6da737759b79cb377d5084e30 src/diffractomorph_pipeline/data/examples/synthetic_minimal/project.yaml +9f0726702033512bd90742a7cf1efb7c47452ad029bcdf2b664a1f0e0d35d290 src/diffractomorph_pipeline/data/examples/synthetic_minimal/project.yaml 0ed9c39d3b27be5dfc0c93c2095aac9a40c6687e16646c85c166137eed7ce0b6 src/diffractomorph_pipeline/data/schema/project_manifest_v1.json 3b0ebb4318ab55aef723f5c506faa68f90a8c695ea484435d80309340d9f526e src/diffractomorph_pipeline/empirical_fit.py 5f6dca7167e2305655a7e1b557f84fdfb35b1a52b2b41e37a522a05629e9e017 src/diffractomorph_pipeline/extract.py @@ -54,12 +54,12 @@ bbdb5a876e5c6d926bc4e78ba5c503f226eec4467f63eb2162e976153be1cbd6 src/diffractom cb7ba9fdba100f8fa58bacbd116e40d268a9b18f1be0fecf0c6848699a022871 src/diffractomorph_pipeline/forward/registry.py f2da49d5de89b5a1f5677c3f8335300bb2d6f9164b62db56172f5cbd6a31e785 src/diffractomorph_pipeline/forward/surface.py d838843c3a5da9ad9f6e52018a54713c9965dce3ca2fb28721ccb6619f137519 src/diffractomorph_pipeline/forward/surface_ode.py -8b9ea72ad8825b49e8e392fdc7a5e3a7dfe93f43e361a68d85ef38f573a1ff13 src/diffractomorph_pipeline/ingest.py +4e719b32015f38a13c9e1061724eacf25218995c24d8348decc979a5596356f4 src/diffractomorph_pipeline/ingest.py b031fd13e224b5c71b5a093271816392876117184c84135e59e4ad8a81845441 src/diffractomorph_pipeline/io/__init__.py 9dab71331cbac58ed797c910857a5d65995e1504851fce7945cefada76fb9e90 src/diffractomorph_pipeline/io/base.py -714fbafd4ca92ab038ce953163cd1f9336a3bb002da167e891c124a1089b2707 src/diffractomorph_pipeline/io/paqxos_rtf.py +96e24067d399d922ab5d5c4db92c6b280dd68911d63eb59b3d00ea79ba3bbb2d src/diffractomorph_pipeline/io/paqxos_rtf.py 56c9b981c6f2230686b4994b214e61262458d3a76708e46165e734cda996acbf src/diffractomorph_pipeline/io/tidy_csv.py -b84f07cbf199059b43d83c6c7785ad1e1d59189a0f86c9f0f799b2c7e3825f57 src/diffractomorph_pipeline/kinetics.py +3d0e39e3cb3e26046fec44e1e80bc0880ddc0d6238983e96e7fb91468b936f60 src/diffractomorph_pipeline/kinetics.py 7199d8c79c935fc9d2ed0b99dd442f8028ea40ca5858a2ced952bde14f39c08b src/diffractomorph_pipeline/model.py 5d0b4fb28f6e7bdf05f438f59588dd366d8f50f5e1dba7d7fcb832f707a58d5d src/diffractomorph_pipeline/noise.py 2cad4e84ba25f8b59634f56efa391c923d46a5d4554a99240bf872de1c47c90f src/diffractomorph_pipeline/noise_filter.py @@ -70,15 +70,16 @@ fe9cc0ddd41f988e26bc86a2a3091c2a06e33737349abfc51da096093b35a3c8 src/diffractom 777c2ddc9a7f20def45982d424771fd8cd36e6a1ed7e80b90918ad1d7250d495 src/diffractomorph_pipeline/optics/mie_candidate.py 68a3282d814a5458c37578e8811447846dc2bceb259804be729c1dc3191f326c src/diffractomorph_pipeline/optics/standards.py a05b30fc6e1526a7db5996d3eec8245ff487713c68f00625d6cd1b5d52563c94 src/diffractomorph_pipeline/plot_styles.py -6449b66d738d1891bf00a42351819ac5865a0d7a31a8ae2bd3b5669d03f3478d src/diffractomorph_pipeline/processing/__init__.py -59106077fee9aa0e91a5add36832783ea163c6d973c8baa2db3acb7b9fa37103 src/diffractomorph_pipeline/processing/aggregate.py +6ae6a6cca7cf5c2f1296308bbead7fa22ee0958576bf792d9c32b4a181bce9eb src/diffractomorph_pipeline/processing/__init__.py +569f3a480553039e66df1eacafce2af432fc55baa3d325131bb0d727bf4dcd1f src/diffractomorph_pipeline/processing/aggregate.py 7396b330612464191ede4d61ea40b58eda10b05d8fae15534a3ae2a13a6538d8 src/diffractomorph_pipeline/processing/artifacts.py 6b06ffa13586a43efa952b4277cb85ed61757668adc41358bbd4bf6ddb0812f8 src/diffractomorph_pipeline/processing/matched_extent.py 85fc03e83ce13823e2dd5c1e0324757f204e2051d832d38512d56b375f7c5625 src/diffractomorph_pipeline/psd.py f81037d4b7763cfb569b819a81c978f9fa89e1c51d5a7e8574f562cd81703703 src/diffractomorph_pipeline/solubility.py 554a4d94b1c09b69a376e78df45bd7eb9e46ea56d6f5824bf7f6bfb06d2a6e57 src/diffractomorph_pipeline/study/__init__.py ff83e5387a722dee57ba72fde230cf058ceeb2fbb4a70d2a67cb99a945f17d9b src/diffractomorph_pipeline/study/aggregate.py -e96e649f877835f346fb5428dfa72049fdca9010a233201e32381c3270b8edd2 src/diffractomorph_pipeline/study/manifest.py +438fa01f3d5716f6ea75d513c6cdb882f0ddb2230fd3e308feb6746574da1991 src/diffractomorph_pipeline/study/manifest.py 2abda7cb13f378a0eebef4e38da3a8df73851bd6591ba8328902c9fc48b1ad27 src/diffractomorph_pipeline/utils.py 0bd1e1dd7a139389de0a29c3ab98a62591acc9ea26e0d2cfe766b9c2a0ce79b4 src/diffractomorph_pipeline/uv_qc.py -cf1b45c675f2518bb5fc76e4368884ed6094d44fcadc7f764f14a0b6db8dff1d tests/public/test_public_release.py +0f7c85f170c0ad97e37f5295db61c8c36437433ff1f6a39ac0b2cd85ddde7bcc tests/public/test_ingest_validation.py +0cccb4deb414ce71e36ab075af3799b37b0afb165c011e2914027fbbd8d4960d tests/public/test_public_release.py diff --git a/docs/scientific_contracts.md b/docs/scientific_contracts.md index 4f82803..68a9080 100644 --- a/docs/scientific_contracts.md +++ b/docs/scientific_contracts.md @@ -7,7 +7,14 @@ quantities, and forward predictions. These categories are not interchangeable. `processing.fit_aggregate_kww` is the manuscript-authoritative particle-side endpoint. It sums the explicitly selected measured detector channels, applies -the one-sided upward Hampel repair, and fits a free-amplitude KWW curve. The +an explicitly declared acquisition-start policy, applies the one-sided upward +Hampel repair, and fits a free-amplitude KWW curve. The released behavior uses +`start_boundary.policy: first_frame`. The optional +`concordant_early_maximum` policy requires a declared acquisition variable, +early search window, maximum startup interval, minimum concordant increase, and +angular-pattern cosine threshold; its selected original frame, elapsed time, and reason +are written with each run-level result, and elapsed time is re-zeroed at the +selected frame. The default profile selects every measured channel and retains the stored reference in the signal. A reference-adjusted fit is available only through the explicit `reference_mode="reference_adjusted"` sensitivity setting. @@ -22,12 +29,30 @@ weight independent units equally. Overall and per-endpoint contributing run and independent-unit counts are stored separately so missing values cannot inflate the reported replication for an endpoint. +An explicit optional start profile has this form: + +```yaml +start_boundary: + policy: concordant_early_maximum + acquisition_variable: copt + search_frames: 3 + maximum_time_min: 1.0 + minimum_relative_increase: 0.20 + minimum_spectral_cosine: 0.995 +``` + +These values define an analysis profile; they are not universal instrument +defaults. The concordant start policy and `processing.correct_artifacts` are +alternative startup treatments. The aggregate CLI therefore requires +`--artifact-correction off` when a manifest declares both. + ## Artifact processing `processing.correct_artifacts` requires a named Copt-like acquisition variable. It records startup removal, synchronized interpolation, isolated-channel median replacement, and gap re-zeroing in a frame-level ledger. The separate aggregate -Hampel repair remains separate in both code and provenance. +Hampel repair remains separate from acquisition-start selection in both code +and provenance. ## Matched-extent q3 diff --git a/src/diffractomorph_pipeline/cli.py b/src/diffractomorph_pipeline/cli.py index c145efb..cccfbba 100644 --- a/src/diffractomorph_pipeline/cli.py +++ b/src/diffractomorph_pipeline/cli.py @@ -168,18 +168,30 @@ def aggregate_kww_main(argv=None): project = load_manifest(args.manifest) analysis_profile = project.require_profile("analysis") required_analysis = { - "channel_set", "stored_reference_subtraction", "upward_hampel", + "channel_set", "stored_reference_subtraction", "start_boundary", "upward_hampel", "tau_bounds_min", "beta_bounds", } missing_analysis = sorted(required_analysis - set(analysis_profile.parameters)) if missing_analysis: parser.error("analysis profile missing: " + ", ".join(missing_analysis)) - config = AggregateKWWConfig.from_profile(analysis_profile.parameters) + try: + config = AggregateKWWConfig.from_profile(analysis_profile.parameters) + except ValueError as exc: + parser.error(str(exc)) artifact_parameters = analysis_profile.parameters.get("artifact_correction") artifact_config = ( ArtifactCorrectionConfig.from_profile(artifact_parameters) if artifact_parameters is not None else None ) + if ( + config.start_policy != "first_frame" + and artifact_config is not None + and args.artifact_correction != "off" + ): + parser.error( + "concordant acquisition-start selection and artifact_correction are " + "alternative startup treatments; rerun with --artifact-correction off" + ) args.output_dir.mkdir(parents=True, exist_ok=True) rows = [] for spec in project.runs: @@ -203,7 +215,10 @@ def aggregate_kww_main(argv=None): run = corrected.run corrected.ledger.to_csv(args.output_dir / f"{spec.run_id}_frame_ledger.csv", index=False) correction_status = "applied" - result = fit_aggregate_kww(run, config) + try: + result = fit_aggregate_kww(run, config) + except KeyError as exc: + parser.error(f"run {spec.run_id!r}: {exc}") condition = spec.metadata.get("condition") if condition is None: parser.error(f"measurement run {spec.run_id!r} metadata requires condition") @@ -215,11 +230,38 @@ def aggregate_kww_main(argv=None): by_run = pd.DataFrame(rows) values = ("tau_min", "beta", "mean_relax_min", "t50_min", "optical_decay_depth_pct", "i0_fit") summary = summarize_hierarchy(by_run, value_columns=values) + unit_start_counts = ( + by_run.assign(_started_late=by_run["start_index"].gt(0).astype(int)) + .groupby(["condition", "independent_unit_id"], as_index=False) + .agg(n_runs_started_late=("_started_late", "sum")) + ) + independent_units = summary.independent_units.merge( + unit_start_counts, + on=["condition", "independent_unit_id"], + how="left", + validate="one_to_one", + ) + condition_start_counts = ( + unit_start_counts.assign( + _unit_started_late=unit_start_counts["n_runs_started_late"].gt(0).astype(int), + ) + .groupby("condition", as_index=False) + .agg( + n_runs_started_late=("n_runs_started_late", "sum"), + n_independent_units_started_late=("_unit_started_late", "sum"), + ) + ) + conditions = summary.conditions.merge( + condition_start_counts, + on="condition", + how="left", + validate="one_to_one", + ) by_run.to_csv(args.output_dir / "aggregate_kww_by_run.csv", index=False) - summary.independent_units.to_csv( + independent_units.to_csv( args.output_dir / "aggregate_kww_by_independent_unit.csv", index=False, ) - summary.conditions.to_csv(args.output_dir / "aggregate_kww_by_condition.csv", index=False) + conditions.to_csv(args.output_dir / "aggregate_kww_by_condition.csv", index=False) print( f"runs={len(by_run)} independent_units={summary.independent_units['independent_unit_id'].nunique()} " f"conditions={len(summary.conditions)} reference_mode={config.reference_mode}" diff --git a/src/diffractomorph_pipeline/data/examples/synthetic_minimal/project.yaml b/src/diffractomorph_pipeline/data/examples/synthetic_minimal/project.yaml index df6410d..fcf3e29 100644 --- a/src/diffractomorph_pipeline/data/examples/synthetic_minimal/project.yaml +++ b/src/diffractomorph_pipeline/data/examples/synthetic_minimal/project.yaml @@ -23,6 +23,8 @@ profiles: parameters: channel_set: all_measured stored_reference_subtraction: false + start_boundary: + policy: first_frame upward_hampel: half_window_frames: 3 threshold_mad: 4.0 diff --git a/src/diffractomorph_pipeline/kinetics.py b/src/diffractomorph_pipeline/kinetics.py index c8578a5..360954b 100644 --- a/src/diffractomorph_pipeline/kinetics.py +++ b/src/diffractomorph_pipeline/kinetics.py @@ -1,8 +1,8 @@ """Scattering-intensity dissolution kinetics — total angular signal ΣI(t) → KWW fit + shape-robust rate. The **overall scattering intensity** workflow (the LD dissolution readout the project settled on): sum -the ring-channel intensities into one signal ``ΣI(t)`` (total angular scattering, ∝ particle area), clean -upward optical glitches (bubbles), and fit the free-amplitude KWW / stretched exponential that a +the ring-channel intensities into one signal ``ΣI(t)`` (total angular scattering, ∝ particle area), screen +brief upward acquisition excursions, and fit the free-amplitude KWW / stretched exponential that a polydisperse population's superposed single-particle decays produce (see :mod:`empirical_fit`). The **meaningful readouts** (per the optics analysis — the ΣI β is a scattering-weighted, uniformity- @@ -11,7 +11,7 @@ - ``mean_relax_min`` = ``⟨t⟩ = (τ/β)·Γ(1/β)`` — the β-decoupled "how fast" timescale (or ``t50``), - ``beta`` — decay **heterogeneity** (β→1 uniform/single-exponential; β<1 a spread of dissolution times), - ``depth`` — extent of the decay, -- ``i0`` — back-extrapolated t=0 signal (recovers the start before the ~30 s injection→first-frame delay). +- ``i0`` — fitted signal at the analysis time origin. Compare rates across conditions with ``⟨t⟩`` / ``t50``, **not** raw ``τ``. @@ -39,9 +39,8 @@ def despike_upward(t_min, y, *, z=4.0, half=3, return_mask=False): """Remove UPWARD glitch spikes from a monotone-ish decaying LD signal (``Copt`` or the total angular ``ΣI``), on the full absolute time grid. - Bubbles / obscuration hits raise the signal above the local trend, but dissolution only *lowers* - the undissolved material — so an upward excursion is never real. A Hampel filter (rolling-median - baseline over ``±half`` frames, robust MAD, **upward tail only**), flagged frames interpolated from + A Hampel filter identifies brief upward excursions from a locally decreasing trajectory + (rolling-median baseline over ``±half`` frames, robust MAD, **upward tail only**). Flagged frames are interpolated from their good neighbours so time stays aligned with UV. A flagged point must also be **locally rising** (a contiguous above-trend run must *contain* a rising frame): the steep monotone *start* of a fast run sits above the local median without being a spike, and the rising-edge guard keeps that leading diff --git a/src/diffractomorph_pipeline/processing/__init__.py b/src/diffractomorph_pipeline/processing/__init__.py index fb4fba6..647f53c 100644 --- a/src/diffractomorph_pipeline/processing/__init__.py +++ b/src/diffractomorph_pipeline/processing/__init__.py @@ -5,12 +5,14 @@ AggregateKWWResult, aggregate_signal, fit_aggregate_kww, + select_aggregate_start, ) from .artifacts import ArtifactCorrectionConfig, ArtifactCorrectionResult, correct_artifacts from .matched_extent import MatchedExtentConfig, matched_q3_extent __all__ = [ "AggregateKWWConfig", "AggregateKWWResult", "aggregate_signal", "fit_aggregate_kww", + "select_aggregate_start", "ArtifactCorrectionConfig", "ArtifactCorrectionResult", "correct_artifacts", "MatchedExtentConfig", "matched_q3_extent", ] diff --git a/src/diffractomorph_pipeline/processing/aggregate.py b/src/diffractomorph_pipeline/processing/aggregate.py index a6a1d12..fe1b77d 100644 --- a/src/diffractomorph_pipeline/processing/aggregate.py +++ b/src/diffractomorph_pipeline/processing/aggregate.py @@ -11,6 +11,7 @@ ReferenceMode = Literal["raw_measured", "reference_adjusted"] +StartPolicy = Literal["first_frame", "concordant_early_maximum"] @dataclass(frozen=True) @@ -19,13 +20,45 @@ class AggregateKWWConfig: reference_mode: ReferenceMode = "raw_measured" upward_hampel_z: float = 4.0 upward_hampel_half_window: int = 3 + start_policy: StartPolicy = "first_frame" + start_acquisition_variable: str = "copt" + start_search_frames: int = 3 + start_maximum_time_min: float | None = None + start_minimum_relative_increase: float | None = None + start_minimum_spectral_cosine: float | None = None tau_bounds_min: tuple[float, float] = (0.05, 500.0) beta_bounds: tuple[float, float] = (0.2, 3.0) + def validate_start_boundary(self) -> None: + """Validate the optional start policy before any run is processed.""" + if self.start_policy == "first_frame": + return + if self.start_policy != "concordant_early_maximum": + raise ValueError(f"unsupported start policy: {self.start_policy!r}") + if self.start_search_frames < 2: + raise ValueError("start_search_frames must be at least 2") + if self.start_maximum_time_min is None or self.start_maximum_time_min < 0: + raise ValueError("start_maximum_time_min must be declared and nonnegative") + if ( + self.start_minimum_relative_increase is None + or self.start_minimum_relative_increase < 0 + ): + raise ValueError("start_minimum_relative_increase must be declared and nonnegative") + if ( + self.start_minimum_spectral_cosine is None + or not 0 <= self.start_minimum_spectral_cosine <= 1 + ): + raise ValueError( + "start_minimum_spectral_cosine must be declared between 0 and 1" + ) + @classmethod def from_profile(cls, parameters) -> "AggregateKWWConfig": """Build from an explicit manifest analysis profile.""" - required = {"channel_set", "stored_reference_subtraction", "upward_hampel", "tau_bounds_min", "beta_bounds"} + required = { + "channel_set", "stored_reference_subtraction", "start_boundary", + "upward_hampel", "tau_bounds_min", "beta_bounds", + } missing = sorted(required - set(parameters)) if missing: raise ValueError("analysis profile missing: " + ", ".join(missing)) @@ -45,14 +78,46 @@ def from_profile(cls, parameters) -> "AggregateKWWConfig": if parameters["stored_reference_subtraction"] else "raw_measured" ) - return cls( + start = parameters["start_boundary"] + policy = str(start.get("policy", "first_frame")) + if policy not in ("first_frame", "concordant_early_maximum"): + raise ValueError( + "start_boundary.policy must be 'first_frame' or " + "'concordant_early_maximum'" + ) + if policy == "concordant_early_maximum": + required_start = { + "acquisition_variable", "search_frames", "maximum_time_min", + "minimum_relative_increase", "minimum_spectral_cosine", + } + missing_start = sorted(required_start - set(start)) + if missing_start: + raise ValueError("start_boundary missing: " + ", ".join(missing_start)) + config = cls( channel_ids=channel_ids, reference_mode=reference_mode, upward_hampel_z=float(hampel["threshold_mad"]), upward_hampel_half_window=int(hampel["half_window_frames"]), + start_policy=policy, + start_acquisition_variable=str(start.get("acquisition_variable", "copt")), + start_search_frames=int(start.get("search_frames", 3)), + start_maximum_time_min=( + float(start["maximum_time_min"]) + if "maximum_time_min" in start else None + ), + start_minimum_relative_increase=( + float(start["minimum_relative_increase"]) + if "minimum_relative_increase" in start else None + ), + start_minimum_spectral_cosine=( + float(start["minimum_spectral_cosine"]) + if "minimum_spectral_cosine" in start else None + ), tau_bounds_min=tuple(float(value) for value in parameters["tau_bounds_min"]), beta_bounds=tuple(float(value) for value in parameters["beta_bounds"]), ) + config.validate_start_boundary() + return config @dataclass(frozen=True) @@ -64,6 +129,9 @@ class AggregateKWWResult: upward_repair_mask: np.ndarray channel_ids: tuple[str, ...] reference_mode: ReferenceMode + start_index: int + selected_elapsed_time_min: float + start_reason: str fit: dict config: dict observable: str = "angle-integrated detector signal" @@ -74,6 +142,10 @@ def to_row(self) -> dict: "run_id": self.run_id, "independent_unit_id": self.independent_unit_id, "reference_mode": self.reference_mode, + "start_policy": self.config["start_policy"], + "start_index": self.start_index, + "selected_elapsed_time_min": self.selected_elapsed_time_min, + "start_reason": self.start_reason, "n_channels": len(self.channel_ids), "n_upward_repairs": int(self.upward_repair_mask.sum()), **self.fit, @@ -101,13 +173,81 @@ def aggregate_signal(run: Run, config: AggregateKWWConfig | None = None) -> tupl return signal, tuple(selected) +def select_aggregate_start( + run: Run, + aggregate: np.ndarray, + channel_ids: tuple[str, ...], + config: AggregateKWWConfig, +) -> tuple[int, str]: + """Select the first analysis frame under an explicit acquisition-start policy. + + ``first_frame`` preserves the released behavior. ``concordant_early_maximum`` + recognizes incomplete initial circulation only when the same later frame is + the early maximum of both total angular signal and a declared acquisition + variable within the declared startup interval, both increases meet or exceed + the declared minimum, and the angular pattern remains shape-preserving. If + those conditions are not met, frame 0 remains the start. The caller re-zeros + time after selection. + """ + if config.start_policy == "first_frame": + return 0, "first_frame" + config.validate_start_boundary() + + aggregate = np.asarray(aggregate, dtype=float) + if aggregate.shape != run.time_min.shape: + raise ValueError("aggregate must contain one value per run frame") + acquisition = run.acquisition_variable(config.start_acquisition_variable) + n_early = min(config.start_search_frames, aggregate.size) + if not ( + np.all(np.isfinite(aggregate[:n_early])) + and np.all(np.isfinite(acquisition[:n_early])) + ): + return 0, "nonfinite_early_start_variable" + signal_index = int(np.argmax(aggregate[:n_early])) + acquisition_index = int(np.argmax(acquisition[:n_early])) + if signal_index == 0 or signal_index != acquisition_index: + return 0, "no_concordant_early_maximum" + selected_time = float(run.time_min[signal_index] - run.time_min[0]) + if selected_time > config.start_maximum_time_min: + return 0, "early_maximum_after_startup_interval" + + baseline_signal = float(aggregate[0]) + baseline_acquisition = float(acquisition[0]) + if baseline_signal <= 0 or baseline_acquisition <= 0: + return 0, "nonpositive_start_variable" + signal_increase = float(aggregate[signal_index] / baseline_signal - 1.0) + acquisition_increase = float(acquisition[signal_index] / baseline_acquisition - 1.0) + minimum = config.start_minimum_relative_increase + if signal_increase < minimum or acquisition_increase < minimum: + return 0, "early_maximum_below_minimum_increase" + + indices = [run.channel_ids.index(channel) for channel in channel_ids] + patterns = np.asarray(run.signal[:, indices], dtype=float) + if config.reference_mode == "reference_adjusted": + if run.stored_reference is None: + raise ValueError("reference_adjusted start selection requires a stored reference") + patterns = patterns - run.stored_reference[indices][None, :] + first_pattern = patterns[0] + peak_pattern = patterns[signal_index] + denominator = float(np.linalg.norm(first_pattern) * np.linalg.norm(peak_pattern)) + cosine = float(np.dot(first_pattern, peak_pattern) / denominator) if denominator else float("nan") + if not np.isfinite(cosine) or cosine < config.start_minimum_spectral_cosine: + return 0, "early_maximum_changed_angular_pattern" + return signal_index, "concordant_early_maximum" + + def fit_aggregate_kww(run: Run, config: AggregateKWWConfig | None = None) -> AggregateKWWResult: """Fit the all-measured-channel aggregate with a free-amplitude KWW descriptor.""" config = config or AggregateKWWConfig() signal, selected = aggregate_signal(run, config) + start_index, start_reason = select_aggregate_start(run, signal, selected, config) + original_start_time = float(run.time_min[start_index]) + selected_elapsed_time = original_start_time - float(run.time_min[0]) + analysis_time = run.time_min[start_index:] - original_start_time + analysis_signal = signal[start_index:] time_min, cleaned, repaired = kinetics.despike_upward( - run.time_min, - signal, + analysis_time, + analysis_signal, z=config.upward_hampel_z, half=config.upward_hampel_half_window, return_mask=True, @@ -126,6 +266,9 @@ def fit_aggregate_kww(run: Run, config: AggregateKWWConfig | None = None) -> Agg upward_repair_mask=repaired, channel_ids=selected, reference_mode=config.reference_mode, + start_index=start_index, + selected_elapsed_time_min=selected_elapsed_time, + start_reason=start_reason, fit=fit, config=asdict(config), ) diff --git a/tests/public/test_public_release.py b/tests/public/test_public_release.py index ec0487f..ab1fc7a 100644 --- a/tests/public/test_public_release.py +++ b/tests/public/test_public_release.py @@ -8,11 +8,16 @@ from pathlib import Path import numpy as np +import pandas as pd import pytest +import yaml import diffractomorph_pipeline as dfm from diffractomorph_pipeline.assay import AssayCalibration -from diffractomorph_pipeline.processing import AggregateKWWConfig +from diffractomorph_pipeline.model import Run, RunProvenance +from diffractomorph_pipeline.processing import ( + AggregateKWWConfig, fit_aggregate_kww, select_aggregate_start, +) from diffractomorph_pipeline.study import bundled_example_manifest, load_manifest @@ -54,6 +59,167 @@ def test_public_analysis_profile_is_explicit(): config = AggregateKWWConfig.from_profile(project.require_profile("analysis").parameters) assert config.channel_ids is None assert config.reference_mode == "raw_measured" + assert config.start_policy == "first_frame" + + +def _synthetic_run(signal, copt, time_min=None): + signal = np.asarray(signal, dtype=float) + if time_min is None: + time_min = np.arange(signal.shape[0], dtype=float) + return Run( + signal=signal, + channel_ids=tuple(f"ch{index + 1}" for index in range(signal.shape[1])), + time_min=np.asarray(time_min, dtype=float), + acquisition={"copt": np.asarray(copt, dtype=float)}, + provenance=RunProvenance( + run_id="synthetic-start-test", source_path="synthetic.csv", + adapter="tidy_csv", sample_id="compound-x", + ), + run_kind="measurement", + ) + + +def test_concordant_early_maximum_selects_and_rezeros_shape_preserving_rise(): + pattern = np.array([4.0, 3.0, 2.0, 1.0]) + scales = np.array([1.0, 1.5, 1.35, 1.0, 0.7, 0.5, 0.4]) + run = _synthetic_run(scales[:, None] * pattern, [1.0, 1.45, 1.3, 1.0, 0.8, 0.6, 0.5]) + config = AggregateKWWConfig( + start_policy="concordant_early_maximum", + start_search_frames=3, + start_maximum_time_min=2.0, + start_minimum_relative_increase=0.2, + start_minimum_spectral_cosine=0.995, + ) + result = fit_aggregate_kww(run, config) + assert result.start_index == 1 + assert result.start_reason == "concordant_early_maximum" + assert result.time_min[0] == 0.0 + assert result.selected_elapsed_time_min == 1.0 + + +def test_concordant_early_maximum_rejects_shape_changing_rise(): + run = _synthetic_run( + [[4, 3, 2, 1], [1, 2, 3, 8], [3, 2, 1.5, 0.8], [2, 1.5, 1, 0.5], + [1.5, 1, 0.8, 0.4], [1.0, 0.8, 0.6, 0.3], [0.8, 0.6, 0.4, 0.2]], + [1.0, 1.5, 1.3, 1.0, 0.8, 0.6, 0.5], + ) + config = AggregateKWWConfig( + start_policy="concordant_early_maximum", + start_search_frames=3, + start_maximum_time_min=2.0, + start_minimum_relative_increase=0.2, + start_minimum_spectral_cosine=0.995, + ) + aggregate = run.signal.sum(axis=1) + index, reason = select_aggregate_start(run, aggregate, run.channel_ids, config) + assert index == 0 + assert reason == "early_maximum_changed_angular_pattern" + + +def test_first_frame_policy_preserves_released_start_behavior(): + pattern = np.array([4.0, 3.0, 2.0, 1.0]) + scales = np.array([1.0, 1.5, 1.35, 1.0, 0.7, 0.5, 0.4]) + run = _synthetic_run(scales[:, None] * pattern, scales) + result = fit_aggregate_kww(run, AggregateKWWConfig(start_policy="first_frame")) + assert result.start_index == 0 + assert result.selected_elapsed_time_min == 0.0 + assert result.start_reason == "first_frame" + + +def test_concordant_start_requires_declared_acquisition_variable(): + run = _synthetic_run([[4, 3], [6, 4.5], [5, 4]], [1.0, 1.5, 1.3]) + config = AggregateKWWConfig( + start_policy="concordant_early_maximum", + start_acquisition_variable="transmission", + start_maximum_time_min=2.0, + start_minimum_relative_increase=0.2, + start_minimum_spectral_cosine=0.995, + ) + with pytest.raises(KeyError, match="transmission"): + fit_aggregate_kww(run, config) + + +def test_concordant_start_rejects_nonfinite_early_acquisition(): + run = _synthetic_run([[4, 3], [6, 4.5], [5, 4]], [1.0, np.nan, 1.3]) + config = AggregateKWWConfig( + start_policy="concordant_early_maximum", + start_maximum_time_min=2.0, + start_minimum_relative_increase=0.2, + start_minimum_spectral_cosine=0.995, + ) + aggregate = run.signal.sum(axis=1) + index, reason = select_aggregate_start(run, aggregate, run.channel_ids, config) + assert index == 0 + assert reason == "nonfinite_early_start_variable" + + +def test_concordant_start_profile_requires_complete_explicit_contract(): + project = load_manifest(bundled_example_manifest()) + parameters = dict(project.require_profile("analysis").parameters) + parameters["start_boundary"] = { + "policy": "concordant_early_maximum", + "acquisition_variable": "copt", + "search_frames": 3, + "maximum_time_min": 1.0, + "minimum_relative_increase": 0.2, + "minimum_spectral_cosine": 0.995, + } + config = AggregateKWWConfig.from_profile(parameters) + assert config.start_policy == "concordant_early_maximum" + assert config.start_maximum_time_min == 1.0 + assert config.start_minimum_relative_increase == 0.2 + assert config.start_minimum_spectral_cosine == 0.995 + + incomplete = dict(parameters) + incomplete["start_boundary"] = dict(parameters["start_boundary"]) + incomplete["start_boundary"].pop("maximum_time_min") + with pytest.raises(ValueError, match="maximum_time_min"): + AggregateKWWConfig.from_profile(incomplete) + + missing_boundary = dict(parameters) + missing_boundary.pop("start_boundary") + with pytest.raises(ValueError, match="start_boundary"): + AggregateKWWConfig.from_profile(missing_boundary) + + invalid = dict(parameters) + invalid["start_boundary"] = dict(parameters["start_boundary"]) + invalid["start_boundary"]["maximum_time_min"] = -1 + with pytest.raises(ValueError, match="nonnegative"): + AggregateKWWConfig.from_profile(invalid) + + +def test_concordant_start_rejects_candidate_after_startup_interval(): + pattern = np.array([4.0, 3.0, 2.0, 1.0]) + scales = np.array([1.0, 1.5, 1.3]) + run = _synthetic_run( + scales[:, None] * pattern, + [1.0, 1.5, 1.3], + time_min=[0.0, 1.2, 1.4], + ) + config = AggregateKWWConfig( + start_policy="concordant_early_maximum", + start_search_frames=3, + start_maximum_time_min=1.0, + start_minimum_relative_increase=0.2, + start_minimum_spectral_cosine=0.995, + ) + aggregate = run.signal.sum(axis=1) + index, reason = select_aggregate_start(run, aggregate, run.channel_ids, config) + assert index == 0 + assert reason == "early_maximum_after_startup_interval" + + +def test_start_selection_rejects_misaligned_aggregate(): + run = _synthetic_run([[4, 3], [6, 4.5], [5, 4]], [1.0, 1.5, 1.3]) + config = AggregateKWWConfig( + start_policy="concordant_early_maximum", + start_search_frames=3, + start_maximum_time_min=1.0, + start_minimum_relative_increase=0.2, + start_minimum_spectral_cosine=0.995, + ) + with pytest.raises(ValueError, match="one value per run frame"): + select_aggregate_start(run, np.array([7.0, 10.5]), run.channel_ids, config) def test_python_310_public_config_entrypoint_imports(): @@ -117,6 +283,83 @@ def test_commands_do_not_claim_withheld_qc_defaults(): assert missing_qc_reference.value.code == 2 +def test_aggregate_cli_writes_start_provenance_and_summary_counts(tmp_path): + from diffractomorph_pipeline import cli + + output = tmp_path / "aggregate" + cli.aggregate_kww_main([ + str(bundled_example_manifest()), "--output-dir", str(output), + "--artifact-correction", "off", + ]) + by_run = pd.read_csv(output / "aggregate_kww_by_run.csv") + by_unit = pd.read_csv(output / "aggregate_kww_by_independent_unit.csv") + by_condition = pd.read_csv(output / "aggregate_kww_by_condition.csv") + assert { + "start_policy", "start_index", "selected_elapsed_time_min", "start_reason", + }.issubset(by_run.columns) + assert by_run.loc[0, "start_index"] == 0 + assert by_unit.loc[0, "n_runs_started_late"] == 0 + assert by_condition.loc[0, "n_runs_started_late"] == 0 + assert by_condition.loc[0, "n_independent_units_started_late"] == 0 + + +def test_aggregate_cli_rejects_two_startup_treatments(tmp_path): + from diffractomorph_pipeline import cli + + source_manifest = bundled_example_manifest() + payload = yaml.safe_load(source_manifest.read_text()) + payload["data_root"] = str(source_manifest.parent) + analysis = payload["profiles"]["analysis"]["parameters"] + analysis["start_boundary"] = { + "policy": "concordant_early_maximum", + "acquisition_variable": "copt", + "search_frames": 3, + "maximum_time_min": 1.0, + "minimum_relative_increase": 0.2, + "minimum_spectral_cosine": 0.995, + } + analysis["artifact_correction"] = { + "acquisition_variable": "copt", + "synchronized_intensity_z": 4.0, + "synchronized_acquisition_z": 1.5, + "synchronized_half_window": 2, + "isolated_spike_mad": 5.0, + "gap_threshold_min": 2.0, + } + manifest = tmp_path / "two-startup-treatments.yaml" + manifest.write_text(yaml.safe_dump(payload, sort_keys=False)) + with pytest.raises(SystemExit) as stopped: + cli.aggregate_kww_main([ + str(manifest), "--output-dir", str(tmp_path / "output"), + ]) + assert stopped.value.code == 2 + + +def test_aggregate_cli_reports_missing_start_acquisition_cleanly(tmp_path): + from diffractomorph_pipeline import cli + + source_manifest = bundled_example_manifest() + payload = yaml.safe_load(source_manifest.read_text()) + payload["data_root"] = str(source_manifest.parent) + start = payload["profiles"]["analysis"]["parameters"]["start_boundary"] + start.update({ + "policy": "concordant_early_maximum", + "acquisition_variable": "transmission", + "search_frames": 3, + "maximum_time_min": 1.0, + "minimum_relative_increase": 0.2, + "minimum_spectral_cosine": 0.995, + }) + manifest = tmp_path / "missing-start-acquisition.yaml" + manifest.write_text(yaml.safe_dump(payload, sort_keys=False)) + with pytest.raises(SystemExit) as stopped: + cli.aggregate_kww_main([ + str(manifest), "--output-dir", str(tmp_path / "output"), + "--artifact-correction", "off", + ]) + assert stopped.value.code == 2 + + def test_diagnostic_bundle_excludes_paths_and_research_values(tmp_path): from diffractomorph_pipeline import cli