read_ms includes dtr != 1

src/juobs_reader.jl

@@ -259,20 +259,20 @@ function read_md(path::String)
 end
 
 @doc raw"""
-read_ms(path::String)
+read_ms(path::String; id::Union{Int64, Nothing}=nothing, dtr::Int64=1)  
 
-Reads openQCD ms dat files at a given path. This method return: 
+Reads openQCD ms dat files at a given path. This method return YData: 
     t(t): flow time values
-    Wsl(icfg, x0, t): the time-slice sums of the densities of the Wilson plaquette action
     Ysl(icfg, x0, t): the time-slice sums of the densities of the Yang-Mills action 
-    Qsl(icfg, x0, t): the time-slice sums of the densities of the topological charge
+    vtr: vector that contains trajectory number
+    id: ensmble id
 
 Examples:
 ```@example    
-t, W, Y, Q = read_ms(path)
+Y = read_ms(path)
 ```
 """
-function read_ms(path::String; id::Union{Int64, Nothing}=nothing)  
+function read_ms(path::String; id::Union{Int64, Nothing}=nothing, dtr::Int64=1)  
     if isnothing(id)
         bname = basename(path)
         m = findfirst(r"[A-Z][0-9]{3}r[0-9]{3}", bname)
@@ -290,27 +290,38 @@ function read_ms(path::String; id::Union{Int64, Nothing}=nothing)
 
     ntr = Int32((fsize - 3*4 - 8) / datsize)
 
-    vntr = Vector{Int32}(undef, ntr)
+    if mod(ntr, dtr) != 0
+        error("ntr / dtr must be exact")
+    end
+
+    vntr = Vector{Int32}(undef, div(ntr, dtr))
     
     # x0, t, cfg 
-    Wsl = Array{Float64}(undef, ntr, tvals, nn + 1) 
-    Ysl = Array{Float64}(undef, ntr, tvals, nn + 1)
-    Qsl = Array{Float64}(undef, ntr, tvals, nn + 1)
+    Wsl = Array{Float64}(undef, div(ntr, dtr), tvals, nn + 1) 
+    Ysl = Array{Float64}(undef, div(ntr, dtr), tvals, nn + 1)
+    Qsl = Array{Float64}(undef, div(ntr, dtr), tvals, nn + 1)
 
+    k = 0
     for itr = 1:ntr
-        vntr[itr] = read(data, Int32)
+        tmp = read(data, Int32)
+        if mod(itr, dtr) == 0
+            k += 1
+            vntr[k] = tmp
+        end
+
         for iobs = 1:3
             for inn = 0:nn
-                tmp = Vector{Float64}(undef, tvals)
-                read!(data, tmp)
-                if iobs == 1
-                    Wsl[itr, :, inn + 1] = tmp
-                elseif iobs == 2
-                    Ysl[itr, :, inn + 1] = tmp
-                elseif iobs == 3
-                    Qsl[itr, :, inn + 1] = tmp
+                tmp2 = Vector{Float64}(undef, tvals)
+                read!(data, tmp2)
+                if mod(itr, dtr) == 0
+                    if iobs == 1
+                        Wsl[k, :, inn + 1] = tmp2
+                    elseif iobs == 2
+                        Ysl[k, :, inn + 1] = tmp2
+                    elseif iobs == 3
+                        Qsl[k, :, inn + 1] = tmp2
+                    end
                 end
-
             end
         end
     end