fit_routine added (juobs_tool) and small modification in juobs_linalg

src/juobs_linalg.jl

@@ -88,13 +88,13 @@ It returns a vector array eff_en where each entry  eff_en[t] contains the first
 function energies(evals::Vector{Vector{uwreal}})
     time = length(evals) 
     n = length(evals[1])
-    eff_en = Vector{Vector}(undef, n)
+    eff_en = Array{Array{uwreal}}(undef, n)
     aux_en = Vector{uwreal}(undef, time-1)
     for i in 1:n
         #aux_en = Vector{uwreal}(undef, time-1)
         for t in 1:time-1
             ratio = evals[t][i] / evals[t+1][i]
-            aux_en[t] = log(sqrt(ratio * ratio))
+            aux_en[t] = 0.5*log(ratio * ratio)
         end
         uwerr.(aux_en)
         eff_en[i] = copy(aux_en)

src/juobs_tools.jl

@@ -140,6 +140,35 @@ Computes the results of a linear interpolation/extrapolation in the y axis
 """
 y_lin_fit(par::Vector{uwreal}, x::Union{uwreal, Float64}) = par[1] + par[2] * x
 
+@doc raw"""
+fit_routine(model::Function, ydata::Array{uwreal}, param::Int64=3)
+
+Given a model function with a number param of parameters and an array of uwreal,
+this function fit ydata with the given model and print  fit information
+The method return an array upar with the best fit parameters with their errors.
+'''@example
+@. model(x,p) = p[1] + p[2] * exp(-(p[3]-p[1])*x)
+fit_routine(model, ydata, param=3)
+"""
+function fit_routine(model::Function, ydata::Array{uwreal}, param::Int64=3)
+    yval = value.(ydata)
+    yer = err.(ydata)
+    xdata = collect(1:length(ydata))
+    chisq = gen_chisq(model, xdata, yer)
+    fit = curve_fit(model, xdata, yval, 1.0 ./ yer.^2, fill(0.5, param))
+    upar, chi_exp = fit_error(chisq, coef(fit), ydata)
+    for i = 1:length(upar)
+        uwerr(upar[i])
+        print("\n Fit parameter: ", i, ": ")
+        details(upar[i])
+    end
+    println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", dof(fit),")")
+    return upar
+end
+function gen_chisq(f::Function, x::Vector{Int64}, err::Vector{Float64})
+    chisq(par, dat) = sum((dat .- f(x,par)).^2 ./err.^2)
+    return chisq
+end
 #TODO: add combined fits
 #=
 using LsqFit