update changes

src/juobs.jl

@@ -7,9 +7,9 @@ include("juobs_reader.jl")
 include("juobs_tools.jl")
 include("juobs_obs.jl")
 
-export read_mesons, read_ms1
+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 
+export meff, dec_const_pcvc, comp_t0
 
 end # module

src/juobs_obs.jl

@@ -93,3 +93,191 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu
 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)
+end
+#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)
+        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_reader.jl

@@ -213,4 +213,70 @@ function read_ms1(path::String; v::String="1.2")
     [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
+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 plaque
+tte action
+Ysl(x0, t, icfg): the time-slice sums of the densities of the Yang-Mills ac
+tion 
+Qsl(x0, t, icfg): the time-slice sums of the densities of the topological c
+harge
+
+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)
+
+    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_types.jl

@@ -138,6 +138,14 @@ mutable struct Corr
     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")