Merge branch 'merging' into alejandro

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, 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]
+const 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)
+const M_error = [0.05, 0.05, 0.07, 0.4, 1.2, 0.011, 0.12, 0.6]
 #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]
 const ZM_error = [35, 27, 33, 38, 45] .* 1e-4
 const ZM_tm_data = [2.6047, 2.6181, 2.6312, 2.6339, 2.6127]
 const ZM_tm_error = [42, 33, 42, 48, 55] .* 1e-4
-# 1808.09236
+#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]
@@ -17,41 +17,63 @@ 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]
-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
+#2101.02694
+# const b_values2 = [3.40, 3.46, 3.55, 3.70]
+const t0_data = [2.846, 3.634, 5.140, 8.564, 14.064]
+const t0_error = [8, 13, 21, 24, 63] .* 1e-3
+# from AS analysis 
+const t0_ph_value = [0.4118]
+const t0_ph_error = ones(1,1) .* 25e-4
+#2101.10969
+const RM_data = [2.335, 1.869, 1.523, 1.267, 1.149]
+const RM_err  = [31, 19, 14, 16, 18] .* 1e-3
+#1906.03445 LCP1 for heavy quarks 
+const BA_MINUS_BP_data = [-0.324, -0.265, -0.196, -0.119, -0.073]
+const BA_MINUS_BP_err  = [17, 14, 14, 14, 12] .*1e-3 
+const ZDATA = [0.8981, 0.9468, 1.0015, 1.0612, 1.0971]
+const ZERR  = [35, 35, 30, 17, 18 ] .* 1e-4
 
+#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 = 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)
+# same for zm wilson
+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]]
+zz = cobs([2.270073, 0.121644, -0.464575], CC, "zm")
+z_covar(b) = value(Mrat) * (zz[1] + zz[2]*(b-3.79) + zz[3]*(b-3.79)^2)
+
+## taking into account correlations in r_m
+#Crm = [[3.699760, -2.198586, -1.476913e-3] [-2.198586, 1.306512, 8.776569e-4] [-1.476913e-3, 8.776569e-4, 5.895710e-7] ]
+#c_i = cobs([-7.86078, 5.49175, -0.54078], Crm, "c_i in rm")
+#rm(b::Float64) = 1.0 +  0.004630*(6/b)^4 * ((1. +c_i[1]*(6/b)^2 + c_i[2]*(6/b)^4) / (1. + c_i[3]*(6/b)^2))
+
 ##
 #Cov matrices
 const C1 = zeros(5, 5)
 const C2 = zeros(5, 5)
-const C3 = zeros(4, 4)
-const C4 = zeros(7, 7)
+const C3 = zeros(5, 5)
+const C4 = zeros(8, 8)
 const C5 = zeros(5, 5)
 const C6 = zeros(5, 5)
 const C7 = zeros(5, 5)
+const C8 = zeros(5, 5)
+const C9 = zeros(5, 5)
+const C10 = zeros(5, 5)
 
-for i = 1:7
+for i = 1:8
     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
         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
+        C7[i, i] = ZV_err[i] ^ 2       
+        C8[i, i] = RM_err[i] ^2 
+        C9[i,i]  = BA_MINUS_BP_err[i] ^2
+        C10[i,i]  = ZERR[i] ^2
+        C3[i, i] = t0_error[i] ^ 2
     end
 
 end
@@ -63,12 +85,18 @@ 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")
+const ZS_over_ZP = cobs(ZS_over_ZP_data, C6, "ZS/ZP")
+const RM = cobs(RM_data, C8, "rm")
+const BA_BP = cobs(BA_MINUS_BP_data, C9, "ba-bp")
+const ZZ = cobs(ZDATA, C10, "Z")
 
-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]
-zv(beta::Float64) = Zv[b_values .== beta][1]
-zs_over_zp(beta::Float64) = ZS_over_ZP[b_values .== beta][1]
+zm(beta::Float64)         = ZM[b_values .== beta][1]
+zm_tm(beta::Float64)      = ZM_tm[b_values .== beta][1]
+t0(beta::Float64)         = t0_[b_values .== beta][1]
+a(beta::Float64)          = a_[b_values .== beta][1]
+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
+rm(beta::Float64)         = RM[b_values .== beta][1]
+ba_bp(beta::Float64)      = BA_BP[b_values .== beta][1]
+Z(beta::Float64)          = ZZ[b_values .== beta][1]
\ No newline at end of file

src/juobs.jl

@@ -12,6 +12,6 @@ include("juobs_obs.jl")
 export read_mesons, read_mesons_correction, read_ms1, read_ms, read_md, truncate_data!, concat_data!
 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_obs_TSM, corr_sym, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine, global_fit_routine, bayesian_av, pvalue
-export meff, mpcac, dec_const, dec_const_pcvc, comp_t0
+export meff, mpcac, dec_const, dec_const_w, dec_const_pcvc, comp_t0
 
 end # module

src/juobs_linalg.jl

@@ -62,7 +62,7 @@ function get_matrix(corr_diag::Vector{Vector{uwreal}}, corr_upper::Vector{Vector
     if d != (n*(n-1) / 2)
         error("Error get_matrix")
     end
-    res = Vector{Matrix}(undef, time) # array with all the matrices at each time step.
+    res = Vector{Matrix{uwreal}}(undef, time) # array with all the matrices at each time step.
     aux = Matrix{uwreal}(undef, n, n)
     for t in 1:time  
         count=0
@@ -180,10 +180,10 @@ matrices = [corr_diag, corr_upper]
 Julia>
 ```
 """
-function getall_eigvals(a::Vector{Matrix{uwreal}}, t0::Int64; iter=30 )
+function getall_eigvals(a::Vector{Matrix{uwreal}}, t0::Int64; iter=5 )
     n = length(a)
     res = Vector{Vector{uwreal}}(undef, n)
-    [res[i] = uweigvals(a[i], a[t0]) for i=1:n]
+    [res[i] = uweigvals(a[i], a[t0], iter=iter) for i=1:n]
     return res
 end
 
@@ -237,7 +237,7 @@ function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 5, diag::Bool=t
 end
 function uweigvals_dim_geq3(a::Matrix{uwreal}; iter = 5, diag::Bool=true)
     n = size(a,1)
-    a = tridiag_reduction(a)
+    #a = tridiag_reduction(a)
     for k in 1:iter
         q_k, r_k = qr(a)
         a = uwdot(r_k, q_k)
@@ -281,7 +281,7 @@ mat_array = get_matrix(diag, upper_diag)
 evecs = getall_eigvecs(mat_array, 5)
 ```
 """
-function getall_eigvecs(a::Vector{Matrix}, delta_t; iter=5)
+function getall_eigvecs(a::Vector{Matrix{uwreal}}, 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]
@@ -315,20 +315,21 @@ res = uweigvecs(a)  ##eigenvectors in column
 res1 = uweigvecs(a,b) ## generalised eigenvectors in column 
 ```
 """
-function uweigvecs(a::Matrix{uwreal}; iter = 5)
+function uweigvecs(a::Matrix{uwreal}; iter = 10)
     n = size(a,1)
-    if n > 2
-        a, q = tridiag_reduction(a, q_mat=true)
-    end
+    # 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
-    n > 2 ? (return  uwdot(q,u)) : (return u)
+    #n > 2 ? (return  uwdot(q,u)) : (return u)
+    return u
 end
-function uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 5)
+function uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 10)
     c  = uwdot(invert(b), a)
     #return uweigvecs(c, iter=iter)
     n = size(c,1)
@@ -691,9 +692,15 @@ end
 """
 Return the transpose of a matrix or vector of uwreals
 """
-function transpose(a::Matrix{uwreal})
+
+function Base.transpose(a::Matrix{uwreal})
     res = similar(a)
     n = size(a, 1)
+    if size(a,1) == 1
+        return reshape(a, :,1)
+    elseif size(a,2) == 1
+        return reshape(a, 1,:)    
+    end
     for i in 1:n-1
         for j in i+1:n
             temp = a[i,j]
@@ -705,7 +712,7 @@ function transpose(a::Matrix{uwreal})
     res[n,n] = a[n,n]
     return res
 end
-function transpose(vec::Vector{uwreal}) 
+function Base.transpose(vec::Vector{uwreal}) 
     res = deepcopy(vec)
     return reshape(res,1,:)
 end
