text:
function extra_doses(vacc_data,N_vacc)
    dose_cum = zeros(size(vacc_data))
    for age in 1:6
        dose_cum[:,age] = cumsum(vacc_data[:,age])
    end
    above_N = trues(size(dose_cum))
    for time in 1:size(dose_cum)[1]
        above_N[time,:] .= (dose_cum[time,:] .> N_vacc)
    end
    vacc_extra = vacc_data .- (vacc_data .* (1 .- above_N))  # set doses past full coverage time in agegroup to zero
    return(vacc_extra)
end

"""
Function to generate vaccination rates: The GoK plan taken more generally, not in Phases
"""

function vacc_rate_data(doses_daily,model_days,pre_vacc_period,age_dist_vacc,N_kenya,county,scenario,rapid)
    N = Array(N_kenya[:,county])
    N = vcat(sum(N[1:4]),sum(N[5:10]),sum(N[11:12]),sum(N[13:14]),sum(N[15:16]),N[17])  # reduce age categories to six
    N_vacc = vcat(N[1]*2/20,N[2:end])  # only 2/20 in the first age group need to be vaccinated

    if scenario == 3 || scenario == 6
        doses_daily = doses_daily .* (25/15.9) # Increase factor: Use 25 million doses instead of the 15.9 million baseline doses
    end

    if scenario == 4 || scenario == 7
        doses_daily = doses_daily .* (35/15.9) # Increase factor: Use 35 million doses instead of the 15.9 million baseline doses
    end

    dose_interval = 8*7 # number of days between 1st and second doses
    doses_daily = doses_daily/2  # not all available doses are given have to reserve one doese for each dose given

    vacc_1_data = zeros(model_days,6)  # To hold number of dose 1 vaccinations per day in each age group ## VACC DATA SHOULD BE TIME (ROWS) x AGE GROUP (COLS)
    vacc_2_data = zeros(model_days,6)  # To hold number of dose 2 vaccinations per day in each age group ## VACC DATA SHOULD BE TIME (ROWS) x AGE GROUP (COLS)

    ds = sum(doses_daily[(pre_vacc_period+1): (pre_vacc_period+dose_interval)]) #/2  # doses that wont be issued at the tail end of vaccination period, recover them by giving earlier to keep total doses

    if rapid
       ds = ds/(180 - 56)  # per day addional dose 1 doses 56 days before end of vaccination period
    else
       ds = ds/(model_days - pre_vacc_period - 56)  # per day addional dose 1 doses 56 days before end of vaccination period
    end

    for time in (pre_vacc_period+1): (pre_vacc_period+dose_interval)
      vacc_1_data[time,:] .= (doses_daily[time] + ds).*age_dist_vacc
    end

    for time in (pre_vacc_period+dose_interval+1): model_days
       if rapid
          if  time <= (pre_vacc_period+180 - 56)
              vacc_1_data[time,:] .= (doses_daily[time] + ds).*age_dist_vacc
          end
       else
          if  time <= (model_days - 56)
              vacc_1_data[time,:] .= (doses_daily[time] + ds).*age_dist_vacc
          end
       end
       vacc_2_data[time,:] .= sum(vacc_1_data[time .- dose_interval,:]).*age_dist_vacc
    end

    # limit the number of doses to not exceed population numbers

    vacc_1_extra = KenyaCoVaccines.extra_doses(vacc_1_data,N_vacc)
    vacc_1_data_ed = vacc_1_data

    age = 6
    while age > 1  # push extra doses down the age groups
        vacc_1_data_ed[:,(age-1)] .+= vacc_1_extra[:,age]
        vacc_1_data_ed[:,age] .-= vacc_1_extra[:,age]
        vacc_1_extra = KenyaCoVaccines.extra_doses(vacc_1_data_ed,N_vacc)
        age = age -1
    end

    for age in 1:5
        vacc_1_data_ed[:,(age+1)] .+= vacc_1_extra[:,age]
        vacc_1_data_ed[:,age] .-= vacc_1_extra[:,age]
        vacc_1_extra = KenyaCoVaccines.extra_doses(vacc_1_data_ed,N_vacc)
    end

    # now check if any extra doses and set to zero (overalocation in the county) assume not transfer between counties ATM
    vacc_1_extra = KenyaCoVaccines.extra_doses(vacc_1_data_ed,N_vacc)
    vacc_1_data_ed .-= vacc_1_extra
    vacc_1_data = vacc_1_data_ed

    vacc_2_extra = KenyaCoVaccines.extra_doses(vacc_2_data,N_vacc)
    vacc_2_data_ed = vacc_2_data

    age = 6
    while age > 1  # push extra doses down the age groups
        vacc_2_data_ed[:,(age-1)] .+= vacc_2_extra[:,age]
        vacc_2_data_ed[:,age] .-= vacc_2_extra[:,age]
        vacc_2_extra = KenyaCoVaccines.extra_doses(vacc_2_data_ed,N_vacc)
        age = age -1
    end

    for age in 1:5
        vacc_2_data_ed[:,(age+1)] .+= vacc_2_extra[:,age]
        vacc_2_data_ed[:,age] .-= vacc_2_extra[:,age]
        vacc_2_extra = KenyaCoVaccines.extra_doses(vacc_2_data_ed,N_vacc)
    end

    # now check if any extra doses and set to zero (overalocation in the county) assume not transfer between counties ATM
    vacc_2_extra = KenyaCoVaccines.extra_doses(vacc_2_data_ed,N_vacc)
    vacc_2_data_ed .-= vacc_2_extra
    vacc_2_data = vacc_2_data_ed

    return vacc_1_data,vacc_2_data
end

"""
  Function to generate target age distribution of vaccinees in a given county
"""

function age_dist(county,total_doses,N_kenya,model_days,pre_vacc_period,rapid)

    pop_county =  N_kenya[:,county]
    pop_kenya =  sum(N_kenya, dims=2)

    pop_kenya_over_18 = sum(vcat((2/5)*pop_kenya[4],pop_kenya[5:end]))  # 2/5 is age 18 and 19 among 15-19 age category
    pop_county_over_18 = sum(vcat((2/5)*pop_county[4],pop_county[5:end]))

    vacc_dist_county = pop_county_over_18/pop_kenya_over_18 # proportion od vaccine to a county is proportional to its population above 18 years: assumption

    doses_per_day = zeros(model_days)
    doses_per_day[(pre_vacc_period+1):model_days]  .= total_doses/(model_days-pre_vacc_period) # use vaccination rate (flat rate, using 15.91 million doses planned for Phase I,II and III
    doses_per_day_county = doses_per_day.*vacc_dist_county #

    if rapid
        doses_per_day = zeros(model_days) # to simulate some rapid vaccination over two months
        doses_per_day[(pre_vacc_period+1):(pre_vacc_period+180)] .= total_doses/180 # compress all doses assigned to over six months
        doses_per_day_county = doses_per_day.*vacc_dist_county
    end

    N50  = sum(Array(N_kenya[:,county])[11:17])  # number >50years

    number_doses_over_50 = N50/(model_days - pre_vacc_period - 56)  # The number of doses (dose 1 ) per day to cover 100% of this age group
    prop_over_50 = min(1,number_doses_over_50/(doses_per_day_county[pre_vacc_period+1]/2)) # prioritise vaccination on 50+, this is the fraction of doses that will go to 50+

    if rapid
        number_doses_over_50 = N50/(180 - 56)  # The number of doses (dose 1 ) per day to cover 100% of this age group
        prop_over_50 = min(1,number_doses_over_50/(doses_per_day_county[pre_vacc_period+1]/2)) # prioritise vaccination on 50+, this is the fraction of doses that will go to 50+
    end

    age_dist_vacc = zeros(6)  # age distrbution of vaccinees
    age_dist_vacc[1] =  (2/32)*(1-prop_over_50)  # make 18-50 age group represent remaining doses. min(1,#) above is used because there are some counties wehre doses assigned is not able to cover beyond the pop aged >50: two age years 18 and 19 out of 32 vaccinated
    age_dist_vacc[2] = (30/32)*(1-prop_over_50)
    age_dist_vacc[3:6] .= prop_over_50/4
    sum(age_dist_vacc)

    return age_dist_vacc,doses_per_day_county
end

# Vacc predictions functions

"""

Predict and plot out vaccination impact in a county given initial condition and vaccination parameter inputs
"""

