API

API#

This page documents all of the available components of the MPoL package. If you do not see something that you think should be documented, please raise an issue.

Utilities#

mpol.utils.torch2npy(tensor)[source]#: Make a copy of a PyTorch tensor on the CPU in numpy format, e.g. for plotting

mpol.utils.ground_cube_to_packed_cube(ground_cube)[source]#

Converts a Ground Cube to a Packed Visibility Cube for visibility-plane work. See Units and Conventions for more details.

Parameters:

ground_cube – a previously initialized Ground Cube object (cube (3D torch tensor
`` (of shape) –

Returns:

3D image cube of shape (nchan, npix, npix); The resulting: array after applying torch.fft.fftshift to the input arg; i.e Returns a Packed Visibility Cube.

Return type:

torch.double

mpol.utils.packed_cube_to_ground_cube(packed_cube) → torch.Tensor[source]#

Converts a Packed Visibility Cube to a Ground Cube for visibility-plane work. See Units and Conventions for more details.

Parameters:

packed_cube – a previously initialized Packed Cube object (cube (3D torch tensor
`` (of shape) –

Returns:

3D image cube of shape (nchan, npix, npix); The resulting: array after applying torch.fft.fftshift to the input arg; i.e Returns a Ground Cube.

Return type:

torch.double

mpol.utils.sky_cube_to_packed_cube(sky_cube)[source]#

Converts a Sky Cube to a Packed Image Cube for image-plane work. See Units and Conventions for more details.

Parameters:

sky_cube – a previously initialized Sky Cube object with RA increasing to the left (cube (3D torch tensor of shape (nchan, npix, npix)))

Returns:

3D image cube of shape (nchan, npix, npix); The resulting: array after applying torch.fft.fftshift to the torch.flip() of the RA axis; i.e Returns a Packed Image Cube.

Return type:

torch.double

mpol.utils.packed_cube_to_sky_cube(packed_cube)[source]#

Converts a Packed Image Cube to a Sky Cube for image-plane work. See Units and Conventions for more details.

Parameters:

packed_cube – a previously initialized Packed Image Cube object (cube (3D torch
`` (tensor of shape) –

Returns:

3D image cube of shape (nchan, npix, npix); The resulting: array after applying torch.fft.fftshift to the torch.flip() of the RA axis; i.e Returns a Sky Cube.

Return type:

torch.double

mpol.utils.get_Jy_arcsec2(T_b, nu=230000000000.0)[source]#

Calculate specific intensity from the brightness temperature, using the Rayleigh-Jeans definition.

Parameters:

T_b – brightness temperature in [$K$]
nu – frequency (in Hz)

Returns:

specific intensity (in [$\mathrm{Jy}\, \mathrm{arcsec}^2]$)

Return type:

float

mpol.utils.log_stretch(x)[source]#

Apply a log stretch to the tensor.

Parameters:: tensor (PyTorch tensor) – input tensor $x$

Returns: $\ln(1 + |x|)$

mpol.utils.loglinspace(start, end, N_log, M_linear=3)[source]#

Return a logspaced array of bin edges, with the first M_linear cells being equal width. There is a one-cell overlap between the linear and logarithmic stretches of the array, since the last linear cell is also the first logarithmic cell, which means the total number of cells is M_linear + N_log - 1.

Parameters:

start (float) – starting cell left edge
end (float) – ending cell right edge
N_log (int) – number of logarithmically spaced bins
M_linear (int) – number of linearly (equally) spaced bins

mpol.utils.fftspace(width, N)[source]#

Delivers a (nearly) symmetric coordinate array that spans $N$ elements (where $N$ is even) from -width to +width, but ensures that the middle point lands on $0$. The array indices go from $0$ to $N -1.$

Parameters:

width (float) – the width of the array
N (int) – the number of elements in the array

Returns:

the fftspace array

Return type:

numpy.float64 1D array

mpol.utils.check_baselines(q, min_feasible_q=1.0, max_feasible_q=100000.0)[source]#

Check if baseline lengths are sensible for expected code unit of [klambda], or if instead they’re being supplied in [lambda].

Parameters:: q (array, unit = $k\lambda$) – Baseline distribution (all values must be non-negative).

min_feasible_qfloat, unit = $k\lambda$, default=1e0: Minimum baseline in code units expected for a dataset. The default value of 1e0 is a conservative value for ALMA, assuming a minimum antenna separation of ~12 m and maximum observing wavelength of 3.6 mm.
max_feasible_qfloat, unit = $k\lambda$, default=1e5: Maximum baseline in code units expected for a dataset. The default value of 1e5 is a conservative value for ALMA, assuming a maximum antenna separation of ~16 km and minimum observing wavelength of 0.3 mm.

mpol.utils.convert_baselines(baselines, freq=None, wle=None)[source]#

Convert baselines in meters to kilolambda. :param baselines: baselines in [m]. :type baselines: float or np.array :param freq: frequencies in [Hz]. :type freq: float or np.array :param wle: wavelengths in [m]. :type wle: float or np.array

Returns:: baselines in [klambda]
Return type:: (1D array nvis)

Notes

If baselines, freq or wle are numpy arrays, their shapes must be broadcast-able.

mpol.utils.broadcast_and_convert_baselines(u, v, chan_freq)[source]#

Convert baselines to kilolambda and broadcast to match shape of channel frequencies. :param u: baseline [m] :type u: 1D array nvis :param v: baseline [m] :type v: 1D array nvis :param chan_freq: frequencies [Hz] :type chan_freq: 1D array nchan

Returns:: (u, v) each of which are (nchan, nvis) arrays of baselines in [klambda]

mpol.utils.get_max_spatial_freq(cell_size, npix)[source]#

Calculate the maximum spatial frequency that the image can represent and still satisfy the Nyquist Sampling theorem.

Parameters:

cell_size (float) – the pixel size in arcseconds
npix (int) – the number of pixels in the image

Returns:

the maximum spatial frequency contained in the image (in kilolambda)

Return type:

max_freq

mpol.utils.get_maximum_cell_size(uu_vv_point)[source]#

Calculate the maximum possible cell_size that will still Nyquist sample the uu or vv point. Note: not q point.

Parameters:: uu_vv_point (float) – a single spatial frequency. Units of [$\mathrm{k}\lambda$].
Returns:: cell_size (in arcsec)

mpol.utils.get_optimal_image_properties(image_width, u, v)[source]#

For an image of desired width, determine the maximum pixel size that ensures Nyquist sampling of the provided spatial frequency points, and the corresponding number of pixels to obtain the desired image width.

Parameters:

image_width (float, unit = arcsec) – Desired width of the image (for a square image of size image_width $\times$ image_width).
u (float or array, unit = $k\lambda$) – u and v spatial frequency points.
v (float or array, unit = $k\lambda$) – u and v spatial frequency points.

Returns:

cell_size (float, unit = arcsec) – Image pixel size required to Nyquist sample.
npix (int) – Number of pixels of cell_size to equal (or slightly exceed) the image width (npix will be rounded up and enforced as even).

Notes

No assumption or correction is made concerning whether the spatial frequency points are projected or deprojected.

mpol.utils.sky_gaussian_radians(l, m, a, delta_l, delta_m, sigma_l, sigma_m, Omega)[source]#

Calculates a 2D Gaussian on the sky plane with inputs in radians. The Gaussian is centered at delta_l, delta_m, has widths of sigma_l, sigma_m, and is rotated Omega degrees East of North.

To evaluate the Gaussian, internally first we translate to center

\[\begin{split}l' = l - \delta_l\\ m' = m - \delta_m\end{split}\]

then rotate coordinates

\[\begin{split}l'' = l' \cos \phi - m' \sin \phi \\ m'' = l' \sin \phi + m' \cos \phi\end{split}\]

and then evaluate the Gaussian

\[f_\mathrm{g}(l,m) = a \exp \left ( - \frac{1}{2} \left [ \left (\frac{l''}{\sigma_l} \right)^2 + \left( \frac{m''}{\sigma_m} \right )^2 \right ] \right )\]

Parameters:

l – units of [radians]
m – units of [radians]
a – amplitude prefactor
delta_l – offset [radians]
delta_m – offset [radians]
sigma_l – width [radians]
sigma_M – width [radians]
Omega – position angle of ascending node [degrees] east of north.

Returns:

2D Gaussian evaluated at input args with peak amplitude $a$

mpol.utils.sky_gaussian_arcsec(x, y, a, delta_x, delta_y, sigma_x, sigma_y, Omega)[source]#

Calculates a Gaussian on the sky plane using inputs in arcsec. This is a convenience wrapper to sky_gaussian_radians() that automatically converts from arcsec to radians.

Parameters:

x – equivalent to l, but in units of [arcsec]
y – equivalent to m, but in units of [arcsec]
a – amplitude prefactor
delta_x – offset [arcsec]
delta_y – offset [arcsec]
sigma_x – width [arcsec]
sigma_y – width [arcsec]
Omega – position angle of ascending node [degrees] east of north.

Returns:

2D Gaussian evaluated at input args with peak amplitude $a$

mpol.utils.fourier_gaussian_lambda_radians(u, v, a, delta_l, delta_m, sigma_l, sigma_m, Omega)[source]#

Calculate the Fourier plane Gaussian $F_\mathrm{g}(u,v)$ corresponding to the Sky plane Gaussian $f_\mathrm{g}(l,m)$ in sky_gaussian_radians(), using analytical relationships. The Fourier Gaussian is parameterized using the sky plane centroid (delta_l, delta_m), widths (sigma_l, sigma_m) and rotation (Omega). Assumes that a was in units of $\mathrm{Jy}/\mathrm{steradian}$.

Parameters:

u – l in units of [lambda]
v – m in units of [lambda]
a – amplitude prefactor, units of $\mathrm{Jy}/\mathrm{steradian}$.
delta_x – offset [radians]
delta_y – offset [radians]
sigma_x – width [radians]
sigma_y – width [radians]
Omega – position angle of ascending node [degrees] east of north.

Returns:

2D Gaussian evaluated at input args

The following is a description of how we derived the analytical relationships. In what follows, all $l$ and $m$ coordinates are assumed to be in units of radians and all $u$ and $v$ coordinates are assumed to be in units of $\lambda$.

We start from Fourier dual relationships in Bracewell’s The Fourier Transform and Its Applications

\[f_0(l, m) \leftrightharpoons F_0(u, v)\]

where the sky-plane and Fourier-plane Gaussians are

\[f_0(l,m) = a \exp \left ( -\pi [l^2 + m^2] \right)\]

and

\[F_0(u,v) = a \exp \left ( -\pi [u^2 + v^2] \right),\]

respectively. The sky-plane Gaussian has a maximum value of $a$.

We will use the similarity, rotation, and shift theorems to turn $f_0$ into a form matching $f_\mathrm{g}$, which simultaneously turns $F_0$ into $F_\mathrm{g}(u,v)$.

The similarity theorem states that (in 1D)

\[f(bl) = \frac{1}{|b|}F\left(\frac{u}{b}\right).\]

First, we scale $f_0$ to include sigmas. Let

\[f_1(l, m) = a \exp \left(-\frac{1}{2} \left [\left(\frac{l}{\sigma_l}\right)^2 + \left( \frac{m}{\sigma_m} \right)^2 \right] \right).\]

i.e., something we might call a normalized Gaussian function. Phrased in terms of $f_0$, $f_1$ is

\[f_1(l, m) = f_0\left ( \frac{l}{\sigma_l \sqrt{2 \pi}},\, \frac{m}{\sigma_m \sqrt{2 \pi}}\right).\]

Therefore, according to the similarity theorem, the equivalent $F_1(u,v)$ is

\[F_1(u, v) = \sigma_l \sigma_m 2 \pi F_0 \left( \sigma_l \sqrt{2 \pi} u,\, \sigma_m \sqrt{2 \pi} v \right),\]

\[F_1(u, v) = a \sigma_l \sigma_m 2 \pi \exp \left ( -2 \pi^2 [\sigma_l^2 u^2 + \sigma_m^2 v^2] \right).\]

Next, we rotate the Gaussian to match the sky plane rotation. A rotation $\Omega$ in the sky plane is carried out in the same direction in the Fourier plane,

\[\begin{split}u' = u \cos \Omega - v \sin \Omega \\ v' = u \sin \Omega + v \cos \Omega\end{split}\]

such that

\[\begin{split}f_2(l, m) = f_1(l', m') \\ F_2(u, v) = F_1(u', m')\end{split}\]

Finally, we translate the sky plane Gaussian by amounts $\delta_l$, $\delta_m$, which corresponds to a phase shift in the Fourier plane Gaussian. The image plane translation is

\[f_3(l,m) = f_2(l - \delta_l, m - \delta_m)\]

According to the shift theorem, the equivalent $F_3(u,v)$ is

\[F_3(u,v) = \exp\left (- 2 i \pi [\delta_l u + \delta_m v] \right) F_2(u,v)\]

We have arrived at the corresponding Fourier Gaussian, $F_\mathrm{g}(u,v) = F_3(u,v)$. The simplified equation is

\[F_\mathrm{g}(u,v) = a \sigma_l \sigma_m 2 \pi \exp \left ( -2 \pi^2 \left [\sigma_l^2 u'^2 + \sigma_m^2 v'^2 \right] - 2 i \pi \left [\delta_l u + \delta_m v \right] \right).\]

N.B. that we have mixed primed ($u'$) and unprimed ($u$) coordinates in the same equation for brevity.

Finally, the same Fourier dual relationship holds

\[f_\mathrm{g}(l,m) \leftrightharpoons F_\mathrm{g}(u,v)\]

mpol.utils.fourier_gaussian_klambda_arcsec(u, v, a, delta_x, delta_y, sigma_x, sigma_y, Omega)[source]#

Calculate the Fourier plane Gaussian $F_\mathrm{g}(u,v)$ corresponding to the Sky plane Gaussian $f_\mathrm{g}(l,m)$ in sky_gaussian_arcsec(), using analytical relationships. The Fourier Gaussian is parameterized using the sky plane centroid (delta_l, delta_m), widths (sigma_l, sigma_m) and rotation (Omega). Assumes that a was in units of $\mathrm{Jy}/\mathrm{arcsec}^2$.

Parameters:

u – l in units of [klambda]
v – m in units of [klambda]
a – amplitude prefactor, units of $\mathrm{Jy}/\mathrm{arcsec}^2$.
delta_x – offset [arcsec]
delta_y – offset [arcsec]
sigma_x – width [arcsec]
sigma_y – width [arcsec]
Omega – position angle of ascending node [degrees] east of north.

Returns:

2D Fourier Gaussian evaluated at input args

Coordinates#

class mpol.coordinates.GridCoords(cell_size: float, npix: int)[source]#

The GridCoords object uses desired image dimensions (via the cell_size and npix arguments) to define a corresponding Fourier plane grid.

Parameters:

cell_size (float) – width of a single square pixel in [arcsec]
npix (int) – number of pixels in the width of the image

The Fourier grid is defined over the domain $[-u,+u]$, $[-v,+v]$, even though for real images, technically we could use an RFFT grid from $[0,+u]$ to $[-v,+v]$. The reason we opt for a full FFT grid in this instance is implementation simplicity.

Images (and their corresponding Fourier transform quantities) are represented as two-dimensional arrays packed as [y, x] and [v, u]. This means that an image with dimensions (npix, npix) will also have a corresponding FFT Fourier grid with shape (npix, npix). The native GridCoords representation assumes the Fourier grid (and thus image) are laid out as one might normally expect an image (i.e., no np.fft.fftshift has been applied).

After the object is initialized, instance variables can be accessed, for example

>>> myCoords = GridCoords(cell_size=0.005, 512)
>>> myCoords.img_ext

Variables:

dl – image-plane cell spacing in RA direction (assumed to be positive) [radians]
dm – image-plane cell spacing in DEC direction [radians]
img_ext – The length-4 list of (left, right, bottom, top) expected by routines like matplotlib.pyplot.imshow in the extent parameter assuming origin='lower'. Units of [arcsec]
packed_x_centers_2D – 2D array of l increasing, with fftshifted applied [arcseconds]. Useful for directly evaluating some function to create a packed cube.
packed_y_centers_2D – 2D array of m increasing, with fftshifted applied [arcseconds]. Useful for directly evaluating some function to create a packed cube.
sky_x_centers_2D – 2D array of l arranged for evaluating a sky image [arcseconds]. l coordinate increases to the left (as on sky).
sky_y_centers_2D – 2D array of m arranged for evaluating a sky image [arcseconds].
du – Fourier-plane cell spacing in East-West direction [$\mathrm{k}\lambda$]
dv – Fourier-plane cell spacing in North-South direction [$\mathrm{k}\lambda$]
u_centers – 1D array of cell centers in East-West direction [$\mathrm{k}\lambda$].
v_centers – 1D array of cell centers in North-West direction [$\mathrm{k}\lambda$].
u_edges – 1D array of cell edges in East-West direction [$\mathrm{k}\lambda$].
v_edges – 1D array of cell edges in North-South direction [$\mathrm{k}\lambda$].
u_bin_min – minimum u edge [$\mathrm{k}\lambda$]
u_bin_max – maximum u edge [$\mathrm{k}\lambda$]
v_bin_min – minimum v edge [$\mathrm{k}\lambda$]
v_bin_max – maximum v edge [$\mathrm{k}\lambda$]
max_grid – maximum spatial frequency enclosed by Fourier grid [$\mathrm{k}\lambda$]
vis_ext – length-4 list of (left, right, bottom, top) expected by routines like matplotlib.pyplot.imshow in the extent parameter assuming origin='lower'. Units of [$\mathrm{k}\lambda$]

check_data_fit(uu: Union[_SupportsArray[dtype], _NestedSequence[_SupportsArray[dtype]], bool, int, float, complex, str, bytes, _NestedSequence[Union[bool, int, float, complex, str, bytes]]], vv: Union[_SupportsArray[dtype], _NestedSequence[_SupportsArray[dtype]], bool, int, float, complex, str, bytes, _NestedSequence[Union[bool, int, float, complex, str, bytes]]]) → None[source]#

Test whether loose data visibilities fit within the Fourier grid defined by cell_size and npix.

Parameters:

uu (np.array) – array of u spatial frequency coordinates. Units of [$\mathrm{k}\lambda$]
vv (np.array) – array of v spatial frequency coordinates. Units of [$\mathrm{k}\lambda$]

Returns:

True if all visibilities fit within the Fourier grid defined by [u_bin_min, u_bin_max, v_bin_min, v_bin_max]. Otherwise an AssertionError is raised on the first violated boundary.

Geometry#

The geometry package provides routines for projecting and de-projecting sky images.

mpol.geometry.flat_to_observer(x, y, omega=None, incl=None, Omega=None)[source]#

Rotate the frame to convert a point in the flat (x,y,z) frame to observer frame (X,Y,Z).

It is assumed that the +Z axis points towards the observer. It is assumed that the: model is flat in the (x,y) frame (like a flat disk), and so the operations involving z are neglected. But the model lives in 3D Cartesian space.

In order,

rotate about the z axis by an amount omega -> x1, y1, z1
rotate about the x1 axis by an amount -incl -> x2, y2, z2
rotate about the z2 axis by an amount Omega -> x3, y3, z3 = X, Y, Z

Inspired by exoplanet/keplerian.py

Parameters:

x (torch tensor) – A tensor representing the x coordinate in the plane of the orbit.
y (torch tensor) – A tensor representing the y coordinate in the plane of the orbit.
omega (torch float tensor) – A tensor representing an argument of periastron [radians] Default 0.0.
incl (torch float tensor) – A tensor representing an inclination value [radians]. Default 0.0.
Omega (torch float tensor) – A tensor representing the position angle of the ascending node in [radians]. Default 0.0

Returns:

Two tensors representing (X, Y) in the observer frame.

mpol.geometry.observer_to_flat(X, Y, omega=None, incl=None, Omega=None)[source]#

Rotate the frame to convert a point in the observer frame (X,Y,Z) to the flat (x,y,z) frame.

It is assumed that the +Z axis points towards the observer. The rotation operations are the inverse of the flat_to_observer() operations.

In order,

inverse rotation about the Z axis by an amount Omega -> x2, y2, z2
inverse rotation about the x2 axis by an amount -incl -> x1, y1, z1
inverse rotation about the z1 axis by an amount omega -> x, y, z

Inspired by exoplanet/keplerian.py

Parameters:

X (torch tensor) – A tensor representing the x coodinate in the plane of the orbit.
Y (torch.tensor) – A tensor representing the y coodinate in the plane of the orbit.
omega (torch float tensor) – A tensor representing an argument of periastron [radians] Default 0.0.
incl (torch float tensor) – A tensor representing an inclination value [radians]. Default 0.0.
Omega (torch float tensor) – A tensor representing the position angle of the ascending node in [radians]. Default 0.0

Returns:

Two tensors representing (x, y) in the flat frame.

Gridding#

The gridding module provides two core classes in mpol.gridding.DataAverager and mpol.gridding.DirtyImager.

mpol.gridding.verify_no_hermitian_pairs(uu, vv, data, test_vis=5, test_channel=0)[source]#

Check that the dataset does not contain Hermitian pairs. Because the sky brightness $I_\nu$ is real, the visibility function $\mathcal{V}$ is Hermitian, meaning that

\[\mathcal{V}(u, v) = \mathcal{V}^*(-u, -v).\]

Most datasets (e.g., those extracted from CASA) will only record one visibility measurement per baseline and not include the duplicate Hermitian pair (to save storage space). This routine attempts to determine if the dataset contains Hermitian pairs or not by choosing one data point at a time and then searching the dataset to see if its Hermitian pair exists. The routine will declare that a dataset contains Hermitian pairs or not after it has searched test_vis number of data points. If 0 Hermitian pairs have been found for all test_vis points, then the dataset will be declared to have no Hermitian pairs. If test_vis Hermitian pairs were found for test_vis points searched, then the dataset will be declared to have Hermitian pairs. If more than 0 but fewer than test_vis Hermitian pairs were found for test_vis points, an error will be raised.

Gridding objects like mpol.gridding.DirtyImager will naturally augment an input dataset to include the Hermitian pairs, so that images of $I_\nu$ produced with the inverse Fourier transform turn out to be real.

Parameters:

uu (numpy array) – array of u spatial frequency coordinates. Units of [$\mathrm{k}\lambda$]
vv (numpy array) – array of v spatial frequency coordinates. Units of [$\mathrm{k}\lambda$]
data (numpy complex) – array of data values
test_vis (int) – the number of data points to search for Hermitian ‘matches’
test_channel (int) – the index of the channel to perform the check

Returns:

True if dataset does contain Hermitian pairs, False if not.

Return type:

boolean

class mpol.gridding.GridderBase(coords=None, uu=None, vv=None, weight=None, data_re=None, data_im=None)[source]#

This class is not designed to be used directly, but rather to be subclassed.

Subclasses will need to implement a _grid_visibilities(self,…) method.

The GridderBase object uses desired image dimensions (via the cell_size and npix arguments) to define a corresponding Fourier plane grid as a GridCoords object. A pre-computed GridCoords can be supplied in lieu of cell_size and npix, but all three arguments should never be supplied at once. For more details on the properties of the grid that is created, see the GridCoords documentation.

Subclasses will accept “loose” ungridded visibility data and store the arrays to the object as instance attributes. The input visibility data should be the set of visibilities over the full $[-u,u]$ and $[-v,v]$ domain, and should not contain Hermitian pairs (an error will be raised, if they are encountered). The visibilities can be 1d for a single continuum channel, or 2d for an image cube. If 1d, visibilities will be converted to 2d arrays of shape (1, nvis). Like the ImageCube class, after construction, GridderBase assumes that you are operating with a multi-channel set of visibilities. These routines will still work with single-channel ‘continuum’ visibilities, they will just have nchan = 1 in the first dimension of most products.

Parameters:

coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
uu (numpy array) – (nchan, nvis) array of u spatial frequency coordinates. Units of [$\mathrm{k}\lambda$]
vv (numpy array) – (nchan, nvis) array of v spatial frequency coordinates. Units of [$\mathrm{k}\lambda$]
weight (2d numpy array) – (nchan, nvis) array of thermal weights. Units of [$1/\mathrm{Jy}^2$]
data_re (2d numpy array) – (nchan, nvis) array of the real part of the visibility measurements. Units of [$\mathrm{Jy}$]
data_im (2d numpy array) – (nchan, nvis) array of the imaginary part of the visibility measurements. Units of [$\mathrm{Jy}$]

property ground_cube#

The visibility FFT cube fftshifted for plotting with imshow.

Returns:: the FFT of the image cube, in sky plane format.
Return type:: (torch.complex tensor, of shape (nchan, npix, npix))

class mpol.gridding.DataAverager(coords=None, uu=None, vv=None, weight=None, data_re=None, data_im=None)[source]#

The DataAverager object uses desired image dimensions (via the cell_size and npix arguments) to define a corresponding Fourier plane grid as a GridCoords object. A pre-computed GridCoords can be supplied in lieu of cell_size and npix, but all three arguments should never be supplied at once. For more details on the properties of the grid that is created, see the GridCoords documentation.

The DataAverager object accepts “loose” ungridded visibility data and stores the arrays to the object as instance attributes. The input visibility data should be the set of visibilities over the full $[-u,u]$ and $[-v,v]$ domain, and should not contain Hermitian pairs (an error will be raised, if they are encountered). (Note that, unlike DirtyImager, this class will not augment the dataset to include Hermitian pairs. This is by design, since Hermitian pairs should not be used in likelihood calculations).

The input visibilities can be 1d for a single continuum channel, or 2d for image cube. If 1d, visibilities will be converted to 2d arrays of shape (1, nvis). Like the ImageCube class, after construction, the DataAverager assumes that you are operating with a multi-channel set of visibilities. These routines will still work with single-channel ‘continuum’ visibilities, they will just have nchan = 1 in the first dimension of most products.

Once the DataAverager object is initialized with loose visibilities, you can average them and export them for use in Regularized Maximum Likelihood imaging with the mpol.gridding.DataAverager.to_pytorch_dataset() routine.

Example:

averager = gridding.DataAverager(
    coords=coords,
    uu=uu,
    vv=vv,
    weight=weight,
    data_re=data_re,
    data_im=data_im,
    )

dset = averager.to_pytorch_dataset()

Parameters:

coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
uu (numpy array) – (nchan, nvis) array of u spatial frequency coordinates. Units of [$\mathrm{k}\lambda$]
vv (numpy array) – (nchan, nvis) array of v spatial frequency coordinates. Units of [$\mathrm{k}\lambda$]
weight (2d numpy array) – (nchan, nvis) array of thermal weights. Units of [$1/\mathrm{Jy}^2$]
data_re (2d numpy array) – (nchan, nvis) array of the real part of the visibility measurements. Units of [$\mathrm{Jy}$]
data_im (2d numpy array) – (nchan, nvis) array of the imaginary part of the visibility measurements. Units of [$\mathrm{Jy}$]

to_pytorch_dataset(check_visibility_scatter=True, max_scatter=1.2)[source]#

Export gridded visibilities to a PyTorch dataset object.

Parameters:

check_visibility_scatter (bool) – whether the routine should check the standard deviation of visibilities in each within each $u,v$ cell ($\mathrm{cell}_{i,j}$) defined by self.coords. Default is True. A RuntimeError will be raised if any cell has a scatter larger than max_scatter.
max_scatter (float) – the maximum allowable standard deviation of visibility values in a given $u,v$ cell ($\mathrm{cell}_{i,j}$) defined by self.coords. Defaults to a factor of 120%.

Returns:

GriddedDataset with gridded visibilities.

class mpol.gridding.DirtyImager(coords=None, uu=None, vv=None, weight=None, data_re=None, data_im=None)[source]#

This class is mainly used for producing diagnostic “dirty” images of the visibility data.

The DirtyImager object uses desired image dimensions (via the cell_size and npix arguments) to define a corresponding Fourier plane grid as a GridCoords object. A pre-computed GridCoords can be supplied in lieu of cell_size and npix, but all three arguments should never be supplied at once. For more details on the properties of the grid that is created, see the GridCoords documentation.

The DirtyImager object accepts “loose” ungridded visibility data and stores the arrays to the object as instance attributes. The input visibility data should be the normal set of visibilities over the full $[-u,u]$ and $[-v,v]$ domain; internally the DirtyImager will automatically augment the dataset to include the complex conjugates, i.e. the ‘Hermitian pairs.’

The input visibilities can be 1d for a single continuum channel, or 2d for image cube. If 1d, visibilities will be converted to 2d arrays of shape (1, nvis). Like the ImageCube class, after construction, the DirtyImager assumes that you are operating with a multi-channel set of visibilities. These routines will still work with single-channel ‘continuum’ visibilities, they will just have nchan = 1 in the first dimension of most products.

Example:

imager = gridding.DirtyImager(
    coords=coords,
    uu=uu,
    vv=vv,
    weight=weight,
    data_re=data_re,
    data_im=data_im,
    )

img, beam = imager.get_dirty_image(weighting="briggs", robust=0.0)

Parameters:

cell_size (float) – width of a single square pixel in [arcsec]
npix (int) – number of pixels in the width of the image
coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
uu (numpy array) – (nchan, nvis) array of u spatial frequency coordinates. Units of [$\mathrm{k}\lambda$]
vv (numpy array) – (nchan, nvis) array of v spatial frequency coordinates. Units of [$\mathrm{k}\lambda$]
weight (2d numpy array) – (nchan, nvis) array of thermal weights. Units of [$1/\mathrm{Jy}^2$]
data_re (2d numpy array) – (nchan, nvis) array of the real part of the visibility measurements. Units of [$\mathrm{Jy}$]
data_im (2d numpy array) – (nchan, nvis) array of the imaginary part of the visibility measurements. Units of [$\mathrm{Jy}$]

get_dirty_beam_area(ntheta=24, single_channel_estimate=True)[source]#

Compute the effective area of the dirty beam for each channel. Assumes that the beam has already been generated by running get_dirty_image(). This is an approximate calculation involving a simple sum over all pixels out to the first null (zero crossing) of the dirty beam. This quantity is designed to approximate the conversion of image units from $[\mathrm{Jy}\,\mathrm{beam}^{-1}]$ to $[\mathrm{Jy}\,\mathrm{arcsec}^{-2}]$, even though units of $[\mathrm{Jy}\,\mathrm{dirty\;beam}^{-1}]$ are technically undefined.

Parameters:

ntheta (int) – number of azimuthal wedges to use for the 1st null calculation. More wedges will result in a more accurate estimate of dirty beam area, but will also take longer.
single_channel_estimate (bool) – If True (the default), use the area estimated from the first channel for all channels in the multi-channel image cube. If False, calculate the beam area for all channels.

Returns:

(1D numpy array float) beam area for each channel in units of $[\mathrm{arcsec}^{2}]$

get_dirty_image(weighting='uniform', robust=None, taper_function=None, unit='Jy/beam', check_visibility_scatter=True, max_scatter=1.2, **beam_kwargs)[source]#

Calculate the dirty image.

Parameters:

weighting (string) – The type of cell averaging to perform. Choices of "natural", "uniform", or "briggs", following CASA tclean. If "briggs", also specify a robust value.
robust (float) – If weighting='briggs', specify a robust value in the range [-2, 2]. robust=-2 approxmately corresponds to uniform weighting and robust=2 approximately corresponds to natural weighting.
taper_function (function reference) – a function assumed to be of the form $f(u,v)$ which calculates a prefactor in the range $[0,1]$ and premultiplies the visibility data. The function must assume that $u$ and $v$ will be supplied in units of $\mathrm{k}\lambda$. By default no taper is applied.
unit (string) – what unit should the image be in. Default is "Jy/beam". If "Jy/arcsec^2", then the effective area of the dirty beam will be used to convert from "Jy/beam" to "Jy/arcsec^2".
check_visibility_scatter (bool) – whether the routine should check the standard deviation of visibilities in each within each $u,v$ cell ($\mathrm{cell}_{i,j}$) defined by self.coords. Default is True. A RuntimeWarning will be raised if any cell has a scatter larger than max_scatter.
max_scatter (float) – the maximum allowable standard deviation of visibility values in a given $u,v$ cell ($\mathrm{cell}_{i,j}$) defined by self.coords. Defaults to a factor of 120%.
**beam_kwargs – all additional keyword arguments passed to get_dirty_beam_area() if unit="Jy/arcsec^2".

Returns:

2-tuple of (image, beam) where image is an (nchan, npix, npix) numpy array of the dirty image cube in units unit. beam is an numpy image cube with a dirty beam (PSF) for each channel. The units of the beam are always Jy/{dirty beam}, i.e., the peak of the beam is normalized to 1.0.

Datasets#

class mpol.datasets.GriddedDataset(*args: Any, **kwargs: Any)[source]#

Parameters:

coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
vis_gridded (torch complex) – the gridded visibility data stored in a “packed” format (pre-shifted for fft)
weight_gridded (torch double) – the weights corresponding to the gridded visibility data, also in a packed format
mask (torch boolean) – a boolean mask to index the non-zero locations of vis_gridded and weight_gridded in their packed format.
nchan (int) – the number of channels in the image (default = 1).
device (torch.device) – the desired device of the dataset. If None, defaults to current device.

After initialization, the GriddedDataset provides the non-zero cells of the gridded visibilities and weights as a 1D vector via the following instance variables. This means that any individual channel information has been collapsed.

Variables:

vis_indexed – 1D complex tensor of visibility data
weight_indexed – 1D tensor of weight values

If you index the output of the Fourier layer in the same manner using self.mask, then the model and data visibilities can be directly compared using a loss function.

classmethod from_image_properties(cell_size: float, npix: int, *, vis_gridded: torch.Tensor, weight_gridded: torch.Tensor, mask: torch.Tensor, nchan: int = 1, device: torch.device | str | None = None)[source]#

Alternative method to instantiate a GriddedDataset object from cell_size and npix.

Parameters:

cell_size (float) – the width of a pixel [arcseconds]
npix (int) – the number of pixels per image side
vis_gridded (torch complex) – the gridded visibility data stored in a “packed” format (pre-shifted for fft)
weight_gridded (torch double) – the weights corresponding to the gridded visibility data, also in a packed format
mask (torch boolean) – a boolean mask to index the non-zero locations of vis_gridded and weight_gridded in their packed format.
nchan (int) – the number of channels in the image (default = 1).

add_mask(mask: Union[_SupportsArray[dtype], _NestedSequence[_SupportsArray[dtype]], bool, int, float, complex, str, bytes, _NestedSequence[Union[bool, int, float, complex, str, bytes]]]) → None[source]#

Apply an additional mask to the data. Only works as a data limiting operation (i.e., mask is more restrictive than the mask already attached to the dataset).

Parameters:: mask (2D numpy or PyTorch tensor) – boolean mask (in packed format) to apply to dataset. Assumes input will be broadcast across all channels.

forward(modelVisibilityCube)[source]#

Parameters:

modelVisibilityCube (complex torch.tensor) – with shape (nchan, npix, npix) to be indexed. In “pre-packed” format, as in output from mpol.fourier.FourierCube.forward()

Returns:

1d torch tensor of indexed model samples collapsed: across cube dimensions.

Return type:

torch complex tensor

property ground_mask: torch.Tensor#

The boolean mask, arranged in ground format.

Returns:: 3D mask cube of shape (nchan, npix, npix)
Return type:: torch.boolean

class mpol.datasets.Dartboard(coords: GridCoords, q_edges: NDArray[floating[Any]] | None = None, phi_edges: NDArray[floating[Any]] | None = None)[source]#

A polar coordinate grid relative to a GridCoords object, reminiscent of a dartboard layout. The main utility of this object is to support splitting a dataset along radial and azimuthal bins for k-fold cross validation.

Parameters:

coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
q_edges (1D numpy array) – an array of radial bin edges to set the dartboard cells in $[\mathrm{k}\lambda]$. If None, defaults to 12 log-linearly radial bins stretching from 0 to the $q_\mathrm{max}$ represented by coords.
phi_edges (1D numpy array) – an array of azimuthal bin edges to set the dartboard cells in [radians], over the domain $[0, \pi]$, which is also implicitly mapped to the domain $[-\pi, \pi]$ to preserve the Hermitian nature of the visibilities. If None, defaults to 8 equal-spaced azimuthal bins stretched from $0$ to $\pi$.

classmethod from_image_properties(cell_size: float, npix: int, q_edges: NDArray[floating[Any]] | None = None, phi_edges: NDArray[floating[Any]] | None = None) → Dartboard[source]#

Alternative method to instantiate a Dartboard object from cell_size and npix.

Parameters:

cell_size (float) – the width of a pixel [arcseconds]
npix (int) – the number of pixels per image side
q_edges (1D numpy array) – an array of radial bin edges to set the dartboard cells in $[\mathrm{k}\lambda]$. If None, defaults to 12 log-linearly radial bins stretching from 0 to the $q_\mathrm{max}$ represented by coords.
phi_edges (1D numpy array) – an array of azimuthal bin edges to set the dartboard cells in [radians], over the domain $[0, \pi]$, which is also implicitly mapped to the domain $[-\pi, \pi]$ to preserve the Hermitian nature of the visibilities. If None, defaults to 8 equal-spaced azimuthal bins stretched from $0$ to $\pi$.

get_polar_histogram(qs: NDArray[floating[Any]], phis: NDArray[floating[Any]]) → NDArray[floating[Any]][source]#

Calculate a histogram in polar coordinates, using the bin edges defined by q_edges and phi_edges during initialization. Data coordinates should include the points for the Hermitian visibilities.

Parameters:

qs – 1d array of q values $[\mathrm{k}\lambda]$
phis – 1d array of datapoint azimuth values [radians] (must be the same length as qs)

Returns:

2d integer numpy array of cell counts, i.e., how many datapoints fell into each dartboard cell.

get_nonzero_cell_indices(qs: NDArray[floating[Any]], phis: NDArray[floating[Any]]) → NDArray[integer[Any]][source]#

Return a list of the cell indices that contain data points, using the bin edges defined by q_edges and phi_edges during initialization. Data coordinates should include the points for the Hermitian visibilities.

Parameters:

qs – 1d array of q values $[\mathrm{k}\lambda]$
phis – 1d array of datapoint azimuth values [radians] (must be the same length as qs)

Returns:

list of cell indices where cell contains at least one datapoint.

build_grid_mask_from_cells(cell_index_list: NDArray[integer[Any]]) → NDArray[np.bool_][source]#

Create a boolean mask of size (npix, npix) (in packed format) corresponding to the vis_gridded and weight_gridded quantities of the GriddedDataset .

Parameters:: cell_index_list (list) – list or iterable containing [q_cell, phi_cell] index pairs to include in the mask.

Returns: (numpy array) 2D boolean mask in packed format.

Images#

The images module provides the core functionality of MPoL via mpol.images.ImageCube.

class mpol.images.BaseCube(coords=None, nchan=1, pixel_mapping=None, base_cube=None)[source]#

A base cube of the same dimensions as the image cube. Designed to use a pixel mapping function $f_\mathrm{map}$ from the base cube values to the ImageCube domain.

\[I = f_\mathrm{map}(b)\]

The base_cube pixel values are set as PyTorch parameters.

Parameters:

cell_size (float) – the width of a pixel [arcseconds]
npix (int) – the number of pixels per image side
coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
nchan (int) – the number of channels in the base cube. Default = 1.
pixel_mapping (torch.nn) – a PyTorch function mapping the base pixel representation to the cube representation. If None, defaults to torch.nn.Softplus(). Output of the function should be in units of [$\mathrm{Jy}\,\mathrm{arcsec}^{-2}$].
base_cube (torch.double tensor, optional) – a pre-packed base cube to initialize the model with. If None, assumes torch.zeros. See Image Cube Packing for FFTs for more information on the expectations of the orientation of the input image.

forward()[source]#

Calculate the image representation from the base_cube using the pixel mapping

\[I = f_\mathrm{map}(b)\]

Returns : an image cube in units of [$\mathrm{Jy}\,\mathrm{arcsec}^{-2}$].

class mpol.images.HannConvCube(nchan, requires_grad=False)[source]#

This convolutional layer convolves an input cube by a small 3x3 filter with shape

\[\begin{split}\begin{bmatrix} 0.0625 & 0.1250 & 0.0625 \\ 0.1250 & 0.2500 & 0.1250 \\ 0.0625 & 0.1250 & 0.0625 \\ \end{bmatrix}\end{split}\]

which is the 2D version of the discretely-sampled response function corresponding to a Hann window, i.e., it is two 1D Hann windows multiplied together. This is a convolutional kernel in the image plane, and so effectively acts as apodization by a Hann window function in the Fourier domain. For more information, see the following Wikipedia articles on Window Functions in general and the Hann Window specifically.

The idea is that this layer would help naturally attenuate high spatial frequency artifacts by baking in a natural apodization in the Fourier plane.

Parameters:

nchan (int) – number of channels
requires_grad (bool) – keep kernel fixed

forward(cube)[source]#

Parameters:

cube (image) – a prepacked
cube –
example (for) –
ImageCube.forward() (from) –

Returns:

the FFT of the image cube, in packed format and of shape (nchan, npix, npix)

Return type:

torch.complex tensor

class mpol.images.ImageCube(coords=None, nchan=1, passthrough=False, cube=None)[source]#

The parameter set is the pixel values of the image cube itself. The pixels are assumed to represent samples of the specific intensity and are given in units of [$\mathrm{Jy}\,\mathrm{arcsec}^{-2}$].

All keyword arguments are required unless noted. The passthrough argument is essential for specifying whether the ImageCube object is the set of root parameters (passthrough==False) or if its simply a passthrough layer (pasthrough==True). In either case, ImageCube is essentially an identity layer, since no transformations are applied to the cube tensor. The main purpose of the ImageCube layer is to provide useful functionality around the cube tensor, such as returning a sky_cube representation and providing FITS writing functionility. In the case of passthrough==False, the ImageCube layer also acts as a container for the trainable parameters.

Parameters:

cell_size (float) – the width of a pixel [arcseconds]
npix (int) – the number of pixels per image side
coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
nchan (int) – the number of channels in the image
passthrough (bool) – if passthrough, assume ImageCube is just a layer as opposed to parameter base.
cube (torch.double tensor, of shape (nchan, npix, npix)) – (optional) a prepacked image cube to initialize the model with in units of [$\mathrm{Jy}\,\mathrm{arcsec}^{-2}$]. If None, assumes starting cube is torch.zeros. See Image Cube Packing for FFTs for more information on the expectations of the orientation of the input image.

forward(cube=None)[source]#

If the ImageCube object was initialized with passthrough=True, the cube argument is required. forward essentially just passes this on as an identity operation.

If the ImageCube object was initialized with passthrough=False, the cube argument is not permitted, and forward passes on the stored nn.Parameter cube as an identity operation.

Parameters:

cube (3D torch tensor of shape (nchan, npix, npix)) – only permitted if
passthrough=True. (the ImageCube object was initialized with) –

Returns: (3D torch.double tensor of shape (nchan, npix, npix)) as identity operation

property sky_cube#

The image cube arranged as it would appear on the sky.

Returns:: 3D image cube of shape (nchan, npix, npix)
Return type:: torch.double

property flux#

The spatially-integrated flux of the image. Returns a ‘spectrum’ with the flux in each channel in units of Jy.

Returns:: a 1D tensor of length (nchan)
Return type:: torch.double

to_FITS(fname='cube.fits', overwrite=False, header_kwargs=None)[source]#

Export the image cube to a FITS file.

Parameters:

fname (str) – the name of the FITS file to export to.
overwrite (bool) – if the file already exists, overwrite?
header_kwargs (dict) – Extra keyword arguments to write to the FITS header.

Returns:

None

Fourier#

The fourier module provides the core functionality of MPoL via mpol.fourier.FourierCube.

class mpol.fourier.FourierCube(coords: GridCoords, persistent_vis: bool = False)[source]#

This layer performs the FFT of an ImageCube and stores the corresponding dense FFT output as a cube. If you are using this layer in a forward-modeling RML workflow, because the FFT of the model is essentially stored as a grid, you will need to make the loss function calculation using a gridded loss function (e.g., mpol.losses.nll_gridded()) and a gridded dataset (e.g., mpol.datasets.GriddedDataset).

Parameters:

coords (GridCoords) – an object already instantiated from the GridCoords class.
persistent_vis (Boolean) – should the visibility cube be stored as part of the module s state_dict? If True, the state of the UV grid will be stored. It is recommended to use False for most applications, since the visibility cube will rarely be a direct parameter of the model.

classmethod from_image_properties(cell_size: float, npix: int, persistent_vis: bool = False) → FourierCube[source]#

Alternative method for instantiating a FourierCube from cell_size and npix

Parameters:

cell_size (float) – the width of an image-plane pixel [arcseconds]
npix (int) – the number of pixels per image side
persistent_vis (Boolean) – should the visibility cube be stored as part of the modules state_dict? If True, the state of the UV grid will be stored. It is recommended to use False for most applications, since the visibility cube will rarely be a direct parameter of the model.

Returns:

instantiated mpol.fourier.FourierCube object.

forward(cube: Tensor) → Tensor[source]#

Perform the FFT of the image cube on each channel.

Parameters:

cube (torch.double tensor) – a prepacked tensor of shape (nchan, npix, npix). For example, an image cube from ImageCube.forward()

Returns:

torch.complex tensor, of shape (nchan, npix, npix). The FFT of the: image cube, in packed format.

property ground_vis: Tensor#

The visibility cube in ground format cube fftshifted for plotting with imshow.

Returns:

the FFT of the: image cube, in sky plane format.

Return type:

(torch.complex tensor, of shape (nchan, npix, npix))

property ground_amp: Tensor#

The amplitude of the cube, arranged in unpacked format corresponding to the FFT of the sky_cube. Array dimensions for plotting given by self.coords.vis_ext.

Returns:: 3D amplitude cube of shape (nchan, npix, npix)
Return type:: torch.double

property ground_phase: Tensor#

The phase of the cube, arranged in unpacked format corresponding to the FFT of the sky_cube. Array dimensions for plotting given by self.coords.vis_ext.

Returns:: 3D phase cube of shape (nchan, npix, npix)
Return type:: torch.double

mpol.fourier.safe_baseline_constant_meters(uu: NDArray[floating[Any]], vv: NDArray[floating[Any]], freqs: NDArray[floating[Any]], coords: GridCoords, uv_cell_frac: float = 0.05) → bool[source]#

This routine determines whether the baselines can safely be assumed to be constant with channel when they converted from meters to units of kilolambda.

The antenna baselines are the same as a function of channel when they are measured in physical distance units, such as meters. However, when these baselines are converted to spatial frequency units, via

\[u = \frac{D}{\lambda},\]

it’s possible that the $u$ and $v$ values of each channel are significantly different if the $\lambda$ values of each channel are significantly different. This routine evaluates whether the maximum change in $u$ or $v$ across channels (when represented in kilolambda) is smaller than some threshold value, calculated as the fraction of a $u,v$ cell defined by coords.

If this function returns True, then it would be safe to proceed with parallelization in the mpol.fourier.NuFFT layer via the coil dimension.

Parameters:

uu (1D np.array) – a 1D array of length nvis array of the u (East-West) spatial frequency coordinate in units of [m]
vv (1D np.array) – a 1D array of length nvis array of the v (North-South) spatial frequency coordinate in units of [m]
freqs (1D np.array) – a 1D array of length nchan of the channel frequencies, in units of [Hz].
coords – a mpol.coordinates.GridCoords object which represents the image and uv-grid dimensions.
uv_cell_frac (float) – the maximum threshold for a change in $u$ or $v$ spatial frequency across channels, measured as a fraction of the $u,v$ cell defined by coords.

Returns:

True if it is safe to assume that the baselines are constant with: channel (at a tolerance of uv_cell_frac.) Otherwise returns False.

Return type:

boolean

mpol.fourier.safe_baseline_constant_kilolambda(uu: NDArray[floating[Any]], vv: NDArray[floating[Any]], coords: GridCoords, uv_cell_frac: float = 0.05) → bool[source]#

This routine determines whether the baselines can safely be assumed to be constant with channel, when the are represented in units of kilolambda.

Compared to mpol.fourier.safe_baseline_constant_meters, this function works with multidimensional arrays of uu and vv that are shape (nchan, nvis) and have units of kilolambda.

If this routine returns True, then it should be safe for the user to either average the baselines across channel or simply choose a single, representative channel. This would enable parallelization in the {class}`mpol.fourier.NuFFT` via the coil dimension.

Parameters:

uu (1D np.array) – a 1D array of length nvis array of the u (East-West) spatial frequency coordinate in units of [m]
vv (1D np.array) – a 1D array of length nvis array of the v (North-South) spatial frequency coordinate in units of [m]
freqs (1D np.array) – a 1D array of length nchan of the channel frequencies, in units of [Hz].
coords – a mpol.coordinates.GridCoords object which represents the image and uv-grid dimensions.
uv_cell_frac (float) – the maximum threshold for a change in $u$ or $v$ spatial frequency across channels, measured as a fraction of the $u,v$ cell defined by coords.

Returns:

True if it is safe to assume that the baselines are constant with: channel (at a tolerance of uv_cell_frac.) Otherwise returns False.

Return type:

boolean

class mpol.fourier.NuFFT(coords: GridCoords, nchan: int = 1)[source]#

This layer translates input from an mpol.images.ImageCube to loose, ungridded samples of the Fourier plane, corresponding to the $u,v$ locations provided. This layer is different than a mpol.Fourier.FourierCube in that, rather than producing the dense cube-like output from an FFT routine, it utilizes the non-uniform FFT or ‘NuFFT’ to interpolate directly to discrete $u,v$ locations. This is implemented using the KbNufft routines of the TorchKbNufft package.

Parameters:

coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
nchan (int) – the number of channels in the mpol.images.ImageCube. Default = 1.

classmethod from_image_properties(cell_size: float, npix: int, nchan: int = 1)[source]#

Instantiate a mpol.fourier.NuFFT object from image properties rather than a mpol.coordinates.GridCoords() instance.

Parameters:

cell_size (float) – the width of an image-plane pixel [arcseconds]
npix (int) – the number of pixels per image side
nchan (int) – the number of channels in the mpol.images.ImageCube. Default = 1.

Returns:

an instance of the mpol.fourier.NuFFT

forward(cube: Tensor, uu, vv, sparse_matrices: bool = False) → Tensor[source]#

Perform the FFT of the image cube for each channel and interpolate to the uu and vv points. This call should automatically take the best parallelization option as indicated by the shape of the uu and vv points. In general, you probably do not want to provide baselines that include Hermitian pairs.

Parameters:

cube (torch.double tensor) – of shape (nchan, npix, npix)). The cube should be a “prepacked” image cube, for example, from mpol.images.ImageCube.forward()
uu (array-like) – array of the u (East-West) spatial frequency coordinate [klambda].
vv (array-like) – array of the v (North-South) spatial frequency coordinate [klambda] (must be the same shape as uu)
sparse_matrices (bool) – If False, use the default table-based interpolation of TorchKbNufft.If True, use TorchKbNuFFT sparse matrices (generally slower but more accurate). Note that sparse matrices are incompatible with multi-channel uu and vv arrays (see below).

Returns:

Fourier samples of shape (nchan, nvis), evaluated: at the uu, vv points

Return type:

torch.complex tensor

Dimensionality: You should consider the dimensionality of your image and your visibility samples when using this method. If your image has multiple channels (nchan > 1), there is the possibility that the $u,v$ sample locations corresponding to each channel may be different. In ALMA/VLA applications, this can arise when continuum observations are taken over significant bandwidth, since the spatial frequency sampled by any pair of antennas is wavelength-dependent

\[u = \frac{D}{\lambda},\]

where $D$ is the projected baseline (measured in, say, meters) and $\lambda$ is the observing wavelength. In this application, the image-plane model could be the same for each channel, or it may vary with channel (necessary if the spectral slope of the source is significant).

On the other hand, with spectral line observations it will usually be the case that the total bandwidth of the observations is small enough such that the $u,v$ sample locations could be considered as the same for each channel. In spectral line applications, the image-plane model usually varies substantially with each channel.

This routine will determine whether the spatial frequencies are treated as constant based upon the dimensionality of the uu and vv input arguments.

If uu and vv have a shape of (nvis), then it will be assumed that
the spatial frequencies can be treated as constant with channel (and will invoke parallelization across the image cube nchan dimension using the ‘coil’ dimension of the TorchKbNufft package).
If the uu and vv have a shape of (nchan, nvis), then it will be
assumed that the spatial frequencies are different for each channel, and the spatial frequencies provided for each channel will be used (and will invoke parallelization across the image cube nchan dimension using the ‘batch’ dimension of the TorchKbNufft package).

Note that there is no straightforward, computationally efficient way to proceed if there are a different number of spatial frequencies for each channel. The best approach is likely to construct uu and vv arrays that have a shape of (nchan, nvis), such that all channels are padded with bogus $u,v$ points to have the same length nvis, and you create a boolean mask to keep track of which points are valid. Then, when this routine returns data points of shape (nchan, nvis), you can use that boolean mask to select only the valid $u,v$ points points.

Interpolation mode: You may choose the type of interpolation mode that KbNufft uses under the hood by changing the boolean value of sparse_matrices. If sparse_matrices=False, this routine will use the default table-based interpolation of TorchKbNufft. If sparse_matrices=True, the routine will calculate sparse matrices (which can be stored for later operations, as in {class}`~mpol.fourier.NuFFTCached`) and use them for the interpolation. This approach is likely to be more accurate but also slower. If Note that as of TorchKbNuFFT version 1.4.0, sparse matrices are not yet available when parallelizing using the ‘batch’ dimension — this will result in a warning. For most applications, we anticipate the accuracy of the table-based interpolation to be sufficiently accurate, but this could change depending on your problem.

class mpol.fourier.NuFFTCached(coords: GridCoords, uu, vv, nchan: int = 1, sparse_matrices: bool = True)[source]#

This layer translates input from an mpol.images.ImageCube directly to loose, ungridded samples of the Fourier plane, directly corresponding to the $u,v$ locations of the data. This layer is different than a mpol.Fourier.FourierCube in that, rather than producing the dense cube-like output from an FFT routine, it utilizes the non-uniform FFT or ‘NuFFT’ to interpolate directly to discrete $u,v$ locations that need not correspond to grid cell centers. This is implemented using the KbNufft routines of the TorchKbNufft package.

Dimensionality: One consideration when using this layer is the dimensionality of your image and your visibility samples. If your image has multiple channels (nchan > 1), there is the possibility that the $u,v$ sample locations corresponding to each channel may be different. In ALMA/VLA applications, this can arise when continuum observations are taken over significant bandwidth, since the spatial frequency sampled by any pair of antennas is wavelength-dependent

\[u = \frac{D}{\lambda},\]

This layer will determine whether the spatial frequencies are treated as constant based upon the dimensionality of the uu and vv input arguments.

If uu and vv have a shape of (nvis), then it will be assumed that the
spatial frequencies can be treated as constant with channel (and will invoke parallelization across the image cube nchan dimension using the ‘coil’ dimension of the TorchKbNufft package).
If the uu and vv have a shape of (nchan, nvis), then it will be
assumed that the spatial frequencies are different for each channel, and the spatial frequencies provided for each channel will be used (and will invoke parallelization across the image cube nchan dimension using the ‘batch’ dimension of the TorchKbNufft package).

Interpolation mode: You may choose the type of interpolation mode that KbNufft uses under the hood by changing the boolean value of sparse_matrices. For repeated evaluations of this layer (as might exist within an optimization loop), sparse_matrices=True is likely to be the more accurate and faster choice. If sparse_matrices=False, this routine will use the default table-based interpolation of TorchKbNufft. Note that as of TorchKbNuFFT version 1.4.0, sparse matrices are not yet available when parallelizing using the ‘batch’ dimension — this will result in a warning.

Parameters:

cell_size (float) – the width of an image-plane pixel [arcseconds]
npix (int) – the number of pixels per image side
coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
nchan (int) – the number of channels in the mpol.images.ImageCube. Default = 1.
uu (array-like) – a length nvis array (not including Hermitian pairs) of the u (East-West) spatial frequency coordinate [klambda]
vv (array-like) – a length nvis array (not including Hermitian pairs) of the v (North-South) spatial frequency coordinate [klambda]

classmethod from_image_properties(cell_size, npix, uu, vv, nchan=1, sparse_matrices=True)[source]#

Instantiate a mpol.fourier.NuFFT object from image properties rather than a mpol.coordinates.GridCoords() instance.

Parameters:

cell_size (float) – the width of an image-plane pixel [arcseconds]
npix (int) – the number of pixels per image side
nchan (int) – the number of channels in the mpol.images.ImageCube. Default = 1.

Returns:

an instance of the mpol.fourier.NuFFT

forward(cube)[source]#

Perform the FFT of the image cube for each channel and interpolate to the uu and vv points set at layer initialization. This call should automatically take the best parallelization option as set by the shape of the uu and vv points.

Parameters:

cube (torch.double tensor) – of shape (nchan, npix, npix)). The cube should be a “prepacked” image cube, for example, from mpol.images.ImageCube.forward()

Returns:

of shape (nchan, nvis), Fourier samples evaluated: corresponding to the uu, vv points set at initialization.

Return type:

torch.complex tensor

mpol.fourier.make_fake_data(image_cube: ImageCube, uu: NDArray[floating[Any]], vv: NDArray[floating[Any]], weight: NDArray[floating[Any]]) → tuple[NDArray[complexfloating[Any, Any]], ...][source]#

Create a fake dataset from a supplied mpol.images.ImageCube. See Making a Mock Dataset for more details on how to prepare a generic image for use in an ImageCube.

The provided visibilities can be 1d for a single continuum channel, or 2d for image cube. If 1d, visibilities will be converted to 2d arrays of shape (1, nvis).

Parameters:

imageCube (ImageCube) – the image layer to put into a fake dataset
uu (numpy array) – array of u spatial frequency coordinates, not including Hermitian pairs. Units of [$\mathrm{k}\lambda$]
vv (numpy array) – array of v spatial frequency coordinates, not including Hermitian pairs. Units of [$\mathrm{k}\lambda$]
weight (2d numpy array) – length array of thermal weights $w_i = 1/\sigma_i^2$. Units of [$1/\mathrm{Jy}^2$]

Returns:

a two tuple of the fake data. The first array is the mock dataset including noise, the second array is the mock dataset without added noise.

Return type:

(2-tuple)

mpol.fourier.get_vis_residuals(model, u_true, v_true, V_true, return_Vmod=False, channel=0)[source]#

Use mpol.fourier.NuFFT to get residuals between gridded model and loose (ungridded) data visiblities at data (u, v) coordinates

Parameters:

model (torch.nn.Module object) – Instance of the mpol.precomposed.SimpleNet class. Contains model visibilities.
u_true (array, unit=[klambda]) – Data u- and v-coordinates
v_true (array, unit=[klambda]) – Data u- and v-coordinates
V_true (array, unit=[Jy]) – Data visibility amplitudes
return_Vmod (bool, default=False) – Whether to return just the residual visibilities, or additionally the loose model visibilities
channel (int, default=0) – Channel (of model) to use to calculate residual visibilities

Returns:

vis_resid – Model loose residual visibility amplitudes of the form Re(V) + 1j * Im(V)

Return type:

array of complex

Precomposed Modules#

For convenience, we provide some “precomposed” modules which may be useful for simple imaging or modeling applications. In general, though, we encourage you to compose your own set of layers if your application requires it. The source code for a precomposed network can provide useful a starting point. We also recommend checking out the PyTorch documentation on modules.

class mpol.precomposed.SimpleNet(coords=None, nchan=1, base_cube=None)[source]#

A basic but fully functional network for RML imaging.

Parameters:

cell_size (float) – the width of a pixel [arcseconds]
npix (int) – the number of pixels per image side
coords (GridCoords) – an object already instantiated from the GridCoords class. If providing this, cannot provide cell_size or npix.
nchan (int) – the number of channels in the base cube. Default = 1.
base_cube – a pre-packed base cube to initialize the model with. If None, assumes torch.zeros.

After the object is initialized, instance variables can be accessed, for example

Variables:

bcube – the BaseCube instance
icube – the ImageCube instance
fcube – the FourierCube instance

For example, you’ll likely want to access the self.icube.sky_model at some point.

The idea is that SimpleNet can save you some keystrokes composing models by connecting the most commonly used layers together.

graph TD subgraph SimpleNet bc(BaseCube) --> HannConvCube HannConvCube --> ImageCube ImageCube --> FourierLayer end FourierLayer --> il([Loss]) ad[[Dataset]] --> il([Loss])

forward()[source]#

Feed forward to calculate the model visibilities. In this step, a BaseCube is fed to a HannConvCube is fed to a ImageCube is fed to a FourierCube to produce the visibility cube.

Returns: 1D complex torch tensor of model visibilities.

Losses#

The following loss functions are available to use in imaging. Many of the definitions follow those in Appendix A of EHT-IV 2019, including the regularization strength, which aspires to be similar across all terms, providing at least a starting point for tuning multiple loss functions.

If you don’t see a loss function you need, it’s easy to write your own directly within your optimization script. If you like it, please consider opening a pull request!

mpol.losses.chi_squared(model_vis, data_vis, weight)[source]#

Compute the $\chi^2$ between the complex data $\boldsymbol{V}$ and model: $M$ visibilities using

\[\chi^2(\boldsymbol{V}|\,\boldsymbol{\theta}) = \sum_i^N \frac{|V_i - M(u_i, v_i |\,\boldsymbol{\theta})|^2}{\sigma_i^2}\]

where $\sigma_i^2 = 1/w_i$. The sum is over all of the provided visibilities. This function is agnostic as to whether the sum should include the Hermitian conjugate visibilities, but be aware that the answer returned will be different between the two cases. We recommend not including the Hermitian conjugates.

Parameters:

model_vis (PyTorch complex) – array tuple of the model representing $\boldsymbol{V}$
data_vis (PyTorch complex) – array of the data values representing $M$
weight (PyTorch real) – array of weight values representing $w_i$

Returns:

the $\chi^2$ likelihood

Return type:

torch.double

mpol.losses.log_likelihood(model_vis, data_vis, weight)[source]#

Compute the log likelihood function $\ln\mathcal{L}$ between the complex data $\boldsymbol{V}$ and model $M$ visibilities using

\[\ln \mathcal{L}(\boldsymbol{V}|\,\boldsymbol{\theta}) = - \left ( N \ln 2 \pi + \sum_i^N \sigma_i^2 + \frac{1}{2} \chi^2(\boldsymbol{V}|\,\boldsymbol{\theta}) \right )\]

where $\chi^2$ is evaluated using mpol.losses.chi_squared().

This function is agnostic as to whether the sum should include the Hermitian conjugate visibilities, but be aware that the normalization of the answer returned will be different between the two cases. Inference of the parameter values should be unaffected. We recommend not including the Hermitian conjugates.

Parameters:

model_vis (PyTorch complex) – array tuple of the model representing $\boldsymbol{V}$
data_vis (PyTorch complex) – array of the data values representing $M$
weight (PyTorch real) – array of weight values representing $w_i$

Returns:

the $\ln\mathcal{L}$ log likelihood

Return type:

torch.double

mpol.losses.nll(model_vis, data_vis, weight)[source]#

Calculate a normalized “negative log likelihood” loss between the complex data $\boldsymbol{V}$ and model $M$ visibilities using

\[L_\mathrm{nll} = \frac{1}{2 N} \chi^2(\boldsymbol{V}|\,\boldsymbol{\theta})\]

where $\chi^2$ is evaluated using mpol.losses.chi_squared(). Visibilities may be any shape as long as all quantities have the same shape. Following EHT-IV 2019, we apply a prefactor $1/(2 N)$, where $N$ is the number of visibilities. The factor of 2 comes in because we must count real and imaginaries in the $\chi^2$ sum. This means that this normalized negative log likelihood loss function will have a minimum value of $L_mathrm{nll}(hat{boldsymbol{theta}}) approx 1$ for a well-fit model (regardless of the number of data points), making it easier to set the prefactor strengths of other regularizers relative to this value.

Note that this function should only be used in an optimization or point estimate situation. If it is used in any situation where uncertainties on parameter values are determined (such as Markov Chain Monte Carlo), it will return the wrong answer. This is because the relative scaling of $L_\mathrm{nll}$ with respect to parameter value is incorrect.

Parameters:

model_vis (PyTorch complex) – array tuple of the model representing $\boldsymbol{V}$
data_vis (PyTorch complex) – array of the data values representing $M$
weight (PyTorch real) – array of weight values representing $w_i$

Returns:

the normalized negative log likelihood likelihood loss

Return type:

torch.double

mpol.losses.chi_squared_gridded(modelVisibilityCube, griddedDataset)[source]#

Calculate the $\chi^2$ (corresponding to chi_squared()) using gridded data and model visibilities.

Parameters:

modelVisibilityCube (torch complex tensor) – torch tensor with shape (nchan, npix, npix) to be indexed by the mask from GriddedDataset. Assumes tensor is “pre-packed,” as in output from mpol.fourier.FourierCube.forward().
griddedDataset – instantiated GriddedDataset object

Returns:

the $\chi^2$ value

Return type:

torch.double

mpol.losses.log_likelihood_gridded(modelVisibilityCube, griddedDataset)[source]#

Calculate the log likelihood function $\ln\mathcal{L}$ (corresponding to log_likelihood()) using gridded data and model visibilities.

Parameters:

modelVisibilityCube (torch complex tensor) – torch tensor with shape (nchan, npix, npix) to be indexed by the mask from GriddedDataset. Assumes tensor is “pre-packed,” as in output from mpol.fourier.FourierCube.forward().
griddedDataset – instantiated GriddedDataset object

Returns:

the $\ln\mathcal{L}$ value

Return type:

torch.double

mpol.losses.nll_gridded(modelVisibilityCube, griddedDataset)[source]#

Calculate a normalized “negative log likelihood” (corresponding to nll()) using gridded data and model visibilities. Function will return the same value regardless of whether Hermitian pairs are included.

Parameters:

vis (torch complex tensor) – torch tensor with shape (nchan, npix, npix) to be indexed by the mask from GriddedDataset. Assumes tensor is “pre-packed,” as in output from mpol.fourier.FourierCube.forward().
griddedDataset – instantiated GriddedDataset object

Returns:

the normalized negative log likelihood likelihood loss

Return type:

torch.double

mpol.losses.entropy(cube, prior_intensity, tot_flux=10)[source]#

Calculate the entropy loss of a set of pixels following the definition in EHT-IV 2019.

Parameters:

cube (any tensor) – pixel values must be positive $I_i > 0$ for all $i$
prior_intensity (any tensor) – the prior value $p$ to calculate entropy against. Could be a single constant or an array the same shape as image.
tot_flux (float) – a fixed normalization factor; the user-defined target total flux density

Returns:

entropy loss

Return type:

torch.double

The entropy loss is calculated as

\[L = \frac{1}{\zeta} \sum_i I_i \; \ln \frac{I_i}{p_i}\]

mpol.losses.TV_image(sky_cube, epsilon=1e-10)[source]#

Calculate the total variation (TV) loss in the image dimension (R.A. and DEC). Following the definition in EHT-IV 2019 Promotes the image to be piecewise smooth and the gradient of the image to be sparse.

Parameters:

sky_cube (any 3D tensor) – the image cube array $I_{lmv}$, where $l$ is R.A. in $ndim=3$, $m$ is DEC in $ndim=2$, and $v$ is the channel (velocity or frequency) dimension in $ndim=1$. Should be in sky format representation.
epsilon (float) – a softening parameter in [$\mathrm{Jy}/\mathrm{arcsec}^2$]. Any pixel-to-pixel variations within each image slice greater than this parameter will have a significant penalty.

Returns:

total variation loss

Return type:

torch.double

\[L = \sum_{l,m,v} \sqrt{(I_{l + 1, m, v} - I_{l,m,v})^2 + (I_{l, m+1, v} - I_{l, m, v})^2 + \epsilon}\]

mpol.losses.TV_channel(cube, epsilon=1e-10)[source]#

Calculate the total variation (TV) loss in the channel dimension. Following the definition in EHT-IV 2019.

Parameters:

cube (any 3D tensor) – the image cube array $I_{lmv}$
epsilon (float) – a softening parameter in [$\mathrm{Jy}/\mathrm{arcsec}^2$]. Any channel-to-channel pixel variations greater than this parameter will have a significant penalty.

Returns:

total variation loss

Return type:

torch.double

\[L = \sum_{l,m,v} \sqrt{(I_{l, m, v + 1} - I_{l,m,v})^2 + \epsilon}\]

mpol.losses.edge_clamp(cube)[source]#

Promote all pixels at the edge of the image to be zero using an $L_2$ norm.

Parameters:: cube (any 3D tensor) – the array and pixel values
Returns:: edge loss
Return type:: torch.double

mpol.losses.sparsity(cube, mask=None)[source]#

Enforce a sparsity prior on the image cube using the $L_1$ norm. Optionally provide a boolean mask to apply the prior to only the True locations. For example, you might want this mask to be True for background regions.

Parameters:

cube (nchan, npix, npix) – tensor image cube
mask (boolean) – tensor array the same shape as cube. The sparsity prior will be applied to those pixels where the mask is True. Default is to apply prior to all pixels.

Returns:

sparsity loss calculated where mask == True

Return type:

torch.double

The sparsity loss calculated as

\[L = \sum_i | I_i |\]

mpol.losses.UV_sparsity(vis, qs, q_max)[source]#

Enforce a sparsity prior for all $q = \sqrt{u^2 + v^2}$ points larger than $q_\mathrm{max}$.

Parameters:

vis (torch.double) – visibility cube of (nchan, npix, npix//2 +1, 2)
qs – numpy array corresponding to visibility coordinates. Dimensionality of (npix, npix//2)
q_max (float) – maximum radial baseline

Returns:

UV sparsity loss above $q_\mathrm{max}$

Return type:

torch.double

mpol.losses.PSD(qs, psd, l)[source]#

Apply a loss function corresponding to the power spectral density using a Gaussian process kernel.

Assumes an image plane kernel of

\[k(r) = exp(-\frac{r^2}{2 \ell^2})\]

The corresponding power spectral density is

\[P(q) = (2 \pi \ell^2) exp(- 2 \pi^2 \ell^2 q^2)\]

Parameters:

qs (torch.double) – the radial UV coordinate (in kilolambda)
psd (torch.double) – the power spectral density cube
l (torch.double) – the correlation length in the image plane (in arcsec)

Returns:

the loss calculated using the power spectral density

Return type:

torch.double

mpol.losses.TSV(sky_cube)[source]#

Calculate the total square variation (TSV) loss in the image dimension (R.A. and DEC). Following the definition in EHT-IV 2019 Promotes the image to be edge smoothed which may be a better reoresentation of the truth image K. Kuramochi et al 2018.

Parameters:: sky_cube (any 3D tensor) – the image cube array $I_{lmv}$, where $l$ is R.A. in $ndim=3$, $m$ is DEC in $ndim=2$, and $v$ is the channel (velocity or frequency) dimension in $ndim=1$. Should be in sky format representation.
Returns:: total square variation loss
Return type:: torch.double

\[L = \sum_{l,m,v} (I_{l + 1, m, v} - I_{l,m,v})^2 + (I_{l, m+1, v} - I_{l, m, v})^2\]

Training and testing#

mpol.training.train_to_dirty_image(model, imager, robust=0.5, learn_rate=100, niter=1000)[source]#

Train against a dirty image of the observed visibilities using a loss function of the mean squared error between the RML model image pixel fluxes and the dirty image pixel fluxes. Useful for initializing a separate RML optimization loop at a reasonable starting image.

Parameters:

model (torch.nn.Module object) – A neural network module; instance of the mpol.precomposed.SimpleNet class.
imager (mpol.gridding.DirtyImager object) – Instance of the mpol.gridding.DirtyImager class.
robust (float, default=0.5) – Robust weighting parameter used to create a dirty image.
learn_rate (float, default=100) – Learning rate for optimization loop
niter (int, default=1000) – Number of iterations for optimization loop

Returns:

model – The input model updated with the state of the training to the dirty image

Return type:

torch.nn.Module object

class mpol.training.TrainTest(imager, optimizer, scheduler=None, regularizers={}, epochs=10000, convergence_tol=1e-05, train_diag_step=None, kfold=None, save_prefix=None, verbose=True)[source]#

Utilities for training and testing an MPoL neural network.

Parameters:

imager (mpol.gridding.DirtyImager object) – Instance of the mpol.gridding.DirtyImager class.
optimizer (torch.optim object) – PyTorch optimizer class for the training loop.
scheduler (torch.optim.lr_scheduler object, default=None) – Scheduler for adjusting learning rate during optimization.
regularizers (nested dict) –
Dictionary of image regularizers to use. For each, a dict of the strength (‘lambda’, float), whether to guess an initial value for lambda (‘guess’, bool), and other quantities needed to compute their loss term.

Example

{"sparsity":{"lambda":1e-3, "guess":False}, "entropy": {"lambda":1e-3, "guess":True, "prior_intensity":1e-10} }
epochs (int) – Number of training iterations, default=10000
convergence_tol (float) – Tolerance for training iteration stopping criterion as assessed by loss function (suggested <= 1e-3)
train_diag_step (int) – Interval at which training diagnostics are output. If None, no diagnostics will be generated.
kfold (int) – The k-fold of the current training set (for diagnostics)
save_prefix (str) – Prefix (path) used for saved figure names. If None, figures won’t be saved
verbose (bool) – Whether to print notification messages

loss_convergence(loss)[source]#

Estimate whether the loss function has converged by assessing its relative change over recent iterations.

Parameters:: loss (array) – Values of loss function over iterations (epochs). If len(loss) < 11, False will be returned, as convergence cannot be adequately assessed.
Return type:: True if the convergence criterion is met, else False.

loss_lambda_guess()[source]#

Set an initial guess for regularizer strengths $\lambda_{x}$ by comparing images generated with different visibility weighting.

The guesses update lambda values in self._regularizers.

loss_eval(vis, dataset, sky_cube=None)[source]#

Parameters:

vis (torch.complex tensor) – Model visibility cube (see mpol.fourier.FourierCube.forward)
dataset (dataset object) – Instance of the mpol.datasets.GriddedDataset class.
sky_cube (torch.double) – MPoL Ground Cube (see mpol.utils.packed_cube_to_ground_cube)

Returns:

loss – Value of loss function

Return type:

torch.double

train(model, dataset)[source]#

Trains a neural network, forward modeling a visibility dataset and evaluating the corresponding model image against the data, using PyTorch with gradient descent.

Parameters:

model (torch.nn.Module object) – A neural network module; instance of the mpol.precomposed.SimpleNet class.
dataset (PyTorch dataset object) – Instance of the mpol.datasets.GriddedDataset class.

Returns:

loss.item() (float) – Value of loss function at end of optimization loop
losses (list of float) – Loss value at each iteration (epoch) in the loop

test(model, dataset)[source]#

Test a model visibility cube against withheld data.

Parameters:

model (torch.nn.Module object) – A neural network module; instance of the mpol.precomposed.SimpleNet class.
dataset (PyTorch dataset object) – Instance of the mpol.datasets.GriddedDataset class.

Returns:

loss.item() – Value of loss function

Return type:

float

property regularizers#: Dict containing regularizers used and their strengths

property train_figure#: (fig, axes) of figure showing training diagnostics

Cross-validation#

class mpol.crossval.CrossValidate(coords, imager, learn_rate=0.3, regularizers={}, epochs=10000, convergence_tol=1e-05, schedule_factor=0.995, start_dirty_image=False, train_diag_step=None, kfolds=5, split_method='dartboard', dartboard_q_edges=None, dartboard_phi_edges=None, split_diag_fig=False, store_cv_diagnostics=False, save_prefix=None, verbose=True, device=None, seed=None)[source]#

Utilities to run a cross-validation loop (implicitly running a training optimization loop), in order to compare MPoL models with different hyperparameter values

Parameters:

coords (mpol.coordinates.GridCoords object) – Instance of the mpol.coordinates.GridCoords class.
imager (mpol.gridding.DirtyImager object) – Instance of the mpol.gridding.DirtyImager class.
learn_rate (float, default=0.3) – Initial learning rate
regularizers (nested dict, default={}) – Dictionary of image regularizers to use. For each, a dict of the strength (‘lambda’, float), whether to guess an initial value for lambda (‘guess’, bool), and other quantities needed to compute their loss term. Example: {“sparsity”:{“lambda”:1e-3, “guess”:False}, “entropy”: {“lambda”:1e-3, “guess”:True, “prior_intensity”:1e-10} }
epochs (int, default=10000) – Number of training iterations
convergence_tol (float, default=1e-5) – Tolerance for training iteration stopping criterion as assessed by loss function (suggested <= 1e-3)
schedule_factor (float, default=0.995) – For the torch.optim.lr_scheduler.ReduceLROnPlateau scheduler, factor to which the learning rate is reduced when learning rate stops decreasing
start_dirty_image (bool, default=False) – Whether to start the RML optimization loop by initializing the model image to a dirty image of the observed data. If False, the optimization loop will start with a blank image.
train_diag_step (int, default=None) – Interval at which training diagnostics are output. If None, no diagnostics will be generated.
kfolds (int, default=5) – Number of k-folds to use in cross-validation
split_method (str, default='dartboard') – Method to split full dataset into train/test subsets
dartboard_q_edges (list of float, default=None, unit=[klambda]) – Radial and azimuthal bin edges of the cells used to split the dataset if split_method`==`dartboard (see datasets.Dartboard)
dartboard_phi_edges (list of float, default=None, unit=[klambda]) – Radial and azimuthal bin edges of the cells used to split the dataset if split_method`==`dartboard (see datasets.Dartboard)
split_diag_fig (bool, default=False) – Whether to generate a diagnostic figure of dataset splitting into train/test sets.
store_cv_diagnostics (bool, default=False) – Whether to store diagnostics of the cross-validation loop.
save_prefix (str, default=None) – Prefix (path) used for saved figure names. If None, figures won’t be saved
verbose (bool, default=True) – Whether to print notification messages.
device (torch.device, default=None) – Which hardware device to perform operations on (e.g., ‘cuda:0’). ‘None’ defaults to current device.
seed (int, default=None) – Seed for random number generator used in splitting data

split_dataset(dataset)[source]#

Split a dataset into training and test subsets.

Parameters:: dataset (PyTorch dataset object) – Instance of the mpol.datasets.GriddedDataset class
Returns:: split_iterator – Iterator that provides a (train, test) pair of GriddedDataset for each k-fold
Return type:: iterator returning tuple

run_crossval(dataset)[source]#

Run a cross-validation loop for a model obtained with a given set of hyperparameters.

Parameters:: dataset (dataset object) – Instance of the mpol.datasets.GriddedDataset class
Returns:: cv_score – Dictionary with mean and standard deviation of cross-validation scores across all k-folds, and all raw scores
Return type:: dict

property model#: For the most recent kfold, trained model (SimpleNet class instance)

property regularizers#: For the most recent kfold, dict containing regularizers used and their strengths

property score#: Dict containing cross-val scores for all k-folds, and mean and standard deviation of these

property split_method#: String of the method used to split the dataset into train/test sets

property train_figure#: For the most recent kfold, (fig, axes) showing training progress

property split_figure#: (fig, axes) of train/test splitting diagnostic figure

property diagnostics#

models, regularizers, loss values

Type:: Dict containing diagnostics of the cross-validation loop across all kfolds

class mpol.crossval.RandomCellSplitGridded(dataset, k=5, seed=None, channel=0)[source]#

Split a GriddedDataset into $k$ subsets. Inherit the properties of the GriddedDataset. This object creates an iterator providing a (train, test) pair of GriddedDataset for each k-fold.

Parameters:

dataset (PyTorch dataset object) – Instance of the mpol.datasets.GriddedDataset class
k (int, default=5) – Number of k-folds (partitions) of dataset
seed (int, default=None) – Seed for PyTorch random number generator used to shuffle data before splitting
channel (int, default=0) – Channel of the dataset to use in determining the splits

Notes

Once initialized, iterate through the datasets like:

>>> split_iterator = crossval.RandomCellSplitGridded(dataset, k)
>>> for (train, test) in split_iterator: # iterate through `k` datasets
>>> ... # working with the n-th slice of `k` datasets
>>> ... # do operations with train dataset
>>> ... # do operations with test dataset

Treats dataset as a single-channel object with all data in channel.

The splitting doesn’t select (preserve) Hermitian pairs of visibilities.

All train splits have the highest 1% of cells by gridded weight

class mpol.crossval.DartboardSplitGridded(gridded_dataset: GriddedDataset, k: int, dartboard: Dartboard | None = None, seed: int | None = None, verbose: bool = True)[source]#

Split a GriddedDataset into $k$ non-overlapping chunks, internally partitioned by a Dartboard. Inherit the properties of the GriddedDataset. This object creates an iterator providing a (train, test) pair of GriddedDataset for each k-fold.

Parameters:

griddedDataset (GriddedDataset) – instance of the gridded dataset
k (int) – the number of subpartitions of the dataset
dartboard (Dartboard) – a pre-initialized Dartboard instance. If dartboard is provided, do not provide q_edges or phi_edges.
q_edges (1D numpy array) – an array of radial bin edges to set the dartboard cells in $[\mathrm{k}\lambda]$. If None, defaults to 12 log-linearly radial bins stretching from 0 to the $q_\mathrm{max}$ represented by coords.
phi_edges (1D numpy array) – an array of azimuthal bin edges to set the dartboard cells in [radians]. If None, defaults to 8 equal-spaced azimuthal bins stretched from $0$ to $\pi$.
seed (int) – (optional) numpy random seed to use for the permutation, for reproducibility

Once initialized, iterate through the datasets like

>>> cv = crossval.DartboardSplitGridded(dataset, k)
>>> for (train, test) in cv: # iterate among k datasets
>>> ... # working with the n-th slice of k datasets
>>> ... # do operations with train dataset
>>> ... # do operations with test dataset

Notes

All train splits have the cells belonging to the shortest dartboard baseline bin.

The number of points in the splits is in general not equal.

classmethod from_dartboard_properties(gridded_dataset: GriddedDataset, k: int, q_edges: NDArray[floating[Any]], phi_edges: NDArray[floating[Any]], seed: int | None = None, verbose: bool = True) → DartboardSplitGridded[source]#

Alternative method to initialize a DartboardSplitGridded object from Dartboard parameters.

Parameters:

griddedDataset (GriddedDataset) – instance of the gridded dataset
k (int) – the number of subpartitions of the dataset q_edges (1D numpy array): an array of radial bin edges to set the dartboard cells in $[\mathrm{k}\lambda]$. If None, defaults to 12 log-linearly radial bins stretching from 0 to the $q_\mathrm{max}$ represented by coords.
phi_edges (1D numpy array) – an array of azimuthal bin edges to set the dartboard cells in [radians]. If None, defaults to 8 equal-spaced azimuthal bins stretched from $0$ to $\pi$. seed (int): (optional) numpy random seed to use for the permutation, for reproducibility
verbose (bool) – whether to print notification messages

Analysis#

mpol.onedim.radialI(icube, geom, chan=0, bins=None)[source]#

Obtain a 1D (radial) brightness profile I(r) from an image cube.

Parameters:

icube (mpol.images.ImageCube object) – Instance of the MPoL images.ImageCube class
geom (dict) –

Dictionary of source geometry. Keys:

”incl”float, unit=[deg]
Inclination

”Omega”float, unit=[deg]
Position angle of the ascending node

”omega”float, unit=[deg]
Argument of periastron

”dRA”float, unit=[arcsec]
Phase center offset in right ascension. Positive is west of north.

”dDec”float, unit=[arcsec]
Phase center offset in declination.
chan (int, default=0) – Channel of the image cube corresponding to the desired image
bins (array, default=None, unit=[arcsec]) – Radial bin edges to use in calculating I(r). If None, bins will span the full image, with widths equal to the hypotenuse of the pixels

Returns:

bin_centers (array, unit=[arcsec]) – Radial coordinates of image at center of bins
Is (array, unit=[Jy / arcsec^2] (if image has these units)) – Azimuthally averaged pixel brightness at rs

mpol.onedim.radialV(fcube, geom, rescale_flux, chan=0, bins=None)[source]#

Obtain the 1D (radial) visibility model V(q) corresponding to a 2D MPoL image.

Parameters:

fcube (~mpol.fourier.FourierCube object) – Instance of the MPoL fourier.FourierCube class
geom (dict) –

Dictionary of source geometry. Keys:

”incl”float, unit=[deg]
Inclination

”Omega”float, unit=[deg]
Position angle of the ascending node

”omega”float, unit=[deg]
Argument of periastron

”dRA”float, unit=[arcsec]
Phase center offset in right ascension. Positive is west of north.

”dDec”float, unit=[arcsec]
Phase center offset in declination
rescale_flux (bool) – If True, the visibility amplitudes are rescaled to account for the difference between the inclined (observed) brightness and the assumed face-on brightness, assuming the emission is optically thick. The source’s integrated (2D) flux is assumed to be: $F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}$. No rescaling would be appropriate in the optically thin limit.
chan (int, default=0) – Channel of the image cube corresponding to the desired image
bins (array, default=None, unit=[klambda]) – Baseline bin edges to use in calculating V(q). If None, bins will span the model baseline distribution, with widths equal to the hypotenuse of the (u, v) coordinates

Returns:

bin_centers (array, unit=:math:[klambda]) – Baselines corresponding to u and v
Vs (array, unit=[Jy]) – Visibility amplitudes at q

Notes

This routine requires the frank package

Plotting#

mpol.plot.get_image_cmap_norm(image, stretch='power', gamma=1.0, asinh_a=0.02, symmetric=False)[source]#

Get a colormap normalization to apply to an image.

imagearray: 2D image array.
stretchstring, default = ‘power’: Transformation to apply to the colormap. ‘power’ is a power law stretch; ‘asinh’ is an arcsinh stretch.
gammafloat, default = 1.0: Index of power law normalization (see matplotlib.colors.PowerNorm). gamma=1.0 yields a linear colormap.
asinh_afloat, default = 0.02: Scale parameter for an asinh stretch.
symmetricbool, default=False: Whether the colormap is symmetric about 0

mpol.plot.get_residual_image(model, u, v, V, weights, robust=0.5)[source]#

Get a dirty image and colormap normalization for residual visibilities, the difference of observed visibilities and an MPoL model sampled at the observed (u,v) points.

Parameters:

model (torch.nn.Module object) – Instance of the mpol.precomposed.SimpleNet class. Contains model visibilities.
u (array, unit=[klambda]) – Data u- and v-coordinates
v (array, unit=[klambda]) – Data u- and v-coordinates
V (array, unit=[Jy]) – Data visibility amplitudes
weights (array, unit=[Jy^-2]) – Data weights
robust (float, default=0.5) – Robust weighting parameter used to create the dirty image of the residual visibilities

Returns:

im_resid (2D image array) – The residual image
norm_resid (Matplotlib colormap normalization) – Symmetric, linear colormap for im_resid

mpol.plot.plot_image(image, extent, cmap='inferno', norm=None, ax=None, clab='I [Jy arcsec$^{-2}$]', xlab='$\\Delta \\alpha \\cos \\delta$ [${}^{\\prime\\prime}$]', ylab='$\\Delta \\delta$ [${}^{\\prime\\prime}$]')[source]#

Wrapper for plt.imshow, with colorbar and colormap normalization.

Parameters:

image (array) – 2D image array.
extent (list, len=4) – x- and y-extents of image: [x-min, x-max, y-min, y-max] (see plt.imshow)
cmap (str, default="inferno) – Matplotlib colormap.
norm (Matplotlib colormap normalization, default=None) – Image colormap norm. If None, a linear normalization is generated with mpol.plot.get_image_cmap_norm
ax (Matplotlib axis instance, default=None) – Axis on which to plot the image. If None, a new figure is created.
clab (str, default=r"Jy arcsec$^{-2}$") – Colorbar axis label
xlab (str, default="RA offset [arcsec]") – Image x-axis label.
ylab (str, default="Dec offset [arcsec]") – Image y-axis label.

Returns:

im (Matplotlib imshow instance) – The plotted image.
cbar (Matplotlib colorbar instance) – Colorbar for the image.

mpol.plot.vis_histogram_fig(dataset, bin_quantity='count', bin_label=None, q_edges=None, phi_edges=None, q_edges1d=None, show_datapoints=False, save_prefix=None)[source]#

Generate a figure with 1d and 2d histograms of (u,v)-plane coverage. Histograms can show different data; see bin_quantity parameter.

Parameters:

dataset (mpol.datasets.GriddedDataset object) –
bin_quantity (str or numpy.ndarray, default='count') –
Which quantity to bin:
- ’count’ bins (u,v) points by their count
- ’weight’ bins points by the data weight (inherited from dataset)
- ’vis_real’ bins points by data Re(V)
- ’vis_imag’ bins points by data Im(V)
- A user-supplied numpy.ndarray to be used as ‘weights’ in np.histogram
bin_label (str, default=None) – Label for 1d histogram y-axis and 2d histogram colorbar.
q_edges (array, optional (default=None), unit=:math:[klambda]) – Radial bin edges for the 1d and 2d histogram. If None, defaults to 12 log-linearly radial bins over [0, 1.1 * maximum baseline in dataset].
phi_edges (array, optional (default=None), unit=[rad]) – Azimuthal bin edges for the 2d histogram. If None, defaults to 16 bins over [-pi, pi]
q_edges1d (array, optional (default=None), unit=:math:[klambda]) – Radial bin edges for a second 1d histogram. If None, defaults to 50 bins equispaced over [0, 1.1 * maximum baseline in dataset].
show_datapoints (bool, default = False) – Whether to overplot the raw visibilities in dataset on the 2d histogram.
save_prefix (string, default = None) – Prefix for saved figure name. If None, the figure won’t be saved

Returns:

fig (Matplotlib .Figure instance) – The generated figure
axes (Matplotlib ~.axes.Axes class) – Axes of the generated figure

Notes

No assumption or correction is made concerning whether the (u,v) distances are projected or deprojected.

mpol.plot.split_diagnostics_fig(splitter, channel=0, save_prefix=None)[source]#

Generate a figure showing (u,v) coverage in train and test sets split from a parent dataset.

Parameters:

splitter (mpol.crossval.RandomCellSplitGridded object) – Iterator that returns a (train, test) pair of GriddedDataset for each iteration.
channel (int, default=0) – Channel (of the datasets in splitter) to use to generate figure
save_prefix (string, default = None) – Prefix for saved figure name. If None, the figure won’t be saved

Returns:

fig (Matplotlib .Figure instance) – The generated figure
axes (Matplotlib ~.axes.Axes class) – Axes of the generated figure

Notes

No assumption or correction is made concerning whether the (u,v) distances are projected or deprojected.

mpol.plot.train_diagnostics_fig(model, losses=None, learn_rates=None, fluxes=None, old_model_image=None, old_model_epoch=None, kfold=None, epoch=None, channel=0, save_prefix=None)[source]#

Figure for model diagnostics at a given model state during an optimization loop.

Plots:

model image
flux of model image
gradient image
difference image between old_model_image and current model image
loss function
learning rate

Parameters:

model (torch.nn.Module object) – A neural network module; instance of the mpol.precomposed.SimpleNet class.
losses (list) – Loss value at each epoch in the training loop
learn_rates (list) – Learning rate at each epoch in the training loop
fluxes (list) – Total flux in model image at each epoch in the training loop
old_model_image (2D image array, default=None) – Model image of a previous epoch for comparison to current image
old_model_epoch (int) – Epoch of old_model_image
kfold (int, default=None) – Current cross-validation k-fold
epoch (int, default=None) – Current training epoch
channel (int, default=0) – Channel (of the datasets in splitter) to use to generate figure
save_prefix (str, default = None) – Prefix for saved figure name. If None, the figure won’t be saved

Returns:

fig (Matplotlib .Figure instance) – The generated figure
axes (Matplotlib ~.axes.Axes class) – Axes of the generated figure

mpol.plot.crossval_diagnostics_fig(cv, title='', save_prefix=None)[source]#

Figure for model diagnostics of a cross-validation run.

Plots:

loss evolution for each k-fold
cross-validation score per k-fold

Parameters:

cv (mpol.crossval.CrossValidate object) – Instance of the CrossValidate class produced by a cross-validation loop
title (str, default="") – Figure super-title
save_prefix (string, default = None) – Prefix for saved figure name. If None, the figure won’t be saved

Returns:

fig (Matplotlib .Figure instance) – The generated figure
axes (Matplotlib ~.axes.Axes class) – Axes of the generated figure

mpol.plot.image_comparison_fig(model, u, v, V, weights, robust=0.5, clean_fits=None, share_cscale=False, xzoom=[None, None], yzoom=[None, None], title='', channel=0, save_prefix=None)[source]#

Figure for comparison of MPoL model image to other image models.

Plots:

dirty image
MPoL model image
MPoL residual visibilities imaged
clean image (if a .fits file is supplied)

Parameters:

model (torch.nn.Module object) – A neural network; instance of the mpol.precomposed.SimpleNet class.
u (array, unit=[klambda]) – Data u- and v-coordinates
v (array, unit=[klambda]) – Data u- and v-coordinates
V (array, unit=[Jy]) – Data visibility amplitudes
weights (array, unit=[Jy^-2]) – Data weights
robust (float, default=0.5) – Robust weighting parameter used to create the dirty image of the observed visibilities and separately of the MPoL residual visibilities
clean_fits (str, default=None) – Path to a clean .fits image
share_cscale (bool, default=False) – Whether the MPoL model image, dirty image and clean image share the same colorscale
xzoom (list of float, default = [None, None]) – X- and y- axis limits to zoom the images to. xzoom and yzoom should both list values in ascending order (e.g. [-2, 3], not [3, -2])
yzoom (list of float, default = [None, None]) – X- and y- axis limits to zoom the images to. xzoom and yzoom should both list values in ascending order (e.g. [-2, 3], not [3, -2])
title (str, default="") – Figure super-title
channel (int, default=0) – Channel of the model to use to generate figure
save_prefix (string, default = None) – Prefix for saved figure name. If None, the figure won’t be saved

Returns:

fig (Matplotlib .Figure instance) – The generated figure
axes (Matplotlib ~.axes.Axes class) – Axes of the generated figure

mpol.plot.vis_1d_fig(model, u, v, V, weights, geom=None, rescale_flux=False, bin_width=20000.0, q_logx=True, title='', channel=0, save_prefix=None)[source]#

Figure for comparison of 1D projected MPoL model visibilities and observed

visibilities.

Plots:

Re(V): observed and MPoL model (projected unless geom is supplied)
Residual Re(V): observed - MPoL model (projected unless geom is supplied)
Im(V): observed and MPoL model (projected unless geom is supplied)
Residual Im(V): observed - MPoL model (projected unless geom is supplied)

Parameters:

model (torch.nn.Module object) – A neural network; instance of the mpol.precomposed.SimpleNet class.
u (array, unit=[klambda]) – Data u- and v-coordinates
v (array, unit=[klambda]) – Data u- and v-coordinates
V (array, unit=[Jy]) – Data visibility amplitudes
weights (array, unit=[Jy^-2]) – Data weights
geom (dict) –

Dictionary of source geometry. If passed in, visibilities will be

deprojected prior to plotting. Keys:

”incl”float, unit=[deg]
Inclination

”Omega”float, unit=[deg]
Position angle of the ascending node

”omega”float, unit=[deg]
Argument of periastron

”dRA”float, unit=[arcsec]
Phase center offset in right ascension. Positive is west of north.

”dDec”float, unit=[arcsec]
Phase center offset in declination.
rescale_flux (bool) –

If True, the visibility amplitudes are rescaled to account
for the difference between the inclined (observed) brightness and the assumed face-on brightness, assuming the emission is optically thick. The source’s integrated (2D) flux is assumed to be: $F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}$. No rescaling would be appropriate in the optically thin limit.
bin_width (float, default=20e3) – Bin size [klambda] for baselines
q_logx (bool, default=True) – Whether to plot visibilities in log-baseline
title (str, default="") – Figure super-title
channel (int, default=0) – Channel of the model to use to generate figure
save_prefix (string, default = None) – Prefix for saved figure name. If None, the figure won’t be saved

Returns:

fig (Matplotlib .Figure instance) – The generated figure
axes (Matplotlib ~.axes.Axes class) – Axes of the generated figure

Notes

This routine requires the frank package

mpol.plot.radial_fig(model, geom, u=None, v=None, V=None, weights=None, dist=None, rescale_flux=False, bin_width=20000.0, q_logx=True, title='', channel=0, save_prefix=None)[source]#

Figure for analysis of 1D (radial) brightness profile of MPoL model image, using a user-supplied geometry.

Plots:

MPoL model image
1D (radial) brightness profile extracted from MPoL image (supply dist to show second x-axis in [AU])
Deprojectd Re(V): binned MPoL model and observations (if u, v, V, weights supplied)
Deprojected Im(V): binned MPoL model and observations (if u, v, V, weights supplied)

Parameters:

model (torch.nn.Module object) – A neural network; instance of the mpol.precomposed.SimpleNet class.
geom (dict) –

Dictionary of source geometry. Used to deproject image and visibilities.

Keys:

”incl”float, unit=[deg]
Inclination

”Omega”float, unit=[deg]
Position angle of the ascending node

”omega”float, unit=[deg]
Argument of periastron

”dRA”float, unit=[arcsec]
Phase center offset in right ascension. Positive is west of north.

”dDec”float, unit=[arcsec]
Phase center offset in declination.
u (array, optional, unit=[klambda], default=None) – Data u- and v-coordinates
v (array, optional, unit=[klambda], default=None) – Data u- and v-coordinates
V (array, optional, unit=[Jy], default=None) – Data visibility amplitudes
weights (array, optional, unit=[Jy^-2], default=None) – Data weights
dist (float, optional, unit = [AU], default = None) – Distance to source, used to show second x-axis for I(r) in [AU]
rescale_flux (bool) – If True, the visibility amplitudes are rescaled to account for the difference between the inclined (observed) brightness and the assumed face-on brightness, assuming the emission is optically thick. The source’s integrated (2D) flux is assumed to be: $F = \cos(i) \int_r^{r=R}{I(r) 2 \pi r dr}$. No rescaling would be appropriate in the optically thin limit.
bin_width (float, default=20e3) – Bin size [klambda] in which to bin observed visibility points
q_logx (bool, default=True) – Whether to plot visibilities in log-baseline
title (str, default="") – Figure super-title
channel (int, default=0) – Channel of the model to use to generate figure
save_prefix (string, default = None) – Prefix for saved figure name. If None, the figure won’t be saved

Returns:

fig (Matplotlib .Figure instance) – The generated figure
axes (Matplotlib ~.axes.Axes class) – Axes of the generated figure

Notes

This routine requires the frank package

API

Contents

API#

Utilities#

Coordinates#

Geometry#

Gridding#

Datasets#

Images#

Fourier#

Precomposed Modules#

Losses#

Training and testing#

Cross-validation#

Analysis#

Plotting#