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

@@ -4,6 +4,7 @@
 #include "ScaleAndMerge.h"
 
 #include <ceres/ceres.h>
+#include <ceres/covariance.h>
 
 #include <algorithm>
 #include <cmath>
@@ -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<double>(r.I)),
           inv_sigma_(SafeInv(sigma_obs, 1.0)),
           delta_phi_rad_(static_cast<double>(r.delta_phi) * M_PI / 180.0),
-          lp_(SafeInv(static_cast<double>(r.rlp), 1.0)), // rlp stores reciprocal Lorentz in this codebase
+          lp_(SafeInv(static_cast<double>(r.rlp), 1.0)),
           half_wedge_rad_(wedge_deg * M_PI / 180.0 / 2.0),
           c1_(std::sqrt(2.0) * SafeInv(static_cast<double>(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<Reflection>& 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<ObsRef> obs;
@@ -211,7 +207,6 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& 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<Reflection>& ob
 
     auto deg2rad = [](double deg) { return deg * M_PI / 180.0; };
 
-    // Mosaicity: either global or per-image (always stored in radians internally)
-    std::vector<double> mosaicity_rad;
-    double mosaicity_rad_global = deg2rad(opt.mosaicity_init_deg);
+    // Mosaicity: always per-image
+    std::vector<double> 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<std::vector<double>> per_hkl_I(nhkl);
         for (const auto& o : obs) {
@@ -306,28 +294,21 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& ob
 
     ceres::Problem problem;
 
-     std::unique_ptr<ceres::LossFunction> loss;
-    if (opt.use_huber_loss) {
-        loss = std::make_unique<ceres::HuberLoss>(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<IntensityResidual, 1, 1, 1, 1, 1>(
             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<ScaleRegularizationResidual, 1, 1>(
@@ -336,7 +317,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& 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<Reflection>& 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<double> min_rad = opt.mosaicity_min_deg ? std::optional<double>(deg2rad(*opt.mosaicity_min_deg)) : std::nullopt;
-        const std::optional<double> max_rad = opt.mosaicity_max_deg ? std::optional<double>(deg2rad(*opt.mosaicity_max_deg)) : std::nullopt;
+        const std::optional<double> min_rad = opt.mosaicity_min_deg
+            ? std::optional<double>(deg2rad(*opt.mosaicity_min_deg)) : std::nullopt;
+        const std::optional<double> max_rad = opt.mosaicity_max_deg
+            ? std::optional<double>(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<Reflection>& 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<int>(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<double>(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<double>(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<double>(o.r->rlp), 1.0);
+
+        const double Ipred = ki * atten * partiality * lp * Itrue;
+        const double resid = (Ipred - static_cast<double>(o.r->I)) / o.sigma;
+        sum_r2 += resid * resid;
+    }
+
+    const double gof2 = (n_obs > n_params) ? sum_r2 / static_cast<double>(n_obs - n_params) : 1.0;
+    out.gof2 = gof2;
+
+    // --- Covariance: σ(I_true) from (J^T W J)^{-1} scaled by GoF² ---
     std::vector<HKLKey> slotToHKL(nhkl);
     for (const auto& kv : hklToSlot) {
         slotToHKL[kv.second] = kv.first;
     }
 
-    // sigma estimate consistent with the optimized forward model
-    std::vector<int> n_per_hkl(nhkl, 0);
-    std::vector<double> 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<double>(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<double>(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<double>(o.r->rlp), 1.0);
-
-        const double Ipred = k * atten * partiality * lp * Itrue;
-        const double r = (Ipred - static_cast<double>(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<double>::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<std::pair<const double*, const double*>> 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<float>(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<double>(n_per_hkl[h] - 1));
-            m.sigma = static_cast<float>(rms / std::sqrt(static_cast<double>(n_per_hkl[h])));
-        } else {
-            m.sigma = std::numeric_limits<float>::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;

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<double> mosaicity_min_deg = 1e-3;
-    std::optional<double> mosaicity_max_deg = 20.0;
+    std::optional<double> 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<double> image_b_factor;
     std::vector<int> 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<double> mosaicity_deg;
+
+    // Goodness-of-fit squared (reduced chi-squared).
+    double gof2 = 1.0;
 };
 
 ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& observations,