Initial commit

Project.toml

+name = "juobs"
+uuid = "3399d769-80ea-4938-a2f8-b7fb0ead32c9"
+authors = ["Javier <j.ugarrio@gmail.com>"]
+version = "0.1.0"
+
+[deps]
+ADerrors = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
+BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

README.md

+# juobs
+
+Analysis code (Julia)
\ No newline at end of file

src/juobs.jl

+module juobs
+using ADerrors, Plots, Statistics
+
+include("juobs_types.jl")
+include("juobs_reader.jl")
+include("juobs_tools.jl")
+
+export read_mesons, read_ms1
+export meff, plat_av, apply_rw, corr_obs
+
+end # module

src/juobs_reader.jl

+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
+
+
+
+function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int, 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
+
+
+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
\ No newline at end of file

src/juobs_tools.jl

+function plat_av(obs::Vector{uwreal}, plat::Vector{Int64})
+    uwerr.(obs)
+    w = 1 ./ err.(obs)[plat[1]:plat[2]].^2
+    av = sum(w .* obs[plat[1]:plat[2]]) / sum(w)
+    return av 
+end
+
+function meff(obs::Vector{uwreal}, plat::Vector{Int64})
+    dim = length(obs)
+    aux = log.(obs[1:dim-1] ./ obs[2:dim])
+    mass = plat_av(aux, plat)
+    uwerr(mass)
+    return mass
+end
+
+function apply_rw(cdata::CData, W::Array{Float64})
+    W1 = W[1, :]
+    W2 = W[2, :]
+    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
+end
+
+function apply_rw(data::Array{Float64}, W::Array{Float64})
+    W1 = W[1, :]
+    W2 = W[2, :]
+
+    data_r = data .* W1 .* W2 / mean(W1 .* W2)
+    return data_r
+end
+
+function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64}, Nothing}=nothing)
+
+    real ? data = cdata.re_data : data = cdata.im_data
+    isnothing(rw) ? data_r = data : data_r = apply_rw(data, rw)
+
+    nt = size(data)[2]
+    obs = Vector{uwreal}(undef, nt)
+    [obs[x0] = uwreal(data_r[:, x0], cdata.id) for x0 = 1:nt]
+    return Corr(obs, cdata)
+end
+
+#function corr_obs for R != 1
+#TODO: vcfg with gaps
+function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing)
+    nr = length(cdata)
+    id = getfield.(cdata, Symbol("id"))
+    vcfg = getfield.(cdata, Symbol("vcfg"))
+    replica = Int64.(maximum.(vcfg))
+    
+    if !all(id .== id[1])
+        println("IDs are not equal")
+        return nothing
+    end
+
+    real ? data = getfield.(cdata, Symbol("re_data")) : data = getfield.(cdata, Symbol("im_data"))
+    isnothing(rw) ? data_r = data : data_r = apply_rw.(data, rw)
+    
+
+    tmp = data_r[1]
+    [tmp = cat(tmp, data_r[k], dims=1) for k = 2:nr]
+
+    nt = size(data[1])[2]
+    obs = Vector{uwreal}(undef, nt)
+
+    [obs[x0] = uwreal(tmp[:, x0], id[1], replica) for x0 = 1:nt]
+
+    return Corr(obs, cdata)
+end    
\ No newline at end of file

src/juobs_types.jl

+ ##############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::Int
+    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}
+        function Corr(a::Vector{uwreal}, b::CData)
+            h = getfield(b, Symbol("header"))
+            mu = [h.mu1, h.mu2]
+            gamma = [gamma_name[h.type1+1], gamma_name[h.type2+1]]
+            return new(a, mu, gamma)
+        end 
+        function Corr(a::Vector{uwreal}, b::Vector{CData})
+            sym = Symbol.(["mu1", "mu2", "type1", "type2"])
+            h = getfield.(b, Symbol("header"))
+            for s in sym
+                if !all(getfield.(h, s) .== getfield(h[1], s))
+                    println("Corr: Parameter missmatch")
+                    return nothing
+                end
+            end
+            mu = [h[1].mu1, h[1].mu2]
+            gamma = [gamma_name[h[1].type1+1], gamma_name[h[1].type2+1]]
+            return new(a, mu, gamma)
+        end 
+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!=Symbol("type1") && k!=Symbol("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