First example of defining C++/Python bindings interface

Define ``cpp_image_visibilities``

Author: timstaley

src/fastimgproto/bindings/__init__.py

@@ -0,0 +1,4 @@
+"""
+Define the Python wrapper-functions which provide an interface to the C++
+implementations.
+"""

src/fastimgproto/bindings/imager.py

@@ -0,0 +1,172 @@
+import astropy.units as u
+import numpy as np
+
+from fastimgproto.gridder.gridder import convolve_to_grid
+from fastimgproto.imager import fft_to_image_plane
+from fastimgproto.gridder.conv_funcs import GaussianSinc
+
+
+class CppKernelFuncs(object):
+    gauss_sinc = 'gauss_sinc'
+
+
+def cpp_image_visibilities(vis, uvw_lambda,
+                           image_size, cell_size,
+                           kernel_func_name=CppKernelFuncs.gauss_sinc,
+                           kernel_trunc_radius=3.0,
+                           kernel_support=3,
+                           kernel_oversampling=None,
+                           normalize=True):
+    """
+    Convenience wrapper over _cpp_image_visibilities.
+
+    Functionality largely mirrors
+    :func:`fastimgproto.imager.image_visibilities`, but the key difference is
+    that instead of passing a callable kernel-function, you must choose
+    ``kernel_func_name`` from a limited selection of kernel-functions
+    implemented in the C++ code. Currently, choices are limited to:
+
+        - ``gauss_sinc``
+
+
+    Performs the following tasks before handing over to C++ bindings:
+    - Checks arguments are of correct type / units
+    - Converts uvw-array from wavelength (lambda) units to pixel units
+
+    Args:
+        vis (numpy.ndarray): Complex visibilities.
+            1d array, shape: `(n_vis,)`.
+        uvw_lambda (numpy.ndarray): UVW-coordinates of visibilities. Units are
+            multiples of wavelength.
+            2d array of ``np.float_``, shape: ``(n_vis, 3)``.
+            Assumed ordering is u,v,w i.e. ``u,v,w = uvw[idx]``
+        image_size (astropy.units.Quantity): Width of the image in pixels.
+            e.g. ``1024 * u.pixel``.
+            NB we assume the pixel ``[image_size//2,image_size//2]``
+            corresponds to the origin in UV-space.
+        cell_size (astropy.units.Quantity): Angular-width of a synthesized pixel
+            in the image to be created, e.g. ``3.5 * u.arcsecond``.
+        kernel_func_name (str): Choice of kernel function from limited C++ selection.
+        kernel_trunc_radius (float): Truncation radius of the kernel to be used.
+        kernel_support (int): Defines the 'radius' of the bounding box within
+            which convolution takes place. `Box width in pixels = 2*support+1`.
+            (The central pixel is the one nearest to the UV co-ordinates.)
+            (This is sometimes known as the 'half-support')
+        kernel_oversampling (int): (Or None). Controls kernel-generation,
+            see :func:`fastimgproto.gridder.gridder.convolve_to_grid` for
+            details.
+        normalize (bool): Whether or not the returned image and beam
+            should be normalized such that the beam peaks at a value of
+            1.0 Jansky. You normally want this to be true, but it may be
+            interesting to check the raw values for debugging purposes.
+
+    Returns:
+        tuple: (image, beam)
+            Tuple of ndarrays representing the image map and beam model.
+            These are 2d arrays of same dtype as ``vis``,
+            (typically ``np._complex``),  shape ``(image_size, image_size)``.
+            Note numpy style index-order, i.e. access like ``image[y,x]``.
+
+    """
+    if kernel_func_name not in (CppKernelFuncs.gauss_sinc,):
+        raise ValueError(
+            "kernel function of type {} not recognised".format(
+                kernel_func_name))
+
+    image_size = image_size.to(u.pix)
+    # Size of a UV-grid pixel, in multiples of wavelength (lambda):
+    grid_pixel_width_lambda = 1.0 / (cell_size.to(u.rad) * image_size)
+    uvw_in_pixels = (uvw_lambda / grid_pixel_width_lambda).value
+    uv_in_pixels = uvw_in_pixels[:, :2]
+
+    # subroutine = _cpp_image_visibilities
+    subroutine = _python_image_visibilities
+    (image, beam) = _python_image_visibilities(
+        vis=vis,
+        uv_pixels=uv_in_pixels,
+        image_size=int(image_size.value),
+        kernel_func_name=kernel_func_name,
+        kernel_trunc_radius=kernel_trunc_radius,
+        kernel_support=kernel_support,
+        kernel_oversampling=kernel_oversampling,
+        normalize=normalize,
+    )
+
+    return image, beam
+
+
+def _cpp_image_visibilities(vis,
+                            uv_pixels,
+                            image_size,
+                            kernel_func_name,
+                            kernel_trunc_radius,
+                            kernel_support,
+                            kernel_oversampling,
+                            normalize=True
+                            ):
+    pass
+    # C++ Bindings here
+
+
+def _python_image_visibilities(vis,
+                               uv_pixels,
+                               image_size,
+                               kernel_func_name,
+                               kernel_trunc_radius,
+                               kernel_support,
+                               kernel_oversampling,
+                               normalize=True
+                               ):
+    """
+    Equivalent Python code for validation of _cpp_image_visibilities
+
+    Args:
+        vis (numpy.ndarray): Complex visibilities.
+            1d array, shape: `(n_vis,)`.
+        uv_pixels (numpy.ndarray): UV-coordinates of visibilities. Units are
+            pixel-widths relative to the grid being sampled onto.
+            2d array of ``np.float_``, shape: ``(n_vis, 2)``.
+            Assumed ordering is u,v i.e. ``u,v = uv[idx]``
+        image_size (int): Width of the image in pixels.
+            NB we assume the pixel ``[image_size//2,image_size//2]``
+            corresponds to the origin in UV-space.
+        kernel_func_name (str): Choice of kernel function from limited C++ selection.
+        kernel_trunc_radius (float): Truncation radius of the kernel to be used.
+        kernel_support (int): Defines the 'radius' of the bounding box within
+            which convolution takes place. `Box width in pixels = 2*support+1`.
+            (The central pixel is the one nearest to the UV co-ordinates.)
+            (This is sometimes known as the 'half-support')
+        kernel_oversampling (int): (Or None). Controls kernel-generation,
+            see :func:`fastimgproto.gridder.gridder.convolve_to_grid` for
+            details.
+        normalize (bool): Whether or not the returned image and beam
+            should be normalized such that the beam peaks at a value of
+            1.0 Jansky. You normally want this to be true, but it may be
+            interesting to check the raw values for debugging purposes.
+
+    Returns:
+        tuple: (image, beam)
+            Tuple of ndarrays representing the image map and beam model.
+            These are 2d arrays of same dtype as ``vis``,
+            (typically ``np._complex``),  shape ``(image_size, image_size)``.
+            Note numpy style index-order, i.e. access like ``image[y,x]``.
+
+    """
+    assert kernel_func_name == CppKernelFuncs.gauss_sinc
+    kernel_func = GaussianSinc(trunc=kernel_trunc_radius)
+
+    vis_grid, sample_grid = convolve_to_grid(kernel_func,
+                                             support=kernel_support,
+                                             image_size=image_size,
+                                             uv=uv_pixels,
+                                             vis=vis,
+                                             oversampling=kernel_oversampling
+                                             )
+    image = fft_to_image_plane(vis_grid)
+    beam = fft_to_image_plane(sample_grid)
+    if normalize:
+        beam_max = np.max(np.real(beam))
+        beam /= beam_max
+        image /= beam_max
+
+    return (image, beam)

