md_val for new version

Julien's most updated code has all valence derivatives with the sign changed, so md_val needs to flip the sign. Also, the reader does not write "_d1, _d2" in the derivatives anymore for the new data files, so md_val does not look for "_d1, _d2". You need to select by hand which correlators are your valence derivatives, looking at the log file e.g.

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.0"
 [deps]
 ADerrors = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
 BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"

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,93 @@ 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_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]
+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]
+
+const ens_db = Dict(
+    #"ens_id"=>[L, T, beta, dtr, vrw, trunc_cnfg]
+    "H101"     => [32, 96, 3.4, 2, "1.2", [1001,1009]],
+    "H102r001" => [32, 96, 3.4, 2, "2.0", [997]],
+    "H102r002" => [32, 96, 3.4, 2, "2.0", [1008]],
+    "H105"     => [32, 96, 3.4, 2, "2.0", [947,1042]],
+    "H105r005" => [32, 96, 3.4, 1, "1.2", [837]],
+    "H400"     => [32, 96, 3.46, 1, "1.2", [505,540]],
+    "D450"     => [64, 128, 3.46, 1, ["2.0"], [1000]],
+    "N202"     => [48, 128, 3.55, 2, "1.2", [899]],
+    "N203"     => [48, 128, 3.55, 1, "2.0", [756,787]],
+    "N200"     => [48, 128, 3.55, 1, "2.0", [856,856]],
+    "D200"     => [64, 128, 3.55, 2, "2.0", [2001]],
+    "E250"     => [96, 192, 3.55, 1, ["2.0"], [1009]],
+    "N300"     => [48, 128, 3.70, 1, "1.2", [1521]],
+    "N302"     => [48, 128, 3.70, 1, "1.2", [2201]],
+    "J303"     => [64, 192, 3.70, 2, "2.0", [1073]],
+    "E300"     => [96, 192, 3.70, 1, ["1.4"], [1139]],
+    "J500"     => [64, 192, 3.85, 2, ["1.2", "1.4", "2.0"], [789,655,431]],
+    "J501"     => [64, 192, 3.85, 1, ["1.2", "1.4", "2.0"], [1635,1142,1150]]
+)
+
+
+const db = Dict(
+    "J501" => ["ts001_chunk1", "ts001_chunk3", "ts190_chunk1", "ts190_chunk2"],
+    "J500" => ["ts001_chunk1", "ts001_chunk2", "ts190_chunk1", "ts190_chunk2"],
+    "E250" => ["ts001", "ts097"],
+    "E300" => ["ts001_chunk2", "ts001_chunk3", "ts190_chunk1"],#, "ts001_chunk1"],
+    "D450" => ["ts001", "ts065_1", "ts065_2"]
+)
+
+const db_c = Dict(
+    "J501" => ["ts001", "ts190"],
+    "J500" => ["ts001_chunk1", "ts190_chunk2"],
+    "E250" => ["ts001"],
+    "E300" => ["ts001_chunk2"],#, "ts001_chunk1"]
+    "D450" => ["ts001_1", "ts001_2"]
+)
+
+const sym_bool = Dict(
+    "J501" => true,
+    "J500" => true,
+    "E250" => false,
+    "E300" => false,
+    "D450" => false
+)
+
+const flag_s = Dict(
+    "J303r003" => [324,325,326],
+    "H105r001" => [100,105,106],
+    "H105r002" => [1],
+    "H105r005" => [254, 255, 256, 257, 259, 260, 261, 264, 265, 266, 269, 280, 282, 283, 284, 285, 286, 287, 288, 289, 291, 299, 301, 313, 314, 315, 316, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342]
+)
 
-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]
+const ensemble_inv = Dict(
+    "H101"     => 1,
+    "H102r001" => 2,
+    "H102r002" => 3,
+    "H105"     => 4,
+    "H105r005" => 5,
+    "H400"     => 6,
\ No newline at end of file
+    "D450"     => 7,
+    "N202"     => 8,
+    "N203"     => 9,
+    "N200"     => 10,
+    "D200"     => 11,
+    "E250"     => 12,
+    "N300"     => 13,
+    "N302"     => 14,
+    "J303"     => 15,
+    "E300"     => 16,
+    "J500"     => 17,
+    "J501"     => 18
+)
\ No newline at end of file

src/juobs.jl

 module juobs
-using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, Optim
+using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, Optim, ForwardDiff
 using Base: Float64
 import Statistics: mean
 
@@ -8,10 +8,14 @@ include("juobs_linalg.jl")
 include("juobs_reader.jl")
 include("juobs_tools.jl")
 include("juobs_obs.jl")
+include("../constants/juobs_const.jl")
+include("../constants/path_csv.jl")
 
-export read_mesons, read_mesons_correction, read_ms1, read_ms, read_md, truncate_data!
+export EnsInfo
+export ens_db, db, db_c, flag_s, sym_bool, ensemble_inv
+export read_mesons, read_mesons_correction, read_mesons_multichunks, read_mesons_correction_multichunks, read_ens_tm, read_ens_tm_sym, read_ens_TSM, read_ens_wil, read_ens_csv, get_corr_tm, get_corr_TSM, get_corr_TSM_multichunks, get_corr_wil, read_ms1, read_ms, read_md, get_YM, get_YM_dYM, 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
-export meff, mpcac, dec_const, dec_const_pcvc, comp_t0
+export corr_obs, corr_obs_TSM, corr_sym, corr_sym_E250, corr_sym_D450, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine, global_fit_routine, bayesian_av, pvalue, fve, model_av, fit_alg
+export meff, mpcac, dec_const, dec_const_w, dec_const_pcvc, comp_t0, get_m, get_m_pbc, get_mpcac, get_f_wil, get_f_wil_pbc, get_f_tm, get_f_tm_pbc, get_w0t0_plat
 
 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})