Skip to content

J
juobs
Commit 619245f8 authored by ale's avatar ale
Browse files
Options

added get_* obs with model av incorporated, w0 routine, model_av update,... 

added get_* obs with model av incorporated, w0 routine, model_av update, fit_alg allows correlated fits, EnsInfo and some const, pvalue with warning and 1.0 safe, read_chunks for TSM. There are duplicities, CHECK. Needs checking and DOCUMENTING
parent 95123bce
Hide whitespace changes
Inline Side-by-side
Showing with 1487 additions and 6 deletions
constants/juobs_const.jl
View file @ 619245f8
...@@ -99,4 +99,43 @@ zv(beta::Float64) = Zv[b_values .== beta][1] ...@@ -99,4 +99,43 @@ zv(beta::Float64) = Zv[b_values .== beta][1]
zs_over_zp(beta::Float64) = ZS_over_ZP[b_values .== beta][1] zs_over_zp(beta::Float64) = ZS_over_ZP[b_values .== beta][1]
rm(beta::Float64) = RM[b_values .== beta][1] rm(beta::Float64) = RM[b_values .== beta][1]
ba_bp(beta::Float64) = BA_BP[b_values .== beta][1] ba_bp(beta::Float64) = BA_BP[b_values .== beta][1]
Z(beta::Float64) = ZZ[b_values .== beta][1] Z(beta::Float64) = ZZ[b_values .== beta][1]
\ No newline at end of file
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"]
)
\ No newline at end of file
src/juobs.jl
View file @ 619245f8
...@@ -9,9 +9,9 @@ include("juobs_reader.jl") ...@@ -9,9 +9,9 @@ include("juobs_reader.jl")
include("juobs_tools.jl") include("juobs_tools.jl")
include("juobs_obs.jl") include("juobs_obs.jl")
export read_mesons, read_mesons_correction, read_ms1, read_ms, read_md, truncate_data!, concat_data! 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 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 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 meff, mpcac, dec_const, dec_const_w, dec_const_pcvc, comp_t0 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 end # module
src/juobs_obs.jl
View file @ 619245f8
...@@ -828,4 +828,838 @@ function t0_guess(t::Vector{Float64}, Ysl::Array{Float64, 3}, plat::Vector{Int64 ...@@ -828,4 +828,838 @@ function t0_guess(t::Vector{Float64}, Ysl::Array{Float64, 3}, plat::Vector{Int64
t0_aux = minimum(abs.(t2E_ax .- 0.3)) t0_aux = minimum(abs.(t2E_ax .- 0.3))
nt0 = findfirst(x-> abs(x - 0.3) == t0_aux, t2E_ax) #index closest value to t0 nt0 = findfirst(x-> abs(x - 0.3) == t0_aux, t2E_ax) #index closest value to t0
return nt0 return nt0
end
function get_m(corr::juobs.Corr, ens::EnsInfo, PS::String;
tm::Union{Vector{Vector{Int64}}, Nothing}=nothing, tM::Union{Vector{Vector{Int64}}, Nothing}=nothing,
pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
corr_d = corr.obs
m_dat = 0.5 .* log.((corr_d[2:end-2] ./ corr_d[3:end-1]) .^ 2)
y0 = corr.y0
T = length(corr_d) - 1 - y0
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(m_dat[Int64(round(T / 2))])
m, syst, m_i, weight, pval = model_av([fit_exp, fit_const], m_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
if pl == true
isnothing(wpm) ? uwerr(m) : uwerr(m, wpm)
isnothing(wpm) ? uwerr.(m_dat) : [uwerr(m_dat[i], wpm) for i in 1:length(m_dat)]
isnothing(wpm) ? uwerr.(m_i) : [uwerr(m_i[i], wpm) for i in 1:length(m_i)]
v = value(m)
e = err(m)
fig = figure("eff")
errorbar(1:length(m_dat), value.(m_dat), err.(m_dat), fmt="x", color="black")
ylabel(L"$am_\mathrm{eff}$")
xlabel(L"$x_0/a$")
ylim(v-7*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/m_",ens.id,"_",PS,"_plat.pdf"))
models = Array{String,1}()
for i in 1:length(tm[1])
for j in 1:length(tM[1])
push!(models, string("[",tm[1][i],",",tM[1][j],"]"))
end
end
for i in 1:length(tm[2])
for j in 1:length(tM[2])
push!(models, string("[",tm[2][i],",",tM[2][j],"]"))
end
end
fig = figure("model av")
subplot(411)
fill_between(1:length(m_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(m_dat), value.(m_dat), err.(m_dat), fmt="x", color="black")
ylabel(L"$am_\mathrm{eff}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+20*e)
ax = gca()
setp(ax.get_xticklabels(),visible=false)
subplot(412)
ylabel(L"$m_i$")
fill_between(1:length(m_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(m_i), value.(m_i), err.(m_i), fmt="x", color="black")
ax = gca()
ax[:set_xlim]([0, length(models)+1])
setp(ax.get_xticklabels(),visible=false)
subplot(413)
ylabel(L"$W_i$")
bar(1:length(m_i), weight, color="green")
ax = gca()
ax[:set_xlim]([0, length(models)+1])
setp(ax.get_xticklabels(),visible=false)
subplot(414)
ylabel("p-value")
bar(1:length(m_i), pval, color="green")
ax = gca()
ax[:set_xlim]([0, length(models)+1])
plt.xticks(collect(1:length(m_i)), models)
xticks(rotation=90)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/m_",ens.id,"_",PS,".pdf"))
#close("all")
end
return m, syst, m_i, weight, pval
end
function get_m_pbc(corr::juobs.Corr, ens::EnsInfo, PS::String;
tm::Union{Vector{Int64}, Nothing}=nothing, tM::Union{Vector{Int64}, Nothing}=nothing,
pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
method::String="cosh")
corr_d = corr.obs
m_dat = 0.5 .* log.((corr_d[2:end-2] ./ corr_d[3:end-1]) .^ 2)
if ens.id == "E250"
global T = 192
global T2 = 192
global Thalf = 96
elseif ens.id == "D450"
global T = 128
global T2 = 128
global Thalf = 64
end
guess = value(m_dat[Int64(round(T / 3))])
isnothing(tm) ? tm = [y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40] : tm=tm
isnothing(tM) ? tM = [Thalf] : tM=tM
aux = [corr_d[i] / corr_d[i+1] for i in 2:length(corr_d)-1]
if ens.id == "E250"
global T = 96
elseif ens.id == "D450"
global T = 64
end
@.fit_fun(x,p) = cosh(p[1] * (x-T)) / cosh(p[1] * (x+1-T))
k1 = 1
aux2 = corr_d
@.fit_fun2(x,p) = p[2] * exp(-p[1] * (x-1)) + p[2] * exp(-p[1] * (T2+1-x))
k2 = 2
if method == "cosh"
m, syst, m_i, weight, pval = model_av(fit_fun, aux, guess, tm=tm, tM=tM, k=k1, wpm=wpm)
elseif method == "corr"
m, syst, m_i, weight, pval = model_av(fit_fun2, aux2, guess, tm=tm, tM=tM, k=k2, wpm=wpm)
end
if pl == true
##TODO
bla = 1
end
return m, syst, m_i, weight, pval
end
function get_mpcac(corr_pp::juobs.Corr, corr_ap::juobs.Corr, ens::EnsInfo, PS::String;
tm::Union{Vector{Vector{Int64}}, Nothing}=nothing, tM::Union{Vector{Vector{Int64}}, Nothing}=nothing,
impr::Bool=true, pl::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
ap_dat = -corr_ap.obs[2:end-1]
pp_dat = corr_pp.obs[2:end-1]
der_ap = (ap_dat[3:end] .- ap_dat[1:end-2]) / 2
if impr == true
ca = ens.ca
der2_pp = pp_dat[1:end-2] + pp_dat[3:end] - 2*pp_dat[2:end-1]
der_ap = der_ap + ca*der2_pp
end
mpcac_dat = der_ap ./ (2*pp_dat[2:end-1])
y0 = corr_pp.y0
T = length(pp_dat) - 1 - y0
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(mpcac_dat[Int64(round(T / 2))])
mpcac, syst, mpcac_i, weight, pval = model_av([fit_exp, fit_const], mpcac_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
if pl == true
isnothing(wpm) ? uwerr(mpcac) : uwerr(mpcac, wpm)
isnothing(wpm) ? uwerr.(mpcac_dat) : [uwerr(mpcac_dat[i], wpm) for i in 1:length(mpcac_dat)]
isnothing(wpm) ? uwerr.(mpcac_i) : [uwerr(mpcac_i[i], wpm) for i in 1:length(mpcac_i)]
v = value(mpcac)
e = err(mpcac)
fig = figure("eff")
errorbar(1:length(mpcac_dat), value.(mpcac_dat), err.(mpcac_dat), fmt="x", color="black")
ylabel(L"$am_\mathrm{pcac}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/mpcac_",ens.id,"_",PS,"_plat.pdf"))
fig = figure("model av")
subplot(411)
fill_between(1:length(mpcac_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(mpcac_dat), value.(mpcac_dat), err.(mpcac_dat), fmt="x", color="black")
ylabel(L"$am_\mathrm{pcac}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+10*e)
subplot(412)
ylabel(L"$mpcac_i$")
fill_between(1:length(mpcac_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(mpcac_i), value.(mpcac_i), err.(mpcac_i), fmt="x", color="black")
subplot(413)
ylabel(L"$W_i$")
bar(1:length(mpcac_i), weight, color="green")
subplot(414)
ylabel("p-value")
bar(1:length(mpcac_i), pval, color="green")
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/mpcac_",ens.id,"_",PS,".pdf"))
#close("all")
end
return mpcac, syst, mpcac_i, weight, pval
end
function get_f_wil(corr_pp::juobs.Corr, corr_ap::juobs.Corr, m::uwreal, ens::EnsInfo, PS::String;
tm::Union{Vector{Vector{Int64}}, Nothing}=nothing, tM::Union{Vector{Vector{Int64}}, Nothing}=nothing,
impr::Bool=true, pl::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
ap_dat = -corr_ap.obs
pp_dat = corr_pp.obs
T = length(pp_dat)
y0 = corr_pp.y0
if impr == true
ca = ens.ca
der_pp = (pp_dat[3:end] .- pp_dat[1:end-2]) / 2
ap_dat = ap_dat[2:end-1] + ca * der_pp
aux = exp.((collect(1:T-2) .- (T-1)/2) .* [m for k = 1:T-2])
else
aux = exp.((collect(0:T-1) .- (T-1)/2) .* [m for k = 1:T])
end
R_dat = ap_dat .* aux ./ [((pp_dat[T-y0])^2)^(1/4) for k = 1:length(ap_dat)]
#f_dat = [sqrt(2) * sqrt(R_dat[i] ^ 2) / sqrt(m) for i in 1:length(R_dat)]
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-1-y0-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(R_dat[Int64(round(T / 2))])
R, syst, R_i, weight, pval = model_av([fit_exp, fit_const], R_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
f = sqrt(2) * sqrt(R^2) / sqrt(m)
if pl == true
isnothing(wpm) ? uwerr(R) : uwerr(R, wpm)
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
isnothing(wpm) ? uwerr.(R_i) : [uwerr(R_i[i], wpm) for i in 1:length(R_i)]
v = value(R)
e = err(R)
fig = figure("eff")
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-20*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
fig = figure("model av")
subplot(411)
fill_between(1:length(R_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+10*e)
subplot(412)
ylabel(L"$R_i$")
fill_between(1:length(R_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_i), value.(R_i), err.(R_i), fmt="x", color="black")
subplot(413)
ylabel(L"$W_i$")
bar(1:length(R_i), weight, color="green")
subplot(414)
ylabel("p-value")
bar(1:length(R_i), pval, color="green")
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,".pdf"))
#close("all")
end
return f, syst, R_i, weight, pval
end
function get_f_wil(corr_ppL::juobs.Corr, corr_ppR::juobs.Corr, corr_apL::juobs.Corr, corr_apR::juobs.Corr,
m::uwreal, ens::EnsInfo, PS::String; tm::Union{Vector{Int64}, Nothing}=nothing,
tM::Union{Vector{Int64}, Nothing}=nothing, impr::Bool=true, pl::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
pp_dat = (corr_ppL.obs + corr_ppR.obs[end:-1:1]) / 2
T = length(corr_ppL.obs)
y0 = corr_ppL.y0
if impr == true
ca = ens.ca
der_ppL = (corr_ppL.obs[3:end] - corr_ppL.obs[1:end-2]) / 2
der_ppR = (corr_ppR.obs[3:end] - corr_ppR.obs[1:end-2]) / 2
apL_dat = -corr_apL.obs[2:end-1] + ca * der_ppL
apR_dat = -corr_apR.obs[2:end-1] + ca * der_ppR
else
apL_dat = -corr_apL.obs[2:end-1]
apR_dat = -corr_apR.obs[2:end-1]
end
f1 = [pp_dat[T - y0] for k = 1:length(apL_dat)]
R_dat = ((apL_dat .* apR_dat ./ f1).^2).^(1/4)
#f_dat = [sqrt(2) * sqrt(R_dat[i] ^ 2) / sqrt(m) for i in 1:length(R_dat)]
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-1-y0-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(R_dat[Int64(round(T / 2))])
R, syst, R_i, weight, pval = model_av([fit_exp, fit_const], R_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
f = sqrt(2) * sqrt(R^2) / sqrt(m)
if pl == true
isnothing(wpm) ? uwerr(R) : uwerr(R, wpm)
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
isnothing(wpm) ? uwerr.(R_i) : [uwerr(R_i[i], wpm) for i in 1:length(R_i)]
v = value(R)
e = err(R)
fig = figure("eff")
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-20*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
fig = figure("model av")
subplot(411)
fill_between(1:length(R_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+10*e)
subplot(412)
ylabel(L"$R_i$")
fill_between(1:length(R_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_i), value.(R_i), err.(R_i), fmt="x", color="black")
subplot(413)
ylabel(L"$W_i$")
bar(1:length(R_i), weight, color="green")
subplot(414)
ylabel("p-value")
bar(1:length(R_i), pval, color="green")
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,".pdf"))
#close("all")
end
return f, syst, R_i, weight, pval
end
function get_f_wil_pbc(corr_pp::juobs.Corr, corr_ap::juobs.Corr, m::uwreal, ens::EnsInfo, PS::String;
tm::Union{Vector{Int64}, Nothing}=nothing, tM::Union{Vector{Int64}, Nothing}=nothing,
impr::Bool=true, pl::Bool=false, method::String="ratio",
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
ap_dat = -corr_ap.obs
pp_dat = corr_pp.obs
if ens.id == "E250"
global T = 192
global Thalf = 96
elseif ens.id == "D450"
global T = 128
global Thalf = 64
end
if impr == true
ca = ens.ca
der_pp = (pp_dat[3:end] .- pp_dat[1:end-2]) / 2
ap_dat = ap_dat[2:end-1] + ca * der_pp
pp_dat = pp_dat[2:end-1]
end
R_dat = ap_dat ./ sqrt.(pp_dat)
isnothing(tm) ? tm = [y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40] : tm=tm
isnothing(tM) ? tM = [Thalf] : tM=tM
@.fit_fun(x,p) = p[1] * (-exp(-value(m) * (x-1)) + exp(-value(m) * (T+1-x))) / sqrt(exp(-value(m) * (x-1)) + exp(-value(m) * (T+1-x)))
@.fit_fun_ap(x,p) = p[1] * (-exp(-value(m) * (x-1)) + exp(-value(m) * (T+1-x)))
@.fit_fun_pp(x,p) = p[1] * (exp(-value(m) * (x-1)) + exp(-value(m) * (T+1-x)))
@.fit_const(x,p) = p[1] * x ^ 0
k = 1
function fit_corr_sim(x,p)
f = [p[1] * (-exp(-value(m) * x[i]) + exp(-value(m) * (T - x[i]))) for i in 1:div(length(x),2)]
g = [p[2] * (exp(-value(m) * x[i]) + exp(-value(m) * (T - x[i]))) for i in div(length(x),2)+1:length(x)]
return [f;g]
end
if method == "ratio"
R, syst, R_i, weight, pval = model_av(fit_fun, R_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
f_i = sqrt.(2 ./ [m for i in 1:length(R_i)]) .* R_i
f = sqrt(2 / m) * R
elseif method == "corr"
cap, syst_ap, cap_i, weight_ap, pval_ap = model_av(fit_fun_ap, ap_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
cpp, syst_pp, cpp_i, weight_pp, pval_pp = model_av(fit_fun_pp, pp_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
f_i = sqrt.(2 ./ [m for i in 1:length(cap_i)]) .* cap_i ./ sqrt.(cpp_i)
f = sqrt(2 / m) * cap / sqrt(cpp)
syst = [syst_ap, syst_pp]
pval = [pval_ap, pval_pp]
weight = [weight_ap, weight_pp]
elseif method == "corr_sim"
pval = Array{Float64,1}()
TIC = Array{Float64,1}()
cap_i = Array{uwreal,1}()
cpp_i = Array{uwreal,1}()
for i in tm
y = [ap_dat[i:end]; pp_dat[i:end]]
x = collect(i:length(ap_dat))
x = [x; x]
uwerr.(y)
W = 1 ./ ADerrors.err.(y) .^ 2
uprm, chi2, chi_exp, pv = fit_alg(fit_corr_sim, x, y, 2, [.5, .5])
push!(cap_i, uprm[1])
push!(cpp_i, uprm[2])
push!(pval, pv)
push!(TIC, chi2 - 2 * chi_exp)
end
weight = exp.(-0.5 * TIC) ./ sum(exp.(-0.5 * TIC))
R_i = cap_i ./ sqrt.(cpp_i)
R_av = sum(R_i .* weight)
syst = sum(R_i .^ 2 .* weight) - R_av ^ 2
f_i = sqrt.(2 ./ [m for i in 1:length(cap_i)]) .* R_i
f = sqrt(2 / m) * R_av
elseif method == "pcac"
if impr == false
error("pcac method implemented only for improved axial current, set impr=true")
end
cpp, syst_pp, cap_i, weight, pval = model_av(fit_fun_pp, pp_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
der_ap = (ap_dat[3:end] .- ap_dat[1:end-2]) / 2
der2_pp = pp_dat[1:end-2] + pp_dat[3:end] - 2*pp_dat[2:end-1]
der_ap = der_ap + ca*der2_pp
mpcac_dat = der_ap ./ (2*pp_dat[2:end-1])
f_dat = [sqrt(8 / m ^ 3 * cpp) for i in 1:length(mpcac_dat)] .* mpcac_dat
f, syst, f_i, weight, pval = model_av(fit_const, f_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
end
if pl == true
##TODO
bla = 1
x = collect(1:length(R_dat))
R_dat = R_dat ./ [(-exp(-value(m) * (x[i]-1)) + exp(-value(m) * (T+1-x[i]))) / sqrt(exp(-value(m) * (x[i]-1)) + exp(-value(m) * (T+1-x[i]))) for i in 1:length(x)]
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
v = value(R_dat[40])
e = err(R_dat[40])
figure()
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-20*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
end
return f, syst, f_i, weight, pval
end
function get_f_tm(corr_pp::juobs.Corr, m::uwreal, ens::EnsInfo, PS::String;
tm::Union{Vector{Vector{Int64}}, Nothing}=nothing, tM::Union{Vector{Vector{Int64}}, Nothing}=nothing,
pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
T = length(corr_pp.obs)
y0 = corr_pp.y0
mu = corr_pp.mu
pp_dat = corr_pp.obs
aux = exp.((collect(0:T-1) .- (T-1)/2 ) .* [m for k in 1:T])
R_dat = pp_dat .* aux ./ [((pp_dat[T-y0])^2)^(1/4) for k = 1:T]
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-1-y0-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(R_dat[Int64(round(T / 2))])
R, syst, R_i, weight, pval = model_av([fit_exp, fit_const], R_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
f = sqrt(2) * (mu[1] + mu[2]) * R / m^1.5
if pl == true
isnothing(wpm) ? uwerr(R) : uwerr(R, wpm)
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
isnothing(wpm) ? uwerr.(R_i) : [uwerr(R_i[i], wpm) for i in 1:length(R_i)]
v = value(R)
e = err(R)
fig = figure("eff")
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-20*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
fig = figure("model av")
subplot(411)
fill_between(1:length(R_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+10*e)
subplot(412)
ylabel(L"$R_i$")
fill_between(1:length(R_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_i), value.(R_i), err.(R_i), fmt="x", color="black")
subplot(413)
ylabel(L"$W_i$")
bar(1:length(R_i), weight, color="green")
subplot(414)
ylabel("p-value")
bar(1:length(R_i), pval, color="green")
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,".pdf"))
#close("all")
end
return f, syst, R_i, weight, pval
end
function get_f_tm(corr_ppL::juobs.Corr, corr_ppR::juobs.Corr, m::uwreal, ens::EnsInfo, PS::String;
tm::Union{Vector{Vector{Int64}}, Nothing}=nothing, tM::Union{Vector{Vector{Int64}}, Nothing}=nothing,
pl::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
T = length(corr_ppL.obs)
y0 = corr_ppL.y0
mu = corr_ppL.mu
ppL_dat = corr_ppL.obs
ppR_dat = corr_ppR.obs
f1 = [ppL_dat[T - y0] for k = 1:T]
R_dat = ((ppL_dat .* ppR_dat ./ f1).^2).^(1/4)
isnothing(tm) ? tm = [[y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40], [i for i in Int(round(T / 3)):Int(round(T / 3))+10]] : tm=tm
isnothing(tM) ? tM = [[T-10,T-15,T-20,T-25,T-30,T-35,T-40], [i for i in Int(round(2 * T / 3)):Int(round(2 * T / 3))+10]] : tM=tM
@.fit_exp(x,p) = p[1] + p[2] * exp(-p[3] * (x-y0)) + p[4] * exp(-p[5] * (T-1-y0-x))
@.fit_const(x,p) = p[1] * x ^ 0
k1 = 5
k2 = 1
guess = value(R_dat[Int64(round(T / 2))])
R, syst, R_i, weight, pval = model_av([fit_exp, fit_const], R_dat, guess, tm=tm, tM=tM, k=[k1,k2], wpm=wpm)
f = sqrt(2) * (mu[1] + mu[2]) * R / m^1.5
if pl == true
isnothing(wpm) ? uwerr(R) : uwerr(R, wpm)
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
isnothing(wpm) ? uwerr.(R_i) : [uwerr(R_i[i], wpm) for i in 1:length(R_i)]
v = value(R)
e = err(R)
fig = figure("eff")
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-20*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
fig = figure("model av")
subplot(411)
fill_between(1:length(R_dat), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
ylim(v-10*e, v+10*e)
subplot(412)
ylabel(L"$R_i$")
fill_between(1:length(R_i), v-e, v+e, color="green", alpha=0.5)
errorbar(1:length(R_i), value.(R_i), err.(R_i), fmt="x", color="black")
subplot(413)
ylabel(L"$W_i$")
bar(1:length(R_i), weight, color="green")
subplot(414)
ylabel("p-value")
bar(1:length(R_i), pval, color="green")
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,".pdf"))
#close("all")
end
return f, syst, R_i, weight, pval
end
function get_f_tm_pbc(corr_pp::juobs.Corr, m::uwreal, ens::EnsInfo, PS::String;
tm::Union{Vector{Int64}, Nothing}=nothing, tM::Union{Vector{Int64}, Nothing}=nothing,
pl::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
pp_dat = corr_pp.obs[2:end]
mu = corr_pp.mu
if ens.id == "E250"
global T = 192
global Thalf = 96
elseif ens.id == "D450"
global T = 128
global Thalf = 64
end
isnothing(tm) ? tm = [y0+10,y0+15,y0+20,y0+25,y0+30,y0+35,y0+40] : tm=tm
isnothing(tM) ? tM = [Thalf] : tM=tM
@.fit_fun_pp(x,p) = p[1] * (exp(-value(m) * (x-1)) + exp(-value(m) * (T+1-x)))
k = 1
cpp, syst, cpp_i, weight, pval = model_av(fit_fun_pp, pp_dat, .5, tm=tm, tM=tM, k=k, wpm=wpm)
f_i = [sqrt(2) * (mu[1] + mu[2]) / m ^ 1.5 for i in 1:length(cpp_i)] .* sqrt.(cpp_i)
f = sqrt(2) * (mu[1] + mu[2]) / m ^ 1.5 * sqrt(cpp)
if pl == true
##TODO
bla = 1
x = collect(1:length(pp_dat))
R_dat = pp_dat #./ [exp(-m * (x[i]-1)) + exp(-m * (T+1-x[i])) for i in 1:length(x)]
vec = cpp .* [exp(-m * (x[i]-1)) + exp(-m * (T+1-x[i])) for i in 1:length(x)]
isnothing(wpm) ? uwerr(cpp) : uwerr(cpp, wpm)
isnothing(wpm) ? uwerr.(R_dat) : [uwerr(R_dat[i], wpm) for i in 1:length(R_dat)]
isnothing(wpm) ? uwerr.(vec) : [uwerr(vec[i], wpm) for i in 1:length(vec)]
v = value.(vec)
e = err.(vec)
#v = value(cpp)
#e = err(cpp)
figure()
errorbar(1:length(R_dat), value.(R_dat), err.(R_dat), fmt="x", color="black")
fill_between(x, v-e, v+e, color="gray", alpha=0.75)
ylabel(L"$R_\mathrm{PS}$")
xlabel(L"$x_0/a$")
#ylim(1e-8, 1e-7)
#xlim(0,length(R_dat))
yscale("log")
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.pdf"))
#@gp x value.(R_dat) err.(R_dat) "w errorbars t ''"
#@gp:- x v-e v+e "w filledcurves t ''"
#@gp:- "set logscale y"
#save(string("/home/asaez/cls_ens/codes/lattA.jl/plots/R_",ens.id,"_",PS,"_plat.gp"))
end
return f, syst, f_i, weight, pval
end
function get_w0t0_plat(path::String, ens::EnsInfo, plat1::Vector{Int64}, plat2::Vector{Int64};
pl::Bool=false,
rw=false, npol::Int64=2, ws::ADerrors.wspace=ADerrors.wsg,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
w0_guess::Union{Float64, Nothing}=nothing,
tau::Union{Float64, Nothing}=nothing)
y0 = 1 ## assumes this is the case, hardcoded, some ensembles will not fulfil !
println("WARNING!: make sure t_src is 1 in this ensemble")
t2YM, tdt2YM, W_obs, t, t_aux = get_YM_dYM(path, ens, rw=rw, w0=w0_guess, tau=tau)
dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1) / 2)
model(x, p) = get_model(x, p, npol)
fmin(x, p) = model(x, p) .- 0.3
function fit_defs(x,W) ## uncorrelated fit
chisq(p,d) = sum((d .- [p[1] for i in 1:length(x)]).^2 .* W)
return chisq
end
## w0
tdt2YM_guess = [plat_av(tdt2YM[:,i], plat1) for i in 1:length(tdt2YM[1,:])]
nt0 = findmin(abs.(value.(tdt2YM_guess[2:end-1]) .- 0.3))[2] + 1
x = t_aux[nt0-dt0:nt0+dt0]
xmax = size(tdt2YM, 1)
T = xmax - 1 - y0
tdt2E_i = Array{uwreal,1}()
syst_i = Array{uwreal,1}()
for j in 1:2*dt0+1
i = Int(round(dt0+0.5))
dat = tdt2YM[:,nt0-dt0-1+j]
if j == i
tdt2E_aux = plat_av(dat, plat1)
uwerr.(dat)
chisq = fit_defs(collect(plat1[1]:plat1[2]), 1 ./ err.(dat[plat1[1]:plat1[2]]) .^ 2)
chi2 = chisq(tdt2E_aux,dat[plat1[1]:plat1[2]])
isnothing(wpm) ? pval = pvalue(chisq,value(chi2),value.(tdt2E_aux),dat[plat1[1]:plat1[2]]) : pval = pvalue(chisq,value(chi2),value.(tdt2E_aux),dat[plat1[1]:plat1[2]],wpm=wpm)
if pl == true
isnothing(wpm) ? uwerr(tdt2E_aux) : uwerr(tdt2E_aux, wpm)
isnothing(wpm) ? uwerr.(dat) : [uwerr(dat[i], wpm) for i in 1:length(dat)]
v = value(tdt2E_aux)
e = err(tdt2E_aux)
fig = figure(string("pvalue=",pval))
errorbar(1:length(dat), value.(dat), err.(dat), fmt="x", color="black")
x_plot = collect(plat1[1]:plat1[2])
fill_between(x_plot, v-e, v+e, color="green", alpha=0.3)
ylabel(L"$t\frac{d}{dt}t^2E(x_0,t)$")
xlabel(L"$x_0$")
ylim(v-5*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/tdt2E_",ens.id,"_plat.pdf"))
end
else
tdt2E_aux = plat_av(dat, plat1)
end
push!(tdt2E_i, tdt2E_aux)
end
if tdt2E_i[end] < 0.30 || tdt2E_i[1] > 0.30
println("WARNING!: extrapolating w0/a")
end
#uwerr.(tdt2E_i)
#uwerr.(syst_i)
#println("tdt2E_i = ", tdt2E_i)
#println("syst_i = ", syst_i)
par, aux = fit_alg(model, x[1:length(tdt2E_i)], tdt2E_i, npol, wpm=wpm)
w0 = root_error(fmin, t_aux[nt0], par)
if pl == true
v = value.(tdt2E_i)
e = err.(tdt2E_i)
uwerr(w0)
fig = figure("tdt2E_vs_t")
errorbar(x[1:length(tdt2E_i)], v, e, fmt="x")
errorbar(value(w0), 0.3, xerr=err(w0), fmt="x")
ylabel(L"$t\frac{d}{dt}t^2\left<E\right>$")
xlabel(L"$t/a^2$")
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/w0_",ens.id,".pdf"))
#close("all")
end
##
## t0
t2YM_guess = [plat_av(t2YM[:,i], plat2) for i in 1:length(t2YM[1,:])]
nt0 = findmin(abs.(value.(t2YM_guess[2:end-1]) .- 0.3))[2] + 1
x = t[nt0-dt0:nt0+dt0]
xmax = size(t2YM, 1)
T = xmax - 1 - y0
t2E_i = Array{uwreal,1}()
syst_i = Array{uwreal,1}()
for j in 1:2*dt0+1
i = Int(round(dt0+0.5))
dat = t2YM[:,nt0-dt0-1+j]
if j == i
t2E_aux = plat_av(dat, plat2)
uwerr.(dat)
chisq = fit_defs(collect(plat2[1]:plat2[2]), 1 ./ err.(dat[plat2[1]:plat2[2]]) .^ 2)
chi2 = chisq(t2E_aux,dat[plat2[1]:plat2[2]])
isnothing(wpm) ? pval = pvalue(chisq,value(chi2),value.(t2E_aux),dat[plat2[1]:plat2[2]]) : pval = pvalue(chisq,value(chi2),value.(t2E_aux),dat[plat2[1]:plat2[2]],wpm=wpm)
if pl == true
isnothing(wpm) ? uwerr(t2E_aux) : uwerr(t2E_aux, wpm)
isnothing(wpm) ? uwerr.(dat) : [uwerr(dat[i], wpm) for i in 1:length(dat)]
v = value(t2E_aux)
e = err(t2E_aux)
fig = figure(string("pvalue= ",pval))
errorbar(1:length(dat), value.(dat), err.(dat), fmt="x", color="black")
x_plot = collect(plat2[1]:plat2[2])
fill_between(x_plot, v-e, v+e, color="green", alpha=0.3)
ylabel(L"$t^2E(x_0,t)$")
xlabel(L"$x_0$")
ylim(v-5*e, v+20*e)
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/t2E_",ens.id,"_plat.pdf"))
end
else
t2E_aux = plat_av(dat, plat2)
end
push!(t2E_i, t2E_aux)
end
if t2E_i[end] < 0.30 || t2E_i[1] > 0.30
println("WARNING!: extrapolating t0/a")
end
#uwerr.(t2E_i)
#uwerr.(syst_i)
#println("t2E_i = ", t2E_i)
#println("syst_i = ", syst_i)
par, aux = fit_alg(model, x[1:length(t2E_i)], t2E_i, npol, wpm=wpm)
t0 = root_error(fmin, t[nt0], par)
if pl == true
v = value.(t2E_i)
e = err.(t2E_i)
uwerr(t0)
fig = figure("t2E_vs_t")
errorbar(x[1:length(t2E_i)], v, e, fmt="x")
errorbar(value(t0), 0.3, xerr=err(t0), fmt="x")
ylabel(L"$t^2\left<E\right>$")
xlabel(L"$t/a^2$")
tight_layout()
#savefig(string("/home/asaez/cls_ens/codes/lattA.jl/plots/t0_",ens.id,".pdf"))
#close("all")
end
##
return w0,t0
end end
\ No newline at end of file
src/juobs_reader.jl
View file @ 619245f8
...@@ -658,7 +658,7 @@ dat2 = read_mesons([path3, path4], "G5", "G5") ...@@ -658,7 +658,7 @@ dat2 = read_mesons([path3, path4], "G5", "G5")
concat_data!(dat, dat2) concat_data!(dat, dat2)
``` ```
""" """
function concat_data!(data1::Vector{CData}, data2::Vector{CData}) function concat_data!(data1::Vector{juobs.CData}, data2::Vector{juobs.CData})
N = length(data1) N = length(data1)
if length(data1) != length(data2) if length(data1) != length(data2)
error("number of correlators do not match") error("number of correlators do not match")
...@@ -667,10 +667,15 @@ function concat_data!(data1::Vector{CData}, data2::Vector{CData}) ...@@ -667,10 +667,15 @@ function concat_data!(data1::Vector{CData}, data2::Vector{CData})
data1[k].vcfg = vcat(data1[k].vcfg, data2[k].vcfg) 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].re_data = vcat(data1[k].re_data, data2[k].re_data)
data1[k].im_data = vcat(data1[k].im_data, data2[k].im_data) data1[k].im_data = vcat(data1[k].im_data, data2[k].im_data)
idx = sortperm(data1[k].vcfg)
data1[k].vcfg = data1[k].vcfg[idx]
data1[k].re_data = data1[k].re_data[idx, :]
data1[k].im_data = data1[k].im_data[idx, :]
end end
return nothing return nothing
end end
function concat_data!(data1::Vector{Vector{CData}}, data2::Vector{Vector{CData}})
function concat_data!(data1::Vector{Vector{juobs.CData}}, data2::Vector{Vector{juobs.CData}})
N = length(data1) N = length(data1)
if length(data1) != length(data2) if length(data1) != length(data2)
error("number of correlators do not match") error("number of correlators do not match")
...@@ -688,7 +693,312 @@ function concat_data!(data1::Vector{Vector{CData}}, data2::Vector{Vector{CData}} ...@@ -688,7 +693,312 @@ function concat_data!(data1::Vector{Vector{CData}}, data2::Vector{Vector{CData}}
data1[k][r].vcfg = data1[k][r].vcfg[idx] data1[k][r].vcfg = data1[k][r].vcfg[idx]
data1[k][r].re_data = data1[k][r].re_data[idx, :] data1[k][r].re_data = data1[k][r].re_data[idx, :]
data1[k][r].im_data = data1[k][r].im_data[idx, :] data1[k][r].im_data = data1[k][r].im_data[idx, :]
end end
end end
return nothing return nothing
end end
function get_YM(path::String, ens::EnsInfo; rw=false, ws::ADerrors.wspace=ADerrors.wsg, w0::Union{Float64, Nothing}=nothing)
path_ms = joinpath(path, ens.id, "gf")
path_ms = filter(x->occursin(".dat", x), readdir(path_ms, join=true))
Y = read_ms.(path_ms, dtr=ens.dtr)
nr = length(Y)
Ysl = getfield.(Y, :obs)
t = getfield.(Y, :t)
t = t[1]
id = getfield.(Y, :id)
replica = size.(Ysl, 1)
L = ens.L
id = ens.id
#T = length(Y[:,1]) - y0
y0 = 1 ## assumes this is the case, hardcoded, some ensembles will not fulfil !
println("WARNING!: make sure t_src is 1 in this ensemble")
#Truncation
if id in keys(ADerrors.wsg.str2id)
n_ws = findfirst(x-> x == ws.str2id[id], ws.map_nob)
if !isnothing(n_ws)
ivrep_ws = ws.fluc[n_ws].ivrep
if length(replica) != length(ivrep_ws)
error("Different number of replicas")
end
for k = 1:length(replica)
if replica[k] > ivrep_ws[k]
println("Automatic truncation in Ysl ", ivrep_ws[k], " / ", replica[k], ". R = ", k)
Ysl[k] = Ysl[k][1:ivrep_ws[k], :, :]
elseif replica[k] < ivrep_ws[k]
error("Automatic truncation failed. R = ", replica[k], "\nTry using truncate_data!")
end
end
replica = size.(Ysl, 1)
end
end
tmp = Ysl[1]
[tmp = cat(tmp, Ysl[k], dims=1) for k = 2:nr]
xmax = size(tmp, 2)
T = xmax - 1 - y0
Y_aux = Matrix{uwreal}(undef, xmax, length(t))
if rw
path_rw = joinpath(path, ens.id, "rwf")
path_rw = filter(x->occursin(".dat", x), readdir(path_rw, 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_2[1],rwf_1[1])]
elseif ens.id in ["E250", "E300", "J500", "J501", "D450"]
rwf = [read_ms1(path_rw[i], v=ens.vrw[i]) for i in 1:length(ens.vrw)]
[Ysl[k] = Ysl[k][1:size(rwf[k],2), :, :] for k in 1:length(Ysl)]
else
rwf = read_ms1.(path_rw, v=ens.vrw)
end
Ysl_r, W = juobs.apply_rw(Ysl, rwf)
tmp_r = Ysl_r[1]
tmp_W = W[1]
[tmp_r = cat(tmp_r, Ysl_r[k], dims=1) for k = 2:nr]
[tmp_W = cat(tmp_W, W[k], dims=1) for k = 2:nr]
W_obs = uwreal(tmp_W, id, replica, collect(1:length(tmp_W)), sum(replica))
WY_aux = Matrix{uwreal}(undef, xmax, length(t))
end
for i = 1:xmax
k = 1
for j = 1:length(t)
if !rw
Y_aux[i, k] = uwreal(tmp[:, i, j], id, replica)
else
WY_aux[i, k] = uwreal(tmp_r[:, i, j], id, replica, collect(1:length(tmp_W)), sum(replica))
Y_aux[i, k] = WY_aux[i, k] / W_obs
end
k = k + 1
end
end
t2YM = similar(Y_aux)
tdt2YM = similar(Y_aux)
for i in 1:length(Y_aux[:,1])
t2YM[i,:] = Y_aux[i,:] .* t .^ 2 ./ L ^ 3
end
for i in 1:length(t2YM[:,1])
if isnothing(w0)
tdt2YM[i,2:end-1] = [(t2YM[i,j+1] - t2YM[i,j-1]) / (t[j+1] - t[j-1]) * t[j] for j in 2:length(t2YM[i,:])-1]
tdt2YM[i,1] = tdt2YM[i,end] = t2YM[i,1]
else
ixm = findmin(abs.(t .- (w0-0.5)))[2]
ixM = findmin(abs.(t[1:end-1] .- (w0+0.5)))[2]
tdt2YM[i,1:end] .= t2YM[i,2]
tdt2YM[i,ixm:ixM] = [(t2YM[i,j+1] - t2YM[i,j-1]) / (t[j+1] - t[j-1]) * t[j] for j in ixm:ixM]
end
end
return t2YM, tdt2YM, W_obs, t
end
function get_YM_dYM(path::String, ens::EnsInfo; rw=false, ws::ADerrors.wspace=ADerrors.wsg, w0::Union{Float64, Nothing}=nothing, tau::Union{Float64, Nothing}=nothing)
path_ms = joinpath(path, ens.id, "gf")
path_ms = filter(x->occursin(".dat", x), readdir(path_ms, join=true))
Y = read_ms.(path_ms, dtr=ens.dtr)
truncate_data!(Y, ens.cnfg)
nr = length(Y)
Ysl = getfield.(Y, :obs)
t = getfield.(Y, :t)
t = t[1]
id = getfield.(Y, :id)
replica = size.(Ysl, 1)
L = ens.L
id = ens.id
#T = length(Y[:,1]) - y0
y0 = 1 ## assumes this is the case, hardcoded, some ensembles will not fulfil !
println("WARNING!: make sure t_src is 1 in this ensemble")
#Truncation
if id in keys(ADerrors.wsg.str2id)
n_ws = findfirst(x-> x == ws.str2id[id], ws.map_nob)
if !isnothing(n_ws)
ivrep_ws = ws.fluc[n_ws].ivrep
if length(replica) != length(ivrep_ws)
error("Different number of replicas")
end
for k = 1:length(replica)
if replica[k] > ivrep_ws[k]
println("Automatic truncation in Ysl ", ivrep_ws[k], " / ", replica[k], ". R = ", k)
Ysl[k] = Ysl[k][1:ivrep_ws[k], :, :]
elseif replica[k] < ivrep_ws[k]
error("Automatic truncation failed. R = ", replica[k], "\nTry using truncate_data!")
end
end
replica = size.(Ysl, 1)
end
end
tmp = Ysl[1]
[tmp = cat(tmp, Ysl[k], dims=1) for k = 2:nr]
xmax = size(tmp, 2)
T = xmax - 1 - y0
Y_aux = Matrix{uwreal}(undef, xmax, length(t))
if rw
path_rw = joinpath(path, ens.id, "rwf")
path_rw = filter(x->occursin(".dat", x), readdir(path_rw, join=true))
if ens.id in ["H105", "J500", "J501"]
rwf = [read_ms1(path_rw[i], v=ens.vrw[i]) for i in 1:length(ens.vrw)]
[Ysl[k] = Ysl[k][1:size(rwf[k],2), :, :] for k in 1:length(Ysl)]
else
rwf = read_ms1.(path_rw, v=ens.vrw)
end
Ysl_r, W = juobs.apply_rw(Ysl, rwf)
tmp_r = Ysl_r[1]
tmp_W = W[1]
[tmp_r = cat(tmp_r, Ysl_r[k], dims=1) for k = 2:nr]
[tmp_W = cat(tmp_W, W[k], dims=1) for k = 2:nr]
W_obs = uwreal(tmp_W, id, replica, collect(1:length(tmp_W)), sum(replica))
WY_aux = Matrix{uwreal}(undef, xmax, length(t))
end
for i = 1:xmax
k = 1
for j = 1:length(t)
if !rw
Y_aux[i, k] = uwreal(tmp[:, i, j], id, replica)
else
WY_aux[i, k] = uwreal(tmp_r[:, i, j], id, replica, collect(1:length(tmp_W)), sum(replica))
Y_aux[i, k] = WY_aux[i, k] / W_obs
end
k = k + 1
end
end
if tau == nothing
t2YM = similar(Y_aux)
tdt2YM = Matrix{uwreal}(undef,size(t2YM)[1],size(t2YM)[2]-2)
else
tau_ind = Int64(div(tau, t[2]-t[1]))
if tau >= 0.0
global t = t[1:end-tau_ind]
elseif tau < 0.0
global t = t[1-tau_ind:end]
end
t2YM = Matrix{uwreal}(undef,size(Y_aux)[1],length(t))
tdt2YM = Matrix{uwreal}(undef,size(t2YM)[1],size(t2YM)[2]-2)
end
for i in 1:length(Y_aux[:,1])
if tau == nothing
t2YM[i,:] = Y_aux[i,:] .* t .^ 2 ./ L ^ 3
else
if tau >= 0.0
t2YM[i,1:end] = Y_aux[i,1+tau_ind:end] .* t .^ 2 ./ L ^ 3
elseif tau < 0.0
t2YM[i,1:end] = Y_aux[i,1:end+tau_ind] .* t .^ 2 ./ L ^ 3
end
end
end
for i in 1:length(t2YM[:,1])
if isnothing(w0)
tdt2YM[i,1:end] = [(t2YM[i,j+1] - t2YM[i,j-1]) / (t[j+1] - t[j-1]) * t[j] for j in 2:length(t2YM[i,:])-1]
global t_aux = t[2:end-1]
else
global t_aux = t[2:end-1]
ixm = findmin(abs.(t_aux .- (w0-0.5)))[2]
ixM = findmin(abs.(t_aux[1:end-1] .- (w0+0.5)))[2]
tdt2YM[i,1:ixm-1] .= t2YM[i,2]
tdt2YM[i,ixm:end] = [(t2YM[i,j+1] - t2YM[i,j-1]) / (t[j+1] - t[j-1]) * t[j] for j in ixm+1:length(t2YM[i,:])-1]
end
end
return t2YM, tdt2YM, W_obs, t, t_aux
end
function read_mesons_multichunks(path::String, ens::String, g1::String="G5", g2::String="G5"; legacy::Bool=false)
idx_ts001 = findall(x->occursin("ts001",x), db[ens])
idx_tsT = findall(x->!occursin("ts001",x), db[ens])
store_cdata_aux = []
for i in db[ens]
aux = filter(x->occursin(i, x), readdir(path, join=true))
data_chunk = read_mesons(aux, g1, g2, legacy=legacy);
push!(store_cdata_aux, data_chunk)
end
#=
if ens == "E300"
for i in 1:length(idx_ts001)-1
concat_data!(store_cdata_aux[idx_ts001[1]][1:2:end], store_cdata_aux[idx_ts001[i+1]])
end
concat_data!(store_cdata_aux[idx_tsT[1]], store_cdata_aux[idx_ts001[1]][2:2:end])
for i in 1:length(idx_tsT)-1
concat_data!(store_cdata_aux[idx_tsT[1]], store_cdata_aux[idx_tsT[i+1]])
end
store_cdata_aux[idx_ts001[1]] = store_cdata_aux[idx_ts001[1]][1:2:end]
else
for i in 1:length(idx_ts001)-1
concat_data!(store_cdata_aux[idx_ts001[1]], store_cdata_aux[idx_ts001[i+1]])
end
for i in 1:length(idx_tsT)-1
concat_data!(store_cdata_aux[idx_tsT[1]], store_cdata_aux[idx_tsT[i+1]])
end
end
=#
for i in 1:length(idx_ts001)-1
concat_data!(store_cdata_aux[idx_ts001[1]], store_cdata_aux[idx_ts001[i+1]])
end
for i in 1:length(idx_tsT)-1
concat_data!(store_cdata_aux[idx_tsT[1]], store_cdata_aux[idx_tsT[i+1]])
end
dat_ts001 = store_cdata_aux[idx_ts001[1]]
dat_ts190 = store_cdata_aux[idx_tsT[1]]
dat = Vector{Vector{juobs.CData}}()
for i in 1:length(dat_ts001)
push!(dat, dat_ts001[i])
push!(dat, dat_ts190[i])
end
return dat
end
function read_mesons_correction_multichunks(path::String, ens::String, g1::String="G5", g2::String="G5"; legacy::Bool=false)
idx_ts001 = findall(x->occursin("ts001",x), db_c[ens])
idx_tsT = findall(x->!occursin("ts001",x), db_c[ens])
store_cdata_aux = []
for i in db_c[ens]
aux = filter(x->occursin(i, x), readdir(path, join=true))
data_chunk = read_mesons_correction(aux, g1, g2, legacy=legacy);
push!(store_cdata_aux, data_chunk)
end
for i in 1:length(idx_ts001)-1
concat_data!(store_cdata_aux[idx_ts001[1]], store_cdata_aux[idx_ts001[i+1]])
end
for i in 1:length(idx_tsT)-1
concat_data!(store_cdata_aux[idx_tsT[1]], store_cdata_aux[idx_tsT[i+1]])
end
dat_ts001 = store_cdata_aux[idx_ts001[1]]
if length(idx_tsT) >= 1
dat_ts190 = store_cdata_aux[idx_tsT[1]]
dat = Vector{Vector{juobs.CData}}()
for i in 1:length(dat_ts001)
push!(dat, dat_ts001[i])
push!(dat, dat_ts190[i])
end
else
dat = dat_ts001
end
return dat
end
\ No newline at end of file
src/juobs_tools.jl
View file @ 619245f8
...@@ -655,6 +655,74 @@ function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, wpm::Union{Dict{Int64 ...@@ -655,6 +655,74 @@ function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, wpm::Union{Dict{Int64
return av return av
end end
function model_av(fun::Function, y::Vector{uwreal}, guess::Float64;
tm::Vector{Int64}, tM::Vector{Int64}, k::Int64,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
pval = Array{Float64,1}()
p_1 = Array{uwreal,1}()
TIC = Array{Float64,1}()
isnothing(wpm) ? uwerr.(y) : [uwerr(y[i], wpm) for i in 1:length(y)]
for i in tm
for j in tM
if i < j
x = collect(i:j)
y_aux = y[i:j]
try
up, chi2, chi_exp, pval_i = fit_alg(fun,x,y_aux,k,guess,wpm=wpm)
push!(pval, pval_i)
push!(TIC, chi2 - 2*chi_exp)
push!(p_1, up[1])
catch e
end
end
end
end
TIC = TIC .- minimum(TIC)
weight = exp.(-0.5 .* TIC) ./ sum(exp.(-0.5 .* TIC))
p_av = sum(p_1 .* weight)
syst2 = sum(p_1 .^ 2 .* weight) - p_av ^ 2
return p_av, syst2, p_1, weight, pval
end
function model_av(fun::Vector{Function}, y::Vector{uwreal}, guess::Float64;
tm::Vector{Vector{Int64}}, tM::Vector{Vector{Int64}}, k::Vector{Int64},
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
pval = Array{Float64,1}()
p_1 = Array{uwreal,1}()
TIC = Array{Float64,1}()
isnothing(wpm) ? uwerr.(y) : [uwerr(y[i], wpm) for i in 1:length(y)]
for ind in 1:length(fun)
f = fun[ind]
for i in tm[ind]
for j in tM[ind]
if i < j
x = collect(i:j)
y_aux = y[i:j]
try
up, chi2, chi_exp, pval_i = fit_alg(f,x,y_aux,k[ind],guess,wpm=wpm)
push!(pval, pval_i)
push!(TIC, chi2 - 2*chi_exp)
push!(p_1, up[1])
catch e
end
end
end
end
end
TIC = TIC .- minimum(TIC)
weight = exp.(-0.5 .* TIC) ./ sum(exp.(-0.5 .* TIC))
p_av = sum(p_1 .* weight)
syst2 = sum(p_1 .^ 2 .* weight) - p_av ^ 2
return p_av, syst2, p_1, weight, pval
end
@doc raw""" @doc raw"""
bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Int64, pl::Bool, data::Bool; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Int64, pl::Bool, data::Bool; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
...@@ -1447,6 +1515,169 @@ function global_fit_routine(models::Vector{Function}, xdata::Vector{Array{Float6 ...@@ -1447,6 +1515,169 @@ function global_fit_routine(models::Vector{Function}, xdata::Vector{Array{Float6
return upar, value(chisq), chi2_exp, dof(fit) return upar, value(chisq), chi2_exp, dof(fit)
end end
function fit_alg(f::Function, x::Union{Vector{Int64}, Vector{Float64}, Matrix{Float64}}, y::Vector{uwreal},
n::Int64, guess::Union{Float64, Vector{Float64}, Nothing}=nothing;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
isnothing(wpm) ? uwerr.(y) : [uwerr(y[i], wpm) for i in 1:length(y)]
W = 1 ./ err.(y) .^ 2
chisq = fit_defs(f,x,W)
#lb = [-Inf for i in 1:n]
#ub = [+Inf for i in 1:n]
if guess == nothing
p0 = [.5 for i in 1:n]
else
p0 = [guess; [1. for i in 1:n-length(guess)]]
#lb[1] = .9 * guess
#ub[1] = 1.1 * guess
end
fit = curve_fit(f,x,value.(y),W,p0)#,lower=lb,upper=ub)
chi2 = sum(fit.resid .^ 2)
isnothing(wpm) ? (up,chi_exp) = fit_error(chisq,coef(fit),y) : (up,chi_exp) = fit_error(chisq,coef(fit),y,wpm)
isnothing(wpm) ? pval = pvalue(chisq,chi2,value.(up),y) : pval = pvalue(chisq,chi2,value.(up),y,wpm=wpm)
return up, chi2, chi_exp, pval
end
function fit_alg(model::Function,x::Array{Float64},y::Array{uwreal},param::Int64,W::Matrix{Float64};
guess::Union{Float64, Vector{Float64}, Nothing}=nothing)
if guess == nothing
p00 = [.5 for i in 1:param]
else
p00 = [guess; [1. for i in 1:param-length(guess)]]
#lb[1] = .9 * guess
#ub[1] = 1.1 * guess
end
chisq_corr(par,dat) = juobs.gen_chisq_correlated(model, x, W, par, dat)
fit = curve_fit(model,x,value.(y),W,p00)
up, chi_exp = fit_error(chisq_corr, coef(fit), y, W=W)
uwerr.(up)
chi2 = sum(fit.resid .^ 2)
pval = pvalue(chisq_corr, sum(fit.resid .^ 2), value.(up), y, W)
doff = dof(fit)
return up, chi2, chi_exp, pval, doff
end
function fit_alg(f::Vector{Function}, x::Vector{Matrix{Float64}}, y::Vector{Vector{uwreal}},
n::Int64, guess::Union{Float64, Vector{Float64}, Nothing}=nothing;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
y_n = y[1]
x_n = x[1]
for i in 2:length(x)
y_n = vcat(y_n, y[i])
x_n = vcat(x_n, x[i])
end
idx = [collect(1:size(x[i],1)) for i in 1:length(x)]
[idx[i] = idx[i] .+ idx[i-1][end] for i in 2:length(idx)]
isnothing(wpm) ? [uwerr.(y[i]) for i in 1:length(y)] : [[uwerr(y[i][j], wpm) for j in 1:length(y[i])] for i in 1:length(y)]
W = [1 ./ err.(y[i]) .^ 2 for i in 1:length(y)]
chisq = (par, y_n) -> sum([sum((y_n[idx[i]] .- f[i](x[i], par)) .^ 2 .* W[i]) for i in 1:length(x)])
min_fun(t) = chisq(t, value.(y_n))
if guess == nothing
p0 = [.5 for i in 1:n]
else
p0 = [guess; [1. for i in 1:n-length(guess)]]
end
sol = optimize(min_fun, p0, Optim.Options(g_tol=1e-8, iterations=1000000))
chi2 = min_fun(sol.minimizer)
isnothing(wpm) ? (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n) : (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n,wpm)
isnothing(wpm) ? pval = pvalue(chisq,chi2,value.(up),y_n) : pval = pvalue(chisq,chi2,value.(up),y_n,wpm=wpm)
return up, chi2, chi_exp, pval
end
##CHECK THIS ONE, NOT WORKING GOOD
function fit_alg(f::Vector{Function}, x::Vector{Matrix{Float64}}, y::Vector{Vector{uwreal}}, W::Vector{Matrix{Float64}},
n::Int64, guess::Union{Float64, Vector{Float64}, Nothing}=nothing;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
y_n = y[1]
x_n = x[1]
for i in 2:length(x)
y_n = vcat(y_n, y[i])
x_n = vcat(x_n, x[i])
end
idx = [collect(1:size(x[i],1)) for i in 1:length(x)]
[idx[i] = idx[i] .+ idx[i-1][end] for i in 2:length(idx)]
isnothing(wpm) ? [uwerr.(y[i]) for i in 1:length(y)] : [[uwerr(y[i][j], wpm) for j in 1:length(y[i])] for i in 1:length(y)]
chisq = (par, y_n) -> sum([juobs.gen_chisq_correlated(f[i], x[i], W[i], par, y_n[i]) for i in 1:length(x)])
min_fun(t) = chisq(t, value.(y_n))
if guess == nothing
p0 = [.5 for i in 1:n]
else
p0 = [guess; [1. for i in 1:n-length(guess)]]
end
sol = optimize(min_fun, p0, Optim.Options(g_tol=1e-8, iterations=1000000))
chi2 = min_fun(sol.minimizer)
W_aux = Matrix{Float64}(undef, length(y_n), length(y_n))
for i in 1:length(y_n)
for j in 1:length(y_n)
if i <= length(y[1])
if j <= length(y[1])
W_aux[i,j] = W[1][i,j]
else
W_aux[i,j] = .0
end
else
if j <= length(y[1])
W_aux[i,j] = .0
else
W_aux[i,j] = W[2][i-length(y[1]),j-length(y[1])]
end
end
end
end
isnothing(wpm) ? (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n,W=W_aux) : (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n,wpm,W=W_aux)
isnothing(wpm) ? pval = pvalue(chisq,chi2,value.(up),y_n,W_aux) : pval = pvalue(chisq,chi2,value.(up),y_n,wpm=wpm,W_aux)
return up, chi2, chi_exp, pval
end
function fit_alg(f::Vector{Function}, x::Vector{Matrix{Float64}}, y::Vector{Vector{uwreal}},
n::Int64, lb::Vector{Float64}, ub::Vector{Float64}, guess::Union{Float64, Vector{Float64}, Nothing}=nothing;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
y_n = y[1]
x_n = x[1]
for i in 2:length(x)
y_n = vcat(y_n, y[i])
x_n = vcat(x_n, x[i])
end
idx = [collect(1:size(x[i],1)) for i in 1:length(x)]
[idx[i] = idx[i] .+ idx[i-1][end] for i in 2:length(idx)]
isnothing(wpm) ? [uwerr.(y[i]) for i in 1:length(y)] : [[uwerr(y[i][j], wpm) for j in 1:length(y[i])] for i in 1:length(y)]
W = [1 ./ err.(y[i]) .^ 2 for i in 1:length(y)]
chisq = (par, y_n) -> sum([sum((y_n[idx[i]] .- f[i](x[i], par)) .^ 2 .* W[i]) for i in 1:length(x)])
min_fun(t) = chisq(t, value.(y_n))
if guess == nothing
p0 = [.5 for i in 1:n]
else
p0 = [guess; [1. for i in 1:n-length(guess)]]
end
sol = optimize(min_fun, lb, ub, p0, Fminbox(NelderMead()), Optim.Options(g_tol=1e-8, iterations=1000000))
chi2 = min_fun(sol.minimizer)
isnothing(wpm) ? (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n) : (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n,wpm)
isnothing(wpm) ? pval = pvalue(chisq,chi2,value.(up),y_n) : pval = pvalue(chisq,chi2,value.(up),y_n,wpm=wpm)
return up, chi2, chi_exp, pval
end
function fit_defs(f::Function,x,W) ## uncorrelated fit
chisq(p,d) = sum((d .- f(x,p)).^2 .* W)
return chisq
end
function gen_chisq(f::Function, x::Array{<:Real}, err::Vector{Float64}) #constrained function gen_chisq(f::Function, x::Array{<:Real}, err::Vector{Float64}) #constrained
chisq(par, dat) = sum((dat .- f(x,par)).^2 ./err.^2) chisq(par, dat) = sum((dat .- f(x,par)).^2 ./err.^2)
return chisq return chisq
...@@ -1728,3 +1959,34 @@ function pvalue(chisq::Function, ...@@ -1728,3 +1959,34 @@ function pvalue(chisq::Function,
return Q return Q
end end
function fve(mpi::uwreal, mk::uwreal, fpi::uwreal, fk::uwreal, ens::EnsInfo)
mm = [6,12,8,6,24,24,0,12,30,24,24,8,24,48,0,6,48,36,24,24]
nn = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]
meta = sqrt(4/3 * mk ^ 2 - 1/3 * mpi ^ 2)
jipi = value(mpi) ^ 2 / (4*pi*fpi) ^ 2
jik = value(mk) ^ 2 / (4*pi*fpi) ^ 2
jieta = value(meta) ^ 2 / (4*pi*fpi) ^ 2
mpiL = value(mpi) * ens.L
mkL = value(mk) * ens.L
metaL = value(meta) * ens.L
lampi = mpiL * sqrt.(nn)
lamk = mkL * sqrt.(nn)
lameta = metaL * sqrt.(nn)
g1pi = sum(4 .* mm ./ lampi .* besselk.(1, lampi))
g1k = sum(4 .* mm ./ lamk .* besselk.(1, lamk))
g1eta = sum(4 .* mm ./ lameta .* besselk.(1, lameta))
fve_mpi = 1/4 * jipi * g1pi - 1/12 * jieta * g1eta
fve_mk = 1/6 * jieta * g1eta
fve_fpi = - jipi * g1pi - 1/2 * jik * g1k
fve_fk = -3/8 * jipi *g1pi - 3/4 * jik * g1k - 3/8 * jieta * g1eta
mpi_infty = mpi / (1+fve_mpi)
mk_infty = mk / (1+fve_mk)
fpi_infty = fpi / (1+fve_fpi)
fk_infty = fk / (1+fve_fk)
return mpi_infty, mk_infty, fpi_infty, fk_infty
end
src/juobs_types.jl
View file @ 619245f8
...@@ -193,6 +193,42 @@ mutable struct YData ...@@ -193,6 +193,42 @@ mutable struct YData
YData(a, b, c, d) = new(a, b, c, d) YData(a, b, c, d) = new(a, b, c, d)
end end
@doc raw"""
EnsInfo(ens_id::String, info::Vector{Any})
Examples:
```@example
id = "H101"
ens = EnsInfo(id, ens_db[id])
```
"""
mutable struct EnsInfo
id::String
L::Int64
T::Int64
beta::Float64
ca::Float64
dtr::Int64
vrw::Union{String, Vector{String}}
cnfg::Vector{Int64}
function EnsInfo(ens_id::String, info::Vector{Any})
id = ens_id
L = info[1]
T = info[2]
beta = info[3]
dtr = info[4]
vrw = info[5]
cnfg = info[6]
p0 = 9.2056
p1 = -13.9847
g2 = 6 ./ beta
ca = - 0.006033 .* g2 .*( 1 .+exp.(p0 .+ p1./g2))
return new(id, L, T, beta, ca, dtr, vrw, cnfg)
end
end
function Base.show(io::IO, a::GHeader) function Base.show(io::IO, a::GHeader)
print(io, "ncorr = ", a.ncorr, "\t") print(io, "ncorr = ", a.ncorr, "\t")
print(io, "nnoise = ", a.nnoise, "\t") print(io, "nnoise = ", a.nnoise, "\t")
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment