diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 8b817a44e..3b99f6e3e 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -1377,6 +1377,8 @@ SearchSpaceGroupResult SearchSpaceGroup( // reason as h_ratio: to report the adopted candidate's whether or not it decided. double r_added = std::numeric_limits::quiet_NaN(); double r_over_floor = std::numeric_limits::quiet_NaN(); + // The candidate's OWN operators, worst R over best R (see max_operator_r_spread). + double r_spread = std::numeric_limits::quiet_NaN(); // Report-only, for the finalist ledger: the subgroup b the tests below compared // this candidate against, whether it survived every consistency test, and - when // it was the highest refusal - why it did not. The selection loop computes all @@ -1629,6 +1631,23 @@ SearchSpaceGroupResult SearchSpaceGroup( if (n_new > 0 && pairs_new >= opt.min_pairs_for_h) r_added = r_new / n_new; } + // Homogeneity of the candidate's OWN operators: the worst-agreeing over the best-agreeing, + // over all of them, parent and added alike (see max_operator_r_spread). + { + double r_lo = std::numeric_limits::infinity(), r_hi = 0.0; + int n_op = 0; + for (const auto &rot : c.pg->rotations) { + if (rot.rot == gemmi::Op::identity().rot) + continue; + const auto &os = operator_score(rot); + if (os.n_pairs < opt.min_pairs_per_operator || !std::isfinite(os.r_stat) + || os.r_stat <= 0.0) + continue; + r_lo = std::min(r_lo, os.r_stat); r_hi = std::max(r_hi, os.r_stat); ++n_op; + } + if (n_op >= 2) + c.r_spread = r_hi / r_lo; + } c.parent_b_used = parent_b; c.r_added = r_added; c.r_over_floor = (std::isfinite(r_added) && std::isfinite(result.merge_r_floor) @@ -1686,16 +1705,28 @@ SearchSpaceGroupResult SearchSpaceGroup( if (r_refused) consistent = false; + // Operator-homogeneity test. A necessary condition like H and the contrast, and the one that + // catches what a MEAN over the added operators cannot: a candidate whose added set mixes real + // operators with false ones. The gates above read summaries - mean H, mean R - so a promotion + // that adds two operators at the noise floor together with four several times worse reads a + // mean between them and passes, while the group it names is a mixture that cannot be a + // symmetry of one crystal. Asking instead that the candidate's operators agree with EACH OTHER + // sees the mixture directly (see max_operator_r_spread). + const bool spread_refused = std::isfinite(c.r_spread) && c.r_spread > opt.max_operator_r_spread; + if (spread_refused) + consistent = false; + if (!consistent) { // Same three sentences the highest refusal gets below, but kept per candidate for the // ledger: short, because the long form is written once for the group the user is told about. - c.why = h_refused ? "H ratio" : r_refused ? "added-operator R" : "merge chi^2"; + c.why = h_refused ? "H ratio" : r_refused ? "added-operator R" + : spread_refused ? "operator spread" : "merge chi^2"; // Record the highest-order refusal so the caller can say WHY it is processing lower. if (c.order > refused_order && c.pg->representative) { refused_order = c.order; refused_pg_hm = c.pg->representative->point_group_hm(); refused_pg_rep = c.pg->representative; - refused_twin_gate = h_refused || r_refused; + refused_twin_gate = h_refused || r_refused || spread_refused; refused_fraction = h_refused ? std::max(0.0, 0.5 - h_added) : std::numeric_limits::quiet_NaN(); if (h_refused) @@ -1727,6 +1758,13 @@ SearchSpaceGroupResult SearchSpaceGroup( + FormatDouble(opt.max_operator_r_over_floor, 1) + "x) - at full resolution the added operator relates reflections that " "disagree far beyond measurement error, i.e. it is not a real symmetry"; + else if (spread_refused) + refused_why = "its own operators disagree among themselves - the worst-agreeing " + "reads " + FormatDouble(c.r_spread, 1) + + "x the intensity-weighted R of the best-agreeing one in the SAME " + "group (bound " + FormatDouble(opt.max_operator_r_spread, 1) + + "x) - so the group pools operators the data support with operators " + "they do not, and cannot be a symmetry of one crystal"; else if (std::isfinite(c.chi2) && std::isfinite(chi2_ref)) refused_why = "merge chi^2 is " + FormatDouble(c.chi2 / chi2_ref, 2) + "x the subgroup's (bound " + FormatDouble(opt.max_merge_chi2_ratio, 2) + ")"; @@ -1822,6 +1860,7 @@ SearchSpaceGroupResult SearchSpaceGroup( / (result.random_pairing_r - global_best_r); } e.r_over_floor = c.r_over_floor; + e.r_spread = c.r_spread; e.adopted = (c.pg == best_pg); e.eligible = c.eligible; e.refused_reason = c.why; @@ -2589,14 +2628,15 @@ std::string FinalistLedgerToText(const SearchSpaceGroupResult& result) { os << " Unrelated reflections (shell-matched random pairing) R " << FormatDouble(result.random_pairing_r, 4) << " - the far end of the contrast scale.\n"; - os << "\n PG ord minCC R_add contrast R/best R/floor | H chi2/best b b/par verdict\n" - " <---- calibrated evidence ----> | <--- diagnostic only, NOT evidence --->\n"; + os << "\n PG ord minCC R_add contrast spread R/best R/floor | H chi2/best b b/par verdict\n" + " <------- calibrated evidence -------> | <--- diagnostic only, NOT evidence --->\n"; for (const auto& e : result.point_group_ledger) { os << " " << std::left << std::setw(8) << e.point_group_hm << std::right << std::setw(3) << e.order << std::setw(8) << LedgerCell(e.min_class_cc, 3) << std::setw(8) << LedgerCell(e.r_added, 4) << std::setw(9) << LedgerCell(e.r_contrast, 2) + << std::setw(7) << LedgerCell(e.r_spread, 2) << std::setw(7) << LedgerCell(e.r_over_best, 2) << std::setw(8) << LedgerCell(e.r_over_floor, 2) << " |" << std::setw(6) << LedgerCell(e.h_ratio, 2) diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 4aba23e2d..1c7bd85f6 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -471,6 +471,26 @@ struct SearchSpaceGroupOptions { // the floor is a noisier reference than a confirmed operator (it can also be inflated by an // indexing ambiguity that mixes hands), and the safe direction here is to stay in P1. double max_operator_r_over_floor = 3.5; + // + // - max_operator_r_spread: the two tests above read a MEAN over the operators a promotion adds, + // so they are blind to an added set that MIXES real operators with false ones - the mean sits + // between the two populations and passes. That is what a pseudo-symmetric metric produces: a + // monoclinic crystal whose lattice is exactly hexagonal (a = c, beta = 120) carries two genuine + // horizontal 2-folds among the six a 622 promotion adds, and their R at the noise floor buys + // the other four a mean well inside the contrast bound. The mixture is invisible to any summary + // of the added set and obvious the moment the candidate's operators are compared with EACH + // OTHER: a point group is the claim that all of its operators are symmetries of one crystal, so + // they must agree with the merge about equally well. + // + // Measured over the 129 rotation datasets of one battery run that adopted a group of order 3 or + // more, as worst-agreeing over best-agreeing operator WITHIN the adopted group: every correct + // adoption sits at 1.0-2.64, the highest being the weak cubic crystal whose near-spindle 3-fold + // is unusually clean (the same crystal that motivated the contrast above), then 2.30 for a + // fine-sliced tetragonal one. The one over-called point group on that run reads 4.87. 3.2 sits + // in that gap with 0.56 of headroom below and 1.67 above - deliberately well clear of the + // 2.0-2.6 spread a genuine group reaches, which is exactly where the old ratio-to-best gate's + // bound of 2.0 sat and why that gate had to go. + double max_operator_r_spread = 3.2; // Above this reduced chi^2 for the best subgroup (chi2_ref), the merged error model is treated as // badly miscalibrated (weak, low-resolution data whose sigmas are far too small): the fixed-sigma @@ -680,6 +700,7 @@ struct PointGroupLedgerEntry { double r_over_best = std::numeric_limits::quiet_NaN(); // ... over the best-agreeing operator anywhere double r_contrast = std::numeric_limits::quiet_NaN(); // ... 0 = unrelated reflections, 1 = that best operator double r_over_floor = std::numeric_limits::quiet_NaN(); // ... over the merge's random-noise R floor + double r_spread = std::numeric_limits::quiet_NaN(); // this group's own worst operator R over its best bool adopted = false; // this is the group the search took bool eligible = false; // it passed every consistency test (so it COULD have been taken) std::string refused_reason; // why it did not, when it was the highest refusal