Rugnux: walk the goniometer rotation scale to its fixed point, decided on the whole sweep
The pass-1 post-refinement fits the rotation scale k only on the frames the stored angles still track, and a rate error is exactly what stops them tracking the rest: on a sweep whose stage turned ~3 % slow the fit read 0.979 over 94 deg, failed its leave-a-fifth-out test and was thrown away, leaving half the frames unscaled. The pass no longer decides. Between the passes, at the pass-2 detector geometry, the lattice is indexed (index-only probe) under the stored angles and under the fitted k and scored on the validation frames of the whole sweep: share of the spots on the lattice beyond the wrong-spindle null. k is adopted only where it scores higher by more than the binomial noise of the two (ValidationEvidencePrefers - the test the beam-centre arms already used, now one function); the run then integrates and post-refines at k (post-refine-only probe), fits again on top of it and repeats until the next k no longer scores better (WalkRotationScale). The stored angles are the first hypothesis. Measured: 1.000 25.1 %, 0.97874 58.9 %, 0.97041 90.0 %, 0.97006 90.7 % (not significant) -> 0.97041 adopted. Probes restore the experiment, the pass-2 geometry, pass-1 mosaicity and the beam-centre-search flag; a probe opens no beam-centre search. Forced pass-1 results get their axis scaled; the header revert drops the scale. GONIOMETER_ROTATION_SCALE reports the adopted k (SUSPECT = adopted). The leave-a-fifth-out figure stays in the log only. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
This commit is contained in:
Claude Opus 5.5
@@ -335,7 +335,7 @@ The space group is determined **after** pass 2, on the geometry the run refined,
|
||||
|
||||
Only pass 2 is written, as the canonical `<prefix>_*` output. Pass 1's merge exists to give the guard something to judge pass 2 against, so it stops short of the parts of the merge that only fill in a file — the correction surfaces, the twinning and radiation-damage analyses, the R-free flags and the amplitudes — and writes no merged files of its own.
|
||||
|
||||
**Goniometer rotation scale (report only).** A stage that turns further than it was commanded to leaves no trace in the file, because the stored $\omega$ values *are* the commanded ones; the excess then presents as the crystal drifting, in this program and in others. The excitation residual already measures it without a new degree of freedom: it rotates by $-\phi\,\mathbf{u}$ with $\mathbf{u}$ an **unnormalised** 3-vector, so $|\mathbf{u}|$ is the factor by which the stage actually turned, and normalising the axis throws it away. It is reported, and warned about beyond 0.5 %, under its own leave-a-fifth-of-the-sweep-out check — a fold that merely soaked up noise cannot raise the flag. It is a detector, not a calibration: nothing corrects the data, and it **under-reads** the true magnitude, because the fit only sees reflections that indexed at the nominal angle and per-frame orientation refinement has already absorbed part of the error.
|
||||
**Goniometer rotation scale.** A stage that turns further than it was commanded to leaves no trace in the file, because the stored $\omega$ values *are* the commanded ones; the excess then presents as the crystal drifting, in this program and in others. The excitation residual already measures it without a new degree of freedom: it rotates by $-\phi\,\mathbf{u}$ with $\mathbf{u}$ an **unnormalised** 3-vector, so $|\mathbf{u}|$ is the factor $k$ by which the stage actually turned, and normalising the axis throws it away. Pass 1 fits $k$ as a single parameter on its rocking events, with the crystal and the axis direction held at their committed values and the angle measured from the centre of the sweep. That fit **under-reads** a real error: it only sees the frames the stored angles still track, and a rate error is exactly what stops them tracking the rest. So it is not acted on directly. Between the passes, at the detector geometry pass 2 runs at, the lattice is indexed under the stored angles and under the fitted $k$, and each is scored on the validation frames spread over the whole sweep, as the share of their spots it puts on the lattice beyond what it puts there at a wrong spindle angle. The fitted $k$ is adopted only where it scores higher by more than the binomial noise of the two scores (z = 3.29); the run then integrates and post-refines at it, fits $k$ again on top of it, and repeats until the next $k$ no longer scores better - the fixed point of the fit. Otherwise the stored angles stand. The adopted $k$ drives every later pass (prediction, integration, scaling and the reported oscillation) and is reported as `GONIOMETER_ROTATION_SCALE`, with `GONIOMETER_ROTATION_SCALE_SUSPECT= TRUE`. `--rotation-scale <k>` asserts a calibration and skips all of this.
|
||||
|
||||
### 7.6 Detector geometry from powder rings
|
||||
|
||||
|
||||
@@ -1080,50 +1080,19 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations,
|
||||
const double k_fit = solve_scale(-1);
|
||||
result.rotation_scale = k_fit;
|
||||
|
||||
// ----- Whether to COMMIT it. A stage fault is rare - 36 of 37 rotation datasets sit at 1.0000
|
||||
// on a direct scan - and a 1 % angle correction applied to a healthy dataset would damage it
|
||||
// silently, so every test below has to pass.
|
||||
// Preconditions: below these the fit is reported but never acted on. Under ~30 deg of sweep k
|
||||
// entangles with the axis direction and 10-20 deg truncations of a perfect dataset wander by
|
||||
// +-0.6 %; a screening wedge must not trigger a correction.
|
||||
constexpr int MIN_SCALE_EVENTS = 5000;
|
||||
constexpr double MIN_SCALE_SWEEP_DEG = 30.0;
|
||||
// T1 significance: 0.5 % is 18 sigma on the between-dataset scatter of healthy stages
|
||||
// (robust sd 2.8e-4) and still 3.5x below the one measured fault.
|
||||
constexpr double ROTATION_SCALE_TOL = 0.005;
|
||||
// T2 relevance: the misorientation the error produces at each end of the sweep. A large k over
|
||||
// a short sweep moves nothing and is not worth correcting.
|
||||
constexpr double MIN_SCALE_END_ERROR_DEG = 0.5;
|
||||
// T3 uniformity: a stage error is a ramp present in EVERY part of the sweep, so dropping any
|
||||
// fifth of it must leave the same k. A second lattice that dominates ONE END of the sweep -
|
||||
// exactly what happens where the primary stops indexing - fakes a k indistinguishable from a
|
||||
// real fault on T1 and T2, and is the reason this test is not optional. It replaces the
|
||||
// hkl-hash split used elsewhere here, which cannot see it: both halves of that split sit at
|
||||
// the same angles, so anything structured in phi survives in both folds.
|
||||
constexpr double MIN_SCALE_JACKKNIFE_FRAC = 0.5;
|
||||
const double end_error_deg = std::fabs(k_fit - 1.0) * sweep_deg / 2.0;
|
||||
const bool enough_data = static_cast<int>(n_events) >= MIN_SCALE_EVENTS
|
||||
&& sweep_deg >= MIN_SCALE_SWEEP_DEG;
|
||||
const bool big_enough = enough_data && std::fabs(k_fit - 1.0) >= ROTATION_SCALE_TOL
|
||||
&& end_error_deg >= MIN_SCALE_END_ERROR_DEG;
|
||||
// Whether the same k comes back with any fifth of the sweep left out, as the smallest share
|
||||
// of the fitted excess the folds keep: a stage error is a ramp present in EVERY part of the
|
||||
// sweep. Reported, not acted on - the fit only sees the frames the angles it was measured at
|
||||
// still track, and a rate error is exactly what stops them tracking the rest, so the part
|
||||
// it sees can be too short to agree with itself. What acts on k is the caller, which walks
|
||||
// it to its fixed point and decides it on the whole sweep (rugnux WalkRotationScale).
|
||||
double jackknife = 1.0;
|
||||
if (big_enough)
|
||||
for (int f = 0; f < 5; ++f)
|
||||
jackknife = std::min(jackknife, (solve_scale(f) - 1.0) / (k_fit - 1.0));
|
||||
result.rotation_scale_suspect = big_enough && jackknife >= MIN_SCALE_JACKKNIFE_FRAC;
|
||||
for (int f = 0; f < 5; ++f)
|
||||
jackknife = std::min(jackknife, (solve_scale(f) - 1.0) / (k_fit - 1.0));
|
||||
logger.Info("Post-refine rotation SCALE: k = {:.5f} over {:.0f} deg of sweep centred on {:.1f} "
|
||||
"deg ({} events): end error {:.2f} deg, leave-a-fifth-out {:.2f} => {}",
|
||||
k_fit, sweep_deg, phi_c * 180.0 / PI, n_events, end_error_deg, jackknife,
|
||||
result.rotation_scale_suspect ? "COMMIT"
|
||||
: !enough_data ? "report only (too little sweep or too few events)"
|
||||
: "reject (kept the stored angles)");
|
||||
if (result.rotation_scale_suspect)
|
||||
logger.Warning("Goniometer rotation scale looks off by {:+.2f} % (fitted {:.5f}): the stage "
|
||||
"appears to have turned {} than the angles stored in the file, which are the "
|
||||
"COMMANDED values. This is a hardware calibration fault, not a data problem - "
|
||||
"left uncorrected it inflates mosaicity, biases the cell and loses "
|
||||
"high-resolution reflections",
|
||||
100.0 * (k_fit - 1.0), k_fit, k_fit > 1.0 ? "further" : "less far");
|
||||
"deg ({} events): end error {:.2f} deg, leave-a-fifth-out {:.2f}",
|
||||
k_fit, sweep_deg, phi_c * 180.0 / PI, n_events,
|
||||
std::fabs(k_fit - 1.0) * sweep_deg / 2.0, jackknife);
|
||||
|
||||
// Assemble the committed geometry.
|
||||
result.distance_after_mm = dist[0];
|
||||
|
||||
@@ -73,8 +73,8 @@ struct PostRefineResult {
|
||||
// the header). Fitted after the joint fit as a single free parameter, with the crystal and the axis
|
||||
// direction held at their committed values. Always the fitted value; 1.0 = header and stage agree.
|
||||
double rotation_scale = 1.0;
|
||||
// Whether the fit passed every test needed to ACT on it: enough sweep and events, a significant and
|
||||
// physically relevant size, and the same k from every fifth of the sweep. Only then is it applied.
|
||||
// Whether the run ACTED on it: set by the caller where the scale, walked to its fixed point, was
|
||||
// adopted on the evidence of the whole sweep (rugnux WalkRotationScale) - never by the fit itself.
|
||||
bool rotation_scale_suspect = false;
|
||||
};
|
||||
|
||||
|
||||
+161
-51
@@ -75,6 +75,48 @@
|
||||
#include <array>
|
||||
#include <map>
|
||||
|
||||
bool ValidationEvidencePrefers(const ValidationSpotEvidence ¤t, const ValidationSpotEvidence &candidate) {
|
||||
const auto excess = [](const ValidationSpotEvidence &e) {
|
||||
return static_cast<double>(e.on_lattice - e.by_chance) / static_cast<double>(std::max<int64_t>(1, e.spots));
|
||||
};
|
||||
const auto rate_var = [](const ValidationSpotEvidence &e) {
|
||||
const double n = static_cast<double>(std::max<int64_t>(1, e.spots));
|
||||
const double p = static_cast<double>(e.on_lattice) / n;
|
||||
return p * (1.0 - p) / n;
|
||||
};
|
||||
return excess(candidate) - excess(current)
|
||||
> SPOT_BUDGET_SIGNIFICANCE_Z * std::sqrt(rate_var(candidate) + rate_var(current));
|
||||
}
|
||||
|
||||
RotationScaleWalk WalkRotationScale(double first_fit,
|
||||
const std::function<ValidationSpotEvidence(float)> &index_at,
|
||||
const std::function<std::optional<double>(float)> &refit_at,
|
||||
int max_rounds) {
|
||||
RotationScaleWalk walk;
|
||||
auto next = static_cast<float>(first_fit);
|
||||
if (next == walk.scale)
|
||||
return walk;
|
||||
const auto score = [](const ValidationSpotEvidence &e) {
|
||||
return 100.0 * static_cast<double>(e.on_lattice - e.by_chance)
|
||||
/ static_cast<double>(std::max<int64_t>(1, e.spots));
|
||||
};
|
||||
walk.evidence = index_at(walk.scale);
|
||||
walk.trail = fmt::format("{:.5f}: {:.1f}%", walk.scale, score(walk.evidence));
|
||||
for (int round = 0; round < max_rounds && next != walk.scale; ++round) {
|
||||
const ValidationSpotEvidence e = index_at(next);
|
||||
walk.trail += fmt::format(", {:.5f}: {:.1f}%", next, score(e));
|
||||
if (!ValidationEvidencePrefers(walk.evidence, e))
|
||||
break;
|
||||
walk.scale = next;
|
||||
walk.evidence = e;
|
||||
const auto fit = refit_at(walk.scale);
|
||||
if (!fit)
|
||||
break;
|
||||
next = static_cast<float>(walk.scale * *fit);
|
||||
}
|
||||
return walk;
|
||||
}
|
||||
|
||||
double MetricViolation(const UnitCell &uc, const gemmi::SpaceGroup &sg) {
|
||||
const double a = uc.a, b = uc.b, c = uc.c;
|
||||
const double d2r = PI / 180.0;
|
||||
@@ -2024,7 +2066,6 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
|
||||
const DiffractionExperiment file_experiment = experiment_;
|
||||
const auto file_mosaicity = prepass_mosaicity_;
|
||||
const auto file_geometry = prepass_detector_geometry_;
|
||||
const auto file_scale = prepass_rotation_scale_;
|
||||
const auto file_result = prepass_result_;
|
||||
|
||||
experiment_ = experiment_before_first_pass;
|
||||
@@ -2036,7 +2077,6 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
|
||||
experiment_.BeamX_pxl(alt_center[0]).BeamY_pxl(alt_center[1]);
|
||||
prepass_mosaicity_.clear();
|
||||
prepass_detector_geometry_.reset();
|
||||
prepass_rotation_scale_.reset();
|
||||
prepass_result_.reset();
|
||||
|
||||
ProcessResult alt;
|
||||
@@ -2113,7 +2153,6 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
|
||||
experiment_ = file_experiment;
|
||||
prepass_mosaicity_ = file_mosaicity;
|
||||
prepass_detector_geometry_ = file_geometry;
|
||||
prepass_rotation_scale_ = file_scale;
|
||||
prepass_result_ = file_result;
|
||||
}
|
||||
// The arm's own pass asks the same question again at its own centre; the run has already
|
||||
@@ -2152,10 +2191,90 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
|
||||
experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2])
|
||||
.PoniRot1_rad(g[3]).PoniRot2_rad(g[4]);
|
||||
}
|
||||
// ... and the goniometer rotation scale, on the same measure-then-re-integrate footing: the angles
|
||||
// in the file are the commanded ones, so a stage that ran fast is a geometry error like any other.
|
||||
// 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;
|
||||
|
||||
// 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
|
||||
// in part (see WalkRotationScale), so it is walked to the fit's fixed point here and adopted
|
||||
// only where the validation frames of the whole sweep prefer it. The walk runs at the detector
|
||||
// geometry the second pass will, so every hypothesis, the stored angles included, is scored
|
||||
// there, and at that geometry only: a probe whose angles lose the lattice must score what it
|
||||
// scores, not start the beam-centre search that a poor first pass otherwise opens. Its passes
|
||||
// only measure - an indexing probe stops once the lattice is scored, a refit stops once the
|
||||
// post-refinement has measured - and they leave nothing behind: the experiment, the geometry
|
||||
// the second pass is to run at, pass 1's mosaicity (a width in degrees fitted against the
|
||||
// stored angles) and whether the run has searched for the beam centre are put back when the
|
||||
// walk ends.
|
||||
std::string rotation_scale_walk;
|
||||
if (!cancelled_ && gonio_snapshot && pass1.post_refine && !config_.rotation_scale) {
|
||||
const DiffractionExperiment before_walk = experiment_;
|
||||
const auto geometry_before_walk = prepass_detector_geometry_;
|
||||
const auto mosaicity_before_walk = prepass_mosaicity_;
|
||||
const bool searched_before_walk = beam_center_searched_;
|
||||
beam_center_searched_ = true;
|
||||
const auto probe = [&](float k, bool index_only) {
|
||||
experiment_ = before_walk;
|
||||
experiment_.Goniometer(ScaleRotation(*before_walk.GetGoniometer(), k));
|
||||
prepass_rotation_scale_ = k;
|
||||
prepass_mosaicity_.clear();
|
||||
indexing_probe_only_ = index_only;
|
||||
postrefine_probe_ = postrefine_probe_only_ = !index_only;
|
||||
ProcessResult r;
|
||||
try {
|
||||
r = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false);
|
||||
} catch (const std::exception &e) {
|
||||
if (IsFatalResourceError(e)) throw;
|
||||
// Angles under which nothing indexes have scored nothing, which is what an empty
|
||||
// result reads as.
|
||||
logger.Info("Rotation scale {:.5f}: the probe pass did not complete ({})", k, e.what());
|
||||
}
|
||||
indexing_probe_only_ = postrefine_probe_ = postrefine_probe_only_ = false;
|
||||
++arm_passes;
|
||||
return r;
|
||||
};
|
||||
const auto index_at = [&](float k) { return probe(k, true).validation_evidence; };
|
||||
const auto refit_at = [&](float k) -> std::optional<double> {
|
||||
const ProcessResult r = probe(k, false);
|
||||
if (!r.post_refine)
|
||||
return std::nullopt;
|
||||
logger.Info("Rotation scale {:.5f}: the post-refinement there fits {:.5f} on top of it, "
|
||||
"held-out residual {:.3e}", k, r.post_refine->rotation_scale,
|
||||
r.post_refine->held_out_before);
|
||||
return r.post_refine->rotation_scale;
|
||||
};
|
||||
constexpr int MAX_ROTATION_SCALE_ROUNDS = 8;
|
||||
const RotationScaleWalk walk = WalkRotationScale(pass1.post_refine->rotation_scale, index_at,
|
||||
refit_at, MAX_ROTATION_SCALE_ROUNDS);
|
||||
experiment_ = before_walk;
|
||||
prepass_detector_geometry_ = geometry_before_walk;
|
||||
prepass_mosaicity_ = mosaicity_before_walk;
|
||||
beam_center_searched_ = searched_before_walk;
|
||||
prepass_rotation_scale_.reset();
|
||||
if (!walk.trail.empty()) {
|
||||
rotation_scale_walk = fmt::format(
|
||||
"goniometer rotation scale walked from the pass-1 fit {:.5f}, validation spots on "
|
||||
"the lattice beyond chance by scale: {} - {}", pass1.post_refine->rotation_scale,
|
||||
walk.trail, walk.scale == 1.0f ? "the stored angles stand"
|
||||
: fmt::format("{:.5f} adopted", walk.scale));
|
||||
logger.Info("Two-pass: {}", rotation_scale_walk);
|
||||
}
|
||||
if (walk.scale != 1.0f) {
|
||||
prepass_rotation_scale_ = walk.scale;
|
||||
pass1.post_refine->rotation_scale = walk.scale;
|
||||
pass1.post_refine->rotation_scale_suspect = true;
|
||||
logger.Warning("Goniometer rotation scale {:.5f} ({:+.2f} %): the stage turned {} than the "
|
||||
"angles stored in the file, which are the COMMANDED values. The second pass "
|
||||
"integrates at the corrected angles; the fault is in the hardware and should "
|
||||
"be fixed there", walk.scale, 100.0 * (walk.scale - 1.0),
|
||||
walk.scale > 1.0f ? "further" : "less far");
|
||||
}
|
||||
}
|
||||
|
||||
// ... and apply the adopted scale, on the same measure-then-re-integrate footing as the geometry.
|
||||
if (prepass_rotation_scale_ && gonio_snapshot) {
|
||||
experiment_.Goniometer(ScaleRotation(*gonio_snapshot, *prepass_rotation_scale_));
|
||||
experiment_.Goniometer(ScaleRotation(*experiment_.GetGoniometer(), *prepass_rotation_scale_));
|
||||
// The pre-pass mosaicity is a width in degrees fitted against the angles the second pass has
|
||||
// just stopped using, and the override can only ever raise the second pass's own estimate (it
|
||||
// takes the larger of the two). Carrying it over would hold the second pass at the rocking
|
||||
@@ -2197,10 +2316,6 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
|
||||
.PoniRot1_rad(g[3]).PoniRot2_rad(g[4]);
|
||||
};
|
||||
|
||||
// The lattice arms below run passes that 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;
|
||||
|
||||
// The metric symmetry, asked of the spots rather than of the tolerance that admitted it. The
|
||||
// Bravais class is chosen by a walk that reads two axes as EQUAL when they agree to a fixed
|
||||
// relative tolerance, and the class carrying that equality is then imposed everywhere below:
|
||||
@@ -2548,6 +2663,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
|
||||
if (const auto current = experiment_.GetGoniometer())
|
||||
restored.Axis(current->GetAxis());
|
||||
experiment_.Goniometer(restored);
|
||||
prepass_rotation_scale_.reset();
|
||||
pass1.post_refine->rotation_scale_suspect = false;
|
||||
}
|
||||
config_.output_prefix = base_prefix;
|
||||
// This pass is the answer whatever it measures - the guard has had its one chance -
|
||||
@@ -2573,6 +2690,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
|
||||
: "the post-refinement committed no geometry change, so this pass reproduces the first";
|
||||
if (!lattice_arm.empty())
|
||||
pass2.pass_decision = lattice_arm + "; " + pass2.pass_decision;
|
||||
if (!rotation_scale_walk.empty())
|
||||
pass2.pass_decision = rotation_scale_walk + "; " + pass2.pass_decision;
|
||||
pass2.geometry_not_converged = geometry_not_converged;
|
||||
return pass2;
|
||||
}
|
||||
@@ -3234,11 +3353,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
|
||||
// The pooled evidence for one candidate lattice: the validation frames' non-ice spots, how many
|
||||
// of them lie on the candidate, and how many a wrong spindle angle still puts there.
|
||||
struct PooledEvidence {
|
||||
int64_t spots = 0;
|
||||
int64_t on_lattice = 0;
|
||||
int64_t by_chance = 0;
|
||||
};
|
||||
using PooledEvidence = ValidationSpotEvidence;
|
||||
|
||||
// Displacements are a fixed fraction of the sweep, so the null is a function of the data alone
|
||||
// and the same file gives the same verdict every time. Sevenths: no crystallographic rotation
|
||||
@@ -3274,19 +3389,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
return static_cast<double>(e.on_lattice - e.by_chance) > SPOT_BUDGET_SIGNIFICANCE_Z * sigma;
|
||||
};
|
||||
|
||||
// The same quantity for COMPARING two lattices: the share of the pooled spots each one puts on
|
||||
// itself over and above what a wrong spindle angle puts there. Each is scored against its own
|
||||
// null, so a denser lattice is not credited for the spots it catches by accident - which is
|
||||
// what makes the comparison fair between a cell and its axis harmonic.
|
||||
auto excess = [](const PooledEvidence &e) {
|
||||
return static_cast<double>(e.on_lattice - e.by_chance)
|
||||
/ static_cast<double>(std::max<int64_t>(1, e.spots));
|
||||
};
|
||||
auto rate_var = [](const PooledEvidence &e) {
|
||||
const double n = static_cast<double>(std::max<int64_t>(1, e.spots));
|
||||
const double p = static_cast<double>(e.on_lattice) / n;
|
||||
return p * (1.0 - p) / n;
|
||||
};
|
||||
// Two lattices are COMPARED on the same quantity, ValidationEvidencePrefers: the share of the
|
||||
// pooled spots each one puts on itself over and above what a wrong spindle angle puts there.
|
||||
// Each is scored against its own null, which is what makes the comparison fair between a cell
|
||||
// and its axis harmonic.
|
||||
|
||||
// How deep into an image's intensity-ordered spot list this lattice is still being seen - the
|
||||
// measured spot budget. Same frames, same per-image path as count_indexed above, but scoring the
|
||||
@@ -4033,10 +4139,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
restore_beam_center(header_x, header_y);
|
||||
if (best.result.has_value())
|
||||
at_header = pooled_evidence(*indexer, *best.result);
|
||||
const double margin = excess(at_measured) - excess(at_header);
|
||||
const double margin_sigma = std::sqrt(rate_var(at_measured) + rate_var(at_header));
|
||||
if (alt.result.has_value() && beats_chance(at_measured)
|
||||
&& margin > SPOT_BUDGET_SIGNIFICANCE_Z * margin_sigma) {
|
||||
&& ValidationEvidencePrefers(at_header, at_measured)) {
|
||||
try_beam_center(measured_x, measured_y);
|
||||
logger.Warning("Beam centre from the background: neither centre indexes a "
|
||||
"validation frame on its own, but the pooled spots put "
|
||||
@@ -4179,11 +4283,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
const PooledEvidence at_header = pooled_evidence(*indexer, *best.result);
|
||||
try_beam_center(measured_x, measured_y);
|
||||
const PooledEvidence at_measured = pooled_evidence(*indexer, *alt.result);
|
||||
const double margin = excess(at_measured) - excess(at_header);
|
||||
const double margin_sigma = std::sqrt(rate_var(at_measured)
|
||||
+ rate_var(at_header));
|
||||
adopted_measured = beats_chance(at_measured)
|
||||
&& margin > SPOT_BUDGET_SIGNIFICANCE_Z * margin_sigma;
|
||||
&& ValidationEvidencePrefers(at_header, at_measured);
|
||||
if (adopted_measured) {
|
||||
logger.Warning("Beam centre check: the larger cell is the measured "
|
||||
"centre's, and it is the one the spots are on - "
|
||||
@@ -4639,6 +4740,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
gemmi::crystal_system_str(prepass_result_->search_result.system),
|
||||
pc.a, pc.b, pc.c, pc.alpha, pc.beta, pc.gamma);
|
||||
best.result = *prepass_result_;
|
||||
// Pass 1's result carries pass 1's goniometer; forcing it whole would put the
|
||||
// uncorrected angles back (the same repair as the supercell re-run in RunAllPasses).
|
||||
if (prepass_rotation_scale_ && best.result->axis)
|
||||
best.result->axis = ScaleRotation(*best.result->axis, *prepass_rotation_scale_);
|
||||
// Pass 1's lattice was FITTED at pass 1's detector distance, and this pass
|
||||
// integrates at the post-refined one. Re-scoring it here is not the same as
|
||||
// re-fitting it: a real-space cell is measured against the distance the spots were
|
||||
@@ -4685,6 +4790,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
/ static_cast<double>(std::max<int64_t>(1, evidence.spots));
|
||||
const double pooled_chance = static_cast<double>(evidence.by_chance)
|
||||
/ static_cast<double>(std::max<int64_t>(1, evidence.spots));
|
||||
result.validation_evidence = evidence;
|
||||
// A pass run only to score a lattice has its score, whatever it is: a lattice that does not
|
||||
// beat chance is a result for the comparison that asked, not a reason to stop the run.
|
||||
if (indexing_probe_only_) {
|
||||
logger.Info("Indexing probe: scheme '{}', {}/{} validation frames, {}/{} validation spots "
|
||||
"= {:.1f}% against {:.1f}% at a wrong spindle angle", best.name, best.score,
|
||||
static_cast<int>(validation.size()), evidence.on_lattice, evidence.spots,
|
||||
100.0 * pooled, 100.0 * pooled_chance);
|
||||
return result;
|
||||
}
|
||||
if (!beats_chance(evidence)) {
|
||||
{
|
||||
// Name the cell and Bravais class that was rejected. The commonest cause is a metric
|
||||
@@ -4821,7 +4936,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
"pass 1 ({}-centred, {:.0f} A^3) - integrating with pass-1's lattice instead",
|
||||
best.result->search_result.centering, v2,
|
||||
prepass_result_->search_result.centering, v1);
|
||||
indexer->ForceRotationIndexerResult(*prepass_result_);
|
||||
RotationIndexerResult forced = *prepass_result_;
|
||||
if (prepass_rotation_scale_ && forced.axis)
|
||||
forced.axis = ScaleRotation(*forced.axis, *prepass_rotation_scale_);
|
||||
indexer->ForceRotationIndexerResult(forced);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -5529,19 +5647,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
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
|
||||
// exactly the ones whose cell and detector steps do NOT pass, because the angle error
|
||||
// is what their residual is made of. A calibration fault is rare (36 of 37 rotation
|
||||
// datasets sit at 1.0000) and applying a 1 % angle correction to a healthy dataset
|
||||
// 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.
|
||||
// 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<float>(pr.rotation_scale);
|
||||
// The goniometer rotation scale in pr is only a fit: RunAllPasses walks it to its
|
||||
// fixed point and decides it on the validation frames of the whole sweep
|
||||
// (WalkRotationScale). It does not ride on pr.ok - the crystals that have a stage
|
||||
// fault are exactly the ones whose cell and detector steps do NOT pass, because the
|
||||
// angle error is what their residual is made of.
|
||||
}
|
||||
}
|
||||
};
|
||||
@@ -5559,8 +5669,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
if (postrefine_probe_only_) {
|
||||
result.processing_time_s = std::chrono::duration<double>(
|
||||
std::chrono::steady_clock::now() - start_time).count();
|
||||
logger.Info("Lattice-arm probe: {} images integrated and post-refined in {:.2f} s "
|
||||
"(the arm reads the held-out residual only, so this pass does not merge)",
|
||||
logger.Info("Probe pass: {} images integrated and post-refined in {:.2f} s "
|
||||
"(the comparison reads the post-refinement only, so this pass does not merge)",
|
||||
result.images_processed, result.processing_time_s);
|
||||
return result;
|
||||
}
|
||||
|
||||
+46
-2
@@ -5,6 +5,7 @@
|
||||
|
||||
#include <array>
|
||||
#include <atomic>
|
||||
#include <functional>
|
||||
#include <optional>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
@@ -186,6 +187,42 @@ struct ProcessConfig {
|
||||
bool finalist_ledger = false; // --finalist-ledger; report-only symmetry evidence table
|
||||
};
|
||||
|
||||
// A rotation lattice's score on the validation frames spread over the whole sweep: their spots (off
|
||||
// the ice rings unless those are indexed), how many lie on the lattice, and how many still do at a
|
||||
// wrong spindle angle - the median over displaced angles, which is what chance and the per-frame
|
||||
// orientation polish give for free.
|
||||
struct ValidationSpotEvidence {
|
||||
int64_t spots = 0;
|
||||
int64_t on_lattice = 0;
|
||||
int64_t by_chance = 0;
|
||||
};
|
||||
|
||||
// Whether `candidate` puts a larger share of its validation spots on its lattice, over and above what
|
||||
// a wrong spindle angle puts there, than `current` does - by more than the binomial noise of the two
|
||||
// shares (SPOT_BUDGET_SIGNIFICANCE_Z). Each is measured against its own null, so a denser lattice is
|
||||
// not credited for the spots it catches by accident.
|
||||
bool ValidationEvidencePrefers(const ValidationSpotEvidence ¤t, const ValidationSpotEvidence &candidate);
|
||||
|
||||
// The goniometer rotation scale - the factor by which the stage turned relative to the angles stored
|
||||
// in the file - walked to the fit's fixed point and decided on the whole sweep. The fit only sees the
|
||||
// frames the angles it was measured at still track, and a stage at the wrong rate is exactly what
|
||||
// stops them tracking the rest, so one fit reads only part of the error. So: index the lattice under
|
||||
// the fitted scale and under the angles in hand; where the fitted one scores better on the validation
|
||||
// frames (ValidationEvidencePrefers), adopt it, fit again there and repeat. The stored angles (k = 1)
|
||||
// are the first hypothesis, and stand unless the evidence moves the run off them.
|
||||
// index_at(k): the validation evidence of the lattice indexed with the stored angles scaled by k
|
||||
// refit_at(k): the scale the post-refinement fits on reflections integrated at k, relative to k;
|
||||
// empty where it fitted none
|
||||
struct RotationScaleWalk {
|
||||
float scale = 1.0f; // adopted; 1 = the stored angles stand
|
||||
ValidationSpotEvidence evidence; // at the adopted scale
|
||||
std::string trail; // every scale tried, with its validation score
|
||||
};
|
||||
RotationScaleWalk WalkRotationScale(double first_fit,
|
||||
const std::function<ValidationSpotEvidence(float)> &index_at,
|
||||
const std::function<std::optional<double>(float)> &refit_at,
|
||||
int max_rounds);
|
||||
|
||||
struct ProcessResult {
|
||||
bool cancelled = false;
|
||||
uint64_t images_processed = 0;
|
||||
@@ -233,6 +270,8 @@ struct ProcessResult {
|
||||
// cell - so two cells can only be compared for volume once this has brought both to primitive.
|
||||
std::optional<char> consensus_centering;
|
||||
bool rotation_lattice_found = false;
|
||||
// The rotation lattice's score on the validation frames; all zero where no lattice was scored.
|
||||
ValidationSpotEvidence validation_evidence;
|
||||
// The metric symmetry the rotation lattice was classified in - the class GeometryRefiner holds the
|
||||
// cell to and the space-group search enumerates under. Empty when no rotation lattice was found.
|
||||
std::optional<LatticeMessage> rotation_lattice_type;
|
||||
@@ -484,8 +523,9 @@ class Rugnux {
|
||||
// describes neither geometry. The two travel together or not at all.
|
||||
std::optional<std::array<float, 5>> prepass_detector_geometry_;
|
||||
|
||||
// Two-pass geometry pre-pass: the goniometer rotation scale the post-refine fitted and flagged as a
|
||||
// stage fault, applied by Run() to the second pass's goniometer. Empty when the fit found nothing.
|
||||
// The goniometer rotation scale the run adopted (WalkRotationScale, in RunAllPasses), applied to the
|
||||
// second pass's goniometer - and, while the walk runs, the scale its current probe pass is at.
|
||||
// Empty where the stored angles stand.
|
||||
std::optional<float> prepass_rotation_scale_;
|
||||
|
||||
// Two-pass geometry pre-pass: pass-1's FULL indexing result (correct lattice + orientation + refined
|
||||
@@ -591,6 +631,10 @@ class Rugnux {
|
||||
// Everything below that point - the merges, the space-group search, the correction surfaces, the
|
||||
// reports - is work no comparison looks at, and on a long sweep it is four fifths of the pass.
|
||||
bool postrefine_probe_only_ = false;
|
||||
// Whether a pass exists only to index: it stops as soon as its first-pass indexing has scored the
|
||||
// lattice on the validation frames (ProcessResult::validation_evidence), before a single image is
|
||||
// integrated. The rotation-scale walk (see RunAllPasses) compares angle models on that score.
|
||||
bool indexing_probe_only_ = false;
|
||||
// The file's detector distance, from before the first pass of the rotation two-pass: a pass whose
|
||||
// distance is still this one asks the post-refinement to test "the header distance is right"
|
||||
// (PostRefineSettings::distance_at_header); a pass that has walked off it does not.
|
||||
|
||||
@@ -59,6 +59,7 @@ ADD_EXECUTABLE(jfjoch_test
|
||||
ResultReportTest.cpp
|
||||
DiagnosticOutputTest.cpp
|
||||
PostRefineTest.cpp
|
||||
RotationScaleWalkTest.cpp
|
||||
RugnuxLargeTest.cpp
|
||||
TestData.h
|
||||
MovingAverageTest.cpp
|
||||
|
||||
@@ -0,0 +1,94 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <cmath>
|
||||
#include <catch2/catch_all.hpp>
|
||||
|
||||
#include "../rugnux/Rugnux.h"
|
||||
|
||||
namespace {
|
||||
// A synthetic sweep whose stage turned `true_scale` times the stored angles. Scored at a scale k,
|
||||
// the validation spots stay on the lattice as long as the angles track the rotation, and the share
|
||||
// that does falls off with the relative rate error; a wrong spindle angle keeps half a percent.
|
||||
struct SyntheticSweep {
|
||||
double true_scale;
|
||||
int64_t spots = 35000;
|
||||
int index_calls = 0;
|
||||
int refit_calls = 0;
|
||||
|
||||
ValidationSpotEvidence IndexAt(float k) {
|
||||
++index_calls;
|
||||
const double error = std::fabs(true_scale / k - 1.0);
|
||||
const double on = 0.9 * std::max(0.0, 1.0 - 30.0 * error);
|
||||
return ValidationSpotEvidence{spots, std::llround(on * spots), std::llround(0.005 * spots)};
|
||||
}
|
||||
|
||||
// The post-refinement at k, relative to k. It reads only 70 % of the error that is left: it
|
||||
// sees only the frames the angles at k still track.
|
||||
std::optional<double> RefitAt(float k) {
|
||||
++refit_calls;
|
||||
return 1.0 + 0.7 * (true_scale / k - 1.0);
|
||||
}
|
||||
|
||||
RotationScaleWalk Walk(double first_fit) {
|
||||
return WalkRotationScale(first_fit, [this](float k) { return IndexAt(k); },
|
||||
[this](float k) { return RefitAt(k); }, 8);
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
TEST_CASE("ValidationEvidencePrefers", "[RotationScale]") {
|
||||
const ValidationSpotEvidence base{10000, 3000, 50};
|
||||
// 1 % more of the spots beyond chance is under the noise of two 30 % shares over 10000 spots
|
||||
// (sqrt(2 * 0.3 * 0.7 / 10000) = 0.65 %, times 3.29); 5 % is well over it.
|
||||
CHECK_FALSE(ValidationEvidencePrefers(base, {10000, 3100, 50}));
|
||||
CHECK(ValidationEvidencePrefers(base, {10000, 3500, 50}));
|
||||
// A candidate is judged against its own null: more spots on the lattice bought by a null that
|
||||
// rose just as much is no gain.
|
||||
CHECK_FALSE(ValidationEvidencePrefers(base, {10000, 3500, 550}));
|
||||
// Never against itself, and nothing that scored nothing wins.
|
||||
CHECK_FALSE(ValidationEvidencePrefers(base, base));
|
||||
CHECK_FALSE(ValidationEvidencePrefers(base, {}));
|
||||
CHECK(ValidationEvidencePrefers({}, base));
|
||||
}
|
||||
|
||||
TEST_CASE("WalkRotationScale_ReachesTheFixedPoint", "[RotationScale]") {
|
||||
// A stage 3 % slow. The first fit reads 70 % of that; each refit at the adopted scale reads 70 %
|
||||
// of what is left, and the walk goes on as long as the validation frames prefer the new scale.
|
||||
SyntheticSweep sweep{0.97};
|
||||
const auto walk = sweep.Walk(sweep.RefitAt(1.0f).value());
|
||||
CHECK(walk.scale == Catch::Approx(0.97).margin(0.001));
|
||||
CHECK(walk.scale != 1.0f);
|
||||
CHECK(walk.evidence.on_lattice > sweep.IndexAt(1.0f).on_lattice);
|
||||
CHECK(sweep.refit_calls > 2);
|
||||
CHECK_FALSE(walk.trail.empty());
|
||||
}
|
||||
|
||||
TEST_CASE("WalkRotationScale_StoredAnglesStand", "[RotationScale]") {
|
||||
SECTION("A healthy stage: a fit off by noise scores no better than the stored angles") {
|
||||
SyntheticSweep sweep{1.0};
|
||||
const auto walk = sweep.Walk(1.0002);
|
||||
CHECK(walk.scale == 1.0f);
|
||||
CHECK(sweep.index_calls == 2); // the stored angles and the fit, nothing more
|
||||
CHECK(sweep.refit_calls == 0);
|
||||
}
|
||||
SECTION("A real but small error the spots cannot resolve beyond their noise") {
|
||||
SyntheticSweep sweep{0.999};
|
||||
sweep.spots = 400;
|
||||
const auto walk = sweep.Walk(0.9993);
|
||||
CHECK(walk.scale == 1.0f);
|
||||
}
|
||||
SECTION("A fit that tracks something other than the rotation scores worse, and is refused") {
|
||||
SyntheticSweep sweep{1.0};
|
||||
const auto walk = sweep.Walk(0.98);
|
||||
CHECK(walk.scale == 1.0f);
|
||||
CHECK(sweep.refit_calls == 0);
|
||||
}
|
||||
SECTION("A fit of exactly one asks for no probe at all") {
|
||||
SyntheticSweep sweep{1.0};
|
||||
const auto walk = sweep.Walk(1.0);
|
||||
CHECK(walk.scale == 1.0f);
|
||||
CHECK(walk.trail.empty());
|
||||
CHECK(sweep.index_calls == 0);
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user