Central Finite Difference Method returns value of type Any

[rvignolo]

Hi!

I was wondering if the following is or is not an issue. This is a MWE:

julia> using FiniteDifferences
julia> const _fdm = FiniteDifferences.central_fdm(5, 1)
julia> @code_warntype _fdm(sin, 1)
Variables
  d::FiniteDifferences.Central{UnitRange{Int64},Array{Float64,1}}
  f::Core.Compiler.Const(sin, false)
  x::Int64

Body::Any
1 ─ %1 = Core.NamedTuple()::Core.Compiler.Const(NamedTuple(), false)
│   %2 = Base.pairs(%1)::Core.Compiler.Const(Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}(), false)
│   %3 = FiniteDifferences.:(var"#_#7")(%2, d, f, x)::Any
└──      return %3

This returns a value of type Any and sometimes I am afraid if this is potentially a performant killer for code where FiniteDifferences is used.

I would appreciate any insight you can give me!

Thanks!

[rvignolo]

I would also want to point out that the code above is type stable when using either FiniteDifferences.forward_fdm or FiniteDifferences.backward_fdm.

[rvignolo]

Hi!

Any ideas?

[willtebbutt]

Hiya, sorry, busy period for me at the minute. Not quite sure what's going on, but this is clearly an issue. Will try and take a look soon.

[rvignolo]

I don't like it, but the following two changes may give you some insight into the underlying problem:

FiniteDifferences.jl/src/methods.jl

Line 238 in f1cdb2a

m.q == 0 && return m, f(x)

to

    m.q == 0 && return m, float(f(x))

and

FiniteDifferences.jl/src/methods.jl

Line 270 in f1cdb2a

dfdx = sum(i->coefs[i] * f(x + ĥ * grid[i]), eachindex(grid)) / ĥ^q

to

    dfdx::Float64 = sum(i->coefs[i] * f(x + ĥ * grid[i]), eachindex(grid)) / ĥ^q

solve the problem.

I have encountered a different case using forward_fdm which also lead to a type instability:

julia> const _fdm2 = FiniteDifferences.forward_fdm(5, 1)
FiniteDifferenceMethod:
  order of method:       5
  order of derivative:   1
  grid:                  0:4
  coefficients:          [-2.0833333333333335, 4.0, -3.0, 1.3333333333333333, -0.25]
julia> @code_warntype _fdm2(t -> t, 1)
Variables
  d::FiniteDifferences.Forward{UnitRange{Int64},Array{Float64,1}}
  f::Core.Compiler.Const(var"#101#102"(), false)
  x::Int64

Body::Union{Float64, Int64, FiniteDifferences.Forward{UnitRange{Int64},Array{Float64,1}}}
1 ─ %1 = Core.NamedTuple()::Core.Compiler.Const(NamedTuple(), false)
│   %2 = Base.pairs(%1)::Core.Compiler.Const(Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}(), false)
│   %3 = FiniteDifferences.:(var"#_#5")(%2, d, f, x)::Union{Float64, Int64, FiniteDifferences.Forward{UnitRange{Int64},Array{Float64,1}}}
└──      return %3

The suggested changes fix the type instabilities but, as I previously said, I don't like those changes. I am just pointing this out to simplify your inspection. I hope this helps!

Thanks and regards!

[oxinabox]

I actually thing they are basically the right way to fix it.
the dfdx is going to be a Float64 because we have the global const TINY defined as a Float64.

I would not use float(f(x)) because I don't just want the floating point version of what ever f(x) returns, I really and actually want Float64

I was thinking of doing

function fdm(
    m::M,
    f,
    x,
    ::Val{true};
    condition=DEFAULT_CONDITION,
    bound=_estimate_bound(f(x), condition),
    eps=(Base.eps(float(bound)) + TINY),
    adapt=m.history.adapt,
    max_step=0.1,
)::Float64 where M<:FiniteDifferenceMethod

as a maybe nicer way to do this

[willtebbutt]

the dfdx is going to be a Float64 because we have the global const TINY defined as a Float64.

@wesselb and I were chatting about this. The TINY value should probably be type-dependent, and assuming it's a Float64 is currently damaging our accuracy. Not sure how this plays with this discussion as I've still not had the time to properly understand what's going on here.

[rvignolo]

function fdm(
m::M,
f,
x,
::Val{true};
condition=DEFAULT_CONDITION,
bound=_estimate_bound(f(x), condition),
eps=(Base.eps(float(bound)) + TINY),
adapt=m.history.adapt,
max_step=0.1,
)::Float64 where M<:FiniteDifferenceMethod

On the other hand, I am not sure we can implement it this way given that fdm returns the tuple (m, dfdx).

[oxinabox]

On the other hand, I am not sure we can implement it this way given that fdm returns the tuple (m, dfdx).

I think we can by using )::Tuple{Any, Float64} or something, but maybe not.
Any equiv way of writing will do anyway.

@wesselb and I were chatting about this. The TINY value should probably be type-dependent, and assuming it's a Float64 is currently damaging our accuracy. Not sure how this plays with this discussion as I've still not had the time to properly understand what's going on he

we should do a bit of a bigger redesign of things before then.
for now i would stick to Float64

[rvignolo]

I think we can by using )::Tuple{Any, Float64} or something, but maybe not.

well, that worked!

[wesselb]

I believe this is resolved with v0.12.0:

julia> using FiniteDifferences

julia> const _fdm = FiniteDifferences.central_fdm(5, 1)

julia> @code_warntype _fdm(sin, 1)
Variables
  m::FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}}
  f::Core.Const(sin)
  x@_3::Int64
  step::Float64
  x@_5::Union{Float64, Int64}

Body::Float64
1 ─      (x@_5 = x@_3)
│        (x@_5 = FiniteDifferences.float(x@_5::Int64))
│   %3 = FiniteDifferences.estimate_step(m, f, x@_5::Float64)::Tuple{Float64, Float64}
│        (step = FiniteDifferences.first(%3))
│   %5 = (m)(f, x@_5::Float64, step)::Float64
└──      return %5

Could anyone confirm this?

[sethaxen]

I can confirm this specific issue is resolved.

[rvignolo]

thank you!