<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Demonstrate full circle wide field imaging
This include prediction of components, inversion, point source fitting. We will compare the output images with the input models, looking for closeness in flux and position.
Step1:  Construct the SKA1-LOW core configuration
Step2:  Use Dask
Step3:  We create the visibility. This just makes the uvw, time, antenna1, antenna2, weight columns in a table
Step4:  Fill in the visibility with exact calculation of a number of point sources
Step5:  Check that the skycoordinate and image coordinate system are consistent by finding the point sources.
Step6:  Make the convolution function
Step7:  Predict the visibility using the different approaches.
Step8:  Make the image using the different approaches. We will evaluate the results using a number of plots
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import os
import sys
sys.path.append(os.path.join('..', '..'))
from data_models.parameters import arl_path
results_dir = arl_path('test_results')
from matplotlib import pylab
pylab.rcParams['figure.figsize'] = (10.0, 10.0)
pylab.rcParams['image.cmap'] = 'rainbow'
from matplotlib import pyplot as plt
import numpy
from astropy.coordinates import SkyCoord
from astropy.time import Time
from astropy import units as u
from astropy.wcs.utils import pixel_to_skycoord
from data_models.polarisation import PolarisationFrame
from wrappers.serial.image.iterators import image_raster_iter
from wrappers.serial.visibility.base import create_visibility
from wrappers.serial.visibility.operations import sum_visibility
from wrappers.serial.visibility.iterators import vis_timeslices, vis_wslices
from wrappers.serial.simulation.configurations import create_named_configuration
from wrappers.serial.skycomponent.operations import create_skycomponent, find_skycomponents, \
    find_nearest_skycomponent, insert_skycomponent
from wrappers.serial.image.operations import show_image, export_image_to_fits, qa_image, smooth_image
from wrappers.serial.simulation.configurations import create_named_configuration
from wrappers.serial.imaging.base import advise_wide_field, create_image_from_visibility, \
    predict_skycomponent_visibility
from wrappers.arlexecute.griddata.kernels import create_awterm_convolutionfunction
from wrappers.arlexecute.griddata.convolution_functions import apply_bounding_box_convolutionfunction
# Use workflows for imaging
from wrappers.arlexecute.execution_support.arlexecute import arlexecute
from workflows.shared.imaging.imaging_shared import imaging_contexts
from workflows.arlexecute.imaging.imaging_arlexecute import predict_list_arlexecute_workflow, \
    invert_list_arlexecute_workflow
import logging
log = logging.getLogger()
log.setLevel(logging.DEBUG)
log.addHandler(logging.StreamHandler(sys.stdout))
pylab.rcParams['figure.figsize'] = (12.0, 12.0)
pylab.rcParams['image.cmap'] = 'rainbow'
Explanation: Demonstrate full circle wide field imaging
This include prediction of components, inversion, point source fitting. We will compare the output images with the input models, looking for closeness in flux and position.
End of explanation
lowcore = create_named_configuration('LOWBD2-CORE')
Explanation: Construct the SKA1-LOW core configuration
End of explanation
arlexecute.set_client(use_dask=True)
Explanation: Use Dask
End of explanation
times = numpy.array([-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0]) * (numpy.pi / 12.0)
frequency = numpy.array([1e8])
channel_bandwidth = numpy.array([1e6])
reffrequency = numpy.max(frequency)
phasecentre = SkyCoord(ra=+15.0 * u.deg, dec=-45.0 * u.deg, frame='icrs', equinox='J2000')
vt = create_visibility(lowcore, times, frequency, channel_bandwidth=channel_bandwidth,
                       weight=1.0, phasecentre=phasecentre, 
                       polarisation_frame=PolarisationFrame('stokesI'))
advice = advise_wide_field(vt, wprojection_planes=1)
Explanation: We create the visibility. This just makes the uvw, time, antenna1, antenna2, weight columns in a table
End of explanation
vt.data['vis'] *= 0.0
npixel=256
model = create_image_from_visibility(vt, npixel=npixel, cellsize=0.001, nchan=1, 
                                     polarisation_frame=PolarisationFrame('stokesI'))
centre = model.wcs.wcs.crpix-1
spacing_pixels = npixel // 8
log.info('Spacing in pixels = %s' % spacing_pixels)
spacing = model.wcs.wcs.cdelt * spacing_pixels
locations = [-3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5]
original_comps = []
# We calculate the source positions in pixels and then calculate the
# world coordinates to put in the skycomponent description
for iy in locations:
    for ix in locations:
        if ix >= iy:
            p = int(round(centre[0] + ix * spacing_pixels * numpy.sign(model.wcs.wcs.cdelt[0]))), \
                int(round(centre[1] + iy * spacing_pixels * numpy.sign(model.wcs.wcs.cdelt[1])))
            sc = pixel_to_skycoord(p[0], p[1], model.wcs)
            log.info("Component at (%f, %f) [0-rel] %s" % (p[0], p[1], str(sc)))
            flux = numpy.array([[100.0 + 2.0 * ix + iy * 20.0]])
            comp = create_skycomponent(flux=flux, frequency=frequency, direction=sc, 
                                       polarisation_frame=PolarisationFrame('stokesI'))
            original_comps.append(comp)
            insert_skycomponent(model, comp)
predict_skycomponent_visibility(vt, original_comps)
        
cmodel = smooth_image(model) 
show_image(cmodel)
plt.title("Smoothed model image")
plt.show()
Explanation: Fill in the visibility with exact calculation of a number of point sources
End of explanation
comps = find_skycomponents(cmodel, fwhm=1.0, threshold=10.0, npixels=5)
plt.clf()
for i in range(len(comps)):
    ocomp, sep = find_nearest_skycomponent(comps[i].direction, original_comps)
    plt.plot((comps[i].direction.ra.value  - ocomp.direction.ra.value)/cmodel.wcs.wcs.cdelt[0], 
             (comps[i].direction.dec.value - ocomp.direction.dec.value)/cmodel.wcs.wcs.cdelt[1], 
             '.', color='r')  
plt.xlabel('delta RA (pixels)')
plt.ylabel('delta DEC (pixels)')
plt.title("Recovered - Original position offsets")
plt.show()
Explanation: Check that the skycoordinate and image coordinate system are consistent by finding the point sources.
End of explanation
wstep = 8.0
nw = int(1.1 * 800/wstep)
    
gcfcf = create_awterm_convolutionfunction(model, nw=110, wstep=8, oversampling=8, 
                                                    support=60,
                                                    use_aaf=True)
    
cf=gcfcf[1]
print(cf.data.shape)
plt.clf()
plt.imshow(numpy.real(cf.data[0,0,0,0,0,:,:]))
plt.title(str(numpy.max(numpy.abs(cf.data[0,0,0,0,0,:,:]))))
plt.show()
    
cf_clipped = apply_bounding_box_convolutionfunction(cf, fractional_level=1e-3)
print(cf_clipped.data.shape)
gcfcf_clipped=(gcfcf[0], cf_clipped)
    
plt.clf()
plt.imshow(numpy.real(cf_clipped.data[0,0,0,0,0,:,:]))
plt.title(str(numpy.max(numpy.abs(cf_clipped.data[0,0,0,0,0,:,:]))))
plt.show()
Explanation: Make the convolution function
End of explanation
contexts = imaging_contexts().keys()
print(contexts)
print(gcfcf_clipped[1])
contexts = ['2d', 'facets', 'timeslice', 'wstack', 'wprojection']
for context in contexts:
    
    print('Processing context %s' % context)
   
    vtpredict_list =[create_visibility(lowcore, times, frequency, channel_bandwidth=channel_bandwidth,
        weight=1.0, phasecentre=phasecentre, polarisation_frame=PolarisationFrame('stokesI'))]
    model_list = [model]
    vtpredict_list = arlexecute.compute(vtpredict_list, sync=True)
    vtpredict_list = arlexecute.scatter(vtpredict_list)
 
    if context == 'wprojection':       
        future = predict_list_arlexecute_workflow(vtpredict_list, model_list, context='2d', gcfcf=[gcfcf_clipped])
    
    elif context == 'facets':
        future = predict_list_arlexecute_workflow(vtpredict_list, model_list, context=context, facets=8)
        
    elif context == 'timeslice':
        future = predict_list_arlexecute_workflow(vtpredict_list, model_list, context=context, vis_slices=vis_timeslices(
            vtpredict, 'auto'))
    elif context == 'wstack':
        future = predict_list_arlexecute_workflow(vtpredict_list, model_list, context=context, vis_slices=31)
    else:
        future = predict_list_arlexecute_workflow(vtpredict_list, model_list, context=context)
        
    vtpredict_list = arlexecute.compute(future, sync=True)
        
    vtpredict = vtpredict_list[0]
    uvdist = numpy.sqrt(vt.data['uvw'][:, 0] ** 2 + vt.data['uvw'][:, 1] ** 2)
    plt.clf()
    plt.plot(uvdist, numpy.abs(vt.data['vis'][:]), '.', color='r', label="DFT")
    
    plt.plot(uvdist, numpy.abs(vtpredict.data['vis'][:]), '.', color='b', label=context)
    plt.plot(uvdist, numpy.abs(vtpredict.data['vis'][:] - vt.data['vis'][:]), '.', color='g', label="Residual")
    plt.xlabel('uvdist')
    plt.ylabel('Amp Visibility')
    plt.legend()
    plt.show()
Explanation: Predict the visibility using the different approaches.
End of explanation
contexts = ['2d', 'facets', 'timeslice', 'wstack', 'wprojection']
for context in contexts:
    targetimage_list = [create_image_from_visibility(vt, npixel=npixel, cellsize=0.001, nchan=1,
                                               polarisation_frame=PolarisationFrame('stokesI'))]
    
    vt_list = [vt]
    print('Processing context %s' % context)
    if context == 'wprojection':
        future = invert_list_arlexecute_workflow(vt_list, targetimage_list, context='2d', gcfcf=[gcfcf_clipped])
    
    elif context == 'facets':
        future = invert_list_arlexecute_workflow(vt_list, targetimage_list, context=context, facets=8)
        
    elif context == 'timeslice':
        future = invert_list_arlexecute_workflow(vt_list, targetimage_list, context=context, vis_slices=vis_timeslices(vt, 'auto'))
    elif context == 'wstack':
        future = invert_list_arlexecute_workflow(vt_list, targetimage_list, context=context, vis_slices=31)
    else:
        future = invert_list_arlexecute_workflow(vt_list, targetimage_list, context=context)
        
    result = arlexecute.compute(future, sync=True)
    targetimage = result[0][0]
    show_image(targetimage)
    plt.title(context)
    plt.show()
    print("Dirty Image %s" % qa_image(targetimage, context="imaging-fits notebook, using processor %s" % context))
    export_image_to_fits(targetimage, '%s/imaging-fits_dirty_%s.fits' % (results_dir, context))
    comps = find_skycomponents(targetimage, fwhm=1.0, threshold=10.0, npixels=5)
    plt.clf()
    for comp in comps:
        distance = comp.direction.separation(model.phasecentre)
        dft_flux = sum_visibility(vt, comp.direction)[0]
        err = (comp.flux[0, 0] - dft_flux) / dft_flux
        plt.plot(distance, err, '.', color='r')
    plt.ylabel('Fractional error of image vs DFT')
    plt.xlabel('Distance from phasecentre (deg)')
    plt.title(
        "Fractional error in %s recovered flux vs distance from phasecentre" %
        context)
    plt.show()
    checkpositions = True
    if checkpositions:
        plt.clf()
        for i in range(len(comps)):
            ocomp, sep = find_nearest_skycomponent(comps[i].direction, original_comps)
            plt.plot(
                (comps[i].direction.ra.value - ocomp.direction.ra.value) /
                targetimage.wcs.wcs.cdelt[0],
                (comps[i].direction.dec.value - ocomp.direction.dec.value) /
                targetimage.wcs.wcs.cdelt[1],
                '.',
                color='r')
        plt.xlabel('delta RA (pixels)')
        plt.ylabel('delta DEC (pixels)')
        plt.title("%s: Position offsets" % context)
        plt.show()
Explanation: Make the image using the different approaches. We will evaluate the results using a number of plots:
The error in fitted versus the radius. The ideal result is a straightline fitted: flux = DFT flux
The offset in RA versus the offset in DEC. The ideal result is a cluster around 0 pixels.
The sampling in w is set to provide 2% decorrelation at the half power point of the primary beam.
End of explanation
<END_TASK>