moved fit routine to their file, added it to the juobs module, fix in...

moved fit routine to their file, added it to the juobs module, fix in fit_routine with combined fits

src/juobs.jl

@@ -12,6 +12,7 @@ include("juobs_tools.jl")
 include("juobs_plots.jl")
 include("juobs_obs.jl")
 include("juobs_pvalue.jl")
+include("juobs_fit.jl")
 include("../constants/juobs_const.jl")
 include("../constants/path_csv.jl")
 #include("../constants/juobs_uwreal_const.jl")

src/juobs_fit.jl

+
+"""  
+    plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, [W::VecOrMat{Float64},] wpm::Union{Dict{Int64},Vector{Float{64}},Dict{String,Vector{Float{64}},Nothing}=nothing} )
+
+It computes plateau average of `obs` in the range `plat`. Effectively this perform a fit to a constant. 
+The field `W` is optional. When absent, `W` is a `Vector{Float64}` constructed from the error of `obs`. 
+When `W` is a `Vector{Float64}` it performs an uncorrelated fit, when `W` is a `Matrix{Float64}`, it perfom a correlated fit
+
+# Potential side effects
+
+When `W` is absent, the function acts with `ADerror.uwerr` on `obs` in order to compute the errors.
+This is the reason why a `wpm` may be given as input.
+
+"""
+function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    isnothing(wpm) ? uwerr.(obs) : [uwerr(obs_aux, wpm) for obs_aux in obs]
+    w = 1 ./ err.(obs)[plat[1]:plat[2]].^2
+    av = sum(w .* obs[plat[1]:plat[2]]) / sum(w)
+    return av 
+end
+
+plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, W::Vector{Float64},wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) = sum(W .* obs[plat[1]:plat[2]]) / sum(W)
+
+function lin_fit(x::Vector{<:Real}, v::Vector{Float64}, e::Vector{Float64})
+    sig2 = e .* e
+    S = sum(1 ./ sig2)
+    Sx = sum(x ./ sig2)
+    Sy = sum(v ./ sig2)
+    Sxy = sum(v .* x ./ sig2)
+    Sxx = sum(x .* x ./sig2)
+    delta = S * Sxx - Sx*Sx
+    par = [Sxx*Sy-Sx*Sxy, S*Sxy-Sx*Sy] ./delta
+    #C = [[Sxx/delta, -Sx/delta], [-Sx/delta,  S/delta]]
+    return par
+end
+@doc raw"""
+    lin_fit(x::Vector{<:Real}, y::Vector{uwreal})
+
+Computes a linear fit of uwreal data points y. This method return uwreal fit parameters and chisqexpected.
+
+```@example
+fitp, csqexp = lin_fit(phi2, m2)
+m2_phys = fitp[1] + fitp[2] * phi2_phys
+```
+"""
+function lin_fit(x::Vector{<:Real}, y::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    isnothing(wpm) ? uwerr.(y) : [uwerr(yaux, wpm) for yaux in y]
+    par = lin_fit(x, value.(y), err.(y))
+    chisq(p, d) = sum((d .- p[1] .- p[2].*x).^2 ./ err.(y) .^2)
+    (fitp, csqexp) = fit_error(chisq, par, y)
+    for i = 1:length(fitp)
+        isnothing(wpm) ? uwerr(fitp[i]) : uwerr(fitp[i], wpm)
+        print("\n Fit parameter: ", i, ": ")
+        details(fitp[i])
+    end
+    println("Chisq / chiexp: ", chisq(par, y), " / ", csqexp, " (dof: ", length(x)-length(par),")")
+    return (fitp, csqexp)
+end
+
+@doc raw"""
+    x_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64})
+
+Computes the results of a linear interpolation/extrapolation in the x axis
+"""
+x_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64}) = (y - par[1]) / par[2]
+@doc raw"""
+    y_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64})
+
+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 """
+    get_chi(model::Function,x::Vector,par,data,W::VecOrMat{Float64})
+
+It computes the chi square of a fit. `model` is assumed to be of the type f(xdata,par). 
+If `W::Vector` the chisquare is computed for an uncorrelated fit
+If `W::Matrix` the chisquare is computed for a correlated fit
+"""
+function get_chi(model::Function,x,par,data,W::Array{Float64,2})
+  return sum((data.-model(x,par))' * W * (data.-model(x,par)))
+end
+
+function get_chi(model::Function,x,par,data,W::Vector{Float64})
+  return sum((data.-model(x,par)).^2 .* W)
+end
+
+
+@doc raw"""
+        fit_routine(model::Function,xdata::AbstractArray{<:Real},ydata::AbstractArray{ADerrors.uwreal},npar; kwargs...)
+        fit_routine(model::Function,xdata::AbstractArray{uwreal},ydata::AbstractArray{ADerrors.uwreal},npar; kwargs...)
+
+It executes a fit of the `ydata` with the fit function `model`. 
+
+# Arguments:
+
+  - `model`: Fit function used to fit the data. the function assumes that its called as `model(xdata,parameters)`.
+  - `xdata`: indipendent variable in the fit. It can be an `AbstractArray{<:Real}`, that is the "default" method, or 
+             a `AbstractArray{uwreal}`. In the latter case, the method builds the function
+``math
+F(q_i) = \begin{cases}
+         model(q_{npar+i},q[1:npar])  \quad \text{if} i<= length(xdata)\\
+         q_i \quad \text{elsewhere}
+       \end{cases}
+``
+             and then call the default method with dummy `xdata` and `vcat(ydata,xdata)` as new `ydata`
+  - `ydata`: Data variable that carries the error 
+  - `npar`: number of parameters used in the fit.
+
+# Keyword parameters:
+
+  - `wpm`: window parameter used by ADerrors to computes the errors. Default to `Dict{Int64,Vector{Float64}}()`.
+           This is used both when computing the error of the data and when computing the covariance matrix through `ADerrors.cov`
+           It accept a dictionary with ensemble id as keys (either `String` or `Int64`) and `Vector{Float64}` as value.
+           See also ADerrors documentation.
+
+  - `W::AbstractVecOrMat{Float64}`: Weight matrix/vector. If `W` is a `Vector` the fit is perfomed uncorrelated ,otherwise is correlated.
+                                    default to `Vector{Float64}()`
+
+  - `corr::Bool`: flag used to generated `W` when not given. Default to `true`. If `false`, `W = [1/y.err^2 for y in ydata]`.
+                  If `true` the covariance matrix of the data is used to compute `W` as `W = C^{-1}`.
+
+  - `C::AbstractMatrix`: covariance matrix used to compute W. Default to `Matrix{Float64}(undef, 0, 0)`.
+                         If `W` is not given, and `corr` is `true`, then it is used to compute `W`.
+                         If need to computed `W` but not given,  then `C = ADerrors.cov(ydata,wpm)`. See also ADerrors documentation.
+                         If available, it will be used to better estimate the fit pvalue and the expected chi-square,
+                         but it will not be computed if not available at this point. 
+
+  - `guess::Vector{Float64}`: Initial guess for the parameter used in the fit routine. If not given, default to `fill(0.5,npar)`.
+                              If `xdata isa Vector{uwreal}`, the mean value of `xdata` are used as guess for the relevant `q_i` parameters.
+  
+  - `logfile`: handle to a logfile. If `nothing`, no log will be written. It can be `IO` file or a custom log structure that has
+               an overloading for `print()` and `println()`.
+
+# Returns
+It returns a `NamedTuple` with names:
+  - `:par`:    fit parameter as `Vector{uwreal}`
+  - `:chi2`:   chisquare,
+  - `:chiexp`: chisquare expected 
+  - `:pval`:   pvalue 
+"""
+function fit_routine(model::Function,
+                     xdata::AbstractArray{<:Real},
+                     ydata::AbstractArray{uwreal},
+                     npar::Int64; 
+                     wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}}=Dict{Int64,Vector{Float64}}(),
+                     corr::Bool = true, 
+                     W::AbstractArray{Float64}=Vector{Float64}(),
+                     guess::AbstractVector{Float64} = fill(0.5,npar),
+                     logfile = nothing,
+                     C::AbstractArray{Float64} = Matrix{Float64}(undef,0,0));
+
+    nobs = length(ydata)
+    if length(xdata)!=nobs
+        error("xdata and ydata must be the of same size")
+    end
+
+    [uwerr(y, wpm) for y in ydata]
+    if length(W) ==0
+        if corr 
+            C = length(C) ==0 ? GammaMethod.cov_AD(ydata,wpm) : C
+            W = LinearAlgebra.pinv(C);
+            W = (W'+W)*0.5
+        else
+            W = [1/y.err^2 for y in ydata]
+        end
+    end
+
+    fit = LsqFit.curve_fit(model,xdata,ADerrors.value.(ydata),W,guess)
+
+    chisq(p,d) = get_chi(model,xdata,p,d,W)
+    chi2 = sum(fit.resid.^2)
+    par = fit_error(chisq,LsqFit.coef(fit),ydata,wpm,W,false)
+    # if the fit is correlated  we have already C, no need to recompute it.
+    chiexp,pval =  if length(C) !=0 
+        GammaMethod.chiexp(chisq,LsqFit.coef(fit),value.(ydata),C,W),juobs.pvalue(chisq,chi2,LsqFit.coef(fit),ydata,W,C=C)
+    else  
+        ADerrors.chiexp(chisq,LsqFit.coef(fit),ydata,wpm,W=W),juobs.pvalue(chisq,chi2,LsqFit.coef(fit),ydata,W)
+    end
+
+    if !isnothing(logfile)
+        uwerr.(par);
+        println(logfile, "juobs.fit_routine output:")
+        println(logfile, "\t\tFit routine in [$(xdata[1]),...,$(xdata[end])]")
+        println(logfile, "\t\tparameters :")
+        for p in par
+            println(logfile,"\t\t\t",p)
+        end
+        println(logfile, "\t\tchi2:       ", chi2)
+        println(logfile, "\t\tchiexp:     ", chiexp)
+        println(logfile, "\t\tpvalue:     ", pval)
+    end
+
+    return (par=par,chi2=chi2,chiexp=chiexp,pval=pval)
+end
+
+"""
+        check_size(M,n)
+
+Check if `M` has sizes all equal to n or if is empty
+"""
+function check_sizes(M,n)
+    if length(M) ==0 || all(size(M) .== n)
+        return true
+    end
+    return false
+end
+
+function fit_routine(model::Function,
+                     xdata::AbstractArray{ADerrors.uwreal},
+                     ydata::AbstractArray{ADerrors.uwreal},
+                     npar::Int64;
+                     W::AbstractArray = Float64[],
+                     wpm::AbstractDict = Dict{Int64,Vector{Float64}}(),
+                     corr::Bool = true,
+                     guess::AbstractArray = fill(0.5,npar),
+                     logfile = nothing,
+                     C::AbstractArray = Float64[])
+    Nx = length(xdata)
+    if (Nx != length(ydata))
+        error("[juobs] Error: xdata and ydata must have the same dimension")
+    end
+    Ndata = 2Nx; ## in total we have 2*lengt(xdata) data points
+    data = vcat(ydata,xdata)
+    if !check_sizes(W,Ndata)
+        error("[juobs] Error: W does not have the correct sizes")
+    end
+    if !check_sizes(C,Ndata)
+        error("[juobs] Error: C does not have the correct sizes")
+    end
+
+    if corr
+        C = length(C) ==0 ? GammaMethod.cov_AD(data) : C
+        W = length(W) ==0 ? LinearAlgebra.pinv(W) : W
+    else
+        W = length(W) ==0 ? [1/yy.err^2 for yy in data] : W
+    end
+
+    f(q,p) = [i<=Nx ? model(q[i],p) : q[i-Nx] for i in 1:Ndata]
+    fit_func(x,p) = f(view(p,npar+1:npar+Nx),view(p,1:npar))
+    new_guess = [guess; value.(xdata)]
+    fit = fit_routine(fit_func,zeros(Ndata),data,npar+Nx,
+                      guess=new_guess, wpm=wpm,corr=corr,
+                      W=W,C=C,logfile=nothing)
+
+    if !isnothing(logfile) 
+        uwerr.(fit.par)
+        println(logfile, "juobs.fit_routine output:")
+        println(logfile, "\t\txdata carries uncertainty")
+        println(logfile, "\t\tFit routine in [$(xdata[1].mean),...,$(xdata[end].mean)]")
+        println(logfile, "\t\tFit parameters:")
+        for i in 1:npar
+            println(logfile,"\t\t\t",fit.par[i])
+        end
+        println(logfile, "\t\tauxiliary fit parameter:")
+        println(logfile, "\t\t\txdata[i] \t q[i]")
+        for i in 1:Nx
+            println(logfile, "\t\t\t",xdata[i]," \t ", fit.par[npar+i])
+        end
+        println(logfile, "\t\tchi2:       ", fit.chi2)
+        println(logfile, "\t\tchiexp:     ", fit.chiexp)
+        println(logfile, "\t\tpvalue:     ", fit.pval)
+    end
+    return fit
+end
+
+
+@doc raw"""
+        fit_routine(model::AbstractArray{Function},
+                    xdata::AbstractArray{<:AbstractArray},
+                    ydata::AbstractArray{<:AbstractArray{ADerrors.uwreal}},npar; kwargs...)
+       
+
+It executes a combined fit of the `ydata` with the fit functions in `model`. The method define the general model `F(X_{mi},p) = model[m](xdata[i],p)`
+and call `fit_routine(F,vcat(xdata),vcat(ydata),npar;kwargs....)`. The `kwargs` parameters are accordingly updated to reflect the new function.
+
+# Returns
+It returns a `NamedTuple` with names:
+  - `:par`:    fit parameter as `Vector{uwreal}`
+  - `:chi2`:   chisquare,
+  - `:chiexp`: chisquare expected 
+  - `:pval`:   pvalue 
+"""
+function fit_routine(models::AbstractArray{Function},
+                     xdata::AbstractArray{<:AbstractArray},
+                     ydata::AbstractArray{<:AbstractArray{ADerrors.uwreal}},
+                     npar::Int64;
+                     W::AbstractArray = Float64[],
+                     C::AbstractArray = Float64[],
+                     wpm::AbstractDict = Dict{Int64,Vector{Float64}}(),
+                     corr::Bool = true,
+                     guess::AbstractArray = fill(0.5,npar),
+                     logfile = nothing,
+                     )
+    Nmodel = length(models)
+    if Nmodel != length(xdata) != length(ydata)
+        error("[juobs] Error: you need the same number of model, xdata and ydata")
+    end
+    if !(length(W) == 0 || length(W)==Nmodel)
+        error("[juobs] Error: You need to pass a W matrix for model or none")
+    end
+    if !(length(C) == 0 || length(C)==Nmodel)
+        error("[juobs] Error: You need to pass a C matrix for model or none")
+    end
+    ndata = length.(xdata)
+    if !all(length.(ydata).== ndata)
+        error("[juobs] Error: Mismatch in xdata and ydata. Make sure that the data corresponds")
+    end
+    Ntotal = sum(ndata);
+    X = vcat(xdata...)
+    Y = vcat(ydata...)
+
+    function make_matrix(M)
+        aux = zeros(Ntotal,Ntotal)
+        ie = 0
+        for m in 1:Nmodel
+            is = ie+1
+            ie += ndata[m];
+            aux[is:ie,is:ie] .= M[m]
+        end
+        return aux
+    end
+
+    WW = length(W)==0 ? W : make_matrix(W);
+    CC = length(C)==0 ? C : make_matrix(C);
+
+    function Model(x,p)
+        res = zeros(Nmodel)
+        res[1] = models[1](x[1:ndata[1]],p)
+        for i in 2:Nmodel
+            rd = (ndata[i-1]+1):ndata[i]
+            rp = (npar[i-1]+1):npar[i]
+            res[i] = models[i](x[rd],p)
+        end
+        return res
+    end
+    fit = fit_routine(Model,X,Y,npar,guess = guess,W=WW,C=CC,corr=corr,logfile=nothing,wpm=wpm)
+
+    if !isnothing(logfile) 
+        uwerr.(fit.par)
+        flag = typeof(xdata) <: AbstractArray{<:AbstractArray{ADerrors.uwreal}}
+        println(logfile, "juobs.fit_routine output:")
+        if flag
+            println(logfile, "\t\txdata carries uncertainty")
+        end
+        println(logfile, "\t\tGlobal fit with $Nmodel functions")
+        for nm in 1:Nmodel
+            println(logfile,"---Model $(nm)---")
+            println(logfile, "\t\tFit routine in [$(xdata[nm][1].mean),...,$(xdata[nm][end].mean)]")
+        end
+        println(logfile, "\t\tFit parameters:")
+        for i in 1:npar
+            println(logfile,"\t\t\t",fit.par[i])
+        end
+        if flag
+            off = npar
+            println(logfile, "\t\tauxiliary fit parameter:")
+            for nm in 1:Nmodel
+                println(logfile,"---Model $(nm)---")
+                println(logfile, "\t\t\txdata[m][i] \t q[i]")
+                for i in eachindex(xdata[nm])
+                    println(logfile, "\t\t\t",xdata[nm][i]," \t ", fit.par[off+i])
+                end
+                off += npar+length(xdata[nm])
+            end
+        end
+        println(logfile, "\t\tchi2:       ", chi2)
+        println(logfile, "\t\tchiexp:     ", chiexp)
+        println(logfile, "\t\tpvalue:     ", pval)
+    end
+    return fit
+end

src/juobs_tools.jl

@@ -664,31 +664,6 @@ function md_val(a::uwreal, obs::Corr, derm::Vector{Corr}; new_version::Bool=fals
     end
 end
 
-"""  
-    plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, [W::VecOrMat{Float64},] wpm::Union{Dict{Int64},Vector{Float{64}},Dict{String,Vector{Float{64}},Nothing}=nothing} )
-
-It computes plateau average of `obs` in the range `plat`. Effectively this perform a fit to a constant. 
-The field `W` is optional. When absent, `W` is a `Vector{Float64}` constructed from the error of `obs`. 
-When `W` is a `Vector{Float64}` it performs an uncorrelated fit, when `W` is a `Matrix{Float64}`, it perfom a correlated fit
-
-# Potential side effects
-
-When `W` is absent, the function acts with `ADerror.uwerr` on `obs` in order to compute the errors.
-This is the reason why a `wpm` may be given as input.
-
-"""
-function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
-    isnothing(wpm) ? uwerr.(obs) : [uwerr(obs_aux, wpm) for obs_aux in obs]
-    w = 1 ./ err.(obs)[plat[1]:plat[2]].^2
-    av = sum(w .* obs[plat[1]:plat[2]]) / sum(w)
-    return av 
-end
-
-plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, W::Vector{Float64},wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) = sum(W .* obs[plat[1]:plat[2]]) / sum(W)
-
-
-
-
 function model_av(fun::Function, y::Vector{uwreal}, guess::Float64; 
     tm::Vector{Int64}, tM::Vector{Int64}, k::Int64, 
     wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) 
@@ -1175,359 +1150,7 @@ function bayesian_av(fun::Array{Function}, y::Array{uwreal}, tmin_array::Array{I
 
 end
 
-function lin_fit(x::Vector{<:Real}, v::Vector{Float64}, e::Vector{Float64})
-    sig2 = e .* e
-    S = sum(1 ./ sig2)
-    Sx = sum(x ./ sig2)
-    Sy = sum(v ./ sig2)
-    Sxy = sum(v .* x ./ sig2)
-    Sxx = sum(x .* x ./sig2)
-    delta = S * Sxx - Sx*Sx
-    par = [Sxx*Sy-Sx*Sxy, S*Sxy-Sx*Sy] ./delta
-    #C = [[Sxx/delta, -Sx/delta], [-Sx/delta,  S/delta]]
-    return par
-end
-@doc raw"""
-    lin_fit(x::Vector{<:Real}, y::Vector{uwreal})
-
-Computes a linear fit of uwreal data points y. This method return uwreal fit parameters and chisqexpected.
-
-```@example
-fitp, csqexp = lin_fit(phi2, m2)
-m2_phys = fitp[1] + fitp[2] * phi2_phys
-```
-"""
-function lin_fit(x::Vector{<:Real}, y::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
-    isnothing(wpm) ? uwerr.(y) : [uwerr(yaux, wpm) for yaux in y]
-    par = lin_fit(x, value.(y), err.(y))
-    chisq(p, d) = sum((d .- p[1] .- p[2].*x).^2 ./ err.(y) .^2)
-    (fitp, csqexp) = fit_error(chisq, par, y)
-    for i = 1:length(fitp)
-        isnothing(wpm) ? uwerr(fitp[i]) : uwerr(fitp[i], wpm)
-        print("\n Fit parameter: ", i, ": ")
-        details(fitp[i])
-    end
-    println("Chisq / chiexp: ", chisq(par, y), " / ", csqexp, " (dof: ", length(x)-length(par),")")
-    return (fitp, csqexp)
-end
-
-@doc raw"""
-    x_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64})
-
-Computes the results of a linear interpolation/extrapolation in the x axis
-"""
-x_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64}) = (y - par[1]) / par[2]
-@doc raw"""
-    y_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64})
-
-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 """
-    get_chi(model::Function,x::Vector,par,data,W::VecOrMat{Float64})
-
-It computes the chi square of a fit. `model` is assumed to be of the type f(xdata,par). 
-If `W::Vector` the chisquare is computed for an uncorrelated fit
-If `W::Matrix` the chisquare is computed for a correlated fit
-"""
-function get_chi(model::Function,x,par,data,W::Array{Float64,2})
-  return sum((data.-model(x,par))' * W * (data.-model(x,par)))
-end
-
-function get_chi(model::Function,x,par,data,W::Vector{Float64})
-  return sum((data.-model(x,par)).^2 .* W)
-end
-
-@doc raw"""
-        fit_routine(model::Function,xdata::AbstractArray{<:Real},ydata::AbstractArray{ADerrors.uwreal},npar; kwargs...)
-        fit_routine(model::Function,xdata::AbstractArray{uwreal},ydata::AbstractArray{ADerrors.uwreal},npar; kwargs...)
-
-It executes a fit of the `ydata` with the fit function `model`. 
-
-# Arguments:
-
-  - `model`: Fit function used to fit the data. the function assumes that its called as `model(xdata,parameters)`.
-  - `xdata`: indipendent variable in the fit. It can be an `AbstractArray{<:Real}`, that is the "default" method, or 
-             a `AbstractArray{uwreal}`. In the latter case, the method builds the function
-``math
-F(q_i) = \begin{cases}
-         model(q_{npar+i},q[1:npar])  \quad \text{if} i<= length(xdata)\\
-         q_i \quad \text{elsewhere}
-       \end{cases}
-``
-             and then call the default method with dummy `xdata` and `vcat(ydata,xdata)` as new `ydata`
-  - `ydata`: Data variable that carries the error 
-  - `npar`: number of parameters used in the fit.
-
-# Keyword parameters:
-
-  - `wpm`: window parameter used by ADerrors to computes the errors. Default to `Dict{Int64,Vector{Float64}}()`.
-           This is used both when computing the error of the data and when computing the covariance matrix through `ADerrors.cov`
-           It accept a dictionary with ensemble id as keys (either `String` or `Int64`) and `Vector{Float64}` as value.
-           See also ADerrors documentation.
-
-  - `W::AbstractVecOrMat{Float64}`: Weight matrix/vector. If `W` is a `Vector` the fit is perfomed uncorrelated ,otherwise is correlated.
-                                    default to `Vector{Float64}()`
-
-  - `corr::Bool`: flag used to generated `W` when not given. Default to `true`. If `false`, `W = [1/y.err^2 for y in ydata]`.
-                  If `true` the covariance matrix of the data is used to compute `W` as `W = C^{-1}`.
-
-  - `C::AbstractMatrix`: covariance matrix used to compute W. Default to `Matrix{Float64}(undef, 0, 0)`.
-                         If `W` is not given, and `corr` is `true`, then it is used to compute `W`.
-                         If need to computed `W` but not given,  then `C = ADerrors.cov(ydata,wpm)`. See also ADerrors documentation.
-                         If available, it will be used to better estimate the fit pvalue and the expected chi-square,
-                         but it will not be computed if not available at this point. 
-
-  - `guess::Vector{Float64}`: Initial guess for the parameter used in the fit routine. If not given, default to `fill(0.5,npar)`.
-                              If `xdata isa Vector{uwreal}`, the mean value of `xdata` are used as guess for the relevant `q_i` parameters.
-  
-  - `logfile`: handle to a logfile. If `nothing`, no log will be written. It can be `IO` file or a custom log structure that has
-               an overloading for `print()` and `println()`.
-
-# Returns
-It returns a `NamedTuple` with names:
-  - `:par`:    fit parameter as `Vector{uwreal}`
-  - `:chi2`:   chisquare,
-  - `:chiexp`: chisquare expected 
-  - `:pval`:   pvalue 
-"""
-function fit_routine(model::Function,
-                     xdata::AbstractArray{<:Real},
-                     ydata::AbstractArray{uwreal},
-                     npar::Int64; 
-                     wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}}=Dict{Int64,Vector{Float64}}(),
-                     corr::Bool = true, 
-                     W::AbstractArray{Float64}=Vector{Float64}(),
-                     guess::AbstractVector{Float64} = fill(0.5,npar),
-                     logfile = nothing,
-                     C::AbstractArray{Float64} = Matrix{Float64}(undef,0,0));
-
-    nobs = length(ydata)
-    if length(xdata)!=nobs
-        error("xdata and ydata must be the of same size")
-    end
-
-    [uwerr(y, wpm) for y in ydata]
-    if length(W) ==0
-        if corr 
-            C = length(C) ==0 ? GammaMethod.cov_AD(ydata,wpm) : C
-            W = LinearAlgebra.pinv(C);
-            W = (W'+W)*0.5
-        else
-            W = [1/y.err^2 for y in ydata]
-        end
-    end
-
-    fit = LsqFit.curve_fit(model,xdata,ADerrors.value.(ydata),W,guess)
-
-    chisq(p,d) = get_chi(model,xdata,p,d,W)
-    chi2 = sum(fit.resid.^2)
-    par = fit_error(chisq,LsqFit.coef(fit),ydata,wpm,W,false)
-    # if the fit is correlated  we have already C, no need to recompute it.
-    chiexp,pval =  if length(C) !=0 
-        GammaMethod.chiexp(chisq,LsqFit.coef(fit),value.(ydata),C,W),juobs.pvalue(chisq,chi2,LsqFit.coef(fit),ydata,W,C=C)
-    else  
-        ADerrors.chiexp(chisq,LsqFit.coef(fit),ydata,wpm,W=W),juobs.pvalue(chisq,chi2,LsqFit.coef(fit),ydata,W)
-    end
-
-    if !isnothing(logfile)
-        uwerr.(par);
-        println(logfile, "juobs.fit_routine output:")
-        println(logfile, "\t\tFit routine in [$(xdata[1]),...,$(xdata[end])]")
-        println(logfile, "\t\tparameters :")
-        for p in par
-            println(logfile,"\t\t\t",p)
-        end
-        println(logfile, "\t\tchi2:       ", chi2)
-        println(logfile, "\t\tchiexp:     ", chiexp)
-        println(logfile, "\t\tpvalue:     ", pval)
-    end
-
-    return (par=par,chi2=chi2,chiexp=chiexp,pval=pval)
-end
-
-"""
-        check_size(M,n)
-
-Check if `M` has sizes all equal to n or if is empty
-"""
-function check_sizes(M,n)
-    if length(M) ==0 || all(size(M) .== n)
-        return true
-    end
-    return false
-end
-
-function fit_routine(model::Function,
-                     xdata::AbstractArray{ADerrors.uwreal},
-                     ydata::AbstractArray{ADerrors.uwreal},
-                     npar::Int64;
-                     W::AbstractArray = Float64[],
-                     wpm::AbstractDict = Dict{Int64,Vector{Float64}}(),
-                     corr::Bool = true,
-                     guess::AbstractArray = fill(0.5,npar),
-                     logfile = nothing,
-                     C::AbstractArray = Float64[])
-    Nx = length(xdata)
-    if (Nx != length(ydata))
-        error("[juobs] Error: xdata and ydata must have the same dimension")
-    end
-    Ndata = 2Nx; ## in total we have 2*lengt(xdata) data points
-    data = vcat(ydata,xdata)
-    if !check_sizes(W,Ndata)
-        error("[juobs] Error: W does not have the correct sizes")
-    end
-    if !check_sizes(C,Ndata)
-        error("[juobs] Error: C does not have the correct sizes")
-    end
-
-    if corr
-        C = length(C) ==0 ? GammaMethod.cov_AD(data) : C
-        W = length(W) ==0 ? LinearAlgebra.pinv(W) : W
-    else
-        W = length(W) ==0 ? [1/yy.err^2 for yy in data] : W
-    end
-
-    f(q,p) = [i<=Nx ? model(q[i],p) : q[i-Nx] for i in 1:Ndata]
-    fit_func(x,p) = f(view(p,npar+1:npar+Nx),view(p,1:npar))
-    new_guess = [guess; value.(xdata)]
-    fit = fit_routine(fit_func,zeros(Ndata),data,npar+Nx,
-                      guess=new_guess, wpm=wpm,corr=corr,
-                      W=W,C=C,logfile=nothing)
-
-    if !isnothing(logfile) 
-        uwerr.(fit.par)
-        println(logfile, "juobs.fit_routine output:")
-        println(logfile, "\t\txdata carries uncertainty")
-        println(logfile, "\t\tFit routine in [$(xdata[1].mean),...,$(xdata[end].mean)]")
-        println(logfile, "\t\tFit parameters:")
-        for i in 1:npar
-            println(logfile,"\t\t\t",fit.par[i])
-        end
-        println(logfile, "\t\tauxiliary fit parameter:")
-        println(logfile, "\t\t\txdata[i] \t q[i]")
-        for i in 1:Nx
-            println(logfile, "\t\t\t",xdata[i]," \t ", fit.par[npar+i])
-        end
-        println(logfile, "\t\tchi2:       ", fit.chi2)
-        println(logfile, "\t\tchiexp:     ", fit.chiexp)
-        println(logfile, "\t\tpvalue:     ", fit.pval)
-    end
-    return fit
-end
-
-
-@doc raw"""
-        fit_routine(model::AbstractArray{Function},
-                    xdata::AbstractArray{<:AbstractArray},
-                    ydata::AbstractArray{<:AbstractArray{ADerrors.uwreal}},npar; kwargs...)
-       
-
-It executes a combined fit of the `ydata` with the fit functions in `model`. The method define the general model `F(X_{mi},p) = model[m](xdata[i],p)`
-and call `fit_routine(F,vcat(xdata),vcat(ydata),npar;kwargs....)`. The `kwargs` parameters are accordingly updated to reflect the new function.
-
-# Returns
-It returns a `NamedTuple` with names:
-  - `:par`:    fit parameter as `Vector{uwreal}`
-  - `:chi2`:   chisquare,
-  - `:chiexp`: chisquare expected 
-  - `:pval`:   pvalue 
-"""
-function fit_routine(models::AbstractArray{Function},
-                     xdata::AbstractArray{<:AbstractArray},
-                     ydata::AbstractArray{<:AbstractArray{ADerrors.uwreal}},
-                     npar::Int64;
-                     W::AbstractArray = Float64[],
-                     C::AbstractArray = Float64[],
-                     wpm::AbstractDict = Dict{Int64,Vector{Float64}}(),
-                     corr::Bool = true,
-                     guess::AbstractArray{<:AbstractArray} = [fill(0.5,npar[i]) for i in eachindex(models)],
-                     logfile = nothing,
-                     )
-    Nmodel = length(models)
-    if Nmodel != length(xdata) != length(ydata)
-        error("[juobs] Error: you need the same number of model, xdata and ydata")
-    end
-    if !(length(W) == 0 || length(W)==Nmodel)
-        error("[juobs] Error: You need to pass a W matrix for model or none")
-    end
-    if !(length(C) == 0 || length(C)==Nmodel)
-        error("[juobs] Error: You need to pass a C matrix for model or none")
-    end
-    if length(guess)!=Nmodel
-        error("[juobs] Error: You need to specify the inital guess for all model or for none")
-    end
-    ndata = length.(xdata)
-    if !all(length.(ydata).== ndata)
-        error("[juobs] Error: Mismatch in xdata and ydata. Make sure that the data corresponds")
-    end
-    Ntotal = sum(ndata);
-    X = vcat(xdata...)
-    Y = vcat(ydata...)
-
-    function make_matrix(M)
-        aux = zeros(Ntotal,Ntotal)
-        ie = 0
-        for m in 1:Nmodel
-            is = ie+1
-            ie += ndata[m];
-            aux[is:ie,is:ie] .= M[m]
-        end
-        return aux
-    end
-
-    WW = length(W)==0 ? W : make_matrix(W);
-    CC = length(C)==0 ? C : make_matrix(C);
 
-    function Model(x,p)
-        res = zeros(Nmodel)
-        res[1] = models[1](x[1:ndata[1]],p)
-        for i in 2:Nmodel
-            rd = (ndata[i-1]+1):ndata[i]
-            rp = (npar[i-1]+1):npar[i]
-            res[i] = models[i](x[rd],p)
-        end
-        return res
-    end
-    Guess = vcat(guess...)
-    fit = fit_routine(Model,X,Y,npar,guess = Guess,W=WW,C=CC,corr=corr,logfile=nothing,wpm=wpm)
-
-    if !isnothing(logfile) 
-        uwerr.(fit.par)
-        flag = typeof(xdata) <: AbstractArray{<:AbstractArray{ADerrors.uwreal}}
-        println(logfile, "juobs.fit_routine output:")
-        if flag
-            println(logfile, "\t\txdata carries uncertainty")
-        end
-        println(logfile, "\t\tGlobal fit with $Nmodel functions")
-        for nm in 1:Nmodel
-            println(logfile,"---Model $(nm)---")
-            println(logfile, "\t\tFit routine in [$(xdata[nm][1].mean),...,$(xdata[nm][end].mean)]")
-        end
-        println(logfile, "\t\tFit parameters:")
-        for i in 1:npar
-            println(logfile,"\t\t\t",fit.par[i])
-        end
-        if flag
-            off = npar
-            pritnln(logfile, "\t\tauxiliary fit parameter:")
-            for nm in 1:Nmodel
-                println(logfile,"---Model $(nm)---")
-                println(logfile, "\t\t\txdata[m][i] \t q[i]")
-                for i in eachindex(xdata[nm])
-                    println(logfile, "\t\t\t",xdata[nm][i]," \t ", fit.par[off+i])
-                end
-                off += npar+length(xdata[nm])
-            end
-        end
-        println(logfile, "\t\tchi2:       ", chi2)
-        println(logfile, "\t\tchiexp:     ", chiexp)
-        println(logfile, "\t\tpvalue:     ", pval)
-    end
-    return fit
-end
 
 function fve(mpi::uwreal, mk::uwreal, fpi::uwreal, fk::uwreal, ens::EnsInfo)
     mm = [6,12,8,6,24,24,0,12,30,24,24,8,24,48,0,6,48,36,24,24]