Fully functional version, but unfortunately slow

Still some work required to avoid allocations in functions:
- uwerr
- basic Math operations
- addobs

Project.toml

@@ -4,7 +4,18 @@ authors = ["Alberto Ramos <alberto.ramos@ific.uv.es>"]
 version = "0.1.0"
 
 [deps]
+BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Nettle = "49dea1ee-f6fa-5aa6-9a11-8816cee7d4b9"
+PGFPlotsX = "8314cec4-20b6-5062-9cdb-752b83310925"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"

src/ADerrors.jl

@@ -11,7 +11,7 @@
 
 module ADerrors
 
-import ForwardDiff, Statistics, FFTW#, BDIO
+import ForwardDiff, Statistics, FFTW, LinearAlgebra, QuadGK, BDIO, Printf
 
 # Include data types
 include("ADerrorsTypes.jl")
@@ -24,7 +24,18 @@ include("ADerrorsMath.jl")
 
 # I/O
 include("ADerrorsIO.jl")
-export err, value, derror, taui, dtaui, window
+
+# Basic tools
+include("ADerrorsTools.jl")
+
+# Root, fit, integral error propagation
+include("ADerrorsUtils.jl")
+
+export err, value, derror, taui, dtaui, window, rho, drho, details
 export uwreal, uwerr
+export cov, trcov, trcorr, neid
+export read_uwreal, write_uwreal
+export addobs, cobs
+export root_error, chiexp, fit_error, int_error
 
 end # module

src/ADerrorsIO.jl

 
+function find_mcid(a::uwreal, mcid::Int64)
+
+    if (length(a.cfd) == 0)
+        return 0
+    else
+        for i in 1:length(a.ids)
+            if (a.ids[i] == mcid)
+                return i
+            end
+        end
+    end
+    
+    return 0
+end
+
 err(a::uwreal)              = a.err
 value(a::uwreal)            = a.mean
 derror(a::uwreal)           = a.derr
-taui(a::uwreal, i::Int64)   = a.cfd[i].taui
-dtaui(a::uwreal, i::Int64)  = a.cfd[i].taui
-window(a::uwreal, i::Int64) = a.cfd[i].iw
+
+function taui(a::uwreal,   mcid::Int64)
+    idx = find_mcid(a, mcid)
+    if (idx == 0)
+        return 0.5
+    else
+        return a.cfd[idx].taui
+    end
+end
+
+function dtaui(a::uwreal,  mcid::Int64)
+    idx = find_mcid(a, mcid)
+    if (idx == 0)
+        return 0.0
+    else
+        return a.cfd[idx].dtaui
+    end
+end
+
+function window(a::uwreal, mcid::Int64)
+    idx = find_mcid(a, mcid)
+    if (idx == 0)
+        return 0
+    else
+        return a.cfd[idx].iw
+    end
+end
+
+function rho(a::uwreal, mcid::Int64)
+    idx = find_mcid(a, mcid)
+    if (idx == 0)
+        return [1.0]
+    else
+        return a.cfd[idx].gamm ./ a.cfd[idx].gamm[1]
+    end
+end
+
+function drho(a::uwreal, mcid::Int64)
+    idx = find_mcid(a, mcid)
+    if (idx == 0)
+        return [0.0]
+    else
+        return a.cfd[idx].drho
+    end
+end
+
+function read_bdio(fb, ws::wspace)
+
+    dfoo = zeros(Float64, 1)
+    ifoo = zeros(Int32, 1)
+    BDIO.BDIO_read(fb, dfoo)
+    BDIO.BDIO_read(fb, ifoo)
+
+    nid::Int32 = ifoo[1]
+
+    p = [false for i in 1:ws.nob+nid]
+    d = zeros(Float64, ws.nob+nid)
+    p[ws.nob+1:end] .= true
+    d[ws.nob+1:end] .= 1.0
+
+    nds  = zeros(Int32, nid)
+    nrep = zeros(Int32, nid)
+    ids  = zeros(Int32, nid)
+    itmp = zeros(Int32, nid)
+
+    BDIO.BDIO_read(fb, nds)
+    BDIO.BDIO_read(fb, nrep)
+    ivrep = zeros(Int32, sum(nrep))
+    BDIO.BDIO_read(fb, ivrep)
+    BDIO.BDIO_read(fb, ids)
+
+    dfl = zeros(Float64, nid)
+    BDIO.BDIO_read(fb, itmp)
+    for i in 1:2
+        BDIO.BDIO_read(fb, dfl)
+    end
+
+    is = 1
+    for i in 1:nid
+        ie = is + nrep[i] - 1
+        if (sum(ivrep[is:ie]) != nds[i])
+            throw("Replica sum does not match number of measurements")
+        end
+
+        dfl = zeros(Float64, nds[i])
+        BDIO.BDIO_read(fb, dfl)
+        add_DB(dfl, convert(Int64, ids[i]), convert(Vector{Int64}, ivrep[is:ie]), ws)
+        
+        is = ie + 1
+    end
+
+    return uwreal(dfoo[1], p, d)
+end
+
+function write_bdio(p::uwreal, fb, iu::Int, ws::wspace)
+
+    BDIO.BDIO_start_record!(fb, BDIO.BDIO_BIN_GENERIC, convert(Int32, iu), true)
+    nid = convert(Int32, unique_ids!(p, ws))
+
+    BDIO.BDIO_write!(fb, [p.mean], true)
+    BDIO.BDIO_write!(fb, [nid], true)
+    for j in 1:nid
+        nd = convert(Int32, ws.fluc[ws.map_ids[p.ids[j]]].nd)
+        BDIO.BDIO_write!(fb, [nd], true)
+    end
+
+    for j in 1:nid
+        nr = convert(Int32, length(ws.fluc[ws.map_ids[p.ids[j]]].ivrep))
+        BDIO.BDIO_write!(fb, [nr], true)
+    end
+
+    for j in 1:nid
+        BDIO.BDIO_write!(fb, convert(Vector{Int32}, ws.fluc[ws.map_ids[p.ids[j]]].ivrep), true)
+    end
+
+    BDIO.BDIO_write!(fb, convert(Vector{Int32}, p.ids), true)
+    
+    for j in 1:nid
+        for i in 1:length(p.prop)
+            if (p.prop[i] && (ws.map_nob[i] == p.ids[j]))
+                (nt,ip) = findmax(ws.fluc[i].ivrep)
+                nt = convert(Int32, div(nt, 2*ws.fluc[i].ibn))
+                BDIO.BDIO_write!(fb, [nt], true)
+                continue
+            end
+        end
+    end
+
+    BDIO.BDIO_write!(fb, zeros(Float64, nid), true)
+    BDIO.BDIO_write!(fb, [DEFAULT_STAU for n in 1:nid], true)
+
+    for j in 1:nid
+        nd = ws.fluc[ws.map_ids[p.ids[j]]].nd
+        dt = zeros(Float64, nd)
+        for i in 1:length(p.prop)
+            if (p.prop[i] && (ws.map_nob[i] == p.ids[j]))
+                dt .= dt .+ p.der[i] .* ws.fluc[i].delta
+            end
+        end
+        BDIO.BDIO_write!(fb, dt, true)
+    end
+
+    BDIO.BDIO_write_hash!(fb)
+
+    return true
+end
+
+"""
+    details(a::uwreal; io::IO=stdout, names::Dict{Int64, String} = Dict{Int64, String}())
+
+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`
+
+## Example
+```@example
+using ADerrors # hide
+a = uwreal(rand(2000),   120)
+b = uwreal([1.2, 0.023], 121)
+c = uwreal([5.2, 0.03],  122)
+d = a + b - c
+uwerr(d)
+
+bnm = Dict{Int64, String}()
+bnm[120] = "Very important ensemble"
+bnm[122] = "H12B K87"
+details(d, bnm)
+```
+"""
+function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String} = Dict{Int64, String}())
+    
+    if (length(a.prop) == 0)
+        print(a.mean)
+        return
+    end
+
+    if (length(a.cfd) > 0) 
+        println(a.mean, " +/- ", a.err)
+        println(" ## Number of error sources: ", length(a.ids))
+
+        n = 0
+        for i in 1:length(a.cfd)
+            if (length(a.cfd[i].gamm) > 0)
+                n = n + 1
+            end
+        end
+        println(" ## Number of MC ids       : ", n)
+        println(" ## Contribution to error  :               Ensemble  [%]     [MC length]")
+
+        
+        truncate_ascii(s::String,n::Int) = s[1:min(sizeof(s),n)]
+        ntrunc = 45
+        v = zeros(length(a.cfd))
+        for i in eachindex(v)
+            v[i] = a.cfd[i].var
+        end
+        ip = sortperm(v, rev=true)
+        for i in 1:length(a.cfd)
+            sid = truncate_ascii(get(names, a.ids[ip[i]], string(a.ids[ip[i]])), ntrunc)
+            if (length(a.cfd[ip[i]].gamm) > 0)
+                idx  = ws.map_ids[a.ids[i]]
+                nd   = ws.fluc[idx].nd
+                Printf.@printf("  #  %45s %6.2f   %10d\n",
+                        sid, 100.0 .* a.cfd[ip[i]].var ./ a.err^2, nd)
+            else
+                Printf.@printf("  #  %45s %6.2f            -\n",
+                        sid, 100.0 .* a.cfd[ip[i]].var ./ a.err^2)
+            end
+        end
+    else
+        print(a.mean, " (Error not available... maybe run uwerr)")
+    end
+end
+
+details(a::uwreal; io::IO=stdout, names::Dict{Int64, String} = Dict{Int64, String}()) = details(a, wsg, io, names)
+
+read_uwreal(fb)  = read_bdio(fb, ADerrors.wsg)
+write_uwreal(p::uwreal, fb, iu::Int) = write_bdio(p::uwreal, fb, iu::Int, ADerrors.wsg)
 

src/ADerrorsTools.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    ADerrorsTools.jl
+### created: Fri Jun 26 21:48:35 2020
+###                               
+
+function cobs(avgs::Vector{Float64}, cov::Array{Float64, 2}, ids::Vector{Int64})
+
+    ch = LinearAlgebra.cholesky(cov)
+
+    n = length(avgs)
+    p = Vector{uwreal}()
+    for j in 1:n
+        new = uwreal([avgs[j], ch.L[j, 1]], ids[1])
+        for i in 2:n
+            new = new + uwreal([0.0, ch.L[j,i]], ids[i])
+        end
+        push!(p, new)
+    end
+    return p
+end
+
+function addobs(a::Vector{uwreal}, der::Vector{Float64}, mvl::Float64)
+
+    n = 0
+    for i in 1:length(a)
+        n = max(n,length(a[i].der))
+    end
+    p = [false for i in 1:n]
+    d = zeros(Float64, n)
+
+    for k in 1:n
+        for j in 1:length(a)
+            if (k <= length(a[j].prop))
+                p[k] = p[k] || a[j].prop[k]
+                if (a[j].prop[k])
+                    d[k] = d[k] + der[j]*a[j].der[k]
+                end
+            end
+        end
+    end
+    
+    return uwreal(mvl, p, d)
+end
+
+function addobs(a::Vector{uwreal}, der::Array{Float64, 2}, mvl::Vector{Float64})
+
+    n = 0
+    for i in 1:length(a)
+        n = max(n,length(a[i].der))
+    end
+
+    x = Vector{uwreal}()
+    for i in 1:size(der, 1)
+        p = [false for i in 1:n]
+        d = zeros(Float64, n)
+        for k in 1:n
+            for j in 1:length(a)
+                if (k <= length(a[j].prop))
+                    p[k] = p[k] || a[j].prop[k]
+                    if (a[j].prop[k])
+                        d[k] = d[k] + der[i,j]*a[j].der[k]
+                    end
+                end
+            end
+        end
+        push!(x, uwreal(mvl[i], p, d))
+    end
+    
+    return x
+end

src/ADerrorsTypes.jl

@@ -32,6 +32,22 @@ mutable struct uwreal
     ids::Array{Int64, 1}
 
     cfd::Vector{cfdata}
+
+    uwreal(a::Float64, b::Float64, c::Float64,
+           d::Array{Bool, 1}, e::Array{Float64, 1},
+           f::Array{Int64, 1}, g::Vector{cfdata}) = new(a, b, c, d, e, f, g) 
+    function uwreal(n::Int64)
+        x = new()
+        x.mean = 0.0
+        x.err  = 0.0
+        x.derr = 0.0
+        x.prop = Vector{Bool}(undef, n)
+        x.der  = Vector{Float64}(undef, n)
+
+        x.prop .= false
+        x.der  .= 0.0
+        return x
+    end
 end
 
 mutable struct fbd
@@ -52,6 +68,10 @@ mutable struct wspace
 end
 
 
+Base.convert(::Type{uwreal}, x::Float64) = uwreal(x, 0.0, 0.0, 
+                                            Vector{Bool}(), Vector{Float64}(), 
+                                            Vector{Int64}(), Vector{cfdata}())
+
 uwreal(v::Float64, n::Int) = uwreal(v, 0.0, 0.0,
                                     [false for i in 1:n], [0.0 for i in 1:n],
                                     Vector{Int64}(), Vector{cfdata}())
