diff --git a/bin/rapthor b/bin/rapthor
index c8e87b4e725c6f2a086a798a855af7cf36174d1f..c574f8166c3c134cc37f30136aea5039a86f61de 100755
--- a/bin/rapthor
+++ b/bin/rapthor
@@ -22,8 +22,7 @@ import sys
if __name__ == '__main__':
- parser = optparse.OptionParser(usage='%prog <parset>',
- version='%%prog v%s' % (version))
+ parser = optparse.OptionParser(usage='%prog <parset>', version='%%prog v%s' % (version))
parser.add_option('-q', help='enable quiet mode', action='store_true', default=False)
parser.add_option('-v', help='enable verbose mode', action='store_true', default=False)
parser.add_option('-r', help='reset one or more operations', action='store_true', default=False)
diff --git a/docs/source/operations.rst b/docs/source/operations.rst
index 8c29809b595b562e0dca2b217e6c72116c15fe23..e2104e6124986812237deeecdfb2b64da230ac63 100644
--- a/docs/source/operations.rst
+++ b/docs/source/operations.rst
@@ -11,7 +11,7 @@ Most of the processing performed by Rapthor is done in "operations," which are s
Calibrate
---------
-This operation calibrates the data using the current sky model. The exact steps done during calibration depend on the strategy, but essentially there are two main parts: a phase-only (scalar) solve on short timescales (the "fast-phase" solve, which corrects for ionospheric errors) and a phase and amplitude (diagonal) solve on long time scales (the "slow-gain" solve, which corrects for beam errors). The slow-gain solve is divided into two parts: an optional amplitude-only solve that constrains all stations to have the same solutions and a phase plus amplitude solve (without the station constraint and usually with a longer solution interval). This calibration strategy is based on the LBA strategy of the `LiLF pipeline <https://github.com/revoltek/LiLF>`_, with the idea that the same strategy can be used for both HBA and LBA. This is similar to the way the calibrator pipeline works in `LINC <https://linc.readthedocs.io/>`_, which has been inspired by LiLF. Lastly, processing of the resulting solutions is done, including smoothing, renormalization, and the generation of a-term images (if 2-D screens are used).
+This operation calibrates the data using the current sky model. The exact steps done during calibration depend on the strategy, but essentially there are two main parts: a phase-only (scalar) solve on short timescales (the "fast-phase" solve, which corrects for ionospheric errors) and a phase and amplitude (diagonal) solve on long time scales (the "slow-gain" solve, which corrects for beam errors). The slow-gain solve is divided into two parts: an optional amplitude-only solve that constrains all stations to have the same solutions and a phase plus amplitude solve (without the station constraint and usually with a longer solution interval). This calibration strategy is based on the LBA strategy of the `LiLF pipeline <https://github.com/revoltek/LiLF>`_, with the idea that the same strategy can be used for both HBA and LBA. This is similar to the way the calibrator pipeline works in `LINC <https://linc.readthedocs.io/>`_, which has been inspired by LiLF. Lastly, processing of the resulting solutions is done, including smoothing and renormalization.
For calibration, Rapthor searches for bright, compact sources (or groups of sources) throughout the field to use as calibrator sources. A target (apparent) flux density is used to ensure that the calibrators are sufficiently bright (set by :term:`target_flux` in the processing strategy). Rapthor then tessellates the full sky model, using the calibrators as the facet centers. This method ensures that each calibration patch (or facet) has a bright calibrator source in it. Despite this designation of calibrators for the tesselation, all sources are used in the calibration (not just the bright sources).
@@ -133,8 +133,8 @@ Primary products:
* In ``visibilities/image_X/sector_Y``, where ``X`` is the cycle number and ``Y`` is the image sector number (only if the :term:`save_visibilities` parameter is set to ``True``):
* ``*.ms`` - measurement sets used as input to WSClean for imaging. Depending on
the value of :term:`dde_method`, some or all of the calibration solutions may be
- preapplied: a value of "none" will preapply all solutions, whereas values of
- "screens" or "facets" will preapply only the full-Jones solutions (if
+ preapplied: a value of "single" will preapply all solutions, whereas a value of
+ "full" will preapply only the full-Jones solutions (if
available), since the direction-dependent solutions in those cases are applied
by WSClean itself. These MS files can be useful for further imaging or self
calibration outside of Rapthor.
diff --git a/docs/source/parset.rst b/docs/source/parset.rst
index e649e9efbdb3be1ad99811b6cf5b53da92c84809..a7951cd6038fcc45f6b5e35bd7c840d27367a097 100644
--- a/docs/source/parset.rst
+++ b/docs/source/parset.rst
@@ -149,6 +149,19 @@ The available options are described below under their respective sections.
calibration patches, in which case the empty patches are removed and the
layout of the remaining patches is set using Voronoi tessellation.
+ dde_mode
+ Mode to use to derive and correct for direction-dependent effects: ``faceting`` or
+ ``hybrid`` (default = ``faceting``). If ``faceting``, Voronoi faceting is used
+ throughout the processing. If ``hybrid``, faceting is used only during the self
+ calibration steps; in the final cycle (done after self calibration has been
+ completed successfully), IDGCal is used during calibration to generate smooth 2-D
+ screens that are then applied by WSClean in the final imaging step.
+
+ .. note::
+
+ The hybrid mode is not yet available; it will be enabled in a future
+ update.
+
.. _parset_calibration_options:
``[calibration]``
@@ -353,21 +366,18 @@ The available options are described below under their respective sections.
Use multiscale cleaning (default = ``True``)?
dde_method
- Method to use to correct for direction-dependent effects during imaging: ``none``,
- ``facets``, or ``screens`` (default = ``facets``). If ``none``, the solutions
- closest to the image centers will be used. If ``facets``, Voronoi faceting is
- used. If ``screens``, smooth 2-D screens are used.
-
- screen_type
- Type of screen to use (default = ``tessellated``), if :term:`dde_method` =
- ``screens``: ``tessellated`` (simple, smoothed Voronoi tessellated screens) or
- ``kl`` (Karhunen-Lo`eve screens).
+ Method to use to correct for direction-dependent effects during imaging:
+ ``single`` or ``full`` (default = ``full``). If ``single``, a single,
+ direction-independent solution (i.e., constant across the image sector) will be
+ applied for each sector. In this case, the solution applied is the one in the
+ direction closest to each sector center. If ``full``, the full,
+ direction-dependent solutions are applied (using either facets or screens).
save_visibilities
Save visibilities used for imaging (default = ``False``). If ``True``, the imaging
MS files will be saved, with the the direction-independent full-Jones solutions,
if available, applied. Note, however, that the direction-dependent solutions will
- not be applied unless :term:`dde_method` = ``none``, in which case the solutions
+ not be applied unless :term:`dde_method` = ``single``, in which case the solutions
closest to the image centers are used.
idg_mode
diff --git a/docs/source/strategy.rst b/docs/source/strategy.rst
index 88298c315dc37394668a73822099da6c55c91501..d51954c0c2af99e9f760261c23ca604365622768 100644
--- a/docs/source/strategy.rst
+++ b/docs/source/strategy.rst
@@ -113,7 +113,7 @@ The following processing parameters can be set for each cycle:
Boolean flag that determines whether the outlier sources (sources that lie outside of any imaging sector region) should be peeled for this cycle. Outliers can only be peeled once (unlike bright sources, see below), as they are not added back for subsequent selfcal cycles. Note that, because they are not imaged, outlier source models do not change during self calibration: however, the solutions they receive may change. To include one or more outlier sources in self calibration, a small imaging sector can be placed on each outlier of interest. The outliers will than be imaging and its model updated with the rest of the field.
peel_bright_sources
- Boolean flag that determines whether the bright sources should be peeled for this cycle (for imaging only). The peeled bright sources are added back before subsequent selfcal cycles are performed (so they are included in the calibration, etc.). Generally, peeling of bright sources is beneficial when using screens but not when using facets.
+ Boolean flag that determines whether the bright sources should be peeled for this cycle (for imaging only). The peeled bright sources are added back before subsequent selfcal cycles are performed (so they are included in the calibration, etc.). Currently, peeling is not supported when screens are used.
max_normalization_delta
Float that sets the maximum allowed fractional delta from unity for the per-station normalization.
diff --git a/docs/source/tips.rst b/docs/source/tips.rst
index 48bc0e9abf21f40772bcc176534782f0973935a5..a04b40ab58ded6085b89cf54e7090e1c25d5b32c 100644
--- a/docs/source/tips.rst
+++ b/docs/source/tips.rst
@@ -76,10 +76,11 @@ Creating a dataset for further self calibration
* Set the :term:`save_visibilities` parameter to ``True`` to save the MS files
used for the imaging. The output MS files will be located in
- ``dir_working/visibilities``. Furthermore, :term:`dde_method` can be set to
- "none" to preapply the solutions in the direction closest to the sector center
- (other values of :term:`dde_method` will not preapply the solutions, since the
- solutions in those cases are applied by WSClean during imaging).
+ ``dir_working/visibilities``. Furthermore, :term:`dde_method` should be set to
+ "single" to apply the solutions in single direction (the direction closest to
+ the sector center). Other values of :term:`dde_method` will not apply the
+ solutions to the visibilities, since the solutions in those cases are applied by
+ WSClean during imaging.
* Define the imaging sectors to cover the targets of interest. Multiple sectors can
be used, and a set of calibrated visibilities will be generated for each sector.
diff --git a/pyproject.toml b/pyproject.toml
index 5c6fb9ebf4cc7d29737975343ee5ce68e3e585d5..61cf85c37608d107d46d45554c70e610c8d51355 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -72,7 +72,6 @@ script-files = [
"rapthor/scripts/combine_h5parms.py",
"rapthor/scripts/concat_ms.py",
"rapthor/scripts/filter_skymodel.py",
- "rapthor/scripts/make_aterm_images.py",
"rapthor/scripts/make_catalog_from_image_cube.py",
"rapthor/scripts/make_image_cube.py",
"rapthor/scripts/make_mosaic.py",
@@ -82,7 +81,6 @@ script-files = [
"rapthor/scripts/normalize_flux_scale.py",
"rapthor/scripts/process_gains.py",
"rapthor/scripts/regrid_image.py",
- "rapthor/scripts/split_h5parms.py",
"rapthor/scripts/subtract_sector_models.py",
]
diff --git a/rapthor/lib/field.py b/rapthor/lib/field.py
index dfb13e0342ee480e5eba47366ab667252c556c0d..5aa92ae567e0b07e3ac7d3c7e5909cd817a034e1 100644
--- a/rapthor/lib/field.py
+++ b/rapthor/lib/field.py
@@ -54,6 +54,7 @@ class Field(object):
self.dd_interval_factor = self.parset['calibration_specific']['dd_interval_factor']
self.h5parm_filename = self.parset['input_h5parm']
self.fulljones_h5parm_filename = self.parset['input_fulljones_h5parm']
+ self.dde_mode = self.parset['dde_mode']
self.fast_smoothnessconstraint = self.parset['calibration_specific']['fast_smoothnessconstraint']
self.fast_smoothnessreffrequency = self.parset['calibration_specific']['fast_smoothnessreffrequency']
self.fast_smoothnessrefdistance = self.parset['calibration_specific']['fast_smoothnessrefdistance']
@@ -70,11 +71,6 @@ class Field(object):
self.stepsigma = self.parset['calibration_specific']['stepsigma']
self.tolerance = self.parset['calibration_specific']['tolerance']
self.dde_method = self.parset['imaging_specific']['dde_method']
- if self.dde_method == 'screens':
- self.use_screens = True
- else:
- self.use_screens = False
- self.screen_type = self.parset['imaging_specific']['screen_type']
self.save_visibilities = self.parset['imaging_specific']['save_visibilities']
self.use_mpi = self.parset['imaging_specific']['use_mpi']
self.parallelbaselines = self.parset['calibration_specific']['parallelbaselines']
@@ -728,7 +724,7 @@ class Field(object):
self.target_flux = target_flux
def update_skymodels(self, index, regroup, target_flux=None, target_number=None,
- calibrator_max_dist_deg=None, final=False):
+ calibrator_max_dist_deg=None, combine_current_and_intial=False):
"""
Updates the source and calibration sky models from the output sector sky model(s)
@@ -749,8 +745,9 @@ class Field(object):
Target number of patches for grouping
calibrator_max_dist_deg : float, optional
Maximum distance in degrees from phase center for grouping
- final : bool, optional
- If True, process as the final pass (combine initial and new sky models)
+ combine_current_and_intial : bool, optional
+ If True, combine the initial and current sky models (needed for the final
+ calibration in order to include potential outlier sources)
"""
# Except for the first iteration, use the results of the previous iteration to
# update the sky models, etc.
@@ -850,11 +847,11 @@ class Field(object):
skymodel_apparent_sky = None
# Concatenate the starting sky model with the new one. This step needs to be
- # done if this iteration is set as a final pass (so that sources outside of
- # imaged areas can be subtracted, since we have to go back to the original
- # input MS files for which no subtraction has been done) or if all sources
- # (and not only the imaged sources) are to be used in calibration
- if final or not self.imaged_sources_only:
+ # done if, e.g., this is the final cycle and sources outside of imaged areas
+ # must be subtracted, since we have to go back to the original input MS files
+ # for which no subtraction has been done) or if all sources (and not only the
+ # imaged sources) are to be used in calibration
+ if combine_current_and_intial or not self.imaged_sources_only:
# Load starting sky model and regroup to one patch per entry to ensure
# any existing patches are removed (otherwise they may propagate to
# the DDE direction determination, leading to unexpected results)
@@ -892,7 +889,7 @@ class Field(object):
# Plot an overview of the field for this cycle, showing the calibration facets
# (patches)
self.log.info('Plotting field overview with calibration patches...')
- if index == 1 or final:
+ if index == 1 or combine_current_and_intial:
# Check the sky model bounds, as they may differ from the sector ones
check_skymodel_bounds = True
else:
@@ -1076,9 +1073,7 @@ class Field(object):
# Compute bounding box for all imaging sectors and store as a
# a semi-colon-separated list of [maxRA; minDec; minRA; maxDec] (we use semi-
# colons as otherwise the workflow parset parser will split the list). Also
- # store the midpoint as [midRA; midDec]. These values are needed for the aterm
- # image generation, so we use the padded polygons to ensure that the final
- # bounding box encloses all of the images *with* padding included.
+ # store the midpoint as [midRA; midDec].
# Note: this is just once, rather than each time the sector borders are
# adjusted, so that the image sizes do not change with iteration (so
# mask images from previous iterations may be used)
@@ -1610,7 +1605,7 @@ class Field(object):
self.update_skymodels(index, step_dict['regroup_model'],
target_flux=target_flux, target_number=target_number,
calibrator_max_dist_deg=calibrator_max_dist_deg,
- final=final)
+ combine_current_and_intial=final)
self.remove_skymodels() # clean up sky models to reduce memory usage
# Always try to peel non-calibrator sources if LBA (we check below whether
diff --git a/rapthor/lib/observation.py b/rapthor/lib/observation.py
index 81aff8119b99244de0e07bcef68702a1ce84c78d..e70c2777a9491115991a8a74098ad492bcca3c51 100644
--- a/rapthor/lib/observation.py
+++ b/rapthor/lib/observation.py
@@ -398,13 +398,6 @@ class Observation(object):
else:
self.parameters[f'solutions_per_direction_{solve_type}'] = [[1] * len(calibrator_fluxes)] * nchunks
- # Set the number of segments to split the h5parm files into for screen fitting.
- # Try to split so that each file gets at least two solutions
- solint_slow_timestep = max(self.parameters['solint_slow_timestep_joint'][0],
- self.parameters['solint_slow_timestep_separate'][0])
- self.parameters['nsplit_fast'] = [max(1, int(self.numsamples / solint_fast_timestep / 2))]
- self.parameters['nsplit_slow'] = [max(1, int(self.numsamples / solint_slow_timestep / 2))]
-
# Set the smoothnessreffrequency for the fast solves, if not set by the user
fast_smoothnessreffrequency = parset['calibration_specific']['fast_smoothnessreffrequency']
if fast_smoothnessreffrequency is None:
@@ -460,7 +453,7 @@ class Observation(object):
def set_imaging_parameters(self, sector_name, cellsize_arcsec, max_peak_smearing, width_ra,
width_dec, solve_fast_timestep, solve_slow_freqstep,
- use_screens):
+ apply_screens):
"""
Sets the imaging parameters
@@ -478,7 +471,7 @@ class Observation(object):
Solution interval in sec for fast solve
solve_slow_freqstep : float
Solution interval in Hz for slow solve
- use_screens : bool
+ apply_screens : bool
If True, use setup appropriate for screens
"""
mean_freq_mhz = self.referencefreq / 1e6
@@ -512,7 +505,7 @@ class Observation(object):
target_timewidth_sec = min(max_timewidth_sec, self.get_target_timewidth(delta_theta_deg,
resolution_deg, peak_smearing_rapthor))
- if use_screens:
+ if apply_screens:
# Ensure we don't average more than the solve time step, as we want to
# preserve the time resolution that matches that of the screens
target_timewidth_sec = min(target_timewidth_sec, solve_fast_timestep)
diff --git a/rapthor/lib/parset.py b/rapthor/lib/parset.py
index 005d7c1e6b74665e348f4cd50135d4e3b23f2ce2..2a7752d6d289b8c7745eb7bbd2a23aa211d4dfce 100644
--- a/rapthor/lib/parset.py
+++ b/rapthor/lib/parset.py
@@ -282,8 +282,7 @@ class Parset:
for opt, valid_values in {
"idg_mode": ("cpu", "gpu", "hybrid"),
- "dde_method": ("none", "screens", "facets"),
- "screen_type": ("kl", "tesselated"),
+ "dde_method": ("single", "full"),
"pol_combine_method": ("link", "join"),
}.items():
if options[opt] not in valid_values:
diff --git a/rapthor/lib/screen.py b/rapthor/lib/screen.py
deleted file mode 100644
index 13f01931b349f17d048c0438f1f5a6cdc57966d0..0000000000000000000000000000000000000000
--- a/rapthor/lib/screen.py
+++ /dev/null
@@ -1,827 +0,0 @@
-"""
-Module that holds screen-related classes and functions
-"""
-import logging
-import numpy as np
-import os
-from rapthor.lib import miscellaneous as misc
-from rapthor.lib import cluster
-from astropy.coordinates import Angle
-from losoto.h5parm import h5parm
-import scipy.interpolate as si
-from astropy.io import fits as pyfits
-import scipy.ndimage as ndimage
-from losoto.operations import stationscreen
-from astropy import wcs
-import lsmtool
-from shapely.geometry import Point
-from scipy.spatial import Voronoi
-import shapely.geometry
-import shapely.ops
-import multiprocessing
-from multiprocessing import Pool, RawArray
-import itertools
-
-
-class Screen(object):
- """
- Master class for a-term screens
-
- Parameters
- ----------
- name : str
- Name of screen
- h5parm_filename : str
- Filename of h5parm containing the input solutions
- skymodel_filename : str
- Filename of input sky model
- ra : float
- RA in degrees of screen center
- dec : float
- Dec in degrees of screen center
- width_ra : float
- Width of screen in RA in degrees, corrected to Dec = 0
- width_dec : float
- Width of screen in Dec in degrees
- solset_name: str, optional
- Name of solset of the input h5parm to use
- phase_soltab_name: str, optional
- Name of the phase soltab of the input h5parm to use
- amplitude_soltab_name: str, optional
- Name of amplitude soltab of the input h5parm to use
- """
- def __init__(self, name, h5parm_filename, skymodel_filename, ra, dec, width_ra, width_dec,
- solset_name='sol000', phase_soltab_name='phase000', amplitude_soltab_name=None):
- self.name = name
- self.log = logging.getLogger('rapthor:{}'.format(self.name))
- self.input_h5parm_filename = h5parm_filename
- self.input_skymodel_filename = skymodel_filename
- self.input_solset_name = solset_name
- self.input_phase_soltab_name = phase_soltab_name
- self.input_amplitude_soltab_name = amplitude_soltab_name
- if self.input_amplitude_soltab_name is not None:
- self.phase_only = False
- else:
- self.phase_only = True
- if type(ra) is str:
- ra = Angle(ra).to('deg').value
- if type(dec) is str:
- dec = Angle(dec).to('deg').value
- self.ra = ra
- self.dec = dec
- width = max(width_ra, width_dec) # force square image until rectangular ones are supported by IDG
- self.width_ra = width
- self.width_dec = width
- self.log_amps = False # sets whether amplitudes are log10 values or not
-
- def fit(self):
- """
- Fits screens to the input solutions
-
- This method should be defined in the subclasses
- """
- pass
-
- def interpolate(self, interp_kind='nearest'):
- """
- Interpolate the slow amplitude values to the fast-phase time and frequency grid
-
- Parameters
- ----------
- interp_kind : str, optional
- Kind of interpolation to use
- """
- if self.phase_only:
- return
-
- if len(self.times_amp) == 1:
- # If only a single time, we just repeat the values as needed
- new_shape = list(self.vals_amp.shape)
- new_shape[0] = self.vals_ph.shape[0]
- new_shape[1] = self.vals_ph.shape[1]
- self.vals_amp = np.resize(self.vals_amp, new_shape)
- else:
- # Interpolate amplitudes (in log space)
- if not self.log_amps:
- logvals = np.log10(self.vals_amp)
- else:
- logvals = self.vals_amp
- if self.vals_amp.shape[0] != self.vals_ph.shape[0]:
- f = si.interp1d(self.times_amp, logvals, axis=0, kind=interp_kind, fill_value='extrapolate')
- logvals = f(self.times_ph)
- if self.vals_amp.shape[1] != self.vals_ph.shape[1]:
- f = si.interp1d(self.freqs_amp, logvals, axis=1, kind=interp_kind, fill_value='extrapolate')
- logvals = f(self.freqs_ph)
- if not self.log_amps:
- self.vals_amp = 10**(logvals)
- else:
- self.vals_amp = logvals
-
- def make_fits_file(self, outfile, cellsize_deg, t_start_index,
- t_stop_index, aterm_type='gain'):
- """
- Makes a FITS data cube and returns the Header Data Unit
-
- Parameters
- ----------
- outfile : str
- Filename of output FITS file
- cellsize_deg : float
- Pixel size of image in degrees
- t_start_index : int
- Index of first time
- t_stop_index : int
- Index of last time
- aterm_type : str, optional
- Type of a-term solutions
- """
- ximsize = int(np.ceil(self.width_ra / cellsize_deg)) # pix
- yimsize = int(np.ceil(self.width_dec / cellsize_deg)) # pix
- misc.make_template_image(outfile, self.ra, self.dec, ximsize=ximsize,
- yimsize=yimsize, cellsize_deg=cellsize_deg, freqs=self.freqs_ph,
- times=self.times_ph[t_start_index:t_stop_index],
- antennas=self.station_names, aterm_type=aterm_type)
- hdu = pyfits.open(outfile, memmap=False)
- return hdu
-
- def make_matrix(self, t_start_index, t_stop_index, freq_ind, stat_ind, cellsize_deg,
- out_dir, ncpu):
- """
- Makes the matrix of values for the given time, frequency, and station indices
-
- This method should be defined in the subclasses, but should conform to the inputs
- below.
-
- Parameters
- ----------
- t_start_index : int
- Index of first time
- t_stop_index : int
- Index of last time
- t_start_index : int
- Index of frequency
- t_stop_index : int
- Index of station
- cellsize_deg : float
- Size of one pixel in degrees
- out_dir : str
- Full path to the output directory (needed for template file generation)
- ncpu : int, optional
- Number of CPUs to use (0 means all)
- """
- pass
-
- def get_memory_usage(self, cellsize_deg):
- """
- Returns memory usage per time slot in GB
-
- This method should be defined in the subclasses, but should conform to the inputs
- below.
-
- Parameters
- ----------
- cellsize_deg : float
- Size of one pixel in degrees
- """
- pass
-
- def write(self, out_dir, cellsize_deg, smooth_pix=0, interp_kind='nearest', ncpu=0):
- """
- Write the a-term screens to a FITS data cube
-
- Parameters
- ----------
- out_dir : str
- Output directory
- cellsize_deg : float
- Size of one pixel in degrees
- smooth_pix : int, optional
- Size of Gaussian in pixels to smooth with
- interp_kind : str, optional
- Kind of interpolation to use
- ncpu : int, optional
- Number of CPUs to use (0 means all)
- """
- self.ncpu = ncpu
-
- # Identify any gaps in time (frequency gaps are not allowed), as we need to
- # output a separate FITS file for each time chunk
- if len(self.times_ph) > 2:
- delta_times = self.times_ph[1:] - self.times_ph[:-1] # time at center of solution interval
- timewidth = np.min(delta_times)
- gaps = np.where(delta_times > timewidth*1.2)
- gaps_ind = gaps[0] + 1
- gaps_ind = np.append(gaps_ind, np.array([len(self.times_ph)]))
- else:
- gaps_ind = np.array([len(self.times_ph)])
-
- # Add additional breaks to gaps_ind to keep memory usage within that available
- if len(self.times_ph) > 2:
- available_mem_gb = cluster.get_available_memory()
- memory_usage = max(1, self.get_memory_usage(cellsize_deg))
- max_ntimes = max(1, int(available_mem_gb / memory_usage))
- check_gaps = True
- while check_gaps:
- check_gaps = False
- g_start = 0
- gaps_ind_copy = gaps_ind.copy()
- for gnum, g_stop in enumerate(gaps_ind_copy):
- if g_stop - g_start > max_ntimes:
- new_gap = g_start + int((g_stop - g_start) / 2)
- gaps_ind = np.insert(gaps_ind, gnum, np.array([new_gap]))
- check_gaps = True
- break
- g_start = g_stop
-
- # Input data are [time, freq, ant, dir, pol] for slow amplitudes
- # and [time, freq, ant, dir] for fast phases (scalarphase).
- # Output data are [RA, DEC, MATRIX, ANTENNA, FREQ, TIME].T
- # Loop over stations, frequencies, and times and fill in the correct
- # matrix values (matrix dimension has 4 elements: real XX, imaginary XX,
- # real YY and imaginary YY)
- outroot = self.name
- outfiles = []
- g_start = 0
- for gnum, g_stop in enumerate(gaps_ind):
- ntimes = g_stop - g_start
- outfile = os.path.join(out_dir, '{0}_{1}.fits'.format(outroot, gnum))
- hdu = self.make_fits_file(outfile, cellsize_deg, g_start, g_stop, aterm_type='gain')
- data = hdu[0].data
- for f, freq in enumerate(self.freqs_ph):
- for s, stat in enumerate(self.station_names):
- data[:, f, s, :, :, :] = self.make_matrix(g_start, g_stop, f, s,
- cellsize_deg, out_dir, self.ncpu)
-
- # Smooth if desired
- if smooth_pix > 0:
- for t in range(ntimes):
- data[t, f, s, :, :, :] = ndimage.gaussian_filter(data[t, f, s, :, :, :],
- sigma=(0, smooth_pix,
- smooth_pix),
- order=0)
-
- # Ensure there are no NaNs in the images, as WSClean will produced uncorrected,
- # uncleaned images if so. We replace NaNs with 1.0 and 0.0 for real and
- # imaginary parts, respectively
- # Note: we iterate over time to reduce memory usage
- for t in range(ntimes):
- for p in range(4):
- if p % 2:
- # Imaginary elements
- nanval = 0.0
- else:
- # Real elements
- nanval = 1.0
- data[t, :, :, p, :, :][~np.isfinite(data[t, :, :, p, :, :])] = nanval
-
- # Write FITS file
- hdu[0].data = data
- hdu.writeto(outfile, overwrite=True)
- outfiles.append(outfile)
- hdu = None
- data = None
-
- # Update start time index before starting next loop
- g_start = g_stop
-
- # Write list of filenames to a text file for later use
- list_file = open(os.path.join(out_dir, '{0}.txt'.format(outroot)), 'w')
- list_file.writelines([o+'\n' for o in outfiles])
- list_file.close()
-
- def process(self, ncpu=0):
- """
- Makes a-term images
-
- Parameters
- ----------
- ncpu : int, optional
- Number of CPUs to use (0 means all)
- """
- self.ncpu = ncpu
-
- # Fit screens to input solutions
- self.fit()
-
- # Interpolate best-fit parameters to common time and frequency grid
- self.interpolate()
-
-
-class KLScreen(Screen):
- """
- Class for KL (Karhunen-Lo`eve) screens
- """
- def __init__(self, name, h5parm_filename, skymodel_filename, ra, dec, width_ra, width_dec,
- solset_name='sol000', phase_soltab_name='phase000', amplitude_soltab_name=None):
- super(KLScreen, self).__init__(name, h5parm_filename, skymodel_filename, ra, dec, width_ra, width_dec,
- solset_name=solset_name, phase_soltab_name=phase_soltab_name,
- amplitude_soltab_name=amplitude_soltab_name)
-
- def fit(self):
- """
- Fits screens to the input solutions
- """
- # Open solution tables
- H = h5parm(self.input_h5parm_filename, readonly=False)
- solset = H.getSolset(self.input_solset_name)
- soltab_ph = solset.getSoltab(self.input_phase_soltab_name)
- if not self.phase_only:
- soltab_amp = solset.getSoltab(self.input_amplitude_soltab_name)
-
- # Set the position of the calibration patches to those of
- # the input sky model, as the patch positions written to the h5parm
- # file by DPPP may be different
- skymod = lsmtool.load(self.input_skymodel_filename)
- source_dict = skymod.getPatchPositions()
- source_positions = []
- for source in soltab_ph.dir:
- radecpos = source_dict[source.strip('[]')]
- source_positions.append([radecpos[0].value, radecpos[1].value])
- source_positions = np.array(source_positions)
- ra_deg = source_positions.T[0]
- dec_deg = source_positions.T[1]
- sourceTable = solset.obj._f_get_child('source')
- vals = [[ra*np.pi/180.0, dec*np.pi/180.0] for ra, dec in zip(ra_deg, dec_deg)]
- sourceTable = list(zip(*(soltab_ph.dir, vals)))
-
- # Now call LoSoTo's stationscreen operation to do the fitting. For the phase
- # screens, we reference the phases to the station with the least amount
- # of flagged solutions, drawn from the first 10 stations (to ensure it is
- # fairly central)
- ref_ind = misc.get_reference_station(soltab_ph, 10)
- adjust_order_amp = True
- screen_order_amp = min(12, max(3, int(np.round(len(source_positions) / 2))))
- adjust_order_ph = True
- screen_order = min(20, len(source_positions)-1)
- misc.remove_soltabs(solset, 'phase_screen000')
- misc.remove_soltabs(solset, 'phase_screen000resid')
- stationscreen.run(soltab_ph, 'phase_screen000', order=screen_order, refAnt=ref_ind,
- scale_order=True, adjust_order=adjust_order_ph, ncpu=self.ncpu)
- soltab_ph_screen = solset.getSoltab('phase_screen000')
- if not self.phase_only:
- misc.remove_soltabs(solset, 'amplitude_screen000')
- misc.remove_soltabs(solset, 'amplitude_screen000resid')
- stationscreen.run(soltab_amp, 'amplitude_screen000', order=screen_order_amp,
- niter=3, scale_order=False, adjust_order=adjust_order_amp,
- ncpu=self.ncpu)
- soltab_amp_screen = solset.getSoltab('amplitude_screen000')
- else:
- soltab_amp_screen = None
-
- # Read in the screen solutions and parameters
- self.vals_ph = soltab_ph_screen.val
- self.times_ph = soltab_ph_screen.time
- self.freqs_ph = soltab_ph_screen.freq
- if not self.phase_only:
- self.log_amps = True
- self.vals_amp = soltab_amp_screen.val
- self.times_amp = soltab_amp_screen.time
- self.freqs_amp = soltab_amp_screen.freq
- self.source_names = soltab_ph_screen.dir
- self.source_dict = solset.getSou()
- self.source_positions = []
- for source in self.source_names:
- self.source_positions.append(self.source_dict[source])
- self.station_names = soltab_ph_screen.ant
- self.station_dict = solset.getAnt()
- self.station_positions = []
- for station in self.station_names:
- self.station_positions.append(self.station_dict[station])
- self.height = soltab_ph_screen.obj._v_attrs['height']
- self.beta_val = soltab_ph_screen.obj._v_attrs['beta']
- self.r_0 = soltab_ph_screen.obj._v_attrs['r_0']
- self.pp = np.array(soltab_ph_screen.obj.piercepoint)
- self.midRA = soltab_ph_screen.obj._v_attrs['midra']
- self.midDec = soltab_ph_screen.obj._v_attrs['middec']
- H.close()
-
- def get_memory_usage(self, cellsize_deg):
- """
- Returns memory usage per time slot in GB
-
- Parameters
- ----------
- cellsize_deg : float
- Size of one pixel in degrees
- """
- ncpu = self.ncpu
- if ncpu == 0:
- ncpu = misc.nproc()
-
- # Make a test array and find its memory usage
- ximsize = int(self.width_ra / cellsize_deg) # pix
- yimsize = int(self.width_dec / cellsize_deg) # pix
- test_array = np.zeros([1, len(self.freqs_ph), len(self.station_names), 4,
- yimsize, ximsize])
- mem_per_timeslot_gb = test_array.nbytes/1024**3 / 10 # include factor of 10 overhead
-
- # Multiply by the number of CPUs, since each gets a copy
- mem_per_timeslot_gb *= ncpu
-
- return mem_per_timeslot_gb
-
- def make_matrix(self, t_start_index, t_stop_index, freq_ind, stat_ind, cellsize_deg,
- out_dir, ncpu):
- """
- Makes the matrix of values for the given time, frequency, and station indices
-
- Parameters
- ----------
- t_start_index : int
- Index of first time
- t_stop_index : int
- Index of last time
- t_start_index : int
- Index of frequency
- t_stop_index : int
- Index of station
- cellsize_deg : float
- Size of one pixel in degrees
- out_dir : str
- Full path to the output directory
- ncpu : int, optional
- Number of CPUs to use (0 means all)
- """
- # Use global variables to avoid serializing the arrays in the multiprocessing calls
- global screen_ph, screen_amp_xx, screen_amp_yy, pp, x, y, var_dict
-
- # Define various parameters
- N_sources = len(self.source_names)
- N_times = t_stop_index - t_start_index
- N_piercepoints = N_sources
- beta_val = self.beta_val
- r_0 = self.r_0
-
- # Make arrays of pixel coordinates for screen
- # We need to convert the FITS cube pixel coords to screen pixel coords. The FITS cube
- # has self.ra, self.dec at (xsize/2, ysize/2)
- ximsize = int(np.ceil(self.width_ra / cellsize_deg)) # pix
- yimsize = int(np.ceil(self.width_dec / cellsize_deg)) # pix
- w = wcs.WCS(naxis=2)
- w.wcs.crpix = [ximsize/2.0, yimsize/2.0]
- w.wcs.cdelt = np.array([-cellsize_deg, cellsize_deg])
- w.wcs.crval = [self.ra, self.dec]
- w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
- w.wcs.set_pv([(2, 1, 45.0)])
-
- x_fits = list(range(ximsize))
- y_fits = list(range(yimsize))
- ra = []
- dec = []
- for xf, yf in zip(x_fits, y_fits):
- x_y = np.array([[xf, yf]])
- ra.append(w.wcs_pix2world(x_y, 0)[0][0])
- dec.append(w.wcs_pix2world(x_y, 0)[0][1])
- xy, _, _ = stationscreen._getxy(ra, dec, midRA=self.midRA, midDec=self.midDec)
- x = xy[0].T
- y = xy[1].T
- Nx = len(x)
- Ny = len(y)
-
- # Select input data and reorder the axes to get axis order of [dir, time]
- # Input data are [time, freq, ant, dir, pol] for slow amplitudes
- # and [time, freq, ant, dir] for fast phases (scalarphase).
- time_axis = 0
- dir_axis = 1
- screen_ph = np.array(self.vals_ph[t_start_index:t_stop_index, freq_ind, stat_ind, :])
- screen_ph = screen_ph.transpose([dir_axis, time_axis])
- if not self.phase_only:
- screen_amp_xx = np.array(self.vals_amp[t_start_index:t_stop_index, freq_ind, stat_ind, :, 0])
- screen_amp_xx = screen_amp_xx.transpose([dir_axis, time_axis])
- screen_amp_yy = np.array(self.vals_amp[t_start_index:t_stop_index, freq_ind, stat_ind, :, 1])
- screen_amp_yy = screen_amp_yy.transpose([dir_axis, time_axis])
-
- # Process phase screens
- pp = self.pp
- val_shape = (Nx, Ny, N_times)
- var_dict = {}
- shared_val = RawArray('d', int(val_shape[0]*val_shape[1]*val_shape[2]))
- screen_type = 'ph'
- with Pool(processes=ncpu, initializer=init_worker, initargs=(shared_val, val_shape)) as pool:
- result = pool.map(calculate_kl_screen_star,
- zip(range(val_shape[2]), itertools.repeat(N_piercepoints),
- itertools.repeat(beta_val), itertools.repeat(r_0),
- itertools.repeat(screen_type)))
- val_phase = np.frombuffer(shared_val, dtype=np.float64).reshape(val_shape).copy()
-
- # Process amplitude screens
- if not self.phase_only:
- # XX amplitudes
- screen_type = 'xx'
- with Pool(processes=ncpu, initializer=init_worker, initargs=(shared_val, val_shape)) as pool:
- result = pool.map(calculate_kl_screen_star,
- zip(range(val_shape[2]), itertools.repeat(N_piercepoints),
- itertools.repeat(beta_val), itertools.repeat(r_0),
- itertools.repeat(screen_type)))
- val_amp_xx = 10**(np.frombuffer(shared_val, dtype=np.float64).reshape(val_shape).copy())
-
- # YY amplitudes
- screen_type = 'yy'
- with Pool(processes=ncpu, initializer=init_worker, initargs=(shared_val, val_shape)) as pool:
- result = pool.map(calculate_kl_screen_star,
- zip(range(val_shape[2]), itertools.repeat(N_piercepoints),
- itertools.repeat(beta_val), itertools.repeat(r_0),
- itertools.repeat(screen_type)))
- val_amp_yy = 10**(np.frombuffer(shared_val, dtype=np.float64).reshape(val_shape).copy())
-
- # Output data are [RA, DEC, MATRIX, ANTENNA, FREQ, TIME].T
- data = np.zeros((N_times, 4, Ny, Nx))
- if self.phase_only:
- data[:, 0, :, :] = np.cos(val_phase.T)
- data[:, 2, :, :] = np.cos(val_phase.T)
- data[:, 1, :, :] = np.sin(val_phase.T)
- data[:, 3, :, :] = np.sin(val_phase.T)
- else:
- data[:, 0, :, :] = val_amp_xx.T * np.cos(val_phase.T)
- data[:, 2, :, :] = val_amp_yy.T * np.cos(val_phase.T)
- data[:, 1, :, :] = val_amp_xx.T * np.sin(val_phase.T)
- data[:, 3, :, :] = val_amp_yy.T * np.sin(val_phase.T)
-
- return data
-
-
-class VoronoiScreen(Screen):
- """
- Class for Voronoi screens
- """
- def __init__(self, name, h5parm_filename, skymodel_filename, ra, dec, width_ra, width_dec,
- solset_name='sol000', phase_soltab_name='phase000', amplitude_soltab_name=None):
- super(VoronoiScreen, self).__init__(name, h5parm_filename, skymodel_filename, ra, dec, width_ra, width_dec,
- solset_name=solset_name, phase_soltab_name=phase_soltab_name, amplitude_soltab_name=amplitude_soltab_name)
- self.data_rasertize_template = None
-
- def fit(self):
- """
- Fitting is not needed: the input solutions are used directly, after
- referencing the phases to a single station
- """
- # Open solution tables
- H = h5parm(self.input_h5parm_filename)
- solset = H.getSolset(self.input_solset_name)
- soltab_ph = solset.getSoltab(self.input_phase_soltab_name)
- if not self.phase_only:
- soltab_amp = solset.getSoltab(self.input_amplitude_soltab_name)
-
- # Input data are [time, freq, ant, dir, pol] for slow amplitudes
- # and [time, freq, ant, dir] for fast phases (scalarphase). We reference
- # the phases to the station with the least amount of flagged solutions,
- # drawn from the first 10 stations (to ensure it is fairly central)
- self.vals_ph = soltab_ph.val
- ref_ind = misc.get_reference_station(soltab_ph, 10)
- vals_ph_ref = self.vals_ph[:, :, ref_ind, :].copy()
- for i in range(len(soltab_ph.ant)):
- # Subtract phases of reference station
- self.vals_ph[:, :, i, :] -= vals_ph_ref
- self.times_ph = soltab_ph.time
- self.freqs_ph = soltab_ph.freq
- if not self.phase_only:
- self.log_amps = False
- self.vals_amp = soltab_amp.val
- self.times_amp = soltab_amp.time
- self.freqs_amp = soltab_amp.freq
- else:
- self.vals_amp = np.ones_like(self.vals_ph)
- self.times_amp = self.times_ph
- self.freqs_amp = self.freqs_ph
-
- self.source_names = soltab_ph.dir
- self.source_dict = solset.getSou()
- self.source_positions = []
- for source in self.source_names:
- self.source_positions.append(self.source_dict[source])
- self.station_names = soltab_ph.ant
- self.station_dict = solset.getAnt()
- self.station_positions = []
- for station in self.station_names:
- self.station_positions.append(self.station_dict[station])
-
- def get_memory_usage(self, cellsize_deg):
- """
- Returns memory usage per time slot in GB
-
- Parameters
- ----------
- cellsize_deg : float
- Size of one pixel in degrees
- """
- # Make a test array and find its memory usage
- ximsize = int(self.width_ra / cellsize_deg) # pix
- yimsize = int(self.width_dec / cellsize_deg) # pix
- test_array = np.zeros([1, len(self.freqs_ph), len(self.station_names), 4,
- yimsize, ximsize])
- mem_per_timeslot_gb = test_array.nbytes/1024**3 * 10 # include factor of 10 overhead
-
- return mem_per_timeslot_gb
-
- def make_matrix(self, t_start_index, t_stop_index, freq_ind, stat_ind, cellsize_deg,
- out_dir, ncpu):
- """
- Makes the matrix of values for the given time, frequency, and station indices
-
- Parameters
- ----------
- t_start_index : int
- Index of first time
- t_stop_index : int
- Index of last time
- t_start_index : int
- Index of frequency
- t_stop_index : int
- Index of station
- cellsize_deg : float
- Size of one pixel in degrees
- out_dir : str
- Full path to the output directory
- ncpu : int, optional
- Number of CPUs to use (0 means all)
- """
- # Make the template that converts polynomials to a rasterized 2-D image.
- # This only needs to be done once
- if self.data_rasertize_template is None:
- self.make_rasertize_template(cellsize_deg, out_dir)
-
- # Fill the output data array
- data = np.zeros((t_stop_index-t_start_index, 4, self.data_rasertize_template.shape[0],
- self.data_rasertize_template.shape[1]))
- for p, poly in enumerate(self.polygons):
- ind = np.where(self.data_rasertize_template == poly.index+1)
- if not self.phase_only:
- val_amp_xx = self.vals_amp[t_start_index:t_stop_index, freq_ind, stat_ind, poly.index, 0]
- val_amp_yy = self.vals_amp[t_start_index:t_stop_index, freq_ind, stat_ind, poly.index, 1]
- else:
- val_amp_xx = self.vals_amp[t_start_index:t_stop_index, freq_ind, stat_ind, poly.index]
- val_amp_yy = val_amp_xx
- val_phase = self.vals_ph[t_start_index:t_stop_index, freq_ind, stat_ind, poly.index]
- for t in range(t_stop_index-t_start_index):
- data[t, 0, ind[0], ind[1]] = val_amp_xx[t] * np.cos(val_phase[t])
- data[t, 2, ind[0], ind[1]] = val_amp_yy[t] * np.cos(val_phase[t])
- data[t, 1, ind[0], ind[1]] = val_amp_xx[t] * np.sin(val_phase[t])
- data[t, 3, ind[0], ind[1]] = val_amp_yy[t] * np.sin(val_phase[t])
-
- return data
-
- def make_rasertize_template(self, cellsize_deg, out_dir):
- """
- Makes the template that is used to fill the output FITS cube
-
- Parameters
- ----------
- cellsize_deg : float
- Size of one pixel in degrees
- out_dir : str
- Full path to the output directory
- """
- temp_image = os.path.join(out_dir, '{}_template.fits'.format(self.name))
- hdu = self.make_fits_file(temp_image, cellsize_deg, 0, 1, aterm_type='gain')
- data = hdu[0].data
- w = wcs.WCS(hdu[0].header)
- RAind = w.axis_type_names.index('RA')
- Decind = w.axis_type_names.index('DEC')
-
- # Get x, y coords for directions in pixels. We use the input calibration sky
- # model for this, as the patch positions written to the h5parm file by DPPP may
- # be different
- skymod = lsmtool.load(self.input_skymodel_filename)
- source_dict = skymod.getPatchPositions()
- source_positions = []
- for source in self.source_names:
- radecpos = source_dict[source.strip('[]')]
- source_positions.append([radecpos[0].value, radecpos[1].value])
- source_positions = np.array(source_positions)
- ra_deg = source_positions.T[0]
- dec_deg = source_positions.T[1]
-
- xy = []
- for RAvert, Decvert in zip(ra_deg, dec_deg):
- ra_dec = np.array([[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]])
- ra_dec[0][RAind] = RAvert
- ra_dec[0][Decind] = Decvert
- xy.append((w.wcs_world2pix(ra_dec, 0)[0][RAind], w.wcs_world2pix(ra_dec, 0)[0][Decind]))
-
- # Get boundary of tessellation region in pixels
- bounds_deg = [self.ra+self.width_ra/2.0, self.dec-self.width_dec/2.0,
- self.ra-self.width_ra/2.0, self.dec+self.width_dec/2.0]
- ra_dec = np.array([[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]])
- ra_dec[0][RAind] = max(bounds_deg[0], np.max(ra_deg)+0.1)
- ra_dec[0][Decind] = min(bounds_deg[1], np.min(dec_deg)-0.1)
- field_minxy = (w.wcs_world2pix(ra_dec, 0)[0][RAind], w.wcs_world2pix(ra_dec, 0)[0][Decind])
- ra_dec[0][RAind] = min(bounds_deg[2], np.min(ra_deg)-0.1)
- ra_dec[0][Decind] = max(bounds_deg[3], np.max(dec_deg)+0.1)
- field_maxxy = (w.wcs_world2pix(ra_dec, 0)[0][RAind], w.wcs_world2pix(ra_dec, 0)[0][Decind])
-
- if len(xy) == 1:
- # If there is only a single direction, just make a single rectangular polygon
- box = [field_minxy, (field_minxy[0], field_maxxy[1]), field_maxxy, (field_maxxy[0], field_minxy[1]), field_minxy]
- polygons = [shapely.geometry.Polygon(box)]
- else:
- # For more than one direction, tessellate
- # Generate array of outer points used to constrain the facets
- nouter = 64
- means = np.ones((nouter, 2)) * np.array(xy).mean(axis=0)
- offsets = []
- angles = [np.pi/(nouter/2.0)*i for i in range(0, nouter)]
- for ang in angles:
- offsets.append([np.cos(ang), np.sin(ang)])
- radius = 2.0*np.sqrt( (field_maxxy[0]-field_minxy[0])**2 + (field_maxxy[1]-field_minxy[1])**2 )
- scale_offsets = radius * np.array(offsets)
- outer_box = means + scale_offsets
-
- # Tessellate and clip
- points_all = np.vstack([xy, outer_box])
- vor = Voronoi(points_all)
- lines = [
- shapely.geometry.LineString(vor.vertices[line])
- for line in vor.ridge_vertices
- if -1 not in line
- ]
- polygons = [poly for poly in shapely.ops.polygonize(lines)]
-
- # Index polygons to directions
- for i, xypos in enumerate(xy):
- for poly in polygons:
- if poly.contains(Point(xypos)):
- poly.index = i
-
- # Rasterize the polygons to an array, with the value being equal to the
- # polygon's index+1
- data_template = np.ones(data[0, 0, 0, 0, :, :].shape)
- data_rasertize_template = np.zeros(data[0, 0, 0, 0, :, :].shape)
- for poly in polygons:
- verts_xy = poly.exterior.xy
- verts = []
- for x, y in zip(verts_xy[0], verts_xy[1]):
- verts.append((x, y))
- poly_raster = misc.rasterize(verts, data_template.copy()) * (poly.index+1)
- filled = np.where(poly_raster > 0)
- data_rasertize_template[filled] = poly_raster[filled]
- zeroind = np.where(data_rasertize_template == 0)
- if len(zeroind[0]) > 0:
- nonzeroind = np.where(data_rasertize_template != 0)
- data_rasertize_template[zeroind] = si.griddata((nonzeroind[0], nonzeroind[1]), data_rasertize_template[nonzeroind],
- (zeroind[0], zeroind[1]), method='nearest')
- self.data_rasertize_template = data_rasertize_template
- self.polygons = polygons
-
-
-def init_worker(shared_val, val_shape):
- """
- Initializer called when a child process is initialized, responsible
- for storing store shared_val and val_shape in var_dict (a global variable).
-
- See https://research.wmz.ninja/articles/2018/03/on-sharing-large-arrays-when-using-pythons-multiprocessing.html
-
- Parameters
- ----------
- shared_val : array
- RawArray to be shared
- val_shape : tuple
- Shape of shared_val array
- """
- global var_dict
-
- var_dict['shared_val'] = shared_val
- var_dict['val_shape'] = val_shape
-
-
-def calculate_kl_screen_star(inputs):
- """
- Simple helper function for pool.map
- """
- return calculate_kl_screen(*inputs)
-
-
-def calculate_kl_screen(k, N_piercepoints, beta_val, r_0, screen_type):
- """
- Calculates screen images
-
- Parameters
- ----------
- k : int
- Time index
- N_piercepoints : int
- Number of pierce points
- beta_val : float
- Power-law index for phase structure function (5/3 =>
- pure Kolmogorov turbulence)
- r_0 : float
- Scale size of phase fluctuations
- screen_type : string
- Type of screen: 'ph' (phase), 'xx' (XX amplitude) or 'yy' (YY amplitude)
- """
- # Use global variables to avoid serializing the arrays in the multiprocessing calls
- global screen_ph, screen_amp_xx, screen_amp_yy, pp, x, y, var_dict
-
- tmp = np.frombuffer(var_dict['shared_val'], dtype=np.float64).reshape(var_dict['val_shape'])
- if screen_type == 'ph':
- inscreen = screen_ph[:, k]
- if screen_type == 'xx':
- inscreen = screen_amp_xx[:, k]
- if screen_type == 'yy':
- inscreen = screen_amp_yy[:, k]
- f = inscreen.reshape(N_piercepoints)
- for i, xi in enumerate(x):
- for j, yi in enumerate(y):
- p = np.array([xi, yi, 0.0])
- d2 = np.sum(np.square(pp - p), axis=1)
- c = -(d2 / ( r_0**2 ))**(beta_val / 2.0) / 2.0
- tmp[i, j, k] = np.dot(c, f)
diff --git a/rapthor/lib/sector.py b/rapthor/lib/sector.py
index 46b41ec4ef0efe11eb24faa230a9db6b24870dee..876904b0d90fdb2fb597b7ce29fd779f1f906165 100644
--- a/rapthor/lib/sector.py
+++ b/rapthor/lib/sector.py
@@ -141,7 +141,7 @@ class Sector(object):
self.reweight = imaging_parameters['reweight']
self.target_fast_timestep = self.field.fast_timestep_sec
self.target_slow_freqstep = self.field.parset['calibration_specific']['slow_freqstep_hz']
- self.use_screens = self.field.use_screens
+ self.apply_screens = self.field.apply_screens
# Set image size based on current sector polygon
if recalculate_imsize or self.imsize is None:
@@ -151,7 +151,7 @@ class Sector(object):
self.imsize = [int(self.width_ra / self.cellsize_deg),
int(self.width_dec / self.cellsize_deg)]
- if self.use_screens:
+ if self.apply_screens:
# IDG does not yet support rectangular images, so ensure image
# is square
self.imsize = [max(self.imsize), max(self.imsize)]
@@ -242,7 +242,7 @@ class Sector(object):
obs.set_imaging_parameters(self.name, self.cellsize_arcsec, max_peak_smearing,
self.width_ra, self.width_dec,
self.target_fast_timestep, self.target_slow_freqstep,
- self.use_screens)
+ self.apply_screens)
# Set BL-dependent averaging parameters
do_bl_averaging = False # does not yet work with IDG
@@ -422,7 +422,7 @@ class Sector(object):
# Save initial polygon, copy of initial polygon (which potentially will be
# altered later for source avoidance), and buffered version of initial polygon
- # (which includes the padding done by WSClean, needed for aterm generation)
+ # (which includes the padding done by WSClean)
self.initial_poly = poly
self.poly = Polygon(poly)
padding_pix = dec_width_pix*(self.wsclean_image_padding - 1.0)
diff --git a/rapthor/operations/calibrate.py b/rapthor/operations/calibrate.py
index 37964409f18d2cbc5d3f4f86ad17d0ee0708766c..bb619866d011c965eb52a2d9ec16c309e527a488 100644
--- a/rapthor/operations/calibrate.py
+++ b/rapthor/operations/calibrate.py
@@ -33,14 +33,9 @@ class CalibrateDD(Operation):
do_joint_solve = True
else:
do_joint_solve = False
- if self.field.dde_method == 'facets':
- use_facets = True
- else:
- use_facets = False
self.parset_parms = {'rapthor_pipeline_dir': self.rapthor_pipeline_dir,
- 'use_screens': self.field.use_screens,
- 'use_facets': use_facets,
+ 'generate_screens': self.field.generate_screens,
'do_slowgain_solve': self.field.do_slowgain_solve,
'do_joint_solve': do_joint_solve,
'use_scalarphase': self.field.use_scalarphase,
@@ -165,23 +160,6 @@ class CalibrateDD(Operation):
sector_bounds_deg = '{}'.format(self.field.sector_bounds_deg)
sector_bounds_mid_deg = '{}'.format(self.field.sector_bounds_mid_deg)
- # Set the number of chunks to split the solution tables into and define
- # the associated filenames
- if self.field.do_slowgain_solve:
- nsplit = sum(self.field.get_obs_parameters('nsplit_slow'))
- else:
- nsplit = sum(self.field.get_obs_parameters('nsplit_fast'))
- split_outh5parm = ['split_solutions_{}.h5'.format(i) for i in
- range(nsplit)]
-
- # Set the root filenames for the a-term images. We save it to an attribute
- # since it is needed in finalize()
- self.output_aterms_root = ['diagonal_aterms_{}'.format(i) for i in
- range(len(split_outh5parm))]
-
- # Set the type of screen to make
- screen_type = self.field.screen_type
-
# Set the DDECal steps depending on whether baseline-dependent averaging is
# activated (and supported) or not. If BDA is used, an "null" step is also
# added to prevent the writing of the BDA data
@@ -261,9 +239,6 @@ class CalibrateDD(Operation):
'slow_datause': slow_datause,
'sector_bounds_deg': sector_bounds_deg,
'sector_bounds_mid_deg': sector_bounds_mid_deg,
- 'split_outh5parm': split_outh5parm,
- 'output_aterms_root': self.output_aterms_root,
- 'screen_type': screen_type,
'combined_h5parms': self.combined_h5parms,
'fast_antennaconstraint': fast_antennaconstraint,
'slow_antennaconstraint': slow_antennaconstraint,
@@ -347,17 +322,6 @@ class CalibrateDD(Operation):
"""
Finalize this operation
"""
- # Get the filenames of the aterm images (for use in the image operation). The files
- # were written by the 'make_aterms' step and the number of them can vary, depending
- # on the node memory, etc.
- if self.field.use_screens:
- self.field.aterm_image_filenames = []
- for aterms_root in self.output_aterms_root:
- with open(os.path.join(self.pipeline_working_dir, aterms_root+'.txt'), 'r') as f:
- self.field.aterm_image_filenames.extend(f.readlines())
- self.field.aterm_image_filenames = [os.path.join(self.pipeline_working_dir, af.strip())
- for af in self.field.aterm_image_filenames]
-
# Copy the solutions (h5parm files) and report the flagged fraction
dst_dir = os.path.join(self.parset['dir_working'], 'solutions', 'calibrate_{}'.format(self.index))
misc.create_directory(dst_dir)
diff --git a/rapthor/operations/image.py b/rapthor/operations/image.py
index 54cb06eec05878d6a454f7d49c6eb348e03be4ec..d039c79ff40cf1ab57bbdd486190a72d7840ab02 100644
--- a/rapthor/operations/image.py
+++ b/rapthor/operations/image.py
@@ -26,8 +26,8 @@ class Image(Operation):
# Initialize various parameters
self.apply_none = None
self.apply_amplitudes = None
- self.use_screens = None
self.apply_fulljones = None
+ self.apply_screens = None
self.make_image_cube = None
self.normalize_flux_scale = None
self.dde_method = None
@@ -50,8 +50,8 @@ class Image(Operation):
self.apply_none = False
if self.apply_amplitudes is None:
self.apply_amplitudes = self.field.apply_amplitudes
- if self.use_screens is None:
- self.use_screens = self.field.use_screens
+ if self.apply_screens is None:
+ self.apply_screens = self.field.apply_screens
if self.apply_fulljones is None:
self.apply_fulljones = self.field.apply_fulljones
if self.make_image_cube is None:
@@ -61,7 +61,7 @@ class Image(Operation):
if self.dde_method is None:
self.dde_method = self.field.dde_method
if self.use_facets is None:
- self.use_facets = True if self.dde_method == 'facets' else False
+ self.use_facets = True if (self.dde_method == 'full' and not self.apply_screens) else False
if self.image_pol is None:
self.image_pol = self.field.image_pol
if self.save_source_list is None:
@@ -79,7 +79,7 @@ class Image(Operation):
'pipeline_working_dir': self.pipeline_working_dir,
'apply_none': self.apply_none,
'apply_amplitudes': self.apply_amplitudes,
- 'use_screens': self.use_screens,
+ 'apply_screens': self.apply_screens,
'apply_fulljones': self.apply_fulljones,
'make_image_cube': self.make_image_cube,
'normalize_flux_scale': self.normalize_flux_scale,
@@ -242,16 +242,13 @@ class Image(Operation):
nnodes_per_subpipeline = max(1, int(nnodes / nsubpipes) - 1)
self.input_parms.update({'mpi_nnodes': [nnodes_per_subpipeline] * nsectors})
self.input_parms.update({'mpi_cpus_per_task': [self.parset['cluster_specific']['cpus_per_task']] * nsectors})
- if self.use_screens:
- # The following parameters were set by the preceding calibrate operation, where
- # aterm image files were generated. They do not need to be set separately for
- # each sector
- self.input_parms.update({'aterm_image_filenames': CWLFile(self.field.aterm_image_filenames).to_json()})
- elif not self.apply_none:
+ if not self.apply_none:
self.input_parms.update({'h5parm': CWLFile(self.field.h5parm_filename).to_json()})
if self.field.fulljones_h5parm_filename is not None:
self.input_parms.update({'fulljones_h5parm': CWLFile(self.field.fulljones_h5parm_filename).to_json()})
- if self.use_facets:
+ if self.field.apply_screens:
+ self.input_parms.update({'idgcal_h5parm': CWLFile(self.field.idgcal_h5parm_filename).to_json()})
+ elif self.use_facets:
# For faceting, we need inputs for making the ds9 facet region files
self.input_parms.update({'skymodel': CWLFile(self.field.calibration_skymodel_file).to_json()})
ra_mid = []
@@ -360,7 +357,7 @@ class Image(Operation):
shutil.copy(src_filename, dst_filename)
# The output ds9 region file, if made
- if self.field.dde_method == 'facets':
+ if self.use_facets:
dst_dir = os.path.join(self.parset['dir_working'], 'regions', 'image_{}'.format(self.index))
misc.create_directory(dst_dir)
region_filename = '{}_facets_ds9.reg'.format(sector.name)
@@ -422,7 +419,7 @@ class ImageInitial(Image):
# Set parameters as needed
self.apply_none = True
self.apply_amplitudes = False
- self.use_screens = False
+ self.apply_screens = False
self.apply_fulljones = False
self.use_facets = False
self.save_source_list = True
@@ -530,7 +527,7 @@ class ImageNormalize(Image):
# No calibration has yet been done, so set various flags as needed
self.apply_none = True
self.use_facets = False
- self.use_screens = False
+ self.apply_screens = False
super().set_parset_parameters()
def set_input_parameters(self):
diff --git a/rapthor/pipeline/execution/examples/hba_nl/L667526.config.json b/rapthor/pipeline/execution/examples/hba_nl/L667526.config.json
index b0f1f8af0dfab025cc66c2d4fa5ff288fb451444..99f12528b1cc175ae02cf82f47127e507366d8c0 100644
--- a/rapthor/pipeline/execution/examples/hba_nl/L667526.config.json
+++ b/rapthor/pipeline/execution/examples/hba_nl/L667526.config.json
@@ -7,7 +7,7 @@
"calibration": {
},
"imaging": {
- "dde_method": "facets"
+ "dde_method": "full"
},
"cluster": {
"use_container": "True",
diff --git a/rapthor/pipeline/execution/examples/hba_nl/hba_nl.json b/rapthor/pipeline/execution/examples/hba_nl/hba_nl.json
index 5ed0f8337d8e27b9b3e6f5d5ace66cbc5ab3ea6c..b49026011a79d7493087ad5b8ab773e7379fdb34 100644
--- a/rapthor/pipeline/execution/examples/hba_nl/hba_nl.json
+++ b/rapthor/pipeline/execution/examples/hba_nl/hba_nl.json
@@ -14,7 +14,7 @@
"calibration": {
},
"imaging": {
- "dde_method": "facets"
+ "dde_method": "full"
},
"cluster": {
"use_container": "True",
diff --git a/rapthor/pipeline/parsets/calibrate_pipeline.cwl b/rapthor/pipeline/parsets/calibrate_pipeline.cwl
index bb287f015e045e6e22d36b24a46d70703f53b406..765f56bb17a068d0daab47a12b63ef91dbdb889e 100644
--- a/rapthor/pipeline/parsets/calibrate_pipeline.cwl
+++ b/rapthor/pipeline/parsets/calibrate_pipeline.cwl
@@ -10,8 +10,7 @@ doc: |
further unconstrained slow gain calibration to correct for station-to-station
differences. Steps (2) and (3) are skipped if the calibration is phase-only.
This calibration scheme works for both HBA and LBA data. The final products of
- this workflow are solution tables (h5parm files), plots, and a-term screens (FITS
- files).
+ this workflow are solution tables (h5parm files) and plots.
requirements:
ScatterFeatureRequirement: {}
@@ -250,24 +249,6 @@ inputs:
The mid point of the boundary of all imaging sectors in degrees (length = 1).
type: string
- - id: split_outh5parm
- label: Output split solution tables
- doc: |
- The filenames of the output split h5parm solution tables (length = n_obs * n_split).
- type: string[]
-
- - id: output_aterms_root
- label: Output root for a-terms
- doc: |
- The root names of the output a-term images (length = n_obs * n_split).
- type: string[]
-
- - id: screen_type
- label: Type of screen for a-terms
- doc: |
- The screen type to use to derive the a-term images (length = 1).
- type: string
-
- id: max_threads
label: Max number of threads
doc: |
@@ -553,12 +534,6 @@ outputs:
outputSource:
- plot_fast_phase_solutions/plots
type: File[]
-{% if use_screens %}
- - id: diagonal_aterms
- outputSource:
- - merge_aterm_files/output
- type: File[]
-{% endif %}
{% if do_slowgain_solve %}
- id: slow_phase_plots
outputSource:
@@ -1079,11 +1054,6 @@ steps:
- id: outh5parm
source: combined_h5parms
- id: mode
-{% if use_screens %}
- valueFrom: 'p1p2a2'
- - id: reweight
- valueFrom: 'True'
-{% elif use_facets %}
{% if apply_diagonal_solutions %}
valueFrom: 'p1p2a2_diagonal'
{% else %}
@@ -1091,11 +1061,6 @@ steps:
{% endif %}
- id: reweight
valueFrom: 'False'
-{% else %}
- valueFrom: 'p1p2a2_diagonal'
- - id: reweight
- valueFrom: 'False'
-{% endif %}
- id: calibrator_names
source: calibrator_patch_names
- id: calibrator_fluxes
@@ -1116,55 +1081,6 @@ steps:
out:
- id: adjustedh5parm
-{% if use_screens %}
-# start use_screens
-
- - id: split_h5parms
- label: Split solution table
- doc: |
- This step splits the final solution table in time, to enable multi-node
- parallel processing of the a-term images.
- run: {{ rapthor_pipeline_dir }}/steps/split_h5parms.cwl
- in:
- - id: inh5parm
- source: adjust_h5parm_sources/adjustedh5parm
- - id: outh5parms
- source: split_outh5parm
- - id: soltabname
- valueFrom: 'gain000'
- out:
- - id: splith5parms
-
- - id: make_aterms
- label: Make a-term images
- doc: |
- This step makes a-term images from the split final solution tables.
- run: {{ rapthor_pipeline_dir }}/steps/make_aterm_images.cwl
- in:
- - id: h5parm
- source: split_h5parms/splith5parms
- - id: soltabname
- valueFrom: 'gain000'
- - id: screen_type
- source: screen_type
- - id: skymodel
- source: calibration_skymodel_file
- - id: outroot
- source: output_aterms_root
- - id: sector_bounds_deg
- source: sector_bounds_deg
- - id: sector_bounds_mid_deg
- source: sector_bounds_mid_deg
- - id: ncpu
- source: max_threads
- scatter: [h5parm, outroot]
- scatterMethod: dotproduct
- out:
- - id: output_images
-
-{% endif %}
-# end use_screens
-
{% else %}
# start not do_slowgain_solve
@@ -1181,68 +1097,5 @@ steps:
out:
- id: adjustedh5parm
-{% if use_screens %}
-# start use_screens
-
- - id: split_h5parms
- label: Split solution table
- doc: |
- This step splits the final solution table in time, to enable multi-node
- parallel processing of the a-term images.
- run: {{ rapthor_pipeline_dir }}/steps/split_h5parms.cwl
- in:
- - id: inh5parm
- source: adjust_h5parm_sources/adjustedh5parm
- - id: outh5parms
- source: split_outh5parm
- - id: soltabname
- valueFrom: 'phase000'
- out:
- - id: splith5parms
-
- - id: make_aterms
- label: Make a-term images
- doc: |
- This step makes a-term images from the split final solution tables.
- run: {{ rapthor_pipeline_dir }}/steps/make_aterm_images.cwl
- in:
- - id: h5parm
- source: split_h5parms/splith5parms
- - id: soltabname
- valueFrom: 'phase000'
- - id: screen_type
- source: screen_type
- - id: skymodel
- source: calibration_skymodel_file
- - id: outroot
- source: output_aterms_root
- - id: sector_bounds_deg
- source: sector_bounds_deg
- - id: sector_bounds_mid_deg
- source: sector_bounds_mid_deg
- - id: ncpu
- source: max_threads
- scatter: [h5parm, outroot]
- scatterMethod: dotproduct
- out:
- - id: output_images
-
-{% endif %}
-# end use_screens
-
{% endif %}
# end do_slowgain_solve / not do_slowgain_solve
-
-{% if use_screens %}
-
- - id: merge_aterm_files
- in:
- - id: input
- source:
- - make_aterms/output_images
- out:
- - id: output
- run: {{ rapthor_pipeline_dir }}/steps/merge_array_files.cwl
- label: merge_aterm_files
-
-{% endif %}
diff --git a/rapthor/pipeline/parsets/image_pipeline.cwl b/rapthor/pipeline/parsets/image_pipeline.cwl
index 7dda15859f6dcaa5af63aca6df6369323d8245c5..b4015f4a8c7dc468f79b907e53b961e769f9af9b 100644
--- a/rapthor/pipeline/parsets/image_pipeline.cwl
+++ b/rapthor/pipeline/parsets/image_pipeline.cwl
@@ -171,16 +171,6 @@ inputs:
type: int[]
{% endif %}
-{% if use_screens %}
-# start use_screens
- - id: aterm_image_filenames
- label: Filenames of a-terms
- doc: |
- The filenames of the a-term images (length = 1, with n_aterms subelements).
- type: File[]
-
-{% else %}
-# start not use_screens
{% if not apply_none %}
- id: h5parm
label: Filename of h5parm
@@ -199,6 +189,15 @@ inputs:
type: File
{% endif %}
+{% if apply_screens %}
+ - id: idgcal_h5parm
+ label: Filename of h5parm
+ doc: |
+ The filename of the h5parm file with the IDGCal screeen solutions
+ (length = 1).
+ type: File
+{% endif %}
+
{% if use_facets %}
# start use_facets
- id: skymodel
@@ -279,9 +278,6 @@ inputs:
{% endif %}
# end use_facets / not use_facets
-{% endif %}
-# end use_screens / not use_screens
-
- id: channels_out
label: Number of channels
doc: |
@@ -607,12 +603,6 @@ steps:
- id: mpi_nnodes
source: mpi_nnodes
{% endif %}
-{% if use_screens %}
-# start use_screens
- - id: aterm_image_filenames
- source: aterm_image_filenames
-{% else %}
-# start not use_screens
{% if not apply_none %}
- id: h5parm
source: h5parm
@@ -621,6 +611,10 @@ steps:
- id: fulljones_h5parm
source: fulljones_h5parm
{% endif %}
+{% if apply_screens %}
+ - id: idgcal_h5parm
+ source: idgcal_h5parm
+{% endif %}
{% if use_facets %}
# start use_facets
- id: skymodel
@@ -649,8 +643,6 @@ steps:
{% endif %}
{% endif %}
# end use_facets / not use_facets
-{% endif %}
-# end use_screens / not use_screens
- id: channels_out
source: channels_out
- id: deconvolution_channels
@@ -717,28 +709,6 @@ steps:
- id: normalize_h5parm
source: normalize_h5parm
{% endif %}
-{% if use_screens %}
-# start use_screens
- scatter: [obs_filename, prepare_filename, concat_filename, starttime, ntimes,
- image_freqstep, image_timestep, previous_mask_filename, mask_filename,
- phasecenter, ra, dec, image_name, cellsize_deg, wsclean_imsize,
- vertices_file, region_file,
-{% if use_mpi %}
- mpi_cpus_per_task, mpi_nnodes,
-{% endif %}
-{% if make_image_cube %}
- image_cube_name,
-{% endif %}
-{% if normalize_flux_scale %}
- output_source_catalog,
- normalize_h5parm,
-{% endif %}
- channels_out, deconvolution_channels, fit_spectral_pol, wsclean_niter,
- wsclean_nmiter, robust, min_uv_lambda, max_uv_lambda, mgain, do_multiscale,
- taper_arcsec, local_rms_strength, wsclean_mem, auto_mask,
- idg_mode, threshisl, threshpix, dd_psf_grid]
-{% else %}
-# start not use_screens
scatter: [obs_filename, prepare_filename, concat_filename, starttime, ntimes,
image_freqstep, image_timestep, previous_mask_filename, mask_filename,
phasecenter, ra, dec, image_name, cellsize_deg, wsclean_imsize,
@@ -764,8 +734,6 @@ steps:
wsclean_nmiter, robust, min_uv_lambda, max_uv_lambda, mgain, do_multiscale,
taper_arcsec, local_rms_strength, wsclean_mem, auto_mask,
idg_mode, threshisl, threshpix, dd_psf_grid]
-{% endif %}
-# end use_screens / not use_screens
scatterMethod: dotproduct
diff --git a/rapthor/pipeline/parsets/image_sector_pipeline.cwl b/rapthor/pipeline/parsets/image_sector_pipeline.cwl
index 0036458e9d89ea874b8b1c1e1d3640a423b8fb7c..6b3799141207db9a227c68c2b6e0ac592e4f7bbd 100644
--- a/rapthor/pipeline/parsets/image_sector_pipeline.cwl
+++ b/rapthor/pipeline/parsets/image_sector_pipeline.cwl
@@ -138,17 +138,6 @@ inputs:
type: int
{% endif %}
-{% if use_screens %}
-# start use_screens
- - id: aterm_image_filenames
- label: Filenames of a-terms
- doc: |
- The filenames of the a-term images (length = 1, with n_aterms subelements).
- type: File[]
-
-{% else %}
-# start not use_screens
-
{% if not apply_none %}
- id: h5parm
label: Filename of h5parm
@@ -167,6 +156,15 @@ inputs:
type: File
{% endif %}
+{% if apply_screens %}
+ - id: idgcal_h5parm
+ label: Filename of h5parm
+ doc: |
+ The filename of the h5parm file with the IDGCal screeen solutions
+ (length = 1).
+ type: File
+{% endif %}
+
{% if use_facets %}
# start use_facets
- id: skymodel
@@ -244,9 +242,6 @@ inputs:
{% endif %}
# end use_facets / not use_facets
-{% endif %}
-# end use_screens / not use_screens
-
- id: channels_out
label: Number of channels
doc: |
@@ -512,8 +507,8 @@ steps:
This step uses DP3 to prepare the input data for imaging. This involves
averaging, phase shifting, and optionally the application of the
calibration solutions.
-{% if use_screens or use_facets %}
-# start use_screens or use_facets
+{% if apply_screens or use_facets %}
+# start apply_screens or use_facets
{% if apply_fulljones %}
run: {{ rapthor_pipeline_dir }}/steps/prepare_imaging_data_fulljones.cwl
{% else %}
@@ -521,7 +516,7 @@ steps:
{% endif %}
{% else %}
-# start not use_screens and not use_facets
+# start not apply_screens and not use_facets
{% if apply_none %}
run: {{ rapthor_pipeline_dir }}/steps/prepare_imaging_data_no_dde_no_cal.cwl
{% else %}
@@ -537,7 +532,7 @@ steps:
{% endif %}
{% endif %}
-# end use_screens or use_facets / not use_screens and not use_facets
+# end apply_screens or use_facets / not apply_screens and not use_facets
{% if max_cores is not none %}
hints:
@@ -564,7 +559,7 @@ steps:
source: phasecenter
- id: numthreads
source: max_threads
-{% if use_screens or use_facets %}
+{% if apply_screens or use_facets %}
{% if apply_fulljones %}
- id: h5parm
source: fulljones_h5parm
@@ -657,8 +652,8 @@ steps:
doc: |
This step makes an image using WSClean. Direction-dependent effects
can be corrected for using a-term images or facet-based corrections.
-{% if use_screens %}
-# start use_screens
+{% if apply_screens %}
+# start apply_screens
{% if use_mpi %}
run: {{ rapthor_pipeline_dir }}/steps/wsclean_mpi_image_screens.cwl
@@ -667,7 +662,7 @@ steps:
{% endif %}
{% else %}
-# start not use_screens
+# start not apply_screens
{% if use_facets %}
# start use_facets
@@ -679,7 +674,7 @@ steps:
{% endif %}
{% else %}
-# start not use_facets and not use_screens (i.e., use no_dde)
+# start not use_facets and not apply_screens (i.e., use no_dde)
{% if use_mpi %}
run: {{ rapthor_pipeline_dir }}/steps/wsclean_mpi_image_no_dde.cwl
@@ -691,7 +686,7 @@ steps:
# end use no_dde
{% endif %}
-# end use_screens / not use_screens
+# end apply_screens / not apply_screens
{% if max_cores is not none %}
hints:
@@ -706,14 +701,14 @@ steps:
source: image_name
- id: mask
source: premask/maskimg
+{% if apply_screens %}
+ - id: h5parm
+ source: idgcal_h5parm
+{% endif %}
{% if use_mpi %}
- id: nnodes
source: mpi_nnodes
{% endif %}
-{% if use_screens %}
- - id: aterm_images
- source: aterm_image_filenames
-{% endif %}
{% if use_facets %}
- id: h5parm
source: h5parm
diff --git a/rapthor/pipeline/steps/make_aterm_images.cwl b/rapthor/pipeline/steps/make_aterm_images.cwl
deleted file mode 100644
index 97aa8143dfa145c8fbc180d80841fbfac6ed7275..0000000000000000000000000000000000000000
--- a/rapthor/pipeline/steps/make_aterm_images.cwl
+++ /dev/null
@@ -1,92 +0,0 @@
-cwlVersion: v1.2
-class: CommandLineTool
-baseCommand: [make_aterm_images.py]
-label: Make a-term images
-doc: |
- This tool makes FITS a-term images from the input solution table. The images
- are suitable for use with WSClean+IDG.
-
-requirements:
- - class: InlineJavascriptRequirement
- - class: InitialWorkDirRequirement
- listing:
- - entry: $(inputs.h5parm)
- writable: true
-
-arguments:
- - '--smooth_deg=0.1'
-
-inputs:
- - id: h5parm
- label: Input solution table
- doc: |
- The filename of the input h5parm file.
- type: File
- inputBinding:
- position: 0
- - id: soltabname
- label: Input soltab name
- doc: |
- The name of the input soltab.
- type: string
- inputBinding:
- prefix: --soltabname=
- separate: false
- - id: screen_type
- label: Screen type
- doc: |
- The type of the screen to generate.
- type: string
- inputBinding:
- prefix: --screen_type=
- separate: false
- - id: outroot
- label: Output root name
- doc: |
- The root of the filenames of the output images.
- type: string
- inputBinding:
- prefix: --outroot=
- separate: false
- - id: skymodel
- label: Input sky model
- doc: |
- The filename of the input sky model file.
- type: File
- inputBinding:
- prefix: --skymodel=
- separate: false
- - id: sector_bounds_deg
- label: Bounds of imaging sectors
- doc: |
- The global bounds of the imaging sectors in deg.
- type: string
- inputBinding:
- prefix: --bounds_deg=
- separate: false
- - id: sector_bounds_mid_deg
- label: Midpoint of imaging bounds
- doc: |
- The global midpoint of the imaging sectors in deg.
- type: string
- inputBinding:
- prefix: --bounds_mid_deg=
- separate: false
- - id: ncpu
- label: Number of CPUs
- doc: |
- The number of CPUs / cores to use.
- type: int
- inputBinding:
- prefix: --ncpu=
- separate: false
-
-outputs:
- - id: output_images
- type: File[]
- outputBinding:
- glob: '$(inputs.outroot)*'
-
-hints:
- - class: DockerRequirement
- dockerPull: 'astronrd/rapthor'
diff --git a/rapthor/pipeline/steps/split_h5parms.cwl b/rapthor/pipeline/steps/split_h5parms.cwl
deleted file mode 100644
index 7c78c98109ef1e5d89913b7716d060bd67d7e0bc..0000000000000000000000000000000000000000
--- a/rapthor/pipeline/steps/split_h5parms.cwl
+++ /dev/null
@@ -1,49 +0,0 @@
-cwlVersion: v1.2
-class: CommandLineTool
-baseCommand: [split_h5parms.py]
-label: Splits an h5parm
-doc: |
- This tool splits the input h5parm in time to produce multiple
- output h5parms.
-
-requirements:
- InlineJavascriptRequirement: {}
-
-inputs:
- - id: inh5parm
- label: Input solution table
- doc: |
- The filename of the input h5parm file.
- type: File
- inputBinding:
- position: 0
- - id: outh5parms
- label: Output solution tables
- doc: |
- The filenames of the output h5parm files.
- type: string[]
- inputBinding:
- position: 1
- itemSeparator: ","
- - id: soltabname
- label: Solution table name
- doc: |
- The name of the solution table to split.
- type: string
- inputBinding:
- prefix: --soltabname=
- separate: false
-
-outputs:
- - id: splith5parms
- label: Output solution tables
- doc: |
- The filenames of the output h5parm files. The value is taken from the input
- parameter "outh5parms"
- type: File[]
- outputBinding:
- glob: $(inputs.outh5parms)
-
-hints:
- - class: DockerRequirement
- dockerPull: 'astronrd/rapthor'
diff --git a/rapthor/pipeline/steps/wsclean_image_screens.cwl b/rapthor/pipeline/steps/wsclean_image_screens.cwl
index c3e9c98de248026d56d20b8dcf6c1e946434565c..bdcf0e41db9e05cc3abea6ea14a7bc3546c63b86 100644
--- a/rapthor/pipeline/steps/wsclean_image_screens.cwl
+++ b/rapthor/pipeline/steps/wsclean_image_screens.cwl
@@ -3,19 +3,22 @@ class: CommandLineTool
baseCommand: [wsclean]
label: Make an image
doc: |
- This tool makes an image using WSClean+IDG with a-term corrections.
+ This tool makes an image using WSClean+IDG with a-term corrections generated by
+ IDGCal.
requirements:
- class: InitialWorkDirRequirement
listing:
- entryname: aterm_plus_beam.cfg
- # Note: WSClean requires that the aterm image filenames be input as part of an
- # aterm config file (and not directly on the command line). Therefore, a config
- # file is made here that contains the filenames defined in the aterm_images
- # input parameter. Also, the required beam parameters are set here
+ # Note: WSClean requires that the h5parm filename be input as part of an aterm
+ # config file (and not directly on the command line). Therefore, a config file is
+ # made here that contains the filename defined in the h5parm input parameter.
+ # Also, the required beam parameters are set here
entry: |
- aterms = [diagonal, beam]
- diagonal.images = [$(inputs.aterm_images.map( (e,i) => (e.path) ).join(' '))]
+ aterms = [idgcalsolutions, beam]
+ idgcalsolutions.type = h5parm
+ idgcalsolutions.files = [$(inputs.h5parm)]
+ idgcalsolutions.update_interval = 8
beam.differential = true
beam.update_interval = 120
beam.usechannelfreq = true
@@ -75,14 +78,11 @@ inputs:
type: File
inputBinding:
prefix: -fits-mask
- - id: aterm_images
- label: Filenames of aterm files
+ - id: h5parm
+ label: h5parm filename
doc: |
- The filenames of the a-term image files. These filenames are not used directly in the
- WSClean call (they are read by WSClean from the aterm config file, defined in the
- requirements section above), hence the value is set to "null" (which results in
- nothing being added to the command for this input).
- type: File[]
+ The filename of the h5parm containing the solutions to apply to correct for DDEs.
+ type: File
inputBinding:
valueFrom: null
- id: wsclean_imsize
diff --git a/rapthor/pipeline/steps/wsclean_mpi_image_screens.cwl b/rapthor/pipeline/steps/wsclean_mpi_image_screens.cwl
index b616fcc053b0f783a0342f7be13cc2a0e5aacb2f..8233fdc1421dc9cf00f37d9cb1bb2cd65a7ddc64 100644
--- a/rapthor/pipeline/steps/wsclean_mpi_image_screens.cwl
+++ b/rapthor/pipeline/steps/wsclean_mpi_image_screens.cwl
@@ -14,13 +14,15 @@ requirements:
- class: InitialWorkDirRequirement
listing:
- entryname: aterm_plus_beam.cfg
- # Note: WSClean requires that the aterm image filenames be input as part of an
- # aterm config file (and not directly on the command line). Therefore, a config
- # file is made here that contains the filenames defined in the aterm_images
- # input parameter. Also, the required beam parameters are set here
+ # Note: WSClean requires that the h5parm filename be input as part of an aterm
+ # config file (and not directly on the command line). Therefore, a config file is
+ # made here that contains the filename defined in the h5parm input parameter.
+ # Also, the required beam parameters are set here
entry: |
- aterms = [diagonal, beam]
- diagonal.images = [$(inputs.aterm_images.map( (e,i) => (e.path) ).join(' '))]
+ aterms = [idgcalsolutions, beam]
+ idgcalsolutions.type = h5parm
+ idgcalsolutions.files = [$(inputs.h5parm)]
+ idgcalsolutions.update_interval = 8
beam.differential = true
beam.update_interval = 120
beam.usechannelfreq = true
@@ -73,14 +75,8 @@ inputs:
type: File
inputBinding:
prefix: -fits-mask
- - id: aterm_images
- label: Filenames of aterm files
- doc: |
- The filenames of the a-term image files. These filenames are not used directly in the
- WSClean call (they are read by WSClean from the aterm config file, defined in the
- requirements section above), hence the value is set to "null" (which results in
- nothing being added to the command for this input).
- type: File[]
+ - id: h5parm
+ type: File
inputBinding:
valueFrom: null
- id: wsclean_imsize
diff --git a/rapthor/process.py b/rapthor/process.py
index 6fac18abe9f82a6ae104f97bd502426660b5930b..68a36714b8f55ba171b7332a0b83cb2606aed7e2 100644
--- a/rapthor/process.py
+++ b/rapthor/process.py
@@ -107,6 +107,14 @@ def run(parset_file, logging_level='info'):
log.info("Using a data fraction of {0:.2f}".format(parset['final_data_fraction']))
if field.make_quv_images:
log.info("Stokes I, Q, U, and V images will be made")
+ if field.dde_mode == 'hybrid':
+ log.info("Screens will be used for calibration and imaging (since dde_mode = "
+ "'hybrid' and this is the final iteration)")
+ if final_step['peel_outliers']:
+ # Currently, when screens are used peeling cannot be done
+ log.warning("Peeling of outliers is currently not supported when using "
+ "screens. Peeling will be skipped")
+ final_step['peel_outliers'] = False
# Set the data chunking to match the longest solution interval set in
# the strategy
@@ -157,6 +165,9 @@ def run_steps(field, steps, final=False):
# Calibrate
if field.do_calibrate:
+ # Set whether screens should be generated
+ field.generate_screens = True if (field.dde_mode == 'hybrid' and final) else False
+
if field.peel_non_calibrator_sources:
# Predict and subtract non-calibrator sources before calibration
op = PredictNC(field, cycle_number)
@@ -174,7 +185,8 @@ def run_steps(field, steps, final=False):
op.run()
# Predict and subtract the sector models
- if field.do_predict:
+ # Note: DD predict is not yet supported when screens are used
+ if field.do_predict and not field.generate_screens:
op = PredictDD(field, cycle_number)
op.run()
@@ -183,6 +195,9 @@ def run_steps(field, steps, final=False):
# Set the Stokes polarizations for imaging
field.image_pol = 'IQUV' if (field.make_quv_images and final) else 'I'
+ # Set whether screens should be applied
+ field.apply_screens = True if (field.dde_mode == 'hybrid' and final) else False
+
op = Image(field, cycle_number)
op.run()
diff --git a/rapthor/scripts/make_aterm_images.py b/rapthor/scripts/make_aterm_images.py
deleted file mode 100755
index d97ac31e70190319f61c81af51a81695aa68a329..0000000000000000000000000000000000000000
--- a/rapthor/scripts/make_aterm_images.py
+++ /dev/null
@@ -1,128 +0,0 @@
-#!/usr/bin/env python3
-"""
-Script to make a-term images from solutions
-"""
-from argparse import ArgumentParser, RawTextHelpFormatter
-import os
-from rapthor.lib.screen import KLScreen, VoronoiScreen
-from losoto.h5parm import h5parm
-
-
-def main(h5parmfile, soltabname='phase000', screen_type='tessellated', outroot='',
- bounds_deg=None, bounds_mid_deg=None, skymodel=None,
- solsetname='sol000', padding_fraction=1.4, cellsize_deg=0.2,
- smooth_deg=0, interp_kind='nearest', ncpu=0):
- """
- Make a-term FITS images
-
- Parameters
- ----------
- h5parmfile : str
- Filename of h5parm
- soltabname : str, optional
- Name of soltab to use. If "gain" is in the name, phase and amplitudes are used
- screen_type : str, optional
- Kind of screen to use: 'tessellated' (simple Voronoi tessellation) or 'kl'
- (Karhunen-Lo`eve transform)
- outroot : str, optional
- Root of filename of output FITS file (root+'_0.fits')
- bounds_deg : list, optional
- List of [maxRA, minDec, minRA, maxDec] for image bounds
- bounds_mid_deg : list, optional
- List of [RA, Dec] for midpoint of image bounds
- skymodel : str, optional
- Filename of calibration sky model (needed for patch positions)
- solsetname : str, optional
- Name of solset
- padding_fraction : float, optional
- Fraction of total size to pad with (e.g., 0.2 => 20% padding all around)
- cellsize_deg : float, optional
- Cellsize of output image
- smooth_deg : float, optional
- Size of smoothing kernel in degrees to apply
- interp_kind : str, optional
- Kind of interpolation to use. Can be any supported by scipy.interpolate.interp1d
- ncpu : int, optional
- Number of CPUs to use (0 means all)
-
- Returns
- -------
- result : dict
- Dict with list of FITS files
- """
- if 'gain' in soltabname:
- # We have scalarphase and XX+YY amplitudes
- soltab_amp = soltabname.replace('gain', 'amplitude')
- soltab_ph = soltabname.replace('gain', 'phase')
- else:
- # We have scalarphase only
- soltab_amp = None
- soltab_ph = soltabname
-
- if type(bounds_deg) is str:
- bounds_deg = [float(f.strip()) for f in bounds_deg.strip('[]').split(';')]
- if type(bounds_mid_deg) is str:
- bounds_mid_deg = [float(f.strip()) for f in bounds_mid_deg.strip('[]').split(';')]
- if padding_fraction is not None:
- padding_fraction = float(padding_fraction)
- padding_ra = (bounds_deg[2] - bounds_deg[0]) * (padding_fraction - 1.0)
- padding_dec = (bounds_deg[3] - bounds_deg[1]) * (padding_fraction - 1.0)
- bounds_deg[0] -= padding_ra
- bounds_deg[1] -= padding_dec
- bounds_deg[2] += padding_ra
- bounds_deg[3] += padding_dec
- cellsize_deg = float(cellsize_deg)
- smooth_deg = float(smooth_deg)
- smooth_pix = smooth_deg / cellsize_deg
- if screen_type == 'kl':
- # No need to smooth KL screens
- smooth_pix = 0.0
-
- # Check whether we just have one direction. If so, force screen_type to 'tessellated'
- # as it can handle this case and KL screens can't
- H = h5parm(h5parmfile)
- solset = H.getSolset(solsetname)
- soltab = solset.getSoltab(soltab_ph)
- source_names = soltab.dir[:]
- if len(source_names) == 1:
- screen_type = 'tessellated'
- H.close()
-
- # Fit screens and make a-term images
- width_deg = bounds_deg[3] - bounds_deg[1] # Use Dec difference and force square images
- rootname = os.path.basename(outroot)
- if screen_type == 'kl':
- screen = KLScreen(rootname, h5parmfile, skymodel, bounds_mid_deg[0], bounds_mid_deg[1],
- width_deg, width_deg, solset_name=solsetname,
- phase_soltab_name=soltab_ph, amplitude_soltab_name=soltab_amp)
- elif screen_type == 'tessellated':
- screen = VoronoiScreen(rootname, h5parmfile, skymodel, bounds_mid_deg[0], bounds_mid_deg[1],
- width_deg, width_deg, solset_name=solsetname,
- phase_soltab_name=soltab_ph, amplitude_soltab_name=soltab_amp)
- screen.process(ncpu=ncpu)
- outdir = os.path.dirname(outroot)
- screen.write(outdir, cellsize_deg, smooth_pix=smooth_pix, interp_kind=interp_kind, ncpu=ncpu)
-
-
-if __name__ == '__main__':
- descriptiontext = "Make a-term images from solutions.\n"
-
- parser = ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
- parser.add_argument('h5parmfile', help='Filename of input h5parm')
- parser.add_argument('--soltabname', help='Name of soltab', type=str, default='phase000')
- parser.add_argument('--screen_type', help='Type of screen', type=str, default='tessellated')
- parser.add_argument('--outroot', help='Root of output images', type=str, default='')
- parser.add_argument('--bounds_deg', help='Bounds list in deg', type=str, default=None)
- parser.add_argument('--bounds_mid_deg', help='Bounds mid list in deg', type=str, default=None)
- parser.add_argument('--skymodel', help='Filename of sky model', type=str, default=None)
- parser.add_argument('--solsetname', help='Solset name', type=str, default='sol000')
- parser.add_argument('--padding_fraction', help='Padding fraction', type=float, default=1.4)
- parser.add_argument('--cellsize_deg', help='Cell size in deg', type=float, default=0.2)
- parser.add_argument('--smooth_deg', help='Smooth scale in degree', type=float, default=0.0)
- parser.add_argument('--ncpu', help='Number of CPUs to use', type=int, default=0)
- args = parser.parse_args()
- main(args.h5parmfile, soltabname=args.soltabname, screen_type=args.screen_type,
- outroot=args.outroot, bounds_deg=args.bounds_deg,
- bounds_mid_deg=args.bounds_mid_deg, skymodel=args.skymodel,
- solsetname=args.solsetname, padding_fraction=args.padding_fraction,
- cellsize_deg=args.cellsize_deg, smooth_deg=args.smooth_deg, ncpu=args.ncpu)
diff --git a/rapthor/scripts/split_h5parms.py b/rapthor/scripts/split_h5parms.py
deleted file mode 100755
index be9dc105f22788bd3cb11ecf90816c3c5a581755..0000000000000000000000000000000000000000
--- a/rapthor/scripts/split_h5parms.py
+++ /dev/null
@@ -1,161 +0,0 @@
-#!/usr/bin/env python3
-"""
-Script to split an h5parm
-"""
-from argparse import ArgumentParser, RawTextHelpFormatter
-from losoto.h5parm import h5parm
-import os
-import shutil
-import numpy as np
-from rapthor.lib import miscellaneous as misc
-
-
-def main(inh5parm, outh5parms, soltabname='phase000', insolset='sol000'):
- """
- Combines two h5parms
-
- Parameters
- ----------
- inh5parm : str
- Filename of h5parm to split
- outh5parms : str
- Filenames of split h5parms as comma-separated string
- soltabname : str, optional
- Name of soltab to use. If "gain" is in the name, phase and amplitudes are used
- insolset : str, optional
- Name of solset to split
- """
- output_h5parm_list = misc.string2list(outh5parms)
- nchunks = len(output_h5parm_list)
- if nchunks == 1:
- # If there is only one output file, just copy the input
- shutil.copy(inh5parm, output_h5parm_list[0])
- return
-
- # Read input table
- h5 = h5parm(inh5parm, readonly=True)
- solset = h5.getSolset(insolset)
- if 'gain' in soltabname:
- soltab_amp = solset.getSoltab(soltabname.replace('gain', 'amplitude'))
- soltab_ph = solset.getSoltab(soltabname.replace('gain', 'phase'))
- else:
- soltab_ph = solset.getSoltab(soltabname)
- pointingNames = []
- antennaNames = []
- pointingDirections = []
- antennaPositions = []
- ants = solset.getAnt()
- sous = solset.getSou()
- for k, v in list(sous.items()):
- if k not in pointingNames:
- pointingNames.append(k)
- pointingDirections.append(v)
- for k, v in list(ants.items()):
- if k not in antennaNames:
- antennaNames.append(k)
- antennaPositions.append(v)
-
- # Read in times
- times_fast = soltab_ph.time
- if 'gain' in soltabname:
- times_slow = soltab_amp.time
-
- # Identify any gaps in time and put initial breaks there. We use median()
- # instead of min() to find the solution interval (timewidth) because the
- # division in time used during calibration to allow processing on multiple
- # nodes can occasionally result in a few smaller solution intervals
- delta_times = times_fast[1:] - times_fast[:-1] # time at center of solution interval
- timewidth = np.median(delta_times)
- gaps = np.where(delta_times > timewidth*1.2)
- gaps_ind = gaps[0] + 1
- gaps_ind = np.append(gaps_ind, np.array([len(times_fast)]))
-
- # Add additional breaks to reach the desired number of chunks
- if len(gaps_ind) >= nchunks:
- gaps_ind = gaps_ind[:nchunks]
-
- while len(gaps_ind) < nchunks:
- # Find the largest existing gap
- g_num_largest = 0
- g_size_largest = 0
- g_start = 0
- for g_num, g_stop in enumerate(gaps_ind):
- if g_stop - g_start > g_size_largest:
- g_num_largest = g_num
- g_size_largest = g_stop - g_start
- g_start = g_stop
-
- # Now split largest gap into two equal parts
- if g_num_largest == 0:
- g_start = 0
- else:
- g_start = gaps_ind[g_num_largest-1]
- g_stop = gaps_ind[g_num_largest]
- new_gap = g_start + int((g_stop - g_start) / 2)
- gaps_ind = np.insert(gaps_ind, g_num_largest, np.array([new_gap]))
-
- gaps_sec = []
- for i, gind in enumerate(gaps_ind):
- if i == nchunks-1:
- gaps_sec.append(times_fast[-1])
- else:
- gaps_sec.append(times_fast[gind])
-
- # Fill the output files
- for i, outh5file in enumerate(output_h5parm_list):
- if os.path.exists(outh5file):
- os.remove(outh5file)
- outh5 = h5parm(outh5file, readonly=False)
- solsetOut = outh5.makeSolset('sol000')
-
- # Store phases
- if i == 0:
- startval = times_fast[0]
- else:
- startval = gaps_sec[i-1]
- if i == nchunks-1:
- endval = times_fast[-1]
- else:
- endval = gaps_sec[i] - 0.5 # subtract 0.5 sec to ensure "[)" range
- soltab_ph.setSelection(time={'min': startval, 'max': endval, 'step': 1})
- solsetOut.makeSoltab('phase', 'phase000', axesNames=soltab_ph.getAxesNames(),
- axesVals=[soltab_ph.getAxisValues(a) for a in soltab_ph.getAxesNames()],
- vals=soltab_ph.getValues(retAxesVals=False),
- weights=soltab_ph.getValues(weight=True, retAxesVals=False))
- soltab_ph.clearSelection()
-
- # Store amps
- if 'gain' in soltabname:
- if i == 0:
- startval = times_slow[0]
- else:
- startval = gaps_sec[i-1]
- if i == nchunks-1:
- endval = times_slow[-1]
- else:
- endval = gaps_sec[i] - 0.5 # subtract 0.5 sec to ensure "[)" range
- soltab_amp.setSelection(time={'min': startval, 'max': endval, 'step': 1})
- solsetOut.makeSoltab('amplitude', 'amplitude000', axesNames=soltab_amp.getAxesNames(),
- axesVals=[soltab_amp.getAxisValues(a) for a in soltab_amp.getAxesNames()],
- vals=soltab_amp.getValues(retAxesVals=False),
- weights=soltab_amp.getValues(weight=True, retAxesVals=False))
- soltab_amp.clearSelection()
-
- # Store metadata
- sourceTable = solsetOut.obj._f_get_child('source')
- antennaTable = solsetOut.obj._f_get_child('antenna')
- antennaTable.append(list(zip(*(antennaNames, antennaPositions))))
- sourceTable.append(list(zip(*(pointingNames, pointingDirections))))
- outh5.close()
-
-
-if __name__ == '__main__':
- descriptiontext = "Combine two h5parms.\n"
-
- parser = ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
- parser.add_argument('inh5parm', help='name of input h5parm')
- parser.add_argument('outh5parms', help='name of output h5parm')
- parser.add_argument('--soltabname', help='name of the soltab', type=str, default='phase000')
- args = parser.parse_args()
-
- main(args.inh5parm, args.outh5parms, soltabname=args.soltabname)
diff --git a/rapthor/settings/defaults.json b/rapthor/settings/defaults.json
index 2ee1b245aaf442dcbc53b7b8903c1b2a2c3a3849..664016df546c0b5e77715a5b981e624b1cdc96b6 100644
--- a/rapthor/settings/defaults.json
+++ b/rapthor/settings/defaults.json
@@ -17,7 +17,8 @@
"final_data_fraction": 1.0,
"input_h5parm": null,
"input_fulljones_h5parm": null,
- "facet_layout": null
+ "facet_layout": null,
+ "dde_mode": "faceting"
},
"calibration": {
"use_included_skymodels": false,
@@ -57,8 +58,7 @@
"taper_arcsec": 0.0,
"local_rms_strength": 1.0,
"do_multiscale_clean": true,
- "dde_method": "facets",
- "screen_type": "kl",
+ "dde_method": "full",
"save_visibilities": false,
"idg_mode": "cpu",
"mem_gb": 0,
diff --git a/rapthor/settings/defaults.parset b/rapthor/settings/defaults.parset
index 178d4835dc452a72ecfbf116e559e83be82ba7d7..f958cc1c8f0eab14d21864cb179ac03d02d37fac 100644
--- a/rapthor/settings/defaults.parset
+++ b/rapthor/settings/defaults.parset
@@ -73,6 +73,16 @@ input_fulljones_h5parm = None
# between cycles)
facet_layout = None
+# Mode to use to derive and correct for direction-dependent effects: faceting or
+# hybrid (default = faceting). If "faceting", Voronoi faceting is used
+# throughout the processing. If "hybrid", faceting is used only during the self
+# calibration steps; in the final cycle (done after self calibration has been
+# completed successfully), IDGCal is used during calibration to generate smooth 2-D
+# screens that are then applied by WSClean in the final imaging step.
+#
+# Note: the hybrid mode is not yet available; it will be enabled in a future update.
+dde_mode = faceting
+
[calibration]
# If one of the included sky models (see rapthor/skymodels) is within 2 *
@@ -140,24 +150,22 @@ local_rms_strength = 0.8
do_multiscale_clean = True
# Method to use to correct for direction-dependent effects during imaging:
-# "none", "facets", or "screens" (default = facets). If "none", the solutions
-# closest to the image centers will be used. If "facets", Voronoi faceting is
-# used. If "screens", smooth 2-D screens are used; the type of screen to use can
-# be specified with screen_type: "tessellated" (simple, smoothed tessellated
-# screens) or "kl" (Karhunen-Lo`eve screens) (default = kl)
-dde_method = facets
-screen_type = kl
+# single or full (default = full). If "single", a single, direction-independent
+# solution (i.e., constant across the image sector) will be applied for each
+# sector. In this case, the solution applied is the one in the direction closest
+# to each sector center. If "full", the full, direction-dependent solutions
+# (either facets or screens) are applied
+dde_method = full
# Save visibilities used for imaging (default = False). If True, the imaging MS
# files will be saved, with the the direction-independent full-Jones solutions,
# if available, applied. Note, however, that the direction-dependent solutions
-# will not be applied unless dde_method = "none", in which case the solutions
+# will not be applied unless dde_method = "single", in which case the solutions
# closest to the image centers are used
save_visibilities = False
# IDG (image domain gridder) mode to use in WSClean (default = cpu). The mode
-# can be "cpu" or "hybrid". Note that IDG is only used when dde_method = "none"
-# or "screens"
+# can be "cpu" or "hybrid". Note that IDG is only used when dde_mode = "hybrid"
idg_mode = cpu
# Maximum memory in GB (per node) to use for WSClean jobs (default = 0 = 100%)
diff --git a/test/resources/rapthor_complete.parset.template b/test/resources/rapthor_complete.parset.template
index ad24cadfe3efa08c204612a34077014927eaab17..a7384a41506473ce9c46ad8ac3800ae29f3b445b 100644
--- a/test/resources/rapthor_complete.parset.template
+++ b/test/resources/rapthor_complete.parset.template
@@ -17,6 +17,7 @@ final_data_fraction = 0.8
input_h5parm = /my/data/dirdependent.h5
input_fulljones_h5parm = /my/data/fulljones.h5
facet_layout = /my/data/facets.reg
+dde_mode = hybrid
[calibration]
use_included_skymodels = True
@@ -55,8 +56,7 @@ max_uv_lambda = 10000.0
taper_arcsec = 2.0
local_rms_strength = 0.5
do_multiscale_clean = False
-dde_method = screens
-screen_type = kl
+dde_method = single
save_visibilities = True
idg_mode = gpu
mem_gb = 128.0
diff --git a/test/resources/rapthor_complete.parset_dict.template b/test/resources/rapthor_complete.parset_dict.template
index ff03e7c200442f11236bca7b45d9d38e6efef24a..f4205950c4b2f092b898cd2a5bcb3722453b695a 100644
--- a/test/resources/rapthor_complete.parset_dict.template
+++ b/test/resources/rapthor_complete.parset_dict.template
@@ -56,7 +56,7 @@
'imaging_specific': {'apply_diagonal_solutions': False,
'cellsize_arcsec': 2.25,
'dd_psf_grid': [2, 2],
- 'dde_method': 'screens',
+ 'dde_method': 'single',
'do_multiscale_clean': False,
'grid_center_dec': 35.50875555555555,
'grid_center_ra': 220.25785,
@@ -74,7 +74,6 @@
'reweight': False,
'robust': -0.3,
'save_visibilities': True,
- 'screen_type': 'kl',
'sector_center_dec_list': [35.50875555555555,
37.36579444444445],
'sector_center_ra_list': [220.25785, 213.3468083333333],
@@ -91,4 +90,5 @@
'mss': ['${input_ms}'],
'regroup_input_skymodel': False,
'selfcal_data_fraction': 0.4,
- 'strategy': 'selfcal'}
+ 'strategy': 'selfcal',
+ 'dde_mode': 'hybrid'}
diff --git a/test/resources/rapthor_minimal.parset_dict.template b/test/resources/rapthor_minimal.parset_dict.template
index 28a241ecca51a1de441e98b040ca28ccb632bee5..09107bbb7a6724eb0450b10db27b863a63430f69 100644
--- a/test/resources/rapthor_minimal.parset_dict.template
+++ b/test/resources/rapthor_minimal.parset_dict.template
@@ -56,7 +56,7 @@
'imaging_specific': {'apply_diagonal_solutions': True,
'cellsize_arcsec': 1.25,
'dd_psf_grid': [0, 0],
- 'dde_method': 'facets',
+ 'dde_method': 'full',
'do_multiscale_clean': True,
'grid_center_dec': None,
'grid_center_ra': None,
@@ -74,7 +74,6 @@
'reweight': False,
'robust': -0.5,
'save_visibilities': False,
- 'screen_type': 'kl',
'sector_center_dec_list': [],
'sector_center_ra_list': [],
'sector_width_dec_deg_list': [],
@@ -90,4 +89,5 @@
'mss': ['${input_ms}'],
'regroup_input_skymodel': True,
'selfcal_data_fraction': 0.2,
- 'strategy': 'selfcal'}
+ 'strategy': 'selfcal',
+ 'dde_mode': 'faceting'}
diff --git a/test/resources/test.parset b/test/resources/test.parset
index 4153a636785dfc1e49c2da5268b3cd617b1d886d..8b46a9632c305ad9e7446759ca4e3f073d4f5db5 100644
--- a/test/resources/test.parset
+++ b/test/resources/test.parset
@@ -26,7 +26,6 @@ grid_center_dec = +57d24m39.0s
grid_nsectors_ra = 1
idg_mode = cpu
reweight = True
-use_screens = True
use_mpi = True
[cluster]
diff --git a/test/scripts/compare_workflow_results.py b/test/scripts/compare_workflow_results.py
index 44b6026dcc1d922b2b8ce5a398d189f4d74d90c7..bc85188a4f09fcf1e710740c12deda033da02f72 100755
--- a/test/scripts/compare_workflow_results.py
+++ b/test/scripts/compare_workflow_results.py
@@ -284,7 +284,7 @@ def compare_results(dcmp, atol, rtol, verbosity):
left, right, rtol=rtol, verbosity=verbosity
):
agree = False
- if ("sky" in ext or ext == ".txt") and "diagonal_aterms" not in root:
+ if ("sky" in ext or ext == ".txt"):
if not skymodel_files_are_similar(
left, right, atol=atol, rtol=rtol, verbosity=verbosity
):
diff --git a/test/test_cwl.py b/test/test_cwl.py
index 87dfc2919bde725d0465aee26eddda3a9757249a..e171ea2ca2860908942f264a14885129dc3282d4 100644
--- a/test/test_cwl.py
+++ b/test/test_cwl.py
@@ -15,8 +15,8 @@ def generate_and_validate(tmp_path, operation, parms, templ, sub_templ=None):
the same way that `rapthor.lib.operation.Operation.setup()` does this.
Validate the workflow file using `cwltool`.
"""
- if parms.get("use_facets") and parms.get("use_screens"):
- pytest.skip("'use_facets' and 'use_screens' cannot be enabled both")
+ if parms.get("use_facets") and parms.get("apply_screens"):
+ pytest.skip("'use_facets' and 'apply_screens' cannot be enabled both")
pipeline_working_dir = tmp_path / "pipelines" / operation
pipeline_working_dir.mkdir(parents=True, exist_ok=True)
parset = pipeline_working_dir / "pipeline_parset.cwl"
@@ -65,7 +65,7 @@ def test_concatenate_workflow(
generate_and_validate(tmp_path, operation, parms, templ)
-@pytest.mark.parametrize("use_screens", (False, True))
+@pytest.mark.parametrize("generate_screens", (False, True))
@pytest.mark.parametrize("use_facets", (False, True))
@pytest.mark.parametrize("do_slowgain_solve", (False, True))
@pytest.mark.parametrize("do_joint_solve", (False, True))
@@ -74,7 +74,7 @@ def test_concatenate_workflow(
@pytest.mark.parametrize("max_cores", (None, 8))
def test_calibrate_workflow(
tmp_path,
- use_screens,
+ generate_screens,
use_facets,
do_slowgain_solve,
do_joint_solve,
@@ -90,7 +90,7 @@ def test_calibrate_workflow(
operation = "calibrate"
templ = rapthor.lib.operation.env_parset.get_template("calibrate_pipeline.cwl")
parms = {
- "use_screens": use_screens,
+ "generate_screens": generate_screens,
"use_facets": use_facets,
"do_slowgain_solve": do_slowgain_solve,
"do_joint_solve": do_joint_solve,
@@ -176,7 +176,7 @@ def test_predict_nc_workflow(tmp_path, max_cores, apply_solutions, apply_amplitu
@pytest.mark.parametrize("apply_amplitudes", (False, True))
-@pytest.mark.parametrize("use_screens", (False, True))
+@pytest.mark.parametrize("apply_screens", (False, True))
@pytest.mark.parametrize("use_facets", (False, True))
@pytest.mark.parametrize("peel_bright_sources", (False, True))
@pytest.mark.parametrize("max_cores", (None, 8))
@@ -186,7 +186,7 @@ def test_predict_nc_workflow(tmp_path, max_cores, apply_solutions, apply_amplitu
def test_image_workflow(
tmp_path,
apply_amplitudes,
- use_screens,
+ apply_screens,
use_facets,
peel_bright_sources,
max_cores,
@@ -206,7 +206,7 @@ def test_image_workflow(
)
parms = {
"apply_amplitudes": apply_amplitudes,
- "use_screens": use_screens,
+ "apply_screens": apply_screens,
"use_facets": use_facets,
"peel_bright_sources": peel_bright_sources,
"max_cores": max_cores,
diff --git a/test/test_screen.py b/test/test_screen.py
deleted file mode 100644
index 55a5f6049ffa8d47e4a61bf482e0a6b6f1f17b90..0000000000000000000000000000000000000000
--- a/test/test_screen.py
+++ /dev/null
@@ -1,50 +0,0 @@
-import unittest
-import os
-import requests
-from rapthor.lib.screen import KLScreen
-
-class TestScreen(unittest.TestCase):
- @classmethod
- def downloadms(self, filename):
- url = 'https://support.astron.nl/software/ci_data/rapthor/' + filename
- r = requests.get(url)
- f = open('resources/' + filename, 'wb')
- f.write(r.content)
- f.close()
-
- @classmethod
- def setUpClass(self):
- cwd = os.getcwd()
- if not cwd.endswith('test'):
- raise SystemExit('Please run this test from the test directory!')
- testh5name = 'split_solutions_0.h5'
- if (not os.path.exists('resources/' + testh5name)):
- print('downloading h5 file')
- self.downloadms(testh5name)
- else:
- print('h5 file found')
-
- self.screen = KLScreen('testscreen', 'resources/split_solutions_0.h5', 'resources/calibration_skymodel.txt', 0.0, 0.0, 1.0, 1.0)
-
- @classmethod
- def tearDownClass(self):
- pass
-
- def test_screen(self):
- self.assertEqual(self.screen.name, 'testscreen')
- self.assertEqual(self.screen.width_ra, 1.0)
-
- def test_fit(self):
- self.screen.ncpu = 1
- self.screen.fit()
- self.assertAlmostEqual(self.screen.midDec, 57.2615, delta=1.0E-4)
-
-
-def suite():
- suite = unittest.TestSuite()
- suite.addTest(TestScreen('test_screen'))
- suite.addTest(TestScreen('test_fit'))
- return suite
-
-if __name__ == '__main__':
- unittest.main()