Testing MPI version for Imaging and deconvolution demonstration

This script makes a fake data set and then deconvolves it. Finally the full and residual visibility are plotted.


In [1]:
%lsmagic


Out[1]:
Available line magics:
%alias  %alias_magic  %autoawait  %autocall  %automagic  %autosave  %bookmark  %cat  %cd  %clear  %colors  %conda  %config  %connect_info  %cp  %debug  %dhist  %dirs  %doctest_mode  %ed  %edit  %env  %gui  %hist  %history  %killbgscripts  %ldir  %less  %lf  %lk  %ll  %load  %load_ext  %loadpy  %logoff  %logon  %logstart  %logstate  %logstop  %ls  %lsmagic  %lx  %macro  %magic  %man  %matplotlib  %mkdir  %more  %mv  %notebook  %page  %pastebin  %pdb  %pdef  %pdoc  %pfile  %pinfo  %pinfo2  %pip  %popd  %pprint  %precision  %prun  %psearch  %psource  %pushd  %pwd  %pycat  %pylab  %qtconsole  %quickref  %recall  %rehashx  %reload_ext  %rep  %rerun  %reset  %reset_selective  %rm  %rmdir  %run  %save  %sc  %set_env  %store  %sx  %system  %tb  %time  %timeit  %unalias  %unload_ext  %who  %who_ls  %whos  %xdel  %xmode

Available cell magics:
%%!  %%HTML  %%SVG  %%bash  %%capture  %%debug  %%file  %%html  %%javascript  %%js  %%latex  %%markdown  %%perl  %%prun  %%pypy  %%python  %%python2  %%python3  %%ruby  %%script  %%sh  %%svg  %%sx  %%system  %%time  %%timeit  %%writefile

Automagic is ON, % prefix IS NOT needed for line magics.

In [2]:
from mpi4py import MPI
print(MPI.COMM_WORLD.size)


1

In [6]:
%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.testing_support import create_test_image
from wrappers.serial.simulation.configurations import create_named_configuration
from wrappers.serial.imaging.base import create_image_from_visibility
from workflows.mpi.imaging.imaging_mpi import predict_list_mpi_workflow, invert_list_mpi_workflow, deconvolve_list_mpi_workflow,weight_list_mpi_workflow

from data_models.polarisation import PolarisationFrame

import logging

log = logging.getLogger()
log.setLevel(logging.DEBUG)
log.addHandler(logging.StreamHandler(sys.stdout))

In [7]:
pylab.rcParams['figure.figsize'] = (12.0, 12.0)
pylab.rcParams['image.cmap'] = 'rainbow'

Construct LOW core configuration


In [8]:
lowcore = create_named_configuration('LOWBD2', rmax=400.0)


create_configuration_from_file: Maximum radius 400.0 m includes 146 antennas/stations

In [9]:
print(lowcore.xyz)


