diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index da6c0905..cd5295cb 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -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 `_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 `_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/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index c4e23584..6e701d74 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -309,12 +309,15 @@ void WriteMmcifReflections(const std::vector &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) - << " # = 2 in a 20 deg cone about eigenvector " << (i + 1) << "\n"; + << " # = 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"; diff --git a/image_analysis/scale_merge/AnisotropyAnalysis.cpp b/image_analysis/scale_merge/AnisotropyAnalysis.cpp index 84f50b76..90a5b962 100644 --- a/image_analysis/scale_merge/AnisotropyAnalysis.cpp +++ b/image_analysis/scale_merge/AnisotropyAnalysis.cpp @@ -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 &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 ops = sg ? sg->operations().sym_ops : std::vector{gemmi::Op::identity()}; std::vector 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); + // 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(cells.entry_cell.size()); out.ok = true; return out; @@ -1155,7 +1173,7 @@ AnisotropyResult AnalyzeAnisotropy(const std::vector &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 &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 &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 ( = " << 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 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" diff --git a/image_analysis/scale_merge/AnisotropyAnalysis.h b/image_analysis/scale_merge/AnisotropyAnalysis.h index 7d650975..a340d0ee 100644 --- a/image_analysis/scale_merge/AnisotropyAnalysis.h +++ b/image_analysis/scale_merge/AnisotropyAnalysis.h @@ -80,6 +80,10 @@ struct AnisotropyResult { // --- 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}; + // True when 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 diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index c7589b36..0dd65bf1 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -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));