Accept Unitful numbers as input

[SeSodesa]

Currently, when attempting to use unitful numbers from Unitful.jl with the functions in this library, one gets MethodErrors related to the fact that Unitful.Quantity{...} <: Number instead of Unitful.Quantity{...} <: Real:

ERROR: MethodError: no method matching (::FiniteDifferences.AdaptedFiniteDifferenceMethod{…})(::var"#35#36", ::Quantity{…})
The object of type `FiniteDifferences.AdaptedFiniteDifferenceMethod{3, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{5, 3}}` exists, but no method is defined for this combination of argument types when trying to treat it as a callable object.

Closest candidates are:
  (::FiniteDifferences.AdaptedFiniteDifferenceMethod)(::TF, ::Real) where TF
   @ FiniteDifferences ~/.julia/packages/FiniteDifferences/IPGFN/src/methods.jl:191
  (::FiniteDifferences.AdaptedFiniteDifferenceMethod{P, Q})(::TF, ::Real, ::Real) where {P, Q, TF}
   @ FiniteDifferences ~/.julia/packages/FiniteDifferences/IPGFN/src/methods.jl:238

At this point I am not sure whether simply relaxing the argument type restriction of functions like FiniteDifferences.central_fdm from Real to Number would suffice (with the possible risk of introducing silent semantic errors), or whether separate methods, that explicitly handle the units in some way (e.g. first stripping them and then reattaching the proper ones to the results) should be defined for Unitful numbers.

Regardless, it would be very nice if this library supported unitful numbers out of the box.

[SeSodesa]

Yep, changing all ::Real and ::AbstractFloat restrictions to ::Number in methods.jl does not come without issues:

import Unitful: @u_str
import Interpolations as Ilib
import FiniteDifferences as FDlib

tt = (0.0 : 0.01 : 1.0) .* u"s"

fn = 3u"J / s" .* tt

ii = Ilib.interpolate((tt,),fn, ComposedFunction(Ilib.Gridded, Ilib.Linear)())

cdm = FDlib.central_fdm(5,1)

dd = cdm(ii,1u"s")

This gives the error

ERROR: LoadError: BoundsError: attempt to access 101-element interpolate(((0.0:0.01:1.0) s,), ::StepRangeLen{Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}}}, Int64}, Gridded(Linear())) with element type Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}} at index [1.0072319884886336 s]
Stacktrace:
  [1] throw_boundserror(A::Interpolations.GriddedInterpolation{…}, I::Tuple{…})
    @ Base ./essentials.jl:14
  [2] GriddedInterpolation
    @ ~/.julia/packages/Interpolations/dR5oF/src/gridded/indexing.jl:3 [inlined]
  [3] #11
    @ ./broadcast.jl:353 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/StaticArrays/MTYJP/src/broadcast.jl:135 [inlined]
  [5] __broadcast
    @ ~/.julia/packages/StaticArrays/MTYJP/src/broadcast.jl:123 [inlined]
  [6] _broadcast
    @ ~/.julia/packages/StaticArrays/MTYJP/src/broadcast.jl:119 [inlined]
  [7] copy
    @ ~/.julia/packages/StaticArrays/MTYJP/src/broadcast.jl:60 [inlined]
  [8] materialize
    @ ./broadcast.jl:872 [inlined]
  [9] _eval_function(m::FiniteDifferences.UnadaptedFiniteDifferenceMethod{…}, f::Interpolations.GriddedInterpolation{…}, x::Unitful.Quantity{…}, step::Float64)
    @ FiniteDifferences ~/FiniteDifferences.jl/src/methods.jl:249
 [10] _estimate_magnitudes(m::FiniteDifferences.UnadaptedFiniteDifferenceMethod{…}, f::Interpolations.GriddedInterpolation{…}, x::Unitful.Quantity{…})
    @ FiniteDifferences ~/FiniteDifferences.jl/src/methods.jl:378
 [11] estimate_step(m::FiniteDifferences.AdaptedFiniteDifferenceMethod{…}, f::Interpolations.GriddedInterpolation{…}, x::Unitful.Quantity{…})
    @ FiniteDifferences ~/FiniteDifferences.jl/src/methods.jl:365
 [12] (::FiniteDifferences.AdaptedFiniteDifferenceMethod{…})(f::Interpolations.GriddedInterpolation{…}, x::Unitful.Quantity{…})
    @ FiniteDifferences ~/FiniteDifferences.jl/src/methods.jl:193
 [13] top-level scope
    @ ~/AneurysmPaper2025.jl/scratchpad.jl:19
 [14] include(fname::String)
    @ Main ./sysimg.jl:38
 [15] top-level scope
    @ REPL[13]:1

Although this seems like it might just be an issue with interpolation bounds. Let me try another case.

[SeSodesa]

Yeah, it was due to the interpolation without extrapolation. The code

import Unitful: @u_str
import Interpolations as Ilib
import FiniteDifferences as FDlib

tt = (0.0 : 0.01 : 1.1) .* u"s"

fn = 3u"J / s" .* tt

ii = Ilib.interpolate((tt,),fn, ComposedFunction(Ilib.Gridded, Ilib.Linear)())

cdm = FDlib.central_fdm(5,1)

dd = cdm(ii,1u"s")

produces the error

