bla

constants/juobs_uwreal_const.jl

@@ -5,7 +5,7 @@ let
   C = zeros(8, 8)
   M_values = [1869.65, 2010.26, 1968.34, 2112.2, 2980.3, 3096.916, 5279.65, 5366.3] #MD, MD*, MDs, MDs*, \eta_c, J/\psi , MB, MBs(MeV)
   M_error = [0.05, 0.05, 0.07, 0.4, 1.2, 0.011, 0.12, 0.6]
-  M_ = Ref{Vector{uwreal}}()
+  M_ = Ref{Vector{ADerrors.uwreal}}()
   for i = 1:8
     C[i, i] = M_error[i] ^ 2
   end
@@ -20,7 +20,7 @@ end
 let 
   #1802.05243 taking into account correlations in zm_tm
   C = [0.375369e-6 0.429197e-6 -0.186896e-5; 0.429197e-6 0.268393e-4 0.686776e-4; -0.186896e-5 0.686776e-4 0.212386e-3]
-  z = Ref{Vector{uwreal}}();
+  z = Ref{Vector{ADerrors.uwreal}}();
   global function ZP(beta::Float64)
     if !isassigned(z)
       z.x =  ADerrors.cobs([0.348629, 0.020921, 0.070613], C, "z")
@@ -28,8 +28,8 @@ let
     return z.x[1] + z.x[2] *(beta-3.79) + z.x[3]*(beta-3.79)^2
   end
   CC = [0.164635e-4 0.215658e-4 -0.754203e-4; 0.215658e-4 0.121072e-2 0.308890e-2; -0.754203e-4 0.308890e-2 0.953843e-2]
-  zm = Ref{Vector{uwreal}}()
-  Mrat_ = Ref{uwreal}()
+  zm = Ref{Vector{ADerrors.uwreal}}()
+  Mrat_ = Ref{ADerrors.uwreal}()
 
   global function Mrat()
     if !isassigned(Mrat_)
@@ -61,7 +61,7 @@ let C = zeros(5, 5)
   for i in 1:5
     C[i, i] = ZM_error[i] ^ 2
   end
-  ZM = Ref{Vector{uwreal}}()
+  ZM = Ref{Vector{ADerrors.uwreal}}()
   global function zm(beta::Float64)
     if !isassigned(ZM)
       ZM.x =  ADerrors.cobs(ZM_data, C, "ZM values")
@@ -77,7 +77,7 @@ let C = zeros(5, 5)
   for i in 1:5
     C[i, i] = ZM_tm_error[i] ^ 2
   end
-  ZM_tm = Ref{Vector{uwreal}}()
+  ZM_tm = Ref{Vector{ADerrors.uwreal}}()
   
   global function zm_tm(beta::Float64)
     if !isassigned(ZM_tm)
@@ -96,7 +96,7 @@ let C = zeros(5, 5)
     C[i, i] = t0_error[i] ^ 2
   end
   
-  t0_ = Ref{Vector{uwreal}}()
+  t0_ = Ref{Vector{ADerrors.uwreal}}()
   global function t0(beta::Float64)
     if !isassigned(beta)
       t0_.x = ADerrors.cobs(t0_data, C, "t0")
@@ -104,7 +104,7 @@ let C = zeros(5, 5)
     return t0_.x[beta_to_idx[beta]];
   end
 
-  t0_ph_ = Ref{uwreal}()
+  t0_ph_ = Ref{ADerrors.uwreal}()
   global function t0_ph()
     if !isassigned(t0_ph)
       t0_ph_.x = ADerrors.uwreal([0.4118,0.25e-4],"sqrt(8 t0) (fm)") 
@@ -126,7 +126,7 @@ let C = zeros(5, 5)
   for i in 1:5
     C[i, i] = ZA_err[i] ^ 2
   end
-  Za = Ref{Vector{uwreal}}()
+  Za = Ref{Vector{ADerrors.uwreal}}()
 
   global function za(beta::Float64)
     if !isassigned(Za)
@@ -144,7 +144,7 @@ let C = zeros(5, 5)
   for i in 1:5
     C[i, i] = ZV_err[i] ^ 2 
   end
-  Zv = Ref{Vector{uwreal}}()
+  Zv = Ref{Vector{ADerrors.uwreal}}()
 
   global function zv(beta::Float64)
     if !isassigned(Zv)
