Adam Ginsburg's blog - Articles by Adam Ginsburg (Adam.G.Ginsburg@gmail.com)

Articles by Adam Ginsburg (Adam.G.Ginsburg@gmail.com)

CASA table locking with MPI

A bunch of CASA MPI runs have race conditioned themselves into a corner, apparently:

2021-11-07 23:33:10     WARN    MPICommandServer::command_request_handler_service::SynthesisImagerVi2::CubeMajorCycle::MPIServer-5 (file src/code/synthesis/ImagerObjects/CubeMajorCycleAlgorithm.cc, line 336) Exception for chan range [104, 129] ---   Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G351.77_B6_spw1_12M_sio.contcube.model/table.lock
##################################
#############################
Exception: Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G351.77_B6_spw1_12M_sio.contcube.model/table.lock
2021-11-07 23:33:10     WARN    MPICommandServer::command_request_handler_service::SynthesisImagerVi2::CubeMajorCycle::MPIServer-10 (file src/code/synthesis/ImagerObjects/CubeMajorCycleAlgorithm.cc, line 336)        Exception for chan range [234, 258] ---   Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G351.77_B6_spw1_12M_sio.contcube.model/table.lock
##################################
#############################
Exception: Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G351.77_B6_spw1_12M_sio.contcube.model/table.lock
2021-11-07 23:33:10     WARN    MPICommandServer::command_request_handler_service::SynthesisImagerVi2::CubeMajorCycle::MPIServer-9 (file src/code/synthesis/ImagerObjects/CubeMajorCycleAlgorithm.cc, line 336) Exception for chan range [208, 233] ---   Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G351.77_B6_spw1_12M_sio.contcube.model/table.lock
##################################
#############################
Exception: Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G351.77_B6_spw1_12M_sio.contcube.model/table.lock
2021-11-07 23:33:10     WARN    MPICommandServer::command_request_handler_service::SynthesisImagerVi2::CubeMajorCycle::MPIServer-2 (file src/code/synthesis/ImagerObjects/CubeMajorCycleAlgorithm.cc, line 336) Exception for chan range [26, 51] ---   Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G351.77_B6_spw1_12M_sio.contcube.model/table.lock
##################################
#############################
Exception: Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G351.77_B6_spw1_12M_sio.contcube.model/table.lock
2021-11-07 23:33:10     WARN    MPICommandServer::command_request_handler_service::SynthesisImagerVi2::CubeMajorCycle::MPIServer-15 (file src/code/synthesis/ImagerObjects/CubeMajorCycleAlgorithm.cc, line 336)        Exception for chan range [359, 383] ---   Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G351.77_B6_spw1_12M_sio.contcube.model/table.lock
##################################
#############################
Exception: Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G351.77_B6_spw1_12M_sio.contcube.model/table.lock

This shows up in ~10 runs that I started around the same time, but I don't usually see it.

This looks like a serious bug having to do with using a startmodel: there is a line earlier copying the startmodel to the model:

2021-11-07 23:27:44     INFO    SIImageStore::setModelImageOne  Copying input model /blue/adamginsburg/adamginsburg/almaimf/workdir//G351.77_B6_spw1_12M_sio.contcube.model to /blue/adamginsburg/adamginsburg/almaimf/workdir/G351.77_B6_spw1_12M_sio.model

but MPI seems to have gotten mixed up and is sill trying to write to the contcube.model

CASA doesn't like VaryingResolutionSpectralCubes

