diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index cc6c37e28..9315f3193 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -114,13 +114,22 @@ namespace { // Only the leading term of the regularized incomplete beta is kept. It is exact as T -> 0, which is // where a condition is claimed, and dropping the (1-T)^n_control factor only ever UNDER-states the // evidence, which is the safe direction for a test that has to clear a bound. - double AbsenceEvidence(double sum_u, int n_absent, int n_control) { - if (n_absent <= 0 || n_control <= 0) - return 0.0; - const double a = n_absent, b = n_control; - const double T = std::max(sum_u / (sum_u + b), 1e-300); - return -(a * std::log(T) + std::lgamma(a + b) - std::lgamma(a + 1) - std::lgamma(b)); - } + // No measurement can place a merged intensity at exactly zero, so sum_u is floored at a + // thousandth of the control mean per absent reflection. Without it, a zone whose absences all + // merged non-positive gives sum_u = 0 exactly, T clamps to the epsilon below, and EACH absent + // reflection is worth ~690 nats. That was harmless while the value only had to clear a bound of + // 20, but it is now summed across zones and ranks the candidates, and it inverted the ranking + // outright: a zone of 2 absences that were never measurable scored 1378 nats where a genuine zone + // of 6 absences at 1% of its row scores 22, so the candidate claiming a screw on an UNMEASURED + // row beat the one whose rows are actually dead, by 60x. + // + // The constant is deliberately an order of magnitude below the precision any real merge reaches - + // a thousandth of the row mean needs I/sigma ~ 1000 against that row, where ISa tops out near 40 - + // so it can only ever remove the singularity, never suppress evidence a measurement could have + // produced. Measured over the realistic range it changes no genuine zone at all (22.0, 34.7 and + // 65.4 nats to four figures) and takes the unmeasurable ones to 13 and 36. + constexpr double MIN_U_PER_ABSENT_REFLECTION = 1e-3; + // The reciprocal-space row {1,0,0} written the way a crystallographer names the zone: h00. std::string RowLabel(std::array row) { @@ -306,6 +315,15 @@ namespace { } } +double AbsenceEvidence(double sum_u, int n_absent, int n_control) { + if (n_absent <= 0 || n_control <= 0) + return 0.0; + const double a = n_absent, b = n_control; + const double u = std::max(sum_u, a * MIN_U_PER_ABSENT_REFLECTION); + const double T = u / (u + b); + return -(a * std::log(T) + std::lgamma(a + b) - std::lgamma(a + 1) - std::lgamma(b)); +} + SearchSpaceGroupResult SearchSpaceGroup( const std::vector& merged, const SearchSpaceGroupOptions& opt) { @@ -1202,8 +1220,11 @@ SearchSpaceGroupResult SearchSpaceGroup( // The screw part is the SUM over the zones, not the group's single gating number: each axial row // is a separate condition tested on its own reflections, so the log-likelihoods add. Summing is // what makes an extra condition pay its own way - a group claiming one screw more gains that - // zone's evidence when the row really is dead and LOSES it when the row is not, where a count - // could only ever go up. Pooling the rows into one Beta instead would make three genuine screws + // zone's evidence when the row really is dead and gains little when the row is not, where a count + // could only ever go up. (It gains LITTLE, not nothing: a zone whose absences were never + // measurable still scores a bounded positive value, see MIN_U_PER_ABSENT_REFLECTION. Requiring a + // minimum number of absences per zone before it may contribute, as the centring class does + // through min_absent_observed, would close that too - it needs its own battery.) Pooling the rows into one Beta instead would make three genuine screws // read as weaker evidence than two whenever the third row is measured less deeply, which is a // property of the pooling, not of the crystal. auto absence_evidence = [](const SpaceGroupCandidateScore& s) { diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 45fe7a83b..de2582942 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -82,6 +82,15 @@ struct SpaceGroupCandidateScore { bool centering_untested = false; }; +// How unlikely a predicted-absent class would be if the condition producing it did not exist, in +// nats: the -log lower tail of Beta(n_absent, n_control) at T = sum_u / (sum_u + n_control), where +// sum_u is the absent intensities expressed in units of their control class's mean. Used for a SCREW +// against the rest of its own axial row and for a CENTRING against the present class, so the two are +// in the same units and add. The control's own strength cancels, which is the property a count of +// absences does not have. Declared here because it is the statistic the candidate ranking is built +// on, and the comments below name it. +double AbsenceEvidence(double sum_u, int n_absent, int n_control); + struct SearchSpaceGroupOptions { // Lattice (metric) symmetry from LatticeSearch. When set, the point-group search is limited to // the Sohncke subgroups of this system's holohedry - the metric is an upper bound on the diff --git a/tests/SearchSpaceGroupTest.cpp b/tests/SearchSpaceGroupTest.cpp index b9cd6a678..9863be0d1 100644 --- a/tests/SearchSpaceGroupTest.cpp +++ b/tests/SearchSpaceGroupTest.cpp @@ -490,3 +490,27 @@ TEST_CASE("SearchSpaceGroup does not promote triclinic data on a pseudo-orthorho CHECK(result.best_space_group->number == 1); CHECK(result.point_group_order == 1); } + +// A zone whose predicted absences were never measurable must not outscore a zone that is genuinely +// dead. sum_u is a sum of max(0, E^2)/row_mean, so it is EXACTLY zero when every absent reflection in +// the zone merged non-positive - and the Beta tail then diverges, worth ~690 nats per reflection. That +// was harmless while the number only had to clear a bound; it is now summed across zones and ranks the +// candidates, so it made a candidate claiming a screw on an UNMEASURED row beat one whose rows are +// actually dead. The evidence is scored through the same entry point for both kinds of absence. +TEST_CASE("AbsenceEvidence does not reward a zone that was never measurable", "[SearchSpaceGroup]") { + // 2 absences that all merged non-positive, against a control of 8... + const double unmeasurable = AbsenceEvidence(0.0, 2, 8); + // ...against a genuinely dead zone: 6 absences at 1% of their row's mean, same control. + const double genuine = AbsenceEvidence(0.06, 6, 8); + + CHECK(std::isfinite(unmeasurable)); + CHECK(unmeasurable < genuine); // the ordering that was inverted + CHECK(unmeasurable < 20.0); // and it does not clear min_screw_absence_evidence + + // The floor is far below any real measurement, so a genuine zone is untouched by it. + CHECK(genuine == Catch::Approx(22.0).margin(0.2)); + CHECK(AbsenceEvidence(0.22, 22, 8) == Catch::Approx(65.4).margin(0.3)); + + // More dead reflections still means more evidence, which is the property the sum relies on. + CHECK(AbsenceEvidence(0.0, 6, 8) > AbsenceEvidence(0.0, 2, 8)); +}