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

Space-group search: read the presence cut on counting significance

Browse Source 
present_i_over_sigma is meant to say "this reflection is really there".
The merged I/sigma is not that quantity. The merge carries the error
model's intensity-proportional term, sigma^2 = a*sigma0^2 + (b*I)^2, so
I/sigma saturates - at ISa*sqrt(n) for a reflection observed n times, and
at ISa exactly for one observed once. Above that knee it stops rising
with the intensity.

Measured over the battery's 76 search merges: the median I/sigma of the
top E^2 decile sits at 0.81 of ISa, and on the weakest merge the decile
profile reads 1.58 1.61 1.63 1.64 1.64 1.65 1.65 against an ISa of 1.70 -
flat to three digits, log-log slope against E^2 of 0.008. Reflections an
order of magnitude apart in real intensity read the same number. So one
constant demands between 2.16 and 11.91 in counting significance
depending on the crystal, and the fraction of a merge clearing it runs
from 3.1% to 98.8%.

Convert it once instead. For a reflection observed once
sigma_counting^2 = sigma^2 - (b*I)^2, so I/sigma_counting >= T is exactly
t >= T/sqrt(1+(T/ISa)^2). That cut lies strictly below ISa for every ISa,
so it is always reachable by the reflection the ceiling binds hardest,
and it is within 1% of T on any merge with ISa >= 21, so a healthy merge
is left where it was. Multiplicity is taken as 1 deliberately rather than
estimated, for the same reason.

The quantile fallback goes with it. Its comment claimed to be inert on
the battery; it fires on 4 of 76 arms, and on one of them it relaxes the
cut BELOW the nominal value because that merge is bimodal. A cap on a
statistic is not a substitute for measuring the right statistic.

Two comments in the tree still described a merged sigma floored at b|I|,
removed in 72efb75a8 before rc-161. Both corrected; that stale claim is
what made this look like a different bug.

Battery, both arms, against a control floor of exactly zero: space groups
35/38 before and after, the same three crystals differing from the
reference, all 380 shells identical and all 38 merged mmCIFs
byte-identical. It is live rather than dead code - the pass_cc population
moves on 45 of 47 arms and 212 of 314 operator correlations move, with 0
of 314 verdicts changed. The weakest crystal's single 2-fold goes from CC
0.489 on 1131 pairs to 0.597 on 4822, widening its margin over the
threshold by 57%.

