Merge pull request #22 from rffscghg/page-discontinuity

WIP discontinuity discrepancy in PAGE

README.md

@@ -67,10 +67,14 @@ MimiIWG.run_scc_mcs(MODEL,
     output_dir = nothing,  # Output directory. If unspecified, a directory with the following name will be created: "output/MODEL yyyy-mm-dd HH-MM-SS SC-$gas MC$trials"
     save_trials = false,   # Whether to save all of the input data sampled for each trial of the Monte Carlo Simulation. If true, values get saved to "output_dir/trials.csv"
     tables = true   # Whether to save a series of summary tables in the output folder; these include statistics such as percentiles and std errors of the SCC values.
+    drop_discontinuities = false # PAGE specific see below
+    save_md = false # Whether to save the global undiscounted marginal damages from each run of the simulation in a subdirectory "output/marginal_damages"
 )
 ```
 Note that the Monte Carlo Simulations are run across all five of the USG socioeconomics scenarios.
 
+If the `drop_discontinuities` optional keyword argument equals `true`, then outliers from the PAGE model (runs where discontinuity damages are triggered in different timesteps in the base and perturbed models) will not contribute to summary statistics. An additional folder "discontinuity_mismatch" contains files identifying in which runs the discrepencies occured. This is a PAGE-specific function that has no bearing on the results of other models, and will mostly be used for internal more advanced cases.
+
 ## Summary of modifications made by the IWG
 
 ### Socioeconomics scenarios

src/core/utils.jl

@@ -40,8 +40,9 @@ function write_scc_values(values, output_dir, perturbation_years, discount_rates
 end
 
 # helper function for computing percentile tables at the end of a monte carlo simulation
-function make_percentile_tables(output_dir, gas, discount_rates, perturbation_years)
+function make_percentile_tables(output_dir, gas, discount_rates, perturbation_years, drop_discontinuities=false)
     scc_dir = "$output_dir/SC-$gas"     # folder with output from the MCS runs
+    drop_discontinuities ? disc_dir = joinpath(output_dir, "discontinuity_mismatch/") : nothing
     tables = "$output_dir/Tables/Percentiles"   # folder to save TSD tables to
     mkpath(tables)
 
@@ -55,7 +56,11 @@ function make_percentile_tables(output_dir, gas, discount_rates, perturbation_ye
             write(f, "Scenario,1st,5th,10th,25th,50th,Avg,75th,90th,95th,99th\n")
             for fn in filter(x -> endswith(x, "$dr.csv"), results)  # Get the results files for this discount rate
                 scenario = split(fn)[1] # get the scenario name
-                d = readdlm(joinpath(scc_dir, fn), ',')[2:end, idx] # just keep 2020 values
+                d = readdlm(joinpath(scc_dir, fn), ',')[2:end, idx]
+                if drop_discontinuities
+                    disc_idx = convert(Array{Bool}, readdlm(joinpath(disc_dir, fn), ',')[2:end, idx])
+                    d = d[map(!, disc_idx)]    
+                end
                 filter!(x->!isnan(x), d)
                 values = [pct == :avg ? Int(round(mean(d))) : Int(round(quantile(d, pct))) for pct in pcts]
                 write(f, "$scenario,", join(values, ","), "\n")
@@ -66,8 +71,9 @@ function make_percentile_tables(output_dir, gas, discount_rates, perturbation_ye
 end
 
 # helper funcion for computing std error tables at the end of a monte carlo simulation
-function make_stderror_tables(output_dir, gas, discount_rates, perturbation_years)
+function make_stderror_tables(output_dir, gas, discount_rates, perturbation_years, drop_discontinuities=false)
     scc_dir = "$output_dir/SC-$gas"     # folder with output from the MCS runs
+    drop_discontinuities ? disc_dir = joinpath(output_dir, "discontinuity_mismatch/") : nothing
     tables = "$output_dir/Tables/Std Errors"   # folder to save the tables to
     mkpath(tables)
 
@@ -79,7 +85,11 @@ function make_stderror_tables(output_dir, gas, discount_rates, perturbation_year
             write(f, "Scenario,Mean,SE\n")
             for fn in filter(x -> endswith(x, "$dr.csv"), results)  # Get the results files for this discount rate
                 scenario = split(fn)[1] # get the scenario name
-                d = readdlm(joinpath(scc_dir, fn), ',')[2:end, idx] # just keep 2020 values
+                d = readdlm(joinpath(scc_dir, fn), ',')[2:end, idx] 
+                if drop_discontinuities
+                    disc_idx = convert(Array{Bool}, readdlm(joinpath(disc_dir, fn), ',')[2:end, idx])
+                    d = d[map(!, disc_idx)]    
+                end
                 filter!(x->!isnan(x), d)
                 write(f, "$scenario, $(round(mean(d), digits=2)), $(round(sem(d), digits=2)) \n")
             end 
@@ -88,10 +98,11 @@ function make_stderror_tables(output_dir, gas, discount_rates, perturbation_year
     nothing  
 end
 
-# helper function for computing a summary table. Reports average values for all discount rates and years, and high impact value (95th pct) for 3%.
-function make_summary_table(output_dir, gas, discount_rates, perturbation_years)
+# helper function for computing a summary table. Reports average values for all discount rates and years, and high impact values (95th pct)
+function make_summary_table(output_dir, gas, discount_rates, perturbation_years, drop_discontinuities=false)
 
     scc_dir = "$output_dir/SC-$gas"     # folder with output from the MCS runs
+    drop_discontinuities ? disc_dir = joinpath(output_dir, "discontinuity_mismatch/") : nothing
     tables = "$output_dir/Tables"   # folder to save the table to
     mkpath(tables)
 
@@ -102,12 +113,17 @@ function make_summary_table(output_dir, gas, discount_rates, perturbation_years)
     data[2:end, 1] = perturbation_years
 
     for (j, dr) in enumerate(discount_rates)
-        vals = Matrix{Float64}(undef, 0, length(perturbation_years))
+        vals = Matrix{Union{Missing, Float64}}(undef, 0, length(perturbation_years))
         for scenario in scenarios
-            vals = vcat(vals, readdlm(joinpath(scc_dir, "$(string(scenario)) $dr.csv"), ',')[2:end, :])
+            curr_vals = convert(Array{Union{Missing, Float64}}, readdlm(joinpath(scc_dir, "$(string(scenario)) $dr.csv"), ',')[2:end, :])
+            if drop_discontinuities
+                disc_idx = convert(Array{Bool}, readdlm(joinpath(disc_dir, "$(string(scenario)) $dr.csv"), ',')[2:end, :])
+                curr_vals[disc_idx] .= missing
+            end
+            vals = vcat(vals, curr_vals)
         end
-        data[2:end, j+1] = mean(vals, dims=1)[:]
-        data[2:end, j+1+length(discount_rates)] = [quantile(vals[2:end, y], .95) for y in 1:length(perturbation_years)]
+        data[2:end, j+1] = mapslices(x -> mean(skipmissing(x)), vals, dims=1)[:]
+        data[2:end, j+1+length(discount_rates)] = [quantile(skipmissing(vals[2:end, y]), .95) for y in 1:length(perturbation_years)]
     end
 
     table = joinpath(tables, "Summary Table.csv")

src/montecarlo/PAGE_mcs.jl

@@ -206,7 +206,7 @@ function page_post_trial_func(mcs::SimulationInstance, trialnum::Int, ntimesteps
     base, marginal = mcs.models 
 
     # Unpack the payload object 
-    discount_rates, discount_factors, gas, perturbation_years, SCC_values, SCC_values_domestic, md_values = Mimi.payload(mcs)
+    discount_rates, discount_factors, discontinuity_mismatch, gas, perturbation_years, SCC_values, SCC_values_domestic, md_values = Mimi.payload(mcs)
 
     DF = discount_factors[rate_num]
     td_base = base[:EquityWeighting, :td_totaldiscountedimpacts]
@@ -222,6 +222,9 @@ function page_post_trial_func(mcs::SimulationInstance, trialnum::Int, ntimesteps
         perturb_marginal_page_emissions!(base, marginal, gas, pyear)
         run(marginal)
 
+        # Stores `true` if the base and marginal models trigger the discontinuity damages in different timesteps (0 otherwise)
+        discontinuity_mismatch[trialnum, j, scenario_num, rate_num] = base[:Discontinuity, :occurdis_occurrencedummy] != marginal[:Discontinuity, :occurdis_occurrencedummy]
+
         td_marginal = marginal[:EquityWeighting, :td_totaldiscountedimpacts]   
         pulse_size = gas == :CO2 ? 100_000 : 1          
         scc = ((td_marginal / UDFT_base[idx]) - (td_base / UDFT_base[idx])) / pulse_size * page_inflator

src/montecarlo/run_scc_mcs.jl

@@ -1,14 +1,15 @@
 """
     run_scc_mcs(model::model_choice; 
-        gas::Union{Symbol, Nothing} = nothing,
-        trials::Int = 10000,
-        perturbation_years::Vector{Int} = _default_perturbation_years,
-        discount_rates::Vector{Float64} = _default_discount_rates, 
-        domestic::Bool = false,
-        output_dir::String = nothing, 
-        save_trials::Bool = false,
-        save_md::Bool = false,
-        tables::Bool = true)
+                gas::Union{Symbol, Nothing} = nothing,
+                trials::Int = 10000,
+                perturbation_years::Vector{Int} = _default_perturbation_years,
+                discount_rates::Vector{Float64} = _default_discount_rates, 
+                domestic::Bool = false,
+                output_dir::Union{String, Nothing} = nothing, 
+                save_trials::Bool = false,
+                tables::Bool = true,
+                drop_discontinuities::Bool = false,
+                save_md::Bool = false)
 
 Run the Monte Carlo simulation used by the IWG for calculating a distribution of SCC values for the 
 Mimi model `model_choice` and the specified number of trials `trials`. The SCC is calculated for all 
@@ -22,9 +23,16 @@ equals `true`, then SCC values will also be calculated using only domestic damag
 Output files will be saved in the `output_dir`. If none is provided, it will default to "./output/". 
 A new sub directory will be created each time this function is called, with the following name: "yyyy-mm-dd HH-MM-SS MODEL SC-\$gas MC\$trials".
 
+Several keyword arguments allow for the following:
 If `tables` equals `true`, then a set of summary statistics tables will also be saved in the output folder.
-If `save_trials` equals `true`, then a file with all of the sampled input trial data will also be saved in
-the output folder. If `save_md` equals `true`, then global undiscounted marginal damages from each run of 
+
+If `save_trials` equals `true`, then a file with all of the sampled input trial data will also be saved in the output folder. 
+
+If `drop_discontinuities` equals `true`, then outliers from the PAGE model (runs where discontinuity damages are triggered
+in different timesteps in the base and perturbed models) will not contribute to summary statistics. An additional folder
+"discontinuity_mismatch" contains files identifying in which runs the discrepencies occured.
+
+If `save_md` equals `true`, then global undiscounted marginal damages from each run of 
 the simulation will be saved in a subdirectory "output/marginal_damages".
 """
 function run_scc_mcs(model::model_choice; 
@@ -35,8 +43,9 @@ function run_scc_mcs(model::model_choice;
     domestic::Bool = false,
     output_dir::Union{String, Nothing} = nothing, 
     save_trials::Bool = false,
-    save_md::Bool = false,
-    tables::Bool = true)
+    tables::Bool = true,
+    drop_discontinuities::Bool = false,
+    save_md::Bool = false)
 
     # Check the gas
     if gas === nothing
@@ -125,6 +134,12 @@ function run_scc_mcs(model::model_choice;
         perturbation_years = model_years[_first_idx : _last_idx]  # figure out which years of the model's time index we need to use to cover all desired perturbation years
     end
 
+    # For each run, this array will store whether there is a discrepency between the base and marginal models triggering the discontinuity damages in different timesteps
+    if model == PAGE
+        discontinuity_mismatch = Array{Bool, 4}(undef, trials, length(perturbation_years), length(scenarios), length(discount_rates))
+        push!(payload, discontinuity_mismatch)
+    end
+
     # Make an array to hold all calculated scc values
     SCC_values = Array{Float64, 4}(undef, trials, length(perturbation_years), length(scenarios), length(discount_rates))
     if domestic 
@@ -188,9 +203,34 @@ function run_scc_mcs(model::model_choice;
             SCC_values_domestic = new_domestic_values
         end
 
-        perturbation_years = all_years
+        # reset perturbation years to all user requested years, unless model is PAGE, for which this is done below
+        if model != PAGE
+            perturbation_years = all_years
+        end
     end
 
+    # Save the information about which runs have a discrepency between base/marginal models of the discontinuity damages
+    if model == PAGE
+        # access the computed saved values from the simulation instance; it's the third item in the payload object for PAGE
+        discontinuity_mismatch = Mimi.payload(sim_results)[3] 
+
+        if _need_to_interpolate
+            new_discontinuity_mismatch = Array{Bool}(undef, trials, length(all_years), length(scenarios), length(discount_rates))
+            for i in 1:trials, j in 1:length(scenarios), k in 1:length(discount_rates)
+                new_discontinuity_mismatch[i, :, j, k] = convert(Array{Bool}, _interpolate(discontinuity_mismatch[i, :, j, k], perturbation_years, all_years) .> 0)
+            end
+            discontinuity_mismatch = new_discontinuity_mismatch
+            perturbation_years = all_years 
+        end
+
+        disc_dir = joinpath(output_dir, "discontinuity_mismatch/")
+        # has the same 4-D array structure as the SCC values, so can use the same function to save them to files
+        write_scc_values(discontinuity_mismatch, disc_dir, perturbation_years, discount_rates)
+
+        disc_sum = dropdims(sum(discontinuity_mismatch, dims=(1,3)), dims=(1,3)) # summary table of how many occured (rows are perturbation years, columns are discount rates)
+        writedlm(joinpath(disc_dir, "discontinuity_summary.csv"), disc_sum, ',') 
+    end
+
     # Save the SCC values
     scc_dir = joinpath(output_dir, "SC-$gas/")
     write_scc_values(SCC_values, scc_dir, perturbation_years, discount_rates)
@@ -201,9 +241,9 @@ function run_scc_mcs(model::model_choice;
 
     # Build the stats tables
     if tables
-        make_percentile_tables(output_dir, gas, discount_rates, perturbation_years)
-        make_stderror_tables(output_dir, gas, discount_rates, perturbation_years)
-        make_summary_table(output_dir, gas, discount_rates, perturbation_years)
+        make_percentile_tables(output_dir, gas, discount_rates, perturbation_years, drop_discontinuities)
+        make_stderror_tables(output_dir, gas, discount_rates, perturbation_years, drop_discontinuities)
+        make_summary_table(output_dir, gas, discount_rates, perturbation_years, drop_discontinuities)
     end
 
     nothing

test/test_PAGE.jl

@@ -19,6 +19,11 @@ using DelimitedFiles
     MimiIWG.run_scc_mcs(PAGE, trials=2, output_dir = tmp_dir, domestic=true)
     rm(tmp_dir, recursive=true)
 
+    # test the drop discontinuities flag set to `true`
+    tmp_dir = joinpath(@__DIR__, "tmp")
+    MimiIWG.run_scc_mcs(PAGE, trials=2, output_dir = tmp_dir, domestic=true, drop_discontinuities = true)
+    rm(tmp_dir, recursive=true)
+
 end
 
 @testset "Deterministic SCC validation" begin