text:
using Mimi
using Distributions
using CSVFiles
using DataFrames
using CSV
using Random

# standard output folder:
global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/"

# set default; interannual global mean standard deviation
global set_gvarsd = 0.11294 # default from https://github.com/jkikstra/climvar
# in case you so desire, set one default; interannual regional mean standard deviation [normally read from rvarsd_regionalvariabilitystandarddeviation.csv]
# global set_rvarsd = 0.29175 # default from mean of mean of SDs: https://github.com/jkikstra/climvar/blob/master/IlyasPastVariability.ipynb
# global v_multiplier = 1 # introduced for sensitivity analysis for variability; default = 1 -> variability according to empirical results.

# set default monte carlo parameters to the full set of probability distributions
global use_only_varMC = false
global use_no_varMC = false
global use_linear = false
global use_logburke = true
global use_logpopulation = false

# set only when doing specific variability sensitivity analysis
global globallyset_use_interannvar = false

# simple: evaluate the differences between variability and non-variability
function overview_deterministic()

    include("getpagefunction.jl")
    include("utils/mctools.jl")

    df_overview = DataFrame(ModelName = String[], SCC = Float64[], Tot_discountedimpacts = Float64[],
                                                                DiscontCost_ann = Float64[],
                                                                DiscontDays_ann = Int64[],
                                                                # BurkeDamCons_ann = Float64[],
                                                                BurkeDamGDP_ann = Float64[],
                                                                # NonMarketCons_ann = Float64[],
                                                                NonMarketGDP_ann = Float64[])

    global use_annual = false
    global use_interannvar = false

    m_PAGEICE = getpage()
    run(m_PAGEICE)
    scc_PAGEICE = compute_scc(m_PAGEICE, year=2020)

    global use_annual = true
    global use_interannvar = false

    m_nonvar = getpage()
    run(m_nonvar)
    scc_nonvar = compute_scc(m_nonvar, year=2020)

    global use_annual = true
    global use_interannvar = true
    global v_multiplier = 1

    Random.seed!(2110);
    m_var = getpage()
    run(m_var)
    scc_var = compute_scc(m_nonvar, year=2020)

    # NB rcons_per_cap_MarketRemainConsumption_ann is not a sum.
    push!(df_overview, ["PAGE-ICE", scc_PAGEICE, m_PAGEICE[:EquityWeighting, :td_totaldiscountedimpacts],
                                    m_PAGEICE[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann_sum],
                                    m_PAGEICE[:Discontinuity, :occurdis_occurrencedummy_ann_sum],
                                    # m_PAGEICE[:MarketDamagesBurke, :rcons_per_cap_MarketRemainConsumption_ann],
                                    m_PAGEICE[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann_sum],
                                    # m_PAGEICE[:NonMarketDamages, :rcons_per_cap_MarketRemainConsumption_ann],
                                    m_PAGEICE[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann_sum]])
    push!(df_overview, ["PAGE-ANN", scc_nonvar, m_nonvar[:EquityWeighting, :td_totaldiscountedimpacts_ann],
                                    m_nonvar[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann_sum],
                                    m_nonvar[:Discontinuity, :occurdis_occurrencedummy_ann_sum],
                                    # m_nonvar[:MarketDamagesBurke, :rcons_per_cap_MarketRemainConsumption_ann],
                                    m_nonvar[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann_sum],
                                    # m_nonvar[:NonMarketDamages, :rcons_per_cap_MarketRemainConsumption_ann],
                                    m_nonvar[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann_sum]])
    push!(df_overview, ["PAGE-VAR", scc_var, m_var[:EquityWeighting, :td_totaldiscountedimpacts_ann],
                                    m_var[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann_sum],
                                    m_var[:Discontinuity, :occurdis_occurrencedummy_ann_sum],
                                    # m_var[:MarketDamagesBurke, :rcons_per_cap_MarketRemainConsumption_ann],
                                    m_var[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann_sum],
                                    # m_var[:NonMarketDamages, :rcons_per_cap_MarketRemainConsumption_ann],
                                    m_var[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann_sum]])

    showall(df_overview)

end

# explore the sensitivity of the SCC for changes in the amplitude of variability for a deterministic run --- only for PAGE-VAR, ovbiously
function var_sensitivity_deterministic(scenario::String = "NDCs")
    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/VarSensitivityTests/deterministic/"
    df_varmultiplier_mybday = DataFrame(Multiplier = Float64[], SCC_2020 = Float64[], NonMarketCost = Float64[], MarketBurkeCost = Float64[], DiscontCost = Float64[])
    df_varmultiplier_vbday = DataFrame(Multiplier = Float64[], SCC_2020 = Float64[], NonMarketCost = Float64[], MarketBurkeCost = Float64[], DiscontCost = Float64[])
    df_varmultiplier_mbday = DataFrame(Multiplier = Float64[], SCC_2020 = Float64[], NonMarketCost = Float64[], MarketBurkeCost = Float64[], DiscontCost = Float64[])
    df_varmultiplier_rbday = DataFrame(Multiplier = Float64[], SCC_2020 = Float64[], NonMarketCost = Float64[], MarketBurkeCost = Float64[], DiscontCost = Float64[])
    df_varmultiplier_lbday = DataFrame(Multiplier = Float64[], SCC_2020 = Float64[], NonMarketCost = Float64[], MarketBurkeCost = Float64[], DiscontCost = Float64[])
    mybday = 21101996
    vbday = 28041993
    mbday = 20041995
    rbday = 3111994
    lbday = 2111995

    global use_annual = true
    global use_interannvar = true
    global globallyset_use_interannvar = true

    for i = collect(range(0, stop = 3, length = 31))
        global v_multiplier = i
        println("var multiplier = ", v_multiplier)

        # for random seed 1.
        Random.seed!(mybday)
        include("getpagefunction.jl")
        include("utils/mctools.jl")
        Random.seed!(mybday)
        m_varmultiplier = getpage(scenario)
        Random.seed!(mybday)
        run(m_varmultiplier)
        Random.seed!(mybday)
        global randomdraw = mybday
        scc_varmultiplier = compute_scc(m_varmultiplier, year=2020)

        push!(df_varmultiplier_mybday, [i, scc_varmultiplier, m_varmultiplier[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann_sum], m_varmultiplier[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann_sum], m_varmultiplier[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann_sum]])

        # for random seed 2.
        Random.seed!(vbday)
        include("getpagefunction.jl")
        include("utils/mctools.jl")
        Random.seed!(vbday)
        m_varmultiplier = getpage(scenario)
        Random.seed!(vbday)
        run(m_varmultiplier)
        Random.seed!(vbday)
        global randomdraw = vbday
        scc_varmultiplier = compute_scc(m_varmultiplier, year=2020)

        push!(df_varmultiplier_vbday, [i, scc_varmultiplier, m_varmultiplier[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann_sum], m_varmultiplier[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann_sum], m_varmultiplier[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann_sum]])

        # for random seed 3.
        Random.seed!(mbday)
        include("getpagefunction.jl")
        include("utils/mctools.jl")
        Random.seed!(mbday)
        m_varmultiplier = getpage(scenario)
        Random.seed!(mbday)
        run(m_varmultiplier)
        Random.seed!(mbday)
        global randomdraw = mbday
        scc_varmultiplier = compute_scc(m_varmultiplier, year=2020)

        push!(df_varmultiplier_mbday, [i, scc_varmultiplier, m_varmultiplier[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann_sum], m_varmultiplier[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann_sum], m_varmultiplier[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann_sum]])

        # for random seed 4.
        Random.seed!(rbday)
        include("getpagefunction.jl")
        include("utils/mctools.jl")
        Random.seed!(rbday)
        m_varmultiplier = getpage(scenario)
        Random.seed!(rbday)
        run(m_varmultiplier)
        Random.seed!(rbday)
        global randomdraw = rbday
        scc_varmultiplier = compute_scc(m_varmultiplier, year=2020)

        push!(df_varmultiplier_rbday, [i, scc_varmultiplier, m_varmultiplier[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann_sum], m_varmultiplier[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann_sum], m_varmultiplier[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann_sum]])

        # for random seed 5.
        Random.seed!(lbday)
        include("getpagefunction.jl")
        include("utils/mctools.jl")
        Random.seed!(lbday)
        m_varmultiplier = getpage(scenario)
        Random.seed!(lbday)
        run(m_varmultiplier)
        Random.seed!(lbday)
        global randomdraw = lbday
        scc_varmultiplier = compute_scc(m_varmultiplier, year=2020)

        push!(df_varmultiplier_lbday, [i, scc_varmultiplier, m_varmultiplier[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann_sum], m_varmultiplier[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann_sum], m_varmultiplier[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann_sum]])

    end
    # showall(df_varmultiplier)
    CSV.write(string(dir_output, "df_varmultiplier_mybday.csv"), df_varmultiplier_mybday, delim = ',')
    CSV.write(string(dir_output, "df_varmultiplier_vbday.csv"), df_varmultiplier_vbday, delim = ',')
    CSV.write(string(dir_output, "df_varmultiplier_rbday.csv"), df_varmultiplier_mbday, delim = ',')
    CSV.write(string(dir_output, "df_varmultiplier_mbday.csv"), df_varmultiplier_rbday, delim = ',')
    CSV.write(string(dir_output, "df_varmultiplier_lbday.csv"), df_varmultiplier_lbday, delim = ',')
end

# var_sensitivity_deterministic()

############# to be done!! #############
#  explore the sensitivity of the SCC for changes in the amplitude of variability for a deterministic run --- only for PAGE-VAR, ovbiously
function var_sensitivity_properMC(numMCruns::Int64 = 2000, scenario::String = "NDCs")
    ######
    ######
    ######
    ######
    ######
    ######
    ###### everything in function below this needs to be rewritten to useful MC stuff.
    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/VarSensitivityTests/montecarlo/"
    global use_no_varMC = false # needs to be full ramge, without variability fluctuating, as we push it ourselves
    df_varmultiplier = DataFrame(Multiplier = Float64[], SCC_2020 = Float64[], NonMarketCost = Float64[], MarketBurkeCost = Float64[], DiscontCost = Float64[])

    global use_annual = true
    global use_interannvar = true
    global globallyset_use_interannvar = true

    for i = collect(range(0, stop = 3, length = 31))
        global v_multiplier = i
        println("var multiplier = ", v_multiplier)

        Random.seed!(2110);
        include("getpagefunction.jl")
        include("utils/mctools.jl")

        m_varmultiplier = getpage(scenario)
        run(m_varmultiplier)
        scc_varmultiplier = compute_scc(m_varmultiplier, year=2020)

        # do monte carlo runs for every iteration here, to get an idea of the range, rather than just one draw with an arbitrary randseed.

        push!(df_varmultiplier, [i, scc_varmultiplier, m_varmultiplier[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann_sum], m_varmultiplier[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann_sum], m_varmultiplier[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann_sum]])
    end
    showall(df_varmultiplier)
    CSV.write(string(dir_output, "df_varmultiplier.csv"), df_varmultiplier, delim = ',')
end
# explore the sensitivity of random.seed by doing several runs. deterministic --- only for PAGE-VAR, to illustrate the point, but also only one affected by it in deterministic runs
function var_sensitivity_randomseed(scenario::String = "NDCs")
    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/RandomSeedSensitivity/page-var/"
    df_varmultiplier = DataFrame(RandomSeed = Int64[], SCC_2020 = Float64[], NonMarketCost = Float64[], MarketBurkeCost = Float64[], DiscontCost = Float64[])

    global use_annual = true
    global use_interannvar = true
    global globallyset_use_interannvar = false

    i = 0
    while i < 50000
        # println("random seed = ", i)

        Random.seed!(i);
        include("getpagefunction.jl")
        include("utils/mctools.jl")

        m_varmultiplier = getpage(scenario)
        run(m_varmultiplier)
        scc_varmultiplier = compute_scc(m_varmultiplier, year=2020)

        # do monte carlo runs for every iteration here, to get an idea of the range, rather than just one draw with an arbitrary randseed.

        push!(df_varmultiplier, [i, scc_varmultiplier, m_varmultiplier[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann_sum], m_varmultiplier[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann_sum], m_varmultiplier[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann_sum]])
        i = i + 1
    end
    showall(df_varmultiplier)
    CSV.write(string(dir_output, "df_varmultiplier.csv"), df_varmultiplier, delim = ',')
end

# explore the sensitivity of random.seed by doing several runs. deterministic --- only for PAGE-VAR, to illustrate the point, but also only one affected by it in deterministic runs
function var_sensitivity_randomseed_deterministic_combo(numSeedRuns::Int64 = 100, scenario::String = "NDCs")
    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/RandomSeedSensitivity_Determinstic_Combo/page-var/"
    df_varmultiplier = DataFrame(RandomSeed = Float64[], mean_SCC = Float64[], sd_SCC = Float64[])

    # GOAL: create what is in var_sensitivity_deterministic but then with means and error bars

    global use_annual = true
    global use_interannvar = true
    global globallyset_use_interannvar = true
    include("getpagefunction.jl")
    include("utils/mctools.jl")

    for i = collect(range(0, stop = 3, length = 51))
        global v_multiplier = i
        println("var multiplier = ", v_multiplier)
        ii = 1
        scc_array = zeros(numSeedRuns)

        while ii < numSeedRuns
            # println("random seed = ", ii)

            Random.seed!(ii);
            include("getpagefunction.jl")
            include("utils/mctools.jl")

            m = getpage(scenario)
            run(m)
            scc_varmultiplier = compute_scc(m, year=2020)

            scc_array[ii] = scc_varmultiplier

            ii = ii + 1
        end
        push!(df_varmultiplier, [i, mean(scc_array), std(scc_array)])
    end
    showall(df_varmultiplier)
    CSV.write(string(dir_output, "df_varmultiplier.csv"), df_varmultiplier, delim = ',')
end

# do a single run with variability
function singlerun_var(scenario::String = "NDCs")

    global use_annual = true
    global use_interannvar = true

    Random.seed!(2110);

    include("getpagefunction.jl")
    include("utils/mctools.jl")

    m_var = getpage(scenario)
    run(m_var)
    # scc_var = compute_scc(m_var, year=2020)

    # println(m_var[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann_sum])
    # explore(m_var)
end

# do Monte Carlo runs for PAGE-ICE, PAGE-ANN, and PAGE-VAR
function do_variability_MCs(numMCruns::Int64 = 2000, scenario::String = "NDCs", calc_scc::Bool=true)

    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    include("compute_scc.jl")

    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEICE"
    global use_annual = false
    global use_interannvar = false
    global v_multiplier = 0

    Random.seed!(2110);
    do_monte_carlo_runs(numMCruns, dir_output)
    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEICE/scc"
    if calc_scc
        get_scc_mcs(numMCruns, 2020, dir_output)
    end

    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    include("compute_scc.jl")

    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEANN"
    global use_annual = true
    global use_interannvar = false
    global v_multiplier = 0

    Random.seed!(2110);
    do_monte_carlo_runs(numMCruns, dir_output)
    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEANN/scc"
    if calc_scc
        get_scc_mcs(numMCruns, 2020, dir_output)
    end

    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    include("compute_scc.jl")

    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR"
    global use_annual = true
    global use_interannvar = true
    global v_multiplier = 1

    Random.seed!(2110);
    do_monte_carlo_runs(numMCruns, dir_output)
    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR/scc"
    if calc_scc
        get_scc_mcs(numMCruns, 2020, dir_output)
    end

end

# do Monte Carlo runs for multiple different scenarios
function do_var_scenario_MCs(numMCruns::Int64=100)
    global use_annual = true
    global use_interannvar = true
    randomdraw = rand(1:1000000000000)

    ## THESE TWO SCENARIOS SEEM NOT TO BE WORKING:
    # global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR/1_5C"
    # global scenario = "1.5 degC Target"
    # include("getpagefunction.jl")
    # include("utils/mctools.jl")
    # include("mcs.jl")
    # Random.seed!(2110);
    # get_scc_mcs(numMCruns, 2020, dir_output, scenario)
    #
    #
    # global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR/2_0C"
    # global scenario = "2 degC Target"
    # include("getpagefunction.jl")
    # include("utils/mctools.jl")
    # include("mcs.jl")
    # Random.seed!(2110);
    # get_scc_mcs(numMCruns, 2020, dir_output, scenario=scenario)

    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR/onlyvarMC/2_5C"
    global scenario = "2.5 degC Target"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    Random.seed!(randomdraw);
    get_scc_mcs(numMCruns, 2020, dir_output, scenario)

    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR/onlyvarMC/NDC"
    global scenario = "NDCs"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    Random.seed!(randomdraw);
    get_scc_mcs(numMCruns, 2020, dir_output, scenario)

    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR/onlyvarMC/BAU"
    global scenario = "BAU"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    Random.seed!(randomdraw);
    get_scc_mcs(numMCruns, 2020, dir_output, scenario)

    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR/onlyvarMC/RCP2_6_SSP1"
    global scenario = "RCP2.6 & SSP1"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    Random.seed!(randomdraw);
    get_scc_mcs(numMCruns, 2020, dir_output, scenario)

    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR/onlyvarMC/RCP4_5_SSP2"
    global scenario = "RCP4.5 & SSP2"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    Random.seed!(randomdraw);
    get_scc_mcs(numMCruns, 2020, dir_output, scenario)

    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR/onlyvarMC/RCP8_5_SSP5"
    global scenario = "RCP8.5 & SSP5"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    Random.seed!(randomdraw);
    get_scc_mcs(numMCruns, 2020, dir_output, scenario)
end
# do Monte Carlo runs for multiple different scenarios
function scenario_pathways()
    global use_annual = true
    global use_interannvar = true
    randomdraw = rand(1:1000000000000)

    ## THESE TWO SCENARIOS SEEM NOT TO BE WORKING:
    # global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR/1_5C"
    # global scenario = "1.5 degC Target"
    # include("getpagefunction.jl")
    # include("utils/mctools.jl")
    # include("mcs.jl")
    # Random.seed!(2110);
    # get_scc_mcs(numMCruns, 2020, dir_output, scenario)
    #
    #
    # global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcPAGEVAR/2_0C"
    # global scenario = "2 degC Target"
    # include("getpagefunction.jl")
    # include("utils/mctools.jl")
    # include("mcs.jl")
    # Random.seed!(2110);
    # get_scc_mcs(numMCruns, 2020, dir_output, scenario=scenario)

    global scenario = "2.5 degC Target"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    Random.seed!(randomdraw)
    m_25 = getpage(scenario)
    run(m_25)

    global scenario = "NDCs"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    Random.seed!(randomdraw)
    m_ndc = getpage(scenario)
    run(m_ndc)

    global scenario = "BAU"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    Random.seed!(randomdraw)
    m_bau = getpage(scenario)
    run(m_bau)

    global scenario = "RCP2.6 & SSP1"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    Random.seed!(randomdraw)
    m_26_1 = getpage(scenario)
    run(m_26_1)

    global scenario = "RCP4.5 & SSP2"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    Random.seed!(randomdraw)
    m_45_2 = getpage(scenario)
    run(m_45_2)

    global scenario = "RCP8.5 & SSP5"
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    Random.seed!(randomdraw)
    m_85_5 = getpage(scenario)
    run(m_85_5)

    ## Using DataFrames
    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/scenarioIllustration/"
    df_temp = DataFrame(year = Int64[], sc25 = Float64[], scndc = Float64[], scbau = Float64[], sc261 = Float64[], sc452 = Float64[], sc855 = Float64[])
    # df_emis = DataFrame(year = Int64[], sc25 = Float64[], scndc = Float64[], scbau = Float64[], sc261 = Float64[], sc452 = Float64[], sc855 = Float64[])
    df_gdp = DataFrame(year = Int64[], sc25 = Float64[], scndc = Float64[], scbau = Float64[], sc261 = Float64[], sc452 = Float64[], sc855 = Float64[])
    df_imp = DataFrame(year = Int64[], sc25 = Float64[], scndc = Float64[], scbau = Float64[], sc261 = Float64[], sc452 = Float64[], sc855 = Float64[])

    for year in 2015:2300
        yr = year - 2015 + 1
        # save file with annual temperatures
        push!(df_temp, [year, m_25[:ClimateTemperature, :rt_g_globaltemperature_ann][yr], m_ndc[:ClimateTemperature, :rt_g_globaltemperature_ann][yr], m_bau[:ClimateTemperature, :rt_g_globaltemperature_ann][yr], m_26_1[:ClimateTemperature, :rt_g_globaltemperature_ann][yr], m_45_2[:ClimateTemperature, :rt_g_globaltemperature_ann][yr], m_85_5[:ClimateTemperature, :rt_g_globaltemperature_ann][yr] ])
        # # save file with annual emissions
        # push!(df_emis, [year, m_25[:EquityWeighting, :cons_percap_consumption_ann[yr]], m_ndc[:EquityWeighting, :cons_percap_consumption_ann[yr]], m_bau[:EquityWeighting, :cons_percap_consumption_ann[yr]], m_26_1[:EquityWeighting, :cons_percap_consumption_ann[yr]], m_45_2[:EquityWeighting, :cons_percap_consumption_ann[yr]], m_85_5[:EquityWeighting, :cons_percap_consumption_ann[yr]] ])
        # save file with annual GDP
        push!(df_gdp, [year, m_25[:EquityWeighting, :cons_percap_consumption_ann][yr], m_ndc[:EquityWeighting, :cons_percap_consumption_ann][yr], m_bau[:EquityWeighting, :cons_percap_consumption_ann][yr], m_26_1[:EquityWeighting, :cons_percap_consumption_ann][yr], m_45_2[:EquityWeighting, :cons_percap_consumption_ann][yr], m_85_5[:EquityWeighting, :cons_percap_consumption_ann][yr] ])
        # save file with annual total impacts
        push!(df_imp, [year, m_25[:EquityWeighting, :te_totaleffect_ann_yr][yr], m_ndc[:EquityWeighting, :te_totaleffect_ann_yr][yr], m_bau[:EquityWeighting, :te_totaleffect_ann_yr][yr], m_26_1[:EquityWeighting, :te_totaleffect_ann_yr][yr], m_45_2[:EquityWeighting, :te_totaleffect_ann_yr][yr], m_85_5[:EquityWeighting, :te_totaleffect_ann_yr][yr] ])
    end
    CSV.write(string(dir_output, "df_temp.csv"), df_temp)
    # CSV.write(string(dir_output, "df_emis.csv"), df_emis)
    CSV.write(string(dir_output, "df_gdp.csv"), df_gdp)
    CSV.write(string(dir_output, "df_imp.csv"), df_imp)

end

