Merging master into jtb (documentation branch)

.gitignore

@@ -1,4 +1,6 @@
+# ignore jupyter temp files
 .ipynb_checkpoints/
 Controller_Plans.ipynb
-./examples/.ipynb_checkpoints/
 XSpect_XES_examples.ipynb
+*.swp
+XSpect/*.swp

.readthedocs.yaml

@@ -0,0 +1,31 @@
+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Set the OS, Python version and other tools you might need
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.12"
+    # You can also specify other tool versions:
+    # nodejs: "19"
+    # rust: "1.64"
+    # golang: "1.19"
+
+# Build documentation in the "docs/" directory with Sphinx
+sphinx:
+  configuration: docs/conf.py
+# Optionally build your docs in additional formats such as PDF and ePub
+# formats:
+#    - pdf
+#    - epub
+
+# Optional but recommended, declare the Python requirements required
+# to build your documentation
+# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html
+python:
+    install:
+    - requirements: docs/requirements.txt

README.md

@@ -5,7 +5,7 @@
 Clone the repository:
 
    ```bash
-   git clone https://github.com/lg345/XSpecT.git
+   git clone https://github.com/lg345/XSpect.git
    ```
 
 Ensure that you have the necessary dependencies installed on your system. 

XSpect/.ipynb_checkpoints/XSpect_Analysis-checkpoint.py

@@ -263,6 +263,7 @@ def purge_all_keys(self,keys_to_keep):
             List of keys to retain.
         """
 
+        keys_to_keep = set(keys_to_keep)  # Remove duplicates by converting to a set
         new_dict = {attr: value for attr, value in self.__dict__.items() if attr in keys_to_keep}
         self.__dict__ = new_dict
 
@@ -397,11 +398,15 @@ def purge_keys(self,run,keys):
             setattr(run, detector_key, None)
             run.update_status(f"Purged key to save room: {detector_key}")
 
-    def reduce_detector_shots(self, run, detector_key,reduction_function=np.sum,  purge=True):
+    def reduce_detector_shots(self, run, detector_key,reduction_function=np.sum,  purge=True,new_key=False):
         detector = getattr(run, detector_key)
         reduced_data=reduction_function(detector,axis=0)
         run.update_status(f"Reduced detector by shots: {detector_key} with number of shots: {np.shape(detector)}")
-        setattr(run, f"{detector_key}_summed", reduced_data)
+        if new_key:
+            target_key=f"{detector_key}_summed"
+        else:
+            target_key=detector_key
+        setattr(run, target_key, reduced_data)
         if purge:
             setattr(run, detector_key,None)
             run.update_status(f"Purged key to save room: {detector_key}")
@@ -431,11 +436,16 @@ def reduce_detector_spatial(self, run, detector_key, shot_range=[0, None], rois=
         if combine:
 
             roi_combined = [rois[0][0], rois[-1][1]]  # Combined ROI spanning the first and last ROI
-            mask = np.zeros(detector.shape[2], dtype=bool)
+            mask = np.zeros(detector.shape[-1], dtype=bool)
             for roi in rois:
                 mask[roi[0]:roi[1]] = True
-            masked_data = detector[shot_range[0]:shot_range[1], :, :][:, :, mask]
-            reduced_data = reduction_function(masked_data, axis=2)
+            if detector.ndim==3:
+                masked_data = detector[shot_range[0]:shot_range[1], :, :][:, :, mask]
+            elif detector.ndim==2:
+                masked_data = detector[:, mask]
+            elif detector.ndim==1:
+                masked_data = detector[mask]
+            reduced_data = reduction_function(masked_data, axis=-1)
             roi_indices = ', '.join([f"{roi[0]}-{roi[1]}" for roi in rois])
             run.update_status(f"Spatially reduced detector: {detector_key} with combined ROI indices: {roi_indices}")
             setattr(run, f"{detector_key}_ROI_1", reduced_data)
@@ -452,7 +462,7 @@ def reduce_detector_spatial(self, run, detector_key, shot_range=[0, None], rois=
             setattr(run, detector_key,None)
             #delattr(run, detector_key)
             #del run.detector_key
-            run.update_status(f"Purged key to save room: {detector_key}")
+            run.update_status(f"Purged key after spatial reduction to save room: {detector_key}")
 
     def time_binning(self,run,bins,lxt_key='lxt_ttc',fast_delay_key='encoder',tt_correction_key='time_tool_correction'):
         """
@@ -471,10 +481,16 @@ def time_binning(self,run,bins,lxt_key='lxt_ttc',fast_delay_key='encoder',tt_cor
         tt_correction_key : str, optional
             The key for the time tool correction data (default is 'time_tool_correction').
         """
+        #print( str(getattr(run,fast_delay_key)))
+        #print( str(getattr(run,tt_correction_key)))
+        #print( str(getattr(run,fast_delay_key)+getattr(run,tt_correction_key)))
         if lxt_key==None:
-            run.delays = 0+ getattr(run,fast_delay_key)  + getattr(run,tt_correction_key)
+            delays = np.array(getattr(run,fast_delay_key)).flatten()  + np.array(getattr(run,tt_correction_key)).flatten()
         else:
-            run.delays = getattr(run,lxt_key)*1.0e12 + getattr(run,fast_delay_key)  + getattr(run,tt_correction_key)
+            delays = np.array(getattr(run,lxt_key)).flatten()*1.0e12 + np.array((getattr(run,fast_delay_key))).flatten()  + np.array(getattr(run,tt_correction_key)).flatten()
+        delays=np.array(getattr(run,fast_delay_key)).flatten()*1.0e12+np.array(getattr(run,tt_correction_key)).flatten()
+        #print(str(delays)) 
+        run.delays=delays
         run.time_bins=bins
         run.timing_bin_indices=np.digitize(run.delays, bins)[:]
         run.update_status('Generated timing bins from %f to %f in %d steps.' % (np.min(bins),np.max(bins),len(bins)))
@@ -550,7 +566,6 @@ def reduce_detector_temporal(self, run, detector_key, timing_bin_key_indices,ave
         """
         detector = getattr(run, detector_key)
         indices = getattr(run, timing_bin_key_indices)
-        #print(len(detector.shape))
         expected_length = len(run.time_bins)+1
         if len(detector.shape) < 2:
             reduced_array = np.zeros((expected_length))
@@ -592,7 +607,7 @@ def patch_pixels(self,run,detector_key,  mode='average', patch_range=4, deg=1, p
             self.patch_pixel(run,detector_key,pixel,mode,patch_range,deg,poly_range,axis=axis)
 
 
-    def patch_pixel(self, run, detector_key, pixel, mode='average', patch_range=4, deg=1, poly_range=6,axis=1):
+    def patch_pixel(self, run, detector_key, pixel, mode='average', patch_range=4, deg=1, poly_range=6, axis=1):
         """
         EPIX detector pixel patching.
         TODO: extend to patch regions instead of per pixel.
@@ -603,9 +618,9 @@ def patch_pixel(self, run, detector_key, pixel, mode='average', patch_range=4, d
         pixel : integer
             Pixel point to be patched
         mode : string
-            Determined which mode to use for patching the pixel. Averaging works well.
+            Determines which mode to use for patching the pixel. Averaging works well.
         patch_range : integer
-            pixels away from the pixel to be patched to be used for patching. Needed if multiple pixels in a row are an issue.
+            Pixels away from the pixel to be patched to be used for patching. Needed if multiple pixels in a row are an issue.
         deg : integer
             Degree of polynomial if polynomial patching is used.
         poly_range : integer
@@ -616,28 +631,62 @@ def patch_pixel(self, run, detector_key, pixel, mode='average', patch_range=4, d
             The original data with the new patch values.
         """
         data = getattr(run, detector_key)
-        if mode == 'average':
-            neighbor_values = data[:, pixel - patch_range:pixel + patch_range + 1, :]
-            new_val=np.sum(neighbor_values, axis=1) / neighbor_values.shape[1]
 
-            if axis==1:
-                data[:, pixel, :] = new_val
-            elif axis==2:
-                data[:,:,pixel]=new_val
-        elif mode == 'polynomial':
-            patch_x = np.arange(pixel - patch_range - poly_range, pixel + patch_range + poly_range + 1, 1)
+        def get_neighbor_values(data, pixel, patch_range, axis):
+            axis_slice = [slice(None)] * data.ndim
+            start_index = max(pixel - patch_range, 0)
+            end_index = min(pixel + patch_range + 1, data.shape[axis])
+            axis_slice[axis] = slice(start_index, end_index)
+            return data[tuple(axis_slice)]
+
+        def patch_value_average(data, pixel, patch_range, axis):
+            neighbor_values = get_neighbor_values(data, pixel, patch_range, axis)
+            neighbor_values = np.moveaxis(neighbor_values, axis, 0)
+            new_val = np.mean(neighbor_values, axis=0)
+            return new_val
+
+        def patch_value_polynomial(data, pixel, patch_range, poly_range, deg, axis):
+            patch_x = np.arange(pixel - patch_range - poly_range, pixel + patch_range + poly_range + 1)
             patch_range_weights = np.ones(len(patch_x))
-            patch_range_weights[pixel - patch_range - poly_range:pixel + patch_range + poly_range] = 0.001
-            coeffs = np.polyfit(patch_x, data[pixel - patch_range - poly_range:pixel + patch_range + poly_range + 1, :], deg,
-                                w=patch_range_weights)
-            data[pixel, :] = np.polyval(coeffs, pixel)
+            patch_range_weights[patch_range:-patch_range] = 0.001
+
+            neighbor_values = get_neighbor_values(data, pixel, patch_range + poly_range, axis)
+            neighbor_values = np.moveaxis(neighbor_values, axis, 0)
+
+            new_vals = []
+            for idx in range(neighbor_values.shape[1]): 
+                ys = neighbor_values[:, idx]
+                coeffs = np.polyfit(patch_x, ys, deg, w=patch_range_weights)
+                new_vals.append(np.polyval(coeffs, pixel))
+            return np.array(new_vals)
+
+        def patch_value_interpolate(data, pixel, patch_range, poly_range, axis):
+            patch_x = np.arange(pixel - patch_range - poly_range, pixel + patch_range + poly_range + 1)
+            neighbor_values = get_neighbor_values(data, pixel, patch_range + poly_range, axis)
+            neighbor_values = np.moveaxis(neighbor_values, axis, 0)
+
+            new_vals = []
+            for idx in range(neighbor_values.shape[1]):
+                ys = neighbor_values[:, idx]
+                interp_func = interp1d(patch_x, ys, kind='quadratic')
+                new_vals.append(interp_func(pixel))
+            return np.array(new_vals)
+
+        if mode == 'average':
+            new_val = patch_value_average(data, pixel, patch_range, axis)
+        elif mode == 'polynomial':
+            new_val = patch_value_polynomial(data, pixel, patch_range, poly_range, deg, axis)
         elif mode == 'interpolate':
-            patch_x = np.arange(pixel - patch_range - poly_range, pixel + patch_range + poly_range + 1, 1)
-            interp = interp1d(patch_x, data[pixel - patch_range - poly_range:pixel + patch_range + poly_range + 1, :],
-                              kind='quadratic')
-            data[pixel, :] = interp(pixel)
-        setattr(run,detector_key,data)
-        run.update_status('Detector %s pixel %d patched. Old value.'%(detector_key, pixel ))
+            new_val = patch_value_interpolate(data, pixel, patch_range, poly_range, axis)
+        else:
+            raise ValueError(f"Unsupported mode: {mode}")
+
+        patch_slice = [slice(None)] * data.ndim
+        patch_slice[axis] = pixel
+        data[tuple(patch_slice)] = new_val
+
+        setattr(run, detector_key, data)
+        run.update_status(f"Detector {detector_key} pixel {pixel} patched. Old value.")
 
     def patch_pixels_1d(self,run,detector_key,  mode='average', patch_range=4, deg=1, poly_range=6):
         """
@@ -700,7 +749,7 @@ def patch_pixel_1d(self, run, detector_key, pixel, mode='average', patch_range=4
                               kind='quadratic')
             data[pixel, :] = interp(pixel)
         setattr(run,detector_key,data)
-        run.update_status('Detector %s pixel %d patched.'%(detector_key, pixel ))
+        run.update_status('Detector %s pixel %d patched in mode %s'%(detector_key, pixel,mode ))
 
 
 
@@ -750,57 +799,7 @@ def make_energy_axis(self, run,energy_axis_length, A, R,  mm_per_pixel=0.05, d=0
 
         setattr(run,self.xes_line+'_energy',xaxis[::-1])
         run.update_status('XES energy axis generated for %s'%(self.xes_line))
-    def reduce_detector_spatial(self, run, detector_key, shot_range=[0, None], rois=[[0, None]], reduction_function=np.sum,  purge=True, combine=True,adu_cutoff=3.0):
-        """
-        Reduce spatial dimensions of detector data by combining or applying reduction functions over regions of interest (ROIs).
 
-        Parameters
-        ----------
-        run : object
-            The spectroscopy run instance.
-        detector_key : str
-            The key corresponding to the detector data.
-        shot_range : list of int, optional
-            The range of shots to consider for reduction (default is [0, None]).
-        rois : list of lists of int, optional
-            The regions of interest for spatial reduction (default is [[0, None]]).
-        reduction_function : function, optional
-            The function to use for reduction (default is np.sum).
-        purge : bool, optional
-            Whether to purge the original detector data after reduction (default is True).
-        combine : bool, optional
-            Whether to combine ROIs into a single reduced dataset (default is True).
-        adu_cutoff : float, optional
-            The ADU cutoff value for filtering (default is 3.0).
-        """
-        detector = getattr(run, detector_key)
-        if combine:
-
-            roi_combined = [rois[0][0], rois[-1][1]]  # Combined ROI spanning the first and last ROI
-            mask = np.zeros(detector.shape[2], dtype=bool)
-            for roi in rois:
-                mask[roi[0]:roi[1]] = True
-            masked_data = detector[shot_range[0]:shot_range[1], :, :][:, :, mask]
-            #masked_data = masked_data * (masked_data > adu_cutoff)
-            #Note the adu filtering should be handled at the controller level
-            reduced_data = reduction_function(masked_data, axis=2)
-            roi_indices = ', '.join([f"{roi[0]}-{roi[1]}" for roi in rois])
-            run.update_status(f"Spatially reduced detector: {detector_key} with combined ROI indices: {roi_indices}")
-            setattr(run, f"{detector_key}_ROI_1", reduced_data)
-        else:
-            for idx, roi in enumerate(rois):
-                data_chunk = detector[shot_range[0]:shot_range[1], roi[0]:roi[1]]
-                reduced_data = reduction_function(data_chunk, **kwargs)
-            if roi[1] is None:
-                roi[1] = detector.shape[1] - 1
-                run.update_status(f"Spatially reduced detector: {detector_key} with ROI: {roi[0]}, {roi[1]}")
-                setattr(run, f"{detector_key}_ROI_{idx+1}", reduced_data)
-        if purge:
-            #pass
-            setattr(run, detector_key,None)
-            #delattr(run, detector_key)
-            #del run.detector_key
-            run.update_status(f"Purged key to save room: {detector_key}")
     def reduce_det_scanvar(self, run, detector_key, scanvar_key, scanvar_bins_key):
         """
         Reduce detector data by binning according to an arbitrary scan variable.