md_sea and md_val rework

src/juobs_obs.jl

@@ -461,7 +461,7 @@ 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, ws::ADerrors.wspace=ADerrors.wsg,
-    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, rw_info::Bool=false)
 
     Ysl = Y.obs
     t = Y.t
@@ -496,6 +496,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
     if !isnothing(rw)
         Ysl_r, W = apply_rw(Ysl, rw)
         W_obs = uwreal(W, id)
+        WY_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
     end
 
     for i = 1:xmax
@@ -504,7 +505,8 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
             if isnothing(rw)
                 Y_aux[i, k] = uwreal(Ysl[:, i, j], id)
             else
-                Y_aux[i, k] = uwreal(Ysl_r[:, i, j], id) / W_obs
+                WY_aux[i, k] = uwreal(Ysl_r[:, i, j], id)
+                Y_aux[i, k] = WY_aux[i, k] / W_obs
             end
             k = k + 1
         end
@@ -546,7 +548,11 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
         title(string(L"$t/a^2 = $", t[nt0]))
         display(gcf())
     end
+    if rw_info && !isnothing(rw)
+        (t0, WY_aux, W_obs)
+    else
         return t0
+    end
 end
 
 function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false, 
@@ -601,6 +607,7 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
         [tmp_r = cat(tmp_r, Ysl_r[k], dims=1) for k = 2:nr]
         [tmp_W = cat(tmp_W, W[k], dims=1) for k = 2:nr]
         W_obs = uwreal(tmp_W, id, replica)
+        WY_aux = Matrix{uwreal}(undef, xmax, 2*dt0+1)
     end
     for i = 1:xmax
         k = 1
@@ -608,7 +615,8 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
             if isnothing(rw)
                 Y_aux[i, k] = uwreal(tmp[:, i, j], id, replica)
             else
-                Y_aux[i, k] = uwreal(tmp_r[:, i, j], id, replica) / W_obs
+                WY_aux[i, k] = uwreal(tmp_r[:, i, j], id, replica)
+                Y_aux[i, k] = WY_aux[i, k] / W_obs
             end
             k = k + 1
         end
@@ -650,7 +658,11 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
         title(string(L"$t/a^2 = $", t[nt0]))
         display(gcf())
     end
+    if rw_info && !isnothing(rw)
+        (t0, WY_aux, W_obs)
+    else
         return t0
+    end
 end
 
 function t0_guess(t::Vector{Float64}, Ysl::Array{Float64, 3}, plat::Vector{Int64}, L::Int64)

src/juobs_tools.jl

@@ -81,25 +81,30 @@ corr_pp = corr_obs.(data)
 corr_pp_r = corr_obs.(data, rw=[rw1, rw2])
 ```
 """
-function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1)
+function corr_obs(cdata::CData; real::Bool=true, rw::Union{Array{Float64, 2}, Nothing}=nothing, L::Int64=1, rw_info::Bool=false)
 
     real ? data = cdata.re_data ./ L^3 : data = cdata.im_data ./ L^3
     nt = size(data)[2]
-    obs = Vector{uwreal}(undef, nt)
 
     if isnothing(rw)
-        [obs[x0] = uwreal(data[:, x0], cdata.id) for x0 = 1:nt]
+        obs = [uwreal(data[:, x0], cdata.id) for x0 = 1:nt]
     else
         data_r, W = apply_rw(data, rw)
+        ow = [uwreal(data_r[:, x0], cdata.id) for x0 = 1:nt]
         W_obs = uwreal(W, cdata.id)
-        [obs[x0] = uwreal(data_r[:, x0], cdata.id) / W_obs for x0 = 1:nt]
+        obs = [ow[x0] / W_obs for x0 = 1:nt]
     end
+
+    if rw_info && !isnothing(rw)
+        return (Corr(obs, cdata), ow, W_obs)
+    else
         return Corr(obs, cdata)
+    end
 end
 
 #function corr_obs for R != 1
 #TODO: vcfg with gaps
-function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1)
+function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1, rw_info::Bool=false)
     nr = length(cdata)
     id = getfield.(cdata, :id)
     vcfg = getfield.(cdata, :vcfg)
@@ -112,12 +117,11 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
     real ? data = getfield.(cdata, :re_data) ./ L^3 : data = getfield.(cdata, :im_data) ./ L^3
 
     nt = size(data[1])[2]
-    obs = Vector{uwreal}(undef, nt)
 
     if isnothing(rw)
         tmp = data[1]
         [tmp = cat(tmp, data[k], dims=1) for k = 2:nr]
-        [obs[x0] = uwreal(tmp[:, x0], id[1], replica) for x0 = 1:nt]
+        obs = [uwreal(tmp[:, x0], id[1], replica) for x0 = 1:nt]
     else
         data_r, W = apply_rw(data, rw)
         tmp = data_r[1]
@@ -125,11 +129,15 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
         [tmp = cat(tmp, data_r[k], dims=1) for k = 2:nr]
         [tmp_W = cat(tmp_W, W[k], dims=1) for k = 2:nr]
 
+        ow = [uwreal(tmp[:, x0], id[1], replica) for x0 = 1:nt]
         W_obs = uwreal(tmp_W, id[1], replica)
-        [obs[x0] = uwreal(tmp[:, x0], id[1], replica) / W_obs for x0 = 1:nt]
+        obs = [ow[x0] / W_obs for x0 = 1:nt]
     end
-
+    if rw_info && !isnothing(rw)
+        return (Corr(obs, cdata), ow, W_obs)
+    else
         return Corr(obs, cdata)
+    end
 end    
 @doc raw"""
     corr_sym(corrL::Corr, corrR::Corr, parity::Int64=1)
@@ -202,7 +210,7 @@ m_mdl, m_mds = md_sea(m, [md1, md2], ADerrors.wsg)
 m_shifted = m + 2 * dml * m_mdl + dms * m_mds
 ```
 """
-function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADerrors.wsg)
+function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::uwreal, w::Union{uwreal, Nothing}=nothing, ws::ADerrors.wspace=ADerrors.wsg)
     nid = neid(a)
     p = findall(t-> t==1, a.prop)
 
@@ -233,22 +241,40 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADer
         md_aux = cat(md_aux, md[k][:, 1:ivrep[k]], dims=2)
     end
     
-    fluc_md = md_aux .- mean(md_aux, dims=2)
-    uwerr(a)
-    fluc_obs = mchist(a, id)
+    nrw = size(md_aux, 1)
+    uwerr(ow)
+    ow_data = ow.mean .+ mchist(ow, id)
+
+    if isnothing(w)
+        der = derivative(a, ow)
+        d1 = der * (ow * uwreal(md_aux[1, :], id, ivrep) - uwreal(ow_data .* md_aux[1, :], id, ivrep))
         
-    nrw = size(fluc_md, 1)
         if nrw == 1
-        der1 = uwreal(-fluc_md[1, :] .* fluc_obs, id, ivrep)
-        return (der1, der1)
+            return (d1, d1)
         elseif nrw == 2
-        der1 = uwreal(-fluc_md[1, :] .* fluc_obs, id, ivrep)
-        der2 = uwreal(-fluc_md[2, :] .* fluc_obs, id, ivrep)
-        return (der1, der2)
+            d2 = der * (ow * uwreal(md_aux[2, :], id, ivrep) - uwreal(ow_data .* md_aux[2, :], id, ivrep))
+            return (d1, d2)
+        end
     else
-        return nothing
+        uwerr(w)
+        der = derivative(a, ow) * w.mean
+        w_data = w.mean .+ mchist(w, id)
+        d1 = der * (ow * uwreal(w_data .* md_aux[1, :], id, ivrep) / w^2 - uwreal(ow_data .* md_aux[1, :], id, ivrep) / w)
+        
+        if nrw == 1
+            return (d1, d1)
+        elseif nrw == 2
+            d2 = der * (ow * uwreal(w_data .* md_aux[2, :], id, ivrep) / w^2 - uwreal(ow_data .* md_aux[2, :], id, ivrep) / w)
+            return (d1, d2)
+        end
     end
 
+    return nothing
+end
+
+function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ow::Array{uwreal}, w::Union{uwreal, Nothing}=nothing, ws::ADerrors.wspace=ADerrors.wsg)
+    d = [md_sea(a, md, ow_, w, ws) for ow_ in ow]
+    return (sum(getindex.(d, 1)), sum(getindex.(d, 2)))
 end
 
 @doc raw"""
@@ -298,8 +324,34 @@ function md_val(a::uwreal, obs::Corr, derm::Vector{Corr})
     end
 
     corr = getfield(obs, :obs)
+    prop = getfield.(corr, :prop)
     
+    if all(count.(prop) .== 1)
         der = [derivative(a, corr[k]) for k = 1:length(corr)]
+    elseif all(count.(prop) .== 2)
+        corr_der = getfield.(corr, :der)
+        n = findall.(t-> t==1, prop)
+        n = vcat(n'...)
+        if all(n[:, 1] .== n[1, 1]) # find ow and w
+            n_w = n[1, 1]
+            n_ow = n[:, 2]
+        else
+            n_w = n[1, 2]
+            n_ow = n[:, 1]
+        end
+        w_mean = 1 / getindex.(corr_der, n_ow)[1]
+        #= Extra information (ow, w)
+        ws = ADerrors.wsg
+        ow_mean = - getindex.(corr_der, n_w) .* w_mean^2
+        w_data = w_mean .+ ws.fluc[n_w].delta
+        fluc_ow = getfield.(ws.fluc[n_ow], :delta)
+        fluc_ow = vcat(fluc_ow'...)
+        ow_data = ow_mean .+ ws.fluc[n_ow]
+        =#
+        der = a.der[n_ow] * w_mean
+    else
+        return nothing
+    end
     derm1, derm2 = derm
     return (sum(der .* derm1.obs), sum(der .* derm2.obs))
 end