workflows/arlexecute/imaging/imaging_arlexecute.py
def deconvolve_list_arlexecute_workflow(dirty_list, psf_list, model_imagelist, prefix='', mask=None, **kwargs):
    """Create a graph for deconvolution, adding to the model

    :param dirty_list:
    :param psf_list:
    :param model_imagelist:
    :param prefix: Informative prefix to log messages
    :param mask: Mask for deconvolution
    :param kwargs: Parameters for functions in components
    :return: graph for the deconvolution
    """
    nchan = len(dirty_list)
    # Number of moments. 1 is the sum.
    nmoment = get_parameter(kwargs, "nmoment", 1)
    
    if get_parameter(kwargs, "use_serial_clean", False):
        from workflows.serial.imaging.imaging_serial import deconvolve_list_serial_workflow
        return arlexecute.execute(deconvolve_list_serial_workflow, nout=nchan) \
            (dirty_list, psf_list, model_imagelist, prefix=prefix, mask=mask, **kwargs)
    
    def deconvolve(dirty, psf, model, facet, gthreshold, msk=None):
        if prefix == '':
            lprefix = "facet %d" % facet
        else:
            lprefix = "%s, facet %d" % (prefix, facet)
        
        if nmoment > 0:
            moment0 = calculate_image_frequency_moments(dirty)
            this_peak = numpy.max(numpy.abs(moment0.data[0, ...])) / dirty.data.shape[0]
        else:
            ref_chan = dirty.data.shape[0] // 2
            this_peak = numpy.max(numpy.abs(dirty.data[ref_chan, ...]))
        
        if this_peak > 1.1 * gthreshold:
            kwargs['threshold'] = gthreshold
            result, _ = deconvolve_cube(dirty, psf, prefix=lprefix, mask=msk, **kwargs)
            
            if result.data.shape[0] == model.data.shape[0]:
                result.data += model.data
            return result
        else:
            return copy_image(model)
    
    deconvolve_facets = get_parameter(kwargs, 'deconvolve_facets', 1)
    deconvolve_overlap = get_parameter(kwargs, 'deconvolve_overlap', 0)
    deconvolve_taper = get_parameter(kwargs, 'deconvolve_taper', None)
    if deconvolve_facets > 1 and deconvolve_overlap > 0:
        deconvolve_number_facets = (deconvolve_facets - 2) ** 2
    else:
        deconvolve_number_facets = deconvolve_facets ** 2
    
    scattered_channels_facets_model_list = \
        [arlexecute.execute(image_scatter_facets, nout=deconvolve_number_facets)(m, facets=deconvolve_facets,
                                                                                 overlap=deconvolve_overlap,
                                                                                 taper=deconvolve_taper)
         for m in model_imagelist]
    scattered_facets_model_list = [
        arlexecute.execute(image_gather_channels, nout=1)([scattered_channels_facets_model_list[chan][facet]
                                                           for chan in range(nchan)])
        for facet in range(deconvolve_number_facets)]
    
    # Scatter the separate channel images into deconvolve facets and then gather channels for each facet.
    # This avoids constructing the entire spectral cube.
    # i.e. SCATTER BY FACET then SCATTER BY CHANNEL
    dirty_list_trimmed = arlexecute.execute(remove_sumwt, nout=nchan)(dirty_list)
    scattered_channels_facets_dirty_list = \
        [arlexecute.execute(image_scatter_facets, nout=deconvolve_number_facets)(d, facets=deconvolve_facets,
                                                                                 overlap=deconvolve_overlap,
                                                                                 taper=deconvolve_taper)
         for d in dirty_list_trimmed]
    scattered_facets_dirty_list = [
        arlexecute.execute(image_gather_channels, nout=1)([scattered_channels_facets_dirty_list[chan][facet]
                                                           for chan in range(nchan)])
        for facet in range(deconvolve_number_facets)]
    
    psf_list_trimmed = arlexecute.execute(remove_sumwt, nout=nchan)(psf_list)
    
    def extract_psf(psf, facets):
        spsf = create_empty_image_like(psf)
        cx = spsf.shape[3] // 2
        cy = spsf.shape[2] // 2
        wx = spsf.shape[3] // facets
        wy = spsf.shape[2] // facets
        xbeg = cx - wx // 2
        xend = cx + wx // 2
        ybeg = cy - wy // 2
        yend = cy + wy // 2
        spsf.data = psf.data[..., ybeg:yend, xbeg:xend]
        spsf.wcs.wcs.crpix[0] -= xbeg
        spsf.wcs.wcs.crpix[1] -= ybeg
        return spsf
    
    psf_list_trimmed = [arlexecute.execute(extract_psf)(p, deconvolve_facets) for p in psf_list_trimmed]
    psf_centre = arlexecute.execute(image_gather_channels, nout=1)([psf_list_trimmed[chan]
                                                                                   for chan in range(nchan)])
    
    # Work out the threshold. Need to find global peak over all dirty_list images
    threshold = get_parameter(kwargs, "threshold", 0.0)
    fractional_threshold = get_parameter(kwargs, "fractional_threshold", 0.1)
    nmoment = get_parameter(kwargs, "nmoment", 1)
    use_moment0 = nmoment > 0
    
    # Find the global threshold. This uses the peak in the average on the frequency axis since we
    # want to use it in a stopping criterion in a moment clean
    global_threshold = arlexecute.execute(threshold_list, nout=1)(scattered_facets_dirty_list, threshold,
                                                                  fractional_threshold,
                                                                  use_moment0=use_moment0, prefix=prefix)
    
    facet_list = numpy.arange(deconvolve_number_facets).astype('int')
    if mask is None:
        scattered_results_list = [
            arlexecute.execute(deconvolve, nout=1)(d, psf_centre, m, facet,
                                                   global_threshold)
            for d, m, facet in zip(scattered_facets_dirty_list, scattered_facets_model_list, facet_list)]
    else:
        mask_list = \
            arlexecute.execute(image_scatter_facets, nout=deconvolve_number_facets)(mask,
                                                                                    facets=deconvolve_facets,
                                                                                    overlap=deconvolve_overlap)
        scattered_results_list = [
            arlexecute.execute(deconvolve, nout=1)(d, psf_centre, m, facet,
                                                   global_threshold, msk)
            for d, m, facet, msk in
            zip(scattered_facets_dirty_list, scattered_facets_model_list, facet_list, mask_list)]
    
    # We want to avoid constructing the entire cube so we do the inverse of how we got here:
    # i.e. SCATTER BY CHANNEL then GATHER BY FACET
    # Gather the results back into one image, correcting for overlaps as necessary. The taper function is is used to
    # feather the facets together
    # gathered_results_list = arlexecute.execute(image_gather_facets, nout=1)(scattered_results_list,
    #                                                                         deconvolve_model_imagelist,
    #                                                                         facets=deconvolve_facets,
    #                                                                         overlap=deconvolve_overlap,
    #                                                                         taper=deconvolve_taper)
    # result_list = arlexecute.execute(image_scatter_channels, nout=nchan)(gathered_results_list, subimages=nchan)
    
    scattered_channel_results_list = [arlexecute.execute(image_scatter_channels, nout=nchan)(scat, subimages=nchan)
                                      for scat in scattered_results_list]
    
    # The structure is now [[channels] for facets]. We do the reverse transpose to the one above.
    result_list = [arlexecute.execute(image_gather_facets, nout=1)([scattered_channel_results_list[facet][chan]
                                                                    for facet in range(deconvolve_number_facets)],
                                                                   model_imagelist[chan], facets=deconvolve_facets,
                                                                   overlap=deconvolve_overlap)
                   for chan in range(nchan)]
    
    return arlexecute.optimize(result_list)
