Added examples

Reviewed-on: #1

README.md

@@ -4,34 +4,64 @@ This pacakge is used to generate region masks separating any focused X-ray flux
 
 ## Installation
 This package is to be used with Python 3.x.x
-```python
-pip install git+https://github.com/andrey-rrousan/nuwavdet
+```bash
+pip install git+http://heagit.cosmos.ru/nustar/nuwavdet.git
 ```
 
-## Main use
+To update the package to the current version one should delete the previous version
+```bash
+pip uninstall nuwavdet
+```
 
-To use the package in your project, import it in by writing
+And simply repeat the intallation procedure again from the repository.
+
+## Installation verification
+If the installation was successful the package can be used with the following import:
 
 ```python
 from nuwavdet import nuwavdet as nw
 ```
 
-The main functionality of the pacakge is presented with a single function
+To verify the installation we suggest running a simple script:
+
 ```python
-process(obs_path, thresh)
+from nuwavdet import nuwavdet as nw
+
+print(nw.binary_array(2))
 ```
+
+The output of the script should be
+
+```bash
+[[False False]
+ [False  True]
+ [ True False]
+ [ True  True]]
+```
+
+## Main use
+
+The main functionality of the package is presented with a single function
+```python
+nw.process(obs_path, thresh)
+```
+
 Inputs are string with path to the _cl.evt file to use and a tuple of thresholds, e.g.
 ```python
