Support for `Unitful.jl`'s units?

singularitti commented 10 months ago

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 commented 9 months ago

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 commented 9 months ago

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 commented 9 months ago

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

extract the units of ydata and p0
extract unitless values for ydata and p0
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
use this wrapped model to call curve_fit
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)?

JuliaNLSolvers / LsqFit.jl

Support for `Unitful.jl`'s units? #250