Create a graph for deconvolution, adding to the model

    :param dirty_list:
    :param psf_list:
    :param model_imagelist:
    :param prefix: Informative prefix to log messages
    :param mask: Mask for deconvolution
    :param kwargs: Parameters for functions in components
    :return: graph for the deconvolution
    
Create a graph for deconvolution, adding to the model
Create a graph for deconvolution adding to the model
def deconvolve_list_arlexecute_workflow(dirty_list, psf_list, model_imagelist, prefix='', mask=None, **kwargs):
    nchan = len(dirty_list)
    nmoment = get_parameter(kwargs, "nmoment", 1)
    if get_parameter(kwargs, "use_serial_clean", False):
        from workflows.serial.imaging.imaging_serial import deconvolve_list_serial_workflow
        return arlexecute.execute(deconvolve_list_serial_workflow, nout=nchan) \
            (dirty_list, psf_list, model_imagelist, prefix=prefix, mask=mask, **kwargs)
    def deconvolve(dirty, psf, model, facet, gthreshold, msk=None):
        if prefix == '':
            lprefix = "facet %d" % facet
        else:
            lprefix = "%s, facet %d" % (prefix, facet)
        if nmoment > 0:
            moment0 = calculate_image_frequency_moments(dirty)
            this_peak = numpy.max(numpy.abs(moment0.data[0, ...])) / dirty.data.shape[0]
        else:
            ref_chan = dirty.data.shape[0] // 2
            this_peak = numpy.max(numpy.abs(dirty.data[ref_chan, ...]))
        if this_peak > 1.1 * gthreshold:
            kwargs['threshold'] = gthreshold
            result, _ = deconvolve_cube(dirty, psf, prefix=lprefix, mask=msk, **kwargs)
            if result.data.shape[0] == model.data.shape[0]:
                result.data += model.data
            return result
        else:
            return copy_image(model)
    deconvolve_facets = get_parameter(kwargs, 'deconvolve_facets', 1)
    deconvolve_overlap = get_parameter(kwargs, 'deconvolve_overlap', 0)
    deconvolve_taper = get_parameter(kwargs, 'deconvolve_taper', None)
    if deconvolve_facets > 1 and deconvolve_overlap > 0:
        deconvolve_number_facets = (deconvolve_facets - 2) ** 2
    else:
        deconvolve_number_facets = deconvolve_facets ** 2
    scattered_channels_facets_model_list = \
        [arlexecute.execute(image_scatter_facets, nout=deconvolve_number_facets)(m, facets=deconvolve_facets,
                                                                                 overlap=deconvolve_overlap,
                                                                                 taper=deconvolve_taper)
         for m in model_imagelist]
    scattered_facets_model_list = [
        arlexecute.execute(image_gather_channels, nout=1)([scattered_channels_facets_model_list[chan][facet]
                                                           for chan in range(nchan)])
        for facet in range(deconvolve_number_facets)]
    dirty_list_trimmed = arlexecute.execute(remove_sumwt, nout=nchan)(dirty_list)
    scattered_channels_facets_dirty_list = \
        [arlexecute.execute(image_scatter_facets, nout=deconvolve_number_facets)(d, facets=deconvolve_facets,
                                                                                 overlap=deconvolve_overlap,
                                                                                 taper=deconvolve_taper)
         for d in dirty_list_trimmed]
    scattered_facets_dirty_list = [
        arlexecute.execute(image_gather_channels, nout=1)([scattered_channels_facets_dirty_list[chan][facet]
                                                           for chan in range(nchan)])
        for facet in range(deconvolve_number_facets)]
    psf_list_trimmed = arlexecute.execute(remove_sumwt, nout=nchan)(psf_list)
    def extract_psf(psf, facets):
        spsf = create_empty_image_like(psf)
        cx = spsf.shape[3] // 2
        cy = spsf.shape[2] // 2
        wx = spsf.shape[3] // facets
        wy = spsf.shape[2] // facets
        xbeg = cx - wx // 2
        xend = cx + wx // 2
        ybeg = cy - wy // 2
        yend = cy + wy // 2
        spsf.data = psf.data[..., ybeg:yend, xbeg:xend]
        spsf.wcs.wcs.crpix[0] -= xbeg
        spsf.wcs.wcs.crpix[1] -= ybeg
        return spsf
    psf_list_trimmed = [arlexecute.execute(extract_psf)(p, deconvolve_facets) for p in psf_list_trimmed]
    psf_centre = arlexecute.execute(image_gather_channels, nout=1)([psf_list_trimmed[chan]
                                                                                   for chan in range(nchan)])
    threshold = get_parameter(kwargs, "threshold", 0.0)
    fractional_threshold = get_parameter(kwargs, "fractional_threshold", 0.1)
    nmoment = get_parameter(kwargs, "nmoment", 1)
    use_moment0 = nmoment > 0
    global_threshold = arlexecute.execute(threshold_list, nout=1)(scattered_facets_dirty_list, threshold,
                                                                  fractional_threshold,
                                                                  use_moment0=use_moment0, prefix=prefix)
    facet_list = numpy.arange(deconvolve_number_facets).astype('int')
    if mask is None:
        scattered_results_list = [
            arlexecute.execute(deconvolve, nout=1)(d, psf_centre, m, facet,
                                                   global_threshold)
            for d, m, facet in zip(scattered_facets_dirty_list, scattered_facets_model_list, facet_list)]
    else:
        mask_list = \
            arlexecute.execute(image_scatter_facets, nout=deconvolve_number_facets)(mask,
                                                                                    facets=deconvolve_facets,
                                                                                    overlap=deconvolve_overlap)
        scattered_results_list = [
            arlexecute.execute(deconvolve, nout=1)(d, psf_centre, m, facet,
                                                   global_threshold, msk)
            for d, m, facet, msk in
            zip(scattered_facets_dirty_list, scattered_facets_model_list, facet_list, mask_list)]
    scattered_channel_results_list = [arlexecute.execute(image_scatter_channels, nout=nchan)(scat, subimages=nchan)
                                      for scat in scattered_results_list]
    result_list = [arlexecute.execute(image_gather_facets, nout=1)([scattered_channel_results_list[facet][chan]
                                                                    for facet in range(deconvolve_number_facets)],
                                                                   model_imagelist[chan], facets=deconvolve_facets,
                                                                   overlap=deconvolve_overlap)
                   for chan in range(nchan)]
    return arlexecute.optimize(result_list)
