diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 6dcfe571..744206f0 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -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 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. diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 19d7989c..403e1202 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -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 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((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 Ecc = shell_normalised(pass_cc); std::unordered_map 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); diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 47210ccc..141b137a 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -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 / (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 diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index bd6f93ad..dfa27de4 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2653,11 +2653,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (sum / static_cast(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