mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 67 Wiki Activity

rugnux: index long finely-spaced axes (de-novo two-pass long-axis rescue)

Browse Source 
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) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-17 20:54:03 +02:00
co-authored by Claude Opus 4.8
parent 799e1e15ec
commit afae1dec2f
1 changed files with 126 additions and 72 deletions
Download Patch File Download Diff File Expand all files Collapse all files
rugnux/Rugnux.cpp
+126 -72
View File
@@ -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<int> &ordinals) {
auto ri = std::make_unique<RotationIndexer>(experiment_, *indexer_pool);
auto feed_scheme = [&](IndexerThreadPool &pool, const std::vector<int> &ordinals) {
auto ri = std::make_unique<RotationIndexer>(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<std::unique_ptr<RotationIndexer>> ris;
ris.reserve(schemes.size());
for (const auto &[name, ordinals] : schemes) {
if (cancelled_)
break;
ris.push_back(feed_scheme(ordinals));
}
std::vector<std::future<void>> 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<RotationIndexerResult> 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<RotationIndexerResult> 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<int>(validation.size()));
// A later scheme wins if it indexes clearly more frames (>10%).
const bool clearly_more = static_cast<float>(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<float>(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<std::unique_ptr<RotationIndexer>> ris;
ris.reserve(schemes.size());
for (const auto &[name, ordinals] : schemes) {
if (cancelled_)
break;
ris.push_back(feed_scheme(pool, ordinals));
}
std::vector<std::future<void>> 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<int>(validation.size()));
// A later scheme wins if it indexes clearly more frames (>10%).
const bool clearly_more = static_cast<float>(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<float>(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<double>(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<std::unique_ptr<RotationIndexer>> cris;
cris.reserve(schemes.size());
for (const auto &[name, ordinals] : schemes) {
if (cancelled_) break;
cris.push_back(feed_scheme(coarse_pool, ordinals));
}
std::vector<std::future<void>> 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<RotationIndexerResult> 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<int>(validation.size()));
best.name, best.score, static_cast<int>(validation.size()));
}
}