###
### "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
###                               

"""
    derivative(a::uwreal, b::uwreal)

Determines the derivate of the observable `a` with respect to observable `b` (i.e. ``\frac{{\rm d} a}{{\rm d} b}). Observable `b` must depend on only a single ensemble (typically `b` is just a primary observable). 
```@example
using ADerrors # hide

# Put some data 
a = uwreal(1.2 .+ randn(2000), "Ensemble test A", ["R0", "Rep 1", "Last rep"], [1000, 550, 450])
b = uwreal(0.8 .+ randn(2034), "Ensemble test B", [1000, 584, 450])

d = sin(a+b)
println("Derivative of d w.r.t a:      ", derivative(d, a))
println("Derivative of d w.r.t log(a): ", derivative(d, log(a)))
```
"""
function derivative(a::uwreal, p::uwreal)

    if count(p.prop .== true) != 1
        error("I do not know how to compute this")
    end
    
    for i in 1:min(length(a.der),length(p.der))
        if (p.prop[i] && a.prop[i])
            return res = a.der[i]/p.der[i]
        end 
    end 
    
    return 0.0


end


"""
    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 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 
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]

# Produce observables with ensemble ID 
# [1, 2001, 32]. Do error analysis
p = cobs(avg, Mcov, "GF beta function parameters")
uwerr.(p)

# Check central values are ok
avg2 = value.(p)
println("Better be zero: ", sum((avg.-avg2).^2))

# Check that the covariance is ok
Mcov2 = cov(p)
println("Better be zero: ", sum((Mcov.-Mcov2).^2))
```
"""
function cobs(avgs::Vector{Float64}, cov::Array{Float64, 2}, ids::Vector{Int64})


    n = length(avgs)
    p = Vector{uwreal}(undef, n)
    try
        ch = LinearAlgebra.cholesky(cov)
        
        for j in 1:n
            p[j] = uwreal([avgs[j], ch.L[j, 1]], ids[1])
            for i in 2:n
                p[j] = p[j] + uwreal([0.0, ch.L[j,i]], ids[i])
            end
        end

        return p
    catch
        usvt = LinearAlgebra.svd(cov)
        SS = similar(usvt.S)
        for i in eachindex(SS)
            if (usvt.S[i] > 1.0e-6)
                SS[i] = usvt.S[i]
            else
                SS[i] = 0.0
            end
        end
        A    = usvt.U * Diagonal(sqrt.(usvt.S))

        for j in 1:n
            p[j] = uwreal([avgs[j], A[j, 1]], ids[1])
            for i in 2:n
                p[j] = p[j] + uwreal([0.0, A[j,i]], ids[i])
            end
        end
    end
    
    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
    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}(undef, size(der, 1))
    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
        x[i] = uwreal(mvl[i], p, d)
    end
    
    return x
end