From 0ade5c255d92f2677c8c059cb3a68f93e5efb0f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Offringa?= <offringa@gmail.com> Date: Thu, 27 Jul 2023 14:11:50 +0000 Subject: [PATCH] Add lofar hba wideband strategy --- data/strategies/lofar-hba-wideband.lua | 155 +++++++++++++++++++++++++ data/strategies/lofar-lba-wideband.lua | 51 ++++---- 2 files changed, 181 insertions(+), 25 deletions(-) create mode 100644 data/strategies/lofar-hba-wideband.lua diff --git a/data/strategies/lofar-hba-wideband.lua b/data/strategies/lofar-hba-wideband.lua new file mode 100644 index 00000000..3717966c --- /dev/null +++ b/data/strategies/lofar-hba-wideband.lua @@ -0,0 +1,155 @@ +--[[ + This is a HBA AOFlagger strategy for combined subbands, version 2023-07-27 + It is designed for the LINC pipeline. + 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 = 3 -- 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() + + for ipol, polarization in ipairs(flag_polarizations) do + local pol_data = input:convert_to_polarization(polarization) + local converted_data + local converted_copy + + for _, representation in ipairs(flag_representations) do + converted_data = pol_data:convert_to_complex(representation) + converted_copy = converted_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 + if exclude_original_flags then + aoflagger.sumthreshold_masked( + converted_data, + converted_copy, + sumthr_level, + sumthr_level * transient_threshold_factor, + true, + true + ) + else + aoflagger.sumthreshold(converted_data, sumthr_level, sumthr_level * transient_threshold_factor, true, true) + end + + -- Do timestep & channel flagging + local chdata = converted_data:copy() + aoflagger.threshold_timestep_rms(converted_data, 3.5) + aoflagger.threshold_channel_rms(chdata, 3.0 * threshold_factor, true) + converted_data:join_mask(chdata) + + -- High pass filtering steps + converted_data:set_visibilities(converted_copy) + if exclude_original_flags then + converted_data:join_mask(converted_copy) + end + + aoflagger.low_pass_filter(converted_data, 21, 31, 2.5, 5.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(converted_data, "Fit #" .. i, i - 1) + + local tmp = converted_copy - converted_data + tmp:set_mask(converted_data) + converted_data = tmp + + aoflagger.visualize(converted_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( + converted_data, + converted_copy, + base_threshold, + base_threshold * transient_threshold_factor, + true, + true + ) + else + aoflagger.sumthreshold(converted_data, base_threshold, base_threshold * transient_threshold_factor, true, true) + end + end -- end of complex representation iteration + + if exclude_original_flags then + converted_data:join_mask(converted_copy) + 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 + converted_data = converted_data:convert_to_complex("complex") + end + input:set_polarization_data(polarization, converted_data) + else + input:join_mask(converted_data) + end + + aoflagger.visualize(converted_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 diff --git a/data/strategies/lofar-lba-wideband.lua b/data/strategies/lofar-lba-wideband.lua index 07125597..01b4c132 100644 --- a/data/strategies/lofar-lba-wideband.lua +++ b/data/strategies/lofar-lba-wideband.lua @@ -45,11 +45,12 @@ function execute(input) for ipol, polarization in ipairs(flag_polarizations) do local pol_data = input:convert_to_polarization(polarization) - local original_data + local converted_data + local converted_copy for _, representation in ipairs(flag_representations) do - data = pol_data:convert_to_complex(representation) - original_data = data:copy() + converted_data = pol_data:convert_to_complex(representation) + converted_copy = converted_data:copy() for i = 1, iteration_count - 1 do local threshold_factor = threshold_factor_step ^ (iteration_count - i) @@ -57,61 +58,61 @@ function execute(input) local sumthr_level = threshold_factor * base_threshold if exclude_original_flags then aoflagger.sumthreshold_masked( - data, - original_data, + converted_data, + converted_copy, sumthr_level, sumthr_level * transient_threshold_factor, true, true ) else - aoflagger.sumthreshold(data, sumthr_level, sumthr_level * transient_threshold_factor, true, true) + aoflagger.sumthreshold(converted_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) + local chdata = converted_data:copy() + aoflagger.threshold_timestep_rms(converted_data, 3.5) aoflagger.threshold_channel_rms(chdata, 3.0 * threshold_factor, true) - data:join_mask(chdata) + converted_data:join_mask(chdata) -- High pass filtering steps - data:set_visibilities(original_data) + converted_data:set_visibilities(converted_copy) if exclude_original_flags then - data:join_mask(original_data) + converted_data:join_mask(converted_copy) end - aoflagger.low_pass_filter(data, 21, 31, 1.0, 1.0) + aoflagger.low_pass_filter(converted_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) + aoflagger.visualize(converted_data, "Fit #" .. i, i - 1) - local tmp = original_data - data - tmp:set_mask(data) - data = tmp + local tmp = converted_copy - converted_data + tmp:set_mask(converted_data) + converted_data = tmp - aoflagger.visualize(data, "Residual #" .. i, i + iteration_count) + aoflagger.visualize(converted_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, + converted_data, + converted_copy, base_threshold, base_threshold * transient_threshold_factor, true, true ) else - aoflagger.sumthreshold(data, base_threshold, base_threshold * transient_threshold_factor, true, true) + aoflagger.sumthreshold(converted_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) + converted_data:join_mask(converted_copy) end -- Helper function used below @@ -126,14 +127,14 @@ function execute(input) if contains(inpPolarizations, polarization) then if input:is_complex() then - data = data:convert_to_complex("complex") + converted_data = converted_data:convert_to_complex("complex") end - input:set_polarization_data(polarization, data) + input:set_polarization_data(polarization, converted_data) else - input:join_mask(data) + input:join_mask(converted_data) end - aoflagger.visualize(data, "Residual #" .. iteration_count, 2 * iteration_count) + aoflagger.visualize(converted_data, "Residual #" .. iteration_count, 2 * iteration_count) aoflagger.set_progress(ipol, #flag_polarizations) end -- end of polarization iterations -- GitLab