diff --git a/common/ResolutionShells.h b/common/ResolutionShells.h index 4e0dc4445..739553bdb 100644 --- a/common/ResolutionShells.h +++ b/common/ResolutionShells.h @@ -12,7 +12,7 @@ class ResolutionShells { const float one_over_dmin2, one_over_dmax2; const int32_t nshells; public: - // Shells of equal volume in 1/d^2 between d_min (exclusive) and d_max (inclusive). d_max may be + // Shells of equal width in 1/d^2 between d_min (exclusive) and d_max (inclusive). d_max may be // infinite, which is "no low-resolution bound": the shells then start at 1/d^2 = 0. ResolutionShells(float d_min, float d_max, int32_t nshells); [[nodiscard]] std::optional GetShell(float d) const; diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index 7712f640b..bbb21f1fe 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -281,7 +281,7 @@ The error model is refined as $\sigma_\mathrm{corr}^2 = a\,\sigma^2 + (b\,\langl $a$ and $b$ are **reported in XDS's convention**, which is $\sigma^2 = a(\sigma_0^2 + b I^2)$ with $\mathrm{ISa}=1/\sqrt{ab}$, so the printed pair can be read straight against a `CORRECT.LP`. The internal fit keeps the form above; only the report converts, as $b_\mathrm{XDS} = b^2/a$. Note that $a$ is the same in both conventions and that the two ISa expressions are the same number, $1/\sqrt{a\cdot b^2/a} = 1/b$ — so the rotation log prints **two** ISa, the whole-range $1/b$ (XDS's meaning) and the strong-reflection asymptote beside it, which can only ever be the more optimistic of the two. The mmCIF follows the same split: `_reflns.jfjoch_diffrn_ISa` is the whole-range value, directly comparable with a `CORRECT.LP`, and the asymptote is written separately as `_reflns.jfjoch_diffrn_ISa_asymptotic`, with `_reflns.jfjoch_error_model_a` and `_b` alongside so the number can be re-derived. Note that a file written before this change carries the *asymptote* under the plain `ISa` name. A third, unrelated $b$ appears in the space-group search (§13.1); it is fitted with the $\sigma^2$ coefficient held at 1 and its gate constants are calibrated in that convention. -**Anomalous signal-to-noise (SigAno).** The strength of the anomalous signal is reported per shell and overall as $\mathrm{SigAno}=\langle|\Delta I|\rangle / \langle\sigma(\Delta I)\rangle$, where $\Delta I = I(+)-I(-)$ over acentric reflections measured in both Bijvoet hands and $\sigma(\Delta I)=\sqrt{\sigma_+^2+\sigma_-^2}$. It is computed from the **full-multiplicity** inverse-variance $I(+)/I(-)$ split (the same one written to the output), i.e. from all observations rather than a half-set. For pure noise $\mathrm{SigAno}$ approaches the half-normal value $\sqrt{2/\pi}\approx0.8$, and it rises above $1$ once a real anomalous difference is present. A half-set anomalous correlation ("$\mathrm{CC}_\mathrm{anom}$") is **not** reported. Its two half estimates $\Delta I_0,\Delta I_1$ are complementary partitions of one observation pool ($\Delta I_0+\Delta I_1=2\,\Delta I_\mathrm{full}$), and subtracting the two Bijvoet hands cancels the large common intensity that keeps $\mathrm{CC}_{1/2}$ non-negative, leaving the small anomalous signal against the per-half split noise; below an anomalous signal-to-noise of $1$ per half that correlation tends towards $-1$ rather than $0$. $\mathrm{SigAno}$ has no such floor. It is emitted only when an anomalous split was made, using the standard PDBx items `_reflns.pdbx_absDiff_over_sigma_anomalous` (overall) and `_reflns_shell.pdbx_absDiff_over_sigma_anomalous` (per shell), and appears as the `SigAno` column of the printed merge-statistics table. +**Anomalous signal-to-noise (SigAno).** The strength of the anomalous signal is reported per shell and overall as $\mathrm{SigAno}=\langle|\Delta I|\rangle / \langle\sigma(\Delta I)\rangle$, where $\Delta I = I(+)-I(-)$ over acentric reflections measured in both Bijvoet hands and $\sigma(\Delta I)=\sqrt{\sigma_+^2+\sigma_-^2}$. It is computed from the **full-multiplicity** inverse-variance $I(+)/I(-)$ split (the same one written to the output), i.e. from all observations rather than a half-set. For pure noise $\mathrm{SigAno}$ approaches the half-normal value $\sqrt{2/\pi}\approx0.8$, and it rises above $1$ once a real anomalous difference is present. A half-set anomalous correlation ($\mathrm{CC}_\mathrm{anom}$) is reported beside it: $\Delta I$ is formed once per half-set and the two are correlated over the acentric pairs where **both** hands split into two non-empty halves, per shell and overall as one correlation rather than a mean of shells. Unlike SigAno it is not a ratio against the error model, so an optimistic $\sigma$ cannot inflate it. It has no floor either: subtracting the two Bijvoet hands cancels the large common intensity that keeps $\mathrm{CC}_{1/2}$ non-negative, so on data with little anomalous signal and about two observations per mate it goes strongly negative. That is a property of the statistic and is reported as measured. It agrees with AIMLESS's `CCanom` and `phenix.merging_statistics`' `cc_anom`; XDS's `Anomal Corr` is a different quantity and is not comparable with it. It is emitted only when an anomalous split was made, using the standard PDBx items `_reflns.pdbx_absDiff_over_sigma_anomalous` (overall) and `_reflns_shell.pdbx_absDiff_over_sigma_anomalous` (per shell), and appears as the `SigAno` column of the printed merge-statistics table. ### 10.6 Rotation datasets: combining partials into fulls (3D integration) diff --git a/image_analysis/geom_refinement/PowderCalibration.cpp b/image_analysis/geom_refinement/PowderCalibration.cpp index 97f5d4782..bf6738b40 100644 --- a/image_analysis/geom_refinement/PowderCalibration.cpp +++ b/image_analysis/geom_refinement/PowderCalibration.cpp @@ -424,6 +424,17 @@ void WritePoniFile(const std::string &path, const DiffractionExperiment &experim if (!f) throw JFJochException(JFJochExceptionCategory::FileWriteError, "Cannot write " + path); + // A PONI states the detector in five numbers - two offsets, a distance and three rotations - and + // has no field for a mirrored or quarter-turned IMAGE. The fits are orientation-aware, so writing + // one here would state a geometry that is not the one that was measured, silently and in the file + // a downstream program trusts. Refusing is the only honest option; the JSON beside it carries the + // full geometry. + if (!geom.GetOrientation().IsIdentity()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "A PONI file cannot state this detector's image orientation (a mirror " + "or a quarter turn), so writing one would mis-state the geometry that " + "was fitted; the .json beside it carries the full geometry"); + const double pixel_m = geom.GetPixelSize_mm() * 1e-3; // pyFAI's axis convention is the trap: Poni1 (and pixel1) is the SLOW axis - rows, our y - and // Poni2 the FAST axis - columns, our x - both in metres from the detector origin. A transposed PONI diff --git a/image_analysis/geom_refinement/RingOptimizer.cpp b/image_analysis/geom_refinement/RingOptimizer.cpp index 8aadcdcd3..6df7cdf61 100644 --- a/image_analysis/geom_refinement/RingOptimizer.cpp +++ b/image_analysis/geom_refinement/RingOptimizer.cpp @@ -148,6 +148,13 @@ void ComputeUncertainty(ceres::Problem &problem, const ceres::Solver::Summary &s DiffractionGeometry RingOptimizer::Run(const std::vector &input, RingFitUncertainty *unc) { + // Cleared on the way in. Run writes *unc only where the solve is usable and the covariance + // conditions, so on any other path it would keep whatever the PREVIOUS fit left there - and + // CalibrateFromSpots runs this twice, once with the tilt free and again with it pinned. A failed + // second fit would then hand back the declined first fit's sigmas, valid flag and all. + if (unc) + *unc = RingFitUncertainty{}; + // Initial guess for the parameters double center_x = reference.GetBeamX_pxl(); double center_y = reference.GetBeamY_pxl(); diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 1b8061f83..5e2a8ac7a 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -755,10 +755,12 @@ namespace { } std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in) { + // NaN, not 0.0, where the denominator is empty: Measured() below turns that into a dash, and a + // printed 0.0% reads as a measured total failure rather than as a quantity nobody could form. double completeness = in.possible_unique_reflections > 0 - ? static_cast(in.unique_reflections) / in.possible_unique_reflections * 100.0 : 0.0; + ? static_cast(in.unique_reflections) / in.possible_unique_reflections * 100.0 : NAN; double multiplicity = in.unique_reflections > 0 - ? static_cast(in.total_observations) / in.unique_reflections : 0.0; + ? static_cast(in.total_observations) / in.unique_reflections : NAN; output << fmt::format("{:8d} {:8d} {:8d} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", in.total_observations, diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index d586b39ff..8b45d5640 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -229,9 +229,13 @@ public: // Outlier rejection (driven by ScalingSettings::GetOutlierRejectNsigma) reports its count. [[nodiscard]] size_t RejectedCount() const { return reject_count; } - // image_id is the image's stable identity (its index in the outcomes vector). The CC1/2 half-set - // is a deterministic hash of it, so the split is reproducible run-to-run and independent of the - // order (or threading) of AddImage calls - not a draw from a shared RNG in call order. + // image_id is the image's stable identity (its index in the outcomes vector). It seeds the CC1/2 + // half-set PREFERENCE, which the balancing rule in AddImage then overrides where a reflection's + // two halves would come out uneven - so the half an observation lands in is a function of the + // order AddImage is called in. CALL IT SERIALLY, in image order: the mutex makes concurrent calls + // safe, not reproducible. The rotation path ranks instead (RotationScaleMerge::AssignHalvesByRank) + // because it holds every observation at once; a streaming accumulator cannot, and there is no + // second, device-side walk of the stills observations for this to disagree with. void AddImage(const IntegrationOutcome& outcome, int64_t image_id); // d_min_override, when set, is the effective high-resolution limit for the shell table (used for diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 219b46aa7..0382fdd9e 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -451,6 +451,15 @@ struct SearchSpaceGroupOptions { // 4_1/4_3 conditions of a crystal that has no screw at all, whose predicted-absent class is as // strong as its control row) read -3.8 to -38.5. The bound sits in that gap - loose enough that // three well-measured dead axial reflections clear it, tight enough that two do not. + // + // That last sentence is no longer true, and the reason is worth recording rather than re-tuning + // blind. MIN_U_PER_ABSENT_REFLECTION floors sum_u at 1e-3 per absent reflection, which caps a + // zone's evidence at -a*log(a*1e-3) + log(a!): 13.1 nats at two absences, 19.2 at three, 25.3 at + // four. A zone of three therefore CANNOT clear 20 however dead its reflections are, so the gate + // is in practice "four absences in one zone" - a count, which is the measure the whole statistic + // exists to avoid. A thin monoclinic sweep reaching only three 0k0-odd reflections loses its + // 2(1) to this. Lowering the floor to 7.7e-4 would let three through; that is a recalibration + // and needs its own battery. double min_screw_absence_evidence = 20.0; // Workers for the operator-correlation stage, which is the bulk of the search: one pass over the diff --git a/reader/MiniCBF.cpp b/reader/MiniCBF.cpp index 4cc2f836c..33ad22213 100644 --- a/reader/MiniCBF.cpp +++ b/reader/MiniCBF.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -19,11 +20,23 @@ namespace minicbf { namespace { // One capture group, first match, or nothing. The headers are a few kB, so a regex per field is -// cheap and keeps each rule next to the thing it reads. +// cheap and keeps each rule next to the thing it reads - but COMPILING one is not: measured at ~60 us +// here, and a header parse runs 26 of them. Parsing happens once per image read across a whole sweep, +// so the same two dozen patterns were being recompiled thousands of times. They are compiled once and +// kept; std::map never invalidates a reference, so the pointer outlives the lock. std::optional Match(const std::string &text, const char *pattern) { + static std::mutex cache_mutex; + static std::map cache; + const std::regex *re; + { + const std::lock_guard lock(cache_mutex); + auto it = cache.find(pattern); + if (it == cache.end()) + it = cache.emplace(pattern, std::regex(pattern)).first; + re = &it->second; + } std::smatch m; - const std::regex re(pattern); - if (!std::regex_search(text, m, re) || m.size() < 2) + if (!std::regex_search(text, m, *re) || m.size() < 2) return {}; return m[1].str(); } diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 55f9b3a20..78065bf30 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -354,19 +354,28 @@ std::string RenderResultReport(const std::string &output_prefix, Key(os, "FRIEDELS_LAW", experiment.GetScalingSettings().GetMergeFriedel() ? "TRUE" : "FALSE"); Key(os, "UNIQUE_REFLECTIONS", o.unique_reflections); Key(os, "TOTAL_OBSERVATIONS", o.total_observations); - Key(os, "COMPLETENESS", o.possible_unique_reflections > 0 - ? fmt::format("{:.1f}", 100.0 * o.unique_reflections / o.possible_unique_reflections) - : std::string("nan")); - Key(os, "MULTIPLICITY", o.unique_reflections > 0 - ? fmt::format("{:.2f}", static_cast(o.total_observations) / o.unique_reflections) - : std::string("nan")); - Key(os, "I_OVER_SIGMA", fmt::format("{:.2f}", o.mean_i_over_sigma)); - Key(os, "R_MEAS", fmt::format("{:.4f}", o.r_meas)); - Key(os, "CC_HALF", fmt::format("{:.4f}", o.cc_half)); - Key(os, "SIGANO", fmt::format("{:.3f}", o.abs_diff_over_sigma_anomalous)); + // One rule for every quantity here: a run that did not measure it writes NO key, rather than + // the word "nan" or a zero that reads as a measured absence. SIGANO is the common case - a + // Friedel-merged run splits no Bijvoet pair - and it sat one line from CC_ANOM, which already + // did this, reporting the same missing quantity two different ways. + if (o.possible_unique_reflections > 0) + Key(os, "COMPLETENESS", + fmt::format("{:.1f}", 100.0 * o.unique_reflections / o.possible_unique_reflections)); + if (o.unique_reflections > 0) + Key(os, "MULTIPLICITY", + fmt::format("{:.2f}", static_cast(o.total_observations) / o.unique_reflections)); + if (std::isfinite(o.mean_i_over_sigma)) + Key(os, "I_OVER_SIGMA", fmt::format("{:.2f}", o.mean_i_over_sigma)); + if (std::isfinite(o.r_meas)) + Key(os, "R_MEAS", fmt::format("{:.4f}", o.r_meas)); + if (std::isfinite(o.cc_half)) + Key(os, "CC_HALF", fmt::format("{:.4f}", o.cc_half)); + if (std::isfinite(o.abs_diff_over_sigma_anomalous)) + Key(os, "SIGANO", fmt::format("{:.3f}", o.abs_diff_over_sigma_anomalous)); if (std::isfinite(o.cc_anom)) Key(os, "CC_ANOM", fmt::format("{:.4f}", o.cc_anom)); - Key(os, "WILSON_B", fmt::format("{:.2f}", result.merge_statistics.wilson_b)); + if (std::isfinite(result.merge_statistics.wilson_b)) + Key(os, "WILSON_B", fmt::format("{:.2f}", result.merge_statistics.wilson_b)); // The error model in XDS's convention, so the numbers are directly comparable with a CORRECT.LP. Key(os, "ERROR_MODEL_A", fmt::format("{:.4f}", result.error_model_a)); Key(os, "ERROR_MODEL_B", fmt::format("{:.4e}", result.error_model_b)); @@ -521,8 +530,15 @@ std::string RenderResultReport(const std::string &output_prefix, // run applied is isotropic, so on an anisotropic crystal the outer shells are complete in count // and empty in signal along the weak directions; this says how far the crystal actually reaches // where it reaches furthest. - Key(os, "ANISOTROPY_D_MIN_BEST", - fmt::format("{:.2f}", *std::min_element(an.d_min_axis, an.d_min_axis + 3))); + // Skipping the non-finite entries, as the warning block below does: a cone too sparse to + // cross the threshold leaves its limit NaN, and every comparison against NaN is false, so a + // bare min_element returns element 0 and prints nan even where the other two are measured. + const double *best_axis = nullptr; + for (const double *p = an.d_min_axis; p != an.d_min_axis + 3; ++p) + if (std::isfinite(*p) && (!best_axis || *p < *best_axis)) + best_axis = p; + if (best_axis) + Key(os, "ANISOTROPY_D_MIN_BEST", fmt::format("{:.2f}", *best_axis)); Key(os, "ANISOTROPY_SHAPE", AnisotropyShapeCode(an.shape)); Key(os, "ANISOTROPY_SHAPE_INTERCEPT", fmt::format("{:.3f}", an.shape_intercept)); Key(os, "ANISOTROPY_SHAPE_INTERCEPT_Z", fmt::format("{:.1f}", an.shape_intercept_z)); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index d2a13ff61..a75631ee9 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -4918,27 +4918,37 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b "absences were never predicted, so a P1 merge of this run would be " "missing whole centring classes. Re-run without -S to get one."); if (config_.write_p1_crosscheck && p1_integration_complete && is_rotation) { - const auto determined_sg = experiment_.GetSpaceGroupNumber(); - const gemmi::SpaceGroup *determined = gemmi::find_spacegroup_by_number( - static_cast(determined_sg.value_or(1))); - // scale_and_merge reports the error model of whatever it merged; this merge is not - // the run's answer, so put the determined group's model back afterwards. + // The whole GROUP, not its number: the search can adopt a non-reference setting, and + // a number only ever names the reference one - so restoring through the number would + // silently re-seat P 1 1 2(1) as P 1 2(1) 1 for everything written after this block, + // the unmerged MTZ included, leaving two files that describe one dataset in two + // different settings. + const auto determined_group = experiment_.GetGemmiSpaceGroup(); + const gemmi::SpaceGroup *determined = + determined_group.has_value() ? &*determined_group : nullptr; + // scale_and_merge reports the error model AND the CC1/2 resolution fit of whatever it + // merged; this merge is not the run's answer, so both go back afterwards. The fit in + // particular is what the report calls the number to quote, and a P1 merge of the same + // observations has n_ops times the unique reflections at a fraction of the + // multiplicity, so its fall-off is not this crystal's. const double em_isa = result.error_model_isa; const double em_isa_asymptotic = result.error_model_isa_asymptotic; const double em_a = result.error_model_a; const double em_b = result.error_model_b; + const auto res_fit = result.resolution_fit_A; // Both the merge and the MTZ read the group from the experiment, so it is set for // the whole of it and restored after. experiment_.SpaceGroupNumber(1); const auto p1 = scale_and_merge("P1 cross-check", false); const std::string path = config_.output_prefix + "_P1.mtz"; WriteMtzReflections(p1.merged, *result.consensus_cell, experiment_, path); - experiment_.SpaceGroupNumber(determined_sg); + experiment_.SetSpaceGroup(determined_group); result.error_model_isa = em_isa; result.error_model_isa_asymptotic = em_isa_asymptotic; result.error_model_a = em_a; result.error_model_b = em_b; - if (determined_sg.value_or(1) > 1) + result.resolution_fit_A = res_fit; + if (determined != nullptr && determined->number > 1) logger.Info("P1 cross-check dataset written to {} ({} unique reflections): the " "same observations merged in P1 instead of {}, so a wrong space group " "can be recovered from without reprocessing. It is not the result of " diff --git a/rugnux/WriteModel.h b/rugnux/WriteModel.h index 3600daabc..356a47952 100644 --- a/rugnux/WriteModel.h +++ b/rugnux/WriteModel.h @@ -13,7 +13,7 @@ class Logger; // Write `placed` - the atomic model as ValidateAgainstModel left it, re-fractionalized into the data -// cell and moved by the rigid-body step - to `output_prefix`_model.cif. Everything the input carried +// cell and moved by the rigid-body step - to `output_prefix`_model.cif, and to `_model.pdb` as well where a PDB can express the cell. Everything the input carried // comes with it: chains, residues, ligands, waters, B-factors, occupancies, anisotropic Us. // // `cell` and `space_group` must be the ones the reflection files beside it are written in, i.e. the diff --git a/tests/SearchSpaceGroupTest.cpp b/tests/SearchSpaceGroupTest.cpp index 9863be0d1..d04324061 100644 --- a/tests/SearchSpaceGroupTest.cpp +++ b/tests/SearchSpaceGroupTest.cpp @@ -240,7 +240,7 @@ TEST_CASE("SearchSpaceGroup finds a screw axis whose absent class is weak only w // itself and symmetry cannot fill it in. Counting the class therefore measures the geometry of the // sweep, not the strength of the evidence, and a count gate refused a monoclinic crystal its 2_1 for // having six 0k0-odd reflections rather than eight, every one of them measured at a thousandth of the -// row beside them. The class is judged by ScrewAbsenceEvidence instead, which reads the contrast +// row beside them. The class is judged by ScrewZoneEvidence instead, which reads the contrast // against the row - so few-but-decisive is accepted and many-but-marginal is not. TEST_CASE("SearchSpaceGroup weighs a screw's absences by evidence, not by how many were recorded") { const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("P 1 21 1");