src/fastimgproto/imager.py

@@ -12,16 +12,57 @@ def image_visibilities(vis, uvw_lambda,
                        kernel_func, kernel_support,
                        kernel_oversampling,
                        normalize=True):
+    """
+    Args:
+        vis (numpy.ndarray): Complex visibilities.
+            1d array, shape: `(n_vis,)`.
+        uvw_lambda (numpy.ndarray): UVW-coordinates of visibilities. Units are
+            multiples of wavelength.
+            2d array of ``np.float_``, shape: ``(n_vis, 3)``.
+            Assumed ordering is u,v,w i.e. ``u,v,w = uvw[idx]``
+        image_size (astropy.units.Quantity): Width of the image in pixels.
+            e.g. ``1024 * u.pixel``.
+            NB we assume the pixel ``[image_size//2,image_size//2]``
+            corresponds to the origin in UV-space.
+        cell_size (astropy.units.Quantity): Angular-width of a synthesized pixel
+            in the image to be created, e.g. ``3.5 * u.arcsecond``.
+        kernel_func (callable): Callable object,
+            (e.g. :class:`.conv_funcs.Pillbox`,)
+            that returns a convolution
+            co-efficient for a given distance in pixel-widths.
+        kernel_support (int): Defines the 'radius' of the bounding box within
+            which convolution takes place. `Box width in pixels = 2*support+1`.
+            (The central pixel is the one nearest to the UV co-ordinates.)
+            (This is sometimes known as the 'half-support')
+        kernel_oversampling (int): (Or None). Controls kernel-generation,
+            see :func:`fastimgproto.gridder.gridder.convolve_to_grid` for
+            details.
+        normalize (bool): Whether or not the returned image and beam
+            should be normalized such that the beam peaks at a value of
+            1.0 Jansky. You normally want this to be true, but it may be
+            interesting to check the raw values for debugging purposes.
+
+    Returns:
+        tuple: (image, beam)
+            Tuple of ndarrays representing the image map and beam model.
+            These are 2d arrays of same dtype as ``vis``,
+            (typically ``np._complex``),  shape ``(image_size, image_size)``.
+            Note numpy style index-order, i.e. access like ``image[y,x]``.
+    """
 
     image_size = image_size.to(u.pix)
