Support for Unitful.jl's units?

[singularitti]

Would you consider supporting fitting x- and y-data which has units that are defined by Unitful.jl?
For example,

using LsqFit
@. model(x, p) = p[1]*exp(-x*p[2])
xdata = range(0, stop=10, length=20) .* u"m"
ydata = model(xdata, [1.0u"kg" 2.0u"m^-1"]) + 0.01u"g" * randn(length(xdata))
p0 = [0.5u"kg", 0.5u"m^-1"]
fit = curve_fit(model, xdata, ydata, p0)

I'm having the following stacktrace:

ERROR: MethodError: no method matching (Quantity{Float64})(::Float64)

Closest candidates are:
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:790
  (::Type{<:Unitful.AbstractQuantity})(::Union{Dates.Day, Dates.Hour, Dates.Microsecond, Dates.Millisecond, Dates.Minute, Dates.Nanosecond, Dates.Second, Dates.Week})
   @ Unitful ~/.julia/packages/Unitful/J4AJj/src/dates.jl:60
  (::Type{<:Unitful.AbstractQuantity})(::Dates.CompoundPeriod)
   @ Unitful ~/.julia/packages/Unitful/J4AJj/src/dates.jl:135
  ...

Stacktrace:
 [1] alloc_DF(x::Vector{Quantity{Float64}}, F::Vector{Quantity{Float64, 𝐌, Unitful.FreeUnits{(kg,), 𝐌, nothing}}})
   @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/objective_types/abstract.jl:19
 [2] lmfit(f::LsqFit.var"#18#20"{…}, p0::Vector{…}, wt::Vector{…}; autodiff::Symbol, kwargs::@Kwargs{})
   @ LsqFit ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:71
 [3] lmfit(f::Function, p0::Vector{Quantity{Float64}}, wt::Vector{Quantity{Float64, 𝐌, Unitful.FreeUnits{(kg,), 𝐌, nothing}}})
   @ LsqFit ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:54
 [4] curve_fit(model::typeof(model), xdata::StepRangeLen{…}, ydata::Vector{…}, p0::Vector{…}; inplace::Bool, kwargs::@Kwargs{})
   @ LsqFit ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:140
 [5] curve_fit(model::Function, xdata::StepRangeLen{…}, ydata::Vector{…}, p0::Vector{…})
   @ LsqFit ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:123

The problem seems to be in this line:

alloc_DF(x, F) = eltype(x)(NaN) .* vec(F) .* vec(x)'

Where x has eltype of Quantity{Float64}, which cannot apply onto NaN. I guess this can be changed since we are only allocating an array, where NaN does not really matter here.

Versions:

• LsqFit v0.15.0
• Julia: v1.10.0-rc1

Similar: #248 #143

[chriswaudby]

My understanding is that Measurements is designed to work nicely with whatever package is used for units: https://juliaphysics.github.io/Measurements.jl/dev/usage/#Error-Propagation-of-Numbers-with-Units. So it ought to be possible to create a package extension to handle Unitful quantities - whether with uncertainties or not.

The easiest way would be a lightweight function, like I did for Measurements in PR #248, that removes the Quantity wrapper and calls LsqFit again with the underlying data. If that data has uncertainties, they should then get handled correctly by the Measurements package extension. However, something I'm less sure about would be how to handle units within the fitting parameters themselves - would it be safe to simply strip the quantities from these, or is there some conversion magic that would be lost by this?

[singularitti]

Hi @chriswaudby, I was doing it in my package: stripping units before fitting, and re-adding units to the fitted parameters.

would it be safe to simply strip the quantities from these, or is there some conversion magic that would be lost by this

My rule of thumb is: use whatever the input parameters' units are (they may be different from each other). However, consider that the units of xdata and ydata might differ from the parameters' units; therefore, we should consider a fitting technique that is agnostic to the units.

In the above example, model with p0 will give results of dimension "mass," supposedly in u"kg". However, the ydata could be in different units, such as u"g" or u"lbs". Therefore, we should ensure that the ydata and the values produced by the model are in the same unit system.

[chriswaudby]

This might get complicated! Perhaps what we need to do is:

1. extract the units of ydata and p0
2. extract unitless values for ydata and p0
3. create a wrapper function for model:
   1. restores units to the input parameters
   2. call the original Unitful model
   3. convert the output to match the units of ydata
   4. return an output stripped of these units
4. use this wrapped model to call curve_fit
5. return the fitted parameters with the original units restored

Does this seem sensible?

What to do about weights - is it sufficient just to strip these and use numerical values (provided all weights have the same unit)?