@@ -741,25 +748,37 @@ 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
+function make_positive_def(A::Matrix; eps::Float64=10^(-14))
+    vals, vecs = eigen(Symmetric(A))
+    idx = findall(x-> x<0.0, vals)
+    vals[idx] .= eps
+    return  Symmetric(vecs * Diagonal(vals) * vecs')
+    #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)
+function invert_covar_matrix(A::Matrix; eps::Float64=10^(-9))
     F = svd(A)
-    cov_inv = F.V * Diagonal(1 ./F.S) * F.U'
+    s_diag = F.S
+    for i in 1:length(s_diag)
+        if s_diag[i] < eps
+            s_diag[i] = 0.0
+        else
+            s_diag[i] = 1 / s_diag[i]
+        end
+    end
+    cov_inv = F.V * Diagonal(s_diag) * F.U'
     return cov_inv  
 end
 
@@ -769,7 +788,80 @@ function more_penrose_inv(A::Matrix)
 end
 
 
+########################
+# OPERATOR OVERLOADING
+########################
 
+function Base.:*(x::uwreal, y::Array{Any})
+    N = size(y, 1)
+    return fill(x, N) .* y
+end  
+Base.:*(x::Array{Any}, y::uwreal) = Base.:*(y,x)
+function Base.:*(x::uwreal, y::Array{Float64})
+    N = size(y, 1)
+    return fill(x, N) .* y
+end
+Base.:*(x::Array{Float64}, y::uwreal) = Base.:*(y,x)
+function Base.:*(x::uwreal, y::Array{uwreal})
+    N = size(y, 1)
+    return fill(x, N) .* y
+end
+Base.:*(x::Array{uwreal}, y::uwreal) = Base.:*(y,x)
+function Base.:+(x::uwreal, y::Vector{uwreal})
+    N = size(y, 1)
+    return fill(x, N) .+ y
+end
+function Base.:+(x::Float64, y::Vector{Float64})
+    N = size(y, 1)
+    return fill(x, N) .+ y
+end
+function Base.:+(x::Float64, y::Vector{uwreal})
+    N = size(y, 1)
+    return fill(x, N) .+ y
+end
+function Base.:-(x::Float64, y::Vector{Any})
+    N = size(y, 1)
+    return fill(x, N) .- y
+end
+function Base.:-(x::Float64, y::Vector{Float64})
+    N = size(y, 1)
+    return fill(x, N) .- y
+end
+function Base.:-(x::Float64, y::Vector{uwreal})
+    N = size(y, 1)
+    return fill(x, N) .- y
+end
+
+function Base.length(uw::uwreal)
+    return length(value(uw))
+end
+function Base.iterate(uw::uwreal)
+    return iterate(uw, 1)
+end
+function Base.iterate(uw::uwreal, state::Int64)
+    if state > length(uw)
+        return nothing
+    else
+        return uw[state], state +1
+    end
+end
+function Base.iterate(uw::Vector{uwreal})
+    return iterate(uw, 1)
+end
+function Base.iterate(uw::Vector{uwreal}, state::Int64)
+    if state > length(uw)
+        return nothing
+    else
+        return uw[state], state +1
+    end
+end
+function Base.getindex(uw::uwreal, ii::Int64)
+    idx = getindex(value(uw), ii)
+    return  uw # uwreal([getindex(value(uw), kwargs...), err(uw)], " ")
+end
+function Base.abs2(uw::uwreal)
+    return (uw^2)^0.5
+end
 """
 Base.:*(x::uwreal, y::Matrix{uwreal})
 

src/juobs_reader.jl

@@ -330,21 +330,130 @@ function read_rw(path::String; v::String="1.2")
     return r_data
 end
 
+@doc raw"""
+read_rw_openQCD2(path::String; print_info::Bool=false)
+
+This function reads the reweighting factors generated with openQCD version 2.#.
+The flag print_info if set to true print additional information for debugging
+"""
+function read_rw_openQCD2(path::String; print_info::Bool=false)
+    
+    data = open(path, "r")
+    nrw = read(data, Int32) 
+    nrw = Int32(nrw / 2)
+    
+    nfct = Array{Int32}(undef, nrw)
+    read!(data, nfct)       
+    
+    nsrc = Array{Int32}(undef, nrw)
+    read!(data, nsrc)        
+    null = read(data, Int32)  
+    if null !== Int32(0)   
+        error("In openQCD 2.0 this Int32 should be a zero.")
+    end
+
+    data_array = Array{Array{Float64}}(undef, nrw)
+    [data_array[k] = Array{Float64}(undef, 0) for k in 1:nrw]
+    vcfg = Vector{Int32}(undef, 0) 
+    while !eof(data)
+
+        push!(vcfg, read(data, Int32)) 
+        if print_info
+            println("\n ######## cnfg: ", vcfg[end])
+        end
+
+        for k in 1:nrw
+            read_array_rwf_dat_openQCD2(data)
+            tmp_rw, n = read_array_rwf_dat_openQCD2(data)
+
+            tmp_nfct=1.0
+            for j in 1:n[1]
+                tmp_nfct *= mean((exp.(.-tmp_rw[j])))
+            end
+            push!(data_array[k], tmp_nfct)
+        end
+    end
+
+    return permutedims(hcat(data_array...), (2,1))
+end
+
+
+function read_array_rwf_dat_openQCD2(data::IOStream; print_info::Bool=false)
+    
+    d = read(data, Int32) 
+    n = Array{Int32}(undef, d)
+    read!(data, n)   
+    size = read(data, Int32) 
+    m = prod(n)
+    
+    if print_info
+        println("d: ", d)
+        println("n: ", n)
+        println("size: ", size)
+        println("m: ", m)
+    end
+
+    if size == 4
+        types = Int32
+    elseif size == 8
+        types = Float64
+    elseif size == 16
+        types = Float64x2
+    else
+        error("No type with size=$(size) supported")
+    end
+
+    tmp_data = Array{types}(undef, m)
+    read!(data, tmp_data)
+
+    res = parse_array_openQCD2(d, n, tmp_data, quadprec=true) 
+
+    return res, n
+end
+
+function parse_array_openQCD2(d, n, dat; quadprec=true)
+    
+    if d != 2
+        error("dat must be a two-dimensional array")
+    end
+    res = Vector{Vector{Float64}}(undef, 0)
+    
+    for k in range(1,n[1])
+        tmp = dat[(k-1)*n[2]+1:k*n[2]]
+        if quadprec
+            tmp2 = Vector{Float64}(undef, 0)
+            for j in range(start=1,step=2,stop=length(tmp))
+                push!(tmp2, tmp[j])
+            end
+            push!(res, tmp2)
+        else
+            push!(res, tmp)
+        end
+    end
+
+    return res
+end
+
 @doc raw"""
     read_ms1(path::String; v::String="1.2")
 
 Reads openQCD ms1 dat files at a given path. This method returns a matrix `W[irw, icfg]` that contains the reweighting factors, where
 `irw` is the `rwf` index and icfg the configuration number.
-The function is compatible with the output files of openQCD v=1.2, 1.4 and 1.6. Version can be specified as argument.
+The function is compatible with the output files of openQCD v=1.2, 1.4, 1.6 and 2.0. Version can be specified as argument.
 
 Examples:
 ```@example    
 read_ms1(path)
 read_ms1(path, v="1.4")
 read_ms1(path, v="1.6")
+read_ms1(path, v="2.0")
 ```
 """
 function read_ms1(path::String; v::String="1.2")
+
+    if v == "2.0"
+        return read_rw_openQCD2(path)
+    end
     r_data = read_rw(path, v=v)
     nrw = length(r_data)
     ncnfg = size(r_data[1])[3]
@@ -417,6 +526,8 @@ function read_ms(path::String; id::Union{String, Nothing}=nothing, dtr::Int64=1
     ntr = Int32((fsize - 3*4 - 8) / datsize)
 
     if mod(ntr, dtr) != 0
+        println("ntr = ", ntr)
+        println("dtr = ", dtr)
         error("ntr / dtr must be exact")
     end
 
@@ -573,7 +684,11 @@ function concat_data!(data1::Vector{Vector{CData}}, data2::Vector{Vector{CData}}
             data1[k][r].vcfg = vcat(data1[k][r].vcfg, data2[k][r].vcfg)
             data1[k][r].re_data = vcat(data1[k][r].re_data, data2[k][r].re_data)
             data1[k][r].im_data = vcat(data1[k][r].im_data, data2[k][r].im_data)
+            idx = sortperm(data1[k][r].vcfg)
+            data1[k][r].vcfg = data1[k][r].vcfg[idx]
+            data1[k][r].re_data = data1[k][r].re_data[idx, :] 
+            data1[k][r].im_data = data1[k][r].im_data[idx, :] 
         end
     end
     return nothing
-end
+end
\ No newline at end of file

src/juobs_types.jl

@@ -182,7 +182,7 @@ mutable struct Corr
         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)
+    Corr(o::Vector{uwreal}, k::Vector{Float64}, m::Vector{Float64}, g::Vector{String}, y::Int64, th1::Vector{Float64}, th2::Vector{Float64} ) = new(o, k, m, g, y, th1, th2)
 end
 
 mutable struct YData