Disable melt_snow by default

This is unreliable, and it's disabled in master.

cal_tools/cal_tools/enums.py

@@ -26,4 +26,11 @@ class BadPixels(Enum):
     OVERSCAN                 = 0b001000000000000000000 # bit 19
     NON_SENSITIVE            = 0b010000000000000000000 # bit 20
     NON_LIN_RESPONSE_REGION  = 0b100000000000000000000 # bit 21
+    
+    
+class SnowResolution(Enum):
+    """ An Enum specifying how to resolve snowy pixels
+    """
+    NONE = "none"
+    INTERPOLATE = "interpolate"
 

cal_tools/cal_tools/step_timing.py

+from collections import defaultdict
+from time import perf_counter
+
+import numpy as np
+
+class StepTimer:
+    def __init__(self):
+        self.steps = defaultdict(list)
+        self.t0 = None
+        self.t_latest = None
+
+    def start(self):
+        """Store a start timestamp"""
+        t = perf_counter()
+        if self.t0 is None:
+            self.t0 = t
+        self.t_latest = t
+
+    def done_step(self, step_name):
+        """Record a step timing, since the last call to done_step() or start()
+        """
+        t = perf_counter()
+        self.steps[step_name].append(t - self.t_latest)
+        print(f"{step_name}: {t - self.t_latest:.01f} s")
+        self.t_latest = t
+
+    def timespan(self):
+        """Wall time from 1st call to start() to last start() or done_step()"""
+        return self.t_latest - self.t0
+
+    def print_summary(self):
+        """Show mean & std for each step"""
+        for step, data in self.steps.items():
+            data = np.asarray(data)
+            print(f'{step}: {data.mean():.01f} +- {data.std():.02f}s')

cal_tools/cython/agipdalgs.pyx

+import numpy as np
+cimport numpy as cnp
+cimport cython
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def histogram(cnp.ndarray[cnp.float32_t, ndim=2] data, range=(0,1), int bins=20,
+           cnp.ndarray[cnp.float32_t, ndim=2] weights=None):
+    """
+    Calculate N histograms along axis 0 with equidistant binning.
+    N is a data size along axis 1.
+    Return histograms and bin edges.
+    """
+
+    cdef cnp.ndarray[cnp.float32_t, ndim=2] ret
+    cdef double min, max 
+    min = range[0] 
+    max = range[1] 
+         
+    ret = np.zeros((bins,data.shape[1]), dtype=np.float32)
+    cdef double bin_width = (max - min) / bins
+    cdef double x
+
+    if weights is not None:
+        assert (<object>weights).shape == (<object>data).shape,\
+            "data & weights must have matching shape"
+
+    for j in xrange(data.shape[1]):
+        for i in xrange(data.shape[0]):
+            x = (data[i,j] - min) / bin_width
+            if 0.0 <= x < bins:
+                if weights is None: 
+                    ret[<int>x,j] += 1.0
+                else: 
+                    ret[<int>x,j] += weights[i,j]
+    return ret, np.linspace(min, max, bins+1)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def histogram2d(cnp.ndarray[float, ndim=1] data_x, cnp.ndarray[float, ndim=1] data_y,
+                range=((0,1), (0,1)), bins=(20, 20)):
+    """
+    Calculate 2D histogram with equidistant binning.
+    Return histogram and bin edges.
+    """
+
+    cdef cnp.ndarray[cnp.int32_t, ndim=2] ret
+    cdef double min_x, max_x, min_y, max_y
+    cdef int bins_x, bins_y
+    min_x = range[0][0]
+    max_x = range[0][1]
+    min_y = range[1][0]
+    max_y = range[1][1]
+    bins_x = bins[0]
+    bins_y = bins[1]
+
+    ret = np.zeros((bins_x, bins_y), dtype=np.int32)
+
+    cdef double bin_width_x = (max_x - min_x) / bins_x
+    cdef double bin_width_y = (max_y - min_y) / bins_y
+    cdef double x, y
+
+    assert (<object>data_x).shape == (<object>data_y).shape,\
+        "data_x & data_y must have matching shape"
+
+    for i in xrange(data_x.shape[0]):
+        x = (data_x[i] - min_x) / bin_width_x
+        y = (data_y[i] - min_y) / bin_width_y
+        if x >= 0.0 and x < bins_x and y>=0.0 and y<bins_y:
+            ret[<int>x, <int>y] += 1
+
+    return ret, np.linspace(min_x, max_x, bins_x+1), np.linspace(min_y, max_y, bins_y+1)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def gain_choose(cnp.ndarray[cnp.uint8_t, ndim=3] a, cnp.ndarray[cnp.float32_t, ndim=4] choices):
+    """Specialised fast equivalent of np.choose(), to select data for a per-pixel gain stage"""
+    cdef int i, j, k
+    cdef cnp.uint8_t v
+    cdef cnp.ndarray[cnp.float32_t, ndim=3] out
+    out = np.zeros_like(a, dtype=np.float32)
+    
+    assert (<object>choices).shape == (3,) + (<object>a).shape
+    
+    with nogil:
+        for i in range(a.shape[0]):
+            for j in range(a.shape[1]):
+                for k in range(a.shape[2]):
+                    v = a[i, j, k]
+                    out[i, j, k] = choices[v, i, j, k]
+    
+    return out
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def gain_choose_int(cnp.ndarray[cnp.uint8_t, ndim=3] a, cnp.ndarray[cnp.int32_t, ndim=4] choices):
+    """Specialised fast equivalent of np.choose(), to select data for a per-pixel gain stage"""
+    cdef int i, j, k
+    cdef cnp.uint8_t v
+    cdef cnp.ndarray[cnp.int32_t, ndim=3] out
+    out = np.zeros_like(a, dtype=np.int32)
+    
+    assert (<object>choices).shape == (3,) + (<object>a).shape
+    
+    with nogil:
+        for i in range(a.shape[0]):
+            for j in range(a.shape[1]):
+                for k in range(a.shape[2]):
+                    v = a[i, j, k]
+                    out[i, j, k] = choices[v, i, j, k]
+    
+    return out
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def sum_and_count_in_range_asic(cnp.ndarray[float, ndim=4] arr, float lower, float upper):
+    """
+    Return the sum & count of values where lower <= x <= upper,
+    across axes 2 & 3 (pixels within an ASIC, as reshaped by AGIPD correction code).
+    Specialised function for performance.
+    """
+
+    cdef int i, j, k, l, m
+    cdef float value
+    cdef cnp.ndarray[unsigned long long, ndim=2] count
+    cdef cnp.ndarray[double, ndim=2] sum_
+    
+    # Drop axes -2 & -1 (pixel dimensions within each ASIC)
+    out_shape = arr[:, :, 0, 0].shape
+    count = np.zeros(out_shape, dtype=np.uint64)
+    sum_ = np.zeros(out_shape, dtype=np.float64)
+    
+    with nogil:
+        for i in range(arr.shape[0]):
+            for k in range(arr.shape[1]):
+                for l in range(arr.shape[2]):
+                    for m in range(arr.shape[3]):
+                        value = arr[i, k, l, m]
+                        if lower <= value <= upper:
+                            sum_[i, k] += value
+                            count[i, k] += 1
+    return sum_, count
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def sum_and_count_in_range_cell(cnp.ndarray[float, ndim=4] arr, float lower, float upper):
+    """
+    Return the sum & count of values where lower <= x <= upper,
+    across axes 0 & 1 (memory cells in the same row, as reshaped by AGIPD correction code).
+    Specialised function for performance.
+    """
+
+    cdef int i, j, k, l, m,
+    cdef float value
+    cdef cnp.ndarray[unsigned long long, ndim=2] count
+    cdef cnp.ndarray[double, ndim=2] sum_
+    
+    # Drop axes 0 & 1
+    out_shape = arr[0, 0, :, :].shape
+    count = np.zeros(out_shape, dtype=np.uint64)
+    sum_ = np.zeros(out_shape, dtype=np.float64)
+    
+    
+    with nogil:
+        for i in range(arr.shape[0]):
+            for k in range(arr.shape[1]):
+                for l in range(arr.shape[2]):
+                    for m in range(arr.shape[3]):
+                        value = arr[i, k, l, m]
+                        if lower <= value <= upper:
+                            sum_[l, m] += value
+                            count[l, m] += 1
+
+    return sum_, count

