From 80e10f51b31a48d37d71f0710d29f1624aa0adea Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 25 Aug 2026 19:13:16 +0200 Subject: [PATCH] rugnux: measure and report diffraction anisotropy rugnux now says whether a dataset's fall-off is direction-dependent, and by how much. It corrects nothing and truncates nothing: no intensity is changed, no reflection is dropped on a directional criterion, and the written files do not depend on direction at all. Two quantities, because they are not the same thing. The anisotropic deltaB is the range of the principal components of the anisotropy tensor - a rate of fall-off. The diffraction limit along each principal direction is where in a 20 degree cone falls through 2 - where signal actually runs out. One battery case has only 0.28 A between its directional limits and a 58x ratio in cone , so reporting either alone would miss it. The tensor is a Laue-constrained deviatoric ADP tensor fitted on INTENSITIES with no positivity cut, by weighted Gauss-Newton over 12 shells x 60 directions with a free constant per shell. Fitting amplitudes after a positivity cut, which is what xtriage and ctruncate do, destroys about 40% of the measured anisotropy - the cut keeps only the positive noise excursions in whichever direction has died, and that is the direction carrying the signal. Against the same 38 merged files rugnux reads 1.24x xtriage's eigenvalue spread and 1.61x ctruncate's; on strong near-isotropic data all three agree to a few percent, and they diverge exactly where a direction has died. The verdict is gated three ways - not detected, detected, or cannot determine - against the dataset's own systematic floor, measured in the tensor directions its Laue symmetry forbids. The floor cannot be measured on merged reflections, which have exact Laue symmetry by construction, so the floor is taken from the unmerged observations and the verdict is "cannot determine" without them. Triclinic has no forbidden subspace and always returns cannot determine. A cubic crystal returns exactly zero, because that is its symmetry and not a measurement. A second axis reports the resolution signature: a genuine Debye-Waller fall-off is linear through the origin in s^2, and a deficit that is flat is something else. Magnitude alone had promoted a crystal that is 68% not a Debye-Waller B into the top five of this battery; it now reads not detected with the caution attached. Following Sheriff & Hendrickson (1987) Acta Cryst. A43, 118-121 for the tensor and Popov & Bourenkov (2003) Acta Cryst. D59, 1145-1153 for the estimator. The directional limits are written as jfjoch_ local mmCIF items rather than _reflns.pdbx_aniso_diffraction_limit_*, whose dictionary definition is explicitly the ellipsoid fitted to a diffraction cut-off surface - a construction rugnux does not perform. The generic anisotropic B tensor items are written. Changes no existing number; only REPORT_VERSION moves, 1 to 2. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P --- docs/ACKNOWLEDGEMENT.md | 15 + docs/CHANGELOG.md | 1 + docs/CPU_DATA_ANALYSIS.md | 28 +- image_analysis/WriteReflections.cpp | 42 + .../scale_merge/AnisotropyAnalysis.cpp | 1357 +++++++++++++++++ .../scale_merge/AnisotropyAnalysis.h | 154 ++ image_analysis/scale_merge/CMakeLists.txt | 2 + image_analysis/scale_merge/Merge.h | 7 + rugnux/ResultReport.cpp | 48 +- rugnux/Rugnux.cpp | 22 + rugnux/rugnux_cli.cpp | 16 + tests/AnisotropyAnalysisTest.cpp | 153 ++ tests/CMakeLists.txt | 1 + 13 files changed, 1842 insertions(+), 4 deletions(-) create mode 100644 image_analysis/scale_merge/AnisotropyAnalysis.cpp create mode 100644 image_analysis/scale_merge/AnisotropyAnalysis.h create mode 100644 tests/AnisotropyAnalysisTest.cpp diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 4afbc131..4211a4e9 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -101,6 +101,21 @@ enumerated from a cell. D. W. Moreau, H. Atakisi and R. E. Thorne, "Ice in biomo cryocrystallography" (2021), Acta Cryst. D77, 540-554 [doi:10.1107/S2059798321001170](https://doi.org/10.1107/S2059798321001170). +**Diffraction anisotropy** — the description of the overall fall-off by a single anisotropic +displacement tensor, its symmetry constraints, and the fact that only its deviatoric part is +determined (the isotropic part being degenerate with the overall scale) are Sheriff and Hendrickson's. +The estimator fits that tensor to the observed intensity distribution, taking sigma(I) into account, +in the sense of Popov and Bourenkov. The directional diffraction limits - in a cone about +each principal direction, and the reporting of the anisotropic deltaB as the range of the principal +components - follow AIMLESS. rugnux reports these; it corrects no intensity and removes no reflection +on a directional criterion. S. Sheriff and W. A. Hendrickson, "Description of overall anisotropy in +diffraction from macromolecular crystals" (1987), Acta Cryst. A43, 118-121 +[doi:10.1107/S010876738709977X](https://doi.org/10.1107/S010876738709977X); A. N. Popov and +G. P. Bourenkov, "Choice of data-collection parameters based on statistic modelling" (2003), Acta +Cryst. D59, 1145-1153 [doi:10.1107/S0907444903008163](https://doi.org/10.1107/S0907444903008163); +P. R. Evans and G. N. Murshudov, "How good are my data and what is the resolution?" (2013), Acta +Cryst. D69, 1204-1214 [doi:10.1107/S0907444913000061](https://doi.org/10.1107/S0907444913000061). + **Data-quality statistics** follow the established conventions rather than any one program: R_meas and R_pim, CC1/2 and CC\*, and the reporting of I/sigma(I). K. Diederichs and P. A. Karplus, "Improved R-factors for diffraction data analysis in macromolecular crystallography" (1997), Nat. Struct. Biol. diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 99720192..db161376 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -5,6 +5,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * The PCIe driver DKMS package builds for the kernel it is being installed for instead of the running one, so a module built while a kernel update is being applied loads after the reboot. * The PCIe driver builds on RHEL 9.5 and later, and on their CentOS Stream, Rocky and AlmaLinux equivalents, where the `vm_flags` kernel interface was backported into the 5.14 kernel. +* rugnux reports diffraction anisotropy: the anisotropic deltaB and the diffraction limit along each * rugnux: `--mode scale` reports the Wilson B-factor estimate instead of `WILSON_B= nan`. * rugnux: `--export-unmerged` writes the integrated observations as `_unmerged.mtz`, an unmerged MTZ readable by aimless, pointless, careless and `iotbx.merging_statistics`, in `--mode mx` and `--mode scale` alike; each rotation reflection's partials are summed into one full, and `--export-unmerged-partials` writes one row per image instead. Lattice-centring absences are not written; screw and glide absences are. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 675d74b0..22d22f4a 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -37,6 +37,8 @@ The methods draw on, and in places reimplement, solutions from: - P. Evans, "Scaling and assessment of data quality", *Acta Cryst.* **D62** (2006), 72-82, and P. R. Evans, *Acta Cryst.* **D67** (2011), 282-292 (POINTLESS: operator-by-operator point-group scoring, and the axial-zone screw-absence test). - A. G. W. Leslie & H. R. Powell, "Processing diffraction data with MOSFLM" (2007), NATO Science Series II **245**, 41-51 (post-refinement practice: what is refined per image and what over a wedge). - D. W. Moreau, H. Atakisi & R. E. Thorne, "Ice in biomolecular cryocrystallography", *Acta Cryst.* **D77** (2021), 540-554 (measured hexagonal-ice ring positions, used by the ice-ring score, the ice flagging and the ice calibrant). +- S. Sheriff & W. A. Hendrickson, "Description of overall anisotropy in diffraction from macromolecular crystals", *Acta Cryst.* **A43** (1987), 118-121, and A. N. Popov & G. P. Bourenkov, *Acta Cryst.* **D59** (2003), 1145-1153 (the overall anisotropic B tensor, its symmetry constraints, and its estimation from the observed intensities). +- 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). - 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). @@ -958,7 +960,31 @@ Merging applies an optional per-observation median-based $N\sigma$ cut (`--rejec By default the reported/written high-resolution limit is trimmed where $\mathrm{CC}_{1/2}$ falls off: a logistic is fitted to $\mathrm{CC}_{1/2}(s)$, and the limit is set **one reported-shell width past** the point where the fit crosses 0.30 — deliberately "one shell too far", so weak-but-real data below the crossing are kept rather than discarded. The extension is measured over the range that is actually kept, not the full measured range, so a detector reaching far past where the crystal diffracts cannot inflate it. `--scaling-high-resolution` overrides the limit and `--resolution-cutoff off` disables it. -### 13.5 Practical notes and limitations +### 13.5 Diffraction anisotropy + +How fast the intensity falls off with resolution can depend on direction. rugnux measures that, reports it, and does nothing else with it: no intensity is corrected, no reflection is removed on a directional criterion, and the merged data and the written files do not depend on direction at all. + +**The tensor.** A deviatoric anisotropic displacement tensor is fitted to the merged intensities as + +$$\ln \langle I(\mathbf{s})\rangle = c(\text{shell}) - \tfrac{1}{2}\,\mathbf{s}^\mathsf{T} B\, \mathbf{s},\qquad \mathbf{s} = \text{reciprocal-space vector},\ |\mathbf{s}| = 1/d$$ + +with one free constant per resolution shell, so every isotropic feature — the Wilson curve, an ice ring, a noise floor, a scaling error — is absorbed exactly and only the $\ell = 2$ angular part drives the tensor. For an isotropic $B$ this reduces to the ordinary Wilson plot, so $B$ here is the ordinary crystallographic ($B = 8\pi^2 U$) $B$, directly comparable with phenix.xtriage's `B_cart`, ctruncate's anisotropic $B$ eigenvalues and AIMLESS's anisotropic $\Delta B$. Only the deviatoric part is fitted: the isotropic part is degenerate with the overall scale. The tensor is constrained to the directions the Laue class allows — five free deviatoric parameters in triclinic, three in monoclinic, two in orthorhombic, one in tetragonal, trigonal and hexagonal, and **none at all in cubic**, where symmetry forces $\Delta B$ to be exactly zero. + +The fit is on **intensities, with no positivity cut**. Fitting amplitudes, or dropping non-positive intensities as an amplitude-based tool must, loses roughly 40% of the signal: in a direction that has died half the merged intensities are negative, so a positivity cut keeps only the positive noise excursions and flattens the fall-off exactly where the anisotropy is largest. + +**Two different quantities are reported, and they are not interchangeable.** $\Delta B$ (the range of the principal components) is a *rate*; the diffraction limit along each principal direction — where $\langle I/\sigma(I)\rangle$ in a 20° cone about that direction falls through 2 — is where the signal actually runs out. A crystal can have a large $\Delta B$ and almost no spread in directional limit, or the reverse. + +**The resolution signature.** A genuine Debye–Waller $B$ makes the directional deficit a straight line through the origin in $s^2$. The per-shell $\ell = 2$ amplitude is therefore fitted against $s^2$ and the curve is classified: *linear* (a real $B$), *flat* (a deficit that does not follow $\exp(-\tfrac12 \mathbf{s}^\mathsf{T} B \mathbf{s})$ at all, so the fitted $\Delta B$ describes the data with the wrong functional form and may be an **under**-estimate), or *convex* (a deficit that grows faster than $s^2$, which a $B$ cannot do). The verdict is re-derived at 8 and at 16 shells, and reported as undetermined if it moves. + +**The verdict, and what it is measured against.** Whether an anisotropy is real is not decided against a counting-statistics error bar. Real data carry systematic error far larger than counting error, and gating on the latter reports anisotropy on datasets that have none. Instead the data set measures its own systematic error: in the tensor directions the Laue class *forbids*, the true tensor is exactly zero whatever the crystal is, so whatever is measured there is systematic. That measurement needs the unmerged observations — a merge has exact Laue symmetry by construction, and the forbidden directions are identically zero in it — so it is made on the scaled, rocking-curve-assembled observations. The counting part is subtracted, the counting error of the directions actually being tested is added back, and the ratio of $\Delta B$ to the result is banded: below 2 not established, 2–3.5 marginal, above 3.5 established, above 5 strong. + +The report says **NOT DETECTED**, **DETECTED**, or **CANNOT DETERMINE**, and the third is a real answer rather than an evasion. It is returned when the Laue class is triclinic (no forbidden direction exists, so there is no internal measurement of the systematic error and no substitute for it), when the observed rotation range is under about 90° (a lab-fixed systematic then reaches several tensor directions instead of one), when the merged data are at the noise floor, when the scale model carried no dose term (an uncorrected dose ramp manufactures anisotropy that no significance test can see through), or when no unmerged observations were available. The smallest $\Delta B$ that could have been established on the data set is reported with the verdict; it is set by the systematic error rather than by counting, so it does **not** improve with more reflections or a longer exposure. + +A too-high space-group assignment is the one failure mode that is silent: real anisotropy is then pushed into the directions used to measure the systematic error, which inflates the floor and biases the answer towards reporting none. A caution says so wherever the Laue class leaves a single free direction. + +Everything lands in `_report.txt` section 9 (`ANISOTROPY_*` keys), in the printed statistics, and in the merged mmCIF: the eigen-decomposition of the tensor as the standard `_reflns.pdbx_aniso_B_tensor_*` items (relative to the weakest direction, since only the deviatoric part is determined), and the directional limits, the shape and the verdict under the `_reflns.jfjoch_aniso_*` local prefix. + +### 13.6 Practical notes and limitations - **Bragg integration is profile-fitted by default** (per-shell Gaussian profile, Kabsch extraction; §9.3), with plain box summation available as a fallback (`--integrator boxsum`). The profiles are built per frame from that frame's strong spots, which suits fast-feedback and serial/streaming use; a profile shared across many frames (as in full offline workflows) is not currently formed. - **Space-group symmetry** beyond centering absences is not enforced during prediction/integration unless the space group is supplied and used downstream. diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 2ebfc046..c4e23584 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -274,6 +274,48 @@ void WriteMmcifReflections(const std::vector &reflections, if (std::isfinite(statistics.radiation_damage_delta_b)) out << "_reflns.jfjoch_radiation_damage_relative_B " << Fmt(statistics.radiation_damage_delta_b, 2) << " # relative-B first->last over the run (A^2); + = high-res fades with dose\n"; + // Diffraction anisotropy. The eigen-decomposition of the anisotropy tensor has standard PDBx + // items; the eigenvalues there must be non-negative, and only the deviatoric part of the tensor + // is determined at all (its isotropic part is degenerate with the overall scale), so they are + // written relative to the weakest direction - eigenvalue_3 is 0 by construction and + // eigenvalue_1 is the anisotropic deltaB. The eigenvectors are in the PDB orthogonalisation + // convention, which is the one gemmi (and hence rugnux) uses throughout. + // + // The directional diffraction LIMITS are deliberately NOT written as + // _reflns.pdbx_aniso_diffraction_limit_*: the dictionary defines those as the semi-axes of an + // ellipsoid fitted to a diffraction cut-off surface, which is a different construction from the + // one below and one rugnux does not perform - it cuts nothing on a directional criterion. They + // go under the jfjoch local prefix with their own definition instead. + const auto &an = statistics.anisotropy; + if (an.n_reflections > 0 && an.n_cells > 0 && std::isfinite(an.delta_b)) { + out << "_reflns.pdbx_orthogonalization_convention pdb\n"; + for (int i = 0; i < 3; ++i) { + out << "_reflns.pdbx_aniso_B_tensor_eigenvalue_" << (i + 1) << " " + << Fmt(an.eigenvalue[i] - an.eigenvalue[2], 2) << "\n"; + for (int j = 0; j < 3; ++j) + out << "_reflns.pdbx_aniso_B_tensor_eigenvector_" << (i + 1) << "_ortho[" << (j + 1) + << "] " << Fmt(an.eigenvector[i][j], 4) << "\n"; + } + out << "_reflns.jfjoch_aniso_delta_B " << Fmt(an.delta_b, 2) + << " # range of the principal components (A^2), fitted on intensities\n"; + if (std::isfinite(an.delta_b_linear)) + out << "_reflns.jfjoch_aniso_delta_B_linear " << Fmt(an.delta_b_linear, 2) + << " # deltaB implied by the s^2 slope alone: what a Debye-Waller B accounts for\n"; + out << "_reflns.jfjoch_aniso_shape " << AnisotropyShapeCode(an.shape) + << " # resolution signature of the directional deficit\n"; + if (std::isfinite(an.floor)) + out << "_reflns.jfjoch_aniso_floor " << Fmt(an.floor, 3) + << " # deltaB this data set's own systematic error could manufacture (A^2)\n"; + if (std::isfinite(an.significance)) + out << "_reflns.jfjoch_aniso_significance " << Fmt(an.significance, 2) + << " # deltaB(linear) / floor\n"; + out << "_reflns.jfjoch_aniso_verdict " << AnisotropyVerdictCode(an.verdict) << "\n"; + for (int i = 0; i < 3; ++i) + if (std::isfinite(an.d_min_axis[i])) + out << "_reflns.jfjoch_aniso_d_min_" << (i + 1) << " " + << Fmt(an.d_min_axis[i], 2) + << " # = 2 in a 20 deg cone about eigenvector " << (i + 1) << "\n"; + } out << "#\n"; // Per-batch relative-B curve (the radiation-damage monitor, rotation): one relative Debye-Waller B diff --git a/image_analysis/scale_merge/AnisotropyAnalysis.cpp b/image_analysis/scale_merge/AnisotropyAnalysis.cpp new file mode 100644 index 00000000..84f50b76 --- /dev/null +++ b/image_analysis/scale_merge/AnisotropyAnalysis.cpp @@ -0,0 +1,1357 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "AnisotropyAnalysis.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "HKLKey.h" + +#include "gemmi/eig3.hpp" +#include "gemmi/scaling.hpp" + +namespace { + + using Tensor = gemmi::SMat33; + + // ---------------------------------------------------------------- tuned constants + // The binning is the one the shape test was validated on: 12 equal-count resolution shells x 60 + // directions on a hemisphere, a cell entering the fit once it holds 12 independent reflections. + // Shell count, direction count and the cell-population floor were each varied by a factor of two + // either way without moving a verdict, so none of the three is sharp. + constexpr int N_SHELLS = 12; + constexpr int N_DIRECTIONS = 60; + constexpr int MIN_CELL_REFLECTIONS = 12; + constexpr int MIN_CELLS = 40; // fewer than this and no tensor is fitted at all + constexpr int MIN_SHAPE_SHELLS = 5; // fewer than this and the s^2 signature is not fitted + constexpr int MIN_CELLS_PER_SHAPE_SHELL = 8; + + // The whole diagnostic is refused below this: at ~ 1 neither the tensor nor its shape is + // interpretable, and one battery dataset at 0.65 produced numbers that mean nothing. + constexpr double MIN_MEAN_I_OVER_SIGMA = 1.0; + + // Below ~90 deg of OBSERVED rotation a lab-fixed systematic reaches a second and a third tensor + // direction (canonical correlation 0.74-0.85 against 0.05 at a full sweep), the anisotropy such a + // systematic manufactures grows 5-9x, and the power to detect a real 10 A^2 falls to <= 0.25. + constexpr double MIN_ROTATION_DEG = 90.0; + + // Symmetry-expanded reflections are capped by striding the unique list: the estimator is a fit of a + // few tens of cell means and does not improve past this, while the memory does grow. + constexpr size_t MAX_ENTRIES = 4000000; + + // The unmerged arm strides whole unique reflections down to this many. The systematic-error floor is + // set by the systematic, not by counting: over a 12x range in the number of unique reflections and a + // 10x range in the smallest anisotropy the gate can establish moves by less than 2x, so + // more than this buys nothing and costs time and memory. + constexpr int MAX_CLUSTERS = 20000; + + // The cone half-angle of the directional diffraction limits, and the the limit is read + // at. Both are AIMLESS's shipped defaults (ANALYSIS CONE 20, ISIGMINIMUM 1.5 -> here 2.0, the value + // this project already uses for a resolution statement). + constexpr double CONE_HALF_ANGLE_DEG = 20.0; + constexpr double CONE_I_OVER_SIGMA = 2.0; + + // The s^2 signature. FLAT is declared at c0 > +0.70 with c0/sigma(c0) > 6: over a 38-crystal battery + // that fires on six of the seven crystals whose deficit does not follow a Debye-Waller B, with no + // false positive among the other 24, and the false-positive count stays 0 anywhere between z > 5 and + // z > 8. A pure-B simulation on the same reflection lists returns |c0| <= 0.18, so the machinery + // cannot invent an intercept of this size. The curvature arm uses the same z for symmetry; it has no + // false-positive calibration of its own, and it is reported rather than acted on. + constexpr double FLAT_INTERCEPT = 0.70; + constexpr double SHAPE_Z = 6.0; + // CONVEX asks for the s^4 term to carry more than half the deficit at the resolution limit, not + // merely to be significant: a mild positive curvature is common and is not a different mechanism. + constexpr double CONVEX_SHARE = 0.5; + + // Gate bands. The denominator is this dataset's own systematic error scale, not a counting-statistics + // error bar: against known ground truth a counting-statistics gate calls "real" in 15-58% of clean + // isotropic datasets, because the symmetry-forbidden directions are determined from contrasts BETWEEN + // symmetry mates of one reflection (which share their true |F|^2, so Wilson scatter cancels) while + // the symmetry-allowed ones can only be determined BETWEEN different reflections (which do not). + // 3.5 is a 2-5% false-positive rate pooled over six artefact classes and 6-14% in the worst of them; + // 5.0 is 0-0.5% and 0-3%. + constexpr double GATE_MARGINAL = 2.0; + constexpr double GATE_ESTABLISHED = 3.5; + constexpr double GATE_STRONG = 5.0; + + // E[dB] / sigma_per_component for a zero-mean Gaussian tensor confined to a subspace: the eigenvalues + // of a noisily estimated tensor repel, so dB has a strictly positive expectation even when the tensor + // is exactly zero. Per-subspace Monte-Carlo medians, indexed by the dimension of the subspace; + // K_ALLOWED is for the symmetry-allowed subspaces, K_FORBIDDEN for the symmetry-forbidden ones (the + // two differ at the same dimension because they are different subspaces of the deviatoric space). + constexpr double K_ALLOWED[6] = {0.0, 0.824, 1.583, 2.075, 0.0, 2.872}; + constexpr double K_FORBIDDEN[6] = {0.0, 0.0, 1.66, 2.13, 2.53, 2.87}; + + // ---------------------------------------------------------------- tensor algebra + double FrobDot(const Tensor &a, const Tensor &b) { + return a.u11 * b.u11 + a.u22 * b.u22 + a.u33 * b.u33 + + 2.0 * (a.u12 * b.u12 + a.u13 * b.u13 + a.u23 * b.u23); + } + + Tensor Scale(const Tensor &a, double s) { + return {a.u11 * s, a.u22 * s, a.u33 * s, a.u12 * s, a.u13 * s, a.u23 * s}; + } + + Tensor Deviatoric(const Tensor &a) { return a.added_kI(-a.trace() / 3.0); } + + // s^T E s from the cell-averaged outer product q = <(sx^2, sy^2, sz^2, sx sy, sx sz, sy sz)>. + double QuadForm(const Tensor &e, const std::array &q) { + return e.u11 * q[0] + e.u22 * q[1] + e.u33 * q[2] + + 2.0 * (e.u12 * q[3] + e.u13 * q[4] + e.u23 * q[5]); + } + + // Gram-Schmidt in the Frobenius inner product; directions that the earlier ones already span drop out. + std::vector Orthonormalise(const std::vector &raw) { + std::vector out; + for (Tensor t : raw) { + for (const auto &o : out) + t = t - Scale(o, FrobDot(t, o)); + const double n = std::sqrt(FrobDot(t, t)); + if (n > 1e-9) + out.push_back(Scale(t, 1.0 / n)); + } + return out; + } + + // The cell with its lengths and angles snapped to what the crystal system requires. + gemmi::UnitCell MetricIdealCell(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg) { + if (!sg) + return cell; + double a = cell.a, b = cell.b, c = cell.c; + double alpha = cell.alpha, beta = cell.beta, gamma = cell.gamma; + switch (sg->crystal_system()) { + case gemmi::CrystalSystem::Triclinic: + break; + case gemmi::CrystalSystem::Monoclinic: + if (sg->monoclinic_unique_axis() == 'a') { beta = 90.0; gamma = 90.0; } + else if (sg->monoclinic_unique_axis() == 'c') { alpha = 90.0; beta = 90.0; } + else { alpha = 90.0; gamma = 90.0; } + break; + case gemmi::CrystalSystem::Orthorhombic: + alpha = beta = gamma = 90.0; + break; + case gemmi::CrystalSystem::Tetragonal: + a = b = 0.5 * (a + b); + alpha = beta = gamma = 90.0; + break; + case gemmi::CrystalSystem::Trigonal: + case gemmi::CrystalSystem::Hexagonal: + if (sg->ext == 'R') { + a = b = c = (a + b + c) / 3.0; + alpha = beta = gamma = (alpha + beta + gamma) / 3.0; + } else { + a = b = 0.5 * (a + b); + alpha = beta = 90.0; + gamma = 120.0; + } + break; + case gemmi::CrystalSystem::Cubic: + a = b = c = (a + b + c) / 3.0; + alpha = beta = gamma = 90.0; + break; + } + gemmi::UnitCell out; + out.set(a, b, c, alpha, beta, gamma); + return out; + } + + // The symmetry-allowed anisotropy directions, deviatoric, in the Cartesian frame. + // + // gemmi::adp_symmetry_constraints gives the standard ADP symmetry constraints in the crystal-axis + // (B*) parameterisation - 6 / 4 / 3 / 2 / 2 / 1 vectors for triclinic / monoclinic / orthorhombic / + // tetragonal / trigonal-hexagonal / cubic. Transforming each to Cartesian is a bijection on symmetric + // matrices, and the isotropic direction is always allowed, so removing the trace leaves exactly + // 5 / 3 / 2 / 1 / 1 / 0 free deviatoric parameters. + // + // Following Sheriff & Hendrickson (1987) Acta Cryst. A43, 118-121 + std::vector AllowedDeviatoricBasis(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg) { + std::vector raw; + // The constraints are a statement about the IDEAL metric, and rugnux writes the unconstrained + // refined cell: a beta of 89.97 deg in a C2 crystal leaves the isotropic tensor outside the span + // of the monoclinic constraints, the trace is then not fully removed, and the count comes out + // one too high. A cell whose angle is a hundredth of a degree off does not create a fourth + // anisotropy direction, so the basis is built on the idealised metric. + const gemmi::UnitCell ideal = MetricIdealCell(cell, sg); + for (const gemmi::Vec6 &v : gemmi::adp_symmetry_constraints(sg)) { + const Tensor b_star{v[0], v[1], v[2], v[3], v[4], v[5]}; + raw.push_back(Deviatoric(b_star.transformed_by(ideal.orth.mat))); + } + return Orthonormalise(raw); + } + + std::vector FullDeviatoricBasis() { + const std::vector raw{{1, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0}, + {0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}}; + std::vector dev; + dev.reserve(raw.size()); + for (const auto &t : raw) + dev.push_back(Deviatoric(t)); + return Orthonormalise(dev); + } + + // What is left of the deviatoric space once the symmetry-allowed directions are taken out. In those + // directions the true tensor is exactly zero by symmetry, whatever the crystal is - so whatever is + // measured there is this dataset's own systematic error, measured on this dataset. + std::vector ForbiddenBasis(const std::vector &full, const std::vector &allowed) { + std::vector raw; + for (Tensor t : full) { + for (const auto &a : allowed) + t = t - Scale(a, FrobDot(t, a)); + raw.push_back(t); + } + return Orthonormalise(raw); + } + + // ---------------------------------------------------------------- binning + // A Fibonacci spiral on the hemisphere: near-uniform, and no pole or seam where a lattice direction + // could pile up. + std::vector HemisphereDirections(int n) { + std::vector d; + d.reserve(n); + const double golden = gemmi::pi() * (1.0 + std::sqrt(5.0)); + for (int i = 0; i < n; ++i) { + const double z = (i + 0.5) / n; + const double r = std::sqrt(std::max(0.0, 1.0 - z * z)); + const double phi = golden * (i + 0.5); + d.emplace_back(r * std::cos(phi), r * std::sin(phi), z); + } + return d; + } + + struct Cells { + std::vector mu, se, s2; + std::vector> q; + std::vector shell, n; + // The cluster structure, kept so the covariance can be made robust to it: one entry per + // measurement that went into a cell, grouped by the unique reflection it belongs to. + std::vector entry_cell, entry_cluster; + std::vector entry_value; + int n_shells = 0; + }; + + struct Usable { + std::array hkl; + double value; // epsilon-corrected merged intensity, negatives kept + double s2; + int shell = 0; + }; + + std::vector SelectReflections(const std::vector &merged, + const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg, + int n_shells) { + const gemmi::GroupOps gops = sg ? sg->operations() : gemmi::GroupOps{}; + std::vector u; + u.reserve(merged.size()); + for (const auto &r : merged) { + if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || !(r.d > 0.0f)) + continue; + const std::array hkl{r.h, r.k, r.l}; + const double eps = sg ? std::max(1, gops.epsilon_factor(hkl)) : 1; + u.push_back({hkl, r.I / eps, 1.0 / (static_cast(r.d) * r.d), 0}); + } + if (u.empty()) + return u; + // Equal-count shells in s^2. Every symmetry copy of a reflection has the same |s|, so equal-count + // over the expanded set and over the unique set are the same edges. + std::vector order(u.size()); + std::iota(order.begin(), order.end(), 0); + std::sort(order.begin(), order.end(), [&](int a, int b) { return u[a].s2 < u[b].s2; }); + for (size_t i = 0; i < order.size(); ++i) + u[order[i]].shell = std::min(n_shells - 1, static_cast(i * n_shells / order.size())); + return u; + } + + // Keep the cells that hold enough measurements, renumber them, and carry the entry list over. + void FinishCells(Cells &c, const std::vector &sum, const std::vector &sumsq, + const std::vector &s2sum, const std::vector> &qsum, + const std::vector &count, const std::vector &raw_cell, + const std::vector &raw_cluster, const std::vector &raw_value) { + std::vector remap(count.size(), -1); + for (size_t cid = 0; cid < count.size(); ++cid) { + if (count[cid] < MIN_CELL_REFLECTIONS) + continue; + const double n = count[cid]; + const double mean = sum[cid] / n; + const double var = std::max(0.0, (sumsq[cid] - n * mean * mean) / (n - 1.0)); + remap[cid] = static_cast(c.mu.size()); + c.mu.push_back(mean); + c.se.push_back(std::sqrt(var / n)); + c.s2.push_back(s2sum[cid] / n); + c.q.push_back({qsum[cid][0] / n, qsum[cid][1] / n, qsum[cid][2] / n, + qsum[cid][3] / n, qsum[cid][4] / n, qsum[cid][5] / n}); + c.shell.push_back(static_cast(cid) / N_DIRECTIONS); + c.n.push_back(count[cid]); + } + // A cell whose measurements happen to agree exactly has no error bar of its own; give it the + // typical one rather than an infinite weight. + std::vector positive; + for (double s : c.se) + if (s > 0.0) + positive.push_back(s); + if (!positive.empty()) { + std::nth_element(positive.begin(), positive.begin() + positive.size() / 2, positive.end()); + const double median = positive[positive.size() / 2]; + for (double &s : c.se) + if (!(s > 0.0)) + s = median; + } + for (size_t e = 0; e < raw_cell.size(); ++e) + if (remap[raw_cell[e]] >= 0) { + c.entry_cell.push_back(remap[raw_cell[e]]); + c.entry_cluster.push_back(raw_cluster[e]); + c.entry_value.push_back(raw_value[e]); + } + } + + Cells BuildCells(const std::vector &refl, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup *sg, int n_shells) { + Cells c; + c.n_shells = n_shells; + const std::vector ops = sg ? sg->operations().sym_ops : std::vector{gemmi::Op::identity()}; + const auto dirs = HemisphereDirections(N_DIRECTIONS); + // Friedel needs no separate operation: -s folds onto the same hemisphere direction as s. + const size_t stride = std::max(1, (refl.size() * ops.size() + MAX_ENTRIES - 1) / MAX_ENTRIES); + + const int n_cells_total = n_shells * N_DIRECTIONS; + std::vector sum(n_cells_total, 0.0), sumsq(n_cells_total, 0.0), s2sum(n_cells_total, 0.0); + std::vector> qsum(n_cells_total, {0, 0, 0, 0, 0, 0}); + std::vector count(n_cells_total, 0); + std::vector raw_cell, raw_cluster; + std::vector raw_value; + + std::vector> hits; // (cell, s) of this reflection's copies + for (size_t i = 0; i < refl.size(); i += stride) { + const auto &r = refl[i]; + hits.clear(); + for (const auto &op : ops) { + const gemmi::Op::Miller m = op.apply_to_hkl(r.hkl); + gemmi::Vec3 s = cell.frac.mat.left_multiply(gemmi::Vec3(m[0], m[1], m[2])); + if (s.z < 0) + s = gemmi::Vec3(-s.x, -s.y, -s.z); + const double len = s.length(); + if (!(len > 0)) + continue; + const gemmi::Vec3 unit(s.x / len, s.y / len, s.z / len); + int best = 0; + double best_dot = -1.0; + for (int j = 0; j < N_DIRECTIONS; ++j) { + const double dot = std::fabs(unit.dot(dirs[j])); + if (dot > best_dot) { best_dot = dot; best = j; } + } + hits.emplace_back(r.shell * N_DIRECTIONS + best, s); + } + // A reflection contributes to a cell ONCE however many of its symmetry copies land there, so + // a cell's standard error counts independent measurements. + std::sort(hits.begin(), hits.end(), [](const auto &a, const auto &b) { return a.first < b.first; }); + const int32_t cluster = static_cast(i); + int previous = -1; + for (const auto &[cid, s] : hits) { + if (cid == previous) + continue; + previous = cid; + sum[cid] += r.value; + sumsq[cid] += r.value * r.value; + s2sum[cid] += r.s2; + qsum[cid][0] += s.x * s.x; qsum[cid][1] += s.y * s.y; qsum[cid][2] += s.z * s.z; + qsum[cid][3] += s.x * s.y; qsum[cid][4] += s.x * s.z; qsum[cid][5] += s.y * s.z; + count[cid] += 1; + raw_cell.push_back(static_cast(cid)); + raw_cluster.push_back(cluster); + raw_value.push_back(static_cast(r.value)); + } + } + + FinishCells(c, sum, sumsq, s2sum, qsum, count, raw_cell, raw_cluster, raw_value); + return c; + } + + // ---------------------------------------------------------------- the tensor fit + struct TensorFit { + bool ok = false; + Tensor b{0, 0, 0, 0, 0, 0}; + Eigen::MatrixXd cov; // nb x nb, cluster-robust + double chi2red = NAN; + }; + + // mu_cell = K_shell * exp(-1/2 s^T B s), fitted on the INTENSITY scale by weighted Gauss-Newton with + // weights 1/SE^2, one free K per resolution shell profiled alongside the tensor. Because K is free per + // shell, every isotropic feature - the Wilson curve, an ice ring, a noise floor, a scaling error - is + // absorbed exactly, and only the l=2 angular part drives the tensor. Nothing is censored: a cell whose + // mean is zero or negative is kept, which is what keeps the fit honest in a direction that has died. + // + // The covariance is a cluster-robust sandwich, clustering by unique reflection. That is not a detail: + // symmetry mates of one reflection share their true |F|^2, so a symmetry-forbidden tensor direction - + // determined from contrasts WITHIN a cluster - carries none of the Wilson scatter that a + // symmetry-allowed direction carries, and a covariance that assumed independent cells would get the + // ratio of the two backwards by a factor of several. + // + // Following Popov & Bourenkov (2003) Acta Cryst. D59, 1145-1153 + TensorFit FitTensor(const Cells &c, const std::vector &basis, bool want_covariance) { + TensorFit fit; + const int nb = static_cast(basis.size()); + const int nc = static_cast(c.mu.size()); + if (nb == 0 || nc < MIN_CELLS) + return fit; + + // Shells that carry a positive mean intensity: ln K is undefined for one that does not, and a + // shell with no signal at all cannot say anything about direction either. + std::vector shell_num(c.n_shells, 0.0), shell_den(c.n_shells, 0.0); + for (int i = 0; i < nc; ++i) { + const double w = 1.0 / (c.se[i] * c.se[i]); + shell_num[c.shell[i]] += w * c.mu[i]; + shell_den[c.shell[i]] += w; + } + std::vector slot(c.n_shells, -1); + int ns = 0; + std::vector ln_k; + for (int s = 0; s < c.n_shells; ++s) + if (shell_den[s] > 0.0 && shell_num[s] / shell_den[s] > 0.0) { + slot[s] = ns++; + ln_k.push_back(std::log(shell_num[s] / shell_den[s])); + } + std::vector keep; + for (int i = 0; i < nc; ++i) + if (slot[c.shell[i]] >= 0) + keep.push_back(i); + if (static_cast(keep.size()) < MIN_CELLS || ns == 0) + return fit; + + const int n = static_cast(keep.size()); + const int np = ns + nb; + Eigen::MatrixXd G(n, nb); + Eigen::VectorXd mu(n), se(n); + std::vector row_slot(n); + for (int i = 0; i < n; ++i) { + const int ci = keep[i]; + for (int k = 0; k < nb; ++k) + G(i, k) = -0.5 * QuadForm(basis[k], c.q[ci]); + mu(i) = c.mu[ci]; + se(i) = c.se[ci]; + row_slot[i] = slot[c.shell[ci]]; + } + + Eigen::VectorXd beta = Eigen::VectorXd::Zero(np); + for (int s = 0; s < ns; ++s) + beta(s) = ln_k[s]; + + auto predict = [&](const Eigen::VectorXd &b, Eigen::VectorXd &pred) { + for (int i = 0; i < n; ++i) { + double e = b(row_slot[i]); + for (int k = 0; k < nb; ++k) + e += G(i, k) * b(ns + k); + pred(i) = std::exp(std::clamp(e, -60.0, 60.0)); + } + }; + auto chi2_of = [&](const Eigen::VectorXd &pred) { + double s = 0.0; + for (int i = 0; i < n; ++i) { + const double r = (mu(i) - pred(i)) / se(i); + s += r * r; + } + return s; + }; + + Eigen::VectorXd pred(n), resid(n); + predict(beta, pred); + double chi2 = chi2_of(pred); + Eigen::MatrixXd J(n, np); + for (int iter = 0; iter < 60; ++iter) { + J.setZero(); + for (int i = 0; i < n; ++i) { + const double f = -pred(i) / se(i); + J(i, row_slot[i]) = f; + for (int k = 0; k < nb; ++k) + J(i, ns + k) = f * G(i, k); + resid(i) = (mu(i) - pred(i)) / se(i); + } + const Eigen::MatrixXd A = J.transpose() * J; + const Eigen::VectorXd grad = J.transpose() * resid; + const Eigen::VectorXd step = + A.ldlt().solve(-grad); // Gauss-Newton on chi2/2 + if (!step.allFinite()) + break; + // Halve the step until chi2 falls; a fit that cannot improve at all is converged. + double t = 1.0; + bool improved = false; + Eigen::VectorXd trial(np), trial_pred(n); + for (int back = 0; back < 20; ++back) { + trial = beta + t * step; + predict(trial, trial_pred); + const double trial_chi2 = chi2_of(trial_pred); + if (trial_chi2 < chi2) { + const bool converged = chi2 - trial_chi2 < 1e-10 * std::max(1.0, chi2); + beta = trial; + pred = trial_pred; + chi2 = trial_chi2; + improved = !converged; + break; + } + t *= 0.5; + } + if (!improved) + break; + } + + Tensor b{0, 0, 0, 0, 0, 0}; + for (int k = 0; k < nb; ++k) + b = b + Scale(basis[k], beta(ns + k)); + fit.b = b; + fit.chi2red = chi2 / std::max(1, n - np); + fit.ok = true; + if (!want_covariance) + return fit; + + // Sandwich: A^-1 B A^-1 with A the Gauss-Newton normal matrix and B the sum of outer products of + // per-cluster scores. A cell mean is the mean of its reflections, so a reflection's share of the + // cell's score is its own share of the cell's residual. + for (int i = 0; i < n; ++i) { + const double f = -pred(i) / se(i); + J(i, row_slot[i]) = f; + for (int k = 0; k < nb; ++k) + J(i, ns + k) = f * G(i, k); + } + const Eigen::MatrixXd A = J.transpose() * J; + std::vector cell_row(c.mu.size(), -1); + for (int i = 0; i < n; ++i) + cell_row[keep[i]] = i; + Eigen::MatrixXd B = Eigen::MatrixXd::Zero(np, np); + Eigen::VectorXd score = Eigen::VectorXd::Zero(np); + int32_t current = -1; + auto flush = [&]() { + if (current >= 0) + B.noalias() += score * score.transpose(); + score.setZero(); + }; + for (size_t e = 0; e < c.entry_cell.size(); ++e) { + if (c.entry_cluster[e] != current) { + flush(); + current = c.entry_cluster[e]; + } + const int i = cell_row[c.entry_cell[e]]; + if (i < 0) + continue; + const int ci = keep[i]; + const double w = (c.entry_value[e] - pred(i)) * pred(i) + / (c.n[ci] * c.se[ci] * c.se[ci]); + score(row_slot[i]) += w; + for (int k = 0; k < nb; ++k) + score(ns + k) += w * G(i, k); + } + flush(); + const Eigen::MatrixXd ainv = A.ldlt().solve(Eigen::MatrixXd::Identity(np, np)); + const Eigen::MatrixXd full = ainv * B * ainv; + fit.cov = full.bottomRightCorner(nb, nb); + return fit; + } + + struct Eigen3 { + double value[3]; + double vec[3][3]; + }; + + Eigen3 Decompose(const Tensor &b) { + Eigen3 out{}; + double d[3]; + const gemmi::Mat33 v = gemmi::eigen_decomposition(b, d); + int idx[3] = {0, 1, 2}; + std::sort(idx, idx + 3, [&](int a, int c) { return d[a] > d[c]; }); + for (int n = 0; n < 3; ++n) { + out.value[n] = d[idx[n]]; + for (int j = 0; j < 3; ++j) + out.vec[n][j] = v[j][idx[n]]; // eigenvectors are the columns + } + return out; + } + + double DeltaB(const Tensor &b) { + double d[3]; + gemmi::eigen_decomposition(b, d); + return *std::max_element(d, d + 3) - *std::min_element(d, d + 3); + } + + // The dB that counting statistics alone would produce in a subspace: E[dB] = k * sigma, with sigma the + // per-component standard error and k the constant for that subspace's dimension. + double CountingNull(const Eigen::MatrixXd &cov, double k) { + if (cov.rows() == 0) + return 0.0; + return k * std::sqrt(std::max(0.0, cov.trace() / cov.rows())); + } + + // ---------------------------------------------------------------- the s^2 signature + struct ShapeFit { + int shells = 0; + double c0 = NAN, c0_err = NAN, c1 = NAN, slope_through_origin = NAN; + double curvature = NAN, curvature_err = NAN, curvature_share = NAN; + double residual = NAN; + bool flat = false, convex = false; + }; + + // One shell's l=2 amplitude: mu_cell = K exp(A u_cell), u = -1/2 (s^T bhat s)/|s|^2, K profiled out, + // fitted on the intensity scale so no cell is dropped. A genuine tensor gives A = |B_dev| s^2, i.e. a + // straight line through the origin; the fit makes no assumption at all about how A depends on s. + bool FitShellAmplitude(const std::vector &u, const std::vector &mu, + const std::vector &se, double &a_out, double &a_err) { + auto chi2 = [&](double a) { + double num = 0.0, den = 0.0; + std::vector g(u.size()); + for (size_t i = 0; i < u.size(); ++i) { + g[i] = std::exp(std::clamp(a * u[i], -60.0, 60.0)); + const double w = 1.0 / (se[i] * se[i]); + num += mu[i] * g[i] * w; + den += g[i] * g[i] * w; + } + const double k = den > 0.0 ? num / den : 0.0; + double s = 0.0; + for (size_t i = 0; i < u.size(); ++i) { + const double r = (mu[i] - k * g[i]) / se[i]; + s += r * r; + } + return s; + }; + // A coarse scan brackets the minimum, golden section refines it. A is |B_dev| s^2, so a few + // hundred covers any crystal; the search is on a smooth one-parameter curve. + double lo = -300.0, hi = 300.0, best = 0.0, best_chi2 = chi2(0.0); + for (int i = 0; i <= 120; ++i) { + const double a = lo + (hi - lo) * i / 120.0; + const double v = chi2(a); + if (v < best_chi2) { best_chi2 = v; best = a; } + } + lo = best - 5.0; + hi = best + 5.0; + constexpr double INV_PHI = 0.6180339887498949; + double x1 = hi - INV_PHI * (hi - lo), x2 = lo + INV_PHI * (hi - lo); + double f1 = chi2(x1), f2 = chi2(x2); + for (int i = 0; i < 60; ++i) { + if (f1 < f2) { hi = x2; x2 = x1; f2 = f1; x1 = hi - INV_PHI * (hi - lo); f1 = chi2(x1); } + else { lo = x1; x1 = x2; f1 = f2; x2 = lo + INV_PHI * (hi - lo); f2 = chi2(x2); } + } + a_out = 0.5 * (lo + hi); + // One-parameter Gauss-Newton error bar, scaled by the residual so a shell the model does not + // describe reports a large one. + double num = 0.0, den = 0.0; + std::vector g(u.size()); + for (size_t i = 0; i < u.size(); ++i) { + g[i] = std::exp(std::clamp(a_out * u[i], -60.0, 60.0)); + const double w = 1.0 / (se[i] * se[i]); + num += mu[i] * g[i] * w; + den += g[i] * g[i] * w; + } + const double k = den > 0.0 ? num / den : 0.0; + double jtj = 0.0, ss = 0.0; + for (size_t i = 0; i < u.size(); ++i) { + const double dr = k * g[i] * u[i] / se[i]; + jtj += dr * dr; + const double r = (mu[i] - k * g[i]) / se[i]; + ss += r * r; + } + if (!(jtj > 0.0) || u.size() < 3) + return false; + a_err = std::sqrt(ss / (u.size() - 2) / jtj); + return std::isfinite(a_out) && std::isfinite(a_err) && a_err > 0.0; + } + + // Weighted least squares of `y` on the columns of `x`; returns the parameters and the covariance. + bool Wls(const Eigen::MatrixXd &x, const Eigen::VectorXd &y, const Eigen::VectorXd &w, + Eigen::VectorXd &p, Eigen::MatrixXd &cov, double &chi2) { + const Eigen::MatrixXd xtw = x.transpose() * w.asDiagonal(); + const Eigen::MatrixXd normal = xtw * x; + const Eigen::LDLT ldlt(normal); + if (ldlt.info() != Eigen::Success) + return false; + p = ldlt.solve(xtw * y); + const Eigen::VectorXd r = y - x * p; + chi2 = r.transpose() * w.asDiagonal() * r; + cov = ldlt.solve(Eigen::MatrixXd::Identity(x.cols(), x.cols())); + return p.allFinite(); + } + + ShapeFit FitShape(const Cells &c, const Tensor &tensor) { + ShapeFit out; + const double norm = std::sqrt(FrobDot(tensor, tensor)); + if (!(norm > 0.0)) + return out; + const Tensor bhat = Scale(tensor, 1.0 / norm); + + std::vector xs, as, errs; + for (int s = 0; s < c.n_shells; ++s) { + std::vector u, mu, se; + double s2sum = 0.0; + for (size_t i = 0; i < c.mu.size(); ++i) { + if (c.shell[i] != s) + continue; + u.push_back(-0.5 * QuadForm(bhat, c.q[i]) / std::max(c.s2[i], 1e-12)); + mu.push_back(c.mu[i]); + se.push_back(c.se[i]); + s2sum += c.s2[i]; + } + if (static_cast(u.size()) < MIN_CELLS_PER_SHAPE_SHELL) + continue; + const auto [lo, hi] = std::minmax_element(u.begin(), u.end()); + if (*hi - *lo < 1e-6) + continue; + double a = 0.0, err = 0.0; + if (!FitShellAmplitude(u, mu, se, a, err)) + continue; + xs.push_back(s2sum / u.size()); + as.push_back(a); + errs.push_back(err); + } + out.shells = static_cast(xs.size()); + if (out.shells < MIN_SHAPE_SHELLS) + return out; + + const int n = out.shells; + Eigen::VectorXd y(n), w(n); + for (int i = 0; i < n; ++i) { + y(i) = as[i]; + w(i) = 1.0 / (errs[i] * errs[i]); + } + Eigen::VectorXd p; + Eigen::MatrixXd cov; + double chi2 = 0.0; + Eigen::MatrixXd x1(n, 1), x2(n, 2), xq(n, 2); + for (int i = 0; i < n; ++i) { + x1(i, 0) = xs[i]; + x2(i, 0) = 1.0; + x2(i, 1) = xs[i]; + xq(i, 0) = xs[i]; + xq(i, 1) = xs[i] * xs[i]; + } + if (!Wls(x1, y, w, p, cov, chi2)) + return out; + out.slope_through_origin = p(0); + out.residual = chi2 / std::max(1, n - 1); + if (!Wls(x2, y, w, p, cov, chi2)) + return out; + out.c0 = p(0); + out.c0_err = std::sqrt(std::max(0.0, cov(0, 0))); + out.c1 = p(1); + if (Wls(xq, y, w, p, cov, chi2)) { + out.curvature = p(1); + out.curvature_err = std::sqrt(std::max(0.0, cov(1, 1))); + // How much of the deficit at the resolution limit the s^4 term carries. + const double s2max = *std::max_element(xs.begin(), xs.end()); + const double total = p(0) * s2max + p(1) * s2max * s2max; + out.curvature_share = total != 0.0 ? p(1) * s2max * s2max / total : NAN; + } + out.flat = out.c0 > FLAT_INTERCEPT && out.c0_err > 0.0 && out.c0 / out.c0_err > SHAPE_Z; + out.convex = !out.flat && out.curvature > 0.0 && out.curvature_err > 0.0 + && out.curvature / out.curvature_err > SHAPE_Z + && out.curvature_share > CONVEX_SHARE; + return out; + } + + // ---------------------------------------------------------------- directional limits + // in a cone about each principal direction, read where it falls through 2.0. This uses no + // model of the fall-off at all, which is why it and the tensor are reported side by side. + // + // Following Evans & Murshudov (2013) Acta Cryst. D69, 1204-1214 + void ConeLimits(const std::vector &merged, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup *sg, const Eigen3 &axes, int n_shells, + double (&d_min)[3]) { + const std::vector ops = sg ? sg->operations().sym_ops + : std::vector{gemmi::Op::identity()}; + std::vector s2; + s2.reserve(merged.size()); + for (const auto &r : merged) + if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f && r.d > 0.0f) + s2.push_back(1.0 / (static_cast(r.d) * r.d)); + if (s2.size() < static_cast(n_shells * MIN_CELL_REFLECTIONS)) + return; + std::vector sorted = s2; + std::sort(sorted.begin(), sorted.end()); + std::vector edge(n_shells + 1); + for (int i = 0; i <= n_shells; ++i) + edge[i] = sorted[std::min(sorted.size() - 1, sorted.size() * i / n_shells)]; + + const double cos_cone = std::cos(CONE_HALF_ANGLE_DEG * gemmi::pi() / 180.0); + std::vector> isig_sum(3, std::vector(n_shells, 0.0)); + std::vector> isig_n(3, std::vector(n_shells, 0)); + size_t iu = 0; + for (const auto &r : merged) { + if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || !(r.d > 0.0f)) + continue; + const double this_s2 = s2[iu++]; + int shell = 0; + while (shell + 1 < n_shells && this_s2 > edge[shell + 1]) + ++shell; + bool in_cone[3] = {false, false, false}; + for (const auto &op : ops) { + const gemmi::Op::Miller m = op.apply_to_hkl({{r.h, r.k, r.l}}); + const gemmi::Vec3 s = cell.frac.mat.left_multiply(gemmi::Vec3(m[0], m[1], m[2])); + const double len = s.length(); + if (!(len > 0)) + continue; + for (int n = 0; n < 3; ++n) { + const double dot = (s.x * axes.vec[n][0] + s.y * axes.vec[n][1] + s.z * axes.vec[n][2]) / len; + if (std::fabs(dot) >= cos_cone) + in_cone[n] = true; + } + } + for (int n = 0; n < 3; ++n) + if (in_cone[n]) { + isig_sum[n][shell] += r.I / r.sigma; + isig_n[n][shell] += 1; + } + } + for (int n = 0; n < 3; ++n) { + double last_s2 = NAN, last_isig = NAN, limit_s2 = NAN; + for (int s = 0; s < n_shells; ++s) { + if (isig_n[n][s] < MIN_CELL_REFLECTIONS) + continue; + const double mid = 0.5 * (edge[s] + edge[s + 1]); + const double isig = isig_sum[n][s] / isig_n[n][s]; + if (isig >= CONE_I_OVER_SIGMA) { + last_s2 = mid; + last_isig = isig; + limit_s2 = edge[s + 1]; + } else if (std::isfinite(last_s2)) { + // Linear interpolation of in s^2 between the last shell above the threshold + // and the first below it. + const double f = (last_isig - CONE_I_OVER_SIGMA) / (last_isig - isig); + limit_s2 = last_s2 + f * (mid - last_s2); + break; + } else { + break; + } + } + if (std::isfinite(limit_s2) && limit_s2 > 0.0) + d_min[n] = 1.0 / std::sqrt(limit_s2); + } + } + + const char *Band(double g) { + if (g < GATE_MARGINAL) return "not established"; + if (g < GATE_ESTABLISHED) return "marginal"; + if (g < GATE_STRONG) return "established"; + return "strong"; + } + + // ---------------------------------------------------------------- the systematic-error floor + // The same cells, built from UNMERGED observations at the Miller index each was measured at, with + // the unique reflection they reduce to as the cluster. Nothing is symmetry-expanded and nothing is + // de-duplicated: two observations of the same reflection recorded on different frames are two + // measurements, and the difference between them is exactly the signal this arm is after. + Cells BuildObservationCells(const std::vector &obs, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup *sg, int n_shells) { + Cells c; + c.n_shells = n_shells; + const HKLKeyGenerator key(/*merge_friedel=*/true, sg ? *sg : *gemmi::find_spacegroup_by_number(1)); + + // Order the observations by the unique reflection they belong to, so the cluster-robust sum + // below can walk them one cluster at a time; then stride whole clusters if there are more than + // the fit needs. The estimate is flat in the number of unique reflections over a factor of 12, + // so striding costs nothing and bounds the work. + std::vector> order; + order.reserve(obs.size()); + for (size_t i = 0; i < obs.size(); ++i) { + const auto &o = obs[i]; + if (!std::isfinite(o.I) || !std::isfinite(o.sigma) || o.sigma <= 0.0f || !(o.d > 0.0f)) + continue; + order.emplace_back(key(o.h, o.k, o.l).pack(), static_cast(i)); + } + if (order.empty()) + return c; + std::sort(order.begin(), order.end()); + int32_t n_clusters = 1; + for (size_t i = 1; i < order.size(); ++i) + if (order[i].first != order[i - 1].first) + ++n_clusters; + const int cluster_stride = std::max(1, n_clusters / MAX_CLUSTERS); + + const auto dirs = HemisphereDirections(N_DIRECTIONS); + const int n_cells_total = n_shells * N_DIRECTIONS; + std::vector sum(n_cells_total, 0.0), sumsq(n_cells_total, 0.0), s2sum(n_cells_total, 0.0); + std::vector> qsum(n_cells_total, {0, 0, 0, 0, 0, 0}); + std::vector count(n_cells_total, 0); + std::vector raw_cell, raw_cluster; + std::vector raw_value; + + std::vector kept; + std::vector cluster_of; + uint64_t previous_key = 0; + int32_t seen = -1; + for (const auto &[k, idx] : order) { + if (seen < 0 || k != previous_key) { + previous_key = k; + ++seen; + } + if (seen % cluster_stride != 0) + continue; + kept.push_back(idx); + cluster_of.push_back(seen / cluster_stride); + } + if (kept.size() < static_cast(MIN_CELLS * MIN_CELL_REFLECTIONS)) + return c; + + // Equal-count shells over the observations that are kept. + std::vector s2(kept.size()); + for (size_t i = 0; i < kept.size(); ++i) { + const auto &o = obs[kept[i]]; + s2[i] = 1.0 / (static_cast(o.d) * o.d); + } + std::vector by_s2(kept.size()); + std::iota(by_s2.begin(), by_s2.end(), 0); + std::sort(by_s2.begin(), by_s2.end(), [&](int32_t a, int32_t b) { return s2[a] < s2[b]; }); + std::vector shell(kept.size()); + for (size_t i = 0; i < by_s2.size(); ++i) + shell[by_s2[i]] = std::min(n_shells - 1, static_cast(i * n_shells / by_s2.size())); + + raw_cell.reserve(kept.size()); + raw_value.reserve(kept.size()); + raw_cluster.reserve(kept.size()); + for (size_t i = 0; i < kept.size(); ++i) { + const auto &o = obs[kept[i]]; + gemmi::Vec3 s = cell.frac.mat.left_multiply(gemmi::Vec3(o.h, o.k, o.l)); + if (s.z < 0) + s = gemmi::Vec3(-s.x, -s.y, -s.z); + const double len = s.length(); + if (!(len > 0)) + continue; + const gemmi::Vec3 unit(s.x / len, s.y / len, s.z / len); + int best = 0; + double best_dot = -1.0; + for (int j = 0; j < N_DIRECTIONS; ++j) { + const double dot = std::fabs(unit.dot(dirs[j])); + if (dot > best_dot) { best_dot = dot; best = j; } + } + const int cid = shell[i] * N_DIRECTIONS + best; + sum[cid] += o.I; + sumsq[cid] += static_cast(o.I) * o.I; + s2sum[cid] += s2[i]; + qsum[cid][0] += s.x * s.x; qsum[cid][1] += s.y * s.y; qsum[cid][2] += s.z * s.z; + qsum[cid][3] += s.x * s.y; qsum[cid][4] += s.x * s.z; qsum[cid][5] += s.y * s.z; + count[cid] += 1; + raw_cell.push_back(cid); + raw_cluster.push_back(cluster_of[i]); + raw_value.push_back(o.I); + } + FinishCells(c, sum, sumsq, s2sum, qsum, count, raw_cell, raw_cluster, raw_value); + return c; + } + + // The systematic error scale this dataset carries, and the part of it that is merely counting + // noise, both measured in the tensor directions the Laue class forbids. + struct SystematicFloor { + bool ok = false; + double sigma_excess = NAN; // per tensor component, A^2 + int n_observations = 0; + }; + + SystematicFloor MeasureSystematicFloor(const std::vector &obs, + const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg, + const std::vector &full, + const std::vector &forbidden) { + SystematicFloor out; + const int q = static_cast(forbidden.size()); + if (q == 0 || obs.empty()) + return out; + const Cells cells = BuildObservationCells(obs, cell, sg, N_SHELLS); + if (cells.mu.empty()) + return out; + const TensorFit fit = FitTensor(cells, full, true); + if (!fit.ok || fit.cov.rows() != static_cast(full.size())) + return out; + + Eigen::MatrixXd basis_map(static_cast(full.size()), q); + for (int j = 0; j < q; ++j) + for (size_t k = 0; k < full.size(); ++k) + basis_map(static_cast(k), j) = FrobDot(full[k], forbidden[j]); + Tensor forb{0, 0, 0, 0, 0, 0}; + Eigen::VectorXd theta(static_cast(full.size())); + for (size_t k = 0; k < full.size(); ++k) + theta(static_cast(k)) = FrobDot(fit.b, full[k]); + const Eigen::VectorXd theta_forb = basis_map.transpose() * theta; + for (int j = 0; j < q; ++j) + forb = forb + Scale(forbidden[j], theta_forb(j)); + const Eigen::MatrixXd cov_forb = basis_map.transpose() * fit.cov * basis_map; + + const double k_q = K_FORBIDDEN[q]; + if (!(k_q > 0.0)) + return out; + const double sigma_sys = DeltaB(forb) / k_q; + const double sigma_stat = CountingNull(cov_forb, k_q) / k_q; + // Take the counting part out: what is left is the systematic alone. + out.sigma_excess = std::sqrt(std::max(0.0, sigma_sys * sigma_sys - sigma_stat * sigma_stat)); + out.n_observations = static_cast(cells.entry_cell.size()); + out.ok = true; + return out; + } +} + +std::vector ScaledObservations(const std::vector &outcomes, + bool rotation, double min_partiality) { + // Per-image scale, indexed the way the outcomes are. + std::vector g(outcomes.size(), 0.0); + for (size_t i = 0; i < outcomes.size(); ++i) + if (outcomes[i].image_scale_g.has_value() && std::isfinite(*outcomes[i].image_scale_g) + && *outcomes[i].image_scale_g > 0.0f) + g[i] = *outcomes[i].image_scale_g; + + std::vector out; + // The same acceptance the merge itself uses (Merge.cpp): everything in an IntegrationOutcome came + // out of the integrator, and Reflection::observed is not set by the _process.h5 reader, so testing + // it would silently empty this on the --mode scale path. + auto usable = [](const Reflection &r) { + return std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f && r.d > 0.0f + && std::isfinite(r.rlp) && r.rlp > 0.0f && std::isfinite(r.partiality); + }; + + if (!rotation) { + // A still is a whole measurement of its reflection, and consecutive stills are different + // crystals, so there is nothing to assemble. + for (size_t i = 0; i < outcomes.size(); ++i) { + if (!(g[i] > 0.0)) + continue; + for (const auto &r : outcomes[i].reflections) { + if (!usable(r) || !(r.partiality >= min_partiality)) + continue; + const double corr = r.rlp / (r.partiality * g[i]); + out.push_back({r.h, r.k, r.l, static_cast(r.I * corr), + static_cast(r.sigma * corr), r.d}); + } + } + return out; + } + + // A rotation reflection is integrated image by image, so it arrives as a run of partials over + // consecutive frames. Assemble each run into one full the way the 3D combine and the unmerged + // export do - same raw hkl, frames no further apart than the combine's own gap, parts added with + // their variances in quadrature - because a single partial divided by its own partiality carries + // the rocking-curve model's error as well as its intensity, and that error is a function of the + // reflection's direction relative to the spindle, which is precisely the direction this + // diagnostic is measuring. + constexpr float MAX_FRAME_GAP = 2.0f; + struct Part { + const Reflection *r; + size_t image; + }; + std::vector parts; + for (size_t i = 0; i < outcomes.size(); ++i) { + if (!(g[i] > 0.0)) + continue; + for (const auto &r : outcomes[i].reflections) + if (usable(r)) + parts.push_back({&r, i}); + } + std::sort(parts.begin(), parts.end(), [](const Part &a, const Part &b) { + return std::tie(a.r->h, a.r->k, a.r->l, a.r->image_number) + < std::tie(b.r->h, b.r->k, b.r->l, b.r->image_number); + }); + for (size_t i = 0; i < parts.size();) { + size_t j = i + 1; + while (j < parts.size() && parts[j].r->h == parts[i].r->h && parts[j].r->k == parts[i].r->k + && parts[j].r->l == parts[i].r->l + && parts[j].r->image_number - parts[j - 1].r->image_number <= MAX_FRAME_GAP) + ++j; + double sum_p = 0.0, sum_I = 0.0, sum_var = 0.0, p_g = 0.0; + for (size_t m = i; m < j; ++m) { + const Reflection &r = *parts[m].r; + sum_p += r.partiality; + sum_I += static_cast(r.I) * r.rlp; + sum_var += static_cast(r.sigma) * r.sigma * r.rlp * r.rlp; + p_g += r.partiality * g[parts[m].image]; + } + const Reflection &first = *parts[i].r; + i = j; + // The scale is the rocking curve's own weighted mean over the frames the event spans; the + // partiality divisor only corrects an event the sweep cut short, since a complete one sums + // to 1 by construction. + if (!(sum_p >= min_partiality) || !(p_g > 0.0)) + continue; + const double scale = p_g / sum_p; // partiality-weighted mean G + const double corr = 1.0 / (sum_p * scale); + out.push_back({first.h, first.k, first.l, static_cast(sum_I * corr), + static_cast(std::sqrt(sum_var) * corr), first.d}); + } + return out; +} + +const char *AnisotropyShapeCode(AnisotropyShape shape) { + switch (shape) { + case AnisotropyShape::Linear: return "LINEAR"; + case AnisotropyShape::Flat: return "FLAT"; + case AnisotropyShape::Convex: return "CONVEX"; + default: return "UNDETERMINED"; + } +} + +const char *AnisotropyVerdictCode(AnisotropyVerdict verdict) { + switch (verdict) { + case AnisotropyVerdict::NotDetected: return "NOT_DETECTED"; + case AnisotropyVerdict::Detected: return "DETECTED"; + default: return "CANNOT_DETERMINE"; + } +} + +AnisotropyResult AnalyzeAnisotropy(const std::vector &merged, + const std::vector &unmerged, + const gemmi::UnitCell &cell, + const gemmi::SpaceGroup *space_group, + const AnisotropyRunInfo &run) { + AnisotropyResult result; + if (merged.empty() || !cell.is_crystal()) + return result; + + const auto allowed = AllowedDeviatoricBasis(cell, space_group); + const auto full = FullDeviatoricBasis(); + const auto forbidden = ForbiddenBasis(full, allowed); + const int p = static_cast(allowed.size()); + const int q = static_cast(forbidden.size()); + result.n_free_parameters = p; + + if (p == 0) { + // Cubic. The anisotropy tensor is forced isotropic by symmetry, so its deviatoric part is exactly + // zero - there is nothing to fit and nothing to gate. + result.n_reflections = static_cast(merged.size()); + result.delta_b = 0.0; + result.delta_b_linear = 0.0; + result.delta_b_flat = 0.0; + result.eigenvalue[0] = result.eigenvalue[1] = result.eigenvalue[2] = 0.0; + result.verdict = AnisotropyVerdict::NotDetected; + result.band = "no anisotropy is possible in this Laue class"; + return result; + } + + double isig_sum = 0.0; + int isig_n = 0; + for (const auto &r : merged) + if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) { + isig_sum += r.I / r.sigma; + ++isig_n; + } + const double mean_isig = isig_n > 0 ? isig_sum / isig_n : 0.0; + + const auto refl = SelectReflections(merged, cell, space_group, N_SHELLS); + if (refl.size() < static_cast(MIN_CELLS * MIN_CELL_REFLECTIONS)) + return result; + + const Cells cells = BuildCells(refl, cell, space_group, N_SHELLS); + const TensorFit constrained = FitTensor(cells, allowed, true); + if (!constrained.ok) + return result; + result.n_reflections = static_cast(refl.size()); + result.n_cells = static_cast(cells.mu.size()); + + const Eigen3 axes = Decompose(constrained.b); + result.delta_b = axes.value[0] - axes.value[2]; + for (int n = 0; n < 3; ++n) { + result.eigenvalue[n] = axes.value[n]; + for (int j = 0; j < 3; ++j) + result.eigenvector[n][j] = axes.vec[n][j]; + } + double d_min = 0.0; + for (const auto &r : merged) + if (r.d > 0.0f && (d_min == 0.0 || r.d < d_min)) + d_min = r.d; + if (d_min > 0.0) + result.fold_weakening = std::exp(result.delta_b / (2.0 * d_min * d_min)); + ConeLimits(merged, cell, space_group, axes, N_SHELLS, result.d_min_axis); + { + double lo = INFINITY, hi = -INFINITY; + for (double d : result.d_min_axis) + if (std::isfinite(d)) { lo = std::min(lo, d); hi = std::max(hi, d); } + if (std::isfinite(lo) && std::isfinite(hi)) + result.d_min_spread = hi - lo; + } + + // --- the resolution signature, and whether the binning decides it --- + const ShapeFit shape = FitShape(cells, constrained.b); + result.shape_shells = shape.shells; + result.shape_intercept = shape.c0; + result.shape_intercept_z = shape.c0_err > 0.0 ? shape.c0 / shape.c0_err : NAN; + result.shape_slope = shape.c1; + result.shape_curvature_z = shape.curvature_err > 0.0 ? shape.curvature / shape.curvature_err : NAN; + result.shape_curvature_share = shape.curvature_share; + result.shape_residual = shape.residual; + if (shape.shells >= MIN_SHAPE_SHELLS) { + result.shape = shape.flat ? AnisotropyShape::Flat + : shape.convex ? AnisotropyShape::Convex : AnisotropyShape::Linear; + if (std::isfinite(shape.slope_through_origin) && shape.slope_through_origin != 0.0) { + result.delta_b_linear = result.delta_b * shape.c1 / shape.slope_through_origin; + result.delta_b_flat = result.delta_b - result.delta_b_linear; + } + // Rebin at 8 and at 16 shells and refit the tensor each time, so the direction the amplitude is + // measured along is re-derived too. A verdict that moves is not a measurement. + result.shape_stable = true; + for (int alt : {8, 16}) { + const auto alt_refl = SelectReflections(merged, cell, space_group, alt); + const Cells alt_cells = BuildCells(alt_refl, cell, space_group, alt); + const TensorFit alt_fit = FitTensor(alt_cells, allowed, false); + if (!alt_fit.ok) { + result.shape_stable = false; + break; + } + const ShapeFit alt_shape = FitShape(alt_cells, alt_fit.b); + if (alt_shape.shells < MIN_SHAPE_SHELLS || alt_shape.flat != shape.flat) { + result.shape_stable = false; + break; + } + } + if (!result.shape_stable) + result.shape = AnisotropyShape::Undetermined; + } + // The gate is applied to the part of delta_b that behaves as a Debye-Waller B, when that could be + // separated; to delta_b itself when it could not. + // A FLAT crystal can fit a negative linear slope; a negative anisotropy is not a smaller one, it is + // no established Debye-Waller component at all. + const double magnitude = std::max(0.0, std::isfinite(result.delta_b_linear) ? result.delta_b_linear + : result.delta_b); + + // --- the gate --- + auto refuse = [&](std::string why) { + result.verdict = AnisotropyVerdict::CannotDetermine; + result.refusal = std::move(why); + }; + if (mean_isig < MIN_MEAN_I_OVER_SIGMA) { + refuse("the merged data are at the noise floor ( below 1), where neither the tensor " + "nor its resolution signature means anything"); + } else if (q == 0) { + // Triclinic. Every other Laue class leaves directions in which the true tensor is zero by + // symmetry, and measuring those is how this dataset's systematic error is established. Triclinic + // leaves none, and three candidate substitutes (a degree-4 spherical-harmonic block, a + // low/high-resolution split, a floor borrowed from similar datasets) were built and measured + // against known truth: the first two are blind to real anisotropy, as a denominator must be, but + // neither delivers a usable scale, and the borrowed floor spans three orders of magnitude. + refuse("the Laue class is triclinic, which leaves no symmetry-forbidden tensor direction, so this " + "dataset carries no measurement of its own systematic error"); + } else if (std::isfinite(run.observed_rotation_deg) && run.observed_rotation_deg < MIN_ROTATION_DEG) { + refuse("the observed rotation range is too short for the tensor to be separated from a lab-fixed " + "systematic"); + } else if (std::isfinite(run.observed_rotation_deg) && !run.dose_term_in_scale_model) { + // A dose ramp is not a significance failure - it is a confident false detection. A simulated + // isotropic crystal given a 25 A^2 relative-B ramp over 180 deg returns dB ~ 9 A^2 with a + // significance no test can see through, because the ramp really is a smooth quadratic function of + // direction once the sweep geometry is folded in. + refuse("the scale model carried no dose term, and an uncorrected dose ramp manufactures anisotropy " + "that no significance test can see through"); + } else if (unmerged.empty()) { + // Merged intensities have exact Laue symmetry by construction, so the symmetry-forbidden tensor + // directions - the only place a dataset can measure its own systematic error - are identically + // zero there whatever the crystal carries. Gating on what is left would be gating on counting + // statistics, which against known ground truth calls "real" in 15-58% of clean isotropic datasets. + refuse("the systematic-error scale can only be measured on unmerged observations, and none were " + "available for this merge"); + } else { + const SystematicFloor sys = MeasureSystematicFloor(unmerged, cell, space_group, full, forbidden); + result.n_observations = sys.n_observations; + if (!sys.ok) { + refuse("the symmetry-forbidden tensor directions of the unmerged observations could not be " + "measured, so this dataset's systematic-error scale is unknown"); + } else { + const double k_p = K_ALLOWED[p]; + result.sigma_systematic = sys.sigma_excess; + const double allowed_null = CountingNull(constrained.cov, k_p); + result.floor = std::sqrt(k_p * result.sigma_systematic * k_p * result.sigma_systematic + + allowed_null * allowed_null); + if (result.floor > 0.0) { + result.significance = magnitude / result.floor; + result.detection_limit = GATE_ESTABLISHED * result.floor; + result.band = Band(result.significance); + result.verdict = result.significance > GATE_ESTABLISHED ? AnisotropyVerdict::Detected + : AnisotropyVerdict::NotDetected; + } else { + refuse("the systematic-error floor could not be measured"); + } + } + } + + // --- cautions: measured things that bias the verdict, in either direction --- + if (result.shape == AnisotropyShape::Flat) + result.cautions.emplace_back( + "the directional deficit does not follow exp(-1/2 s^T B s), so the fitted deltaB is a fit " + "of the wrong functional form and may be an under-estimate"); + if (result.shape == AnisotropyShape::Undetermined && shape.shells >= MIN_SHAPE_SHELLS) + result.cautions.emplace_back( + "the resolution signature changed when the shells were rebinned, so it is reported as " + "undetermined rather than as a measurement"); + if (std::isfinite(run.observed_rotation_deg) && run.observed_rotation_deg < 180.0) + result.cautions.emplace_back( + "the observed rotation range is under 180 deg, which leaves a second tensor direction " + "reachable by a lab-fixed systematic and reduces the power of the test"); + if (std::isfinite(run.radiation_damage_relative_b) && std::fabs(run.radiation_damage_relative_b) > 10.0) + result.cautions.emplace_back( + "the run carries substantial radiation damage; the scale model removes its average " + "monotone part, but a non-monotone dose ramp needs --relative-b to be taken out too"); + if (p == 1) + result.cautions.emplace_back( + "this Laue class leaves a single free anisotropy direction; if the space group has been " + "assigned too high a symmetry, real anisotropy is pushed into the directions used to " + "measure the systematic error and the test is biased towards reporting none"); + return result; +} + +std::string AnisotropyToText(const AnisotropyResult &result) { + std::ostringstream os; + if (result.n_reflections == 0) + return os.str(); + os << std::fixed; + os << "Diffraction anisotropy\n"; + if (result.n_free_parameters == 0) { + os << " The Laue class forces the anisotropy tensor to be isotropic: deltaB = 0 exactly, with no\n" + << " free parameter to fit. This is symmetry, not a measurement.\n"; + return os.str(); + } + os << std::setprecision(2); + os << " Anisotropic deltaB (range of principal components) = " << result.delta_b << " A^2" + << " [" << result.n_free_parameters << " free direction" + << (result.n_free_parameters == 1 ? "" : "s") << " in this Laue class]\n"; + os << " Principal components (A^2, relative):"; + for (double v : result.eigenvalue) + os << " " << (v - result.eigenvalue[2]); + os << "\n"; + if (std::isfinite(result.fold_weakening)) + os << " Strongest / weakest direction at the resolution limit: " << result.fold_weakening << "x\n"; + if (std::isfinite(result.d_min_axis[0]) || std::isfinite(result.d_min_axis[2])) { + os << " d_min along the principal directions ( = " << CONE_I_OVER_SIGMA + << " in a " << CONE_HALF_ANGLE_DEG << " deg cone):"; + for (double d : result.d_min_axis) { + if (std::isfinite(d)) os << " " << d; + else os << " -"; + } + os << " A\n"; + } + if (result.shape_shells >= MIN_SHAPE_SHELLS) { + os << " Resolution signature of the deficit over " << result.shape_shells << " shells: " + << AnisotropyShapeCode(result.shape) << " (intercept " << std::showpos << result.shape_intercept + << std::noshowpos << ", z = " << result.shape_intercept_z << "; slope " << result.shape_slope + << " A^2; through-origin chi2/dof " << result.shape_residual << ")\n"; + if (result.shape == AnisotropyShape::Flat) + os << " The deficit does not follow exp(-1/2 s^T B s): " << result.delta_b_flat + << " A^2 of the deltaB above does not behave as a Debye-Waller B. The fall-off is being\n" + << " described by the wrong functional form, so the number may be an UNDER-estimate,\n" + << " not an over-estimate.\n"; + else if (result.shape == AnisotropyShape::Convex) + os << " The deficit grows faster than s^2, which a Debye-Waller B cannot do.\n"; + } + if (result.verdict == AnisotropyVerdict::CannotDetermine) { + os << " => CANNOT DETERMINE: " << result.refusal << ".\n"; + } else { + // The magnitude judged is the Debye-Waller part of the deltaB, which is what the floor is a + // floor on; on a crystal whose deficit is not a Debye-Waller fall-off the two differ. + const double judged = std::max(0.0, std::isfinite(result.delta_b_linear) ? result.delta_b_linear + : result.delta_b); + os << std::setprecision(2) + << (result.verdict == AnisotropyVerdict::Detected ? " => DETECTED (" : " => NOT DETECTED (") + << result.band << "): a Debye-Waller deltaB of " << judged + << " A^2 against this dataset's own systematic-error floor of " << result.floor + << " A^2, ratio " << result.significance << ".\n"; + } + if (std::isfinite(result.detection_limit)) + os << " Below about " << result.detection_limit << " A^2 nothing could be established on this\n" + << " dataset. That limit is set by the systematic error, not by counting, so it does NOT\n" + << " improve with more reflections or a longer exposure.\n"; + for (const auto &c : result.cautions) + os << " Note: " << c << "\n"; + os << " Reported and not corrected: no intensity is changed and no reflection is removed.\n"; + return os.str(); +} diff --git a/image_analysis/scale_merge/AnisotropyAnalysis.h b/image_analysis/scale_merge/AnisotropyAnalysis.h new file mode 100644 index 00000000..7d650975 --- /dev/null +++ b/image_analysis/scale_merge/AnisotropyAnalysis.h @@ -0,0 +1,154 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "../../common/Reflection.h" +#include "../IntegrationOutcome.h" +#include "gemmi/symmetry.hpp" +#include "gemmi/unitcell.hpp" + +// Diffraction-anisotropy diagnostic. It DESCRIBES and it REPORTS: no intensity is corrected, no +// reflection is removed, and the written reflection file does not depend on direction in any way. +// Two quantities are reported, and they are not the same thing: +// +// - the anisotropic deltaB, the range of the principal components of the anisotropy tensor (the +// AIMLESS quantity), which is a RATE of fall-off; +// - the diffraction limit along each principal direction, from in a cone about that +// direction (the AIMLESS construction), which is where the signal actually runs out. +// +// A crystal can have a large deltaB and almost no spread in directional limit, or the reverse. +// +// CONVENTION. The tensor is fitted as ln = c(shell) - 1/2 s^T B s, s = 1/d in A^-1, B in A^2. +// For an isotropic B this is the ordinary Wilson plot, so B here is the ordinary crystallographic +// (structure-factor-level, B = 8 pi^2 U) B - the same scale as phenix.xtriage's B_cart, ctruncate's +// anisotropic B eigenvalues and AIMLESS's anisotropic deltaB, and directly comparable with all three. +// Only the DEVIATORIC part is fitted: the isotropic part is degenerate with the overall scale and with +// the Wilson curve, and one free constant per resolution shell absorbs it exactly. +// +// The tensor is fitted on INTENSITIES, with no positivity cut. Fitting amplitudes, or dropping +// non-positive intensities as the amplitude-based tools must, loses roughly 40% of the signal: in a +// direction that has died, half the merged intensities are negative, so a positivity cut keeps only +// the positive noise excursions and flattens the fall-off exactly where the anisotropy is largest. +// +// TWO INPUTS, because they answer different halves of the question. The tensor, its resolution +// signature and the directional limits are measured on the MERGED reflections. The error bar the +// verdict is gated on cannot be: merged data has exact Laue symmetry by construction, so the +// symmetry-forbidden tensor directions - where the true tensor is zero whatever the crystal is, and +// which are therefore this dataset's own measurement of its systematic error - are identically zero +// there. That measurement needs the UNMERGED, scaled observations, which still carry the differences +// between symmetry mates. Without them the verdict is "cannot determine" and says so. +// +// Method credits are at the algorithms in the .cpp. + +enum class AnisotropyShape { + Undetermined, // too few shells, or the verdict moved when the binning changed + Linear, // the deficit follows exp(-1/2 s^T B s): a Debye-Waller B + Flat, // the deficit does not follow exp(-1/2 s^T B s) - see the note on Flat below + Convex // the deficit grows faster than s^2 +}; + +enum class AnisotropyVerdict { + CannotDetermine, // this dataset cannot answer the question - `refusal` says why + NotDetected, // anisotropy is not established above this dataset's own systematic error + Detected // anisotropy is established; `band` grades it +}; + +const char *AnisotropyShapeCode(AnisotropyShape shape); // e.g. "LINEAR" +const char *AnisotropyVerdictCode(AnisotropyVerdict verdict); // e.g. "NOT_DETECTED" + +struct AnisotropyResult { + // "Could not decide" is a counter, not a flag: every renderer tests n_reflections and falls silent + // when it is zero (the TwinningAnalysis idiom). n_free_parameters is the deviatoric degrees of + // freedom the Laue class allows - 5 triclinic, 3 monoclinic, 2 orthorhombic, 1 tetragonal/trigonal/ + // hexagonal, 0 cubic - so 0 means "no anisotropy is possible here", not "nothing was measured", and + // a cubic crystal reports n_reflections > 0 with n_cells = 0 and deltaB exactly zero. + int n_reflections = 0; + int n_cells = 0; + int n_free_parameters = 0; + int n_observations = 0; // unmerged observations the systematic-error floor was measured on + + // --- the tensor --------------------------------------------------------------------------- + double delta_b = NAN; // B_max - B_min of the deviatoric tensor, A^2 + double eigenvalue[3] = {NAN, NAN, NAN}; // descending; deviatoric, so they sum to zero + double eigenvector[3][3] = {{NAN, NAN, NAN}, {NAN, NAN, NAN}, {NAN, NAN, NAN}}; // Cartesian rows + double fold_weakening = NAN; // exp(delta_b / 2 d_min^2): strong/weak intensity ratio at d_min + + // --- the diffraction limits (model-free) -------------------------------------------------- + // Highest resolution at which in a 20 deg cone about principal axis n is still 2.0. + double d_min_axis[3] = {NAN, NAN, NAN}; + double d_min_spread = NAN; // max - min of the above + + // --- the resolution signature of the deficit ---------------------------------------------- + // The l=2 directional deficit A of each shell against s^2, fitted as A = c0 + c1 s^2. A genuine + // Debye-Waller B gives a straight line through the origin (c0 = 0, c1 = |B_dev|). + int shape_shells = 0; + double shape_intercept = NAN; // c0, dimensionless + double shape_intercept_z = NAN; // c0 / sigma(c0) + double shape_slope = NAN; // c1, A^2 + double shape_curvature_z = NAN; // significance of an s^4 term + double shape_curvature_share = NAN; // fraction of the deficit that term carries at the limit + double shape_residual = NAN; // chi2/dof of the through-origin line + AnisotropyShape shape = AnisotropyShape::Undetermined; + bool shape_stable = false; // the verdict survived rebinning at 8 and 16 shells + // deltaB implied by the s^2 slope alone - what a Debye-Waller B accounts for - and the rest. + // The first can exceed delta_b when the fitted curve passes below the origin. + double delta_b_linear = NAN; + double delta_b_flat = NAN; + + // --- the gate ----------------------------------------------------------------------------- + // Everything here is measured from this dataset alone. sigma_systematic is the per-component + // systematic error scale left in the symmetry-FORBIDDEN tensor directions of the unmerged + // observations after their own counting noise has been taken out; floor is what that systematic, + // plus counting noise, would manufacture in the symmetry-ALLOWED directions; significance is + // delta_b_linear / floor. + double sigma_systematic = NAN; + double floor = NAN; + double significance = NAN; + double detection_limit = NAN; // the smallest delta_b this dataset could establish, A^2 + AnisotropyVerdict verdict = AnisotropyVerdict::CannotDetermine; + std::string band; // "not established" / "marginal" / "established" / "strong" + std::string refusal; // why, when the verdict is CannotDetermine + std::vector cautions; // things that bias the verdict, each measured +}; + +// One scaled, unmerged observation: the Miller index it was measured at (NOT reduced to the +// asymmetric unit - the whole point is that symmetry mates are separate measurements) and its +// intensity on the merge's own scale. +struct AnisotropyObservation { + int32_t h = 0, k = 0, l = 0; + float I = NAN; + float sigma = NAN; + float d = NAN; +}; + +// Build those from the per-image integration outcomes, on the scale rugnux itself applied. On rotation +// data each reflection's partials are first assembled into one full, as the 3D combine and the +// unmerged export do; on stills each reflection is already a whole measurement. Images with no fitted +// per-image scale are left out, and so is an event that caught less than min_partiality of its rocking +// curve - the combine does not merge one either. +std::vector ScaledObservations(const std::vector &outcomes, + bool rotation, double min_partiality = 0.5); + +// What the caller knows about the run and the merge that the reflections alone do not say. +struct AnisotropyRunInfo { + // Rotation actually covered by the processed images, in degrees (NOT the header's nominal sweep - + // one battery dataset spans 179 deg against a 360 deg header). NaN for stills / not known. + double observed_rotation_deg = NAN; + // Whether the scale model that produced these intensities carried a dose (decay) term. Mandatory + // on rotation data: see the note at the refusal in the .cpp. + bool dose_term_in_scale_model = true; + // The run's measured radiation damage (MergeStatistics::radiation_damage_delta_b), for the caution. + double radiation_damage_relative_b = NAN; +}; + +AnisotropyResult AnalyzeAnisotropy(const std::vector &merged, + const std::vector &unmerged, + const gemmi::UnitCell &cell, + const gemmi::SpaceGroup *space_group, + const AnisotropyRunInfo &run = {}); + +std::string AnisotropyToText(const AnisotropyResult &result); diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index 338b1290..5c9dd34b 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -3,6 +3,8 @@ ADD_LIBRARY(JFJochScaleMerge SearchSpaceGroup.h TwinningAnalysis.cpp TwinningAnalysis.h + AnisotropyAnalysis.cpp + AnisotropyAnalysis.h Merge.cpp Merge.h ScaleOnTheFly.cpp diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 4635a92a..0d85450e 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -14,6 +14,7 @@ #include "../../common/Reflection.h" #include "../IntegrationOutcome.h" +#include "AnisotropyAnalysis.h" #include "HKLKey.h" struct MergeStatisticsShell { @@ -104,6 +105,12 @@ struct MergeStatistics { // 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. SweepQuality sweep_quality; + + // Diffraction anisotropy (AnalyzeAnisotropy): the anisotropy tensor, the diffraction limit along + // each of its principal directions, and whether either is established above this dataset's own + // systematic error. Report-only - nothing is corrected and no reflection is removed. Empty + // (n_reflections = 0) when the diagnostic did not run. + AnisotropyResult anisotropy; }; diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 6157dfce..e5168e7b 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1,6 +1,7 @@ // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include #include #include @@ -8,6 +9,7 @@ #include "../common/GitInfo.h" #include "../common/time_utc.h" +#include "../image_analysis/scale_merge/AnisotropyAnalysis.h" #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" #include "../image_analysis/scale_merge/TwinningAnalysis.h" @@ -17,7 +19,7 @@ namespace { // The version of this file format. Bumped when a key is renamed or removed, a table column moves, // or a reason code changes meaning - a consumer can gate on it. - constexpr int REPORT_VERSION = 1; + constexpr int REPORT_VERSION = 2; const char *BANNER = " ******************************************************************************"; @@ -263,11 +265,51 @@ std::string RenderResultReport(const std::string &output_prefix, } os << " ----------- ----------- --------- -------- -------------------- -------- ------ ------ --------\n"; - // --------------------------------------------------------------- 9. WARNINGS + // ---------------------------------------------------------- 9. DIFFRACTION ANISOTROPY + const auto &an = result.merge_statistics.anisotropy; + if (merged && an.n_reflections > 0) { + Section(os, "9. DIFFRACTION ANISOTROPY"); + os << " How much the fall-off depends on direction, and whether that is established above this\n" + << " data set's own systematic error. Nothing here corrects an intensity or removes a\n" + << " reflection: the merged data and the written files do not depend on direction at all.\n" + << " ANISOTROPY_DELTA_B is the range of the principal components of the anisotropy tensor,\n" + << " on the ordinary crystallographic B scale (the same scale as phenix.xtriage's B_cart and\n" + << " ctruncate's anisotropic B), fitted on intensities with nothing dropped.\n\n"; + Key(os, "ANISOTROPY_VERDICT", AnisotropyVerdictCode(an.verdict)); + Key(os, "ANISOTROPY_FREE_DIRECTIONS", an.n_free_parameters); + Key(os, "ANISOTROPY_DELTA_B", fmt::format("{:.2f}", an.delta_b)); + Key(os, "ANISOTROPY_DELTA_B_LINEAR", fmt::format("{:.2f}", an.delta_b_linear)); + Key(os, "ANISOTROPY_PRINCIPAL_B", fmt::format("{:.2f} {:.2f} {:.2f}", + an.eigenvalue[0] - an.eigenvalue[2], + an.eigenvalue[1] - an.eigenvalue[2], 0.0)); + Key(os, "ANISOTROPY_FOLD_WEAKENING", fmt::format("{:.1f}", an.fold_weakening)); + Key(os, "ANISOTROPY_D_MIN_PRINCIPAL", fmt::format("{:.2f} {:.2f} {:.2f}", an.d_min_axis[0], + an.d_min_axis[1], an.d_min_axis[2])); + Key(os, "ANISOTROPY_D_MIN_SPREAD", fmt::format("{:.2f}", an.d_min_spread)); + Key(os, "ANISOTROPY_SHAPE", AnisotropyShapeCode(an.shape)); + Key(os, "ANISOTROPY_SHAPE_INTERCEPT", fmt::format("{:.3f}", an.shape_intercept)); + Key(os, "ANISOTROPY_SHAPE_INTERCEPT_Z", fmt::format("{:.1f}", an.shape_intercept_z)); + Key(os, "ANISOTROPY_SHAPE_SLOPE", fmt::format("{:.2f}", an.shape_slope)); + Key(os, "ANISOTROPY_SHAPE_RESIDUAL", fmt::format("{:.1f}", an.shape_residual)); + Key(os, "ANISOTROPY_SIGMA_SYSTEMATIC", fmt::format("{:.3f}", an.sigma_systematic)); + Key(os, "ANISOTROPY_FLOOR", fmt::format("{:.3f}", an.floor)); + Key(os, "ANISOTROPY_SIGNIFICANCE", fmt::format("{:.2f}", an.significance)); + Key(os, "ANISOTROPY_DETECTION_LIMIT", fmt::format("{:.2f}", an.detection_limit)); + os << "\n" << AnisotropyToText(an) << "\n"; + if (an.verdict == AnisotropyVerdict::Detected && an.d_min_spread > 0.5) + warnings.emplace_back(fmt::format( + "Diffraction is anisotropic (deltaB {:.1f} A^2; the diffraction limit runs from " + "{:.2f} to {:.2f} A depending on direction) - refinement and map interpretation " + "should allow for it; no intensity has been corrected for it here", + an.delta_b, *std::max_element(an.d_min_axis, an.d_min_axis + 3), + *std::min_element(an.d_min_axis, an.d_min_axis + 3))); + } + + // --------------------------------------------------------------- 10. WARNINGS if (result.cancelled) warnings.emplace_back(fmt::format("Processing was cancelled after {} images - this report " "describes an incomplete run", result.images_processed)); - Section(os, "9. WARNINGS"); + Section(os, "10. WARNINGS"); os << " Everything that needs a person's attention, one line each, marked so a script can find\n" << " them with a single grep for \"WARNING:\".\n\n"; Key(os, "WARNING_COUNT", warnings.size()); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 211949d5..af5aef6d 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -47,6 +47,7 @@ #include "../image_analysis/scale_merge/SearchSpaceGroup.h" #include "../image_analysis/geom_refinement/PostRefine.h" #include "../image_analysis/lattice_search/LatticeSearch.h" +#include "../image_analysis/scale_merge/AnisotropyAnalysis.h" #include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../image_analysis/scale_merge/HKLKey.h" #include "../image_analysis/scale_merge/ScaleOnTheFly.h" @@ -2815,6 +2816,27 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (std::isfinite(wilson.b) && wilson.b > 0.0) logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)", wilson.b, wilson.correlation, wilson.n_shells); + + // Diffraction anisotropy. Report-only, like the twinning and Wilson analyses above: it + // describes the data and corrects nothing. The tensor and its resolution signature come from + // the merged intensities; the error bar the verdict is gated on has to come from the + // unmerged observations, because a merge has exact Laue symmetry by construction and the + // tensor directions the symmetry forbids - the only place a dataset measures its own + // systematic error - are identically zero in it. + if (result.consensus_cell) { + AnisotropyRunInfo aniso_run; + if (sm.statistics.sweep_quality.measured && sm.statistics.sweep_quality.sweep_deg > 0.0f) + aniso_run.observed_rotation_deg = sm.statistics.sweep_quality.sweep_deg; + aniso_run.dose_term_in_scale_model = + experiment_.GetScalingSettings().GetCorrectionSurfaces(); + aniso_run.radiation_damage_relative_b = sm.statistics.radiation_damage_delta_b; + sm.statistics.anisotropy = AnalyzeAnisotropy( + sm.merged, + ScaledObservations(indexer->GetIntegrationOutcome(), + experiment_.IsRotationIndexing()), + *result.consensus_cell, twin_sg, aniso_run); + stats_text << AnisotropyToText(sm.statistics.anisotropy) << "\n"; + } } stats_text << sm.statistics; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index ef677a21..70bebdc3 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -36,6 +36,7 @@ #include "../image_analysis/scale_merge/StillsPartialityRefine.h" #include "../image_analysis/scale_merge/RotationScaleMerge.h" #include "../image_analysis/scale_merge/ResolutionCutoff.h" +#include "../image_analysis/scale_merge/AnisotropyAnalysis.h" #include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" #include "Rugnux.h" @@ -1552,6 +1553,21 @@ static int RunRugnux(int argc, char **argv) { const auto twinning = AnalyzeTwinning(merged_reflections, twin_sg); std::cout << std::endl << TwinningAnalysisToText(twinning) << std::endl; + // Diffraction anisotropy, exactly as the full pipeline computes it (Rugnux.cpp) - report-only, + // nothing is corrected. The stored reflections have just been re-scaled above, so they carry + // this merge's own per-image scale. + if (experiment.GetUnitCell()) { + AnisotropyRunInfo aniso_run; + if (merged_statistics.sweep_quality.measured && merged_statistics.sweep_quality.sweep_deg > 0.0f) + aniso_run.observed_rotation_deg = merged_statistics.sweep_quality.sweep_deg; + aniso_run.dose_term_in_scale_model = experiment.GetScalingSettings().GetCorrectionSurfaces(); + aniso_run.radiation_damage_relative_b = merged_statistics.radiation_damage_delta_b; + merged_statistics.anisotropy = AnalyzeAnisotropy(merged_reflections, + ScaledObservations(reflections, is_rotation), + *experiment.GetUnitCell(), twin_sg, aniso_run); + std::cout << AnisotropyToText(merged_statistics.anisotropy) << std::endl; + } + // Unmerged observations (--export-unmerged), from the integrated observations rather than the // merged ones: the partiality and the per-image scale are left for the reading program, which // fits a scale model of its own. diff --git a/tests/AnisotropyAnalysisTest.cpp b/tests/AnisotropyAnalysisTest.cpp new file mode 100644 index 00000000..8bd68d0f --- /dev/null +++ b/tests/AnisotropyAnalysisTest.cpp @@ -0,0 +1,153 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include + +#include "../image_analysis/scale_merge/AnisotropyAnalysis.h" +#include "gemmi/symmetry.hpp" +#include "gemmi/unitcell.hpp" + +namespace { + gemmi::UnitCell Cell(double a, double b, double c, double al, double be, double ga) { + gemmi::UnitCell out; + out.set(a, b, c, al, be, ga); + return out; + } + + // A synthetic merge: Wilson-distributed intensities with an isotropic fall-off and a known + // deviatoric anisotropy on top, over every hkl inside the resolution limit. The "structure + // factor" is a hash of the index, so the set is reproduced bit for bit. + std::vector SyntheticMerge(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg, + double d_min, double b_iso, + const gemmi::SMat33 &b_dev) { + const gemmi::GroupOps ops = sg->operations(); + const gemmi::ReciprocalAsu asu(sg); + std::vector out; + const int hmax = static_cast(cell.a / d_min) + 1; + const int kmax = static_cast(cell.b / d_min) + 1; + const int lmax = static_cast(cell.c / d_min) + 1; + for (int h = -hmax; h <= hmax; ++h) + for (int k = -kmax; k <= kmax; ++k) + for (int l = -lmax; l <= lmax; ++l) { + const gemmi::Op::Miller hkl{{h, k, l}}; + if ((h == 0 && k == 0 && l == 0) || !asu.is_in(hkl) || ops.is_systematically_absent(hkl)) + continue; + const double d2 = cell.calculate_1_d2(hkl); + if (d2 <= 0 || 1.0 / std::sqrt(d2) < d_min) + continue; + // Wilson draw from a hash of the index: deterministic, and spanning a realistic range. + const uint32_t seed = static_cast(h * 73856093 ^ k * 19349663 ^ l * 83492791); + const double u = ((seed * 2654435761u) >> 8) / static_cast(1 << 24); + const double wilson = -std::log(std::max(1e-6, u)); + // The tensor is applied in the Cartesian frame, s = F^T h. + const gemmi::Vec3 s = cell.frac.mat.left_multiply(gemmi::Vec3(h, k, l)); + const double aniso = -0.5 * b_dev.r_u_r(s); + MergedReflection r; + r.h = h; r.k = k; r.l = l; + r.d = static_cast(1.0 / std::sqrt(d2)); + r.I = static_cast(1000.0 * wilson * std::exp(-0.5 * b_iso * d2 + aniso)); + r.sigma = static_cast(0.02 * std::fabs(r.I) + 1.0); + out.push_back(r); + } + return out; + } +} + +// The number of free deviatoric anisotropy parameters is fixed by the Laue class alone. This is the +// self-test of the constraint basis: 5 / 3 / 2 / 1 / 1 / 0 for triclinic / monoclinic / orthorhombic / +// tetragonal / trigonal-hexagonal / cubic, and nothing else is possible. +TEST_CASE("Anisotropy free-parameter count", "[anisotropy]") { + struct Case { + const char *space_group; + gemmi::UnitCell cell; + int free_directions; + }; + const std::vector cases{ + {"P 1", Cell(51, 62, 73, 84.0, 95.0, 103.0), 5}, + {"P 1 21 1", Cell(51, 62, 73, 90.0, 95.0, 90.0), 3}, + {"C 1 2 1", Cell(91, 62, 73, 90.0, 105.0, 90.0), 3}, + {"P 21 21 21", Cell(51, 62, 73, 90.0, 90.0, 90.0), 2}, + {"I 2 2 2", Cell(51, 62, 73, 90.0, 90.0, 90.0), 2}, + {"P 43 21 2", Cell(79, 79, 38, 90.0, 90.0, 90.0), 1}, // lysozyme, the field's test specimen + {"P 31 2 1", Cell(62, 62, 91, 90.0, 90.0, 120.0), 1}, + {"R 3 :H", Cell(78, 78, 33, 90.0, 90.0, 120.0), 1}, + {"P 63", Cell(62, 62, 91, 90.0, 90.0, 120.0), 1}, + {"I 2 3", Cell(78, 78, 78, 90.0, 90.0, 90.0), 0}, + {"F 4 3 2", Cell(78, 78, 78, 90.0, 90.0, 90.0), 0}, + }; + for (const auto &c : cases) { + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(c.space_group); + REQUIRE(sg != nullptr); + std::vector merged(1); + merged[0].h = 1; merged[0].k = 0; merged[0].l = 0; + merged[0].I = 1.0f; merged[0].sigma = 1.0f; merged[0].d = 10.0f; + const auto result = AnalyzeAnisotropy(merged, {}, c.cell, sg); + INFO(c.space_group); + CHECK(result.n_free_parameters == c.free_directions); + } +} + +// A refined cell need not obey its space group's metric constraints exactly, and rugnux writes the +// unconstrained refined cell. An angle a hundredth of a degree off 90 must not create an extra +// anisotropy direction. +TEST_CASE("Anisotropy free-parameter count with an off-metric refined cell", "[anisotropy]") { + std::vector merged(1); + merged[0].h = 1; merged[0].k = 0; merged[0].l = 0; + merged[0].I = 1.0f; merged[0].sigma = 1.0f; merged[0].d = 10.0f; + CHECK(AnalyzeAnisotropy(merged, {}, Cell(51, 62, 73, 90.02, 105.0, 89.97), + gemmi::find_spacegroup_by_name("C 1 2 1")).n_free_parameters == 3); + CHECK(AnalyzeAnisotropy(merged, {}, Cell(51.0, 62.0, 73.0, 89.98, 90.03, 90.01), + gemmi::find_spacegroup_by_name("I 2 2 2")).n_free_parameters == 2); + CHECK(AnalyzeAnisotropy(merged, {}, Cell(78.01, 77.99, 78.02, 90.01, 89.99, 90.0), + gemmi::find_spacegroup_by_name("I 2 3")).n_free_parameters == 0); +} + +// A cubic crystal has no free deviatoric parameter, so its deltaB is exactly zero by symmetry - not +// small, not measured, zero - and the verdict is a statement about symmetry rather than about data. +TEST_CASE("Anisotropy is exactly zero in a cubic Laue class", "[anisotropy]") { + const gemmi::UnitCell cell = Cell(78, 78, 78, 90, 90, 90); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("I 2 3"); + const auto merged = SyntheticMerge(cell, sg, 2.5, 20.0, {8.0, -4.0, -4.0, 0.0, 0.0, 0.0}); + REQUIRE(merged.size() > 1000); + const auto result = AnalyzeAnisotropy(merged, {}, cell, sg); + CHECK(result.n_free_parameters == 0); + CHECK(result.delta_b == 0.0); + CHECK(result.verdict == AnisotropyVerdict::NotDetected); +} + +// The tensor itself: put a known deviatoric B into a tetragonal merge and read it back. The +// tetragonal Laue class leaves one free direction, along c*, and its magnitude is what deltaB means. +TEST_CASE("Anisotropy tensor is recovered from a synthetic merge", "[anisotropy]") { + const gemmi::UnitCell cell = Cell(79, 79, 38, 90, 90, 90); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 43 21 2"); + // Uniaxial about c, deltaB = B_zz - B_xx = 15 A^2. + const gemmi::SMat33 b_dev{-5.0, -5.0, 10.0, 0.0, 0.0, 0.0}; + const auto merged = SyntheticMerge(cell, sg, 2.0, 20.0, b_dev); + REQUIRE(merged.size() > 2000); + const auto result = AnalyzeAnisotropy(merged, {}, cell, sg); + REQUIRE(result.n_free_parameters == 1); + REQUIRE(result.n_cells > 0); + CHECK(result.delta_b == Catch::Approx(15.0).margin(1.5)); + // c* is the weak direction here, so the largest principal value points along z. + CHECK(std::fabs(result.eigenvector[0][2]) == Catch::Approx(1.0).margin(0.05)); + // A pure Debye-Waller fall-off is a straight line through the origin in s^2. + CHECK(result.shape == AnisotropyShape::Linear); + // With no unmerged observations the systematic-error scale cannot be measured, and the verdict + // says so rather than falling back on a counting-statistics error bar. + CHECK(result.verdict == AnisotropyVerdict::CannotDetermine); + CHECK_FALSE(result.refusal.empty()); +} + +// An isotropic merge must not produce a tensor, whatever the Laue class allows. +TEST_CASE("Anisotropy of an isotropic synthetic merge is small", "[anisotropy]") { + const gemmi::UnitCell cell = Cell(79, 79, 38, 90, 90, 90); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 43 21 2"); + const auto merged = SyntheticMerge(cell, sg, 2.0, 20.0, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); + REQUIRE(merged.size() > 2000); + const auto result = AnalyzeAnisotropy(merged, {}, cell, sg); + REQUIRE(result.n_cells > 0); + CHECK(std::fabs(result.delta_b) < 1.0); +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 9892c0b5..0757fd7c 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -89,6 +89,7 @@ ADD_EXECUTABLE(jfjoch_test SyntheticMergedReflections.h XDSPluginTest.cpp MergeScaleTest.cpp + AnisotropyAnalysisTest.cpp RfreeFlagsTest.cpp FrenchWilsonTest.cpp ReindexAmbiguityTest.cpp