[[-1697.451395    1786.90252634 -4220.53729446]
 [-1899.245706    1676.22129979 -4000.47217072]
 [-1749.332463    1743.04062531 -4133.32761149]
 [-1980.112908    1741.06827262 -4129.40602467]
 [-1797.563795    1711.65424735 -4070.92274514]
 [-1856.918727    1634.47055507 -3917.46005486]
 [-1854.091115    1685.47915605 -4018.87936901]
 [-1841.933631    1735.45400747 -4118.24330093]
 [-1710.69259     1686.89797425 -4021.70037499]
 [-1998.747741    1692.42716624 -4032.69394958]
 [-1942.270247    1691.73930513 -4031.32629   ]
 [-1654.987395    1654.73474441 -3957.75091061]
 [-1834.982938    1762.89064318 -4172.79497894]
 [-1792.469112    1750.71903274 -4148.59442532]
 [-1882.541622    1709.46473429 -4066.569383  ]
 [-2062.126102    1714.44766716 -4076.47684221]
 [-1789.25535     1675.75003965 -3999.53517422]
 [-1660.667531    1684.05418174 -4016.04612295]
 [-1885.643275    1739.10263634 -4125.49779196]
 [-1776.101351    1727.22766226 -4101.88703409]
 [-1650.104471    1708.87873739 -4065.40425786]
 [-1619.872693    1741.2904703  -4129.84781558]
 [-1744.337431    1765.3622405  -4177.70920321]
 [-1888.296497    1651.32622307 -3950.97382027]
 [-1714.541629    1648.98111578 -3946.31109344]
 [-1936.228874    1638.66040531 -3925.79064478]
 [-1821.529984    1657.8986803  -3964.04169691]
 [-1801.313183    1797.67154979 -4241.94911424]
 [-1891.074282    1770.02558984 -4186.98124131]
 [-1755.697704    1660.18837802 -3968.59425408]
 [-1710.119593    1753.92886297 -4154.97646236]
 [-1753.953802    1631.41687847 -3911.38849475]
 [-1906.906416    1791.85635547 -4230.38688733]
 [-1611.165642    1685.29445809 -4018.51213798]
 [-1928.804836    1749.452451   -4146.07610786]
 [-1996.368561    1769.92596458 -4186.78315854]
 [-1744.723581    1712.97008163 -4073.5389904 ]
 [-1969.811433    1717.41739031 -4082.38147945]
 [-1751.599535    1799.33719149 -4245.26087414]
 [-2060.494309    1689.04551013 -4025.97027481]
 [-1927.833234    1708.77351497 -4065.19504636]
 [-1599.332025    1703.44973006 -4054.60987829]
 [-2064.542652    1745.66110089 -4138.53784723]
 [-1768.764508    1697.33174237 -4042.44561376]
 [-1909.609464    1810.43139608 -4267.31924467]
 [-2030.573942    1668.48188608 -3985.08405939]
 [-1652.570845    1774.18801544 -4195.25730342]
 [-2033.937181    1728.8725953  -4105.1576194 ]
 [-1818.802024    1619.00348139 -3886.7072019 ]
 [-1613.13376     1764.06487547 -4175.12967999]
 [-1793.029651    1643.42505856 -3935.26410328]
 [-1858.488238    1786.22194151 -4219.18410214]
 [-1970.920056    1666.16476328 -3980.47697351]
 [-1701.150957    1704.66929775 -4057.03471873]
 [-1584.596058    1723.99208605 -4095.45380687]
 [-2020.608791    1755.07735019 -4157.25997502]
 [-1784.422252    1775.50552883 -4197.8768872 ]
 [-1946.580174    1778.18980905 -4203.21398436]
 [-1837.162815    1702.94153005 -4053.59943505]
 [-1678.92867     1734.62398497 -4116.5929849 ]
 [-1991.448268    1791.10189136 -4228.88680242]
 [-1953.642975    1801.4119678  -4249.38610765]
 [-1928.630446    1729.68134973 -4106.76564859]
 [-1983.264387    1642.49709036 -3933.4190439 ]
 [-1899.295532    1621.75660357 -3892.18117607]
 [-1843.689989    1813.21586096 -4272.85553685]
 [-2014.330746    1712.49658298 -4072.59754313]
 [-1709.471859    1629.13949389 -3906.86041953]
 [-1693.72692     1666.1882703  -3980.52371202]
 [-1861.191285    1665.57876638 -3979.31184838]
 [-1656.469711    1755.16634096 -4157.43691347]
 [-1620.146735    1665.83062765 -3979.81261877]
 [-1785.119813    1815.25593636 -4276.91177528]
 [-2042.233169    1770.55225955 -4188.02840747]
 [-1721.679169    1728.63080866 -4104.67688017]
 [-2068.628364    1597.71337726 -3844.37654164]
 [-2184.622723    1667.47947876 -3983.09099428]
 [-1998.722828    1825.14064745 -4296.56533566]
 [-2038.757823    1812.87612817 -4272.18005338]
 [-1582.951808    1628.10854178 -3904.81059943]
 [-1939.517374    1837.44826271 -4321.03630475]
 [-2120.771017    1757.98214899 -4163.03552456]
 [-1764.952838    1841.09521248 -4328.28745726]
 [-2016.996424    1611.99670387 -3872.77577551]
 [-2196.506166    1688.45335667 -4024.7929087 ]
 [-2039.654686    1836.38596835 -4318.92416753]
 [-2169.077087    1707.21925181 -4062.10473804]
 [-1982.753673    1573.70876991 -3796.64869248]
 [-2118.354468    1685.43550022 -4018.79256905]
 [-2154.091991    1635.51158091 -3919.52990435]
 [-1715.326385    1604.50523337 -3857.88064438]
 [-1688.906278    1804.54344174 -4255.61235055]
 [-1913.309026    1576.24137397 -3801.68421515]
 [-1961.615096    1606.72720837 -3862.29854993]
 [-1623.049085    1595.01062733 -3839.00272155]
 [-1886.951201    1599.59058261 -3848.10894902]
 [-2123.312131    1618.11413126 -3884.93892601]
 [-2105.163099    1734.98386658 -4117.30852982]
 [-1661.576851    1586.36171513 -3821.80627377]
 [-1765.799876    1584.24160388 -3817.59090177]
 [-2060.843089    1787.05923982 -4220.84888453]
 [-1820.271883    1840.9536105  -4328.00591305]
 [-1968.204553    1854.46400284 -4354.86833819]
 [-2107.841233    1642.9980141  -3934.41501989]
 [-2179.602778    1775.49377554 -4197.8535184 ]
 [-2035.79319     1650.2269893  -3948.7882372 ]
 [-1859.098604    1564.81695182 -3778.96928008]
 [-1552.707574    1642.67059445 -3933.76401837]
 [-2123.411782    1661.25794819 -3970.72085767]
 [-1817.357077    1583.45635648 -3816.02961109]
 [-1611.477053    1643.87672974 -3936.16215147]
 [-2084.722083    1817.90999359 -4282.18878071]
 [-2136.241914    1816.46207166 -4279.30990838]
 [-2101.874599    1801.11029412 -4248.78629632]
 [-2173.100517    1755.00626927 -4157.11864634]
 [-1883.724983    1854.49702464 -4354.93399474]
 [-1551.374735    1684.53775502 -4017.0076014 ]
 [-1717.929781    1577.6881762  -3804.56086119]
 [-2182.841452    1736.15586037 -4119.63878009]
 [-2059.460424    1620.59748232 -3889.87651999]
 [-1800.790013    1558.88478179 -3767.17447285]
 [-1998.598264    1594.93059148 -3838.84358799]
 [-2103.792891    1782.09925413 -4210.98705069]
 [-1859.385102    1832.79890606 -4311.79208803]
 [-1616.783496    1613.44462535 -3875.65464694]
 [-1838.50811     1601.09615262 -3851.10244179]
 [-1579.675765    1603.8145739  -3856.50742088]
 [-1575.652335    1661.53947306 -3971.28060755]
 [-2221.319392    1704.54336711 -4056.78433353]
 [-2234.161981    1760.02558215 -4167.09843916]
 [-1886.041881    1506.16128566 -3662.34546922]
 [-2154.664987    1603.95505663 -3856.7867397 ]
 [-1720.881957    1531.28359693 -3712.29562527]
 [-2252.659793    1660.70721206 -3969.62584074]
 [-2267.769453    1641.24953819 -3930.93856251]
 [-2050.292485    1535.39621013 -3720.47264644]
 [-2221.842563    1615.02183614 -3878.79058152]
 [-2238.210324    1725.04878324 -4097.55481537]
 [-1830.311773    1527.98365622 -3705.73442349]
 [-2215.539605    1644.71682642 -3937.83249778]
 [-2244.338892    1685.21889971 -4018.36190686]
 [-1925.279664    1516.90288447 -3683.70276119]
 [-2255.400209    1596.34773012 -3841.66125455]
 [-2278.070928    1708.74944855 -4065.1471956 ]
 [-2293.267784    1752.89231459 -4152.91551532]
 [-2309.884674    1734.08388287 -4115.51911139]]

