diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 1deaa9b4..27e47e7f 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -33,8 +33,8 @@ packages is linked or vendored, with the single exception of GEMMI (see **[XDS](https://xds.mr.mpg.de/)** — rotation geometry and notation, the reciprocal Lorentz and partiality treatment, the maximum-likelihood mosaicity estimate, the `MINPK` criterion for rejecting -a reflection whose predicted profile is not cleanly its own, and the intensity-based test for a -centred lattice. W. Kabsch, "XDS" (2010), Acta Cryst. D66, 125-132 +a reflection whose predicted profile is not cleanly its own, the intensity-based test for a +centred lattice, and the scaling correction surfaces indexed by image number and detector region. W. Kabsch, "XDS" (2010), Acta Cryst. D66, 125-132 [doi:10.1107/S0907444909047337](https://doi.org/10.1107/S0907444909047337); W. Kabsch, "Integration, scaling, space-group assignment and post-refinement" (2010), Acta Cryst. D66, 133-144 [doi:10.1107/S0907444909047374](https://doi.org/10.1107/S0907444909047374). diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 50c4f58e..bf2805f8 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.161 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. +* Scaling: rotation merging now also fits a correction surface over **rotation angle × detector position**, which reaches the absorption left behind when the diffracting volume drifts through the beam; like the other surfaces it is cross-validated and applied only where it improves reflections it was not fitted on. * Scaling: the rotation correction surfaces (decay, absorption, detector-plane modulation) no longer pull each cell's factor towards zero in proportion to that cell's own noise, which had put a spurious centre-to-edge ramp on runs whose outer resolution shells carry no signal. * rugnux: A goniometer that turned further than the angles stored in the file — which are the commanded ones — is now measured and **corrected in the second rotation pass**, or set by hand with `--rotation-scale `; the fitted correction is applied only when it exceeds 0.5 %, moves each end of the sweep by at least 0.5°, and comes out the same with any fifth of the sweep left out. * rugnux: Every `mx` and `scale` run now writes `_report.txt`, a results report modelled on XDS's `CORRECT.LP` — `KEY= value` lines, fixed-width tables and `WARNING:` sentences covering indexing, geometry post-refinement, the space-group decision, merging, twinning, radiation damage and the stretches of the sweep over which the crystal delivered much less than the rest of the run. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index ee854bd0..2a12663e 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -747,11 +747,12 @@ The combine groups each reflection's partials into rocking events (contiguous ru The fulls are then re-scaled in the XDS sense — a per-image scale refit directly on the complete reflections under the unity partiality model — and merged (§10.4). Because every merged observation is now a counting-statistics-limited full rather than a partiality-divided slice, the error model reaches a far higher asymptotic $I/\sigma$. -After scale-fulls, three **correction surfaces** are fitted on the combined fulls (rotation path, **on by default**; disable all with `--no-scaling-corrections`), each an alternating multiplicative refinement of the per-full scale against the merged reference: +After scale-fulls, four **correction surfaces** are fitted on the combined fulls (rotation path, **on by default**; disable all with `--no-scaling-corrections`), each an alternating multiplicative refinement of the per-full scale against the merged reference: - **Decay.** Radiation damage weakens later frames more at higher resolution — a resolution×time (Debye–Waller) systematic the resolution-flat per-image scale cannot capture. A single global relative-$B$ rate is fitted, $\ln(I_\mathrm{ref}/I_\mathrm{obs}) = 2\,(\mathrm{d}B/\mathrm{d}n)\,(n-\bar n)\,s^2$ (frame $n$, $s^2 = 1/4d^2$), and folded into the scale. It engages only when the total relative-$B$ over the run exceeds a physical floor (2 Ų); below that the decay is negligible and "correcting" it only spreads symmetry equivalents (same $s^2$, different frames). An optional **per-batch relative-$B$** (`--relative-b[=deg]`, off unless requested; 10°-of-rotation batches by default) extends the single global rate to a smooth $B(n)$ curve — the same $s^2$-weighted decay fit solved independently over short frame batches, curvature-penalized so it cannot over-fit and cross-validated like the surfaces below — for crystals whose decay is non-linear in dose. Its cross-validation splits on **ASU-group parity**, not the frame parity the surfaces below use: a per-batch parameter owns whole frames and so cannot be scored on a held-out frame, whereas splitting the symmetry equivalents tests whether a batch's $B$ generalises to reflections it was not fitted on. - **Absorption.** A smooth multiplicative factor over the diffracted-beam direction expressed in the goniometer (crystal) frame: each full's predicted detector position gives the lab diffracted direction, de-rotated by the spindle so a fixed crystal-frame direction is sampled at many rotation angles and its grid cell is well-determined. Negligible at hard X-rays / thin crystals; it matters at low photon energy. - **Modulation** (detector-plane flat-field). A smooth multiplicative factor over where each reflection lands on the detector (predicted $x,y$): symmetry-equivalents land at different positions as the crystal rotates, over-determining the surface. It absorbs detector-response and geometric systematics that inflate $R_\mathrm{meas}$. +- **Time-dependent absorption.** The same surface as *Absorption*, but over (rotation angle × detector position) instead of the crystal-frame direction alone. The two agree while the illuminated volume stays put — the incident path then depends only on the spindle angle, which the per-image scale already takes, and the exit path is fixed in the crystal frame. Once the diffracting volume drifts through the beam the exit path becomes a function of the spindle angle as well, and nothing time-independent describes it. Fitted last, on 12 rotation bins × a 10×10 equal-occupancy detector grid, so the two time-independent surfaces get first claim on what they can explain. Each surface is **cross-validated**: fitted on even-numbered frames and kept only if it improves the held-out odd-frame agreement by a clear margin (and vice versa), scored by a **σ-independent, $R_\mathrm{meas}$-like** fractional agreement $\sum|I_s-I_\mathrm{ref}|/\sum|I_\mathrm{ref}|$ rather than a studentized $\chi^2$ — so a surface cannot "pass" by reshaping the sigmas instead of tightening the intensities. A surface fitted to noise where its systematic is absent does not generalize and is discarded — a correction never adds scatter. diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 08ce4a96..cd303138 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1510,6 +1510,61 @@ void RotationScaleMerge::RefineAbsorption(int n_iter, int n_groups) { ApplyCellSurface(cell, NB * NB, n_iter, n_groups, "Absorption (goniometer-frame 8x8)"); } +void RotationScaleMerge::RefineAbsorptionTime(int n_iter, int n_groups) { + // Time-dependent absorption: the same multiplicative surface as RefineAbsorption, but indexed by + // (rotation, detector position) rather than by the de-rotated crystal-frame direction alone. For a rigid + // absorber illuminating a fixed volume the two are the same model - the incident path depends on the + // spindle angle only, which the per-frame scale already takes, and the exit path is then fixed in the + // crystal frame. What breaks that factorisation is the diffracting volume moving through the sweep: a + // crystal larger than the beam, a mis-centred loop, ice building up. The exit path then depends on the + // spindle angle AND the direction, and no time-independent surface reaches it. XDS parameterises it the + // same way, spending most of its correction budget there. + // Following Kabsch (2010) Acta Cryst. D66, 133-144 + if (fulls.empty() || n_frames <= 0) + return; + // 12 time bins x a 10 x 10 detector grid = 1200 factors. Held out on unseen frames over 34 rotation + // datasets and scored PER RESOLUTION SHELL - the gate's own whole-range ratio is lowered by any + // resolution-dependent scale without a reflection getting tighter, so it cannot be trusted to choose a + // grid. On that score the median gain only moves between 3.1 % and 4.1 % as the grid goes from 216 to + // 2400 cells, so the grid is second order; 12 x 10 is where the net across the set is largest and where + // fewest crystals are hurt, 3 of the 22 it engages on. The time axis wants to be coarse either way - + // the crystal drifts through the beam slowly. + constexpr int NT = 12; // time bins across the sweep + constexpr int ND = 10; // detector bins per axis + auto usable = [&](const Obs &o) { + return o.group >= 0 && o.corr > 0.0f && std::isfinite(o.corr) && o.partiality >= min_partiality + && std::isfinite(o.px) && std::isfinite(o.py); + }; + // Equal-OCCUPANCY detector bins. An equal-width grid starves the edges and the corners, where a + // half-empty cell is exactly where a free surface over-fits; equal occupancy also spends the grid + // where the reflections are, which on a rotation dataset is the ring the lattice actually fills. + std::vector vx, vy; + vx.reserve(fulls.size()); + vy.reserve(fulls.size()); + for (const Obs &o : fulls) + if (usable(o)) { vx.push_back(o.px); vy.push_back(o.py); } + if (static_cast(vx.size()) < NT * ND * ND) + return; + std::sort(vx.begin(), vx.end()); + std::sort(vy.begin(), vy.end()); + std::vector ex(ND - 1), ey(ND - 1); + for (int i = 1; i < ND; ++i) { + ex[i - 1] = vx[vx.size() * i / ND]; + ey[i - 1] = vy[vy.size() * i / ND]; + } + const int frames_per_bin = std::max(1, n_frames / NT); + std::vector cell(fulls.size(), -1); + for (size_t i = 0; i < fulls.size(); ++i) { + const Obs &o = fulls[i]; + if (!usable(o)) continue; + const int it = std::min(NT - 1, o.frame / frames_per_bin); + const int ix = static_cast(std::upper_bound(ex.begin(), ex.end(), o.px) - ex.begin()); + const int iy = static_cast(std::upper_bound(ey.begin(), ey.end(), o.py) - ey.begin()); + cell[i] = (it * ND + ix) * ND + iy; + } + ApplyCellSurface(cell, NT * ND * ND, n_iter, n_groups, "Absorption (time x detector 12x10x10)"); +} + void RotationScaleMerge::RefineModulation(int n_iter, int n_groups) { // Detector-plane modulation (flat-field): a smooth multiplicative factor over the position where a // reflection lands on the detector (predicted px, py). It absorbs detector-response and geometric @@ -3009,6 +3064,8 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) { RefineAbsorption(absorption_iter, n_groups); if (!for_search && modulation_iter > 0) RefineModulation(modulation_iter, n_groups); + if (!for_search && absorption_iter > 0) + RefineAbsorptionTime(absorption_iter, n_groups); // last: the static surfaces get first claim #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 9e0345b8..e0b1edcb 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -287,6 +287,11 @@ private: // tests whether a batch's B generalises to reflections it was not fit on). Opt-in (--relative-b). void RefineRelativeB(int n_groups); void RefineAbsorption(int n_iter, int n_groups); + // Time-dependent absorption: the same cross-validated surface, indexed by (rotation, detector position) + // instead of by the crystal-frame direction alone. RefineAbsorption's parameterisation is the whole + // model only while the illuminated volume stays put; once the crystal drifts through the beam the exit + // path becomes a function of the spindle angle too, and nothing time-independent reaches it. + void RefineAbsorptionTime(int n_iter, int n_groups); // Detector-plane modulation (flat-field): the same cross-validated surface fit as absorption, but the // cell is the predicted detector position (px, py) instead of the goniometer-frame direction. Corrects // detector-response / geometric systematics that vary with where a reflection lands; because it lives