mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 2 Packages Projects Releases 71 Wiki Activity

Anisotropy: say what was measured, and stop crying wolf

Browse Source 
Four reporting changes and a restored changelog line. The gate boundaries are
unchanged - see below - and no merged number moves: .mtz and .hkl are
byte-identical on two crystals, and REPORT_VERSION stays 3 because everything
here is an added key or prose.

A directional diffraction limit that is really the edge of the measured data
was printed as if it were the crystal's. One crystal reported 1.09 / 0.99 /
0.99 A against a 0.994 A cutoff, so two of its three limits were the cutoff.
A censored limit now prints with a "<" and sets ANISOTROPY_D_MIN_CENSORED. It
fires on 4 of 28 crystals, and on two of those all three limits are censored,
which means their reported directional spread was binning noise.

The caution about a too-high symmetry assignment fired on every crystal with
one deviatoric degree of freedom - every tetragonal, trigonal and hexagonal
case, 17 of 38. A caution that fires on 45% of runs is noise. It is now
conditional on the signature that actually indicates it, a large
symmetry-forbidden signal with a gate that established nothing, and fires on
2 of 38: the crystal already on record as space-group unstable, and one
other. The forbidden-direction z is now reported so the reader can see it;
its battery median is 6.5.

The verdict line quoted deltaB_linear while the headline reported deltaB, and
the two differ - deltaB_linear is the larger on 12 of 31 crystals. It now
writes out the arithmetic, says which number the gate tests, and says which
to act on. The observation count the floor was measured on is reported, and
the mmCIF verdict item carries its vocabulary like its neighbours.

On the boundaries: they were calibrated against XDS-derived unmerged data,
while production measures the floor from rugnux's own observations, and those
disagree by 0.12x to 16.2x with 7 of 29 verdicts flipping. Re-measured on
production input - by merging each cubic crystal in proper subgroups of its
own Laue class, where the true anisotropy is exactly zero and every deltaB
that comes back is manufactured - the false-positive rate is 24% at 2.0, 10%
at 3.5 and 5% at 5.0. Each sits inside the published band read as its
worst-case figure rather than its average. 3.5 maximises power minus
false-positive rate across the grid; moving to 4.0 costs a real detection and
reduces the false-positive rate by nothing. So the boundaries stay, and the
"strong" band at 5.0 is a 5% statement on real data rather than the 0.5% its
derivation suggested.

