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

Anisotropy verdict: the directional limits are a second hypothesis beside the tensor

Browse Source 
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 (<I/sigma> = 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 <I/sigma> 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'
<I/sigma> 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) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-10-11 02:21:32 +02:00
co-authored by Claude Opus 5.5
parent e34a5fb88f
commit fe181ca042
5 changed files with 149 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/WriteReflections.cpp
+1 -1
View File
@@ -474,7 +474,7 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &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) << " "
image_analysis/scale_merge/AnisotropyAnalysis.cpp
+75 -5
View File
@@ -787,7 +787,8 @@ 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], 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<gemmi::Op> ops = sg ? sg->operations().sym_ops
: std::vector<gemmi::Op>{gemmi::Op::identity()};
std::vector<double> s2;
@@ -805,6 +806,7 @@ namespace {
const double cos_cone = std::cos(CONE_HALF_ANGLE_DEG * gemmi::pi() / 180.0);
std::vector<std::vector<double>> isig_sum(3, std::vector<double>(n_shells, 0.0));
std::vector<std::vector<double>> isig_sum2(3, std::vector<double>(n_shells, 0.0));
std::vector<std::vector<int>> isig_n(3, std::vector<int>(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]);
// <I/sigma(I)> 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<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, 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<MergedReflection> &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 <I/sigma(I)> 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();
}
image_analysis/scale_merge/AnisotropyAnalysis.h
+16
View File
@@ -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
// <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;
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
rugnux/ResultReport.cpp
+34 -16
View File
@@ -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) {
tests/AnisotropyAnalysisTest.cpp
+23
View File
@@ -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*: <I/sigma> = 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) {