Trying to get how to use it - reading Documentation unfortunately doesn’t help.
The use case - a model described by a formula, and a vector of (t, y) value pairs representing e.g. experimental data. Like the example problem in the docs of the LsqFit.jl package.
Can you show what you tried so far?
I’ve tried to read the whole package documentation and analyze the example, spend all in all probably an hour, made different attempts to write an example of my own (which errored, and which I probably can’t reproduce now), and failed. The example problem in the documentation has a different structure from the (quite common) problem I’d like to be able to solve. I could find e.g. no explanation why ‘NonlinearFunction’ should be there.
There is an obvious mismatch between me and the docs, and you may consider the problem is not the documentation. On the other side, I had no such problems with the docs of some other SciML (sub-) packages, including Problem Type 1: Solving Nonlinear Systems of Equations from the same package.
using NonlinearSolve
function nlls!(du, u, p)
du[1] = 2u[1] - 2
du[2] = u[1] - 4u[2]
du[3] = 0
end
u0 = [0.0, 0.0]
prob = NonlinearLeastSquaresProblem(
NonlinearFunction(nlls!, resid_prototype = zeros(3)), u0)
solve(prob)
Here’s an example. What did you try? I’ll need a bit more feedback from you to figure out where the disconnect is here, like see what you misinterpreted from some example.
E.g. this:
function nlls!(du, u, p)
du[1] = 2u[1] - 2
du[2] = (u[1] - 4u[2])^2 + 0.1
du[3] = 0
end
ERROR: LinearAlgebra.LAPACKException(2)
@avikpal maybe we should default to pivoted QR?
The function nlls! has a root du.==0. A typical least squares problem doesn’t. Why such an example?
to get an error, I’ve changed the function nlls!
du[2] = (u[1] - 4u[2])^2 + 0.1
Probably, change the linear solve default right?
I use OSQP for things like that. See, for instance, this notebook for a restricted LS: FinancialTheoryMSc/Ch08_PerfEval.ipynb at master · PaulSoderlind/FinancialTheoryMSc · GitHub
Paul thank you for the suggestion. Actually there are quite a few packages out there, as discussed in this thread as well as here.
For most cases, LsqFit.jl just does the job, however I got some issues with it, and @ChrisRackauckas recommended to avoid it.
What I need is just a package which I could use for simple curve fitting. If there is a SciML package doing what I need, that’d be my default choice.
There’s ongoing effort for that: Migrate to SciML and extend? · Issue #35 · pjabardo/CurveFit.jl · GitHub. We’re close.
Your original problem here just had something in the defaults, which we can handle today. Change the default QR to ColumnNorm by avik-pal · Pull Request #492 · SciML/LinearSolve.jl · GitHub May take a bit of discussion though.
OK, now I got it.
using NonlinearSolve, Plots
function curve_fit(model, u0, xdata, ydata, p)
data = (xdata, ydata, p)
function lossfn!(du, u, data)
(xs, ys, p) = data
du .= model.(xs, Ref(u), Ref(p)) .- ys
return nothing
end
prob = NonlinearLeastSquaresProblem(
NonlinearFunction(lossfn!, resid_prototype = similar(ydata)), u0, data)
sol = solve(prob)
u = sol.u
fit = model.(xdata, Ref(u), Ref(p))
return (;sol, fit)
end
define a model:
function expmodel(x, u, t₀=0)
y0 = u[1]
a = u[2]
τ = u[3]
return y0 + a * exp(-(x-t₀)/τ)
end
synthesize data:
y0, a, τ, t₀ = 0.0, 1.0, 2.0, 0
xs = 0:0.1:3;
ys0 = expmodel.(xs, Ref([y0, a, τ]), Ref(t₀));
ys = @. ys0 + 0.02*randn()
Solve, plot:
(;sol, fit) = curve_fit(expmodel, [y0, a, τ ], xs, ys, t₀)
plot(xs, [ys, fit])
