From 86fe70c135e084a9d1b0709ac610e060f696fa51 Mon Sep 17 00:00:00 2001 From: Jerome Kieffer Date: Fri, 15 Nov 2024 09:55:08 +0100 Subject: [PATCH 01/13] Implement the sort function in cython --- src/pyFAI/ext/CSR_common.pxi | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/src/pyFAI/ext/CSR_common.pxi b/src/pyFAI/ext/CSR_common.pxi index 847f8af80..4ecb541de 100644 --- a/src/pyFAI/ext/CSR_common.pxi +++ b/src/pyFAI/ext/CSR_common.pxi @@ -29,11 +29,15 @@ __author__ = "Jérôme Kieffer" __contact__ = "Jerome.kieffer@esrf.fr" -__date__ = "06/09/2024" +__date__ = "15/11/2024" __status__ = "stable" __license__ = "MIT" +from libcpp cimport bool +from libcpp.algorithm cimport sort +from cython cimport floating + import cython from cython.parallel import prange import numpy @@ -42,6 +46,23 @@ from .preproc import preproc from ..containers import Integrate1dtpl, Integrate2dtpl, ErrorModel +cdef struct float4: + float s0 + float s1 + float s2 + float s3 + +cdef bool inline cmp(float4 a, float4 b) noexcept nogil: + return True if a.s0 Date: Fri, 15 Nov 2024 10:07:16 +0100 Subject: [PATCH 02/13] Update docstring --- src/pyFAI/ext/CSR_common.pxi | 274 ++++++++++++++++++++++++++++++++++- 1 file changed, 270 insertions(+), 4 deletions(-) diff --git a/src/pyFAI/ext/CSR_common.pxi b/src/pyFAI/ext/CSR_common.pxi index 4ecb541de..7891aaa95 100644 --- a/src/pyFAI/ext/CSR_common.pxi +++ b/src/pyFAI/ext/CSR_common.pxi @@ -319,8 +319,8 @@ cdef class CsrIntegrator(object): :type absorption: ndarray :param normalization_factor: divide the valid result by this value :param bool weighted_average: set to False to use an unweighted mean (similar to legacy) instead of the weighted average. WIP - :return: positions, pattern, weighted_histogram and unweighted_histogram - :rtype: Integrate1dtpl 4-named-tuple of ndarrays + :return: namedtuple with "position intensity sigma signal variance normalization count std sem norm_sq" + :rtype: Integrate1dtpl named-tuple of ndarrays """ cdef: index_t i, j, idx = 0 @@ -506,8 +506,274 @@ cdef class CsrIntegrator(object): :param error_model: set to "poissonian" to use signal as variance (minimum 1), "azimuthal" to use the variance in a ring. :param normalization_factor: divide the valid result by this value - :return: positions, pattern, weighted_histogram and unweighted_histogram - :rtype: Integrate1dtpl 4-named-tuple of ndarrays + :return: namedtuple with "position intensity sigma signal variance normalization count std sem norm_sq" + :rtype: Integrate1dtpl named-tuple of ndarrays + + """ + error_model = ErrorModel.parse(error_model) + cdef: + index_t i, j, c, bad_pix, idx = 0 + acc_t acc_sig = 0.0, acc_var = 0.0, acc_norm = 0.0, acc_count = 0.0, coef = 0.0, acc_norm_sq=0.0 + acc_t sig, norm, count, var + acc_t delta1, delta2, b, x, omega_A, omega_B, aver, std, chauvenet_cutoff, omega2_A, omega2_B, w + data_t empty + acc_t[::1] sum_sig = numpy.empty(self.output_size, dtype=acc_d) + acc_t[::1] sum_var = numpy.empty(self.output_size, dtype=acc_d) + acc_t[::1] sum_norm = numpy.empty(self.output_size, dtype=acc_d) + acc_t[::1] sum_norm_sq = numpy.empty(self.output_size, dtype=acc_d) + acc_t[::1] sum_count = numpy.empty(self.output_size, dtype=acc_d) + data_t[::1] merged = numpy.empty(self.output_size, dtype=data_d) + data_t[::1] stda = numpy.empty(self.output_size, dtype=data_d) + data_t[::1] sema = numpy.empty(self.output_size, dtype=data_d) + data_t[:, ::1] preproc4 + bint do_azimuthal_variance = error_model == ErrorModel.AZIMUTHAL + bint do_hybrid_variance = error_model == ErrorModel.HYBRID + assert weights.size == self.input_size, "weights size" + empty = dummy if dummy is not None else self.empty + #Call the preprocessor ... + preproc4 = preproc(weights.ravel(), + dark=dark, + flat=flat, + solidangle=solidangle, + polarization=polarization, + absorption=absorption, + mask=self.cmask if self.check_mask else None, + dummy=dummy, + delta_dummy=delta_dummy, + normalization_factor=normalization_factor, + empty=self.empty, + split_result=4, + variance=variance, + dtype=data_d, + error_model=error_model, + out=self.preprocessed) + with nogil: + # Integrate once + for i in prange(self.output_size, schedule="guided"): + acc_sig = 0.0 + acc_var = 0.0 + acc_norm = 0.0 + acc_norm_sq = 0.0 + acc_count = 0.0 + for j in range(self._indptr[i], self._indptr[i + 1]): + coef = self._data[j] + if coef == 0.0: + continue + idx = self._indices[j] + sig = preproc4[idx, 0] + var = preproc4[idx, 1] + norm = preproc4[idx, 2] + count = preproc4[idx, 3] + if isnan(sig) or isnan(var) or isnan(norm) or isnan(count) or (norm == 0.0) or (count == 0.0): + continue + w = coef * norm + if do_azimuthal_variance or (do_hybrid_variance and cycle): + if acc_norm == 0.0: + # Initialize the accumulator with data from the pixel + acc_sig = coef * sig + #Variance remains at 0 + acc_norm = w + acc_norm_sq = w*w + acc_count = coef * count + else: + omega_A = acc_norm + omega_B = coef * norm # ω_i = c_i * norm_i + omega2_A = acc_norm_sq + omega2_B = omega_B*omega_B + acc_norm = omega_A + omega_B + acc_norm_sq = omega2_A + omega2_B + # omega3 = acc_norm * omega_A * omega2_B + # VV_{AUb} = VV_A + ω_b^2 * (b-) * (b-) + b = sig / norm + delta1 = acc_sig/omega_A - b + acc_sig = acc_sig + coef * sig + delta2 = acc_sig / acc_norm - b + acc_var = acc_var + omega2_B * delta1 * delta2 + acc_count = acc_count + coef * count + else: + acc_sig = acc_sig + coef * sig + acc_var = acc_var + coef * coef * var + acc_norm = acc_norm + w + acc_norm_sq = acc_norm_sq + w*w + acc_count = acc_count + coef * count + + if (acc_norm_sq > 0.0): + aver = acc_sig / acc_norm + std = sqrt(acc_var / acc_norm_sq) + # sem = sqrt(acc_var) / acc_norm + else: + aver = NAN; + std = NAN; + # sem = NAN; + + # cycle for sigma-clipping + for c in range(cycle): + # Sigma-clip + if (acc_norm == 0.0) or (acc_count == 0.0): + break + + #This is for Chauvenet's criterion + acc_count = max(3.0, acc_count) + chauvenet_cutoff = max(cutoff, sqrt(2.0*log(acc_count/sqrt(2.0*M_PI)))) + # Discard outliers + bad_pix = 0 + for j in range(self._indptr[i], self._indptr[i + 1]): + coef = self._data[j] + if coef == 0.0: + continue + idx = self._indices[j] + sig = preproc4[idx, 0] + var = preproc4[idx, 1] + norm = preproc4[idx, 2] + count = preproc4[idx, 3] + if isnan(sig) or isnan(var) or isnan(norm) or isnan(count) or (norm == 0.0) or (count == 0.0): + continue + # Check validity (on cnt, i.e. s3) and normalisation (in s2) value to avoid zero division + x = sig / norm + if fabs(x - aver) > chauvenet_cutoff * std: + preproc4[idx, 3] = NAN; + bad_pix = bad_pix + 1 + if bad_pix == 0: + if do_hybrid_variance: + c = cycle #enforce the leave of the loop + else: + break + + #integrate again + acc_sig = 0.0 + acc_var = 0.0 + acc_norm = 0.0 + acc_norm_sq = 0.0 + acc_count = 0.0 + + for j in range(self._indptr[i], self._indptr[i + 1]): + coef = self._data[j] + if coef == 0.0: + continue + idx = self._indices[j] + sig = preproc4[idx, 0] + var = preproc4[idx, 1] + norm = preproc4[idx, 2] + count = preproc4[idx, 3] + if isnan(sig) or isnan(var) or isnan(norm) or isnan(count) or (norm == 0.0) or (count == 0.0): + continue + w = coef * norm + if do_azimuthal_variance or (do_hybrid_variance and c+1) * (b-) + b = sig / norm + delta1 = acc_sig/omega_A - b + acc_sig = acc_sig + coef * sig + delta2 = acc_sig / acc_norm - b + acc_var = acc_var + omega2_B * delta1 * delta2 + acc_count = acc_count + coef * count + else: + acc_sig = acc_sig + coef * sig + acc_var = acc_var + coef * coef * var + acc_norm = acc_norm + w + acc_norm_sq = acc_norm_sq + w*w + acc_count = acc_count + coef * count + if (acc_norm_sq > 0.0): + aver = acc_sig / acc_norm + std = sqrt(acc_var / acc_norm_sq) + # sem = sqrt(acc_var) / acc_norm # Not needed yet + else: + aver = NAN; + std = NAN; + # sem = NAN; + if bad_pix == 0: + break + + #collect things ... + sum_sig[i] = acc_sig + sum_var[i] = acc_var + sum_norm[i] = acc_norm + sum_norm_sq[i] = acc_norm_sq + sum_count[i] = acc_count + if (acc_norm_sq > 0.0): + merged[i] = aver # calculated previously + stda[i] = std # sqrt(acc_var / acc_norm_sq) # already calculated previously + sema[i] = sqrt(acc_var) / acc_norm # not calculated previously since it was not needed + else: + merged[i] = empty + stda[i] = empty + sema[i] = empty + + #"position intensity error signal variance normalization count" + return Integrate1dtpl(self.bin_centers, + numpy.asarray(merged),numpy.asarray(sema) , + numpy.asarray(sum_sig),numpy.asarray(sum_var), + numpy.asarray(sum_norm), numpy.asarray(sum_count), + numpy.asarray(stda), numpy.asarray(sema), numpy.asarray(sum_norm_sq)) + + + def medfilt( self, + weights, + dark=None, + dummy=None, + delta_dummy=None, + variance=None, + dark_variance=None, + flat=None, + solidangle=None, + polarization=None, + absorption=None, + bint safe=True, + error_model=ErrorModel.NO, + data_t normalization_factor=1.0, + double quant_min=0.5, + double quant_max=0.5, + ): + """Perform a median filter/quantile averaging in azimuthal space + + Else, the error is propagated like Poisson or pre-defined variance, no azimuthal variance for now. + + Integration is performed using the CSR representation of the look-up table on all + arrays: signal, variance, normalization and count + + All data are duplicated, sorted and the relevant values (i.e. within [quant_min..quant_max]) + are averaged like in integrate_ng + + :param weights: input image + :type weights: ndarray + :param dark: array with the dark-current value to be subtracted (if any) + :type dark: ndarray + :param dummy: value for dead pixels (optional) + :type dummy: float + :param delta_dummy: precision for dead-pixel value in dynamic masking + :type delta_dummy: float + :param variance: the variance associate to the image + :type variance: ndarray + :param dark_variance: the variance associate to the dark + :type dark_variance: ndarray + :param flat: array with the dark-current value to be divided by (if any) + :type flat: ndarray + :param solidAngle: array with the solid angle of each pixel to be divided by (if any) + :type solidAngle: ndarray + :param polarization: array with the polarization correction values to be divided by (if any) + :type polarization: ndarray + :param absorption: Apparent efficiency of a pixel due to parallax effect + :type absorption: ndarray + :param safe: set to True to save some tests + :param error_model: set to "poissonian" to use signal as variance (minimum 1), "azimuthal" to use the variance in a ring. + :param normalization_factor: divide the valid result by this value + + :return: namedtuple with "position intensity sigma signal variance normalization count std sem norm_sq" + :rtype: Integrate1dtpl named-tuple of ndarrays """ error_model = ErrorModel.parse(error_model) cdef: From 6b7f99f8694f52dfcd43685fdfb25794a916db95 Mon Sep 17 00:00:00 2001 From: Jerome Kieffer Date: Fri, 15 Nov 2024 14:42:39 +0100 Subject: [PATCH 03/13] inline functions --- src/pyFAI/ext/CSR_common.pxi | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pyFAI/ext/CSR_common.pxi b/src/pyFAI/ext/CSR_common.pxi index 7891aaa95..9e862fc83 100644 --- a/src/pyFAI/ext/CSR_common.pxi +++ b/src/pyFAI/ext/CSR_common.pxi @@ -52,10 +52,10 @@ cdef struct float4: float s2 float s3 -cdef bool inline cmp(float4 a, float4 b) noexcept nogil: +cdef inline bool cmp(float4 a, float4 b) noexcept nogil: return True if a.s0 Date: Fri, 15 Nov 2024 15:53:10 +0100 Subject: [PATCH 04/13] Implement cython medfilt --- src/pyFAI/ext/CSR_common.pxi | 246 ++++++++-------------------- src/pyFAI/ext/meson.build | 5 +- src/pyFAI/ext/splitPixelFullCSR.pyx | 4 +- 3 files changed, 77 insertions(+), 178 deletions(-) diff --git a/src/pyFAI/ext/CSR_common.pxi b/src/pyFAI/ext/CSR_common.pxi index 9e862fc83..19ac09980 100644 --- a/src/pyFAI/ext/CSR_common.pxi +++ b/src/pyFAI/ext/CSR_common.pxi @@ -38,6 +38,7 @@ from libcpp cimport bool from libcpp.algorithm cimport sort from cython cimport floating +import os import cython from cython.parallel import prange import numpy @@ -45,18 +46,25 @@ import numpy from .preproc import preproc from ..containers import Integrate1dtpl, Integrate2dtpl, ErrorModel +# cdef Py_ssize_t MAX_THREADS = 8 +# try: +# MAX_THREADS = min(MAX_THREADS, len(os.sched_getaffinity(os.getpid()))) # Limit to the actual number of threads +# except Exception: +# MAX_THREADS = min(MAX_THREADS, os.cpu_count() or 1) -cdef struct float4: + +cdef struct float4_t: float s0 float s1 float s2 float s3 +float4_d = numpy.dtype([('s0','f4'),('s1','f4'),('s2','f4'),('s3','f4')]) -cdef inline bool cmp(float4 a, float4 b) noexcept nogil: +cdef inline bool cmp(float4_t a, float4_t b) noexcept nogil: return True if a.s0) * (b-) - b = sig / norm - delta1 = acc_sig/omega_A - b - acc_sig = acc_sig + coef * sig - delta2 = acc_sig / acc_norm - b - acc_var = acc_var + omega2_B * delta1 * delta2 - acc_count = acc_count + coef * count - else: - acc_sig = acc_sig + coef * sig - acc_var = acc_var + coef * coef * var - acc_norm = acc_norm + w - acc_norm_sq = acc_norm_sq + w*w - acc_count = acc_count + coef * count - - if (acc_norm_sq > 0.0): - aver = acc_sig / acc_norm - std = sqrt(acc_var / acc_norm_sq) - # sem = sqrt(acc_var) / acc_norm - else: - aver = NAN; - std = NAN; - # sem = NAN; - - # cycle for sigma-clipping - for c in range(cycle): - # Sigma-clip - if (acc_norm == 0.0) or (acc_count == 0.0): - break - - #This is for Chauvenet's criterion - acc_count = max(3.0, acc_count) - chauvenet_cutoff = max(cutoff, sqrt(2.0*log(acc_count/sqrt(2.0*M_PI)))) - # Discard outliers - bad_pix = 0 - for j in range(self._indptr[i], self._indptr[i + 1]): - coef = self._data[j] - if coef == 0.0: - continue - idx = self._indices[j] - sig = preproc4[idx, 0] - var = preproc4[idx, 1] - norm = preproc4[idx, 2] - count = preproc4[idx, 3] - if isnan(sig) or isnan(var) or isnan(norm) or isnan(count) or (norm == 0.0) or (count == 0.0): - continue - # Check validity (on cnt, i.e. s3) and normalisation (in s2) value to avoid zero division - x = sig / norm - if fabs(x - aver) > chauvenet_cutoff * std: - preproc4[idx, 3] = NAN; - bad_pix = bad_pix + 1 - if bad_pix == 0: - if do_hybrid_variance: - c = cycle #enforce the leave of the loop - else: - break - - #integrate again - acc_sig = 0.0 - acc_var = 0.0 - acc_norm = 0.0 - acc_norm_sq = 0.0 - acc_count = 0.0 - - for j in range(self._indptr[i], self._indptr[i + 1]): - coef = self._data[j] - if coef == 0.0: - continue - idx = self._indices[j] - sig = preproc4[idx, 0] - var = preproc4[idx, 1] - norm = preproc4[idx, 2] - count = preproc4[idx, 3] - if isnan(sig) or isnan(var) or isnan(norm) or isnan(count) or (norm == 0.0) or (count == 0.0): - continue - w = coef * norm - if do_azimuthal_variance or (do_hybrid_variance and c+1) * (b-) - b = sig / norm - delta1 = acc_sig/omega_A - b - acc_sig = acc_sig + coef * sig - delta2 = acc_sig / acc_norm - b - acc_var = acc_var + omega2_B * delta1 * delta2 - acc_count = acc_count + coef * count - else: - acc_sig = acc_sig + coef * sig - acc_var = acc_var + coef * coef * var - acc_norm = acc_norm + w - acc_norm_sq = acc_norm_sq + w*w - acc_count = acc_count + coef * count - if (acc_norm_sq > 0.0): - aver = acc_sig / acc_norm - std = sqrt(acc_var / acc_norm_sq) - # sem = sqrt(acc_var) / acc_norm # Not needed yet - else: - aver = NAN; - std = NAN; - # sem = NAN; - if bad_pix == 0: - break - + # Duplicate the input data and populate the large work-array + for i in range(npix): # NOT faster in parallel ! + j = self._indices[i] + sig = preproc4[j,0] + var = preproc4[j,1] + nrm = preproc4[j,2] + weight = self._data[i] + element.s0 = sig/nrm # average signal + element.s1 = sig * weight # weighted raw signal + element.s2 = var * weight * weight # weighted raw variance + element.s3 = nrm * weight # weighted raw normalization + work[i] = element + for idx in prange(self.output_size, schedule="guided"): + start = self._indptr[idx] + stop = self._indptr[idx+1] + acc_sig = acc_var = acc_norm = acc_norm_sq = 0.0 + sort_float4(work[start:stop]) + + cumsum = 0.0 + for i in range(start, stop): + cumsum = cumsum + work[i].s3 + work[i].s0 = cumsum + qmin = quant_min * cumsum + qmax = quant_max * cumsum + + element.s0 = 0.0 + for i in range(start, stop): + former_element = element + element = work[i] + if (former_element.s0 >= qmin) and (element.s0 <= qmax): + acc_sig = acc_sig + element.s1 + acc_var = acc_var + element.s2 + acc_norm = acc_norm + element.s3 + acc_norm_sq = acc_norm_sq + element.s3*element.s3 + cnt = cnt + 1.0 #collect things ... - sum_sig[i] = acc_sig - sum_var[i] = acc_var - sum_norm[i] = acc_norm - sum_norm_sq[i] = acc_norm_sq - sum_count[i] = acc_count - if (acc_norm_sq > 0.0): - merged[i] = aver # calculated previously - stda[i] = std # sqrt(acc_var / acc_norm_sq) # already calculated previously - sema[i] = sqrt(acc_var) / acc_norm # not calculated previously since it was not needed + sum_sig[idx] = acc_sig + sum_var[idx] = acc_var + sum_norm[idx] = acc_norm + sum_norm_sq[idx] = acc_norm_sq + sum_count[idx] = cnt + if (acc_norm_sq): + merged[idx] = acc_sig/acc_var + stda[idx] = sqrt(acc_var / acc_norm_sq) + sema[idx] = sqrt(acc_var) / acc_norm else: - merged[i] = empty - stda[i] = empty - sema[i] = empty + merged[idx] = empty + stda[idx] = empty + sema[idx] = empty - #"position intensity error signal variance normalization count" + #"position intensity sigma signal variance normalization count std sem norm_sq" return Integrate1dtpl(self.bin_centers, numpy.asarray(merged),numpy.asarray(sema) , numpy.asarray(sum_sig),numpy.asarray(sum_var), diff --git a/src/pyFAI/ext/meson.build b/src/pyFAI/ext/meson.build index 0e26b8c8e..140e2c87e 100644 --- a/src/pyFAI/ext/meson.build +++ b/src/pyFAI/ext/meson.build @@ -83,8 +83,9 @@ py.extension_module('watershed', 'watershed.pyx', py.extension_module('_tree', '_tree.pyx', dependencies : py_dep, install: true, subdir: 'pyFAI/ext') -py.extension_module('sparse_utils', 'sparse_utils.pyx', #Deprecated - dependencies : py_dep, install: true, subdir: 'pyFAI/ext') +py.extension_module('sparse_utils', 'sparse_utils.pyx', + override_options : ['cython_language=cpp'], + dependencies : [py_dep, omp], install: true, subdir: 'pyFAI/ext') py.extension_module('inpainting', 'inpainting.pyx', dependencies : py_dep, install: true, subdir: 'pyFAI/ext') diff --git a/src/pyFAI/ext/splitPixelFullCSR.pyx b/src/pyFAI/ext/splitPixelFullCSR.pyx index e52d37017..c358f04f0 100644 --- a/src/pyFAI/ext/splitPixelFullCSR.pyx +++ b/src/pyFAI/ext/splitPixelFullCSR.pyx @@ -2,7 +2,7 @@ #cython: embedsignature=True, language_level=3, binding=True #cython: boundscheck=False, wraparound=False, cdivision=True, initializedcheck=False, ## This is for developping -## cython: profile=True, warn.undeclared=True, warn.unused=True, warn.unused_result=False, warn.unused_arg=True +##cython: profile=True, warn.undeclared=True, warn.unused=True, warn.unused_result=False, warn.unused_arg=True # # Project: Fast Azimuthal Integration # https://github.com/silx-kit/pyFAI @@ -35,7 +35,7 @@ Sparse matrix represented using the CompressedSparseRow. __author__ = "Jérôme Kieffer" __contact__ = "Jerome.kieffer@esrf.fr" -__date__ = "04/10/2023" +__date__ = "15/11/2024" __status__ = "stable" __license__ = "MIT" From 5f61cf8d379b584771714920ac9f2a75aa93add1 Mon Sep 17 00:00:00 2001 From: Jerome Kieffer Date: Mon, 18 Nov 2024 11:22:22 +0100 Subject: [PATCH 05/13] Use same API as cython --- src/pyFAI/engines/CSR_engine.py | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/src/pyFAI/engines/CSR_engine.py b/src/pyFAI/engines/CSR_engine.py index 178f8903c..4a79814dd 100644 --- a/src/pyFAI/engines/CSR_engine.py +++ b/src/pyFAI/engines/CSR_engine.py @@ -26,7 +26,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "14/11/2024" +__date__ = "18/11/2024" __status__ = "development" from collections.abc import Iterable @@ -395,7 +395,9 @@ def medfilt(self, data, dark=None, dummy=None, delta_dummy=None, variance=None, dark_variance=None, flat=None, solidangle=None, polarization=None, absorption=None, safe=True, error_model=None, - normalization_factor=1.0, quantile=0.5 + normalization_factor=1.0, + quant_min=0.5, + quant_max=0.5, ): """ Perform a median-filter/quantile mean in azimuthal space. @@ -425,17 +427,12 @@ def medfilt(self, data, dark=None, dummy=None, delta_dummy=None, :param safe: Unused in this implementation :param error_model: Enum or str, "azimuthal" or "poisson" :param normalization_factor: divide raw signal by this value - :param quantile: which percentile/100 use for cutting out quantil. - can be a 2-tuple to specify a region to average out. - By default, takes the median + :param quant_min: start percentile/100 to use. Use 0.5 for the median (default). 0<=quant_min<=1 + :param quant_max: stop percentile/100 to use. Use 0.5 for the median (default). 0<=quant_max<=1 + :return: namedtuple with "position intensity error signal variance normalization count" """ - if isinstance(quantile, Iterable): - q_start = min(quantile) - q_stop = max(quantile) - else: - q_stop = q_start = quantile indptr = self._csr.indptr indices = self._csr.indices @@ -484,7 +481,7 @@ def medfilt(self, data, dark=None, dummy=None, delta_dummy=None, upper = numpy.cumsum(tmp["norm"]) last = upper[-1] lower = numpy.concatenate(([0],upper[:-1])) - mask = numpy.logical_and(upper>=q_start*last, lower<=q_stop*last) + mask = numpy.logical_and(upper>=quant_min*last, lower<=quant_max*last) tmp = tmp[mask] cnt[i] = tmp.size signal[i] = tmp["sig"].sum(dtype=numpy.float64) From 36576089102411c9486178e10432c84c293d596b Mon Sep 17 00:00:00 2001 From: Jerome Kieffer Date: Tue, 19 Nov 2024 09:16:27 +0100 Subject: [PATCH 06/13] Correct calculation & fix test --- src/pyFAI/ext/CSR_common.pxi | 39 ++++++++++++++------------ src/pyFAI/test/test_medfilt_engine.py | 40 +++++++++++++++++++++++---- 2 files changed, 57 insertions(+), 22 deletions(-) diff --git a/src/pyFAI/ext/CSR_common.pxi b/src/pyFAI/ext/CSR_common.pxi index 19ac09980..9a83af563 100644 --- a/src/pyFAI/ext/CSR_common.pxi +++ b/src/pyFAI/ext/CSR_common.pxi @@ -29,7 +29,7 @@ __author__ = "Jérôme Kieffer" __contact__ = "Jerome.kieffer@esrf.fr" -__date__ = "15/11/2024" +__date__ = "19/11/2024" __status__ = "stable" __license__ = "MIT" @@ -787,23 +787,24 @@ cdef class CsrIntegrator(object): """ error_model = ErrorModel.parse(error_model) cdef: - index_t i, j, c, bad_pix, npix = self._indices.shape[0], idx = 0, start, stop + index_t i, j, c, bad_pix, npix = self._indices.shape[0], idx = 0, start, stop, cnt=0 acc_t acc_sig = 0.0, acc_var = 0.0, acc_norm = 0.0, acc_count = 0.0, coef = 0.0, acc_norm_sq=0.0 - acc_t cumsum, cnt = 0.0, qmin, qmax + acc_t cumsum = 0.0 + data_t qmin, qmax data_t empty, sig, var, nrm, weight, nrm2 - acc_t[::1] sum_sig = numpy.empty(self.output_size, dtype=acc_d) - acc_t[::1] sum_var = numpy.empty(self.output_size, dtype=acc_d) - acc_t[::1] sum_norm = numpy.empty(self.output_size, dtype=acc_d) - acc_t[::1] sum_norm_sq = numpy.empty(self.output_size, dtype=acc_d) - acc_t[::1] sum_count = numpy.empty(self.output_size, dtype=acc_d) - data_t[::1] merged = numpy.empty(self.output_size, dtype=data_d) - data_t[::1] stda = numpy.empty(self.output_size, dtype=data_d) - data_t[::1] sema = numpy.empty(self.output_size, dtype=data_d) + acc_t[::1] sum_sig = numpy.zeros(self.output_size, dtype=acc_d) + acc_t[::1] sum_var = numpy.zeros(self.output_size, dtype=acc_d) + acc_t[::1] sum_norm = numpy.zeros(self.output_size, dtype=acc_d) + acc_t[::1] sum_norm_sq = numpy.zeros(self.output_size, dtype=acc_d) + index_t[::1] sum_count = numpy.zeros(self.output_size, dtype=index_d) + data_t[::1] merged = numpy.zeros(self.output_size, dtype=data_d) + data_t[::1] stda = numpy.zeros(self.output_size, dtype=data_d) + data_t[::1] sema = numpy.zeros(self.output_size, dtype=data_d) data_t[:, ::1] preproc4 bint do_azimuthal_variance = error_model == ErrorModel.AZIMUTHAL bint do_hybrid_variance = error_model == ErrorModel.HYBRID float4_t element, former_element - float4_t[::1] work = numpy.empty(npix, dtype=float4_d) + float4_t[::1] work = numpy.zeros(npix, dtype=float4_d) assert weights.size == self.input_size, "weights size" empty = dummy if dummy is not None else self.empty @@ -828,11 +829,11 @@ cdef class CsrIntegrator(object): with nogil: # Duplicate the input data and populate the large work-array for i in range(npix): # NOT faster in parallel ! + weight = self._data[i] j = self._indices[i] sig = preproc4[j,0] var = preproc4[j,1] nrm = preproc4[j,2] - weight = self._data[i] element.s0 = sig/nrm # average signal element.s1 = sig * weight # weighted raw signal element.s2 = var * weight * weight # weighted raw variance @@ -842,12 +843,15 @@ cdef class CsrIntegrator(object): start = self._indptr[idx] stop = self._indptr[idx+1] acc_sig = acc_var = acc_norm = acc_norm_sq = 0.0 + cnt = 0 + cumsum = 0.0 + sort_float4(work[start:stop]) - cumsum = 0.0 for i in range(start, stop): cumsum = cumsum + work[i].s3 work[i].s0 = cumsum + qmin = quant_min * cumsum qmax = quant_max * cumsum @@ -855,12 +859,13 @@ cdef class CsrIntegrator(object): for i in range(start, stop): former_element = element element = work[i] - if (former_element.s0 >= qmin) and (element.s0 <= qmax): + if (qmin<=former_element.s0) and (element.s0 <= qmax): acc_sig = acc_sig + element.s1 acc_var = acc_var + element.s2 acc_norm = acc_norm + element.s3 acc_norm_sq = acc_norm_sq + element.s3*element.s3 - cnt = cnt + 1.0 + cnt = cnt + 1 + #collect things ... sum_sig[idx] = acc_sig sum_var[idx] = acc_var @@ -868,7 +873,7 @@ cdef class CsrIntegrator(object): sum_norm_sq[idx] = acc_norm_sq sum_count[idx] = cnt if (acc_norm_sq): - merged[idx] = acc_sig/acc_var + merged[idx] = acc_sig/acc_norm stda[idx] = sqrt(acc_var / acc_norm_sq) sema[idx] = sqrt(acc_var) / acc_norm else: diff --git a/src/pyFAI/test/test_medfilt_engine.py b/src/pyFAI/test/test_medfilt_engine.py index 4bc472d08..a594b2e15 100644 --- a/src/pyFAI/test/test_medfilt_engine.py +++ b/src/pyFAI/test/test_medfilt_engine.py @@ -32,7 +32,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "14/11/2024" +__date__ = "19/11/2024" import unittest import numpy @@ -50,7 +50,7 @@ class TestMedfilt(unittest.TestCase): @classmethod def setUpClass(cls)->None: super(TestMedfilt, cls).setUpClass() - cls.method = ["full", "csr", "python"] + cls.method = ("full", "csr", "python") cls.img = fabio.open(UtilsTest.getimage("mock.tif")).data cls.ai = load({ "dist": 0.1, "poni1":0.03, @@ -67,23 +67,53 @@ def tearDownClass(cls)->None: super(TestMedfilt, cls).tearDownClass() cls.method = cls.img =cls.ai =cls.npt =None - def test_python(self): - print(self.ai) + method = list(self.method) + method[-1] = "python" method = tuple(self.method) ref = self.ai.integrate1d(self.img, self.npt, unit="2th_rad", method=method, error_model="poisson") + # print(ref.method) + engine = self.ai.engines[ref.method].engine + obt = engine.medfilt(self.img, + solidangle=self.ai.solidAngleArray(), + quant_min=0,quant_max=1, # taking all Like this it works like a normal mean + error_model="poisson") + + self.assertTrue(numpy.allclose(ref.radial, obt.position), "radial matches") + self.assertTrue(numpy.allclose(ref.sum_signal, obt.signal), "signal matches") + self.assertTrue(numpy.allclose(ref.sum_variance, obt.variance), "variance matches") + self.assertTrue(numpy.allclose(ref.sum_normalization, obt.normalization), "normalization matches") + self.assertTrue(numpy.allclose(ref.sum_normalization2, obt.norm_sq), "norm_sq matches") + + self.assertTrue(numpy.allclose(ref.intensity, obt.intensity), "intensity matches") + self.assertTrue(numpy.allclose(ref.sigma, obt.sigma), "sigma matches") + self.assertTrue(numpy.allclose(ref.std, obt.std), "std matches") + self.assertTrue(numpy.allclose(ref.sem, obt.sem), "sem matches") + + def test_cython(self): + # print(self.ai) + + method = list(self.method) + # method[0] = "no" + method[-1] = "cython" + method = tuple(method) + ref = self.ai.integrate1d(self.img, self.npt, unit="2th_rad", method=method, error_model="poisson") print(ref.method) engine = self.ai.engines[ref.method].engine + print(engine) obt = engine.medfilt(self.img, solidangle=self.ai.solidAngleArray(), - quantile=(0,1), # taking all Like this it works like a normal mean + quant_min=0,quant_max=1, # taking all Like this it works like a normal mean error_model="poisson") + # print(ref.count-obt.count) + # print() self.assertTrue(numpy.allclose(ref.radial, obt.position), "radial matches") self.assertTrue(numpy.allclose(ref.sum_signal, obt.signal), "signal matches") self.assertTrue(numpy.allclose(ref.sum_variance, obt.variance), "variance matches") self.assertTrue(numpy.allclose(ref.sum_normalization, obt.normalization), "normalization matches") self.assertTrue(numpy.allclose(ref.sum_normalization2, obt.norm_sq), "norm_sq matches") + # self.assertTrue(numpy.allclose(ref.count, obt.count), "count matches") # not valid with pixel splitting self.assertTrue(numpy.allclose(ref.intensity, obt.intensity), "intensity matches") self.assertTrue(numpy.allclose(ref.sigma, obt.sigma), "sigma matches") From d278eb5a958a1b83470689ceea6716fd77528813 Mon Sep 17 00:00:00 2001 From: Jerome Kieffer Date: Tue, 19 Nov 2024 09:47:47 +0100 Subject: [PATCH 07/13] Issue with pre-commit ... --- src/pyFAI/engines/CSR_engine.py | 2 +- src/pyFAI/ext/splitPixelFullCSR.pyx | 2 +- src/pyFAI/test/test_all.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pyFAI/engines/CSR_engine.py b/src/pyFAI/engines/CSR_engine.py index 4a79814dd..5b30b61a9 100644 --- a/src/pyFAI/engines/CSR_engine.py +++ b/src/pyFAI/engines/CSR_engine.py @@ -26,7 +26,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "18/11/2024" +__date__ = "19/11/2024" __status__ = "development" from collections.abc import Iterable diff --git a/src/pyFAI/ext/splitPixelFullCSR.pyx b/src/pyFAI/ext/splitPixelFullCSR.pyx index c358f04f0..7148f1b26 100644 --- a/src/pyFAI/ext/splitPixelFullCSR.pyx +++ b/src/pyFAI/ext/splitPixelFullCSR.pyx @@ -35,7 +35,7 @@ Sparse matrix represented using the CompressedSparseRow. __author__ = "Jérôme Kieffer" __contact__ = "Jerome.kieffer@esrf.fr" -__date__ = "15/11/2024" +__date__ = "19/11/2024" __status__ = "stable" __license__ = "MIT" diff --git a/src/pyFAI/test/test_all.py b/src/pyFAI/test/test_all.py index d738acab3..9ccc5c896 100644 --- a/src/pyFAI/test/test_all.py +++ b/src/pyFAI/test/test_all.py @@ -32,7 +32,7 @@ __contact__ = "jerome.kieffer@esrf.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "14/11/2024" +__date__ = "19/11/2024" import sys import unittest From ba2d955226cae8ddd7bd027892654a0f6f11cd47 Mon Sep 17 00:00:00 2001 From: Jerome Kieffer Date: Tue, 19 Nov 2024 09:50:41 +0100 Subject: [PATCH 08/13] revert to previous version of pre-commit --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index abe722cf4..6415cc8d2 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -2,7 +2,7 @@ # See https://pre-commit.com/hooks.html for more hooks repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v5.0.0 + rev: v4.6.0 hooks: - id: trailing-whitespace - id: end-of-file-fixer From 0fd06d9b079f8579d4a8555713b0e4ef8d982f21 Mon Sep 17 00:00:00 2001 From: Jerome Kieffer Date: Tue, 19 Nov 2024 10:28:53 +0100 Subject: [PATCH 09/13] precommit issue --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6415cc8d2..abe722cf4 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -2,7 +2,7 @@ # See https://pre-commit.com/hooks.html for more hooks repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v4.6.0 + rev: v5.0.0 hooks: - id: trailing-whitespace - id: end-of-file-fixer From e15a3e407c27df358d9c63a760717a74121327eb Mon Sep 17 00:00:00 2001 From: Jerome Kieffer Date: Tue, 19 Nov 2024 10:35:17 +0100 Subject: [PATCH 10/13] work on pre-commit --- src/pyFAI/containers.py | 2 +- src/pyFAI/ext/sparse_utils.pyx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pyFAI/containers.py b/src/pyFAI/containers.py index b39887b49..265254da8 100644 --- a/src/pyFAI/containers.py +++ b/src/pyFAI/containers.py @@ -30,7 +30,7 @@ __contact__ = "valentin.valls@esrf.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "12/06/2024" +__date__ = "19/11/2024" __status__ = "development" from collections import namedtuple diff --git a/src/pyFAI/ext/sparse_utils.pyx b/src/pyFAI/ext/sparse_utils.pyx index 837eae2d3..708ebd37f 100644 --- a/src/pyFAI/ext/sparse_utils.pyx +++ b/src/pyFAI/ext/sparse_utils.pyx @@ -33,7 +33,7 @@ __author__ = "Jerome Kieffer" __contact__ = "Jerome.kieffer@esrf.fr" -__date__ = "21/08/2024" +__date__ = "19/11/2024" __status__ = "stable" __license__ = "MIT" From cd2343ae661e100ef3e08bb6987119cd2ebb9b0c Mon Sep 17 00:00:00 2001 From: Jerome Kieffer Date: Tue, 19 Nov 2024 10:38:35 +0100 Subject: [PATCH 11/13] fix trailing spaces in doc --- doc/source/publications.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/source/publications.rst b/doc/source/publications.rst index e77f2464e..6ff386de2 100644 --- a/doc/source/publications.rst +++ b/doc/source/publications.rst @@ -1,5 +1,5 @@ :Author: Jérôme Kieffer -:Date: 07/01/2021 +:Date: 19/11/2024 :Keywords: List of publications @@ -34,11 +34,11 @@ Publications about pyFAI * *Application of signal separation to diffraction image compression and serial crystallography*; Jérôme Kieffer, Julien Orlans, Nicolas Coquelle, Samuel Debionne, Shibom Basu, Alejandro Homs, Gianluca Santonia and Daniele De Sanctis; - `Accepted `_ in **J. Applied Crystallography** (2024); + `Accepted `_ in **J. Applied Crystallography** (2024); In depth explainaion of sigma-clipping background assessment and error models. The latest paper should be the cited in publications using pyFAI. There are already 1400 publications referring to pyFAI, some of them in the most -prestigious scientific journals (Nature, PNAS, ...) and -40 other `applications `_ +prestigious scientific journals (Nature, PNAS, ...) and +40 other `applications `_ using pyFAI as a library. From 8ac0124fbb1fb7819982cd162e21b6fe4881142a Mon Sep 17 00:00:00 2001 From: Edgar Gutierrez <112327909+EdgarGF93@users.noreply.github.com> Date: Tue, 19 Nov 2024 10:49:23 +0100 Subject: [PATCH 12/13] regression --- doc/source/publications.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/publications.rst b/doc/source/publications.rst index 6ff386de2..87e2acb10 100644 --- a/doc/source/publications.rst +++ b/doc/source/publications.rst @@ -1,5 +1,5 @@ :Author: Jérôme Kieffer -:Date: 19/11/2024 +:Date: 19/11/2024 :Keywords: List of publications From 7e313f5e1c85ab835127dd361b9d1e26662fd881 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 19 Nov 2024 09:49:51 +0000 Subject: [PATCH 13/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- doc/source/publications.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/publications.rst b/doc/source/publications.rst index 87e2acb10..6ff386de2 100644 --- a/doc/source/publications.rst +++ b/doc/source/publications.rst @@ -1,5 +1,5 @@ :Author: Jérôme Kieffer -:Date: 19/11/2024 +:Date: 19/11/2024 :Keywords: List of publications