Merge pull request #7 from mlubin/nanmath

use NaNMath package to return NaNs instead of throwing DomainErrors

Author: mlubin

REQUIRE

@@ -1,3 +1,4 @@
 julia 0.3
 DualNumbers
+NaNMath
 Graphs

src/ReverseDiffSparse.jl

@@ -3,6 +3,12 @@ module ReverseDiffSparse
 import Calculus
 using DualNumbers
 using Base.Meta
+# Override basic math functions to return NaN instead of throwing errors.
+# This is what NLP solvers expect, and
+# sometimes the results aren't needed anyway,
+# because the code may compute derivatives wrt constants.
+import NaNMath: sin, cos, tan, asin, acos, acosh, atanh, log, log2, log10, lgamma, log1p, pow
+
 if isdir(Pkg.dir("ArrayViews"))
     eval(Expr(:import,:ArrayViews))
     const subarr = ArrayViews.view

src/revmode.jl

@@ -244,7 +244,11 @@ function forwardpass(x::ExprNode, expr_out)
         for i in 2:length(x.ex.args)
             push!(values, forwardpass(x.ex.args[i], expr_out))
         end
-        fcall = Expr(:call, x.ex.args[1], values...)
+        if x.ex.args[1] == :(^) # Use NaNMath.pow instead of ^
+            fcall = Expr(:call, :pow, values...)
+        else
+            fcall = Expr(:call, x.ex.args[1], values...)
+        end
         push!(expr_out.args, :( $(x.value) = $fcall ))
         return x.value
     elseif isexpr(x.ex, :curly)
@@ -284,11 +288,6 @@ forwardpass(x, expr_out) = :(forwardvalue($x, __placevalues, __placeindex_in))
 forwardvalue(x::Placeholder, placevalues, placeindex_in) = placevalues[placeindex_in[getplaceindex(x)]]
 forwardvalue(x, placevalues, placeindex_in) = float(x)
 
-# better to return NaNs than throw DomainErrors.
-# sometimes the results aren't needed anyway,
-# because the code may compute derivatives wrt constants.
-log(x) = x <= 0 ? NaN : Base.log(x)
-
 function revpass(x::ExprNode, expr_out)
     @assert isexpr(expr_out, :block)
     # compute the partial drivative wrt. each expression down the graph
@@ -342,12 +341,12 @@ function revpass(x::ExprNode, expr_out)
             if k == 2 # base
                 exponent = getvalue(p.ex.args[3])
                 push!(expr_out.args, 
-                  :( $(x.deriv) += $(p.deriv)*$exponent*$(x.value)^($exponent-1) ))
+                  :( $(x.deriv) += $(p.deriv)*$exponent*pow($(x.value),$exponent-1) ))
             else
                 @assert k == 3
                 base = getvalue(p.ex.args[2])
                 push!(expr_out.args,
-                  :( $(x.deriv) += $(p.deriv)*$base^($(x.value))*log($base) ))
+                  :( $(x.deriv) += $(p.deriv)*pow($base,$(x.value))*log($base) ))
             end
         elseif f == :(/)
             if k == 2 # numerator
@@ -358,7 +357,7 @@ function revpass(x::ExprNode, expr_out)
                 @assert k == 3 # denominator
                 numer = getvalue(p.ex.args[2])
                 push!(expr_out.args,
-                  :( $(x.deriv) += -1*$(p.deriv)*$numer*($(x.value))^(-2) ))
+                  :( $(x.deriv) += -1*$(p.deriv)*$numer*pow($(x.value),-2) ))
             end
         else
             # try one of the derivative rules

test/test_grad.jl

@@ -205,6 +205,13 @@ fval = fg([-3.0],out)
 @test_approx_eq fval (-3)^2
 @test_approx_eq out[1] 2*-3
 
+y = -2
+ex = @processNLExpr y^x[1]
+fg = genfgrad_simple(ex)
+fval = fg([0.3],out)
+@test isnan(fval)
+
+
 # zeros in products
 ex = @processNLExpr prod{ x[i], i = 1:2 }
 fg = genfgrad_simple(ex)