diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 1e717dde5..1817ea052 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -134,6 +134,15 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotat outcome.symmetry.centering = result->search_result.centering; outcome.symmetry.niggli_class = result->search_result.niggli_class; outcome.symmetry.crystal_system = result->search_result.system; + if (const auto &tw = result->twin_overlap) { + // t is a vector of the lattice, so in a centred cell's basis its coordinates are multiples + // of 1/2 or 1/3: in sixths they are integers, and the test is then modulo 6n. + const auto &l = result->lattice; + outcome.twin_overlap_w = {static_cast(std::lround(6.0f * (tw->t * l.Astar()))), + static_cast(std::lround(6.0f * (tw->t * l.Bstar()))), + static_cast(std::lround(6.0f * (tw->t * l.Cstar())))}; + outcome.twin_overlap_n = 6 * tw->n; + } } return outcome; } @@ -587,7 +596,13 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, // caller's callback binds it to the right image (GPU-resident buffer, host buffer, or the assembled // FPGA image read straight on the CPU). auto integration_start_time = std::chrono::steady_clock::now(); - const std::vector integrated = integrate(prediction.GetReflections(), nrefl, msg.number); + std::vector integrated = integrate(prediction.GetReflections(), nrefl, msg.number); + if (const int n = outcome.twin_overlap_n; n > 0) { + const auto &w = outcome.twin_overlap_w; + std::erase_if(integrated, [&](const Reflection &r) { + return (r.h * w[0] + r.k * w[1] + r.l * w[2]) % n == 0; + }); + } msg.integrated_reflections = integrated.size(); auto integration_end_time = std::chrono::steady_clock::now(); msg.integration_time_s = std::chrono::duration(integration_end_time - integration_start_time).count(); diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index a5eb15a28..f27efe47e 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -85,6 +85,11 @@ class IndexAndRefine { .crystal_system = gemmi::CrystalSystem::Triclinic }; bool beam_center_updated = false; + // The rotation result's twin overlap (RotationIndexerResult::TwinOverlap) in the candidate's + // basis, in sixths: reflections with h*w[0] + k*w[1] + l*w[2] divisible by n (6 times the + // overlap's index) are left out of integration. + std::array twin_overlap_w{}; + int twin_overlap_n = 0; explicit IndexingOutcome(const DiffractionExperiment& experiment_ref) : experiment(experiment_ref) {} diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index bb177c12c..e077d828b 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -244,6 +244,7 @@ void RotationIndexer::RunIndexing() { std::unique_lock ul(m); if (!axis_) return; + twin_overlap_.reset(); const KeyDigest key = Digest(InputKey()); std::optional kept; @@ -886,6 +887,7 @@ std::optional RotationIndexer::GetLattice() const { .axis = axis_, .unconstrained = unconstrained_, .tilt_walk = tilt_walk_, + .twin_overlap = twin_overlap_, }; } @@ -903,6 +905,7 @@ void RotationIndexer::ForceResult(const RotationIndexerResult &result) { axis_ = result.axis; unconstrained_ = result.unconstrained; tilt_walk_ = result.tilt_walk; + twin_overlap_ = result.twin_overlap; } bool RotationIndexer::AccumulationFull() const { @@ -912,6 +915,7 @@ bool RotationIndexer::AccumulationFull() const { void RotationIndexer::ForceLattice(const CrystalLattice &lattice) { indexed_lattice = lattice; + twin_overlap_.reset(); const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1(); search_result_ = LatticeSearchResult{ .niggli_class = 0, // Since Niggli class was not searched for, we don't know which one diff --git a/image_analysis/rotation_indexer/RotationIndexer.h b/image_analysis/rotation_indexer/RotationIndexer.h index 49bf1e0c7..fe83d261c 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.h +++ b/image_analysis/rotation_indexer/RotationIndexer.h @@ -37,6 +37,15 @@ struct RotationIndexerResult { bool refused = false; }; std::optional tilt_walk; + // Twinning by reticular merohedry (rugnux's twin-composite step): a second domain whose lattice + // shares the reflections h of this one with h . t divisible by n, t a vector of this lattice (real + // space, in the frame of `lattice`, so it holds through any change of basis). Those reflections + // carry both domains' intensity, and are not integrated. + struct TwinOverlap { + Coord t; + int n = 0; + }; + std::optional twin_overlap; }; // What a memo verification throws where a reused result differs from its recomputation - its own @@ -92,6 +101,7 @@ private: std::vector extra_lattices_; std::shared_ptr unconstrained_; std::optional tilt_walk_; + std::optional twin_overlap_; // only ever forced, see ForceResult IndexerThreadPool &indexer_; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 0f2806303..cfa6392aa 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -404,6 +404,9 @@ struct ProcessResult { // number a few dozen against the merge's tens of thousands, and a bulk statistic therefore cannot // see them at all - which is why the two-pass quality guard counts them separately. int64_t search_merge_axial_reflections = 0; + // The index n of a twin overlap the pass left out of integration (RotationIndexerResult:: + // TwinOverlap), 0 where there was none: its merges then hold (n - 1)/n of the lattice's nodes. + int twin_overlap_index = 0; // How much signal that merge holds: reflections at I/sigma >= 2, and how many reflections it holds // at all. This is what the two-pass quality guard compares the passes on - see RefinedPassIsWorse // for why it is these and not CC1/2. diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index d2d8e1b5e..cc44a332b 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -1553,14 +1553,24 @@ auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass cand.lattice = cand.search_result.conventional; cand.unconstrained.reset(); cand.extra_lattices.clear(); + // The overlap, the reflections on both sub-lattices, holds both domains' intensity: the + // stronger domain's reflections h on the weaker's sub-lattice are those with h . t divisible + // by n, t the weaker sub-lattice's form taken as a vector of the standing cell. + cand.twin_overlap = RotationIndexerResult::TwinOverlap{ + .t = forms[bj][0] * prim.Vec0() + forms[bj][1] * prim.Vec1() + forms[bj][2] * prim.Vec2(), + .n = n}; const int score = CountIndexed(*indexer, cand); const auto uc = cand.lattice.GetUnitCell(); const bool adopt = score > static_cast(validation.size()) / 2; logger.Info("Twin-composite hypothesis: the stronger domain's cell, {}-centred {} {:.2f} {:.2f} " - "{:.2f} {:.1f} {:.1f} {:.1f}, indexes {}/{} validation frames - {}", + "{:.2f} {:.1f} {:.1f} {:.1f}, indexes {}/{} validation frames - {} (the 1/{} of its " + "reflections the other domain overlaps, h.({:.2f},{:.2f},{:.2f}) divisible by {}, is left out " + "of integration)", cand.search_result.centering, gemmi::crystal_system_str(cand.search_result.system), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, score, - static_cast(validation.size()), adopt ? "adopting it" : "keeping the standing cell"); + static_cast(validation.size()), adopt ? "adopting it" : "keeping the standing cell", + n, cand.twin_overlap->t * cand.lattice.Astar(), cand.twin_overlap->t * cand.lattice.Bstar(), + cand.twin_overlap->t * cand.lattice.Cstar(), n); if (adopt) { best = FirstPass{std::move(cand), score, fmt::format("{} (twin domain)", best.name), std::abs(sub.CalcVolume())}; diff --git a/rugnux/RugnuxImagePass.cpp b/rugnux/RugnuxImagePass.cpp index f0d159506..63def847a 100644 --- a/rugnux/RugnuxImagePass.cpp +++ b/rugnux/RugnuxImagePass.cpp @@ -557,6 +557,7 @@ void Rugnux::ImagePass(PipelineLocals &p) { }; result.rotation_lattice_found = true; result.rotation_lattice_type = end_msg.rotation_lattice_type; + result.twin_overlap_index = rot->twin_overlap ? rot->twin_overlap->n : 0; // What the equality this class asserts is worth on a REFINED cell. The indexer refines the // same candidate a second time with nothing held (its pseudo-symmetry guard) and hands that // cell over whenever the constrained one won; it comes back in the primitive reduced diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index 1ad46836d..3aa1ed4e5 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -122,9 +122,10 @@ namespace { // of reflections at a given Miller index the same count on both sides. A pass that re-indexes at // the refined geometry can settle on a centred setting of the same lattice the other pass held // primitive; the cells are then different bases and only the counts taken over all reflections - // survive the change. + // survive the change. Nor is it the same count where one pass left a twin overlap out. bool SameSetting(const ProcessResult &pass1, const ProcessResult &pass2) { if (pass1.consensus_centering != pass2.consensus_centering + || pass1.twin_overlap_index != pass2.twin_overlap_index || !pass1.consensus_cell || !pass2.consensus_cell) return false; const UnitCell &c1 = *pass1.consensus_cell, &c2 = *pass2.consensus_cell; @@ -193,8 +194,14 @@ std::string RefinedPassIsWorse(const ProcessResult &pass1, const ProcessResult & // pass on the crystal's own lattice and the run went back to the supercell, while that merge // held more reflections at <= -2 than at >= 2. On a real superstructure the extra classes hold // signal, and it is counted. - const auto signal = [](const ProcessResult &r) { - return static_cast(r.search_merge_strong_reflections - r.search_merge_negative_reflections); + // A pass that left a twin overlap of index n out of integration merged (n - 1)/n of the + // lattice's nodes, and its counts are taken up to the whole lattice. + const auto whole_lattice = [](const ProcessResult &r) { + return r.twin_overlap_index > 0 ? r.twin_overlap_index / (r.twin_overlap_index - 1.0) : 1.0; + }; + const auto signal = [&](const ProcessResult &r) { + return static_cast(r.search_merge_strong_reflections - r.search_merge_negative_reflections) + * whole_lattice(r); }; // ... and the loss has to be a loss rather than a dilution. The wider disk pulls weaker // reflections into the merge and the error model is then refitted over them, so a pass can hold @@ -208,8 +215,8 @@ std::string RefinedPassIsWorse(const ProcessResult &pass1, const ProcessResult & // the smaller lattice is scaled by the index n, which is an integer - the same lattice, whose // volume two refinements put a percent or two apart, is index 1 and compares as it always did. const auto coverage_does_not_grow = [&] { - double n1 = static_cast(pass1.search_merge_reflections); - double n2 = static_cast(pass2.search_merge_reflections); + double n1 = static_cast(pass1.search_merge_reflections) * whole_lattice(pass1); + double n2 = static_cast(pass2.search_merge_reflections) * whole_lattice(pass2); const double v1 = PrimitiveCellVolume(pass1), v2 = PrimitiveCellVolume(pass2); if (v1 > 0.0 && v2 > 0.0) { if (v1 >= v2) diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 3a99f3c2d..b970745e1 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -395,6 +395,15 @@ TEST_CASE("RefinedPassIsWorse", "[process]") { with_rows.search_merge_axial_reflections = 16; lost_rows.search_merge_axial_reflections = 8; CHECK(!RefinedPassIsWorse(with_rows, lost_rows).empty()); + + // A refined pass that left a twin overlap of index 3 out of integration merged two thirds of the + // nodes: two thirds of the signal is no loss, and its axial rows are not the header pass's rows. + ProcessResult twin = pass(cell, 'P', 34000, 700, 132000); + twin.twin_overlap_index = 3; + twin.search_merge_axial_reflections = 10; + CHECK(RefinedPassIsWorse(with_rows, twin).empty()); + twin.search_merge_strong_reflections = 28000; + CHECK(!RefinedPassIsWorse(with_rows, twin).empty()); } // The beam-centre arbitration between a first pass at the file's centre and one at the measured