Merge pull request #263 from texadactyl/master

Improve CPU utilisation in turbo_seti use of blank_dc function

Author: texadactyl

VERSION-HISTORY.md

@@ -4,8 +4,9 @@ This file is a version history of turbo_seti amendments, beginning with version
 
 | `YYYY_MM_DD` | `Version` | `Contents` |
 | :--: | :--: | :-- |
+| 2021-07-18 | 2.1.7 | Create a turbo_seti clone of blank_dc that is optional and uses a different strategy (issue #262). |
 | 2021-07-15 | 2.1.6 | Calculate normalized value inside hitsearch kernel on GPU-mode. |
-| 2021-07-16 | 2.1.5 | Failed to pass the gpu_id from find_doppler.py to data_handler.py. (issue #254). |
+| 2021-07-16 | 2.1.5 | Failed to pass the gpu_id from find_doppler.py to data_handler.py (issue #254). |
 | 2021-07-15 | 2.1.4 | Add GPU device selection with cli argument gpu_id. (issue #254). |
 | 2021-07-15 | 2.1.3 | Diagnose out of range time steps with correct messages (issue #256). |
 | | | Also, stop catching exceptions in seti_event.py which causes a cascade in tracebacks. |

setup.py

@@ -9,7 +9,7 @@
 
 from setuptools import setup, find_packages
 
-__version__ = "2.1.6"
+__version__ = "2.1.7"
 
 with open("turbo_seti/find_doppler/turbo_seti_version.py", "w") as fh:
     fh.write("TURBO_SETI_VERSION = '{}'\n".format(__version__))

test/test_turbo_seti.py

@@ -218,6 +218,14 @@ def test_turboSETI_entry_point():
     h5_5 = os.path.join(TESTDIR, OFFNIL_H5)
     args = [h5_5, "-P", "y", "-s", str(MIN_SNR), "-M", str(MAX_DRIFT), "-o", TESTDIR, ]
     seti_event.main(args)
+    print("\n===== test_turboSETI_entry_point 6 =====")
+    h5_5 = os.path.join(TESTDIR, VOYAH5)
+    args = [h5_5, "--blank_dc", "y", "-s", str(MIN_SNR), "-M", str(MAX_DRIFT), "-o", TESTDIR, ]
+    seti_event.main(args)
+    print("\n===== test_turboSETI_entry_point 7 =====")
+    h5_5 = os.path.join(TESTDIR, VOYAH5)
+    args = [h5_5, "--blank_dc", "n", "-s", str(MIN_SNR), "-M", str(MAX_DRIFT), "-o", TESTDIR, ]
+    seti_event.main(args)
 
 
 def test_make_waterfall_plots():
@@ -350,3 +358,4 @@ def test_flipx_kernel(kernels):
 if __name__ == "__main__":
     print("Please run: pytest test_turbo_seti.py")
     test_find_doppler_voyager()
+    test_turboSETI_entry_point()

turbo_seti/find_doppler/data_handler.py

@@ -11,6 +11,7 @@
 from blimpy import Waterfall
 from blimpy.io import sigproc
 
+from .helper_functions import cut_the_mid_spike
 from .kernels import Kernels
 
 logger = logging.getLogger('data_handler')
@@ -249,7 +250,7 @@ def __init__(self, filename, f_start=None, f_stop=None, t_start=None, t_stop=Non
         self.shoulder_size = 0
         self.tdwidth = self.fftlen + self.shoulder_size * self.tsteps
 
-    def load_data(self):
+    def load_data(self, flag_blank_dc=True):
         r"""
         Read the spectra and drift indices from file.
 
@@ -261,11 +262,16 @@ def load_data(self):
         self.fil_file.read_data(f_start=self.f_start, f_stop=self.f_stop)
 
         #Blanking DC bin.
-        if self.n_coarse_chan is not None:
-            n_coarse_chan = self.n_coarse_chan
+        if flag_blank_dc:
+            logger.debug("blank_dc is enabled.")
+            if self.n_coarse_chan is not None:
+                n_coarse_chan = self.n_coarse_chan
+            else:
+                n_coarse_chan = int(self.fil_file.calc_n_coarse_chan())
+            cut_the_mid_spike(self.fil_file.data, n_coarse_chan)
         else:
-            n_coarse_chan = int(self.fil_file.calc_n_coarse_chan())
-        self.fil_file.blank_dc(n_coarse_chan)
+            logger.debug("blank_dc is disabled.")
+
 
         dim_time = self.fil_file.data.shape[0]
         if dim_time < 2:

turbo_seti/find_doppler/find_doppler.py

@@ -73,11 +73,13 @@ class FindDoppler:
         Append output DAT & LOG files? (True/False)
     log_level_int : int, optional
         Python logging threshold level (INFO, DEBUG, or WARNING)
+    blank_dc : bool, optional
+        Use blank_dc() for spike smoothing.
 
     """
     def __init__(self, datafile, max_drift=4.0, min_drift=0.00001, snr=25.0, out_dir='./', coarse_chans=None,
                  obs_info=None, flagging=False, n_coarse_chan=None, kernels=None, gpu_backend=False, gpu_id=0,
-                 precision=2, append_output=False, log_level_int=logging.INFO):
+                 precision=2, append_output=False, log_level_int=logging.INFO, blank_dc=True):
 
         print(version_announcements)
 
@@ -113,10 +115,11 @@ def __init__(self, datafile, max_drift=4.0, min_drift=0.00001, snr=25.0, out_dir
         self.status = True
         self.flagging = flagging
         self.append_output = append_output
+        self.flag_blank_dc = blank_dc
         self.parms = 'datafile={}, max_drift={}, min_drift={}, snr={}, out_dir={}, coarse_chans={}' \
                         .format(datafile, max_drift, min_drift, snr, out_dir, coarse_chans) \
-                    + ', flagging={}, n_coarse_chan={}, kernels={}, gpu_id={}, gpu_backend={}' \
-                        .format(flagging, self.n_coarse_chan, kernels, gpu_id, gpu_backend) \
+                    + ', flagging={}, n_coarse_chan={}, kernels={}, gpu_id={}, gpu_backend={}, blank_dc={}' \
+                        .format(flagging, self.n_coarse_chan, kernels, gpu_id, gpu_backend, blank_dc) \
                     + ', precision={}, append_output={}, log_level_int={}, obs_info={}' \
                         .format(precision, append_output, log_level_int, obs_info)
         if min_drift < 0 or max_drift < 0:
@@ -166,7 +169,6 @@ def search(self, n_partitions=1, progress_bar='n'):
 
         # As of 2.1.0, add max_drift_rate and obs_length to FileWriter header input
         header_in['max_drift_rate'] = self.max_drift
-        #header_in['obs_length'] was already set in data_handler.py DATAH __init__
 
         wfilename = filename_in.split('/')[-1].replace('.h5', '').replace('.fits', '').replace('.fil', '')
         path_log = '{}/{}.log'.format(self.out_dir.rstrip('/'), wfilename)
@@ -198,7 +200,7 @@ def search(self, n_partitions=1, progress_bar='n'):
 
         # Run serial version
         if n_partitions == 1:
-            sched = Scheduler(load_the_data, [ (l, self.kernels.precision) for l in self.data_handle.data_list ])
+            sched = Scheduler(load_the_data, [ (l, self.kernels.precision, self.flag_blank_dc) for l in self.data_handle.data_list ])
             for dl in self.data_handle.data_list:
                 search_coarse_channel(dl, self, dataloader=sched, filewriter=filewriter, logwriter=logwriter)
         # Run Parallel version via dask
@@ -216,15 +218,15 @@ def search(self, n_partitions=1, progress_bar='n'):
         t1 = time.time()
         self.last_logwriter(path_log, '\n===== Search time: {:.2f} minutes'.format((t1 - t0)/60.0))
 
-def load_the_data(d, precision):
+def load_the_data(d, precision, flag_blank_dc):
     data_obj = DATAH5(d['filename'],
                   f_start=d['f_start'],
                   f_stop=d['f_stop'],
                   coarse_chan=d['coarse_chan'],
                   n_coarse_chan=d['n_coarse_chan'],
                   gpu_backend=False,
                   precision=precision)
-    spectra, drift_indices = data_obj.load_data()
+    spectra, drift_indices = data_obj.load_data(flag_blank_dc=flag_blank_dc)
     data_obj.close()
 
     return (data_obj, spectra, drift_indices)
@@ -277,7 +279,7 @@ def search_coarse_channel(data_dict, find_doppler_instance, dataloader=None, log
     if dataloader:
         data_obj, spectra, drift_indices = dataloader.get()
     else:
-        data_obj, spectra, drift_indices = load_the_data(d, fd.kernels.precision)
+        data_obj, spectra, drift_indices = load_the_data(d, fd.kernels.precision, fd.flag_blank_dc)
 
     fileroot_out = filename_in.split('/')[-1].replace('.h5', '').replace('.fits', '').replace('.fil', '')
     if logwriter is None:

turbo_seti/find_doppler/helper_functions.py

@@ -1,6 +1,38 @@
 #!/usr/bin/env python
 
 import numpy as np
+import logging
+
+
+def cut_the_mid_spike(data_array, n_coarse_chan):
+    """ Cut the mid-point spike in coarse channels.
+
+    Removes the DC spike in the centre of each coarse channel bin.
+
+    Parameters
+    ----------
+    data_array : ndarray
+        Full data array.
+    n_coarse_chan : int
+        Number of coarse channels.
+    """
+    logger = logging.getLogger('cut_the_mid_spike')
+    if not type(n_coarse_chan) != "int":
+        logger.error("Number of coarse channels is not an integer, no action taken!")
+        return
+    if n_coarse_chan < 1:
+        logger.error = "Number of coarse channels < 1, no action taken!"
+        return
+
+    n_fine_chan = data_array.shape[-1]
+    n_fine_chan_per_coarse_chan = int(n_fine_chan / n_coarse_chan) # ratio of fine channels to coarse channels
+
+    mid_chan = int(n_fine_chan_per_coarse_chan / 2)
+
+    for ii in range(n_coarse_chan):
+        ss = ii * n_fine_chan_per_coarse_chan
+        # Replace the mid point value with the neighbour's value.
+        data_array[..., ss + mid_chan] = data_array[..., ss + mid_chan + 1]
 
 
 def chan_freq(header, fine_channel, tdwidth, ref_frame):

turbo_seti/find_doppler/seti_event.py

@@ -54,6 +54,8 @@ def main(args=None):
                    help='Use a progress bar with dask? (y/n)')
     p.add_argument('-g', '--gpu', dest='flag_gpu', type=str, default='n',
                    help='Compute on the GPU? (y/n)')
+    p.add_argument('-z', '--blank_dc', dest='flag_blank_dc', type=str, default='y',
+                   help='Smooth out the DC spike? (y/n)')
     p.add_argument('-d', '--gpu_id', dest='gpu_id', type=int, default=0,
                    help='Use which GPU device? (0,1,...)')
     p.add_argument('-P', '--profile', dest='flag_profile', type=str, default='n',
@@ -132,6 +134,7 @@ def exec_proc(args):
                                   n_coarse_chan=args.n_coarse_chan,
                                   gpu_backend=(args.flag_gpu == "y"),
                                   gpu_id=args.gpu_id,
+                                  blank_dc=(args.flag_blank_dc == "y"),
                                   precision=1 if args.flag_single_precision == "y" else 2,
                                   log_level_int=log_level_int)