Debugged the pipeline after running unreduced data.

The changes include:
* Core allocation of computationally heavy steps
* Added RFI strategies, since these may not exist within the container image
* Updated python scripts to reflect changes in the generic pipeline implementation

rfistrategies/lofar-default.lua

+--[[
+ This is the LOFAR strategy. It is almost equal to the
+ generic "minimal" AOFlagger strategy, version 2020-06-14.
+ Author: André Offringa
+]]--
+
+function execute(input)
+
+  --
+  -- Generic settings
+  --
+
+  local base_threshold = 1.0  -- lower means more sensitive detection
+  -- How to flag complex values, options are: phase, amplitude, real, imaginary, complex
+  local representation = "amplitude"
+  local iteration_count = 3  -- how many iterations to perform?
+  local threshold_factor_step = 2.0 -- How much to increase the sensitivity each iteration?
+  local frequency_resize_factor = 3.0 -- Amount of "extra" smoothing in frequency direction
+  local transient_threshold_factor = 1.0 -- decreasing this value makes detection of transient RFI more aggressive
+
+  --
+  -- End of generic settings
+  --
+
+  local inpPolarizations = input:get_polarizations()
+
+  input:clear_mask()
+
+  -- For collecting statistics. Note that this is done after clear_mask(),
+  -- so that the statistics ignore any flags in the input data. 
+  local copy_of_input = input:copy()
+  
+  for ipol,polarization in ipairs(inpPolarizations) do
+ 
+    local data = input:convert_to_polarization(polarization)
+
+    data = data:convert_to_complex(representation)
+    local original_data = data:copy()
+
+    for i=1,iteration_count-1 do
+      local threshold_factor = threshold_factor_step^(iteration_count-i)
+
+      local sumthr_level = threshold_factor * base_threshold
+      aoflagger.sumthreshold(data, sumthr_level, sumthr_level*transient_threshold_factor, true, true)
+ 
+      -- Do timestep & channel flagging
+      local chdata = data:copy()
+      aoflagger.threshold_timestep_rms(data, 3.5)
+      aoflagger.threshold_channel_rms(chdata, 3.0 * threshold_factor, true)
+      data:join_mask(chdata)
+
+      -- High pass filtering steps
+      data:set_visibilities(original_data)
+      local resized_data = aoflagger.downsample(data, 3, frequency_resize_factor, true)
+      aoflagger.low_pass_filter(resized_data, 21, 31, 2.5, 5.0)
+      aoflagger.upsample(resized_data, data, 3, frequency_resize_factor)
+
+      local tmp = original_data - data
+      tmp:set_mask(data)
+      data = tmp
+
+      aoflagger.set_progress((ipol-1)*iteration_count+i, #inpPolarizations*iteration_count )
+    end -- end of iterations
+
+    aoflagger.sumthreshold(data, base_threshold, base_threshold*transient_threshold_factor, true, true)
+
+    if input:is_complex() then
+      data = data:convert_to_complex("complex")
+    end
+    input:set_polarization_data(polarization, data)
+
+    aoflagger.set_progress(ipol, #inpPolarizations )
+  end -- end of polarization iterations
+
+  aoflagger.scale_invariant_rank_operator(input, 0.2, 0.2)
+  aoflagger.threshold_timestep_rms(input, 4.0)
+  input:flag_nans()
+end
+

rfistrategies/lofar-lba-wideband.lua

+--[[
+ This is a LBA AOFlagger strategy for combined subbands, version 2021-03-30
+ Author: André Offringa
+]]--
+
+aoflagger.require_min_version("3.1")
+
+function execute(input)
+
+  --
+  -- Generic settings
+  --
+
+  -- What polarizations to flag? Default: input:get_polarizations() (=all that are in the input data)
+  -- Other options are e.g.:
+  -- { 'XY', 'YX' } to flag only XY and YX, or
+  -- { 'I', 'Q' } to flag only on Stokes I and Q
+  local flag_polarizations = input:get_polarizations()
+
+  local base_threshold = 1.2  -- lower means more sensitive detection
+  -- How to flag complex values, options are: phase, amplitude, real, imaginary, complex
+  -- May have multiple values to perform detection multiple times
+  local flag_representations = { "amplitude" }
+  local iteration_count = 5  -- how many iterations to perform?
+  local threshold_factor_step = 2.0 -- How much to increase the sensitivity each iteration?
+  -- If the following variable is true, the strategy will consider existing flags
+  -- as bad data. It will exclude flagged data from detection, and make sure that any existing
+  -- flags on input will be flagged on output. If set to false, existing flags are ignored.
+  local exclude_original_flags = true
+  local transient_threshold_factor = 1.0 -- decreasing this value makes detection of transient RFI more aggressive
+ 
+  --
+  -- End of generic settings
+  --
+
+  local inpPolarizations = input:get_polarizations()
+
+  if(not exclude_original_flags) then
+    input:clear_mask()
+  end
+  -- For collecting statistics. Note that this is done after clear_mask(),
+  -- so that the statistics ignore any flags in the input data. 
+  local copy_of_input = input:copy()
+
+  aoflagger.normalize_bandpass(input)
+  
+  for ipol,polarization in ipairs(flag_polarizations) do
+ 
+    local pol_data = input:convert_to_polarization(polarization)
+    local original_data
+
+    for _,representation in ipairs(flag_representations) do
+
+      data = pol_data:convert_to_complex(representation)
+      original_data = data:copy()
+
+      for i=1,iteration_count-1 do
+        local threshold_factor = math.pow(threshold_factor_step, iteration_count-i)
+
+        local sumthr_level = threshold_factor * base_threshold
+        if(exclude_original_flags) then
+          aoflagger.sumthreshold_masked(data, original_data, sumthr_level, sumthr_level*transient_threshold_factor, true, true)
+        else
+          aoflagger.sumthreshold(data, sumthr_level, sumthr_level*transient_threshold_factor, true, true)
+        end
+
+        -- Do timestep & channel flagging
+        local chdata = data:copy()
+        aoflagger.threshold_timestep_rms(data, 3.5)
+        aoflagger.threshold_channel_rms(chdata, 3.0 * threshold_factor, true)
+        data:join_mask(chdata)
+
+        -- High pass filtering steps
+        data:set_visibilities(original_data)
+        if(exclude_original_flags) then
+          data:join_mask(original_data)
+        end
+
+        aoflagger.low_pass_filter(data, 21, 31, 1.0, 1.0)
+
+        -- In case this script is run from inside rfigui, calling
+        -- the following visualize function will add the current result
+        -- to the list of displayable visualizations.
+        -- If the script is not running inside rfigui, the call is ignored.
+        aoflagger.visualize(data, "Fit #"..i, i-1)
+
+        local tmp = original_data - data
+        tmp:set_mask(data)
+        data = tmp
+
+        aoflagger.visualize(data, "Residual #"..i, i+iteration_count)
+        aoflagger.set_progress((ipol-1)*iteration_count+i, #flag_polarizations*iteration_count )
+      end -- end of iterations
+
+      if(exclude_original_flags) then
+        aoflagger.sumthreshold_masked(data, original_data, base_threshold, base_threshold*transient_threshold_factor, true, true)
+      else
+        aoflagger.sumthreshold(data, base_threshold, base_threshold*transient_threshold_factor, true, true)
+      end
+    end -- end of complex representation iteration
+
+    if(exclude_original_flags) then
+      data:join_mask(original_data)
+    end
+
+    -- Helper function used below
+    function contains(arr, val)
+      for _,v in ipairs(arr) do
+        if v == val then return true end
+      end
+      return false
+    end
+
+    if contains(inpPolarizations, polarization) then
+      if input:is_complex() then
+        data = data:convert_to_complex("complex")
+      end
+      input:set_polarization_data(polarization, data)
+    else
+      input:join_mask(data)
+    end
+
+    aoflagger.visualize(data, "Residual #"..iteration_count, 2*iteration_count)
+    aoflagger.set_progress(ipol, #flag_polarizations )
+  end -- end of polarization iterations
+
+  if(exclude_original_flags) then
+    aoflagger.scale_invariant_rank_operator_masked(input, copy_of_input, 0.2, 0.2)
+  else
+    aoflagger.scale_invariant_rank_operator(input, 0.2, 0.2)
+  end
+
+  aoflagger.threshold_timestep_rms(input, 4.0)
+
+  if input:is_complex() and input:has_metadata() then
+    -- This command will calculate a few statistics like flag% and stddev over
+    -- time, frequency and baseline and write those to the MS. These can be
+    -- visualized with aoqplot.
+    aoflagger.collect_statistics(input, copy_of_input)
+  end
+  input:flag_nans()
+end
+

scripts/TargetListToCoords.py

@@ -27,8 +27,8 @@ def plugin_main(**kwargs):
     
     # read in the catalogue to get source_id, RA, and DEC
     t = Table.read(target_file, format='csv')
-    RA_val = t['RA_LOTSS'].data[0]
-    DEC_val = t['DEC_LOTSS'].data[0]
+    RA_val = t['RA'].data[0]
+    DEC_val = t['DEC'].data[0]
     Source_id = t['Source_id'].data[0]
     if str(Source_id)[0:1] == 'I':
         pass

scripts/check_Ateam_separation.py

-#!/usr/bin/python
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 
 #

scripts/skynet.py

@@ -4,7 +4,7 @@ import subprocess
 import shutil
 import sys
 import glob
-import regex as re
+import re
 from astropy.coordinates import SkyCoord
 from astropy.table import Table
 

steps/Ateamclipper.cwl

 class: CommandLineTool
 cwlVersion: v1.2
-id: check_ateam_separation
+id: Ateamclipper
 label: Ateamclipper
 
 baseCommand:
@@ -36,8 +36,13 @@ hints:
     listing:
       - entry: $(inputs.msin)
         writable: true
+  - class: InplaceUpdateRequirement
+    inplaceUpdate: true
   - class: DockerRequirement
     dockerPull: vlbi-cwl
   - class: InlineJavascriptRequirement
+  - class: ResourceRequirement
+    coresMin: 12
 
 stdout: Ateamclipper.log
+stderr: Ateamclipper_err.log

steps/aoflagger.cwl

@@ -25,7 +25,7 @@ inputs:
       doc: Input data Column
     - id: memoryperc
       type: int?
-      default: 15
+      default: 10
       inputBinding:
         position: 0
         prefix: aoflagger.memoryperc=
@@ -72,6 +72,8 @@ requirements:
       - entry: $(inputs.msin)
         writable: true
   - class: ShellCommandRequirement
+  - class: ResourceRequirement
+    coresMin: 6
 
 hints:
   DockerRequirement:

steps/delay_solve.cwl

@@ -26,6 +26,10 @@ outputs:
       type: File[]
       outputBinding:
         glob: merged_selfcal*.h5
+    - id: images
+      type: File[]
+      outputBinding:
+        glob: image*.png
     - id: logfile
       type: File[]
       outputBinding:

steps/download_cats.cwl

@@ -7,28 +7,28 @@ baseCommand:
     - downloadCats.py
 
 inputs:
-    - id: msin # mapfile_in = kwargs['mapfile_in']
+    - id: msin
       type: Directory[]
       inputBinding:
         position: 0
     - id: lotss_radius
       type: float?
-      default: 2.0
+      default: 1.5
     - id: lbcs_radius
       type: float?
-      default: 2.0
+      default: 1.5
+    - id: im_radius
+      type: float?
+      default: 1.24
     - id: bright_limit_Jy
       type: float
       default: 5.0
     - id: match_tolerance
       type: float?
       default: 5.0
-    - id: subtract_limit
-      type: float?
-      default: 0.5
     - id: image_limit_Jy
       type: float?
-      default: 0.05
+      default: 0.01
     - id: continue_no_lotss
       type: boolean?
       default: true
@@ -47,7 +47,7 @@ inputs:
       default: "image_catalogue.csv"
     - id: delay_catalogue_name
       type: string?
-      default: "delay_catalogue.csv"
+      default: "delay_calibrators.csv"
     - id: subtract_catalogue_name
       type: string?
       default: "subtract_sources.csv"
@@ -65,18 +65,18 @@ outputs:
     #  type: File
     #  outputBinding:
     #    glob: $(inputs.lbcs_skymodel_name)
-    - id: best_delay_catalogue
-      type: File
-      outputBinding:
-        glob: best_delay_*.csv
-    #- id: image_catalogue
+    #- id: best_delay_catalogue
     #  type: File
     #  outputBinding:
-    #    glob: $(inputs.image_catalogue_name)
-    #- id: delay_catalogue
+    #    glob: best_delay_*.csv
+    #- id: image_catalogue
     #  type: File
     #  outputBinding:
-    #    glob: $(inputs.delay_catalogue_name)
+    #    glob: $(inputs.image_catalogue_name)
+    - id: delay_catalogue
+      type: File
+      outputBinding:
+        glob: $(inputs.delay_catalogue_name)
     #- id: subtract_catalogue
     #  type: File
     #  outputBinding:

steps/dp3_concat.cwl

@@ -89,6 +89,8 @@ hints:
     listing:
       - entry: $(inputs.msin)
         writable: true
+  - class: ResourceRequirement
+    coresMin: 6
 
 stdout: dp3_concat.log
 stderr: dp3_concat_err.log

steps/dp3_prep_target.cwl

@@ -70,6 +70,8 @@ arguments:
 
 requirements:
   - class: InlineJavascriptRequirement
+  - class: ResourceRequirement
+    coresMin: 6
 #  - class: InitialWorkDirRequirement
 #    listing:
 #      - entry: $(inputs.msin)

steps/predict.cwl

@@ -109,6 +109,8 @@ outputs:
 hints:
   - class: DockerRequirement
     dockerPull: vlbi-cwl:latest
+  - class: ResourceRequirement
+    coresMin: 6
 
 stdout: predict_ateam.log
 stderr: predict_ateam_err.log

workflows/delay-calibration.cwl

@@ -50,7 +50,7 @@ steps:
           source: phasesol
       out:
         - id: parset
-        - id: best_delay_cats
+        - id: delay_calibrators
         - id: logdir
         - id: msout
         - id: initial_flags
@@ -85,7 +85,7 @@ steps:
         - id: msin
           source: sort-concatenate-flag/msout
         - id: delay_calibrator
-          source: setup/best_delay_cats
+          source: setup/delay_calibrators
         - id: configfile
           source: configfile
         - id: selfcal
@@ -102,6 +102,7 @@ steps:
       out:
         - id: msout
         - id: solutions
+        - id: pictures
         - id: phaseup_flags
         - id: logdir
 #        - id: summary_file
@@ -129,13 +130,21 @@ outputs:
     type: Directory[]
 
   - id: delay_cat
-    outputSource: setup/best_delay_cats
+    outputSource: setup/delay_calibrators
     type: File
 
   - id: logs
     outputSource: store_logs/dir
     type: Directory
 
+  - id: pictures
+    outputSource: phaseup/pictures
+    type: File[]
+
+  - id: solutions
+    outputSource: phaseup/solutions
+    type: File[]
+
 #  - id: summary_file
 #    outputSource: phaseup-concat.cwl
 #    type: File

workflows/phaseup-concat.cwl

@@ -204,8 +204,6 @@ steps:
     in:
       - id: msin
         source: delay_cal_model/msout
-        #source: phaseup_concatenate/msout
-        #valueFrom: $(self[0])
       - id: configfile
         source: configfile
       - id: selfcal
@@ -214,6 +212,7 @@ steps:
         source: h5merger
     out:
       - id: h5parm
+      - id: images
       - id: logfile
     run: ../steps/delay_solve.cwl
     label: delay_solve
@@ -262,7 +261,7 @@ steps:
         linkMerge: merge_flattened
         source:
           - prep_delay/logfile
-          - concat_logfiles_phaseup/output 
+          - concat_logfiles_phaseup/output
           - sort_concatenate/logfile
           - concat_logfiles_phaseup/output
           - delay_cal_model/logfile
@@ -288,6 +287,9 @@ outputs:
   - id: logdir
     outputSource: save_logfiles/dir
     type: Directory
+  - id: pictures
+    type: File[]
+    outputSource: delay_solve/images
 #  - id: summary_file
 #    type: File
 #    outputSource: summary/summary_file

workflows/setup.cwl

@@ -46,7 +46,7 @@ steps:
         - id: msin
           source: msin
       out:
-        - id: best_delay_catalogue
+        - id: delay_catalogue
         - id: logfile
       run: ../steps/download_cats.cwl
       label: download_cats
@@ -154,8 +154,8 @@ outputs:
     - id: logdir
       outputSource: save_logfiles/dir
       type: Directory
-    - id: best_delay_cats
-      outputSource: download_cats/best_delay_catalogue
+    - id: delay_calibrators
+      outputSource: download_cats/delay_catalogue
       type: File
     - id: parset
       outputSource: dp3_make_parset/parset