<jupyter_start><jupyter_text>Pipeline processing using serial workflows.This notebook demonstrates the continuum imaging and ICAL pipelines. These are based on RASCIL functions wrapped     up as SDP workflows using the serial class.<jupyter_code>%matplotlib inline

import os
import sys

sys.path.append(os.path.join('..', '..'))

results_dir = '/tmp/'

from matplotlib import pylab
pylab.rcParams['figure.figsize'] = (12.0, 12.0)
pylab.rcParams['image.cmap'] = 'rainbow'

import numpy

from astropy.coordinates import SkyCoord
from astropy import units as u

from matplotlib import pyplot as plt

from rascil.data_models.polarisation import PolarisationFrame

from rascil.processing_components.visibility import  convert_blockvisibility_to_visibility
from rascil.processing_components.image import show_image, export_image_to_fits, qa_image
from rascil.processing_components.simulation import create_low_test_image_from_gleam
from rascil.processing_components.imaging import  create_image_from_visibility, advise_wide_field

from rascil.workflows.serial.imaging.imaging_serial import invert_list_serial_workflow, \
    predict_list_serial_workflow, deconvolve_list_serial_workflow
from rascil.workflows.serial.simulation.simulation_serial import simulate_list_serial_workflow, \
    corrupt_list_serial_workflow
from rascil.workflows.serial.pipelines.pipeline_serial import continuum_imaging_list_serial_workflow

import pprint

pp = pprint.PrettyPrinter()

import logging

def init_logging():
    log = logging.getLogger()
    logging.basicConfig(filename='%s/imaging-pipeline.log' % results_dir,
                        filemode='a',
                        format='%(asctime)s,%(msecs)d %(name)s %(levelname)s %(message)s',
                        datefmt='%H:%M:%S',
                        level=logging.INFO)
log = logging.getLogger()
mpl_logger = logging.getLogger("matplotlib") 
mpl_logger.setLevel(logging.WARNING) 

logging.info("Starting imaging-pipeline")
pylab.rcParams['figure.figsize'] = (12.0, 12.0)
pylab.rcParams['image.cmap'] = 'Greys'<jupyter_output><empty_output><jupyter_text>We make the visibility. The parameter rmax determines the distance of the furthest antenna/stations used. All over parameters are determined from this number.<jupyter_code>nfreqwin=7
ntimes=5
rmax=300.0
frequency=numpy.linspace(1e8,1.2e8,nfreqwin)
channel_bandwidth=numpy.array(nfreqwin*[frequency[1]-frequency[0]])
times = numpy.linspace(-numpy.pi/3.0, numpy.pi/3.0, ntimes)
phasecentre=SkyCoord(ra=+30.0 * u.deg, dec=-60.0 * u.deg, frame='icrs', equinox='J2000')

bvis_list=simulate_list_serial_workflow('LOWBD2',
                                         frequency=frequency, 
                                         channel_bandwidth=channel_bandwidth,
                                         times=times,
                                         phasecentre=phasecentre,
                                         order='frequency',
                                        rmax=rmax, format='blockvis')
vis_list = [convert_blockvisibility_to_visibility(bv) for bv in bvis_list]

print('%d elements in vis_list' % len(vis_list))
wprojection_planes=1
advice_low=advise_wide_field(vis_list[0], guard_band_image=8.0, delA=0.02,
                             wprojection_planes=wprojection_planes)

advice_high=advise_wide_field(vis_list[-1], guard_band_image=8.0, delA=0.02,
                              wprojection_planes=wprojection_planes)

vis_slices = advice_low['vis_slices']
npixel=advice_high['npixels2']
cellsize=min(advice_low['cellsize'], advice_high['cellsize'])<jupyter_output><empty_output><jupyter_text>Now make a graph to fill with a model drawn from GLEAM<jupyter_code>gleam_model = [create_low_test_image_from_gleam(npixel=npixel,
                                                               frequency=[frequency[f]],
                                                               channel_bandwidth=[channel_bandwidth[f]],
                                                               cellsize=cellsize,
                                                               phasecentre=phasecentre,
                                                               polarisation_frame=PolarisationFrame("stokesI"),
                                                               flux_limit=1.0,
                                                               applybeam=True)
                     for f, freq in enumerate(frequency)]
log.info('About to make GLEAM model')
log.info('About to run predict to get predicted visibility')
predicted_vislist = predict_list_serial_workflow(vis_list, gleam_model,  
                                                context='wstack', vis_slices=vis_slices)
