doc updated (rwf rework)

src/juobs_obs.jl

@@ -430,23 +430,26 @@ function get_model(x, p, n)
     return s
 end
 @doc raw"""
-    comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Matrix{Float64}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)
 
-    comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Union{Vector{Matrix{Float64}}, Nothing}=nothing, npol::Int64=2, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, info::Bool=false)
 
 Computes `t0` using the energy density of the action `Ysl`(Yang-Mills action).
 `t0` is computed in the plateau `plat`.
 A polynomial interpolation in `t` is performed to find `t0`, where `npol` is the degree of the polynomial (linear fit by default)
 
-The flag `pl` allows to show the plot.
+The flag `pl` allows to show the plot. 
+
+The flag `info` provides extra output that contains information about the primary observables. The function returns the primary observables ``<WY>`` and ``<W>``
+(it returns the observable <Y> if rw=nothing)
 
 ```@example
 #Single replica
 Y = read_ms(path)
 rw = read_ms(path_rw)
 
-t0 = comp_t0(Y, [38, 58], L=32)
-t0_r = comp_t0(Y, [38, 58], L=32, rw=rw)
+t0, Yobs = comp_t0(Y, [38, 58], L=32, info=true)
+t0_r, WYobs, Wobs = comp_t0(Y, [38, 58], L=32, rw=rw, info=true)
 
 #Two replicas
 Y1 = read_ms(path1)

src/juobs_tools.jl

@@ -54,9 +54,9 @@ function check_corr_der(obs::Corr, derm::Vector{Corr})
 end
 
 @doc raw"""
-    corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1)
+    corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1, info::Bool=false)
 
-    corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1)
+    corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, info::Bool=false)
 
 Creates a `Corr` struct with the given `CData` struct `cdata` (`read_mesons`) for a single replica.
 An array of `CData` can be passed as argument for multiple replicas.
@@ -65,6 +65,9 @@ The flag `real` select the real or imaginary part of the correlator.
 If `rw` is specified, the method applies reweighting. `rw` is passed as a matrix of Float64 (`read_ms1`)
 The correlator can be normalized with the volume factor if `L` is fixed.
 
+The flag `info` provides extra output that contains information about the primary observables. The function returns the primary observables ``<WO>`` and ``<W>``
+(it returns the observable <O> if rw=nothing)
+
 ```@example
 #Single replica
 data = read_mesons(path, "G5", "G5")
@@ -72,6 +75,12 @@ rw = read_ms1(path_rw)
 corr_pp = corr_obs.(data)
 corr_pp_r = corr_obs.(data, rw=rw)
 
+#Single replica + Info
+data = read_mesons(path, "G5", "G5")
+rw = read_ms1(path_rw)
+corr_pp, O = corr_obs(data[1], info=true)
+corr_pp_r, WO, W = corr_obs(data[1], rw=rw, info=true)
+
 #Two replicas
 data = read_mesons([path_r1, path_r2], "G5", "G5")
 rw1 = read_ms1(path_rw1)
@@ -174,17 +183,22 @@ function corr_sym(corrL::Corr, corrR::Corr, parity::Int64=1)
 end
 #TODO: VECTORIZE, uwreal?
 @doc raw"""
-    md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADerrors.wsg)
+    md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::uwreal, w::Union{uwreal, Nothing}=nothing, ws::ADerrors.wspace=ADerrors.wsg)
+
+    md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::Array{uwreal}, w::Union{uwreal, Nothing}=nothing, ws::ADerrors.wspace=ADerrors.wsg)
 
 Computes the derivative of an observable A with respect to the sea quark masses.
 
 ``\frac{d <A>}{dm(sea)} = \sum_i \frac{\partial <A>}{\partial <O_i>}  \frac{d <O_i>}{d m(sea)}``
 
 
-``\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}>)>``
+``\frac{d <O_i>}{dm(sea)} = <O_i> <\frac{\partial S}{\partial m}> - <O_i \frac{\partial S}{\partial m}>``
 
-where ``O_i`` are primary observables 
+where ``O_i`` are primary observables. The function returns  ``\frac{\partial <A>}{\partial <O_i>}  \frac{d <O_i>}{d m(sea)}``, where
+the primary observable ``O_i`` is specified as an input argument.
+
+`ow` is a primary observable or a vector of primary observables that contains ``<OW>`` (``<O>`` if reweighting is not applied) and 
+`w` is a primary obserable that contains the reweighting factors ``<W>`` (`w`=nothing if reweighting is not applied).
 
 `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]`
@@ -198,9 +212,9 @@ data = read_mesons(path, "G5", "G5")
 md = read_md(path_md)
 rw = read_ms1(path_rw)
 
-corr_pp = corr_obs.(data, rw=rw)
+corr_pp, wpp, w = corr_obs.(data[1], rw=rw, info=true)
 m = meff(corr_pp[1], plat)
-m_mdl, m_mds = md_sea(m, [md], ADerrors.wsg)
+m_mdl, m_mds = md_sea(m, [md], wpp, w, ADerrors.wsg)
 m_shifted = m + 2 * dml * m_mdl + dms * m_mds
 
 #Two replicas
@@ -208,9 +222,9 @@ data = read_mesons([path_r1, path_r2], "G5", "G5")
 md1 = read_md(path_md1)
 md2 = read_md(path_md2)
 
-corr_pp = corr_obs.(data)
+corr_pp, pp = corr_obs(data[1], info=true)
 m = meff(corr_pp[1], plat)
-m_mdl, m_mds = md_sea(m, [md1, md2], ADerrors.wsg)
+m_mdl, m_mds = md_sea(m, [md1, md2], pp, ADerrors.wsg)
 m_shifted = m + 2 * dml * m_mdl + dms * m_mds
 ```
 """