# do deterministic runs for PAGE-ICE, PAGE-ANN, and PAGE-VAR and output various variables for analysis (using DataFrame)
function compare_deterministic_ICEANNVAR(scenario::String = "NDCs")
    # for saving the results
    df_ICE_global = DataFrame(ModelName = String[], year = String[], SCC = Float64[], g_temp = Float64[],
                                td_totdiscountimpact = Float64[])
    df_ICE_regional = DataFrame(ModelName = String[], year = String[],
                                r_temp = Float64[], gdprem_marketburke = Float64[],
                                consremain_marketburke = Float64[], gdprem_nonmarket = Float64[],
                                consremain_nonmarket = Float64[], consremain_discon = Float64[],
                                impact_equity = Float64[], adaptation_equity = Float64[]) # NOTE; no easy way to do SCC yet regionally.
    df_ANN_global = DataFrame(ModelName = String[], year = String[], SCC = Float64[], g_temp = Float64[],
                                td_totdiscountimpact = Float64[])
    df_ANN_regional = DataFrame(ModelName = String[], year = String[],
                                r_temp = Float64[], gdprem_marketburke = Float64[],
                                consremain_marketburke = Float64[], gdprem_nonmarket = Float64[],
                                consremain_nonmarket = Float64[], consremain_discon = Float64[],
                                impact_equity = Float64[], adaptation_equity = Float64[])
    df_VAR_global = DataFrame(ModelName = String[], year = String[], SCC = Float64[], g_temp = Float64[],
                                td_totdiscountimpact = Float64[])
    df_VAR_regional = DataFrame(ModelName = String[], year = String[],
                                r_temp = Float64[], gdprem_marketburke = Float64[],
                                consremain_marketburke = Float64[], gdprem_nonmarket = Float64[],
                                consremain_nonmarket = Float64[], consremain_discon = Float64[],
                                impact_equity = Float64[], adaptation_equity = Float64[])

    # timestep PAGE-ICE:
    global use_annual = false
    global use_interannvar = false
    global v_multiplier = 0

    include("getpagefunction.jl")
    include("utils/mctools.jl")

    m_ICE = getpage(scenario)
    run(m_ICE)

    # get intermediary variables
        ### variables I want: summary variables of ClimateTemperature MarketDamagesBurke, NonMarketDamages, Discontinuity, EquityWeighting
        ###   which are [annual versions]:
        ###     Climate:        rt_g_globaltemperature_ann**, rtl_realizedtemperature_ann
        ###     Burke:          rgdp_per_cap_MarketRemainGDP_ann, rcons_per_cap_MarketRemainConsumption_ann
        ###     NonMarket:      rgdp_per_cap_NonMarketRemainGDP_ann, rcons_per_cap_NonMarketRemainConsumption_ann
        ###     Discontinuity:  rcons_per_cap_DiscRemainConsumption_ann
        ###     Equity:         td_totaldiscountedimpacts**, addt_equityweightedimpact_discountedaggregated_ann, aact_equityweightedadaptation_discountedaggregated_ann
        ###   **not regional

    var1 = m_ICE[:ClimateTemperature, :rt_g_globaltemperature]
    var2 = m_ICE[:ClimateTemperature, :rtl_realizedtemperature]
    var3 = m_ICE[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP]
    var4 = m_ICE[:MarketDamagesBurke, :rcons_per_cap_MarketRemainConsumption]
    var5 = m_ICE[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP]
    var6 = m_ICE[:NonMarketDamages, :rcons_per_cap_NonMarketRemainConsumption]
    var7 = m_ICE[:Discontinuity, :rcons_per_cap_DiscRemainConsumption]
    var8 = m_ICE[:EquityWeighting, :addt_equityweightedimpact_discountedaggregated]
    var9 = m_ICE[:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated]
    var10 = m_ICE[:EquityWeighting, :td_totaldiscountedimpacts_ts]

    # create timeseries of SCC:
    t = 1
    for y in [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]
        scc_ICE = compute_scc(m_ICE, year=y)
        push!(df_ICE_global, ["PAGE-ICE", string(y), scc_ICE, var1[t], var10[t]])

        for r = 1:1
            push!(df_ICE_regional, ["PAGE-ICE", string(y),
                            var2[t,r], var3[t,r],
                            var4[t,r], var5[t,r],
                            var6[t,r], var7[t,r],
                            var8[t,r], var9[t,r]])
        end
        t = t + 1
    end

    # PAGE-ANN:
    global use_annual = true
    global use_interannvar = false
    global v_multiplier = 0

    include("getpagefunction.jl")
    include("utils/mctools.jl")

    m_ANN = getpage(scenario)
    run(m_ANN)

    # get intermediary variables
    var1 = m_ANN[:ClimateTemperature, :rt_g_globaltemperature_ann]
    var2 = m_ANN[:ClimateTemperature, :rtl_realizedtemperature_ann]
    var3 = m_ANN[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann]
    var4 = m_ANN[:MarketDamagesBurke, :rcons_per_cap_MarketRemainConsumption_ann]
    var5 = m_ANN[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann]
    var6 = m_ANN[:NonMarketDamages, :rcons_per_cap_NonMarketRemainConsumption_ann]
    var7 = m_ANN[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann]
    var8 = m_ANN[:EquityWeighting, :addt_equityweightedimpact_discountedaggregated_ann]
    var9 = m_ANN[:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated_ann]
    var10 = m_ANN[:EquityWeighting, :td_totaldiscountedimpacts_ann_yr]

    # create timeseries of SCC:
    for y in [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]
        yr = y - 2015 + 1

        scc_ANN = compute_scc(m_ANN, year=y)
        push!(df_ANN_global, ["PAGE-ANN", string(y), scc_ANN, var1[yr], var10[yr]])
    end
    for y in collect(2015:2300)
        yr = y - 2015 + 1
        for r = 1:1
            push!(df_ANN_regional, ["PAGE-ANN", string(y),
                            var2[yr,r], var3[yr,r],
                            var4[yr,r], var5[yr,r],
                            var6[yr,r], var7[yr,r],
                            var8[yr,r], var9[yr,r]])
        end
    end

    # PAGE-VAR:
    global use_annual = true
    global use_interannvar = true
    global v_multiplier = 1

    Random.seed!(2110);

    include("getpagefunction.jl")
    include("utils/mctools.jl")

    m_VAR = getpage(scenario)
    run(m_VAR)

    # get intermediary variables
    var1 = m_VAR[:ClimateTemperature, :rt_g_globaltemperature_ann]
    var2 = m_VAR[:ClimateTemperature, :rtl_realizedtemperature_ann]
    var3 = m_VAR[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP_ann]
    var4 = m_VAR[:MarketDamagesBurke, :rcons_per_cap_MarketRemainConsumption_ann]
    var5 = m_VAR[:NonMarketDamages, :rgdp_per_cap_NonMarketRemainGDP_ann]
    var6 = m_VAR[:NonMarketDamages, :rcons_per_cap_NonMarketRemainConsumption_ann]
    var7 = m_VAR[:Discontinuity, :rcons_per_cap_DiscRemainConsumption_ann]
    var8 = m_VAR[:EquityWeighting, :addt_equityweightedimpact_discountedaggregated_ann]
    var9 = m_VAR[:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated_ann]
    var10 = m_VAR[:EquityWeighting, :td_totaldiscountedimpacts_ann_yr]

    # create timeseries of SCC:
    for y in [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300]
        yr = y - 2015 + 1
        scc_VAR = compute_scc(m_ANN, year=y)
        push!(df_VAR_global, ["PAGE-VAR", string(y), scc_VAR, var1[yr], var10[yr]])
    end
    for y in collect(2015:2300)
        yr = y - 2015 + 1
        for r = 1:1
            push!(df_VAR_regional, ["PAGE-VAR", string(y),
                            var2[yr,r], var3[yr,r],
                            var4[yr,r], var5[yr,r],
                            var6[yr,r], var7[yr,r],
                            var8[yr,r], var9[yr,r]])
        end
    end

    # showall(df_ICE_global)
    # showall(df_ICE_regional)
    # showall(df_ANN_global)
    # showall(df_ANN_regional)
    # showall(df_VAR_global)
    # showall(df_VAR_regional)
    global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/ICEvANNvVAR/"
    CSV.write(string(dir_output, scenario, "/df_ICE_global.csv"), df_ICE_global, delim = ',')
    CSV.write(string(dir_output, scenario, "/df_ICE_regional.csv"), df_ICE_regional, delim = ',')
    CSV.write(string(dir_output, scenario, "/df_ANN_global.csv"), df_ANN_global, delim = ',')
    CSV.write(string(dir_output, scenario, "/df_ANN_regional.csv"), df_ANN_regional, delim = ',')
    CSV.write(string(dir_output, scenario, "/df_VAR_global.csv"), df_VAR_global, delim = ',')
    CSV.write(string(dir_output, scenario, "/df_VAR_regional.csv"), df_VAR_regional, delim = ',')

    # explore(m_ICE)
    # explore(m_ANN)
    # explore(m_VAR)
end

# do deterministic runs for PAGE-ICE, PAGE-ANN, and PAGE-VAR and output various variables for analysis (using DataFrame)
function compare_stochastic_ICEANNVAR_fullMCoptions(numMCruns::Int64 = 10000, scenario::String = "NDCs", onlyvar::Bool=false, novar::Bool=false)
    global use_only_varMC = onlyvar
    global use_no_varMC = novar

    ### PAGE-ICE ###
    if use_only_varMC
        global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/ICEANNVARmcAugust_onlyvarMC/page-ice/"
    elseif use_no_varMC
        global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/ICEANNVARmcAugust_novarMC/page-ice/"
    else
        global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/ICEANNVARmcAugust_fullMC/page-ice/"
    end
    global use_annual = false
    global use_interannvar = false
    # global v_multiplier = 0
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    Random.seed!(21101996);
    do_monte_carlo_runs(numMCruns, dir_output)
    # get_scc_mcs(numMCruns, 2020, dir_output, scenario)

    ### PAGE-ANN ###
    if use_only_varMC
        global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/ICEANNVARmcAugust_onlyvarMC/page-ann/"
    elseif use_no_varMC
        global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/ICEANNVARmcAugust_novarMC/page-ann/"
    else
        global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/ICEANNVARmcAugust_fullMC/page-ann/"
    end
    global use_annual = true
    global use_interannvar = false
    # global v_multiplier = 0
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    Random.seed!(21101996);
    do_monte_carlo_runs(numMCruns, dir_output)
    # get_scc_mcs(numMCruns, 2020, dir_output, scenario)

    ### PAGE-VAR ###
    if use_only_varMC
        global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/ICEANNVARmcAugust_onlyvarMC/page-var/"
    elseif use_no_varMC
        global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/ICEANNVARmcAugust_novarMC/page-var/"
    else
        global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/ICEANNVARmcAugust_fullMC/page-var/"
    end
    global use_annual = true
    global use_interannvar = true
    # global v_multiplier = 1
    include("getpagefunction.jl")
    include("utils/mctools.jl")
    include("mcs.jl")
    Random.seed!(21101996);
    do_monte_carlo_runs(numMCruns, dir_output)
    # get_scc_mcs(numMCruns, 2020, dir_output, scenario)

end

# function explore_discontinuity(scenario::String = "NDCs")
#     global use_annual = true
#     global use_interannvar = true
#     global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/mcDiscontuinty"
#     include("getpagefunction.jl")
#     include("utils/mctools.jl")
#     include("mcs.jl")
#     Random.seed!(2110);
#     get_scc_mcs(numMCruns, 2020, dir_output, scenario)
# end

##########################################################
# choose what to run:
##########################################################
# global use_linear = true
# global use_logburke = false
# global use_logpopulation = false
# println("Deterministic, for full linear interpolation:")
# overview_deterministic()
# global use_linear = false
# global use_logburke = true
# global use_logpopulation = false
# println("\n Deterministic, for full logarithmic interpolation:")
# overview_deterministic()
# global use_linear = false
# global use_logburke = false
# global use_logpopulation = true
# println("\n Deterministic, for full partial logarithmic (socioeconomics) interpolation:")
# overview_deterministic()

# # test pulse_size
# global use_linear = false
# global use_logburke = true
# global use_logpopulation = false
# println("\n\n Deterministic, for Burke log interpolation:")
#
# global use_annual = false
# global use_interannvar = false
# println("PAGE-ICE (timestep)")
# include("getpagefunction.jl")
# include("utils/mctools.jl")
# m_nonvar = getpage()
# run(m_nonvar)
# for pulse = [0.00001,  0.01,  1., 50.,  500.,  2000., 5000., 7500., 10000., 50000., 100000., 1000000., 10000000., 100000000.]
#     scc_nonvar = compute_scc(m_nonvar, year=2020, pulse_size = pulse)
#     println(string("SCC for pulse_size ", pulse, " = ", scc_nonvar))
# end
#
# global use_annual = true
# global use_interannvar = false
# println("PAGE-ANN (annual, no variability)")
# include("getpagefunction.jl")
# include("utils/mctools.jl")
# m_nonvar = getpage()
# run(m_nonvar)
# for pulse = [0.00001,  0.01,  1., 50.,  500.,  2000., 5000., 7500., 10000., 50000., 100000., 1000000., 10000000., 100000000.]
#     scc_nonvar = compute_scc(m_nonvar, year=2020, pulse_size = pulse)
#     println(string("SCC for pulse_size ", pulse, " = ", scc_nonvar))
# end
#
# global use_annual = true
# global use_interannvar = true
# println("PAGE-VAR (annual, with variability)")
# include("getpagefunction.jl")
# include("utils/mctools.jl")
# m_nonvar = getpage()
# run(m_nonvar)
# for pulse = [0.00001,  0.01,  1., 50.,  500.,  2000., 5000., 7500., 10000., 50000., 100000., 1000000., 10000000., 100000000.]
#     scc_nonvar = compute_scc(m_nonvar, year=2020, pulse_size = pulse)
#     println(string("SCC for pulse_size ", pulse, " = ", scc_nonvar))
# end

# var_sensitivity_deterministic()
# var_sensitivity_properMC()
# var_sensitivity_randomseed()
# var_sensitivity_randomseed_deterministic_combo(10)

### compare ICE v ANN v VAR with 3 different selected sets of MC parameters:
# ## only var:
# compare_stochastic_ICEANNVAR_fullMCoptions(5000, "NDCs", true, false)
# ## no var, all other uncertainties:
# compare_stochastic_ICEANNVAR_fullMCoptions(5000, "NDCs", false, true)
# ## full set of uncertainties:
# compare_stochastic_ICEANNVAR_fullMCoptions(5000, "NDCs", false, false)

# singlerun_var("NDCs")

# do_variability_MCs(10000, "NDCs", true)

# do_var_scenario_MCs(2000) # NOTE: still to make appropriate changes in `mcs.jl`
#
# for sc in ["2.5 degC Target", "NDCs", "BAU", "RCP2.6 & SSP1", "RCP4.5 & SSP2", "RCP8.5 & SSP5"]
#     compare_deterministic_ICEANNVAR(sc)
# end

