fit_routine with Optim

Project.toml

@@ -7,7 +7,6 @@ version = "0.1.0"
 ADerrors = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
 BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
-LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"

src/juobs.jl

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

src/juobs_obs.jl

@@ -519,7 +519,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
     
     model(x, p) = get_model(x, p, npol)
 
-    par = fit_routine(model, x, t2E, npol)
+    par, ch2 = fit_routine(model, x, t2E, npol)
     fmin(x, p) = model(x, p) .- 0.3
     t0 = root_error(fmin, t[nt0], par)
     if pl
@@ -631,7 +631,7 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
     
     model(x, p) = get_model(x, p, npol)
 
-    par = fit_routine(model, x, t2E, npol)
+    par, ch2 = fit_routine(model, x, t2E, npol)
     fmin(x, p) = model(x, p) .- 0.3
     t0 = root_error(fmin, t[nt0], par)
     if pl

src/juobs_tools.jl

@@ -560,13 +560,56 @@ function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}
     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
+    #for i = 1:length(upar)
+    #    isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm)
+    #    print("\n Fit parameter: ", i, ": ")
+    #    details(upar[i])
+    #end
+    chis2_corrected = (length(yval) - param) * chisq(coef(fit), ydata) / chi_exp
     println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", length(yval) - param,")")
-    return upar
+    println("Chisq corrected: ", chis2_corrected)
+    return upar, value(chis2_corrected)
+end
+
+function fit_routine(model::Vector{Function}, xdata::Vector{Array{Float64, N}} where N, ydata::Vector{Array{uwreal, N}} where N, param::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    if !(length(model) == length(xdata) == length(ydata))
+        error("Dimension mismatch")
+    end
+    N = length(model)
+
+    dat = ydata[1]
+    idx = Vector{Vector{Int64}}(undef, N)
+    e = Vector{Vector{Float64}}(undef, N)
+    j = 1
+    for i = 1:N
+        if isnothing(wpm)
+            uwerr.(ydata[i])
+        else
+            [uwerr(yaux, wpm) for yaux in ydata[i]]
+        end
+        
+        e[i] = err.(ydata[i])
+        if i > 1 
+            dat = vcat(dat, ydata[i])
+        end
+        stp = j + length(ydata[i]) - 1
+        idx[i] = collect(j:stp)
+        j = stp + 1
+    end
+
+    chisq = (par, dat) -> sum([sum((dat[idx[i]] .- model[i](xdata[i], par)).^2 ./e[i].^2) for i=1:N])
+    min_fun(t) = chisq(t, value.(dat))
+    
+    p = fill(0.5, param)
+    sol = optimize(min_fun, p, LBFGS()) 
+    
+    (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq, sol.minimizer, dat) : fit_error(chisq, sol.minimizer, dat, wpm)
+    
+    println("Chisq / chiexp: ", min_fun(sol.minimizer), " / ", chi2_exp, " (dof: ", length(dat) - param,")")
+    chis2_corrected = (length(dat) - param) * min_fun(sol.minimizer) / chi2_exp
+    println("Chisq corrected: ", chis2_corrected)
+    
+    return upar, chis2_corrected
 end
 
 function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; 
@@ -617,30 +660,35 @@ 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(fit.param, dat[1:Nalpha*Ndata]), LevenbergMarquardt())
+        sol = optimize(min_fun_cov, vcat(fit.param, dat[1:Nalpha*Ndata]), method=LBFGS())
         
         (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full_cov, sol.minimizer, data) : fit_error(chisq_full_cov, sol.minimizer, data, wpm)
         println("Chisq / chiexp: ", min_fun_cov(sol.minimizer), " / ", chi2_exp, " (dof: ", length(ydata) - param,")")
+        chis2_corrected = (length(ydata) - param) * min_fun_cov(sol.minimizer) / chi2_exp
+        println("Chisq corrected: ", chis2_corrected)
     else
         chisq_full(p, d) = get_chi2(model, d, ddat, p, Nalpha)
         min_fun(t) = chisq_full(t, dat)
-        sol = optimize(min_fun, vcat(fit.param, dat[1:Nalpha*Ndata]), LevenbergMarquardt())
+        sol = optimize(min_fun, vcat(fit.param, dat[1:Nalpha*Ndata]), method=LBFGS())
         
         (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full, sol.minimizer, data) : fit_error(chisq_full, sol.minimizer, data, wpm)
         println("Chisq / chiexp: ", min_fun(sol.minimizer), " / ", chi2_exp, " (dof: ", length(ydata) - param,")")
-
+        chis2_corrected = (length(ydata) - param) * min_fun(sol.minimizer) / chi2_exp
+        println("Chisq corrected: ", chis2_corrected)
     end
 
     #### chisq_full, min_fun out of conditional ->
     #### COMPILER WARNING ** incremental compilation may be fatally broken for this module **
 
     # Info
+    #=
     for i = 1:length(upar)
         isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm)
         print("\n Fit parameter: ", i, ": ")
         details(upar[i])
     end
-    return upar
+    =#
+    return upar, chis2_corrected
     
 end