From fe181ca04232ccc982035decc142cbc55ce250d4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 11 Oct 2026 02:21:32 +0200 Subject: [PATCH] Anisotropy verdict: the directional limits are a second hypothesis beside the tensor The verdict gated only deltaB_linear, the Debye-Waller part of the tensor, against the data set's own systematic floor. A deficit that does not follow s^2 (a FLAT signature) can stop the signal far earlier in one direction and still read NOT_DETECTED: one open-arm set reported directional limits of 3.18/3.18/2.14 A beside NOT_DETECTED. The cone limits ( = 2 in 20 deg cones) were computed but never consulted. Now each statement is decided on its own evidence and either establishes anisotropy: - tensor: unchanged gate; - directional limits: at the best direction's limit the weakest cone's is below the best cone's by more than 3.5 of their joint standard error (both read linearly in s^2 between shell centres, as the limits are), the weakest limit is measured (not censored at the data edge), and the limits are at least 0.25 A apart - the spread the report already calls "noticeable". Why the magnitude term: with the standard errors alone (the counting statistics of the cones) the weakest cone is "significantly" below the best on 172 of 216 battery sets, down to limits 0.5% apart (measured with a replica of ConeLimits on the a897b08 battery's merged data). The cones differ systematically (multiplicity and coverage depend on direction), so significance alone made DETECTED near-universal. A refusal (CANNOT_DETERMINE: dose ramp, short sweep, triclinic, no unmerged data) is not overruled - its reasons bias the cones too - but now reports the measured directional limits. Report: ANISOTROPY_DETECTED_BY (TENSOR/DIRECTIONAL/BOTH/NONE); developer keys for the two cones' and errors; the summary row names the test that carried the verdict, and CANNOT_DETERMINE rows give the directional limits; the mmCIF verdict comment says so. Report only: no intensity, flag or decision downstream reads the verdict (the French-Wilson prior uses the tensor regardless). Measured (battery replica, 225 sets with limits): directional statement established on 95; 23 NOT_DETECTED -> DETECTED (directional ratio 1.14-1.59); of the 17 sets with ratio > 1.2 and an uncensored weakest limit, 11 flip, 5 stay CANNOT_DETERMINE with limits reported, 1 (sparse, noisy cones) stays NOT_DETECTED. ANISOTROPY warnings (spread > 0.5 A) appear on 6 more sets. Confirmed with rugnux runs on the ci tier, those 17 and the private arm: p.hkl and MTZ data identical to the base (e34a5fb88) on every set compared. Co-Authored-By: Claude Opus 5.5 (1M context) --- image_analysis/WriteReflections.cpp | 2 +- .../scale_merge/AnisotropyAnalysis.cpp | 80 +++++++++++++++++-- .../scale_merge/AnisotropyAnalysis.h | 16 ++++ rugnux/ResultReport.cpp | 50 ++++++++---- tests/AnisotropyAnalysisTest.cpp | 23 ++++++ 5 files changed, 149 insertions(+), 22 deletions(-) diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index fba15332c..e1c1bce98 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -474,7 +474,7 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_reflns.jfjoch_aniso_significance " << Fmt(an.significance, 2) << " # deltaB(linear) / floor\n"; out << "_reflns.jfjoch_aniso_verdict " << AnisotropyVerdictCode(an.verdict) - << " # DETECTED / NOT_DETECTED / CANNOT_DETERMINE, at deltaB(linear)/floor > 3.5\n"; + << " # DETECTED / NOT_DETECTED / CANNOT_DETERMINE: deltaB(linear)/floor > 3.5, or the directional limits\n"; for (int i = 0; i < 3; ++i) if (std::isfinite(an.d_min_axis[i])) out << "_reflns.jfjoch_aniso_d_min_" << (i + 1) << " " diff --git a/image_analysis/scale_merge/AnisotropyAnalysis.cpp b/image_analysis/scale_merge/AnisotropyAnalysis.cpp index ce8546b5f..8a983222b 100644 --- a/image_analysis/scale_merge/AnisotropyAnalysis.cpp +++ b/image_analysis/scale_merge/AnisotropyAnalysis.cpp @@ -787,7 +787,8 @@ 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], bool (&censored)[3]) { + double (&d_min)[3], bool (&censored)[3], double &weak_isig, double &weak_isig_se, + double &best_isig, double &best_isig_se) { const std::vector ops = sg ? sg->operations().sym_ops : std::vector{gemmi::Op::identity()}; std::vector s2; @@ -805,6 +806,7 @@ namespace { 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_sum2(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) { @@ -830,6 +832,7 @@ namespace { for (int n = 0; n < 3; ++n) if (in_cone[n]) { isig_sum[n][shell] += r.I / r.sigma; + isig_sum2[n][shell] += (r.I / r.sigma) * (r.I / r.sigma); isig_n[n][shell] += 1; } } @@ -864,6 +867,44 @@ namespace { censored[n] = !crossed; } } + + // The directional-limit statement: the weakest and the best direction's cones, both read at the + // best direction's limit (linear in s^2 between shell centres, as the limits themselves are). It + // needs the weakest limit to be measured - a censored one is the edge of the data, and the + // spread built on it only a bound. + int best = -1, weak = -1; + for (int n = 0; n < 3; ++n) { + if (!std::isfinite(d_min[n])) + continue; + if (best < 0 || d_min[n] < d_min[best]) best = n; + if (weak < 0 || d_min[n] > d_min[weak]) weak = n; + } + if (best < 0 || weak == best || censored[weak]) + return; + const double best_s2 = 1.0 / (d_min[best] * d_min[best]); + // of cone n at s^2 and its standard error. + auto at = [&](int n, double &mean, double &se) { + double prev_x = NAN, prev_m = NAN, prev_e = NAN; + for (int sh = 0; sh < n_shells; ++sh) { + const int k = isig_n[n][sh]; + if (k < MIN_CELL_REFLECTIONS) + continue; + const double x = 0.5 * (edge[sh] + edge[sh + 1]); + const double m = isig_sum[n][sh] / k; + const double e = std::sqrt(std::max(0.0, (isig_sum2[n][sh] - k * m * m) / (k - 1)) / k); + if (x >= best_s2) { + const double f = std::isfinite(prev_x) ? (best_s2 - prev_x) / (x - prev_x) : 1.0; + mean = std::isfinite(prev_x) ? prev_m + f * (m - prev_m) : m; + se = std::isfinite(prev_x) ? prev_e + f * (e - prev_e) : e; + return; + } + prev_x = x; prev_m = m; prev_e = e; + } + mean = prev_m; + se = prev_e; + }; + at(weak, weak_isig, weak_isig_se); + at(best, best_isig, best_isig_se); } const char *Band(double g) { @@ -1267,7 +1308,9 @@ 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, result.d_min_censored); + ConeLimits(merged, cell, space_group, axes, N_SHELLS, result.d_min_axis, result.d_min_censored, + result.weak_cone_i_over_sigma, result.weak_cone_i_over_sigma_se, + result.best_cone_i_over_sigma, result.best_cone_i_over_sigma_se); { double lo = INFINITY, hi = -INFINITY; for (double d : result.d_min_axis) @@ -1412,6 +1455,27 @@ AnisotropyResult AnalyzeAnisotropy(const std::vector &merged, "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"); + + // --- the second hypothesis: the directional limits --- + // The tensor gate judges the part of the fall-off that follows a Debye-Waller B. A deficit that + // does not - flat in s^2, say - can stop the signal far earlier in one direction and still pass + // under it. The cones measure where the signal stops with no model at all, so either statement + // establishes anisotropy on its own evidence. The cones' own standard errors are counting + // statistics, and on the open battery the weakest cone sits below the best one by more than 3.5 of + // them on four sets in five, down to limits 0.5% apart: the cones differ systematically + // (multiplicity and coverage depend on direction too) far beyond their counting noise. So the + // difference must be established AND large enough to be the noticeable spread the report already + // grades by, rather than significant alone. + result.tensor_detected = result.verdict == AnisotropyVerdict::Detected; + result.directional_detected = std::isfinite(result.weak_cone_i_over_sigma) + && std::isfinite(result.best_cone_i_over_sigma) + && result.best_cone_i_over_sigma - result.weak_cone_i_over_sigma + > GATE_ESTABLISHED * std::hypot(result.weak_cone_i_over_sigma_se, result.best_cone_i_over_sigma_se) + && result.d_min_spread >= ANISOTROPY_SPREAD_NOTICEABLE_A; + // Only a tensor that could be tested and was not established is overruled: where the tensor was + // refused, the reasons (a dose ramp, a short sweep) bias the cones just as much. + if (result.directional_detected && result.verdict == AnisotropyVerdict::NotDetected) + result.verdict = AnisotropyVerdict::Detected; return result; } @@ -1486,8 +1550,8 @@ std::string AnisotropyToText(const AnisotropyResult &result) { 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"; + if (!result.refusal.empty()) { + os << " => tensor: 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. Which of @@ -1497,7 +1561,7 @@ std::string AnisotropyToText(const AnisotropyResult &result) { 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.tensor_detected ? " => tensor: DETECTED (" : " => tensor: NOT DETECTED (") << 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"; @@ -1519,6 +1583,12 @@ std::string AnisotropyToText(const AnisotropyResult &result) { << " improve with more reflections or a longer exposure.\n"; for (const auto &c : result.cautions) os << " Note: " << c << "\n"; + if (std::isfinite(result.weak_cone_i_over_sigma)) + os << " => directional limits: " << (result.directional_detected ? "DETECTED" : "NOT DETECTED") + << ": at the best direction's limit is " << result.best_cone_i_over_sigma << " +/- " + << result.best_cone_i_over_sigma_se << " along it and " << result.weak_cone_i_over_sigma << " +/- " + << result.weak_cone_i_over_sigma_se << " along the weakest direction.\n"; + os << " Verdict: " << AnisotropyVerdictCode(result.verdict) << "\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 index 47021eb69..021eda984 100644 --- a/image_analysis/scale_merge/AnisotropyAnalysis.h +++ b/image_analysis/scale_merge/AnisotropyAnalysis.h @@ -58,6 +58,9 @@ enum class AnisotropyVerdict { Detected // anisotropy is established; `band` grades it }; +// The spread of the directional limits the report calls noticeable rather than small, in A. +constexpr double ANISOTROPY_SPREAD_NOTICEABLE_A = 0.25; + const char *AnisotropyShapeCode(AnisotropyShape shape); // e.g. "LINEAR" const char *AnisotropyVerdictCode(AnisotropyVerdict verdict); // e.g. "NOT_DETECTED" @@ -86,6 +89,12 @@ struct AnisotropyResult { // 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 + // in the weakest and the best direction's cones at the best direction's limit, and + // their standard errors. NaN where the weakest limit is censored or no two limits differ. + double weak_cone_i_over_sigma = NAN; + double weak_cone_i_over_sigma_se = NAN; + double best_cone_i_over_sigma = NAN; + double best_cone_i_over_sigma_se = NAN; // --- 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 @@ -119,6 +128,13 @@ struct AnisotropyResult { double floor = NAN; double significance = NAN; double detection_limit = NAN; // the smallest delta_b this dataset could establish, A^2 + // Two hypotheses, each decided on its own evidence, and the verdict is DETECTED when either holds: + // the tensor gate above, and the directional limits - at the best direction's limit the weakest + // cone is below the best one by more than 3.5 of their joint standard error, and the limits are + // at least ANISOTROPY_SPREAD_NOTICEABLE_A apart. The directional limits overrule NOT_DETECTED, + // never a refusal (CANNOT_DETERMINE), whose reasons bias the cones as much as the tensor. + bool tensor_detected = false; + bool directional_detected = false; AnisotropyVerdict verdict = AnisotropyVerdict::CannotDetermine; std::string band; // "not established" / "marginal" / "established" / "strong" std::string refusal; // why, when the verdict is CannotDetermine diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 55ab95b5a..db38ff87e 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -2170,6 +2170,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, if (merged && an.n_reflections > 0) { Add(s, KeyEnum("ANISOTROPY_VERDICT", AnisotropyVerdictCode(an.verdict), {"DETECTED", "NOT_DETECTED", "CANNOT_DETERMINE"})); + Add(s, KeyEnum("ANISOTROPY_DETECTED_BY", + an.tensor_detected && an.directional_detected ? "BOTH" + : an.tensor_detected ? "TENSOR" : an.directional_detected ? "DIRECTIONAL" : "NONE", + {"TENSOR", "DIRECTIONAL", "BOTH", "NONE"})); Add(s, KeyReal("ANISOTROPY_DELTA_B", an.delta_b, "{:.2f}")); // The one rule again: a direction the fit could not limit prints "-", and a triplet // with no measured direction at all writes no key - never the word "nan". @@ -2220,6 +2224,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyReal("ANISOTROPY_FLOOR", an.floor, "{:.3f}", true)); Add(s, KeyReal("ANISOTROPY_SIGNIFICANCE", an.significance, "{:.2f}", true)); Add(s, KeyReal("ANISOTROPY_DETECTION_LIMIT", an.detection_limit, "{:.2f}", true)); + Add(s, KeyReal("ANISOTROPY_WEAK_CONE_I_SIGMA", an.weak_cone_i_over_sigma, "{:.2f}", true)); + Add(s, KeyReal("ANISOTROPY_WEAK_CONE_I_SIGMA_SE", an.weak_cone_i_over_sigma_se, "{:.2f}", true)); + Add(s, KeyReal("ANISOTROPY_BEST_CONE_I_SIGMA", an.best_cone_i_over_sigma, "{:.2f}", true)); + Add(s, KeyReal("ANISOTROPY_BEST_CONE_I_SIGMA_SE", an.best_cone_i_over_sigma_se, "{:.2f}", true)); const double *worst = nullptr, *best = nullptr; for (const double *p = an.d_min_axis; p != an.d_min_axis + 3; ++p) { @@ -2237,7 +2245,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, // Severity in the one unit a user feels - how far apart the diffraction limits are - and // never in the confidence of the detection. "DETECTED (strong)" meant statistically // established and was read as severely anisotropic on every second dataset in the corpus. - const char *severity = an.d_min_spread < 0.25 ? "measurable but small" + const char *severity = an.d_min_spread < ANISOTROPY_SPREAD_NOTICEABLE_A ? "measurable but small" : an.d_min_spread < 0.5 ? "noticeable" : "STRONG"; Add(s, Blank()); @@ -2261,8 +2269,11 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " removed; only the French-Wilson amplitudes (F) take the direction into account.", severity, an.delta_b))); } else if (an.verdict == AnisotropyVerdict::CannotDetermine) { - Add(s, Prose(" Anisotropy: cannot be determined in this Laue class - there is no\n" - " symmetry-forbidden direction to calibrate the fit's own noise against.")); + std::string p = " Anisotropy: the tensor cannot be tested - " + an.refusal + "."; + if (worst && best && worst != best) + p += fmt::format("\n Diffraction reaches {:.2f} A along {} and {:.2f} A along {}.", + *best, along(best), *worst, along(worst)); + Add(s, Prose(p)); } else { Add(s, Prose(" Anisotropy: below this data set's own noise floor. No directional statement can\n" " be made, and no intensity was corrected for direction.")); @@ -2739,23 +2750,30 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, if (an.n_reflections == 0) untested("Anisotropy", {PathologyCode::ANISOTROPY}); else { - if (an.verdict != AnisotropyVerdict::Detected) + const double *lo = nullptr, *hi = nullptr; + for (const double *p = an.d_min_axis; p != an.d_min_axis + 3; ++p) { + if (!std::isfinite(*p)) + continue; + if (!lo || *p < *lo) lo = p; + if (!hi || *p > *hi) hi = p; + } + const std::string limits = lo && hi && lo != hi + ? fmt::format(", d_min {:.2f}-{:.2f} A by direction", *lo, *hi) : std::string(); + if (an.verdict == AnisotropyVerdict::CannotDetermine) + checked_row("Anisotropy", {PathologyCode::ANISOTROPY}, + "the tensor cannot be tested on this data set" + limits); + else if (an.verdict != AnisotropyVerdict::Detected) checked_row("Anisotropy", {PathologyCode::ANISOTROPY}, "not established above this data set's own noise"); else { - const double *lo = nullptr, *hi = nullptr; - for (const double *p = an.d_min_axis; p != an.d_min_axis + 3; ++p) { - if (!std::isfinite(*p)) - continue; - if (!lo || *p < *lo) lo = p; - if (!hi || *p > *hi) hi = p; - } - const char *sev = an.d_min_spread < 0.25 ? "small" + const char *sev = an.d_min_spread < ANISOTROPY_SPREAD_NOTICEABLE_A ? "small" : an.d_min_spread < 0.5 ? "noticeable" : "STRONG"; - checked_row("Anisotropy", {PathologyCode::ANISOTROPY}, lo && hi - ? fmt::format("{}: deltaB {:.1f} A^2, d_min {:.2f}-{:.2f} A by direction", - sev, an.delta_b, *lo, *hi) - : fmt::format("{}: deltaB {:.1f} A^2", sev, an.delta_b)); + const char *by = an.tensor_detected && an.directional_detected + ? "the tensor and the directional limits" + : an.tensor_detected ? "the tensor" : "the directional limits"; + checked_row("Anisotropy", {PathologyCode::ANISOTROPY}, + fmt::format("{}: deltaB {:.1f} A^2{} (established by {})", + sev, an.delta_b, limits, by)); } } if (result.tncs.measurable) { diff --git a/tests/AnisotropyAnalysisTest.cpp b/tests/AnisotropyAnalysisTest.cpp index fef49078c..400ea81ee 100644 --- a/tests/AnisotropyAnalysisTest.cpp +++ b/tests/AnisotropyAnalysisTest.cpp @@ -157,6 +157,29 @@ TEST_CASE("Anisotropy of an isotropic synthetic merge is small", "[anisotropy]") CHECK(std::fabs(result.delta_b) < 1.0); } +// The directional limits are a second hypothesis beside the tensor: the signal stops earlier along c* +// than along a* and b*, and the cones see it with no model of the fall-off. Without unmerged +// observations the tensor is refused, and a refusal stands - the cones overrule only NOT_DETECTED - +// but the directional statement is still made and reported. +TEST_CASE("Anisotropy directional limits establish a difference the cones measure", "[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"); + // B 70 A^2 along a*/b* and 100 A^2 along c*: = 2 near 2.4 A and 2.8 A respectively. + const auto merged = SyntheticMerge(cell, sg, 2.0, 80.0, {-10.0, -10.0, 20.0, 0.0, 0.0, 0.0}); + const auto result = AnalyzeAnisotropy(merged, {}, cell, sg); + REQUIRE(result.n_cells > 0); + REQUIRE(std::isfinite(result.weak_cone_i_over_sigma)); + CHECK(result.d_min_spread > ANISOTROPY_SPREAD_NOTICEABLE_A); + CHECK(result.weak_cone_i_over_sigma < result.best_cone_i_over_sigma); + CHECK(result.directional_detected); + CHECK_FALSE(result.tensor_detected); + CHECK(result.verdict == AnisotropyVerdict::CannotDetermine); + + // The same fall-off in every direction: the limits agree and nothing is established. + const auto isotropic = AnalyzeAnisotropy(SyntheticMerge(cell, sg, 2.0, 80.0, {0, 0, 0, 0, 0, 0}), {}, cell, sg); + CHECK_FALSE(isotropic.directional_detected); +} + namespace { // One partial of reflection (h,k,l) on image `frame`, with the fields ScaledObservations reads. Reflection Partial(int h, int k, int l, float frame, float I, float partiality) {