Repository: dask/dask

Issue:
Dropping NumPy 1.10
Increasingly we discover issues with PRs after merging them on NumPy 1.10, which we then either need to fix ourselves or hope the submitter of the PR will kindly fix the issue. While these issues can normally be solved by skipping, adding compat functions, etc., it starts to raise the question: how long do we want to support NumPy 1.10. As of January NumPy 1.14 was relased and has since had 2 patch releases. The list of NumPy 1.15 issue is not too long either (roughly two thirds have already been resolved). To me this is a pretty good case for dropping, but am interested to hear other thoughts.

If we do feel strongly about keeping NumPy 1.10, could we discuss adding it as part of the PR builds to cutdown on strain caused by the issues mentioned above?

Discussion:

Patch:
diff --git a/.gitignore b/.gitignore
index cb1fc67ff..5b3080424 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,6 @@
 *.pyc
 *.egg-info
+dask-worker-space/
 docs/build
 build/
 dist/
diff --git a/.travis.yml b/.travis.yml
index 8ffe4f6ed..a1edcbfd9 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -35,7 +35,7 @@ jobs:
 
     - env:
       - PYTHON=3.4
-      - NUMPY=1.10.4
+      - NUMPY=1.12.1
       - PANDAS=0.19.1
       - *test_and_lint
       - *no_coverage
@@ -45,7 +45,7 @@ jobs:
 
     - env:
       - PYTHON=3.5
-      - NUMPY=1.12.1
+      - NUMPY=1.11.0
       - PANDAS=0.19.2
       - *test_and_lint
       - *no_coverage
diff --git a/dask/array/__init__.py b/dask/array/__init__.py
index 4e43f29b6..67c697118 100644
--- a/dask/array/__init__.py
+++ b/dask/array/__init__.py
@@ -6,7 +6,7 @@ from .core import (Array, block, concatenate, stack, from_array, store,
                    from_delayed, asarray, asanyarray,
                    broadcast_arrays, broadcast_to)
 from .routines import (take, choose, argwhere, where, coarsen, insert,
-                       ravel, roll, unique, squeeze, topk, ptp, diff, ediff1d,
+                       ravel, roll, unique, squeeze, ptp, diff, ediff1d,
                        bincount, digitize, histogram, cov, array, dstack,
                        vstack, hstack, compress, extract, round, count_nonzero,
                        flatnonzero, nonzero, around, isin, isnull, notnull,
@@ -37,7 +37,8 @@ from .reductions import (sum, prod, mean, std, var, any, all, min, max, vnorm,
                          argmin, argmax,
                          nansum, nanmean, nanstd, nanvar, nanmin,
                          nanmax, nanargmin, nanargmax,
-                         cumsum, cumprod)
+                         cumsum, cumprod,
+                         topk, argtopk)
 from .percentile import percentile
 with ignoring(ImportError):
     from .reductions import nanprod, nancumprod, nancumsum
diff --git a/dask/array/chunk.py b/dask/array/chunk.py
index d38fe660b..2879f38e5 100644
--- a/dask/array/chunk.py
+++ b/dask/array/chunk.py
@@ -61,15 +61,14 @@ nanargmax = keepdims_wrapper(np.nanargmax)
 any = keepdims_wrapper(np.any)
 all = keepdims_wrapper(np.all)
 nansum = keepdims_wrapper(np.nansum)
+nanprod = keepdims_wrapper(np.nanprod)
 
 try:
-    from numpy import nanprod, nancumprod, nancumsum
+    from numpy import nancumprod, nancumsum
 except ImportError:  # pragma: no cover
-    nanprod = npcompat.nanprod
     nancumprod = npcompat.nancumprod
     nancumsum = npcompat.nancumsum
 
-nanprod = keepdims_wrapper(nanprod)
 nancumprod = keepdims_wrapper(nancumprod)
 nancumsum = keepdims_wrapper(nancumsum)
 
@@ -176,16 +175,48 @@ except ImportError:  # pragma: no cover
     broadcast_to = npcompat.broadcast_to
 
 
