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

rugnux: report anomalous signal-to-noise SigAno (stats table + mmCIF)

Browse Source 
Compute SigAno = <|I(+)-I(-)|>/<sigma(I(+)-I(-))> per report shell and overall
from the full-multiplicity I(+)/I(-) Bijvoet split, store it on
MergeStatisticsShell, and surface it in the console merge-statistics table (new
SigAno column) and the mmCIF via the standard wwPDB PDBx items
_reflns.pdbx_absDiff_over_sigma_anomalous (overall) and
_reflns_shell.pdbx_absDiff_over_sigma_anomalous (per shell) - the dictionary's
home for this quantity, so no jfjoch_ local item is needed (unlike ISa).

SigAno approaches ~0.8 for pure noise and rises above 1 with real anomalous
signal. A half-set anomalous CC is deliberately not used: its two half estimates
are complementary partitions of one observation pool (dI0 + dI1 = 2*dI_full), so
differencing the two Bijvoet hands cancels the large common intensity that keeps
CC1/2 non-negative and, once the anomalous SNR per half drops below 1, drives the
correlation towards -1 rather than 0 - misrepresenting a weak-but-real signal.

Emitted only when an anomalous split was made (has_anom guard), so a non-anomalous
merge keeps its stats block and shell-loop columns byte-for-byte. Export-neutral:
the I(+)/I(-)/IMEAN accumulation is unchanged. On low-energy S-SAD test data SigAno
tracks XDS (13 keV null ~0.89, 6 keV ~1.44, 5 keV ~1.5-1.7) and never goes
spuriously negative.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-17 20:54:03 +02:00
co-authored by Claude Opus 4.8
parent 1b5eb80ba4
commit 7fe12773de
4 changed files with 45 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+19 -2
View File
@@ -1447,11 +1447,16 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
// 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;
// SigAno = <|I(+)-I(-)|>/<sigma(I(+)-I(-))> accumulated on the standard report shells (num=sum|dI|,
// den=sum sigma(dI)); mmCIF pdbx_absDiff_over_sigma_anomalous. NaN'd below when there is no split.
std::vector<double> sig_num(n_shells, 0.0), sig_den(n_shells, 0.0);
double sig_num_all = 0.0, sig_den_all = 0.0;
size_t sig_n = 0;
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(-)
struct AnomAcc { double swI[2] = {}; double sw[2] = {}; int32_t h = 0, k = 0, l = 0; float d = NAN; }; // [hand] 0=I(+) 1=I(-)
std::unordered_map<uint64_t, AnomAcc> anom;
anom.reserve(result.merged.size() * 2 + 1);
for (const auto &o : fulls) {
@@ -1464,7 +1469,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
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.h = ak.h; a.k = ak.k; a.l = ak.l; a.d = o.d;
a.swI[hand] += w * static_cast<double>(I_corr); a.sw[hand] += w;
}
for (const auto &[fkey, a] : anom) {
@@ -1474,6 +1479,13 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
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;
if (a.sw[0] > 0.0 && a.sw[1] > 0.0 && a.d > 0.0f) {
const double abs_di = std::fabs(a.swI[0] / a.sw[0] - a.swI[1] / a.sw[1]);
const double sig_di = std::sqrt(1.0 / a.sw[0] + 1.0 / a.sw[1]);
const auto shell = shells.GetShell(a.d);
if (shell && *shell >= 0 && *shell < n_shells) { sig_num[*shell] += abs_di; sig_den[*shell] += sig_di; }
sig_num_all += abs_di; sig_den_all += sig_di; ++sig_n;
}
}
}
@@ -1491,6 +1503,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
ss.cc_half = sa[s].cc_half.GetCC();
ss.cc_ref = NAN;
ss.r_meas = rmeas_den[s] > 0.0 ? rmeas_num[s] / rmeas_den[s] : NAN;
ss.abs_diff_over_sigma_anomalous = sig_den[s] > 0.0 ? sig_num[s] / sig_den[s] : NAN;
}
auto &overall = out.overall;
overall.d_min = sd_min; overall.d_max = sd_max;
@@ -1505,6 +1518,10 @@ 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;
overall.abs_diff_over_sigma_anomalous = sig_den_all > 0.0 ? sig_num_all / sig_den_all : NAN;
if (std::isfinite(overall.abs_diff_over_sigma_anomalous))
logger.Info("Anomalous signal SigAno = {:.2f} ({} acentric pairs)",
overall.abs_diff_over_sigma_anomalous, sig_n);
// 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).