apply_rw updated with gaps in measurements.

corr_obs_TMS updated with gaps in measurements and automatic detection of the gaps

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]
@@ -54,7 +54,7 @@ z_covar(b) = value(Mrat) * (zz[1] + zz[2]*(b-3.79) + zz[3]*(b-3.79)^2)
 const C1 = zeros(5, 5)
 const C2 = zeros(5, 5)
 const C3 = zeros(5, 5)
-const C4 = zeros(7, 7)
+const C4 = zeros(8, 8)
 const C5 = zeros(5, 5)
 const C6 = zeros(5, 5)
 const C7 = zeros(5, 5)
@@ -62,7 +62,7 @@ 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

src/juobs_linalg.jl

@@ -692,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]
@@ -706,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

src/juobs_obs.jl

@@ -672,21 +672,29 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
         e = err.(t2E)
 
         figure()
-        errorbar(x, v, e, fmt="x")
-        
-        errorbar(value(t0), 0.3, xerr=err(t0), fmt="x")
-        ylabel(L"$t^2E$")
+        errorbar(x, v, e, fmt="d", capsize=2, mfc="none", color="black")
+        errorbar(value(t0), 0.3, xerr=err(t0), capsize=2, color="red", fmt="d")
+        xarr = Float64.(range(5.0, 5.2, length=100))
+        yarr = model(xarr, par); uwerr.(yarr)
+        fill_between(xarr, value.(yarr) .- err.(yarr), value.(yarr) .+ err.(yarr), color="royalblue", alpha=0.2 )
+        ylabel(L"$t^2\langle E(x_0, t)\rangle$")
         xlabel(L"$t/a^2$")
+        tight_layout()
+        savefig("/Users/alessandroconigli/Desktop/t0_interp.pdf")
         display(gcf())
 
         figure()
         t_pl = dt0 + 1
-        plot(collect(1:xmax), value.(Y_aux[:, t_pl]) * t[nt0]^2 / L^3, "x") 
-        fill_between(collect(plat[1]:plat[2]), v[t_pl]+e[t_pl], v[t_pl]-e[t_pl], alpha=0.75, color="green")
-        ylabel(L"$t^2E(x0, t)$")
+        yyy = Y_aux[:, t_pl] * t[nt0]^2 / L^3 ; uwerr.(yyy)
+        # plot(collect(1:xmax), value.(Y_aux[:, t_pl]) * t[nt0]^2 / L^3, "x") 
+        errorbar(collect(1:xmax), value.(yyy), err.(yyy), fmt="d", mfc="none", color="black", capsize=2) 
+        fill_between(collect(plat[1]:plat[2]), v[t_pl]+e[t_pl], v[t_pl]-e[t_pl], alpha=0.75, color="royalblue")
+        ylabel(L"$t^2\langle E(x_0, t)\rangle$")
         xlabel(L"$x_0/a$")
-        ylim(v[t_pl]-20*e[t_pl], v[t_pl]+20*e[t_pl])
+        ylim(v[t_pl]-30*e[t_pl], v[t_pl]+50*e[t_pl])
         title(string(L"$t/a^2 = $", t[nt0]))
+        tight_layout()
+        savefig("/Users/alessandroconigli/Desktop/t0_plateau.pdf")
         display(gcf())
     end
     if info && !isnothing(rw)

src/juobs_reader.jl

@@ -518,4 +518,63 @@ function truncate_data!(data::Vector{Vector{CData}}, nc::Vector{Int64})
         end
     end
     return nothing
-end
\ No newline at end of file
+end
+
+@doc raw"""
+
+    concat_data!(data1::Vector{CData}, data2::Vector{CData})
+
+    concat_data!(data1::Vector{Vector{CData}}, data2::Vector{Vector{CData}})
+
+Concatenates the output of 2 calls to `read_mesons` over configurations. 
+Both data must have the same number of replicas and correlators. 
+The output is saved in the first argument, so if you want to concatenate 
+3 data sets: `concat_data!(data1, data2); concat_data!(data1, data3)`
+
+Examples:
+```@example    
+#Single replica
+dat = read_mesons(path, "G5", "G5")
+dat2 = read_mesons(path2, "G5", "G5")
+concat_data!(dat, dat2)
+
+#Two replicas
+dat = read_mesons([path1, path2], "G5", "G5")
+dat2 = read_mesons([path3, path4], "G5", "G5")
+concat_data!(dat, dat2)
+```
+"""
+function concat_data!(data1::Vector{CData}, data2::Vector{CData})
+    N = length(data1)
+    if length(data1) != length(data2) 
+        error("number of correlators do not match")
+    end
+    for k = 1:N
+        data1[k].vcfg = vcat(data1[k].vcfg, data2[k].vcfg)
+        data1[k].re_data = vcat(data1[k].re_data, data2[k].re_data)
+        data1[k].im_data = vcat(data1[k].im_data, data2[k].im_data)
+    end
+    return nothing
+end
+function concat_data!(data1::Vector{Vector{CData}}, data2::Vector{Vector{CData}})
+    N = length(data1)
+    if length(data1) != length(data2) 
+        error("number of correlators do not match")
+    end
+    R = length(data1[1])
+    if length(data1[1]) != length(data2[1])
+        error("number of replicas do not match")
+    end
+    for k = 1:N
+        for r = 1:R
+            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

src/juobs_tools.jl

 #TODO: apply_rw with gaps
-function apply_rw(data::Array{Float64}, W::Matrix{Float64})
-    nc = size(data, 1)
-    W1 = W[1, 1:nc]
-    W2 = W[2, 1:nc]
+function apply_rw(data::Array{Float64}, W::Matrix{Float64}; vcfg::Union{Nothing, Vector{Int32}}=nothing)
+    nc =  isnothing(vcfg) ? collect(1:size(data, 1)) : vcfg
+    W1 = W[1, nc]
+    W2 = W[2, nc]
 
     data_r = data .* W1 .* W2
     return (data_r, W1 .* W2)
@@ -20,6 +20,16 @@ function apply_rw(data::Vector{<:Array{Float64}}, W::Vector{Matrix{Float64}})
     data_r = [data[k] .* rw[k] for k=1:length(data)]
     return (data_r, rw)
 end
+# rw with gaps
+function apply_rw(data::Vector{<:Array{Float64}}, W::Vector{Matrix{Float64}}, vcfg::Vector{Vector{Int32}})
+    nc = size.(W, 2)
+
+    rw1 = [W[k][1, 1:nc[k]] for k=1:length(W)]
+    rw2 = [W[k][2, 1:nc[k]] for k=1:length(W)]
+    rw = [rw1[k] .* rw2[k] for k =1:length(W)]
+    data_r = [data[k] .* rw[k][vcfg[k]] for k=1:length(data)]
+    return (data_r, rw)
+end
 
 function check_corr_der(obs::Corr, derm::Vector{Corr})   
     g1 = Vector{String}(undef, 0)
@@ -218,7 +228,7 @@ corr_pp = corr_obs.(data)
 corr_pp_r = corr_obs.(data, rw=[rw1, rw2])
 ```
 """
-function corr_obs_TSM(cdata_sl::CData, cdata_corr::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1, info::Bool=false, idm_sl::Union{Vector{Int64},Nothing}=nothing, idm_corr::Union{Vector{Int64},Nothing}=nothing, nms::Int64=Int64(maximum(cdata_sl.vcfg)))
+function corr_obs_TSM(cdata_sl::CData, cdata_corr::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1, info::Bool=false, idm_sl::Union{Vector{Int64},Nothing}=nothing, idm_corr::Union{Vector{Int64},Nothing}=nothing, nms::Union{Int64,Nothing}=nothing)
     # cdata1 is sloppy, cdata2 is correction 
     if cdata_sl.id != cdata_corr.id
         error("Error: cdata_sl id != cdata_corr id")
@@ -227,43 +237,45 @@ function corr_obs_TSM(cdata_sl::CData, cdata_corr::CData; real::Bool=true, rw::U
         error("Error: cdata_sl header != cdata_corr header")
     end
 
-    data1 = real ? cdata_sl.re_data   ./ L^3 : cdata_sl.im_data   ./ L^3
-    data2 = real ? cdata_corr.re_data ./ L^3 : cdata_corr.im_data ./ L^3
+    id = getfield(cdata_sl, :id)
 
-    nt1 = size(data1, 2) # 1
-    nt2 = size(data2, 2)
-    id_sl   = cdata_sl.id
-    id_corr = cdata_corr.id 
-    # added automatic idsm detection - to be tested 
-    nms_vec_sl   = cdata_sl.vcfg
-    nms_vec_corr = cdata_corr.vcfg
+    vcfg_sl   = getfield(cdata_sl, :vcfg)
+    vcfg_corr = getfield(cdata_corr, :vcfg)
+    nms = isnothing(nms) ? Int64(maximum(vcfg_sl)) : nms # assuming vcfg_sl >= vcfg_corr
+    
+    idm_sl = isnothing(idm_sl) ? Int64.(vcfg_sl) : idm_sl
+    idm_corr = isnothing(idm_corr) ? Int64.(vcfg_corr) : idm_corr
+    
 
-    delta_cnfg_sl   = nms_vec_sl[2] - nms_vec_sl[1]
-    delta_cnfg_corr = nms_vec_corr[2] - nms_vec_corr[1]
+    data1 = real ? cdata_sl.re_data   ./ L^3 : cdata_sl.im_data   ./ L^3
+    data2 = real ? cdata_corr.re_data ./ L^3 : cdata_corr.im_data ./ L^3
 
-    idm_sl   = isnothing(idm_sl) ? collect(1:delta_cnfg_sl:nms_vec_sl[end]) : idm_sl
-    idm_corr = isnothing(idm_corr) ? collect(1:delta_cnfg_corr:nms_vec_corr[end]) : idm_corr
+    nt = size(data1, 2) # 1
+    # added automatic idsm detection  
+    # nms_vec_sl   = cdata_sl.vcfg
+    # nms_vec_corr = cdata_corr.vcfg
+ 
+    # delta_cnfg_sl   = nms_vec_sl[2] - nms_vec_sl[1]
+    # delta_cnfg_corr = nms_vec_corr[2] - nms_vec_corr[1]
+ 
+    # idm_sl   = isnothing(idm_sl) ? collect(1:delta_cnfg_sl:nms_vec_sl[end]) : idm_sl
+    # idm_corr = isnothing(idm_corr) ? collect(1:delta_cnfg_corr:nms_vec_corr[end]) : idm_corr
     # end section to be tested 
     if isnothing(rw)
-        idm_sl   = isnothing(idm_sl)   ? collect(1:nms) : idm_sl
-        idm_corr = isnothing(idm_corr) ? collect(1:nms) : idm_corr
-        obs1 = [uwreal(data1[:, x0], id_sl, idm_sl, nms) for x0 = 1:nt1]
-        obs2 = [uwreal(data2[:, x0], id_corr, idm_corr, nms) for x0 = 1:nt2]
-
+        obs1 = [uwreal(data1[:, x0], id, idm_sl, nms) for x0 = 1:nt]
+        obs2 = [uwreal(data2[:, x0], id, idm_corr, nms) for x0 = 1:nt]
     else
-        idm_sl   = isnothing(idm_sl)   ? collect(1:nms) : idm_sl
-        idm_corr = isnothing(idm_corr) ? collect(1:nms) : idm_corr
 
-        data1_r, W_sl = apply_rw(data1, rw)
-        data2_r, W_corr = apply_rw(data2, rw)
-        ow1 = [uwreal(data1_r[:, x0], id_sl, idm_sl, nms) for x0 = 1:nt1]
-        ow2 = [uwreal(data2_r[:, x0], id_corr, idm_corr, nms) for x0 = 1:nt2]
+        data1_r, W_sl = apply_rw(data1, rw, vcfg=vcfg_sl)
+        data2_r, W_corr = apply_rw(data2, rw, vcfg=vcfg_corr)
+        ow1 = [uwreal(data1_r[:, x0], id, idm_sl, nms) for x0 = 1:nt]
+        ow2 = [uwreal(data2_r[:, x0], id, idm_corr, nms) for x0 = 1:nt]
 
-        W_obs_sl   = uwreal(W_sl, id_sl, idm_sl, nms)
-        W_obs_corr = uwreal(W_corr, id_corr, idm_corr, nms)
+        W_obs_sl   = uwreal(W_sl, id, idm_sl, nms)
+        W_obs_corr = uwreal(W_corr, id, idm_corr, nms)
 
-        obs1 = [ow1[x0] / W_obs_sl for x0 = 1:nt1]
-        obs2 = [ow2[x0] / W_obs_corr for x0 = 1:nt2]
+        obs1 = [ow1[x0] / W_obs_sl for x0 = 1:nt]
+        obs2 = [ow2[x0] / W_obs_corr for x0 = 1:nt]
     end
 
     if info && !isnothing(rw)
@@ -318,26 +330,25 @@ function corr_obs_TSM(cdata1::Array{CData, 1}, cdata2::Array{CData, 1}; real::Bo
     data2 = real ? getfield.(cdata2, :re_data) ./ L^3 : getfield.(cdata2, :im_data) ./ L^3
 
     nt = size(data1[1], 2)
-
     if isnothing(rw)
         tmp1 = cat(data1..., dims=1)
         tmp2 = cat(data2..., dims=1)
 
         obs1 = [uwreal(tmp1[:, x0], id[1], replica_sl, idm_sl, nms) for x0 = 1:nt]
-        obs2 = [uwreal(tmp2[:, x0], id[1], replica_corr, idm_corr, nms) for x0 = 1:nt]
+        obs2 = [uwreal(tmp2[:, x0], id[1], replica_sl, idm_corr, nms) for x0 = 1:nt]
 
     else
-        data1_r, W = apply_rw(data1, rw)
-        data2_r, W = apply_rw(data2, rw)
+        data1_r, W = apply_rw(data1, rw, vcfg_sl)
+        data2_r, W = apply_rw(data2, rw, vcfg_corr)
 
         tmp1 = cat(data1_r..., dims=1)
         tmp2 = cat(data2_r..., dims=1)
         tmp_W = cat(W..., dims=1)
 
         ow1 = [uwreal(tmp1[:, x0], id[1], replica_sl, idm_sl, nms) for x0 = 1:nt]
-        ow2 = [uwreal(tmp2[:, x0], id[1], replica_corr, idm_corr, nms) for x0 = 1:nt]
+        ow2 = [uwreal(tmp2[:, x0], id[1], replica_sl, idm_corr, nms) for x0 = 1:nt]
 
-        W_obs = uwreal(tmp_W, id[1], replica_sl, idm_sl, nms  )
+        W_obs = uwreal(tmp_W, id[1], replica_sl, idm_sl, nms)
 
         obs1 = [ow1[x0] / W_obs for x0 = 1:nt]
         obs2 = [ow2[x0] / W_obs for x0 = 1:nt]
@@ -645,7 +656,7 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
     
     for INDEX in tmin_array ## vary tmin
         for j in tmax_array ## vary tmax
-             try
+            try
                 x = [i for i in INDEX+1:1:j]	
 	            yy = y[INDEX+1:1:j]
 	            Ncut = total - length(x)
@@ -672,7 +683,7 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
                 push!(p1, up[1])
                 push!(mods,string("[", INDEX+1, ",", j, "]"))
                 # println("chi2_vs_exp " , chi2 / dof(fit),  " chi2/dof ", sum(fit.resid.^2)/dof(fit))
-            catch e
+            catch 
                 @warn string(":/ Negative window for error propagation at tmin = ", INDEX, ", tmax = ", j, "; skipping that point")
             end
         end
@@ -1448,7 +1459,6 @@ function get_covar(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},D
     tmp_fluc = mchist.(obs, ens)
     fluc = mapreduce(permutedims, vcat, tmp_fluc)' 
     covar = fluc' * fluc ./ size(fluc)[1]^2
-    #@warn("should double check it. no justification to have a square in the Ncnfg")
     return covar
 end
 function get_covar_multi_id(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)