ENH: missingdata: Add skipna=, keepdims= parameters to methods

Also fix some memory leaks, improve some type resolution code. The
methods still have some issues with array subtypes that needs working
through.

7 files changed, 148 insertions, 229 deletions

diff --git a/doc/release/2.0.0-notes.rst b/doc/release/2.0.0-notes.rst
index d1ebd1b75..4150b1ee7 100644
--- a/doc/release/2.0.0-notes.rst
+++ b/doc/release/2.0.0-notes.rst
@@ -25,7 +25,7 @@ list of things that do and do not work with NA values:
 What works with NA:
     * Basic indexing and slicing, as well as full boolean mask indexing.
     * All element-wise ufuncs.
-    * UFunc.reduce methods, with a new skipna parameter.
+    * All UFunc.reduce methods, with a new skipna parameter.
     * The nditer object.
     * Array methods:
        + ndarray.clip, ndarray.min, ndarray.max, ndarray.sum, ndarray.prod,
@@ -40,8 +40,6 @@ What doesn't work with NA:
       mechanism instead of the newer nditer.
     * Struct dtypes, which will have corresponding struct masks with
       one mask value per primitive field of the struct dtype.
-    * UFunc.reduce of multi-dimensional arrays, with skipna=True and a ufunc
-      that doesn't have an identity.
     * UFunc.accumulate, UFunc.reduceat.
     * Ufuncs calls with both NA masks and a where= mask at the same time.
     * np.logical_and, np.logical_or, np.all, and np.any don't satisfy the
@@ -59,7 +57,7 @@ Differences with R:
     * np.isna(nan) is False, but R's is.na(nan) is TRUE. This is because
       NumPy's NA is treated independently of the underlying data type.
     * Boolean indexing, where the result is compressed to just
-      the elements with true in the mask, raises if the booelan mask
+      the elements with true in the mask, raises if the boolean mask
       has an NA value in it. This is because that value could be either
       True or False, meaning the count of the output array is actually
       NA. R treats this case in a manner inconsistent with the NA model,
diff --git a/numpy/core/fromnumeric.py b/numpy/core/fromnumeric.py
index f3930eaa9..f374951c9 100644
--- a/numpy/core/fromnumeric.py
+++ b/numpy/core/fromnumeric.py
@@ -16,6 +16,7 @@ import multiarray as mu
 import umath as um
 import numerictypes as nt
 from numeric import asarray, array, asanyarray, concatenate
+import _methods
 _dt_ = nt.sctype2char
 
 import types
@@ -1482,12 +1483,12 @@ def sum(a, axis=None, dtype=None, out=None, skipna=False, keepdims=False):
         try:
             sum = a.sum
         except AttributeError:
