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.

6 files changed, 146 insertions, 225 deletions

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()