The changelog entry for the diagnostic itself was committed with only its
first line, having been extracted by a script that took one line of a
four-line bullet. Restored.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P
This commit is contained in:
leonarski_f
2026-08-25 23:56:08 +02:00
co-authored by Claude Opus 5
parent 48008e1447
commit 0da593b32a
5 changed files with 101 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/CHANGELOG.md
+8
View File
@@ -10,7 +10,15 @@ This is an UNSTABLE release. It includes many experimental features, as well as
* rugnux: the space-group search prints the twin-law disagreement H for every operator it tested and the H ratio of the point group it adopted, on every run, instead of only when that ratio refuses a promotion.
* rugnux: `--mode scale` works on a `_process.h5` whose space group came from re-seating the lattice; a file written before this stops with a message naming the two cells instead of failing inside the merge.
* rugnux: on rotation data the integration signal radius is set from how wide the crystal's own spots are, measured in the pre-scan, instead of the fixed 4 px; `--adaptive-integration-radius=off` restores the fixed radius, and an explicit `--integration-radius` still overrides both.
* rugnux: a directional diffraction limit that is the edge of the measured data rather than the crystal's own limit is marked as such - with a `<` in the report and in `ANISOTROPY_D_MIN_CENSORED`, and in the mmCIF - so `ANISOTROPY_D_MIN_SPREAD` is not read as a measurement when it is a lower bound.
* rugnux: the anisotropy verdict line names which of `ANISOTROPY_DELTA_B` and `ANISOTROPY_DELTA_B_LINEAR` it is quoting, says which of the two to act on, and says why the second can be the larger.
* rugnux: the anisotropy caution about a too-high symmetry assignment now fires only where that is actually indicated - the symmetry-forbidden tensor directions far above their own counting noise together with a gate that established nothing - instead of on every tetragonal, trigonal and hexagonal data set.
* rugnux: the anisotropy section reports `ANISOTROPY_N_OBSERVATIONS`, the unmerged observations the systematic-error floor was measured on, and `ANISOTROPY_FORBIDDEN_Z`, that floor against its own counting noise.
* rugnux: the mmCIF item `_reflns.jfjoch_aniso_verdict` carries its vocabulary as a comment, as the other `jfjoch_aniso_` items already did.
* rugnux reports diffraction anisotropy: the anisotropic deltaB and the diffraction limit along each
principal direction, with a NOT DETECTED / DETECTED / CANNOT DETERMINE verdict against the data
set's own systematic error, in `<prefix>_report.txt` section 9 and in the mmCIF. No intensity is
corrected and no reflection is removed.
* rugnux: `--mode scale` reports the Wilson B-factor estimate instead of `WILSON_B= nan`.
* rugnux: `--export-unmerged` writes the integrated observations as `<prefix>_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.
image_analysis/WriteReflections.cpp
+5 -2
View File
@@ -309,12 +309,15 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
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";
out << "_reflns.jfjoch_aniso_verdict " << AnisotropyVerdictCode(an.verdict)
<< " # DETECTED / NOT_DETECTED / CANNOT_DETERMINE, at deltaB(linear)/floor > 3.5\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)
<< " # <I/sigma(I)> = 2 in a 20 deg cone about eigenvector " << (i + 1) << "\n";
<< " # <I/sigma(I)> = 2 in a 20 deg cone about eigenvector " << (i + 1)
<< (an.d_min_censored[i] ? "; at the edge of the measured data, so a bound" : "")
<< "\n";
}
out << "#\n";
image_analysis/scale_merge/AnisotropyAnalysis.cpp
+70 -16
View File
@@ -83,6 +83,14 @@ namespace {
constexpr double GATE_ESTABLISHED = 3.5;
constexpr double GATE_STRONG = 5.0;
// The forbidden subspace of a real dataset is dominated by systematics, not by counting noise, so
// this ratio is above 1 on essentially every dataset: measured over the battery its median is 6.5
// and only two of 22 exceed 20, so "it is above its counting null" would fire on almost all of
// them. This is the upper tail. Together with p == 1 and a gate that established nothing it picks
// out two datasets in the battery - one of them the crystal whose space-group determination is
// itself known to be unstable - where the unconditional p == 1 caution fired on 17 of 38.
constexpr double FORBIDDEN_Z_SUSPECT = 20.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;
@@ -757,7 +765,7 @@ namespace {
// Following Evans & Murshudov (2013) Acta Cryst. D69, 1204-1214
void ConeLimits(const std::vector<MergedReflection> &merged, const gemmi::UnitCell &cell,
const gemmi::SpaceGroup *sg, const Eigen3 &axes, int n_shells,
double (&d_min)[3]) {
double (&d_min)[3], bool (&censored)[3]) {
const std::vector<gemmi::Op> ops = sg ? sg->operations().sym_ops
: std::vector<gemmi::Op>{gemmi::Op::identity()};
std::vector<double> s2;
@@ -805,6 +813,7 @@ namespace {
}
for (int n = 0; n < 3; ++n) {
double last_s2 = NAN, last_isig = NAN, limit_s2 = NAN;
bool crossed = false;
for (int s = 0; s < n_shells; ++s) {
if (isig_n[n][s] < MIN_CELL_REFLECTIONS)
continue;
@@ -819,13 +828,19 @@ namespace {
// 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);
crossed = true;
break;
} else {
break;
}
}
if (std::isfinite(limit_s2) && limit_s2 > 0.0)
if (std::isfinite(limit_s2) && limit_s2 > 0.0) {
d_min[n] = 1.0 / std::sqrt(limit_s2);
// <I/sigma> never fell through the threshold: the outermost shell that carries this cone
// is still above it, so the value is the edge of the measured data and not this crystal's
// own limit along that direction.
censored[n] = !crossed;
}
}
}
@@ -944,6 +959,7 @@ namespace {
struct SystematicFloor {
bool ok = false;
double sigma_excess = NAN; // per tensor component, A^2
double forbidden_z = NAN; // that scatter over its own counting null
int n_observations = 0;
};
@@ -982,6 +998,8 @@ namespace {
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));
if (sigma_stat > 0.0)
out.forbidden_z = sigma_sys / sigma_stat;
out.n_observations = static_cast<int>(cells.entry_cell.size());
out.ok = true;
return out;
@@ -1155,7 +1173,7 @@ AnisotropyResult AnalyzeAnisotropy(const std::vector<MergedReflection> &merged,
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);
ConeLimits(merged, cell, space_group, axes, N_SHELLS, result.d_min_axis, result.d_min_censored);
{
double lo = INFINITY, hi = -INFINITY;
for (double d : result.d_min_axis)
@@ -1244,6 +1262,7 @@ AnisotropyResult AnalyzeAnisotropy(const std::vector<MergedReflection> &merged,
} else {
const SystematicFloor sys = MeasureSystematicFloor(unmerged, cell, space_group, full, forbidden);
result.n_observations = sys.n_observations;
result.forbidden_z = sys.forbidden_z;
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");
@@ -1282,11 +1301,20 @@ AnisotropyResult AnalyzeAnisotropy(const std::vector<MergedReflection> &merged,
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)
// R2's Laue-class caution, fired on its measured signature rather than on the Laue class alone.
// Too high a symmetry assignment pushes real anisotropy into the directions the systematic error is
// measured in, so sigma_sys is inflated by the signal itself and the gate collapses (measured: a
// 42x inflation of sigma_sys, a gate of 106 falling to 3.1). The signature is a forbidden subspace
// far above its own counting noise TOGETHER with a gate that has established nothing. p == 1 on its
// own is not the signature - it is true of every tetragonal, trigonal and hexagonal dataset, 17 of
// 38 in the battery, most of them correctly assigned.
if (p == 1 && result.verdict == AnisotropyVerdict::NotDetected
&& std::isfinite(result.forbidden_z) && result.forbidden_z > FORBIDDEN_Z_SUSPECT)
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");
"the space group is worth checking before this verdict is believed: the tensor directions "
"this Laue class forbids carry far more scatter than their own counting noise, yet nothing "
"was established in the single direction it allows - which is what a too-high symmetry "
"assignment looks like, real anisotropy being measured as systematic error");
return result;
}
@@ -1312,13 +1340,25 @@ std::string AnisotropyToText(const AnisotropyResult &result) {
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])) {
bool any_censored = false;
os << " d_min along the principal directions (<I/sigma(I)> = " << 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 << " -";
for (int n = 0; n < 3; ++n) {
if (!std::isfinite(result.d_min_axis[n])) {
os << " -";
} else if (result.d_min_censored[n]) {
os << " <" << result.d_min_axis[n];
any_censored = true;
} else {
os << " " << result.d_min_axis[n];
}
}
os << " A\n";
if (any_censored)
os << " \"<\" marks a direction in which <I/sigma(I)> had not fallen to "
<< CONE_I_OVER_SIGMA << " where the measured data\n"
<< " end: that figure is the edge of the data, not this crystal's limit along that\n"
<< " direction, and the spread between the directions is therefore a lower bound.\n";
}
if (result.shape_shells >= MIN_SHAPE_SHELLS) {
os << " Resolution signature of the deficit over " << result.shape_shells << " shells: "
@@ -1337,14 +1377,28 @@ std::string AnisotropyToText(const AnisotropyResult &result) {
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);
// floor on; on a crystal whose deficit is not a Debye-Waller fall-off the two differ. Which of
// the two numbers is which is said here rather than left to the reader: deltaB_linear exceeds
// the deltaB on roughly a third of a real battery, and a verdict line that quoted the larger
// one without saying so read as a contradiction of the headline above.
const bool linear = std::isfinite(result.delta_b_linear);
const double judged = std::max(0.0, 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";
<< result.band << "): ratio " << result.significance << " = "
<< (linear ? "deltaB_linear " : "deltaB ") << judged
<< " A^2 / this dataset's own systematic-error floor of " << result.floor << " A^2.\n";
if (linear) {
os << " The gate is applied to deltaB_linear, the part of the deltaB that follows\n"
<< " exp(-1/2 s^T B s)";
if (result.delta_b_linear > result.delta_b)
os << " - it exceeds the deltaB here because the fitted line passes below\n"
<< " the origin";
else if (result.delta_b_linear < 0.0)
os << " - the fitted linear slope is negative here, which is no established\n"
<< " Debye-Waller component at all rather than a small one, so it enters as zero";
os << ". The anisotropy to quote is the deltaB of " << result.delta_b << " A^2 above.\n";
}
}
if (std::isfinite(result.detection_limit))
os << " Below about " << result.detection_limit << " A^2 nothing could be established on this\n"
image_analysis/scale_merge/AnisotropyAnalysis.h
+9
View File
@@ -80,6 +80,10 @@ struct AnisotropyResult {
// --- the diffraction limits (model-free) --------------------------------------------------
// Highest resolution at which <I/sigma(I)> in a 20 deg cone about principal axis n is still 2.0.
double d_min_axis[3] = {NAN, NAN, NAN};
// True when <I/sigma(I)> in that cone had not yet fallen to 2.0 at the edge of the measured data.
// The value returned is then where the DATA stop, not where the signal does: it is a bound on the
// limit, not a measurement of it, and d_min_spread built from it is a lower bound.
bool d_min_censored[3] = {false, false, false};
double d_min_spread = NAN; // max - min of the above
// --- the resolution signature of the deficit ----------------------------------------------
@@ -106,6 +110,11 @@ struct AnisotropyResult {
// plus counting noise, would manufacture in the symmetry-ALLOWED directions; significance is
// delta_b_linear / floor.
double sigma_systematic = NAN;
// The scatter measured in the symmetry-FORBIDDEN directions over its own counting null. It is
// normally well above 1 (a real dataset's forbidden subspace is dominated by systematics, not by
// counting noise); it is large together with a gate that establishes nothing when a too-high Laue
// class has pushed real anisotropy into those directions - R2's "space group suspect" signature.
double forbidden_z = NAN;
double floor = NAN;
double significance = NAN;
double detection_limit = NAN; // the smallest delta_b this dataset could establish, A^2
rugnux/ResultReport.cpp
+9 -1
View File
@@ -290,7 +290,10 @@ std::string RenderResultReport(const std::string &output_prefix,
<< " 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";
<< " ctruncate's anisotropic B), fitted on intensities with nothing dropped;\n"
<< " ANISOTROPY_SIGNIFICANCE gates ANISOTROPY_DELTA_B_LINEAR, the part of it that follows\n"
<< " exp(-1/2 s^T B s), which is not the same number. A 1 in ANISOTROPY_D_MIN_CENSORED marks\n"
<< " a direction whose limit is the edge of the measured data rather than the crystal's own.\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));
@@ -301,13 +304,18 @@ std::string RenderResultReport(const std::string &output_prefix,
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_CENSORED", fmt::format("{} {} {}", an.d_min_censored[0] ? 1 : 0,
an.d_min_censored[1] ? 1 : 0,
an.d_min_censored[2] ? 1 : 0));
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_N_OBSERVATIONS", an.n_observations);
Key(os, "ANISOTROPY_SIGMA_SYSTEMATIC", fmt::format("{:.3f}", an.sigma_systematic));
Key(os, "ANISOTROPY_FORBIDDEN_Z", fmt::format("{:.1f}", an.forbidden_z));
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));