Skip to content

J
juobs
Commit 9a0d2137 authored by Javier's avatar Javier
Browse files
Options

read_ms + comp_t0 modified 

indexes permuted in output read_ms
comp_t0 fixes rwf, plot plateau added, L mandatory
parent 55c47ca9
Hide whitespace changes
Inline Side-by-side
Showing with 45 additions and 30 deletions
src/juobs_obs.jl
View file @ 9a0d2137
......@@ -102,9 +102,9 @@ dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::
#t0
@doc raw"""
comp_t0(y::YData, plat::Vector{Int64}; pl::Bool=false, L::Int64=1)
comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing)
comp_t0(Y::Vector{YData}, plat::Vector{Int64}; pl::Bool=false, L::Int64=1)
comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing)
Computes t0 using the energy density of the action Ysl(Yang-Mills action).
t0 is computed in the plateau plat.
......@@ -131,26 +131,25 @@ 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)
function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, 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]
xmax = size(Ysl, 2)
nt0 = t0_guess(t, Ysl, plat, L)
Y = Matrix{uwreal}(undef, xmax, 3)
Y_aux = 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)
Y_aux[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
t2E = [plat_av(Y_aux[:, j], plat) for j=1:3] .* x.^2 / L^3
par, chi2exp = lin_fit(x, t2E)
......@@ -167,15 +166,25 @@ function comp_t0(Y::YData, plat::Vector{Int64}; pl::Bool=false, L::Int64=1, rw::
errorbar(value(t0), 0.3, xerr=err(t0), fmt="x")
ylabel(L"$t^2E$")
xlabel(L"$t/a^2$")
display(gcf())
display(gcf())
figure()
plot(collect(1:xmax), value.(Y_aux[:, 2]) * t[nt0]^2 / L^3, "x")
fill_between(collect(plat[1]:plat[2]), v[2]+e[2], v[2]-e[2], alpha=0.75, color="green")
ylabel(L"$t^2E(x0, t)$")
xlabel(L"$x_0/a$")
title(string(L"$t/a^2 = $", t[nt0]))
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)
function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing)
nr = length(Y)
Ysl = getfield.(Y, :Ysl)
t = getfield.(Y, :t)
t = t[1]
id = getfield.(Y, :id)
vtr = getfield.(Y, :vtr)
replica = length.(vtr)
......@@ -188,23 +197,21 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; pl::Bool=false, L::Int64
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)
[tmp = cat(tmp, Ysl[k], dims=1) for k = 2:nr]
nt0 = t0_guess(t, tmp, plat, L)
xmax = size(tmp, 2)
Y = Matrix{uwreal}(undef, xmax, 3)
Y_aux = 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)
Y_aux[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)
x = t[nt0-1:nt0+1]
t2E = [plat_av(Y_aux[:, j], plat) for j=1:3] .* x.^2 / L^3
par, chi2exp = lin_fit(x, t2E)
t0 = x_lin_fit(par, 0.3)
......@@ -220,12 +227,20 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; pl::Bool=false, L::Int64
errorbar(value(t0), 0.3, xerr=err(t0), fmt="x")
ylabel(L"$t^2E$")
xlabel(L"$t/a^2$")
display(gcf())
display(gcf())
figure()
plot(collect(1:xmax), value.(Y_aux[:, 2]) * t[nt0]^2 / L^3, "x")
fill_between(collect(plat[1]:plat[2]), v[2]+e[2], v[2]-e[2], alpha=0.75, color="green")
ylabel(L"$t^2E(x0, t)$")
xlabel(L"$x_0/a$")
title(string(L"$t/a^2 = $", t[nt0]))
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
t2E_ax = t.^2 .* mean(mean(Ysl[:, plat[1]:plat[2], :], dims=2), dims=1)[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
......
src/juobs_reader.jl
View file @ 9a0d2137
......@@ -245,9 +245,9 @@ 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
Wsl(icfg, x0, t): the time-slice sums of the densities of the Wilson plaquette action
Ysl(icfg, x0, t): the time-slice sums of the densities of the Yang-Mills action
Qsl(icfg, x0, t): the time-slice sums of the densities of the topological charge
Examples:
```@example
......@@ -275,9 +275,9 @@ function read_ms(path::String; id::Union{Int64, Nothing}=nothing)
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)
Wsl = Array{Float64}(undef, ntr, tvals, nn + 1)
Ysl = Array{Float64}(undef, ntr, tvals, nn + 1)
Qsl = Array{Float64}(undef, ntr, tvals, nn + 1)
for itr = 1:ntr
vntr[itr] = read(data, Int32)
......@@ -286,11 +286,11 @@ function read_ms(path::String; id::Union{Int64, Nothing}=nothing)
tmp = Vector{Float64}(undef, tvals)
read!(data, tmp)
if iobs == 1
Wsl[:, inn + 1, itr] = tmp
Wsl[itr, :, inn + 1] = tmp
elseif iobs == 2
Ysl[:, inn + 1, itr] = tmp
Ysl[itr, :, inn + 1] = tmp
elseif iobs == 3
Qsl[:, inn + 1, itr] = tmp
Qsl[itr, :, inn + 1] = tmp
end
end
......
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