ERROR: LoadError: DimensionError: 1.2873703414413103e-13 J s⁻³⁶⁰²⁸⁷⁹⁷⁰¹⁸⁹⁶³⁹⁷ᐟ⁴⁵⁰³⁵⁹⁹⁶²⁷³⁷⁰⁴⁹⁶ and 3.21842585360327e-14 J s⁻⁹⁰⁰⁷¹⁹⁹²⁵⁴⁷⁴⁰⁹⁹ᐟ¹¹²⁵⁸⁹⁹⁹⁰⁶⁸⁴²⁶²⁴ are not dimensionally compatible.
Stacktrace:
 [1] +(x::Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻¹²⁶¹⁰⁰⁷⁸⁹⁵⁶⁶³⁷³⁸⁹ᐟ⁴⁵⁰³⁵⁹⁹⁶²⁷³⁷⁰⁴⁹⁶, Unitful.FreeUnits{(J, s⁻³⁶⁰²⁸⁷⁹⁷⁰¹⁸⁹⁶³⁹⁷ᐟ⁴⁵⁰³⁵⁹⁹⁶²⁷³⁷⁰⁴⁹⁶), 𝐋² 𝐌 𝐓⁻¹²⁶¹⁰⁰⁷⁸⁹⁵⁶⁶³⁷³⁸⁹ᐟ⁴⁵⁰³⁵⁹⁹⁶²⁷³⁷⁰⁴⁹⁶, nothing}}, y::Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻³¹⁵²⁵¹⁹⁷³⁹¹⁵⁹³⁴⁷ᐟ¹¹²⁵⁸⁹⁹⁹⁰⁶⁸⁴²⁶²⁴, Unitful.FreeUnits{(J, s⁻⁹⁰⁰⁷¹⁹⁹²⁵⁴⁷⁴⁰⁹⁹ᐟ¹¹²⁵⁸⁹⁹⁹⁰⁶⁸⁴²⁶²⁴), 𝐋² 𝐌 𝐓⁻³¹⁵²⁵¹⁹⁷³⁹¹⁵⁹³⁴⁷ᐟ¹¹²⁵⁸⁹⁹⁹⁰⁶⁸⁴²⁶²⁴, nothing}})
   @ Unitful ~/.julia/packages/Unitful/sxiHA/src/quantities.jl:137
 [2] _compute_step_acc(m::FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}}, ∇f_magnitude::Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻⁶, Unitful.FreeUnits{(J, s⁻⁴), 𝐋² 𝐌 𝐓⁻⁶, nothing}}, f_error::Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}})
   @ FiniteDifferences ~/FiniteDifferences.jl/src/methods.jl:404
 [3] estimate_step(m::FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}}, f::Interpolations.GriddedInterpolation{Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}}, 1, StepRangeLen{Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}}}, Int64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, Int64}}}, x::Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}})
   @ FiniteDifferences ~/FiniteDifferences.jl/src/methods.jl:369
 [4] (::FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}})(f::Interpolations.GriddedInterpolation{Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}}, 1, StepRangeLen{Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋² 𝐌 𝐓⁻², Unitful.FreeUnits{(J,), 𝐋² 𝐌 𝐓⁻², nothing}}}, Int64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, Int64}}}, x::Unitful.Quantity{Int64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}})
   @ FiniteDifferences ~/FiniteDifferences.jl/src/methods.jl:193
 [5] top-level scope

So now we get funky dimension mismatches related to powers of units being almost the same, but not quite. Maybe a more correct approach is to just strip units for now.

[SeSodesa]

Ok, so _compute_step_acc seems to be behaving strangely:

m.f_error_mult = 10.0
m.factor = 1.0
∇f_magnitude = 10.0
C₁ = f_error * m.f_error_mult * m.factor = 2.220446049250313e-15
C₂ = ∇f_magnitude * m.∇f_magnitude_mult = 5.365079365079365
P = 7
Q = 5
step = ((Q / (P - Q)) * (C₁ / C₂)) ^ (1 / P) = 0.007231988488633704
acc = C₁ * step ^ -Q + C₂ * step ^ (P - Q) = 0.00039284356077575253
m.f_error_mult = 10.0
m.factor = 1.0
∇f_magnitude = 10.0
C₁ = f_error * m.f_error_mult * m.factor = 2.220446049250313e-15
C₂ = ∇f_magnitude * m.∇f_magnitude_mult = 5.365079365079365
P = 7
Q = 5
step = ((Q / (P - Q)) * (C₁ / C₂)) ^ (1 / P) = 0.007231988488633704
acc = C₁ * step ^ -Q + C₂ * step ^ (P - Q) = 0.00039284356077575253
m.f_error_mult = 1.4999999999999998
m.factor = 1.0
∇f_magnitude = 0.0008081353250244043 J s⁻⁴
C₁ = f_error * m.f_error_mult * m.factor = 6.661338147750938e-16 J
C₂ = ∇f_magnitude * m.∇f_magnitude_mult = 4.489640694580024e-5 J s⁻⁴
P = 5
Q = 1
step = ((Q / (P - Q)) * (C₁ / C₂)) ^ (1 / P) = 0.0051743759610719765 s³⁶⁰²⁸⁷⁹⁷⁰¹⁸⁹⁶³⁹⁷ᐟ⁴⁵⁰³⁵⁹⁹⁶²⁷³⁷⁰⁴⁹⁶

If quantities are involved, the units only show up at around iteration 3, and it is the Pth root that is the issue at that point with the computation of acc failing because of dimension mismatch. Changing the root operation to use // instead of / allows us to proceed further, but I'll need to investigate a bit more.