diff --git a/CHANGELOG.md b/CHANGELOG.md index 0551ca1..2fdde6f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Added more detailed instructions for installing analysisUtils (#340). - Included initial documentation (#348). - Add badges to README.md (#358). +- Add more control over multiscale clean scales (#362). ### Changed diff --git a/docs/docs_phangsPipeline/keys/config_key.rst b/docs/docs_phangsPipeline/keys/config_key.rst index 723144b..4aef570 100644 --- a/docs/docs_phangsPipeline/keys/config_key.rst +++ b/docs/docs_phangsPipeline/keys/config_key.rst @@ -121,6 +121,17 @@ well as more complicated joint array combinations (e.g. ``12m+7m``). attempt to detect and assign data to arrays on its own. - ``clean_scales_arcsec``: These are the scales used for multiscale cleaning of data from this array. The units are arcseconds. +- ``clean_scales_beam``: These are the scales used for multiscale + cleaning of data from this array. The units are in terms of the synthesized beam. +- ``clean_scales_auto``: If True, then the pipeline will automatically + determine the scales to use for multiscale cleaning. This is determined as ``[0, 1*beam]``, + and then multiplying the largest scale by ``clean_scales_auto_factor`` until + ``clean_scales_max_las_fraction`` * largest angular scale in the measurement set is reached. + If False, then the user must provide ``clean_scales_arcsec`` and/or ``clean_scales_beam``. +- ``clean_scales_auto_factor``: Factor to multiply clean scales by, if automatically setting them. + Defaults to 2, and cannot be less than 1. +- ``clean_scales_max_las_fraction``: Maximum fraction of the LAS in the measurement set to consider + clean scales up to. Defaults to 1, and cannot be greater than 1. - ``requires``: By default, we require ALL arrays that make up the configuration, but this can be changed to an OR if you only need one of a certain combination, e.g. diff --git a/phangsPipeline/clean_call.py b/phangsPipeline/clean_call.py index 1f5138e..28251a6 100644 --- a/phangsPipeline/clean_call.py +++ b/phangsPipeline/clean_call.py @@ -1,12 +1,16 @@ """ This is a dummy CleanCall class for dry run only, or to be inheritted by casaImagingRoutines.CleanCall. """ - +import copy import logging +import os import re +import analysisUtils as au import numpy as np +from .casaStuff import imhead + logger = logging.getLogger(__name__) #region class CleanCall @@ -141,23 +145,135 @@ def set_reffreq_ghz(self, value=None): self.clean_params['reffreq']='' return() - def set_multiscale_arcsec(self, scales=[]): + def set_multiscale_clean_scales(self, + key_handler, + config: str, + imaging_method: str="tclean", + ): """ - Set the scales for deconvolution in acseconds. Requires that a + Set the scales for deconvolution in arcseconds. Requires that a cell size already be defined so that these can be translated into pixel units. + + Args: + key_handler (handlerTemplate.HandlerTemplate): KeyHandler object + config (str): Configuration name + imaging_method (str, optional): Should be one of 'tclean', 'sdintimaging'. + Defaults to 'tclean'. """ cell_in_pix = self.get_cell_in_arcsec() if cell_in_pix is None: return() + clean_scales_auto = key_handler.get_clean_scales_auto_for_config(config) + if clean_scales_auto: + + logger.info("Setting clean scales automatically") + + # Get parameters for the automatic clean scale handling + clean_scales_auto_factor = key_handler.get_clean_scales_auto_factor_for_config(config) + clean_scales_max_las_fraction = key_handler.get_clean_scales_max_las_fraction_for_config(config) + + # Calculate LAS and beam size + vis = self.get_param("vis") + las = au.estimateMRS(vis) + beam = self.estimate_synthesised_beam(imaging_method=imaging_method) + + # Start with scales of 0 and the beam + scales = [0, float(beam)] + + # Keep multiplying the scale by the set factor, until we reach the fraction of LAS + while scales[-1] < clean_scales_max_las_fraction * las: + scale = scales[-1] * clean_scales_auto_factor + scales.append(float(scale)) + + # Trim off the final scale, since we'll have overshot + scales = scales[:-1] + + logger.info("Will use scales (arcsec): "+str(scales)) + + else: + + logger.info("Setting clean scales manually:") + + clean_scales_arcsec = key_handler.get_clean_scales_arcsec_for_config(config) + if clean_scales_arcsec is None: + clean_scales_arcsec = [] + + logger.info("clean_scales_arcsec: "+str(clean_scales_arcsec)) + + clean_scales_beam = key_handler.get_clean_scales_beam_for_config(config) + if clean_scales_beam is None: + clean_scales_beam = [] + + logger.info("clean_scales_beam: " + str(clean_scales_beam)) + + if len(clean_scales_arcsec) + len(clean_scales_beam) == 0: + raise ValueError("At least one of clean_scales_arcsec, clean_scales_beam must be defined") + + # If we have beam scales defined, convert to arcsec + if len(clean_scales_beam) > 0: + beam = self.estimate_synthesised_beam(imaging_method=imaging_method) + + clean_scales_beam = np.asarray(clean_scales_beam) * beam + + # Combine these lists, make sure they're unique + scales = list(clean_scales_arcsec) + list(clean_scales_beam) + scales = [float(s) for s in np.unique(scales)] + scales_in_pix = [] for this_scale_arcsec in scales: this_scale_pix = this_scale_arcsec/cell_in_pix scales_in_pix.append(this_scale_pix) + # Sort, round, make sure unique scales_in_pix.sort() - self.clean_params['scales'] = [int(t) for t in scales_in_pix] + scales_in_pix = [np.round(t) for t in scales_in_pix] + scales_in_pix = [int(s) for s in np.unique(scales_in_pix)] + + self.clean_params['scales'] = copy.deepcopy(scales_in_pix) + + return True + + def estimate_synthesised_beam( + self, + imaging_method: str = "tclean", + ): + """Estimate the synthesised beam, either from an existing image or from the visibility data + + Args: + imaging_method (str, optional): either 'tclean' or 'sdintimaging'. Defaults to 'tclean'. + """ + + im_name = self.get_param("imagename") + + # Get out expected image name + if imaging_method == "tclean": + im = f"{im_name}.image" + elif imaging_method == "sdintimaging": + im = f"{im_name}.int.cube" + else: + raise ValueError("Unknown imaging method: " + str(imaging_method)) + + # If available, pull from the existing image + if os.path.exists(im): + + logger.info("Getting beam from existing image") + + # Get average of bmaj/bmin + bmaj = imhead(imagename=im, mode="get", hdkey="bmaj") + bmin = imhead(imagename=im, mode="get", hdkey="bmin") + beam = (bmaj["value"] + bmin["value"]) / 2 + + # Otherwise, estimate using visibilities + else: + + logger.info("Estimating beam from visibilities") + + vis = self.get_param("vis") + beam = au.estimateSynthesizedBeam(vis) + + return beam def set_round_uvtaper_arcsec(self, taper=0.0): """ diff --git a/phangsPipeline/handlerImaging.py b/phangsPipeline/handlerImaging.py index 721f58e..e0084b6 100644 --- a/phangsPipeline/handlerImaging.py +++ b/phangsPipeline/handlerImaging.py @@ -678,18 +678,25 @@ def task_assign_multiscales( self, clean_call=None, config=None, + imaging_method : str='tclean', ): """ Look up the appropriate scales for multi-scale clean associated with this config. + + Args: + imaging_method (str, optional): Either 'tclean' + or 'sdintimaging'. Defaults to 'tclean'. """ if clean_call is None: logger.warning("Require a clean_call object. Returning.") return () - scales_to_clean = self._kh.get_clean_scales_for_config(config=config) - clean_call.set_multiscale_arcsec(scales=scales_to_clean) + clean_call.set_multiscale_clean_scales(key_handler=self._kh, + config=config, + imaging_method=imaging_method, + ) return () @@ -1326,7 +1333,10 @@ def recipe_phangsalma_imaging( # Look up the angular scales to clean for this config - self.task_assign_multiscales(config=config, clean_call=clean_call) + self.task_assign_multiscales(config=config, + clean_call=clean_call, + imaging_method=imaging_method, + ) # Run a multiscale clean until it converges. diff --git a/phangsPipeline/handlerImagingChunked.py b/phangsPipeline/handlerImagingChunked.py index 5f4ad5b..2ab62ee 100644 --- a/phangsPipeline/handlerImagingChunked.py +++ b/phangsPipeline/handlerImagingChunked.py @@ -699,9 +699,11 @@ def task_initialize_clean_call( elif stage == 'multiscale': clean_call.set_param('deconvolver', 'multiscale') - # Set the multiscale to use: - scales_to_clean = self._kh.get_clean_scales_for_config(config=self.config) - clean_call.set_multiscale_arcsec(scales=scales_to_clean) + # Set the multiscale scales to use + clean_call.set_multiscale_clean_scales(key_handler=self._kh, + config=self.config, + imaging_method=self.imaging_method, + ) elif stage == 'singlescale': clean_call.set_param('deconvolver', 'hogbom') diff --git a/phangsPipeline/handlerKeys.py b/phangsPipeline/handlerKeys.py index 2a3cfe2..27821ba 100644 --- a/phangsPipeline/handlerKeys.py +++ b/phangsPipeline/handlerKeys.py @@ -880,7 +880,7 @@ def _read_all_keys(self): self._distance_dict = key_readers.batch_read( key_list=self._distance_keys, reader_function=key_readers.read_distance_key, key_dir=self._key_dir) - + self._window_dict = key_readers.batch_read( key_list=self._window_keys, reader_function=key_readers.read_window_key, key_dir=self._key_dir) @@ -1351,7 +1351,7 @@ def get_vfield_dir_for_target(self, target=None, changeto=False): changeto is true, then change directory to that location. """ return self._get_dir_for_target(target=target, changeto=changeto, vfield=True) - + def get_cleanmask_dir_for_target(self, target=None, changeto=False): """ Return the release working directory given a target name. If @@ -1739,7 +1739,7 @@ def get_distance_for_target(self, target=None): distance = self._distance_dict[target_name]['distance'] return distance - + def get_window_for_target(self, target=None): """ Get the velocity window (in km/s) associated with a target. If the @@ -2704,7 +2704,7 @@ def get_interf_config_for_feather_config( return feather_config_dict[feather_config]['interf_config'] - def get_clean_scales_for_config( + def get_clean_scales_arcsec_for_config( self, config=None, ): @@ -2716,12 +2716,92 @@ def get_clean_scales_for_config( if config is None: return None - if config in self._config_dict['interf_config'].keys(): - this_dict = self._config_dict['interf_config'][config] - else: + clean_scales_arcsec = self._config_dict.get("interf_config", {}).get(config, {}).get("clean_scales_arcsec", []) + + return clean_scales_arcsec + + def get_clean_scales_beam_for_config( + self, + config=None, + ): + """ + Return the angular scales as multiples of the beam + used for multiscale clean for an interferometric configuration. + """ + + if config is None: + return None + + clean_scales_beam = self._config_dict.get("interf_config", {}).get(config, {}).get("clean_scales_beam", []) + + return clean_scales_beam + + def get_clean_scales_auto_for_config( + self, + config=None, + ): + """ + Return whether we are automatically setting clean + scales for multiscale clean for an interferometric configuration. + """ + + if config is None: + return False + + clean_scales_auto = ( + self._config_dict.get("interf_config", {}) + .get(config, {}) + .get("clean_scales_auto", False) + ) + + return clean_scales_auto + + def get_clean_scales_auto_factor_for_config( + self, + config=None, + ): + """ + Return the multiplicative factor for automatic clean scale + """ + + if config is None: return None - return this_dict['clean_scales_arcsec'] + clean_scales_auto_factor = ( + self._config_dict.get("interf_config", {}) + .get(config, {}) + .get("clean_scales_auto_factor", 2) + ) + + # If we have a value less than 1, then we won't converge + if clean_scales_auto_factor <= 1: + logger.warning("clean_scales_auto_factor should not be smaller than 1. Will set to 1.1") + clean_scales_auto_factor = 1.1 + + return clean_scales_auto_factor + + def get_clean_scales_max_las_fraction_for_config( + self, + config=None, + ): + """ + Return the maximum fraction of the LAS for automatic clean scale + """ + + if config is None: + return None + + clean_scales_max_las_fraction = ( + self._config_dict.get("interf_config", {}) + .get(config, {}) + .get("clean_scales_max_las_fraction", 1) + ) + + if clean_scales_max_las_fraction > 1: + logger.warning("clean_scales_max_las_fraction should not be larger than 1. Will set to 1") + clean_scales_max_las_fraction = 1 + + return clean_scales_max_las_fraction def get_ang_res_dict(self, config=None, product=None, ): @@ -2797,7 +2877,7 @@ def get_derived_kwargs(self, config=None, product=None, if product not in self._derived_dict[config].keys(): return {} - + if kwarg_type not in self._derived_dict[config][product].keys(): return {} @@ -2901,10 +2981,33 @@ def print_configs( logger.info("... " + this_config) this_arrays = self._config_dict['interf_config'][this_config]['array_tags'] this_other_config = self._config_dict['interf_config'][this_config]['feather_config'] - scales_for_clean = self._config_dict['interf_config'][this_config]['clean_scales_arcsec'] + + # Get out various clean scales + scales_for_clean_arcsec = self.get_clean_scales_arcsec_for_config(this_config) + if scales_for_clean_arcsec is None: + scales_for_clean_arcsec = [] + + scales_for_clean_beam = self.get_clean_scales_beam_for_config(this_config) + if scales_for_clean_beam is None: + scales_for_clean_beam = [] + + scales_for_clean_auto = self.get_clean_scales_auto_for_config(this_config) + logger.info("... ... includes arrays " + str(this_arrays)) logger.info("... ... maps to feather config " + str(this_other_config)) - logger.info("... ... clean these scales in arcsec " + str(scales_for_clean)) + + if scales_for_clean_auto: + logger.info("... ... automatically set clean scales") + else: + + # Crash out if we don't have any clean scales defined + if len(scales_for_clean_beam) + len(scales_for_clean_arcsec) == 0: + raise ValueError("At least one of clean_scales_arcsec, clean_scales_beam must be defined") + + if len(scales_for_clean_arcsec) > 0: + logger.info("... ... clean these scales in arcsec: " + str(scales_for_clean_arcsec)) + if len(scales_for_clean_beam) > 0: + logger.info("... ... clean these scales as multiples of the beam: " + str(scales_for_clean_beam)) if 'feather_config' in self._config_dict: logger.info("Feather Configurations") diff --git a/phangsPipeline/utilsKeyReaders.py b/phangsPipeline/utilsKeyReaders.py index 52ed44f..c7fe71b 100644 --- a/phangsPipeline/utilsKeyReaders.py +++ b/phangsPipeline/utilsKeyReaders.py @@ -740,7 +740,11 @@ def read_config_key(fname='', existing_dict=None, delim=None): 'res_max_pc': 0.0, 'res_step_factor': 1.0, 'res_list': [], - 'clean_scales_arcsec': [] + 'clean_scales_arcsec': [], + 'clean_scales_beam': [], + 'clean_scales_auto': False, + 'clean_scales_auto_factor': 2.0, + 'clean_scales_max_las_fraction': 1.0, } if this_type == "feather_config":