Merge pull request #46 from peter-madigan/upgrade/implement-unique-ch…

…annel-thresholds

Revert "Revert "Upgrade/implement unique channel thresholds""

.gitattributes

@@ -1,2 +1,3 @@
 *.h5 filter=lfs diff=lfs merge=lfs -text
 *.npy filter=lfs diff=lfs merge=lfs -text
+*.npz filter=lfs diff=lfs merge=lfs -text

cli/simulate_pixels.py

@@ -15,6 +15,7 @@
 from tqdm import tqdm
 
 from larndsim import consts
+from larndsim.cuda_dict import CudaDict
 
 LOGO = """
   _                      _            _
@@ -32,7 +33,8 @@ def run_simulation(input_filename,
                    output_filename='',
                    response_file='../larndsim/response_44.npy',
                    bad_channels=None,
-                   n_tracks=None):
+                   n_tracks=None,
+                   pixel_thresholds_file=None):
     """
     Command-line interface to run the simulation of a pixelated LArTPC
 
@@ -42,6 +44,7 @@ def run_simulation(input_filename,
             layout and connection details.
         detector_properties (str): path of the YAML file containing
             the detector properties
+        pixel_thresholds_file (str): path to npz file containing pixel thresholds
         output_filename (str): path of the HDF5 output file. If not specified
             the output is added to the input file.
         response_file: path of the Numpy array containing the pre-calculated
@@ -74,6 +77,14 @@ def run_simulation(input_filename,
     from larndsim import quenching, drifting, detsim, pixels_from_track, fee
     RangePop()
 
+    RangePush("load_pixel_thresholds")
+    if pixel_thresholds_file is not None:
+        print("Pixel thresholds file:", pixel_thresholds_file)
+        pixel_thresholds_lut = CudaDict.load(pixel_thresholds_file, 256)
+    else:
+        pixel_thresholds_lut = CudaDict(cp.array([fee.DISCRIMINATION_THRESHOLD]), 1, 1)
+    RangePop()
+
     RangePush("load_hd5_file")
     # First of all we load the edep-sim output
     # For this sample we need to invert $z$ and $y$ axes
@@ -207,9 +218,9 @@ def run_simulation(input_filename,
 
             RangePush("unique_pix")
             shapes = neighboring_pixels.shape
-            joined = neighboring_pixels.reshape(shapes[0]*shapes[1])
+            joined = neighboring_pixels.reshape(shapes[0] * shapes[1])
             unique_pix = cp.unique(joined)
-            unique_pix = unique_pix[(unique_pix != -1),:]
+            unique_pix = unique_pix[(unique_pix != -1)]
             RangePop()
 
             if not unique_pix.shape[0]:
@@ -282,15 +293,20 @@ def run_simulation(input_filename,
 
             current_fractions = cp.zeros((pixels_signals.shape[0], fee.MAX_ADC_VALUES, track_pixel_map.shape[1]))
             rng_states = create_xoroshiro128p_states(TPB * BPG, seed=ievd)
-
-            fee.get_adc_values[BPG,TPB](pixels_signals,
-                                        pixels_tracks_signals,
-                                        time_ticks,
-                                        integral_list,
-                                        adc_ticks_list,
-                                        0,
-                                        rng_states,
-                                        current_fractions)
+            pixel_thresholds_lut.tpb = TPB
+            pixel_thresholds_lut.bpg = BPG
+            orig_shape = unique_pix.shape
+            pixel_thresholds = pixel_thresholds_lut[unique_pix.ravel()].reshape(orig_shape)
+
+            fee.get_adc_values[BPG, TPB](pixels_signals,
+                                         pixels_tracks_signals,
+                                         time_ticks,
+                                         integral_list,
+                                         adc_ticks_list,
+                                         0,
+                                         rng_states,
+                                         current_fractions,
+                                         pixel_thresholds)
 
             adc_list = fee.digitize(integral_list)
             adc_event_ids = np.full(adc_list.shape, unique_eventIDs[0]) # FIXME: only works if looping on a single event

larndsim/cuda_dict.py

@@ -0,0 +1,216 @@
+from numba import cuda
+import cupy as cp
+from math import ceil
+
+
+_EMPTY_KEY = 2**31-1
+
+
+class CudaDict(object):
+    '''
+        A numba implementation of a static hash table that lives on the GPU.
+        Based on this project
+        https://github.com/nosferalatu/SimpleGPUHashTable.git
+
+        Initialization is preformed via::
+
+            cd = CudaDict(default=cp.array([0.]), tpb=256, bpg=256)
+            cd[keys] = values
+
+        Lookup is available via::
+
+            values = cd[keys]
+
+        ``keys`` are expected to be a 1D array of integer key values. When
+        initializing, ``keys`` must be unique.
+
+        Removal of items can be performed via::
+
+            del cd[keys]
+
+    '''
+
+    def __init__(self, default, tpb, bpg):
+        self.tpb = tpb
+        self.bpg = bpg
+        self.default = default
+        self._hashtable_keys = cp.full(1, _EMPTY_KEY, dtype=cp.int32)
+        self._hashtable_values = cp.full(1, default, dtype=default.dtype)
+
+    def keys(self):
+        mask = self._hashtable_keys != _EMPTY_KEY
+        return self._hashtable_keys[mask]
+
+    def values(self):
+        mask = self._hashtable_keys != _EMPTY_KEY
+        return self._hashtable_values[mask]
+
+    def items(self):
+        return self.keys(), self.values()
+
+    def __getitem__(self, key):
+        values = cp.full(key.shape[0], self.default, dtype=self._hashtable_values.dtype)
+        cuda_hashtable_lookup[self.tpb, self.bpg](
+            key, values, self._hashtable_keys, self._hashtable_values, self.default)
+        return values
+
+    def __setitem__(self, key, value):
+        if len(self) == 0:
+            self._hashtable_keys = cp.full(key.shape[0] + 1, _EMPTY_KEY, dtype=cp.int32)
+            self._hashtable_values = cp.full(key.shape[0] + 1, self.default, dtype=value.dtype)
+        else:
+            raise NotImplementedError('Trying to update CudaDict, not yet supported')
+        cuda_hashtable_insert[self.tpb, self.bpg](
+            key, value, self._hashtable_keys, self._hashtable_values)
+
+    def __delitem__(self, key):
+        if len(self) == 0:
+            return
+        cuda_hashtable_delete[self.tpb, self.bpg](
+            key, self._hashtable_key_arr, self._hashtable_value_arr)
+
+    def __len__(self):
+        return len(self._hashtable_keys) - 1
+
+    def contains(self, key):
+        exists = cp.zeros(key.shape[0], dtype=bool)
+        if len(self) == 0:
+            return exists
+        cuda_hashtable_exists[self.tpb, self.bpg](
+            key, exists, self._hashtable_keys)
+        return exists
+
+    @staticmethod
+    def load(filename, tpb):
+        data = cp.load(filename)
+        keys = data['keys']
+        values = data['values']
+        default = data['default']
+        bpg = ceil(len(keys)/tpb)
+        cd = CudaDict(default=default, tpb=tpb, bpg=bpg)
+        cd[keys] = values
+        return cd
+
+    @staticmethod
+    def save(filename, cdict):
+        mask = cdict._hashtable_keys != _EMPTY_KEY
+        keys = cdict._hashtable_keys[mask]
+        values = cdict._hashtable_values[mask]
+        default = cdict.default
+        data = dict(
+            keys=keys,
+            values=values,
+            default=default
+        )
+        cp.savez_compressed(filename, **data)
+
+
+_HASH_CONSTANT = 257
+
+
+@cuda.jit(device=True)
+def cuda_hashtable_encode(key_arr, entries):
+    '''
+    Encodes keys into a integer number between 0 and entries-1
+    Args:
+        key_arr (: obj: `numpy.ndarray`): list of keys
+        entries (int): maximum integer value to return
+    '''
+    return ((key_arr * _HASH_CONSTANT) % entries)
+
+
+@cuda.jit
+def cuda_hashtable_insert(key_arr, value_arr, hashtable_key_arr, hashtable_value_arr):
+    '''
+    Inserts keys into a hashtable
+    Args:
+        key_arr (: obj: `numpy.ndarray`): list of keys to insert
+        value_arr (: obj: `numpy.ndarray`): list of values to insert
+        hashtable_key_arr (: obj: `numpy.ndarray`): list of unique keys in hashtable
+        hashtable_value_arr (: obj: `numpy.ndarray`): list of values in hashtable
+    '''
+    ikey = cuda.grid(1)
+
+    if ikey < key_arr.shape[0]:
+
+        encoding = cuda_hashtable_encode(key_arr[ikey], hashtable_key_arr.shape[0])
+        while True:
+            cuda.atomic.compare_and_swap(hashtable_key_arr[encoding:encoding+1], _EMPTY_KEY, key_arr[ikey])
+            if hashtable_key_arr[encoding] == key_arr[ikey]:
+                hashtable_value_arr[encoding] = value_arr[ikey]
+                break
+            encoding = (encoding + 1) % hashtable_key_arr.shape[0]
+
+
+@cuda.jit
+def cuda_hashtable_lookup(key_arr, value_arr, hashtable_key_arr, hashtable_value_arr,
+                          default):
+    '''
+    Fetches values from hashtable
+    Args:
+        key_arr (: obj: `numpy.ndarray`): list of keys to lookup
+        value_arr (: obj: `numpy.ndarray`): output array
+        hashtable_key_arr (: obj: `numpy.ndarray`): list of unique keys in hashtable
+        hashtable_value_arr (: obj: `numpy.ndarray`): list of values in hashtable
+        default (: obj: `numpy.ndarray`): singleton array to use if key is not present
+    '''
+    ikey = cuda.grid(1)
+
+    if ikey < key_arr.shape[0]:
+
+        encoding = cuda_hashtable_encode(key_arr[ikey], hashtable_key_arr.shape[0])
+        while True:
+            if hashtable_key_arr[encoding] == key_arr[ikey]:
+                value_arr[ikey] = hashtable_value_arr[encoding]
+                break
+            elif hashtable_key_arr[encoding] == _EMPTY_KEY:
+                value_arr[ikey] = default[0]
+                break
+            encoding = (encoding + 1) % hashtable_key_arr.shape[0]
+
+
+@cuda.jit
+def cuda_hashtable_exists(key_arr, exists_arr, hashtable_key_arr):
+    '''
+    Checks if a key is present in the hashtable
+    Args:
+        key_arr (: obj: `numpy.ndarray`): list of keys to remove
+        exists_arr (: obj: `numpy.ndarray`): output array
+        hashtable_key_arr (: obj: `numpy.ndarray`): list of unique keys in hashtable
+    '''
+    ikey = cuda.grid(1)
+
+    if ikey < key_arr.shape[0]:
+
+        encoding = cuda_hashtable_encode(key_arr[ikey], hashtable_key_arr.shape[0])
+        while True:
+            if hashtable_key_arr[encoding] == key_arr[ikey]:
+                exists_arr[ikey] = True
+                break
+            if hashtable_key_arr[encoding] == _EMPTY_KEY:
+                exists_arr[ikey] = False
+                break
+            encoding = (encoding + 1) % hashtable_key_arr.shape[0]
+
+
+@cuda.jit
+def cuda_hashtable_delete(key_arr, hashtable_key_arr, hashtable_value_arr):
+    '''
+    Removes a key from the hashtable
+    Args:
+        key_arr (: obj: `numpy.ndarray`): list of keys to remove
+        hashtable_key_arr (: obj: `numpy.ndarray`): list of unique keys in hashtable
+        hashtable_value_arr (: obj: `numpy.ndarray`): list of values in hashtable
+    '''
+    ikey = cuda.grid(1)
+
+    if ikey < key_arr.shape[0]:
+
+        encoding = cuda_hashtable_encode(key_arr[ikey], hashtable_key_arr.shape[0])
+        while True:
+            if hashtable_key_arr[encoding] == key_arr[ikey]:
+                hashtable_key_arr[encoding] == _EMPTY_KEY
+                break
+            if hashtable_key_arr[encoding] == _EMPTY_KEY:
+                break
+            encoding = (encoding + 1) % hashtable_key_arr.shape[0]

larndsim/detector_properties/module0-far_field.yaml

@@ -10,5 +10,7 @@ response_sampling: 0.1 # us
 tile_map: 
   - [[7,5,3,1],[8,6,4,2]]
   - [[16,14,12,10],[15,13,11,9]]
+module_to_io_groups:
+  1: [1,2,3,4]
 time_padding: 190 # us
 time_window: 189.1 # us

larndsim/detsim.py

@@ -232,11 +232,11 @@ def get_pixel_coordinates(pixel_id):
     """
     Returns the coordinates of the pixel center given the pixel ID
     """
-    px, py, plane_id = id2pixel(pixel_id)
+    i_x, i_y, plane_id = id2pixel(pixel_id)
 
     this_border = tpc_borders[int(plane_id)]
-    pix_x = px * consts.pixel_pitch + this_border[0][0]
-    pix_y = py * consts.pixel_pitch + this_border[1][0]
+    pix_x = i_x * consts.pixel_pitch + this_border[0][0]
+    pix_y = i_y * consts.pixel_pitch + this_border[1][0]
 
     return pix_x,pix_y
 
@@ -276,9 +276,8 @@ def tracks_current(signals, pixels, tracks, response):
         signals (:obj:`numpy.ndarray`): empty 3D array with dimensions S x P x T,
             where S is the number of track segments, P is the number of pixels, and T is
             the number of time ticks. The output is stored here.
-        pixels (:obj:`numpy.ndarray`): 3D array with dimensions S x P x 2, where S is
-            the number of track segments, P is the number of pixels and the third dimension
-            contains the two pixel ID numbers.
+        pixels (:obj:`numpy.ndarray`): 2D array with dimensions S x P , where S is
+            the number of track segments, P is the number of pixels and contains the pixel ID number.
         tracks (:obj:`numpy.ndarray`): 2D array containing the detector segments.
         response (:obj:`numpy.ndarray`): 3D array containing the tabulated response.
     """
@@ -288,8 +287,9 @@ def tracks_current(signals, pixels, tracks, response):
     if itrk < signals.shape[0] and ipix < signals.shape[1] and it < signals.shape[2]:
         t = tracks[itrk]
         pID = pixels[itrk][ipix]
+        pID_x, pID_y, pID_plane = id2pixel(pID)
 
-        if pID >= 0:
+        if pID_x >= 0 and pID_y >= 0:
 
             # Pixel coordinates
             x_p, y_p = get_pixel_coordinates(pID)