diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index 6ef913795..b34ce7469 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -43,6 +43,21 @@ constexpr double BEAM_BOUND_PXL = 15.0; // +-2.86 deg box the solver works in, so the box reaches everywhere the gate accepts. constexpr double ANGLE_BOUND_DEG = 1.0; +// How far ONE step of the joint fit may move the detector distance or a cell length. This was a cap +// on the whole fit, and as a cap it was the opposite of its job: the header is most worth correcting +// when it is most wrong, and it threw away fits that were right (measured on a crystal whose refused +// 1.39 % move landed on the deposited cell to 0.06 %). What it is still for - a fit corrupted by +// something every cross-validation fold shares, a second lattice most often, which the relative "the +// held-out residual improved" gate cannot see - is kept by making it a TRUST REGION on the step: a +// solve reaches at most this far, and a fit that wants more takes another step from there. A move of +// more than one step is then ratified by re-indexing at it, which is what actually tells a second +// lattice from a wrong header (see PostRefineResult::large_move). +constexpr double STEP_BOUND = 0.01; + +// ... and the walk is bounded, so a fit that never settles still terminates. Five steps is the +-5 % +// box the solver was always given on the distance, i.e. as far as this fit has ever been trusted. +constexpr int MAX_STEPS = 5; + // One integrated partial, flattened across all images. Kept as narrow as the sort and the event // split allow: on a large cell this array is gigabytes, and the scatter and every level of the // per-bucket sort move all of it. The goniometer angle is not stored - it is a function of the @@ -491,14 +506,14 @@ PostRefineResult PostRefineRotationGeometry(const std::vector(np + ne) : 0.0; }; - // One problem, both residual families, every block seeded at nominal. The detector tilt is - // declared and held constant so that a non-zero rot1/rot2 still acts on the observed side. - auto solve_joint = [&](Subset s, double bm[2], double ds[1], double rv[3], + // One problem, both residual families, every block seeded at the geometry it is handed - + // nominal for the first step of the walk below, the last step's answer after that. The + // detector tilt is declared and held constant so that a non-zero rot1/rot2 still acts on + // the observed side. + auto solve_joint = [&](Subset s, bool wide, double bm[2], double ds[1], double rv[3], double q0[3], double q1[3], double q2[3]) { - bm[0] = beam_x0; bm[1] = beam_y0; ds[0] = dist0; - for (int j = 0; j < 3; ++j) { - rv[j] = ax0[j]; q0[j] = p0_0[j]; q1[j] = p1_0[j]; q2[j] = p2_0[j]; - } + const double ds_seed = ds[0]; + const double q1_seed[3] = {q1[0], q1[1], q1[2]}; ceres::Problem p; for (size_t oi = 0; oi < n_obs; ++oi) { const Partial *pp = obs[oi]; @@ -519,20 +534,30 @@ PostRefineResult PostRefineRotationGeometry(const std::vector= static_cast(settings.min_events)) { double pos_nom = 0.0, exc_nom = 0.0, pos_ref = 0.0, exc_ref = 0.0; const double cv_nom = joint_cost(VAL, beam, dist, axv, p0, p1, p2, pos_nom, exc_nom); - double bm_f[2], ds_f[1], rv_f[3], q0_f[3], q1_f[3], q2_f[3]; - const bool convJ = solve_joint(FIT, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f); - const double cv_ref = joint_cost(VAL, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f, pos_ref, exc_ref); - // Commit only for a small, credible move: distance < 1 % (a calibrated header needs - // < ~0.6 %), each cell length < 1 %, each free cell angle < 1 deg (ANGLE_BOUND_DEG), and - // the beam inside the bound measured from whichever centre anything already believes is - // nearer. A larger move is the red flag for an unreliable fit - typically a second - // lattice whose spots bias every cross-validation fold identically, so the relative "it - // improved" gate is blind to it. The absolute size of the move discriminates a genuine - // header correction from that failure far better than the absolute residual, which real - // marginal (noisy / iced) data shares with it. Names the test it failed, or nullptr. - const auto out_of_bounds = [&](const double bm[2], const double ds[1], - const double q1c[3], const double q2c[3]) -> const char * { - if (std::fabs(ds[0] - dist0) >= 0.01 * dist0) - return "the distance moved more than 1 %"; - double len_shift = 0.0, ang_shift = 0.0; - for (int j = 0; j < 3; ++j) { + // How far a geometry has moved from nominal, in units of one step. + const auto beyond_a_step = [&](const double ds[1], const double q1c[3]) { + double len_shift = 0.0; + for (int j = 0; j < 3; ++j) len_shift = std::max(len_shift, std::fabs(q1c[j] - p1_0[j]) / std::max(1e-9, p1_0[j])); - ang_shift = std::max(ang_shift, std::fabs(q2c[j] - p2_0[j])); + return std::fabs(ds[0] - dist0) >= STEP_BOUND * dist0 || len_shift >= STEP_BOUND; + }; + + // Ask the fit the way it has always been asked: ONE solve, seeded at nominal, in the + // wide box. A move that stays inside the bound is the fit this code has always + // committed and is committed unchanged - which is every crystal on the corpus but a + // handful, so nothing else about their runs moves either. + double bm_f[2] = {beam_x0, beam_y0}, ds_f[1] = {dist0}; + double rv_f[3], q0_f[3], q1_f[3], q2_f[3]; + for (int j = 0; j < 3; ++j) { + rv_f[j] = ax0[j]; q0_f[j] = p0_0[j]; q1_f[j] = p1_0[j]; q2_f[j] = p2_0[j]; + } + bool convJ = solve_joint(FIT, /*wide=*/true, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f); + double cv_ref = convJ ? joint_cost(VAL, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f, + pos_ref, exc_ref) + : cv_nom; + int steps = 1; + + // A larger move is no longer refused for its size - but it is not taken as one solve + // asserted it either. It is re-fitted as a WALK: each step seeded where the last one + // arrived, each bounded by the trust region, and each required to lower the held-out + // residual, so the size is arrived at by a sequence of checked steps. The walk stops + // when a step stops paying, which is the same cross-validated quantity the commit gate + // decides on rather than a new rule, and a walk that never stops is refused below. + const bool walked = convJ && beyond_a_step(ds_f, q1_f); + if (walked) { + bm_f[0] = beam_x0; bm_f[1] = beam_y0; ds_f[0] = dist0; + for (int j = 0; j < 3; ++j) { + rv_f[j] = ax0[j]; q0_f[j] = p0_0[j]; q1_f[j] = p1_0[j]; q2_f[j] = p2_0[j]; } - if (len_shift >= 0.01) - return "a cell length moved more than 1 %"; + convJ = false; cv_ref = cv_nom; steps = 0; + for (int step = 0; step < MAX_STEPS; ++step) { + double bm_s[2] = {bm_f[0], bm_f[1]}, ds_s[1] = {ds_f[0]}; + double rv_s[3], q0_s[3], q1_s[3], q2_s[3]; + for (int j = 0; j < 3; ++j) { + rv_s[j] = rv_f[j]; q0_s[j] = q0_f[j]; q1_s[j] = q1_f[j]; q2_s[j] = q2_f[j]; + } + if (!solve_joint(FIT, /*wide=*/false, bm_s, ds_s, rv_s, q0_s, q1_s, q2_s)) + break; + double pos_s = 0.0, exc_s = 0.0; + const double cv_s = joint_cost(VAL, bm_s, ds_s, rv_s, q0_s, q1_s, q2_s, + pos_s, exc_s); + if (step > 0 && !(cv_s < cv_ref)) + break; + bm_f[0] = bm_s[0]; bm_f[1] = bm_s[1]; ds_f[0] = ds_s[0]; + for (int j = 0; j < 3; ++j) { + rv_f[j] = rv_s[j]; q0_f[j] = q0_s[j]; q1_f[j] = q1_s[j]; q2_f[j] = q2_s[j]; + } + cv_ref = cv_s; pos_ref = pos_s; exc_ref = exc_s; + convJ = true; + ++steps; + } + } + + // Commit only for a credible move: each free cell angle under 1 deg (ANGLE_BOUND_DEG) + // and the beam inside the bound measured from whichever centre anything already + // believes is nearer. The distance and the cell lengths are no longer bounded here - + // they are bounded one step at a time by the walk above, and a move of more than one + // step is ratified by re-indexing at where it arrived (result.large_move), which is the + // discriminator the size of the move was standing in for. Names the test it failed, or + // nullptr. + const auto out_of_bounds = [&](const double bm[2], const double q2c[3]) -> const char * { + double ang_shift = 0.0; + for (int j = 0; j < 3; ++j) + ang_shift = std::max(ang_shift, std::fabs(q2c[j] - p2_0[j])); if (ang_shift >= ANGLE_BOUND_DEG * PI / 180.0) return "a cell angle moved more than 1 deg"; const double from_nominal = std::hypot(bm[0] - beam_x0, bm[1] - beam_y0); @@ -594,10 +667,17 @@ PostRefineResult PostRefineRotationGeometry(const std::vector= MAX_STEPS ? "the fit used every step it was given and never settled" : !(cv_ref < 0.98 * cv_nom) ? "the held-out residual did not improve enough" : !(pos_ref < pos_nom) ? "the held-out positional residual did not improve" : !(exc_ref < exc_nom) ? "the held-out excitation residual did not improve" - : out_of_bounds(bm_f, ds_f, q1_f, q2_f); + : out_of_bounds(bm_f, q2_f); // What the half that was held out earned the right to fit is re-fitted on all of it, and // that second solve moves the geometry - so the bounds are asked again of what will @@ -607,10 +687,19 @@ PostRefineResult PostRefineRotationGeometry(const std::vector #include +#include #include #include "../../common/DiffractionGeometry.h" @@ -40,6 +41,15 @@ struct PostRefineResult { // together - so these two are always equal. Both are kept because the report names them separately. bool cell_refined = false; // the joint fit passed cross-validation (crystal half) bool detector_refined = false; // the joint fit passed cross-validation (detector half) + // The fit moved the distance or a cell length by more than one step of its trust region (see + // STEP_BOUND in PostRefine.cpp), so it walked rather than settling beside the geometry it started + // from. Such a move is not refused for its size - it is RATIFIED by re-integrating and re-indexing + // at it, which a second lattice does not survive and a wrong header does (see RunAllPasses). + bool large_move = false; + // Which test refused the fit and what it wanted, where it was refused; empty where it committed. + // Nothing else says it - the geometry that comes out of a refusal is the nominal one whatever + // refused it, so a silent refusal cannot be told from a fit that never wanted to move. + std::string refused_reason; // GONIOMETER ROTATION SCALE: the factor by which the stage actually turned relative to the angle // stored in the file (which is the COMMANDED value, hence a stage calibration error is invisible in // the header). Fitted after the joint fit as a single free parameter, with the crystal and the axis diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index c23ec23d1..465e81d68 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -457,6 +457,12 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyText("POSTREFINE_BEAM_CENTRE", fmt::format("{:.2f} {:.2f} -> {:.2f} {:.2f}", pr.beam_x_before_px, pr.beam_y_before_px, pr.beam_x_after_px, pr.beam_y_after_px))); + // Which test refused the fit, where one did. Every refusal leaves the same geometry - the + // header's - so without this the report cannot tell a fit that was thrown away from one that + // never wanted to move, and a refusal is the one case where the number worth reading is the + // one that was thrown away. + if (!pr.refused_reason.empty()) + Add(s, KeyText("POSTREFINE_REFUSED", pr.refused_reason)); Add(s, KeyReal("GONIOMETER_ROTATION_SCALE", pr.rotation_scale, "{:.5f}")); Add(s, KeyBool("GONIOMETER_ROTATION_SCALE_SUSPECT", pr.rotation_scale_suspect)); if (pr.rotation_scale_suspect) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 1fc8b52cc..4b5adf964 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1969,12 +1969,68 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { logger.Info("Rotation two-pass geometry post-refinement: second pass (refined geometry, canonical) -> {}_*", base_prefix); config_.output_prefix = base_prefix; // the refined pass is the canonical result (no _02 suffix) + // The canonical pass measures the post-refinement again where the pre-pass moved the geometry + // (it applies nothing), so that the walk below has something to decide on. + postrefine_probe_ = prepass_detector_geometry_.has_value(); auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + + // Walk the geometry to the fit's fixed point. The post-refinement is fitted to reflections + // integrated at the very geometry it is correcting, so one fit reaches only as far as those + // reflections see: on a header whose distance is several percent long, the pre-pass absorbs the + // error into an inflated cell, and from that self-consistent pair the fit finds a fraction of it + // and stops. The canonical pass has just re-integrated and re-indexed at what it committed, so + // asking the same fit there is the next step of the same walk - and that re-indexing is also the + // RATIFICATION the old size bound stood in for: a second lattice does not index better at a + // geometry it did not come from, a wrong header does, and a pass that lost the crystal is thrown + // out by the quality guard below, which returns the run to the header geometry. + // + // What starts the walk is a move of more than one step: that is the fit saying the geometry is + // somewhere else entirely rather than a fraction of a percent away, and it leaves every run + // whose fit settles beside its header at the two passes it always had. Once the run IS walking + // it keeps walking while the fit still commits, because the steps do not get smaller in + // proportion to what is left: measured on a header 3.9 % long, the first round took a third of + // the error and each of the next took a third of the rest. The walk is bounded rather than + // stopped at a size, because there is no size that separates the last step of a real walk from + // the first step of a runaway - and it does not run away: started BELOW the optimum the same fit + // turns back up, so it settles around the answer instead of past it. + constexpr int MAX_GEOMETRY_ROUNDS = 4; + std::optional> adopted_geometry; + int geometry_rounds = 0; + while (!cancelled_ && geometry_rounds < MAX_GEOMETRY_ROUNDS && prepass_detector_geometry_ + && pass2.post_refine && pass2.post_refine->detector_refined + && (pass2.post_refine->large_move || geometry_rounds > 0)) { + adopted_geometry = *prepass_detector_geometry_; + const auto &g = *adopted_geometry; + logger.Info("Two-pass: the post-refinement at the adopted geometry still moves it " + "(distance {:.3f} -> {:.3f} mm) - re-integrating and re-indexing at what it " + "asks for (round {} of at most {})", pass2.post_refine->distance_before_mm, + pass2.post_refine->distance_after_mm, geometry_rounds + 1, + MAX_GEOMETRY_ROUNDS); + experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]) + .PoniRot1_rad(g[3]).PoniRot2_rad(g[4]); + pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + ++geometry_rounds; + } + postrefine_probe_ = false; + // The post-refinement is measured by pass 1 and consumed by pass 2, so carry it onto whichever - // result is returned - it is a result of the run, not of the pass that happened to fit it. + // result is returned - it is a result of the run, not of the pass that happened to fit it. Where + // the walk took rounds they are ONE geometry decision, so what the report names is where the + // geometry started and where the run left it, not the first step of the way there. + if (geometry_rounds > 0 && pass1.post_refine && adopted_geometry) { + pass1.post_refine->beam_x_after_px = (*adopted_geometry)[0]; + pass1.post_refine->beam_y_after_px = (*adopted_geometry)[1]; + pass1.post_refine->distance_after_mm = (*adopted_geometry)[2]; + } pass2.post_refine = pass1.post_refine; - pass2.pass_number = 2; - pass2.pass_count = 2; + pass2.pass_number = 2 + geometry_rounds; + pass2.pass_count = 2 + geometry_rounds; + if (geometry_rounds > 0) + pass2.pass_decision = fmt::format( + "post-refined geometry adopted over {} rounds: the fit was still walking at each of " + "them, and the detector distance ended at {:.3f} mm against the header's {:.3f}", + geometry_rounds + 1, (*adopted_geometry)[2], + pass1.post_refine ? pass1.post_refine->distance_before_mm : 0.0); // r2 = r1 + 2 is also the inner edge of the background ring, so widening the signal disk pushes // that ring out into the neighbours; a reflection whose ring is left with five or fewer clean @@ -4901,14 +4957,23 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // apply the detector geometry from this integration. That is all the pre-pass is for: everything // the second pass consumes from it is fitted here, and the space group is left to the second pass // (see the search gate below). - if (geometry_prepass && rsm && experiment_.IsRotationIndexing()) { + // postrefine_probe_ brings the post-refinement into a CANONICAL pass as well, where the + // pre-pass moved the geometry: the fit it committed was made on reflections integrated at the + // geometry BEFORE the move, and whether it has finished walking can only be asked of + // reflections integrated at the geometry after it. The probe measures - it applies nothing and + // leaves the smoothed mosaicity and the pre-pass lattice alone, which stay pass 1's for every + // pass that consumes them; RunAllPasses decides what to do with the answer. + if ((geometry_prepass || postrefine_probe_) && rsm && experiment_.IsRotationIndexing()) { const auto &outcomes = indexer->GetIntegrationOutcome(); - prepass_mosaicity_.assign(outcomes.size(), NAN); - for (size_t o = 0; o < outcomes.size(); ++o) - if (outcomes[o].mosaicity_deg) prepass_mosaicity_[o] = *outcomes[o].mosaicity_deg; + if (geometry_prepass) { + prepass_mosaicity_.assign(outcomes.size(), NAN); + for (size_t o = 0; o < outcomes.size(); ++o) + if (outcomes[o].mosaicity_deg) prepass_mosaicity_[o] = *outcomes[o].mosaicity_deg; + } if (result.consensus_cell.has_value()) { if (const auto rot = indexer->FinalizeRotationIndexing(); rot && rot->axis) { - prepass_result_ = *rot; // full result (lattice+geom+axis), for the supercell fallback + if (geometry_prepass) + prepass_result_ = *rot; // full result (lattice+geom+axis), for the supercell fallback PostRefineSettings prs; prs.crystal_system = rot->search_result.system; prs.num_threads = static_cast(config_.nthreads); @@ -4940,14 +5005,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b rot->geom.GetPoniRot1_rad(), rot->geom.GetPoniRot2_rad()}; logger.Info("Two-pass: post-refined detector geometry (distance {:.3f} mm, " - "beam {:.2f},{:.2f} px, tilt {:.4f},{:.4f} deg) will drive the " - "second integration pass", + "beam {:.2f},{:.2f} px, tilt {:.4f},{:.4f} deg) - {}", pr.distance_after_mm, pr.beam_x_after_px, pr.beam_y_after_px, rot->geom.GetPoniRot1_rad() * 180.0 / PI, - rot->geom.GetPoniRot2_rad() * 180.0 / PI); + rot->geom.GetPoniRot2_rad() * 180.0 / PI, + geometry_prepass ? "it will drive the second integration pass" + : "what the fit still asks for at the geometry " + "this pass ran at"); } else { - logger.Info("Two-pass: geometry post-refine committed no detector change " - "- the second pass reproduces the first"); + logger.Info("Two-pass: geometry post-refine committed no detector change{}", + geometry_prepass ? " - the second pass reproduces the first" : ""); } // DECISION POINT for the goniometer rotation scale. It does not ride on pr.ok - the // scale is its own cross-validated fit and the crystals that have a stage fault are @@ -4957,7 +5024,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // would silently damage it, so the asymmetry is deliberate: committed only when the // fit is both cross-validated and outside the tolerance. A manual --rotation-scale is // already on the goniometer and is left alone. - if (pr.rotation_scale_suspect && !config_.rotation_scale.has_value()) + // Only the pre-pass fits the rotation scale. It is applied from the goniometer + // the file came with, so a scale measured again on angles that have already been + // corrected once is not a correction that can be applied on top of that one. + if (pr.rotation_scale_suspect && !config_.rotation_scale.has_value() && geometry_prepass) prepass_rotation_scale_ = static_cast(pr.rotation_scale); } } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 642ff9f29..3fd58f2ba 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -501,6 +501,11 @@ class Rugnux { // of itself as soon as its search merge is in, which is where it learns that its report and its // files are about to be superseded. std::optional quality_guard_pass1_; + // Whether a canonical pass also MEASURES the post-refinement (it never applies it). The fit the + // pre-pass committed was made on reflections integrated at the geometry before the move, so + // whether it has finished walking is a question only reflections integrated after the move can + // answer; RunAllPasses turns this on for the passes that answer it. + bool postrefine_probe_ = false; // Pre-scan: read a spread sample of frames and take three things off them - the shadow of the beam // stop and its holder, added to the pixel mask (config_.detect_beam_stop), the beam centre