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

stills fixed cost: geometry refinement and pre-scan do the same work in less wall time

Browse Source 
Four changes to the per-run fixed cost of a stills run (Step 2 geometry refinement, Step 3 pre-scan),
none of which changes what is computed:

- Module placement: the four held-out fits and the placement on every strong frame are independent, and
  one Ceres fit keeps only ~4.6-6.5 of 32 cores busy (perf, JF16M, 4000 frames / 150888 spots). They
  now run side by side on a fifth of the threads each; the placement fit is computed speculatively and
  used only where the held-out residuals call for it, as before. Baseline (quiet): held-out 12.3 s +
  placement 5.1 s on JF16M, 16.4 + 6.5 s on JF4M, 19.0 + 7.6 s on MPCCD.
- Module scan (+-25 px grid per module ahead of the placement) was one thread: every grid shift is now
  scored on the workers and the grid then walked in its old order, so the first best still wins
  (3.3 s on JF16M, ~5 s on JF4M/MPCCD before). New test ScanModuleShifts_FindsShiftedModule checks the
  shift is found and that 1 and 4 threads give the same answer.
- Stills geometry trial: each hypothesis (centre, distance, frame sample) is indexed once. The distance
  grid ran twice on the same centre whenever the file's centre was kept (40 of ~275 hypotheses,
  ~15% of the trial: 6.2 s on JF16M, 6.7 s on the flat JF8M), and the final here/there check repeated
  the first trial. Repeats always gave the same count (checked on three logs).
- The saturation histogram (40 frames, 0.1-2.2 s) reads only frames and the mask: it is taken beside
  the geometry refinement and the limit decided where it was, after it.

Plus TIMING marks for the sub-steps (geometry refinement sample, module scan, module placement fits).

Outputs identical to the baseline binary (bb7c9c) where checked: JF16M int16 15 runs (geometry, all
32 module shifts, held-out residuals, cell), JF4M int16, MPCCD 236 (indexing rate, CC1/2, R_split,
R-free with --model identical), PAL 90 and LCLS 222 (identical reports).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-10-21 17:48:09 +02:00
1 parent 987e812e23
commit 96322a8ec9
6 files changed
+171 -64
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

