don't propagate NaN/Inf if Partials are zero

Author: jrevels

.travis.yml

@@ -10,10 +10,10 @@ script:
     - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
     - if (julia -e 'VERSION < v"0.5" && exit(1)'); then
         julia -e 'include(joinpath(JULIA_HOME, Base.DATAROOTDIR, "julia", "build_sysimg.jl")); build_sysimg(force=true)';
-        julia -e 'Pkg.clone(pwd()); Pkg.build("ForwardDiff"); Pkg.test("ForwardDiff"; coverage=true)';
+        julia -e 'Pkg.clone(pwd()); Pkg.build("ForwardDiff"); Pkg.checkout("DiffBase"); Pkg.test("ForwardDiff"; coverage=true)';
         julia -O3 -e 'include(joinpath(Pkg.dir("ForwardDiff"), "test/SIMDTest.jl"))';
       else
-        julia -e 'Pkg.clone(pwd()); Pkg.build("ForwardDiff"); Pkg.test("ForwardDiff"; coverage=true)';
+        julia -e 'Pkg.clone(pwd()); Pkg.build("ForwardDiff"); Pkg.checkout("DiffBase"); Pkg.test("ForwardDiff"; coverage=true)';
       fi
 after_success:
     - julia -e 'cd(Pkg.dir("ForwardDiff")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'

src/partials.jl

@@ -75,18 +75,29 @@ Base.convert{N,T}(::Type{Partials{N,T}}, partials::Partials{N,T}) = partials
 @inline @compat(Base.:+){N}(a::Partials{N}, b::Partials{N}) = Partials(add_tuples(a.values, b.values))
 @inline @compat(Base.:-){N}(a::Partials{N}, b::Partials{N}) = Partials(sub_tuples(a.values, b.values))
 @inline @compat(Base.:-)(partials::Partials) = Partials(minus_tuple(partials.values))
-@inline @compat(Base.:*)(partials::Partials, x::Real) = Partials(scale_tuple(partials.values, x))
 @inline @compat(Base.:*)(x::Real, partials::Partials) = partials*x
-@inline @compat(Base.:/)(partials::Partials, x::Real) = Partials(div_tuple_by_scalar(partials.values, x))
 
-@inline function _mul_partials{N}(a::Partials{N}, b::Partials{N}, afactor, bfactor)
-    return Partials(mul_tuples(a.values, b.values, afactor, bfactor))
+# NaN/Inf-safe methods #
+#----------------------#
+
+@inline function @compat(Base.:*)(partials::Partials, x::Real)
+    x = ifelse((isnan(x) || isinf(x)) && iszero(partials), one(x), x)
+    return Partials(scale_tuple(partials.values, x))
+end
+
+@inline function @compat(Base.:/)(partials::Partials, x::Real)
+    x = ifelse(x == zero(x) && iszero(partials), one(x), x)
+    return Partials(div_tuple_by_scalar(partials.values, x))
+end
+
+@inline function _mul_partials{N}(a::Partials{N}, b::Partials{N}, x_a, x_b)
+    x_a = ifelse((isnan(x_a) || isinf(x_a)) && iszero(a), one(x_a), x_a)
+    x_b = ifelse((isnan(x_b) || isinf(x_b)) && iszero(b), one(x_b), x_b)
+    return Partials(mul_tuples(a.values, b.values, x_a, x_b))
 end
 
 @inline function _div_partials(a::Partials, b::Partials, aval, bval)
-    afactor = inv(bval)
-    bfactor = -aval/(bval*bval)
-    return _mul_partials(a, b, afactor, bfactor)
+    return _mul_partials(a, b, inv(bval), -(aval / (bval*bval)))
 end
 
 # edge cases where N == 0 #

test/JacobianTest.jl

@@ -99,11 +99,11 @@ for f in DiffBase.ARRAY_TO_ARRAY_FUNCS
         out = ForwardDiff.jacobian(f, X, cfg)
         @test_approx_eq out j
 
-        out = similar(X, length(X), length(X))
+        out = similar(X, length(v), length(X))
         ForwardDiff.jacobian!(out, f, X, cfg)
         @test_approx_eq out j
 
