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

Sweep ends, overloads and ice flags: one walk over the reflections instead of two

Browse Source 
Before the merge, the integrated reflections were walked twice - once per frame for what the
blank sweep ends are decided on, once more for the overload counts and the ice flags - and
each walk reads the whole reflection set from memory. The frame's sweep-end tally is now a
function of its own (TallySweepEnd), and the merge step makes it in the same per-frame walk
that counts the overloads and flags the ice, then adds the counts up over the frames still in,
so the trimmed and left-out frames count nothing, as before. Per frame the sums are the same
sums in the same order; the totals are integers. TrimBlankSweepEnds keeps its one-argument
form for the callers that walk the frames for nothing else.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
This commit is contained in:
leonarski_fandClaude Fable 5.1 committed 2026-10-26 20:58:10 +01:00
1 parent 093ac1612c
commit dfb1536fe2
3 files changed
+70 -35
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

image_analysis/scale_merge/RotationScaleMerge.cpp
+28 -13
View File
@@ -408,23 +408,38 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment,
// the scaling and the sweep ledger to judge as they judge any other. Only from the ends inward: a frame
// inside the sweep without beam is news, and is left to the ledger too. A frame with nothing integrated
// on it (missing from disk, not processed) is walked past and not counted.
std::pair<int, int> TrimBlankSweepEnds(std::vector<IntegrationOutcome> &outcomes, size_t nthreads) {
namespace {
constexpr double BLANK_TAIL_SIGMA = 3.0;
}
SweepEndTally TallySweepEnd(const IntegrationOutcome &outcome) {
SweepEndTally t;
for (const auto &r : outcome.reflections) {
if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) {
t.upper += r.I > BLANK_TAIL_SIGMA * r.sigma;
t.lower += r.I < -BLANK_TAIL_SIGMA * r.sigma;
}
if (std::isfinite(r.bkg)) { t.bkg_sum += r.bkg; ++t.bkg_n; }
}
return t;
}
std::pair<int, int> TrimBlankSweepEnds(std::vector<IntegrationOutcome> &outcomes, size_t nthreads) {
std::vector<SweepEndTally> tally(outcomes.size());
ParallelFor(static_cast<int>(outcomes.size()), nthreads, [&](int f) { tally[f] = TallySweepEnd(outcomes[f]); });
return TrimBlankSweepEnds(outcomes, tally);
}
std::pair<int, int> TrimBlankSweepEnds(std::vector<IntegrationOutcome> &outcomes,
const std::vector<SweepEndTally> &tally) {
const int n = static_cast<int>(outcomes.size());
std::vector<int64_t> upper(n, 0), lower(n, 0);
std::vector<double> bkg(n, 0.0);
ParallelFor(n, nthreads, [&](int f) {
double bkg_sum = 0.0;
int64_t bkg_n = 0;
for (const auto &r : outcomes[f].reflections) {
if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) {
upper[f] += r.I > BLANK_TAIL_SIGMA * r.sigma;
lower[f] += r.I < -BLANK_TAIL_SIGMA * r.sigma;
}
if (std::isfinite(r.bkg)) { bkg_sum += r.bkg; ++bkg_n; }
}
if (bkg_n > 0) bkg[f] = bkg_sum / bkg_n;
});
for (int f = 0; f < n; ++f) {
upper[f] = tally[f].upper;
lower[f] = tally[f].lower;
if (tally[f].bkg_n > 0) bkg[f] = tally[f].bkg_sum / tally[f].bkg_n;
}
std::vector<double> excess, tails, background;
for (int f = 0; f < n; ++f)
if (!outcomes[f].reflections.empty()) {
image_analysis/scale_merge/RotationScaleMerge.h
+11 -1
View File
@@ -67,7 +67,17 @@ template <class T> using FillVector = std::vector<T, DefaultInitAllocator<T>>;
// Stills use the per-image ScaleOnTheFly (fixed partiality) instead.
//
// Before any of it: TrimBlankSweepEnds clears the reflections of the frames at either end of the sweep
// that the beam never reached, and returns how many it cleared at the start and at the end.
// that the beam never reached, and returns how many it cleared at the start and at the end. It decides
// on each frame's SweepEndTally, one walk over the frame's reflections in their order; a caller walking
// them anyway for something else hands the tallies over, the other form makes them on all threads.
struct SweepEndTally {
int64_t upper = 0, lower = 0; // reflections standing above / below zero by more than the tail sigma
double bkg_sum = 0.0;
int64_t bkg_n = 0;
};
SweepEndTally TallySweepEnd(const IntegrationOutcome &outcome);
std::pair<int, int> TrimBlankSweepEnds(std::vector<IntegrationOutcome> &outcomes,
const std::vector<SweepEndTally> &tally);
std::pair<int, int> TrimBlankSweepEnds(std::vector<IntegrationOutcome> &outcomes, size_t nthreads);
class RotationScaleMerge {
rugnux/RugnuxScaleMerge.cpp
+31 -21
View File
@@ -507,8 +507,31 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
// and what a model's low-resolution terms rest on. The ice flags are set in the same sweep. Integer
// counts, one image per task: the totals do not depend on the split.
TimingMark ends_mark(logger, "sweep ends, overloads, ice rings");
// One walk over the reflections, per frame, for all three: what the sweep ends are decided on, the
// overload counts and the ice flags. The counts are added up below over the frames still in.
auto &outcomes = indexer->GetIntegrationOutcome();
struct FrameCounts { int64_t all = 0, low = 0, over = 0, over_low = 0, on_ice = 0; };
std::vector<SweepEndTally> tally(outcomes.size());
std::vector<FrameCounts> counts(outcomes.size());
ParallelFor(static_cast<int>(outcomes.size()), config_.nthreads, [&](int o) {
tally[o] = TallySweepEnd(outcomes[o]);
FrameCounts &c = counts[o];
for (auto &r : outcomes[o].reflections) {
const bool low_res = r.d >= OVERLOAD_LOW_RES_A;
++c.all;
c.low += low_res;
if (r.overloaded) {
++c.over;
c.over_low += low_res;
}
if (flag_ice) {
r.on_ice_ring = IsOnIceRing(r.d, ice_width);
c.on_ice += r.on_ice_ring;
}
}
});
if (experiment_.IsRotationIndexing()) {
const auto [at_start, at_end] = TrimBlankSweepEnds(indexer->GetIntegrationOutcome(), config_.nthreads);
const auto [at_start, at_end] = TrimBlankSweepEnds(outcomes, tally);
result.blank_start_frames = at_start;
result.blank_end_frames = at_end;
if (at_start + at_end > 0)
@@ -518,7 +541,6 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
// Stills given a cell: only the frames whose refined cell belongs to its population are merged
// (see GivenCellPopulation) - a second crystal form, fitted to the same cell, is left out.
if (!experiment_.IsRotationIndexing() && stills_user_cell_) {
auto &outcomes = indexer->GetIntegrationOutcome();
const double given = CrystalLattice(*stills_user_cell_).CalcVolume();
std::vector<double> rel;
for (const auto &o : outcomes)
@@ -550,25 +572,13 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
}
size_t ice_flagged = 0, ice_total = 0;
{
auto &outcomes = indexer->GetIntegrationOutcome();
std::atomic<int64_t> n_all{0}, n_low{0}, n_over{0}, n_over_low{0}, n_ice{0};
ParallelFor(static_cast<int>(outcomes.size()), config_.nthreads, [&](int o) {
int64_t all = 0, low = 0, over = 0, over_low = 0, on_ice = 0;
for (auto &r : outcomes[o].reflections) {
const bool low_res = r.d >= OVERLOAD_LOW_RES_A;
++all;
low += low_res;
if (r.overloaded) {
++over;
over_low += low_res;
}
if (flag_ice) {
r.on_ice_ring = IsOnIceRing(r.d, ice_width);
on_ice += r.on_ice_ring;
}
}
n_all += all; n_low += low; n_over += over; n_over_low += over_low; n_ice += on_ice;
});
int64_t n_all = 0, n_low = 0, n_over = 0, n_over_low = 0, n_ice = 0;
for (size_t o = 0; o < outcomes.size(); ++o) {
if (outcomes[o].reflections.empty()) // a frame left out above counts nothing
continue;
const FrameCounts &c = counts[o];
n_all += c.all; n_low += c.low; n_over += c.over; n_over_low += c.over_low; n_ice += c.on_ice;
}
result.integrated_reflections += n_all;
result.integrated_reflections_low_res += n_low;
result.overloaded_reflections += n_over;