Skip to content

A
ADerrors.jl
Commit fabaca28 authored by Alberto Ramos's avatar Alberto Ramos
Browse files
Options

Adds support to tag ensembles with a String. 

- Native BDIO format also adds this option
parent 6cd322d4
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 482 additions and 228 deletions
Manifest.toml
View file @ fabaca28
This diff is collapsed.
docs/src/tutorial.md
View file @ fabaca28
......@@ -7,7 +7,7 @@ contains variables with uncertainties. In very simple cases this is
just a central value and an error
```@repl gs
using ADerrors
a = uwreal([1.0, 0.1], 1) # 1.0 +/- 0.1
a = uwreal([1.0, 0.1], "Var with error") # 1.0 +/- 0.1
```
But it also supports the case where the uncertainty comes from Monte
Carlo measurements
......@@ -23,10 +23,10 @@ for i in 2:1000
end
end
b = uwreal(x.^2, 200)
c = uwreal(x.^4, 200)
b = uwreal(x.^2, "Random walk ensemble in [-1,1]")
c = uwreal(x.^4, "Random walk ensemble in [-1,1]")
```
Correlations between variables are treated consistently. Each call to `uwreal` contains an ensemble `ID` tag. In the previous examples `a` has been measured on ensemble `ID` 1, while both `b` and `c` have been measured on ensemble `ID` 200. Monte Carlo measurements with the same `ID` tag are assumed to be measured on the same configurations: `b` and `c` above are statistically correlated, while `a` will be uncorrelated with both `b` and `c`.
Correlations between variables are treated consistently. Each call to `uwreal` contains an ensemble `ID` tag. In the previous examples `a` has been measured on ensemble `ID` `Var with error`, while both `b` and `c` have been measured on ensemble `ID` `Random walk ensemble in [-1,1]`. Monte Carlo measurements with the same `ID` tag are assumed to be measured on the same configurations: `b` and `c` above are statistically correlated, while `a` will be uncorrelated with both `b` and `c`.
One can perform operations with `uwreal` variables as if there were normal floats
```@repl gs
......@@ -45,7 +45,7 @@ A call to `uwerr` will apply the ``\Gamma``-method to MC ensembles. `ADerrors.jl
println("Details on variable d: ")
details(d)
```
Here we can see that there are two ensembles contributing to the uncertainty in `d`. We recognize this as `d` being a function of both `a` (measured on ensemble 1), and `b` and `c`, measured on ensemble 200. Most of the error in `d` comes in fact from ensemble 200.
Here we can see that there are two ensembles contributing to the uncertainty in `d`. We recognize this as `d` being a function of both `a` (measured on ensemble `Var with error`), and `b` and `c`, measured on ensemble `Random walk ensemble in [-1,1]`. Most of the error in `d` comes in fact from ensemble `Random walk ensemble in [-1,1]`.
Note that one does not need to use `uwerr` on a variable unless one is interested in the error on that variable. For example, continuing with the previous example
```@repl gs
......@@ -84,8 +84,7 @@ using LsqFit, ADerrors
xdata = range(0, stop=10, length=20);
ydata = Vector{uwreal}(undef, length(xdata));
for i in eachindex(ydata)
# Note that we assign to point i the ID 1000*i
ydata[i] = uwreal([model(xdata[i], [1.0 2.0]) + 0.01*getindex(randn(1),1), 0.01], 1000*i)
ydata[i] = uwreal([model(xdata[i], [1.0 2.0]) + 0.01*getindex(randn(1),1), 0.01], "Point "*string(i))
uwerr(ydata[i]) # We will need the errors for the weights of the fit
end
p0 = [0.5, 0.5];
......@@ -142,9 +141,9 @@ for i in 2:10000
end
# x^2 only measured on odd configurations
x2 = uwreal(x[1:2:9999].^2, 1001, collect(1:2:9999), 10000)
x2 = uwreal(x[1:2:9999].^2, 1001, collect(1:2:9999), "Random walk in [-1,1]")
# x^4 only measured on even configurations
x4 = uwreal(x[2:2:10000].^4, 1001, collect(2:2:10000), 10000)
x4 = uwreal(x[2:2:10000].^4, 1001, collect(2:2:10000), "Random walk in [-1,1]")
```
the variables `x2` and `x4` are normal `uwreal` variables, and one can work with them as with any other variable. Note however that the automatic determination of the summation window can fail in these cases, because the autocorrelation function with missing measurements can be very complicated. For example
```@repl gaps
......@@ -153,9 +152,9 @@ uwerr(rat); # Use automatic window: might fail!
```
In this case `uwerr` complains because with the automatically chosen window the variance for ensemble with ID 1001 is negative. Let's have a look at the normalized autocorrelation function
```@repl gaps
iw = window(rat, 1001)
r = rho(rat, 1001);
dr = drho(rat, 1001);
iw = window(rat, "Random walk in [-1,1]")
r = rho(rat, "Random walk in [-1,1]");
dr = drho(rat, "Random walk in [-1,1]");
plot(collect(1:100),
r[1:100],
yerr = dr[1:100],
......@@ -166,7 +165,7 @@ savefig("rat_cf.png") # hide
The normalized autocorrelation function is oscillating, and the automatically chosen window is 2, producing a negative variance! We better fix the window to 50 for this case
```@repl gaps
wpm = Dict{Int64, Vector{Float64}}()
wpm[1001] = [50.0, -1.0, -1.0, -1.0]
wpm["Random walk in [-1,1]"] = [50.0, -1.0, -1.0, -1.0]
uwerr(rat, wpm) # repeat error analysis with our choice (window=50)
println("Ratio: ", rat)
```
......@@ -175,9 +174,9 @@ Note however that it is very difficult to anticipate which observables will have
```@repl gaps
prod = x2*x4
uwerr(prod)
iw = window(prod, 1001)
r = rho(prod, 1001);
dr = drho(prod, 1001);
iw = window(prod, "Random walk in [-1,1]")
r = rho(prod, "Random walk in [-1,1]");
dr = drho(prod, "Random walk in [-1,1]");
plot(collect(1:2*iw),
r[1:2*iw],
yerr = dr[1:2*iw],
......@@ -197,11 +196,13 @@ All observables are stored in a single `BDIO` record of type `BDIO_BIN_GENERIC`.
- ndata (`Vector{Int32}(neid)`): The length of each ensemble.
- nrep (`Vector{Int32}(neid)`): The number of replica of each enemble.
- vrep (`Vector{Int32}`): The replica vector for each ensemble.
- ids (`Vector{Int32}(neid)`): The ensemble `ID`'s.
- ids (`Vector{Int32}(neid)`): The ensemble numeric `ID`'s.
- nt (`Vector{Int32}(neid)`): Half the largest replica of each ensemble.
- zero (`Vector{Float64}(neid)`): just `neid` zeros.
- four (`Vector{Float64}(neid)`): just `neid` fours.
- delta (`Vector{Float64}`): The fluctuations for each ensemble.
- name (NULL terminated `String`): A description of the observable.
- `ID` tags: A list of `neid` tuples `(Int34, String)` that maps each numeric `ID` to an ensemble tag. All strings are NULL terminated.
!!! alert
Obviously this weird format is what it is for some legacy reasons, but it is strongly encouraged that new implementations respect this standard with all its weirdness.
......
src/ADerrorsCF.jl
View file @ fabaca28
This diff is collapsed.
src/ADerrorsIO.jl
View file @ fabaca28
......@@ -74,9 +74,9 @@ function derror(a::uwreal)
end
"""
taui(a::uwreal, id::Int64)
taui(a::uwreal, mcid)
Returns the value of tauint for the ensemble `id`. It is assumed that `uwerr` has been run on the variable and that `id` contributes to the observable `a`. Otherwise an error message is printed.
Returns the value of tauint for the ensemble `mcid`. It is assumed that `uwerr` has been run on the variable and that `mcid` contributes to the observable `a`. Otherwise an error message is printed. `mcid` can be either an `Int64` (the proper ensemble ID), or a `String` (the ensemble tag).
```@example
using ADerrors # hide
# Generate some correlated data
......@@ -90,9 +90,9 @@ for i in 2:1000
end
end
a = uwreal(x.^2, 666)
a = uwreal(x.^2, "Some simple ensemble")
uwerr(a)
println("Error analysis result: ", a, " (tauint = ", taui(a, 666), ")")
println("Error analysis result: ", a, " (tauint = ", taui(a, "Some simple ensemble"), ")")
```
"""
function taui(a::uwreal, mcid::Int64)
......@@ -105,9 +105,9 @@ function taui(a::uwreal, mcid::Int64)
end
"""
dtaui(a::uwreal, id::Int64)
dtaui(a::uwreal, mcid)
Returns an estimate on the error of tauint for the ensemble `id`. It is assumed that `uwerr` has been run on the variable and that `id` contributes to the observable `a`. Otherwise an error message is printed.
Returns an estimate on the error of tauint for the ensemble `mcid`. It is assumed that `uwerr` has been run on the variable and that `mcid` contributes to the observable `a`. Otherwise an error message is printed.
```@example
using ADerrors # hide
# Generate some correlated data
......@@ -121,10 +121,10 @@ for i in 2:1000
end
end
a = uwreal(x.^2, 666)
a = uwreal(x.^2, "Some simple ensemble")
uwerr(a)
println("Error analysis result: ", a,
" (tauint = ", taui(a, 666), " +/- ", dtaui(a, 666), ")")
" (tauint = ", taui(a, "Some simple ensemble"), " +/- ", dtaui(a, "Some simple ensemble"), ")")
```
"""
function dtaui(a::uwreal, mcid::Int64)
......@@ -137,9 +137,9 @@ function dtaui(a::uwreal, mcid::Int64)
end
"""
window(a::uwreal, id::Int64)
window(a::uwreal, mcid)
Returns the summation window for the ensemble `id`. It is assumed that `uwerr` has been run on the variable and that `id` contributes to the observable `a`. Otherwise an error message is printed.
Returns the summation window for the ensemble `mcid`. It is assumed that `uwerr` has been run on the variable and that `mcid` contributes to the observable `a`. Otherwise an error message is printed.
```@example
using ADerrors # hide
# Generate some correlated data
......@@ -153,10 +153,10 @@ for i in 2:1000
end
end
a = uwreal(x.^2, 666)
a = uwreal(x.^2, "Some simple ensemble")
uwerr(a)
println("Error analysis result: ", a,
" (window = ", window(a, 666), ")")
" (window = ", window(a, "Some simple ensemble"), ")")
```
"""
function window(a::uwreal, mcid::Int64)
......@@ -169,9 +169,9 @@ function window(a::uwreal, mcid::Int64)
end
"""
rho(a::uwreal, id::Int64)
rho(a::uwreal, mcid)
Returns the normalized autocorrelation function of `a` for the ensemble `id`. It is assumed that `uwerr` has been run on the variable and that `id` contributes to the observable `a`. Otherwise an error message is printed.
Returns the normalized autocorrelation function of `a` for the ensemble `mcid`. It is assumed that `uwerr` has been run on the variable and that `mcid` contributes to the observable `a`. Otherwise an error message is printed.
```@example
using ADerrors # hide
# Generate some correlated data
......@@ -185,9 +185,9 @@ for i in 2:1000
end
end
a = uwreal(x.^2, 666)
a = uwreal(x.^2, "Some simple ensemble")
uwerr(a)
v = rho(a, 666)
v = rho(a, "Some simple ensemble")
for i in 1:length(v)
println(i, " ", v[i])
end
......@@ -203,9 +203,9 @@ function rho(a::uwreal, mcid::Int64)
end
"""
rho(a::uwreal, id::Int64)
drho(a::uwreal, mcid)
Returns an estimate of the error on the normalized autocorrelation function of `a` for the ensemble `id`. It is assumed that `uwerr` has been run on the variable and that `id` contributes to the observable `a`. Otherwise an error message is printed.
Returns an estimate of the error on the normalized autocorrelation function of `a` for the ensemble `mcid`. It is assumed that `uwerr` has been run on the variable and that `mcid` contributes to the observable `a`. Otherwise an error message is printed.
```@example
using ADerrors # hide
# Generate some correlated data
......@@ -219,10 +219,10 @@ for i in 2:1000
end
end
a = uwreal(x.^2, 666)
a = uwreal(x.^2, "Some simple ensemble")
uwerr(a)
v = rho(a, 666)
dv = drho(a, 666)
v = rho(a, "Some simple ensemble")
dv = drho(a, "Some simple ensemble")
for i in 1:length(v)
println(i, " ", v[i], " +/- ", dv[i])
end
......@@ -278,15 +278,33 @@ function read_bdio(fb, ws::wspace, mapids::Dict{Int64, Int64})
dfl = zeros(Float64, nds[i])
BDIO.BDIO_read(fb, dfl)
id = convert(Int64, ids[i])
add_DB(dfl, get(mapids, id, id), convert(Vector{Int64}, ivrep[is:ie]), ws)
add_DB(dfl, get(mapids, id, id), convert(Vector{Int64}, ivrep[is:ie]), ws, false)
is = ie + 1
end
if BDIO.BDIO_eor(fb)
for i in 1:nid
id = convert(Int64, ids[i])
add_maps(id, ws)
get_name_from_id(id, ws)
end
else
name = BDIO.BDIO_read_str(fb)
for i in 1:nid
BDIO.BDIO_read(fb, ifoo)
str = BDIO.BDIO_read_str(fb)
id = get_id_from_name(str, ws)
add_maps(id, ws)
end
end
return uwreal(dfoo[1], p, d)
end
function write_bdio(p::uwreal, fb, iu::Int, ws::wspace)
function write_bdio(p::uwreal, fb, iu::Int, ws::wspace; name="NO NAME")
BDIO.BDIO_start_record!(fb, BDIO.BDIO_BIN_GENERIC, convert(Int32, iu), true)
nid = convert(Int32, unique_ids!(p, ws))
......@@ -329,6 +347,12 @@ function write_bdio(p::uwreal, fb, iu::Int, ws::wspace)
BDIO.BDIO_write!(fb, dt, true)
end
BDIO.BDIO_write!(fb, name*"\0")
for i in 1:nid
BDIO.BDIO_write!(fb, [convert(Int32, p.ids[i])])
BDIO.BDIO_write!(fb, get_name_from_id(p.ids[i], ws)*"\0")
end
BDIO.BDIO_write_hash!(fb)
return true
......@@ -341,24 +365,21 @@ Write out a detailed information on the error of `a`.
## Arguments
Optionally one can pass as a keyword argument (`io`) the `IO` stream to write to and a dictionary (`names`) that translates ensemble `ID`s into more human friendly `Strings`
Optionally one can pass as a keyword argument (`io`) the `IO` stream to write to.
## Example
```@example
using ADerrors # hide
a = uwreal(rand(2000), 120)
b = uwreal([1.2, 0.023], 121)
c = uwreal([5.2, 0.03], 122)
a = uwreal(rand(2000), "Ensemble A12")
b = uwreal([1.2, 0.023], "Ensemble XYZ")
c = uwreal([5.2, 0.03], "Ensemble RRR")
d = a + b - c
uwerr(d)
bnm = Dict{Int64, String}()
bnm[120] = "Very important ensemble"
bnm[122] = "H12B K87"
details(d, bnm)
details(d)
```
"""
function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String} = Dict{Int64, String}())
function details(a::uwreal, ws::wspace, io::IO=stdout)
if (length(a.prop) == 0)
print(a.mean)
......@@ -390,7 +411,7 @@ function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String
for i in 1:length(a.cfd)
idx = ws.map_ids[a.ids[ip[i]]]
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_name_from_id(a.ids[ip[i]], ws), ntrunc)
if (ws.fluc[idx].nd > 1)
Printf.@printf(io, " # %45s %6.2f %s\n",
sid, 100.0 .* a.cfd[ip[i]].var ./ a.err^2, sndt)
......@@ -404,7 +425,7 @@ function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String
end
end
details(a::uwreal; io::IO=stdout, names::Dict{Int64, String} = Dict{Int64, String}()) = details(a, wsg, io, names)
details(a::uwreal; io::IO=stdout) = details(a, wsg, io)
"""
read_uwreal(fb[, map_ids::Dict{Int64, Int64}])
......
src/ADerrorsTools.jl
View file @ fabaca28
......@@ -10,9 +10,10 @@
###
"""
cobs(avgs::Vector{Float64}, Mcov::Array{Float64, 2}, str::String)
cobs(avgs::Vector{Float64}, Mcov::Array{Float64, 2}, ids::Vector{Int64})
Returns a vector of `uwreal` such that their mean values are `avgs[:]` and their covariance is `Mcov[:,:]`. In order to construct these observables `n=length(avgs)` ensemble ID are used. These have to be specified in the vector `ids[:]`.
Returns a vector of `uwreal` such that their mean values are `avgs[:]` and their covariance is `Mcov[:,:]`. In order to construct these observables `n=length(avgs)` ensemble ID are used. These are generated either from the string `str` or have to be specified in the vector `ids[:]`.
```@example
using ADerrors # hide
# Put some average values and covariance
......@@ -23,7 +24,7 @@ Mcov = [0.478071 -0.176116 0.0135305
# Produce observables with ensemble ID
# [1, 2001, 32]. Do error analysis
p = cobs(avg, Mcov, [1, 2001, 32])
p = cobs(avg, Mcov, "GF beta function parameters")
uwerr.(p)
# Check central values are ok
......@@ -50,6 +51,19 @@ function cobs(avgs::Vector{Float64}, cov::Array{Float64, 2}, ids::Vector{Int64})
return p
end
function cobs(avgs::Vector{Float64}, cov::Array{Float64, 2}, str::String)
n = length(avgs)
ids = Vector{Int64}(undef, n)
for i in 1:n
tstr = Printf.@sprintf "%s %8.8d" str i
ids[i] = get_id_from_name(tstr)
end
return cobs(avgs, cov, ids)
end
function addobs(a::Vector{uwreal}, der::Vector{Float64}, mvl::Float64)
n = 0
......
src/ADerrorsTypes.jl
View file @ fabaca28
......@@ -73,6 +73,11 @@ mutable struct wspace
map_nob::Array{Int64, 1} # The id that corresponds to each ob
map_ids::Dict{Int64, Int64} # For each id, a fluc index
id2str::Dict{Int64, String} # Ensemble id for each String
str2id::Dict{String, Int64} # Ensemble String for each id
newid::Int64
end
......
src/ADerrorsUtils.jl
View file @ fabaca28
......@@ -17,9 +17,9 @@ Returns `x` such that `fnew(x, data) = 0` (i.e. a root of the function). The err
using ADerrors # hide
# First define some arbitrary data
data = Vector{uwreal}(undef, 3)
data[1] = uwreal([1.0, 0.2], 120)
data[2] = uwreal([1.2, 0.023], 121)
data[3] = uwreal(rand(1000), 122)
data[1] = uwreal([1.0, 0.2], "Var A")
data[2] = uwreal([1.2, 0.023], "Var B")
data[3] = uwreal(rand(1000), "White noise ensemble")
# Now define a function
f(x, p) = x + p[1]*x + cos(p[2]*x+p[3])
......@@ -33,7 +33,8 @@ println("Root: ", x)
# Check
z = f(x, data)
uwerr(z)
println("Better be zero (with zero error): ", z)
print("Better be zero (with zero error): ")
details(z)
```
"""
function root_error(fnew::Function, x::Float64,
......@@ -59,7 +60,7 @@ function root_error(fnew::Function, x::Float64,
return addobs(data, -gw ./ grad[1], xv[1])
end
function chiexp(hess::Array{Float64, 2}, data::Vector{uwreal}, W::Vector{Float64})
function chiexp(hess::Array{Float64, 2}, data::Vector{uwreal}, W::Vector{Float64}, wpm::Dict{Int64,Vector{Float64}})
m = length(data)
n = size(hess, 1) - m
......@@ -77,10 +78,10 @@ function chiexp(hess::Array{Float64, 2}, data::Vector{uwreal}, W::Vector{Float64
Px[i,i] = W[i] + Px[i,i]
end
return trcov(Px, data)
return trcov(Px, data, wpm)
end
function chiexp(hess::Array{Float64, 2}, data::Vector{uwreal}, W::Array{Float64, 2})
function chiexp(hess::Array{Float64, 2}, data::Vector{uwreal}, W::Array{Float64, 2}, wpm::Dict{Int64,Vector{Float64}})
m = length(data)
n = size(hess, 1) - m
......@@ -95,7 +96,7 @@ function chiexp(hess::Array{Float64, 2}, data::Vector{uwreal}, W::Array{Float64,
hi = LinearAlgebra.pinv(maux)
Px = W - hm' * hi * hm
return trcov(Px, data)
return trcov(Px, data, wpm)
end
@doc raw"""
......@@ -153,8 +154,9 @@ println("chi^2 / chi_exp^2: ", chisq(xp, value.(dt)), " / ", chiexp(chisq, xp, d
"""
function chiexp(chisq::Function,
xp::Vector{Float64},
data::Vector{uwreal};
W = Vector{Float64}())
data::Vector{uwreal},
wpm::Dict{Int64,Vector{Float64}};
W::Vector{Float64} = Vector{Float64}())
n = length(xp) # Number of fit parameters
m = length(data) # Number of data
......@@ -182,7 +184,7 @@ function chiexp(chisq::Function,
Ww = zeros(Float64, m)
for i in 1:m
if (data[i].err == 0.0)
uwerr(data[i])
uwerr(data[i], wpm)
if (data[i].err == 0.0)
error("Zero error in fit data")
end
......@@ -193,15 +195,28 @@ function chiexp(chisq::Function,
Ww = W
end
cse = chiexp(hess, data, Ww)
cse = chiexp(hess, data, Ww, wpm)
end
return cse
end
chiexp(chisq::Function,
xp::Vector{Float64},
data::Vector{uwreal};
W::Vector{Float64} = Vector{Float64}()) =
chiexp(chisq, xp, data, Dict{Int64,Vector{Float64}}(), W)
chiexp(chisq::Function,
xp::Vector{Float64},
data::Vector{uwreal},
wpm::Dict{String,Vector{Float64}};
W::Vector{Float64} = Vector{Float64}()) =
chiexp(chisq, xp, data, dict_names_to_id(wpm), W)
@doc raw"""
fit_error(chisq::Function, xp::Vector{Float64}, data::Vector{uwreal};
fit_error(chisq::Function, xp::Vector{Float64}, data::Vector{uwreal}[, wpm];
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)``.
......@@ -215,6 +230,7 @@ Given a ``\chi^2(p, d)``, function of the fit parameters `p[:]` and the data `d[
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`.
- `wpm`: `Dict{Int64,Vector{Float64}}` or `Dict{String,Vector{Float64}}`. The criteria to determine the summation window. See the documentation on `uwerr` function for more details.
- `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)``.
......@@ -239,7 +255,7 @@ 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])
dt = cobs(vs[:,1], sig, "Data points")
# Define the chi^2
chisq(p, d) = sum( (d .- p[1]) .^ 2 ./ dx .^2 )
......@@ -254,14 +270,18 @@ xp = [sum(value.(dt) ./ dx)/sum(1.0 ./ dx)]
(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)
println(" *** FIT RESULTS ***")
print("Fit parameter: ")
details.(fitp)
println("chi^2 / chi_exp^2: ", chisq(xp, value.(dt)), " / ", csqexp, " (dof: ", npt-1, ")")
```
"""
function fit_error(chisq::Function,
xp::Vector{Float64},
data::Vector{uwreal};
W = Vector{Float64}(), chi_exp = true)
data::Vector{uwreal},
wpm::Dict{Int64,Vector{Float64}};
W::Vector{Float64} = Vector{Float64}(),
chi_exp::Bool = true)
n = length(xp) # Number of fit parameters
m = length(data) # Number of data
......@@ -301,7 +321,7 @@ function fit_error(chisq::Function,
Ww = zeros(Float64, m)
for i in 1:m
if (data[i].err == 0.0)
uwerr(data[i])
uwerr(data[i], wpm)
if (data[i].err == 0.0)
error("Zero error in fit data")
end
......@@ -312,12 +332,87 @@ function fit_error(chisq::Function,
Ww = W
end
cse = chiexp(hess, data, Ww)
cse = chiexp(hess, data, Ww, wpm)
end
return param, cse
end
function fit_error(chisq::Function,
xp::Vector{Float64},
data::Vector{uwreal},
wpm::Dict{Int64,Vector{Float64}},
W::Array{Float64, 2};
chi_exp::Bool = true)
n = length(xp) # Number of fit parameters
m = length(data) # Number of data
xav = Vector{Float64}(undef, n+m)
for i in 1:n
xav[i] = xp[i]
end
for i in n+1:n+m
xav[i] = data[i-n].mean
end
ccsq(x::Vector) = chisq(x[1:n], x[n+1:n+m])
if (n+m < 4)
cfg = ForwardDiff.HessianConfig(ccsq, xav, Chunk{1}());
else
cfg = ForwardDiff.HessianConfig(ccsq, xav, Chunk{4}());
end
hess = Array{Float64}(undef, n+m, n+m)
ForwardDiff.hessian!(hess, ccsq, xav, cfg)
hm = view(hess, 1:n, 1:n)
sm = view(hess, 1:n, n+1:n+m)
hinv = LinearAlgebra.pinv(hm)
grad = - hinv * sm
param = addobs(data, grad, xp)
if (!chi_exp)
return param
end
cse = 0.0
if (m-n > 0)
cse = chiexp(hess, data, W, wpm)
end
return param, cse
end
fit_error(chisq::Function,
xp::Vector{Float64},
data::Vector{uwreal};
W::Vector{Float64} = Vector{Float64}(),
chi_exp::Bool = true) =
fit_error(chisq, xp, data, Dict{Int64,Vector{Float64}}(), W=W, chi_exp=chi_exp)
fit_error(chisq::Function,
xp::Vector{Float64},
data::Vector{uwreal},
W::Array{Float64,2};
chi_exp::Bool = true) =
fit_error(chisq, xp, data, Dict{Int64,Vector{Float64}}(), W=W, chi_exp=chi_exp)
fit_error(chisq::Function,
xp::Vector{Float64},
data::Vector{uwreal},
wpm::Dict{String,Vector{Float64}};
W::Vector{Float64} = Vector{Float64}(),
chi_exp::Bool = true) =
fit_error(chisq, xp, data, dict_names_to_id(wpm), W=W, chi_exp=chi_exp)
fit_error(chisq::Function,
xp::Vector{Float64},
data::Vector{uwreal},
wpm::Dict{String,Vector{Float64}},
W::Array{Float64,2};
chi_exp::Bool = true) =
fit_error(chisq, xp, data, dict_names_to_id(wpm), W=W, chi_exp=chi_exp)
@doc raw"""
int_error(fint::Function, a, b, p::Vector{uwreal})
......@@ -334,11 +429,11 @@ using ADerrors
# 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]
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])
p = cobs(avg, Mcov, "Beta function fit parameters")
g1s = uwreal([2.6723, 0.0064], 4)
g2s = 11.31
......
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