rugnux: grade anisotropy on the ratio of the directional limits, not their spread in A

The small / noticeable / STRONG grade (and the ANISOTROPY warning) was read on the spread of
the directional diffraction limits in A, 0.25 and 0.5 A, so the same anisotropy graded
differently at every resolution: a 0.5 A spread is a 12% loss of resolution at 4 A and a
33% one at 1.5 A, and the battery had 4 A crystals losing 15% of their reach graded STRONG
beside 1 A crystals losing a third of it graded noticeable. The grade is now read on the
ratio of the weakest limit to the best (new key ANISOTROPY_D_MIN_RATIO; 1.125 and 1.25,
the former spreads read at 2 A), and so is the directional hypothesis of the verdict. The
headline prose quotes the factor.

ANISOTROPY_DELTA_B is unchanged. Its difference from phenix.xtriage's B_cart (a quarter
larger on a typical crystal, several-fold on one whose weakest direction stops diffracting
inside the measured range) was traced to xtriage's positivity cut on amplitudes and its
single Wilson curve over the written range - restricting the model-free fall-off of <I> in
cones to I > 0 reproduces xtriage's number, and xtriage's range grows as the data are
truncated at the weak direction's limit - while rugnux's tensor follows the model-free
fall-off; the report and the docs now say what to expect from amplitude-based tools. The
directional limits reproduce the ellipsoids of STARANISO-processed depositions of the same
crystals to about 0.1 A per axis. The log's strongest/weakest factor is capped at 1000x as
the report's key already was.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK

docs/CPU_DATA_ANALYSIS_DECISIONS.md

@@ -102,6 +102,10 @@ with one free constant per resolution shell, so every isotropic feature — the
 
 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.
 
+That is not a small effect, and it is the whole of the difference from phenix.xtriage's `B_cart`. Measured over the battery against the model-free fall-off of $\langle I\rangle$ (negatives kept) in 20° cones about the principal directions: the tensor above reproduces that fall-off; xtriage's eigenvalue range is a quarter lower on a typical crystal (median ratio 1.36), and four to eight times lower on a crystal whose weakest direction stops diffracting inside the measured range — its positivity cut floors the dead direction at the noise's positive half, and its single Wilson curve over 11.2 Å to $d_\min$ averages the deficit over the shells beyond that direction's limit, so the value also shrinks as the written range extends past it (truncating one such data set at the fitted resolution raised xtriage's range by a third; a shell earlier, it doubled). Restricting the same cone fall-off to positive intensities reproduces xtriage's number. xtriage's own verdict, range over mean $B \ge 0.5$ with mean $B > 20$ (the Read et al. 2011 cutoffs), is blind on low-resolution data for a second reason: the mean $B$ it divides by is 100–180 Ų there, so a crystal whose weakest direction reaches 5.6 Å against 3.7 Å along the best scores "not significantly anisotropic".
+
+**The grade** — *measurable but small*, *noticeable*, *STRONG* (the warning) — is read on the directional limits, as the ratio of the weakest to the best (1.125 and 1.25), not as their spread in Å: the same spread is a 12% loss of resolution at 4 Å and a 33% one at 1.5 Å, and a spread graded every 4 Å crystal strong and a 1 Å crystal losing a third of its reach noticeable. The two ratios are the former 0.25 Å and 0.5 Å spreads read at 2 Å, so the bands keep their meaning where most data sit. The limits are the number to compare between programs: on the open-arm depositions that were processed with STARANISO, rugnux's limits reproduce the deposited ellipsoid to about 0.1 Å per axis and its axis ratio to a few per cent.
+
 **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, read by interpolation in $s^2$ over equal-count shells — 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. Where $\langle I/\sigma(I)\rangle$ never falls through 2 in a direction, the limit returned is the **edge of the measured data** rather than the crystal's own; such a direction is marked — with a `<` in the report, a 1 in `ANISOTROPY_D_MIN_CENSORED`, and a note in the mmCIF — so the spread is not read as a measurement when it is a lower bound. Both constants are AIMLESS's defaults (cone half-angle 20°, $\langle I/\sigma\rangle = 2$). The 2 is a per-direction diagnostic level only: the dataset-wide resolution cut (§13.4) is CC$_{1/2}$-based, and the two criteria are not interchangeable.
 
 **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.

docs/RUGNUX_REPORT.md

@@ -1094,13 +1094,14 @@ spread in directional limit, or the reverse.
 | `ANISOTROPY_D_MIN_PRINCIPAL` | Diffraction limit (Å) along each principal direction — where ⟨I/σ(I)⟩ in a 20° cone about it falls through 2 |
 | `ANISOTROPY_D_MIN_CENSORED` | One flag per direction. `1` means ⟨I/σ(I)⟩ never fell through 2, so the limit is the **edge of the measured data**, a bound and not a measurement. The prose marks it with a `<` |
 | `ANISOTROPY_D_MIN_SPREAD` | Range of the three limits — itself a lower bound if any is censored |
+| `ANISOTROPY_D_MIN_RATIO` | The weakest limit over the best. The grade is read on it — small below 1.125, noticeable below 1.25, **STRONG** from 1.25 (the `WARNING:` line) — a ratio rather than a spread in Å, so that it means the same at every resolution |
 | `ANISOTROPY_SHAPE` | `LINEAR` (a real Debye–Waller *B*) \| `FLAT` (the deficit does not follow a *B* at all, so Δ*B* may be an **under**-estimate) \| `CONVEX` (grows faster than a *B* can) \| `UNDETERMINED` (the verdict moved on rebinning) |
 | `ANISOTROPY_FLOOR`, `ANISOTROPY_SIGNIFICANCE` | The data set's own systematic-error floor (Ų) and Δ*B*<sub>linear</sub> over it. Banded: below 2 not established, 2–3.5 marginal, above 3.5 established, above 5 strong |
 | `ANISOTROPY_DETECTION_LIMIT` | The smallest Δ*B* that could have been established on these data. It is set by systematic error, not by counting, so it does **not** improve with more reflections or a longer exposure |
 | `ANISOTROPY_N_OBSERVATIONS`, `ANISOTROPY_FORBIDDEN_Z`, `ANISOTROPY_SIGMA_SYSTEMATIC` | The unmerged observations the floor was measured on, that measurement against its own counting noise, and the floor before the counting part is added back |
 
-The default report carries `ANISOTROPY_VERDICT`, `ANISOTROPY_DELTA_B`, `ANISOTROPY_D_MIN_PRINCIPAL`
-and `ANISOTROPY_D_MIN_SPREAD`; the rest of this table — the detection gate's own parameters — is
+The default report carries `ANISOTROPY_VERDICT`, `ANISOTROPY_DELTA_B`, `ANISOTROPY_D_MIN_PRINCIPAL`,
+`ANISOTROPY_D_MIN_SPREAD` and `ANISOTROPY_D_MIN_RATIO`; the rest of this table — the detection gate's own parameters — is
 written with `--developer`.
 
 **`CANNOT_DETERMINE` is a real answer, not an evasion.** The verdict is not measured against counting
@@ -1113,8 +1114,19 @@ reason: a **triclinic** Laue class (no forbidden direction exists), an observed
 (`--no-scaling-corrections`), or no unmerged observations. A **cubic** Laue class is different again
 — symmetry forces Δ*B* to be exactly zero, and the run says that rather than reporting a measurement.
 
-Where anisotropy is detected and the directional limits differ by more than 0.5 Å, a `WARNING:` line
-says so, since refinement and map interpretation should allow for it.
+Where anisotropy is detected and the weakest directional limit is at least 1.25 times the best
+(`ANISOTROPY_D_MIN_RATIO`), a `WARNING:` line says so, since refinement and map interpretation
+should allow for it.
+
+**Comparing Δ*B* with other programs.** Expect phenix.xtriage (and any amplitude-based tool) to
+report a smaller anisotropic *B*: it converts to amplitudes and drops the non-positive intensities,
+and fits one Wilson curve over the whole range, which flattens the fall-off exactly where a direction
+is weak. Measured over the battery against the model-free fall-off of ⟨I⟩ in cones about the
+principal directions, rugnux's Δ*B* follows that fall-off, while xtriage's eigenvalue range sits a
+quarter lower on a typical crystal and four to eight times lower on one whose weakest direction stops
+diffracting inside the measured range; its value also shrinks the further the written data extend
+past that direction's limit. The directional limits are the comparison to make — they reproduce the
+ellipsoid of a STARANISO-processed deposition of the same crystal to within about 0.1 Å per axis.
 
 
 ## Model validation

image_analysis/scale_merge/AnisotropyAnalysis.cpp

@@ -1344,8 +1344,10 @@ AnisotropyResult AnalyzeAnisotropy(const std::vector<MergedReflection> &merged,
         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))
+        if (std::isfinite(lo) && std::isfinite(hi)) {
             result.d_min_spread = hi - lo;
+            result.d_min_ratio = hi / lo;
+        }
     }
 
     // --- the resolution signature, and whether the binning decides it ---