# scenario_pathways()

## for visualising output in easy-to-use (but limited for annual) GUI (for model "m")
# explore(m)

max_stars_repo_name: jkikstra/PAGE-VAR

max_issues_repo_name: jkikstra/PAGE-VAR

max_forks_repo_name: jkikstra/PAGE-VAR


text:
# ---
# title: 1295. Find Numbers with Even Number of Digits
# id: problem1295
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/find-numbers-with-even-number-of-digits/description/>
# hidden: true
# ---
#
# Given an array `nums` of integers, return how many of them contain an **even
# number** of digits.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [12,345,2,6,7896]
#     Output: 2
#     Explanation: 12 contains 2 digits (even number of digits).
#     345 contains 3 digits (odd number of digits).
#     2 contains 1 digit (odd number of digits).
#     6 contains 1 digit (odd number of digits).
#     7896 contains 4 digits (even number of digits).
#     Therefore only 12 and 7896 contain an even number of digits.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [555,901,482,1771]
#     Output: 1
#     Explanation:
#     Only 1771 contains an even number of digits.
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 500`
#   * `1 <= nums[i] <= 10^5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
import Base.findall

"""

Replacement for `findall` for Julia versions < 1.3

$(SIGNATURES)

# Extended help

Exported if Julia VERSION.minor < 3
"""
function findall(t::Union{AbstractString,Regex}, s::AbstractString; overlap::Bool=false)
    found = UnitRange{Int}[]
    i, e = firstindex(s), lastindex(s)
    while true
        r = findnext(t, s, i)
        isnothing(r) && return found
        push!(found, r)
        j = overlap || isempty(r) ? first(r) : last(r)
        j > e && return found
        @inbounds i = nextind(s, j)
    end
end

export
    findall


text:
module TestStatisticsWithCollect

using LazyGroupBy
using Statistics
using Test
using Transducers: MapCat

≅(a, b) = isequal(a, b) || a ≈ b

dataset = [
    "mean(_)" => g -> mean.(g),
    "mean(x -> 2x, _)" => g -> mean.(x -> 2x, g),
    "var(_)" => g -> var.(g),
    "std(_)" => g -> std.(g),
]

@testset "$label" for (label, reducer) in dataset
    xs0 = [0, 7, 3, 1, 5, 9, 4, 3, 0, 5]
    @testset for idx in [1:1, 1:2, :]
        xs = xs0[idx]
        g = grouped(isodd, xs)
        d = collect.(g)
        actual = sort!(collect(reducer(g)); by = first)
        desired = sort!(collect(zip(keys(d), reducer(values(d)))); by = first)
        @test first.(actual) == first.(desired)
        @test last.(actual) ≅ last.(desired)
    end
end

dataset_with_mapping = [
    "var(x -> 2x, _)" => (g -> var.(x -> 2x, g), g -> var(2 .* g)),
    "std(x -> 2x, _)" => (g -> std.(x -> 2x, g), g -> std(2 .* g)),
    "mean(MapCat(x -> -1:x), _)" =>
        (g -> mean.(MapCat(x -> -1:x), g), g -> mean(y for x in g for y in -1:x)),
]

@testset "$label" for (label, (reducer, mapped)) in dataset_with_mapping
    xs0 = [0, 7, 3, 1, 5, 9, 4, 3, 0, 5]
    @testset for idx in [1:1, 1:2, :]
        xs = xs0[idx]
        g = grouped(isodd, xs)
        d = collect.(g)
        actual = sort!(collect(reducer(g)); by = first)
        desired = sort!(collect(zip(keys(d), mapped.(values(d)))); by = first)
        @test first.(actual) == first.(desired)
        @test last.(actual) ≅ last.(desired)
    end
end

end  # module


text:
function binarysearch(a::Vector, item, isless=isless)
    leftindex = 1
    rightindex = length(a)
    while leftindex <= rightindex
        middleindex = (leftindex + rightindex) ÷ 2
        if isless(a[middleindex], item)
            leftindex = middleindex + 1
        elseif isless(item, a[middleindex])
            rightindex = middleindex - 1
        else
            return middleindex
        end
    end
    return false
end


text:
struct HsAtom end
struct HsContext end
struct HsJIT end
struct NativeFunctionObj end
struct NativeFunctionSignature
    arg::Cstring
    res::Cstring
    _ccall::Cstring  # can't name this field `ccall` as that is a built-in in julia
end

struct CRectArray
    data::Ptr{Nothing}
    shape_ptr::Ptr{Int64}
    strides_ptr::Ptr{Int64}
end

"""
    TaggedUnion{Tuple{A, B, ...}}

Represents a tagged union over types `A`, `B` etc.
Must have a first field `tag::UInt64`
and a second field `payload` which must be some isbits type that can hold the largest
element of the union, which you can (e.g.) declare as a custom `primitive`.
This is required as Julia doens't directly support Unions in the mapping to-from C
so we store the data as arbitary bits then force it to be reinterpretted based on the tag
"""
abstract type TaggedUnion{T<:Tuple} end
function bust_union(x::TaggedUnion{U}) where U
    T = U.parameters[Int(x.tag) + 1]
    return bust_union(force_reinterpret(T, x.payload))
end
bust_union(x) = x  # not a union, leave it alone

"Forces reinterpretting `raw` as a `T`. If sizeof(raw) <: size(T) high bits will be filled with junk."
function force_reinterpret(::Type{T}, raw) where T
    isbits(raw) || throw(ArgumentError("Can only reinterpret from a isbits type"))
    isbitstype(T) || throw(ArgumentError("Can only reinterpret into a isbits type"))
    return unsafe_load(Ptr{T}(Base.pointer_from_objref(Ref(raw))))
end

"Holds data for CLit, big enough to hold the larges of its union members, such as a Float64"
primitive type CLitPayload 64 end
struct CLit <: TaggedUnion{Tuple{Int64, Int32, UInt8, Float64, Float32, UInt32, UInt64}}
    tag::UInt64
    payload::CLitPayload
end
CLit(x) = CLit(CLit_tag(x), force_reinterpret(CLitPayload, x))

CLit_tag(::Int64) = 0x0000_0000_0000_0000
CLit_tag(::Int32) = 0x0000_0000_0000_0001
CLit_tag(::UInt8) = 0x0000_0000_0000_0002
CLit_tag(::Float64) = 0x0000_0000_0000_0003
CLit_tag(::Float32) = 0x0000_0000_0000_0004
CLit_tag(::UInt32) = 0x0000_0000_0000_0005
CLit_tag(::UInt64) = 0x0000_0000_0000_0006

"Holds data for CAtom, big enough to hold the larges of its union members, which is the CRectArray"
primitive type CAtomPayload 3*64 end
struct CAtom <: TaggedUnion{Tuple{CLit, CRectArray}}
    tag::UInt64
    payload::CAtomPayload
end
CAtom(x) = CAtom(0x0000_0000_0000_0000, force_reinterpret(CAtomPayload, CLit(x)))
CAtom(x::AbstractArray) = throw(DomainError(typeof(x), "Arrays not yet supported"))
function CAtom(atm::Ptr{HsAtom})
    result = Ref{CAtom}()
    success = to_CAtom(atm, result)
    iszero(success) && throw_from_dex()
    return result[]
end

max_stars_repo_name: justinchiu/dex-lang

max_issues_repo_name: justinchiu/dex-lang

max_forks_repo_name: justinchiu/dex-lang


text:
export totalistic
function totalistic(code::Int, steps::Int, colors=3, initial_state::Vector{UInt8}=UInt8[ 1 ])
    totalistic(UInt64(code), Int64(steps), UInt8(colors), initial_state)
end
function totalistic(code::UInt64, steps::Int64, colors::UInt8=0x3, initial_state::Vector{UInt8}=UInt8[ 1 ])
    code_count = 3 * (colors - 1) + 1
    code_map = zeros(UInt8, code_count)
    for i in 1:code_count
        code_map[i] = (code ÷ Int64(colors)^(i-1)) % Int64(colors)
    end

    current_state = initial_state
    uns = 0x0
    history = AutomataColorHistory(colors)

    push!(history, current_state, uns)

    for i = 2:steps
        next_state = zeros(UInt8, length(current_state) + 2)
        for i in 1:length(next_state)
            next_state[i] = code_map[state_at(current_state, uns, i-2) + state_at(current_state, uns, i-1) + state_at(current_state, uns, i) + 1]
        end
        uns = code_map[uns * 3 + 1]
        current_state = next_state

        push!(history, next_state, uns)
    end

    return history
end


text:
module test_72_ResultDictWithMatrixOfPlots

using ModiaResult
using Unitful
ModiaResult.@usingModiaPlot

tr = [0.0, 15.0]
t0 = (tr, tr, ModiaResult.Independent)
t1 = (0.0  : 0.1 : 15.0)
t2 = (0.0  : 0.1 : 3.0)
t3 = (5.0  : 0.3 : 9.5)
t4 = (11.0 : 0.1 : 15.0)

sigA1 = 0.9*sin.(t2)u"m"
sigA2 =     cos.(t3)u"m"
sigA3 = 1.1*sin.(t4)u"m"
R2    = [[0.4 * cos(t), 0.5 * sin(t), 0.3 * cos(t)] for t in t2]u"m"
R4    = [[0.2 * cos(t), 0.3 * sin(t), 0.2 * cos(t)] for t in t4]u"m"

sigA  = ([t2,t3,t4], [sigA1,sigA2,sigA3 ], ModiaResult.Continuous)
sigB  = ([t1]      , [0.7*sin.(t1)u"m/s"], ModiaResult.Continuous)
sigC  = ([t3]      , [sin.(t3)u"N*m"]    , ModiaResult.Clocked)
r     = ([t2,t4]   , [R2,R4]             , ModiaResult.Continuous)

result = ModiaResult.ResultDict("time" => t0,
                                "sigA" => sigA,
                                "sigB" => sigB,
                                "sigC" => sigC,
                                "r"    => r,
                                defaultHeading = "Signals from test_72_ResultDictWithMatrixOfPlots",
                                hasOneTimeSignal = false)

println("\n... test_72_ResultDictWithMatrixOfPlots.jl:\n")
ModiaResult.printResultInfo(result)

plot(result, [("sigA", "sigB", "sigC"), "r[2:3]"])

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
# TODO is this really a reasonable default way of doing nulls?

import CategoricalArrays.categorical

const ReferenceType{J} = AbstractPrimitive{T} where {J<:Integer,T<:Union{J,Union{J,Missing}}}

"""
    DictEncoding{P<:ArrowVector,J} <: ArrowVector{J}
"""
struct DictEncoding{J,R<:ReferenceType,P<:ArrowVector} <: ArrowVector{J}
    refs::R
    pool::P
end
export DictEncoding

function DictEncoding{J}(refs::R, pool::P) where {J,R<:ReferenceType,P<:ArrowVector}
    DictEncoding{J,R,P}(refs, pool)
end

function DictEncoding{J}(data::Vector{UInt8}, refs_idx::Integer, len::Integer, pool::P
                        ) where {J,P<:ArrowVector}
    refs = Primitive{Int32}(data, refs_idx, len)
    DictEncoding{J}(refs, pool)
end
function DictEncoding{Union{J,Missing}}(data::Vector{UInt8}, refs_bmask_idx::Integer,
                                        refs_values_idx::Integer,
                                        len::Integer, pool::P) where {J,P<:ArrowVector}
    refs = NullablePrimitive{Int32}(data, refs_bmask_idx, refs_values_idx, len)
    DictEncoding{Union{J,Missing}}(refs, pool)
end

function DictEncoding(data::Vector{UInt8}, refs_idx::Integer, pool_idx::Integer,
                      x::CategoricalArray{J,1,U}) where {J,U}
    refs = Primitive{Int32}(data, refs_idx, getrefs(x))
    pool = Primitive{J}(data, pool_idx, getlevels(x))
    DictEncoding{J}(refs, pool)
end
function DictEncoding(data::Vector{UInt8}, refs_bmask_idx::Integer, refs_values_idx::Integer,
                      pool_idx::Integer, x::CategoricalArray{Union{J,Missing},1,U}) where {J,U}
    refs = NullablePrimitive{Int32}(data, refs_bmask_idx, refs_values_idx, getrefs(x))
    pool = Primitive{J}(data, pool_idx, getlevels(x))
    DictEncoding{J}(refs, pool)
end

function DictEncoding(data::Vector{UInt8}, i::Integer, x::CategoricalArray{J,1,U}
                     ) where {J,U}
    refs = Primitive{Int32}(data, i, getrefs(x))
    pool = createpool(data, i+refsbytes(x), x)
    DictEncoding{T}(refs, pool)
