Create Measurements package extension

[chriswaudby]

This PR creates a simple extension for the Measurements package, allowing data with uncertainties to be passed directly for fitting, e.g.:

using LsqFit, Measurements
x = [1.0, 2.0, 3.0]
y = [1 ± 0.1, 2.2 ± 0.2, 2.9 ± 0.3]
model(x, p) = p[1] * x
p0 = [1.]
curve_fit(model, x, y, p0)

The extension creates a simple wrapper for LsqFit.curve_fit that dispatches on ydata::AbstractArray{Measurement{T}} where T. Values and uncertainties are extracted from the input, then weights are calculated as the inverse square of the uncertainties and passed to the main curve_fit function.

A unit test checks that the extension gives the same results as when weights are passed manually.

An __init__ function is defined in LsqFit.jl to load the extension for Julia versions <1.9 (see discussion at https://discourse.julialang.org/t/package-extensions-for-julia-1-9/93397).

[chriswaudby]

Hi @pkofod, haven't made a pull request before so don't know the etiquette, but it's been a few weeks since I submitted this - are you able to take a look? Let me know if there's anything else I should be doing!

[pkofod]

Thanks, not it's not bad etiquette to bump :) Sorry for not seeing this.

I think sometimes there can be corner cases with automatic differentiation where things can go wrong, so if we are to merge this PR, I think the tests have to be significantly extended to cover the various ways which the user-facing functions can be called: different combinations of inputs, different choices for AD, etc.

[TheFibonacciEffect]

This works very well, thanks a lot :D

using LsqFit, Plots, Measurements
@. gauss(x,p) = p[1] + p[2] * exp(-(x-p[3])^2/p[4]^2)
x = range(-1,1,100)
yerr = gauss(x,[0,1,0,0.5]) + 0.1*randn(length(x)) .± 0.1
fit = curve_fit(gauss, x,yerr,Float64[0,1,0,1])
begin
	scatter(x,yerr)
	plot!(x,gauss(x,fit.param);ribbon = gauss(x,fit.param.+stderror(fit)) - gauss(x,fit.param.-stderror(fit)))
end
savefig("fit.png")

[donovaly]

This would be a helpful extension. Why is this PR stuck?

The PR has a test, the request to write more tests is a bit abstract because it is not clear what exactly should be tested.

[chriswaudby]

I've just added an additional method for compatibility with Jacobians, and some tests for autodiff - matching those for the main curve_fit method. Hopefully this addresses the original concern?

[pkofod]

The PR has a test, the request to write more tests is a bit abstract because it is not clear what exactly should be tested.

I respectfully disagree. Since the package supports various ways of calculating derivatives by automatic differentiation it is not just a given that everything will work because one version of curve_fit works. There is such a large surface to cover in terms of potentially problematically typed inner functions and calls to dependency code and so on that need to work with measurements and autodiff. Just to take one example.

I've just added an additional method for compatibility with Jacobians, and some tests for autodiff - matching those for the main curve_fit method. Hopefully this addresses the original concern?

I had a quick look and yes that appears to cover exactly what I was looking for.

[donovaly]

There is such a large surface to cover in terms of potentially problematically typed inner functions and calls to dependency code and so on that need to work with measurements and autodiff.

The point is that you write requiring a lot of background knowledge about the history of LsqFit etc. Therefore it was not clear to me what exactly you mean. For example despite being a Ph.D. physicist I never hear of the term automatic differentiation yet. And I don't understand what is meant by "automatic".

However, I am happy that this PR move an because it is a cool feature.

[pkofod]

The point is that you write requiring a lot of background knowledge about the history of LsqFit etc. Therefore it was not clear to me what exactly you mean. For example despite being a Ph.D. physicist I never hear of the term automatic differentiation yet. And I don't understand what is meant by "automatic".

@chriswaudby seemed to understand it just fine. It has nothing to do with having a Ph.D. in physics or computer science or anything else. It has to do about project management and making sure PRs do not add code that is insufficiently tested. It also has nothing to do with the history of LsqFit, but only the features you can find here now. Automatic differentiation, well you can just google it. It's a quite old method that allows you to calculate exact derivatives instead of finite differences.

I appreciate the interest, but there is more to package development than adding some code that appears to work. I hope you appreciate that if this code is merged without sufficient tests and contains bugs that could have been avoided, I'm the one who receives the complaints. Code is proposed and reviewed in an iterative process and is merged when the reviewer approves. Again, the PR author seemed to understand it just fine which is what is important here.

[chriswaudby]

Ok, to get the discussion back on track - I appreciate the nudge on this PR from @donovaly to write more tests, which picked up that explicit jacobians weren't originally being handled correctly. That's now fixed and tests are covering AD, so are there any other issues outstanding?

[pkofod]

Is the sigma here not a standard deviation? Then I think you're using the weight from c.f. the discussion in the other active thread about weights being wrong :) For the vector format you need inverse standard deviations, for the matrix input it's a variance-covariance matrix.

[donovaly]

@chriswaudby seemed to understand it just fine. It has nothing to do with having a Ph.D. in physics or computer science or anything else. It has to do about project management and making sure PRs do not add code that is insufficiently tested.

The PR was stuck and he did not reply for a long time to you. I therefore gave (unrequested) feedback, that you should more clearly say what you want to be improved in a PR. Exactly what test for what purpose or so. I came here as a complete outsider, wondered why the PR was stuck and did not understood what exactly you wanted. My personal experience as PR reviewer and maintainer is that when things are not clear enough, PRs got stuck or you have to take over and implement the missing things. We do this in our spare time and when I feel things get too complicated for me or I don't understand, I will not continue. I speak here for myself and this is also the experiences I made with the PRs I reviewed.

However, I am happy that the PR is alive and I look forward to get it merged soon.

Regarding the terminology, I wanted to give feedback because also in your other replies to me you refer to things that newcomers might not know (package history, maintainer history, technical terms). You can take this as criticism, but it was meant just as feedback. For example I was involved for many years as maintainer in LyX and made often the mistake to expect too much LaTeX background knowledge instead of only that portion necessary for the particular PR.
Here, of course I can google what is meant by "automatic differentiation", but it is not good to be forced to do. Is it necessary to understand the whole concept of automatic differentiation or only a portion, if only a portion, what portion, what corner case, etc. to get the PR merged or not, that was my point. For example for you automatic differentiation is so normal that you abbreviate it directly as "AD". And for me I felt bad that I never every heard from it yet.

And we don't know each other, so I think it is no offense to give just feedback. It is sad that you feel offended by this.

[pkofod]

Ci is run on v1.6 here and I think maybe the reason why the tests fail because extensions were not included in Julia then. Maybe we need restrict Julia to be >= v1.9 since we're adding an extension here and adjust the CI actions accordingly.

[aplavin]

This PR's implementation seems to ignore correlation between measurements.
The point of Measurements.jl, the only reason they use a very involved implementation (compared to a much simpler one as in Uncertain.jl), is tracking correlation between variables. Ignoring them for fitting can easily be surprising for the user if they even catch it – otherwise it directly leads to incorrect results.

[TheFibonacciEffect]

I think maybe the covariance matrix might have to be used here instead, but I dont know if how the weights would need to be adjusted for that.

[TheFibonacciEffect]

Maybe a monte carlo simulation of a fit could be used as a test case? That way, we know that the uncertainties are taken into account correctly.
Possibly something like this:

using LsqFit
using LinearAlgebra
using Distributions
# create datapoints with noise and correalation
N = 10
x = range(0, stop=10, length=N)
corr = 0.5^2 * diagm(0 => 2*ones(N), 1 => ones(N-1), -1 => ones(N-1))

# define model function
@. model(x, p) = p[1] * sin(p[2] * x)

p0 = [1.0, 1.0]
p = []
for i in 1:10
    y = model(x, p0) + rand(MvNormal(corr))
    fit = curve_fit(model, x, y, p0)
    push!(p, fit.param)
end
p_mean = mean(p, dims=1)
p_cov = cov(p)

which samples from a distribution according to a given correalation matrix and makes a fit. The same result should then be reproduced by this code with Measurements.jl with a given covariance matrix, at least approximately.