BUG: core: Fix things so scipy trunk passes all tests (but one)

With this patch, the latest scipy trunk (7087), built against NumPy
1.5.1, passes all tests when run against the numpy trunk.  The single
failing test, test_imresize, fails because it tests all float types,
and the new 'half' type lacks the precision to pass that test.

8 files changed, 112 insertions, 64 deletions

diff --git a/numpy/core/src/multiarray/convert_datatype.c b/numpy/core/src/multiarray/convert_datatype.c
index 45a20b5de..a2ce0ccf3 100644
--- a/numpy/core/src/multiarray/convert_datatype.c
+++ b/numpy/core/src/multiarray/convert_datatype.c
@@ -1007,31 +1007,50 @@ static int min_scalar_type_num(char *valueptr, int type_num,
             break;
         }
         /*
-         * Complex types may be demoted to float types if the
-         * imaginary part is zero.
+         * The code to demote complex to float is disabled for now,
+         * as forcing complex by adding 0j is probably desireable.
          */
         case NPY_CFLOAT: {
             npy_cfloat value = *(npy_cfloat *)valueptr;
+            /*
             if (value.imag == 0) {
-                return min_scalar_type_num((char *)&value.real, NPY_FLOAT, is_small_unsigned);
+                return min_scalar_type_num((char *)&value.real,
+                                            NPY_FLOAT, is_small_unsigned);
             }
+            */
             break;
         }
         case NPY_CDOUBLE: {
             npy_cdouble value = *(npy_cdouble *)valueptr;
+            /*
             if (value.imag == 0) {
-                return min_scalar_type_num((char *)&value.real, NPY_DOUBLE, is_small_unsigned);
+                return min_scalar_type_num((char *)&value.real,
+                                            NPY_DOUBLE, is_small_unsigned);
             }
-            /* TODO: Check overflow values as for float case */
-            return NPY_CFLOAT;
+            */
+            if (value.real > -3.4e38 && value.real < 3.4e38 &&
+                     value.imag > -3.4e38 && value.imag < 3.4e38) {
+                return NPY_CFLOAT;
+            }
+            break;
         }
         case NPY_CLONGDOUBLE: {
             npy_cdouble value = *(npy_cdouble *)valueptr;
+            /*
             if (value.imag == 0) {
-                return min_scalar_type_num((char *)&value.real, NPY_LONGDOUBLE, is_small_unsigned);
+                return min_scalar_type_num((char *)&value.real,
+                                            NPY_LONGDOUBLE, is_small_unsigned);
+            }
+            */
+            if (value.real > -3.4e38 && value.real < 3.4e38 &&
+                     value.imag > -3.4e38 && value.imag < 3.4e38) {
+                return NPY_CFLOAT;
             }
-            /* TODO: Check overflow values as for float case */
-            return NPY_CFLOAT;
+            else if (value.real > -1.7e308 && value.real < 1.7e308 &&
+                     value.imag > -1.7e308 && value.imag < 1.7e308) {
+                return NPY_CDOUBLE;
+            }
+            break;
         }
     }
 
diff --git a/numpy/core/src/multiarray/item_selection.c b/numpy/core/src/multiarray/item_selection.c
index cc8e6aa8a..293b26326 100644
--- a/numpy/core/src/multiarray/item_selection.c
+++ b/numpy/core/src/multiarray/item_selection.c
@@ -1819,6 +1819,7 @@ PyArray_Nonzero(PyArrayObject *self)
     /* Build an iterator with coordinates, in C order */
     iter = NpyIter_New(self, NPY_ITER_READONLY|
                              NPY_ITER_COORDS|
+                             NPY_ITER_ZEROSIZE_OK|
                              NPY_ITER_REFS_OK,
                         NPY_CORDER, NPY_NO_CASTING,
                         NULL, 0, NULL, 0);
@@ -1828,31 +1829,33 @@ PyArray_Nonzero(PyArrayObject *self)
         return NULL;
     }
 
-    /* Get the pointers for inner loop iteration */
-    iternext = NpyIter_GetIterNext(iter, NULL);
-    if (iternext == NULL) {
-        NpyIter_Deallocate(iter);
-        Py_DECREF(ret);
-        return NULL;
-    }
-    getcoords = NpyIter_GetGetCoords(iter, NULL);
-    if (getcoords == NULL) {
-        NpyIter_Deallocate(iter);
-        Py_DECREF(ret);
-        return NULL;
-    }
-    dataptr = NpyIter_GetDataPtrArray(iter);
-    innersizeptr = NpyIter_GetInnerLoopSizePtr(iter);
+    if (NpyIter_GetIterSize(iter) != 0) {
+        /* Get the pointers for inner loop iteration */
+        iternext = NpyIter_GetIterNext(iter, NULL);
+        if (iternext == NULL) {
+            NpyIter_Deallocate(iter);
+            Py_DECREF(ret);
+            return NULL;
+        }
+        getcoords = NpyIter_GetGetCoords(iter, NULL);
+        if (getcoords == NULL) {
+            NpyIter_Deallocate(iter);
+            Py_DECREF(ret);
+            return NULL;
+        }
+        dataptr = NpyIter_GetDataPtrArray(iter);
+        innersizeptr = NpyIter_GetInnerLoopSizePtr(iter);
 
-    coords = (npy_intp *)PyArray_DATA(ret);
+        coords = (npy_intp *)PyArray_DATA(ret);
 
-    /* Get the coordinates for each non-zero element */
-    do {
-        if (nonzero(*dataptr, self)) {
-            getcoords(iter, coords);
-            coords += ndim;
-        }
-    } while(iternext(iter));
+        /* Get the coordinates for each non-zero element */
+        do {
+            if (nonzero(*dataptr, self)) {
+                getcoords(iter, coords);
+                coords += ndim;
+            }
+        } while(iternext(iter));
+    }
 
     NpyIter_Deallocate(iter);
 
diff --git a/numpy/core/src/multiarray/lowlevel_strided_loops.c.src b/numpy/core/src/multiarray/lowlevel_strided_loops.c.src
index 24265c3cf..70dd3437b 100644
--- a/numpy/core/src/multiarray/lowlevel_strided_loops.c.src
+++ b/numpy/core/src/multiarray/lowlevel_strided_loops.c.src
@@ -753,6 +753,8 @@ NPY_NO_EXPORT PyArray_StridedTransferFn *
 #    define _CONVERT_FN(x) npy_halfbits_to_doublebits(x)
 #  elif @is_half2@
 #    define _CONVERT_FN(x) (x)
+#  elif @is_bool2@
+#    define _CONVERT_FN(x) ((npy_bool)!npy_half_iszero(x))
 #  else
 #    define _CONVERT_FN(x) ((_TYPE2)npy_half_to_float(x))
 #  endif
