Merge branch 'md_val'

src/juobs.jl

@@ -10,7 +10,7 @@ include("juobs_obs.jl")
 
 export read_mesons, read_ms1, read_ms, read_md, truncate_data!
 export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs
-export corr_obs, md_sea, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
+export corr_obs, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
 export meff, dec_const_pcvc, comp_t0
 
 end # module

src/juobs_reader.jl

-function read_global_header(path::String)
+function read_GHeader(path::String)
     data = open(path, "r")
     g_header = zeros(Int32, 4)
     read!(data, g_header)
@@ -7,77 +7,89 @@ function read_global_header(path::String)
     return a
 end
 
-function read_CHeader(path::String)
-    gh = read_global_header(path)
+function read_CHeader(path::String; legacy::Bool=false)
+    gh = read_GHeader(path)
     
     data = open(path, "r")
     seek(data, gh.hsize)
     
-    aux_f = zeros(Float64, 6)
-    aux_i = zeros(Int32, 4)
-    theta = zeros(Float64, 6)
     a = Vector{CHeader}(undef, gh.ncorr)
-    for k = 1:gh.ncorr
-        read!(data, aux_f)
-        read!(data, theta)
-
-        qs1 = read(data, Int32)
-        if qs1 != 0
-            qn1 = read(data, Int32)
-            qeps1 = read(data, Float64)
-            q1 = Sm(qs1, qn1, qeps1, 1)
-        else
-            q1 = Sm(qs1, 1)
-        end
+    if !legacy
+        aux_f = zeros(Float64, 6)
+        aux_i = zeros(Int32, 4)
+        theta = zeros(Float64, 6)
+
+        for k = 1:gh.ncorr
+            read!(data, aux_f)
+            read!(data, theta)
+
+            qs1 = read(data, Int32)
+            if qs1 != 0
+                qn1 = read(data, Int32)
+                qeps1 = read(data, Float64)
+                q1 = Sm(qs1, qn1, qeps1, 1)
+            else
+                q1 = Sm(qs1, 1)
+            end
 
-        qs2 = read(data, Int32)
-        if qs2 != 0
-            qn2 = read(data, Int32)
-            qeps2 = read(data, Float64)
-            q2 = Sm(qs2, qn2, qeps2, 1)
-        else
-            q2 = Sm(qs2, 1)
-        end
+            qs2 = read(data, Int32)
+            if qs2 != 0
+                qn2 = read(data, Int32)
+                qeps2 = read(data, Float64)
+                q2 = Sm(qs2, qn2, qeps2, 1)
+            else
+                q2 = Sm(qs2, 1)
+            end
 
 
-        gs1 = read(data, Int32)
-        if gs1 != 0 && gs1 != 3 && gs1 != 4
-            gn1 = read(data, Int32)
-            geps1 = read(data, Float64)
-            g1 = Sm(gs1, gn1, geps1, 2)
-        elseif gs1 == 3 || gs1 == 4
-            g1 = Sm(gs1, q1.niter, q1.eps, 2)
-        else
-            g1 = Sm(gs1, 2)
+            gs1 = read(data, Int32)
+            if gs1 != 0 && gs1 != 3 && gs1 != 4
+                gn1 = read(data, Int32)
+                geps1 = read(data, Float64)
+                g1 = Sm(gs1, gn1, geps1, 2)
+            elseif gs1 == 3 || gs1 == 4
+                g1 = Sm(gs1, q1.niter, q1.eps, 2)
+            else
+                g1 = Sm(gs1, 2)
+            end
+            
+            gs2 = read(data, Int32)
+            if gs2 != 0 && gs2 != 3 && gs2 != 4
+                gn2 = read(data, Int32)
+                geps2 = read(data, Float64)
+                g2 = Sm(gs2, gn2, geps2, 2)
+            elseif gs1 == 3 || gs1 == 4
+                g2 = Sm(gs2, q2.niter, q2.eps, 2)
+            else
+                g2 = Sm(gs2, 2)
+            end
+            
+            read!(data, aux_i)
+            a[k] = CHeader(aux_f, aux_i, theta, [q1, q2, g1, g2])
         end
