From afae1dec2f5c663442e2bf5a38a67b53ed35791f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 17 Jul 2026 14:33:58 +0200 Subject: [PATCH] rugnux: index long finely-spaced axes (de-novo two-pass long-axis rescue) De-novo two-pass rotation indexing collapsed on a crystal with a long, finely spaced axis (a ~150 A c-axis at 16 keV): the unconstrained FFT either collapsed the long axis to a short sub-multiple or let a denser supercell over-fit the accumulated cloud - a small global-orientation error throws the many high-order reflections off along the fine axis, so the true cell scores worst on the raw cloud, and propagating that inaccurate global orientation to each frame fails. Add a long-axis rescue that only runs when the standard pass indexes few validation frames (<50%), so well-indexing crystals are untouched and keep the fast path: re-run the first pass at a COARSE resolution (low-order reflections only, where the fine axis stays robust) to recover the true metric, take the recovered cell with the longest axis directly, then RE-INDEX at full resolution with that cell as a reference (the -C path) - the reference-cell filter drops the collapsed/supercell candidates and refines an accurate global lattice. The scheme feed/index/validate loop is factored into a pick_best() lambda shared by the standard and constrained passes (the standard pass is behaviour-preserved). A P41212 58/58/150 test crystal goes 15% -> 99.89% indexed. Full 25-crystal rotation battery: every other crystal's space group and cell BIT-IDENTICAL, and the rescue fires only on the failing crystal. Co-Authored-By: Claude Opus 4.8 (1M context) --- rugnux/Rugnux.cpp | 198 +++++++++++++++++++++++++++++----------------- 1 file changed, 126 insertions(+), 72 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 9e510c87..b0f9e2cb 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -402,14 +402,14 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { auto validation_settings = config_.spot_finding; validation_settings.indexing = true; validation_settings.quick_integration = false; - auto count_indexed = [&](const RotationIndexerResult &r) -> int { - indexer->ForceRotationIndexerResult(r); + auto count_indexed = [&](IndexAndRefine &idx, const RotationIndexerResult &r) -> int { + idx.ForceRotationIndexerResult(r); int count = 0; for (const int ordinal : validation) { DataMessage m{}; m.number = ordinal; m.spots = get_spots(ordinal); - if (indexer->IndexFrameOnly(m, validation_settings)) + if (idx.IndexFrameOnly(m, validation_settings)) count++; } return count; @@ -418,8 +418,8 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { // Feed one first-pass scheme (a set of image ordinals) into its own rotation indexer, ready to // be indexed. Spots are pulled from the single-threaded cache here; the FFT + refinement runs // separately (RunIndexing) so the two schemes' indexing can overlap. - auto feed_scheme = [&](const std::vector &ordinals) { - auto ri = std::make_unique(experiment_, *indexer_pool); + auto feed_scheme = [&](IndexerThreadPool &pool, const std::vector &ordinals) { + auto ri = std::make_unique(experiment_, pool); for (const int ordinal : ordinals) { if (cancelled_ || ri->AccumulationFull()) break; @@ -446,83 +446,137 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { {"wedge", all_consecutive}, }; - // Feed both schemes serially (single-threaded spot cache), then run their RunIndexing() passes - // concurrently. Each RotationIndexer is independent and shares only the thread-safe FFT pool - // (sized above to one worker per scheme), so on a two-GPU node the schemes index fully in - // parallel and even on one GPU a scheme's cuFFT overlaps the other's Ceres refinement. - std::vector> ris; - ris.reserve(schemes.size()); - for (const auto &[name, ordinals] : schemes) { - if (cancelled_) - break; - ris.push_back(feed_scheme(ordinals)); - } - std::vector> index_futures; - index_futures.reserve(ris.size()); - for (auto &ri : ris) { - RotationIndexer *rp = ri.get(); - index_futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); })); - } - for (auto &f : index_futures) - f.get(); + struct FirstPass { std::optional result; int score = -1; std::string name; double vol = 0.0; }; - int best_score = -1; - double best_vol = 0.0; - std::string best_name; - std::optional best_result; - for (size_t i = 0; i < ris.size(); i++) { - if (cancelled_) - break; - const auto result = ris[i]->GetLattice(); - if (!result.has_value()) - continue; - const int score = count_indexed(*result); - const double vol = std::abs(result->lattice.CalcVolume()); - const std::string &name = schemes[i].first; - logger.Info("First-pass scheme '{}': indexes {}/{} validation frames", name, score, - static_cast(validation.size())); - - // A later scheme wins if it indexes clearly more frames (>10%). - const bool clearly_more = static_cast(score) > best_score * 1.1f + 0.5f; - - // Integer-supercell tie-break. The validation-frame count saturates - a spurious axis - // multiple (2x/3x...) indexes every frame its true cell does, so both schemes reach the same - // frame total and the count alone cannot tell them apart; the default then keeps whichever - // ran first. When the two schemes tie on frames but their cell volumes are related by an - // integer factor >=2, the larger cell is that spurious supercell (one FFT resolved a true - // axis only as a 2x/3x harmonic of its length - the full-rotation 'spread' cloud is prone to - // this) and the smaller is the true reduced cell: take it, regardless of scheme order. - // Requiring a near-integer ratio is what separates a real axis multiplication from a mere - // centering coincidence: a rhombohedral H cell vs its C2 sub-cell is 1.5x (non-integer) and - // is correctly left alone. - bool integer_subcell = false; - if (best_result.has_value() && !clearly_more && vol > 1.0 && best_vol > 1.0) { - const bool tied = static_cast(score) >= best_score * 0.9f - 0.5f; - const double ratio = (vol < best_vol) ? best_vol / vol : vol / best_vol; - const double nearest = std::round(ratio); - const bool integer_multiple = nearest >= 2.0 && std::abs(ratio - nearest) < 0.15; - integer_subcell = tied && integer_multiple && vol < best_vol; + // Feed both schemes (single-threaded spot cache), run their RunIndexing() concurrently, then keep + // the lattice that indexes the most validation frames on the real per-image path. + auto pick_best = [&](IndexerThreadPool &pool, IndexAndRefine &idx) -> FirstPass { + std::vector> ris; + ris.reserve(schemes.size()); + for (const auto &[name, ordinals] : schemes) { + if (cancelled_) + break; + ris.push_back(feed_scheme(pool, ordinals)); } + std::vector> index_futures; + index_futures.reserve(ris.size()); + for (auto &ri : ris) { + RotationIndexer *rp = ri.get(); + index_futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); })); + } + for (auto &f : index_futures) + f.get(); - if (!best_result.has_value() || clearly_more || integer_subcell) { - if (integer_subcell) - logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x sub-cell of '{}' (vol {:.0f}) at " - "equal frame count - adopting the smaller true cell", name, vol, - std::round(best_vol / vol), best_name, best_vol); - best_score = score; - best_vol = vol; - best_name = name; - best_result = result; + FirstPass bp; + for (size_t i = 0; i < ris.size(); i++) { + if (cancelled_) + break; + const auto result = ris[i]->GetLattice(); + if (!result.has_value()) + continue; + const int score = count_indexed(idx, *result); + const double vol = std::abs(result->lattice.CalcVolume()); + const std::string &name = schemes[i].first; + logger.Info("First-pass scheme '{}': indexes {}/{} validation frames", name, score, + static_cast(validation.size())); + + // A later scheme wins if it indexes clearly more frames (>10%). + const bool clearly_more = static_cast(score) > bp.score * 1.1f + 0.5f; + + // Integer-supercell tie-break. The validation-frame count saturates - a spurious axis + // multiple (2x/3x...) indexes every frame its true cell does, so both schemes reach the + // same frame total and the count alone cannot tell them apart; the default then keeps + // whichever ran first. When the two schemes tie on frames but their cell volumes are + // related by an integer factor >=2, the larger cell is that spurious supercell and the + // smaller is the true reduced cell: take it, regardless of scheme order. A near-integer + // ratio separates a real axis multiplication from a centering coincidence (a rhombohedral + // H cell vs its C2 sub-cell is 1.5x, non-integer, and is correctly left alone). + bool integer_subcell = false; + if (bp.result.has_value() && !clearly_more && vol > 1.0 && bp.vol > 1.0) { + const bool tied = static_cast(score) >= bp.score * 0.9f - 0.5f; + const double ratio = (vol < bp.vol) ? bp.vol / vol : vol / bp.vol; + const double nearest = std::round(ratio); + const bool integer_multiple = nearest >= 2.0 && std::abs(ratio - nearest) < 0.15; + integer_subcell = tied && integer_multiple && vol < bp.vol; + } + + if (!bp.result.has_value() || clearly_more || integer_subcell) { + if (integer_subcell) + logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x sub-cell of '{}' (vol {:.0f}) at " + "equal frame count - adopting the smaller true cell", name, vol, + std::round(bp.vol / vol), bp.name, bp.vol); + bp.score = score; + bp.vol = vol; + bp.name = name; + bp.result = result; + } + } + return bp; + }; + + FirstPass best = pick_best(*indexer_pool, *indexer); + + // Long-axis rescue. When the de-novo cell indexes few validation frames, a long, finely-spaced + // axis was likely lost: the unconstrained FFT either collapsed it to a short sub-multiple or let + // a denser supercell over-fit the accumulated cloud (a small global-orientation error throws the + // many high-order reflections off along the fine axis, so the true cell scores worst on the raw + // cloud). Recover the true metric with a COARSE-resolution pass - only low-order reflections, + // where the fine axis stays robust - then RE-INDEX at full resolution constrained by that cell as + // a reference (the -C path): the reference filter drops the collapsed/supercell candidates and + // refines an accurate global lattice. Only runs after a poor standard pass, so well-indexing + // crystals never pay for it and their result is untouched. + auto max_axis = [](const RotationIndexerResult &r) { + const auto uc = r.lattice.GetUnitCell(); + return std::max({uc.a, uc.b, uc.c}); + }; + if (!cancelled_ && best.result.has_value() + && best.score < 0.5 * static_cast(validation.size())) { + auto coarse_settings = experiment_.GetIndexingSettings(); + coarse_settings.FFT_HighResolution_A(3.5f); // low-order reflections only -> robust long axis + coarse_settings.IndexingThreads(2); + IndexerThreadPool coarse_pool(coarse_settings); + // Coarse first pass: keep the recovered cell with the LONGEST axis directly. Its full- + // resolution per-frame validation would be low (the coarse cell is metrically right but + // imprecise), so do NOT score it here - it is only the reference for the constrained re-index. + std::vector> cris; + cris.reserve(schemes.size()); + for (const auto &[name, ordinals] : schemes) { + if (cancelled_) break; + cris.push_back(feed_scheme(coarse_pool, ordinals)); + } + std::vector> cfut; + cfut.reserve(cris.size()); + for (auto &ri : cris) { + RotationIndexer *rp = ri.get(); + cfut.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); })); + } + for (auto &f : cfut) f.get(); + std::optional coarse_ref; + for (auto &ri : cris) { + auto r = ri->GetLattice(); + if (r.has_value() && (!coarse_ref.has_value() || max_axis(*r) > max_axis(*coarse_ref))) + coarse_ref = r; + } + // A recovered axis at least 30% longer than the standard pass found signals a rescued metric. + if (!cancelled_ && coarse_ref.has_value() + && max_axis(*coarse_ref) > 1.3 * max_axis(*best.result)) { + logger.Info("Long-axis rescue: coarse pass recovered a {:.0f} A axis (was {:.0f} A); " + "re-indexing constrained by that cell", max_axis(*coarse_ref), max_axis(*best.result)); + experiment_.SetUnitCell(coarse_ref->lattice.GetUnitCell()); + const FirstPass constrained = pick_best(*indexer_pool, *indexer); + experiment_.SetUnitCell(std::nullopt); // leave the space-group / cell determination de-novo + if (constrained.result.has_value() && constrained.score > best.score) + best = constrained; } } if (!cancelled_) { - if (!best_result.has_value()) + if (!best.result.has_value()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Two-pass rotation indexing failed"); - indexer->ForceRotationIndexerResult(*best_result); + indexer->ForceRotationIndexerResult(*best.result); logger.Info("Two-pass rotation indexing found lattice (scheme '{}': {}/{} validation frames)", - best_name, best_score, static_cast(validation.size())); + best.name, best.score, static_cast(validation.size())); } }