end
function DictEncoding(data::Vector{UInt8}, i::Integer, x::CategoricalArray{Union{J,Missing},1,U}
                     ) where {J,U}
    refs = NullablePrimitive{Int32}(data, i, getrefs(x))
    pool = createpool(data, i+refsbytes(x), x)
    DictEncoding{J}(refs, pool)
end

function DictEncoding(::Type{<:Array}, x::CategoricalArray)
    b = Vector{UInt8}(totalbytes(x))
    DictEncoding(b, 1, x)
end

function DictEncoding(x::CategoricalArray{J,1,U}) where {J,U}
    refs = arrowformat(getrefs(x))
    pool = arrowformat(getlevels(x))
    DictEncoding{J}(refs, pool)
end

DictEncoding(v::AbstractVector) = DictEncoding(CategoricalArray(v))

DictEncoding{J,P}(d::DictEncoding{J,P}) where {J,P} = DictEncoding{J,P}(d.refs, d.pool)
DictEncoding{J}(d::DictEncoding{J}) where J = DictEncoding{J}(d.refs, d.pool)
DictEncoding{J}(d::DictEncoding{T}) where {J,T} = DictEncoding{J}(convert(AbstractVector{J}, d[:]))
DictEncoding(d::DictEncoding{J}) where J = DictEncoding{J}(d)

copy(d::DictEncoding) = DictEncoding(d)

referencetype(d::DictEncoding{J,R,P}) where {J,K,R<:ReferenceType{K},P} = K
export referencetype

length(d::DictEncoding) = length(d.refs)

references(d::DictEncoding) = d.refs
levels(d::DictEncoding) = d.pool
export references, levels

function createpool(data::Vector{UInt8}, i::Integer, x::CategoricalArray{J,1,U}) where {J,U}
    Primitive{J}(data, i, getlevels(x))
end
function createpool(data::Vector{UInt8}, i::Integer, x::CategoricalArray{T,1,U}
                   ) where {J<:AbstractString,U,T<:Union{J,Union{J,Missing}}}
    List{J}(data, i, getlevels(x))
end

# both defined to avoid method ambiguity
isnull(d::DictEncoding, i::Integer) = isnull(d.refs, i)
isnull(d::DictEncoding, idx::AbstractVector{<:Integer}) = isnull(d.refs, idx)

unsafe_isnull(d::DictEncoding, i::Integer) = unsafe_isnull(d.refs, i)
usnafe_isnull(d::DictEncoding, idx::AbstractVector{<:Integer}) = unsafe_isnull(d.refs, idx)

# helper function for categorical
function _getrefsvec(d::DictEncoding{J,R}, idx::AbstractVector{<:Integer}
                    ) where {J,K,R<:ArrowVector{K}}
    convert(Vector{K}, d.refs[idx] .+ one(K))
end
function _getrefsvec(d::DictEncoding{J,R}, idx::AbstractVector{<:Integer}
                    ) where {J,K,R<:ArrowVector{Union{K,Missing}}}
    K[ismissing(x) ? zero(K) : x+1 for x ∈ d.refs]
end

# for now this always transfers the entire pool
function categorical(d::DictEncoding{J,R}, idx::AbstractVector{<:Integer}
                    ) where {L,J<:Union{L,Union{L,Missing}},K,R<:ReferenceType{K}}
    p = CategoricalPool{L,K}(convert(Vector, d.pool))
    refs = _getrefsvec(d, idx)
    CategoricalArray{J,1}(refs, p)
end
categorical(d::DictEncoding) = categorical(d, 1:length(d))

convert(::Type{CategoricalArray}, d::DictEncoding{J,R}) where {J,K,R<:ReferenceType{K}} = categorical(d)

@inline function getindex(d::DictEncoding{J}, i::Integer)::J where J
    @boundscheck checkbounds(d, i)
    unsafe_isnull(d, i) ? missing : d.pool[d.refs[i]+1]
end
function getindex(d::DictEncoding{J}, idx::AbstractVector{<:Integer}) where J
    @boundscheck checkbounds(d, idx)
    @inbounds o = J[getindex(d, i) for i ∈ idx]
    o
end
function getindex(d::DictEncoding{J}, idx::AbstractVector{Bool}) where J
    @boundscheck checkbounds(d, idx)
    @inbounds o = J[getindex(d, i) for i ∈ 1:length(d) if idx[i]]
    o
end
getindex(d::DictEncoding, ::Colon) = categorical(d)

nullcount(d::DictEncoding{Union{J,Missing}}) where J = nullcount(d.refs)

#====================================================================================================
    utilities specific to DictEncoding
====================================================================================================#
getrefs(x::CategoricalArray) = convert(Vector{Int32}, x.refs) .- Int32(1)
function getrefs(x::CategoricalArray{Union{J,Missing},1,U}) where {J,U}
    refs = Vector{Union{Int32,Missing}}(length(x))
    for i ∈ 1:length(x)
        x.refs[i] == 0 ? (refs[i] = missing) : (refs[i] = x.refs[i] - 1)
    end
    refs
end

getlevels(x::CategoricalArray) = x.pool.index

refsbytes(len::Integer) = padding(sizeof(Int32)*len)
refsbytes(::Type{Union{J,Missing}}, len::Integer) where J = bitmaskbytes(len) + refsbytes(len)
refsbytes(x::AbstractVector) = refsbytes(length(x))
refsbytes(::Type{Union{J,Missing}}, x::AbstractVector) where J = refsbytes(Union{J,Missing}, length(x))
refsbytes(x::AbstractVector{Union{J,Missing}}) where J = refsbytes(Union{J,Missing}, length(x))

totalbytes(x::CategoricalArray) = refsbytes(x) + totalbytes(getlevels(x))
function totalbytes(x::CategoricalArray{Union{J,Missing},1,U}) where {J,U}
    refsbytes(Union{J,Missing}, x) + totalbytes(getlevels(x))
end


text:
abstract type AbstractArchitecture end

array_type(::KA.CPU) = Array
array_type(::KA.CUDADevice) = CUDA.CuArray

# norm
xnorm(x::AbstractVector) = norm(x, 2)
xnorm(x::CUDA.CuVector) = CUBLAS.nrm2(x)

# Array initialization
xzeros(S, n) = fill!(S(undef, n), zero(eltype(S)))
xones(S, n) = fill!(S(undef, n), one(eltype(S)))


text:
using PQEcolaPoint
using Dates

# setup DTE -
T = 36.0
DTE = range(T, stop=14.0,step=-1.0) |> collect

# setup underlying range -
Sₒ = 79.74
IV = 0.4456
σ = Sₒ*(IV)*(√(T/365))
L = Sₒ - 2.96*σ
U = Sₒ + 2.96*σ
S = range(L,stop=U, length=10) |> collect

# setup IV array -
IV_array = range(0.1, stop = 1.2, step = 0.1) |> collect

# setup strike coordinates -
K = range(50.0, stop = 125.0, step = 5.0) |> collect

# what is the contract?
short_put_contract = PutContractModel()
short_put_contract.ticker = "MU"
short_put_contract.expiration_date = Date(2022, 05, 06)
short_put_contract.strike_price = 75.0
short_put_contract.premium = 0.0
short_put_contract.number_of_contracts = 1
short_put_contract.direction = 1

# compute the grid -
lattice_model = projection(short_put_contract, S, K, IV_array, DTE; number_of_levels = 80)


text:
using FFTW
using Statistics: mean
using LinearAlgebra: norm

import GPUifyLoops: @launch, @loop, @unroll, @synchronize
@hascuda using CuArrays

using Oceananigans.Operators