We create the visibility. This just makes the uvw, time, antenna1, antenna2, weight columns in a table


In [10]:
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(lowcore, times, frequency, channel_bandwidth=channel_bandwidth,
                       weight=1.0, phasecentre=phasecentre, polarisation_frame=PolarisationFrame('stokesI'))
print(vt)


create_visibility: 10585 rows, 0.001 GB
Visibility:
	Number of visibilities: 10585
	Number of channels: 1
	Frequency: [1.e+08]
	Number of polarisations: 1
	Visibility shape: (10585, 1)
	Polarisation Frame: stokesI
	Phasecentre: <SkyCoord (ICRS): (ra, dec) in deg
    (15., -45.)>
	Configuration: LOWBD2

Plot the synthesized uv coverage.


In [11]:
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.show()


findfont: Matching :family=sans-serif:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to DejaVu Sans ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/DejaVuSans.ttf') with score of 0.050000

Read the venerable test image, constructing an image


In [12]:
m31image = create_test_image(frequency=frequency, cellsize=0.0005)
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)
print("My image:")
print(m31image)
print(m31image.data)
#print(m31image.wcs)
#print(m31image.polarisation_frame)
sumwt = numpy.zeros([m31image.nchan, m31image.npol])
fig=show_image(m31image)


import_image_from_fits: created >f4 image of shape (256, 256), size 0.000 (GB)
import_image_from_fits: Max, min in /Users/montsefarreras/Documents/repo/tmp/algorithm-reference-library/data/models/M31.MOD = 1.006458, 0.000000
replicate_image: replicating shape (256, 256) to (1, 1, 256, 256)
My image:
Image:
	Shape: (1, 1, 256, 256)
	WCS: WCS Keywords

