From 76b3956e118ea29dbfe0c5708319ca7a567bf9cd Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 6 Oct 2026 23:10:26 +0200 Subject: [PATCH] French-Wilson once, on the written merge; none for the P1 cross-check Every full-statistics rotation merge computed French-Wilson amplitudes - the search candidates, the pinned low/high pair, the final, and the P1 cross-check - although only the written merge's amplitudes are used, and those were recomputed anyway with the anisotropic prior. RotationScaleMerge gets SetFrenchWilson(); Rugnux turns it off for its merges and makes the amplitudes once after the anisotropy analysis (anisotropic prior when a tensor was fitted, isotropic otherwise). --mode scale keeps the engine's own amplitudes. p_P1.mtz is now intensities only (owner decision): the F/SIGF/F(+)/F(-) columns are omitted; every intensity column is bit-identical to before. French-Wilson was 29% of the CPU samples in the large-cell tail (8tyy P1 cross-check window, ~84 core-s). p.mtz md5 unchanged on myob, cytc, thau, kdp, 8tyy (GPU build). Wall: 8tyy 173.9 -> 150.9 s; cytc/thau -0.2..-0.3 s, myob unchanged (3 interleaved repeats). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- docs/RUGNUX_INTEGRATION.md | 4 +++- docs/RUGNUX_TUTORIAL.md | 5 +++-- image_analysis/WriteReflections.cpp | 21 ++++++++++++------- image_analysis/WriteReflections.h | 4 +++- .../scale_merge/RotationScaleMerge.cpp | 2 +- .../scale_merge/RotationScaleMerge.h | 5 +++++ rugnux/Rugnux.cpp | 15 ++++++++++--- 7 files changed, 40 insertions(+), 16 deletions(-) diff --git a/docs/RUGNUX_INTEGRATION.md b/docs/RUGNUX_INTEGRATION.md index 88b6362d7..294d8439a 100644 --- a/docs/RUGNUX_INTEGRATION.md +++ b/docs/RUGNUX_INTEGRATION.md @@ -62,7 +62,9 @@ H K L IMEAN SIGIMEAN I(+) SIGI(+) I(-) SIGI(-) F SIGF F(+) SIGF(+) F(-) SIGF(-) H H H J Q K M K M F Q G L G L I ``` -`F` is the French–Wilson amplitude. The header carries the determined space group, the refined cell +`F` is the French–Wilson amplitude. `_P1.mtz` carries the intensity columns only (no `F` +columns): it is there to be re-merged or re-searched, and amplitudes are made from the group finally +chosen. The header carries the determined space group, the refined cell and the wavelength, on a dataset of its own behind the reserved `HKL_base`, which is where the MTZ format puts them. Older Rugnux wrote the data on dataset 0, the id reserved for `HKL_base`, and CCP4's `mtzinfo` then reported its 1.54187 Å (Cu Kα) default instead of the real wavelength — every diff --git a/docs/RUGNUX_TUTORIAL.md b/docs/RUGNUX_TUTORIAL.md index ce6e29912..04d73c34d 100644 --- a/docs/RUGNUX_TUTORIAL.md +++ b/docs/RUGNUX_TUTORIAL.md @@ -277,8 +277,9 @@ Absorption is corrected only by the empirical surface a rotation run fits and ke `_unmerged_partials.mtz`, one row per image, instead of summing. - `_P1.mtz` — the **P1 cross-check dataset**: the same observations merged in P1 instead of the space group the run determined, so a wrong call can be re-merged, re-solved or re-refined without - processing the images again. It is a full merge, not the degraded one the search itself runs on, and - it is what `rugnux --mode scale -S P1` would make from a `_process.h5`. **Every rotation run that + processing the images again. It holds intensities only, without French–Wilson amplitudes. It is a + full merge, not the degraded one the search itself runs on, and + its intensities are what `rugnux --mode scale -S P1` would make from a `_process.h5`. **Every rotation run that determines its own space group writes it** — including one that determined P1, where it simply repeats the merged output — so a script harvesting results can always expect the file rather than having to reproduce the search's decision to know whether it exists. It is **not** written when `-S` diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 013b74f15..049ac5e5a 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -460,7 +460,8 @@ void WriteMmcifReflections(const std::vector &reflections, void WriteMtzReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, - const std::string &filename) { + const std::string &filename, + bool amplitudes) { gemmi::Mtz mtz; // Optional but recommended metadata @@ -495,9 +496,11 @@ void WriteMtzReflections(const std::vector &reflections, mtz.add_column("I(-)", 'K', dataset_id, -1, false); mtz.add_column("SIGI(-)", 'M', dataset_id, -1, false); } - mtz.add_column("F", 'F', dataset_id, -1, false); // French-Wilson amplitude - mtz.add_column("SIGF", 'Q', dataset_id, -1, false); - if (has_anom) { + if (amplitudes) { + mtz.add_column("F", 'F', dataset_id, -1, false); // French-Wilson amplitude + mtz.add_column("SIGF", 'Q', dataset_id, -1, false); + } + if (amplitudes && has_anom) { mtz.add_column("F(+)", 'G', dataset_id, -1, false); mtz.add_column("SIGF(+)", 'L', dataset_id, -1, false); mtz.add_column("F(-)", 'G', dataset_id, -1, false); @@ -506,7 +509,7 @@ void WriteMtzReflections(const std::vector &reflections, mtz.add_column("FreeR_flag", 'I', dataset_id, -1, false); mtz.nreflections = static_cast(out_rows.size()); - mtz.data.reserve(out_rows.size() * (has_anom ? 16 : 8)); + mtz.data.reserve(out_rows.size() * 16); for (const auto& row : out_rows) { mtz.data.push_back(static_cast(row.h)); mtz.data.push_back(static_cast(row.k)); @@ -519,9 +522,11 @@ void WriteMtzReflections(const std::vector &reflections, mtz.data.push_back(row.Im); mtz.data.push_back(row.sIm); } - mtz.data.push_back(row.Fmean); - mtz.data.push_back(row.sFmean); - if (has_anom) { + if (amplitudes) { + mtz.data.push_back(row.Fmean); + mtz.data.push_back(row.sFmean); + } + if (amplitudes && has_anom) { mtz.data.push_back(row.Fp); mtz.data.push_back(row.sFp); mtz.data.push_back(row.Fm); diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 697cbe87c..0b6187eb9 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -37,10 +37,12 @@ void WriteMmcifReflections(const std::vector &reflections, const std::string &filename, size_t nthreads); +// amplitudes: write the French-Wilson F columns; off for a merge that has none (the P1 cross-check). void WriteMtzReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, - const std::string &filename); + const std::string &filename, + bool amplitudes = true); // SHELX HKLF-4 text file, FORMAT(3I4,2F8.2) h k l I sigma(I), scaled so the largest value fits the // field, ended by a 0 0 0 record. nthreads: workers for the row formatting, as for the mmCIF. diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index ff59811dc..b0a3dde74 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -5814,7 +5814,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool } // French-Wilson amplitudes for IMEAN and (now that they are attached) each Bijvoet hand. - if (full_stats) { + if (full_stats && french_wilson) { FrenchWilsonOptions fw_opts; fw_opts.num_threads = static_cast(nthreads); ApplyFrenchWilson(result.merged, Group(), fw_opts); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index e4a312a91..a5ea82261 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -147,6 +147,10 @@ public: // it is a copy as long as the fulls, wanted only by a caller that writes them. void SetExportScaledFulls(bool on) { export_scaled_fulls = on; } + // Whether a full-statistics Run() gives the merge French-Wilson amplitudes (on by default). Rugnux + // turns it off and makes the amplitudes once, on the merge it writes, with the anisotropic prior. + void SetFrenchWilson(bool on) { french_wilson = on; } + // Where the engine logs from now on (by default the logger it was made with) - a merge run beside // other work logs into a held buffer, so the run's log keeps one order. void SetLogger(Logger &to) { logger = &to; } @@ -467,6 +471,7 @@ private: #endif bool write_back_per_frame_scale = true; // see SetWriteBackPerFrameScale bool export_scaled_fulls = false; // see SetExportScaledFulls + bool french_wilson = true; // see SetFrenchWilson // --- helpers (each a flat pass; see the .cpp) --- // Turn the per-frame mean background under the reflections (accumulated by the ingest fill loop) into diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index a533a4904..43442378c 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -7055,6 +7055,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // fulls, and the search merges, the pre-pass and the P1 cross-check never write a .hkl. rsm->SetExportScaledFulls(!for_search && !geometry_prepass && write_files && config_.write_merged && export_scaled_fulls); + // Amplitudes are made once, on the merge that is written, where the anisotropy tensor + // their prior uses is known (see below) - not on every merge the tail makes on its way. + rsm->SetFrenchWilson(false); auto r = already_run ? std::move(*already_run) : rsm->Run(for_search, /*full_stats=*/!geometry_prepass, /*measure_cc_before_corrections=*/measure_cc); @@ -8872,6 +8875,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (p1_crosscheck && rsm) { rsm->SetWriteBackPerFrameScale(false); rsm->SetExportScaledFulls(false); + rsm->SetFrenchWilson(false); rsm->SetLogger(p1_log); p1_ahead = std::async(std::launch::async, [&rsm] { return rsm->Run(/*for_search=*/false, /*full_stats=*/true, @@ -9142,14 +9146,19 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (anisotropy.valid()) { sm.statistics.anisotropy = anisotropy.get(); stats_text << AnisotropyToText(sm.statistics.anisotropy) << "\n"; - // The amplitudes are the one place the tensor is used: French-Wilson again, with the - // Wilson prior of each reflection taken along its own direction. The intensities and + // The amplitudes are the one place the tensor is used: French-Wilson with the Wilson + // prior of each reflection taken along its own direction. The intensities and // everything measured on them above are left as they are. FrenchWilsonOptions fw_opts; fw_opts.num_threads = config_.nthreads; fw_opts.anisotropy_hkl = AnisotropyTensorHKL(sm.statistics.anisotropy, gemmi::UnitCell(*result.consensus_cell)); ApplyFrenchWilson(sm.merged, experiment_.GetSpaceGroupOrP1(), fw_opts); + } else if (rsm) { + // No tensor: the isotropic prior. The rotation merge leaves its amplitudes to here. + FrenchWilsonOptions fw_opts; + fw_opts.num_threads = config_.nthreads; + ApplyFrenchWilson(sm.merged, experiment_.GetSpaceGroupOrP1(), fw_opts); } } @@ -9636,7 +9645,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const std::string path = config_.output_prefix + "_P1.mtz"; WriteMtzReflections(p1.merged, - *result.consensus_cell, experiment_, path); + *result.consensus_cell, experiment_, path, /*amplitudes=*/false); experiment_.SetSpaceGroup(determined_group); result.error_model_isa = em_isa; result.error_model_isa_asymptotic = em_isa_asymptotic;