Scaling: correct the absorption that changes as the crystal turns
RefineAbsorption indexes its surface by the diffracted direction de-rotated into the crystal frame, deliberately without a time axis, and RefineModulation indexes its by detector position, also without one. Nothing is indexed by (rotation, detector position), so the part of the absorption that changes as the crystal turns has no parameter at all. For a rigid absorber illuminating a fixed volume that is the right model: the incident path is a function of the spindle angle alone and the per-image scale takes it, and the exit path is then fixed in the crystal frame. What breaks the factorisation is the diffracting volume moving - a crystal larger than the beam, a mis-centred loop, ice building up. The exit path then depends on the spindle angle as well as the direction, and no time-independent surface reaches it. Measured on 34 rotation datasets, on XDS's own uncorrected intensities, as what is left after the crystal-frame absorption and detector modulation surfaces have taken what they can. The cross-validation gate lets the surface engage on 22 of them. Scored per resolution shell - the gate's whole-range ratio is lowered by any resolution-dependent scale without a reflection getting tighter, so the honest readout is each shell's own ratio, which a per-shell scale leaves unchanged - the median engaged crystal gains 4.1 %, the set gains 117 % summed against 19 % of damage, and 3 of the 22 are hurt. The surface has to be smooth in rotation angle to be absorption at all, and it is: the lag-1 autocorrelation of the fitted factor along the time axis runs +0.32 to +0.71 on the crystals it engages, against -0.08 for the same surface with its time bins shuffled. Where it is not smooth it is fitting something else, and says so - on a sweep whose beam was obstructed for a 70 deg wedge the autocorrelation is +0.16 and the profile is a cliff at the wedge, not a turn. Two null controls. Assign every observation a random cell and the gate refuses it (-1.2 % to -3.8 %). Keep the detector bin and shuffle only the time bin - a surface that cannot contain any time-dependent information - and the gate refuses that too, at +0.04 %, -0.60 % and +0.33 % on three crystals. Against the real surface's +3.7 % to +14.8 % on the same three. 12 time bins x a 10 x 10 detector grid = 1200 factors. On the per-shell score the median gain moves only between 3.1 % and 4.1 % across grids from 216 to 2400 cells, so the grid is second order; 12 x 10 has the largest net and the fewest crystals hurt. Equal-occupancy detector bins, not equal width: an equal-width grid starves the edges and the corners, and a starved cell is where a free surface over-fits. Fitted last, so the two time-independent surfaces get first claim on what they can explain. QUALIFICATION, measured after this was written: the "33 better / 0 worse" above is overall R_meas, which is a ratio of sums across every shell and is therefore lowered by any resolution-dependent scale without a reflection getting tighter - the same property that let the estimator bias pass its own gate. Scored per resolution shell instead, this surface HURTS 6 of 18 crystals under the acceptance gate as it currently stands, because that gate shares the defect and admits the surface where it should not. With a per-shell gate the surface is refused on exactly those crystals and its net over the chain goes from +86.3 to +159.5 per cent with none worse. The correction is right; the gate that decides where to apply it is the next commit's problem, not this one's. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Full 38-crystal rotation battery against its own matched baseline - the same binary with the corrected estimator and without this surface: R_meas better 33 / worse 0, -78.7 R_meas_lo better 24 / worse 2, -20.7 CC1/2 better 8 / worse 0, +9.6 ISa better 30 / worse 2, +99.15 space groups unchanged at 35/38 The low-energy datasets gain most, which is what absorption should do: at 5-6 keV one crystal goes ISa 24.68 -> 37.42 and another 14.21 -> 22.08, while the same protein measured at 13 keV moves 13.41 -> 14.83. Taken with the estimator fix it precedes, against a clean baseline: R_meas_lo better 25 / worse 3 summed -26.0, ISa better 31 / worse 2 summed +113.1, outer-shell CC1/2 +83.0, observations +25 880 on 36 crystals of 38, CC1/2 flat at -1.4 and no space group moved. That last number is the point of the pair: the estimator fix alone costs CC1/2 -11.4, because the ramp it removes was partly standing in for this correction. One cost, predicted in advance and still unexplained: outer-shell CC1/2 falls on three of the four low-energy crystals, by 15.8 points on the worst, while every other statistic on those same crystals improves. The fourth goes up. On 5000-9000 That outer-shell fall has since been attributed, and it is not this surface: with the merge's 6-sigma outlier rejection turned off, the sign flips on every crystal that lost, +20.5 and +20.7 where it read -15.8 and -8.0. The surface removes most of the deviants in sample - it is fitted on all the data and applied to it, with no robustness of its own - so the merge's cut stops firing and the survivors land in a shell whose multiplicity is about three. Last-shell CC1/2 is largely made by that cut: one crystal's baseline goes 10.4 to 92.7 purely by dropping 19 per cent of the shell. A second qualification, measured after the numbers above were taken: overall R_meas is a ratio of sums across every shell, so any resolution-dependent scale lowers it without a reflection getting tighter - the same property that let the estimator bias pass its own gate. Scored per resolution shell instead, this surface hurts 6 of 18 crystals under the acceptance gate as it stands, because that gate shares the defect and admits the surface where it should not. Under a per-shell gate it is refused on exactly those crystals and its net over the correction chain goes from +86.3 to +159.5 per cent with none worse. The correction is right; where to apply it is the gate's problem. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5
@@ -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).
|
||||
|
||||
@@ -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 <k>`; 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 `<prefix>_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.
|
||||
|
||||
@@ -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.
|
||||
|
||||
|
||||
@@ -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<float> 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<int>(vx.size()) < NT * ND * ND)
|
||||
return;
|
||||
std::sort(vx.begin(), vx.end());
|
||||
std::sort(vy.begin(), vy.end());
|
||||
std::vector<float> 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<int32_t> 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<int>(std::upper_bound(ex.begin(), ex.end(), o.px) - ex.begin());
|
||||
const int iy = static_cast<int>(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.
|
||||
|
||||
@@ -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
|
||||
|
||||
Reference in New Issue
Block a user