Number of WCS axes: 4
CTYPE : 'RA---SIN'  'DEC--SIN'  'STOKES'  'FREQ'  
CRVAL : 15.0  -45.0  1.0  100000000.0  
CRPIX : 128.0  128.0  1.0  1.0  
PC1_1 PC1_2 PC1_3 PC1_4  : 1.0  0.0  0.0  0.0  
PC2_1 PC2_2 PC2_3 PC2_4  : 0.0  1.0  0.0  0.0  
PC3_1 PC3_2 PC3_3 PC3_4  : 0.0  0.0  1.0  0.0  
PC4_1 PC4_2 PC4_3 PC4_4  : 0.0  0.0  0.0  1.0  
CDELT : -0.02864788975654116  0.02864788975654116  1.0  100000.0  
NAXIS : 0  0
	Polarisation frame: stokesI

[[[[0. 0. 0. ... 0. 0. 0.]
   [0. 0. 0. ... 0. 0. 0.]
   [0. 0. 0. ... 0. 0. 0.]
   ...
   [0. 0. 0. ... 0. 0. 0.]
   [0. 0. 0. ... 0. 0. 0.]
   [0. 0. 0. ... 0. 0. 0.]]]]
findfont: Matching :family=STIXGeneral:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to STIXGeneral ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/STIXGeneral.ttf') with score of 0.050000
findfont: Matching :family=STIXGeneral:style=italic:variant=normal:weight=normal:stretch=normal:size=10.0 to STIXGeneral ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/STIXGeneralItalic.ttf') with score of 0.050000
findfont: Matching :family=STIXGeneral:style=normal:variant=normal:weight=bold:stretch=normal:size=10.0 to STIXGeneral ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/STIXGeneralBol.ttf') with score of 0.000000
findfont: Matching :family=STIXNonUnicode:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to STIXNonUnicode ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/STIXNonUni.ttf') with score of 0.050000
findfont: Matching :family=STIXNonUnicode:style=italic:variant=normal:weight=normal:stretch=normal:size=10.0 to STIXNonUnicode ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/STIXNonUniIta.ttf') with score of 0.050000
findfont: Matching :family=STIXNonUnicode:style=normal:variant=normal:weight=bold:stretch=normal:size=10.0 to STIXNonUnicode ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/STIXNonUniBol.ttf') with score of 0.000000
findfont: Matching :family=STIXSizeOneSym:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to STIXSizeOneSym ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/STIXSizOneSymReg.ttf') with score of 0.050000
findfont: Matching :family=STIXSizeTwoSym:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to STIXSizeTwoSym ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/STIXSizTwoSymReg.ttf') with score of 0.050000
findfont: Matching :family=STIXSizeThreeSym:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to STIXSizeThreeSym ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/STIXSizThreeSymReg.ttf') with score of 0.050000
findfont: Matching :family=STIXSizeFourSym:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to STIXSizeFourSym ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/STIXSizFourSymReg.ttf') with score of 0.050000
findfont: Matching :family=STIXSizeFiveSym:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to STIXSizeFiveSym ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/STIXSizFiveSymReg.ttf') with score of 0.050000
findfont: Matching :family=cmsy10:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to cmsy10 ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/cmsy10.ttf') with score of 0.050000
findfont: Matching :family=cmr10:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to cmr10 ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/cmr10.ttf') with score of 0.050000
findfont: Matching :family=cmtt10:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to cmtt10 ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/cmtt10.ttf') with score of 0.050000
findfont: Matching :family=cmmi10:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to cmmi10 ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/cmmi10.ttf') with score of 0.050000
findfont: Matching :family=cmb10:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to cmb10 ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/cmb10.ttf') with score of 0.050000
findfont: Matching :family=cmss10:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to cmss10 ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/cmss10.ttf') with score of 0.050000
findfont: Matching :family=cmex10:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to cmex10 ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/cmex10.ttf') with score of 0.050000
findfont: Matching :family=DejaVu Sans:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to DejaVu Sans ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/DejaVuSans.ttf') with score of 0.050000
findfont: Matching :family=DejaVu Sans:style=italic:variant=normal:weight=normal:stretch=normal:size=10.0 to DejaVu Sans ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/DejaVuSans-Oblique.ttf') with score of 0.150000
findfont: Matching :family=DejaVu Sans:style=normal:variant=normal:weight=bold:stretch=normal:size=10.0 to DejaVu Sans ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/DejaVuSans-Bold.ttf') with score of 0.000000
findfont: Matching :family=DejaVu Sans Mono:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to DejaVu Sans Mono ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/DejaVuSansMono.ttf') with score of 0.050000
findfont: Matching :family=DejaVu Sans Display:style=normal:variant=normal:weight=normal:stretch=normal:size=10.0 to DejaVu Sans Display ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/DejaVuSansDisplay.ttf') with score of 0.050000

