1D analysis

.github/workflows/tests.yml

@@ -17,6 +17,7 @@ jobs:
         # in __init__. below we'll reinstall for the tests.
         run: |
           pip install astropy
+          pip install frank
           pip install .
       - name: Cache/Restore the .mpol folder cache
         uses: actions/cache@v3

mypy.ini

@@ -14,5 +14,8 @@ ignore_missing_imports = True
 [mypy-torchkbnufft.*]
 ignore_missing_imports = True
 
+[mypy-frank.*]
+ignore_missing_imports = True
+
 [mypy-fast_histogram.*]
 ignore_missing_imports = True

setup.py

@@ -33,6 +33,7 @@ def get_version(rel_path):
         "astropy",
         "tensorboard",
         "mypy",
+        "frank>=1.2.1",
     ],
     "docs": [
         "sphinx>=2.3.0",
@@ -55,10 +56,13 @@ def get_version(rel_path):
         "pyro-ppl",
         "arviz[all]"
     ],
+    "analysis": [
+        "frank>=1.2.1",
+    ],
 }
 
 EXTRA_REQUIRES["dev"] = (
-    EXTRA_REQUIRES["test"] + EXTRA_REQUIRES["docs"] + ["pylint", "black", "pre-commit"]
+    EXTRA_REQUIRES["test"] + EXTRA_REQUIRES["docs"] + EXTRA_REQUIRES["analysis"] + ["pylint", "black", "pre-commit"]
 )
 
 

src/mpol/geometry.py

@@ -28,7 +28,6 @@ def flat_to_observer(x, y, omega=None, incl=None, Omega=None):
     Returns:
         Two tensors representing ``(X, Y)`` in the observer frame.
     """
-
     # Rotation matrices result in a *clockwise* rotation of the axes, as defined using the righthand rule.
     # For example, looking down the z-axis, a positive angle will rotate the x,y axes clockwise.
     # A vector in the coordinate system will appear as though it has been rotated counter-clockwise.
@@ -70,7 +69,7 @@ def flat_to_observer(x, y, omega=None, incl=None, Omega=None):
 
 
 def observer_to_flat(X, Y, omega=None, incl=None, Omega=None):
-    """Rotate the from to convert a point in the observer frame (X,Y,Z) to the flat (x,y,z) frame.
+    """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 :func:`~mpol.geometry.flat_to_observer` operations.
 
@@ -90,9 +89,8 @@ def observer_to_flat(X, Y, omega=None, incl=None, Omega=None):
         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.
+        Two tensors representing ``(x, y)`` in the flat frame.
     """
-
     # Rotation matrices result in a *clockwise* rotation of the axes, as defined using the righthand rule.
     # For example, looking down the z-axis, a positive angle will rotate the x,y axes clockwise.
     # A vector in the coordinate system will appear as though it has been rotated counter-clockwise.

src/mpol/images.py

@@ -315,4 +315,4 @@ def to_FITS(self, fname="cube.fits", overwrite=False, header_kwargs=None):
         hdul = fits.HDUList([hdu])
         hdul.writeto(fname, overwrite=overwrite)
 
-        hdul.close()
+        hdul.close()

src/mpol/onedim.py

@@ -0,0 +1,170 @@
+import numpy as np 
+from mpol.utils import torch2npy
+
+def radialI(icube, geom, chan=0, bins=None):
+    r"""
+    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`
+    """
+
+    # projected Cartesian pixel coordinates [arcsec]
+    xx, yy = icube.coords.sky_x_centers_2D, icube.coords.sky_y_centers_2D
+
+    # shift image center to source center
+    xc, yc = xx - geom["dRA"], yy - geom["dDec"]
+
+    # deproject image
+    cos_PA = np.cos(geom["Omega"] * np.pi / 180)
+    sin_PA = np.sin(geom["Omega"] * np.pi / 180)
+    xd = xc * cos_PA - yc * sin_PA
+    yd = xc * sin_PA + yc * cos_PA
+    xd /= np.cos(geom["incl"] * np.pi / 180)
+
+    # deprojected radial coordinates
+    rr = np.ravel(np.hypot(xd, yd))
+
+    if bins is None:
+        # choose sensible bin size and range
+        step = np.hypot(icube.coords.cell_size, icube.coords.cell_size)
+        bins = np.arange(0.0, np.max((abs(xc.ravel()), abs(yc.ravel()))), step)
+
+    bin_counts, bin_edges = np.histogram(a=rr, bins=bins, weights=None)
+
+    # cumulative binned brightness in each annulus
+    Is, _ = np.histogram(a=rr, bins=bins, 
+                         weights=torch2npy(icube.sky_cube[chan]).ravel()
+                         )
+
+    # mask empty bins
+    mask = (bin_counts == 0)
+    Is = np.ma.masked_where(mask, Is)
+
+    # average binned brightness in each annulus
+    Is /= bin_counts
+
+    bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
+
+    return bin_centers, Is
+
+
+def radialV(fcube, geom, rescale_flux, chan=0, bins=None):
+    r"""
+    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 and weights 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:
+            :math:`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=[k\lambda]
+        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:[`k\lambda`]
+        Baselines corresponding to `u` and `v`
+    Vs : array, unit=[Jy] 
+        Visibility amplitudes at `q`
+
+    Notes
+    -----
+    This routine requires the `frank <https://github.com/discsim/frank>`_ package
+    """
+    from frank.geometry import apply_phase_shift, deproject
+
+    # projected model (u,v) points [k\lambda]
+    uu, vv = fcube.coords.sky_u_centers_2D, fcube.coords.sky_v_centers_2D
+
+    # visibilities 
+    V = torch2npy(fcube.ground_cube[chan]).ravel()
+
+    # phase-shift the visibilities
+    Vp = apply_phase_shift(uu.ravel() * 1e3, vv.ravel() * 1e3, V, geom["dRA"], 
+                           geom["dDec"], inverse=True)
+
+    # deproject the (u,v) points
+    up, vp, _ = deproject(uu.ravel() * 1e3, vv.ravel() * 1e3, geom["incl"], 
+                          geom["Omega"])
+
+    # if the source is optically thick, rescale the deprojected V(q)
+    if rescale_flux: 
+        Vp.real /= np.cos(geom["incl"] * np.pi / 180)
+
+    # convert back to [k\lambda]
+    up /= 1e3
+    vp /= 1e3
+
+    # deprojected baselines
+    qq = np.hypot(up, vp) 
+
+    if bins is None:
+        # choose sensible bin size and range
+        step = np.hypot(fcube.coords.du, fcube.coords.dv) / 2
+        bins = np.arange(0.0, max(qq), step)
+
+    bin_counts, bin_edges = np.histogram(a=qq, bins=bins, weights=None)
+
+    # cumulative binned visibility amplitude in each annulus
+    Vs, _ = np.histogram(a=qq, bins=bins, weights=Vp)
+
+    # mask empty bins
+    mask = (bin_counts == 0)
+    Vs = np.ma.masked_where(mask, Vs)
+
+    # average binned visibility amplitude in each annulus
+    Vs /= bin_counts
+
+    bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
+
+    return bin_centers, Vs

