Bug fixes and doc updates (#40)

* Fix lint issues on long lines in docstrings

* Fix various bugs

* Bug fixes (#38)

* Fix #29 add unit decorator to init method

* Fix #30 sign change on u, v coords on import

* Fix lint issues on long lines in docstrings

* Change normalisation of idft to no of visibilities from number of pixels

* Temporary fix for to_map for RHESSI

* Update .rtd-environment.yml so doc build works on RTD

* Update tutorial and readme

* Update tutorial add images

* Add RHESSI time parsing

* More travis updates

* Update README fix badges

* Bug fixes and documentation updates

README.rst

@@ -1,7 +1,7 @@
 XRAYVISION  - X-RAY VIsibility Synthesis ImagiNg
 ================================================
 
-|Powered By| |Build Status| |Doc Status|
+|Powered By| |Build Status| |Doc Status| |Python Versions|
 
 XRAYVISION is an open-source Python library for Fourier or synthesis imaging of X-Rays. The most
 common usage of this technique is radio interferometry however there have been a number of solar
@@ -48,12 +48,16 @@ follow our `Code of Conduct`_.
     :target: http://www.sunpy.org
     :alt: Powered by SunPy Badge
 
-.. |Build Status| image:: https://travis-ci.org/samaloney/xrayvision.svg?branch=master
+.. |Build Status| image:: https://travis-ci.org/sunpy/xrayvision.svg?branch=master
     :target: https://travis-ci.org/sunpy/xrayvision
     :alt: Travis-CI build status
 
-.. |Doc Status|  image:: https://readthedocs.org/projects/xrayvision/badge/?version=latest
+.. |Doc Status| image:: https://readthedocs.org/projects/xrayvision/badge/?version=latest
     :target: http://xrayvision.readthedocs.io/en/latest/?badge=latest
     :alt: Documentation Status
 
+.. |Python Versions| image:: https://img.shields.io/badge/python-3.6-blue.svg
+    :target: https://www.python.org/downloads/release/python-360/
+    :alt: Python Versions
+
 .. _Code of Conduct: http://docs.sunpy.org/en/stable/coc.html

docs/xrayvision/introduction.rst

@@ -1,17 +1,128 @@
 Introduction
 ============
 
+Below is a figure inspired by `Dale Gary`_, which summarises the problem XRAYVISION tries to solve.
+The top row shows a source map, the point spread function (PSF) or dirty beam and the convolution of
+two, the dirty map (left to right). The bottom row shows the corresponding visibilities,
+notice the convolution is replaced by multiplication in frequency space. The problem is given the
+measured visibilities in f can the original map a be recovered?
+
+.. plot::
+
+    import numpy as np
+
+    from astropy import units as u
+    from astropy.convolution import convolve, Gaussian2DKernel
+    from matplotlib import pyplot as plt
+    from matplotlib.colors import LogNorm
+
+    from xrayvision import transform
+    from xrayvision import clean
+
+
+    def make_data(noisy=False):
+        g1 = Gaussian2DKernel(stddev=2, x_size=65, y_size=65).array
+        g2 = Gaussian2DKernel(stddev=6, x_size=65, y_size=65).array
+
+        temp = np.zeros((65, 65))
+        temp[50,50] = 400.0
+
+        data = convolve(temp, g1)
+
+        temp = np.zeros((65, 65))
+        temp[20,20] =2600.0
+
+        data = data + convolve(temp, g2);
+
+        if noisy:
+            data = data + np.random.normal(loc=0.0, scale=0.005, size=(65, 65));
+
+        return data
+
+    data = make_data()
+
+    full_uv = transform.generate_uv(65)
+
+    uu = transform.generate_uv(33)
+    vv = transform.generate_uv(33)
+
+    uu, vv = np.meshgrid(uu, vv)
+
+    # uv coordinates require unit
+    uv = np.array([uu, vv]).reshape(2, 33**2)/u.arcsec
+
+    full_vis = transform.dft_map(data, uv)
+
+    res = transform.idft_map(full_vis, (33, 33), uv)
+    # assert np.allclose(data, res)
+
+    # Generate log spaced radial u, v sampeling
+    half_log_space = np.logspace(np.log10(np.abs(uv[uv != 0]).value.min()), np.log10(np.abs(uv.value.max())), 10)
+
+    theta = np.linspace(0, 2*np.pi, 31)
+    theta = theta[np.newaxis,:]
+    theta = np.repeat(theta, 10, axis=0)
+
+    r = half_log_space
+    r = r[:, np.newaxis]
+    r = np.repeat(r, 31, axis=1)
+
+    x = r * np.sin(theta)
+    y = r * np.cos(theta)
+
+
+    sub_uv = np.vstack([x.flatten(), y.flatten()])
+    sub_uv = np.hstack([sub_uv, np.zeros((2,1))])/u.arcsec
+
+    sub_vis = transform.dft_map(data, sub_uv)
+
+    psf1 = transform.idft_map(np.full(sub_vis.size, 1), (65, 65), sub_uv)
+
+    sub_res = transform.idft_map(sub_vis, (65, 65), sub_uv)
+
+    xp = np.round(x * 33 + 33/2 - 0.5 + 16).astype(int)
+    yp = np.round(y * 33 + 33/2 - 0.5 + 16).astype(int)
+
+    sv = np.zeros((65, 65))
+    sv[:,:] = 0
+    sv[xp.flatten(), yp.flatten()] = 1
+
+    v = np.pad(np.abs(full_vis.reshape(33, 33))**2, ((16, 16),(16,16)), mode='constant', constant_values=0.11)
+
+    s_v = v*sv
+    s_v[s_v == 0] = 0.11
+
+    def im_plot(axis, data, *, text, label, **imshow_kwargs):
+        axis.text(0.5, 0.85, text, fontsize=15, color='white', transform=axis.transAxes,
+        horizontalalignment='center', verticalalignment='bottom', fontweight='bold')
+        axis.text(0.05, 0.9, label, fontsize=14, color='white', transform=axis.transAxes)
+        axis.imshow(data, origin='lower', **imshow_kwargs)
+        axis.axis('off')
+
+    f, (r1, r2) = plt.subplots(2, 3, figsize=(12,8))
+
+    im_plot(r1[0], data, text=r'$I(l, m)$', label='a )')
+    im_plot(r1[1], psf1, text=r'$B(l, m)$', label='b )')
+    im_plot(r1[2], sub_res, text=r'$I(l, m) *B(l, m)$', label='c )')
+    im_plot(r2[0], v, text=r'$V(u, v)$', norm=LogNorm(0.1), label='d )')
+    im_plot(r2[1], np.ones((65,65)), text=r'$S(u, v)$', label='e )', extent=(-1, 1, -1, 1))
+    r2[1].plot(x.flatten(), y.flatten(), 'w.', ms=2.5)
+    im_plot(r2[2], s_v, text=r'$S(u,v)V(u, v)$', label='f )', extent=(-1, 1, -1, 1), norm=LogNorm(0.1))
+
+    f.subplots_adjust(hspace=0.05, wspace=0.025)
+    plt.show()
+
 Theory
 ------
 Synthesis imaging relies upon describing the amplitude of some quantity on the sky (radio flux or
 x-ray photon flux) in terms of complex visibilities as:
 
-.. math:: I(x,y) = \int^{\infty}_{-\infty}\int^{\infty}_{-\infty}V(u, v)e^{-2 i \pi(ux+vy}) du dv
+.. math:: I(l,m) = \int^{\infty}_{-\infty}\int^{\infty}_{-\infty}V(u, v)e^{-2 i \pi(ul+vm}) du dv
    :label: ifft
 
 and the complex visibilties are given by:
 
-.. math:: V(u,v) = \int^{\infty}_{-\infty}\int^{\infty}_{-\infty}I(x, y)e^{2 i \pi(ux+vy}) du dv
+.. math:: V(u,v) = \int^{\infty}_{-\infty}\int^{\infty}_{-\infty}I(l, m)e^{2 i \pi(ux+vy}) dl dm
    :label: fft
 
 In the case where the :math:`u, v` plane is fully sampled the amplitude can be retrieved by simple
@@ -31,7 +142,7 @@ applying the convolution theorem
 where :math:`B = \mathscr{F}^{-1} S` is the point spread function (PSF) also known as the dirty
 beam given by
 
-.. math:: B(x, y) = \sum_{i} e^{-2 i \pi(u_{i}x+v_{i}y)}w_{i}.
+.. math:: B(l, m) = \sum_{i} e^{-2 i \pi(u_{i}l+v_{i}m)}w_{i}.
 
 So the problem is to deconvolve the effects of the PSF or diry beam :math:`B` from the dirty
 image :math:`I^{D}` to obtain the true image :math:`I`.
@@ -41,6 +152,8 @@ Implementation
 In reality the integrals above must be turned into summations over finite coordinates so :eq:`ifft`
 can be written as
 
-.. math:: I(x_i, y_j) = \sum_{k=0}^{N} e^{2 \pi i ( x_i u_k + y_i v_k)}
+.. math:: I(l_i, m_j) = \sum_{k=0}^{N} e^{2 \pi i ( l_i u_k + m_i v_k)}
+
+where :math:`x_i`
 
-where :math:`x_i`
+.. _Dale Gary: https://web.njit.edu/~gary/728/Lecture6.html

docs/xrayvision/tutorial.rst

@@ -12,6 +12,7 @@ Begin by importing the necessary libraries
     import numpy as np
 
     from astropy import units as u
+    from sunpy.map import Map
     from xrayvision.visibility import RHESSIVisibility
     from xrayvision import clean, SAMPLE_RHESSI_VISIBILITIES
 
@@ -22,16 +23,27 @@ as visibilities in a fits file using the RHESSI IDL software.
 
     rhessi_vis = RHESSIVisibility.from_fits_file(SAMPLE_RHESSI_VISIBILITIES)
 
-At this stage we can create a map with a back projection or inverse transform also known as the dirty map
-of the visibility data. The size of the map and the physical pixel size can be specified as
-to the `to_map` function of the visibility.
+At this stage we can create a map with a back projection or inverse transform also known as the
+dirty map of the visibility data. The size of the map and the physical pixel size can be specified
+as to the :meth:`~xrayvision.visibility.Visibility.to_map` function of the visibility.
 
 .. code:: python
 
     rhessi_map = rhessi_vis.to_map(shape=(65, 65), pixel_size=[4., 4.] * u.arcsec)
     rhessi_map.peek()
 
-.. image:: tutorial/bproj.png
+.. plot::
+
+    import numpy as np
+
+    from astropy import units as u
+    from xrayvision.visibility import RHESSIVisibility
+    from xrayvision import clean, SAMPLE_RHESSI_VISIBILITIES
+
+    rhessi_vis = RHESSIVisibility.from_fits_file(SAMPLE_RHESSI_VISIBILITIES)
+
+    rhessi_map = rhessi_vis.to_map(shape=(65, 65), pixel_size=[4., 4.] * u.arcsec)
+    rhessi_map.plot()
 
 The artifacts due to the under sampling of the u, v plane are clear. The main goal
 of this library is to provide a number image reconstruction methods such as CLEAN.
@@ -49,4 +61,64 @@ inputs.
     rhessi_map.data[:] = clean_map
     rhessi_map.peek()
 
-.. image:: tutorial/clean.png
+.. plot::
+
+    import numpy as np
+
+    from astropy import units as u
+    from xrayvision.visibility import RHESSIVisibility
+    from xrayvision import clean, SAMPLE_RHESSI_VISIBILITIES
+
+    rhessi_vis = RHESSIVisibility.from_fits_file(SAMPLE_RHESSI_VISIBILITIES)
+
+    rhessi_map = rhessi_vis.to_map(shape=(65, 65), pixel_size=[4., 4.] * u.arcsec)
+
+    rhessi_vis.vis = np.ones(rhessi_vis.vis.shape)
+    dirty_beam = rhessi_vis.to_image(shape=(65*3, 65*3), pixel_size=[4., 4.] * u.arcsec)
+
+    clean_data, residuals = clean.clean(dirty_map = rhessi_map.data, dirty_beam = dirty_beam,
+                                       gain=0.05, niter=1000, clean_beam_width = 1.0)
+
+    rhessi_map.data[:] = clean_data + residuals
+    rhessi_map.plot()
+
+XRAYVISION also implements a multi-scale version of clean :meth:`~xrayvision.clean.ms_clean`
+
+.. code:: python
+
+    rhessi_vis.vis = np.ones(rhessi_vis.vis.shape)
+    dirty_beam = rhessi_vis.to_image(shape=(65*3, 65*3), pixel_size=[4., 4.] * u.arcsec)
+
+    ms_clean_data, ms_residuals = clean.ms_clean(dirty_map = rhessi_map.data,
+                                                 dirty_beam = dirty_beam,
+                                                 gain=0.05, niter=1000,
+                                                 clean_beam_width = 1.0,
+                                                 scales=[1,2,4])
+
+    rhessi_map.data[:] = clean_data + residuals
+    rhessi_map.plot()
+
+.. plot::
+
+    import numpy as np
+
+    from astropy import units as u
+    from sunpy.map import Map
+    from xrayvision.visibility import RHESSIVisibility
+    from xrayvision import clean, SAMPLE_RHESSI_VISIBILITIES
+
+    rhessi_vis = RHESSIVisibility.from_fits_file(SAMPLE_RHESSI_VISIBILITIES)
+
+    rhessi_map = rhessi_vis.to_map(shape=(65, 65), pixel_size=[4., 4.] * u.arcsec)
+
+    rhessi_vis.vis = np.ones(rhessi_vis.vis.shape)
+    dirty_beam = rhessi_vis.to_image(shape=(65*3, 65*3), pixel_size=[4., 4.] * u.arcsec)
+
+    ms_clean_data, ms_residuals = clean.ms_clean(dirty_map = rhessi_map.data,
+                                                       dirty_beam = dirty_beam,
+                                                       gain=0.05, niter=1000,
+                                                       clean_beam_width = 1.0,
+                                                       scales=[1,2,4])
+
+    rhessi_map.data[:] = ms_clean_data + ms_residuals
+    rhessi_map.plot()

xrayvision/clean.py

@@ -15,6 +15,12 @@
 
 __all__ = ['clean', 'ms_clean']
 
+import logging
+import sys
+
+logging.basicConfig(stream=sys.stdout, level=logging.WARNING)
+logger = logging.getLogger(__name__)
+
 
 def clean(dirty_map, dirty_beam, clean_beam_width=4.0, gain=0.1, thres=0.01, niter=5000):
     r"""
@@ -90,7 +96,7 @@ def clean(dirty_map, dirty_beam, clean_beam_width=4.0, gain=0.1, thres=0.01, nit
         # imax = imax * max_beam
         model[mx, my] += gain * imax
 
-        print(f"Iter: {i}, strength: {imax}, location: {mx, my}")
+        logger.info(f"Iter: {i}, strength: {imax}, location: {mx, my}")
 
         offset = map_center[0] - mx, map_center[1] - my
         shifted_beam_center = int(beam_center[0] + offset[0]), int(beam_center[1] + offset[1])
@@ -104,11 +110,11 @@ def clean(dirty_map, dirty_beam, clean_beam_width=4.0, gain=0.1, thres=0.01, nit
         dirty_map = np.subtract(dirty_map, comp)
 
         if dirty_map.max() <= thres:
-            print("Threshold reached")
+            logger.info("Threshold reached")
             # break
         # el
         if np.abs(dirty_map.min()) > dirty_map.max():
-            print("Largest residual negative")
+            logger.info("Largest residual negative")
             break
 
     else:
@@ -121,7 +127,7 @@ def clean(dirty_map, dirty_beam, clean_beam_width=4.0, gain=0.1, thres=0.01, nit
         model = signal.convolve2d(model, clean_beam, mode='same')
 
         # clean_beam = clean_beam * (1/clean_beam.max())
-        dirty_map = dirty_map / clean_beam.sum()
+        dirty_map = dirty_map * clean_beam.sum()
 
     return model, dirty_map
 
@@ -221,7 +227,7 @@ def ms_clean(dirty_map, dirty_beam, scales=None,
         # Adjust for the max of scaled beam
         strength = strength / max_scaled_dirty_beams[max_scale]
 
-        print(f"Iter: {i}, max scale: {max_scale}, strength: {strength}")
+        logger.info(f"Iter: {i}, max scale: {max_scale}, strength: {strength}")
 
         # Loop gain and scale dependent bias
         strength = strength * scale_biases[max_scale] * gain
@@ -260,23 +266,24 @@ def ms_clean(dirty_map, dirty_beam, scales=None,
 
         # End max(res(a)) or niter
         if scaled_residuals[:, :, max_scale].max() <= thres:
-            print("Threshold reached")
+            logger.info("Threshold reached")
             # break
 
         # Largest scales largest residual is negative
         if np.abs(scaled_residuals[:, :, 0].min()) > scaled_residuals[:, :, 0].max():
-            print("Max scale residual negative")
+            logger.info("Max scale residual negative")
             break
 
     else:
-        print("Max iterations reached")
+        logger.info("Max iterations reached")
 
     # Convolve model with clean beam  B_G * I^M
     if clean_beam_width != 0.0:
         clean_beam = Gaussian2DKernel(stddev=clean_beam_width, x_size=dirty_beam.shape[1],
                                       y_size=dirty_beam.shape[0]).array
 
         model = signal.convolve2d(model, clean_beam, mode='same')  # noqa
+        return model, scaled_residuals.sum(axis=2) * clean_beam.sum()
 
     # Add residuals B_G * I^M + I^R
     return model, scaled_residuals.sum(axis=2)