From 626afc49ddd719b57d137ec77cc56840b8179165 Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Sun, 25 Oct 2026 20:37:20 +0100
Subject: [PATCH] Refinement likelihood: the searches over D alone with the
 variance taken once

EstimateLikelihoodParameters' golden-section searches over D at a fixed beta - four per bin on the free
reflections, and the one on the working reflections the D/FOM change added, 50 evaluations over 95 % of
the data - evaluated -ln L reading every row through the index and taking log(var) at every evaluation.
The rows are now gathered per search with their variance and its logarithm; the sum is the same
expressions in the same order, so D, beta and everything after are the same to the bit (refined model and
report identical on the seven timing sets). refit.ml 200 -> 140 ms per cycle on the 1.1 A set, one thread per
trial slot.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit 23ad3bfc254d10af310f17dbb8f85616f666a783)
---
 .../structure_refinement/RefinementXray.cpp   | 47 +++++++++++++++++--
 1 file changed, 43 insertions(+), 4 deletions(-)

diff --git a/image_analysis/structure_refinement/RefinementXray.cpp b/image_analysis/structure_refinement/RefinementXray.cpp
index 8337c8c68..1ce87e07d 100644
--- a/image_analysis/structure_refinement/RefinementXray.cpp
+++ b/image_analysis/structure_refinement/RefinementXray.cpp
@@ -59,6 +59,43 @@ double BinNegLogLikelihood(const RefinementReflections &refl, const std::vector<
     return sum;
 }
 
+// Rows `idx` at a fixed beta, gathered for searches over D alone: the variance and its logarithm taken once rather
+// than at every evaluation. FixedVarianceNegLogLikelihood() at D is BinNegLogLikelihood() at D and that beta, to
+// the bit: the same expressions in the same order, only log(var) hoisted.
+struct FixedVarianceRows {
+    std::vector<double> fo, fc, var, log_var;
+    std::vector<char> centric;
+    FixedVarianceRows(const RefinementReflections &refl, const std::vector<double> &fc_all,
+                      const std::vector<size_t> &idx, double beta, bool with_sigma) {
+        for (size_t i : idx) {
+            fo.push_back(refl.fo[i]);
+            fc.push_back(fc_all[i]);
+            var.push_back(Variance(refl, i, beta, with_sigma));
+            log_var.push_back(std::log(var.back()));
+            centric.push_back(refl.centric[i]);
+        }
+    }
+};
+
+// refine_math::AmplitudeNegLogLikelihood() summed over the rows at D, with log(var) given.
+double FixedVarianceNegLogLikelihood(const FixedVarianceRows &rows, double d) {
+    double sum = 0;
+    for (size_t k = 0; k < rows.fo.size(); k++) {
+        const double fo = rows.fo[k], fc = rows.fc[k], var = rows.var[k];
+        if (!rows.centric[k]) {
+            const double x = 2.0 * d * fo * fc / var;
+            double log_i0, ratio;
+            refine_math::BesselI0I1(x, log_i0, ratio);
+            sum += rows.log_var[k] + (fo * fo + d * d * fc * fc) / var - log_i0;
+        } else {
+            const double x = d * fo * fc / var;
+            const double log_cosh = x + std::log(1.0 + std::exp(-2.0 * x)) - 0.6931471805599453;
+            sum += 0.5 * rows.log_var[k] + (fo * fo + d * d * fc * fc) / (2.0 * var) - log_cosh;
+        }
+    }
+    return sum;
+}
+
 // The minimum of f over [lo, hi] by golden-section search.
 template <typename F>
 double GoldenMinimum(F f, double lo, double hi, int iterations) {
@@ -304,7 +341,8 @@ void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::
         // beta: the spread of the free amplitudes about D Fc, measurement error included, by maximum likelihood
         // (alternating one-dimensional minimisations over ln D and ln beta)
         for (int round = 0; round < 4; round++) {
-            const double ld = GoldenMinimum([&](double x) { return BinNegLogLikelihood(refl, fc, idx, std::exp(x), beta, false); },
+            const FixedVarianceRows at_beta(refl, fc, idx, beta, false);
+            const double ld = GoldenMinimum([&](double x) { return FixedVarianceNegLogLikelihood(at_beta, std::exp(x)); },
                                             std::log(d) - 2.0, std::log(std::min(d * 7.4, 3.0)), 40);
             d = std::exp(ld);
             const double lb = GoldenMinimum([&](double x) { return BinNegLogLikelihood(refl, fc, idx, d, std::exp(x), false); },
@@ -315,9 +353,10 @@ void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::
         // Where the error is a large part of the spread, the variance counts it twice and D falls towards 0 -
         // the high-resolution D and figure of merit REFMAC 5.8 reports.
         // Following Murshudov et al. (2011) Acta Cryst. D67, 355-367 (D from the working set, the variance from the free)
-        if (!work_in_bin[b].empty())
-            d = GoldenMinimum([&](double x) { return BinNegLogLikelihood(refl, fc, work_in_bin[b], x, beta, true); },
-                              0.0, 3.0, 50);
+        if (!work_in_bin[b].empty()) {
+            const FixedVarianceRows working(refl, fc, work_in_bin[b], beta, true);
+            d = GoldenMinimum([&](double x) { return FixedVarianceNegLogLikelihood(working, x); }, 0.0, 3.0, 50);
+        }
         bin_d[b] = d;
         bin_beta[b] = beta;
     });