@@ -164,7 +164,7 @@ let C = zeros(5, 5)
   for i in 1:5
     C[i, i] = ZS_over_ZP_err[i] ^ 2
   end
-  ZS_over_ZP = Ref{Vector{uwreal}}()
+  ZS_over_ZP = Ref{Vector{ADerrors.uwreal}}()
   global function zs_over_zp(beta::Float64)
     if !isassigned(ZS_over_ZP_data)
       ZS_over_ZP.x = ADerrors.cobs(ZS_over_ZP_data, C, "ZS/ZP")
@@ -181,7 +181,7 @@ let C = zeros(5, 5)
   for i in 1:5
     C[i, i] = RM_err[i] ^2 
   end
-  RM  = Ref{Vector{uwreal}}()
+  RM  = Ref{Vector{ADerrors.uwreal}}()
   global function rm(beta::Float64)
     if !isassigned(rm)
       RM.x = ADerrors.cobs(RM_data, C, "rm")
@@ -199,7 +199,7 @@ let C = zeros(5, 5)
   for i in 1:5
     C[i,i]  = BA_MINUS_BP_err[i] ^2
   end
-  BA_BP = Ref{Vector{uwreal}}()
+  BA_BP = Ref{Vector{ADerrors.uwreal}}()
   global function ba_bp(beta::Float64)
     if !isassigned(BA_BP)
       BA_BP.x =  ADerrors.cobs(BA_MINUS_BP_data, C, "ba-bp")
@@ -215,7 +215,7 @@ let C = zeros(5, 5)
   for i in 1:5
     C[i,i]  = ZERR[i] ^2
   end
-  ZZ = Ref{Vector{uwreal}}();
+  ZZ = Ref{Vector{ADerrors.uwreal}}();
   global function Z(beta::Float64)
     if !isassinged(ZZ)
       ZZ.x =ADerrors.cobs(ZDATA, C, "Z")

src/juobs_tools.jl

@@ -1218,85 +1218,83 @@ 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, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
-
-    fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, covar::Bool=false)
-
-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.
-The flag `wpm` is an optional array of Float64 of lenght 4. The first three paramenters specify the criteria to determine
-the summation windows:
 
