diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 869b72efa..f2915b50b 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -12,6 +12,7 @@ * `rugnux` measures the beam centre on every run and indexes with it when the file's value indexes nothing, refines only the detector-tilt component the data determine - a beam-centre error is no longer reported as a tilt - and places a detector swung out on a 2theta arm where the file says it stands. * `rugnux` accepts a spot into the per-frame geometry refinement where the rotation the exposure could have supplied accounts for the miss, on data collected at 0.5 degrees per image or coarser. * `rugnux` writes the unmerged MTZ by default, with a P1 merge beside it, a batch header for every image the observations span, and events kept to the same `--min-captured-fraction` as the merge, so a wrong space group can be re-merged in a scaling program without reprocessing. +* `rugnux` decides which partials belong to one rocking event from an angle rather than a frame count, so a coarsely sliced sweep no longer joins two crossings of the Ewald sphere into a single full reflection. * `rugnux` writes reflection files in the conventions downstream programs read: `FreeR_flag` is 0 for the test set and 1 for the working set - it was the other way round - the merged and P1 MTZ carry the reserved `HKL_base` dataset so a CCP4 program reads the wavelength instead of falling back to 1.54187 A, and the merged mmCIF marks the free set as `_refln.status` = `f`. * The rugnux results report carries the space groups the data cannot separate and the enantiomorph state, the model's verdict and what it was allowed to decide, the detector geometry measured and what a single sweep cannot determine, the resolution the CC1/2 fit reached, which reciprocal axis each anisotropic diffraction limit belongs to, and twinning measured before and after the space group was decided; `REPORT_VERSION` is 6, and `SPACE_GROUP_ENANTIOMORPH= DETERMINED_FROM_MODEL` is now `ASSUMED_FROM_MODEL`. * `rugnux --mode calibration` writes `.json` beside the `.poni`, holding the geometry as a `jfjoch_broker` `dataset_settings` body, and refuses a fit that is not a measurement - no `.poni`, a non-zero exit, `converged` recorded in the `.json`; `--no-refine-tilt` holds the detector tilt at the file's value instead of zeroing it. diff --git a/image_analysis/IntegrationOutcome.h b/image_analysis/IntegrationOutcome.h index 75b2a09cd..c0d8855fa 100644 --- a/image_analysis/IntegrationOutcome.h +++ b/image_analysis/IntegrationOutcome.h @@ -3,6 +3,8 @@ #pragma once +#include + #include "../common/Reflection.h" #include "../common/CrystalLattice.h" #include "../common/DiffractionGeometry.h" @@ -16,4 +18,16 @@ struct IntegrationOutcome { std::optional image_scale_cc_n; std::optional image_scale_g; std::optional image_scale_wedge_deg; -}; \ No newline at end of file +}; + +// How far apart, in frames, two partials of one raw hkl may sit and still belong to the same rocking +// event. The quantity being bridged is an ANGLE - a reflecting range, a tenth to half a degree - so +// spelling it as a frame count makes the bridge grow with the slicing: at 1 deg per frame the two +// frames that were always allowed are 2 deg of dead rotation, as wide as a whole event, and two +// genuine Ewald crossings fuse into one "full". Half a degree of bridge instead, floored at one frame +// (never cut an event at its own neighbours) and capped at the two frames that were always used. For +// any wedge of 0.25 deg or less the quotient is at least two and the cap returns the literal 2.0f, so +// finely sliced data is bridged exactly as before. +inline float RockingEventFrameGap(float wedge_deg) { + return std::min(2.0f, std::max(1.0f, 0.5f / wedge_deg)); +} \ No newline at end of file diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 8332288e4..e14cc2a4b 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -598,7 +598,7 @@ float DetectorY(const Reflection &r) { return std::isfinite(r.observed_y) ? r.ob // Sum each rocking event into one full observation. A rotation reflection is integrated image by // image, so it arrives here as a run of partials over consecutive frames; the run is cut where the -// 3D combine cuts it - same raw hkl, frames no further apart than MAX_FRAME_GAP - so the exported +// 3D combine cuts it - same raw hkl, frames no further apart than RockingEventFrameGap - so the exported // file and rugnux's own merge see exactly the same events. // The parts are added, plainly, with their variances in quadrature, which is what every other // rotation program writes as a full. Nothing is divided by the partiality: FRACTIONCALC carries the @@ -630,8 +630,9 @@ float DetectorY(const Reflection &r) { return std::isfinite(r.observed_y) ? r.ob // mosaicity, RotationScaleMerge::SmoothMosaicity), so clamping it to 1 would hide the spread and buy // a reading program nothing. std::vector SumRockingEvents(const std::vector &outcomes, - double min_partiality, double min_captured_fraction) { - constexpr float MAX_FRAME_GAP = 2.0f; // == RotationScaleMerge's: what makes one rocking event + double min_partiality, double min_captured_fraction, + float wedge_deg) { + const float max_frame_gap = RockingEventFrameGap(wedge_deg); // == RotationScaleMerge's // The sort key travels with the part instead of being read back through the pointer, the way the // merge's own ingest sort carries it (RotationScaleMerge's SortKey): there are millions of parts @@ -661,7 +662,7 @@ std::vector SumRockingEvents(const std::vector & size_t j = i + 1; while (j < parts.size() && parts[j].h == parts[i].h && parts[j].k == parts[i].k && parts[j].l == parts[i].l - && parts[j].image_number - parts[j - 1].image_number <= MAX_FRAME_GAP) + && parts[j].image_number - parts[j - 1].image_number <= max_frame_gap) ++j; double sum_p = 0.0, sum_I = 0.0, sum_var = 0.0, sum_var_bkg = 0.0; @@ -815,7 +816,8 @@ void WriteUnmergedMtzReflections(const std::vector &outcomes if (scanning && sum_partials) { for (const auto &r : SumRockingEvents(outcomes, experiment.GetScalingSettings().GetMinPartiality(), - experiment.GetScalingSettings().GetMinCapturedFraction())) + experiment.GetScalingSettings().GetMinCapturedFraction(), + wedge_deg)) add_row(r); } else { for (const auto &outcome : outcomes) diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index a4ba669c7..6ef913795 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -260,11 +260,11 @@ PostRefineResult PostRefineRotationGeometry(const std::vector=2-frame events carry an // unbiased phi_obs (a single-frame centroid is just the frame centre). - constexpr float MAX_FRAME_GAP = 2.0f; + 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) + && pts[j].img - pts[j - 1].img <= max_frame_gap) ++j; return j; }; diff --git a/image_analysis/scale_merge/AnisotropyAnalysis.cpp b/image_analysis/scale_merge/AnisotropyAnalysis.cpp index 893389bf1..c490f11dd 100644 --- a/image_analysis/scale_merge/AnisotropyAnalysis.cpp +++ b/image_analysis/scale_merge/AnisotropyAnalysis.cpp @@ -1029,7 +1029,7 @@ namespace { std::vector ScaledObservations(const std::vector &outcomes, bool rotation, const gemmi::SpaceGroup *space_group, - double min_partiality) { + float wedge_deg, double min_partiality) { // Per-image scale, indexed the way the outcomes are. std::vector g(outcomes.size(), 0.0); for (size_t i = 0; i < outcomes.size(); ++i) @@ -1070,7 +1070,7 @@ std::vector ScaledObservations(const std::vector ScaledObservations(const std::vector ScaledObservations(const std::vector &outcomes, bool rotation, const gemmi::SpaceGroup *space_group = nullptr, + float wedge_deg = 0.0f, double min_partiality = 0.5); // What the caller knows about the run and the merge that the reflections alone do not say. diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 0d879ddd9..f61da9be2 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -91,7 +91,6 @@ namespace { // guarantees the modulation costs data is the CC dip) constexpr double SWEEP_HARMONIC_CC_DIP = 0.8; // ... and cost real signal at its trough - constexpr float MAX_FRAME_GAP = 2.0f; // a rocking event is a run of frames no more apart than this constexpr double CHI2_1_MEDIAN = 0.454936; // A post-scale-fulls correction surface (decay / absorption) is applied only if its held-out // cross-validation improvement exceeds this fraction of the held-out scatter. A margin (not just >0) @@ -228,6 +227,8 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, merge_friedel = s.GetMergeFriedel(); capture_uncertainty_coeff = s.GetCaptureUncertaintyCoeff(); min_captured_fraction = s.GetMinCapturedFraction(); + if (const auto gon = x.GetGoniometer(); gon.has_value()) + max_frame_gap = RockingEventFrameGap(gon->GetWedge_deg()); min_cc_for_image = s.GetMinCCForImage(); search_min_zeta = s.GetSearchMinZeta(); reject_nsigma = s.GetOutlierRejectNsigma(); @@ -874,7 +875,7 @@ void RotationScaleMerge::SmoothMosaicityAndPartiality() { while (i < hi) { int j = i + 1; while (j < hi && partials[perm[j]].image_number - - partials[perm[j - 1]].image_number <= MAX_FRAME_GAP) + - partials[perm[j - 1]].image_number <= max_frame_gap) ++j; if (j - i >= 2) { const float f0 = partials[perm[i]].image_number; @@ -2458,7 +2459,7 @@ void RotationScaleMerge::Combine() { float last_frame = partials[ev[i]].image_number; while (kk < ev.size()) { const float frame = partials[ev[kk]].image_number; - if (frame - last_frame > MAX_FRAME_GAP) break; + if (frame - last_frame > max_frame_gap) break; last_frame = frame; ++kk; } @@ -3850,7 +3851,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st if (gpu_active_ && observation_dump_path.empty()) { // The smoothed corr is already resident (scaling + smooth-G ran on the device, no round-trip). const int nf = gpu_->Combine(rawrun_group.data(), min_partiality, capture_uncertainty_coeff, - min_captured_fraction); + min_captured_fraction, max_frame_gap); g_full.assign(n_frames, 1.0); if (scale_fulls && nf > 0) { diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 3ee86b5ee..13e78a5c9 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -154,6 +154,10 @@ private: bool merge_friedel = true; double capture_uncertainty_coeff = 0.0; double min_captured_fraction = 0.0; + // RockingEventFrameGap of this run's oscillation width, taken once so that the two places that + // cut events on the CPU and the GPU kernel that cuts them on the device all compare the same + // float (RotationScaleMergeGPU.cu documents that bit-parity contract). + float max_frame_gap = 2.0f; // Drop a frame's observations entirely when the frame disagrees with the merged reference below this // correlation (--min-image-cc). A mis-centred or off-crystal frame still produces spots, still diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 8b778ec7d..e6b2821f8 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -275,8 +275,6 @@ namespace { __device__ __forceinline__ double Dmin(double a, double b) { return (b < a) ? b : a; } __device__ __forceinline__ float Fmax(float a, float b) { return (a < b) ? b : a; } - constexpr float COMBINE_MAX_FRAME_GAP = 2.0f; // == RotationScaleMerge::MAX_FRAME_GAP - // A partial is usable for the combine iff its corr and (I, sigma) are finite with corr>0, sigma>0. __device__ __forceinline__ bool CombineUsable(int i, const float *I, const float *sigma, const float *corr) { @@ -289,6 +287,9 @@ namespace { struct CombineParams { int n_runs; double min_partiality, capture_uncertainty_coeff, min_captured_fraction; + // RotationScaleMerge::max_frame_gap, passed in rather than recomputed here so that the host + // and the device compare the same float and the combine stays bit-identical on both paths. + float max_frame_gap; const float *__restrict__ I, *__restrict__ sigma, *__restrict__ corr, *__restrict__ partiality, *__restrict__ bkg, *__restrict__ var_bkg, *__restrict__ image_number, *__restrict__ d, *__restrict__ px, *__restrict__ py; @@ -328,7 +329,7 @@ namespace { while (probe < hi && !CombineUsable(p.perm[probe], p.I, p.sigma, p.corr)) ++probe; if (probe >= hi) break; const float img = p.image_number[p.perm[probe]]; - if (img - last_img > COMBINE_MAX_FRAME_GAP) break; + if (img - last_img > p.max_frame_gap) break; last_img = img; ev_end = probe; ++probe; @@ -1024,7 +1025,8 @@ void RotationScaleMergeGPU::SetRawRuns(int n_runs, int n_perm, const int32_t *pe } int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_partiality, - double capture_uncertainty_coeff, double min_captured_fraction) { + double capture_uncertainty_coeff, double min_captured_fraction, + float max_frame_gap) { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; CudaCheck(cudaMemcpy(d.rr_group.get(), rawrun_group, size_t(d.n_runs) * sizeof(int32_t), @@ -1035,6 +1037,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti p.min_partiality = min_partiality; p.capture_uncertainty_coeff = capture_uncertainty_coeff; p.min_captured_fraction = min_captured_fraction; + p.max_frame_gap = max_frame_gap; p.I = d.I.get(); p.sigma = d.sigma.get(); p.corr = d.corr.get(); p.partiality = d.partiality.get(); p.bkg = d.bkg.get(); p.var_bkg = d.var_bkg.get(); p.image_number = d.image_number.get(); p.d = d.d_obs.get(); p.px = d.px_obs.get(); p.py = d.py_obs.get(); diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index e821ae916..50bf0c23f 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -118,13 +118,13 @@ public: // Combine the resident partials (reading the current resident corr) into fulls on the device, // mirroring RotationScaleMerge::Combine: one thread per raw-hkl run splits its usable partials into - // rocking events (frame gap <= 2), pools background, seeds F, does 3 de-biased Poisson reweights and + // rocking events (frame gap <= max_frame_gap, the host's RockingEventFrameGap), pools background, seeds F, does 3 de-biased Poisson reweights and // adds the capture-uncertainty term. rawrun_group (length n_runs) is the current space group's ASU // id per raw hkl (it becomes the full's group). Deterministic: fulls are emitted in raw-run-major, // event order (a count pass -> host prefix sum -> emit-at-offset), matching the CPU path. Returns the // number of fulls (call GetFulls with buffers of that length). int Combine(const int32_t *rawrun_group, double min_partiality, double capture_uncertainty_coeff, - double min_captured_fraction); + double min_captured_fraction, float max_frame_gap); // Download the combined fulls SoA (length = Combine()'s return). The working corr is downloaded // separately by GetFullsCorr (it is only meaningful after ScaleFulls; otherwise the caller sets it). diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index d2f717c31..72a07c51e 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -4810,7 +4810,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b sm.statistics.anisotropy = AnalyzeAnisotropy( sm.merged, ScaledObservations(indexer->GetIntegrationOutcome(), - experiment_.IsRotationIndexing(), twin_sg), + experiment_.IsRotationIndexing(), twin_sg, + experiment_.GetGoniometer() + ? experiment_.GetGoniometer()->GetWedge_deg() : 0.0f), *result.consensus_cell, twin_sg, aniso_run); stats_text << AnisotropyToText(sm.statistics.anisotropy) << "\n"; } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index d672d9a80..a403ca728 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1770,7 +1770,10 @@ static int RunRugnux(int argc, char **argv) { aniso_run.dose_term_in_scale_model = experiment.GetScalingSettings().GetCorrectionSurfaces(); aniso_run.radiation_damage_relative_b = merged_statistics.radiation_damage_delta_b; merged_statistics.anisotropy = AnalyzeAnisotropy(merged_reflections, - ScaledObservations(reflections, is_rotation, twin_sg), + ScaledObservations(reflections, is_rotation, twin_sg, + experiment.GetGoniometer() + ? experiment.GetGoniometer()->GetWedge_deg() + : 0.0f), *experiment.GetUnitCell(), twin_sg, aniso_run); std::cout << AnisotropyToText(merged_statistics.anisotropy) << std::endl; }