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.popgrw_populationgrowth_ann[yr, r] = p.popgrw_populationgrowth[t, r]
            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 || use_logburke || use_logpopulation || use_logwherepossible
                    v.popgrw_populationgrowth_ann[yr, r] = p.popgrw_populationgrowth[t, r] * fraction_timestep + p.popgrw_populationgrowth[t - 1, r] * (1 - fraction_timestep)
                else
                    error("NO INTERPOLATION METHOD SELECTED! Specify linear or logarithmic interpolation.")
                end

            end
        end
    end
end

function calc_equityweighting(p, v, d, t, annual_year)
    # 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))

    if p.discfix_fixediscountrate != 0.
        v.df_utilitydiscountfactor_ann[yr] = (1 + p.discfix_fixediscountrate / 100)^(-(p.y_year_ann[yr] - p.y_year_0))
    end

    for r in d.region

        ## Gas Costs Accounting
        # Weighted costs (Page 23 of Hope 2009)
        v.wtct_percap_weightedcosts_ann[yr, r] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (p.cons_percap_consumption_ann[yr, r]^(1 - p.emuc_utilityconvexity) - (p.cons_percap_consumption_ann[yr, r] - p.tct_percap_totalcosts_total_ann[yr, r] < 0.01 * p.cons_percap_consumption_0[1] ? 0.01 * p.cons_percap_consumption_0[1] : p.cons_percap_consumption_ann[yr, r] - p.tct_percap_totalcosts_total_ann[yr, r])^(1 - p.emuc_utilityconvexity))

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

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

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

        # Discount rate calculations
        v.dr_discountrate_ann[yr, r] = p.ptp_timepreference + p.emuc_utilityconvexity * (p.grwnet_realizedgdpgrowth_ann[yr, r] - v.popgrw_populationgrowth_ann[yr, r])

        if p.discfix_fixediscountrate != 0.
            v.dr_discountrate_ann[yr, r] = p.discfix_fixediscountrate
        end

        v.yp_yearsperiod_ann[yr] = 1 # every timestep is 1 year long.

        if is_first(t)
            v.dfc_consumptiondiscountrate_ann[yr, r] = (1 + v.dr_discountrate_ann[yr, r] / 100)^(-v.yp_yearsperiod_ann[yr])
        else
            v.dfc_consumptiondiscountrate_ann[yr, r] = v.dfc_consumptiondiscountrate_ann[yr - 1, r] * (1 + v.dr_discountrate_ann[yr, r] / 100)^(-v.yp_yearsperiod_ann[yr])
        end

        # Discounted costs
        if p.equity_proportion == 0
            v.pcdt_partiallyweighted_discounted_ann[yr, r] = v.pct_partiallyweighted_ann[yr, r] * v.dfc_consumptiondiscountrate_ann[yr, r]
            v.wacdt_partiallyweighted_discounted_ann[yr, r] = p.act_adaptationcosts_total_ann[yr, r] * v.dfc_consumptiondiscountrate_ann[yr, r]
        else
            v.pcdt_partiallyweighted_discounted_ann[yr, r] = v.pct_partiallyweighted_ann[yr, r] * v.df_utilitydiscountfactor_ann[yr]
            v.wacdt_partiallyweighted_discounted_ann[yr, r] = v.wact_partiallyweighted_ann[yr, r] * v.df_utilitydiscountfactor_ann[yr]
        end

        v.pcdat_partiallyweighted_discountedaggregated_ann[yr, r] = v.pcdt_partiallyweighted_discounted_ann[yr, r]

        # calculate the total damages due to impacts
        v.currentdam_currentdamages_ann[yr, r] = (p.cons_percap_aftercosts_ann[yr, r] - p.rcons_percap_dis_ann[yr, r]) * p.pop_population_ann[yr, r]
        v.currentdampc_percapitacurrentdamages_ann[yr, r] = v.currentdam_currentdamages_ann[yr, r] / p.pop_population_ann[yr, r]
        v.damshare_currentdamagesshare_ann[yr, r] = 100 * v.currentdampc_percapitacurrentdamages_ann[yr, r] / p.cons_percap_aftercosts_ann[yr, r]

        # if damages including GDP losses exceed current consumption levels, calculate the share for equity weighting based on the convergence system
        v.eqwshare_shareofweighteddamages_noconvergence_ann[yr,r] = (v.currentdam_currentdamages_ann[yr, r] + p.lgdp_gdploss_ann[yr,r]) / (p.cons_percap_aftercosts_ann[yr, r] * p.pop_population_ann[yr, r])

        if p.use_convergence == 1.
            if v.eqwshare_shareofweighteddamages_noconvergence_ann[yr,r] > p.eqwboundn_maxshareofweighteddamages_neighbourhood
                v.eqwshare_shareofweighteddamages_ann[yr,r] = p.eqwboundn_maxshareofweighteddamages_neighbourhood - 0.5 * v.eqwaux1_weighteddamages_auxiliary1 +
                                v.eqwaux1_weighteddamages_auxiliary1 * exp(v.eqwaux2_weighteddamages_auxiliary2 *
                                                                        (v.eqwshare_shareofweighteddamages_noconvergence_ann[yr,r] - p.eqwboundn_maxshareofweighteddamages_neighbourhood)) /
                                                                    (1 + exp(v.eqwaux2_weighteddamages_auxiliary2 *
                                                                            (v.eqwshare_shareofweighteddamages_noconvergence_ann[yr,r] - p.eqwboundn_maxshareofweighteddamages_neighbourhood)))
                # for very large excess damages, exp(...) becomes infite and the convergence system returns NaN. In this case, simply set eqwshare to the upper bound
                if isnan(v.eqwshare_shareofweighteddamages_ann[yr,r])
                    v.eqwshare_shareofweighteddamages_ann[yr,r] = p.eqwbound_maxshareofweighteddamages
                end

                v.excdam_excessdamages_ann[yr, r] =  max(0, v.currentdam_currentdamages_ann[yr, r] + p.lgdp_gdploss_ann[yr,r] - v.eqwshare_shareofweighteddamages_ann[yr,r] * p.cons_percap_aftercosts_ann[yr, r] * p.pop_population_ann[yr, r])
            else
                v.excdam_excessdamages_ann[yr, r] =  0
            end
        else
            # use hard boundaries if use_convergence is not set to 1
            v.eqwshare_shareofweighteddamages_ann[yr,r] = min(v.eqwshare_shareofweighteddamages_noconvergence_ann[yr,r], p.eqwbound_maxshareofweighteddamages)
            v.excdam_excessdamages_ann[yr, r] =  max(0, v.currentdam_currentdamages_ann[yr, r] + p.lgdp_gdploss_ann[yr,r] - v.eqwshare_shareofweighteddamages_ann[yr,r] * p.cons_percap_aftercosts_ann[yr, r] * p.pop_population_ann[yr, r])
        end

        ## Equity weighted impacts (end of page 28, Hope 2009)
        if p.lossinc_includegdplosses == 0. && p.equity_proportion == 0.
            v.wit_partiallyweighted_ann[yr, r] = (p.cons_percap_aftercosts_ann[yr, r]  - p.rcons_percap_dis_ann[yr, r]) * p.pop_population_ann[yr, r]
            v.widt_partiallyweighted_discounted_ann[yr, r] = v.wit_partiallyweighted_ann[yr, r] * v.dfc_consumptiondiscountrate_ann[yr, r]
        elseif p.lossinc_includegdplosses == 0. && p.equity_proportion == 1.
            v.wit_partiallyweighted_ann[yr, r] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (p.cons_percap_aftercosts_ann[yr, r]^(1 - p.emuc_utilityconvexity) - p.rcons_percap_dis_ann[yr, r]^(1 - p.emuc_utilityconvexity)) * p.pop_population_ann[yr, r]
            v.widt_partiallyweighted_discounted_ann[yr, r] = v.wit_partiallyweighted_ann[yr, r] * v.df_utilitydiscountfactor_ann[yr]
        elseif p.lossinc_includegdplosses == 1. && p.equity_proportion == 0.
            v.wit_partiallyweighted_ann[yr, r] = (p.cons_percap_aftercosts_ann[yr, r]  - p.rcons_percap_dis_ann[yr, r]) * p.pop_population_ann[yr, r] + p.lgdp_gdploss_ann[yr, r]
            v.widt_partiallyweighted_discounted_ann[yr, r] = v.wit_partiallyweighted_ann[yr, r] * v.dfc_consumptiondiscountrate_ann[yr, r]
        elseif p.lossinc_includegdplosses == 1. && p.equity_proportion == 1.
            if v.excdam_excessdamages_ann[yr, r] == 0
                v.wit_partiallyweighted_ann[yr, r] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (p.cons_percap_aftercosts_ann[yr, r]^(1 - p.emuc_utilityconvexity) - p.rcons_percap_dis_ann[yr, r]^(1 - p.emuc_utilityconvexity)) * p.pop_population_ann[yr, r]
                v.widt_partiallyweighted_discounted_ann[yr, r] = v.wit_partiallyweighted_ann[yr, r] * v.df_utilitydiscountfactor_ann[yr]
            else
                v.wit_partiallyweighted_ann[yr, r] = ((p.cons_percap_consumption_0[1]^p.emuc_utilityconvexity) / (1 - p.emuc_utilityconvexity)) * (p.cons_percap_aftercosts_ann[yr, r]^(1 - p.emuc_utilityconvexity) - ((1 - v.eqwshare_shareofweighteddamages_ann[yr,r]) * p.cons_percap_aftercosts_ann[yr, r])^(1 - p.emuc_utilityconvexity)) * p.pop_population_ann[yr, r]
                v.widt_partiallyweighted_discounted_ann[yr, r] = v.wit_partiallyweighted_ann[yr, r] * v.df_utilitydiscountfactor_ann[yr]  + v.excdam_excessdamages_ann[yr, r] *  v.dfc_consumptiondiscountrate_ann[yr, r]
            end
        end

        v.excdampv_excessdamagespresvalue_ann[yr, r] = v.excdam_excessdamages_ann[yr, r] *  v.dfc_consumptiondiscountrate_ann[yr, r]

        v.addt_equityweightedimpact_discountedaggregated_ann[yr, r] = v.widt_partiallyweighted_discounted_ann[yr, r]
        v.aact_equityweightedadaptation_discountedaggregated_ann[yr, r] = v.wacdt_partiallyweighted_discounted_ann[yr, r]

    end

    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, :])

    v.td_totaldiscountedimpacts_ann = min(v.addt_gt_equityweightedimpact_discountedglobal_ann, p.civvalue_civilizationvalue)
    v.td_totaldiscountedimpacts_ann_yr[yr] = v.td_totaldiscountedimpacts_ann # for saving as a timeseries for possible post-analysis

    # Total effect of climate change
    v.te_totaleffect_ann = min(v.td_totaldiscountedimpacts_ann + v.tpc_totalaggregatedcosts_ann + v.tac_totaladaptationcosts_ann, p.civvalue_civilizationvalue)
    v.te_totaleffect_ann_yr[yr] = v.te_totaleffect_ann # for saving as a timeseries for possible post-analysis
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 = Parameter(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 = Parameter(index=[year, region], unit="\$/person")
    act_adaptationcosts_total = Parameter(index=[time, region], unit="\$million")
    act_adaptationcosts_total_ann = Parameter(index=[year, region], unit="\$million")
    act_percap_adaptationcosts = Parameter(index=[time, region], unit="\$/person")
    act_percap_adaptationcosts_ann = Parameter(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 = Parameter(index=[year, region], unit="\$/person")
    cons_percap_aftercosts = Parameter(index=[time, region], unit="\$/person")
    cons_percap_aftercosts_ann = Parameter(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>
    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

    ###############################################
    # Growth Effects - additional variables and parameters
    ###############################################
    # additional paramters and variables for growth effects and boundaries
    grwnet_realizedgdpgrowth = Parameter(index=[time, region], unit="%/year")
    grwnet_realizedgdpgrowth_ann = Parameter(index=[year, region], unit="%/year")                  # interpolated in GDP
    lgdp_gdploss =  Parameter(index=[time, region], unit="\$M")
    lgdp_gdploss_ann =  Parameter(index=[year, region], unit="\$M")                                  # interpolated in GDP
    lossinc_includegdplosses = Parameter(unit="none", default=1.)
    excdam_excessdamages = Variable(index=[time,region], unit="\$million")
    excdam_excessdamages_ann = Variable(index=[year,region], unit="\$million")
    excdampv_excessdamagespresvalue = Variable(index=[time,region], unit="\$million")
    excdampv_excessdamagespresvalue_ann = Variable(index=[year,region], unit="\$million")

    # convergence system for equity weighting threshold
    use_convergence = Parameter(unit="none", default=1.)
    eqwshare_shareofweighteddamages = Variable(index=[time,region], unit="none")
    eqwshare_shareofweighteddamages_ann = Variable(index=[year,region], unit="none")
    eqwshare_shareofweighteddamages_noconvergence = Variable(index=[time,region], unit="none")
    eqwshare_shareofweighteddamages_noconvergence_ann = Variable(index=[year,region], unit="none")
    eqwbound_maxshareofweighteddamages = Parameter(unit="none", default=0.99)
    eqwboundn_maxshareofweighteddamages_neighbourhood = Parameter(unit="none", default=0.9)
    eqwaux1_weighteddamages_auxiliary1 = Variable(unit="none")
    eqwaux2_weighteddamages_auxiliary2 = Variable(unit="none")

    currentdam_currentdamages = Variable(index=[time,region], unit="\$million")
    currentdam_currentdamages_ann = Variable(index=[year,region], unit="\$million")
    damshare_currentdamagesshare = Variable(index=[time,region], unit="%GDP")
    damshare_currentdamagesshare_ann = Variable(index=[year,region], unit="%GDP")
    currentdampc_percapitacurrentdamages = Variable(index=[time,region], unit="\$/person")
    currentdampc_percapitacurrentdamages_ann = Variable(index=[year,region], unit="\$/person")

    discfix_fixediscountrate = Parameter(unit="none", default=0.) # override the discount rates with something exogenous
    ###############################################

    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

            v.eqwaux1_weighteddamages_auxiliary1 = 2 * (p.eqwbound_maxshareofweighteddamages - p.eqwboundn_maxshareofweighteddamages_neighbourhood)
            v.eqwaux2_weighteddamages_auxiliary2 = 4 / v.eqwaux1_weighteddamages_auxiliary1

        end

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

        if p.discfix_fixediscountrate != 0.
            v.df_utilitydiscountfactor[tt] = (1 + p.discfix_fixediscountrate / 100)^(-(p.y_year[tt] - p.y_year_0))
        end

        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.grwnet_realizedgdpgrowth[tt, rr] - p.popgrw_populationgrowth[tt, rr])

            if p.discfix_fixediscountrate != 0.
                v.dr_discountrate[tt, rr] = p.discfix_fixediscountrate
            end

            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]
            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]
            end

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

            # calculate the total damages due to impacts
            v.currentdam_currentdamages[tt, rr] = (p.cons_percap_aftercosts[tt, rr] - p.rcons_percap_dis[tt, rr]) * p.pop_population[tt, rr]
            v.currentdampc_percapitacurrentdamages[tt, rr] = v.currentdam_currentdamages[tt, rr] / p.pop_population[tt, rr]
            v.damshare_currentdamagesshare[tt, rr] = 100 * v.currentdampc_percapitacurrentdamages[tt, rr] / p.cons_percap_aftercosts[tt, rr]

            # if damages including GDP losses exceed current consumption levels, calculate the share for equity weighting based on the convergence system
            v.eqwshare_shareofweighteddamages_noconvergence[tt,rr] = (v.currentdam_currentdamages[tt, rr] + p.lgdp_gdploss[tt,rr]) / (p.cons_percap_aftercosts[tt, rr] * p.pop_population[tt, rr])

            if p.use_convergence == 1.
                if v.eqwshare_shareofweighteddamages_noconvergence[tt,rr] > p.eqwboundn_maxshareofweighteddamages_neighbourhood
                    v.eqwshare_shareofweighteddamages[tt,rr] = p.eqwboundn_maxshareofweighteddamages_neighbourhood - 0.5 * v.eqwaux1_weighteddamages_auxiliary1 +
                                    v.eqwaux1_weighteddamages_auxiliary1 * exp(v.eqwaux2_weighteddamages_auxiliary2 *
                                                                            (v.eqwshare_shareofweighteddamages_noconvergence[tt,rr] - p.eqwboundn_maxshareofweighteddamages_neighbourhood)) /
                                                                        (1 + exp(v.eqwaux2_weighteddamages_auxiliary2 *
                                                                                (v.eqwshare_shareofweighteddamages_noconvergence[tt,rr] - p.eqwboundn_maxshareofweighteddamages_neighbourhood)))
                    # for very large excess damages, exp(...) becomes infite and the convergence system returns NaN. In this case, simply set eqwshare to the upper bound
                    if isnan(v.eqwshare_shareofweighteddamages[tt,rr])
                        v.eqwshare_shareofweighteddamages[tt,rr] = p.eqwbound_maxshareofweighteddamages
                    end

                    v.excdam_excessdamages[tt, rr] =  max(0, v.currentdam_currentdamages[tt, rr] + p.lgdp_gdploss[tt,rr] - v.eqwshare_shareofweighteddamages[tt,rr] * p.cons_percap_aftercosts[tt, rr] * p.pop_population[tt, rr])
                else
                    v.excdam_excessdamages[tt, rr] =  0
                end
            else
                # use hard boundaries if use_convergence is not set to 1
                v.eqwshare_shareofweighteddamages[tt,rr] = min(v.eqwshare_shareofweighteddamages_noconvergence[tt,rr], p.eqwbound_maxshareofweighteddamages)
                v.excdam_excessdamages[tt, rr] =  max(0, v.currentdam_currentdamages[tt, rr] + p.lgdp_gdploss[tt,rr] - v.eqwshare_shareofweighteddamages[tt,rr] * p.cons_percap_aftercosts[tt, rr] * p.pop_population[tt, rr])
            end

            ## Equity weighted impacts (end of page 28, Hope 2009)
            if p.lossinc_includegdplosses == 0. && p.equity_proportion == 0.
                v.wit_partiallyweighted[tt, rr] = (p.cons_percap_aftercosts[tt, rr]  - p.rcons_percap_dis[tt, rr]) * p.pop_population[tt, rr]
                v.widt_partiallyweighted_discounted[tt, rr] = v.wit_partiallyweighted[tt, rr] * v.dfc_consumptiondiscountrate[tt, rr]
            elseif p.lossinc_includegdplosses == 0. && p.equity_proportion == 1.
                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]
            elseif p.lossinc_includegdplosses == 1. && p.equity_proportion == 0.
                v.wit_partiallyweighted[tt, rr] = (p.cons_percap_aftercosts[tt, rr]  - p.rcons_percap_dis[tt, rr]) * p.pop_population[tt, rr] + p.lgdp_gdploss[tt, rr]
                v.widt_partiallyweighted_discounted[tt, rr] = v.wit_partiallyweighted[tt, rr] * v.dfc_consumptiondiscountrate[tt, rr]
            elseif p.lossinc_includegdplosses == 1. && p.equity_proportion == 1.
                if v.excdam_excessdamages[tt, rr] == 0
                    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]
                else
                    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) - ((1 - v.eqwshare_shareofweighteddamages[tt,rr]) * p.cons_percap_aftercosts[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]  + v.excdam_excessdamages[tt, rr] *  v.dfc_consumptiondiscountrate[tt, rr]
                end
            end

            v.excdampv_excessdamagespresvalue[tt, rr] = v.excdam_excessdamages[tt, rr] *  v.dfc_consumptiondiscountrate[tt, rr]

            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]

        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, :])

        v.td_totaldiscountedimpacts = min(v.addt_gt_equityweightedimpact_discountedglobal, p.civvalue_civilizationvalue)

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

        # 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)
            end
        else
            for annual_year = (gettime(tt - 1) + 1):(gettime(tt))
                calc_equityweighting(p, v, d, tt, annual_year)
            end
        end

    end
