text:
"""
    function flatten_chain(chn::Chains)

Provide a vector of all MCMC draws as `NamedTuple` objects.
"""
function flatten_chain(chn::Chains)
    param_names = names(chn)
    [get(chn[i, :, k], param_names; flatten = true) for i = 1:size(chn, 1), k = 1:size(chn, 3)][:]
end

"""
    function flatten_chain(chn::Chains,n::Integer)

Provide a vector of `n` sampled MCMC draws as `NamedTuple` objects.
"""
function flatten_chain(chn::Chains, n::Integer)
    sampled_chn = sample(chn, n)
    param_names = names(chn)
    [get(sampled_chn[i, :, k], param_names; flatten = true) for i = 1:size(sampled_chn, 1), k = 1:size(sampled_chn, 3)][:]
end

"""
    function flatten_chain(model::KenyaCoVaccines.CoVAreaModel)

Provide a vector of all MCMC draws as `NamedTuple` objects from a `KenyaCoVaccines.CoVAreaModel` object
"""
function flatten_chain(model::KenyaCoVaccines.CoVAreaModel)
    flatten_chain(model.MCMC_results.chain)
end

"""
    function flatten_chain(model::KenyaCoVaccines.CoVAreaModel,n::Integer)

Provide a vector of `n` sampled MCMC draws as `NamedTuple` objects from a `KenyaCoVaccines.CoVAreaModel` object
"""
function flatten_chain(model::KenyaCoVaccines.CoVAreaModel, n::Integer)
    flatten_chain(model.MCMC_results.chain, n)
end

"""
    project_daily_infections_all(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date; ω = 1 / 180)

Return total infections as a timeseries, this sums over age and vaccine status. Simulation use the parameters in `θ` with waning immunity rate `ω`.
"""
function project_daily_infections_all(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date; ω = 1 / 180)
    @unpack PCR_cases, sero_cases, init_sero, baseline_sero_array, PCR_array, sero_sensitivity, sero_specificity, p_symp, p_sus, hosp_rate_by_age, N, M_BB, prob, α, αP, γA, γM, relative_testing_rate, alpha_variant_time, delta_variant_time, model_start_time, model_end_time, janstarttime, VE_acquisition, VE_infectiousness, VE_severe_disease_risk, M_county_ho, M_county_school, M_county_work, M_county_other = model
    @unpack β₀, β_home, β_school, β_other, β_work, ϵ, χ, p_test, E₀, inc_R_αvar, time_scale_αvar, mid_point_αvar, inc_R_δvar, time_scale_δvar, mid_point_δvar, init_scale = θ  #
    # Vaccine efficacies
    ve_ac1 = VE_acquisition[1]
    ve_ac2 = VE_acquisition[2]
    ve_ac3 = VE_acquisition[3]
    ve_inf1 = VE_infectiousness[1]
    ve_inf2 = VE_infectiousness[2]
    ve_inf3 = VE_infectiousness[3]
    ve_dis1 = VE_severe_disease_risk[1]
    ve_dis2 = VE_severe_disease_risk[2]
    ve_dis3 = VE_severe_disease_risk[3]
    #Set transmission parameters --- SCHOOLS START SHUT IN THIS MODEL
    p = [β₀, β_home, 0.0, β_other, β_work, α, ϵ, αP, γA, γM, ω, inc_R_αvar, time_scale_αvar, alpha_variant_time + mid_point_αvar * 30, inc_R_δvar, time_scale_δvar, delta_variant_time + mid_point_δvar * 30, init_scale, ve_ac1, ve_ac2, ve_ac3, ve_inf1, ve_inf2, ve_inf3, ve_dis1, ve_dis2, ve_dis3]
    u0 = KenyaCoVaccines.create_initial_conditions(θ, model)
    #Set up the callback for lockdown
    schoolsclosure_times = [Float64((Date(2021, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 12, 23) - Date(2020, 12, 1)).value)]
    function affect_schoolsclosed!(integrator)
        integrator.p[3] = 0.0 #Schools shut
    end
    schoolsclosed_cb = PresetTimeCallback(schoolsclosure_times, affect_schoolsclosed!, save_positions = (false, false))

    #Set up callback for reopening Schools
    schoolsreopening_times = [Float64((Date(2021, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 11) - Date(2020, 12, 1)).value)]
    function affect_schoolsopening!(integrator)
        integrator.p[3] = β_school # Schools open #
    end
    schoolsreopening_cb = PresetTimeCallback(schoolsreopening_times, affect_schoolsopening!, save_positions = (false, false))

    cb = CallbackSet(schoolsclosed_cb, schoolsreopening_cb)

    #Solve for daily incidence by age (not currently broken down by severity of disease)
    sim_duration = (projection_date - Date(2020, 12, 1)).value
    sol = solve(prob, BS3(); abstol = 1e-3,
        tspan = (0, sim_duration),
        callback = cb,
        u0 = u0, p = p,
        saveat = 1,
        # isoutofdomain = (u, p, t) -> any(x -> x < 0, u),
        verbose = true)
    ι = clamp!(KenyaCoVaccines.get_incidence_time_array(sol, 10), 0.0, Inf)  # Sum over age and vaccine doses
    return ι
end

"""
    project_daily_infections_by_age_and_vaccine(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date; ω = 1 / 180)

Return total infections as a timeseries, this sums over age and vaccine status. Simulation use the parameters in `θ` with waning immunity rate `ω`.
"""
function project_daily_infections_by_age_and_vaccine(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date; ω = 1 / 180)
    @unpack PCR_cases, sero_cases, init_sero, baseline_sero_array, PCR_array, sero_sensitivity, sero_specificity, p_symp, p_sus, hosp_rate_by_age, N, M_BB, prob, α, αP, γA, γM, relative_testing_rate, alpha_variant_time, delta_variant_time, model_start_time, model_end_time, janstarttime, VE_acquisition, VE_infectiousness, VE_severe_disease_risk, M_county_ho, M_county_school, M_county_work, M_county_other = model
    @unpack β₀, β_home, β_school, β_other, β_work, ϵ, χ, p_test, E₀, inc_R_αvar, time_scale_αvar, mid_point_αvar, inc_R_δvar, time_scale_δvar, mid_point_δvar, init_scale = θ  #
    # Vaccine efficacies
    ve_ac1 = VE_acquisition[1]
    ve_ac2 = VE_acquisition[2]
    ve_ac3 = VE_acquisition[3]
    ve_inf1 = VE_infectiousness[1]
    ve_inf2 = VE_infectiousness[2]
    ve_inf3 = VE_infectiousness[3]
    ve_dis1 = VE_severe_disease_risk[1]
    ve_dis2 = VE_severe_disease_risk[2]
    ve_dis3 = VE_severe_disease_risk[3]
    #Set transmission parameters --- SCHOOLS START SHUT IN THIS MODEL
    p = [β₀, β_home, 0.0, β_other, β_work, α, ϵ, αP, γA, γM, ω, inc_R_αvar, time_scale_αvar, alpha_variant_time + mid_point_αvar * 30, inc_R_δvar, time_scale_δvar, delta_variant_time + mid_point_δvar * 30, init_scale, ve_ac1, ve_ac2, ve_ac3, ve_inf1, ve_inf2, ve_inf3, ve_dis1, ve_dis2, ve_dis3]
    u0 = KenyaCoVaccines.create_initial_conditions(θ, model)
    #Set up the callback for lockdown
    schoolsclosure_times = [Float64((Date(2021, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 12, 23) - Date(2020, 12, 1)).value)]
    function affect_schoolsclosed!(integrator)
        integrator.p[3] = 0.0 #Schools shut
    end
    schoolsclosed_cb = PresetTimeCallback(schoolsclosure_times, affect_schoolsclosed!, save_positions = (false, false))

    #Set up callback for reopening Schools
    schoolsreopening_times = [Float64((Date(2021, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 11) - Date(2020, 12, 1)).value)]
    function affect_schoolsopening!(integrator)
        integrator.p[3] = β_school # Schools open #
    end
    schoolsreopening_cb = PresetTimeCallback(schoolsreopening_times, affect_schoolsopening!, save_positions = (false, false))

    cb = CallbackSet(schoolsclosed_cb, schoolsreopening_cb)

    #Solve for daily incidence by age (not currently broken down by severity of disease)
    sim_duration = (projection_date - Date(2020, 12, 1)).value
    sol = solve(prob, AutoTsit5(Rosenbrock23()); abstol = 1e-3,
        tspan = (0, sim_duration),
        callback = cb,
        u0 = u0, p = p,
        saveat = 1,
        # isoutofdomain = (u, p, t) -> any(x -> x < -1e-6, u),
        verbose = true)

    ι_first = clamp!(KenyaCoVaccines.get_incidence_time_array_by_age_dose(sol, 9), 0.0, Inf)  # Sum over age and vaccine doses
    ι_all = clamp!(KenyaCoVaccines.get_incidence_time_array_by_age_dose(sol, 10), 0.0, Inf)
    return ι_first, ι_all
end

"""
    function project_daily_infections_by_age_and_vaccine_immune_escape(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date,immune_esc_cb; ω = 1 / 180)

Project daily infections, broken down by age and vaccination status, in a structure that can incoporate immune escape variants.

The specific effect of the immune escape variant is encoded in the
"""
function project_daily_infections_by_age_and_vaccine_immune_escape(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date, immune_esc_cb, cb_vac; ω = 1 / 180, num_vac_cats = 3)
    @unpack PCR_cases, sero_cases, init_sero, baseline_sero_array, PCR_array, sero_sensitivity, sero_specificity, p_symp, p_sus, hosp_rate_by_age, N, M_BB, prob, α, αP, γA, γM, relative_testing_rate, alpha_variant_time, delta_variant_time, model_start_time, model_end_time, janstarttime, VE_acquisition, VE_infectiousness, VE_severe_disease_risk, M_county_ho, M_county_school, M_county_work, M_county_other = model
    @unpack β₀, β_home, β_school, β_other, β_work, ϵ, χ, p_test, E₀, inc_R_αvar, time_scale_αvar, mid_point_αvar, inc_R_δvar, time_scale_δvar, mid_point_δvar, init_scale = θ  #
    # Vaccine efficacies
    ve_ac1 = VE_acquisition[1]
    ve_ac2 = VE_acquisition[2]
    ve_ac3 = VE_acquisition[3]
    ve_inf1 = VE_infectiousness[1]
    ve_inf2 = VE_infectiousness[2]
    ve_inf3 = VE_infectiousness[3]
    ve_dis1 = VE_severe_disease_risk[1]
    ve_dis2 = VE_severe_disease_risk[2]
    ve_dis3 = VE_severe_disease_risk[3]
    #Set transmission parameters --- SCHOOLS START SHUT IN THIS MODEL
    p = [β₀, β_home, 0.0, β_other, β_work, α, ϵ, αP, γA, γM, ω, inc_R_αvar, time_scale_αvar, alpha_variant_time + mid_point_αvar * 30, inc_R_δvar, time_scale_δvar, delta_variant_time + mid_point_δvar * 30, init_scale, ve_ac1, ve_ac2, ve_ac3, ve_inf1, ve_inf2, ve_inf3, ve_dis1, ve_dis2, ve_dis3, 1.0, 0.0, 0.0]
    u0 = KenyaCoVaccines.create_initial_conditions(θ, model; num_vac_cats = num_vac_cats)
    #Set up the callback for lockdown
    schoolsclosure_times = [Float64((Date(2021, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 12, 23) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 12, 23) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 12, 23) - Date(2020, 12, 1)).value)]
    function affect_schoolsclosed!(integrator)
        integrator.p[3] = 0.0 #Schools shut
    end
    schoolsclosed_cb = PresetTimeCallback(schoolsclosure_times, affect_schoolsclosed!, save_positions = (false, false))

    #Set up callback for reopening Schools
    schoolsreopening_times = [Float64((Date(2021, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 11) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 10, 11) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 10, 11) - Date(2020, 12, 1)).value)]
    function affect_schoolsopening!(integrator)
        integrator.p[3] = β_school # Schools open
    end
    schoolsreopening_cb = PresetTimeCallback(schoolsreopening_times, affect_schoolsopening!, save_positions = (false, false))
    cb = CallbackSet(schoolsclosed_cb, schoolsreopening_cb, immune_esc_cb, cb_vac)

    #Solve for daily incidence by age (not currently broken down by severity of disease)
    sim_duration = (projection_date - Date(2020, 12, 1)).value
    sol = solve(prob, CVODE_BDF(linear_solver = :GMRES);# abstol = 1e-15,reltol = 1e-15,
        tspan = (0, sim_duration),
        callback = cb,
        u0 = u0, p = p,
        saveat = 1,
        # isoutofdomain = (u, p, t) -> any(x -> x < -1e-6, u),
        verbose = true,
        force_dtmin = true)
    # global solves
    # solves += 1
    # println(solves)
    # sol = solve(prob, Tsit5();tspan = (0, sim_duration),AutoTsit5(Rodas5())
    #     callback = cb,
    #     u0 = u0, p = p,
    #     saveat = 1,
    #     # isoutofdomain = (u, p, t) -> any(x -> x < -1e-6, u),
    #     verbose = true)

    ι_first = clamp!(KenyaCoVaccines.get_incidence_time_array_by_age_dose(sol, 9), 0.0, Inf)  # Sum over age and vaccine doses
    ι_all = clamp!(KenyaCoVaccines.get_incidence_time_array_by_age_dose(sol, 10), 0.0, Inf)
    if size(ι_first, 1) < 912
        ι_first = vcat(ι_first, zeros(912 - size(ι_first, 1), 6, 3))
    end
    if size(ι_all, 1) < 912
        ι_all = vcat(ι_first, zeros(912 - size(ι_first, 1), 6, 3))
    end

    return ι_first, ι_all
