Tighten and tidy API, docs page per module and better tests

* Optonal args are now keyword only
* Use of qantity aware type annotations
* Doc refaactor - page per submodule
* Fix bug in transform code
* Add fft equivalence test

changelog/58.breaking.rst

@@ -0,0 +1,2 @@
+Optional parameters are now keyword only for the :mod:`xrayvision.transform`, :mod:`xrayvision.imaging` and :mod:`xrayvision.visibility` modules.
+Remove ``natural`` keyword in favour of ``scheme`` keyword which can be either 'natural' or 'uniform'.

changelog/58.bugfix.rst

@@ -0,0 +1 @@
+Fix a bug where the x, y dimensions were not being treated consistently in :mod:`xrayvision.transform`.

changelog/58.docs.rst

@@ -0,0 +1 @@
+Add per module reference pages, switch to documenting types using type annotations.

docs/conf.py

@@ -44,6 +44,9 @@
 # Add any paths that contain templates here, relative to this directory.
 # templates_path = ['_templates']
 
+# Set automodapi to generate files inside the generated directory
+automodapi_toctreedirnm = "generated/api"
+
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
 # This pattern also affects html_static_path and html_extra_path.
@@ -61,14 +64,17 @@
 default_role = 'obj'
 
 # Disable having a separate return type row
-napoleon_use_rtype = False
+napoleon_use_rtype = True
 
 # Disable google style docstrings
 napoleon_google_docstring = False
 
 # until sphinx-gallery / sphinx is fixed https://github.com/sphinx-doc/sphinx/issues/12300
 suppress_warnings = ["config.cache"]
 
+autodoc_typehints = "description"
+autoclass_content = "init"
+
 # -- Options for intersphinx extension ---------------------------------------
 
 # Example configuration for intersphinx: refer to the Python standard library.

docs/reference/clean.rst

@@ -0,0 +1,8 @@
+.. _clean:
+
+Clean ('xrayvision.clean')
+**************************
+
+The ``clean`` submodule contains clean imaging methods.
+
+.. automodapi:: xrayvision.clean

docs/reference/imaging.rst

@@ -0,0 +1,8 @@
+.. _imaging:
+
+Imaging ('xrayvision.imaging')
+******************************
+
+The ``imaging`` submodule contains functions to make map and images from visibilities.
+
+.. automodapi:: xrayvision.imaging

docs/reference/index.rst

@@ -1,23 +1,15 @@
-.. _reference:
-
 *********
 Reference
 *********
 
-.. automodapi:: xrayvision.visibility
-
-.. automodapi:: xrayvision.imaging
-
-.. automodapi:: xrayvision.transform
-
-.. automodapi:: xrayvision.clean
-
-.. automodapi:: xrayvision.mem
-
-.. automodapi:: xrayvision.utils
-
-
 .. toctree::
-   :maxdepth: 2
+   :maxdepth: 1
+
+   clean
+   imaging
+   mem
+   transform
+   utils
+   visibility
 
    ../whatsnew/index

docs/reference/mem.rst

@@ -0,0 +1,8 @@
+.. _mem:
+
+MEM ('xrayvision.mem')
+**********************
+
+The ``mem`` submodule contains the Maximum Entropy methods.
+
+.. automodapi:: xrayvision.mem

docs/reference/transform.rst

@@ -0,0 +1,8 @@
+.. _transform:
+
+Transform ('xrayvision.transform')
+**********************************
+
+The ``transform`` submodule forward and reverse transforms.
+
+.. automodapi:: xrayvision.transform

docs/reference/utils.rst

@@ -0,0 +1,8 @@
+.. _utils:
+
+Utils ('xrayvision.utils')
+**************************
+
+The ``utils`` submodule contains utility functions.
+
+.. automodapi:: xrayvision.utils

docs/reference/visibility.rst

@@ -0,0 +1,8 @@
+.. _visibility:
+
+Visibility ('xrayvision.visibility')
+************************************
+
+The ``visibility`` submodule contains generic visibility classes.
+
+.. automodapi:: xrayvision.visibility

docs/tutorials/introduction.rst

@@ -41,10 +41,10 @@ measured visibilities in f can the original map a be recovered?
 
     data = make_data()
 
-    full_uv = transform.generate_uv(65)
+    full_uv = transform.generate_uv(65*u.pix)
 
-    uu = transform.generate_uv(33)
-    vv = transform.generate_uv(33)
+    uu = transform.generate_uv(33*u.pix)
+    vv = transform.generate_uv(33*u.pix)
 
     uu, vv = np.meshgrid(uu, vv)
 
@@ -53,7 +53,7 @@ measured visibilities in f can the original map a be recovered?
 
     full_vis = transform.dft_map(data, u=uv[0,:], v=uv[1,:])
 
-    res = transform.idft_map(full_vis, u=uv[0,:], v=uv[1,:], shape=(33, 33))
+    res = transform.idft_map(full_vis, u=uv[0,:], v=uv[1,:], weights=1/33**2, shape=(33, 33)*u.pix)
     # assert np.allclose(data, res)
 
     # Generate log spaced radial u, v sampeling
@@ -77,9 +77,10 @@ measured visibilities in f can the original map a be recovered?
     sub_vis = transform.dft_map(data, u=sub_uv[0,:], v=sub_uv[1,:])
 
     psf1 = transform.idft_map(np.full(sub_vis.size, 1), u=sub_uv[0,:], v=sub_uv[1,:],
-                              shape=(65, 65))
+                              weights=1/sub_vis.size, shape=(65, 65)*u.pix)
 
-    sub_res = transform.idft_map(sub_vis, u=sub_uv[0,:], v=sub_uv[1,:], shape=(65, 65))
+    sub_res = transform.idft_map(sub_vis, u=sub_uv[0,:], v=sub_uv[1,:],
+                                 weights=1/sub_vis.size, shape=(65, 65)*u.pix)
 
     xp = np.round(x * 33 + 33/2 - 0.5 + 16).astype(int)
     yp = np.round(y * 33 + 33/2 - 0.5 + 16).astype(int)

examples/rhessi.py

@@ -58,15 +58,15 @@
 # Lets have a look at the point spread function (PSF) or dirty beam
 
 psf_map = vis_psf_map(vis, shape=(101, 101)*apu.pixel,
-                      pixel_size=1.5*apu.arcsec,
-                      natural=False)
+                      pixel_size=1.5*apu.arcsec/apu.pixel,
+                      scheme='uniform')
 
 ###############################################################################
 # We can now make an image using the back projection algorithm essentially and
 # inverse Fourier transform of the visibilities.
 
 backproj_map = vis_to_map(vis, shape=[101, 101]*apu.pixel,
-                          pixel_size=1.5*apu.arcsec)
+                          pixel_size=1.5*apu.arcsec/apu.pix)
 
 ###############################################################################
 # Back projection contain many artifact due to the incomplete sampling of the u-v
@@ -78,15 +78,14 @@
 # vis_59 = Visibility(vis=vis_data_59['obsvis']*apu.Unit('ph/cm*s'), u=vis_data_59['u']/apu.arcsec,
 #                     v=vis_data_59['v']/apu.arcsec, offset=vis_data_59['xyoffset'][0]*apu.arcsec)
 
-clean_map, model_map, residual_map = vis_clean(vis, shape=[101, 101]*apu.pixel,
-                                               pixel=[1.5, 1.5]*apu.arcsec,
-                                               clean_beam_width=10*apu.arcsec,
-                                               niter=100)
+clean_map, model_map, residual_map = vis_clean(vis, shape=[101, 101] * apu.pixel,
+                                               pixel_size=[1.5, 1.5] * apu.arcsec / apu.pix,
+                                               clean_beam_width=10 * apu.arcsec, niter=100)
 
 ###############################################################################
 # MEM
 
-mem_map = mem(vis, shape=[129, 129]*apu.pixel, pixel=[2, 2]*apu.arcsec)
+mem_map = mem(vis, shape=[129, 129]*apu.pixel, pixel=[2, 2]*apu.arcsec/apu.pix)
 mem_map.plot()
 
 

examples/stix.py

@@ -25,35 +25,36 @@
 stix_data = pickle.load(urllib.request.urlopen("https://pub099.cs.technik.fhnw.ch/demo/stix_vis.pkl"))
 
 time_range, energy_range, offset, stix_vis = stix_data
+stix_vis.phase_centre = [0, 0] * apu.arcsec
 stix_vis.offset = offset
 
 ###############################################################################
 # Lets have a look at the point spread function (PSF) or dirty beam
 
 psf_map = vis_psf_map(stix_vis, shape=(129, 129)*apu.pixel,
-                      pixel_size=2*apu.arcsec,
-                      natural=False)
+                      pixel_size=2*apu.arcsec/apu.pix,
+                      scheme='uniform')
 psf_map.plot()
 
 ###############################################################################
 # Back projection
 
 backproj_map = vis_to_map(stix_vis, shape=(129, 129)*apu.pixel,
-                          pixel_size=2*apu.arcsec, natural=False)
+                          pixel_size=2*apu.arcsec/apu.pix, scheme='uniform')
 backproj_map.plot()
 
 ###############################################################################
 # Clean
 
-clean_map, model_map, resid_map = vis_clean(stix_vis, shape=[129, 129]*apu.pixel,
-                                            pixel=[2, 2]*apu.arcsec, niter=100,
-                                            clean_beam_width=20*apu.arcsec)
+clean_map, model_map, resid_map = vis_clean(stix_vis, shape=[129, 129] * apu.pixel,
+                                            pixel_size=[2, 2] * apu.arcsec / apu.pix, clean_beam_width=20 * apu.arcsec,
+                                            niter=100)
 clean_map.plot()
 
 ###############################################################################
 # MEM
 
-mem_map = mem(stix_vis, shape=[129, 129]*apu.pixel, pixel=[2, 2]*apu.arcsec)
+mem_map = mem(stix_vis, shape=[129, 129]*apu.pixel, pixel=[2, 2]*apu.arcsec/apu.pix)
 mem_map.plot()
 
 ###############################################################################

setup.cfg

@@ -96,6 +96,9 @@ ignore-words-list =
     process,
     technik
 
+[mypy]
+disable_error_code = import-untyped
+
 [coverage:run]
 omit =
   xrayvision/_sunpy_init*

xrayvision/__init__.py

@@ -1,11 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
 try:
-    from xrayvision.version import __version__
+    from xrayvision.version import __version__ # type: ignore
 except ImportError:
     __version__ = "unknown"
-__all__ = []
-
-from pkg_resources import resource_filename
-
-SAMPLE_RHESSI_VISIBILITIES = resource_filename('xrayvision', 'data/hsi_visibili_20131028_0156_20131028_0200_6_12.fits')

xrayvision/clean.py

@@ -11,7 +11,7 @@
 import numpy as np
 from astropy.convolution import Gaussian2DKernel
 from scipy import signal
-from scipy.ndimage.interpolation import shift
+from scipy.ndimage import shift
 from sunpy.map.map_factory import Map
 
 from xrayvision.imaging import vis_psf_image, vis_to_map
@@ -58,7 +58,7 @@
     """
 
 
-def clean(dirty_map, dirty_beam, pixel=None, clean_beam_width=4.0,
+def clean(dirty_map, dirty_beam, pixel_size=None, clean_beam_width=4.0,
           gain=0.1, thres=0.01, niter=5000):
     r"""
     Clean the image using Hogbom's original method.
@@ -74,15 +74,15 @@ def clean(dirty_map, dirty_beam, pixel=None, clean_beam_width=4.0,
         The dirty map to be cleaned 2D
     dirty_beam : `numpy.ndarray`
         The dirty beam or point spread function (PSF) 2D must
-    pixel :
+    pixel_size :
         Size of a pixel
     """
     # Ensure both beam and map are even/odd on same axes
     # if not [x % 2 == 0 for x in dirty_map.shape] == [x % 2 == 0 for x in dirty_beam.shape]:
     #     raise ValueError('')
     pad = [0 if x % 2 == 0 else 1 for x in dirty_map.shape]
 
-    # Assume beam, map center is in middle
+    # Assume beam, map phase_centre is in middle
     beam_center = (dirty_beam.shape[0] - 1)/2.0, (dirty_beam.shape[1] - 1)/2.0
     map_center = (dirty_map.shape[0] - 1)/2.0, (dirty_map.shape[1] - 1)/2.0
 
@@ -132,28 +132,28 @@ def clean(dirty_map, dirty_beam, pixel=None, clean_beam_width=4.0,
 
     if clean_beam_width != 0.0:
         # Convert from FWHM to StDev   FWHM = sigma*(8ln2)**0.5 = 2.3548200450309493
-        x_stdev = (clean_beam_width/pixel[0] / 2.3548200450309493).value
-        y_stdev = (clean_beam_width/pixel[1] / 2.3548200450309493).value
+        x_stdev = (clean_beam_width / pixel_size[0] / 2.3548200450309493).value
+        y_stdev = (clean_beam_width / pixel_size[1] / 2.3548200450309493).value
         clean_beam = Gaussian2DKernel(x_stdev, y_stdev, x_size=dirty_beam.shape[1],
                                       y_size=dirty_beam.shape[0]).array
         # Normalise beam
         clean_beam = clean_beam / clean_beam.max()
 
         # Convolve clean beam with model and scale
         clean_map = (signal.convolve2d(model, clean_beam/clean_beam.sum(), mode='same')
-                     / (pixel[0]*pixel[1]))
+                     / (pixel_size[0] * pixel_size[1]))
 
         # Scale residual map with model and scale
-        dirty_map = dirty_map / clean_beam.sum() / (pixel[0] * pixel[1])
+        dirty_map = dirty_map / clean_beam.sum() / (pixel_size[0] * pixel_size[1])
         return clean_map+dirty_map, model, dirty_map
 
     return model+dirty_map, model, dirty_map
 
 
-clean.__doc__ += __common_clean_doc__
+clean.__doc__ += __common_clean_doc__  # type: ignore
 
 
-def vis_clean(vis, shape, pixel, clean_beam_width=4.0, niter=5000, map=True, gain=0.1, **kwargs):
+def vis_clean(vis, shape, pixel_size, clean_beam_width=4.0, niter=5000, map=True, gain=0.1, **kwargs):
     r"""
     Clean the visibilities using Hogbom's original method.
 
@@ -165,24 +165,24 @@ def vis_clean(vis, shape, pixel, clean_beam_width=4.0, niter=5000, map=True, gai
         The visibilities to clean
     shape :
         Size of map
-    pixel :
+    pixel_size :
         Size of pixel
     map : `boolean` optional
         Return an `sunpy.map.Map` by default or data only if `False`
     """
 
-    dirty_map = vis_to_map(vis, shape=shape, pixel_size=pixel, **kwargs)
+    dirty_map = vis_to_map(vis, shape=shape, pixel_size=pixel_size, **kwargs)
     dirty_beam_shape = [x.value*3 + 1 if x.value*3 % 2 == 0 else x.value*3 for x in shape]*shape.unit
-    dirty_beam = vis_psf_image(vis, shape=dirty_beam_shape, pixel_size=pixel, **kwargs)
-    clean_map, model, residual = clean(dirty_map.data, dirty_beam.value, pixel=pixel, gain=gain,
-                                       clean_beam_width=clean_beam_width, niter=niter)
+    dirty_beam = vis_psf_image(vis, shape=dirty_beam_shape, pixel_size=pixel_size, **kwargs)
+    clean_map, model, residual = clean(dirty_map.data, dirty_beam.value, pixel_size=pixel_size,
+                                       clean_beam_width=clean_beam_width, gain=gain, niter=niter)
     if not map:
         return clean_map, model, residual
 
     return [Map((data, dirty_map.meta)) for data in (clean_map, model, residual)]
 
 
-vis_clean.__doc__ += __common_clean_doc__
+vis_clean.__doc__ += __common_clean_doc__  # type: ignore
 
 
 __common_ms_clean_doc__ = r"""
@@ -332,7 +332,7 @@ def ms_clean(dirty_map, dirty_beam, pixel, scales=None,
     else:
         logger.info("Max iterations reached")
 
-    # Convolve model with clean beam  B_G * I^M
+    # Convolve model with clean beam B_G * I^M
     if clean_beam_width != 0.0:
         x_stdev = (clean_beam_width/pixel[0] / (2.0 * np.sqrt(2.0 * np.log(2.0)))).value
         y_stdev = (clean_beam_width/pixel[1] / (2.0 * np.sqrt(2.0 * np.log(2.0)))).value
@@ -352,7 +352,7 @@ def ms_clean(dirty_map, dirty_beam, pixel, scales=None,
     return model, scaled_residuals.sum(axis=2)
 
 
-ms_clean.__doc__ += __common_ms_clean_doc__
+ms_clean.__doc__ += __common_ms_clean_doc__  # type: ignore
 
 
 def vis_ms_clean(vis, shape, pixel, scales=None, clean_beam_width=4.0,