end


text:
module DreamStars

using ..Ahorn, Maple

placements = Maple.DreamStars

function Ahorn.canFgBg(effect::Maple.DreamStars)
    return true, true
end

end

ext: jl


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

export Maybe, MaybeVector

"""
Haskell-style [option type](https://en.wikipedia.org/wiki/Option_type) primarily used to allow optional inputs/outputs
for functions.

For example, in `src/Examples/docs_examples.jl`, we have `draw_cooketriplet(filename::Maybe{AbstractString} = nothing)`.
If the function is called with a `filename` argument, then it is passed as a string, and `Vis.save(filename)` saves an
image of the drawn example. If there is no argument, then the save call is dispatched to the empty `function
Vis.save(::Nothing) end`. This allows the user to run the example without saving a file.
"""
const Maybe{T} = Union{T, Nothing}

"""
Extension of `Maybe{T}` which allows things like `[nothing, 1, 1.0, nothing] isa MaybeVector{Real}`
"""
const MaybeVector{T} = Vector{<:Maybe{T}}


text:
mutable struct LatinGenerator <: AbstractModelGenerator
    N::Int
    p::Float64
    A::Union{Nothing,Matrix{Int}}
    function LatinGenerator(N,p)
        new(N,p,nothing)
    end
end

"""
    fill_with_generator!(cpmodel::CPModel, gen::GraphColoringGenerator)
Fill a CPModel with the variables and constraints generated. We fill it directly instead of
creating temporary files for efficiency purpose.

This generator create graps for the NQueens problem.

"""

function proper_move(A::Matrix{Int},x::Int,y::Int,z::Int,v::Int)
    z2 = A[x,y]
    if z2==v
        A[x,y]=z
    else
        y2 = rand(findall(k->k==v, A[x,:]))
        x2 = rand(findall(k->k==v, A[:,y]))
        A[x,y] = z
        A[x,y2]= z2
        A[x2,y]= z2
        proper_move(A,x2,y2,v,z2)
    end
end

function fill_with_generator!(cpmodel::CPModel, gen::LatinGenerator; seed=nothing)
    N = gen.N
    p = gen.p
    cpmodel.limit.numberOfSolutions = 1
    if isnothing(gen.A)
        gen.A = Matrix{Int}(undef,N,N)
        for i in 1:N gen.A[i,:]= [(i+j-2)%N + 1 for j in 1:N] end
    end
    A = gen.A
    if !isnothing(seed)
        Random.seed!(seed)
    end
    for i in 1:N^3
        x,y,z = rand(1:N,3)
        proper_move(A,x,y,z,z)
    end
    B = copy(A)
    n = floor(Int,p*N^2)
    indicies = shuffle(1:N^2)[1:n]
    for x in indicies
        i = div(x-1,N) + 1
        j = (x-1)%N + 1
        B[i,j] = 0
    end

    puzzle = Matrix{SeaPearl.AbstractIntVar}(undef, N,N)
    for i = 1:N
        for j in 1:N
            puzzle[i,j] = SeaPearl.IntVar(1, N, "puzzle_"*string(i)*","*string(j), cpmodel.trailer)
            SeaPearl.addVariable!(cpmodel, puzzle[i,j]; branchable=true)
            if B[i,j]>0 SeaPearl.addConstraint!(cpmodel,SeaPearl.EqualConstant(puzzle[i,j], B[i,j], cpmodel.trailer)) end
        end
    end
    for i in 1:N
        SeaPearl.addConstraint!(cpmodel, SeaPearl.AllDifferent(puzzle[i,:], cpmodel.trailer))
        SeaPearl.addConstraint!(cpmodel, SeaPearl.AllDifferent(puzzle[:,i], cpmodel.trailer))
    end
    return nothing
