text:
"""
```
Particle
```

The `Particle` type contains the values and weight of a given vector of parameters.

### Fields
- `weights::Vector{Float64}`: The weights of a given particle (from the perspective of importance sampling). The incremental weight that updates the weights at each step φ_n is given by p(Y|θ^i_{n-1})^(φ_n-φ_{n-1}).
- `keys::Vector{Symbol}`: The key of the parameter corresponding to a given value.
- `value::Vector{Float64}`: The matrix of particles (n_params by n_parts)
- `loglh::Float64`: The log-likelihood of the Particle
- `logpost::Float64`: The log-posterior of the Particle
- `old_loglh::Float64`: The log-likelihood of the Particle evaluated at the old data (only non-zero during time tempering)
- `old_logpost::Float64`: The log-posterior of the Particle evaluated at the old data (only non-zero during time tempering)
- `accept::Float64`: Acceptance rate of Metropolis-Hastings steps in the mutation of the particle in the previous period
"""
mutable struct Particle
    weight::Float64
    keys::Vector{Symbol}
    value::Vector{Float64}
    loglh::Float64
    logpost::Float64
    old_loglh::Float64
    accept::Float64
end

"""
```
ParticleCloud
```

The `ParticleCloud` type contains all of the relevant information for a given cloud of particles
in the SMC algorithm. Information for a single iteration is stored at any given time (and thus
the final output will be the final cloud of particles, of which only the particle values will be saved).

### Fields
- `particles::Vector{Particle}`: The vector of particles (which contain weight, keys, and value information)
- `tempering_schedule::Vector{Float64}`: The vector of ϕ_ns (tempering factors)
- `ESS::Vector{Float64}`: The vector of effective sample sizes (resample if ESS falls under the threshold)
- `stage_index::Int`: The current iteration index of the algorithm
- `n_Φ::Int`: The total number of stages of in the fixed tempering schedule
    (if the algorithm is run with an adaptive ϕ schedule then this is used to calibrate the ϕ_prop)
- `resamples::Int`: The number of times the particle population was resampled
- `c::Float64`: The mutation step size
- `accept::Float64`: The average acceptance rate of mutation steps
- `total_sampling_time::Float64`: The total amount of time that the smc algorith took to execute
"""
mutable struct ParticleCloud
    particles::Vector{Particle}
    tempering_schedule::Vector{Float64}
    ESS::Vector{Float64}
    stage_index::Int
    n_Φ::Int
    resamples::Int
    c::Float64
    accept::Float64
    total_sampling_time::Float64
end

# Implementing a sensible showing option for clouds within DataFrames
function Base.show(io::IO, cloud::ParticleCloud)
    @printf io "Particle Cloud: N_parts=%i, N_stages=%i" length(cloud) cloud.stage_index
end

"""
```
function ParticleCloud(m::AbstractDSGEModel, n_parts::Int)
```
Easier constructor for ParticleCloud, which initializes the weights to be equal, and everything else in the Particle object etc. to be empty.
"""
function ParticleCloud(m::AbstractDSGEModel, n_parts::Int)
    return ParticleCloud([Particle(1/n_parts,[m.parameters[i].key for i in 1:length(m.parameters)],
                         zeros(length(m.parameters)),0.,0.,0.,false) for n in 1:n_parts],
                         zeros(1),zeros(1),1,0,0,0.,0.25, 0.)
end

"""
```
Cloud
```

The `Cloud` type contains all of the relevant information for a given cloud of particles
in the SMC algorithm. Information for a single iteration is stored at any given time (and thus
the final output will be the final cloud of particles, of which only the particle values will be saved).

### Fields
- `particles::Matrix{Float64}`: The vector of particles (which contain weight, keys, and value information)
- `tempering_schedule::Vector{Float64}`: The vector of ϕ_ns (tempering factors)
- `ESS::Vector{Float64}`: The vector of effective sample sizes (resample if ESS falls under the threshold)
- `stage_index::Int`: The current iteration index of the algorithm
- `n_Φ::Int`: The total number of stages of in the fixed tempering schedule
    (if the algorithm is run with an adaptive ϕ schedule then this is used to calibrate the ϕ_prop)
- `resamples::Int`: The number of times the particle population was resampled
- `c::Float64`: The mutation step size
- `accept::Float64`: The average acceptance rate of mutation steps
- `total_sampling_time::Float64`: The total amount of time that the smc algorith took to execute
"""
mutable struct Cloud
    particles::Matrix{Float64}
    tempering_schedule::Vector{Float64}
    ESS::Vector{Float64}
    stage_index::Int
    n_Φ::Int
    resamples::Int
    c::Float64
    accept::Float64
    total_sampling_time::Float64
end

"""
```
function Cloud(m::AbstractDSGEModel, n_parts::Int)
```
Easier constructor for Cloud, which initializes the weights to be equal, and everything else in the Particle object etc. to be empty.
"""
function Cloud(m::AbstractDSGEModel, n_parts::Int)
    return Cloud(Matrix{Float64}(undef, n_parts, n_parameters(m) + 5),
                 zeros(1),zeros(1),1,0,0,0.,0.25, 0.)
end

"""
```
function Cloud(c::ParticleCloud)
function Cloud(c::Cloud)
```
Returns a Cloud type.
"""
function Cloud(cloud::ParticleCloud)
    return Cloud(hcat(Matrix{Float64}(get_vals(cloud)'), get_loglh(cloud), get_logprior(cloud),
                      get_old_loglh(cloud), get_accept(cloud), get_weights(cloud)),
                 cloud.tempering_schedule, cloud.ESS, cloud.stage_index, cloud.n_Φ,
                 cloud.resamples, cloud.c, cloud.accept, cloud.total_sampling_time)
end
function Cloud(cloud::Cloud)
    return cloud
end

"""
```
function ParticleCloud(cloud::Cloud, para_symbols::Vector{Symbol})
```
Returns a ParticleCloud in the place of a Cloud.
"""
function ParticleCloud(cloud::Union{SMC.Cloud, DSGE.Cloud}, para_symbols::Vector{Symbol})
    N = size(cloud.particles, 2)
    return ParticleCloud([Particle(cloud.particles[k, ind_weight(N)], para_symbols,
                                   cloud.particles[k, 1:ind_para_end(N)],
                                   cloud.particles[k, ind_loglh(N)],
                                   cloud.particles[k, ind_logpost(N)],
                                   cloud.particles[k, ind_old_loglh(N)],
                                   cloud.particles[k, ind_accept(N)]) for k in 1:length(cloud)],
                         cloud.tempering_schedule, cloud.ESS, cloud.stage_index, cloud.n_Φ,
                         cloud.resamples, cloud.c, cloud.accept, cloud.total_sampling_time)
end

"""
Find correct indices for accessing columns of cloud array.
"""
ind_para_end(N::Int)  = N-5
ind_loglh(N::Int)     = N-4
ind_logpost(N::Int)   = N-3
ind_logprior(N::Int)  = N-3 # TODO: Fix logprior/logpost shenanigans
ind_old_loglh(N::Int) = N-2
ind_accept(N::Int)    = N-1
ind_weight(N::Int)    = N

"""
```
function get_weights(c::ParticleCloud)
function get_weights(c::Cloud)
function get_weights(c::Matrix{Float64})
```
Returns Vector{Float64}(n_parts) of weights of particles in cloud.
"""
function get_weights(c::ParticleCloud)
    return map(p -> p.weight, c.particles)
end
@inline function get_weights(c::Matrix{Float64})
    return c[:, ind_weight(size(c,2))]
end
@inline function get_weights(c::Cloud)
    return c.particles[:, ind_weight(size(c.particles,2))]
end

"""
```
function get_vals(c::ParticleCloud)
function get_vals(c::Matrix{Float64})
function get_vals(c::Cloud)
```
Returns Matrix{Float64}(n_params, n_parts) of parameter values in particle cloud.
"""
function get_vals(c::ParticleCloud)
    return hcat(map(p -> p.value, c.particles)...)
end
@inline function get_vals(c::Matrix{Float64}; transpose::Bool = true)
    return transpose ? Matrix{Float64}(c[:, 1:ind_para_end(size(c, 2))]') :
                                       c[:, 1:ind_para_end(size(c, 2))]
end
@inline function get_vals(c::Cloud; transpose::Bool = true)
    return transpose ? Matrix{Float64}(c.particles[:, 1:ind_para_end(size(c.particles,2))]') :
                                       c.particles[:, 1:ind_para_end(size(c.particles,2))]
end

"""
```
function get_loglh(c::ParticleCloud)
function get_loglh(c::Matrix{Float64})
function get_loglh(c::Cloud)
```
Returns Vector{Float64}(n_parts) of log-likelihood of particles in cloud.
"""
function get_loglh(c::ParticleCloud)
    return map(p -> p.loglh, c.particles)
end
@inline function get_loglh(c::Matrix{Float64})
    return c[:, ind_loglh(size(c,2))]
end
@inline function get_loglh(c::Cloud)
    return c.particles[:, ind_loglh(size(c.particles,2))]
end

"""
```
function cloud_isempty(c::ParticleCloud)
function cloud_isempty(c::Matrix{Float64})
function cloud_isempty(c::Cloud)
```
Check if cloud has no particles.
"""
function cloud_isempty(c::ParticleCloud)
    return isempty(c.particles)
end
function cloud_isempty(c::Matrix{Float64})
    return isempty(c)
end
function cloud_isempty(c::Cloud)
    println("TODO: 'ISEMPTY' IS INCOMPATIBLE WITH INITIALIZATION")
    return isempty(c.particles)
end

"""
```
function get_old_loglh(c::ParticleCloud)
function get_old_loglh(c::Matrix{Float64})
function get_old_loglh(c::Cloud)
```
Returns Vector{Float64}(n_parts) of old log-likelihood of particles in cloud.
"""
function get_old_loglh(c::ParticleCloud)
    return map(p -> p.old_loglh, c.particles)
end
@inline function get_old_loglh(c::Matrix{Float64})
    return c[:, ind_old_loglh(size(c,2))]
end
@inline function get_old_loglh(c::Cloud)
    return c.particles[:, ind_old_loglh(size(c.particles,2))]
end

"""
```
function get_logpost(c::ParticleCloud)
function get_logpost(c::Matrix{Float64})
function get_logpost(c::Cloud)
```
Returns Vector{Float64}(n_parts) of log-posterior of particles in cloud.
"""
function get_logpost(c::ParticleCloud)
    return map(p -> p.logpost + p.loglh, c.particles)
end
@inline function get_logpost(c::Matrix{Float64})
    return c[:, ind_loglh(size(c,2))] .+ c[:, ind_logprior(size(c,2))]
end
@inline function get_logpost(c::Cloud)
    return c.particles[:, ind_loglh(size(c.particles,2))] .+
        c.particles[:, ind_logprior(size(c.particles,2))]
end

"""
```
function get_logprior(c::ParticleCloud)
function get_logprior(c::Matrix{Float64})
function get_logprior(c::Cloud)
```
Returns Vector{Float64}(n_parts) of log-prior of particles in cloud.
"""
function get_logprior(c::ParticleCloud)
    #TODO: Fix logpost/logprior confusion
    return map(p -> p.logpost, c.particles)
end
@inline function get_logprior(c::Matrix{Float64})
    return c[:, ind_logprior(size(c,2))]
end
@inline function get_logprior(c::Cloud)
    return c.particles[:, ind_logprior(size(c.particles,2))]
end

"""
```
function get_accept(c::ParticleCloud)
function get_accept(c::Matrix{Float64})
function get_accept(c::Cloud)
```
Returns Vector{Float64}(n_parts) of old log-likelihood of particles in cloud.
"""
function get_accept(c::ParticleCloud)
    return map(p -> p.accept, c.particles)
end
@inline function get_accept(c::Matrix{Float64})
    return c[:, ind_accept(size(c,2))]
end
@inline function get_accept(c::Cloud)
    return c.particles[:, ind_accept(size(c.particles,2))]
end

"""
```
function get_likeliest_particle_value(c::ParticleCloud)
function get_likeliest_particle_value(c::Matrix{Float64})
function get_likeliest_particle_value(c::Cloud)
```
Return parameter vector of particle with highest log-likelihood.
"""
function get_likeliest_particle_value(c::ParticleCloud)
    return c.particles[indmax(get_loglh(c))].value
end
@inline function get_likeliest_particle_value(c::Matrix{Float64})
    return c[indmax(get_loglh(c)), 1:ind_para_end(size(c,2))]
end
@inline function get_likeliest_particle_value(c::Cloud)
    return c.particles[indmax(get_loglh(c)), 1:ind_para_end(size(c.particles,2))]
end

"""
```
function update_draws!(c::ParticleCloud, draws::Matrix{Float64})
function update_draws!(c::ParticleCloud, draws::Array{Particle,1})
@inline function update_draws!(c::Matrix{Float64}, draws::Matrix{Float64})
```
Update parameter draws in cloud.
"""
function update_draws!(c::ParticleCloud, draws::Matrix{Float64})
    @assert size(draws) == (length(c.particles[1]), length(c))
    for (i,p) in enumerate(c.particles)
        update_val!(p, draws[:,i])
    end
end
function update_draws!(c::ParticleCloud, draws::Array{Particle,1})
    for (i,d) in enumerate(draws)
        c.particles[i] = d
    end
end
@inline function update_draws!(c::Cloud, draws::Matrix{Float64})
    I, J     = size(draws)
    n_parts  = length(c)
    n_params = ind_para_end(size(c.particles, 2))
    if (I, J) == (n_parts, n_params)
        for i = 1:I, j=1:J
            c.particles[i, j] = draws[i, j]
        end
    elseif (I, J) == (n_params, n_parts)
        for i = 1:I, j=1:J
            c.particles[j, i] = draws[i, j]
        end
    else
        throw(error("update_draws!(c::Cloud, draws::Matrix): Draws are incorrectly sized!"))
    end
end

"""
```
function update_weights!(c::ParticleCloud, incweight::Vector{Float64})
@inline function update_weights!(c::Matrix{Float64}, incweight::Vector{Float64})
@inline function update_weights!(c::Cloud, weights::Vector{Float64})
```
Update weights in cloud.
"""
function update_weights!(c::ParticleCloud, incweight::Vector{Float64})
    new_weights = get_weights(c) .* incweight
    for (p,w) in zip(c.particles, new_weights)
        update_weight!(p, w)
    end
end
@inline function update_weights!(c::Matrix{Float64}, incweight::Vector{Float64})
    @assert size(c, 1) == length(incweight) "Dimensional mismatch in inc. weights"
    N = ind_weight(size(c,2))
    for i=1:length(incweight)
        c[i, N] *= incweight[i]
    end
end
@inline function update_weights!(c::Cloud, weights::Vector{Float64})
    update_weights!(c.particles, weights)
end

"""
```
@inline function set_weights!(c::Cloud, weights::Vector{Float64})
```
Set weights to specific values. Contrast to update_weights, which multiplies
existing weights by provided incremental weights.
"""
function set_weights!(c::ParticleCloud, weights::Vector{Float64})
    for (p,w) in zip(c.particles, weights)
        update_weight!(p,w)
    end
end
@inline function set_weights!(c::Cloud, weights::Vector{Float64})
    @assert length(c) == length(weights) "Dimensional mismatch in set_weights"
    N = ind_weight(size(c.particles,2))
    for i=1:length(c)
        c.particles[i, N] = weights[i]
    end
end

"""
```
function update_loglh!(c::ParticleCloud, incweight::Vector{Float64})
@inline function update_loglh!(c::Matrix{Float64}, incweight::Vector{Float64})
@inline function update_loglh!(c::Cloud, loglh::Vector{Float64})
```
Update log-likelihood in cloud.
"""
function update_loglh!(c::ParticleCloud, loglh::Vector{Float64})
    for (p,l) in zip(c.particles,loglh)
        p.loglh = l
    end
end
@inline function update_loglh!(c::Matrix{Float64}, loglh::Vector{Float64})
    @assert size(c,1) == length(loglh) "Dimensional mismatch"
    N = ind_loglh(size(c,2))
    for i=1:length(loglh)
        c[i, N] = loglh[i]
    end
end
@inline function update_loglh!(c::Cloud, loglh::Vector{Float64})
    update_loglh!(c.particles, loglh)
end

"""
```
function update_logpost!(c::ParticleCloud, incweight::Vector{Float64})
@inline function update_logpost!(c::Matrix{Float64}, incweight::Vector{Float64})
@inline function update_logpost!(c::Cloud, logpost::Vector{Float64})
```
Update log-posterior in cloud.
"""
function update_logpost!(c::ParticleCloud, logpost::Vector{Float64})
    for (p,l) in zip(c.particles,logpost)
        p.logpost = l
    end
end
@inline function update_logpost!(c::Matrix{Float64}, logpost::Vector{Float64})
    @assert size(c, 1) == length(logpost) "Dimensional mismatch"
    N = ind_logpost(size(c,2))
    for i=1:length(logpost)
        c[i, N] = logpost[i]
    end
end
@inline function update_logpost!(c::Cloud, logpost::Vector{Float64})
    update_logpost!(c.particles, logpost)
end

"""
```
function update_old_loglh!(c::ParticleCloud, incweight::Vector{Float64})
@inline function update_old_loglh!(c::Matrix{Float64}, incweight::Vector{Float64})
@inline function update_old_loglh!(c::Cloud, old_loglh::Vector{Float64})
```
Update log-likelihood in cloud.
"""
function update_old_loglh!(c::ParticleCloud, old_loglh::Vector{Float64})
    for (p,l) in zip(c.particles,old_loglh)
        p.old_loglh = l
    end
end
@inline function update_old_loglh!(c::Matrix{Float64}, old_loglh::Vector{Float64})
    @assert size(c, 1) == length(old_loglh) "Dimensional mismatch"
    N = ind_old_loglh(size(c,2))
    for i=1:length(old_loglh)
        c[i, N] = old_loglh[i]
    end
end
@inline function update_old_loglh!(c::Cloud, old_loglh::Vector{Float64})
    update_old_loglh!(c.particles, old_loglh)
end

"""
```
function normalize_weights!(c::ParticleCloud)
function normalize_weights!(c::Matrix{Float64})
@inline function normalize_weights!(c::Cloud)
```
Normalize weights in cloud.
"""
function normalize_weights!(c::ParticleCloud)
    new_weights = get_weights(c) / sum(get_weights(c))
    for (p,w) in zip(c.particles, new_weights)
        update_weight!(p,w)
    end
end
@inline function normalize_weights!(c::Matrix{Float64})
    sum_weights = sum(get_weights(c))
    c[:, ind_weight(size(c,2))] /= sum_weights
end
@inline function normalize_weights!(c::Cloud)
    normalize_weights!(c.particles)
end

"""
```
function reset_weights!(c::ParticleCloud)
function reset_weights!(c::Matrix{Float64})
function reset_weights!(c::Cloud)
```
Uniformly reset weights of all particles to 1/n_parts.
"""
function reset_weights!(c::ParticleCloud)
    map(p->update_weight!(p,1/length(c)), c.particles)
end
@inline function reset_weights!(c::Matrix{Float64})
    n_parts = size(c, 1)
    c[:, ind_weight(size(c,2))] .= 1.0 / n_parts
end
function reset_weights!(c::Cloud)
    reset_weights!(c.particles)
end

"""
```
function update_mutation!(p::Particle, para::Vector{Float64},
                          like::Float64, post::Float64, old_like::Float64, accept::Bool)
```
Deprecated syntax. Accept should be a Float64, not a Bool.
"""
function update_mutation!(p::Particle, para::Vector{Float64},
                          like::Float64, post::Float64, old_like::Float64, accept::Bool)
    p.value = para
    p.loglh = like
    p.logpost = post
    p.old_loglh = old_like
    p.accept = accept
end

"""
```
function update_mutation!(p::Particle, para::Vector{Float64},
                          like::Float64, post::Float64, old_like::Float64, accept::Float64)
function update_mutation!(p::Vector{Float64}, para::Vector{Float64},
                          like::Float64, post::Float64, old_like::Float64, accept::Float64)
```
Update a particle's parameter vector, log-likelihood, log-posteriod, old log-likelihood, and acceptance rate at the end of mutation.
"""
function update_mutation!(p::Particle, para::Vector{Float64},
                          like::Float64, post::Float64, old_like::Float64, accept::Float64)
    p.value = para
    p.loglh = like
    p.logpost = post
    p.old_loglh = old_like
    p.accept = accept
end
@inline function update_mutation!(p::Vector{Float64}, para::Vector{Float64}, like::Float64,
                                  post::Float64, old_like::Float64, accept::Float64)
    N = length(p)
    p[1:ind_para_end(N)] = para
    p[ind_loglh(N)]      = like
    p[ind_logpost(N)]    = post
    p[ind_old_loglh(N)]  = old_like
    p[ind_accept(N)]     = accept
end

"""
```
function update_cloud!(cloud::ParticleCloud, new_particles::Matrix{Float64})
function update_cloud!(cloud::Cloud, new_particles::Matrix{Float64})
```
Updates cloud values with those of new particles in particle array.
"""
function update_cloud!(cloud::ParticleCloud, new_particles::Matrix{Float64})
    M, N = size(new_particles)
    if M == length(cloud.particles)
        for k=1:length(cloud.particles)
            cloud.particles[k].value     = new_particles[k, 1:ind_para_end(N)]
            cloud.particles[k].loglh     = new_particles[k, ind_loglh(N)]
            cloud.particles[k].logpost   = new_particles[k, ind_logpost(N)]
            cloud.particles[k].old_loglh = new_particles[k, ind_old_loglh(N)]
            cloud.particles[k].accept    = new_particles[k, ind_accept(N)]
        end
    elseif N == length(cloud.particles)
        for k=1:length(cloud.particles)
            cloud.particles[k].value     = new_particles[1:ind_para_end(M), k]
            cloud.particles[k].loglh     = new_particles[ind_loglh(M), k]
            cloud.particles[k].logpost   = new_particles[ind_logpost(M), k]
            cloud.particles[k].old_loglh = new_particles[ind_old_loglh(M), k]
            cloud.particles[k].accept    = new_particles[ind_accept(M), k]
        end
    else
        throw(error("update_cloud!(c::ParticleCloud, draws): draws are incorrect size!"))
    end
end
function update_cloud!(cloud::Cloud, new_particles::Matrix{Float64})
    I, J = size(new_particles)
    if I == length(cloud)
        cloud.particles = new_particles
    elseif J == length(cloud)
        for k = 1:length(cloud)
            cloud.particles[k,:] = new_particles[:, k]
        end
    else
        throw(error("update_cloud!(c::Cloud, draws): draws are incorrect size!"))
    end
end

"""
```
function get_cloud_array(c::ParticleCloud)
```
Returns a particle array in the place of a ParticleCloud.
"""
function get_cloud_array(c::ParticleCloud)
    return hcat(Matrix{Float64}(get_vals(c)'), get_loglh(c), get_logprior(c),
                get_old_loglh(c), get_accept(c), get_weights(c))
end

"""
```
function update_val!(p::Particle, val::Vector{Float64})
@inline function update_val!(p::Vector{Float64}, val::Vector{Float64})
```
Update parameter vector of particle.
"""
function update_val!(p::Particle, val::Vector{Float64})
    p.value = val
end
@inline function update_val!(p::Vector{Float64}, val::Vector{Float64})
    @assert ind_para_end(length(p)) == length(val) "Parameter vector length is wrong!"
    p[1:length(val)] = val
end

"""
```
function update_weight!(p::Particle, weight::Vector{Float64})
@inline function update_weight!(p::Vector{Float64}, weight::Vector{Float64})
```
Update weight of particle.
"""
function update_weight!(p::Particle, weight::Float64)
    p.weight = weight
end
@inline function update_weight!(p::Vector{Float64}, weight::Float64)
    p[ind_weight(length(p))] = weight
end

"""
```
function update_acceptance_rate!(c::ParticleCloud)
function update_acceptance_rate!(c::Cloud)
```
Update cloud's acceptance rate with the mean of its particle acceptance rates.
"""
function update_acceptance_rate!(c::ParticleCloud)
    c.accept = mean(map(p -> p.accept, c.particles))
end
function update_acceptance_rate!(c::Cloud)
    c.accept = mean(get_accept(c))
end

@inline function Base.length(c::ParticleCloud)
        return length(c.particles)
end
@inline function Base.length(c::Cloud)
        return size(c.particles, 1)
end
@inline function Base.length(p::Particle)
    return length(p.value)
end

"""
```
function weighted_mean(c::ParticleCloud)
function weighted_mean(c::Cloud)
function weighted_mean(c::Matrix{Float64})
```
Compute weighted mean of particle cloud.
"""
function weighted_mean(c::ParticleCloud)
    return get_vals(c) * get_weights(c)
end
@inline function weighted_mean(c::Matrix{Float64})
    return get_vals(c) * get_weights(c)
end
@inline function weighted_mean(c::Cloud)
    return weighted_mean(c.particles)
end

"""
```
function weighted_quantile(c::ParticleCloud, i::Int64)
@inline function weighted_quantile(c::Matrix{Float64}, i::Int64)
function weighted_quantile(c::Cloud, i::Int64)
```
Compute weighted quantiles of particle cloud for input parameter, indexed by i.
"""
function weighted_quantile(c::ParticleCloud, i::Int64)
    lb = quantile(get_vals(c)[i, :], Weights(get_weights(c)), .05)
    ub = quantile(get_vals(c)[i, :], Weights(get_weights(c)), .95)
    return lb, ub
end
@inline function weighted_quantile(c::Matrix{Float64}, i::Int64)
    @assert i <= ind_end_para(size(c,2)) "Parameter index invalid."
    lb = quantile(c[:, i], Weights(get_weights(c)), .05)
    ub = quantile(c[:, i], Weights(get_weights(c)), .95)
    return lb, ub
end
function weighted_quantile(c::Cloud, i::Int64)
    return weighted_quantile(c.particles, i)
end

"""
```
function weighted_std(c::ParticleCloud)
function weighted_std(c::Cloud)
@inline function weighted_std(c::Matrix{Float64})
```
Compute weighted standard deviation of particle cloud.
"""
function weighted_std(c::ParticleCloud)
    return sqrt.(diag(weighted_cov(c)))
end
@inline function weighted_std(c::Matrix{Float64})
    return sqrt.(diag(weighted_cov(c)))
end
function weighted_std(c::Cloud)
    return weighted_std(c.particles)
end

"""
```
function weighted_cov(c::ParticleCloud)
function weighted_cov(c::Cloud)
@inline function weighted_cov(c::Matrix{Float64})
```
Compute weighted covariance of particle cloud.
"""
function weighted_cov(c::ParticleCloud)
    return cov(Matrix(get_vals(c)'), Weights(get_weights(c)), corrected = false)
end
@inline function weighted_cov(c::Matrix{Float64})
    return cov(get_vals(c; transpose = false), Weights(get_weights(c)), corrected = false)
end
@inline function weighted_cov(c::Cloud)
    return weighted_cov(c.particles)
end

for op in (:(Base.:+),
           :(Base.:-),
           :(Base.:*),
           :(Base.:/),
           :(Base.:^))

    @eval ($op)(p::Particle, q::Particle)   = ($op)(p.value,q.value)
    @eval ($op)(p::Particle, x::Integer)    = ($op)(p.value,x)
    @eval ($op)(p::Particle, x::Number)     = ($op)(p.value,x)
    @eval ($op)(x::Integer, p::Particle)    = ($op)(p.value,x)
    @eval ($op)(x::Number, p::Particle)     = ($op)(p.value,x)
    @eval ($op)(p::Particle, q::Vector{Float64})   = @. ($op)(p.value,q)
    @eval ($op)(p::Vector{Float64}, q::Particle)   = @. ($op)(p,q.value)
end


text:
using Base.Broadcast: broadcastable
using ChainRules
using ChainRulesCore
using FiniteDifferences
using LinearAlgebra
using LinearAlgebra.BLAS
using LinearAlgebra: dot
using Random
using Statistics
using Test

# For testing purposes we use a lot of
using ChainRulesCore: cast, extern, accumulate, accumulate!, store!, @scalar_rule,
    Wirtinger, wirtinger_primal, wirtinger_conjugate,
    Zero, One, Casted, DNE, Thunk, AbstractDifferential

include("test_util.jl")

println("Testing ChainRules.jl")
@testset "ChainRules" begin
    include("helper_functions.jl")
    @testset "rulesets" begin
        @testset "Base" begin
            include(joinpath("rulesets", "Base", "base.jl"))
            include(joinpath("rulesets", "Base", "array.jl"))
            include(joinpath("rulesets", "Base", "mapreduce.jl"))
            include(joinpath("rulesets", "Base", "broadcast.jl"))
        end

        print(" ")

        @testset "Statistics" begin
            include(joinpath("rulesets", "Statistics", "statistics.jl"))
        end

        print(" ")

        @testset "LinearAlgebra" begin
            include(joinpath("rulesets", "LinearAlgebra", "dense.jl"))
            include(joinpath("rulesets", "LinearAlgebra", "structured.jl"))
            include(joinpath("rulesets", "LinearAlgebra", "factorization.jl"))
            include(joinpath("rulesets", "LinearAlgebra", "blas.jl"))
        end

        print(" ")

        @testset "packages" begin
            include(joinpath("rulesets", "packages", "NaNMath.jl"))
            include(joinpath("rulesets", "packages", "SpecialFunctions.jl"))
        end
    end
end

ext: jl


text:
cd(@__DIR__)
ks, ctr, face, simTime = KitBase.initialize("config_1d0f.txt")

using KitBase.JLD2
set = ks
t = 0.0
cd(@__DIR__)
@save "restart.jld2" set ctr t
ks, ctr, face, simTime = KitBase.initialize("restart.jld2")

dt = KitBase.timestep(ks, ctr, 0.0)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :period)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :extra)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :balance)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :mirror)
t = KitBase.solve!(ks, ctr, face, 0.0)
KitBase.evolve!(ks, ctr, face, dt; mode = :roe)

ks, ctr, face, simTime = KitBase.initialize("config_1d0f2s.txt")
dt = KitBase.timestep(ks, ctr, 0.0)
KitBase.update!(ks, ctr, face, dt, zeros(3, 2); bc = :extra)
KitBase.update!(ks, ctr, face, dt, zeros(3, 2); bc = :period)
KitBase.update!(ks, ctr, face, dt, zeros(3, 2); bc = :balance)
KitBase.update!(ks, ctr, face, dt, zeros(3, 2); bc = :mirror)
KitBase.evolve!(ks, ctr, face, dt; mode = :gks)

ks, ctr, face, simTime = KitBase.initialize("config_1d1f1v.txt")
KitBase.reconstruct!(ks, ctr)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :extra)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :period)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :balance)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :mirror)
KitBase.evolve!(ks, ctr, face, dt; mode = :kfvs)
KitBase.evolve!(ks, ctr, face, dt; mode = :kcu)