-            return um.add.reduce(a, axis=axis, dtype=dtype,
+            return _methods._sum(a, axis=axis, dtype=dtype,
                                 out=out, skipna=skipna, keepdims=keepdims)
         # NOTE: Dropping the skipna and keepdims parameters here...
         return sum(axis=axis, dtype=dtype, out=out)
     else:
-        return um.add.reduce(a, axis=axis, dtype=dtype,
+        return _methods._sum(a, axis=axis, dtype=dtype,
                             out=out, skipna=skipna, keepdims=keepdims)
 
 def product (a, axis=None, dtype=None, out=None, skipna=False, keepdims=False):
@@ -1603,7 +1604,8 @@ def any(a, axis=None, out=None, skipna=False, keepdims=False):
     (191614240, 191614240)
 
     """
-    return um.logical_or.reduce(a, axis=axis, out=out, skipna=skipna, keepdims=keepdims)
+    return _methods._any(a, axis=axis, out=out,
+                            skipna=skipna, keepdims=keepdims)
 
 def all(a, axis=None, out=None, skipna=False, keepdims=False):
     """
@@ -1674,7 +1676,8 @@ def all(a, axis=None, out=None, skipna=False, keepdims=False):
     (28293632, 28293632, array([ True], dtype=bool))
 
     """
-    return um.logical_and.reduce(a, axis=axis, out=out, skipna=skipna, keepdims=keepdims)
+    return _methods._all(a, axis=axis, out=out,
+                            skipna=skipna, keepdims=keepdims)
 
 def cumsum (a, axis=None, dtype=None, out=None):
     """
@@ -1873,12 +1876,12 @@ def amax(a, axis=None, out=None, skipna=False, keepdims=False):
         try:
             amax = a.max
         except AttributeError:
-            return um.maximum.reduce(a, axis=axis,
+            return _methods._amax(a, axis=axis,
                                 out=out, skipna=skipna, keepdims=keepdims)
         # NOTE: Dropping the skipna and keepdims parameters
         return amax(axis=axis, out=out)
     else:
-        return um.maximum.reduce(a, axis=axis,
+        return _methods._amax(a, axis=axis,
                             out=out, skipna=skipna, keepdims=keepdims)
 
 def amin(a, axis=None, out=None, skipna=False, keepdims=False):
@@ -1947,12 +1950,12 @@ def amin(a, axis=None, out=None, skipna=False, keepdims=False):
         try:
             amin = a.min
         except AttributeError:
-            return um.minimum.reduce(a, axis=axis,
+            return _methods._amin(a, axis=axis,
                                 out=out, skipna=skipna, keepdims=keepdims)
         # NOTE: Dropping the skipna and keepdims parameters
         return amin(axis=axis, out=out)
     else:
-        return um.minimum.reduce(a, axis=axis,
+        return _methods._amin(a, axis=axis,
                             out=out, skipna=skipna, keepdims=keepdims)
 
 def alen(a):
@@ -2080,11 +2083,11 @@ def prod(a, axis=None, dtype=None, out=None, skipna=False, keepdims=False):
         try:
             prod = a.prod
         except AttributeError:
-            return um.multiply.reduce(a, axis=axis, dtype=dtype,
+            return _methods._prod(a, axis=axis, dtype=dtype,
                                 out=out, skipna=skipna, keepdims=keepdims)
         return prod(axis=axis, dtype=dtype, out=out)
     else:
-        return um.multiply.reduce(a, axis=axis, dtype=dtype,
+        return _methods._prod(a, axis=axis, dtype=dtype,
                             out=out, skipna=skipna, keepdims=keepdims)
 
 def cumprod(a, axis=None, dtype=None, out=None):
@@ -2459,22 +2462,8 @@ def mean(a, axis=None, dtype=None, out=None, skipna=False, keepdims=False):
         except AttributeError:
             pass
 
-    arr = asarray(a)
-
-    # Upgrade bool, unsigned int, and int to float64
-    if dtype is None and arr.dtype.kind in ['b','u','i']:
-        ret = um.add.reduce(arr, axis=axis, dtype='f8',
+    return _methods._mean(a, axis=axis, dtype=dtype,
                             out=out, skipna=skipna, keepdims=keepdims)
-    else:
-        ret = um.add.reduce(arr, axis=axis, dtype=dtype,
-                            out=out, skipna=skipna, keepdims=keepdims)
-    rcount = mu.count_reduce_items(arr, axis=axis,
-                            skipna=skipna, keepdims=keepdims)
-    if isinstance(ret, mu.ndarray):
-        um.true_divide(ret, rcount, out=ret, casting='unsafe')
-    else:
-        ret = ret / float(rcount)
-    return ret
 
 
 def std(a, axis=None, dtype=None, out=None, ddof=0,
@@ -2579,16 +2568,9 @@ def std(a, axis=None, dtype=None, out=None, ddof=0,
         except AttributeError:
             pass
 
-    ret = var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
+    return _methods._std(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
                                 skipna=skipna, keepdims=keepdims)
 
-    if isinstance(ret, mu.ndarray):
-        um.sqrt(ret, out=ret)
-    else:
-        ret = um.sqrt(ret)
-
-    return ret
-
 def var(a, axis=None, dtype=None, out=None, ddof=0,
                             skipna=False, keepdims=False):
     """
@@ -2692,43 +2674,6 @@ def var(a, axis=None, dtype=None, out=None, ddof=0,
         except AttributeError:
             pass
 
-    arr = asarray(a)
-
-    # First compute the mean, saving 'rcount' for reuse later
-    if dtype is None and arr.dtype.kind in ['b','u','i']:
-        arrmean = um.add.reduce(arr, axis=axis, dtype='f8',
-                            skipna=skipna, keepdims=True)
-    else:
-        arrmean = um.add.reduce(arr, axis=axis, dtype=dtype,
-                            skipna=skipna, keepdims=True)
-    rcount = mu.count_reduce_items(arr, axis=axis,
-                            skipna=skipna, keepdims=True)
-    if isinstance(arrmean, mu.ndarray):
-        um.true_divide(arrmean, rcount, out=arrmean, casting='unsafe')
-    else:
-        arrmean = arrmean / float(rcount)
-
-    # arr - arrmean
-    x = arr - arrmean
-
-    # (arr - arrmean) ** 2
-    if arr.dtype.kind == 'c':
-        um.multiply(x, um.conjugate(x), out=x)
-        x = x.real
-    else:
-        um.multiply(x, x, out=x)
-
-    # add.reduce((arr - arrmean) ** 2, axis)
-    ret = um.add.reduce(x, axis=axis, dtype=dtype, out=out,
+    return _methods._var(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
                                 skipna=skipna, keepdims=keepdims)
 
-    # add.reduce((arr - arrmean) ** 2, axis) / (n - ddof)
-    if not keepdims and isinstance(rcount, mu.ndarray):
-        rcount = rcount.squeeze(axis=axis)
-    rcount -= ddof
-    if isinstance(ret, mu.ndarray):
-        um.true_divide(ret, rcount, out=ret, casting='unsafe')
-    else:
-        ret = ret / float(rcount)
-
-    return ret
diff --git a/numpy/core/src/multiarray/methods.c b/numpy/core/src/multiarray/methods.c
index 5cff8f767..afeae9a30 100644
--- a/numpy/core/src/multiarray/methods.c
+++ b/numpy/core/src/multiarray/methods.c
@@ -49,6 +49,54 @@ NpyArg_ParseKeywords(PyObject *keys, const char *format, char **kwlist, ...)
     return ret;
 }
 
+/*
+ * Forwards an ndarray method to a the Python function
+ * numpy.core._methods.<name>(...)
+ */
+static PyObject *
+forward_ndarray_method(PyArrayObject *self, PyObject *args, PyObject *kwds,
+                            const char *name)
+{
+    PyObject *sargs, *ret;
+    PyObject *module_methods, *callable;
+    int i, n;
+
+    /* Get a reference to the function we're calling */
+    module_methods = PyImport_ImportModule("numpy.core._methods");
+    if (module_methods == NULL) {
+        return NULL;
+    }
+    callable = PyDict_GetItemString(PyModule_GetDict(module_methods), name);
+    if (callable == NULL) {
+        Py_DECREF(module_methods);
+        PyErr_Format(PyExc_RuntimeError,
+                "NumPy internal error: could not find function "
+                "numpy.core._methods.%s", name);
+        return NULL;
+    }
+
+    /* Combine 'self' and 'args' together into one tuple */
+    n = PyTuple_GET_SIZE(args);
+    sargs = PyTuple_New(n + 1);
+    if (sargs == NULL) {
+        Py_DECREF(module_methods);
+        return NULL;
+    }
+    Py_INCREF(self);
+    PyTuple_SET_ITEM(sargs, 0, (PyObject *)self);
+    for (i = 0; i < n; ++i) {
+        PyObject *item = PyTuple_GET_ITEM(args, i);
+        Py_INCREF(item);
+        PyTuple_SET_ITEM(sargs, i+1, item);
+    }
+
+    /* Call the function and return */
+    ret = PyObject_Call(callable, sargs, kwds);
+    Py_DECREF(sargs);
+    Py_DECREF(module_methods);
+    return ret;
+}
+
 static PyObject *
 array_take(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
@@ -292,36 +340,17 @@ array_argmin(PyArrayObject *self, PyObject *args, PyObject *kwds)
 static PyObject *
 array_max(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
-    int axis = MAX_DIMS;
-    PyArrayObject *out = NULL;
-    static char *kwlist[] = {"axis", "out", NULL};
-
-    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist,
-                                     PyArray_AxisConverter, &axis,
-                                     PyArray_OutputConverter, &out))
-        return NULL;
-
-    return PyArray_Max(self, axis, out);
+    return forward_ndarray_method(self, args, kwds, "_amax");
 }
 
 static PyObject *
-array_ptp(PyArrayObject *self, PyObject *args, PyObject *kwds)
+array_min(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
-    int axis = MAX_DIMS;
-    PyArrayObject *out = NULL;
-    static char *kwlist[] = {"axis", "out", NULL};
-
-    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&", kwlist,
-                                     PyArray_AxisConverter, &axis,
-                                     PyArray_OutputConverter, &out))
-        return NULL;
-
-    return PyArray_Ptp(self, axis, out);
+    return forward_ndarray_method(self, args, kwds, "_amin");
 }
 
-
 static PyObject *
-array_min(PyArrayObject *self, PyObject *args, PyObject *kwds)
+array_ptp(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
     int axis = MAX_DIMS;
     PyArrayObject *out = NULL;
@@ -332,9 +361,10 @@ array_min(PyArrayObject *self, PyObject *args, PyObject *kwds)
                                      PyArray_OutputConverter, &out))
         return NULL;
 
-    return PyArray_Min(self, axis, out);
+    return PyArray_Ptp(self, axis, out);
 }
 
+
 static PyObject *
 array_swapaxes(PyArrayObject *self, PyObject *args)
 {
@@ -1855,45 +1885,13 @@ _get_type_num_double(PyArray_Descr *dtype1, PyArray_Descr *dtype2)
 static PyObject *
 array_mean(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
-    int axis = MAX_DIMS;
-    PyArray_Descr *dtype = NULL;
-    PyArrayObject *out = NULL;
-    int num;
-    static char *kwlist[] = {"axis", "dtype", "out", NULL};
-
-    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&", kwlist,
-                                     PyArray_AxisConverter, &axis,
-                                     PyArray_DescrConverter2, &dtype,
-                                     PyArray_OutputConverter, &out)) {
-        Py_XDECREF(dtype);
-        return NULL;
-    }
-
-    num = _get_type_num_double(PyArray_DESCR(self), dtype);
-    Py_XDECREF(dtype);
-    return PyArray_Mean(self, axis, num, out);
+    return forward_ndarray_method(self, args, kwds, "_mean");
 }
 
 static PyObject *
 array_sum(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
-    int axis = MAX_DIMS;
-    PyArray_Descr *dtype = NULL;
-    PyArrayObject *out = NULL;
-    int rtype;
-    static char *kwlist[] = {"axis", "dtype", "out", NULL};
-
-    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&", kwlist,
-                                     PyArray_AxisConverter, &axis,
-                                     PyArray_DescrConverter2, &dtype,
-                                     PyArray_OutputConverter, &out)) {
-        Py_XDECREF(dtype);
-        return NULL;
-    }
-
-    rtype = _CHKTYPENUM(dtype);
-    Py_XDECREF(dtype);
-    return PyArray_Sum(self, axis, rtype, out);
+    return forward_ndarray_method(self, args, kwds, "_sum");
 }
 
 
@@ -1922,23 +1920,7 @@ array_cumsum(PyArrayObject *self, PyObject *args, PyObject *kwds)
 static PyObject *
 array_prod(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
-    int axis = MAX_DIMS;
-    PyArray_Descr *dtype = NULL;
-    PyArrayObject *out = NULL;
-    int rtype;
-    static char *kwlist[] = {"axis", "dtype", "out", NULL};
-
-    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&", kwlist,
-                                     PyArray_AxisConverter, &axis,
-                                     PyArray_DescrConverter2, &dtype,
-                                     PyArray_OutputConverter, &out)) {
-        Py_XDECREF(dtype);
-        return NULL;
-    }
-
-    rtype = _CHKTYPENUM(dtype);
-    Py_XDECREF(dtype);
-    return PyArray_Prod(self, axis, rtype, out);
+    return forward_ndarray_method(self, args, kwds, "_prod");
 }
 
 static PyObject *
@@ -2022,50 +2004,14 @@ array_all(PyArrayObject *self, PyObject *args, PyObject *kwds)
 static PyObject *
 array_stddev(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
-    int axis = MAX_DIMS;
-    PyArray_Descr *dtype = NULL;
-    PyArrayObject *out = NULL;
-    int num;
-    int ddof = 0;
-    static char *kwlist[] = {"axis", "dtype", "out", "ddof", NULL};
-
-    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&i", kwlist,
-                                     PyArray_AxisConverter, &axis,
-                                     PyArray_DescrConverter2, &dtype,
-                                     PyArray_OutputConverter, &out,
-                                     &ddof)) {
-        Py_XDECREF(dtype);
-        return NULL;
-    }
-
-    num = _get_type_num_double(PyArray_DESCR(self), dtype);
-    Py_XDECREF(dtype);
-    return __New_PyArray_Std(self, axis, num, out, 0, ddof);
+    return forward_ndarray_method(self, args, kwds, "_std");
 }
 
 
 static PyObject *
 array_variance(PyArrayObject *self, PyObject *args, PyObject *kwds)
 {
-    int axis = MAX_DIMS;
-    PyArray_Descr *dtype = NULL;
-    PyArrayObject *out = NULL;
-    int num;
-    int ddof = 0;
-    static char *kwlist[] = {"axis", "dtype", "out", "ddof", NULL};
-
-    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|O&O&O&i", kwlist,
-                                     PyArray_AxisConverter, &axis,
-                                     PyArray_DescrConverter2, &dtype,
-                                     PyArray_OutputConverter, &out,
-                                     &ddof)) {
-        Py_XDECREF(dtype);
-        return NULL;
-    }
-
-    num = _get_type_num_double(PyArray_DESCR(self), dtype);
-    Py_XDECREF(dtype);
-    return __New_PyArray_Std(self, axis, num, out, 1, ddof);
+    return forward_ndarray_method(self, args, kwds, "_var");
 }
 
 
diff --git a/numpy/core/src/multiarray/reduction.c b/numpy/core/src/multiarray/reduction.c
index d60b07fce..94ca111a9 100644
--- a/numpy/core/src/multiarray/reduction.c
+++ b/numpy/core/src/multiarray/reduction.c
@@ -366,6 +366,13 @@ initialize_reduce_result_noidentity_skipna(
         goto fail;
     }
 
+    /*
+     * Track how many initializations we've done, both to
+     * short circuit completion and to raise an error if
+     * any remained uninitialized.
+     */
+    initialized_countdown = PyArray_SIZE(result);
+
     if (NpyIter_GetIterSize(iter) != 0) {
         NpyIter_IterNextFunc *iternext;
         char **dataptr;
@@ -381,13 +388,6 @@ initialize_reduce_result_noidentity_skipna(
         strideptr = NpyIter_GetInnerStrideArray(iter);
         countptr = NpyIter_GetInnerLoopSizePtr(iter);
 
-        /*
-         * Track how many initializations we've done, both to
-         * short circuit completion and to raise an error if
-         * any remained uninitialized.
-         */
-        initialized_countdown = PyArray_SIZE(result);
-
         if (!needs_api) {
             NPY_BEGIN_THREADS;
         }
@@ -870,7 +870,7 @@ PyArray_ReduceWrapper(PyArrayObject *operand, PyArrayObject *out,
          */
         if (!reorderable && check_nonreorderable_axes(PyArray_NDIM(operand),
                                         axis_flags, funcname) < 0) {
-            return NULL;
+            goto fail;
         }
 
         if (assign_identity(result, !skipna, data) < 0) {
@@ -884,8 +884,7 @@ PyArray_ReduceWrapper(PyArrayObject *operand, PyArrayObject *out,
                             axis_flags, reorderable, skipna,
                             &skip_first_count, funcname);
         if (op_view == NULL) {
-            Py_DECREF(result);
-            return NULL;
+            goto fail;
         }
         if (PyArray_SIZE(op_view) == 0) {
             Py_DECREF(op_view);
@@ -943,9 +942,7 @@ PyArray_ReduceWrapper(PyArrayObject *operand, PyArrayObject *out,
                                op_dtypes,
                                0, NULL, NULL, buffersize);
     if (iter == NULL) {
-        Py_DECREF(result);
-        Py_DECREF(op_dtypes[0]);
-        return NULL;
+        goto fail;
     }
 
     if (NpyIter_GetIterSize(iter) != 0) {
@@ -957,10 +954,7 @@ PyArray_ReduceWrapper(PyArrayObject *operand, PyArrayObject *out,
 
         iternext = NpyIter_GetIterNext(iter, NULL);
         if (iternext == NULL) {
-            Py_DECREF(result);
-            Py_DECREF(op_dtypes[0]);
-            NpyIter_Deallocate(iter);
-            return NULL;
+            goto fail;
         }
         dataptr = NpyIter_GetDataPtrArray(iter);
         strideptr = NpyIter_GetInnerStrideArray(iter);
@@ -1000,6 +994,7 @@ PyArray_ReduceWrapper(PyArrayObject *operand, PyArrayObject *out,
     }
 
     NpyIter_Deallocate(iter);
+    Py_DECREF(op_view);
 
 finish:
     /* Strip out the extra 'one' dimensions in the result */
diff --git a/numpy/core/src/umath/ufunc_object.c b/numpy/core/src/umath/ufunc_object.c
index 9de0e7ef5..46513cda1 100644
--- a/numpy/core/src/umath/ufunc_object.c
+++ b/numpy/core/src/umath/ufunc_object.c
@@ -2536,6 +2536,7 @@ reduce_type_resolver(PyUFuncObject *ufunc, PyArrayObject *arr,
     retcode = ufunc->type_resolver(
                         ufunc, NPY_UNSAFE_CASTING,
                         op, type_tup, dtypes);
+    Py_DECREF(type_tup);
     if (retcode == -1) {
         return -1;
     }
@@ -2806,6 +2807,8 @@ PyUFunc_Reduce(PyUFuncObject *ufunc, PyArrayObject *arr, PyArrayObject *out,
     int buffersize = 0, errormask = 0;
     PyObject *errobj = NULL;
 
+    NPY_UF_DBG_PRINT1("\nEvaluating ufunc %s.reduce\n", ufunc_name);
+
     ndim = PyArray_NDIM(arr);
 
     /* Create an array of flags for reduction */
@@ -3628,7 +3631,7 @@ static PyObject *
 PyUFunc_GenericReduction(PyUFuncObject *ufunc, PyObject *args,
                          PyObject *kwds, int operation)
 {
-    int i, naxes=0;
+    int i, naxes=0, ndim;
     int axes[NPY_MAXDIMS];
     PyObject *axes_in = NULL;
     PyArrayObject *mp, *ret = NULL;
@@ -3712,6 +3715,9 @@ PyUFunc_GenericReduction(PyUFuncObject *ufunc, PyObject *args,
     if (mp == NULL) {
         return NULL;
     }
+
+    ndim = PyArray_NDIM(mp);
+
     /* Check to see that type (and otype) is not FLEXIBLE */
     if (PyArray_ISFLEXIBLE(mp) ||
         (otype && PyTypeNum_ISFLEXIBLE(otype->type_num))) {
@@ -3730,7 +3736,7 @@ PyUFunc_GenericReduction(PyUFuncObject *ufunc, PyObject *args,
     }
     /* Convert 'None' into all the axes */
     else if (axes_in == Py_None) {
-        naxes = PyArray_NDIM(mp);
+        naxes = ndim;
         for (i = 0; i < naxes; ++i) {
             axes[i] = i;
         }
@@ -3753,9 +3759,9 @@ PyUFunc_GenericReduction(PyUFuncObject *ufunc, PyObject *args,
                 return NULL;
             }
             if (axis < 0) {
-                axis += PyArray_NDIM(mp);
+                axis += ndim;
             }
-            if (axis < 0 || axis >= PyArray_NDIM(mp)) {
+            if (axis < 0 || axis >= ndim) {
                 PyErr_SetString(PyExc_ValueError,
                         "'axis' entry is out of bounds");
                 Py_XDECREF(otype);
@@ -3775,11 +3781,13 @@ PyUFunc_GenericReduction(PyUFuncObject *ufunc, PyObject *args,
             return NULL;
         }
         if (axis < 0) {
-            axis += PyArray_NDIM(mp);
+            axis += ndim;
+        }
+        /* Special case letting axis={0 or -1} slip through for scalars */
+        if (ndim == 0 && (axis == 0 || axis == -1)) {
+            axis = 0;
         }
-        /* Special case letting axis=0 slip through for scalars */
-        if ((axis < 0 || axis >= PyArray_NDIM(mp)) &&
-                !(axis == 0 && PyArray_NDIM(mp) == 0)) {
+        else if (axis < 0 || axis >= ndim) {
             PyErr_SetString(PyExc_ValueError,
                     "'axis' entry is out of bounds");
             Py_XDECREF(otype);
@@ -3791,7 +3799,7 @@ PyUFunc_GenericReduction(PyUFuncObject *ufunc, PyObject *args,
     }
 
     /* Check to see if input is zero-dimensional. */
-    if (PyArray_NDIM(mp) == 0) {
+    if (ndim == 0) {
         /*
          * A reduction with no axes is still valid but trivial.
          * As a special case for backwards compatibility in 'sum',
diff --git a/numpy/core/src/umath/ufunc_type_resolution.c b/numpy/core/src/umath/ufunc_type_resolution.c
index de1da5c69..0d6cf19f1 100644
--- a/numpy/core/src/umath/ufunc_type_resolution.c
+++ b/numpy/core/src/umath/ufunc_type_resolution.c
@@ -217,6 +217,9 @@ PyUFunc_SimpleBinaryComparisonTypeResolver(PyUFuncObject *ufunc,
         Py_INCREF(out_dtypes[1]);
     }
     else {
+        PyObject *item;
+        PyArray_Descr *dtype = NULL;
+
         /*
          * If the type tuple isn't a single-element tuple, let the
          * default type resolution handle this one.
@@ -226,14 +229,18 @@ PyUFunc_SimpleBinaryComparisonTypeResolver(PyUFuncObject *ufunc,
                     operands, type_tup, out_dtypes);
         }
 
-        if (!PyArray_DescrCheck(PyTuple_GET_ITEM(type_tup, 0))) {
+        item = PyTuple_GET_ITEM(type_tup, 0);
+
+        if (item == Py_None) {
             PyErr_SetString(PyExc_ValueError,
                     "require data type in the type tuple");
             return -1;
         }
+        else if (!PyArray_DescrConverter(item, &dtype)) {
+            return -1;
+        }
 
-        out_dtypes[0] = ensure_dtype_nbo(
-                            (PyArray_Descr *)PyTuple_GET_ITEM(type_tup, 0));
+        out_dtypes[0] = ensure_dtype_nbo(dtype);
         if (out_dtypes[0] == NULL) {
             return -1;
         }
@@ -314,6 +321,9 @@ PyUFunc_SimpleUnaryOperationTypeResolver(PyUFuncObject *ufunc,
         Py_INCREF(out_dtypes[1]);
     }
     else {
+        PyObject *item;
+        PyArray_Descr *dtype = NULL;
+
         /*
          * If the type tuple isn't a single-element tuple, let the
          * default type resolution handle this one.
@@ -323,14 +333,18 @@ PyUFunc_SimpleUnaryOperationTypeResolver(PyUFuncObject *ufunc,
                     operands, type_tup, out_dtypes);
         }
 
-        if (!PyArray_DescrCheck(PyTuple_GET_ITEM(type_tup, 0))) {
+        item = PyTuple_GET_ITEM(type_tup, 0);
+
+        if (item == Py_None) {
             PyErr_SetString(PyExc_ValueError,
                     "require data type in the type tuple");
             return -1;
         }
+        else if (!PyArray_DescrConverter(item, &dtype)) {
+            return -1;
+        }
 
-        out_dtypes[0] = ensure_dtype_nbo(
-                            (PyArray_Descr *)PyTuple_GET_ITEM(type_tup, 0));
+        out_dtypes[0] = ensure_dtype_nbo(dtype);
         if (out_dtypes[0] == NULL) {
             return -1;
         }
@@ -424,6 +438,9 @@ PyUFunc_SimpleBinaryOperationTypeResolver(PyUFuncObject *ufunc,
         Py_INCREF(out_dtypes[2]);
     }
     else {
+        PyObject *item;
+        PyArray_Descr *dtype = NULL;
+
         /*
          * If the type tuple isn't a single-element tuple, let the
          * default type resolution handle this one.
@@ -433,14 +450,18 @@ PyUFunc_SimpleBinaryOperationTypeResolver(PyUFuncObject *ufunc,
                     operands, type_tup, out_dtypes);
         }
 
-        if (!PyArray_DescrCheck(PyTuple_GET_ITEM(type_tup, 0))) {
+        item = PyTuple_GET_ITEM(type_tup, 0);
+
+        if (item == Py_None) {
             PyErr_SetString(PyExc_ValueError,
                     "require data type in the type tuple");
             return -1;
         }
+        else if (!PyArray_DescrConverter(item, &dtype)) {
+            return -1;
+        }
 
-        out_dtypes[0] = ensure_dtype_nbo(
-                            (PyArray_Descr *)PyTuple_GET_ITEM(type_tup, 0));
+        out_dtypes[0] = ensure_dtype_nbo(dtype);
         if (out_dtypes[0] == NULL) {
             return -1;
         }
@@ -591,8 +612,8 @@ PyUFunc_AdditionTypeResolver(PyUFuncObject *ufunc,
 
     /* Use the default when datetime and timedelta are not involved */
     if (!PyTypeNum_ISDATETIME(type_num1) && !PyTypeNum_ISDATETIME(type_num2)) {
-        return PyUFunc_DefaultTypeResolver(ufunc, casting, operands,
-                    type_tup, out_dtypes);
+        return PyUFunc_SimpleBinaryOperationTypeResolver(ufunc, casting,
+                    operands, type_tup, out_dtypes);
     }
 
     if (type_num1 == NPY_TIMEDELTA) {
@@ -780,8 +801,8 @@ PyUFunc_SubtractionTypeResolver(PyUFuncObject *ufunc,
 
     /* Use the default when datetime and timedelta are not involved */
     if (!PyTypeNum_ISDATETIME(type_num1) && !PyTypeNum_ISDATETIME(type_num2)) {
-        return PyUFunc_DefaultTypeResolver(ufunc, casting, operands,
-                    type_tup, out_dtypes);
+        return PyUFunc_SimpleBinaryOperationTypeResolver(ufunc, casting,
+                    operands, type_tup, out_dtypes);
     }
 
     if (type_num1 == NPY_TIMEDELTA) {
@@ -947,8 +968,8 @@ PyUFunc_MultiplicationTypeResolver(PyUFuncObject *ufunc,
 
     /* Use the default when datetime and timedelta are not involved */
     if (!PyTypeNum_ISDATETIME(type_num1) && !PyTypeNum_ISDATETIME(type_num2)) {
-        return PyUFunc_DefaultTypeResolver(ufunc, casting, operands,
-                    type_tup, out_dtypes);
+        return PyUFunc_SimpleBinaryOperationTypeResolver(ufunc, casting,
+                    operands, type_tup, out_dtypes);
     }
 
     if (type_num1 == NPY_TIMEDELTA) {
diff --git a/numpy/core/tests/test_regression.py b/numpy/core/tests/test_regression.py
index 2bad0e90c..03c5f026b 100644
--- a/numpy/core/tests/test_regression.py
+++ b/numpy/core/tests/test_regression.py
@@ -1626,5 +1626,11 @@ class TestRegression(TestCase):
             a = np.empty((100000000,), dtype='i1')
             del a
 
+    def test_ufunc_reduce_memoryleak(self):
+        a = np.arange(6)
+        acnt = sys.getrefcount(a)
+        res = np.add.reduce(a)
+        assert_equal(sys.getrefcount(a), acnt)
+
 if __name__ == "__main__":
     run_module_suite()