In [13]:
vt = predict_list_mpi_workflow([vt], [m31image], context='2d',use_serial_predict=True)
print(vt)


0: In predict_list_mpi_workflow: 1 elements in vis_list
shift_vis_from_image: shifting phasecentre from image phase centre <SkyCoord (ICRS): (ra, dec) in deg
    (14.95950601, -44.97134495)> to visibility phasecentre <SkyCoord (ICRS): (ra, dec) in deg
    (15., -45.)>
[<data_models.memory_data_models.Visibility object at 0x11a3b5780>]
[<data_models.memory_data_models.Visibility object at 0x1196284a8>]

In [14]:
# To check that we got the prediction right, plot the amplitude of the visibility.
vt=vt[0]
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()


Make the dirty image and point spread function


In [15]:
model = create_image_from_visibility(vt, cellsize=0.001, npixel=256)
dirty, sumwt = invert_list_mpi_workflow([vt], [model], context='2d')[0]
psf, sumwt = invert_list_mpi_workflow([vt], [model], context='2d', dopsf=True)[0]

show_image(dirty)
show_image(psf)
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))


create_image_from_visibility: Parsing parameters to get definition of WCS
create_image_from_visibility: Defining single channel Image at <SkyCoord (ICRS): (ra, dec) in deg
    (15., -45.)>, starting frequency 100000000.0 Hz, and bandwidth 999999.99999 Hz
create_image_from_visibility: uvmax = 253.011682 wavelengths
create_image_from_visibility: Critical cellsize = 0.001976 radians, 0.113228 degrees
create_image_from_visibility: Cellsize          = 0.001000 radians, 0.057296 degrees
0: In invert_list_mpi_workflow: 1 elements in vis_list 1 in model
[(<data_models.memory_data_models.Image object at 0x105569550>, array([[10585.]]))]
0: In invert_list_mpi_workflow: 1 elements in vis_list 1 in model
[(<data_models.memory_data_models.Image object at 0x11910e860>, array([[10585.]]))]
Max, min in dirty image = 33.418636, -7.539496, sumwt = 10585.000000
Max, min in PSF         = 0.999209, -0.024851, sumwt = 10585.000000

Deconvolve using clean


In [18]:
result=deconvolve_list_mpi_workflow([(dirty,sumwt)], [(psf,sumwt)], [model])


