From f2b92e3f4dc22f6b179923a2686cf70d88f1d8b4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 28 Jul 2026 20:03:43 +0200 Subject: [PATCH] rugnux: --search-min-zeta drops badly-measured observations from the symmetry search zeta is the sine of the angle between a reflection's rocking path and the spindle. Near 0 the reflection crosses the Ewald sphere almost tangentially, spends many frames in diffracting position and is measured worst. The de-novo space-group search asks how EQUAL an operator's paired intensities are, so its answer is dominated by whichever reflections are measured worst - and when the spindle lies in a lattice plane, an operator that permutes the two in-plane axes samples a different mixture of measurement qualities than one that only flips signs. That is not a fair comparison, and it can make a real symmetry operator look like a twin law. Measured on a thaumatin set mounted that way (its 4-fold is 88.9 deg from the spindle), the added operators' disagreement is 1.74x the parent's over pairs where both reflections have zeta < 0.85 and 1.003x - i.e. the symmetry is exact - over pairs where both are above it. The search consequently refuses the 422 promotion and merges the crystal in P222, while the same data forced to the right group give CC1/2 99.2% at multiplicity 10.7, matching XDS. With the option the de-novo pass ignores those observations (the final merge keeps everything - there completeness is the point): zeta cut observations ignored H ratio adopted 0 (off) - 1.47 P222 0.5 1620648 1.44 P222 0.7 3006013 1.34 P21212 0.85 4536724 promoted P4212 (correct point group) OFF BY DEFAULT, and it must stay off, because the same cut costs four other crystals their space group (P41212 -> P212121, I23 -> P2, I23 -> I222 twice): at 0.85 it discards 40-80% of all observations, which on a crystal whose geometry is not the problem simply starves the search. Two independent implementations - filtering the pairs that enter the statistic, and filtering the observations that enter the merge - trade exactly the same crystals, so this is a property of the cut and not of where it is applied. Verified bit-identical to the previous binary when off. The companion diagnostic is already there: the run now reports how close a symmetry axis lies to the spindle, which is the geometry that makes this option worth reaching for. Implementation note for anyone tempted by the cheaper route: excluding these observations from the ASU grouping alone does NOT work. The 3D combine selects partials on corr, not on their group, so their intensity still reaches the fulls and the merged intensities are unchanged - measured, the statistic did not move by 0.03 while 67% of observations were nominally excluded. Zeroing corr is what removes an observation from the combine, the merge and the error model alike. Co-Authored-By: Claude Opus 5 (1M context) --- common/ScalingSettings.cpp | 12 +++++++++ common/ScalingSettings.h | 6 +++++ .../scale_merge/RotationScaleMerge.cpp | 25 +++++++++++++++++++ .../scale_merge/RotationScaleMerge.h | 10 ++++++++ rugnux/rugnux_cli.cpp | 8 ++++++ 5 files changed, 61 insertions(+) diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index 9a43f393..617ff2ac 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -89,6 +89,18 @@ double ScalingSettings::GetMinCCForImage() const { return min_cc_for_image; } +double ScalingSettings::GetSearchMinZeta() const { + return search_min_zeta; +} + +ScalingSettings &ScalingSettings::SearchMinZeta(double input) { + if (input < 0.0 || input >= 1.0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Search zeta limit must be in [0,1)"); + search_min_zeta = input; + return *this; +} + ScalingSettings &ScalingSettings::MinCCForImage(double input) { if (input < 0.0 || input > 1.0) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Min CC for image must be between 0 and 1"); diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index abef9dfa..60b98b0d 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -28,6 +28,10 @@ class ScalingSettings { // min_partiality, which gates individual partials; this gates the assembled full. 0 = off (baseline). double min_captured_fraction = 0.0; double min_cc_for_image = 0.0; + + // Exclude observations whose Lorentz geometry |zeta| falls below this from the DE-NOVO space-group + // search merge only (see RotationScaleMerge::search_min_zeta). 0 = off. + double search_min_zeta = 0.0; double outlier_reject_nsigma = 0.0; // per-observation merge outlier rejection (XDS/DIALS-style); 0 = off, e.g. 6 enables // Scale fulls: after the rotation 3D combine, refit a per-frame scale on the combined fulls (XDS @@ -92,6 +96,7 @@ public: ScalingSettings& CaptureUncertaintyCoeff(double input); ScalingSettings& MinCapturedFraction(double input); ScalingSettings& MinCCForImage(double min_cc_for_image); + ScalingSettings& SearchMinZeta(double search_min_zeta); ScalingSettings& OutlierRejectNsigma(double input); ScalingSettings& ScaleFulls(bool input); ScalingSettings& AbsorptionIter(int input); @@ -127,6 +132,7 @@ public: [[nodiscard]] double GetCaptureUncertaintyCoeff() const; [[nodiscard]] double GetMinCapturedFraction() const; [[nodiscard]] double GetMinCCForImage() const; + [[nodiscard]] double GetSearchMinZeta() const; [[nodiscard]] double GetOutlierRejectNsigma() const; [[nodiscard]] bool GetScaleFulls() const; [[nodiscard]] int GetAbsorptionIter() const; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index c8e3f4ae..9b698dbe 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -177,6 +177,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, capture_uncertainty_coeff = s.GetCaptureUncertaintyCoeff(); min_captured_fraction = s.GetMinCapturedFraction(); min_cc_for_image = s.GetMinCCForImage(); + search_min_zeta = s.GetSearchMinZeta(); reject_nsigma = s.GetOutlierRejectNsigma(); reject_outliers = reject_nsigma > 0.0; rfree_fraction = s.GetRfreeFraction(); @@ -1995,6 +1996,30 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, } FinalizePerFrameScale(cc, cc_n, partial_scaled); + // --- 2a. On the de-novo search pass only, drop observations whose Lorentz geometry is poor. + // Zeroing corr is what removes an observation everywhere: the 3D combine, the merge and the + // error model all require corr > 0. Excluding them from the ASU grouping alone is NOT + // enough - the combine selects on corr, so their intensity would still reach the fulls. --- + if (for_search && search_min_zeta > 0.0) { + int64_t n_dropped = 0; + for (auto &o : partials) + if (!(std::isfinite(o.zeta) && o.zeta >= search_min_zeta)) { + if (std::isfinite(o.corr) && o.corr > 0.0f) ++n_dropped; + o.corr = 0.0f; + } +#ifdef JFJOCH_USE_CUDA + if (gpu_active_ && n_dropped > 0) { + std::vector corr(partials.size()); + for (size_t i = 0; i < partials.size(); ++i) corr[i] = partials[i].corr; + gpu_->SetCorr(corr.data()); + } +#endif + if (n_dropped > 0) + logger.Info("Space-group search: ignoring {} observations with |zeta| < {:.2f} " + "(they cross the Ewald sphere near-tangentially and are measured worst)", + n_dropped, search_min_zeta); + } + // --- 2b. Drop frames that do not agree with the merged reference (--min-image-cc). --- if (min_cc_for_image > 0.0) { std::vector reject(n_frames, 0); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index adb37fa9..f93ec27c 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -106,6 +106,16 @@ private: // indexes and still integrates - it just measures something that is not the crystal's diffraction, // and nothing downstream removes it. 0 = off. double min_cc_for_image = 0.0; + + // Exclude observations with |zeta| below this from the DE-NOVO SEARCH merge only (the final merge + // keeps everything). zeta is the sine of the angle between a reflection's rocking path and the + // spindle: near 0 it crosses the Ewald sphere almost tangentially, spends many frames in + // diffracting position and is measured badly. The symmetry search compares how equal an operator's + // paired intensities are, so it is answered by whichever reflections are worst measured - and when + // the spindle lies in a lattice plane, an operator permuting the two in-plane axes samples a + // different mixture of qualities than one that only flips signs, which is not a fair comparison. + // Unlike a bound on I/sigma this is pure geometry, identical in meaning on every dataset. 0 = off. + double search_min_zeta = 0.0; double reject_nsigma = 0.0; bool reject_outliers = false; double rfree_fraction = 0.0; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index cde3d25e..30388b68 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -109,6 +109,7 @@ void print_usage() { std::cout << " --min-captured-fraction rot3d: drop a combined full whose rocking curve was captured below this fraction (edge-of-sweep truncated fulls) (default: 0.7 for rotation, 0 otherwise; 0 = off)" << std::endl; std::cout << " --mosaicity Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl; std::cout << " --reject-outliers Per-observation merge outlier rejection, N sigma from the per-reflection median (default: 6 for rot3d, XDS/DIALS-style; 0 = off)" << std::endl; + std::cout << " --search-min-zeta De-novo space-group search only: ignore observations whose Lorentz geometry |zeta| is below this (0-1). Reflections crossing the Ewald sphere near-tangentially are measured badly and can make a symmetry operator look wrong (default: 0 = use all)" << std::endl; std::cout << " --min-image-cc Per-image CC limit in percent (default: no limit)" << std::endl; std::cout << " --scaling-iterations Number of scaling iterations with no reference data (default: 3)" << std::endl; std::cout << " -z, --reference-mtz Reference MTZ file" << std::endl; @@ -152,6 +153,7 @@ enum { OPT_MAX_SPOTS, OPT_MIN_PARTIALITY, OPT_MIN_IMAGE_CC, + OPT_SEARCH_MIN_ZETA, OPT_SCALING_ITERATIONS, OPT_SCALING_HIGH_RESOLUTION, OPT_RESOLUTION_CUTOFF, @@ -263,6 +265,7 @@ static option long_options[] = { {"min-captured-fraction", required_argument, nullptr, OPT_MIN_CAPTURED_FRACTION}, {"mosaicity", required_argument, nullptr, OPT_MOSAICITY}, {"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC}, + {"search-min-zeta", required_argument, nullptr, OPT_SEARCH_MIN_ZETA}, {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS}, {"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION}, {"resolution-cutoff", required_argument, nullptr, OPT_RESOLUTION_CUTOFF}, @@ -527,6 +530,7 @@ int main(int argc, char **argv) { std::optional capture_uncertainty_arg; // explicit --capture-uncertainty; default depends on rot3d std::optional forced_mosaicity_arg; // diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed double min_image_cc = 0.0; + double search_min_zeta = 0.0; // --search-min-zeta; de-novo search merge only int64_t scaling_iter = 3; std::optional forced_rotation_lattice; std::optional refine_geometry; // --refine-geometry[=N]: stills global geometry-refinement pass @@ -851,6 +855,9 @@ int main(int argc, char **argv) { case OPT_MIN_IMAGE_CC: min_image_cc = parse_double_arg(optarg, "--min-image-cc", logger); break; + case OPT_SEARCH_MIN_ZETA: + search_min_zeta = parse_double_arg(optarg, "--search-min-zeta", logger); + break; case OPT_SCALING_HIGH_RESOLUTION: d_min_scale_merge = atof(optarg); break; @@ -1425,6 +1432,7 @@ int main(int argc, char **argv) { scaling_settings.CaptureUncertaintyCoeff(capture_uncertainty_arg.value_or(rotation_indexing ? 1.0 : 0.0)); scaling_settings.ForcedMosaicity(forced_mosaicity_arg); scaling_settings.MinCCForImage(min_image_cc / 100.0); // --min-image-cc is in percent; the setting is a fraction + scaling_settings.SearchMinZeta(search_min_zeta); scaling_settings.OutlierRejectNsigma( outlier_reject_nsigma.value_or(rotation_indexing ? REJECT_OUTLIERS_DEFAULT_NSIGMA : 0.0));