-def topk(k, x):
-    """ Top k elements of an array
-
-    >>> topk(2, np.array([5, 1, 3, 6]))
-    array([6, 5])
+def topk(a, k, axis, keepdims):
+    """Kernel of topk and argtopk.
+    Extract the k largest elements from a on the given axis.
+    If k is negative, extract the -k smallest elements instead.
+    Note that, unlike in the parent function, the returned elements
+    are not sorted internally.
+    """
+    axis = axis[0]
+    if abs(k) >= a.shape[axis]:
+        return a
+    a = np.partition(a, -k, axis=axis)
+    # return a[-k:] if k>0 else a[:-k], on arbitrary axis
+    return a[[
+        (slice(-k, None) if k > 0 else slice(None, -k))
+        if i == axis else slice(None)
+        for i in range(a.ndim)
+    ]]
+
+
+def topk_postprocess(a, k, axis):
+    """Kernel of topk and argtopk.
+    Post-processes the output of topk, sorting the results internally.
+    """
+    a = np.sort(a, axis=axis)
+    if k > 0:
+        # a = a[::-1] on arbitrary axis
+        a = a[[
+            slice(None, None, -1) if i == axis else slice(None)
+            for i in range(a.ndim)
+        ]]
+    return a
+
+
+def argtopk_preprocess(a, idx):
+    """Kernel of argtopk.
+    Preprocess data, by putting it together with its indexes in a recarray
     """
-    # http://stackoverflow.com/a/23734295/616616 by larsmans
-    k = np.minimum(k, len(x))
-    ind = np.argpartition(x, -k)[-k:]
-    return np.sort(x[ind])[::-1]
+    # np.core.records.fromarrays won't work if a and idx don't have the same shape
+    res = np.recarray(a.shape, dtype=[('values', a.dtype), ('idx', idx.dtype)])
+    res.values = a
+    res.idx = idx
+    return res
 
 
 def arange(start, stop, step, length, dtype):
diff --git a/dask/array/core.py b/dask/array/core.py
index b2c44a4e0..4f0c58e25 100644
--- a/dask/array/core.py
+++ b/dask/array/core.py
@@ -1401,12 +1401,19 @@ class Array(DaskMethodsMixin):
             shape = shape[0]
         return reshape(self, shape)
 
-    def topk(self, k):
+    def topk(self, k, axis=-1, split_every=None):
         """The top k elements of an array.
 
         See ``da.topk`` for docstring"""
-        from .routines import topk
-        return topk(k, self)
+        from .reductions import topk
+        return topk(self, k, axis=axis, split_every=split_every)
+
+    def argtopk(self, k, axis=-1, split_every=None):
+        """The indices of the top k elements of an array.
+
+        See ``da.argtopk`` for docstring"""
+        from .reductions import argtopk
+        return argtopk(self, k, axis=axis, split_every=split_every)
 
     def astype(self, dtype, **kwargs):
         """Copy of the array, cast to a specified type.
diff --git a/dask/array/einsumfuncs.py b/dask/array/einsumfuncs.py
index f896610f3..483335a1d 100644
--- a/dask/array/einsumfuncs.py
+++ b/dask/array/einsumfuncs.py
@@ -193,7 +193,7 @@ def _einsum_kernel(*operands, **kwargs):
     return chunk.reshape(chunk.shape + (1,) * ncontract_inds)
 
 
-_einsum_can_optimize = LooseVersion(np.__version__) >= LooseVersion("1.12.0")
+einsum_can_optimize = LooseVersion(np.__version__) >= LooseVersion("1.12.0")
 
 
 @wraps(np.einsum)
@@ -212,7 +212,7 @@ def einsum(*operands, **kwargs):
     if optimize is None:
         optimize = False
 
-    if _einsum_can_optimize and optimize is not False:
+    if einsum_can_optimize and optimize is not False:
         # Avoid computation of dask arrays within np.einsum_path
         # by passing in small numpy arrays broadcasted
         # up to the right shape
@@ -234,7 +234,7 @@ def einsum(*operands, **kwargs):
     kwargs['kernel_dtype'] = einsum_dtype
     kwargs['ncontract_inds'] = ncontract_inds
 
-    if _einsum_can_optimize:
+    if einsum_can_optimize:
         kwargs['optimize'] = optimize
 
     # Update kwargs with atop parameters
diff --git a/dask/array/ma.py b/dask/array/ma.py
index f656eaf97..67186e1ce 100644
--- a/dask/array/ma.py
+++ b/dask/array/ma.py
@@ -10,8 +10,8 @@ from .core import (concatenate_lookup, tensordot_lookup, map_blocks,
                    asanyarray, atop)
 
 
-if LooseVersion(np.__version__) < '1.11.0':
-    raise ImportError("dask.array.ma requires numpy >= 1.11.0")
+if LooseVersion(np.__version__) < '1.11.2':
+    raise ImportError("dask.array.ma requires numpy >= 1.11.2")
 
 
 @normalize_token.register(np.ma.masked_array)
diff --git a/dask/array/numpy_compat.py b/dask/array/numpy_compat.py
index 8f27e5c91..6e1999e84 100644
--- a/dask/array/numpy_compat.py
+++ b/dask/array/numpy_compat.py
@@ -2,29 +2,9 @@ from __future__ import absolute_import, division, print_function
 
 from distutils.version import LooseVersion
 
-from ..compatibility import builtins
 import numpy as np
 import warnings
 
-try:
-    isclose = np.isclose
-except AttributeError:
-    def isclose(*args, **kwargs):
-        raise RuntimeError("You need numpy version 1.7 or greater to use "
-                           "isclose.")
-
-try:
-    full = np.full
-except AttributeError:
-    def full(shape, fill_value, dtype=None, order=None):
-        """Our implementation of numpy.full because your numpy is old."""
-        if order is not None:
-            raise NotImplementedError("`order` kwarg is not supported upgrade "
-                                      "to Numpy 1.8 or greater for support.")
-        return np.multiply(fill_value, np.ones(shape, dtype=dtype),
-                           dtype=dtype)
-
-
 # Taken from scikit-learn:
 # https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/utils/fixes.py#L84
 try:
@@ -54,79 +34,12 @@ except TypeError:
                                                     np.ma.core._DomainSafeDivide(),
                                                     0, 1)
 
