Merge branch 'alessandro'

constants/juobs_const.jl

 using ADerrors
 #PDG
 const hc = 197.3269804 #MeV fm
-const M_values = [1869.65, 2010.26, 1968.34, 2112.2, 2980.3, 3096.916] #MD, MD*, MDs, MDs*, \eta_c, J/\psi (MeV)
-const M_error = [0.05, 0.05, 0.07, 0.4, 1.2, 0.011]
+const M_values = [1869.65, 2010.26, 1968.34, 2112.2, 2980.3, 3096.916, 5279.65] #MD, MD*, MDs, MDs*, \eta_c, J/\psi , MB (MeV)
+const M_error = [0.05, 0.05, 0.07, 0.4, 1.2, 0.011, 0.12]
 #1802.05243 
 const b_values = [3.40, 3.46, 3.55, 3.70, 3.85]
 const ZM_data = [1.9684, 1.9935, 2.0253, 2.0630, 2.0814]
@@ -12,6 +12,11 @@ const ZM_tm_error = [42, 33, 42, 48, 55] .* 1e-4
 # 1808.09236
 const ZA_data = [0.75642, 0.76169, 0.76979, 0.78378, 0.79667]
 const ZA_err = [72, 93, 43, 47, 47] .*1e-5
+const ZV_data = [0.72259, 0.72845, 0.73886, 0.75574, 0.77074]
+const ZV_err = [74, 89, 46, 48, 49] .* 1e-5
+#2005.01352
+const ZS_over_ZP_data = [1.3497, 1.2914, 1.2317, 1.1709, 1.1343]
+const ZS_over_ZP_err = [83, 64, 48, 23, 25 ] .* 1e-4
 #1608.08900
 const b_values2 = [3.40, 3.46, 3.55, 3.70]
 const t0_data = [2.86, 3.659, 5.164, 8.595]
@@ -19,18 +24,31 @@ const t0_error = [11, 16, 18, 29] .* 1e-3
 const t0_ph_value = [0.413]
 const t0_ph_error = ones(1,1) .* 5e-3
 
+## 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 = cobs([0.348629, 0.020921, 0.070613], C, "z")
+ZP(b) = z[1] + z[2] * (b - 3.79) + z[3] * (b - 3.79)^2
+Mrat = uwreal([.9148, 0.0088], "M/mhad")
+zm_tm_covar(b) = Mrat / ZP(b)
+##
 #Cov matrices
 const C1 = zeros(5, 5)
 const C2 = zeros(5, 5)
 const C3 = zeros(4, 4)
-const C4 = zeros(6, 6)
+const C4 = zeros(7, 7)
 const C5 = zeros(5, 5)
-for i = 1:6
+const C6 = zeros(5, 5)
+const C7 = zeros(5, 5)
+
+for i = 1:7
     C4[i, i] = M_error[i] ^ 2
     if i<=5
         C1[i, i] = ZM_error[i] ^ 2
         C2[i, i] = ZM_tm_error[i] ^ 2
-        C5[i, i] = ZA_err[i]^2
+        C5[i, i] = ZA_err[i] ^ 2
+        C6[i, i] = ZS_over_ZP_err[i] ^ 2
+        C7[i, i] = ZV_err[i] ^ 2        
         if i <= 4
             C3[i, i] = t0_error[i] ^ 2
         end
@@ -44,9 +62,13 @@ const t0_ph = cobs(t0_ph_value, t0_ph_error .^ 2, "sqrt(8 t0) (fm)")
 const a_ = t0_ph ./ sqrt.(8 .* t0_)
 const M = cobs(M_values, C4, "charmed meson masses")
 const Za = cobs(ZA_data, C5, "ZA")
+const Zv = cobs(ZV_data, C7, "ZV")
+const ZS_over_ZP = cobs(ZS_over_ZP_data, C6, "ZS over ZP")
 
 zm(beta::Float64) = ZM[b_values .== beta][1]
 zm_tm(beta::Float64) = ZM_tm[b_values .== beta][1]
 t0(beta::Float64) = t0_[b_values2 .== beta][1]
 a(beta::Float64) = a_[b_values2 .== beta][1]
