"""
    adjust_distributions_mean(distribution_matrix::AbstractArray{T},mean_shift_percentage::Number)

Adjust the means of the distributions in `distribution_matrix` by `mean_shift_percentage`.
"""
function adjust_distributions_mean!(distribution_matrix,mean_shift_percentage) 
    for i in eachindex(distribution_matrix)
        dist = distribution_matrix[i]
        new_mean = mean(dist) + mean_shift_percentage

        distribution_matrix[i] = from_mean(typeof(dist), new_mean)

    end
end

"""
    symmetrize_means(population_sizes::Vector{Int},mixing_matrix::Matrix{<:Distribution})

Symmetrize `mixing_matrix` such that it satisifes the conditions needed for `generate_contact_vectors!(..)`, given `population_sizes`.
"""
function symmetrize_means(population_sizes,mixing_matrix::Matrix{<:Distribution})
    mixing_means = mean.(mixing_matrix)
    symmetrized_means = zeros((length(population_sizes),length(population_sizes)))

    for i in 1:length(population_sizes), j in 1:length(population_sizes)
        symmetrized_means[i,j] = 0.5*(mixing_means[i,j] + population_sizes[j]/population_sizes[i] * mixing_means[j,i])
    end
    return map(((dst,sym_mean),)->from_mean(typeof(dst),sym_mean), zip(mixing_matrix,symmetrized_means))
end

"""
This is a namedtuple with entries `hh,ws,rest` which define matrices of Poisson distributions, with parameters given by IntervalsModel.

See `load_contact_time_distributions()`.
"""
const contact_time_distributions = load_contact_time_distributions()

function shift_contact_distributions(sampler_matrix::AbstractMatrix{T},proportion) where T<:Distributions.PoissonADSampler
    return map(sampler_matrix) do sampler
        μ_old = sampler.μ
        return Poisson(μ_old*proportion)
    end
end

```
Given alphas given in `unemployment_matrix` computes the zero weight for each pair of demographic contacts given individual zero weights.
```
function alpha_matrix(alphas)
    M = zeros(length(alphas),length(alphas))
    for i in 1:length(alphas), j in 1:length(alphas)
        M[i,j] = alphas[i] + alphas[j] - alphas[j]*alphas[i]
    end
    return M
end

const unemployment_matrix = alpha_matrix(
    (
        Y = 0.11,
        M = 0.14,
        O = 0.87
    )
)

"""
Sample initial_workschool_mixing_matrix, which is the workschool distributions symmetrized for the full Canadian population, rather than subsets (as used in the ABM). This is used in IntervalsModel.
"""
@views function ws_sample(age)
    return rand.(Random.default_rng(Threads.threadid()),initial_workschool_mixing_matrix[age,:]) * (rand(Random.default_rng(Threads.threadid())) < (5/7))
end