From 7f1f28f8d3fbf4bf763547614b202e48296746fe Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 8 Feb 2026 17:22:02 +0100 Subject: [PATCH] ScaleAndMerge: Merge intensities taken from scaling - variances taken directly from least-squares (with some magic from Opus 4.6, which I don't yet understand) --- image_analysis/scale_merge/ScaleAndMerge.cpp | 188 +++++++++---------- image_analysis/scale_merge/ScaleAndMerge.h | 12 +- 2 files changed, 96 insertions(+), 104 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index bc85d816..9b4b92e2 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -4,6 +4,7 @@ #include "ScaleAndMerge.h" #include +#include #include #include @@ -77,9 +78,6 @@ inline double SafeInv(double x, double fallback) { return 1.0 / x; } -// Canonicalize HKL according to Gemmi Reciprocal ASU if space group is provided. -// If merge_friedel==true -> Friedel mates collapse (key.is_positive always true). -// If merge_friedel==false -> keep I+ vs I- separate by key.is_positive. inline HKLKey CanonicalizeHKLKey(const Reflection& r, const ScaleMergeOptions& opt) { HKLKey key{}; key.h = r.h; @@ -87,7 +85,6 @@ inline HKLKey CanonicalizeHKLKey(const Reflection& r, const ScaleMergeOptions& o key.l = r.l; key.is_positive = true; - // If no SG provided, we can still optionally separate Friedel mates deterministically. if (!opt.space_group.has_value()) { if (!opt.merge_friedel) { const HKLKey neg{-r.h, -r.k, -r.l, true}; @@ -122,7 +119,7 @@ struct IntensityResidual { : Iobs_(static_cast(r.I)), inv_sigma_(SafeInv(sigma_obs, 1.0)), delta_phi_rad_(static_cast(r.delta_phi) * M_PI / 180.0), - lp_(SafeInv(static_cast(r.rlp), 1.0)), // rlp stores reciprocal Lorentz in this codebase + lp_(SafeInv(static_cast(r.rlp), 1.0)), half_wedge_rad_(wedge_deg * M_PI / 180.0 / 2.0), c1_(std::sqrt(2.0) * SafeInv(static_cast(r.zeta), 1.0)), s2_(s2), @@ -140,7 +137,6 @@ struct IntensityResidual { T partiality = T(1.0); if (partiality_model_) { - // partiality = 0.5 * [ erf( (Δφ+Δω/2) * (sqrt(2)*η/ζ) ) - erf( (Δφ-Δω/2) * (sqrt(2)*η/ζ) ) ] const T arg_plus = T(delta_phi_rad_ + half_wedge_rad_) * (T(c1_) * mosaicity_rad[0]); const T arg_minus = T(delta_phi_rad_ - half_wedge_rad_) * (T(c1_) * mosaicity_rad[0]); partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0); @@ -192,7 +188,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob int img_slot = -1; int hkl_slot = -1; double s2 = 0.0; - double sigma = 0.0; // sanitized sigma for weighting + double sigma = 0.0; }; std::vector obs; @@ -211,7 +207,6 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob if (!std::isfinite(r.I)) continue; - // Need valid zeta/rlp for the model to behave. if (!std::isfinite(r.zeta) || r.zeta <= 0.0f) continue; if (!std::isfinite(r.rlp) || r.rlp == 0.0f) @@ -274,17 +269,10 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob auto deg2rad = [](double deg) { return deg * M_PI / 180.0; }; - // Mosaicity: either global or per-image (always stored in radians internally) - std::vector mosaicity_rad; - double mosaicity_rad_global = deg2rad(opt.mosaicity_init_deg); + // Mosaicity: always per-image + std::vector mosaicity_rad(nimg, deg2rad(opt.mosaicity_init_deg)); - if (opt.mosaicity_per_image) { - mosaicity_rad.assign(nimg, deg2rad(opt.mosaicity_init_deg)); - } else { - mosaicity_rad_global = deg2rad(opt.mosaicity_init_deg); - } - - // Initialize Itrue from per-HKL median of observed intensities (rough; model contains partiality/LP) + // Initialize Itrue from per-HKL median of observed intensities { std::vector> per_hkl_I(nhkl); for (const auto& o : obs) { @@ -306,28 +294,21 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob ceres::Problem problem; - std::unique_ptr loss; - if (opt.use_huber_loss) { - loss = std::make_unique(opt.huber_delta); - } - const bool partiality_model = opt.wedge_deg > 0.0; for (const auto& o : obs) { - double* mosaicity_block = opt.mosaicity_per_image ? &mosaicity_rad[o.img_slot] : &mosaicity_rad_global; - auto* cost = new ceres::AutoDiffCostFunction( new IntensityResidual(*o.r, o.sigma, opt.wedge_deg, o.s2, opt.refine_b_factor, partiality_model)); problem.AddResidualBlock(cost, - loss.get(), + nullptr, // no loss function &log_k[o.img_slot], &b[o.img_slot], - mosaicity_block, + &mosaicity_rad[o.img_slot], &log_Itrue[o.hkl_slot]); } - // Optional Kabsch-like regularization for k: (k - 1)/sigma + // Optional Kabsch-like regularization for k if (opt.regularize_scale_to_one) { for (int i = 0; i < nimg; ++i) { auto* rcost = new ceres::AutoDiffCostFunction( @@ -336,7 +317,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob } } - // Fix gauge freedom: anchor first image scale to 1.0 + // Fix gauge freedom if (opt.fix_first_image_scale && nimg > 0) { log_k[0] = 0.0; problem.SetParameterBlockConstant(&log_k[0]); @@ -356,26 +337,19 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob } } - // Mosaicity refinement + bounds (degrees -> radians) + // Mosaicity refinement + bounds if (!opt.refine_mosaicity) { - if (opt.mosaicity_per_image) { - for (int i = 0; i < nimg; ++i) - problem.SetParameterBlockConstant(&mosaicity_rad[i]); - } else { - problem.SetParameterBlockConstant(&mosaicity_rad_global); - } + for (int i = 0; i < nimg; ++i) + problem.SetParameterBlockConstant(&mosaicity_rad[i]); } else { - const std::optional min_rad = opt.mosaicity_min_deg ? std::optional(deg2rad(*opt.mosaicity_min_deg)) : std::nullopt; - const std::optional max_rad = opt.mosaicity_max_deg ? std::optional(deg2rad(*opt.mosaicity_max_deg)) : std::nullopt; + const std::optional min_rad = opt.mosaicity_min_deg + ? std::optional(deg2rad(*opt.mosaicity_min_deg)) : std::nullopt; + const std::optional max_rad = opt.mosaicity_max_deg + ? std::optional(deg2rad(*opt.mosaicity_max_deg)) : std::nullopt; - if (opt.mosaicity_per_image) { - for (int i = 0; i < nimg; ++i) { - if (min_rad) problem.SetParameterLowerBound(&mosaicity_rad[i], 0, *min_rad); - if (max_rad) problem.SetParameterUpperBound(&mosaicity_rad[i], 0, *max_rad); - } - } else { - if (min_rad) problem.SetParameterLowerBound(&mosaicity_rad_global, 0, *min_rad); - if (max_rad) problem.SetParameterUpperBound(&mosaicity_rad_global, 0, *max_rad); + for (int i = 0; i < nimg; ++i) { + if (min_rad) problem.SetParameterLowerBound(&mosaicity_rad[i], 0, *min_rad); + if (max_rad) problem.SetParameterUpperBound(&mosaicity_rad[i], 0, *max_rad); } } @@ -389,74 +363,94 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob ceres::Solver::Summary summary; ceres::Solve(options, &problem, &summary); + // --- Export per-image results --- for (int i = 0; i < nimg; ++i) { out.image_scale_k[i] = std::exp(log_k[i]); out.image_b_factor[i] = opt.refine_b_factor ? b[i] : 0.0; } - // Export mosaicity (degrees) to result - if (opt.mosaicity_per_image) { - out.mosaicity_deg.resize(nimg); - for (int i = 0; i < nimg; ++i) - out.mosaicity_deg[i] = mosaicity_rad[i] * 180.0 / M_PI; - } else { - out.mosaicity_deg = { mosaicity_rad_global * 180.0 / M_PI }; + out.mosaicity_deg.resize(nimg); + for (int i = 0; i < nimg; ++i) + out.mosaicity_deg[i] = mosaicity_rad[i] * 180.0 / M_PI; + + // --- Compute goodness-of-fit (reduced chi-squared) --- + const int n_obs = static_cast(obs.size()); + // Count free parameters: nhkl log_Itrue + per-image (log_k + b + mosaicity) minus fixed ones + int n_params = nhkl; + for (int i = 0; i < nimg; ++i) { + if (!(opt.fix_first_image_scale && i == 0)) + n_params += 1; // log_k + if (opt.refine_b_factor) + n_params += 1; // b + if (opt.refine_mosaicity) + n_params += 1; // mosaicity } - // Reverse maps for merged output + const double half_wedge_rad = opt.wedge_deg * M_PI / 180.0 / 2.0; + + double sum_r2 = 0.0; + for (const auto& o : obs) { + const int i = o.img_slot; + const int h = o.hkl_slot; + + const double ki = std::exp(log_k[i]); + const double atten = opt.refine_b_factor ? std::exp(-b[i] * o.s2) : 1.0; + const double Itrue = std::exp(log_Itrue[h]); + + const double zeta = static_cast(o.r->zeta); + const double mosa_i = mosaicity_rad[i]; + const double c1 = std::sqrt(2.0) * (mosa_i / zeta); + + const double delta_phi_rad = static_cast(o.r->delta_phi) * M_PI / 180.0; + const double partiality = partiality_model + ? (std::erf((delta_phi_rad + half_wedge_rad) * c1) - std::erf((delta_phi_rad - half_wedge_rad) * c1)) / 2.0 + : 1.0; + + const double lp = SafeInv(static_cast(o.r->rlp), 1.0); + + const double Ipred = ki * atten * partiality * lp * Itrue; + const double resid = (Ipred - static_cast(o.r->I)) / o.sigma; + sum_r2 += resid * resid; + } + + const double gof2 = (n_obs > n_params) ? sum_r2 / static_cast(n_obs - n_params) : 1.0; + out.gof2 = gof2; + + // --- Covariance: σ(I_true) from (J^T W J)^{-1} scaled by GoF² --- std::vector slotToHKL(nhkl); for (const auto& kv : hklToSlot) { slotToHKL[kv.second] = kv.first; } - // sigma estimate consistent with the optimized forward model - std::vector n_per_hkl(nhkl, 0); - std::vector ss_per_hkl(nhkl, 0.0); - - const double half_wedge_rad = opt.wedge_deg * M_PI / 180.0 / 2.0; - - for (const auto& o : obs) { - const int i = o.img_slot; - const int h = o.hkl_slot; - - const double k = std::exp(log_k[i]); - const double atten = opt.refine_b_factor ? std::exp(-b[i] * o.s2) : 1.0; - const double Itrue = std::exp(log_Itrue[h]); - - const double zeta = static_cast(o.r->zeta); - const double mosa_i = opt.mosaicity_per_image ? mosaicity_rad[i] : mosaicity_rad_global; - const double c1 = std::sqrt(2.0) * (mosa_i / zeta); - - const double delta_phi_rad = static_cast(o.r->delta_phi) * M_PI / 180.0; - const double partiality = partiality_model ? - (std::erf((delta_phi_rad + half_wedge_rad) * c1) - std::erf((delta_phi_rad - half_wedge_rad) * c1)) / 2.0 - : 1.0; - - const double lp = SafeInv(static_cast(o.r->rlp), 1.0); - - const double Ipred = k * atten * partiality * lp * Itrue; - const double r = (Ipred - static_cast(o.r->I)); - - n_per_hkl[h] += 1; - ss_per_hkl[h] += r * r; + out.merged.resize(nhkl); + for (int h = 0; h < nhkl; ++h) { + out.merged[h].h = slotToHKL[h].h; + out.merged[h].k = slotToHKL[h].k; + out.merged[h].l = slotToHKL[h].l; + out.merged[h].I = std::exp(log_Itrue[h]); + out.merged[h].sigma = std::numeric_limits::quiet_NaN(); } - out.merged.reserve(nhkl); + ceres::Covariance::Options cov_options; + cov_options.algorithm_type = ceres::SPARSE_QR; + + ceres::Covariance covariance(cov_options); + + std::vector> covariance_blocks; + covariance_blocks.reserve(nhkl); for (int h = 0; h < nhkl; ++h) { - MergedReflection m{}; - m.h = slotToHKL[h].h; - m.k = slotToHKL[h].k; - m.l = slotToHKL[h].l; - m.I = static_cast(std::exp(log_Itrue[h])); + covariance_blocks.emplace_back(&log_Itrue[h], &log_Itrue[h]); + } - if (n_per_hkl[h] >= 2) { - const double rms = std::sqrt(ss_per_hkl[h] / static_cast(n_per_hkl[h] - 1)); - m.sigma = static_cast(rms / std::sqrt(static_cast(n_per_hkl[h]))); - } else { - m.sigma = std::numeric_limits::quiet_NaN(); + if (covariance.Compute(covariance_blocks, &problem)) { + for (int h = 0; h < nhkl; ++h) { + double var_log_I = 0.0; + covariance.GetCovarianceBlock(&log_Itrue[h], &log_Itrue[h], &var_log_I); + + // σ(I) = I * sqrt( var(log I) * GoF² ) + const double Itrue = std::exp(log_Itrue[h]); + out.merged[h].sigma = Itrue * std::sqrt(var_log_I * gof2); } - - out.merged.push_back(m); } return out; diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 848c7509..576b3e4c 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -13,9 +13,6 @@ struct ScaleMergeOptions { bool refine_b_factor = true; - bool use_huber_loss = true; - double huber_delta = 2.0; - int max_num_iterations = 100; double max_solver_time_s = 1.0; @@ -41,13 +38,11 @@ struct ScaleMergeOptions { double wedge_deg = 1.0; // --- Mosaicity (user input in degrees; internally converted to radians) --- - // If true, refine mosaicity. If mosaicity_per_image==true, refine one value per image. bool refine_mosaicity = true; - bool mosaicity_per_image = true; double mosaicity_init_deg = 0.17; // ~0.003 rad std::optional mosaicity_min_deg = 1e-3; - std::optional mosaicity_max_deg = 20.0; + std::optional mosaicity_max_deg = 2.0; // --- Optional: regularize per-image scale k towards 1 (Kabsch-like) --- bool regularize_scale_to_one = false; @@ -68,8 +63,11 @@ struct ScaleMergeResult { std::vector image_b_factor; std::vector image_ids; - // If mosaicity_per_image==true, one value per image (degrees). Otherwise size is 1 (global). + // One mosaicity value per image (degrees). std::vector mosaicity_deg; + + // Goodness-of-fit squared (reduced chi-squared). + double gof2 = 1.0; }; ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& observations,