Correctness issues + annoying behavior

[lrnv]

Hey,

The following behaviors are really annoying:

julia> t = TaylorScalar{Float64, 3}((Inf, 1.0, 0.0))     
TaylorScalar{Float64, 3}((Inf, 1.0, 0.0))

julia> t+1
TaylorScalar{Float64, 3}((Inf, 1.0, 0.0)) ##### OK

julia> t*1
TaylorScalar{Float64, 3}((Inf, NaN, NaN)) ##### Should be Inf, 1.0, 0.0

julia> t*2
TaylorScalar{Float64, 3}((Inf, NaN, NaN)) ##### Should be Inf, 2.0, 0.0

julia> 

The last one even looks like a correctness issue...

but also :

julia> t2 = TaylorScalar{Float64, 3}((0.0, 1.0, 0.0))    
TaylorScalar{Float64, 3}((0.0, 1.0, 0.0))

julia> t2^1
TaylorScalar{Float64, 3}((0.0, NaN, NaN)) ##### Should be 0.0, 1.0, 0.0

julia> t2^2
TaylorScalar{Float64, 3}((0.0, NaN, NaN)) ##### Should not be NaNs...

julia> 

Using TaylorSeries.jl dos not produce these NaNs, and this makes my test cases fail... Do you think this is something that would be fixable in the current state of TaylorDiff ?

[tansongchen]

Thanks for noticing that! I will take a closer look later this week

[lrnv]

So I wrote the following test:

all_equal(a,b) = all(a.value .== b.value)
offenders = (
    TaylorDiff.TaylorScalar{Float64, 4}((Inf, 1.0, 0.0, 0.0)),
    TaylorDiff.TaylorScalar{Float64, 4}((Inf, 0.0, 0.0, 0.0)),
    TaylorDiff.TaylorScalar{Float64, 4}((1.0, 0.0, 0.0, 0.0)),
    TaylorDiff.TaylorScalar{Float64, 4}((1.0, Inf, 0.0, 0.0)),
    TaylorDiff.TaylorScalar{Float64, 4}((0.0, 1.0, 0.0, 0.0)),
    TaylorDiff.TaylorScalar{Float64, 4}((0.0, Inf, 0.0, 0.0))
    # Others ?
)
f_id = (
    :id => x -> x,
    :add0 => x -> x+0,
    :sub0 => x -> x-0,
    :mul1 => x -> x*1,
    :div1 => x -> x/1,
    :pow1 => x -> x^1,
)
for (name_f,f) in f_id, t in offenders
    if !all_equal(f(t),t)
        @show name_f, t, f(t)
    end
end

On the main branch, I have 11 faillures as follows:

(name_f, t, f(t)) = (:mul1, TaylorScalar{Float64, 4}((Inf, 1.0, 0.0, 0.0)), TaylorScalar{Float64, 4}((Inf, NaN, NaN, NaN)))
(name_f, t, f(t)) = (:mul1, TaylorScalar{Float64, 4}((Inf, 0.0, 0.0, 0.0)), TaylorScalar{Float64, 4}((Inf, NaN, NaN, NaN)))
(name_f, t, f(t)) = (:mul1, TaylorScalar{Float64, 4}((1.0, Inf, 0.0, 0.0)), TaylorScalar{Float64, 4}((1.0, Inf, NaN, NaN)))
(name_f, t, f(t)) = (:mul1, TaylorScalar{Float64, 4}((0.0, Inf, 0.0, 0.0)), TaylorScalar{Float64, 4}((0.0, Inf, NaN, NaN)))
(name_f, t, f(t)) = (:div1, TaylorScalar{Float64, 4}((Inf, 1.0, 0.0, 0.0)), TaylorScalar{Float64, 4}((Inf, NaN, NaN, NaN)))
(name_f, t, f(t)) = (:div1, TaylorScalar{Float64, 4}((Inf, 0.0, 0.0, 0.0)), TaylorScalar{Float64, 4}((Inf, NaN, NaN, NaN)))
(name_f, t, f(t)) = (:div1, TaylorScalar{Float64, 4}((1.0, Inf, 0.0, 0.0)), TaylorScalar{Float64, 4}((1.0, Inf, NaN, NaN)))
(name_f, t, f(t)) = (:div1, TaylorScalar{Float64, 4}((0.0, Inf, 0.0, 0.0)), TaylorScalar{Float64, 4}((0.0, Inf, NaN, NaN)))
(name_f, t, f(t)) = (:pow1, TaylorScalar{Float64, 4}((1.0, Inf, 0.0, 0.0)), TaylorScalar{Float64, 4}((1.0, Inf, NaN, NaN)))
(name_f, t, f(t)) = (:pow1, TaylorScalar{Float64, 4}((0.0, 1.0, 0.0, 0.0)), TaylorScalar{Float64, 4}((0.0, 1.0, NaN, NaN)))
(name_f, t, f(t)) = (:pow1, TaylorScalar{Float64, 4}((0.0, Inf, 0.0, 0.0)), TaylorScalar{Float64, 4}((0.0, Inf, NaN, NaN)))

