Make min/max on Dual consistent with isless

Project.toml

@@ -1,6 +1,6 @@
 name = "ForwardDiff"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "1.4.7"
+version = "1.4.8"
 
 [deps]
 CommonSubexpressions = "bbf7d656-a473-5ed7-a52c-81e309532950"

src/dual.jl

@@ -474,7 +474,7 @@ Base.float(d::Dual) = convert(float(typeof(d)), d)
 ###################################
 
 for (M, f, arity) in DiffRules.diffrules(filter_modules = nothing)
-    if (M, f) in ((:Base, :^), (:NaNMath, :pow), (:Base, :/), (:Base, :+), (:Base, :-), (:Base, :sin), (:Base, :cos))
+    if (M, f) in ((:Base, :^), (:NaNMath, :pow), (:Base, :/), (:Base, :+), (:Base, :-), (:Base, :sin), (:Base, :cos), (:Base, :min), (:Base, :max))
         continue  # Skip methods which we define elsewhere.
     elseif !(isdefined(@__MODULE__, M) && isdefined(getfield(@__MODULE__, M), f))
         continue  # Skip rules for methods not defined in the current scope
@@ -543,6 +543,13 @@ end
     end
 )
 
+# min/max #
+#---------#
+
+# Consistent with `isless` (also on ties in value), and propagate NaN like `min`/`max` on floats; see issue #868
+@inline Base.min(x::D, y::D) where {D<:Dual} = ifelse(isnan(x) | (!isnan(y) & !isless(y, x)), x, y)
+@inline Base.max(x::D, y::D) where {D<:Dual} = ifelse(isnan(x) | (!isnan(y) & !isless(x, y)), x, y)
+
 # exponentiation #
 #----------------#
 

test/DualTest.jl

@@ -702,6 +702,57 @@ end
   @test floatmin(Dual{Nothing, ForwardDiff.Dual{Nothing, Float64, 2}, 1}) === Dual{Nothing}(Dual{Nothing}(floatmin(Float64),0.0,0.0),Dual{Nothing}(0.0,0.0,0.0))
 end
 
+@testset "min/max consistent with isless (issue #868)" begin
+    D(v, p...) = Dual{TestTag}(v, p...)
+    one_ = D(1.0, 0.0)
+    one_pos = D(1.0, 1.0)
+    one_neg = D(1.0, -1.0)
+    one_nan = D(1.0, NaN)
+    nan_ = D(NaN, 0.0)
+    nested_one = D(D(1.0, 0.0), D(0.0, 0.0))
+    nested_one_pos = D(D(1.0, 1.0), D(0.0, 0.0))
+    nested_one_innernan = D(D(1.0, NaN), D(0.0, 0.0))
+    nested_one_outernan = D(D(1.0, 0.0), D(NaN, 0.0))
+    nested_nan = D(D(NaN, 0.0), D(0.0, 0.0))
+
+    # (a, b, min(a, b), max(a, b))
+    cases = [
+        (one_pos, 1.0, one_, one_pos),
+        (one_pos, 1, one_, one_pos),
+        (one_neg, 1.0, one_neg, one_),
+        (one_pos, one_neg, one_neg, one_pos),
+        (one_pos, D(2.0, -1.0), one_pos, D(2.0, -1.0)),
+        (one_nan, one_, one_, one_nan),
+        (nan_, one_pos, nan_, nan_),
+        (nan_, 1.0, nan_, nan_),
+        (nested_one_pos, 1.0, nested_one, nested_one_pos),
+        (nested_one_innernan, nested_one, nested_one, nested_one_innernan),
+        (nested_one_outernan, nested_one, nested_one, nested_one_outernan),
+        (nested_nan, nested_one_pos, nested_nan, nested_nan),
+        (nested_nan, 1.0, nested_nan, nested_nan),
+    ]
+    for (a, b, lo, hi) in cases
+        @test min(a, b) === min(b, a) === lo
+        @test max(a, b) === max(b, a) === hi
+        if !isnan(a) && !isnan(b)
+            @test !isless(a, min(a, b)) && !isless(b, min(a, b))
+            @test !isless(max(a, b), a) && !isless(max(a, b), b)
+        end
+    end
+
+    @test ForwardDiff.derivative(t -> min(t, 48.0), 48.0) === ForwardDiff.derivative(t -> min(48.0, t), 48.0) === 0.0
+    @test ForwardDiff.derivative(t -> max(t, 48.0), 48.0) === ForwardDiff.derivative(t -> max(48.0, t), 48.0) === 1.0
+
+    d(f) = t -> ForwardDiff.derivative(f, t)
+    for f in (t -> min(3t - 2, t^2), t -> min(t^2, 3t - 2))
+        @test d(f)(1.0) === d(d(f))(1.0) === 2.0
+    end
+    for f in (t -> max(3t - 2, t^2), t -> max(t^2, 3t - 2))
+        @test d(f)(1.0) === 3.0
+        @test d(d(f))(1.0) === 0.0
+    end
+end
+
 @testset "Integer" begin
     x = Dual(1.0,0)
     @test Int(x) ≡ Integer(x) ≡ convert(Int,x) ≡ convert(Integer,x) ≡ 1