diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index 4909193e1..08c19ca5b 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -50,7 +50,7 @@ Both search merges also drop the **frames whose fitted per-frame scale came out **A confirmed promotion that is refused is decided by the twin-immune zone of its added operators, and only where that is undecided by merging under it.** The zone (§13.2) is the one reading that does not depend on the per-frame scales or on a reference operator: read on the all-observation P1 merge where the two search arms disagree, and on the adopted group's own merge at the remerge, an index-2 promotion whose calibrated zone evidence is acentric by 20 nats stays refused — on all observations, and the first added operator is recorded as the twin law — and one whose zone is centric by 20 nats is taken. Where the zone cannot decide (a higher index, no zone reflections in the shells with signal, within the margin), the merge is asked. Every gate that refuses a *confirmed* higher point group is a **ratio** — the added operators' disagreement against the parent's, or their $R$ against the best-agreeing operator anywhere — and both references are properties of how the crystal was mounted rather than of its symmetry. A 2-fold within a degree of the spindle records its mates on the same detector pixel half a turn later, so it carries no geometry-dependent systematic at all, and by being that clean it makes every other operator look bad against it; the best-agreeing operator can also be one of the operators under test, which collapses the ratio on worse data and raises it on better. So where a promotion is confirmed and then refused, the **merge** is asked instead of the ratio. Both groups are scaled and merged over one pinned resolution range — the adopted group's own cut — and the promotion is taken only if **neither $R_\mathrm{meas}$ nor ISa gets worse**. $R_\mathrm{meas}$ is multiplicity-corrected, so folding non-equivalent reflections together has to inflate it, and the error model is refitted per merge, so ISa says whether the extra multiplicity was bought with systematic disagreement. $CC_{1/2}$ and $\langle I/\sigma\rangle$ cannot arbitrate this: both *rise* on a false promotion too. Two details the comparison turns on: the refused point group's representative is primitive and symmorphic, so the arm merged under it takes the **adopted group's centring** — merged in the bare representative, a centred crystal is handed its centring-absent class as data, which inverts the decision — and the absence stage is asked of the adopted merge, so both arms carry the same absence classes and differ by the added rotations alone. That merge no longer holds the adopted group's screw-absent axial reflections, so the higher group also **keeps the screws the adopted group decided** — among the higher candidates the merge cannot separate, the one that predicts every such axial absence is taken — rather than the lowest-numbered one, which claims none. The cost is two extra merges, and only on a run that records a refusal, which is a few in a hundred. -**On a merge that reads twinned, the search asks the zone itself.** A twin of an index-2 subgroup by exactly the operators a promotion adds defeats every agreement gate once its fraction is high — the added operators then agree like real ones — so where the P1 merge's $\langle|L|\rangle$ lies in the partial-twin band (0.375 to 0.44) and the candidate does not hold every rotation of its lattice (so a twin law outside it is possible), the candidate is refused when the twin-immune zone of one of its index-2 subgroups reads acentric by 20 nats. Below 0.375 something other than a twin compresses the intensities, the zones with them, and they are not read. A $P3_1$ crystal twinned at $\alpha \ge 0.4$ by its 321 law, whose 32 promotion every agreement gate passed, reads −230 nats there. +**On a merge that reads twinned, the search asks the zone itself.** A twin of an index-2 subgroup by exactly the operators a promotion adds defeats every agreement gate once its fraction is high — the added operators then agree like real ones — so where the P1 merge's $\langle|L|\rangle$ lies in the partial-twin band (0.375 to 0.44), the candidate is refused when the twin-immune zone of one of its index-2 subgroups reads acentric by 20 nats. That holds for a candidate that is the whole symmetry of its lattice too: the subgroup's twin law is one of the candidate's own operators, so a crystal in 4 twinned by $k,h,\bar{l}$ is exactly such a twin of 422. Below 0.375 something other than a twin compresses the intensities, the zones with them, and they are not read. A $P3_1$ crystal twinned at $\alpha \ge 0.4$ by its 321 law, whose 32 promotion every agreement gate passed, reads −230 nats there; a body-centred crystal in 4 twinned at about 0.3, promoted to 422 at an $H$ ratio of 1.51, −58 nats. Between two candidates of the same order the one whose own two-folds carry screw absences is preferred, but an axial row whose odd class a lattice centring extinguishes is not counted where the merge reads that centring in its general reflections: on a body-centred crystal every principal row reads extinct whatever the two-folds are, which would hand 222 the tie over 4. **A point group's own operators must agree with each other.** A point group is the claim that all of its operators are symmetries of one crystal, and a mixture of real and false operators has a specific shape: the real ones form a subgroup — the true group's intersection with the candidate — and everything outside it is false. The gate therefore reads the operators sorted by their $R$ against the merge, finds every prefix that closes into a subgroup, and takes the largest factor by which every operator outside such a subgroup reads worse than every operator inside it; above 3.2 the candidate is refused. A genuine group whose operators merely spread over a continuum — anisotropic systematic error spreads a cubic crystal's operators twofold — has no such gap. The earlier form, worst-agreeing over best-agreeing operator, divided by an extreme order statistic and refused a genuine $432$ as soon as one of its operators happened to agree unusually well. diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 59a9c3630..7625c9d39 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -1693,6 +1694,46 @@ SearchSpaceGroupResult SearchSpaceGroup( // row is the crystal's and the other is the twin law that out-correlated it. // Only a row whose axis is a two-fold and nothing more: on a 3-, 4- or 6-fold axis the odd class // is not what the screw extinguishes (a 3_1 keeps l = 3, 9, ... strong), so it says nothing here. + // Nor a row a lattice centring extinguishes, where the merge reads that centring: an I lattice + // kills the odd class of all three principal rows whatever the two-folds are, so on a primitive + // merge of a body-centred crystal every two-fold along an axis read "screwed", and a 222 made of + // a real 2-fold and two twin laws out-scored the crystal's own 4 (a tetragonal twin, 30 nats on + // each row). The centring is read on the general reflections, where no screw reaches. + using Extinction = bool (*)(int, int, int); + std::optional> centrings_read; + const auto row_extinct_by_centring = [&](const AxialRow& row) { + if (!centrings_read) { + const Extinction centrings[] = { + [](int h, int k, int l) { return (h + k + l) % 2 != 0; }, // I + [](int, int k, int l) { return (k + l) % 2 != 0; }, // A + [](int h, int, int l) { return (h + l) % 2 != 0; }, // B + [](int h, int k, int) { return (h + k) % 2 != 0; }, // C + [](int h, int k, int l) { return (h + k) % 2 != 0 || (h + l) % 2 != 0; }}; // F + centrings_read.emplace(); + for (const auto& absent : centrings) { + double sum_absent = 0.0, sum_present = 0.0; + int n_absent = 0, n_present = 0; + for (size_t i = 0; i < n; ++i) { + if (!pass_absence[i] || (H[i] == 0) + (K[i] == 0) + (L[i] == 0) >= 2) + continue; + if (absent(H[i], K[i], L[i])) { + sum_absent += Esq[i]; + ++n_absent; + } else { + sum_present += Esq[i]; + ++n_present; + } + } + // A primitive lattice leaves the class at the mean of the rest (0.97-1.00 measured); a + // real centring at noise, which a weak, high-background body-centred crystal still + // reads at 0.11-0.13 of it. + if (n_absent >= 50 && n_present >= 50 && sum_absent / n_absent < 0.25 * (sum_present / n_present)) + centrings_read->push_back(absent); + } + } + return std::any_of(centrings_read->begin(), centrings_read->end(), + [&](const auto& absent) { return absent(row[0], row[1], row[2]); }); + }; const auto screwed_two_folds = [&](const PointGroupInfo& pg) { const auto is_two_fold = [](const gemmi::Op& op) { return op.rot != gemmi::Op::identity().rot && op.combine(op).rot == gemmi::Op::identity().rot; @@ -1704,7 +1745,7 @@ SearchSpaceGroupResult SearchSpaceGroup( for (const auto& rot : pg.rotations) if (rot.rot != gemmi::Op::identity().rot && rot.apply_to_hkl(hkl) == hkl) (is_two_fold(rot) ? two_fold : higher) = true; - if (two_fold && !higher) + if (two_fold && !higher && !row_extinct_by_centring(row)) total += row_screw_evidence(row); } return total; @@ -1959,10 +2000,12 @@ SearchSpaceGroupResult SearchSpaceGroup( if (l_refused) consistent = false; - // Twin-immune zone test, for a group that is NOT the whole symmetry of its lattice, on a merge - // that reads twinned: the other half of the twin argument. Where a lattice operation lies - // outside the candidate, a twin of one of its index-2 subgroups by exactly the operators the - // candidate adds over it is possible, and at a high twin fraction it defeats every gate above - + // Twin-immune zone test, on a merge that reads twinned: the other half of the twin argument. A + // twin of one of the candidate's index-2 subgroups by exactly the operators the candidate adds + // over it is possible wherever those operators are lattice operations - which they always are, + // the whole symmetry of the lattice included: a tetragonal crystal in 4 twinned by k,h,-l is the + // commonest case of all, and the L-test gate above leaves it in 422 by design (a twin of the + // subgroup closes 0.19-0.30 of the gap). At a high twin fraction such a twin defeats every gate above - // the added operators then agree like real ones (H, R, chi^2 all read "real"), which is the // theorem that no statistic of the operators' agreement separates a near-perfect twin from the // higher group. What it cannot fake are the reflections centric in the candidate but not in the @@ -1975,7 +2018,8 @@ SearchSpaceGroupResult SearchSpaceGroup( // - Measured: a P3_1 crystal twinned at 0.4 or more by its 321 law, whose 32 promotion every // gate above passed (contrast 0.87, H ratio 0.97), reads -230 nats; a genuine P3_121 // crystal twinned by a 622 law, read at its other law's fraction 0.07, +209 (-206 read - // untwinned, which is why the fraction is read). + // untwinned, which is why the fraction is read). A body-centred tetragonal crystal in 4 + // twinned ~0.3 by k,h,-l, promoted to 422 at an H ratio of 1.51, reads -58 nats over 4. // - Only on a merge whose L-test reads twinned (merge_reads_twinned, on the P1 merge before any // promotion): without twinning a twin law cannot pass the gates above, and a crystal that is // not twinned is not asked a twin's question. Below the perfect twin's 0.375 something other @@ -1987,7 +2031,7 @@ SearchSpaceGroupResult SearchSpaceGroup( std::string zone_operators; double zone_evidence = 0.0; if (consistent && !c.pg->rotations.empty() && c.pg->representative && opt.cell.has_value() - && lattice_rotations && !holds_lattice && merge_reads_twinned()) { + && lattice_rotations && merge_reads_twinned()) { const auto zones = AnalyzeTwinImmuneZones(merged, *opt.cell, *c.pg->representative, nullptr, opt.nthreads, twin_fraction_outside(*c.pg)); if (!zones.zones_ambiguous) diff --git a/tests/SearchSpaceGroupTwinTest.cpp b/tests/SearchSpaceGroupTwinTest.cpp index acc533344..33ead7c2f 100644 --- a/tests/SearchSpaceGroupTwinTest.cpp +++ b/tests/SearchSpaceGroupTwinTest.cpp @@ -239,6 +239,29 @@ TEST_CASE("SearchSpaceGroup keeps a partially twinned crystal in its true subgro } } +// A body-centred crystal in 4, twinned near-perfectly by a 2-fold of 422, given its cell as rugnux +// gives it. The agreement gates pass 422 at this fraction and only the twin-immune zone refuses it; +// what remains is a tie of order 4 between 4 and a 222 made of the crystal's 2-fold and two twin laws. +// The I centring extinguishes the odd class of every principal row, and that must not be read as +// three screw axes for the 222. +TEST_CASE("SearchSpaceGroup keeps a near-perfectly twinned body-centred 4 in 4", "[SearchSpaceGroup][twin]") { + SyntheticMergeParams p; + p.true_space_group = "I 4"; + p.twin_supergroup = "I 4 2 2"; + p.twin_fraction = 0.45; + p.multiplicity = 6; + p.seed = 20260727; + const auto merged = jfjoch_test::GenerateSyntheticMerged(p); + + SearchSpaceGroupOptions opt; + opt.merge_friedel = true; + opt.lattice_system = gemmi::CrystalSystem::Tetragonal; + opt.cell = jfjoch_test::SyntheticCellFor(gemmi::get_spacegroup_by_name("I 4")); + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + CHECK(result.point_group_hm == "4"); +} + // A PERFECT (alpha = 0.5) merohedral twin produces intensities that are exactly invariant under the // twin law: I_obs(h) = I_obs(twin h) for every reflection. No intensity statistic can tell it from a // crystal that genuinely has the higher symmetry - the information is not in the data (it takes a