Two things this does not fix, both quantified in the work behind it. The
conversion is monotone, so it corrects the threshold but cannot restore
the ranking the (b*I)^2 term already destroyed - that needs a
counting-only merged sigma and touches the CUDA merge kernel. And one
crystal's gate still admits 98.8% of its merge with a slope of 0.037,
discriminating nothing; the instrument for that is E^2, which needs
min_operator_cc re-derived first.

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-26 12:43:04 +02:00
co-authored by Claude Opus 5
parent 278f8e3d73
commit e332aaf28c
4 changed files with 52 additions and 35 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/CHANGELOG.md
+1
View File
@@ -16,6 +16,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as
* rugnux: the spot-width pre-scan stops once the integration radius it is measuring has settled, instead of always working through the whole sample; the radius it chooses is unchanged.
* rugnux: on stills the integration background ring runs to 14 px instead of 12, so fewer reflections are discarded for want of a background ring and per-shell R_meas improves over most of the signal-bearing range; the gain reverses in shells below <I/sigma> of about 4.
* rugnux: the anisotropy diagnostic no longer spends minutes sorting at the end of a large rotation run; the run finishes around 10% sooner and reports the same numbers.
* rugnux: the space-group search decides whether a reflection is genuinely present from its counting significance rather than from the merged I/sigma, which saturates at the merge's own ISa; the cut now means the same thing on a weak crystal as on a strong one, and the quantile fallback that used to paper over it is gone.
* 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.
image_analysis/scale_merge/SearchSpaceGroup.cpp
+14 -27
View File
@@ -246,26 +246,16 @@ SearchSpaceGroupResult SearchSpaceGroup(
pass_absence[i] = in_range;
}
// A fixed I/sigma cut is a statement about the error model as much as about a reflection: on a weak
// merge, whose error-model sigmas are large for every reflection, NOTHING passes, every operator is
// left with no pairs and the point group collapses to 1. Measured over the rotation battery, the four
// crystals that lose symmetry are exactly the four whose search merge is weakest - just above this
// cut - and that drop below it when the integration background changes.
// So cap the cut at the merge's own I/sigma quantile: the correlation stage always keeps at least
// its strongest quarter. This is a no-op on any merge where the fixed cut already keeps that many.
constexpr double MIN_PRESENT_FRACTION = 0.25;
// present_i_over_sigma is a cut on the reflection's own significance, and the merged I/sigma is not
// that: it carries the error model's (b*I)^2 term, so it saturates at ISa = 1/b for a reflection
// measured once and stops rising with the intensity above that knee. Convert the cut to the
// quantity the merge exports, once, here - see SearchSpaceGroupOptions::merge_isa for the
// derivation and for what it is worth. An unknown ISa (or b = 0, no systematic term) leaves the cut
// exactly where the caller set it.
double present_cut = opt.present_i_over_sigma;
{
std::vector<double> v;
v.reserve(n);
for (size_t i = 0; i < n; ++i)
if (pass_absence[i])
v.push_back(IoverSigma[i]);
if (!v.empty()) {
const size_t k = static_cast<size_t>((1.0 - MIN_PRESENT_FRACTION) * (v.size() - 1));
std::nth_element(v.begin(), v.begin() + k, v.end());
present_cut = std::min(present_cut, v[k]);
}
if (opt.merge_isa > 0.0) {
const double r = opt.present_i_over_sigma / opt.merge_isa;
present_cut = opt.present_i_over_sigma / std::sqrt(1.0 + r * r);
}
// The correlation stage uses only genuinely-present reflections. Near-zero (systematically
@@ -341,10 +331,9 @@ SearchSpaceGroupResult SearchSpaceGroup(
// DEFINE pass_cc, so normalising that over pass_cc would be circular.
//
// The price is a new dependence in place of the old one: the CC now moves with whatever defines
// pass_cc, i.e. with present_i_over_sigma and the MIN_PRESENT_FRACTION cap above. On the battery
// that is inert - the cap fires on no crystal that has an operator to score - but the cap exists
// for weak merges, so a merge on which it fires decides its point group on a differently
// normalised statistic.
// pass_cc, i.e. with present_i_over_sigma and the merge_isa conversion above. That conversion is
// what keeps the dependence harmless - it holds the cut at one counting significance on every
// crystal, so the population this is normalised over means the same thing on all of them.
const std::vector<double> Ecc = shell_normalised(pass_cc);
std::unordered_map<HKLKey, int, HKLKeyHash> key_to_index;
@@ -900,10 +889,8 @@ SearchSpaceGroupResult SearchSpaceGroup(
// present_cut, not the fixed cut: on a merge weak enough that nothing clears the fixed cut,
// screw_violations is identically zero, so every screw axis passes unchallenged, and
// present_strong is zero, so the centering rescue below switches itself off on exactly the
// weak data it exists for. Stage A already caps the cut at the merge's own 75th percentile;
// reusing it here keeps the two stages on one definition. Where the fixed cut is already the
// smaller of the two - any merge with a healthy I/sigma - present_cut EQUALS it and this is
// a no-op.
// weak data it exists for. present_cut is the same cut converted to the counting scale the
// two stages share (see merge_isa); on a healthy merge it is the fixed cut to within 1%.
const bool present = IoverSigma[i] > present_cut &&
(opt.present_e_squared <= 0.0 || Esq[i] > opt.present_e_squared);
image_analysis/scale_merge/SearchSpaceGroup.h
+25 -3
View File
@@ -284,11 +284,33 @@ struct SearchSpaceGroupOptions {
// Signed I/sigma above which a reflection counts as genuinely present. Signed, so negative
// noise is never mistaken for a real reflection. Used both to spot reflections that violate a
// wrongly-assumed absence, and to keep the correlation stage from pairing near-zero reflections.
// For the CORRELATION stage this is an upper bound on the cut rather than the cut itself: merged
// sigma is floored at b|I|, so no reflection can read above ISa = 1/b, and a merge whose ISa fell
// below this value would otherwise contribute nothing at all (see SearchSpaceGroup.cpp).
// Read against the reflection's COUNTING sigma, not against the merged one - see merge_isa, which
// is what converts it into a cut on the quantity the merge actually exports.
double present_i_over_sigma = 3.0;
// ISa = 1/b of the error model the merge handed in was scaled with: the I/sigma one observation
// of an arbitrarily strong reflection can reach. 0 = not known, or b = 0 (no systematic term at
// all), and then present_i_over_sigma is used as it stands.
//
// It is needed because a merged I/sigma is not a measure of how strong a reflection is. The merge
// carries the error model's intensity-proportional term, sigma^2 = a*sigma0^2 + (b*I)^2, so
// I/sigma saturates: at ISa*sqrt(n) for a reflection observed n times, and at ISa exactly for one
// observed once. Above that knee the quantity stops rising with the intensity - measured over the
// rotation battery's space-group search merges, the median I/sigma of the top E^2 decile sits at
// 0.97 of ISa on the weakest merge, and its log-log slope against E^2 over the strong half of the
// merge is 0.008, i.e. flat. An absolute cut on it is therefore a different test on every crystal,
// and on a merge whose ISa is near the cut it is close to unreachable: 3% of reflections clear 3.0
// there, against a battery median of 53%.
//
// So the cut is converted to the counting scale once, using the reflection's own significance.
// For a reflection observed once sigma_counting^2 = sigma^2 - (b*I)^2, hence
// I/sigma_counting = t/sqrt(1-(t/ISa)^2) with t = I/sigma, and "I/sigma_counting >= T" is exactly
// "t >= T/sqrt(1+(T/ISa)^2)". That converted cut lies strictly below ISa for every ISa, so it is
// always reachable, and it is within 1% of T on any merge with ISa >= 21, so it leaves a healthy
// merge where it was. Multiplicity is taken as 1 deliberately and not estimated: the cut has to be
// reachable by the reflection the ceiling binds hardest, which is the one measured once.
double merge_isa = 0.0;
// Extra intensity gate for the systematic-absence test: a reflection also has to reach this
// resolution-normalised intensity E^2 = I / <I>(shell) to count as violating a predicted absence.
// The error model can under-estimate sigma on weak axial reflections and fake a high I/sigma, so a
rugnux/Rugnux.cpp
+12 -5
View File
@@ -2653,11 +2653,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
if (sum / static_cast<double>(hi - lo) < 1.0) {
// Cut the noise-dominated high-res shells. When even the lowest-res shell fails,
// keep that shell alone rather than abandoning the cut. The bound is absolute
// while the merged I/sigma it tests is capped by the merge's own ISa, so a merge
// whose ISa has fallen below 1 cannot satisfy it in ANY shell - and that is
// exactly the merge the search must not be handed whole, since the noise which
// drove ISa down is what the cut exists to remove. Abandoning it feeds back:
// more noise into the search, a weaker operator statistic, and a lost symmetry.
// while the merged I/sigma it tests saturates at the merge's own ISa (times the
// square root of the multiplicity), so a merge whose ISa has collapsed towards 1
// can fail it in every shell - and that is exactly the merge the search must not
// be handed whole, since the noise which drove ISa down is what the cut exists to
// remove. Abandoning it feeds back: more noise into the search, a weaker operator
// statistic, and a lost symmetry.
if (b > 0)
d_min_search = rs[lo].first;
else if (per >= 400)
@@ -2692,6 +2693,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
// indexing. Centering is not constrained here - it is determined from the absences.
if (end_msg.rotation_lattice_type.has_value())
sg_opts.lattice_system = end_msg.rotation_lattice_type->crystal_system;
// The error model this merge was scaled with, so the search's "genuinely present" cut can be
// read on the counting sigma instead of on the merged one (SearchSpaceGroupOptions::merge_isa).
// Set per merge, not once: the two arms below are scaled separately and their ISa differ by up
// to a factor of six. 0 (b not measurable) leaves the cut as the caller set it.
sg_opts.merge_isa = result.error_model_isa;
auto sg_search = SearchSpaceGroup(sm.merged, sg_opts);
// Second opinion from a merge that keeps only well-measured observations (see
@@ -2711,6 +2717,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
rsm->SetSearchMinZeta(0.0);
auto sm_all = scale_and_merge("P1, all observations", true);
rsm->SetSearchMinZeta(zeta);
sg_opts.merge_isa = result.error_model_isa; // this arm's own error model
const auto alt = SearchSpaceGroup(sm_all.merged, sg_opts);
// The order of the point group each arm confirmed. Read from the search, not from the
// space group it went on to pick: Stage B leaves best_space_group unset when no