end

"""
    function solve_model(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date; ω = 1 / 180)

Return `ODESolution` object. Simulation use the parameters in `θ` with waning immunity rate `ω`.
NB: this saves every aspect of the simulation, and, therefore, costs a large amount of memory is not used temporarily.
"""
function solve_model(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date; ω = 1 / 180)
    @unpack PCR_cases, sero_cases, init_sero, baseline_sero_array, PCR_array, sero_sensitivity, sero_specificity, p_symp, p_sus, hosp_rate_by_age, N, M_BB, prob, α, αP, γA, γM, relative_testing_rate, alpha_variant_time, delta_variant_time, model_start_time, model_end_time, janstarttime, VE_acquisition, VE_infectiousness, VE_severe_disease_risk, M_county_ho, M_county_school, M_county_work, M_county_other = model
    @unpack β₀, β_home, β_school, β_other, β_work, ϵ, χ, p_test, E₀, inc_R_αvar, time_scale_αvar, mid_point_αvar, inc_R_δvar, time_scale_δvar, mid_point_δvar, init_scale = θ  #
    # Vaccine efficacies
    ve_ac1 = VE_acquisition[1]
    ve_ac2 = VE_acquisition[2]
    ve_ac3 = VE_acquisition[3]
    ve_inf1 = VE_infectiousness[1]
    ve_inf2 = VE_infectiousness[2]
    ve_inf3 = VE_infectiousness[3]
    ve_dis1 = VE_severe_disease_risk[1]
    ve_dis2 = VE_severe_disease_risk[2]
    ve_dis3 = VE_severe_disease_risk[3]
    #Set transmission parameters --- SCHOOLS START SHUT IN THIS MODEL
    p = [β₀, β_home, 0.0, β_other, β_work, α, ϵ, αP, γA, γM, ω, inc_R_αvar, time_scale_αvar, alpha_variant_time + mid_point_αvar * 30, inc_R_δvar, time_scale_δvar, delta_variant_time + mid_point_δvar * 30, init_scale, ve_ac1, ve_ac2, ve_ac3, ve_inf1, ve_inf2, ve_inf3, ve_dis1, ve_dis2, ve_dis3]
    u0 = KenyaCoVaccines.create_initial_conditions(θ, model)
    #Set up the callback for lockdown
    schoolsclosure_times = [Float64((Date(2021, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 12, 23) - Date(2020, 12, 1)).value)]
    function affect_schoolsclosed!(integrator)
        integrator.p[3] = 0.0 #Schools shut
    end
    schoolsclosed_cb = PresetTimeCallback(schoolsclosure_times, affect_schoolsclosed!, save_positions = (false, false))

    #Set up callback for reopening Schools
    schoolsreopening_times = [Float64((Date(2021, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 11) - Date(2020, 12, 1)).value)]
    function affect_schoolsopening!(integrator)
        integrator.p[3] = β_school # Schools open #
    end
    schoolsreopening_cb = PresetTimeCallback(schoolsreopening_times, affect_schoolsopening!, save_positions = (false, false))

    cb = CallbackSet(schoolsclosed_cb, schoolsreopening_cb)

    #Solve for daily incidence by age (not currently broken down by severity of disease)
    sim_duration = (projection_date - Date(2020, 12, 1)).value
    sol = solve(prob, AutoTsit5(Rosenbrock23()); abstol = 1e-3,
        tspan = (0, sim_duration),
        callback = cb,
        u0 = u0, p = p,
        saveat = 1,
        isoutofdomain = (u, p, t) -> any(x -> x < -1e-6, u),
        verbose = true)

    return sol

end

function solve_model_immune_escape(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date, immune_esc_cb, vac_cb, solver; ω = 1 / 180)
    @unpack PCR_cases, sero_cases, init_sero, baseline_sero_array, PCR_array, sero_sensitivity, sero_specificity, p_symp, p_sus, hosp_rate_by_age, N, M_BB, prob, α, αP, γA, γM, relative_testing_rate, alpha_variant_time, delta_variant_time, model_start_time, model_end_time, janstarttime, VE_acquisition, VE_infectiousness, VE_severe_disease_risk, M_county_ho, M_county_school, M_county_work, M_county_other = model
    @unpack β₀, β_home, β_school, β_other, β_work, ϵ, χ, p_test, E₀, inc_R_αvar, time_scale_αvar, mid_point_αvar, inc_R_δvar, time_scale_δvar, mid_point_δvar, init_scale = θ  #
    # Vaccine efficacies
    ve_ac1 = VE_acquisition[1]
    ve_ac2 = VE_acquisition[2]
    ve_ac3 = VE_acquisition[3]
    ve_inf1 = VE_infectiousness[1]
    ve_inf2 = VE_infectiousness[2]
    ve_inf3 = VE_infectiousness[3]
    ve_dis1 = VE_severe_disease_risk[1]
    ve_dis2 = VE_severe_disease_risk[2]
    ve_dis3 = VE_severe_disease_risk[3]
    #Set transmission parameters --- SCHOOLS START SHUT IN THIS MODEL
    p = [β₀, β_home, 0.0, β_other, β_work, α, ϵ, αP, γA, γM, ω, inc_R_αvar, time_scale_αvar, alpha_variant_time + mid_point_αvar * 30, inc_R_δvar, time_scale_δvar, delta_variant_time + mid_point_δvar * 30, init_scale, ve_ac1, ve_ac2, ve_ac3, ve_inf1, ve_inf2, ve_inf3, ve_dis1, ve_dis2, ve_dis3, 1.0, 0.0, 0.0]
    u0 = KenyaCoVaccines.create_initial_conditions(θ, model)
    #Set up the callback for lockdown
    schoolsclosure_times = [Float64((Date(2021, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 12, 23) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 12, 23) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 12, 23) - Date(2020, 12, 1)).value)]
    function affect_schoolsclosed!(integrator)
        integrator.p[3] = 0.0 #Schools shut
    end
    schoolsclosed_cb = PresetTimeCallback(schoolsclosure_times, affect_schoolsclosed!, save_positions = (false, false))

    #Set up callback for reopening Schools
    schoolsreopening_times = [Float64((Date(2021, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 11) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 10, 11) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 10, 11) - Date(2020, 12, 1)).value)]
    function affect_schoolsopening!(integrator)
        integrator.p[3] = β_school # Schools open #
    end
    schoolsreopening_cb = PresetTimeCallback(schoolsreopening_times, affect_schoolsopening!, save_positions = (false, false))
    cb = CallbackSet(schoolsclosed_cb, schoolsreopening_cb, immune_esc_cb, cb_vac)

    #Solve for daily incidence by age (not currently broken down by severity of disease)
    sim_duration = (projection_date - Date(2020, 12, 1)).value
    sol = solve(prob, solver; abstol = 1e-3, reltol = 1e-3,
        tspan = (0, sim_duration),
        callback = cb,
        u0 = u0, p = p,
        saveat = 1,
        # isoutofdomain = (u, p, t) -> any(x -> x < -1e-6, u),
        verbose = true)

    return sol

end

function solve_model_waning_vac(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date, immune_esc_cb, vac_cb, solver,
    VE_acq, VE_inf, VE_sev; ω=1 / 180)
    @unpack PCR_cases, sero_cases, init_sero, baseline_sero_array, PCR_array, sero_sensitivity, sero_specificity, p_symp, p_sus, hosp_rate_by_age, N, M_BB, prob, α, αP, γA, γM, relative_testing_rate, alpha_variant_time, delta_variant_time, model_start_time, model_end_time, janstarttime, VE_acquisition, VE_infectiousness, VE_severe_disease_risk, M_county_ho, M_county_school, M_county_work, M_county_other = model
    @unpack β₀, β_home, β_school, β_other, β_work, ϵ, χ, p_test, E₀, inc_R_αvar, time_scale_αvar, mid_point_αvar, inc_R_δvar, time_scale_δvar, mid_point_δvar, init_scale = θ  #
    # Vaccine efficacies
    ve_ac1 = VE_acq[1]
    ve_ac2 = VE_acq[2]
    ve_ac3 = VE_acq[3]
    ve_inf1 = VE_inf[1]
    ve_inf2 = VE_inf[2]
    ve_inf3 = VE_inf[3]
    ve_dis1 = VE_sev[1]
    ve_dis2 = VE_sev[2]
    ve_dis3 = VE_sev[3]
    #Set transmission parameters --- SCHOOLS START SHUT IN THIS MODEL
    p = [β₀, β_home, 0.0, β_other, β_work, α, ϵ, αP, γA, γM, ω, inc_R_αvar, time_scale_αvar, alpha_variant_time + mid_point_αvar * 30, inc_R_δvar, time_scale_δvar, delta_variant_time + mid_point_δvar * 30, init_scale, ve_ac1, ve_ac2, ve_ac3, ve_inf1, ve_inf2, ve_inf3, ve_dis1, ve_dis2, ve_dis3, 1.0, 0.0, 0.0]
    u0 = KenyaCoVaccines.create_initial_conditions_with_waning(θ, model)
    #Set up the callback for lockdown
    schoolsclosure_times = [Float64((Date(2021, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 12, 23) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 12, 23) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 12, 23) - Date(2020, 12, 1)).value)]
    function affect_schoolsclosed!(integrator)
        integrator.p[3] = 0.0 #Schools shut
    end
    schoolsclosed_cb = PresetTimeCallback(schoolsclosure_times, affect_schoolsclosed!, save_positions=(false, false))

    #Set up callback for reopening Schools
    schoolsreopening_times = [Float64((Date(2021, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 11) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2022, 10, 11) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2023, 10, 11) - Date(2020, 12, 1)).value)]
    function affect_schoolsopening!(integrator)
        integrator.p[3] = β_school # Schools open #
    end
    schoolsreopening_cb = PresetTimeCallback(schoolsreopening_times, affect_schoolsopening!, save_positions=(false, false))
    cb = CallbackSet(schoolsclosed_cb, schoolsreopening_cb, immune_esc_cb, vac_cb)

    #Solve for daily incidence by age (not currently broken down by severity of disease)
    sim_duration = (projection_date - Date(2020, 12, 1)).value
    sol = solve(prob, solver; #abstol = 1e-3, reltol = 1e-3,
        tspan=(0, sim_duration),
        callback=cb,
        u0=u0, p=p,
        saveat=1,
        # isoutofdomain = (u, p, t) -> any(x -> x < -1e-6, u),
        verbose=true)

    return sol

end

"""
    function project_daily_infections_by_age_and_vaccine_waning_vac(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date,immune_esc_cb; ω = 1 / 180)

Project daily infections, broken down by age and vaccination status, in a structure that can incoporate immune escape variants, and waning immunity.

The specific effect of the immune escape variant is encoded in the `immune_esc_cb` callback object. The specific vaccination programme is encoded in the `cb_vac` callback object.
"""
function project_daily_infections_by_age_and_vaccine_waning_vac(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date, immune_esc_cb, cb_vac, solver,
    VE_acq, VE_inf, VE_sev; ω=1 / 180, data_length::Integer=912)

    sol = solve_model_waning_vac(θ, model, projection_date, immune_esc_cb, cb_vac, solver,
        VE_acq, VE_inf, VE_sev; ω=ω)

    ι_all = clamp!(KenyaCoVaccines.get_incidence_time_array_by_age_dose(sol, 10), 0.0, Inf) #All infections
    ι_dis_weighted = clamp!(KenyaCoVaccines.get_incidence_time_array_by_age_dose(sol, 11), 0.0, Inf) #Weighted by infection risk

    if size(ι_all, 1) < data_length
        ι_all = vcat(ι_all, zeros(912 - size(ι_all, 1), size(ι_all, 2), size(ι_all, 3)))
    end

    if size(ι_dis_weighted, 1) < data_length
        ι_dis_weighted = vcat(ι_dis_weighted, zeros(912 - size(ι_dis_weighted, 1), size(ι_dis_weighted, 2), size(ι_dis_weighted, 3)))
    end

    return ι_all, ι_dis_weighted
end

