read_mesons_correction + corr_obs_TSM added

functions for TSM runs.
read_mesons_correction -> read correction runs (exact + sloppy solutions)
corr_obs_TSM -> create Corr with <O_sloppy> + <O_correction>

src/juobs.jl

@@ -8,9 +8,9 @@ include("juobs_reader.jl")
 include("juobs_tools.jl")
 include("juobs_obs.jl")
 
-export read_mesons, read_ms1, read_ms, read_md, truncate_data!
+export read_mesons, read_mesons_correction, read_ms1, read_ms, read_md, truncate_data!
 export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs
-export corr_obs, corr_sym, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
+export corr_obs, corr_obs_TSM, corr_sym, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
 export meff, mpcac, dec_const, dec_const_pcvc, comp_t0
 
 end # module

src/juobs_reader.jl

@@ -191,6 +191,99 @@ function read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g
     return cdata
 end
 
+function read_mesons_correction(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false,
+    nnoise_trunc::Union{Int64, Nothing}=nothing)
+    t1 = isnothing(g1) ? nothing : findfirst(x-> x==g1, gamma_name) - 1
+    t2 = isnothing(g2) ? nothing : 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 = bname[m[1:4]]
+        #id = parse(Int64, bname[m[2:4]])
+    end
+
+    data = open(path, "r")
+    g_header = read_GHeader(path)
+    c_header = read_CHeader(path, legacy=legacy)
+
+    ncorr = g_header.ncorr
+    tvals = g_header.tvals
+    nnoise = g_header.nnoise
+
+    nnoise_trunc = isnothing(nnoise_trunc) ? nnoise : min(nnoise, nnoise_trunc)
+
+    fsize = filesize(path)
+
+    datsize = 4 + sum(getfield.(c_header, :dsize)) * tvals * nnoise #data_size / ncnfg
+    ncfg = div(fsize - g_header.hsize - sum(getfield.(c_header, :hsize)), datsize) #(total size - header_size) / data_size
+
+    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, div(length(corr_match), 2)) # Modification: total length is divided by 2
+
+    data_re = zeros(div(length(corr_match), 2), ncfg, tvals) # Modification: total length is divided by 2
+    data_im = zeros(div(length(corr_match), 2), ncfg, tvals) # Modification: total length is divided by 2
+    vcfg = Array{Int32}(undef, ncfg)
+    
+    for icfg = 1:ncfg
+        vcfg[icfg] = read(data, Int32)
+        c = 1
+        sgn = +1 # sign it changes at ncorr / 2. O_exact - O_sloppy
+        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[1:nnoise_trunc, :], dims=1)
+
+                    data_re[c, icfg, :] = data_re[c, icfg, :] + sgn * 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[:, 1:nnoise_trunc, :], dims=2)
+                    data_re[c, icfg, :] = data_re[c, icfg, :] + sgn * tmp2[1, 1, :]
+                    data_im[c, icfg, :] = data_im[c, icfg, :] + sgn * tmp2[2, 1, :]
+
+                end
+                c += 1
+            else
+                seek(data, position(data)  + c_header[k].dsize*tvals*nnoise)
+            end
+            if k == div(ncorr, 2)
+                c = 1
+                sgn = -1
+            end
+        end
+    end
+    for c = 1:div(length(corr_match), 2)
+        res[c] = CData(c_header[corr_match[c]], vcfg, data_re[c, :, :], data_im[c, :, :], id)
+    end
+    close(data)
+
+    return res
+end
+
+function read_mesons_correction(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{String, Nothing}=nothing, legacy::Bool=false,
+    nnoise_trunc::Union{Int64, Nothing}=nothing)
+    res = read_mesons_correction.(path, g1, g2, id=id, legacy=legacy, nnoise_trunc=nnoise_trunc)
+    nrep = length(res)
+    ncorr = length(res[1])
+
+    cdata = Vector{Vector{CData}}(undef, ncorr)
+    for icorr = 1:ncorr
+        cdata[icorr] = Vector{CData}(undef, nrep) 
+        for r = 1:nrep
+            cdata[icorr][r] = res[r][icorr]
+        end
+    end
+    return cdata
+end
+
 function read_rw(path::String; v::String="1.2")
     data = open(path, "r")
     nrw = read(data, Int32)

src/juobs_tools.jl

@@ -112,7 +112,6 @@ function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, No
         return Corr(obs, cdata)
     end
 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, L::Int64=1, info::Bool=false)
@@ -151,6 +150,99 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
     else
         return Corr(obs, cdata)
     end
