Commit e3968345 authored by Alberto Ramos's avatar Alberto Ramos

Documentation of all functions complete

parent ff40c6fb
...@@ -802,6 +802,25 @@ An optional parameter `wpm` can be used to choose the summation window for the r ...@@ -802,6 +802,25 @@ An optional parameter `wpm` can be used to choose the summation window for the r
cov(a::Vector{uwreal}) = cov(a::Vector{uwreal}, wsg, empt) cov(a::Vector{uwreal}) = cov(a::Vector{uwreal}, wsg, empt)
cov(a::Vector{uwreal}, wpm::Dict{Int64,Vector{Float64}}) = cov(a::Vector{uwreal}, wsg, wpm::Dict{Int64,Vector{Float64}}) cov(a::Vector{uwreal}, wpm::Dict{Int64,Vector{Float64}}) = cov(a::Vector{uwreal}, wsg, wpm::Dict{Int64,Vector{Float64}})
@doc raw"""
trcov(M::Array{Float64, 2}, a::Vector{uwreal})
Given a vector of `uwreal`, `a[:]` and a two dimensional array `M`, this routine computes ``{\rm tr}(MC)``, where ``C_{ij} = {\rm cov}(a[i], a[j])``.
```@example
using ADerrors, LinearAlgebra # hide
a = uwreal([1.3, 0.01], 1) # 1.3 +/- 0.01
b = uwreal([5.3, 0.23], 2) # 5.3 +/- 0.23
c = uwreal(rand(2000), 3)
x = [a+b+sin(c), a-b+cos(c), c-b/a]
M = [1.0 0.2 0.1
0.2 2.0 0.3
0.1 0.3 1.0]
mcov = cov(x)
d = tr(mcov * M)
println("Better be zero: ", d -trcov(M, x))
"""
trcov(M, a::Vector{uwreal}) = trcov(M, a::Vector{uwreal}, wsg, empt) trcov(M, a::Vector{uwreal}) = trcov(M, a::Vector{uwreal}, wsg, empt)
trcov(M, a::Vector{uwreal}, wpm::Dict{Int64,Vector{Float64}}) = trcov(M, a::Vector{uwreal}, wsg, wpm::Dict{Int64,Vector{Float64}}) trcov(M, a::Vector{uwreal}, wpm::Dict{Int64,Vector{Float64}}) = trcov(M, a::Vector{uwreal}, wsg, wpm::Dict{Int64,Vector{Float64}})
...@@ -813,5 +832,18 @@ trcorr(M, a::Vector{uwreal}, ...@@ -813,5 +832,18 @@ trcorr(M, a::Vector{uwreal},
wpm::Dict{Int64,Vector{Float64}}) = trcorr(M, a::Vector{uwreal}, wsg, wpm::Dict{Int64,Vector{Float64}}, Vector{Float64}()) wpm::Dict{Int64,Vector{Float64}}) = trcorr(M, a::Vector{uwreal}, wsg, wpm::Dict{Int64,Vector{Float64}}, Vector{Float64}())
"""
neid(a::uwreal)
Returns the number of different ensemble ID's contributing to `a`
```@example
using ADerrors # hide
a = uwreal([1.2, 0.2], 12) # a = 1.2 +/- 0.2
b = uwreal([7.2, 0.5], 13) # a = 7.2 +/- 0.5
c = a*b
println("Number od ID contributing to c: ", neid(a))
```
"""
neid(a::uwreal) = ADerrors.unique_ids!(a::uwreal, wsg) neid(a::uwreal) = ADerrors.unique_ids!(a::uwreal, wsg)
...@@ -31,6 +31,7 @@ function err(a::uwreal) ...@@ -31,6 +31,7 @@ function err(a::uwreal)
end end
return a.err return a.err
end end
""" """
value(a::uwreal) value(a::uwreal)
...@@ -43,6 +44,7 @@ println("a has central value: ", value(a)) ...@@ -43,6 +44,7 @@ println("a has central value: ", value(a))
``` ```
""" """
value(a::uwreal) = a.mean value(a::uwreal) = a.mean
""" """
derror(a::uwreal) derror(a::uwreal)
...@@ -83,6 +85,7 @@ uwerr(a) ...@@ -83,6 +85,7 @@ uwerr(a)
println("Error analysis result: ", a, " (tauint = ", taui(a, 666), ")") println("Error analysis result: ", a, " (tauint = ", taui(a, 666), ")")
``` ```
""" """
function taui(a::uwreal, mcid::Int64) function taui(a::uwreal, mcid::Int64)
idx = find_mcid(a, mcid) idx = find_mcid(a, mcid)
if (idx == nothing) if (idx == nothing)
...@@ -115,6 +118,7 @@ println("Error analysis result: ", a, ...@@ -115,6 +118,7 @@ println("Error analysis result: ", a,
" (tauint = ", taui(a, 666), " +/- ", dtaui(a, 666), ")") " (tauint = ", taui(a, 666), " +/- ", dtaui(a, 666), ")")
``` ```
""" """
function dtaui(a::uwreal, mcid::Int64) function dtaui(a::uwreal, mcid::Int64)
idx = find_mcid(a, mcid) idx = find_mcid(a, mcid)
if (idx == nothing) if (idx == nothing)
...@@ -357,8 +361,8 @@ function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String ...@@ -357,8 +361,8 @@ function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String
end end
if (length(a.cfd) > 0) if (length(a.cfd) > 0)
println(a.mean, " +/- ", a.err) println(io, a.mean, " +/- ", a.err)
println(" ## Number of error sources: ", length(a.ids)) println(io, " ## Number of error sources: ", length(a.ids))
n = 0 n = 0
for i in 1:length(a.cfd) for i in 1:length(a.cfd)
...@@ -367,8 +371,8 @@ function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String ...@@ -367,8 +371,8 @@ function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String
n = n + 1 n = n + 1
end end
end end
println(" ## Number of MC ids : ", n) println(io, " ## Number of MC ids : ", n)
println(" ## Contribution to error : Ensemble [%] [MC length]") println(io, " ## Contribution to error : Ensemble [%] [MC length]")
truncate_ascii(s::String,n::Int) = s[1:min(sizeof(s),n)] truncate_ascii(s::String,n::Int) = s[1:min(sizeof(s),n)]
...@@ -383,15 +387,15 @@ function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String ...@@ -383,15 +387,15 @@ function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String
sndt = join(ws.fluc[idx].ivrep, ",") sndt = join(ws.fluc[idx].ivrep, ",")
sid = truncate_ascii(get(names, a.ids[ip[i]], string(a.ids[ip[i]])), ntrunc) sid = truncate_ascii(get(names, a.ids[ip[i]], string(a.ids[ip[i]])), ntrunc)
if (ws.fluc[idx].nd > 1) if (ws.fluc[idx].nd > 1)
Printf.@printf(" # %45s %6.2f %s\n", Printf.@printf(io, " # %45s %6.2f %s\n",
sid, 100.0 .* a.cfd[ip[i]].var ./ a.err^2, sndt) sid, 100.0 .* a.cfd[ip[i]].var ./ a.err^2, sndt)
else else
Printf.@printf(" # %45s %6.2f -\n", Printf.@printf(io, " # %45s %6.2f -\n",
sid, 100.0 .* a.cfd[ip[i]].var ./ a.err^2) sid, 100.0 .* a.cfd[ip[i]].var ./ a.err^2)
end ' end
end end
else else
print(a.mean, " (Error not available... maybe run uwerr)") print(io, a.mean, " (Error not available... maybe run uwerr)")
end end
end end
...@@ -428,13 +432,13 @@ read_uwreal(fb) = read_bdio(fb, ADerrors.wsg) ...@@ -428,13 +432,13 @@ read_uwreal(fb) = read_bdio(fb, ADerrors.wsg)
""" """
write_uwreal(p::uwreal, fb, iu::Int) write_uwreal(p::uwreal, fb, iu::Int)
Given a `BDIO` file handler `fb`, this writes the observable `p` in a BDIO resord with user info `iu`. Given a `BDIO` file handler `fb`, this writes the observable `p` in a BDIO record with user info `iu`.
```@example ```@example
using ADerrors # hide using ADerrors # hide
using BDIO using BDIO
a = uwreal(rand(2000), 12) a = uwreal(rand(2000), 12)
# Create a BDIO file and write a with user info 8. # Create a BDIO file and write observable with user info 8.
fb = BDIO_open("/tmp/foo.bdio", "w", "Test file") fb = BDIO_open("/tmp/foo.bdio", "w", "Test file")
write(a, fb, 8) write(a, fb, 8)
BDIO_close(fb) BDIO_close(fb)
......
...@@ -106,7 +106,7 @@ end ...@@ -106,7 +106,7 @@ end
Given a ``\chi^2(p, d)``, function of the fit parameters `p[:]` and the data `d[:]`, compute the expected value of the ``\chi^2(p, d)``. Given a ``\chi^2(p, d)``, function of the fit parameters `p[:]` and the data `d[:]`, compute the expected value of the ``\chi^2(p, d)``.
### Arguments #### Arguments
- `chisq`: Must be a function of two vectors (i.e. `chisq(p::Vector, d::Vector)`). The function is assumed to have the form - `chisq`: Must be a function of two vectors (i.e. `chisq(p::Vector, d::Vector)`). The function is assumed to have the form
...@@ -117,10 +117,38 @@ where the function ``f_i(p)`` is an arbitrary function of the fit parameters. In ...@@ -117,10 +117,38 @@ where the function ``f_i(p)`` is an arbitrary function of the fit parameters. In
- `data`: A vector of `uwreal`. The data whose fluctuations enter in the evaluation of the `chisq`. - `data`: A vector of `uwreal`. The data whose fluctuations enter in the evaluation of the `chisq`.
- `W`: A matrix. The weights that enter in the evaluation of the `chisq` function. If a vector is passed, the matrix is assumed to be diagonal (i.e. **uncorrelated** fit). If no weights are passed, the routines assumes that `W` is diagonal with entries given by the inverse errors squared of the data (i.w. the `chisq` is weighted with the errors of the data). - `W`: A matrix. The weights that enter in the evaluation of the `chisq` function. If a vector is passed, the matrix is assumed to be diagonal (i.e. **uncorrelated** fit). If no weights are passed, the routines assumes that `W` is diagonal with entries given by the inverse errors squared of the data (i.w. the `chisq` is weighted with the errors of the data).
### Example #### Example
```@example ```@example
using ADerrors # hide using ADerrors, Distributions # hide
# Generate correlated samples with average 0.1
npt = 12
sig = zeros(npt, npt)
dx = zeros(npt)
for i in 1:npt
dx[i] = 0.01*i
sig[i,i] = dx[i]^2
for j in i+1:npt
sig[i,j] = 0.0001 - 0.000005*abs(i-j)
sig[j,i] = 0.0001 - 0.000005*abs(i-j)
end
end
dmv = MvNormal([0.1 for n in 1:npt], sig)
vs = rand(dmv, 1)
# Create the uwreal data that we want to
# fit to a constant
dt = cobs(vs[:,1], sig, [100+n for n in 1:npt])
# Define the chi^2
chisq(p, d) = sum( (d .- p[1]) .^ 2 ./ dx .^2 )
# The result of an uncorrelated fit to a
# constant is the weighted average
xp = [sum(value.(dt) ./ dx)/sum(1.0 ./ dx)]
# Compare chi^2 and expected chi^2
println("chi^2 / chi_exp^2: ", chisq(xp, value.(dt)), " / ", chiexp(chisq, xp, dt))
``` ```
""" """
...@@ -168,6 +196,64 @@ function chiexp(chisq::Function, ...@@ -168,6 +196,64 @@ function chiexp(chisq::Function,
return cse return cse
end end
@doc raw"""
fit_error(chisq::Function, xp::Vector{Float64}, data::Vector{uwreal};
W = Vector{Float64}(), chi_exp = true)
Given a ``\chi^2(p, d)``, function of the fit parameters `p[:]` and the data `d[:]`, this routine return the fit parameters as `uwreal` type and optionally, the expected value of ``\chi^2(p, d)``.
#### Arguments
- `chisq`: Must be a function of two vectors (i.e. `chisq(p::Vector, d::Vector)`). To determine the fit parameters, the function can be arbitrary, but for the determination of the expected ``\chi^2(p, d)`` is assumed to have the form
``\chi^2(p, d) = \sum_{ij} [d_i - f_i(p)]W_{ij}[d_j - f_j(p)]``
where the function ``f_i(p)`` is an arbitrary function of the fit parameters. In simple words, the expected ``\chi^2(p, d)`` is determined assuming that the function ``\chi^2(p, d)`` is quadratic in the data.
- `xp`: A vector of `Float64`. The value of the fit parameters at the minima.
- `data`: A vector of `uwreal`. The data whose fluctuations enter in the evaluation of the `chisq`.
- `W`: A matrix. The weights that enter in the evaluation of the `chisq` function. If a vector is passed, the matrix is assumed to be diagonal (i.e. **uncorrelated** fit). If no weights are passed, the routines assumes that `W` is diagonal with entries given by the inverse errors squared of the data (i.e. the `chisq` is weighted with the errors of the data).
- `chi_exp`: Bool type. If false, do not compute the expected ``\chi^2(p, d)``.
#### Example
```@example
using ADerrors, Distributions # hide
# Generate correlated samples with average 0.1
npt = 12
sig = zeros(npt, npt)
dx = zeros(npt)
for i in 1:npt
dx[i] = 0.01*i
sig[i,i] = dx[i]^2
for j in i+1:npt
sig[i,j] = 0.0001 - 0.000005*abs(i-j)
sig[j,i] = 0.0001 - 0.000005*abs(i-j)
end
end
dmv = MvNormal([0.1 for n in 1:npt], sig)
vs = rand(dmv, 1)
# Create the uwreal data that we want to
# fit to a constant
dt = cobs(vs[:,1], sig, [100+n for n in 1:npt])
# Define the chi^2
chisq(p, d) = sum( (d .- p[1]) .^ 2 ./ dx .^2 )
# The result of an uncorrelated fit to a
# constant is the weighted average
xp = [sum(value.(dt) ./ dx)/sum(1.0 ./ dx)]
# Propagate errors to the fit parameters and
# determine the expected chi^2
(fitp, csqexp) = fit_error(chisq, xp, dt)
uwerr.(fitp)
println("Fit parameter: ", fitp[1])
println("chi^2 / chi_exp^2: ", chisq(xp, value.(dt)), " / ", csqexp)
"""
function fit_error(chisq::Function, function fit_error(chisq::Function,
xp::Vector{Float64}, xp::Vector{Float64},
data::Vector{uwreal}; data::Vector{uwreal};
...@@ -190,8 +276,10 @@ function fit_error(chisq::Function, ...@@ -190,8 +276,10 @@ function fit_error(chisq::Function,
hess = Array{Float64}(undef, n+m, n+m) hess = Array{Float64}(undef, n+m, n+m)
ForwardDiff.hessian!(hess, ccsq, xav, cfg) ForwardDiff.hessian!(hess, ccsq, xav, cfg)
hinv = LinearAlgebra.pinv(hess[1:n,1:n]) hm = view(hess, 1:n, 1:n)
grad = - hinv[1:n,1:n] * hess[1:n,n+1:n+m] sm = view(hess, 1:n, n+1:n+m)
hinv = LinearAlgebra.pinv(hm)
grad = - hinv * sm
param = addobs(data, grad, xp) param = addobs(data, grad, xp)
...@@ -222,6 +310,40 @@ function fit_error(chisq::Function, ...@@ -222,6 +310,40 @@ function fit_error(chisq::Function,
return param, cse return param, cse
end end
@doc raw"""
int_error(fint::Function, a, b, p::Vector{uwreal})
Returns an `uwreal` with the values of
``\int_a^b {\rm fint}(x; p)\, {\rm d}x``
and the errors the parameters `p[:]` and (optionally) the limits of the intgeral (`a`, `b`).
```@example
using ADerrors, QuadGK, ForwardDiff
# Average values and covariance from
# https://inspirehep.net/literature/1477411
# here we try to reproduce the result
# Scale ratio = 21.86(42) Eq.5.6
avg = [16.26, 0.12, -0.0038]
Mcov = [0.478071 -0.176116 0.0135305
-0.176116 0.0696489 -0.00554431
0.0135305 -0.00554431 0.000454180]
p = cobs(avg, Mcov, [1, 2001, 32])
g1s = uwreal([2.6723, 0.0064], 4)
g2s = 11.31
fint(x, p) = - (p[1] + p[2]*x^2 + p[3]*x^4)/x^3
g1 = sqrt(g1s)
g2 = sqrt(g2s)
sint = 2.0*exp(-int_error(fint, g1, g2, p))
uwerr(sint)
print(" From integral evaluation: ")
details(sint)
```
"""
function int_error(fint::Function, function int_error(fint::Function,
a, a,
b, b,
......
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