"""
    function project_daily_infections_first_by_age(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date; ω = 1 / 180)

Return total first (i.e. no reinfections) infections as a timeseries, this sums over vaccine status, but seperates by age. Simulation use the parameters in `θ` with waning immunity rate `ω`.
"""
function project_daily_infections_first_by_age(θ, model::KenyaCoVaccines.CoVAreaModel, projection_date::Date; ω = 1 / 180)
    @unpack PCR_cases, sero_cases, init_sero, baseline_sero_array, PCR_array, sero_sensitivity, sero_specificity, p_symp, p_sus, hosp_rate_by_age, N, M_BB, prob, α, αP, γA, γM, relative_testing_rate, alpha_variant_time, delta_variant_time, model_start_time, model_end_time, janstarttime, VE_acquisition, VE_infectiousness, VE_severe_disease_risk, M_county_ho, M_county_school, M_county_work, M_county_other = model
    @unpack β₀, β_home, β_school, β_other, β_work, ϵ, χ, p_test, E₀, inc_R_αvar, time_scale_αvar, mid_point_αvar, inc_R_δvar, time_scale_δvar, mid_point_δvar, init_scale = θ  #
    # Vaccine efficacies
    ve_ac1 = VE_acquisition[1]
    ve_ac2 = VE_acquisition[2]
    ve_ac3 = VE_acquisition[3]
    ve_inf1 = VE_infectiousness[1]
    ve_inf2 = VE_infectiousness[2]
    ve_inf3 = VE_infectiousness[3]
    ve_dis1 = VE_severe_disease_risk[1]
    ve_dis2 = VE_severe_disease_risk[2]
    ve_dis3 = VE_severe_disease_risk[3]
    #Set transmission parameters --- SCHOOLS START SHUT IN THIS MODEL
    p = [β₀, β_home, 0.0, β_other, β_work, α, ϵ, αP, γA, γM, ω, inc_R_αvar, time_scale_αvar, alpha_variant_time + mid_point_αvar * 30, inc_R_δvar, time_scale_δvar, delta_variant_time + mid_point_δvar * 30, init_scale, ve_ac1, ve_ac2, ve_ac3, ve_inf1, ve_inf2, ve_inf3, ve_dis1, ve_dis2, ve_dis3]
    u0 = KenyaCoVaccines.create_initial_conditions(θ, model)
    #Set up the callback for lockdown
    schoolsclosure_times = [Float64((Date(2021, 3, 19) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 16) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 1) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 12, 23) - Date(2020, 12, 1)).value)]
    function affect_schoolsclosed!(integrator)
        integrator.p[3] = 0.0 #Schools shut
    end
    schoolsclosed_cb = PresetTimeCallback(schoolsclosure_times, affect_schoolsclosed!, save_positions = (false, false))

    #Set up callback for reopening Schools
    schoolsreopening_times = [Float64((Date(2021, 1, 4) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 5, 10) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 7, 26) - Date(2020, 12, 1)).value),
        Float64((Date(2021, 10, 11) - Date(2020, 12, 1)).value)]
    function affect_schoolsopening!(integrator)
        integrator.p[3] = β_school # Schools open #
    end
    schoolsreopening_cb = PresetTimeCallback(schoolsreopening_times, affect_schoolsopening!, save_positions = (false, false))

    cb = CallbackSet(schoolsclosed_cb, schoolsreopening_cb)

    #Solve for daily incidence by age (not currently broken down by severity of disease)
    sim_duration = (projection_date - Date(2020, 12, 1)).value
    sol = solve(prob, BS3(); abstol = 1e-3,
        tspan = (0, sim_duration),
        callback = cb,
        u0 = u0, p = p,
        saveat = 1,
        isoutofdomain = (u, p, t) -> any(x -> x < 0, u),
        verbose = false)
    ι = clamp!(KenyaCoVaccines.get_incidence_time_array_by_age(sol, 9), 0.0, Inf)  # Sum over vaccine doses
    return ι
end

"""
    function perc_pos_to_num_tests_scale(model;lookback=21)

Fit a linear realationship between the % PCR positive (7 day running average) and the number of PCR tests performed on each day (7-day running average),
over the last `lookback` days of available data. Returns a `NamedTuple` with the regression coefficients and the std over error.
"""
function perc_pos_to_num_tests_scale(model::KenyaCoVaccines.CoVAreaModel; lookback = 21)
    prop_pos = model.PCR_cases[:, 1] ./ model.PCR_cases[:, 2]
    prop_pos_mv_av = [mean(prop_pos[(t-3):t+3]) for t = 4:(length(prop_pos)-3)]
    number_pos_mv_av = [mean(model.PCR_cases[:, 1][(t-3):t+3]) for t = 4:(length(model.PCR_cases[:, 1])-3)]
    number_tests_mv_av = [mean(model.PCR_cases[:, 2][(t-3):t+3]) for t = 4:(length(model.PCR_cases[:, 2])-3)]

    idxs = (.~isnan.(prop_pos_mv_av)) .& (prop_pos_mv_av .> 0) .& (1:length(prop_pos_mv_av) .>= length(prop_pos_mv_av) - lookback) .& (number_tests_mv_av .> 0)
    xs = prop_pos_mv_av[idxs]
    ys = number_tests_mv_av[idxs]
    β = [ones(length(xs)) xs] \ ys
    pred_ys = [ones(length(xs)) xs] * β
    std_ys = std(ys .- pred_ys)
    return (β = β, σ = std_ys)
end

function generate_PCR_predictions(ι::AbstractVector, θ, model::KenyaCoVaccines.CoVAreaModel, fit; M_PCR = 30.0, p_test_scale = 5e-4)
    @unpack PCR_cases, sero_cases, init_sero, baseline_sero_array, PCR_array, sero_sensitivity, sero_specificity, p_symp, p_sus, hosp_rate_by_age, N, M_BB, prob, α, αP, γA, γM, relative_testing_rate, alpha_variant_time, delta_variant_time, model_start_time, model_end_time, janstarttime, VE_acquisition, VE_infectiousness, VE_severe_disease_risk, M_county_ho, M_county_school, M_county_work, M_county_other = model
    @unpack β₀, β_home, β_school, β_other, β_work, ϵ, χ, p_test, E₀, inc_R_αvar, time_scale_αvar, mid_point_αvar, inc_R_δvar, time_scale_δvar, mid_point_δvar, init_scale = θ  #

    PCR = KenyaCoVaccines.simple_conv(ι, PCR_array)
    p_PCR_pred_num = p_test_scale .* p_test .* PCR
    p_PCR_pred_denom = p_PCR_pred_num .+ p_test_scale .* ((p_test / χ) .* (sum(N) .- PCR))
    p_PCR_pred = (p_PCR_pred_num ./ p_PCR_pred_denom) .+ 0.001
    _PCR_cases = [PCR_cases[model_start_time:end, :]; fill(-1, length(ι) - size(PCR_cases[model_start_time:end, :], 1), 2)]
    n_tests = _PCR_cases[:, 2]
    for i = 1:length(n_tests)
        if n_tests[i] < 0
            n_tests[i] = round(Int64, fit.β[1] + fit.β[2] * p_PCR_pred[i] + fit.σ * randn())
        end
    end
    clamp!(n_tests, 0, Inf)
    sampled_PCR_pos_tests = [rand(BetaBinomial(n_tests[i], p_PCR_pred[i] * M_PCR, (1.0 - p_PCR_pred[i]) * M_PCR)) for i = 1:length(n_tests)]
    return [sampled_PCR_pos_tests n_tests]
end

function generate_sero_predictions(ι_sero::AbstractVector, θ, age_grp, model::KenyaCoVaccines.CoVAreaModel; seroreversionrate = 0.0)
    @unpack PCR_cases, sero_cases, init_sero, baseline_sero_array, PCR_array, sero_sensitivity, sero_specificity, p_symp, p_sus, hosp_rate_by_age, N, M_BB, prob, α, αP, γA, γM, relative_testing_rate, alpha_variant_time, delta_variant_time, model_start_time, model_end_time, janstarttime, VE_acquisition, VE_infectiousness, VE_severe_disease_risk, M_county_ho, M_county_school, M_county_work, M_county_other = model
    @unpack β₀, β_home, β_school, β_other, β_work, ϵ, χ, p_test, E₀, inc_R_αvar, time_scale_αvar, mid_point_αvar, inc_R_δvar, time_scale_δvar, mid_point_δvar, init_scale = θ  #

    sero_array = vcat(baseline_sero_array[1:30], [(1 - seroreversionrate)^k for k = 1:600])
    rand_sensitivity = rand(Beta(M_BB * sero_sensitivity, M_BB * (1 - sero_sensitivity)))
    unscaled_sero = (init_sero[age_grp] * init_scale) .+ KenyaCoVaccines.simple_conv(ι_sero, sero_array) ./ N[age_grp]
    sero = rand_sensitivity .* unscaled_sero .+ (1 - sero_specificity) .* (1 .- unscaled_sero)

    return sero
end

"""
    function get_variant_prediction(θ, model, n::Integer)

For times `1:n` give the model prediction of frequency of different variants conditional on parameter `θ`.
"""
function get_variant_prediction(θ, model, n::Integer)
    @unpack PCR_cases, sero_cases, init_sero, baseline_sero_array, PCR_array, sero_sensitivity, sero_specificity, p_symp, p_sus, hosp_rate_by_age, N, M_BB, prob, α, αP, γA, γM, relative_testing_rate, alpha_variant_time, delta_variant_time, model_start_time, model_end_time, janstarttime, VE_acquisition, VE_infectiousness, VE_severe_disease_risk, M_county_ho, M_county_school, M_county_work, M_county_other = model
    @unpack β₀, β_home, β_school, β_other, β_work, ϵ, χ, p_test, E₀, inc_R_αvar, time_scale_αvar, mid_point_αvar, inc_R_δvar, time_scale_δvar, mid_point_δvar, init_scale = θ  #

    prop_inc_delta = [KenyaCoVaccines.gen_logistic(t; κ = time_scale_δvar, x₀ = delta_variant_time + mid_point_δvar * 30) for t = 1:n]
    prop_inc_alpha_beta = (1 .- prop_inc_delta) .* [KenyaCoVaccines.gen_logistic(t; κ = time_scale_αvar, x₀ = alpha_variant_time + mid_point_αvar * 30) for t = 1:n]
    prop_inc_wt = (1 .- prop_inc_delta) .* (1 .- [KenyaCoVaccines.gen_logistic(t; κ = time_scale_αvar, x₀ = alpha_variant_time + mid_point_αvar * 30) for t = 1:n])

    return (prop_inc_wt = prop_inc_wt, prop_inc_alpha_beta = prop_inc_alpha_beta, prop_inc_delta = prop_inc_delta)
end

"""
    function weekly_proportions_with_errors(X::AbstractMatrix; lwr = 0.025, upr = 0.975)

Generate a time series of weekly proportions of pos/neg test results.
"""
function weekly_proportions_with_errors(X::AbstractMatrix; lwr = 0.025, upr = 0.975)
    @assert size(X, 2) == 2 "This method only works for two column matrices: first column being positive tests and second column being negative tests"

    weekly_pos_tests = [sum(X[grp, 1]) for grp in Iterators.partition(1:size(X, 1), 7)]
    weekly_neg_tests = [sum(X[grp, 2]) for grp in Iterators.partition(1:size(X, 1), 7)]
    weekly_total_tests = weekly_pos_tests + weekly_neg_tests

    ts = collect(0:7:(7*(length(weekly_total_tests)-1)))
    prop = [weekly_total_tests[i] > 0 ? weekly_pos_tests[i] / weekly_total_tests[i] : 0.0 for i = 1:length(weekly_total_tests)]
    lb = [weekly_total_tests[i] > 0 ? prop[i] - invlogcdf(Beta(0.5 + weekly_pos_tests[i], 0.5 + weekly_neg_tests[i]), log(lwr)) : prop[i] for i = 1:length(weekly_total_tests)]
    ub = [weekly_total_tests[i] > 0 ? invlogcdf(Beta(0.5 + weekly_pos_tests[i], 0.5 + weekly_neg_tests[i]), log(upr)) - prop[i] : 1 - prop[i] for i = 1:length(weekly_total_tests)]
    idxs = weekly_total_tests .> 0
    return (ts = ts[idxs], prop = prop[idxs], lb = lb[idxs], ub = ub[idxs])
end

"""
    function monthly_proportions_with_errors(X::AbstractMatrix, dates::Vector{Date}, date0::Date; lwr = 0.025, upr = 0.975)

Generate a time series of monthly proportions of pos/neg test results.
"""
function monthly_proportions_with_errors(X::AbstractMatrix, dates::Vector{Date}, date0::Date; lwr = 0.025, upr = 0.975)
    @assert size(X, 2) == 2 "This method only works for two column matrices: first column being positive tests and second column being negative tests"

    samplemonths = month.(dates)
    samplyears = year.(dates)

    ts = Float64[]
    monthly_pos_tests = Int64[]
    monthly_neg_tests = Int64[]

    for y = minimum(samplyears):maximum(samplyears)
        for m = minimum(samplemonths):maximum(samplemonths)
            idxs = (samplemonths .== m) .& (samplyears .== y)
            if any(idxs)
                day1 = Date(y, m, 1)
                dayend = lastdayofmonth(day1)
                push!(ts, mean([(day1 - date0).value, (dayend - date0).value]))
                push!(monthly_pos_tests, sum(X[idxs, 1]))
                push!(monthly_neg_tests, sum(X[idxs, 2]))
            end
        end
    end

    monthly_total_tests = monthly_pos_tests + monthly_neg_tests
    idxs = monthly_total_tests .> 0
    prop = [monthly_total_tests[i] > 0 ? monthly_pos_tests[i] / monthly_total_tests[i] : 0.0 for i = 1:length(monthly_total_tests)]
    lb = [monthly_total_tests[i] > 0 ? prop[i] - invlogcdf(Beta(0.5 + monthly_pos_tests[i], 0.5 + monthly_neg_tests[i]), log(lwr)) : prop[i] for i = 1:length(monthly_total_tests)]
    ub = [monthly_total_tests[i] > 0 ? invlogcdf(Beta(0.5 + monthly_pos_tests[i], 0.5 + monthly_neg_tests[i]), log(upr)) - prop[i] : 1 - prop[i] for i = 1:length(monthly_total_tests)]

    return (ts = ts[idxs], prop = prop[idxs], lb = lb[idxs], ub = ub[idxs])
end

"""
    function get_credible_intervals(X::AbstractMatrix; lwr = 0.025, upr = 0.975)

Get credible intervals of a matrix `X` over time where dim 1 is assumed to be time, and dim 2 is the variation over MCMC samples.
"""
function get_credible_intervals(X::AbstractMatrix; lwr = 0.025, upr = 0.975)
    pred = mean(X, dims = 2)
    lb = pred .- [quantile(X[t, :], lwr) for t = 1:size(X, 1)]
    ub = [quantile(X[t, :], upr) for t = 1:size(X, 1)] .- pred
    return (pred = pred, lb = lb, ub = ub)
end