src/mpol/utils.py

@@ -302,7 +302,6 @@ def get_optimal_image_properties(image_width, u, v):
     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 

test/conftest.py

@@ -1,8 +1,9 @@
 import numpy as np
 import pytest
+import torch
 from astropy.utils.data import download_file
 
-from mpol import coordinates, gridding
+from mpol import coordinates, fourier, gridding, images
 from mpol.__init__ import zenodo_record
 
 # We need a fixture which provides mock visibilities of the sort we'd
@@ -116,6 +117,71 @@ def dataset_cont(mock_visibility_data_cont, coords):
     return averager.to_pytorch_dataset()
 
 
+@pytest.fixture(scope="session")
+def mock_1d_archive():
+    # use astropy routines to cache data
+    fname = download_file(
+        f"https://zenodo.org/record/{zenodo_record}/files/mock_disk_1d.npz",
+        cache=True,
+        pkgname="mpol",
+    )
+
+    return np.load(fname, allow_pickle=True)
+
+
+@pytest.fixture
+def mock_1d_image_model(mock_1d_archive):
+    m = mock_1d_archive
+    rtrue = m['rtrue']
+    itrue = m['itrue']
+    i2dtrue = m['i2dtrue']
+    xmax = ymax = m['xmax']
+    geom = m['geometry']
+    geom = geom[()]
+
+    coords = coordinates.GridCoords(cell_size=xmax * 2 / i2dtrue.shape[0], 
+                                    npix=i2dtrue.shape[0])
+
+    # the center of the array is already at the center of the image -->
+    # undo this as expected by input to ImageCube
+    i2dtrue = np.flip(np.fft.fftshift(i2dtrue), 1)
+
+    # pack the numpy image array into an ImageCube
+    packed_cube = np.broadcast_to(i2dtrue, (1, coords.npix, coords.npix)).copy()
+    packed_tensor = torch.from_numpy(packed_cube)
+    cube_true = images.ImageCube(coords=coords, nchan=1, cube=packed_tensor)
+
+    return rtrue, itrue, cube_true, xmax, ymax, geom
+
+
+@pytest.fixture
+def mock_1d_vis_model(mock_1d_archive):
+    m = mock_1d_archive
+    Vtrue = m['vis']
+    Vtrue_dep = m['vis_dep']
+    q_dep = m['baselines_dep']
+    geom = m['geometry']
+    geom = geom[()]
+
+    xmax = m['xmax']
+    i2dtrue = m['i2dtrue']
+
+    coords = coordinates.GridCoords(cell_size=xmax * 2 / i2dtrue.shape[0],
+                                    npix=i2dtrue.shape[0])
+
+    # create a FourierCube
+    packed_cube = np.broadcast_to(Vtrue, (1, len(Vtrue))).copy()
+    packed_tensor = torch.from_numpy(packed_cube)
+    cube_true = fourier.FourierCube(coords=coords)
+
+    # insert the vis tensor into the FourierCube ('vis' would typically be 
+    # populated by taking the FFT of an image)
+    cube_true.ground_cube = packed_tensor
+
+    return cube_true, Vtrue_dep, q_dep, geom
+
+
+
 @pytest.fixture
 def crossvalidation_products(mock_visibility_data):
     # test the crossvalidation with a smaller set of image / Fourier coordinates than normal,