From ffc990a330c7b3af0fe99a4ab0e14fb26226313e Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 3 Nov 2026 06:13:45 +0100 Subject: [PATCH] rugnux: panel study log of held-out residual per panel and shell (RUGNUX_PANEL_STUDY) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UqXYdJKPUFLoeAchjzv3x4 --- rugnux/RugnuxStills.cpp | 90 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 79 insertions(+), 11 deletions(-) diff --git a/rugnux/RugnuxStills.cpp b/rugnux/RugnuxStills.cpp index df4dfa57d..38945d030 100644 --- a/rugnux/RugnuxStills.cpp +++ b/rugnux/RugnuxStills.cpp @@ -116,6 +116,59 @@ GeomRefineFrame ToRefineFrame(const DataMessage &msg, int ordinal) { return f; } +// One held-out fit: the geometry fitted on one half, and its residuals on the other half's spots. +struct HeldOutFit { + double rms = std::numeric_limits::infinity(); + std::vector modules; // the modules as fitted on the fitting half + GeometryRefinerResult eval; // the frozen geometry on the other half + DiffractionGeometry geom; +}; + +// PANEL STUDY (RUGNUX_PANEL_STUDY): the held-out residual of each fold per module and per resolution shell, +// and each module's shift on the two halves, to the log under "PanelStudy". +void LogPanelStudy(Logger &logger, const std::string &label, const std::vector half[2], + const HeldOutFit fits[2]) { + constexpr double CAP_PX = 2.0; + const std::vector edges = {1e9, 4.0, 3.0, 2.5, 2.2, 2.0, 1.8, 1.6, 1.4, 0.0}; + const size_t nm = fits[0].modules.size(); + std::vector m_sq(nm, 0), s_sq(edges.size(), 0); + std::vector m_n(nm, 0), m_hit(nm, 0), s_n(edges.size(), 0), s_hit(edges.size(), 0); + for (int fit = 0; fit < 2; fit++) { + const auto &e = fits[fit].eval; + const auto &frames = half[1 - fit]; + if (e.predicted_px.size() != frames.size()) + continue; + for (size_t i = 0; i < frames.size(); i++) + for (size_t j = 0; j < frames[i].spots.size() && j < e.predicted_px[i].size(); j++) { + const auto &sp = frames[i].spots[j]; + const auto &pr = e.predicted_px[i][j]; + double r = CAP_PX; + if (std::isfinite(pr[0]) && std::isfinite(pr[1])) + r = std::min(CAP_PX, static_cast(std::hypot(pr[0] - sp.x, pr[1] - sp.y))); + const double d = 1.0 / std::max(1e-6f, fits[fit].geom.DetectorToRecip(sp.x, sp.y).Length()); + size_t sh = 1; + while (sh + 1 < edges.size() && d < edges[sh]) + sh++; + s_sq[sh] += r * r; s_n[sh]++; s_hit[sh] += r < 1.0; + if (sp.module >= 0 && static_cast(sp.module) < nm) { + m_sq[sp.module] += r * r; m_n[sp.module]++; m_hit[sp.module] += r < 1.0; + } + } + } + for (size_t sh = 1; sh < edges.size(); sh++) + if (s_n[sh] > 0) + logger.Info("PanelStudy {} shell {:.2f} {:.2f} n {} rms {:.4f} hit1 {:.4f}", label, edges[sh - 1], edges[sh], + s_n[sh], std::sqrt(s_sq[sh] / s_n[sh]), s_hit[sh] / static_cast(s_n[sh])); + for (size_t m = 0; m < nm; m++) + if (m_n[m] > 0) + logger.Info("PanelStudy {} module {} x0 {:.0f} y0 {:.0f} n {} rms {:.4f} hit1 {:.4f} " + "fold0 {:+.3f} {:+.3f} {} fold1 {:+.3f} {:+.3f} {}", label, m, + fits[0].modules[m].x0, fits[0].modules[m].y0, m_n[m], std::sqrt(m_sq[m] / m_n[m]), + m_hit[m] / static_cast(m_n[m]), + fits[0].modules[m].dx, fits[0].modules[m].dy, fits[0].modules[m].spots, + fits[1].modules[m].dx, fits[1].modules[m].dy, fits[1].modules[m].spots); +} + bool SameGeometry(const DiffractionGeometry &a, const DiffractionGeometry &b) { return a.GetBeamX_pxl() == b.GetBeamX_pxl() && a.GetBeamY_pxl() == b.GetBeamY_pxl() && a.GetDetectorDistance_mm() == b.GetDetectorDistance_mm() && a.GetPixelSize_mm() == b.GetPixelSize_mm() @@ -506,32 +559,47 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ }); // Held-out residual of the geometry fitted on half `fit`, scored on every spot of the other half. - auto held_out = [&](const std::vector &modules, bool refine, int fit) { + auto held_out = [&](const std::vector &modules, bool refine, int fit, int min_spots) { + HeldOutFit h; GeometryRefinerSettings sf = gr_settings; sf.num_threads = std::max(1, config_.nthreads / 8); sf.modules = modules; sf.refine_modules = refine; + sf.min_module_spots = min_spots; const auto f = RefineGlobalGeometry(nominal, *cell, half[fit], sf); if (!f.ok) - return std::numeric_limits::infinity(); + return h; DiffractionGeometry g = nominal; g.BeamX_pxl(f.beam_x_px).BeamY_pxl(f.beam_y_px).DetectorDistance_mm(f.distance_mm); GeometryRefinerSettings ev = gr_settings; ev.num_threads = sf.num_threads; ev.refine_beam_distance = false; ev.modules = f.modules; - const auto e = RefineGlobalGeometry(g, f.cell, half[1 - fit], ev); - if (!e.ok) - return std::numeric_limits::infinity(); - return e.capped_rms_px; + h.modules = f.modules; + h.geom = g; + h.eval = RefineGlobalGeometry(g, f.cell, half[1 - fit], ev); + if (h.eval.ok) + h.rms = h.eval.capped_rms_px; + return h; }; - std::future held_fits[2], moved_fits[2]; + const int min_spots = gr_settings.min_module_spots; + std::future held_fits[2], moved_fits[2]; for (int fit = 0; fit < 2; fit++) { - held_fits[fit] = std::async(std::launch::async, held_out, std::cref(centres), false, fit); - moved_fits[fit] = std::async(std::launch::async, held_out, std::cref(start), true, fit); + held_fits[fit] = std::async(std::launch::async, held_out, std::cref(centres), false, fit, min_spots); + moved_fits[fit] = std::async(std::launch::async, held_out, std::cref(start), true, fit, min_spots); + } + HeldOutFit held_f[2] = {held_fits[0].get(), held_fits[1].get()}; + HeldOutFit moved_f[2] = {moved_fits[0].get(), moved_fits[1].get()}; + const double held = (held_f[0].rms + held_f[1].rms) / 2.0; + const double moved = (moved_f[0].rms + moved_f[1].rms) / 2.0; + if (std::getenv("RUGNUX_PANEL_STUDY")) { + LogPanelStudy(logger, "given", half, held_f); + LogPanelStudy(logger, "placed", half, moved_f); + HeldOutFit low_f[2] = {held_out(start, true, 0, 100), held_out(start, true, 1, 100)}; + LogPanelStudy(logger, "placed100", half, low_f); + logger.Info("PanelStudy summary given {:.4f} placed {:.4f} placed100 {:.4f}", held, moved, + (low_f[0].rms + low_f[1].rms) / 2.0); } - const double held = (held_fits[0].get() + held_fits[1].get()) / 2.0; - const double moved = (moved_fits[0].get() + moved_fits[1].get()) / 2.0; const GeometryRefinerResult rm = placement.get(); placement_mark.Stop(); logger.Info("Module placement: {} modules; held-out spot residual (capped rms) {:.3f} px as given, "