Utilities

Utilities#

mpol.utils.torch2npy(t: torch.Tensor) → ndarray[Any, dtype[_ScalarType_co]][source]#

Copy a tensor (potentially on the GPU) to the CPU and convert to a numpy np.ndarray, e.g., for visualization or further analysis with non-PyTorch scientific libraries.

Parameters:

t (torch.Tensor)

Return type:

np.ndarray

mpol.utils.ground_cube_to_packed_cube(ground_cube: torch.Tensor) → torch.Tensor[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) → 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: torch.Tensor) → torch.Tensor[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: torch.Tensor) → torch.Tensor[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: float, nu: float = 230000000000.0) → float[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.loglinspace(start: float, end: float, N_log: int, M_linear: int = 3) → ndarray[Any, dtype[floating[Any]]][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

Returns:

logspaced bin edges

Return type:

np.ndarray

mpol.utils.fftspace(width: float, N: int) → ndarray[Any, dtype[floating[Any]]][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=1000.0, max_feasible_q=100000000.0)[source]#

Check if baseline lengths are sensible for expected code unit of [\(\lambda\)], or if instead they’re being supplied in [\(\mathrm{k}\lambda\)].

Parameters:

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

min_feasible_qfloat, unit = \(\lambda\), default=1e3

Minimum baseline in code units expected for a dataset. The default value of 1e3 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 = \(\lambda\), default=1e8

Maximum baseline in code units expected for a dataset. The default value of 1e8 is a conservative value for ALMA, assuming a maximum antenna separation of ~16 km and minimum observing wavelength of 0.3 mm.

mpol.utils.get_max_spatial_freq(cell_size: float, npix: int) → float[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 [\(\lambda\)]

Return type:

float

mpol.utils.get_maximum_cell_size(uu_vv_point: float) → float[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 [\(\lambda\)].

Returns:

cell_size (in arcsec)

mpol.utils.get_optimal_image_properties(image_width: float, u: torch.Tensor, v: torch.Tensor) → tuple[float, int][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 (torch.Tensor , unit = \(\lambda\)) – u and v baselines.

  • v (torch.Tensor , unit = \(\lambda\)) – u and v baselines.

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

Assumes baselines are as-observed.

mpol.utils.sky_gaussian_radians(l: ndarray[Any, dtype[floating[Any]]], m: ndarray[Any, dtype[floating[Any]]], a: float, delta_l: float, delta_m: float, sigma_l: float, sigma_m: float, Omega: float) → ndarray[Any, dtype[floating[Any]]][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: ndarray[Any, dtype[floating[Any]]], y: ndarray[Any, dtype[floating[Any]]], a: float, delta_x: float, delta_y: float, sigma_x: float, sigma_y: float, Omega: float) → ndarray[Any, dtype[floating[Any]]][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: ndarray[Any, dtype[floating[Any]]], v: ndarray[Any, dtype[floating[Any]]], a: float, delta_l: float, delta_m: float, sigma_l: float, sigma_m: float, Omega: float) → ndarray[Any, dtype[floating[Any]]][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),\]

or

\[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_lambda_arcsec(u: ndarray[Any, dtype[floating[Any]]], v: ndarray[Any, dtype[floating[Any]]], a: float, delta_x: float, delta_y: float, sigma_x: float, sigma_y: float, Omega: float) → ndarray[Any, dtype[floating[Any]]][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 [lambda]

  • v – m in units of [lambda]

  • 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