+    image_size_int = int(image_size.value)
+    if image_size_int != image_size.value:
+        raise ValueError("Please supply an integer-valued image size")
+
     # Size of a UV-grid pixel, in multiples of wavelength (lambda):
     grid_pixel_width_lambda = 1.0 / (cell_size.to(u.rad) * image_size)
     uvw_in_pixels = (uvw_lambda / grid_pixel_width_lambda).value
 
     uv_in_pixels = uvw_in_pixels[:, :2]
     vis_grid, sample_grid = convolve_to_grid(kernel_func,
                                              support=kernel_support,
-                                             image_size=int(image_size.value),
+                                             image_size=image_size_int,
                                              uv=uv_in_pixels,
                                              vis=vis,
                                              oversampling=kernel_oversampling

tests/test_cpp_bindings/__init__.py

@@ -0,0 +1,46 @@
+import numpy as np
+import astropy.units as u
+from fastimgproto.fixtures.data import simple_vis_npz_filepath
+from fastimgproto.imager import image_visibilities
+from fastimgproto.bindings.imager import cpp_image_visibilities, CppKernelFuncs
+from fastimgproto.gridder.conv_funcs import GaussianSinc
+
+
+def test_cpp_imager_results_parity():
+    with open(simple_vis_npz_filepath, 'rb') as f:
+        npz_data_dict = np.load(f)
+        uvw_lambda = npz_data_dict['uvw_lambda']
+        vis = npz_data_dict['vis']
+
+    trunc = 3.0
+    support = 3
+
+    kernel_func = GaussianSinc(trunc=trunc)
+    kernel_func_name = CppKernelFuncs.gauss_sinc
+
+    image_size = 1024 * u.pix
+    cell_size = 3. * u.arcsec
+    py_img, py_beam = image_visibilities(
+        vis=vis,
+        uvw_lambda=uvw_lambda,
+        image_size= image_size,
+        cell_size=cell_size,
+        kernel_func = kernel_func,
+        kernel_support=support,
+        kernel_oversampling=None,
+        normalize=True
+    )
+
+    cpp_img, cpp_beam = cpp_image_visibilities(
+        vis=vis,
+        uvw_lambda=uvw_lambda,
+        image_size=image_size,
+        cell_size=cell_size,
+        kernel_func_name=kernel_func_name,
+        kernel_trunc_radius=trunc,
+        kernel_support=support,
+        kernel_oversampling=None,
+        normalize=True
+    )
+
+    assert (np.abs(cpp_img - py_img)< 1E-5).all()

tests/test_cpp_bindings/test_imager_bindings.py

@@ -27,7 +27,7 @@ def test_fractional_coord_to_oversampled_index_math():
     assert (calculate_oversampled_kernel_indices(
         subpixel_coord=subpix_offset, oversampling=oversampling) == 3).all()
 
-    # OK, now demonstrate values with oversampling of 0.5, which has easy
+    # OK, now demonstrate values with oversampling of 5, which has easy
     # numbers to calculate since 1/5 = 0.2
     oversampling = 5
     io_pairs = np.array([