Added methods for read_meson and read_meson_correnction that take a...

Added methods for read_meson and read_meson_correnction that take a Vector{Tuple{String,String}} as Gamma structure. This method allow to read in one sweep all correlator for any specific gamma structure combination. It is, in a way, complementary to what was already available when one of the two gamma structure was not specified

src/juobs_reader.jl

@@ -196,6 +196,114 @@ function read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g
     return cdata
 end
 
+@doc raw"""
+    read_mesons(path::String, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,  nnoise_trunc::Union{Int64, Nothing}=nothing)
+    
+    read_mesons(path::Vector{String}, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,  nnoise_trunc::Union{Int64, Nothing}=nothing)
+
+
+Read a meson.dat file and return a `Vector{Cdata}` with for different masses and the  Dirac structures specified in `gamma`.
+The Dirac structures are specified into `gamma` in a vector of Tuple as `(Gsrc,Gsnk)`. 
+
+This method is best used when muliple, but specific, Dirac structure are requested. It is up to the user reorganized the output in the preferred way.
+"""
+function read_mesons(path::String, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,
+  nnoise_trunc::Union{Int64, Nothing}=nothing)    
+  
+  type = Vector{Tuple{Int64,Int64}}(undef,length(gamma)) 
+  
+  for i in eachindex(gamma)
+    G1,G2 = gamma[i]
+    type[i] = (findfirst(x-> x==G1, juobs.gamma_name) - 1, findfirst(x-> x==G2, juobs.gamma_name) - 1)
+  end
+
+     
+  if isnothing(id)
+      bname = basename(path)
+      m = findfirst(r"[A-Z][0-9]{3}r[0-9]{3}", bname)
+      id = bname[m[1: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,x.type2) in type, c_header)
+      
+  seek(data, g_header.hsize + sum(c.hsize for c in c_header))
+
+  res = Array{juobs.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[1:nnoise_trunc, :], 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[:, 1:nnoise_trunc, :], 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 eachindex(corr_match)
+      res[c] = juobs.CData(c_header[corr_match[c]], vcfg, data_re[c, :, :], data_im[c, :, :], id)
+  end
+  close(data)
+
+  return res
+end
+
+function read_mesons(path::Vector{String}, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,
+  nnoise_trunc::Union{Int64, Nothing}=nothing)
+  res = [read_mesons(p, gamma, id=id, legacy=legacy, nnoise_trunc=nnoise_trunc) for p in path]
+  nrep = length(res)
+  ncorr = length(res[1])
+
+  cdata = Vector{Vector{juobs.CData}}(undef, ncorr)
+  for icorr = 1:ncorr
+      cdata[icorr] = Vector{juobs.CData}(undef, nrep) 
+      for r = 1:nrep
+          cdata[icorr][r] = res[r][icorr]
+      end
+  end
+  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
@@ -273,6 +381,18 @@ function read_mesons_correction(path::String, g1::Union{String, Nothing}=nothing
     return res
 end
 
+
+@doc raw"""
+    read_mesons_correction(path::String, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,  nnoise_trunc::Union{Int64, Nothing}=nothing)
+
+    read_mesons_correction(path::Vector{String}, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,  nnoise_trunc::Union{Int64, Nothing}=nothing)
+
+
+Read a meson.dat file with correction data and return a `Vector{Cdata}` with for different masses and the  Dirac structures specified in `gamma`. (See also [`read_meson`](@ref))
+The Dirac structures are specified into `gamma` in a vector of Tuple as `(Gsrc,Gsnk)`. 
+
+This method is best used when muliple, but specific, Dirac structure are requested. It is up to the user reorganized the output in the preferred way.
+"""
 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)
@@ -289,6 +409,105 @@ function read_mesons_correction(path::Vector{String}, g1::Union{String, Nothing}
     return cdata
 end
 
+function dev_read_mesons_correction(path::String, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,
+  nnoise_trunc::Union{Int64, Nothing}=nothing)
+  
+  type = Vector{Tuple{Int64,Int64}}(undef,length(gamma)) 
+  
+  for i in eachindex(gamma)
+    G1,G2 = gamma[i]
+    type[i] = (findfirst(x-> x==G1, juobs.gamma_name) - 1, findfirst(x-> x==G2, juobs.gamma_name) - 1)
+  end
+
+  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 = juobs.read_GHeader(path)
+  c_header = juobs.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,x.type2) in type, c_header)
+
+      
+  seek(data, g_header.hsize + sum(c.hsize for c in c_header))
+
+  res = Array{juobs.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] = juobs.CData(c_header[corr_match[c]], vcfg, data_re[c, :, :], data_im[c, :, :], id)
+  end
+  close(data)
+
+  return res
+end
+
+function dev_read_mesons_correction(path::Vector{String}, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,
+nnoise_trunc::Union{Int64, Nothing}=nothing)
+res = [dev_read_mesons_correction(p, gamma, id=id, legacy=legacy, nnoise_trunc=nnoise_trunc) for p in path]
+nrep = length(res)
+ncorr = length(res[1])
+
+cdata = Vector{Vector{juobs.CData}}(undef, ncorr)
+for icorr = 1:ncorr
+    cdata[icorr] = Vector{juobs.CData}(undef, nrep) 
+    for r = 1:nrep
+        cdata[icorr][r] = res[r][icorr]
+    end
+end
+return cdata
+end
+
 @doc raw"""
     read_mesons_chunks(paths::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)