From 1bd39b0563ec3e35dda5219e6fc3fea92aa7f589 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 3 Oct 2026 21:00:25 +0200 Subject: [PATCH 01/13] Space-group search: cut the P1 search merge where a monotone fit of drops under 1 The search merge was cut at the first thin shell whose mean fell under 1. On a merge whose is flat near 1 - its ISa has collapsed, as on a strongly absorbing garnet - the first single shell to dip is decided by noise: the same crystal cut at 1.19 A without -march and 0.85 A with -march=x86-64-v3, and the coarser cut left each glide zone fewer than 20 absences, so Ia-3d became unjudgeable and I4(1)32 was adopted. The cut now comes from a non-increasing (pool-adjacent-violators) fit of the shell means, which moves only as much as its input does; that garnet's profile never falls under 1, so its search sees the full range and adopts Ia-3d under both builds with identical candidate tables. On monotone profiles nothing changes: myob/cytc/thau x10sa keep their cuts and p.mtz byte for byte. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- rugnux/Rugnux.cpp | 34 +++++++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 94a9325af..aaaa11d02 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -7133,11 +7133,39 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b [](const auto &a, const auto &b) { return a.first > b.first; }); // low -> high res const int bins = std::clamp(static_cast(rs.size() / 100), 1, 40); const size_t per = (rs.size() + bins - 1) / static_cast(bins); - for (size_t b = 0; b * per < rs.size(); ++b) { - const size_t lo = b * per, hi = std::min(rs.size(), lo + per); + // The shell means, then their best NON-INCREASING fit (pool adjacent violators, weighted + // by shell size). The cut is where the FIT drops under 1, not where the first shell + // does: on a merge whose is flat near 1 - its ISa has collapsed - the first + // single shell to dip under 1 is a coin toss that moved the cut from 0.85 to 1.19 A + // between two compiler flags on one crystal, and the coarser cut left its glide zones + // too few absences to be judged. The fit moves by as much as its input does. + std::vector mean, weight; + for (size_t lo = 0; lo < rs.size(); lo += per) { + const size_t hi = std::min(rs.size(), lo + per); double sum = 0.0; for (size_t j = lo; j < hi; ++j) sum += rs[j].second; - if (sum / static_cast(hi - lo) < 1.0) { + mean.push_back(sum / static_cast(hi - lo)); + weight.push_back(static_cast(hi - lo)); + } + std::vector block_mean, block_weight; + std::vector block_size; + for (size_t b = 0; b < mean.size(); ++b) { + block_mean.push_back(mean[b]); block_weight.push_back(weight[b]); block_size.push_back(1); + while (block_mean.size() > 1 && block_mean[block_mean.size() - 2] < block_mean.back()) { + const size_t k = block_mean.size() - 1; + const double w = block_weight[k - 1] + block_weight[k]; + block_mean[k - 1] = (block_mean[k - 1] * block_weight[k - 1] + block_mean[k] * block_weight[k]) / w; + block_weight[k - 1] = w; + block_size[k - 1] += block_size[k]; + block_mean.pop_back(); block_weight.pop_back(); block_size.pop_back(); + } + } + std::vector fit; + for (size_t k = 0; k < block_mean.size(); ++k) + fit.insert(fit.end(), block_size[k], block_mean[k]); + for (size_t b = 0; b < fit.size(); ++b) { + const size_t lo = b * per; + if (fit[b] < 1.0) { // Cut the noise-dominated high-res shells. When even the lowest-res shell fails, // keep that shell alone rather than abandoning the cut. The bound is absolute // while the merged I/sigma it tests saturates at the merge's own ISa (times the From 445804fa43211064c0cc10932b448638030984ef Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 3 Oct 2026 21:33:58 +0200 Subject: [PATCH 02/13] rugnux: carry pass 1's lattice into the refined pass; revert to pass 1 itself Two knife-edges let a sub-0.05 px change of the measured beam centre turn a P1 crystal (9qw8) from a pass into a merge with R_meas 60000% / ISa 0: - The refined-geometry pass re-indexes de novo and fell back to pass 1's lattice only when the re-index scored under the 1/6 validation floor. The long-axis rescue lifted a related C-centred 303 A cell to 41/60 - over the floor - while pass 1's lattice indexes 55/60 at the same geometry. Pass 1's lattice is now scored as a hypothesis of its own whenever the re-index found a DIFFERENT lattice (class + primitive volume within 2 %), and is integrated when it indexes more validation frames. The same lattice found again is kept as the re-index refined it, so ordinary crystals are untouched. - The two-pass quality guard's "going back to the header geometry" re-ran pass 1 de novo, a hypothesis nobody had judged: at the centre where pass 1 indexed 56/60 it found 5/60, flipped the axis sign and shipped garbage. It now forces pass 1's whole indexing result, which is what the guard preferred. 9qw8: P1 1.71 A on -march=x86-64-v3 and on no-march GPU builds (both failed before, one catastrophically). myob/cytc/thau p.mtz md5 unchanged (GPU). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- rugnux/Rugnux.cpp | 107 ++++++++++++++++++++++++++++------------------ 1 file changed, 66 insertions(+), 41 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 94a9325af..74e05a89d 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -3097,6 +3097,13 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { // This pass is the answer whatever it measures - the guard has had its one chance - // so it must not conclude from its own search merge that it is going to be re-run. quality_guard_pass1_.reset(); + // What the guard preferred is pass 1 - its lattice, orientation and axis at this very + // geometry - so that is what is re-integrated, not a fresh de-novo index. A de-novo + // re-index here is a new hypothesis nobody judged: measured, at the same centre where + // pass 1 indexed 56/60 validation frames it found 5/60, took the opposite axis sign and + // shipped a merge with R_meas in the hundreds of percent. + if (prepass_result_) + force_rotation_result_ = *prepass_result_; rerun_if_starved_ = true; auto redo = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); rerun_if_starved_ = false; @@ -3105,7 +3112,9 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { redo.pass_count = pass2.pass_count + 1; redo.pass_decision = fmt::format( "header geometry re-adopted: the post-refined pass was worse ({})", worse); + // Still forced: the re-run at the fixed radius is the same pass. rerun_if_starved(redo, redo.pass_decision); + force_rotation_result_.reset(); return redo; } if (pass2.pass_decision.empty()) @@ -5475,55 +5484,71 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // frame in six, several times below the weakest crystal that still merges. It is checked AFTER // the long-axis rescue, so a metric that rescue recovers is never rejected on its pre-rescue // score. - if (best.score < static_cast(validation.size()) / 6) { - // The second (refined-geometry) pass re-indexes DE NOVO, and the post-refined geometry - - // fitted to the lattice pass 1 already found - can tip the blind FFT onto an axis harmonic - // of a long cell that indexes too few frames to clear this floor (measured: a tripled - // ~100 A c-axis at 7/60). Pass 1 found and integrated the true cell at the header geometry, - // so prefer its whole result rather than abort the run on a re-index only the refined - // geometry made fail. This is the floor-side companion of the supercell-collapse guard - // below: that one only reaches a pass-2 lattice that DID clear the floor, whereas here the - // harmonic fell BELOW it and the throw would run before the guard could ever see it. - if (!geometry_prepass && prepass_result_.has_value()) { - const auto &pc = prepass_result_->search_result.conventional.GetUnitCell(); + // The second (refined-geometry) pass re-indexes DE NOVO, and the post-refined geometry - fitted + // to the lattice pass 1 already found - can tip the blind FFT, or the long-axis rescue above, + // onto another lattice: an axis harmonic of a long cell (measured: a tripled ~100 A c-axis at + // 7/60), or a related centred cell that indexes fewer frames than pass 1's lattice does at + // this very geometry (measured: a C-centred 303 A cell at 41/60 against pass 1's triclinic + // cell at 60/60). Pass 1's lattice is therefore carried here as a hypothesis of its own and + // the two are judged on the same validation frames: where the re-index found a different + // lattice and indexes fewer frames than pass 1's, or too few to be this crystal's lattice at + // all, pass 1's whole result is integrated instead. The same lattice found again is kept as + // the re-index refined it at this geometry. This is the floor-side companion of the + // supercell-collapse guard in RunAllPasses, which only sees a pass-2 lattice that is LARGER. + if (!geometry_prepass && prepass_result_.has_value()) { + RotationIndexerResult pass1_lattice = *prepass_result_; + // Pass 1's result carries pass 1's goniometer; forcing it whole would put the + // uncorrected angles back (the same repair as the supercell re-run in RunAllPasses). + if (prepass_rotation_scale_ && pass1_lattice.axis) + pass1_lattice.axis = ScaleRotation(*pass1_lattice.axis, *prepass_rotation_scale_); + // Pass 1's lattice was FITTED at pass 1's detector distance, and this pass + // integrates at the post-refined one. Re-scoring it here is not the same as + // re-fitting it: a real-space cell is measured against the distance the spots were + // seen at, so carrying it across a distance change scales the whole cell by the + // ratio of the two. Measured on a 545 A axis over 225 -> 226 mm: 0.4 %, exactly + // 225/226, shipped as the run's answer. Scale it to the distance it will be used + // at; the orientation is unaffected. + const float d_fit = prepass_result_->geom.GetDetectorDistance_mm(); + const float d_use = experiment_.GetDetectorDistance_mm(); + const bool rescaled = d_fit > 0 && d_use > 0 && d_fit != d_use; + if (rescaled) { + const float f = d_use / d_fit; + const auto &l = pass1_lattice.lattice; + pass1_lattice.lattice = CrystalLattice(l.Vec0() * f, l.Vec1() * f, l.Vec2() * f); + } + const double pass1_vol = std::abs(pass1_lattice.lattice + .ToPrimitive(pass1_lattice.search_result.centering).CalcVolume()); + // Same lattice as in the beam-centre check: same class and primitive volumes within 2 %. + const bool same_lattice = best.result.has_value() + && best.result->search_result.centering == pass1_lattice.search_result.centering + && best.result->search_result.system == pass1_lattice.search_result.system + && std::abs(best.vol - pass1_vol) <= 0.02 * pass1_vol; + const int floor = static_cast(validation.size()) / 6; + const int pass1_score = (best.score < floor || !same_lattice) + ? count_indexed(*indexer, pass1_lattice) : -1; + if (best.score < floor || pass1_score > best.score) { + const auto &pc = pass1_lattice.search_result.conventional.GetUnitCell(); const auto &bc = best.result->search_result.conventional.GetUnitCell(); - logger.Warning("Two-pass: the refined-geometry re-index indexes only {}/{} validation " - "frames ({}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}) - too " - "few for this crystal's lattice; integrating with pass-1's lattice " - "instead ({}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f})", + logger.Warning("Two-pass: the refined-geometry re-index indexes {}/{} validation frames " + "({}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}), pass-1's lattice " + "{}/{} at this geometry ({}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} " + "{:.2f}) - integrating with pass-1's lattice", best.score, static_cast(validation.size()), best.result->search_result.centering, gemmi::crystal_system_str(best.result->search_result.system), bc.a, bc.b, bc.c, bc.alpha, bc.beta, bc.gamma, - prepass_result_->search_result.centering, - gemmi::crystal_system_str(prepass_result_->search_result.system), + pass1_score, static_cast(validation.size()), + pass1_lattice.search_result.centering, + gemmi::crystal_system_str(pass1_lattice.search_result.system), pc.a, pc.b, pc.c, pc.alpha, pc.beta, pc.gamma); - best.result = *prepass_result_; - // Pass 1's result carries pass 1's goniometer; forcing it whole would put the - // uncorrected angles back (the same repair as the supercell re-run in RunAllPasses). - if (prepass_rotation_scale_ && best.result->axis) - best.result->axis = ScaleRotation(*best.result->axis, *prepass_rotation_scale_); - // Pass 1's lattice was FITTED at pass 1's detector distance, and this pass - // integrates at the post-refined one. Re-scoring it here is not the same as - // re-fitting it: a real-space cell is measured against the distance the spots were - // seen at, so carrying it across a distance change scales the whole cell by the - // ratio of the two. Measured on a 545 A axis over 225 -> 226 mm: 0.4 %, exactly - // 225/226, shipped as the run's answer. Scale it to the distance it will be used - // at; the orientation is unaffected. - const float d_fit = prepass_result_->geom.GetDetectorDistance_mm(); - const float d_use = experiment_.GetDetectorDistance_mm(); - if (d_fit > 0 && d_use > 0 && d_fit != d_use) { - const float f = d_use / d_fit; - const auto &l = best.result->lattice; - best.result->lattice = CrystalLattice(l.Vec0() * f, l.Vec1() * f, l.Vec2() * f); + if (rescaled) logger.Info("Two-pass: pass-1's cell was fitted at {:.3f} mm and is used at " "{:.3f} mm - scaled by {:.6f} to the distance it is integrated at", - d_fit, d_use, f); - } - best.score = count_indexed(*indexer, *best.result); // re-score at the refined geometry - best.vol = std::abs(best.result->lattice - .ToPrimitive(best.result->search_result.centering).CalcVolume()); - best.name = "pass-1 lattice (refined-geometry re-index too sparse)"; + d_fit, d_use, d_use / d_fit); + best.result = std::move(pass1_lattice); + best.score = pass1_score; + best.vol = pass1_vol; + best.name = "pass-1 lattice (refined-geometry re-index found less)"; } } From da5ef0c5032031a329943bab69714c1575146d64 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 3 Oct 2026 21:40:42 +0200 Subject: [PATCH 03/13] SearchSpaceGroup: offer every setting of the chosen rotation set, not only on a near-exact cell The non-reference settings of the adopted point group were offered only when the cell hosted their rotations to within ~0.1 deg, while the reference settings - which hold the very same rotations - were offered without asking. A cell refined free after integration a few tenths of a degree off 90 therefore lost every setting but the reference one. On an orthorhombic set whose measured screws lie on a and c and whose b row was never recorded, that left P2(1)2(1)2(1) as the only candidate covering both screws, and it was reported as determined. With P 21 2 21 offered, the two tie, the b-axis screw is reported as undetermined, and the model check uses the setting the data describe. The point-group stage still asks the cell whether a rotation set it adds is hosted; only the setting enumeration within an already chosen set stops asking. Validation: myob/cytc/thau x10sa p.mtz byte-identical (GPU); 7mzt fail -> unscored (b screw undetermined, P 21 21 21 or P 21 2 21); 5cc8 unchanged; [SearchSpaceGroup] 23 cases pass incl. a new section with a cell 0.2-0.3 deg off 90. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- image_analysis/scale_merge/SearchSpaceGroup.cpp | 14 +++++++++----- tests/SearchSpaceGroupTest.cpp | 11 +++++++++++ 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 3b99f6e3e..59eb289eb 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -1900,9 +1900,14 @@ SearchSpaceGroupResult SearchSpaceGroup( // the point group is one only a non-reference setting carries (Stage A's second pass), since // otherwise there is no candidate at all and the group would be lost after being found. These are // alternative namings at the same order, so nothing here can promote the point group; what they - // add is a screw or a centering on the axis the data show it on. Two refusals bound them: a - // setting whose axes the cell does not have is not offered, and one predicting exactly the - // absences a candidate already offered predicts is the same hypothesis under another name. + // add is a screw or a centering on the axis the data show it on. One refusal bounds them: a setting + // predicting exactly the absences a candidate already offered predicts is the same hypothesis + // under another name. The cell is NOT asked again here: a setting of the chosen rotation set holds + // the very rotations the reference candidates hold, so it fits the cell exactly as well as they + // do - and they are offered without asking. Asking only the non-reference ones turned a free- + // refined cell a few tenths of a degree off 90 into a refusal of every setting but the reference + // one, leaving the group with the measured screws on a and c (P 2_1 2 2_1) unoffered and P2_12_12_1 + // the only candidate that covered both. if (opt.cell.has_value() && (opt.enumerate_all_settings || (opt.enumerate_all_rotation_sets && sg_cands.empty()))) { std::vector> signatures; @@ -1910,8 +1915,7 @@ SearchSpaceGroupResult SearchSpaceGroup( signatures.push_back(AbsenceSignature(*c)); for (const auto& sg : gemmi::spacegroup_tables::main) { if (!sg.is_sohncke() || sg.is_reference_setting() || - RotationSetOf(sg) != best_pg->rotation_set || - !CellHostsRotations(*opt.cell, sg.operations())) + RotationSetOf(sg) != best_pg->rotation_set) continue; auto sig = AbsenceSignature(sg); if (std::find(signatures.begin(), signatures.end(), sig) != signatures.end()) diff --git a/tests/SearchSpaceGroupTest.cpp b/tests/SearchSpaceGroupTest.cpp index 88123bdd3..59b912021 100644 --- a/tests/SearchSpaceGroupTest.cpp +++ b/tests/SearchSpaceGroupTest.cpp @@ -444,6 +444,17 @@ TEST_CASE("SearchSpaceGroup names an orthorhombic screw pair on the axes it lies REQUIRE(result.best_space_group.has_value()); CHECK(result.best_space_group->xhm() == "P 2 21 21"); } + + // A cell refined free after integration is a few tenths of a degree off 90. The reference + // candidates hold the same rotations and are offered on it, so the other settings must be too. + SECTION("a free cell a few tenths off 90 still offers it") { + opt.cell = gemmi::UnitCell(40.0, 50.0, 60.0, 90.32, 90.19, 90.19); + opt.enumerate_all_settings = true; + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->xhm() == "P 2 21 21"); + } } // A screw on a row the sweep never recorded is not a group the data refused, it is a question From 42441ac3283ccc359f231dc588ec9cba595a8646 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 3 Oct 2026 21:40:51 +0200 Subject: [PATCH 04/13] Twinning report: warn when a promotion stands on intensities below a perfect twin's; flag compressed twin-immune zone controls A promoted point group whose <|L|> reads below 0.375 (NOT_READABLE) previously raised no TWINNING warning at all; a twinned subgroup whose law the promotion absorbed predicts the same data, so the adopted group is not confirmed. Warn. The twin-immune zone control reading more compressed than a perfect twin's acentric population (0.541) is something no twin fraction produces (overlap or neighbour correlation); genuine symmetry then reads acentric in its zones too (measured: a genuine 622 with its control at 0.528 read its 2-folds at -950 to -2640 nats). Such zones are now marked ambiguous in the text and the zone decision line. Report-only: no decision and no output file but the report changes. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- image_analysis/scale_merge/TwinningAnalysis.cpp | 10 ++++++++++ image_analysis/scale_merge/TwinningAnalysis.h | 9 +++++++++ rugnux/ResultReport.cpp | 11 +++++++++++ rugnux/Rugnux.cpp | 4 +++- 4 files changed, 33 insertions(+), 1 deletion(-) diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index 78ffeb15e..3f7ba6f6e 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -596,6 +596,10 @@ namespace { result.control.evidence_nats - result.control.n * result.control_excess_per_reflection; for (auto& z : result.zones) z.calibrated_evidence_nats = z.evidence_nats - z.n * result.control_excess_per_reflection; + constexpr double PERFECT_TWIN_ACENTRIC_MEAN_ABS_E2_MINUS_1 = 0.541; + result.zones_ambiguous = result.control.n > 0 + && result.control.mean_abs_e2_minus_1 + 3.0 * result.control.standard_error + < PERFECT_TWIN_ACENTRIC_MEAN_ABS_E2_MINUS_1; } } @@ -862,5 +866,11 @@ std::string TwinImmuneZonesToText(const TwinImmuneZoneResult& result) { for (const auto& z : result.zones) row(z); row(result.control); + if (result.zones_ambiguous) + os << " The control reads more compressed than a perfect twin's acentric reflections (0.541), which\n" + " no twin fraction can do: something else - overlapping spots, neighbour correlation - averages\n" + " each reflection with unrelated ones, and it compresses the zones as well. Genuine symmetry then\n" + " reads acentric in its zones too, so here they do not separate a twinned subgroup from the\n" + " adopted group, and no verdict is read from them either way.\n"; return os.str(); } diff --git a/image_analysis/scale_merge/TwinningAnalysis.h b/image_analysis/scale_merge/TwinningAnalysis.h index cd831d2fb..688daf948 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.h +++ b/image_analysis/scale_merge/TwinningAnalysis.h @@ -143,6 +143,15 @@ struct TwinImmuneZoneResult { // nats were the normalisation's, and the verdict rescued a twin law. Calibrated, that zone reads // -115 nats; the zones of genuine promotions keep +140 nats (a 90 deg tetragonal sweep) to +2000. double control_excess_per_reflection = 0.0; + // The control reads more compressed than a PERFECT twin's acentric population (<|E^2-1|> 0.541, by + // more than three standard errors). Twinning cannot do that at any fraction, so something else + // averages each reflection with unrelated ones - overlapping spots of a long cell, neighbour + // correlation - and it compresses the zones as much as the control: genuine symmetry then reads + // acentric in its zones. Measured on the open battery: a 622 crystal with its control at 0.528 + // read its genuine 2-folds at -950 to -2640 nats, and a twinned 4/m crystal promoted to 4/mmm, its + // control at 0.454, read its added operators at -226 to -616. The zones are reported, but no + // verdict is read from them either way. + bool zones_ambiguous = false; double anisotropy_delta_b_A2 = 0.0; // the anisotropy taken out of the intensities before normalising double d_max_A = 0.0; // resolution range of the shells read double d_min_A = 0.0; diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 189295c7b..7cd5edeae 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1663,6 +1663,17 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Warn(doc, PathologyCode::TWINNING, fmt::format( "Twinning is indicated ({}) - refinement against the merged data needs a twin law", TwinningVerdictLine(tw))); + else if (abnormal && tw.laue_class_was_chosen_by_promotion && tw.mean_abs_l < 0.375) + // Below a perfect twin's 0.375 the L-test says nothing about the twin fraction, but it does + // say the intensities are at least as twin-like as a perfect twin's - and a perfect twin of + // a subgroup, its law absorbed by the promotion, predicts exactly the data the promoted + // group does. The promotion is then unconfirmed, however clean its operators look. + Warn(doc, PathologyCode::TWINNING, fmt::format( + "The point group was promoted although the intensities read at least as twin-like as a " + "perfect twin's (<|L|> = {:.3f}, below 0.375): a twin of a subgroup whose law the promotion " + "absorbed predicts the same data, so the adopted group is not confirmed - refine in the " + "subgroup with that twin law too (p_P1.mtz carries the same observations)", + tw.mean_abs_l)); } // ---- radiation damage diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 94a9325af..917f12968 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2103,7 +2103,9 @@ namespace { zone.standard_error, z.control.mean_abs_e2_minus_1, z.control_excess_per_reflection, zone.evidence_nats, zone.calibrated_evidence_nats, z.d_max_A, z.d_min_A, z.tncs_normalised ? ", normalised per pseudo-translation class" : "", - verdict < 0 ? "acentric, the added operators are a twin law or a pseudo-symmetry" + verdict < 0 && z.zones_ambiguous + ? "acentric, but the control is compressed beyond a perfect twin, so genuine symmetry would read so too" + : verdict < 0 ? "acentric, the added operators are a twin law or a pseudo-symmetry" : verdict > 0 ? "centric, the added operators are real" : "undecided"); // Recorded only where the zone DECIDED something. An undecided read decides nothing - the // rule stays as it was - and recording it had the report announce "a twin gate refused this From 769a4dd8debc2c5732e700156005a77070a9ed35 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 3 Oct 2026 22:35:00 +0200 Subject: [PATCH 05/13] rugnux: pass-1 lattice hypothesis compares lattices, not classes The first version judged "the same lattice" by class and centring, so an F-centred cubic re-index (57/60) and pass 1's triclinic primitive of the very same lattice (60/60) counted as different, and the unconstrained primitive won the frame count - 6oel lost its cubic setting (R_meas 36.6 -> 42.1 %). The two are now compared on their Niggli-reduced primitive edges (within 2 %, the battery's lattice-identity test), so a symmetric setting never loses to its own primitive, while 9qw8's C-centred re-index (reduced edges 5.7 % off pass 1's) is still judged against pass 1's lattice on the frames. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- rugnux/Rugnux.cpp | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 74e05a89d..bd60493e4 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -5516,13 +5516,25 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const auto &l = pass1_lattice.lattice; pass1_lattice.lattice = CrystalLattice(l.Vec0() * f, l.Vec1() * f, l.Vec2() * f); } - const double pass1_vol = std::abs(pass1_lattice.lattice - .ToPrimitive(pass1_lattice.search_result.centering).CalcVolume()); - // Same lattice as in the beam-centre check: same class and primitive volumes within 2 %. - const bool same_lattice = best.result.has_value() - && best.result->search_result.centering == pass1_lattice.search_result.centering - && best.result->search_result.system == pass1_lattice.search_result.system - && std::abs(best.vol - pass1_vol) <= 0.02 * pass1_vol; + const CrystalLattice pass1_primitive = pass1_lattice.lattice + .ToPrimitive(pass1_lattice.search_result.centering).NiggliReduce(); + const double pass1_vol = std::abs(pass1_primitive.CalcVolume()); + // The same lattice whatever setting or class it is held in - a centred cell and the + // triclinic primitive of the same points are one lattice, and the re-index's symmetric + // setting must not lose to its own unconstrained primitive on a frame count: the Niggli- + // reduced primitive edges agree within 2 % (the battery's lattice-identity test). + const auto reduced_edges = [](const CrystalLattice &l) { + const auto uc = l.GetUnitCell(); + std::array e = {uc.a, uc.b, uc.c}; + std::sort(e.begin(), e.end()); + return e; + }; + const auto e1 = reduced_edges(pass1_primitive); + const auto e2 = reduced_edges(best.result->lattice + .ToPrimitive(best.result->search_result.centering).NiggliReduce()); + bool same_lattice = true; + for (int i = 0; i < 3; i++) + same_lattice = same_lattice && std::abs(e2[i] - e1[i]) <= 0.02 * e1[i]; const int floor = static_cast(validation.size()) / 6; const int pass1_score = (best.score < floor || !same_lattice) ? count_indexed(*indexer, pass1_lattice) : -1; From f69339ce6387234c8e3fdde14f9918238dd50725 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 3 Oct 2026 23:08:12 +0200 Subject: [PATCH 06/13] RotationScaleMerge: fulls-only per-frame scaling (prototype), pooled over sparse frames Prototype switch --no-scale-partials: skip the per-frame scaling of the partials and let the fulls (scale-fulls) carry the per-frame scale, as XDS does. Without partial scaling a frame whose fulls number fewer than MIN_REFLECTIONS is fitted over the nearest frames on either side that together hold enough (PoolHalfWidth), on the host and in the GPU kernel. Default path unchanged (p.mtz md5 identical on the three profiling sets). Why: on fine-sliced small-molecule sweeps the per-frame partial scale and the partiality model are degenerate within a rocking curve, and the fit swings G 0.23..1.2 with a 180 deg period (XDS's own frame scale: 0.79..0.99). That imprints an hkl-dependent bias common to all equivalents, which R_meas/CC1/2/ISa cannot see but a refinement against the known structure does. And scale-fulls never fitted a frame on such data: a full is filed under one frame, about 8 per frame, below MIN_REFLECTIONS, so every frame kept G = 1. Measured with SHELXL refining the COD structures (R1 >4sig), default -> switch: aspirin 20 keV 0.096 -> 0.062, aspirin 25 keV 0.094 -> 0.046, citric acid 0.161 -> 0.109, HEPES 0.090 -> 0.070 (XDS 0.030-0.038). Proteins lose ISa with the switch (myob 9.1 -> 7.6, cytc 25.8 -> 13.7), so it is not a default; the choice is to be made from the data. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- common/ScalingSettings.cpp | 9 ++++ common/ScalingSettings.h | 3 ++ .../scale_merge/RotationScaleMerge.cpp | 41 +++++++++++++--- .../scale_merge/RotationScaleMerge.h | 7 +++ .../scale_merge/RotationScaleMergeGPU.cu | 48 ++++++++++++------- .../scale_merge/RotationScaleMergeGPU.h | 4 +- rugnux/rugnux_cli.cpp | 9 ++++ tests/RotationScaleWalkTest.cpp | 13 +++++ 8 files changed, 110 insertions(+), 24 deletions(-) diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index 2016afba5..4fabbd44d 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -208,6 +208,15 @@ double ScalingSettings::GetSmoothGDegrees() const { return smooth_g_deg; } +ScalingSettings &ScalingSettings::ScalePartials(bool input) { + scale_partials = input; + return *this; +} + +bool ScalingSettings::GetScalePartials() const { + return scale_partials; +} + ScalingSettings &ScalingSettings::RelativeBDegrees(double input) { if (input < 0) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Relative-B batch width must be non-negative"); diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index 198c6366c..e9b2cd72e 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -98,6 +98,7 @@ class ScalingSettings { // degrees (like XDS DELPHI), converted to an odd frame window from the oscillation step; this keeps // the smoothing physical (independent of frame slicing). 0 = off. A no-op without rot3d. double smooth_g_deg = 0.0; + bool scale_partials = true; // Per-batch relative-B on the rot3d fulls (beyond the single global decay slope): bin frames into // rotation-range batches of this width in degrees and refine one relative Debye-Waller B per batch, so @@ -143,6 +144,7 @@ public: ScalingSettings& IceMinScore(float input); ScalingSettings& IceMinSpotRatio(float input); ScalingSettings& SmoothGDegrees(double input); + ScalingSettings& ScalePartials(bool input); ScalingSettings& RelativeBDegrees(double input); ScalingSettings& RfreeFraction(double input); @@ -182,6 +184,7 @@ public: [[nodiscard]] float GetIceMinScore() const; [[nodiscard]] float GetIceMinSpotRatio() const; [[nodiscard]] double GetSmoothGDegrees() const; + [[nodiscard]] bool GetScalePartials() const; [[nodiscard]] double GetRelativeBDegrees() const; [[nodiscard]] double GetRfreeFraction() const; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 81dabcd22..2b71049b4 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1407,17 +1407,34 @@ void RotationScaleMerge::ReduceScalingGroupMeans(int n_groups, const std::vector }); } +int RotationScaleMerge::PoolHalfWidth(const std::vector &fcount, int f) { + const int n = static_cast(fcount.size()); + int64_t total = fcount[f]; + int h = 0; + while (total < static_cast(MIN_REFLECTIONS) && h < n) { + ++h; + if (f - h >= 0) total += fcount[f - h]; + if (f + h < n) total += fcount[f + h]; + } + return h; +} + template void RotationScaleMerge::FitPerFrameG(const std::vector &obs, const std::vector &fstart, const std::vector &fcount, const std::vector &group_mean_in, bool unity, std::vector &g) { std::vector scaled(fstart.size(), 0); - ParallelFor(static_cast(fstart.size()), nthreads, [&](int f) { + const int n_fr = static_cast(fstart.size()); + ParallelFor(n_fr, nthreads, [&](int f) { + // A full is a whole rocking event filed under one frame, so on a sparse sweep a frame holds a + // handful of them - too few to fit a scale on, and the frame was left unscaled, which on such a + // sweep was every frame. Those frames are fitted over the nearest frames on either side that + // together hold enough; a frame with enough of its own is fitted on its own, as before. + const int h = unity && pool_sparse_fulls ? PoolHalfWidth(fcount, f) : 0; std::vector so; - so.reserve(fcount[f]); - const int lo = fstart[f], hi = fstart[f] + fcount[f]; - for (int i = lo; i < hi; ++i) { + for (int j = std::max(0, f - h); j <= std::min(n_fr - 1, f + h); ++j) + for (int i = fstart[j]; i < fstart[j] + fcount[j]; ++i) { const auto &o = obs[i]; if (o.group < 0) continue; if (o.on_ice) continue; @@ -5607,9 +5624,15 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st "curves span {} frames)", smooth_window, smooth_window * osc_deg, smooth_g_deg, rocking_event_frames_at_start); } + // Without partial scaling every frame keeps G = 1 here and its scale comes from the fulls alone + // (scale-fulls below), as in XDS: a partial's scale and an error of the partiality model are the + // same thing within a rocking curve, and on a sparse fine-sliced sweep the fit takes the one for + // the other. + const bool scale_partials = s.GetScalePartials(); + pool_sparse_fulls = !scale_partials; ScalingLoopOutcome partial_loop; #ifdef JFJOCH_USE_CUDA - if (gpu_active_) { + if (gpu_active_ && scale_partials) { // The scaling loop runs on the GPU one iteration per call, and corr stays RESIDENT across // scaling -> smooth-G -> CC -> combine (and across passes, exactly as the old host round-trip // did). Only the per-frame G/scaled come back after each iteration, for the gauge pin and the @@ -5630,7 +5653,11 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st scaled_on_gpu = true; } #endif - if (!scaled_on_gpu) { + if (!scale_partials) { + std::fill(g_partial.begin(), g_partial.end(), 1.0); + frame_scaled_scratch.assign(n_frames, 1); + partial_loop.converged = true; + } else if (!scaled_on_gpu) { const PartialLoopKey key{x.GetSpaceGroupOrP1().xhm(), merge_friedel, d_min_limit, d_max_limit, min_partiality, smooth_window, scaling_iter}; const auto memo = std::find_if(partial_loop_memos.begin(), partial_loop_memos.end(), @@ -5920,7 +5947,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st std::vector scaled_dev(n_frames); fulls_loop = RunScalingLoop("fulls", g_full, f_count, smooth_window, /*release_restraint=*/true, [&] { - gpu_->ScaleFulls(1, min_partiality); + gpu_->ScaleFulls(1, min_partiality, pool_sparse_fulls); gpu_->GetG(g_dev.data(), scaled_dev.data()); for (int f = 0; f < n_frames; ++f) if (scaled_dev[f]) g_full[f] = g_dev[f]; diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 79278560f..30da42438 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -43,6 +43,10 @@ // Stills use the per-image ScaleOnTheFly (fixed partiality) instead. class RotationScaleMerge { public: + // How many frames on either side a frame's fulls scale is fitted over (used without partial + // scaling): none when the frame holds MIN_REFLECTIONS fulls itself, else the fewest that bring the + // pooled count there. + static int PoolHalfWidth(const std::vector &fcount, int f); struct Result { std::vector merged; MergeStatistics statistics; @@ -213,6 +217,9 @@ private: // this is how far they may go before giving up and saying so. int scaling_iter = 100; bool scale_fulls = true; + // Without partial scaling the fulls carry the whole per-frame scale, and a sparse frame is fitted + // over its neighbours (PoolHalfWidth). + bool pool_sparse_fulls = false; bool refine_decay_b = false; // per-time-block Debye-Waller decay correction (radiation damage) int absorption_iter = 0; // >0: fit a goniometer-frame absorption surface over this many iterations int modulation_iter = 0; // >0: fit a detector-plane modulation (flat-field) surface, this many iterations diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 66aeabacb..06b85d5bd 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -123,16 +123,31 @@ namespace { const int32_t *__restrict__ frame_count, const float *__restrict__ I, const double *__restrict__ inv_sigma, const float *__restrict__ sco_coeff, const uint8_t *__restrict__ sco_ok, - const int32_t *__restrict__ perm, + const int32_t *__restrict__ perm, bool pool, double *__restrict__ g, uint8_t *__restrict__ scaled) { const int f = blockIdx.x; if (f >= n_frames) return; - const int lo = frame_start[f], hi = frame_start[f] + frame_count[f]; + // pool (the fulls): a frame holding fewer than MIN_REFLECTIONS is fitted over the fewest frames + // on either side that bring the count there - RotationScaleMerge::PoolHalfWidth, same rule. + __shared__ int s_h; + if (threadIdx.x == 0) { + long total = frame_count[f]; + int h = 0; + while (pool && total < long(MIN_REFLECTIONS) && h < n_frames) { + ++h; + if (f - h >= 0) total += frame_count[f - h]; + if (f + h < n_frames) total += frame_count[f + h]; + } + s_h = h; + } + __syncthreads(); + const int j_lo = max(0, f - s_h), j_hi = min(n_frames - 1, f + s_h); __shared__ double sh[BLK]; long cnt_local = 0; - for (int i = lo + threadIdx.x; i < hi; i += blockDim.x) - if (sco_ok[perm ? perm[i] : i]) ++cnt_local; + for (int j = j_lo; j <= j_hi; ++j) + for (int i = frame_start[j] + threadIdx.x; i < frame_start[j] + frame_count[j]; i += blockDim.x) + if (sco_ok[perm ? perm[i] : i]) ++cnt_local; const double cnt = BlockReduceSum(double(cnt_local), sh); __shared__ double s_cnt; if (threadIdx.x == 0) s_cnt = cnt; @@ -140,15 +155,16 @@ namespace { if (s_cnt < MIN_REFLECTIONS) return; // leave g[f]/scaled[f] as-is double num = 0.0, den = 0.0; - for (int i = lo + threadIdx.x; i < hi; i += blockDim.x) { - const int a = perm ? perm[i] : i; - if (!sco_ok[a]) continue; - const double coeff = sco_coeff[a]; - const double w = inv_sigma[a]; - const double w2 = w * w; - num += w2 * coeff * double(I[a]); - den += w2 * coeff * coeff; - } + for (int j = j_lo; j <= j_hi; ++j) + for (int i = frame_start[j] + threadIdx.x; i < frame_start[j] + frame_count[j]; i += blockDim.x) { + const int a = perm ? perm[i] : i; + if (!sco_ok[a]) continue; + const double coeff = sco_coeff[a]; + const double w = inv_sigma[a]; + const double w2 = w * w; + num += w2 * coeff * double(I[a]); + den += w2 * coeff * coeff; + } const double tnum = BlockReduceSum(num, sh); __syncthreads(); const double tden = BlockReduceSum(den, sh); if (threadIdx.x == 0) { @@ -1081,7 +1097,7 @@ void RotationScaleMergeGPU::ScalePartials(int iters, double min_partiality, bool d.sco_coeff.get(), d.sco_ok.get()); CudaCheck(cudaGetLastError(), "PrepScaleObsKernel launch"); FitPerFrameGKernel<<s()>>>(d.n_frames, d.frame_start.get(), d.frame_count.get(), - d.I.get(), d.inv_sigma.get(), d.sco_coeff.get(), d.sco_ok.get(), nullptr, d.g.get(), d.scaled.get()); + d.I.get(), d.inv_sigma.get(), d.sco_coeff.get(), d.sco_ok.get(), nullptr, /*pool=*/false, d.g.get(), d.scaled.get()); CudaCheck(cudaGetLastError(), "FitPerFrameGKernel launch"); UpdateCorrKernel<<s()>>>(d.n_obs, d.frame.get(), d.prescaling_corr.get(), d.partiality.get(), d.g.get(), d.scaled.get(), d.corr.get()); @@ -1497,7 +1513,7 @@ void RotationScaleMergeGPU::ResetFullsScale() { CudaCheck(cudaStreamSynchronize(impl_->s()), "reset fulls scale sync"); } -void RotationScaleMergeGPU::ScaleFulls(int iters, double min_partiality) { +void RotationScaleMergeGPU::ScaleFulls(int iters, double min_partiality, bool pool) { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; const int nf = d.n_fulls; @@ -1523,7 +1539,7 @@ void RotationScaleMergeGPU::ScaleFulls(int iters, double min_partiality) { CudaCheck(cudaGetLastError(), "PrepScaleObsKernel launch"); FitPerFrameGKernel<<s()>>>(d.n_frames, d.f_frame_start.get(), d.f_frame_count.get(), d.f_I.get(), d.f_inv_sigma.get(), - d.f_sco_coeff.get(), d.f_sco_ok.get(), d.f_frame_perm.get(), d.g.get(), d.scaled.get()); + d.f_sco_coeff.get(), d.f_sco_ok.get(), d.f_frame_perm.get(), pool, d.g.get(), d.scaled.get()); CudaCheck(cudaGetLastError(), "FitPerFrameGKernel launch"); UpdateCorrKernel<<s()>>>(nf, d.f_frame.get(), d.f_rlp.get(), d.f_partiality.get(), d.g.get(), d.scaled.get(), d.f_corr.get()); diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index 766e2f640..98aa1b18a 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -184,7 +184,9 @@ public: // Run `iters` of the Unity scaling loop on the resident fulls (reduce group means -> per-frame LS G // -> update corr), in place on the fulls' working corr. Requires SetFullsFrameCSR + SetFullsGroups // and ResetFullsScale. - void ScaleFulls(int iters, double min_partiality); + // pool: a frame with fewer than MIN_REFLECTIONS fulls is fitted over its neighbours + // (RotationScaleMerge::PoolHalfWidth). + void ScaleFulls(int iters, double min_partiality, bool pool); // The fulls' counterpart of SmoothCorr: f_corr[i] *= ratio[f_frame[i]] where apply[f]. void SmoothFullsCorr(const uint8_t *apply, const double *ratio); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 4eb457f6b..9cb7403f5 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -157,6 +157,7 @@ void print_usage() { std::cout << " --write-process-h5 Also write the (large) _process.h5 when merging (default: only .mtz/.cif when merging)" << std::endl; std::cout << " --finalist-ledger Report the full-resolution evidence for each space group the search considered, not only the one it adopted (report-only; the decision is unchanged)" << std::endl; std::cout << " --developer Write the full _report.txt: the pipeline-internal keys (the anisotropy gate, the space-group candidate and operator tables, the model-fit null, the sweep internals) and the long explanations, which the default report leaves out. Nothing is computed differently - the same report, rendered in full" << std::endl; + std::cout << " --no-scale-partials rot3d: do not fit per-frame scales on the partials; the per-frame scale comes from the fulls (scale-fulls) alone, as in XDS" << std::endl; std::cout << " --smooth-g[=deg] rot3d: smooth per-frame scale G over a deg-degree rotation range (XDS DELPHI-like) before the combine (default: 5 for rot3d; 0 = off)" << std::endl; std::cout << " --relative-b[=deg] rot3d: fit a per-batch relative-B (beyond the single decay slope) over deg-degree batches; cross-validated (default: 10 deg when bare; off otherwise)" << std::endl; std::cout << " --no-scaling-corrections rot3d: disable the (default-on) decay + absorption + modulation correction surfaces fitted on the fulls after scale-fulls" << std::endl; @@ -282,6 +283,7 @@ enum { OPT_MOSAICITY, OPT_PREDICTION_MOSAICITY, OPT_SMOOTH_G, + OPT_NO_SCALE_PARTIALS, OPT_RELATIVE_B, OPT_NO_SCALING_CORRECTIONS, OPT_NO_EXPECTED_VARIANCE_MERGE, @@ -351,6 +353,7 @@ static option long_options[] = { {"finalist-ledger", no_argument, nullptr, OPT_FINALIST_LEDGER}, {"developer", no_argument, nullptr, OPT_DEVELOPER}, {"smooth-g", optional_argument, nullptr, OPT_SMOOTH_G}, + {"no-scale-partials", no_argument, nullptr, OPT_NO_SCALE_PARTIALS}, {"relative-b", optional_argument, nullptr, OPT_RELATIVE_B}, {"no-scaling-corrections", no_argument, nullptr, OPT_NO_SCALING_CORRECTIONS}, {"no-expected-variance-merge", no_argument, nullptr, OPT_NO_EXPECTED_VARIANCE_MERGE}, @@ -773,6 +776,7 @@ static int RunRugnux(int argc, char **argv) { std::optional beam_x, beam_y, detector_distance_mm, wavelength_A, rot1_rad, rot2_rad, rot3_rad, polarization_factor; bool detector_mirror_y = false; int64_t detector_quarter_turns = 0; + bool no_scale_partials = false; // --no-scale-partials (prototype) std::optional smooth_g_deg_arg; // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise std::optional relative_b_deg_arg; // --relative-b[=deg]; per-batch relative-B width, 0 (off) unless given bool no_scaling_corrections = false; // --no-scaling-corrections: disable rot3d decay+absorption+modulation surfaces @@ -1207,6 +1211,9 @@ static int RunRugnux(int argc, char **argv) { case OPT_DEVELOPER: provenance.developer = true; break; + case OPT_NO_SCALE_PARTIALS: + no_scale_partials = true; + break; case OPT_SMOOTH_G: smooth_g_deg_arg = optarg ? parse_double_arg(optarg, "--smooth-g", logger) : SMOOTH_G_DEFAULT_DEG; break; @@ -1652,6 +1659,7 @@ static int RunRugnux(int argc, char **argv) { outlier_reject_nsigma.value_or(scaling_settings.GetOutlierRejectNsigma())); scaling_settings.ScaleFulls(scale_fulls_arg.value_or(scaling_settings.GetScaleFulls())); scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(scaling_settings.GetSmoothGDegrees())); + if (no_scale_partials) scaling_settings.ScalePartials(false); scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0)); // opt-in only; default off if (no_scaling_corrections) scaling_settings.CorrectionSurfaces(false); @@ -2538,6 +2546,7 @@ static int RunRugnux(int argc, char **argv) { ScalingSettings scaling_settings = RugnuxDefaultScalingSettings(rotation_indexing); scaling_settings.ScaleFulls(scale_fulls); scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(scaling_settings.GetSmoothGDegrees())); + if (no_scale_partials) scaling_settings.ScalePartials(false); scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0)); // opt-in only; default off if (no_scaling_corrections) scaling_settings.CorrectionSurfaces(false); diff --git a/tests/RotationScaleWalkTest.cpp b/tests/RotationScaleWalkTest.cpp index 562d46a24..f01334b3f 100644 --- a/tests/RotationScaleWalkTest.cpp +++ b/tests/RotationScaleWalkTest.cpp @@ -5,6 +5,7 @@ #include #include "../rugnux/Rugnux.h" +#include "../image_analysis/scale_merge/RotationScaleMerge.h" namespace { // A synthetic sweep whose stage turned `true_scale` times the stored angles. Scored at a scale k, @@ -92,3 +93,15 @@ TEST_CASE("WalkRotationScale_StoredAnglesStand", "[RotationScale]") { CHECK(sweep.index_calls == 0); } } + +TEST_CASE("PoolHalfWidth_FitsSparseFramesOverTheirNeighbours", "[RotationScale]") { + // 20 fulls of its own: fitted on its own. + CHECK(RotationScaleMerge::PoolHalfWidth({0, 20, 0}, 1) == 0); + // 5 per frame: two frames either side bring 25 >= 20. + const std::vector five(11, 5); + CHECK(RotationScaleMerge::PoolHalfWidth(five, 5) == 2); + // At the end of the sweep the window grows on the one side there is. + CHECK(RotationScaleMerge::PoolHalfWidth(five, 0) == 3); + // A sweep that never holds enough stops at its length. + CHECK(RotationScaleMerge::PoolHalfWidth({1, 1, 1}, 1) == 3); +} From 20ab92fcaba21293209c67daefb5fc97832e235f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 4 Oct 2026 00:05:36 +0200 Subject: [PATCH 07/13] Build: no FMA contraction; first-pass cell refinement independent of SIMD width The same source built with and without -march=x86-64-v3 gave different results on battery sets (9zmu axis-harmonic supercell arbiter fired in one build only; 8rud resolution cut 1.69 vs 1.70 A; myob_x06da_powder_1 cut 1.07 vs 1.42 A; 7mzt short-axis first pass; lcystine_x10sa_20keV indexed vs no lattice), against the rule that no decision may depend on compiler flags. Two mechanisms, found by building the merged tree four ways (baseline, x86-64-v2, x86-64-v3, x86-64-v3 -mno-fma) with and without -ffp-contract=off and comparing p.mtz: - FMA contraction. GCC contracts a*b+c whenever the target has FMA. With -ffp-contract=off the x86-64-v3 build gives a p.mtz byte-identical to the baseline build on 8 of 9 sets (7mzt, 8rud, 9zmu, insu_I_x06da_5keV_2, myob_x06da_powder_1, myob_x10sa, cytc_x10sa, thau_x10sa_16keV), on both the GPU and the CPU build. Cost: none measurable (user core-s, CPU build, x86-64-v3 vs the same with -ffp-contract=off: 1877/1874, 3285/3253, 2212/2195 on myob/cytc/thau; GPU likewise within noise). Set project-wide for C, C++ and CUDA host code; MSVC does not contract under /fp:precise. - SIMD width. lcystine still differed: baseline and x86-64-v2 (128-bit) agreed, x86-64-v3 with or without FMA (256-bit) disagreed - Eigen's HouseholderQR in the FFT indexer's candidate refinement (PostIndexingRefinement.cpp) reduces column norms over all spots in packets of the target width. Replaced by the 3x3 normal equations summed in spot order in double. All four builds now agree on all nine sets. Changes results of the default x86-64-v3 build (contraction off); lcystine_x10sa_20keV now gives no lattice in every build (its first pass is a knife-edge: 0/60 vs 9/60 validation frames before). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- CMakeLists.txt | 12 +++++++++++ .../indexing/PostIndexingRefinement.cpp | 20 +++++++++++++++---- 2 files changed, 28 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1caea9809..dfc4da29e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -41,6 +41,18 @@ IF (MSVC) ADD_COMPILE_OPTIONS($<$:-Xcompiler=/Zc:preprocessor>) ENDIF() +# No contraction of a*b+c into a fused multiply-add. GCC contracts whenever the target has FMA, so a +# -march=x86-64-v3 build (CI, production) and a baseline build rounded differently, and knife-edge +# decisions - which candidate cell a first pass indexes with, a resolution cut, a symmetry gate - +# came out differently from the same source and data. Off, every x86-64 level computes the same bits +# (measured: identical p.mtz on nine battery sets, no measurable time cost), and aarch64, where FMA is +# always there and GCC contracts by default, follows the same arithmetic. MSVC does not contract under +# its default /fp:precise. Host code of .cu sources gets it through nvcc; device code is not affected. +IF (NOT MSVC) + ADD_COMPILE_OPTIONS($<$:-ffp-contract=off>) + ADD_COMPILE_OPTIONS($<$:-Xcompiler=-ffp-contract=off>) +ENDIF() + SET(JFJOCH_WRITER_ONLY OFF CACHE BOOL "Compile HDF5 writer only") SET(JFJOCH_INSTALL_DRIVER_SOURCE OFF CACHE BOOL "Install kernel driver source (ignored if building writer only; necessary for RPM building)") SET(JFJOCH_USE_CUDA ON CACHE BOOL "Compile Jungfraujoch with CUDA") diff --git a/image_analysis/indexing/PostIndexingRefinement.cpp b/image_analysis/indexing/PostIndexingRefinement.cpp index d8d210eb0..a31bb7de3 100644 --- a/image_analysis/indexing/PostIndexingRefinement.cpp +++ b/image_analysis/indexing/PostIndexingRefinement.cpp @@ -82,7 +82,6 @@ namespace { const unsigned nspots = spots.rows(); const unsigned ncells = scores.rows(); VectorX below{nspots}; - MatrixX3 sel{nspots, 3u}; Mx3 resid{nspots, 3u}; Mx3 miller{nspots, 3u}; M3 cell; @@ -111,9 +110,22 @@ namespace { break; threshold *= cifssr.threshold_contraction; - sel.colwise() = below; - HouseholderQR qr{sel.select(spots, .0f)}; - cell = qr.solve(sel.select(miller, .0f)); + // Least squares spots * cell = miller over the spots below the threshold, by the normal + // equations summed in spot order in double. A QR over all the spots reduces its column + // norms in SIMD packets of the target's width, so its cell - and which candidate a + // borderline first pass indexed with - depended on -march; these sums do not. + Matrix3d ata = Matrix3d::Zero(); + Matrix3d atb = Matrix3d::Zero(); + for (unsigned i = 0; i < nspots; i++) { + if (!below[i]) + continue; + for (int r = 0; r < 3; r++) + for (int c = 0; c < 3; c++) { + ata(r, c) += double(spots(i, r)) * double(spots(i, c)); + atb(r, c) += double(spots(i, r)) * double(miller(i, c)); + } + } + cell = (ata.inverse() * atb).cast(); } resid = CalculateResiduals(spots, cell); From ba92417210e06899810e117e3b9778e912fd3fda Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 4 Oct 2026 00:34:39 +0200 Subject: [PATCH 08/13] RotationScaleMerge: per-frame scale from the fulls alone on sparse sweeps (data-decided) Replaces the prototype --no-scale-partials switch with a rule read off the data: a merge whose sweep holds fewer than 50 rocking events per frame takes its per-frame scale from the fulls alone (scale-fulls, a sparse frame fitted over its neighbours - PoolHalfWidth), and logs that it did. The rocking-event walk already run for the smoothing window now also returns its event count. Within one rocking curve a partial's scale and an error of the partiality model are the same thing; on a sparse fine-sliced sweep the partial fit takes the one for the other and imprints an hkl-dependent bias common to all equivalents. Measured populations: small-molecule sweeps 2.6-23 events per frame, protein sets 84-900; on the proteins the fulls-only scale leaves model R-free unchanged (+-0.002 over the smoke tier's open-arm sets) but lowers ISa, so they keep the partial scale. p.mtz md5 unchanged on the three profiling sets, GPU and CPU builds. SHELXL against the COD structures, R1(>4sig) rc174 -> this: aspirin 20 keV 0.096 -> 0.062, aspirin 25 keV 0.094 -> 0.046, citric acid 0.161 -> 0.109, HEPES 0.090 -> 0.070 (XDS 0.030-0.038); SHELXL's weight a comes off its 0.2 cap on all four. CPU and GPU paths agree. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- common/ScalingSettings.cpp | 9 ---- common/ScalingSettings.h | 3 -- .../scale_merge/RotationScaleMerge.cpp | 42 +++++++++++++------ .../scale_merge/RotationScaleMerge.h | 7 +++- rugnux/rugnux_cli.cpp | 9 ---- 5 files changed, 35 insertions(+), 35 deletions(-) diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index 4fabbd44d..2016afba5 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -208,15 +208,6 @@ double ScalingSettings::GetSmoothGDegrees() const { return smooth_g_deg; } -ScalingSettings &ScalingSettings::ScalePartials(bool input) { - scale_partials = input; - return *this; -} - -bool ScalingSettings::GetScalePartials() const { - return scale_partials; -} - ScalingSettings &ScalingSettings::RelativeBDegrees(double input) { if (input < 0) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Relative-B batch width must be non-negative"); diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index e9b2cd72e..198c6366c 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -98,7 +98,6 @@ class ScalingSettings { // degrees (like XDS DELPHI), converted to an odd frame window from the oscillation step; this keeps // the smoothing physical (independent of frame slicing). 0 = off. A no-op without rot3d. double smooth_g_deg = 0.0; - bool scale_partials = true; // Per-batch relative-B on the rot3d fulls (beyond the single global decay slope): bin frames into // rotation-range batches of this width in degrees and refine one relative Debye-Waller B per batch, so @@ -144,7 +143,6 @@ public: ScalingSettings& IceMinScore(float input); ScalingSettings& IceMinSpotRatio(float input); ScalingSettings& SmoothGDegrees(double input); - ScalingSettings& ScalePartials(bool input); ScalingSettings& RelativeBDegrees(double input); ScalingSettings& RfreeFraction(double input); @@ -184,7 +182,6 @@ public: [[nodiscard]] float GetIceMinScore() const; [[nodiscard]] float GetIceMinSpotRatio() const; [[nodiscard]] double GetSmoothGDegrees() const; - [[nodiscard]] bool GetScalePartials() const; [[nodiscard]] double GetRelativeBDegrees() const; [[nodiscard]] double GetRfreeFraction() const; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 2b71049b4..fef918578 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -36,6 +36,16 @@ namespace { // These mirror the per-image ScaleOnTheFly / Merge rocking-curve physics verbatim so this flat // implementation is numerically identical - see the comments there for the details. constexpr size_t MIN_REFLECTIONS = 20; // per-frame scale needs at least this many + // Per-frame scaling of the PARTIALS needs a frame to carry many rocking events: within one rocking + // curve a change of scale and an error of the partiality model are the same thing, and only events + // caught at different points of their curves tell them apart. On a fine-sliced sparse sweep the fit + // takes one for the other - G swung 0.23..1.2 with a 180 deg period on a small-molecule sweep whose + // own frame scale varies 0.79..0.99 (XDS), and the merge carried an hkl-dependent bias that the + // equivalents cannot see but a refinement against the known structure does (R1 0.096 -> 0.062 + // without it). Below this the per-frame scale comes from the fulls alone, as XDS takes it. Measured + // on the battery: the small-molecule sweeps hold 2.6 to 23 events per frame, the protein sets 84 to + // 900 - and on the proteins the fulls-only scale leaves model R-free unchanged (+-0.002). + constexpr double MIN_EVENTS_PER_FRAME_SCALE_PARTIALS = 50.0; constexpr int64_t MIN_REFLECTIONS_FOR_IMAGE_CC = 20; // below this a frame's CC means nothing // The per-frame scaling loop has settled when the scales moved less than this between two // iterations, as the rms of |log(G_new/G_old)| over the frames fitted in both. @@ -591,7 +601,7 @@ void RotationScaleMerge::Ingest() { // (see the header). Take its answer now and hand them back. The dump path keeps the full Obs // array: the CPU combine is what writes the dump. if (resident_ingest) { - rocking_event_frames_at_ingest = RockingEventFrames(); + rocking_event_frames_at_ingest = RockingEventFrames(&rocking_events_at_ingest); partials_released = true; } #endif @@ -2473,24 +2483,26 @@ namespace { // image order, bridged by the same max_frame_gap - taken as the median over the run's events. It is // the finest stretch of the sweep this file is allowed to speak about: the partials of one event are // welded into a single full, so two frames closer together than this are not separable observations. -int RotationScaleMerge::RockingEventFrames() const { - if (partials_released) +int RotationScaleMerge::RockingEventFrames(int64_t *n_events) const { + if (partials_released) { + if (n_events) *n_events = rocking_events_at_ingest; return rocking_event_frames_at_ingest; + } if (resident_ingest) // Ingest's own call, before the ingest arrays are handed back return RockingEventFramesOver( [&](int i) { return ingest_rock_ok[i] != 0; }, - [&](int i) { return ingest_image_number[i]; }); + [&](int i) { return ingest_image_number[i]; }, n_events); return RockingEventFramesOver( [&](int i) { const Obs &o = partials[i]; return std::isfinite(o.corr) && o.corr > 0.0f && std::isfinite(o.I) && std::isfinite(o.sigma) && o.sigma > 0.0f; }, - [&](int i) { return partials[i].image_number; }); + [&](int i) { return partials[i].image_number; }, n_events); } template -int RotationScaleMerge::RockingEventFramesOver(UsableFn usable, ImgFn img) const { +int RotationScaleMerge::RockingEventFramesOver(UsableFn usable, ImgFn img, int64_t *n_events) const { // Blocks of runs on all threads, each counting into a histogram of its own. The counts are // integers, so the blocks add up to the serial histogram whatever the split. constexpr int RUNS_PER_BLOCK = 16384; @@ -2525,6 +2537,7 @@ int RotationScaleMerge::RockingEventFramesOver(UsableFn usable, ImgFn img) const hist[w] += h[w]; n_event += h[w]; } + if (n_events) *n_events = n_event; int acc = 0; for (int w = 1; w <= n_frames + 1; ++w) if ((acc += hist[w]) * 2 >= n_event) @@ -5616,7 +5629,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st int smooth_window = 1; if (smooth_g_deg > 0.0 && osc_deg > 1e-6) { if (rocking_event_frames_at_start < 0) - rocking_event_frames_at_start = RockingEventFrames(); + rocking_event_frames_at_start = RockingEventFrames(&rocking_events_at_start); smooth_window = std::max(static_cast(std::lround(smooth_g_deg / osc_deg)), SMOOTH_G_MIN_ROCKING_EVENTS * rocking_event_frames_at_start); if (smooth_window % 2 == 0) ++smooth_window; @@ -5624,12 +5637,17 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st "curves span {} frames)", smooth_window, smooth_window * osc_deg, smooth_g_deg, rocking_event_frames_at_start); } - // Without partial scaling every frame keeps G = 1 here and its scale comes from the fulls alone - // (scale-fulls below), as in XDS: a partial's scale and an error of the partiality model are the - // same thing within a rocking curve, and on a sparse fine-sliced sweep the fit takes the one for - // the other. - const bool scale_partials = s.GetScalePartials(); + // Without partial scaling (MIN_EVENTS_PER_FRAME_SCALE_PARTIALS) every frame keeps G = 1 here and + // its scale comes from the fulls alone (scale-fulls below). + if (rocking_event_frames_at_start < 0) + rocking_event_frames_at_start = RockingEventFrames(&rocking_events_at_start); + const double events_per_frame = n_frames > 0 ? static_cast(rocking_events_at_start) / n_frames : 0.0; + const bool scale_partials = events_per_frame >= MIN_EVENTS_PER_FRAME_SCALE_PARTIALS; pool_sparse_fulls = !scale_partials; + if (!scale_partials) + logger.Info("Per-frame scale from the fulls alone: {:.1f} rocking events per frame, under the {:.0f} " + "a partial's scale needs to stay apart from the partiality model", events_per_frame, + MIN_EVENTS_PER_FRAME_SCALE_PARTIALS); ScalingLoopOutcome partial_loop; #ifdef JFJOCH_USE_CUDA if (gpu_active_ && scale_partials) { diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 30da42438..c5261da57 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -276,6 +276,7 @@ private: // full-size copy of the partials, gigabytes on a fine-sliced long axis - straight back. bool partials_released = false; int rocking_event_frames_at_ingest = 0; + int64_t rocking_events_at_ingest = 0; std::vector frame_start, frame_count; // CSR ranges of `partials` per frame std::vector frame_cell_ok; // per-frame cell-consistency mask (1 = kept) std::vector finite_ok; // per-obs AcceptReflection finiteness (immutable; 1 = kept) @@ -306,6 +307,7 @@ private: // RockingEventFrames at the start of a Run, where corr has just been restored to corr_ingested - so // the same for every Run on the same ingest; -1 until the first Run takes it. int rocking_event_frames_at_start = -1; + int64_t rocking_events_at_start = 0; // the events that walk counted // Raw-hkl ordering, built ONCE by Ingest and reused: `perm` lists partial indices sorted by // (raw h,k,l, image_number); each distinct raw hkl is a contiguous run [rawrun_start, +count) of it. @@ -626,9 +628,10 @@ private: // The frames one rocking event spans, as the combine cuts them: the median over the run's events. // Walked on the partials, because the combine also runs on the GPU and a full keeps only the frame // of its peak partial. - [[nodiscard]] int RockingEventFrames() const; + // n_events (if given) receives how many rocking events the walk counted. + [[nodiscard]] int RockingEventFrames(int64_t *n_events = nullptr) const; template - [[nodiscard]] int RockingEventFramesOver(UsableFn usable, ImgFn img) const; + [[nodiscard]] int RockingEventFramesOver(UsableFn usable, ImgFn img, int64_t *n_events) const; // Per-batch delta-CC1/2 on the corrected fulls: measure what keeping each batch of the sweep costs // the merged intensities, convict the batches that cost significantly, slide the conviction's edges // onto the frames that carry it, and turn the result into the disposition ledger. Fills diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 9cb7403f5..4eb457f6b 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -157,7 +157,6 @@ void print_usage() { std::cout << " --write-process-h5 Also write the (large) _process.h5 when merging (default: only .mtz/.cif when merging)" << std::endl; std::cout << " --finalist-ledger Report the full-resolution evidence for each space group the search considered, not only the one it adopted (report-only; the decision is unchanged)" << std::endl; std::cout << " --developer Write the full _report.txt: the pipeline-internal keys (the anisotropy gate, the space-group candidate and operator tables, the model-fit null, the sweep internals) and the long explanations, which the default report leaves out. Nothing is computed differently - the same report, rendered in full" << std::endl; - std::cout << " --no-scale-partials rot3d: do not fit per-frame scales on the partials; the per-frame scale comes from the fulls (scale-fulls) alone, as in XDS" << std::endl; std::cout << " --smooth-g[=deg] rot3d: smooth per-frame scale G over a deg-degree rotation range (XDS DELPHI-like) before the combine (default: 5 for rot3d; 0 = off)" << std::endl; std::cout << " --relative-b[=deg] rot3d: fit a per-batch relative-B (beyond the single decay slope) over deg-degree batches; cross-validated (default: 10 deg when bare; off otherwise)" << std::endl; std::cout << " --no-scaling-corrections rot3d: disable the (default-on) decay + absorption + modulation correction surfaces fitted on the fulls after scale-fulls" << std::endl; @@ -283,7 +282,6 @@ enum { OPT_MOSAICITY, OPT_PREDICTION_MOSAICITY, OPT_SMOOTH_G, - OPT_NO_SCALE_PARTIALS, OPT_RELATIVE_B, OPT_NO_SCALING_CORRECTIONS, OPT_NO_EXPECTED_VARIANCE_MERGE, @@ -353,7 +351,6 @@ static option long_options[] = { {"finalist-ledger", no_argument, nullptr, OPT_FINALIST_LEDGER}, {"developer", no_argument, nullptr, OPT_DEVELOPER}, {"smooth-g", optional_argument, nullptr, OPT_SMOOTH_G}, - {"no-scale-partials", no_argument, nullptr, OPT_NO_SCALE_PARTIALS}, {"relative-b", optional_argument, nullptr, OPT_RELATIVE_B}, {"no-scaling-corrections", no_argument, nullptr, OPT_NO_SCALING_CORRECTIONS}, {"no-expected-variance-merge", no_argument, nullptr, OPT_NO_EXPECTED_VARIANCE_MERGE}, @@ -776,7 +773,6 @@ static int RunRugnux(int argc, char **argv) { std::optional beam_x, beam_y, detector_distance_mm, wavelength_A, rot1_rad, rot2_rad, rot3_rad, polarization_factor; bool detector_mirror_y = false; int64_t detector_quarter_turns = 0; - bool no_scale_partials = false; // --no-scale-partials (prototype) std::optional smooth_g_deg_arg; // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise std::optional relative_b_deg_arg; // --relative-b[=deg]; per-batch relative-B width, 0 (off) unless given bool no_scaling_corrections = false; // --no-scaling-corrections: disable rot3d decay+absorption+modulation surfaces @@ -1211,9 +1207,6 @@ static int RunRugnux(int argc, char **argv) { case OPT_DEVELOPER: provenance.developer = true; break; - case OPT_NO_SCALE_PARTIALS: - no_scale_partials = true; - break; case OPT_SMOOTH_G: smooth_g_deg_arg = optarg ? parse_double_arg(optarg, "--smooth-g", logger) : SMOOTH_G_DEFAULT_DEG; break; @@ -1659,7 +1652,6 @@ static int RunRugnux(int argc, char **argv) { outlier_reject_nsigma.value_or(scaling_settings.GetOutlierRejectNsigma())); scaling_settings.ScaleFulls(scale_fulls_arg.value_or(scaling_settings.GetScaleFulls())); scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(scaling_settings.GetSmoothGDegrees())); - if (no_scale_partials) scaling_settings.ScalePartials(false); scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0)); // opt-in only; default off if (no_scaling_corrections) scaling_settings.CorrectionSurfaces(false); @@ -2546,7 +2538,6 @@ static int RunRugnux(int argc, char **argv) { ScalingSettings scaling_settings = RugnuxDefaultScalingSettings(rotation_indexing); scaling_settings.ScaleFulls(scale_fulls); scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(scaling_settings.GetSmoothGDegrees())); - if (no_scale_partials) scaling_settings.ScalePartials(false); scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0)); // opt-in only; default off if (no_scaling_corrections) scaling_settings.CorrectionSurfaces(false); From 1baf926066d4bd9f617a0a84c20a046088b09ddf Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 4 Oct 2026 00:48:35 +0200 Subject: [PATCH 09/13] RotationScaleMerge: refit the error model about the merge's own mean The error model was fitted once, on deviations from a mean that weighs every full by its COUNTING variance, and the outlier test's median took the same weights. Where equivalents disagree beyond counting statistics, the low-count observations dominate that centre: on a strongly absorbing crystal (YAG, Ia-3d, equivalents spread over a factor 100 after scaling) the mean of (0 4 0) sat at 21k among observations from 11k to 1.9M, a ran into its bound (100), b came out at 800% internally, and the six-sigma test about the biased median removed 63% of the observations - the strong ones. Now, after the first fit, the model is refitted about the model-weighted mean (each full weighted by the variance the fitted model gives it at the reflection's mean) until a and b settle, and the rejection median takes the same weights. Where counting statistics are right nothing moves. Following Blessing (1997) J. Appl. Cryst. 30, 421-426. GPU path: the refitted means are uploaded (SetEmMean). Measured (SHELXL R1(>4sigma) against the COD model, sm-a's harness; rc174 scaling): YAG 0.556 -> 0.127 (XDS 0.083 merged), rejected 8055 -> 53, normalised deviations calibrated (median |z| 0.62-0.69 in every intensity decile); aspirin 20 keV 0.0964 -> 0.0958; citric acid 0.161 -> 0.159; HEPES 0.0903 -> 0.0899; L-cystine 25 keV 0.1425 -> 0.1456; aspirin 25 keV 0.094 -> 0.106 (fixed-model R1 0.107 -> 0.173): its strong equivalents split into two frame-dependent populations from the per-frame partial scaling (sm-a's dq-smallmol), which the old under-sized sigmas happened to cut; with that scaling fixed (f69339ce6 + pooling, --no-scale-partials) this change is neutral to better on every small molecule (aspirin 20 .0618 -> .0586, aspirin 25 .0456 -> .0454, citric .1093 -> .1006, HEPES .0703 -> .0697, YAG .649 -> .222; SHELXL GooF ~1.1). => ship together with the scaling fix. Proteins (GPU full runs): CC1/2 and R_meas unchanged to 0.002; ISa myob 9.06 -> 8.20, thau 52.5 -> 47.5, cytc 25.8 -> 25.4, lyso 29.4 -> 29.4 (still above XDS's 5.2 / 44.5 / 31.8 / 28.3 except cytc). CPU build gives the same statistics as the GPU build on aspirin 20 keV and myob. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- .../scale_merge/RotationScaleMerge.cpp | 111 +++++++++++++++++- .../scale_merge/RotationScaleMergeGPU.cu | 8 ++ .../scale_merge/RotationScaleMergeGPU.h | 3 + 3 files changed, 118 insertions(+), 4 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 81dabcd22..edccd2d8c 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -4448,6 +4448,104 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool resid * resid / factor, o.d}); } } + fit_error_model(samples); + + // ---- The error model refitted about the mean the merge takes. ---- + // The fit above is made about a mean that weighs each full by its COUNTING variance. Where the + // equivalents of a reflection agree to within counting statistics that is the merge's own centre, + // and the model it gives is the model the merge applies. Where they do not - a systematic error + // the counting variance does not know about - it is the wrong centre, and in the worst way: the + // full with the fewest counts has the smallest variance, so the mean, and every deviation the + // model is fitted on, is pulled toward the observations that read lowest. Measured on a strongly + // absorbing crystal whose equivalents spread over a factor 100: the mean of (0 4 0) sat at 21k + // among observations from 11k to 1.9M, the deviations of the rest read as hundreds of counting + // sigmas, a ran into its bound and the six-sigma test about the median taken with the same + // weights deleted 63% of the observations - the strong ones (SHELXL R1 0.50 on that merge). + // So refit about the mean the MERGE takes - each full weighted by the variance the fitted model + // gives it, evaluated at the reflection's mean - until the model stops moving; the outlier + // test's median below takes the same weights. Where counting statistics were right the first fit + // is the answer and nothing moves; where b dominates, the weights even out, as the model says. + // Following Blessing (1997) J. Appl. Cryst. 30, 421-426: the centre an outlier is judged from + // is weighted as the merge weighs. + if (error_model_active) { + // Fulls in group order, so each pass over a group's members is one contiguous walk and the + // groups split across threads without every thread reading every full. + std::vector gstart(n_groups + 1, 0); + for (int i = 0; i < n_full; ++i) + if (mf.group[i] >= 0 && cnt[mf.group[i]] >= 2) ++gstart[mf.group[i] + 1]; + for (int g = 0; g < n_groups; ++g) gstart[g + 1] += gstart[g]; + std::vector member(gstart[n_groups]); + { + std::vector fill(gstart.begin(), gstart.end() - 1); + for (int i = 0; i < n_full; ++i) + if (mf.group[i] >= 0 && cnt[mf.group[i]] >= 2) member[fill[mf.group[i]]++] = i; + } + std::vector next_mean; + std::vector next; + constexpr int REFIT_ROUNDS = 10; + for (int round = 0; round < REFIT_ROUNDS; ++round) { + const double a = error_model_a, b = error_model_b; + const auto model_var = [&](int i, double mean) { + const double sc = static_cast(mf.sigma[i]) * mf.corr[i]; + return a * counting_variance(fulls[i], mean, sc * sc) + (b * mean) * (b * mean); + }; + std::vector> part(ThreadsForWork(member.size(), nthreads)); + const int nt = static_cast(part.size()); + next_mean.assign(n_groups, NAN); + // Each group's mean and samples come from its own members in index order, so the result + // does not depend on how the groups were split. + ParallelFor(nt, nt, [&](int t) { + const int g0 = static_cast(static_cast(n_groups) * t / nt); + const int g1 = static_cast(static_cast(n_groups) * (t + 1) / nt); + for (int g = g0; g < g1; ++g) { + if (gstart[g + 1] - gstart[g] < 2 || !std::isfinite(em_mean[g])) continue; + double sw = 0.0, swI = 0.0, swh[2] = {0.0, 0.0}, swIh[2] = {0.0, 0.0}; + int nh[2] = {0, 0}; + for (int q = gstart[g]; q < gstart[g + 1]; ++q) { + const int i = member[q]; + const double v = model_var(i, em_mean[g]); + if (!(v > 0.0)) continue; + const double I_corr = static_cast(mf.I[i]) * mf.corr[i]; + sw += 1.0 / v; swI += I_corr / v; + if (merge_friedel && group_has_hands[g]) { + swh[obs_hand[i]] += 1.0 / v; swIh[obs_hand[i]] += I_corr / v; nh[obs_hand[i]]++; + } + } + if (!(sw > 0.0)) continue; + const double mean = swI / sw; + next_mean[g] = mean; + // The samples about the hand's own mean where it has two of its own, as above. + for (int q = gstart[g]; q < gstart[g + 1]; ++q) { + const int i = member[q]; + const double v = model_var(i, em_mean[g]); + if (!(v > 0.0)) continue; + const int hh = obs_hand[i]; + const bool on_hand = merge_friedel && group_has_hands[g] && nh[hh] >= 2 && swh[hh] > 0.0; + const double centre = on_hand ? swIh[hh] / swh[hh] : mean; + const double factor = 1.0 - (1.0 / v) / (on_hand ? swh[hh] : sw); + if (factor < 0.05) continue; + const double sc = static_cast(mf.sigma[i]) * mf.corr[i]; + const double resid = static_cast(mf.I[i]) * mf.corr[i] - centre; + part[t].push_back({counting_variance(fulls[i], mean, sc * sc), centre * centre, + resid * resid / factor, mf.d[i]}); + } + } + }); + for (int g = 0; g < n_groups; ++g) + if (!std::isfinite(next_mean[g])) next_mean[g] = em_mean[g]; + next.clear(); + for (auto &v : part) next.insert(next.end(), v.begin(), v.end()); + em_mean.swap(next_mean); + samples.swap(next); + fit_error_model(samples); + if (std::fabs(error_model_a - a) <= 1e-3 * a + && std::fabs(error_model_b - b) <= 1e-3 * std::max(b, 1e-6)) + break; + } +#ifdef JFJOCH_USE_CUDA + if (use_gpu_merge) gpu_->SetEmMean(em_mean.data()); +#endif + } // Per-group outlier-rejection median of I*corr (host both paths - a per-group median is awkward on // the GPU; cheap here, cnt >= 3 filter from the em pass). Fed to the merge accumulate. @@ -4510,15 +4608,21 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // Inverse-variance WEIGHTED median: a frame the crystal barely diffracted on is scaled up by // 1/G together with its sigma, so a plain median lets two such observations outvote one // well-measured one - and the test below then rejects the well-measured one against its - // own small sigma. - std::vector> iv(start[n_sets]); // (I*corr, 1/(sigma*corr)^2) + // own small sigma. The variance is the model's, at the reflection's mean - the weight the + // merge gives the observation (see the centre refit above). + std::vector> iv(start[n_sets]); // (I*corr, 1/model variance) { std::vector fill(start.begin(), start.end() - 1); for (int i = 0; i < n_full; ++i) { const int g = mf.group[i]; if (g < 0) continue; const float sc = mf.sigma[i] * mf.corr[i]; - const std::pair v{mf.I[i] * mf.corr[i], sc > 0.0f ? 1.0f / (sc * sc) : 0.0f}; + const double mean = std::isfinite(em_mean[g]) ? em_mean[g] : static_cast(mf.I[i]) * mf.corr[i]; + const double mv = error_model_active + ? error_model_a * counting_variance(fulls[i], mean, static_cast(sc) * sc) + + (error_model_b * mean) * (error_model_b * mean) + : static_cast(sc) * sc; + const std::pair v{mf.I[i] * mf.corr[i], mv > 0.0 ? static_cast(1.0 / mv) : 0.0f}; if (cnt[g] >= 3) iv[fill[g]++] = v; if (!pair_needed.empty() && pair_needed[pair_of_group[g]]) iv[fill[n_groups + pair_of_group[g]]++] = v; @@ -4545,7 +4649,6 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool reject_median[g] = cnt[g] >= 3 ? set_median[g] : !pair_needed.empty() ? set_median[n_groups + pair_of_group[g]] : NAN; } - fit_error_model(samples); } // The full's sigma under the error model, with the variance evaluated at intensity I_for_b. diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 8a0e97f85..60a2e8d76 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -1110,6 +1110,14 @@ void RotationScaleMergeGPU::SetFrameCellOk(const uint8_t *frame_cell_ok) { // The per-group inv-var mean (em_mean) + the per-full leverage-corrected error-model samples over the // resident+scaled fulls. Stashes the filter context for the later MergeAccum/MergeRmeas calls. +void RotationScaleMergeGPU::SetEmMean(const double *em_mean) { + DeviceGuard guard(impl_->device, impl_->available); + auto &d = *impl_; + if (d.n_groups > 0) + CopyAndWait(d.m_em_mean.get(), em_mean, size_t(d.n_groups) * sizeof(double), cudaMemcpyHostToDevice, + impl_->s(), "ul em_mean"); +} + void RotationScaleMergeGPU::MergeEmSamples(bool for_search, double min_partiality, const uint8_t *hand, const uint8_t *has_hands, double *em_mean_out, int32_t *cnt_out, double *s2_out, diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index 766e2f640..17b56c520 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -103,6 +103,9 @@ public: // half-set weights multiplied by it. Requires MergeEmSamples first (em_mean resident). // reject_var_add (n_groups) widens the pooled cut by the shell's own measured Bijvoet // variance; null leaves the plain n-sigma test. + // Replace the per-group means MergeEmSamples left on the device (n_groups values): the merge's + // model sigmas are evaluated at them. + void SetEmMean(const double *em_mean); void MergeAccum(double error_model_a, double error_model_b, bool error_model_active, bool reject_outliers, double reject_nsigma, const float *reject_median, const float *reject_var_add, From 8394b5988d8a067b5ad001c3e38401d084b86132 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 4 Oct 2026 02:02:15 +0200 Subject: [PATCH 10/13] Integration follows the measured spot footprint where a spot outgrows the r1 disk The integrator's r1 disk and r2..r3 background ring are fixed in pixels and chosen from spots near the beam. On small-molecule data at 20-25 keV a spot's standard deviation grows from ~1 px near the beam to ~5 px at the edge (radially from parallax/obliquity, tangentially from the crystal's azimuthal spread), so the r1 = 4 disk holds a quarter of the flux there, the background ring a third of it, and the in-disk second moments the Gaussian is built from saturate near r1^2/4. On top of that, the profile/summation runaway guard sent 20-30% of these reflections - the strong, wide ones - back to the truncated r1 box sum. - SpotFootprint: every pre-scan spot (width frames) is measured with a window that follows it (3 sigma, iterated, re-centred), radially and tangentially; the medians per distance-from-beam bin become BraggIntegrationSettings::Footprint. Installed only where some bin outgrows r1, and on the adaptive side like the radius (pre-pass without; the starvation guard falls back to the settings without it). - BraggStencil: where 3 sigma > r1 the background ring starts at 3 sigma along and across the radius, the summation region is the r1 disk plus the 3-sigma footprint ellipse (so the guard's fallback is a complete intensity), and the per-reflection Gaussian takes the footprint widths. Compact spots keep the stencil bit for bit. Both engines build it from the same header. SHELXL against COD (R1 / fixed-XDS-model R1(F)): citric acid .101/.230 -> .077/.055, HEPES .070/.179 -> .048/.050, aspirin 20 keV .059/.070 -> .052/.061, aspirin 25 keV unchanged, L-cystine 25 keV unchanged (.145 -> .144). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- common/BraggIntegrationSettings.cpp | 13 ++ common/BraggIntegrationSettings.h | 16 +++ .../BraggIntegrationEngine.cpp | 15 ++- .../BraggIntegrationEngineCPU.cpp | 26 +++- .../BraggIntegrationEngineGPU.cu | 24 +++- .../bragg_integration/BraggStencil.h | 103 +++++++++++++- .../bragg_integration/CMakeLists.txt | 2 + .../bragg_integration/SpotFootprint.cpp | 126 ++++++++++++++++++ .../bragg_integration/SpotFootprint.h | 53 ++++++++ rugnux/Rugnux.cpp | 68 +++++++++- tests/BraggStencilTest.cpp | 67 ++++++++++ tests/CMakeLists.txt | 2 + tests/SpotFootprintTest.cpp | 68 ++++++++++ 13 files changed, 560 insertions(+), 23 deletions(-) create mode 100644 image_analysis/bragg_integration/SpotFootprint.cpp create mode 100644 image_analysis/bragg_integration/SpotFootprint.h create mode 100644 tests/SpotFootprintTest.cpp diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index ae8086b26..03b1a92bf 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -204,3 +204,16 @@ BraggIntegrationSettings &BraggIntegrationSettings::FlightPath(FlightPathMedium FlightPathMedium BraggIntegrationSettings::GetFlightPath() const { return flight_path; } + +BraggIntegrationSettings &BraggIntegrationSettings::Footprint(const SpotFootprint &input) { + if (input.sigma_rad.size() != input.sigma_tan.size() + || input.sigma_rad.size() > static_cast(SpotFootprint::MAX_BINS) + || (!input.empty() && !(input.bin_px > 0.0f))) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Invalid spot footprint table"); + footprint = input; + return *this; +} + +const SpotFootprint &BraggIntegrationSettings::GetFootprint() const { + return footprint; +} diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index 4e4b5645e..3555ccbdc 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -4,6 +4,7 @@ #pragma once #include +#include // Spot-intensity extraction method used by the Bragg integration engine. ProfileGaussian (default) // profile-fits with a measured-width Gaussian (Kabsch-style) - more accurate intensities than the @@ -34,6 +35,18 @@ enum class FlightPathMedium { Air, Helium, Vacuum }; // the max_hkl default in broker/jfjoch_api.yaml, so an omitting client and an omitting config agree. constexpr int BRAGG_ONLINE_DEFAULT_MAX_HKL = 100; +// The spot footprint measured on the data (image_analysis/bragg_integration/SpotFootprint.h): the +// radial and tangential standard deviations of a spot [px], tabulated against the distance from the +// beam centre in bins of bin_px. Empty until measured; the integrator then uses the fixed r1/r2/r3 +// stencil alone. +struct SpotFootprint { + static constexpr int MAX_BINS = 16; + float bin_px = 0.0f; + std::vector sigma_rad; + std::vector sigma_tan; + [[nodiscard]] bool empty() const { return sigma_rad.empty(); } +}; + class BraggIntegrationSettings { IntegratorMode integrator_mode = IntegratorMode::ProfileGaussian; float r_1 = 4; @@ -148,6 +161,7 @@ class BraggIntegrationSettings { // threshold fitted between two points rather than a physical criterion. The honest design is an // explicit option, a stated default, and the assumption printed in the report. FlightPathMedium flight_path = FlightPathMedium::Air; + SpotFootprint footprint; public: BraggIntegrationSettings& R1(float input); @@ -165,6 +179,7 @@ public: BraggIntegrationSettings& Overlap(OverlapMode input); BraggIntegrationSettings& OverlapMinPeak(float input); BraggIntegrationSettings& FlightPath(FlightPathMedium input); + BraggIntegrationSettings& Footprint(const SpotFootprint &input); [[nodiscard]] IntegratorMode GetIntegrator() const; @@ -186,4 +201,5 @@ public: [[nodiscard]] float GetOverlapMinPeak() const; // See flight_path: air unless the user said otherwise, which no file can say for them. [[nodiscard]] FlightPathMedium GetFlightPath() const; + [[nodiscard]] const SpotFootprint &GetFootprint() const; }; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index a0a8b263d..5a6e85db0 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -96,6 +96,19 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe stencil.bw_sigma = static_cast(bw_sigma); stencil.k_sigma = settings.GetStencilKSigma(); stencil.max_grow = bragg_engine::MAX_STENCIL_GROW_OVER_R3 * r3; + // The measured spot footprint (SpotFootprint.h). It moves the ring only where a spot outgrows the + // r1 disk, so a pattern of compact spots integrates exactly as without it. + static_assert(SpotFootprint::MAX_BINS <= BRAGG_FOOTPRINT_MAX_BINS); + stencil.r1 = settings.GetR1(); + const auto &fp = settings.GetFootprint(); + if (!empirical && !fp.empty()) { + stencil.fp_n = static_cast(fp.sigma_rad.size()); + stencil.fp_bin_px = fp.bin_px; + for (int i = 0; i < stencil.fp_n; ++i) { + stencil.fp_sigma_rad[i] = fp.sigma_rad[i]; + stencil.fp_sigma_tan[i] = fp.sigma_tan[i]; + } + } // Robust background ring, one estimator or the other (see BraggIntegrationSettings): a high-side // sigma-clip (rugnux --background-clip, the default) or, when the clip is switched off, a @@ -137,7 +150,7 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe r_max = std::hypot(std::max(beam_x, static_cast(xpixel) - beam_x), std::max(beam_y, static_cast(ypixel) - beam_y)); bkg_radial_built = bkg_radial || bkg_radial_auto; - const float grow_max = bkg_radial_built ? BraggStencilGrow_px(static_cast(r_max), stencil) + const float grow_max = bkg_radial_built ? BraggStencilMaxGrow_px(static_cast(r_max), stencil) : 0.0f; n_kern = static_cast(std::lround(grow_max)) + 1; // Every row must fit: the last one is built at grow = n_kern - 1, which rounding can put just diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 9c5d509d9..f7a46ac8d 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -224,7 +224,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, for (int x = x0; x <= x1; ++x) { const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y); const int32_t px = img[y * W + x]; - if (d.signal < r1_sq) { + if (BraggStencilInSignal(st, x - r.predicted_x, y - r.predicted_y, d.signal, r1_sq)) { // A pixel a nearer neighbour owns carries that neighbour's flux, so in Exclude // mode it leaves the disk entirely - the sum, the pixel count the background is // subtracted with, and the all-or-nothing validity gate alike. The box sum only @@ -519,11 +519,12 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, int Rf = R; const std::vector *Pvec = &shell_P[sh]; + const BraggStencil st = MakeBraggStencil(predicted[i].predicted_x, predicted[i].predicted_y, stencil); if (!empirical) { const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y; const double Rpx = std::hypot(rx, ry); const double tan2t = Rpx / F_px; - const double s2t = shell_sigma2[sh].tan; + double s2t = shell_sigma2[sh].tan; double s2r = s2t, ux = 1.0, uy = 0.0; bool elong = false; if (use_ellipse) { @@ -539,10 +540,18 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, elong = true; } } + // A spot that outgrows the r1 disk has widths the disk cannot measure - the moments above + // saturate near r1^2/4 - so the measured footprint sets them instead (SpotFootprint.h). + if (st.fp_s2t > 0.0f) { + ux = st.ux; uy = st.uy; + s2t = std::max(s2t, st.fp_s2t); + s2r = std::max(s2r, st.fp_s2r); + elong = true; + } // Build the Gaussian per reflection, centred on the sub-pixel predicted position and (when - // needed) radially elongated, on a grid grown to hold the streak. + // needed) elongated, on a grid grown to hold the wider of its two axes. const double fx = predicted[i].predicted_x - rh.cx, fy = predicted[i].predicted_y - rh.cy; - Rf = elong ? std::min(3 * R, static_cast(std::ceil(r2 + 2.0 * std::sqrt(s2r)))) : R; + Rf = elong ? std::min(3 * R, static_cast(std::ceil(r2 + 2.0 * std::sqrt(std::max(s2r, s2t))))) : R; // The grid's image rows, asked for now (as in pass A): building the profile below takes // long enough for them to arrive before the fit reads them. const int gx0 = std::max(0, rh.cx - Rf), gx1 = std::min(W - 1, rh.cx + Rf); @@ -582,7 +591,14 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (Pp <= 0.0) continue; const int x = rh.cx + dx, y = rh.cy + dy; if (x < 0 || y < 0 || x >= W || y >= H) continue; - const bool in_disk = dx * dx + dy * dy < r1_sq; + // The region the summation seed was taken over (pass A), so the guard below + // compares the two over the same pixels. + const bool in_disk = st.fp_s2t > 0.0f + ? BraggStencilInSignal(st, x - predicted[i].predicted_x, y - predicted[i].predicted_y, + (x - predicted[i].predicted_x) * (x - predicted[i].predicted_x) + + (y - predicted[i].predicted_y) * (y - predicted[i].predicted_y), + r1_sq) + : dx * dx + dy * dy < r1_sq; p_grid += Pp; p_peak = std::max(p_peak, Pp); if (in_disk) m_all += Pp; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 4ff0d148b..3f4e3d4ca 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -220,7 +220,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, // The window is the bounding box of the outer ellipse, so nearly half of it is neither the // signal disk nor the background ring. Reading the pixel only once it is known to be wanted // keeps those slots from fetching a cache line for nothing. - if (d.signal < p.r1_sq) { + if (BraggStencilInSignal(st, (float) x - cx, (float) y - cy, d.signal, p.r1_sq)) { // A pixel a nearer neighbour owns carries that neighbour's flux; see the CPU engine. ++l_nd; if (p.overlap) { @@ -625,6 +625,7 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x if (!ok_a[i]) { if (threadIdx.x == 0) ok_o[i] = 0; return; } const int cx = cx_a[i], cy = cy_a[i]; + const BraggStencil st = MakeBraggStencil(px_x[i], px_y[i], p.stencil); // As in learn_profile: one thread per block, not one per thread. Two double-precision divisions // on a card whose double throughput is a sixty-fourth of its single is worth doing once. __shared__ int s_sh; @@ -642,7 +643,7 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x __syncthreads(); } else { const int si = use_shell ? sh : N_SHELL; // N_SHELL is the global slot - const float s2t = sigma2_t[si]; + float s2t = sigma2_t[si]; const float rx = px_x[i] - p.beam_x, ry = px_y[i] - p.beam_y; const float Rpx = sqrtf(rx * rx + ry * ry); const float tan2t = Rpx / p.F_px; @@ -655,7 +656,14 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x const float radial_extra = fmaxf(sigma2_r[si] - s2t, sbw * sbw + p.c_radial * tan2t * tan2t); if (Rpx > 1e-6f && radial_extra > 0.25f) { ux = rx / Rpx; uy = ry / Rpx; s2r = s2t + radial_extra; elong = true; } } - const int Rf = elong ? min(3 * p.R, (int) ceilf(p.r2 + 2.0f * sqrtf(s2r))) : p.R; + // A spot that outgrows the r1 disk takes the measured footprint's widths (see the CPU engine). + if (st.fp_s2t > 0.0f) { + ux = st.ux; uy = st.uy; + s2t = fmaxf(s2t, st.fp_s2t); + s2r = fmaxf(s2r, st.fp_s2r); + elong = true; + } + const int Rf = elong ? min(3 * p.R, (int) ceilf(p.r2 + 2.0f * sqrtf(fmaxf(s2r, s2t)))) : p.R; const int Gf = 2 * Rf + 1; if (threadIdx.x == 0) { s_Rf = Rf; s_Gf = Gf; } __syncthreads(); @@ -697,7 +705,11 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x const int dx = k % Gf - Rf, dy = k / Gf - Rf; const int x = cx + dx, y = cy + dy; if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue; - const bool in_disk = (float) (dx * dx + dy * dy) < p.r1_sq; + // The region the summation seed was taken over (boxsum), as on the CPU. + const float ddx = (float) x - px_x[i], ddy = (float) y - px_y[i]; + const bool in_disk = st.fp_s2t > 0.0f + ? BraggStencilInSignal(st, ddx, ddy, ddx * ddx + ddy * ddy, p.r1_sq) + : (float) (dx * dx + dy * dy) < p.r1_sq; l_pgrid += Pp; l_ppeak = fmaxf(l_ppeak, Pp); if (in_disk) l_mall += Pp; @@ -819,11 +831,11 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment // Shared radial window of the background curve in boxsum. The stencil spans r0 +- the radial // semi-axis of the outer ellipse, widest at the far corner of the detector, and the window // covers twice that so a reflection anywhere on the detector fits. - rad_w = 2 * (static_cast(std::ceil(r3 + BraggStencilGrow_px(static_cast(r_max), stencil))) + 1) + 1; + rad_w = 2 * (static_cast(std::ceil(r3 + BraggStencilMaxGrow_px(static_cast(r_max), stencil))) + 1) + 1; // The pixels one reflection can mark: the box mark_mask walks, at the widest aperture on the // detector. Run() weighs that against the frame to decide how to clear the mask afterwards. - const int mark_half = static_cast(std::ceil(r2 + BraggStencilGrow_px(static_cast(r_max), stencil))) + 1; + const int mark_half = static_cast(std::ceil(r2 + BraggStencilMaxGrow_px(static_cast(r_max), stencil))) + 1; mask_box_px = static_cast(2 * mark_half + 1) * static_cast(2 * mark_half + 1); boxsum_shared_bytes = static_cast(rad_w) * (sizeof(unsigned long long) + sizeof(int)); diff --git a/image_analysis/bragg_integration/BraggStencil.h b/image_analysis/bragg_integration/BraggStencil.h index d85866f2f..bc49cbd06 100644 --- a/image_analysis/bragg_integration/BraggStencil.h +++ b/image_analysis/bragg_integration/BraggStencil.h @@ -65,13 +65,27 @@ #define BRAGG_STENCIL_HD inline #endif +// The measured spot footprint (SpotFootprint.h): the radial and tangential standard deviations of a +// spot [px], tabulated against distance from the beam centre in bins of fp_bin_px. Held here, by value, +// because both engines build every stencil from these parameters - the GPU gets the table with the +// kernel arguments and needs no buffer of its own. +constexpr int BRAGG_FOOTPRINT_MAX_BINS = 16; +// A spot reaches about this many of its standard deviations: the background ring starts there and the +// profile is fitted out to it. +constexpr float BRAGG_FOOTPRINT_NSIGMA = 3.0f; + // Fixed per-experiment inputs to the stencil law (mirrors BraggIntegrationEngine's members). struct BraggStencilParams { float beam_x = 0.0f, beam_y = 0.0f; + float r1 = 4.0f; float r2 = 6.0f, r3 = 10.0f; float bw_sigma = 0.0f; // radial streak per pixel of radius (bandwidth sigma, dimensionless) float k_sigma = 0.0f; // radial sigmas to push the ring out by; 0 = the old circular stencil float max_grow = 0.0f; // hard cap on the radial growth [px]; <= 0 means uncapped + int fp_n = 0; // bins of the measured footprint; 0 = none measured + float fp_bin_px = 0.0f; + float fp_sigma_rad[BRAGG_FOOTPRINT_MAX_BINS] = {}; + float fp_sigma_tan[BRAGG_FOOTPRINT_MAX_BINS] = {}; }; // One reflection's stencil, in its own radial/tangential frame. @@ -79,7 +93,15 @@ struct BraggStencil { float ux = 1.0f, uy = 0.0f; // unit vector beam -> reflection float r0 = 0.0f; // distance from the beam centre [px] float grow = 0.0f; // radial growth of the ring [px] (0 = circular) + float grow_tan = 0.0f; // tangential growth of the ring [px] (footprint only) float q_in = 0.0f, q_out = 0.0f; // radial shrink coefficients (0 = circular) + float qt_in = 0.0f, qt_out = 0.0f; // tangential shrink coefficients (0 = circular) + // The measured footprint where the spot outgrows the r1 disk (variances [px^2]), else 0: the + // profile fit takes its Gaussian at least this wide. + float fp_s2r = 0.0f, fp_s2t = 0.0f; + // 1 / (BRAGG_FOOTPRINT_NSIGMA * sigma)^2 along and across the radius: the footprint ellipse a wide + // spot is SUMMED over, beside the r1 disk (BraggStencilInSignal). + float core_r = 0.0f, core_t = 0.0f; float ex_in = 0.0f, ey_in = 0.0f; // axis-aligned half-extent of the inner (r2) ellipse float ex_out = 0.0f, ey_out = 0.0f; // axis-aligned half-extent of the outer (r3) ellipse }; @@ -94,6 +116,38 @@ BRAGG_STENCIL_HD float BraggStencilGrow_px(float Rpx, const BraggStencilParams & return (p.max_grow > 0.0f && grow > p.max_grow) ? p.max_grow : grow; } +// The footprint at Rpx: linear between bin centres, flat beyond the first and the last. +BRAGG_STENCIL_HD void BraggFootprintAt(float Rpx, const BraggStencilParams &p, float &s_rad, float &s_tan) { + const float t = Rpx / p.fp_bin_px - 0.5f; + int i = (int) floorf(t); + float f = t - (float) i; + if (i < 0) { i = 0; f = 0.0f; } + if (i >= p.fp_n - 1) { i = p.fp_n - 1; f = 0.0f; } + const int j = i + 1 < p.fp_n ? i + 1 : i; + s_rad = p.fp_sigma_rad[i] + f * (p.fp_sigma_rad[j] - p.fp_sigma_rad[i]); + s_tan = p.fp_sigma_tan[i] + f * (p.fp_sigma_tan[j] - p.fp_sigma_tan[i]); +} + +// Growth of the ring beyond r2 that the footprint asks for along one axis: none while the spot fits +// the r1 disk - there the disk resolves the spot and the stencil stays what it always was - and +// otherwise out to BRAGG_FOOTPRINT_NSIGMA of it, capped as the bandwidth growth is. +BRAGG_STENCIL_HD float BraggFootprintGrow_px(float sigma, const BraggStencilParams &p) { + const float g = BRAGG_FOOTPRINT_NSIGMA * sigma - p.r2; + if (!(g > 0.0f)) return 0.0f; + return (p.max_grow > 0.0f && g > p.max_grow) ? p.max_grow : g; +} + +// The largest growth any reflection on a detector reaching r_max can get, radially - what the +// buffers that hold a ring are sized by. +BRAGG_STENCIL_HD float BraggStencilMaxGrow_px(float r_max, const BraggStencilParams &p) { + float g = BraggStencilGrow_px(r_max, p); + for (int i = 0; i < p.fp_n; ++i) { + const float s = fmaxf(p.fp_sigma_rad[i], p.fp_sigma_tan[i]); + if (BRAGG_FOOTPRINT_NSIGMA * s > p.r1) g = fmaxf(g, BraggFootprintGrow_px(s, p)); + } + return g; +} + BRAGG_STENCIL_HD BraggStencil MakeBraggStencil(float px_x, float px_y, const BraggStencilParams &p) { BraggStencil s; const float rx = px_x - p.beam_x, ry = px_y - p.beam_y; @@ -105,24 +159,44 @@ BRAGG_STENCIL_HD BraggStencil MakeBraggStencil(float px_x, float px_y, const Bra } s.grow = BraggStencilGrow_px(r0, p); - const float grow = s.grow; - if (grow == 0.0f) { // the common case, and the only one the online path can reach + // A spot wider than the r1 disk can resolve keeps its background ring clear of itself: the ring + // starts BRAGG_FOOTPRINT_NSIGMA of the measured footprint out, radially and tangentially apart. + if (p.fp_n > 0) { + float s_rad, s_tan; + BraggFootprintAt(r0, p, s_rad, s_tan); + if (BRAGG_FOOTPRINT_NSIGMA * fmaxf(s_rad, s_tan) > p.r1) { + s.fp_s2r = s_rad * s_rad; + s.fp_s2t = s_tan * s_tan; + s.core_r = 1.0f / (BRAGG_FOOTPRINT_NSIGMA * BRAGG_FOOTPRINT_NSIGMA * s.fp_s2r); + s.core_t = 1.0f / (BRAGG_FOOTPRINT_NSIGMA * BRAGG_FOOTPRINT_NSIGMA * s.fp_s2t); + s.grow = fmaxf(s.grow, BraggFootprintGrow_px(s_rad, p)); + s.grow_tan = BraggFootprintGrow_px(s_tan, p); + } + } + const float grow = s.grow, grow_t = s.grow_tan; + if (grow == 0.0f && grow_t == 0.0f) { // the common case, and the only one the online path can reach s.ex_in = p.r2; s.ey_in = p.r2; s.ex_out = p.r3; s.ey_out = p.r3; return s; } - const float a_in = p.r2 + grow, a_out = p.r3 + grow; + const float a_in = p.r2 + grow, a_out = p.r3 + grow; // semi-axes along the radius + const float b_in = p.r2 + grow_t, b_out = p.r3 + grow_t; // and across it const float si = p.r2 / a_in, so = p.r3 / a_out; s.q_in = 1.0f - si * si; s.q_out = 1.0f - so * so; + if (grow_t > 0.0f) { + const float ti = p.r2 / b_in, to = p.r3 / b_out; + s.qt_in = 1.0f - ti * ti; + s.qt_out = 1.0f - to * to; + } // Axis-aligned half-extents of an ellipse with semi-axes (a along u, b across it). const float ux2 = s.ux * s.ux, uy2 = s.uy * s.uy; - s.ex_in = sqrtf(a_in * a_in * ux2 + p.r2 * p.r2 * uy2); - s.ey_in = sqrtf(a_in * a_in * uy2 + p.r2 * p.r2 * ux2); - s.ex_out = sqrtf(a_out * a_out * ux2 + p.r3 * p.r3 * uy2); - s.ey_out = sqrtf(a_out * a_out * uy2 + p.r3 * p.r3 * ux2); + s.ex_in = sqrtf(a_in * a_in * ux2 + b_in * b_in * uy2); + s.ey_in = sqrtf(a_in * a_in * uy2 + b_in * b_in * ux2); + s.ex_out = sqrtf(a_out * a_out * ux2 + b_out * b_out * uy2); + s.ey_out = sqrtf(a_out * a_out * uy2 + b_out * b_out * ux2); return s; } @@ -183,5 +257,20 @@ BRAGG_STENCIL_HD BraggStencilDist BraggStencilDistances(const BraggStencil &s, f const float rad2 = d.rad * d.rad; d.inner = d.signal - s.q_in * rad2; d.outer = d.signal - s.q_out * rad2; + if (s.qt_in != 0.0f) { // grown tangentially as well (the footprint); tan^2 = d2 - rad^2 + const float tan2 = d.signal - rad2; + d.inner -= s.qt_in * tan2; + d.outer -= s.qt_out * tan2; + } return d; } + +// The signal region: the r1 disk, and for a spot that outgrows it the whole footprint ellipse too, so +// the summation - the seed of the profile fit, and what the fit falls back to - holds the spot rather +// than its core. Without a footprint this is the plain disk test, bit for bit. +BRAGG_STENCIL_HD bool BraggStencilInSignal(const BraggStencil &s, float ddx, float ddy, float d2, float r1_sq) { + if (d2 < r1_sq) return true; + if (!(s.fp_s2t > 0.0f)) return false; + const float rad = ddx * s.ux + ddy * s.uy, tn = -ddx * s.uy + ddy * s.ux; + return rad * rad * s.core_r + tn * tn * s.core_t < 1.0f; +} diff --git a/image_analysis/bragg_integration/CMakeLists.txt b/image_analysis/bragg_integration/CMakeLists.txt index 7978c0439..aa8b3cd80 100644 --- a/image_analysis/bragg_integration/CMakeLists.txt +++ b/image_analysis/bragg_integration/CMakeLists.txt @@ -3,6 +3,8 @@ ADD_LIBRARY(JFJochBraggIntegration STATIC BraggIntegrationEngine.h BraggIntegrationEngineCPU.cpp BraggIntegrationEngineCPU.h + SpotFootprint.cpp + SpotFootprint.h Regression.h CalcISigma.cpp CalcISigma.h diff --git a/image_analysis/bragg_integration/SpotFootprint.cpp b/image_analysis/bragg_integration/SpotFootprint.cpp new file mode 100644 index 000000000..bfa6d0257 --- /dev/null +++ b/image_analysis/bragg_integration/SpotFootprint.cpp @@ -0,0 +1,126 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "SpotFootprint.h" + +#include +#include + +namespace { + +// The window around a spot is BRAGG_FOOTPRINT_NSIGMA-like: three of its standard deviations, at least +// a few pixels and at most this many, which is wider than any spot the integrator could hold. +constexpr float WINDOW_NSIGMA = 3.0f; +constexpr float WINDOW_MIN_PX = 3.0f; +constexpr float WINDOW_MAX_PX = 24.0f; +// The background is the median of an elliptical annulus between these multiples of the window. +constexpr float BKG_INNER = 1.5f; +constexpr float BKG_OUTER = 2.2f; +constexpr int ITERATIONS = 8; + +inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MAX; } + +float median_of(std::vector &v) { + const size_t m = v.size() / 2; + std::nth_element(v.begin(), v.begin() + static_cast(m), v.end()); + return v[m]; +} + +} // namespace + +void MeasureFootprintSpots(const int32_t *img, int width, int height, float beam_x, float beam_y, + const std::vector &x, const std::vector &y, + std::vector &out) { + const int half = static_cast(std::ceil(BKG_OUTER * WINDOW_MAX_PX)) + 1; + std::vector ring; + for (size_t s = 0; s < x.size(); ++s) { + const float rx = x[s] - beam_x, ry = y[s] - beam_y; + const float r = std::sqrt(rx * rx + ry * ry); + if (!(r > 1.0f)) continue; + const float ux = rx / r, uy = ry / r; + const int ix = static_cast(std::lround(x[s])), iy = static_cast(std::lround(y[s])); + if (ix - half < 0 || iy - half < 0 || ix + half >= width || iy + half >= height) continue; + + // Start from a compact spot at the prediction; each round re-centres on the signal and takes + // the window to three of the widths just measured. + float cx = x[s], cy = y[s]; + float s2r = 1.0f, s2t = 1.0f; + bool ok = true; + for (int it = 0; it < ITERATIONS && ok; ++it) { + const float wr = std::clamp(WINDOW_NSIGMA * std::sqrt(s2r), WINDOW_MIN_PX, WINDOW_MAX_PX); + const float wt = std::clamp(WINDOW_NSIGMA * std::sqrt(s2t), WINDOW_MIN_PX, WINDOW_MAX_PX); + const float reach = BKG_OUTER * std::max(wr, wt); + const int x0 = static_cast(std::floor(cx - reach)), x1 = static_cast(std::ceil(cx + reach)); + const int y0 = static_cast(std::floor(cy - reach)), y1 = static_cast(std::ceil(cy + reach)); + if (x0 < 0 || y0 < 0 || x1 >= width || y1 >= height) { ok = false; break; } + + ring.clear(); + for (int py = y0; py <= y1; ++py) + for (int px = x0; px <= x1; ++px) { + const float dx = px - cx, dy = py - cy; + const float rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; + const float e = rad * rad / (wr * wr) + tn * tn / (wt * wt); + const int32_t v = img[static_cast(py) * width + px]; + if (e >= BKG_INNER * BKG_INNER && e < BKG_OUTER * BKG_OUTER && valid(v)) + ring.push_back(static_cast(v)); + } + if (ring.size() < 10) { ok = false; break; } + const double bkg = median_of(ring); + + double w = 0.0, mr = 0.0, mt = 0.0, m2r = 0.0, m2t = 0.0; + for (int py = y0; py <= y1 && ok; ++py) + for (int px = x0; px <= x1; ++px) { + const float dx = px - cx, dy = py - cy; + const float rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; + if (rad * rad / (wr * wr) + tn * tn / (wt * wt) >= 1.0f) continue; + const int32_t v = img[static_cast(py) * width + px]; + if (!valid(v)) { ok = false; break; } + const double net = static_cast(v) - bkg; + w += net; + mr += net * rad; + mt += net * tn; + m2r += net * rad * rad; + m2t += net * tn * tn; + } + if (!ok || !(w > 0.0)) { ok = false; break; } + const double cr = mr / w, ct = mt / w; + s2r = static_cast(std::max(0.25, m2r / w - cr * cr)); + s2t = static_cast(std::max(0.25, m2t / w - ct * ct)); + cx += static_cast(cr * ux - ct * uy); + cy += static_cast(cr * uy + ct * ux); + } + if (ok) + out.push_back({r, std::sqrt(s2r), std::sqrt(s2t)}); + } +} + +SpotFootprint FootprintFromSpots(const std::vector &spots, float r_max) { + SpotFootprint fp; + if (!(r_max > 0.0f)) return fp; + const float bin = r_max / FOOTPRINT_BINS; + std::vector> rad(FOOTPRINT_BINS), tan(FOOTPRINT_BINS); + for (const auto &s : spots) { + const int b = std::clamp(static_cast(s.r_px / bin), 0, FOOTPRINT_BINS - 1); + rad[b].push_back(s.sigma_rad); + tan[b].push_back(s.sigma_tan); + } + std::vector filled; + std::vector mr(FOOTPRINT_BINS, 0.0f), mt(FOOTPRINT_BINS, 0.0f); + for (int b = 0; b < FOOTPRINT_BINS; ++b) + if (static_cast(rad[b].size()) >= FOOTPRINT_MIN_SPOTS_PER_BIN) { + mr[b] = median_of(rad[b]); + mt[b] = median_of(tan[b]); + filled.push_back(b); + } + if (filled.empty()) return fp; + fp.bin_px = bin; + for (int b = 0; b < FOOTPRINT_BINS; ++b) { + // The nearest bin that has enough spots; the inner one on a tie. + int best = filled.front(); + for (int f : filled) + if (std::abs(f - b) < std::abs(best - b)) best = f; + fp.sigma_rad.push_back(mr[best]); + fp.sigma_tan.push_back(mt[best]); + } + return fp; +} diff --git a/image_analysis/bragg_integration/SpotFootprint.h b/image_analysis/bragg_integration/SpotFootprint.h new file mode 100644 index 000000000..d973481a2 --- /dev/null +++ b/image_analysis/bragg_integration/SpotFootprint.h @@ -0,0 +1,53 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +// ============================================================================= +// SpotFootprint - how far a spot really reaches, measured on the data +// ============================================================================= +// +// The integrator's r1 disk and r2..r3 background ring are fixed in pixels, chosen from spots near the +// beam. Away from it a spot can grow several times wider: radially from the sensor's parallax and the +// obliquity of the incidence, tangentially from the crystal's own spread (an azimuthal mosaic spread +// rotates the diffracted beam about the incident one, smearing the spot along the ring). Measured on +// small-molecule data at 20 keV, the standard deviation goes from 1 px near the beam to 5 px at the +// edge, so the fixed disk holds a quarter of the flux there and the background ring a third of it - +// an intensity loss that grows with resolution and that no scaling can see, because symmetry mates +// share their resolution. +// +// The widths the integrator learns itself cannot follow: they are second moments taken inside the r1 +// disk, so they saturate near r1^2/4. Here each strong spot is measured with a window that follows +// the spot instead - three of its own standard deviations, iterated - separately along the radius and +// across it, and the widths are tabulated against the distance from the beam centre. The table goes +// to the integrator through BraggIntegrationSettings; where it says a spot outgrows the r1 disk, the +// background ring is moved clear of the spot and the profile is fitted at the measured width +// (BraggStencil.h). Compact spots leave the integration exactly as it was. +// ============================================================================= + +#include +#include + +#include "../../common/BraggIntegrationSettings.h" + +// One spot's measured widths, and its distance from the beam centre [px]. +struct FootprintSpot { + float r_px = 0.0f; + float sigma_rad = 0.0f; + float sigma_tan = 0.0f; +}; + +// Measure the spots at the given (sub-pixel) positions on one preprocessed frame (INT32_MIN masked, +// INT32_MAX saturated). A spot whose window holds an unreadable pixel, or no signal above its +// background, is left out. Results are appended in the order of the positions. +void MeasureFootprintSpots(const int32_t *img, int width, int height, float beam_x, float beam_y, + const std::vector &x, const std::vector &y, + std::vector &out); + +// Tabulate the spots: the median widths in bins of distance from the beam centre out to r_max. A bin +// holding too few spots takes the nearest bin that has enough; with no such bin the table is empty. +SpotFootprint FootprintFromSpots(const std::vector &spots, float r_max); + +// Spots per distance bin a table needs before it is believed. +constexpr int FOOTPRINT_MIN_SPOTS_PER_BIN = 20; +constexpr int FOOTPRINT_BINS = 12; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 720f41689..0800758b0 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -7,6 +7,8 @@ #include "WriteModel.h" #include "SpindleCuspLoss.h" #include "SpotWidth.h" +#include "../image_analysis/bragg_integration/SpotFootprint.h" +#include "../image_analysis/bragg_integration/BraggStencil.h" // BRAGG_FOOTPRINT_NSIGMA #include "HotPixels.h" #include "../image_analysis/SensorAbsorption.h" #include "DiagnosticOutput.h" @@ -833,6 +835,7 @@ std::optional Rugnux::RebinAndRefit(const CalibrationResult & return best; } + void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, RugnuxObserver *observer) { Logger logger("Rugnux"); @@ -975,7 +978,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const auto find_spots = [&](PreScanWorker &w, CompressedImage &image, int image_idx, bool for_beam_center, std::vector &out, bool for_width, std::vector &curves, - std::vector &spot_q) { + std::vector &spot_q, std::vector &footprint) { try { // As in the image loops: a frame the device cannot decode goes to the host decoder. ImageStatistics stats; @@ -1010,10 +1013,23 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru } } // The width goes first: it wants the whole spot list, to know which spots are isolated. - if (for_width) + if (for_width) { MeasureSpotFluxCurves(*w.preprocessed, static_cast(prescan_x.GetXPixelsNumConv()), static_cast(prescan_x.GetYPixelsNumConv()), prescan_x.GetDiffractionGeometry(), spots, curves); + // The footprint is measured on every spot, not only the isolated ones the width takes: + // it is the far, wide spots it is for (SpotFootprint.h). + std::vector fx, fy; + for (const auto &spot : spots) { + const Coord c = spot.RawCoord(); + fx.push_back(c.x); + fy.push_back(c.y); + } + const auto &g = prescan_x.GetDiffractionGeometry(); + MeasureFootprintSpots(w.preprocessed->data(), static_cast(prescan_x.GetXPixelsNumConv()), + static_cast(prescan_x.GetYPixelsNumConv()), + g.GetBeamX_pxl(), g.GetBeamY_pxl(), fx, fy, footprint); + } if (!for_beam_center) return; // The strongest of a crowded frame: the symmetry is over-determined either way, and the @@ -1037,6 +1053,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // The spots the bandwidth is measured on: the width's, and more where those are too few // (spot_width::BANDWIDTH_POOL_SPOTS). std::vector shape_curves; + // Every spot's measured widths, for the spot footprint (SpotFootprint.h). + std::vector footprint_spots; // Images the width ended up being measured on, for the log: the tiers below stop as soon as the // answer has settled, so this is a property of the crystal and worth reporting. size_t width_images = 0; @@ -1102,6 +1120,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru std::vector> spots_of(ordinals.size()); std::vector> curves_of(ordinals.size()); std::vector> spot_q_of(ordinals.size()); + std::vector> footprint_of(ordinals.size()); // A frame joins the pool only if it could be read, as it did when this was a serial loop. std::vector spot_read(ordinals.size(), 0); // Frames read for the bandwidth alone, past the tier the width settled on. Their spots stay @@ -1185,7 +1204,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru if (for_beam_center) spot_read[i] = 1; if (for_shape && !for_beam_center) shape_only[i] = 1; find_spots(w, w.raw_image.image, image_idx, for_beam_center, spots_of[i], - for_width || for_shape, curves_of[i], spot_q_of[i]); + for_width || for_shape, curves_of[i], spot_q_of[i], footprint_of[i]); } })); for (auto &f : futures) f.get(); @@ -1216,6 +1235,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru for (size_t i = 0; i < ordinals.size(); i++) if (!shape_only[i]) prescan_spot_q.insert(prescan_spot_q.end(), spot_q_of[i].begin(), spot_q_of[i].end()); + // In sample order, so the table does not depend on how the workers interleaved. + for (const auto &f : footprint_of) + footprint_spots.insert(footprint_spots.end(), f.begin(), f.end()); // Frame numbering follows the sample order, exactly as the serial read did. for (size_t i = 0; i < ordinals.size(); i++) { @@ -1278,6 +1300,43 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru "integration radii r1={:.1f} r2={:.1f} r3={:.2f}", *r80, spot_width::D_REF_A, width_curves.size(), width_images, r1, r2, r3); } + + // How far the spots reach away from the beam, where the radius chosen above at 5 A no + // longer holds them (SpotFootprint.h). It changes the integration only where a spot + // outgrows the r1 disk. + const auto &g = experiment_.GetDiffractionGeometry(); + const float bx = g.GetBeamX_pxl(), by = g.GetBeamY_pxl(); + const float W = static_cast(experiment_.GetXPixelsNum()); + const float H = static_cast(experiment_.GetYPixelsNum()); + const SpotFootprint fp = FootprintFromSpots(footprint_spots, + std::hypot(std::max(bx, W - bx), std::max(by, H - by))); + // It acts only where a spot outgrows the r1 disk (BraggStencil.h), so a pattern of compact + // spots is left on exactly the settings - and the passes - it had without it. + const float r1_now = experiment_.GetBraggIntegrationSettings().GetR1(); + bool outgrows = false; + for (size_t b = 0; b < fp.sigma_rad.size(); ++b) + outgrows |= BRAGG_FOOTPRINT_NSIGMA * std::max(fp.sigma_rad[b], fp.sigma_tan[b]) > r1_now; + if (fp.empty()) { + logger.Info("Spot footprint: too few spots to tabulate ({})", footprint_spots.size()); + } else if (!outgrows) { + logger.Info("Spot footprint from {} spots: every spot fits the r1={:.1f} disk", footprint_spots.size(), + r1_now); + } else { + std::string table; + for (size_t b = 0; b < fp.sigma_rad.size(); ++b) + table += fmt::format(" {:.0f}:{:.1f}/{:.1f}", (b + 0.5f) * fp.bin_px, + fp.sigma_rad[b], fp.sigma_tan[b]); + logger.Info("Spot footprint from {} spots, sigma radial/tangential [px] by distance from " + "the beam [px]:{}", footprint_spots.size(), table); + // On the adaptive side, like the radius: the geometry pre-pass integrates without it, + // and a canonical pass whose wider rings the neighbours starve falls back to the + // settings before it (bragg_before_adaptive_, see RunAllPasses). + BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); + if (!bragg_before_adaptive_) + bragg_before_adaptive_ = bis; + bis.Footprint(fp); + experiment_.ImportBraggIntegrationSettings(bis); + } } // The X-ray bandwidth, read off the shapes of the same spots (spot_width::EstimateBandwidth). A @@ -1605,8 +1664,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // The extra frames go to the beam centre alone; the width and the powder // rings have already been measured on the projection's own sample. std::vector unused_q; + std::vector unused_fp; find_spots(w, w.raw_image.image, image_idx, true, extra_spots[i], - false, width_curves, unused_q); + false, width_curves, unused_q, unused_fp); } } })); diff --git a/tests/BraggStencilTest.cpp b/tests/BraggStencilTest.cpp index 4d37aa86a..c6b9f9389 100644 --- a/tests/BraggStencilTest.cpp +++ b/tests/BraggStencilTest.cpp @@ -158,3 +158,70 @@ TEST_CASE("BraggStencil_KernelIndexInRange", "[Integration][portable]") { } } } + +namespace { + +// A footprint growing linearly from sigma 1 px at the beam to `edge` px at 800 px, radial and tangential +// alike unless told otherwise. +BraggStencilParams FootprintParams(float edge_rad, float edge_tan) { + BraggStencilParams p = Params(0.0f, 0.0f); + p.r1 = 4.0f; + p.fp_n = 8; + p.fp_bin_px = 100.0f; + for (int i = 0; i < p.fp_n; ++i) { + const float t = (i + 0.5f) / p.fp_n; + p.fp_sigma_rad[i] = 1.0f + t * (edge_rad - 1.0f); + p.fp_sigma_tan[i] = 1.0f + t * (edge_tan - 1.0f); + } + return p; +} + +} // namespace + +// Where the footprint fits the r1 disk (3 sigma <= r1) the stencil is the circular one, bit for bit: +// compact spots integrate exactly as without a footprint. +TEST_CASE("BraggStencil_FootprintInsideDiskIsInert", "[Integration]") { + const BraggStencilParams p = FootprintParams(1.3f, 1.3f); // 3 sigma < 4 everywhere + const BraggStencilParams none = Params(0.0f, 0.0f); + for (float py = 0.0f; py < 800.0f; py += 37.0f) + for (float px = 0.0f; px < 800.0f; px += 41.0f) { + const BraggStencil s = MakeBraggStencil(px, py, p); + const BraggStencil n = MakeBraggStencil(px, py, none); + REQUIRE(s.fp_s2r == 0.0f); + REQUIRE(s.grow == 0.0f); + REQUIRE(s.grow_tan == 0.0f); + for (int dy = -12; dy <= 12; ++dy) + for (int dx = -12; dx <= 12; ++dx) { + const auto d = BraggStencilDistances(s, static_cast(dx), static_cast(dy)); + const auto e = BraggStencilDistances(n, static_cast(dx), static_cast(dy)); + REQUIRE(d.inner == e.inner); + REQUIRE(d.outer == e.outer); + } + } +} + +// A spot wider than the disk pushes the ring to 3 sigma along each axis separately: a pixel at 3 sigma +// along the radius (or across it) is no longer background, one just beyond the grown ring's inner edge +// is. +TEST_CASE("BraggStencil_FootprintGrowsRingAlongEachAxis", "[Integration]") { + const BraggStencilParams p = FootprintParams(3.0f, 5.0f); + const float px = 400.0f + 700.0f, py = 400.0f; // on +x, 700 px out: radial = x, tangential = y + const BraggStencil s = MakeBraggStencil(px, py, p); + float sr, st; + BraggFootprintAt(700.0f, p, sr, st); + REQUIRE(s.fp_s2r == sr * sr); + REQUIRE(s.fp_s2t == st * st); + REQUIRE_THAT(s.grow, Catch::Matchers::WithinAbs(3.0f * sr - p.r2, 1e-4)); + REQUIRE_THAT(s.grow_tan, Catch::Matchers::WithinAbs(3.0f * st - p.r2, 1e-4)); + const float ar = p.r2 + s.grow, at = p.r2 + s.grow_tan; + // Just inside the inner ellipse along each axis: not background. + REQUIRE(BraggStencilDistances(s, 0.98f * ar, 0.0f).inner < p.r2 * p.r2); + REQUIRE(BraggStencilDistances(s, 0.0f, 0.98f * at).inner < p.r2 * p.r2); + // Just outside: background ring. + REQUIRE(BraggStencilDistances(s, 1.02f * ar, 0.0f).inner >= p.r2 * p.r2); + REQUIRE(BraggStencilDistances(s, 0.0f, 1.02f * at).inner >= p.r2 * p.r2); + // The bounding box holds the outer ellipse. + REQUIRE(s.ex_out >= p.r3 + s.grow - 1e-3f); + REQUIRE(s.ey_out >= p.r3 + s.grow_tan - 1e-3f); + REQUIRE(BraggStencilMaxGrow_px(1000.0f, p) >= std::max(s.grow, s.grow_tan)); +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c6d6df0bb..816808b1b 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -17,6 +17,7 @@ ADD_EXECUTABLE(jfjoch_portable_test EXCLUDE_FROM_ALL AzimuthalIntegrationTest.cpp CalcBraggPredictionTest.cpp BraggStencilTest.cpp + SpotFootprintTest.cpp BraggIntegrationEngineCompressedImageTest.cpp BraggIntegrationEngineCPUTest.cpp WriteReflectionsTest.cpp @@ -79,6 +80,7 @@ ADD_EXECUTABLE(jfjoch_test BraggIntegrationEngineCompressedImageTest.cpp BraggIntegrationEngineCPUTest.cpp BraggStencilTest.cpp + SpotFootprintTest.cpp LossyFilterTest.cpp ImageBufferTest.cpp PixelMaskTest.cpp diff --git a/tests/SpotFootprintTest.cpp b/tests/SpotFootprintTest.cpp new file mode 100644 index 000000000..46699b50e --- /dev/null +++ b/tests/SpotFootprintTest.cpp @@ -0,0 +1,68 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#include +#include + +#include "../image_analysis/bragg_integration/SpotFootprint.h" + +// Spots drawn as Gaussians elongated along and across the radius are measured back at their widths, +// whatever their azimuth, and tabulated by distance from the beam. +TEST_CASE("SpotFootprint_MeasuresRadialAndTangentialWidths", "[Integration][portable]") { + const int W = 1200, H = 1200; + const float bx = 600.0f, by = 600.0f; + std::vector img(static_cast(W) * H, 10); // flat background + std::vector xs, ys; + // Width grows with distance: sigma_rad = 1 + r/200, sigma_tan = 1 + r/100. + for (int k = 0; k < 400; ++k) { + const float r = 60.0f + 480.0f * static_cast(k % 20) / 20.0f; + const float phi = 0.61f * static_cast(k); + const float x = bx + r * std::cos(phi), y = by + r * std::sin(phi); + const float ux = std::cos(phi), uy = std::sin(phi); + const float sr = 1.0f + r / 200.0f, st = 1.0f + r / 100.0f; + for (int py = static_cast(y) - 30; py <= static_cast(y) + 30; ++py) + for (int px = static_cast(x) - 30; px <= static_cast(x) + 30; ++px) { + if (px < 0 || py < 0 || px >= W || py >= H) continue; + const float dx = px - x, dy = py - y; + const float rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; + img[static_cast(py) * W + px] += static_cast( + std::lround(2000.0f * std::exp(-rad * rad / (2 * sr * sr) - tn * tn / (2 * st * st)))); + } + xs.push_back(x); + ys.push_back(y); + } + // Overlapping spots are not what this checks: keep those far from every other one. + std::vector kx, ky; + for (size_t i = 0; i < xs.size(); ++i) { + bool alone = true; + for (size_t j = 0; j < xs.size(); ++j) + if (i != j && std::hypot(xs[i] - xs[j], ys[i] - ys[j]) < 45.0f) alone = false; + if (alone) { kx.push_back(xs[i]); ky.push_back(ys[i]); } + } + REQUIRE(kx.size() > 40); + std::vector spots; + MeasureFootprintSpots(img.data(), W, H, bx, by, kx, ky, spots); + REQUIRE(spots.size() > 30); + for (const auto &s : spots) { + CHECK_THAT(s.sigma_rad, Catch::Matchers::WithinRel(1.0f + s.r_px / 200.0f, 0.12f)); + CHECK_THAT(s.sigma_tan, Catch::Matchers::WithinRel(1.0f + s.r_px / 100.0f, 0.12f)); + } +} + +TEST_CASE("SpotFootprint_TableFillsSparseBinsFromNeighbours", "[Integration][portable]") { + std::vector spots; + for (int i = 0; i < FOOTPRINT_MIN_SPOTS_PER_BIN; ++i) { + spots.push_back({50.0f, 1.0f, 1.5f}); // bin 0 + spots.push_back({1150.0f, 3.0f, 4.0f}); // last bin + } + const SpotFootprint fp = FootprintFromSpots(spots, 1200.0f); + REQUIRE(fp.sigma_rad.size() == static_cast(FOOTPRINT_BINS)); + REQUIRE(fp.sigma_rad.front() == 1.0f); + REQUIRE(fp.sigma_tan.back() == 4.0f); + REQUIRE(fp.sigma_rad[2] == 1.0f); // nearer the first bin + REQUIRE(fp.sigma_rad[FOOTPRINT_BINS - 3] == 3.0f); // nearer the last + REQUIRE(FootprintFromSpots({}, 1200.0f).empty()); +} From a6378800b85d2605f5d0a67dd18f32c6eb92e935 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 4 Oct 2026 02:42:09 +0200 Subject: [PATCH 11/13] Footprint integration: GPU/CPU parity sections, docs, changelog BraggIntegrationEngineGPU_MatchesCPU gains three footprint sections (spaced, crowded under overlap exclude, with the radial background correction): both engines classify the summation ellipse, the grown ring and the footprint Gaussian alike. Integration chapter and changelog describe the measured footprint. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- docs/CHANGELOG.md | 1 + docs/CPU_DATA_ANALYSIS_INTEGRATION.md | 2 ++ tests/BraggIntegrationEngineGPUTest.cpp | 32 +++++++++++++++++++++++-- 3 files changed, 33 insertions(+), 2 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 9a34391b6..8cbe704c3 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -4,6 +4,7 @@ ### 1.0.0-rc.174 * Rugnux reads Rigaku d*TREK SMV images (Saturn CCD), including detector 2theta, image orientation and encoded pixel overflows. +* Rugnux integrates spots that grow wider than the integration disk away from the beam (typical of small molecules at high X-ray energy) over their measured footprint. ### 1.0.0-rc.173 diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index cf9aed468..4a55123ab 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -106,6 +106,8 @@ Only the ring moves. The signal disk $r_1$ stays circular, deliberately: it sets What a circular $r_1$ loses is flux, and that loss is **not** a function of resolution alone: measured per reflection, it carries a directional component worth several Ų with a definite principal axis, on top of the isotropic part. Nor is there anything in the merge to absorb it. There is **no per-shell scale**, and there cannot usefully be one: every scale in §10 is fitted against a reference built from a reflection's own symmetry equivalents, and equivalents share $s^2$ exactly, so any function of $s^2$ lies in the exact null space of the whole scaling model — a per-shell parameter would have zero residual to fit against. (XDS and DIALS have the same null space, for the same reason.) The isotropic part of the loss is instead degenerate with the overall Wilson $B$ and is silently reported as part of it, so **the reported `WILSON_B` / `_reflns.B_iso_Wilson_estimate` carries an $r_1$-dependent contribution**: measured across a constant-ring-area radius sweep it falls monotonically as the disk grows, by 0.5 Ų on sharp strong data and by up to ~10 Ų on weak wide-spot data. What this costs the *data* is much less than what it costs the flux, because most of the loss is matched by a proportional $\sigma$: it moves no CC$_{1/2}$ and no $R_\text{meas}$, and — to within a few hundredths of an ångström — no resolution cut. +**Measured spot footprint (automatic).** The radii above are chosen from spots at 5 Å, which at high X-ray energy sit close to the beam. Away from it a spot can grow several times wider — radially from the sensor's parallax and the obliquity of the incidence, tangentially from the crystal's azimuthal spread, which rotates the diffracted beam about the incident one and smears the spot along its ring. On small-molecule data at 20–25 keV the standard deviation grows from ~1 px near the beam to ~5 px at the detector edge: the $r_1 = 4$ disk holds a quarter of the flux there, the $6\ldots13$ px ring a third of it, and the profile widths learned inside $r_1$ (§9.3) saturate near $r_1^2/4$. So the pre-scan measures every spot it finds with a window that follows the spot — three of its own standard deviations, iterated and re-centred — separately along and across the radius, and tabulates the median widths $\sigma_ ho,\sigma_ au$ against the distance from the beam. Wherever $3\max(\sigma_ ho,\sigma_ au)>r_1$ the integrator then (i) starts the background ring at $3\sigma$ along each axis, (ii) sums the reflection over the $r_1$ disk **and** the $3\sigma$ footprint ellipse, so the summation — the profile fit's seed and its fallback — holds the spot rather than its core, and (iii) builds the per-reflection Gaussian at the measured widths on a grid grown to hold them. Where every spot fits the disk nothing is installed and the integration is unchanged bit for bit, which is the case for compact protein spots; like the measured radius, the footprint applies to the canonical pass and not to the geometry pre-pass, and a canonical pass whose wider rings the neighbours starve falls back to the settings without it. Judged by refining the published structures with SHELXL, it removes the intensity loss that grew with resolution on the small-molecule sets (rugnux/model intensity in the outermost shell 0.81–0.91 → 0.98–1.02). + ### 9.2 Box summation (seed and fallback) Let: diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index 2f5c45425..55d160c6e 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -146,10 +146,16 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, float stencil_k = 0.0f, float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f, OverlapMode overlap = OverlapMode::Off, float companion_dx = 0.0f, - bool clip_spots = false, float odd_partiality = 1.0f) { - const DiffractionExperiment experiment = + bool clip_spots = false, float odd_partiality = 1.0f, + const SpotFootprint *footprint = nullptr) { + DiffractionExperiment experiment = MakeExperiment(mode, bandwidth_fwhm, clip_nsigma, radial, DetJF(2), stencil_k, r1, r2, r3, overlap); + if (footprint) { + BraggIntegrationSettings settings = experiment.GetBraggIntegrationSettings(); + settings.Footprint(*footprint); + experiment.ImportBraggIntegrationSettings(settings); + } const size_t width = experiment.GetXPixelsNum(); const size_t height = experiment.GetYPixelsNum(); const size_t npixel = experiment.GetPixelsNum(); @@ -254,6 +260,28 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, 0.04f, 0.0f, false, 120, 4.0f, 6.0f, 8.0f, 12.0f); } SECTION("ProfileEmpirical") { CompareCpuVsGpu(IntegratorMode::ProfileEmpirical, std::nullopt); } + // A measured footprint wider than the r1 disk (SpotFootprint.h): the summation ellipse, the ring + // grown along and across the radius and the footprint-width Gaussian are all reflection-dependent + // geometry the two engines have to classify alike - spaced so the rings stay clear, and crowded + // so the grown neighbour mask and the ring gate come into play. + SpotFootprint fp; + fp.bin_px = 100.0f; + for (int b = 0; b < 8; ++b) { + fp.sigma_rad.push_back(1.5f + 0.2f * b); + fp.sigma_tan.push_back(1.8f + 0.3f * b); + } + SECTION("ProfileGaussian footprint") { + CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 120, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 1.0f, &fp); + } + SECTION("ProfileGaussian footprint crowded exclude") { + CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 40, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Exclude, 0.0f, false, 1.0f, &fp); + } + SECTION("ProfileGaussian footprint radial") { + CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, true, 120, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 1.0f, &fp); + } // Overlap treatment: companions 4 px apart put each reflection's centre inside its neighbour's // signal disk, so the owner map, the excluded pixels and the profile fraction the two modes act on // all have to come out the same in both engines - the ownership atomic in particular is settled by From 9c0a8c740cb543a4778547e9bef9441e45503dfe Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 4 Oct 2026 11:09:25 +0200 Subject: [PATCH 12/13] Battery: drop 7UDI from the open arm Its deposited sweep lacks 40 deg (frames 561-720 are not in the archive; the depositors processed it as two sweeps), so it tests gap handling rather than processing. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- docs/EXTERNAL_TEST_DATA.md | 1 - tools/battery/open.json | 1 - 2 files changed, 2 deletions(-) diff --git a/docs/EXTERNAL_TEST_DATA.md b/docs/EXTERNAL_TEST_DATA.md index 7db2cd403..210890797 100644 --- a/docs/EXTERNAL_TEST_DATA.md +++ b/docs/EXTERNAL_TEST_DATA.md @@ -131,7 +131,6 @@ the table below; the repositories themselves are cited in | [7RJI](https://www.rcsb.org/structure/7RJI) | IRRMC [10.18430/M37RJI](https://doi.org/10.18430/M37RJI) | LNLS W01B-MX2 | 1.71 | H 3 2 | 83.0 83.0 124.8 90.0 90.0 120.0 | PILATUS 2M | BthTX-II variant b, from Bothrops jararacussu venom, complexed with stearic acid | | [7T5T](https://www.rcsb.org/structure/7T5T) | SBGrid [10.15785/sbgrid/864](https://doi.org/10.15785/sbgrid/864) | SSRL BL9-2 | 1.35 | P 42 21 2 | 95.3 95.3 104.9 90.0 90.0 90.0 | PILATUS 6M | Structure of Thauera sp. K11 CapP | | [7TCD](https://www.rcsb.org/structure/7TCD) | IRRMC [10.18430/m37tcd](https://doi.org/10.18430/m37tcd) | SLS X06SA | 1.70 | C 1 2 1 | 138.5 47.9 78.1 90.0 107.6 90.0 | Dectris Eiger 16M | LOV2-DARPIN fusion: D13 | -| [7UDI](https://www.rcsb.org/structure/7UDI) | Zenodo [10.5281/zenodo.10022358](https://doi.org/10.5281/zenodo.10022358) | CLSI 08B1-1 | 2.24 | P 41 | 66.7 66.7 129.6 90.0 90.0 90.0 | PILATUS3 S 6M | Full-length dimer of DNA-Damage Response Protein C from Deinococcus radiodurans | | [7YZX](https://www.rcsb.org/structure/7YZX) | IRRMC [10.18430/M37YZX](https://doi.org/10.18430/M37YZX) | Diamond I24 | 1.90 | P 63 2 2 | 169.4 169.4 141.8 90.0 90.0 120.0 | PILATUS3 6M | ScpA from Streptococcus pyogenes, D783A mutant. | | [8A1A](https://www.rcsb.org/structure/8A1A) | IRRMC [10.18430/M38A1A](https://doi.org/10.18430/M38A1A) | SLS X06SA | 2.05 | P 65 | 191.9 191.9 122.4 90.0 90.0 120.0 | Dectris Eiger 16M | Structure of a leucinostatin derivative determined by host lattice display : L1F11V1 construct | | [8AGQ](https://www.rcsb.org/structure/8AGQ) | IRRMC [10.18430/M38AGQ](https://doi.org/10.18430/M38AGQ) | SLS X06DA | 1.09 | C 1 2 1 | 89.9 55.4 54.8 90.0 113.5 90.0 | PILATUS 2MF | Crystal structure of anthocyanin-related GSTF8 from Populus trichocarpa in complex with (-)-catechin and glutathione | diff --git a/tools/battery/open.json b/tools/battery/open.json index 5da289926..1f10a48f8 100644 --- a/tools/battery/open.json +++ b/tools/battery/open.json @@ -198,7 +198,6 @@ {"id": "8v2t", "input": "8v2t/data/MJ504_2_3_00001.cbf", "wavelength": 1.1, "ref": {"sg": "P 42 21 2", "sgno": 94, "cell": [60.906, 60.906, 92.728, 90.0, 90.0, 90.0], "dmin": 1.402}, "tags": ["cbf", "tetragonal"]}, {"id": "5jk4", "input": "5jk4/dioxyhipeg/dioxyhipeg_1_001.img", "wavelength": 0.933, "ref": {"sg": "P 1 21 1", "sgno": 4, "cell": [37.729, 77.856, 56.287, 90.0, 102.08, 90.0], "dmin": 1.1}, "tags": ["smv", "monoclinic"], "pinned": true}, {"id": "6cs9", "input": "6cs9/568/HBD2D7_2_1_001.img", "wavelength": 0.9537, "ref": {"sg": "P 1 21 1", "sgno": 4, "cell": [32.871, 25.538, 40.17, 90.0, 98.64, 90.0], "dmin": 1.85}, "tags": ["smv", "monoclinic"]}, - {"id": "7udi", "input": "7udi/data/MJ7121_0001.cbf", "wavelength": 0.9795, "ref": {"sg": "P 41", "sgno": 76, "cell": [66.698, 66.698, 129.581, 90.0, 90.0, 90.0], "dmin": 2.24}, "tags": ["cbf", "tetragonal"]}, {"id": "9lxl", "input": "9lxl/6/1_T10Y-FFp-PNPI_6_master.h5", "wavelength": 0.97919, "ref": {"sg": "P 41 21 2", "sgno": 92, "cell": [76.757, 76.757, 225.338, 90.0, 90.0, 90.0], "dmin": 2.191}, "tags": ["h5", "tetragonal"]}, {"id": "9s02", "input": "9s02/PYCR1-D11_1_master.h5", "wavelength": 0.97625, "ref": {"sg": "P 21 21 2", "sgno": 18, "cell": [163.674, 88.042, 116.731, 90.0, 90.0, 90.0], "dmin": 1.65}, "tags": ["h5", "orthorhombic"]}, {"id": "6gvk", "input": "6gvk/212_3_0001.cbf", "wavelength": 0.97915, "ref": {"sg": "C 1 2 1", "sgno": 5, "cell": [105.58, 59.52, 42.4, 90.0, 113.5, 90.0], "dmin": 1.55}, "tags": ["cbf", "monoclinic"], "pinned": true} From 15fb2543ca99e1efdba2f2a70ba6badffdee82c5 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 4 Oct 2026 11:26:58 +0200 Subject: [PATCH 13/13] Revert "RotationScaleMerge: refit the error model about the merge's own mean" This reverts commit 1baf92606. The refit fixed one small-molecule set but cost a split crystal its screw axis and moved ISa down across the battery; the rc173 error model stays. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- .../scale_merge/RotationScaleMerge.cpp | 111 +----------------- .../scale_merge/RotationScaleMergeGPU.cu | 8 -- .../scale_merge/RotationScaleMergeGPU.h | 3 - 3 files changed, 4 insertions(+), 118 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 7885d767c..b1371ee4e 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -4478,104 +4478,6 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool resid * resid / factor, o.d}); } } - fit_error_model(samples); - - // ---- The error model refitted about the mean the merge takes. ---- - // The fit above is made about a mean that weighs each full by its COUNTING variance. Where the - // equivalents of a reflection agree to within counting statistics that is the merge's own centre, - // and the model it gives is the model the merge applies. Where they do not - a systematic error - // the counting variance does not know about - it is the wrong centre, and in the worst way: the - // full with the fewest counts has the smallest variance, so the mean, and every deviation the - // model is fitted on, is pulled toward the observations that read lowest. Measured on a strongly - // absorbing crystal whose equivalents spread over a factor 100: the mean of (0 4 0) sat at 21k - // among observations from 11k to 1.9M, the deviations of the rest read as hundreds of counting - // sigmas, a ran into its bound and the six-sigma test about the median taken with the same - // weights deleted 63% of the observations - the strong ones (SHELXL R1 0.50 on that merge). - // So refit about the mean the MERGE takes - each full weighted by the variance the fitted model - // gives it, evaluated at the reflection's mean - until the model stops moving; the outlier - // test's median below takes the same weights. Where counting statistics were right the first fit - // is the answer and nothing moves; where b dominates, the weights even out, as the model says. - // Following Blessing (1997) J. Appl. Cryst. 30, 421-426: the centre an outlier is judged from - // is weighted as the merge weighs. - if (error_model_active) { - // Fulls in group order, so each pass over a group's members is one contiguous walk and the - // groups split across threads without every thread reading every full. - std::vector gstart(n_groups + 1, 0); - for (int i = 0; i < n_full; ++i) - if (mf.group[i] >= 0 && cnt[mf.group[i]] >= 2) ++gstart[mf.group[i] + 1]; - for (int g = 0; g < n_groups; ++g) gstart[g + 1] += gstart[g]; - std::vector member(gstart[n_groups]); - { - std::vector fill(gstart.begin(), gstart.end() - 1); - for (int i = 0; i < n_full; ++i) - if (mf.group[i] >= 0 && cnt[mf.group[i]] >= 2) member[fill[mf.group[i]]++] = i; - } - std::vector next_mean; - std::vector next; - constexpr int REFIT_ROUNDS = 10; - for (int round = 0; round < REFIT_ROUNDS; ++round) { - const double a = error_model_a, b = error_model_b; - const auto model_var = [&](int i, double mean) { - const double sc = static_cast(mf.sigma[i]) * mf.corr[i]; - return a * counting_variance(fulls[i], mean, sc * sc) + (b * mean) * (b * mean); - }; - std::vector> part(ThreadsForWork(member.size(), nthreads)); - const int nt = static_cast(part.size()); - next_mean.assign(n_groups, NAN); - // Each group's mean and samples come from its own members in index order, so the result - // does not depend on how the groups were split. - ParallelFor(nt, nt, [&](int t) { - const int g0 = static_cast(static_cast(n_groups) * t / nt); - const int g1 = static_cast(static_cast(n_groups) * (t + 1) / nt); - for (int g = g0; g < g1; ++g) { - if (gstart[g + 1] - gstart[g] < 2 || !std::isfinite(em_mean[g])) continue; - double sw = 0.0, swI = 0.0, swh[2] = {0.0, 0.0}, swIh[2] = {0.0, 0.0}; - int nh[2] = {0, 0}; - for (int q = gstart[g]; q < gstart[g + 1]; ++q) { - const int i = member[q]; - const double v = model_var(i, em_mean[g]); - if (!(v > 0.0)) continue; - const double I_corr = static_cast(mf.I[i]) * mf.corr[i]; - sw += 1.0 / v; swI += I_corr / v; - if (merge_friedel && group_has_hands[g]) { - swh[obs_hand[i]] += 1.0 / v; swIh[obs_hand[i]] += I_corr / v; nh[obs_hand[i]]++; - } - } - if (!(sw > 0.0)) continue; - const double mean = swI / sw; - next_mean[g] = mean; - // The samples about the hand's own mean where it has two of its own, as above. - for (int q = gstart[g]; q < gstart[g + 1]; ++q) { - const int i = member[q]; - const double v = model_var(i, em_mean[g]); - if (!(v > 0.0)) continue; - const int hh = obs_hand[i]; - const bool on_hand = merge_friedel && group_has_hands[g] && nh[hh] >= 2 && swh[hh] > 0.0; - const double centre = on_hand ? swIh[hh] / swh[hh] : mean; - const double factor = 1.0 - (1.0 / v) / (on_hand ? swh[hh] : sw); - if (factor < 0.05) continue; - const double sc = static_cast(mf.sigma[i]) * mf.corr[i]; - const double resid = static_cast(mf.I[i]) * mf.corr[i] - centre; - part[t].push_back({counting_variance(fulls[i], mean, sc * sc), centre * centre, - resid * resid / factor, mf.d[i]}); - } - } - }); - for (int g = 0; g < n_groups; ++g) - if (!std::isfinite(next_mean[g])) next_mean[g] = em_mean[g]; - next.clear(); - for (auto &v : part) next.insert(next.end(), v.begin(), v.end()); - em_mean.swap(next_mean); - samples.swap(next); - fit_error_model(samples); - if (std::fabs(error_model_a - a) <= 1e-3 * a - && std::fabs(error_model_b - b) <= 1e-3 * std::max(b, 1e-6)) - break; - } -#ifdef JFJOCH_USE_CUDA - if (use_gpu_merge) gpu_->SetEmMean(em_mean.data()); -#endif - } // Per-group outlier-rejection median of I*corr (host both paths - a per-group median is awkward on // the GPU; cheap here, cnt >= 3 filter from the em pass). Fed to the merge accumulate. @@ -4638,21 +4540,15 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // Inverse-variance WEIGHTED median: a frame the crystal barely diffracted on is scaled up by // 1/G together with its sigma, so a plain median lets two such observations outvote one // well-measured one - and the test below then rejects the well-measured one against its - // own small sigma. The variance is the model's, at the reflection's mean - the weight the - // merge gives the observation (see the centre refit above). - std::vector> iv(start[n_sets]); // (I*corr, 1/model variance) + // own small sigma. + std::vector> iv(start[n_sets]); // (I*corr, 1/(sigma*corr)^2) { std::vector fill(start.begin(), start.end() - 1); for (int i = 0; i < n_full; ++i) { const int g = mf.group[i]; if (g < 0) continue; const float sc = mf.sigma[i] * mf.corr[i]; - const double mean = std::isfinite(em_mean[g]) ? em_mean[g] : static_cast(mf.I[i]) * mf.corr[i]; - const double mv = error_model_active - ? error_model_a * counting_variance(fulls[i], mean, static_cast(sc) * sc) - + (error_model_b * mean) * (error_model_b * mean) - : static_cast(sc) * sc; - const std::pair v{mf.I[i] * mf.corr[i], mv > 0.0 ? static_cast(1.0 / mv) : 0.0f}; + const std::pair v{mf.I[i] * mf.corr[i], sc > 0.0f ? 1.0f / (sc * sc) : 0.0f}; if (cnt[g] >= 3) iv[fill[g]++] = v; if (!pair_needed.empty() && pair_needed[pair_of_group[g]]) iv[fill[n_groups + pair_of_group[g]]++] = v; @@ -4679,6 +4575,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool reject_median[g] = cnt[g] >= 3 ? set_median[g] : !pair_needed.empty() ? set_median[n_groups + pair_of_group[g]] : NAN; } + fit_error_model(samples); } // The full's sigma under the error model, with the variance evaluated at intensity I_for_b. diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 257886bc3..0ecf08707 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -1126,14 +1126,6 @@ void RotationScaleMergeGPU::SetFrameCellOk(const uint8_t *frame_cell_ok) { // The per-group inv-var mean (em_mean) + the per-full leverage-corrected error-model samples over the // resident+scaled fulls. Stashes the filter context for the later MergeAccum/MergeRmeas calls. -void RotationScaleMergeGPU::SetEmMean(const double *em_mean) { - DeviceGuard guard(impl_->device, impl_->available); - auto &d = *impl_; - if (d.n_groups > 0) - CopyAndWait(d.m_em_mean.get(), em_mean, size_t(d.n_groups) * sizeof(double), cudaMemcpyHostToDevice, - impl_->s(), "ul em_mean"); -} - void RotationScaleMergeGPU::MergeEmSamples(bool for_search, double min_partiality, const uint8_t *hand, const uint8_t *has_hands, double *em_mean_out, int32_t *cnt_out, double *s2_out, diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index 349343b1a..98aa1b18a 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -103,9 +103,6 @@ public: // half-set weights multiplied by it. Requires MergeEmSamples first (em_mean resident). // reject_var_add (n_groups) widens the pooled cut by the shell's own measured Bijvoet // variance; null leaves the plain n-sigma test. - // Replace the per-group means MergeEmSamples left on the device (n_groups values): the merge's - // model sigmas are evaluated at them. - void SetEmMean(const double *em_mean); void MergeAccum(double error_model_a, double error_model_b, bool error_model_active, bool reject_outliers, double reject_nsigma, const float *reject_median, const float *reject_var_add,