aboutsummaryrefslogtreecommitdiff
path: root/julia/process_hkl.jl
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2024-02-22 17:49:24 +0100
committerThomas White <taw@physics.org>2024-02-22 17:49:24 +0100
commit5e8c1b3f964432a16f180cb64f199f3d02218f3f (patch)
tree57dab1cf5f016318ba8b2886e947e524d8629a58 /julia/process_hkl.jl
parentfa7c1e4aa507489f1c5ae9d9c107beaca5a705d3 (diff)
julia/process_hkl.jl: Abstract stable running mean/variance
Diffstat (limited to 'julia/process_hkl.jl')
-rw-r--r--julia/process_hkl.jl35
1 files 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