def deconvolve_list_arlexecute_workflow ( dirty_list , psf_list , model_imagelist , prefix = '' , mask = None , ** kwargs ) : nchan = len ( dirty_list ) nmoment = get_parameter ( kwargs , "nmoment" , 1 ) if get_parameter ( kwargs , "use_serial_clean" , False ) : from workflows . serial . imaging . imaging_serial import deconvolve_list_serial_workflow return arlexecute . execute ( deconvolve_list_serial_workflow , nout = nchan ) ( dirty_list , psf_list , model_imagelist , prefix = prefix , mask = mask , ** kwargs ) def deconvolve ( dirty , psf , model , facet , gthreshold , msk = None ) : if prefix == '' : lprefix = "facet %d" % facet else : lprefix = "%s, facet %d" % ( prefix , facet ) if nmoment > 0 : moment0 = calculate_image_frequency_moments ( dirty ) this_peak = numpy . max ( numpy . abs ( moment0 . data [ 0 , ... ] ) ) / dirty . data . shape [ 0 ] else : ref_chan = dirty . data . shape [ 0 ] // 2 this_peak = numpy . max ( numpy . abs ( dirty . data [ ref_chan , ... ] ) ) if this_peak > 1.1 * gthreshold : kwargs [ 'threshold' ] = gthreshold result , _ = deconvolve_cube ( dirty , psf , prefix = lprefix , mask = msk , ** kwargs ) if result . data . shape [ 0 ] == model . data . shape [ 0 ] : result . data += model . data return result else : return copy_image ( model ) deconvolve_facets = get_parameter ( kwargs , 'deconvolve_facets' , 1 ) deconvolve_overlap = get_parameter ( kwargs , 'deconvolve_overlap' , 0 ) deconvolve_taper = get_parameter ( kwargs , 'deconvolve_taper' , None ) if deconvolve_facets > 1 and deconvolve_overlap > 0 : deconvolve_number_facets = ( deconvolve_facets - 2 ) ** 2 else : deconvolve_number_facets = deconvolve_facets ** 2 scattered_channels_facets_model_list = [ arlexecute . execute ( image_scatter_facets , nout = deconvolve_number_facets ) ( m , facets = deconvolve_facets , overlap = deconvolve_overlap , taper = deconvolve_taper ) for m in model_imagelist ] scattered_facets_model_list = [ arlexecute . execute ( image_gather_channels , nout = 1 ) ( [ scattered_channels_facets_model_list [ chan ] [ facet ] for chan in range ( nchan ) ] ) for facet in range ( deconvolve_number_facets ) ] dirty_list_trimmed = arlexecute . execute ( remove_sumwt , nout = nchan ) ( dirty_list ) scattered_channels_facets_dirty_list = [ arlexecute . execute ( image_scatter_facets , nout = deconvolve_number_facets ) ( d , facets = deconvolve_facets , overlap = deconvolve_overlap , taper = deconvolve_taper ) for d in dirty_list_trimmed ] scattered_facets_dirty_list = [ arlexecute . execute ( image_gather_channels , nout = 1 ) ( [ scattered_channels_facets_dirty_list [ chan ] [ facet ] for chan in range ( nchan ) ] ) for facet in range ( deconvolve_number_facets ) ] psf_list_trimmed = arlexecute . execute ( remove_sumwt , nout = nchan ) ( psf_list ) def extract_psf ( psf , facets ) : spsf = create_empty_image_like ( psf ) cx = spsf . shape [ 3 ] // 2 cy = spsf . shape [ 2 ] // 2 wx = spsf . shape [ 3 ] // facets wy = spsf . shape [ 2 ] // facets xbeg = cx - wx // 2 xend = cx + wx // 2 ybeg = cy - wy // 2 yend = cy + wy // 2 spsf . data = psf . data [ ... , ybeg : yend , xbeg : xend ] spsf . wcs . wcs . crpix [ 0 ] -= xbeg spsf . wcs . wcs . crpix [ 1 ] -= ybeg return spsf psf_list_trimmed = [ arlexecute . execute ( extract_psf ) ( p , deconvolve_facets ) for p in psf_list_trimmed ] psf_centre = arlexecute . execute ( image_gather_channels , nout = 1 ) ( [ psf_list_trimmed [ chan ] for chan in range ( nchan ) ] ) threshold = get_parameter ( kwargs , "threshold" , 0.0 ) fractional_threshold = get_parameter ( kwargs , "fractional_threshold" , 0.1 ) nmoment = get_parameter ( kwargs , "nmoment" , 1 ) use_moment0 = nmoment > 0 global_threshold = arlexecute . execute ( threshold_list , nout = 1 ) ( scattered_facets_dirty_list , threshold , fractional_threshold , use_moment0 = use_moment0 , prefix = prefix ) facet_list = numpy . arange ( deconvolve_number_facets ) . astype ( 'int' ) if mask is None : scattered_results_list = [ arlexecute . execute ( deconvolve , nout = 1 ) ( d , psf_centre , m , facet , global_threshold ) for d , m , facet in zip ( scattered_facets_dirty_list , scattered_facets_model_list , facet_list ) ] else : mask_list = arlexecute . execute ( image_scatter_facets , nout = deconvolve_number_facets ) ( mask , facets = deconvolve_facets , overlap = deconvolve_overlap ) scattered_results_list = [ arlexecute . execute ( deconvolve , nout = 1 ) ( d , psf_centre , m , facet , global_threshold , msk ) for d , m , facet , msk in zip ( scattered_facets_dirty_list , scattered_facets_model_list , facet_list , mask_list ) ] scattered_channel_results_list = [ arlexecute . execute ( image_scatter_channels , nout = nchan ) ( scat , subimages = nchan ) for scat in scattered_results_list ] result_list = [ arlexecute . execute ( image_gather_facets , nout = 1 ) ( [ scattered_channel_results_list [ facet ] [ chan ] for facet in range ( deconvolve_number_facets ) ] , model_imagelist [ chan ] , facets = deconvolve_facets , overlap = deconvolve_overlap ) for chan in range ( nchan ) ] return arlexecute . optimize ( result_list )
Create a graph for deconvolution, adding to the model
Create a graph for deconvolution adding to the model
"""Create a graph for deconvolution, adding to the model

    :param dirty_list:
    :param psf_list:
    :param model_imagelist:
    :param prefix: Informative prefix to log messages
    :param mask: Mask for deconvolution
    :param kwargs: Parameters for functions in components
    :return: graph for the deconvolution
    """ # Number of moments. 1 is the sum. # Scatter the separate channel images into deconvolve facets and then gather channels for each facet. # This avoids constructing the entire spectral cube. # i.e. SCATTER BY FACET then SCATTER BY CHANNEL # Work out the threshold. Need to find global peak over all dirty_list images # Find the global threshold. This uses the peak in the average on the frequency axis since we # want to use it in a stopping criterion in a moment clean # We want to avoid constructing the entire cube so we do the inverse of how we got here: # i.e. SCATTER BY CHANNEL then GATHER BY FACET # Gather the results back into one image, correcting for overlaps as necessary. The taper function is is used to # feather the facets together # gathered_results_list = arlexecute.execute(image_gather_facets, nout=1)(scattered_results_list, #                                                                         deconvolve_model_imagelist, #                                                                         facets=deconvolve_facets, #                                                                         overlap=deconvolve_overlap, #                                                                         taper=deconvolve_taper) # result_list = arlexecute.execute(image_scatter_channels, nout=nchan)(gathered_results_list, subimages=nchan) # The structure is now [[channels] for facets]. We do the reverse transpose to the one above.
graph for the deconvolution
graph for the deconvolution
dirty_list
None
psf_list
None
model_imagelist
None
prefix
Informative prefix to log messages
Informative prefix to log messages
mask
Mask for deconvolution
Mask for deconvolution
kwargs
Parameters for functions in components
Parameters for functions in components