add integer compression error estimation script

stixdcpy/integer_compression.py

@@ -13,17 +13,35 @@
 """
 
 import numpy as np
+from pprint import pprint
 
 MAX_STORED_INTEGER = 1e8
 
 
 
 class Compression(object):
-    def __init__(self, k,m, with_stat_error=True):
-        self.lut=Compression.get_error_lut(k,m, with_stat_error)
+    def __init__(self,s, k,m, include_stat_error=True):
+        self.skm=(s,k,m)
+        self.lut=Compression.get_error_lut(s,k,m, include_stat_error)
+
+    def get_errors(self,counts:np.array):
+        """
+            calculate errors for counts
+        """
+        fv=np.vectorize(self.get_error, otypes=[np.float])
+        #otypes must be specified
+        return fv(counts)
+
+
 
     def get_error(self,counts):
-        return self.lut[int(counts)]
+        if counts==0:
+            return 0
+        try:
+            return self.lut[counts]
+        except KeyError:
+            s,k,m=self.skm
+            raise Exception(f'Failed to error of {counts}! Could the compression scheme ({s=},{k=},{m=}) wrong?')
 
     @staticmethod
     def decompress(x,S, K, M):
@@ -54,7 +72,7 @@ def decompress(x,S, K, M):
 
         x0 = 1 << (M + 1)
         if x < x0:
-            return None, x
+            return  x,0
         mask1 = (1 << M) - 1
         mask2 = (1 << M)
         mantissa1 = x & mask1
@@ -68,12 +86,12 @@ def decompress(x,S, K, M):
         #error of a flat distribution
 
         if mean > MAX_STORED_INTEGER:
-            return error, float(mean)
+            return float(mean),error
 
-        return error, sign * mean
+        return  sign * mean, error
 
     @staticmethod
-    def get_error_lut(k,  m, with_stat_error=True, s=0):
+    def get_error_lut(s, k,  m, include_stat_error=True):
         """
         Creates a lookup table for error calculation based on k and m values.
 
@@ -86,14 +104,15 @@ def get_error_lut(k,  m, with_stat_error=True, s=0):
         """
         res = {}
         for i in range(256):
-            err, val = Compression.decompress(i,s, k, m)
+            val,err = Compression.decompress(i,s, k, m)
             if err is not None:
-                res[val] = np.sqrt(err ** 2 + val) if with_stat_error else err
+                res[val] = np.sqrt(err ** 2 + np.abs(val)) if include_stat_error else err
         return res
 
-    #error_luts = {'035': make_lut(3, 5), '044': make_lut(4, 4), '053': make_lut(5, 3)}
 
-#if __name__ == '__main__':
-#    from pprint import pprint
-#    pprint(error_luts)
+if __name__ == '__main__':
+    from pprint import pprint
+    c=Compression(0,4,4)
+    error_luts=c.get_error_lut(0,4,4)
+    pprint(error_luts)
 

stixdcpy/science.py

@@ -13,8 +13,11 @@
 from stixdcpy.logger import logger
 from stixdcpy import io as sio
 from stixdcpy.net import FitsQuery as freq
+from stixdcpy import integer_compression as sic
+
 from stixdcpy import net as net
 from stixdcpy import instrument as inst
+
 from pathlib import PurePath
 from datetime import datetime as dt
 from datetime import timedelta as td
@@ -44,7 +47,6 @@ def __init__(self, fname=None, request_id=None):
         self.energies = []
         self.corrected = None
         # self.read_data()
-        self.skm = {}
 
     @property
     def url(self):
@@ -69,13 +71,7 @@ def filename(self):
         return self.fname
 
 
-    #def init_skm(self):
-    #    try:
-    #        self.skm={'counts':self.hdul['CONTROL'].data['compression_scheme_counts_skm'],
-    #                'triggers':self.hdul['CONTROL'].data['compression_scheme_triggers_skm']
-    #                }
-    #    except (KeyError,TypeError, ValueError):
-    #        pass
+
 
 
 
@@ -108,10 +104,9 @@ def read_fits(self, light_time_correction=True):
         self.timedel = self.data['timedel']
         self.time = self.data['time']
 
-        if self.is_time_bin_shifted(self.T0_unix) and len(self.timedel) > 1:
+        if self.is_time_bin_shifted() and len(self.timedel) > 1:
             self.timedel = self.timedel[:-1]
             self.time = self.time[1:]
-            logger.info('Shifted time bins have been corrected automatically!')
             if self.data_type == 'PixelData':
                 self.counts = self.counts[1:, :, :, :]
                 self.triggers = self.triggers[1:, :]
@@ -120,6 +115,7 @@ def read_fits(self, light_time_correction=True):
                 self.counts = self.counts[1:, :]
                 self.triggers = self.triggers[1:]
                 # self.rcr = self.rcr[1:]
+            logger.info('Shifted time bins have been corrected automatically!')
 
         self.request_id = self.hdul['CONTROL'].data['request_id']
 
@@ -167,7 +163,7 @@ def read_fits(self, light_time_correction=True):
             self.count_rates = self.counts / self.timedel[:, None]
             self.trigger_rates = self.triggers / self.timedel
 
-    def is_time_bin_shifted(self, unix_time):
+    def is_time_bin_shifted(self):
         """
             Time bins are shifted in the data collected before 2021-12-09 due a bug in the flight software
 
