Change way maximum and minimum deal with nans.

Add ufuncs fmax and fmin.

In the following, a complex number is considered a nan if the real or imaginary
part is nan. This means that there are many different complex numbers that are
nans and this effects the nan values returned by the maximum, minimum, fmax, and
fmin ufuncs.

The maximum and minimum ufuncs are the same as before unless nans are involved.
In the case of nans, if both values being compared are nans, then the first is
returned, otherwise the nan value is returned. This has the effect of
propagating nans.

The fmax and fmin ufuncs return the same values as maximum and minimum if
neither value being compared is nan.  In the case of nans, if both values being
compared are nans, then the first is returned, otherwise the non-nan value is
returned. This has the effect that nans are ignored unless both values are nan.

3 files changed, 169 insertions, 24 deletions

diff --git a/numpy/core/code_generators/generate_umath.py b/numpy/core/code_generators/generate_umath.py
index d831ceda9..a7ebd93fb 100644
--- a/numpy/core/code_generators/generate_umath.py
+++ b/numpy/core/code_generators/generate_umath.py
@@ -316,6 +316,16 @@ defdict = {
           TD(noobj),
           TD(O, f='_npy_ObjectMin')
           ),
+'fmax' :
+    Ufunc(2, 1, None,
+          "",
+          TD(inexact)
+          ),
+'fmin' :
+    Ufunc(2, 1, None,
+          "",
+          TD(inexact)
+          ),
 'bitwise_and' :
     Ufunc(2, 1, One,
           docstrings.get('numpy.core.umath.bitwise_and'),
diff --git a/numpy/core/src/umathmodule.c.src b/numpy/core/src/umathmodule.c.src
index af66c256b..cc413085d 100644
--- a/numpy/core/src/umathmodule.c.src
+++ b/numpy/core/src/umathmodule.c.src
@@ -654,7 +654,7 @@ static void
 {
     UNARY_LOOP {
         const @s@@type@ in1 = *(@s@@type@ *)ip1;
-        *((@s@@type@ *)op) = 1.0/in1;
+        *((@s@@type@ *)op) = (@s@@type@)(1.0/in1);
     }
 }
 
@@ -924,6 +924,9 @@ U@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(f
  *  #C = F, , L#
  */
 
+#define ONE 1.0@c@
+#define ZERO 0.0@c@
+
 /**begin repeat1
  * Arithmetic
  * # kind = add, subtract, multiply, divide#
@@ -991,7 +994,7 @@ static void
 
 /**begin repeat1
  * #kind = maximum, minimum#
- * #OP =  >, <#
+ * #OP =  >=, <=#
  **/
 static void
 @TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
@@ -1000,7 +1003,23 @@ static void
     BINARY_LOOP {
         const @type@ in1 = *(@type@ *)ip1;
         const @type@ in2 = *(@type@ *)ip2;
-        *((@type@ *)op) = in1 @OP@ in2 ? in1 : in2;
+        *((@type@ *)op) = (in1 @OP@ in2 || isnan(in1)) ? in1 : in2;
+    }
+}
+/**end repeat1**/
+
+/**begin repeat1
+ * #kind = fmax, fmin#
+ * #OP =  >=, <=#
+ **/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+    /*  */
+    BINARY_LOOP {
+        const @type@ in1 = *(@type@ *)ip1;
+        const @type@ in2 = *(@type@ *)ip2;
+        *((@type@ *)op) = (in1 @OP@ in2 || isnan(in2)) ? in1 : in2;
     }
 }
 /**end repeat1**/
@@ -1045,7 +1064,7 @@ static void
 {
     UNARY_LOOP {
         const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op) = 1.0/in1;
+        *((@type@ *)op) = ONE/in1;
     }
 }
 
@@ -1053,7 +1072,7 @@ static void
 @TYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
 {
     OUTPUT_LOOP {
-        *((@type@ *)op) = 1;
+        *((@type@ *)op) = ONE;
     }
 }
 