end


text:
#=

Code for Sparsification.

Right now, just implements Spielman-Srivastava

=#

"""
    as = sparsify(a; ep=0.5)

Apply Spielman-Srivastava sparsification: sampling by effective resistances.
`ep` should be less than 1.
"""
function sparsify(a; ep=0.3, matrixConcConst=4.0, JLfac=4.0)

  if ep > 1
    @warn "Calling sparsify with ep > 1 can produce a disconnected graph."
  end

  f = approxchol_lap(a,tol=1e-2);

  n = size(a,1)
  k = round(Int, JLfac*log(n)) # number of dims for JL

  U = wtedEdgeVertexMat(a)
  m = size(U,1)
  R = randn(Float64, m,k)
  UR = U'*R;

  V = zeros(n,k)
  for i in 1:k
    V[:,i] = f(UR[:,i])
  end

  (ai,aj,av) = findnz(triu(a))
  prs = zeros(size(av))
  for h in 1:length(av)
      i = ai[h]
      j = aj[h]
      prs[h] = min(1,av[h]* (norm(V[i,:]-V[j,:])^2/k) * matrixConcConst *log(n)/ep^2)
  end

  ind = rand(Float64,size(prs)) .< prs

  as = sparse(ai[ind],aj[ind],av[ind]./prs[ind],n,n)
  as = as + as'

  return as

end


text:
using Base.Test
using MPI

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

dst = mod(rank+1, size)
src = mod(rank-1, size)

N = 32

send_mesg = Array(Float64, N)
recv_mesg = Array(Float64, N)
recv_mesg_expected = Array(Float64, N)
fill!(send_mesg, Float64(rank))
fill!(recv_mesg_expected, Float64(src))

rreq = MPI.Irecv!(recv_mesg, src,  src+32, comm)
sreq = MPI.Isend(send_mesg, dst, rank+32, comm)

(ind,stats) = MPI.Waitsome!([sreq, rreq])
(done, ind, stats) = MPI.Testany!([sreq, rreq])
arr = [sreq,rreq]
if done
    (onedone,stats) = MPI.Test!(arr[ind])
    @test onedone
end

MPI.Finalize()


text:
const USE_GPU = true
using BenchmarkTools
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2);
else
    @init_parallel_stencil(Threads, Float64, 2);
end

@parallel_indices (ix,iy) function diffusion2D_step!(T2, T, Ci, lam, dt, _dx, _dy)
    tx  = @threadIdx().x + 1
    ty  = @threadIdx().y + 1
    T_l = @sharedMem(eltype(T), (@blockDim().x+2, @blockDim().y+2))
    T_l[tx,ty] = T[ix,iy]
    if (ix>1 && ix<size(T2,1) && iy>1 && iy<size(T2,2))
        if (@threadIdx().x == 1)             T_l[tx-1,ty] = T[ix-1,iy] end
        if (@threadIdx().x == @blockDim().x) T_l[tx+1,ty] = T[ix+1,iy] end
        if (@threadIdx().y == 1)             T_l[tx,ty-1] = T[ix,iy-1] end
        if (@threadIdx().y == @blockDim().y) T_l[tx,ty+1] = T[ix,iy+1] end
        @sync_threads()
        T2[ix,iy] = T_l[tx,ty] + dt*(Ci[ix,iy]*(
                    - ((-lam*(T_l[tx+1,ty] - T_l[tx,ty])*_dx) - (-lam*(T_l[tx,ty] - T_l[tx-1,ty])*_dx))*_dx
                    - ((-lam*(T_l[tx,ty+1] - T_l[tx,ty])*_dy) - (-lam*(T_l[tx,ty] - T_l[tx,ty-1])*_dy))*_dy
                    ));
    end
    return
end

function diffusion2D()
# Physics
lam      = 1.0;                                          # Thermal conductivity
c0       = 2.0;                                          # Heat capacity
lx, ly   = 1.0, 1.0;                                     # Length of computational domain in dimension x and y

# Numerics
nx, ny   = 512*32, 512*32;                               # Number of gridpoints in dimensions x and y
nt       = 100;                                          # Number of time steps
dx       = lx/(nx-1);                                    # Space step in x-dimension
dy       = ly/(ny-1);                                    # Space step in y-dimension
_dx, _dy = 1.0/dx, 1.0/dy;

# Array initializations
T   = @zeros(nx, ny);
T2  = @zeros(nx, ny);
Ci  = @zeros(nx, ny);

# Initial conditions
Ci .= 1/c0;                                              # 1/Heat capacity
T  .= 1.7;
T2 .= T;                                                 # Assign also T2 to get correct boundary conditions.
#
# GPU launch parameters
threads = (32, 8)
blocks  = (nx, ny) .÷ threads
#
# Time loop
dt   = min(dx^2,dy^2)/lam/maximum(Ci)/4.1;               # Time step for 2D Heat diffusion
for it = 1:nt
    if (it == 11) GC.gc(); global t_tic=time(); end      # Start measuring time.
    @parallel blocks threads shmem=prod(threads.+2)*sizeof(Float64) diffusion2D_step!(T2, T, Ci, lam, dt, _dx, _dy);
    T, T2 = T2, T;
end
time_s = time() - t_tic

# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*sizeof(eltype(T));           # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it  = time_s/(nt-10);                                  # Execution time per iteration [s]
T_eff = A_eff/t_it;                                      # Effective memory throughput [GB/s]
println("time_s=$time_s t_it=$t_it T_eff=$T_eff");

# Performance
A_eff = (2*1+1)*1/1e9*nx*ny*sizeof(eltype(T));           # Effective main memory access per iteration [GB] (Lower bound of required memory access: T has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
t_it = @belapsed begin @parallel $blocks $threads shmem=prod($threads.+2)*sizeof(Float64) diffusion2D_step!($T2, $T, $Ci, $lam, $dt, $_dx, $_dy); end
println("Benchmarktools (min): t_it=$t_it T_eff=$(A_eff/t_it)");

end

