Merge branch 'javier'

Project.toml

@@ -8,6 +8,6 @@ ADerrors = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
 BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

src/juobs.jl

 module juobs
-using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit
+using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, Optim
 import Statistics: mean
 
 include("juobs_types.jl")
@@ -8,7 +8,7 @@ include("juobs_reader.jl")
 include("juobs_tools.jl")
 include("juobs_obs.jl")
 
-export read_mesons, read_ms1, read_ms, read_md
+export read_mesons, read_ms1, read_ms, read_md, truncate_data!
 export get_matrix, uwgevp_tot, energies, uwdot
 export corr_obs, md_sea, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
 export meff, dec_const_pcvc, comp_t0

src/juobs_obs.jl

@@ -14,12 +14,13 @@ corr_pp = corr_obs.(data)
 m = meff(corr_pp[1], [50, 60], pl=false)
 ```
 """
-function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, mu::Union{Vector{Float64}, Nothing}=nothing)     
+function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bool=false, mu::Union{Vector{Float64}, Nothing}=nothing,
+    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)     
     dim = length(corr)                                                                                    
     aux = 0.5 .* log.((corr[2:dim-2] ./ corr[3:dim-1]).^2)     
-    mass = plat_av(aux, plat)                                                                            
-    uwerr(mass)                       
+    mass = plat_av(aux, plat, wpm)                                                                            
     if pl == true
+        isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm)                       
         x = 1:length(aux)
         y = value.(aux)
         dy = err.(aux)
@@ -65,7 +66,8 @@ m = meff(corr_pp[1], [50, 60], pl=false)
 f = dec_const_pcvc(corr_pp[1], [50, 60], m, pl=false)
 ```
 """
-function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false)
+function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu::Vector{Float64}, y0::Int64 ; pl::Bool=true, data::Bool=false,
+    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
     """
     compute the decay constant when the source is far from the boundaries
     """
@@ -74,14 +76,21 @@ function dec_const_pcvc(corr::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, mu
                                                                  
     aux = exp.((collect(1:dim) .- y0 ) .* [m for k in 1:dim])
     R = ((aux .* corr_pp).^2).^0.25
-    R_av = plat_av(R, plat)
+    R_av = plat_av(R, plat, wpm)
     f = sqrt(2) * (mu[1] + mu[2]) *R_av / m^1.5
-    uwerr(f)
     if pl == true
-        uwerr(R_av)
+        if isnothing(wpm)
+            uwerr(f)
+            uwerr(R_av)
+            uwerr.(R)
+        else
+            uwerr(f, wpm)
+            uwerr(R_av, wpm)
+            uwerr.(R, wpm)
+        end
+    
         v = value(R_av)
         e = err(R_av)                                
-        uwerr.(R)
         figure()
         fill_between(plat[1]:plat[2], v-e, v+e, color="green", alpha=0.75)                                            
         errorbar(1:length(R), value.(R), err.(R), fmt="x", color="black")
@@ -138,10 +147,31 @@ t0_r = comp_t0(Y, [38, 58], L=32, rw=[rw1, rw2], pl=true)
 
 ```
 """
-function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2)
+function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, 
+    rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, ws::ADerrors.wspace=ADerrors.wsg,
+    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+
     Ysl = Y.Ysl
     t = Y.t
     id = Y.id
+    replica = size.([Ysl], 1)
+
+    #Truncation
+    n_ws = findfirst(x-> x == id, ws.map_nob)
+    if !isnothing(n_ws)
+        ivrep_ws = ws.fluc[n_ws].ivrep
+
+        if length(ivrep_ws) != 1
+            error("Different number of replicas")
+        end
+
+        if replica[1] > ivrep_ws[1]
+            println("Automatic truncation in Ysl ", ivrep_ws[1], " / ", replica[1], ". R = 1")
+            Ysl = Ysl[1:ivrep_ws[1], :, :]
+        elseif replica[1] < ivrep_ws[1]
+            error("Automatic truncation failed. R = 1\nTry using truncate_data!")
+        end
+    end
 
     Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw)
     xmax = size(Ysl, 2)
@@ -157,7 +187,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Un
         end
     end
     x = t[nt0-dt0:nt0+dt0]
-    t2E = [plat_av(Y_aux[:, j], plat) for j=1:2*dt0+1] .* x.^2 / L^3
+    t2E = [plat_av(Y_aux[:, j], plat, wpm) for j=1:2*dt0+1] .* x.^2 / L^3
     
     model(x, p) = get_model(x, p, npol)
 
@@ -165,8 +195,14 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Un
     fmin(x, p) = model(x, p) .- 0.3
     t0 = root_error(fmin, t[nt0], par)
     if pl
-        uwerr(t0)
-        uwerr.(t2E)
+        if isnothing(wpm)
+            uwerr(t0)
+            uwerr.(t2E)
+        else
+            uwerr(t0, wpm)
+            uwerr.(t2E, wpm)
+        end
+
         v = value.(t2E)
         e = err.(t2E)
 
@@ -186,23 +222,43 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Un
         xlabel(L"$x_0/a$")
         title(string(L"$t/a^2 = $", t[nt0]))
         display(gcf())
-        end
+    end
     return t0
 end
 
-function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2)
+function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, 
+    rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, ws::ADerrors.wspace=ADerrors.wsg,
+    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+
     nr = length(Y)
     Ysl = getfield.(Y, :Ysl)
     t = getfield.(Y, :t)
     t = t[1]
     id = getfield.(Y, :id)
