Create Measurements package extension

Project.toml

@@ -10,16 +10,24 @@ NLSolversBase = "d41bc354-129a-5804-8e4c-c37616107c6c"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 
+[weakdeps]
+Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
+
+[extensions]
+MeasurementsExt = "Measurements"
+
 [compat]
 Distributions = "0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
 ForwardDiff = "0.10"
+Measurements = "2.10"
 NLSolversBase = "7.5"
 StatsAPI = "1"
 julia = "1.6"
 
 [extras]
+Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "StableRNGs"]
+test = ["Test", "StableRNGs", "Measurements"]

ext/MeasurementsExt.jl

@@ -0,0 +1,35 @@
+module MeasurementsExt
+
+using LsqFit
+isdefined(Base, :get_extension) ? (using Measurements) : (using ..Measurements)
+
+function LsqFit.curve_fit(
+    model,
+    xdata::AbstractArray,
+    ydata::AbstractArray{Measurement{T}} where T,
+    p0::AbstractArray;
+    inplace=false,
+    kwargs...,
+)
+    y = Measurements.value.(ydata)
+    ye = Measurements.uncertainty.(ydata)
+    wt = ye .^ -2
+    curve_fit(model, xdata, y, wt, p0; inplace=inplace, kwargs...)
+end
+
+function LsqFit.curve_fit(
+    model,
+    jacobian,
+    xdata::AbstractArray,
+    ydata::AbstractArray{Measurement{T}} where T,
+    p0::AbstractArray;
+    inplace=false,
+    kwargs...,
+)
+    y = Measurements.value.(ydata)
+    ye = Measurements.uncertainty.(ydata)
+    wt = ye .^ -2
+    curve_fit(model, jacobian, xdata, y, wt, p0; inplace=inplace, kwargs...)
+end
+
+end # module

src/LsqFit.jl

@@ -28,6 +28,16 @@ using StatsAPI: coef, confint, dof, nobs, rss, stderror, weights, residuals, vco
 
 import Base.summary
 
+# 
+function __init__()
+    @static if !isdefined(Base, :get_extension)
+        @require Measurements="eff96d63-e80a-5855-80a2-b1b0885c5ab7" begin
+             include("../ext/MeasurementsExt.jl")
+         end
+    end
+end
+
+
 include("geodesic.jl")
 include("levenberg_marquardt.jl")
 include("curve_fit.jl")

test/measurements.jl

@@ -0,0 +1,52 @@
+using LsqFit
+using Measurements
+using StableRNGs
+
+@testset "MeasurementsExt" begin
+    model(x, p) = p[1] .* exp.(-x .* p[2])
+    ptrue = [10, 0.3]
+    x = LinRange(0, 2, 50)
+    y0 = model(x, ptrue)
+    σ = rand(1:5, 50)
+    y = y0 .± σ
+    wt = σ .^ -2
+    fit0 = curve_fit(model, x, y0, wt, ptrue) # fit to data using weights
+    fit1 = curve_fit(model, x, y, ptrue) # fit to data using Measurements
+    @test coef(fit0) ≈ coef(fit1)
+
+    # some example data
+    rng = StableRNG(125)
+    x = range(0, stop = 10, length = 20)
+    y0 = model(x, [1.0, 2.0]) + 0.01 * randn(rng, length(x))
+    σ = 0.01
+    y = y0 .± σ
+    wt = ones(length(x)) * σ .^ -2
+    p0 = [0.5, 0.5]
+
+    for ad in (:finite, :forward, :forwarddiff)
+        fit0 = curve_fit(model, x, y0, wt, p0; autodiff = ad)
+        fit1 = curve_fit(model, x, y, p0; autodiff = ad)
+        for fit in (fit0, fit1)
+            @test norm(fit.param - [1.0, 2.0]) < 0.05
+            @test fit.converged
+
+            # can also get error estimates on the fit parameters
+            errors = margin_error(fit, 0.05)
+            @test norm(errors - [0.021, 0.093]) < 0.01
+        end
+    end
+    # if your model is differentiable, it can be faster and/or more accurate
+    # to supply your own jacobian instead of using the finite difference
+    function jacobian_model(x, p)
+        J = Array{Float64}(undef, length(x), length(p))
+        J[:, 1] = exp.(-x .* p[2])             #dmodel/dp[1]
+        J[:, 2] = -x .* p[1] .* J[:, 1]           #dmodel/dp[2]
+        J
+    end
+    jacobian_fit0 = curve_fit(model, jacobian_model, x, y0, wt, p0;show_trace=true)
+    jacobian_fit1 = curve_fit(model, jacobian_model, x, y, p0;show_trace=true)
+    for jacobian_fit in (jacobian_fit0, jacobian_fit1)
+        @test norm(jacobian_fit.param - [1.0, 2.0]) < 0.05
+        @test jacobian_fit.converged
+    end
+end

test/runtests.jl

@@ -4,7 +4,7 @@
 using LsqFit, Test, LinearAlgebra
 using NLSolversBase
 
-my_tests = ["curve_fit.jl", "levenberg_marquardt.jl", "curve_fit_inplace.jl", "geodesic.jl"]
+my_tests = ["curve_fit.jl", "levenberg_marquardt.jl", "curve_fit_inplace.jl", "geodesic.jl", "measurements.jl"]
 
 println("Running tests:")