#TODO: apply_rw with gaps function apply_rw(data::Array{Float64}, W::Matrix{Float64}) nc = size(data, 1) W1 = W[1, 1:nc] W2 = W[2, 1:nc] data_r = data .* W1 .* W2 / mean(W1 .* W2) return data_r end function apply_rw(data::Vector{<:Array{Float64}}, W::Vector{Matrix{Float64}}) if length(W) != length(data) error("Lenghts must match") end nc = size.(data, 1) rw1 = [W[k][1, 1:nc[k]] for k=1:length(W)] rw2 = [W[k][2, 1:nc[k]] for k=1:length(W)] rw1_cat = rw1[1] rw2_cat = rw2[1] for k = 2:length(W) rw1_cat = cat(rw1_cat, rw1[k], dims=1) rw2_cat = cat(rw2_cat, rw2[k], dims=1) end rw_mean = mean(rw1_cat .* rw2_cat) data_r = [data[k] .* rw1[k].* rw2[k] / rw_mean for k=1:length(data)] return data_r 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) corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array{Float64, 2}, 1}, Nothing}=nothing, L::Int64=1) 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. 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. ```@example #Single replica data = read_mesons(path, "G5", "G5") rw = read_ms1(path_rw) corr_pp = corr_obs.(data) corr_pp_r = corr_obs.(data, rw=rw) #Two replicas data = read_mesons([path_r1, path_r2], "G5", "G5") rw1 = read_ms1(path_rw1) rw2 = read_ms1(path_rw2) 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) real ? data = cdata.re_data ./ L^3 : data = cdata.im_data ./ L^3 data_r = isnothing(rw) ? data : apply_rw(data, rw) nt = size(data)[2] obs = Vector{uwreal}(undef, nt) [obs[x0] = uwreal(data_r[:, x0], cdata.id) for x0 = 1:nt] return Corr(obs, cdata) 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) nr = length(cdata) id = getfield.(cdata, :id) vcfg = getfield.(cdata, :vcfg) replica = Int64.(maximum.(vcfg)) if !all(id .== id[1]) error("IDs are not equal") end real ? data = getfield.(cdata, :re_data) ./ L^3 : data = getfield.(cdata, :im_data) ./ L^3 data_r = isnothing(rw) ? data : apply_rw(data, rw) tmp = data_r[1] [tmp = cat(tmp, data_r[k], dims=1) for k = 2:nr] nt = size(data[1])[2] obs = Vector{uwreal}(undef, nt) [obs[x0] = uwreal(tmp[:, x0], id[1], replica) for x0 = 1:nt] return Corr(obs, cdata) end @doc raw""" corr_sym(corrL::Corr, corrR::Corr, parity::Int64=1) Computes the symmetrized correlator using the left correlador `corrL` and the right correlator `corrR`. The source position of `corrR` must be `T - 1 - y0`, where `y0` is the source position of `corrL`. ```@example pp_sym = corr_sym(ppL, ppR, +1) a0p_sym = corr_sym(a0pL, a0pR, -1) ``` """ function corr_sym(corrL::Corr, corrR::Corr, parity::Int64=1) T = length(corrL.obs) sym = [:kappa, :mu, :gamma] if corrL.y0 != T - 1 - corrR.y0 error("Corr: Parameter mismatch") end for s in sym if getfield(corrL, s) != getfield(corrR, s) error("Corr: Parameter mismatch") end end if abs(parity) != 1 error("incorrect value of parity (+- 1)") end res = (corrL.obs[1:end] + parity * corrR.obs[end:-1:1]) / 2 return Corr(res, corrL.kappa, corrL.mu, corrL.gamma, corrL.y0) end #TODO: VECTORIZE, uwreal? @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. ``\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}>)>`` 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) 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 = read_mesons([path_r1, path_r2], "G5", "G5") md1 = read_md(path_md1) md2 = read_md(path_md2) 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 ``` """ 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 error("Error: neid > 1") end id = ws.map_nob[p] if !all(id .== id[1]) error("ids do not match") end id = ws.id2str[id[1]] ivrep = getfield.(ws.fluc[p], :ivrep) ivrep1 = fill(ivrep[1], length(ivrep)) if !all(ivrep .== ivrep1) error("ivreps do not match") end ivrep = ivrep[1] if length(md) != length(ivrep) error("Nr obs != Nr md") end #md_aux as a Matrix + Automatic truncation 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) uwerr(a) fluc_obs = mchist(a, id) nrw = size(fluc_md, 1) if nrw == 1 der1 = uwreal(-fluc_md[1, :] .* fluc_obs, id, ivrep) return (der1, der1) elseif nrw == 2 der1 = uwreal(-fluc_md[1, :] .* fluc_obs, id, ivrep) der2 = uwreal(-fluc_md[2, :] .* fluc_obs, id, ivrep) return (der1, der2) else return nothing 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 = [derivative(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_aux, wpm) for obs_aux in obs] w = 1 ./ err.(obs)[plat[1]:plat[2]].^2 av = sum(w .* obs[plat[1]:plat[2]]) / sum(w) return av end function lin_fit(x::Vector{<:Real}, v::Vector{Float64}, e::Vector{Float64}) sig2 = e .* e S = sum(1 ./ sig2) Sx = sum(x ./ sig2) Sy = sum(v ./ sig2) Sxy = sum(v .* x ./ sig2) Sxx = sum(x .* x ./sig2) delta = S * Sxx - Sx*Sx par = [Sxx*Sy-Sx*Sxy, S*Sxy-Sx*Sy] ./delta #C = [[Sxx/delta, -Sx/delta], [-Sx/delta, S/delta]] return par end @doc raw""" lin_fit(x::Vector{<:Real}, y::Vector{uwreal}) Computes a linear fit of uwreal data points y. This method return uwreal fit parameters and chisqexpected. ```@example fitp, csqexp = lin_fit(phi2, m2) m2_phys = fitp[1] + fitp[2] * phi2_phys ``` """ function lin_fit(x::Vector{<:Real}, y::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) isnothing(wpm) ? uwerr.(y) : [uwerr(yaux, wpm) for yaux in y] par = lin_fit(x, value.(y), err.(y)) chisq(p, d) = sum((d .- p[1] .- p[2].*x).^2 ./ err.(y) .^2) (fitp, csqexp) = fit_error(chisq, par, y) for i = 1:length(fitp) isnothing(wpm) ? uwerr(fitp[i]) : uwerr(fitp[i], wpm) print("\n Fit parameter: ", i, ": ") details(fitp[i]) end println("Chisq / chiexp: ", chisq(par, y), " / ", csqexp, " (dof: ", length(x)-length(par),")") return (fitp, csqexp) end @doc raw""" x_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64}) Computes the results of a linear interpolation/extrapolation in the x axis """ x_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64}) = (y - par[1]) / par[2] @doc raw""" y_lin_fit(par::Vector{uwreal}, y::Union{uwreal, Float64}) Computes the results of a linear interpolation/extrapolation in the y axis """ y_lin_fit(par::Vector{uwreal}, x::Union{uwreal, Float64}) = par[1] + par[2] * x @doc raw""" fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, covar::Bool=false) Given a model function with a number param of parameters and an array of `uwreal`, this function fit ydata with the given `model` and print fit information The method return an array `upar` with the best fit parameters with their errors. The flag `wpm` is an optional array of Float64 of lenght 4. The first three paramenters specify the criteria to determine the summation windows: - `vp[1]`: The autocorrelation function is summed up to ``t = round(vp[1])``. - `vp[2]`: The sumation window is determined using U. Wolff poposal with ``S_\tau = wpm[2]`` - `vp[3]`: The autocorrelation function ``\Gamma(t)`` is summed up a point where its error ``\delta\Gamma(t)`` is a factor `vp[3]` times larger than the signal. An additional fourth parameter `vp[4]`, tells ADerrors to add a tail to the error with ``\tau_{exp} = wpm[4]``. Negative values of `wpm[1:4]` are ignored and only one component of `wpm[1:3]` needs to be positive. If the flag `covar`is set to true, `fit_routine` takes into account covariances between x and y for each data point. ```@example @. model(x,p) = p[1] + p[2] * exp(-(p[3]-p[1])*x) @. model2(x,p) = p[1] + p[2] * x[:, 1] + (p[3] + p[4] * x[:, 1]) * x[:, 2] fit_routine(model, xdata, ydata, param=3) fit_routine(model, xdata, ydata, param=3, covar=true) ``` """ function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) isnothing(wpm) ? uwerr.(ydata) : [uwerr(yaux, wpm) for yaux in ydata] yval = value.(ydata) yer = err.(ydata) # Generate chi2 + solver chisq = gen_chisq(model, xdata, yer) fit = curve_fit(model, xdata, yval, 1.0 ./ yer.^2, fill(0.5, param)) (upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata) : fit_error(chisq, coef(fit), ydata, wpm) #Info for i = 1:length(upar) isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm) print("\n Fit parameter: ", i, ": ") details(upar[i]) end println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", length(yval) - param,")") return upar end function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, covar::Bool=false) Nalpha = size(xdata, 2) # number of x-variables Ndata = size(xdata, 1) # number of datapoints if isnothing(wpm) uwerr.(ydata) uwerr.(xdata) else [uwerr(yaux, wpm) for yaux in ydata] [uwerr(xaux, wpm) for xaux in xdata] end yval = value.(ydata) yer = err.(ydata) xval = value.(xdata) xer = err.(xdata) dat = Vector{Float64}(undef, Ndata * (Nalpha+1)) ddat = Vector{Float64}(undef, Ndata * (Nalpha+1)) data = Vector{uwreal}(undef, Ndata * (Nalpha+1)) for i = 1:Nalpha dat[(i-1)*Ndata+1:i*Ndata] = xval[:, i] ddat[(i-1)*Ndata+1:i*Ndata] = xer[:, i] data[(i-1)*Ndata+1:i*Ndata] = xdata[:, i] end dat[Nalpha*Ndata+1:end] = yval ddat[Nalpha*Ndata+1:end] = yer data[Nalpha*Ndata+1:end] = ydata # Guess fit = curve_fit(model, xval, yval, 1.0 ./ yer.^2, fill(0.5, param)) # Generate chi2 + solver if covar aux = Vector{Vector{uwreal}}(undef, Ndata) for k = 1:Ndata aux[k] = Vector{uwreal}(undef, Nalpha+1) for i = 1:Nalpha aux[k][i] = xdata[k, i] end aux[k][Nalpha+1] = ydata[k] end C = isnothing(wpm) ? [ADerrors.cov(aux[k]) for k = 1:Ndata] : [ADerrors.cov(aux[k], wpm) for k = 1:Ndata] chisq_full_cov(p, d) = get_chi2_cov(model, d, C, p, Nalpha) min_fun_cov(t) = chisq_full_cov(t, dat) sol = optimize(min_fun_cov, vcat(fit.param, dat[1:Nalpha*Ndata]), LevenbergMarquardt()) (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full_cov, sol.minimizer, data) : fit_error(chisq_full_cov, sol.minimizer, data, wpm) println("Chisq / chiexp: ", min_fun_cov(sol.minimizer), " / ", chi2_exp, " (dof: ", length(ydata) - param,")") else chisq_full(p, d) = get_chi2(model, d, ddat, p, Nalpha) min_fun(t) = chisq_full(t, dat) sol = optimize(min_fun, vcat(fit.param, dat[1:Nalpha*Ndata]), LevenbergMarquardt()) (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full, sol.minimizer, data) : fit_error(chisq_full, sol.minimizer, data, wpm) println("Chisq / chiexp: ", min_fun(sol.minimizer), " / ", chi2_exp, " (dof: ", length(ydata) - param,")") end #### chisq_full, min_fun out of conditional -> #### COMPILER WARNING ** incremental compilation may be fatally broken for this module ** # Info for i = 1:length(upar) isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm) print("\n Fit parameter: ", i, ": ") details(upar[i]) end return upar end function gen_chisq(f::Function, x::Array{<:Real}, err::Vector{Float64}) #constrained chisq(par, dat) = sum((dat .- f(x,par)).^2 ./err.^2) return chisq end function get_chi2(f::Function, data, ddata, par, Nalpha) #full chi2 = 0.0 Ndata = div(length(data), Nalpha+1) Npar = length(par) - Ndata * Nalpha p = par[1:Npar] for k = 1:Ndata xx = [par[Npar + k + (i-1)*Ndata] for i = 1:Nalpha] Cinv = zeros(Nalpha+1, Nalpha+1) [Cinv[i, i] = 1 / ddata[k + (i-1)*Ndata]^2 for i = 1:Nalpha+1] xx = [par[Npar + k + (i-1)*Ndata] for i = 1:Nalpha] delta = [data[k + (i-1)*Ndata] - xx[i] for i = 1:Nalpha] yy = f(xx', p) push!(delta, data[k + Nalpha*Ndata] - yy[1]) chi2 += delta' * Cinv * delta end return chi2 end function get_chi2_cov(f::Function, data, C, par, Nalpha) # full + cov chi2 = 0.0 Ndata = div(length(data), Nalpha+1) Npar = length(par) - Ndata * Nalpha p = par[1:Npar] for k = 1:Ndata if det(C[k]) / prod(diag(C[k])) > 1e-6 Cinv = inv(C[k]) else Cinv = zeros(Nalpha+1, Nalpha+1) [Cinv[i, i] = 1 / C[k][i, i] for i = 1:Nalpha+1] end xx = [par[Npar + k + (i-1)*Ndata] for i = 1:Nalpha] delta = [data[k + (i-1)*Ndata] - xx[i] for i = 1:Nalpha] yy = f(xx', p) push!(delta, data[k + Nalpha*Ndata] - yy[1]) chi2 += delta' * Cinv * delta end return chi2 end