diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index a6b48cd3..4835739d 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -173,7 +173,7 @@ SearchSpaceGroupResult SearchSpaceGroup( // Flatten the reflections and mark which ones each stage may use. The correlation stage drops // weak reflections; the absence stage must keep them - that is where the screw-axis signal is. std::vector H(n), K(n), L(n); - std::vector I(n), IoverSigma(n); + std::vector I(n), Sigma(n), IoverSigma(n); std::vector key(n); std::vector pass_absence(n, 0), pass_cc(n, 0); @@ -181,6 +181,7 @@ SearchSpaceGroupResult SearchSpaceGroup( const auto& r = merged[i]; H[i] = r.h; K[i] = r.k; L[i] = r.l; I[i] = r.I; + Sigma[i] = std::isfinite(r.sigma) && r.sigma > 0 ? r.sigma : 0.0; key[i] = Canonicalize(r.h, r.k, r.l, opt.merge_friedel); const bool finite = std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0 && @@ -250,14 +251,18 @@ SearchSpaceGroupResult SearchSpaceGroup( holohedry = HolohedryRotationSet(opt.lattice_system.value()); const auto point_groups = EnumeratePointGroups(holohedry); - // 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) + // Reduced chi^2 of the intensities merged under a point group's rotations - how well its symmetry + // equivalents agree RELATIVE TO THEIR ERRORS. A real point group gives ~1; a false operator forces + // non-equivalent reflections together, so they disagree by many sigma and chi^2 blows up. This is + // more sensitive than R-meas to a strong pseudo-symmetry (where the intensities still correlate well + // - high operator CC - but not within their errors). Inverse-variance weighted mean per orbit, over + // the present (pass_cc) reflections. + auto chi2_under = [&](const std::vector& rotations) -> double { + struct Acc { double sw = 0.0, swI = 0.0; int n = 0; }; + std::unordered_map grp; std::vector rep(n); for (size_t i = 0; i < n; ++i) { - if (!pass_cc[i]) + if (!pass_cc[i] || !(Sigma[i] > 0.0)) continue; HKLKey best = key[i]; for (const auto& op : rotations) { @@ -268,30 +273,30 @@ SearchSpaceGroupResult SearchSpaceGroup( } rep[i] = best; auto& g = grp[best]; - g.first += 1; g.second += I[i]; + const double w = 1.0 / (Sigma[i] * Sigma[i]); + g.sw += w; g.swI += w * I[i]; g.n += 1; } - std::unordered_map absdev; + double chi2 = 0.0; + long dof = 0; for (size_t i = 0; i < n; ++i) { - if (!pass_cc[i]) + if (!pass_cc[i] || !(Sigma[i] > 0.0)) 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) + if (g.n < 2) continue; - num += std::sqrt(static_cast(g.first) / (g.first - 1)) * absdev[k]; - den += g.second; + const double mean = g.swI / g.sw, dev = I[i] - mean; + chi2 += dev * dev / (Sigma[i] * Sigma[i]); } - return den > 0 ? num / den : std::numeric_limits::quiet_NaN(); + for (const auto& [k, g] : grp) + if (g.n >= 2) + dof += g.n - 1; + return dof > 0 ? chi2 / static_cast(dof) : 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; }; + // Operator-CC-confirmed candidates, each with its merge chi^2; chi2_ref = the most consistent. + struct PGCand { const PointGroupInfo* pg; int order; double min_cc; double chi2; }; std::vector pg_cands; - double rmeas_ref = std::numeric_limits::infinity(); + double chi2_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(); @@ -302,22 +307,22 @@ 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 double ch = pg.rotations.empty() ? std::numeric_limits::quiet_NaN() + : chi2_under(pg.rotations); + pg_cands.push_back({&pg, static_cast(pg.rotations.size()) + 1, min_cc, ch}); + if (!pg.rotations.empty() && std::isfinite(ch)) + chi2_ref = std::min(chi2_ref, ch); } // 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. + // chi^2 is not inflated past max_merge_chi2_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; + const bool consistent = c.pg->rotations.empty() || !std::isfinite(c.chi2) || + !std::isfinite(chi2_ref) || c.chi2 <= chi2_ref * opt.max_merge_chi2_ratio; if (!consistent) continue; if (c.order > best_pg_order || (c.order == best_pg_order && c.min_cc > best_pg_min_cc)) { diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 89b87a48..82109b9c 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -71,12 +71,15 @@ 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; + // Per-operator CC alone cannot tell a real weak operator from a false strong one (a noisy crystal's + // genuine 2-fold can score below a pseudo-symmetric crystal's near-perfect false one). So a point + // group is also required to be SELF-CONSISTENT: merging the intensities under it must not inflate the + // reduced chi^2 (within-orbit scatter / sigma^2) beyond this factor times the most-consistent + // candidate. A false operator forces non-equivalent reflections together so they disagree by many + // sigma and chi^2 blows up; a real one leaves it ~flat even when the operator CC is only moderate. + // 2.0 separates the test set: a true point group's chi^2 stays within ~1.6x the best subgroup + // (the cubic merge of a noisy crystal is the worst real case), while a false operator gives >=2.3x. + double max_merge_chi2_ratio = 2.0; // --- Stage B: space group (screw axes / centering) --- bool determine_space_group = true; // false: stop at the symmorphic representative