-# functions copied from numpy
-try:
-    from numpy import broadcast_to, nanprod, nancumsum, nancumprod
-except ImportError:  # pragma: no cover
-    # these functions should arrive in numpy v1.10 to v1.12.  Until then,
-    # they are duplicated here
+if LooseVersion(np.__version__) < '1.12.0':
+    # These functions were added in numpy 1.12.0. For previous versions they
+    # are duplicated here
 
     # See https://github.com/numpy/numpy/blob/master/LICENSE.txt
     # or NUMPY_LICENSE.txt within this directory
-    def _maybe_view_as_subclass(original_array, new_array):
-        if type(original_array) is not type(new_array):
-            # if input was an ndarray subclass and subclasses were OK,
-            # then view the result as that subclass.
-            new_array = new_array.view(type=type(original_array))
-            # Since we have done something akin to a view from original_array, we
-            # should let the subclass finalize (if it has it implemented, i.e., is
-            # not None).
-            if new_array.__array_finalize__:
-                new_array.__array_finalize__(original_array)
-        return new_array
-
-    def _broadcast_to(array, shape, subok, readonly):
-        shape = tuple(shape) if np.iterable(shape) else (shape,)
-        array = np.array(array, copy=False, subok=subok)
-        if not shape and array.shape:
-            raise ValueError('cannot broadcast a non-scalar to a scalar array')
-        if builtins.any(size < 0 for size in shape):
-            raise ValueError('all elements of broadcast shape must be non-'
-                             'negative')
-        broadcast = np.nditer(
-            (array,), flags=['multi_index', 'zerosize_ok', 'refs_ok'],
-            op_flags=['readonly'], itershape=shape, order='C').itviews[0]
-        result = _maybe_view_as_subclass(array, broadcast)
-        if not readonly and array.flags.writeable:
-            result.flags.writeable = True
-        return result
-
-    def broadcast_to(array, shape, subok=False):
-        """Broadcast an array to a new shape.
-
-        Parameters
-        ----------
-        array : array_like
-            The array to broadcast.
-        shape : tuple
-            The shape of the desired array.
-        subok : bool, optional
-            If True, then sub-classes will be passed-through, otherwise
-            the returned array will be forced to be a base-class array (default).
-
-        Returns
-        -------
-        broadcast : array
-            A readonly view on the original array with the given shape. It is
-            typically not contiguous. Furthermore, more than one element of a
-            broadcasted array may refer to a single memory location.
-
-        Raises
-        ------
-        ValueError
-            If the array is not compatible with the new shape according to NumPy's
-            broadcasting rules.
-
-        Examples
-        --------
-        >>> x = np.array([1, 2, 3])
-        >>> np.broadcast_to(x, (3, 3))  # doctest: +SKIP
-        array([[1, 2, 3],
-               [1, 2, 3],
-               [1, 2, 3]])
-        """
-        return _broadcast_to(array, shape, subok=subok, readonly=True)
-
     def _replace_nan(a, val):
         """
         If `a` is of inexact type, make a copy of `a`, replace NaNs with
@@ -168,75 +81,6 @@ except ImportError:  # pragma: no cover
         np.copyto(a, val, where=mask)
         return a, mask
 
-    def nanprod(a, axis=None, dtype=None, out=None, keepdims=0):
-        """
-        Return the product of array elements over a given axis treating Not a
-        Numbers (NaNs) as zero.
-
-        One is returned for slices that are all-NaN or empty.
-
-        .. versionadded:: 1.10.0
-
-        Parameters
-        ----------
-        a : array_like
-            Array containing numbers whose sum is desired. If `a` is not an
-            array, a conversion is attempted.
-        axis : int, optional
-            Axis along which the product is computed. The default is to compute
-            the product of the flattened array.
-        dtype : data-type, optional
-            The type of the returned array and of the accumulator in which the
-            elements are summed.  By default, the dtype of `a` is used.  An
-            exception is when `a` has an integer type with less precision than
-            the platform (u)intp. In that case, the default will be either
-            (u)int32 or (u)int64 depending on whether the platform is 32 or 64
-            bits. For inexact inputs, dtype must be inexact.
-        out : ndarray, optional
-            Alternate output array in which to place the result.  The default
-            is ``None``. If provided, it must have the same shape as the
-            expected output, but the type will be cast if necessary.  See
-            `doc.ufuncs` for details. The casting of NaN to integer can yield
-            unexpected results.
-        keepdims : bool, optional
-            If True, the axes which are reduced are left in the result as
-            dimensions with size one. With this option, the result will
-            broadcast correctly against the original `arr`.
-
-        Returns
-        -------
-        y : ndarray or numpy scalar
-
-        See Also
-        --------
-        :func:`numpy.prod` : Product across array propagating NaNs.
-        isnan : Show which elements are NaN.
-
-        Notes
-        -----
-        Numpy integer arithmetic is modular. If the size of a product exceeds
-        the size of an integer accumulator, its value will wrap around and the
-        result will be incorrect. Specifying ``dtype=double`` can alleviate
-        that problem.
-
-        Examples
-        --------
-        >>> np.nanprod(1)
-        1
-        >>> np.nanprod([1])
-        1
-        >>> np.nanprod([1, np.nan])
-        1.0
-        >>> a = np.array([[1, 2], [3, np.nan]])
-        >>> np.nanprod(a)
-        6.0
-        >>> np.nanprod(a, axis=0)
-        array([ 3.,  2.])
-
-        """
-        a, mask = _replace_nan(a, 1)
-        return np.prod(a, axis=axis, dtype=dtype, out=out, keepdims=keepdims)
-
     def nancumsum(a, axis=None, dtype=None, out=None):
         """
         Return the cumulative sum of array elements over a given axis treating Not a
