Commit 4e0996ee authored by Alberto Ramos's avatar Alberto Ramos

Added main code for basica analysis.

Added also analysis of run SU3_PBC_bt5.96_L24x24x24x48
parent b102f136
## Results of analysis of file: data/SU3_PBC_bt5.96_L24x24x24x48/SU3_PBC_bt5.96_L24x24x24x48.bdio
- Number of thermalization steps: 1000
- Number of measurements: 10000
- flow measurements for times: 0.0 - 3.0
##
# Observable <exp(-dH)>
- value: 0.9999944130766438 +/- 1.1212491302990926e-5
- tauint: 0.20025032769951406 +/- 0.010131855128562722
# Observable <P>
- value: 1.7674737643415093 +/- 1.3696790349384248e-5
- tauint: 3.6602965209536578 +/- 0.5559307314300315
# Observable t0 plaquette
- value: 2.5376415250592386 +/- 0.0041083924549615725
- tauint: 2.7303027666101665 +/- 0.9479279130978608
# Observable w0 plaquette
- value: 1.597321379304086 +/- 0.0012407878980924313
- tauint: 2.7203815563944858 +/- 0.9445611237721703
# Observable t0 clover
- value: 2.771930894344794 +/- 0.004275783803034508
- tauint: 2.8329858832313772 +/- 0.9989400205301594
# Observable w0 clover
- value: 1.6679427847316557 +/- 0.0012743822359628887
- tauint: 2.8327189004176665 +/- 0.9988480244606983
###
### "THE BEER-WARE LICENSE":
### Alberto Ramos wrote this file. As long as you retain this
### notice you can do whatever you want with this stuff. If we meet some
### day, and you think this stuff is worth it, you can buy me a beer in
### return. <alberto.ramos@cern.ch>
###
### file: ana.jl
### created: Tue Nov 9 17:28:04 2021
###
using Printf, Plots, BDIO, ADerrors, ArgParse
function parse_commandline()
s = ArgParseSettings()
@add_arg_table s begin
"-i"
help = "input file"
required = true
arg_type = String
end
@add_arg_table s begin
"-o"
help = "output directory"
required = true
arg_type = String
end
return parse_args(s)
end
function get_msm_info(fb)
nthm = 0
nmsm = 0
nflw = 0
flw_ns = 0
eps = 0.0
BDIO_seek!(fb, 0)
while BDIO_seek!(fb)
if BDIO_get_uinfo(fb) == 2
nthm = round(Int64, BDIO_get_len(fb)/8)
end
if BDIO_get_uinfo(fb) == 3
nmsm = round(Int64, BDIO_get_len(fb)/8)
end
if BDIO_get_uinfo(fb) == 12
nflw = round(Int64, BDIO_get_len(fb)/8-2)
foo = zeros(2)
BDIO_read(fb, foo)
flw_ns = round(Int64, foo[2])
eps = foo[1]
end
end
return nthm, nmsm, round(Int64, nflw/(flw_ns+1)), eps, flw_ns
end
function read_record(dt::Vector, fb, uinfo)
BDIO_seek!(fb, 0)
while BDIO_seek!(fb)
if BDIO_get_uinfo(fb) == uinfo
BDIO_read(fb, dt)
break
end
end
return nothing
end
function read_record(dt::Array{Float64, 2}, fb, uinfo)
BDIO_seek!(fb, 0)
while BDIO_seek!(fb)
if BDIO_get_uinfo(fb) == uinfo
foo = Vector{Float64}(undef, size(dt,1))
for i in 1:size(dt,2)
BDIO_read(fb, foo)
for j in 1:length(foo)
dt[j,i] = foo[j]
end
end
break
end
end
return nothing
end
function ana_record(fb, nmsm, uinfo)
dt = Vector{Float64}(undef, nmsm)
read_record(dt, fb, uinfo)
if uinfo == 3
obs = uwreal(exp.(-dt), "Ensemble")
else
obs = uwreal(dt, "Ensemble")
end
uwerr(obs)
return obs
end
function ana_flow_record(fb, nflw, flw_ns, eps, uinfo, tsqet::Union{Nothing, Float64} = nothing)
dt = Array{Float64, 2}(undef, flw_ns+1, nflw)
read_record(dt, fb, uinfo)
for i in 1:size(dt, 1)
t = (i-1)*eps
for j in 1:size(dt,2)
dt[i,j] = dt[i,j]*t^2
end
end
obs = Vector{uwreal}(undef, flw_ns+1)
for i in 1:flw_ns+1
obs[i] = uwreal(dt[i,:], "Flow")
end
uwerr.(obs)
if tsqet == nothing
return obs
else
# Now determine t0 scale
id = 0
d = 100.0
for i in 1:length(obs)
if abs(value(obs[i])-tsqet) < d
d = abs(value(obs[i])-tsqet)
id = i
end
end
ta = (id-2)*eps
tb = (id-1)*eps
B = (obs[id-1]-obs[id])/(ta-tb)
A = obs[id-1] - B*ta
t0 = (tsqet-A)/B
# w0 scale
dobs = Vector{uwreal}(undef, length(obs)-2)
for i in 2:length(obs)-1
dobs[i-1] = i*(obs[i+1]-obs[i-1])/2
end
id = 0
d = 100.0
for i in 1:length(dobs)
if abs(value(dobs[i])-tsqet) < d
d = abs(value(dobs[i])-tsqet)
id = i
end
end
ta = (id-1)*eps
tb = (id)*eps
B = (obs[id-1]-obs[id])/(ta-tb)
A = obs[id-1] - B*ta
w0 = sqrt((tsqet-A)/B)
return obs, t0, w0
end
return nothing
end
function make_plot(obs)
ens = ensembles(obs)
l = grid(2,1, heights=[0.3,0.7])
pl1 = plot(mchist(obs, ens[1]) .+ value(obs), xlabel="MC time",
ylabel="", label="")
r = rho(obs, ens[1])
dr = drho(obs, ens[1])
w = max(20, 2*window(obs, ens[1]))
pl2 = plot(r[1:w], yerr=dr[1:w], xlabel="W", ylabel="AC function",
seriestype=:scatter, label="")
pl2 = vline!([window(obs, ens[1])], label="")
pl = plot(pl1, pl2, layout=l)
return pl
end
parsed_args = parse_commandline()
fname = parsed_args["i"]
dout = parsed_args["o"]
fb = BDIO_open(fname, "r")
nthm, nmsm, nflw, eps, flw_ns = get_msm_info(fb)
fout = open(dout*"/results.txt", "w")
println(fout, "\n## Results of analysis of file: "*fname)
println(fout, " - Number of thermalization steps: ", nthm)
println(fout, " - Number of measurements: ", nmsm)
println(fout, " - flow measurements for times: ", 0.0, " - ", flw_ns*eps)
println(fout, "## ")
println(fout, " # Observable <exp(-dH)>")
obs = ana_record(fb, nmsm, 3)
println(fout, " - value: ", obs)
println(fout, " - tauint: ", taui(obs, "Ensemble"), " +/- ", dtaui(obs, "Ensemble"))
pl = make_plot(obs)
savefig(pl, dout*"/IamOne.pdf")
println(fout, " # Observable <P>")
obs = ana_record(fb, nmsm, 5)
println(fout, " - value: ", obs)
println(fout, " - tauint: ", taui(obs, "Ensemble"), " +/- ", dtaui(obs, "Ensemble"))
pl = make_plot(obs)
savefig(pl, dout*"/plaq.pdf")
println(fout, " # Observable t0 plaquette")
obs, t0, w0 = ana_flow_record(fb, nflw, flw_ns, eps, 12, 0.3)
uwerr(t0)
println(fout, " - value: ", t0)
println(fout, " - tauint: ", taui(t0, "Flow"), " +/- ", dtaui(t0, "Flow"))
pl = make_plot(t0)
savefig(pl, dout*"/t0_pl.pdf")
println(fout, " # Observable w0 plaquette")
uwerr(w0)
println(fout, " - value: ", w0)
println(fout, " - tauint: ", taui(w0, "Flow"), " +/- ", dtaui(w0, "Flow"))
pl = make_plot(w0)
savefig(pl, dout*"/w0_pl.pdf")
println(fout, " # Observable t0 clover")
obs, t0, w0 = ana_flow_record(fb, nflw, flw_ns, eps, 13, 0.3)
uwerr(t0)
println(fout, " - value: ", t0)
println(fout, " - tauint: ", taui(t0, "Flow"), " +/- ", dtaui(t0, "Flow"))
pl = make_plot(t0)
savefig(pl, dout*"/t0_cl.pdf")
println(fout, " # Observable w0 clover")
uwerr(w0)
println(fout, " - value: ", w0)
println(fout, " - tauint: ", taui(w0, "Flow"), " +/- ", dtaui(w0, "Flow"))
pl = make_plot(w0)
savefig(pl, dout*"/w0_cl.pdf")
BDIO_close!(fb)
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