-- `vp[1]`: The autocorrelation function is summed up to ``t = round(vp[1])``.
+@doc """
+    get_chi(model::Function,x::Vector,par,data,W::VecOrMat{Float64})
 
-- `vp[2]`: The sumation window is determined using U. Wolff poposal with ``S_\tau = wpm[2]``
+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::Vector,par,data,W::Array{Float64,2})
+  return sum((data.-model(x,par))' * W * (data.-model(x,par)))
+end
 
-- `vp[3]`: The autocorrelation function ``\Gamma(t)`` is summed up a point where its error ``\delta\Gamma(t)`` is a factor `vp[3]` times larger than the signal.
+function get_chi(model::Function,x::Vector,par,data,W::Vector{Float64})
+  return sum((data.-model(x,par)).^2 .* W)
+end
 
-An additional fourth parameter `vp[4]`, tells ADerrors to add a tail to the error with ``\tau_{exp} = wpm[4]``.
-Negative values of `wpm[1:4]` are ignored and only one component of `wpm[1:3]` needs to be positive.
-If the flag `covar`is set to true, `fit_routine` takes into account covariances between x and y for each data point.
-```@example
-@. model(x,p) = p[1] + p[2] * exp(-(p[3]-p[1])*x)
-@. model2(x,p) = p[1] + p[2] * x[:, 1] + (p[3] + p[4] * x[:, 1]) * x[:, 2]
-fit_routine(model, xdata, ydata, param=3)
-fit_routine(model, xdata, ydata, param=3, covar=true)
-```
+@doc raw"""
+    fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, npar::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}}=Dict{Int64,Vector{Float64}}(),corr::Bool = true, W::VecOrMat{Float64}=Vector{Float64}(),guess::Vector{Float64} = fill(0.5,npar));
+
+It executes a fit of the `ydata` with the fit function `model`.
+
+# parameters
+  - `npar::Int64`: number of fit parameters.
+  - `wpm`: Windows parameter that ADerrors uses to computes the errors. This is used both when computing the error of the data and when computing the covariance matrix through `ADerrors.cov`
+  - `corr::Bool`: if `true` (default) it exectues the correlated fit, this flag is overridden if `W` is passed as a `Vector`. When `corr=true`, and `W` is not given, it computes the covarinace matrix through `ADerrors.cov`, then it is inverted and symmetrized. 
+  - `W::VecOrMat{Float64}`: Weight matrix/vector. When given, the flag `corr` is overridden. If `W` is a `Vector` the fit is perfomed uncorrelated,otherwise is correlated. 
+  - `guess::Vector{Float64}`: Initial guess used in the fit routine.
+  - `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::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; info::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
-    inv_cov::Union{Matrix{Float64}, Nothing}=nothing)
-
-    isnothing(wpm) ? uwerr.(ydata) : [uwerr(yaux, wpm) for yaux in ydata]
-    
-    yval = value.(ydata)
-    yer = err.(ydata)
-    
-    ## CODE HERE THE UPDATED VERSION OF CORRELATED FITS
-    if isnothing(inv_cov)
-        @info("Uncorrelated fit")
-        chisq = gen_chisq(model, xdata, yer)
-        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)
+function fit_routine(model::Function,
+  xdata::Array{<:Real},
+  ydata::Array{uwreal}, 
+  npar::Int64; 
+  wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}}=Dict{Int64,Vector{Float64}}(),
+  corr::Bool = true, 
+  W::VecOrMat{Float64}=Vector{Float64}(),
+  guess::Vector{Float64} = fill(0.5,npar),
+  logfile = nothing);
+
+  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 
+      W = ADerrors.cov(ydata,wpm)
+      W = LinearAlgebra.pinv(W);
+      W = (W'+W)*0.5
     else
-        @info("Correlated fit")
-        #covar_inv = inv(get_covariance(ydata)) 
-        #covar_inv = (covar_inv + covar_inv')/2
-        #covar_inv = inv(Symmetric(get_covariance(ydata)))
-        chisq(par, dat) = gen_chisq_correlated(model, xdata, covar_inv, par, dat)
-        fit = curve_fit(model, xdata, yval, covar_inv, fill(0.5, param))        
-        #println(chisq(coef(fit), ydata))
-        (upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata, W=covar_inv) : fit_error(chisq, coef(fit), ydata, wpm, W=covar_inv)
-    end
-    chi2_fit_res = sum(fit.resid.^2 )
-    # compute and print single point contribution to chi2
-    # println("\n")
-    # for i in 1:length(fit.resid)
-        # println((fit.resid[i])^2)
-    # end
-    # println("\n")
-
-    # println("chi2 from fit residual = ", chi2_fit_res)
-    # println("chi2 from chi2 function = ", chisq(coef(fit), ydata))
-
-
-    #Info
-    if info
-        for i = 1:length(upar)
-            isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm)
-            print("\n Fit parameter: ", i, ": ")
-            details(upar[i])
-        end
+      W = [1/y.err^2 for y in ydata]
     end
-    #chis2_corrected = (length(yval) - param) * chisq(coef(fit), ydata) / chi_exp
-    chis2_corrected = (length(yval) - param) * chi2_fit_res / chi_exp
+  end
 
-    #println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", length(yval) - param,")")
-    println("Chisq / chiexp: ", chi2_fit_res, " / ", chi_exp, " (dof: ", length(yval) - param,")")
-    println("Chisq corrected: ", chis2_corrected)
-    println("Return: params, chi2, chi2exp")
-    return upar, chi2_fit_res,  chi_exp
+  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,chiexp = fit_error(chisq,LsqFit.coef(fit),ydata,wpm,W)
+  pval = juobs.pvalue(chisq,chi2,LsqFit.coef(fit),ydata,W)
+
+  if !isnothing(logfile) 
+    println(logfile, "Fit routine in [$(xdata[1]),...,$(xdata[2])]")
+    println(logifle, "parameters :", par...)
+    println(logfile, "chi2:       ", chi2)
+    println(logfile, "chiexp:     ", chiexp)
+    println(logfile, "pvalue:     ", pval)
+  end
+
+  return (par=par,chi2=chi2,chiexp=chiexp,pval=pval)
 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,