text:
using Mimi
using Distributions
using CSVFiles
using DataFrames

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

function getsim()

        ############################################################################
        # Define random variables (RVs)
        ############################################################################

        if use_only_varMC
            mcs = @defsim begin
                # Climate Variability Parameters
                gvarsd_globalvariabilitystandarddeviation = Normal(0.11294, 0.021644) # ranges from https://github.com/jkikstra/climvar
                rvarsd_regionalvariabilitystandarddeviation["EU"] = Normal(0.35251, 0.019147)
                rvarsd_regionalvariabilitystandarddeviation["USA"] = Normal(0.39171, 0.030440)
                rvarsd_regionalvariabilitystandarddeviation["OECD"] = Normal(0.35272, 0.042962)
                rvarsd_regionalvariabilitystandarddeviation["USSR"] = Normal(0.39679, 0.030228)
                rvarsd_regionalvariabilitystandarddeviation["China"] = Normal(0.26032, 0.059394)
                rvarsd_regionalvariabilitystandarddeviation["SEAsia"] = Normal(0.20009, 0.039896)
                rvarsd_regionalvariabilitystandarddeviation["Africa"] = Normal(0.18202, 0.050171)
                rvarsd_regionalvariabilitystandarddeviation["LatAmerica"] = Normal(0.15618, 0.054156)
                # #The folllowing RVs are in more than one component.  For clarity they are
                # #set here as opposed to below within the blocks of RVs separated by component
                # #so that they are not set more than once.
                # save_savingsrate = TriangularDist(15, 15, 15) # components: MarketDamages, MarketDamagesBurke, NonMarketDamages. GDP, SLRDamages
                #
                # wincf_weightsfactor_sea["USA"] = TriangularDist(.8, .8, .8) # components: SLRDamages, Discountinuity (weights for market and nonmarket are non-stochastic and uniformly 1)
                # wincf_weightsfactor_sea["OECD"] = TriangularDist(.8, .8, .8)
                # wincf_weightsfactor_sea["USSR"] = TriangularDist(.4, .4, .4)
                # wincf_weightsfactor_sea["China"] = TriangularDist(.8, .8, .8)
                # wincf_weightsfactor_sea["SEAsia"] = TriangularDist(.8, .8, .8)
                # wincf_weightsfactor_sea["Africa"] = TriangularDist(.6, .6, .6)
                # wincf_weightsfactor_sea["LatAmerica"] = TriangularDist(.6, .6, .6)
                #
                # automult_autonomouschange = TriangularDist(0.65, 0.65, 0.65)  #components: AdaptationCosts, AbatementCosts
                #
                # #The following RVs are divided into blocks by component
                #
                # # CO2cycle
                # air_CO2fractioninatm = TriangularDist(62, 62, 62)
                # res_CO2atmlifetime = TriangularDist(70, 70, 70)
                # #ccf_CO2feedback = TriangularDist(0, 0, 0) # only usable if lb <> ub
                # ccfmax_maxCO2feedback = TriangularDist(20, 20, 20)
                # stay_fractionCO2emissionsinatm = TriangularDist(0.3,0.3,0.3)
                # ce_0_basecumCO2emissions=TriangularDist(2040000, 2040000, 2040000)
                # a1_percentco2oceanlong=TriangularDist(23.0,	23.0, 23.0)
                # a2_percentco2oceanshort=TriangularDist(26.6, 26.6, 26.6)
                # a3_percentco2land=TriangularDist(27.0, 27.0, 27.0)
                # t1_timeco2oceanlong=TriangularDist(312.5, 312.5, 312.5)
                # t2_timeco2oceanshort=TriangularDist(34.9, 34.9, 34.9)
                # t3_timeco2land=TriangularDist(4.3, 4.3, 4.3)
                # rt_g0_baseglobaltemp=TriangularDist(0.946, 0.946, 0.946)
                # SiBCASA Permafrost
                # perm_sib_af = TriangularDist(1.42609149897258, 2.32504747848815, 1.87556948873036)
                # perm_sib_sens_c_co2 = TriangularDist(28191.1555428869, 35688.3253432574, 31939.7404430722)
                # perm_sib_lag_c_co2 = TriangularDist(35.4926669856915, 87.8949041341782, 61.6937855599349)
                # perm_sib_pow_c_co2 = TriangularDist(0.107020247715729, 0.410961185142816, 0.258990716429273)
                # perm_sib_sens_c_ch4 = TriangularDist(1240.3553299183, 3348.11995329232, 2294.23764160531)
                # perm_sib_lag_c_ch4 = TriangularDist(75.1943160023131, 337.382510123922, 206.288413063117)
                # perm_sib_pow_c_ch4 = TriangularDist(-0.108779283732708, 0.610889007954489, 0.25105486211089)
                #
                # # SulphateForcing
                # d_sulphateforcingbase = TriangularDist(-0.4, -0.4, -0.4)
                # ind_slopeSEforcing_indirect = TriangularDist(-0.4, -0.4, -0.4)
                #
                # # ClimateTemperature
                # frt_warminghalflife = TriangularDist(20, 20, 20)        # from PAGE-ICE v6.2 documentation
                # tcr_transientresponse = TriangularDist(1.8, 1.8, 1.8)   # from PAGE-ICE v6.2 documentation
                # alb_emulator_rand = TriangularDist(-1., 1., 0.)
                # ampf_amplification["EU"] = TriangularDist(1.23, 1.23, 1.23)
                # ampf_amplification["USA"] = TriangularDist(1.32, 1.32, 1.32)
                # ampf_amplification["OECD"] = TriangularDist(1.21, 1.21, 1.21)
                # ampf_amplification["USSR"] = TriangularDist(1.64, 1.64, 1.64)
                # ampf_amplification["China"] = TriangularDist(1.21, 1.21, 1.21)
                # ampf_amplification["SEAsia"] = TriangularDist(1.04, 1.04, 1.04)
                # ampf_amplification["Africa"] = TriangularDist(1.22, 1.22, 1.22)
                # ampf_amplification["LatAmerica"] = TriangularDist(1.04, 1.04, 1.04)
                #
                # # SeaLevelRise
                # s0_initialSL = TriangularDist(0.19, 0.19, 0.19)                             # taken from PAGE-ICE v6.20 default
                # sltemp_SLtemprise = TriangularDist(1.5, 1.5, 1.5)                            # median sensitivity to GMST changes
                # sla_SLbaselinerise = TriangularDist(1., 1., 1.)                           # asymptote for pre-industrial
                # sltau_SLresponsetime = TriangularDist(386.0, 386.0, 386.0)                      # mean taken from the default; fat-tailed distribution of time constant T_sl, sea level response time, from mode=362, mean = 386
                #
                # # GDP
                # isat0_initialimpactfxnsaturation = TriangularDist(20, 20, 20)
                #
                # # MarketDamages
                # tcal_CalibrationTemp = TriangularDist(3., 3., 3.)
                # iben_MarketInitialBenefit = TriangularDist(.1, .1, .1)
                # W_MarketImpactsatCalibrationTemp = TriangularDist(.5, .5, .5)
                # pow_MarketImpactExponent = TriangularDist(2, 2, 2)
                # ipow_MarketIncomeFxnExponent = TriangularDist(-.1, -.1, -.1)
                #
                # # MarketDamagesBurke
                # impf_coeff_lin = TriangularDist(-0.00829990966469437, -0.00829990966469437, -0.00829990966469437)
                # impf_coeff_quadr = TriangularDist(-0.000500003403703578, -0.000500003403703578, -0.000500003403703578)
                #
                # rtl_abs_0_realizedabstemperature["EU"] = TriangularDist(10.1222005657257, 10.1222005657257, 10.1222005657257)
                # rtl_abs_0_realizedabstemperature["USA"] = TriangularDist(13.4286201975287, 13.4286201975287, 13.4286201975287)
                # rtl_abs_0_realizedabstemperature["OECD"] = TriangularDist(12.0633541623396, 12.0633541623396, 12.0633541623396)
                # rtl_abs_0_realizedabstemperature["USSR"] = TriangularDist(7.11321295839747, 7.11321295839747, 7.11321295839747)
                # rtl_abs_0_realizedabstemperature["China"] = TriangularDist(15.0129576117242, 15.0129576117242, 15.0129576117242)
                # rtl_abs_0_realizedabstemperature["SEAsia"] = TriangularDist(24.9499789823412, 24.9499789823412, 24.9499789823412)
                # rtl_abs_0_realizedabstemperature["Africa"] = TriangularDist(21.892251349428, 21.892251349428, 21.892251349428)
                # rtl_abs_0_realizedabstemperature["LatAmerica"] = TriangularDist(21.1203990194037, 21.1203990194037, 21.1203990194037)
                #
                # # NonMarketDamages
                # tcal_CalibrationTemp = TriangularDist(3., 3., 3.)
                # iben_NonMarketInitialBenefit = TriangularDist(.05, .05, .05)
                # w_NonImpactsatCalibrationTemp = TriangularDist(.5, .5, .5)
                # pow_NonMarketExponent = TriangularDist(2, 2, 2)
                # ipow_NonMarketIncomeFxnExponent = TriangularDist(0, 0, 0)
                #
                # # SLRDamages
                # scal_calibrationSLR = TriangularDist(.5, .5, .5)
                # #iben_SLRInitialBenefit = TriangularDist(0, 0, 0) # only usable if lb <> ub
                # W_SatCalibrationSLR = TriangularDist(1, 1, 1)
                # pow_SLRImpactFxnExponent = TriangularDist(.7, .7, .7)
                # ipow_SLRIncomeFxnExponent = TriangularDist(-.3, -.3, -.3)
                #
                # # Discountinuity
                # rand_discontinuity = Uniform(.5, .5)
                # tdis_tolerabilitydisc = TriangularDist(1.5, 1.5, 1.5)
                # pdis_probability = TriangularDist(20, 20, 20)
                # wdis_gdplostdisc = TriangularDist(3, 3, 3)
                # ipow_incomeexponent = TriangularDist(-.1, -.1, -.1)
                # distau_discontinuityexponent = TriangularDist(20, 20, 20)
                #
                # # EquityWeighting
                # civvalue_civilizationvalue = TriangularDist(5e10, 5e10, 5e10)
                # ptp_timepreference = TriangularDist(1,1,1)
                # emuc_utilityconvexity = TriangularDist(1,1,1)
                #
                # # AbatementCosts
                # AbatementCostParametersCO2_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-22,-22,-22)
                # AbatementCostParametersCH4_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-30,-30,-30)
                # AbatementCostParametersN2O_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-7.0,-7.0,-7.0)
                # AbatementCostParametersLin_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(0,0,0)
                #
                # AbatementCostParametersCO2_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(20,20,20)
                # AbatementCostParametersCH4_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(10,10,10)
                # AbatementCostParametersN2O_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(10,10,10)
                # AbatementCostParametersLin_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(10,10,10)
                #
                # AbatementCostParametersCO2_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-200,-200,-200)
                # AbatementCostParametersCH4_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-4000,-4000,-4000)
                # AbatementCostParametersN2O_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-7000,-7000,-7000)
                # AbatementCostParametersLin_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-200,-200,-200)
                #
                # AbatementCostParametersCO2_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(70,70,70)
                # AbatementCostParametersCH4_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(50,50,50)
                # AbatementCostParametersN2O_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(50,50,50)
                # AbatementCostParametersLin_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(70,70,70)
                #
                # AbatementCostParametersCO2_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(400,400,400)
                # AbatementCostParametersCH4_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(6000,6000,6000)
                # AbatementCostParametersN2O_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(20000,20000,20000)
                # AbatementCostParametersLin_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(300,300,300)
                #
                # AbatementCostParametersCO2_ies_InitialExperienceStockofCutbacks = TriangularDist(150000,150000,150000)
                # AbatementCostParametersCH4_ies_InitialExperienceStockofCutbacks = TriangularDist(2000,2000,2000)
                # AbatementCostParametersN2O_ies_InitialExperienceStockofCutbacks = TriangularDist(50,50,50)
                # AbatementCostParametersLin_ies_InitialExperienceStockofCutbacks = TriangularDist(2000,2000,2000)
                #
                # #the following variables need to be set, but set the same in all 4 abatement cost components
                # #note that for these regional variables, the first region is the focus region (EU), which is set in the preceding code, and so is always one for these variables
                #
                # emitf_uncertaintyinBAUemissfactor["USA"] = TriangularDist(1.0,1.0,1.0)
                # emitf_uncertaintyinBAUemissfactor["OECD"] = TriangularDist(1.0,1.0,1.0)
                # emitf_uncertaintyinBAUemissfactor["USSR"] = TriangularDist(1.0,1.0,1.0)
                # emitf_uncertaintyinBAUemissfactor["China"] = TriangularDist(1.0,1.0,1.0)
                # emitf_uncertaintyinBAUemissfactor["SEAsia"] = TriangularDist(1.0,1.0,1.0)
                # emitf_uncertaintyinBAUemissfactor["Africa"] = TriangularDist(1.0,1.0,1.0)
                # emitf_uncertaintyinBAUemissfactor["LatAmerica"] = TriangularDist(1.0,1.0,1.0)
                #
                # q0f_negativecostpercentagefactor["USA"] = TriangularDist(1.0,1.0,1.0)
                # q0f_negativecostpercentagefactor["OECD"] = TriangularDist(1.0,1.0,1.0)
                # q0f_negativecostpercentagefactor["USSR"] = TriangularDist(0.7,0.7,0.7)
                # q0f_negativecostpercentagefactor["China"] = TriangularDist(0.7,0.7,0.7)
                # q0f_negativecostpercentagefactor["SEAsia"] = TriangularDist(0.7,0.7,0.7)
                # q0f_negativecostpercentagefactor["Africa"] = TriangularDist(0.7,0.7,0.7)
                # q0f_negativecostpercentagefactor["LatAmerica"] = TriangularDist(0.7,0.7,0.7)
                #
                # cmaxf_maxcostfactor["USA"] = TriangularDist(1.0,1.0,1.0)
                # cmaxf_maxcostfactor["OECD"] = TriangularDist(1.2,1.2,1.2)
                # cmaxf_maxcostfactor["USSR"] = TriangularDist(0.7,0.7,0.7)
                # cmaxf_maxcostfactor["China"] = TriangularDist(1.0,1.0,1.0)
                # cmaxf_maxcostfactor["SEAsia"] = TriangularDist(1.2,1.2,1.2)
                # cmaxf_maxcostfactor["Africa"] = TriangularDist(1.2,1.2,1.2)
                # cmaxf_maxcostfactor["LatAmerica"] = TriangularDist(0.7,0.7,0.7)
                #
                # q0propmult_cutbacksatnegativecostinfinalyear = TriangularDist(0.7,0.7,0.7)
                # qmax_minus_q0propmult_maxcutbacksatpositivecostinfinalyear = TriangularDist(1.3,1.3,1.3)
                # c0mult_mostnegativecostinfinalyear = TriangularDist(0.8,0.8,0.8)
                # curve_below_curvatureofMACcurvebelowzerocost = TriangularDist(0.45,0.45,0.45)
                # curve_above_curvatureofMACcurveabovezerocost = TriangularDist(0.4,0.4,0.4)
                # cross_experiencecrossoverratio = TriangularDist(0.2,0.2,0.2)
                # learn_learningrate = TriangularDist(0.2,0.2,0.2)
                #
                # # AdaptationCosts
                # AdaptiveCostsSeaLevel_cp_costplateau_eu = TriangularDist(0.02, 0.02, 0.02)
                # AdaptiveCostsSeaLevel_ci_costimpact_eu = TriangularDist(0.001, 0.001, 0.001)
                # AdaptiveCostsEconomic_cp_costplateau_eu = TriangularDist(0.01, 0.01, 0.01)
                # AdaptiveCostsEconomic_ci_costimpact_eu = TriangularDist(0.003, 0.003, 0.003)
                # AdaptiveCostsNonEconomic_cp_costplateau_eu = TriangularDist(0.02, 0.02, 0.02)
                # AdaptiveCostsNonEconomic_ci_costimpact_eu = TriangularDist(0.005, 0.005, 0.005)
                #
                # cf_costregional["USA"] = TriangularDist(0.8, 0.8, 0.8)
                # cf_costregional["OECD"] = TriangularDist(0.8, 0.8, 0.8)
                # cf_costregional["USSR"] = TriangularDist(0.4, 0.4, 0.4)
                # cf_costregional["China"] = TriangularDist(0.8, 0.8, 0.8)
                # cf_costregional["SEAsia"] = TriangularDist(0.8, 0.8, 0.8)
                # cf_costregional["Africa"] = TriangularDist(0.6, 0.6, 0.6)
                # cf_costregional["LatAmerica"] = TriangularDist(0.6, 0.6, 0.6)

                ############################################################################
                # Indicate which parameters to save for each model run
                ############################################################################
                    # add some annual variables here:
                save(EquityWeighting.td_totaldiscountedimpacts,
                     EquityWeighting.tpc_totalaggregatedcosts,
                     EquityWeighting.tac_totaladaptationcosts,
                     EquityWeighting.te_totaleffect,
                     CO2Cycle.c_CO2concentration,
                     TotalForcing.ft_totalforcing,
                     ClimateTemperature.rt_g_globaltemperature,
                     SeaLevelRise.s_sealevel,
                     SLRDamages.rgdp_per_cap_SLRRemainGDP,
                     MarketDamagesBurke.rgdp_per_cap_MarketRemainGDP,
                     NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP,
                     MarketDamagesBurke.isat_per_cap_ImpactperCapinclSaturationandAdaptation_sum, # added: burke damages (timestep)
                     MarketDamagesBurke.isat_per_cap_ImpactperCapinclSaturationandAdaptation_ann_sum, # added: burke damages (annual)
                     Discontinuity.isat_per_cap_DiscImpactperCapinclSaturation_sum, # added: discontinuity damages (timestep)
                     Discontinuity.isat_per_cap_DiscImpactperCapinclSaturation_ann_sum, # added: discontinuity damages (annual)
                     #Discontinuity.rgdp_per_cap_NonMarketRemainGDP, # this is a parameter in Discontinuity and therefore does not change from NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP.
                     Discontinuity.rcons_per_cap_DiscRemainConsumption,
                     Discontinuity.rcons_per_cap_DiscRemainConsumption_sum, # NEW
                     # # add these new variables to the RV bit below as well.
                     EquityWeighting.td_totaldiscountedimpacts_ann, # NEW
                     Discontinuity.rcons_per_cap_DiscRemainConsumption_ann_sum, # NEW
                     Discontinuity.occurdis_occurrencedummy_sum, # NEW
                     Discontinuity.occurdis_occurrencedummy_ann_sum, # NEW
                     MarketDamagesBurke.rgdp_per_cap_MarketRemainGDP_ann_sum, # NEW
                     MarketDamagesBurke.rgdp_per_cap_MarketRemainGDP_sum, # NEW
                     NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP_ann_sum, # NEW
                     NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP_sum, # NEW
                     NonMarketDamages.isat_per_cap_ImpactperCapinclSaturationandAdaptation_ann_sum_nonmarket) # NEW

                 end
        elseif use_no_varMC
            mcs = @defsim begin
                gvarsd_globalvariabilitystandarddeviation = Normal(0.11294, 0.) # ranges from https://github.com/jkikstra/climvar
                rvarsd_regionalvariabilitystandarddeviation["EU"] = Normal(0.35251, 0.)
                rvarsd_regionalvariabilitystandarddeviation["USA"] = Normal(0.39171, 0.)
                rvarsd_regionalvariabilitystandarddeviation["OECD"] = Normal(0.35272, 0.)
                rvarsd_regionalvariabilitystandarddeviation["USSR"] = Normal(0.39679, 0.)
                rvarsd_regionalvariabilitystandarddeviation["China"] = Normal(0.26032, 0.)
                rvarsd_regionalvariabilitystandarddeviation["SEAsia"] = Normal(0.20009, 0.)
                rvarsd_regionalvariabilitystandarddeviation["Africa"] = Normal(0.18202, 0.)
                rvarsd_regionalvariabilitystandarddeviation["LatAmerica"] = Normal(0.15618, 0.)
                #The folllowing RVs are in more than one component.  For clarity they are
                #set here as opposed to below within the blocks of RVs separated by component
                #so that they are not set more than once.
                save_savingsrate = TriangularDist(10, 20, 15) # components: MarketDamages, MarketDamagesBurke, NonMarketDamages. GDP, SLRDamages

                wincf_weightsfactor_sea["USA"] = TriangularDist(.6, 1, .8) # components: SLRDamages, Discountinuity (weights for market and nonmarket are non-stochastic and uniformly 1)
                wincf_weightsfactor_sea["OECD"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor_sea["USSR"] = TriangularDist(.2, .6, .4)
                wincf_weightsfactor_sea["China"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor_sea["SEAsia"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor_sea["Africa"] = TriangularDist(.4, .8, .6)
                wincf_weightsfactor_sea["LatAmerica"] = TriangularDist(.4, .8, .6)

                automult_autonomouschange = TriangularDist(0.5, 0.8, 0.65)  #components: AdaptationCosts, AbatementCosts

                #The following RVs are divided into blocks by component

                # CO2cycle
                air_CO2fractioninatm = TriangularDist(57, 67, 62)
                res_CO2atmlifetime = TriangularDist(50, 100, 70)
                #ccf_CO2feedback = TriangularDist(0, 0, 0) # only usable if lb <> ub
                ccfmax_maxCO2feedback = TriangularDist(10, 30, 20)
                stay_fractionCO2emissionsinatm = TriangularDist(0.25,0.35,0.3)
                ce_0_basecumCO2emissions=TriangularDist(1830000, 2240000, 2040000)
                a1_percentco2oceanlong=TriangularDist( 4.3,	41.6, 23.0)
                a2_percentco2oceanshort=TriangularDist(23.1, 30.1, 26.6)
                a3_percentco2land=TriangularDist(11.4, 42.5, 27.0)
                t1_timeco2oceanlong=TriangularDist(248.9, 376.2, 312.5)
                t2_timeco2oceanshort=TriangularDist(25.9, 43.9, 34.9)
                t3_timeco2land=TriangularDist(2.8, 5.7, 4.3)
                rt_g0_baseglobaltemp=TriangularDist(0.903, 0.989, 0.946)

                # SiBCASA Permafrost
                perm_sib_af = TriangularDist(1.42609149897258, 2.32504747848815, 1.87556948873036)
                perm_sib_sens_c_co2 = TriangularDist(28191.1555428869, 35688.3253432574, 31939.7404430722)
                perm_sib_lag_c_co2 = TriangularDist(35.4926669856915, 87.8949041341782, 61.6937855599349)
                perm_sib_pow_c_co2 = TriangularDist(0.107020247715729, 0.410961185142816, 0.258990716429273)
                perm_sib_sens_c_ch4 = TriangularDist(1240.3553299183, 3348.11995329232, 2294.23764160531)
                perm_sib_lag_c_ch4 = TriangularDist(75.1943160023131, 337.382510123922, 206.288413063117)
                perm_sib_pow_c_ch4 = TriangularDist(-0.108779283732708, 0.610889007954489, 0.25105486211089)

                # SulphateForcing
                d_sulphateforcingbase = TriangularDist(-0.8, -0.2, -0.4)
                ind_slopeSEforcing_indirect = TriangularDist(-0.8, 0, -0.4)

                # ClimateTemperature
                frt_warminghalflife = TriangularDist(10, 55, 20)        # from PAGE-ICE v6.2 documentation
                tcr_transientresponse = TriangularDist(0.8, 2.7, 1.8)   # from PAGE-ICE v6.2 documentation
                alb_emulator_rand = TriangularDist(-1., 1., 0.)
                ampf_amplification["EU"] = TriangularDist(1.05, 1.53, 1.23)
                ampf_amplification["USA"] = TriangularDist(1.16, 1.54, 1.32)
                ampf_amplification["OECD"] = TriangularDist(1.14, 1.31, 1.21)
                ampf_amplification["USSR"] = TriangularDist(1.41, 1.9, 1.64)
                ampf_amplification["China"] = TriangularDist(1, 1.3, 1.21)
                ampf_amplification["SEAsia"] = TriangularDist(0.84, 1.15, 1.04)
                ampf_amplification["Africa"] = TriangularDist(0.99, 1.42, 1.22)
                ampf_amplification["LatAmerica"] = TriangularDist(0.9, 1.18, 1.04)

                # SeaLevelRise
                s0_initialSL = TriangularDist(0.17, 0.21, 0.19)                             # taken from PAGE-ICE v6.20 default
                sltemp_SLtemprise = TriangularDist(0.7, 3., 1.5)                            # median sensitivity to GMST changes
                sla_SLbaselinerise = TriangularDist(0.5, 1.5, 1.)                           # asymptote for pre-industrial
                sltau_SLresponsetime = Gamma(16.0833333333333333, 24.)                      # fat-tailed distribution of time constant T_sl, sea level response time, from mode=362, mean = 386

                # GDP
                isat0_initialimpactfxnsaturation = TriangularDist(15, 25, 20)

                # MarketDamages
                tcal_CalibrationTemp = TriangularDist(2.5, 3.5, 3.)
                iben_MarketInitialBenefit = TriangularDist(0, .3, .1)
                W_MarketImpactsatCalibrationTemp = TriangularDist(.2, .8, .5)
                pow_MarketImpactExponent = TriangularDist(1.5, 3, 2)
                ipow_MarketIncomeFxnExponent = TriangularDist(-.3, 0, -.1)

                # MarketDamagesBurke
                impf_coeff_lin = TriangularDist(-0.0139791885347898, -0.0026206307945989, -0.00829990966469437)
                impf_coeff_quadr = TriangularDist(-0.000599999506482576, -0.000400007300924579, -0.000500003403703578)

                rtl_abs_0_realizedabstemperature["EU"] = TriangularDist(6.76231496767033, 13.482086163781, 10.1222005657257)
                rtl_abs_0_realizedabstemperature["USA"] = TriangularDist(9.54210085883826, 17.3151395362191, 13.4286201975287)
                rtl_abs_0_realizedabstemperature["OECD"] = TriangularDist(9.07596053028087, 15.0507477943984, 12.0633541623396)
                rtl_abs_0_realizedabstemperature["USSR"] = TriangularDist(3.01320548016903, 11.2132204366259, 7.11321295839747)
                rtl_abs_0_realizedabstemperature["China"] = TriangularDist(12.2330402806912, 17.7928749427573, 15.0129576117242)
                rtl_abs_0_realizedabstemperature["SEAsia"] = TriangularDist(23.3863348263352, 26.5136231383473, 24.9499789823412)
                rtl_abs_0_realizedabstemperature["Africa"] = TriangularDist(20.1866940491107, 23.5978086497453, 21.892251349428)
                rtl_abs_0_realizedabstemperature["LatAmerica"] = TriangularDist(19.4846849750102, 22.7561130637973, 21.1203990194037)

                # NonMarketDamages
                tcal_CalibrationTemp = TriangularDist(2.5, 3.5, 3.)
                iben_NonMarketInitialBenefit = TriangularDist(0, .2, .05)
                w_NonImpactsatCalibrationTemp = TriangularDist(.1, 1, .5)
                pow_NonMarketExponent = TriangularDist(1.5, 3, 2)
                ipow_NonMarketIncomeFxnExponent = TriangularDist(-.2, .2, 0)

                # SLRDamages
                scal_calibrationSLR = TriangularDist(0.45, 0.55, .5)
                #iben_SLRInitialBenefit = TriangularDist(0, 0, 0) # only usable if lb <> ub
                W_SatCalibrationSLR = TriangularDist(.5, 1.5, 1)
                pow_SLRImpactFxnExponent = TriangularDist(.5, 1, .7)
                ipow_SLRIncomeFxnExponent = TriangularDist(-.4, -.2, -.3)

                # Discountinuity
                rand_discontinuity = Uniform(0, 1)
                tdis_tolerabilitydisc = TriangularDist(1, 2, 1.5)
                pdis_probability = TriangularDist(10, 30, 20)
                wdis_gdplostdisc = TriangularDist(1, 5, 3)
                ipow_incomeexponent = TriangularDist(-.3, 0, -.1)
                distau_discontinuityexponent = TriangularDist(10, 30, 20)

                # EquityWeighting
                civvalue_civilizationvalue = TriangularDist(1e10, 1e11, 5e10)
                ptp_timepreference = TriangularDist(0.1,2,1)
                emuc_utilityconvexity = TriangularDist(0.5,2,1)

                # AbatementCosts
                AbatementCostParametersCO2_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,6.0,-22)
                AbatementCostParametersCH4_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-67,6.0,-30)
                AbatementCostParametersN2O_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-20,6.0,-7.0)
                AbatementCostParametersLin_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,50,0)

                AbatementCostParametersCO2_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,40,20)
                AbatementCostParametersCH4_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)
                AbatementCostParametersN2O_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)
                AbatementCostParametersLin_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)

                AbatementCostParametersCO2_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400,-100,-200)
                AbatementCostParametersCH4_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-8000,-1000,-4000)
                AbatementCostParametersN2O_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-15000,0,-7000)
                AbatementCostParametersLin_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400,-100,-200)

                AbatementCostParametersCO2_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60,80,70)
                AbatementCostParametersCH4_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35,70,50)
                AbatementCostParametersN2O_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35,70,50)
                AbatementCostParametersLin_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60,80,70)

                AbatementCostParametersCO2_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100,700,400)
                AbatementCostParametersCH4_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(3000,10000,6000)
                AbatementCostParametersN2O_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(2000,60000,20000)
                AbatementCostParametersLin_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100,600,300)

                AbatementCostParametersCO2_ies_InitialExperienceStockofCutbacks = TriangularDist(100000,200000,150000)
                AbatementCostParametersCH4_ies_InitialExperienceStockofCutbacks = TriangularDist(1500,2500,2000)
                AbatementCostParametersN2O_ies_InitialExperienceStockofCutbacks = TriangularDist(30,80,50)
                AbatementCostParametersLin_ies_InitialExperienceStockofCutbacks = TriangularDist(1500,2500,2000)

                #the following variables need to be set, but set the same in all 4 abatement cost components
                #note that for these regional variables, the first region is the focus region (EU), which is set in the preceding code, and so is always one for these variables

                emitf_uncertaintyinBAUemissfactor["USA"] = TriangularDist(0.8,1.2,1.0)
                emitf_uncertaintyinBAUemissfactor["OECD"] = TriangularDist(0.8,1.2,1.0)
                emitf_uncertaintyinBAUemissfactor["USSR"] = TriangularDist(0.65,1.35,1.0)
                emitf_uncertaintyinBAUemissfactor["China"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["SEAsia"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["Africa"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["LatAmerica"] = TriangularDist(0.5,1.5,1.0)

                q0f_negativecostpercentagefactor["USA"] = TriangularDist(0.75,1.5,1.0)
                q0f_negativecostpercentagefactor["OECD"] = TriangularDist(0.75,1.25,1.0)
                q0f_negativecostpercentagefactor["USSR"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["China"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["SEAsia"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["Africa"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["LatAmerica"] = TriangularDist(0.4,1.0,0.7)

                cmaxf_maxcostfactor["USA"] = TriangularDist(0.8,1.2,1.0)
                cmaxf_maxcostfactor["OECD"] = TriangularDist(1.0,1.5,1.2)
                cmaxf_maxcostfactor["USSR"] = TriangularDist(0.4,1.0,0.7)
                cmaxf_maxcostfactor["China"] = TriangularDist(0.8,1.2,1.0)
                cmaxf_maxcostfactor["SEAsia"] = TriangularDist(1,1.5,1.2)
                cmaxf_maxcostfactor["Africa"] = TriangularDist(1,1.5,1.2)
                cmaxf_maxcostfactor["LatAmerica"] = TriangularDist(0.4,1.0,0.7)

                q0propmult_cutbacksatnegativecostinfinalyear = TriangularDist(0.3,1.2,0.7)
                qmax_minus_q0propmult_maxcutbacksatpositivecostinfinalyear = TriangularDist(1,1.5,1.3)
                c0mult_mostnegativecostinfinalyear = TriangularDist(0.5,1.2,0.8)
                curve_below_curvatureofMACcurvebelowzerocost = TriangularDist(0.25,0.8,0.45)
                curve_above_curvatureofMACcurveabovezerocost = TriangularDist(0.1,0.7,0.4)
                cross_experiencecrossoverratio = TriangularDist(0.1,0.3,0.2)
                learn_learningrate = TriangularDist(0.05,0.35,0.2)

                # AdaptationCosts
                AdaptiveCostsSeaLevel_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
                AdaptiveCostsSeaLevel_ci_costimpact_eu = TriangularDist(0.0005, 0.002, 0.001)
                AdaptiveCostsEconomic_cp_costplateau_eu = TriangularDist(0.005, 0.02, 0.01)
                AdaptiveCostsEconomic_ci_costimpact_eu = TriangularDist(0.001, 0.008, 0.003)
                AdaptiveCostsNonEconomic_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
                AdaptiveCostsNonEconomic_ci_costimpact_eu = TriangularDist(0.002, 0.01, 0.005)

                cf_costregional["USA"] = TriangularDist(0.6, 1, 0.8)
                cf_costregional["OECD"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["USSR"] = TriangularDist(0.2, 0.6, 0.4)
                cf_costregional["China"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["SEAsia"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["Africa"] = TriangularDist(0.4, 0.8, 0.6)
                cf_costregional["LatAmerica"] = TriangularDist(0.4, 0.8, 0.6)

                ############################################################################
                # Indicate which parameters to save for each model run
                ############################################################################
                    # add some annual variables here:
                save(EquityWeighting.td_totaldiscountedimpacts,
                     EquityWeighting.tpc_totalaggregatedcosts,
                     EquityWeighting.tac_totaladaptationcosts,
                     EquityWeighting.te_totaleffect,
                     CO2Cycle.c_CO2concentration,
                     TotalForcing.ft_totalforcing,
                     ClimateTemperature.rt_g_globaltemperature,
                     SeaLevelRise.s_sealevel,
                     SLRDamages.rgdp_per_cap_SLRRemainGDP,
                     MarketDamagesBurke.rgdp_per_cap_MarketRemainGDP,
                     NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP,
                     MarketDamagesBurke.isat_per_cap_ImpactperCapinclSaturationandAdaptation_sum, # added: burke damages (timestep)
                     MarketDamagesBurke.isat_per_cap_ImpactperCapinclSaturationandAdaptation_ann_sum, # added: burke damages (annual)
                     Discontinuity.isat_per_cap_DiscImpactperCapinclSaturation_sum, # added: discontinuity damages (timestep)
                     Discontinuity.isat_per_cap_DiscImpactperCapinclSaturation_ann_sum, # added: discontinuity damages (annual)
                     #Discontinuity.rgdp_per_cap_NonMarketRemainGDP, # this is a parameter in Discontinuity and therefore does not change from NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP.
                     Discontinuity.rcons_per_cap_DiscRemainConsumption,
                     Discontinuity.rcons_per_cap_DiscRemainConsumption_sum, # NEW
                     # # add these new variables to the RV bit below as well.
                     EquityWeighting.td_totaldiscountedimpacts_ann, # NEW
                     Discontinuity.rcons_per_cap_DiscRemainConsumption_ann_sum, # NEW
                     Discontinuity.occurdis_occurrencedummy_sum, # NEW
                     Discontinuity.occurdis_occurrencedummy_ann_sum, # NEW
                     MarketDamagesBurke.rgdp_per_cap_MarketRemainGDP_ann_sum, # NEW
                     MarketDamagesBurke.rgdp_per_cap_MarketRemainGDP_sum, # NEW
                     NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP_ann_sum, # NEW
                     NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP_sum, # NEW
                     NonMarketDamages.isat_per_cap_ImpactperCapinclSaturationandAdaptation_ann_sum_nonmarket) # NEW
            end
        else
            mcs = @defsim begin
                gvarsd_globalvariabilitystandarddeviation = Normal(0.11294, 0.021644) # ranges from https://github.com/jkikstra/climvar
                rvarsd_regionalvariabilitystandarddeviation["EU"] = Normal(0.35251, 0.019147)
                rvarsd_regionalvariabilitystandarddeviation["USA"] = Normal(0.39171, 0.030440)
                rvarsd_regionalvariabilitystandarddeviation["OECD"] = Normal(0.35272, 0.042962)
                rvarsd_regionalvariabilitystandarddeviation["USSR"] = Normal(0.39679, 0.030228)
                rvarsd_regionalvariabilitystandarddeviation["China"] = Normal(0.26032, 0.059394)
                rvarsd_regionalvariabilitystandarddeviation["SEAsia"] = Normal(0.20009, 0.039896)
                rvarsd_regionalvariabilitystandarddeviation["Africa"] = Normal(0.18202, 0.050171)
                rvarsd_regionalvariabilitystandarddeviation["LatAmerica"] = Normal(0.15618, 0.054156)

                save_savingsrate = TriangularDist(10, 20, 15) # components: MarketDamages, MarketDamagesBurke, NonMarketDamages. GDP, SLRDamages

                wincf_weightsfactor_sea["USA"] = TriangularDist(.6, 1, .8) # components: SLRDamages, Discountinuity (weights for market and nonmarket are non-stochastic and uniformly 1)
                wincf_weightsfactor_sea["OECD"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor_sea["USSR"] = TriangularDist(.2, .6, .4)
                wincf_weightsfactor_sea["China"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor_sea["SEAsia"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor_sea["Africa"] = TriangularDist(.4, .8, .6)
                wincf_weightsfactor_sea["LatAmerica"] = TriangularDist(.4, .8, .6)

                automult_autonomouschange = TriangularDist(0.5, 0.8, 0.65)  #components: AdaptationCosts, AbatementCosts

                #The following RVs are divided into blocks by component

                # CO2cycle
                air_CO2fractioninatm = TriangularDist(57, 67, 62)
                res_CO2atmlifetime = TriangularDist(50, 100, 70)
                #ccf_CO2feedback = TriangularDist(0, 0, 0) # only usable if lb <> ub
                ccfmax_maxCO2feedback = TriangularDist(10, 30, 20)
                stay_fractionCO2emissionsinatm = TriangularDist(0.25,0.35,0.3)
                ce_0_basecumCO2emissions=TriangularDist(1830000, 2240000, 2040000)
                a1_percentco2oceanlong=TriangularDist( 4.3,	41.6, 23.0)
                a2_percentco2oceanshort=TriangularDist(23.1, 30.1, 26.6)
                a3_percentco2land=TriangularDist(11.4, 42.5, 27.0)
                t1_timeco2oceanlong=TriangularDist(248.9, 376.2, 312.5)
                t2_timeco2oceanshort=TriangularDist(25.9, 43.9, 34.9)
                t3_timeco2land=TriangularDist(2.8, 5.7, 4.3)
                rt_g0_baseglobaltemp=TriangularDist(0.903, 0.989, 0.946)

                # SiBCASA Permafrost
                perm_sib_af = TriangularDist(1.42609149897258, 2.32504747848815, 1.87556948873036)
                perm_sib_sens_c_co2 = TriangularDist(28191.1555428869, 35688.3253432574, 31939.7404430722)
                perm_sib_lag_c_co2 = TriangularDist(35.4926669856915, 87.8949041341782, 61.6937855599349)
                perm_sib_pow_c_co2 = TriangularDist(0.107020247715729, 0.410961185142816, 0.258990716429273)
                perm_sib_sens_c_ch4 = TriangularDist(1240.3553299183, 3348.11995329232, 2294.23764160531)
                perm_sib_lag_c_ch4 = TriangularDist(75.1943160023131, 337.382510123922, 206.288413063117)
                perm_sib_pow_c_ch4 = TriangularDist(-0.108779283732708, 0.610889007954489, 0.25105486211089)

                # SulphateForcing
                d_sulphateforcingbase = TriangularDist(-0.8, -0.2, -0.4)
                ind_slopeSEforcing_indirect = TriangularDist(-0.8, 0, -0.4)

                # ClimateTemperature
                frt_warminghalflife = TriangularDist(10, 55, 20)        # from PAGE-ICE v6.2 documentation
                tcr_transientresponse = TriangularDist(0.8, 2.7, 1.8)   # from PAGE-ICE v6.2 documentation
                alb_emulator_rand = TriangularDist(-1., 1., 0.)
                ampf_amplification["EU"] = TriangularDist(1.05, 1.53, 1.23)
                ampf_amplification["USA"] = TriangularDist(1.16, 1.54, 1.32)
                ampf_amplification["OECD"] = TriangularDist(1.14, 1.31, 1.21)
                ampf_amplification["USSR"] = TriangularDist(1.41, 1.9, 1.64)
                ampf_amplification["China"] = TriangularDist(1, 1.3, 1.21)
                ampf_amplification["SEAsia"] = TriangularDist(0.84, 1.15, 1.04)
                ampf_amplification["Africa"] = TriangularDist(0.99, 1.42, 1.22)
                ampf_amplification["LatAmerica"] = TriangularDist(0.9, 1.18, 1.04)

                # SeaLevelRise
                s0_initialSL = TriangularDist(0.17, 0.21, 0.19)                             # taken from PAGE-ICE v6.20 default
                sltemp_SLtemprise = TriangularDist(0.7, 3., 1.5)                            # median sensitivity to GMST changes
                sla_SLbaselinerise = TriangularDist(0.5, 1.5, 1.)                           # asymptote for pre-industrial
                sltau_SLresponsetime = Gamma(16.0833333333333333, 24.)                      # fat-tailed distribution of time constant T_sl, sea level response time, from mode=362, mean = 386

                # GDP
                isat0_initialimpactfxnsaturation = TriangularDist(15, 25, 20)

                # MarketDamages
                tcal_CalibrationTemp = TriangularDist(2.5, 3.5, 3.)
                iben_MarketInitialBenefit = TriangularDist(0, .3, .1)
                W_MarketImpactsatCalibrationTemp = TriangularDist(.2, .8, .5)
                pow_MarketImpactExponent = TriangularDist(1.5, 3, 2)
                ipow_MarketIncomeFxnExponent = TriangularDist(-.3, 0, -.1)

                # MarketDamagesBurke
                impf_coeff_lin = TriangularDist(-0.0139791885347898, -0.0026206307945989, -0.00829990966469437)
                impf_coeff_quadr = TriangularDist(-0.000599999506482576, -0.000400007300924579, -0.000500003403703578)

                rtl_abs_0_realizedabstemperature["EU"] = TriangularDist(6.76231496767033, 13.482086163781, 10.1222005657257)
                rtl_abs_0_realizedabstemperature["USA"] = TriangularDist(9.54210085883826, 17.3151395362191, 13.4286201975287)
                rtl_abs_0_realizedabstemperature["OECD"] = TriangularDist(9.07596053028087, 15.0507477943984, 12.0633541623396)
                rtl_abs_0_realizedabstemperature["USSR"] = TriangularDist(3.01320548016903, 11.2132204366259, 7.11321295839747)
                rtl_abs_0_realizedabstemperature["China"] = TriangularDist(12.2330402806912, 17.7928749427573, 15.0129576117242)
                rtl_abs_0_realizedabstemperature["SEAsia"] = TriangularDist(23.3863348263352, 26.5136231383473, 24.9499789823412)
                rtl_abs_0_realizedabstemperature["Africa"] = TriangularDist(20.1866940491107, 23.5978086497453, 21.892251349428)
                rtl_abs_0_realizedabstemperature["LatAmerica"] = TriangularDist(19.4846849750102, 22.7561130637973, 21.1203990194037)

                # NonMarketDamages
                tcal_CalibrationTemp = TriangularDist(2.5, 3.5, 3.)
                iben_NonMarketInitialBenefit = TriangularDist(0, .2, .05)
                w_NonImpactsatCalibrationTemp = TriangularDist(.1, 1, .5)
                pow_NonMarketExponent = TriangularDist(1.5, 3, 2)
                ipow_NonMarketIncomeFxnExponent = TriangularDist(-.2, .2, 0)

                # SLRDamages
                scal_calibrationSLR = TriangularDist(0.45, 0.55, .5)
                #iben_SLRInitialBenefit = TriangularDist(0, 0, 0) # only usable if lb <> ub
                W_SatCalibrationSLR = TriangularDist(.5, 1.5, 1)
                pow_SLRImpactFxnExponent = TriangularDist(.5, 1, .7)
                ipow_SLRIncomeFxnExponent = TriangularDist(-.4, -.2, -.3)

                # Discountinuity
                rand_discontinuity = Uniform(0, 1)
                tdis_tolerabilitydisc = TriangularDist(1, 2, 1.5)
                pdis_probability = TriangularDist(10, 30, 20)
                wdis_gdplostdisc = TriangularDist(1, 5, 3)
                ipow_incomeexponent = TriangularDist(-.3, 0, -.1)
                distau_discontinuityexponent = TriangularDist(10, 30, 20)

                # EquityWeighting
                civvalue_civilizationvalue = TriangularDist(1e10, 1e11, 5e10)
                ptp_timepreference = TriangularDist(0.1,2,1)
                emuc_utilityconvexity = TriangularDist(0.5,2,1)

                # AbatementCosts
                AbatementCostParametersCO2_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,6.0,-22)
                AbatementCostParametersCH4_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-67,6.0,-30)
                AbatementCostParametersN2O_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-20,6.0,-7.0)
                AbatementCostParametersLin_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,50,0)

                AbatementCostParametersCO2_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,40,20)
                AbatementCostParametersCH4_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)
                AbatementCostParametersN2O_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)
                AbatementCostParametersLin_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)

                AbatementCostParametersCO2_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400,-100,-200)
                AbatementCostParametersCH4_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-8000,-1000,-4000)
                AbatementCostParametersN2O_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-15000,0,-7000)
                AbatementCostParametersLin_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400,-100,-200)

                AbatementCostParametersCO2_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60,80,70)
                AbatementCostParametersCH4_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35,70,50)
                AbatementCostParametersN2O_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35,70,50)
                AbatementCostParametersLin_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60,80,70)

                AbatementCostParametersCO2_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100,700,400)
                AbatementCostParametersCH4_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(3000,10000,6000)
                AbatementCostParametersN2O_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(2000,60000,20000)
                AbatementCostParametersLin_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100,600,300)

                AbatementCostParametersCO2_ies_InitialExperienceStockofCutbacks = TriangularDist(100000,200000,150000)
                AbatementCostParametersCH4_ies_InitialExperienceStockofCutbacks = TriangularDist(1500,2500,2000)
                AbatementCostParametersN2O_ies_InitialExperienceStockofCutbacks = TriangularDist(30,80,50)
                AbatementCostParametersLin_ies_InitialExperienceStockofCutbacks = TriangularDist(1500,2500,2000)

                #the following variables need to be set, but set the same in all 4 abatement cost components
                #note that for these regional variables, the first region is the focus region (EU), which is set in the preceding code, and so is always one for these variables

                emitf_uncertaintyinBAUemissfactor["USA"] = TriangularDist(0.8,1.2,1.0)
                emitf_uncertaintyinBAUemissfactor["OECD"] = TriangularDist(0.8,1.2,1.0)
                emitf_uncertaintyinBAUemissfactor["USSR"] = TriangularDist(0.65,1.35,1.0)
                emitf_uncertaintyinBAUemissfactor["China"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["SEAsia"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["Africa"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["LatAmerica"] = TriangularDist(0.5,1.5,1.0)

                q0f_negativecostpercentagefactor["USA"] = TriangularDist(0.75,1.5,1.0)
                q0f_negativecostpercentagefactor["OECD"] = TriangularDist(0.75,1.25,1.0)
                q0f_negativecostpercentagefactor["USSR"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["China"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["SEAsia"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["Africa"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["LatAmerica"] = TriangularDist(0.4,1.0,0.7)

                cmaxf_maxcostfactor["USA"] = TriangularDist(0.8,1.2,1.0)
                cmaxf_maxcostfactor["OECD"] = TriangularDist(1.0,1.5,1.2)
                cmaxf_maxcostfactor["USSR"] = TriangularDist(0.4,1.0,0.7)
                cmaxf_maxcostfactor["China"] = TriangularDist(0.8,1.2,1.0)
                cmaxf_maxcostfactor["SEAsia"] = TriangularDist(1,1.5,1.2)
                cmaxf_maxcostfactor["Africa"] = TriangularDist(1,1.5,1.2)
                cmaxf_maxcostfactor["LatAmerica"] = TriangularDist(0.4,1.0,0.7)

                q0propmult_cutbacksatnegativecostinfinalyear = TriangularDist(0.3,1.2,0.7)
                qmax_minus_q0propmult_maxcutbacksatpositivecostinfinalyear = TriangularDist(1,1.5,1.3)
                c0mult_mostnegativecostinfinalyear = TriangularDist(0.5,1.2,0.8)
                curve_below_curvatureofMACcurvebelowzerocost = TriangularDist(0.25,0.8,0.45)
                curve_above_curvatureofMACcurveabovezerocost = TriangularDist(0.1,0.7,0.4)
                cross_experiencecrossoverratio = TriangularDist(0.1,0.3,0.2)
                learn_learningrate = TriangularDist(0.05,0.35,0.2)

                # AdaptationCosts
                AdaptiveCostsSeaLevel_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
                AdaptiveCostsSeaLevel_ci_costimpact_eu = TriangularDist(0.0005, 0.002, 0.001)
                AdaptiveCostsEconomic_cp_costplateau_eu = TriangularDist(0.005, 0.02, 0.01)
                AdaptiveCostsEconomic_ci_costimpact_eu = TriangularDist(0.001, 0.008, 0.003)
                AdaptiveCostsNonEconomic_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
                AdaptiveCostsNonEconomic_ci_costimpact_eu = TriangularDist(0.002, 0.01, 0.005)

                cf_costregional["USA"] = TriangularDist(0.6, 1, 0.8)
                cf_costregional["OECD"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["USSR"] = TriangularDist(0.2, 0.6, 0.4)
                cf_costregional["China"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["SEAsia"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["Africa"] = TriangularDist(0.4, 0.8, 0.6)
                cf_costregional["LatAmerica"] = TriangularDist(0.4, 0.8, 0.6)
                ############################################################################
                # Indicate which parameters to save for each model run
                ############################################################################
                    # add some annual variables here:
                save(EquityWeighting.td_totaldiscountedimpacts,
                     EquityWeighting.tpc_totalaggregatedcosts,
                     EquityWeighting.tac_totaladaptationcosts,
                     EquityWeighting.te_totaleffect,
                     CO2Cycle.c_CO2concentration,
                     TotalForcing.ft_totalforcing,
                     ClimateTemperature.rt_g_globaltemperature,
                     SeaLevelRise.s_sealevel,
                     SLRDamages.rgdp_per_cap_SLRRemainGDP,
                     MarketDamagesBurke.rgdp_per_cap_MarketRemainGDP,
                     NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP,
                     MarketDamagesBurke.isat_per_cap_ImpactperCapinclSaturationandAdaptation_sum, # added: burke damages (timestep)
                     MarketDamagesBurke.isat_per_cap_ImpactperCapinclSaturationandAdaptation_ann_sum, # added: burke damages (annual)
                     Discontinuity.isat_per_cap_DiscImpactperCapinclSaturation_sum, # added: discontinuity damages (timestep)
                     Discontinuity.isat_per_cap_DiscImpactperCapinclSaturation_ann_sum, # added: discontinuity damages (annual)
                     #Discontinuity.rgdp_per_cap_NonMarketRemainGDP, # this is a parameter in Discontinuity and therefore does not change from NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP.
                     Discontinuity.rcons_per_cap_DiscRemainConsumption,
                     Discontinuity.rcons_per_cap_DiscRemainConsumption_sum, # NEW
                     # # add these new variables to the RV bit below as well.
                     EquityWeighting.td_totaldiscountedimpacts_ann, # NEW
                     Discontinuity.rcons_per_cap_DiscRemainConsumption_ann_sum, # NEW
                     Discontinuity.occurdis_occurrencedummy_sum, # NEW
                     Discontinuity.occurdis_occurrencedummy_ann_sum, # NEW
                     MarketDamagesBurke.rgdp_per_cap_MarketRemainGDP_ann_sum, # NEW
                     MarketDamagesBurke.rgdp_per_cap_MarketRemainGDP_sum, # NEW
                     NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP_ann_sum, # NEW
                     NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP_sum, # NEW
                     NonMarketDamages.isat_per_cap_ImpactperCapinclSaturationandAdaptation_ann_sum_nonmarket) # NEW
                 end #defsim
        end #if-elif-else
end

#Reformat the RV results into the format used for testing
function reformat_RV_outputs(samplesize::Int; output_path::String = joinpath(@__DIR__, "../output/"))
    #create vectors to hold results of Monte Carlo runs
    td=zeros(samplesize);
    tpc=zeros(samplesize);
    tac=zeros(samplesize);
    te=zeros(samplesize);
    ft=zeros(samplesize);
    rt_g=zeros(samplesize);
    s=zeros(samplesize);
    c_co2concentration=zeros(samplesize);
    rgdppercap_slr=zeros(samplesize);
    rgdppercap_market=zeros(samplesize);
    rgdppercap_nonmarket=zeros(samplesize);
    rconspercap_disc=zeros(samplesize);
    rconspercap_disc_sum=zeros(samplesize); # NEW
    td_ann = zeros(samplesize); # NEW
    rdgppercap_disc_ann_sum = zeros(samplesize); # NEW
    occurdis_occurrencedummy_sum = zeros(samplesize); # NEW
    occurdis_occurrencedummy_ann_sum = zeros(samplesize); # NEW
    rgdp_per_cap_MarketRemainGDP_ann_sum = zeros(samplesize); # NEW
    rgdp_per_cap_MarketRemainGDP_sum = zeros(samplesize); # NEW
    rgdp_per_cap_NonMarketRemainGDP_ann_sum = zeros(samplesize); # NEW
    rgdp_per_cap_NonMarketRemainGDP_sum = zeros(samplesize); # NEW
    isat_per_cap_ImpactperCapinclSaturationandAdaptation_sum = zeros(samplesize); # added: burke damages (timestep)
    isat_per_cap_ImpactperCapinclSaturationandAdaptation_ann_sum = zeros(samplesize); # added: burke damages (annual)
    isat_per_cap_DiscImpactperCapinclSaturation_sum = zeros(samplesize); # added: discontinuity damages (timestep)
    isat_per_cap_DiscImpactperCapinclSaturation_ann_sum = zeros(samplesize); # added: discontinuity damages (annual)
    isat_per_cap_ImpactperCapinclSaturationandAdaptation_ann_sum_nonmarket = zeros(samplesize); # added: nonmarket damages (annual)

    #load raw data
    #no filter
    td      = load_RV("td_totaldiscountedimpacts"; output_path = output_path)
    tpc     = load_RV("tpc_totalaggregatedcosts"; output_path = output_path)
    tac     = load_RV("tac_totaladaptationcosts"; output_path = output_path)
    te      = load_RV("te_totaleffect"; output_path = output_path)
    td_ann  = load_RV("td_totaldiscountedimpacts_ann"; output_path = output_path) # NEW

    occurdis_occurrencedummy_sum = load_RV("occurdis_occurrencedummy_sum"; output_path = output_path) # NEW
    occurdis_occurrencedummy_ann_sum = load_RV("occurdis_occurrencedummy_ann_sum"; output_path = output_path) # NEW
    burke_impact_per_cap_sum = load_RV("isat_per_cap_ImpactperCapinclSaturationandAdaptation_sum"; output_path = output_path) # added: burke damages (timestep)
    burke_impact_per_cap_ann_sum = load_RV("isat_per_cap_ImpactperCapinclSaturationandAdaptation_ann_sum"; output_path = output_path) # added: burke damages (annual)
    discontinuity_impact_per_cap_sum = load_RV("isat_per_cap_DiscImpactperCapinclSaturation_sum"; output_path = output_path) # added: discontinuity damages (timestep)
    discontinuity_impact_per_cap_ann_sum = load_RV("isat_per_cap_DiscImpactperCapinclSaturation_ann_sum"; output_path = output_path) # added: discontinuity damages (annual)
    nonmmarket_impact_per_cap_ann_sum = load_RV("isat_per_cap_ImpactperCapinclSaturationandAdaptation_ann_sum_nonmarket"; output_path = output_path) # added: burke damages (annual)

    #time index
    c_co2concentration = load_RV("c_CO2concentration"; output_path = output_path)
    ft      = load_RV("ft_totalforcing"; output_path = output_path)
    rt_g    = load_RV("rt_g_globaltemperature"; output_path = output_path)
    s       = load_RV("s_sealevel"; output_path = output_path)

    #region index
    rgdppercap_slr          = load_RV("rgdp_per_cap_SLRRemainGDP"; output_path = output_path)
    rgdppercap_market       = load_RV("rgdp_per_cap_MarketRemainGDP"; output_path = output_path)
    rgdppercap_nonmarket    = load_RV("rgdp_per_cap_NonMarketRemainGDP"; output_path = output_path)
    rconspercap_disc         = load_RV("rcons_per_cap_DiscRemainConsumption"; output_path = output_path)
    rconspercap_disc_sum         = load_RV("rcons_per_cap_DiscRemainConsumption_sum"; output_path = output_path) # NEW
    rconspercap_disc_ann_sum = load_RV("rcons_per_cap_DiscRemainConsumption_ann_sum"; output_path = output_path) # NEW
    rgdppercap_market_sum       = load_RV("rgdp_per_cap_MarketRemainGDP_sum"; output_path = output_path) # NEW
    rgdppercap_market_ann_sum       = load_RV("rgdp_per_cap_MarketRemainGDP_ann_sum"; output_path = output_path) # NEW
    rgdppercap_nonmarket_sum    = load_RV("rgdp_per_cap_NonMarketRemainGDP_sum"; output_path = output_path) # NEW
    rgdppercap_nonmarket_ann_sum    = load_RV("rgdp_per_cap_NonMarketRemainGDP_ann_sum"; output_path = output_path) # NEW

    #resave data
    df=DataFrame(td=td,
            tpc=tpc,
            tac=tac,
            te=te,
            c_co2concentration=c_co2concentration,
            ft=ft,
            rt_g=rt_g,
            sealevel=s,
            rgdppercap_slr=rgdppercap_slr,
            rgdppercap_market=rgdppercap_market,
            rgdppercap_nonmarket=rgdppercap_nonmarket,
            rconspercap_di=rconspercap_disc,
            td_ann=td_ann,
            rconspercap_disc_sum=rconspercap_disc_sum,
            rconspercap_di_ann=rconspercap_disc_ann_sum,
            occurdis_occurrencedummy_sum=occurdis_occurrencedummy_sum,
            occurdis_occurrencedummy_ann_sum=occurdis_occurrencedummy_ann_sum,
            burke_impact_per_cap_sum=burke_impact_per_cap_sum,
            burke_impact_per_cap_ann_sum=burke_impact_per_cap_ann_sum,
            discontinuity_impact_per_cap_sum=discontinuity_impact_per_cap_sum,
            discontinuity_impact_per_cap_ann_sum=discontinuity_impact_per_cap_ann_sum,
            rgdppercap_market_sum=rgdppercap_market_sum,
            rgdppercap_market_ann_sum=rgdppercap_market_ann_sum,
            rgdppercap_nonmarket_sum=rgdppercap_nonmarket_sum,
            rgdppercap_nonmarket_ann_sum=rgdppercap_nonmarket_ann_sum,
            nonmmarket_impact_per_cap_ann_sum=nonmmarket_impact_per_cap_ann_sum)

    save(joinpath(output_path, "mimipagemontecarlooutput.csv"),df)
end

function do_monte_carlo_runs(samplesize::Int, output_path::String = joinpath(@__DIR__, "../output"), scenario::String = "NDCs")

    # get simulation
    mcs = getsim()

    # get a model
    m = getpage(scenario)
    run(m)

    # Generate trial data for all RVs and save to a file
    generate_trials!(mcs, samplesize, filename = joinpath(output_path, "trialdata.csv"))

    # set model
    set_models!(mcs, m)

    # Run trials 1:samplesize, and save results to the indicated directory, one CSV file per RV
    run_sim(mcs, output_dir = output_path)

    # reformat outputs for testing and analysis
    reformat_RV_outputs(samplesize, output_path=output_path)
end

function get_scc_mcs(samplesize::Int, year::Int, output_path::String = joinpath(@__DIR__, "../output"), scenario::String = "NDCs";
                     eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing,
                     pulse_size::Union{Float64, Nothing} = 2000.)
    # Setup the marginal model
    m = getpage(scenario)
    mm = compute_scc_mm(m, year=year, eta=eta, prtp=prtp, pulse_size=pulse_size)[:mm]

    # Setup SCC calculation and place for results
    scc_results = zeros(samplesize)

    function my_scc_calculation(mcs::Simulation, trialnum::Int, ntimesteps::Int, tup::Union{Tuple, Nothing})
        base, marginal = mcs.models
        if use_annual
            scc_results[trialnum] = (marginal[:EquityWeighting, :td_totaldiscountedimpacts_ann] - base[:EquityWeighting, :td_totaldiscountedimpacts_ann]) / pulse_size # CHANGED: made to new annual variable
        else
            scc_results[trialnum] = (marginal[:EquityWeighting, :td_totaldiscountedimpacts] - base[:EquityWeighting, :td_totaldiscountedimpacts]) / pulse_size
        end
    end

    # Setup MC simulation
    mcs = getsim()
    set_models!(mcs, [mm.base, mm.marginal])
    generate_trials!(mcs, samplesize, filename = joinpath(output_path, "scc_trials.csv"))

    # Run it!
    run_sim(mcs, output_dir=output_path, post_trial_func=my_scc_calculation)

    writedlm(joinpath(output_path, "scc.csv"), scc_results, ',')

end

### To run the general model in MC mode:
# do_monte_carlo_runs(2000, dir_output) # 500 iterations takes about 40 seconds...

### To calculate the SSC by means of MC sampling
# include("mcs.jl")
# do_monte_carlo_runs(100)
# include("compute_scc.jl")
# get_scc_mcs(100, 2020, dir_output)

max_stars_repo_name: jkikstra/PAGE-VAR

max_issues_repo_name: jkikstra/PAGE-VAR

max_forks_repo_name: jkikstra/PAGE-VAR


text:
using Manifolds, LinearAlgebra, PGFPlotsX, Colors, Distributions, Contour, Random

#
# Settings
#
dark_mode = false

line_offset_brightness = 0.25
patch_opacity = 1.0
geo_opacity = dark_mode ? 0.66 : 0.5
geo_line_width = 35
mesh_line_width = 5
mesh_opacity = dark_mode ? 0.5 : 0.7
logo_colors = [(77, 100, 174), (57, 151, 79), (202, 60, 50), (146, 89, 163)] # Julia colors

rgb_logo_colors = map(x -> RGB(x ./ 255...), logo_colors)
rgb_logo_colors_bright = map(x -> RGB( (1 + line_offset_brightness) .* x ./ 255...), logo_colors)
rgb_logo_colors_dark = map(x -> RGB( (1 - line_offset_brightness) .* x ./ 255...), logo_colors)

out_file_prefix = dark_mode ? "logo-dark" : "logo"
out_file_ext = ".svg"

#
# Helping functions
#
polar_to_cart(r, θ) = (r * cos(θ), r * sin(θ))

cart_to_polar(x, y) = (hypot(x, y), atan(y, x))

normal_coord_to_vector(M, x, rθ, B) = get_vector(M, x, collect(polar_to_cart(rθ...)), B)

normal_coord_to_point(M, x, rθ, B) = exp(M, x, normal_coord_to_vector(M, x, rθ, B))

function plot_normal_coord!(
    ax,
    M,
    x,
    B,
    rs,
    θs;
    ncirc = 9,
    options = Dict(),
    kwargs...,
)
    for r in rs[2:end-1]
        push!(ax, Plot3(
            options,
            Coordinates(map(θ -> Tuple(normal_coord_to_point(M, x, [r, θ], B)), θs))
        ))
    end
    for θ in range(0, 2π; length = ncirc)
        push!(ax, Plot3(
            options,
            Coordinates(map(r -> Tuple(normal_coord_to_point(M, x, [r, θ], B)), rs))
        ))
    end
    return ax
end

function plot_patch!(ax, M, x, B, r, θs; options = Dict())
    push!(ax, Plot3(
        options,
        Coordinates(map(θ -> Tuple(normal_coord_to_point(M, x, [r, θ], B)), θs))
    ))
    return ax
end

function plot_geodesic!(ax, M, x, y; n = 100, options = Dict())
    γ = shortest_geodesic(M, x, y)
    T = range(0, 1; length = n)
    push!(ax, Plot3(options, Coordinates(Tuple.(γ(T)))))
    return ax
end

#
# Prepare document
#
resize!(PGFPlotsX.CUSTOM_PREAMBLE, 0)
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\pgfplotsset{scale=6.0}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usetikzlibrary{arrows.meta}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\pgfplotsset{roundcaps/.style={line cap=round}}")
push!(
    PGFPlotsX.CUSTOM_PREAMBLE,
    raw"\pgfplotsset{meshlinestyle/.style={dash pattern=on 0pt off 2.5\pgflinewidth, line cap=round}}",
)
if dark_mode
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\pagecolor{black}")
end
S = Sphere(2)

center = normalize([1, 1, 1])
x, y, z = eachrow(Matrix{Float64}(I, 3, 3))
γ1 = shortest_geodesic(S, center, z)
γ2 = shortest_geodesic(S, center, x)
γ3 = shortest_geodesic(S, center, y)
p1 = γ1(1)
p2 = γ2(1)
p3 = γ3(1)

#
# Setup Axes
if dark_mode
    tp = @pgf Axis({
        axis_lines = "none",
        axis_equal,
        view = "{135}{35}",
        zmin = -0.05, zmax = 1.0, xmin = 0.0, xmax = 1.0, ymin = 0.0, ymax = 1.0
    })
else
    tp = @pgf Axis({
        axis_lines = "none",
        axis_equal,
        view = "{135}{35}",
        zmin = -0.05, zmax = 1.0, xmin = 0.0, xmax = 1.0, ymin = 0.0, ymax = 1.0
    })
end
rs = range(0, π / 5; length = 6)
θs = range(0, 2π; length = 100)

#
# Plot manifold patches
patch_colors = rgb_logo_colors[2:end]
patch_colors_line = dark_mode ? rgb_logo_colors_bright[2:end] : rgb_logo_colors_dark[2:end]
base_points = [p1, p2, p3]
basis_vectors = [log(S, p1, p2), log(S, p2, p1), log(S, p3, p1)]
for i in eachindex(base_points)
    x = base_points[i]
    B = DiagonalizingOrthonormalBasis(basis_vectors[i])
    basis = get_basis(S, x, B)
    optionsP = @pgf {fill = patch_colors[i], draw = "none", opacity = patch_opacity}
    plot_patch!(tp, S, x, basis, π / 5, θs; options = optionsP)
    optionsL = @pgf {
        meshlinestyle,
        color = dark_mode ? "white" : "black",
        line_width = mesh_line_width,
        opacity = mesh_opacity,
    }
    plot_normal_coord!(tp, S, x, basis, rs, θs; options = optionsL)
end

#
# Plot geodesics
options = @pgf {
    opacity=geo_opacity,
    no_markers,
    roundcaps,
    line_width = geo_line_width,
    color = dark_mode ? "white" : "black"}
plot_geodesic!(tp, S, base_points[1], base_points[2]; options = options)
plot_geodesic!(tp, S, base_points[1], base_points[3]; options = options)
plot_geodesic!(tp, S, base_points[2], base_points[3]; options = options)

#
# Export Logo.
out_file = "$(out_file_prefix)$(out_file_ext)"
pgfsave(out_file, tp)

ext: jl


text:
struct MultipleShooting{T, N, Ny, Nx, Nθ, M}
    simulations
    list_N
end

function MultipleShooting{T}(f::Function, g::Function,
                             x0_list::Vector{T}, y::Vector{T},
                             k0_list::Vector{Int}, θ,
                             loss=L2DistLoss())
    # Get lengths
    M = length(x0_list)
    N = length(y)
    Ny = length(y[1])
    Nx = length(x0_list[1])
    Nθ = length(θ)
    # Simulate
    k0 = k0_list[1]
    grid = [k0_list; k0+N]
    list_N = Vector{Int}(M)
    simulations = Vector{Union{Future,OneShootSimulation}}(M)
    for i =1:M
        x0 = x0_list[i]; k0 = k0_list[i];
        yi = y[grid[i]:grid[i+1]-1]; list_N[i] = length(yi);
        simulations[i] = OneShootSimulation(f, g, x0, yi, k0, θ, loss)

    end
    MultipleShooting{T, N, Ny, Nx, Nθ, M}(simulations, list_N)
end

function new_simulation!{T, N, Ny, Nx, Nθ, M}(
        ms::MultipleShooting{T, N, Ny, Nx, Nθ, M}, x0_list::Vector{T}, θ)
    for i = 1:M
        x0 = x0_list[i]
        ms.simulations[i] = new_simulation!(ms.simulations[i], x0, θ)
    end
    return ms
end

function cost_function{T, N, Ny, Nx, Nθ, M}(
        ms::MultipleShooting{T, N, Ny, Nx, Nθ, M})
    cost = 0
    for i = 1:M
        cost += cost_function(ms.simulations[i])
    end
    return cost
end

function derivatives_θ!{T, N, Ny, Nx, Nθ, M}(
        dvec::Vector{Float64},
        ms::MultipleShooting{T, N, Ny, Nx, Nθ, M},
        accumulat=false, dtype="gradient",
        p::Vector{Float64}=Float64[])
    variable="θ"
    # Set initial values to zero if not accumulating
    if !accumulat
        fill!(dvec, 0.0)
    end
    # Make all the computations on remote instances
    for i = 1:M
        if dtype == "gradient"
            gradient!(dvec, ms.simulations[i],
                      variable, true)
        elseif dtype == "hessian_approx"
            hessian_approx!(dvec, ms.simulations[i], p,
                           variable, true)
        end
    end
    return dvec
end

function derivatives_x0!{T, N, Ny, Nx, Nθ, M}(
        dvec::Vector{Vector{Float64}},
        ms::MultipleShooting{T, N, Ny, Nx, Nθ, M},
        accumulat=false, dtype="gradient",
        p::Vector{Vector{Float64}}=Vector{Float64}[])
    variable="x0"
    # Set initial values to zero if not accumulating
    if !accumulat
        for i = 1:M
            fill!(dvec[i], 0.0)
        end
    end
    # Make all the computations on remote instances
    for i = 1:M
        if dtype == "gradient"
            gradient!(dvec[i], ms.simulations[i],
                      variable, false)
        elseif dtype == "hessian_approx"
            hessian_approx!(dvec[i], ms.simulations[i],
                            p[i], variable, false)
        end
    end
    return dvec
end

function gradient!{T, N, Ny, Nx, Nθ, M}(
        grad, ms::MultipleShooting{T, N, Ny, Nx, Nθ, M},
        variable="θ", accumulat=false)
    if variable == "θ"
        derivatives_θ!(grad, ms, accumulat, "gradient")
    elseif variable == "x0"
        derivatives_x0!(grad, ms, accumulat, "gradient")
    end
end

function hessian_approx!{T, N, Ny, Nx, Nθ, M}(
        hessp, ms::MultipleShooting{T, N, Ny, Nx, Nθ, M}, p,
        variable="θ", accumulat=false)
    if variable == "θ"
        derivatives_θ!(hessp, ms, accumulat, "hessian_approx", p)
    elseif variable == "x0"
        derivatives_x0!(hessp, ms, accumulat, "hessian_approx", p)
    end
end

function constr!{T, N, Ny, Nx, Nθ, M}(
        c::Vector{Vector{Float64}},
        ms::MultipleShooting{T, N, Ny, Nx, Nθ, M})
    # Get x
    for i = 1:M-1
        copy!(c[i], ms.simulations[i].x)
    end
    # subtract x0
    for i = 2:M
        Base.LinAlg.axpy!(-1, ms.simulations[i].x0, c[i-1])
    end
    return c
end

function constr_jac!{T, N, Ny, Nx, Nθ, M}(
        jac::Vector{Matrix{Float64}},
        ms::MultipleShooting{T, N, Ny, Nx, Nθ, M},
        variable="θ")
    # Get x
    for i = 1:M-1
        if variable == "x0"
            jac[i] .= get_dxdx0(ms.simulations[i])
        elseif variable == "θ"
            jac[i] .= get_dxdθ(ms.simulations[i])
        end
    end
    return jac
end


text:
function _extract_reaction_phrase_table(phrase::String)::Dict{String,Float64}

    # initialize -
    reaction_phrase_table = Dict{String,Float64}()

    # compnents -
    if (occursin("+",phrase) == true)

        # ok, so we have some +'s
        component_array = split(phrase,"+")
        for component in component_array

            # split around the * -
            inner_component_array = split(component,"*")

            # if len(inner_component_array) == 1, then no *
            # if len(inner_component_array) == 2, then we have a st coeff -
            if (length(inner_component_array) == 1)

                # we have a single species -
                key = string(last(inner_component_array))
                value = 1.0
                reaction_phrase_table[key] = value

            elseif (length(inner_component_array) == 2)

                # we have a species w/a coefficient -
                key = string(last(inner_component_array))
                value = parse(Float64, first(inner_component_array))
                reaction_phrase_table[key] = value
            else
                # error state -
                throw(DimensionMismatch("component array larger than 2"))
            end
        end
    else

        # single species, possibly w/*
        # split around the * -
        phrase_component_array = split(phrase,"*")
        if (length(phrase_component_array) == 1)

            # we have a single species -
            key = string(last(phrase_component_array))
            value = 1.0
            reaction_phrase_table[key] = value

        elseif (length(phrase_component_array) == 2)

            # we have a species w/a coefficient -
            key = string(last(phrase_component_array))
            value = parse(Float64, first(phrase_component_array))
            reaction_phrase_table[key] = value

        else
            # error state -
            throw(DimensionMismatch("phrase component array larger than 2"))
        end
    end

    # return -
    return reaction_phrase_table
end

function _extract_stoichiometric_coefficient(phrase::String, speciesSymbol::String)::Float64

    # initialize -
    stoichiometric_coefficient = 0.0 # default is 0, species in NOT involved

    # get the reaction phrase table -
    reaction_phrase_table = _extract_reaction_phrase_table(phrase)

    # ok, look up the species in the phrase table, if its there we have a match - no = 0
    if (haskey(reaction_phrase_table,speciesSymbol) == true)
        return reaction_phrase_table[speciesSymbol]
    end

    # return -
    return stoichiometric_coefficient
end

# == PUBLIC METHODS BELOW HERE ======================================================================================== #

"""
    generate_stoichiometric_matrix(intermediate_dictionary::Dict{String,Any})::VLResult

Generate a stoichiometric matrix based on the biochemical model reaction network.

Input arguments:
`intermediate_dictionary::Dict{String,Any}` - data dictionary containing the master reaction table and molecular species participating in the reactions.

Output arguments:
`VLResult::VLResult` - concrete data type holding the generated stoichiometric matrix.

"""
function generate_stoichiometric_matrix(intermediate_dictionary::Dict{String,Any})::VLResult

    # initialize -

    try

        # check: do we have the ir_master_reaction_table_key?
        if (haskey(intermediate_dictionary, ir_master_reaction_table_key) == false)
            throw(ArgumentError("intermediate_dictionary is missing the ir_master_reaction_table_key"))
        end

        # check: do we have the ir_list_of_molecular_species_key?
        if (haskey(intermediate_dictionary, ir_list_of_molecular_species_key) == false)
            throw(ArgumentError("intermediate_dictionary is missing the ir_list_of_molecular_species_key"))
        end

        # get reaction data frame, and species list -
        reaction_data_frame = intermediate_dictionary[ir_master_reaction_table_key]
        list_of_species = intermediate_dictionary[ir_list_of_molecular_species_key]

        # ok, so we need to initialize some space for the stm -
        number_of_species = length(list_of_species)
        number_of_reactions = size(reaction_data_frame,1)
        stoichiometric_matrix = Array{Float64,2}(undef, number_of_species, number_of_reactions)

        # process each species -
        for (species_index,species_symbol) in enumerate(list_of_species)

            # process each reaction -
            for reaction_index = 1:number_of_reactions

                # get the left and right phrases -
                left_phrase = reaction_data_frame[reaction_index,:left_phrase]
                right_phrase = reaction_data_frame[reaction_index,:right_phrase]

                # check: is this species symbol in the left or right phrase?
                left_st_coeff = _extract_stoichiometric_coefficient(left_phrase, species_symbol)
                right_st_coeff = _extract_stoichiometric_coefficient(right_phrase, species_symbol)

                # overall coeff -
                overall_st_coeff = right_st_coeff - left_st_coeff

                # grab -
                stoichiometric_matrix[species_index,reaction_index] = overall_st_coeff
            end
        end

        # return -
        return VLResult(stoichiometric_matrix)
    catch error
        return VLResult(error)
    end
end

# == PUBLIC METHODS ABOVE HERE ======================================================================================== #

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
# Note: these function live in the Resampling module, which extends the SamplingConstraints methods

# Add uncertainties to the time series values
n_points = 40
x_uncertain = [UncertainValue(Normal, x, rand(Uniform(0.1, 0.8))) for x in rand(n_points)]
y_uncertain = [UncertainValue(Normal, y, rand(Uniform(0.1, 0.8))) for y in rand(n_points)]
x = UncertainValueDataset(x_uncertain)
y = UncertainValueDataset(y_uncertain)

time_uncertain = [UncertainValue(Normal, i, 1) for i = 1:length(x)];
time_certain = [CertainValue(i) for i = 1:length(x)];
timeinds_x = UncertainIndexDataset(time_uncertain)
timeinds_y = UncertainIndexDataset(time_certain)

X = UncertainIndexValueDataset(timeinds_x, x)
Y = UncertainIndexValueDataset(timeinds_y, y);

# Truncate each of the indices for x at 0.8 their standard deviation around the mean
constraints_x_inds = TruncateStd(0.8)

# Truncate each of the indices for y at 1.5 their standard deviation around the mean
constraints_y_inds = TruncateStd(1.5)

# Truncate each of the values of x at the 20th percentile range
constraints_x_vals = [TruncateQuantiles(0.4, 0.6) for i = 1:length(x)];

# Truncate each of the values of x at the 80th percentile range
constraints_y_vals = [TruncateQuantiles(0.1, 0.9) for i = 1:length(x)];

@test constrain(X.indices, ConstrainedValueResampling(constraints_x_inds)) isa ConstrainedUncertainIndexDataset
@test constrain(X.values, ConstrainedValueResampling(constraints_x_vals)) isa ConstrainedUncertainValueDataset
@test constrain(X, ConstrainedValueResampling(constraints_x_inds), ConstrainedValueResampling(constraints_x_vals)) isa UncertainIndexValueDataset

idxval_resampling = ConstrainedIndexValueResampling((constraints_x_inds, constraints_x_vals))
@test constrain(X, idxval_resampling) isa UncertainIndexValueDataset
@test constrain(X, constraints_x_inds, constraints_x_vals) isa UncertainIndexValueDataset


text:
module MyExample

using ForwardDiff
# Write your package code here.

include("extra_file.jl")

export my_f, derivative_of_my_f
end

max_stars_repo_name: filippigianluca/MyExample

max_issues_repo_name: filippigianluca/MyExample


text:
#  Copyright 2016,2017 Eric S. Tellez <eric.tellez@infotec.mx>
#
#    Licensed under the Apache License, Version 2.0 (the "License");
#    you may not use this file except in compliance with the License.
#    You may obtain a copy of the License at
#
#        http://www.apache.org/licenses/LICENSE-2.0
#
#    Unless required by applicable law or agreed to in writing, software
#    distributed under the License is distributed on an "AS IS" BASIS,
#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#    See the License for the specific language governing permissions and
#    limitations under the License.

using SimilaritySearch
using Dates

import SimilaritySearch:
    search

import Base:
    push!

export Knr, optimize!

mutable struct Knr{T, D} <: Index
    db::Vector{T}
    dist::D
    refs::Vector{T}
    k::Int
    ksearch::Int
    minmatches::Int
    invindex::Vector{Vector{Int32}}
end

function Knr(db::Vector{T}, dist::D, refs::Vector{T}, k::Int, minmatches::Int=1) where {T,D}
    @info "Knr> refs=$(typeof(db)), k=$(k), numrefs=$(length(refs)), dist=$(dist)"
    invindex = [Vector{Int32}(undef, 0) for i in 1:length(refs)]
    seqindex = Sequential(refs, dist)

    pc = round(Int, length(db) / 20)
    for i=1:length(db)
        if i % pc == 0
            @info "Knr> advance $(round(i/length(db), digits=4)), now: $(now())"
        end

        res = search(seqindex, db[i], KnnResult(k))
        for p in res
            push!(invindex[p.objID], i)
        end
    end

    Knr(db, dist, refs, k, k, minmatches, invindex)
end

function Knr(db::Array{T,1}, dist::D; numrefs::Int=1024, k::Int=7, minmatches::Int=1, tournamentsize::Int=3) where {T,D}
    # refs = rand(db, numrefs)
    refs = [db[x] for x in select_tournament(db, dist, numrefs, tournamentsize)]
    Knr(db, dist, refs, k, minmatches)
end

"""
    search(index::Knr{T, D}, q::T) where {T,D}

Solves the search specified by `q`and `res` using `index`
"""
function search(index::Knr{T,D}, q::T, res::Result) where {T,D}
    dz = zeros(Int16, length(index.db))
    # M = BitArray(length(index.db))
    seqindex = Sequential(index.refs, index.dist)
    kres = search(seqindex, q, KnnResult(index.ksearch))

    for p in kres
        @inbounds for objID in index.invindex[p.objID]
            c = dz[objID] + 1
            dz[objID] = c

            if c == index.minmatches
                d = index.dist(q, index.db[objID])
                push!(res, objID, d)
            end
        end
    end

    return res
end

function push!(index::Knr{T, D}, obj::T) where {T,D}
    push!(index.db, obj)
    seqindex = Sequential(index.refs, index.dist)
    res = search(seqindex, obj, KnnResult(index.k))
    for p in res
        push!(index.invindex[p.objID], length(index.db))
    end
    return length(index.db)
end

function optimize!(index::Knr{T, D}; recall::Float64=0.9, k::Int=1, numqueries::Int=128, use_distances::Bool=false) where {T,D}
    @info "Knr> optimizing index for recall=$(recall)"
    perf = Performance(index.db, index.dist; numqueries=numqueries, expected_k=k)
    index.minmatches = 1
    index.ksearch = 1
    p = probe(perf, index, use_distances=use_distances)

    while p.recall < recall && index.ksearch < length(index.refs)
        index.ksearch += 1
        @info "Knr> opt step ksearch=$(index.ksearch), performance $(p)"
        p = probe(perf, index, use_distances=use_distances)

    end
    @info "Knr> reached performance $(p)"
    return index
end


text:
# This file contains functions for the Modified Directional Distance Function DEA model
"""
    ModifiedDDFDEAModel
An data structure representing a modified directional distance function DEA model.
"""
struct ModifiedDDFDEAModel <: AbstractTechnicalDEAModel
    n::Int64
    m::Int64
    s::Int64
    Gx::Symbol
    Gy::Symbol
    rts::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    betax::Vector
    betay::Vector
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deamddf(X, Y; Gx, Gy)
Compute data envelopment analysis modified directional distance function model for inputs
`X` and outputs `Y`, using directions `Gx` and `Gy`.

# Direction specification:

The directions `Gx` and `Gy` can be one of the following symbols.
- `:Ones`: use ones.
- `:Observed`: use observed values.
- `:Mean`: use column means.

Alternatively, a vector or matrix with the desired directions can be supplied.

# Optional Arguments
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `slack=true`: computes input and output slacks.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [2; 4; 8; 12; 6; 14; 14; 9.412];

julia> Y = [1; 5; 8; 9; 3; 7; 9; 2.353] ;

julia> deamddf(X, Y, Gx = :Ones, Gy = :Ones)
Modified DDF DEA Model
DMUs = 8; Inputs = 1; Outputs = 1
Returns to Scale = CRS
Gx = Ones; Gy = Ones
────────────────────────────────────────────
   efficiency   βx      βy  slackX1  slackY1
────────────────────────────────────────────
1       1.5    0.0   1.5        0.0      0.0
2       0.0    0.0   0.0        0.0      0.0
3       2.0    0.0   2.0        0.0      0.0
4       6.0    0.0   6.0        0.0      0.0
5       4.5    0.0   4.5        0.0      0.0
6      10.5    0.0  10.5        0.0      0.0
7       8.5    0.0   8.5        0.0      0.0
8       9.412  0.0   9.412      0.0      0.0
────────────────────────────────────────────
```
"""
function deamddf(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
    Gx::Union{Symbol, Matrix, Vector}, Gy::Union{Symbol, Matrix, Vector},
    rts::Symbol = :CRS, slack::Bool = true,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix,Vector,Nothing} = nothing,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::ModifiedDDFDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nrefx != nrefy
        throw(DimensionMismatch("number of rows in Xref and Yref ($nrefx, $nrefy) are not equal"));
    end
    if m != mref
        throw(DimensionMismatch("number of columns in X and Xref ($m, $mref) are not equal"));
    end
    if s != sref
        throw(DimensionMismatch("number of columns in Y and Yref ($s, $sref) are not equal"));
    end

    # Build or get user directions
    if typeof(Gx) == Symbol
        Gxsym = Gx

        if Gx == :Ones
            Gx = ones(size(X))
        elseif Gx == :Observed
            Gx = X
        elseif Gx == :Mean
            Gx = repeat(mean(X, dims = 1), size(X, 1))
        else
            throw(ArgumentError("Invalid `Gx`"));
        end

    else
        Gxsym = :Custom
    end

    if typeof(Gy) == Symbol
        Gysym = Gy

        if Gy == :Ones
            Gy = ones(size(Y))
        elseif Gy == :Observed
            Gy = Y
        elseif Gy == :Mean
            Gy = repeat(mean(Y, dims = 1), size(Y, 1))
        else
            throw(ArgumentError("Invalid `Gy`"));
        end

    else
        Gysym = :Custom
    end

    nGx, mGx = size(Gx, 1), size(Gx, 2)
    nGy, sGy = size(Gy, 1), size(Gy, 2)

    if (size(Gx, 1) != size(X, 1)) | (size(Gx, 2) != size(X, 2))
        throw(DimensionMismatch("size of Gx and X ($(size(Gx)), $(size(X))) are not equal"));
    end
    if (size(Gy, 1) != size(Y, 1)) | (size(Gy, 2) != size(Y, 2))
        throw(DimensionMismatch("size of Gy and Y ($(size(Gy)), $(size(Y))) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Compute efficiency for each DMU
    n = nx
    nref = nrefx

    effi = zeros(n)
    betaxi = zeros(n)
    betayi = zeros(n)
    lambdaeff = spzeros(n, nref)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]

        # Directions to use
        Gx0 = Gx[i,:]
        Gy0 = Gy[i,:]

        # Solve if any direction is different from zero
        if any(Gx0 .!= 0) | any(Gy0 .!= 0)
            # Create the optimization model
            deamodel = newdeamodel(optimizer)
            set_silent(deamodel)

            @variable(deamodel, betax >= 0)
            @variable(deamodel, betay >= 0)
            @variable(deamodel, lambda[1:nref] >= 0)

            @objective(deamodel, Max, betax + betay)

            @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) <= x0[j] - betax * Gx0[j])
            @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) >= y0[j] + betay * Gy0[j])

            # Add return to scale constraints
            if rts == :CRS
                # No contraint to add for constant returns to scale
            elseif rts == :VRS
                @constraint(deamodel, sum(lambda) == 1)
            else
                throw(ArgumentError("`rts` must be :CRS or :VRS"));
            end

            # Optimize and return results
            JuMP.optimize!(deamodel)

            effi[i]  = JuMP.objective_value(deamodel)
            betaxi[i] = JuMP.value(betax)
            betayi[i] = JuMP.value(betay)
            lambdaeff[i,:] = JuMP.value.(lambda)

            # Check termination status
            if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
                @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
            end
        else
            effi[i]  = 0.0
            betaxi[i] = 0.0
            betayi[i] = 0.0
            lambdaeff[i,:] .= 0.0
            lambdaeff[i,i] = 1.0
        end

    end

    # Get first-stage X and Y targets
    Xtarget = X .- betaxi .* Gx
    Ytarget = Y .+ betayi .* Gy

    # Compute slacks
    if slack == true
        # Compute slacks if more than one input or more than 1 output
        if (m > 1) || (s > 1)
            # Use additive model with X and Y targets to get slacks
            slacksmodel = deaadd(Xtarget, Ytarget, :Ones, rts = rts, Xref = Xref, Yref = Yref, optimizer = optimizer)
            slackX = slacks(slacksmodel, :X)
            slackY = slacks(slacksmodel, :Y)

            # Get second-stage X and Y targets
            Xtarget = Xtarget - slackX
            Ytarget = Ytarget + slackY
        else
            slackX = zeros(n, m)
            slackY = zeros(n, s)

            if typeof(Xtarget) <: AbstractVector    Xtarget = Xtarget[:,:]  end
            if typeof(Ytarget) <: AbstractVector    Ytarget = Ytarget[:,:]  end
        end
    else
        if typeof(Xtarget) <: AbstractVector    Xtarget = Xtarget[:,:]  end
        if typeof(Ytarget) <: AbstractVector    Ytarget = Ytarget[:,:]  end

        slackX = Array{Float64}(undef, 0, 0)
        slackY = Array{Float64}(undef, 0, 0)
    end

    return ModifiedDDFDEAModel(n, m, s, Gxsym, Gysym, rts, names, effi, betaxi, betayi, slackX, slackY, lambdaeff, Xtarget, Ytarget)

end

function Base.show(io::IO, x::ModifiedDDFDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    eff = efficiency(x)
    dmunames = names(x)

    betax = x.betax
    betay = x.betay
    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)
    hasslacks = ! isempty(slackX)

    if !compact
        print(io, "Modified DDF DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Returns to Scale = ", string(x.rts))
        print(io, "\n")
        print(io, "Gx = ", string(x.Gx), "; Gy = ", string(x.Gy))
        print(io, "\n")

        if hasslacks == true
            show(io, CoefTable(hcat(eff, betax, betay, slackX, slackY), ["efficiency"; "βx"; "βy"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
        else
            show(io, CoefTable(hcat(eff, betax, betay), ["efficiency"; "βx"; "βy"], dmunames))
        end
    end

end

function efficiency(model::ModifiedDDFDEAModel, type::Symbol)::Vector

    if type == :X return model.betax end
    if type == :Y return model.betay end

    throw(ArgumentError("$(typeof(model)) has no efficiency $(type)"));

end


text:
# Copyright 2020 Paul Melis (paul.melis@surf.nl)
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

mutable struct HalfEdge
    source  ::UInt32        # Vertex index
    target  ::UInt32        # Vertex index
    face    ::UInt32        # Face index

    index   ::UInt32        # Should be same value for sibling edge

    next    ::HalfEdge
    prev    ::HalfEdge

    # Can be nothing for a boundary edge
    sibling ::Union{HalfEdge, Nothing}

    HalfEdge(s, t, f, i) = new(s, t, f, i)
end

Base.show(io::IO, he::HalfEdge) = print(io,
    "HalfEdge([", he.index, "] ", he.source, " -> ", he.target, ", f:", he.face, ", ",
    he.sibling != nothing ? "s" : "!s",
    ")")

const VertexPair = Tuple{UInt32, UInt32}

function build(vertices, loop_start, loop_total, loops)

    # Temporary index
    half_edges = Dict{VertexPair, HalfEdge}()

    # For each vertex an arbitrary half-edge starting from that vertex
    vertex_start_edges = Dict{UInt32, HalfEdge}()
    # For each face an arbitrary half-edge making up that face
    face_start_edges = Dict{UInt32, HalfEdge}()

    # For each polygon edge (not half-edge) one of the half-edges making up
    # that polygon edge
    edges = HalfEdge[]
    num_edges = 0

    num_vertices = trunc(UInt32, length(vertices) / 3)
    num_faces = length(loop_start)

    for fi = 1:num_faces

        ls = loop_start[fi]
        nv = loop_total[fi]
        loop = loops[ls:ls+nv-1]

        first_he = nothing
        last_halfedge = nothing

        for i = 1:nv
            cur = i
            next = i % nv + 1

            A = loop[cur]
            B = loop[next]

            key = (A,B)
            rkey = (B,A)

            if haskey(half_edges, rkey)
                he2 = half_edges[rkey]
                he = HalfEdge(A, B, fi, he2.index)
                he2.sibling = he
                he.sibling = he2
            else
                num_edges += 1
                he = HalfEdge(A, B, fi, num_edges)
                he.sibling = nothing
                push!(edges, he)
            end

            @assert !haskey(half_edges, key)
            half_edges[key] = he

            if i == 1
                face_start_edges[fi] = he
                first_he = he
            end

            if !haskey(vertex_start_edges, A)
                vertex_start_edges[A] = he
            end

            if last_halfedge != nothing
                last_halfedge.next = he
                he.prev = last_halfedge
            end

            last_halfedge = he
        end

        last_halfedge.next = first_he
        first_he.prev = last_halfedge

    end

    # XXX this will trigger on models with unconnected vertices, but those
    # shouldn't really hurt later on
    #@assert length(vertex_start_edges) == num_vertices

    return num_vertices, num_faces, num_edges, vertex_start_edges, face_start_edges, edges

end

max_stars_repo_name: paulmelis/blender-julia-test

max_issues_repo_name: paulmelis/blender-julia-test

max_forks_repo_name: paulmelis/blender-julia-test


text:
using QuantileRegressions, CSV, DataFrames
@testset "npqreg" begin
	mycyle = CSV.read("../examples/mcycle.csv", DataFrame)
	x = mycyle.times
	y = mycyle.accel

	for t = range(0.01,0.99,length=99)
		res = QuantileRegressions.npqreg(Array(y), Array(x), t; h=2)
		@test all(!isnan, res[2])
		@test all(!isnan, res[3])
		xrange = sort(unique(x))
		res = QuantileRegressions.npqreg(Array(y), Array(x), t; xrange=xrange, h=0.5)
	end
	res = QuantileRegressions.npqreg([y[1]], [x[1]], 0.4; xrange=[x[1]], h=2)
	@test !isnan(res[2][1])
	@test isnan(res[3][1])
	res = QuantileRegressions.npqreg([y[1]], [x[1]], 0.8; xrange=[x[1]], h=0.0001)
	@test !isnan(res[2][1])
	@test isnan(res[3][1])
	res = QuantileRegressions.npqreg([y[1]], [x[1]], 0.8; xrange=[x[1].+1], h=0.1)
	@test isnan(res[2][1])
	@test isnan(res[3][1])

	@test all(QuantileRegressions.npqreg([1.0,2.0,4.0,5.0], [0.0,0.0,6.0,12.0], 0.3; h=0.001, xrange=[0.0,6.0,12.0])[2] .≈ [1.0, 4.0, 5.0])

	yrand = randn(100)
	# Integer xrange gives same results as float xrange
	@test all(QuantileRegressions.npqreg(yrand, collect(1:100), 0.4; xrange=1:1:10) .== QuantileRegressions.npqreg(yrand, collect(1:100), 0.4; xrange=1.0:1.0:10.0))

	# Integer x gives same results as float x
	@test all(QuantileRegressions.npqreg(yrand, collect(1:100), 0.4) .== QuantileRegressions.npqreg(yrand, collect(1.0:100.0), 0.4))

	# Combinations
	@test all(QuantileRegressions.npqreg(yrand, collect(1.0:100.0), 0.4; xrange=1:1:10) .== QuantileRegressions.npqreg(yrand, collect(1:100), 0.4; xrange=1.0:1.0:10.0))
	@test all(QuantileRegressions.npqreg(yrand, collect(1:100), 0.4; xrange=1.0:1.0:10.0) .== QuantileRegressions.npqreg(yrand, collect(1:100), 0.4; xrange=1.0:1.0:10.0))
end