diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index ec434941..e407cd56 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -271,6 +271,54 @@ SearchSpaceGroupResult SearchSpaceGroup( return op_cache.emplace(rk, score_operator(op)).first->second; }; + // Conjugate rotations (symmetry-equivalent within the point group) relate symmetry-equivalent + // reflection sets, so on real data their CCs cluster; a noisy crystal can push one class member + // below min_operator_cc while the class is unmistakably present (e.g. one cubic 3-fold at 0.48 + // among siblings at 0.53-0.66). Judge each conjugacy class by its mean CC, not its weakest + // member, so a genuine high-symmetry point group is not lost to one marginal operator. chi2_under + // (below) remains the safety net against a truly false promotion. Returns {all classes present, + // worst class-mean CC}. + auto point_group_present = [&](const std::vector& rots) -> std::pair { + const size_t m = rots.size(); + std::vector cls(m, -1); + int n_cls = 0; + for (size_t i = 0; i < m; ++i) { + if (cls[i] >= 0) + continue; + cls[i] = n_cls; + for (size_t j = i + 1; j < m; ++j) + if (cls[j] < 0) + for (const auto& p : rots) + if ((p * rots[i] * p.inverse()).rot == rots[j].rot) { + cls[j] = n_cls; + break; + } + ++n_cls; + } + bool ok = true; + double worst_mean = 1.0; + for (int c = 0; c < n_cls; ++c) { + // Average over the class members that actually have enough pairs to score; a single + // low-multiplicity / degenerate (NaN) operator in an otherwise strong class is skipped, + // not allowed to veto the class. The class must still have at least one scored member. + double sum_cc = 0.0; + int n_valid = 0; + for (size_t i = 0; i < m; ++i) + if (cls[i] == c) { + const auto& s = operator_score(rots[i]); + if (s.n_pairs >= opt.min_pairs_per_operator && std::isfinite(s.cc)) { + sum_cc += s.cc; + ++n_valid; + } + } + const double mean_cc = n_valid > 0 ? sum_cc / n_valid : 0.0; + worst_mean = std::min(worst_mean, mean_cc); + if (n_valid == 0 || mean_cc < opt.min_operator_cc) + ok = false; + } + return {ok, worst_mean}; + }; + std::optional holohedry; if (opt.lattice_system.has_value()) holohedry = HolohedryRotationSet(opt.lattice_system.value()); @@ -319,41 +367,44 @@ SearchSpaceGroupResult SearchSpaceGroup( }; // 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; }; + struct PGCand { const PointGroupInfo* pg; int order; double min_class_cc; double chi2; }; std::vector pg_cands; 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(); - for (const auto& op : pg.rotations) { - const auto& s = operator_score(op); - all_present = all_present && s.present; - min_cc = std::min(min_cc, s.cc); - } - if (!all_present) + const auto [present, min_class_cc] = point_group_present(pg.rotations); + if (!present) continue; 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}); + pg_cands.push_back({&pg, static_cast(pg.rotations.size()) + 1, min_class_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 - // 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. + // chi^2 is not inflated past the miscalibration-widened bound below; ties -> higher min class 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) { + // Widen the chi^2 ratio bound by the error-model miscalibration. When even the best candidate's + // reduced chi^2 (chi2_ref) is far above 1 (weak data / uncorrected decay), the ratio grows with + // point-group order for genuine high symmetry too, so the fixed bound wrongly rejects it. Adding + // log10(chi2_ref) lets the bound exceed max_merge_chi2_ratio by ~the order of magnitude of the + // miscalibration - enough to keep the true group (F432 weak data: chi2_ref~150, ratio 3.2 vs + // bound ~4.2) while STILL rejecting a false operator whose ratio is far worse than its subgroup's + // (a twin/pseudo R32 on calibrated data: chi2_ref~2.6, ratio 4.2 vs bound ~2.4). The bound stays + // a per-candidate chi^2 test, so chi^2 still arbitrates between candidates - it is never bypassed. + const double allowed_ratio = opt.max_merge_chi2_ratio + std::log10(std::max(1.0, chi2_ref)); const bool consistent = c.pg->rotations.empty() || !std::isfinite(c.chi2) || - !std::isfinite(chi2_ref) || c.chi2 <= chi2_ref * opt.max_merge_chi2_ratio; + !std::isfinite(chi2_ref) || c.chi2 <= chi2_ref * allowed_ratio; if (!consistent) continue; - if (c.order > best_pg_order || (c.order == best_pg_order && c.min_cc > best_pg_min_cc)) { + if (c.order > best_pg_order || (c.order == best_pg_order && c.min_class_cc > best_pg_min_cc)) { best_pg = c.pg; best_pg_order = c.order; - best_pg_min_cc = c.min_cc; + best_pg_min_cc = c.min_class_cc; } } diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index c512525f..021e38ea 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -79,6 +79,9 @@ struct SearchSpaceGroupOptions { // 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. + // This bound is the calibrated-error baseline; SearchSpaceGroup widens it by log10(chi2_ref) so a + // badly miscalibrated error model (weak data / uncorrected decay, where the ratio grows with + // point-group order for real symmetry too) does not spuriously reject a genuine high-symmetry group. double max_merge_chi2_ratio = 2.0; // --- Stage B: space group (screw axes / centering) ---