diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index 6c2bce937..ff43c94bc 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -53,7 +53,7 @@ constexpr double ANGLE_BOUND_DEG = 1.0; // 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; +constexpr double STEP_BOUND = POSTREFINE_STEP_BOUND; // ... 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. @@ -424,6 +424,20 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations, const UnitCell r0 = reference_latt.GetUnitCell(); const double beam_x0 = nominal_geom.GetBeamX_pxl(), beam_y0 = nominal_geom.GetBeamY_pxl(); const double dist0 = nominal_geom.GetDetectorDistance_mm(); + // The beam the fit refines is the PONI, the foot of the normal from the sample. On a detector + // swung out on a 2theta arm the direct beam lies a distance-proportional way from it, so at a + // fixed direct beam the PONI moves with the distance (7 px per mm at 35 deg on a 0.1 mm pixel), + // and a bound on the PONI caps the distance instead of the beam. What the bound is about is + // where the beam hits the detector, which is also what a measured centre measures, so it is + // applied there: direct beam = PONI + beam_per_mm * distance. Zero on an untilted detector. + double beam_per_mm[2] = {0.0, 0.0}; + { + const auto [dbx, dby] = nominal_geom.GetDirectBeam_pxl(); + if (std::isfinite(dbx) && std::isfinite(dby)) { + beam_per_mm[0] = (dbx - beam_x0) / dist0; + beam_per_mm[1] = (dby - beam_y0) / dist0; + } + } // Deterministic split of the reflections into a fit half and a held-out half. Avalanche-mix the // hkl hash so the split bit is decorrelated from the LSB - a plain h+k+l parity collides with the @@ -670,11 +684,15 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations, // the bound and is then refused for being there. p.SetParameterLowerBound(ds, 0, wide ? dist0 * 0.95 : ds_seed * (1.0 - STEP_BOUND)); p.SetParameterUpperBound(ds, 0, wide ? dist0 * 1.05 : ds_seed * (1.0 + STEP_BOUND)); + // The PONI box reaches every PONI whose direct beam the gate below can accept, at any + // distance the walk can reach (MAX_STEPS steps). for (int j = 0; j < 2; ++j) { const double b0 = (j == 0) ? beam_x0 : beam_y0; - const double m = settings.measured_beam_px ? (*settings.measured_beam_px)[j] : b0; - p.SetParameterLowerBound(bm, j, std::min(b0, m) - BEAM_BOUND_PXL); - p.SetParameterUpperBound(bm, j, std::max(b0, m) + BEAM_BOUND_PXL); + const double m = settings.measured_beam_px + ? (*settings.measured_beam_px)[j] - beam_per_mm[j] * dist0 : b0; + const double reach = std::fabs(beam_per_mm[j]) * dist0 * MAX_STEPS * STEP_BOUND; + p.SetParameterLowerBound(bm, j, std::min(b0, m) - BEAM_BOUND_PXL - reach); + p.SetParameterUpperBound(bm, j, std::max(b0, m) + BEAM_BOUND_PXL + reach); } for (int j = 0; j < 3; ++j) { p.SetParameterLowerBound(rv, j, ax0[j] - 0.05); @@ -770,6 +788,11 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations, bool convJ = solve_joint(FIT, /*wide=*/true, /*fix_distance=*/false, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f, distance_corr); const bool held_converged = held_fit.valid() && held_fit.get(); + if (convJ) { + result.distance_free_mm = ds_f[0]; + result.beam_x_free_px = bm_f[0]; + result.beam_y_free_px = bm_f[1]; + } double exc_ref_se = NAN; double cv_ref = convJ ? joint_cost(VAL, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f, pos_ref, exc_ref, nullptr, &exc_ref_se) @@ -876,16 +899,19 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations, // 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 * { + const auto out_of_bounds = [&](const double bm[2], const double ds[1], + 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); + const double db[2] = {bm[0] + beam_per_mm[0] * ds[0], bm[1] + beam_per_mm[1] * ds[0]}; + const double from_nominal = std::hypot(db[0] - beam_x0 - beam_per_mm[0] * dist0, + db[1] - beam_y0 - beam_per_mm[1] * dist0); const double from_measured = settings.measured_beam_px - ? std::hypot(bm[0] - (*settings.measured_beam_px)[0], - bm[1] - (*settings.measured_beam_px)[1]) + ? std::hypot(db[0] - (*settings.measured_beam_px)[0], + db[1] - (*settings.measured_beam_px)[1]) : std::numeric_limits::infinity(); if (std::min(from_nominal, from_measured) >= BEAM_BOUND_PXL) return "the beam moved further than the bound from every centre anything believes"; @@ -949,7 +975,7 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations, ? "the held-out positional residual did not fall by more than its noise" : !(change.exc_n < 2 || change.exc <= change.exc_se) ? "the held-out excitation residual rose by more than its noise" - : out_of_bounds(bm_f, q2_f); + : out_of_bounds(bm_f, ds_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 @@ -971,7 +997,7 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations, solve_joint(ALL, !walked, hold_distance, beam, dist, axv, p0, p1, p2); // commit: re-fit on all data ds_log[0] = dist[0]; bm_log[0] = beam[0]; bm_log[1] = beam[1]; for (int j = 0; j < 3; ++j) { q1_log[j] = p1[j]; q2_log[j] = p2[j]; } - refused = out_of_bounds(beam, p2); + refused = out_of_bounds(beam, dist, p2); if (refused != nullptr) { beam[0] = beam_x0; beam[1] = beam_y0; dist[0] = dist0; for (int j = 0; j < 3; ++j) { diff --git a/image_analysis/geom_refinement/PostRefine.h b/image_analysis/geom_refinement/PostRefine.h index 7d4c40d20..2f7d8c360 100644 --- a/image_analysis/geom_refinement/PostRefine.h +++ b/image_analysis/geom_refinement/PostRefine.h @@ -31,6 +31,10 @@ // Where there are more reflections than the fit's caps, the strongest by I/sigma are the ones kept; // inside the fit every reflection carries the same weight. +// How far one step of the joint fit may move the distance or a cell length, relative (see STEP_BOUND in +// PostRefine.cpp, where it is explained). +constexpr double POSTREFINE_STEP_BOUND = 0.01; + struct PostRefineResult { bool ok = false; DiffractionGeometry geom; // refined (distance / beam left as configured) @@ -53,6 +57,11 @@ struct PostRefineResult { // ever - because freeing it did not lower the held-out excitation residual by more than that // residual's noise (see the joint fit in PostRefine.cpp). The report names it. bool distance_held = false; + // Where the fit with the distance free went, held or not (NaN where it did not converge): the + // detector move these data ask for at the geometry they were integrated at. + double distance_free_mm = std::numeric_limits::quiet_NaN(); + double beam_x_free_px = std::numeric_limits::quiet_NaN(); + double beam_y_free_px = std::numeric_limits::quiet_NaN(); // The held-out residual at the geometry the fit STARTED from - what the integration this fit was // handed actually realises there, before the fit moves anything - and its standard error over the // held-out values. Two passes at two geometries are compared on this, not on the fit's own diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index a8bf7753a..ef411059c 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -1820,7 +1820,14 @@ class Rugnux { std::optional gonio_snapshot; std::array header_geometry{}; RotationScaleOutcome scale_walk; // where the walk has ended before they start + int distance_passes = 0; // the pre-passes the fixed-point distance ran (FixedPointDistance) }; + // The detector distance as the fixed point of the integration, where pass 1 held the header but its + // free fit asked for a move of a step or more (see RunAllPasses). Returns the geometry to start the + // refined passes from, or nothing where the header stands; `passes` counts the pre-passes it ran. + std::optional> FixedPointDistance(RugnuxObserver *observer, const ProcessResult &pass1, + const DiffractionExperiment &first_pass_start, + int &passes); // The refined passes of the rotation two-pass, from the second pass on. Empty where they ran beside // the rotation-scale walk and it overturned them. // After the refined passes: where the supercell probe's class stands clear, integrate and merge the diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index 56745c722..e6b1bea90 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -627,12 +627,21 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { experiment_.GetDetectorDistance_mm(), experiment_.GetPoniRot1_rad(), experiment_.GetPoniRot2_rad()}; + int distance_passes = 0; + if (const auto g = FixedPointDistance(observer, pass1, experiment_before_first_pass, distance_passes)) { + prepass_detector_geometry_ = g; + pass1.post_refine->distance_after_mm = (*g)[2]; + pass1.post_refine->beam_x_after_px = (*g)[0]; + pass1.post_refine->beam_y_after_px = (*g)[1]; + pass1.post_refine->distance_held = false; + } + if (cancelled_) return pass1; if (prepass_detector_geometry_) { const auto &g = *prepass_detector_geometry_; experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]) .PoniRot1_rad(g[3]).PoniRot2_rad(g[4]); } - RefinedPassStart start{gonio_snapshot, header_geometry, {}}; + RefinedPassStart start{gonio_snapshot, header_geometry, {}, distance_passes}; // The goniometer rotation scale: the angles in the file are the commanded ones, so a stage that // turned at the wrong rate is a geometry error like any other - but one pass 1's fit reads only @@ -713,6 +722,138 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); } +// The detector distance as the fixed point of the integration. Pass 1 integrated at the header, and its +// post-refinement left the distance there - the rocking angles did not ask for the move the free fit +// wanted, or the held-out residual did not pay for it - but a fit is measured on reflections integrated at the geometry it starts from, and those +// carry that geometry with them: which reflections are predicted, and on which frames their rocking +// centroids are taken. Measured on a detector swung out on a 2theta arm with a header 4 % short, the fit +// asked for 104.4 mm from data integrated at 103 and for 108.6 mm from data integrated at 109, each a +// third of the way to where the other started, so a walk that moves to what each round asks for crawls +// and stalls short; spot-finder centroids in place of the integrated ones moved neither number. +// +// So the distance is asked for as what it is, the point where the data integrated there ask for no +// move: a few pre-passes at distances along the free fit's direction (the direct beam held where the free +// fit put it), the move each one asks for against the distance it was integrated at, and the zero of the +// straight line through them. That distance is integrated once more and adopted only for what it +// REALISES against the header pass: a lower held-out residual beyond its noise, or validation spots the +// lattice takes better - the two tests a round of the geometry walk has to pass. A header that is right +// asks for less than a step and never gets here. +std::optional> Rugnux::FixedPointDistance(RugnuxObserver *observer, const ProcessResult &pass1, + const DiffractionExperiment &first_pass_start, + int &passes) { + if (cancelled_ || config_.fast_live || !pass1.post_refine) + return std::nullopt; + const PostRefineResult &held = *pass1.post_refine; + const double d0 = experiment_.GetDetectorDistance_mm(); + if (prepass_detector_geometry_ && std::fabs((*prepass_detector_geometry_)[2] - d0) > 1e-3 * d0) + return std::nullopt; // pass 1 moved the distance itself; the walk takes it from there + const double move = held.distance_free_mm - d0; + if (!std::isfinite(move) || std::fabs(move) < POSTREFINE_STEP_BOUND * d0) + return std::nullopt; + Logger logger("Rugnux"); + if (observer) + observer->OnStep("Distance check"); + + // The direct beam the free fit put the beam at, held there at every distance tried: on a swung + // detector the PONI moves with the distance at a fixed direct beam (zero per mm on an untilted one). + double per_mm[2] = {0.0, 0.0}; + const auto [dbx, dby] = experiment_.GetDiffractionGeometry().GetDirectBeam_pxl(); + if (std::isfinite(dbx) && std::isfinite(dby)) { + per_mm[0] = (dbx - experiment_.GetBeamX_pxl()) / d0; + per_mm[1] = (dby - experiment_.GetBeamY_pxl()) / d0; + } + const double direct_beam[2] = {held.beam_x_free_px + per_mm[0] * held.distance_free_mm, + held.beam_y_free_px + per_mm[1] * held.distance_free_mm}; + const auto geometry_at = [&](double d) { + return std::array{static_cast(direct_beam[0] - per_mm[0] * d), + static_cast(direct_beam[1] - per_mm[1] * d), static_cast(d), + experiment_.GetPoniRot1_rad(), experiment_.GetPoniRot2_rad()}; + }; + // A pre-pass at a distance, on a copy of the run as pass 1 started it: it leaves nothing on the run + // but its first-pass memo. + const auto prepass_at = [&](double d) { + Rugnux probe(*this); + probe.experiment_ = first_pass_start; + if (bragg_adaptive_) + probe.experiment_.ImportBraggIntegrationSettings(*bragg_adaptive_); + const auto g = geometry_at(d); + probe.experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]); + probe.prepass_mosaicity_.clear(); + probe.prepass_detector_geometry_.reset(); + probe.prepass_result_.reset(); + ProcessResult r; + try { + r = probe.RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/true); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Info("Distance check: the pre-pass at {:.3f} mm did not complete ({})", d, e.what()); + } + for (const auto &m : probe.first_pass_memo_) + AddFirstPassMemo(m); + ++passes; + return r; + }; + const auto asked = [](const ProcessResult &r) { + if (!r.post_refine) return std::numeric_limits::quiet_NaN(); + return std::isfinite(r.post_refine->distance_free_mm) ? r.post_refine->distance_free_mm + : r.post_refine->distance_after_mm; + }; + + // Distances along the free fit's move, one move apart, until the move asked for changes sign. + std::vector> points{{d0, held.distance_free_mm}}; + std::string text = fmt::format("{:.2f} -> {:.2f}", d0, held.distance_free_mm); + for (int k = 1; k <= 4 && !cancelled_; ++k) { + const double d = d0 + k * move; + const double a = asked(prepass_at(d)); + text += fmt::format(", {:.2f} -> {:.2f}", d, a); + if (!std::isfinite(a)) + break; + points.emplace_back(d, a); + if ((a - d) * move <= 0.0) + break; + } + if (cancelled_) + return std::nullopt; + // The straight line through (distance, asked - distance), least squares, and its zero. + double sx = 0, sy = 0, sxx = 0, sxy = 0; + for (const auto &[d, a] : points) { + sx += d; sy += a - d; sxx += d * d; sxy += d * (a - d); + } + const double n = static_cast(points.size()); + const double slope = (n * sxy - sx * sy) / (n * sxx - sx * sx); + const double intercept = (sy - slope * sx) / n; + // A fixed point is taken only where the distances tried bracket it - the move asked for changed + // sign - so it is interpolated, never extrapolated. Where every distance asks for more of the same + // move the data do not pin the distance at all (measured: a crystal at 3 A on a detector 10 deg out + // asked for 5.5 mm more from every distance between 183 and 205 mm), and the header stands. + const double last_move = points.back().second - points.back().first; + if (points.size() < 2 || !(slope < 0.0) || last_move * move > 0.0) { + logger.Info("Distance check: the moves asked for at each distance integrated ({}) do not bracket a " + "fixed point - the header stands", text); + if (observer) + observer->OnNote("the header distance stands"); + return std::nullopt; + } + const double d_star = std::clamp(-intercept / slope, std::min(d0, points.back().first), + std::max(d0, points.back().first)); + const ProcessResult at_star = prepass_at(d_star); + const bool residual_fell = at_star.post_refine && HeldOutResidualFell(held, *at_star.post_refine); + const bool evidence_prefers = ValidationEvidencePrefers(pass1.validation_evidence, at_star.validation_evidence); + const bool adopt = residual_fell || evidence_prefers; + logger.Info("Distance check: integrated at -> asked for (mm): {}; fixed point {:.3f} mm (asked minus integrated " + "falls {:.2f} mm per mm); there the held-out residual is {:.3e} against the header's {:.3e} and " + "the validation spots on the lattice {} against {} - {}", text, d_star, -slope, + at_star.post_refine ? at_star.post_refine->held_out_before : NAN, held.held_out_before, + at_star.validation_evidence.on_lattice, pass1.validation_evidence.on_lattice, + adopt ? "adopted" : "the header stands"); + if (observer) + observer->OnNote(adopt ? fmt::format("detector distance {:.2f} mm, where the data ask for no move", d_star) + : std::string("the header distance stands")); + if (!adopt) + return std::nullopt; + return geometry_at(d_star); +} + // The rotation-scale walk of RunAllPasses, run beside the refined passes on a copy of the run (`check` // is what the two share): it starts once their first pass has indexed, calls them off (output_cancel) // before a refit, and says at its end whether it overturns them. @@ -918,7 +1059,7 @@ std::optional Rugnux::RefinedPasses(RugnuxObserver *observer, Pro RotationScaleOutcome scale_walk = start.scale_walk; // The probe passes below and the lattice arms after them exist ONLY to measure - see there - so // the count of passes the run made has to include them even though they write nothing. - int arm_passes = 0; + int arm_passes = start.distance_passes; // ... and apply the adopted scale, on the same measure-then-re-integrate footing as the geometry. if (prepass_rotation_scale_ && gonio_snapshot) {