LsqFit restored

Manifest.toml

@@ -153,6 +153,12 @@ version = "1.0.2"
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
+[[Distributions]]
+deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
+git-tree-sha1 = "0fc424e725eaec6ea3e9fa8df773bee18a1ab503"
+uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
+version = "0.24.14"
+
 [[DocStringExtensions]]
 deps = ["LibGit2", "Markdown", "Pkg", "Test"]
 git-tree-sha1 = "50ddf44c53698f5e784bbebb3f4b21c5807401b1"
@@ -454,6 +460,12 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 [[Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
+[[LsqFit]]
+deps = ["Distributions", "ForwardDiff", "LinearAlgebra", "NLSolversBase", "OptimBase", "Random", "StatsBase"]
+git-tree-sha1 = "0c33987800fbc37edf3a5fd94520e6b06783a63b"
+uuid = "2fda8390-95c7-5789-9bda-21331edee243"
+version = "0.12.0"
+
 [[MKL_jll]]
 deps = ["IntelOpenMP_jll", "Libdl", "Pkg"]
 git-tree-sha1 = "eb540ede3aabb8284cb482aa41d00d6ca850b1f8"
@@ -543,6 +555,12 @@ git-tree-sha1 = "8a8208d3a1b97994d15ebcdae235d332b4a69e78"
 uuid = "429524aa-4258-5aef-a3af-852621145aeb"
 version = "1.2.3"
 
+[[OptimBase]]
+deps = ["NLSolversBase", "Printf", "Reexport"]
+git-tree-sha1 = "4c26a757fbb5b1893b7df19a44e21762d8f8e470"
+uuid = "87e2bd06-a317-5318-96d9-3ecbac512eee"
+version = "2.0.1"
+
 [[Opus_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "f9d57f4126c39565e05a2b0264df99f497fc6f37"
@@ -560,6 +578,12 @@ git-tree-sha1 = "1b556ad51dceefdbf30e86ffa8f528b73c7df2bb"
 uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
 version = "8.42.0+4"
 
+[[PDMats]]
+deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"]
+git-tree-sha1 = "95a4038d1011dfdbde7cecd2ad0ac411e53ab1bc"
+uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
+version = "0.10.1"
+
 [[PGFPlotsX]]
 deps = ["ArgCheck", "DataStructures", "Dates", "DefaultApplication", "DocStringExtensions", "MacroTools", "Parameters", "Requires", "Tables"]
 git-tree-sha1 = "1adde3d07cce96b6a3bb88572612db4bd9d6153b"
@@ -671,6 +695,18 @@ git-tree-sha1 = "cfbac6c1ed70c002ec6361e7fd334f02820d6419"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
 version = "1.1.2"
 
+[[Rmath]]
+deps = ["Random", "Rmath_jll"]
+git-tree-sha1 = "86c5647b565873641538d8f812c04e4c9dbeb370"
+uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
+version = "0.6.1"
+
+[[Rmath_jll]]
+deps = ["Libdl", "Pkg"]
+git-tree-sha1 = "d76185aa1f421306dec73c057aa384bad74188f0"
+uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
+version = "0.2.2+1"
+
 [[Roots]]
 deps = ["Printf"]
 git-tree-sha1 = "8f743e4f4368d1d753f3806bf635899dad6b4847"
@@ -734,12 +770,22 @@ git-tree-sha1 = "7bab7d4eb46b225b35179632852b595a3162cb61"
 uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 version = "0.33.2"
 
+[[StatsFuns]]
+deps = ["Rmath", "SpecialFunctions"]
+git-tree-sha1 = "3b9f665c70712af3264b61c27a7e1d62055dafd1"
+uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+version = "0.9.6"
+
 [[StructArrays]]
 deps = ["DataAPI", "Tables"]
 git-tree-sha1 = "ad1f5fd155426dcc879ec6ede9f74eb3a2d582df"
 uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 version = "0.4.2"
 
+[[SuiteSparse]]
+deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
+uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+
 [[TableTraits]]
 deps = ["IteratorInterfaceExtensions"]
 git-tree-sha1 = "b1ad568ba658d8cbb3b892ed5380a6f3e781a81e"

Project.toml

@@ -8,6 +8,7 @@ ADerrors = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
 BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

src/juobs.jl

 module juobs
-using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, Optim
+using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, Optim
 import Statistics: mean
 
 include("juobs_types.jl")

src/juobs_tools.jl

@@ -277,18 +277,15 @@ function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}
     
     # Generate chi2 + solver
     chisq = gen_chisq(model, xdata, yer)
-    min_fun(t) = chisq(t, yval)
-    sol = optimize(min_fun, fill(0.5, param), method=LBFGS())
-    par = Optim.minimizer(sol)
-
-    # Info
-    (upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, par, ydata) : fit_error(chisq, par, ydata, wpm)
+    fit = curve_fit(model, xdata, yval, 1.0 ./ yer.^2, fill(0.5, param))
+    (upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata) : fit_error(chisq, coef(fit), ydata, wpm)
+    #Info
     for i = 1:length(upar)
         isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm)
         print("\n Fit parameter: ", i, ": ")
         details(upar[i])
     end
-    println("Chisq / chiexp: ", sol.minimum, " / ", chi_exp, " (dof: ", length(yval) - param,")")
+    println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", length(yval) - param,")")
     return upar
 end
 
@@ -324,10 +321,7 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
     data[Nalpha*Ndata+1:end] = ydata
 
     # Guess
-    chisq = gen_chisq(model, xval, yer)
-    min_fun_cons(t) = chisq(t, yval)
-    sol_cons = optimize(min_fun_cons, fill(0.5, param), method=LBFGS())
-    par_cons = Optim.minimizer(sol_cons)
+    fit = curve_fit(model, xval, yval, 1.0 ./ yer.^2, fill(0.5, param))
 
     # Generate chi2 + solver
     if covar
@@ -343,13 +337,13 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
         C = isnothing(wpm) ? [ADerrors.cov(aux[k]) for k = 1:Ndata] : [ADerrors.cov(aux[k], wpm) for k = 1:Ndata]
         chisq_full_cov(p, d) = get_chi2_cov(model, d, C, p, Nalpha)
         min_fun_cov(t) = chisq_full_cov(t, dat)
-        sol = optimize(min_fun_cov, vcat(par_cons, dat[1:Nalpha*Ndata]), method=LBFGS())
+        sol = optimize(min_fun_cov, vcat(fit.param, dat[1:Nalpha*Ndata]), method=LBFGS())
         
         (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full_cov, Optim.minimizer(sol), data) : fit_error(chisq_full_cov, Optim.minimizer(sol), data, wpm)
     else
         chisq_full(p, d) = get_chi2(model, d, ddat, p, Nalpha)
         min_fun(t) = chisq_full(t, dat)
-        sol = optimize(min_fun, vcat(par_cons, dat[1:Nalpha*Ndata]), method=LBFGS())
+        sol = optimize(min_fun, vcat(fit.param, dat[1:Nalpha*Ndata]), method=LBFGS())
         (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full, Optim.minimizer(sol), data) : fit_error(chisq_full, Optim.minimizer(sol), data, wpm)
     end