Format files using DocumentFormat

docs/make.jl

@@ -1,15 +1,15 @@
 using Documenter, MimiMooreEtAlAgricultureImpacts
 
 makedocs(
-    modules = [MimiMooreEtAlAgricultureImpacts],
-    sitename = "Moore et al. Agriculture Documentation",
-    pages = [
+    modules=[MimiMooreEtAlAgricultureImpacts],
+    sitename="Moore et al. Agriculture Documentation",
+    pages=[
         "Home" => "index.md",
         "Reference" => "reference.md"
     ],
-    format = Documenter.HTML(prettyurls = get(ENV, "JULIA_NO_LOCAL_PRETTY_URLS", nothing) === nothing)
+    format=Documenter.HTML(prettyurls=get(ENV, "JULIA_NO_LOCAL_PRETTY_URLS", nothing) === nothing)
 )
 
 deploydocs(
-    repo = "github.com/rffscghg/MimiMooreEtAlAgricultureImpacts.jl.git",
+    repo="github.com/rffscghg/MimiMooreEtAlAgricultureImpacts.jl.git",
 )

examples/main_scc.jl

@@ -18,4 +18,4 @@ for gtap in MimiMooreEtAlAgricultureImpacts.gtaps
     write(f, "$gtap,$ag_scc\n")
 end
 
-close(f)
+close(f)

src/MimiMooreEtAlAgricultureImpacts.jl

@@ -1,4 +1,4 @@
-module MimiMooreEtAlAgricultureImpacts 
+module MimiMooreEtAlAgricultureImpacts
 
 using Mimi
 using Interpolations
@@ -13,4 +13,4 @@ include("core/get_model.jl")
 include("core/get_ag_scc.jl")
 include("core/mcs.jl")
 
-end
+end

src/core/AgricultureComponent.jl

@@ -1,4 +1,4 @@
-using Interpolations
+using Interpolations
 using Mimi
 
 # Moore et al Agriculture component (with linear interpolation between gtap temperature points)
@@ -7,26 +7,26 @@ using Mimi
     fund_regions = Index()
 
     gdp90 = Parameter(index=[fund_regions])
-    income = Parameter(index=[time,fund_regions])
+    income = Parameter(index=[time, fund_regions])
     pop90 = Parameter(index=[fund_regions], unit="million")
-    population = Parameter(index=[time,fund_regions], unit="million")
+    population = Parameter(index=[time, fund_regions], unit="million")
 
-    agrish = Variable(index=[time,fund_regions])     # agricultural share of the economy
+    agrish = Variable(index=[time, fund_regions])     # agricultural share of the economy
     agrish0 = Parameter(index=[fund_regions])        # initial share 
-    agel = Parameter(default = 0.31)            # elasticity
+    agel = Parameter(default=0.31)            # elasticity
 
-    agcost = Variable(index=[time,fund_regions])     # This is the main damage variable
+    agcost = Variable(index=[time, fund_regions])     # This is the main damage variable
 
     temp = Parameter(index=[time], unit="degC")              # Moore et al uses global temperature (original FUND ImpactAgriculture component uses regional temperature)
 
     # Moore additions:
 
-    AgLossGTAP = Variable(index=[time,fund_regions]) # Moore's fractional loss (intermediate variable for calculating agcost)
+    AgLossGTAP = Variable(index=[time, fund_regions]) # Moore's fractional loss (intermediate variable for calculating agcost)
 
     gtap_df = Parameter(index=[fund_regions, 3])  # three temperature data points per region
 
-    floor_on_damages = Parameter{Bool}(default = true)
-    ceiling_on_benefits = Parameter{Bool}(default = false)
+    floor_on_damages = Parameter{Bool}(default=true)
+    ceiling_on_benefits = Parameter{Bool}(default=false)
 
     function run_timestep(p, v, d, t)
         for r in d.fund_regions
@@ -38,8 +38,8 @@ using Mimi
             # Interpolate for p.temp, using the three gtap welfare points with the additional origin (0,0) point
             impact = linear_interpolate([0, p.gtap_df[r, :]...], collect(0:3), p.temp[t])
             impact = p.floor_on_damages ? max(-100, impact) : impact
-            impact = p.ceiling_on_benefits ? min(100, impact)  : impact
-            v.AgLossGTAP[t, r] = - impact / 100 # We take the negative to go from impact to loss
+            impact = p.ceiling_on_benefits ? min(100, impact) : impact
+            v.AgLossGTAP[t, r] = -impact / 100 # We take the negative to go from impact to loss
 
             # Calculate total cost for the ag sector based on the percent loss
             v.agcost[t, r] = p.income[t, r] * v.agrish[t, r] * v.AgLossGTAP[t, r]

src/core/get_ag_scc.jl

@@ -17,35 +17,35 @@ allowed to exceed 100% of the size of the agricultural sector in each region.
 If `ceiling_on_benefits` = true, then the agricultural benefits in each timestep will not
 be allowed to exceed 100% of the size of the agricultural sector in each region.
 """
-function get_ag_scc(gtap::String; 
-    prtp::Float64 = 0.03, 
-    horizon::Int = _default_horizon,
-    floor_on_damages::Bool = true,
-    ceiling_on_benefits::Bool = false)
+function get_ag_scc(gtap::String;
+    prtp::Float64=0.03,
+    horizon::Int=_default_horizon,
+    floor_on_damages::Bool=true,
+    ceiling_on_benefits::Bool=false)
 
     horizon in years ? nothing : error("Invalid value: $horizon for `horizon`, must be within the model years.")
 
     # Run base model
-    base_m = get_model(gtap, floor_on_damages = floor_on_damages, ceiling_on_benefits = ceiling_on_benefits)
+    base_m = get_model(gtap, floor_on_damages=floor_on_damages, ceiling_on_benefits=ceiling_on_benefits)
     run(base_m)
 
     # Run model with pulse in 2020
-    pulse_m = get_model(gtap, pulse=true, floor_on_damages = floor_on_damages, ceiling_on_benefits = ceiling_on_benefits)
+    pulse_m = get_model(gtap, pulse=true, floor_on_damages=floor_on_damages, ceiling_on_benefits=ceiling_on_benefits)
     run(pulse_m)
 
     # calculate SCC 
     base_damages = dropdims(sum(base_m[:Agriculture, :agcost], dims=2), dims=2)
     pulse_damages = dropdims(sum(pulse_m[:Agriculture, :agcost], dims=2), dims=2)
-    marginal_damages = (pulse_damages - base_damages) * 10^9 / 10^9 * 12/44  # 10^9 for billions of dollars; /10^9 for Gt pulse; 12/44 to go from $/ton C to $/ton CO2
+    marginal_damages = (pulse_damages - base_damages) * 10^9 / 10^9 * 12 / 44  # 10^9 for billions of dollars; /10^9 for Gt pulse; 12/44 to go from $/ton C to $/ton CO2
 
     start_idx = findfirst(isequal(pulse_year), years)
     end_idx = findfirst(isequal(horizon), years)
 
     # Implement discounting as a 10-year step function as described by Delevane
-    discount_factor = [(1 + prtp) ^ (-1 * t * 10) for t in 0:end_idx-start_idx]
+    discount_factor = [(1 + prtp)^(-1 * t * 10) for t in 0:end_idx-start_idx]
     npv = marginal_damages[start_idx:end_idx] .* 10 .* discount_factor  # multiply by 10 so that value of damages is used for all 10 years
 
-    ag_scc = sum(npv) 
+    ag_scc = sum(npv)
 
     return ag_scc
-end
+end

src/core/get_model.jl

@@ -22,29 +22,29 @@ to exceed 100% of the size of the agricultural sector in each region.
 If `ceiling_on_benefits` = true, then the agricultural benefits in each timestep will not be
 allowed to exceed 100% of the size of the agricultural sector in each region.
 """
-function get_model( gtap::String; 
-                    pulse::Bool=false,
-                    floor_on_damages::Bool = true,
-                    ceiling_on_benefits::Bool = false
-                )
+function get_model(gtap::String;
+    pulse::Bool=false,
+    floor_on_damages::Bool=true,
+    ceiling_on_benefits::Bool=false
+)
 
     gtap in gtaps ? nothing : error("Unknown GTAP dataframe specification: \"$gtap\". Must be one of the following: $gtaps")
 
     # Read in the USG2 socioeconomics data 
-    usg2_population = Array{Float64, 2}(readdlm(joinpath(fund_datadir, "usg2_population.csv"),',')[2:end, 2:end])   # Saved from SCCinputs.rdata from Delavane
-    usg2_income = Array{Float64, 2}(readdlm(joinpath(fund_datadir, "usg2_income.csv"),',')[2:end, 2:end])   # Saved from SCCinputs.rdata from Delavane
-    
+    usg2_population = Array{Float64,2}(readdlm(joinpath(fund_datadir, "usg2_population.csv"), ',')[2:end, 2:end])   # Saved from SCCinputs.rdata from Delavane
+    usg2_income = Array{Float64,2}(readdlm(joinpath(fund_datadir, "usg2_income.csv"), ',')[2:end, 2:end])   # Saved from SCCinputs.rdata from Delavane
+
     # Read in DICE temperature pathway
     dice_temp_file = pulse ? "dice_temp_pulse.csv" : "dice_temp.csv"
-    dice_temp = readdlm(joinpath(dice_datadir, dice_temp_file), Float64)[:]      
-
-    params = Dict{Tuple, Any}([
-        (:Agriculture, :population) =>  usg2_population[2:end, :],     # 2000:10:2300
-        (:Agriculture, :income) =>      usg2_income[2:end, :],         # 2000:10:2300
-        (:Agriculture, :pop90) =>       usg2_population[1, :],         # 1990 is the first row
-        (:Agriculture, :gdp90) =>       usg2_income[1, :],             # 1990 is the first row
-        (:Agriculture, :temp) =>        dice_temp,
-        (:Agriculture, :agrish0) =>     Array{Float64, 1}(readdlm(joinpath(fund_datadir, "agrish0.csv"), ',', skipstart=1)[:,2]),
+    dice_temp = readdlm(joinpath(dice_datadir, dice_temp_file), Float64)[:]
+
+    params = Dict{Tuple,Any}([
+        (:Agriculture, :population) => usg2_population[2:end, :],     # 2000:10:2300
+        (:Agriculture, :income) => usg2_income[2:end, :],         # 2000:10:2300
+        (:Agriculture, :pop90) => usg2_population[1, :],         # 1990 is the first row
+        (:Agriculture, :gdp90) => usg2_income[1, :],             # 1990 is the first row
+        (:Agriculture, :temp) => dice_temp,
+        (:Agriculture, :agrish0) => Array{Float64,1}(readdlm(joinpath(fund_datadir, "agrish0.csv"), ',', skipstart=1)[:, 2]),
         (:Agriculture, :gtap_df_all) => gtap_df_all
     ])
 

src/core/mcs.jl

@@ -5,14 +5,14 @@ function get_probdists_gtap_df(n=1000)
 
     # For each region and temperature point we construct an interpolation where the x values are between 0 and 1
     # and the y values are the values from the three scenarios.
-    dists = [LinearInterpolation([0.,0.5,1.], [lowDF[r,temp],midDF[r,temp], highDF[r,temp]]) for r in 1:16, temp in 1:3]
+    dists = [LinearInterpolation([0., 0.5, 1.], [lowDF[r, temp], midDF[r, temp], highDF[r, temp]]) for r in 1:16, temp in 1:3]
 
     # We only sample one set of random numbers, as we want perfect correlation between all the individual
     # parameter values.
     samples = rand(TriangularDist(0., 1., 0.5), n)
 
     # Now evaluate the interpolated function we created above with the samples from the triangular distributions
-    sample_stores = [Mimi.SampleStore(dists[r,temp].(samples)) for r in 1:16, temp in 1:3]
+    sample_stores = [Mimi.SampleStore(dists[r, temp].(samples)) for r in 1:16, temp in 1:3]
 
     return sample_stores
 end

src/core/utils.jl

@@ -13,9 +13,9 @@ const fund_datadir = joinpath(@__DIR__, "../../data/FUND params")
 const dice_datadir = joinpath(@__DIR__, "../../data/DICE climate output")
 
 # Moore et al uses regions in alphabetical order; need to be conscious of switching the regional ordering for running with FUND parameters
-alpha_order = ["ANZ","CAM","CAN","CEE","CHI","FSU","JPK","MDE","NAF","SAM","SAS","SEA","SIS","SSA","USA","WEU"]
+alpha_order = ["ANZ", "CAM", "CAN", "CEE", "CHI", "FSU", "JPK", "MDE", "NAF", "SAM", "SAS", "SEA", "SIS", "SSA", "USA", "WEU"]
 alpha_order[[4, 9, 10]] = ["EEU", "MAF", "LAM"]     # three regions named slightly different: CEU, NAF, CAM --> EEU, MAF, LAM
-const fund_regions = ["USA","CAN","WEU","JPK","ANZ","EEU","FSU","MDE","CAM","LAM","SAS","SEA","CHI","MAF","SSA","SIS"]
+const fund_regions = ["USA", "CAN", "WEU", "JPK", "ANZ", "EEU", "FSU", "MDE", "CAM", "LAM", "SAS", "SEA", "CHI", "MAF", "SSA", "SIS"]
 const switch_region_indices = [findfirst(isequal(region), alpha_order) for region in fund_regions]
 
 # Returns the Moore gtap data points (16 regions x 3 points) in the FUND regional order
@@ -24,17 +24,17 @@ function get_gtap_df(gtap::String)
     gtap_dir = joinpath(@__DIR__, "../../data/GTAP DFs")   # The five welfare dataframes from Fran Moore are in this folder
     gtap_data = Array(readdlm(joinpath(gtap_dir, "$gtap.csv"), ',', skipstart=1)')
     return gtap_data[switch_region_indices, :]
-end 
+end
 
 const gtap_df_all = reshape(reduce(hcat, [get_gtap_df(gtap) for gtap in gtaps]), (16, 3, 5))
 
 # helper function for linear interpolation
-function linear_interpolate(values::AbstractArray, original_domain::AbstractArray, new_domain::Union{AbstractArray, Number})
+function linear_interpolate(values::AbstractArray, original_domain::AbstractArray, new_domain::Union{AbstractArray,Number})
     # Build the interpolation object with linear interpolation between the provided points, and extrapolation beyond the points
-    itp = extrapolate(interpolate((original_domain,), values, Gridded(Linear())), Line())    
+    itp = extrapolate(interpolate((original_domain,), values, Gridded(Linear())), Line())
 
     # Get the interpolated values for the point(s) in new_domain
-    if new_domain isa Number 
+    if new_domain isa Number
         return itp(convert(Float64, new_domain)) # itp(x) returns a Ratios.SimpleRatio is x is just a single number, needs to be converted to a Float
     elseif new_domain isa Array
         return itp(new_domain)   # itp([x1, x2, etc]) returns an Array
@@ -43,6 +43,6 @@ end
 
 # TODO: implement quadratic interpolation correctly
 # function quadratic_interpolate(values, original_domain, new_domain)
-    # itp = interpolate(values, BSpline(Quadratic(Free())), OnGrid()) # TODO: only works if original x values are 1,2,3,etc.
-    # return [itp[i] for i in new_domain]
+# itp = interpolate(values, BSpline(Quadratic(Free())), OnGrid()) # TODO: only works if original x values are 1,2,3,etc.
+# return [itp[i] for i in new_domain]
 # end

test/test_api.jl

@@ -19,4 +19,4 @@ end
 
 @testitem "Test the ceiling on benefits" begin
     @test MimiMooreEtAlAgricultureImpacts.get_ag_scc("lowDF", prtp=0.03, ceiling_on_benefits=false) < MimiMooreEtAlAgricultureImpacts.get_ag_scc("lowDF", prtp=0.03, ceiling_on_benefits=true) # ceiling on benefits is only binding in the "lowDF" scenario
-end
+end

test/test_validation.jl

@@ -6,12 +6,12 @@
     results = readdlm(joinpath(@__DIR__, "../data/validation/ag_scc.csv"), ',')
     mimi_sccs = Vector{Any}()
     global i = 2
-    for gtap in MimiMooreEtAlAgricultureImpacts.gtaps 
+    for gtap in MimiMooreEtAlAgricultureImpacts.gtaps
         for dr in [0.025, 0.03, 0.05]
             println(gtap, dr)
-            mimi_scc = MimiMooreEtAlAgricultureImpacts.get_ag_scc(gtap, prtp = dr)
+            mimi_scc = MimiMooreEtAlAgricultureImpacts.get_ag_scc(gtap, prtp=dr)
             r_scc = results[i, 3]
-            @test mimi_scc ≈ r_scc atol=0.034
+            @test mimi_scc ≈ r_scc atol = 0.034
             push!(mimi_sccs, mimi_scc)
             println(mimi_scc)
             println(r_scc)
@@ -21,4 +21,4 @@
 
     # comparison = hcat(results, ["Mimi_scc", mimi_sccs...])
     # writedlm(joinpath(@__DIR__, "scc_comparison.csv"), comparison, ",")
-end
+end