tests/ModuleRefinementTest.cpp
+76 -19
View File
@@ -9,27 +9,14 @@
#include "../image_analysis/geom_refinement/GeometryRefiner.h"
#include "../image_analysis/bragg_prediction/BraggPrediction.h"
// A detector of four modules (the quadrants around the beam), each displaced by a known rigid move: the
// bundle adjustment with module placement must find the moves, and the fit without it must be worse.
TEST_CASE("RefineGlobalGeometry_ModulePlacement") {
DiffractionExperiment exp;
exp.IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(150);
const DiffractionGeometry geom = exp.GetDiffractionGeometry();
namespace {
// Moves that hold none of what the gauge takes away: zero in sum, in turn and outwards.
std::vector<ModuleShift> truth(4);
const double centre[4][2] = {{500, 500}, {1500, 500}, {500, 1500}, {1500, 1500}};
const double move[4][3] = {{2.0, -2.0, 0.2}, {2.0, 2.0, -0.2}, {-2.0, -2.0, 0.1}, {-2.0, 2.0, -0.1}};
for (int m = 0; m < 4; m++) {
truth[m].x0 = centre[m][0];
truth[m].y0 = centre[m][1];
truth[m].dx = move[m][0];
truth[m].dy = move[m][1];
truth[m].rot_rad = move[m][2] * PI / 180.0;
}
// Frames of a 50/60/70 A orthorhombic crystal in random orientations on a detector of four modules (the
// quadrants around the beam at 1000, 1000), each displaced by the rigid move `truth` gives it: the spots are
// where the image shows them, each indexed as predicted.
std::vector<GeomRefineFrame> FramesOnMovedModules(const DiffractionExperiment &exp, const CrystalLattice &latt,
const std::vector<ModuleShift> &truth) {
auto module_of = [](double x, double y) { return (x < 1000 ? 0 : 1) + (y < 1000 ? 0 : 2); };
const CrystalLattice latt(50, 60, 70, 90, 90, 90);
BraggPredictionSettings settings{.high_res_A = 2.0, .ewald_dist_cutoff = 0.001};
std::mt19937 rng(7);
std::uniform_real_distribution<float> u(-1.0f, 1.0f);
@@ -62,6 +49,31 @@ TEST_CASE("RefineGlobalGeometry_ModulePlacement") {
if (frame.spots.size() >= 20)
frames.push_back(frame);
}
return frames;
}
} // namespace
// A detector of four modules (the quadrants around the beam), each displaced by a known rigid move: the
// bundle adjustment with module placement must find the moves, and the fit without it must be worse.
TEST_CASE("RefineGlobalGeometry_ModulePlacement") {
DiffractionExperiment exp;
exp.IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(150);
const DiffractionGeometry geom = exp.GetDiffractionGeometry();
// Moves that hold none of what the gauge takes away: zero in sum, in turn and outwards.
std::vector<ModuleShift> truth(4);
const double centre[4][2] = {{500, 500}, {1500, 500}, {500, 1500}, {1500, 1500}};
const double move[4][3] = {{2.0, -2.0, 0.2}, {2.0, 2.0, -0.2}, {-2.0, -2.0, 0.1}, {-2.0, 2.0, -0.1}};
for (int m = 0; m < 4; m++) {
truth[m].x0 = centre[m][0];
truth[m].y0 = centre[m][1];
truth[m].dx = move[m][0];
truth[m].dy = move[m][1];
truth[m].rot_rad = move[m][2] * PI / 180.0;
}
const CrystalLattice latt(50, 60, 70, 90, 90, 90);
const std::vector<GeomRefineFrame> frames = FramesOnMovedModules(exp, latt, truth);
REQUIRE(frames.size() > 30);
GeometryRefinerSettings st;
@@ -92,3 +104,48 @@ TEST_CASE("RefineGlobalGeometry_ModulePlacement") {
CHECK(placed.beam_y_px == Catch::Approx(1000).margin(0.1));
CHECK(placed.distance_mm == Catch::Approx(150).margin(0.05));
}
// The coarse search ahead of the placement: a module off by several pixels is found near where it is, and the
// search gives the same answer on any number of threads.
TEST_CASE("ScanModuleShifts_FindsShiftedModule") {
DiffractionExperiment exp;
exp.IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(150);
const DiffractionGeometry geom = exp.GetDiffractionGeometry();
std::vector<ModuleShift> truth(4), given(4);
const double centre[4][2] = {{500, 500}, {1500, 500}, {500, 1500}, {1500, 1500}};
for (int m = 0; m < 4; m++) {
truth[m].x0 = given[m].x0 = centre[m][0];
truth[m].y0 = given[m].y0 = centre[m][1];
}
truth[1].dx = 6.0;
truth[1].dy = -4.0;
const CrystalLattice latt(50, 60, 70, 90, 90, 90);
std::vector<GeomRefineFrame> frames = FramesOnMovedModules(exp, latt, truth);
REQUIRE(frames.size() > 30);
// The module's spots did not index where it was given: the search sees them as unindexed.
for (auto &f: frames)
for (auto &s: f.spots)
if (s.module == 1) {
s.h = s.k = s.l = 0;
f.unindexed.push_back(s);
}
for (auto &f: frames)
std::erase_if(f.spots, [](const GeomRefineSpot &s) { return s.module == 1; });
const auto one = ScanModuleShifts(geom, frames, given, 25.0, 1);
const auto many = ScanModuleShifts(geom, frames, given, 25.0, 4);
REQUIRE(one.size() == 4);
// Coarse: a spot is on its lattice within 0.15 of an index, here ~4 px of the spacing, so the search lands
// within that of the true shift and the placement takes it from there.
CHECK(one[1].dx == Catch::Approx(6.0).margin(4.0));
CHECK(one[1].dy == Catch::Approx(-4.0).margin(4.0));
CHECK(one[1].on_lattice_at_best > 2 * one[1].on_lattice_at_zero);
for (int m = 0; m < 4; m++) {
CHECK(many[m].dx == one[m].dx);
CHECK(many[m].dy == one[m].dy);
CHECK(many[m].on_lattice_at_best == one[m].on_lattice_at_best);
CHECK(many[m].on_lattice_at_zero == one[m].on_lattice_at_zero);
}
}