@@ -60,6 +80,12 @@ uwreal(v::Float64, prop::Vector{Bool}, der::Vector{Float64}) = uwreal(v, 0.0, 0.
                                     Vector{Int64}(), Vector{cfdata}())
 
 function Base.show(io::IO, a::uwreal)
+    
+    if (length(a.prop) == 0)
+        print(a.mean)
+        return
+    end
+    
     if (length(a.cfd) > 0) 
         print(a.mean, " +/- ", a.err)
     else

src/ADerrorsUtils.jl

+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    ADerrorsUtils.jl
+### created: Fri Jun 26 19:30:27 2020
+###                               
+
+function root_error(fnew::Function, x::Float64,
+                    data::Vector{uwreal}) 
+
+    function fvec(x::Vector)
+        return fnew(x[1], x[2:end])
+    end
+    xv = zeros(Float64, length(data)+1)
+    xv[1] = x
+    for i in 1:length(data)
+        xv[i+1] = data[i].mean
+    end
+    grad = ForwardDiff.gradient(fvec, xv)
+
+    return addobs(data, -grad[2:end]./grad[1], x)
+end
+
+function chiexp(hess::Array{Float64, 2}, data::Vector{uwreal}, W::Vector{Float64})
+
+    m = length(data)
+    n = size(hess, 1) - m
+
+    hm = view(hess, 1:n, n+1:n+m)
+    sm = hm * LinearAlgebra.Diagonal(1.0 ./ sqrt.(W))
+
+    Px = LinearAlgebra.Symmetric(LinearAlgebra.Diagonal(W)) -
+        LinearAlgebra.Symmetric(hm' *
+                                LinearAlgebra.inv(LinearAlgebra.Symmetric(sm * sm')) *
+                                hm)
+
+    return trcov(Px, data)
+end
+    
+function chiexp(hess::Array{Float64, 2}, data::Vector{uwreal}, W::Array{Float64, 2})
+
+    m = length(data)
+    n = size(hess, 1) - m
+    
+    Lm = LinearAlgebra.cholesky(LinearAlgebra.Symmetric(W))
+    Li = LinearAlgebra.inv(Lm.L)
+    
+    hm = view(hess, 1:n, n+1:n+m)
+    sm = hm * Li'
+    
+    Px = LinearAlgebra.Symmetric(W) -
+        LinearAlgebra.Symmetric(hm' *
+                                LinearAlgebra.inv(LinearAlgebra.Symmetric(sm * sm')) *
+                                hm)
+    
+    return trcov(Px, data)
+end
+
+function chiexp(chisq::Function,
+                xp::Vector{Float64}, 
+                data::Vector{uwreal};
+                W = Vector{Float64}())
+
+    n = length(xp)   # Number of fit parameters
+    m = length(data) # Number of data
+    ccsq(x::Vector) = chisq(x[1:n], x[n+1:end])
+
+    xav = zeros(Float64, 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
+    hess = ForwardDiff.hessian(ccsq, xav)
+
+    cse = 0.0
+    if (m-n > 0)
+        if (length(W) == 0)
+            Ww = zeros(Float64, m)
+            for i in 1:m
+                if (data[i].err == 0.0)
+                    uwerr(data[i])
+                    if (data[i].err == 0.0)
+                        error("Zero error in fit data")
+                    end
+                end
+                Ww[i] = 1.0 / data[i].err^2
+            end
+        else
+            Ww = W
+        end
+
+        cse = chiexp(hess, data, Ww)        
+    end
+
+    return cse
+end
+
+function fit_error(chisq::Function,
+                   xp::Vector{Float64}, 
+                   data::Vector{uwreal};
+                   W = Vector{Float64}(), chi_exp = true)
+
+    n = length(xp)   # Number of fit parameters
+    m = length(data) # Number of data
+    ccsq(x::Vector) = chisq(x[1:n], x[n+1:end])
+
+    xav = zeros(Float64, 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
+    hess = ForwardDiff.hessian(ccsq, xav)
+
+    hinv = LinearAlgebra.pinv(hess[1:n,1:n])
+    grad = - hinv[1:n,1:n] * hess[1:n,n+1:n+m]
+    
+    param = addobs(data, grad, xp)
+
+    if (!chi_exp)
+        return param
+    end
+    
+    cse = 0.0
+    if (m-n > 0)
+        if (length(W) == 0)
+            Ww = zeros(Float64, m)
+            for i in 1:m
+                if (data[i].err == 0.0)
+                    uwerr(data[i])
+                    if (data[i].err == 0.0)
+                        error("Zero error in fit data")
+                    end
+                end
+                Ww[i] = 1.0 / data[i].err^2
+            end
+        else
+            Ww = W
+        end
+
+        cse = chiexp(hess, data, Ww)        
+    end
+
+    return param, cse
+end
+
+function int_error(fint::Function,
+                   a,
+                   b,
+                   param::Vector{uwreal})
+
+    # First basic integral
+    xp = zeros(Float64, length(param))
+    for i in 1:length(xp)
+        xp[i] = param[i].mean
+    end
+    f(x) = fint(x, xp)
+    (val, foo) = QuadGK.quadgk(f, a, b)
+
+    # Functions that return the derivatives
+    # with respect to the parameters as functions
+    function fp(n::Int64)
+        function fd(x)
+            fax(pa) = fint(x,[xp[1:n-1];[pa];xp[n+1:end]])
+            return ForwardDiff.derivative(fax, xp[n])
+        end
+        return fd
+    end
+
+    grad = zeros(Float64, length(param))
+    for i in 1:length(param)
+        (grad[i], foo) = QuadGK.quadgk(fp(i), a, b)
+    end
+
+    return addobs(param, grad, val)
+end
+
+function int_error(fint::Function,
+                   a::uwreal,
+                   b::Float64,
+                   param::Vector{uwreal})
+
+    # First basic integral
+    xp = zeros(Float64, length(param))
+    for i in 1:length(xp)
+        xp[i] = param[i].mean
+    end
+    f(x) = fint(x, xp)
+    (val, foo) = QuadGK.quadgk(f, a.mean, b)
+
+    # Functions that return the derivatives
+    # with respect to the parameters as functions
+    function fp(n::Int64)
+        function fd(x)
+            fax(pa) = fint(x,[xp[1:n-1];[pa];xp[n+1:end]])
+            return ForwardDiff.derivative(fax, xp[n])
+        end
+        return fd
+    end
+
+    grad = zeros(Float64, length(param)+1)
+    for i in 1:length(param)
+        (grad[i], foo) = QuadGK.quadgk(fp(i), a.mean, b)
+    end
+    grad[end] = fint(a.mean, xp)
+
+    return addobs([param; a], grad, val)
+end
+
+function int_error(fint::Function,
+                   a::uwreal,
+                   b::uwreal,
+                   param::Vector{uwreal})
+
+    # First basic integral
+    xp = zeros(Float64, length(param))
+    for i in 1:length(xp)
+        xp[i] = param[i].mean
+    end
+    f(x) = fint(x, xp)
+    (val, foo) = QuadGK.quadgk(f, a.mean, b.mean)
+
+    # Functions that return the derivatives
+    # with respect to the parameters as functions
+    function fp(n::Int64)
+        function fd(x)
+            fax(pa) = fint(x,[xp[1:n-1];[pa];xp[n+1:end]])
+            return ForwardDiff.derivative(fax, xp[n])
+        end
+        return fd
+    end
+
+    grad = zeros(Float64, length(param)+2)
+    for i in 1:length(param)
+        (grad[i], foo) = QuadGK.quadgk(fp(i), a.mean, b.mean)
+    end
+    grad[end-1] = fint(a.mean, xp)
+    grad[end]   = fint(b.mean, xp)
+
+    return addobs([param; a; b], grad, val)
+end
+
+function int_error(fint::Function,
+                   a::Float64,
+                   b::uwreal,
+                   param::Vector{uwreal})
+
+    # First basic integral
+    xp = zeros(Float64, length(param))
+    for i in 1:length(xp)
+        xp[i] = param[i].mean
+    end
+    f(x) = fint(x, xp)
+    (val, foo) = QuadGK.quadgk(f, a, b.mean)
+
+    # Functions that return the derivatives
+    # with respect to the parameters as functions
+    function fp(n::Int64)
+        function fd(x)
+            fax(pa) = fint(x,[xp[1:n-1];[pa];xp[n+1:end]])
+            return ForwardDiff.derivative(fax, xp[n])
+        end
+        return fd
+    end
+
+    grad = zeros(Float64, length(param)+1)
+    for i in 1:length(param)
+        (grad[i], foo) = QuadGK.quadgk(fp(i), a, b.mean)
+    end
+    grad[end]   = fint(b.mean, xp)
+
+    return addobs([param; b], grad, val)
+end