-process('D:\\Data\\obs_cl.evt', (3, 2))
+nw.process('D:\\Data\\obs_cl.evt', (3, 2))
 ```
 
-Outputs of the function are:
-1. dictionary with some metadata and properties of the observation after mask generation procedure.
-2. region array with mask in DET1 coordinate frame. Note that this mask is for numpy mask application so 1 corresponds to masked pixel and 0 otherwise.
-3. custom bad pixel table with flagged pixels in RAW coordinates. It can be exported as fits file or each separate table can be acessed directly.
-4. array with the sum of wavelet planes used in the processing.
+The detailed script description of the data extraction with the script is presented in the examples folder of the repository.
 
-Metadata about the observation file:
+The function nw.process returns severl python objects:
+1. python-dictionary with some metadata and properties of the observation after mask generation procedure.
+2. region array with mask in DET1 coordinate frame. Note that this mask is for numpy mask application so  True (1) corresponds to masked pixel and False (0) otherwise.
+3. custom bad pixel table with flagged pixels in RAW coordinates. It can be exported as fits file for further application to the nupipeline as fpma_userbpfile or fpmb_userbpfile. 
+4. array with the sum of wavelet planes for potential alternative applications.
+
+Metadata about the observation returned by the nw.process is:
+
+Observation metadata: 
 
 1. OBS_ID
 2. Detector
@@ -47,14 +77,8 @@ Useful algorythm-related data:
 
 ## Other uses
 
-You can process the cl.evt file by creating an Observation class object:
+Other possbile usecases are shown in the examples folder.
 
-```python
-obs = nw.Observation(path_to_evt_file)
-```
+## Contact information
 
-Additionally, the energy band in KeV to get events from can be passed as an argument. The default value is [3,20].
-
-```python
-obs = nuwavsource.Observation(path_to_evt_file,E_borders=[E_min,E_max])
-```
+If you have any questions or issues with the code, feel free to contact Andrey Mukhin: amukhin@cosmos.ru

examples/1_save_results.py

@@ -0,0 +1,54 @@
+from nuwavdet import nuwavdet as nw
+
+
+OBS_PATH = r'.//path_to_obs//nu<obsid><DET>01_cl.evt'
+THRESH = (3, 2)
+
+SAVE_BADPIX_PATH = r'.//out//badpix.fits'
+SAVE_REGION_PATH = r'.//out//region.fits'
+SAVE_WAVSUM_PATH = r'.//out//wavsum.fits'
+
+METADATA_PATH = r'.//out//metadata.csv'
+METADATA_FITS_PATH = r'.//out//metadata.fits'
+
+
+if __name__ == '__main__':
+    # PROCESS THE OBSERVATION WITH GIVEN THRESHOLD
+    result, region, region_raw, wav_sum = nw.process(OBS_PATH, thresh=THRESH)
+
+    # SAVE THE REGION BAD PIXEL FILES TO THE FITS FILE WITH NUPIPELINE
+    # COMPATIBLE FORMAT AND HEADERS.
+    region_raw.writeto(SAVE_BADPIX_PATH)
+
+    # SAVE REGION MASK AS A FITS IMAGE
+    nw.save_region(region, SAVE_REGION_PATH, overwrite=False)
+    # Note that the Python script uses numpy masked array with
+    # True (1) as as masked and False (0) as unmasked pixel.
+    # nw.save_region transfers the numpy masked array to
+    # conventional format with 1 for unmasked and 0 for masked pixel.
+    # However, if mask is used in the Python you need to transfer it back with
+    # numpy.logical_not(mask).
+
+    # SAVE WAVSUM ARRAY AS A FITS IMAGE
+    nw.fits.writeto(SAVE_WAVSUM_PATH, wav_sum, overwrite=False)
+
+    # SAVE METADATA
+    # WE SUGGEST SAVING ALL THE METADATA FOR SEVERAL OBSERVATIONS
+    # IN ONE FILE.
+
+    # CREATE CSV FILE TO SAVE DATA
+    # IF FILE ALREADY EXISTS YOU SHOULD REMOVE THIS BLOCK FROM YOUR CODE
+    table = {
+        'obs_id': [], 'detector': [], 'ra': [], 'dec': [],
+        'lon': [], 'lat': [], 't_start': [], 'exposure': [],
+        'count_rate': [], 'remaining_area': [], 'cash_stat': [],
+        'cash_stat_full': []
+    }
+    out_table = nw.DataFrame(table)
+    out_table.to_csv(METADATA_PATH)
+
+    # SAVE DATA TO CREATED CSV
+    nw.DataFrame(result, index=[0]).to_csv(METADATA_PATH, mode='a', header=False)
+
+    # TRANSFORM THE CSV TO FITS-TABLE
+    nw.csv_to_table(METADATA_PATH, METADATA_FITS_PATH)

examples/2_directory_processing.py

@@ -0,0 +1,33 @@
+from nuwavdet import nuwavdet as nw
+
+INPUT_FOLDER = r'path_to_directory'
+OUTPUT_FOLDER = r'.//Output'
+
+
+if __name__ == '__main__':
+    # BEGIN PROCESSING ALL THE OBSERVATIONS INSIDE THE FOLDER
+    nw.process_folder(input_folder=INPUT_FOLDER,
+                      start_new_file='y',
+                      fits_folder=OUTPUT_FOLDER,
+                      thresh=(3, 2),
+                      cpu_num=10
+                      )
+
+    # IF THE PROCESSING WAS INTERRUPTED YOU CAN CONTINUE IT WITH THE SAME CODE
+    # BY CHANGING THE start_new_file TO 'n'.
+
+    # THE STRUCTURE OF THE FOLDER IS
+    # OUTPUT_FOLDER
+    # __overview.csv                csv-table with observations metadata
+    # __overvies.fits               fits-table with the same metadata
+    # __overview_skipped.csv        csv-table with the skipped observations
+    # __Region                      folder for region mask images
+    # ____<obsid><DET>_region.fits
+    # __Region_raw                  folder for region masks in RAW coordinates
+    # ____<obsid><DET>_reg_raw.fits
+    # __Wav_sum                     folder for sum of wavelet layers
+    # ____<obsid><DET>_wav_sum.fits
+
+    # Note nw.process_folder uses multiprocessing with cpu_num cores.
+    # The number of cores can be manually chosen or automatically
+    # detected if cpu_num = 0. 

examples/3_wavelet.py

@@ -0,0 +1,23 @@
+from nuwavdet import nuwavdet as nw
+
+OBS_PATH = r'.//path_to_obs//nu<obsid><DET>01_cl.evt'
+THRESH = (3, 2)
+
+
+if __name__ == '__main__':
+    # CREATE THE OBSERVATION CLASS OBJECT
+    obs = nw.Observation(OBS_PATH)
+
+    # CALCULATE THE WAVLET LAYERS WITH GIVEN THRESHOLD
+    wav_layers = obs.wavdecomp(mode='atrous', occ_coeff=True, thresh=THRESH)
+
+    # ALL THE LAYERS CAN BE ACCESSED AS AN ELEMENT OF wav_layers VARIABLE
+    # wav_layers[0] for the 1st wavelet layer
+    # wav_layers[4] for 5th wavelet layer
+    # wav_layers[-1] for the last wavelet layer
+    # wav_layers[2:5] for the list of the layers from 3 to 5
+    # wav_layers[[1, 3, 5]] for the list of layers 2, 4 and 6
+
+    # To calculate the sum of wavelet layers one should use sum() method
+    # wav_layers[2:7].sum(0) returns a sum of layers from 3 to 7
+    # wav_layers[[1, 3, 5]].sum(0) returns a sum of layers 2, 4 and 6.

examples/4_cstat.py

@@ -0,0 +1,23 @@
+from nuwavdet import nuwavdet as nw
+import numpy as np
+
+
+OBS_PATH = r'.//path_to_obs//nu<obsid><DET>01_cl.evt'
+MASK_PATH = r'.//path_to_mask//<obsid><DET>.fits'
+
+
+if __name__ == '__main__':
+    # CREATE THE OBSERVATION CLASS OBJECT
+    obs = nw.Observation(OBS_PATH)
+
+    # READ THE REGION MASK FILE
+    region = nw.fits.getdata(MASK_PATH)
+
+    # TRANSFORM REGION MASK DATA TO NUMPY MASK DATA (SEE 1_save_results.py).
+    region = np.logical_not(region.astype(bool))
+
+    # CREATE MASKED ARRAY CLASS OBJECT
+    masked_data = np.ma.masked_array(obs, mask=region)
+
+    # CALCULATE THE CSTAT ON THE MASKED DATA
+    print(nw.сstat(masked_data.mean(), masked_data))

nuwavdet/nuwavdet.py

@@ -30,6 +30,14 @@ def get_link_list(folder: str, sort_list: bool = True) -> list[str]:
         return np.array(links)
 
 
+def csv_to_table(csv_path, fits_path):
+    """
+    Transform the csv table to fits table with astropy.
+    """
+    csv_file = read_csv(csv_path, index_col=0, dtype={'obs_id': str})
+    Table.from_pandas(csv_file).write(fits_path, overwrite=True)
+
+
 def binary_array(num: int) -> list[list[bool]]:
     """
     Returns list of all possible combinations of num of bool values.
@@ -199,6 +207,7 @@ def count_binning(array, count_per_bin: int = 2):
 
 def cstat(expected, data: list, count_per_bin: int = 2) -> float:
     _data = data.flatten()
+    if type(data) is np.ma.masked_array:
         _data = _data[_data.mask == False]
     _expected = expected
     c_stat = 0
@@ -235,7 +244,7 @@ class Observation:
         resized_coeff = (coeff).reshape(2, 2).repeat(180, 0).repeat(180, 1)
         return resized_coeff
 
-    def get_data(self, file, E_borders=[3, 20]):
+    def get_data(self, file, E_borders=[3, 20], generate_mask=True):
         """
         Returns masked array with DET1 image data for given energy band.
         Mask is created from observations badpix tables and to mask the border and gaps. 
@@ -248,10 +257,12 @@ class Observation:
         data_mask = data[np.logical_not(idx_mask)]
         build_hist = lambda array: np.histogram2d(array['DET1Y'], array['DET1X'], 360, [[0, 360], [0, 360]])[0]
         output = build_hist(data_output)
+        if generate_mask:
             mask = build_hist(data_mask)
             mask = np.logical_or(mask, add_borders(output))
             mask = np.logical_or(mask, self.get_bad_pix(file))
             return output, mask
+        return output
 
     def get_bad_pix(self, file, threshold=0.9):
         """
@@ -276,7 +287,7 @@ class Observation:
         correction_poiss = np.random.poisson(corr*array, corr.shape)
         return array + correction_poiss
 
-    def wavdecomp(self, mode='gauss', thresh=False, occ_coeff=False):
+    def wavdecomp(self, mode='gauss', thresh=0, occ_coeff=False):
         """
         Performs a wavelet decomposition of image.
         """
@@ -330,7 +341,7 @@ class Observation:
         """
         Returns a hdu_list with positions of masked pixels in RAW coordinates.
         """
-        x_region, y_region = np.where(region)
+        y_region, x_region = np.where(region)
         hdus = []
         for i in range(4):
             current_dir = os.path.dirname(__file__)
@@ -390,7 +401,6 @@ def process(obs_path, thresh):
     Arguments: path to the file of interest and threshold,
     e.g. process('D:\\Data\\obs_cl.evt', (3, 2))
     """
-    bin_num = 6
 
     table = {
         'obs_id': [], 'detector': [], 'ra': [], 'dec': [],
@@ -405,18 +415,19 @@ def process(obs_path, thresh):
                              dec=obs.dec*u.deg,
                              frame='fk5').transform_to('galactic')
         lon, lat = sky_coord.l.value, sky_coord.b.value
-        rem_signal, rem_area, poiss_comp, rms = np.zeros((4, 2**bin_num))
+        useful_bin_num = 6
+        rem_signal, rem_area, poiss_comp, rms = np.zeros((4, 2**useful_bin_num))
         region = np.zeros(obs.data.shape, dtype=bool)
         region_raw = -1
         rem_region = np.logical_and(region, np.logical_not(obs.data.mask))
         masked_obs = np.ma.masked_array(obs.data, mask=region)
-        good_lvl = np.zeros(bin_num, dtype=bool)
-        good_idx = 0
+        good_lvl = np.zeros(useful_bin_num, dtype=bool)
         if obs.exposure > 1000:
             wav_obs = obs.wavdecomp('atrous', thresh, occ_coeff=True)
             wav_sum = wav_obs[2:-1].sum(0)
             occ_coeff = obs.get_coeff()
-            binary_arr = binary_array(bin_num)
+            binary_arr = binary_array(useful_bin_num)
+            good_idx = len(binary_arr) - 1
 
             for idx, lvl in enumerate(binary_arr):
                 try:
@@ -434,12 +445,9 @@ def process(obs_path, thresh):
                 rms[idx] = np.sqrt(((masked_obs-masked_obs.mean())**2).mean())
 
             for idx in range(len(poiss_comp)):
-                if ((poiss_comp[idx] < poiss_comp[good_idx]) and
-                    (poiss_comp[idx] < poiss_comp[-1] + 0.05) and
-                    (rem_area[idx] > rem_area[-1])):
+                if ((poiss_comp[idx] < poiss_comp[-1] + 0.05) and
+                    (rem_area[idx] > rem_area[good_idx])):
                     good_idx = idx
-            if good_idx == 0:
-                good_idx = len(binary_arr) - 1
             good_lvl = binary_arr[good_idx]
 
             try:
@@ -483,7 +491,7 @@ def process(obs_path, thresh):
                         ]
 
         for key, value in zip(table.keys(), to_table):
-            table[key] = [value]
+            table[key] = value
         return table, region.astype(int), region_raw, wav_sum
     except TypeError:
         return obs_path, -1, -1, -1
@@ -494,7 +502,7 @@ def _process_multi(args):
 
 
 def process_folder(input_folder=None, start_new_file=None, fits_folder=None,
-                   thresh=None):
+                   thresh=None, cpu_num=0):
     """
     Generates a fits-table of parameters, folder with mask images in DET1 and
     BADPIX tables in RAW for all observations in given folder.
@@ -586,6 +594,14 @@ def process_folder(input_folder=None, start_new_file=None, fits_folder=None,
     # START PROCESSING
     print('Started processing...')
     num = 0
+    if cpu_num == 0:
+        cpu_num = cpu_count()
+    elif cpu_num < 0:
+        raise ValueError('cpu_num must be a positive integer')
+    elif cpu_num > cpu_count():
+        print('Chosen cpu_num exceed the number of CPU cores. Using cpu_count() instead.')
+        cpu_num = cpu_count()
+
     for group_idx in range(len(obs_list)//group_size+1):
         print(f'Started group {group_idx}')
         group_list = obs_list[group_size*group_idx:min(group_size*(group_idx+1), len(obs_list))]
@@ -593,7 +609,7 @@ def process_folder(input_folder=None, start_new_file=None, fits_folder=None,
             os.stat(file).st_size/2**20
             for file in group_list
         ]).max()
-        process_num = (cpu_count() if max_size < 50 else (cpu_count()//2 if max_size < 200 else (cpu_count()//4 if max_size < 1000 else 1)))
+        process_num = (cpu_num if max_size < 50 else (cpu_num//2 if max_size < 200 else cpu_num//4 if max_size < 1000 else cpu_num//8))
         print(f"Max file size in group is {max_size:.2f}Mb, create {process_num} processes")
         with get_context('spawn').Pool(processes=process_num) as pool:
             packed_args = map(lambda _: (_, thresh), group_list)
@@ -604,14 +620,14 @@ def process_folder(input_folder=None, start_new_file=None, fits_folder=None,
                     num += 1
                     continue
 
-                obs_name = str(result['obs_id'][0])+result['detector'][0]
-                if result['exposure'][0] < 1000:
+                obs_name = str(result['obs_id'])+result['detector']
+                if result['exposure'] < 1000:
                     print(f'{num:>3} {obs_name} is skipped. Exposure < 1000')
-                    DataFrame(result).to_csv(os.path.join(fits_folder, 'overview_skipped.csv'), mode='a', header=False)
+                    DataFrame(result, index=[0]).to_csv(os.path.join(fits_folder, 'overview_skipped.csv'), mode='a', header=False)
                     num += 1
                     continue
 
-                DataFrame(result).to_csv(os.path.join(fits_folder, 'overview.csv'), mode='a', header=False)
+                DataFrame(result, index=[0]).to_csv(os.path.join(fits_folder, 'overview.csv'), mode='a', header=False)
                 save_region(region, os.path.join(region_folder, f'{obs_name}_region.fits'), overwrite=True)
                 region_raw.writeto(os.path.join(region_raw_folder, f'{obs_name}_reg_raw.fits'), overwrite=True)
                 fits.writeto(os.path.join(wav_sum_folder, f'{obs_name}_wav_sum.fits'), wav_sum, overwrite=True)

setup.py

@@ -5,10 +5,10 @@ with open("README.md", "r") as fh:
 
 setuptools.setup(
     name="nuwavdet",
-    version="0.1.0",
+    version="0.1.1",
     author="Andrey Mukhin",
-    author_email="amukhin@phystech.edu",
-    description="A package for source exclusion in NuStar observation data using wavelet decomposition",
+    author_email="amukhin@cosmos.ru",
+    description="A package for source exclusion in NuSTAR observation data using wavelet decomposition",
     long_description=long_description,
     long_description_content_type="text/markdown",
     url="https://github.com/andrey-rrousan/nuwavdet",