Fix array_finalize for beams

.travis.yml

@@ -96,8 +96,6 @@ matrix:
         # versions of Python, we can vary Python and Numpy versions at the same
         # time.
 
-        - os: linux
-          env: PYTHON_VERSION=3.5 NUMPY_VERSION=1.12
         - os: linux
           env: PYTHON_VERSION=3.6 NUMPY_VERSION=1.13
         - os: linux

CHANGES.rst

@@ -1,5 +1,8 @@
 0.3 (unreleased)
 ----------------
+ - Set mult/div for convolution/deconvolution in `Beam` and `Beams`.
+   The `==` and `!=` operators also work with `Beams` now.
+   (https://github.com/radio-astro-tools/radio-beam/pull/75)
  - Added common beam operations to `Beams`.
    (https://github.com/radio-astro-tools/radio-beam/pull/67)
  - Fix PA usage for plotting and kernel routines.

radio_beam/beam.py

@@ -348,6 +348,9 @@ def __mul__(self, other):
     def __sub__(self, other):
         return self.deconvolve(other)
 
+    def __truediv__(self, other):
+        return self.deconvolve(other)
+
     def deconvolve(self, other, failure_returns_pointlike=False):
         """
         Deconvolve a beam from another

radio_beam/conftest.py

@@ -2,8 +2,19 @@
 # by importing them here in conftest.py they are discoverable by py.test
 # no matter how it is invoked within the source tree.
 
-from astropy.tests.pytest_plugins import *
+from astropy.version import version as astropy_version
+if astropy_version < '3.0':
+    # With older versions of Astropy, we actually need to import the pytest
+    # plugins themselves in order to make them discoverable by pytest.
+    from astropy.tests.pytest_plugins import *
+else:
+    # As of Astropy 3.0, the pytest plugins provided by Astropy are
+    # automatically made available when Astropy is installed. This means it's
+    # not necessary to import them here, but we still need to import global
+    # variables that are used for configuration.
+    from astropy.tests.plugins.display import PYTEST_HEADER_MODULES, TESTED_VERSIONS
 
+from astropy.tests.helper import enable_deprecations_as_exceptions
 ## Uncomment the following line to treat all DeprecationWarnings as
 ## exceptions
 # enable_deprecations_as_exceptions()

radio_beam/multiple_beams.py

@@ -8,6 +8,7 @@
 
 from .beam import Beam, _to_area, SIGMA_TO_FWHM
 from .commonbeam import commonbeam
+from .utils import InvalidBeamOperationError
 
 
 class Beams(u.Quantity):
@@ -16,7 +17,8 @@ class Beams(u.Quantity):
     """
 
     def __new__(cls, major=None, minor=None, pa=None,
-                areas=None, default_unit=u.arcsec, meta=None):
+                areas=None, default_unit=u.arcsec, meta=None,
+                beams=None):
         """
         Create a new set of Gaussian beams
 
@@ -34,12 +36,20 @@ def __new__(cls, major=None, minor=None, pa=None,
             Gaussian beam.
         default_unit : :class:`~astropy.units.Unit`
             The unit to impose on major, minor if they are specified as floats
+        beams : List of :class:`~radio_beam.Beam` objects
+            List of individual `Beam` objects. The resulting `Beams` object will
+            have major and minor axes in degrees.
         """
 
         # improve to some kwargs magic later
 
         # error checking
 
+        if beams is not None:
+            major = [beam.major.to(u.deg).value for beam in beams] * u.deg
+            minor = [beam.major.to(u.deg).value for beam in beams] * u.deg
+            pa = [beam.pa.to(u.deg).value for beam in beams] * u.deg
+
         # ... given an area make a round beam assuming it is Gaussian
         if areas is not None:
             rad = np.sqrt(areas / (2 * np.pi)) * u.deg
@@ -143,8 +153,10 @@ def __array_finalize__(self, obj):
                 self._set_unit(unit)
 
         if isinstance(obj, Beams):
+            # Multiplication and division should change the area,
+            # but not the PA or major/minor ratio
             self.major = obj.major
-            self.minajor = obj.minor
+            self.minor = obj.minor
             self.pa = obj.pa
             self.meta = obj.meta
 
@@ -325,3 +337,60 @@ def common_beam(self, includemask=None, method='pts', **kwargs):
     def __iter__(self):
         for i in range(len(self)):
             yield self[i]
+
+    def __mul__(self, other):
+        # Other must be a single beam. Assume multiplying is convolving
+        # as set of beams with a given beam
+        if not isinstance(other, Beam):
+            raise InvalidBeamOperationError("Multiplication is defined as a "
+                                            "convolution of the set of beams "
+                                            "with a given beam. Must be "
+                                            "multiplied with a Beam object.")
+
+        return Beams(beams=[beam * other for beam in self])
+
+    def __truediv__(self, other):
+        # Other must be a single beam. Assume dividing is deconvolving
+        # as set of beams with a given beam
+        if not isinstance(other, Beam):
+            raise InvalidBeamOperationError("Division is defined as a "
+                                            "deconvolution of the set of beams"
+                                            " with a given beam. Must be "
+                                            "divided by a Beam object.")
+
+        return Beams(beams=[beam - other for beam in self])
+
+    def __add__(self, other):
+        raise InvalidBeamOperationError("Addition of a set of Beams "
+                                        "is not defined.")
+
+    def __sub__(self, other):
+        raise InvalidBeamOperationError("Addition of a set of Beams "
+                                        "is not defined.")
+
+    def __eq__(self, other):
+        # other should be a single beam, or a another Beams object
+        if isinstance(other, Beam):
+            return np.array([beam == other for beam in self])
+        elif isinstance(other, Beams):
+            # These should have the same size.
+            if not self.size == other.size:
+                raise InvalidBeamOperationError("Beams objects must have the "
+                                                "same shape to test "
+                                                "equality.")
+
+            return np.all([beam == other_beam for beam, other_beam in
+                           zip(self, other)])
+        else:
+            raise InvalidBeamOperationError("Must test equality with a Beam"
+                                            " or Beams object.")
+
+    def __ne__(self, other):
+        eq_out = self.__eq__(other)
+
+        # If other is a Beam, will get array back
+        if isinstance(eq_out, np.ndarray):
+            return ~eq_out
+        # If other is a Beams, will get boolean back
+        else:
+            return not eq_out

radio_beam/tests/test_beam.py

@@ -247,6 +247,10 @@ def test_conv_deconv():
     assert beam2 == beam3 - beam1
     assert beam1 == beam3 - beam2
 
+    # Dividing should give the same thing
+    assert beam2 == beam3 / beam1
+    assert beam1 == beam3 / beam2
+
     assert beam3 == beam1 * beam2
 
 

radio_beam/tests/test_beams.py

@@ -9,6 +9,7 @@
 from ..multiple_beams import Beams
 from ..beam import Beam
 from ..commonbeam import common_2beams, common_manybeams_mve
+from ..utils import InvalidBeamOperationError
 
 from .test_beam import data_path
 
@@ -62,6 +63,132 @@ def test_beams_from_fits_bintable():
     assert (beams.pa.value == bintable.data['BPA']).all()
 
 
+def test_beams_from_list_of_beam():
+
+    beams, majors = symm_beams_for_tests()[:2]
+
+    new_beams = Beams(beams=[beam for beam in beams])
+
+    assert beams == new_beams
+
+
+def test_beams_equality_beams():
+
+    beams, majors = symm_beams_for_tests()[:2]
+
+    assert beams == beams
+
+    assert not beams != beams
+
+    abeams, amajors = asymm_beams_for_tests()[:2]
+
+    assert not (beams == abeams)
+
+    assert beams != abeams
+
+
+def test_beams_equality_beam():
+
+    # Test whether all are equal to a single beam
+    beams = Beams([1.] * 5 * u.arcsec)
+
+    beam = Beam(1 * u.arcsec)
+
+    assert np.all(beams == beam)
+
+    assert not np.any(beams != beam)
+
+
+@pytest.mark.xfail(raises=InvalidBeamOperationError, strict=True)
+def test_beams_equality_fail():
+
+    # Test whether all are equal to a single beam
+    beams = Beams([1.] * 5 * u.arcsec)
+
+    beams == 2
+
+
+@pytest.mark.xfail(raises=InvalidBeamOperationError, strict=True)
+def test_beams_notequality_fail():
+
+    # Test whether all are equal to a single beam
+    beams = Beams([1.] * 5 * u.arcsec)
+
+    beams != 2
+
+
+@pytest.mark.xfail(raises=InvalidBeamOperationError, strict=True)
+def test_beams_equality_fail_shape():
+
+    # Test whether all are equal to a single beam
+    beams = Beams([1.] * 5 * u.arcsec)
+
+    assert np.all(beams == beams[1:])
+
+@pytest.mark.xfail(raises=InvalidBeamOperationError, strict=True)
+def test_beams_add_fail():
+
+    # Test whether all are equal to a single beam
+    beams = Beams([1.] * 5 * u.arcsec)
+
+    beams + 2
+
+
+@pytest.mark.xfail(raises=InvalidBeamOperationError, strict=True)
+def test_beams_sub_fail():
+
+    # Test whether all are equal to a single beam
+    beams = Beams([1.] * 5 * u.arcsec)
+
+    beams - 2
+
+
+@pytest.mark.xfail(raises=InvalidBeamOperationError, strict=True)
+def test_beams_mult_fail():
+
+    # Test whether all are equal to a single beam
+    beams = Beams([1.] * 5 * u.arcsec)
+
+    beams * 2
+
+
+@pytest.mark.xfail(raises=InvalidBeamOperationError, strict=True)
+def test_beams_div_fail():
+
+    # Test whether all are equal to a single beam
+    beams = Beams([1.] * 5 * u.arcsec)
+
+    beams / 2
+
+
+def test_beams_mult_convolution():
+
+    beams, majors = asymm_beams_for_tests()[:2]
+
+    beam = Beam(1 * u.arcsec)
+
+    conv_beams = beams * beam
+
+    individ_conv_beams = [beam_i.convolve(beam) for beam_i in beams]
+    new_beams = Beams(beams=individ_conv_beams)
+
+    assert conv_beams == new_beams
+
+
+def test_beams_div_deconvolution():
+
+    beams, majors = asymm_beams_for_tests()[:2]
+
+    beam = Beam(0.25 * u.arcsec)
+
+    deconv_beams = beams / beam
+
+    individ_deconv_beams = [beam_i.deconvolve(beam) for beam_i in beams]
+    new_beams = Beams(beams=individ_deconv_beams)
+
+    assert deconv_beams == new_beams
+
+
 def test_indexing():
 
     beams, majors = symm_beams_for_tests()[:2]
@@ -91,6 +218,12 @@ def test_indexing():
     assert np.all(beams[mask].major.value == majors[mask].value)
 
 
+# def test_beams_mult_convolution():
+
+#     beams, majors = symm_beams_for_tests()[:2]
+
+
+
 def test_average_beams():
 
     beams, majors = symm_beams_for_tests()[:2]

radio_beam/utils.py

@@ -8,6 +8,10 @@ class BeamError(Exception):
     pass
 
 
+class InvalidBeamOperationError(Exception):
+    pass
+
+
 def deconvolve(beam, other, failure_returns_pointlike=False):
     """
     Deconvolve a beam from another

setup.cfg

@@ -7,7 +7,7 @@ all_files = 1
 upload-dir = docs/_build/html
 show-response = 1
 
-[pytest]
+[tool:pytest]
 minversion = 2.2
 norecursedirs = build docs/_build
 doctest_plus = enabled