text:
function interpolate_parameters_equityweighting(p, v, d, t)

    # interpolation of parameters, see notes in run_timestep for why these parameters
    if is_first(t)
        for annual_year = 2015:(gettime(t))
            yr = annual_year - 2015 + 1
            for r in d.region

                v.pop_population_ann[yr, r] = p.pop_population[t, r]

                v.tct_percap_totalcosts_total_ann[yr, r] = p.tct_percap_totalcosts_total[t, r]
                v.act_adaptationcosts_total_ann[yr, r] = p.act_adaptationcosts_total[t, r]
                v.act_percap_adaptationcosts_ann[yr, r] = p.act_percap_adaptationcosts[t, r]

                v.cons_percap_consumption_ann[yr, r] = p.cons_percap_consumption[t, r]
                v.cons_percap_aftercosts_ann[yr, r] = p.cons_percap_aftercosts[t, r]

                v.grw_gdpgrowthrate_ann[yr, r] = p.grw_gdpgrowthrate[t, r]
                v.popgrw_populationgrowth_ann[yr, r] = p.popgrw_populationgrowth[t, r]

                v.yagg_periodspan_ann[yr] = 1 # p.yagg_periodspan[t] --- needs REVIEW

            end
        end
    else
        for annual_year = (gettime(t - 1) + 1):(gettime(t))
            yr = annual_year - 2015 + 1
            frac = annual_year - gettime(t - 1)
            fraction_timestep = frac / ((gettime(t)) - (gettime(t - 1)))

            for r in d.region
                if use_linear
                    v.pop_population_ann[yr, r] = p.pop_population[t, r] * fraction_timestep + p.pop_population[t - 1, r] * (1 - fraction_timestep)

                    v.tct_percap_totalcosts_total_ann[yr, r] = p.tct_percap_totalcosts_total[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)
                    v.act_adaptationcosts_total_ann[yr, r] = p.act_adaptationcosts_total[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)
                    v.act_percap_adaptationcosts_ann[yr, r] = p.act_percap_adaptationcosts[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)

                    v.cons_percap_consumption_ann[yr, r] = p.cons_percap_consumption[t, r] * fraction_timestep + p.cons_percap_consumption[t - 1, r] * (1 - fraction_timestep)
                    v.cons_percap_aftercosts_ann[yr, r] = p.cons_percap_aftercosts[t, r] * fraction_timestep + p.cons_percap_aftercosts[t - 1, r] * (1 - fraction_timestep)

                    v.grw_gdpgrowthrate_ann[yr, r] = p.grw_gdpgrowthrate[t, r] * fraction_timestep + p.grw_gdpgrowthrate[t - 1, r] * (1 - fraction_timestep)
                    v.popgrw_populationgrowth_ann[yr, r] = p.popgrw_populationgrowth[t, r] * fraction_timestep + p.popgrw_populationgrowth[t - 1, r] * (1 - fraction_timestep)
                elseif use_logburke
                    ## fully linear (partially because everything except for pop_population causes imaginary numbers (due to negative numbers)).
                    # linear
                    v.tct_percap_totalcosts_total_ann[yr, r] = p.tct_percap_totalcosts_total[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)
                    v.act_adaptationcosts_total_ann[yr, r] = p.act_adaptationcosts_total[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)
                    v.act_percap_adaptationcosts_ann[yr, r] = p.act_percap_adaptationcosts[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)

                    v.cons_percap_consumption_ann[yr, r] = p.cons_percap_consumption[t, r] * fraction_timestep + p.cons_percap_consumption[t - 1, r] * (1 - fraction_timestep)
                    v.cons_percap_aftercosts_ann[yr, r] = p.cons_percap_aftercosts[t, r] * fraction_timestep + p.cons_percap_aftercosts[t - 1, r] * (1 - fraction_timestep)

                    v.grw_gdpgrowthrate_ann[yr, r] = p.grw_gdpgrowthrate[t, r] * fraction_timestep + p.grw_gdpgrowthrate[t - 1, r] * (1 - fraction_timestep)
                    v.popgrw_populationgrowth_ann[yr, r] = p.popgrw_populationgrowth[t, r] * fraction_timestep + p.popgrw_populationgrowth[t - 1, r] * (1 - fraction_timestep)

                    v.pop_population_ann[yr, r] = p.pop_population[t, r] * fraction_timestep + p.pop_population[t - 1, r] * (1 - fraction_timestep)

                elseif use_logpopulation
                    # linear
                    v.tct_percap_totalcosts_total_ann[yr, r] = p.tct_percap_totalcosts_total[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)
                    v.act_adaptationcosts_total_ann[yr, r] = p.act_adaptationcosts_total[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)
                    v.act_percap_adaptationcosts_ann[yr, r] = p.act_percap_adaptationcosts[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)

                    v.cons_percap_consumption_ann[yr, r] = p.cons_percap_consumption[t, r] * fraction_timestep + p.cons_percap_consumption[t - 1, r] * (1 - fraction_timestep)
                    v.cons_percap_aftercosts_ann[yr, r] = p.cons_percap_aftercosts[t, r] * fraction_timestep + p.cons_percap_aftercosts[t - 1, r] * (1 - fraction_timestep)

                    v.grw_gdpgrowthrate_ann[yr, r] = p.grw_gdpgrowthrate[t, r] * fraction_timestep + p.grw_gdpgrowthrate[t - 1, r] * (1 - fraction_timestep)
                    v.popgrw_populationgrowth_ann[yr, r] = p.popgrw_populationgrowth[t, r] * fraction_timestep + p.popgrw_populationgrowth[t - 1, r] * (1 - fraction_timestep)

                    # log
                    v.pop_population_ann[yr, r] = p.pop_population[t, r]^fraction_timestep * p.pop_population[t - 1, r]^(1 - fraction_timestep)
                elseif use_logwherepossible
                    # linear
                    v.tct_percap_totalcosts_total_ann[yr, r] = p.tct_percap_totalcosts_total[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)
                    v.act_adaptationcosts_total_ann[yr, r] = p.act_adaptationcosts_total[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)
                    v.act_percap_adaptationcosts_ann[yr, r] = p.act_percap_adaptationcosts[t, r] * fraction_timestep + p.tct_percap_totalcosts_total[t - 1, r] * (1 - fraction_timestep)

                    v.cons_percap_consumption_ann[yr, r] = p.cons_percap_consumption[t, r] * fraction_timestep + p.cons_percap_consumption[t - 1, r] * (1 - fraction_timestep)
                    v.cons_percap_aftercosts_ann[yr, r] = p.cons_percap_aftercosts[t, r] * fraction_timestep + p.cons_percap_aftercosts[t - 1, r] * (1 - fraction_timestep)

                    v.grw_gdpgrowthrate_ann[yr, r] = p.grw_gdpgrowthrate[t, r] * fraction_timestep + p.grw_gdpgrowthrate[t - 1, r] * (1 - fraction_timestep)
                    v.popgrw_populationgrowth_ann[yr, r] = p.popgrw_populationgrowth[t, r] * fraction_timestep + p.popgrw_populationgrowth[t - 1, r] * (1 - fraction_timestep)

                    # log
                    v.pop_population_ann[yr, r] = p.pop_population[t, r]^fraction_timestep * p.pop_population[t - 1, r]^(1 - fraction_timestep)
                else
                    error("NO INTERPOLATION METHOD SELECTED! Specify linear or logarithmic interpolation.")
                end

                v.yagg_periodspan_ann[yr] = 1
            end
        end
    end
end

function calc_equityweighting(p, v, d, t, annual_year, r)
    # setting the year for entry in lists.
    yr = annual_year - 2015 + 1 # + 1 because of 1-based indexing in Julia

    v.df_utilitydiscountfactor_ann[yr] = (1 + p.ptp_timepreference / 100)^(-(p.y_year_ann[yr] - p.y_year_0))

    for rr in d.region
        ## Gas Costs Accounting
        # Weighted costs (Page 23 of Hope 2009)
        v.wtct_percap_weightedcosts_ann[yr, rr] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (v.cons_percap_consumption_ann[yr, rr]^(1 - p.emuc_utilityconvexity) - (v.cons_percap_consumption_ann[yr, rr] - v.tct_percap_totalcosts_total_ann[yr, rr])^(1 - p.emuc_utilityconvexity))

        # Add these into consumption
        v.eact_percap_weightedadaptationcosts_ann[yr, rr] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (v.cons_percap_consumption_ann[yr, rr]^(1 - p.emuc_utilityconvexity) - (v.cons_percap_consumption_ann[yr, rr] - v.act_percap_adaptationcosts_ann[yr, rr])^(1 - p.emuc_utilityconvexity))

        # Do partial weighting
        if p.equity_proportion == 0
            v.pct_percap_partiallyweighted_ann[yr, rr] = v.tct_percap_totalcosts_total_ann[yr, rr]
            v.wact_percap_partiallyweighted_ann[yr, rr] = v.act_percap_adaptationcosts_ann[yr, rr]
        else
            v.pct_percap_partiallyweighted_ann[yr, rr] = (1 - p.equity_proportion) * v.tct_percap_totalcosts_total_ann[yr, rr] + p.equity_proportion * v.wtct_percap_weightedcosts_ann[yr, rr]
            v.wact_percap_partiallyweighted_ann[yr, rr] = (1 - p.equity_proportion) * v.act_percap_adaptationcosts_ann[yr, rr] + p.equity_proportion * v.eact_percap_weightedadaptationcosts_ann[yr, rr]
        end

        v.pct_partiallyweighted_ann[yr, rr] = v.pct_percap_partiallyweighted_ann[yr, rr] * v.pop_population_ann[yr, rr]
        v.wact_partiallyweighted_ann[yr, rr] = v.wact_percap_partiallyweighted_ann[yr, rr] * v.pop_population_ann[yr, rr]

        # Discount rate calculations
        v.dr_discountrate_ann[yr, rr] = p.ptp_timepreference + p.emuc_utilityconvexity * (v.grw_gdpgrowthrate_ann[yr, rr] - v.popgrw_populationgrowth_ann[yr, rr])
        v.yp_yearsperiod_ann[yr] = 1 # every timestep is 1 year long.
        if yr == 1
            v.dfc_consumptiondiscountrate_ann[yr, rr] = (1 + v.dr_discountrate_ann[yr, rr] / 100)^(-v.yp_yearsperiod[TimestepIndex(1)])
        else
            v.dfc_consumptiondiscountrate_ann[yr, rr] = v.dfc_consumptiondiscountrate_ann[yr - 1, rr] * (1 + v.dr_discountrate_ann[yr, rr] / 100)^(-v.yp_yearsperiod_ann[yr])
        end

        # Discounted costs
        if p.equity_proportion == 0
            v.pcdt_partiallyweighted_discounted_ann[yr, rr] = v.pct_partiallyweighted_ann[yr, rr] * v.dfc_consumptiondiscountrate_ann[yr, rr]
            v.wacdt_partiallyweighted_discounted_ann[yr, rr] = v.act_adaptationcosts_total_ann[yr, rr] * v.dfc_consumptiondiscountrate_ann[yr, rr]
            v.wit_partiallyweighted_ann[yr, rr] = (v.cons_percap_aftercosts_ann[yr, rr] - p.rcons_percap_dis_ann[yr, rr]) * v.pop_population_ann[yr, rr] # equivalent to emuc = 0
            v.widt_partiallyweighted_discounted_ann[yr, rr] = v.wit_partiallyweighted_ann[yr, rr] * v.dfc_consumptiondiscountrate_ann[yr] # apply Ramsey discounting
        else
            v.pcdt_partiallyweighted_discounted_ann[yr, rr] = v.pct_partiallyweighted_ann[yr, rr] * v.df_utilitydiscountfactor_ann[yr]
            v.wacdt_partiallyweighted_discounted_ann[yr, rr] = v.wact_partiallyweighted_ann[yr, rr] * v.df_utilitydiscountfactor_ann[yr]
            ## Equity weighted impacts (end of page 28, Hope 2009)
            v.wit_partiallyweighted_ann[yr, rr] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (v.cons_percap_aftercosts_ann[yr, rr]^(1 - p.emuc_utilityconvexity) - p.rcons_percap_dis_ann[yr, rr]^(1 - p.emuc_utilityconvexity)) * v.pop_population_ann[yr, rr]
            v.widt_partiallyweighted_discounted_ann[yr, rr] = v.wit_partiallyweighted_ann[yr, rr] * v.df_utilitydiscountfactor_ann[yr]
        end
        v.pcdat_partiallyweighted_discountedaggregated_ann[yr, rr] = v.pcdt_partiallyweighted_discounted_ann[yr, rr] * v.yagg_periodspan_ann[yr]

        v.addt_equityweightedimpact_discountedaggregated_ann[yr, rr] = v.widt_partiallyweighted_discounted_ann[yr, rr] * v.yagg_periodspan_ann[yr]
        v.aact_equityweightedadaptation_discountedaggregated_ann[yr, rr] = v.wacdt_partiallyweighted_discounted_ann[yr, rr] * v.yagg_periodspan_ann[yr]
    end

    # some new annual variables that sum over the regions
    if isequal(r, 8)
        v.pct_g_partiallyweighted_global_ann[yr] = sum(v.pct_partiallyweighted_ann[yr, :])
        v.pcdt_g_partiallyweighted_discountedglobal_ann[yr] = sum(v.pcdt_partiallyweighted_discounted_ann[yr, :])

        v.tpc_totalaggregatedcosts_ann = v.tpc_totalaggregatedcosts_ann + sum(v.pcdat_partiallyweighted_discountedaggregated_ann[yr, :])

        v.addt_gt_equityweightedimpact_discountedglobal_ann = v.addt_gt_equityweightedimpact_discountedglobal_ann + sum(v.addt_equityweightedimpact_discountedaggregated_ann[yr, :])

        v.tac_totaladaptationcosts_ann = v.tac_totaladaptationcosts_ann + sum(v.aact_equityweightedadaptation_discountedaggregated_ann[yr, :])

        # without civilisation value
        # v.td_totaldiscountedimpacts_ann = v.addt_gt_equityweightedimpact_discountedglobal_ann
        # with civilisation value
        v.td_totaldiscountedimpacts_ann = min(v.addt_gt_equityweightedimpact_discountedglobal_ann, p.civvalue_civilizationvalue)

        # Total effect of climate change
        # without civilisation value
        # v.te_totaleffect_ann = v.td_totaldiscountedimpacts_ann + v.tpc_totalaggregatedcosts_ann + v.tac_totaladaptationcosts_ann
        # with civilisation value
        v.te_totaleffect_ann = min(v.td_totaldiscountedimpacts_ann + v.tpc_totalaggregatedcosts_ann + v.tac_totaladaptationcosts_ann, p.civvalue_civilizationvalue)

        # for analysis / debugging / code analysis
        v.td_totaldiscountedimpacts_ann_yr[yr] = v.td_totaldiscountedimpacts_ann
        v.te_totaleffect_ann_yr[yr] = v.te_totaleffect_ann
    end