I get this from trying to run makemask:

 makemask(mode='copy',
         inpimage=lineimagename+".image",
         inpmask=lineimagename+".image:mask0",
         output=lineimagename+".mask")

 2021-11-09 02:27:10     INFO    makemask::::casa        ##########################################
 2021-11-09 02:27:10     INFO    makemask::::casa        ##### Begin Task: makemask           #####
 2021-11-09 02:27:10     INFO    makemask::::casa        makemask( mode='copy', inpimage='/blue/adamginsburg/adamginsburg/almaimf/workdir/G353.41_B6_spw5_12M_spw5.image', inpmask='/blue/adamginsburg/adamginsburg/almaimf/workdir/G353.41_B6_spw5_12M_spw5.image:mask0', output='/blue/adamginsburg/adamginsburg/almaimf/workdir/G353.41_B6_spw5_12M_spw5.mask', overwrite=False, inpfreqs=[], outfreqs=[] )
 2021-11-09 02:27:11     INFO    makemask::ImageFactory::createImage     Created Paged image '__tmp_outputmask_22199' of shape [1080, 900, 1, 480] with float valued pixels.
 2021-11-09 02:27:11     INFO    makemask::::casa        Summing all T/F mask in inpmask and normalized to 1 for mask
 2021-11-09 02:27:51     SEVERE  image::regrid   Exception Reported: Exception: An image with multiple beams cannot be regridded along the spectral axis. You may wish to convolve all channels to a common resolution and retry.
 2021-11-09 02:27:51     SEVERE  image::regrid+  ... thrown by std::shared_ptr<casa6core::ImageInterface<T> > casa::ImageRegridder<T>::regrid() const [with T = float] at File: src/code/imageanalysis/ImageAnalysis/ImageRegridder.tcc, line: 102
 2021-11-09 02:27:52     INFO            deleted table /blue/adamginsburg/adamginsburg/almaimf/workdir/G353.41_spw5_12M_B6/_tmp_copy___tmp_frominmask_22199
 2021-11-09 02:27:56     SEVERE  image::calc (file src/tools/image/image_cmpt.cc, line 553)      Exception Reported: ImageExprParse: '__tmp_fromTFmask_22199' is an unknown lattice or image or it is an unqualified region
 2021-11-09 02:27:56     SEVERE  image::calc (file src/tools/image/image_cmpt.cc, line 553)+     Scanned so far: iif ("__tmp_fromTFmask_22199"
 2021-11-09 02:27:56     SEVERE  makemask::::casa        Task makemask raised an exception of class RuntimeError with the following message: *** Error (2), in mode copy: *** ImageExprParse: '__tmp_fromTFmask_22199' is an unknown lattice or image or it is an unqualified region
 2021-11-09 02:27:56     SEVERE  makemask::::casa+       Scanned so far: iif ("__tmp_fromTFmask_22199"
 2021-11-09 02:27:56     INFO    makemask::::casa        Task makemask complete. Start time: 2021-11-08 21:27:09.614534 End time: 2021-11-08 21:27:56.399654
 2021-11-09 02:27:56     INFO    makemask::::casa        ##### End Task: makemask             #####
 2021-11-09 02:27:56     INFO    makemask::::casa        ##########################################


2021-11-09 02:27:51     SEVERE  image::regrid   Exception Reported: Exception: An image with multiple beams cannot be regridded along the spectral axis. You may wish to convolve all channels to a common
 resolution and retry.
2021-11-09 02:27:51     SEVERE  image::regrid+  ... thrown by std::shared_ptr<casa6core::ImageInterface<T> > casa::ImageRegridder<T>::regrid() const [with T = float] at File: src/code/imageanalysis/Imag
eAnalysis/ImageRegridder.tcc, line: 102
2021-11-09 02:27:56     SEVERE  image::calc (file src/tools/image/image_cmpt.cc, line 553)      Exception Reported: ImageExprParse: '__tmp_fromTFmask_22199' is an unknown lattice or image or it is an un
qualified region
2021-11-09 02:27:56     SEVERE  image::calc (file src/tools/image/image_cmpt.cc, line 553)+     Scanned so far: iif ("__tmp_fromTFmask_22199"
2021-11-09 02:27:56     SEVERE  makemask::::casa        Task makemask raised an exception of class RuntimeError with the following message: *** Error (2), in mode copy: *** ImageExprParse: '__tmp_fromTF
mask_22199' is an unknown lattice or image or it is an unqualified region
2021-11-09 02:27:56     SEVERE  makemask::::casa+       Scanned so far: iif ("__tmp_fromTFmask_22199"

I'm pretty sure the complaint is caused by the top line: there are multiple beams in the image. This is _absurd_, though, because the makemask command should not care at all about beams, there is no spectral regridding happening.

This might be new to CASA 6.4, not sure.

EDIT 2021-11-11 15:51: The workaround is easy and maybe better:

ia.open(lineimagename+".image")
shape = ia.shape()
csys = ia.coordsys().torecord()
ia.close()

ia.fromshape(outfile=lineimagename+".mask", shape=shape, csys=csys, type='f')

Walkthrough of PSF characterization

2021/11/06
For ALMA-IMF, we want to characterize the PSF shapes and the recovered beam sizes.

We showed the beams and a few features of the PSFs, but they were only described tersely in the paper. This document goes through what we did in a little more detail
Front matter (imports)
```
%matplotlib inline
import pylab as pl

# for display within notebook
pl.rcParams['figure.facecolor'] = 'w'

from spectral_cube import SpectralCube
from astropy import units as u

# Debuggy loading: this approach needed b/c the code isn't properly packaged
import imp
import sys
sys.path.append('/orange/adamginsburg/ALMA_IMF/reduction/analysis')
import imstats
imp.reload(imstats)
from imstats import get_psf_secondpeak
```
N2H+

I chose to do this only for an N$_2$H+ example; the N$_2$H+ line cubes seem to have the worst-behaved beams in our samples.
```
basename = '/orange/adamginsburg/ALMAIMF_Images/G327.29/B3/linecubes_12m/G327.29_B3_spw0_12M_n2hp'
img = SpectralCube.read(f'{basename}.image')
psf = SpectralCube.read(f'{basename}.psf', format='casa_image')
```
The code get_psf_secondpeak does all of the work, we’re just going to talk about what gets plotted.
```
psffn = f'{basename}.psf'
(psf_secondpeak, psf_secondpeak_loc, psf_sidelobe1_fraction, epsilon, firstnull, r_sidelobe,
 (rr, pixscale, cutout, beam, fullbeam, view, bmfit_residual)) = \
        get_psf_secondpeak(psffn, show_image=True, min_radial_extent=2.5*u.arcsec,
                   max_radial_extent=5*u.arcsec,
                           specslice=slice(50,51),
                  )
```
```
INFO: first_sidelobe_ind=44, first_min_ind = 43 [imstats]
```
The plot is showing:
- in red contours, the 10th, 50th, and 90th percentile of the Gaussian synthesized beam
- in grayscale, the residual of the dirty beam minus the synthesized beam
- in green dashed (faint, outermost ellipse), the location of the first sidelobe peak. The first sidelobe is defined as the first peak outside of the first null.
- in blue dashed (inner ellipse), the location of the first null. This may not be a true null: it is the first minimum of the radial profile
Next, we show the radial profile. The radial profile is an elliptical profile, which is why the elliptical synthesized beam has a single profile.

The X-axis is calculated as:

\begin{equation*} r = \left(C^2 (x \cos \theta + y \sin \theta)^2 + (x \sin \theta - y \cos \theta)^2\right)^{1/2} \end{equation*}

where $C = \theta_{min} / \theta_{maj}$. So the radial axis is scaled to the major axis of the beam; the true beam is smaller along one axis than shown in these plots.
```
# radial profile
pl.figure(figsize=(12,6))
ax2=pl.subplot(1,2,1)
ax1=pl.subplot(1,2,2)


ax2.set_title('G327 N2H+')
rr_inds = rr.ravel().argsort()
sorted_synth = (fullbeam.array.ravel()/fullbeam.array.max())[rr_inds]
pixscale = pixscale.to(u.arcsec)
ax2.plot(pixscale.value*rr.ravel()[rr_inds], sorted_synth, '-', label='Synth')
ax2.plot(pixscale.value*rr.ravel(),
         cutout.value.ravel()/cutout.max().value, '.', label='Dirty', alpha=0.75, markersize=2)
ax2.axvline(firstnull.value, linestyle='--', color='b')
ax2.set_xlim(0, firstnull.value*2)# rr[view].max())
ax2.text(firstnull.value*1.4, 0.9, f'$\epsilon={epsilon:0.2f}$', fontsize=14)
ax2.text(firstnull.value*1.4, 0.85, f'$\\theta_{{maj}}={beam.major.value:0.2f}$"', fontsize=14)
ax2.text(firstnull.value*1.4, 0.8, f'$\\theta_{{min}}={beam.minor.value:0.2f}$"', fontsize=14)

ax2.axvline(psf_secondpeak_loc, linestyle=':', color='k')
ax2.axvline(r_sidelobe.value, linestyle=':', color='g')

ax2.plot([0, beam.major.value/2, beam.major.value/2],
         [0.5, 0.5, 0],
         color='k', alpha=0.5)


ax1.set_title('G327 N2H+')
ax1.plot(pixscale.value*rr.ravel()[rr_inds], sorted_synth, '-', label='Synth')
ax1.plot(pixscale.value*rr.ravel(),
         cutout.value.ravel()/cutout.max().value, '.', label='Dirty', alpha=0.75, markersize=2)
ax1.axvline(firstnull.value, linestyle='--', color='b')
ax1.axvline(psf_secondpeak_loc, linestyle=':', color='k')
ax1.axvline(r_sidelobe.value, linestyle=':', color='g')

ax1.set_xlim(0, firstnull.value*2)# rr[view].max())
ax1.set_ylim(-0.035, 0.2)

ax1.set_xlabel("Radius along PSF (\")")
ax2.set_xlabel("Radius along PSF (\")")
ax1.set_ylabel("Normalized PSF profile")
```
```
Text(0, 0.5, 'Normalized PSF profile')
```
The right plot is a vertical zoom-in of the left plot.

The light blue line shows the synthesized beam.

The orange dots are the individual pixel values from the dirty beam.

The vertical dotted black line shows the peak of the residual - it’s not used for anything, it just highlights where the “shelf” is most prominent (and it is very prominent in this example).

The blue vertical dashed line shows the location of the first null.

The green vertical dotted line shows the location of the first sidelobe.

$\epsilon$ is the ratio of the integral of the synthesized beam to the integral of the dirty beam out to the first null. It is the correction factor to apply to the residuals (multiply the residuals by this number) to convert the residuals from Jy/dirtybeam to Jy/synthesizedbeam.

$\theta_{min}$ and $\theta_{maj}$ are the FWHM widths of the beam along the major and minor axes. The grey right-angle line is drawn at $y=0.5$ out to $\theta = \theta_{maj}/2$ (the half-width half-max)

CASA MPI problems

There are several CASA MPI errors I encounter regularly.

2021-10-17 15:56:06     SEVERE  task_tclean::SynthesisImagerVi2::runCubeGridding (file src/code/synthesis/ImagerObjects/SynthesisImagerVi2.cc, line 1579)       remainder rank 7 failed
master 1 init 1
2021-10-17 15:59:21     SEVERE  tclean::::casa  Task tclean raised an exception of class RuntimeError with the following message: Error in making PSF : One or more  of the cube section failed in de/grid
ding. Return values for the sections: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

 Traceback (most recent call last):
   File "/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-12/lib/py/lib/python3.6/site-packages/casashell/private/init_system.py", line 238, in __evprop__
     exec(stmt)
   File "<string>", line 1, in <module>
   File "/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-12/lib/py/lib/python3.6/site-packages/casashell/private/init_system.py", line 175, in execfile
     newglob = run_path( filename, init_globals=globals )
   File "/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-12/lib/py/lib/python3.6/runpy.py", line 263, in run_path
     pkg_name=pkg_name, script_name=fname)
   File "/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-12/lib/py/lib/python3.6/runpy.py", line 96, in _run_module_code
     mod_name, mod_spec, pkg_name, script_name)
   File "/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-12/lib/py/lib/python3.6/runpy.py", line 85, in _run_code
     exec(code, run_globals)
   File "/orange/adamginsburg/ALMA_IMF/reduction/reduction/line_imaging.py", line 935, in <module>
     **impars
   File "/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-12/lib/py/lib/python3.6/site-packages/casatasks/tclean.py", line 1660, in __call__
     task_result = _tclean_t( _pc.document['vis'], _pc.document['selectdata'], _pc.document['field'], _pc.document['spw'], _pc.document['timerange'], _pc.document['uvrange'], _pc.document['antenna'], _pc.document['scan'], _pc.document['observation'], _pc.document['intent'], _pc.document['datacolumn'], _pc.document['imagename'], _pc.document['imsize'], _pc.document['cell'], _pc.document['phasecenter'], _pc.document['stokes'], _pc.document['projection'], _pc.document['startmodel'], _pc.document['specmode'], _pc.document['reffreq'], _pc.document['nchan'], _pc.document['start'], _pc.document['width'], _pc.document['outframe'], _pc.document['veltype'], _pc.document['restfreq'], _pc.document['interpolation'], _pc.document['perchanweightdensity'], _pc.document['gridder'], _pc.document['facets'], _pc.document['psfphasecenter'], _pc.document['wprojplanes'], _pc.document['vptable'], _pc.document['mosweight'], _pc.document['aterm'], _pc.document['psterm'], _pc.document['wbawp'], _pc.document['conjbeams'], _pc.document['cfcache'], _pc.document['usepointing'], _pc.document['computepastep'], _pc.document['rotatepastep'], _pc.document['pointingoffsetsigdev'], _pc.document['pblimit'], _pc.document['normtype'], _pc.document['deconvolver'], _pc.document['scales'], _pc.document['nterms'], _pc.document['smallscalebias'], _pc.document['restoration'], _pc.document['restoringbeam'], _pc.document['pbcor'], _pc.document['outlierfile'], _pc.document['weighting'], _pc.document['robust'], _pc.document['noise'], _pc.document['npixels'], _pc.document['uvtaper'], _pc.document['niter'], _pc.document['gain'], _pc.document['threshold'], _pc.document['nsigma'], _pc.document['cycleniter'], _pc.document['cyclefactor'], _pc.document['minpsffraction'], _pc.document['maxpsffraction'], _pc.document['interactive'], _pc.document['usemask'], _pc.document['mask'], _pc.document['pbmask'], _pc.document['sidelobethreshold'], _pc.document['noisethreshold'], _pc.document['lownoisethreshold'], _pc.document['negativethreshold'], _pc.document['smoothfactor'], _pc.document['minbeamfrac'], _pc.document['cutthreshold'], _pc.document['growiterations'], _pc.document['dogrowprune'], _pc.document['minpercentchange'], _pc.document['verbose'], _pc.document['fastnoise'], _pc.document['restart'], _pc.document['savemodel'], _pc.document['calcres'], _pc.document['calcpsf'], _pc.document['psfcutoff'], _pc.document['parallel'] )
   File "/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-12/lib/py/lib/python3.6/site-packages/casatasks/private/task_tclean.py", line 364, in tclean
     imager.makePSF()
   File "/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-12/lib/py/lib/python3.6/site-packages/casatasks/private/imagerhelpers/imager_base.py", line 344, in makePSF
     self.makePSFCore()
   File "/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-12/lib/py/lib/python3.6/site-packages/casatasks/private/imagerhelpers/imager_base.py", line 496, in makePSFCore
     self.SItool.makepsf()
   File "/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-12/lib/py/lib/python3.6/site-packages/casatools/synthesisimager.py", line 70, in makepsf
     return self._swigobj.makepsf()
   File "/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-12/lib/py/lib/python3.6/site-packages/casatools/__casac__/synthesisimager.py", line 322, in makepsf
     return _synthesisimager.synthesisimager_makepsf(self)
 RuntimeError: Error in making PSF : One or more  of the cube section failed in de/gridding. Return values for the sections: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

That one appears to be caused by:

2021-10-17 15:53:34     WARN    MPICommandServer::command_request_handler_service::SynthesisImagerVi2::CubeMajorCycle::MPIServer-7 (file src/code/synthesis/ImagerObjects/CubeMajorCycleAlgorithm.cc, line
 A) Exception for chan range [1921, 1924] ---   Error in making PSF : Interpolate1D::operator() data has repeated x values
##################################
#############################
Exception: Error in making PSF : Interpolate1D::operator() data has repeated x values

This happened after a lot of apparently successful cleaning, as an .image was created. However, the failure caused CASA to exit.

This tclean command:

2021-10-17 15:19:43     INFO    tclean::::casa  tclean( vis='/blue/adamginsburg/adamginsburg/almaimf/workdir/G008.67_B3_spw2_12M.concat.ms', selectdata=True, field='G008.67', spw='', timerange='', uvrange='', antenna='', scan='', observation='', intent='', datacolumn='corrected', imagename='/blue/adamginsburg/adamginsburg/almaimf/workdir/G008.67_B3_spw2_12M_spw2', imsize=[2880, 2250], cell=['0.08arcsec', '0.08arcsec'], phasecenter='ICRS 271.5877979623041deg -21.620789662367244deg', stokes='I', projection='SIN', startmodel='', specmode='cube', reffreq='', nchan=-1, start='', width='', outframe='LSRK', veltype='radio', restfreq=[], interpolation='linear', perchanweightdensity=True, gridder='mosaic', facets=1, psfphasecenter='', wprojplanes=1, vptable='', mosweight=True, aterm=True, psterm=False, wbawp=True, conjbeams=False, cfcache='', usepointing=False, computepastep=360.0, rotatepastep=360.0, pointingoffsetsigdev=[], pblimit=0.05, normtype='flatnoise', deconvolver='multiscale', scales=[0, 4, 8, 16, 32], nterms=2, smallscalebias=0.5, restoration=True, restoringbeam='', pbcor=False, outlierfile='', weighting='briggsbwtaper', robust=0.0, noise='1.0Jy', npixels=0, uvtaper=[''], niter=5000000, gain=0.1, threshold='0.0168Jy', nsigma=0.0, cycleniter=-1, cyclefactor=2.0, minpsffraction=0.05, maxpsffraction=0.8, interactive=0, usemask='user', mask='', pbmask=0.1, sidelobethreshold=3.0, noisethreshold=5.0, lownoisethreshold=1.5, negativethreshold=0.0, smoothfactor=1.0, minbeamfrac=0.3, cutthreshold=0.01, growiterations=75, dogrowprune=True, minpercentchange=-1.0, verbose=False, fastnoise=True, restart=True, savemodel='none', calcres=False, calcpsf=True, psfcutoff=0.35, parallel=True )

using concatenated data caused the failure. It was run as an MPI job with 32 cores and 128 GB RAM. The same issue recurred for the next MS (spw3), but not for spw0 or spw1 - yet, they're still running.

This: "Exception for chan range [1921, 1924] --- Error in making PSF : Interpolate1D::operator() data has repeated x values" suggests to me that there's a problem with the gridder thinking there are more channels than there actually are; tclean's behavior is not consistent between operating on concatenated data sets and on lists of data sets. This issue may be a manifestation of tclean misunderstanding the grid when operating with MPI enabled. This could also prove not to be an MPI error, but one that is only encountered if MPI is enabled b/c MPI gets to channel 1921 while non-MPI never gets there before the job is canceled for time.

W51 spw2 had similar:

WARN    MPICommandServer::command_request_handler_service::SynthesisImagerVi2::CubeMajorCycle::MPIServer-21 (file src/code/synthesis/ImagerObjects/CubeMajorCycleAlgorithm.cc, line 336)        Exception for chan range [1269, 1277] ---   FilebufIO::readBlock - incorrect number of bytes read for file /blue/adamginsburg/adamginsburg/almaimf/workdir/W51-E_B3_spw2_12M_spw2.sumwt/table.f0
Exception: FilebufIO::readBlock - incorrect number of bytes read for file /blue/adamginsburg/adamginsburg/almaimf/workdir/W51-E_B3_spw2_12M_spw2.sumwt/table.f0
WARN    MPICommandServer::command_request_handler_service::SynthesisImagerVi2::CubeMajorCycle::MPIServer-9 (file src/code/synthesis/ImagerObjects/CubeMajorCycleAlgorithm.cc, line 336) Exception for chan range [1314, 1322] ---   FilebufIO::readBlock - incorrect number of bytes read for file /blue/adamginsburg/adamginsburg/almaimf/workdir/W51-E_B3_spw2_12M_spw2.sumwt/table.f0
Exception: FilebufIO::readBlock - incorrect number of bytes read for file /blue/adamginsburg/adamginsburg/almaimf/workdir/W51-E_B3_spw2_12M_spw2.sumwt/table.f0

I have not found solutions to any of these; the workaround appears to be just to not use MPI.

EDIT: here's another one

2021-10-22 22:36:53     WARN    MPICommandServer::command_request_handler_service::SynthesisImagerVi2::CubeMajorCycle::MPIServer-2 (file src/code/synthesis/ImagerObjects/CubeMajorCycleAlgorithm.cc, line 336) Exception for chan range [55, 109] ---   Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G327.29_B6_spw1_12M_sio.contcube.model/table.lock
##################################
#############################
Exception: Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G327.29_B6_spw1_12M_sio.contcube.model/table.lock
2021-10-22 22:36:54     WARN    MPICommandServer::command_request_handler_service::SynthesisImagerVi2::CubeMajorCycle::MPIServer-4 (file src/code/synthesis/ImagerObjects/CubeMajorCycleAlgorithm.cc, line 336) Exception for chan range [165, 219] ---   Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G327.29_B6_spw1_12M_sio.contcube.model/table.lock
##################################
#############################
Exception: Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G327.29_B6_spw1_12M_sio.contcube.model/table.lock
2021-10-22 22:36:57     WARN    MPICommandServer::command_request_handler_service::SynthesisImagerVi2::CubeMajorCycle::MPIServer-1 (file src/code/synthesis/ImagerObjects/CubeMajorCycleAlgorithm.cc, line 336) Exception for chan range [0, 54] ---   Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G327.29_B6_spw1_12M_sio.contcube.model/table.lock
##################################
#############################
Exception: Setting masked pixels to zero for input startmodel : Error (Resource deadlock avoided) when acquiring lock on /blue/adamginsburg/adamginsburg/almaimf/workdir/G327.29_B6_spw1_12M_sio.contcube.model/table.lock

Confetti artifacts in clean

We have noted a “confetti” problem in some of the ALMA-IMF cube images. It is called confetti because one of its manifestations is individual bright pixels in the model image.

This notebook shows where confetti is - and is not a problem.

https://github.com/ALMA-IMF/notebooks/blob/master/ConfettiQA.ipynb

Preliminaries: skip this

from spectral_cube import SpectralCube
import warnings
warnings.simplefilter('ignore')

%matplotlib inline
import pylab as pl
from astropy.visualization import simple_norm

from astropy import units as u

G008.67 Continuum

This is the G008.67 continuum image in Band 3.

It has “confetti” pieces in the model image, but there is no problem: no artificial point sources are introduced.

contmodel = SpectralCube.read('/orange/adamginsburg/ALMA_IMF/2017.1.01355.L/imaging_results/G008.67_B3_uid___A001_X1296_X1c1_continuum_merged_12M_robust-0.5_selfcal5_finaliter.model.tt0', format='casa_image')
contimage = SpectralCube.read('/orange/adamginsburg/ALMA_IMF/2017.1.01355.L/imaging_results/G008.67_B3_uid___A001_X1296_X1c1_continuum_merged_12M_robust-0.5_selfcal5_finaliter.image.tt0', format='casa_image')
contresid = SpectralCube.read('/orange/adamginsburg/ALMA_IMF/2017.1.01355.L/imaging_results/G008.67_B3_uid___A001_X1296_X1c1_continuum_merged_12M_robust-0.5_selfcal5_finaliter.residual.tt0', format='casa_image')

# zoom in
contslc = (0, slice(950,1050), slice(1900,2050))

pl.figure(figsize=(15,15))
pl.imshow(contmodel[contslc].value, norm=simple_norm(contmodel[contslc].value, max_percent=99.995, stretch='asinh'), interpolation='none', origin='lower')
_=pl.title("model")

pl.figure(figsize=(15,15))
im = contimage[0][contslc[1:]].value
pl.imshow(im, norm=simple_norm(im, stretch='log', max_percent=99.995), interpolation='none',  origin='lower')
_=pl.title("image")

pl.figure(figsize=(15,15))
im = contresid[0][contslc[1:]].value
pl.imshow(im, norm=simple_norm(im, stretch='linear'), interpolation='none',  origin='lower')
_=pl.title("Residual")

G008.67 Cube

The G008.67 cube does have a problem!

cube = SpectralCube.read('/orange/adamginsburg/ALMA_IMF/2017.1.01355.L/imaging_results/G008.67_B3_spw0_12M_n2hp.image', format='casa_image').with_spectral_unit(u.km/u.s, velocity_convention='radio')
resid = SpectralCube.read('/orange/adamginsburg/ALMA_IMF/2017.1.01355.L/imaging_results/G008.67_B3_spw0_12M_n2hp.residual', format='casa_image').with_spectral_unit(u.km/u.s, velocity_convention='radio')
model = SpectralCube.read('/orange/adamginsburg/ALMA_IMF/2017.1.01355.L/imaging_results/G008.67_B3_spw0_12M_n2hp.model', format='casa_image').with_spectral_unit(u.km/u.s, velocity_convention='radio')

xslc = slice(1620,1654)
yslc = slice(709,743)
cubeslice = (slice(None), yslc, xslc)
cube1slice = (98, yslc, xslc)

Confetti: a false (?) point source

In the image below, we see a new source introduced in the 35.140 km/s panel.

There is no comparable source in either of the adjacent panels!

That means either there is a genuinely pointlike source with linewidth much narrower than 0.77 km/s, or this source is fake.

specpix = np.arange(94,94+9)
panels = cube[cubeslice].plot_channel_maps(3,3,specpix, fig_smallest_dim_inches=15, norm=simple_norm(cube[cube1slice].value))

It is caused by one hot (bright) pixel surrounded by many faint pixels

panels = model[cubeslice].plot_channel_maps(3,3,specpix, fig_smallest_dim_inches=15, norm=simple_norm(model[cube1slice].value, stretch='log', max_percent=99.93))

There is no sign of it in the residual - it is apparently canceled out between positive and negative components added to the model

panels = resid[cubeslice].plot_channel_maps(3,3,specpix, fig_smallest_dim_inches=15)

The manifestation in the spectra is a very sharp peak at just one pixel.

pl.figure(figsize=(14,6))
xx,yy=1638,725
cube[:,yy,xx].quicklook()
resid[:,yy,xx].quicklook()
modspec = model[:,yy,xx]
pl.plot(cube.spectral_axis, modspec*cube.pixels_per_beam, label='model')

Pixels to either side look fine

pl.figure(figsize=(14,14))
ii = 1
for dx in (-1,0,1):
    for dy in (-1,0,1):
        ax = pl.subplot(3, 3, ii)
        ii +=1
        cube[:,yy+dy,xx+dx].quicklook()
        resid[:,yy+dy,xx+dx].quicklook()
        modspec = model[:,yy+dy,xx+dx]
        ax.plot(cube.spectral_axis, modspec*cube.pixels_per_beam, label='model')
        ax.set_title(f"dx={dx} dy={dy}")

When is confetti not a problem?

Here is an example you might think is “confetti”, but is a real source

xc, yc = 1812, 905
xslc = slice(xc-12,xc+12)
yslc = slice(yc-12,yc+12)
cubeslice = (slice(94,94+9), yslc, xslc)
cube1slice = (94, yslc, xslc)

specpix = np.arange(9)
panels = cube[cubeslice].plot_channel_maps(3,3,specpix, fig_smallest_dim_inches=15, norm=simple_norm(cube[cube1slice].value, min_cut=cube[cubeslice].min().value, max_cut=cube[cubeslice].max().value, stretch='linear'))

panels = model[cubeslice].plot_channel_maps(3,3,specpix, fig_smallest_dim_inches=15, norm=simple_norm(model[cube1slice].value, stretch='asinh', max_percent=99.9, min_cut=model[cubeslice].min().value, max_cut=model[cubeslice].max().value, ))

panels = resid[cubeslice].plot_channel_maps(3,3,specpix, fig_smallest_dim_inches=15)

pl.figure(figsize=(14,14))
ii = 1
for dx in (-1,0,1):
    for dy in (-1,0,1):
        ax = pl.subplot(3, 3, ii)
        ii +=1
        cube[:,yc+dy,xc+dx].quicklook()
        resid[:,yc+dy,xc+dx].quicklook()
        modspec = model[:,yc+dy,xc+dx]
        ax.plot(cube.spectral_axis, modspec*cube.pixels_per_beam, label='model')
        ax.set_title(f"dx={dx} dy={dy}")

Conclusion

“Confetti” is a problem, but it is not trivial to find in the model data alone. It exhibits no features in the residual.

Key features to look out for are:

Narrow-line features (single channel)
Extremely bright pointlike peaks in the model (note that this is not independent evidence of a problem! All point sources occupy single pixels in the model image by definition!)
Bright pointlike peaks sitting on top of slightly extended negative backgrounds in the model (this is how the model is “canceled out” in the residual)

casa's plotms in a notebook

CASA's plotms requires an X11 frontend to work properly, so you either need to run it with xvfb or with X.

I want to just make UV data plots in jupyter notebooks.

This is harder than it sounds, particularly for concatenated data sets.

First step is that you need to assemble the metadata:

msmd.open(scvis)
scsum = msmd.summary()
msmd.close()

The resulting scsum is a nested dictionary containing a ton of information. We need only some of it.

To get out some flattened information, such as the "data description IDs" associated with data that are actually in my concatenated MS (in my case, 8 of 11 observation blocks in the metadata are not present in the data), this tool gets and flattens the resulting information:

def walk_summary(summary, key):
    if isinstance(summary, dict):
        try:
            if key in summary:
                return summary[key]
            else:
                result = [walk_summary(summary[kk], key) for kk in summary]
                try:
                    result = np.unique([x for x in result if x is not None])
                    if len(result) > 0:
                        return result
                except:
                    return [x for x in result if x is not None]
        except (KeyError,TypeError):
            return

Then we grab the usable observation and data description IDs:

obsids = [key for key in scsum if 'observationID' in key and scsum[key] != {}]
ddids = walk_summary(scsum, 'data description IDs')

Now we have to load the actual data. I needed to create new versions of the mstool each time because I was getting weird errors where ms.close() followed by ms.open(fn) was not properly resetting, which resulted in ms.getdata returning the same data for two different MSes - BAD!

ms6 = mstool()
ms6.open(scvis)
scdata_all = {}
for ddid in ddids:
    ms6.selectinit(ddid, reset=True)
    ms6.selectinit(ddid)
    scdata_all[ddid] = ms6.getdata(['amplitude', 'weight', 'uvdist', 'axis_info', 'flag', 'corrected_amplitude'], ifraxis=True)
ms6.close()

Then plotting is a hassle, but at least is not insane - it's just normal array manipulation:

pl.figure(figsize=(12,12))
colors = itertools.cycle(pl.rcParams["axes.prop_cycle"].by_key()["color"])
for key in scdata_all: # the keys are SPW numbers from the multi-MS file
    weight = scdata_all[key]['weight']
    uvd = scdata_all[key]['uvdist']
    pl.semilogy(uvd.ravel(), weight.reshape(weight.shape[0], uvd.size).T, ',', color=next(colors), label=key, alpha=0.5)
pl.legend(loc='best')
pl.xlabel("UV Distance (m)")
pl.ylabel("Weight")

That is basically equivalent to plotms(vis, xaxis='uvdist', yaxis='weight', coloraxis='spw'). It's just slightly more flexible, since you can also make waterfall plots, etc.

Selecting along various axes, like antenna, gets tricky. You'd rather not go back to ms.getdata because it is _slow_.

The "identifying baselines" step that you can do in plotms can be achieved through normal boolean array selection, with a result that's somewhat easier to read than CASA's logger:

highwt = (scdata_all[56]['weight'] > 1e5)
scdata_all[56]['axis_info']['ifr_axis']['ifr_name'][highwt.any(axis=(0,2))]
# array(['DA41-DV12', 'DA42-DV12', 'DA43-DV12', 'DA44-DV12', 'DA45-DV12',
#  'DA46-DV12', 'DA47-DV12', 'DA48-DV12', 'DA49-DV12', 'DA50-DV12',
#  'DA51-DV12', 'DA52-DV12', 'DA53-DV12', 'DA54-DV12', 'DA55-DV12',
#  'DA56-DV12', 'DA57-DV12', 'DA58-DV12', 'DA60-DV12', 'DA61-DV12',
#  'DA62-DV12', 'DA63-DV12', 'DA64-DV12', 'DA65-DV12', 'DV01-DV12',
#  'DV02-DV12', 'DV03-DV12', 'DV04-DV12', 'DV05-DV12', 'DV06-DV12',
#  'DV07-DV12', 'DV08-DV12', 'DV09-DV12', 'DV10-DV12', 'DV11-DV12',
#  'DV12-DV14', 'DV12-DV16', 'DV12-DV17', 'DV12-DV19', 'DV12-DV20',
#  'DV12-DV21', 'DV12-DV22', 'DV12-DV23', 'DV12-DV24', 'DV12-DV25'],
# dtype='<U16')

in this case, the output shows that DV12 is clearly the antenna that's overweighted.

Handling flags is tricky. The 'weight' array has shape [2,nbaseline,ntime], but I don't know if that 2 refers to frequency or polarization, since I have two of each and there is no information about this in the axis_info dictionary. So... I'm doing the conservative thing and ignoring any row of frequency or polarization if _any_ of the data are flagged.

pl.figure(figsize=(12,12))
colors = itertools.cycle(pl.rcParams["axes.prop_cycle"].by_key()["color"])
for key in scdata_all:
    weight = scdata_all[key]['weight']
    uvd = scdata_all[key]['uvdist']
    flag = scdata_all[key]['flag']
    okflag = ~np.any(flag, axis=(0,1))
    pl.semilogy(uvd[okflag], weight.reshape(weight.shape[0], uvd.size).T[okflag.ravel(), :], ',', color=next(colors), label=key, alpha=0.5)
pl.legend(loc='best')
_=pl.xlabel("UV Distance (m)")
_=pl.ylabel("Weight")

And, indeed, doing this revealed that the super-high-weight antenna was already flagged out.

This has been a demo of how to do some stuff with CASA to replicate plotms

Errors

Along the way, I hit a ton of errors I didn't understand, so I need a record of what they mean.

"Data shape varies, selecting first data desc id only"

This is the worst. It means that you have a multi-ms file (mms) that has been produced by concat. It is therefore doing something explicitly wrong and returning only the first data description ID.

Data description IDs are super counterintuitive. They refer, afaict, to SPW numbers. If you have an MMS with unmerged SPWs (same frequency, but you didn't merge them because of tolerances or something), each SPW will correspond to a single observation ID.

The solution is to explicitly ms.selectinit, i.e.:

ms.open(vis)
ms.selectinit(datadescid) # note: "reset=True" appears to _override_ the selection
ms.select(...) # use this as needed
data = ms.getdata(...)
ms.close()

but datadescid is _not_ the observation ID. It's just... some number you need to get from the metadata.

Mismatched data shapes when comparing two MSes

I got into this whole thing because I wanted to compare the data in two MSes to determine where we went wrong (we did go wrong).

Unfortunately, somewhere along the way, the autocorrelations got dropped. I don't know how it happened because I never did it (but I wish I had earlier), but that means that now the data shapes are mismatched.

The solution was to use ms.select and explicitly downselect to the baselines from MS #1 in MS #2 using MS #1's axis_info: {'ifr_number': scdata['axis_info']['ifr_axis']['ifr_number']}:

ms3 = mstool()
ms3.open(newvis)
ms3.selectinit(80) # 80? what the hell kind of datadescid is that?!
ms3.select({'ifr_number': scdata['axis_info']['ifr_axis']['ifr_number']})
newdata = ms3.getdata(['amplitude', 'uvdist', 'axis_info', 'corrected_amplitude'], ifraxis=True)
ms3.close()
newdata.keys()

That downselected the baselines to match the baselines present in scdata, which was MS 1.

The wrong data are showing up when I open a new MS

I was comparing several MSes, and sometimes the data looked the same for MSes I was dead certain were different. I was right, and the problem was apparently reusing the ms tool.

I didn't show the import statement above, but it was this:

from casatools import ms, msmetadata
mstool = ms
ms = ms()
msmd = msmetadata()

If I tried to reuse ms as normal, like so:

ms.open(vis1)
ms.selectinit(80)
data1 = ms.getdata(...)
ms.close()

ms.open(vis2)
ms.selectinit(80)
data2 = ms.getdata(...)
ms.close()

then I can't say exactly what data1 and data2 are, but they weren't the same.

So, instead, I created a new mstool instance each time:

ms = mstool()
ms.open(vis1)
ms.selectinit(80)
data1 = ms.getdata(...)
ms.close()

ms = mstool()
ms.open(vis2)
ms.selectinit(80)
data2 = ms.getdata(...)
ms.close()

Put all together...

Here's some code to loop over a list of MSes and plot both amp and weight vs. uvdist for the field W43-MM3.

Of course this crashed immediately because it exceeded the available memory.

for vis in mses:
    msmd.open(vis)
    vissum = msmd.summary()
    msmd.close()

    vobsids = [key for key in vissum if 'observationID' in key and vissum[key] != {}]
    vddids = walk_summary(vissum, 'data description IDs')
    fids = np.where(vissum['fields'] == 'W43-MM3')[0]
    print(f"obsids: {vobsids}, ddids: {vddids}")

    ms = mstool()
    ms.open(vis)
    msdata_all = {}
    for ddid in vddids:
        ms.selectinit(ddid, reset=True)
        ms.selectinit(ddid)
        ms.select({'field_id': fids})
        msdata_all[ddid] = ms.getdata(['amplitude', 'weight', 'uvdist', 'axis_info', 'flag', 'corrected_amplitude'], ifraxis=True)
    ms.close()

    colors = itertools.cycle(pl.rcParams["axes.prop_cycle"].by_key()["color"])
    pl.figure(figsize=(12,12))

    for key in msdata_all:
        weight = msdata_all[key]['weight']
        uvd = msdata_all[key]['uvdist']
        flag = msdata_all[key]['flag']
        okflag = ~np.any(flag, axis=(0,1))
        pl.semilogy(uvd[okflag], weight.reshape(weight.shape[0], uvd.size).T[okflag.ravel(), :], ',', color=next(colors), label=key, alpha=0.5)
    pl.legend(loc='best')
    pl.set_title(vis)
    _=pl.xlabel("UV Distance (m)")
    _=pl.ylabel("Weight")

    colors = itertools.cycle(pl.rcParams["axes.prop_cycle"].by_key()["color"])
    pl.figure(figsize=(12,12))

    for key in msdata_all:
        amp = msdata_all[key]['amplitude']
        uvd = msdata_all[key]['uvdist']
        flag = msdata_all[key]['flag']
        okflag = ~np.any(flag, axis=(0,1))
        pl.plot(uvd[okflag], amp.reshape(amp.shape[0]*amp.shape[1], uvd.size).T[okflag.ravel(), :], ',', color=next(colors), label=key, alpha=0.5)
    pl.legend(loc='best')
    pl.set_title(vis)
    _=pl.xlabel("UV Distance (m)")
    _=pl.ylabel("Amplitude")

Reading a getsf catalog

To read a getsf catalog, you can do:

from astropy.io import ascii
cat = ascii.read('getsf.cat', data_start=0, format='commented_header', header_start=110, comment=\"!\")

though header_start may be different; I've seen 110 and 120 at least.

Determing ALMA Array Configurations

ALMA's named array configurations are not accessible through CASA. They are available from https://almascience.eso.org/observing/observing-configuration-schedule/prior-cycle-observing-and-configuration-schedule

This snippet will let you grab the tables and create a mapping from date to array config name.

However, it's a HUGE waste of time, because these data are stored directly in the MS!

This is the right way:

tb.open(vis+'/ASDM_EXECBLOCK')
tb.getcol('configName')
# array(['C43-3'], dtype='<U16')

This is the hacky, reconstructed from the website, bad way:

import requests
from bs4 import BeautifulSoup
from astropy.table import Table, vstack
from astropy.io import ascii
from astropy.time import Time

url = "https://almascience.eso.org/observing/observing-configuration-schedule/prior-cycle-observing-and-configuration-schedule"

response = requests.get(url)
response.raise_for_status()
soup = BeautifulSoup(response.text)
tables = soup.findAll('table', class_="grid listing")

def clean_lines(soup, badrows=['Long Baseline Campaign',
                               'February Maintenance Period',
                               'End of Cycle',
                               'Engineering/Software',
                               'Engineering/Software Time',
                                   ]):
    for tr in soup('tr'):
        if any(bad in str(tr) for bad in badrows):
             _=tr.extract()
    return soup

tables = soup.findAll('table', class_="grid listing")
tables = [ascii.read(str(clean_lines(tbl)), format='html') for tbl in tables]
tables[3] = tables[3][:-1]
tables[3]['Block'] = tables[3]['Block'].astype('int')
tables[4]['maximumrecoverablescale2(")'] = tables[4]['maximumrecoverablescale2(")'].astype('str')

stacked = vstack(tables)
start_times = Time(stacked['Start date'])
end_times = Time(stacked['End date'])

Then, say you have a list of measurement sets mses, you can look up the array configuration for each date and field. The choice of fields[0] here only makes sense if your data had a single target; it's a bad choice if there are multiple targets in a scheduling block:

import json
from casatools import msmetadata
msmd = msmetadata()

results = {}

for vis in mses:
    msmd.open(mses)
    obstime = Time(msmd.timerangeforobs(0)['begin']['m0']['value'], format='mjd')
    fieldnames = np.array(msmd.fieldnames())
    fields = np.unique(fieldnames[msmd.fieldsforintent('OBSERVE_TARGET#ON_SOURCE')])
    msmd.close()

    array_config = stacked[(obstime > start_times) & (obstime < end_times)]['Approx\xa0Config.']
    if fields[0] in results:
        results[fields[0]][obstime.strftime('%Y-%m-%d')] = array_config[0]
    else:
        results[fields[0]] = {obstime.strftime('%Y-%m-%d'): array_config[0]}

with open('array_configurations.json', 'w') as fh:
    json.dump(results, fh)

CASA: Logfiles & tclean from pysynthesisimager

2020/05/09

A few brief notes I don't want to forget, so I'm posting them here:

The name of the log file being used by the currently active CASA session is in casalog.logfile().

I'm trying to get the pysynthesisimager version of tclean to work. This is, in theory, "exactly the same thing as tclean", except you (the user) have the ability to control individual steps. The process is documented in the tclean documentation under examples.

However, it is not a one-to-one mapping. tclean is not implemented using the exact script quoted in the examples. The differences are the variety that shouldn't matter but inevitably do.

I'm working on a translation of tclean's input parameters to pysynthesisimager's, but it's not complete: https://gist.github.com/keflavich/a9953c88cbe098b5e9bb4bdac6485f09

That was dumb anyway, because https://github.com/radio-astro/casa/blob/master/gcwrap/python/scripts/task_tclean.py does everything I really needed.

My version saves the image, residual, and model at each major cycle.

CASA Cross-Calibration Application

2019/08/22

On several independent occasions, I've had data where the continuum is bright enough to self-calibrate, but the lines aren't. In these cases, I want to apply the continuum-derived selfcal solutions to the line data.

In all of these cases, I have several independent execution blocks that I've split and downsampled before combining to make the continuum image.

I got stuck trying to make this happen, but we reconvened on the topic in early 2021, and Manuel Fernandez and Roberto Galvan-Madrid and I came up with this tool to match spectral windows in a concatenated data set to the appropriate windows in the original: https://github.com/ALMA-IMF/reduction/blob/master/misc/applycal_tests.py

« Page 2 / 7 »

Recent Posts