By the way, I checked on my old PR #25, and there only the power function has issue:

(name_f, t, f(t)) = (:pow1, TaylorScalar{Float64, 4}((Inf, 1.0, 0.0, 0.0)), TaylorScalar{Float64, 4}((Inf, NaN, NaN, NaN)))
(name_f, t, f(t)) = (:pow1, TaylorScalar{Float64, 4}((Inf, 0.0, 0.0, 0.0)), TaylorScalar{Float64, 4}((Inf, NaN, NaN, NaN)))
(name_f, t, f(t)) = (:pow1, TaylorScalar{Float64, 4}((1.0, Inf, 0.0, 0.0)), TaylorScalar{Float64, 4}((1.0, Inf, NaN, NaN)))
(name_f, t, f(t)) = (:pow1, TaylorScalar{Float64, 4}((0.0, 1.0, 0.0, 0.0)), TaylorScalar{Float64, 4}((0.0, NaN, NaN, NaN)))
(name_f, t, f(t)) = (:pow1, TaylorScalar{Float64, 4}((0.0, Inf, 0.0, 0.0)), TaylorScalar{Float64, 4}((0.0, NaN, NaN, NaN)))

I do not recall what was the problem with this PR, why didnt it got merged at the time ? Maybe you could reconsider ?

Version of the code as a `@testset` to directly include it in your test suite

@testset "Correct limiting behaviors" begin
    # Test included related to the discussion in https://github.com/JuliaDiff/TaylorDiff.jl/issues/61
    # In particular, no nans should be produced ! 
    all_equal(a,b) = all(a.value .== b.value)
    offenders = (
        TaylorDiff.TaylorScalar{Float64, 4}((Inf, 1.0, 0.0, 0.0)),
        TaylorDiff.TaylorScalar{Float64, 4}((Inf, 0.0, 0.0, 0.0)),
        TaylorDiff.TaylorScalar{Float64, 4}((1.0, 0.0, 0.0, 0.0)),
        TaylorDiff.TaylorScalar{Float64, 4}((1.0, Inf, 0.0, 0.0)),
        TaylorDiff.TaylorScalar{Float64, 4}((0.0, 1.0, 0.0, 0.0)),
        TaylorDiff.TaylorScalar{Float64, 4}((0.0, Inf, 0.0, 0.0))
        # Others ?
    )
    f_id = (
        :id => x -> x,
        :add0 => x -> x+0,
        :sub0 => x -> x-0,
        :mul1 => x -> x*1,
        :div1 => x -> x/1,
        :pow1 => x -> x^1,
    )
    for (name_f,f) in f_id, t in offenders
        @test all_equal(f(t),t)
    end
end

[tansongchen]

These tests are added to package tests in c90fc69, and they passed

[lrnv]

@tansongchen Waouh that i really nice thanks a lot. I will come back If i need more as I just said in another message.

[tansongchen]

oh I just realized I misunderstood == and your allequal... in fact there's still 3 cases not passing. I need to read the source code of TaylorSeries.jl since they are really comprehensive at this

[lrnv]

Indeed the test you added in c90fc69 uses == which is not enough here as it only checks the first value.

Also the test added does only check identities, but we also have issues with other things like

2 * TaylorScalar{Float64, 3}((Inf, 1.0, 0.0))
TaylorScalar{Float64, 3}((0.0, 1.0, 0.0))^2

which are still producing NaNs for the moment. The product one should work with the new product with scalar implementation, but the square does not yet.

[tansongchen]

Update: I added specializations to Base.literal_pow which lowers x ^ 1 to x and x ^ 2 to x * x etc., which partially alleviates your problem.

I also introduced a new macro @to_static which can define nonlinear functions on polynomials without any metaprogramming tricks.

However, when I look at the source code of TaylorSeries.jl, I couldn't exactly understand how they implemented special case handling for the case where zeroth coefficient is 0. Therefore I couldn't translate to equivalent stuff at TaylorDiff.jl.

I raised an issue there JuliaDiff/TaylorSeries.jl#370 , hopefully they will explan :)

[lrnv]

Hope you'll find the right thing to do. Keep in mind that the examples I proposed here are only examples of a bigger, not fully understood by myself, issue. Might indeed be related to these elementary function shortcuts that TaylorSeries is using, but might be something deeper :/