-za(beta::Float64) = Za[b_values .== beta][1]
\ No newline at end of file
+za(beta::Float64) = Za[b_values .== beta][1]
+zv(beta::Float64) = Zv[b_values .== beta][1]
+zs_over_zp(beta::Float64) = ZS_over_ZP[b_values .== beta][1]
\ No newline at end of file

gevp tools/tools.jl

@@ -138,11 +138,11 @@ function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.C
         evec = Array{Array{}}(undef, length(aux_mat))
         for i = 1:length(aux_mat)
             mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
-            evalues, eigvec = uwgevp_tot(mat_obs, tnot, evec=true, iter=50)
-            res_en[i] = infoEn(energies(evalues), aux_mat[i])
+            eigvec = uwgevp_tot(mat_obs, delta_t=3)
+            #res_en[i] = infoEn(energies(evalues), aux_mat[i])
             evec[i] = eigvec
         end
-        return res_en, evec
+        return  evec
     end
 end
 function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot13::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, tot33::Array{Array{juobs.Corr}}, tot23::Array{Array{juobs.Corr}}, evec::Bool=false; tnot::Int64=2)
@@ -160,12 +160,12 @@ function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.C
         evec = Array{Array{}}(undef, length(aux_mat))
         for i = 1:length(aux_mat)
             mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
-            evalues, eigvec = uwgevp_tot(mat_obs, tnot, evec=true)
+            eigvec = uwgevp_tot(mat_obs, delta_t=3)
             println("in comp_energy, i is = ", i)
-            res_en[i] = infoEn(energies(evalues), aux_mat[i])
+            #res_en[i] = infoEn(energies(evalues), aux_mat[i])
             evec[i] = eigvec
         end
-        return res_en, evec
+        return  evec
     end
 end
 mutable struct infoEn
@@ -180,13 +180,22 @@ mutable struct infoEn
     end
 end
 function pseudo_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal, mu::Array{Float64})
-    Rn = [exp(mass * (t)/2) * evec[t][1] for t =2:length(evec)]
+    Rn = [exp(mass * (t)/2) * evec[t][1] for t =1:length(evec)]
     aux = sqrt(2)*sum(mu)  / mass^1.5 
     return uwdot(aux,Rn)
 end
-function vec_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal)
-    Rn = [exp(mass * t/2) * evec[t][1] for t =1:length(evec)]
-    return uwdot(Rn , 1/ mass)
+function vec_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal, deg::Bool)
+    if !deg
+        Rn = [exp(mass * t/2) * evec[t][1] for t =1:length(evec)]
+        return uwdot(Rn , za(beta[iens]) * sqrt(2/ mass))
+    else
+        Rn = [exp(mass * t/2) * evec[t][4] for t =1:length(evec)]
+        println(za(beta[iens]))
+        println(mass)
+        return uwdot(Rn , za(beta[iens]) * sqrt(2/ mass))
+        #return uwdot(Rn , 1/mass)# * sqrt(2/ mass))
+end
+
 end
 function extract_dec_const(mat_elem::uwreal, mass::uwreal, mu::Array{Float64})
     return sqrt(2)*sum(mu)*mat_elem / mass^1.5

src/.vscode/settings.json

+{}
\ No newline at end of file

src/juobs.jl

 module juobs
+using Base: Float64
 using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, LeastSquaresOptim
 import Statistics: mean
 
@@ -9,8 +10,8 @@ include("juobs_tools.jl")
 include("juobs_obs.jl")
 
 export read_mesons, read_ms1, read_ms, read_md, truncate_data!
-export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs
-export corr_obs, corr_sym, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
+export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs, hess_reduce, uwcholesky, transpose, tridiag_reduction, make_positive_def, invert_covar_matrix
+export corr_obs, corr_sym, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine, global_fit_routine
 export meff, mpcac, dec_const, dec_const_pcvc, comp_t0
 
 end # module

src/juobs_linalg.jl

@@ -126,7 +126,8 @@ function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}} }; wp
             ratio = evals[t][i] / evals[t+1][i]
             aux_en[t] = 0.5*log(ratio * ratio)
         end
-        isnothing(wpm) ? uwerr.(aux_en) : uwerr.(aux_en, wpm)
+        #uwerr.(aux_en)
+        #isnothing(wpm) ? uwerr.(aux_en) : [uwerr(a, wpm) for a in aux_en]
         eff_en[i] = copy(aux_en)
     end
     return eff_en