ks, ctr, face, simTime = KitBase.initialize("config_1d1f3v.txt")
KitBase.reconstruct!(ks, ctr)
KitBase.update!(ks, ctr, face, dt, zeros(5); bc = :extra)
KitBase.update!(ks, ctr, face, dt, zeros(5); bc = :period)
KitBase.update!(ks, ctr, face, dt, zeros(5); bc = :balance)
KitBase.update!(ks, ctr, face, dt, zeros(5); bc = :mirror)
KitBase.evolve!(ks, ctr, face, dt; mode = :kfvs)
KitBase.evolve!(ks, ctr, face, dt; mode = :kcu)

ks, ctr, face, simTime = KitBase.initialize("config_1d2f.txt")
KitBase.reconstruct!(ks, ctr)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :extra)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :period)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :balance)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :mirror)
KitBase.evolve!(ks, ctr, face, dt; mode = :kfvs)
KitBase.evolve!(ks, ctr, face, dt; mode = :kcu)

gas = KitBase.Gas()
ks1 = KitBase.SolverSet(ks.set, ks.pSpace, ks.vSpace, gas, ks.ib, ks.outputFolder)
KitBase.init_ptc!(ks1, ctr, mode = :soa)
KitBase.init_ptc!(ks1, ctr, mode = :aos)

ks, ctr, face, simTime = KitBase.initialize("config_1d2f2s.txt")
KitBase.reconstruct!(ks, ctr)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :extra)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :period)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :balance)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :mirror)
KitBase.evolve!(ks, ctr, face, dt; mode = :kfvs)
KitBase.evolve!(ks, ctr, face, dt; mode = :kcu)

ks, ctr, face, simTime = KitBase.initialize("config_1d4f2s.txt")
KitBase.reconstruct!(ks, ctr)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :extra)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :period)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :balance)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :mirror)
KitBase.evolve!(ks, ctr, face, dt; mode = :kfvs)
KitBase.evolve!(ks, ctr, face, dt; mode = :kcu)

ks, ctr, face, simTime = KitBase.initialize("config_1d3f2s.txt")
KitBase.reconstruct!(ks, ctr)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :extra)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :period)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :balance)
KitBase.update!(ks, ctr, face, dt, zeros(3); bc = :mirror)
KitBase.evolve!(ks, ctr, face, dt; mode = :kfvs)
KitBase.evolve!(ks, ctr, face, dt; mode = :kcu)


text:
# plot first 5000 time points (before stabilization)
include("../Model_scripts/Functions.jl");
# plot the number of infected & uninfected hosts
plot_population_numbers(
    hcat(labels...),
    "c=4.0\n \nNumber of infected and uninfected hosts",
    100,
    5000,
    "$figure_directory/supplemental/supp_2b.png"
)
# plot the average weighted R0
plot_reproductive_rate(5, 5000, "$figure_directory/supplemental/supp_2d.png")
# plot the average V0
plot_vertical_reproductive_ratio(1, 5000, "$figure_directory/supplemental/supp_2f.png")
# plot the average weighted horizontal transmission
plot_horizontal_and_virulence(1, 5000, "$figure_directory/supplemental/supp_2h.png")
# plot the evenness through time
plot_evenness(1, 5000, "$figure_directory/supplemental/supp_2j.png")


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Deprecations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


text:
#=
real symmetric or complex Hermitian positive semidefinite matrices (svec scaled triangle)
W ⪰ 0
self-dual
=#

mutable struct PosSemidefTri{C <: RealCompF} <: Cache
    oa_s::Vector{AE}
    d::Int
    w_temp::Vector{RealF}
    W_temp::Matrix{C}
    PosSemidefTri{C}() where {C <: RealCompF} = new{C}()
end

function MOIPajarito.Cones.create_cache(
    oa_s::Vector{AE},
    cone::Hypatia.PosSemidefTriCone{RealF, C},
    ::Optimizer,
) where {C <: RealCompF}
    cache = PosSemidefTri{C}()
    cache.oa_s = oa_s
    dim = MOI.dimension(cone)
    d = cache.d = svec_side(C, dim)
    cache.w_temp = zeros(RealF, dim)
    cache.W_temp = zeros(C, d, d)
    return cache
end

function MOIPajarito.Cones.add_init_cuts(cache::PosSemidefTri{C}, opt::Optimizer) where {C}
    # add variable bounds W_ii ≥ 0
    d = cache.d
    w = cache.oa_s
    JuMP.@constraint(opt.oa_model, [i in 1:d], w[svec_idx(C, i, i)] >= 0)
    opt.use_init_fixed_oa || return

    # add cuts (1, 1, ±2) on (W_ii, W_jj, W_ij), ∀i != j
    # (a linearization of W_ii * W_jj ≥ W_ij^2)
    # initial OA polyhedron is the dual cone of diagonally dominant matrices
    for j in 1:d, i in 1:(j - 1)
        _add_init_cuts(i, j, cache, opt)
    end
    return
end

# real: cuts on (w_ii, w_jj, w_ij) are (1, 1, ±rt2), ∀i != j
function _add_init_cuts(i::Int, j::Int, cache::PosSemidefTri{RealF}, opt::Optimizer)
    w = cache.oa_s
    w_ii = w[svec_idx(RealF, i, i)]
    w_jj = w[svec_idx(RealF, j, j)]
    w_ij = w[svec_idx(RealF, i, j)]
    JuMP.@constraints(opt.oa_model, begin
        w_ii + w_jj + rt2 * w_ij >= 0
        w_ii + w_jj - rt2 * w_ij >= 0
    end)
    return
end

# complex: cuts on smat (w_ii, w_jj, wr_ij, wi_ij) are (1, 1, ±1, ±1), ∀i != j
function _add_init_cuts(i::Int, j::Int, cache::PosSemidefTri{CompF}, opt::Optimizer)
    w = cache.oa_s
    w_ii = w[svec_idx(CompF, i, i)]
    w_jj = w[svec_idx(CompF, j, j)]
    ij_idx = svec_idx(CompF, i, j)
    wr_ij = w[ij_idx]
    wi_ij = w[ij_idx + 1]
    JuMP.@constraints(opt.oa_model, begin
        w_ii + w_jj + wr_ij + wi_ij >= 0
        w_ii + w_jj + wr_ij - wi_ij >= 0
        w_ii + w_jj - wr_ij + wi_ij >= 0
        w_ii + w_jj - wr_ij - wi_ij >= 0
    end)
    return
end

function MOIPajarito.Cones.get_subp_cuts(
    z::Vector{RealF},
    cache::PosSemidefTri,
    opt::Optimizer,
)
    # strengthened cuts from eigendecomposition are λᵢ * rᵢ * rᵢ'
    R = cache.W_temp
    svec_to_smat!(R, z, rt2)
    F = eigen(Hermitian(R, :U), 1e-8, Inf) # TODO tune
    isempty(F.values) && return AE[]
    R_eig = F.vectors * Diagonal(sqrt.(F.values))
    return _get_psd_cuts(R_eig, cache.oa_s, cache, opt)
end

function MOIPajarito.Cones.get_sep_cuts(
    s::Vector{RealF},
    cache::PosSemidefTri,
    opt::Optimizer,
)
    Ws = cache.W_temp
    svec_to_smat!(Ws, s, rt2)
    F = eigen(Hermitian(Ws, :U), -Inf, -opt.tol_feas)
    isempty(F.values) && return AE[]
    return _get_psd_cuts(F.vectors, cache.oa_s, cache, opt)
end


text:
# Compute the base SF index and verify same answer with culnerability model
SF = findstochasticflexibility(m, num_pts = 100, seed = 42)
@test abs(findstochasticflexibility(m, num_pts = 100, use_vulnerability_model = true, seed = 42) - SF) <= 1e-5
@test_logs (:warn, "Flexibility index hasn't yet been computed. Setting use_flexibility_index = false.") abs(findstochasticflexibility(m, num_pts = 100, use_flexibility_index = true, seed = 42) - SF) <= 1e-5

# Compute the SF with all positive values and verify both methods yield the same solution
SF2 = findstochasticflexibility(m, num_pts = 100, only_positive = true, seed = 42)
@test abs(findstochasticflexibility(m, num_pts = 100, use_vulnerability_model = true, only_positive = true, seed = 42) - SF2) <= 1e-5


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(
        full_bridge_optimizer,
        (MOI.AbstractOptimizer, Type{Float64}),
    )
    Base.precompile(
        unbridged_function,
        (LazyBridgeOptimizer, MOI.AbstractScalarFunction),
    )   # time: 0.1381832
    return
end

_precompile_()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

eltype(::Type{<:AbstractSet{T}}) where {T} = @isdefined(T) ? T : Any
sizehint!(s::AbstractSet, n) = nothing

copy!(dst::AbstractSet, src::AbstractSet) = union!(empty!(dst), src)

## set operations (union, intersection, symmetric difference)

"""
    union(s, itrs...)
    ∪(s, itrs...)

Construct the union of sets. Maintain order with arrays.

# Examples
```jldoctest
julia> union([1, 2], [3, 4])
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> union([1, 2], [2, 4])
3-element Array{Int64,1}:
 1
 2
 4

julia> union([4, 2], 1:2)
3-element Array{Int64,1}:
 4
 2
 1

julia> union(Set([1, 2]), 2:3)
Set{Int64} with 3 elements:
  2
  3
  1
```
"""
function union end

_in(itr) = x -> x in itr

union(s, sets...) = union!(emptymutable(s, promote_eltype(s, sets...)), s, sets...)
union(s::AbstractSet) = copy(s)

const ∪ = union

"""
    union!(s::Union{AbstractSet,AbstractVector}, itrs...)

Construct the union of passed in sets and overwrite `s` with the result.
Maintain order with arrays.

# Examples
```jldoctest
julia> a = Set([1, 3, 4, 5]);

julia> union!(a, 1:2:8);

julia> a
Set{Int64} with 5 elements:
  7
  4
  3
  5
  1
```
"""
function union!(s::AbstractSet, sets...)
    for x in sets
        union!(s, x)
    end
    return s
end

max_values(::Type) = typemax(Int)
max_values(T::Union{map(X -> Type{X}, BitIntegerSmall_types)...}) = 1 << (8*sizeof(T))
# saturated addition to prevent overflow with typemax(Int)
max_values(T::Union) = max(max_values(T.a), max_values(T.b), max_values(T.a) + max_values(T.b))
max_values(::Type{Bool}) = 2
max_values(::Type{Nothing}) = 1

function union!(s::AbstractSet{T}, itr) where T
    haslength(itr) && sizehint!(s, length(s) + length(itr))
    for x in itr
        push!(s, x)
        length(s) == max_values(T) && break
    end
    return s
end

