Migrate A-Team-Clipping towards DP3 clipper step

pyproject.toml

@@ -50,7 +50,6 @@ dynamic = ["version"]
 [tool.setuptools]
 packages = []
 script-files = [
-    "scripts/Ateamclipper.py",
     "scripts/BLsmooth.py",
     "scripts/add_missing_stations.py",
     "scripts/blank_image_reg.py",

scripts/Ateamclipper.py

-#!/usr/bin/env python
-
-## changelog
-# W.Williams 2014/11/03  add - to give input/output statistics per channel
-# W.Williams 2014/11/03  fix - statistics per correlation
-# A.Drabent 2019/07/24   write fraction of flagged data into output file (for prefactor/LINC)
-
-import numpy
-import pyrap.tables as pt
-import os
-import sys
-
-msname = str(sys.argv[1])
-
-cliplevelhba = 5.0
-cliplevellba = 50.0
-
-t = pt.table(msname, readonly=False)
-data = t.getcol('MODEL_DATA')
-flag = t.getcol('FLAG')
-freq_tab= pt.table(msname + '::SPECTRAL_WINDOW')
-freq    = freq_tab.getcol('REF_FREQUENCY')
-
-if freq[0] > 100e6:
- cliplevel = cliplevelhba
-if freq[0] < 100e6:
- cliplevel = cliplevellba
-
-
-print('------------------------------')
-print('SB Frequency [MHz]', freq[0]/1e6)
-for chan in range(0,numpy.size(data[0,:,0])):
-  print('chan %i : %.5f%% input XX flagged' %( chan, 100.*numpy.sum(flag[:,chan,0] == True)/numpy.size(flag[:,chan,0]) ))
-  print('chan %i : %.5f%% input YY flagged' %( chan, 100.*numpy.sum(flag[:,chan,3] == True)/numpy.size(flag[:,chan,3]) ))
-input_flags_xx = 100. * numpy.sum(flag[:,:,0] == True)/numpy.size(flag[:,:,0])
-input_flags_yy = 100. * numpy.sum(flag[:,:,3] == True)/numpy.size(flag[:,:,3])
-print('Total : %.5f%% input XX flagged' %(   input_flags_xx ))
-print('Total : %.5f%% input YY flagged' %(   input_flags_yy ))
-print('')
-print('Cliplevel used [Jy]', cliplevel)
-print('\n\n')
-
-for pol in range(0,numpy.size(data[0,0,:])):
- for chan in range(0,numpy.size(data[0,:,0])):
-  print('Doing polarization,chan', pol, chan)
-  idx = numpy.where(abs(data[:,chan,pol]) > cliplevel)
-  flag[idx,chan,0] = True
-  flag[idx,chan,1] = True
-  flag[idx,chan,2] = True
-  flag[idx,chan,3] = True 
-
-print('')
-for chan in range(0,numpy.size(data[0,:,0])):
-  print('chan %i : %.5f%% output XX flagged' %( chan, 100.*numpy.sum(flag[:,chan,0] == True)/numpy.size(flag[:,chan,0]) ))
-  print('chan %i : %.5f%% output YY flagged' %( chan, 100.*numpy.sum(flag[:,chan,3] == True)/numpy.size(flag[:,chan,3]) ))
-output_flags_xx = 100. * numpy.sum(flag[:,:,0] == True)/numpy.size(flag[:,:,0])
-output_flags_yy = 100. * numpy.sum(flag[:,:,3] == True)/numpy.size(flag[:,:,3])
-print('Total : %.5f%% output XX flagged' %(   output_flags_xx ))
-print('Total : %.5f%% output YY flagged' %(   output_flags_yy ))
-print('')
-t.putcol('FLAG', flag)
-t.close()
-freq_tab.close()
-
-os.system('echo ' + str(freq[0]) + ' ' + str(output_flags_xx - input_flags_xx) + ' ' + str(output_flags_yy - input_flags_yy) + ' >> Ateamclipper.txt')

scripts/plot_Ateamclipper.py

@@ -2,48 +2,92 @@
 # -* coding: utf-8 -*-
 
 """
-Adds phases (=0) and amplitudes (=1) to any missing station if they appear in an h5parm, but not in a particular soltab. 
-
-Created on Tue Jul 24 2019
+Created on Wed Mar 12, 2025
 
 @author: Alexander Drabent
 """
 
-import argparse
+import numpy
+import multiprocessing
 import matplotlib as mpl
-mpl.use('Agg')
-import matplotlib
+
+mpl.use("Agg")
 import matplotlib.pyplot as plt
-import numpy
+import casacore.tables as ctab
+
+
+def process_table(table):
+    tt = ctab.table(table, readonly=True)
+    freqs = tt.getcol("Frequency")
+    flag_perc = tt.getcol("Percentage")
+    flags_per_freq = {}
+    for i, freq in enumerate(freqs):
+        flags_per_freq[freq] = flag_perc[i]
+    tt.close()
+    return flags_per_freq
+
+
+def main(flagfreq_before, flagfreq_after, outfile="Ateamclipper.png"):
+
+    ## reading in all flagging percentages
+    pool = multiprocessing.Pool(processes=int(multiprocessing.cpu_count() / 2))
+    flag_perc_before = {
+        freq: flag_perc
+        for entry in pool.map(process_table, flagfreq_before)
+        for freq, flag_perc in entry.items()
+    }
+    flag_perc_after = {
+        freq: flag_perc
+        for entry in pool.map(process_table, flagfreq_after)
+        for freq, flag_perc in entry.items()
+    }
+
+    freq_list = list(flag_perc_before.keys() | flag_perc_after.keys())
+    flag_perc_diff = {
+        freq: flag_perc_after[freq] - flag_perc_before[freq] for freq in freq_list
+    }
+
+    perc_list = [flag_perc_diff[freq] for freq in freq_list]
 
-def main(txtfile = 'Ateamclipper.txt', outfile = 'Ateamclipper.png'):
-
-    frac_list_xx = []
-    frac_list_yy = []
-    freq_list    = []
-    with open(txtfile, 'r') as infile:
-        for line in infile:
-            freq_list.append(float(line.split()[0]))
-            frac_list_xx.append(float(line.split()[1]))
-            frac_list_yy.append(float(line.split()[2]))
-    
     # Plot the amount of clipped data vs. frequency potentially contaminated by the A-team
-    plt.scatter(numpy.array(freq_list) / 1e6, numpy.array(frac_list_xx), marker = '.', s = 10)
-    plt.xlabel('frequency [MHz]')
-    plt.ylabel('A-team clipping fraction [%]')
+    plt.scatter(numpy.array(freq_list) / 1e6, numpy.array(perc_list), marker=".", s=10)
+    plt.xlabel("frequency [MHz]")
+    plt.ylabel("A-Team clipping fraction [%]")
     plt.savefig(outfile)
-    return(0)
+    return 0
+
 
 if __name__ == "__main__":
-    parser = argparse.ArgumentParser(description='Adds phases and amplitudes to any missing station if they appear in an h5parm, but not in a particular soltab.')
+    import argparse
 
-    parser.add_argument('txtfile', type=str,
-                        help='Input text file containing frequency and flag fraction of the XX and YY polarization.')
-    parser.add_argument('outfile', type=str,
-                        help='Output image file containing frequency and flag fraction of the XX and YY polarization.')
+    parser = argparse.ArgumentParser(
+        description="Compare the flagging percentage of two given sets/lists of flagfreq tables provided by the DP3 count step"
+    )
 
+    parser.add_argument(
+        "--flagfreq_before",
+        type=str,
+        nargs="+",
+        help="One or more flagfreq tables before running flagging",
+    )
+    parser.add_argument(
+        "--flagfreq_after",
+        type=str,
+        nargs="+",
+        help="One or more flagfreq tables after  running flagging",
+    )
+    parser.add_argument(
+        "--outfile",
+        type=str,
+        nargs="?",
+        help="Output file name for the plot",
+        default="Ateamclipper.png",
+    )
 
     args = parser.parse_args()
 
-    main(txtfile = args.txtfile, outfile = args.outfile)
-
+    main(
+        flagfreq_before=args.flagfreq_before,
+        flagfreq_after=args.flagfreq_after,
+        outfile=args.outfile,
+    )

steps/Ateamclipper.cwl

-class: CommandLineTool
-cwlVersion: v1.2
-id: check_ateam_separation
-baseCommand:
-  - Ateamclipper.py
-inputs:
-  - id: msin
-    type:
-      - Directory
-      - type: array
-        items: Directory
-    inputBinding:
-      position: 0
-    doc: Input measurement set
-outputs:
-  - id: msout
-    doc: Output MS
-    type: Directory
-    outputBinding:
-      glob: $(inputs.msin.basename)
-  - id: logfile
-    type: File[]
-    outputBinding:
-      glob: Ateamclipper.log
-  - id: output
-    type: File
-    outputBinding:
-      glob: Ateamclipper.txt
-label: Ateamclipper
-hints:
-  - class: InitialWorkDirRequirement
-    listing:
-      - entry: $(inputs.msin)
-        writable: true
-  - class: InplaceUpdateRequirement
-    inplaceUpdate: true
-  - class: DockerRequirement
-    dockerPull: astronrd/linc
-  - class: InlineJavascriptRequirement
-  - class: ResourceRequirement
-    coresMin: 8
-stdout: Ateamclipper.log

steps/filter_predict.cwl → steps/clipper.cwl

 class: CommandLineTool
 cwlVersion: v1.2
-id: predict
-baseCommand:
-  - DP3
+id: Ateamclipper
+label: clip A-team
+doc: |
+    Simulates data for the A-team sources based off a skymodel,
+    and flags the visibilities of the input MeasuremenSet where
+    the model data exceeds the threshold for LBA (50 janskys)
+    or HBA (5 janksys).
+
+baseCommand: DP3
+arguments:
+  - steps=[filter,count1,clipper,count2]
+  - msout=.
+  - count1.save=true
+  - count2.save=true
+
 inputs:
-  - id: max_dp3_threads
-    type: int?
-    inputBinding:
-      position: 0
-      prefix: numthreads=
-      separate: false
   - id: msin
     type: Directory
     inputBinding:
       position: 0
       prefix: msin=
       separate: false
-    doc: Input Measurement Set
-  - default: DATA
-    id: msin_datacolumn
+      shellQuote: false
+    doc: Input data in MeasurementSet format.
+  - id: msin_datacolumn
     type: string?
+    default: DATA
     inputBinding:
       position: 0
       prefix: msin.datacolumn=
       separate: false
-    doc: Input data Column
-  - default: MODEL_DATA
-    id: msout_datacolumn
+      shellQuote: false
+    doc: |
+        Data column of the MeasurementSet
+        from which input data is read.
+
+  - id: storagemanager
     type: string?
+    default: ""
     inputBinding:
-      position: 0
-      prefix: msout.datacolumn=
+       prefix: msout.storagemanager=
+       separate: false
+  - id: databitrate
+    type: int?
+    default: 0
+    inputBinding:
+       prefix: msout.storagemanager.databitrate=
+       separate: false
+
+  - id: operation
+    type:
+      type: enum
+      symbols:
+        - replace
+    doc: |
+      Type of operation to be performed on clipped data.
+    inputBinding:
+      prefix: clipper.operation=
       separate: false
-  - id: sources_db
+
+  - id: sourcedb
     type:
       - File
       - Directory
     inputBinding:
       position: 0
-      prefix: predict.sourcedb=
+      prefix: clipper.sourcedb=
       separate: false
   - default: null
     id: sources
     type: string[]?
     inputBinding:
       position: 0
-      prefix: predict.sources=
+      prefix: clipper.sources=
       separate: false
       itemSeparator: ','
       valueFrom: "[$(self.join(','))]"
-  - default: false
-    id: usebeammodel
-    type: boolean?
-    inputBinding:
-      position: 0
-      prefix: predict.usebeammodel=True
+    doc: |
+        Labels of the skymodel patches to
+        use to simulate visibilities.
   - id: usechannelfreq
-    default: true
+    default: false
     type: boolean?
     inputBinding:
       valueFrom: $(!self)
       position: 0
-      prefix: predict.usechannelfreq=False
-  - default: false
-    id: onebeamperpatch
+      prefix: clipper.usechannelfreq=False
+      separate: false
+  - id: usebeammodel
     type: boolean?
+    default: true
     inputBinding:
       position: 0
-      prefix: predict.onebeamperpatch=True
-  - default: null
-    id: filter_baselines
+      prefix: clipper.usebeammodel=True
+      shellQuote: false
+    doc: |
+        Determines whether to use the beam model.
+  - id: beamproximitylimit
+    doc: |
+         Specified in arcseconds, and if non-zero, sources that are near each other are clustered
+         and the beam is only calculated for each cluster.
+    type: int?
+    default: 2000
+    inputBinding:
+      position: 0
+      prefix: clipper.beamproximitylimit=
+      separate: false
+
+  - id: filter_baselines
     type: string?
     inputBinding:
       position: 0
       prefix: filter.baseline=
       separate: false
       valueFrom: $(self)
-  - default: false
-    id: filter_remove
+  - id: filter_remove
+    default: false
     type: boolean?
     inputBinding:
       position: 0
       prefix: filter.remove=True
-  - id: writefullresflag
-    type: boolean?
-    default: false
-    inputBinding:
-       prefix: msout.writefullresflag=True
-  - default: default
-    id: beammode
-    type: string?
-    inputBinding:
-      position: 0
-      prefix: predict.beammode=
-      separate: false
-  - id: overwrite
-    type: boolean?
-    default: false
-    inputBinding:
-       prefix: msout.overwrite=True
-  - id: storagemanager
-    type: string?
-    default: ""
-    inputBinding:
-       prefix: msout.storagemanager=
-       separate: false
-  - id: databitrate
+
+  - id: max_dp3_threads
     type: int?
     inputBinding:
-       prefix: msout.storagemanager.databitrate=
-       separate: false
+      prefix: numthreads=
+      separate: false
+    doc: The number of threads per DP3 process.
+
 outputs:
   - id: msout
-    doc: Output Measurement Set
+    doc: Output data in MeasurementSet format.
     type: Directory
     outputBinding:
       glob: $(inputs.msin.basename)
   - id: logfile
     type: File[]
     outputBinding:
-      glob: 'filter_predict*.log'
-arguments:
-  - steps=[filter,predict,count]
-  - msout=.
+      glob: clipper*.log
+    doc: |
+        The files containing the stdout
+        and stderr from the step.
+  - id: flagfreq_before
+    type: Directory
+    outputBinding:
+      glob: '$(inputs.msin.nameroot.split(".")[0])_count1.flagfreq'
+    doc: A MS file containing flagging fraction statistics before clipping.
+  - id: flagfreq_after
+    type: Directory
+    outputBinding:
+      glob: '$(inputs.msin.nameroot.split(".")[0])_count2.flagfreq'
+    doc: A MS file containing flagging fraction statistics after clipping.
+
 requirements:
+  - class: InplaceUpdateRequirement
+    inplaceUpdate: true 
   - class: InitialWorkDirRequirement
     listing:
       - entry: $(inputs.msin)
         writable: true
-  - class: InplaceUpdateRequirement
-    inplaceUpdate: true
   - class: InlineJavascriptRequirement
   - class: ResourceRequirement
     coresMin: $(inputs.max_dp3_threads)
 hints:
   - class: DockerRequirement
     dockerPull: astronrd/linc
-stdout: filter_predict.log
-stderr: filter_predict_err.log
+
+stdout: clipper.log
+stderr: clipper_err.log

steps/plot_Ateamclipper.cwl

@@ -4,15 +4,22 @@ id: plot_Ateamclipper
 baseCommand:
   - plot_Ateamclipper.py
 inputs:
-  - id: clipper_output
-    type: File?
+  - id: flagfreq_before
+    type: Directory[]
     inputBinding:
       position: 1
+      prefix: '--flagfreq_before'
+  - id: flagfreq_after
+    type: Directory[]
+    inputBinding:
+      position: 1
+      prefix: '--flagfreq_after'
   - id: outputimage
     type: string?
     default: Ateamclipper.png
     inputBinding:
       position: 2
+      prefix: '--outfile'
 outputs:
   - id: output_imag
     doc: Output image
@@ -22,4 +29,4 @@ outputs:
 label: plot_Ateamclipper
 hints:
   - class: DockerRequirement
-    dockerPull: astronrd/linc
+    dockerPull: astronrd/linc
\ No newline at end of file

workflows/linc_target/dp3_prep_targ.cwl

@@ -105,11 +105,6 @@ outputs:
     outputSource:
       - concat_logfiles_dp3/output
     type: File
-  - id: predict_logfile
-    outputSource:
-      - concat_logfiles_predict/output
-    type: File
-    pickValue: all_non_null
   - id: clipper_logfile
     outputSource:
       - concat_logfiles_clipper/output
@@ -117,14 +112,19 @@ outputs:
     pickValue: all_non_null
   - id: msout
     outputSource:
-      - Ateamclipper/msout
+      - clipper/msout
       - dp3_execute/msout
     type: Directory
     pickValue: first_non_null
-  - id: clipper_output
+  - id: clipper_flags_before
     outputSource:
-      - Ateamclipper/output
-    type: File
+      - clipper/flagfreq_before
+    type: Directory
+    pickValue: all_non_null
+  - id: clipper_flags_after
+    outputSource:
+      - clipper/flagfreq_after
+    type: Directory
     pickValue: all_non_null
   - id: instrument_tables
     outputSource:
@@ -198,26 +198,11 @@ steps:
       - id: output
     run: ../../steps/concatenate_files.cwl
     label: concat_logfiles_dp3
-  - id: concat_logfiles_predict
-    in:
-      - id: file_list
-        source:
-          - predict/logfile
-        pickValue: all_non_null
-      - id: file_prefix
-        default: predict_targ
-      - id: execute
-        source: clipAteam
-    out:
-      - id: output
-    run: ../../steps/concatenate_files.cwl
-    when: $(inputs.execute)
-    label: concat_logfiles_predict
   - id: concat_logfiles_clipper
     in:
       - id: file_list
         source:
-          - Ateamclipper/logfile
+          - clipper/logfile
         pickValue: all_non_null
       - id: file_prefix
         default: Ateamclipper
@@ -263,53 +248,45 @@ steps:
       - id: logfile
     run: ../../steps/dp3_prep_target.cwl
     label: DP3.Execute
-  - id: predict
+  - id: clipper
     in:
-      - id: max_dp3_threads
-        source: max_dp3_threads
       - id: msin
         source: dp3_execute/msout
       - id: msin_datacolumn
         default: DATA
-      - id: msout_datacolumn
-        default: MODEL_DATA
-      - id: sources_db
+      - id: storagemanager
+        default: Dysco
+      - id: databitrate
+        default: 0
+      - id: operation
+        default: replace
+      - id: sourcedb
         source: skymodel
       - id: sources
         source:
           - clip_sources
       - id: usebeammodel
         default: true
-      - id: storagemanager
-        default: Dysco
-      - id: databitrate
-        default: 0
+      - id: usechannelfreq
+        default: false
+      - id: beamproximitylimit
+        default: 2000
       - id: filter_baselines
         source: process_baselines_target
-      - id: usechannelfreq
+      - id: filter_remove
         default: false
+      - id: max_dp3_threads
+        source: max_dp3_threads
       - id: execute
         source: clipAteam
     out:
       - id: msout
       - id: logfile
-    run: ../../steps/filter_predict.cwl
-    when: $(inputs.execute)
-  - id: Ateamclipper
-    in:
-      - id: msin
-        source:
-          - predict/msout
-        pickValue: all_non_null
-      - id: execute
-        source: clipAteam
-    out:
-      - id: msout
-      - id: logfile
-      - id: output
-    run: ../../steps/Ateamclipper.cwl
+      - id: flagfreq_before
+      - id: flagfreq_after
+    run: ../../steps/clipper.cwl
     when: $(inputs.execute)
-    label: Ateamclipper
+
 requirements:
   - class: InlineJavascriptRequirement
   - class: StepInputExpressionRequirement

workflows/linc_target/prep.cwl

@@ -167,7 +167,6 @@ outputs:
       - concat_logfiles_stationlist/output
       - concat_logfiles_RMextract/output
       - concat_logfiles_prep_targ/output
-      - concat_logfiles_predict_targ/output
       - concat_logfiles_clipper_targ/output
       - concat_logfiles_plot_demix/output
     type: File[]
@@ -476,10 +475,10 @@ steps:
       - id: prep_flags_out
       - id: initial_flags_out
       - id: prep_logfile
-      - id: predict_logfile
       - id: clipper_logfile
       - id: msout
-      - id: clipper_output
+      - id: clipper_flags_before
+      - id: clipper_flags_after
       - id: instrument_tables
     run: ./dp3_prep_targ.cwl
     label: dp3_prep_target
@@ -487,15 +486,17 @@ steps:
       - msin
   - id: plot_Ateamclipper
     in:
-      - id: clipper_output
-        source: concat_logfiles_clipper_output/output
+      - id: flagfreq_before
+        source: dp3_prep_target/clipper_flags_before
+      - id: flagfreq_after
+        source: dp3_prep_target/clipper_flags_after
       - id: execute
         source: clipAteam
     out:
       - id: output_imag
     run: ../../steps/plot_Ateamclipper.cwl
     when: $(inputs.execute)
-    label: concat_logfiles_clipper_output
+    label: plot_Ateamclipper
   - id: plot_demix
     in:
       - id: instrument_tables
@@ -508,23 +509,6 @@ steps:
     run: ../../steps/plot_demix.cwl
     label: plot_demix_solutions
     when: $(inputs.demix)
-  - id: concat_logfiles_clipper_output
-    in:
-      - id: file_list
-        linkMerge: merge_flattened
-        source:
-          - dp3_prep_target/clipper_output
-      - id: file_prefix
-        default: Ateamclipper
-      - id: file_suffix
-        default: txt
-      - id: execute
-        source: clipAteam
-    out:
-      - id: output
-    run: ../../steps/concatenate_files.cwl
-    when: $(inputs.execute)
-    label: concat_logfiles_clipper_output
   - id: concat_logfiles_prep_targ
     in:
       - id: file_list
@@ -537,21 +521,6 @@ steps:
       - id: output
     run: ../../steps/concatenate_files.cwl
     label: concat_logfiles_prep_target
-  - id: concat_logfiles_predict_targ
-    in:
-      - id: file_list
-        linkMerge: merge_flattened
-        source:
-          - dp3_prep_target/predict_logfile
-      - id: file_prefix
-        default: predict_targ
-      - id: execute
-        source: clipAteam
-    out:
-      - id: output
-    run: ../../steps/concatenate_files.cwl
-    when: $(inputs.execute)
-    label: concat_logfiles_predict_targ
   - id: concat_logfiles_clipper_targ
     in:
       - id: file_list