-        out = DiffBase.JacobianResult(X)
+        out = DiffBase.DiffResult(similar(v, length(v)), similar(v, length(v), length(X)))
         ForwardDiff.jacobian!(out, f, X, cfg)
         @test_approx_eq DiffBase.value(out) v
         @test_approx_eq DiffBase.jacobian(out) j

test/PartialsTest.jl

@@ -9,10 +9,10 @@ samerng() = MersenneTwister(1)
 for N in (0, 3), T in (Int, Float32, Float64)
     println("  ...testing Partials{$N,$T}")
 
-    VALUES = ntuple(n -> rand(T), Val{N})
+    VALUES = (rand(T,N)...)
     PARTIALS = Partials{N,T}(VALUES)
 
-    VALUES2 = ntuple(n -> rand(T), Val{N})
+    VALUES2 = (rand(T,N)...)
     PARTIALS2 = Partials{N,T}(VALUES2)
 
     ##############################
@@ -70,7 +70,7 @@ for N in (0, 3), T in (Int, Float32, Float64)
     @test hash(PARTIALS, hash(1)) == hash(copy(PARTIALS), hash(1))
     @test hash(PARTIALS, hash(1)) == hash(copy(PARTIALS), hash(1))
 
-    const TMPIO = IOBuffer()
+    TMPIO = IOBuffer()
     write(TMPIO, PARTIALS)
     seekstart(TMPIO)
     @test read(TMPIO, typeof(PARTIALS)) == PARTIALS
@@ -84,8 +84,8 @@ for N in (0, 3), T in (Int, Float32, Float64)
     # Conversion/Promotion #
     ########################
 
-    const WIDE_T = widen(T)
-    const WIDE_PARTIALS = convert(Partials{N,WIDE_T}, PARTIALS)
+    WIDE_T = widen(T)
+    WIDE_PARTIALS = convert(Partials{N,WIDE_T}, PARTIALS)
 
     @test typeof(WIDE_PARTIALS) == Partials{N,WIDE_T}
     @test WIDE_PARTIALS == PARTIALS
@@ -101,8 +101,8 @@ for N in (0, 3), T in (Int, Float32, Float64)
     @test (PARTIALS - PARTIALS).values == map(v -> v - v, VALUES)
     @test getfield(-(PARTIALS), :values) == map(-, VALUES)
 
-    const X = rand()
-    const Y = rand()
+    X = rand()
+    Y = rand()
 
     @test X * PARTIALS == PARTIALS * X
     @test (X * PARTIALS).values == map(v -> X * v, VALUES)
@@ -111,6 +111,19 @@ for N in (0, 3), T in (Int, Float32, Float64)
     if N > 0
         @test ForwardDiff._mul_partials(PARTIALS, PARTIALS2, X, Y).values == map((a, b) -> (X * a) + (Y * b), VALUES, VALUES2)
         @test ForwardDiff._div_partials(PARTIALS, PARTIALS2, X, Y) == ForwardDiff._mul_partials(PARTIALS, PARTIALS2, inv(Y), -X/(Y^2))
+
+        ZEROS = Partials((zeros(T, N)...))
+
+        @test (NaN * ZEROS).values == ZEROS.values
+        @test (Inf * ZEROS).values == ZEROS.values
+        @test (ZEROS / 0).values == ZEROS.values
+
+        @test ForwardDiff._mul_partials(ZEROS, ZEROS, X, NaN).values == ZEROS.values
+        @test ForwardDiff._mul_partials(ZEROS, ZEROS, NaN, X).values == ZEROS.values
+        @test ForwardDiff._mul_partials(ZEROS, ZEROS, X, Inf).values == ZEROS.values
+        @test ForwardDiff._mul_partials(ZEROS, ZEROS, Inf, X).values == ZEROS.values
+        @test ForwardDiff._mul_partials(ZEROS, ZEROS, Inf, NaN).values == ZEROS.values
+        @test ForwardDiff._mul_partials(ZEROS, ZEROS, NaN, Inf).values == ZEROS.values
     end
 end