function predict_vacc_county_new(name::String,linelist_data_with_pos_neg,serology_data,
                    M_Kenya_ho,M_Kenya_other,M_Kenya_school,M_Kenya_work,N_kenya;
                    PCR_array,relative_testing_rate,
                    rapid, scenario,
                    prob_sus_w,
                    prob_sus,
                    prob_symp,
                    H_rate_by_age,
                    p_ID,p_IH,
                    deaths_data,
                    priors,
                    T_to_full_efficacy,
                    vacc_end_date = Date(2022,6,30), # date vaccination program ends
                    startdate = Date(2020,12,31),
                    enddate = Date(2022,6,30),
                    obs_factor)

    # load all_IFR_fits,ICU_fit,hosp_fit here, or pass in to the function as arguments

    # Load the model object for the county
    fit_dict = load("modelfits/$(name)_model.jld2")
    model = fit_dict[first(keys(fit_dict))]
    up_county = uppercase(model.areaname)

    #Set times, all times after model_start_time are in reference to the model_start_time
    model_start_time = model.model_start_time # start in 1st Dec 2020
    model_end_time = model.model_end_time   # end of model run for fitting
    marchendpoint = model.marchendpoint # 10th March when roudn 3 serosurvey data ends
    alpha_variant_time = model.alpha_variant_time  # earliest time of alpha variant introduction
    delta_variant_time = model.delta_variant_time # earliest time of delta variant introduction
    janstarttime = model.janstarttime  # time t start including observations in likelihood
    vacc_end_point =  (vacc_end_date - model.startdate).value # End of model run
    pre_vacc_period = (Date(2021,3,6) - model.startdate).value # vaccination starts on 6th March 2021
    model_end_point =  (enddate - startdate).value # End of model run

    # Variants data
    kenya_variants = CSV.File("data/kenya_gisaid_variants.csv") |> DataFrame
    kenya_variants.date = [Date(d,DateFormat("dd/mm/yyyy")) for d in kenya_variants.date]
    kenya_variants.rel_time = [(d - Date(2020,2,20)).value for d in kenya_variants.date]

    idxs_alpha = kenya_variants.variant .== "Alpha"
    idxs_beta = kenya_variants.variant .== "Beta"
    idxs_delta = kenya_variants.variant .== "Delta"

    #Daily prediction of variant prevalence using linear interpolation
    interp_alpha = DataInterpolations.LinearInterpolation(kenya_variants.perc_sequences[idxs_alpha]./100,kenya_variants.rel_time[idxs_alpha])
    alpha_perc = [interp_alpha(t) for t in 1:maximum(kenya_variants.rel_time[idxs_alpha])]
    interp_beta = DataInterpolations.LinearInterpolation(kenya_variants.perc_sequences[idxs_beta]./100,kenya_variants.rel_time[idxs_beta])
    beta_perc = [interp_beta(t) for t in 1:maximum(kenya_variants.rel_time[idxs_beta])]
    interp_delta = DataInterpolations.LinearInterpolation(kenya_variants.perc_sequences[idxs_delta]./100,kenya_variants.rel_time[idxs_delta])
    delta_perc = [interp_delta(t) for t in 1:maximum(kenya_variants.rel_time[idxs_delta])]

    delta_perc = delta_perc[model_start_time:end]  # pick data from 1st Dec 2020
    alpha_perc = alpha_perc[model_start_time:end]
    beta_perc  = beta_perc[model_start_time:end]

    n = model_end_point

    delta_perc = [delta_perc;ones(n-length(delta_perc))]
    alpha_beta_perc = alpha_perc .+ beta_perc
    alpha_beta_perc = [alpha_beta_perc;zeros(n-length(alpha_beta_perc))]
    other_perc = 1.0 .- alpha_beta_perc .- delta_perc

    # fitted parameters for county
    MCMCdraws = get(model.MCMC_results.chain,[:β₀,:β_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])
    prior = priors[countynames.==model.areaname][1]

    # Data sets to populate and save after loop
    crit_dis_inc_samples  = zeros(model_end_point,6,2000) # 6 age groups, 2000 samples from HMC
    severe_dis_inc_samples = zeros(model_end_point,6,2000)
    asymp_dis_inc_samples = zeros(model_end_point,6,2000)
    mild_dis_inc_samples = zeros(model_end_point,6,2000)
    death_inc_samples = zeros(model_end_point,6,2000)

    for k in 1:2000
        α= model.α  # 1/mean latent period
        ϵ= MCMCdraws.ϵ[k] #0.5    # Relative infectiousness of asymptomatic cases
        αP= model.αP # rate of movement from Pre symptomatic state
        γA = model.γA #Recovery rate from Asymptomatic state
        γM= model.γM  #Recovery rate from Mild disease

        relative_testing_rate = KenyaCoVaccines.relative_testing_rate
        PCR_array= KenyaCoVaccines.PCR_array #Relative Sensitivity of PCR test on days post-infection
        PCR_sensitivity = model.PCR_sensitivity #Base sensitivity for RT-PCR test
        PCR_specificity = model.PCR_specificity #Base specificity for RT_PCR test
        baseline_sero_array = KenyaCoVaccines.rel_sero_array_26days #Relative sensitivity of serology test on days post-infection
        sero_sensitivity = model.sero_sensitivity#Base sensitivity for serology test
        sero_specificity= model.sero_specificity #Base specificity for serology test
        M_BB= model.M_BB

        ω= model.ω # Waning of natural protection, rate per day
        σω= model.p_sus_w # relative suceptibility after first episode: 0.16 for all ages for now.
        σ = model.p_sus
        δ = model.p_symp
        υ = model.hosp_rate_by_age

        β₀ = MCMCdraws.β₀[k]
        β_home = MCMCdraws.β_home[k]
        β_school = MCMCdraws.β_school[k]
        β_other = MCMCdraws.β_other[k]
        β_work = MCMCdraws.β_work[k]

        #Set variance scalers
        clustering_factor_PCR = 0.5
        M_PCR = 30.
        T = eltype(β₀)

        inc_R_αvar= MCMCdraws.inc_R_αvar[k]
        time_scale_αvar= MCMCdraws.time_scale_αvar[k]
        mid_point_αvar= MCMCdraws.mid_point_αvar[k]
        inc_R_δvar= MCMCdraws.inc_R_δvar[k]
        time_scale_δvar= MCMCdraws.time_scale_δvar[k]
        mid_point_δvar= MCMCdraws.mid_point_δvar[k]
        init_scale = MCMCdraws.init_scale[k]

        p = convert.(T,[β₀,β_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])  # start with schools closed
        p = vcat(p,model.VE_acquisition,model.VE_infectiousness,model.VE_severe_disease_risk) # add vaccine efficacies to the paramter vector

        init_sero = model.init_sero
        γV = γM

        N= model.N
        E₀ = MCMCdraws.E₀[k]

        θ = (β₀=β₀,β_home=β_home,β_school=β_school,β_other=β_other,β_work=β_work,E₀=E₀,ϵ=ϵ,init_scale=init_scale)
        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)

        immune_escape_variant_time = Float64((Date(2022,1,1) - model.startdate).value) #  December 1 2021, time to implement variant immunity escape

        function immune_escape_variant!(integrator)
            # Reduce vaccine efficacies by 50%
            integrator.p[19:27] *= 0.5

            # Start with 'variant_E' exposed number of individuals
            variant_E = 1
            prop_vacc = sum(integrator.u[:,1:8,:],dims=2)[:,1,:] ./ sum(integrator.u[:,1:8,:],dims=[2,3])[:]
            integrator.u[:,2,:] .= variant_E .* prop_vacc # number exposed in each age group and vaccinatin status
            integrator.u[:,1,:] .= integrator.u[:,1,:] .- integrator.u[:,2,:] #substract E from S

            # move 50% from R and W -> S
            loss_susceptible_R = 0.5 .* integrator.u[:,7,:]
            loss_susceptible_W = 0.5 .* integrator.u[:,8,:]
            integrator.u[:,1,:] .+= loss_susceptible_R
            integrator.u[:,7,:] .-= loss_susceptible_R
            integrator.u[:,1,:] .+= loss_susceptible_W
            integrator.u[:,8,:] .-= loss_susceptible_W
        end
        immune_escape_variant_cb = PresetTimeCallback([immune_escape_variant_time],immune_escape_variant!,save_positions=(false,false))

        # Update callbacks based on scenario
        if  scenario==1 || scenario>4  # For Scenario 4-6
         cb = CallbackSet(cb, immune_escape_variant_cb)
        end

        county_model = KenyaCoVaccines.CoVAreaModel(name,KenyaCoVaccines.ll_onegroup_newvariant_infboost,prior;
                                PCR_cases = linelist_data_with_pos_neg,
                                sero_cases = serology_data,
                                average_sero_init = zeros(6),
                                deaths = deaths_data,
                                pop_data= N_kenya,
                                M_county_ho=M_Kenya_ho,
                                M_county_other=M_Kenya_other,
                                M_county_school=M_Kenya_school,
                                M_county_work=M_Kenya_work,
                                rapid = rapid,
                                scenario= scenario,
                                startdate=startdate,
                                enddate=enddate)

        #Solve dynamical system
        sol = solve(county_model.prob, BS3();tspan = (0,model_end_point),callback = cb,u0=u0, p=p,saveat = 1) #isoutofdomain=(u,p,t)->any(x->x<0,u)

        # Number of unscalled (for testing) cases over time by severity
        ι = KenyaCoVaccines.get_incidence_time_array_by_age_dose(sol,10)

        # Number of tests over time by age
        PCR_cases = sum(linelist_data_with_pos_neg.cases[:,linelist_data_with_pos_neg.areas .== up_county,:,:],dims = [2])[:,1,:,:]
        p_vec = [[1,2,3,4],[5,6,7,8,9,10],[11,12],[13,14],[15,16],[17]]
        PCR_cases_collaped =  zeros(Int64,size(PCR_cases,1),6,2)
        for a in 1:6
            PCR_cases_collaped[:,a:a,:] .= sum(PCR_cases[:,p_vec[a][1]:p_vec[a][end],:], dims=2)
        end
        PCR_cases = PCR_cases_collaped
        num_tests =  PCR_cases[model_start_time:end,:,1] .+ PCR_cases[model_start_time:end,:,2]
        num_tests_age = zeros(size(ι,1),6)
        for a in 1:6
            num_tests_age[:,a] = vcat(num_tests[:,a],fill(mean(num_tests[:,a]),(length(ι)-length(num_tests[:,a]))))  # use mean number of teste for forcast periods, could be changed
        end

        # scalling parameters
        p_test = MCMCdraws.p_test[k]
        χ = MCMCdraws.χ[k]
        p_test_scale = 5e-4

        # Apply vaccine efficacy against disease
        v = (1 .- model.VE_severe_disease_risk') .* υ  # efficacy

        l_inc = [inc .* v for inc in ι]
        cases_wt = similar(ι)
        cases_αβ = similar(ι)
        cases_δ  = similar(ι)
        crit_incidence = similar(ι)
        sev_incidence = similar(ι)

        for i in 1:size(ι,1)
            cases_wt[i] = sum(l_inc[i] .* other_perc[i],dims=2) # sum over dose
            cases_αβ[i] = sum(l_inc[i] .* alpha_beta_perc[i], dims=2)
            cases_δ[i]  = sum(l_inc[i] .* delta_perc[i], dims=2)

            PCR_num_wt = cases_wt[i] .* p_test_scale .* p_test
            PCR_denom_wt = PCR_num_wt .+ p_test_scale*((p_test/χ)*(N .- cases_wt[i]))
            p_PCR_pred_wt = PCR_num_wt./PCR_denom_wt
            cases_wt[i]  =  p_PCR_pred_wt .* num_tests_age[i,:]

            PCR_num_αβ = cases_αβ[i] .* p_test_scale .* p_test
            PCR_denom_αβ = PCR_num_αβ .+ p_test_scale*((p_test/χ)*(N .- cases_αβ[i]))
            p_PCR_pred_αβ = PCR_num_αβ./PCR_denom_αβ
            cases_αβ[i]  =  p_PCR_pred_αβ .* num_tests_age[i,:]

            PCR_num_δ = cases_δ[i] .* p_test_scale .* p_test
            PCR_denom_δ = PCR_num_wt .+ p_test_scale*((p_test/χ)*(N .- cases_δ[i]))
            p_PCR_pred_δ = PCR_num_δ./PCR_denom_δ
            cases_δ[i]  =  p_PCR_pred_δ .* num_tests_age[i,:]

            # need to edit structure on ICU and hosp fit to include kth sample fit rathert than point est (single value)
            crit_incidence[i] = ICU_fit[1]*(cases_wt[i] .+ ICU_fit[2].*cases_αβ[i] .+ ICU_fit[3].*cases_δ[i])
            sev_incidence[i] = hosp_fit[1]*(cases_wt[i] .+ hosp_fit[2].*cases_αβ[i] .+ hosp_fit[3].*cases_δ[i])
        end

        crit_incidence =  Matrix(VectorOfArray([inc[:] for inc in crit_incidence])')  # convert to a matrix
        sev_incidence =  Matrix(VectorOfArray([inc[:] for inc in sev_incidence])')

        for a in 1:6 # convolve with PCR array to get right timing of PCR detectable incidence
            crit_incidence[:,a] = KenyaCoVaccines.simple_conv(crit_incidence[:,a],KenyaCoVaccines.PCR_array)
            sev_incidence[:,a] = KenyaCoVaccines.simple_conv(sev_incidence[:,a],KenyaCoVaccines.PCR_array)
        end

        # Saving severe and critical disease incidence
        crit_dis_inc_samples[:,:,k]   .= crit_incidence
        severe_dis_inc_samples[:,:,k] .= sev_incidence

        # Asymptomatic cases
        asymp_dis_inc = [inc .* (ones(3)'.* (1 .- δ )) for inc in ι]
        asymp_dis_inc = Matrix(VectorOfArray([sum(inc,dims=2)[:] for inc in asymp_dis_inc])')  #sum over dose
        asymp_dis_inc_samples[:,:,k] .= asymp_dis_inc

        # Mild disease
        symp_dis_inc  =  [inc .* (ones(3)'.* δ) for inc in ι]
        symp_dis_inc =  Matrix(VectorOfArray([sum(inc,dims=2)[:] for inc in symp_dis_inc])')  #sum over dose
        for a in 1:6
            symp_dis_inc[:,a] =  KenyaCoVaccines.simple_conv(symp_dis_inc[:,a],p_IH)  # convolve with time to hsop function to get right timing of Hospitalisation
        end
        mild_dis_inc = symp_dis_inc .- sev_incidence .- crit_incidence
        mild_dis_inc_samples[:,:,k] .= mild_dis_inc

        # Number of unscalled deaths by age group over time
        IFR = all_IFR_fits[[nm[1] .== name for nm in all_IFR_fits]][1].IFR # EDIT to include HMC samples, so there shoudl be a 'k' index on the IFR
        d = (1 .- VE_severe_disease_risk') .* IFR # efficacy multiplied by risk of death by age, assumes only those with severe disease die

        death_inc =  [inc .* d for inc in ι]
        death_inc =  Matrix(VectorOfArray([sum(inc,dims=2)[:] for inc in death_inc])')  #sum over dose
        for a in 1:6
            death_inc[:,a] =  KenyaCoVaccines.simple_conv(death_inc[:,a],p_ID)  # convolve with time to death function to get right timing of deaths
        end
        death_inc = death_inc ./ obs_factor

    end # end loop over HMC chains

    # pick data from Jan for CEA analysis
    asymp_dis_inc_samples = asymp_dis_inc_samples[janstarttime:end,:,:]
    crit_dis_inc_samples = crit_dis_inc_samples[janstarttime:end,:,:]
    severe_dis_inc_samples = severe_dis_inc_samples[janstarttime:end,:,:]
    mild_dis_inc_samples = mild_dis_inc_samples[janstarttime:end,:,:]
    death_inc_samples = death_inc_samples[janstarttime:end,:,:]

    if rapid  # becasue we now have rapid scenario we save them in different directory: I have created this filing structure in the repository
        rapid_text = "rapid"
    else
        rapid_text = "non_rapid"
    end

    @save("scenario_$(scenario)/$(rapid_text)/mild_cases/severe_cases_$(name).jld2",mild_dis_inc_samples)
    @save("scenario_$(scenario)/$(rapid_text)/severe_cases/severe_cases_$(name).jld2",severe_dis_inc_samples)
    @save("scenario_$(scenario)/$(rapid_text)/critical_cases/critical_cases_$(name).jld2",crit_dis_inc_samples)
    @save("scenario_$(scenario)/$(rapid_text)/asymp_cases/asymp_cases_$(name).jld2",asymp_dis_inc_samples)
    @save("scenario_$(scenario)/$(rapid_text)/deaths/deaths_$(name).jld2",death_inc_samples)

    # Number of vaccines administered over time, this is the same for all 2000 chains
    doses_Kenya_baseline = 0.
    if scenario==1
        doses_Kenya = 0.
    else
        doses_Kenya = (1.25 + 9.76 +4.9) .* 1000000  #Number of doses planned for Kenya over the ynext one year under a 30% cover plan
    end

    age_dist_vacc, doses_daily = KenyaCoVaccines.age_dist([name],doses_Kenya,N_kenya,vacc_end_point,pre_vacc_period,rapid)
    vacc_rate_1,vacc_rate_2 = KenyaCoVaccines.vacc_rate_data(doses_daily,vacc_end_point,pre_vacc_period,age_dist_vacc,N_kenya,name,scenario,rapid)
    num_doses = sum(vacc_rate_1, dims=2)[:] .+ sum(vacc_rate_2, dims=2)[:]
    num_doses = num_doses[janstarttime:end]
    complete_doses_age = vacc_rate_2[janstarttime:end,:]  # number with complete doses

    age_dist_vacc, doses_daily = KenyaCoVaccines.age_dist([name],doses_Kenya_baseline,N_kenya,vacc_end_point,pre_vacc_period,rapid)
    vacc_rate_1,vacc_rate_2 = KenyaCoVaccines.vacc_rate_data(doses_daily,vacc_end_point,pre_vacc_period,age_dist_vacc,N_kenya,name,baseline_scenario,rapid)
    num_doses_baseline = sum(vacc_rate_1, dims=2)[:] .+ sum(vacc_rate_2, dims=2)[:]
    num_doses_baseline = num_doses_baseline[janstarttime:end]
    complete_doses_age_baseline = vacc_rate_2[janstarttime:end,:]

    @save("scenario_$(scenario)/$(rapid_text)/doses/doses_$(name).jld2",num_doses)
    @save("scenario_$(scenario)/$(rapid_text)/doses_age/doses_$(name).jld2",complete_doses_age)

    return nothing
end

"""
OLDER FUNTION THAT WAS USING POINT ESTIMATES
Predict and plot out vaccination impact in a county given initial condition and vaccination parameter inputs
"""

function predict_vacc_county(name::String,linelist_data_with_pos_neg,serology_data,
                    M_Kenya_ho,M_Kenya_other,M_Kenya_school,M_Kenya_work,N_kenya;
                    PCR_array,relative_testing_rate,
                    all_param_fits_newvariant,
                    rapid, scenario,baseline_scenario,countynames,
                    prob_sus_w,
                    prob_sus,
                    prob_symp,
                    H_rate_by_age,
                    av_sero_init,p_ID,p_IH, all_IFR_fits,ICU_fit,hosp_fit,
                    deaths_data,
                    priors,
                    T_to_full_efficacy,
                    VE_acquisition,  #Vaccine efficacy against acquisition by dose
                    VE_infectiousness, #Vaccine efficacy against transmission by dose
                    VE_severe_disease_risk,  #Vaccine efficacy against disease by dose
                    vacc_end_date = Date(2022,6,30), # date vaccination program ends
                    startdate = Date(2020,12,31),
                    enddate = Date(2022,6,30),
                    obs_factor)

    up_county = uppercase(name)

    # fitted parameters for county
    fit_p = all_param_fits_newvariant[[fit[1] .== name for fit in all_param_fits_newvariant]][1][2]
    prior = priors[countynames.==name][1]

    α= 1/3.1  # 1/mean latent period
    ϵ= fit_p[6] #0.5    # Relative infectiousness of asymptomatic cases
    αP= 1/1.9  # rate of movement from Pre symptomatic state
    γA = 1/2.4 #Recovery rate from Asymptomatic state
    γM= 1/2.4 #Recovery rate from Mild disease

    relative_testing_rate = KenyaCoVaccines.relative_testing_rate
    PCR_array= KenyaCoVaccines.PCR_array #Relative Sensitivity of PCR test on days post-infection
    PCR_sensitivity = 1. #Base sensitivity for RT-PCR test
    PCR_specificity = 0.995 #Base specificity for RT_PCR test
    baseline_sero_array = KenyaCoVaccines.rel_sero_array_26days #Relative sensitivity of serology test on days post-infection
    sero_sensitivity = 0.825 #Base sensitivity for serology test
    sero_specificity= 0.992 #Base specificity for serology test
    M_BB= 40.74

    ω= 1/180 # Waning of natural protection, rate per day
    σω= prob_sus_w  # relative suceptibility after first episode: 0.16 for all ages for now.
    σ = prob_sus
    δ = prob_symp
    υ = H_rate_by_age

    β₀ = fit_p[1]
    β_home = fit_p[2]
    β_school = fit_p[3]
    β_other = fit_p[4]
    β_work = fit_p[5]

    #Set variance scalers
    clustering_factor_PCR = 0.5
    M_PCR = 30.
    T = eltype(β₀)

    inc_R_αvar= fit_p[15]
    time_scale_αvar= fit_p[16]
    mid_point_αvar= fit_p[17]
    inc_R_δvar= fit_p[18]
    time_scale_δvar= fit_p[19]
    mid_point_δvar= fit_p[20]
    init_scale = fit_p[21]

    #Set times, all times after model_start_time are in reference to the model_start_time
    model_start_time = (startdate - Date(2020,2,20)).value # start in 1st Dec 2020
    model_end_time = (Date(2021,9,24) - startdate).value  # end of model run for fitting
    marchendpoint = (Date(2021,3,10) - startdate).value    # 10th March when roudn 3 serosurvey data ends
    alpha_variant_time = (Date(2021,2,1) - startdate).value  # earliest time of alpha variant introduction
    delta_variant_time = (Date(2021,4,1) - startdate).value  # earliest time of delta variant introduction
    janstarttime = (Date(2021,1,1) - startdate).value # time t start including observations in likelihood
    vacc_end_point =  (vacc_end_date - startdate).value # End of model run
    pre_vacc_period = (Date(2021,3,6) - startdate).value # vaccination starts on 6th March 2021
    model_end_point =  (enddate - startdate).value # End of model run

    p = convert.(T,[β₀,β_home,β_school,β_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])
    p = vcat(p,VE_acquisition,VE_infectiousness,VE_severe_disease_risk) # add vaccine efficacies to the paramter vector

    sero_cases = sum(serology_data.sero[:,serology_data.areas .== up_county,:,:],dims = [2])[:,1,:,:]
    if sum(sero_cases[(model_start_time:(model_start_time+marchendpoint)),:,:]) == 0
       init_sero = av_sero_init  # replace with the average
    else
        p_vec = [[1,2,3,4],[5,6,7,8,9,10],[11,12],[13,14],[15,16],[17]]  # reduce sero cases to 6 age groups
        sero_cases_collaped =  zeros(Int64,size(sero_cases,1),6,2)
        for a in 1:6
            sero_cases_collaped[:,a:a,:] .= sum(sero_cases[:,p_vec[a][1]:p_vec[a][end],:], dims=2)
        end
        sero_cases = sero_cases_collaped
        u_sero = sero_cases[(model_start_time:(model_start_time+marchendpoint)),:,:]
        u_sero = sum(u_sero, dims=1)[1,:,:]
        init_sero = u_sero[:,1] ./ sum(u_sero,dims=2)[:]
    end

    if isnan(init_sero[1])  # action if first value is isnan: replace with next succeding available value
        a = 1
        while isnan(init_sero[1])
            init_sero[1] = init_sero[a+1]
            a += 1
        end
    end

    for a in 2:6
        if isnan(init_sero[a])
            init_sero[a] = init_sero[a-1] # use the seroprevalence in the preceding age groups
        end
    end

    γV = γM

    N= Array(N_kenya[:,name])  #Area's population size
    N = vcat(sum(N[1:4]),sum(N[5:10]),sum(N[11:12]),sum(N[13:14]),sum(N[15:16]),N[17]) # reduce to six age categories

    E_vec = fit_p[9:14]
    u0 = KenyaCoVaccines.create_initial_conditions(init_sero,init_scale,E_vec,N,α,ω,δ,αP,υ,γA,γM,γV) #NOTE_THAT γV = γM

    #Set up the callback for lockdown
    schoolsclosed_time = Float64((Date(2021,1,1) - startdate).value)
    function affect_schoolsclosed!(integrator)
        integrator.p[3] = 0. #Schools shut
    end
    schoolsclosed_cb = PresetTimeCallback([schoolsclosed_time],affect_schoolsclosed!,save_positions=(false,false))

    #Set up callback for reopening Schools
    schoolsreopening_time = Float64((Date(2021,1,6) - startdate).value)
    function affect_schoolsopening!(integrator)
        integrator.p[3] = β_school # Schools open #
    end
    schoolsreopening_cb = PresetTimeCallback([schoolsreopening_time],affect_schoolsopening!,save_positions=(false,false))

    cb = CallbackSet(schoolsclosed_cb,schoolsreopening_cb)
    cb_baseline  = CallbackSet(schoolsclosed_cb,schoolsreopening_cb)

    immune_escape_variant_time = Float64((Date(2021,12,1) - startdate).value) #  December 1 2021, time to implement variant immunity escape

    function immune_escape_variant!(integrator)
        # Reduce vaccine efficacies by 50%
        integrator.p[19:27] *= 0.5

        # Start with 'variant_E' exposed number of individuals
        variant_E = 1
        prop_vacc = sum(integrator.u[:,1:8,:],dims=2)[:,1,:] ./ sum(integrator.u[:,1:8,:],dims=[2,3])[:]
        integrator.u[:,2,:] .= variant_E .* prop_vacc # number exposed in each age group and vaccinatin status
        integrator.u[:,1,:] .= integrator.u[:,1,:] .- integrator.u[:,2,:] #substract E from S

        # move 50% from R and W -> S
        loss_susceptible_R = 0.5 .* integrator.u[:,7,:]
        loss_susceptible_W = 0.5 .* integrator.u[:,8,:]
        integrator.u[:,1,:] .+= loss_susceptible_R
        integrator.u[:,7,:] .-= loss_susceptible_R
        integrator.u[:,1,:] .+= loss_susceptible_W
        integrator.u[:,8,:] .-= loss_susceptible_W
    end
    immune_escape_variant_cb = PresetTimeCallback([immune_escape_variant_time],immune_escape_variant!,save_positions=(false,false))

    # Update callbacks based on scenario
    if  scenario==1 || scenario>4  # For Scenario 4-6
     cb = CallbackSet(cb, immune_escape_variant_cb)
    end

    if  baseline_scenario==1 || baseline_scenario>4  # For Scenario 4-6
     cb_baseline = CallbackSet(cb_baseline, immune_escape_variant_cb)
    end

    # Vaccination rates and vaccine uptake
    county_model_baseline = KenyaCoVaccines.CoVAreaModel_proj(name,KenyaCoVaccines.ll_onegroup_newvariant_infboost,prior;
                            PCR_cases = linelist_data_with_pos_neg,
                            sero_cases = serology_data,
                            average_sero_init = av_sero_init,
                            deaths = deaths_data,
                            pop_data= N_kenya,
                            M_county_ho=M_Kenya_ho,
                            M_county_other=M_Kenya_other,
                            M_county_school=M_Kenya_school,
                            M_county_work=M_Kenya_work,
                            rapid = rapid,
                            scenario= baseline_scenario,
                            time_to_full_efficacy = T_to_full_efficacy,
                            hosp_rate_by_age = H_rate_by_age,
                            p_symp = prob_symp,
                            p_sus = prob_sus,
                            p_sus_w=prob_sus_w,
                            startdate=startdate,
                            enddate=enddate)

    county_model = KenyaCoVaccines.CoVAreaModel_proj(name,KenyaCoVaccines.ll_onegroup_newvariant_infboost,prior;
                            PCR_cases = linelist_data_with_pos_neg,
                            sero_cases = serology_data,
                            average_sero_init = av_sero_init,
                            deaths = deaths_data,
                            pop_data= N_kenya,
                            M_county_ho=M_Kenya_ho,
                            M_county_other=M_Kenya_other,
                            M_county_school=M_Kenya_school,
                            M_county_work=M_Kenya_work,
                            rapid = rapid,
                            scenario= scenario,
                            time_to_full_efficacy = T_to_full_efficacy,
                            hosp_rate_by_age = H_rate_by_age,
                            p_symp = prob_symp,
                            p_sus = prob_sus,
                            p_sus_w=prob_sus_w,
                            startdate=startdate,
                            enddate=enddate)

    #Solve dynamical system
    sol = solve(county_model.prob, BS3();tspan = (0,model_end_point),callback = cb,u0=u0, p=p,saveat = 1) #isoutofdomain=(u,p,t)->any(x->x<0,u)
    sol_baseline = solve(county_model_baseline.prob, BS3();tspan = (0,model_end_point),callback = cb_baseline,u0=u0, p=p, saveat = 1) # isoutofdomain=(u,p,t)->any(x->x<0,u)

    # Number of unscalled (for testing) cases over time by severity
    ι = KenyaCoVaccines.get_incidence_time_array_by_age_dose(sol,10)
    ι_baseline = KenyaCoVaccines.get_incidence_time_array_by_age_dose(sol_baseline,10)

    # Variants data
    kenya_variants = CSV.File("data/kenya_gisaid_variants.csv") |> DataFrame
    kenya_variants.date = [Date(d,DateFormat("dd/mm/yyyy")) for d in kenya_variants.date]
    kenya_variants.rel_time = [(d - Date(2020,2,20)).value for d in kenya_variants.date]

    idxs_alpha = kenya_variants.variant .== "Alpha"
    idxs_beta = kenya_variants.variant .== "Beta"
    idxs_delta = kenya_variants.variant .== "Delta"

    # Number of tests over time by age
    PCR_cases = sum(linelist_data_with_pos_neg.cases[:,linelist_data_with_pos_neg.areas .== up_county,:,:],dims = [2])[:,1,:,:]
    p_vec = [[1,2,3,4],[5,6,7,8,9,10],[11,12],[13,14],[15,16],[17]]
    PCR_cases_collaped =  zeros(Int64,size(PCR_cases,1),6,2)
    for a in 1:6
        PCR_cases_collaped[:,a:a,:] .= sum(PCR_cases[:,p_vec[a][1]:p_vec[a][end],:], dims=2)
    end
    PCR_cases = PCR_cases_collaped
    num_tests =  PCR_cases[model_start_time:end,:,1] .+ PCR_cases[model_start_time:end,:,2]
    num_tests_age = zeros(size(ι,1),6)
    for a in 1:6
        num_tests_age[:,a] = vcat(num_tests[:,a],fill(mean(num_tests[:,a]),(length(ι)-length(num_tests[:,a]))))  # use mean number of teste for forcast periods
    end

    #Daily prediction of variant prevalence using linear interpolation
    interp_alpha = DataInterpolations.LinearInterpolation(kenya_variants.perc_sequences[idxs_alpha]./100,kenya_variants.rel_time[idxs_alpha])
    alpha_perc = [interp_alpha(t) for t in 1:maximum(kenya_variants.rel_time[idxs_alpha])]
    interp_beta = DataInterpolations.LinearInterpolation(kenya_variants.perc_sequences[idxs_beta]./100,kenya_variants.rel_time[idxs_beta])
    beta_perc = [interp_beta(t) for t in 1:maximum(kenya_variants.rel_time[idxs_beta])]
    interp_delta = DataInterpolations.LinearInterpolation(kenya_variants.perc_sequences[idxs_delta]./100,kenya_variants.rel_time[idxs_delta])
    delta_perc = [interp_delta(t) for t in 1:maximum(kenya_variants.rel_time[idxs_delta])]

    delta_perc = delta_perc[model_start_time:end]  # pick data from 1st Dec 2020
    alpha_perc = alpha_perc[model_start_time:end]
    beta_perc  = beta_perc[model_start_time:end]

    n = size(ι ,1)

    delta_perc = [delta_perc;ones(n-length(delta_perc))]
    alpha_beta_perc = alpha_perc .+ beta_perc
    alpha_beta_perc = [alpha_beta_perc;zeros(n-length(alpha_beta_perc))]
    other_perc = 1.0 .- alpha_beta_perc .- delta_perc

    # scalling parameters
    p_test = fit_p[8]
    χ = fit_p[7]
    p_test_scale = 5e-4

    # Apply vaccine efficacy against disease
    υ = (1 .- VE_severe_disease_risk') .* ones(6)  # efficacy
    if baseline_scenario==0
        υ_baseline = ones(3)' .* ones(6) # No vaccine efficacy
    else
        υ_baseline = (1 .- VE_severe_disease_risk') .* ones(6)  # in case of a baselise scenarion with vaccination
    end

    l_inc = [inc .* υ for inc in ι]
    cases_wt = similar(l_inc)
    cases_αβ = similar(l_inc)
    cases_δ  = similar(l_inc)
    crit_incidence = similar(l_inc)
    sev_incidence = similar(l_inc)

    l_baseline_inc = [inc .* υ_baseline for inc in ι_baseline]
    cases_baseline_wt = similar(l_baseline_inc)
    cases_baseline_αβ = similar(l_baseline_inc)
    cases_baseline_δ  = similar(l_baseline_inc)
    crit_incidence_baseline = similar(l_baseline_inc)
    sev_incidence_baseline = similar(l_baseline_inc)

    for i in 1:size(l_inc,1)
        cases_wt[i] = sum(l_inc[i] .* other_perc[i],dims=2) # sum over dose
        cases_αβ[i] = sum(l_inc[i] .* alpha_beta_perc[i], dims=2)
        cases_δ[i]  = sum(l_inc[i] .* delta_perc[i], dims=2)

        PCR_num_wt = cases_wt[i] .* p_test_scale .* p_test
        PCR_denom_wt = PCR_num_wt .+ p_test_scale*((p_test/χ)*(N .- cases_wt[i]))
        p_PCR_pred_wt = PCR_num_wt./PCR_denom_wt
        cases_wt[i]  =  p_PCR_pred_wt .* num_tests_age[i,:]

        PCR_num_αβ = cases_αβ[i] .* p_test_scale .* p_test
        PCR_denom_αβ = PCR_num_αβ .+ p_test_scale*((p_test/χ)*(N .- cases_αβ[i]))
        p_PCR_pred_αβ = PCR_num_αβ./PCR_denom_αβ
        cases_αβ[i]  =  p_PCR_pred_αβ .* num_tests_age[i,:]

        PCR_num_δ = cases_δ[i] .* p_test_scale .* p_test
        PCR_denom_δ = PCR_num_wt .+ p_test_scale*((p_test/χ)*(N .- cases_δ[i]))
        p_PCR_pred_δ = PCR_num_δ./PCR_denom_δ
        cases_δ[i]  =  p_PCR_pred_δ .* num_tests_age[i,:]

        crit_incidence[i] = ICU_fit[1]*(cases_wt[i] .+ ICU_fit[2].*cases_αβ[i] .+ ICU_fit[3].*cases_δ[i])
        sev_incidence[i] = hosp_fit[1]*(cases_wt[i] .+ hosp_fit[2].*cases_αβ[i] .+ hosp_fit[3].*cases_δ[i])

        # baseline  iteration
        cases_baseline_wt[i] = sum(l_baseline_inc[i] .* other_perc[i],dims=2) # sum over dose
        cases_baseline_αβ[i] = sum(l_baseline_inc[i] .* alpha_beta_perc[i], dims=2)
        cases_baseline_δ[i]  = sum(l_baseline_inc[i] .* delta_perc[i], dims=2)

        PCR_num_wt = cases_baseline_wt[i] .* p_test_scale .* p_test
        PCR_denom_wt = PCR_num_wt .+ p_test_scale*((p_test/χ)*(N .- cases_baseline_wt[i]))
        p_PCR_pred_wt = PCR_num_wt./PCR_denom_wt
        cases_baseline_wt[i]  =  p_PCR_pred_wt .* num_tests_age[i,:]

        PCR_num_αβ = cases_baseline_αβ[i] .* p_test_scale .* p_test
        PCR_denom_αβ = PCR_num_αβ .+ p_test_scale*((p_test/χ)*(N .- cases_baseline_αβ[i]))
        p_PCR_pred_αβ = PCR_num_αβ./PCR_denom_αβ
        cases_baseline_αβ[i]  =  p_PCR_pred_αβ .* num_tests_age[i,:]

        PCR_num_δ = cases_baseline_δ[i] .* p_test_scale .* p_test
        PCR_denom_δ = PCR_num_wt .+ p_test_scale*((p_test/χ)*(N .- cases_baseline_δ[i]))
        p_PCR_pred_δ = PCR_num_δ./PCR_denom_δ
        cases_baseline_δ[i]  =  p_PCR_pred_δ .* num_tests_age[i,:]

        crit_incidence_baseline[i] = ICU_fit[1]*(cases_baseline_wt[i] .+ ICU_fit[2].*cases_baseline_αβ[i] .+ ICU_fit[3].*cases_baseline_δ[i])
        sev_incidence_baseline[i] = hosp_fit[1]*(cases_baseline_wt[i] .+ hosp_fit[2].*cases_baseline_αβ[i] .+ hosp_fit[3].*cases_baseline_δ[i])
    end

    crit_incidence =  Matrix(VectorOfArray([inc[:] for inc in crit_incidence])')  # convert to a matrix
    sev_incidence =  Matrix(VectorOfArray([inc[:] for inc in sev_incidence])')
    crit_incidence_baseline=  Matrix(VectorOfArray([inc[:] for inc in crit_incidence_baseline])')  # convert to a matrix
    sev_incidence_baseline =  Matrix(VectorOfArray([inc[:] for inc in sev_incidence_baseline])')

    for a in 1:6 # convolve with PCR array to get right timing of PCR detectable incidence
        crit_incidence[:,a] = KenyaCoVaccines.simple_conv(crit_incidence[:,a],KenyaCoVaccines.PCR_array)
        sev_incidence[:,a] = KenyaCoVaccines.simple_conv(sev_incidence[:,a],KenyaCoVaccines.PCR_array)
        crit_incidence_baseline[:,a] = KenyaCoVaccines.simple_conv(crit_incidence_baseline[:,a],KenyaCoVaccines.PCR_array)
        sev_incidence_baseline[:,a] = KenyaCoVaccines.simple_conv(sev_incidence_baseline[:,a],KenyaCoVaccines.PCR_array)
    end

    # Saving severe and critical disease incidence
    crit_dis_inc   = crit_incidence[janstarttime:end,:]
    severe_dis_inc = sev_incidence[janstarttime:end,:]
    crit_dis_inc_baseline   = crit_incidence_baseline[janstarttime:end,:]
    severe_dis_inc_baseline = sev_incidence_baseline[janstarttime:end,:]

    @save("scenario_$(scenario)/severe_cases/severe_cases_$(name).jld2",severe_dis_inc)
    @save("scenario_$(baseline_scenario)/severe_cases/severe_cases_$(name).jld2",severe_dis_inc_baseline)
    @save("scenario_$(scenario)/critical_cases/critical_cases_$(name).jld2",crit_dis_inc)
    @save("scenario_$(baseline_scenario)/critical_cases/critial_cases_$(name).jld2",crit_dis_inc_baseline)

    # Asymptomatic cases
    asymp_dis_inc = [inc .* (ones(3)'.* (1 .- prob_symp)) for inc in ι]
    asymp_dis_inc = Matrix(VectorOfArray([sum(inc,dims=2)[:] for inc in asymp_dis_inc])')  #sum over dose
    asymp_dis_inc = asymp_dis_inc[janstarttime:end,:]

    asymp_dis_inc_baseline  =  [inc .* (ones(3)'.* (1 .- prob_symp)) for inc in ι_baseline]
    asymp_dis_inc_baseline = Matrix(VectorOfArray([sum(inc,dims=2)[:] for inc in asymp_dis_inc_baseline])')  #sum over dose
    asymp_dis_inc_baseline = asymp_dis_inc_baseline[janstarttime:end,:]

    @save("scenario_$(scenario)/asymp_cases/asymp_cases_$(name).jld2",asymp_dis_inc)
    @save("scenario_$(baseline_scenario)/asymp_cases/asymp_cases_$(name).jld2",asymp_dis_inc_baseline)

    # Mild disease
    symp_dis_inc  =  [inc .* (ones(3)'.* prob_symp) for inc in ι]
    symp_dis_inc =  Matrix(VectorOfArray([sum(inc,dims=2)[:] for inc in symp_dis_inc])')  #sum over dose
    for a in 1:6
        symp_dis_inc[:,a] =  KenyaCoVaccines.simple_conv(symp_dis_inc[:,a],p_IH)  # convolve with time to hsop function to get right timing of Hospitalisation
    end
    symp_dis_inc = symp_dis_inc[janstarttime:end,:]
    mild_dis_inc = symp_dis_inc .- severe_dis_inc .- crit_dis_inc

    symp_dis_inc_baseline =  [inc .* (ones(3)'.* prob_symp) for inc in ι_baseline]
    symp_dis_inc_baseline =  Matrix(VectorOfArray([sum(inc,dims=2)[:] for inc in symp_dis_inc_baseline])')  #sum over dose
    for a in 1:6
        symp_dis_inc_baseline[:,a] =  KenyaCoVaccines.simple_conv(symp_dis_inc_baseline[:,a],p_IH)  # convolve with time to hsop function to get right timing of Hospitalisation
    end
    symp_dis_inc_baseline = symp_dis_inc_baseline[janstarttime:end,:]
    mild_dis_inc_baseline = symp_dis_inc_baseline .- severe_dis_inc_baseline .- crit_dis_inc_baseline

    @save("scenario_$(scenario)/mild_cases/mild_cases_$(name).jld2",mild_dis_inc)
    @save("scenario_$(baseline_scenario)/mild_cases/mild_cases_$(name).jld2",mild_dis_inc_baseline)

    # Number of unscalled deaths by age group over time
    IFR = all_IFR_fits[[nm[1] .== name for nm in all_IFR_fits]][1].IFR
    υ = (1 .- VE_severe_disease_risk') .* IFR # efficacy multiplied by risk of death by age, assumes only those with severe disease die
    υ_baseline = ones(3)' .* IFR # No vaccine efficacy

    death_inc =  [inc .* υ for inc in ι]
    death_inc =  Matrix(VectorOfArray([sum(inc,dims=2)[:] for inc in death_inc])')  #sum over dose
    for a in 1:6
        death_inc[:,a] =  KenyaCoVaccines.simple_conv(death_inc[:,a],p_ID)  # convolve with time to death function to get right timing of deaths
    end
    death_inc = death_inc[janstarttime:end,:] ./ obs_factor

    death_inc_baseline =  [inc .* υ for inc in ι_baseline]
    death_inc_baseline =  Matrix(VectorOfArray([sum(inc,dims=2)[:] for inc in death_inc_baseline])')  #sum over dose
    for a in 1:6
        death_inc_baseline[:,a] =  KenyaCoVaccines.simple_conv(death_inc_baseline[:,a],p_ID)  # convolve with time to death function to get right timing of deaths
    end
    death_inc_baseline = death_inc_baseline[janstarttime:end,:] ./ obs_factor

    @save("scenario_$(scenario)/deaths/deaths_$(name).jld2",death_inc)
    @save("scenario_$(baseline_scenario)/deaths/deaths_$(name).jld2",death_inc_baseline)

    # Number of vaccines administered over time
    # number of doses planned to cover about 30% of pop>18y

    doses_Kenya_baseline = 0.
    if scenario==1
        doses_Kenya = 0.
    else
        doses_Kenya = (1.25 + 9.76 +4.9) .* 1000000  #Number of doses planned for Kenya over the ynext one year under a 30% cover plan
    end

    age_dist_vacc, doses_daily = KenyaCoVaccines.age_dist([name],doses_Kenya,N_kenya,vacc_end_point,pre_vacc_period,rapid)
    vacc_rate_1,vacc_rate_2 = KenyaCoVaccines.vacc_rate_data(doses_daily,vacc_end_point,pre_vacc_period,age_dist_vacc,N_kenya,name,scenario,rapid)
    num_doses = sum(vacc_rate_1, dims=2)[:] .+ sum(vacc_rate_2, dims=2)[:]
    num_doses = num_doses[janstarttime:end]
    complete_doses_age = vacc_rate_2[janstarttime:end,:]  # number with complete doses

    age_dist_vacc, doses_daily = KenyaCoVaccines.age_dist([name],doses_Kenya_baseline,N_kenya,vacc_end_point,pre_vacc_period,rapid)
    vacc_rate_1,vacc_rate_2 = KenyaCoVaccines.vacc_rate_data(doses_daily,vacc_end_point,pre_vacc_period,age_dist_vacc,N_kenya,name,baseline_scenario,rapid)
    num_doses_baseline = sum(vacc_rate_1, dims=2)[:] .+ sum(vacc_rate_2, dims=2)[:]
    num_doses_baseline = num_doses_baseline[janstarttime:end]
    complete_doses_age_baseline = vacc_rate_2[janstarttime:end,:]

    @save("scenario_$(scenario)/doses/doses_$(name).jld2",num_doses)
    @save("scenario_$(baseline_scenario)/doses/doses_$(name).jld2",num_doses_baseline)

    @save("scenario_$(scenario)/doses_age/doses_$(name).jld2",complete_doses_age)
    @save("scenario_$(baseline_scenario)/doses_age/doses_$(name).jld2",complete_doses_age_baseline)

    # Plot vaccination impact for county

    ι = KenyaCoVaccines.get_incidence_time_array(sol,10)
    ι_baseline = KenyaCoVaccines.get_incidence_time_array(sol_baseline,10)
    p_test = fit_p[8]
    χ = fit_p[7]
    p_test_scale = 5e-4

    PCR = KenyaCoVaccines.simple_conv(ι,PCR_array)
    PCR_cases = sum(linelist_data_with_pos_neg.cases[:,linelist_data_with_pos_neg.areas .== up_county,:,:],dims = [2,3])[:,1,1,:]
    PCR_cases_pos = PCR_cases[model_start_time:end,1]

    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
    num_tests =  PCR_cases[model_start_time:end,1] + PCR_cases[model_start_time:end,2]

    # Might later change number of tests function from mean
    num_tests = vcat(num_tests,fill(mean(num_tests),(length(p_PCR_pred)-length(num_tests))))
    pred = num_tests .* p_PCR_pred

    PCR = KenyaCoVaccines.simple_conv(ι_baseline,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
    pred_baseline = num_tests .* p_PCR_pred

    pcr_pos_ma = KenyaCoVaccines.sma(PCR_cases_pos,7)  # add 7-day movinf average to the plot

    cover = 0
    if scenario == 2 || scenario == 5
        cover =30
    end
    if scenario == 3 || scenario == 6
        cover = 50
    end
    if scenario == 4 || scenario == 7
        cover = 70
    end

    baseline_cover = 0
    if baseline_scenario == 2 || baseline_scenario == 5
        baseline_cover =30
    end
    if baseline_scenario == 3 || baseline_scenario == 6
        baseline_cover = 50
    end
    if baseline_scenario == 4 || baseline_scenario == 7
        baseline_cover = 70
    end

    xticks_monthly = [(Date(y,m,1) - Date(2021,1,1)).value for m = 1:12,y = 2021:2022][1:end]
    xticklabs_monthly = [monthname(m)[1:3]*"-"*string(y-2000) for m = 1:12,y = 2021:2022][1:end]

    plt = plot(pred[janstarttime:end],
            lab = "$(cover)% coverage",
            lw=1,
            legend = :topright,
            title = "PCR prevalence $(name)",
            xticks = (xticks_monthly,xticklabs_monthly),
            ylabel = "Daily cases",
            color=:blue,
            size = (800,600),dpi = 300)
    plot!(plt,pred_baseline[janstarttime:end],color=:red,lab = "Baseline:$(baseline_cover)% coverage",lw=1)
    plot!(plt,pcr_pos_ma[janstarttime:end],lw=1,color=:green,lab="7-day MA PCR prevalence")
    scatter!(plt,PCR_cases_pos[janstarttime:end],color=:green,lab="Observed PCR prevalence")

    savefig(plt,"scenario_$(scenario)/plots/vaccination_impact/predicted_PCR_incidence_$(name).png")

    # Save scalled PCR incidence for later plotting Kenyawide model fit
    if  baseline_scenario==0
        pred_baseline_fit =  pred_baseline[janstarttime:model_end_time]
        @save("scenario_$(baseline_scenario)/scaled_PCR_pred/scaled_PCR_pred_$(name).jld2",pred_baseline_fit)
    end

    return nothing
end

"""

Pedict vaccination impact of multiple counties using the `predict_vacc_county` method.
"""
function predict_vacc_for_multiple_counties(namelist::Vector{String},
                    linelist_data_with_pos_neg,serology_data,
                    M_Kenya_ho,M_Kenya_other,M_Kenya_school,M_Kenya_work,N_kenya;
                    PCR_array,
                    relative_testing_rate,
                    all_param_fits_newvariant,
                    rapid,
                    scenario,
                    baseline_scenario,
                    countynames,
                    prob_sus_w,
                    prob_sus,
                    prob_symp,
                    H_rate_by_age,
                    av_sero_init,
                    p_ID,
                    p_IH,
                    all_IFR_fits,
                    ICU_fit,
                    hosp_fit,
                    deaths_data,
                    priors,
                    T_to_full_efficacy,
                    VE_acquisition,  #Vaccine efficacy against acquisition by dose
                    VE_infectiousness, #Vaccine efficacy against transmission by dose
                    VE_severe_disease_risk,  #Vaccine efficacy against disease by dose
                    vacc_end_date, # date vaccination program ends
                    startdate,
                    enddate,
                    obs_factor)

    for name in namelist
        println("Starting projection for $(name)")
        pred = KenyaCoVaccines.predict_vacc_county(name,linelist_data_with_pos_neg,serology_data,
                            M_Kenya_ho,M_Kenya_other,M_Kenya_school,M_Kenya_work,N_kenya;
                            PCR_array=PCR_array,
                            relative_testing_rate=relative_testing_rate,
                            all_param_fits_newvariant = all_param_fits_newvariant,
                            rapid=rapid,
                            scenario=scenario,
                            baseline_scenario=baseline_scenario,
                            countynames=countynames,
                            prob_sus_w=prob_sus_w,
                            prob_sus=prob_sus,
                            prob_symp=prob_symp,
                            H_rate_by_age=H_rate_by_age,
                            av_sero_init=av_sero_init,
                            p_ID=p_ID,
                            p_IH=p_IH,
                            all_IFR_fits=all_IFR_fits,
                            ICU_fit=ICU_fit,
                            hosp_fit=hosp_fit,
                            deaths_data = deaths_data,
                            priors=priors,
                            T_to_full_efficacy=T_to_full_efficacy,
                            VE_acquisition=VE_acquisition,  #Vaccine efficacy against acquisition by dose
                            VE_infectiousness=VE_infectiousness, #Vaccine efficacy against transmission by dose
                            VE_severe_disease_risk=VE_severe_disease_risk,  #Vaccine efficacy against disease by dose
                            vacc_end_date = vacc_end_date, # date vaccination program ends
                            startdate = startdate,
                            enddate = enddate,
                            obs_factor=obs_factor)
    end
    return nothing
end

function sma(a::Vector, n::Int)  # moving average function
    vals = zeros(size(a,1) - (n-1))

    for i in 1:size(a,1) - (n-1)
      vals[i] = mean(a[i:i+(n-1)])
    end
    return vcat(missings(n-1),vals)
end


text:
export
    equinox_approx,
    equinox

# Calculate the times of solstice and equinox events for Earth

# Meeus-1998 Table 27.A

equinox_approx_1000 = Dict{String,Array{Float64,1}}(
    "spring" => [1721139.29189, 365242.13740,  0.06134,  0.00111, -0.00071],
    "summer" => [1721233.25401, 365241.72562, -0.05323,  0.00907, -0.00025],
    "autumn" => [1721325.70455, 365242.49558, -0.11677, -0.00297,  0.00074],
    "winter" => [1721414.39987, 365242.88257, -0.00769, -0.00933, -0.00006]
    )

# Meeus-1998 Table 27.B

equinox_approx_3000 = Dict{String,Array{Float64,1}}(
    "spring" => [2451623.80984, 365242.37404,  0.05169, -0.00411, -0.00057],
    "summer" => [2451716.56767, 365241.62603,  0.00325,  0.00888, -0.00030],
    "autumn" => [2451810.21715, 365242.01767, -0.11575,  0.00337,  0.00078],
    "winter" => [2451900.05952, 365242.74049, -0.06223, -0.00823,  0.00032]
    )

# Meeus-1998 Table 27.C

equinox_terms = Array[
    [485, deg2rad(324.96),  deg2rad(  1934.136)],
    [203, deg2rad(337.23),  deg2rad( 32964.467)],
    [199, deg2rad(342.08),  deg2rad(    20.186)],
    [182, deg2rad( 27.85),  deg2rad(445267.112)],
    [156, deg2rad( 73.14),  deg2rad( 45036.886)],
    [136, deg2rad(171.52),  deg2rad( 22518.443)],
    [ 77, deg2rad(222.54),  deg2rad( 65928.934)],
    [ 74, deg2rad(296.72),  deg2rad(  3034.906)],
    [ 70, deg2rad(243.58),  deg2rad(  9037.513)],
    [ 58, deg2rad(119.81),  deg2rad( 33718.147)],
    [ 52, deg2rad(297.17),  deg2rad(   150.678)],
    [ 50, deg2rad( 21.02),  deg2rad(  2281.226)],
    [ 45, deg2rad(247.54),  deg2rad( 29929.562)],
    [ 44, deg2rad(325.15),  deg2rad( 31555.956)],
    [ 29, deg2rad( 60.93),  deg2rad(  4443.417)],
    [ 18, deg2rad(155.12),  deg2rad( 67555.328)],
    [ 17, deg2rad(288.79),  deg2rad(  4562.452)],
    [ 16, deg2rad(198.04),  deg2rad( 62894.029)],
    [ 14, deg2rad(199.76),  deg2rad( 31436.921)],
    [ 12, deg2rad( 95.39),  deg2rad( 14577.848)],
    [ 12, deg2rad(287.11),  deg2rad( 31931.756)],
    [ 12, deg2rad(320.81),  deg2rad( 34777.259)],
    [  9, deg2rad(227.73),  deg2rad(  1222.114)],
    [  8, deg2rad( 15.45),  deg2rad( 16859.074)]
    ]

"""
    equinox_approx(yr, season)

Returns the approximate time of a solstice or equinox event.

The year must be in the range -1000...3000. Within that range the the error from
the precise instant is at most 2.16 minutes.

Parameters:
    yr     : year -1000 ... 3000
    season : one of ("spring", "summer", "autumn", "winter")

Returns:
    Julian Day of the event in dynamical time
"""
function equinox_approx(yr, season)
    if yr < -1000 || yr > 3000
        error("year is out of range")
    end
    yr = floor(Integer, yr)
    if yr > -1000 && yr <= 1000
        Y = yr / 1000.0
        tbl = equinox_approx_1000
    else
        Y = (yr - 2000) / 1000.0
        tbl = equinox_approx_3000
    end
    jd = polynomial(tbl[season], Y)
    T = jd_to_jcent(jd)
    W = deg2rad(35999.373 * T - 2.47)
    delta_lambda = 1 + 0.0334 * cos(W) + 0.0007 * cos(2 * W)
    S = 0.0
    for v in equinox_terms
        S = S + v[1] * cos(v[2] + v[3] * T)
    end
    jd += (0.00001 * S) / delta_lambda
    return jd
end

equinox_circle = Dict{AbstractString, Float64}(
    "spring" => 0.0,
    "summer" => pi * 0.5,
    "autumn" => pi,
    "winter" => pi * 1.5
    )

equinox_k_sun_motion = 365.25 / (pi * 2)

"""
    equinox(jd, season, delta)

Return the precise moment of an equinox or solstice event on Earth.

Parameters:
    jd     : Julian Day of an approximate time of the event in dynamical time
    season : one of ("spring", "summer", "autumn", "winter")
    delta  : the required precision in days. Times accurate to a second are
        reasonable when using the VSOP model.

Returns:
    Julian Day of the event in dynamical time
"""
function equinox(jd::Float64, season, delta)
    # TODO this fails
    # if we knew that the starting approximate time was close enough
    # to the actual time, we could pull nut_in_lon() and the aberration
    # out of the loop and save some calculating.

    #   calls: sun_dimension3(jd) which calls vsop87d_dimension(jd, "Earth") to get (L, B, R)
    circ = equinox_circle[lowercase(season)]

    for k in 1:20
        jd0 = jd
        L, B, R = sun_dimension3(jd)
        L += nut_in_lon(jd) + aberration_low(R)
        L, B = vsop_to_fk5(jd, L, B)
        jd += 58 * sin(circ - L)
        if abs(jd - jd0) < delta
           return jd
        end
    end
    error("bailout")
end

ext: jl


text:
using FunctionalCollections
using Test

const PAM = PersistentArrayMap

@testset "Persistent Array Maps" begin

    @testset "construction" begin
        @test length(PAM{Int, Int}().kvs) == 0
        @test length(PAM((1, 1), (2, 2)).kvs) == 2

        @test length(PAM((1, 1))) == 1
        @test length(PAM((1, 1), (2, 2))) == 2
    end

    @testset "accessing" begin
        m = PAM((1, "one"), (2, "two"), (3, "three"))
        @test m[1] == "one"
        @test m[2] == "two"
        @test m[3] == "three"

        @test get(m, 1) == "one"
        @test get(m, 1, "foo") == "one"
        @test get(m, 4, "foo") == "foo"
    end

    @testset "haskey" begin
        m = PAM((1, "one"), (2, "two"), (3, "three"))
        @test haskey(m, 1)
        @test haskey(m, 2)
        @test haskey(m, 3)
        @test !haskey(m, 4)
    end

    @testset "assoc" begin
        m = PAM{Int, String}()
        @test assoc(m, 1, "one")[1] == "one"
        @test try m[1]; false catch e true end

        m = PAM((1, "one"))
        @test assoc(m, 1, "foo")[1] == "foo"
    end

    @testset "dissoc" begin
        m = PAM((1, "one"))
        m = dissoc(m, 1)
        @test try m[1]; false catch e true end

        m = PAM((1, "one"), (2, "two"))
        @test dissoc(m, 1) == PAM((2, "two"))
    end

    @testset "iterating" begin
        m = PAM((1, "one"), (2, "two"), (3, "three"))
        @test [v for (k, v) in m] == ["one", "two", "three"]
    end

    @testset "isempty" begin
        @test isempty(PAM{Int, Int}())
        @test !isempty(PAM((1, "one")))
    end

    @testset "equality" begin
        @test PAM((1, "one")) == PAM((1, "one"))
        @test PAM((1, "one"), (2, "two")) == PAM((2, "two"), (1, "one"))
        @test isequal(PAM((1, "one")), PAM((1, "one")))

        @test PAM((1, "one")) != (PAM((2, "two")))
    end

    @testset "kwargs construction" begin
        @test PAM(x=1, y=2, z=3) == PAM((:x, 1), (:y, 2), (:z, 3))
    end

    @testset "map" begin
        m = PAM((1, 1), (2, 2), (3, 3))

        @test map((kv) -> (kv[1], kv[2]+1), m) == PAM((1, 2), (2, 3), (3, 4))
    end
end

const PHM = PersistentHashMap

@testset "Persistent Hash Maps" begin

    @testset "constructor" begin
        hashmap = PHM{Int, Int}()
        @test length(hashmap) == 0
        @test length(PHM((1, 1), (2, 2), (3, 3))) == 3
        @test length(PHM(x=1, y=2, z=3)) == 3
    end

    @testset "equality" begin
        @test PHM{Int, Int}() == PHM{Int, Int}()
        @test PHM{Int, Int}() == PHM{String, String}()

        m1 = PHM{Int, Int}()
        m2 = PHM{Int, Int}()
        @test assoc(m1, 1, 100) == assoc(m2, 1, 100)
        @test assoc(m1, 1, 100) != (assoc(m2, 1, 200))
        @test assoc(m1, 1, 100) != (assoc(m2, 2, 100))

        m3 = PHM([(1 => 10), (2 => 20), (3 => 30)])
        m4 = PHM((3, 30), (2, 20), (1, 10))
        @test m3 == m4
        @test m3 != (m1)

        @test m3 == Dict(1 => 10, 2 => 20, 3 => 30)
    end

    @testset "assoc" begin
        m = PHM{Int, String}()
        @test assoc(m, 1, "one")[1] == "one"
        @test try m[1]; false catch e true end

        m = PHM{Int, String}()
        m = assoc(m, 1, "one")
        @test assoc(m, 1, "foo")[1] == "foo"
    end

    @testset "covariance" begin
        m = PHM{Any, Any}()
        @test assoc(m, "foo", "bar") == (Dict("foo" => "bar"))
    end

    @testset "dissoc" begin
        m = PAM((1, "one"))
        m = dissoc(m, 1)
        @test try m[1]; false catch e true end

        m = PHM((1, "one"), (2, "two"))
        @test dissoc(m, 1) == PHM((2, "two"))
    end

    @testset "get" begin
        m = PHM{Int, String}()
        @test get(m, 1, "default") == "default"
        m = assoc(m, 1, "one")
        @test get(m, 1, "default") == "one"
        m = assoc(m, 1, "newone")
        @test get(m, 1, "default") == "newone"
        @test try get(m, 2); false catch e true end
    end

    @testset "haskey" begin
        m = PHM{Int, String}()
        @test !haskey(m, 1)
        m = assoc(m, 1, "one")
        @test haskey(m, 1)
        @test !haskey(m, 2)
    end

    @testset "haskey dissoc" begin
        m = PHM{Int, String}()
        m = assoc(m, 1, "one")
        m = dissoc(m, 1)
        @test !haskey(m, 1)
    end

    @testset "map" begin
        m = PHM((1, 1), (2, 2), (3, 3))
        @test map((kv) -> (kv[1], kv[2]+1), m) == PHM((1, 2), (2, 3), (3, 4))
    end

    @testset "filter" begin
        @test filter((kv) -> iseven(kv[2]), PHM((1, 1), (2, 2))) == PHM((2, 2))
    end

    @testset "merge" begin
        @test merge(PHM((1, 1), (2, 2)), PHM((2, 3), (3, 4))) ==
                PHM((1,1), (2, 3), (3, 4))
    end
end


text:
struct Input

    y

    Z
    T
    d
    c
    R
end

struct Options
    constant_in_time::Bool
end

Base.@kwdef mutable struct Sizes
    n::Integer = 0
    p::Integer = 0
    m::Integer = 0
    r::Integer = 0

end

struct Problem
    input::Input
    options::Options
    sizes::Sizes
end

Base.@kwdef mutable struct MidTerm

end

Base.@kwdef mutable struct Output
end

max_stars_repo_path: Modelos Estruturais/Kalman/src/types.jl

max_issues_repo_path: Modelos Estruturais/Kalman/src/types.jl

max_forks_repo_path: Modelos Estruturais/Kalman/src/types.jl


text:
import CpuId

const REGISTER_SIZE = CpuId.simdbytes()
const REGISTER_COUNT = CpuId.cpufeature(CpuId.AVX512F) ? 32 : 16
const FP256 = (CpuId.cpufeature(CpuId.FP256)) # Is AVX2 fast?
const CACHELINE_SIZE = (CpuId.cachelinesize())
const CACHE_SIZE = CpuId.cachesize()
const NUM_CORES = CpuId.cpucores()
const FMA = (CpuId.cpufeature(CpuId.FMA3))
const AVX2 = (CpuId.cpufeature(CpuId.AVX2))
const AVX512F = (CpuId.cpufeature(CpuId.AVX512F))
const AVX512ER = (CpuId.cpufeature(CpuId.AVX512ER))
const AVX512PF = (CpuId.cpufeature(CpuId.AVX512PF))
const AVX512VL = (CpuId.cpufeature(CpuId.AVX512VL))
const AVX512BW = (CpuId.cpufeature(CpuId.AVX512BW))
const AVX512DQ = (CpuId.cpufeature(CpuId.AVX512DQ))
const AVX512CD = (CpuId.cpufeature(CpuId.AVX512CD))
const SIMD_NATIVE_INTEGERS = (CpuId.cpufeature(CpuId.AVX2))

const FMA3 = FMA


text:
# syntax: proto2
# COV_EXCL_START
using ProtoBuf
import ProtoBuf.meta

const AclEntryProto_AclEntryScopeProto = (;[
    Symbol("ACCESS") => Int32(0),
    Symbol("DEFAULT") => Int32(1),
]...)

const AclEntryProto_AclEntryTypeProto = (;[
    Symbol("USER") => Int32(0),
    Symbol("GROUP") => Int32(1),
    Symbol("MASK") => Int32(2),
    Symbol("OTHER") => Int32(3),
]...)

const AclEntryProto_FsActionProto = (;[
    Symbol("NONE") => Int32(0),
    Symbol("EXECUTE") => Int32(1),
    Symbol("WRITE") => Int32(2),
    Symbol("WRITE_EXECUTE") => Int32(3),
    Symbol("READ") => Int32(4),
    Symbol("READ_EXECUTE") => Int32(5),
    Symbol("READ_WRITE") => Int32(6),
    Symbol("PERM_ALL") => Int32(7),
]...)

mutable struct AclEntryProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function AclEntryProto(; kwargs...)
        obj = new(meta(AclEntryProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct AclEntryProto
const __meta_AclEntryProto = Ref{ProtoMeta}()
function meta(::Type{AclEntryProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_AclEntryProto)
            __meta_AclEntryProto[] = target = ProtoMeta(AclEntryProto)
            req = Symbol[:_type,:scope,:permissions]
            allflds = Pair{Symbol,Union{Type,String}}[:_type => Int32, :scope => Int32, :permissions => Int32, :name => AbstractString]
            meta(target, AclEntryProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_AclEntryProto[]
    end
end
function Base.getproperty(obj::AclEntryProto, name::Symbol)
    if name === :_type
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :scope
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :permissions
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct AclStatusProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function AclStatusProto(; kwargs...)
        obj = new(meta(AclStatusProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct AclStatusProto
const __meta_AclStatusProto = Ref{ProtoMeta}()
function meta(::Type{AclStatusProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_AclStatusProto)
            __meta_AclStatusProto[] = target = ProtoMeta(AclStatusProto)
            req = Symbol[:owner,:group,:sticky]
            allflds = Pair{Symbol,Union{Type,String}}[:owner => AbstractString, :group => AbstractString, :sticky => Bool, :entries => Base.Vector{AclEntryProto}, :permission => FsPermissionProto]
            meta(target, AclStatusProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_AclStatusProto[]
    end
end
function Base.getproperty(obj::AclStatusProto, name::Symbol)
    if name === :owner
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :group
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :sticky
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :entries
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AclEntryProto}
    elseif name === :permission
        return (obj.__protobuf_jl_internal_values[name])::FsPermissionProto
    else
        getfield(obj, name)
    end
end

mutable struct ModifyAclEntriesRequestProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ModifyAclEntriesRequestProto(; kwargs...)
        obj = new(meta(ModifyAclEntriesRequestProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ModifyAclEntriesRequestProto
const __meta_ModifyAclEntriesRequestProto = Ref{ProtoMeta}()
function meta(::Type{ModifyAclEntriesRequestProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ModifyAclEntriesRequestProto)
            __meta_ModifyAclEntriesRequestProto[] = target = ProtoMeta(ModifyAclEntriesRequestProto)
            req = Symbol[:src]
            allflds = Pair{Symbol,Union{Type,String}}[:src => AbstractString, :aclSpec => Base.Vector{AclEntryProto}]
            meta(target, ModifyAclEntriesRequestProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ModifyAclEntriesRequestProto[]
    end
end
function Base.getproperty(obj::ModifyAclEntriesRequestProto, name::Symbol)
    if name === :src
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :aclSpec
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AclEntryProto}
    else
        getfield(obj, name)
    end
end

mutable struct ModifyAclEntriesResponseProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ModifyAclEntriesResponseProto(; kwargs...)
        obj = new(meta(ModifyAclEntriesResponseProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ModifyAclEntriesResponseProto
const __meta_ModifyAclEntriesResponseProto = Ref{ProtoMeta}()
function meta(::Type{ModifyAclEntriesResponseProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ModifyAclEntriesResponseProto)
            __meta_ModifyAclEntriesResponseProto[] = target = ProtoMeta(ModifyAclEntriesResponseProto)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, ModifyAclEntriesResponseProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ModifyAclEntriesResponseProto[]
    end
end

mutable struct RemoveAclRequestProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RemoveAclRequestProto(; kwargs...)
        obj = new(meta(RemoveAclRequestProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RemoveAclRequestProto
const __meta_RemoveAclRequestProto = Ref{ProtoMeta}()
function meta(::Type{RemoveAclRequestProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RemoveAclRequestProto)
            __meta_RemoveAclRequestProto[] = target = ProtoMeta(RemoveAclRequestProto)
            req = Symbol[:src]
            allflds = Pair{Symbol,Union{Type,String}}[:src => AbstractString]
            meta(target, RemoveAclRequestProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RemoveAclRequestProto[]
    end
end
function Base.getproperty(obj::RemoveAclRequestProto, name::Symbol)
    if name === :src
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct RemoveAclResponseProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RemoveAclResponseProto(; kwargs...)
        obj = new(meta(RemoveAclResponseProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RemoveAclResponseProto
const __meta_RemoveAclResponseProto = Ref{ProtoMeta}()
function meta(::Type{RemoveAclResponseProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RemoveAclResponseProto)
            __meta_RemoveAclResponseProto[] = target = ProtoMeta(RemoveAclResponseProto)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, RemoveAclResponseProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RemoveAclResponseProto[]
    end
end

mutable struct RemoveAclEntriesRequestProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RemoveAclEntriesRequestProto(; kwargs...)
        obj = new(meta(RemoveAclEntriesRequestProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RemoveAclEntriesRequestProto
const __meta_RemoveAclEntriesRequestProto = Ref{ProtoMeta}()
function meta(::Type{RemoveAclEntriesRequestProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RemoveAclEntriesRequestProto)
            __meta_RemoveAclEntriesRequestProto[] = target = ProtoMeta(RemoveAclEntriesRequestProto)
            req = Symbol[:src]
            allflds = Pair{Symbol,Union{Type,String}}[:src => AbstractString, :aclSpec => Base.Vector{AclEntryProto}]
            meta(target, RemoveAclEntriesRequestProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RemoveAclEntriesRequestProto[]
    end
end
function Base.getproperty(obj::RemoveAclEntriesRequestProto, name::Symbol)
    if name === :src
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :aclSpec
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AclEntryProto}
    else
        getfield(obj, name)
    end
end

mutable struct RemoveAclEntriesResponseProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RemoveAclEntriesResponseProto(; kwargs...)
        obj = new(meta(RemoveAclEntriesResponseProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RemoveAclEntriesResponseProto
const __meta_RemoveAclEntriesResponseProto = Ref{ProtoMeta}()
function meta(::Type{RemoveAclEntriesResponseProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RemoveAclEntriesResponseProto)
            __meta_RemoveAclEntriesResponseProto[] = target = ProtoMeta(RemoveAclEntriesResponseProto)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, RemoveAclEntriesResponseProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RemoveAclEntriesResponseProto[]
    end
end

mutable struct RemoveDefaultAclRequestProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RemoveDefaultAclRequestProto(; kwargs...)
        obj = new(meta(RemoveDefaultAclRequestProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RemoveDefaultAclRequestProto
const __meta_RemoveDefaultAclRequestProto = Ref{ProtoMeta}()
function meta(::Type{RemoveDefaultAclRequestProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RemoveDefaultAclRequestProto)
            __meta_RemoveDefaultAclRequestProto[] = target = ProtoMeta(RemoveDefaultAclRequestProto)
            req = Symbol[:src]
            allflds = Pair{Symbol,Union{Type,String}}[:src => AbstractString]
            meta(target, RemoveDefaultAclRequestProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RemoveDefaultAclRequestProto[]
    end
end
function Base.getproperty(obj::RemoveDefaultAclRequestProto, name::Symbol)
    if name === :src
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct RemoveDefaultAclResponseProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RemoveDefaultAclResponseProto(; kwargs...)
        obj = new(meta(RemoveDefaultAclResponseProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RemoveDefaultAclResponseProto
const __meta_RemoveDefaultAclResponseProto = Ref{ProtoMeta}()
function meta(::Type{RemoveDefaultAclResponseProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RemoveDefaultAclResponseProto)
            __meta_RemoveDefaultAclResponseProto[] = target = ProtoMeta(RemoveDefaultAclResponseProto)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, RemoveDefaultAclResponseProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RemoveDefaultAclResponseProto[]
    end
end

mutable struct SetAclRequestProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SetAclRequestProto(; kwargs...)
        obj = new(meta(SetAclRequestProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SetAclRequestProto
const __meta_SetAclRequestProto = Ref{ProtoMeta}()
function meta(::Type{SetAclRequestProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SetAclRequestProto)
            __meta_SetAclRequestProto[] = target = ProtoMeta(SetAclRequestProto)
            req = Symbol[:src]
            allflds = Pair{Symbol,Union{Type,String}}[:src => AbstractString, :aclSpec => Base.Vector{AclEntryProto}]
            meta(target, SetAclRequestProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SetAclRequestProto[]
    end
end
function Base.getproperty(obj::SetAclRequestProto, name::Symbol)
    if name === :src
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :aclSpec
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AclEntryProto}
    else
        getfield(obj, name)
    end
end

mutable struct SetAclResponseProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SetAclResponseProto(; kwargs...)
        obj = new(meta(SetAclResponseProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SetAclResponseProto
const __meta_SetAclResponseProto = Ref{ProtoMeta}()
function meta(::Type{SetAclResponseProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SetAclResponseProto)
            __meta_SetAclResponseProto[] = target = ProtoMeta(SetAclResponseProto)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, SetAclResponseProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SetAclResponseProto[]
    end
end

mutable struct GetAclStatusRequestProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function GetAclStatusRequestProto(; kwargs...)
        obj = new(meta(GetAclStatusRequestProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct GetAclStatusRequestProto
const __meta_GetAclStatusRequestProto = Ref{ProtoMeta}()
function meta(::Type{GetAclStatusRequestProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_GetAclStatusRequestProto)
            __meta_GetAclStatusRequestProto[] = target = ProtoMeta(GetAclStatusRequestProto)
            req = Symbol[:src]
            allflds = Pair{Symbol,Union{Type,String}}[:src => AbstractString]
            meta(target, GetAclStatusRequestProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_GetAclStatusRequestProto[]
    end
end
function Base.getproperty(obj::GetAclStatusRequestProto, name::Symbol)
    if name === :src
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct GetAclStatusResponseProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function GetAclStatusResponseProto(; kwargs...)
        obj = new(meta(GetAclStatusResponseProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct GetAclStatusResponseProto
const __meta_GetAclStatusResponseProto = Ref{ProtoMeta}()
function meta(::Type{GetAclStatusResponseProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_GetAclStatusResponseProto)
            __meta_GetAclStatusResponseProto[] = target = ProtoMeta(GetAclStatusResponseProto)
            req = Symbol[:result]
            allflds = Pair{Symbol,Union{Type,String}}[:result => AclStatusProto]
            meta(target, GetAclStatusResponseProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_GetAclStatusResponseProto[]
    end
end
function Base.getproperty(obj::GetAclStatusResponseProto, name::Symbol)
    if name === :result
        return (obj.__protobuf_jl_internal_values[name])::AclStatusProto
    else
        getfield(obj, name)
    end
end

export AclEntryProto_AclEntryScopeProto, AclEntryProto_AclEntryTypeProto, AclEntryProto_FsActionProto, AclEntryProto, AclStatusProto, ModifyAclEntriesRequestProto, ModifyAclEntriesResponseProto, RemoveAclRequestProto, RemoveAclResponseProto, RemoveAclEntriesRequestProto, RemoveAclEntriesResponseProto, RemoveDefaultAclRequestProto, RemoveDefaultAclResponseProto, SetAclRequestProto, SetAclResponseProto, GetAclStatusRequestProto, GetAclStatusResponseProto
# COV_EXCL_STOP


text:
import Base: @propagate_inbounds

const IsingSpin = Int64
const IsingDistribution = IsingSpin[-1,1]
const IsingConf = Array{IsingSpin}

const IsingTc = 1/(1/2*log(1+sqrt(2)))

"""
Famous Ising model on a cubic lattice.

    IsingModel(; dims, L)

Create Ising model on `dims`-dimensional cubic lattice
with linear system size `L`.
"""
@with_kw_noshow struct IsingModel{C<:AbstractCubicLattice} <: Model # noshow because we override it below
    L::Int
    dims::Int
    l::C = choose_lattice(IsingModel, dims, L)
    neighs::Matrix{Int} = neighbors_lookup_table(l)
end

function choose_lattice(::Type{IsingModel}, dims::Int, L::Int)
    if dims == 1
        return Chain(L)
    elseif dims == 2
        return SquareLattice(L)
    else
        return CubicLattice(dims, L)
    end
end

"""
    IsingModel(params::Dict)
    IsingModel(params::NamedTuple)

Create an Ising model with (keyword) parameters as specified in the
dictionary/named tuple `params`.
"""
IsingModel(params::Dict{Symbol, T}) where T = IsingModel(; params...)
IsingModel(params::NamedTuple) = IsingModel(; params...)

# convenience
@inline Base.ndims(m::IsingModel) = m.dims

# cosmetics
import Base.summary
import Base.show
Base.summary(model::IsingModel) = "$(model.dims)D-Ising model"
Base.show(io::IO, model::IsingModel) = print(io, "$(model.dims)D-Ising model, L=$(model.L) ($(model.l.sites) sites)")
Base.show(io::IO, m::MIME"text/plain", model::IsingModel) = print(io, model)

# implement `Model` interface
@inline nsites(m::IsingModel) = nsites(m.l)

# implement `MC` interface
Base.rand(::Type{MC}, m::IsingModel) = rand(IsingDistribution, fill(m.L, ndims(m))...)

@propagate_inbounds function propose_local(mc::MC, m::IsingModel, i::Int, conf::IsingConf)
    delta_E = 2. * conf[i] * sum(conf[m.neighs[:,i]])
    return delta_E, nothing
end

@propagate_inbounds function accept_local!(mc::MC, m::IsingModel, i::Int, conf::IsingConf, delta_i, delta_E::Float64)
    conf[i] *= -1
    nothing
end

# optimized for 2D case
@propagate_inbounds function propose_local(mc::MC, m::IsingModel{SquareLattice}, i::Int, conf::IsingConf)
    neighs = m.neighs
    @inbounds delta_E = 2. * conf[i] * (conf[neighs[1, i]] + conf[neighs[2, i]] +
                              + conf[neighs[3, i]] + conf[neighs[4, i]])
    return delta_E, nothing
end

# optional functions/methods
"""
    global_move(mc::MC, m::IsingModel, conf::IsingConf) -> accepted::Bool

Constructs a Wolff cluster spinflip for configuration `conf`.
Returns wether a cluster spinflip has been performed (any spins have been flipped).
"""
function global_move(mc::MC, m::IsingModel, conf::IsingConf)
    N = nsites(m)
    neighs = m.neighs
    beta = mc.p.beta

    cluster = Array{Int, 1}()
    tocheck = Array{Int, 1}()

    s = rand(1:N)
    push!(tocheck, s)
    push!(cluster, s)

    while !isempty(tocheck)
        cur = pop!(tocheck)
        @inbounds for n in neighs[:,cur]

            @inbounds if conf[cur] == conf[n] && !(n in cluster) && rand() < (1 - exp(- 2.0 * beta))
                push!(tocheck, n)
                push!(cluster, n)
            end

        end
    end

    for spin in cluster
        conf[spin] *= -1
    end

    return length(cluster)>1
end

@inline function prepare_observables(mc::MC, m::IsingModel)
    obs = Dict{String,Observable}()
    obs["confs"] = Observable(typeof(mc.conf), "Configurations")

    obs["E"] = Observable(Float64, "Total energy")
    obs["E2"] = Observable(Float64, "Total energy squared")
    obs["e"] = Observable(Float64, "Energy (per site)")

    obs["M"] = Observable(Float64, "Total magnetization")
    obs["M2"] = Observable(Float64, "Total magnetization squared")
    obs["m"] = Observable(Float64, "Magnetization (per site)")

    obs["χ"] = Observable(Float64, "Susceptibility")

    obs["C"] = Observable(Float64, "Specific Heat")

    return obs
end

"""
    energy(mc::MC, m::IsingModel, conf::IsingConf)

Calculate energy of Ising configuration `conf` for Ising model `m`.
"""
function energy(mc::MC, m::IsingModel, conf::IsingConf)
    E = 0.0
    for n in 1:ndims(m)
        @inbounds @simd for i in 1:nsites(m)
            E -= conf[i]*conf[m.neighs[n,i]]
        end
    end
    return E
end

"""
    energy(mc::MC, m::IsingModel{SquareLattice}, conf::IsingConf)

Calculate energy of Ising configuration `conf` for 2D Ising model `m`.
This method is a faster variant of the general method for the
square lattice case. (It is roughly twice as fast in this case.)
"""
function energy(mc::MC, m::IsingModel{SquareLattice}, conf::IsingConf)
    neighs = m.neighs
    E = 0.0
    @inbounds @simd for i in 1:nsites(m)
        E -= conf[i]*conf[neighs[1,i]] + conf[i]*conf[neighs[2,i]]
    end
    return E
end

@inline function measure_observables!(mc::MC, m::IsingModel, obs::Dict{String,Observable}, conf::IsingConf)
    N = nsites(m)
    push!(obs["confs"], conf)

    # energie
    E = energy(mc, m, conf)
    E2 = E^2
    push!(obs["E"], E)
    push!(obs["E2"], E2)
    push!(obs["e"], E/N)

    # magnetization
    M::Float64 = abs(sum(conf))
    M2 = M^2
    push!(obs["M"], M)
    push!(obs["M2"], M2)
    push!(obs["m"], M/N)
    nothing
end

"""
    measure_observables!(mc::MC, m::IsingModel, obs::Dict{String,Observable}, conf::IsingConf)

Calculates magnetic susceptibility and specific heat and updates corresponding `Observable` objects in `obs`.

See also [`prepare_observables`](@ref) and [`measure_observables!`](@ref).
"""
@inline function finish_observables!(mc::MC, m::IsingModel, obs::Dict{String,Observable})
    N = nsites(m)
    beta = mc.p.beta

    # specific heat
    E = mean(obs["E"])
    E2 = mean(obs["E2"])
    push!(obs["C"], beta*beta*(E2/N - E*E/N))

    # susceptibility
    M = mean(obs["M"])
    M2 = mean(obs["M2"])
    push!(obs["χ"], beta*(M2/N - M*M/N))
    nothing
end

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC