<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Imaging and deconvolution demonstration
This script makes a fake data set and then deconvolves it. Finally the full and residual visibility are plotted.
Step1:  Construct LOW core configuration
Step2:  We create the visibility. This just makes the uvw, time, antenna1, antenna2, weight columns in a table
Step3:  Plot the synthesized uv coverage.
Step4:  Read the venerable test image, constructing an image
Step5:  Make the dirty image and point spread function
Step6:  Deconvolve using clean
Step7:  Predict the visibility of the model
Step8:  Now we will plot the original visibility and the residual visibility.
<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'] = (8.0, 8.0)
pylab.rcParams['image.cmap'] = 'rainbow'
import numpy
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.wcs.utils import pixel_to_skycoord
from matplotlib import pyplot as plt
from processing_components.image.iterators import image_raster_iter
from wrappers.serial.visibility.base import create_visibility
from wrappers.serial.skycomponent.operations import create_skycomponent
from wrappers.serial.image.operations import show_image, export_image_to_fits
from wrappers.serial.image.deconvolution import deconvolve_cube, restore_cube
from wrappers.serial.visibility.iterators import vis_timeslice_iter
from wrappers.serial.simulation.configurations import create_named_configuration
from wrappers.serial.simulation.testing_support import create_test_image
from wrappers.serial.imaging.base import create_image_from_visibility
from wrappers.serial.imaging.base import advise_wide_field
from workflows.serial.imaging.imaging_serial import invert_list_serial_workflow, predict_list_serial_workflow
from data_models.polarisation import PolarisationFrame
import logging
log = logging.getLogger()
log.setLevel(logging.DEBUG)
log.addHandler(logging.StreamHandler(sys.stdout))
mpl_logger = logging.getLogger("matplotlib") 
mpl_logger.setLevel(logging.WARNING) 
pylab.rcParams['figure.figsize'] = (12.0, 12.0)
pylab.rcParams['image.cmap'] = 'rainbow'
Explanation: Imaging and deconvolution demonstration
This script makes a fake data set and then deconvolves it. Finally the full and residual visibility are plotted.
End of explanation
lowr3 = create_named_configuration('LOWBD2', rmax=750.0)
print(lowr3.xyz)
Explanation: Construct LOW core configuration
End of explanation
times = numpy.zeros([1])
frequency = numpy.array([1e8])
channel_bandwidth = numpy.array([1e6])
phasecentre = SkyCoord(ra=+15.0 * u.deg, dec=-45.0 * u.deg, frame='icrs', equinox='J2000')
vt = create_visibility(lowr3, times, frequency, channel_bandwidth=channel_bandwidth,
                       weight=1.0, phasecentre=phasecentre, polarisation_frame=PolarisationFrame('stokesI'))
advice = advise_wide_field(vt, guard_band_image=3.0, delA=0.1, facets=1, wprojection_planes=1, 
                           oversampling_synthesised_beam=4.0)
cellsize = advice['cellsize']
Explanation: We create the visibility. This just makes the uvw, time, antenna1, antenna2, weight columns in a table
End of explanation
plt.clf()
plt.plot(vt.data['uvw'][:,0], vt.data['uvw'][:,1], '.', color='b')
plt.plot(-vt.data['uvw'][:,0], -vt.data['uvw'][:,1], '.', color='b')
plt.xlim([-400.0, 400.0])
plt.ylim([-400.0, 400.0])
plt.show()
Explanation: Plot the synthesized uv coverage.
End of explanation
m31image = create_test_image(frequency=frequency, cellsize=cellsize)
nchan, npol, ny, nx = m31image.data.shape
m31image.wcs.wcs.crval[0] = vt.phasecentre.ra.deg
m31image.wcs.wcs.crval[1] = vt.phasecentre.dec.deg
m31image.wcs.wcs.crpix[0] = float(nx // 2)
m31image.wcs.wcs.crpix[1] = float(ny // 2)
fig=show_image(m31image)
vt = predict_list_serial_workflow([vt], [m31image], context='2d')[0]
# To check that we got the prediction right, plot the amplitude of the visibility.
uvdist=numpy.sqrt(vt.data['uvw'][:,0]**2+vt.data['uvw'][:,1]**2)
plt.clf()
plt.plot(uvdist, numpy.abs(vt.data['vis']), '.')
plt.xlabel('uvdist')
plt.ylabel('Amp Visibility')
plt.show()
Explanation: Read the venerable test image, constructing an image
End of explanation
model = create_image_from_visibility(vt, cellsize=cellsize, npixel=512)
dirty, sumwt = invert_list_serial_workflow([vt], [model], context='2d')[0]
psf, sumwt = invert_list_serial_workflow([vt], [model], context='2d', dopsf=True)[0]
show_image(dirty)
print("Max, min in dirty image = %.6f, %.6f, sumwt = %f" % (dirty.data.max(), dirty.data.min(), sumwt))
print("Max, min in PSF         = %.6f, %.6f, sumwt = %f" % (psf.data.max(), psf.data.min(), sumwt))
export_image_to_fits(dirty, '%s/imaging_dirty.fits'%(results_dir))
export_image_to_fits(psf, '%s/imaging_psf.fits'%(results_dir))
Explanation: Make the dirty image and point spread function
End of explanation
comp, residual = deconvolve_cube(dirty, psf, niter=10000, threshold=0.001, fractional_threshold=0.001,
                                 window_shape='quarter', gain=0.7, scales=[0, 3, 10, 30])
restored = restore_cube(comp, psf, residual)
# Show the results
fig=show_image(comp)
plt.title('Solution')
fig=show_image(residual)
plt.title('Residual')
fig=show_image(restored)
plt.title('Restored')
Explanation: Deconvolve using clean
End of explanation
vtmodel = create_visibility(lowr3, times, frequency, channel_bandwidth=channel_bandwidth,
                            weight=1.0, phasecentre=phasecentre, 
                            polarisation_frame=PolarisationFrame('stokesI'))
vtmodel=predict_list_serial_workflow([vtmodel], [comp], context='2d')[0]
Explanation: Predict the visibility of the model
End of explanation
uvdist=numpy.sqrt(vt.data['uvw'][:,0]**2+vt.data['uvw'][:,1]**2)
plt.clf()
plt.plot(uvdist, numpy.abs(vt.data['vis'][:]-vtmodel.data['vis'][:]), '.', color='r', 
         label='Residual')
plt.plot(uvdist, numpy.abs(vt.data['vis'][:]), '.', color='b', label='Original')
plt.xlabel('uvdist')
plt.ylabel('Amp Visibility')
plt.legend()
plt.show()
Explanation: Now we will plot the original visibility and the residual visibility.
End of explanation
<END_TASK>