@@ -173,7 +174,7 @@ corr_upper = [corr_pa[1]]
 matrices = [corr_diag, corr_upper]
 
 ## solve the GEVP
-evals = getall_eigvals(matrices, 5) #where t_0=5
+#evals = getall_eigvals(matrices, 5) #where t_0=5
 
 
 Julia>
@@ -213,23 +214,46 @@ res = uweigvals(a)  ##eigenvalues
 res1 = uweigvals(a,b) ## generalised eigenvalues
 ```
 """
-function uweigvals(a::Matrix{uwreal}; iter = 30)
+function uweigvals(a::Matrix{uwreal}; iter::Int64=5, diag::Bool=true)
     n = size(a,1)
-    for k in 1:iter
-        q_k, r_k = qr(a)
-        a = uwdot(r_k, q_k)
+    if n == 2 
+        return uweigvals_dim_geq3(a, iter=iter, diag=diag) #uweigvals_dim2(a) 
+    else
+        return uweigvals_dim_geq3(a, iter=iter, diag=diag)
     end
-    return get_diag(a)
 end
-function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)
+function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 5, diag::Bool=true)
     c = uwdot(invert(b),a)
+    n = size(c,1)
+    if n == 2
+        return uweigvals_dim_geq3(c, iter=iter, diag=diag) #uweigvals_dim2(c) 
+    else
+        return uweigvals_dim_geq3(c, iter=iter, diag=diag)
+    end
+    #for k in 1:iter
+    #    q_k, r_k = qr(c)
+    #    c = uwdot(r_k, q_k)
+    #end
+end
+function uweigvals_dim_geq3(a::Matrix{uwreal}; iter = 5, diag::Bool=true)
+    n = size(a,1)
+    a = tridiag_reduction(a)
     for k in 1:iter
-        q_k, r_k = qr(c)
-        c = uwdot(r_k, q_k)
+        q_k, r_k = qr(a)
+        a = uwdot(r_k, q_k)
     end
-    return get_diag(c)
+    diag == true ? (return get_diag(a)) : (return a)
+end
+function uweigvals_dim2(a::Matrix{uwreal})
+    res = Vector{uwreal}(undef,2)
+    aux_sum = a[1,1] + a[2,2]
+    delta = ADerrors.sqrt((aux_sum)^2 - 4*(a[1,1]*a[2,2]- a[1,2]*a[2,1]))
+    res[1] = (aux_sum + delta) / 2
+    res[2] = (aux_sum - delta) / 2
+    return res
 end
 
+
 @doc raw"""
     getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30 )
 
@@ -257,7 +281,7 @@ mat_array = get_matrix(diag, upper_diag)
 evecs = getall_eigvecs(mat_array, 5)
 ```
 """
-function getall_eigvecs(a::Vector{Matrix{uwreal}}, delta_t; iter=30)
+function getall_eigvecs(a::Vector{Matrix}, delta_t; iter=5)
     n = length(a)-delta_t
     res = Vector{Matrix{uwreal}}(undef, n)
     [res[i] = uweigvecs(a[i+delta_t], a[i]) for i=1:n]
@@ -291,18 +315,22 @@ res = uweigvecs(a)  ##eigenvectors in column
 res1 = uweigvecs(a,b) ## generalised eigenvectors in column 
 ```
 """
-function uweigvecs(a::Matrix{uwreal}; iter = 30)
+function uweigvecs(a::Matrix{uwreal}; iter = 5)
     n = size(a,1)
+    if n > 2
+        a, q = tridiag_reduction(a, q_mat=true)
+    end
     u = idty(n)
     for k in 1:iter
         q_k, r_k = qr(a)
         a = uwdot(r_k, q_k)
         u = uwdot(u, q_k)
     end
-    return  u
+    n > 2 ? (return  uwdot(q,u)) : (return u)
 end
-function uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30)
+function uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 5)
     c  = uwdot(invert(b), a)
+    #return uweigvecs(c, iter=iter)
     n = size(c,1)
     u = idty(n)
     for k in 1:iter
