NaNs inputs to function when using one-sided differences

[tpapp]

Adapted one-sided differences gives a NaN input while receiving finite values from the function.

(jl_yIAkgW) pkg> st
Status `/tmp/jl_yIAkgW/Project.toml`
  [26cc04aa] FiniteDifferences v0.12.32 `https://github.com/JuliaDiff/FiniteDifferences.jl.git#main`

julia> using FiniteDifferences

julia> function f(x::Real)
       A = L = 1.0
       x2 = abs2(x)
       @show x, x2
       @assert x2 ≤ 1
       if x2 == 1
       y = oftype(A + L + x, Inf)
       x < 0 ? -y : y
       else
       A + L * x / √(1 - x2)
       end
       end
f (generic function with 1 method)

julia> forward_fdm(5, 1)(f, -1.0)
(x, x2) = (-1.0, 1.0)
(x, x2) = (-0.9959053413840389, 0.9918274489972589)
(x, x2) = (-0.9918106827680776, 0.9836884304528802)
(x, x2) = (-0.9877160241521165, 0.9755829443668643)
(x, x2) = (-0.9836213655361552, 0.9675109907392107)
(x, x2) = (-0.9795267069201941, 0.9594725695699198)
(x, x2) = (-0.9754320483042329, 0.9514676808589914)
(x, x2) = (NaN, NaN)
ERROR: AssertionError: x2 ≤ 1
Stacktrace:
  [1] f(x::Float64)
    @ Main ./REPL[10]:5
  [2] newf
    @ ~/.julia/packages/StaticArrays/85pEu/src/broadcast.jl:186 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/StaticArrays/85pEu/src/broadcast.jl:135 [inlined]
  [4] __broadcast
    @ ~/.julia/packages/StaticArrays/85pEu/src/broadcast.jl:123 [inlined]
  [5] _broadcast
    @ ~/.julia/packages/StaticArrays/85pEu/src/broadcast.jl:119 [inlined]
  [6] copy
    @ ~/.julia/packages/StaticArrays/85pEu/src/broadcast.jl:60 [inlined]
  [7] materialize
    @ ./broadcast.jl:903 [inlined]
  [8] _eval_function(m::FiniteDifferences.AdaptedFiniteDifferenceMethod{…}, f::typeof(f), x::Float64, step::Float64)
    @ FiniteDifferences ~/.julia/packages/FiniteDifferences/IPGFN/src/methods.jl:249
  [9] AdaptedFiniteDifferenceMethod
    @ ~/.julia/packages/FiniteDifferences/IPGFN/src/methods.jl:240 [inlined]
 [10] (::FiniteDifferences.AdaptedFiniteDifferenceMethod{…})(f::typeof(f), x::Float64)
    @ FiniteDifferences ~/.julia/packages/FiniteDifferences/IPGFN/src/methods.jl:194
 [11] top-level scope
    @ REPL[11]:1
Some type information was truncated. Use `show(err)` to see complete types.

[wesselb]

@tpapp Might the issue be that the f(1) and f(-1) are not finite?

julia> function f(x::Real)
       A = L = 1.0
       x2 = abs2(x)
       @show x, x2
       @assert x2 ≤ 1
       if x2 == 1
       y = oftype(A + L + x, Inf)
       x < 0 ? -y : y
       else
       A + L * x / √(1 - x2)
       end
       end
f (generic function with 1 method)

julia> f(1)
(x, x2) = (1, 1)
Inf

julia> f(-1)
(x, x2) = (-1, 1)
-Inf

[tpapp]

Thanks, it is user error then. Is there a way to catch this in the finite differencing method though? Eg check that internal arguments remain isfinite.

[wesselb]

@tpapp It should be fairly simple to add a check whether the function evaluations are finite or not. I suppose that could help. :)