Skip to content

J
juobs
Commit 8aa2abdb authored by Antonino D'Anna's avatar Antonino D'Anna
Browse files
Options

changes in fit_routines plus tests

parent bf1b33cf
Show whitespace changes
Inline Side-by-side
Showing with 230 additions and 738 deletions
dev/read_meson_vec_gamma.jl deleted 100644 → 0
View file @ bf1b33cf
using juobs
function dev_read_mesons(path::String, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,
nnoise_trunc::Union{Int64, Nothing}=nothing)
type = Vector{Tuple{Int64,Int64}}(undef,length(gamma))
for i in eachindex(gamma)
G1,G2 = gamma[i]
type[i] = (findfirst(x-> x==G1, juobs.gamma_name) - 1, findfirst(x-> x==G2, juobs.gamma_name) - 1)
end
if isnothing(id)
bname = basename(path)
m = findfirst(r"[A-Z][0-9]{3}r[0-9]{3}", bname)
id = bname[m[1:4]]
end
data = open(path, "r")
g_header = juobs.read_GHeader(path)
c_header = juobs.read_CHeader(path, legacy=legacy)
ncorr = g_header.ncorr
tvals = g_header.tvals
nnoise = g_header.nnoise
nnoise_trunc = isnothing(nnoise_trunc) ? nnoise : min(nnoise, nnoise_trunc)
fsize = filesize(path)
datsize = 4 + sum(getfield.(c_header, :dsize)) * tvals * nnoise #data_size / ncnfg
ncfg = div(fsize - g_header.hsize - sum(getfield.(c_header, :hsize)), datsize) #(total size - header_size) / data_size
corr_match = findall(x-> (x.type1,x.type2) in type, c_header)
seek(data, g_header.hsize + sum(c.hsize for c in c_header))
res = Array{juobs.CData}(undef, length(corr_match))
data_re = Array{Float64}(undef, length(corr_match), ncfg, tvals)
data_im = zeros(length(corr_match), ncfg, tvals)
vcfg = Array{Int32}(undef, ncfg)
for icfg = 1:ncfg
vcfg[icfg] = read(data, Int32)
c=1
for k = 1:ncorr
if k in corr_match
if c_header[k].is_real == 1
tmp = Array{Float64}(undef, tvals*nnoise)
read!(data, tmp)
tmp2 = reshape(tmp, (nnoise, tvals))
tmp2 = mean(tmp2[1:nnoise_trunc, :], dims=1)
data_re[c, icfg, :] = tmp2[1, :]
elseif c_header[k].is_real == 0
tmp = Array{Float64}(undef, 2*tvals*nnoise)
read!(data, tmp)
tmp2 = reshape(tmp, (2, nnoise, tvals))
tmp2 = mean(tmp2[:, 1:nnoise_trunc, :], dims=2)
data_re[c, icfg, :] = tmp2[1, 1, :]
data_im[c, icfg, :] = tmp2[2, 1, :]
end
c += 1
else
seek(data, position(data) + c_header[k].dsize*tvals*nnoise)
end
end
end
for c in eachindex(corr_match)
res[c] = juobs.CData(c_header[corr_match[c]], vcfg, data_re[c, :, :], data_im[c, :, :], id)
end
close(data)
return res
end
function dev_read_mesons(path::Vector{String}, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,
nnoise_trunc::Union{Int64, Nothing}=nothing)
res = [dev_read_mesons(p, gamma, id=id, legacy=legacy, nnoise_trunc=nnoise_trunc) for p in path]
nrep = length(res)
ncorr = length(res[1])
cdata = Vector{Vector{juobs.CData}}(undef, ncorr)
for icorr = 1:ncorr
cdata[icorr] = Vector{juobs.CData}(undef, nrep)
for r = 1:nrep
cdata[icorr][r] = res[r][icorr]
end
end
return cdata
end
function dev_read_mesons_correction(path::String, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,
nnoise_trunc::Union{Int64, Nothing}=nothing)
type = Vector{Tuple{Int64,Int64}}(undef,length(gamma))
for i in eachindex(gamma)
G1,G2 = gamma[i]
type[i] = (findfirst(x-> x==G1, juobs.gamma_name) - 1, findfirst(x-> x==G2, juobs.gamma_name) - 1)
end
if isnothing(id)
bname = basename(path)
m = findfirst(r"[A-Z][0-9]{3}r[0-9]{3}", bname)
id = bname[m[1:4]]
#id = parse(Int64, bname[m[2:4]])
end
data = open(path, "r")
g_header = juobs.read_GHeader(path)
c_header = juobs.read_CHeader(path, legacy=legacy)
ncorr = g_header.ncorr
tvals = g_header.tvals
nnoise = g_header.nnoise
nnoise_trunc = isnothing(nnoise_trunc) ? nnoise : min(nnoise, nnoise_trunc)
fsize = filesize(path)
datsize = 4 + sum(getfield.(c_header, :dsize)) * tvals * nnoise #data_size / ncnfg
ncfg = div(fsize - g_header.hsize - sum(getfield.(c_header, :hsize)), datsize) #(total size - header_size) / data_size
corr_match = findall(x-> (x.type1,x.type2) in type, c_header)
seek(data, g_header.hsize + sum(c.hsize for c in c_header))
res = Array{juobs.CData}(undef, div(length(corr_match), 2)) # Modification: total length is divided by 2
data_re = zeros(div(length(corr_match), 2), ncfg, tvals) # Modification: total length is divided by 2
data_im = zeros(div(length(corr_match), 2), ncfg, tvals) # Modification: total length is divided by 2
vcfg = Array{Int32}(undef, ncfg)
for icfg = 1:ncfg
vcfg[icfg] = read(data, Int32)
c = 1
sgn = +1 # sign it changes at ncorr / 2. O_exact - O_sloppy
for k = 1:ncorr
if k in corr_match
if c_header[k].is_real == 1
tmp = Array{Float64}(undef, tvals*nnoise)
read!(data, tmp)
tmp2 = reshape(tmp, (nnoise, tvals))
tmp2 = mean(tmp2[1:nnoise_trunc, :], dims=1)
data_re[c, icfg, :] = data_re[c, icfg, :] + sgn * tmp2[1, :]
elseif c_header[k].is_real == 0
tmp = Array{Float64}(undef, 2*tvals*nnoise)
read!(data, tmp)
tmp2 = reshape(tmp, (2, nnoise, tvals))
tmp2 = mean(tmp2[:, 1:nnoise_trunc, :], dims=2)
data_re[c, icfg, :] = data_re[c, icfg, :] + sgn * tmp2[1, 1, :]
data_im[c, icfg, :] = data_im[c, icfg, :] + sgn * tmp2[2, 1, :]
end
c += 1
else
seek(data, position(data) + c_header[k].dsize*tvals*nnoise)
end
if k == div(ncorr, 2)
c = 1
sgn = -1
end
end
end
for c = 1:div(length(corr_match), 2)
res[c] = juobs.CData(c_header[corr_match[c]], vcfg, data_re[c, :, :], data_im[c, :, :], id)
end
close(data)
return res
end
function dev_read_mesons_correction(path::Vector{String}, gamma::Vector{Tuple{String,String}}; id::Union{String, Nothing}=nothing, legacy::Bool=false,
nnoise_trunc::Union{Int64, Nothing}=nothing)
res = [dev_read_mesons_correction(p, gamma, id=id, legacy=legacy, nnoise_trunc=nnoise_trunc) for p in path]
nrep = length(res)
ncorr = length(res[1])
cdata = Vector{Vector{juobs.CData}}(undef, ncorr)
for icorr = 1:ncorr
cdata[icorr] = Vector{juobs.CData}(undef, nrep)
for r = 1:nrep
cdata[icorr][r] = res[r][icorr]
end
end
return cdata
end
\ No newline at end of file
src/juobs_tools.jl
View file @ 8aa2abdb
...@@ -685,10 +685,10 @@ function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, wpm::Union{Dict{Int64 ...@@ -685,10 +685,10 @@ function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, wpm::Union{Dict{Int64
end end
plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, W::Vector{Float64},wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) = sum(W .* obs[plat[1]:plat[2]]) / sum(W) plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, W::Vector{Float64},wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) = sum(W .* obs[plat[1]:plat[2]]) / sum(W)
function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, W::Matrix{Float64},wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
o = obs[plat[1]:plat[2]]
return sum([o[i]*W[i,j]*o[j] for i in axes(W,1), j in axes(W,2)])/sum([o[i]*W[i,j] for i in axes(W,1), j in axes(W,2)])
end
function model_av(fun::Function, y::Vector{uwreal}, guess::Float64; function model_av(fun::Function, y::Vector{uwreal}, guess::Float64;
tm::Vector{Int64}, tM::Vector{Int64}, k::Int64, tm::Vector{Int64}, tM::Vector{Int64}, k::Int64,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing) wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
...@@ -827,7 +827,6 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, ...@@ -827,7 +827,6 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
Ncut = total - length(x) Ncut = total - length(x)
dy = err.(yy) dy = err.(yy)
W = 1 ./ dy .^2 W = 1 ./ dy .^2
p00 = [0.5 for i in 1:1:k] p00 = [0.5 for i in 1:1:k]
chisq = gen_chisq(fun,x,dy) chisq = gen_chisq(fun,x,dy)
fit = curve_fit(fun,x,value.(yy),W,p00) fit = curve_fit(fun,x,value.(yy),W,p00)
...@@ -835,19 +834,12 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, ...@@ -835,19 +834,12 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
isnothing(wpm) ? uwerr(up[1]) : uwerr(up[1],wpm) isnothing(wpm) ? uwerr(up[1]) : uwerr(up[1],wpm)
chi2 = sum(fit.resid.^2) # * dof(fit) / chi_exp chi2 = sum(fit.resid.^2) # * dof(fit) / chi_exp
Q = pvalue(chisq, sum(fit.resid.^2), value.(up), yy, wpm=wpm, W=W, nmc=10000) Q = pvalue(chisq, sum(fit.resid.^2), value.(up), yy, wpm=wpm, W=W, nmc=10000)
# if chi2 / dof(fit) > 5. || abs(value(up[1])) > 2. || err(up[1])/value(up[1])*100 >= 10.
# error()
# end
# println(up[1], " ", up[2], " ", up[3])
push!(Pval, Q) push!(Pval, Q)
# push!(AIC, chi2 + 2*k + 2*Ncut) # AIC criteria
push!(AIC, chi2 - 2*chi_exp) # TIC criteria push!(AIC, chi2 - 2*chi_exp) # TIC criteria
push!(chi2chi2exp, chi2 / dof(fit)) #substitute chi_exp with dof(fit) for chi^2/dof push!(chi2chi2exp, chi2 / dof(fit)) #substitute chi_exp with dof(fit) for chi^2/dof
push!(p1, up[1]) push!(p1, up[1])
push!(mods,string("[", INDEX+1, ",", j, "]")) push!(mods,string("[", INDEX+1, ",", j, "]"))
# println("chi2_vs_exp " , chi2 / dof(fit), " chi2/dof ", sum(fit.resid.^2)/dof(fit))
catch catch
@warn string(":/ Negative window for error propagation at tmin = ", INDEX, ", tmax = ", j, "; skipping that point") @warn string(":/ Negative window for error propagation at tmin = ", INDEX, ", tmax = ", j, "; skipping that point")
end end
...@@ -1248,13 +1240,13 @@ function get_chi(model::Function,x,par,data,W::Vector{Float64}) ...@@ -1248,13 +1240,13 @@ function get_chi(model::Function,x,par,data,W::Vector{Float64})
return sum((data.-model(x,par)).^2 .* W) return sum((data.-model(x,par)).^2 .* W)
end end
@doc raw""" @doc """
fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, npar::Int64; fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, npar::Int64;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}}=Dict{Int64,Vector{Float64}}(), wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}}=Dict{Int64,Vector{Float64}}(),
corr::Bool = true, corr::Bool = true,
W::VecOrMat{Float64}=Vector{Float64}(), W::VecOrMat{Float64}=Vector{Float64}(),
guess::Vector{Float64} = fill(0.5,npar), guess::Vector{Float64} = fill(0.5,npar),
logfile logfile,
C::AbstractMatrix{Float64} = Matrix{Float64}()); C::AbstractMatrix{Float64} = Matrix{Float64}());
It executes a fit of the `ydata` with the fit function `model`. It executes a fit of the `ydata` with the fit function `model`.
...@@ -1283,14 +1275,14 @@ It returns a NamedTuple with names: ...@@ -1283,14 +1275,14 @@ It returns a NamedTuple with names:
""" """
function fit_routine(model::Function, function fit_routine(model::Function,
xdata::AbstractArray{<:Real}, xdata::AbstractArray{<:Real},
ydata::Vector{uwreal}, ydata::AbstractArray{uwreal},
npar::Int64; npar::Int64;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}}=Dict{Int64,Vector{Float64}}(), wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}}=Dict{Int64,Vector{Float64}}(),
corr::Bool = true, corr::Bool = true,
W::AbstractVecOrMat{Float64}=Vector{Float64}(), W::AbstractArray{Float64}=Vector{Float64}(),
guess::AbstractVector{Float64} = fill(0.5,npar), guess::AbstractVector{Float64} = fill(0.5,npar),
logfile = nothing, logfile = nothing,
C::AbstractMatrix{Float64} = Matrix{Float64}(undef,0,0)); C::AbstractArray{Float64} = Matrix{Float64}(undef,0,0));
nobs = length(ydata) nobs = length(ydata)
if length(xdata)!=nobs if length(xdata)!=nobs
...@@ -1324,7 +1316,10 @@ function fit_routine(model::Function, ...@@ -1324,7 +1316,10 @@ function fit_routine(model::Function,
uwerr.(par); uwerr.(par);
println(logfile, "juobs.fit_routine output:") println(logfile, "juobs.fit_routine output:")
println(logfile, "\t\tFit routine in [$(xdata[1]),...,$(xdata[end])]") println(logfile, "\t\tFit routine in [$(xdata[1]),...,$(xdata[end])]")
println(logfile, "\t\tparameters :", par...) println(logfile, "\t\tparameters :")
for p in par
println(logfile,"\t\t\t",p)
end
println(logfile, "\t\tchi2: ", chi2) println(logfile, "\t\tchi2: ", chi2)
println(logfile, "\t\tchiexp: ", chiexp) println(logfile, "\t\tchiexp: ", chiexp)
println(logfile, "\t\tpvalue: ", pval) println(logfile, "\t\tpvalue: ", pval)
...@@ -1333,510 +1328,172 @@ function fit_routine(model::Function, ...@@ -1333,510 +1328,172 @@ function fit_routine(model::Function,
return (par=par,chi2=chi2,chiexp=chiexp,pval=pval) return (par=par,chi2=chi2,chiexp=chiexp,pval=pval)
end end
#= function fit_routine(model::Vector{Function}, xdata::Vector{Array{Float64, N}} where N, ydata::Vector{Array{uwreal, N}} where N, param::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, """
correlated_fit::Bool=false) check_size(M,n)
if !(length(model) == length(xdata) == length(ydata))
error("Dimension mismatch")
end
N = length(model)
dat = ydata[1] Check if `M` has sizes all equal to n or, if is empty
idx = Vector{Vector{Int64}}(undef, N) """
e = Vector{Vector{Float64}}(undef, N) function check_sizes(M,n)
j = 1 if length(M) ==0 || all(size(M) .== n)
for i = 1:N return true
if isnothing(wpm)
uwerr.(ydata[i])
else
# println(yaux for yaux in ydata[i])
for yaux in ydata[i]
[uwerr(yaux[k], wpm) for k in 1:length(yaux)]
end
# [[uwerr(yaux[k], wpm) for k in eachindex(yaux)] for yaux in ydata[i]]
end end
return false
end
e[i] = err.(ydata[i]) function fit_routine(model::Function,
if i > 1 xdata::AbstractArray{ADerrors.uwreal},
dat = vcat(dat, ydata[i]) ydata::AbstractArray{ADerrors.uwreal},
end npar::Int64;
stp = j + length(ydata[i]) - 1 W::AbstractArray = Float64[],
idx[i] = collect(j:stp) wpm::AbstractDict = Dict{Int64,Vector{Float64}}(),
j = stp + 1 corr::Bool = true,
end guess::AbstractArray = fill(0.5,npar),
if !correlated_fit logfile = nothing,
chisq = (par, dat) -> sum([sum((dat[idx[i]] .- model[i](xdata[i], par)).^2 ./e[i].^2) for i=1:N]) C::AbstractArray = Float64[])
else Nx = length(xdata)
inv_cov_tot = juobs.inv_covar(vcat(ydata...)) #[inv_covar_multi_id(dat[idx[ii]]) for ii in 1:N] if (Nx != length(ydata))
function chisq(par, dat) error("[juobs] Error: xdata and ydata must have the same dimension")
res = Vector{uwreal}(undef,0)
for ii in 1:N
res = vcat(res, dat[idx[ii]] .- model[ii](xdata[ii],par))
end end
chi = reshape(res, 1,:) * inv_cov_tot*res Ndata = 2Nx; ## in total we have 2*lengt(xdata) data points
return chi[1] data = vcat(ydata,xdata)
if !check_sizes(W,Ndata)
error("[juobs] Error: W does not have the correct sizes")
end end
if !check_sizes(C,Ndata)
error("[juobs] Error: C does not have the correct sizes")
end end
min_fun(t) = chisq(t, value.(dat)) if corr
p = fill(0.5, param) C = length(C) ==0 ? GammaMethod.cov_AD(data) : C
#println(chisq(p,dat)) W = length(W) ==0 ? LinearAlgebra.pinv(W) : W
sol = optimize(min_fun, p, LBFGS())
#println(chisq(sol.minimizer,dat))
if !correlated_fit
@info("Uncorrelated fit ")
(upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq, sol.minimizer, dat) : fit_error(chisq, sol.minimizer, dat, wpm)
else else
@info("Correlated fit ") W = length(W) ==0 ? [1/yy.err^2 for yy in data] : W
#inv_cov_merge = zeros(size(inv_cov_tot[1],1)*N, size(inv_cov_tot[1],1)*N )
#count = 1
#for k in 1:N
# dim_tmp = size(inv_cov_tot[k], 1)
# inv_cov_merge[count:dim_tmp+count-1, count:dim_tmp+count-1] = inv_cov_tot[k]
# count +=dim_tmp
#end
(upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq, sol.minimizer, dat, W=inv_cov_tot) : fit_error(chisq, sol.minimizer, dat, wpm, W=inv_cov_tot)
end end
println("Chisq / chiexp: ", min_fun(sol.minimizer), " / ", chi2_exp, " (dof: ", length(dat) - param,")") f(q,p) = [i<=Nx ? model(q[i],p) : q[i-Nx] for i in 1:Ndata]
chis2_corrected = (length(dat) - param) * min_fun(sol.minimizer) / chi2_exp fit_func(x,p) = f(view(p,npar+1:npar+Nx),view(p,1:npar))
println("Chisq corrected: ", chis2_corrected) new_guess = [guess; value.(xdata)]
fit = fit_routine(fit_func,zeros(Ndata),data,npar+Nx,
guess=new_guess, wpm=wpm,corr=corr,
W=W,C=C,logfile=nothing)
@info("Returns: params, chi2, chi2_exp\n") if !isnothing(logfile)
return upar, min_fun(sol.minimizer) , chi2_exp uwerr.(fit.par)
println(logfile, "juobs.fit_routine output:")
println(logfile, "\t\txdata carries uncertainty")
println(logfile, "\t\tFit routine in [$(xdata[1].mean),...,$(xdata[end].mean)]")
println(logfile, "\t\tFit parameters:")
for i in 1:npar
println(logfile,"\t\t\t",fit.par[i])
end
println(logfile, "\t\tauxiliary fit parameter:")
println(logfile, "\t\t\txdata[i] \t q[i]")
for i in 1:Nx
println(logfile, "\t\t\t",xdata[i]," \t ", fit.par[npar+i])
end
println(logfile, "\t\tchi2: ", fit.chi2)
println(logfile, "\t\tchiexp: ", fit.chiexp)
println(logfile, "\t\tpvalue: ", fit.pval)
end
return fit
end end
function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; function fit_routine(models::AbstractArray{Function},
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, covar::Bool=false) xdata::AbstractArray{<:AbstractArray},
ydata::AbstractArray{<:AbstractArray{ADerrors.uwreal}},
Nalpha = size(xdata, 2) # number of x-variables npar::AbstractArray{Int64};
Ndata = size(xdata, 1) # number of datapoints W::AbstractArray = Float64[],
if isnothing(wpm) C::AbstractArray = Float64[],
uwerr.(ydata) wpm::AbstractDict = Dict{Int64,Vector{Float64}}(),
uwerr.(xdata) corr::Bool = true,
else guess::AbstractArray{<:AbstractArray} = [fill(0.5,npar[i]) for i in eachindex(models)],
[uwerr(yaux, wpm) for yaux in ydata] logfile = nothing,
[uwerr(xaux, wpm) for xaux in xdata] )
end Nmodel = length(models)
if Nmodel != length(xdata) != length(ydata)
yval = value.(ydata) error("[juobs] Error: you need the same number of model, xdata and 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 end
if !(length(W) == 0 || length(W)==Nmodel)
C = isnothing(wpm) ? [ADerrors.cov(aux[k]) for k = 1:Ndata] : [ADerrors.cov(aux[k], wpm) for k = 1:Ndata] error("[juobs] Error: You need to pass a W matrix for model or none")
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]), method=LBFGS())
(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,")")
chis2_corrected = (length(ydata) - param) * min_fun_cov(sol.minimizer) / chi2_exp
println("Chisq corrected: ", chis2_corrected)
return upar, min_fun_cov(sol.minimizer) / chi2_exp, chi2_exp
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]), method=LBFGS())
(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,")")
chis2_corrected = (length(ydata) - param) * min_fun(sol.minimizer) / chi2_exp
println("Chisq corrected: ", chis2_corrected)
return upar, min_fun(sol.minimizer) / chi2_exp, chi2_exp
end end
if !(length(C) == 0 || length(C)==Nmodel)
#### chisq_full, min_fun out of conditional -> error("[juobs] Error: You need to pass a C matrix for model or none")
#### 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
end
function global_fit_routine(models::Vector{Function}, xdata::Vector{Array{Float64, 1}} , ydata::Vector{Array{uwreal, 1}}, n_par::Int64; guess_param::Vector{Float64}=fill(0.5, n_par), log::Bool=true,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
inv_cov::Union{Matrix{Float64}, Nothing}=nothing)
if !(length(models) == length(ydata))
error("Dimension mismatch")
end
N_sets = length(models) # number of datasets
x_flat = xdata[1] # unrolling xdata
y_flat = ydata[1] # unrolling ydata
chunk_idx = Vector{Vector{Int64}}(undef, N_sets)
j = 1
for i = 1:N_sets
isnothing(wpm) ? uwerr.(ydata[i]) : [uwerr(yaux, wpm) for yaux in ydata[i]]
if i > 1
y_flat = vcat(y_flat, ydata[i])
x_flat =vcat(x_flat, xdata[i])
end
stp = j + length(ydata[i]) - 1
chunk_idx[i] = collect(j:stp)
j = stp + 1
end
function flat_model(x_flat::Vector{Float64}, param::Vector{Float64})
fx = similar(x_flat)
for ii in 1:N_sets
mod = models[ii]
fx[chunk_idx[ii]] = mod(x_flat[chunk_idx[ii]], param)
end
return fx
end
function flat_chisq(param, data)
if isnothing(inv_cov)
chi = 0.0
for ii in 1:N_sets
mod = models[ii]
chi += sum((data[chunk_idx[ii]] .- mod(x_flat[chunk_idx[ii]], param)).^2 ./ err.(y_flat[chunk_idx[ii]]).^2 )
end
return chi
else
# take pieces of covariance and build correlated chi2 piece by piec
# chi =0.0
# for ii in 1:N_sets
# mod = models[ii]
# cov_inv = juovs.inv_covar(ydata[ii])
# res = ydata[ii] .- mod(xdata[ii], param)
# chi += (reshape(res, 1,:) * cov_inv * res)[1]
# end
# return chi
# take whole covariance and multiply for whole vector of dat-mod(x,param)
res = Vector{uwreal}(undef, 0)
for ii in 1:N_sets
mod = models[ii]
res =vcat(res, data[chunk_idx[ii]] .- mod(x_flat[chunk_idx[ii]], param))
end end
chi = reshape(res, 1,:) * inv_cov * res if length(npar) != Nmodel
return chi[1] error("[juobs] Error: You need to specify the number of parameter for every model")
end end
if length(guess)!=Nmodel
error("[juobs] Error: You need to specify the inital guess for all model or for none")
end end
ndata = length.(xdata)
if isnothing(inv_cov) if !all(length.(ydata).== ndata)
@info("Uncorrelated fit") error("[juobs] Error: Mismatch in xdata and ydata. Make sure that the data corresponds")
fit = curve_fit(flat_model, x_flat, value.(y_flat), 1.0 ./ err.(y_flat).^2, guess_param )
(upar, chi2_exp) = isnothing(wpm) ? fit_error(flat_chisq, coef(fit) , y_flat) : fit_error(flat_chisq, coef(fit) , y_flat, wpm)
else
@info("Correlated fit")
fit = curve_fit(flat_model, x_flat, value.(y_flat), inv_cov, guess_param )
(upar, chi2_exp) = isnothing(wpm) ? fit_error(flat_chisq, coef(fit) , y_flat, W=inv_cov) : fit_error(flat_chisq, coef(fit) , y_flat, wpm, W=inv_cov)
end
println("chi2 from fit residue: ", sum(fit.resid.^2))
chisq = flat_chisq(coef(fit), y_flat)
if log
println("Converged param: ", fit.converged)
println("Chisq / chiexp: ", chisq, " / ", chi2_exp, " (dof: ", dof(fit),")")
chis2_corrected = (dof(fit)) * chisq / chi2_exp
println("Chisq corrected: ", chis2_corrected)
end
return upar, value(chisq), chi2_exp, dof(fit)
end
function fit_alg(f::Function, x::Union{Vector{Int64}, Vector{Float64}, Matrix{Float64}}, y::Vector{uwreal},
n::Int64, guess::Union{Float64, Vector{Float64}, Nothing}=nothing;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
isnothing(wpm) ? uwerr.(y) : [uwerr(y[i], wpm) for i in 1:length(y)]
W = 1 ./ err.(y) .^ 2
chisq = fit_defs(f,x,W)
#lb = [-Inf for i in 1:n]
#ub = [+Inf for i in 1:n]
if guess == nothing
p0 = [.5 for i in 1:n]
else
p0 = [guess; [1. for i in 1:n-length(guess)]]
#lb[1] = .9 * guess
#ub[1] = 1.1 * guess
end
fit = curve_fit(f,x,value.(y),W,p0)#,lower=lb,upper=ub)
chi2 = sum(fit.resid .^ 2)
isnothing(wpm) ? (up,chi_exp) = fit_error(chisq,coef(fit),y) : (up,chi_exp) = fit_error(chisq,coef(fit),y,wpm)
isnothing(wpm) ? pval = pvalue(chisq,chi2,value.(up),y) : pval = pvalue(chisq,chi2,value.(up),y,wpm=wpm)
return up, chi2, chi_exp, pval
end
function fit_alg(model::Function,x::Array{Float64},y::Array{uwreal},param::Int64,W::Matrix{Float64};
guess::Union{Float64, Vector{Float64}, Nothing}=nothing,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
if guess == nothing
p00 = [.5 for i in 1:param]
else
p00 = [guess; [1. for i in 1:param-length(guess)]]
#lb[1] = .9 * guess
#ub[1] = 1.1 * guess
end end
Ntotal = sum(ndata);
X = vcat(xdata...)
Y = vcat(ydata...)
chisq_corr(par,dat) = juobs.gen_chisq_correlated(model, x, W, par, dat) function make_matrix(M)
fit = curve_fit(model,x,value.(y),W,p00) aux = zeros(Ntotal,Ntotal)
up, chi_exp = fit_error(chisq_corr, coef(fit), y, wpm, W=W) ie = 0
uwerr.(up) for m in 1:Nmodel
chi2 = sum(fit.resid .^ 2) is = ie+1
pval = pvalue(chisq_corr, sum(fit.resid .^ 2), value.(up), y, W) ie += ndata[m];
doff = dof(fit) aux[is:ie,is:ie] .= M[m]
return up, chi2, chi_exp, pval, doff
end
function fit_alg(f::Vector{Function}, x::Vector{Matrix{Float64}}, y::Vector{Vector{uwreal}},
n::Int64, guess::Union{Float64, Vector{Float64}, Nothing}=nothing;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
y_n = y[1]
x_n = x[1]
for i in 2:length(x)
y_n = vcat(y_n, y[i])
x_n = vcat(x_n, x[i])
end end
idx = [collect(1:size(x[i],1)) for i in 1:length(x)] return aux
[idx[i] = idx[i] .+ idx[i-1][end] for i in 2:length(idx)]
isnothing(wpm) ? [uwerr.(y[i]) for i in 1:length(y)] : [[uwerr(y[i][j], wpm) for j in 1:length(y[i])] for i in 1:length(y)]
W = [1 ./ err.(y[i]) .^ 2 for i in 1:length(y)]
chisq = (par, y_n) -> sum([sum((y_n[idx[i]] .- f[i](x[i], par)) .^ 2 .* W[i]) for i in 1:length(x)])
min_fun(t) = chisq(t, value.(y_n))
if guess == nothing
p0 = [.5 for i in 1:n]
else
p0 = [guess; [1. for i in 1:n-length(guess)]]
end end
sol = optimize(min_fun, p0, Optim.Options(g_tol=1e-8, iterations=1000000))
chi2 = min_fun(sol.minimizer)
isnothing(wpm) ? (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n) : (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n,wpm)
isnothing(wpm) ? pval = pvalue(chisq,chi2,value.(up),y_n) : pval = pvalue(chisq,chi2,value.(up),y_n,wpm=wpm)
return up, chi2, chi_exp, pval
end
##CHECK THIS ONE, NOT WORKING GOOD WW = length(W)==0 ? W : make_matrix(W);
function fit_alg(f::Vector{Function}, x::Vector{Matrix{Float64}}, y::Vector{Vector{uwreal}}, W::Vector{Matrix{Float64}}, CC = length(C)==0 ? C : make_matrix(C);
n::Int64, guess::Union{Float64, Vector{Float64}, Nothing}=nothing;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
y_n = y[1] function Model(x,p)
x_n = x[1] res = zeros(Nmodel)
for i in 2:length(x) res[1] = models[1](x[1:ndata[1]],p[1:npar[1]])
y_n = vcat(y_n, y[i]) for i in 2:Nmodel
x_n = vcat(x_n, x[i]) rd = (ndata[i-1]+1):ndata[i]
rp = (npar[i-1]+1):npar[i]
res[i] = models[i](x[rd],p[rp])
end end
idx = [collect(1:size(x[i],1)) for i in 1:length(x)] return res
[idx[i] = idx[i] .+ idx[i-1][end] for i in 2:length(idx)]
isnothing(wpm) ? [uwerr.(y[i]) for i in 1:length(y)] : [[uwerr(y[i][j], wpm) for j in 1:length(y[i])] for i in 1:length(y)]
chisq = (par, y_n) -> sum([juobs.gen_chisq_correlated(f[i], x[i], W[i], par, y_n[i]) for i in 1:length(x)])
min_fun(t) = chisq(t, value.(y_n))
if guess == nothing
p0 = [.5 for i in 1:n]
else
p0 = [guess; [1. for i in 1:n-length(guess)]]
end end
sol = optimize(min_fun, p0, Optim.Options(g_tol=1e-8, iterations=1000000)) Guess = vcat(guess...)
chi2 = min_fun(sol.minimizer) fit = fit_routine(Model,X,Y,sum(npar),guess = Guess,W=WW,C=CC,corr=corr,logfile=nothing,wpm=wpm)
W_aux = Matrix{Float64}(undef, length(y_n), length(y_n)) if !isnothing(logfile)
for i in 1:length(y_n) uwerr.(fit.par)
for j in 1:length(y_n) flag = typeof(xdata) <: AbstractArray{<:AbstractArray{ADerrors.uwreal}}
if i <= length(y[1]) println(logfile, "juobs.fit_routine output:")
if j <= length(y[1]) if flag
W_aux[i,j] = W[1][i,j] println(logfile, "\t\txdata carries uncertainty")
else
W_aux[i,j] = .0
end
else
if j <= length(y[1])
W_aux[i,j] = .0
else
W_aux[i,j] = W[2][i-length(y[1]),j-length(y[1])]
end
end end
println(logfile, "\t\tGlobal fit with $Nmodel functions")
for nm in 1:Nmodel
println(logfile,"---Model $(nm)---")
println(logfile, "\t\tFit routine in [$(xdata[nm][1].mean),...,$(xdata[nm][end].mean)]")
println(logfile, "\t\tFit parameters:")
for i in 1:npar[nm]
off = nm ==1 ? 0 : npar[nm-1]
println(logfile,"\t\t\t",fit.par[off+i])
end end
end end
if flag
isnothing(wpm) ? (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n,W=W_aux) : (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n,wpm,W=W_aux) pritnln(logfile, "\t\tauxiliary fit parameter:")
isnothing(wpm) ? pval = pvalue(chisq,chi2,value.(up),y_n,W_aux) : pval = pvalue(chisq,chi2,value.(up),y_n,wpm=wpm,W_aux) println(logfile, "\t\t\txdata[m][i] \t q[i]")
return up, chi2, chi_exp, pval for nm in 1:Nmodel, i in 1:npar[nm]
end off = nm==1 ? 0 : npar[end]+npar[i-1]
println(logfile, "\t\t\t",xdata[nm][i]," \t ", fit.par[off+i])
function fit_alg(f::Vector{Function}, x::Vector{Matrix{Float64}}, y::Vector{Vector{uwreal}},
n::Int64, lb::Vector{Float64}, ub::Vector{Float64}, guess::Union{Float64, Vector{Float64}, Nothing}=nothing;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
y_n = y[1]
x_n = x[1]
for i in 2:length(x)
y_n = vcat(y_n, y[i])
x_n = vcat(x_n, x[i])
end end
idx = [collect(1:size(x[i],1)) for i in 1:length(x)]
[idx[i] = idx[i] .+ idx[i-1][end] for i in 2:length(idx)]
isnothing(wpm) ? [uwerr.(y[i]) for i in 1:length(y)] : [[uwerr(y[i][j], wpm) for j in 1:length(y[i])] for i in 1:length(y)]
W = [1 ./ err.(y[i]) .^ 2 for i in 1:length(y)]
chisq = (par, y_n) -> sum([sum((y_n[idx[i]] .- f[i](x[i], par)) .^ 2 .* W[i]) for i in 1:length(x)])
min_fun(t) = chisq(t, value.(y_n))
if guess == nothing
p0 = [.5 for i in 1:n]
else
p0 = [guess; [1. for i in 1:n-length(guess)]]
end end
sol = optimize(min_fun, lb, ub, p0, Fminbox(NelderMead()), Optim.Options(g_tol=1e-8, iterations=1000000)) println(logfile, "\t\tchi2: ", chi2)
chi2 = min_fun(sol.minimizer) println(logfile, "\t\tchiexp: ", chiexp)
isnothing(wpm) ? (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n) : (up,chi_exp) = fit_error(chisq,sol.minimizer,y_n,wpm) println(logfile, "\t\tpvalue: ", pval)
isnothing(wpm) ? pval = pvalue(chisq,chi2,value.(up),y_n) : pval = pvalue(chisq,chi2,value.(up),y_n,wpm=wpm)
return up, chi2, chi_exp, pval
end
function fit_defs(f::Function,x,W) ## uncorrelated fit
chisq(p,d) = sum((d .- f(x,p)).^2 .* W)
return chisq
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 gen_chisq_correlated(f::Function, x::Array{<:Real}, covar_inv::Matrix, par, dat )
chi2 = 0.0
#delta = dat - ffit_routine(x, par)
#chi2 = (reshape(delta, 1,:) * covar_inv * delta )[1]
for i in eachindex(dat)
for j in eachindex(dat)
#println(f(x,par)[i])
#println(f(x,par)[j])
#println(f(x[j,:],par))
#println(f(x[i,:],par))
chi2 = chi2 + (dat[i] - f(x,par)[i]) * covar_inv[i,j] * (dat[j] - f(x,par)[j])
end
end
#chisq(par, dat) = sum(
# transpose((dat .- f(x,par))) * covar_inv * (dat .- f(x,par))
#)
return chi2
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 end
return chi2 return fit
end 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
function get_covar(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
isnothing(wpm) ? uwerr.(obs) : [uwerr(yaux, wpm) for yaux in obs]
ens = vcat(ensembles.(obs)...)
tmp_fluc = mchist.(obs, ens)
fluc = mapreduce(permutedims, vcat, tmp_fluc)'
covar = fluc' * fluc ./ size(fluc)[1]^2
return covar
end
function get_covar_multi_id(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
isnothing(wpm) ? uwerr.(obs) : [uwerr(yaux, wpm) for yaux in obs]
id_set = getindex.(ensembles.(obs), 2) # hardcoded assuming that ens label is always the second entry in each ensembles(obs[i])
cov_tot = zeros(length(id_set), length(id_set))
count = 1
for (k,id) in enumerate(unique(id_set))
idx = findall(x-> x==id, id_set )
fluc_aux = mchist.(obs[idx], fill(id, length(idx)))
fluc = mapreduce(permutedims, vcat, fluc_aux)'
covar_set = fluc' * fluc ./ size(fluc)[1]^2
dim_tmp = size(covar_set,1)
cov_tot[count:dim_tmp+count-1, count:dim_tmp+count-1] = covar_set
count +=dim_tmp
end
return cov_tot
end
function inv_covar(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
cov = get_covar(obs, wpm=wpm)
cov_inv = inv(cov)
return (cov_inv' + cov_inv) / 2.
end
function inv_covar_multi_id(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
cov = get_covar_multi_id(obs, wpm=wpm)
cov_inv = inv(cov)
return (cov_inv' + cov_inv) / 2.
end =#
@doc raw""" @doc raw"""
pvalue(chisq::Function, pvalue(chisq::Function,
......
src/juobs_types.jl
View file @ 8aa2abdb
##############UTILITY NAMES##############
##############UTILITY NAMES##############
#= STARTING AT 0 #= STARTING AT 0
noise_name=['Z2', 'GAUSS', 'U1', 'Z4'] noise_name=['Z2', 'GAUSS', 'U1', 'Z4']
gamma_name = ['G0', 'G1', 'G2', 'G3', gamma_name = ['G0', 'G1', 'G2', 'G3',
......
test/test_fit_routines.jl 0 → 100644
View file @ 8aa2abdb
using Revise, juobs, ADerrors
true_model(x) = 0.5 + 0.02*x +0.001*x^2
xdata = collect(0.0:0.1:1)
ydata = let
yaux = true_model.(xdata)
[uwreal(randn(1000)*0.0003.+y,"fit_test") for y in yaux]
end
uwerr.(ydata)
@.model(x,p) = p[1] + p[2]*x + p[3]*x^2
println("true model: 0.5 + 0.02x + 0.0001x^2")
println("fit model: p[1]+ p[2]x + p[3]x^2 ")
fit = fit_routine(model,xdata,ydata,3,logfile=stdout,corr=false);
xdata_s = [uwreal(randn(1000)*0.0003.+x,"xdata") for x in xdata]
fit2 = fit_routine(model,xdata_s,ydata,3,logfile=stdout,corr=false)
@.model2(x,p) = p[1] + (p[2] + p[3]*x)^2
Cy = [y.err^2 for y in ydata]
Cy = [i==j ? Cy[i] : 0 for i in eachindex(ydata), j in eachindex(ydata)]
W = [C==0 ? 0 : 1/C for C in Cy ]
fit3 = fit_routine([model,model2],[xdata,xdata],[ydata,ydata],[3,3],corr=false,logfile=stdout, W=[W,W],C=[Cy,Cy])
true
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment