Skip to content

Modification of equality checking between two Duals not synchronized between master and last release #726

Closed
@clguillot

Description

@clguillot

Hi everyone,

It seems that the modifications required by issue #481, which should have been in the code starting from release 0.10.33 (according to the release notes) are indeed included in the branch master, but are not present in the last release 0.10.38.
In particular, the following

ForwardDiff.jl/src/dual.jl

Lines 392 to 401 in 0a35a80

for pred in BINARY_PREDICATES
@eval begin
@define_binary_dual_op(
Base.$(pred),
$(pred)(value(x), value(y)),
$(pred)(value(x), y),
$(pred)(x, value(y))
)
end
end

should have been replaced for the operator == by

ForwardDiff.jl/src/dual.jl

Lines 391 to 407 in 3acc8a6

for pred in [:<, :>]
predeq = Symbol(pred, :(=))
@eval begin
@define_binary_dual_op(
Base.$(pred),
$(pred)(value(x), value(y)) || (value(x) == value(y) && $(pred)(partials(x), partials(y))),
$(pred)(value(x), y) || (value(x) == y && $(pred)(partials(x), zero(partials(x)))),
$(pred)(x, value(y)) || (x == value(y) && $(pred)(zero(partials(y)), partials(y))),
)
@define_binary_dual_op(
Base.$(predeq),
$(pred)(value(x), value(y)) || (value(x) == value(y) && $(predeq)(partials(x), partials(y))),
$(pred)(value(x), y) || (value(x) == y && $(predeq)(partials(x), zero(partials(x)))),
$(pred)(x, value(y)) || (x == value(y) && $(predeq)(zero(partials(y)), partials(y))),
)
end
end

As a consequence of this modification, equality between 2 Dual with the same tags should check both the equality of the value and the partials, hence avoiding some branching that would have been taken otherwise.
In particular, without this fix, we obtain

julia> f(x) = exp(complex(0, x))
f (generic function with 1 method)

julia> df(x) = ForwardDiff.derivative(f, x)
df (generic function with 1 method)

julia> d2f(x) = ForwardDiff.derivative(df, x)
d2f (generic function with 1 method)

julia> d2f(0.0)
0.0 + 0.0im

which is obviously wrong, but

julia> d2f(1e-15)
-1.0 - 1.0e-15im

which is right.
This can be explained by looking at the code of exp in complex.jl

function exp(z::Complex)
    zr, zi = reim(z)
    if isnan(zr)
        Complex(zr, zi==0 ? zi : zr)
    elseif !isfinite(zi)
        if zr == Inf
            Complex(-zr, oftype(zr,NaN))
        elseif zr == -Inf
            Complex(-zero(zr), copysign(zero(zi), zi))
        else
            Complex(oftype(zr,NaN), oftype(zi,NaN))
        end
    else
        er = exp(zr)
        if iszero(zi)
            Complex(er, zi)
        else
            s, c = sincos(zi)
            Complex(er * c, er * s)
        end
    end
end

Here, the branching if iszero(zi) is taken because iszero will only check that value(zi) is zero when it should also check that partials(zi) is zero (which won't be the case in this example).

Is there a reason why the modification that should have been included from release 0.10.33 onward can be found in the master branch but not in the last release ?

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions