Skip to content

J
juobs
Commit 54e7b8a7 authored by Javier's avatar Javier
Browse files
Options

Merge branch 'javier'

parents 55c47ca9 15e017ec
Hide whitespace changes
Inline Side-by-side
Showing with 200 additions and 72 deletions
src/juobs.jl
View file @ 54e7b8a7
module juobs
using ADerrors, PyPlot, Statistics, LaTeXStrings, LinearAlgebra, LsqFit
using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit
import Statistics: mean
include("juobs_types.jl")
include("juobs_linalg.jl")
......@@ -9,7 +10,7 @@ include("juobs_obs.jl")
export read_mesons, read_ms1, read_ms, read_md
export get_matrix, uwgevp_tot, energies, uwdot
export corr_obs, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
export corr_obs, md_sea, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
export meff, dec_const_pcvc, comp_t0
end # module
src/juobs_linalg.jl
View file @ 54e7b8a7
......@@ -56,8 +56,7 @@ function get_matrix(corr_diag::Vector{Vector{uwreal}}, corr_upper::Vector{Vector
n = length(corr_diag) # matrix dimension
d = length(corr_upper) # upper elements
if d != (n*(n-1) / 2)
println("Error get_matrix")
return nothing
error("Error get_matrix")
end
res = Vector{Matrix}(undef, time) # array with all the matrices at each time step.
aux = Matrix{uwreal}(undef, n, n)
......@@ -388,8 +387,7 @@ function uwdot(a::Matrix{uwreal}, b::Matrix{uwreal})
k,p = size(b)
c = Matrix{uwreal}(undef, n, p)
if m != k
println("Error uwdot")
return nothing
error("Error uwdot")
end
for i = 1:n
for j = 1:p
......
src/juobs_obs.jl
View file @ 54e7b8a7
......@@ -100,15 +100,22 @@ dec_const_pcvc(corr::Corr, plat::Vector{Int64}, m::uwreal; pl::Bool=true, data::
dec_const_pcvc(corr.obs, plat, m, corr.mu, corr.y0, pl=pl, data=data)
#t0
function get_model(x, p, n)
s = 0.0
for k = 1:n
s = s .+ p[k] .* x.^(k-1)
end
return s
end
@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, npol::Int64=2)
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, npol::Int64=2)
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.
A polynomial interpolation in t is performed to find t0, where npol is the degree of the polynomial (linear fit by default)
The flag pl allows to show the plot.
......@@ -131,30 +138,32 @@ 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, npol::Int64=2)
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)
dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1)/ 2)
Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
for i = 1:xmax
k = 1
for j = nt0-1:nt0+1
Y[i, k] = uwreal(Ysl[i, j, :], id)
for j = nt0-dt0:nt0+dt0
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
x = t[nt0-dt0:nt0+dt0]
t2E = [plat_av(Y_aux[:, j], plat) for j=1:2*dt0+1] .* x.^2 / L^3
par, chi2exp = lin_fit(x, t2E)
model(x, p) = get_model(x, p, npol)
t0 = x_lin_fit(par, 0.3)
par = fit_routine(model, x, t2E, npol)
fmin(x, p) = model(x, p) .- 0.3
t0 = root_error(fmin, t[nt0], par)
if pl
uwerr(t0)
uwerr.(t2E)
......@@ -167,47 +176,56 @@ 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()
t_pl = dt0 + 1
plot(collect(1:xmax), value.(Y_aux[:, t_pl]) * t[nt0]^2 / L^3, "x")
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)$")
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, npol::Int64=2)
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)
if !all(id .== id[1])
println("IDs are not equal")
return nothing
error("IDs are not equal")
end
Ysl = isnothing(rw) ? Ysl : apply_rw.(Ysl, rw)
xmax = size(mean(Ysl[1], dims=3))[1]
Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw)
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)
dt0 = iseven(npol) ? Int64(npol / 2) : Int64((npol+1) / 2)
Y_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
for i = 1:xmax
k = 1
for j = nt0-1:nt0+1
Y[i, k] = uwreal(tmp[i, j, :], id[1], replica)
for j = nt0-dt0:nt0+dt0
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)
par, chi2exp = lin_fit(x, t2E)
x = t[nt0-dt0:nt0+dt0]
t2E = [plat_av(Y_aux[:, j], plat) for j=1:2*dt0+1] .* x.^2 / L^3
model(x, p) = get_model(x, p, npol)
t0 = x_lin_fit(par, 0.3)
par = fit_routine(model, x, t2E, npol)
fmin(x, p) = model(x, p) .- 0.3
t0 = root_error(fmin, t[nt0], par)
if pl
uwerr(t0)
uwerr.(t2E)
......@@ -220,12 +238,21 @@ 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()
t_pl = dt0 + 1
plot(collect(1:xmax), value.(Y_aux[:, t_pl]) * t[nt0]^2 / L^3, "x")
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)$")
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 @ 54e7b8a7
......@@ -168,8 +168,7 @@ function read_rw(path::String; v::String="1.2")
nfct = ones(Int32, nrw)
nfct_inheader = 0
else
println("Error: Version not supported")
return nothing
error("Version not supported")
end
nsrc = Array{Int32}(undef, nrw)
read!(data, nsrc)
......@@ -225,6 +224,8 @@ read_md(path::String)
Reads openQCD pbp.dat files at a given path. This method returns a matrix md[irw, icfg] that contains the derivatives dS/dm, where
md[irw=1] = dS/dm_l and md[irw=2] = dS/dm_s
Seff = -tr(log(D+m))
dSeff/ dm = -tr((D+m)^-1)
Examples:
```@example
......@@ -236,7 +237,7 @@ function read_md(path::String)
nrw = length(r_data)
ncnfg = size(r_data[1])[3]
md = zeros(Float64, nrw, ncnfg)
[md[k, :] = prod(mean(r_data[k], dims=2), dims=1) for k = 1:nrw]
[md[k, :] = prod(.-mean(r_data[k], dims=2), dims=1) for k = 1:nrw]
return md
end
......@@ -245,9 +246,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 +276,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 +287,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
......
src/juobs_tools.jl
View file @ 54e7b8a7
#TODO: apply_rw with gaps
function apply_rw(cdata::CData, W::Array{Float64, 2})
nc = maximum(cdata.vcfg)
function apply_rw(data::Array{Float64}, W::Matrix{Float64})
nc = size(data, 1)
W1 = W[1, 1:nc]
W2 = W[2, 1:nc]
cdata_r = copy(cdata)
cdata_r.re_data = cdata.re_data .* W1 .* W2 / mean(W1 .* W2)
cdata_r.im_data = cdata.im_data .* W1 .* W2 / mean(W1 .* W2)
return cdata_r
data_r = data .* W1 .* W2 / mean(W1 .* W2)
return data_r
end
function apply_rw(data::Array{Float64}, W::Array{Float64, 2})
nc = size(data)[1]
W1 = W[1, 1:nc]
W2 = W[2, 1:nc]
function apply_rw(data::Vector{<:Array{Float64}}, W::Vector{Matrix{Float64}})
if length(W) != length(data)
error("Lenghts must match")
end
nc = size.(data, 1)
data_r = data .* W1 .* W2 / mean(W1 .* W2)
rw1 = [W[k][1, 1:nc[k]] for k=1:length(W)]
rw2 = [W[k][2, 1:nc[k]] for k=1:length(W)]
rw1_cat = rw1[1]
rw2_cat = rw2[1]
for k = 2:length(W)
rw1_cat = cat(rw1_cat, rw1[k], dims=1)
rw2_cat = cat(rw2_cat, rw2[k], dims=1)
end
rw_mean = mean(rw1_cat .* rw2_cat)
data_r = [data[k] .* rw1[k].* rw2[k] / rw_mean for k=1:length(data)]
return data_r
end
@doc raw"""
......@@ -52,7 +62,7 @@ corr_pp_r = corr_obs.(cdata, rw=[rw1, rw2])
function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1)
real ? data = cdata.re_data ./ L^3 : data = cdata.im_data ./ L^3
isnothing(rw) ? data_r = data : data_r = apply_rw(data, rw)
data_r = isnothing(rw) ? data : apply_rw(data, rw)
nt = size(data)[2]
obs = Vector{uwreal}(undef, nt)
......@@ -69,13 +79,11 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
replica = Int64.(maximum.(vcfg))
if !all(id .== id[1])
println("IDs are not equal")
return nothing
error("IDs are not equal")
end
real ? data = getfield.(cdata, :re_data) ./ L^3 : data = getfield.(cdata, :im_data) ./ L^3
isnothing(rw) ? data_r = data : data_r = apply_rw.(data, rw)
data_r = isnothing(rw) ? data : apply_rw(data, rw)
tmp = data_r[1]
[tmp = cat(tmp, data_r[k], dims=1) for k = 2:nr]
......@@ -88,6 +96,100 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
return Corr(obs, cdata)
end
@doc raw"""
md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADerrors.wsg)
Computes the derivative of an observable A with respect to the sea quark masses.
d <A> / dm(sea) = sum_i (d<A> / d<Oi>) * (d<Oi> / dm(sea))
d <O> / dm(sea) = <O> <dS/dm> - <O dS/dm> = - <(O - <O>) (dS/dm - <dS/dm>)>
where <Oi> are primary observables
md is a vector that contains the derivative of the action S with respect
to the sea quark masses for each replica. md[irep][irw, icfg]
md_sea returns a tuple of uwreal observables (dA/dml, dA/dms)|sea,
where ml and ms are the light and strange quark masses.
@example```
#Single replica
data = read_mesons(path, "G5", "G5")
md = read_md(path_md)
corr_pp = corr_obs.(data)
m = meff(corr_pp[1], plat)
m_mdl, m_mds = md_sea(m, [md], ADerrors.wsg)
m_shifted = m + 2 * dml * m_mdl + dms * m_mds
#Two replicas
data_r1 = read_mesons(path_r1, "G5", "G5")
data_r2 = read_mesons(path_r2, "G5", "G5")
md1 = read_md(path_md1)
md2 = read_md(path_md2)
cdata = [[data_r1[k], data_r2[k]] for k=1:length(data_r1)]
corr_pp = corr_obs.(cdata)
m = meff(corr_pp[1], plat)
m_mdl, m_mds = md_sea(m, [md1, md2], ADerrors.wsg)
m_shifted = m + 2 * dml * m_mdl + dms * m_mds
```
"""
#TODO: VECTORIZE
function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADerrors.wsg)
nid = neid(a)
p = findall(t-> t==1, a.prop)
if nid != 1
error("neid > 1")
end
id = ws.map_nob[p]
if !all(id .== id[1])
error("ids do not match")
end
id = id[1]
ivrep = getfield.(ws.fluc[p], :ivrep)
ivrep1 = fill(ivrep[1], length(ivrep))
if !all(ivrep .== ivrep1)
error("ivreps do not match")
end
ivrep = ivrep[1]
if length(md) != length(ivrep)
error("Nr obs != Nr md")
end
md_aux = md[1][:, 1:ivrep[1]]
for k = 2:length(md)
md_aux = cat(md_aux, md[k][:, 1:ivrep[k]], dims=2)
end
fluc_obs = getfield.(ws.fluc[p], :delta)
fluc_md = md_aux .- mean(md_aux, dims=2)
nrw = size(fluc_md, 1)
d = a.der[p]
nobs = sum(a.prop)
if nrw == 1
der1 = uwreal(0)
for k = 1:nobs
der1 = der1 - d[k] * uwreal(fluc_md[1, :] .* fluc_obs[k], id, ivrep)
end
return (der1, der1)
elseif nrw == 2
der1 = uwreal(0)
der2 = uwreal(0)
for k = 1:nobs
der1 = der1 - d[k] * uwreal(fluc_md[1, :] .* fluc_obs[k], id, ivrep)
der2 = der2 - d[k] * uwreal(fluc_md[2, :] .* fluc_obs[k], id, ivrep)
end
return (der1, der2)
else
return nothing
end
end
function plat_av(obs::Vector{uwreal}, plat::Vector{Int64})
uwerr.(obs)
w = 1 ./ err.(obs)[plat[1]:plat[2]].^2
......@@ -127,7 +229,7 @@ function lin_fit(x::Vector{<:Real}, y::Vector{uwreal})
print("\n Fit parameter: ", i, ": ")
details(fitp[i])
end
println("Chisq / chiexp: ", chisq(par, y), " / ", csqexp, " (dof: ", x[end]-length(par),")")
println("Chisq / chiexp: ", chisq(par, y), " / ", csqexp, " (dof: ", length(x)-length(par),")")
return (fitp, csqexp)
end
......@@ -155,6 +257,7 @@ The method return an array upar with the best fit parameters with their errors.
fit_routine(model, ydata, param=3)
"""
function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing )
uwerr.(ydata)
yval = value.(ydata)
yer = err.(ydata)
chisq = gen_chisq(model, xdata, yer)
......@@ -172,7 +275,6 @@ function gen_chisq(f::Function, x::Array{<:Real}, err::Vector{Float64})
chisq(par, dat) = sum((dat .- f(x,par)).^2 ./err.^2)
return chisq
end
#TODO: add combined fits
#=
using LsqFit
@doc raw"""
......
src/juobs_types.jl
View file @ 54e7b8a7
......@@ -127,8 +127,7 @@ mutable struct Corr
h = getfield.(b, :header)
for s in sym
if !all(getfield.(h, s) .== getfield(h[1], s))
println("Corr: Parameter mismatch")
return nothing
error("Corr: Parameter mismatch")
end
end
mu = [h[1].mu1, h[1].mu2]
......
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