diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 264bf1cb9..325be8583 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -139,7 +139,7 @@ and a run that produced a textbook data set read identically for their first thr `SYMMETRY_AMBIGUITY`, `CENTERING_UNTESTED`, `UNUSABLE_MERGE`, `LOW_COMPLETENESS`, `SWEEP_GAPS`, `GONIO_SCALE`, `SPINDLE_CAP`, `ANISOTROPY`, `TWINNING`, `PSEUDO_TRANSLATION`, `LATTICE_TRANSLATION`, `MODEL_HAND`, `MODEL_NOT_VALIDATED`, `CANCELLED`, `RESOLUTION_FIT`, - `FLIGHT_PATH`. `NONE` when nothing fired. A code appears if and only if its + `FLIGHT_PATH`, `GEOMETRY_NOT_CONVERGED`. `NONE` when nothing fired. A code appears if and only if its warning fired, so the flags and the `WARNING:` lines are two renderings of one list — the closed type for machinery, the open sentence for a person. - **`WARNING_COUNT=`** and the **`WARNING:`** lines follow, in the same section. They are what they diff --git a/rugnux/ReportDocument.h b/rugnux/ReportDocument.h index 48b42c917..bc587c577 100644 --- a/rugnux/ReportDocument.h +++ b/rugnux/ReportDocument.h @@ -43,6 +43,7 @@ namespace PathologyCode { // in how exact it is. ANISOTROPY and RESOLUTION_FIT are already shared this way. constexpr const char *LATTICE_TRANSLATION = "LATTICE_TRANSLATION"; // a translation the cell does not declare constexpr const char *HARMONIC_CONTAMINATION = "HARMONIC_CONTAMINATION"; // the beam carries a higher harmonic + constexpr const char *GEOMETRY_NOT_CONVERGED = "GEOMETRY_NOT_CONVERGED"; // the geometry walk ran out of rounds } // A closed type plus an open description. The type is what machinery switches on; the text is what a diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index cbb2ae8ba..0a976d95d 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -385,6 +385,11 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, if (result.pass_count > 1) { Add(s, KeyText("PASS", fmt::format("{} of {}", result.pass_number, result.pass_count))); Add(s, KeyText("PASS_DECISION", result.pass_decision, true)); + if (result.geometry_not_converged) + Warn(doc, PathologyCode::GEOMETRY_NOT_CONVERGED, + "the detector geometry post-refinement was still moving when the run stopped walking " + "it: every round lowered the held-out residual, so the geometry this result was " + "integrated at is the best one reached, not a settled one (see PASS_DECISION)"); } Add(s, Blank()); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 4e6d247b0..a6da5bbb2 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1987,50 +1987,124 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { // 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); + // it keeps walking 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. + // + // But a round is kept only for what it REALISES, not for what the fit predicts: the held-out + // residual of the pass re-integrated and re-indexed at the new geometry has to fall below the + // best round's by more than the cross-validation noise of the two (HeldOutResidualFell). With + // the walk allowed twelve rounds, every round's fit predicted a gain the re-integration did not + // show - the realised residual stayed flat, or passed a minimum and rose again, while the + // distance crept on along the distance/cell-scale degeneracy - so "the fit still commits" never + // ended a walk, and the cap decided where on the drift the run stopped. A round that does not + // pay ends the walk and the run goes back to the best round it reached. So does a fit whose + // committed cell the re-indexing pushes back the same way twice running (ReindexPushesCellBack): + // the fit keeps asking for a cell the indexer at the fit's own geometry does not return, and + // nothing ratifies it. A walk still paying when the rounds run out keeps its best round - the + // last - and says it did not converge; the quality guard below judges it against the header + // pass as it judges every refined pass. + constexpr int MAX_GEOMETRY_ROUNDS = 8; + const auto current_geometry = [&] { + return std::array{experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), + experiment_.GetDetectorDistance_mm(), + experiment_.GetPoniRot1_rad(), experiment_.GetPoniRot2_rad()}; + }; + const auto set_geometry = [&](const std::array &g) { experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]) .PoniRot1_rad(g[3]).PoniRot2_rad(g[4]); + }; + int geometry_rounds = 0; // passes run at a geometry the walk moved to + int walk_passes = 0; // every pass the walk ran, the return to its best round included + int best_round = 0; + std::array best_geometry = current_geometry(); + PostRefineResult best_fit; // the fit measured in the best round's pass: its realised residual + int last_pushback = 0; + bool geometry_not_converged = false; + std::string walk_stop; // why the walk ended, where anything but the fit settling ended it + std::string walk_rounds; // distance / realised residual per round, for the pass decision + if (pass2.post_refine) { + best_fit = *pass2.post_refine; + walk_rounds = fmt::format("{:.3f} mm / {:.3e}", best_geometry[2], best_fit.held_out_before); + } + while (!cancelled_ && prepass_detector_geometry_ + && pass2.post_refine && pass2.post_refine->detector_refined + && (pass2.post_refine->large_move || geometry_rounds > 0)) { + if (geometry_rounds == MAX_GEOMETRY_ROUNDS) { + geometry_not_converged = true; + walk_stop = fmt::format("the fit still moved the geometry after {} rounds, every one of " + "which lowered the realised residual", MAX_GEOMETRY_ROUNDS); + logger.Warning("Two-pass: {} - the walk did not converge; the run keeps its last round " + "and the quality guard judges it against the header pass", walk_stop); + break; + } + const PostRefineResult fit = *pass2.post_refine; + const std::array g = *prepass_detector_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 {})", fit.distance_before_mm, + fit.distance_after_mm, geometry_rounds + 1, MAX_GEOMETRY_ROUNDS); + set_geometry(g); pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); ++geometry_rounds; + ++walk_passes; + const double realised = pass2.post_refine ? pass2.post_refine->held_out_before : NAN; + walk_rounds += fmt::format(", {:.3f} mm / {:.3e}", g[2], realised); + if (!pass2.post_refine || !HeldOutResidualFell(best_fit, *pass2.post_refine)) { + walk_stop = fmt::format( + "round {} did not lower the realised held-out residual by more than its noise " + "({:.3e} against {:.3e}, standard errors {:.1e} and {:.1e})", geometry_rounds, + realised, best_fit.held_out_before, + pass2.post_refine ? pass2.post_refine->held_out_before_se : NAN, + best_fit.held_out_before_se); + logger.Info("Two-pass: {}", walk_stop); + break; + } + best_round = geometry_rounds; + best_geometry = g; + best_fit = *pass2.post_refine; + const int pushback = ReindexPushesCellBack(fit, *pass2.post_refine); + if (pushback != 0 && pushback == last_pushback) { + walk_stop = fmt::format( + "re-indexing pushed the cell length {} back against the fit for the second " + "round running, so nothing ratifies the cell the fit asks for", + "abc"[std::abs(pushback) - 1]); + logger.Info("Two-pass: {} - the walk stops at round {}", walk_stop, geometry_rounds); + break; + } + last_pushback = pushback; } postrefine_probe_ = false; + // The canonical files were written by the last round; where that is not the best one, the best + // one is run again to write them. + if (!cancelled_ && geometry_rounds > best_round) { + logger.Info("Two-pass: going back to round {} of the geometry walk (distance {:.3f} mm), the " + "last one that lowered the realised residual", best_round, best_geometry[2]); + set_geometry(best_geometry); + pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + ++walk_passes; + } // 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. 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]; + if (best_round > 0 && pass1.post_refine) { + pass1.post_refine->beam_x_after_px = best_geometry[0]; + pass1.post_refine->beam_y_after_px = best_geometry[1]; + pass1.post_refine->distance_after_mm = best_geometry[2]; } pass2.post_refine = pass1.post_refine; - pass2.pass_number = 2 + geometry_rounds; - pass2.pass_count = 2 + geometry_rounds; + pass2.pass_number = 2 + walk_passes; + pass2.pass_count = 2 + walk_passes; 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); + "post-refined geometry walked {} round{} and kept round {}, detector distance {:.3f} mm " + "against the header's {:.3f}: {}; distance / realised held-out residual by round: {}", + geometry_rounds, geometry_rounds == 1 ? "" : "s", best_round, best_geometry[2], + pass1.post_refine ? pass1.post_refine->distance_before_mm : 0.0, + walk_stop.empty() ? "the fit at the kept geometry commits no further move" : walk_stop, + walk_rounds); // 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 @@ -2157,6 +2231,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { pass2.pass_decision = (prepass_detector_geometry_ || prepass_rotation_scale_) ? "post-refined geometry adopted" : "the post-refinement committed no geometry change, so this pass reproduces the first"; + pass2.geometry_not_converged = geometry_not_converged; return pass2; } return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index d56907983..91312251a 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -361,6 +361,9 @@ struct ProcessResult { int pass_number = 1; int pass_count = 1; std::string pass_decision; + // The rotation geometry walk used every round it was given and its fit still asked for a move: the + // geometry this result was integrated at is the best round the walk reached, not a settled one. + bool geometry_not_converged = false; // The radiation-damage report (rotation), the same text the log prints. Empty when not measured. std::string radiation_damage_text; diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index 594306c74..12d61a22b 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -793,3 +793,22 @@ TEST_CASE("ResultReport_TranslationalNCS", "[Diagnostics]") { CHECK(lat.find("\nPATHOLOGY_FLAGS= LATTICE_TRANSLATION\n") != std::string::npos); CHECK(lat.find("\nTNCS_SUBLATTICE=") == std::string::npos); } + +TEST_CASE("ResultReport_GeometryNotConverged", "[Diagnostics]") { + DiffractionExperiment x(DetJF(1)); + ProcessResult result; + result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f, + .alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f}; + result.pass_number = 10; + result.pass_count = 10; + result.pass_decision = "post-refined geometry walked 8 rounds and kept round 8"; + + CHECK(RenderResultReport("p", "in.h5", x, result).find("GEOMETRY_NOT_CONVERGED") == std::string::npos); + + result.geometry_not_converged = true; + const auto text = RenderResultReport("p", "in.h5", x, result); + const auto flags = text.find("\nPATHOLOGY_FLAGS= "); + REQUIRE(flags != std::string::npos); + CHECK(text.substr(flags, text.find('\n', flags + 1) - flags).find("GEOMETRY_NOT_CONVERGED") + != std::string::npos); +}