Merge pull request #68 from fund-model/SC-N2O

add gas weight conversion for N2O to N and fix weight conversion for SC-SF6

src/new_marginaldamages.jl

@@ -304,6 +304,33 @@ end
     end
 end
 
+function _gas_normalization(gas::Symbol)
+    if gas == :CO2 
+        return 12/44    # convert from tons CO2 to tons C
+    elseif gas == :CH4
+        return 1
+    elseif gas == :N2O
+        return 28/44    # convert from tons N2O to tons N
+    elseif gas == :SF6 
+        return 1
+    else
+        error("Unknown gas :$gas.")
+    end
+end
+
+function _weight_normalization(gas::Symbol)
+    if gas == :CO2 
+        return 1e6   # convert from tonnes to Mt since component expects Mt
+    elseif gas == :CH4
+        return 1e6   # convert from tonnes to Mt since component expects Mt
+    elseif gas == :N2O
+        return 1e6   # convert from tonnes to Mt since component expects Mt
+    elseif gas == :SF6 
+        return 1e3   # convert from tonnes to kt since component expects Mt
+    else
+        error("Unknown gas :$gas.")
+    end
+end 
 """
 Adds an emissionspulse component to `m`, and sets the additional emissions if a year is specified.
 The size of the marginal emission pulse can be modified with the `pulse_size` keyword argument, in metric 
@@ -317,11 +344,10 @@ function add_marginal_emissions!(m, year::Union{Int, Nothing} = nothing; gas::Sy
     nyears = length(Mimi.time_labels(m))
     addem = zeros(nyears) 
     if year != nothing 
-        # pulse is spread over ten years, and emissions components is in Mt so 
-        # divide by 1e7, and convert from CO2 to C if gas==:CO2 because emissions 
-        # component is in MtC; also use 12/44 to convert from GtCO2 to GtC for CO2 since the 
-        # emissions component expects GtC not GtCO2 but pulse is in CO2
-        addem[getindexfromyear(year):getindexfromyear(year) + 9] .= pulse_size / 1e7 * (gas == :CO2 ? 12/44 : 1)
+        # need to (1) divide by pulse spread over 10 years and weight adjustment 
+        # from tons to expected units and then (2) normalize from entered gas 
+        # units ie. CO2 to expected units ie. C
+        addem[getindexfromyear(year):getindexfromyear(year) + 9] .= pulse_size / (10 * _weight_normalization(gas)) *  _gas_normalization(gas)
     end
     set_param!(m, :emissionspulse, :add, addem)
 
@@ -354,7 +380,10 @@ function perturb_marginal_emissions!(m::Model, year; comp_name::Symbol = :emissi
 
     nyears = length(Mimi.dimension(m, :time))
     new_em = zeros(nyears)
-    new_em[getindexfromyear(year):getindexfromyear(year) + 9] .= pulse_size / 1e7 * (gas == :CO2 ? 12/44 : 1)
+    # need to (1) divide by pulse spread over 10 years and weight adjustment 
+    # from tons to expected units and then  (2) normalize from entered gas 
+    # units ie. CO2 to expected units ie. C
+    new_em[getindexfromyear(year):getindexfromyear(year) + 9] .= pulse_size / (10 * _weight_normalization(gas)) * _gas_normalization(gas)
     emissions[:] = new_em
 
 end

test/create_scc_validation_file.jl

@@ -0,0 +1,67 @@
+# This file can be used to create the validation file for scc files. It will create
+# a validation file containing all possibilities of parameter values defined in the 
+# specs dictionary below produce the same results. A subset of these values are 
+# tested in the "runtests.jl" file deployed by Travis, and the full set can be tested
+# manually using the scc_validation_full.jl file
+
+using MimiFUND
+using DataFrames
+using Query
+using CSVFiles
+using Test
+
+specs = Dict([
+    :gas => [:CO2, :CH4, :N2O, :SF6],
+    :year => [2020, 2055],
+    :eta => [0, 1.45],
+    :prtp => [0.015, 0.03],
+    :equity_weights => [true, false],
+    :equity_weights_normalization_region => [0, 1, 10],
+    :last_year => [2300, 3000],
+    :pulse_size => [1, 1e7]
+])
+
+results = DataFrame(gas = [], year = [], eta = [], prtp = [], equity_weights = [], equity_weights_normalization_region = [], last_year = [], pulse_size = [], SC = [])
+
+for gas in specs[:gas]
+    for year in specs[:year]
+        for eta in specs[:eta]
+            for prtp in specs[:prtp]
+                for equity_weights in specs[:equity_weights]
+                    for equity_weights_normalization_region in specs[:equity_weights_normalization_region]
+                        for last_year in specs[:last_year]
+                            for pulse_size in specs[:pulse_size]
+                                sc = MimiFUND.compute_sc(gas=gas, year=year, eta=eta, prtp=prtp, equity_weights=equity_weights, equity_weights_normalization_region=equity_weights_normalization_region, last_year=last_year, pulse_size=pulse_size)
+                                push!(results, (gas, year, eta, prtp, equity_weights, equity_weights_normalization_region, last_year, pulse_size, sc))
+                            end
+                        end
+                    end
+                end
+            end
+        end
+    end
+end
+
+path = joinpath(@__DIR__, "SC validation data", "deterministic_sc_values_v3-13-0.csv")
+save(path, results)
+
+## Compare v3.11.7 to v3.13.0 ... :CO2 and :CH4 should not change
+
+old_sc = load(joinpath(@__DIR__, "SC validation data", "deterministic_sc_values_v3-11-7.csv")) |> DataFrame 
+new_sc = load(joinpath(@__DIR__, "SC validation data", "deterministic_sc_values_v3-13-0.csv")) |> DataFrame 
+
+filter!(:gas => x -> (x == "CO2" || x == "CH4"), old_sc)
+filter!(:gas => x -> (x == "CO2" || x == "CH4"), new_sc)
+
+# after running scc_validation_all.jl on Mimi v3.12.1 on Lisa Rennels' machine,
+# we note that the differences between results and validation results pass only at
+# a tolerance of 1.0e-1 for pulse_size of 1 and 1e-9 for pulse_size of 1.0e7, so 
+# we mirror that here
+
+old_sc_bigpulse = filter!(:pulse_size => x -> x == 1.0e7, old_sc)
+new_sc_bigpulse = filter!(:pulse_size => x -> x == 1.0e7, new_sc)
+@test all(isapprox.(old_sc_bigpulse[!, :SC], new_sc_bigpulse[!, :SC], atol = 1e-9))
+
+old_sc_smallpulse = filter!(:pulse_size => x -> x == 1.0e7, old_sc)
+new_sc_smallpulse = filter!(:pulse_size => x -> x == 1.0e7, new_sc)
+@test all(isapprox.(old_sc_smallpulse[!, :SC], new_sc_smallpulse[!, :SC], atol = 1e-1))

test/runtests.jl

@@ -120,7 +120,7 @@ result = MimiFUND.compute_sc(year=2050, return_mm = true)
 scch4 = MimiFUND.compute_scch4(year = 2020)
 scn2o = MimiFUND.compute_scn2o(year = 2020)
 scsf6 = MimiFUND.compute_scsf6(year = 2020)
-@test scsf6 < scch4 < scn2o
+@test scch4 < scn2o < scsf6
 
 # Test that modifying the pulse_size keyword changes the values, but not by much
 scch4_2 = MimiFUND.compute_scch4(year = 2020, pulse_size = 1e3)
@@ -131,7 +131,7 @@ scsf6_2 = MimiFUND.compute_scsf6(year = 2020, pulse_size = 1e3)
 @test scn2o != scn2o_2
 @test scn2o ≈ scn2o_2 rtol = 1e-3
 @test scsf6 != scsf6_2
-@test scsf6 ≈ scsf6_2 rtol = 1e-3
+@test scsf6 ≈ scsf6_2 rtol = 1e-1 # TODO Lisa Rennels - ok to increase tolerance here?
 
 # Test monte carlo simulation
 scco2_values = MimiFUND.compute_sc(year = 2020, gas = :CO2, n = 10)
@@ -171,9 +171,9 @@ end #test-mcs testset
 datadir = joinpath(@__DIR__, "SC validation data")
 atol = 1e-2
 
-# Test a subset of all validation configurations against the pre-saved values from MimiFUND v3.11.7
+# Test a subset of all validation configurations against the pre-saved values from MimiFUND v3.13.0
 
-validation_results = load(joinpath(datadir, "deterministic_sc_values.csv")) |> DataFrame
+validation_results = load(joinpath(datadir, "deterministic_sc_values_v3-13-0.csv")) |> DataFrame
 
 for spec in [
     (gas = :CO2, year = 2020, eta = 0., prtp = 0.03, equity_weights = true, equity_weights_normalization_region = 1, last_year = 3000, pulse_size = 1.),
@@ -192,9 +192,10 @@ for spec in [
     @test sc ≈ validation_value[1, :SC] atol = atol
 end
 
-# Test Monte Carlo results with the same configuration and seed
+# Test Monte Carlo results with the same configuration and seed - note here we use
+# FUND 3.11.7 since 3.13.0 only updated SC-N2O and SC-SF6
 sc_mcs = MimiFUND.compute_sc(gas = :CO2, year = 2020, eta = 1.45, prtp = 0.015, n = 25, seed = 350)
-validation_mcs = load(joinpath(datadir, "mcs_sc_values.csv")) |> DataFrame
+validation_mcs = load(joinpath(datadir, "mcs_sc_values_v3-11-7.csv")) |> DataFrame
 @test all(isapprox.(sc_mcs, validation_mcs[!, :SCCO2], atol = atol))
 
 end

test/scc_validation_full.jl

@@ -2,7 +2,7 @@
 # This file can be run manually to check that Social Cost calculations using all possibilities of
 # parameter values defined in the specs dictionary below produce the same results. A subset of
 # these values are tested in the "runtests.jl" file deployed by Travis. The validation of these
-# runs depends on the pre-saved values from MimiFUND v3.11.7 in the "SC validation data" folder.
+# runs depends on the pre-saved values from MimiFUND v3.13.0 in the "SC validation data" folder.
 
 using MimiFUND
 using DataFrames
@@ -44,5 +44,5 @@ for gas in specs[:gas]
     end
 end
 
-validation_results = load(joinpath(datadir, "deterministic_sc_values.csv")) |> DataFrame
-@test all(isapprox.(results[!, :SC], validation_results[!, :SC], atol = 1e-11))
+validation_results = load(joinpath(datadir, "deterministic_sc_values_v3-13-0.csv")) |> DataFrame
+@test all(isapprox.(results[!, :SC], validation_results[!, :SC], atol = 1e-10)