diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index c381f5f1..9d108de5 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -273,16 +273,18 @@ SearchSpaceGroupResult SearchSpaceGroup( s.n_pairs = static_cast(x.size()); s.cc = PearsonCC(x, y); // Sigma-free disagreement over the same pairs (see SpaceGroupOptions::max_operator_h_ratio). - double h_sum = 0.0; - int h_n = 0; + std::vector hv; + hv.reserve(x.size()); for (size_t p = 0; p < x.size(); ++p) { const double denom = x[p] + y[p]; - if (denom > 0.0) { - h_sum += std::fabs(x[p] - y[p]) / denom; - ++h_n; - } + if (denom > 0.0) + hv.push_back(std::fabs(x[p] - y[p]) / denom); + } + if (!hv.empty()) { + const size_t mid = hv.size() / 2; + std::nth_element(hv.begin(), hv.begin() + mid, hv.end()); + s.h_stat = hv[mid]; } - s.h_stat = h_n > 0 ? h_sum / h_n : 0.0; s.present = s.n_pairs >= opt.min_pairs_per_operator && std::isfinite(s.cc) && s.cc >= opt.min_operator_cc; return s; @@ -536,7 +538,7 @@ SearchSpaceGroupResult SearchSpaceGroup( const bool in_parent = std::binary_search(parent->rotation_set.begin(), parent->rotation_set.end(), RotKey(rot)); if (in_parent) { h_par += os.h_stat; ++n_par; pairs_par += os.n_pairs; } - else { h_new += os.h_stat; ++n_new; pairs_new += os.n_pairs; } + else { h_new += os.h_stat; ++n_new; pairs_new += os.n_pairs; } } if (n_new > 0 && n_par > 0 && pairs_new >= opt.min_pairs_for_h && pairs_par >= opt.min_pairs_for_h && h_par > 0.0) { diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index bafa48ca..ad630672 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -30,11 +30,12 @@ struct SpaceGroupOperatorScore { double cc = 0.0; // correlation of I(h) with I(Rh) int n_pairs = 0; // independent reflection pairs the CC was computed from bool present = false; // operator confirmed as a real symmetry of the intensities - // Mean |I1-I2|/(I1+I2) over this operator's pairs - the disagreement the operator implies, with no + // MEDIAN |I1-I2|/(I1+I2) over this operator's pairs - the disagreement the operator implies, with no // sigma in it. Unlike chi^2 and the systematic-b, which are ratios to a merge error model that // drifts with multiplicity, this is a property of the intensities alone; compared against an // operator already confirmed on the same reflections it is what separates a real symmetry from a - // twin law (see max_operator_h_ratio). + // twin law (see max_operator_h_ratio). The median rather than the mean because a merohedral twin + // perturbs EVERY pair while a badly integrated minority perturbs only the tail. double h_stat = 0.0; }; @@ -131,16 +132,20 @@ struct SearchSpaceGroupOptions { // is a non-negligible fraction of I; below the floor the increase is treated as noise, not a twin. double min_systematic_b_for_veto = 0.05; - // Promotion gate on the operator disagreement H = <|I1-I2|/(I1+I2)>, taken as the ratio of the + // Promotion gate on the operator disagreement H = median|I1-I2|/(I1+I2), taken as the ratio of the // operators the promotion ADDS to the operators of the parent group already confirmed on the same // reflections. A real symmetry operator relates equal intensities, so its H matches the parent's // (ratio ~1); a merohedral twin law relates DIFFERENT reflections mixed in proportion alpha, so its - // H is systematically larger. Measured over 27 runs spanning 5 promotion types and 450-1800 images: - // genuine symmetry 0.862-1.219, merohedral twins 1.270-2.084 - a clean gap, and unlike the chi^2 and - // systematic-b ratios (genuine 1.00-3.47 / 1.09-3.89 vs twin 1.35-3.32 / 1.77-4.76, fully - // interleaved) it does not drift with multiplicity, because there is no sigma in it and the parent - // normalisation cancels data quality. For a partial twin over the twin law is (1-2*alpha)/2, so - // the excess over the parent also estimates the twin fraction rather than being a tuned constant. + // H is systematically larger. Unlike the chi^2 and systematic-b ratios (genuine 1.00-3.47 / + // 1.09-3.89 vs twin 1.35-3.32 / 1.77-4.76, fully interleaved) it does not drift with multiplicity, + // because there is no sigma in it and the parent normalisation cancels data quality. + // + // A MEDIAN, not a mean. A merohedral twin mixes every reflection with its twin mate, so it shifts the + // whole distribution of |I1-I2|/(I1+I2); a minority of badly measured reflections shifts only the + // tail. Measured on the same real crystals, moving from the mean to the median leaves genuine + // promotions where they are (1.016 -> 1.013, 1.051 -> 1.067, 1.238 -> 1.231) and pushes every real + // twin UP (1.280 -> 1.622, 1.427 -> 2.010, 1.272 -> 1.447, 1.441 -> 1.522), widening the margin + // around this bound from 2.7% to 17.5%. Battery-neutral: 33 crystals, no point group changed. // // The parent normalisation is what makes this work, and it is not optional. Symmetry-related // reflections never agree exactly on real data - absorption, illumination and partiality differ diff --git a/tests/SearchSpaceGroupTwinTest.cpp b/tests/SearchSpaceGroupTwinTest.cpp index 43f7402b..e9a46e88 100644 --- a/tests/SearchSpaceGroupTwinTest.cpp +++ b/tests/SearchSpaceGroupTwinTest.cpp @@ -149,9 +149,9 @@ namespace { return out; } - // The statistic the promotion is actually decided on: the mean operator disagreement - // H = <|I1-I2|/(I1+I2)> over the operators the promotion ADDS, divided by the mean over the parent - // group's own operators, measured on the same reflections. Mirrors what SearchSpaceGroup computes + // The statistic the promotion is actually decided on: the operator disagreement + // H = median|I1-I2|/(I1+I2) over the operators the promotion ADDS, divided by the same over the + // parent group's own operators, measured on the same reflections. Mirrors what SearchSpaceGroup computes // for the sub -> super step, so a test can report the margin the max_operator_h_ratio bound has. // // The parent normalisation is the whole design, not a detail. An ABSOLUTE per-operator bound cannot @@ -172,6 +172,7 @@ namespace { for (const auto& s : operators) { if (s.n_pairs < 200) // SearchSpaceGroupOptions::min_pairs_for_h continue; + // The median, which is what the promotion is gated on (SearchSpaceGroup). if (parent_ops.count(s.op_triplet_hkl) > 0) { h_parent += s.h_stat; ++n_parent; } else { h_added += s.h_stat; ++n_added; } }