diff --git a/dask/array/reductions.py b/dask/array/reductions.py
index 8598a3519..29e6a9782 100644
--- a/dask/array/reductions.py
+++ b/dask/array/reductions.py
@@ -12,6 +12,7 @@ from toolz import compose, partition_all, get, accumulate, pluck
 
 from . import chunk
 from .core import _concatenate2, Array, atop, lol_tuples, handle_out
+from .creation import arange
 from .ufunc import sqrt
 from .wrap import zeros, ones
 from .numpy_compat import ma_divide, divide as np_divide
@@ -722,3 +723,73 @@ def validate_axis(ndim, axis):
         return axis + ndim
     else:
         return axis
+
+
+def topk(a, k, axis=-1, split_every=None):
+    """Extract the k largest elements from a on the given axis,
+    and return them sorted from largest to smallest.
+    If k is negative, extract the -k smallest elements instead,
+    and return them sorted from smallest to largest.
+
+    This assumes that ``k`` is small.  All results will be returned in a single
+    chunk along the given axis.
+
+    Examples
+    --------
+    >>> import dask.array as da
+    >>> x = np.array([5, 1, 3, 6])
+    >>> d = da.from_array(x, chunks=2)
+    >>> d.topk(2).compute()
+    array([6, 5])
+    >>> d.topk(-2).compute()
+    array([1, 3])
+    """
+    if isinstance(a, int) and isinstance(k, Array):
+        warnings.warn("DeprecationWarning: topk(k, a) has been replaced with topk(a, k)")
+        a, k = k, a
+
+    axis = validate_axis(a.ndim, axis)
+
+    kernel = partial(chunk.topk, k=k)
+    res = reduction(a, kernel, kernel, axis=axis, keepdims=True,
+                    dtype=a.dtype, split_every=split_every)
+    # reduction(keepdims=True) sets shape[axis] to 1. Fix it.
+    chunks = list(res.chunks)
+    chunks[axis] = (abs(k), )
+    res = Array(res.dask, res.name, chunks, res.dtype)
+
+    # Sort result internally
+    return res.map_blocks(chunk.topk_postprocess, k=k, axis=axis, dtype=a.dtype)
+
+
+def argtopk(a, k, axis=-1, split_every=None):
+    """Extract the indices of the k largest elements from a on the given axis,
+    and return them sorted from largest to smallest.
+    If k is negative, extract the indices of the -k smallest elements instead,
+    and return them sorted from smallest to largest.
+
+    This assumes that ``k`` is small.  All results will be returned in a single
+    chunk along the given axis.
+
+    Examples
+    --------
+    >>> import dask.array as da
+    >>> x = np.array([5, 1, 3, 6])
+    >>> d = da.from_array(x, chunks=2)
+    >>> d.argtopk(2).compute()
+    array([3, 0])
+    >>> d.argtopk(-2).compute()
+    array([1, 2])
+    """
+    axis = validate_axis(a.ndim, axis)
+
+    # Convert a to a recarray that contains its index
+    idx = arange(a.shape[axis], chunks=a.chunks[axis], dtype=np.int64)
+    idx = idx[tuple(slice(None) if i == axis else np.newaxis for i in range(a.ndim))]
+    a_rec = a.map_blocks(chunk.argtopk_preprocess, idx,
+                         dtype=[('a', a.dtype), ('idx', idx.dtype)])
+
+    res = topk(a_rec, k, axis=axis, split_every=split_every)
+
+    # Discard values
+    return res['idx']
diff --git a/dask/array/routines.py b/dask/array/routines.py
index 819daba9e..72e88433a 100644
--- a/dask/array/routines.py
+++ b/dask/array/routines.py
@@ -6,7 +6,6 @@ from collections import Iterable
 from distutils.version import LooseVersion
 from functools import wraps, partial
 from numbers import Integral
-from operator import getitem
 
 import numpy as np
 from toolz import concat, sliding_window, interleave
@@ -14,7 +13,7 @@ from toolz import concat, sliding_window, interleave
 from .. import sharedict
 from ..core import flatten
 from ..base import tokenize
-from . import numpy_compat, chunk
+from . import chunk
 from .creation import arange
 from .utils import safe_wraps
 from .wrap import ones
@@ -881,38 +880,6 @@ def squeeze(a, axis=None):
     return a[sl]
 
 