diffusion2D()


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -1.263105417837706 <= q <= 2.4766737340949954)

                     add_NL_constraint(m, :(gelu(0.7635571507341847 + -0.35250618093605723*gelu(0.5051094873509578 + -0.17934163301535788*gelu(-0.331542222824047 + 0.7450661875335634*gelu(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*gelu(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9610191438924418*gelu(0.6781688636512797 + -0.22580923517577522*gelu(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*gelu(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4])))) + 0.1709251594821839*gelu(0.058205314222445015 + 0.7954955842799771*gelu(-0.331542222824047 + 0.7450661875335634*gelu(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*gelu(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9762293687340509*gelu(0.6781688636512797 + -0.22580923517577522*gelu(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*gelu(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))) + gelu(-0.29663194049169395 + 0.013502946151282558*gelu(0.5051094873509578 + -0.17934163301535788*gelu(-0.331542222824047 + 0.7450661875335634*gelu(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*gelu(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9610191438924418*gelu(0.6781688636512797 + -0.22580923517577522*gelu(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*gelu(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4])))) + -0.2459643040595232*gelu(0.058205314222445015 + 0.7954955842799771*gelu(-0.331542222824047 + 0.7450661875335634*gelu(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + -0.10668129539960924*gelu(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))) + -0.9762293687340509*gelu(0.6781688636512797 + -0.22580923517577522*gelu(0.04378671584682481 + 0.6928572029860351*$(x[1]) + 0.8337172801246142*$(x[2]) + 0.1196869024103484*$(x[3]) + -0.7815194586429408*$(x[4])) + 0.2965219706242044*gelu(0.8574467662893772 + -0.06659837723966744*$(x[1]) + -0.15959391344635332*$(x[2]) + -0.6414263171723844*$(x[3]) + 0.9253455779716249*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
export hook!, addmixin!

type Mixin
    hooks::Dict{Symbol, Vector{Function}}

    Mixin() = new(Dict(i => Function[] for i in HOOKS))
end

hook!(m::Mixin, t::Symbol, f::Function) = hook!(m, t, Function[f])
hook!(m::Mixin, t::Symbol, f::Vector{Function}) = t in HOOKS ? for i in f push!(m.hooks[t], i) end : error("No hook called $t")

addmixin!(r::Resource, m::Mixin) = addmixin!(r, Mixin[m])
addmixin!(r::Resource, ms::Vector{Mixin}) = for m in ms
    for (k, v) in m.hooks
        for i in v push!(r.hooks[k], i) end
    end
end


text:
#
println("Fd) Symbolic simplification of Racah expressions including spherical harmonics and Wigner rotation matrices.")
#
using SymEngine
j1  = Basic(:j1);    j2  = Basic(:j2);    j3  = Basic(:j3);    j4  = Basic(:j4);    j5  = Basic(:j5);   j  = Basic(:j)
l1  = Basic(:l1);    l2  = Basic(:l2);    l3  = Basic(:l3);    l4  = Basic(:l4);    l5  = Basic(:l5)
m1  = Basic(:m1);    m2  = Basic(:m2);    m3  = Basic(:m3);    m4  = Basic(:m4);    m5  = Basic(:m5)
p1  = Basic(:p1);    p2  = Basic(:p2);    p3  = Basic(:p3);    p4  = Basic(:p4);    p5  = Basic(:p5);   p  = Basic(:p)
q1  = Basic(:q1);    q2  = Basic(:q2);    q3  = Basic(:q3);    q4  = Basic(:q4);    q5  = Basic(:q5);   q  = Basic(:q)

alpha  = Basic(:alpha);    beta  = Basic(:beta);    gamma  = Basic(:gamma)
alpha1 = Basic(:alpha1);   beta1 = Basic(:beta1);   gamma1 = Basic(:gamma1)
theta  = Basic(:theta);    phi   = Basic(:phi)
theta1 = Basic(:theta1);   phi1  = Basic(:phi1)
theta2 = Basic(:theta2);   phi2  = Basic(:phi2)
J   = Basic(:J)

if  true
    #
    w1 = RacahAlgebra.DFunction(j1, p1, q1, alpha, beta, gamma)
    w2 = RacahAlgebra.Djpq(j,p,q,beta)
    wa = RacahAlgebra.Ylm(l1,m1,theta1,phi1)
    RacahAlgebra.symmetricForms(w2)
    #
elseif false
    #
end


text:
module PoreMatMod_Test

using Test, PoreMatMod

@testset "correct_missing_Hs" begin
    include("../examples/correct_missing_Hs.jl")
    @test true
end

@testset "disorder_and_guests" begin
    include("../examples/disorder_and_guests.jl")
    @test true
end

@testset "make_hypothetical_MOF" begin
    include("../examples/make_hypothetical_MOF.jl")
    @test true
end

@testset "missing_linker_defect" begin
    include("../examples/missing_linker_defect.jl")
    @test true
end

@testset "replacement_modes" begin
    include("../examples/replacement_modes.jl")
    @test true
end

end


text:
module SQL

using Compat

import ..MySQL: MySQLHandle, MySQLInternalError, MySQLRowIterator, MYSQL_TUPLES, MYSQL_OPT_RECONNECT
import ..MySQL: mysql_connect, mysql_disconnect, mysql_stmt_prepare, mysql_metadata, mysql_execute

import Octo
import .Octo: QueryBuilders
import .QueryBuilders: DefaultArrangement, arrange

include("sql/adapter.jl")
include("sql/connect.jl")
include("sql/prepare.jl")
include("sql/iterate.jl")
include("sql/execute.jl")
include("sql/arrange.jl")
include("sql/table_from.jl")

end


text:
export @spinlock, @lock, LazyInitialized

const var"@lock" = Base.var"@lock"

# a safe way to acquire locks from finalizers, where we can't wait (which switches tasks)
macro spinlock(l, ex)
  quote
    temp = $(esc(l))
    while !trylock(temp)
      ccall(:jl_cpu_pause, Cvoid, ())
      # Temporary solution before we have gc transition support in codegen.
      ccall(:jl_gc_safepoint, Cvoid, ())
      # we can't yield here
    end
    try
      $(esc(ex))
    finally
      unlock(temp)
    end
  end
end

"""
    LazyInitialized{T}()

A thread-safe, lazily-initialized wrapper for a value of type `T`. Initialize and fetch the
value by calling `get!`. The constructor is ensured to only be called once.

This type is intended for lazy initialization of e.g. global structures, without using
`__init__`. It is similar to protecting accesses using a lock, but is much cheaper.

"""
struct LazyInitialized{T}
    # 0: uninitialized
    # 1: initializing
    # 2: initialized
    guard::Threads.Atomic{Int}
    value::Base.RefValue{T}

    LazyInitialized{T}() where {T} =
        new(Threads.Atomic{Int}(0), Ref{T}())
end

function Base.get!(constructor, x::LazyInitialized; hook=nothing)
    while x.guard[] != 2
        initialize!(x, constructor, hook)
    end
    assume(isassigned(x.value)) # to get rid of the check
    x.value[]
end

@noinline function initialize!(x::LazyInitialized, constructor::F1, hook::F2) where {F1, F2}
    status = Threads.atomic_cas!(x.guard, 0, 1)
    if status == 0
        try
          x.value[] = constructor()
          x.guard[] = 2
        catch
          x.guard[] = 0
          rethrow()
        end

        if hook !== nothing
          hook()
        end
    else
        yield()
    end
    return
end


text:
# ---
# title: Colorbar tweaking
# cover: assets/colorbar_tweaking.png
# description: Setting colorbar attributes.
# ---

using AlgebraOfGraphics, CairoMakie
set_aog_theme!() #src

# To tweak the position and appearance of the colorbar, simply use the `colorbar` keyword when plotting. For example

df = (x=rand(100), y=rand(100), z=rand(100))
plt = data(df) * mapping(:x, :y, color=:z)
draw(plt)

#

fg = draw(plt, colorbar=(position=:top, size=25))

# save cover image #src
mkpath("assets") #src
save("assets/colorbar_tweaking.png", fg) #src


text:
#############################################################################
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#
# Author: Louis Luangkesorn
# Date: February 26, 2015
#
# JuMP implementation of the multicommodity transportation model
#   AMPL: A Modeling Language for Mathematical Programming, 2nd ed
#   by Robert Fourer, David Gay, and Brian W. Kernighan
#   4-1 multi.mod and multi.dat
#############################################################################

using Compat

function PrintSolution(status, Trans, ORIG, DEST, PROD)
    println("RESULTS:")
    if status == :Optimal
      for i = 1:length(ORIG)
        for j = 1:length(DEST)
          for p = 1:length(PROD)
            print(" $(PROD[p]) $(ORIG[i]) $(DEST[j]) = $(getvalue(Trans[i,j, p])) \t")
          end
          println()
        end
      end
    else
        println("  No solution")
    end
    println("")
end

using JuMP

# PARAMETERS

orig = ["GARY" "CLEV" "PITT"]
dest = ["FRA" "DET" "LAN" "WIN" "STL" "FRE" "LAF"]
prod = ["bands" "coils" "plate"]
numorig = length(orig)
numdest = length(dest)
numprod = length(prod)

# supply(prod, orig) amounts available at origins
supply = [[ 400    700    800]
          [ 800   1600   1800]
          [ 200    300    300]]

# demand(prod, dest) amounts required at destinations
demand = [[  300   300   100    75   650   225   250]
          [ 500   750   400   250   950   850   500]
          [ 100   100     0    50   200   100   250]]

# limit(orig, dest) of total units from any origin to destination
defaultlimit = float(625)

limit = [[defaultlimit::AbstractFloat for j=1:numdest] for i=1:numorig]

# cost(dest, orig, prod) Shipment cost per unit
cost = reshape([[[  30,   10,   8 ,  10,   11 ,  71,    6];
         [  22,    7,   10,    7,   21,   82,   13];
         [  19,   11,   12,   10,   25,   83,   15]];
        [[  39,   14,   11,   14,   16,   82,    8];
         [  27,    9,   12,    9,   26,   95,   17];
         [  24,   14,   17,   13,   28,   99,   20]];
        [[  41,   15,   12,   16,   17,   86,    8];
         [  29,    9,   13,    9,   28,   99,   18];
         [  26,   14,   17,   13,   31,  104,   20]]],
               7, 3, 3)

#  DECLARE MODEL

multi = Model()

#  VARIABLES

@variable(multi, Trans[1:numorig, 1:numdest, 1:numprod] >= 0)

#  OBJECTIVE

length(cost)

@objective(multi, Max,
              sum{sum{sum{cost[j, i, p] * Trans[i,j, p],
                        i=1:numorig}, j=1:numdest}, p=1:numprod})

#  CONSTRAINTS

# Supply constraint
@constraint(multi, xyconstr[i=1:numorig, p=1:numprod],
               sum{Trans[i,j,p], j=1:numdest} == supply[p,i])

# Demand constraint
@constraint(multi, xyconstr[j=1:numdest, p=1:numprod],
               sum{Trans[i,j,p], i=1:numorig} == demand[p,j])

# Total shipment constraint
@constraint(multi, xyconstr[i=1:numorig, j=1:numdest],
               sum{Trans[i,j,p], p=1:numprod} - limit[i][j] <= 0)
limit[2][3]
status = solve(multi)
if status == :Optimal
  result = Trans
  PrintSolution(status, Trans, orig, dest, prod)
else
  result = status
end


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ b65c67ec-b79f-4f0e-85e6-78ff22b279d4
using HypertextLiteral

# ╔═╡ a8c1e0d2-3604-4e1d-a87c-c8f5b86b79ed
md"""
# MultiCheckBox
"""

# ╔═╡ c8350f43-0d30-45d0-873b-ff56c5801ac1
md"""
## Definition
"""

# ╔═╡ 79c4dc76-efa2-4b7f-99ea-243e3a17f81c
function skip_as_script(m::Module)
	if isdefined(m, :PlutoForceDisplay)
		return m.PlutoForceDisplay
	else
		isdefined(m, :PlutoRunner) && parentmodule(m) == Main
	end
end

# ╔═╡ 499ca710-1a50-4aa1-87d8-d213416e8e30
if skip_as_script(@__MODULE__)
	import Pkg
	Pkg.activate(Base.current_project(@__DIR__))
	Text("Project env active")
end

# ╔═╡ 631c14bf-e2d3-4a24-8ddc-095a3dab80ef
import AbstractPlutoDingetjes.Bonds

# ╔═╡ c38de38d-e900-4309-a9f6-1392af2f245b
subarrays(x) = (
	x[collect(I)]
	for I in Iterators.product(Iterators.repeated([true,false],length(x))...) |> collect |> vec
)

# ╔═╡ 430e2c1a-832f-11eb-024a-13e3989fd7c2
begin
	export MultiCheckBox
    local result = begin
    """
    ```julia
    MultiCheckBox(options::Vector; [default::Vector], [orientation ∈ [:row, :column]], [select_all::Bool])
    ```

    A group of checkboxes - the user can choose which of the `options` to return.
    The value returned via `@bind` is a list containing the currently checked items.

    See also: [`MultiSelect`](@ref).

    `options` can also be an array of pairs `key::Any => value::String`. The `key` is returned via `@bind`; the `value` is shown.

    # Keyword arguments
    - `defaults` specifies which options should be checked initally.
    - `orientation` specifies whether the options should be arranged in `:row`'s `:column`'s.
    - `select_all` specifies whether or not to include a "Select All" checkbox.

    # Examples
    ```julia
    @bind snacks MultiCheckBox(["🥕", "🐟", "🍌"]))

    if "🥕" ∈ snacks
        "Yum yum!"
    end
    ```

    ```julia
    @bind functions MultiCheckBox([sin, cos, tan])

    [f(0.5) for f in functions]
    ```

    ```julia
    @bind snacks MultiCheckBox(["🥕" => "🐰", "🐟" => "🐱", "🍌" => "🐵"]; default=["🥕", "🍌"])
    ```

    ```julia
    @bind animals MultiCheckBox(["🐰", "🐱" , "🐵", "🐘", "🦝", "🐿" , "🐝",  "🐪"]; orientation=:column, select_all=true)
    ```
    """
    struct MultiCheckBox{BT,DT}
        options::AbstractVector{Pair{BT,DT}}
        default::Union{Missing,AbstractVector{BT}}
        orientation::Symbol
        select_all::Bool
    end
    end

    MultiCheckBox(options::AbstractVector{<:Pair{BT,DT}}; default=missing, orientation=:row, select_all=false) where {BT,DT} = MultiCheckBox(options, default, orientation, select_all)

    MultiCheckBox(options::AbstractVector{BT}; default=missing, orientation=:row, select_all=false) where BT = MultiCheckBox{BT,BT}(Pair{BT,BT}[o => o for o in options], default, orientation, select_all)

    const js = read(joinpath(@__DIR__, "..", "assets", "multicheckbox.js"), String)
    const css = read(joinpath(@__DIR__, "..", "assets", "multicheckbox.css"), String)

    function Base.show(io::IO, m::MIME"text/html", mc::MultiCheckBox)
        @assert mc.orientation == :column || mc.orientation == :row "Invalid orientation $(mc.orientation). Orientation should be :row or :column"

        defaults = coalesce(mc.default, [])

		# Old:
		# checked = [k in defaults for (k,v) in mc.options]
		#
		# More complicated to fix https://github.com/JuliaPluto/PlutoUI.jl/issues/106
		defaults_copy = copy(defaults)
		checked = [
			let
				i = findfirst(isequal(k), defaults_copy)
				if i === nothing
					false
				else
					deleteat!(defaults_copy, i)
					true
				end
			end
		for (k,v) in mc.options]

        show(io, m, @htl("""
        <plj-multi-checkbox style="flex-direction: $(mc.orientation);"></plj-multi-checkbox>
        <script type="text/javascript">
            const labels = $([string(v) for (k,v) in mc.options]);
            const values = $(1:length(mc.options));
            const checked = $(checked);
            const includeSelectAll = $(mc.select_all);
            $(HypertextLiteral.JavaScript(js))
        </script>
        <style type="text/css">$(css)</style>
        """))
    end

    Base.get(select::MultiCheckBox) = Bonds.initial_value(select)
    Bonds.initial_value(select::MultiCheckBox{BT,DT}) where {BT,DT} =
        ismissing(select.default) ? BT[] : select.default
    Bonds.possible_values(select::MultiCheckBox) =
        subarrays((string(i) for i in 1:length(select.options)))

    function Bonds.transform_value(select::MultiCheckBox{BT,DT}, val_from_js) where {BT,DT}
        # val_from_js will be a vector of Strings, but let's allow Integers as well, there's no harm in that
        @assert val_from_js isa Vector

        val_nums = (
            v isa Integer ? v : tryparse(Int64, v)
            for v in val_from_js
        )

        BT[select.options[v].first for v in val_nums]
    end

    function Bonds.validate_value(select::MultiCheckBox, val)
        val isa Vector && all(val_from_js) do v
            val_num = v isa Integer ? v : tryparse(Int64, v)
            1 ≤ val_num ≤ length(select.options)
        end
    end
    result
end

# ╔═╡ 8bfaf4c8-557d-433e-a228-aac493746efc
@bind animals MultiCheckBox(["🐰", "🐱" , "🐵", "🐘", "🦝", "🐿" , "🐝",  "🐪"]; orientation=:column, select_all=true)

# ╔═╡ 8e9f3962-d86c-4e07-b5d3-f31ee5361ca2
animals

# ╔═╡ d8613c4f-6936-4fb2-9b9f-acf34377091f
@bind funcs MultiCheckBox([sin, cos, tan]; default=[sin, tan])

# ╔═╡ 92fe974e-cc6c-4387-8c6e-85813b222f25
funcs

# ╔═╡ 1dd8a041-60c0-4b26-8207-3dc9cca0d3eb
[f(0.5) for f in funcs]

# ╔═╡ a8a7e90d-8bbf-4ab6-90a8-24a10885fb0a
@bind animals2 MultiCheckBox(["\"🐰\\\"", "🐱" , "🐵", "🐘", "🦝", "🐿" , "🐝",  "🐪"]; orientation=:column, default=["🐿" , "🐝"])

# ╔═╡ 6123cf6d-fc29-4a8f-a5a1-4366cc6457b6
animals2

# ╔═╡ 60183ad1-4919-4402-83fb-d53b86dda0a6
MultiCheckBox(["🐰 &&\\a \$\$", "🐱" , "🐵", "🐘", "🦝", "🐿" , "🐝",  "🐪"])

# ╔═╡ ad6dcdec-2fc9-45d2-8828-62ac857b4afa
@bind snacks MultiCheckBox(
	["🐱" => "🐝", "🐵" => "🦝", "🐱" => "🐿"];
	default=["🐱", "🐱"]
)

# ╔═╡ abe4c3e0-6e1e-4e26-a4fa-bd60f31c1a4c
snacks

# ╔═╡ Cell order:
# ╟─a8c1e0d2-3604-4e1d-a87c-c8f5b86b79ed
# ╠═8bfaf4c8-557d-433e-a228-aac493746efc
# ╠═8e9f3962-d86c-4e07-b5d3-f31ee5361ca2
# ╠═d8613c4f-6936-4fb2-9b9f-acf34377091f
# ╠═92fe974e-cc6c-4387-8c6e-85813b222f25
# ╠═1dd8a041-60c0-4b26-8207-3dc9cca0d3eb
# ╠═a8a7e90d-8bbf-4ab6-90a8-24a10885fb0a
# ╠═6123cf6d-fc29-4a8f-a5a1-4366cc6457b6
# ╠═60183ad1-4919-4402-83fb-d53b86dda0a6
# ╠═abe4c3e0-6e1e-4e26-a4fa-bd60f31c1a4c
# ╠═ad6dcdec-2fc9-45d2-8828-62ac857b4afa
# ╟─c8350f43-0d30-45d0-873b-ff56c5801ac1
# ╟─79c4dc76-efa2-4b7f-99ea-243e3a17f81c
# ╟─499ca710-1a50-4aa1-87d8-d213416e8e30
# ╠═631c14bf-e2d3-4a24-8ddc-095a3dab80ef
# ╠═b65c67ec-b79f-4f0e-85e6-78ff22b279d4
# ╠═430e2c1a-832f-11eb-024a-13e3989fd7c2
# ╠═c38de38d-e900-4309-a9f6-1392af2f245b


text:
####
#### Flags
####

export StartFlag, StopFlag, FlagPair, FlagSet
export GreedyType, ScopeType
export TokenStream

abstract type AbstractGrepType end
struct ScopeType <: AbstractGrepType end
struct GreedyType <: AbstractGrepType end

abstract type AbstractFlag{GT<:AbstractGrepType} end

"""
    StartFlag(flag::String,
              flag_boundaries_left::Vector{String},
              flag_boundaries_right::Vector{String})

A "start" flag that BetweenFlags looks for to denote
the start position of a given scope.
"""
struct StartFlag{GT} <: AbstractFlag{GT}
  flag::String
  flag_boundaries_left::Vector{String}
  flag_boundaries_right::Vector{String}
  trigger::Dict{Tuple,String}
end
function StartFlag(flag::S,
  flag_boundaries_left=S[],
  flag_boundaries_right=S[];
  grep_type::AbstractGrepType=GreedyType()
  ) where {S<:AbstractString, GT<:AbstractGrepType}
  trigger = Dict{Tuple,String}()
  for left in flag_boundaries_left
    for right in flag_boundaries_right
      trigger[(left,right)] = string(left, flag, right)
    end
  end
  isempty(trigger) && (trigger[("","")] = flag)
  return StartFlag{typeof(grep_type)}(flag, flag_boundaries_left,
        flag_boundaries_right, trigger)
end

"""
    StopFlag(flag::String,
              flag_boundaries_left::Vector{String},
              flag_boundaries_right::Vector{String})

A "stop" flag that BetweenFlags looks for to denote
the stop position of a given scope.
"""
struct StopFlag{GT} <: AbstractFlag{GT}
  flag::String
  flag_boundaries_left::Vector{String}
  flag_boundaries_right::Vector{String}
  trigger::Dict{Tuple,String}
end
function StopFlag(flag::S,
  flag_boundaries_left=S[],
  flag_boundaries_right=S[];
  grep_type::AbstractGrepType=GreedyType()
  ) where {S<:AbstractString}
  trigger = Dict{Tuple,String}()
  for left in flag_boundaries_left
    for right in flag_boundaries_right
      trigger[(left,right)] = string(left, flag, right)
    end
  end
  isempty(trigger) && (trigger[("","")] = flag)
  return StopFlag{typeof(grep_type)}(flag, flag_boundaries_left,
        flag_boundaries_right, trigger)
end

# TODO: resolve inconsistency between
# grep_type of flag and flag pair, as
# these are easily set different from
# each other.
grep_type(::AbstractFlag{T}) where T = T

####
#### FlagPair
####

"""
    FlagPair(start::Flag, stop::Flag)

A flag pair that defines the start and stop of
the substring of interest.

```julia
julia>
using BetweenFlags
# find: ["\\nfunction", " function", ";function"]
start_flag = StartFlag("function",
                  ["\\n", "\\s", ";"],
                  ["\\n", "\\s"])
# find: ["\\nend", " end", ";end"]
stop_flag = StopFlag("end",
                 ["\\n", "\\s", ";"],
                 ["\\n", "\\s", ";"])
flag_pair = FlagPair{ScopeType}(start_flag, stop_flag)
```
"""
struct FlagPair{FPT}
  start::StartFlag
  stop::StopFlag
  ID::String
  function FlagPair{FPT}(
    start::StartFlag,
    stop::StopFlag
    ) where {FPT<:Union{ScopeType,GreedyType}}
    return new{FPT}(start, stop, start.flag*"-"*stop.flag)
  end
end
FlagPair(start, stop) = FlagPair{GreedyType}(start, stop)

grep_type(::FlagPair{T}) where T = T

struct FlagSet{FP}
  flag_pairs::Vector{FP}
end

"""
    TokenStream

A token stream, containing
 - `flag_set` a `FlagSet`
 - `code` a string of code to be tokenized
 - `token_stream` a token stream (dict)
    whose keys are the flag IDs,
    and whose values are vectors of `Int`s
    indicating the nestedness of the flag.
"""
struct TokenStream{FS,C,TS}
  flag_set::FS
  code::C
  token_stream::TS
end


text:
module numder
using LinearAlgebra

include("num_diff.jl")
include("noisy_num_diff.jl")

end # numder

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
Base.@kwdef mutable struct ProblemDesc
    name::String
    params = []
    fn = nothing
    bounds = nothing
    pfront = nothing
    pfront_image = nothing
end

function problem_gui(problem_desc)
    problem_pfront_image_id = nothing

    function ()
        problems = ["ZDT1", "ZDT2", "ZDT3", "ZDT4", "ZDT6",
            "DTLZ1", "DTLZ2", "DTLZ3", "DTLZ4", "DTLZ5", "DTLZ6"]
        @cstatic item_current = "ZDT1" begin
            if CImGui.BeginCombo("Problem", item_current)
                for n = 0:length(problems)-1
                    is_selected = item_current == problems[n+1]
                    if CImGui.Selectable(problems[n+1], is_selected)
                        item_current = problems[n+1]
                        update_problem_desc(problem_desc, problems[n+1])
                    end
                    # set the initial focus when opening the combo (scrolling + for keyboard navigation support in the upcoming navigation branch)
                    is_selected && CImGui.SetItemDefaultFocus()
                end
                CImGui.EndCombo()
            end
        end
        for param_desc in problem_desc.params
            eval(param_gui(problem_desc, param_desc))
        end
        if problem_pfront_image_id !== nothing
            ImGui_ImplOpenGL3_DestroyImageTexture(problem_pfront_image_id)
        end
        if problem_desc.pfront_image !== nothing
            img = problem_desc.pfront_image
            _, img_width, img_height = size(img)
            problem_pfront_image_id = ImGui_ImplOpenGL3_CreateImageTexture(img_width, img_height)
            ImGui_ImplOpenGL3_UpdateImageTexture(problem_pfront_image_id, img, img_width, img_height)
            CImGui.Image(Ptr{Cvoid}(problem_pfront_image_id), (img_width, img_height))
        end
    end
end

function get_problem_desc(problem_name)
    ProblemDesc(name = problem_name,
                params = [ParamDesc(name="D", type="Int",
                                    value=30, default_value=30,
                                    desc="Dimension",
                                    min=2, max=100),
                          ParamDesc(name="n", type="Int",
                                    value=100, default_value=100,
                                    desc="Number of solutions",
                                    min=2, max=1000)])
end

function get_problem_param(problem_desc, param_name)
    for param in problem_desc.params
        if param.name == param_name
            return param.value === nothing ? param.default_value : param.value
        end
    end
end

function update_problem_desc(problem_desc, problem_name, param_name=nothing, param_value=nothing)
    # println("update_problem_desc ", problem_name, " ", param_name, " ", param_value)
    should_update_fn = problem_desc.fn === nothing ? true : false
    if problem_name !== nothing && problem_desc.name != problem_name
        new_problem_desc = get_problem_desc(problem_name)
        problem_desc.name = new_problem_desc.name
        problem_desc.params = new_problem_desc.params
        should_update_fn = true
    else
        for param in problem_desc.params
            if param.name == param_name
                if param.value != param_value
                    should_update_fn = true
                    if param.max !== nothing && param_value > param.max
                        param_value = param.max
                    end
                    if param.min !== nothing && param_value < param.min
                        param_value = param.min
                    end
                    param.value = param_value
                end
            end
        end
    end
    if should_update_fn
        if problem_desc.name in ["ZDT1", "ZDT2", "ZDT3", "ZDT4", "ZDT6"]
            fn, bounds, pfront = eval(Expr(:call, Symbol(problem_desc.name),
                                           get_problem_param(problem_desc, "D"),
                                           get_problem_param(problem_desc, "n")))
            problem_desc.fn = fn
            problem_desc.bounds = bounds
            problem_desc.pfront = pfront
            size = get_figure_size()
            pfront_fig = plot(
                objectives(pfront)[:, 1],
                objectives(pfront)[:, 2],
                seriestype = :scatter,
                legend = false,
                width = size,
                height = size
            )
            pfront_image = plot_image(pfront_fig)
            pfront_image = fit_image(pfront_image, size)
            problem_desc.pfront_image = image_uint8(pfront_image)
        end
    end
end


text:
g = Multigraph(6)
add_edge!(g, 1, 3)
add_edge!(g, 2, 4)
add_edge!(g, 3, 4)
add_edge!(g, 3, 5)
add_edge!(g, 4, 6)
ps = [0//1 for i = 1:6]
v_t = [SpiderType.In, SpiderType.In, SpiderType.X, SpiderType.Z, SpiderType.Out, SpiderType.Out]
zxd = ZXDiagram(g, v_t, ps)
zxg1 = ZXGraph(zxd)
@test outneighbors(zxg1, 1) == inneighbors(zxg1, 1)
@test !ZXCalculus.is_hadamard(zxg1, 2, 4) && !ZXCalculus.is_hadamard(zxg1, 4, 6)
@test add_edge!(zxg1, 1, 1)
@test !add_edge!(zxg1, 2, 4)
@test !add_edge!(zxg1, 7, 8)
@test sum([ZXCalculus.is_hadamard(zxg1, src(e), dst(e)) for e in edges(zxg1.mg)]) == 3
replace!(Rule{:b}(), zxd)
zxg2 = ZXGraph(zxd)
@test !ZXCalculus.is_hadamard(zxg2, 5, 8) && !ZXCalculus.is_hadamard(zxg2, 1, 7)

zxd = ZXDiagram(2)
push_gate!(zxd, Val(:H), 1)
push_gate!(zxd, Val(:CNOT), 2, 1)
zxg = ZXGraph(zxd)

zxg3 = ZXGraph(ZXDiagram(3))
ZXCalculus.add_global_phase!(zxg3, ZXCalculus.Phase(1//4))
ZXCalculus.add_power!(zxg3, 3)
@test ZXCalculus.scalar(zxg3) == Scalar(3, 1//4)
@test degree(zxg3, 1) == indegree(zxg3, 1) == outdegree(zxg3, 1)
@test ZXCalculus.qubit_loc(zxg3, 1) == ZXCalculus.qubit_loc(zxg3, 2)
@test ZXCalculus.column_loc(zxg3, 1) == 1//1
@test ZXCalculus.column_loc(zxg3, 2) == 3//1


text:
"""
Opaque Encoder struct
"""
type OpusEncoder  # type not immutable so that finalizer can be applied
    v::Ptr{Void}
    fs::Int32
    channels::Cint

    """
    Create new OpusEncoder object with given samplerate and channels
    """
    function OpusEncoder(samplerate, channels; application = OPUS_APPLICATION_AUDIO)
        errorptr = Ref{Cint}(0);
        # Create new encoder object with the given samplerate and channel
        encptr = ccall((:opus_encoder_create,libopus), Ptr{Void}, (Int32, Cint, Cint, Ref{Cint}), samplerate, channels, application, errorptr)
        err = errorptr[]
        enc = new(encptr, samplerate, channels)
        if err != OPUS_OK
            error("opus_encoder_create() failed: $(OPUS_ERROR_MESSAGE_STRS[err])")
        end

        # Register finalizer to cleanup this encoder
        finalizer(enc,x -> ccall((:opus_encoder_destroy,libopus),Void,(Ptr{Void},),x.v))
        return enc
    end
end

function encode_frame(enc::OpusEncoder, data::Vector{Float32})
    frame_len = div(length(data),enc.channels)
    if !(frame_len in [120, 240, 480, 960, 1920, 2880])
        error("Invalid frame length of $(length(data)/enc.channels) samples")
    end
    packet = Vector{UInt8}(length(data)*4*enc.channels)

    #print("opus_encode_float: ")
    num_bytes = ccall((:opus_encode_float,libopus), Cint, (Ptr{Void}, Ptr{Float32}, Cint, Ptr{UInt8}, Int32),
                        enc.v, data, frame_len, packet, length(packet))
    if num_bytes < 0
        error("opus_encode_float() failed: $(OPUS_ERROR_MESSAGE_STRS[num_bytes])")
    end
    return copy(packet[1:num_bytes])
end

"""
Array goes in, packets of bytes come out
"""
function encode(enc::OpusEncoder, audio::Array{Float32}; chunksize=960)
    if size(audio, 2) != enc.channels
        error("Audio data must have the same number of channels as encoder!")
    end

    # Do we have multiple channels?  If so, flatten audio
    if size(audio, 2) > 1
        audio = reshape(audio', (prod(size(audio)),))
    end

    packets = Vector{Vector{UInt8}}()
    # Split audio up into chunks of size chunksize
    for chunk_idx in 1:div(length(audio),enc.channels*chunksize)
        chunk = audio[(chunk_idx-1)*chunksize*enc.channels + 1:chunk_idx*chunksize*enc.channels]
        encoded_chunk = encode_frame(enc, chunk)
        push!(packets, encoded_chunk)
    end

    # Encode last chunk, zero-padding if we need to
    leftover_samples = div(length(audio),enc.channels) % chunksize
    if leftover_samples != 0
        chunk = audio[length(audio) - leftover_samples*enc.channels + 1:length(audio)]
        chunk = vcat(chunk, zeros(Float32, chunksize - leftover_samples))
        encoded_chunk = encode_frame(enc, chunk)
        push!(packets, encoded_chunk)
    end
    return packets
end

function save(output::Union{File{format"OPUS"},AbstractString,IO}, audio::Array{Float32}, fs; chunksize=960)
    # Encode the audio into packets
    enc = OpusEncoder(fs, size(audio, 2))
    packets = encode(enc, audio; chunksize=chunksize)

    # Calculate granule positions for each packet.  We mark the first two packets as
    # "header packets" by setting their granulepos to zero, which forces them into
    # their own ogg pages
    granulepos = vcat(0, 0, collect(1:length(packets))*div(chunksize,enc.channels))

    # Save it out into an Ogg file with an OpusHead and OpusTags
    opus_head = OpusHead(fs, size(audio, 2))
    insert!(packets, 1, opus_head)
    opus_tags = OpusTags()
    insert!(packets, 2, opus_tags)
    Ogg.save(output, packets, granulepos)
end

function save(output::Union{File{format"OPUS"},AbstractString,IO}, audio::Array, fs; chunksize=960)
    return save(output, map(Float32, audio), fs, chunksize=chunksize)
end


text:
module LLVMCalls

using Static: Zero, One, True, False, StaticInt, StaticBool, lt, gt
using Base: FastMath
import IfElse: ifelse

using Base: HWReal

export Vec, VecUnroll, MM, Mask

const LLVMCALL = GlobalRef(Base, :llvmcall)
const FloatingTypes = Union{Float32, Float64} # Float16

# const SignedHW = Union{Int8,Int16,Int32,Int64,Int128}
# const UnsignedHW = Union{UInt8,UInt16,UInt32,UInt64,UInt128}
const SignedHW = Union{Int8,Int16,Int32,Int64}
const UnsignedHW = Union{UInt8,UInt16,UInt32,UInt64}
const IntegerTypesHW = Union{SignedHW,UnsignedHW}
const IntegerTypes = Union{StaticInt,IntegerTypesHW}

struct Bit; data::Bool; end # Dummy for Ptr
const Boolean = Union{Bit,Bool}
const NativeTypesExceptBit = Union{Bool,HWReal}
const NativeTypes = Union{NativeTypesExceptBit, Bit}

const _Vec{W,T<:Number} = NTuple{W,Core.VecElement{T}}

abstract type AbstractSIMD{W,T <: Union{<:StaticInt,NativeTypes}} <: Real end
abstract type AbstractSIMDVector{W,T} <: AbstractSIMD{W,T} end
struct VecUnroll{N,W,T,V<:Union{NativeTypes,AbstractSIMD{W,T}}} <: AbstractSIMD{W,T}
  data::Tuple{V,Vararg{V,N}}
  @inline (VecUnroll(data::Tuple{V,Vararg{V,N}})::VecUnroll{N,W,T,V}) where {N,W,T,V<:AbstractSIMD{W,T}} = new{N,W,T,V}(data)
  @inline (VecUnroll(data::Tuple{T,Vararg{T,N}})::VecUnroll{N,T,T}) where {N,T<:NativeTypes} = new{N,1,T,T}(data)
end
# @inline VecUnroll(data::Tuple) = VecUnroll(promote(data))
const VecOrScalar = Union{AbstractSIMDVector,NativeTypes}
const NativeTypesV = Union{AbstractSIMD,NativeTypes,StaticInt}
const IntegerTypesV = Union{AbstractSIMD{<:Any,<:IntegerTypes},IntegerTypesHW}
struct Vec{W,T} <: AbstractSIMDVector{W,T}
  data::NTuple{W,Core.VecElement{T}}
  @inline Vec{W,T}(x::NTuple{W,Core.VecElement{T}}) where {W,T<:NativeTypes} = new{W,T}(x)
  @generated function Vec(x::Tuple{Core.VecElement{T},Vararg{Core.VecElement{T},_W}}) where {_W,T<:NativeTypes}
    W = _W + 1
    vtyp = Expr(:curly, :Vec, W, T)
    Expr(:block, Expr(:meta,:inline), Expr(:(::), Expr(:call, vtyp, :x), vtyp))
  end
end

@generated mask_type(::StaticInt{W}) where {W} = Symbol(:UInt, clamp(nextpow2(W), 8, 128))
@inline data(vu::VecUnroll) = getfield(vu, :data)
abstract type AbstractMask{W,U<:Union{UnsignedHW,UInt128}} <: AbstractSIMDVector{W,Bit} end
struct Mask{W,U} <: AbstractMask{W,U}
  u::U
  @inline function Mask{W,U}(u::Unsigned) where {W,U} # ignores U...
    U2 = mask_type(StaticInt{W}())
    new{W,U2}(u % U2)
  end
end
struct EVLMask{W,U} <: AbstractMask{W,U}
  u::U
  evl::UInt32
  @inline function EVLMask{W,U}(u::Unsigned, evl) where {W,U} # ignores U...
    U2 = mask_type(StaticInt{W}())
    new{W,U2}(u % U2, evl % UInt32)
  end
end
const AnyMask{W} = Union{AbstractMask{W},VecUnroll{<:Any,W,Bit,<:AbstractMask{W}}}
@inline Mask{W}(u::U) where {W,U<:Unsigned} = Mask{W,U}(u)
@inline EVLMask{W}(u::U, i) where {W,U<:Unsigned} = EVLMask{W,U}(u, i)
@inline Mask{1}(b::Bool) = b
@inline EVLMask{1}(b::Bool, i) = b
@inline Mask(m::EVLMask{W,U}) where {W,U} = Mask{W,U}(getfield(m,:u))
# Const prop is good enough; added an @inferred test to make sure.
# Removed because confusion can cause more harm than good.

@inline Base.broadcastable(v::AbstractSIMDVector) = Ref(v)

Vec{W,T}(x::Vararg{NativeTypes,W}) where {W,T<:NativeTypes} = Vec(ntuple(w -> Core.VecElement{T}(x[w]), Val{W}()))
Vec{1,T}(x::Union{Float32,Float64}) where {T<:NativeTypes} = T(x)
Vec{1,T}(x::Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Bool}) where {T<:NativeTypes} = T(x)

@inline Base.length(::AbstractSIMDVector{W}) where W = W
@inline Base.size(::AbstractSIMDVector{W}) where W = (W,)
@inline Base.eltype(::AbstractSIMD{W,T}) where {W,T} = T
@inline Base.conj(v::AbstractSIMDVector) = v # so that things like dot products work.
@inline Base.adjoint(v::AbstractSIMDVector) = v # so that things like dot products work.
@inline Base.transpose(v::AbstractSIMDVector) = v # so that things like dot products work.

# Not using getindex/setindex as names to emphasize that these are generally treated as single objects, not collections.
@generated function extractelement(v::Vec{W,T}, i::I) where {W,I <: IntegerTypesHW,T}
  typ = LLVM_TYPES[T]
  instrs = """
          %res = extractelement <$W x $typ> %0, i$(8sizeof(I)) %1
          ret $typ %res
      """
  call = :($LLVMCALL($instrs, $T, Tuple{_Vec{$W,$T},$I}, data(v), i))
  Expr(:block, Expr(:meta, :inline), call)
end
@generated function insertelement(v::Vec{W,T}, x::T, i::I) where {W,I <: IntegerTypesHW,T}
  typ = LLVM_TYPES[T]
  instrs = """
          %res = insertelement <$W x $typ> %0, $typ %1, i$(8sizeof(I)) %2
          ret <$W x $typ> %res
      """
  call = :(Vec($LLVMCALL($instrs, _Vec{$W,$T}, Tuple{_Vec{$W,$T},$T,$I}, data(v), x, i)))
  Expr(:block, Expr(:meta, :inline), call)
end
@inline (v::AbstractSIMDVector)(i::IntegerTypesHW) = extractelement(v, i - one(i))
@inline (v::AbstractSIMDVector)(i::Integer) = extractelement(v, Int(i) - 1)
Base.@propagate_inbounds (vu::VecUnroll)(i::Integer, j::Integer) = getfield(vu, :data)[j](i)

@inline Base.Tuple(v::Vec{W}) where {W} = ntuple(v, Val{W}())

# Use with care in function signatures; try to avoid the `T` to stay clean on Test.detect_unbound_args

@inline Base.copy(v::AbstractSIMDVector) = v
@inline data(v) = v
@inline data(v::Vec) = getfield(v, :data)

function Base.show(io::IO, v::AbstractSIMDVector{W,T}) where {W,T}
  name = typeof(v)
  print(io, "$(name)<")
  for w ∈ 1:W
    print(io, repr(extractelement(v, w-1)))
    w < W && print(io, ", ")
  end
  print(io, ">")
end
Base.bitstring(m::AbstractMask{W}) where {W} = bitstring(data(m))[end-W+1:end]
function Base.show(io::IO, m::AbstractMask{W}) where {W}
  bits = data(m)
  if m isa EVLMask
    print(io, "EVLMask{$W,Bit}<")
  else
    print(io, "Mask{$W,Bit}<")
  end
  for w ∈ 0:W-1
    print(io, bits & 1)
    bits >>= 1
    w < W-1 && print(io, ", ")
  end
  print(io, ">")
end
function Base.show(io::IO, vu::VecUnroll{N,W,T,V}) where {N,W,T,V}
  println(io, "$(N+1) x $V")
  d = data(vu)
  for n in 1:N+1
    show(io, d[n]);
    n > N || println(io)
  end
end

"""
  The name `MM` type refers to _MM registers such as `XMM`, `YMM`, and `ZMM`.
  `MMX` from the original MMX SIMD instruction set is a [meaningless initialism](https://en.wikipedia.org/wiki/MMX_(instruction_set)#Naming).

  The `MM{W,X}` type is used to represent SIMD indexes of width `W` with stride `X`.
  """
struct MM{W,X,I<:Union{HWReal,StaticInt}} <: AbstractSIMDVector{W,I}
  i::I
  @inline MM{W,X}(i::T) where {W,X,T<:Union{HWReal,StaticInt}} = new{W,X::Int,T}(i)
end
@inline MM(i::MM{W,X}) where {W,X} = MM{W,X}(getfield(i, :i))
@inline MM{W}(i::Union{HWReal,StaticInt}) where {W} = MM{W,1}(i)
@inline MM{W}(i::Union{HWReal,StaticInt}, ::StaticInt{X}) where {W,X} = MM{W,X}(i)
@inline data(i::MM) = getfield(i, :i)

@inline extractelement(i::MM{W,X,I}, j) where {W,X,I} = getfield(i, :i) + (X % I) * (j % I)

function Base.getproperty(::AbstractSIMD, ::Symbol)
  throw(ErrorException("""
  `Base.getproperty` not defined on AbstractSIMD.
  If you wish to work with the data as a tuple, it is recommended to use `Tuple(v)`. Once you have an ordinary tuple, you can access
  individual elements normally. Alternatively, you can index using parenthesis, e.g. `v(1)` indexes the first element.
  Parenthesis are used instead of `getindex`/square brackets because `AbstractSIMD` objects represent a single number, and
  for `x::Number`, `x[1] === x`.

  If you wish to perform a reduction on the collection, the naming convention is prepending the base function with a `v`. These functions
  are not overloaded, because for `x::Number`, `sum(x) === x`. Functions include `vsum`, `vprod`, `vmaximum`, `vminimum`, `vany`, and `vall`.

  If you wish to define a new operation applied to the entire vector, do not define it in terms of operations on the individual eleemnts.
  This will often lead to bad code generation -- bad in terms of both performance, and often silently producing incorrect results!
  Instead, implement them in terms of existing functions defined on `::AbstractSIMD`. Please feel free to file an issue if you would like
  clarification, and especially if you think the function may be useful for others and should be included in `VectorizationBase.jl`.
  """))
end

nextpow2(W::T) where {T<:Base.BitInteger} = one(T) << (T(8sizeof(T)) - leading_zeros(W - one(T)))

@generated function Vec(y::T, x::Vararg{T,_W}) where {T<:NativeTypes, _W}
  W = 1 + _W
  Wfull = nextpow2(W)
  ty = LLVM_TYPES[T]
  init = W == Wfull ? "undef" : "zeroinitializer"
  instrs = ["%v0 = insertelement <$Wfull x $ty> $init, $ty %0, i32 0"]
  Tup = Expr(:curly, :Tuple, T)
  for w ∈ 1:_W
    push!(instrs, "%v$w = insertelement <$Wfull x $ty> %v$(w-1), $ty %$w, i32 $w")
    push!(Tup.args, T)
  end
  push!(instrs, "ret <$Wfull x $ty> %v$_W")
  llvmc = :($LLVMCALL($(join(instrs,"\n")), _Vec{$Wfull,$T}, $Tup))
  push!(llvmc.args, :y)
  for w ∈ 1:_W
    push!(llvmc.args, Expr(:ref, :x, w))
  end
  quote
    $(Expr(:meta,:inline))
    Vec($llvmc)
  end
end
@inline reduce_to_onevec(f::F, vu::VecUnroll) where {F} = ArrayInterface.reduce_tup(f, data(vu))

include("static.jl")
include("llvm_types.jl")
include("lazymul.jl")
include("vector_width.jl")
include("ranges.jl")
include("binary_ops.jl")
include("conversion.jl")
include("masks.jl")
include("intrin_funcs.jl")
include("memory_addr.jl")
include("unary_ops.jl")
include("vbroadcast.jl")
include("vector_ops.jl")
include("nonbroadcastingops.jl")
include("integer_fma.jl")
include("fmap.jl")

end


text:
# Autogenerated wrapper script for x265_jll for i686-w64-mingw32
export libx265, x265

JLLWrappers.@generate_wrapper_header("x265")
JLLWrappers.@declare_library_product(libx265, "libx265.dll")
JLLWrappers.@declare_executable_product(x265)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libx265,
        "bin\\libx265.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        x265,
        "bin\\x265.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
config_file = isempty(ARGS) ? joinpath(@__DIR__, "config", "config.jl") : ARGS[1]
isfile(config_file) || error("Cannot read '$config_file'")
println("Config supplied: '$config_file'")
include(config_file)

using Random
using DelimitedFiles
using Distributions

function process_file(input_file, output_file)
    Random.seed!(0)
    raw = readdlm(input_file, ',')
    num_attributes = length(findall(x -> occursin("@ATTRIBUTE", string(x)), raw[:, 1])) - 2
    id_column = findfirst(x -> occursin("@ATTRIBUTE 'id'", string(x)), raw[:, 1]) - 1
    label_column = findfirst(x -> occursin("@ATTRIBUTE 'outlier'", string(x)), raw[:, 1]) - 1
    data_start_row = findlast(x -> x == "@DATA", raw[:, 1]) + 1
    raw[:, label_column] = map(x -> x == "'yes'" ? :outlier : :inlier, raw[:, label_column])
    data, labels = raw[data_start_row:end, [i for i in 1:size(raw, 2) if i != id_column && i != label_column]], raw[data_start_row:end, label_column]
    data = hcat(data, labels)
    @assert size(data, 2) - 1 == num_attributes
    @assert size(data, 1) == length(labels)
    if ADD_NOISE
        @info "Adding noise."
        data[:, 1:end-1] += rand(NOISE_DIST, size(data, 1), size(data, 2) - 1)
    end
    @info "Saving to '$output_file'."
    writedlm(output_file, data, ',')
    return nothing
end

ADD_NOISE = true
NOISE_DIST = Normal(0, 0.01)
target_versions_semantic = r"withoutdupl_norm_[0-9][0-9].arff"
target_versions_literature = r"withoutdupl_norm"
dataset_dir = normpath(data_raw_input_root)
output_path = normpath(data_input_root)

mkpath(output_path)
@info "Saving processed files to $output_path."

for dataset_class in ["semantic", "literature"]
    for d in raw_data_dirs[dataset_class]
        @info d
        outdir = joinpath(output_path, d)
        isdir(outdir) || mkpath(outdir)
        if dataset_class == "semantic"
            target_files = filter(x -> occursin(target_versions_semantic, x), readdir(joinpath(dataset_dir, dataset_class, d)))
            @assert length(target_files) == 1
            @info "[$(d)] Found $(length(target_files)) files."
            for f in target_files
               process_file(joinpath(dataset_dir, "semantic", d, f), joinpath(outdir, f[1:end-5] * ".csv"))
            end
        else
            target_files = filter(x -> occursin(target_versions_literature, x), readdir(joinpath(dataset_dir, dataset_class, d)))
            if (i = findfirst(x -> occursin("catremoved", x), target_files)) !== nothing
                target_file = target_files[i]
            else
                target_file = first(target_files)
            end
            @info "[$(d)] Resampling '$target_file'."
            process_file(joinpath(dataset_dir, dataset_class, d, target_file), joinpath(outdir, target_file[1:end-5] * ".csv"))
        end
    end
end

@info "Preprocessing done."

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
# Define initialization Function
function MathProgBase.initialize(d::STOev, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad, :Jac, :Hess])
            error("Unsupported feature $feat")
        end
    end
end

# List Available Features
MathProgBase.features_available(d::STOev) = [:Grad, :Jac, :Hess]

# Functions to precompute matrices for linear STO

"Precompute matrices for cost function"
function precompMatJ!(d::STOev, x)

  # Propagate dynamics
  propagateDynamics!(d, x)

  # DEBUG: Timing
  # timeprop = @elapsed propagateDynamics!(d, x)
  # @printf("Time elapsed propag = %.4f mus\n", timeprop*10^6)

  # Compute Matrices S
  computeSlw!(d)

  # DEBUG: Timing
  # timeS = @elapsed computeSlw!(d)
  # @printf("Time elapsed S = %.4f mus\n", timeS*10^6)

end

"Precompute matrices for gradient"
function precompMatGradJ!(d::STOev, x)
  if d.deltafun_prev!= x  # Different point than the previous cost function evaluation
    precompMatJ!(d, x)  # precompute matrices for cost function
    d.deltafun_prev[:] = x

    computeClw!(d)  # precompute matrices for gradient

    # DEBUG: Timing
    # timeC = @elapsed computeClw!(d)  # precompute matrices for gradient
    # @printf("Time elapsed C = %.4f mus\n", timeC*10^6)

  else
    computeClw!(d)  # precompute matrices for gradient

    # DEBUG: Timing
    # timeC = @elapsed computeClw!(d)  # precompute matrices for gradient
    # @printf("Time elapsed C = %.4f mus\n", timeC*10^6)
  end
end

"Precompute matrices for hessian"
function precompMatHessJ!(d::STOev, x)
  if d.deltagrad_prev!= x  # Different point than the previous gradient evaluation
    precompMatGradJ!(d, x)  # precompute matrices for gradient function
    d.deltagrad_prev[:] = x

    computePhilw!(d)  # precompute matrices for hessian

    # DEBUG: Timing
    # timePhi = @elapsed computePhilw!(d)  # precompute matrices for hessian
    # @printf("Time elapsed Phi = %.4f mus\n", timePhi*10^6)

  else
    computePhilw!(d)  # precompute matrices for hessian

    # DEBUG: Timing
    # timePhi = @elapsed computePhilw!(d)  # precompute matrices for hessian
    # @printf("Time elapsed Phi = %.4f mus\n", timePhi*10^6)

  end
end

# Evaluate Cost Function
function MathProgBase.eval_f(d::STOev, x)
  # Increase counter of objective function evaluations
  d.nobjeval += 1

  # Check if the matrices have already been precomputed
  if d.deltafun_prev!= x
    precompMatJ!(d, x)
    d.deltafun_prev[:] = x
  end

  J = (d.x0'*d.S[:,:,1]*d.x0)[1]

  return J
end

# Evaluate Gradient
function MathProgBase.eval_grad_f(d::STOev, grad_f, x)

  # Increase counter of gradient evaluations
  d.ngradeval += 1

  # Check if the matrices have already been precomputed
  if d.deltagrad_prev!= x
    precompMatGradJ!(d, x)
    d.deltagrad_prev[:] = x
  end

  # Update objective and cost function values to check iteration status (DEBUG)
  d.obj = [d.obj; MathProgBase.eval_f(d, x)]
  d.deltaval = [d.deltaval x]

  for i = 1:d.N+1
    grad_f[i] = (d.xpts[:,d.tauIdx[i+1]]'*d.C[:, :, i]*d.xpts[:,d.tauIdx[i+1]])[1]
  end

end

function MathProgBase.eval_hesslag(d::linSTOev, H, x, sigma, mu )

  # Increase counter of Hessian evaluations
  d.nhesseval += 1

  # Check if the matrices have already been precomputed
  if d.deltahess_prev!= x
    precompMatHessJ!(d, x)
    d.deltahess_prev[:] = x
  end

  Htemp = zeros(d.N+1, d.N+1)

  for i = 1:d.N+1
    Htemp[i, i] = 2*(d.xpts[:, d.tauIdx[i+1]]'*d.C[:, :, i]*d.A[:, :, i]*d.xpts[:, d.tauIdx[i+1]])[1]
  end

  for j = 2:d.N+1
      for i = 1:j-1
      Htemp[j, i] = 2*(d.xpts[:, d.tauIdx[j+1]]'*d.C[:, :, j]*d.Phi[:, :, i+1, j+1]*d.A[:, :, i]*d.xpts[:, d.tauIdx[i+1]])[1]
    end
  end

  Htemp *= sigma

  # Assign Elements of Hessian
  H[:] = Htemp[d.IndTril]
end

function MathProgBase.eval_hesslag(d::nlinSTOev, H, x, sigma, mu )

  # Increase counter of Hessian evaluations
  d.nhesseval += 1

  # Check if the matrices have already been precomputed
  if d.deltahess_prev!= x
    precompMatHessJ!(d, x)
    d.deltahess_prev[:] = x
  end

  Htemp = zeros(d.N+1, d.N+1)

  for i = 1:d.N+1
    Htemp[i, i] = 2*(d.xpts[:, d.tauIdx[i+1]]'*d.C[:, :, i]*d.A[:, :, d.tauIdx[i+1]-1]*d.xpts[:, d.tauIdx[i+1]])[1]
  end

  for j = 2:d.N+1
      for i = 1:j-1
      Htemp[j, i] = 2*(d.xpts[:, d.tauIdx[j+1]]'*d.C[:, :, j]*d.Phi[:, :, i+1, j+1]*d.A[:, :, d.tauIdx[i+1]-1]*d.xpts[:, d.tauIdx[i+1]])[1]
    end
  end

  H[:] = sigma*Htemp[d.IndTril]
end

# Constraints
function MathProgBase.eval_g(d::STOev, g, x)
  g[:] = d.Ag*x
end

function MathProgBase.eval_jac_g(d::STOev, J, x)
  J[:] = d.Vg
end

function MathProgBase.jac_structure(d::STOev)
  return d.Ig, d.Jg
end

# Hessian structure for Linear and Nonlinear Systems
function MathProgBase.hesslag_structure(d::STOev)
  return d.Itril, d.Jtril
end

# Linearize Dynamics
function linearizeDyn(nonlin_dyn::Function, nonlin_dyn_deriv::Function, x::Array{Float64,1}, u::Array{Float64,1})

  f = nonlin_dyn(x,u)
  df = nonlin_dyn_deriv(x, u)

  A = [df f-df*x; zeros(1, length(x)+1)]

end

function mergeSortFindIndex(tgrid::Array{Float64, 1}, tau::Array{Float64,1})

  ngrid = length(tgrid)
  N = length(tau)
  tauIdx = Array{Int}(N+2); tauIdx[1] = 1; tauIdx[end]= N + ngrid
  tgridIdx = Array{Int}(ngrid); tgridIdx[1] = 1; tgridIdx[end]= N + ngrid

  # Create merged and sorted time vector with grid and switching times
  ttemp = vcat(tgrid, tau)  # Concatenate grid and tau vector
  tidxtemp = sortperm(ttemp)  # Find permutation vector to sort ttemp
  tvec = ttemp[tidxtemp]    # Create full sorted tvec

  # # Create index of the tau vector elements inside tvec
  for i = 1:N
    tauIdx[i+1] = findfirst(tidxtemp, ngrid + i)
  end

  # # Create index of the tgrid vector elements inside tvec
  for i = 1:ngrid
    tgridIdx[i] = findfirst(tidxtemp, i)
  end

  return tvec, tauIdx, tgridIdx

end

# Propagate dynamic for linear STO
function propagateDynamics!(d::linSTOev, x::Array{Float64,1})

  # Get positive delta
  x = max.(x, 0)

  # Get switching times from delta (tfdelta is the final time we get from the current delta vector)
  tau, d.tfdelta = delta2tau(x, d.t0)

  # Create grid from t0 to tfdelta
  d.tgrid[end] = d.tfdelta
  d.tgrid = min.(d.tgrid, d.tgrid[end])

  # Create merged and sorted time vector with grid and switching times
  d.tvec, d.tauIdx, d.tgridIdx = mergeSortFindIndex(d.tgrid, tau)

  # Get complete delta vector with all intervals
  d.deltacomplete = tau2delta(d.tvec[2:end-1], d.t0, d.tfdelta)

  # The derivative checker in IPOPT will fail in computing the numerical derivatives because of the numerical issues in going to tau formulation and then back to delta formulation. To double check, please uncomment the following line in the case of only 2 points in the grid. The derivative checker should have no errors.

  # d.deltacomplete = x

  # Compute Exponentials required for the computations (for the complete grid)
  Aidx = 1  # Initialize index for current A mode

  # Compute Matrix Exponentials
  for i = 1:d.N+d.ngrid-1  # Iterate over all grid (incl sw times)

    # Verify which mode input applies
    if Aidx <= d.N
      if i>= d.tauIdx[Aidx + 1]
        Aidx += 1
      end
    end

    # Get temporary Matrix relative to current Aidx
    if d.isDiag[Aidx]  # Diagonalizable Matrix -> Compute Fast Exponential
      tempMat = real(d.V[:, :, Aidx]*diagm(exp.(d.D[:,Aidx]*d.deltacomplete[i]))*d.invV[:, :, Aidx])
    else  # Nondiagonalizable Matrix -> Compute Standard Exponential
      # Compute Matrix Exponential
      tempMat = expm([-d.A[:, :, i]'  d.Q;
                      zeros(d.nx, d.nx)  d.A[:, :, i]]*d.deltacomplete[i])
    end

    # Assign \mathcal{E}_i
    d.expMat[:,:,i] = tempMat[d.nx+1:2*d.nx, d.nx+1:2*d.nx]

    # Assign M_k
    d.M[:,:,i] = tempMat[d.nx+1:2*d.nx, d.nx+1:2*d.nx]'*tempMat[1:d.nx, d.nx+1:2*d.nx]

    # Compute next state
    d.xpts[:, i+1] = d.expMat[:,:,i]*d.xpts[:, i]

  end

end

# Propagate dynamics for nonlinear STO
function propagateDynamics!(d::nlinSTOev, x::Array{Float64,1})

  # Get positive delta
  x = max.(x, 0)

  # Get switching times from delta
  tau, d.tfdelta = delta2tau(x, d.t0)

  # Create grid from t0 to tfdelta
  d.tgrid[end] = d.tfdelta
  d.tgrid = min.(d.tgrid, d.tgrid[end])

  # Fit switching times withing the grid
  d.tvec, d.tauIdx, d.tgridIdx = mergeSortFindIndex(d.tgrid, tau)

  # Get complete delta vector with all intervals
  d.deltacomplete = tau2delta(d.tvec[2:end-1], d.t0, d.tfdelta)

  # Propagate Dynamics and Compute Exponentials over the whole grid
  uIdx = 1  # Initialize index for current u

  # Compute Matrix Exponentials
  for i =1:d.N+d.ngrid-1  # Iterate over all grid (incl sw times)

    # Verify which U input applies
    if uIdx <= d.N
      if i>= d.tauIdx[uIdx + 1]
        uIdx += 1
      end
    end

    # Linearize Dynamics
    d.A[:,:,i] = linearizeDyn(d.nonlin_dyn, d.nonlin_dyn_deriv, d.xpts[1:end-1,i], d.uvec[:, uIdx])

    # Compute Matrix Exponential
    tempMat = expm([-d.A[:, :, i]'  d.Q;
                    zeros(d.nx, d.nx)  d.A[:, :, i]]*d.deltacomplete[i])

    # Assign \mathcal{E}_i
    d.expMat[:,:,i] = tempMat[d.nx+1:2*d.nx, d.nx+1:2*d.nx]

    # Assign M_k
    d.M[:,:,i] = tempMat[d.nx+1:2*d.nx, d.nx+1:2*d.nx]'*tempMat[1:d.nx, d.nx+1:2*d.nx]

    # Compute state at the next grid point
    d.xpts[:, i+1] = d.expMat[:, :, i]*d.xpts[:, i]

  end

end

"Low level function for computing S matrices"
function computeSlw!(d::STOev)
  d.S[:,:,d.N+d.ngrid] = d.E
  for i = d.N+d.ngrid-1:-1:1
    d.S[:,:,i] = d.M[:,:,i] + d.expMat[:,:,i]'*d.S[:,:,i+1]*d.expMat[:,:,i]
  end
end

"Low level function for computing C matrices for linear dynamics"
function computeClw!(d::linSTOev)
  for i = 1:d.N+1
    d.C[:, :, i] = d.Q + d.A[:, :, i]'*d.S[:, :, d.tauIdx[i+1]] + d.S[:, :, d.tauIdx[i+1]]*d.A[:, :, i]
  end
end

"Low level function for computing C matrices for nonlinear dynamics"
function computeClw!(d::nlinSTOev)
  for i = 1:d.N+1
    d.C[:, :, i] = d.Q + d.A[:, :, d.tauIdx[i+1]-1]'*d.S[:, :, d.tauIdx[i+1]] + d.S[:, :, d.tauIdx[i+1]]*d.A[:, :, d.tauIdx[i+1]-1]
  end
end

"Low level function for computing Phi matrices"
function computePhilw!(d::STOev)
  for i = 1 : d.N + 2
    d.Phi[:, :, i, i] = eye(d.nx)  # Identity Matrix to Start
    for l = i:d.N + 1  # Iterate over all successive Phi matrices
      d.Phi[:, :, i, l+1] = d.Phi[:, :, i, l]  # Initialize with previous matrix
      for j=d.tauIdx[l]:d.tauIdx[l+1]-1  # Iterate over all points between switching times
        d.Phi[:, :, i, l+1] = d.expMat[:, :, j] * d.Phi[:, :, i, l+1]
      end
    end
  end
end


text:
module Utils

using MosimoBase: Vector3s, Vector3

@inline function acceleration(forces::Vector3s, m::Function)
    N = length(forces)
    a = zeros(Vector3, N)
    @inbounds for i = 1:N
        a[i] = forces[i] / m(i)
    end
    a
end

end