diff --git a/common/IndexingSettings.cpp b/common/IndexingSettings.cpp index 6b3c994c8..8d35f57b1 100644 --- a/common/IndexingSettings.cpp +++ b/common/IndexingSettings.cpp @@ -223,6 +223,15 @@ float IndexingSettings::GetRotationIndexingAngularStride_deg() const { return rotation_indexing_angular_stride_deg; } +IndexingSettings & IndexingSettings::LatticeClassDemoted(bool input) { + lattice_class_demoted = input; + return *this; +} + +bool IndexingSettings::GetLatticeClassDemoted() const { + return lattice_class_demoted; +} + bool IndexingSettings::GetIndexIceRings() const { return index_ice_rings; } diff --git a/common/IndexingSettings.h b/common/IndexingSettings.h index 852c951d1..a779fda7b 100644 --- a/common/IndexingSettings.h +++ b/common/IndexingSettings.h @@ -62,6 +62,14 @@ private: bool blocking_behavior = true; bool index_ice_rings = false; + // Refuse the Bravais class the metric was promoted to on a LENGTH equality, and index in the + // class left when that equality is struck out. Not a user setting: rugnux runs one pass with it + // and one without, and keeps the one whose re-integrated held-out positional residual is lower + // (see RunAllPasses). The promotion is decided against a 3 % relative tolerance on a cell the + // spot positions place to ~0.1 %, so on a small cell a genuinely orthorhombic crystal is + // integrated with a = b imposed and nothing downstream can recover from it. + bool lattice_class_demoted = false; + bool enable_rotation_indexing = false; float rotation_indexing_min_angular_range_deg = 20.0; float rotation_indexing_angular_stride_deg = 0.5; @@ -89,6 +97,7 @@ public: IndexingSettings& RotationIndexingAngularStride_deg(float input); IndexingSettings& BlockingBehavior(bool input); IndexingSettings& MaxExtraLattices(int64_t input); + IndexingSettings& LatticeClassDemoted(bool input); [[nodiscard]] int64_t GetViableCellMinSpots() const; [[nodiscard]] IndexingAlgorithmEnum GetAlgorithm() const; @@ -110,4 +119,5 @@ public: [[nodiscard]] float GetRotationIndexingAngularStride_deg() const; [[nodiscard]] bool GetBlockingBehavior() const; [[nodiscard]] int64_t GetMaxExtraLattices() const; + [[nodiscard]] bool GetLatticeClassDemoted() const; }; diff --git a/image_analysis/lattice_search/LatticeSearch.cpp b/image_analysis/lattice_search/LatticeSearch.cpp index b22a331c3..232c5d05f 100644 --- a/image_analysis/lattice_search/LatticeSearch.cpp +++ b/image_analysis/lattice_search/LatticeSearch.cpp @@ -476,6 +476,19 @@ std::optional LatticeSearchForClass(const CrystalLattice &L return SearchCharacters(L, dist_tolerance, angle_tolerance, &only); } + +bool ClassImposesLengthEquality(gemmi::CrystalSystem system) { + switch (system) { + case gemmi::CrystalSystem::Tetragonal: + case gemmi::CrystalSystem::Trigonal: + case gemmi::CrystalSystem::Hexagonal: + case gemmi::CrystalSystem::Cubic: + return true; + default: + return false; + } +} + double MetricDeviation(const LatticeSearchResult &sr) { const UnitCell uc = sr.conventional.GetUnitCell(); double angle = 0.0, length = 0.0; diff --git a/image_analysis/lattice_search/LatticeSearch.h b/image_analysis/lattice_search/LatticeSearch.h index 06f3dead8..396219423 100644 --- a/image_analysis/lattice_search/LatticeSearch.h +++ b/image_analysis/lattice_search/LatticeSearch.h @@ -48,6 +48,14 @@ std::optional LatticeSearchForClass(const CrystalLattice& L double dist_tolerance = LATTICE_SEARCH_DIST_TOLERANCE, double angle_tolerance_deg = LATTICE_SEARCH_ANGLE_TOLERANCE_DEG); +// Whether the conventional metric of this class asserts an equality between cell LENGTHS (a = b, and +// for cubic a = b = c). Those are the promotions the relative tolerance above admits far more widely +// than the spot positions resolve; every other class is carried by its angles alone. The walk run +// with dist_tolerance = 0 grants no such equality, and so returns the class the metric carries on +// its angles alone - the OTHER hypothesis about any of these promotions, which rugnux runs as its +// own pass and judges on what it realises (see RunAllPasses). +bool ClassImposesLengthEquality(gemmi::CrystalSystem system); + // How far the matched lattice's metric sits from the IDEAL metric of its class, as a fraction of the // tolerance that admitted it: 0 is exactly the class, 1 is the edge of what the walk accepts. The // walk computes this and throws it away - a cell 0.02 deg from being C-centred monoclinic and one diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 90845254a..568f21297 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -160,8 +160,18 @@ void RotationIndexer::RunIndexing() { // 60 validation frames and an F-cubic one relabelled trigonal indexed 0 of 60, where the same // frames index at 36/60 and 51/60 without a group. The group is applied where it belongs - to // the scaling and the merge. + const bool demote_class = experiment.GetIndexingSettings().GetLatticeClassDemoted(); auto build_sr = [&](const CrystalLattice &cand) -> LatticeSearchResult { auto ls = LatticeSearch(cand); + // The other hypothesis about a promotion that rests on a = b: walk the characters again + // granting NO length equality at all, and take the class the metric carries on its angles + // alone. Not the same as striking one class out - the tetragonal-P lattice also has a + // C-centred orthorhombic character, which is a different (doubled) cell and rests on the + // very equality in question. The caller runs this as its own pass and keeps whichever pass + // realises the lower held-out positional residual. The hexagonal re-expression below is a + // promotion of the same kind, so it is not applied here. + if (demote_class && ClassImposesLengthEquality(ls.system)) + return LatticeSearch(cand, /*dist_tolerance=*/0.0); if (!IsHexagonalSystem(ls.system) && IsMetricallyHexagonal(ls.primitive_reduced)) { ls.conventional = HexagonalConventional(ls.primitive_reduced); ls.system = gemmi::CrystalSystem::Hexagonal; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 4b7b32632..fa13e6200 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2086,6 +2086,106 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { .PoniRot1_rad(g[3]).PoniRot2_rad(g[4]); }; + // The lattice arms below run passes that exist ONLY to measure - see there - so the + // count of passes the run made has to include them even though they write nothing. + int arm_passes = 0; + + // The metric symmetry, asked of the spots rather than of the tolerance that admitted it. The + // Bravais class is chosen by a walk that reads two axes as EQUAL when they agree to a fixed + // relative tolerance, and the class carrying that equality is then imposed everywhere below: + // the cell is refined with a = b, reflections are predicted from it and integrated at the + // predicted positions. On a small cell the tolerance is far wider than the positions resolve - + // a 5.8 A axis pair 2.4 % apart is read as tetragonal, and the mean the constraint substitutes + // mispredicts the outermost reflections by a good ten pixels against an integration disc of + // four - so a genuinely orthorhombic crystal is integrated as tetragonal and every decision + // downstream, up to the space group, is read off the wreckage. + // + // So the equality is measured rather than assumed. The indexer already refines the same + // candidate a second time with nothing held, and pass 1 reports how far that free refinement + // leaves the two axes apart (lattice_length_split). A relative split of eps displaces a + // reflection at radius r by eps/2 * r pixels, so the question "do these data resolve it" has + // an answer in the units the integration works in: where that displacement at the far corner + // of the detector stays inside the integration disc, imposing the equality moves nothing out + // of its own box and the higher symmetry is kept with no pass run. Where it does not, both + // hypotheses are run - the promoted class, and the class the same walk carries when it is + // granted no length equality at all - and whichever REALISES the lower held-out residual on a + // re-integrated, re-indexed pass is kept, the + // promoted class on a tie, because it is the one that asserts more. Both arms are probe-only + // passes (postrefine_probe_only_): the residual is all the comparison reads and it is in hand + // before the scaling engine is built, so neither arm pays for a merge, a space-group search or + // a set of correction surfaces, and only the class that wins is ever run as a canonical pass. + std::string lattice_arm; // the arm decision, for the pass decision text + if (!cancelled_ && pass1.lattice_length_split && pass1.rotation_lattice_type) { + const double split = *pass1.lattice_length_split; + const double dx = std::max(experiment_.GetBeamX_pxl(), + static_cast(experiment_.GetXPixelsNum()) - experiment_.GetBeamX_pxl()); + const double dy = std::max(experiment_.GetBeamY_pxl(), + static_cast(experiment_.GetYPixelsNum()) - experiment_.GetBeamY_pxl()); + const double r_corner_px = std::hypot(dx, dy); + const double miss_px = 0.5 * split * r_corner_px; + const double disc_px = experiment_.GetBraggIntegrationSettings().GetR1(); + if (miss_px > disc_px) { + logger.Info("Metric symmetry: {} was promoted on a = b, but the free refinement leaves " + "the two axes {:.2f} % apart - {:.1f} px of misprediction at the far corner " + "against an integration disc of {:.1f} px. Both classes are re-integrated " + "and judged on the held-out residual each realises", + gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), + 100.0 * split, miss_px, disc_px); + // The arms are compared on the post-refinement's held-out residual, so they have to + // measure one: a canonical pass only does where the pre-pass moved the geometry. + const bool probe_was = postrefine_probe_; + postrefine_probe_ = true; + postrefine_probe_only_ = true; + const auto promoted_arm = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); + ++arm_passes; + IndexingSettings demoted_settings = experiment_.GetIndexingSettings(); + demoted_settings.LatticeClassDemoted(true); + experiment_.ImportIndexingSettings(demoted_settings); + const auto demoted_arm = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); + ++arm_passes; + postrefine_probe_only_ = false; + postrefine_probe_ = probe_was; + const double realised_promoted = promoted_arm.post_refine ? promoted_arm.post_refine->held_out_before : NAN; + const double realised_demoted = demoted_arm.post_refine ? demoted_arm.post_refine->held_out_before : NAN; + const auto class_text = [](const ProcessResult &r) { + return r.rotation_lattice_type + ? fmt::format("{}-centred {}", r.rotation_lattice_type->centering, + gemmi::crystal_system_str(r.rotation_lattice_type->crystal_system)) + : std::string("no lattice"); + }; + if (!cancelled_ && promoted_arm.post_refine && demoted_arm.post_refine + && HeldOutResidualFell(*promoted_arm.post_refine, *demoted_arm.post_refine)) { + logger.Warning("Metric symmetry: re-integrated as {} the held-out residual is {:.3e} " + "against {:.3e} as {} - the equality is not this crystal's, and the " + "run goes on without it", class_text(demoted_arm), realised_demoted, + realised_promoted, class_text(promoted_arm)); + lattice_arm = fmt::format( + "the metric symmetry was promoted to {} on axes the free refinement leaves " + "{:.2f} % apart; re-integrated as {} the held-out residual is {:.3e} against " + "{:.3e}, so the equality was dropped", class_text(promoted_arm), 100.0 * split, + class_text(demoted_arm), realised_demoted, realised_promoted); + } else if (!cancelled_) { + IndexingSettings promoted_settings = experiment_.GetIndexingSettings(); + promoted_settings.LatticeClassDemoted(false); + experiment_.ImportIndexingSettings(promoted_settings); + logger.Info("Metric symmetry: re-integrated as {} the held-out residual is {:.3e}, " + "which {} does not beat ({:.3e}) - the higher symmetry stands", + class_text(promoted_arm), realised_promoted, class_text(demoted_arm), + realised_demoted); + lattice_arm = fmt::format( + "{} was run as well and realised a held-out residual of {:.3e} against {:.3e} " + "as {}, so the metric symmetry stands", class_text(demoted_arm), + realised_demoted, realised_promoted, class_text(promoted_arm)); + } + } else { + logger.Info("Metric symmetry: {} asserts a = b and the free refinement leaves the two " + "axes {:.2f} % apart - {:.1f} px at the far corner, inside the integration " + "disc of {:.1f} px, so these data do not resolve the difference", + gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system), + 100.0 * split, miss_px, disc_px); + } + } + logger.Info("Rotation two-pass geometry post-refinement: second pass (refined geometry, canonical) -> {}_*", base_prefix); auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); @@ -2205,8 +2305,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { pass1.post_refine->distance_held = false; } pass2.post_refine = pass1.post_refine; - pass2.pass_number = 2 + walk_passes; - pass2.pass_count = 2 + walk_passes; + pass2.pass_number = 2 + arm_passes + walk_passes; + pass2.pass_count = 2 + arm_passes + walk_passes; if (geometry_rounds > 0) pass2.pass_decision = fmt::format( "post-refined geometry walked {} round{} and kept round {}, detector distance {:.3f} mm " @@ -2341,6 +2441,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { pass2.pass_decision = (prepass_detector_geometry_ || prepass_rotation_scale_) ? "post-refined geometry adopted" : "the post-refinement committed no geometry change, so this pass reproduces the first"; + if (!lattice_arm.empty()) + pass2.pass_decision = lattice_arm + "; " + pass2.pass_decision; pass2.geometry_not_converged = geometry_not_converged; return pass2; } @@ -4869,6 +4971,28 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b }; result.rotation_lattice_found = true; result.rotation_lattice_type = end_msg.rotation_lattice_type; + // 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 was seeded from the primitive reduced cell + // and comes back in that basis, so the class's own reindex re-expresses it in the + // conventional one and a and b can be read off it. Measured here, decided in RunAllPasses, + // which is where there is a pass to compare it against. + if (rot->unconstrained && ClassImposesLengthEquality(rot->search_result.system)) { + const UnitCell free_uc = + rot->unconstrained->lattice.Multiply(rot->search_result.reindex).GetUnitCell(); + double split = std::fabs(free_uc.a - free_uc.b) / (0.5 * (free_uc.a + free_uc.b)); + if (rot->search_result.system == gemmi::CrystalSystem::Cubic) + split = std::max(split, std::fabs(free_uc.b - free_uc.c) + / (0.5 * (free_uc.b + free_uc.c))); + if (std::isfinite(split)) { + result.lattice_length_split = split; + logger.Info("Rotation lattice: {}-centred {} asserts a = b; the free refinement of " + "the same candidate leaves them {:.2f} % apart ({:.4f} / {:.4f} A)", + rot->search_result.centering, + gemmi::crystal_system_str(rot->search_result.system), + 100.0 * split, free_uc.a, free_uc.b); + } + } // The consensus cell below is this lattice's conventional cell, so this is the centring // it is written in - the pair the volume comparisons in Run() need. result.consensus_centering = rot->search_result.centering; @@ -5214,8 +5338,20 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // before the engine is built rather than beside it - the gathered observations and their // rocking events are gigabytes on a long axis, and beside the engine's own arrays they set // the pass's memory high-water mark. - if (postrefine_probe_ && !geometry_prepass && is_rotation) + if (postrefine_probe_ && !geometry_prepass && is_rotation) { post_refine_geometry(); + // A pass the lattice arms run to compare two Bravais classes has now produced the only + // thing the comparison reads - the held-out residual of the fit just measured. Stop here + // rather than scale, search, fit surfaces and merge for a result nothing looks at. + if (postrefine_probe_only_) { + result.processing_time_s = std::chrono::duration( + std::chrono::steady_clock::now() - start_time).count(); + logger.Info("Lattice-arm probe: {} images integrated and post-refined in {:.2f} s " + "(the arm reads the held-out residual only, so this pass does not merge)", + result.images_processed, result.processing_time_s); + return result; + } + } if (is_rotation) { if (experiment_.GetRefineRotationWedgeInScaling() || rot_ss.GetRotationWedgeForScaling().has_value()) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 6e606e086..12fa96e32 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -232,6 +232,13 @@ struct ProcessResult { // The metric symmetry the rotation lattice was classified in - the class GeometryRefiner holds the // cell to and the space-group search enumerates under. Empty when no rotation lattice was found. std::optional rotation_lattice_type; + // Where that class asserts a = b: how far apart the FREE (triclinic) refinement of the same + // candidate - which the indexer runs anyway, as its pseudo-symmetry guard - leaves those two + // axes, relative. The class was decided against a fixed relative tolerance on an UNREFINED FFT + // cell; this is the same equality measured on a refined one, and it is what RunAllPasses turns + // into pixels to ask whether this crystal's spot positions resolve the equality at all. Empty + // where the class asserts no length equality, or where the free refinement is not available. + std::optional lattice_length_split; MeanProcessingTime mean_processing_time{}; std::optional written_master_path; std::string merge_statistics_text; // populated when scaling/merging ran @@ -555,6 +562,12 @@ class Rugnux { // whether it has finished walking is a question only reflections integrated after the move can // answer; RunAllPasses turns this on for the passes that answer it. bool postrefine_probe_ = false; + // Whether such a pass exists ONLY for that measurement. The lattice-class arms (see RunAllPasses) + // read the held-out residual the probe returns and nothing else, so a pass run for one of them + // stops as soon as the probe has measured - which is before the scaling engine is built. + // Everything below that point - the merges, the space-group search, the correction surfaces, the + // reports - is work no comparison looks at, and on a long sweep it is four fifths of the pass. + bool postrefine_probe_only_ = false; // The file's detector distance, from before the first pass of the rotation two-pass: a pass whose // distance is still this one asks the post-refinement to test "the header distance is right" // (PostRefineSettings::distance_at_header); a pass that has walked off it does not. diff --git a/tests/LatticeSearchTest.cpp b/tests/LatticeSearchTest.cpp index 56d9094e2..fb16da942 100644 --- a/tests/LatticeSearchTest.cpp +++ b/tests/LatticeSearchTest.cpp @@ -639,3 +639,32 @@ TEST_CASE("SymmetrizeMetric puts a cell onto the metric its group fixes") { CHECK(turned_fixed.Vec0() * turned.Vec0() == Catch::Approx(turned_fixed.Vec0().Length() * turned.Vec0().Length()).epsilon(1e-4)); } + +TEST_CASE("LatticeSearch - a pseudo-tetragonal cell is tetragonal at the walk's tolerance and " + "orthorhombic with no length equality granted") { + // a and b 2.4 % apart: inside LATTICE_SEARCH_DIST_TOLERANCE, far outside what spot positions + // resolve. The two hypotheses rugnux carries forward are exactly these two answers. + const double a = 5.795, b = 5.933, c = 12.362; + CrystalLattice L(Coord(a, 0, 0), Coord(0, b, 0), Coord(0, 0, c)); + + const auto promoted = LatticeSearch(L); + CHECK(promoted.system == gemmi::CrystalSystem::Tetragonal); + CHECK(promoted.centering == 'P'); + CHECK(ClassImposesLengthEquality(promoted.system)); + + const auto free_class = LatticeSearch(L, /*dist_tolerance=*/0.0); + CHECK(free_class.system == gemmi::CrystalSystem::Orthorhombic); + CHECK(free_class.centering == 'P'); + const auto uc = free_class.conventional.GetUnitCell(); + check_uc(uc, a, b, c, 90.0, 90.0, 90.0, 1e-2, 1e-2); +} + +TEST_CASE("ClassImposesLengthEquality names the classes whose metric asserts a = b") { + CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Tetragonal)); + CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Trigonal)); + CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Hexagonal)); + CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Cubic)); + CHECK_FALSE(ClassImposesLengthEquality(gemmi::CrystalSystem::Orthorhombic)); + CHECK_FALSE(ClassImposesLengthEquality(gemmi::CrystalSystem::Monoclinic)); + CHECK_FALSE(ClassImposesLengthEquality(gemmi::CrystalSystem::Triclinic)); +}