diff --git a/numpy/core/src/multiarray/new_iterator.c.src b/numpy/core/src/multiarray/new_iterator.c.src
index 8c302e161..bc819b150 100644
--- a/numpy/core/src/multiarray/new_iterator.c.src
+++ b/numpy/core/src/multiarray/new_iterator.c.src
@@ -274,7 +274,8 @@ npyiter_check_casting(npy_intp niter, PyArrayObject **op,
 static int
 npyiter_fill_axisdata(NpyIter *iter, npy_uint32 flags, char *op_itflags,
                     char **op_dataptr,
-                    npy_uint32 *op_flags, npy_intp **op_axes);
+                    npy_uint32 *op_flags, npy_intp **op_axes,
+                    int output_scalars);
 static void
 npyiter_replace_axisdata(NpyIter *iter, npy_intp iiter,
                       PyArrayObject *op,
@@ -429,7 +430,7 @@ NpyIter_MultiNew(npy_intp niter, PyArrayObject **op_in, npy_uint32 flags,
 
     /* Fill in the AXISDATA arrays and set the ITERSIZE field */
     if (!npyiter_fill_axisdata(iter, flags, op_itflags, op_dataptr,
-                                                op_flags, op_axes)) {
+                                        op_flags, op_axes, output_scalars)) {
         NpyIter_Deallocate(iter);
         return NULL;
     }
@@ -2985,7 +2986,8 @@ npyiter_shape_string(npy_intp n, npy_intp *vals, char *ending)
 static int
 npyiter_fill_axisdata(NpyIter *iter, npy_uint32 flags, char *op_itflags,
                     char **op_dataptr,
-                    npy_uint32 *op_flags, npy_intp **op_axes)
+                    npy_uint32 *op_flags, npy_intp **op_axes,
+                    int output_scalars)
 {
     npy_uint32 itflags = NIT_ITFLAGS(iter);
     npy_intp idim, ndim = NIT_NDIM(iter);
@@ -3188,8 +3190,12 @@ npyiter_fill_axisdata(NpyIter *iter, npy_uint32 flags, char *op_itflags,
     /* Go through and check for operands for broadcasting and reduction */
     for (iiter = 0; iiter < niter; ++iiter) {
         if (op[iiter] != NULL) {
-            /* If broadcasting is disallowed for this operand */
-            if (op_flags[iiter]&NPY_ITER_NO_BROADCAST) {
+            /*
+             * If broadcasting is disallowed for this operand,
+             * unless scalars are output, which means all operands are scalar
+             * and no broadcasting errors could occur
+             */
+            if ((op_flags[iiter]&NPY_ITER_NO_BROADCAST) && !output_scalars) {
                 npy_intp *axes;
 
                 axes = op_axes ? op_axes[iiter] : NULL;
diff --git a/numpy/core/src/umath/ufunc_object.c b/numpy/core/src/umath/ufunc_object.c
index f4dc1185c..1dc8354ee 100644
--- a/numpy/core/src/umath/ufunc_object.c
+++ b/numpy/core/src/umath/ufunc_object.c
@@ -2717,33 +2717,34 @@ PyUFunc_ReductionOp(PyUFuncObject *self, PyArrayObject *arr,
      * unit for UFUNC_REDUCE, or return the zero-sized output array
      * for UFUNC_ACCUMULATE.
      */
-    if (PyArray_DIM(op[1], axis) == 0) {
-        if (operation == UFUNC_REDUCE) {
-            if (self->identity == PyUFunc_None) {
-                PyErr_Format(PyExc_ValueError,
-                             "zero-size array to %s.%s "
-                             "without identity", ufunc_name, opname);
+    if (operation == UFUNC_REDUCE && PyArray_DIM(op[1], axis) == 0) {
+        if (self->identity == PyUFunc_None) {
+            PyErr_Format(PyExc_ValueError,
+                         "zero-size array to %s.%s "
+                         "without identity", ufunc_name, opname);
+            goto fail;
+        }
+        if (self->identity == PyUFunc_One) {
+            PyObject *obj = PyInt_FromLong((long) 1);
+            if (obj == NULL) {
                 goto fail;
             }
-            if (self->identity == PyUFunc_One) {
-                PyObject *obj = PyInt_FromLong((long) 1);
-                if (obj == NULL) {
-                    goto fail;
-                }
-                PyArray_FillWithScalar(op[0], obj);
-                Py_DECREF(obj);
-            } else {
-                PyObject *obj = PyInt_FromLong((long) 0);
-                if (obj == NULL) {
-                    goto fail;
-                }
-                PyArray_FillWithScalar(op[0], obj);
-                Py_DECREF(obj);
+            PyArray_FillWithScalar(op[0], obj);
+            Py_DECREF(obj);
+        } else {
+            PyObject *obj = PyInt_FromLong((long) 0);
+            if (obj == NULL) {
+                goto fail;
             }
+            PyArray_FillWithScalar(op[0], obj);
+            Py_DECREF(obj);
         }
 
         goto finish;
     }
+    else if (PyArray_SIZE(op[0]) == 0) {
+        goto finish;
+    }
 
     /* Only allocate an inner iterator if it's necessary */
     if (!PyArray_ISALIGNED(op[1]) || !PyArray_ISALIGNED(op[0]) ||
diff --git a/numpy/core/tests/test_numeric.py b/numpy/core/tests/test_numeric.py
index 2cb2ff532..743528487 100644
--- a/numpy/core/tests/test_numeric.py
+++ b/numpy/core/tests/test_numeric.py
@@ -657,12 +657,16 @@ class TestTypes(TestCase):
 
         f64 = float64(0)
         c64 = complex64(0)
-        # Scalars do not coerce to complex if the value is real
-        assert_equal(res_type(c64,array([f64])), np.dtype(float64))
+        ## Scalars do not coerce to complex if the value is real
+        #assert_equal(res_type(c64,array([f64])), np.dtype(float64))
         # But they do if the value is complex
         assert_equal(res_type(complex64(3j),array([f64])),
                                                     np.dtype(complex128))
 
+        # Scalars do coerce to complex even if the value is real
+        # This is so "a+0j" can be reliably used to make something complex.
+        assert_equal(res_type(c64,array([f64])), np.dtype(complex128))
+
         ctx.__exit__()
 
 
@@ -671,12 +675,17 @@ class TestTypes(TestCase):
 
         f64 = float64(0)
         c64 = complex64(0)
-        # Scalars do not coerce to complex if the value is real
-        assert_equal(np.result_type(c64,array([f64])), np.dtype(float64))
+        ## Scalars do not coerce to complex if the value is real
+        #assert_equal(np.result_type(c64,array([f64])), np.dtype(float64))
         # But they do if the value is complex
         assert_equal(np.result_type(complex64(3j),array([f64])),
                                                     np.dtype(complex128))
 
+        # Scalars do coerce to complex even if the value is real
+        # This is so "a+0j" can be reliably used to make something complex.
+        assert_equal(np.result_type(c64,array([f64])), np.dtype(complex128))
+
+
     def can_cast(self):
         assert_(np.can_cast(np.int32, np.int64))
         assert_(np.can_cast(np.float64, np.complex))
diff --git a/numpy/core/tests/test_regression.py b/numpy/core/tests/test_regression.py
index 8439c9b80..7c0b24953 100644
--- a/numpy/core/tests/test_regression.py
+++ b/numpy/core/tests/test_regression.py
@@ -1471,5 +1471,12 @@ class TestRegression(TestCase):
         a = np.arange(10000, dtype='f')
         assert_equal(a.sum(dtype='d'), a.astype('d').sum())
 
+    def test_ufunc_casting_out(self):
+        a = np.array(1.0, dtype=np.float32)
+        b = np.array(1.0, dtype=np.float64)
+        c = np.array(1.0, dtype=np.float32)
+        np.add(a, b, out=c)
+        assert_equal(c, 2.0)
+
 if __name__ == "__main__":
     run_module_suite()
diff --git a/numpy/testing/utils.py b/numpy/testing/utils.py
index 6f18b5468..1a03f98ea 100644
--- a/numpy/testing/utils.py
+++ b/numpy/testing/utils.py
@@ -421,7 +421,7 @@ def assert_almost_equal(actual,desired,decimal=7,err_msg='',verbose=True):
         usecomplex = False
 
     msg = build_err_msg([actual, desired], err_msg, verbose=verbose,
-                         header='Arrays are not almost equal')
+             header=('Arrays are not almost equal to %d decimals' % decimal))
 
     if usecomplex:
         if iscomplexobj(actual):
@@ -616,7 +616,8 @@ def assert_array_compare(comparison, x, y, err_msg='', verbose=True,
                                 names=('x', 'y'))
             if not cond :
                 raise AssertionError(msg)
-    except ValueError:
+    except ValueError as e:
+        header = 'error during assertion:\n%s\n\n%s' % (e, header)
         msg = build_err_msg([x, y], err_msg, verbose=verbose, header=header,
                             names=('x', 'y'))
         raise ValueError(msg)
@@ -771,7 +772,7 @@ def assert_array_almost_equal(x, y, decimal=6, err_msg='', verbose=True):
             z = z.astype(float_) # handle object arrays
         return around(z, decimal) <= 10.0**(-decimal)
     assert_array_compare(compare, x, y, err_msg=err_msg, verbose=verbose,
-                         header='Arrays are not almost equal')
+             header=('Arrays are not almost equal to %d decimals' % decimal))
 
 def assert_array_less(x, y, err_msg='', verbose=True):
     """