-        
-        gs2 = read(data, Int32)
-        if gs2 != 0 && gs2 != 3 && gs2 != 4
-            gn2 = read(data, Int32)
-            geps2 = read(data, Float64)
-            g2 = Sm(gs2, gn2, geps2, 2)
-        elseif gs1 == 3 || gs1 == 4
-            g2 = Sm(gs2, q2.niter, q2.eps, 2)
-        else
-            g2 = Sm(gs2, 2)
+    else
+        aux_f = zeros(Float64, 4)
+        aux_i = zeros(Int32, 4)
+        for k = 1:gh.ncorr
+            read!(data, aux_f)
+            read!(data, aux_i)
+            a[k] = CHeader(aux_f, aux_i)
         end
-        
-        read!(data, aux_i)
-        a[k] = CHeader(aux_f, aux_i, theta, [q1, q2, g1, g2])
     end
     close(data)
     return a
 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::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing, legacy::Bool=false)
 
-    read_mesons(path::Vector{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, legacy::Bool=false)
 
 This function 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.
 ADerrors id can be specified as argument. If is not specified, the `id` is fixed according to the ensemble name (example: "H400"-> id = 400)
 
+*For the old version (without smearing, distance preconditioning and theta) set legacy=true
 Examples:
 ```@example
 read_mesons(path)
@@ -85,13 +97,12 @@ 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")
+read_mesons(path, "G5_d2", "G5_d2", legacy=true)
 ```
 """
-function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing)
-    
-    isnothing(g1) ? t1=nothing : t1 = findfirst(x-> x==g1, gamma_name) - 1
-    isnothing(g2) ? t2=nothing : t2 = findfirst(x-> x==g2, gamma_name) - 1
+function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing, legacy::Bool=false)    
+    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)
@@ -99,8 +110,8 @@ function read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union
     end
 
     data = open(path, "r")
-    g_header = read_global_header(path)
-    c_header = read_CHeader(path)
+    g_header = read_GHeader(path)
+    c_header = read_CHeader(path, legacy=legacy)
 
     ncorr = g_header.ncorr
     tvals = g_header.tvals
@@ -159,8 +170,8 @@ 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)
+function read_mesons(path::Vector{String}, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing, legacy::Bool=false)
+    res = read_mesons.(path, g1, g2, id=id, legacy=legacy)
     nrep = length(res)
     ncorr = length(res[1])
 

src/juobs_tools.jl

@@ -28,6 +28,54 @@ function apply_rw(data::Vector{<:Array{Float64}}, W::Vector{Matrix{Float64}})
     data_r = [data[k] .* rw1[k].* rw2[k] / rw_mean for k=1:length(data)]
     return data_r
 end
+
+function dobsdp(a::uwreal, p::uwreal) # Compute da / dp
+    
+    if count(p.prop .== true) != 1
+        error("I do not know how to compute this")
+    end
+    
+    for i = 1:min(length(a.der), length(p.der))
+        if (p.prop[i] && a.prop[i])
+            return a.der[i] / p.der[i]
+        end 
+    end 
+    
+    return 0.0
+end
+
+function check_corr_der(obs::Corr, derm::Vector{Corr})   
+    g1 = Vector{String}(undef, 0)
+    g2 = Vector{String}(undef, 0)
+
+    for d in derm
+        aux = [d.gamma[1], d.gamma[2]]
+        push!(g1, aux[1][1:end-3])
+        push!(g2, aux[2][1:end-3])
+    end
+
+    h = copy(derm)
+    push!(h, obs)
+    
+    
+    if any(getfield.(h, :y0) .!= getfield(h[1], :y0))
+        return false
+    end
+    for s in [:kappa, :mu]
+        for k = 1:2
+            if any(getindex.(getfield.(h, s), k) .!= getindex(getfield(h[1], s), k))
+                return false
+            end
+        end
+    end
+    #gamma check
+    if any(g1 .!= obs.gamma[1]) || any(g2 .!= obs.gamma[2])
+        return false
+    end
+
+    return true
+end
+
 @doc raw"""
     corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1)
 
@@ -99,23 +147,27 @@ end
 
 Computes the derivative of an observable A with respect to the sea quark masses.
 
-``d <A> / dm(sea) = \sum_i (d<A> / d<O_i>) * (d<O_i> / dm(sea))`` 
+``\frac{d <A>}{dm(sea)} = \sum_i \frac{\partial <A>}{\partial <O_i>}  \frac{d <O_i>}{d m(sea)}``
+
 
