Replace asarray_of_arrays() with XyzComponents

Author: mgeier

sfs/mono/source.py

@@ -53,7 +53,7 @@ def point(omega, x0, n0, grid, c=None):
     """
     k = util.wavenumber(omega, c)
     x0 = util.asarray_1d(x0)
-    grid = util.asarray_of_arrays(grid)
+    grid = util.XyzComponents(grid)
 
     r = np.linalg.norm(grid - x0)
     return 1 / (4*np.pi) * np.exp(-1j * k * r) / r
@@ -65,7 +65,7 @@ def point_velocity(omega, x0, n0, grid, c=None):
     """
     k = util.wavenumber(omega, c)
     x0 = util.asarray_1d(x0)
-    grid = util.asarray_of_arrays(grid)
+    grid = util.XyzComponents(grid)
     offset = grid - x0
     r = np.linalg.norm(offset)
     v = point(omega, x0, n0, grid, c=c)
@@ -108,7 +108,7 @@ def point_modal(omega, x0, n0, grid, L, N=None, deltan=0, c=None):
     """
     k = util.wavenumber(omega, c)
     x0 = util.asarray_1d(x0)
-    x, y, z = util.asarray_of_arrays(grid)
+    x, y, z = util.XyzComponents(grid)
 
     if N is None:
         # determine maximum modal order per dimension
@@ -178,7 +178,7 @@ def point_modal_velocity(omega, x0, n0, grid, L, N=None, deltan=0, c=None):
     """
     k = util.wavenumber(omega, c)
     x0 = util.asarray_1d(x0)
-    x, y, z = util.asarray_of_arrays(grid)
+    x, y, z = util.XyzComponents(grid)
 
     if N is None:
         # determine maximum modal order per dimension
@@ -251,7 +251,7 @@ def line(omega, x0, n0, grid, c=None):
     k = util.wavenumber(omega, c)
     x0 = util.asarray_1d(x0)
     x0 = x0[:2]  # ignore z-component
-    grid = util.asarray_of_arrays(grid)
+    grid = util.XyzComponents(grid)
 
     r = np.linalg.norm(grid[:2] - x0)
     p = -1j/4 * special.hankel2(0, k * r)
@@ -265,7 +265,7 @@ def line_velocity(omega, x0, n0, grid, c=None):
     k = util.wavenumber(omega, c)
     x0 = util.asarray_1d(x0)
     x0 = x0[:2]  # ignore z-component
-    grid = util.asarray_of_arrays(grid)
+    grid = util.XyzComponents(grid)
 
     offset = grid[:2] - x0
     r = np.linalg.norm(offset)
@@ -296,7 +296,7 @@ def line_dipole(omega, x0, n0, grid, c=None):
     x0 = util.asarray_1d(x0)
     x0 = x0[:2]  # ignore z-component
     n0 = n0[:2]
-    grid = util.asarray_of_arrays(grid)
+    grid = util.XyzComponents(grid)
     dx = grid[:2] - x0
     r = np.linalg.norm(dx)
     p = 1j*k/4 * special.hankel2(1, k * r) * np.inner(dx, n0) / r
@@ -325,7 +325,7 @@ def plane(omega, x0, n0, grid, c=None):
     k = util.wavenumber(omega, c)
     x0 = util.asarray_1d(x0)
     n0 = util.asarray_1d(n0)
-    grid = util.asarray_of_arrays(grid)
+    grid = util.XyzComponents(grid)
     return np.exp(-1j * k * np.inner(grid - x0, n0))
 
 

sfs/plot.py

@@ -145,7 +145,7 @@ def loudspeaker_2d(x0, n0, a0=0.5, size=0.08, show_numbers=False, grid=None,
 
 def _visible_secondarysources_2d(x0, n0, grid):
     """Determine secondary sources which lie within `grid`."""
-    grid = util.asarray_of_arrays(grid)
+    grid = util.XyzComponents(grid)
     x, y = grid[:2]
     idx = np.where((x0[:, 0] > x.min()) & (x0[:, 0] < x.max()) &
                    (x0[:, 1] > y.min()) & (x0[:, 1] < x.max()))
@@ -242,7 +242,7 @@ def soundfield(p, grid, xnorm=None, cmap='coolwarm_clip', colorbar=True,
 
     """
     p = np.asarray(p)
-    grid = util.asarray_of_arrays(grid)
+    grid = util.XyzComponents(grid)
 
     # normalize sound field wrt xnorm
     if xnorm is not None:

sfs/util.py

@@ -98,20 +98,6 @@ def asarray_of_rows(a, **kwargs):
     return result
 
 
-def asarray_of_arrays(a, **kwargs):
-    """Convert the input to an array consisting of arrays.
-
-    A one-dimensional :class:`numpy.ndarray` with `dtype=object` is
-    returned, containing the elements of `a` as :class:`numpy.ndarray`
-    (whose `dtype` and other options can be specified with `**kwargs`).
-
-    """
-    result = np.empty(len(a), dtype=object)
-    for i, element in enumerate(a):
-        result[i] = np.asarray(element, **kwargs)
-    return result
-
-
 def strict_arange(start, stop, step=1, endpoint=False, dtype=None, **kwargs):
     """Like :func:`numpy.arange`, but compensating numeric errors.
 
@@ -202,7 +188,7 @@ def normalize(p, grid, xnorm):
 
 def probe(p, grid, x):
     """Determine the value at position `x` in the sound field `p`."""
-    grid = asarray_of_arrays(grid)
+    grid = XyzComponents(grid)
     x = asarray_1d(x)
     r = np.linalg.norm(grid - x)
     idx = np.unravel_index(r.argmin(), r.shape)
@@ -228,7 +214,7 @@ def displacement(v, omega):
         d(x, t) = \int_0^t v(x, t) dt
 
     """
-    v = asarray_of_arrays(v)
+    v = XyzComponents(v)
     return v / (1j * omega)
 
 
@@ -246,3 +232,91 @@ def db(x, power=False):
     """
     with np.errstate(divide='ignore'):
         return 10 if power else 20 * np.log10(np.abs(x))
+
+
+class XyzComponents(np.ndarray):
+    """See __init__()."""
+
+    def __init__(self, arg, **kwargs):
+        """Triple (or pair) of arrays: x, y, and optionally z.
+
+        Instances of this class can be used to store coordinate grids
+        (either regular grids like in :func:`sfs.util.xyz_grid` or
+        arbitrary point clouds) or vector fields (e.g. particle
+        velocity).
+
+        This class is a subclass of :class:`numpy.ndarray`.  It is
+        one-dimensional (like a plain :class:`list`) and has a length of
+        3 (or 2, if no z-component is available).  It uses
+        `dtype=object` in order to be able to store other
+        `numpy.ndarray`\s of arbitrary shapes but also scalars, if
+        needed.  Because it is a NumPy array subclass, it can be used in
+        operations with scalars and "normal" NumPy arrays, as long as
+        they have a compatible shape.  Like any NumPy array, instances
+        of this class are iterable and can be used, e.g., in for-loops
+        and tuple unpacking.  If slicing or broadcasting leads to an
+        incompatible shape, a plain `numpy.ndarray` with `dtype=object`
+        is returned.
+
+        Parameters
+        ----------
+        arg : triple or pair of array_like
+            The values to be used as X, Y and Z arrays.  Z is optional.
+        **kwargs
+            All further arguments are forwarded to
+            :func:`numpy.asarray`.
+
+        """
+        # This method does nothing, it's only here for the documentation!
+
+    def __new__(cls, arg, **kwargs):
+        # object arrays cannot be created and populated in a single step:
+        obj = np.ndarray.__new__(cls, len(arg), dtype=object)
+        for i, component in enumerate(arg):
+            obj[i] = np.asarray(component, **kwargs)
+        return obj
+
+    def __array_finalize__(self, obj):
+        if self.ndim == 0:
+            pass  # this is allowed, e.g. for np.inner()
+        elif self.ndim > 1 or len(self) not in (2, 3):
+            raise ValueError("XyzComponents can only have 2 or 3 components")
+
+    def __array_prepare__(self, obj, context=None):
+        if obj.ndim == 1 and len(obj) in (2, 3):
+            return obj.view(XyzComponents)
+        return obj
+
+    def __array_wrap__(self, obj, context=None):
+        if obj.ndim != 1 or len(obj) not in (2, 3):
+            return obj.view(np.ndarray)
+        return obj
+
+    def __getitem__(self, index):
+        if isinstance(index, slice):
+            start, stop, step = index.indices(len(self))
+            if start == 0 and stop in (2, 3) and step == 1:
+                return np.ndarray.__getitem__(self, index)
+        # Slices other than xy and xyz are "downgraded" to ndarray
+        return np.ndarray.__getitem__(self.view(np.ndarray), index)
+
+    def __repr__(self):
+        return 'XyzComponents(\n' + ',\n'.join(
+            '    {0}={1}'.format(name, repr(data).replace('\n', '\n      '))
+            for name, data in zip('xyz', self)) + ')'
+
+    def make_property(index, doc):
+
+        def getter(self):
+            return self[index]
+
+        def setter(self, value):
+            self[index] = value
+
+        return property(getter, setter, doc=doc)
+
+    x = make_property(0, doc='x-component.')
+    y = make_property(1, doc='y-component.')
+    z = make_property(2, doc='z-component (optional).')
+
+    del make_property