rotation: per-frame refinement restrained toward the sweep (empirical Bayes)
An integrated rotation frame is fitted twice, each fit independent of the other frames. The first fit is free and sets the frame's residual scale. The second starts again from the sweep's geometry with a Gaussian prior on how far the orientation and the beam centre move. The prior's two widths are measured once per run, on the geometry pre-pass, as the rms of the seven-frame running median of the frames' free shifts. This is a pass-level input, like the pre-pass mosaicity. A beam centre that does not wander in the pre-pass is therefore held near the sweep, and the data decide this per run without a switch. The orientation prior acts on the rotation away from the sweep orientation and not on the angle-axis parameters. Applied to the parameters it was up to 2.4x too tight, depending on how the crystal sits in the laboratory. Both fits use XtalOptimizer's back-rotating form with the spindle given, so the beam-centre component along the spindle, which on a single frame nearly trades against the orientation, gets the same restraint it gets over a sweep. The frames the first pass scores lattices on, the pre-pass and stills keep the free fit. The residual scale counts 3N components and has no parameter correction (measured as immaterial). Squashes the sparse-eb prototype (2a54995fd, 8fb05679a, 51ae67467) onto twin-subgroup-check, simplified as described above. On sparse small molecules SHELXL R1 improved by a geometric mean of 3.943 -> 3.768 over 11 well-measured sets. Protein REFMAC R-free and anomalous peaks stayed within noise. The rotation-scale walk is unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK (cherry picked from commit 045ae107d1c6b6ab6dc40773924426134dee69bf) (cherry picked from commit 1f3c11f61c693e6eb5debbbb881245c69efe7472)
This commit is contained in:
leonarski_f committed
1 parent
cc1553a767
commit
37c3661b33
13 files changed
+355
-13
No files matched your search
@@ -356,6 +356,8 @@ Refinement is performed in stages with decreasing acceptance tolerance for inclu
|
||||
|
||||
The loose first stage necessarily admits some spots that are not reflections of this lattice — the fraction of *randomly* placed spots inside a fractional-Miller tolerance $t$ is $\tfrac{4}{3}\pi t^3$, i.e. 11 % at $t=0.3$ — and an unweighted fit lets them pull the orientation. Each residual is therefore weighted by how strong its spot is **for its resolution**: the frame's spots are cut into equal-count resolution shells and each intensity is divided by its shell median, mapped to $w^2=r/(1+r)$. The shell normalisation is what makes this safe — genuine high-resolution spots are legitimately weaker and carry the cell and distance information, so an un-normalised intensity weight would suppress exactly the spots the fit needs. The weight is a property of the spot and never of the current residual, so it does not depend on how far the geometry is from convergence.
|
||||
|
||||
**A rotation frame's refinement is restrained toward the sweep.** A frame is integrated on the orientation and beam centre refined against its own spots, starting from the sweep's geometry. Where the crystal really moves the refined orientation follows it, and integrating on it is what keeps such a sweep together. Where a frame has few or poorly placed spots the same fit can instead jump on its own, by up to a few degrees, which drops strong reflections out of prediction on the frame they peak on. So after the geometry pre-pass, which refines its frames freely, each integrated frame is fitted twice: freely, for the scale $\sigma^2$ of its own weighted misses, and then from the sweep again with a Gaussian prior on how far it moves, minimising $\sum_i w_i^2\|\mathbf{r}_i\|^2/\sigma^2 + \sum_k \theta_k^2/\tau_k^2$ over the orientation and beam-centre shifts $\theta_k$ (the orientation's as the three components of the rotation from the sweep's orientation, so that its width means the same in every direction however the crystal sits in the laboratory). Both fits of an integrated frame are made with the spindle given, as over a sweep, so the beam centre's component along the spindle - which on a single frame nearly trades against the orientation - carries the same weak restraint toward the start that it carries in the sweep refinement. The widths $\tau$ - one for the orientation, one for the beam centre - are measured once per run on the pre-pass frames, as the root mean square of the seven-frame running median of their free shifts (the track's smooth part, which a jump on one frame does not reach). There is no gate: the data of each frame decide how far it moves, a frame with many well placed spots following a moving crystal and a sparse one staying near the sweep. Each frame is still refined on its own spots, independently of the others. The frames the first pass scores candidate lattices on and the pre-pass keep the free refinement, so the lattice and the geometry the run refines do not depend on the prior, and stills are not affected.
|
||||
|
||||
**The rotation chain commits its best round, not the round it stops on.** After the winning candidate is selected, the refinement is run again — solve, re-accumulate the reciprocal-space cloud under the refined geometry, solve again — up to twenty times, and the loop stops on a test of the detector-tilt step. The chain is a trajectory and its last point is not always its best one: every solve ends by fitting only the spots inside its **tightest** gate, so a cell with a direction the data barely constrain (which is what a free cell whose metric is near a Bravais class has) can slide along it, pulling a core of spots tighter while the periphery falls out of the fit altogether. Measured on such a chain, the spots inside the tight gate rise over seventeen rounds while the spots inside the widest gate peak at round three and fall away — and round three is the round that merges at ISa 11.0 against 6.3 and $R_\mathrm{meas}$ 0.148 against 0.213. Nothing the run consulted could see it: indexed fraction, validation frames, validation spots and the tight-gate count all prefer the overfitted end.
|
||||
|
||||
So every round is scored on the **widest** gate — the population the first pass selects on and the last pass does not fit, which makes it the one a converged solve is not optimising — and the best-scoring round is committed. Two conditions keep that from acting on noise. The score is a *count* of spots, so a lead of fewer than $\sqrt{\text{count}}$ of them leaves the last round standing. And the round taken has to be the less distorted lattice as well as the better-fitting one: the lattice search (§6) is re-asked every round to **measure** how far the cell sits from the ideal metric of the class it matches (imposing that class is measured fatal — the snap puts almost everything outside the refinement's own gate), and an earlier round is taken only when it matched the *same* class and sits closer to it. Same class is a precondition and not a precaution: the deviation is a fraction of whichever class's tolerance admitted it, so two classes' deviations are not the same quantity, and a round that matched no class reports zero, which means "nothing was asserted" rather than "undistorted". A chain that has settled scores its rounds within a spot or two of each other and a symmetry-constrained solve holds its distortion at zero throughout, so the rule fires on neither: measured over 914 chains, an earlier round scores higher on 44 % of them and the committed round changes on 1 dataset in 54.
|
||||
|
||||
+1
-1
@@ -268,7 +268,7 @@ For very long axes the direction grid's angular resolution binds well below 500
|
||||
| `--fft-min-unit-cell <num>` | Shortest cell axis the FFT search accepts, Å (default: 10). A candidate with a shorter axis is discarded — but on rotation data a row shorter than the floor, down to 3 Å, is recognised under the multiple of it the search reports past the floor, so a small-molecule cell is indexed without this flag. `-C` lowers the floor on its own to cover the cell given. On stills a crystal below the floor cannot be indexed unless this is lowered |
|
||||
| `--min-indexed-spots <num>` | Spots a frame must have on the lattice before it counts as indexed (default: 6, minimum 4). It sets the reported indexing rate and the count the rotation first pass ranks candidate lattices by; whether the run has a lattice at all is decided on the pooled spots instead ([§4.1](CPU_DATA_ANALYSIS_INDEXING.md#41-indexed-spot-decision-inlier-test)), and integration is not gated by it |
|
||||
| `-S, --space-group <num\|symbol>` | Space group number (`96`) or Hermann-Mauguin symbol (`P43212`) — for indexing and scaling. On serial stills whose group admits two indexings, it also lets the crystals be put on one from the data ([the indexing ambiguity](RUGNUX_ADVANCED.md#the-indexing-ambiguity)) |
|
||||
| `-r, --geom-refine <txt>` | Geometry refinement: `none` \| `orientation` \| `beam_and_lattice` (default) \| `flex` (try all three per image, keep whichever indexes the most spots; alias `multi`) |
|
||||
| `-r, --geom-refine <txt>` | Geometry refinement: `none` \| `orientation` \| `beam_and_lattice` (default; on rotation data restrained toward the sweep's geometry, [§7.4](CPU_DATA_ANALYSIS_INDEXING.md#74-multi-stage-tightening-of-inlier-tolerance)) \| `flex` (try all three per image, keep whichever indexes the most spots; alias `multi`) |
|
||||
| `-R, --two-pass-rotation[=num]` | Two-pass offline rotation indexing (default for goniometer data; optional first-pass image count, default 100) |
|
||||
| `--force-rotation-lattice <vec>` | Force rotation lattice (9 floats, Å), skipping the first pass |
|
||||
| `--rotation-no-postrefine` | Rotation: disable the default-on two-pass geometry post-refine (see [Rotation data](RUGNUX_TUTORIAL.md#rotation-data)) |
|
||||
|
||||
@@ -7,6 +7,7 @@
|
||||
#include <cstdlib>
|
||||
#include <unordered_map>
|
||||
#include "IndexAndRefine.h"
|
||||
#include <Eigen/Dense>
|
||||
|
||||
#include "bragg_integration/CalcISigma.h"
|
||||
#include "bragg_integration/SystematicAbsence.h"
|
||||
@@ -418,9 +419,35 @@ int CountIndexedSpots(const DiffractionGeometry &geom, const CrystalLattice &lat
|
||||
}
|
||||
return n;
|
||||
}
|
||||
|
||||
// One rotation frame refined through the back-rotating form of XtalOptimizer, the spindle given: the lattice
|
||||
// (and the restraint's centre) is turned back to the sweep's zero angle, refined against the frame's spots
|
||||
// back-rotated by the frame's own angle, and turned forward again - the identity on the orientation - but
|
||||
// with the axis known the fit holds the beam centre's spindle-parallel component to the start with the
|
||||
// restraint XtalOptimizer applies over a sweep. Over a sweep that component is a gauge; on one frame it
|
||||
// nearly is, trading against the crystal's orientation, and a frame with few spots spends its freedom there.
|
||||
bool XtalOptimizerOnSpindle(XtalOptimizerData &data, const std::vector<SpotToSave> &spots, int64_t frame,
|
||||
const GoniometerAxis &gon) {
|
||||
const float angle_deg = gon.GetAngle_deg(frame) + gon.GetWedge_deg() / 2.0f;
|
||||
XtalOptimizerData d = data;
|
||||
d.axis = GoniometerAxis(gon.GetName(), gon.GetAngle_deg(frame), gon.GetIncrement_deg(), gon.GetAxis(),
|
||||
gon.GetHelicalStep()).ScreeningWedge(gon.GetScreeningWedge());
|
||||
d.latt = data.latt.Multiply(gon.GetTransformationAngle(angle_deg));
|
||||
if (data.prior_lattice)
|
||||
d.prior_lattice = data.prior_lattice->Multiply(gon.GetTransformationAngle(angle_deg));
|
||||
if (!XtalOptimizer(d, std::span(&spots, 1)))
|
||||
return false;
|
||||
data.geom = d.geom;
|
||||
data.beam_corr_x = d.beam_corr_x;
|
||||
data.beam_corr_y = d.beam_corr_y;
|
||||
data.fit_residual_var = d.fit_residual_var;
|
||||
data.latt = d.latt.Multiply(gon.GetTransformationAngle(-angle_deg));
|
||||
return true;
|
||||
}
|
||||
} // namespace
|
||||
|
||||
void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::IndexingOutcome &outcome) {
|
||||
void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::IndexingOutcome &outcome,
|
||||
bool integrating) {
|
||||
if (!outcome.lattice_candidate)
|
||||
return;
|
||||
|
||||
@@ -462,6 +489,9 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In
|
||||
if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Trigonal)
|
||||
data.crystal_system = gemmi::CrystalSystem::Hexagonal;
|
||||
|
||||
const bool rotation_integration = integrating && experiment.IsRotationIndexing();
|
||||
const CrystalLattice start_lattice = data.latt;
|
||||
const float start_beam_x = data.geom.GetBeamX_pxl(), start_beam_y = data.geom.GetBeamY_pxl();
|
||||
switch (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm()) {
|
||||
case GeomRefinementAlgorithmEnum::None:
|
||||
break;
|
||||
@@ -470,13 +500,38 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In
|
||||
XtalOptimizerRotationOnly(data, msg.spots, 0.1);
|
||||
XtalOptimizerRotationOnly(data, msg.spots, 0.05);
|
||||
break;
|
||||
case GeomRefinementAlgorithmEnum::BeamCenter:
|
||||
if (XtalOptimizer(data, msg.spots)) {
|
||||
case GeomRefinementAlgorithmEnum::BeamCenter: {
|
||||
const XtalOptimizerData start = data;
|
||||
// An integrated rotation frame is refined about the spindle (XtalOptimizerOnSpindle), both fits.
|
||||
const auto fit = [&](XtalOptimizerData &d) {
|
||||
return rotation_integration
|
||||
? XtalOptimizerOnSpindle(d, msg.spots, msg.number, *experiment.GetGoniometer())
|
||||
: XtalOptimizer(d, msg.spots);
|
||||
};
|
||||
bool refined = fit(data);
|
||||
// Rotation, once the run has its prior (SetRefinementPrior): the same fit again from the start,
|
||||
// with a Gaussian restraint on how far the frame moves from the sweep's geometry, weighed against
|
||||
// the frame's own spots at the residual scale the free fit just found. A frame with few or noisy
|
||||
// spots then stays near the sweep and a well determined one follows a crystal that moves
|
||||
// (empirical Bayes: the widths are the run's, the evidence each frame's own).
|
||||
if (refined && rotation_integration && refinement_prior_) {
|
||||
const double sigma_res = std::sqrt(data.fit_residual_var);
|
||||
XtalOptimizerData restrained = start;
|
||||
restrained.prior_lattice = start.latt;
|
||||
restrained.prior_beam_x = start.geom.GetBeamX_pxl();
|
||||
restrained.prior_beam_y = start.geom.GetBeamY_pxl();
|
||||
restrained.orientation_prior_weight = sigma_res / refinement_prior_->orientation_rad;
|
||||
restrained.beam_prior_weight = sigma_res / refinement_prior_->beam_px;
|
||||
refined = fit(restrained);
|
||||
data = refined ? restrained : start;
|
||||
}
|
||||
if (refined) {
|
||||
outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl())
|
||||
.BeamY_pxl(data.geom.GetBeamY_pxl());
|
||||
outcome.beam_center_updated = true;
|
||||
}
|
||||
break;
|
||||
}
|
||||
case GeomRefinementAlgorithmEnum::Flex: {
|
||||
// Try all three refinements per image and keep whichever indexes the most spots. Beam+cell
|
||||
// refinement helps some stills but diverges on sparse spot lists (few spots, long axes),
|
||||
@@ -514,6 +569,25 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In
|
||||
|
||||
outcome.lattice_candidate = data.latt;
|
||||
|
||||
// The cell is held on rotation, so the refined vectors are the sweep's turned by one rotation R,
|
||||
// A_refined = R A_start for the matrix A of the three real-space vectors.
|
||||
if (rotation_integration && collect_frame_refinement_) {
|
||||
const auto matrix = [](const CrystalLattice &l) {
|
||||
Eigen::Matrix3d m;
|
||||
const Coord v[3] = {l.Vec0(), l.Vec1(), l.Vec2()};
|
||||
for (int i = 0; i < 3; i++)
|
||||
m.col(i) << v[i].x, v[i].y, v[i].z;
|
||||
return m;
|
||||
};
|
||||
const Eigen::AngleAxisd aa(Eigen::Matrix3d(matrix(data.latt) * matrix(start_lattice).inverse()));
|
||||
const Eigen::Vector3d w = aa.angle() * aa.axis();
|
||||
const FrameRefinement r{.rotation = Coord(w.x(), w.y(), w.z()),
|
||||
.beam_dx = data.geom.GetBeamX_pxl() - start_beam_x,
|
||||
.beam_dy = data.geom.GetBeamY_pxl() - start_beam_y};
|
||||
const std::unique_lock ul(frame_refinement_mutex_);
|
||||
frame_refinement_.insert_or_assign(msg.number, r);
|
||||
}
|
||||
|
||||
if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Monoclinic)
|
||||
outcome.lattice_candidate->ReorderMonoclinic();
|
||||
|
||||
@@ -873,6 +947,48 @@ void IndexAndRefine::ProbeSupercellFrame(const DataMessage &msg, BraggPrediction
|
||||
supercell_probe_by_frame_.insert_or_assign(msg.number, ProbedFrame{frame, trend, prim});
|
||||
}
|
||||
|
||||
std::optional<RefinementPrior> IndexAndRefine::EstimateRefinementPrior() {
|
||||
const std::unique_lock ul(frame_refinement_mutex_);
|
||||
std::vector<int64_t> frames;
|
||||
std::vector<FrameRefinement> track;
|
||||
for (const auto &[f, r]: frame_refinement_) {
|
||||
frames.push_back(f);
|
||||
track.push_back(r);
|
||||
}
|
||||
if (track.size() < 2 * REFINEMENT_NEIGHBOURS + 1)
|
||||
return std::nullopt;
|
||||
const auto median = [](std::vector<float> &v) {
|
||||
std::ranges::nth_element(v, v.begin() + v.size() / 2);
|
||||
return v[v.size() / 2];
|
||||
};
|
||||
// The track's smooth part: each component's median over the frames within REFINEMENT_NEIGHBOURS
|
||||
// images, the frame itself included. A crystal that moves moves it; a frame whose fit jumped on its
|
||||
// own does not. Its mean square about the sweep is the width the crystal's real motion has.
|
||||
double sum_w = 0.0, sum_b = 0.0;
|
||||
std::vector<float> c[5];
|
||||
for (size_t i = 0; i < track.size(); i++) {
|
||||
for (auto &v: c)
|
||||
v.clear();
|
||||
for (size_t j = i > REFINEMENT_NEIGHBOURS ? i - REFINEMENT_NEIGHBOURS : 0;
|
||||
j < track.size() && frames[j] <= frames[i] + REFINEMENT_NEIGHBOURS; j++) {
|
||||
if (frames[j] < frames[i] - REFINEMENT_NEIGHBOURS)
|
||||
continue;
|
||||
c[0].push_back(track[j].rotation.x);
|
||||
c[1].push_back(track[j].rotation.y);
|
||||
c[2].push_back(track[j].rotation.z);
|
||||
c[3].push_back(track[j].beam_dx);
|
||||
c[4].push_back(track[j].beam_dy);
|
||||
}
|
||||
for (int k = 0; k < 3; k++)
|
||||
sum_w += std::pow(median(c[k]), 2);
|
||||
for (int k = 3; k < 5; k++)
|
||||
sum_b += std::pow(median(c[k]), 2);
|
||||
}
|
||||
const double n = static_cast<double>(track.size());
|
||||
return RefinementPrior{.orientation_rad = std::max(std::sqrt(sum_w / (3.0 * n)), MIN_PRIOR_ORIENTATION_RAD),
|
||||
.beam_px = std::max(std::sqrt(sum_b / (2.0 * n)), MIN_PRIOR_BEAM_PX)};
|
||||
}
|
||||
|
||||
SupercellTrend IndexAndRefine::GetSupercellTrend() {
|
||||
const std::unique_lock ul(supercell_probe_mutex_);
|
||||
SupercellTrend sum{};
|
||||
@@ -885,7 +1001,8 @@ SupercellTrend IndexAndRefine::GetSupercellTrend() {
|
||||
|
||||
std::optional<IndexAndRefine::IndexingOutcome>
|
||||
IndexAndRefine::DetermineRefineAnalyze(DataMessage &msg, const SpotFindingSettings &spot_finding_settings,
|
||||
int64_t *spots_on_lattice, FrameIndexerCalls *frame_calls) {
|
||||
bool integrating, int64_t *spots_on_lattice,
|
||||
FrameIndexerCalls *frame_calls) {
|
||||
if (!indexer_ || !spot_finding_settings.indexing)
|
||||
return std::nullopt;
|
||||
|
||||
@@ -900,7 +1017,7 @@ IndexAndRefine::DetermineRefineAnalyze(DataMessage &msg, const SpotFindingSettin
|
||||
return std::nullopt;
|
||||
|
||||
if (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm() != GeomRefinementAlgorithmEnum::None)
|
||||
RefineGeometryIfNeeded(msg, outcome);
|
||||
RefineGeometryIfNeeded(msg, outcome, integrating);
|
||||
|
||||
if (!outcome.lattice_candidate.has_value())
|
||||
return std::nullopt;
|
||||
@@ -927,7 +1044,7 @@ void IndexAndRefine::ProcessImage(DataMessage &msg,
|
||||
const BraggIntegrateFn &integrate,
|
||||
const BraggIntegrateFn &probe_integrate,
|
||||
FrameIndexerCalls *frame_calls) {
|
||||
auto outcome = DetermineRefineAnalyze(msg, spot_finding_settings, nullptr, frame_calls);
|
||||
auto outcome = DetermineRefineAnalyze(msg, spot_finding_settings, true, nullptr, frame_calls);
|
||||
if (outcome && spot_finding_settings.quick_integration)
|
||||
QuickPredictAndIntegrate(msg, spot_finding_settings, prediction, integrate, *outcome, probe_integrate);
|
||||
}
|
||||
@@ -937,13 +1054,13 @@ bool IndexAndRefine::IndexFrameOnly(DataMessage &msg, const SpotFindingSettings
|
||||
// The per-frame verdict, not "was there anything to integrate": this is what the rotation first pass
|
||||
// scores candidate lattices on, and on rotation a frame can be integrable without indexing on its
|
||||
// own (AnalyzeIndexing), which would score a sparse frame for every candidate alike.
|
||||
return DetermineRefineAnalyze(msg, spot_finding_settings, nullptr, frame_calls).has_value()
|
||||
return DetermineRefineAnalyze(msg, spot_finding_settings, false, nullptr, frame_calls).has_value()
|
||||
&& msg.indexing_result.value_or(false);
|
||||
}
|
||||
|
||||
int64_t IndexAndRefine::CountSpotsOnLattice(DataMessage &msg, const SpotFindingSettings &spot_finding_settings) {
|
||||
int64_t spots_on_lattice = 0;
|
||||
DetermineRefineAnalyze(msg, spot_finding_settings, &spots_on_lattice);
|
||||
DetermineRefineAnalyze(msg, spot_finding_settings, false, &spots_on_lattice);
|
||||
return spots_on_lattice;
|
||||
}
|
||||
|
||||
|
||||
@@ -24,6 +24,7 @@
|
||||
#include "scale_merge/ScaleOnTheFly.h"
|
||||
#include "scale_merge/ScalingResult.h"
|
||||
#include "IntegrationOutcome.h"
|
||||
#include "geom_refinement/XtalOptimizer.h" // RefinementPrior
|
||||
|
||||
// Integrates the predicted reflections off whatever image the caller holds: the preprocessed GPU/CPU
|
||||
// buffer on the WithoutFPGA path (GPU when available), or the assembled detector image read straight,
|
||||
@@ -199,11 +200,26 @@ class IndexAndRefine {
|
||||
// Returns the outcome (ready for integration) when the frame indexes, nullopt otherwise. Both the
|
||||
// real per-image ProcessImage and the first-pass scheme validation go through this, so they cannot
|
||||
// diverge.
|
||||
// `integrating`: the frame is integrated with the outcome (ProcessImage), rather than only asked
|
||||
// whether it indexes - see RefineGeometryIfNeeded.
|
||||
std::optional<IndexingOutcome> DetermineRefineAnalyze(DataMessage &msg,
|
||||
const SpotFindingSettings &spot_finding_settings,
|
||||
bool integrating,
|
||||
int64_t *spots_on_lattice = nullptr,
|
||||
FrameIndexerCalls *frame_calls = nullptr);
|
||||
void RefineGeometryIfNeeded(DataMessage &msg, IndexingOutcome &outcome);
|
||||
void RefineGeometryIfNeeded(DataMessage &msg, IndexingOutcome &outcome, bool integrating);
|
||||
|
||||
// Rotation: what the per-frame refinement did to each integrated frame - the rotation taking the
|
||||
// sweep's lattice to the refined one (angle-axis, radians) and the beam centre's shift - by image
|
||||
// number.
|
||||
struct FrameRefinement {
|
||||
Coord rotation;
|
||||
float beam_dx = 0.0f, beam_dy = 0.0f;
|
||||
};
|
||||
bool collect_frame_refinement_ = false;
|
||||
std::map<int64_t, FrameRefinement> frame_refinement_;
|
||||
std::optional<RefinementPrior> refinement_prior_;
|
||||
std::mutex frame_refinement_mutex_;
|
||||
void QuickPredictAndIntegrate(DataMessage &msg,
|
||||
const SpotFindingSettings &spot_finding_settings,
|
||||
BraggPrediction &prediction,
|
||||
@@ -271,6 +287,20 @@ public:
|
||||
void ProbeSupercell(std::vector<int64_t> image_numbers);
|
||||
[[nodiscard]] SupercellProbe GetSupercellProbe();
|
||||
[[nodiscard]] SupercellTrend GetSupercellTrend();
|
||||
|
||||
// Rotation: the integrated frames (ProcessImage) are refined under this prior; the frames only asked
|
||||
// whether they index keep the free refinement.
|
||||
void SetRefinementPrior(std::optional<RefinementPrior> prior) { refinement_prior_ = prior; }
|
||||
// Rotation: remember what the refinement did to each integrated frame, for EstimateRefinementPrior.
|
||||
void CollectFrameRefinement() { collect_frame_refinement_ = true; }
|
||||
// The prior's widths from the collected frames (see the .cpp); empty for a run with too few
|
||||
// refined frames to say.
|
||||
std::optional<RefinementPrior> EstimateRefinementPrior();
|
||||
static constexpr int REFINEMENT_NEIGHBOURS = 3;
|
||||
// Floors on the widths, far below what a crystal is ever measured to drift: they only keep the
|
||||
// prior's weight finite on a perfectly still track.
|
||||
static constexpr double MIN_PRIOR_ORIENTATION_RAD = 1e-5;
|
||||
static constexpr double MIN_PRIOR_BEAM_PX = 1e-3;
|
||||
// The primitive lattice the parity classes are indexed in (the highest probed frame's; only its metric
|
||||
// is meant to be read), so a class can be turned into the cell it would double to.
|
||||
[[nodiscard]] std::optional<CrystalLattice> GetSupercellProbePrimitive();
|
||||
|
||||
@@ -141,6 +141,44 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
|
||||
break;
|
||||
}
|
||||
|
||||
// The orientation restraint's centre, in the same parameters as latt_vec0.
|
||||
if (data.prior_lattice && data.orientation_prior_weight > 0.0) {
|
||||
double lengths[3], beta_prior;
|
||||
switch (data.crystal_system) {
|
||||
case gemmi::CrystalSystem::Hexagonal:
|
||||
LatticeToRodriguesAndLengths_Hex(*data.prior_lattice, problem.orientation_prior, lengths);
|
||||
break;
|
||||
case gemmi::CrystalSystem::Monoclinic:
|
||||
LatticeToRodriguesLengthsBeta_Mono(*data.prior_lattice, problem.orientation_prior, lengths,
|
||||
beta_prior);
|
||||
break;
|
||||
default:
|
||||
LatticeToRodriguesAndLengths_GS(*data.prior_lattice, problem.orientation_prior, lengths);
|
||||
break;
|
||||
}
|
||||
problem.orientation_prior_weight = data.orientation_prior_weight;
|
||||
// The rotation a change dw of the angle-axis parameters w makes, R(w + dw) = R(J dw) R(w):
|
||||
// J = I + (1 - cos t)/t^2 [w]x + (t - sin t)/t^3 [w]x^2, t = |w|. Along w it is 1, across it
|
||||
// 2 sin(t/2)/t, which is 0.64 for a crystal turned by t = 180 deg - a restraint on dw alone
|
||||
// would be that much tighter across the axis than the width it was given.
|
||||
const Eigen::Vector3d w0(problem.orientation_prior[0], problem.orientation_prior[1],
|
||||
problem.orientation_prior[2]);
|
||||
const double t = w0.norm();
|
||||
Eigen::Matrix3d J = Eigen::Matrix3d::Identity();
|
||||
if (t > 1e-8) {
|
||||
Eigen::Matrix3d W;
|
||||
W << 0, -w0.z(), w0.y(), w0.z(), 0, -w0.x(), -w0.y(), w0.x(), 0;
|
||||
J += (1.0 - std::cos(t)) / (t * t) * W + (t - std::sin(t)) / (t * t * t) * W * W;
|
||||
}
|
||||
for (int i = 0; i < 3; i++)
|
||||
for (int k = 0; k < 3; k++)
|
||||
problem.orientation_prior_jacobian[3 * i + k] = J(i, k);
|
||||
}
|
||||
if (data.refine_beam_center && data.beam_prior_weight > 0.0) {
|
||||
problem.priors.push_back({XtalRefinePrior::Block::Beam, 1.0, 0.0, data.prior_beam_x, data.beam_prior_weight});
|
||||
problem.priors.push_back({XtalRefinePrior::Block::Beam, 0.0, 1.0, data.prior_beam_y, data.beam_prior_weight});
|
||||
}
|
||||
|
||||
// The spindle. `rocking_spindle` is the fallback for a caller that holds one frame and so
|
||||
// passes no `axis` to back-rotate by: the back-rotation is the identity there either way
|
||||
// (angle_rad is zero and an AngleAxisRotator of a zero angle-axis ignores the vector, so the
|
||||
@@ -456,6 +494,7 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
|
||||
// failed refinement leaves data untouched rather than committing half a fit.
|
||||
if (!summary.IsSolutionUsable())
|
||||
return false;
|
||||
data.fit_residual_var = 2.0 * summary.final_cost / (3.0 * static_cast<double>(problem.residuals.size()));
|
||||
|
||||
// A solve can walk the basis flat and still report success - the volume clamp in XtalResidual
|
||||
// keeps every step finite on the way. A coplanar cell has no reciprocal basis, so that is a
|
||||
|
||||
@@ -12,6 +12,13 @@
|
||||
#include "../common/SpotToSave.h"
|
||||
#include "gemmi/symmetry.hpp"
|
||||
|
||||
// A Gaussian prior on how far one rotation frame's refinement moves from the sweep's geometry: the width
|
||||
// for each orientation component (radians) and each beam-centre coordinate (pixels).
|
||||
struct RefinementPrior {
|
||||
double orientation_rad = 0.0;
|
||||
double beam_px = 0.0;
|
||||
};
|
||||
|
||||
struct XtalOptimizerData {
|
||||
DiffractionGeometry geom;
|
||||
CrystalLattice latt;
|
||||
@@ -55,7 +62,19 @@ struct XtalOptimizerData {
|
||||
std::optional<Coord> rocking_spindle;
|
||||
float rocking_wedge_deg = 0.0f;
|
||||
|
||||
// A Gaussian restraint of the orientation toward prior_lattice's and of the beam centre toward
|
||||
// (prior_beam_x, prior_beam_y): the residual weight * (parameter - prior) on each component of the
|
||||
// rotation from the prior orientation (radians) and on each beam coordinate, beside the spots' own.
|
||||
// Zero weight, no restraint.
|
||||
std::optional<CrystalLattice> prior_lattice;
|
||||
double prior_beam_x = 0.0, prior_beam_y = 0.0;
|
||||
double orientation_prior_weight = 0.0;
|
||||
double beam_prior_weight = 0.0;
|
||||
|
||||
// output
|
||||
// The last solve's mean squared (weighted) miss per component. Meant for an unrestrained fit: a
|
||||
// restraint's cost is counted in.
|
||||
double fit_residual_var = 0.0;
|
||||
std::optional<double> beam_corr_x;
|
||||
std::optional<double> beam_corr_y;
|
||||
|
||||
|
||||
@@ -121,7 +121,25 @@ namespace {
|
||||
}
|
||||
|
||||
void AddPriors(const XtalRefineProblem &p, const double *x, double &cost, Eigen::VectorXd *g,
|
||||
Eigen::MatrixXd *H, int beam_tangent, int rot_tangent) {
|
||||
Eigen::MatrixXd *H, int beam_tangent, int rot_tangent, int orientation_tangent) {
|
||||
if (p.orientation_prior_weight > 0.0) {
|
||||
const double w = p.orientation_prior_weight;
|
||||
const double *J = p.orientation_prior_jacobian;
|
||||
double d[3], r[3];
|
||||
for (int k = 0; k < 3; k++)
|
||||
d[k] = x[OFF_P0 + k] - p.orientation_prior[k];
|
||||
for (int i = 0; i < 3; i++) {
|
||||
r[i] = w * (J[3 * i] * d[0] + J[3 * i + 1] * d[1] + J[3 * i + 2] * d[2]);
|
||||
cost += 0.5 * r[i] * r[i];
|
||||
}
|
||||
if (g && orientation_tangent >= 0)
|
||||
for (int k = 0; k < 3; k++)
|
||||
for (int i = 0; i < 3; i++) {
|
||||
(*g)[orientation_tangent + k] += w * J[3 * i + k] * r[i];
|
||||
for (int l = 0; l < 3; l++)
|
||||
(*H)(orientation_tangent + k, orientation_tangent + l) += w * w * J[3 * i + k] * J[3 * i + l];
|
||||
}
|
||||
}
|
||||
for (const auto &prior: p.priors) {
|
||||
const bool on_beam = prior.block == XtalRefinePrior::Block::Beam;
|
||||
const double *v = x + (on_beam ? OFF_BEAM : OFF_ROT);
|
||||
@@ -231,7 +249,7 @@ namespace {
|
||||
cost = s.cost;
|
||||
if (g)
|
||||
ToTangent<LANES>(s, map, *g, *H);
|
||||
AddPriors(p, x, cost, g, H, map[0], map[2]);
|
||||
AddPriors(p, x, cost, g, H, map[0], map[2], map[OBS]);
|
||||
return AllFinite(cost, g, H);
|
||||
}
|
||||
|
||||
@@ -271,7 +289,7 @@ namespace {
|
||||
cost = s.cost;
|
||||
if (g)
|
||||
ToTangent<R_LANES>(s, map, *g, *H);
|
||||
AddPriors(p, x, cost, g, H, map[0], -1);
|
||||
AddPriors(p, x, cost, g, H, map[0], -1, map[R_OBS]);
|
||||
return AllFinite(cost, g, H);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -54,6 +54,13 @@ struct XtalRefineProblem {
|
||||
double latt_vec2_lower[3] = {-kNoBound, -kNoBound, -kNoBound}, latt_vec2_upper[3] = {kNoBound, kNoBound, kNoBound};
|
||||
|
||||
std::vector<XtalRefinePrior> priors;
|
||||
// An isotropic restraint orientation_prior_weight * J (latt_vec0 - orientation_prior) on the rotation
|
||||
// that takes the prior orientation to the refined one; zero weight, none. J (row-major) maps a change
|
||||
// of the angle-axis parameters at orientation_prior to the rotation it makes, so the restraint has
|
||||
// the same width in every direction whatever the crystal's orientation in the laboratory.
|
||||
double orientation_prior[3] = {0, 0, 0};
|
||||
double orientation_prior_jacobian[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1};
|
||||
double orientation_prior_weight = 0.0;
|
||||
|
||||
LMOptions options;
|
||||
};
|
||||
|
||||
@@ -1178,6 +1178,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
// change the processing, so a run with one merges exactly what a run without one does.
|
||||
if (!config_.reference_data.empty())
|
||||
indexer->ReferenceIntensities(config_.reference_data);
|
||||
// The geometry pre-pass refines its frames freely and measures them for the prior every later pass
|
||||
// refines its integrated frames under.
|
||||
if (geometry_prepass)
|
||||
indexer->CollectFrameRefinement();
|
||||
else
|
||||
indexer->SetRefinementPrior(refinement_prior_);
|
||||
// Second pass of the rotation two-pass: predict with the smoothed mosaicity fitted in the pre-pass.
|
||||
if (!geometry_prepass && !prepass_mosaicity_.empty())
|
||||
indexer->SetPredictionMosaicityOverride(prepass_mosaicity_);
|
||||
|
||||
@@ -44,6 +44,7 @@
|
||||
#include "../image_analysis/scale_merge/SearchSpaceGroup.h" // SearchSpaceGroupResult
|
||||
#include "../image_analysis/scale_merge/CrystalSetting.h" // MetricViolation, MAX_METRIC_VIOLATION
|
||||
#include "../image_analysis/geom_refinement/PostRefine.h" // PostRefineResult
|
||||
#include "../image_analysis/geom_refinement/XtalOptimizer.h" // RefinementPrior
|
||||
#include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" // BeamCenterEstimate
|
||||
#include "../image_analysis/geom_refinement/GeometryRefiner.h" // GeomRefineFrame
|
||||
#include "../image_analysis/rotation_indexer/RotationIndexer.h" // RotationIndexerResult
|
||||
@@ -1358,6 +1359,9 @@ class Rugnux {
|
||||
// pre-pass at the file's centre merges the pre-pass subset as well where that arm follows it
|
||||
// (ProcessResult::prepass_subset_merge); the arm itself does not.
|
||||
bool measured_centre_arm_ = false;
|
||||
// Rotation: the prior the passes after the geometry pre-pass refine each integrated frame under,
|
||||
// measured on that pre-pass's own per-frame refinement (IndexAndRefine::EstimateRefinementPrior).
|
||||
std::optional<RefinementPrior> refinement_prior_;
|
||||
// Whether a geometry pre-pass integrates only the pre-pass subset (PrepassSubset): the beam-centre
|
||||
// arm at the measured centre, first.
|
||||
bool prepass_subset_only_ = false;
|
||||
|
||||
@@ -506,6 +506,14 @@ void Rugnux::ImagePass(PipelineLocals &p) {
|
||||
const auto stopped = [&] { return config_.fast_live ? result.cancelled : static_cast<bool>(cancelled_); };
|
||||
Note(fmt::format("{} images, {:.0f} per second", result.images_processed,
|
||||
result.images_processed / std::max(1e-6, result.image_loop_time_s)));
|
||||
if (geometry_prepass && indexer && experiment_.IsRotationIndexing() && !stopped()) {
|
||||
refinement_prior_ = indexer->EstimateRefinementPrior();
|
||||
if (refinement_prior_)
|
||||
logger.Info("Per-frame refinement prior: the pre-pass's refined orientation wanders {:.4f} deg "
|
||||
"and its beam centre {:.3f} px from the sweep's (rms of the {}-frame running median)",
|
||||
refinement_prior_->orientation_rad * 180.0 / M_PI, refinement_prior_->beam_px,
|
||||
2 * IndexAndRefine::REFINEMENT_NEIGHBOURS + 1);
|
||||
}
|
||||
|
||||
// A pre-pass that stopped at loop_end is a sweep of the frames it read, and that is all the merge and
|
||||
// the post-refinement are handed. Left in, the empty tail is a stretch of frames every per-frame
|
||||
|
||||
@@ -386,6 +386,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
|
||||
prepass_rotation_scale_.reset();
|
||||
prepass_result_.reset();
|
||||
force_rotation_result_.reset();
|
||||
refinement_prior_.reset();
|
||||
bragg_adaptive_.reset();
|
||||
beam_center_same_lattice_alternative_.reset();
|
||||
lattice_arm_check_.reset();
|
||||
|
||||
@@ -741,3 +741,94 @@ TEST_CASE("XtalOptimizer_rotation_axis_direction_only") {
|
||||
|
||||
// --- helpers for lattice sanity tests ---
|
||||
#include <Eigen/Dense>
|
||||
|
||||
TEST_CASE("XtalOptimizer_prior_shrinks_toward_its_centre") {
|
||||
DiffractionExperiment exp_i;
|
||||
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(200);
|
||||
|
||||
const CrystalLattice true_latt(40, 50, 80, 90, 90, 90);
|
||||
BraggPredictionSettings prediction_settings{.high_res_A = 1.5, .ewald_dist_cutoff = 0.001};
|
||||
BraggPrediction prediction;
|
||||
const auto count = prediction.Calc(exp_i, true_latt, prediction_settings);
|
||||
std::vector<SpotToSave> spots;
|
||||
for (int i = 0; i < count; ++i) {
|
||||
const auto &refl = prediction.GetReflections().at(i);
|
||||
spots.push_back(SpotToSave{refl.predicted_x, refl.predicted_y});
|
||||
}
|
||||
|
||||
// The centre the refinement starts from and is restrained to: the crystal turned by 0.2 deg and the
|
||||
// beam centre one pixel off, with the data all pointing back at the truth.
|
||||
const CrystalLattice centre = true_latt.Multiply(RotMatrix(0.2f * static_cast<float>(M_PI) / 180.0f, Coord(0, 1, 0)));
|
||||
const auto refined = [&](double orientation_weight, double beam_weight) {
|
||||
XtalOptimizerData data{};
|
||||
data.max_iterations = 50;
|
||||
data.refine_unit_cell = false;
|
||||
data.crystal_system = gemmi::CrystalSystem::Orthorhombic;
|
||||
data.latt = centre;
|
||||
data.geom.BeamX_pxl(1001).BeamY_pxl(1000).DetectorDistance_mm(200);
|
||||
data.prior_lattice = centre;
|
||||
data.prior_beam_x = 1001;
|
||||
data.prior_beam_y = 1000;
|
||||
data.orientation_prior_weight = orientation_weight;
|
||||
data.beam_prior_weight = beam_weight;
|
||||
REQUIRE(XtalOptimizer(data, spots));
|
||||
return data;
|
||||
};
|
||||
const auto moved_deg = [&](const XtalOptimizerData &d) {
|
||||
const auto r = d.latt.Vec0() - centre.Vec0();
|
||||
return r.Length() / centre.Vec0().Length() * 180.0 / M_PI;
|
||||
};
|
||||
|
||||
// No weight: the data alone, which reach the truth.
|
||||
const auto free_fit = refined(0.0, 0.0);
|
||||
CHECK(moved_deg(free_fit) == Catch::Approx(0.2).margin(0.02));
|
||||
CHECK(free_fit.geom.GetBeamX_pxl() == Catch::Approx(1000.0).margin(0.05));
|
||||
|
||||
// A weight far above the data's: the centre.
|
||||
const auto held = refined(1e3, 1e3);
|
||||
CHECK(moved_deg(held) < 0.005);
|
||||
CHECK(held.geom.GetBeamX_pxl() == Catch::Approx(1001.0).margin(0.01));
|
||||
|
||||
// A weight of the data's own order: in between.
|
||||
const auto between = refined(0.3, 3e-3);
|
||||
CHECK(moved_deg(between) > moved_deg(held));
|
||||
CHECK(moved_deg(between) < moved_deg(free_fit));
|
||||
}
|
||||
|
||||
TEST_CASE("XtalOptimizer_prior_same_whatever_the_basis_is_turned_by") {
|
||||
// The same crystal written with a and b negated - a basis turned by 180 deg about c - is the same
|
||||
// physics and must be restrained the same: the restraint acts on the rotation, not on the angle-axis
|
||||
// parameters, which across the axis change by 1/0.64 as much for the same rotation at 180 deg.
|
||||
DiffractionExperiment exp_i;
|
||||
exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(200);
|
||||
|
||||
const CrystalLattice true_latt(40, 50, 80, 90, 90, 90);
|
||||
BraggPredictionSettings prediction_settings{.high_res_A = 1.5, .ewald_dist_cutoff = 0.001};
|
||||
BraggPrediction prediction;
|
||||
const auto count = prediction.Calc(exp_i, true_latt, prediction_settings);
|
||||
std::vector<SpotToSave> spots;
|
||||
for (int i = 0; i < count; ++i) {
|
||||
const auto &refl = prediction.GetReflections().at(i);
|
||||
spots.push_back(SpotToSave{refl.predicted_x, refl.predicted_y});
|
||||
}
|
||||
|
||||
const CrystalLattice centre = true_latt.Multiply(RotMatrix(0.2f * static_cast<float>(M_PI) / 180.0f, Coord(1, 1, 0).Normalize()));
|
||||
const CrystalLattice flipped(-centre.Vec0(), -centre.Vec1(), centre.Vec2());
|
||||
const auto moved_deg = [&](const CrystalLattice &start) {
|
||||
XtalOptimizerData data{};
|
||||
data.max_iterations = 50;
|
||||
data.refine_unit_cell = false;
|
||||
data.refine_beam_center = false;
|
||||
data.crystal_system = gemmi::CrystalSystem::Orthorhombic;
|
||||
data.latt = start;
|
||||
data.geom.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(200);
|
||||
data.prior_lattice = start;
|
||||
data.orientation_prior_weight = 3.0;
|
||||
REQUIRE(XtalOptimizer(data, spots));
|
||||
const auto r = data.latt.Vec2() - centre.Vec2();
|
||||
return r.Length() / centre.Vec2().Length() * 180.0 / M_PI;
|
||||
};
|
||||
const double a = moved_deg(centre), b = moved_deg(flipped);
|
||||
CHECK(a > 0.01);
|
||||
CHECK(b == Catch::Approx(a).epsilon(0.02));
|
||||
}
|
||||
Reference in new issue
Block a user