+end   
+
+function corr_obs_TSM(cdata1::CData, cdata2::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1, info::Bool=false)
+    if cdata1.id != cdata2.id
+        error("Error: cdata1 id != cdata2 id")
+    end
+    if cdata1.header != cdata2.header # Base.:(==) and Base.:(!=) are redifined in juobs_types.jl
+        error("Error: cdata1 header != cdata2 header")
+    end
+
+    data1 = real ? cdata1.re_data ./ L^3 : cdata1.im_data ./ L^3
+    data2 = real ? cdata2.re_data ./ L^3 : cdata2.im_data ./ L^3
+
+    nt = size(data1, 2)
+    id = cdata1.id
+
+    if isnothing(rw)
+        obs1 = [uwreal(data1[:, x0], id) for x0 = 1:nt]
+        obs2 = [uwreal(data2[:, x0], id) for x0 = 1:nt]
+    else
+        data1_r, W = apply_rw(data1, rw)
+        data2_r, W = apply_rw(data2, rw)
+
+        ow1 = [uwreal(data1_r[:, x0], id) for x0 = 1:nt]
+        ow2 = [uwreal(data2_r[:, x0], id) for x0 = 1:nt]
+
+        W_obs = uwreal(W, id)
+
+        obs1 = [ow1[x0] / W_obs for x0 = 1:nt]
+        obs2 = [ow2[x0] / W_obs for x0 = 1:nt]
+    end
+
+    if info && !isnothing(rw)
+        return (Corr(obs1 + obs2, cdata1), ow1, ow2, W_obs)
+    elseif info && isnothing(rw)
+        return (Corr(obs1 + obs2, cdata1), obs1, obs2)
+    else
+        return Corr(obs1 + obs2, cdata1)
+    end
+
+end
+function corr_obs_TSM(cdata1::Array{CData, 1}, cdata2::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, info::Bool=false)
+    if any(getfield.(cdata1, :id) .!= getfield.(cdata2, :id))
+        error("Error: cdata1 id != cdata2 id")
+    end
+    if any(getfield.(cdata1, :header) .!= getfield.(cdata2, :header)) # Base.:(==) and Base.:(!=) are redifined in juobs_types.jl
+        error("Error: cdata1 header != cdata2 header")
+    end
+    
+    id = getfield.(cdata1, :id)
+    vcfg = getfield.(cdata1, :vcfg)
+    replica = Int64.(maximum.(vcfg))
+    
+    if !all(id .== id[1])
+        error("IDs are not equal")
+    end
+
+    data1 = real ? getfield.(cdata1, :re_data) ./ L^3 : getfield.(cdata1, :im_data) ./ L^3
+    data2 = real ? getfield.(cdata2, :re_data) ./ L^3 : getfield.(cdata2, :im_data) ./ L^3
+
+    nt = size(data1[1], 2)
+
+    if isnothing(rw)
+        tmp1 = cat(data1..., dims=1)
+        tmp2 = cat(data2..., dims=1)
+
+        obs1 = [uwreal(tmp1[:, x0], id[1], replica) for x0 = 1:nt]
+        obs2 = [uwreal(tmp2[:, x0], id[1], replica) for x0 = 1:nt]
+
+    else
+        data1_r, W = apply_rw(data1, rw)
+        data2_r, W = apply_rw(data2, rw)
+
+        tmp1 = cat(data1_r..., dims=1)
+        tmp2 = cat(data2_r..., dims=1)
+        tmp_W = cat(W..., dims=1)
+
+        ow1 = [uwreal(tmp1[:, x0], id[1], replica) for x0 = 1:nt]
+        ow2 = [uwreal(tmp2[:, x0], id[1], replica) for x0 = 1:nt]
+
+        W_obs = uwreal(tmp_W, id[1], replica)
+
+        obs1 = [ow1[x0] / W_obs for x0 = 1:nt]
+        obs2 = [ow2[x0] / W_obs for x0 = 1:nt]
+
+    end
+    if info && !isnothing(rw)
+        return (Corr(obs1 + obs2, cdata), ow1, ow2, W_obs)
+    elseif info && isnothing(rw)
+        return (Corr(obs1 + obs2, cdata), obs1, obs2)
+    else
+        return Corr(obs1 + obs2, cdata)
+    end
 end    
 @doc raw"""
     corr_sym(corrL::Corr, corrR::Corr, parity::Int64=1)

src/juobs_types.jl

@@ -127,6 +127,15 @@ mutable struct CHeader
         return a
     end
 end
+function Base.:(==)(a::CHeader, b::CHeader)
+    for s in [:k1, :k2, :mu1, :mu2, :dp1, :dp2, :type1, :type2, :x0, :is_real, :theta1, :theta2]
+        if getfield(a, s) != getfield(b, s)
+            return false
+        end
+    end
+    return true
+end
+#Base.:(!=)(a::CHeader, b::CHeader) = !(a == b)
 
 mutable struct CData
     header::CHeader