From 5e8c1b3f964432a16f180cb64f199f3d02218f3f Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Thu, 22 Feb 2024 17:49:24 +0100
Subject: julia/process_hkl.jl: Abstract stable running mean/variance

---
 julia/process_hkl.jl | 35 +++++++++++++++++++++++------------
 1 file changed, 23 insertions(+), 12 deletions(-)

diff --git a/julia/process_hkl.jl b/julia/process_hkl.jl
index 6d0dbd40..1d0bbdea 100644
--- a/julia/process_hkl.jl
+++ b/julia/process_hkl.jl
@@ -1,6 +1,21 @@
 using CrystFEL
 
 
+macro addmeasurement(measurement, weight,
+                     mean, sumweight, wksp)
+    return quote
+        delta = $(esc(measurement)) - $(esc(mean))
+        newsumweight = $(esc(sumweight)) + $(esc(weight))
+        R = delta * $(esc(weight)) / newsumweight
+        $(esc(mean)) += R
+        $(esc(wksp)) += $(esc(sumweight)) * delta * R
+        $(esc(sumweight)) = newsumweight
+    end
+end
+
+cstddev(nmeas, work1, work2) = sqrt(work2/work1)/sqrt(nmeas)
+
+
 function mergereflections(correction, crystalrefls, sym)
 
     merged = RefList{MergedReflection}(sym)
@@ -11,17 +26,10 @@ function mergereflections(correction, crystalrefls, sym)
 
             indices = asymmetricindices(sym, refl.indices)
             model_version = get!(merged, indices)
-
-            w = 1.0
-            mean = model_version.intensity
-            sumweight = model_version.temp1
-            M2 = model_version.temp2
-            temp = w + sumweight
-            delta = correction(refl.intensity, cr) - mean
-            R = delta * w / temp
-            model_version.intensity = mean + R
-            model_version.temp1 = temp
-            model_version.temp2 = M2 + sumweight * delta * R
+            @addmeasurement(correction(refl.intensity, cr), 1.0,
+                            model_version.intensity,
+                            model_version.temp1,
+                            model_version.temp2)
             model_version.nmeasurements += 1
 
         end
@@ -29,8 +37,11 @@ function mergereflections(correction, crystalrefls, sym)
 
     for refl in merged
         if refl.nmeasurements > 1
-            refl.sigintensity = sqrt(refl.temp2/refl.temp1)/sqrt(refl.nmeasurements)
+            refl.sigintensity = cstddev(refl.nmeasurements, refl.temp1, refl.temp2)
         else
+            # Cannot delete a reflection from a list (especially not
+            # while iterating), but setting nmeasurements to zero
+            # prevents it from being written to the output file.
             refl.nmeasurements = 0
         end
     end
-- 
cgit v1.2.3