end

@defcomp EquityWeighting begin

    region = Index()
    year = Index()

    # Basic information
    y_year = Parameter(index=[time], unit="year")
    y_year_0 = Parameter(unit="year")
    y_year_ann = Parameter(index=[year], unit="year")

    # Impacts across all gases
    pop_population = Parameter(index=[time, region], unit="million person")
    pop_population_ann = Variable(index=[year, region], unit="million person")

    # Total and Per-Capita Abatement and Adaptation Costs
    tct_percap_totalcosts_total = Parameter(index=[time, region], unit="\$/person")
    tct_percap_totalcosts_total_ann = Variable(index=[year, region], unit="\$/person")
    act_adaptationcosts_total = Parameter(index=[time, region], unit="\$million")
    act_adaptationcosts_total_ann = Variable(index=[year, region], unit="\$million")
    act_percap_adaptationcosts = Parameter(index=[time, region], unit="\$/person")
    act_percap_adaptationcosts_ann = Variable(index=[year, region], unit="\$/person")

    # Consumption
    cons_percap_consumption_0 = Parameter(index=[region], unit="\$/person")
    cons_percap_consumption = Parameter(index=[time, region], unit="\$/person") # Called "CONS_PER_CAP"
    cons_percap_consumption_ann = Variable(index=[year, region], unit="\$/person")
    cons_percap_aftercosts = Parameter(index=[time, region], unit="\$/person")
    cons_percap_aftercosts_ann = Variable(index=[year, region], unit="\$/person")

    # Calculation of weighted costs
    emuc_utilityconvexity = Parameter(unit="none", default=1.1666666666666667)

    wtct_percap_weightedcosts = Variable(index=[time, region], unit="\$/person")
    wtct_percap_weightedcosts_ann = Variable(index=[year, region], unit="\$/person")
    eact_percap_weightedadaptationcosts = Variable(index=[time, region], unit="\$/person")
    eact_percap_weightedadaptationcosts_ann = Variable(index=[year, region], unit="\$/person")
    wact_percap_partiallyweighted = Variable(index=[time, region], unit="\$/person")
    wact_percap_partiallyweighted_ann = Variable(index=[year, region], unit="\$/person")
    wact_partiallyweighted = Variable(index=[time, region], unit="\$million")
    wact_partiallyweighted_ann = Variable(index=[year, region], unit="\$million")

    # Amount of equity weighting variable (0, (0, 1), or 1)
    equity_proportion = Parameter(unit="fraction", default=1.0)

    pct_percap_partiallyweighted = Variable(index=[time, region], unit="\$/person")
    pct_percap_partiallyweighted_ann = Variable(index=[year, region], unit="\$/person")
    pct_partiallyweighted = Variable(index=[time, region], unit="\$million")
    pct_partiallyweighted_ann = Variable(index=[year, region], unit="\$million")
    pct_g_partiallyweighted_global = Variable(index=[time], unit="\$million")
    pct_g_partiallyweighted_global_ann = Variable(index=[year], unit="\$million")

    # Discount rates
    ptp_timepreference = Parameter(unit="%/year", default=1.0333333333333334) # <0.1,1, 2>
    grw_gdpgrowthrate = Parameter(index=[time, region], unit="%/year")
    grw_gdpgrowthrate_ann = Variable(index=[year, region], unit="%/year")
    popgrw_populationgrowth = Parameter(index=[time, region], unit="%/year")
    popgrw_populationgrowth_ann = Variable(index=[year, region], unit="%/year")

    dr_discountrate = Variable(index=[time, region], unit="%/year")
    dr_discountrate_ann = Variable(index=[year, region], unit="%/year")
    yp_yearsperiod = Variable(index=[time], unit="year") # defined differently from yagg
    yp_yearsperiod_ann = Variable(index=[year], unit="year")
    dfc_consumptiondiscountrate = Variable(index=[time, region], unit="1/year")
    dfc_consumptiondiscountrate_ann = Variable(index=[year, region], unit="1/year")

    df_utilitydiscountfactor = Variable(index=[time], unit="fraction")
    df_utilitydiscountfactor_ann = Variable(index=[year], unit="fraction")

    # Discounted costs
    pcdt_partiallyweighted_discounted = Variable(index=[time, region], unit="\$million")
    pcdt_partiallyweighted_discounted_ann = Variable(index=[year, region], unit="\$million")
    pcdt_g_partiallyweighted_discountedglobal = Variable(index=[time], unit="\$million")
    pcdt_g_partiallyweighted_discountedglobal_ann = Variable(index=[year], unit="\$million")

    pcdat_partiallyweighted_discountedaggregated = Variable(index=[time, region], unit="\$million")
    pcdat_partiallyweighted_discountedaggregated_ann = Variable(index=[year, region], unit="\$million")
    tpc_totalaggregatedcosts = Variable(unit="\$million")
    tpc_totalaggregatedcosts_ann = Variable(unit="\$million")

    wacdt_partiallyweighted_discounted = Variable(index=[time, region], unit="\$million")
    wacdt_partiallyweighted_discounted_ann = Variable(index=[year, region], unit="\$million")

    # Equity weighted impact totals
    rcons_percap_dis = Parameter(index=[time, region], unit="\$/person")
    rcons_percap_dis_ann = Parameter(index=[year, region], unit="\$/person")

    wit_partiallyweighted = Variable(index=[time, region], unit="\$million")
    wit_partiallyweighted_ann = Variable(index=[year, region], unit="\$million")
    widt_partiallyweighted_discounted = Variable(index=[time, region], unit="\$million")
    widt_partiallyweighted_discounted_ann = Variable(index=[year, region], unit="\$million")

    yagg_periodspan = Parameter(index=[time], unit="year")
    yagg_periodspan_ann = Variable(index=[year], unit="year")

    addt_equityweightedimpact_discountedaggregated = Variable(index=[time, region], unit="\$million")
    addt_equityweightedimpact_discountedaggregated_ann = Variable(index=[year, region], unit="\$million")
    addt_gt_equityweightedimpact_discountedglobal = Variable(unit="\$million")
    addt_gt_equityweightedimpact_discountedglobal_ann = Variable(unit="\$million")

    civvalue_civilizationvalue = Parameter(unit="\$million", default=6.1333333333333336e10) # Called "CIV_VALUE"
    td_totaldiscountedimpacts = Variable(unit="\$million")
    td_totaldiscountedimpacts_ts = Variable(index=[time], unit="\$million") # for analysis
    td_totaldiscountedimpacts_ann = Variable(unit="\$million")
    td_totaldiscountedimpacts_ann_yr = Variable(index=[year], unit="\$million") # for analysis

    aact_equityweightedadaptation_discountedaggregated = Variable(index=[time, region], unit="\$million")
    aact_equityweightedadaptation_discountedaggregated_ann = Variable(index=[year, region], unit="\$million")
    tac_totaladaptationcosts = Variable(unit="\$million")
    tac_totaladaptationcosts_ann = Variable(unit="\$million")

    # Final result: total effect of climate change
    te_totaleffect = Variable(unit="\$million")
    te_totaleffect_ann = Variable(unit="\$million")
    te_totaleffect_ann_yr = Variable(index=[year], unit="\$million") # for analysis

    function run_timestep(p, v, d, tt)

        # interpolate the parameters that require interpolation:
        interpolate_parameters_equityweighting(p, v, d, tt)

        if is_first(tt)
            v.tpc_totalaggregatedcosts = 0
            v.addt_gt_equityweightedimpact_discountedglobal = 0
            v.tac_totaladaptationcosts = 0
            v.te_totaleffect = 0
        end

        v.df_utilitydiscountfactor[tt] = (1 + p.ptp_timepreference / 100)^(-(p.y_year[tt] - p.y_year_0))

        for rr in d.region

            ## Gas Costs Accounting
            # Weighted costs (Page 23 of Hope 2009)
            v.wtct_percap_weightedcosts[tt, rr] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (p.cons_percap_consumption[tt, rr]^(1 - p.emuc_utilityconvexity) - (p.cons_percap_consumption[tt, rr] - p.tct_percap_totalcosts_total[tt, rr] < 0.01 * p.cons_percap_consumption_0[1] ? 0.01 * p.cons_percap_consumption_0[1] : p.cons_percap_consumption[tt, rr] - p.tct_percap_totalcosts_total[tt, rr])^(1 - p.emuc_utilityconvexity))

            # Add these into consumption
            v.eact_percap_weightedadaptationcosts[tt, rr] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (p.cons_percap_consumption[tt, rr]^(1 - p.emuc_utilityconvexity) - (p.cons_percap_consumption[tt, rr] - p.act_percap_adaptationcosts[tt, rr] < 0.01 * p.cons_percap_consumption_0[1] ? 0.01 * p.cons_percap_consumption_0[1] : p.cons_percap_consumption[tt, rr] - p.act_percap_adaptationcosts[tt, rr])^(1 - p.emuc_utilityconvexity))

            # Do partial weighting
            if p.equity_proportion == 0
                v.pct_percap_partiallyweighted[tt, rr] = p.tct_percap_totalcosts_total[tt, rr]
                v.wact_percap_partiallyweighted[tt, rr] = p.act_percap_adaptationcosts[tt, rr]
            else
                v.pct_percap_partiallyweighted[tt, rr] = (1 - p.equity_proportion) * p.tct_percap_totalcosts_total[tt, rr] + p.equity_proportion * v.wtct_percap_weightedcosts[tt, rr]
                v.wact_percap_partiallyweighted[tt, rr] = (1 - p.equity_proportion) * p.act_percap_adaptationcosts[tt, rr] + p.equity_proportion * v.eact_percap_weightedadaptationcosts[tt, rr]
            end

            v.pct_partiallyweighted[tt, rr] = v.pct_percap_partiallyweighted[tt, rr] * p.pop_population[tt, rr]
            v.wact_partiallyweighted[tt, rr] = v.wact_percap_partiallyweighted[tt, rr] * p.pop_population[tt, rr]

            # Discount rate calculations
            v.dr_discountrate[tt, rr] = p.ptp_timepreference + p.emuc_utilityconvexity * (p.grw_gdpgrowthrate[tt, rr] - p.popgrw_populationgrowth[tt, rr])
            if is_first(tt)
                v.yp_yearsperiod[TimestepIndex(1)] = p.y_year[TimestepIndex(1)] - p.y_year_0
            else
                v.yp_yearsperiod[tt] = p.y_year[tt] - p.y_year[tt - 1]
            end

            if is_first(tt)
                v.dfc_consumptiondiscountrate[TimestepIndex(1), rr] = (1 + v.dr_discountrate[TimestepIndex(1), rr] / 100)^(-v.yp_yearsperiod[TimestepIndex(1)])
            else
                v.dfc_consumptiondiscountrate[tt, rr] = v.dfc_consumptiondiscountrate[tt - 1, rr] * (1 + v.dr_discountrate[tt, rr] / 100)^(-v.yp_yearsperiod[tt])
            end

            # Discounted costs
            if p.equity_proportion == 0
                v.pcdt_partiallyweighted_discounted[tt, rr] = v.pct_partiallyweighted[tt, rr] * v.dfc_consumptiondiscountrate[tt, rr]
                v.wacdt_partiallyweighted_discounted[tt, rr] = p.act_adaptationcosts_total[tt, rr] * v.dfc_consumptiondiscountrate[tt, rr]
                v.wit_partiallyweighted[tt, rr] = (p.cons_percap_aftercosts[tt, rr] - p.rcons_percap_dis[tt, rr]) * p.pop_population[tt, rr] # equivalent to emuc = 0
                v.widt_partiallyweighted_discounted[tt, rr] = v.wit_partiallyweighted[tt, rr] * v.dfc_consumptiondiscountrate[tt] # apply Ramsey discounting
            else
                v.pcdt_partiallyweighted_discounted[tt, rr] = v.pct_partiallyweighted[tt, rr] * v.df_utilitydiscountfactor[tt]
                v.wacdt_partiallyweighted_discounted[tt, rr] = v.wact_partiallyweighted[tt, rr] * v.df_utilitydiscountfactor[tt]
                ## Equity weighted impacts (end of page 28, Hope 2009)
                v.wit_partiallyweighted[tt, rr] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (p.cons_percap_aftercosts[tt, rr]^(1 - p.emuc_utilityconvexity) - p.rcons_percap_dis[tt, rr]^(1 - p.emuc_utilityconvexity)) * p.pop_population[tt, rr]

                v.widt_partiallyweighted_discounted[tt, rr] = v.wit_partiallyweighted[tt, rr] * v.df_utilitydiscountfactor[tt]
            end

            v.pcdat_partiallyweighted_discountedaggregated[tt, rr] = v.pcdt_partiallyweighted_discounted[tt, rr] * p.yagg_periodspan[tt]

            v.addt_equityweightedimpact_discountedaggregated[tt, rr] = v.widt_partiallyweighted_discounted[tt, rr] * p.yagg_periodspan[tt]
            v.aact_equityweightedadaptation_discountedaggregated[tt, rr] = v.wacdt_partiallyweighted_discounted[tt, rr] * p.yagg_periodspan[tt]

            # calculate  for this specific year
            if is_first(tt)
                v.tpc_totalaggregatedcosts_ann = 0
                v.addt_gt_equityweightedimpact_discountedglobal_ann = 0
                v.tac_totaladaptationcosts_ann = 0
                v.te_totaleffect_ann = 0
                for annual_year = 2015:(gettime(tt))
                    calc_equityweighting(p, v, d, tt, annual_year, rr)
                end
            else
                for annual_year = (gettime(tt - 1) + 1):(gettime(tt))
                    calc_equityweighting(p, v, d, tt, annual_year, rr)
                end
            end
        end

        v.pct_g_partiallyweighted_global[tt] = sum(v.pct_partiallyweighted[tt, :])
        v.pcdt_g_partiallyweighted_discountedglobal[tt] = sum(v.pcdt_partiallyweighted_discounted[tt, :])
        v.tpc_totalaggregatedcosts = v.tpc_totalaggregatedcosts + sum(v.pcdat_partiallyweighted_discountedaggregated[tt, :])

        v.addt_gt_equityweightedimpact_discountedglobal = v.addt_gt_equityweightedimpact_discountedglobal + sum(v.addt_equityweightedimpact_discountedaggregated[tt, :])

        v.tac_totaladaptationcosts = v.tac_totaladaptationcosts + sum(v.aact_equityweightedadaptation_discountedaggregated[tt, :])

        # without  civilisation value
        # v.td_totaldiscountedimpacts = v.addt_gt_equityweightedimpact_discountedglobal
        # with civilisation value
        v.td_totaldiscountedimpacts = min(v.addt_gt_equityweightedimpact_discountedglobal, p.civvalue_civilizationvalue)

        v.td_totaldiscountedimpacts_ts[tt] = v.td_totaldiscountedimpacts # for analysis

        # Total effect of climate change
        # without civilisation value
        # v.te_totaleffect = v.td_totaldiscountedimpacts + v.tpc_totalaggregatedcosts + v.tac_totaladaptationcosts
        # with civilisation value
        v.te_totaleffect = min(v.td_totaldiscountedimpacts + v.tpc_totalaggregatedcosts + v.tac_totaladaptationcosts, p.civvalue_civilizationvalue)
    end
end


text:
using Weave
using NeXLMatrixCorrection

let start_dir = pwd()
    cd(@__DIR__)
    outpath = normpath(joinpath(@__DIR__, "..", "docs", "src"))
    @show outpath
    if !isdirpath(outpath)
        mkpath(outpath)
    end

    weave("example.jmd", out_path=joinpath(outpath,"example.md"), doctype="github")
    weave("coatingthickness.jmd", out_path=joinpath(outpath,"coatingthickness.md"), doctype="github")
    weave("testagainstpap.jmd", out_path=joinpath(outpath,"testagainstpap.md"), doctype="github")
    weave("testagainstheinrich.jmd", out_path=joinpath(outpath,"testagainstheinrich.md"), doctype="github")

    cd(start_dir)
end

ext: jl

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
############################################################################################
"""
$(TYPEDEF)
HMC - Container used throughout sampling process.

# Fields
$(TYPEDFIELDS)
"""
mutable struct HMC{
    R<:ℓObjectiveResult,D<:AbstractDifferentiableTune,C<:KineticEnergy,T<:BaytesCore.UpdateBool
} <: MCMCKernel
    "Target result stored for caching."
    result::R
    "Differentiation tuning struct."
    diff::D
    "Energy used for Hamiltonian."
    energy::C
    "Tuning struct for discretization steps"
    stepnumber::StepNumberTune{T}
    function HMC(
        result::R, diff::D, energy::C, stepnumber::StepNumberTune{T}
    ) where {
        R<:ℓObjectiveResult,
        D<:AbstractDifferentiableTune,
        C<:KineticEnergy,
        T<:BaytesCore.UpdateBool,
    }
        return new{R,D,C,T}(result, diff, energy, stepnumber)
    end
end
function HMC(sym::S; kwargs...) where {S<:Union{Symbol,NTuple{k,Symbol} where k}}
    return MCMC(HMC, sym; kwargs...)
end

function update!(kernel::HMC, objective::Objective, up::BaytesCore.UpdateTrue)
    ## Update log-target result with current (latent) data
    kernel.diff = update(kernel.diff, objective)
    kernel.result = ModelWrappers.log_density_and_gradient(objective, kernel.diff)
    return nothing
end
function update!(kernel::HMC, objective::Objective, up::BaytesCore.UpdateFalse)
    return nothing
end

############################################################################################
#Trajectory
"""
$(TYPEDEF)
Representation of a trajectory, ie a Hamiltonian with a discrete integrator that also checks for divergence.

# Fields
$(TYPEDFIELDS)
"""
struct TrajectoryHMC{S<:Hamiltonian,F<:AbstractFloat}
    "Hamiltonian."
    H::S
    "Log density of negative log energy at initial point."
    ℓH₀::F
    "Stepsize for leapfrog."
    ϵ::F
    function TrajectoryHMC(H::S, ℓH₀::F, ϵ::F) where {S<:Hamiltonian,F<:AbstractFloat}
        return new{S,F}(H, ℓH₀, ϵ)
    end
end

function move(trajectory::T, phasepoint::PhasePoint) where {T<:TrajectoryHMC}
    @unpack H, ϵ = trajectory
    return leapfrog(H, phasepoint, ϵ)
end

function checkfinite(trajectory::TrajectoryHMC, phasepoint::PhasePoint)
    return checkfinite(
        trajectory.ℓH₀, ℓdensity(trajectory.H, phasepoint), phasepoint.result
    )
end

############################################################################################
"""
$(SIGNATURES)
Propagate forward one HMC step.

# Examples
```julia
```

"""
function propagate(
    _rng::Random.AbstractRNG, kernel::HMC, tune::MCMCTune, objective::Objective
)
    @unpack ϵ = tune.stepsize
    ## Set new DiffObjective
    diff = DiffObjective(objective, kernel.diff)
    ## If needed, update kinetic energy and proposal leapfrog steps
    update!(kernel.energy, tune.proposal)
    update!(kernel.stepnumber, ϵ)
    ## Fill new Hamiltonian with current Energy and log-posterior adjusted for (latent/new) data
    H = Hamiltonian(kernel.energy, diff)
    ## Create new trajectory and phasepoint, and evaluate Hamiltonian at phasepoint
    phasepoint = PhasePoint(kernel.result, rand_ρ(_rng, H.K))
    ℓH₀ = ℓdensity(H, phasepoint)
    trajectory = TrajectoryHMC(H, ℓH₀, ϵ)
    ## leapfrog and check for divergence
    divergent = false
    for _ in Base.OneTo(kernel.stepnumber.steps)
        phasepoint = move(trajectory, phasepoint)
        if !checkfinite(trajectory, phasepoint)
            divergent = true
            break
        end
    end
    ## Pack container and return output
    accept_statistic = BaytesCore.AcceptStatistic(_rng, ℓdensity(H, phasepoint) - ℓH₀)
    return phasepoint.result,
    divergent, accept_statistic,
    DiagnosticsHMC(tune.stepsize.ϵ, kernel.stepnumber.steps)
end

############################################################################################
"""
$(SIGNATURES)
Tune initial step size.

# Examples
```julia
```

"""
function get_acceptrate(
    _rng::Random.AbstractRNG, kernel::HMC, objective::Objective, Σ::M
) where {M<:AbstractMatrix}
    ## Update Hamiltonian with current Σ and log-posterior adjusted for (latent/new) data
    H = Hamiltonian(kernel.energy, DiffObjective(objective, kernel.diff))
    ## Create new trajectory and phasepoint, and evaluate Hamiltonian at phasepoint
    phasepoint = PhasePoint(kernel.result, rand_ρ(_rng, H.K))
    H₀ = ℓdensity(H, phasepoint)
    ## Function of stepsize
    return function acceptrate(ϵ::T) where {T<:Real}
        # Make leapfrog step
        phasepointᵖ = leapfrog(H, phasepoint, ϵ)
        # Calculate acceptance rate
        return exp(ℓdensity(H, phasepointᵖ) - H₀)
    end
end

############################################################################################
#export
export HMC


text:
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
#                                           ANGULAR COEFFICIENTS                                          #
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
export Epsilon, Eta, mλ, RSphCA, CSphCA, SphCB, SphPCD, SphPCG0, SphPCG, ClebschGordanG, ClebschGordan,
GauntGF, GauntG, GGauntG, GauntWF, GauntW
###########################################################################################################
############################################ SIGN COEFFICIENTS ############################################
Epsilon(x, y) = Times(Nonnegative(x), Nonnegative(y))
Eta(x) = Power((-1), ((abs(x)-x)//2)) # Spherical spinors sign coefficient to convert to Bethe phase
Eta(x, y, z) = Sign(x, y, z)
mλ(m, λ) = m + λ - 1//2 # Nonrelativistic angular momentum quamtum number m from mj
###########################################################################################################
############################### SPHERICAL HARMONICS PRODUCT COEFFICIENTS A ################################(25)
################# A_{σ}(m, mü)
function RSphCA(σ::Int, m1::Int, m2::Int)
    res1 = Sqrt(2-Abs(Eta(m1, m2, m1-m2)))*KroneckerDelta(σ, Epsilon(m1, m2)*Abs(m1-m2))
    res2 = Eta(m1, m2, m1+m2)*KroneckerDelta(σ, Epsilon(m1, m2)*Abs(m1+m2))
    RSqrt(2) * (res1 + res2)
end
RSphCA(σ::arb, m1::arb, m2::arb) = RSphCA(
    Int(convert(Float64, σ)),
    Int(convert(Float64, m1)),
    Int(convert(Float64, m2))
)
RSphCA(σ::Real, m1::Real, m2::Real) = RSphCA(Int(σ), Int(m1), Int(m2))

CSphCA(σ::Int, m1::Int, m2::Int) = σ == m1 - m2 ? One(ArbField(sprec)) : Zero(ArbField(sprec))
CSphCA(σ::arb, m1::arb, m2::arb) = CSphCA(
    Int(convert(Float64, σ)),
    Int(convert(Float64, m1)),
    Int(convert(Float64, m2))
)
CSphCA(σ::Real, m1::Real, m2::Real) = CSphCA(Int(σ), Int(m1), Int(m2))
###########################################################################################################
################################### SPHERICAL HARMONICS COEFFICIENTS B ####################################(25)
################# B(k,l,m)
function SphCB(k::Int, l::Int, λ::Int)
    km2 = 2 * k
    lm2 = 2 * l

    res1 = (((-1)^k)//Power(2,l))*RSqrt(Binomial(l, λ) * Binomial(l + λ, λ))
    res2 = Binomial(k + λ, k) * Binomial(lm2 - km2, l - k) * Binomial(l - k, l - λ - km2)
    res1 * res2
end
SphCB(k::arb, l::arb, m::arb) = SphCB(
    Int(convert(Float64, k)),
    Int(convert(Float64, l)),
    Int(convert(Float64, m))
)
SphCB(k::Real, l::Real, λ::Real) = SphCB(Int(k), Int(l), Int(λ))
###########################################################################################################
############################### SPHERICAL HARMONICS PRODUCT COEFFICIENTS D ################################(25)
################# D(l,lambda,beta)
function SphPCD(l::Int, λ::Int, β::Int)
    lmβd2 = (l - β) / 2 # Exponentiation below line yielding a complex result when mod(l - b, 2) != 0
    #res1 = ((-1)^lmβd2) / (2^l)
    res1 = (EvenQ(lmβd2) ? 1 : -1) / (2^l)
    res2 = Sqrt((2 * l + 1) * Binomial(l + λ, l) / (2 * Binomial(l, λ)))
    res3 = Binomial(l, floor(Int, lmβd2)) * Binomial(l + β, β - λ)
    res1 * res2 * res3
end
SphPCD(l::arb, λ::arb, β::arb) = β < 0 ? 0 : SphPCD(
    UInt(convert(Float64, l)),
    UInt(convert(Float64, λ)),
    UInt(convert(Float64, β))
)
SphPCD(l::Real, λ::Real, β::Real) = β < 0 ? 0 : SphPCD(Int(l), Int(λ), Int(β))
###########################################################################################################
############################### SPHERICAL HARMONICS PRODUCT COEFFICIENTS G ################################(25)
################# G(alpha, beta, l1, lambda1, l2, lambda2, Lambda)
function SphPCG0(α::Int, β::Int, l1::Int, λ1::Int, l2::Int, λ2::Int, Λ::Int)
    res = zeros(RF, Λ + 1)
    for s in 0 : Λ
        res[s + 1] = ((-1)^(s)) * Binomial(Λ, s) * SphPCD(l1, λ1, α + 2 * Λ - 2 * s)
    end
    sum(res) * SphPCD(l2, λ2, β)
end
#################
function SphPCG(q::Int, α::Int, β::Int, l1::Int, λ1::Int, l2::Int, λ2::Int, Λ::Int)
    SphPCG0(α, β, l1, λ1, l2, λ2, Λ) * GBinomial(q, α + 2 * Λ - λ1, β - λ2)
end
###########################################################################################################
####################################### CLEBSCH-GORDAN COEFFICIENTS #######################################(32)
################# I. I. Guseinov formulation
function ClebschGordanG(l1::Int, m1::Int, l2::Int, m2::Int, L::Int, M::Int)
    if L < M || m1 + m2 != M
        return Zero(ArbField(sprec))
    else
        phs = (-1)^((m1 + abs(m1) + m2 + abs(m2) + M + abs(M)) / 2)
        Lt2 = 2 * L + 1
        l1t2 = 2 * l1 + 1
        l2t2 = 2 * l2 + 1

        lL1 = l1 + l2 + L + 1
        lL2 = l1 + l2 - L
        lL3 = l1 - l2 + L
        lL4 = l2 - l1 + L
        res1 = (Lt2^2)*Binomial(lL1, lL2)*Binomial(2*L, L+M)
        res2 = l1t2*l2t2*Binomial(lL1, lL3)*Binomial(lL1, lL4)*Binomial(2*l1, l1+m1)*Binomial(2*l2, l2+m2)

        lowl = max(0, l1-m1-L+M, l2+m2-L-M)
        upl = min(l1-m1, l2+m2, l1+l2-L)

        if lowl > upl
            lowl = upl = 0
        end

        res3 = zeros(RF, upl-lowl+1)
        a1 = zeros(RF, upl-lowl+1)
        a2 = zeros(RF, upl-lowl+1)
        a3 = zeros(RF, upl-lowl+1)
        a1[1] = Binomial(lL2,lowl)
        a2[1] = Binomial(L-M,l1-m1-lowl)
        a3[1] = Binomial(L+M,l2+m2-lowl)

        for s in lowl : upl-1
            s1 = s-lowl+1
            a1[s1+1] = ((lL2-s) // (s+1))*a1[s1]
            a2[s1+1] = ((l1-m1-s) // ((L-M)-(l1-m1-s)+1))*a2[s1]
            a3[s1+1] = ((l2+m2-s) // ((L+M)-(l2+m2-s)+1))*a3[s1]
        end

        for s in lowl : upl
            s1 = s-lowl+1
            res3[s1] = ((-1)^s)*a1[s1]*a2[s1]*a3[s1]
        end
        res = phs * KroneckerDelta(M,m1 + m2) * rsqrt(res2) * sqrt(res1) * sum(res3)
        return res
    end
end
################# Using Wigner3j symbol via WignerSymbol package
function ClebschGordan(l1, m1, l2, m2, L, M)
    return clebschgordan(BigFloat, l1, m1, l2, m2, L, M)
end
###########################################################################################################
########################################### GAUNT COEFFICIENTS ############################################(37)
################# If "," is used the given value is defauls. If ";" is used the value assigned.
function GauntGF(l1::Int, m1::Int, l2::Int, m2::Int, L::Int; M::Int = m1-m2)
    lowl = max(0, L-m1-l2)
    upl = min(l1-abs(m1), L-M, L-m1+l2)

    if EvenQ(l1 + l2 + L) == true && l1 - l2 <= L <= l1 + l2 && upl >= lowl
        γ = (l1 + l2 + L) // 2

        phs = (-1)^(γ - (l2-m2) + (abs(m1)+abs(m2)+abs(M)) // 2)
        res1 = (Binomial(2*γ-l1-l2, γ-l1)*Binomial(γ, L)) // ((2*γ+1)Binomial(2*γ, 2*L))
        res2 = (2*l1+1)*(2*l2+1)*Binomial(l1+l2+M, l1+m1)*Binomial(2*L+l1+l2+M, l1+l2+M)
        res3 = Binomial(l1+l2-M, l1-m1)*Binomial(2*L+l1+l2+M, l1+l2-M)*Binomial(2*L, L-M)*Binomial(2*L+2*M, L+M)
        res4 = Binomial(l1+l2+M, l1+l2-L) // Binomial(l1+l2+M, l1+m1)

        bin1 = zeros(RF, upl-lowl+1)
        bin2 = zeros(RF, upl-lowl+1)
        bin3 = zeros(RF, upl-lowl+1)
        bin4 = zeros(RF, upl-lowl+1)

        bin1[1] = Binomial(l1+m1+lowl,lowl)
        bin2[1] = Binomial(l1+l2-L, l1-m1-lowl)
        bin3[1] = Binomial(l2-m2+L-M-lowl, L-M-lowl)
        bin4[1] = ((-1)^lowl) * bin1[1] * bin2[1] * bin3[1]

        for s in lowl : upl-1
            s1 = s-lowl+1
            bin1[s1+1] = ((l1+m1+s+1) // (s+1))*bin1[s1]
            bin2[s1+1] = ((l1-m1-s) // (l2-L+m1+s+1))*bin2[s1]
            bin3[s1+1] = ((L-M-s) // (l2-m2+L-M-s))*bin3[s1]
            bin4[s1+1] = ((-1)^(s+1)) * bin1[s1+1] * bin2[s1+1] * bin3[s1+1]
        end
        res5 = res1 * (Sqrt(res2 // res3)) * res4
        res = phs * res5 * sum(bin4)
        return res
    else
        return Zero(ArbField(sprec))
    end
end

function GauntG(l1::Int, m1::Int, l2::Int, m2::Int, L::Int)
    M = m1-m2
    if l1 >= l2
        if m1 >= 0 && m2 >= 0
            GauntGF(l1, m1, l2, m2, L)
        elseif m1 < 0 && m2 >= 0
            GauntGF(l2, m2, l1, m1, L)
        elseif m1 >= 0 && m2 < 0
            Sqrt((2*l1+1) // (2*L+1))*GauntGF(L, M, l2, -m2, l1)
        elseif m1 < 0 && m2 < 0
            GauntGF(l1, -m1, l2, -m2, L)
        else
            Zero(ArbField(sprec))
        end
    else l1 < l2
        if m1 >= 0 && m2 >= 0 || m1 < 0 && m2 >= 0
            GauntGF(l2, m2, l1, m1, L)
        elseif  m1 >= 0 && m2 < 0
            Sqrt((2*l1+1) // (2*L+1))*GauntGF(L, M, l2, -m2, l1)
        elseif m1 < 0 && m2 < 0
            GauntGF(l1, -m1, l2, -m2, L)
        else
            Zero(ArbField(sprec))
        end
    end
end

function GGauntG(l1::Int, m1::Int, l2::Int, m2::Int, L::Int, M::Int)
    if abs(M) == abs(m1 - m2)
        GauntG(l1, m1, l2, m2, L)
    elseif abs(M) == abs(m1 + m2)
        GauntG(l1, m1, l2, -m2, L)
    else
        return Zero(ArbField(sprec))
    end
end
################# Using Wigner3j symbol via WignerSymbol package
function GauntWF(l1::Int, m1::Int, l2::Int, m2::Int, L::Int)
    #phs = (-1)^(((l1 + l2 + L) // 2) - (l2- m1))
    res1 = Sqrt((2*l1 + 1) * (2*l2 + 1)) // (2*L + 1)
    res2 = ClebschGordan(l1, m1, l2, -m2, L, m1 - m2)
    res3 = ClebschGordan(l1, 0, l2, 0, L, 0)
    return res1 * res2 * res3
end

function GauntW(l1::Int, m1::Int, l2::Int, m2::Int, L::Int)
    M = m1 - m2
    if l1 >= l2 && m1 >= 0 && m2 >= 0
        GauntWF(l1, m1, l2, m2, L)
    elseif l1 >= l2 && m1 < 0 && m2 >= 0 || l1 < l2 && m1 >= 0 && m2 >= 0 || l1 < l2 && m1 < 0 && m2 >= 0
        GauntWF(l2, m2, l1, m1, L)
    elseif l1 >= l2 && m1 >= 0 && m2 < 0 || l1 < l2 && m1 >= 0 && m2 < 0
        sqrt((2*l1+1) // (2*L+1))*GauntWF(L, M, l2, -m2, l1)
    elseif l1 >= l2 && m1 < 0 && m2 < 0  || l1 < l2 && m1 < 0 && m2 < 0
        GauntWF(l1, -m1, l2, -m2, L)
    else
        Zero(ArbField(sprec))
    end
end
#For Sage_math results (GGauntG = GauntW) * Sqrt((2 * L + 1) // (4 * Pi(RF)))


text:
using MixedModeDebugger

function foo(x)
    z = 0
    for ii in 1:10
        z+=ii
    end
    y = sum((x, z, 11))
    return y
end
@breakpoint sum((1,2,3))
@run_mixedmode(Main.foo(10))


text:
import Base: @pure

@pure function printf_impl(fmt::Type{PrintArgument}, len, i=1)
  assert_argument_number(len, 1)
  quote
    Base.print(args[$i])
  end
end

@pure assert_argument_number(len, required) =
  @assert len == required ARGN_ERR

@pure function printf_impl(fmt::Type{PrintValue{V}}, len, i=1) where V
  assert_argument_number(len, 0)
  let s = String(V)
    if s == ""
      quote end
    else
      quote
        Base.print($s)
      end
    end
  end
end

@pure function printf_impl(fmt::Type{PrintBoth{A,B}}, len, i=1) where {A,B}
  l1 = arglength(A)
  l2 = arglength(B)
  assert_argument_number(len, l1+l2)

  quote
    $(printf_impl(A, l1, i))
    $(printf_impl(B, l2, i+l1))
  end
end

@pure function printf_impl(fmt::Type{String}, len, i=1) where {A,B}
  :?
end

@pure function arglength(fmt)
  a = Nothing
  b = fmt
  result = 0
  while b <: PrintBoth
    a <: PrintArgument && (result += 1)
    a = b.parameters[1]
    b = b.parameters[2]
  end
  a <: PrintArgument && (result += 1)
  b <: PrintArgument && (result += 1)
  result
end


text:
include("basic_setup.jl")

using Plots
using Statistics
using LaTeXStrings

grid_size = 200
num_samples = 2^18
θs = range(0,1,length=3)
grid = range(0.0, 2.0, length=grid_size)

gen = get_gen(5.0, 5.0, 1.0, 0.1)
system, initial = Trajectories.generate_configuration(gen, duration=2.0)

p = plot(initial)
plot!(p, system.response)

resampled_samples = zeros((grid_size, num_samples, length(θs)))
for (i, θ) in enumerate(θs)
    println("θ=$θ")
    chain = Trajectories.chain(system, θ)
    samples, acceptance = Trajectories.generate_mcmc_samples(initial, chain, 2^8, num_samples)

    for (j, s) in enumerate(samples)
        for (k, x) in enumerate(s(grid))
            resampled_samples[k, j, i] = x[1]
        end
    end
end

x = grid
y = reshape(mean(resampled_samples, dims=2), (length(x), :))
yerror = reshape(std(resampled_samples, dims=2), (length(x), :))

pyplot()
using Plots.PlotMeasures
p1 = plot(x, y[:,end],
    ribbon=yerror[:,end],
    fillalpha=0.2,
    fill_z=1.0,
    # ylim=(70,130),
    linewidth=0,
    label=nothing,
    c=:viridis,
)
plot!(p1, initial, linewidth=0.9, label="orig. signal")
plot!(p1, x, y,
    fillalpha=0.4,
    line_z=θs',
    fill_z=θs',
    label=hcat(("θ=$θ" for θ in θs)...),
    linewidth=2.5,
    ylabel=L"S_t",
    xlabel=L"t",
    title="signals",
    legend=true,
    colorbar=false,
    c=:viridis,
)

p2 = plot(system.response,
    c=:blue,
    title="response",
    legend=:false,
    ylabel=L"X_t",
    # ylim=(70,130)
)

p = plot(p2, p1, layout=Plots.grid(2,1), size=(600,600))
pos = p.o.axes[3].get_position()
pos.y1 *= 0.5
p.o.axes[3].set_position(pos)
PyPlot.show(p.o)
p.o

savefig(p, projectdir("plots", "signal_example.pdf"))


text:
@defcomp climatesf6cycle begin
    # Global SF6 emissions in kt of SF6
    globsf6 = Parameter(index=[time])

    #Atmospheric SF6 concentrations
    acsf6 = Variable(index=[time])

    #SF6 pre industrial
    sf6pre = Parameter(default = 0.04)

    #SF6 decay
    sf6decay = Variable()

    lifesf6 = Parameter(default = 3200.0)

    function run_timestep(p, v, d, t)

        if is_first(t)
            v.sf6decay = 1.0 / p.lifesf6

            v.acsf6[t] = p.sf6pre
        else
            # Calculate SF6 concentrations
            v.acsf6[t] = p.sf6pre + (v.acsf6[t - 1] - p.sf6pre) * (1 - v.sf6decay) + p.globsf6[t] / 25.1

            if v.acsf6[t] < 0
                error("sf6 atmospheric concentration out of range")
            end
        end
    end
end


text:
export GraphSpace

#######################################################################################
# Basic space definition
#######################################################################################
struct GraphSpace{G} <: DiscreteSpace
    graph::G
    s::Vector{Vector{Int}}
end

"""
    GraphSpace(graph::AbstractGraph)
Create a `GraphSpace` instance that is underlined by an arbitrary graph from
[LightGraphs.jl](https://github.com/JuliaGraphs/LightGraphs.jl).
The position type for this space is `Int`, use [`GraphAgent`](@ref) for convenience.
The underlying graph can be altered using [`add_node!`](@ref) and [`rem_node!`](@ref).

`GraphSpace` represents a space where each node (i.e. position) of a graph can hold an arbitrary
amount of agents, and each agent can move between the nodes of the graph.
An example of its usage can be found in [SIR model for the spread of COVID-19](@ref).
If you want to model social networks, where each agent is equivalent with a node of
a graph, you're better of using `nothing` (or other spaces) as the model space, and using
a graph from LightGraphs.jl directly in the model parameters, as shown in the
[Social networks with LightGraphs.jl](@ref) integration example.
"""
function GraphSpace(graph::G) where {G<:AbstractGraph}
    agent_positions = [Int[] for i in 1:LightGraphs.nv(graph)]
    return GraphSpace{G}(graph, agent_positions)
end

function Base.show(io::IO, s::GraphSpace)
    print(io, "GraphSpace with $(nv(s.graph)) positions and $(ne(s.graph)) edges")
end

#######################################################################################
# Agents.jl space API
#######################################################################################
random_position(model::ABM{<:GraphSpace}) = rand(model.rng, 1:nv(model))

function remove_agent_from_space!(
    agent::A,
    model::ABM{<:GraphSpace,A},
) where {A<:AbstractAgent}
    agentpos = agent.pos
    ids = ids_in_position(agentpos, model)
    splice!(ids, findfirst(a -> a == agent.id, ids))
    return model
end

function move_agent!(
    agent::A,
    pos::ValidPos,
    model::ABM{<:GraphSpace,A},
) where {A<:AbstractAgent}
    oldpos = agent.pos
    ids = ids_in_position(oldpos, model)
    splice!(ids, findfirst(a -> a == agent.id, ids))
    agent.pos = pos
    push!(ids_in_position(agent.pos, model), agent.id)
    return agent
end

function add_agent_to_space!(
    agent::A,
    model::ABM{<:DiscreteSpace,A},
) where {A<:AbstractAgent}
    push!(ids_in_position(agent.pos, model), agent.id)
    return agent
end

# The following two is for the discrete space API:
ids_in_position(n::Integer, model::ABM{<:GraphSpace}) = model.space.s[n]
# NOTICE: The return type of `ids_in_position` must support `length` and `isempty`!

#######################################################################################
# Neighbors
#######################################################################################
function nearby_ids(pos::Int, model::ABM{<:GraphSpace}, args...; kwargs...)
    np = nearby_positions(pos, model, args...; kwargs...)
    # This call is faster than reduce(vcat, ..), or Iterators.flatten
    vcat(model.space.s[pos], model.space.s[np]...)
end

function nearby_ids(agent::A, model::ABM{<:GraphSpace,A}, args...; kwargs...) where {A<:AbstractAgent}
    all = nearby_ids(agent.pos, model, args...; kwargs...)
    filter!(i -> i ≠ agent.id, all)
end

function nearby_positions(
    position::Integer,
    model::ABM{<:GraphSpace};
    neighbor_type::Symbol = :default,
)
    @assert neighbor_type ∈ (:default, :all, :in, :out)
    neighborfn = if neighbor_type == :default
        LightGraphs.neighbors
    elseif neighbor_type == :in
        LightGraphs.inneighbors
    elseif neighbor_type == :out
        LightGraphs.outneighbors
    else
        LightGraphs.all_neighbors
    end
    neighborfn(model.space.graph, position)
end

#######################################################################################
# Mutable graph functions
#######################################################################################
export rem_node!, add_node!, add_edge!

"""
    rem_node!(model::ABM{<: GraphSpace}, n::Int)
Remove node (i.e. position) `n` from the model's graph. All agents in that node are killed.

**Warning:** LightGraphs.jl (and thus Agents.jl) swaps the index of the last node with
that of the one to be removed, while every other node remains as is. This means that
when doing `rem_node!(n, model)` the last node becomes the `n`-th node while the previous
`n`-th node (and all its edges and agents) are deleted.
"""
function rem_node!(model::ABM{<:GraphSpace}, n::Int)
    for id in copy(ids_in_position(n, model))
        kill_agent!(model[id], model)
    end
    V = nv(model)
    success = LightGraphs.rem_vertex!(model.space.graph, n)
    n > V && error("Node number exceeds amount of nodes in graph!")
    s = model.space.s
    s[V], s[n] = s[n], s[V]
    pop!(s)
end

"""
    add_node!(model::ABM{<: GraphSpace})
Add a new node (i.e. possible position) to the model's graph and return it.
You can connect this new node with existing ones using [`add_edge!`](@ref).
"""
function add_node!(model::ABM{<:GraphSpace})
    add_vertex!(model.space.graph)
    push!(model.space.s, Int[])
    return nv(model)
end

"""
    add_edge!(model::ABM{<: GraphSpace}, n::Int, m::Int)
Add a new edge (relationship between two positions) to the graph.
Returns a boolean, true if the operation was succesful.
"""
LightGraphs.add_edge!(model, n, m) = add_edge!(model.space.graph, n, m)


text:
###

println("High level method write_number()...")

data = [0,1,2,3,4,5,6,7]
write_number(d, data)
println("$data\n")
sleep(5)

data = [8,9,10,11,12,13,14,15]
write_number(d, data)
println("$data\n")
sleep(5)

data = [6,6,6,nothing,9,9,9,nothing]
write_number(d, data)
println("$data\n")
sleep(5)

write_number(d, data, 1)
println("$data with dot at 1\n")
sleep(3)

write_number(d, data, 2)
println("$data with dot at 2\n")
sleep(3)

write_number(d, data, 3)
println("$data with dot at 3\n")
sleep(3)

write_number(d, data, 4)
println("$data with dot at 4\n")
sleep(3)

data = 31415926
write_number(d, data, 8)
println("$data with dot at 8\n")
sleep(5)

println("High level method write_symbols()...")

data = [0b01001001, 0b01001001, 0b01001001, 0b01001001, 0b01001001, 0b01001001, 0b01001001, 0b01001001]
write_symbols(d, data)
println("$data\n")
sleep(5)

data = ['A', 'I', 'L', 'U', 'J']
write_symbols(d, data)
println("$data\n")
sleep(5)

data = "Hello Hi"
write_symbols(d, data)
println("$data\n")
sleep(5)


text:
using DecomposedMDPSolver
using Flux
using Test

Na = 5
Np = 10
solutions = [(x) -> rand(Na) for i=1:Np ]

## Weights network
base = Chain(Dense(4, 32, relu), Dense(32, Na))
attn = Chain(Dense(4, 32, relu), Dense(32, Np+1), softmax)
a2t_model = A2TNetwork(base, attn, solutions)
@test size(a2t_model(rand(4))) == (Na, 1)

## Constant Weights
base = Chain(Dense(4, 32, relu), Dense(32, Na))
attn = Chain(ConstantLayer(Np+1), softmax)
a2t_model = A2TNetwork(base, attn, solutions)
@test size(a2t_model(rand(4))) == (Na, 1)

## Version that was failing -- local approx policy eval

sols = rand(100)
function val(s)
    sols[1] = rand()
    return sols[1]
end
base = Chain(Dense(2, 32, relu), Dense(32, 1, σ))
attn = Chain(Dense(2, 32, relu), Dense(32, 2, exp))
solutions = [val]

model = A2TNetwork(base, attn, solutions)

S, G = rand(2, 100), rand(1,100)
data = Flux.Data.DataLoader(S, G, batchsize=32, shuffle = true)
opt = ADAM()
Flux.train!((x, y) -> Flux.mse(model(x), y), Flux.params(model), data, opt)


text:
using AbstractPlotting
using MakieGallery
using Test
using StaticArrays, GeometryBasics
using Observables
using GeometryBasics: Pyramid
using PlotUtils
using MeshIO, FileIO, AbstractPlotting.MakieLayout

# Minimal sanity checks for MakieLayout
@testset "Layoutables constructors" begin
    scene, layout = layoutscene()
    ax = layout[1, 1] = LAxis(scene)
    cb = layout[1, 2] = LColorbar(scene)
    gl2 = layout[2, :] = MakieLayout.GridLayout()
    bu = gl2[1, 1] = LButton(scene)
    sl = gl2[1, 2] = LSlider(scene)

    scat = scatter!(ax, rand(10))
    le = gl2[1, 3] = LLegend(scene, [scat], ["scatter"])

    to = gl2[1, 4] = LToggle(scene)
    te = layout[0, :] = LText(scene, "A super title")
    me = layout[end + 1, :] = LMenu(scene, options=["one", "two", "three"])
    tb = layout[end + 1, :] = LTextbox(scene)
    @test true
end

@testset "Minimal AbstractPlotting tests" begin

    include("conversions.jl")
    include("quaternions.jl")
    include("projection_math.jl")
    include("shorthands.jl")
    include("liftmacro.jl")

    @testset "basic functionality" begin
        scene = scatter(rand(4))
        @test scene[Axis].ticks.title_gap[] == 3
        scene[Axis].ticks.title_gap = 4
        @test scene[Axis].ticks.title_gap[] == 4
        @test scene[Axis].tickmarks.length[] == (1, 1)

        scene = scatter([22.0, 28.0])
        AbstractPlotting.update_limits!(scene)
        lim = scene.data_limits[]
        @test lim.origin[1] <= 1 && lim.widths[1] >= 1 && lim.origin[1]+lim.widths[1] >= 2
        @test lim.origin[2] <= 22 && lim.widths[2] >= 6 && lim.origin[2]+lim.widths[2] >= 28
    end
end

# run software only tests...
include("no_backend_tests.jl")

# test statistical recipes
include("statistical_tests.jl")


text:
@testset "WreathProducts" begin
   S_3 = PermutationGroup(3)
   S_2 = PermutationGroup(2)
   b = perm"(1,2,3)"
   a = perm"(1,2)"

   @testset "Constructors" begin
      @test Groups.WreathProduct(S_2, S_3) isa AbstractAlgebra.Group
      B3 = Groups.WreathProduct(S_2, S_3)
      @test B3 isa Groups.WreathProduct
      @test B3 isa WreathProduct{3, Generic.PermGroup{Int}, Generic.PermGroup{Int}}

      aa = Groups.DirectPowerGroupElem((a^0 ,a, a^2))

      @test Groups.WreathProductElem(aa, b) isa AbstractAlgebra.GroupElem
      x = Groups.WreathProductElem(aa, b)
      @test x isa Groups.WreathProductElem
      @test x isa
         Groups.WreathProductElem{3, perm{Int}, perm{Int}}

      @test B3.N == Groups.DirectPowerGroup(S_2, 3)
      @test B3.P == S_3

      @test B3(aa, b) == Groups.WreathProductElem(aa, b)
      w = B3(aa, b)
      @test B3(w) == w
      @test B3(b) == Groups.WreathProductElem(B3.N(), b)
      @test B3(aa) == Groups.WreathProductElem(aa, S_3())

      @test B3((a^0 ,a, a^2), b) isa WreathProductElem

      @test B3((a^0 ,a, a^2), b) == B3(aa, b)
   end

   @testset "Types" begin
      B3 = Groups.WreathProduct(S_2, S_3)

      @test elem_type(B3) == Groups.WreathProductElem{3, perm{Int}, perm{Int}}

      @test parent_type(typeof(B3())) == Groups.WreathProduct{3, parent_type(typeof(B3.N.group())), Generic.PermGroup{Int}}

      @test parent(B3()) == Groups.WreathProduct(S_2,S_3)
      @test parent(B3()) == B3
   end

   @testset "Basic operations on WreathProductElem" begin
      aa = Groups.DirectPowerGroupElem((a^0 ,a, a^2))
      B3 = Groups.WreathProduct(S_2, S_3)
      g = B3(aa, b)

      @test g.p == b
      @test g.n == DirectPowerGroupElem(aa.elts)

      h = deepcopy(g)
      @test h == g
      @test !(g === h)

      g = B3(Groups.DirectPowerGroupElem((a ,a, a^2)), g.p)

      @test g.n[1] == parent(g.n[1])(a)
      @test g != h

      @test hash(g) != hash(h)
   end

   @testset "Group arithmetic" begin
      B4 = Groups.WreathProduct(AdditiveGroup(GF(3)), PermutationGroup(4))

      x = B4((0,1,2,0), perm"(1,2,3)(4)")
      @test inv(x) == B4((1,0,2,0), perm"(1,3,2)(4)")

      y = B4((1,0,1,2), perm"(1,4)(2,3)")
      @test inv(y) == B4((1,2,0,2), perm"(1,4)(2,3)")

      @test x*y == B4((0,2,0,2), perm"(1,3,4)(2)")

      @test y*x == B4((1,2,2,2), perm"(1,4,2)(3)")

      @test inv(x)*y == B4((2,1,2,2), perm"(1,2,4)(3)")

      @test y*inv(x) == B4((1,2,1,0), perm"(1,4,3)(2)")

      @test (x*y)^6 == ((x*y)^2)^3

   end

   @testset "Iteration" begin
      B3_a = Groups.WreathProduct(AdditiveGroup(GF(3)), S_3)
      @test order(B3_a) == 3^3*6
      @test collect(B3_a) isa Vector{
      WreathProductElem{3, AddGrpElem{AbstractAlgebra.gfelem{Int}}, perm{Int}}}

      B3_m = Groups.WreathProduct(MultiplicativeGroup(GF(3)), S_3)
      @test order(B3_m) == 2^3*6
      @test collect(B3_m) isa Vector{
      WreathProductElem{3, MltGrpElem{AbstractAlgebra.gfelem{Int}}, perm{Int}}}

      @test length(Set([B3_a, B3_m, B3_a])) == 2

      Wr = WreathProduct(PermutationGroup(2),PermutationGroup(4))

      elts = collect(Wr)
      @test elts isa Vector{Groups.WreathProductElem{4, perm{Int}, perm{Int}}}
      @test order(Wr) == 2^4*factorial(4)

      @test length(elts) == order(Wr)
      @test all([g*inv(g) == Wr() for g in elts])
      @test all(inv(g*h) == inv(h)*inv(g) for g in elts for h in elts)
   end

   @testset "Misc" begin
      B3_a = Groups.WreathProduct(AdditiveGroup(GF(3)), S_3)
      @test string(B3_a) == "Wreath Product of The additive group of Finite field F_3 by Permutation group over 3 elements"

      @test string(B3_a(perm"(1,3)")) == "([0,0,0]≀(1,3))"

      B3_m = Groups.WreathProduct(MultiplicativeGroup(GF(3)), S_3)
      @test string(B3_m) == "Wreath Product of The multiplicative group of Finite field F_3 by Permutation group over 3 elements"

      @test string(B3_m(perm"(1,3)")) == "([1,1,1]≀(1,3))"

   end

end


text:
# Module-related types and auxiliary functions

export
    CuModule, CuModuleFile

include(joinpath("module", "jit.jl"))

"""
    CuModule(data, options::Dict{CUjit_option,Any})
    CuModuleFile(path, options::Dict{CUjit_option,Any})

Create a CUDA module from a data, or a file containing data. The data may be PTX code, a
CUBIN, or a FATBIN.

The `options` is an optional dictionary of JIT options and their respective value.
"""
mutable struct CuModule
    handle::CUmodule
    ctx::CuContext

    function CuModule(data, options::Dict{CUjit_option,Any}=Dict{CUjit_option,Any}())
        handle_ref = Ref{CUmodule}()

        options[JIT_ERROR_LOG_BUFFER] = Vector{UInt8}(undef, 1024*1024)
        @debug begin
            options[JIT_INFO_LOG_BUFFER] = Vector{UInt8}(undef, 1024*1024)
            options[JIT_LOG_VERBOSE] = true
            "JIT compiling code" # FIXME: remove this useless message
        end
        optionKeys, optionVals = encode(options)

        # FIXME: maybe all CUDA API calls need to run under retry_reclaim?
        #        that would require a redesign of the memory pool,
        #        so maybe do so when we replace it with CUDA 11.2's pool.
        res = GC.@preserve data retry_reclaim(isequal(ERROR_OUT_OF_MEMORY)) do
            unsafe_cuModuleLoadDataEx(handle_ref, pointer(data),
                                      length(optionKeys),
                                      optionKeys, optionVals)
        end
        if res == ERROR_NO_BINARY_FOR_GPU ||
           res == ERROR_INVALID_IMAGE ||
           res == ERROR_INVALID_PTX
            options = decode(optionKeys, optionVals)
            throw(CuError(res, options[JIT_ERROR_LOG_BUFFER]))
        elseif res != SUCCESS
            throw_api_error(res)
        end

        @debug begin
            options = decode(optionKeys, optionVals)
            if isempty(options[JIT_INFO_LOG_BUFFER])
                """JIT info log is empty"""
            else
                """JIT info log:
                   $(options[JIT_INFO_LOG_BUFFER])"""
            end
        end

        ctx = CuCurrentContext()
        obj = new(handle_ref[], ctx)
        finalizer(unsafe_unload!, obj)
        return obj
    end
end

function unsafe_unload!(mod::CuModule)
    if isvalid(mod.ctx)
        cuModuleUnload(mod)
    end
end

Base.unsafe_convert(::Type{CUmodule}, mod::CuModule) = mod.handle

Base.:(==)(a::CuModule, b::CuModule) = a.handle == b.handle
Base.hash(mod::CuModule, h::UInt) = hash(mod.handle, h)

CuModuleFile(path) = CuModule(read(path, String))

include(joinpath("module", "function.jl"))
include(joinpath("module", "global.jl"))
include(joinpath("module", "linker.jl"))


text:
using OCToolbox
using Frechet_Evolution

# should probably test ODE evolution, non ODE Evolution
# TODO work out a good test case, Rabi?

#
# H_d = [1 0; 0 1]
# H_ti = nothing
# H_pc = nothing
# time_steps = 20
# tau = 0.1
#
# muh = Frechet_Time_Evolution(H_d)


text:
using CSyntax

@static if Sys.isapple()
    const VERSION_MAJOR = 4
    const VERSION_MINOR = 1
end

include(joinpath(@__DIR__, "glutils.jl"))

# init window
width, height = fb_width, fb_height = 640, 480
window = startgl(width, height)

glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)

# vertex data
points = GLfloat[ 0.0, 0.75, 0.0,
                  0.5, 0.25, 0.0,
                 -0.5, 0.25, 0.0,
                  0.5,-0.25, 0.0,
                  0.0,-0.75, 0.0,
                 -0.5,-0.25, 0.0]

# create buffers located in the memory of graphic card
points_vbo = GLuint(0)
@c glGenBuffers(1, &points_vbo)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(points), points, GL_STATIC_DRAW)

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)

# load and compile shaders from file
vert_shader = createshader(joinpath(@__DIR__, "tess.vert"), GL_VERTEX_SHADER)
tesc_shader = createshader(joinpath(@__DIR__, "tess.tesc"), GL_TESS_CONTROL_SHADER)
tese_shader = createshader(joinpath(@__DIR__, "tess.tese"), GL_TESS_EVALUATION_SHADER)
frag_shader = createshader(joinpath(@__DIR__, "tess.frag"), GL_FRAGMENT_SHADER)

# link program
shader_prog = createprogram(vert_shader, tesc_shader, tese_shader, frag_shader)

inner_tess_factor_loc = glGetUniformLocation(shader_prog, "tess_fac_inner")
outer_tess_factor_loc = glGetUniformLocation(shader_prog, "tess_fac_outer")

glEnable(GL_CULL_FACE)
glCullFace(GL_BACK)
glFrontFace(GL_CW)
glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)
glPatchParameteri(GL_PATCH_VERTICES, 3)

# set background color to gray
glClearColor(0.2, 0.2, 0.2, 1.0)

# render
let
inner_tess_factor = GLfloat(1)
outer_tess_factor = GLfloat(4)
qWasDown = false
aWasDown = false
wWasDown = false
sWasDown = false
updatefps = FPSCounter()
while !GLFW.WindowShouldClose(window)
    updatefps(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, GLFW.GetFramebufferSize(window)...)
    # drawing
    glUseProgram(shader_prog)
    glBindVertexArray(vao)
    glDrawArrays(GL_PATCHES, 0, 6)
    # check and call events
    GLFW.PollEvents()
    if GLFW.GetKey(window, GLFW.KEY_Q)
        if !qWasDown
            inner_tess_factor += 1
            println("inner tessellation factor = $inner_tess_factor.")
            glUniform1f(inner_tess_factor_loc, inner_tess_factor)
            qWasDown = true
        end
    else
        qWasDown = false
    end
    if GLFW.GetKey(window, GLFW.KEY_A)
        if !aWasDown
            inner_tess_factor -= 1
            println("inner tessellation factor = $inner_tess_factor.")
            glUniform1f(inner_tess_factor_loc, inner_tess_factor)
            aWasDown = true
        end
    else
        aWasDown = false
    end
    if GLFW.GetKey(window, GLFW.KEY_W)
        if !wWasDown
            outer_tess_factor += 1
            println("outer tessellation factor = $outer_tess_factor.")
            glUniform1f(outer_tess_factor_loc, outer_tess_factor)
            wWasDown = true
        end
    else
        wWasDown = false
    end
    if GLFW.GetKey(window, GLFW.KEY_S)
        if !sWasDown
            outer_tess_factor -= 1
            println("outer tessellation factor = $outer_tess_factor.")
            glUniform1f(outer_tess_factor_loc, outer_tess_factor)
            sWasDown = true
        end
    else
        sWasDown = false
    end
    # swap the buffers
    GLFW.SwapBuffers(window)
end
end # let
GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
function inex02()
#
#  test polynomial suggested by Goedecker
#
    inex(6)
end


text:
abstract type AbstractMaterial end

@inline function scatter(mat::AbstractMaterial, ray_in::Ray{T}, hitrecord::HitRecord{T}, rng::MersenneTwister) where {T<:Real}
     error("Scatter for $(typeof(mat)) has not been defined.")
end

struct Lambertian{T<:Real} <: AbstractMaterial
    albedo::Vec{3, T}
end

albedo(lambertian::Lambertian{T}) where {T<:Real} = lambertian.albedo

@inline function scatter(lambertian::Lambertian{T}, ray_in::Ray{T}, hitrecord::HitRecord{T}, rng::MersenneTwister) where {T<:Real}
    color = albedo(lambertian)
    facing_normal = normal(hitrecord)
    new_dir = normalize(facing_normal + random_unit_vector3(rng, T))
    ray_out = Ray(point(hitrecord), new_dir) # direction could be the normal?
    absorbed = false
    (color=color, ray_out=ray_out, absorbed=absorbed)
end

struct NormalShading <: AbstractMaterial
end

@inline function scatter(::NormalShading, ray_in::Ray{T}, hitrecord::HitRecord{T}, rng::MersenneTwister) where {T<:Real}
    color = (((normal(hitrecord) .* (front_face(hitrecord) ? 1 : -1)) .+ 1) ./ 2)
    ray_out = Ray(point(hitrecord), direction(ray_in)) # direction could be the normal?
    absorbed = true
    (color=color, ray_out=ray_out, absorbed=absorbed)
end


text:
using Bridge, LinearAlgebra, Test, Statistics
x = [rand(5) for i in 1:10]
m, sc, m2 = Bridge.iterateall(Bridge.MeanCov(x))
@test norm(cov(x, corrected=true) - sc*m2) < eps(100.0)
@test norm(mean(x) - m) < eps(100.0)

y = [rand() for i in 1:10]
m, sc, m2 = Bridge.iterateall(Bridge.MeanCov(y))
@test norm(var(y, corrected=true) - sc*m2) < eps(100.0)
@test norm(mean(y) - m) < eps(100.0)

@inferred Bridge.iterateall(Bridge.MeanCov(1:1000))
@inferred Bridge.iterateall(Bridge.MeanCov(x))


text:
using FileIO
using Printf
include("config.jl")
include("calc_loader.jl")
include("utils/utils.jl")

function save_result(X, Xr, runtime, method, reqrank, cnt, fix)
    tensor_size = size(X)
    LSerror = norm(X-Xr)
    if method in non_negative_methods
        KLerror = get_dkl(X, Xr)
    else
        KLerror = missing
    end
    result = Dict(
                  "tensor_size" => tensor_size,
                  "runtime" => runtime,
                  "LSerror" => LSerror,
                  "KLerror" => KLerror,
                  "method" => method,
                  "reqrank" => reqrank
                 )
    size_str =  get_tensorinfo(tensor_size)
    reqrank_str =  get_tensorinfo(reqrank)
    if fix == "size"
        savepath = "../result/synthetic/size$size_str/reqrank$reqrank_str"
    elseif fix == "rank"
        savepath = "../result/synthetic/reqrank$reqrank_str/size$size_str"
    else
        error("improper $sz")
    end
    mkpath(savepath)
    save("$savepath/$method$cnt.jld2", result)
    return LSerror, KLerror
end

function get_random_tensor( tensor_size, random_tenosr_dist)
    if random_tenosr_dist == "uniform"
        X = rand(tensor_size...)
    elseif random_tenosr_dist == "gauss"
        X = abs.(randn(tensor_size...))
    end
    return X
end

function run_fix_size(tensor_size, reqranks)
    X = get_random_tensor(tensor_size, random_tenosr_dist)
    println("input tensor shape $tensor_size")
    println("__________________________________________")
    @printf("%5s %20s %8s %10s %10s %10s \n", "cnt", "reqrank", "method", "runtime", "LSerror", "KLerror")
    println("__________________________________________")
    for reqrank in reqranks
        for method in methods
            for cnt = 1:rep_max
                runtime = @elapsed begin
                    Xr = calc(X, method, reqrank)
                end
                LSerror, KLerror = save_result(X, Xr, runtime, method, reqrank, cnt, "size")
                @printf("%5s %20s %8s %10f %10f %10s \n", cnt, reqrank, method, runtime, LSerror, KLerror)
            end
        end
        println("__________________________________________")
    end
end

function run_fix_rank(tensor_sizes, reqrank)
    println("fixed reqrank $reqrank")
    println("__________________________________________")
    @printf("%5s %20s %8s %10s %10s %10s \n", "cnt", "tensor_size", "method", "runtime", "LSerror", "KLerror")
    println("__________________________________________")
    for tensor_size in tensor_sizes
        X = get_random_tensor(tensor_size, random_tenosr_dist)
        for method in methods
            for cnt = 1:rep_max
                runtime = @elapsed begin
                    Xr = calc(X, method, reqrank)
                end
                LSerror, KLerror = save_result(X, Xr, runtime, method, reqrank, cnt, "rank")
                @printf("%5s %20s %8s %10f %10f %10s \n", cnt, tensor_size, method, runtime, LSerror, KLerror)
            end
        end
        println("__________________________________________")
    end
end
d = 5
tensor_size = [30 for i=1:d]
reqranks = []
for n = 2:1:15
     reqrank = [n for i=1:d]
     push!(reqranks, reqrank)
end
run_fix_size(tensor_size, reqranks)

reqrank = [10,10,10]
tensor_sizes = [[30,30,30], [50,50,50], [60,60,60], [80,80,80],
               [100,100,100], [120,120,120], [150,150,150], [200,200,200], [300,300,300], [400,400,400], [500,500,500]]
run_fix_rank(tensor_sizes, reqrank)

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
using Accessors, Documenter, Literate

inputdir = joinpath(@__DIR__, "..", "examples")
outputdir = joinpath(@__DIR__, "src", "examples")
mkpath(outputdir)
mkpath(joinpath(outputdir, "examples"))
for filename in readdir(inputdir)
    inpath = joinpath(inputdir, filename)
    cp(inpath, joinpath(outputdir, "examples", filename), force = true)
    Literate.markdown(inpath, outputdir; documenter = true)
end
cp(joinpath(@__DIR__, "..", "README.md"), joinpath(@__DIR__, "src", "index.md"))

makedocs(
    modules = [Accessors],
    sitename = "Accessors.jl",
    pages = [
        "Home" => "index.md",
        "Tutorials" => ["Getting started" => "getting_started.md",],
        "Explanation" => ["Lenses" => "lenses.md",],
        "Reference" => ["Docstrings" => "docstrings.md"],
        "How-to guides" => [
            "Custom Optics" => "examples/custom_optics.md",
            "Custom Macros" => "examples/custom_macros.md",
        ],
        hide("examples/specter.md"),
        hide("internals.md"),
    ],
    strict = true,  # to exit with non-zero code on error
)

deploydocs(;
    repo = "github.com/JuliaObjects/Accessors.jl.git",
    push_preview=true
)


text:
using Test
@info "The SimpleRationalFunction type is now part of the SimplePolynomials package."
@test true


text:
using LSODA

function rhs!(t, x, ydot, data)
	ydot[1]=1.0E4 * x[2] * x[3] - .04E0 * x[1]
	ydot[3]=3.0E7 * x[2] * x[2]
	ydot[2]=-ydot[1] - ydot[3]
  nothing
end

y0 = [1.,0.,0.]
println("\n####################################\n--> Use of a old wrapper for speed comparison")
tspan = [0., 0.4]
@time LSODA.lsoda_0(rhs!, y0, tspan, reltol= 1e-4,abstol = Vector([1.e-6,1.e-10,1.e-6]))
@time LSODA.lsoda_0(rhs!, y0, tspan, reltol= 1e-4,abstol = Vector([1.e-6,1.e-10,1.e-6]))

# case with a vector
println("\n####################################\n--> Use of a vector of times where output is required")
tspan = [4 .* 10.0^k for k=-1:10]
res = LSODA.lsoda(rhs!, y0, tspan, reltol= 1e-4,abstol = Vector([1.e-6,1.e-10,1.e-6]))
res = @time LSODA.lsoda(rhs!, y0, tspan, reltol= 1e-4,abstol = Vector([1.e-6,1.e-10,1.e-6]))
println(res)

# #case where we don't have to declare a new context
println("\n####################################\n--> Use of a lsoda_evolve!")
# y0 = [1.,0.,0.]
# tspan = [4.*10.0^k for k=-1:10]
# ctx, _ = LSODA.lsoda(rhs!, y0, tspan[1:2], reltol= 1e-4,abstol = Vector([1.e-6,1.e-10,1.e-6]))
# @time for k=2:length(tspan)
# 	LSODA.lsoda_evolve!(ctx,y0,tspan[k-1:k])
# 	@printf("at t = %12.4e y= %14.6e %14.6e %14.6e\n",tspan[k],y0[1], y0[2], y0[3])
# end
# lsoda_free(ctx)


text:
using DelimitedFiles

function read_data(
                    problem::Int64,
                    filename::String
                    )

    # Reads the execution data
    data = readdlm(filename)

    # Print data
    println("Dimension: $(data[problem, 1])")
    println("N. functions: $(data[problem, 2])")
    println("Status: $(data[problem, 3])")
    println("True value: $(data[problem, 4])")
    println("N. iterations: $(data[problem, 5])")
    println("N. evaluations: $(data[problem, 6])")
    println("N. eval. / N. functions: $(data[problem, 7])")
    println("Stationarity: $(data[problem, 8])")
    println("δ: $(data[problem, 9])")
    println("Δ: $(data[problem, 10])")
    println("Index: $(data[problem, 11])")
    println("F. value: $(data[problem, 12])")
    println("Solution: $(data[problem, 13:(12+data[problem, 1])])")

end

function print_sol(
                    filename::String
                    )

    # Reads the execution data
    data = readdlm(filename)
    n = size(data)[1]

    # Print data
    for i = 1:n
        println("$(i) - Solution: $(data[i, 13:(12+data[i, 1])])")
    end

end

filename = "../data_files/mw_uncons_test_beta_1.dat"

print_sol(filename)

#read_data(4, filename)


text:
@assert begin
    res = PropBiology.simulation(
        [0.0, 1.0] ./ 4,
        [0.2, 0.1] ./ 4,
        9000.0 / 60.0;
        γ = [deg2rad(45.0), deg2rad(15.0)] .+ deg2rad(15.0),
        V∞ = 10.0
    )

    true
end

max_stars_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/PropBiology


text:
table = dbc_table(
    # using the same table as in the above example
    table_header * table_body,
    bordered=true,
    dark=true,
    hover=true,
    responsive=true,
    striped=true,
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# AUTO GENERATED FILE - DO NOT EDIT

export paper

"""
    paper(;kwargs...)
    paper(children::Any;kwargs...)
    paper(children_maker::Function;kwargs...)

A Paper component.
Renders white or dark background depending on color scheme. For more information, see: https://mantine.dev/core/paper/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Paper content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `padding` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Predefined padding value from theme.spacing or number for padding in px
- `radius` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Predefined border-radius value from theme.radius or number for border-radius in px
- `shadow` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Predefined box-shadow from theme.shadows (xs, sm, md, lg, xl) or any valid css box-shadow property
- `withBorder` (Bool; optional): Adds 1px border with theme.colors.gray[2] color in light color scheme and theme.colors.dark[6] in dark color scheme
"""
function paper(; kwargs...)
        available_props = Symbol[:children, :id, :className, :padding, :radius, :shadow, :withBorder]
        wild_props = Symbol[]
        return Component("paper", "Paper", "dash_mantine_components", available_props, wild_props; kwargs...)
end

paper(children::Any; kwargs...) = paper(;kwargs..., children = children)
paper(children_maker::Function; kwargs...) = paper(children_maker(); kwargs...)


text:
"""
    exclusivenearestneighbors

Exclusively pair points between two provided lists.

Finds the closest of all pairs, removes both, and iterates on the remainder until all pairs have been created.

NOTE: Assumes that your data is far away from the type's maximum value.
"""
function exclusivenearestneighbors(points1::Vector{T}, points2::Vector{T}) where T <: DataEntity
    points2larger = length(points2) > length(points1)

    (points1, points2) = points2larger ?
        (points2, copy(points1)) :
        (points1, copy(points2))

    coordinates1 = extractcoordinates(points1)
    coordinates1 = ndims(coordinates1) == 1 ?
        make2d(coordinates1) :
        coordinates1

    neighbors1 = T[]
    neighbors2 = T[]
    distances = eltype(coordinates1)[]

    neighbortree = BallTree(coordinates1)
    coordinates1used = falses(length(points1))

    while !isempty(points2)
        coordinates2 = extractcoordinates(points2)
        coordinates2used = falses(length(points2))
        inearestneighbors, nearestneighbor_distances = nn(neighbortree, coordinates2, i -> coordinates1used[i])
        while any(nearestneighbor_distances .< typemax(eltype(coordinates1)))
            mindist, mini = findmin(nearestneighbor_distances)
            nearestneighbor_distances[mini] = typemax(eltype(coordinates1))
            coordinates2used[mini] && continue
            inearestneighbor = inearestneighbors[mini]
            coordinates1used[inearestneighbor] && continue
            coordinates2used[mini] = coordinates1used[inearestneighbor] = true
            push!(neighbors1, points1[inearestneighbor])
            push!(neighbors2, points2[mini])
            push!(distances, mindist)
        end
        deleteat!(points2, coordinates2used)
    end

    if points2larger
        neighbors2, neighbors1, distances
    else
        neighbors1, neighbors2, distances
    end
end

"""
    montecarloaffinity()

Evaluate the probability of chance association across all molecules by randomizing their neighbors.
"""
function montecarloaffinity(molecules1::Vector{T}, molecules2::Vector{T}, ch1_neighbors::Vector{T},
                            ch2_neighbors::Vector{T}, distances, localradius, iterations) where T <: DataEntity
    coordinates1, coordinates2 = extractcoordinates.([molecules1, molecules2])

    percentileranks = ones(length(ch1_neighbors))

    neighborcoordinates1, neighborcoordinates2 = extractcoordinates.([ch1_neighbors, ch2_neighbors])

    centerxcoordinates = mean([neighborcoordinates1[1, :]'; neighborcoordinates2[1, :]'], dims = 1)
    centerycoordinates = mean([neighborcoordinates1[2, :]'; neighborcoordinates2[2, :]'], dims = 1)
    centerzcoordinates = mean([neighborcoordinates1[3, :]'; neighborcoordinates2[3, :]'], dims = 1)
    centercoordinates = [centerxcoordinates; centerycoordinates; centerzcoordinates]

    neighbor1tree, neighbor2tree = BallTree.([coordinates1, coordinates2])
    ineighbors1, ineighbors2 = inrange.([neighbor1tree, neighbor2tree], Ref(centercoordinates), localradius, true)
    nlocalmolecules1, nlocalmolecules2 = [length.(x) for x ∈ [ineighbors1, ineighbors2]]

    unique_nlocalmolecules1 = filter(x -> x > 0, unique(nlocalmolecules1))
    unique_nlocalmolecules2 = filter(x -> x > 0, unique(nlocalmolecules2))

    # precompute for the possible sets of distances, as there are typically many fewer combinations than their are potential pairs.
    localcounts = [(i,j) for i in unique_nlocalmolecules1 for j in unique_nlocalmolecules2]
    localmindistances = pmap(localmontecarlo_mindistance, localcounts, (localradius for i = eachindex(localcounts)), (iterations for i = eachindex(localcounts)), on_error = identity)
    localmindistancesdict = Dict{Tuple{Int, Int}, Vector{Float64}}(zip(localcounts, localmindistances))
    mindistances = [first(x) == 0 || last(x) == 0 ? repeat([0.0], iterations) : localmindistancesdict[x] for x ∈ zip(nlocalmolecules1, nlocalmolecules2)]
    percentileranks = [count(mindistances[i] .≤ distances[i]) / iterations for i ∈ eachindex(ch1_neighbors)]

    return percentileranks
end

"""
    localmontecarlo()

Evaluate the probabilty of chance association of a molecule pair by randomizing its neighbors.
"""
function localmontecarlo_mindistance(nlocalmolecules1_2, radius, iterations) where T <: DataEntity
    nlocalmolecules1 = first(nlocalmolecules1_2)
    nlocalmolecules2 = last(nlocalmolecules1_2)
    (nlocalmolecules1 == 0 || nlocalmolecules2 == 0) && return 1.0

    randomcoordinates1 = [randomcoordinates2d(nlocalmolecules1, radius) for i ∈ 1:iterations]
    randomcoordinates2 = [randomcoordinates2d(nlocalmolecules2, radius) for i ∈ 1:iterations]
    randomtrees = nlocalmolecules1 > nlocalmolecules2 ? KDTree.(randomcoordinates1) : KDTree.(randomcoordinates2)
    randomcoordinates = nlocalmolecules1 > nlocalmolecules2 ? randomcoordinates2 : randomcoordinates1
    mindistance = (nn(randomtrees[i], randomcoordinates[i]) |> last for i ∈ 1:iterations) .|> minimum
end


text:
"""
    normalize_name(file::AbstractString)

Given a filename, normalize it, stripping out extensions.  E.g. the file path
`"foo/libfoo.tar.gz"` would get mapped to `"libfoo"`.
"""
function normalize_name(file::AbstractString)
    file = basename(file)
    idx = findfirst(isequal('.'), file)
    if idx !== nothing
        file = file[1:prevind(file, idx)]
    end
    # Strip -123, which is a common thing for libraries on windows
    idx = findlast(isequal('-'), file)
    if idx !== nothing && all(isnumeric, file[nextind(file, idx):end])
        file = file[1:prevind(file, idx)]
    end
    return file
end

"""
    filter_object_files(files)

Given a list of files, filter out any that cannot be opened by `readmeta()`
from `ObjectFile`.
"""
function filter_object_files(files)
    return filter(files) do f
        try
            readmeta(f) do oh
                return true
            end
        catch e
            # If it was just a MagicMismatch, then return false for this file
            if isa(e, ObjectFile.MagicMismatch)
                return false
            end

            # If it was an EOFError (e.g. this was an empty file) then return false
            if isa(e, EOFError)
                return false
            end

            # If something else wrong then rethrow that error and pass it up
            rethrow(e)
            return false
        end
    end
end

"""
    match_files(state::WizardState, prefix::Prefix,
                platform::Platform, files::Vector; silent::Bool = false)

Inspects all binary files within a prefix, matching them with a given list of
`files`, complaining if there are any files that are not properly matched and
returning the set of normalized names that were not matched, or an empty set if
all names were properly matched.
"""
function match_files(state::WizardState, prefix::Prefix,
                     platform::Platform, files::Vector; silent::Bool = false)
    # Collect all executable/library files
    prefix_files = collapse_symlinks(collect_files(prefix))
    prefix_files = filter_object_files(prefix_files)
    # Check if we can load them as an object file
    prefix_files = filter(prefix_files) do f
        readmeta(f) do oh
            if !is_for_platform(oh, platform)
                if !silent
                    @warn("Skipping binary `$f` with incorrect platform")
                end
                return false
            end
            return true
        end
    end

    norm_prefix_files = Set(map(normalize_name, prefix_files))
    norm_files = Set(map(normalize_name, files))
    d = setdiff(norm_files, norm_prefix_files)
    if !isempty(d)
        if !silent
            @warn("Could not find correspondences for $(join(d, ' '))")
        end
    end
    return d
end

"""
    edit_script(state::WizardState, script::AbstractString)

For consistency (and security), use the sandbox for editing a script, launching
`vi` within an interactive session to edit a buildscript.
"""
function edit_script(state::WizardState, script::AbstractString)
    # Create a temporary directory to hold the script inside of
    temp_prefix() do prefix
        # Write the script out to a file
        path = joinpath(prefix, "script")
        open(path, "w") do f
            write(f, script)
        end

        # Launch a sandboxed vim editor
        ur = preferred_runner()(
            prefix.path,
            cwd = "/workspace/",
            platform = Linux(:x86_64))
        run_interactive(ur, `/usr/bin/vi /workspace/script`;
                        stdin=state.ins, stdout=state.outs, stderr=state.outs)

        # Once the user is finished, read the script back in
        script = String(read(path))
    end
    return script
end

"""
    yn_prompt(state::WizardState, question::AbstractString, default = :y)

Perform a `[Y/n]` or `[y/N]` question loop, using `default` to choose between
the prompt styles, and looping until a proper response (e.g. `"y"`, `"yes"`,
`"n"` or `"no"`) is received.
"""
function yn_prompt(state::WizardState, question::AbstractString, default = :y)
    @assert default in (:y, :n)
    ynstr = default == :y ? "[Y/n]" : "[y/N]"
    while true
        print(state.outs, question, " ", ynstr, ": ")
        answer = lowercase(strip(readline(state.ins)))
        if isempty(answer)
            return default
        elseif answer == "y" || answer == "yes"
            return :y
        elseif answer == "n" || answer == "no"
            return :n
        else
            println(state.outs, "Unrecognized answer. Answer `y` or `n`.")
        end
    end
end

function script_for_dep(dep, install_dir)
    # If dep is a pair, peel it to override install_dir
    if isa(dep, Pair)
        dep, install_dir = dep
    end

    # Since remote dependencies are most common, we default plain strings to
    # this in order to keep the build scripts small
    if isa(dep, String)
        dep = RemoteBuildDependency(dep, nothing)
    end

    if isa(dep, InlineBuildDependency)
        script = dep.script
    elseif isa(dep, RemoteBuildDependency)
        script = dep.script === nothing ? String(HTTP.get(dep.url).body) :
            dep.script
    elseif isa(dep, TarballDependency)
        script = """
        prefix = Prefix(ARGS[1])
        install($(repr(dep.url)), $(repr(dep.hash)); prefix=prefix, force=true)
        download_info = Dict(platform_key_abi() => ($(repr(dep.url)), $(repr(dep.hash))))
        """
    end
    return script, install_dir
end

"""
    setup_workspace(build_path::String, src_paths::Vector,
                    src_hashes::Vector, platform::Platform,
                    extra_env::Dict{String, String};
                    verbose::Bool = false, tee_stream::IO = stdout,
                    downloads_dir = nothing)

Sets up a workspace within `build_path`, creating the directory structure
needed by further steps, unpacking the source within `build_path`, and defining
the environment variables that will be defined within the sandbox environment.

This method returns the `Prefix` to install things into, and the runner
that can be used to launch commands within this workspace.
"""
function setup_workspace(build_path::AbstractString, src_paths::Vector,
                         src_hashes::Vector, dependencies::Vector,
                         platform::Platform,
                         extra_env::Dict{String, String} =
                             Dict{String, String}();
                         verbose::Bool = false,
                         tee_stream::IO = stdout,
                         downloads_dir = nothing)
    # Use a random nonce to make detection of paths in embedded binary easier
    nonce = randstring()
    mkdir(joinpath(build_path, nonce))

    # We now set up two directories, one as a source dir, one as a dest dir
    srcdir = joinpath(build_path, nonce, "srcdir")
    destdir = joinpath(build_path, nonce, "destdir")
    metadir = joinpath(build_path, nonce, "metadir")
    mkdir(srcdir); mkdir(destdir); mkdir(metadir)

    # Create a runner to work inside this workspace with the nonce built-in
    ur = preferred_runner()(
        joinpath(build_path, nonce),
        cwd = "/workspace/srcdir",
        platform = platform,
        extra_env = merge(extra_env,
            Dict(
                "prefix" => "/workspace/destdir",
                "WORKSPACE" => "/workspace",
                "PS1" => "sandbox:\${PWD//\$WORKSPACE/\\\\\$WORKSPACE}\\\$ "
             )
        ),
        verbose = verbose,
        workspaces = [
            metadir => "/meta",
        ],
    )

    # For each source path, unpack it
    for (src_path, src_hash) in zip(src_paths, src_hashes)
        if isdir(src_path)
            # Chop off the `.git` at the end of the src_path
            repo_dir = joinpath(srcdir, basename(src_path)[1:end-4])
            LibGit2.with(LibGit2.clone(src_path, repo_dir)) do repo
                LibGit2.checkout!(repo, src_hash)
            end
        else
            # Extract with host tools because it is _much_ faster on e.g. OSX.
            # If this becomes a compatibility problem, we'll just have to install
            # our own `tar` and `unzip` through BP as dependencies for BB.
            tar_extensions = [".tar", ".tar.gz", ".tgz", ".tar.bz", ".tar.bz2",
                              ".tar.xz", ".tar.Z", ".txz"]
            cd(joinpath(build_path, nonce, "srcdir")) do
                if any(endswith(src_path, ext) for ext in tar_extensions)
                    if verbose
                        @info "Extracting tarball $(basename(src_path))..."
                    end
                    tar_flags = verbose ? "xvof" : "xof"
                    run(`tar $(tar_flags) $(src_path)`)
                elseif endswith(src_path, ".zip")
                    if verbose
                        @info "Extracting zipball $(basename(src_path))..."
                    end
                    run(`unzip -q $(src_path)`)
               else
                   if verbose
                       @info "Copying in $(basename(src_path))..."
                   end
                   cp(src_path, basename(src_path))
                end
            end
        end
    end

    # If we haven't been given a downloads_dir, just dump it all into `srcdir`:
    if downloads_dir == nothing
        downloads_dir = srcdir
    end

    # For each dependency, install it into the prefix
    for dep in dependencies
        script, install_dir = script_for_dep(dep, destdir)
        m = Module(:__anon__)
        Core.eval(m, quote
            using BinaryProvider

            # Force the script to download for this platform.
            platform_key() = $platform
            platform_key_abi() = $platform

            # We don't want any deps files being written out
            function write_deps_file(args...; kwargs...) end

            # Override @__DIR__ to return the destination directory,
            # so that things get installed into there.  This is a protection
            # against older scripts that ignore `ARGS[1]`, which is set below.
            # Eventually we should be able to forgo this skullduggery.
            macro __DIR__(args...); return $destdir; end

            # Override install() to cache in our downloads directory, to
            # ignore platforms, and to be verbose if we want it to be.
            function install(url, hash; kwargs...)
                BinaryProvider.install(url, hash;
                    tarball_path=joinpath($downloads_dir, string(hash, "-", basename(url))),
                    ignore_platform=true,
                    verbose=$verbose,
                    kwargs...,
                )
            end

            ARGS = [$install_dir]
        end)
        include_string(m, script)
    end

    # Return the prefix and the runner
    return Prefix(destdir), ur
end

"""
    Change the script. This will invalidate all platforms to make sure we later
    verify that they still build with the new script.
"""
function change_script!(state, script)
    state.history = script
    empty!(state.validated_platforms)
end

function print_wizard_logo(outs)
    logo = raw"""

            o      `.
           o*o      \'-_                 00000000: 01111111 $.
             \\      \;"".     ,;.--*    00000001: 01000101 $E
              \\     ,\''--.--'/         00000003: 01001100 $L
              :=\--<' `""  _   |         00000003: 01000110 $F
              ||\\     `" / ''--         00000004: 00000010  .
              `/_\\,-|    |              00000005: 00000001  .
                  \\/     L
                   \\ ,'   \
                 _/ L'   `  \
                /  /    /   /          Julia Binzard
               /  /    |    \          JuliaPackaging/BinaryBuilder.jl
              "_''--_-''---__=;

    """

    blue    = "\033[34m"
    red     = "\033[31m"
    green   = "\033[32m"
    magenta = "\033[35m"
    normal  = "\033[0m\033[0m"

    logo = replace(logo, " o " => " $(green)o$(normal) ")
    logo = replace(logo, "o*o" => "$(red)o$(blue)*$(magenta)o$(normal)")

    logo = replace(logo, ".--*" => "$(red).$(green)-$(magenta)-$(blue)*$(normal)")

    logo = replace(logo, "\$." => " $(blue).$(normal)")
    logo = replace(logo, "\$E" => " $(red)E$(normal)")
    logo = replace(logo, "\$L" => " $(green)L$(normal)")
    logo = replace(logo, "\$F" => " $(magenta)F$(normal)")

    print(outs, logo)
end


text:
"""
po()

Returns information about the platform ontology.

### Examples

```
julia> po()
"10.15.6\n"

julia> po()
Apple clang version 11.0.3 . . .
```
"""
function po()
    # anIOBuffer = IOBuffer();
    # versioninfo(anIOBuffer, verbose=true);
    # println("System Version Info: ", String(take!(anIOBuffer)))
    # println("*\n")
    println("Clang: ", read(`clang --version`, String))
    println("*\n")
    println("Load paths: ", Base.load_path())
    println("*\n")
    #println("Jupyter kernels: ", run(`jupyter kernelspec list`))
    println("*\n")
    println("OS X: ", read(`sw_vers -productVersion`, String))
    println("*\n")
    println("Julia: ", read(`julia -v`, String))
    println("*\n")
    #println("Jupyter: ", read(`jupyter --version`, String))
    println("*\n")
    println("GitKraken: ", read(`/Applications/GitKraken.app/Contents/MacOS/GitKraken --version`, String))
    println("*\n")
    println("git: ", read(`git --version`, String))
    println("*\n")
    println("Atom: ", read(`apm -v`, String))
    println("*\n")
    println("Atom, Julia-client: ", read(`apm show julia-client`, String))
    println("*\n")
    println("Atom, Ink: ", read(`apm show ink`, String))
    println("*\n")
    println("Python: ", read(`python --version`, String))
    println("*\n")
    println("Clang: ", read(`clang --version`, String))
    println("*\n")
    println("VS Code: ", read(`code --version`, String))
    println("*\n")
    println("NOTE: for additional information use: versioninfo(verbose=true)")

end