Type Conversion Error with Measurements

[lstagner]

using FiniteDifferences
using Measurements

f(x) = cos(x + (1.0 ± 0.1)*sin(x))
m = central_fdm(5,1)

m(f, 0.0)

MethodError: no method matching Float64(::Measurement{Float64})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:772
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...

Stacktrace:
 [1] _eval_function(m::FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}}, f::typeof(f), x::Float64, step::Measurement{Float64})
   @ FiniteDifferences ~/.julia/packages/FiniteDifferences/8zye5/src/methods.jl:249
 [2] (::FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}})(f::typeof(f), x::Float64, step::Measurement{Float64})
   @ FiniteDifferences ~/.julia/packages/FiniteDifferences/8zye5/src/methods.jl:240
 [3] (::FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}})(f::typeof(f), x::Float64)
   @ FiniteDifferences ~/.julia/packages/FiniteDifferences/8zye5/src/methods.jl:194
 [4] top-level scope
   @ In[6]:7
 [5] eval
   @ ./boot.jl:368 [inlined]
 [6] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1428

[oxinabox]

Quick workaround: turn-off adaption.

julia> m2 = central_fdm(5, 1; adapt=0);

julia> m2(f, 0.0)
0.0 ± 0.00014

---

There are a few places in the code,
namely

FiniteDifferences.jl/src/methods.jl

Lines 246 to 250 in 93cd547

	function _eval_function(
	m::FiniteDifferenceMethod, f::TF, x::T, step::Real,
	) where {TF<:Function,T<:AbstractFloat}
	return f.(x .+ T(step) .* m.grid)
	end

and

FiniteDifferences.jl/src/methods.jl

Line 263 in 93cd547

return sum(fs .* _coefs) ./ T(step)^Q

where we use T(step) to convert the step into the same type as the primal value x.
This was done when chasing down allocations from type-instabilities.
Mathematically, I don't think this is required.
We just need x+step to be the same type as x.
However, I don't thing step::Measurement{Float64} satisfies that: x + step gives a result that is a Measurement{Float64}.

But the question then is where is our step coming from?
That is coming out of FiniteDifferences.estimate_step

julia> FiniteDifferences.estimate_step(m, f, 0.0)
(0.00097 ± 2500.0, 4.3e-13 ± 1.1e-6)

Something is off somewhere in the adaption logic of

FiniteDifferences.jl/src/methods.jl

Lines 362 to 426 in 93cd547

	function estimate_step(
	m::AdaptedFiniteDifferenceMethod{P,Q}, f::TF, x::T,
	) where {P,Q,TF<:Function,T<:AbstractFloat}
	∇f_magnitude, f_magnitude = _estimate_magnitudes(m.bound_estimator, f, x)
	if ∇f_magnitude == 0.0 || f_magnitude == 0.0
	step, acc = _compute_step_acc_default(m, x)
	else
	step, acc = _compute_step_acc(m, ∇f_magnitude, eps(f_magnitude))
	end
	return _limit_step(m, x, step, acc)
	end
	function _estimate_magnitudes(
	m::FiniteDifferenceMethod{P,Q}, f::TF, x::T,
	) where {P,Q,TF<:Function,T<:AbstractFloat}
	step = first(estimate_step(m, f, x))
	fs = _eval_function(m, f, x, step)
	# Estimate magnitude of `∇f` in a neighbourhood of `x`.
	∇fs = SVector{3}(
	_compute_estimate(m, fs, x, step, m.coefs_neighbourhood[1]),
	_compute_estimate(m, fs, x, step, m.coefs_neighbourhood[2]),
	_compute_estimate(m, fs, x, step, m.coefs_neighbourhood[3])
	)
	∇f_magnitude = maximum(maximum.(abs, ∇fs))
	# Estimate magnitude of `f` in a neighbourhood of `x`.
	f_magnitude = maximum(maximum.(abs, fs))
	return ∇f_magnitude, f_magnitude
	end
	function _compute_step_acc_default(m::FiniteDifferenceMethod, x::T) where {T<:AbstractFloat}
	# Compute a default step size using a heuristic and [`DEFAULT_CONDITION`](@ref).
	return _compute_step_acc(m, m.condition, eps(T))
	end
	function _compute_step_acc(
	m::FiniteDifferenceMethod{P,Q}, ∇f_magnitude::Real, f_error::Real,
	) where {P,Q}
	# Set the step size by minimising an upper bound on the error of the estimate.
	C₁ = f_error * m.f_error_mult * m.factor
	C₂ = ∇f_magnitude * m.∇f_magnitude_mult
	step = (Q / (P - Q) * (C₁ / C₂))^(1 / P)
	# Estimate the accuracy of the method.
	acc = C₁ * step^(-Q) + C₂ * step^(P - Q)
	return step, acc
	end
	function _limit_step(
	m::FiniteDifferenceMethod, x::T, step::Real, acc::Real,
	) where {T<:AbstractFloat}
	# First, limit the step size based on the maximum range.
	step_max = m.max_range / maximum(abs.(m.grid))
	if step > step_max
	step = step_max
	acc = NaN
	end
	# Second, prevent very large step sizes, which can occur for high-order methods or
	# slowly-varying functions.
	step_default, _ = _compute_step_acc_default(m, x)
	step_max_default = 1000step_default
	if step > step_max_default
	step = step_max_default
	acc = NaN
	end
	return step, acc
	end

because it is yielding a step that is of the same type as f(x) rather than the type of x.
@wesselb might have ideas
but otherwise this is going to take a fair bit of digging to uncover.