Skip to content

J
juobs
Commit d6bd07c8 authored by Alessandro 's avatar Alessandro
Browse files
Options

corr_obs_TSM now supports different statistics between sloppy and correction... 

corr_obs_TSM now supports different statistics between sloppy and correction without neglecting the correlation. Extra flags idm_corr, idm_sl, nms have been added to allow this feature.
Documentation added
parent 85034d61
Hide whitespace changes
Inline Side-by-side
Showing with 367 additions and 116 deletions
constants/juobs_const.jl
View file @ d6bd07c8
...@@ -23,14 +23,22 @@ const t0_data = [2.86, 3.659, 5.164, 8.595] ...@@ -23,14 +23,22 @@ const t0_data = [2.86, 3.659, 5.164, 8.595]
const t0_error = [11, 16, 18, 29] .* 1e-3 const t0_error = [11, 16, 18, 29] .* 1e-3
const t0_ph_value = [0.413] const t0_ph_value = [0.413]
const t0_ph_error = ones(1,1) .* 5e-3 const t0_ph_error = ones(1,1) .* 5e-3
#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
## taking into account correlations in zm_tm ## 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]] 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") 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 ZP(b) = z[1] + z[2] * (b - 3.79) + z[3] * (b - 3.79)^2
Mrat = uwreal([.9148, 0.0088], "M/mhad") Mrat = uwreal([.9148, 0.0088], "M/mhad")
zm_tm_covar(b) = Mrat / ZP(b) zm_tm_covar(b) = Mrat / ZP(b)
## 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 #Cov matrices
const C1 = zeros(5, 5) const C1 = zeros(5, 5)
...@@ -40,7 +48,7 @@ const C4 = zeros(7, 7) ...@@ -40,7 +48,7 @@ const C4 = zeros(7, 7)
const C5 = zeros(5, 5) const C5 = zeros(5, 5)
const C6 = zeros(5, 5) const C6 = zeros(5, 5)
const C7 = zeros(5, 5) const C7 = zeros(5, 5)
const C8 = zeros(5, 5)
for i = 1:7 for i = 1:7
C4[i, i] = M_error[i] ^ 2 C4[i, i] = M_error[i] ^ 2
if i<=5 if i<=5
...@@ -48,7 +56,8 @@ for i = 1:7 ...@@ -48,7 +56,8 @@ for i = 1:7
C2[i, i] = ZM_tm_error[i] ^ 2 C2[i, i] = ZM_tm_error[i] ^ 2
C5[i, i] = ZA_err[i] ^ 2 C5[i, i] = ZA_err[i] ^ 2
C6[i, i] = ZS_over_ZP_err[i] ^ 2 C6[i, i] = ZS_over_ZP_err[i] ^ 2
C7[i, i] = ZV_err[i] ^ 2 C7[i, i] = ZV_err[i] ^ 2
C8[i, i] = RM_err[i] ^2
if i <= 4 if i <= 4
C3[i, i] = t0_error[i] ^ 2 C3[i, i] = t0_error[i] ^ 2
end end
...@@ -64,11 +73,13 @@ const M = cobs(M_values, C4, "charmed meson masses") ...@@ -64,11 +73,13 @@ const M = cobs(M_values, C4, "charmed meson masses")
const Za = cobs(ZA_data, C5, "ZA") const Za = cobs(ZA_data, C5, "ZA")
const Zv = cobs(ZV_data, C7, "ZV") 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 over ZP")
const RM = cobs(RM_data, C8, "rm")
zm(beta::Float64) = ZM[b_values .== beta][1] zm(beta::Float64) = ZM[b_values .== beta][1]
zm_tm(beta::Float64) = ZM_tm[b_values .== beta][1] zm_tm(beta::Float64) = ZM_tm[b_values .== beta][1]
t0(beta::Float64) = t0_[b_values2 .== beta][1] t0(beta::Float64) = t0_[b_values2 .== beta][1]
a(beta::Float64) = a_[b_values2 .== beta][1] a(beta::Float64) = a_[b_values2 .== beta][1]
za(beta::Float64) = Za[b_values .== beta][1] za(beta::Float64) = Za[b_values .== beta][1]
zv(beta::Float64) = Zv[b_values .== beta][1] 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]
\ No newline at end of file rm(beta::Float64) = RM[b_values .== beta][1]
\ No newline at end of file
src/juobs.jl
View file @ d6bd07c8
...@@ -12,6 +12,6 @@ include("juobs_obs.jl") ...@@ -12,6 +12,6 @@ include("juobs_obs.jl")
export read_mesons, read_mesons_correction, read_ms1, read_ms, read_md, truncate_data! export read_mesons, read_mesons_correction, read_ms1, read_ms, read_md, truncate_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 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 meff, mpcac, dec_const, dec_const_w, dec_const_pcvc, comp_t0
end # module end # module
src/juobs_linalg.jl
View file @ d6bd07c8
...@@ -180,7 +180,7 @@ matrices = [corr_diag, corr_upper] ...@@ -180,7 +180,7 @@ matrices = [corr_diag, corr_upper]
Julia> 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) n = length(a)
res = Vector{Vector{uwreal}}(undef, n) res = Vector{Vector{uwreal}}(undef, n)
[res[i] = uweigvals(a[i], a[t0]) for i=1:n] [res[i] = uweigvals(a[i], a[t0]) for i=1:n]
...@@ -237,7 +237,7 @@ function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 5, diag::Bool=t ...@@ -237,7 +237,7 @@ function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 5, diag::Bool=t
end end
function uweigvals_dim_geq3(a::Matrix{uwreal}; iter = 5, diag::Bool=true) function uweigvals_dim_geq3(a::Matrix{uwreal}; iter = 5, diag::Bool=true)
n = size(a,1) n = size(a,1)
a = tridiag_reduction(a) #a = tridiag_reduction(a)
for k in 1:iter for k in 1:iter
q_k, r_k = qr(a) q_k, r_k = qr(a)
a = uwdot(r_k, q_k) a = uwdot(r_k, q_k)
...@@ -315,20 +315,21 @@ res = uweigvecs(a) ##eigenvectors in column ...@@ -315,20 +315,21 @@ res = uweigvecs(a) ##eigenvectors in column
res1 = uweigvecs(a,b) ## generalised 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) n = size(a,1)
if n > 2 # if n > 2
a, q = tridiag_reduction(a, q_mat=true) # a, q = tridiag_reduction(a, q_mat=true)
end # end
u = idty(n) u = idty(n)
for k in 1:iter for k in 1:iter
q_k, r_k = qr(a) q_k, r_k = qr(a)
a = uwdot(r_k, q_k) a = uwdot(r_k, q_k)
u = uwdot(u, q_k) u = uwdot(u, q_k)
end end
n > 2 ? (return uwdot(q,u)) : (return u) #n > 2 ? (return uwdot(q,u)) : (return u)
return u
end 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) c = uwdot(invert(b), a)
#return uweigvecs(c, iter=iter) #return uweigvecs(c, iter=iter)
n = size(c,1) n = size(c,1)
...@@ -741,25 +742,37 @@ function idty(n::Int64) ...@@ -741,25 +742,37 @@ function idty(n::Int64)
return res return res
end end
function make_positive_def(A::Matrix) function make_positive_def(A::Matrix; eps::Float64=10^(-14))
B = 0.5 * (A + A') vals, vecs = eigen(Symmetric(A))
U, S, V = svd(B) idx = findall(x-> x<0.0, vals)
H = V * Diagonal(S) * V' vals[idx] .= eps
A_hat = 0.5 * (B + H) return Symmetric(vecs * Diagonal(vals) * vecs')
A_hat = 0.5 * ( A_hat + A_hat') #B = 0.5 * (A + A')
k = 0 #U, S, V = svd(B)
while !isposdef(A_hat) #H = V * Diagonal(S) * V'
mineig = minimum(eigvals(A_hat)) #A_hat = 0.5 * (B + H)
eps = 2.220446049250313e-16 #A_hat = 0.5 * ( A_hat + A_hat')
A_hat = A_hat + (-mineig*k^2 .+ eps*Matrix(I, size(A))) #k = 0
k+=1 #while !isposdef(A_hat)
end # mineig = minimum(eigvals(A_hat))
return 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 end
function invert_covar_matrix(A::Matrix) function invert_covar_matrix(A::Matrix; eps::Float64=10^(-9))
F = svd(A) 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 return cov_inv
end end
......
src/juobs_obs.jl
View file @ d6bd07c8
...@@ -23,9 +23,9 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo ...@@ -23,9 +23,9 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo
if pl if pl
isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm) isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm)
x = 1:length(aux) x = 1:length(aux)
y = value.(aux) y = abs.(value.(aux))
dy = err.(aux) dy = err.(aux)
v = value(mass) v = abs(value(mass))
e = err(mass) e = err(mass)
figure() figure()
...@@ -33,7 +33,9 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo ...@@ -33,7 +33,9 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo
errorbar(x, y, dy, fmt="x", color="black") errorbar(x, y, dy, fmt="x", color="black")
ylabel(L"$m_\mathrm{eff}$") ylabel(L"$m_\mathrm{eff}$")
xlabel(L"$x_0$") xlabel(L"$x_0$")
ylim(v-100*e, v+100*e) ylim(v-30*e, v+30*e)
# ylim(0.05, 0.075)
#xlim(64,length(y))
if !isnothing(kappa) if !isnothing(kappa)
title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2])) title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
...@@ -92,10 +94,12 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca: ...@@ -92,10 +94,12 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca:
corr_a0p = -a0p[1:end] corr_a0p = -a0p[1:end]
corr_pp = pp[1:end] corr_pp = pp[1:end]
der_a0p = (corr_a0p[3:end] .- corr_a0p[1:end-2]) / 2
if ca != 0.0 if ca != 0.0
der_a0p = (corr_a0p[3:end] .- corr_a0p[1:end-2]) / 2
der2_pp = (corr_pp[1:end-4] - 2*corr_pp[3:end-2] + corr_pp[5:end])/4 der2_pp = (corr_pp[1:end-4] - 2*corr_pp[3:end-2] + corr_pp[5:end])/4
der_a0p = der_a0p[2:end-1] + ca * der2_pp der_a0p = der_a0p[2:end-1] + ca * der2_pp
else
der_a0p = (corr_a0p[4:end-1] .- corr_a0p[2:end-3]) / 2
end end
#der_a0p = (corr_a0p[3:end] .- corr_a0p[1:end-2]) / 2 #der_a0p = (corr_a0p[3:end] .- corr_a0p[1:end-2]) / 2
...@@ -132,7 +136,7 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca: ...@@ -132,7 +136,7 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca:
if !data if !data
return mass return mass
else else
return (mass, aux) return (mass, aux)
end end
end end
...@@ -199,16 +203,17 @@ f_ratio = dec_const(corr_a0pL, corr_a0pR, corr_ppL, corr_ppR, [50, 60], m, pl=tr ...@@ -199,16 +203,17 @@ f_ratio = dec_const(corr_a0pL, corr_a0pR, corr_ppL, corr_ppR, [50, 60], m, pl=tr
function dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, function dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false,
kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
corr_a0p = -a0p corr_a0p = -a0p[2:end-1]
corr_pp = pp corr_pp = pp[2:end-1]
T = length(corr_a0p)
if ca != 0.0 if ca != 0.0
der_pp = (corr_pp[3:end] - corr_pp[1:end-2]) / 2 der_pp = (corr_pp[3:end] - corr_pp[1:end-2]) / 2
corr_a0p = corr_a0p[2:end-1] + ca * der_pp corr_a0p = corr_a0p[2:end-1] + ca * der_pp
aux = exp.((collect(1:T-2) .- (T-1)/2) .* [m for k = 1:T-2]) T = length(corr_a0p)
aux = exp.((collect(1:T) .- (T)/2) .* [m for k = 1:T])
else else
aux = exp.((collect(0:T-1) .- T/2) .* [m for k = 1:T]) T = length(corr_a0p)
aux = exp.((collect(1:T) .- T/2) .* [m for k = 1:T])
end end
R = corr_a0p .* aux ./ [((corr_pp[T-y0])^2)^(1/4) for k = 1:length(corr_a0p)] R = corr_a0p .* aux ./ [((corr_pp[T-y0])^2)^(1/4) for k = 1:length(corr_a0p)]
...@@ -369,7 +374,68 @@ end ...@@ -369,7 +374,68 @@ end
function dec_const(vv::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false, function dec_const(vv::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
return dec_const(vv.obs, plat, m, vv.y0+1, kappa=vv.kappa, pl=pl, data=data, wpm=wpm) return dec_const(vv.obs, plat, m, vv.y0, kappa=vv.kappa, pl=pl, data=data, wpm=wpm)
end
# this can be used also for pseudoscalar decay constants with wilson fermions
function dec_const_w(vv::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false,
kappa::Union{Vector{Float64}, Nothing}=nothing, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
corr_vv = vv[2:end-1]
corr_pp = pp[2:end-1]
T = length(corr_vv)
if ca != 0.0
der_pp = (corr_pp[3:end] - corr_pp[1:end-2]) / 2
corr_vv = corr_vv[2:end-1] + ca* der_pp
aux = exp.((collect(1:T-2) .- y0 ) .* fill(m, T-2))
else
#corr_vv = corr_vv[2:end-1]
aux = exp.((collect(1:T) .- y0 ) .* fill(m, T))
end
R = (aux .* corr_vv)
R_av = plat_av(R, plat, wpm)
f = sqrt(2 / m) * sqrt(R_av)
if pl
R .*= sqrt.(2 ./ [m for i in 1:length(R)])
uwerr.(R)
R_av *= sqrt(2/m)
isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
isnothing(wpm) ? uwerr(f) : uwerr(f, wpm)
x = 1:length(R)
y = value.(R)
dy = err.(R)
v = value(R_av)
e = err(R_av)
figure()
lbl = string(L"$af = $", sprint(show, f))
fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75, label=L"$R$")
errorbar(x, y, dy, fmt="x", color="black", label=lbl)
legend()
xlim(y0, y0 +plat[2] )
ylim(v-10*e, v+10*e)
ylabel(L"$af_{vec}$")
#ylabel(L"$R_\mathrm{av}$")
xlabel(L"$x_0$")
#title(L"$f_{B^*}$")
if !isnothing(kappa)
title(string( L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
end
display(gcf())
#t = "ps_decay_$(kappa[1])_$(kappa[2]).pdf"
#savefig(joinpath("/Users/ale/Il mio Drive/phd/secondment/3pf test/analysis/plots",t))
end
if !data
return f
else
return f, R
end
end
function dec_const_w(vv::Corr, pp::Corr, plat::Vector{Int64}, m::uwreal; ca::Float64=0.0, pl::Bool=true, data::Bool=false,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
return dec_const_w(vv.obs, pp.obs, plat, m, vv.y0, ca=ca, kappa=vv.mu, pl=pl, data=data, wpm=wpm)
end end
@doc raw""" @doc raw"""
...@@ -424,6 +490,8 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu ...@@ -424,6 +490,8 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu
errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black", label=lbl) errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black", label=lbl)
ylabel(L"$R$") ylabel(L"$R$")
xlabel(L"$x_0$") xlabel(L"$x_0$")
xlim(y0,y0+plat[2])
ylim(v-10*e, v+10*e)
legend() legend()
title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2])) title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
display(gcf()) display(gcf())
...@@ -612,6 +680,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, ...@@ -612,6 +680,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
fill_between(collect(plat[1]:plat[2]), v[t_pl]+e[t_pl], v[t_pl]-e[t_pl], alpha=0.75, color="green") 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)$") ylabel(L"$t^2E(x0, t)$")
xlabel(L"$x_0/a$") xlabel(L"$x_0/a$")
ylim(v[t_pl]-20*e[t_pl], v[t_pl]+20*e[t_pl])
title(string(L"$t/a^2 = $", t[nt0])) title(string(L"$t/a^2 = $", t[nt0]))
display(gcf()) display(gcf())
end end
......
src/juobs_tools.jl
View file @ d6bd07c8
...@@ -90,17 +90,19 @@ corr_pp = corr_obs.(data) ...@@ -90,17 +90,19 @@ corr_pp = corr_obs.(data)
corr_pp_r = corr_obs.(data, rw=[rw1, rw2]) corr_pp_r = corr_obs.(data, rw=[rw1, rw2])
``` ```
""" """
function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1, info::Bool=false) function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1, info::Bool=false, idm::Union{Vector{Int64},Nothing}=nothing, nms::Int64=Int64(maximum(cdata.vcfg)))
real ? data = cdata.re_data ./ L^3 : data = cdata.im_data ./ L^3 real ? data = cdata.re_data ./ L^3 : data = cdata.im_data ./ L^3
nt = size(data)[2] nt = size(data)[2]
if isnothing(rw) if isnothing(rw)
obs = [uwreal(data[:, x0], cdata.id) for x0 = 1:nt] idm = isnothing(idm) ? collect(1:nms) : idm
obs = [uwreal(data[:, x0], cdata.id, idm, nms) for x0 = 1:nt]
else else
idm = isnothing(idm) ? collect(1:nms) : idm
data_r, W = apply_rw(data, rw) data_r, W = apply_rw(data, rw)
ow = [uwreal(data_r[:, x0], cdata.id) for x0 = 1:nt] ow = [uwreal(data_r[:, x0], cdata.id, idm, nms) for x0 = 1:nt]
W_obs = uwreal(W, cdata.id) W_obs = uwreal(W, cdata.id, idm, nms)
obs = [ow[x0] / W_obs for x0 = 1:nt] obs = [ow[x0] / W_obs for x0 = 1:nt]
end end
...@@ -152,46 +154,103 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array ...@@ -152,46 +154,103 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
end end
end end
function corr_obs_TSM(cdata1::CData, cdata2::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1, info::Bool=false) @doc raw"""
if cdata1.id != cdata2.id 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)))
error("Error: cdata1 id != cdata2 id")
corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, info::Bool=false)
Creates a `Corr` struct for Truncated Solver Method (TSM) with the given `CData` structs `cdata_sl` and `cdata_ex` for a single replica.
Two arrays of `CData` can be passed as argument for multiple replicas.
The flag `real` select the real or imaginary part of the correlator.
If `rw` is specified, the method applies reweighting. `rw` is passed as a matrix of Float64 (`read_ms1`)
The correlator can be normalized with the volume factor if `L` is fixed.
The flag `info` provides extra output that contains information about the primary observables. The function returns the primary observables ``<WO>`` and ``<W>``
(it returns the observable <O> if rw=nothing)
The flags `idm_sl`, `idm_corr` and `nms` can be used if the observable is not measure in every configuration. In particular:
`idm_sl::Vector{Int64}` If given, label the configurations where the sloppy observables is measured. Not required if sloppy is full statistics.
`idm_corr::Vector{Int64}` If given, label the configurations where the correction observables is measured. Not required if correction is full statistics.
`nms::Int64` Total number of measurements in the ensembles. By default, nms=maximum(cdata_sl.vcfg) it is assumed that the sloppy observable is full statistics, while the exact observable is not.
Typically, only `idm_corr` has to be passed, if the exact observable is not full statistics.
By playing with these three flags you ensure that the Gamma method is still working even if the configurations between sloppy and correction do not match.
```@example
#Single replica
data_sloppy = read_mesons(path_sl, "G5", "G5")
data_correction = read_mesons_correction(path_corr, "G5", "G5")
rw = read_ms1(path_rw)
corr_pp = corr_obs_TSM.(data_sloppy, data_correction)
corr_pp_r = corr_obs_TSM.(data_sloppy, data_correction rw=rw)
#Single replica + Info
data_sloppy = read_mesons(path_sl, "G5", "G5")
data_correction = read_mesons_correction(path_corr, "G5", "G5")
rw = read_ms1(path_rw)
corr_pp, O = corr_obs_TSM(data_sloppy, data_correction, info=true)
corr_pp_r, WO, W = corr_obs(data_sloppy, data_correction,, rw=rw, info=true)
#Two replicas
data = read_mesons([path_r1, path_r2], "G5", "G5")
rw1 = read_ms1(path_rw1)
rw2 = read_ms1(path_rw2)
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)))
# cdata1 is sloppy, cdata2 is correction
if cdata_sl.id != cdata_corr.id[1:4]
error("Error: cdata_sl id != cdata_corr id")
end end
if cdata1.header != cdata2.header # Base.:(==) and Base.:(!=) are redifined in juobs_types.jl if cdata_sl.header != cdata_corr.header # Base.:(==) and Base.:(!=) are redifined in juobs_types.jl
error("Error: cdata1 header != cdata2 header") error("Error: cdata_sl header != cdata_corr header")
end end
data1 = real ? cdata1.re_data ./ L^3 : cdata1.im_data ./ L^3 data1 = real ? cdata_sl.re_data ./ L^3 : cdata_sl.im_data ./ L^3
data2 = real ? cdata2.re_data ./ L^3 : cdata2.im_data ./ L^3 data2 = real ? cdata_corr.re_data ./ L^3 : cdata_corr.im_data ./ L^3
nt = size(data1, 2)
id = cdata1.id
nt1 = size(data1, 2) # 1
nt2 = size(data2, 2)
id_sl = cdata_sl.id
id_corr = cdata_corr.id
# println(id_sl, " ", id_corr)
if isnothing(rw) if isnothing(rw)
obs1 = [uwreal(data1[:, x0], id) for x0 = 1:nt] idm_sl = isnothing(idm_sl) ? collect(1:nms) : idm_sl
obs2 = [uwreal(data2[:, x0], id) for x0 = 1:nt] 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]
else else
data1_r, W = apply_rw(data1, rw) idm_sl = isnothing(idm_sl) ? collect(1:nms) : idm_sl
data2_r, W = apply_rw(data2, rw) idm_corr = isnothing(idm_corr) ? collect(1:nms) : idm_corr
ow1 = [uwreal(data1_r[:, x0], id) for x0 = 1:nt] data1_r, W_sl = apply_rw(data1, rw)
ow2 = [uwreal(data2_r[:, x0], id) for x0 = 1:nt] 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]
W_obs = uwreal(W, id) W_obs_sl = uwreal(W_sl, id_sl, idm_sl, nms)
W_obs_corr = uwreal(W_corr, id_corr, idm_corr, nms)
obs1 = [ow1[x0] / W_obs for x0 = 1:nt] obs1 = [ow1[x0] / W_obs_sl for x0 = 1:nt1]
obs2 = [ow2[x0] / W_obs for x0 = 1:nt] obs2 = [ow2[x0] / W_obs_corr for x0 = 1:nt2]
end end
if info && !isnothing(rw) if info && !isnothing(rw)
return (Corr(obs1 + obs2, cdata1), ow1, ow2, W_obs) return (Corr(obs1 + obs2, cdata_sl), ow1, ow2, W_obs_sl)
elseif info && isnothing(rw) elseif info && isnothing(rw)
return (Corr(obs1 + obs2, cdata1), obs1, obs2) return (Corr(obs1 + obs2, cdata_sl), obs1, obs2)
else else
return Corr(obs1 + obs2, cdata1) return Corr(obs1 + obs2, cdata_sl)
end end
end end
function corr_obs_TSM(cdata1::Array{CData, 1}, cdata2::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, info::Bool=false) function corr_obs_TSM(cdata1::Array{CData, 1}, cdata2::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, info::Bool=false, idm::Union{Vector{Int64},Nothing}=nothing, nms::Union{Int64, Nothing}=nothing)
if any(getfield.(cdata1, :id) .!= getfield.(cdata2, :id)) if any(getfield.(cdata1, :id) .!= getfield.(cdata2, :id))
error("Error: cdata1 id != cdata2 id") error("Error: cdata1 id != cdata2 id")
end end
...@@ -237,11 +296,11 @@ function corr_obs_TSM(cdata1::Array{CData, 1}, cdata2::Array{CData, 1}; real::Bo ...@@ -237,11 +296,11 @@ function corr_obs_TSM(cdata1::Array{CData, 1}, cdata2::Array{CData, 1}; real::Bo
end end
if info && !isnothing(rw) if info && !isnothing(rw)
return (Corr(obs1 + obs2, cdata), ow1, ow2, W_obs) return (Corr(obs1 + obs2, cdata1), ow1, ow2, W_obs)
elseif info && isnothing(rw) elseif info && isnothing(rw)
return (Corr(obs1 + obs2, cdata), obs1, obs2) return (Corr(obs1 + obs2, cdata1), obs1, obs2)
else else
return Corr(obs1 + obs2, cdata) return Corr(obs1 + obs2, cdata1)
end end
end end
@doc raw""" @doc raw"""
...@@ -325,24 +384,24 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::uwreal, w::Union{uwr ...@@ -325,24 +384,24 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::uwreal, w::Union{uwr
p = findall(t-> t==1, a.prop) p = findall(t-> t==1, a.prop)
if nid != 1 if nid != 1
error("Error: neid > 1") @info("Error: neid > 1")
end end
id = ws.map_nob[p] id = ws.map_nob[p]
if !all(id .== id[1]) if !all(id .== id[1])
error("ids do not match") @info("ids do not match")
end end
id = ws.id2str[id[1]] id = ws.id2str[id[1]]
ivrep = getfield.(ws.fluc[p], :ivrep) ivrep = getfield.(ws.fluc[p], :ivrep)
ivrep1 = fill(ivrep[1], length(ivrep)) ivrep1 = fill(ivrep[1], length(ivrep))
if !all(ivrep .== ivrep1) if !all(ivrep .== ivrep1)
error("ivreps do not match") @info("ivreps do not match")
end end
ivrep = ivrep[1] ivrep = ivrep[1]
if length(md) != length(ivrep) if length(md) != length(ivrep)
error("Nr obs != Nr md") warning("Nr obs != Nr md")
end end
#md_aux as a Matrix + Automatic truncation #md_aux as a Matrix + Automatic truncation
...@@ -518,7 +577,7 @@ tmax_array = [80,81,82,83,84,85] ...@@ -518,7 +577,7 @@ tmax_array = [80,81,82,83,84,85]
(average, systematics) = bayesian_av([fun1,fun2],x,tmin_array,tmax_array,[k1,k2]) (average, systematics) = bayesian_av([fun1,fun2],x,tmin_array,tmax_array,[k1,k2])
``` ```
""" """
function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Int64, pl::Bool=false, data::Bool=false; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Int64; pl::Bool=false, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
weight_model = Array{Float64,1}() weight_model = Array{Float64,1}()
AIC = Array{Float64,1}() AIC = Array{Float64,1}()
...@@ -531,7 +590,7 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, ...@@ -531,7 +590,7 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
end end
total = length(y) total = length(y)
isnothing(wpm) ? uwerr.(y) : for i in 1:length(y) uwerr(y[i],wpm) end isnothing(wpm) ? uwerr.(y) : [uwerr(y[i],wpm) for i in 1:length(y) ]
for INDEX in tmin_array ## vary tmin for INDEX in tmin_array ## vary tmin
for j in tmax_array ## vary tmax for j in tmax_array ## vary tmax
...@@ -551,8 +610,8 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, ...@@ -551,8 +610,8 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
push!(AIC, chi2 + 2*k + 2*Ncut) push!(AIC, chi2 + 2*k + 2*Ncut)
push!(chi2chi2exp, chi2 / dof(fit)) push!(chi2chi2exp, chi2 / dof(fit))
push!(p1, up[1]) push!(p1, up[1])
push!(mods,string("[", INDEX+1, ",", j, "]")) push!(mods,string("[", INDEX+1, ",", j, "]"))
println("chi2_vs_exp " , chi2 / dof(fit), " chi2/dof ", sum(fit.resid.^2)/dof(fit))
catch e catch e
@warn string(":/ Negative window for error propagation at tmin = ", INDEX, ", tmax = ", j, "; skipping that point") @warn string(":/ Negative window for error propagation at tmin = ", INDEX, ", tmax = ", j, "; skipping that point")
end end
...@@ -566,7 +625,6 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, ...@@ -566,7 +625,6 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
p1_mean = sum(p1 .* weight_model)/sum(weight_model) ; isnothing(wpm) ? uwerr(p1_mean) : uwerr(p1_mean,wpm) p1_mean = sum(p1 .* weight_model)/sum(weight_model) ; isnothing(wpm) ? uwerr(p1_mean) : uwerr(p1_mean,wpm)
weight_model = weight_model ./ sum(weight_model) weight_model = weight_model ./ sum(weight_model)
systematic_err = sqrt(sum(p1 .^ 2 .* weight_model) - (sum(p1 .* weight_model)) ^ 2) systematic_err = sqrt(sum(p1 .^ 2 .* weight_model) - (sum(p1 .* weight_model)) ^ 2)
if pl if pl
x = 1:length(p1) x = 1:length(p1)
y = value.(p1) y = value.(p1)
...@@ -588,17 +646,18 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, ...@@ -588,17 +646,18 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
xlabel(L"model") xlabel(L"model")
display(gcf()) display(gcf())
close()
end end
if !data if !data
return (p1_mean, systematic_err) return (p1_mean, systematic_err)
else else
return (p1_mean, systematic_err, p1, weight_model) return (p1_mean, systematic_err, p1, weight_model, AIC .+ offset)
end end
end end
function bayesian_av(fun1::Function, fun2::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k1::Int64, k2::Int64, pl::Bool=false, data::Bool=false; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) function bayesian_av(fun1::Function, fun2::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k1::Int64, k2::Int64; pl::Bool=false, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
weight_model = Array{Float64,1}() weight_model = Array{Float64,1}()
AIC = Array{Float64,1}() AIC = Array{Float64,1}()
...@@ -849,24 +908,29 @@ fit_routine(model, xdata, ydata, param=3) ...@@ -849,24 +908,29 @@ fit_routine(model, xdata, ydata, param=3)
fit_routine(model, xdata, ydata, param=3, covar=true) fit_routine(model, xdata, ydata, param=3, covar=true)
``` ```
""" """
function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; info::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, correlated_fit::Bool=false) function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; info::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
inv_cov::Union{Matrix{Float64}, Nothing}=nothing)
isnothing(wpm) ? uwerr.(ydata) : [uwerr(yaux, wpm) for yaux in ydata] isnothing(wpm) ? uwerr.(ydata) : [uwerr(yaux, wpm) for yaux in ydata]
yval = value.(ydata) yval = value.(ydata)
yer = err.(ydata) yer = err.(ydata)
if !correlated_fit ## CODE HERE THE UPDATED VERSION OF CORRELATED FITS
if isnothing(inv_cov)
@info("Uncorrelated fit")
chisq = gen_chisq(model, xdata, yer) chisq = gen_chisq(model, xdata, yer)
fit = curve_fit(model, xdata, yval, 1.0 ./ yer.^2, fill(0.5, param)) fit = curve_fit(model, xdata, yval, 1.0 ./ yer.^2, fill(0.5, param))
(upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata) : fit_error(chisq, coef(fit), ydata, wpm) (upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata) : fit_error(chisq, coef(fit), ydata, wpm)
else else
covar_inv = make_positive_def(invert_covar_matrix(cov(ydata))) @info("Correlated fit")
#covar_inv = inv(get_covariance(ydata))
#covar_inv = (covar_inv + covar_inv')/2 #covar_inv = (covar_inv + covar_inv')/2
chisq(par,dat) = gen_chisq_correlated(model, xdata, covar_inv, par, dat) #covar_inv = inv(Symmetric(get_covariance(ydata)))
chisq(par, dat) = gen_chisq_correlated(model, xdata, covar_inv, par, dat)
fit = curve_fit(model, xdata, yval, covar_inv, fill(0.5, param)) fit = curve_fit(model, xdata, yval, covar_inv, fill(0.5, param))
println(chisq(coef(fit), ydata)) #println(chisq(coef(fit), ydata))
(upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata, W=covar_inv) : fit_error(chisq, coef(fit), ydata, wpm, W=covar_inv) (upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata, W=covar_inv) : fit_error(chisq, coef(fit), ydata, wpm, W=covar_inv)
end end
chi2_fit_res = sum(fit.resid.^2 ) chi2_fit_res = sum(fit.resid.^2 )
println("\n") println("\n")
...@@ -896,7 +960,9 @@ function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal} ...@@ -896,7 +960,9 @@ function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}
return upar, chi2_fit_res / chi_exp return upar, chi2_fit_res / chi_exp
end end
function fit_routine(model::Vector{Function}, xdata::Vector{Array{Float64, N}} where N, ydata::Vector{Array{uwreal, N}} where N, param::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) function fit_routine(model::Vector{Function}, xdata::Vector{Array{Float64, N}} where N, ydata::Vector{Array{uwreal, N}} where N, param::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
correlated_fit::Bool=false)
if !(length(model) == length(xdata) == length(ydata)) if !(length(model) == length(xdata) == length(ydata))
error("Dimension mismatch") error("Dimension mismatch")
end end
...@@ -921,20 +987,48 @@ function fit_routine(model::Vector{Function}, xdata::Vector{Array{Float64, N}} w ...@@ -921,20 +987,48 @@ function fit_routine(model::Vector{Function}, xdata::Vector{Array{Float64, N}} w
idx[i] = collect(j:stp) idx[i] = collect(j:stp)
j = stp + 1 j = stp + 1
end end
if !correlated_fit
chisq = (par, dat) -> sum([sum((dat[idx[i]] .- model[i](xdata[i], par)).^2 ./e[i].^2) for i=1:N])
else
inv_cov_tot = juobs.inv_covar(vcat(ydata...)) #[inv_covar_multi_id(dat[idx[ii]]) for ii in 1:N]
function chisq(par, dat)
res = Vector{uwreal}(undef,0)
for ii in 1:N
res = vcat(res, dat[idx[ii]] .- model[ii](xdata[ii],par))
end
chi = reshape(res, 1,:) * inv_cov_tot*res
return chi[1]
end
end
chisq = (par, dat) -> sum([sum((dat[idx[i]] .- model[i](xdata[i], par)).^2 ./e[i].^2) for i=1:N])
min_fun(t) = chisq(t, value.(dat)) min_fun(t) = chisq(t, value.(dat))
p = fill(0.5, param) p = fill(0.5, param)
#println(chisq(p,dat))
sol = optimize(min_fun, p, LBFGS()) sol = optimize(min_fun, p, LBFGS())
#println(chisq(sol.minimizer,dat))
(upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq, sol.minimizer, dat) : fit_error(chisq, sol.minimizer, dat, wpm) if !correlated_fit
@info("Uncorrelated fit ")
(upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq, sol.minimizer, dat) : fit_error(chisq, sol.minimizer, dat, wpm)
else
@info("Correlated fit ")
#inv_cov_merge = zeros(size(inv_cov_tot[1],1)*N, size(inv_cov_tot[1],1)*N )
#count = 1
#for k in 1:N
# dim_tmp = size(inv_cov_tot[k], 1)
# inv_cov_merge[count:dim_tmp+count-1, count:dim_tmp+count-1] = inv_cov_tot[k]
# count +=dim_tmp
#end
(upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq, sol.minimizer, dat, W=inv_cov_tot) : fit_error(chisq, sol.minimizer, dat, wpm, W=inv_cov_tot)
end
println("Chisq / chiexp: ", min_fun(sol.minimizer), " / ", chi2_exp, " (dof: ", length(dat) - param,")") println("Chisq / chiexp: ", min_fun(sol.minimizer), " / ", chi2_exp, " (dof: ", length(dat) - param,")")
chis2_corrected = (length(dat) - param) * min_fun(sol.minimizer) / chi2_exp chis2_corrected = (length(dat) - param) * min_fun(sol.minimizer) / chi2_exp
println("Chisq corrected: ", chis2_corrected) println("Chisq corrected: ", chis2_corrected)
return upar, chis2_corrected @info("fit_routine is returning par, chi2 / chi2_exp\n")
return upar, min_fun(sol.minimizer) / chi2_exp
end end
function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3;
...@@ -1017,7 +1111,10 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal} ...@@ -1017,7 +1111,10 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
end end
function global_fit_routine(models::Vector{Function}, xdata::Vector{Array{Float64, 1}}, ydata::Vector{Array{uwreal, 1}}, n_par::Int64; guess_param::Vector{Float64}=fill(0.5, n_par), log::Bool=true, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) function global_fit_routine(models::Vector{Function}, xdata::Vector{Array{Float64, 1}} , ydata::Vector{Array{uwreal, 1}}, n_par::Int64; guess_param::Vector{Float64}=fill(0.5, n_par), log::Bool=true,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
inv_cov::Union{Matrix{Float64}, Nothing}=nothing)
if !(length(models) == length(ydata)) if !(length(models) == length(ydata))
error("Dimension mismatch") error("Dimension mismatch")
end end
...@@ -1046,28 +1143,52 @@ function global_fit_routine(models::Vector{Function}, xdata::Vector{Array{Float6 ...@@ -1046,28 +1143,52 @@ function global_fit_routine(models::Vector{Function}, xdata::Vector{Array{Float6
return fx return fx
end end
function flat_chisq(param, data) function flat_chisq(param, data)
chi = 0.0 if isnothing(inv_cov)
for ii in 1:N_sets chi = 0.0
mod = models[ii] for ii in 1:N_sets
# uncorrelated chi2 mod = models[ii]
chi += sum((data[chunk_idx[ii]] .- mod(x_flat[chunk_idx[ii]], param)).^2 ./ err.(y_flat[chunk_idx[ii]]).^2 ) chi += sum((data[chunk_idx[ii]] .- mod(x_flat[chunk_idx[ii]], param)).^2 ./ err.(y_flat[chunk_idx[ii]]).^2 )
# correlated chi2 end
return chi
else
# take pieces of covariance and build correlated chi2 piece by piec
# chi =0.0
# for ii in 1:N_sets
# mod = models[ii]
# cov_inv = juovs.inv_covar(ydata[ii])
# res = ydata[ii] .- mod(xdata[ii], param)
# chi += (reshape(res, 1,:) * cov_inv * res)[1]
# end
# return chi
# take whole covariance and multiply for whole vector of dat-mod(x,param)
res = Vector{uwreal}(undef, 0)
for ii in 1:N_sets
mod = models[ii]
res =vcat(res, data[chunk_idx[ii]] .- mod(x_flat[chunk_idx[ii]], param))
end
chi = reshape(res, 1,:) * inv_cov * res
return chi[1]
end end
return chi
end end
fit = curve_fit(flat_model, x_flat, value.(y_flat), 1.0 ./ err.(y_flat).^2, guess_param )
#uncorrelated chi2exp if isnothing(inv_cov)
(upar, chi2_exp) = isnothing(wpm) ? fit_error(flat_chisq, coef(fit) , y_flat) : fit_error(flat_chisq, coef(fit) , y_flat, wpm) @info("Uncorrelated fit")
#correlated chi2exp fit = curve_fit(flat_model, x_flat, value.(y_flat), 1.0 ./ err.(y_flat).^2, guess_param )
#(upar, chi2_exp) = isnothing(wpm) ? fit_error(flat_chisq, coef(fit) , y_flat, W=cov(y_flat)^(-1)) : fit_error(flat_chisq, coef(fit) , y_flat, wpm) (upar, chi2_exp) = isnothing(wpm) ? fit_error(flat_chisq, coef(fit) , y_flat) : fit_error(flat_chisq, coef(fit) , y_flat, wpm)
chisq = value(flat_chisq(coef(fit), y_flat)) else
@info("Correlated fit")
fit = curve_fit(flat_model, x_flat, value.(y_flat), inv_cov, guess_param )
(upar, chi2_exp) = isnothing(wpm) ? fit_error(flat_chisq, coef(fit) , y_flat, W=inv_cov) : fit_error(flat_chisq, coef(fit) , y_flat, wpm, W=inv_cov)
end
println("chi2 from fit residue: ", sum(fit.resid.^2))
chisq = flat_chisq(coef(fit), y_flat)
if log if log
println("Converged param: ", fit.converged) println("Converged param: ", fit.converged)
println("Chisq / chiexp: ", chisq, " / ", chi2_exp, " (dof: ", dof(fit),")") println("Chisq / chiexp: ", chisq, " / ", chi2_exp, " (dof: ", dof(fit),")")
chis2_corrected = (dof(fit)) * chisq / chi2_exp chis2_corrected = (dof(fit)) * chisq / chi2_exp
println("Chisq corrected: ", chis2_corrected) println("Chisq corrected: ", chis2_corrected)
end end
return upar, chisq, chi2_exp, dof(fit) return upar, value(chisq), chi2_exp, dof(fit)
end 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
...@@ -1078,8 +1199,8 @@ function gen_chisq_correlated(f::Function, x::Array{<:Real}, covar_inv::Matrix, ...@@ -1078,8 +1199,8 @@ function gen_chisq_correlated(f::Function, x::Array{<:Real}, covar_inv::Matrix,
chi2 = 0.0 chi2 = 0.0
#delta = dat - f(x, par) #delta = dat - f(x, par)
#chi2 = (reshape(delta, 1,:) * covar_inv * delta )[1] #chi2 = (reshape(delta, 1,:) * covar_inv * delta )[1]
for i in 1:length(dat) for i in eachindex(dat)
for j in 1:length(dat) for j in eachindex(dat)
#println(f(x,par)[i]) #println(f(x,par)[i])
#println(f(x,par)[j]) #println(f(x,par)[j])
#println(f(x[j,:],par)) #println(f(x[j,:],par))
...@@ -1138,4 +1259,41 @@ function get_chi2_cov(f::Function, data, C, par, Nalpha) # full + cov ...@@ -1138,4 +1259,41 @@ function get_chi2_cov(f::Function, data, C, par, Nalpha) # full + cov
chi2 += delta' * Cinv * delta chi2 += delta' * Cinv * delta
end end
return chi2 return chi2
end
function get_covar(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
isnothing(wpm) ? uwerr.(obs) : [uwerr(yaux, wpm) for yaux in obs]
ens = vcat(ensembles.(obs)...)
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)
isnothing(wpm) ? uwerr.(obs) : [uwerr(yaux, wpm) for yaux in obs]
id_set = getindex.(ensembles.(obs), 2) # hardcoded assuming that ens label is always the second entry in each ensembles(obs[i])
cov_tot = zeros(length(id_set), length(id_set))
count = 1
for (k,id) in enumerate(unique(id_set))
idx = findall(x-> x==id, id_set )
fluc_aux = mchist.(obs[idx], fill(id, length(idx)))
fluc = mapreduce(permutedims, vcat, fluc_aux)'
covar_set = fluc' * fluc ./ size(fluc)[1]^2
dim_tmp = size(covar_set,1)
cov_tot[count:dim_tmp+count-1, count:dim_tmp+count-1] = covar_set
count +=dim_tmp
end
return cov_tot
end
function inv_covar(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
cov = get_covar(obs, wpm=wpm)
cov_inv = inv(cov)
return (cov_inv' + cov_inv) / 2.
end
function inv_covar_multi_id(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
cov = get_covar_multi_id(obs, wpm=wpm)
cov_inv = inv(cov)
return (cov_inv' + cov_inv) / 2.
end end
\ No newline at end of file
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