"""
    intersect(s, itrs...)
    ∩(s, itrs...)

Construct the intersection of sets.
Maintain order with arrays.

# Examples
```jldoctest
julia> intersect([1, 2, 3], [3, 4, 5])
1-element Array{Int64,1}:
 3

julia> intersect([1, 4, 4, 5, 6], [4, 6, 6, 7, 8])
2-element Array{Int64,1}:
 4
 6

julia> intersect(Set([1, 2]), BitSet([2, 3]))
Set{Int64} with 1 element:
  2
```
"""
intersect(s::AbstractSet, itr, itrs...) = intersect!(intersect(s, itr), itrs...)
intersect(s) = union(s)
intersect(s::AbstractSet, itr) = mapfilter(_in(s), push!, itr, emptymutable(s))

const ∩ = intersect

"""
    intersect!(s::Union{AbstractSet,AbstractVector}, itrs...)

Intersect all passed in sets and overwrite `s` with the result.
Maintain order with arrays.
"""
function intersect!(s::AbstractSet, itrs...)
    for x in itrs
        intersect!(s, x)
    end
    return s
end
intersect!(s::AbstractSet, s2::AbstractSet) = filter!(_in(s2), s)
intersect!(s::AbstractSet, itr) =
    intersect!(s, union!(emptymutable(s, eltype(itr)), itr))

"""
    setdiff(s, itrs...)

Construct the set of elements in `s` but not in any of the iterables in `itrs`.
Maintain order with arrays.

# Examples
```jldoctest
julia> setdiff([1,2,3], [3,4,5])
2-element Array{Int64,1}:
 1
 2
```
"""
setdiff(s::AbstractSet, itrs...) = setdiff!(copymutable(s), itrs...)
setdiff(s) = union(s)

"""
    setdiff!(s, itrs...)

Remove from set `s` (in-place) each element of each iterable from `itrs`.
Maintain order with arrays.

# Examples
```jldoctest
julia> a = Set([1, 3, 4, 5]);

julia> setdiff!(a, 1:2:6);

julia> a
Set{Int64} with 1 element:
  4
```
"""
function setdiff!(s::AbstractSet, itrs...)
    for x in itrs
        setdiff!(s, x)
    end
    return s
end
function setdiff!(s::AbstractSet, itr)
    for x in itr
        delete!(s, x)
    end
    return s
end

"""
    symdiff(s, itrs...)

Construct the symmetric difference of elements in the passed in sets.
When `s` is not an `AbstractSet`, the order is maintained.
Note that in this case the multiplicity of elements matters.

# Examples
```jldoctest
julia> symdiff([1,2,3], [3,4,5], [4,5,6])
3-element Array{Int64,1}:
 1
 2
 6

julia> symdiff([1,2,1], [2, 1, 2])
2-element Array{Int64,1}:
 1
 2

julia> symdiff(unique([1,2,1]), unique([2, 1, 2]))
0-element Array{Int64,1}
```
"""
symdiff(s, sets...) = symdiff!(emptymutable(s, promote_eltype(s, sets...)), s, sets...)
symdiff(s) = symdiff!(copy(s))

"""
    symdiff!(s::Union{AbstractSet,AbstractVector}, itrs...)

Construct the symmetric difference of the passed in sets, and overwrite `s` with the result.
When `s` is an array, the order is maintained.
Note that in this case the multiplicity of elements matters.
"""
function symdiff!(s::AbstractSet, itrs...)
    for x in itrs
        symdiff!(s, x)
    end
    return s
end

function symdiff!(s::AbstractSet, itr)
    for x in itr
        x in s ? delete!(s, x) : push!(s, x)
    end
    return s
end

## non-strict subset comparison

const ⊆ = issubset
function ⊇ end
"""
    issubset(a, b) -> Bool
    ⊆(a, b) -> Bool
    ⊇(b, a) -> Bool

Determine whether every element of `a` is also in `b`, using [`in`](@ref).

# Examples
```jldoctest
julia> issubset([1, 2], [1, 2, 3])
true

julia> [1, 2, 3] ⊆ [1, 2]
false

julia> [1, 2, 3] ⊇ [1, 2]
true
```
"""
issubset, ⊆, ⊇

const FASTIN_SET_THRESHOLD = 70

function issubset(l, r)
    if haslength(r) && (isa(l, AbstractSet) || !hasfastin(r))
        rlen = length(r) # conditions above make this length computed only when needed
        # check l for too many unique elements
        if isa(l, AbstractSet) && length(l) > rlen
            return false
        end
        # when `in` would be too slow and r is big enough, convert it to a Set
        # this threshold was empirically determined (cf. #26198)
        if !hasfastin(r) && rlen > FASTIN_SET_THRESHOLD
            return issubset(l, Set(r))
        end
    end
    for elt in l
        elt in r || return false
    end
    return true
end

"""
    hasfastin(T)

Determine whether the computation `x ∈ collection` where `collection::T` can be considered
as a "fast" operation (typically constant or logarithmic complexity).
The definition `hasfastin(x) = hasfastin(typeof(x))` is provided for convenience so that instances
can be passed instead of types.
However the form that accepts a type argument should be defined for new types.
"""
hasfastin(::Type) = false
hasfastin(::Union{Type{<:AbstractSet},Type{<:AbstractDict},Type{<:AbstractRange}}) = true
hasfastin(x) = hasfastin(typeof(x))

⊇(l, r) = r ⊆ l

## strict subset comparison

function ⊊ end
function ⊋ end
"""
    ⊊(a, b) -> Bool
    ⊋(b, a) -> Bool

Determines if `a` is a subset of, but not equal to, `b`.

# Examples
```jldoctest
julia> (1, 2) ⊊ (1, 2, 3)
true

julia> (1, 2) ⊊ (1, 2)
false
```
"""
⊊, ⊋

⊊(l::AbstractSet, r) = length(l) < length(r) && l ⊆ r
⊊(l, r) = Set(l) ⊊ r
⊋(l, r) = r ⊊ l

function ⊈ end
function ⊉ end
"""
    ⊈(a, b) -> Bool
    ⊉(b, a) -> Bool

Negation of `⊆` and `⊇`, i.e. checks that `a` is not a subset of `b`.

# Examples
```jldoctest
julia> (1, 2) ⊈ (2, 3)
true

julia> (1, 2) ⊈ (1, 2, 3)
false
```
"""
⊈, ⊉

⊈(l, r) = !⊆(l, r)
⊉(l, r) = r ⊈ l

## set equality comparison

"""
    issetequal(a, b) -> Bool

Determine whether `a` and `b` have the same elements. Equivalent
to `a ⊆ b && b ⊆ a` but more efficient when possible.

# Examples
```jldoctest
julia> issetequal([1, 2], [1, 2, 3])
false

julia> issetequal([1, 2], [2, 1])
true
```
"""
issetequal(l::AbstractSet, r::AbstractSet) = l == r
issetequal(l::AbstractSet, r) = issetequal(l, Set(r))

function issetequal(l, r::AbstractSet)
    if haslength(l)
        # check r for too many unique elements
        length(l) < length(r) && return false
    end
    return issetequal(Set(l), r)
end

function issetequal(l, r)
    haslength(l) && return issetequal(l, Set(r))
    haslength(r) && return issetequal(r, Set(l))
    return issetequal(Set(l), Set(r))
end

## set disjoint comparison
"""
    isdisjoint(v1, v2) -> Bool

Return whether the collections `v1` and `v2` are disjoint, i.e. whether
their intersection is empty.

!!! compat "Julia 1.5"
    This function requires at least Julia 1.5.
"""
function isdisjoint(l, r)
    function _isdisjoint(l, r)
        hasfastin(r) && return !any(in(r), l)
        hasfastin(l) && return !any(in(l), r)
        haslength(r) && length(r) < FASTIN_SET_THRESHOLD &&
            return !any(in(r), l)
        return !any(in(Set(r)), l)
    end
    if haslength(l) && haslength(r) && length(r) < length(l)
        return _isdisjoint(r, l)
    end
    _isdisjoint(l, r)
end

## partial ordering of sets by containment

==(l::AbstractSet, r::AbstractSet) = length(l) == length(r) && l ⊆ r
# convenience functions for AbstractSet
# (if needed, only their synonyms ⊊ and ⊆ must be specialized)
<( l::AbstractSet, r::AbstractSet) = l ⊊ r
<=(l::AbstractSet, r::AbstractSet) = l ⊆ r

## filtering sets

filter(pred, s::AbstractSet) = mapfilter(pred, push!, s, emptymutable(s))

# it must be safe to delete the current element while iterating over s:
unsafe_filter!(pred, s::AbstractSet) = mapfilter(!pred, delete!, s, s)

# TODO: delete mapfilter in favor of comprehensions/foldl/filter when competitive
function mapfilter(pred, f, itr, res)
    for x in itr
        pred(x) && f(res, x)
    end
    res
end

max_stars_repo_name: Dictino/julia

max_issues_repo_name: Dictino/julia

max_forks_repo_name: Dictino/julia


text:
# remove ++++++
# using CollegeStratData, Test

# csd = CollegeStratData;

# function school_groups_test(sg); # :: csd.AbstractSchoolGroups)
#     @testset "School groups" begin
#         sGroupV = csd.s_groups(sg);
#         @test length(sGroupV) == csd.n_school(sg)
#         for j = 1 : csd.n_school(sg)
#             @test csd.s_idx(sg, sGroupV[j]) == j
#         end
# 	end
# end

# @testset "SchoolGroups" begin
#     school_groups_test(csd.SchoolGroups3());
# end

# ------------


text:
using NicePipes
using Test

@static if Sys.iswindows()
    const LE = "\r\n"
else
    const LE = "\n"
end

function test_show(x, show_x)
    io = IOBuffer()
    show(io, x)
    output = String(take!(io))
    output = replace(output, LE*"\e[1A" => "")
    @test output == show_x
end

@testset "NicePipes.jl" begin
    a = ["foo", "bar"]
    show_a = [
        "2-element $(Vector{String}):",
        " \"foo\"",
        " \"bar\"",
    ]
    test_show((a | @grep foo), show_a[2])
    test_show((a | @grep -iv FoO), join(show_a[[1, 3]], LE))
    test_show((3 | @grep 4), "No matches found!\n")
    test_show((a | @sed "/foo/d"), join(show_a[[1, 3]], LE))
    test_show((a | @sed raw"s/f\(o\+\)/b\1/g"), show_a[1] * LE * " \"boo\"" * LE * show_a[3])
end