-    vtr = getfield.(Y, :vtr)
-    replica = length.(vtr)
+    replica = size.(Ysl, 1)
 
     if !all(id .== id[1])
         error("IDs are not equal")
     end
+    #Truncation
+    n_ws = findfirst(x-> x == id[1], ws.map_nob)
+    if !isnothing(n_ws)
+        ivrep_ws = ws.fluc[n_ws].ivrep
+
+        if length(replica) != length(ivrep_ws)
+            error("Different number of replicas")
+        end
 
+        for k = 1:length(replica)
+            if replica[k] > ivrep_ws[k]
+                println("Automatic truncation in Ysl ", ivrep_ws[k], " / ", replica[k], ". R = ", k)
+                Ysl[k] = Ysl[k][1:ivrep_ws[k], :, :]
+            elseif replica[k] < ivrep_ws[k]
+                error("Automatic truncation failed. R = ", replica[k], "\nTry using truncate_data!")
+            end
+        end
+    end
+    
     Ysl = isnothing(rw) ? Ysl : apply_rw(Ysl, rw)
     tmp = Ysl[1]
     [tmp = cat(tmp, Ysl[k], dims=1) for k = 2:nr]
@@ -219,7 +275,7 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
         end
     end
     x = t[nt0-dt0:nt0+dt0]
-    t2E = [plat_av(Y_aux[:, j], plat) for j=1:2*dt0+1] .* x.^2 / L^3
+    t2E = [plat_av(Y_aux[:, j], plat, wpm) for j=1:2*dt0+1] .* x.^2 / L^3
     
     model(x, p) = get_model(x, p, npol)
 
@@ -227,8 +283,14 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
     fmin(x, p) = model(x, p) .- 0.3
     t0 = root_error(fmin, t[nt0], par)
     if pl
-        uwerr(t0)
-        uwerr.(t2E)
+        if isnothing(wpm)
+            uwerr(t0)
+            uwerr.(t2E)
+        else
+            uwerr(t0, wpm)
+            uwerr.(t2E, wpm)
+        end
+
         v = value.(t2E)
         e = err.(t2E)
 

src/juobs_reader.jl

@@ -71,6 +71,7 @@ end
 
 @doc raw"""
 read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
+read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
 
 This faction read a mesons dat file at a given path and returns a vector of CData structures for different masses and Dirac structures.
 Dirac structures g1 and/or g2 can be passed as string arguments in order to filter correaltors.
@@ -83,6 +84,7 @@ read_mesons(path, "G5")
 read_mesons(path, nothing, "G5")
 read_mesons(path, "G5", "G5")
 read_mesons(path, "G5", "G5", id=1)
+read_mesons([path1, path2], "G5", "G5")
 ```
 """
 function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
@@ -156,6 +158,21 @@ function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union
     return res
 end
 
+function read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
+    res = read_mesons.(path, g1, g2, id=id)
+    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)
@@ -242,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)
@@ -273,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
@@ -301,3 +329,60 @@ function read_ms(path::String; id::Union{Int64, Nothing}=nothing)
     t = Float64.(0:nn) .* dn .* eps
     return YData(vntr, t, Ysl, id)
 end
+
+@doc raw"""
+truncate_data!(data::YData, nc::Int64)
+
+truncate_data!(data::Vector{YData}, nc::Vector{Int64})
+
+truncate_data!(data::Vector{CData}, nc::Int64)
+
+truncate_data!(data::Vector{Vector{CData}}, nc::Vector{Int64})
+
+Truncates the output of read_mesons and read_ms taking the first nc configurations.
+
+Examples:
+```@example    
+#Single replica
+dat = read_mesons(path, "G5", "G5")
+Y = read_ms(path)
+truncate_data!(dat, nc)
+truncate_data!(Y, nc)
+
+#Two replicas
+dat = read_mesons([path1, path2], "G5", "G5")
+Y = read_ms.([path1_ms, path2_ms])
+truncate_data!(dat, [nc1, nc2])
+truncate_data!(Y, [nc1, nc2])
+```
+"""
+function truncate_data!(data::YData, nc::Int64)
+    data.vtr = data.vtr[1:nc]
+    data.Ysl = data.Ysl[1:nc, :, :]
+    return nothing
+end
+function truncate_data!(data::Vector{YData}, nc::Vector{Int64})
+    truncate_data!.(data, nc)
+    return nothing
+end
+function truncate_data!(data::Vector{CData}, nc::Int64)
+    N = length(data)
+    for k = 1:N
+        data[k].vcfg = data[k].vcfg[1:nc]
+        data[k].re_data = data[k].re_data[1:nc, :]
+        data[k].im_data = data[k].im_data[1:nc, :]
+    end
+    return nothing
+end
+function truncate_data!(data::Vector{Vector{CData}}, nc::Vector{Int64})
+    N = length(data)
+    R = length(data[1])
+    for k = 1:N
+        for r = 1:R
+            data[k][r].vcfg = data[k][r].vcfg[1:nc[r]]
+            data[k][r].re_data = data[k][r].re_data[1:nc[r], :]
+            data[k][r].im_data = data[k][r].im_data[1:nc[r], :]
+        end
+    end
+    return nothing
+end
\ No newline at end of file