diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index dd00c673..64361b47 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -135,6 +135,9 @@ void WriteMmcifReflections(const std::vector &reflections, ? 100.0 * static_cast(s.unique_reflections) / s.possible_unique_reflections : 0.0; }; if (!statistics.shells.empty()) { const auto &ov = statistics.overall; + // Anomalous signal-to-noise (SigAno) is written only when an anomalous split was made, so a + // non-anomalous merge keeps its previous stats block / shell-loop columns unchanged. + const bool has_anom = std::isfinite(ov.abs_diff_over_sigma_anomalous); out << "_reflns.d_resolution_high " << Fmt(ov.d_min, 2) << "\n"; out << "_reflns.d_resolution_low " << Fmt(ov.d_max, 2) << "\n"; out << "_reflns.number_obs " << ov.unique_reflections << "\n"; @@ -144,6 +147,9 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_reflns.pdbx_netI_over_sigmaI " << Fmt(ov.mean_i_over_sigma, 2) << "\n"; out << "_reflns.pdbx_Rrim_I_all " << Fmt(ov.r_meas, 4) << "\n"; out << "_reflns.pdbx_CC_half " << Fmt(ov.cc_half, 4) << "\n"; + if (has_anom) + out << "_reflns.pdbx_absDiff_over_sigma_anomalous " << Fmt(ov.abs_diff_over_sigma_anomalous, 3) + << " # SigAno = <|dano|>/\n"; out << "_reflns.jfjoch_diffrn_ISa " << CifStr(isa) << " # asymptotic I/sigma (Diederichs)\n"; // Dataset-wide isotropic Wilson B-factor estimate (standard PDBx item), analogous to XDS's // "WILSON LINE ... B=". Emitted only when the log-linear fit succeeded. @@ -171,13 +177,18 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_reflns_shell.meanI_over_sigI_obs\n"; out << "_reflns_shell.pdbx_Rrim_I_all\n"; out << "_reflns_shell.pdbx_CC_half\n"; + if (has_anom) + out << "_reflns_shell.pdbx_absDiff_over_sigma_anomalous\n"; for (const auto &s : statistics.shells) { if (s.unique_reflections == 0) continue; out << Fmt(s.d_min, 2) << " " << Fmt(s.d_max, 2) << " " << s.total_observations << " " << s.unique_reflections << " " << Fmt(mult(s), 2) << " " << Fmt(compl_pct(s), 1) << " " - << Fmt(s.mean_i_over_sigma, 2) << " " << Fmt(s.r_meas, 4) << " " << Fmt(s.cc_half, 4) << "\n"; + << Fmt(s.mean_i_over_sigma, 2) << " " << Fmt(s.r_meas, 4) << " " << Fmt(s.cc_half, 4); + if (has_anom) + out << " " << Fmt(s.abs_diff_over_sigma_anomalous, 3); + out << "\n"; } out << "#\n"; } diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index dd5aae6d..518fee26 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -886,7 +886,7 @@ std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in) { double multiplicity = in.unique_reflections > 0 ? static_cast(in.total_observations) / in.unique_reflections : 0.0; - output << fmt::format("{:8d} {:8d} {:8d} {:7.1f}% {:7.1f} {:8.1f} {:7.1f}% {:7.1f}% {:7.1f}%", + output << fmt::format("{:8d} {:8d} {:8d} {:7.1f}% {:7.1f} {:8.1f} {:7.1f}% {:7.1f}% {:7.1f}% {:8.2f}", in.total_observations, in.unique_reflections, in.possible_unique_reflections, @@ -895,17 +895,18 @@ std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in) { in.mean_i_over_sigma, in.r_meas*100.0, in.cc_half*100.0, - in.cc_ref*100.0); + in.cc_ref*100.0, + in.abs_diff_over_sigma_anomalous); return output; } std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) { output << std::endl; - output << fmt::format(" {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s}", - "d_min", "N_obs", "N_uniq", "N_possib", "Compl", "Mult", "", "R_meas", "CC1/2", "CCref") + output << fmt::format(" {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", + "d_min", "N_obs", "N_uniq", "N_possib", "Compl", "Mult", "", "R_meas", "CC1/2", "CCref", "SigAno") << std::endl; - output << fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->7s} {:->8s} {:->8s} {:->8s} {:->8s}", - "", "", "", "", "", "", "", "", "", "") << std::endl; + output << fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->7s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", + "", "", "", "", "", "", "", "", "", "", "") << std::endl; for (const auto &sh: in.shells) { if (sh.unique_reflections == 0) continue; @@ -913,8 +914,8 @@ std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) { output << sh; output << std::endl; } - output << fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->7s} {:->8s} {:->8s} {:->8s} {:->8s}", - "", "", "", "", "", "", "", "", "", "") << std::endl; + output << fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->7s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", + "", "", "", "", "", "", "", "", "", "", "") << std::endl; output << fmt::format(" {:>8s} ", "Overall"); output << in.overall; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index e0968e37..2636a1f7 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -46,6 +46,11 @@ struct MergeStatisticsShell { // Redundancy-independent merging R-factor (Diederichs & Karplus 1997), computed over the // observations that enter the merge: R_meas = sum_hkl sqrt(n/(n-1)) sum_i|I_i-| / sum I_i. double r_meas = NAN; + + // Anomalous signal-to-noise (XDS "SigAno" / mmCIF pdbx_absDiff_over_sigma_anomalous): + // <|I(+)-I(-)|> / over acentric reflections measured in both hands. NaN when + // there is no anomalous split (e.g. Friedel-merged with no mates, or the stills path). + double abs_diff_over_sigma_anomalous = NAN; }; struct MergeStatistics { diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 7bd73187..c7ca1e42 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -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 anom_export; + // SigAno = <|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 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 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(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(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(a.swI[0] / a.sw[0]); ex.sIp = static_cast(1.0 / std::sqrt(a.sw[0])); } if (a.sw[1] > 0.0) { ex.Im = static_cast(a.swI[1] / a.sw[1]); ex.sIm = static_cast(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).