requirements.txt

@@ -4,6 +4,7 @@ git+file:///gpfs/exfel/sw/calsoft/git/pyDetLib@2.5.3-2.7.0#subdirectory=lib
 astcheck == 0.2.5
 astsearch == 0.1.3
 dill == 0.3.0
+extra_data == 1.2.0
 extra_geom == 0.8.0
 fabio == 0.9.0
 gitpython == 3.1.0
@@ -28,6 +29,7 @@ python-dateutil == 2.8.1
 pyyaml == 5.3
 pyzmq == 19.0.0
 scikit-learn == 0.22.2.post1
+sharedmem == 0.3.7 
 sphinx == 1.4.5
 tabulate == 0.8.6
 unittest-xml-reporting == 3.0.2

setup.py

@@ -2,8 +2,18 @@ from setuptools import setup
 from setuptools.command.install import install
 from subprocess import check_call, check_output
 from distutils.command.build import build
+
+from Cython.Distutils import build_ext
+from distutils.extension import Extension
 import sys
+import numpy
 
+extensions = [Extension("cal_tools.cython.agipdalgs",
+                        ['cal_tools/cython/agipdalgs.pyx'],
+                        include_dirs=[numpy.get_include()],
+                        extra_compile_args=['-fopenmp', '-march=native'],
+                        extra_link_args=['-fopenmp'], ),
+              ]
 
 class PostInstallCommand(install):
     """Post-installation for installation mode."""
@@ -61,6 +71,7 @@ setup(
     cmdclass={
         'build' : PreInstallCommand,
         'install': PostInstallCommand,
+        'build_ext': build_ext
     },
     url='',
     license='(c) European XFEL GmbH 2018',
@@ -72,6 +83,7 @@ setup(
                   'xfel-calibrate = xfel_calibrate.calibrate:run',                  
               ],              
           },
+    ext_modules=extensions
     
 )