@@ -1493,14 +1495,14 @@ AnisotropyResult AnalyzeAnisotropy(const std::vector<MergedReflection> &merged,
     // 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
+    // difference must be established AND large enough to be the noticeable ratio 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;
+            && result.d_min_ratio >= ANISOTROPY_LIMIT_RATIO_NOTICEABLE;
     // 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)
@@ -1566,7 +1568,11 @@ std::string AnisotropyToText(const AnisotropyResult &result) {
     for (double v : result.eigenvalue)
         os << " " << (v - result.eigenvalue[2]);
     os << "\n";
-    if (std::isfinite(result.fold_weakening))
+    // Capped as the report caps it: past the weak direction's own limit the number is the Debye-Waller
+    // factor extrapolated into a region with no data, and seven figures of it state nothing physical.
+    if (result.fold_weakening > 1000.0)
+        os << "  Strongest / weakest direction at the resolution limit: > 1000x\n";
+    else if (std::isfinite(result.fold_weakening))
         os << "  Strongest / weakest direction at the resolution limit: " << result.fold_weakening << "x\n";
     // Every axis, not just the outer two: ConeLimits fills each independently, so a crystal whose
     // first and third directions fall in cones too sparse to measure - but whose second is well

image_analysis/scale_merge/AnisotropyAnalysis.h

@@ -58,8 +58,14 @@ 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;
+// The ratio of the weakest to the best directional limit (d_weak / d_best) the report calls noticeable
+// rather than small, and strong rather than noticeable. A ratio and not a spread in A, because the
+// same spread is a 12% loss of resolution at 4 A and a 33% one at 1.5 A, and a grade that read
+// "strong" on every 4 A crystal and "noticeable" on a 1 A one losing a third of its reach was
+// grading the resolution, not the anisotropy. 1.125 and 1.25 are the former 0.25 A and 0.5 A
+// spreads read at 2 A, where most of the battery sits, so the bands keep their meaning there.
+constexpr double ANISOTROPY_LIMIT_RATIO_NOTICEABLE = 1.125;
+constexpr double ANISOTROPY_LIMIT_RATIO_STRONG = 1.25;
 
 const char *AnisotropyShapeCode(AnisotropyShape shape);       // e.g. "LINEAR"
 const char *AnisotropyVerdictCode(AnisotropyVerdict verdict); // e.g. "NOT_DETECTED"
@@ -89,6 +95,7 @@ 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
+    double d_min_ratio = NAN;             // max / min of the above: what the bands above are read on
     // <I/sigma(I)> 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;
@@ -130,8 +137,8 @@ struct AnisotropyResult {
     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,
+    // cone is below the best one by more than 3.5 of their joint standard error, and the limits differ
+    // by at least ANISOTROPY_LIMIT_RATIO_NOTICEABLE. 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;

rugnux/ResultReport.cpp

@@ -2901,8 +2901,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                             ? fmt::format("{:.2f}", an.d_min_axis[i]) : std::string("-"));
                 Add(s, KeyText("ANISOTROPY_D_MIN_PRINCIPAL", principal));
             }
-            if (std::isfinite(an.d_min_spread))
+            if (std::isfinite(an.d_min_spread)) {
                 Add(s, KeyReal("ANISOTROPY_D_MIN_SPREAD", an.d_min_spread, "{:.2f}"));
+                Add(s, KeyReal("ANISOTROPY_D_MIN_RATIO", an.d_min_ratio, "{:.3f}"));
+            }
 
             // The gate's internals. Correct, and none of it changes what a user does next; the verdict
             // key above already carries the decision they were computed to make.
@@ -2958,18 +2960,19 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                        ? ReciprocalAxisLabel(*result.consensus_cell, an.eigenvector[n])
                        : fmt::format("principal direction {}", n + 1);
             };
-            // 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 < ANISOTROPY_SPREAD_NOTICEABLE_A ? "measurable but small"
-                                 : an.d_min_spread < 0.5  ? "noticeable"
-                                                          : "STRONG";
+            // Severity in the one unit a user feels - how far apart the diffraction limits are, as a
+            // ratio, so that it means the same at 4 A as at 1 A - 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_ratio < ANISOTROPY_LIMIT_RATIO_NOTICEABLE ? "measurable but small"
+                                 : an.d_min_ratio < ANISOTROPY_LIMIT_RATIO_STRONG     ? "noticeable"
+                                                                                      : "STRONG";
             Add(s, Blank());
             if (an.verdict == AnisotropyVerdict::Detected && worst && best && worst != best) {
                 std::string p = fmt::format(
                         " Anisotropy: {}. Diffraction reaches {:.2f} A along {} and {:.2f} A along {}\n"
-                        " (deltaB {:.1f} A^2).", severity, *best, along(best), *worst, along(worst),
-                        an.delta_b);
+                        " (a factor {:.2f} in resolution; deltaB {:.1f} A^2).", severity, *best, along(best),
+                        *worst, along(worst), an.d_min_ratio, an.delta_b);
                 // "Censored" is survival-analysis jargon for good news, and it read as an accusation.
                 // The whole content of the flag, in plain words:
                 const int worst_i = static_cast<int>(worst - an.d_min_axis);