@@ -342,12 +370,12 @@ b = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
 eval, evec = uweigen(a)   
 eval1, evec1 = uweigvecs(a,b) 
 ```
-"""
-function uweigen(a::Matrix{uwreal}; iter = 30)
-    return uweigvals(a, iter=iter), uweigvecs(a, iter=iter)
+""" 
+function uweigen(a::Matrix{uwreal}; iter = 5, diag::Bool=true)
+    return uweigvals(a, iter=iter,diag=diag), uweigvecs(a, iter=iter)
 end
-function uweigen(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30)
-    return uweigvals(a, b, iter=iter), uweigvecs(a, b, iter=iter)
+function uweigen(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 5, diag::Bool=true)
+    return uweigvals(a, b, iter=iter, diag=diag), uweigvecs(a, b, iter=iter)
 end
 
 function norm_evec(evec::Matrix{uwreal}, ctnot::Matrix{uwreal})
@@ -379,13 +407,32 @@ Given an input matrix of uwreals, this function returns the inverse of that matr
 backward substitution.
 """
 function invert(a::Matrix{uwreal})
-    #q,r = qr(hess_red(a))
     q,r = qr(a)
     r_inv = upper_inversion(r)
     return uwdot(r_inv, transpose(q))
 end
 
 
+"""
+Performs a Cholesky decomposition of A, which must 
+    be a symmetric and positive definite matrix. The function
+    returns the lower variant triangular matrix, L.
+"""
+function uwcholesky(A::Matrix{uwreal})
+    n = size(A,1)
+    L = fill(uwreal(0), n, n)
+    for i in 1:n 
+        for k in 1:i
+            tmp_sum = sum(L[i,j] * L[k,j] for j in 1:k)
+            if i==k 
+                L[i,k] = sqrt(((A[i,i] - tmp_sum)))
+            else
+                L[i,k] = (1.0 / L[k,k] * (A[i,k]-tmp_sum))
+            end
+        end
+    end
+    return L 
+end
 """
 qr(A::Matrix{uwreal})
 
@@ -441,13 +488,6 @@ function make_householder(a::Vector{uwreal})
     v[1] = uwreal(1.0) 
     H = idty(n)
     H = H -  (2.0 / uwdot(v,v)) * uwdot(reshape(v, :,1), transpose(v)) #this last term is a vector product  
-    for i in 1:n
-        for j in 1:n
-            if H[i,j]!= 0 &&  H[i,j]!=1 
-
-            end
-        end
-    end
     return H
 end
 
@@ -462,28 +502,99 @@ function make_householder(a::Vector{Float64})
     H = H .-  (2.0 / dot(v,v)) * mul!(y, v, transpose(v))
     return H
 end
-
+"""
+This function computes a vector u that generates a Householder reflection H = I - uu* satisfying Ha=nu*e1.
+Accumulating this transformations any matrix can
+be reduced to upper Hessenberg form.
+"""
+function house_gen(a::Vector{uwreal})
+    u = deepcopy(a)
+    nu = uwnorm(a)
+    if nu == 0 
+        u[1] = sqrt(2)
+    end
+    if u[1] != 0 
+        rho = u[1] / sqrt(u[1]*u[1])
+    else
+        rho = 1
+    end
+    u = uwdot(rho/nu, u) 
+    u[1] += 1
+    u =  uwdot(u , 1/sqrt(u[1]))
+    nu = -(rho) * nu
+    return u, nu 
+end
 """
 hess_red(a::Matrix{uwreal})
 
-Given a matrix of uwreal variables, this function returns the Hessenberg form of the matrix.
+Given a matrix of uwreal variables, this function returns the Hessenberg form of the matrix by applying householder transformations.
 """
-function hess_red(a::Matrix{uwreal})
-    n = size(a, 1)
-    res = Matrix{uwreal}(undef, n, n)
-    res = a
-    for j in 1:n-2
-        mi = idty(n)
-        mi_inv = idty(n)
-        for k in j+2:n
-            mi[k,j+1] = res[k,j] / res[j+1, j]
-            mi_inv[k,j+1] = -mi[k,j+1]
+function hess_reduce(a::Matrix{uwreal}; q_mat::Bool=false)
+    n = size(a,1)
+    H = deepcopy(a)
+    Q = idty(n)
+    
+    for k in 1:n-2
+        u, H[k+1,k] = house_gen(H[k+1:n,k])
+        Q[k+1:n,k] = u
+
+        v = uwdot(u,H[k+1:n, k+1:n])
+        H[k+1:n, k+1:n] .= H[k+1:n, k+1:n] .- uwdot(reshape(u,length(u),1),v)
+        H[k+2:n,k] .= 0.0 
+        
+        v = uwdot(H[1:n,k+1:n], u)
+        H[1:n,k+1:n] .= H[1:n,k+1:n] .- uwdot(v,reshape(u,1,length(u)))
+    end
+    
+    if q_mat 
+        Id = idty(n)
+        for k  in  reverse(1:n-2)
+            u = Q[k+1:n, k]
+            v = uwdot(u, Q[k+1:n,k+1:n])
+            Q[k+1:n,k+1:n] .= Q[k+1:n,k+1:n] .- uwdot(reshape(u,length(u),1),v)
+            Q[:,k] = Id[:,k]
         end
-        res = uwdot(mi_inv, uwdot(res, mi) )
+        return H,Q
+    else
+        return H
     end
-    return res
+    
 end
+"""
 
+
+"""
+function tridiag_reduction(a::Matrix{uwreal}; q_mat::Bool=false)
+    n = size(a,1)
+    H = deepcopy(a)
+    Q = idty(n)
+    for k in 1:n-2
+        u, H[k+1,k] = house_gen(H[k+1:n,k])
+        Q[k+1:n,k] = u
+
+        v = uwdot(u,H[k+1:n, k+1:n])
+        H[k+1:n, k+1:n] .= H[k+1:n, k+1:n] .- uwdot(reshape(u,length(u),1),v)
+        H[k+2:n,k] .= 0.0 
+        
+        v = uwdot(H[1:n,k+1:n], u)
+        H[1:n,k+1:n] .= H[1:n,k+1:n] .- uwdot(v,reshape(u,1,length(u)))
+        H[k,k+2:n] .= 0.0 
+    end
+    
+    if q_mat 
+        Id = idty(n)
+        for k  in  reverse(1:n-2)
+            u = Q[k+1:n, k]
+            v = uwdot(u, Q[k+1:n,k+1:n])
+            Q[k+1:n,k+1:n] .= Q[k+1:n,k+1:n] .- uwdot(reshape(u,length(u),1),v)
+            Q[:,k] = Id[:,k]
+        end
+        return H,Q
+    else
+        return H
+    end
+    
+end
 """
 upper_inversion(u::Matrix{uwreal})
 
@@ -556,7 +667,7 @@ function uwdot(a::Matrix{uwreal}, b::Matrix{uwreal})
     k,p = size(b)
     c = Matrix{uwreal}(undef, n, p)
     if m != k
-        error("Error uwdot")
+        error("Error uwdot: in a*b the size of a is ", size(a), " while the size of b is ", size(b) )
     end
     for i = 1:n
         for j = 1:p
@@ -573,7 +684,7 @@ It returns the norm of a vector of uwreal
 """
 function uwnorm(a::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
     norm = sqrt(uwdot(a,a))
-    isnothing(wpm) ? uwerr(norm) : uwerr(norm, wpm)
+    #isnothing(wpm) ? uwerr(norm) : uwerr(norm, wpm)
     return norm 
 end
 
@@ -630,8 +741,34 @@ function idty(n::Int64)
     return res
 end
 
+function make_positive_def(A::Matrix)
+    B = 0.5 * (A + A')
+    U, S, V = svd(B)
+    H = V * Diagonal(S) * V'
+    A_hat = 0.5 * (B + H)
+    A_hat = 0.5 * ( A_hat + A_hat')
+    k = 0
+    while !isposdef(A_hat)
+        mineig = minimum(eigvals(A_hat))
+        eps =  2.220446049250313e-16
+        A_hat = A_hat + (-mineig*k^2 .+ eps*Matrix(I, size(A)))
+        k+=1
+    end
+    return A_hat
+end
+
+function invert_covar_matrix(A::Matrix)
+    F = svd(A)
+    cov_inv = F.V * Diagonal(1 ./F.S) * F.U'
+    return cov_inv  
+end
+
+function more_penrose_inv(A::Matrix)
+    F = svd(A)
+
+end
+
 
-# operator overloading
 
 """
 Base.:*(x::uwreal, y::Matrix{uwreal})

src/juobs_tools.jl

@@ -300,7 +300,7 @@ function md_val(a::uwreal, obs::Corr, derm::Vector{Corr})
 end
 
 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, wpm)
+    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 
@@ -329,7 +329,7 @@ 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.(y, wpm)
+    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)
@@ -382,24 +382,51 @@ fit_routine(model, xdata, ydata, param=3)
 fit_routine(model, xdata, ydata, param=3, covar=true)
 ```
 """
-function 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)
-    #isnothing(wpm) ? uwerr.(ydata) : uwerr.(ydata, wpm)
+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, correlated_fit::Bool=false)
+    isnothing(wpm) ? uwerr.(ydata) : [uwerr(yaux, wpm) for yaux in ydata]
     
     yval = value.(ydata)
     yer = err.(ydata)
     
-    # Generate chi2 + solver
-    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)
+    if !correlated_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)
+    else
+        covar_inv = make_positive_def(invert_covar_matrix(cov(ydata))) 
+        #covar_inv = (covar_inv + covar_inv')/2
+        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 )
+    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
-    for i = 1:length(upar)
-        isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm)
-        print("\n Fit parameter: ", i, ": ")
-        details(upar[i])
+    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
     end
-    println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", length(yval) - param,")")
-    return upar
+    #chis2_corrected = (length(yval) - param) * chisq(coef(fit), ydata) / chi_exp
+    chis2_corrected = (length(yval) - param) * chi2_fit_res / chi_exp
+
+    #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)
+    return upar, chi2_fit_res / chi_exp
 end
 
 function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; 
@@ -411,8 +438,8 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
         uwerr.(ydata)
         uwerr.(xdata)
     else
-        uwerr.(ydata, wpm)
-        uwerr.(xdata, wpm)
+        [uwerr(yaux, wpm) for yaux in ydata] 
+        [uwerr(xaux, wpm) for xaux in xdata]
     end
 
     yval = value.(ydata)
@@ -454,6 +481,10 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
         
         (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)
+        return upar, min_fun_cov(sol.minimizer) / chi2_exp
+
     else
         chisq_full(p, d) = get_chi2(model, d, ddat, p, Nalpha)
         min_fun(t) = chisq_full(t, dat)
@@ -461,26 +492,98 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
         
         (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)
+        return upar, min_fun(sol.minimizer) / chi2_exp
     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
+    #for i = 1:length(upar)
+    #    isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm)
+    #    print("\n Fit parameter: ", i, ": ")
+    #    details(upar[i])
+    #end
     
 end
 
+function global_fit_routine(models::Vector{Function}, xdata::Vector{Array{Float64, 1}}, ydata::Vector{Array{uwreal, 1}}, n_par::Int64; guess_param::Vector{Float64}=fill(0.5, n_par), log::Bool=true, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    if !(length(models) == length(ydata))
+        error("Dimension mismatch")
+    end
+    N_sets = length(models)             # number of datasets
+    x_flat = xdata[1]                   # unrolling xdata 
+    y_flat = ydata[1]                   # unrolling ydata
+
+    chunk_idx = Vector{Vector{Int64}}(undef, N_sets)
+    j = 1
+    for i = 1:N_sets
+        isnothing(wpm) ? uwerr.(ydata[i]) :  [uwerr(yaux, wpm) for yaux in ydata[i]]
+        if i > 1 
+            y_flat = vcat(y_flat, ydata[i])
+            x_flat =vcat(x_flat, xdata[i])
+        end
+        stp = j + length(ydata[i]) - 1
+        chunk_idx[i] = collect(j:stp)
+        j = stp + 1
+    end
+    function flat_model(x_flat::Vector{Float64}, param::Vector{Float64})
+        fx = similar(x_flat)
+        for ii in 1:N_sets
+            mod = models[ii]
+            fx[chunk_idx[ii]] = mod(x_flat[chunk_idx[ii]], param)
+        end
+        return fx
+    end
+    function flat_chisq(param, data)
+        chi = 0.0
+        for ii in 1:N_sets
+            mod = models[ii]
+            # uncorrelated chi2
+            chi += sum((data[chunk_idx[ii]] .- mod(x_flat[chunk_idx[ii]], param)).^2 ./ err.(y_flat[chunk_idx[ii]]).^2 )        
+            # correlated chi2
+        end
+        return chi
+    end
+    fit = curve_fit(flat_model, x_flat, value.(y_flat), 1.0 ./ err.(y_flat).^2, guess_param )
+    #uncorrelated chi2exp
+    (upar, chi2_exp) = isnothing(wpm) ? fit_error(flat_chisq, coef(fit) , y_flat) : fit_error(flat_chisq, coef(fit) , y_flat, wpm)
+    #correlated chi2exp
+    #(upar, chi2_exp) = isnothing(wpm) ? fit_error(flat_chisq, coef(fit) , y_flat, W=cov(y_flat)^(-1)) : fit_error(flat_chisq, coef(fit) , y_flat, wpm)
+    chisq = value(flat_chisq(coef(fit), y_flat))  
+    if log
+        println("Converged param: ", fit.converged)
+        println("Chisq / chiexp: ", chisq, " / ", chi2_exp, " (dof: ", dof(fit),")")
+        chis2_corrected = (dof(fit)) * chisq / chi2_exp
+        println("Chisq corrected: ", chis2_corrected)
+    end
+    return upar, chisq, chi2_exp, dof(fit)
+end
+
 function gen_chisq(f::Function, x::Array{<:Real}, err::Vector{Float64}) #constrained
     chisq(par, dat) = sum((dat .- f(x,par)).^2 ./err.^2)
     return chisq
 end
+function gen_chisq_correlated(f::Function, x::Array{<:Real}, covar_inv::Matrix, par, dat )
+    chi2 = 0.0
+    #delta = dat - f(x, par)
+    #chi2 = (reshape(delta, 1,:) * covar_inv * delta )[1]
+    for i in 1:length(dat)
+        for j in 1:length(dat)
+            #println(f(x,par)[i])
+            #println(f(x,par)[j])
+            #println(f(x[j,:],par))
+            #println(f(x[i,:],par))
+            chi2 = chi2 +  (dat[i] - f(x,par)[i]) * covar_inv[i,j] * (dat[j] - f(x,par)[j])
+        end
+    end
+    #chisq(par, dat) = sum(
+    #    transpose((dat .- f(x,par))) * covar_inv * (dat .- f(x,par))
+    #) 
+    return chi2
+end
 
 function get_chi2(f::Function, data, ddata, par, Nalpha) #full
     chi2 = 0.0

src/juobs_types.jl

@@ -145,13 +145,17 @@ mutable struct Corr
     mu::Vector{Float64}
     gamma::Vector{String}
     y0::Int64
+    theta1::Vector{Float64}
+    theta2::Vector{Float64}
     function Corr(a::Vector{uwreal}, b::CData)
         h = getfield(b, :header)
         kappa = [h.k1, h.k2]
         mu = [h.mu1, h.mu2]
         gamma = [gamma_name[h.type1+1], gamma_name[h.type2+1]]
         y0 = Int64(h.x0)
-        return new(a, kappa, mu, gamma, y0)
+        theta1 = h.theta1
+        theta2 = h.theta2
+        return new(a, kappa, mu, gamma, y0, theta1, theta2)
     end 
     function Corr(a::Vector{uwreal}, b::Vector{CData})
         sym = [:k1, :k2, :mu1, :mu2, :type1, :type2, :x0]
@@ -165,7 +169,9 @@ mutable struct Corr
         mu = [h[1].mu1, h[1].mu2]
         gamma = [gamma_name[h[1].type1+1], gamma_name[h[1].type2+1]]
         y0 = Int64(h[1].x0)
-        return new(a, kappa, mu, gamma, y0)
+        theta1 = h[1].theta1
+        theta2 = h[1].theta2
+        return new(a, kappa, mu, gamma, y0, theta1, theta2)
     end 
     Corr(o::Vector{uwreal}, k::Vector{Float64}, m::Vector{Float64}, g::Vector{String}, y::Int64) = new(o, k, m, g, y)
 end