threshold_list: using refchan 0 , sub_image 0, peak = 33.418636,
threshold_list : Global peak = 33.418636, sub-image clean threshold will be 3.341864
deconvolve_cube facet 0: PSF support = +/- 128 pixels
deconvolve_cube facet 0: PSF shape (1, 1, 256, 256)
deconvolve_cube facet 0: Multi-scale clean of each polarisation and channel separately
deconvolve_cube facet 0: Processing pol 0, channel 0
msclean facet 0: Peak of PSF = 0.9992088831836738 at (128, 128)
msclean facet 0: Peak of Dirty = 33.418636 Jy/beam at (148, 124) 
msclean facet 0: Coupling matrix =
 [[1.         0.96508872 0.59971887 0.12718816]
 [0.96508872 0.93218613 0.58573945 0.12654143]
 [0.59971887 0.58573945 0.42406774 0.11731114]
 [0.12718816 0.12654143 0.11731114 0.06923701]]
msclean facet 0: Max abs in dirty Image = 33.445095 Jy/beam
msclean facet 0: Start of minor cycle
msclean facet 0: This minor cycle will stop at 100 iterations or peak < 3.341864 (Jy/beam)
msclean facet 0: Timing for setup: 0.158 (s) for dirty shape (256, 256), PSF shape (256, 256) , scales [0, 3, 10, 30]
msclean facet 0: Minor cycle 0, peak [23.04812494 23.11827695 23.66695334 17.98439773] at [137, 128, 3]
msclean facet 0: Minor cycle 10, peak [5.45804726 5.45115393 5.29728783 3.2850186 ] at [166, 138, 3]
msclean facet 0: At iteration 12, absolute value of peak 2.846833 is below stopping threshold 3.341864
msclean facet 0: End of minor cycle
msclean facet 0: Timing for clean: 0.058 (s) for dirty shape (256, 256), PSF shape (256, 256) , scales [0, 3, 10, 30], 13 iterations, time per clean 4.492 (ms)

In [16]:
comp, residual = deconvolve_cube(dirty, psf, niter=1000, threshold=0.001, fracthresh=0.01, 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')


deconvolve_cube : Cleaning inner quarter of each sky plane
deconvolve_cube : PSF support = +/- 128 pixels
deconvolve_cube : PSF shape (1, 1, 256, 256)
deconvolve_cube : Multi-scale clean of each polarisation and channel separately
deconvolve_cube : Processing pol 0, channel 0
msclean : Peak of PSF = 0.9992088831836738 at (128, 128)
msclean : Peak of Dirty = 33.418636 Jy/beam at (148, 124) 
msclean : Coupling matrix =
 [[1.         0.96508872 0.59971887 0.12718816]
 [0.96508872 0.93218613 0.58573945 0.12654143]
 [0.59971887 0.58573945 0.42406774 0.11731114]
 [0.12718816 0.12654143 0.11731114 0.06923701]]
msclean : Max abs in dirty Image = 33.445095 Jy/beam
msclean : Start of minor cycle
msclean : This minor cycle will stop at 1000 iterations or peak < 0.334451 (Jy/beam)
msclean : Timing for setup: 0.219 (s) for dirty shape (256, 256), PSF shape (256, 256) , scales [0, 3, 10, 30]
msclean : Minor cycle 0, peak [23.04812494 23.11827695 23.66695334 17.98439773] at [137, 128, 3]
msclean : Minor cycle 100, peak [1.25715144 1.26278188 1.24960649 0.43864615] at [163, 120, 3]
msclean : At iteration 140, absolute value of peak 0.299661 is below stopping threshold 0.334451
msclean : End of minor cycle
msclean : Timing for clean: 0.362 (s) for dirty shape (256, 256), PSF shape (256, 256) , scales [0, 3, 10, 30], 141 iterations, time per clean 2.564 (ms)
restore_cube: psfwidth = Parameter('x_stddev', value=1.880165227887132)
Out[16]:
Text(0.5,1,'Restored')
findfont: Matching :family=sans-serif:style=normal:variant=normal:weight=normal:stretch=normal:size=12.0 to DejaVu Sans ('/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/DejaVuSans.ttf') with score of 0.050000

Predict the visibility of the model


In [ ]:
vtmodel = create_visibility(lowcore, 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]

Now we will plot the original visibility and the residual visibility.


In [ ]:
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()

In [ ]: