From 2590c1ccb37edc4b39751f740f9343f6d7b09afc Mon Sep 17 00:00:00 2001
From: Alix Damman <ald@plan.be>
Date: Mon, 19 Aug 2019 09:16:12 +0200
Subject: [PATCH 1/3] moved imports into functions instead of on top of the
 module tests/common.py to avoid import loop problems

---
 larray/tests/common.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/larray/tests/common.py b/larray/tests/common.py
index e3126d5fe..084234d35 100644
--- a/larray/tests/common.py
+++ b/larray/tests/common.py
@@ -12,8 +12,6 @@
 except ImportError:
     xw = None
 
-from larray import LArray, isnan, aslarray, Metadata
-
 
 TESTDATADIR = os.path.dirname(__file__)
 
@@ -38,6 +36,7 @@ def inputpath(relpath):
 
 def assert_equal_factory(test_func):
     def assert_equal(a, b):
+        from larray.core.array import LArray
         if isinstance(a, LArray) and isinstance(b, LArray) and a.axes != b.axes:
             raise AssertionError("axes differ:\n%s\n\nvs\n\n%s" % (a.axes.info, b.axes.info))
         if not isinstance(a, (np.ndarray, LArray)):
@@ -57,6 +56,7 @@ def assert_equal(a, b):
 
 def assert_larray_equal_factory(test_func, convert=True, check_axes=False):
     def assert_equal(a, b):
+        from larray.core.array import aslarray
         if convert:
             a = aslarray(a)
             b = aslarray(b)
@@ -88,6 +88,7 @@ def equal(a, b):
 
 
 def nan_equal(a, b):
+    from larray.core.ufuncs import isnan
     return (a == b) | (isnan(a) & isnan(b))
 
 
@@ -121,6 +122,7 @@ def tmp_path(tmpdir, fname):
 
 @pytest.fixture
 def meta():
+    from larray.core.metadata import Metadata
     title = 'test array'
     description = 'Array used for testing'
     author = 'John Cleese'

From f72af2045d0853f9bb55b3bb16e53061ea4cb1e7 Mon Sep 17 00:00:00 2001
From: Alix Damman <ald@plan.be>
Date: Fri, 12 Apr 2019 15:59:24 +0200
Subject: [PATCH 2/3] - fix #234 : implemented support for np.__array_ufunc__

- merged npufuncs.py module into ufuncs.py

- updated setup.cfg config file to skip all functions derived from Numpy and for which copied doctests fail
---
 doc/source/changes/version_0_32.rst.inc |   2 +-
 larray/__init__.py                      |  18 +-
 larray/core/array.py                    |   7 +-
 larray/core/npufuncs.py                 | 135 --------
 larray/core/ufuncs.py                   | 397 ++++++++++++++++++------
 larray/tests/test_array.py              |  73 -----
 larray/tests/test_numpy_funcs.py        |  87 ++++++
 setup.cfg                               |  35 ++-
 8 files changed, 438 insertions(+), 316 deletions(-)
 delete mode 100644 larray/core/npufuncs.py
 create mode 100644 larray/tests/test_numpy_funcs.py

diff --git a/doc/source/changes/version_0_32.rst.inc b/doc/source/changes/version_0_32.rst.inc
index bbb2f6750..4898269dc 100644
--- a/doc/source/changes/version_0_32.rst.inc
+++ b/doc/source/changes/version_0_32.rst.inc
@@ -49,7 +49,7 @@ New features
 Miscellaneous improvements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-* improved something.
+* made it possible to pass LArray objects to Numpy ufuncs directly (e.g. np.sqrt(ndtest(3)), closes :issue:`234`).
 
 
 Fixes
diff --git a/larray/__init__.py b/larray/__init__.py
index ab65cf548..38446e900 100644
--- a/larray/__init__.py
+++ b/larray/__init__.py
@@ -12,15 +12,15 @@
 from larray.core.session import Session, local_arrays, global_arrays, arrays
 from larray.core.constants import nan, inf, pi, e, euler_gamma
 from larray.core.metadata import Metadata
-from larray.core.ufuncs import wrap_elementwise_array_func, maximum, minimum, where
-from larray.core.npufuncs import (sin, cos, tan, arcsin, arccos, arctan, hypot, arctan2, degrees,
-                                  radians, unwrap, sinh, cosh, tanh, arcsinh, arccosh, arctanh,
-                                  angle, real, imag, conj,
-                                  round, around, round_, rint, fix, floor, ceil, trunc,
-                                  exp, expm1, exp2, log, log10, log2, log1p, logaddexp, logaddexp2,
-                                  i0, sinc, signbit, copysign, frexp, ldexp,
-                                  convolve, clip, sqrt, absolute, fabs, sign, fmax, fmin, nan_to_num,
-                                  real_if_close, interp, isnan, isinf, inverse)
+from larray.core.ufuncs import (wrap_elementwise_array_func, maximum, minimum, where,
+                                sin, cos, tan, arcsin, arccos, arctan, hypot, arctan2, degrees,
+                                radians, unwrap, sinh, cosh, tanh, arcsinh, arccosh, arctanh,
+                                angle, real, imag, conj,
+                                round, around, round_, rint, fix, floor, ceil, trunc,
+                                exp, expm1, exp2, log, log10, log2, log1p, logaddexp, logaddexp2,
+                                i0, sinc, signbit, copysign, frexp, ldexp,
+                                convolve, clip, sqrt, absolute, fabs, sign, fmax, fmin, nan_to_num,
+                                real_if_close, interp, isnan, isinf, inverse)
 
 from larray.inout.misc import from_lists, from_string
 from larray.inout.pandas import from_frame, from_series
diff --git a/larray/core/array.py b/larray/core/array.py
index 17612fa91..e880be02b 100644
--- a/larray/core/array.py
+++ b/larray/core/array.py
@@ -57,6 +57,7 @@
 from larray.core.constants import nan
 from larray.core.metadata import Metadata
 from larray.core.expr import ExprNode
+from larray.core.ufuncs import SupportsNumpyUfuncs
 from larray.core.group import (Group, IGroup, LGroup, remove_nested_groups, _to_key, _to_keys,
                                _translate_sheet_name, _translate_group_key_hdf)
 from larray.core.axis import Axis, AxisReference, AxisCollection, X, _make_axis
@@ -700,7 +701,7 @@ def _handle_meta(meta, title):
     return Metadata(meta)
 
 
-class LArray(ABCLArray):
+class LArray(ABCLArray, SupportsNumpyUfuncs):
     r"""
     A LArray object represents a multidimensional, homogeneous array of fixed-size items with labeled axes.
 
@@ -5862,7 +5863,7 @@ def eq(self, other, rtol=0, atol=0, nans_equal=False):
             if not nans_equal:
                 return self == other
             else:
-                from larray.core.npufuncs import isnan
+                from larray.core.ufuncs import isnan
 
                 def general_isnan(a):
                     if np.issubclass_(a.dtype.type, np.inexact):
@@ -6625,7 +6626,7 @@ def clip(self, minval=None, maxval=None, out=None):
          a1   0   1   2
          a2   2   2   2
         """
-        from larray.core.npufuncs import clip
+        from larray.core.ufuncs import clip
         return clip(self, minval, maxval, out)
 
     @deprecate_kwarg('transpose', 'wide')
diff --git a/larray/core/npufuncs.py b/larray/core/npufuncs.py
deleted file mode 100644
index 4203f8d85..000000000
--- a/larray/core/npufuncs.py
+++ /dev/null
@@ -1,135 +0,0 @@
-# numpy ufuncs -- this module is excluded from pytest
-# http://docs.scipy.org/doc/numpy/reference/routines.math.html
-
-import numpy as np
-
-from larray.core.ufuncs import broadcastify
-
-
-# Trigonometric functions
-
-sin = broadcastify(np.sin)
-cos = broadcastify(np.cos)
-tan = broadcastify(np.tan)
-arcsin = broadcastify(np.arcsin)
-arccos = broadcastify(np.arccos)
-arctan = broadcastify(np.arctan)
-hypot = broadcastify(np.hypot)
-arctan2 = broadcastify(np.arctan2)
-degrees = broadcastify(np.degrees)
-radians = broadcastify(np.radians)
-unwrap = broadcastify(np.unwrap)
-# deg2rad = broadcastify(np.deg2rad)
-# rad2deg = broadcastify(np.rad2deg)
-
-# Hyperbolic functions
-
-sinh = broadcastify(np.sinh)
-cosh = broadcastify(np.cosh)
-tanh = broadcastify(np.tanh)
-arcsinh = broadcastify(np.arcsinh)
-arccosh = broadcastify(np.arccosh)
-arctanh = broadcastify(np.arctanh)
-
-# Rounding
-
-# TODO: add examples for round, floor, ceil and trunc
-# all 3 are equivalent, I am unsure I should support around and round_
-round = broadcastify(np.round)
-around = broadcastify(np.around)
-round_ = broadcastify(np.round_)
-rint = broadcastify(np.rint)
-fix = broadcastify(np.fix)
-floor = broadcastify(np.floor)
-ceil = broadcastify(np.ceil)
-trunc = broadcastify(np.trunc)
-
-# Sums, products, differences
-
-# prod = broadcastify(np.prod)
-# sum = broadcastify(np.sum)
-# nansum = broadcastify(np.nansum)
-# cumprod = broadcastify(np.cumprod)
-# cumsum = broadcastify(np.cumsum)
-
-# cannot use a simple wrapped ufunc because those ufuncs do not preserve
-# shape or dimensions so labels are wrong
-# diff = broadcastify(np.diff)
-# ediff1d = broadcastify(np.ediff1d)
-# gradient = broadcastify(np.gradient)
-# cross = broadcastify(np.cross)
-# trapz = broadcastify(np.trapz)
-
-# Exponents and logarithms
-
-# TODO: add examples for exp and log
-exp = broadcastify(np.exp)
-expm1 = broadcastify(np.expm1)
-exp2 = broadcastify(np.exp2)
-log = broadcastify(np.log)
-log10 = broadcastify(np.log10)
-log2 = broadcastify(np.log2)
-log1p = broadcastify(np.log1p)
-logaddexp = broadcastify(np.logaddexp)
-logaddexp2 = broadcastify(np.logaddexp2)
-
-# Other special functions
-
-i0 = broadcastify(np.i0)
-sinc = broadcastify(np.sinc)
-
-# Floating point routines
-
-signbit = broadcastify(np.signbit)
-copysign = broadcastify(np.copysign)
-frexp = broadcastify(np.frexp)
-ldexp = broadcastify(np.ldexp)
-
-# Arithmetic operations
-
-# add = broadcastify(np.add)
-# reciprocal = broadcastify(np.reciprocal)
-# negative = broadcastify(np.negative)
-# multiply = broadcastify(np.multiply)
-# divide = broadcastify(np.divide)
-# power = broadcastify(np.power)
-# subtract = broadcastify(np.subtract)
-# true_divide = broadcastify(np.true_divide)
-# floor_divide = broadcastify(np.floor_divide)
-# fmod = broadcastify(np.fmod)
-# mod = broadcastify(np.mod)
-# modf = broadcastify(np.modf)
-# remainder = broadcastify(np.remainder)
-
-# Handling complex numbers
-
-angle = broadcastify(np.angle)
-real = broadcastify(np.real)
-imag = broadcastify(np.imag)
-conj = broadcastify(np.conj)
-
-# Miscellaneous
-
-convolve = broadcastify(np.convolve)
-clip = broadcastify(np.clip)
-# square = broadcastify(np.square)
-absolute = broadcastify(np.absolute)
-fabs = broadcastify(np.fabs)
-sign = broadcastify(np.sign)
-fmax = broadcastify(np.fmax)
-fmin = broadcastify(np.fmin)
-nan_to_num = broadcastify(np.nan_to_num)
-real_if_close = broadcastify(np.real_if_close)
-
-# TODO: add examples for functions below
-sqrt = broadcastify(np.sqrt)
-isnan = broadcastify(np.isnan)
-isinf = broadcastify(np.isinf)
-inverse = broadcastify(np.linalg.inv)
-
-# XXX: create a new LArray method instead ?
-# TODO: should appear in the API doc if it actually works with LArrays,
-#       which I have never tested (and I doubt is the case).
-#       Might be worth having specific documentation if it works well.
-#       My guess is that we should rather make a new LArray method for that one.
-interp = broadcastify(np.interp)
diff --git a/larray/core/ufuncs.py b/larray/core/ufuncs.py
index 6a3b8c8a7..bce68ea41 100644
--- a/larray/core/ufuncs.py
+++ b/larray/core/ufuncs.py
@@ -1,9 +1,10 @@
-# numpy ufuncs
+# numpy ufuncs -- this module is excluded from pytest
 # http://docs.scipy.org/doc/numpy/reference/routines.math.html
+# https://docs.scipy.org/doc/numpy/reference/ufuncs.html#available-ufuncs
 
-import numpy as np
+import numbers
 
-from larray.core.array import LArray, make_args_broadcastable
+import numpy as np
 
 
 def wrap_elementwise_array_func(func):
@@ -50,6 +51,7 @@ def wrap_elementwise_array_func(func):
            F  -0.37  -0.61  -0.87
     """
     def wrapper(*args, **kwargs):
+        from larray.core.array import LArray, make_args_broadcastable
         raw_bcast_args, raw_bcast_kwargs, res_axes = make_args_broadcastable(args, kwargs)
 
         # We pass only raw numpy arrays to the ufuncs even though numpy is normally meant to handle those cases itself
@@ -78,95 +80,182 @@ def wrapper(*args, **kwargs):
     return wrapper
 
 
-# TODO: rename to wrap_numpy_func
-def broadcastify(func):
-    wrapper = wrap_elementwise_array_func(func)
-    # update documentation by inserting a warning message after the short description of the numpy function
-    # (otherwise the description of ufuncs given in the corresponding API 'autosummary' tables will always
-    #  start with 'larray specific variant of ...' without giving a meaningful description of what does the ufunc)
-    if func.__doc__.startswith('\n'):
-        # docstring starts with short description
-        end_signature = 1
-        end_short_desc = func.__doc__.find('\n\n')
-    else:
-        # docstring starts with signature
-        end_signature = func.__doc__.find('\n\n') + 2
-        end_short_desc = func.__doc__.find('\n\n', end_signature)
-    short_desc = func.__doc__[:end_short_desc]
-    numpy_doc = func.__doc__[end_short_desc:]
-    ident = ' ' * (len(short_desc[end_signature:]) - len(short_desc[end_signature:].lstrip()))
-    wrapper.__doc__ = '{short_desc}' \
-                      '\n\n{ident}larray specific variant of ``numpy.{fname}``.' \
-                      '\n\n{ident}Documentation from numpy:' \
-                      '{numpy_doc}'.format(short_desc=short_desc, ident=ident, fname=func.__name__, numpy_doc=numpy_doc)
-    # set __qualname__ explicitly (all these functions are supposed to be top-level function in the ufuncs module)
-    wrapper.__qualname__ = func.__name__
-    # we should not copy __module__
-    return wrapper
-
-
-where = broadcastify(np.where)
-where.__doc__ = r"""
-where(condition, x, y)
-
-    Return elements, either from `x` or `y`, depending on `condition`.
-
-    Parameters
-    ----------
-    condition : boolean LArray
-        When True, yield `x`, otherwise yield `y`.
-    x, y : LArray
-        Values from which to choose.
-
-    Returns
-    -------
-    out : LArray
-        If both `x` and `y` are specified, the output array contains
-        elements of `x` where `condition` is True, and elements from
-        `y` elsewhere.
-
-    Examples
-    --------
-    >>> from larray import LArray
-    >>> arr = LArray([[10, 7, 5, 9],
-    ...               [5, 8, 3, 7],
-    ...               [6, 2, 0, 9],
-    ...               [9, 10, 5, 6]], "a=a0..a3;b=b0..b3")
-    >>> arr
-    a\b  b0  b1  b2  b3
-     a0  10   7   5   9
-     a1   5   8   3   7
-     a2   6   2   0   9
-     a3   9  10   5   6
-
-    Simple use
-
-    >>> where(arr <= 5, 0, arr)
-    a\b  b0  b1  b2  b3
-     a0  10   7   0   9
-     a1   0   8   0   7
-     a2   6   0   0   9
-     a3   9  10   0   6
-
-    With broadcasting
+class SupportsNumpyUfuncs(object):
+    """
+    Base class for larray types that support ufuncs.
+    Used by LArray.
 
-    >>> mean_by_col = arr.mean('a')
-    >>> mean_by_col
-    b   b0    b1    b2    b3
-       7.5  6.75  3.25  7.75
-    >>> # for each column, set values below the mean value to the mean value
-    >>> where(arr < mean_by_col, mean_by_col, arr)
-    a\b    b0    b1    b2    b3
-     a0  10.0   7.0   5.0   9.0
-     a1   7.5   8.0  3.25  7.75
-     a2   7.5  6.75  3.25   9.0
-     a3   9.0  10.0   5.0  7.75
-"""
+    Notes
+    -----
+    See https://docs.scipy.org/doc/numpy/reference/arrays.classes.html#numpy.class.__array_ufunc__
+    for more details about __array_ufunc__ and
+    https://docs.scipy.org/doc/numpy/reference/generated/numpy.lib.mixins.NDArrayOperatorsMixin.html
+    for an example.
+    """
+    _HANDLED_TYPES = (np.ndarray, np.generic, numbers.Number, bytes, str)
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        """
+        Override the behavior of NumPy’s ufuncs.
+
+        Parameters
+        ----------
+        ufunc: callable
+            Ufunc object that was called.
+        method: str
+            String indicating which Ufunc method was called
+            (one of "__call__", "reduce", "reduceat", "accumulate", "outer", "inner").
+        inputs: tuple
+            Input arguments to the ufunc.
+        kwargs: dict
+            Dictionary containing the optional input arguments of the ufunc.
+            If given, any out arguments, both positional and keyword, are passed as a tuple in kwargs.
+        """
+        out = kwargs.get('out', ())
+        for x in inputs + out:
+            if not isinstance(x, self._HANDLED_TYPES + (SupportsNumpyUfuncs,)):
+                return NotImplemented
+
+        if out:
+            if len(out) > 1:
+                raise TypeError("Passing an iterable for the argument 'out' is not supported")
+            kwargs['out'] = out[0]
+
+        if ufunc.signature is not None:
+            # In regular ufuncs, the elementary function is limited to element-by-element operations,
+            # whereas the generalized version (gufuncs) supports "sub-array" by "sub-array" operations.
+            raise NotImplementedError('{} not supported: larray objects do not directly implement '
+                                      'generalized ufuncs.'.format(ufunc))
+
+        if method != '__call__':
+            raise NotImplemented
+
+        wrapped_ufunc = wrap_elementwise_array_func(ufunc)
+        return wrapped_ufunc(*inputs, **kwargs)
+
+
+def broadcastify(npfunc, copy_numpy_doc=True):
+    wrapper = wrap_elementwise_array_func(npfunc)
+    if copy_numpy_doc:
+        # copy meaningful attributes (numpy ufuncs do not have __annotations__ nor __qualname__)
+        wrapper.__name__ = npfunc.__name__
+        # update documentation by inserting a warning message after the short description of the numpy function
+        # (otherwise the description of npfuncs given in the corresponding API 'autosummary' tables will always
+        #  start with 'larray specific variant of ...' without giving a meaningful description of what does the npfunc)
+        if npfunc.__doc__.startswith('\n'):
+            # docstring starts with short description
+            end_signature = 1
+            end_short_desc = npfunc.__doc__.find('\n\n')
+        else:
+            # docstring starts with signature
+            end_signature = npfunc.__doc__.find('\n\n') + 2
+            end_short_desc = npfunc.__doc__.find('\n\n', end_signature)
+        short_desc = npfunc.__doc__[:end_short_desc]
+        numpy_doc = npfunc.__doc__[end_short_desc:]
+        ident = ' ' * (len(short_desc[end_signature:]) - len(short_desc[end_signature:].lstrip()))
+        wrapper.__doc__ = '{short_desc}' \
+                          '\n\n{ident}larray specific variant of ``numpy.{fname}``.' \
+                          '\n\n{ident}Documentation from numpy:' \
+                          '{numpy_doc}'.format(short_desc=short_desc, ident=ident,
+                                               fname=npfunc.__name__, numpy_doc=numpy_doc)
+        # set __qualname__ explicitly (all these functions are supposed to be top-level function in the npfuncs module)
+        wrapper.__qualname__ = npfunc.__name__
+        # we should not copy __module__
+    return wrapper
 
-maximum = broadcastify(np.maximum)
-maximum.__doc__ = r"""
-maximum(x1, x2, out=None, dtype=None)
 
+# list of available Numpy ufuncs
+# https://docs.scipy.org/doc/numpy/reference/ufuncs.html#available-ufuncs
+
+# Trigonometric functions
+
+sin = broadcastify(np.sin)
+cos = broadcastify(np.cos)
+tan = broadcastify(np.tan)
+arcsin = broadcastify(np.arcsin)
+arccos = broadcastify(np.arccos)
+arctan = broadcastify(np.arctan)
+arctan2 = broadcastify(np.arctan2)
+hypot = broadcastify(np.hypot)
+degrees = broadcastify(np.degrees)
+radians = broadcastify(np.radians)
+unwrap = broadcastify(np.unwrap)
+# deg2rad = broadcastify(np.deg2rad)
+# rad2deg = broadcastify(np.rad2deg)
+
+# Hyperbolic functions
+
+sinh = broadcastify(np.sinh)
+cosh = broadcastify(np.cosh)
+tanh = broadcastify(np.tanh)
+arcsinh = broadcastify(np.arcsinh)
+arccosh = broadcastify(np.arccosh)
+arctanh = broadcastify(np.arctanh)
+
+# Rounding
+
+rint = broadcastify(np.rint)
+# TODO: fix fails because of its 'out' argument
+fix = broadcastify(np.fix)
+# TODO: add examples for round, floor, ceil and trunc
+floor = broadcastify(np.floor)
+ceil = broadcastify(np.ceil)
+trunc = broadcastify(np.trunc)
+
+# Exponents and logarithms
+
+# TODO: add examples for exp and log
+exp = broadcastify(np.exp)
+exp2 = broadcastify(np.exp2)
+expm1 = broadcastify(np.expm1)
+log = broadcastify(np.log)
+log10 = broadcastify(np.log10)
+log2 = broadcastify(np.log2)
+log1p = broadcastify(np.log1p)
+logaddexp = broadcastify(np.logaddexp)
+logaddexp2 = broadcastify(np.logaddexp2)
+
+# Floating functions
+
+signbit = broadcastify(np.signbit)
+copysign = broadcastify(np.copysign)
+# TODO: understand why frexp fails
+frexp = broadcastify(np.frexp)
+ldexp = broadcastify(np.ldexp)
+
+# Arithmetic operations
+
+# add = broadcastify(np.add)
+# reciprocal = broadcastify(np.reciprocal)
+# positive = broadcastify(np.positive)
+# negative = broadcastify(np.negative)
+# multiply = broadcastify(np.multiply)
+# divide = broadcastify(np.divide)
+# power = broadcastify(np.power)
+# subtract = broadcastify(np.subtract)
+# true_divide = broadcastify(np.true_divide)
+# floor_divide = broadcastify(np.floor_divide)
+# fmod = broadcastify(np.fmod)
+# mod = broadcastify(np.mod)
+# modf = broadcastify(np.modf)
+# remainder = broadcastify(np.remainder)
+# divmod = broadcastify(np.divmod)
+
+# Handling complex numbers
+
+conj = broadcastify(np.conj)
+
+# Miscellaneous
+
+sqrt = broadcastify(np.sqrt)
+square = broadcastify(np.square)
+absolute = broadcastify(np.absolute)
+fabs = broadcastify(np.fabs)
+sign = broadcastify(np.sign)
+
+
+def maximum(x1, x2, out=None, dtype=None):
+    r"""
     Element-wise maximum of array elements.
 
     Compare two arrays and returns a new array containing the element-wise
@@ -203,7 +292,7 @@ def broadcastify(func):
 
     Examples
     --------
-    >>> from larray import LArray
+    >>> from larray import LArray, zeros_like
     >>> arr1 = LArray([[10, 7, 5, 9],
     ...                [5, 8, 3, 7]], "a=a0,a1;b=b0..b3")
     >>> arr2 = LArray([[6, 2, 9, 0],
@@ -231,12 +320,12 @@ def broadcastify(func):
     a\b  b0  b1  b2  b3
      a0  10   7   9   9
      a1   6   8   9   7
-"""
+    """
+    return np.maximum(x1, x2, out=out, dtype=dtype)
 
-minimum = broadcastify(np.minimum)
-minimum.__doc__ = r"""
-minimum(x1, x2, out=None, dtype=None)
 
+def minimum(x1, x2, out=None, dtype=None):
+    r"""
     Element-wise minimum of array elements.
 
     Compare two arrays and returns a new array containing the element-wise
@@ -301,4 +390,124 @@ def broadcastify(func):
     a\b  b0  b1  b2  b3
      a0   6   2   5   0
      a1   5   2   3   0
+    """
+    return np.minimum(x1, x2, out=out, dtype=dtype)
+
+fmax = broadcastify(np.fmax)
+fmin = broadcastify(np.fmin)
+isnan = broadcastify(np.isnan)
+isinf = broadcastify(np.isinf)
+
+
+# --------------------------------
+# numpy funcs which are not ufuncs
+
+
+# Rounding
+
+# all 3 are equivalent, I am unsure I should support around and round_
+around = broadcastify(np.around)
+round_ = broadcastify(np.round_)
+round = broadcastify(np.round)
+
+# Sums, products, differences
+
+# prod = broadcastify(np.prod)
+# sum = broadcastify(np.sum)
+# nansum = broadcastify(np.nansum)
+# cumprod = broadcastify(np.cumprod)
+# cumsum = broadcastify(np.cumsum)
+
+# cannot use a simple wrapped ufunc because those ufuncs do not preserve
+# shape or dimensions so labels are wrong
+# diff = broadcastify(np.diff)
+# ediff1d = broadcastify(np.ediff1d)
+# gradient = broadcastify(np.gradient)
+# cross = broadcastify(np.cross)
+# trapz = broadcastify(np.trapz)
+
+# Other special functions
+
+i0 = broadcastify(np.i0)
+sinc = broadcastify(np.sinc)
+
+# Handling complex numbers
+
+angle = broadcastify(np.angle)
+real = broadcastify(np.real)
+imag = broadcastify(np.imag)
+
+# Miscellaneous
+
+# put clip here even it is a ufunc because it ends up within a recursion loop with LArray.clip
+clip = broadcastify(np.clip)
+
+where = broadcastify(np.where, False)
+where.__doc__ = r"""
+    Return elements, either from `x` or `y`, depending on `condition`.
+
+    If only `condition` is given, return ``condition.nonzero()``.
+
+    Parameters
+    ----------
+    condition : boolean LArray
+        When True, yield `x`, otherwise yield `y`.
+    x, y : LArray
+        Values from which to choose.
+
+    Returns
+    -------
+    out : LArray
+        If both `x` and `y` are specified, the output array contains
+        elements of `x` where `condition` is True, and elements from
+        `y` elsewhere.
+
+    Examples
+    --------
+    >>> from larray import LArray
+    >>> arr = LArray([[10, 7, 5, 9],
+    ...               [5, 8, 3, 7],
+    ...               [6, 2, 0, 9],
+    ...               [9, 10, 5, 6]], "a=a0..a3;b=b0..b3")
+    >>> arr
+    a\b  b0  b1  b2  b3
+     a0  10   7   5   9
+     a1   5   8   3   7
+     a2   6   2   0   9
+     a3   9  10   5   6
+
+    Simple use
+
+    >>> where(arr <= 5, 0, arr)
+    a\b  b0  b1  b2  b3
+     a0  10   7   0   9
+     a1   0   8   0   7
+     a2   6   0   0   9
+     a3   9  10   0   6
+
+    With broadcasting
+
+    >>> mean_by_col = arr.mean('a')
+    >>> mean_by_col
+    b   b0    b1    b2    b3
+       7.5  6.75  3.25  7.75
+    >>> # for each column, set values below the mean value to the mean value
+    >>> where(arr < mean_by_col, mean_by_col, arr)
+    a\b    b0    b1    b2    b3
+     a0  10.0   7.0   5.0   9.0
+     a1   7.5   8.0  3.25  7.75
+     a2   7.5  6.75  3.25   9.0
+     a3   9.0  10.0   5.0  7.75
 """
+
+nan_to_num = broadcastify(np.nan_to_num)
+real_if_close = broadcastify(np.real_if_close)
+convolve = broadcastify(np.convolve)
+inverse = broadcastify(np.linalg.inv)
+
+# XXX: create a new LArray method instead ?
+# TODO: should appear in the API doc if it actually works with LArrays,
+#       which I have never tested (and I doubt is the case).
+#       Might be worth having specific documentation if it works well.
+#       My guess is that we should rather make a new LArray method for that one.
+interp = broadcastify(np.interp)
diff --git a/larray/tests/test_array.py b/larray/tests/test_array.py
index 0615f987e..41c63a565 100644
--- a/larray/tests/test_array.py
+++ b/larray/tests/test_array.py
@@ -4324,79 +4324,6 @@ def test_open_excel(tmpdir):
         os.remove(fpath)
 
 
-def test_ufuncs(small_array):
-    raw = small_array.data
-
-    # simple one-argument ufunc
-    assert_array_equal(exp(small_array), np.exp(raw))
-
-    # with out=
-    la_out = zeros(small_array.axes)
-    raw_out = np.zeros(raw.shape)
-
-    la_out2 = exp(small_array, la_out)
-    raw_out2 = np.exp(raw, raw_out)
-
-    # FIXME: this is not the case currently
-    # self.assertIs(la_out2, la_out)
-    assert_array_equal(la_out2, la_out)
-    assert raw_out2 is raw_out
-
-    assert_array_equal(la_out, raw_out)
-
-    # with out= and broadcasting
-    # we need to put the 'a' axis first because array numpy only supports that
-    la_out = zeros([Axis([0, 1, 2], 'a')] + list(small_array.axes))
-    raw_out = np.zeros((3,) + raw.shape)
-
-    la_out2 = exp(small_array, la_out)
-    raw_out2 = np.exp(raw, raw_out)
-
-    # self.assertIs(la_out2, la_out)
-    # XXX: why is la_out2 transposed?
-    assert_array_equal(la_out2.transpose(X.a), la_out)
-    assert raw_out2 is raw_out
-
-    assert_array_equal(la_out, raw_out)
-
-    sex, lipro = small_array.axes
-
-    low = small_array.sum(sex) // 4 + 3
-    raw_low = raw.sum(0) // 4 + 3
-    high = small_array.sum(sex) // 4 + 13
-    raw_high = raw.sum(0) // 4 + 13
-
-    # LA + scalars
-    assert_array_equal(small_array.clip(0, 10), raw.clip(0, 10))
-    assert_array_equal(clip(small_array, 0, 10), np.clip(raw, 0, 10))
-
-    # LA + LA (no broadcasting)
-    assert_array_equal(clip(small_array, 21 - small_array, 9 + small_array // 2),
-                       np.clip(raw, 21 - raw, 9 + raw // 2))
-
-    # LA + LA (with broadcasting)
-    assert_array_equal(clip(small_array, low, high),
-                       np.clip(raw, raw_low, raw_high))
-
-    # where (no broadcasting)
-    assert_array_equal(where(small_array < 5, -5, small_array),
-                       np.where(raw < 5, -5, raw))
-
-    # where (transposed no broadcasting)
-    assert_array_equal(where(small_array < 5, -5, small_array.T),
-                       np.where(raw < 5, -5, raw))
-
-    # where (with broadcasting)
-    result = where(small_array['P01'] < 5, -5, small_array)
-    assert result.axes.names == ['sex', 'lipro']
-    assert_array_equal(result, np.where(raw[:, [0]] < 5, -5, raw))
-
-    # round
-    small_float = small_array + 0.6
-    rounded = round(small_float)
-    assert_array_equal(rounded, np.round(raw + 0.6))
-
-
 def test_diag():
     # 2D -> 1D
     a = ndtest((3, 3))
diff --git a/larray/tests/test_numpy_funcs.py b/larray/tests/test_numpy_funcs.py
new file mode 100644
index 000000000..9d923fe22
--- /dev/null
+++ b/larray/tests/test_numpy_funcs.py
@@ -0,0 +1,87 @@
+from __future__ import absolute_import, division, print_function
+
+import numpy as np
+
+from larray.tests.common import assert_array_equal
+from larray.tests.test_array import small_array
+from larray import Axis, zeros, exp, clip, where, round, cos
+
+
+def test_unary():
+    args = [0, np.zeros(2), zeros((2, 2))]
+    for a in args:
+        assert_array_equal(a + 1, np.cos(a))
+        assert_array_equal(a + 1, cos(a))
+
+
+def test_ufuncs(small_array):
+    raw = small_array.data
+
+    # simple one-argument ufunc
+    assert_array_equal(exp(small_array), np.exp(raw))
+
+    # with out=
+    la_out = zeros(small_array.axes)
+    raw_out = np.zeros(raw.shape)
+
+    la_out2 = exp(small_array, la_out)
+    raw_out2 = np.exp(raw, raw_out)
+
+    # FIXME: this is not the case currently
+    # self.assertIs(la_out2, la_out)
+    assert_array_equal(la_out2, la_out)
+    assert raw_out2 is raw_out
+
+    assert_array_equal(la_out, raw_out)
+
+    # with out= and broadcasting
+    # we need to put the 'a' axis first because array numpy only supports that
+    la_out = zeros([Axis([0, 1, 2], 'a')] + list(small_array.axes))
+    raw_out = np.zeros((3,) + raw.shape)
+
+    la_out2 = exp(small_array, la_out)
+    raw_out2 = np.exp(raw, raw_out)
+
+    # self.assertIs(la_out2, la_out)
+    # XXX: why is la_out2 transposed?
+    assert_array_equal(la_out2.transpose('a'), la_out)
+    assert raw_out2 is raw_out
+
+    assert_array_equal(la_out, raw_out)
+
+    sex, lipro = small_array.axes
+
+    low = small_array.sum(sex) // 4 + 3
+    raw_low = raw.sum(0) // 4 + 3
+    high = small_array.sum(sex) // 4 + 13
+    raw_high = raw.sum(0) // 4 + 13
+
+    # LA + scalars
+    assert_array_equal(small_array.clip(0, 10), raw.clip(0, 10))
+    assert_array_equal(clip(small_array, 0, 10), np.clip(raw, 0, 10))
+
+    # LA + LA (no broadcasting)
+    assert_array_equal(clip(small_array, 21 - small_array, 9 + small_array // 2),
+                       np.clip(raw, 21 - raw, 9 + raw // 2))
+
+    # LA + LA (with broadcasting)
+    assert_array_equal(clip(small_array, low, high),
+                       np.clip(raw, raw_low, raw_high))
+
+    # where (no broadcasting)
+    assert_array_equal(where(small_array < 5, -5, small_array),
+                       np.where(raw < 5, -5, raw))
+
+    # where (transposed no broadcasting)
+    assert_array_equal(where(small_array < 5, -5, small_array.T),
+                       np.where(raw < 5, -5, raw))
+
+    # where (with broadcasting)
+    result = where(small_array['P01'] < 5, -5, small_array)
+    assert result.axes.names == ['sex', 'lipro']
+    assert_array_equal(result, np.where(raw[:, [0]] < 5, -5, raw))
+
+    # round
+    small_float = small_array + 0.6
+    rounded = round(small_float)
+    assert_array_equal(rounded, np.round(raw + 0.6))
diff --git a/setup.cfg b/setup.cfg
index 3aa0093e2..acc070d8f 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -7,10 +7,43 @@ testpaths = larray
 #   and many of those doctests are failing
 # - deselect LArray.astype since doctests fails for Python 3.6 and numpy >= 1.17
 addopts = -v --doctest-modules
-          --ignore=larray/core/npufuncs.py
           --ignore=larray/ipfp
           --ignore=larray/inout/xw_reporting.py
           --deselect larray/core/array.py::larray.core.array.LArray.astype
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.angle
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.arccos
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.arctan
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.cos
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.cosh
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.frexp
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.ldexp
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.log10
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.log2
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.signbit
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.sin
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.sinh
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.tan
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.tanh
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.interp
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.sinc
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.unwrap
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.absolute
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.exp
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.fmax
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.fmin
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.inverse
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.isinf
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.isnan
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.log
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.sqrt
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.square
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.around
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.convolve
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.i0
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.imag
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.real
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.real_if_close
+          --deselect larray/core/ufuncs.py::larray.core.ufuncs.nan_to_num
           --pep8
           #--cov
 # E122: continuation line missing indentation or outdented

From da29162a97e58ffe3b084ec3627abf537d0e38c7 Mon Sep 17 00:00:00 2001
From: Alix Damman <ald@plan.be>
Date: Thu, 16 May 2019 08:59:14 +0200
Subject: [PATCH 3/3] set and updated documentation and doctests of functions
 derived in the ufuncs.py module

---
 larray/core/ufuncs.py | 943 +++++++++++++++++++++++++++++++++++++++---
 setup.cfg             |  10 -
 2 files changed, 883 insertions(+), 70 deletions(-)

diff --git a/larray/core/ufuncs.py b/larray/core/ufuncs.py
index bce68ea41..0ccf9b819 100644
--- a/larray/core/ufuncs.py
+++ b/larray/core/ufuncs.py
@@ -6,6 +6,8 @@
 
 import numpy as np
 
+from larray.tests.common import needs_python36
+
 
 def wrap_elementwise_array_func(func):
     r"""
@@ -197,18 +199,261 @@ def broadcastify(npfunc, copy_numpy_doc=True):
 rint = broadcastify(np.rint)
 # TODO: fix fails because of its 'out' argument
 fix = broadcastify(np.fix)
-# TODO: add examples for round, floor, ceil and trunc
-floor = broadcastify(np.floor)
-ceil = broadcastify(np.ceil)
-trunc = broadcastify(np.trunc)
+
+
+def floor(x, out=None, dtype=None, where=True):
+    r"""
+    Return the floor of the input, element-wise.
+
+    The floor of the scalar `x` is the largest integer `i`, such that
+    `i <= x`.  It is often denoted as :math:`\lfloor x \rfloor`.
+
+    Parameters
+    ----------
+    x : LArray
+        Input array.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+    where : array_like, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    y : LArray or scalar
+        The floor of each element in `x`.
+        This is a scalar if `x` is a scalar.
+
+    See Also
+    --------
+    ceil, trunc, rint
+
+    Notes
+    -----
+    Some spreadsheet programs calculate the "floor-towards-zero", in other
+    words ``floor(-2.5) == -2``.  LArray instead uses the definition of
+    `floor` where `floor(-2.5) == -3`.
+
+    Examples
+    --------
+    >>> from larray import LArray
+    >>> arr = LArray([-1.7, -1.5, -0.2, 0.2, 1.5, 1.7, 2.0], 'a=a0..a6')
+    >>> floor(arr)
+    a    a0    a1    a2   a3   a4   a5   a6
+       -2.0  -2.0  -1.0  0.0  1.0  1.0  2.0
+    """
+    return np.floor(x, out=out, dtype=dtype, where=where)
+
+
+def ceil(x, out=None, dtype=None, where=True):
+    r"""
+    Return the ceiling of the input, element-wise.
+
+    The ceil of the scalar `x` is the smallest integer `i`, such that
+    `i >= x`.  It is often denoted as :math:`\lceil x \rceil`.
+
+    Parameters
+    ----------
+    x : LArray
+        Input array.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+    where : array_like, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    y : LArray or scalar
+        The ceiling of each element in `x`, with `float` dtype.
+        This is a scalar if `x` is a scalar.
+
+    See Also
+    --------
+    floor, trunc, rint
+
+    Examples
+    --------
+    >>> from larray import LArray
+    >>> arr = LArray([-1.7, -1.5, -0.2, 0.2, 1.5, 1.7, 2.0], 'a=a0..a6')
+    >>> ceil(arr)
+    a    a0    a1    a2   a3   a4   a5   a6
+       -1.0  -1.0  -0.0  1.0  2.0  2.0  2.0
+    """
+    return np.ceil(x, out=out, dtype=dtype, where=where)
+
+
+def trunc(x, out=None, dtype=None, where=True):
+    r"""
+    Return the truncated value of the input, element-wise.
+
+    The truncated value of the scalar `x` is the nearest integer `i` which
+    is closer to zero than `x` is. In short, the fractional part of the
+    signed number `x` is discarded.
+
+    Parameters
+    ----------
+    x : LArray
+        Input array.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+    where : array_like, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    y : LArray or scalar
+        The truncated value of each element in `x`.
+        This is a scalar if `x` is a scalar.
+
+    See Also
+    --------
+    ceil, floor, rint
+
+    Examples
+    --------
+    >>> from larray import LArray
+    >>> arr = LArray([-1.7, -1.5, -0.2, 0.2, 1.5, 1.7, 2.0], 'a=a0..a6')
+    >>> trunc(arr)
+    a    a0    a1    a2   a3   a4   a5   a6
+       -1.0  -1.0  -0.0  0.0  1.0  1.0  2.0
+    """
+    return np.trunc(x, out=out, dtype=dtype, where=where)
+
 
 # Exponents and logarithms
 
-# TODO: add examples for exp and log
-exp = broadcastify(np.exp)
+
+def exp(x, out=None, dtype=None, where=True):
+    r"""
+    Calculate the exponential of all elements in the input array.
+
+    Parameters
+    ----------
+    x : LArray
+        Input array.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+    where : array_like, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    out : LArray or scalar
+        Output array, element-wise exponential of `x`.
+        This is a scalar if `x` is a scalar.
+
+    See Also
+    --------
+    expm1 : Calculate ``exp(x) - 1`` for all elements in the array.
+    exp2  : Calculate ``2**x`` for all elements in the array.
+
+    Notes
+    -----
+    The irrational number ``e`` is also known as Euler's number.  It is
+    approximately 2.718281, and is the base of the natural logarithm,
+    ``ln`` (this means that, if :math:`x = \ln y = \log_e y`,
+    then :math:`e^x = y`. For real input, ``exp(x)`` is always positive.
+
+    For complex arguments, ``x = a + ib``, we can write
+    :math:`e^x = e^a e^{ib}`.  The first term, :math:`e^a`, is already
+    known (it is the real argument, described above).  The second term,
+    :math:`e^{ib}`, is :math:`\cos b + i \sin b`, a function with
+    magnitude 1 and a periodic phase.
+
+    References
+    ----------
+    .. [1] Wikipedia, "Exponential function",
+           https://en.wikipedia.org/wiki/Exponential_function
+    .. [2] M. Abramovitz and I. A. Stegun, "Handbook of Mathematical Functions
+           with Formulas, Graphs, and Mathematical Tables," Dover, 1964, p. 69,
+           http://www.math.sfu.ca/~cbm/aands/page_69.htm
+
+    Examples
+    --------
+    >>> from larray import LArray, e
+    >>> exp(LArray([-1, 0, 1])) / LArray([1/e, 1, e])
+    {0}*    0    1    2
+          1.0  1.0  1.0
+    """
+    return np.exp(x, out=out, dtype=dtype, where=where)
+
+
 exp2 = broadcastify(np.exp2)
 expm1 = broadcastify(np.expm1)
-log = broadcastify(np.log)
+
+
+def log(x, out=None, dtype=None, where=True):
+    r"""
+    Natural logarithm, element-wise.
+
+    The natural logarithm `log` is the inverse of the exponential function,
+    so that `log(exp(x)) = x`. The natural logarithm is logarithm in base
+    `e`.
+
+    Parameters
+    ----------
+    x : LArray
+        Input array.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+    where : array_like, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    y : LArray
+        The natural logarithm of `x`, element-wise.
+        This is a scalar if `x` is a scalar.
+
+    See Also
+    --------
+    log10, log2, log1p
+
+    Notes
+    -----
+    Logarithm is a multivalued function: for each `x` there is an infinite
+    number of `z` such that `exp(z) = x`. The convention is to return the
+    `z` whose imaginary part lies in `[-pi, pi]`.
+
+    For real-valued input data types, `log` always returns real output. For
+    each value that cannot be expressed as a real number or infinity, it
+    yields ``nan`` and sets the `invalid` floating point error flag.
+
+    For complex-valued input, `log` is a complex analytical function that
+    has a branch cut `[-inf, 0]` and is continuous from above on it. `log`
+    handles the floating-point negative zero as an infinitesimal negative
+    number, conforming to the C99 standard.
+
+    References
+    ----------
+    .. [1] M. Abramowitz and I.A. Stegun, "Handbook of Mathematical Functions",
+           10th printing, 1964, pp. 67. http://www.math.sfu.ca/~cbm/aands/
+    .. [2] Wikipedia, "Logarithm". https://en.wikipedia.org/wiki/Logarithm
+
+    Examples
+    --------
+    >>> from larray import LArray, e
+    >>> log(LArray([1, e, e**2, 0]))
+    {0}*    0    1    2     3
+          0.0  1.0  2.0  -inf
+    """
+    return np.log(x, out=out, dtype=dtype, where=where)
+
+
 log10 = broadcastify(np.log10)
 log2 = broadcastify(np.log2)
 log1p = broadcastify(np.log1p)
@@ -219,7 +464,6 @@ def broadcastify(npfunc, copy_numpy_doc=True):
 
 signbit = broadcastify(np.signbit)
 copysign = broadcastify(np.copysign)
-# TODO: understand why frexp fails
 frexp = broadcastify(np.frexp)
 ldexp = broadcastify(np.ldexp)
 
@@ -247,11 +491,220 @@ def broadcastify(npfunc, copy_numpy_doc=True):
 
 # Miscellaneous
 
-sqrt = broadcastify(np.sqrt)
-square = broadcastify(np.square)
-absolute = broadcastify(np.absolute)
-fabs = broadcastify(np.fabs)
-sign = broadcastify(np.sign)
+
+def sqrt(x, out=None, dtype=None, where=True):
+    r"""
+    Return the non-negative square-root of an array, element-wise.
+
+    Parameters
+    ----------
+    x : LArray
+        The values whose square-roots are required.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+    where : array_like, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    y : LArray
+        An array of the same shape as `x`, containing the positive
+        square-root of each element in `x`.  If any element in `x` is
+        complex, a complex array is returned (and the square-roots of
+        negative reals are calculated).  If all of the elements in `x`
+        are real, so is `y`, with negative elements returning ``nan``.
+        If `out` was provided, `y` is a reference to it.
+        This is a scalar if `x` is a scalar.
+
+    Notes
+    -----
+    *sqrt* has--consistent with common convention--as its branch cut the
+    real "interval" [`-inf`, 0), and is continuous from above on it.
+    A branch cut is a curve in the complex plane across which a given
+    complex function fails to be continuous.
+
+    Examples
+    --------
+    >>> from larray import LArray, inf
+    >>> sqrt(LArray([1, 4, 9], 'a=a0..a2'))
+    a   a0   a1   a2
+       1.0  2.0  3.0
+
+    >>> sqrt(LArray([4, -1, -3+4J], 'a=a0..a2'))
+    a      a0  a1      a2
+       (2+0j)  1j  (1+2j)
+
+    >>> sqrt(LArray([4, -1, inf], 'a=a0..a2'))
+    a   a0   a1   a2
+       2.0  nan  inf
+    """
+    return np.sqrt(x, out=out, dtype=dtype, where=where)
+
+
+def square(x, out=None, dtype=None, where=True):
+    r"""
+    Return the element-wise square of the input.
+
+    Parameters
+    ----------
+    x : LArray
+        The values whose square are required.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+    where : array_like, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    out : LArray or scalar
+        Element-wise `x*x`, of the same shape and dtype as `x`.
+        This is a scalar if `x` is a scalar.
+
+    See Also
+    --------
+    sqrt
+
+    Examples
+    --------
+    >>> from larray import LArray, nan, inf
+    >>> square(LArray([2., 3., nan, inf], 'a=a0..a3'))
+    a   a0   a1   a2   a3
+       4.0  9.0  nan  inf
+    """
+    return np.square(x, out=out, dtype=dtype, where=where)
+
+
+def absolute(x, out=None, dtype=None, where=True):
+    r"""
+    Calculate the absolute value element-wise.
+
+    Parameters
+    ----------
+    x : LArray
+        Input array.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+    where : array_like, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    absolute : LArray
+        An array containing the absolute value of
+        each element in `x`.  For complex input, ``a + ib``, the
+        absolute value is :math:`\sqrt{ a^2 + b^2 }`.
+        This is a scalar if `x` is a scalar.
+
+    See Also
+    --------
+    fabs
+
+    Examples
+    --------
+    >>> from larray import LArray, nan, inf
+    >>> absolute(LArray([1, -1, -3+4J, nan, -inf, inf], 'a=a0..a5'))
+    a   a0   a1   a2   a3   a4   a5
+       1.0  1.0  5.0  nan  inf  inf
+    """
+    return np.absolute(x, out=out, dtype=dtype, where=where)
+
+
+def fabs(x, out=None, dtype=None, where=True):
+    r"""
+    Compute the absolute values element-wise.
+
+    This function returns the absolute values (positive magnitude) of the
+    data in `x`. Complex values are not handled, use `absolute` to find the
+    absolute values of complex data.
+
+    Parameters
+    ----------
+    x : LArray
+        Input array.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+    where : array_like, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    y : LArray or scalar
+        The absolute values of `x`, the returned values are always floats.
+        This is a scalar if `x` is a scalar.
+
+    See Also
+    --------
+    absolute : Absolute values including `complex` types.
+
+    Examples
+    --------
+    >>> from larray import LArray, nan, inf
+    >>> fabs(LArray([1, -1, nan, -inf, inf], 'a=a0..a4'))
+    a   a0   a1   a2   a3   a4
+       1.0  1.0  nan  inf  inf
+    """
+    return np.fabs(x, out=out, dtype=dtype, where=where)
+
+
+def sign(x, out=None, dtype=None, where=True):
+    r"""
+    Returns an element-wise indication of the sign of a number.
+
+    The `sign` function returns ``-1 if x < 0, 0 if x==0, 1 if x > 0``.
+    nan is returned for nan inputs.
+
+    For complex inputs, the `sign` function returns
+    ``sign(x.real) + 0j if x.real != 0 else sign(x.imag) + 0j``.
+
+    complex(nan, 0) is returned for complex nan inputs.
+
+    Parameters
+    ----------
+    x : LArray
+        Input array.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+    where : array_like, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    y : LArray or scalar
+        The sign of `x`.
+        This is a scalar if `x` is a scalar.
+
+    Notes
+    -----
+    There is more than one definition of sign in common use for complex
+    numbers.  The definition used here is equivalent to :math:`x/\sqrt{x*x}`
+    which is different from a common alternative, :math:`x/|x|`.
+
+    Examples
+    --------
+    >>> from larray import LArray, nan, inf
+    >>> sign(LArray([-5., 4.5, 0, -inf, inf, nan], 'a=a0..a5'))
+    a    a0   a1   a2    a3   a4   a5
+       -1.0  1.0  0.0  -1.0  1.0  nan
+    >>> np.sign(LArray([5, 5-2J, 2J], 'a=a0..a2'))
+    a      a0      a1      a2
+       (1+0j)  (1+0j)  (1+0j)
+    """
+    return np.sign(x, out=out, dtype=dtype, where=where)
 
 
 def maximum(x1, x2, out=None, dtype=None):
@@ -284,6 +737,10 @@ def maximum(x1, x2, out=None, dtype=None):
     --------
     minimum :
         Element-wise minimum of two arrays, propagates NaNs.
+    fmax :
+        Element-wise maximum of two arrays, ignores NaNs.
+    fmin :
+        Element-wise minimum of two arrays, ignores NaNs.
 
     Notes
     -----
@@ -292,34 +749,34 @@ def maximum(x1, x2, out=None, dtype=None):
 
     Examples
     --------
-    >>> from larray import LArray, zeros_like
-    >>> arr1 = LArray([[10, 7, 5, 9],
-    ...                [5, 8, 3, 7]], "a=a0,a1;b=b0..b3")
-    >>> arr2 = LArray([[6, 2, 9, 0],
-    ...                [9, 10, 5, 6]], "a=a0,a1;b=b0..b3")
+    >>> from larray import LArray, nan
+    >>> arr1 = LArray([[8.0, 7.0, 5.0, 6.0],
+    ...                [5.0, 8.0, 3.0, nan]], "a=a0,a1;b=b0..b3")
+    >>> arr2 = LArray([[10.0, 2.0, 9.0, nan],
+    ...                [9.0, 10.0, 5.0, nan]], "a=a0,a1;b=b0..b3")
     >>> arr1
-    a\b  b0  b1  b2  b3
-     a0  10   7   5   9
-     a1   5   8   3   7
+    a\b   b0   b1   b2   b3
+     a0  8.0  7.0  5.0  6.0
+     a1  5.0  8.0  3.0  nan
     >>> arr2
-    a\b  b0  b1  b2  b3
-     a0   6   2   9   0
-     a1   9  10   5   6
+    a\b    b0    b1   b2   b3
+     a0  10.0   2.0  9.0  nan
+     a1   9.0  10.0  5.0  nan
 
     >>> maximum(arr1, arr2)
-    a\b  b0  b1  b2  b3
-     a0  10   7   9   9
-     a1   9  10   5   7
+    a\b    b0    b1   b2   b3
+     a0  10.0   7.0  9.0  nan
+     a1   9.0  10.0  5.0  nan
 
     With broadcasting
 
     >>> arr2['a0']
-    b  b0  b1  b2  b3
-        6   2   9   0
+    b    b0   b1   b2   b3
+       10.0  2.0  9.0  nan
     >>> maximum(arr1, arr2['a0'])
-    a\b  b0  b1  b2  b3
-     a0  10   7   9   9
-     a1   6   8   9   7
+    a\b    b0   b1   b2   b3
+     a0  10.0  7.0  9.0  nan
+     a1  10.0  8.0  9.0  nan
     """
     return np.maximum(x1, x2, out=out, dtype=dtype)
 
@@ -354,6 +811,10 @@ def minimum(x1, x2, out=None, dtype=None):
     --------
     maximum :
         Element-wise maximum of two arrays, propagates NaNs.
+    fmax :
+        Element-wise maximum of two arrays, ignores NaNs.
+    fmin :
+        Element-wise minimum of two arrays, ignores NaNs.
 
     Notes
     -----
@@ -362,41 +823,279 @@ def minimum(x1, x2, out=None, dtype=None):
 
     Examples
     --------
-    >>> from larray import LArray
-    >>> arr1 = LArray([[10, 7, 5, 9],
-    ...                [5, 8, 3, 7]], "a=a0,a1;b=b0..b3")
-    >>> arr2 = LArray([[6, 2, 9, 0],
-    ...                [9, 10, 5, 6]], "a=a0,a1;b=b0..b3")
+    >>> from larray import LArray, nan
+    >>> arr1 = LArray([[8.0, 7.0, 5.0, 6.0],
+    ...                [5.0, 8.0, 3.0, nan]], "a=a0,a1;b=b0..b3")
+    >>> arr2 = LArray([[10.0, 2.0, 9.0, nan],
+    ...                [9.0, 10.0, 5.0, nan]], "a=a0,a1;b=b0..b3")
     >>> arr1
-    a\b  b0  b1  b2  b3
-     a0  10   7   5   9
-     a1   5   8   3   7
+    a\b   b0   b1   b2   b3
+     a0  8.0  7.0  5.0  6.0
+     a1  5.0  8.0  3.0  nan
     >>> arr2
-    a\b  b0  b1  b2  b3
-     a0   6   2   9   0
-     a1   9  10   5   6
+    a\b    b0    b1   b2   b3
+     a0  10.0   2.0  9.0  nan
+     a1   9.0  10.0  5.0  nan
 
     >>> minimum(arr1, arr2)
-    a\b  b0  b1  b2  b3
-     a0   6   2   5   0
-     a1   5   8   3   6
+    a\b   b0   b1   b2   b3
+     a0  8.0  2.0  5.0  nan
+     a1  5.0  8.0  3.0  nan
 
     With broadcasting
 
     >>> arr2['a0']
-    b  b0  b1  b2  b3
-        6   2   9   0
+    b    b0   b1   b2   b3
+       10.0  2.0  9.0  nan
     >>> minimum(arr1, arr2['a0'])
-    a\b  b0  b1  b2  b3
-     a0   6   2   5   0
-     a1   5   2   3   0
+    a\b   b0   b1   b2   b3
+     a0  8.0  2.0  5.0  nan
+     a1  5.0  2.0  3.0  nan
     """
     return np.minimum(x1, x2, out=out, dtype=dtype)
 
-fmax = broadcastify(np.fmax)
-fmin = broadcastify(np.fmin)
-isnan = broadcastify(np.isnan)
-isinf = broadcastify(np.isinf)
+
+def fmax(x1, x2, out=None, dtype=None):
+    r"""
+    Element-wise maximum of array elements.
+
+    Compare two arrays and returns a new array containing the element-wise
+    maxima. If one of the elements being compared is a NaN, then the
+    non-nan element is returned. If both elements are NaNs then the first
+    is returned.  The latter distinction is important for complex NaNs,
+    which are defined as at least one of the real or imaginary parts being
+    a NaN. The net effect is that NaNs are ignored when possible.
+
+    Parameters
+    ----------
+    x1, x2 : LArray
+        The arrays holding the elements to be compared.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+
+    Returns
+    -------
+    y : LArray or scalar
+        The maximum of `x1` and `x2`, element-wise.
+        This is a scalar if both `x1` and `x2` are scalars.
+
+    See Also
+    --------
+    fmin :
+        Element-wise minimum of two arrays, ignores NaNs.
+    maximum :
+        Element-wise maximum of two arrays, propagates NaNs.
+    minimum :
+        Element-wise minimum of two arrays, propagates NaNs.
+
+    Notes
+    -----
+    The fmax is equivalent to ``np.where(x1 >= x2, x1, x2)`` when neither
+    x1 nor x2 are NaNs, but it is faster.
+
+    Examples
+    --------
+    >>> from larray import LArray, nan
+    >>> arr1 = LArray([[8.0, 7.0, 5.0, 6.0],
+    ...                [5.0, 8.0, 3.0, nan]], "a=a0,a1;b=b0..b3")
+    >>> arr2 = LArray([[10.0, 2.0, 9.0, nan],
+    ...                [9.0, 10.0, 5.0, nan]], "a=a0,a1;b=b0..b3")
+    >>> arr1
+    a\b   b0   b1   b2   b3
+     a0  8.0  7.0  5.0  6.0
+     a1  5.0  8.0  3.0  nan
+    >>> arr2
+    a\b    b0    b1   b2   b3
+     a0  10.0   2.0  9.0  nan
+     a1   9.0  10.0  5.0  nan
+
+    >>> fmax(arr1, arr2)
+    a\b    b0    b1   b2   b3
+     a0  10.0   7.0  9.0  6.0
+     a1   9.0  10.0  5.0  nan
+
+    With broadcasting
+
+    >>> arr2['a0']
+    b    b0   b1   b2   b3
+       10.0  2.0  9.0  nan
+    >>> fmax(arr1, arr2['a0'])
+    a\b    b0   b1   b2   b3
+     a0  10.0  7.0  9.0  6.0
+     a1  10.0  8.0  9.0  nan
+    """
+    return np.fmax(x1, x2, out=out, dtype=dtype)
+
+
+def fmin(x1, x2, out=None, dtype=None):
+    r"""
+    Element-wise minimum of array elements.
+
+    Compare two arrays and returns a new array containing the element-wise
+    minima. If one of the elements being compared is a NaN, then the
+    non-nan element is returned. If both elements are NaNs then the first
+    is returned.  The latter distinction is important for complex NaNs,
+    which are defined as at least one of the real or imaginary parts being
+    a NaN. The net effect is that NaNs are ignored when possible.
+
+    Parameters
+    ----------
+    x1, x2 : LArray
+        The arrays holding the elements to be compared.
+    out : LArray, optional
+        An array into which the result is stored.
+    dtype : data-type, optional
+        Overrides the dtype of the output array.
+
+    Returns
+    -------
+    y : LArray or scalar
+        The minimum of `x1` and `x2`, element-wise.
+        This is a scalar if both `x1` and `x2` are scalars.
+
+    See Also
+    --------
+    fmax :
+        Element-wise maximum of two arrays, ignores NaNs.
+    maximum :
+        Element-wise maximum of two arrays, propagates NaNs.
+    minimum :
+        Element-wise minimum of two arrays, propagates NaNs.
+
+    Notes
+    -----
+    The fmin is equivalent to ``np.where(x1 <= x2, x1, x2)`` when neither
+    x1 nor x2 are NaNs, but it is faster.
+
+    Examples
+    --------
+    >>> from larray import LArray, nan
+    >>> arr1 = LArray([[8.0, 7.0, 5.0, 6.0],
+    ...                [5.0, 8.0, 3.0, nan]], "a=a0,a1;b=b0..b3")
+    >>> arr2 = LArray([[10.0, 2.0, 9.0, nan],
+    ...                [9.0, 10.0, 5.0, nan]], "a=a0,a1;b=b0..b3")
+    >>> arr1
+    a\b   b0   b1   b2   b3
+     a0  8.0  7.0  5.0  6.0
+     a1  5.0  8.0  3.0  nan
+    >>> arr2
+    a\b    b0    b1   b2   b3
+     a0  10.0   2.0  9.0  nan
+     a1   9.0  10.0  5.0  nan
+
+    >>> fmin(arr1, arr2)
+    a\b   b0   b1   b2   b3
+     a0  8.0  2.0  5.0  6.0
+     a1  5.0  8.0  3.0  nan
+
+    With broadcasting
+
+    >>> arr2['a0']
+    b    b0   b1   b2   b3
+       10.0  2.0  9.0  nan
+    >>> fmin(arr1, arr2['a0'])
+    a\b   b0   b1   b2   b3
+     a0  8.0  2.0  5.0  6.0
+     a1  5.0  2.0  3.0  nan
+    """
+    return np.fmin(x1, x2, out=out, dtype=dtype)
+
+
+def isnan(x, out=None, where=True):
+    r"""
+    Test element-wise for NaN and return result as a boolean array.
+
+    Parameters
+    ----------
+    x : LArray
+        Input array.
+    out : LArray, optional
+        An array into which the result is stored.
+    where : LArray, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    y : boolean LArray or bool
+        True where ``x`` is NaN, false otherwise.
+        This is a scalar if `x` is a scalar.
+
+    See Also
+    --------
+    isinf
+
+    Notes
+    -----
+    NumPy uses the IEEE Standard for Binary Floating-Point for Arithmetic
+    (IEEE 754). This means that Not a Number is not equivalent to infinity.
+
+    Examples
+    --------
+    >>> from larray import LArray, nan, inf
+    >>> arr1 = LArray([[-inf, 7.0, 5.0, inf],
+    ...                [5.0, 8.0, 3.0, nan]], "a=a0,a1;b=b0..b3")
+    >>> arr1
+    a\b    b0   b1   b2   b3
+     a0  -inf  7.0  5.0  inf
+     a1   5.0  8.0  3.0  nan
+    >>> isnan(arr1)
+    a\b     b0     b1     b2     b3
+     a0  False  False  False  False
+     a1  False  False  False   True
+    """
+    return np.isnan(x, out=out, where=where)
+
+
+def isinf(x, out=None, where=True):
+    r"""
+    Test element-wise for positive or negative infinity.
+
+    Returns a boolean array of the same shape as `x`, True where ``x ==
+    +/-inf``, otherwise False.
+
+    Parameters
+    ----------
+    x : LArray
+        Input array.
+    out : LArray, optional
+        An array into which the result is stored.
+    where : LArray, optional
+        Values of True indicate to calculate the ufunc at that position, values
+        of False indicate to leave the value in the output alone.
+
+    Returns
+    -------
+    y : boolean LArray or bool
+        True where ``x`` is positive or negative infinity, false otherwise.
+        This is a scalar if `x` is a scalar.
+
+    See Also
+    --------
+    isnan
+
+    Notes
+    -----
+    NumPy uses the IEEE Standard for Binary Floating-Point for Arithmetic
+    (IEEE 754). This means that Not a Number is not equivalent to infinity.
+
+    Examples
+    --------
+    >>> from larray import LArray, nan, inf
+    >>> arr1 = LArray([[-inf, 7.0, 5.0, inf],
+    ...                [5.0, 8.0, 3.0, nan]], "a=a0,a1;b=b0..b3")
+    >>> arr1
+    a\b    b0   b1   b2   b3
+     a0  -inf  7.0  5.0  inf
+     a1   5.0  8.0  3.0  nan
+    >>> isinf(arr1)
+    a\b     b0     b1     b2     b3
+     a0   True  False  False   True
+     a1  False  False  False  False
+    """
+    return np.isinf(x, out=out, where=where)
 
 
 # --------------------------------
@@ -408,7 +1107,76 @@ def minimum(x1, x2, out=None, dtype=None):
 # all 3 are equivalent, I am unsure I should support around and round_
 around = broadcastify(np.around)
 round_ = broadcastify(np.round_)
-round = broadcastify(np.round)
+
+
+def round(a, decimals=0, out=None):
+    r"""
+    Round an array to the given number of decimals.
+
+    Parameters
+    ----------
+    a : LArray
+        Input array.
+    decimals : int, optional
+        Number of decimal places to round to (default: 0).
+        If decimals is negative, it specifies the number of positions to
+        the left of the decimal point.
+    out : LArray, optional
+        An array into which the result is stored.
+
+    Returns
+    -------
+    rounded_array : LArray
+        An array of the same type as `a`, containing the rounded values.
+        Unless `out` was specified, a new array is created.
+        A reference to the result is returned.
+
+        The real and imaginary parts of complex numbers are rounded
+        separately.  The result of rounding a float is a float.
+
+    See Also
+    --------
+    ceil, fix, floor, rint, trunc
+
+    Notes
+    -----
+    For values exactly halfway between rounded decimal values, NumPy
+    rounds to the nearest even value. Thus 1.5 and 2.5 round to 2.0,
+    -0.5 and 0.5 round to 0.0, etc. Results may also be surprising due
+    to the inexact representation of decimal fractions in the IEEE
+    floating point standard [1]_ and errors introduced when scaling
+    by powers of ten.
+
+    References
+    ----------
+    .. [1] "Lecture Notes on the Status of IEEE 754", William Kahan,
+           https://people.eecs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF
+    .. [2] "How Futile are Mindless Assessments of
+           Roundoff in Floating-Point Computation?", William Kahan,
+           https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf
+
+    Examples
+    --------
+    >>> from larray import LArray
+    >>> round(LArray([0.37, 1.64], 'a=a0,a1'))
+    a   a0   a1
+       0.0  2.0
+    >>> round(LArray([0.37, 1.64], 'a=a0,a1'), decimals=1)
+    a   a0   a1
+       0.4  1.6
+    >>> round(LArray([.5, 1.5, 2.5, 3.5, 4.5], 'a=a0..a4')) # rounds to nearest even value
+    a   a0   a1   a2   a3   a4
+       0.0  2.0  2.0  4.0  4.0
+    >>> round(LArray([1, 2, 3, 11], 'a=a0..a3'), decimals=1) # array of ints is returned
+    a  a0  a1  a2  a3
+        1   2   3  11
+    >>> round(LArray([1, 2, 3, 11], 'a=a0..a3'), decimals=-1)
+    a  a0  a1  a2  a3
+        0   0   0  10
+    """
+    func = broadcastify(np.round, False)
+    return func(a, decimals, out)
+
 
 # Sums, products, differences
 
@@ -442,8 +1210,9 @@ def minimum(x1, x2, out=None, dtype=None):
 # put clip here even it is a ufunc because it ends up within a recursion loop with LArray.clip
 clip = broadcastify(np.clip)
 
-where = broadcastify(np.where, False)
-where.__doc__ = r"""
+
+def where(condition, x=None, y=None):
+    r"""
     Return elements, either from `x` or `y`, depending on `condition`.
 
     If only `condition` is given, return ``condition.nonzero()``.
@@ -498,12 +1267,66 @@ def minimum(x1, x2, out=None, dtype=None):
      a1   7.5   8.0  3.25  7.75
      a2   7.5  6.75  3.25   9.0
      a3   9.0  10.0   5.0  7.75
-"""
+    """
+    func = broadcastify(np.where, False)
+    return func(condition, x, y)
 
 nan_to_num = broadcastify(np.nan_to_num)
 real_if_close = broadcastify(np.real_if_close)
 convolve = broadcastify(np.convolve)
-inverse = broadcastify(np.linalg.inv)
+
+
+@needs_python36
+def inverse(a):
+    r"""
+    Compute the (multiplicative) inverse of a matrix.
+
+    Given a square matrix `a`, return the matrix `ainv` satisfying
+    ``dot(a, ainv) = dot(ainv, a) = eye(a.shape[0])``.
+
+    Parameters
+    ----------
+    a : Larray
+        Matrix(ces) to be inverted.
+
+    Returns
+    -------
+    ainv : LArray
+        Inverse of the matrix(ces) `a`.
+
+    Raises
+    ------
+    LinAlgError
+        If `a` is not square or inversion fails.
+
+    Examples
+    --------
+    >>> from larray import LArray, eye
+    >>> arr = LArray([[1., 2.],
+    ...               [3., 4.]], "a=a0,a1;b=b0,b1")
+    >>> arr
+    a\b   b0   b1
+     a0  1.0  2.0
+     a1  3.0  4.0
+    >>> arr_inv = inverse(arr)
+    >>> (arr @ arr_inv).equals(eye(arr.a, arr.b), atol=1.e-10)
+    True
+    >>> (arr_inv @ arr).equals(eye(arr.a, arr.b), atol=1.e-10)
+    True
+
+    Inverses of several matrices can be computed at once:
+
+    >>> arr = LArray([[[1., 2.], [3., 4.]],
+    ...               [[1., 3.], [3., 5.]]], "a=a0,a1;b=b0,b1;c=c0,c1")
+    >>> arr_inv = inverse(arr)
+    >>> (arr_inv['a0'] @ arr['a0']).equals(eye(arr.b, arr.c), atol=1.e-10)
+    True
+    >>> (arr_inv['a1'] @ arr['a1']).equals(eye(arr.b, arr.c), atol=1.e-10)
+    True
+    """
+    func = broadcastify(np.linalg.inv, False)
+    return func(a)
+
 
 # XXX: create a new LArray method instead ?
 # TODO: should appear in the API doc if it actually works with LArrays,
diff --git a/setup.cfg b/setup.cfg
index acc070d8f..33acd9ef1 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -27,16 +27,6 @@ addopts = -v --doctest-modules
           --deselect larray/core/ufuncs.py::larray.core.ufuncs.interp
           --deselect larray/core/ufuncs.py::larray.core.ufuncs.sinc
           --deselect larray/core/ufuncs.py::larray.core.ufuncs.unwrap
-          --deselect larray/core/ufuncs.py::larray.core.ufuncs.absolute
-          --deselect larray/core/ufuncs.py::larray.core.ufuncs.exp
-          --deselect larray/core/ufuncs.py::larray.core.ufuncs.fmax
-          --deselect larray/core/ufuncs.py::larray.core.ufuncs.fmin
-          --deselect larray/core/ufuncs.py::larray.core.ufuncs.inverse
-          --deselect larray/core/ufuncs.py::larray.core.ufuncs.isinf
-          --deselect larray/core/ufuncs.py::larray.core.ufuncs.isnan
-          --deselect larray/core/ufuncs.py::larray.core.ufuncs.log
-          --deselect larray/core/ufuncs.py::larray.core.ufuncs.sqrt
-          --deselect larray/core/ufuncs.py::larray.core.ufuncs.square
           --deselect larray/core/ufuncs.py::larray.core.ufuncs.around
           --deselect larray/core/ufuncs.py::larray.core.ufuncs.convolve
           --deselect larray/core/ufuncs.py::larray.core.ufuncs.i0