@@ -178,8 +174,7 @@ def is_time_bin_shifted(self, unix_time):
                 is_shifted: bool
                     True if time bin is shifted else False
         """
-
-        return (unix_time < sdt.utc2unix('2021-12-09T14:00:00'))
+        return (self.T0_unix< sdt.utc2unix('2021-12-09T14:00:00'))
 
     @classmethod
     def from_sdc(cls, request_id, level='L1'):
@@ -274,6 +269,7 @@ def __init__(self, fname, request_id, ltc=False):
         self.pixel_count_rates = None
         self.correct_pixel_count_rates = None
         self.read_fits(light_time_correction=ltc)
+        self.pixel_counts_comp_stat_err= None
         self.make_spectra()
 
 
@@ -294,15 +290,46 @@ def make_spectra(self, pixel_counts=None):
 
         self.pixel_total_counts = np.sum(self.pixel_counts, axis=(0, 3))
 
+    def compute_errors_from_counts(self, skm=None):
+        """
+            compute errors using counts
+        Args
+            include_stat_error: include statistical errors in the final results
+
+        """
+
+        if skm is None:
+            try:
+                skm=self.hdul['CONTROL'].data['compression_scheme_counts_skm']
+                s,k,m=skm.flatten()
+            except (KeyError,TypeError, ValueError):
+                raise Exception("Couldn't find SKM found in the FITS file!")
+        else:
+            s,k,m=skm
+        comp = sic.Compression(s,k,m, include_stat_error=True)
+        self.pixel_counts_comp_stat_err=comp.get_errors(self.counts)
+
+
     @property
     def pixel_counts_error(self):
-        try:
-            counts_err=self.hdul['data'].data['counts_err']
-        except KeyError:
-            counts_err=self.hdul['data'].data['counts_comp_err']
-        return error_computation(counts_err,
+        """
+            compute or read counts errors from FITS file
+        """
+
+        if self.pixel_counts_comp_stat_err is not None:
+            #if already computed
+            return self.pixel_counts_comp_stat_err
+        else:
+            try:
+                counts_err=self.hdul['data'].data['counts_err']
+            except KeyError:
+                counts_err=self.hdul['data'].data['counts_comp_err']
+                #in some FITS files compression errors are zeros
+            return error_computation(counts_err,
                                  self.pixel_counts)
 
+
+
 
 
     def correct_dead_time(self):
@@ -351,7 +378,6 @@ def correct(triggers, counts_arr, counts_err_arr, time_bins):
 
             count_rate = counts_arr / time_bins
             count_rate_err = counts_err_arr / time_bins
-            # print(counts_arr.shape)
             for det in range(32):
                 trig_idx = inst.detector_id_to_trigger_index(det)
                 nin = photons_in[:, trig_idx]
@@ -362,13 +388,16 @@ def correct(triggers, counts_arr, counts_err_arr, time_bins):
 
             #live_ratio=np.zeros_like(live_ratio)+1
             #disable live time correction
+            live_ratio=live_ratio[:, :, None, None]
 
-            corrected_rate = count_rate / live_ratio[:, :, None, None]
+            corrected_rate = count_rate / live_ratio
 
-            corrected_rate_err = count_rate_err / live_ratio[:, :, None, None]
+            corrected_rate_err =count_rate_err/ live_ratio
 
             corrected_counts = corrected_rate * time_bins
             corrected_counts_err = corrected_rate_err * time_bins
+            #errors of live ratio not taken into account yet
+
 
             return {
                 'corrected_rates': corrected_rate,
@@ -465,8 +494,7 @@ def get_sum_counts(self, start_utc=None, end_utc=None) :
                 + sum_counts['bottom_err']**2 )
 
         sum_counts['total'] = sum_counts['top']+sum_counts['bottom']+sum_counts['small']
-        sum_counts['total_err'] = np.sqrt(sum_counts['top_err']**2 
-                + sum_counts['bottom_err']**2 + sum_counts['small_err']**2)
+        sum_counts['total_err'] = np.sqrt(sum_counts['big_err']**2  + sum_counts['small_err']**2)
 
         return sum_counts
 
@@ -617,7 +645,6 @@ def __init__(self, l1sig: PixelData, l1bkg: PixelData):
         ])
         # set counts beyond the signal energy range to 0
         self.subtracted_counts = (self.l1sig.counts - self.pixel_bkg_counts)
-        #print(self.l1sig.inversed_energy_bin_mask)
         self.subtracted_counts *= self.l1sig.energy_bin_mask
 
         # Dead time correction needs to be included in the future
@@ -789,7 +816,6 @@ def correct_dead_time(self) -> dict:
         def correct(triggers, counts_arr, time_bins, num_detectors):
             ''' correct dead time using triggers '''
             tau_conv_const = 1e-6
-
             photons_in = triggers / (time_bins * num_detectors -
                                      TRIG_TAU * triggers)
             # photon rate approximated using triggers
@@ -880,7 +906,8 @@ def fits_time_corrections(primary_header, tstart, tend):
     primary_header.set('DATE_EAR', date_ear)
     primary_header.set('DATE_SUN', date_sun)
     primary_header.set('MJDREF', Time(tstart).mjd)
-    # Note that MJDREF in these files is NOT 1979-01-01 but rather the observation start time! Affects time calculation later, but further processing of the file corrects this
+    # Note that MJDREF in these files is NOT 1979-01-01 but rather the 
+    # observation start time! Affects time calculation later, but further processing of the file corrects this
 
     return primary_header
 
@@ -982,7 +1009,8 @@ def time_select_indices(tstart, tend, primary_header, data_table, factor=1.):
 
 
 def spec_fits_crop(fitsfile, tstart, tend, outfilename=None):
-    """ Crop a STIX science data product (L1A, L1, or L4) to contain only the data within a given time interval. Create a new FITS file containing this data. The new file will be of the same processing level as the input file.
+    """ Crop a STIX science data product (L1A, L1, or L4) to contain only the data within a given time interval. 
+        Create a new FITS file containing this data. The new file will be of the same processing level as the input file.
 
     Inputs:
     fitsfile : str