Skip to content

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

update changes

parent ef788cef
Hide whitespace changes
Inline Side-by-side
Showing with 0 additions and 710 deletions
src/juobs_obsmod.jl deleted 100644 → 0
View file @ ef788cef
@doc raw"""
meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false )
meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false)
Computes effective mass for a given correlator corr at a given plateau plat.
Correlator can be passed as an Corr struct or Vector{uwreal}.
The flags pl and data allow to show the plots and return data as an extra result.
```@example
data = read_mesons(path, "G5", "G5")
corr_pp = corr_obs.(data)
m = meff(corr_pp[1], [50, 60], pl=false)
```
"""
function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, mu::Union{Vector{Float64}, Nothing}=nothing)
dim = length(corr)
aux = 0.5 .* log.((corr[2:dim-2] ./ corr[3:dim-1]).^2)
mass = plat_av(aux, plat)
uwerr(mass)
if pl == true
x = 1:length(aux)
y = value.(aux)
dy = err.(aux)
v = value(mass)
e = err(mass)
figure()
fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)
errorbar(x, y, dy, fmt="x", color="black")
ylabel(L"$m_\mathrm{eff}$")
xlabel(L"$x_0$")
if !isnothing(mu)
title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
end
display(gcf())
end
if data == false
return mass
else
return (mass, aux)
end
end
meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false) =
meff(corr.obs, plat, pl=pl, data=data, mu=corr.mu)
## Decay constants
@doc raw"""
dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false)meff(corr::Corr, plat::Vector{Int64}; pl::Bool=true, data::Bool=false)
dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false)
Computes decay constant using the PCVC relation for twisted mass fermions. The decay constant is computed in the plateau plat.
Correlator can be passed as an Corr struct or Vector{uwreal}. If it is passed as a uwreal vector, vector of twisted masses mu and source position y0
must be specified.
The flags pl and data allow to show the plots and return data as an extra result.
```@example
data = read_mesons(path, "G5", "G5")
corr_pp = corr_obs.(data)
m = meff(corr_pp[1], [50, 60], pl=false)
f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)
```
"""
function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false)
"""
compute the decay constant when the source is far from the boundaries
"""
corr_pp = corr[2:end-1]
dim = length(corr_pp)
aux = exp.((collect(1:dim) .- y0 ) .* [m for k in 1:dim])
R = (aux .* corr_pp).^0.5
R_av = plat_av(R, plat)
f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5
uwerr(f)
if pl == true
uwerr(R_av)
v = value(R_av)
e = err(R_av)
uwerr.(R)
figure()
fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)
errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black")
ylabel(L"$R$")
xlabel(L"$x_0$")
title(string(L"$\mu_1 = $", mu[1], L" $\mu_2 = $", mu[2]))
display(gcf())
end
if data == false
return f
else
return (f, R)
end
end
dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::Bool=false) =
dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data)
#t0
@doc raw"""
comp_t0(y::YData, plat::Vector{Int64}; pl::Bool=false, L::Int64=1)
comp_t0(Y::Vector{YData}, plat::Vector{Int64}; pl::Bool=false, L::Int64=1)
Computes t0 using the energy density of the action Ysl(Yang-Mills action).
t0 is computed in the plateau plat.
A (linear) interpolation in t is performed to find t0.
The flag pl allows to show the plot.
```@example
#Single replica
Y = read_ms(path)
rw = read_ms(path_rw)
t0 = comp_t0(Y, [38, 58], L=32)
t0_r = comp_t0(Y, [38, 58], L=32, rw=rw)
#Two replicas
Y1 = read_ms(path1)
Y2 = read_ms(path2)
rw1 = read_ms(path_rw1)
rw2 = read_ms(path_rw2)
t0 = comp_t0([Y1, Y2], [38, 58], L=32, pl=true)
t0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true)
```
"""
function comp_t0(Y::YData, plat::Vector{Int64}; pl::Bool=false, L::Int64=1, rw::Union{Matrix{Float64}, Nothing}=nothing)
Ysl = Y.Ysl
t = Y.t
id = Y.id
Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw)
xmax = size(mean(Ysl, dims=3))[1]
nt0 = t0_guess(t, Ysl, plat, L)
Y = Matrix{uwreal}(undef, xmax, 3)
for i = 1:xmax
k = 1
for j = nt0-1:nt0+1
Y[i, k] = uwreal(Ysl[i, j, :], id)
k = k + 1
end
end
x = t[nt0-1:nt0+1]
t2E = [plat_av(Y[:, j], plat) for j=1:3] .* x.^2 / L^3
par, chi2exp = lin_fit(x, t2E)
t0 = x_lin_fit(par, 0.3)
if pl
uwerr(t0)
uwerr.(t2E)
v = value.(t2E)
e = err.(t2E)
figure()
errorbar(x, v, e, fmt="x")
errorbar(value(t0), 0.3, xerr=err(t0), fmt="x")
ylabel(L"$t^2E$")
xlabel(L"$t/a^2$")
display(gcf())
end
return t0
end
function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; pl::Bool=false, L::Int64=1, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing)
nr = length(Y)
Ysl = getfield.(Y, :Ysl)
t = getfield.(Y, :t)
id = getfield.(Y, :id)
vtr = getfield.(Y, :vtr)
replica = length.(vtr)
if !all(id .== id[1])
println("IDs are not equal")
return nothing
end
Ysl = isnothing(rw) ? Ysl : apply_rw.(Ysl, rw)
xmax = size(mean(Ysl[1], dims=3))[1]
tmp = Ysl[1]
[tmp = cat(tmp, Ysl[k], dims=3) for k = 2:nr]
nt0 = t0_guess(t[1], tmp, plat, L)
Y = Matrix{uwreal}(undef, xmax, 3)
for i = 1:xmax
k = 1
for j = nt0-1:nt0+1
Y[i, k] = uwreal(tmp[i, j, :], id[1], replica)
k = k + 1
end
end
x = t[1][nt0-1:nt0+1]
t2E = [plat_av(Y[:, j], plat) for j=1:3] .* x.^2 / L^3
println(t2E)
par, chi2exp = lin_fit(x, t2E)
t0 = x_lin_fit(par, 0.3)
if pl
uwerr(t0)
uwerr.(t2E)
v = value.(t2E)
e = err.(t2E)
figure()
errorbar(x, v, e, fmt="x")
errorbar(value(t0), 0.3, xerr=err(t0), fmt="x")
ylabel(L"$t^2E$")
xlabel(L"$t/a^2$")
display(gcf())
end
return t0
end
function t0_guess(t::Vector{Float64}, Ysl::Array{Float64, 3}, plat::Vector{Int64}, L::Int64)
t2E_ax = t.^2 .* mean(mean(Ysl[plat[1]:plat[2], :, :], dims=1), dims=3)[1, :, 1] / L^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
return nt0
end
\ No newline at end of file
src/juobs_readermod.jl deleted 100644 → 0
View file @ ef788cef
function read_global_header(path::String)
data = open(path, "r")
g_header = zeros(Int32, 4)
read!(data, g_header)
close(data)
a = GHeader(g_header)
return a
end
function read_CHeader(path::String)
gh = read_global_header(path)
data = open(path, "r")
seek(data, gh.hsize)
aux_f = zeros(Float64, 6)
aux_i = zeros(Int32, 4)
theta = zeros(Float64, 6)
a = Vector{CHeader}(undef, gh.ncorr)
for k = 1:gh.ncorr
read!(data, aux_f)
read!(data, theta)
qs1 = read(data, Int32)
if qs1 != 0
qn1 = read(data, Int32)
qeps1 = read(data, Float64)
q1 = Sm(qs1, qn1, qeps1, 1)
else
q1 = Sm(qs1, 1)
end
qs2 = read(data, Int32)
if qs2 != 0
qn2 = read(data, Int32)
qeps2 = read(data, Float64)
q2 = Sm(qs2, qn2, qeps2, 1)
else
q2 = Sm(qs2, 1)
end
gs1 = read(data, Int32)
if gs1 != 0 && gs1 != 3 && gs1 != 4
gn1 = read(data, Int32)
geps1 = read(data, Float64)
g1 = Sm(gs1, gn1, geps1, 2)
elseif gs1 == 3 || gs1 == 4
g1 = Sm(gs1, q1.niter, q1.eps, 2)
else
g1 = Sm(gs1, 2)
end
gs2 = read(data, Int32)
if gs2 != 0 && gs2 != 3 && gs2 != 4
gn2 = read(data, Int32)
geps2 = read(data, Float64)
g2 = Sm(gs2, gn2, geps2, 2)
elseif gs1 == 3 || gs1 == 4
g2 = Sm(gs2, q2.niter, q2.eps, 2)
else
g2 = Sm(gs2, 2)
end
read!(data, aux_i)
a[k] = CHeader(aux_f, aux_i, theta, [q1, q2, g1, g2])
end
close(data)
return a
end
@doc raw"""
read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
This faction read a mesons dat file at a given path and returns a vector of CData structures for different masses and Dirac structures.
Dirac structures g1 and/or g2 can be passed as string arguments in order to filter correaltors.
ADerrors id can be specified as argument. If is not specified, the id is fixed according to the ensemble name (example: "H400"-> id = 400)
Examples:
```@example
read_mesons(path)
read_mesons(path, "G5")
read_mesons(path, nothing, "G5")
read_mesons(path, "G5", "G5")
read_mesons(path, "G5", "G5", id=1)
```
"""
function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
isnothing(g1) ? t1=nothing : t1 = findfirst(x-> x==g1, gamma_name) - 1
isnothing(g2) ? t2=nothing : t2 = findfirst(x-> x==g2, gamma_name) - 1
if isnothing(id)
bname = basename(path)
m = findfirst(r"[A-Z][0-9]{3}r[0-9]{3}", bname)
id = parse(Int64, bname[m[2:4]])
end
data = open(path, "r")
g_header = read_global_header(path)
c_header = read_CHeader(path)
ncorr = g_header.ncorr
tvals = g_header.tvals
nnoise = g_header.nnoise
fsize = filesize(path)
datsize = 4 + sum(c.dsize for c in c_header)*tvals*nnoise #datasize / ncnfg
ncfg = Int32((fsize-g_header.hsize-sum(c.hsize for c in c_header)) / datsize)
corr_match = findall(x-> (x.type1==t1 || isnothing(t1)) && (x.type2==t2 || isnothing(t2)), c_header)
seek(data, g_header.hsize + sum(c.hsize for c in c_header))
res = Array{CData}(undef, length(corr_match))
data_re = Array{Float64}(undef, length(corr_match), ncfg, tvals)
data_im = zeros(length(corr_match), ncfg, tvals)
vcfg = Array{Int32}(undef, ncfg)
for icfg = 1:ncfg
vcfg[icfg] = read(data, Int32)
c=1
for k = 1:ncorr
if k in corr_match
if c_header[k].is_real == 1
tmp = Array{Float64}(undef, tvals*nnoise)
read!(data, tmp)
tmp2 = reshape(tmp, (nnoise, tvals))
tmp2 = mean(tmp2, dims=1)
data_re[c, icfg, :] = tmp2[1, :]
elseif c_header[k].is_real == 0
tmp = Array{Float64}(undef, 2*tvals*nnoise)
read!(data, tmp)
tmp2 = reshape(tmp, (2, nnoise, tvals))
tmp2 = mean(tmp2, dims=2)
data_re[c, icfg, :] = tmp2[1, 1, :]
data_im[c, icfg, :] = tmp2[2, 1, :]
end
c += 1
else
seek(data, position(data) + c_header[k].dsize*tvals*nnoise)
end
end
end
for c in 1:length(corr_match)
res[c] = CData(c_header[corr_match[c]], vcfg, data_re[c, :, :], data_im[c, :, :], id)
end
close(data)
return res
end
@doc raw"""
read_ms1(path::String; v::String="1.2")
Reads openQCD ms1 dat files at a given path. This method returns a matrix W[irw, icfg] that contains the reweighting factors, where
irw is the rwf index and icfg the configuration number.
The function is compatible with the output files of openQCD v=1.2, 1.4 and 1.6. Version can be specified as argument.
Examples:
```@example
read_ms1(path)
read_ms1(path, v="1.4")
read_ms1(path, v="1.6")
```
"""
function read_ms1(path::String; v::String="1.2")
data = open(path, "r")
nrw = read(data, Int32)
if v == "1.4" || v =="1.6"
nfct = Array{Int32}(undef, nrw)
read!(data, nfct)
nfct_inheader = 1
elseif v=="1.2"
nfct = ones(Int32, nrw)
nfct_inheader = 0
else
println("Error: Version not supported")
return nothing
end
nsrc = Array{Int32}(undef, nrw)
read!(data, nsrc)
glob_head_size = 4 + 4*nrw*(nfct_inheader + 1)
datsize = 4 + 2*8*sum(nsrc .* nfct)
fsize = filesize(path)
ncnfg = Int32((fsize - glob_head_size)/datsize)
r_data = Array{Array{Float64}}(undef, nrw)
[r_data[k] = zeros(Float64, nfct[k], nsrc[k], ncnfg) for k = 1:nrw]
vcfg = Array{Int32}(undef, ncnfg)
for icfg = 1:ncnfg
vcfg[icfg] = read(data, Int32)
for irw = 1:nrw
for ifct = 1:nfct[irw]
tmp = zeros(Float64, nsrc[irw])
seek(data, position(data) + 8 * nsrc[irw])
read!(data, tmp)
r_data[irw][ifct, :, icfg] = tmp[:]
end
end
end
W = zeros(Float64, nrw, ncnfg)
[W[k, :] = prod(mean(exp.(.-r_data[k]), dims=2), dims=1) for k = 1:nrw]
close(data)
return W
end
@doc raw"""
read_ms(path::String)
Reads openQCD ms dat files at a given path. This method return:
t(t): flow time values
Wsl(x0, t, icfg): the time-slice sums of the densities of the Wilson plaquette action
Ysl(x0, t, icfg): the time-slice sums of the densities of the Yang-Mills action
Qsl(x0, t, icfg): the time-slice sums of the densities of the topological charge
Examples:
```@example
t, W, Y, Q = read_ms(path)
```
"""
function read_ms(path::String; id::Union{Int64, Nothing}=nothing)
if isnothing(id)
bname = basename(path)
m = findfirst(r"[A-Z][0-9]{3}r[0-9]{3}", bname)
id = parse(Int64, bname[m[2:4]])
end
data = open(path, "r")
dn = read(data, Int32)
nn = read(data, Int32)
tvals = read(data, Int32)
eps = read(data, Float64)
fsize = filesize(path)
datsize=4 + 3*8*(nn + 1) * tvals # measurement size of each cnfg
ntr = Int32((fsize - 3*4 - 8) / datsize)
vntr = Vector{Int32}(undef, ntr)
# x0, t, cfg
Wsl = Array{Float64}(undef, tvals, nn + 1, ntr)
Ysl = Array{Float64}(undef, tvals, nn + 1, ntr)
Qsl = Array{Float64}(undef, tvals, nn + 1, ntr)
for itr = 1:ntr
vntr[itr] = read(data, Int32)
for iobs = 1:3
for inn = 0:nn
tmp = Vector{Float64}(undef, tvals)
read!(data, tmp)
if iobs == 1
Wsl[:, inn + 1, itr] = tmp
elseif iobs == 2
Ysl[:, inn + 1, itr] = tmp
elseif iobs == 3
Qsl[:, inn + 1, itr] = tmp
end
end
end
end
close(data)
t = Float64.(0:nn) .* dn .* eps
return YData(vntr, t, Ysl, id)
end
src/juobs_typesmod.jl deleted 100644 → 0
View file @ ef788cef
##############UTILITY NAMES##############
#= STARTING AT 0
noise_name=['Z2', 'GAUSS', 'U1', 'Z4']
gamma_name = ['G0', 'G1', 'G2', 'G3',
'invalid', 'G5', '1', 'G0G1', 'G0G2', 'G0G3',
'G0G5', 'G1G2', 'G1G3', 'G1G5', 'G2G3', 'G2G5', 'G3G5']
# quark and gauge Smearing types
qs_name=['Local', 'Wuppertal', '3D Gradient Flow', 'Gradient Flow']
gs_name=['Local', 'APE', '3D Wilson Flow', 'Quark 3D Gradient Flow', 'Quark Gradient Flow']
=#
const noise_name=["Z2", "GAUSS", "U1", "Z4"]
const gamma_name = ["G0", "G1", "G2", "G3",
"invalid", "G5", "1", "G0G1", "G0G2", "G0G3",
"G0G5", "G1G2", "G1G3", "G1G5", "G2G3", "G2G5", "G3G5"]
const qs_name=["Local", "Wuppertal", "3D Gradient Flow", "Gradient Flow"]
const gs_name=["Local", "APE", "3D Wilson Flow", "Quark 3D Gradient Flow", "Quark Gradient Flow"]
mutable struct GHeader
ncorr::Int32
nnoise::Int32
tvals::Int32
noisetype::Int32
hsize::Int32 #headersize
GHeader(a, b, c, d) = new(a, b, c, d, 4*4)
function GHeader(x::Vector{Int32})
a = new(x[1], x[2], x[3], x[4], 4*4)
return a
end
end
mutable struct Sm
type::Int32
niter::Int32
eps::Float64
qg::Int32 #1->fermionic 2->gluonic
Sm(t,q) = new(t, 0, 0.0, q)
Sm(t, n, e, q) = new(t, n, e, q)
end
mutable struct CHeader
k1::Float64
k2::Float64
mu1::Float64
mu2::Float64
dp1::Float64
dp2::Float64
theta1::Vector{Float64}
theta2::Vector{Float64}
q1::Sm
q2::Sm
g1::Sm
g2::Sm
type1::Int32
type2::Int32
x0::Int32
is_real::Int32
hsize::Int32 #headersize
dsize::Int32 #datasize / (nnoise * T * ncfg)
function CHeader(aux_f::Vector{Float64}, aux_i::Vector{Int32}, theta::Vector{Float64}, sm_par::Vector{Sm})
a = new()
a.k1 = aux_f[1]
a.k2 = aux_f[2]
a.mu1 = aux_f[3]
a.mu2 = aux_f[4]
a.dp1 = aux_f[5]
a.dp2 = aux_f[6]
a.type1 = aux_i[1]
a.type2 = aux_i[2]
a.x0 = aux_i[3]
a.is_real = aux_i[4]
a.theta1 = theta[1:3]
a.theta2 = theta[4:6]
a.q1 = sm_par[1]
a.q2 = sm_par[2]
a.g1 = sm_par[3]
a.g2 = sm_par[4]
a.hsize = 8*12 + 4*8 #without smearing parameters niter, neps
if a.q1.type != 0
a.hsize += 8+4
end
if a.q2.type != 0
a.hsize += 8+4
end
if a.g1.type != 0 && a.g1.type != 3 && a.g1.type != 4
a.hsize += 8+4
end
if a.g2.type != 0 && a.g2.type != 3 && a.g2.type != 4
a.hsize += 8+4
end
a.dsize = 16 - 8* a.is_real
return a
end
end
mutable struct CData
header::CHeader
vcfg::Array{Int32}
re_data::Array{Float64}
im_data::Array{Float64}
id::Int64
CData(a, b, c, d, e) = new(a, b, c, d, e)
end
Base.copy(a::CData) = CData(a.header, a.vcfg, a.re_data, a.im_data, a.id)
mutable struct Corr
obs::Vector{uwreal}
mu::Vector{Float64}
gamma::Vector{String}
y0::Int64
function Corr(a::Vector{uwreal}, b::CData)
h = getfield(b, :header)
mu = [h.mu1, h.mu2]
gamma = [gamma_name[h.type1+1], gamma_name[h.type2+1]]
y0 = Int64(h.x0)
return new(a, mu, gamma, y0)
end
function Corr(a::Vector{uwreal}, b::Vector{CData})
sym = [:mu1, :mu2, :type1, :type2, :x0]
h = getfield.(b, :header)
for s in sym
if !all(getfield.(h, s) .== getfield(h[1], s))
println("Corr: Parameter mismatch")
return nothing
end
end
mu = [h[1].mu1, h[1].mu2]
gamma = [gamma_name[h[1].type1+1], gamma_name[h[1].type2+1]]
y0 = Int64(h[1].x0)
return new(a, mu, gamma, y0)
end
end
mutable struct YData
vtr::Vector{Int32}
t::Vector{Float64}
Ysl::Array{Float64, 3}
id::Int64
YData(a, b, c, d) = new(a, b, c, d)
end
function Base.show(io::IO, a::GHeader)
print(io, "ncorr = ", a.ncorr, "\t")
print(io, "nnoise = ", a.nnoise, "\t")
print(io, "tvals = ", a.tvals, "\t")
print(io, "noisetype = ", noise_name[a.noisetype + 1], "\t")
print(io, "\n")
end
function Base.show(io::IO, a::CHeader)
fnames = fieldnames(CHeader)
for k in fnames
f = getfield(a, k)
if k != :type1 && k!= :type2
print(io, "$k = $f\t")
else
print(io, "$k = ", gamma_name[f + 1], "\t")
end
end
print(io, "\n")
end
function Base.show(io::IO, a::Sm)
print(io, "\n")
if a.qg == 1
print(io, "type = ", qs_name[a.type+1], "\t")
elseif a.qg == 2
print(io, "type = ", gs_name[a.type+1], "\t")
end
print(io, "niter = ", a.niter, "\t")
print(io, "eps = ", a.eps, "\t")
print(io, "\n")
end
function Base.show(io::IO, a::CData)
print(io, a.header)
end
\ No newline at end of file
src/juobsmod.jl deleted 100644 → 0
View file @ ef788cef
module juobs
using ADerrors, PyPlot, Statistics, LaTeXStrings, LinearAlgebra, LsqFit
include("juobs_types.jl")
include("juobs_linalg.jl")
include("juobs_reader.jl")
include("juobs_tools.jl")
include("juobs_obs.jl")
export read_mesons, read_ms1, read_ms
export get_matrix, uwgevp_tot, energies, uwdot
export corr_obs, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
export meff, dec_const_pcvc, comp_t0
end # module
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