From dd49d4e2f2e717b3a8f706a098066f1d36c58244 Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Sun, 25 Oct 2026 13:27:47 +0100
Subject: [PATCH] Refinement: the likelihood's D and beta fitted bin by bin on
 the refinement's threads

EstimateLikelihoodParameters fits each resolution bin on its own and writes only that bin's values, so
the bins run under ParallelFor; the result is the same, bit for bit. It was the largest single stage of
a refinement (107 calls, one core): 5reo 2.42 -> 0.56 s, 8dyz 2.46 -> 0.57 s, 9yl4 3.54 -> 0.88 s;
whole refinement 5reo 6.3 -> 4.8 s, 8dyz 6.9 -> 5.3 s (RTX 5080, -N 8, loaded box).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit 85aa34d3031c8e36be4462fca544bc801af45d5a)
(cherry picked from commit d84ba55252adf455d5b4ba62052dfbafa6290a64)
---
 image_analysis/structure_refinement/Refinement.cpp     | 2 +-
 image_analysis/structure_refinement/RefinementXray.cpp | 7 ++++---
 image_analysis/structure_refinement/RefinementXray.h   | 2 +-
 3 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp
index ea8771af8..560d3ff93 100644
--- a/image_analysis/structure_refinement/Refinement.cpp
+++ b/image_analysis/structure_refinement/Refinement.cpp
@@ -208,7 +208,7 @@ struct Context {
         }
         {
             const RefinementStage s4("refit.ml");
-            EstimateLikelihoodParameters(refl, fmodel, rows);
+            EstimateLikelihoodParameters(refl, fmodel, rows, nthreads);
         }
         engine->SetRowModel(rows);
         return ComputeRFactors(refl, fmodel);
diff --git a/image_analysis/structure_refinement/RefinementXray.cpp b/image_analysis/structure_refinement/RefinementXray.cpp
index ccbe8e3ab..61120304d 100644
--- a/image_analysis/structure_refinement/RefinementXray.cpp
+++ b/image_analysis/structure_refinement/RefinementXray.cpp
@@ -229,7 +229,7 @@ RefinementReflections PrepareRefinementReflections(const std::vector<MergedRefle
 }
 
 void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::vector<std::complex<float>> &fmodel,
-                                  XrayRowModel &rows) {
+                                  XrayRowModel &rows, size_t nthreads) {
     const size_t n = refl.Size();
     std::vector<double> fc(n);
     for (size_t i = 0; i < n; i++)
@@ -248,7 +248,8 @@ void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::
     const size_t bins = std::clamp<size_t>(free.size() / FREE_PER_BIN, 1, MAX_BINS);
 
     std::vector<double> centre(bins), bin_d(bins), bin_beta(bins);
-    for (size_t b = 0; b < bins; b++) {
+    // each bin is fitted on its own
+    ParallelFor(static_cast<int>(bins), nthreads, [&](int b) {
         const std::vector<size_t> idx(free.begin() + b * free.size() / bins, free.begin() + (b + 1) * free.size() / bins);
         double s2 = 0, sfo = 0, sfc2 = 0, sfofc = 0;
         for (size_t i : idx) {
@@ -271,7 +272,7 @@ void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::
         }
         bin_d[b] = d;
         bin_beta[b] = beta;
-    }
+    });
     rows.d.resize(n);
     rows.var.resize(n);
     for (size_t i = 0; i < n; i++) {
diff --git a/image_analysis/structure_refinement/RefinementXray.h b/image_analysis/structure_refinement/RefinementXray.h
index a0a459f21..fc37a9e85 100644
--- a/image_analysis/structure_refinement/RefinementXray.h
+++ b/image_analysis/structure_refinement/RefinementXray.h
@@ -58,7 +58,7 @@ struct XrayRowModel {
 // epsilon beta + sigma^2. `fmodel` is k_total (F + fmask), without D.
 // Following Lunin & Skovoroda (1995) Acta Cryst. A51, 880-887, and Afonine et al. (2012) Acta Cryst. D68, 352-367
 void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::vector<std::complex<float>> &fmodel,
-                                  XrayRowModel &rows);
+                                  XrayRowModel &rows, size_t nthreads = 1);
 
 // R-work and R-free of k_total (F + fmask) against the observations.
 struct RFactors {