diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index c4dba77a5..d6b970040 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -468,6 +468,14 @@ struct EndMessage { // not the same as every image being clean. Offline (rugnux) only; the broker does not merge. std::vector sweep_quality; std::vector sweep_quality_reasons; + + // Per-image disposition: what became of the image's observations in the merged data, as an index + // into frame_disposition_codes (FrameDisposition in image_analysis/scale_merge/Merge.h - merged, + // downgraded, rejected). Written to /entry/MX/frameDisposition with the vocabulary alongside it in + // /entry/MX/frameDispositionCodes. The sweep-quality code above says what was SEEN over a stretch; + // this says what was DONE about it. Offline (rugnux) only; the broker does not merge. + std::vector frame_disposition; + std::vector frame_disposition_codes; }; struct MetadataMessage { diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 659d707a7..a79c1654c 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -211,6 +211,24 @@ K. Diederichs, "Linking crystallographic model and data quality" (2012), Science P. A. Karplus, "Better models by discarding data?" (2013), Acta Cryst. D69, 1215-1222 [doi:10.1107/S0907444913001121](https://doi.org/10.1107/S0907444913001121). +**The frame disposition** reports what keeping each stretch of a rotation sweep costs the merged +intensities as delta-CC1/2, the change in the overall CC1/2 when a group of images is left out, +measured in the sigma-tau form so that no random half-dataset split is involved. The statistic, its +standard error going as the inverse square root of the reflection count, and the reading of it (a +group whose delta-CC1/2 is positive or near zero is not evidence of harm) are taken from its authors, +whose XDSCC12 is the reference implementation. +G. Assmann, W. Brehm and K. Diederichs, "Identification of rogue datasets in serial crystallography" +(2016), J. Appl. Cryst. 49, 1021-1028 +[doi:10.1107/S1600576716005471](https://doi.org/10.1107/S1600576716005471); G. M. Assmann, M. Wang and +K. Diederichs, "Making a difference in multi-data-set crystallography: simple and deterministic +data-scaling/selection methods" (2020), Acta Cryst. D76, 636-652 +[doi:10.1107/S2059798320006348](https://doi.org/10.1107/S2059798320006348). That the same statistic +belongs at scaling, applied to groups of images rather than to whole datasets, follows +[DIALS](https://dials.github.io/) (`dials.scale`, delta-CC1/2 image-group filtering): +J. Beilsten-Edmands, G. Winter, R. Gildea et al., "Scaling diffraction data in the DIALS software +package: algorithms and new approaches for multi-crystal scaling" (2020), Acta Cryst. D76, 385-399 +[doi:10.1107/S2059798320003198](https://doi.org/10.1107/S2059798320003198). + **Uncertainty conventions** follow the IUCr Commission on Crystallographic Nomenclature: D. Schwarzenbach, S. C. Abrahams, H. D. Flack et al., "Statistical descriptors in crystallography: Report of the IUCr Subcommittee on Statistical Descriptors" (1989), Acta Cryst. A45, 63-75 diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 3ff5cec6f..37480f7c4 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.169 +* rugnux reports what became of every frame of a rotation sweep - how many were merged, downgraded and rejected, in frames and in degrees - and what keeping each flagged stretch costs the merged intensities (delta-CC1/2); nothing is excluded on the strength of it. * The viewer draws its spot markers over a black outline, so they stay visible on a light colour map and are no longer mistaken for a magenta bad pixel or the coral beam stop. * The viewer has a font size of its own - View -> Font size, or Ctrl+plus and Ctrl+minus - at 100%, 125% or 150% of whatever size the desktop asks for, remembered across restarts. * The viewer's dataset plots label their axes in full, instead of rendering every label as "...". diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 78843cc1c..fdbd25646 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -68,6 +68,7 @@ The methods draw on, and in places reimplement, solutions from: - K. Röttger, A. Endriss, J. Ihringer, S. Doyle & W. F. Kuhs, "Lattice constants and thermal expansion of H2O and D2O ice Ih between 10 and 265 K", *Acta Cryst.* **B50** (1994), 644-648 (the ice Ih cell the ring positions below 1.522 Å are calculated from). - S. Sheriff & W. A. Hendrickson, "Description of overall anisotropy in diffraction from macromolecular crystals", *Acta Cryst.* **A43** (1987), 118-121 (the overall anisotropic B tensor and its symmetry constraints), and A. N. Popov & G. P. Bourenkov, "Choice of data-collection parameters based on statistic modelling", *Acta Cryst.* **D59** (2003), 1145-1153 (the sigma-aware estimation of the anisotropy of the observed intensity distribution, part of that paper's statistic modelling). - P. R. Evans & G. N. Murshudov, "How good are my data and what is the resolution?", *Acta Cryst.* **D69** (2013), 1204-1214 (AIMLESS: the anisotropic deltaB as the range of the principal components, and diffraction limits from a cone about each principal direction). +- G. Assmann, W. Brehm & K. Diederichs, "Identification of rogue datasets in serial crystallography", *J. Appl. Cryst.* **49** (2016), 1021-1028, and G. M. Assmann, M. Wang & K. Diederichs, *Acta Cryst.* **D76** (2020), 636-652 (XDSCC12: the sigma-tau CC1/2 and the delta-CC1/2 the frame disposition reports). - K. Diederichs & P. A. Karplus, *Nat. Struct. Biol.* **4** (1997), 269-275, and P. A. Karplus & K. Diederichs, *Science* **336** (2012), 1030-1033 (R_meas / R_pim, CC1/2 and CC\*). - IUCr Commission on Crystallographic Nomenclature, "Statistical descriptors in crystallography", *Acta Cryst.* **A45** (1989), 63-75, and *Acta Cryst.* **A51** (1995), 565-569 (uncertainty conventions). diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index 23230649c..867c9c654 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -378,6 +378,8 @@ Each surface is **cross-validated**: fitted on even-numbered frames and kept onl Each batch's $B$ is fitted on **resolution-shell means**, not on single observations: $\ln(I_\mathrm{ref}/I_\mathrm{obs})$ of one weak observation is unbounded and biased downwards — the observation appears in the response and in its own weight, and the logarithm needs $I_\mathrm{obs} > 0$, which keeps only the upward half of the noise — and on decayed data that bias grows with dose until it reverses the sign of the answer. The shells are laid inside the range the run actually diffracted to, and the fit carries an intercept as well as a slope, so a batch that is merely *dimmer* than the run (an attenuated beam, a mis-fitted frame scale) is not reported as damage. A batch whose shells are too weak to fit, or whose solved value reaches the bound the smoothing solve clamps to, is reported as **absent** rather than as a number. The first→last headline is reported only where a straight line describes the curve: radiation damage is progressive, so a curve that dips and recovers is a disturbance rather than dose, and is left to the sweep-quality report (`docs/RUGNUX_REPORT.md`) to name. +**Frame disposition (rotation, report-only).** After the correction surfaces are fitted, and on the corrected fulls, rugnux measures **Δ*CC*₁⁄₂** — the overall *CC*₁⁄₂ of the merged data with a group of images minus the *CC*₁⁄₂ without it, over the reflections that group touches — for each 10° batch of the sweep and for each stretch the sweep-quality diagnostic flagged. It is computed in the σ-τ form (no random half-dataset split, so the answer is the same every run) with each reflection's error variance taken from the **observed** scatter of its own observations rather than from the error model's σ's: a bad stretch claims the same σ's as a good one, so an error-model estimate would read a batch that adds noise as a batch that adds precision, inverting the sign of the measurement. The leave-one-out is a subtraction of the group's own $(n, \sum I, \sum I^2)$ from the per-reflection totals, so measuring a group costs one pass over its observations rather than a re-merge, and reflections the group holds the only observations of are excluded from both sides — the published statistic's own restriction. Its standard error is reported beside it as that of a single *CC*₁⁄₂ on the same reflection count, $(1-CC^2)/\sqrt{n_\mathrm{refl}-3}$, and a Δ*CC*₁⁄₂ smaller than that says nothing; a group whose Δ*CC*₁⁄₂ is positive or near zero is not evidence of harm. **Nothing is excluded on the strength of it** — the per-frame `merged` / `downgraded` / `rejected` ledger in `docs/RUGNUX_REPORT.md` records what the scaling and merging already did, and Δ*CC*₁⁄₂ measures what that cost. + ### 10.7 R-free test-set flags A fraction of the unique reflections (`rfree_fraction`, default 0.05) is flagged as a **free (test) set**, written to the output (MTZ `FreeR_flag`, mmCIF `_refln.status` = `f`, a text-HKL column) for model validation (§14) and for downstream refinement. The flag is a pure function of the reflection's **Friedel-merged (Laue) ASU index**, which gives three properties: diff --git a/docs/HDF5.md b/docs/HDF5.md index d514a6122..ca8b2a57b 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -352,6 +352,7 @@ In legacy/VDS mode these live in the data files and are linked/virtual-stacked i | `imageScaleCC` | | on-the-fly scaling correlation coefficient | | `imageScaleMosaicity` | deg | scaling-model mosaicity | | `sweepQuality` | | why this image's stretch of the sweep was flagged — see below | +| `frameDisposition` | | what became of this image's observations in the merged data — see below | **Per-image lattices:** `latticeIndexed` `[n_images, 9]` (Å) — the real-space lattice (flattened 3×3); `latticeIndexedExtra` `[n_images, max_extra_lattices, 9]` (Å) — additional orientation @@ -402,13 +403,21 @@ not, and any other value is a **1-based index into `sweepQualityReasons`**, a st beside it that carries the whole vocabulary, so the codes can be read without this source. The vocabulary is closed and stable — a code is never renamed and never reused — and currently reads `no_diffraction`, `crystal_out_of_beam`, `weak_diffraction`, `loss_of_centring`, `radiation_damage`; -[the rugnux results report](RUGNUX_REPORT.md#sweep-quality-and-the-reason-vocabulary) defines what each one -means. Both datasets are **absent** unless the sweep-quality diagnostic ran, which needs scaling and -merging; their absence therefore means "not looked for", *not* "every image clean". Written by the -offline `rugnux` path only — the broker does not merge — and not carried on the CBOR stream, in the -same way as the other offline-only fields (`space_group_number`, the refined geometry). Nothing is -excluded from processing on the strength of it. The condensed, dataset-wide form of the same finding -is in `_report.txt`. +[the rugnux results report](RUGNUX_REPORT.md#sweep-quality-the-disposition-and-their-vocabularies) +defines what each one means. Both datasets are **absent** unless the sweep-quality diagnostic ran, +which needs scaling and merging; their absence therefore means "not looked for", *not* "every image +clean". Written by the offline `rugnux` path only — the broker does not merge — and not carried on the +CBOR stream, in the same way as the other offline-only fields (`space_group_number`, the refined +geometry). Nothing is excluded from processing on the strength of it. The condensed, dataset-wide form +of the same finding is in `_report.txt`. + +**Frame disposition.** `frameDisposition` `[n_images]` (`uint8`) says what became of the image's +observations: a **0-based index into `frameDispositionCodes`**, written beside it, which reads +`merged`, `downgraded`, `rejected`. Where `sweepQuality` says what was *seen* over a stretch, this is +the ledger of what the scaling and merging already did with it — a `rejected` image contributed +nothing to the merged intensities, because it was never scaled or because an always-on guard dropped +it. It is not a decision of its own: no image is excluded on the strength of either array. Same +availability rule as `sweepQuality`: absent means the diagnostic never ran. CrystFEL can read the spots directly with: diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index d19c07fdf..9464dddee 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -173,31 +173,45 @@ inversion centre, so a reader who knows their sample is a protein reads this key run. Where no glide was found the two keys read the same, deliberately: greppability is the point, and a key that appears only sometimes has to be tested for before it can be read. -## Sweep quality and the reason vocabulary +## Sweep quality, the disposition, and their vocabularies Section 4 lists the stretches of the sweep over which the crystal delivered much less than the rest of the run — the feedback a beamline control system needs to tell an operator that a crystal should -be recentred or recollected. Nothing is excluded on the strength of it; the frames still carry -signal, and this is a message for the beamline, not a filter. +be recentred or recollected — says what became of each of them, and measures what keeping each one +costs the merged intensities. Nothing is excluded on the strength of any of it. ``` SWEEP_QUALITY_STATUS= COMPUTED -SWEEP_QUALITY_COUNT= 1 +SWEEP_QUALITY_COUNT= 2 SWEEP_QUALITY_REASONS= no_diffraction crystal_out_of_beam weak_diffraction loss_of_centring radiation_damage +SWEEP_DISPOSITIONS= merged downgraded rejected +FRAMES_MERGED= 1712 +FRAMES_DOWNGRADED= 52 +FRAMES_REJECTED= 36 +FRAMES_REJECTED_PCT= 2.00 +ROTATION_REJECTED_DEG= 7.2 SWEEP_ROTATION= 360.0 FLUX_PEAK_TO_TROUGH= 1.03 SCALE_MODULATION_PEAK_TO_TROUGH= 1.00 - FIRST_IMAGE LAST_IMAGE N_IMAGES ROTATION REASON SEVERITY SCALE CC INDEXED - ----------- ----------- --------- -------- -------------------- -------- ------ ------ -------- - 500 600 101 10.1 crystal_out_of_beam 0.83 0.12 0.30 0.02 - ----------- ----------- --------- -------- -------------------- -------- ------ ------ -------- + FIRST_IMAGE LAST_IMAGE N_IMAGES ROTATION REASON SEVERITY SCALE CC INDEXED DISPOSITION DELTA_CC_HALF DELTA_CC_HALF_SE + ----------- ----------- --------- -------- -------------------- -------- ------ ------ -------- ----------- ------------- ---------------- + 500 600 101 10.1 crystal_out_of_beam 0.83 0.12 0.30 0.02 downgraded +0.0004 0.0031 + 612 630 19 1.9 no_diffraction 0.98 0.02 0.00 0.00 rejected - - + ----------- ----------- --------- -------- -------------------- -------- ------ ------ -------- ----------- ------------- ---------------- ``` -Of those keys, only `SWEEP_QUALITY_COUNT` is in the default report; `SWEEP_QUALITY_STATUS`, -`SWEEP_QUALITY_REASONS`, `SWEEP_ROTATION`, `FLUX_PEAK_TO_TROUGH` and -`SCALE_MODULATION_PEAK_TO_TROUGH` appear with `--developer` (the default report states in prose -whether the diagnostic ran). +Of those keys, `SWEEP_QUALITY_COUNT` and the five disposition keys (`FRAMES_MERGED`, +`FRAMES_DOWNGRADED`, `FRAMES_REJECTED`, `FRAMES_REJECTED_PCT`, `ROTATION_REJECTED_DEG`) are in the +default report; `SWEEP_QUALITY_STATUS`, `SWEEP_QUALITY_REASONS`, `SWEEP_DISPOSITIONS`, +`SWEEP_ROTATION`, `FLUX_PEAK_TO_TROUGH` and `SCALE_MODULATION_PEAK_TO_TROUGH` appear with +`--developer` (the default report states in prose whether the diagnostic ran). + +The three frame counts **partition the sweep** — every processed image is exactly one of them and +they add up to the frame count — so `FRAMES_REJECTED_PCT` is the answer to "how much of this +experiment was useless". It is reported beside `ROTATION_REJECTED_DEG` on purpose: a percentage of +*frames* moves when the same experiment is re-sliced, and a percentage of the *rotation* does not. +The prose headline above the table states both. `SWEEP_QUALITY_STATUS` distinguishes **`COMPUTED`** (the diagnostic ran; a count of 0 means the sweep was clean throughout) from **`NOT_COMPUTED`** (it did not run — no scaling and merging, or stills @@ -222,10 +236,55 @@ source image is `start + ordinal * stride`); `ROTATION` — the width of the ran `SEVERITY` — the fraction of the run's typical diffracting power missing over the range, 0 (as good as the run) to 1 (nothing at all); `SCALE` and `CC` — the range's mean per-image scale and CC-to-merge relative to the run median; `INDEXED` — the fraction of the range's frames that were -scaled at all. Every range also appears as a `WARNING:` sentence in the SUMMARY. +scaled at all; `DISPOSITION` — what became of it; `DELTA_CC_HALF` and `DELTA_CC_HALF_SE` — what +keeping it costs the merged intensities, and how precisely that is known. Every range also appears as +a `WARNING:` sentence in the SUMMARY, with the same cost in words. -The same finding is written **per image** into the `_process.h5` as `/entry/MX/sweepQuality`, when -one is written — see [HDF5](HDF5.md#entry-mx-spot-finding-and-indexing-cxi-style). +### What the disposition means + +**Nothing in this section is a filter.** No observation is excluded on the strength of the sweep +quality or of Δ*CC*₁⁄₂; the disposition is a **ledger of what the scaling and merging already did**, +and Δ*CC*₁⁄₂ is a measurement of what that cost. + +| Disposition | Meaning | +|-------------|---------| +| `merged` | The frame's observations are in the merged data at their own weight. | +| `downgraded` | They are in the merged data, but over a stretch the run itself flagged, carried at the reduced weight the frame's own scale and sigmas give it. Nothing extra is subtracted: for weak-but-consistent data that reduced weight *is* the honest weight, and a second, invented per-frame weight would double-count with the σ's. | +| `rejected` | Nothing of the frame reached the merge — the frame was never scaled at all, or an earlier always-on guard dropped it because its scale had collapsed to an unusable number. To a user asking how much of the experiment was useless these are the same answer, and the `REASON` column separates them. | + +`DELTA_CC_HALF` is **Δ*CC*₁⁄₂**: the overall *CC*₁⁄₂ of the merged data **with** the range minus the +*CC*₁⁄₂ **without** it, evaluated over the reflections the range touches. Negative means keeping the +range makes the merged intensities worse. It is computed in the σ-τ form — no random half-dataset +split, so the same input gives the same answer every run — with each reflection's error variance taken +from the observed scatter of its own observations, not from the error model's σ's (a bad stretch claims +the same σ's as a good one, so an error-model estimate would read a stretch that adds noise as one that +adds precision). It is a *CC*₁⁄₂ over **the range's own reflections**, not over the whole dataset — a +range that touches a few hundred reflections can carry a large Δ*CC*₁⁄₂ without the dataset's headline +*CC*₁⁄₂ moving by anything like as much. `DELTA_CC_HALF_SE` is the standard error of a *CC*₁⁄₂ on that +many reflections, in the same units, and a Δ*CC*₁⁄₂ smaller than it says nothing. The same quantity is +logged per 10° batch over the whole sweep, on the same batches as the radiation-damage curve. It is +measured last, after the per-frame scale, the decay slope and the per-batch relative-*B* have been +fitted, so it is the cost of what the corrections could **not** remove. + +**What Δ*CC*₁⁄₂ cannot do**, because the report must not imply otherwise: + +- it says nothing about the **cause**: a shutter fault and a crystal that slipped have the identical + signature, both integrating background, so the cause comes from the `REASON` column and never from + the Δ*CC*₁⁄₂ itself; +- a **second lattice entering** is invisible to it: those spots were never integrated, so they are not + in the merged intensities it measures; +- a **centring drift that is pure attenuation** reads ≈ 0. That is the right answer, not a blind spot: + the data are weak but consistent, the σ's already say so, and their Δ*CC*₁⁄₂ is the evidence that + discarding them would cost completeness for nothing; +- it is attributed to the frame carrying a rocking event's **peak** partial, so a Δ*CC*₁⁄₂ on a range + narrower than one rocking event may be carried by a neighbouring frame; +- `loss_of_centring` needs ≥ 350° of sweep to be named at all. On a 90° sweep the same drift is still + *detected*, only as `crystal_out_of_beam` or `weak_diffraction` — "cause not determined" here means + this sweep cannot determine it, not that it is undeterminable. + +The same finding is written **per image** into the `_process.h5` as `/entry/MX/sweepQuality` and +`/entry/MX/frameDisposition`, when one is written — see +[HDF5](HDF5.md#entry-mx-spot-finding-and-indexing-cxi-style). ## Translational pseudo-symmetry diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 5e2a8ac7a..97cd7161d 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -742,6 +742,15 @@ const char *SweepQualityReasonText(SweepQualityReason reason) { return "unknown"; } +const char *FrameDispositionCode(FrameDisposition disposition) { + switch (disposition) { + case FrameDisposition::Merged: return "merged"; + case FrameDisposition::Downgraded: return "downgraded"; + case FrameDisposition::Rejected: return "rejected"; + } + return "unknown"; +} + namespace { // A quantity a run did not measure prints as a dash, not as "nan": CCref has nothing to compare // against unless a reference was given, CCanom needs both Bijvoet mates split in two, and a thin diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index c2da284bf..f33146976 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -51,7 +51,8 @@ struct MergeStatisticsShell { // Why a stretch of the sweep came out much weaker than the rest of the run. Report-only: nothing is // excluded on the strength of it. The codes are the field's own words - "crystal rotating out of the -// beam" (HKL-2000 manual), "loss of centring during crystal rotation" (autoPROC). +// beam" (HKL-2000 manual), "loss of centring during crystal rotation" (autoPROC). The reason names what +// was SEEN over a stretch; the disposition below says what became of its frames. enum class SweepQualityReason { NoDiffraction, // the range recorded essentially no diffraction from the indexed lattice CrystalOutOfBeam, // frames were lost: the range gets a per-image scale far less often than the run @@ -63,6 +64,19 @@ enum class SweepQualityReason { const char *SweepQualityReasonCode(SweepQualityReason reason); // machine-readable, e.g. "out_of_beam" const char *SweepQualityReasonText(SweepQualityReason reason); // for a sentence, e.g. "out of beam" +// What became of a frame's observations. The three codes partition the sweep - every processed image is +// exactly one of them - so the percentages add to 100 and "5% of the frames were useless" has a meaning. +// This is a LEDGER of what the scaling and merging already did, not a decision of its own. +enum class FrameDisposition { + Merged, // its observations are in the merged data at their own weight + Downgraded, // merged, but over a flagged stretch: carried at the reduced weight its own scale and + // sigmas give it, which for weak-but-consistent data is the only honest weight there is + Rejected // nothing of it reached the merge - an earlier per-frame guard dropped it, or it + // recorded nothing to drop in the first place +}; + +const char *FrameDispositionCode(FrameDisposition disposition); // "merged", "downgraded", "rejected" + struct SweepQualityRange { // Inclusive, in processed-image ordinals - the numbering of _image.dat and of every other // per-image array rugnux writes. With -s/--stride the source image is start + ordinal * stride. @@ -77,6 +91,13 @@ struct SweepQualityRange { float mean_relative_cc = 1.0f; // / run median float indexed_fraction = 1.0f; // frames in the range that got a per-image scale at all float relative_b = NAN; // mean of the radiation-damage monitor's per-batch B over the range + // What became of the range, and what the merged data say about keeping it: delta-CC1/2 is the + // overall CC1/2 with the range minus the CC1/2 without it, over the reflections the range touches + // (negative = keeping it makes the merged intensities worse), beside the standard error of a CC1/2 + // on that many reflections. NaN where the range carries too few reflections to be measured. + FrameDisposition disposition = FrameDisposition::Downgraded; + float delta_cc_half = NAN; + float delta_cc_half_se = NAN; }; // Sweep-quality diagnostic (rotation only). `measured` separates "the run is clean" from "this did not @@ -87,6 +108,22 @@ struct SweepQuality { float flux_peak_to_trough = 1.0f; // the incident-flux proxy, over the whole run float modulation_peak_to_trough = 1.0f; // depth of a DIAGNOSED once-per-revolution modulation (1 = none) std::vector ranges; + + // The disposition of the whole sweep. The three counts add up to the frame count; `rejected_deg` is + // the same rejection in degrees, because a percentage of FRAMES can be moved by re-slicing the same + // experiment and a percentage of the rotation cannot. Rejected = never reached the merge; nothing + // is rejected on the strength of the delta-CC1/2 below. + int frames_merged = 0; + int frames_downgraded = 0; + int frames_rejected = 0; + float rejected_deg = 0.0f; + std::vector frame_disposition; // FrameDisposition per processed image + + // What keeping each batch of the sweep costs the merged intensities (NaN where a batch carries too + // few reflections to be measured). Same batches as the radiation-damage monitor's relative-B curve, + // so the two read side by side. + std::vector delta_cc_half_batch; + float delta_cc_half_batch_deg = 0.0f; }; struct MergeStatistics { @@ -117,7 +154,8 @@ struct MergeStatistics { double radiation_damage_batch_deg = 0.0; // Stretches of the sweep over which the crystal delivered much less than the rest of the run - // (MeasureSweepQuality). Report-only - no observation is dropped because of it. + // (MeasureSweepQuality), what keeping each one costs (MeasureBatchDeltaCCHalf) and what became of + // every frame of the sweep. Report-only - no observation is dropped because of any of it. SweepQuality sweep_quality; // Diffraction anisotropy (AnalyzeAnisotropy): the anisotropy tensor, the diffraction limit along diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index f18e20aa8..8562ae472 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -69,6 +69,11 @@ namespace { // which keeps every frame. constexpr double SEARCH_MIN_SCALE_RATIO = 0.1; + // Rotation width per batch of the radiation-damage curve - and of the delta-CC1/2 ledger, which + // deliberately reuses it rather than choosing a granularity of its own. In degrees rather than + // images so that re-slicing the same experiment does not change the unit. + constexpr double MONITOR_BATCH_DEG = 10.0; + // --- Sweep-quality diagnostic (MeasureSweepQuality) --- // A stretch is reported only when BOTH per-frame channels are down: the scale (how much the crystal // diffracted) and the CC to merge (whether what it diffracted is still usable). The CC channel is what @@ -1492,7 +1497,6 @@ void RotationScaleMerge::MeasureRadiationDamageB(int n_groups) { const double osc_deg = gon ? std::fabs(gon->GetIncrement_deg()) : 0.0; if (!(osc_deg > 1e-6)) return; - constexpr double MONITOR_BATCH_DEG = 10.0; // rotation width per batch for the damage curve const int frames_per_batch = std::max(1, static_cast(std::lround(MONITOR_BATCH_DEG / osc_deg))); const int n_batch = std::max(1, (n_frames + frames_per_batch - 1) / frames_per_batch); if (n_batch < 2) // a single batch is a global Wilson-B (degenerate with scale) - nothing to monitor @@ -1901,6 +1905,181 @@ void RotationScaleMerge::MeasureSweepQuality(const std::vector &partial } } +namespace { + // --- delta-CC1/2 (MeasureBatchDeltaCCHalf) --- + constexpr int DELTA_CC_MIN_REFLECTIONS = 50; // fewer than this and the pair of CC1/2 means nothing + + // CC1/2 in the sigma-tau form (Assmann, Brehm & Diederichs (2016) J. Appl. Cryst. 49, 1021-1028): + // the variance of the unique reflections' mean intensities against the error variance of a half + // dataset's mean, accumulated one reflection at a time. No half-set split, so no random numbers - + // the same answer every run, and more precise than the split estimate. + struct SigmaTau { + double n = 0, sum_mean = 0, sum_mean2 = 0, sum_eps = 0; + void Add(double mean, double eps) { + n += 1; sum_mean += mean; sum_mean2 += mean * mean; sum_eps += eps; + } + [[nodiscard]] double CC() const { + if (n < 4) return NAN; + const double var_y = (sum_mean2 - sum_mean * sum_mean / n) / (n - 1); + const double eps = sum_eps / n; + return var_y + eps / 2 > 0.0 ? (var_y - eps / 2) / (var_y + eps / 2) : NAN; + } + }; + + // A reflection's mean intensity and the error variance of a half dataset's mean of it, from the + // OBSERVED scatter of its own observations. Observed, and not the error model's sigmas, on purpose: + // a batch of garbage claims the same sigmas as a good one, so an error-model estimate would read a + // batch that adds noise as a batch that adds precision, and the sign of the test would invert. + void ReflectionMoments(int n, double s1, double s2, double &mean, double &eps) { + mean = s1 / n; + eps = 2.0 * std::max(0.0, (s2 - s1 * s1 / n) / (n - 1)) / n; + } +} + +// What keeping each stretch of the sweep costs the merged intensities. REPORT-ONLY: no observation is +// dropped on the strength of it, exactly as for MeasureSweepQuality above. +// +// delta-CC1/2 of a group of images is the overall CC1/2 with the group minus the CC1/2 without it, +// evaluated over the reflections the group touches (Assmann, Brehm & Diederichs (2016) J. Appl. Cryst. +// 49, 1021-1028; Assmann, Wang & Diederichs (2020) Acta Cryst. D76, 636-652 - XDSCC12). Negative means +// keeping the group makes the merged data worse. It is the only instrument here that measures harm to +// the quantity the run actually ships: the channels that DETECT and NAME a bad stretch (the per-image +// scale, the per-image CC, the once-per-revolution harmonic, the relative-B curve) cannot say whether +// removing it would pay. What a bad stretch MEANS is never in this number; the reason codes name that. +// +// It runs LAST, on the corrected fulls, so it measures data the per-frame scale, the decay slope and +// the per-batch relative-B have already downgraded - harm the model has already removed must not be +// counted a second time as if it were still there. +// +// Two properties of the statistic set the rest of this, and both are the authors': its standard error +// goes as 1/sqrt(n_refl), so the curve is measured at the batch the scaling model already forms +// (MONITOR_BATCH_DEG) rather than per frame; and a group whose delta-CC1/2 is positive or near zero is +// not evidence of anything - "these may just be weak, and rejection without good reason may ultimately +// reduce the completeness". +void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { + sweep_quality.delta_cc_half_batch.clear(); + sweep_quality.delta_cc_half_batch_deg = 0.0f; + const auto gon = x.GetGoniometer(); + const double osc_deg = gon ? std::fabs(gon->GetIncrement_deg()) : 0.0; + const int n_full = static_cast(fulls.size()); + if (!(osc_deg > 1e-6) || n_full == 0 || n_groups <= 0) + return; + const int per_batch = std::max(1, static_cast(std::lround(MONITOR_BATCH_DEG / osc_deg))); + const int n_batch = (n_frames + per_batch - 1) / per_batch; + if (n_batch < 3) // with one or two batches "the rest of the run" is not a reference + return; + + // The observations, once: scaled intensity, ASU group (-1 = not in this merge, the same test the + // merge itself applies), and the frame order that makes any frame range a contiguous span. + std::vector obs_I(n_full); + std::vector obs_group(n_full), order(n_full), fstart(n_frames + 1, 0); + for (int i = 0; i < n_full; ++i) { + const Obs &o = fulls[i]; + obs_I[i] = o.I * o.corr; + obs_group[i] = UsableFull(o) ? o.group : -1; + ++fstart[o.frame + 1]; + } + for (int f = 0; f < n_frames; ++f) fstart[f + 1] += fstart[f]; + { + std::vector fill(fstart.begin(), fstart.end() - 1); + for (int i = 0; i < n_full; ++i) order[fill[fulls[i].frame]++] = i; + } + + // Per-reflection totals, and the per-reflection contribution of one frame range - which is all a + // leave-one-out needs, because a count, a sum and a sum of squares subtract. + std::vector tot_n(n_groups, 0), d_n(n_groups, 0), touched; + std::vector tot_s1(n_groups, 0.0), tot_s2(n_groups, 0.0), + d_s1(n_groups, 0.0), d_s2(n_groups, 0.0); + for (int i = 0; i < n_full; ++i) { + const int g = obs_group[i]; + if (g < 0) continue; + tot_n[g] += 1; tot_s1[g] += obs_I[i]; tot_s2[g] += obs_I[i] * obs_I[i]; + } + + // delta-CC1/2 of the frames [f0, f1], and its standard error. `se` is the standard error of a + // SINGLE CC1/2 on that many reflections - single, not of the difference: the two CC1/2 share most + // of their data, so their difference is the more precise of the two and this errs toward not + // claiming a harm the data cannot resolve. + double delta_cc = 0.0, se = 0.0; + auto measure = [&](int f0, int f1) { + touched.clear(); + for (int f = f0; f <= f1; ++f) + for (int j = fstart[f]; j < fstart[f + 1]; ++j) { + const int i = order[j], g = obs_group[i]; + if (g < 0) continue; + if (d_n[g] == 0) touched.push_back(g); + d_n[g] += 1; d_s1[g] += obs_I[i]; d_s2[g] += obs_I[i] * obs_I[i]; + } + // Only reflections that are still measurable WITHOUT the range count, which is the published + // statistic's own restriction: a range holding the only observations of the reflections it + // touches contributes to neither CC1/2 and reads ~0. + SigmaTau with, without; + for (int g : touched) { + const int n_without = tot_n[g] - d_n[g]; + if (tot_n[g] >= 2 && n_without >= 2) { + double mean, eps; + ReflectionMoments(tot_n[g], tot_s1[g], tot_s2[g], mean, eps); + with.Add(mean, eps); + ReflectionMoments(n_without, tot_s1[g] - d_s1[g], tot_s2[g] - d_s2[g], mean, eps); + without.Add(mean, eps); + } + d_n[g] = 0; d_s1[g] = 0.0; d_s2[g] = 0.0; + } + const double cc_with = with.CC(), cc_without = without.CC(); + if (with.n < DELTA_CC_MIN_REFLECTIONS || !std::isfinite(cc_with) || !std::isfinite(cc_without)) + return false; + delta_cc = cc_with - cc_without; + se = (1.0 - cc_with * cc_with) / std::sqrt(with.n - 3); + return true; + }; + + // The curve over the whole sweep, on the same batches as the radiation-damage monitor so the two + // read side by side, and then what keeping each flagged stretch costs. + sweep_quality.delta_cc_half_batch.assign(n_batch, NAN); + sweep_quality.delta_cc_half_batch_deg = static_cast(MONITOR_BATCH_DEG); + for (int b = 0; b < n_batch; ++b) + if (measure(b * per_batch, std::min(n_frames, (b + 1) * per_batch) - 1)) + sweep_quality.delta_cc_half_batch[b] = static_cast(delta_cc); + for (auto &r : sweep_quality.ranges) + if (measure(r.first_image, r.last_image)) { + r.delta_cc_half = static_cast(delta_cc); + r.delta_cc_half_se = static_cast(se); + } +} + +// Give every frame and every range its disposition, and count the sweep. A frame is merged if its +// observations reached the merge and rejected if none did - because an earlier per-frame guard dropped +// it, or because it recorded nothing to drop in the first place, which are the same thing to a user +// asking how much of the experiment was useless. Of the merged frames, the ones inside a flagged range +// are downgraded: they are in the merge, carried at the reduced weight their own scale and sigmas give +// them, which is the only honest weight weak data can have. +void RotationScaleMerge::CountSweepDisposition() { + const auto gon = x.GetGoniometer(); + const double osc_deg = gon ? std::fabs(gon->GetIncrement_deg()) : 0.0; + auto &sq = sweep_quality; + sq.frame_disposition.assign(n_frames, static_cast(FrameDisposition::Merged)); + for (int f = 0; f < n_frames; ++f) + if (!frame_in_merge[f]) + sq.frame_disposition[f] = static_cast(FrameDisposition::Rejected); + for (auto &r : sq.ranges) { + bool all_out = true; + for (int f = r.first_image; f <= r.last_image; ++f) + if (sq.frame_disposition[f] != static_cast(FrameDisposition::Rejected)) { + sq.frame_disposition[f] = static_cast(FrameDisposition::Downgraded); + all_out = false; + } + r.disposition = all_out ? FrameDisposition::Rejected : FrameDisposition::Downgraded; + } + sq.frames_merged = sq.frames_downgraded = sq.frames_rejected = 0; + for (int f = 0; f < n_frames; ++f) + switch (static_cast(sq.frame_disposition[f])) { + case FrameDisposition::Merged: ++sq.frames_merged; break; + case FrameDisposition::Downgraded: ++sq.frames_downgraded; break; + case FrameDisposition::Rejected: ++sq.frames_rejected; break; + } + sq.rejected_deg = static_cast(sq.frames_rejected * osc_deg); +} + void RotationScaleMerge::RefineRelativeB(int n_groups) { // RefineDecay removes the AVERAGE radiation-damage falloff as a single global relative-B slope, but the // relative scattering power drifts NON-monotonically across a run (absorption path as the crystal @@ -2506,6 +2685,7 @@ bool RotationScaleMerge::DropCollapsedFullScales(bool from_staging) { for (int f = 0; f < n_frames; ++f) if (std::isfinite(g_frame[f]) && g_frame[f] < g_floor) { worst = std::max(worst, fitted_median / g_frame[f]); + frame_in_merge[f] = 0; // out of the merge, so out of the disposition ledger's denominator ++n_frames_dropped; } // With no frame below the floor there is nothing for the walk below to zero, and it is a pass over @@ -2528,6 +2708,15 @@ bool RotationScaleMerge::DropCollapsedFullScales(bool from_staging) { return n_dropped > 0; } +bool RotationScaleMerge::UsableFull(const Obs &o) const { + if (o.group < 0) return false; + if (!frame_cell_ok[o.frame]) return false; + if (!(o.corr > 0.0f) || !std::isfinite(o.corr)) return false; + if (o.partiality < min_partiality) return false; + const float I_corr = o.I * o.corr, sigma_corr = o.sigma * o.corr; + return std::isfinite(I_corr) && std::isfinite(sigma_corr) && sigma_corr > 0.0f; +} + void RotationScaleMerge::Combine() { fulls.clear(); g_full.assign(n_frames, 1.0); @@ -2824,15 +3013,9 @@ namespace { RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool for_search, bool fulls_resident, bool full_stats) { // A full is usable for the merge / error model if it passes AddImage's filters (with the current - // ice context). group >= 0 already encodes "not absent and passes AcceptReflection". + // ice context). auto usable_merge = [&](const Obs &o) { - if (o.group < 0) return false; - if (!frame_cell_ok[o.frame]) return false; - if (!(o.corr > 0.0f) || !std::isfinite(o.corr)) return false; - if (for_search && o.on_ice) return false; - if (o.partiality < min_partiality) return false; - const float I_corr = o.I * o.corr, sigma_corr = o.sigma * o.corr; - return std::isfinite(I_corr) && std::isfinite(sigma_corr) && sigma_corr > 0.0f; + return UsableFull(o) && !(for_search && o.on_ice); }; // The em-stats / samples / merge-accumulate / R_meas reductions run on the resident, scaled fulls @@ -3749,6 +3932,10 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st } } const std::vector partial_scaled = frame_scaled_scratch; + // The disposition ledger's denominator; the per-frame drops below take frames out of it. + frame_in_merge.assign(n_frames, 0); + for (int f = 0; f < n_frames; ++f) + frame_in_merge[f] = partial_scaled[f] && frame_cell_ok[f] ? 1 : 0; // --- 1b. Guard the per-frame partial scales. Unconditional: the smooth-G window below used to be // where this lived, so --smooth-g 0, a dataset with no oscillation width, and every caller @@ -3757,6 +3944,8 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st std::vector apply; std::vector ratio; if (DropCollapsedScales(frame_scaled_scratch, g_partial, apply, ratio)) { + for (int f = 0; f < n_frames; ++f) + if (apply[f]) frame_in_merge[f] = 0; bool applied_on_gpu = false; #ifdef JFJOCH_USE_CUDA if (gpu_active_) { @@ -3920,6 +4109,8 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st ++n_rejected; } if (n_rejected > 0) { + for (int f = 0; f < n_frames; ++f) + if (reject[f]) frame_in_merge[f] = 0; bool reject_on_gpu = false; #ifdef JFJOCH_USE_CUDA if (gpu_active_) { @@ -4050,7 +4241,8 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st if (refine_surfaces) MeasureRadiationDamageB(n_groups); // Sweep-quality diagnostic, on the per-frame scale the partial scaling just fitted (the flux is - // already out of it) and the per-frame CC computed above. Report-only; drops nothing. + // already out of it) and the per-frame CC computed above. It drops nothing itself; the ranges it + // segments are the ledger 4c below decides on. if (refine_surfaces) MeasureSweepQuality(partial_scaled, cc, cc_n); // CC1/2 of these same fulls with the surfaces not yet folded in, for the caller's cross-pass quality @@ -4081,6 +4273,15 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st RefineModulation(modulation_iter, n_groups); if (refine_surfaces && absorption_iter > 0) RefineAbsorptionTime(absorption_iter, n_groups); // last: the static surfaces get first claim + + // --- 4c. What keeping each stretch of the sweep costs the merged intensities (delta-CC1/2), and + // what became of every frame of the run. LAST, on the corrected fulls: harm the corrections + // have already removed must not be counted a second time as if it were still there. + // Report-only, like the two diagnostics above - nothing here changes an observation. --- + if (refine_surfaces && sweep_quality.measured) { + MeasureBatchDeltaCCHalf(n_groups); + CountSweepDisposition(); + } #ifdef JFJOCH_USE_CUDA // The corrections mutate the host fulls' corr; when the merge runs on the resident (GPU) fulls, push // the corrected corr back to the device so the merge reads it. diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 1d7879617..52cd15e36 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -303,10 +303,16 @@ private: std::vector rad_damage_b_batch; // per-batch relative-B curve (A^2) double rad_damage_batch_deg = 0.0; // rotation width per batch (deg) - // Sweep-quality diagnostic (MeasureSweepQuality; report-only, copied into the result statistics by - // MergeAndStats). Empty and not measured until it runs. + // Sweep-quality diagnostic (MeasureSweepQuality, then the delta-CC1/2 and the disposition; + // copied into the result statistics by MergeAndStats). Empty and not measured until it runs. SweepQuality sweep_quality; + // Frames whose observations reach this pass's merge at all (1 = in). Built by Run(), where every + // per-frame drop is decided. A full cannot answer this: the combine attributes a rocking event to + // its PEAK frame, so on fine slicing most frames own no full at all while their measurements sit + // inside one. + std::vector frame_in_merge; + // Working per-group arrays (sized to the current group count; reused). std::vector group_h, group_k, group_l; @@ -377,6 +383,11 @@ private: void Combine(); // partials -> fulls (CPU) + // Is this full in the merge? group >= 0 already encodes "not absent and passes AcceptReflection"; + // the rest is the frame's cell-consistency mask, a usable scale and a usable sigma. The P1 search + // pass adds its own ice test on top of this (see MergeAndStats). + [[nodiscard]] bool UsableFull(const Obs &o) const; + // Drop the fulls of any frame whose scale collapsed toward zero. The fulls are scaled with the Unity // model, so their corr IS 1/G and a collapsed G multiplies every intensity on that frame without // bound. Covers the CPU and GPU scaling paths alike; `from_staging` says the fulls' frame and corr @@ -420,11 +431,17 @@ private: // any decay correction and store the first->last relative-B change + the per-batch curve on this object // (copied into the result statistics by MergeAndStats, then printed / logged / written to the mmCIF). void MeasureRadiationDamageB(int n_groups); - // Sweep-quality diagnostic (report-only): find the contiguous stretches of the sweep over which the - // crystal delivered much less than the rest of the run, and say what each one looks like. Reads the + // Sweep-quality diagnostic: find the contiguous stretches of the sweep over which the crystal + // delivered much less than the rest of the run, and say what each one looks like. Reads the // per-frame scale (with the incident flux already divided out) and the per-frame CC to merge. void MeasureSweepQuality(const std::vector &partial_scaled, const std::vector &cc, const std::vector &cc_n); + // Per-batch delta-CC1/2 on the corrected fulls: what keeping each stretch of the sweep costs the + // merged intensities. REPORT-ONLY, like MeasureSweepQuality - no observation is dropped on the + // strength of it. See the .cpp for the statistic. + void MeasureBatchDeltaCCHalf(int n_groups); + // Give every frame and every flagged range its disposition, and count the sweep. + void CountSweepDisposition(); // Per-batch relative-B, applied after RefineDecay: the single decay slope removes the average // radiation-damage falloff, but the relative scattering power drifts NON-monotonically across a run // (absorption path, crystal slippage, dose bursts). Refine one relative Debye-Waller B per batch diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index fb6f67488..e65edee4e 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -596,6 +596,14 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen master_file->ReadElement("/entry/MX/sweepQualityReasons", i) .value_or("")); } + dataset->frame_disposition = master_file->ReadOptVector("/entry/MX/frameDisposition"); + if (master_file->Exists("/entry/MX/frameDispositionCodes")) { + const auto dim = master_file->GetDimension("/entry/MX/frameDispositionCodes"); + for (size_t i = 0; i < (dim.empty() ? 0 : dim[0]); i++) + dataset->frame_disposition_codes.push_back( + master_file->ReadElement("/entry/MX/frameDispositionCodes", i) + .value_or("")); + } } if (master_file->Exists("/entry/image")) dataset->max_value = master_file->ReadOptVector("/entry/image/max_value"); diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index cc50376eb..9d82b2b5a 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -76,6 +76,12 @@ struct JFJochReaderDataset { std::vector sweep_quality; std::vector sweep_quality_reasons; + // Per-image disposition from /entry/MX/frameDisposition: a 0-based index into + // frame_disposition_codes (merged, downgraded, rejected). What was DONE with the image, where + // sweep_quality above says what was SEEN. + std::vector frame_disposition; + std::vector frame_disposition_codes; + // Maps this dataset's image index -> the original/collected image number it came from. // Empty means identity (image i == original image i). Lets a dataset be a subset (or strided // selection) of the truly collected images: reprocessing snapshots over a sub-range, and (in diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index cbb70b313..c23ec23d1 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1011,6 +1011,16 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, // ---- sweep quality { const auto &sq = result.merge_statistics.sweep_quality; + const int n_sweep_frames = sq.frames_merged + sq.frames_downgraded + sq.frames_rejected; + auto pct = [&](int n) { return n_sweep_frames > 0 ? 100.0 * n / n_sweep_frames : 0.0; }; + // A delta-CC1/2 or its standard error that could not be measured prints as a dash - a + // range too thin to judge is not a range that was judged and came out at zero. + auto cc_cell = [](float v) { + return std::isfinite(v) ? fmt::format("{:+.4f}", v) : std::string("-"); + }; + auto se_cell = [](float v) { + return std::isfinite(v) ? fmt::format("{:.4f}", v) : std::string("-"); + }; Add(s, KeyInt("SWEEP_QUALITY_COUNT", static_cast(sq.ranges.size()))); Add(s, KeyEnum("SWEEP_QUALITY_STATUS", sq.measured ? "COMPUTED" : "NOT_COMPUTED", {"COMPUTED", "NOT_COMPUTED"}, true)); @@ -1020,14 +1030,35 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, codes += (codes.empty() ? "" : " ") + std::string(SweepQualityReasonCode(static_cast(r))); Add(s, KeyText("SWEEP_QUALITY_REASONS", codes, true)); + std::string dispositions; + for (int d = 0; d <= static_cast(FrameDisposition::Rejected); ++d) + dispositions += (dispositions.empty() ? "" : " ") + + std::string(FrameDispositionCode(static_cast(d))); + Add(s, KeyText("SWEEP_DISPOSITIONS", dispositions, true)); } if (sq.measured) { + Add(s, KeyInt("FRAMES_MERGED", sq.frames_merged)); + Add(s, KeyInt("FRAMES_DOWNGRADED", sq.frames_downgraded)); + Add(s, KeyInt("FRAMES_REJECTED", sq.frames_rejected)); + Add(s, KeyReal("FRAMES_REJECTED_PCT", pct(sq.frames_rejected), "{:.2f}")); + Add(s, KeyReal("ROTATION_REJECTED_DEG", sq.rejected_deg, "{:.1f}")); Add(s, KeyReal("SWEEP_ROTATION", sq.sweep_deg, "{:.1f}", true)); Add(s, KeyReal("FLUX_PEAK_TO_TROUGH", sq.flux_peak_to_trough, "{:.2f}", true)); Add(s, KeyReal("SCALE_MODULATION_PEAK_TO_TROUGH", sq.modulation_peak_to_trough, "{:.2f}", true)); } Add(s, Blank()); + if (sq.measured) + // Frames AND degrees, because a percentage of frames can be moved by re-slicing the + // same experiment - 5% of 0.05 deg frames is not 5% of the experiment - and a percentage + // of the rotation cannot. + Add(s, Prose(fmt::format( + " Sweep: {} frames / {:.1f} deg; {} merged ({:.1f}%), {} downgraded ({:.1f}%),\n" + " {} rejected ({:.1f}% of frames, {:.1f} deg).", + n_sweep_frames, sq.sweep_deg, + sq.frames_merged, pct(sq.frames_merged), + sq.frames_downgraded, pct(sq.frames_downgraded), + sq.frames_rejected, pct(sq.frames_rejected), sq.rejected_deg))); if (sq.ranges.empty()) { // The empty table used to print its header and two rules around nothing on every // clean run - a table that says "no problem" by being blank. @@ -1039,20 +1070,29 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, Prose(" Stretches of the sweep over which the crystal delivered much less than the rest of\n" " the run. SEVERITY is the fraction of the run's typical diffracting power missing over\n" " the range; SCALE and CC are relative to the run median; INDEXED is the fraction of the\n" - " range's frames that were scaled at all. Nothing is excluded on the strength of this.\n")); + " range's frames that were scaled at all. DELTA_CC_HALF is what keeping the range costs\n" + " the merged intensities - the overall CC1/2 with it minus the CC1/2 without it, over the\n" + " reflections it touches - beside the standard error of a CC1/2 on that many reflections.\n" + " DISPOSITION says what became of the range: REJECTED means nothing of it reached the\n" + " merge - the frames were never scaled, or an earlier guard dropped them - and DOWNGRADED\n" + " means it is in the merge at the reduced weight its own scale and sigmas give it. Nothing\n" + " is excluded on the strength of DELTA_CC_HALF: it is a measurement, not a filter.\n")); ReportEntry e; e.kind = ReportEntry::Kind::Table; e.table.columns = {"FIRST_IMAGE", "LAST_IMAGE", "N_IMAGES", "ROTATION", "REASON", - "SEVERITY", "SCALE", "CC", "INDEXED"}; + "SEVERITY", "SCALE", "CC", "INDEXED", "DISPOSITION", + "DELTA_CC_HALF", "DELTA_CC_HALF_SE"}; e.table.text_header = - " FIRST_IMAGE LAST_IMAGE N_IMAGES ROTATION REASON SEVERITY SCALE CC INDEXED\n" - " ----------- ----------- --------- -------- -------------------- -------- ------ ------ --------"; + " FIRST_IMAGE LAST_IMAGE N_IMAGES ROTATION REASON SEVERITY SCALE CC INDEXED DISPOSITION DELTA_CC_HALF DELTA_CC_HALF_SE\n" + " ----------- ----------- --------- -------- -------------------- -------- ------ ------ -------- ----------- ------------- ----------------"; for (const auto &r : sq.ranges) { e.table.text_rows.push_back(fmt::format( - " {:11d} {:11d} {:9d} {:8.1f} {:<20} {:8.2f} {:6.2f} {:6.2f} {:8.2f}", + " {:11d} {:11d} {:9d} {:8.1f} {:<20} {:8.2f} {:6.2f} {:6.2f} {:8.2f} {:<11} {:>13} {:>16}", r.first_image, r.last_image, r.last_image - r.first_image + 1, r.rotation_deg, SweepQualityReasonCode(r.reason), r.severity, - r.mean_relative_scale, r.mean_relative_cc, r.indexed_fraction)); + r.mean_relative_scale, r.mean_relative_cc, r.indexed_fraction, + FrameDispositionCode(r.disposition), cc_cell(r.delta_cc_half), + se_cell(r.delta_cc_half_se))); std::vector row; row.push_back(KeyInt("", r.first_image).value); row.push_back(KeyInt("", r.last_image).value); @@ -1063,12 +1103,25 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, row.push_back(KeyReal("", r.mean_relative_scale, "{:.2f}").value); row.push_back(KeyReal("", r.mean_relative_cc, "{:.2f}").value); row.push_back(KeyReal("", r.indexed_fraction, "{:.2f}").value); + row.push_back(KeyText("", FrameDispositionCode(r.disposition)).value); + row.push_back(KeyText("", cc_cell(r.delta_cc_half)).value); + row.push_back(KeyText("", se_cell(r.delta_cc_half_se)).value); e.table.cells.push_back(std::move(row)); + std::string cost; + if (r.disposition == FrameDisposition::Rejected) + cost = " - nothing of them reached the merge"; + else if (std::isfinite(r.delta_cc_half)) + cost = fmt::format(" - kept at their own reduced weight; removing them would change " + "CC1/2 by {:+.4f} +/- {:.4f} over the reflections they touch", + -r.delta_cc_half, r.delta_cc_half_se); + else + cost = " - kept at their own reduced weight; too few reflections to measure what " + "they cost"; Warn(doc, PathologyCode::SWEEP_GAPS, fmt::format( - "Frames {}-{} {} ({:.1f} deg, scale {:.2f} and CC {:.2f} of the run, {:.0f}% scaled)", + "Frames {}-{} {} ({:.1f} deg, scale {:.2f} and CC {:.2f} of the run, {:.0f}% scaled){}", r.first_image, r.last_image, SweepQualityReasonText(r.reason), r.rotation_deg, r.mean_relative_scale, r.mean_relative_cc, - 100.0 * r.indexed_fraction)); + 100.0 * r.indexed_fraction, cost)); } Add(s, std::move(e)); } @@ -1574,11 +1627,17 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, if (std::isfinite(db)) row("Radiation damage", fmt::format("relative B {:+.2f} A^2 over the sweep", db)); const auto &sq = result.merge_statistics.sweep_quality; - if (sq.measured) - row("Sweep", sq.ranges.empty() - ? fmt::format("no degraded ranges in {:.1f} deg", sq.sweep_deg) - : fmt::format("{} degraded range(s) in {:.1f} deg", sq.ranges.size(), - sq.sweep_deg)); + if (sq.measured) { + const int n_sweep = sq.frames_merged + sq.frames_downgraded + sq.frames_rejected; + row("Sweep", sq.frames_rejected > 0 + ? fmt::format("{} of {} frames rejected ({:.1f} deg of {:.1f}), " + "{} degraded range(s)", sq.frames_rejected, n_sweep, + sq.rejected_deg, sq.sweep_deg, sq.ranges.size()) + : (sq.ranges.empty() + ? fmt::format("no degraded ranges in {:.1f} deg", sq.sweep_deg) + : fmt::format("{} degraded range(s) in {:.1f} deg, none rejected", + sq.ranges.size(), sq.sweep_deg))); + } } if (!f.empty()) { f.pop_back(); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index d1673031c..1fc8b52cc 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -5988,16 +5988,26 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // Per-image form of the sweep-quality ranges, for the _process.h5: one code per image, 0 where // the image is in no flagged range, plus the vocabulary the codes index. Only filled when the - // diagnostic ran, so absent datasets mean "not looked for" rather than "all clean". + // diagnostic ran, so absent datasets mean "not looked for" rather than "all clean". The + // disposition rides alongside as a second per-image array: what the reason SAW and what became + // of the frame are different statements, and a consumer needs both. if (sm.statistics.sweep_quality.measured) { + const auto &sq = sm.statistics.sweep_quality; end_msg.sweep_quality.assign(end_msg.max_image_number, 0); - for (const auto &r : sm.statistics.sweep_quality.ranges) + for (const auto &r : sq.ranges) for (int64_t i = std::max(0, r.first_image); i <= r.last_image && i < static_cast(end_msg.sweep_quality.size()); ++i) end_msg.sweep_quality[i] = static_cast(r.reason) + 1; for (int r = 0; r <= static_cast(SweepQualityReason::RadiationDamage); ++r) end_msg.sweep_quality_reasons.emplace_back( SweepQualityReasonCode(static_cast(r))); + end_msg.frame_disposition.assign( + sq.frame_disposition.begin(), + sq.frame_disposition.begin() + + std::min(end_msg.max_image_number, sq.frame_disposition.size())); + for (int d = 0; d <= static_cast(FrameDisposition::Rejected); ++d) + end_msg.frame_disposition_codes.emplace_back( + FrameDispositionCode(static_cast(d))); } { // Stride rather than take the head: the merged list is ordered by hkl, so the first N @@ -6086,30 +6096,48 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } // Sweep-quality report: the stretches of the sweep over which the crystal delivered much - // less than the rest of the run, and what each one looks like. Nothing is excluded because - // of it - the frames still carry signal, and this is a message for the beamline, not a - // filter. Frame numbers are processed-image ordinals, inclusive, as in _image.dat. + // less than the rest of the run, what each one looks like, what became of it and what + // keeping it costs the merged intensities. Nothing is excluded because of it - the frames + // still carry signal, and this is a message for the beamline, not a filter. Frame numbers + // are processed-image ordinals, inclusive, as in _image.dat. const auto &sq = sm.statistics.sweep_quality; if (sq.measured) { + const int n_sweep = sq.frames_merged + sq.frames_downgraded + sq.frames_rejected; + auto pct = [&](int n) { return n_sweep > 0 ? 100.0 * n / n_sweep : 0.0; }; std::ostringstream os; os << fmt::format("Sweep quality over {:.0f} deg (per-image scale with the incident flux, " "which varied {:.2f}x, already divided out):\n", sq.sweep_deg, sq.flux_peak_to_trough); + os << fmt::format(" {} frames: {} merged ({:.1f}%), {} downgraded ({:.1f}%), " + "{} rejected ({:.1f}%, {:.1f} deg)\n", + n_sweep, sq.frames_merged, pct(sq.frames_merged), + sq.frames_downgraded, pct(sq.frames_downgraded), + sq.frames_rejected, pct(sq.frames_rejected), sq.rejected_deg); if (sq.ranges.empty()) { os << " no stretch of the sweep is materially worse than the run"; } else { - os << " frames rotation diagnosis severity scale CC indexed\n"; + os << " frames rotation diagnosis severity scale CC indexed disposition dCC1/2\n"; for (const auto &r : sq.ranges) - os << fmt::format(" {:<17s} {:6.1f} deg {:<16s} {:5.2f} {:5.2f} {:5.2f} {:4.0f}%\n", + os << fmt::format(" {:<17s} {:6.1f} deg {:<16s} {:5.2f} {:5.2f} {:5.2f} {:4.0f}% {:<11s} {}\n", fmt::format("{}-{}", r.first_image, r.last_image), r.rotation_deg, SweepQualityReasonText(r.reason), r.severity, r.mean_relative_scale, - r.mean_relative_cc, 100.0 * r.indexed_fraction); + r.mean_relative_cc, 100.0 * r.indexed_fraction, + FrameDispositionCode(r.disposition), + std::isfinite(r.delta_cc_half) + ? fmt::format("{:+.4f} +/- {:.4f}", r.delta_cc_half, + r.delta_cc_half_se) + : std::string("-")); os << " => severity is the fraction of the run's typical diffracting power missing " - "over the range"; + "over the range; dCC1/2 is what keeping it costs the merged intensities"; if (sq.modulation_peak_to_trough >= 1.05f) os << fmt::format("\n => once-per-revolution modulation of the per-image scale: " "{:.1f}x peak to trough", sq.modulation_peak_to_trough); } + if (!sq.delta_cc_half_batch.empty()) { + os << fmt::format("\n => delta-CC1/2 per {:.0f} deg:", sq.delta_cc_half_batch_deg); + for (float d : sq.delta_cc_half_batch) + os << (std::isfinite(d) ? fmt::format(" {:+.3f}", d) : std::string(" -")); + } logger.Info("{}", os.str()); } } diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index f71b675d1..4f5821e64 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -20,12 +20,20 @@ namespace { .reason = SweepQualityReason::CrystalOutOfBeam, .severity = 0.83f, .rotation_deg = 5.0f, .mean_relative_scale = 0.12f, .mean_relative_cc = 0.30f, - .indexed_fraction = 0.015f}); + .indexed_fraction = 0.015f, + .disposition = FrameDisposition::Rejected, + .delta_cc_half = -0.0312f, .delta_cc_half_se = 0.0060f}); out.ranges.push_back(SweepQualityRange{.first_image = 400, .last_image = 499, .reason = SweepQualityReason::RadiationDamage, .severity = 0.41f, .rotation_deg = 10.0f, .mean_relative_scale = 0.55f, .mean_relative_cc = 0.80f, - .indexed_fraction = 0.62f}); + .indexed_fraction = 0.62f, + .disposition = FrameDisposition::Downgraded, + .delta_cc_half = 0.0011f, .delta_cc_half_se = 0.0042f}); + out.frames_merged = 450; + out.frames_downgraded = 100; + out.frames_rejected = 50; + out.rejected_deg = 5.0f; return out; } @@ -35,6 +43,13 @@ namespace { out.emplace_back(SweepQualityReasonCode(static_cast(r))); return out; } + + std::vector DispositionVocabulary() { + std::vector out; + for (int d = 0; d <= static_cast(FrameDisposition::Rejected); ++d) + out.emplace_back(FrameDispositionCode(static_cast(d))); + return out; + } } TEST_CASE("SweepQuality_ReasonVocabulary", "[Diagnostics]") { @@ -48,6 +63,12 @@ TEST_CASE("SweepQuality_ReasonVocabulary", "[Diagnostics]") { CHECK(codes[2] == "weak_diffraction"); CHECK(codes[3] == "loss_of_centring"); CHECK(codes[4] == "radiation_damage"); + + // The disposition vocabulary is an interface for the same reasons, and the per-image HDF5 array + // indexes it from 0. + CHECK(std::string(FrameDispositionCode(FrameDisposition::Merged)) == "merged"); + CHECK(std::string(FrameDispositionCode(FrameDisposition::Downgraded)) == "downgraded"); + CHECK(std::string(FrameDispositionCode(FrameDisposition::Rejected)) == "rejected"); } TEST_CASE("ResultReport_Render", "[Diagnostics]") { @@ -100,12 +121,26 @@ TEST_CASE("ResultReport_Render", "[Diagnostics]") { CHECK(dev_text.find("\nSWEEP_QUALITY_STATUS= COMPUTED\n") != std::string::npos); CHECK(dev_text.find("\nSWEEP_QUALITY_REASONS= no_diffraction crystal_out_of_beam weak_diffraction " "loss_of_centring radiation_damage\n") != std::string::npos); + CHECK(dev_text.find("\nSWEEP_DISPOSITIONS= merged downgraded rejected\n") != std::string::npos); + + // The disposition headline: frames AND degrees, so that re-slicing cannot move the percentage + // without the degrees moving with it. + CHECK(text.find("\nFRAMES_MERGED= 450\n") != std::string::npos); + CHECK(text.find("\nFRAMES_DOWNGRADED= 100\n") != std::string::npos); + CHECK(text.find("\nFRAMES_REJECTED= 50\n") != std::string::npos); + CHECK(text.find("\nFRAMES_REJECTED_PCT= 8.33\n") != std::string::npos); + CHECK(text.find("\nROTATION_REJECTED_DEG= 5.0\n") != std::string::npos); + CHECK(text.find("600 frames / 60.0 deg; 450 merged (75.0%), 100 downgraded (16.7%)") + != std::string::npos); // One table row per range, with the reason code verbatim. CHECK(text.find(" 100 149 50 5.0 crystal_out_of_beam ") != std::string::npos); CHECK(text.find(" 400 499 100 10.0 radiation_damage ") != std::string::npos); + // ... and its disposition beside what keeping it costs the merged intensities. + CHECK(text.find("rejected -0.0312 0.0060") != std::string::npos); + CHECK(text.find("downgraded +0.0011 0.0042") != std::string::npos); // ... and one plain-English WARNING line per range, greppable by the marker alone. They sit in // the SUMMARY at the top of the file now, which is what the verdict is composed from. @@ -115,6 +150,9 @@ TEST_CASE("ResultReport_Render", "[Diagnostics]") { CHECK(text.find("VERDICT=") < text.find("SWEEP_QUALITY_COUNT=")); CHECK(text.find("\n# WARNING: Frames 100-149 out of beam (5.0 deg,") != std::string::npos); CHECK(text.find("\n# WARNING: Frames 400-499 radiation damage (10.0 deg,") != std::string::npos); + CHECK(text.find("nothing of them reached the merge") != std::string::npos); + CHECK(text.find("kept at their own reduced weight; removing them would change CC1/2 by -0.0011") + != std::string::npos); } TEST_CASE("ResultReport_RenderEmpty", "[Diagnostics]") { @@ -184,6 +222,15 @@ TEST_CASE("SweepQuality_HDF5RoundTrip", "[HDF5][Full][Diagnostics]") { static_cast(SweepQualityReason::CrystalOutOfBeam) + 1, 0, static_cast(SweepQualityReason::RadiationDamage) + 1}; + // What was DONE about them: the out-of-beam pair is dropped, the damaged tail stays in at its own + // reduced weight, the rest merges. Indexed from 0, unlike the reason codes above. + const std::vector expected_disposition{ + static_cast(FrameDisposition::Merged), + static_cast(FrameDisposition::Merged), + static_cast(FrameDisposition::Rejected), + static_cast(FrameDisposition::Rejected), + static_cast(FrameDisposition::Merged), + static_cast(FrameDisposition::Downgraded)}; { RegisterHDF5Filter(); @@ -194,6 +241,8 @@ TEST_CASE("SweepQuality_HDF5RoundTrip", "[HDF5][Full][Diagnostics]") { end_message.max_image_number = x.GetImageNum(); end_message.sweep_quality = expected; end_message.sweep_quality_reasons = ReasonVocabulary(); + end_message.frame_disposition = expected_disposition; + end_message.frame_disposition_codes = DispositionVocabulary(); FileWriter writer(start_message); std::vector image(x.GetPixelsNum(), 42); @@ -214,6 +263,8 @@ TEST_CASE("SweepQuality_HDF5RoundTrip", "[HDF5][Full][Diagnostics]") { CHECK(dataset->sweep_quality == expected); // The vocabulary travels with the codes, so a consumer can name them without this source. CHECK(dataset->sweep_quality_reasons == ReasonVocabulary()); + CHECK(dataset->frame_disposition == expected_disposition); + CHECK(dataset->frame_disposition_codes == DispositionVocabulary()); } REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 30eb8d1b8..5416cb5e8 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -1285,4 +1285,10 @@ void NXmx::EndResultVectors(const EndMessage &end) { if (!hdf5_file->Exists("/entry/MX/sweepQualityReasons")) hdf5_file->SaveVector("/entry/MX/sweepQualityReasons", end.sweep_quality_reasons); } + // Per-image disposition, same shape: frameDisposition[i] indexes frameDispositionCodes from 0. + if (!end.frame_disposition.empty() && !end.frame_disposition_codes.empty()) { + SaveVectorIfMissing(*hdf5_file, "/entry/MX/frameDisposition", end.frame_disposition); + if (!hdf5_file->Exists("/entry/MX/frameDispositionCodes")) + hdf5_file->SaveVector("/entry/MX/frameDispositionCodes", end.frame_disposition_codes); + } }