-def topk(k, x):
-    """ The top k elements of an array
-
-    Returns the k greatest elements of the array in sorted order.  Only works
-    on arrays of a single dimension.
-
-    This assumes that ``k`` is small.  All results will be returned in a single
-    chunk.
-
-    Examples
-    --------
-
-    >>> x = np.array([5, 1, 3, 6])
-    >>> d = from_array(x, chunks=2)
-    >>> d.topk(2).compute()
-    array([6, 5])
-    """
-    if x.ndim != 1:
-        raise ValueError("Topk only works on arrays of one dimension")
-
-    token = tokenize(k, x)
-    name = 'chunk.topk-' + token
-    dsk = {(name, i): (chunk.topk, k, key)
-           for i, key in enumerate(x.__dask_keys__())}
-    name2 = 'topk-' + token
-    dsk[(name2, 0)] = (getitem, (np.sort, (np.concatenate, list(dsk))),
-                       slice(-1, -k - 1, -1))
-    chunks = ((k,),)
-
-    return Array(sharedict.merge((name2, dsk), x.dask), name2, chunks, dtype=x.dtype)
-
-
 @wraps(np.compress)
 def compress(condition, a, axis=None):
     if axis is None:
@@ -995,9 +962,9 @@ def notnull(values):
     return ~isnull(values)
 
 
-@wraps(numpy_compat.isclose)
+@wraps(np.isclose)
 def isclose(arr1, arr2, rtol=1e-5, atol=1e-8, equal_nan=False):
-    func = partial(numpy_compat.isclose, rtol=rtol, atol=atol, equal_nan=equal_nan)
+    func = partial(np.isclose, rtol=rtol, atol=atol, equal_nan=equal_nan)
     return elemwise(func, arr1, arr2, dtype='bool')
 
 
diff --git a/dask/array/utils.py b/dask/array/utils.py
index 7e35865f7..3ac9b43e2 100644
--- a/dask/array/utils.py
+++ b/dask/array/utils.py
@@ -1,6 +1,5 @@
 from __future__ import absolute_import, division, print_function
 
-from distutils.version import LooseVersion
 import difflib
 import functools
 import math
@@ -14,23 +13,9 @@ from ..local import get_sync
 from ..sharedict import ShareDict
 
 
-if LooseVersion(np.__version__) >= '1.10.0':
-    _allclose = np.allclose
-else:
-    def _allclose(a, b, **kwargs):
-        if kwargs.pop('equal_nan', False):
-            a_nans = np.isnan(a)
-            b_nans = np.isnan(b)
-            if not (a_nans == b_nans).all():
-                return False
-            a = a[~a_nans]
-            b = b[~b_nans]
-        return np.allclose(a, b, **kwargs)
-
-
 def allclose(a, b, equal_nan=False, **kwargs):
     if getattr(a, 'dtype', None) != 'O':
-        return _allclose(a, b, equal_nan=equal_nan, **kwargs)
+        return np.allclose(a, b, equal_nan=equal_nan, **kwargs)
     if equal_nan:
         return (a.shape == b.shape and
                 all(np.isnan(b) if np.isnan(a) else a == b
diff --git a/dask/array/wrap.py b/dask/array/wrap.py
index accee3d64..1306c781d 100644
--- a/dask/array/wrap.py
+++ b/dask/array/wrap.py
@@ -12,7 +12,6 @@ except ImportError:
 
 from ..base import tokenize
 from .core import Array, normalize_chunks
-from .numpy_compat import full
 
 
 def wrap_func_shape_as_first_arg(func, *args, **kwargs):
@@ -68,4 +67,4 @@ w = wrap(wrap_func_shape_as_first_arg)
 ones = w(np.ones, dtype='f8')
 zeros = w(np.zeros, dtype='f8')
 empty = w(np.empty, dtype='f8')
-full = w(full)
+full = w(np.full)
diff --git a/docs/source/array-api.rst b/docs/source/array-api.rst
index b253f28ec..f6fa498c4 100644
--- a/docs/source/array-api.rst
+++ b/docs/source/array-api.rst
@@ -22,6 +22,7 @@ Top level user functions:
    arctanh
    argmax
    argmin
+   argtopk
    argwhere
    around
    array
@@ -341,7 +342,6 @@ Other functions
 .. autofunction:: from_array
 .. autofunction:: from_delayed
 .. autofunction:: store
-.. autofunction:: topk
 .. autofunction:: coarsen
 .. autofunction:: stack
 .. autofunction:: concatenate
@@ -362,6 +362,7 @@ Other functions
 .. autofunction:: arctanh
 .. autofunction:: argmax
 .. autofunction:: argmin
+.. autofunction:: argtopk
 .. autofunction:: argwhere
 .. autofunction:: around
 .. autofunction:: array
diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst
index d55f3a042..301db78af 100644
--- a/docs/source/changelog.rst
+++ b/docs/source/changelog.rst
@@ -12,6 +12,12 @@ Array
 - Add ``piecewise`` for Dask Arrays (:pr:`3350`) `John A Kirkham`_
 - Fix handling of ``nan`` in ``broadcast_shapes`` (:pr:`3356`) `John A Kirkham`_
 - Add ``isin`` for dask arrays (:pr:`3363`). `Stephan Hoyer`_
+- Overhauled ``topk`` for Dask Arrays: faster algorithm, particularly for large k's; added support
+  for multiple axes, recursive aggregation, and an option to pick the bottom k elements instead.
+  (:pr:`3395`) `Guido Imperiale`_
+- The ``topk`` API has changed from topk(k, array) to the more conventional topk(array, k).
+  The legacy API still works but is now deprecated. (:pr:`2965`) `Guido Imperiale`_
+- New function ``argtopk`` for Dask Arrays (:pr:`3396`) `Guido Imperiale`_
 
 DataFrame
 +++++++++
@@ -1027,6 +1033,7 @@ Other
 - There is also a gitter chat room and a stackoverflow tag
 
 
+.. _`Guido Imperiale`: https://github.com/crusaderky
 .. _`John A Kirkham`: https://github.com/jakirkham
 .. _`Matthew Rocklin`: https://github.com/mrocklin
 .. _`Jim Crist`: https://github.com/jcrist
diff --git a/setup.py b/setup.py
index 75ccaa59f..8daeba180 100755
--- a/setup.py
+++ b/setup.py
@@ -8,9 +8,9 @@ import versioneer
 # NOTE: These are tested in `continuous_integration/travis/test_imports.sh` If
 # you modify these, make sure to change the corresponding line there.
 extras_require = {
-  'array': ['numpy >= 1.10.4', 'toolz >= 0.7.3'],
+  'array': ['numpy >= 1.11.0', 'toolz >= 0.7.3'],
   'bag': ['cloudpickle >= 0.2.1', 'toolz >= 0.7.3', 'partd >= 0.3.8'],
-  'dataframe': ['numpy >= 1.10.4', 'pandas >= 0.19.0', 'toolz >= 0.7.3',
+  'dataframe': ['numpy >= 1.11.0', 'pandas >= 0.19.0', 'toolz >= 0.7.3',
                 'partd >= 0.3.8', 'cloudpickle >= 0.2.1'],
   'distributed': ['distributed >= 1.21'],
   'delayed': ['toolz >= 0.7.3'],


Test patch:
diff --git a/dask/array/tests/test_array_core.py b/dask/array/tests/test_array_core.py
index adf24fe13..211a2b849 100644
--- a/dask/array/tests/test_array_core.py
+++ b/dask/array/tests/test_array_core.py
@@ -286,8 +286,6 @@ def test_stack_promote_type():
     assert_eq(res, np.stack([i, f]))
 
 
-@pytest.mark.skipif(LooseVersion(np.__version__) < '1.10.0',
-                    reason="NumPy doesn't yet support stack")
 def test_stack_rechunk():
     x = da.random.random(10, chunks=5)
     y = da.random.random(10, chunks=4)
diff --git a/dask/array/tests/test_linalg.py b/dask/array/tests/test_linalg.py
index eafd86144..18d9e5473 100644
--- a/dask/array/tests/test_linalg.py
+++ b/dask/array/tests/test_linalg.py
@@ -500,10 +500,6 @@ def test_norm_1dim(shape, chunks, axis, norm, keepdims):
 
     a_r = np.linalg.norm(a, ord=norm, axis=axis, keepdims=keepdims)
     d_r = da.linalg.norm(d, ord=norm, axis=axis, keepdims=keepdims)
-
-    # Fix a type mismatch on NumPy 1.10.
-    a_r = a_r.astype(float)
-
     assert_eq(a_r, d_r)
 
 
diff --git a/dask/array/tests/test_reductions.py b/dask/array/tests/test_reductions.py
index 2dda9e1d5..ed6a2e6a0 100644
--- a/dask/array/tests/test_reductions.py
+++ b/dask/array/tests/test_reductions.py
@@ -1,21 +1,13 @@
 from __future__ import absolute_import, division, print_function
 
 import pytest
-pytest.importorskip('numpy')
+np = pytest.importorskip('numpy')
 
 import dask.array as da
 from dask.array.utils import assert_eq as _assert_eq, same_keys
 from dask.core import get_deps
 from dask.context import set_options
 
-import numpy as np
-# temporary until numpy functions migrated
-try:
-    from numpy import nanprod
-except ImportError:  # pragma: no cover
-    import dask.array.numpy_compat as npcompat
-    nanprod = npcompat.nanprod
-
 
 def assert_eq(a, b):
     _assert_eq(a, b, equal_nan=True)
@@ -56,7 +48,7 @@ def test_reductions_1D(dtype):
     reduction_1d_test(da.all, a, np.all, x, False)
 
     reduction_1d_test(da.nansum, a, np.nansum, x)
-    reduction_1d_test(da.nanprod, a, nanprod, x)
+    reduction_1d_test(da.nanprod, a, np.nanprod, x)
     reduction_1d_test(da.nanmean, a, np.mean, x)
     reduction_1d_test(da.nanvar, a, np.var, x)
     reduction_1d_test(da.nanstd, a, np.std, x)
@@ -127,7 +119,7 @@ def test_reductions_2D(dtype):
     reduction_2d_test(da.all, a, np.all, x, False)
 
     reduction_2d_test(da.nansum, a, np.nansum, x)
-    reduction_2d_test(da.nanprod, a, nanprod, x)
+    reduction_2d_test(da.nanprod, a, np.nanprod, x)
     reduction_2d_test(da.nanmean, a, np.mean, x)
     reduction_2d_test(da.nanvar, a, np.nanvar, x, False)  # Difference in dtype algo
     reduction_2d_test(da.nanstd, a, np.nanstd, x, False)  # Difference in dtype algo
@@ -207,7 +199,7 @@ def test_reductions_2D_nans():
     reduction_2d_test(da.all, a, np.all, x, False, False)
 
     reduction_2d_test(da.nansum, a, np.nansum, x, False, False)
-    reduction_2d_test(da.nanprod, a, nanprod, x, False, False)
+    reduction_2d_test(da.nanprod, a, np.nanprod, x, False, False)
     reduction_2d_test(da.nanmean, a, np.nanmean, x, False, False)
     with pytest.warns(None):  # division by 0 warning
         reduction_2d_test(da.nanvar, a, np.nanvar, x, False, False)
@@ -287,7 +279,7 @@ def test_nan():
     assert_eq(np.nanstd(x, axis=0), da.nanstd(d, axis=0))
     assert_eq(np.nanargmin(x, axis=0), da.nanargmin(d, axis=0))
     assert_eq(np.nanargmax(x, axis=0), da.nanargmax(d, axis=0))
-    assert_eq(nanprod(x), da.nanprod(d))
+    assert_eq(np.nanprod(x), da.nanprod(d))
 
 
 @pytest.mark.skipif(np.__version__ < '1.13.0', reason='nanmax/nanmin for object dtype')
@@ -422,3 +414,65 @@ def test_array_cumreduction_out(func):
     x = da.ones((10, 10), chunks=(4, 4))
     func(x, axis=0, out=x)
     assert_eq(x, func(np.ones((10, 10)), axis=0))
+
+
+@pytest.mark.parametrize('npfunc,daskfunc', [
+    (np.sort, da.topk),
+    (np.argsort, da.argtopk),
+])
+@pytest.mark.parametrize('split_every', [None, 2])
+def test_topk_argtopk1(npfunc, daskfunc, split_every):
+    # Test data
+    k = 5
+    a = da.random.random(1000, chunks=250)
+    b = da.random.random((10, 20, 30), chunks=(4, 8, 8))
+    c = da.from_array(np.array([(1, 'Hello'), (2, 'World')], dtype=[('foo', int), ('bar', '<U5')]),
+                      chunks=1)
+
+    # 1-dimensional arrays
+    # top 5 elements, sorted descending
+    assert_eq(npfunc(a)[-k:][::-1],
+              daskfunc(a, k, split_every=split_every))
+    # bottom 5 elements, sorted ascending
+    assert_eq(npfunc(a)[:k],
+              daskfunc(a, -k, split_every=split_every))
+
+    # n-dimensional arrays
+    # also testing when k > chunk
+    # top 5 elements, sorted descending
+    assert_eq(npfunc(b, axis=0)[-k:, :, :][::-1, :, :],
+              daskfunc(b, k, axis=0, split_every=split_every))
+    assert_eq(npfunc(b, axis=1)[:, -k:, :][:, ::-1, :],
+              daskfunc(b, k, axis=1, split_every=split_every))
+    assert_eq(npfunc(b, axis=-1)[:, :, -k:][:, :, ::-1],
+              daskfunc(b, k, axis=-1, split_every=split_every))
+    with pytest.raises(ValueError):
+        daskfunc(b, k, axis=3, split_every=split_every)
+
+    # bottom 5 elements, sorted ascending
+    assert_eq(npfunc(b, axis=0)[:k, :, :],
+              daskfunc(b, -k, axis=0, split_every=split_every))
+    assert_eq(npfunc(b, axis=1)[:, :k, :],
+              daskfunc(b, -k, axis=1, split_every=split_every))
+    assert_eq(npfunc(b, axis=-1)[:, :, :k],
+              daskfunc(b, -k, axis=-1, split_every=split_every))
+    with pytest.raises(ValueError):
+        daskfunc(b, -k, axis=3, split_every=split_every)
+
+    # structured arrays
+    assert_eq(npfunc(c, axis=0)[-1:][::-1],
+              daskfunc(c, 1, split_every=split_every))
+    assert_eq(npfunc(c, axis=0)[:1],
+              daskfunc(c, -1, split_every=split_every))
+
+
+def test_topk_argtopk2():
+    a = da.random.random((10, 20, 30), chunks=(4, 8, 8))
+
+    # Support for deprecated API for topk
+    with pytest.warns(UserWarning):
+        assert_eq(da.topk(a, 5), da.topk(5, a))
+
+    # As Array methods
+    assert_eq(a.topk(5, axis=1, split_every=2), da.topk(a, 5, axis=1, split_every=2))
+    assert_eq(a.argtopk(5, axis=1, split_every=2), da.argtopk(a, 5, axis=1, split_every=2))
diff --git a/dask/array/tests/test_routines.py b/dask/array/tests/test_routines.py
index 882941d35..7b8443644 100644
--- a/dask/array/tests/test_routines.py
+++ b/dask/array/tests/test_routines.py
@@ -11,6 +11,7 @@ np = pytest.importorskip('numpy')
 import dask.array as da
 from dask.utils import ignoring
 from dask.array.utils import assert_eq, same_keys
+from dask.array.einsumfuncs import einsum_can_optimize
 
 
 def test_array():
@@ -418,27 +419,6 @@ def test_ediff1d(shape, to_end, to_begin):
     assert_eq(da.ediff1d(a, to_end, to_begin), np.ediff1d(x, to_end, to_begin))
 
 
-def test_topk():
-    x = np.array([5, 2, 1, 6])
-    d = da.from_array(x, chunks=2)
-
-    e = da.topk(2, d)
-
-    assert e.chunks == ((2,),)
-    assert_eq(e, np.sort(x)[-1:-3:-1])
-    assert same_keys(da.topk(2, d), e)
-
-
-def test_topk_k_bigger_than_chunk():
-    x = np.array([5, 2, 1, 6])
-    d = da.from_array(x, chunks=2)
-
-    e = da.topk(3, d)
-
-    assert e.chunks == ((3,),)
-    assert_eq(e, np.array([6, 5, 2]))
-
-
 def test_bincount():
     x = np.array([2, 1, 5, 2, 1])
     d = da.from_array(x, chunks=2)
@@ -469,8 +449,6 @@ def test_bincount_raises_informative_error_on_missing_minlength_kwarg():
         assert False
 
 
-@pytest.mark.skipif(LooseVersion(np.__version__) < '1.10.0',
-                    reason="NumPy doesn't yet support nd digitize")
 def test_digitize():
     x = np.array([2, 4, 5, 6, 1])
     bins = np.array([1, 2, 3, 4, 5])
@@ -1356,6 +1334,8 @@ def test_einsum(einsum_signature):
               da.einsum(einsum_signature, *da_inputs))
 
 
+@pytest.mark.skipif(not einsum_can_optimize,
+                    reason="np.einsum(optimize) unavailable")
 @pytest.mark.parametrize('optimize_opts', [
     (True, False),
     ('greedy', False),
diff --git a/dask/array/tests/test_sparse.py b/dask/array/tests/test_sparse.py
index dc7519f5a..028bcfef4 100644
--- a/dask/array/tests/test_sparse.py
+++ b/dask/array/tests/test_sparse.py
@@ -9,7 +9,7 @@ from dask.array.utils import assert_eq
 
 sparse = pytest.importorskip('sparse')
 
-if LooseVersion(np.__version__) < '1.11.0':
+if LooseVersion(np.__version__) < '1.11.2':
     pytestmark = pytest.mark.skip
 
 
diff --git a/dask/dataframe/tests/test_dataframe.py b/dask/dataframe/tests/test_dataframe.py
index 6c447ecd2..35e290f7c 100644
--- a/dask/dataframe/tests/test_dataframe.py
+++ b/dask/dataframe/tests/test_dataframe.py
@@ -1,6 +1,7 @@
 import sys
-from operator import add
+from distutils.version import LooseVersion
 from itertools import product
+from operator import add
 
 import pandas as pd
 import pandas.util.testing as tm
@@ -331,24 +332,25 @@ def test_cumulative():
     assert_eq(ddf.cummin(axis=1), df.cummin(axis=1))
     assert_eq(ddf.cummax(axis=1), df.cummax(axis=1))
 
-    # testing out parameter
-    np.cumsum(ddf, out=ddf_out)
-    assert_eq(ddf_out, df.cumsum())
-    np.cumprod(ddf, out=ddf_out)
-    assert_eq(ddf_out, df.cumprod())
-    ddf.cummin(out=ddf_out)
-    assert_eq(ddf_out, df.cummin())
-    ddf.cummax(out=ddf_out)
-    assert_eq(ddf_out, df.cummax())
-
-    np.cumsum(ddf, out=ddf_out, axis=1)
-    assert_eq(ddf_out, df.cumsum(axis=1))
-    np.cumprod(ddf, out=ddf_out, axis=1)
-    assert_eq(ddf_out, df.cumprod(axis=1))
-    ddf.cummin(out=ddf_out, axis=1)
-    assert_eq(ddf_out, df.cummin(axis=1))
-    ddf.cummax(out=ddf_out, axis=1)
-    assert_eq(ddf_out, df.cummax(axis=1))
+    # testing out parameter if out parameter supported
+    if LooseVersion(np.__version__) >= '1.13.0':
+        np.cumsum(ddf, out=ddf_out)
+        assert_eq(ddf_out, df.cumsum())
+        np.cumprod(ddf, out=ddf_out)
+        assert_eq(ddf_out, df.cumprod())
+        ddf.cummin(out=ddf_out)
+        assert_eq(ddf_out, df.cummin())
+        ddf.cummax(out=ddf_out)
+        assert_eq(ddf_out, df.cummax())
+
+        np.cumsum(ddf, out=ddf_out, axis=1)
+        assert_eq(ddf_out, df.cumsum(axis=1))
+        np.cumprod(ddf, out=ddf_out, axis=1)
+        assert_eq(ddf_out, df.cumprod(axis=1))
+        ddf.cummin(out=ddf_out, axis=1)
+        assert_eq(ddf_out, df.cummin(axis=1))
+        ddf.cummax(out=ddf_out, axis=1)
+        assert_eq(ddf_out, df.cummax(axis=1))
 
     assert_eq(ddf.a.cumsum(), df.a.cumsum())
     assert_eq(ddf.a.cumprod(), df.a.cumprod())