function ∇²_ppn!(grid::RegularCartesianGrid, f, ∇²f)
    @loop for k in (1:grid.Nz; blockIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds ∇²f[i, j, k] = ∇²_ppn(grid, f, i, j, k)
            end
        end
    end

    @synchronize
end

function mixed_fft_commutes(N)
    A = rand(N, N, N)
    Ã1 = FFTW.dct(FFTW.rfft(A, [1, 2]), 3)
    Ã2 = FFTW.rfft(FFTW.dct(A, 3), [1, 2])
    Ã1 ≈ Ã2
end

function mixed_ifft_commutes(N)
    A = rand(N, N, N)

    Ã1 = FFTW.dct(FFTW.rfft(A, [1, 2]), 3)
    Ã2 = FFTW.rfft(FFTW.dct(A, 3), [1, 2])

    A11 = FFTW.irfft(FFTW.idct(Ã1, 3), N, [1, 2])
    A12 = FFTW.idct(FFTW.irfft(Ã1, N, [1, 2]), 3)
    A21 = FFTW.irfft(FFTW.idct(Ã2, 3), N, [1, 2])
    A22 = FFTW.idct(FFTW.irfft(Ã2, N, [1, 2]), 3)
    A ≈ A11 && A ≈ A12 && A ≈ A21 && A ≈ A22
end

function fftw_planner_works(ft, Nx, Ny, Nz, planner_flag)
    g = RegularCartesianGrid(ft, (Nx, Ny, Nz), (100, 100, 100))
    RHS = CellField(Complex{ft}, CPU(), g)
    solver = PoissonSolver(g, RHS, FFTW.ESTIMATE)
    true  # Just making sure our PoissonSolver does not error/crash.
end

function poisson_ppn_planned_div_free_cpu(ft, Nx, Ny, Nz, planner_flag)
    g = RegularCartesianGrid(ft, (Nx, Ny, Nz), (100, 100, 100))

    RHS = CellField(Complex{ft}, CPU(), g)
    RHS_orig = CellField(Complex{ft}, CPU(), g)
    ϕ = CellField(Complex{ft}, CPU(), g)
    ∇²ϕ = CellField(Complex{ft}, CPU(), g)

    RHS.data .= rand(Nx, Ny, Nz)
    RHS.data .= RHS.data .- mean(RHS.data)

    RHS_orig.data .= copy(RHS.data)

    solver = PoissonSolver(g, RHS, planner_flag)

    solve_poisson_3d_ppn_planned!(solver, g, RHS, ϕ)
    ∇²_ppn!(g, ϕ, ∇²ϕ)

    ∇²ϕ.data ≈ RHS_orig.data
end

function poisson_ppn_planned_div_free_gpu(ft, Nx, Ny, Nz)
    grid = RegularCartesianGrid(ft, (Nx, Ny, Nz), (100, 100, 100))

    RHS = CellField(Complex{ft}, GPU(), grid)
    RHS_orig = CellField(Complex{ft}, GPU(), grid)
    ϕ = CellField(Complex{ft}, GPU(), grid)
    ∇²ϕ = CellField(Complex{ft}, GPU(), grid)

    solver = init_poisson_solver(GPU(), grid, RHS)

    RHS.data .= CuArray(rand(Nx, Ny, Nz))
    RHS.data .= RHS.data .- mean(RHS.data)

    RHS_orig.data .= copy(RHS.data)

    # Performing the permutation [a, b, c, d, e, f] -> [a, c, e, f, d, b] in the z-direction in preparation to calculate
    # the DCT in the Poisson solver.
    RHS.data .= cat(RHS.data[:, :, 1:2:Nz], RHS.data[:, :, Nz:-2:2]; dims=3)

    Tx, Ty = 16, 16
    Bx, By, Bz = floor(Int, Nx/Tx), floor(Int, Ny/Ty), Nz  # Blocks in grid
    solve_poisson_3d_ppn_gpu_planned!(Tx, Ty, Bx, By, Bz, solver, grid, RHS, ϕ)

    # Undoing the permutation made above to complete the IDCT.
    ϕ.data .= CuArray{eltype(ϕ.data)}(reshape(permutedims(cat(ϕ.data[:, :, 1:Int(Nz/2)], ϕ.data[:, :, end:-1:Int(Nz/2)+1]; dims=4), (1, 2, 4, 3)), Nx, Ny, Nz))
    @. ϕ.data = real(ϕ.data)

    ∇²_ppn!(grid, ϕ, ∇²ϕ)

    Array(∇²ϕ.data) ≈ Array(RHS_orig.data)
end

"""
    Test that the Poisson solver can recover an analytic solution. In this test, we are trying to see if the solver
    can recover the solution ``\\Psi(x, y, z) = cos(\\pi m_z z / L_z) sin(2\\pi m_y y / L_y) sin(2\\pi m_x x / L_x)``
    by giving it the source term or right hand side (RHS), which is ``f(x, y, z) = \\nabla^2 \\Psi(x, y, z) =
    -((\\pi m_z / L_z)^2 + (2\\pi m_y / L_y)^2 + (2\\pi m_x/L_x)^2) \\Psi(x, y, z)``.
"""
function poisson_ppn_recover_sine_cosine_solution(ft, Nx, Ny, Nz, Lx, Ly, Lz, mx, my, mz)
    grid = RegularCartesianGrid(ft, (Nx, Ny, Nz), (Lx, Ly, Lz))

    RHS = CellField(Complex{ft}, CPU(), grid)
    ϕ = CellField(Complex{ft}, CPU(), grid)

    solver = PoissonSolver(grid, RHS)

    xC, yC, zC = grid.xC, grid.yC, grid.zC
    xC = reshape(xC, (Nx, 1, 1))
    yC = reshape(yC, (1, Ny, 1))
    zC = reshape(zC, (1, 1, Nz))

    Ψ(x, y, z) = cos(π*mz*z/Lz) * sin(2π*my*y/Ly) * sin(2π*mx*x/Lx)
    f(x, y, z) = -((mz*π/Lz)^2 + (2π*my/Ly)^2 + (2π*mx/Lx)^2) * Ψ(x, y, z)

    @. RHS.data = f(xC, yC, zC)
    solve_poisson_3d_ppn_planned!(solver, grid, RHS, ϕ)

    error = norm(ϕ.data - Ψ.(xC, yC, zC)) / √(Nx*Ny*Nz)

    @info "Error (ℓ²-norm), $ft, N=($Nx, $Ny, $Nz), m=($mx, $my, $mz): $error"

    isapprox(real.(ϕ.data),  Ψ.(xC, yC, zC); rtol=5e-2)
end


text:
"""
    destination(origin::Position, distance::Real, bearing::Real, units::String="kilometers")

Take a Point or a Position and calculate the location of a destination point given a distance in
degrees, radians, miles, or kilometers; and bearing in degrees.
The destination is calculated using the [Haversine formula](http://en.wikipedia.org/wiki/Haversine_formula) to account for global curvature.

# Examples
```jldoctest
julia> point = Point([-75, 38])
Point([-75.0, 38.0])

julia> destination(point, 100, 0)
Point([-75.0, 38.8993])

julia> destination(point, 100, 45)
Point([-74.1859, 38.6331])

julia> destination(point, 50, 0, "miles")
Point([-75.0, 38.7237])
```
"""
function destination(origin::Position, distance::Real, bearing::Real, units::String="kilometers")
    lon1 = deg2rad(origin[1])
    lat1 = deg2rad(origin[2])
    bearingRad = deg2rad(bearing)
    radians = length_to_radians(distance, units)

    lat2 = asin(sin(lat1) * cos(radians) + cos(lat1) * sin(radians) * cos(bearingRad))
    lon2 = lon1 + atan(sin(bearingRad) * sin(radians) * cos(lat1), cos(radians) - sin(lat1) * sin(lat2))

    return Point([rad2deg(lon2), rad2deg(lat2)])
end

destination(origin::Point, distance::Real, bearing::Real, units::String="kilometers") = destination(origin.coordinates, distance, bearing, units)

"""
    rhumb_destination(origin::Position, distance::Real, bearing::Real, units::String="kilometers")

Take a Point or a Position and return the destination Point having travelled the given distance along a Rhumb line from the
origin Point with the (varant) given bearing.

# Examples
```jldoctest
julia> point = Point([-75, 38])
Point([-75.0, 38.0])

julia> rhumb_destination(point, 100, 0)
Point([-75.0, 38.8993])

julia> rhumb_destination(point, 100, 45)
Point([-74.1895, 38.6359])

julia> rhumb_destination(point, 50, 0, "miles")
Point([-75.0, 38.7237])
```
"""
function rhumb_destination(origin::Position, distance::Real, bearing::Real, units::String="kilometers")
    negative::Bool = distance < 0
    distanceinMeters = convert_length(abs(distance), units, "meters")
    if negative
        distanceinMeters = - abs(distanceinMeters)
    end

    dest = calculateRhumbDestination(origin, distanceinMeters, bearing)

    dest[1] += (dest[1] - origin[1] > 180) ? -360 : (origin[1] - dest[1] > 180) ? 360 : 0

    return Point(dest)
end

rhumb_destination(origin::Point, distance::Real, bearing::Real, units::String="kilometers") = rhumb_destination(origin.coordinates, distance, bearing, units)

function calculateRhumbDestination(origin::Position, distance::Real, bearing::Real)
    Δ = distance / earth_radius
    λ1 = origin[1] * pi /180
    ϕ1 = deg2rad(origin[2])
    θ = deg2rad(bearing)

    Δϕ = Δ * cos(θ)
    ϕ2 = ϕ1 + Δϕ

    if abs(ϕ2 > pi / 2)
        ϕ2 = ϕ2 > 0 ? pi - ϕ2 : -pi - ϕ2
    end

    Δψ = log(tan(ϕ2 / 2 + pi / 4) / tan(ϕ1 / 2 + pi / 4))
    q = abs(Δψ) > 10e-12 ? Δϕ / Δψ : cos(ϕ1)

    Δλ = Δ * sin(θ) / q
    λ2 = λ1 + Δλ

    return [((λ2 * 180 / pi) + 540) % 360 - 180, ϕ2 * 180 / pi]
end


text:
using ArcadeLearningEnvironment
using ArgParse
using Cambrian
using CartesianGeneticProgramming
using Dates
using IICGP
using Random

import Cambrian.mutate # function extension

s = ArgParseSettings()
@add_arg_table! s begin
    "--cfg"
    help = "configuration script"
    default = "cfg/dualcgpga_atari_pooling.yaml"
    "--game"
    help = "game rom name"
    default = "gravitar"
    "--seed"
    help = "random seed for evolution"
    arg_type = Int
    default = 0
    "--out"
    help = "output directory"
    arg_type = String
    default = dirname(@__DIR__)
    "--ind"
    help = "individual for evaluation"
    arg_type = String
    default = ""
end

args = parse_args(ARGS, s)
const rom = args["game"]
const seed = args["seed"]
const resdir = args["out"]
Random.seed!(seed)
mcfg, ecfg, ccfg, reducer, bootstrap = IICGP.dualcgp_config(args["cfg"], rom)
const max_frames = mcfg["max_frames"]
const stickiness = mcfg["stickiness"]
const grayscale = mcfg["grayscale"]
const downscale = mcfg["downscale"]
const logid = mcfg["id"]

function play_atari(
    encoder::CGPInd,
    reducer::Reducer,
    controller::CGPInd,
    seed::Int64,
    lck::ReentrantLock;
    rom=rom,
    max_frames=max_frames,
    grayscale=grayscale,
    downscale=downscale,
    stickiness=stickiness
)
    # Random.seed!(seed)
    mt = MersenneTwister(seed)
    game = Game(rom, seed, lck=lck)
    IICGP.reset!(reducer) # zero buffers
    reward = 0.0
    frames = 0
    prev_action = Int32(0)
    while ~game_over(game.ale)
        if rand(mt) > stickiness || frames == 0
            s = get_state(game, grayscale, downscale)
            output = IICGP.process(encoder, reducer, controller, s)
            action = game.actions[argmax(output)]
        else
            action = prev_action
        end
        reward += act(game.ale, action)
        frames += 1
        prev_action = action
        if frames > max_frames
            break
        end
    end
    close!(game)
    [reward]
end

# Extend Cambrian mutate function
function mutate(ind::CGPInd, ind_type::String)
    if ind_type == "encoder"
        return goldman_mutate(ecfg, ind, init_function=IPCGPInd)
    elseif ind_type == "controller"
        return goldman_mutate(ccfg, ind)
    end
end

lck = ReentrantLock()
fit(e::CGPInd, c::CGPInd, seed::Int64) = play_atari(e, reducer, c, seed, lck)

# Create an evolution framework
evo = IICGP.DualCGPGAEvo(mcfg, ecfg, ccfg, fit, logid, resdir)
#encoder_init_function=IPCGPInd, game=game

# Run evolution
init_backup(logid, resdir, args["cfg"])
run!(evo)


text:
using Fire, Distributed

"
	julia abcmk_cli.jl rates --samples 661 --gamNeg -1000,-200 --gL 1,10 --gH 400,1000 --rho 0.001 --theta 0.001 --solutions 1000000 --output rates.jld2 --dac 1,2,4,5,10,20,50,100,200,400,500,661,925,1000  --nthreads 7

Function to solve fixation and polymorphic rates analitically. The function will create N random models from prior values. Use the arguments to defined the input range for each parameter.

If rho and/or theta are set to nothing, the function will input random values given the range 0.0005:0.0005:0.01. Otherwise you can fix the values.
If gL is set to nothing, the function will not account the role of the weakly selected alleles in the estimation.

The function returns a HDF5 file containing models solved and rates. The rates will be used to compute summary statistics required at ABC.

Please check the documentation to get more info about models parameters or detailed arguments description https://jmurga.github.io/Analytical.jl/dev/cli/ to check model
"
@main function rates(;ne::Int64=1000, samples::Int64=500, gamNeg::String="-1000,-200", gL::String="5,10", gH::String="400,1000",dac::String="2,4,5,10,20,50,200,500,700",shape::Float64=0.184,rho::String="nothing",theta::String="nothing",solutions::Int64=100000,output::String="/home/jmurga/rates.jld2",scheduler::String="local",nthreads::Int64=1)

	tmpNeg    = parse.(Int,split(gamNeg,","))
	tmpStrong = parse.(Int,split(gH,","))
	dac       = parse.(Int,split(dac,","))

	if (gL == "nothing")
		tmpWeak = nothing
	else
		tmpWeak = parse.(Int,split(gL,","))
		tmpWeak = collect(tmpWeak[1]:tmpWeak[2])
	end

	if (rho == "nothing")
		rho = nothing
	else
		rho = parse(Float64,rho)
	end

	if (theta == "nothing")
		theta = nothing
	else
		theta = parse(Float64,theta)
	end

	if scheduler == "slurm"
		@eval using ClusterManagers
		@eval addprocs_slurm($nthreads)
	elseif scheduler == "htcondor"
		@eval using ClusterManagers
		@eval addprocs_htc($nthreads)
	else
		@eval addprocs($nthreads)
	end

	@eval @everywhere using Analytical, ParallelUtilities
	@eval adap = Analytical.parameters(N=$ne,n=$samples,dac=$dac,al=$shape)

	@eval convolutedSamples = Analytical.binomialDict()
	@eval Analytical.binomOp!($adap,$convolutedSamples.bn);
	@eval Analytical.rates(param = $adap,convolutedSamples=$convolutedSamples,gH=collect($tmpStrong[1]:$tmpStrong[2]),gL=$tmpWeak,gamNeg=collect($tmpNeg[1]:$tmpNeg[2]),iterations = $solutions,rho=$rho,theta=$theta,shape=$adap.al,output=$output,scheduler=$scheduler);

	for i in workers()
		rmprocs(i)
	end
end

"
	julia abcmk_cli.jl joinRates --analysisFolder --output
"
@main function joinRates(;analysisFolder::String="<folder>",samples::Int64=661,output::String="<name>")

	@eval using JLD2, CSV ,DataFrames, OrderedCollections

	searchdir(path) = filter(x->occursin("jld2",x), readdir(path))
	files = analysisFolder .* "/" .* searchdir(analysisFolder)

	j = jldopen.(files)

	models = []
	n = []
	s = []
	dsdn = []
	dac = j[1]["1000/" * string(samples)]["dac"]

	for i in eachindex(files)
		tmp = j[i]["1000/" * string(samples)]

		push!(models,tmp["models"])
		push!(n,tmp["neut"])
		push!(s,tmp["sel"])
		push!(dsdn,tmp["dsdn"])

	end

	neut = OrderedDict{Int,Array}()
	sel = OrderedDict{Int,Array}()
	o(x,k) = x[k]
	for x in keys(n[1])
		neut[x] = vcat(o.(n,x)...)
		sel[x] = vcat(o.(s,x)...)
	end
	models = vcat(models...)
	dsdn = vcat(dsdn...)

	# Writting HDF5 file
	JLD2.jldopen(output, "a+") do file
		file[string(1000)* "/" * string(samples) * "/models"] = models;
		file[string(1000)* "/" * string(samples) * "/neut"]   = neut;
		file[string(1000)* "/" * string(samples) * "/sel"]    = sel;
		file[string(1000)* "/" * string(samples) * "/dsdn"]   = dsdn;
		file[string(1000)* "/" * string(samples) * "/dac"]    = dac;
	end;

	rm.(files)
end

"
	julia abcmk_cli.jl parseData --analysisFolder analysis/ --geneList analysis/dnaVipsList.txt

Function to parse polymorphic and divergence data from Uricchio et. al (2019) and Murga-Moreno et al (2019). Please input a path to create a new analysis folder. You can filter the dataset using a file containing a list of Ensembl IDs.

The function returns files containing raw polymorphic and divergence data, parsed SFS and parsed divegence required to estimate summary statistics.

Please check the documentation to get more info https://jmurga.github.io/Analytical.jl/dev/cli/
"
@main function parseData(;analysisFolder::String="<folder>",dataset::String="tgp",geneList::String="false",bins::String="false")

	@eval using Analytical, DataFrames, CSV

	run(`mkdir -p $analysisFolder`)

	dataset = lowercase(dataset)
	data    = analysisFolder * "/" * dataset * ".txt"

	download("https://raw.githubusercontent.com/jmurga/Analytical.jl/master/data/"* dataset * ".txt",data)

	# Check if bins or genelist are defined
	@eval if $geneList != "false"
		@eval gList = CSV.read($geneList,DataFrame,header=false) |> Array
	else
		gList = nothing
	end

	@eval if $bins != "false"
		@eval binsSize = parse(Int,$bins)
	else
		binsSize = nothing
	end

	# Parsing TGP data
	if dataset == "tgp"
		@eval α,sfs, divergence = Analytical.parseSfs(sampleSize=661,data=$data,geneList=$gList,bins=$binsSize)
	elseif occursin("zi",dataset)
		@eval α,sfs, divergence = Analytical.parseSfs(sampleSize=154,data=$data,geneList=$gList,bins=$binsSize,isolines=true)
	elseif occursin("ral",dataset)
		@eval α,sfs, divergence = Analytical.parseSfs(sampleSize=160,data=$data,geneList=$gList,bins=$binsSize,isolines=true)
	end
	# Writting data to folder
	@eval sName = $analysisFolder * "/sfs.tsv"
	@eval dName = $analysisFolder * "/div.tsv"

	@eval CSV.write($sName,DataFrame($sfs,:auto),delim='\t',header=false)
	@eval CSV.write($dName,DataFrame($divergence',:auto),delim='\t',header=false)
end

"
	julia abcmk_cli.jl summaries --analysisFolder analysis/ --rates analysis/rates.jld2 --samples 661 --dac 2,4,5,10,20,50,200,661,925 --summstatSize 1000000 --nthreads 7

Estimate summary statistics from analytical rates. You must provide a path containing the parsed SFS and divergence file.

The function returns files containing bootstrapped datasets (alphas.txt) and summary statistics (summstat.txt)

Check the documentation to get more info https://jmurga.github.io/Analytical.jl/dev/cli"
@main function summaries(;analysisFolder::String="<folder>",rates::String="rates.jld2",ne::Int64=1000, samples::Int64=500,dac::String="2,4,5,10,20,50,200,500,700",summstatSize::Int64=100000,replicas::Int64=100,bootstrap::String="true")

	@eval  using Analytical, JLD2, DataFrames, CSV, ProgressMeter

	@eval h5file    = jldopen($rates)

	@eval adap      = Analytical.parameters(N=$ne,n=$samples,dac =parse.(Int,split($dac,",")))

	@eval if $bootstrap == "true"
		@eval summstat  = Analytical.summaryStatsFromRates(param=$adap,rates=$h5file,analysisFolder=$analysisFolder,summstatSize=$summstatSize,replicas=$replicas,bootstrap=true)
	else
		@eval summstat  = Analytical.summaryStatsFromRates(param=$adap,rates=$h5file,analysisFolder=$analysisFolder,summstatSize=$summstatSize,replicas=$replicas,bootstrap=false)
	end
end

"ABCreg inference.

The function returns posterior distributions from ABC inference. Each posterior file contains information about alpha_w, alpha_s, alpha, gamNeg and shape parameter. The number of posterior distributions will depend on the number of bootstrap replicas.

Check the documentation to get more info https://jmurga.github.io/Analytical.jl/dev/cli
"
@main function abcInference(;analysisFolder::String="<folder>",S::Int64=9,tol::Float64=0.01,ABCreg::String="/home/jmurga/ABCreg/src/reg")

	@eval using Analytical
	@eval Analytical.ABCreg(analysisFolder=$analysisFolder,S=$S,tol=$tol,abcreg=$ABCreg)

end

#="Plot Maximum a posterior distribution"
@main function plotMap(;analysisFolder::String="<folder>")
	try
		@eval using Analytical, RCall, GZip, DataFrames, CSV

		@eval Analytical.sourcePlotMapR(script=$analysisFolder)
		@eval Analytical.plotMap(analysisFolder=$analysisFolder)
		@eval RCall.endEmbeddedR()
	catch
		println("Please install R, ggplot2, data.table and locfit in your system before execute this function")
	end
end=#


text:
module RulesBIFMOutTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_rules

@testset "rules:BIFM:out" begin

     @testset "Belief Propagation: (m_in::MultivariateNormalDistributionsFamily, m_zprev::ProdFinal{<:MultivariateNormalDistributionsFamily}, m_znext::MultivariateNormalDistributionsFamily, meta::BIFMMeta)" begin

        meta = BIFMMeta([2.0 0; 0 1], # A
                        [3.0 0; 0 2], # B
                        [4.0 0; 0 3], # C
                        [5.0 0; 0 4], # H
                        [6.0 0; 0 5], # BHBt
                        [1.0, 2],     # ξz
                        [7.0 0; 0 6], # Λz
                        [2.0, 3],     # ξztilde
                        [8.0 0; 0 7], # Λztilde
                        [3.0, 4],     # μu
                        [9.0 0; 0 8]  # Σu
                        )

        @test_rules [ with_float_conversions = true, float32_atol = 1e-2, bigfloat_atol = 1e-8 ] BIFM(:out, Marginalisation) [
            (input = (m_in = MvNormalMeanPrecision([1,2], [2 0; 0 1]),     m_zprev = ProdFinal(MvNormalMeanPrecision([1,2], [1 0; 0 2])), m_znext = MvNormalMeanPrecision([1,2], [1 0; 0 2]), meta = meta),
                output = ProdFinal(MvNormalMeanCovariance([-280, -126], [107680.0 0.0; 0.0 3829.4999999999986]))),
            (input = (m_in = MvNormalMeanPrecision([3,4], [6 0; 0 1]),     m_zprev = ProdFinal(MvNormalMeanPrecision([1,6], [2 0; 0 2])), m_znext = MvNormalMeanPrecision([8,2], [2 0; 0 2]), meta = meta),
                output = ProdFinal(MvNormalMeanCovariance([-280, -462], [53887.99999999999 0.0; 0.0 3829.4999999999986]))),
            (input = (m_in = MvNormalMeanPrecision([5,6], [2 0; 0 5]),     m_zprev = ProdFinal(MvNormalMeanPrecision([6,2], [1 0; 0 1])), m_znext = MvNormalMeanPrecision([1,9], [1 0; 0 1]), meta = meta),
                output = ProdFinal(MvNormalMeanCovariance([-1872, -130.8], [107680.0 0.0; 0.0 7614.0])))
        ]

    end

end

end


text:
abstract type AbstractTurnover end

@columns struct Turnover{gV0,gV1,gX,gDads,gXads,FED,FVD,FVP,FV,FXD,FXP,FX} <: AbstractTurnover
    γ_V_0::gV0     | 1/h       | "Max. reserve and structure turnover rate"
    γ_V_1::gV1     | mol/m^3   | "Half-saturation constant for reserve and structure turnover rate"
    γ_X::gX        | 1/h       | "Enzyme turnover rate"
    γ_D_ads::gDads | 1/h       | "Adsorbed monomer turnover rate"
    γ_X_ads::gXads | 1/h       | "Adsorbed enzyme turnover rate"
    f_ED::FED      | _         | "Fraction of decayed reserve recycling to monomers"
    f_VD::FVD      | _         | "Fraction of decayed biomass recycling to monomer fractions"
    f_VP::FVP      | _         | "Fraction of decayed biomass recycling to polymer fractions"
    f_V::FV        | _         | "Fraction of decayed biomass recycling to monomers"
    f_XD::FXD      | _         | "Fraction of decayed enzyme recycling to monomer fractions"
    f_XP::FXP      | _         | "Fraction of decayed enzyme recycling to polymer fractions"
    f_X::FX        | _         | "Fraction of decayed enzyme recycling to polymers"
end

for fn in fieldnames(Turnover)
    @eval $fn(p::Turnover) = p.$fn
end

biomass_turnover!(J_B::Vector{Float64}, p::AbstractTurnover, B::Vector{Float64}) = begin
    J_B = γ_V_0(p).*B./(γ_V_1(p) .+ B)
end

biomass_turnover!(J_B::Vector{Float64}, p::AbstractTurnover, B::Vector{Float64}, t::String) = begin
    J_B = γ_V_0(p).*B
end

enzyme_decay!(J_X::Vector{Float64}, p::AbstractTurnover, X::Vector{Float64}) = begin
    J_X = γ_X(p).*X
end

adsorbed_monomer_decay!(J_D_ads::Vector{Float64}, p::AbstractTurnover, D_ads::Vector{Float64}) = begin
    J_D_ads = γ_D_ads(p).*D_ads
end

adsorbed_enzyme_decay!(J_X_ads::Vector{Float64}, p::AbstractTurnover, X_ads::Vector{Float64}) = begin
    J_X_ads = γ_X_ads(p).*X_ads
end

reserve_recycling!(J_ED::Vector{Float64}, p::AbstractTurnover, E::Vector{Float64}) = begin
    J_ED = f_ED(p).*sum(γ_V_0(p).*E./(γ_V_1(p) .+ E), dims=1)
end

reserve_recycling!(J_ED::Vector{Float64}, p::AbstractTurnover, E::Vector{Float64}, t::String) = begin
    J_ED = f_ED(p).*sum(γ_V_0(p).*E, dims=1)
end

biomass_recycling!(J_VD::Vector{Float64}, p::AbstractTurnover, V::Vector{Float64}) = begin
    J_VD = f_VD(p).*sum(f_V(p).*γ_V_0(p).*V./(γ_V_1(p) .+ V), dims=1)
    J_VP = f_VP(p).*sum((1 .- f_V(p)).*sum(γ_V_0(p).*V./(γ_V_1(p) .+ V), dims=1), dims=1)
    return J_VD, J_VP
end

biomass_recycling!(J_VD::Vector{Float64}, p::AbstractTurnover, V::Vector{Float64}, t::String) = begin
    J_VD = f_V(p).*f_VD(p).*sum(γ_V_0(p).*V, dims=1)
    J_VP = f_VP(p).*sum((1 .- f_V(p)).*sum(γ_V_0(p).*V, dims=1), dims=1)
    return J_VD, J_VP
end

enzyme_recycling!(J_XD::Vector{Float64}, p::AbstractTurnover, X::Vector{Float64}) = begin
    J_XD = (1 .- f_X(p)).*f_XD(p).*sum(γ_X(p).*X, dims=1)
    J_XP = f_X(p).*f_XP(p).*sum(γ_X(p).*X, dims=1)
    return J_XD, J_XP
end

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
using PrincipalMomentAnalysisApp
using Test
using PrincipalMomentAnalysis
using LinearAlgebra
using DataFrames
using CSV

import PrincipalMomentAnalysisApp: JobGraph, process_thread, process_step
using PrincipalMomentAnalysisApp.Schedulers

include("utils.jl")

@testset "PrincipalMomentAnalysisApp" begin
    include("test_content.jl")
    include("test_app.jl")
end


text:
using LinearAlgebra

import CPLEX

include("LPBenchmark.jl")

function run_cplex(fname::String, fsol::String="")

    cplex = CPLEX.Optimizer()
    bridged = MOIB.full_bridge_optimizer(cplex, Float64)

    # Read file
    load_problem!(bridged, fname)

    # MOI parameters
    MOI.set(bridged, MOI.Silent(), false)
    MOI.set(bridged, MOI.TimeLimitSec(), 10_000)
    MOI.set(bridged, MOI.RawParameter("CPXPARAM_Threads"), 1)  # Single threads
    MOI.set(bridged, MOI.RawParameter("CPXPARAM_SolutionType"), 2)  # No crossover
    MOI.set(bridged, MOI.RawParameter("CPXPARAM_LPMethod"), 4)  # barrier

    # Solve instance
    run(bridged)
    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    run_cplex(ARGS[1])
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Deprecation warnings.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

#                     Introduced in SatelliteToolbox v0.4
# ==============================================================================

@deprecate parse_gfc(filename) create_gravity_model_coefs(parse_icgem(filename))
@deprecate satellite_orbit_compute_f(a, e, i, M) M_to_f(e, M)
@deprecate satellite_orbit_compute_f(a, e, i, M, tol) M_to_f(e, M, tol)

#                     Introduced in SatelliteToolbox v0.5
# ==============================================================================

@deprecate legendre_fully_normalized(ϕ::Number, n_max::Number, ph_term::Bool)         legendre_fully_normalized(ϕ, n_max, -1, ph_term)
@deprecate legendre_schmidt_quasi_normalized(ϕ::Number, n_max::Number, ph_term::Bool) legendre_schmidt_quasi_normalized(ϕ, n_max, -1, ph_term)
@deprecate legendre_conventional(ϕ::Number, n_max::Number, ph_term::Bool)             legendre_conventional(ϕ, n_max, -1, ph_term)
@deprecate legendre(::Type{Val{:conv}}, ϕ::Number, n_max::Number, ph_term::Bool)      legendre(Val{:conv}, ϕ, n_max, -1, ph_term)
@deprecate legendre(::Type{Val{:schmidt}}, ϕ::Number, n_max::Number, ph_term::Bool)   legendre(Val{:schmidt}, ϕ, n_max, -1, ph_term)
@deprecate legendre(::Type{Val{:full}}, ϕ::Number, n_max::Number, ph_term::Bool)      legendre(Val{:full}, ϕ, n_max, -1, ph_term)
@deprecate legendre(ϕ::Number, n_max::Number, ph_term::Bool)                          legendre(ϕ, n_max, -1, ph_term)

@deprecate dlegendre_fully_normalized(ϕ::Number, n_max::Number, ph_term::Bool)         dlegendre_fully_normalized(ϕ, n_max, -1, ph_term)
@deprecate dlegendre_schmidt_quasi_normalized(ϕ::Number, n_max::Number, ph_term::Bool) dlegendre_schmidt_quasi_normalized(ϕ, n_max, -1, ph_term)
@deprecate dlegendre_conventional(ϕ::Number, n_max::Number, ph_term::Bool)             dlegendre_conventional(ϕ, n_max, -1, ph_term)
@deprecate dlegendre(::Type{Val{:conv}}, ϕ::Number, n_max::Number, ph_term::Bool)      dlegendre(Val{:conv}, ϕ, n_max, -1, ph_term)
@deprecate dlegendre(::Type{Val{:schmidt}}, ϕ::Number, n_max::Number, ph_term::Bool)   dlegendre(Val{:schmidt}, ϕ, n_max, -1, ph_term)
@deprecate dlegendre(::Type{Val{:full}}, ϕ::Number, n_max::Number, ph_term::Bool)      dlegendre(Val{:full}, ϕ, n_max, -1, ph_term)
@deprecate dlegendre(ϕ::Number, n_max::Number, ph_term::Bool)                          dlegendre(ϕ, n_max, -1, ph_term)


text:
# # struct GeneralizedEstimatingEquation <: GeneralizedEstimatingEquationsModel
# #     resp
# #     μ̂
# #     β̂
# #     id
# #     corstr
# #     ϕ
# #     vcov
# #     iter
# #     converged
# # end

# # mutable struct GeneralizedEstimatingEquation{G<:GeeResp, L<:LinPred}
# #     rr::G
# #     pp::L
# # end

# """
#     GeeResp
# The response vector and various derived vectors in a GEE model.
# """
# struct GeeResp{V<:AbstractArray{AbstractFloat,1}, D<:UnivariateDistribution, L<:Link}
#     "`y`: response vector"
#     y::V
#     d::D
#     "`eta`: the linear predictor"
#     eta::V
#     "`mu`: mean response"
#     mu::V
# end

# mutable struct GeePred{T<:BlasReal} <: LinPred
#     X::Matrix{T}
#     beta0::Vector{T}
#     iter::Real
#     converged::Bool
# end

# # step 1: generate initial estimates of β̂
# # This needs to fit the GLM and extract the GeneralizedLinearModel object
# function initialfit(::Type{GeneralizedEstimatingEquation},
#     f::FormulaTerm,
#     data,
#     d::UnivariateDistribution = Normal(),
#     l::Link = canonicallink(d))
#     return glm(f, data, d, l).model
# end

# mutable struct GeeDens
#     yᵢ
#     Xᵢ
#     μᵢ
#     rᵢ
#     β
#     iter
#     converged
# end

# function setup(id, y, X, μ, r, β)
#     uid = unique(id)
#     GeeResp(
#         [y[id .== i] for i in uid],
#         [X[id .== i, :] for i in uid],
#         [μ[id .== i] for i in uid],
#         [r[id .== i] for i in uid],
#         β,
#         0,
#         false
#     )
# end

# # step 2: using current estimates of α (function of working covariance) and ϕ, compute Vᵢ and update β̂
# # slide 34 http://biostat.jhsph.edu/~jleek/teaching/2011/754/lecture7updated.pdf
# function updateβ(β, D, V, r)
#     K = size(D, 1)
#     β + sum([D'[j] * V[j]^-1 * D[j] for j in 1:K])^-1 * sum([D'[j] * V[j]^-1 * r[j] for j in 1:K])
# end

# # step 3: using new β̂, obtain updated α and ϕ using standardized residuals
# # repeat until convergance
# function fit(::Type{GeneralizedEstimatingEquation},
#     f::FormulaTerm,
#     data,
#     d::UnivariateDistribution,
#     l::Link,
#     id,
#     corstr,
#     maxit = 20,
#     tolerance = 0.00001)

#     sch = schema(f, data)
#     f = apply_schema(f, sch)
#     y, X = modelcols(f, data)

#     # step 1: generate initial estimates of β̂
#     initial = initialfit(GeneralizedEstimatingEquation, f, data, d, l)

#     resp = GEE.setup(
#         id,
#         y,
#         X,
#         initial.rr.mu,
#         initial.rr.wrkresid,
#         initial.pp.beta0
#     )

#     n = length(unique(id))         # total number of cluster
#     m = maximum(size.(resp.Xᵢ, 1)) # number of observations within cluster

#     # step 2: compute α and ϕ
#     # step 3: using current estimates of α (function of working covariance) and ϕ, β̂
#     # Fisher scoring loop, repeat until convergence
#     local ϕ̂, D̂, V̂, R̂
#     while resp.iter <= maxit
#         resp.iter += 1
#         e = pearson.(d, resp.yᵢ, resp.μᵢ)
#         ϕ̂ = ϕ(e)
#         α̂ = updateα(corstr, n, m, e, ϕ̂)
#         R̂ = corstr(m, α̂)
#         V̂ = [V(R̂, μ) for μ in resp.μᵢ]
#         D̂ = [D(d, x, μ) for (x, μ) in zip(resp.Xᵢ, resp.μᵢ)]
#         βnew = GEE.updateβ(resp.β, D̂, V̂, resp.rᵢ)
#         if maximum(abs.((βnew .- resp.β) ./ (resp.β .+ eps(0.0)))) < tolerance
#             resp.converged = true
#             break
#         else
#             resp.β = βnew
#             resp.μᵢ = [X*resp.β for X in resp.Xᵢ]
#             resp.rᵢ = [y - μ for (y, μ) in zip(resp.yᵢ, resp.μᵢ)]
#         end
#     end

#     GeneralizedEstimatingEquation(
#         initial,
#         collect(Iterators.flatten(resp.μᵢ)),
#         resp.β,
#         id,
#         R̂,
#         ϕ̂,
#         sandwich(bread.(D̂, V̂), meat.(D̂, V̂, resp.rᵢ)),
#         resp.iter,
#         resp.converged
#     )
# end

# coef(m::GeneralizedEstimatingEquation) = m.β

# function StatsBase.coeftable(m::GeneralizedEstimatingEquation)
#     co = coef(m)
#     se = stderror(m)
#     z = co ./ se
#     pvalue = ccdf.(Chisq(1), abs2.(z))
#     CoefTable(
#         hcat(co, se, z, pvalue),
#         ["Coef.", "Std. Error", "z", "Pr(>|z|)"],
#         coefnames(m),
#         4, # pvalcol
#         3, # teststatcol
#     )
# end


text:
@testset "transform" begin

let
    pmb = PassManagerBuilder()
    dispose(pmb)
end

PassManagerBuilder() do pmb
    optlevel!(pmb, 0)
    sizelevel!(pmb, 0)

    unit_at_a_time!(pmb, false)
    unroll_loops!(pmb, false)
    simplify_libcalls!(pmb, false)
    inliner!(pmb, 0)

    Context() do ctx
    LLVM.Module("SomeModule", ctx) do mod
        FunctionPassManager(mod) do fpm
            populate!(fpm, pmb)
        end
        ModulePassManager() do mpm
            populate!(mpm, pmb)
        end
    end
    end
end

Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
ModulePassManager() do pm
    aggressive_dce!(pm)
    bit_tracking_dce!(pm)
    alignment_from_assumptions!(pm)
    cfgsimplification!(pm)
    dead_store_elimination!(pm)
    scalarizer!(pm)
    merged_load_store_motion!(pm)
    gvn!(pm)
    ind_var_simplify!(pm)
    instruction_combining!(pm)
    jump_threading!(pm)
    licm!(pm)
    loop_deletion!(pm)
    loop_idiom!(pm)
    loop_rotate!(pm)
    loop_reroll!(pm)
    loop_unroll!(pm)
    loop_unswitch!(pm)
    mem_cpy_opt!(pm)
    partially_inline_lib_calls!(pm)
    lower_switch!(pm)
    promote_memory_to_register!(pm)
    reassociate!(pm)
    sccp!(pm)
    scalar_repl_aggregates!(pm)
    scalar_repl_aggregates!(pm, 1)
    scalar_repl_aggregates_ssa!(pm)
    simplify_lib_calls!(pm)
    tail_call_elimination!(pm)
    constant_propagation!(pm)
    demote_memory_to_register!(pm)
    verifier!(pm)
    correlated_value_propagation!(pm)
    early_cse!(pm)
    lower_expect_intrinsic!(pm)
    type_based_alias_analysis!(pm)
    scoped_no_alias_aa!(pm)
    basic_alias_analysis!(pm)

    if LLVM.version() < v"8.0"
        bbvectorize!(pm)
    end
    loop_vectorize!(pm)
    slpvectorize!(pm)

    argument_promotion!(pm)
    constant_merge!(pm)
    dead_arg_elimination!(pm)
    function_attrs!(pm)
    function_inlining!(pm)
    always_inliner!(pm)
    global_dce!(pm)
    global_optimizer!(pm)
    ipconstant_propagation!(pm)
    prune_eh!(pm)
    ipsccp!(pm)
    strip_dead_prototypes!(pm)
    strip_symbols!(pm)
    internalize!(pm)
    internalize!(pm, true)
    internalize!(pm, false)
    internalize!(pm, ["SomeFunction", "SomeOtherFunction"])

    if :NVPTX in LLVM.libllvm_targets
        nvvm_reflect!(pm)
    end
end
end
end

end


text:
__precompile__()

# This file creates the Manpower module, holding all the types, definitions, and
#   methods for the manpower planning project.

module ManpowerPlanning
    using DataFrames
    using DataStructures
    using Dates
    using Distributions
    using LightGraphs
    using MetaGraphs
    using Random
    using ResumableFunctions
    using SimJulia
    using SQLite
    using StatsBase
    using XLSX

    version = v"2.7.0"

    export versionMP
    versionMP() = @info string( "Running version ", version,
        " of ManpowerPlanning module in Julia v", VERSION )

    # include( "Functions/XLSXfix.jl")

    types = [
        "attrition",
        "basenode",
        "attribute",
        "condition",
        "transition",
        "recruitment",
        "retirement",
        "compoundnode",
        "manpowersimulation",
        "multirunsimulation",
        "multirunreport",
        "subpopulation"
    ]

    rootPath = Base.source_path()
    rootPath = rootPath isa Nothing ? "" : dirname( rootPath )
    typePath = joinpath( rootPath, "Types" )
    funcPath = joinpath( rootPath, "Functions" )
    privPath = joinpath( funcPath, "private" )
    repPath = joinpath( funcPath, "reports" )
    repPrivPath = joinpath( repPath, "private" )
    simPath = joinpath( funcPath, "sim" )
    simPrivPath = joinpath( simPath, "private" )

    # The types
    map( mpType -> include( joinpath( typePath, mpType * ".jl" ) ), types )

    # Some union type aliases
    include( joinpath( typePath, "typeAliases.jl" ) )

    # The functions
    map( mpType -> include( joinpath( funcPath, mpType * ".jl" ) ), types )
    include( joinpath( funcPath, "reports.jl" ) )

    # Deprecated functions
    # include( joinpath( funcPath, "deprecated.jl" ) )
end  # module ManpowerPlanning