mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases 67 Wiki Activity

rugnux: write I(+)/I(-) and F(+)/F(-) by default (rotation)

Browse Source 
A rotation merge now keeps each acentric reflection's Bijvoet split - the
inverse-variance I(+)/I(-) from the same scaled fulls - even when the merge is
Friedel-averaged (the default, no -A). IMEAN stays the primary intensity and is
bit-identical to before; the split is purely additive and scaled
non-anomalously, so a weak anomalous signal is preserved in the output without
having to reprocess with -A. French-Wilson now also fills F(+)/F(-) from the two
hands (one pass, shared Wilson prior).

Writers: the default (Friedel-merged) MTZ and mmCIF now carry IMEAN plus
I(+)/I(-) and F/F(+)/F(-) whenever a reflection has an anomalous split (the
stills path, which computes none, keeps the previous columns). A missing mate or
a centric is emitted as the CCP4 missing-value flag (NaN).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-17 11:52:31 +02:00
co-authored by Claude Opus 4.8
parent b945444372
commit b647cf454c
4 changed files with 146 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files
common/Reflection.h
+12
View File
@@ -44,6 +44,18 @@ struct MergedReflection {
bool rfree_flag = false;
float F = NAN; // French-Wilson amplitude |F| (filled by ApplyFrenchWilson at end of merge)
float sigmaF = NAN; // its sigma
// Anomalous (Bijvoet) split of this reflection's own observations, kept even when the merge is
// Friedel-averaged (I above is the Friedel mean). Lets I(+)/I(-) be written and CCano reported by
// default without scaling anomalously; NaN when a hand was not measured or for centrics.
float I_plus = NAN;
float sigma_plus = NAN;
float I_minus = NAN;
float sigma_minus = NAN;
// French-Wilson amplitudes of the two hands (filled by ApplyFrenchWilson from I_plus/I_minus).
float F_plus = NAN;
float sigmaF_plus = NAN;
float F_minus = NAN;
float sigmaF_minus = NAN;
};
image_analysis/WriteReflections.cpp
+49 -1
View File
@@ -6,6 +6,7 @@
#include "scale_merge/HKLKey.h"
#include "scale_merge/TwinningAnalysis.h"
#include <algorithm>
#include <cmath>
#include <map>
#include <tuple>
@@ -188,8 +189,16 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
out << "_refln.index_l\n";
out << "_refln.intensity_meas\n";
out << "_refln.intensity_sigma\n";
out << "_refln.pdbx_I_plus\n";
out << "_refln.pdbx_I_plus_sigma\n";
out << "_refln.pdbx_I_minus\n";
out << "_refln.pdbx_I_minus_sigma\n";
out << "_refln.F_meas_au\n";
out << "_refln.F_meas_sigma_au\n";
out << "_refln.pdbx_F_plus\n";
out << "_refln.pdbx_F_plus_sigma\n";
out << "_refln.pdbx_F_minus\n";
out << "_refln.pdbx_F_minus_sigma\n";
out << "_refln.status_free\n";
out << "_refln.status\n";
@@ -199,8 +208,16 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &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.I_plus, 4) << " "
<< std::setw(14) << Fmt(r.sigma_plus, 4) << " "
<< std::setw(14) << Fmt(r.I_minus, 4) << " "
<< std::setw(14) << Fmt(r.sigma_minus, 4) << " "
<< std::setw(14) << Fmt(r.F, 4) << " "
<< std::setw(14) << Fmt(r.sigmaF, 4) << " "
<< std::setw(14) << Fmt(r.F_plus, 4) << " "
<< std::setw(14) << Fmt(r.sigmaF_plus, 4) << " "
<< std::setw(14) << Fmt(r.F_minus, 4) << " "
<< std::setw(14) << Fmt(r.sigmaF_minus, 4) << " "
<< (r.rfree_flag ? 1 : 0) << " "
<< "o" // 'o' = observed
<< "\n";
@@ -235,25 +252,56 @@ void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
// reflection with the standard CCP4 anomalous layout (IMEAN + I(+)/I(-), and the same split for
// the French-Wilson amplitude), which aimless / ctruncate / mtz2sca / ANODE read directly.
if (experiment.GetScalingSettings().GetMergeFriedel()) {
// Friedel-merged output: IMEAN is the primary intensity. I(+)/I(-) are written by default too
// when the merge kept a Bijvoet split (rotation always does; the split is scaled non-anomalously),
// so a weak anomalous signal is preserved without reprocessing. A reflection with only one mate,
// or a centric, gets a missing value (NaN) for the absent hand. When no reflection has an
// anomalous split (e.g. the stills path, which does not compute one) the columns are omitted.
const bool has_anom = std::any_of(reflections.begin(), reflections.end(),
[](const MergedReflection& r){ return std::isfinite(r.I_plus) || std::isfinite(r.I_minus); });
mtz.add_column("H", 'H', dataset_id, -1, false);
mtz.add_column("K", 'H', dataset_id, -1, false);
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);
if (has_anom) {
mtz.add_column("I(+)", 'K', dataset_id, -1, false);
mtz.add_column("SIGI(+)", 'M', dataset_id, -1, false);
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) {
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);
mtz.add_column("SIGF(-)", 'L', dataset_id, -1, false);
}
mtz.add_column("FreeR_flag", 'I', dataset_id, -1, false);
mtz.nreflections = static_cast<int>(reflections.size());
mtz.data.reserve(reflections.size() * 8);
mtz.data.reserve(reflections.size() * (has_anom ? 16 : 8));
for (const auto& r : reflections) {
mtz.data.push_back(static_cast<float>(r.h));
mtz.data.push_back(static_cast<float>(r.k));
mtz.data.push_back(static_cast<float>(r.l));
mtz.data.push_back(r.I);
mtz.data.push_back(r.sigma);
if (has_anom) {
mtz.data.push_back(r.I_plus);
mtz.data.push_back(r.sigma_plus);
mtz.data.push_back(r.I_minus);
mtz.data.push_back(r.sigma_minus);
}
mtz.data.push_back(r.F);
mtz.data.push_back(r.sigmaF);
if (has_anom) {
mtz.data.push_back(r.F_plus);
mtz.data.push_back(r.sigmaF_plus);
mtz.data.push_back(r.F_minus);
mtz.data.push_back(r.sigmaF_minus);
}
mtz.data.push_back(r.rfree_flag ? 1.0f : 0.0f);
}
mtz.write_to_file(filename);
image_analysis/scale_merge/FrenchWilson.cpp
+28 -20
View File
@@ -59,16 +59,24 @@ Posterior integrate_posterior(double I, double sigma, double sigma_wilson, bool
void ApplyFrenchWilson(std::vector<MergedReflection> &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<float>(std::sqrt(ip));
r.sigmaF = (ip > 0.0 && std::isfinite(r.sigma)) ? static_cast<float>(r.sigma / (2.0 * std::sqrt(ip)))
: NAN;
// Naive amplitude sqrt(max(I,0)) for a missing / strong / untrusted intensity; NaN in -> NaN out
// (a missing Bijvoet hand stays missing). Fills one (F, sigmaF) pair.
auto naive_one = [](float I, float sigma, float &F, float &sigF) {
if (!std::isfinite(I)) { F = NAN; sigF = NAN; return; }
const double ip = std::max(I, 0.0f);
F = static_cast<float>(std::sqrt(ip));
sigF = (ip > 0.0 && std::isfinite(sigma)) ? static_cast<float>(sigma / (2.0 * std::sqrt(ip))) : NAN;
};
// The mean intensity and each measured hand share the reflection's Wilson prior, so fill all three.
auto naive_all = [&](MergedReflection &r) {
naive_one(r.I, r.sigma, r.F, r.sigmaF);
naive_one(r.I_plus, r.sigma_plus, r.F_plus, r.sigmaF_plus);
naive_one(r.I_minus, r.sigma_minus, r.F_minus, r.sigmaF_minus);
};
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number);
if (sg == nullptr || merged.empty()) {
for (auto &r : merged) naive(r);
for (auto &r : merged) naive_all(r);
return;
}
const gemmi::GroupOps gops = sg->operations();
@@ -80,7 +88,7 @@ void ApplyFrenchWilson(std::vector<MergedReflection> &merged, int32_t space_grou
d_max = std::max(d_max, r.d);
}
if (!(d_min < d_max && d_min > 0.0f)) {
for (auto &r : merged) naive(r);
for (auto &r : merged) naive_all(r);
return;
}
@@ -110,22 +118,22 @@ void ApplyFrenchWilson(std::vector<MergedReflection> &merged, int32_t space_grou
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;
}
// French-Wilson |F| for one intensity of reflection r (its mean, or one Bijvoet hand); the shell
// Wilson prior, epsilon and centric flag are the reflection's, shared by all three.
auto fw_one = [&](const MergedReflection &r, float I, float sigma, float &F, float &sigF) {
if (!std::isfinite(I) || !std::isfinite(sigma) || sigma <= 0.0f) { naive_one(I, sigma, F, sigF); return; }
// Strong reflections: the FW correction is negligible, <|F|> = sqrt(I).
if (r.I > opts.strong_cutoff * r.sigma) {
naive(r);
continue;
}
if (I > opts.strong_cutoff * sigma) { naive_one(I, sigma, F, sigF); return; }
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<float>(post.mean_F);
r.sigmaF = static_cast<float>(std::sqrt(std::max(0.0, post.mean_I - post.mean_F * post.mean_F)));
const Posterior post = integrate_posterior(I, sigma, sigma_wilson, centric, opts.integration_points);
F = static_cast<float>(post.mean_F);
sigF = static_cast<float>(std::sqrt(std::max(0.0, post.mean_I - post.mean_F * post.mean_F)));
};
for (auto &r : merged) {
fw_one(r, r.I, r.sigma, r.F, r.sigmaF);
fw_one(r, r.I_plus, r.sigma_plus, r.F_plus, r.sigmaF_plus);
fw_one(r, r.I_minus, r.sigma_minus, r.F_minus, r.sigmaF_minus);
}
}
image_analysis/scale_merge/RotationScaleMerge.cpp
+57 -1
View File
@@ -11,6 +11,7 @@
#include <future>
#include <limits>
#include <random>
#include <unordered_map>
#include <gemmi/reciproc.hpp>
#include <gemmi/symmetry.hpp>
@@ -1340,7 +1341,8 @@ 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));
// French-Wilson (F, and F(+)/F(-) from the anomalous split) is deferred until after the anomalous
// accumulator below has attached I(+)/I(-), so the two hands get their amplitudes in one pass.
if (reject_count > 0)
logger.Info("Merge outlier rejection: dropped {} observations", reject_count);
@@ -1438,6 +1440,43 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
rmeas_num_all += factor * r.sum_abs_dev; rmeas_den_all += r.sum_I;
}
// ---- Anomalous split (always, even when the merge is Friedel-averaged): for each acentric
// reflection, the inverse-variance I(+)/I(-) from the SAME scaled fulls (so it never touches the
// Friedel-merged IMEAN, scaling or error model above). Lets I(+)/I(-) and F(+)/F(-) be written by
// default without scaling anomalously; a reflection with only one mate, or a centric, is left
// without the split. Skipped for the P1 search pass, which has no use for it. ----
struct AnomExport { float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN; };
std::unordered_map<uint64_t, AnomExport> anom_export;
if (!for_search) {
const int sg_num = x.GetSpaceGroupNumber().value_or(1);
const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, sg_num);
const gemmi::GroupOps gops = gemmi::find_spacegroup_by_number(sg_num)->operations();
struct AnomAcc { double swI[2] = {}; double sw[2] = {}; int32_t h = 0, k = 0, l = 0; }; // [hand] 0=I(+) 1=I(-)
std::unordered_map<uint64_t, AnomAcc> anom;
anom.reserve(result.merged.size() * 2 + 1);
for (const auto &o : fulls) {
if (!usable_merge(o)) continue;
if (rejected_obs[&o - fulls.data()]) continue; // outlier-rejected in the merge (CPU path)
const HKLKey ak = anom_keygen(o.h, o.k, o.l);
const int hand = ak.plus ? 0 : 1;
const float I_corr = o.I * o.corr;
const float sigma_corr = corrected_sigma(I_corr, o.sigma * o.corr, o.group);
if (!(sigma_corr > 0.0f) || !std::isfinite(sigma_corr)) continue;
const double w = 1.0 / (static_cast<double>(sigma_corr) * sigma_corr);
AnomAcc &a = anom[HKLKey{ak.h, ak.k, ak.l, true}.pack()];
a.h = ak.h; a.k = ak.k; a.l = ak.l;
a.swI[hand] += w * static_cast<double>(I_corr); a.sw[hand] += w;
}
for (const auto &[fkey, a] : anom) {
// Centrics have I(+)=I(-) by symmetry; leave them without an anomalous split.
if (gops.is_reflection_centric(gemmi::Op::Miller{a.h, a.k, a.l})) continue;
AnomExport ex;
if (a.sw[0] > 0.0) { ex.Ip = static_cast<float>(a.swI[0] / a.sw[0]); ex.sIp = static_cast<float>(1.0 / std::sqrt(a.sw[0])); }
if (a.sw[1] > 0.0) { ex.Im = static_cast<float>(a.swI[1] / a.sw[1]); ex.sIm = static_cast<float>(1.0 / std::sqrt(a.sw[1])); }
if (std::isfinite(ex.Ip) || std::isfinite(ex.Im)) anom_export[fkey] = ex;
}
}
MergeStatistics &out = result.statistics;
out.shells.resize(n_shells);
for (int s = 0; s < n_shells; ++s) {
@@ -1466,6 +1505,23 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
overall.cc_half = cc_half_overall.GetCC();
overall.cc_ref = NAN;
overall.r_meas = rmeas_den_all > 0.0 ? rmeas_num_all / rmeas_den_all : NAN;
// Attach the per-reflection anomalous split so the writer can emit I(+)/I(-) by default (each merged
// reflection maps to its Friedel-ASU key; in an anomalous merge both mates map to the same key).
if (!anom_export.empty()) {
const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, x.GetSpaceGroupNumber().value_or(1));
for (auto &r : result.merged) {
const HKLKey ak = anom_keygen(r.h, r.k, r.l);
const auto it = anom_export.find(HKLKey{ak.h, ak.k, ak.l, true}.pack());
if (it == anom_export.end()) continue;
r.I_plus = it->second.Ip; r.sigma_plus = it->second.sIp;
r.I_minus = it->second.Im; r.sigma_minus = it->second.sIm;
}
}
// French-Wilson amplitudes for IMEAN and (now that they are attached) each Bijvoet hand.
ApplyFrenchWilson(result.merged, x.GetSpaceGroupNumber().value_or(1));
logger.Info("Merge complete ({} unique reflections)", result.merged.size());
return result;
}