md_val added

src/juobs.jl

@@ -9,7 +9,7 @@ include("juobs_obs.jl")
 
 export read_mesons, read_ms1, read_ms, read_md
 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 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_tools.jl

@@ -18,6 +18,54 @@ function apply_rw(data::Array{Float64}, W::Array{Float64, 2})
     data_r = data .* W1 .* W2 / mean(W1 .* W2)
     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)
 
@@ -88,72 +136,71 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
     return Corr(obs, cdata)
 end    
 
+#TODO: VECTORIZE
 @doc raw"""
 md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADerrors.wsg)
 
 Computes the derivative of an observable A with respect to the sea quark masses.
-d <A> / dm(sea) = sum_i (d<A> / d<Oi>) * (d<Oi> / dm(sea)) 
-d <O> / dm(sea) = <O> <dS/dm> - <O dS/dm> = - <(O - <O>) (dS/dm - <dS/dm>)>  
-where <Oi> 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]
+``\frac{d <A>}{dm(sea)} = \sum_i \frac{\partial <A>}{\partial <O_i>}  \frac{d <O_i>}{d m(sea)}``
 
-md_sea returns a tuple of uwreal observables (dA/dml, dA/dms)|sea, 
-where ml and ms are the light and strange quark masses.
 
-@example```
+``\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_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
 
 #Two replicas
-data_r1 = read_mesons(path_r1, "G5", "G5")
-data_r2 = read_mesons(path_r2, "G5", "G5")
+data = read_mesons([path_r1, path_r2], "G5", "G5")
 md1 = read_md(path_md1)
 md2 = read_md(path_md2)
 
-
-cdata = [[data_r1[k], data_r2[k]] for k=1:length(data_r1)]
-
-corr_pp = corr_obs.(cdata)
+corr_pp = corr_obs.(data)
 m = meff(corr_pp[1], plat)
 m_mdl, m_mds = md_sea(m, [md1, md2], ADerrors.wsg)
 m_shifted = m + 2 * dml * m_mdl + dms * m_mds
 ```
 """
-#TODO: VECTORIZE
 function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADerrors.wsg)
     nid = neid(a)
     p = findall(t-> t==1, a.prop)
 
     if nid != 1
-        println("Error: neid > 1")
-        return nothing
+        error("Error: neid > 1")
     end
 
     id = ws.map_nob[p]
     if !all(id .== id[1])
-        println("ids do not match")
-        return nothing
+        error("ids do not match")
     end
     id = id[1]
 
     ivrep = getfield.(ws.fluc[p], :ivrep)
     ivrep1 = fill(ivrep[1], length(ivrep))
     if !all(ivrep .== ivrep1)
-        println("ivreps do not match")
-        return nothing
+        error("ivreps do not match")
     end
     ivrep = ivrep[1]
 
     if length(md) != length(ivrep)
-        println("Nr obs != Nr md")
-        return nothing
+        error("Nr obs != Nr md")
     end
 
     md_aux = md[1][:, 1:ivrep[1]]
@@ -186,6 +233,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})
     uwerr.(obs)
     w = 1 ./ err.(obs)[plat[1]:plat[2]].^2

src/juobs_types.jl

@@ -141,18 +141,20 @@ 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))
@@ -160,10 +162,11 @@ mutable struct Corr
                 return nothing
             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