From 63542ec9e052b147da54a656c9eed6ec35efe431 Mon Sep 17 00:00:00 2001
From: Wolfgang Brehm <wolfgang.brehm@gmail.com>
Date: Thu, 1 Aug 2019 13:10:51 +0200
Subject: add stable inline mean variance function

---
 libcrystfel/src/utils.h | 21 ++++++++++++++++++++-
 1 file changed, 20 insertions(+), 1 deletion(-)

(limited to 'libcrystfel/src/utils.h')

diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
index cd21bd93..44e6af6c 100644
--- a/libcrystfel/src/utils.h
+++ b/libcrystfel/src/utils.h
@@ -35,6 +35,7 @@
 
 #include <math.h>
 #include <complex.h>
+#include <float.h>
 #include <string.h>
 #include <stdlib.h>
 #include <pthread.h>
@@ -268,7 +269,25 @@ static inline struct quaternion invalid_quaternion(void)
 	return quat;
 }
 
-
+/* function to compute mean and variance stably
+ * \param x value
+ * \param w weight
+ * \param sumw pointer to accumulator variable for the sum of weights
+ * \param mean pointer to online mean value
+ * \param M2 pointer to online variance times sumw
+ */
+static inline void mean_variance(const double x,
+                                 const double w,
+                                 double* sumw, double* mean, double* M2)
+{
+	if (w<DBL_MIN) return;
+	const double temp  = w + *sumw;
+	const double delta = x - *mean;
+	const double R     = delta * w / temp;
+	*mean += R;
+	*M2   += *sumw*delta*R;
+	*sumw  = temp;
+}
 
 #ifdef __cplusplus
 }
-- 
cgit v1.2.3