Perform fancy indexing on csr_matrices

Rather than expanding the full csr_matrix into a numpy array, then
performing fancy indexing on that, do the reverse: perform the
fancy indexing on the csr_matrix, and then expand it into a numpy
array.

This has the potential to save huge amounts of computing time, since
we do not have to create the entire dense array before performing
the indexing.

For the multiruby test case, this did not show any change to the
performance. However, for a monolithic Eiger example from Josh Ward,
it showed a 5x speed improvement.

Signed-off-by: Patrick Avery <[email protected]>

Author: psavery

hexrd/imageseries/load/framecache.py

@@ -209,7 +209,14 @@ def _load_framelist_if_needed(self):
 
     def __getitem__(self, key):
         self._load_framelist_if_needed()
-        return self._framelist[key].toarray()
+        if isinstance(key, tuple):
+            # Extract only what we need from the sparse array
+            # using fancy indexing before we convert it to a
+            # numpy array.
+            sparse_frame = self._framelist[key[0]][*key[1:]]
+        else:
+            sparse_frame = self._framelist[key]
+        return sparse_frame.toarray()
 
     def __iter__(self):
         return ImageSeriesIterator(self)

hexrd/imageseries/load/hdf5.py

@@ -64,6 +64,11 @@ def __del__(self):
             warnings.warn("HDF5ImageSeries could not close h5 file")
 
     def __getitem__(self, key):
+        if isinstance(key, tuple):
+            # FIXME: we do not yet support fancy indexing here.
+            # Fully expand the array then apply the fancy indexing.
+            return self[key[0]][*key[1:]]
+
         if self._ndim == 2:
             if key != 0:
                 raise IndexError(

hexrd/instrument/hedm_instrument.py

@@ -1723,7 +1723,7 @@ def pull_spots(self, plane_data, grain_params,
                         contains_signal = False
                         for i_frame in frame_indices:
                             contains_signal = contains_signal or np.any(
-                                ome_imgser[i_frame][ii, jj] > threshold
+                                ome_imgser[i_frame, ii, jj] > threshold
                             )
                         compl.append(contains_signal)
                         patch_output.append((ii, jj, frame_indices))
@@ -1792,7 +1792,7 @@ def pull_spots(self, plane_data, grain_params,
                         contains_signal = False
                         patch_data_raw = []
                         for i_frame in frame_indices:
-                            tmp = ome_imgser[i_frame][ijs[0], ijs[1]]
+                            tmp = ome_imgser[i_frame, ijs[0], ijs[1]]
                             contains_signal = contains_signal or np.any(
                                 tmp > threshold
                             )