some changes

constants/juobs_uwreal_const.jl

@@ -132,7 +132,7 @@ let C = zeros(5, 5)
     if !isassigned(Za)
       Za.x = ADerrors.cobs(ZA_data, C, "ZA")
     end
-    return za.x[beta_to_idx[beta]];
+    return Za.x[beta_to_idx[beta]];
   end
 end
 
@@ -150,7 +150,7 @@ let C = zeros(5, 5)
     if !isassigned(Zv)
       Zv.x = ADerrors.cobs(ZV_data, C, "ZV")
     end
-    return zv.x[beta_to_idx[beta]];
+    return Zv.x[beta_to_idx[beta]];
   end
 
 end
@@ -217,9 +217,9 @@ let C = zeros(5, 5)
   end
   ZZ = Ref{Vector{ADerrors.uwreal}}();
   global function Z(beta::Float64)
-    if !isassinged(ZZ)
+    if !isassigned(ZZ)
       ZZ.x =ADerrors.cobs(ZDATA, C, "Z")
     end
-    return ZZ[beta_to_idx[beta]];
+    return ZZ.x[beta_to_idx[beta]];
   end
 end

src/juobs_gammaMethod.jl

@@ -286,7 +286,7 @@ module GammaMethod
   cov_min(vobs::Vector{ADerrors.uwreal}) = cov_min(vobs,ADerrors.wsg) 
 
   @doc raw""" 
-      cov_AD(vobs::Vector{ADerrors.uwrela} [ws::ADerrors.wspace])
+      cov_AD(vobs::Vector{ADerrors.uwrela} [ws::ADerrors.wspace]; biased::Bool=true, info::Bool=false)
 
   It computes the unbiased covariance matrix of `vobs` "à la ADerrors".
 
@@ -297,9 +297,12 @@ module GammaMethod
   This method is not stable, in the sense that the covarince between two `ADerrors.uwreal` 
   depend also on the other `ADerrors.uwreal` included in `vobs`
   
+  # KEYWORD ARGUMENTS
+  - `bias::Bool`: if `true` (default), it computes the covariance matrix including the bias as `ADerrors`
+  - `info::Bool`: if `true` return the integration window used to compute the covariance matrix as a Dict{Int64,Int64}(id=> iw)
   See also [`cov_min`](@ref) [`cov_pe`](@ref)
   """
-  function cov_AD(vobs::Vector{ADerrors.uwreal},ws::ADerrors.wspace; biased::Bool=false, info::Bool=false)
+  function cov_AD(vobs::Vector{ADerrors.uwreal},ws::ADerrors.wspace; biased::Bool=true, info::Bool=false)
     ids = ADerrors.unique_ids_multi(vobs, ws)
     nid = length(ids)
     iw  = zeros(Int64,nid) 
@@ -342,7 +345,7 @@ module GammaMethod
     for n1 in 1:nobs-1, n2 in n1+1:nobs
       cov[n2,n1] = cov[n1,n2]
     end
-    return info ? (cov, iw) : cov
+    return info ? (cov, Dict(ids.=>iw)) : cov
   end
 
   cov_AD(vobs::Vector{ADerrors.uwreal};biased::Bool=false) = cov_AD(vobs,ADerrors.wsg,biased=biased) 

src/juobs_tools.jl

@@ -1235,18 +1235,30 @@ function get_chi(model::Function,x::Vector,par,data,W::Vector{Float64})
 end
 
 @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),logfile);
+    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
+    C::AbstractMatrix{Float64} = Matrix{Float64}());
 
 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. 
+  - `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::AbstractVecOrMat{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()`.
-
+  - `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()`.
+  - `C::AbstractMatrix`: covariance matrix used to compute W. 
+     If W is computed by using a different covariance matrix that what given by ADerrors, it is advised to pass it to the function
+     to have a better chiexp and pvalue. If W and C are not given but the fit is correlated, then C is computed using ADerrors
 
 # Returns
 It returns a NamedTuple with names:
@@ -1261,9 +1273,10 @@ function fit_routine(model::Function,
   npar::Int64; 
   wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}}=Dict{Int64,Vector{Float64}}(),
   corr::Bool = true, 
-  W::VecOrMat{Float64}=Vector{Float64}(),
+  W::AbstractVecOrMat{Float64}=Vector{Float64}(),
   guess::Vector{Float64} = fill(0.5,npar),
-  logfile = nothing,);
+  logfile = nothing,
+  C::AbstractMatrix{Float64} = Matrix{Float64}());
 
   nobs = length(ydata)
   if length(xdata)!=nobs
@@ -1273,8 +1286,8 @@ function fit_routine(model::Function,
   [uwerr(y, wpm) for y in ydata]
   if length(W) ==0
     if corr 
-      C = ADerrors.cov(ydata,wpm)
-      W = LinearAlgebra.pinv(W);
+      C = length(C) ==0 ? ADerrors.cov(ydata,wpm) : C
+      W = LinearAlgebra.pinv(C);
       W = (W'+W)*0.5
     else
       W = [1/y.err^2 for y in ydata]
@@ -1286,13 +1299,14 @@ function fit_routine(model::Function,
   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)
-
-  chiexp = ADerrors.chiexp(chisq,LsqFit.coef(fit),ydata,wpm,W=W)
-
-  pval = juobs.pvalue(chisq,chi2,LsqFit.coef(fit),ydata,W)
+ # if the fit is 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)
-
     println(logfile, "juobs.fit_routine output:") 
     println(logfile, "\t\tFit routine in [$(xdata[1]),...,$(xdata[2])]")
     println(logfile, "\t\tparameters :", par...)