-`` d <O> / dm(sea) = <O> <dS/dm> - <O dS/dm> = - <(O - <O>) (dS/dm - <dS/dm>)>`` 
+``\frac{d <O_i>}{dm(sea)} = <O_i> <\frac{\partial S}{\partial m}> - <O_i \frac{\partial S}{\partial m}> 
+= - <(O_i - <O_i>) (\frac{\partial S}{\partial m} - <\frac{\partial S}{\partial m}>)>``
 
 where ``O_i`` are primary observables 
 
-md is a vector that contains the derivative of the action S with respect
-to the sea quark masses for each replica. md[irep][irw, icfg]
+`md` is a vector that contains the derivative of the action ``S`` with respect
+to the sea quark masses for each replica. `md[irep][irw, icfg]`
 
-md_sea returns a tuple of uwreal observables (dA/dml, dA/dms)|sea, 
-where ml and ms are the light and strange quark masses.
+`md_sea` returns a tuple of uwreal observables ``(dA/dm_l, dA/dm_s)|_sea``, 
+where ``m_l`` and ``m_s`` are the light and strange quark masses.
 
 ```@example
 #Single replica
 data = read_mesons(path, "G5", "G5")
 md = read_md(path_md)
-corr_pp = corr_obs.(data)
+rw = read_ms1(path_rw)
+
+corr_pp = corr_obs.(data, rw=rw)
 m = meff(corr_pp[1], plat)
 m_mdl, m_mds = md_sea(m, [md], ADerrors.wsg)
 m_shifted = m + 2 * dml * m_mdl + dms * m_mds
@@ -136,7 +188,7 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADer
     p = findall(t-> t==1, a.prop)
 
     if nid != 1
-        error("neid > 1")
+        error("Error: neid > 1")
     end
 
     id = ws.map_nob[p]
@@ -186,6 +238,59 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADer
     end
 end
 
+@doc raw"""
+    md_val(a::uwreal, obs::Corr, derm::Vector{Corr})
+
+Computes the derivative of an observable A with respect to the valence quark masses.
+
+``\frac{d <A>}{dm(val)} = \sum_i \frac{\partial <A>}{\partial <O_i>}  \frac{d <O_i>}{d m(val)}``
+
+``\frac{d <O_i>}{dm(val)} = <\frac{\partial O_i}{\partial m(val)}>``
+
+where ``O_i`` are primary observables 
+
+`md` is a vector that contains the derivative of the action ``S`` with respect
+to the sea quark masses for each replica. `md[irep][irw, icfg]`
+
+`md_val` returns a tuple of `uwreal` observables ``(dA/dm_1, dA/dm_2)|_val``, 
+where ``m_1`` and ``m_2`` are the correlator masses.
+
+```@example
+data = read_mesons(path, "G5", "G5", legacy=true)
+data_d1 = read_mesons(path, "G5_d1", "G5_d1", legacy=true)
+data_d2 = read_mesons(path, "G5_d2", "G5_d2", legacy=true)
+
+rw = read_ms1(path_rw)
+
+corr_pp = corr_obs.(data, rw=rw)
+corr_pp_d1 = corr_obs.(data_d1, rw=rw)
+corr_pp_d2 = corr_obs.(data_d2, rw=rw)
+derm = [[corr_pp_d1[k], corr_pp_d2[k]] for k = 1:length(pp_d1)]
+
+m = meff(corr_pp[1], plat)
+m_md1, m_md2 = md_val(m, corr_pp[1], derm[1])
+m_shifted = m + 2 * dm1 * m_md1 + dm2 * m_md2
+```
+"""
+function md_val(a::uwreal, obs::Corr, derm::Vector{Corr})
+    nid = neid(a)
+    if nid != 1
+        error("Error: neid > 1")
+    end
+    if length(derm) != 2
+        error("Error: length derm != 2")
+    end
+    if !check_corr_der(obs, derm)
+        error("Corr parameters does not match")
+    end
+
+    corr = getfield(obs, :obs)
+
+    der = [dobsdp(a, corr[k]) for k = 1:length(corr)]
+    derm1, derm2 = derm
+    return (sum(der .* derm1.obs), sum(der .* derm2.obs))
+end
+
 function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
     isnothing(wpm) ? uwerr.(obs) : uwerr.(obs, wpm)
     w = 1 ./ err.(obs)[plat[1]:plat[2]].^2

src/juobs_types.jl

@@ -12,7 +12,13 @@
 const noise_name=["Z2", "GAUSS", "U1", "Z4"]
 const gamma_name = ["G0", "G1", "G2", "G3",
 "invalid", "G5", "1", "G0G1", "G0G2", "G0G3",
-"G0G5", "G1G2", "G1G3", "G1G5", "G2G3", "G2G5", "G3G5"]
+"G0G5", "G1G2", "G1G3", "G1G5", "G2G3", "G2G5", "G3G5", 
+"G0_d1", "G1_d1", "G2_d1", "G3_d1",
+"invalid", "G5_d1", "1_d1", "G0G1_d1", "G0G2_d1", "G0G3_d1",
+"G0G5_d1", "G1G2_d1", "G1G3_d1", "G1G5_d1", "G2G3_d1", "G2G5_d1", "G3G5_d1",
+"G0_d2", "G1_d2", "G2_d2", "G3_d2",
+"invalid", "G5_d2", "1_d2", "G0G1_d2", "G0G2_d2", "G0G3_d2",
+"G0G5_d2", "G1G2_d2", "G1G3_d2", "G1G5_d2", "G2G3_d2", "G2G5_d2", "G3G5_d2"]
 const qs_name=["Local", "Wuppertal", "3D Gradient Flow", "Gradient Flow"]
 const gs_name=["Local", "APE", "3D Wilson Flow", "Quark 3D Gradient Flow", "Quark Gradient Flow"]
 
@@ -96,7 +102,30 @@ mutable struct CHeader
         
         return a
     end
+    function CHeader(aux_f::Vector{Float64}, aux_i::Vector{Int32})
+        a = new()
+        a.k1 = aux_f[1]
+        a.k2 = aux_f[2]
+        a.mu1 = aux_f[3]
+        a.mu2 = aux_f[4]
+        a.dp1 = 0.0
+        a.dp2 = 0.0
+        a.type1 = aux_i[1]
+        a.type2 = aux_i[2]
+        a.x0 = aux_i[3]
+        a.is_real = aux_i[4]
+
+        a.theta1 = zeros(3)
+        a.theta2 = zeros(3)
+        a.q1 = Sm(0, 1)
+        a.q2 = Sm(0, 1)
+        a.g1 = Sm(0, 2)
+        a.g2 = Sm(0, 2)
+        a.hsize = 8*4 + 4*4
+        a.dsize = 16 - 8* a.is_real 
 
+        return a
+    end
 end
 
 mutable struct CData
@@ -112,28 +141,31 @@ Base.copy(a::CData) = CData(a.header, a.vcfg, a.re_data, a.im_data, a.id)
 
 mutable struct Corr
     obs::Vector{uwreal}
+    kappa::Vector{Float64}
     mu::Vector{Float64}
     gamma::Vector{String}
     y0::Int64
     function Corr(a::Vector{uwreal}, b::CData)
         h = getfield(b, :header)
+        kappa = [h.k1, h.k2]
         mu = [h.mu1, h.mu2]
         gamma = [gamma_name[h.type1+1], gamma_name[h.type2+1]]
         y0 = Int64(h.x0)
-        return new(a, mu, gamma, y0)
+        return new(a, kappa, mu, gamma, y0)
     end 
     function Corr(a::Vector{uwreal}, b::Vector{CData})
-        sym = [:mu1, :mu2, :type1, :type2, :x0]
+        sym = [:k1, :k2, :mu1, :mu2, :type1, :type2, :x0]
         h = getfield.(b, :header)
         for s in sym
             if !all(getfield.(h, s) .== getfield(h[1], s))
                 error("Corr: Parameter mismatch")
             end
         end
+        kappa = [h[1].k1, h[1].k2]
         mu = [h[1].mu1, h[1].mu2]
         gamma = [gamma_name[h[1].type1+1], gamma_name[h[1].type2+1]]
         y0 = Int64(h[1].x0)
-        return new(a, mu, gamma, y0)
+        return new(a, kappa, mu, gamma, y0)
     end 
 end