diff --git a/common/Reflection.h b/common/Reflection.h index 37afcf9c..3b74a695 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -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; }; diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 6f54e0d6..dd00c673 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -6,6 +6,7 @@ #include "scale_merge/HKLKey.h" #include "scale_merge/TwinningAnalysis.h" +#include #include #include #include @@ -188,8 +189,16 @@ void WriteMmcifReflections(const std::vector &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 &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 &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(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(r.h)); mtz.data.push_back(static_cast(r.k)); mtz.data.push_back(static_cast(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); diff --git a/image_analysis/scale_merge/FrenchWilson.cpp b/image_analysis/scale_merge/FrenchWilson.cpp index 3e5a1cde..bbb7c169 100644 --- a/image_analysis/scale_merge/FrenchWilson.cpp +++ b/image_analysis/scale_merge/FrenchWilson.cpp @@ -59,16 +59,24 @@ Posterior integrate_posterior(double I, double sigma, double sigma_wilson, bool 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; + // 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(std::sqrt(ip)); + sigF = (ip > 0.0 && std::isfinite(sigma)) ? static_cast(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 &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 &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(post.mean_F); - r.sigmaF = static_cast(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(post.mean_F); + sigF = static_cast(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); } } diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 35180df5..7bd73187 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -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 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 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(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(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(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; + } + } + 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; }