comp_t0 added

read_ms now returns YData
comp_t0 computes t0 R!=1 + rwf

src/juobs.jl

@@ -10,6 +10,6 @@ include("juobs_obs.jl")
 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

@@ -98,3 +98,135 @@ 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)
+
+    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)
+        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

@@ -229,7 +229,12 @@ Examples:
 t, W, Y, Q = read_ms(path)
 ```
 """
-function read_ms(path::String)  
+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)
@@ -240,27 +245,27 @@ function read_ms(path::String)
     fsize = filesize(path)
     datsize=4 + 3*8*(nn + 1) * tvals # measurement size of each cnfg
 
-    ncfg = Int32((fsize - 3*4 - 8) / datsize)
+    ntr = Int32((fsize - 3*4 - 8) / datsize)
 
-    vcfg = Vector{Int32}(undef, ncfg)
+    vntr = Vector{Int32}(undef, ntr)
     
     # x0, t, cfg 
-    Wsl = Array{Float64}(undef, tvals, nn + 1, ncfg) 
-    Ysl = Array{Float64}(undef, tvals, nn + 1, ncfg)
-    Qsl = Array{Float64}(undef, tvals, nn + 1, ncfg)
+    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 icfg = 1:ncfg
-        vcfg[icfg] = read(data, Int32)
+    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, icfg] = tmp
+                    Wsl[:, inn + 1, itr] = tmp
                 elseif iobs == 2
-                    Ysl[:, inn + 1, icfg] = tmp
+                    Ysl[:, inn + 1, itr] = tmp
                 elseif iobs == 3
-                    Qsl[:, inn + 1, icfg] = tmp
+                    Qsl[:, inn + 1, itr] = tmp
                 end
 
             end
@@ -268,5 +273,5 @@ function read_ms(path::String)
     end
     close(data)
     t = Float64.(0:nn) .* dn .* eps
-    return (t, Wsl, Ysl, Qsl)
+    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")