log.info('About to run corrupt to get corrupted visibility')
corrupted_vislist = corrupt_list_serial_workflow(predicted_vislist, phase_error=1.0)<jupyter_output><empty_output><jupyter_text>Get the LSM. This is currently blank.<jupyter_code>model_list = [create_image_from_visibility(vis_list[f],
                                                     npixel=npixel,
                                                     frequency=[frequency[f]],
                                                     channel_bandwidth=[channel_bandwidth[f]],
                                                     cellsize=cellsize,
                                                     phasecentre=phasecentre,
                                                     polarisation_frame=PolarisationFrame("stokesI"))
               for f, freq in enumerate(frequency)]
dirty_list = invert_list_serial_workflow(predicted_vislist, model_list, 
                                  context='wstack',
                                  vis_slices=vis_slices, dopsf=False)
psf_list = invert_list_serial_workflow(predicted_vislist, model_list, 
                                context='wstack',
                                vis_slices=vis_slices, dopsf=True)<jupyter_output><empty_output><jupyter_text>Create and execute graphs to make the dirty image and PSF<jupyter_code>log.info('About to run invert to get dirty image')
dirty = dirty_list[0][0]
show_image(dirty, cm='Greys', vmax=1.0, vmin=-0.1)
plt.show()

log.info('About to run invert to get PSF')
psf = psf_list[0][0]
show_image(psf, cm='Greys', vmax=0.1, vmin=-0.01)
plt.show()<jupyter_output><empty_output><jupyter_text>Now deconvolve using msclean<jupyter_code>log.info('About to run deconvolve')

deconvolved = \
    deconvolve_list_serial_workflow(dirty_list, psf_list, model_imagelist=model_list, 
                            deconvolve_facets=8, deconvolve_overlap=16, deconvolve_taper='tukey',
                            scales=[0, 3, 10],
                            algorithm='msclean', niter=1000, 
                            fractional_threshold=0.1,
                            threshold=0.1, gain=0.1, psf_support=64)
    
centre=nfreqwin // 2
show_image(deconvolved[centre], cm='Greys', vmax=0.1, vmin=-0.01)
plt.show()
continuum_imaging_list = \
    continuum_imaging_list_serial_workflow(predicted_vislist, 
                                            model_imagelist=model_list, 
                                            context='wstack', vis_slices=vis_slices, 
                                            scales=[0, 3, 10], algorithm='mmclean', 
                                            nmoment=3, niter=1000, 
                                            fractional_threshold=0.1,
                                            threshold=0.1, nmajor=5, gain=0.25,
                                            deconvolve_facets = 8, deconvolve_overlap=16, 
                                            deconvolve_taper='tukey', psf_support=64)

log.info('About to run continuum imaging')

centre=nfreqwin // 2
deconvolved = continuum_imaging_list[0][centre]
residual = continuum_imaging_list[1][centre]
restored = continuum_imaging_list[2][centre]

f=show_image(deconvolved, title='Clean image - no selfcal', cm='Greys', 
             vmax=0.1, vmin=-0.01)
print(qa_image(deconvolved, context='Clean image - no selfcal'))

plt.show()

f=show_image(restored, title='Restored clean image - no selfcal', 
             cm='Greys', vmax=1.0, vmin=-0.1)
print(qa_image(restored, context='Restored clean image - no selfcal'))
plt.show()
export_image_to_fits(restored, '%s/imaging-dask_continuum_imaging_restored.fits' 
                     %(results_dir))

f=show_image(residual[0], title='Residual clean image - no selfcal', cm='Greys', 
             vmax=0.1, vmin=-0.01)
print(qa_image(residual[0], context='Residual clean image - no selfcal'))
plt.show()
export_image_to_fits(residual[0], '%s/imaging-dask_continuum_imaging_residual.fits' 
                     %(results_dir))
for chan in range(nfreqwin):
    residual = continuum_imaging_list[1][chan]
    show_image(residual[0], title='Channel %d' % chan, cm='Greys', 
               vmax=0.1, vmin=-0.01)
    plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0942215267687493, 'fraction_numerical': 0.016318318989605455, 'mean_word_length': 2.3067257945306725, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 57, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '362516', 'license': 'Apache-2.0', 'n_tokens_mistral': 2718, 'n_tokens_neox': 2549, 'n_words': 450, 'path': 'examples/notebooks/imaging-pipelines_serial.ipynb', 'repo_name': 'SKA-ScienceDataProcessor/rascil', 'size': '8947'}