"""
    function get_credible_intervals(X::Vector{AbstractArray}; lwr = 0.025, upr = 0.975)

Converts the vector of arrays into a matrix then applies `get_credible_intervals`.
"""
function get_credible_intervals(X::AbstractVector; lwr = 0.025, upr = 0.975)
    X = Matrix(RecursiveArrayTools.vecvec_to_mat(X)')
    return get_credible_intervals(X; lwr = lwr, upr = upr)
end

"""
    function mvav_cols(X::AbstractMatrix)

Return 7-day central moving average for each column in Matrix `X`.
"""
function mvav_cols(X::AbstractMatrix)
    _X = zeros(size(X))
    for j = 1:size(X, 2)
        _X[:, j] .= [X[1:3, j]
            [mean(X[(t-3):(t+3), j]) for t = 4:(size(X, 1)-3)]
            X[(end-2):end, j]]
    end
    return _X
end

"""
    function smooth_vector(X)

Return 7-day central moving average for vector `X`.
"""
function smooth_vector(X::AbstractVector)
    return [mean(X[1:3])
        [mean(X[(t-3):(t+3)]) for t = 4:(length(X)-3)]
        mean(X[(end-2):end])]
end

function vecmat_to_array(X)
    _X = zeros(size(X[1], 1), size(X[1], 2), length(X))
    for k = 1:length(X)
        _X[:, :, k] .= X[k]
    end
    return permutedims(_X, [1, 3, 2])
end

function gather_infection_type_for_econ_analysis(model, projection_date::Date;
    mort_scale::Float64, mort_age, mort_var,
    H_ICU_to_death_fit, H_hosp, H_var_ab, H_var_delta,
    veff1 = 0.24, veff2 = 0.187, veff_death1 = 0.2, veff_death2 = 0.1)

    infections_by_age_dose = model |> flatten_chain .|> θ -> project_daily_infections_by_age_and_vaccine(θ, model, Date(2023, 6, 1))
    age_symptoms = model.p_symp

    uncertaintyinfections_by_age_unvac_first = Matrix.([vecvec_to_mat([infs[1][:, a, 1] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_unvac_reinf = Matrix.([vecvec_to_mat([infs[2][:, a, 1] .- infs[1][:, a, 1] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_1dose_first = Matrix.([vecvec_to_mat([infs[1][:, a, 2] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_1dose_reinf = Matrix.([vecvec_to_mat([infs[2][:, a, 2] .- infs[1][:, a, 2] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_2dose_first = Matrix.([vecvec_to_mat([infs[1][:, a, 3] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_2dose_reinf = Matrix.([vecvec_to_mat([infs[2][:, a, 3] .- infs[1][:, a, 3] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])

    θs = model |> flatten_chain
    n = size(uncertaintyinfections_by_age_unvac_first[1], 1)
    var_frqs = θs .|> θ -> get_variant_prediction(θ, model, n)
    wt_frqs = Matrix(vecvec_to_mat([data.prop_inc_wt for data in var_frqs])')
    alpha_beta_frqs = Matrix(vecvec_to_mat([data.prop_inc_alpha_beta for data in var_frqs])')
    delta_frqs = Matrix(vecvec_to_mat([data.prop_inc_delta for data in var_frqs])')

    var_deadliness = wt_frqs .+ mort_var[1] .* alpha_beta_frqs .+ mort_var[2] .* delta_frqs
    var_severe_risk = H_hosp .* wt_frqs .+ H_var_ab .* mort_var[1] .* alpha_beta_frqs .+ H_var_delta .* mort_var[2] .* delta_frqs
    var_critical_risk = H_ICU_to_death_fit .* var_deadliness

    all_infs = [uncertaintyinfections_by_age_unvac_first[a] + uncertaintyinfections_by_age_unvac_reinf[a] +
                uncertaintyinfections_by_age_1dose_first[a] + uncertaintyinfections_by_age_1dose_reinf[a] +
                uncertaintyinfections_by_age_2dose_first[a] + uncertaintyinfections_by_age_2dose_reinf[a]
                for a = 1:6]

    temp_deadly_infs = [uncertaintyinfections_by_age_unvac_first[a] + 0.15 * uncertaintyinfections_by_age_unvac_reinf[a] +
                        veff_death1 * uncertaintyinfections_by_age_1dose_first[a] + veff_death1 * 0.15 * uncertaintyinfections_by_age_1dose_reinf[a] +
                        veff_death2 * uncertaintyinfections_by_age_2dose_first[a] + veff_death2 * 0.15 * uncertaintyinfections_by_age_2dose_reinf[a]
                        for a = 1:6] .* mort_scale .* mort_age

    temp_sev_crit_infs = [uncertaintyinfections_by_age_unvac_first[a] + 0.15 * uncertaintyinfections_by_age_unvac_reinf[a] +
                          veff1 * uncertaintyinfections_by_age_1dose_first[a] + veff1 * 0.15 * uncertaintyinfections_by_age_1dose_reinf[a] +
                          veff2 * uncertaintyinfections_by_age_2dose_first[a] + veff2 * 0.15 * uncertaintyinfections_by_age_2dose_reinf[a]
                          for a = 1:6] .* mort_scale .* mort_age

    asymptomatic_infs = vecmat_to_array(all_infs .* (1.0 .- age_symptoms))
    deadly_infs = vecmat_to_array(temp_deadly_infs .|> x -> x .* var_deadliness)
    severe_infs = vecmat_to_array(temp_sev_crit_infs .|> x -> x .* var_severe_risk)
    critical_infs = vecmat_to_array(temp_sev_crit_infs .|> x -> x .* var_critical_risk)
    mild_infections = vecmat_to_array((all_infs .* age_symptoms)) .- deadly_infs .- severe_infs .- critical_infs

    return (; asymptomatic_infs, mild_infections, severe_infs, critical_infs, deadly_infs)

end

function gather_infection_type_for_econ_analysis_waning_vac(model, projection_date::Date, cb_imm_esc, cb_vac, solver, VE_acqs, VE_infs, VE_sevs, VE_deaths;
    mort_scale::Float64, mort_age, mort_var, H_ICU_to_death_fit, H_hosp, H_var_ab, H_var_delta, VE_waning=0.0)

    #Solve the infection predictions for each parameter set
    θs = model |> flatten_chain
    # infections_by_age_dose = θs .|> θ -> project_daily_infections_by_age_and_vaccine_waning_vac(θ, model, projection_date, cb_imm_esc, cb_vac, solver)
    infections_by_age_dose = map((θ, VE_acq, VE_inf, VE_sev) -> project_daily_infections_by_age_and_vaccine_waning_vac(θ, model, Date(2023, 6, 1), cb_imm_esc, cb_vac, solver, VE_acq, VE_inf, VE_sev),
        θs, VE_acqs, VE_infs, VE_sevs)
    #Gather information about size of data arrays
    age_symptoms = model.p_symp
    n_t = size(infections_by_age_dose[1][2], 1)
    n_a = size(infections_by_age_dose[1][2], 2)
    n_v = size(infections_by_age_dose[1][2], 3)
    #Solve the variant frequency predictions
    n = size(infections_by_age_dose[1][1], 1)
    var_frqs = θs .|> θ -> get_variant_prediction(θ, model, n)
    wt_frqs = Matrix(vecvec_to_mat([data.prop_inc_wt for data in var_frqs])')
    alpha_beta_frqs = Matrix(vecvec_to_mat([data.prop_inc_alpha_beta for data in var_frqs])')
    delta_frqs = Matrix(vecvec_to_mat([data.prop_inc_delta for data in var_frqs])')
    #Weight by variant severity, this is weighted average of risk per infection by variant frequency
    var_deadliness = wt_frqs .+ mort_var[1] .* alpha_beta_frqs .+ mort_var[2] .* delta_frqs
    var_severe_risk = H_hosp .* wt_frqs .+ H_var_ab .* mort_var[1] .* alpha_beta_frqs .+ H_var_delta .* mort_var[2] .* delta_frqs
    var_critical_risk = H_ICU_to_death_fit .* var_deadliness
    #Solve number of asymptomatic infections
    symptomatic_infs = permutedims(vecmat_to_array(infections_by_age_dose .|> x -> sum(x[1], dims=3)[:, :] .* age_symptoms'), [1, 3, 2])
    asymptomatic_infs = permutedims(vecmat_to_array(infections_by_age_dose .|> x -> sum(x[1], dims=3)[:, :] .* (1.0 .- age_symptoms)'), [1, 3, 2])
    #Gather infections weighted by reduction in risk because a reinfection (from simulations)
    weighted_infs = [infs[2][t, a, v] for t = 1:n_t, a = 1:n_a, v = 1:n_v, infs in infections_by_age_dose]
    #Sensitivity over possible vaccine effectiveness range
    VE_sevs_mat = 1 .- [vecvec_to_mat(VE_sevs) vecvec_to_mat(VE_sevs)[:, end] VE_waning .* vecvec_to_mat(VE_sevs)[:, end]]
    VE_deaths_mat = 1 .- [vecvec_to_mat(VE_deaths) vecvec_to_mat(VE_deaths)[:, end] VE_waning .* vecvec_to_mat(VE_deaths)[:, end]]
    #Fully weight infections by risk from all factors
    #Indices: t - day, a - age group, v - vaccination status, s - sample over MCMC draws of underlying parameters
    @tullio severe_infs[t, a, s] := mort_scale * weighted_infs[t, a, v, s] * mort_age[a] * VE_sevs_mat[s, v] * var_severe_risk[t, s]
    @tullio critical_infs[t, a, s] := mort_scale * weighted_infs[t, a, v, s] * mort_age[a] * VE_sevs_mat[s, v] * var_critical_risk[t, s]
    @tullio deadly_infs[t, a, s] := mort_scale * weighted_infs[t, a, v, s] * mort_age[a] * VE_deaths_mat[s, v] * var_deadliness[t, s]
    #Left over symptomatic infections are mild
    mild_infections = symptomatic_infs .- deadly_infs .- severe_infs .- critical_infs

    return (; asymptomatic_infs, mild_infections, severe_infs, critical_infs, deadly_infs)

end

function gather_infection_type_for_econ_analysis_immune_escape(model, projection_date::Date, cb_imm_esc, cb_vac;
    mort_scale::Float64, mort_age, mort_var,
    H_ICU_to_death_fit, H_hosp, H_var_ab, H_var_delta,
    veff1 = 0.24, veff2 = 0.187, veff_death1 = 0.2, veff_death2 = 0.1)

    infections_by_age_dose = model |> flatten_chain .|> θ -> project_daily_infections_by_age_and_vaccine_immune_escape(θ, model, projection_date, cb_imm_esc, cb_vac)
    age_symptoms = model.p_symp

    uncertaintyinfections_by_age_unvac_first = Matrix.([vecvec_to_mat([infs[1][:, a, 1] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_unvac_reinf = Matrix.([vecvec_to_mat([infs[2][:, a, 1] .- infs[1][:, a, 1] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_1dose_first = Matrix.([vecvec_to_mat([infs[1][:, a, 2] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_1dose_reinf = Matrix.([vecvec_to_mat([infs[2][:, a, 2] .- infs[1][:, a, 2] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_2dose_first = Matrix.([vecvec_to_mat([infs[1][:, a, 3] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_2dose_reinf = Matrix.([vecvec_to_mat([infs[2][:, a, 3] .- infs[1][:, a, 3] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])

    θs = model |> flatten_chain
    n = size(uncertaintyinfections_by_age_unvac_first[1], 1)
    var_frqs = θs .|> θ -> get_variant_prediction(θ, model, n)
    wt_frqs = Matrix(vecvec_to_mat([data.prop_inc_wt for data in var_frqs])')
    alpha_beta_frqs = Matrix(vecvec_to_mat([data.prop_inc_alpha_beta for data in var_frqs])')
    delta_frqs = Matrix(vecvec_to_mat([data.prop_inc_delta for data in var_frqs])')

    var_deadliness = wt_frqs .+ mort_var[1] .* alpha_beta_frqs .+ mort_var[2] .* delta_frqs
    var_severe_risk = H_hosp .* wt_frqs .+ H_var_ab .* mort_var[1] .* alpha_beta_frqs .+ H_var_delta .* mort_var[2] .* delta_frqs
    var_critical_risk = H_ICU_to_death_fit .* var_deadliness

    all_infs = [uncertaintyinfections_by_age_unvac_first[a] + uncertaintyinfections_by_age_unvac_reinf[a] +
                uncertaintyinfections_by_age_1dose_first[a] + uncertaintyinfections_by_age_1dose_reinf[a] +
                uncertaintyinfections_by_age_2dose_first[a] + uncertaintyinfections_by_age_2dose_reinf[a]
                for a = 1:6]

    temp_deadly_infs = [uncertaintyinfections_by_age_unvac_first[a] + 0.15 * uncertaintyinfections_by_age_unvac_reinf[a] +
                        veff_death1 * uncertaintyinfections_by_age_1dose_first[a] + veff_death1 * 0.15 * uncertaintyinfections_by_age_1dose_reinf[a] +
                        veff_death2 * uncertaintyinfections_by_age_2dose_first[a] + veff_death2 * 0.15 * uncertaintyinfections_by_age_2dose_reinf[a]
                        for a = 1:6] .* mort_scale .* mort_age

    temp_sev_crit_infs = [uncertaintyinfections_by_age_unvac_first[a] + 0.15 * uncertaintyinfections_by_age_unvac_reinf[a] +
                          veff1 * uncertaintyinfections_by_age_1dose_first[a] + veff1 * 0.15 * uncertaintyinfections_by_age_1dose_reinf[a] +
                          veff2 * uncertaintyinfections_by_age_2dose_first[a] + veff2 * 0.15 * uncertaintyinfections_by_age_2dose_reinf[a]
                          for a = 1:6] .* mort_scale .* mort_age

    asymptomatic_infs = vecmat_to_array(all_infs .* (1.0 .- age_symptoms))
    deadly_infs = vecmat_to_array(temp_deadly_infs .|> x -> x .* var_deadliness)
    severe_infs = vecmat_to_array(temp_sev_crit_infs .|> x -> x .* var_severe_risk)
    critical_infs = vecmat_to_array(temp_sev_crit_infs .|> x -> x .* var_critical_risk)
    mild_infections = vecmat_to_array((all_infs .* age_symptoms)) .- deadly_infs .- severe_infs .- critical_infs

    return (; asymptomatic_infs, mild_infections, severe_infs, critical_infs, deadly_infs)

end

function gather_infection_type_for_econ_analysis_waning_vac(model, projection_date::Date, cb_imm_esc, cb_vac;
    mort_scale::Float64, mort_age, mort_var,
    H_ICU_to_death_fit, H_hosp, H_var_ab, H_var_delta,
    veff1 = 0.24, veff2 = 0.187, veff_death1 = 0.2, veff_death2 = 0.1)

    infections_by_age_dose = model |> flatten_chain .|> θ -> project_daily_infections_by_age_and_vaccine_immune_escape(θ, model, projection_date, cb_imm_esc, cb_vac)
    age_symptoms = model.p_symp

    uncertaintyinfections_by_age_unvac_first = Matrix.([vecvec_to_mat([infs[1][:, a, 1] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_unvac_reinf = Matrix.([vecvec_to_mat([infs[2][:, a, 1] .- infs[1][:, a, 1] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_1dose_first = Matrix.([vecvec_to_mat([infs[1][:, a, 2] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_1dose_reinf = Matrix.([vecvec_to_mat([infs[2][:, a, 2] .- infs[1][:, a, 2] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_2dose_first = Matrix.([vecvec_to_mat([infs[1][:, a, 3] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])
    uncertaintyinfections_by_age_2dose_reinf = Matrix.([vecvec_to_mat([infs[2][:, a, 3] .- infs[1][:, a, 3] for infs in infections_by_age_dose])' for a = 1:size(infections_by_age_dose[1][1], 2)])

    θs = model |> flatten_chain
    n = size(uncertaintyinfections_by_age_unvac_first[1], 1)
    var_frqs = θs .|> θ -> get_variant_prediction(θ, model, n)
    wt_frqs = Matrix(vecvec_to_mat([data.prop_inc_wt for data in var_frqs])')
    alpha_beta_frqs = Matrix(vecvec_to_mat([data.prop_inc_alpha_beta for data in var_frqs])')
    delta_frqs = Matrix(vecvec_to_mat([data.prop_inc_delta for data in var_frqs])')

    var_deadliness = wt_frqs .+ mort_var[1] .* alpha_beta_frqs .+ mort_var[2] .* delta_frqs
    var_severe_risk = H_hosp .* wt_frqs .+ H_var_ab .* mort_var[1] .* alpha_beta_frqs .+ H_var_delta .* mort_var[2] .* delta_frqs
    var_critical_risk = H_ICU_to_death_fit .* var_deadliness

    all_infs = [uncertaintyinfections_by_age_unvac_first[a] + uncertaintyinfections_by_age_unvac_reinf[a] +
                uncertaintyinfections_by_age_1dose_first[a] + uncertaintyinfections_by_age_1dose_reinf[a] +
                uncertaintyinfections_by_age_2dose_first[a] + uncertaintyinfections_by_age_2dose_reinf[a]
                for a = 1:6]

    temp_deadly_infs = [uncertaintyinfections_by_age_unvac_first[a] + 0.15 * uncertaintyinfections_by_age_unvac_reinf[a] +
                        veff_death1 * uncertaintyinfections_by_age_1dose_first[a] + veff_death1 * 0.15 * uncertaintyinfections_by_age_1dose_reinf[a] +
                        veff_death2 * uncertaintyinfections_by_age_2dose_first[a] + veff_death2 * 0.15 * uncertaintyinfections_by_age_2dose_reinf[a]
                        for a = 1:6] .* mort_scale .* mort_age

    temp_sev_crit_infs = [uncertaintyinfections_by_age_unvac_first[a] + 0.15 * uncertaintyinfections_by_age_unvac_reinf[a] +
                          veff1 * uncertaintyinfections_by_age_1dose_first[a] + veff1 * 0.15 * uncertaintyinfections_by_age_1dose_reinf[a] +
                          veff2 * uncertaintyinfections_by_age_2dose_first[a] + veff2 * 0.15 * uncertaintyinfections_by_age_2dose_reinf[a]
                          for a = 1:6] .* mort_scale .* mort_age

    asymptomatic_infs = vecmat_to_array(all_infs .* (1.0 .- age_symptoms))
    deadly_infs = vecmat_to_array(temp_deadly_infs .|> x -> x .* var_deadliness)
    severe_infs = vecmat_to_array(temp_sev_crit_infs .|> x -> x .* var_severe_risk)
    critical_infs = vecmat_to_array(temp_sev_crit_infs .|> x -> x .* var_critical_risk)
    mild_infections = vecmat_to_array((all_infs .* age_symptoms)) .- deadly_infs .- severe_infs .- critical_infs

    return (; asymptomatic_infs, mild_infections, severe_infs, critical_infs, deadly_infs)

end

function add_infs(infs1, infs2)
    return (asymptomatic_infs = infs1.asymptomatic_infs + infs2.asymptomatic_infs,
        mild_infections = infs1.mild_infections + infs2.mild_infections,
        severe_infs = infs1.severe_infs + infs2.severe_infs,
        critical_infs = infs1.critical_infs + infs2.critical_infs,
        deadly_infs = infs1.deadly_infs + infs2.deadly_infs)
end


text:
"""
`MCIntegrationVI(;ϵ::T=1e-5,nMC::Integer=1000,optimizer::Optimizer=Adam(α=0.1))`

Variational Inference solver by approximating gradients via MC Integration.

**Keyword arguments**

    - `ϵ::T` : convergence criteria, which can be user defined
    - `nMC::Int` : Number of samples per data point for the integral evaluation
    - `optimizer::Optimizer` : Optimizer used for the variational updates. Should be an Optimizer object from the [GradDescent.jl]() package. Default is `Adam()`
"""
mutable struct MCIntegrationVI{T<:Real,N} <: NumericalVI{T}
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    nMC::Int64 #Number of samples for MC Integrations
    Stochastic::Bool #Use of mini-batches
    nSamples::Int64 #Number of samples of the data
    nMinibatch::Int64 #Size of mini-batches
    ρ::T #Stochastic Coefficient
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    vi_opt::NTuple{N,NVIOptimizer}
    MBIndices::Vector #Indices of the minibatch
    xview::SubArray{T,2,Matrix{T}}
    yview::SubArray
    function MCIntegrationVI{T}(ϵ::T,nMC::Integer,optimizer::Opt,Stochastic::Bool,nSamplesUsed::Integer=1) where {T<:Real,Opt<:Optimizer}
        return new{T,1}(ϵ,0,nMC,Stochastic,1,nSamplesUsed)
    end
end

function MCIntegrationVI(;ϵ::T=1e-5,nMC::Integer=1000,optimizer::Optimizer=Momentum(η=0.01)) where {T<:Real}
    MCIntegrationVI{T}(ϵ,nMC,optimizer,false,1)
end

"""
`MCIntegrationSVI(;ϵ::T=1e-5,nMC::Integer=1000,optimizer::Optimizer=Adam(α=0.1))`

Stochastic Variational Inference solver by approximating gradients via Monte Carlo integration

**Argument**

    -`nMinibatch::Integer` : Number of samples per mini-batches

**Keyword arguments**

    - `ϵ::T` : convergence criteria, which can be user defined
    - `nMC::Int` : Number of samples per data point for the integral evaluation
    - `optimizer::Optimizer` : Optimizer used for the variational updates. Should be an Optimizer object from the [GradDescent.jl]() package. Default is `Adam()`
"""
function MCIntegrationSVI(nMinibatch::Integer;ϵ::T=1e-5,nMC::Integer=200,optimizer::Optimizer=Momentum(η=0.001)) where {T<:Real}
    MCIntegrationVI{T}(ϵ,nMC,0,optimizer,true,nMinibatch)
end

function compute_grad_expectations!(model::VGP{T,L,<:MCIntegrationVI}) where {T,L}
    raw_samples = randn(model.inference.nMC,model.nLatent)
    samples = similar(raw_samples)
    for i in 1:model.nSample
        samples .= raw_samples.*[sqrt(model.Σ[k][i,i]) for k in 1:model.nLatent]' .+ [model.μ[k][i] for k in 1:model.nLatent]'
        grad_samples(model,samples,i)
    end
end

function compute_grad_expectations!(model::SVGP{T,L,<:MCIntegrationVI}) where {T,L}
    raw_samples = randn(model.inference.nMC,model.nLatent)
    samples = similar(raw_samples)
    Σ = opt_diag.(model.κ.*model.Σ,model.κ)
    μ = model.κ.*model.μ
    for i in 1:model.nSample
        samples .= raw_samples.*[sqrt(Σ[k][i]) for k in 1:model.nLatent]' .+ [μ[k][i] for k in 1:model.nLatent]'
        grad_samples(model,samples,i)
    end
end

function expec_log_likelihood(l::Likelihood,i::MCIntegrationVI,y,μ,diag_cov) where {T,L}
    raw_samples = randn(i.nMC,i.nLatent)
    samples = similar(raw_samples)
    loglike = 0.0
    for i in 1:model.nSample
        samples .= raw_samples.*[sqrt(diag_cov[k][i]) for k in 1:i.nLatent]' .+ [μ[k][i] for k in 1:i.nLatent]'
        loglike += mean(mapslices(f->logpdf(l,y[i],f),samples,dims=2))
    end
    return loglike
end

ext: jl


text:
"""
    unsafe_clone(feature::Feature)

Duplicate feature.

The newly created feature is owned by the caller, and will have its own
reference to the OGRFeatureDefn.
"""
unsafe_clone(feature::Feature)::Feature = Feature(GDAL.ogr_f_clone(feature.ptr))

"""
    destroy(feature::Feature)

Destroy the feature passed in.

The feature is deleted, but within the context of the GDAL/OGR heap. This is
necessary when higher level applications use GDAL/OGR from a DLL and they want
to delete a feature created within the DLL. If the delete is done in the calling
application the memory will be freed onto the application heap which is
inappropriate.
"""
function destroy(feature::Feature)::Nothing
    GDAL.ogr_f_destroy(feature.ptr)
    feature.ptr = C_NULL
    return nothing
end

"""
    setgeom!(feature::Feature, geom::AbstractGeometry)

Set feature geometry.

This method updates the features geometry, and operate exactly as
SetGeometryDirectly(), except that this method does not assume ownership of the
passed geometry, but instead makes a copy of it.

### Parameters
* `feature`: the feature on which new geometry is applied to.
* `geom`: the new geometry to apply to feature.

### Returns
`OGRERR_NONE` if successful, or `OGR_UNSUPPORTED_GEOMETRY_TYPE` if the geometry
type is illegal for the `OGRFeatureDefn` (checking not yet implemented).
"""
function setgeom!(feature::Feature, geom::AbstractGeometry)::Feature
    result = GDAL.ogr_f_setgeometry(feature.ptr, geom.ptr)
    @ogrerr result "OGRErr $result: Failed to set feature geometry."
    return feature
end

"""
    getgeom(feature::Feature)

Returns a clone of the geometry corresponding to the feature.
"""
function getgeom(feature::Feature)::IGeometry
    result = GDAL.ogr_f_getgeometryref(feature.ptr)
    return if result == C_NULL
        IGeometry()
    else
        IGeometry(GDAL.ogr_g_clone(result))
    end
end

function unsafe_getgeom(feature::Feature)::Geometry
    result = GDAL.ogr_f_getgeometryref(feature.ptr)
    return if result == C_NULL
        Geometry()
    else
        Geometry(GDAL.ogr_g_clone(result))
    end
end

"""
    nfield(feature::Feature)

Fetch number of fields on this feature.

This will always be the same as the field count for the OGRFeatureDefn.
"""
nfield(feature::Feature)::Integer = GDAL.ogr_f_getfieldcount(feature.ptr)

"""
    getfielddefn(feature::Feature, i::Integer)

Fetch definition for this field.

### Parameters
* `feature`: the feature on which the field is found.
* `i`: the field to fetch, from 0 to GetFieldCount()-1.

### Returns
an handle to the field definition (from the `FeatureDefn`). This is an
internal reference, and should not be deleted or modified.
"""
getfielddefn(feature::Feature, i::Integer)::IFieldDefnView =
    IFieldDefnView(GDAL.ogr_f_getfielddefnref(feature.ptr, i))

"""
    findfieldindex(feature::Feature, name::Union{AbstractString, Symbol})

Fetch the field index given field name.

### Parameters
* `feature`: the feature on which the field is found.
* `name`: the name of the field to search for.

### Returns
the field index, or `nothing` if no matching field is found.

### Remarks
This is a cover for the `OGRFeatureDefn::GetFieldIndex()` method.
"""
function findfieldindex(
    feature::Feature,
    name::Union{AbstractString,Symbol},
)::Union{Integer,Nothing}
    i = GDAL.ogr_f_getfieldindex(feature.ptr, name)
    return if i == -1
        nothing
    else
        i
    end
end

"""
    isfieldset(feature::Feature, i::Integer)

Test if a field has ever been assigned a value or not.

### Parameters
* `feature`: the feature that owned the field.
* `i`: the field to fetch, from 0 to GetFieldCount()-1.
"""
isfieldset(feature::Feature, i::Integer)::Bool =
    Bool(GDAL.ogr_f_isfieldset(feature.ptr, i))

"""
    unsetfield!(feature::Feature, i::Integer)

Clear a field, marking it as unset.

### Parameters
* `feature`: the feature that owned the field.
* `i`: the field to fetch, from 0 to GetFieldCount()-1.
"""
function unsetfield!(feature::Feature, i::Integer)::Feature
    GDAL.ogr_f_unsetfield(feature.ptr, i)
    return feature
end

"""
    isfieldnull(feature::Feature, i::Integer)

Test if a field is null.

### Parameters
* `feature`: the feature that owned the field.
* `i`: the field to test, from 0 to GetFieldCount()-1.

### Returns
`true` if the field is null, otherwise `false`.

### References
* https://gdal.org/development/rfc/rfc67_nullfieldvalues.html
"""
isfieldnull(feature::Feature, i::Integer)::Bool =
    Bool(GDAL.ogr_f_isfieldnull(feature.ptr, i))

"""
    isfieldsetandnotnull(feature::Feature, i::Integer)

Test if a field is set and not null.

### Parameters
* `feature`: the feature that owned the field.
* `i`: the field to test, from 0 to GetFieldCount()-1.

### Returns
`true` if the field is set and not null, otherwise `false`.

### References
* https://gdal.org/development/rfc/rfc67_nullfieldvalues.html
"""
isfieldsetandnotnull(feature::Feature, i::Integer)::Bool =
    Bool(GDAL.ogr_f_isfieldsetandnotnull(feature.ptr, i))

"""
    setfieldnull!(feature::Feature, i::Integer)

Clear a field, marking it as null.

### Parameters
* `feature`: the feature that owned the field.
* `i`: the field to set to null, from 0 to GetFieldCount()-1.

### References
* https://gdal.org/development/rfc/rfc67_nullfieldvalues.html
"""
function setfieldnull!(feature::Feature, i::Integer)::Feature
    GDAL.ogr_f_setfieldnull(feature.ptr, i)
    return feature
end

# """
#     OGR_F_GetRawFieldRef(OGRFeatureH hFeat,
#                          int iField) -> OGRField *
# Fetch an handle to the internal field value given the index.
# ### Parameters
# * `hFeat`: handle to the feature on which field is found.
# * `iField`: the field to fetch, from 0 to GetFieldCount()-1.
# ### Returns
# the returned handle is to an internal data structure, and should not be freed,
# or modified.
# """
# function getrawfieldref(arg1::Ptr{OGRFeatureH},arg2::Integer)
#     ccall((:OGR_F_GetRawFieldRef,libgdal),Ptr{OGRField},(Ptr{OGRFeatureH},
#           Cint),arg1,arg2)
# end

"""
    asint(feature::Feature, i::Integer)

Fetch field value as integer.

### Parameters
* `feature`: the feature that owned the field.
* `i`: the field to fetch, from 0 to GetFieldCount()-1.
"""
asint(feature::Feature, i::Integer)::Int32 =
    GDAL.ogr_f_getfieldasinteger(feature.ptr, i)

"""
    asint64(feature::Feature, i::Integer)

Fetch field value as integer 64 bit.

### Parameters
* `feature`: the feature that owned the field.
* `i`: the field to fetch, from 0 to GetFieldCount()-1.
"""
asint64(feature::Feature, i::Integer)::Int64 =
    GDAL.ogr_f_getfieldasinteger64(feature.ptr, i)

"""
    asdouble(feature::Feature, i::Integer)

Fetch field value as a double.

### Parameters
* `feature`: the feature that owned the field.
* `i`: the field to fetch, from 0 to GetFieldCount()-1.
"""
asdouble(feature::Feature, i::Integer)::Float64 =
    GDAL.ogr_f_getfieldasdouble(feature.ptr, i)

"""
    asstring(feature::Feature, i::Integer)

Fetch field value as a string.

### Parameters
* `feature`: the feature that owned the field.
* `i`: the field to fetch, from 0 to GetFieldCount()-1.
"""
asstring(feature::Feature, i::Integer)::String =
    GDAL.ogr_f_getfieldasstring(feature.ptr, i)

"""
    asintlist(feature::Feature, i::Integer)

Fetch field value as a list of integers.

### Parameters
* `hFeat`: handle to the feature that owned the field.
* `iField`: the field to fetch, from 0 to GetFieldCount()-1.
* `pnCount`: an integer to put the list count (number of integers) into.

### Returns
the field value. This list is internal, and should not be modified, or freed.
Its lifetime may be very brief. If *pnCount is zero on return the returned
pointer may be NULL or non-NULL.
"""
function asintlist(feature::Feature, i::Integer)::Vector{Int32}
    n = Ref{Cint}()
    ptr = GDAL.ogr_f_getfieldasintegerlist(feature.ptr, i, n)
    return (n.x == 0) ? Int32[] : unsafe_wrap(Vector{Int32}, ptr, n.x)
end

"""
    asint64list(feature::Feature, i::Integer)

Fetch field value as a list of 64 bit integers.

### Parameters
* `hFeat`: handle to the feature that owned the field.
* `iField`: the field to fetch, from 0 to GetFieldCount()-1.
* `pnCount`: an integer to put the list count (number of integers) into.

### Returns
the field value. This list is internal, and should not be modified, or freed.
Its lifetime may be very brief. If *pnCount is zero on return the returned
pointer may be NULL or non-NULL.
"""
function asint64list(feature::Feature, i::Integer)::Vector{Int64}
    n = Ref{Cint}()
    ptr = GDAL.ogr_f_getfieldasinteger64list(feature.ptr, i, n)
    return (n.x == 0) ? Int64[] : unsafe_wrap(Vector{Int64}, ptr, n.x)
end

"""
    asdoublelist(feature::Feature, i::Integer)

Fetch field value as a list of doubles.

### Parameters
* `hFeat`: handle to the feature that owned the field.
* `iField`: the field to fetch, from 0 to GetFieldCount()-1.
* `pnCount`: an integer to put the list count (number of doubles) into.

### Returns
the field value. This list is internal, and should not be modified, or freed.
Its lifetime may be very brief. If *pnCount is zero on return the returned
pointer may be NULL or non-NULL.
"""
function asdoublelist(feature::Feature, i::Integer)::Vector{Float64}
    n = Ref{Cint}()
    ptr = GDAL.ogr_f_getfieldasdoublelist(feature.ptr, i, n)
    return (n.x == 0) ? Float64[] : unsafe_wrap(Vector{Float64}, ptr, n.x)
end

"""
    asstringlist(feature::Feature, i::Integer)

Fetch field value as a list of strings.

### Parameters
* `hFeat`: handle to the feature that owned the field.
* `iField`: the field to fetch, from 0 to GetFieldCount()-1.

### Returns
the field value. This list is internal, and should not be modified, or freed.
Its lifetime may be very brief.
"""
asstringlist(feature::Feature, i::Integer)::Vector{String} =
    GDAL.ogr_f_getfieldasstringlist(feature.ptr, i)

"""
    asbinary(feature::Feature, i::Integer)

Fetch field value as binary.

### Parameters
* `hFeat`: handle to the feature that owned the field.
* `iField`: the field to fetch, from 0 to GetFieldCount()-1.

### Returns
the field value. This list is internal, and should not be modified, or freed.
Its lifetime may be very brief.
"""
function asbinary(feature::Feature, i::Integer)::Vector{UInt8}
    n = Ref{Cint}()
    ptr = GDAL.ogr_f_getfieldasbinary(feature.ptr, i, n)
    return (n.x == 0) ? UInt8[] : unsafe_wrap(Vector{UInt8}, ptr, n.x)
end

"""
    asdatetime(feature::Feature, i::Integer)

Fetch field value as date and time. Currently this method only works for
OFTDate, OFTTime and OFTDateTime fields.

### Parameters
* `hFeat`: handle to the feature that owned the field.
* `iField`: the field to fetch, from 0 to GetFieldCount()-1.

### Returns
`true` on success or `false` on failure.
"""
function asdatetime(feature::Feature, i::Integer)::DateTime
    pyr = Ref{Cint}()
    pmth = Ref{Cint}()
    pday = Ref{Cint}()
    phr = Ref{Cint}()
    pmin = Ref{Cint}()
    psec = Ref{Cint}()
    ptz = Ref{Cint}()
    result = Bool(
        GDAL.ogr_f_getfieldasdatetime(
            feature.ptr,
            i,
            pyr,
            pmth,
            pday,
            phr,
            pmin,
            psec,
            ptz,
        ),
    )
    (result == false) && error("Failed to fetch datetime at index $i")
    return DateTime(pyr[], pmth[], pday[], phr[], pmin[], psec[])
end

# """
#     OGR_F_GetFieldAsDateTimeEx(OGRFeatureH hFeat,
#                                int iField,
#                                int * pnYear,
#                                int * pnMonth,
#                                int * pnDay,
#                                int * pnHour,
#                                int * pnMinute,
#                                float * pfSecond,
#                                int * pnTZFlag) -> int
# Fetch field value as date and time.
# ### Parameters
# * `hFeat`: handle to the feature that owned the field.
# * `iField`: the field to fetch, from 0 to GetFieldCount()-1.
# * `pnYear`: (including century)
# * `pnMonth`: (1-12)
# * `pnDay`: (1-31)
# * `pnHour`: (0-23)
# * `pnMinute`: (0-59)
# * `pfSecond`: (0-59 with millisecond accuracy)
# * `pnTZFlag`: (0=unknown, 1=localtime, 100=GMT, see data model for details)
# ### Returns
# `true` on success or `false` on failure.
# """
# function getfieldasdatetimeex(hFeat::Ptr{OGRFeatureH},iField::Integer,pnYear,
#                               pnMonth,pnDay,pnHour,pnMinute,pfSecond,pnTZFlag)
#     ccall((:OGR_F_GetFieldAsDateTimeEx,libgdal),Cint,(Ptr{OGRFeatureH},Cint,
#           Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cfloat},
#           Ptr{Cint}),hFeat,iField,pnYear,pnMonth,pnDay,pnHour,pnMinute,
#           pfSecond,pnTZFlag)
# end

function getdefault(feature::Feature, i::Integer)::Union{String,Nothing}
    return getdefault(getfielddefn(feature, i))
end

getfield(feature::Feature, i::Nothing)::Missing = missing

const _FETCHFIELD = Dict{OGRFieldType,Function}(
    OFTInteger => asint,
    OFTIntegerList => asintlist,
    OFTReal => asdouble,
    OFTRealList => asdoublelist,
    OFTString => asstring,
    OFTStringList => asstringlist,
    OFTBinary => asbinary,
    OFTDate => asdatetime,
    OFTTime => asdatetime,
    OFTDateTime => asdatetime,
    OFTInteger64 => asint64,
    OFTInteger64List => asint64list,
)

"""
    getfield(feature, i)

When the field is unset, it will return `nothing`. When the field is set but
null, it will return `missing`.

### References
* https://gdal.org/development/rfc/rfc53_ogr_notnull_default.html
* https://gdal.org/development/rfc/rfc67_nullfieldvalues.html
"""
function getfield(feature::Feature, i::Integer)
    return if !isfieldset(feature, i)
        nothing
    elseif isfieldnull(feature, i)
        missing
    else
        _fieldtype = getfieldtype(getfielddefn(feature, i))
        try
            _fetchfield = _FETCHFIELD[_fieldtype]
            _fetchfield(feature, i)
        catch e
            if e isa KeyError
                error(
                    "$_fieldtype not implemented in getfield, please report an issue to https://github.com/yeesian/ArchGDAL.jl/issues",
                )
            else
                rethrow(e)
            end
        end
    end
end

function getfield(feature::Feature, name::Union{AbstractString,Symbol})
    return getfield(feature, findfieldindex(feature, name))
end

"""
    setfield!(feature::Feature, i::Integer, value)
    setfield!(feature::Feature, i::Integer, value::DateTime, tzflag::Int = 0)

Set a feature's `i`-th field to `value`.

The following types for `value` are accepted: `Int32`, `Int64`, `Float64`,
`AbstractString`, or a `Vector` with those in it, as well as `Vector{UInt8}`.
For `DateTime` values, an additional keyword argument `tzflag` is accepted
(0=unknown, 1=localtime, 100=GMT, see data model for details).

OFTInteger, OFTInteger64 and OFTReal fields will be set directly. OFTString
fields will be assigned a string representation of the value, but not
necessarily taking into account formatting constraints on this field. Other
field types may be unaffected.

### Parameters
* `feature`: handle to the feature that owned the field.
* `i`: the field to fetch, from 0 to GetFieldCount()-1.
* `value`: the value to assign.
"""
function setfield! end

function setfield!(feature::Feature, i::Integer, value::Nothing)::Feature
    unsetfield!(feature, i)
    return feature
end

function setfield!(feature::Feature, i::Integer, value::Missing)::Feature
    if isnullable(getfielddefn(feature, i))
        setfieldnull!(feature, 1)
    else
        setfield!(feature, i, getdefault(feature, i))
    end
    return feature
end

function setfield!(feature::Feature, i::Integer, value::Int32)::Feature
    GDAL.ogr_f_setfieldinteger(feature.ptr, i, value)
    return feature
end

function setfield!(feature::Feature, i::Integer, value::Int16)::Feature
    GDAL.ogr_f_setfieldinteger(feature.ptr, i, value)
    return feature
end

function setfield!(feature::Feature, i::Integer, value::Bool)::Feature
    GDAL.ogr_f_setfieldinteger(feature.ptr, i, value)
    return feature
end

function setfield!(feature::Feature, i::Integer, value::Int64)::Feature
    GDAL.ogr_f_setfieldinteger64(feature.ptr, i, value)
    return feature
end

function setfield!(feature::Feature, i::Integer, value::Float32)::Feature
    GDAL.ogr_f_setfielddouble(feature.ptr, i, value)
    return feature
end

function setfield!(feature::Feature, i::Integer, value::Float64)::Feature
    GDAL.ogr_f_setfielddouble(feature.ptr, i, value)
    return feature
end

function setfield!(feature::Feature, i::Integer, value::AbstractString)::Feature
    GDAL.ogr_f_setfieldstring(feature.ptr, i, value)
    return feature
end

function setfield!(feature::Feature, i::Integer, value::Vector{Int32})::Feature
    GDAL.ogr_f_setfieldintegerlist(feature.ptr, i, length(value), value)
    return feature
end

function setfield!(feature::Feature, i::Integer, value::Vector{Int64})::Feature
    GDAL.ogr_f_setfieldinteger64list(feature.ptr, i, length(value), value)
    return feature
end

function setfield!(
    feature::Feature,
    i::Integer,
    value::Vector{Float64},
)::Feature
    GDAL.ogr_f_setfielddoublelist(feature.ptr, i, length(value), value)
    return feature
end

function setfield!(
    feature::Feature,
    i::Integer,
    value::Vector{T},
)::Feature where {T<:AbstractString}
    GDAL.ogr_f_setfieldstringlist(feature.ptr, i, value)
    return feature
end

# """
#     OGR_F_SetFieldRaw(OGRFeatureH hFeat,
#                       int iField,
#                       OGRField * psValue) -> void
# Set field.
# ### Parameters
# * `hFeat`: handle to the feature that owned the field.
# * `iField`: the field to fetch, from 0 to GetFieldCount()-1.
# * `psValue`: handle on the value to assign.
# """
# function setfieldraw(arg1::Ptr{OGRFeatureH},arg2::Integer,arg3)
#     ccall((:OGR_F_SetFieldRaw,libgdal),Void,(Ptr{OGRFeatureH},Cint,
#           Ptr{OGRField}),arg1,arg2,arg3)
# end

function setfield!(feature::Feature, i::Integer, value::Vector{UInt8})::Feature
    GDAL.ogr_f_setfieldbinary(feature.ptr, i, sizeof(value), value)
    return feature
end

function setfield!(
    feature::Feature,
    i::Integer,
    dt::DateTime,
    tzflag::Int = 0,
)::Feature
    GDAL.ogr_f_setfielddatetime(
        feature.ptr,
        i,
        Dates.year(dt),
        Dates.month(dt),
        Dates.day(dt),
        Dates.hour(dt),
        Dates.minute(dt),
        Dates.second(dt),
        tzflag,
    )
    return feature
end

"""
    ngeom(feature::Feature)

Fetch number of geometry fields on this feature.

This will always be the same as the geometry field count for OGRFeatureDefn.
"""
ngeom(feature::Feature)::Integer = GDAL.ogr_f_getgeomfieldcount(feature.ptr)

"""
    getgeomdefn(feature::Feature, i::Integer)

Fetch definition for this geometry field.

### Parameters
* `feature`: the feature on which the field is found.
* `i`: the field to fetch, from 0 to GetGeomFieldCount()-1.

### Returns
The field definition (from the OGRFeatureDefn). This is an
internal reference, and should not be deleted or modified.
"""
function getgeomdefn(feature::Feature, i::Integer)::IGeomFieldDefnView
    return IGeomFieldDefnView(GDAL.ogr_f_getgeomfielddefnref(feature.ptr, i))
end

"""
    findgeomindex(feature::Feature, name::Union{AbstractString, Symbol} = "")

Fetch the geometry field index given geometry field name.

### Parameters
* `feature`: the feature on which the geometry field is found.
* `name`: the name of the geometry field to search for. (defaults to \"\")

### Returns
the geometry field index, or -1 if no matching geometry field is found.

### Remarks
This is a cover for the `OGRFeatureDefn::GetGeomFieldIndex()` method.
"""
function findgeomindex(
    feature::Feature,
    name::Union{AbstractString,Symbol} = "",
)::Integer
    return GDAL.ogr_f_getgeomfieldindex(feature.ptr, name)
end

"""
    getgeom(feature::Feature, i::Integer)

Returns a clone of the feature geometry at index `i`.

### Parameters
* `feature`: the feature to get geometry from.
* `i`: geometry field to get.
"""
function getgeom(feature::Feature, i::Integer)::IGeometry
    result = GDAL.ogr_f_getgeomfieldref(feature.ptr, i)
    return if result == C_NULL
        IGeometry()
    else
        IGeometry(GDAL.ogr_g_clone(result))
    end
end

function unsafe_getgeom(feature::Feature, i::Integer)::Geometry
    result = GDAL.ogr_f_getgeomfieldref(feature.ptr, i)
    return if result == C_NULL
        Geometry()
    else
        Geometry(GDAL.ogr_g_clone(result))
    end
end

function getgeom(
    feature::Feature,
    name::Union{AbstractString,Symbol},
)::IGeometry
    i = findgeomindex(feature, name)
    return if i == -1
        IGeometry()
    else
        getgeom(feature, i)
    end
end

function unsafe_getgeom(
    feature::Feature,
    name::Union{AbstractString,Symbol},
)::Geometry
    i = findgeomindex(feature, name)
    return if i == -1
        Geometry()
    else
        unsafe_getgeom(feature, i)
    end
end

"""
    setgeom!(feature::Feature, i::Integer, geom::AbstractGeometry)

Set feature geometry of a specified geometry field.

This function updates the features geometry, and operate exactly as
SetGeometryDirectly(), except that this function does not assume ownership of
the passed geometry, but instead makes a copy of it.

### Parameters
* `feature`: the feature on which to apply the geometry.
* `i`: geometry field to set.
* `geom`: the new geometry to apply to feature.

### Returns
`OGRERR_NONE` if successful, or `OGR_UNSUPPORTED_GEOMETRY_TYPE` if the geometry
type is illegal for the `OGRFeatureDefn` (checking not yet implemented).
"""
function setgeom!(feature::Feature, i::Integer, geom::AbstractGeometry)::Feature
    result = GDAL.ogr_f_setgeomfield(feature.ptr, i, geom.ptr)
    @ogrerr result "OGRErr $result: Failed to set feature geometry"
    return feature
end

"""
    getfid(feature::Feature)

Get feature identifier.

### Returns
feature id or `OGRNullFID` (`-1`) if none has been assigned.
"""
getfid(feature::Feature) = GDAL.ogr_f_getfid(feature.ptr)::Integer

"""
    setfid!(feature::Feature, i::Integer)

Set the feature identifier.

### Parameters
* `feature`: handle to the feature to set the feature id to.
* `i`: the new feature identifier value to assign.

### Returns
On success OGRERR_NONE, or on failure some other value.
"""
function setfid!(feature::Feature, i::Integer)::Feature
    result = GDAL.ogr_f_setfid(feature.ptr, i)
    @ogrerr result "OGRErr $result: Failed to set FID $i"
    return feature
end

"""
    setfrom!(feature1::Feature, feature2::Feature, forgiving::Bool = false)
    setfrom!(feature1::Feature, feature2::Feature, indices::Vector{Cint},
        forgiving::Bool = false)

Set one feature from another.

### Parameters
* `feature1`: handle to the feature to set to.
* `feature2`: handle to the feature from which geometry, and field values
    will be copied.
* `indices`: indices of the destination feature's fields stored at the
    corresponding index of the source feature's fields. A value of `-1` should
    be used to ignore the source's field. The array should not be NULL and be
    as long as the number of fields in the source feature.
* `forgiving`: `true` if the operation should continue despite lacking output
    fields matching some of the source fields.

### Returns
OGRERR_NONE if the operation succeeds, even if some values are not transferred,
otherwise an error code.
"""
function setfrom! end

function setfrom!(
    feature1::Feature,
    feature2::Feature,
    forgiving::Bool = false,
)::Feature
    result = GDAL.ogr_f_setfrom(feature1.ptr, feature2.ptr, forgiving)
    @ogrerr result "OGRErr $result: Failed to set feature"
    return feature1
end

function setfrom!(
    feature1::Feature,
    feature2::Feature,
    indices::Vector{Cint},
    forgiving::Bool = false,
)::Feature
    result = GDAL.ogr_f_setfromwithmap(
        feature1.ptr,
        feature2.ptr,
        forgiving,
        indices,
    )
    @ogrerr result "OGRErr $result: Failed to set feature with map"
    return feature1
end

"""
    getstylestring(feature::Feature)

Fetch style string for this feature.
"""
getstylestring(feature::Feature)::String =
    GDAL.ogr_f_getstylestring(feature.ptr)

"""
    setstylestring!(feature::Feature, style::AbstractString)

Set feature style string.

This method operate exactly as `setstylestringdirectly!()` except that
it doesn't assume ownership of the passed string, but makes a copy of it.
"""
function setstylestring!(feature::Feature, style::AbstractString)::Feature
    GDAL.ogr_f_setstylestring(feature.ptr, style)
    return feature
end

"""
    getstyletable(feature::Feature)

Fetch style table for this feature.
"""
getstyletable(feature::Feature)::StyleTable =
    StyleTable(GDAL.ogr_f_getstyletable(feature.ptr))

"""
    setstyletable!(feature::Feature, styletable::StyleTable)

Set the style table for this feature.
"""
function setstyletable!(feature::Feature, styletable::StyleTable)::Feature
    GDAL.ogr_f_setstyletable(feature.ptr, styletable.ptr)
    return feature
end

"""
    getnativedata(feature::Feature)

Returns the native data for the feature.

The native data is the representation in a "natural" form that comes from the
driver that created this feature, or that is aimed at an output driver. The
native data may be in different format, which is indicated by
GetNativeMediaType().

Note that most drivers do not support storing the native data in the feature
object, and if they do, generally the NATIVE_DATA open option must be passed at
dataset opening.

The "native data" does not imply it is something more performant or powerful
than what can be obtained with the rest of the API, but it may be useful in
round-tripping scenarios where some characteristics of the underlying format
are not captured otherwise by the OGR abstraction.
"""
getnativedata(feature::Feature)::String = GDAL.ogr_f_getnativedata(feature.ptr)

"""
    setnativedata!(feature::Feature, data::AbstractString)

Sets the native data for the feature.

The native data is the representation in a "natural" form that comes from the
driver that created this feature, or that is aimed at an output driver. The
native data may be in different format, which is indicated by
GetNativeMediaType().
"""
function setnativedata!(feature::Feature, data::AbstractString)::Feature
    GDAL.ogr_f_setnativedata(feature.ptr, data)
    return feature
end

"""
    getmediatype(feature::Feature)

Returns the native media type for the feature.

The native media type is the identifier for the format of the native data. It
follows the IANA RFC 2045 (see https://en.wikipedia.org/wiki/Media_type),
e.g. \"application/vnd.geo+json\" for JSON.
"""
getmediatype(feature::Feature)::String =
    GDAL.ogr_f_getnativemediatype(feature.ptr)

"""
    setmediatype!(feature::Feature, mediatype::AbstractString)

Sets the native media type for the feature.

The native media type is the identifier for the format of the native data. It
follows the IANA RFC 2045 (see https://en.wikipedia.org/wiki/Media_type),
e.g. \"application/vnd.geo+json\" for JSON.
"""
function setmediatype!(feature::Feature, mediatype::AbstractString)::Feature
    GDAL.ogr_f_setnativemediatype(feature.ptr, mediatype)
    return feature
end

"""
    fillunsetwithdefault!(feature::Feature; notnull = true,
        options = StringList(C_NULL))

Fill unset fields with default values that might be defined.

### Parameters
* `feature`: handle to the feature.
* `notnull`: if we should fill only unset fields with a not-null constraint.
* `papszOptions`: unused currently. Must be set to `NULL`.

### References
* https://gdal.org/development/rfc/rfc53_ogr_notnull_default.html
"""
function fillunsetwithdefault!(
    feature::Feature;
    notnull::Bool = true,
    options = StringList(C_NULL),
)::Feature
    GDAL.ogr_f_fillunsetwithdefault(feature.ptr, notnull, options)
    return feature
end

"""
    validate(feature::Feature, flags::Integer, emiterror::Bool)

Validate that a feature meets constraints of its schema.

The scope of test is specified with the nValidateFlags parameter.

Regarding `OGR_F_VAL_WIDTH`, the test is done assuming the string width must be
interpreted as the number of UTF-8 characters. Some drivers might interpret the
width as the number of bytes instead. So this test is rather conservative (if it
fails, then it will fail for all interpretations).

### Parameters
* `feature`: handle to the feature to validate.
* `flags`: `OGR_F_VAL_ALL` or combination of `OGR_F_VAL_NULL`,
    `OGR_F_VAL_GEOM_TYPE`, `OGR_F_VAL_WIDTH` and
    `OGR_F_VAL_ALLOW_NULL_WHEN_DEFAULT` with `|` operator
* `emiterror`: `true` if a `CPLError()` must be emitted when a check fails

### Returns
`true` if all enabled validation tests pass.

### References
* https://gdal.org/development/rfc/rfc53_ogr_notnull_default.html
"""
validate(feature::Feature, flags::FieldValidation, emiterror::Bool)::Bool =
    Bool(GDAL.ogr_f_validate(feature.ptr, flags, emiterror))


text:
using Pipelines
using JobSchedulers, Dates

scheduler_start()
sleep(2)

prog = CmdProgram(
    inputs = ["PATH", "SHELL", "INPUT1", "INPUT2", "INPUT3"],
    # outputs = "OUTPUT_FILE",
    cmd = pipeline(`SHELL PATH -c INPUT1`,
                   `SHELL PATH -c INPUT2`,
                   `SHELL PATH -c INPUT3`
                   ),
)

inputs = Dict(
    "SHELL" => `nix develop`,
    "PATH" => ENV["PRJ_ROOT"] * "/devshell",
    "INPUT1" => `std run //tenzir//entrypoints:config-vast-prod`,
    "INPUT2" => `std run //cliche//entrypoints:example add 1 2`,
    "INPUT3" => `std run //makes//entrypoints:scriptEnv`,
)
#outputs = "OUTPUT_FILE" => "out" # save output to file

program_job = Job(prog, inputs; touch_run_id_file = false)

submit!(program_job)

while JobSchedulers.DataFrames.nrow(queue()) > 0
    sleep(2)
    all_queue()
    queue(all = true)
end
# run(`cat out.txt`)

max_stars_repo_name: GTrunSec/DevSecOps-Cells-Lab

max_issues_repo_name: GTrunSec/DevSecOps-cells

max_forks_repo_name: GTrunSec/DevSecOps-Cells-Lab


text:
#=##############################################################################
# DESCRIPTION
    Functions for database manipulation

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jun 2020
  * License   : MIT
=###############################################################################

"""
    `new_database(path::String; index_file::String=DEF_NDXFL, mkpath_optargs=[], prompt=true, v_lvl=0)`

Creates a new airfoil database in the given path.
"""
function new_database(path::String;                          # Path of new database
                        # DATABASE OPTIONS
                        index_file::String=DEF_NDXFL,        # Name of indexing file
                        # OUTPUT OPTIONS
                        mkpath_optargs=[], prompt=true, v_lvl=0)

    # Case that path already exists
    if isdir(path)

        remove = true # Determine whether to remove it or not
        if prompt
            inp = ""
            while !(inp in ["y", "n"])
              print("\t"^v_lvl*"Directory $path already exists. Remove? (y/n) ")
              inp = readline()[1:end]
            end
            remove = inp=="y"
        end

        if remove; rm(path; recursive=true); end;
    end

    # Create path if it doesn't already exists
    mkpath(path; mkpath_optargs...)

    # Create subdirectories
    for dir in [DIR_XY, DIR_CL, DIR_CD, DIR_CM, DIR_XUP, DIR_XLO]
        mkpath(joinpath(path, dir); mkpath_optargs...)
    end

    # Create default indexing file
    f = open(joinpath(path, index_file), "w")

    for (coli, (key, val)) in enumerate(LBLS)
        print(f, val[1])
        print(f, coli==length(LBLS) ? "\n" : ",")
    end

    close(f)

    return joinpath(path, index_file)
end

"""
    `new_entry(; database_path::String=$(def_database), index_file::String=$(DEF_NDXFL), lbls...)`

Adds a new entry to the database under `database_path` with indexing file
`index_file`. `lbls` are the labels of the entry (or columns of database), as
shown above (Fields is the name of the label):

$(_lbls2prettystring())

For instance, here is a dummy example on how to add an airfoil with only the
required fields: `db.new_entry(; airfoilname="NACA 0012", xyfile="naca0012.csv)`
"""
function new_entry(; database_path::String=def_database,
                        index_file::String=DEF_NDXFL,
                        check_duplicates::Bool=true,
                        lbls...)

    # Check that all required fields where given
    for field in RQRD_FIELDS
        if findfirst(x->x==field, keys(lbls))==nothing
            error("Required field $field was not given!")
        end
    end

    # Test that all fields are the correct type
    for (field, val) in LBLS
        if findfirst(x->x==field, keys(lbls)) != nothing
            try
                lbls[findfirst(x->x==field, keys(lbls))] :: val[3]
            catch e
                error("Expected type $(val[3]) under field $(field);"*
                        " got $(lbls[findfirst(x->x==field, keys(lbls))])"*
                        " (type $(typeof(lbls[findfirst(x->x==field, keys(lbls))])))."*
                        " Error: $(e)")
            end
        end
    end

    # Add missing optional fields
    this_lbls = [ (key,
                    findfirst(x->x==key, keys(lbls))==nothing ? val[2] : lbls[findfirst(x->x==key, keys(lbls))])
                   for (key, val) in LBLS]

    # Create default indexing file
    f = open(joinpath(database_path, index_file), "a")

    for (coli, key) in enumerate(keys(LBLS))
        val = this_lbls[findfirst(x->x[1]==key, this_lbls)][2]
        print(f, val)
        print(f, coli==length(LBLS) ? "\n" : ",")
    end

    close(f)
end

"""
    `new_entry(polar::AirfoilPrep.Polar; database_path::String=$(def_database), warn=true, index_file::String=$(DEF_NDXFL), lbls...)`

Adds a new entry to the database under `database_path` with indexing file
`index_file`. `lbls` are the labels of the entry (or columns of database), as
shown above (Fields is the name of the label).
"""
function new_entry(polar::ap.Polar; database_path::String=def_database,
                                    ext::String=".csv", warn=true, lbls...)

    # Get airfoil name
    if findfirst(x->x==:airfoilname, keys(lbls))==nothing
        error("Required field $(:airfoilname) was not given!")
    end

    rflname = replace("$(lbls[findfirst(x->x==:airfoilname, keys(lbls))])", " "=>"")

    # Add fields from polar
    this_lbls = Any[(key, val) for (key, val) in lbls]

    for (field, fun) in [ (:re, ap.get_Re), (:ma, ap.get_Ma),
                          (:npanels, ap.get_npanels), (:ncrit, ap.get_ncrit),
                          (:diff, (args...)->PTNL_I)]

        if findfirst(x->x==field, keys(lbls))==nothing # Check if field wasn't prescribed by user
            push!(this_lbls, (field, fun(polar)))
        end

    end

    # Build file names
    re = "$(ceil(Int, this_lbls[findfirst(x->x[1]==:re, this_lbls)][2]))"
    ma = "$(replace("$(this_lbls[findfirst(x->x[1]==:ma, this_lbls)][2])", "."=>"p"))"
    npanels = "$(replace("$(this_lbls[findfirst(x->x[1]==:npanels, this_lbls)][2])", "."=>"p"))"
    ncrit = "$(replace("$(this_lbls[findfirst(x->x[1]==:ncrit, this_lbls)][2])", "."=>"p"))"
    diff = "$(replace("$(this_lbls[findfirst(x->x[1]==:diff, this_lbls)][2])", "."=>"p"))"

    xyfile = rflname*"-npanels"*npanels*"-"*diff*ext
    suff = "-re"*re*"-ma"*ma*"-ncrit"*ncrit*"-"*diff*ext
    clfile = rflname*"-Cl"*suff
    cdfile = rflname*"-Cd"*suff
    cmfile = rflname*"-Cm"*suff
    cmfile = rflname*"-Cm"*suff
    xupsepfile = rflname*"-xupsep"*suff
    xlosepfile = rflname*"-xlosep"*suff

    # Create files and add them to the entry
    for (field, fun, file, dir, headers) in [ (:xyfile, ap.get_geometry, xyfile, DIR_XY, ["x", "y"]),
                                              (:clfile, ap.get_cl, clfile, DIR_CL, ["alpha", "cl"]),
                                              (:cdfile, ap.get_cd, cdfile, DIR_CD, ["alpha", "cd"]),
                                              (:cmfile, ap.get_cm, cmfile, DIR_CM, ["alpha", "cm"]),
                                              (:xupsepfile, ap.get_xsepup, xupsepfile, DIR_XUP, ["alpha", "xsep"]),
                                              (:xlosepfile, ap.get_xseplo, xlosepfile, DIR_XLO, ["alpha", "xsep"]),
                                            ]

        if findfirst(x->x==field, keys(lbls))==nothing # Check if field wasn't prescribed by user

            col1, col2 = fun(polar)

            fname = joinpath(database_path, dir, file)

            if isfile(fname) && warn
                @warn("Overwriting file $fname.")
            end

            # Create file
            f = open(fname, "w")

            println(f, headers[1], ",", headers[2])
            for rowi in 1:length(col1)
                println(f, col1[rowi], ",", col2[rowi])
            end

            close(f)

            # Add field
            push!(this_lbls, (field, file))
        end

    end

    return new_entry(; database_path=database_path, this_lbls...)
end