text:
@testset "Lanczos - Exponentiate full" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T,(n,n)) .- one(T)/2
            A = (A+A')/2
            V = one(A)
            W = zero(A)
            alg = Lanczos(orth = orth, krylovdim = n, maxiter = 1, tol = 10*n*eps(real(T)))
            for k = 1:n
                W[:,k], =  @inferred exponentiate(A, 1, view(V,:,k), alg)
            end
            @test W ≈ exp(A)
        end
    end
end

@testset "Lanczos - Exponentiate iteratively" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T,(N,N)) .- one(T)/2
            A = (A+A')/2
            v = rand(T,N)
            t = rand(complex(T))
            alg = Lanczos(orth = orth, krylovdim = n, maxiter = 30, tol = 10*N*eps(real(T))*abs(t))
            w, info =  @inferred exponentiate(A, t, v, alg)
            @test info.converged > 0
            @test w ≈ exp(t*A)*v
        end
    end
end


text:
include("Graph.jl")
include("SimMetric.jl")

using SparseArrays
using LinearAlgebra
import Laplacians

function getSimSelectorByS(S)
    n = size(S, 1)
    f(u, k) = begin
        ret = Vector{Int32}()
        A = Vector{Tuple{Int32, Float64}}()
        for i = 1 : n
            if i != u
                push!(A, (i, S[u, i]))
            end
        end
        sort!(A, by=x->-x[2])
        for i = 1 : k
            push!(ret, A[i][1])
        end
        return ret
    end
    return f
end

function RoleSimSelector(G)
    S = RoleSim(G)
    return getSimSelectorByS(S)
end

function ForestSimSelector(G)
    S = ForestSim(G)
    return getSimSelectorByS(S)
end

function ForestSimApproxSelector(G)
    S = approxForest(G)
    A = Vector{Tuple{Int32, Float64}}()
    foreach(i -> push!(A, (i, S[i])), 1 : G.n)
    sort!(A, by=x->x[2])
    pos = zeros(Int32, G.n)
    foreach(i -> pos[A[i][1]] = i, 1 : G.n)

    sim(x, y) = min(S[x], S[y]) / max(S[x], S[y])

    f(x, k) = begin
        l = pos[x]-1
        r = pos[x]+1
        ret = Vector{Int32}()
        for rep = 1 : k
            if l == 0
                push!(ret, A[r][1])
                r += 1
            elseif r > G.n
                push!(ret, A[l][1])
                l -= 1
            elseif sim(x, A[l][1]) > sim(x, A[r][1])
                push!(ret, A[l][1])
                l -= 1
            else
                push!(ret, A[r][1])
                r += 1
            end
        end
        return ret
    end

    return f
end

function StructSimSelector(G; maxk = 3, wi = 0.5)
    n = G.n
    d = degrees(G)
    B, dk = calcBins(G, maxk)
    maxb = size(B, 3)

    f(x, k) = begin
        S = zeros(n)
        foreach(i -> S[i] = min(d[i], d[x]) / max(d[i], d[x]), 1 : n)
        for k = 1 : maxk
            for i = 1 : n
                delta = 0.0
                foreach(z -> delta += min(B[k, i, z], B[k, x, z]), 1 : maxb)
                delta /= max(dk[k, i], dk[k, x])
                S[i] = (1.0 - wi) * S[i] + wi * delta
            end
        end
        A = Vector{Tuple{Int32, Float64}}()
        for i = 1 : n
            if i != x
                push!(A, (i, S[i]))
            end
        end
        sort!(A, by=x->-x[2])
        ret = Vector{Int32}()
        for i = 1 : k
            push!(ret, A[i][1])
        end
        return ret
    end

    return f
end


text:
using SimpleSDMLayers
using Test

global anyerrors = false

tests = [
   "construction" => "core/construction.jl",
   "lat./lon." => "core/latlon.jl",
   "lat./lon. conversion" => "core/coordconvert.jl",
   "iteration" => "core/iteration.jl",
   "tiling" => "core/tiling.jl",
   "setindex" => "core/setindex.jl",
   "overloads" => "overloads.jl",
   "clipping" => "operations/clip.jl",
   "rescale" => "operations/rescale.jl",
   "mosaic" => "operations/mosaic.jl",
   "coarsen" => "operations/coarsen.jl",
   "subsetting" => "operations/subsetting.jl",
   "generated" => "generated.jl",
   "ascii" => "data/ascii.jl",
   "import" => "data/dataread.jl",
   "worldclim" => "data/worldclim.jl",
   "earthenv" => "data/earthenv.jl",
   "chelsa" => "data/chelsa.jl",
   "plotting" => "extensions/plots.jl",
   "GBIF" => "extensions/gbif.jl",
]

for test in tests
   try
      include(test.second)
      println("\033[1m\033[32m✓\033[0m\t$(test.first)")
   catch e
      global anyerrors = true
      println("\033[1m\033[31m×\033[0m\t$(test.first)")
      println("\033[1m\033[38m→\033[0m\ttest/$(test.second)")
      showerror(stdout, e, backtrace())
      println()
      break
   end
end

if anyerrors
   throw("Tests failed")
end


text:
# import TensorOperations: checked_similar_from_indices, scalar, isblascontractable, add!, contract!, unsafe_contract!

# Add support for TensorOperations.jl
function checked_similar_from_indices(tC, ::Type{T}, p1::IndexTuple{N₁}, p2::IndexTuple{N₂},
        t::AbstractTensorMap{S}) where {T,S<:IndexSpace,N₁,N₂}

    cod = ProductSpace{S,N₁}(space.(Ref(t), p1))
    dom = ProductSpace{S,N₂}(dual.(space.(Ref(t), p2)))
    A = similarstoragetype(t, T)

    if tC !== nothing && (tC isa TensorMap{S,N₁,N₂}) && eltype(tC) == T &&
        codomain(tC) == cod && domain(tC) == dom && storagetype(tC) == A
        G = sectortype(S)
        if G === Trivial
            return tC::TensorMap{S,N₁,N₂,G,A,Nothing,Nothing}
        else
            F₁ = fusiontreetype(G, StaticLength(N₁))
            F₂ = fusiontreetype(G, StaticLength(N₂))
            return tC::TensorMap{S,N₁,N₂,G,SectorDict{G,A},F₁,F₂}
        end
    else
        tC = similar(t, T, cod←dom)
        return tC
    end
end

function checked_similar_from_indices(tC, ::Type{T}, oindA::IndexTuple, oindB::IndexTuple,
    p1::IndexTuple{N₁}, p2::IndexTuple{N₂}, tA::AbstractTensorMap{S},
    tB::AbstractTensorMap{S}) where {T, S<:IndexSpace,N₁,N₂}

    spaces = (space.(Ref(tA), oindA)..., space.(Ref(tB), oindB)...)
    cod = ProductSpace{S,N₁}(getindex.(Ref(spaces), p1))
    dom = ProductSpace{S,N₂}(dual.(getindex.(Ref(spaces), p2)))
    A = similarstoragetype(tA, T)
    @assert A === similarstoragetype(tB, T)

    if tC !== nothing && (tC isa TensorMap{S,N₁,N₂}) && eltype(tC) == T &&
        codomain(tC) == cod && domain(tC) == dom && storagetype(tC) == A
        G = sectortype(S)
        if G === Trivial
            return tC::TensorMap{S,N₁,N₂,G,A,Nothing,Nothing}
        else
            F₁ = fusiontreetype(G, StaticLength(N₁))
            F₂ = fusiontreetype(G, StaticLength(N₂))
            return tC::TensorMap{S,N₁,N₂,G,SectorDict{G,A},F₁,F₂}
        end
    else
        tC = similar(tA, T, cod←dom)
        return tC
    end
end

function cached_permuteind(sym::Symbol, t::TensorMap{S},
                            p1::IndexTuple{N₁},  p2::IndexTuple{N₂}=()) where {S,N₁,N₂}
    cod = ProductSpace{S,N₁}(map(n->space(t, n), p1))
    dom = ProductSpace{S,N₂}(map(n->dual(space(t, n)), p2))

    # share data if possible
    if p1 === codomainind(t) && p2 === domainind(t)
        return t
    elseif isa(t, TensorMap) && sectortype(S) === Trivial
        stridet = i->stride(t[], i)
        sizet = i->size(t[], i)
        canfuse1, d1, s1 = TensorOperations._canfuse(sizet.(p1), stridet.(p1))
        canfuse2, d2, s2 = TensorOperations._canfuse(sizet.(p2), stridet.(p2))
        if canfuse1 && canfuse2 && s1 == 1 && (d2 == 1 || s2 == d1)
            return TensorMap(reshape(t.data, dim(cod), dim(dom)), cod, dom)
        end
    end
    # general case
    @inbounds begin
        tp = TO.cached_similar_from_indices(sym, eltype(t), p1, p2, t, :N)
        return permuteind!(tp, t, p1, p2)
    end
end

function cached_permuteind(sym::Symbol, t::AdjointTensorMap{S},
                            p1::IndexTuple{N₁},  p2::IndexTuple{N₂}=()) where {S,N₁,N₂}

    p1′ = adjointtensorindices(t, p2)
    p2′ = adjointtensorindices(t, p1)
    adjoint(cached_permuteind(sym, adjoint(t), p1′, p2′))
end

scalar(t::AbstractTensorMap{S}) where {S<:IndexSpace} =
    dim(codomain(t)) == dim(domain(t)) == 1 ?
        first(blocks(t))[2][1,1] : throw(SpaceMismatch())

function add!(α, tsrc::AbstractTensorMap{S}, β, tdst::AbstractTensorMap{S,N₁,N₂},
     p1::IndexTuple{N₁}, p2::IndexTuple{N₂}) where {S,N₁,N₂}
    # TODO: check Frobenius-Schur indicators!, and add fermions!
    @boundscheck begin
        all(i->space(tsrc, p1[i]) == space(tdst, i), 1:N₁) ||
            throw(SpaceMismatch("tsrc = $(codomain(tsrc))←$(domain(tsrc)),
            tdst = $(codomain(tdst))←$(domain(tdst)), p1 = $(p1), p2 = $(p2)"))
        all(i->space(tsrc, p2[i]) == space(tdst, N₁+i), 1:N₂) ||
            throw(SpaceMismatch("tsrc = $(codomain(tsrc))←$(domain(tsrc)),
            tdst = $(codomain(tdst))←$(domain(tdst)), p1 = $(p1), p2 = $(p2)"))
    end

    G = sectortype(S)
    if G === Trivial
        cod = codomain(tsrc)
        dom = domain(tsrc)
        n = length(cod)
        pdata = (p1..., p2...)
        axpby!(α, permutedims(tsrc[], pdata), β, tdst[])
    elseif FusionStyle(G) isa Abelian
        K = Threads.nthreads()
        if K > 1
            let iterator = fusiontrees(tsrc)
                Threads.@threads for k = 1:K
                    counter = 0
                    for (f1,f2) in iterator
                        counter += 1
                        if mod1(counter, K) == k
                            _addabelianblock!(α, tsrc, β, tdst, p1, p2, f1, f2)
                        end
                    end
                end
            end
        else # debugging is easier this way
            for (f1,f2) in fusiontrees(tsrc)
                _addabelianblock!(α, tsrc, β, tdst, p1, p2, f1, f2)
            end
        end
    else
        cod = codomain(tsrc)
        dom = domain(tsrc)
        n = length(cod)
        pdata = (p1...,p2...)
        if iszero(β)
            fill!(tdst, β)
        elseif β != 1
            mul!(tdst, β, tdst)
        end
        for (f1,f2) in fusiontrees(tsrc)
            for ((f1′,f2′), coeff) in permute(f1, f2, p1, p2)
                @inbounds for i in p2
                    if i <= n && !isdual(cod[i])
                        b = f1.uncoupled[i]
                        coeff *= frobeniusschur(b) #*fermionparity(b)
                    end
                end
                @inbounds for i in p1
                    if i > n && isdual(dom[i-n])
                        b = f2.uncoupled[i-n]
                        coeff /= frobeniusschur(b) #*fermionparity(b)
                    end
                end
                @inbounds axpy!(α*coeff, permutedims(tsrc[f1,f2], pdata), tdst[f1′,f2′])
            end
        end
    end
    return tdst
end

function _addabelianblock!(α, tsrc::AbstractTensorMap{S},
                            β, tdst::AbstractTensorMap{S,N₁,N₂},
                            p1::IndexTuple{N₁}, p2::IndexTuple{N₂},
                            f1::FusionTree, f2::FusionTree) where {S,N₁,N₂}
    cod = codomain(tsrc)
    dom = domain(tsrc)
    n = length(cod)
    (f1′,f2′), coeff = first(permute(f1, f2, p1, p2))
    @inbounds for i in p2
        if i <= n && !isdual(cod[i])
            b = f1.uncoupled[i]
            coeff *= frobeniusschur(b) #*fermionparity(b)
        end
    end
    @inbounds for i in p1
        if i > n && isdual(dom[i-n])
            b = f2.uncoupled[i-n]
            coeff /= frobeniusschur(b) #*fermionparity(b)
        end
    end
    pdata = (p1...,p2...)
    @inbounds axpby!(α*coeff, permutedims(tsrc[f1,f2], pdata), β, tdst[f1′,f2′])
end

function trace!(α, tsrc::AbstractTensorMap{S}, β, tdst::AbstractTensorMap{S,N₁,N₂},
                p1::IndexTuple{N₁}, p2::IndexTuple{N₂},
                q1::IndexTuple{N₃}, q2::IndexTuple{N₃}) where {S,N₁,N₂,N₃}
    # TODO: check Frobenius-Schur indicators!, and  add fermions!
    @boundscheck begin
        all(i->space(tsrc, p1[i]) == space(tdst, i), 1:N₁) ||
            throw(SpaceMismatch("trace: tsrc = $(codomain(tsrc))←$(domain(tsrc)),
                    tdst = $(codomain(tdst))←$(domain(tdst)), p1 = $(p1), p2 = $(p2)"))
        all(i->space(tsrc, p2[i]) == space(tdst, N₁+i), 1:N₂) ||
            throw(SpaceMismatch("trace: tsrc = $(codomain(tsrc))←$(domain(tsrc)),
                    tdst = $(codomain(tdst))←$(domain(tdst)), p1 = $(p1), p2 = $(p2)"))
        all(i->space(tsrc, q1[i]) == dual(space(tsrc, q2[i])), 1:N₃) ||
            throw(SpaceMismatch("trace: tsrc = $(codomain(tsrc))←$(domain(tsrc)),
                    q1 = $(p1), q2 = $(q2)"))
    end

    G = sectortype(S)
    if G === Trivial
        cod = codomain(tsrc)
        dom = domain(tsrc)
        n = length(cod)
        pdata = (p1..., p2...)
        TO._trace!(α, tsrc[], β, tdst[], pdata, q1, q2)
    # elseif FusionStyle(G) isa Abelian
    # TODO: is it worth multithreading Abelian case for traces?
    else
        cod = codomain(tsrc)
        dom = domain(tsrc)
        n = length(cod)
        pdata = (p1...,p2...)
        if iszero(β)
            fill!(tdst, β)
        elseif β != 1
            mul!(tdst, β, tdst)
        end
        r1 = (p1..., q1...)
        r2 = (p2..., q2...)
        for (f1,f2) in fusiontrees(tsrc)
            for ((f1′,f2′), coeff) in permute(f1, f2, r1, r2)
                f1′′, g1 = split(f1′, StaticLength(N₁))
                f2′′, g2 = split(f2′, StaticLength(N₂))
                if g1 == g2
                    @inbounds for i in r2
                        if i <= n && !isdual(cod[i])
                            b = f1.uncoupled[i]
                            coeff *= frobeniusschur(b) #*fermionparity(b)
                        end
                    end
                    @inbounds for i in r1
                        if i > n && isdual(dom[i-n])
                            b = f2.uncoupled[i-n]
                            coeff /= frobeniusschur(b) #*fermionparity(b)
                        end
                    end
                    coeff *= dim(g1.coupled)/dim(g1.uncoupled[1])
                    TO._trace!(α*coeff, tsrc[f1,f2], true, tdst[f1′′,f2′′], pdata, q1, q2)
                end
            end
        end
    end
    return tdst
end

# TODO: contraction with either A or B a rank (1,1) tensor does not require to
# permute the fusion tree and should therefore be special cased. This will speed
# up MPS algorithms
function contract!(α, A::AbstractTensorMap{S}, B::AbstractTensorMap{S},
                    β, C::AbstractTensorMap{S},
                    oindA::IndexTuple{N₁}, cindA::IndexTuple,
                    oindB::IndexTuple{N₂}, cindB::IndexTuple,
                    p1::IndexTuple, p2::IndexTuple,
                    syms::Union{Nothing, NTuple{3,Symbol}} = nothing) where {S,N₁,N₂}

    # check if it is beneficial to chagne the role of A and B
    hsp = has_shared_permuteind
    ipC = TupleTools.invperm((p1..., p2...))
    oindAinC = TupleTools.getindices(ipC, ntuple(n->n, StaticLength(N₁)))
    oindBinC = TupleTools.getindices(ipC, ntuple(n->n+N₁, StaticLength(N₂)))
    memcost1 = dim(A)*(!hsp(A, oindA, cindA))
                + dim(B)*(!hsp(B, cindB, oindB))
                + dim(C)*(!hsp(C, oindAinC, oindBinC))
    memcost2 = dim(B)*(!hsp(B, oindB, cindB))
                + dim(A)*(!hsp(B, cindA, oindA))
                + dim(C)*(!hsp(C, oindBinC, oindAinC))

    if memcost1 > memcost2
        p1′ = map(n->ifelse(n>N₁, n-N₁, n+N₂), p1)
        p2′ = map(n->ifelse(n>N₁, n-N₁, n+N₂), p2)
        return contract!(α, B, A, β, C, oindB, cindB, oindA, cindA, p1′, p2′, syms)
    end

    if syms === nothing
        A′ = permuteind(A, oindA, cindA)
        B′ = permuteind(B, cindB, oindB)
    else
        A′ = cached_permuteind(syms[1], A, oindA, cindA)
        B′ = cached_permuteind(syms[2], B, cindB, oindB)
    end
    ipC = TupleTools.invperm((p1..., p2...))
    oindAinC = TupleTools.getindices(ipC, ntuple(n->n, StaticLength(N₁)))
    oindBinC = TupleTools.getindices(ipC, ntuple(n->n+N₁, StaticLength(N₂)))
    if has_shared_permuteind(C, oindAinC, oindBinC)
        C′ = permuteind(C, oindAinC, oindBinC)
        mul!(C′, A′, B′, α, β)
    else
        if syms === nothing
            C′ = A′*B′
        else
            p1′ = ntuple(identity, StaticLength(N₁))
            p2′ = N₁ .+ ntuple(identity, StaticLength(N₂))
            TC = eltype(C)
            C′ = TO.cached_similar_from_indices(syms[3], TC, oindA, oindB, p1′, p2′, A, B, :N, :N)
            mul!(C′, A′, B′)
        end
        add!(α, C′, β, C, p1, p2)
    end
    return C
end

# # Compatibility layer for working with the `@tensor` macro from TensorOperations
function TO.checked_similar_from_indices(C, T::Type,
    p1::IndexTuple, p2::IndexTuple, A::AbstractTensorMap, CA::Symbol = :N)
    if CA == :N
        checked_similar_from_indices(C, T, p1, p2, A)
    else
        p1 = adjointtensorindices(A, p1)
        p2 = adjointtensorindices(A, p2)
        checked_similar_from_indices(C, T, p1, p2, adjoint(A))
    end
end

function TO.checked_similar_from_indices(C, T::Type,
    poA::IndexTuple, poB::IndexTuple,
    p1::IndexTuple, p2::IndexTuple,
    A::AbstractTensorMap, B::AbstractTensorMap, CA::Symbol = :N, CB::Symbol = :N)

    if CA == :N && CB == :N
        checked_similar_from_indices(C, T, poA, poB, p1, p2, A, B)
    elseif CA == :C && CB == :N
        poA = adjointtensorindices(A, poA)
        checked_similar_from_indices(C, T, poA, poB, p1, p2, adjoint(A), B)
    elseif CA == :N && CB == :C
        poB = adjointtensorindices(B, poB)
        checked_similar_from_indices(C, T, poA, poB, p1, p2, A, adjoint(B))
    else
        poA = adjointtensorindices(A, poA)
        poB = adjointtensorindices(B, poB)
        checked_similar_from_indices(C, T, poA, poB, p1, p2, adjoint(A), adjoint(B))
    end
end

TO.scalar(t::AbstractTensorMap) = scalar(t)

function TO.add!(α, tsrc::AbstractTensorMap{S}, CA::Symbol, β,
    tdst::AbstractTensorMap{S,N₁,N₂}, p1::IndexTuple, p2::IndexTuple) where {S,N₁,N₂}

    if CA == :N
        p = (p1..., p2...)
        pl = TupleTools.getindices(p, codomainind(tdst))
        pr = TupleTools.getindices(p, domainind(tdst))
        add!(α, tsrc, β, tdst, pl, pr)
    else
        p = adjointtensorindices(tsrc, (p1..., p2...))
        pl = TupleTools.getindices(p, codomainind(tdst))
        pr = TupleTools.getindices(p, domainind(tdst))
        add!(α, adjoint(tsrc), β, tdst, pl, pr)
    end
    return tdst
end

function TO.trace!(α, tsrc::AbstractTensorMap{S}, CA::Symbol, β,
    tdst::AbstractTensorMap{S,N₁,N₂}, p1::IndexTuple, p2::IndexTuple,
    q1::IndexTuple, q2::IndexTuple) where {S,N₁,N₂}

    if CA == :N
        p = (p1..., p2...)
        pl = TupleTools.getindices(p, codomainind(tdst))
        pr = TupleTools.getindices(p, domainind(tdst))
        trace!(α, tsrc, β, tdst, pl, pr, q1, q2)
    else
        p = adjointtensorindices(tsrc, (p1..., p2...))
        pl = TupleTools.getindices(p, codomainind(tdst))
        pr = TupleTools.getindices(p, domainind(tdst))
        q1 = adjointtensorindices(tsrc, q1)
        q2 = adjointtensorindices(tsrc, q2)
        trace!(α, adjoint(tsrc), β, tdst, pl, pr, q1, q2)
    end
    return tdst
end

function TO.contract!(α,
    tA::AbstractTensorMap{S}, CA::Symbol,
    tB::AbstractTensorMap{S}, CB::Symbol,
    β, tC::AbstractTensorMap{S,N₁,N₂},
    oindA::IndexTuple, cindA::IndexTuple,
    oindB::IndexTuple, cindB::IndexTuple,
    p1::IndexTuple, p2::IndexTuple,
    syms::Union{Nothing, NTuple{3,Symbol}} = nothing) where {S,N₁,N₂}

    p = (p1..., p2...)
    pl = ntuple(n->p[n], StaticLength(N₁))
    pr = ntuple(n->p[N₁+n], StaticLength(N₂))
    if CA == :N && CB == :N
        contract!(α, tA, tB, β, tC, oindA, cindA, oindB, cindB, pl, pr, syms)
    elseif CA == :N && CB == :C
        oindB = adjointtensorindices(tB, oindB)
        cindB = adjointtensorindices(tB, cindB)
        contract!(α, tA, tB', β, tC, oindA, cindA, oindB, cindB, pl, pr, syms)
    elseif CA == :C && CB == :N
        oindA = adjointtensorindices(tA, oindA)
        cindA = adjointtensorindices(tA, cindA)
        contract!(α, tA', tB, β, tC, oindA, cindA, oindB, cindB, pl, pr, syms)
    elseif CA == :C && CB == :C
        oindA = adjointtensorindices(tA, oindA)
        cindA = adjointtensorindices(tA, cindA)
        oindB = adjointtensorindices(tB, oindB)
        cindB = adjointtensorindices(tB, cindB)
        contract!(α, tA', tB', β, tC, oindA, cindA, oindB, cindB, pl, pr, syms)
    else
        error("unknown conjugation flags: $CA and $CB")
    end
    return tC
end


text:
# optimized methods for ReshapedArrays
function cartesian(a::CartesianIndex{N}, b::CartesianIndex{N}) where N
    CartesianIndices(UnitRange.(a.I, b.I))
end

# column iteration
@inline function iterate(iter::ContigCartIterator, state=first(iter.columnrange))
    state === CartesianIndex() && return nothing
    isless(last(iter.columnrange), state) && return nothing
    v = iterate(iter.arrayrange, state)
    newstate = v === nothing ? CartesianIndex() : v[2]
    ReshapedIndex(state), newstate
end

function _contiguous_iterator(W::ArrayIndexingWrapper{AA}, ::IndexCartesian) where AA<:ReshapedArray
    fi, li = firstlast(W)
    A = parent(W)
    ax = axes(parent(A))
    f = Base.ind2sub_rs(ax, A.mi, fi)
    l = Base.ind2sub_rs(ax, A.mi, li)
    c = ContigCartIterator(CartesianIndices(inds(parent(A))),
                           cartesian(CartesianIndex(f), CartesianIndex(l)))
end

import Base: ==
function ==(a::C, b::C) where C<:ContigCartIterator
    a.arrayrange == b.arrayrange && a.columnrange == b.columnrange
end

# Branching implementation
# @inline isless(I1::CartesianIndex{N}, I2::CartesianIndex{N}) where N = _isless(I1.I, I2.I)
# @inline function _isless(I1::NTuple{N,Int}, I2::NTuple{N,Int}) where N
#     i1, i2 = I1[N], I2[N]
#     isless(i1, i2) && return true
#     isless(i2, i1) && return false
#     _isless(Base.front(I1), Base.front(I2))
# end
# _isless(::Tuple{}, ::Tuple{}) = false

# Select implementation
# @inline isless(I1::CartesianIndex{N}, I2::CartesianIndex{N}) where N = _isless(0, I1.I, I2.I)
# @inline function _isless(ret, I1::NTuple{N,Int}, I2::NTuple{N,Int}) where N
#     newret = ifelse(ret==0, icmp(I1[N], I2[N]), ret)
#     _isless(newret, Base.front(I1), Base.front(I2))
# end
# _isless(ret, ::Tuple{}, ::Tuple{}) = ifelse(ret==1, true, false)
# icmp(a, b) = ifelse(isless(a,b), 1, ifelse(a==b, 0, -1))


text:
@recipe function f(ddf::DateDataFrame)
    # label --> names(ddf)
    return ddf.timestamp, collect(eachcol(ddf.df))
end

@recipe function f(ddf_vec::Vector{<:DateDataFrame})
    ts_vec = timestamp.(ddf_vec)
    @assert all([ts_vec[1]] .== ts_vec[2:end])
    ts = ddf_vec[1].timestamp
    return ts, [ddf.df for ddf in ddf_vec]
end


text:
#TODO: Fix the cross-tab table
function createquantiles!(df::AbstractDataFrame, source::Symbol; target=Symbol(:q, source),
  groupfields::Union{Vector{Symbol}, Nothing} = [:fisyr])

  if target ∉ names(df)
    df[target] = Vector{Union{Missing,eltype(df[source])}}(undef, size(df,1))
  end

  if isnothing(groupfields) || length(groupfields)==0 #no grouping needed
    F = ecdf(collect(skipmissing(df[source])))

    if isfinite(F(0.0)) #we can shortcut this loop if there are no valid ecdf values
      df[target] .= (f::MFloat64->ismissing(f) ? missing : F(f)).(df[source])
    end
  else #group the data using the provided fields
    for subdf::SubDataFrame ∈ groupby(df, :fisyr)
      F = ecdf(collect(skipmissing(subdf[source])))

      if isfinite(F(0.0)) #we can shortcut this loop if there are no valid ecdf values
        subdf[target] .= (f::MFloat64->ismissing(f) ? missing : F(f)).(subdf[source])
      end
  end
end

  return nothing
end

#creates a series of percentile intervals given the cutoffs
function createcutoffintervals(thresholds::Vector{Float64})
  intervals::Vector{Tuple{Float64,Float64}} = Vector{Tuple{Float64,Float64}}()

  sizehint!(intervals, length(thresholds))

  #make sure we are dealing with quantiles
  @assert minimum(thresholds) > 0.0
  @assert maximum(thresholds) == 1.0

  for i ∈ 1:length(thresholds)
    if i == 1
      push!(intervals, (0.0,thresholds[1]))
    else
      push!(intervals, (thresholds[i-1],thresholds[i]))
    end
  end

  return intervals
end

#holds all of the info to make a cross-tab table
struct CrossTabTable
  crosstabcolumns::Vector{Vector{Union{MFloat64,Int}}}

  rowshortnames::Vector{Symbol}
  #rowstatfields::Vector{Symbol}
  rowindex::Dict

  colshortnames::Vector{Symbol}
  #colstatfields::Vector{Symbol}
  colindex::Dict

  rowfield::Symbol
  colfield::Symbol
  tabfield::Symbol
  suffix::String
end

#a helper function to aggregate quantile columns for crosstab
function crosstabaggregate(target::AbstractVector; aggfunction::Function=median)::MFloat64
  local vcomplete::T where T <: Union{typeof(target), Vector{MFloat64}}

  if Missing <: eltype(target)
    vcomplete =  target[(!ismissing).(target)] #don't trust how missings are handled
  else
    vcomplete=target
  end

  if length(vcomplete) > 0
    return aggfunction(vcomplete)
  else
    return missing
  end
end

#makes fourblocker histograms of returns
function tabbedhistograms(df::AbstractDataFrame;
    outputpath::String = GRAPH_PATH,
    rowfield::Symbol=:qlreturn3yr,
    rowthresholds::Vector{Float64} = CROSSTAB_CUTOFFS,
    colfield = :qpprogramexpense,
    colthresholds::Vector{Float64} = CROSSTAB_CUTOFFS,
    tabfield::Symbol = :qpcontributions,
    suffix::String = "",
    shortnameindex::Dict = SHORT_NAME_INDEX,
    graphsuffix::String = GRAPH_SUFFIX)

  local plotnames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()

  rowcutoffintervals::Vector{Tuple{Float64,Float64}} = createcutoffintervals(rowthresholds)
  colcutoffintervals::Vector{Tuple{Float64,Float64}} = createcutoffintervals(colthresholds)

  graphindex::Dict = Dict{Tuple{Int,Int}, Int}()

  #iterate over the thresholds
  for r ∈ 1:length(rowcutoffintervals)
    for c ∈ 1:length(colcutoffintervals)
      dfsub::SubDataFrame =view(df, ((frow::MFloat64,fcol::MFloat64)->
        (!ismissing(frow)) && (!ismissing(fcol)) &&
        (colcutoffintervals[c][1] < fcol) && (fcol < colcutoffintervals[c][2]) &&
        (rowcutoffintervals[r][1] < frow) && (frow < rowcutoffintervals[r][2])
        ).(df[rowfield], df[colfield]), :)

      #index the location of the graph in the vector
      graphindex[(r,c)] = length(plots) + 1
      push!(plotnames,
        "$(shortnameindex[rowfield]): $(rowcutoffintervals[r][1])-$(rowcutoffintervals[r][2])
        $(shortnameindex[colfield]): $(colcutoffintervals[c][1])-$(colcutoffintervals[c][2])")

      #write the grpah
      push!(plots, plot(dfsub,
        x=tabfield,
        Guide.title(plotnames[end]),
        Guide.xlabel("contribution percentile"), Guide.ylabel("density"),
        Geom.line, Stat.histogram(bincount=50, density=true),
        Coord.cartesian(xmin=0.0, xmax=1.0)))

    end
  end

  finalplot::PlotContainer = title(gridstack(
    [plots[graphindex[(1,1)]] plots[graphindex[(1,2)]]
    plots[graphindex[(2,1)]] plots[graphindex[(2,2)]]]),
    "Normalized Contribution Percentile Histograms")

  outname::String = ("tabhist-$(shortnameindex[rowfield])X$(
    shortnameindex[colfield])BY$(shortnameindex[tabfield])$(suffix)_$(graphsuffix)")

  draw(PDF("$outputpath\\$(outname).pdf", 9inch, 7inch), finalplot)

  return nothing
end

#makes a histogram of betas
function plotβtab!(df::DataFrame;
  tabfield::Symbol = :b1lRet3yrCon,#:b1lRet3yrCon,
  suffix::String = "",
  outputpath::String = GRAPH_PATH,
  thresholds::Vector{Float64} = CROSSTAB_CUTOFFS,
  groupfield = :pprogramexpense_avg,
  shortnameindex::Dict = SHORT_NAME_INDEX,
  makequantiles::Bool = true,
  graphsuffix::String = GRAPH_SUFFIX
  )

  local plotnames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()

  cutoffintervals::Vector{Tuple{Float64,Float64}} = createcutoffintervals(thresholds)

  #make quintiles and classifers
  df = deepcopy(df[completecases(df[[groupfield, tabfield]]),:])
  makequantiles && createquantiles!(df, groupfield, groupfields = nothing)
  makequantiles && createquantiles!(df, tabfield, groupfields = nothing)
  tabfield = Symbol(:q, tabfield)
  groupfield = Symbol(:q, groupfield)
  df[:classification] = Vector{MSymbol}(undef, size(df,1))
  intervalnames = [:lowexp, :highexp]

  #given how much of a hack the below is...
  @assert length(cutoffintervals) == 2
  df[:classification] .= (f::Float64 ->
    f < thresholds[1] ? intervalnames[1] : intervalnames[2]).(df[groupfield])
  #=for c ∈ 1:length(cutoffintervals)
    dfsub::SubDataFrame =view(df, (f::Float64->
      (cutoffintervals[c][1] < f) && (f < cutoffintervals[c][2])
      ).(df[groupfield]), :)

    #add the classifier
    dfsub[:classification] .= intervalnames[c]

    #index the location of the graph in the vector
    #=push!(plotnames,
      "$(shortnameindex[groupfield]): $(cutoffintervals[c][1])-$(cutoffintervals[c][2])")

    #write the grpah
    push!(plots, plot(dfsub,
      x=tabfield,
      Guide.title(plotnames[end]),
      Guide.xlabel("contribution beta percentile"), Guide.ylabel("density"),
      Geom.line, Stat.histogram(bincount=25, density=false),
      Coord.cartesian(xmin=0.0, xmax=1.0)))=#

  end=#

  #=finalplot::PlotContainer = title(vstack(
    [plots[1],plots[2]]), "Normalized Contribution Percentile Histograms")=#

  finalplot = plot(df,
    x=tabfield, color=:classification,
    Guide.title("Histogram of beta percentile by governance category"),
    Guide.xlabel("contribution beta percentile"),
    Guide.ylabel("count"),
    Geom.line, Stat.histogram(bincount=25, density=true),
    Coord.cartesian(xmin=0.0, xmax=1.0))

  outname::String = (
    "tab-$(shortnameindex[groupfield])BY$(shortnameindex[tabfield])$(suffix)_$(graphsuffix)")

  draw(PDF("$outputpath\\$(outname).pdf", 9inch, 7inch), finalplot)

  return nothing
end

#creates the cross-tabulation table
function CrossTabTable(df::DataFrame;rowfield::Symbol=:qlreturn3yr,
  rowthresholds::Vector{Float64} = CROSSTAB_CUTOFFS,
  colfield = :qpprogramexpense,
  colthresholds::Vector{Float64} = CROSSTAB_CUTOFFS,
  yearrange::UnitRange = CROSSTAB_YEAR_RANGE,
  groupfield = :fisyr,
  tabfield::Symbol = :qpcontributions,
  aggfunction::Function = crosstabaggregate,
  suffix::String = "",
  shortnameindex::Dict = SHORT_NAME_INDEX)

  #build up indices and naming conventions
  local rowthresholdfields::Vector{Symbol} = (f::Float64 ->
    Symbol(rowfield, "p",Int(round(f*100)))).(rowthresholds)
  #local rowstatfields::Vector{Symbol} = [rowthresholdfields; :total; :N]
  local rowshortnames = (f::Float64 ->
    Symbol(shortnameindex[rowfield], "p",Int(round(f*100)))).(rowthresholds)
  rowshortnames = [rowshortnames; :total; :N]
  local nrows = length(rowshortnames)
  local rowindex::Dict = Dict(rowshortnames[r]=>r for r ∈ 1:length(rowshortnames))

  local colthresholdfields::Vector{Symbol} = (f::Float64 ->
    Symbol(colfield, "p",Int(round(f*100)))).(colthresholds)
  #local colstatfields::Vector{Symbol} = [colthresholdfields; :total; :N]
  local colshortnames = (f::Float64 ->
    Symbol(shortnameindex[colfield], "p",Int(round(f*100)))).(colthresholds)
  colshortnames = [colshortnames; :total; :N]
  local ncols::Int = length(colshortnames)
  local colindex::Dict = Dict(colshortnames[r]=>r for r ∈ 1:length(colshortnames))

  #this is what will actually hold the results
  crosstabcolumns::Vector{Vector{Union{MFloat64,Int}}} =
      (i->Vector{Union{MFloat64,Int}}(undef, nrows)).(1:ncols)

  #create a dataframe to work with
  dfcross::DataFrame = df[(y::MInt -> y ∈ yearrange).(
    df[:fisyr]), [rowfield, colfield, tabfield, groupfield]]

  #draw the graphs
  tabbedhistograms(dfcross,
    rowfield=rowfield,
    rowthresholds=rowthresholds,
    colfield=colfield,
    colthresholds=colthresholds,
    tabfield=tabfield,
    suffix=suffix,
    shortnameindex=shortnameindex)

  dfcross = dfcross[completecases(dfcross),:]
  #CSV.write("$WORKING_PATH\\crossdftest.csv", dfcross)

  #holds the top and bottom cutoff for each interval
  rowcutoffintervals::Vector{Tuple{Float64,Float64}} = createcutoffintervals(rowthresholds)
  colcutoffintervals::Vector{Tuple{Float64,Float64}} = createcutoffintervals(colthresholds)

  #now create the main table (We will need to make the row totals later)
  #NOTE: assumes that the cross tab columns come first!
  for c ∈ 1:length(colcutoffintervals)
    dfsub::SubDataFrame =view(dfcross, (f::MFloat64-> (!ismissing(f)) &&
      (colcutoffintervals[c][1] < f) && (f < colcutoffintervals[c][2])).(dfcross[colfield]), :)
    crosstabcolumns[c][rowindex[:total]] = aggfunction(dfsub[tabfield])
    crosstabcolumns[c][rowindex[:N]] = size(dfsub,1)

    for r ∈ 1:length(rowcutoffintervals)
      dfssub::SubDataFrame =view(dfsub, (f::MFloat64-> (!ismissing(f)) &&
        (rowcutoffintervals[r][1] < f) && (f < rowcutoffintervals[r][2])).(dfsub[rowfield]), :)
      crosstabcolumns[c][r] = aggfunction(dfssub[tabfield])
    end
  end

  #now make the rowwise subtotals
  for  r ∈ 1:length(rowcutoffintervals)
    dfsub::SubDataFrame =view(dfcross, (f::MFloat64-> (!ismissing(f)) &&
      (rowcutoffintervals[r][1] < f) && (f < rowcutoffintervals[r][2])).(dfcross[rowfield]), :)
    crosstabcolumns[colindex[:total]][r] = aggfunction(dfsub[tabfield])
    crosstabcolumns[colindex[:N]][r] = size(dfsub,1)
  end

  println(crosstabcolumns)

  #now make the grand total
  crosstabcolumns[colindex[:total]][rowindex[:total]] = aggfunction(dfcross[tabfield])
  crosstabcolumns[colindex[:N]][rowindex[:N]] = size(dfcross,1)

  return CrossTabTable(
    crosstabcolumns,

    rowshortnames,
    rowindex,

    colshortnames,
    colindex,

    rowfield,
    colfield,
    tabfield,
    suffix)
end

function writecrosstabtable(tab::CrossTabTable,
  outputpath::String = TABLE_PATH,
  decimals::Int = DECIMALS_RETURN,
  shortnameindex::Dict = SHORT_NAME_INDEX,
  graphsuffix::String = GRAPH_SUFFIX
  )

  local ncols::Int = length(tab.crosstabcolumns) #number of columns excluding the year
  local nrows::Int = length(tab.crosstabcolumns[1])

  #collect the keys reverse the dicitonaries for notational convenience
  local colindexkeys::Vector{Symbol} = collect(keys(tab.colindex))
  local rowindexkeys::Vector{Symbol} = collect(keys(tab.rowindex))

  #local revcolindex::Dict = Dict(c=>tab.colindex[colindexkeys[c]]  for c ∈ 1:ncols)
  #local revrowindex::Dict = Dict(r=>tab.rowindex[rowindexkeys[r]]  for r ∈ 1:nrows)

  #make the column names and sizes
  local colNames::Vector{Vector{String}} = Vector{Vector{String}}()
  push!(colNames, (s->string(s)).(tab.colshortnames))

  local widthColNames::Vector{Vector{Int}} = Vector{Vector{Int}}()
  push!(widthColNames, ones(Int, ncols))

  #make the row names
  local descRowNames::Vector{String} = (s->string(s)).(tab.rowshortnames)

  #need to print the descriptive rows
  descContent::Vector{Vector{String}} = (
    (i::Int)->Vector{String}(undef, ncols)).(1:nrows)

  #NOTE: tab has the data in COLUMNS while descContet absorbs data in ROWS (flipped indices)
  for r::Int ∈ 1:nrows
    for c::Int ∈ 1:ncols

      #Ns are ints
      if tab.rowshortnames[r] == :N || tab.colshortnames[c] == :N
        #lookup the collumn number from the column index for the pth column symbol
        descContent[r][c] = num2Str(tab.crosstabcolumns[c][r], Ints=true)
      else
        descContent[r][c] = num2Str(tab.crosstabcolumns[c][r], decimals)
      end
    end
  end

  outtabletitle::String = ("crosstab-$(shortnameindex[tab.rowfield])X$(
    shortnameindex[tab.colfield])BY$(shortnameindex[tab.tabfield])$(tab.suffix)")
  outtable::String = texTable(titleCaption=outtabletitle,
      caption="""See Tex File""", #caption
      colNames= colNames, #colNames
      descRowNames=descRowNames, #descRowNames
      descContent = descContent, #descContent
      nakedTable=true,
      widthColNames=widthColNames
    )

  writeNakedTable(outtable, path=outputpath, outName="$(outtabletitle)_$(graphsuffix).tex")

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using RecordedArrays: return_type

@testset "$E" for E in (
    StaticEntry,
    StaticEntry{Int},
    StaticEntry{Int,Float64},
    DynamicEntry,
    DynamicEntry{Int},
    DynamicEntry{Int,Float64},
)
    c = DiscreteClock(3)
    e1 = E(1, c)
    e2 = E(0x2, currenttime(c))
    @test e1 isa return_type(E, 1, c)
    @test e1 isa return_type(E, Int, c)
    @test e2 isa return_type(E, 0x2, c)
    @test e2 isa return_type(E, typeof(0x2), c)

    for t in c
        if E <: DynamicEntry
            store!(e1, t, c)
            t != 3 && store!(e2, t + 1, c)
        end
        t == 2 && del!(e1, c)
        t == 3 && del!(e2, c)
    end

    @test e1 == deepcopy(e1)
    @test e2 == deepcopy(e2)
    @test e1 != e2

    @test gettime(e1, -1) == 0
    @test gettime(e2, -1) == 0
    if E <: StaticEntry
        @test getts(e1) == [0, 2]
        @test getts(e2) == [0, 3]
        @test tspan(e1) == 2
        @test tspan(e2) == 3
        @test getvs(e1) == [1, 1]
        @test getvs(e2) == [2, 2]
        @test gettime(e1, 0) == 1
        @test gettime(e2, 0) == 2
        @test gettime(e1, 1) == 1
        @test gettime(e2, 1) == 2
        @test gettime(e1, 2) == 1
        @test gettime(e2, 2) == 2
        @test gettime(e1, 3) == 0
        @test gettime(e2, 3) == 2
        @test gettime(e1, 1:2) == [1, 1]
        @test gettime(e2, 1:2) == [2, 2]
        @test gettime(e1, 1:3) == [1, 1, 0]
        @test gettime(e2, 1:3) == [2, 2, 2]
        @testset "gettime($alg)" for alg in (LinearSearch(), BinarySearch())
            @test gettime(alg, e1, 0) == 1
            @test gettime(alg, e2, 0) == 2
            @test gettime(alg, e1, 1) == 1
            @test gettime(alg, e2, 1) == 2
            @test gettime(alg, e1, 2) == 1
            @test gettime(alg, e2, 2) == 2
            @test gettime(alg, e1, 3) == 0
            @test gettime(alg, e2, 3) == 2
            @test gettime(alg, e1, 1:2) == [1, 1]
            @test gettime(alg, e2, 1:2) == [2, 2]
            @test gettime(alg, e1, 2:3) == [1, 0]
            @test gettime(alg, e2, 2:3) == [2, 2]
            @test gettime(alg, e1, 1:3) == [1, 1, 0]
            @test gettime(alg, e2, 1:3) == [2, 2, 2]
            @test gettime(alg, e1, -1:3) == [0, 1, 1, 1, 0]
            @test gettime(alg, e2, -1:3) == [0, 2, 2, 2, 2]
        end
    elseif E <: DynamicEntry
        @test getts(e1) == 0:3
        @test getts(e2) == 0:2
        @test tspan(e1) == 3
        @test tspan(e2) == 2
        @test getvs(e1) == [1, 1, 2, 3]
        @test getvs(e2) == [2, 2, 3]
        @test gettime(e1, 0) == 1
        @test gettime(e2, 0) == 2
        @test gettime(e1, 1) == 1
        @test gettime(e2, 1) == 2
        @test gettime(e1, 2) == 2
        @test gettime(e2, 2) == 3
        @test gettime(e1, 3) == 3
        @test gettime(e2, 3) == 3
        @test gettime(e1, 1:2) == [1, 2]
        @test gettime(e2, 1:2) == [2, 3]
        @test gettime(e1, 1:3) == [1, 2, 3]
        @test gettime(e2, 1:3) == [2, 3, 3]
        @test gettime(e1, 1.5:2.5) == [1, 2]
        @test gettime(e2, 1.5:2.5) == [2, 3]
        @test gettime(e1, 1.5:3.5) == [1, 2, 3]
        @test gettime(e2, 1.5:3.5) == [2, 3, 3]
        @testset "gettime($alg)" for alg in (LinearSearch(), BinarySearch())
            @test gettime(alg, e1, 0) == 1
            @test gettime(alg, e2, 0) == 2
            @test gettime(alg, e1, 1) == 1
            @test gettime(alg, e2, 1) == 2
            @test gettime(alg, e1, 2) == 2
            @test gettime(alg, e2, 2) == 3
            @test gettime(alg, e1, 3) == 3
            @test gettime(alg, e2, 3) == 3
            @test gettime(alg, e1, 1:2) == [1, 2]
            @test gettime(alg, e2, 1:2) == [2, 3]
            @test gettime(alg, e1, 1:3) == [1, 2, 3]
            @test gettime(alg, e2, 1:3) == [2, 3, 3]
            @test gettime(alg, e1, 1.5:2.5) == [1, 2]
            @test gettime(alg, e2, 1.5:2.5) == [2, 3]
            @test gettime(alg, e1, 1.5:3.5) == [1, 2, 3]
            @test gettime(alg, e2, 1.5:3.5) == [2, 3, 3]
            @test gettime(alg, e1, -1:3) == [0, 1, 1, 2, 3]
            @test gettime(alg, e2, -1:3) == [0, 2, 2, 3, 3]
        end
    end
end


text:
module Atomistic

using AtomsBase
using InteratomicPotentials
using DFTK
using NBodySimulator
using PeriodicTable
using Plots
using PyCall
using Random
using StaticArrays
using Unitful
using UnitfulAtomic
using UnitfulRecipes

import Base: @kwdef, Fix2
import InteratomicPotentials: ArbitraryPotential
import DFTK: Mixing
import NBodySimulator: Body, BoundaryConditions, NullThermostat, SimulationResult, Thermostat
import Plots: Plot

include("unit_convention.jl")

# API
include("api/exceptions.jl")
export MolecularDynamicsResult, get_system, get_time_range, temperature, reference_temperature, kinetic_energy, potential_energy, total_energy, rdf
include("api/md_result.jl")
export MolecularDynamicsSimulator, simulate
include("api/md_simulator.jl")

# Integrations
export DynamicSystem
include("integrations/atomsbase_integration.jl")
export DFTKPotential, analyze_convergence
include("integrations/dftk_integration.jl")

# Implementations
export NBSimulator, NBSResult
include("implementations/nbodysimulator.jl")

# Analysis
export plot_temperature, plot_temperature!, plot_energy, plot_energy!, plot_rdf
include("analysis/plotting.jl")
export write_nbs_animation, write_ase_trajectory
include("analysis/visualization.jl")

end


text:
#///////////////////////////////////////
#// File Name: pets_test.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2020/11/06
#// Description: Test code for src/pets.jl
#///////////////////////////////////////

using LinearAlgebra
using Statistics
using Random

@testset "PETS Test" begin

    f_stochastic(x, u, rng, use_true_model=false) = x + u + rand(rng, length(x))
    c(k, x, u) = sum(abs.(u));
    h(x) = 1.0;
    N = 20;

    problem = FiniteHorizonGenerativeOptimalControlProblem(f_stochastic, c, h, N);

    # CrossEntropyDirectOptimizationSolver test
    μ_init_array = [zeros(2) for ii = 1 : N];
    Σ_init_array = [Matrix(1.0I, 2, 2) for ii = 1 : N];
    direct_solver = CrossEntropyDirectOptimizationSolver(μ_init_array, Σ_init_array,
                                                         num_control_samples=20, num_trajectory_samples=100,
                                                         num_elite=5, iter_max=20, smoothing_factor=0.1);

    @test direct_solver.N == N
    @test direct_solver.iter_current == 0;
    @test direct_solver.μ_array == μ_init_array
    @test direct_solver.Σ_array == Σ_init_array

    direct_solver.iter_current = 10;
    direct_solver.μ_array = [ones(2) for ii = 1 : N];
    direct_solver.Σ_array = [Matrix(0.1I, 2, 2) for ii = 1 : N];
    initialize!(direct_solver)

    @test direct_solver.iter_current == 0
    @test direct_solver.μ_array == μ_init_array
    @test direct_solver.Σ_array == Σ_init_array

    # compute_cost_serial test
    rng = MersenneTwister(1234);
    control_sequence_array = [[rand(rng, 2) for tt = 1 : N] for ii = 1 : direct_solver.num_control_samples];
    x_init = zeros(2);
    cost_array_serial = compute_cost_serial(direct_solver, problem, x_init, control_sequence_array, rng);
    # compute_cost test
    cost_array = compute_cost(direct_solver, problem, x_init, control_sequence_array, rng);
    @test all(cost_array .== cost_array_serial) # note that this holds because c is not dependent on x. Otherwise they'd be different.

    rng = MersenneTwister(1234);
    @test length(cost_array) == direct_solver.num_control_samples
    for ii = 1 : direct_solver.num_control_samples
        x = x_init;
        cost = 0.0;
        for tt = 1 : direct_solver.N
            cost += problem.c(tt - 1, x, control_sequence_array[ii][tt])
            x = problem.f_stochastic(x, control_sequence_array[ii][tt], rng)
        end
        cost += problem.h(x)
        @test cost ≈ cost_array[ii]
    end

    # get_elite_samples test
    control_sequence_elite_array = get_elite_samples(direct_solver, control_sequence_array, cost_array);

    @test length(control_sequence_elite_array) == direct_solver.num_elite
    elite_idx_array = [pair[1] for pair in sort(collect(zip(1:1:direct_solver.num_control_samples, cost_array)), by=x->x[2])[1 : direct_solver.num_elite]]
    @test control_sequence_elite_array == control_sequence_array[elite_idx_array]

    # compute_new_distribution test
    μ_new_array, Σ_new_array = compute_new_distribution(direct_solver, control_sequence_elite_array);

    @test length(μ_new_array) == direct_solver.N;
    @test all([length(μ) == length(x_init) for μ in μ_new_array])
    @test length(Σ_new_array) == direct_solver.N;
    @test all([size(Σ) == (length(x_init), length(x_init)) for Σ in Σ_new_array])
    for tt = 1 : direct_solver.N
        @test μ_new_array[tt] ≈ (1.0 - direct_solver.smoothing_factor).*mean([control_sequence_elite_array[ii][tt] for ii = 1 : direct_solver.num_elite]) +
                                direct_solver.smoothing_factor.*direct_solver.μ_array[tt]
        @test Σ_new_array[tt] ≈ (1.0 - direct_solver.smoothing_factor).*Diagonal(var([control_sequence_elite_array[ii][tt] for ii = 1 : direct_solver.num_elite])) +
                                direct_solver.smoothing_factor.*direct_solver.Σ_array[tt]
    end

    # step! test
    rng = MersenneTwister(1234);
    step!(direct_solver, problem, x_init, rng, false, false);
    @test direct_solver.iter_current == 1;

    # solve! test
    control_array, ~ = solve!(direct_solver, problem, x_init, rng, use_true_model=false, verbose=false);

    @test direct_solver.iter_current == direct_solver.iter_max

end


text:
using Gtk, Graphics
using ALifeBenchmark
import Random
using Profile
using Serialization

model = GebModel(size=20)

c = @GtkCanvas()
win = GtkWindow(c, "Geb", 1000, 800)
showall(win)

function paint(ctx, h, w, model::GebModel)
    set_source_rgb(ctx, 0.5, 0.5, 0.5)

    for i in 1:(model.size-1)
        move_to(ctx, 0, h/model.size * i)
        line_to(ctx, w, h/model.size * i)
    end
    for i in 1:(model.size-1)
        move_to(ctx, w/model.size * i, 0)
        line_to(ctx, w/model.size * i, h)
    end
    stroke(ctx)

    for organism in model.organisms
        paint(ctx, h, w, organism)
    end
end

function paint(ctx, h, w, organism::GebOrganism)
    x, y = organism.coordinates

    x, y = x / model.size, y / model.size
    x, y = x*w, y*h

    set_source_rgb(ctx, 0, 0, min(1.0, organism.age / 10))

    circle(ctx, x, y, 3)
    fill(ctx)

    move_to(ctx, x, y)
    line_to(ctx, x + 10*cosd(organism.direction), y + 10*sind(organism.direction))
    stroke(ctx)
end

function paint()
    @guarded draw(c) do widget
        ctx = getgc(c)
        h = height(c)
        w = width(c)

        rectangle(ctx, 0, 0, w, h)
        set_source_rgb(ctx, 1, 1, 1)
        fill(ctx)

        paint(ctx, h, w, model)
    end
    show(c)
end

function run(m, n)
    Random.seed!(105)

    n_epochs = 0

    for _ in 1:m
        @time for _ in 1:n
            execute!(model)
            n_epochs += 1

            if n_epochs % 500 == 0
                println(n_epochs)
            end
        end

        println(model)

        for _ in 1:200
            execute!(model)
            paint()
            sleep(0.1)
            n_epochs += 1
        end
    end
end

run(100, 0)

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# This file constains testset for Writer

@testset "WriterTestSet" begin
    @info "Running WriterTestSet ..."

    # Preliminaries
    randdirname() = string(uuid4())
    randfilename() = join([string(uuid4()), ".jld2"], "")
    testdir = tempdir()

    # Writer construction
    writer = Writer(input=Inport(3), buflen=10)
    writer = Writer(input=Inport(3), buflen=10, path=joinpath(testdir, randfilename()))
    path = joinpath(testdir, randfilename())
    writer = Writer(input=Inport(3), path=path)
    @test typeof(writer.trigger) == Inpin{Float64}
    @test typeof(writer.handshake) == Outpin{Bool}
    @test isa(writer.input, Inport)
    @test length(writer.input) == 3
    @test size(writer.timebuf) == (64,)
    @test size(writer.databuf) == (3, 64)
    @test writer.plugin === nothing
    @test writer.callbacks === nothing
    @test typeof(writer.sinkcallback) <: Callback

    # Reading and writing into Writer
    writer = Writer(input=Inport())
    open(writer)
    for t in 1 : 10
        write!(writer, t, sin(t))
    end
    close(writer)
    data = read(writer, flatten=false)
    for (t, u) in data
        @test sin(t) == u
    end
    data = fread(writer.file.path, flatten=false)
    for (t, u) in data
        @test sin(t) == u
    end

    # Moving/Copying Writer file.
    filename = randfilename()
    dirnames = map(1:3) do i
        path = joinpath(testdir, randdirname())
        ispath(path) || mkdir(path)
        path
    end
    paths = map(dname -> joinpath(dname, filename), dirnames)
    w = Writer(input=Inport(), path=paths[1])
    mv(w, dirnames[2], force=true)
    @test w.file.path == paths[2]
    cp(w, dirnames[3], force=true)
    @test isfile(paths[3])

    # Driving Writer
    writer = Writer(input=Inport(3), buflen=10)
    open(writer)
    oport, iport, trg, hnd, comptask, outtask = prepare(writer)
    for t in 1 : 100
        put!(trg, t)
        put!(oport, ones(3)*t)
        take!(hnd)
        @test read(writer.timebuf) == t
        @test [read(pin.links[1].buffer) for pin in oport] == ones(3) * t
    end
    close(writer)
    t, x = read(writer, flatten=true)
    @test t == collect(1 : 100)
    @test x == [collect(1:100) collect(1:100) collect(1:100)]
    put!(trg, NaN)
    sleep(0.1)
    @test istaskdone(comptask)

    @info "Done WriterTestSet."
end  # testset


text:
using StrTables, HTML_Entities

@static VERSION < v"0.7.0-DEV" ? (using Base.Test) : (using Test)

# Test the functions lookupname, matches, longestmatches, completions
# Check that characters from all 3 tables (BMP, non-BMP, string) are tested

const def = HTML_Entities.default

@testset "HTML_Entities" begin
@testset "lookupname" begin
    @test lookupname(def, SubString("My name is Spock", 12)) == ""
    @test lookupname(def, "foobar") == ""
    @test lookupname(def, "nle")    == "\u2270"
    @test lookupname(def, "Pscr")   == "\U1d4ab"
    @test lookupname(def, "lvnE")   == "\u2268\ufe00"
end

@testset "matches" begin
    @test isempty(matches(def, ""))
    @test isempty(matches(def, "\u201f"))
    @test isempty(matches(def, SubString("This is \u201f", 9)))
    for (chrs, exp) in (("\u2270", ["nle", "nleq"]),
                        ("\U1d4ab", ["Pscr"]),
                        ("\U1d51e", ["afr"]),
                        ("\u2268\ufe00", ["lvertneqq", "lvnE"]))
        res = matches(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end

@testset "longestmatches" begin
    @test isempty(longestmatches(def, "\u201f abcd"))
    @test isempty(longestmatches(def, SubString("This is \U201f abcd", 9)))
    for (chrs, exp) in (("\u2270 abcd", ["nle", "nleq"]),
                        ("\U1d4ab abcd", ["Pscr"]),
                        ("\u2268\ufe00 silly", ["lvertneqq", "lvnE"]))
        res = longestmatches(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end

@testset "completions" begin
    @test isempty(completions(def, "ScottPaulJones"))
    @test isempty(completions(def, SubString("My name is Scott", 12)))
    for (chrs, exp) in (("and", ["and", "andand", "andd", "andslope", "andv"]),
                        ("um", ["umacr", "uml"]))
        res = completions(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end
end