@@ -1071,7 +1090,7 @@ static void
 {
     UNARY_LOOP {
         const @type@ in1 = *(@type@ *)ip1;
-        const @type@ tmp = (in1 > 0) ? in1 : -in1;
+        const @type@ tmp = in1 > 0 ? in1 : -in1;
         /* add 0 to clear -0.0 */
         *((@type@ *)op) = tmp + 0;
     }
@@ -1138,6 +1157,9 @@ static void
 
 #define @TYPE@_true_divide @TYPE@_divide
 
+#undef ONE
+#undef ZERO
+
 /**end repeat**/
 
 
@@ -1155,6 +1177,11 @@ static void
  * #c = f, , l#
  */
 
+#define CGE(xr,xi,yr,yi) (xr > yr || (xr == yr && xi >= yi))
+#define CLE(xr,xi,yr,yi) (xr < yr || (xr == yr && xi <= yi))
+#define ONE 1.0@c@
+#define ZERO 0.0@c@
+
 /**begin repeat1
  * arithmetic
  * #kind = add, subtract#
@@ -1348,8 +1375,8 @@ static void
 @CTYPE@_ones_like(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(data))
 {
     OUTPUT_LOOP {
-        ((@type@ *)op)[0] = 1;
-        ((@type@ *)op)[1] = 0;
+        ((@type@ *)op)[0] = ONE;
+        ((@type@ *)op)[1] = ZERO;
     }
 }
 
@@ -1380,29 +1407,30 @@ static void
         const @type@ in1r = ((@type@ *)ip1)[0];
         const @type@ in1i = ((@type@ *)ip1)[1];
         if (in1r > 0) {
-            ((@type@ *)op)[0] = 1;
+            ((@type@ *)op)[0] = ONE;
         }
         else if (in1r < 0) {
-            ((@type@ *)op)[0] = -1;
+            ((@type@ *)op)[0] = -ONE;
         }
         else {
             if (in1i > 0) {
-                ((@type@ *)op)[0] = 1;
+                ((@type@ *)op)[0] = ONE;
             }
             else if (in1i < 0) {
-                ((@type@ *)op)[0] = -1;
+                ((@type@ *)op)[0] = -ONE;
             }
             else {
-                ((@type@ *)op)[0] = 0;
+                ((@type@ *)op)[0] = ZERO;
             }
         }
-        ((@type@ *)op)[1] = 0;
+        ((@type@ *)op)[1] = ZERO;
     }
 }
 
 /**begin repeat1
  * #kind = maximum, minimum#
- * #OP = >, <#
+ * #OP1 = CGE, CLE#
+ * #OP2 = CLE, CGE#
  */
 static void
 @CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
@@ -1412,7 +1440,31 @@ static void
         const @type@ in1i = ((@type@ *)ip1)[1];
         const @type@ in2r = ((@type@ *)ip2)[0];
         const @type@ in2i = ((@type@ *)ip2)[1];
-        if (in1r @OP@ in2r || ((in1r == in2r) && (in1i @OP@ in2i))) {
+        if (@OP1@(in1r, in1i, in2r, in2i) || isnan(in1r) || isnan(in1i)) {
+            ((@type@ *)op)[0] = in1r;
+            ((@type@ *)op)[1] = in1i;
+        }
+        else {
+            ((@type@ *)op)[0] = in2r;
+            ((@type@ *)op)[1] = in2i;
+        }
+    }
+}
+/**end repeat1**/
+
+/**begin repeat1
+ * #kind = fmax, fmin#
+ * #OP1 = CGE, CLE#
+ */
+static void
+@CTYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+    BINARY_LOOP {
+        const @type@ in1r = ((@type@ *)ip1)[0];
+        const @type@ in1i = ((@type@ *)ip1)[1];
+        const @type@ in2r = ((@type@ *)ip2)[0];
+        const @type@ in2i = ((@type@ *)ip2)[1];
+        if (@OP1@(in1r, in1i, in2r, in2i) || isnan(in2r) || isnan(in2i)) {
             ((@type@ *)op)[0] = in1r;
             ((@type@ *)op)[1] = in1i;
         }
@@ -1425,6 +1477,12 @@ static void
 /**end repeat1**/
 
 #define @CTYPE@_true_divide @CTYPE@_divide
+
+#undef CGE
+#undef CLE
+#undef ONE
+#undef ZERO
+
 /**end repeat**/
 
 /*
@@ -1447,6 +1505,7 @@ OBJECT_@kind@(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)
 }
 /**end repeat**/
 
+static void
 OBJECT_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
 {
     PyObject *zero = PyInt_FromLong(0);
@@ -1464,6 +1523,16 @@ OBJECT_sign(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
  */
 
 
+/*
+ *****************************************************************************
+ **                            SETUP UFUNCS                                 **
+ *****************************************************************************
+ */
+
+#include "__umath_generated.c"
+#include "ufuncobject.c"
+#include "__ufunc_api.c"
+
 static PyUFuncGenericFunction frexp_functions[] = {
 #ifdef HAVE_FREXPF
     FLOAT_frexp,
@@ -1485,7 +1554,6 @@ static char frexp_signatures[] = {
 #endif
 };
 
-
 static PyUFuncGenericFunction ldexp_functions[] = {
 #ifdef HAVE_LDEXPF
     FLOAT_ldexp,
@@ -1507,10 +1575,6 @@ static char ldexp_signatures[] = {
 };
 
 
-#include "__umath_generated.c"
-#include "ufuncobject.c"
-#include "__ufunc_api.c"
-
 static double
 pinf_init(void)
 {
diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
index ee2893f18..34df54563 100644
--- a/numpy/core/tests/test_umath.py
+++ b/numpy/core/tests/test_umath.py
@@ -47,12 +47,83 @@ class TestExpm1(TestCase):
 
 class TestMaximum(TestCase):
     def test_reduce_complex(self):
-        assert_equal(ncu.maximum.reduce([1,2j]),1)
-        assert_equal(ncu.maximum.reduce([1+3j,2j]),1+3j)
+        assert_equal(np.maximum.reduce([1,2j]),1)
+        assert_equal(np.maximum.reduce([1+3j,2j]),1+3j)
+
+    def test_float_nans(self):
+        nan = np.nan
+        arg1 = np.array([0,   nan, nan])
+        arg2 = np.array([nan, 0,   nan])
+        out  = np.array([nan, nan, nan])
+        assert_equal(np.maximum(arg1, arg2), out)
+
+    def test_complex_nans(self):
+        nan = np.nan
+        for cnan in [nan, nan*1j, nan + nan*1j] :
+            arg1 = np.array([0, cnan, cnan], dtype=np.complex)
+            arg2 = np.array([cnan, 0, cnan], dtype=np.complex)
+            out  = np.array([nan, nan, nan], dtype=np.complex)
+            assert_equal(np.maximum(arg1, arg2), out)
 
 class TestMinimum(TestCase):
     def test_reduce_complex(self):
-        assert_equal(ncu.minimum.reduce([1,2j]),2j)
+        assert_equal(np.minimum.reduce([1,2j]),2j)
+        assert_equal(np.minimum.reduce([1+3j,2j]),2j)
+
+    def test_float_nans(self):
+        nan = np.nan
+        arg1 = np.array([0,   nan, nan])
+        arg2 = np.array([nan, 0,   nan])
+        out  = np.array([nan, nan, nan])
+        assert_equal(np.minimum(arg1, arg2), out)
+
+    def test_complex_nans(self):
+        nan = np.nan
+        for cnan in [nan, nan*1j, nan + nan*1j] :
+            arg1 = np.array([0, cnan, cnan], dtype=np.complex)
+            arg2 = np.array([cnan, 0, cnan], dtype=np.complex)
+            out  = np.array([nan, nan, nan], dtype=np.complex)
+            assert_equal(np.minimum(arg1, arg2), out)
+
+class TestFmax(TestCase):
+    def test_reduce_complex(self):
+        assert_equal(np.fmax.reduce([1,2j]),1)
+        assert_equal(np.fmax.reduce([1+3j,2j]),1+3j)
+
+    def test_float_nans(self):
+        nan = np.nan
+        arg1 = np.array([0,   nan, nan])
+        arg2 = np.array([nan, 0,   nan])
+        out  = np.array([0,   0,   nan])
+        assert_equal(np.fmax(arg1, arg2), out)
+
+    def test_complex_nans(self):
+        nan = np.nan
+        for cnan in [nan, nan*1j, nan + nan*1j] :
+            arg1 = np.array([0, cnan, cnan], dtype=np.complex)
+            arg2 = np.array([cnan, 0, cnan], dtype=np.complex)
+            out  = np.array([0,    0, nan], dtype=np.complex)
+            assert_equal(np.fmax(arg1, arg2), out)
+
+class TestFmin(TestCase):
+    def test_reduce_complex(self):
+        assert_equal(np.fmin.reduce([1,2j]),2j)
+        assert_equal(np.fmin.reduce([1+3j,2j]),2j)
+
+    def test_float_nans(self):
+        nan = np.nan
+        arg1 = np.array([0,   nan, nan])
+        arg2 = np.array([nan, 0,   nan])
+        out  = np.array([0,   0,   nan])
+        assert_equal(np.fmin(arg1, arg2), out)
+
+    def test_complex_nans(self):
+        nan = np.nan
+        for cnan in [nan, nan*1j, nan + nan*1j] :
+            arg1 = np.array([0, cnan, cnan], dtype=np.complex)
+            arg2 = np.array([cnan, 0, cnan], dtype=np.complex)
+            out  = np.array([0,    0, nan], dtype=np.complex)
+            assert_equal(np.fmin(arg1, arg2), out)
 
 class TestFloatingPoint(TestCase):
     def test_floating_point(self):