@@ -2998,12 +3001,17 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                          " data set's own systematic error. ANISOTROPY_DELTA_B is the range of the principal\n"
                          " components of the anisotropy tensor, on the ordinary crystallographic B scale (the same\n"
                          " scale as phenix.xtriage's B_cart and ctruncate's anisotropic B), fitted on intensities\n"
-                         " with nothing dropped; ANISOTROPY_SIGNIFICANCE gates ANISOTROPY_DELTA_B_LINEAR, the part\n"
-                         " of it that follows exp(-1/2 s^T B s), which is not the same number. A 1 in\n"
-                         " ANISOTROPY_D_MIN_CENSORED marks a direction whose limit is the edge of the measured\n"
-                         " data rather than the crystal's own. ANISOTROPY_FOLD_WEAKENING is evaluated at the\n"
-                         " GLOBAL d_min, which extrapolates past the weak direction's own limit, and is capped\n"
-                         " at 1000 here for that reason.", true));
+                         " with nothing dropped. Expect those amplitude-based tools to return less: they drop the\n"
+                         " non-positive intensities and fit one Wilson curve over the whole range, which flattens\n"
+                         " the fall-off wherever a direction is weak - by a quarter on a typical crystal, several-\n"
+                         " fold on one whose weakest direction stops diffracting inside the measured range.\n"
+                         " ANISOTROPY_SIGNIFICANCE gates ANISOTROPY_DELTA_B_LINEAR, the part of the deltaB that\n"
+                         " follows exp(-1/2 s^T B s), which is not the same number. The grade (small, noticeable,\n"
+                         " STRONG) is read on ANISOTROPY_D_MIN_RATIO, the weakest directional limit over the best,\n"
+                         " so that it means the same at every resolution. A 1 in ANISOTROPY_D_MIN_CENSORED marks a\n"
+                         " direction whose limit is the edge of the measured data rather than the crystal's own.\n"
+                         " ANISOTROPY_FOLD_WEAKENING is evaluated at the GLOBAL d_min, which extrapolates past the\n"
+                         " weak direction's own limit, and is capped at 1000 here for that reason.", true));
             if (fitted_resolution && best_axis)
                 Add(s, Prose(fmt::format(
                         " Three resolutions, three questions. INCLUDE_RESOLUTION_RANGE ({:.2f} A) is what was\n"
@@ -3013,7 +3021,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                         " one without saying which it is, is.",
                         result.merge_statistics.overall.d_min, *fitted_resolution, *best_axis), true));
 
-            if (an.verdict == AnisotropyVerdict::Detected && an.d_min_spread > 0.5) {
+            if (an.verdict == AnisotropyVerdict::Detected && an.d_min_ratio >= ANISOTROPY_LIMIT_RATIO_STRONG) {
                 if (worst && best && worst != best)
                     Warn(doc, PathologyCode::ANISOTROPY, fmt::format(
                             "Diffraction is anisotropic (deltaB {:.1f} A^2; the diffraction limit is {:.2f} A "
@@ -3698,8 +3706,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                         checked_row("Anisotropy", {PathologyCode::ANISOTROPY},
                                     "not established above this data set's own noise");
                     else {
-                        const char *sev = an.d_min_spread < ANISOTROPY_SPREAD_NOTICEABLE_A ? "small"
-                                        : an.d_min_spread < 0.5 ? "noticeable" : "STRONG";
+                        const char *sev = an.d_min_ratio < ANISOTROPY_LIMIT_RATIO_NOTICEABLE ? "small"
+                                        : an.d_min_ratio < ANISOTROPY_LIMIT_RATIO_STRONG ? "noticeable" : "STRONG";
                         const char *by = an.tensor_detected && an.directional_detected
                                          ? "the tensor and the directional limits"
                                          : an.tensor_detected ? "the tensor" : "the directional limits";

tests/AnisotropyAnalysisTest.cpp

@@ -194,7 +194,7 @@ TEST_CASE("Anisotropy directional limits establish a difference the cones measur
     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.d_min_ratio > ANISOTROPY_LIMIT_RATIO_NOTICEABLE);
     CHECK(result.weak_cone_i_over_sigma < result.best_cone_i_over_sigma);
     CHECK(result.directional_detected);
     CHECK_FALSE(result.tensor_detected);