diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index cc1535a8..a6b48cd3 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -250,11 +250,48 @@ SearchSpaceGroupResult SearchSpaceGroup( holohedry = HolohedryRotationSet(opt.lattice_system.value()); const auto point_groups = EnumeratePointGroups(holohedry); - // Choose the largest point group all of whose operators are confirmed (ties -> higher min CC). - const PointGroupInfo* best_pg = nullptr; - int best_pg_order = 0; - double best_pg_min_cc = -2.0; + // R-meas of the intensities merged under a point group's rotations - how well its symmetry + // equivalents agree. A real point group keeps this low; a false operator forces non-equivalent + // reflections together and inflates it. Computed over the present (pass_cc) reflections. + auto rmeas_under = [&](const std::vector& rotations) -> double { + std::unordered_map, HKLKeyHash> grp; // orbit rep -> (n, sum I) + std::vector rep(n); + for (size_t i = 0; i < n; ++i) { + if (!pass_cc[i]) + continue; + HKLKey best = key[i]; + for (const auto& op : rotations) { + const auto m = op.apply_to_hkl(gemmi::Op::Miller{{H[i], K[i], L[i]}}); + const HKLKey k2 = Canonicalize(m[0], m[1], m[2], opt.merge_friedel); + if (std::make_tuple(k2.h, k2.k, k2.l) < std::make_tuple(best.h, best.k, best.l)) + best = k2; + } + rep[i] = best; + auto& g = grp[best]; + g.first += 1; g.second += I[i]; + } + std::unordered_map absdev; + for (size_t i = 0; i < n; ++i) { + if (!pass_cc[i]) + continue; + const auto& g = grp[rep[i]]; + if (g.first >= 2) + absdev[rep[i]] += std::fabs(I[i] - g.second / g.first); + } + double num = 0, den = 0; + for (const auto& [k, g] : grp) { + if (g.first < 2) + continue; + num += std::sqrt(static_cast(g.first) / (g.first - 1)) * absdev[k]; + den += g.second; + } + return den > 0 ? num / den : std::numeric_limits::quiet_NaN(); + }; + // Operator-CC-confirmed candidates, each with its merge R-meas; rmeas_ref = the most consistent. + struct PGCand { const PointGroupInfo* pg; int order; double min_cc; double rmeas; }; + std::vector pg_cands; + double rmeas_ref = std::numeric_limits::infinity(); for (const auto& pg : point_groups) { bool all_present = true; double min_cc = pg.rotations.empty() ? 1.0 : std::numeric_limits::infinity(); @@ -265,12 +302,28 @@ SearchSpaceGroupResult SearchSpaceGroup( } if (!all_present) continue; + const double rm = pg.rotations.empty() ? std::numeric_limits::quiet_NaN() + : rmeas_under(pg.rotations); + pg_cands.push_back({&pg, static_cast(pg.rotations.size()) + 1, min_cc, rm}); + if (!pg.rotations.empty() && std::isfinite(rm)) + rmeas_ref = std::min(rmeas_ref, rm); + } - const int order = static_cast(pg.rotations.size()) + 1; - if (order > best_pg_order || (order == best_pg_order && min_cc > best_pg_min_cc)) { - best_pg = &pg; - best_pg_order = order; - best_pg_min_cc = min_cc; + // Choose the largest point group that is both operator-confirmed AND self-consistent (its merge + // R-meas is not inflated past max_rmeas_ratio x the most-consistent candidate; ties -> higher min + // CC). Identity (no operators) is always consistent, so it stays the P1 fallback. + const PointGroupInfo* best_pg = nullptr; + int best_pg_order = 0; + double best_pg_min_cc = -2.0; + for (const auto& c : pg_cands) { + const bool consistent = c.pg->rotations.empty() || !std::isfinite(c.rmeas) || + !std::isfinite(rmeas_ref) || c.rmeas <= rmeas_ref * opt.max_rmeas_ratio; + if (!consistent) + continue; + if (c.order > best_pg_order || (c.order == best_pg_order && c.min_cc > best_pg_min_cc)) { + best_pg = c.pg; + best_pg_order = c.order; + best_pg_min_cc = c.min_cc; } } diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index acdb0286..89b87a48 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -71,6 +71,13 @@ struct SearchSpaceGroupOptions { double min_operator_cc = 0.5; int min_pairs_per_operator = 20; + // Per-operator CC alone cannot tell a real weak operator from a false moderate one (a noisy + // crystal's genuine 2-fold can score below a pseudo-symmetric crystal's false one). So a point + // group is also required to be SELF-CONSISTENT: merging the intensities under it must not inflate + // R-meas beyond this factor times the most-consistent (lowest-R-meas) candidate. A false operator + // forces non-equivalent reflections together and blows R-meas up; a real one leaves it ~flat. + double max_rmeas_ratio = 1.5; + // --- Stage B: space group (screw axes / centering) --- bool determine_space_group = true; // false: stop at the symmorphic representative