md_sea added + fix md_reader

md_reader sign problem fixed

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, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
+export corr_obs, md_sea, 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

@@ -225,6 +225,8 @@ read_md(path::String)
 
 Reads openQCD  pbp.dat files at a given path. This method returns a matrix md[irw, icfg] that contains the derivatives dS/dm, where
 md[irw=1] = dS/dm_l and md[irw=2] = dS/dm_s
+Seff = -tr(log(D+m))
+dSeff/ dm = -tr((D+m)^-1)
 
 Examples:
 ```@example    
@@ -236,7 +238,7 @@ function read_md(path::String)
     nrw = length(r_data)
     ncnfg = size(r_data[1])[3]
     md = zeros(Float64, nrw, ncnfg)
-    [md[k, :] = prod(mean(r_data[k], dims=2), dims=1) for k = 1:nrw]
+    [md[k, :] = prod(.-mean(r_data[k], dims=2), dims=1) for k = 1:nrw]
     return md
 end
 

src/juobs_tools.jl

@@ -88,6 +88,104 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
     return Corr(obs, cdata)
 end    
 
+@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]
+
+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```
+#Single replica
+data = read_mesons(path, "G5", "G5")
+md = read_md(path_md)
+corr_pp = corr_obs.(data)
+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")
+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)
+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
+    end
+
+    id = ws.map_nob[p]
+    if !all(id .== id[1])
+        println("ids do not match")
+        return nothing
+    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
+    end
+    ivrep = ivrep[1]
+
+    if length(md) != length(ivrep)
+        println("Nr obs != Nr md")
+        return nothing
+    end
+
+    md_aux = md[1][:, 1:ivrep[1]]
+    for k = 2:length(md)
+        md_aux = cat(md_aux, md[k][:, 1:ivrep[k]], dims=2)
+    end
+
+    fluc_obs = getfield.(ws.fluc[p], :delta)
+    fluc_md = md_aux .- mean(md_aux, dims=2)
+    
+    nrw = size(fluc_md, 1)
+    d = a.der[p]
+    nobs = sum(a.prop)
+    if nrw == 1
+        der1 = uwreal(0)
+        for k = 1:nobs
+            der1 = der1 - d[k] * uwreal(fluc_md[1, :] .* fluc_obs[k], id, ivrep)
+        end
+        return (der1, der1)
+    elseif nrw == 2
+        der1 = uwreal(0)
+        der2 = uwreal(0)
+        for k = 1:nobs
+            der1 = der1 - d[k] * uwreal(fluc_md[1, :] .* fluc_obs[k], id, ivrep)
+            der2 = der2 - d[k] * uwreal(fluc_md[2, :] .* fluc_obs[k], id, ivrep)
+        end
+        return (der1, der2)
+    else
+        return nothing
+    end
+end
+
 function plat_av(obs::Vector{uwreal}, plat::Vector{Int64})
     uwerr.(obs)
     w = 1 ./ err.(obs)[plat[1]:plat[2]].^2