bug flag_strange fixed, added a few constants and reader for TSM and generic...

bug flag_strange fixed, added a few constants and reader for TSM and generic wil & tm corrs, needs checking but should not clash with anything anyway

constants/juobs_const.jl

@@ -140,6 +140,14 @@ const db_c = Dict(
     "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],

src/juobs.jl

@@ -12,10 +12,10 @@ include("../constants/juobs_const.jl")
 include("juobs_types.jl")
 
 export EnsInfo
-export ens_db, db, db_c, flag_s
+export ens_db, db, db_c, flag_s, sym_bool
 export read_mesons, read_mesons_correction, read_mesons_multichunks, read_mesons_correction_multichunks, 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, pvalue, fve, model_av, fit_alg
+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_reader.jl

@@ -1002,3 +1002,182 @@ function read_mesons_correction_multichunks(path::String, ens::String, g1::Strin
 
     return dat
 end
+
+function get_corr_TSM_multichunks(path::String, ens::EnsInfo; info=false)
+    path = joinpath(path, ens.id)
+    path_sl = joinpath.(path, "sloppy")
+    path_c = joinpath.(path, "correc")
+    path_rw = joinpath(path, "rwf_def")
+    path_rw = filter(x->occursin(".dat", x), readdir(path_rw, join=true))
+    if length(path_rw) == 0
+        path_rw = joinpath(path, "rwf")
+        global path_rw = filter(x->occursin(".dat", x), readdir(path_rw, join=true))
+    end
+    
+    pp_dat = read_mesons_multichunks(path_sl, ens.id, "G5", "G5")
+    ap_dat = read_mesons_multichunks(path_sl, ens.id, "G5", "G0G5")
+    pp_dat_c = read_mesons_correction_multichunks(path_c, ens.id, "G5", "G5")
+    ap_dat_c = read_mesons_correction_multichunks(path_c, ens.id, "G5", "G0G5")
+
+    rwf = [read_ms1(path_rw[i], v=ens.vrw[i]) for i in 1:length(ens.vrw)]
+    cnfg_rw = size.(rwf,2)
+    cnfg_trunc_ts001 = [findall(pp_dat[1][i].vcfg .< cnfg_rw[i])[end] for i in 1:length(cnfg_rw)] ## rwf missing some configs at the end
+    cnfg_trunc_ts001_c = [findall(pp_dat_c[1][i].vcfg .< cnfg_rw[i])[end] for i in 1:length(cnfg_rw)]
+    cnfg_trunc_ts190 = [findall(pp_dat[2][i].vcfg .< cnfg_rw[i])[end] for i in 1:length(cnfg_rw)] ## rwf missing some configs at the end
+    cnfg_trunc_ts190_c = [findall(pp_dat_c[2][i].vcfg .< cnfg_rw[i])[end] for i in 1:length(cnfg_rw)]
+    truncate_data!(pp_dat[1:2:end], cnfg_trunc_ts001)
+    truncate_data!(ap_dat[1:2:end], cnfg_trunc_ts001)
+    truncate_data!(pp_dat_c[1:2:end], cnfg_trunc_ts001_c)
+    truncate_data!(ap_dat_c[1:2:end], cnfg_trunc_ts001_c)
+    truncate_data!(pp_dat[2:2:end], cnfg_trunc_ts190) 
+    truncate_data!(ap_dat[2:2:end], cnfg_trunc_ts190)
+    truncate_data!(pp_dat_c[2:2:end], cnfg_trunc_ts190_c)
+    truncate_data!(ap_dat_c[2:2:end], cnfg_trunc_ts190_c)
+
+    if sym_bool[ens.id] == true
+        pp = corr_obs_TSM.(pp_dat[1:length(ap_dat_c)], pp_dat_c[1:length(ap_dat_c)], rw=rwf, L=ens.L, info=info, replica_sl=ens.cnfg, nms=sum(ens.cnfg))
+        ap = corr_obs_TSM.(ap_dat[1:length(ap_dat_c)], ap_dat_c[1:length(ap_dat_c)], rw=rwf, L=ens.L, info=info, replica_sl=ens.cnfg, nms=sum(ens.cnfg))
+    
+        if ens.id == "E250"
+            pp_sym = [corr_sym_E250(pp[i], pp[i+1], +1) for i in 1:2:length(ap)]
+            ap_sym = [corr_sym_E250(ap[i], ap[i+1], -1) for i in 1:2:length(ap)]
+        elseif ens.id == "D450"
+            pp_sym = [corr_sym_D450(pp_ts001[i], pp_tsT[i], +1) for i in 1:length(pp_ts001)]
+            ap_sym = [corr_sym_D450(ap_ts001[i], ap_tsT[i], -1) for i in 1:length(pp_ts001)]
+        else
+            pp_sym = [corr_sym(pp[i], pp[i+1], +1) for i in 1:2:length(ap)]
+            ap_sym = [corr_sym(ap[i], ap[i+1], -1) for i in 1:2:length(ap)]
+        end
+    else
+        pp_ts001 = corr_obs_TSM.(pp_dat[1:2:length(ap_dat)], pp_dat_c[1:length(ap_dat_c)], rw=rwf, L=ens.L, info=info, replica_sl=ens.cnfg, nms=sum(ens.cnfg))
+        ap_ts001 = corr_obs_TSM.(ap_dat[1:2:length(ap_dat)], ap_dat_c[1:length(ap_dat_c)], rw=rwf, L=ens.L, info=info, replica_sl=ens.cnfg, nms=sum(ens.cnfg))
+        pp_tsT = corr_obs.(pp_dat[2:2:length(ap_dat)], rw=rwf, L=ens.L, info=info, replica=ens.cnfg, nms=sum(ens.cnfg))
+        ap_tsT = corr_obs.(ap_dat[2:2:length(ap_dat)], rw=rwf, L=ens.L, info=info, replica=ens.cnfg, nms=sum(ens.cnfg))
+    
+        if ens.id == "E250"
+            pp_sym = [corr_sym_E250(pp_ts001[i], pp_tsT[i], +1) for i in 1:length(pp_ts001)]
+            ap_sym = [corr_sym_E250(ap_ts001[i], ap_tsT[i], -1) for i in 1:length(pp_ts001)]
+        elseif ens.id == "D450"
+            pp_sym = [corr_sym_D450(pp_ts001[i], pp_tsT[i], +1) for i in 1:length(pp_ts001)]
+            ap_sym = [corr_sym_D450(ap_ts001[i], ap_tsT[i], -1) for i in 1:length(pp_ts001)]
+        else
+            pp_sym = [corr_sym(pp_ts001[i], pp_tsT[i], +1) for i in 1:length(pp_ts001)]
+            ap_sym = [corr_sym(ap_ts001[i], ap_tsT[i], -1) for i in 1:length(pp_ts001)]
+        end
+    end
+    
+    return pp_sym, ap_sym
+end
+
+function get_corr_TSM(path::String, ens::EnsInfo, g1::String, g2::String; rw=false, info=false, legacy=false, fs=false)
+    path_rw = joinpath(path, ens.id, "rwf_def")
+    path_rw = filter(x->occursin(".dat", x), readdir(path_rw, join=true))
+    if length(path_rw) == 0
+        path_rw = joinpath(path, ens.id, "rwf")
+        global path_rw = filter(x->occursin(".dat", x), readdir(path_rw, join=true))
+    end
+    path_sl = joinpath(path, ens.id, "sloppy")
+    path_sl = filter(x->occursin(".mesons.dat", x), readdir(path_sl, join=true))
+    path_c = joinpath(path, ens.id, "correc")
+    path_c = filter(x->occursin(".mesons.dat", x), readdir(path_c, join=true))
+
+    rwf = read_ms1.(path_rw, v=ens.vrw)
+    dat_sl = read_mesons([path_sl[i] for i in 1:length(path_sl)], g1, g2, legacy=legacy, id=ens.id)
+    dat_c = read_mesons_correction([path_c[i] for i in 1:length(path_c)], g1, g2, legacy=legacy, id=ens.id)
+
+    rw ? corr = [corr_obs_TSM(dat_sl[i], dat_c[i], L=ens.L, rw=rwf, info=info) for i in 1:length(dat_sl)] : corr = [corr_obs_TSM(dat_sl[i], dat_c[i], L=ens.L, info=info) for i in 1:length(dat_sl)]
+
+    if info == false
+        return corr
+    else
+        pp = getindex.(corr,1)
+        ppw = getindex.(corr,2)
+        w = getindex.(corr,3)
+        return pp, ppw, w[1]
+    end
+end
+
+function get_corr_wil(path::String, ens::EnsInfo, g1::String, g2::String; rw=false, info=false, legacy=false, fs=false)
+    path_rw = joinpath(path, ens.id, "rwf_def")
+    path_rw = filter(x->occursin(".dat", x), readdir(path_rw, join=true))
+    if length(path_rw) == 0
+        path_rw = joinpath(path, ens.id, "rwf")
+        global path_rw = filter(x->occursin(".dat", x), readdir(path_rw, join=true))
+    end
+    path = joinpath(path, ens.id, "wil")
+    path = filter(x->occursin(".mesons.dat", x), readdir(path, join=true))
+
+    if ens.id == "D200"
+        #rwf_1 = read_ms1.([path_rw[1]], v=ens.vrw)
+        #rwf_2 = read_ms1.([path_rw[2]], v=ens.vrw)
+        #rwf = [hcat(rwf_1[1],rwf_2[1])]
+        rwf = read_ms1.(path_rw, v=ens.vrw)
+        dat = read_mesons([path[1]], g1, g2, legacy=legacy, id=ens.id)
+        dat_2 = read_mesons([path[2]], g1, g2, legacy=legacy, id=ens.id)
+        concat_data!(dat,dat_2)
+    else
+        rwf = read_ms1.(path_rw, v=ens.vrw)
+        dat = read_mesons([path[i] for i in 1:length(path)], g1, g2, legacy=legacy, id=ens.id)
+        truncate_data!(dat,ens.cnfg)
+    end
+
+    rw ? corr = [corr_obs(dat[i], L=ens.L, rw=rwf, info=info, flag_strange=fs) for i in 1:length(dat)] : corr = [corr_obs(dat[i], L=ens.L, info=info) for i in 1:length(dat)]
+
+    if info == false
+        return corr
+    else
+        pp = getindex.(corr,1)
+        ppw = getindex.(corr,2)
+        w = getindex.(corr,3)
+        return pp, ppw, w[1]
+    end
+end
+
+function get_corr_tm(path::String, ens::EnsInfo, g1::String, g2::String; rw=false, info=false, legacy=false, fs=false)
+    path_rw = joinpath(path, ens.id, "rwf_def")
+    path_rw = filter(x->occursin(".dat", x), readdir(path_rw, join=true))
+    if length(path_rw) == 0
+        path_rw = joinpath(path, ens.id, "rwf")
+        global path_rw = filter(x->occursin(".dat", x), readdir(path_rw, join=true))
+    end
+    path = joinpath(path, ens.id, "tm")
+    path = filter(x->occursin(".mesons.dat", x), readdir(path, join=true))
+
+    if ens.id == "J303"
+        rwf = read_ms1.(path_rw, v=ens.vrw)
+        dat = read_mesons([path[1]], g1, g2, legacy=legacy, id=ens.id)
+        dat_2 = read_mesons([path[2]], g1, g2, legacy=legacy, id=ens.id)
+        concat_data!(dat,dat_2)
+    elseif ens.id == "D200"
+        #rwf_1 = read_ms1.([path_rw[1]], v=ens.vrw)
+        #rwf_2 = read_ms1.([path_rw[2]], v=ens.vrw)
+        #rwf = [hcat(rwf_2[1],rwf_1[1])]
+        rwf = read_ms1.(path_rw, v=ens.vrw)
+        dat = read_mesons([path[1]], g1, g2, legacy=legacy, id=ens.id)
+        dat_2 = read_mesons([path[2]], g1, g2, legacy=legacy, id=ens.id)
+        dat_3 = read_mesons([path[3]], g1, g2, legacy=legacy, id=ens.id)
+        concat_data!(dat,dat_3)
+        concat_data!(dat,dat_2)
+    elseif ens.id == "N300"
+        rwf = read_ms1.(path_rw, v=ens.vrw)
+        dat_1 = read_mesons([path[1]], g1, g2, legacy=legacy, id=ens.id)
+        dat = read_mesons([path[2]], g1, g2, legacy=legacy, id=ens.id)
+        truncate_data!(dat,[199])
+        concat_data!(dat,dat_1)
+    else
+        rwf = read_ms1.(path_rw, v=ens.vrw)
+        dat = read_mesons([path[i] for i in 1:length(path)], g1, g2, legacy=legacy, id=ens.id)
+        truncate_data!(dat,ens.cnfg)
+    end
+
+    rw ? corr = [corr_obs(dat[i], L=ens.L, rw=rwf, info=info) for i in 1:length(dat)] : corr = [corr_obs(dat[i], L=ens.L, info=info) for i in 1:length(dat)]
+
+    if info == false
+        return corr
+    else
+        pp = getindex.(corr,1)
+        ppw = getindex.(corr,2)
+        w = getindex.(corr,3)
+        return pp, ppw, w[1]
+    end
+end
\ No newline at end of file

src/juobs_tools.jl

@@ -429,6 +429,25 @@ function corr_sym(corrL::Corr, corrR::Corr, parity::Int64=1)
     res = (corrL.obs[1:end] + parity * corrR.obs[end:-1:1]) / 2
     return Corr(res, corrL.kappa, corrL.mu, corrL.gamma, corrL.y0, corrL.theta1, corrL.theta2)
 end
+
+function corr_sym_E250(corr1::Corr, corr2::Corr, parity::Int64=1)
+    aux = [corr2.obs[97:end]; corr2.obs[1:96]]
+    corr2_sym = Corr(aux, corr2.kappa, corr2.mu, corr2.gamma, corr1.y0, corr2.theta1, corr2.theta2)
+
+    corr = [corr1.obs[1:3]; (corr1.obs[4:98] .+ parity * corr2_sym.obs[192:-1:98]) / 2]
+
+    return Corr(corr, corr1.kappa, corr1.mu, corr1.gamma, corr1.y0, corr1.theta1, corr1.theta2)
+end
+
+function corr_sym_D450(corr1::Corr, corr2::Corr, parity::Int64=1)
+    aux = [corr2.obs[65:end]; corr2.obs[1:64]]
+    corr2_sym = Corr(aux, corr2.kappa, corr2.mu, corr2.gamma, corr1.y0, corr2.theta1, corr2.theta2)
+
+    corr = [corr1.obs[1:3]; (corr1.obs[4:66] .+ parity * corr2_sym.obs[128:-1:66]) / 2]
+
+    return Corr(corr, corr1.kappa, corr1.mu, corr1.gamma, corr1.y0, corr1.theta1, corr1.theta2)
+end
+
 #TODO: VECTORIZE, uwreal?
 @doc raw"""
     md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::uwreal, w::Union{uwreal, Nothing}=nothing, ws::ADerrors.wspace=ADerrors.wsg, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)