leonarski_fandClaude Opus 5
e61454a27b
The twin-law report block dereferenced the adopted space group without the guard every other use of it carries; a de-novo run whose search refused a point group but left no adopted group would crash there, after the merge. The joint post-refinement's distance-correlation array was read on every run that holds the header distance, but only written when the solve is usable, so an unusable first solve printed four indeterminate values as correlations. The merge-degradation gate ran neither of its two branches when the data confirm more than one operator and the best of them is itself close to the random-pairing end: the contrast has no scale to read on such a merge, but the candidate was then promoted on no merge-degradation evidence at all. It now falls back to the noise-floor ratio, as it does when only one operator is confirmed. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT
1174 lines
76 KiB
C++
1174 lines
76 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "../../common/ParallelFor.h"
|
|
#include "PostRefine.h"
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <limits>
|
|
#include <map>
|
|
#include <memory>
|
|
#include <numeric>
|
|
#include <string>
|
|
|
|
#include "../../common/JFJochMath.h" // PI
|
|
#include "XtalResidual.h" // the positional detector<->reciprocal residual, and the cell it is parameterised by
|
|
#include "LatticeReduction.h"
|
|
#include "ceres/ceres.h"
|
|
#include "ceres/rotation.h"
|
|
|
|
namespace {
|
|
|
|
// How far the joint fit may move the beam centre away from a value something else already
|
|
// believes. The bound is not about how large a real correction can be - it is the backstop 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. The absolute size of the move is what
|
|
// separates that from a genuine header correction.
|
|
//
|
|
// It is measured from whichever of the nominal centre and the run's own MEASUREMENT of the centre
|
|
// (PostRefineSettings::measured_beam_px) is nearer. Anchored on the nominal centre alone it caps the
|
|
// correction at exactly the header's own error, which is the opposite of its job: the header is most
|
|
// worth correcting when it is most wrong. Measured on a rotation crystal whose header centre is 22.1 px
|
|
// out - a value the same run's beam-centre check had already placed to +-0.30 px and then discarded -
|
|
// the old bound rejected a fit that improved the held-out positional residual nine-fold.
|
|
constexpr double BEAM_BOUND_PXL = 15.0;
|
|
|
|
// How far the joint fit may move a cell ANGLE. Only monoclinic (beta) and triclinic leave one free,
|
|
// and those are the two systems whose conditioning is weakest, so the angles need the same absolute
|
|
// backstop the lengths and the distance already have rather than being left to the solver's box. The
|
|
// largest genuine angle move anywhere on the corpus is 0.35 deg, on a monoclinic beta, measured
|
|
// against the same run's own second pass; one degree is about twice that, and still well inside the
|
|
// +-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 - PostRefinePartial in the header, under its
|
|
// working name (see there for why it is kept this narrow).
|
|
using Partial = PostRefinePartial;
|
|
|
|
// A rocking event: one reflection's intensity-weighted centroid over the frames it spans.
|
|
struct Event {
|
|
double phi_obs; // rad, intensity-weighted rocking centroid
|
|
double weight; // sqrt(sum I / sum sigma)
|
|
int h, k, l;
|
|
};
|
|
|
|
// Distance-INDEPENDENT Ewald excitation residual with the whole crystal free. The observed rocking
|
|
// centroid phi_obs says where the reflection actually crossed the Ewald sphere, which fixes the
|
|
// ABSOLUTE size of the reciprocal lattice with the detector never entering - that is what breaks the
|
|
// distance <-> cell-scale degeneracy the spot positions alone leave open, and it is why the two can be
|
|
// fitted together. The parameter blocks are XtalResidual's crystal half (rotation axis, orientation,
|
|
// cell lengths, cell angles), so one problem can share them between the two residuals.
|
|
// On the Ewald sphere <=> |p|^2 + 2 p_z/lambda == 0.
|
|
struct JointExcitationResidual {
|
|
JointExcitationResidual(double lambda, double angle_rad, double weight,
|
|
int h, int k, int l, gemmi::CrystalSystem symmetry)
|
|
: inv_lambda(1.0 / lambda), angle_rad(angle_rad), weight(weight),
|
|
h(h), k(k), l(l), symmetry(symmetry) {}
|
|
template<typename T>
|
|
bool operator()(const T *const axis, const T *const p0, const T *const p1, const T *const p2,
|
|
T *residual) const {
|
|
Eigen::Matrix<T, 3, 1> bxc, cxa, axb;
|
|
T invV;
|
|
XtalResidual::ReciprocalBasis(p1, p2, symmetry, bxc, cxa, axb, invV);
|
|
const Eigen::Matrix<T, 3, 1> unrot = (bxc * T(h) + cxa * T(k) + axb * T(l)) * invV;
|
|
const T recip_unrot[3] = {unrot[0], unrot[1], unrot[2]};
|
|
T p_ref[3];
|
|
const AngleAxisRotator<T> rot_p0(p0);
|
|
rot_p0.Rotate(recip_unrot, p_ref);
|
|
const T aa[3] = {T(-angle_rad) * axis[0], T(-angle_rad) * axis[1], T(-angle_rad) * axis[2]};
|
|
T p_lab[3];
|
|
ceres::AngleAxisRotatePoint(aa, p_ref, p_lab);
|
|
const T zeta = p_lab[0] * p_lab[0] + p_lab[1] * p_lab[1] + p_lab[2] * p_lab[2]
|
|
+ T(2.0) * p_lab[2] * T(inv_lambda);
|
|
residual[0] = T(weight) * zeta * T(0.5) / T(inv_lambda);
|
|
return true;
|
|
}
|
|
const double inv_lambda, angle_rad, weight;
|
|
const double h, k, l;
|
|
const gemmi::CrystalSystem symmetry;
|
|
};
|
|
|
|
// GONIOMETER ROTATION SCALE k: the same Ewald excitation residual, but with the crystal and the axis
|
|
// DIRECTION already committed by the joint fit, so the single free quantity is how far the stage actually
|
|
// turned per unit of commanded angle. Two differences from that fit matter:
|
|
// * the angle is measured from the CENTRE of the sweep, not from the goniometer's zero. The reference
|
|
// orientation is the one rotation indexing fitted against the commanded angles, so it has already
|
|
// absorbed the MEAN angle error; only the part that varies across the sweep is left to fit. Scaling the
|
|
// absolute angle instead - which is what reading k off the length of the fitted axis vector does - asks
|
|
// the fit to also produce a constant offset it has no parameter for, and the least-squares compromise
|
|
// shrinks k towards 1 by var(phi) / (var(phi) + phi_centre^2): exactly a factor of four for the common
|
|
// case of a sweep starting at zero.
|
|
// * e_mid is the committed reciprocal vector already turned to the sweep centre, so nothing but k is
|
|
// free.
|
|
struct RotationScaleResidual {
|
|
RotationScaleResidual(double lambda, double dangle_rad, const double u[3], const double e_mid[3])
|
|
: inv_lambda(1.0 / lambda), dangle_rad(dangle_rad),
|
|
ux(u[0]), uy(u[1]), uz(u[2]), ex(e_mid[0]), ey(e_mid[1]), ez(e_mid[2]) {}
|
|
template<typename T>
|
|
bool operator()(const T *const k, T *residual) const {
|
|
const T a = T(-dangle_rad) * k[0];
|
|
const T aa[3] = {a * T(ux), a * T(uy), a * T(uz)};
|
|
const T p_ref[3] = {T(ex), T(ey), T(ez)};
|
|
T p_lab[3];
|
|
ceres::AngleAxisRotatePoint(aa, p_ref, p_lab);
|
|
const T zeta = p_lab[0] * p_lab[0] + p_lab[1] * p_lab[1] + p_lab[2] * p_lab[2]
|
|
+ T(2.0) * p_lab[2] * T(inv_lambda);
|
|
residual[0] = zeta * T(0.5) / T(inv_lambda);
|
|
return true;
|
|
}
|
|
const double inv_lambda, dangle_rad, ux, uy, uz, ex, ey, ez;
|
|
};
|
|
|
|
} // namespace
|
|
|
|
PostRefineObservations GatherPostRefineObservations(std::vector<IntegrationOutcome> &outcomes,
|
|
size_t nthreads, bool release_reflections,
|
|
Logger &logger) {
|
|
PostRefineObservations out;
|
|
try {
|
|
nthreads = std::max<size_t>(1, nthreads);
|
|
|
|
// Count first, then fill. Growing one vector by push_back over tens of millions of
|
|
// reflections copies the whole thing every time it doubles - several gigabytes of pure
|
|
// copying - and the counts are cheap to take. The h range comes out of the same sweep: the
|
|
// sort below buckets by h and needs to know how many buckets that is, and this pass already
|
|
// reads every reflection.
|
|
const int n_out = static_cast<int>(outcomes.size());
|
|
std::vector<size_t> counts(n_out, 0);
|
|
std::vector<int> h_lo_of(n_out), h_hi_of(n_out);
|
|
ParallelChunks(n_out, nthreads, [&](int lo, int hi) {
|
|
for (int o = lo; o < hi; o++) {
|
|
size_t keep = 0;
|
|
int lmin = std::numeric_limits<int>::max(), lmax = std::numeric_limits<int>::min();
|
|
for (const auto &r : outcomes[o].reflections)
|
|
if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) {
|
|
keep++;
|
|
lmin = std::min(lmin, r.h);
|
|
lmax = std::max(lmax, r.h);
|
|
}
|
|
counts[o] = keep;
|
|
h_lo_of[o] = lmin;
|
|
h_hi_of[o] = lmax;
|
|
}
|
|
});
|
|
int h_lo = std::numeric_limits<int>::max(), h_hi = std::numeric_limits<int>::min();
|
|
size_t n_pts = 0;
|
|
for (int o = 0; o < n_out; o++) {
|
|
n_pts += counts[o];
|
|
h_lo = std::min(h_lo, h_lo_of[o]);
|
|
h_hi = std::max(h_hi, h_hi_of[o]);
|
|
}
|
|
const int H = (h_lo <= h_hi) ? (h_hi - h_lo + 1) : 1;
|
|
|
|
// The bucket histogram, in a pass of its own: the fill below scatters each partial straight
|
|
// into its bucket, so the gathered-order copy of the whole array that used to be scattered
|
|
// from - gigabytes on a large cell, alive beside the bucketed one - no longer exists at all.
|
|
// Its chunks are the outcome chunks ParallelChunks makes, so the scatter below has to be
|
|
// split the same way.
|
|
const int nt = static_cast<int>(std::clamp<size_t>(nthreads, 1, std::max(1, n_out)));
|
|
const int chunk = (n_out + nt - 1) / nt;
|
|
std::vector<std::vector<int32_t>> hist(nt, std::vector<int32_t>(H, 0));
|
|
ParallelChunks(n_out, nthreads, [&](int lo, int hi) {
|
|
std::vector<int32_t> &h_count = hist[lo / chunk];
|
|
for (int o = lo; o < hi; o++)
|
|
for (const auto &r : outcomes[o].reflections)
|
|
if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f)
|
|
h_count[r.h - h_lo]++;
|
|
});
|
|
|
|
// Where each chunk writes each bucket: buckets in order, and within a bucket the chunks in
|
|
// index order, so the scatter below is stable and lands in original (per-image) order inside
|
|
// a bucket - which is what makes the tie order below reproducible at every thread count.
|
|
out.bstart.assign(H + 1, 0);
|
|
{
|
|
int32_t acc = 0;
|
|
for (int b = 0; b < H; ++b) {
|
|
out.bstart[b] = acc;
|
|
for (int t = 0; t < nt; ++t) { const int32_t c = hist[t][b]; hist[t][b] = acc; acc += c; }
|
|
}
|
|
out.bstart[H] = acc;
|
|
}
|
|
|
|
// A vector of n partials VALUE-initialises them: on a large cell that is gigabytes of zeroing
|
|
// on one thread, and it is that one thread which first touches every page - which on a
|
|
// multi-socket machine leaves the whole array on its node, so every pass that follows runs at
|
|
// one node's memory bandwidth. new[] leaves the partials untouched, so the parallel scatter is
|
|
// the first touch and each page lands on the node of the thread that filled it.
|
|
std::unique_ptr<Partial[]> pts(new Partial[n_pts]);
|
|
ParallelChunks(n_out, nthreads, [&](int lo, int hi) {
|
|
std::vector<int32_t> fill = hist[lo / chunk];
|
|
for (int o = lo; o < hi; o++) {
|
|
for (const auto &r : outcomes[o].reflections) {
|
|
if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f) continue;
|
|
const float ox = std::isfinite(r.observed_x) ? r.observed_x : NAN;
|
|
const float oy = std::isfinite(r.observed_y) ? r.observed_y : NAN;
|
|
pts[fill[r.h - h_lo]++] = Partial{static_cast<int16_t>(r.h), static_cast<int16_t>(r.k),
|
|
static_cast<int16_t>(r.l), r.image_number,
|
|
r.I, r.sigma, ox, oy};
|
|
}
|
|
// The caller said it will never read this image's reflections again, so hand the
|
|
// vector back the moment it is consumed: the payload shrinks image by image as it
|
|
// moves into the (much narrower) partials instead of doubling up in full.
|
|
if (release_reflections)
|
|
ReflectionVector().swap(outcomes[o].reflections);
|
|
}
|
|
});
|
|
logger.Info("Post-refine: {} partials gathered", n_pts);
|
|
|
|
// Sort each bucket, on its own thread and the largest first (the tail of this is one bucket,
|
|
// so it should be the big one). h is the leading key, so the sorted array is the buckets laid
|
|
// end to end. Sorting the whole thing in one pass moved every partial through every level of
|
|
// a comparison sort, on one thread, over tens of millions of reflections.
|
|
//
|
|
// Not a total order: two partials of one reflection on one image still tie, as they did
|
|
// before this was bucketed. What makes the result reproducible is the scatter above rather
|
|
// than the comparator - each bucket is laid out chunk by chunk and a chunk covers the images
|
|
// in index order, so a bucket reaches std::sort in global per-image order whatever the
|
|
// thread count. Ties therefore resolve the same way on every run and at every -N; they are
|
|
// simply not resolved by rank.
|
|
const auto part_less = [](const Partial &a, const Partial &b) {
|
|
if (a.h != b.h) return a.h < b.h;
|
|
if (a.k != b.k) return a.k < b.k;
|
|
if (a.l != b.l) return a.l < b.l;
|
|
return a.img < b.img;
|
|
};
|
|
out.order.resize(H);
|
|
std::iota(out.order.begin(), out.order.end(), 0);
|
|
std::sort(out.order.begin(), out.order.end(), [&](int a, int b) {
|
|
return (out.bstart[a + 1] - out.bstart[a]) > (out.bstart[b + 1] - out.bstart[b]);
|
|
});
|
|
ParallelFor(H, nthreads, [&](int oi) {
|
|
const int b = out.order[oi];
|
|
std::sort(pts.get() + out.bstart[b], pts.get() + out.bstart[b + 1], part_less);
|
|
});
|
|
out.pts = std::move(pts);
|
|
out.n_pts = n_pts;
|
|
} catch (...) {
|
|
out = PostRefineObservations{}; // empty; the fit refuses it (see the header)
|
|
}
|
|
return out;
|
|
}
|
|
|
|
PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations,
|
|
const GoniometerAxis &axis,
|
|
const DiffractionGeometry &nominal_geom,
|
|
const CrystalLattice &reference_latt,
|
|
const PostRefineSettings &settings,
|
|
Logger &logger) {
|
|
PostRefineResult result;
|
|
result.geom = nominal_geom;
|
|
result.cell = reference_latt.GetUnitCell();
|
|
result.distance_before_mm = nominal_geom.GetDetectorDistance_mm();
|
|
result.distance_after_mm = nominal_geom.GetDetectorDistance_mm();
|
|
try {
|
|
const double wedge_half = axis.GetWedge_deg() / 2.0;
|
|
const double lambda = nominal_geom.GetWavelength_A();
|
|
const Coord ax = axis.GetAxis();
|
|
|
|
// The goniometer angle is a function of the image number alone, so it is recomputed where it
|
|
// is used rather than carried through the gathered array: eight bytes per partial cost more
|
|
// in the fill, the scatter and every level of the sort than the multiply-add that rebuilds
|
|
// them.
|
|
const auto angle_rad = [&](float img) { return (axis.GetAngle_deg(img) + wedge_half) * PI / 180.0; };
|
|
|
|
const size_t nthreads = std::max(1, settings.num_threads);
|
|
const size_t n_pts = observations.n_pts;
|
|
if (n_pts < static_cast<size_t>(settings.min_events)) return result;
|
|
const int H = static_cast<int>(observations.order.size());
|
|
const std::vector<int32_t> &bstart = observations.bstart;
|
|
const std::vector<int> &order = observations.order;
|
|
const Partial *const pts = observations.pts.get();
|
|
|
|
// Split into rocking events (same raw hkl, adjacent frames). Only >=2-frame events carry an
|
|
// unbiased phi_obs (a single-frame centroid is just the frame centre).
|
|
const float max_frame_gap = RockingEventFrameGap(axis.GetWedge_deg());
|
|
const auto run_end = [&](size_t i, size_t end) {
|
|
size_t j = i + 1;
|
|
while (j < end && pts[j].h == pts[i].h && pts[j].k == pts[i].k && pts[j].l == pts[i].l
|
|
&& pts[j].img - pts[j - 1].img <= max_frame_gap)
|
|
++j;
|
|
return j;
|
|
};
|
|
// The event the partials [i, j) make, or false where their intensities cannot place a
|
|
// centroid. Counting the events and writing them both walk the buckets, and both build the
|
|
// event this way.
|
|
const auto make_event = [&](size_t i, size_t j, Event &out) {
|
|
double sumI = 0, sumIphi = 0, sumSig = 0;
|
|
for (size_t m = i; m < j; ++m) {
|
|
const double Ipos = std::max(0.0, static_cast<double>(pts[m].I));
|
|
sumI += Ipos; sumIphi += Ipos * angle_rad(pts[m].img); sumSig += pts[m].sigma;
|
|
}
|
|
if (!(sumI > 0.0 && sumSig > 0.0)) return false;
|
|
out = Event{sumIphi / sumI, std::sqrt(sumI / sumSig), pts[i].h, pts[i].k, pts[i].l};
|
|
return true;
|
|
};
|
|
// An event never crosses an h boundary - h is the leading sort key - so the buckets can be
|
|
// walked independently, and laying their events out in bucket order gives exactly the order
|
|
// the serial walk produced. Counting first also sizes the array in one go, in place of a
|
|
// push_back that grew a gigabyte by doubling.
|
|
std::vector<int32_t> ev_count(H, 0);
|
|
std::vector<size_t> ev_frames(H, 0);
|
|
ParallelFor(H, nthreads, [&](int oi) {
|
|
const int b = order[oi];
|
|
const size_t end = bstart[b + 1];
|
|
int c = 0;
|
|
Event ev;
|
|
for (size_t i = bstart[b]; i < end; ) {
|
|
const size_t j = run_end(i, end);
|
|
if (j - i >= 2 && make_event(i, j, ev)) ++c;
|
|
i = j;
|
|
}
|
|
ev_count[b] = c;
|
|
});
|
|
std::vector<int32_t> ev_start(H + 1, 0);
|
|
for (int b = 0; b < H; ++b) ev_start[b + 1] = ev_start[b] + ev_count[b];
|
|
const size_t n_events = ev_start[H];
|
|
std::unique_ptr<Event[]> events(new Event[n_events]);
|
|
ParallelFor(H, nthreads, [&](int oi) {
|
|
const int b = order[oi];
|
|
const size_t end = bstart[b + 1];
|
|
int at = ev_start[b];
|
|
size_t frames = 0;
|
|
for (size_t i = bstart[b]; i < end; ) {
|
|
const size_t j = run_end(i, end);
|
|
if (j - i >= 2 && make_event(i, j, events[at])) { frames += j - i; ++at; }
|
|
i = j;
|
|
}
|
|
ev_frames[b] = frames;
|
|
});
|
|
size_t event_frames = 0;
|
|
for (int b = 0; b < H; ++b) event_frames += ev_frames[b];
|
|
// Frames per event is the phi_obs sampling: near 2 the reflections barely rock, so the angle
|
|
// this refinement is fitted to is under-determined. It is a geometry count, so unlike an
|
|
// intensity-weighted width it cannot be inflated by noise.
|
|
logger.Info("Post-refine: {} multi-frame rocking events ({:.1f} frames per event)", n_events,
|
|
n_events == 0 ? 0.0 : static_cast<double>(event_frames) / n_events);
|
|
if (static_cast<int>(n_events) < settings.min_events) return result;
|
|
|
|
|
|
// The rotation-scale fit further down is a single scalar whose whole point is how the residual
|
|
// varies ALONG the sweep, so it keeps every event. The cap below ranks by I/sigma, and on the
|
|
// crystals that have a stage fault the strong events sit in the middle of the sweep - the part
|
|
// that still indexes - so a capped set would leave the ends unrepresented in exactly the fit that
|
|
// has to see them. Both sets come out of the one array by selecting on indices instead: with
|
|
// the weights in the same places nth_element takes the same decisions it would take on the
|
|
// events themselves, so the selection is the same one in the same order and the whole list no
|
|
// longer has to be duplicated to survive it.
|
|
constexpr size_t MAX_EVENTS = 20000;
|
|
std::vector<int32_t> selected(n_events);
|
|
std::iota(selected.begin(), selected.end(), 0);
|
|
if (selected.size() > MAX_EVENTS) {
|
|
std::nth_element(selected.begin(), selected.begin() + MAX_EVENTS, selected.end(),
|
|
[&](int32_t a, int32_t b) { return events[a].weight > events[b].weight; });
|
|
selected.resize(MAX_EVENTS);
|
|
}
|
|
|
|
// ---- GEOMETRY REFINEMENT: ONE JOINT fit of the crystal (orientation, cell, rotation axis) and
|
|
// the detector (distance, beam centre) against BOTH residuals at once:
|
|
// * the positional detector<->reciprocal residual at each partial's observed spot, and
|
|
// * the distance-INDEPENDENT Ewald excitation residual at each rocking centroid phi_obs.
|
|
// It replaces a two-step fit that scaled the whole cell by ONE scalar against phi_obs and then
|
|
// read the distance off that scaled cell. The two were separated because the positional residual
|
|
// is degenerate with the cell scale - which is true of the positions ALONE, and is exactly what
|
|
// the excitation residual breaks, so the degeneracy that motivated the split is already resolved
|
|
// inside the same problem. Splitting it cost accuracy twice over: the first pass frees the whole
|
|
// lattice against a frozen distance, so the distortion it absorbs is ANISOTROPIC and no single
|
|
// scale can undo it; and whatever bias is left in that scale goes straight into the distance,
|
|
// which is only ever determined relative to the cell. Measured over three wavelengths of one
|
|
// crystal, that scale's bias changed SIGN with the wavelength and the distance error followed it
|
|
// with a slope of 1.3 and a correlation of 0.99, while the rocking centroids on their own fixed
|
|
// the cell volume to 0.06 %.
|
|
// Detector tilt is held fixed (gauge-coupled to the orientation on a single crystal).
|
|
// Committed only if it lowers a HELD-OUT (deterministic split-half) residual and the move stays
|
|
// within the bounds below - otherwise the geometry is left at nominal ("quit when things go wrong").
|
|
if (settings.refine_geometry) {
|
|
const gemmi::CrystalSystem sys =
|
|
(settings.crystal_system == gemmi::CrystalSystem::Trigonal) ? gemmi::CrystalSystem::Hexagonal
|
|
: settings.crystal_system;
|
|
const double ax0[3] = {ax.x, ax.y, ax.z};
|
|
const double lambda_l = lambda;
|
|
const double rot3 = nominal_geom.GetPoniRot3_rad();
|
|
// Same for every observation, so taken once here rather than per residual.
|
|
const double cos_rot3 = std::cos(rot3), sin_rot3 = std::sin(rot3);
|
|
const DetectorOrientation orientation = nominal_geom.GetOrientation();
|
|
const double pixel_mm = nominal_geom.GetPixelSize_mm();
|
|
double det_rot[2] = {nominal_geom.GetPoniRot1_rad(), nominal_geom.GetPoniRot2_rad()};
|
|
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();
|
|
|
|
// 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
|
|
// lattice centering condition (e.g. an I-centred lattice has h+k+l even for EVERY present
|
|
// reflection, so a parity split would leave the validation half empty). A reflection's rocking
|
|
// event and its spot positions carry the same hkl, so both residual families split together.
|
|
auto is_val = [](int h, int k, int l) {
|
|
unsigned u = static_cast<unsigned>(h) * 2654435761u + static_cast<unsigned>(k) * 2246822519u
|
|
+ static_cast<unsigned>(l) * 3266489917u;
|
|
u ^= u >> 15; u *= 2246822519u; u ^= u >> 13;
|
|
return (u & 1u) != 0u;
|
|
};
|
|
enum Subset { FIT, VAL, ALL };
|
|
auto in = [&](int h, int k, int l, Subset s) {
|
|
return s == ALL || (is_val(h, k, l) == (s == VAL)); };
|
|
|
|
// The crystal as XtalResidual's parameter blocks - orientation (angle-axis), cell lengths,
|
|
// cell angles - in the per-system parameterisation that residual's B matrix reads. Seeded from
|
|
// the lattice rotation indexing settled on, and free from here on: all of it, not one scale.
|
|
double p0[3] = {0, 0, 0}, p1[3] = {0, 0, 0}, p2[3] = {0, 0, 0};
|
|
double beta = r0.beta;
|
|
switch (sys) {
|
|
case gemmi::CrystalSystem::Tetragonal:
|
|
LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1);
|
|
p1[0] = (p1[0] + p1[1]) / 2.0; break;
|
|
case gemmi::CrystalSystem::Cubic:
|
|
LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1);
|
|
p1[0] = (p1[0] + p1[1] + p1[2]) / 3.0; break;
|
|
case gemmi::CrystalSystem::Hexagonal:
|
|
LatticeToRodriguesAndLengths_Hex(reference_latt, p0, p1); break;
|
|
case gemmi::CrystalSystem::Monoclinic:
|
|
LatticeToRodriguesLengthsBeta_Mono(reference_latt, p0, p1, beta);
|
|
p2[0] = beta; break;
|
|
case gemmi::CrystalSystem::Orthorhombic:
|
|
LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1); break;
|
|
default:
|
|
LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1);
|
|
p2[0] = r0.alpha * PI / 180.0; p2[1] = r0.beta * PI / 180.0; p2[2] = r0.gamma * PI / 180.0; break;
|
|
}
|
|
const double p0_0[3] = {p0[0], p0[1], p0[2]};
|
|
const double p1_0[3] = {p1[0], p1[1], p1[2]};
|
|
const double p2_0[3] = {p2[0], p2[1], p2[2]};
|
|
// Only monoclinic (beta) and triclinic read p2; for every other system it is a block the
|
|
// residual never touches, so it is held constant rather than left with a zero Jacobian column.
|
|
const bool p2_free = (sys == gemmi::CrystalSystem::Monoclinic
|
|
|| sys == gemmi::CrystalSystem::Triclinic);
|
|
|
|
// The observed spot positions, one entry per partial that produced a centroid. Count first,
|
|
// then fill, exactly as the partial gather above does and for the same reason: this walks the
|
|
// same tens of millions of partials, and a pointer vector grown by push_back copies itself
|
|
// every time it doubles. new[] rather than a sized vector so the array is not zeroed on one
|
|
// thread before the parallel fill overwrites it. The fill lands in the order the serial loop
|
|
// produced, so the selection below sees the same sequence it always did.
|
|
const int n_obs_chunks = static_cast<int>(std::clamp<size_t>(nthreads, 1,
|
|
std::max<size_t>(1, n_pts)));
|
|
const size_t obs_chunk = (n_pts + n_obs_chunks - 1) / n_obs_chunks;
|
|
std::vector<size_t> obs_offset(n_obs_chunks + 1, 0);
|
|
const auto keep_obs = [&](size_t i) {
|
|
return std::isfinite(pts[i].obs_x) && std::isfinite(pts[i].obs_y);
|
|
};
|
|
ParallelChunks(static_cast<int>(n_pts), nthreads, [&](int lo, int hi) {
|
|
size_t keep = 0;
|
|
for (int i = lo; i < hi; ++i) if (keep_obs(i)) keep++;
|
|
obs_offset[static_cast<size_t>(lo) / obs_chunk + 1] = keep;
|
|
});
|
|
for (int c = 0; c < n_obs_chunks; ++c) obs_offset[c + 1] += obs_offset[c];
|
|
size_t n_obs = obs_offset[n_obs_chunks];
|
|
std::unique_ptr<const Partial *[]> obs(new const Partial *[n_obs]);
|
|
ParallelChunks(static_cast<int>(n_pts), nthreads, [&](int lo, int hi) {
|
|
size_t at = obs_offset[static_cast<size_t>(lo) / obs_chunk];
|
|
for (int i = lo; i < hi; ++i) if (keep_obs(i)) obs[at++] = &pts[i];
|
|
});
|
|
constexpr size_t MAX_OBS = 20000;
|
|
if (n_obs > MAX_OBS) {
|
|
std::nth_element(obs.get(), obs.get() + MAX_OBS, obs.get() + n_obs,
|
|
[](const Partial *a, const Partial *b) {
|
|
return a->I / std::max(1e-9, static_cast<double>(a->sigma))
|
|
> b->I / std::max(1e-9, static_cast<double>(b->sigma)); });
|
|
n_obs = MAX_OBS;
|
|
}
|
|
result.obs_used = static_cast<int>(n_obs);
|
|
|
|
// ===== The joint fit =====
|
|
// Cost of a whole geometry over one subset, the two residual families kept apart so the log
|
|
// can say which of them moved. Each is a mean per residual VALUE, and the number they are
|
|
// summed over below is the same weighting the solver itself applies.
|
|
// se_out, where asked for: the standard error of the returned mean over the residual values
|
|
// it averages.
|
|
// exc_se_out, where asked for: the same standard error for the excitation mean alone.
|
|
auto joint_cost = [&](Subset s, const double bm[2], const double ds[1], const double rv[3],
|
|
const double q0[3], const double q1[3], const double q2[3],
|
|
double &pos_out, double &exc_out, double *se_out = nullptr,
|
|
double *exc_se_out = nullptr) {
|
|
double cp = 0.0, ce = 0.0, c4 = 0.0, c4e = 0.0;
|
|
size_t np = 0, ne = 0;
|
|
for (size_t oi = 0; oi < n_obs; ++oi) {
|
|
const Partial *pp = obs[oi];
|
|
if (!in(pp->h, pp->k, pp->l, s)) continue;
|
|
XtalResidual r(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, cos_rot3, sin_rot3,
|
|
angle_rad(pp->img), pp->h, pp->k, pp->l, sys, orientation);
|
|
double res[3] = {0, 0, 0};
|
|
r(bm, ds, det_rot, rv, q0, q1, q2, res);
|
|
for (int j = 0; j < 3; ++j) {
|
|
cp += res[j] * res[j];
|
|
c4 += res[j] * res[j] * res[j] * res[j];
|
|
}
|
|
np += 3;
|
|
}
|
|
for (const int32_t i : selected) {
|
|
const Event &ev = events[i];
|
|
if (!in(ev.h, ev.k, ev.l, s)) continue;
|
|
JointExcitationResidual r(lambda_l, ev.phi_obs, settings.excitation_weight,
|
|
ev.h, ev.k, ev.l, sys);
|
|
double res = 0.0;
|
|
r(rv, q0, q1, q2, &res);
|
|
ce += res * res;
|
|
c4 += res * res * res * res;
|
|
c4e += res * res * res * res;
|
|
++ne;
|
|
}
|
|
pos_out = np ? cp / np : 0.0;
|
|
exc_out = ne ? ce / ne : 0.0;
|
|
const double n = static_cast<double>(np + ne);
|
|
const double mean = (np + ne) ? (cp + ce) / n : 0.0;
|
|
if (se_out)
|
|
*se_out = (np + ne) > 1 ? std::sqrt(std::max(0.0, c4 / n - mean * mean) / (n - 1.0))
|
|
: std::numeric_limits<double>::quiet_NaN();
|
|
if (exc_se_out)
|
|
*exc_se_out = ne > 1 ? std::sqrt(std::max(0.0, c4e / ne - exc_out * exc_out) / (ne - 1.0))
|
|
: std::numeric_limits<double>::quiet_NaN();
|
|
return mean;
|
|
};
|
|
|
|
// 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.
|
|
// fix_distance holds the distance block at its seed, which makes the solve the nested
|
|
// hypothesis "the header distance is right" of the free one. distance_corr, where given
|
|
// (four values), receives the free fit's own formal picture of the distance: sd(distance)
|
|
// over the distance, and its correlation with each cell length, from the covariance of
|
|
// the solution.
|
|
auto solve_joint = [&](Subset s, bool wide, bool fix_distance, double bm[2], double ds[1],
|
|
double rv[3], double q0[3], double q1[3], double q2[3],
|
|
double *distance_corr = nullptr) {
|
|
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];
|
|
if (!in(pp->h, pp->k, pp->l, s)) continue;
|
|
p.AddResidualBlock(new ceres::AutoDiffCostFunction<XtalResidual, 3, 2, 1, 2, 3, 3, 3, 3>(
|
|
new XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, cos_rot3, sin_rot3,
|
|
angle_rad(pp->img), pp->h, pp->k, pp->l, sys, orientation)),
|
|
new ceres::CauchyLoss(0.02), bm, ds, det_rot, rv, q0, q1, q2);
|
|
}
|
|
for (const int32_t i : selected) {
|
|
const Event &ev = events[i];
|
|
if (!in(ev.h, ev.k, ev.l, s)) continue;
|
|
p.AddResidualBlock(new ceres::AutoDiffCostFunction<JointExcitationResidual, 1, 3, 3, 3, 3>(
|
|
new JointExcitationResidual(lambda_l, ev.phi_obs, settings.excitation_weight,
|
|
ev.h, ev.k, ev.l, sys)),
|
|
new ceres::CauchyLoss(0.02), rv, q0, q1, q2);
|
|
}
|
|
if (p.NumResidualBlocks() == 0) return false;
|
|
p.SetParameterBlockConstant(det_rot);
|
|
if (fix_distance) p.SetParameterBlockConstant(ds);
|
|
if (!p2_free) p.SetParameterBlockConstant(q2);
|
|
// Two boxes on the distance and the cell lengths. WIDE is the one this fit has always
|
|
// had, around nominal: it is what the first solve is asked in, so a fit that settles
|
|
// inside the bound is the fit this code has always committed. Otherwise the box is one
|
|
// STEP around wherever this solve starts - the trust region the walk below steps in.
|
|
// The beam, the axis and the cell angles keep their own absolute bounds around NOMINAL
|
|
// either way, since a box that travelled with the seed would let them drift a whole
|
|
// bound per step. Each of those has to reach everywhere the gate below would accept, or
|
|
// the gate is never the thing that decides - a fit pinned at a box face lands exactly ON
|
|
// 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));
|
|
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);
|
|
}
|
|
for (int j = 0; j < 3; ++j) {
|
|
p.SetParameterLowerBound(rv, j, ax0[j] - 0.05);
|
|
p.SetParameterUpperBound(rv, j, ax0[j] + 0.05);
|
|
p.SetParameterLowerBound(q1, j, wide ? 0.97 * p1_0[j]
|
|
: (1.0 - STEP_BOUND) * q1_seed[j]);
|
|
p.SetParameterUpperBound(q1, j, wide ? 1.03 * p1_0[j]
|
|
: (1.0 + STEP_BOUND) * q1_seed[j]);
|
|
}
|
|
if (p2_free)
|
|
for (int j = 0; j < 3; ++j) {
|
|
p.SetParameterLowerBound(q2, j, p2_0[j] - 0.05);
|
|
p.SetParameterUpperBound(q2, j, p2_0[j] + 0.05);
|
|
}
|
|
ceres::Solver::Options o; o.linear_solver_type = ceres::DENSE_QR; o.max_num_iterations = 60;
|
|
o.num_threads = std::max(1, settings.num_threads); o.logging_type = ceres::LoggingType::SILENT;
|
|
ceres::Solver::Summary sum; ceres::Solve(o, &p, &sum);
|
|
if (distance_corr && !fix_distance && sum.IsSolutionUsable()) {
|
|
std::fill(distance_corr, distance_corr + 4, std::numeric_limits<double>::quiet_NaN());
|
|
ceres::Covariance::Options co;
|
|
co.algorithm_type = ceres::DENSE_SVD;
|
|
co.null_space_rank = -1; // a symmetry-constrained cell leaves unused length components
|
|
ceres::Covariance cov(co);
|
|
if (cov.Compute({{ds, ds}, {ds, q1}, {q1, q1}}, &p)) {
|
|
double cdd[1], cdq[3], cqq[9];
|
|
cov.GetCovarianceBlock(ds, ds, cdd);
|
|
cov.GetCovarianceBlock(ds, q1, cdq);
|
|
cov.GetCovarianceBlock(q1, q1, cqq);
|
|
// Ceres' covariance is (J^T J)^-1 for unit-variance residuals; the residuals'
|
|
// own variance, from the fit's cost, puts it in real units.
|
|
const double sigma2 = 2.0 * sum.final_cost
|
|
/ std::max(1, sum.num_residuals_reduced - sum.num_effective_parameters_reduced);
|
|
distance_corr[0] = std::sqrt(std::max(0.0, cdd[0]) * sigma2) / ds[0];
|
|
for (int j = 0; j < 3; ++j)
|
|
distance_corr[1 + j] = cdq[j] / std::sqrt(std::max(1e-300, cdd[0] * cqq[4 * j]));
|
|
}
|
|
}
|
|
return sum.IsSolutionUsable();
|
|
};
|
|
|
|
double beam[2] = {beam_x0, beam_y0}, dist[1] = {dist0};
|
|
double axv[3] = {ax0[0], ax0[1], ax0[2]};
|
|
bool commit = false;
|
|
if (n_obs >= static_cast<size_t>(settings.min_events)) {
|
|
double pos_nom = 0.0, exc_nom = 0.0, pos_ref = 0.0, exc_ref = 0.0;
|
|
double cv_nom_se = 0.0;
|
|
const double cv_nom = joint_cost(VAL, beam, dist, axv, p0, p1, p2, pos_nom, exc_nom,
|
|
&cv_nom_se);
|
|
result.held_out_before = cv_nom;
|
|
result.held_out_before_se = cv_nom_se;
|
|
|
|
// 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]));
|
|
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];
|
|
}
|
|
// Only written when the joint solve is usable, and read unconditionally below.
|
|
double distance_corr[4] = {NAN, NAN, NAN, NAN};
|
|
bool convJ = solve_joint(FIT, /*wide=*/true, /*fix_distance=*/false, bm_f, ds_f, rv_f,
|
|
q0_f, q1_f, q2_f, distance_corr);
|
|
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)
|
|
: cv_nom;
|
|
int steps = 1;
|
|
|
|
// Do these data determine the distance at all? At a detector far enough away that
|
|
// no reflection reaches more than a few degrees of 2theta, a longer distance and a
|
|
// larger cell move every spot the same way to first order (the difference is of
|
|
// order sin^2 theta of the spot's own position). The joint fit then finds a
|
|
// distance/cell pair that fits its own spot positions a little better than the
|
|
// header, commits it, and the pass re-integrated there finds the next pair - a
|
|
// walk along the degenerate direction that the realised residual never ratifies.
|
|
// So the same fit is asked once more with the distance HELD at the header, every
|
|
// other block as free as before - the nested hypothesis "the header distance is
|
|
// right" - and the two are compared on the one residual family that can tell them
|
|
// apart. The excitation residual never involves the detector, so it is blind to the
|
|
// distance itself; what it sees is the cell scale, and a held fit at a wrong header
|
|
// distance is forced into a wrong cell scale by the spot positions, which the
|
|
// rocking angles then refuse (measured: a header 1.4 % long leaves the held fit's
|
|
// excitation residual seventeen times the free fit's). Where freeing the distance
|
|
// lowers the held-out excitation residual below the held fit's by more than the
|
|
// standard error of that residual - the same noise a round of the geometry walk has
|
|
// to beat (HeldOutResidualFell) - the free fit is committed exactly as before. Where
|
|
// it does not, the held fit is the candidate from here on: what the walk, the gate
|
|
// and the commit below see. The positional residual is deliberately not consulted
|
|
// for this. Its in-fit gain along the degenerate direction is the one re-integration
|
|
// erases, and pooled with the excitation family it either drowns a decisive
|
|
// excitation gain in its own noise (measured: a 54 % excitation gain read as 9 %
|
|
// pooled against a 9 % noise) or lends the degenerate direction a gain that is not
|
|
// there.
|
|
double bm_h[2] = {beam_x0, beam_y0}, ds_h[1] = {dist0};
|
|
double rv_h[3], q0_h[3], q1_h[3], q2_h[3];
|
|
for (int j = 0; j < 3; ++j) {
|
|
rv_h[j] = ax0[j]; q0_h[j] = p0_0[j]; q1_h[j] = p1_0[j]; q2_h[j] = p2_0[j];
|
|
}
|
|
// A run that has already walked off the header has refuted that hypothesis by
|
|
// re-integrating - the rounds of a walk are judged on what they realise - so the fit
|
|
// is not asked to hold a distance that is no longer the header's.
|
|
const bool convH = settings.distance_at_header && convJ
|
|
&& solve_joint(FIT, /*wide=*/true, /*fix_distance=*/true, bm_h, ds_h,
|
|
rv_h, q0_h, q1_h, q2_h);
|
|
double pos_h = 0.0, exc_h = 0.0, cv_h = cv_nom, exc_h_se = NAN;
|
|
if (convH)
|
|
cv_h = joint_cost(VAL, bm_h, ds_h, rv_h, q0_h, q1_h, q2_h, pos_h, exc_h, nullptr,
|
|
&exc_h_se);
|
|
const double noise = std::hypot(exc_ref_se, exc_h_se);
|
|
const bool hold_distance = convH && !(std::isfinite(noise) && exc_ref < exc_h - noise);
|
|
if (settings.distance_at_header)
|
|
logger.Info("Post-refine GEOM distance: free {:.3f} mm gives held-out excitation {:.3e} +- "
|
|
"{:.1e} (positional {:.3e}); held at the header {:.3f} mm gives {:.3e} +- "
|
|
"{:.1e} (positional {:.3e}); free minus held {:+.3e} against the noise {:.1e} "
|
|
"=> {}; the free fit's formal sd(distance) {:.4f} %, corr(distance, cell "
|
|
"a/b/c) {:.4f} {:.4f} {:.4f}",
|
|
ds_f[0], exc_ref, exc_ref_se, pos_ref, dist0, exc_h, exc_h_se, pos_h,
|
|
exc_ref - exc_h, noise,
|
|
hold_distance ? "the rocking angles do not ask for it, the header distance is held"
|
|
: convH ? "the free distance pays"
|
|
: "the held fit did not converge, the free distance stands",
|
|
100.0 * distance_corr[0], distance_corr[1], distance_corr[2],
|
|
distance_corr[3]);
|
|
if (hold_distance) {
|
|
bm_f[0] = bm_h[0]; bm_f[1] = bm_h[1]; ds_f[0] = ds_h[0];
|
|
for (int j = 0; j < 3; ++j) {
|
|
rv_f[j] = rv_h[j]; q0_f[j] = q0_h[j]; q1_f[j] = q1_h[j]; q2_f[j] = q2_h[j];
|
|
}
|
|
cv_ref = cv_h; pos_ref = pos_h; exc_ref = exc_h;
|
|
}
|
|
|
|
// 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];
|
|
}
|
|
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, hold_distance, 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);
|
|
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::numeric_limits<double>::infinity();
|
|
if (std::min(from_nominal, from_measured) >= BEAM_BOUND_PXL)
|
|
return "the beam moved further than the bound from every centre anything believes";
|
|
return nullptr;
|
|
};
|
|
|
|
// The two residual families are asked separately as well as together. Pooled, the
|
|
// positional values outnumber the excitation ones about three to one where both caps
|
|
// saturate, and the excitation residual is the only one that identifies the cell SCALE -
|
|
// so a pooled mean can improve on the strength of the positions alone while the one
|
|
// quantity the cell is committed for has got worse. Neither may degrade.
|
|
const char *refused =
|
|
!convJ ? "the fit did not converge"
|
|
// A walk that used every step it was given never settled: it stopped because it
|
|
// ran out of steps, not because it had arrived, and where that happens it is
|
|
// walking away rather than to something. This is the runaway the size of the
|
|
// move used to stand in for, tested where it can actually be seen - measured on
|
|
// a weak 7 A crystal whose fit asked for another 1 % at every one of five steps
|
|
// and 4.3 % in all, against the 1.2 % its cell is actually out by.
|
|
: steps >= 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, 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
|
|
// actually be committed, not only of the half that passed the gate. Refused here, the
|
|
// run keeps its header geometry exactly as the other refusals leave it.
|
|
double ds_log[1] = {ds_f[0]}, bm_log[2] = {bm_f[0], bm_f[1]};
|
|
double q1_log[3], q2_log[3];
|
|
for (int j = 0; j < 3; ++j) { q1_log[j] = q1_f[j]; q2_log[j] = q2_f[j]; }
|
|
if (refused == nullptr) {
|
|
// A walk is re-fitted from where it arrived, one step wide, so the re-fit is its
|
|
// last step taken over all the data rather than the whole walk taken again on it.
|
|
// Everything else is re-fitted from nominal in the wide box, exactly as before.
|
|
if (walked) {
|
|
beam[0] = bm_f[0]; beam[1] = bm_f[1]; dist[0] = ds_f[0];
|
|
for (int j = 0; j < 3; ++j) {
|
|
axv[j] = rv_f[j]; p0[j] = q0_f[j]; p1[j] = q1_f[j]; p2[j] = q2_f[j];
|
|
}
|
|
}
|
|
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);
|
|
if (refused != nullptr) {
|
|
beam[0] = beam_x0; beam[1] = beam_y0; dist[0] = dist0;
|
|
for (int j = 0; j < 3; ++j) {
|
|
axv[j] = ax0[j]; p0[j] = p0_0[j]; p1[j] = p1_0[j]; p2[j] = p2_0[j];
|
|
}
|
|
}
|
|
}
|
|
commit = (refused == nullptr);
|
|
result.distance_held = commit && hold_distance;
|
|
// Name which test refused it. Several different things reject here and the geometry that
|
|
// comes out is the same in all of them, so a run that silently keeps its header geometry
|
|
// says nothing about whether the fit was bad, the improvement too small, or the move too
|
|
// large for the bound - which is the one case where the number worth reading is the one
|
|
// that was thrown away. The report carries the same sentence, so a refusal is visible
|
|
// without the log.
|
|
const std::string verdict = commit ? fmt::format("COMMIT ({} step{}{})", steps,
|
|
steps == 1 ? "" : "s",
|
|
hold_distance ? ", distance held at the header"
|
|
: "")
|
|
: "reject (" + std::string(refused) + ")";
|
|
// The cell the log names is the EFFECTIVE one - what the residual's B matrix builds
|
|
// from the blocks - not the blocks themselves, whose unused components a high-symmetry
|
|
// system leaves at whatever the seed happened to put there. The geometry beside the
|
|
// verdict is the one the verdict is about: the all-data re-fit where it committed, the
|
|
// half-data fit where it did not. The residual pair stays the split-half gate's, since
|
|
// the fit that is committed has no held-out half of its own.
|
|
double len_nom[3], ang_nom[3], len_fit[3], ang_fit[3];
|
|
EffectiveCellFromParams(sys, p1_0, p2_0, len_nom, ang_nom);
|
|
EffectiveCellFromParams(sys, q1_log, q2_log, len_fit, ang_fit);
|
|
result.cell_before = UnitCell{static_cast<float>(len_nom[0]), static_cast<float>(len_nom[1]),
|
|
static_cast<float>(len_nom[2]),
|
|
static_cast<float>(ang_nom[0] * 180.0 / PI),
|
|
static_cast<float>(ang_nom[1] * 180.0 / PI),
|
|
static_cast<float>(ang_nom[2] * 180.0 / PI)};
|
|
logger.Info("Post-refine GEOM (joint crystal + detector): dist {:.3f} -> {:.3f} mm, beam "
|
|
"({:.2f},{:.2f}) -> ({:.2f},{:.2f}), cell {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} "
|
|
"{:.2f} -> {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}, held-out positional "
|
|
"{:.3e} -> {:.3e}, excitation {:.3e} -> {:.3e} => {}",
|
|
dist0, ds_log[0], beam_x0, beam_y0, bm_log[0], bm_log[1],
|
|
len_nom[0], len_nom[1], len_nom[2],
|
|
ang_nom[0] * 180.0 / PI, ang_nom[1] * 180.0 / PI, ang_nom[2] * 180.0 / PI,
|
|
len_fit[0], len_fit[1], len_fit[2],
|
|
ang_fit[0] * 180.0 / PI, ang_fit[1] * 180.0 / PI, ang_fit[2] * 180.0 / PI,
|
|
pos_nom, pos_ref, exc_nom, exc_ref, verdict);
|
|
if (!commit)
|
|
result.refused_reason = fmt::format(
|
|
"{} - it wanted {:.3f} mm and cell {:.3f} {:.3f} {:.3f}", refused,
|
|
ds_log[0], len_fit[0], len_fit[1], len_fit[2]);
|
|
// A move of more than one step is a walk, not a correction next to where it started,
|
|
// and it is the one the old bound refused outright. It is committed here and ratified
|
|
// where it can be: by re-indexing at it (see RunAllPasses).
|
|
// A held distance has nothing to walk: the walk is the distance finding the rest of an
|
|
// error the rocking angles have just said they do not see.
|
|
result.large_move = commit && !hold_distance && beyond_a_step(dist, p1);
|
|
} else {
|
|
logger.Info("Post-refine GEOM: only {} positional observations - the joint fit needs the "
|
|
"spot positions as well as the rocking angles, so nothing is refined", n_obs);
|
|
}
|
|
result.cell_refined = commit;
|
|
result.detector_refined = commit;
|
|
|
|
const double axlen = std::sqrt(axv[0]*axv[0] + axv[1]*axv[1] + axv[2]*axv[2]);
|
|
const double axdev = std::acos(std::clamp((axv[0]*ax0[0]+axv[1]*ax0[1]+axv[2]*ax0[2])
|
|
/ std::max(1e-9, axlen), -1.0, 1.0)) * 180.0 / PI;
|
|
// The committed crystal, as a lattice again. Where nothing was committed this is the lattice
|
|
// rotation indexing handed in, unchanged - not its round trip through the symmetry-constrained
|
|
// parameterisation, which would move the cell for a fit that was refused.
|
|
double eff_len[3], eff_ang[3];
|
|
EffectiveCellFromParams(sys, p1, p2, eff_len, eff_ang);
|
|
const CrystalLattice committed_latt = commit
|
|
? AngleAxisAndCellToLattice(p0, eff_len, eff_ang[0], eff_ang[1], eff_ang[2])
|
|
: reference_latt;
|
|
const Coord As = committed_latt.Astar(), Bs = committed_latt.Bstar(), Cs = committed_latt.Cstar();
|
|
if (commit) {
|
|
result.cell = UnitCell{static_cast<float>(eff_len[0]), static_cast<float>(eff_len[1]),
|
|
static_cast<float>(eff_len[2]),
|
|
static_cast<float>(eff_ang[0] * 180.0 / PI),
|
|
static_cast<float>(eff_ang[1] * 180.0 / PI),
|
|
static_cast<float>(eff_ang[2] * 180.0 / PI)};
|
|
logger.Info("Post-refine GEOM: committed cell {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}, "
|
|
"rotation axis moved {:.3f} deg (the cell is the fit's own; the second pass "
|
|
"re-indexes at the refined detector geometry)", result.cell.a, result.cell.b,
|
|
result.cell.c, result.cell.alpha, result.cell.beta, result.cell.gamma, axdev);
|
|
}
|
|
|
|
// ===== Goniometer rotation SCALE k, its own one-parameter fit on the same rocking events =====
|
|
// The angles stored in the file are the COMMANDED ones, so a stage that turned k times as far
|
|
// is invisible in the header. Nothing else here can represent it: the cell scale, the axis
|
|
// direction, the distance and the beam are all orthogonal to a rotation MAGNITUDE error. Fitted
|
|
// after step A so the cell scale and the axis direction are fixed at their committed values and
|
|
// k is the only free quantity.
|
|
const double u[3] = {axv[0] / axlen, axv[1] / axlen, axv[2] / axlen};
|
|
double phi_c = 0.0, phi_lo = events[0].phi_obs, phi_hi = events[0].phi_obs;
|
|
for (size_t e = 0; e < n_events; ++e) {
|
|
phi_c += events[e].phi_obs;
|
|
phi_lo = std::min(phi_lo, events[e].phi_obs);
|
|
phi_hi = std::max(phi_hi, events[e].phi_obs);
|
|
}
|
|
phi_c /= static_cast<double>(n_events);
|
|
const double sweep_deg = (phi_hi - phi_lo) * 180.0 / PI;
|
|
// The committed reciprocal vector turned to the sweep centre. The
|
|
// angle then enters the fit measured FROM that centre. A constant crystal missetting about the
|
|
// spindle is k with a slope in phi, so measuring the angle from the goniometer's zero instead
|
|
// lets a missetting leak into k with gain <phi>/<phi^2> - which depends only on where the sweep
|
|
// happens to sit. On a short sweep starting near zero that gain is enormous: a 0.14 deg
|
|
// missetting on a 10 deg wedge fakes 1.4 % of k. Referred to the sweep centre the leak is
|
|
// identically zero at any width, and no parameter has to be added to get it.
|
|
// The residual is closed-form in k, so this is a one-parameter minimisation rather than a
|
|
// solver problem. A rotation preserves length, so |p_lab| = |e_mid| whatever k is, and only
|
|
// the z component moves; Rodrigues gives it exactly:
|
|
//
|
|
// r(k) = C + A cos(a k) - B sin(a k) = C + R cos(a k + psi)
|
|
// C = lambda |e|^2 / 2 + u_z (u.e), A = e_z - u_z (u.e), B = (u x e)_z, a = phi_obs - phi_c
|
|
//
|
|
// which is the same function the residual functor computes, to the last bit. Handing 8 million
|
|
// one-parameter residual blocks to Ceres instead cost tens of millions of allocations and a
|
|
// dense factorisation per iteration, for a fit that a scan over a bounded interval settles.
|
|
// Coefficients are computed in double and stored narrowed: their rounding perturbs the
|
|
// minimiser by ~1e-10, and k is carried downstream as a float.
|
|
struct ScaleTerm { float a, C, R, psi; };
|
|
std::vector<ScaleTerm> terms(n_events);
|
|
std::vector<int> fifth_of(n_events);
|
|
const double aa_c[3] = {-phi_c * u[0], -phi_c * u[1], -phi_c * u[2]};
|
|
ParallelChunks(static_cast<int>(n_events), nthreads, [&](int lo, int hi) {
|
|
for (int e = lo; e < hi; ++e) {
|
|
const Coord ec = As * static_cast<float>(events[e].h)
|
|
+ Bs * static_cast<float>(events[e].k)
|
|
+ Cs * static_cast<float>(events[e].l);
|
|
const double p[3] = {ec.x, ec.y, ec.z};
|
|
double em[3];
|
|
ceres::AngleAxisRotatePoint(aa_c, p, em);
|
|
const double ue = u[0] * em[0] + u[1] * em[1] + u[2] * em[2];
|
|
const double e2 = em[0] * em[0] + em[1] * em[1] + em[2] * em[2];
|
|
const double C = 0.5 * lambda_l * e2 + u[2] * ue;
|
|
const double A = em[2] - u[2] * ue;
|
|
const double B = u[0] * em[1] - u[1] * em[0];
|
|
terms[e] = ScaleTerm{static_cast<float>(events[e].phi_obs - phi_c),
|
|
static_cast<float>(C), static_cast<float>(std::hypot(A, B)),
|
|
static_cast<float>(std::atan2(B, A))};
|
|
fifth_of[e] = std::clamp(static_cast<int>(
|
|
5.0 * (events[e].phi_obs - phi_lo) / std::max(1e-9, phi_hi - phi_lo)), 0, 4);
|
|
}
|
|
});
|
|
|
|
// Robust-loss scale from the scatter the events actually have: it varies by more than a decade
|
|
// between datasets, so any fixed constant is either inert or throws away real data. Taken once,
|
|
// over every event, so the all-data fit and every jackknife fold share it.
|
|
const auto residual_at = [&](const ScaleTerm &t, double k) {
|
|
return static_cast<double>(t.C)
|
|
+ static_cast<double>(t.R) * std::cos(static_cast<double>(t.a) * k + t.psi);
|
|
};
|
|
double rms = 0.0;
|
|
for (const auto &t : terms) {
|
|
const double r = residual_at(t, 1.0);
|
|
rms += r * r;
|
|
}
|
|
rms = std::sqrt(rms / static_cast<double>(terms.size()));
|
|
const double huber_delta = std::max(1e-12, 2.0 * rms);
|
|
const double huber_d2 = huber_delta * huber_delta;
|
|
|
|
// Ceres minimises half the sum of the loss applied to the SQUARED residual, so that is what is
|
|
// reproduced here. One pass yields the five per-fifth partial sums, which serve the all-data
|
|
// fit and all five leave-a-fifth-out folds together.
|
|
// Each chunk folds into its own slot and the slots are summed in chunk order, so the
|
|
// sum does not depend on which worker finishes first: the same events always add up in
|
|
// the same sequence, and the fit is reproducible run to run.
|
|
const int n_terms = static_cast<int>(terms.size());
|
|
const int cost_nt = static_cast<int>(std::max<size_t>(1, std::min(nthreads,
|
|
static_cast<size_t>(n_terms))));
|
|
const int cost_chunk = (n_terms + cost_nt - 1) / cost_nt;
|
|
// Several k are always wanted at once (the grid below asks for 101), and they all sweep the
|
|
// same event list, so sweep it ONCE and evaluate every k on each event while it is still in
|
|
// registers. The per-thread accumulator is one slot per (k, fifth) - 4 kB for the grid, small
|
|
// enough to stay in L1 - against re-reading the whole term array once per k. Each (k, fifth)
|
|
// still receives its events in the same order and the chunks are still summed in chunk order,
|
|
// so the sums are the ones a k-at-a-time loop produced, bit for bit.
|
|
const auto cost_grid = [&](const std::vector<double> &ks) {
|
|
const int nk = static_cast<int>(ks.size());
|
|
std::vector<std::vector<std::array<double, 5>>> per_chunk(
|
|
cost_nt, std::vector<std::array<double, 5>>(nk));
|
|
ParallelChunks(n_terms, nthreads, [&](int lo, int hi) {
|
|
std::vector<std::array<double, 5>> acc(nk);
|
|
for (int e = lo; e < hi; ++e) {
|
|
const ScaleTerm &term = terms[e];
|
|
const int fifth = fifth_of[e];
|
|
for (int g = 0; g < nk; ++g) {
|
|
const double r = residual_at(term, ks[g]);
|
|
const double s2 = r * r;
|
|
acc[g][fifth] += (s2 <= huber_d2) ? s2
|
|
: (2.0 * huber_delta * std::sqrt(s2) - huber_d2);
|
|
}
|
|
}
|
|
per_chunk[lo / cost_chunk] = std::move(acc);
|
|
});
|
|
std::vector<std::array<double, 5>> total(nk);
|
|
for (const auto &acc : per_chunk)
|
|
for (int g = 0; g < nk; ++g)
|
|
for (int j = 0; j < 5; ++j) total[g][j] += acc[g][j];
|
|
return total;
|
|
};
|
|
const auto cost_by_fifth = [&](double k) { return cost_grid({k})[0]; };
|
|
|
|
// Scan the interval Ceres was bounded to, then close in. No event's phase can move by more than
|
|
// a fraction of a period over an interval this narrow, so the objective has no structure the
|
|
// grid could step over; the refinement is only there to place the minimum precisely.
|
|
constexpr int SCALE_GRID = 101;
|
|
constexpr double SCALE_K_LO = 0.95, SCALE_K_HI = 1.05;
|
|
std::vector<double> grid_k(SCALE_GRID);
|
|
for (int g = 0; g < SCALE_GRID; ++g)
|
|
grid_k[g] = SCALE_K_LO + (SCALE_K_HI - SCALE_K_LO) * g / (SCALE_GRID - 1);
|
|
const std::vector<std::array<double, 5>> grid = cost_grid(grid_k);
|
|
auto solve_scale = [&](int drop_fifth) {
|
|
const auto total = [&](const std::array<double, 5> &f) {
|
|
double t = 0.0;
|
|
for (int j = 0; j < 5; ++j)
|
|
if (j != drop_fifth) t += f[j];
|
|
return t;
|
|
};
|
|
int best = 0;
|
|
for (int g = 1; g < SCALE_GRID; ++g)
|
|
if (total(grid[g]) < total(grid[best])) best = g;
|
|
const double step = (SCALE_K_HI - SCALE_K_LO) / (SCALE_GRID - 1);
|
|
double a = std::max(SCALE_K_LO, SCALE_K_LO + step * (best - 1));
|
|
double b = std::min(SCALE_K_HI, SCALE_K_LO + step * (best + 1));
|
|
// Golden section: the objective is smooth but its curvature jumps wherever an event
|
|
// crosses the Huber knee, which a derivative method would have to cope with.
|
|
constexpr double INV_PHI = 0.6180339887498949;
|
|
double c = b - INV_PHI * (b - a), d = a + INV_PHI * (b - a);
|
|
double fc = total(cost_by_fifth(c)), fd = total(cost_by_fifth(d));
|
|
// The fit is narrowed to a float before it is applied (Rugnux.h prepass_rotation_scale_),
|
|
// so bracketing it below that type's epsilon, 6e-8, only buys about ten more full
|
|
// passes over the events for a digit that cannot survive being stored.
|
|
while (b - a > 1e-7) {
|
|
if (fc < fd) { b = d; d = c; fd = fc; c = b - INV_PHI * (b - a); fc = total(cost_by_fifth(c)); }
|
|
else { a = c; c = d; fc = fd; d = a + INV_PHI * (b - a); fd = total(cost_by_fifth(d)); }
|
|
}
|
|
return 0.5 * (a + b);
|
|
};
|
|
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;
|
|
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;
|
|
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");
|
|
|
|
// Assemble the committed geometry.
|
|
result.distance_after_mm = dist[0];
|
|
result.beam_x_before_px = beam_x0; result.beam_x_after_px = beam[0];
|
|
result.beam_y_before_px = beam_y0; result.beam_y_after_px = beam[1];
|
|
result.events_used = static_cast<int>(selected.size());
|
|
result.ok = commit;
|
|
if (!result.ok)
|
|
logger.Info("Post-refine GEOM: the joint fit did not pass cross-validation - geometry "
|
|
"left at nominal");
|
|
return result;
|
|
}
|
|
|
|
return result; // refine_geometry is the only supported mode; nothing refined otherwise
|
|
} catch (...) {
|
|
result.ok = false;
|
|
return result;
|
|
}
|
|
}
|
|
|
|
|
|
bool HeldOutResidualFell(const PostRefineResult &before, const PostRefineResult &after) {
|
|
const double noise = std::hypot(before.held_out_before_se, after.held_out_before_se);
|
|
return std::isfinite(noise) && after.held_out_before < before.held_out_before - noise;
|
|
}
|
|
|
|
int ReindexPushesCellBack(const PostRefineResult &fit, const PostRefineResult &next) {
|
|
if (!fit.cell_refined)
|
|
return 0;
|
|
const double from[3] = {fit.cell_before.a, fit.cell_before.b, fit.cell_before.c};
|
|
const double to[3] = {fit.cell.a, fit.cell.b, fit.cell.c};
|
|
const double again[3] = {next.cell_before.a, next.cell_before.b, next.cell_before.c};
|
|
int axis = -1;
|
|
double largest = 0.0;
|
|
for (int j = 0; j < 3; ++j)
|
|
if (from[j] > 0.0 && std::fabs(to[j] - from[j]) / from[j] > largest) {
|
|
largest = std::fabs(to[j] - from[j]) / from[j];
|
|
axis = j;
|
|
}
|
|
if (axis < 0 || !(again[axis] > 0.0))
|
|
return 0;
|
|
const double move = to[axis] - from[axis];
|
|
const double back = again[axis] - to[axis];
|
|
if (!(move * back < 0.0))
|
|
return 0;
|
|
return move > 0.0 ? axis + 1 : -(axis + 1);
|
|
}
|