From 1d9709725d81e782ab5b0e1e6b5a117e57cbbe58 Mon Sep 17 00:00:00 2001 From: Bojan Nikolic Date: Wed, 30 Jul 2025 15:05:02 +0000 Subject: [PATCH 1/2] Save on standard deviation calculation Reuse means, and reuse squared deviations from the raw mean. The numerical stability of this approach should be very good as long as clipped mean is not very far from the raw mean. --- bdsf/functions.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/bdsf/functions.py b/bdsf/functions.py index 719ffa1..75fd38c 100644 --- a/bdsf/functions.py +++ b/bdsf/functions.py @@ -2143,7 +2143,8 @@ def bstat(indata, mask, kappa_npixbeam): maxiter = 200 converge_num = 1e-6 m_raw = numpy.mean(skpix) - r_raw = numpy.std(skpix, ddof=1) + d2 = (skpix-m_raw)**2 + r_raw = (d2.sum()/(ct-1))**0.5 while (c1 >= c2) and (iter < maxiter): npix = skpix.size @@ -2156,13 +2157,18 @@ def bstat(indata, mask, kappa_npixbeam): kappa = 3.0 lastct = ct medval = skpix[npix//2] - sig = numpy.std(skpix) + if iter == 0: + sig = r_raw + else: + m_new = numpy.mean(skpix) + sig = (d2.sum()/ct - (m_raw-m_new)**2)**0.5 wsm1, wsm2 = numpy.searchsorted(skpix, [medval - kappa*sig, medval + kappa*sig]) ct = wsm2 - wsm1 if ct > 0: skpix = skpix[wsm1:wsm2] + d2 = d2[wsm1:wsm2] c1 = abs(ct - lastct) c2 = converge_num * lastct From 34070da568b52045cefdcf1ff8900329789b16f9 Mon Sep 17 00:00:00 2001 From: Bojan Nikolic Date: Wed, 30 Jul 2025 15:39:14 +0000 Subject: [PATCH 2/2] Adapt also the last sigma calculation --- bdsf/functions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bdsf/functions.py b/bdsf/functions.py index 75fd38c..4e82d52 100644 --- a/bdsf/functions.py +++ b/bdsf/functions.py @@ -2176,7 +2176,7 @@ def bstat(indata, mask, kappa_npixbeam): mean = numpy.mean(skpix) median = skpix[skpix.size//2] - sigma = numpy.std(skpix, ddof=1) + sigma = (d2.sum()/(ct-1) - (m_raw-mean)**2 * ct/(ct-1) )**0.5 mode = 2.5*median - 1.5*mean if sigma > 0.0: skew_par = abs(mean - median)/sigma