From 87af0f419bd167bece314ea2c9b15e722cb8c7eb Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 12 Jul 2026 21:12:26 +0200 Subject: [PATCH] scale_merge: compute French-Wilson amplitudes at end of merge, write them alongside intensities Merged reflections carried only intensities; a naive sqrt(max(I,0)) turns every weak or negative measurement into a biased (or zero) amplitude. Add a French-Wilson step so the merged reflections themselves carry proper Bayesian amplitudes. New FrenchWilson.{h,cpp}: ApplyFrenchWilson() fills F and sigmaF on each MergedReflection with the posterior mean |F| given I and its sigma under the Wilson prior: - correct centric vs acentric prior (gemmi is_reflection_centric); - epsilon (symmetry-enhancement) multiplicity: the shell mean is and the prior mean for a reflection is eps * ; - numerically stable log-shift integration of the posterior; - strong reflections (I > 4 sigma) short-circuit to sqrt(I) where the FW bias is negligible; - unusable I/sigma falls back to sqrt(max(I,0)). It runs as the last step of the merge routine (both MergeOnTheFly and RotationScaleMerge), so F/sigmaF are part of the merged result and downstream consumers (file writing AND, later, map calculation) all use the same amplitudes rather than each recomputing FW and risking divergence. The writers just emit the fields; extended additively: - MTZ: F (F) + SIGF (Q) columns; - mmCIF: _refln.F_meas_au / _refln.F_meas_sigma_au; - HKL text: F sigmaF appended (rfree flag kept in place for back-compat). Existing intensity columns are unchanged, so current readers keep working while phenix.refine (etc.) can now refine against the amplitudes. Verified on lyso (1.3 A): all 173 negative-intensity reflections get F>0, strong reflections have F/sqrt(I)=1.000, no NaN, and the mmCIF _refln loop is well-formed. Co-Authored-By: Claude Opus 4.8 (1M context) --- common/Reflection.h | 2 + image_analysis/WriteReflections.cpp | 13 +- image_analysis/scale_merge/CMakeLists.txt | 2 + image_analysis/scale_merge/FrenchWilson.cpp | 131 ++++++++++++++++++ image_analysis/scale_merge/FrenchWilson.h | 24 ++++ image_analysis/scale_merge/Merge.cpp | 2 + .../scale_merge/RotationScaleMerge.cpp | 2 + 7 files changed, 174 insertions(+), 2 deletions(-) create mode 100644 image_analysis/scale_merge/FrenchWilson.cpp create mode 100644 image_analysis/scale_merge/FrenchWilson.h diff --git a/common/Reflection.h b/common/Reflection.h index 4f7b4251..bf78ee49 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -43,6 +43,8 @@ struct MergedReflection { float sigma_half[2] = {NAN, NAN}; float d = 0.0; bool rfree_flag = false; + float F = NAN; // French-Wilson amplitude |F| (filled by ApplyFrenchWilson at end of merge) + float sigmaF = NAN; // its sigma }; diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index c2e10d7d..b01bc0a7 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -181,6 +181,8 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_refln.index_l\n"; out << "_refln.intensity_meas\n"; out << "_refln.intensity_sigma\n"; + out << "_refln.F_meas_au\n"; + out << "_refln.F_meas_sigma_au\n"; out << "_refln.status_free\n"; out << "_refln.status\n"; @@ -190,6 +192,8 @@ void WriteMmcifReflections(const std::vector &reflections, << std::setw(5) << r.l << " " << std::setw(14) << Fmt(r.I, 4) << " " << std::setw(14) << Fmt(r.sigma, 4) << " " + << std::setw(14) << Fmt(r.F, 4) << " " + << std::setw(14) << Fmt(r.sigmaF, 4) << " " << (r.rfree_flag ? 1 : 0) << " " << "o" // 'o' = observed << "\n"; @@ -224,13 +228,15 @@ void WriteMtzReflections(const std::vector &reflections, mtz.add_column("L", 'H', dataset_id, -1, false); mtz.add_column("IMEAN", 'J', dataset_id, -1, false); mtz.add_column("SIGIMEAN", 'Q', 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); mtz.add_column("FreeR_flag", 'I', dataset_id, -1, false); // Number of rows mtz.nreflections = reflections.size(); // Allocate data table - mtz.data.reserve(reflections.size() * 6); + mtz.data.reserve(reflections.size() * 8); for (const auto& r : reflections) { mtz.data.push_back(static_cast(r.h)); @@ -238,6 +244,8 @@ void WriteMtzReflections(const std::vector &reflections, mtz.data.push_back(static_cast(r.l)); mtz.data.push_back(r.I); mtz.data.push_back(r.sigma); + mtz.data.push_back(r.F); + mtz.data.push_back(r.sigmaF); mtz.data.push_back(r.rfree_flag ? 1 : 0); } @@ -254,7 +262,8 @@ void WriteHKLReflections(const std::vector &reflections, for (const auto &r: reflections) hkl_file << r.h << " " << r.k << " " << r.l << " " << r.I << " " << r.sigma << " " - << (r.rfree_flag ? 1 : 0) << std::endl; + << (r.rfree_flag ? 1 : 0) << " " + << r.F << " " << r.sigmaF << std::endl; hkl_file.close(); } diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index ea3e655a..9eb12b27 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -15,6 +15,8 @@ ADD_LIBRARY(JFJochScaleMerge HKLKey.h RfreeFlags.cpp RfreeFlags.h + FrenchWilson.cpp + FrenchWilson.h ScalingResult.h ScalingResult.cpp) TARGET_LINK_LIBRARIES(JFJochScaleMerge Ceres::ceres Eigen3::Eigen JFJochCommon) diff --git a/image_analysis/scale_merge/FrenchWilson.cpp b/image_analysis/scale_merge/FrenchWilson.cpp new file mode 100644 index 00000000..0f3a5d8c --- /dev/null +++ b/image_analysis/scale_merge/FrenchWilson.cpp @@ -0,0 +1,131 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "FrenchWilson.h" + +#include +#include +#include +#include + +#include "../../common/ResolutionShells.h" +#include "gemmi/symmetry.hpp" + +namespace { + +struct Posterior { + double mean_I; // (posterior mean true intensity) + double mean_F; // <|F|> (posterior mean amplitude) +}; + +// Posterior moments of the true intensity J >= 0 given a measurement I +/- sigma and the Wilson +// prior with mean sigma_wilson. Integrated numerically over J in [0, I + 8 sigma] with a log-shift +// so the exponentials never overflow/underflow. acentric: p(J) ~ exp(-J/S); centric: +// p(J) ~ exp(-J/2S)/sqrt(J). +Posterior integrate_posterior(double I, double sigma, double sigma_wilson, bool centric, int npts) { + const double inv_2s2 = 1.0 / (2.0 * sigma * sigma); + const double j_max = std::max(I, 0.0) + 8.0 * sigma; + const double dj = j_max / npts; + + std::vector logw(npts); + double max_logw = -std::numeric_limits::infinity(); + for (int i = 0; i < npts; ++i) { + const double j = (i + 0.5) * dj; + const double diff = I - j; + const double log_prior = centric ? (-j / (2.0 * sigma_wilson) - 0.5 * std::log(j)) + : (-j / sigma_wilson); + logw[i] = log_prior - diff * diff * inv_2s2; + max_logw = std::max(max_logw, logw[i]); + } + + double sum_w = 0, sum_wI = 0, sum_wF = 0; + for (int i = 0; i < npts; ++i) { + const double j = (i + 0.5) * dj; + const double w = std::exp(logw[i] - max_logw); + if (!std::isfinite(w)) + continue; + sum_w += w; + sum_wI += w * j; + sum_wF += w * std::sqrt(j); + } + if (sum_w <= 0.0) { + const double j = std::max(I, 0.0); + return {j, std::sqrt(j)}; + } + return {sum_wI / sum_w, sum_wF / sum_w}; +} + +} // namespace + +void ApplyFrenchWilson(std::vector &merged, int32_t space_group_number, + const FrenchWilsonOptions &opts) { + auto naive = [](MergedReflection &r) { + const double ip = std::max(r.I, 0.0f); + r.F = static_cast(std::sqrt(ip)); + r.sigmaF = (ip > 0.0 && std::isfinite(r.sigma)) ? static_cast(r.sigma / (2.0 * std::sqrt(ip))) + : NAN; + }; + + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); + if (sg == nullptr || merged.empty()) { + for (auto &r : merged) naive(r); + return; + } + const gemmi::GroupOps gops = sg->operations(); + + float d_min = std::numeric_limits::max(), d_max = 0.0f; + for (const auto &r : merged) + if (std::isfinite(r.d) && r.d > 0.0f) { + d_min = std::min(d_min, r.d); + d_max = std::max(d_max, r.d); + } + if (!(d_min < d_max && d_min > 0.0f)) { + for (auto &r : merged) naive(r); + return; + } + + // Wilson mean intensity per resolution shell. + ResolutionShells shells(d_min * 0.999f, d_max * 1.001f, opts.num_shells); + std::vector shell_sum(opts.num_shells, 0.0); + std::vector shell_count(opts.num_shells, 0); + double global_sum = 0.0; + int global_count = 0; + auto epsilon = [&](const MergedReflection &r) { + return std::max(1, gops.epsilon_factor_without_centering({{r.h, r.k, r.l}})); + }; + for (const auto &r : merged) { + if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f) + continue; + const double i_over_eps = r.I / epsilon(r); + global_sum += i_over_eps; + ++global_count; + if (const auto s = shells.GetShell(r.d)) { + shell_sum[*s] += i_over_eps; + ++shell_count[*s]; + } + } + const double global_mean = global_count > 0 ? std::max(global_sum / global_count, 1e-10) : 1.0; + std::vector shell_mean(opts.num_shells, global_mean); + for (int s = 0; s < opts.num_shells; ++s) + if (shell_count[s] >= opts.min_reflections_per_shell) + shell_mean[s] = std::max(shell_sum[s] / shell_count[s], 1e-10); + + for (auto &r : merged) { + if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f) { + naive(r); + continue; + } + // Strong reflections: the FW correction is negligible, <|F|> = sqrt(I). + if (r.I > opts.strong_cutoff * r.sigma) { + naive(r); + continue; + } + const auto s = shells.GetShell(r.d); + const double sigma_wilson = epsilon(r) * (s ? shell_mean[*s] : global_mean); + const bool centric = gops.is_reflection_centric({{r.h, r.k, r.l}}); + const Posterior post = integrate_posterior(r.I, r.sigma, sigma_wilson, !centric, + opts.integration_points); + r.F = static_cast(post.mean_F); + r.sigmaF = static_cast(std::sqrt(std::max(0.0, post.mean_I - post.mean_F * post.mean_F))); + } +} diff --git a/image_analysis/scale_merge/FrenchWilson.h b/image_analysis/scale_merge/FrenchWilson.h new file mode 100644 index 00000000..9ea38350 --- /dev/null +++ b/image_analysis/scale_merge/FrenchWilson.h @@ -0,0 +1,24 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "../../common/Reflection.h" + +struct FrenchWilsonOptions { + int num_shells = 20; // resolution shells for the Wilson mean + int min_reflections_per_shell = 20; // below this, fall back to the overall mean for the shell + int integration_points = 400; // posterior integration steps (weak reflections only) + double strong_cutoff = 4.0; // I/sigma above which <|F|> = sqrt(I) (FW bias negligible) +}; + +// Fill F and sigmaF on every merged reflection with the French-Wilson estimate of the amplitude: +// the posterior mean |F| given the measured intensity I and its sigma under the Wilson prior. The +// prior uses the resolution-shell mean intensity, the correct centric/acentric form, and the +// reflection's epsilon (symmetry-enhancement) multiplicity. Strong reflections reduce to sqrt(I); +// reflections with an unusable I/sigma fall back to sqrt(max(I,0)) with propagated sigma. +void ApplyFrenchWilson(std::vector &merged, int32_t space_group_number, + const FrenchWilsonOptions &opts = {}); diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 9bfa9b62..405a9b26 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -18,6 +18,7 @@ #include "../../common/Definitions.h" #include "HKLKey.h" #include "RfreeFlags.h" +#include "FrenchWilson.h" namespace { // Deterministic CC1/2 half-set assignment: a splitmix64 bit-mix of the image's stable index. @@ -472,6 +473,7 @@ std::vector MergeOnTheFly::ExportReflections() { } AssignRfreeFlags(out, space_group_number, scaling_settings.GetRfreeFraction()); + ApplyFrenchWilson(out, space_group_number); return out; } diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 71ec34d9..af6db35d 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -17,6 +17,7 @@ #include "HKLKey.h" #include "RfreeFlags.h" +#include "FrenchWilson.h" #include "ResolutionCutoff.h" #include "../../common/CorrelationCoefficient.h" #include "../../common/CrystalLattice.h" @@ -1300,6 +1301,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool result.merged, d_min_limit, resolution_cutoff_method, resolution_cc_target, for_search, logger); AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction); + ApplyFrenchWilson(result.merged, x.GetSpaceGroupNumber().value_or(1)); if (reject_count > 0) logger.Info("Merge outlier rejection: dropped {} observations", reject_count);