interval_overlap_sampling.jl

# Set the underlying parameters for the intervals model
# This is the distribution of start_times
const startdist = Normal(sparam...)


function coverage!(cov,S_j,E_j)
    if E_j < S_j
        push!(cov,Interval(0,E_j))
        push!(cov,Interval(S_j,durmax))
    else
        push!(cov,Interval(S_j,E_j))
    end
end

"""

    tot_dur_sample(n::Integer,durlist::Vector{Integer})

Returns the mean of the distribution of total durations, given a set of contact durations `durlist`, and a number of samples `n`. 
Note that this is different from `tot_dur_sample!` in that it returns only a mean and does not mutate it's arguments. This could be reducible to that one but not without losing some performance I think. Used in `err_hh` where we only need a mean. 
"""
function tot_dur_sample(n,durlist)
    if isempty(durlist)
        return 0
    end
    total_dur= 0
    numcontact = length(durlist)

    int_list = Vector{Interval{Int,Closed,Closed}}()
    sizehint!(int_list,numcontact*2)
    start_matrix = trunc.(Int,rand(rng,startdist,(numcontact,n)))
    @inbounds for i in 1:n  
        empty!(int_list)
        @inbounds for j in 1:numcontact
            S_j = start_matrix[j,i] % durmax
            E_j = (S_j + durlist[j]) % durmax
            coverage!(int_list,S_j,E_j)
        end
        union!(int_list)
        total_dur += mapreduce(Intervals.span,+,int_list)
    end
    return total_dur/n 
end


"""

    tot_dur_sample(sample_list::Vector{Integer},durlist::Vector{Integer})

Mutates `sample_list` to contain samples from the distribution of total durations, given a set of contact durations `durlist`.

This is different from `tot_dur_sample` because it modifies it's arguments. Used in `err_ws` and `err_rest` where we need a distribution.
"""
function tot_dur_sample!(sample_list,durlist)
    if isempty(durlist)
        sample_list .= 0.0
        return
    end
    numcontact = length(durlist)
    n = length(sample_list)
    int_list = Vector{Interval{Int,Closed,Closed}}()
    sizehint!(int_list,numcontact*2)
    start_matrix = trunc.(Int,rand(rng,startdist,(numcontact,n)))
    for i in 1:n  
        empty!(int_list)
        for j in 1:numcontact
            S_j = start_matrix[j,i] % durmax
            E_j = (S_j + durlist[j]) % durmax
            coverage!(int_list,S_j,E_j)
        end
        union!(int_list)
        sample_list[i] = mapreduce(Intervals.span,+,int_list)
    end
end