text:
function haversine(p1, p2; R=6371.0)
    lon1, lat1 = p1
    lon2, lat2 = p2
    φ1 = lat1 * π/180.0
    φ2 = lat2 * π/180.0
    Δφ = (lat2-lat1) * π/180.0
    Δλ = (lon2-lon1) * π/180.0
    a = sin(Δφ/2.0)^2.0 + cos(φ1)*cos(φ2) * sin(Δλ)^2.0
    c = 2.0 * atan(sqrt(a), sqrt(1.0-a))
    return R*c
end

"""
    slidingwindow(L::LT, f::FT, d::IT) where {LT <: SimpleSDMLayer, FT <: Function, IT <: Number}

This function will replace the value at any cell by applying the function `f` to
the array of cells values that are within a distance `d` (in kilometers) from
the focal cell. This is, for example, useful to use an average to smooth out the
layers. The distance is estimated using the haversine distance, assuming that
the radius of the Earth is 6371.0 km. This means that the size of the window
will vary a little bit across latitudes, but this is far better than using a
number of cells, which would have dramatic consequences near the poles.

It *always* returns a `SimpleSDMResponse`, and the cells containing `nothing`
will also not contain a value in the output. This is *different* from the
behavior of `coarsen`, which tends to expand the area of the layer in which we
have data.

This function is currently relatively slow. Performance improvements will arrive
at some point.
"""
function slidingwindow(layer::LT, f::FT, d::IT) where {LT <: SimpleSDMLayer, FT <: Function, IT <: Number}
    return_type = typeof(f(collect(layer)[1:min(3, length(layer))]))
    newgrid = convert(Matrix{Union{Nothing,return_type}}, fill(nothing, size(layer)))
    N = SimpleSDMResponse(newgrid, layer)
    pixels = []
    for lat in latitudes(layer)
        for lon in longitudes(layer)
            if !isnothing(layer[lon,lat])
                push!(pixels, (lon, lat) => layer[lon,lat])
            end
        end
    end

    for p1 in pixels
        ok = filter(p2 -> haversine(p2.first, p1.first) < 100.0, pixels)
        val = [p2.second for p2 in ok]
        N[p1.first...] = f(val)
    end

    #internal_types = unique(typeof.(N.grid))
    #N.grid = convert(Matrix{Union{internal_types...}}, N.grid)
    N = typeof(layer) <: SimpleSDMPredictor ? convert(SimpleSDMPredictor, N) : N

    return N
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 1ca14906-eca1-11ea-23f6-472ed97d75aa
begin
	using Statistics
	using Images
	using FFTW
	using Plots
	using DSP
	using ImageFiltering
	using PlutoUI
end

# ╔═╡ 42ed52ba-ed34-11ea-26b5-05379824cbc0
md"""
# Convolutions with various kernels
"""

# ╔═╡ 4c13d558-ee15-11ea-2ed9-c5fb90d93881
kernel = Kernel.gaussian((2, 2))

# ╔═╡ 673f7ac0-ee16-11ea-35d0-cf3da430b843
sum(kernel)

# ╔═╡ 9c90feb8-ec79-11ea-2870-31be5cedff43
md"""
# Function definitions
"""

# ╔═╡ 84e6a57c-edfc-11ea-01a0-157f1df77518
function show_colored_kernel(kernel)
	to_rgb(x) = RGB(max(-x, 0), max(x, 0), 0)
	to_rgb.(kernel) / maximum(abs.(kernel))
end

# ╔═╡ 9424b46a-ee16-11ea-1819-f17ce53e9997
show_colored_kernel(kernel)

# ╔═╡ 68f2afec-eca2-11ea-0758-2f22c7afdd94
function decimate(arr, ratio=5)
	return arr[1:ratio:end, 1:ratio:end]
end

# ╔═╡ aa3b9bd6-ed35-11ea-1bdc-33861bdbd29a
function shrink_image(image, ratio=5)
	(height, width) = size(image)
	new_height = height ÷ ratio - 1
	new_width = width ÷ ratio - 1
	list = [
		mean(image[
			ratio * i:ratio * (i + 1),
			ratio * j:ratio * (j + 1),
		])
		for j in 1:new_width
		for i in 1:new_height
	]
	reshape(list, new_height, new_width)
end

# ╔═╡ 6d39fea8-ed3c-11ea-3d7c-3f62ca91ce23
begin
	large_image = load("tom_in_bowtie.jpg")
	image = shrink_image(large_image, 7)
end

# ╔═╡ 2f446dcc-ee15-11ea-0e78-931ff507b5e5
size(image)

# ╔═╡ 14d5b144-ee18-11ea-0080-c187f068c168
image

# ╔═╡ 160eb236-eca1-11ea-1dbe-47ad61cc9397
function rgb_to_float(color)
    return mean([color.r, color.g, color.b])
end

# ╔═╡ fa3c5074-eca0-11ea-2d2d-bb6bcdeb834c
function fourier_spectrum_magnitudes(img)
    grey_values = rgb_to_float.(img)
    spectrum = fftshift(fft(grey_values))
	return abs.(spectrum)
end

# ╔═╡ e40d807e-ed3a-11ea-2340-7f98bd5d04a2
function plot_1d_fourier_spectrum(img, dims=1)
	spectrum = fourier_spectrum_magnitudes(img)
	plot(centered(mean(spectrum, dims=1)[1:end]))
end

# ╔═╡ beb6b4b0-eca1-11ea-1ece-e3c9931c9c13
function heatmap_2d_fourier_spectrum(img)
	heatmap(log.(fourier_spectrum_magnitudes(img)))
end

# ╔═╡ 18045956-ee18-11ea-3e34-612133e2e39c
heatmap_2d_fourier_spectrum(image)

# ╔═╡ 58f4754e-ed31-11ea-0464-5bfccf397966
function clamp_at_boundary(M, i, j)
	return M[
		clamp(i, 1, size(M, 1)),
		clamp(j, 1, size(M, 2)),
	]
end

# ╔═╡ f28af11e-ed31-11ea-2b46-7dff147ccb48
function rolloff_boundary(M, i, j)
	if (1 ≤ i ≤ size(M, 1)) && (1 ≤ j ≤ size(M, 2))
		return M[i, j]
	else
		return 0 * M[1, 1]
	end
end

# ╔═╡ 572cf620-ecb2-11ea-0019-21666a30d9d2
function convolve(M, kernel, M_index_func=clamp_at_boundary)
    height = size(kernel, 1)
    width = size(kernel, 2)

    half_height = height ÷ 2
    half_width = width ÷ 2

    new_image = similar(M)

    # (i, j) loop over the original image
    @inbounds for i in 1:size(M, 1)
        for j in 1:size(M, 2)
            # (k, l) loop over the neighbouring pixels
			new_image[i, j] = sum([
				kernel[k, l] * M_index_func(M, i - k, j - l)
				for k in -half_height:-half_height + height - 1
				for l in -half_width:-half_width + width - 1
			])
        end
    end

    return new_image
end

# ╔═╡ 5afed4ea-ee18-11ea-1aa4-abca154b3793
conv_image = convolve(image, kernel)

# ╔═╡ 6340c0f8-ee18-11ea-1765-45f4bc140670
heatmap_2d_fourier_spectrum(conv_image)

# ╔═╡ 587092e4-ecb2-11ea-18fc-ad5e9778fb30
box_blur(n) = centered(ones(n, n) ./ (n^2))

# ╔═╡ 991cb9b8-ecb8-11ea-3f80-5d95b2200259
function gauss_blur(n, sigma=0.25)
	kern = gaussian((n, n), sigma)
	return kern / sum(kern)
end

# ╔═╡ Cell order:
# ╟─42ed52ba-ed34-11ea-26b5-05379824cbc0
# ╠═6d39fea8-ed3c-11ea-3d7c-3f62ca91ce23
# ╠═2f446dcc-ee15-11ea-0e78-931ff507b5e5
# ╠═4c13d558-ee15-11ea-2ed9-c5fb90d93881
# ╠═9424b46a-ee16-11ea-1819-f17ce53e9997
# ╠═673f7ac0-ee16-11ea-35d0-cf3da430b843
# ╠═14d5b144-ee18-11ea-0080-c187f068c168
# ╠═18045956-ee18-11ea-3e34-612133e2e39c
# ╠═5afed4ea-ee18-11ea-1aa4-abca154b3793
# ╠═6340c0f8-ee18-11ea-1765-45f4bc140670
# ╟─9c90feb8-ec79-11ea-2870-31be5cedff43
# ╟─1ca14906-eca1-11ea-23f6-472ed97d75aa
# ╟─84e6a57c-edfc-11ea-01a0-157f1df77518
# ╟─68f2afec-eca2-11ea-0758-2f22c7afdd94
# ╟─aa3b9bd6-ed35-11ea-1bdc-33861bdbd29a
# ╟─160eb236-eca1-11ea-1dbe-47ad61cc9397
# ╟─fa3c5074-eca0-11ea-2d2d-bb6bcdeb834c
# ╟─e40d807e-ed3a-11ea-2340-7f98bd5d04a2
# ╟─beb6b4b0-eca1-11ea-1ece-e3c9931c9c13
# ╟─58f4754e-ed31-11ea-0464-5bfccf397966
# ╟─f28af11e-ed31-11ea-2b46-7dff147ccb48
# ╟─572cf620-ecb2-11ea-0019-21666a30d9d2
# ╟─587092e4-ecb2-11ea-18fc-ad5e9778fb30
# ╟─991cb9b8-ecb8-11ea-3f80-5d95b2200259

ext: jl

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
using Juriba
using Documenter

DocMeta.setdocmeta!(Juriba, :DocTestSetup, :(using Juriba); recursive=true)

makedocs(;
    modules=[Juriba],
    authors="Caspar <c.w.bruenech@pm.me> and contributors",
    repo="https://github.com/casparwb/Juriba.jl/blob/{commit}{path}#{line}",
    sitename="Juriba.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://casparwb.github.io/Juriba.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/casparwb/Juriba.jl",
)


text:
module ShaleDrillingLikelihood_DynamicDrillModel_ComparativeStatics

DOBTIME = false
const DOPROFILE = false

using ShaleDrillingLikelihood
using Test
using SparseArrays
using BenchmarkTools
using Calculus
using Random
using InteractiveUtils
using Distributions
using Dates
using StatsBase
using Optim
using CountPlus
using LinearAlgebra
# using Profile
# using PProf
# using Juno
# if DOPROFILE
#     using ProfileView
# end

using Optim: minimizer, minimum, maximize, maximizer, maximum

using Base.Iterators: product, OneTo

using ShaleDrillingLikelihood: ObservationDrill,
    SimulationDraw,
    idx_ρ,
    actionspace,
    ichars_sample,
    j1chars,
    tptr,
    _x, _y,
    zchars,
    ztrange,
    InitialDrilling,
    _i,
    j1_sample,
    j1_range,
    jtstart,
    exploratory_terminal,
    ssprime,
    total_leases,
    zero_out_small_probs,
    lrvar,
    lrstd,
    lrmean,
    simulate,
    theta_royalty_ρ,
    revenue,
    cost,
    extend,
    theta_ρ,
    vw_revenue,
    vw_cost,
    vw_extend,
    DataDynamicDrill,
    max_states,
    _D,
    simloglik!,
    simloglik_drill_data!,
    SimulationDraws,
    DCDPTmpVars,
    _model,
    theta_drill_ρ,
    update!,
    value_function,
    EV, dEV,
    _x, _y,
    end_ex0, end_ex1, end_lrn, end_inf,
    state_idx,
    j2ptr, tptr

println("print to keep from blowing up")

SEED = 1234

@testset "simulate all data" begin

    Random.seed!(SEED)

    discount = ((0x1.006b55c832502p+0)^12 / 1.125) ^ (1/4)  # real discount rate

    # set up coefs
    θρ = 0.0
    αψ = 0.33
    αg = 0.56

    # set up coefs
    #            θρ  ψ        g    x_r     κ₁   κ₂
    θ_royalty = [θρ, 1.0,    0.1, -0.3,  -0.6, 0.6]
    #                αψ, αg, γx   σ2η, σ2u   # η is iwt, u is iw
    θ_produce = vcat(αψ, αg, 0.2, 0.3, 0.4)
    #            drill  ext    α0  αg  αψ  θρ

    # model
    rwrd_c = DrillReward(DrillingRevenue(Constrained(;log_ogip=αg, α_ψ=αψ),NoTrend(),NoTaxes()), DrillingCost_constant(), ExtensionCost_Constant())
    rwrd_u = UnconstrainedProblem(rwrd_c)
    rwrd = rwrd_u

    # parms
    num_i = 250

    # grid sizes
    nψ =  13
    nz = 15

    # simulations
    M = 250

    # observations
    num_zt = 150          # drilling
    obs_per_well = 10:20  # pdxn

    # roylaty
    royalty_rates = [1/8, 3/16, 1/4,]
    num_royalty_rates = length(royalty_rates)
    nk_royalty = length(θ_royalty) - (length(royalty_rates)-1)-2

    # geology
    ogip_dist = Normal(4.68,0.31)

    # price parameters
    zrho = 0.8                # rho
    zmean = 1.33*(1-zrho)     # mu
    zvar = 0.265^2*(1-zrho^2) # sigsq

    # state space, psi-space
    # wp = LeasedProblem(_dmax,4,5,3,2)
    # wp = LeasedProblem(8, 8, 30, 20, 8)
    wp = LeasedProblem(8, 8, 12, 12, 8)
    ψs = range(-4.5; stop=4.5, length=nψ)

    # simulate price process
    zprocess = AR1process(zmean, zrho, zvar)
    zvec = simulate(zprocess, num_zt)
    ztrng = range(Date(2003,10); step=Dates.Quarter(1), length=num_zt)
    _zchars = ExogTimeVars(tuple.(zvec), ztrng)

    # price grid2
    zmin = min(minimum(zvec), lrmean(zprocess)-3*lrstd(zprocess))
    zmax = max(maximum(zvec), lrmean(zprocess)+3*lrstd(zprocess))
    zrng = range(zmin, zmax; length=nz)
    zs = tuple(zrng)
    @test prod(length.(zs)) == nz
    ztrans = sparse(zero_out_small_probs(tauchen_1d(zprocess, zrng), 1e-4))
    @test all(isfinite.(ztrans))

    # random vars
    log_ogip = rand(ogip_dist, num_i)
    Xroyalty = vcat(log_ogip', randn(nk_royalty-1, num_i))
    u,v = randn(num_i), randn(num_i)

    # construct royalty data
    data_roy = DataRoyalty(u, v, Xroyalty, θ_royalty, num_royalty_rates)
    royrates_endog = royalty_rates[_y(data_roy)]
    royrates_exog = sample(royalty_rates, num_i)

    # model
    ddm = DynamicDrillModel(rwrd_u, discount, wp, zs, ztrans, ψs, false)
    ddm_opts = (minmaxleases=0:0, nper_initial=0:0, nper_development=60:60, tstart=1:75, xdomain=0:0)

    theta1  = [-5.5, -2.0, -2.8, αg, αψ, θρ]
    theta2  = [-4.5, -2.0, -2.8, αg, αψ, θρ]
    theta3  = [-5.5, -1.0, -2.8, αg, αψ, θρ]
    theta4  = [-5.5, -2.0, -2.5, αg, αψ, θρ]
    theta5  = [-5.5, -2.0, -2.8, αg+0.1, αψ, θρ]
    theta6  = [-5.5, -2.0, -2.8, αg, αψ+0.3, θρ]
    theta7  = [-5.5, -2.0, -2.8, αg, αψ, θρ+0.5]

    thetas = (theta1,theta2,theta3,theta4,theta5,theta6,theta7)

    @test _nparm(rwrd) == length(theta1)

    # construct ichars for drilling
    _ichars = [gr for gr in zip(log_ogip, royrates_exog)]

    function set_seed_and_data(thet)
        Random.seed!(SEED)
        return DataDrill(u, v, _zchars, _ichars, ddm, thet; ddm_opts...)
    end

    function last_state_dist(data::DataDrill)
        x = _x(data)
        _tptr = tptr(data)
        _j2ptr = j2ptr(data)
        first_tptr = _j2ptr[2]
        last_obs = _tptr[first_tptr:end] .- 1
        last_x = x[ last_obs ]
        return sort(countmap(last_x))
    end

    # using Juno
    # Juno.@enter DataDrill(u, v, _zchars, _ichars, ddm, theta1; ddm_opts...)

    datas = map(set_seed_and_data, thetas)
    @show ys = map(d -> sum(_y(d)), datas)
    xterms = map(last_state_dist, datas)

    for x in xterms
        @show x
    end

    @test ys[1] < ys[2] # cost down ⟹ drill up
    @test ys[1] > ys[3] # ext cost down ⟹ drill down
    @test ys[1] < ys[4] # α₀ up ⟹ drill up
    @test ys[1] < ys[5] # αg up ⟹ drill up b/c E(gᵢ) > 0
    @test ys[1] < ys[6] # αψ up ⟹ more exploration, more exhaustion

    expired = end_ex0(wp) + 1
    exhausted = end_inf(wp)
    @test xterms[1][expired] < xterms[3][expired] # cheap expiration ⟹ more expire

    @test xterms[1][expired] > xterms[6][expired] # αψ ↑ ⟹ LESS expired
    @test xterms[1][exhausted] < xterms[6][exhausted] # αψ ↑ ⟹ MORE exhausted

    @show xterms[1][expired]   < xterms[7][expired] # ρ ↑ ⟹ MORE expired
    @show xterms[1][exhausted] , xterms[7][exhausted] # ρ ↑ ⟹ MORE exhausted

    #TODO: test payoff(d=0) == 0 ∀ sidx ≠ end_ex1(wp)

    # # states
    # end_ex1(wp), end_ex0(wp)
    #
    # sort(countmap(_x(data_drill_n_con)))
    # # count of where we are before dirlling first
    # # frequency of choices before expiration
    # sort(countmap([x for (x,y) in zip(_x(data_drill_n_con), _y(data_drill_n_con)) if x <= end_ex0(wp) && y > 0]))
    #
    # sort(countmap(_x(data_drill_n_con)[tptr(data_drill_n_con)[j2ptr(data_drill_n_con)[2]:end].-1]))
    #
    # # frequency of choices before expiration
    # sort(countmap([y for (x,y) in zip(_x(data_drill_n_con), _y(data_drill_n_con)) if x <= end_ex0(wp)]))
    # # infill drilling choices
    # sort(countmap([y for (x,y) in zip(_x(data_drill_n_con), _y(data_drill_n_con)) if x > end_lrn(wp)]))
    # # frequency of states just before where D==2 at end_lrn(wp)+2
    # sort(countmap(_x(data_drill_n_con)[findall(x -> x==end_lrn(wp)+2, _x(data_drill_n_con)).-1]))

end

end # module


text:
# dualtensormap.jl

immutable DualTensorMap{S,P,T} <: AbstractTensorMap{S,P,T}
    map::AbstractTensorMap{S,P,T}
end

# transposition behavior of AbstractTensorMap objects
Base.transpose(A::AbstractTensorMap) = DualTensorMap(A)
Base.transpose(A::DualTensorMap) = A.map

# properties
domain(A::DualTensorMap) = dual(codomain(A.map))
codomain(A::DualTensorMap) = dual(domain(A.map))

# comparison of DualTensorMap objects
==(A::DualTensorMap,B::DualTensorMap)=A.lmap==B.lmap

# multiplication
A_mul_B!{S,P}(Y::AbstractTensor{S,P},A::DualTensorMap{S,P},X::AbstractTensor{S,P})=At_mul_B!(Y,A.map,X)
At_mul_B!{S,P}(Y::AbstractTensor{S,P},A::DualTensorMap{S,P},X::AbstractTensor{S,P})=A_mul_B!(Y,A.map,X)


text:
# Script to run a test case for understanding how M scales with β #

# Add paths
loc_gm_package = "../src/"
if !(loc_gm_package in LOAD_PATH)
    push!(LOAD_PATH, "../src/")
end

# Required packages
using LinearAlgebra

using GraphicalModelLearning
using LightGraphs

# Function definitions
using Printf, SparseArrays
include("beta_scaling.jl")

using DelimitedFiles

# parameters
FLAG_create_struct_gm = true
FLAG_regular_random_gm = false
FLAG_lattice_gm = true
FLAG_weak_impurity = false
sampling_regime = M_regime()
learning_algo = NLP()

if FLAG_lattice_gm
    # text file to save adjacency matrix in
    file_adj_matrix_gm = "adj_matrix_ferro_lattice_gm_A_M.txt"
    # name of picture to save graphical model in
    file_plot_gm = "ferro_lattice_gm_A_M.eps"
    # File to save final results
    file_M_opt_gm = "M_opt_Ferro_Lattice_A_M_tol.txt"
elseif FLAG_regular_random_gm
    # text file to save adjacency matrix in
    file_adj_matrix_gm = "adj_matrix_spin_glass_gm.txt"
    # name of picture to save graphical model in
    file_plot_gm = "spin_glass_gm3.eps"
    # File to save final results
    file_M_opt_gm = "M_opt_M_Regime_Regular_Random_B.txt"
end

# Simple test
N = 16
d = 4
α = 0.4
#β_array = [0.8+0.2*i for i=1:11]
β_array = [1.1,1.5,1.9,2.3]

# Create and plot the initial graphical graphical model
β = copy(β_array[4])
adj_matrix, struct_adj_matrix = ferro_lattice(N,α,β,FLAG_weak_impurity)
m = n = Int(sqrt(N))
open(file_adj_matrix_gm, "w") do io
    writedlm(io, adj_matrix)
end;

## Start the complexity studies
τ = α/2
L = 45
M_factor = 0.05

M_opt = Array{Int64,1}(undef,length(β_array))

let
    M_guess = 9000
    for i = 1:length(β_array)
        # define β
        β = copy(β_array[i])
        println(β)
        println(M_guess)

        # Create the adjacency matrix
        adj_matrix, struct_adj_matrix = ferro_lattice(N,α,β,FLAG_weak_impurity)

        # Get the optimum number of samples
        M_opt[i] = get_M_opt_glauber_dynamics(adj_matrix, RPLE(), learning_algo, sampling_regime, τ, L, M_guess, M_factor)

        # Update the guess of M_opt
        @printf("beta=%f, M_opt=%d \n", β, M_opt[i])
        M_guess = Int(floor(0.8*copy(M_opt[i])))
    end
end

writedlm(file_M_opt_gm,M_opt)


text:
"""
	clip(cmd0::String="", arg1=nothing; kwargs...)

Reads (length,azimuth) pairs from file and plot a windclip diagram.

Full option list at [`psclip`]($(GMTdoc)psclip.html)

Parameters
----------

- **C** | **end_clip_path** :: [Type => Bool]

    Mark end of existing clip path. No input file is needed.
    ($(GMTdoc)psclip.html#c)
- $(GMT.opt_J)

- **A** | **inc** :: [Type => Str or []]

    By default, geographic line segments are connected as great circle arcs. To connect them as straight lines, use **A**
    ($(GMTdoc)psclip.html#a)
- $(GMT.opt_B)
- $(GMT.opt_Jz)
- **N** | **invert** :: [Type => Bool]

    Invert the sense of the test, i.e., clip regions where there is data coverage.
    ($(GMTdoc)psclip.html#n)
- $(GMT.opt_P)
- $(GMT.opt_R)
- **T** | **clip_limits** :: [Type => Bool]

    Rather than read any input files, simply turn on clipping for the current map region.
    ($(GMTdoc)psclip.html#t)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function clip(cmd0::String="", arg1=nothing; first=true, kwargs...)

	gmt_proggy = (IamModern[1]) ? "clip "  : "psclip "
	length(kwargs) == 0 && return monolitic(gmt_proggy, cmd0, arg1)

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode
	K, O = set_KO(first)		# Set the K O dance

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12c/12c")
	cmd, = parse_common_opts(d, cmd, [:UVXY :JZ :c :e :f :g :p :t :yx :params], first)
	cmd  = parse_these_opts(cmd, d, [[:A :straight_lines], [:C :end_clip_path], [:N :invert], [:T :clip_limits]])

	# If file name sent in, read it and compute a tight -R if this was not provided
	cmd, arg1, = read_data(d, cmd0, cmd, arg1, opt_R)

	return finish_PS_module(d, gmt_proggy * cmd, "", K, O, true, arg1)
end

# ---------------------------------------------------------------------------------------------------
clip!(cmd0::String="", arg1=nothing; kw...) = clip(cmd0, arg1; first=false, kw...)
clip(arg1, cmd0::String=""; kw...)  = clip(cmd0, arg1; first=true, kw...)
clip!(arg1, cmd0::String=""; kw...) = clip(cmd0, arg1; first=false, kw...)

psclip  = clip			# Alias
psclip! = clip!			# Alias


text:
#= ######################################################
includes
###################################################### =#

include("./models2D.jl")

#= ######################################################
inside functions
###################################################### =#

function get_epsilon_effective(m, mat, p_hole; epsilon_ink = 0.1)
    # calculate epsilon effective
    # hole is part with epsilon effective / encl is enclosure
    p_encl = collect(1:m.npar)
    filter!(x->x != p_hole, p_encl)
    temp_encl = 0
    temp_hole = 1000 # Kelvin
    temp = zeros(m.nelem,1)
    set_bc_part!(m, temp, p_hole, temp_hole)
    set_bc_part!(m, temp, p_encl, temp_encl)
    epsilon_hole = 1.0
    n_epsilon_encl = round(Integer,1.0/epsilon_ink)
    epsilon_encl = Matrix{Float64}(undef,n_epsilon_encl,2)
    epsilon_encl[:,1] = [epsilon_ink*i for i =1:n_epsilon_encl]
    epsilon = zeros(m.nelem,1)
    set_bc_part!(m, epsilon, p_hole, epsilon_hole)
    for i = 1:n_epsilon_encl
        set_bc_part!(m, epsilon, p_encl, epsilon_encl[i,1])
        # calculation of Qp
        Qp, Ga = tempsolver(m, mat, temp, epsilon)
        area_hole = get_area_of_part(m, p_hole)
        G_hole = sum(Ga[m.elem2par[p_hole].first:m.elem2par[p_hole].last])
        Qp_hole = sum(Qp[m.elem2par[p_hole].first:m.elem2par[p_hole].last])
        # calculation of epsilon effective
        epsilon_eff = 1.0 - (G_hole / (area_hole * sigma * (temp_hole^4 - temp_encl^4)))
        epsilon_encl[i,2] = epsilon_eff
    end
    return epsilon_encl
end


text:
export AbstractExporter,
    InMemoryExporter,
    ConsoleExporter,
    EXPORT_SUCCESS,
    EXPORT_FAILURE,
    export!,
    shut_down!,
    force_flush!

@enum ExportResult begin
    EXPORT_SUCCESS
    EXPORT_FAILURE
end

"""
An `AbstractExporter` is to export a collection of [`AbstractSpan`](@ref)s and
(or) [`Metric`](@ref)s. Each method should have the following interfaces
implemented:

  - [`force_flush!(::AbstractExporter)`](@ref)
  - [`shut_down!(::AbstractExporter)`](@ref)
  - [`export!(::AbstractExporter, collection)`](@ref)
"""
abstract type AbstractExporter end

force_flush!(::AbstractExporter) = true

shut_down!(e::AbstractExporter) = force_flush!(e)

#####

"""
    InMemoryExporter(;pool=[], is_shut_down=Ref(false))

Simply store all `export!`ed elements into the `pool`.
"""
Base.@kwdef struct InMemoryExporter <: AbstractExporter
    pool::Vector = []
    is_shut_down::Ref{Bool} = Ref(false)
end

Base.show(io::IO, e::InMemoryExporter) =
    print(io, "InMemoryExporter with $(length(e.pool)) elements")
Base.empty!(se::InMemoryExporter) = empty!(se.pool)
Base.length(se::InMemoryExporter) = length(se.pool)
Base.iterate(e::InMemoryExporter, args...) = iterate(e.pool, args...)

function force_flush!(se::InMemoryExporter)
    empty!(se)
    true
end

function shut_down!(se::InMemoryExporter)
    force_flush!(se)
    se.is_shut_down[] = true
    true
end

function export!(e::InMemoryExporter, x)
    if e.is_shut_down[]
        EXPORT_FAILURE
    else
        append!(e.pool, x)
        EXPORT_SUCCESS
    end
end

#####

"""
    ConsoleExporter(;io=stdout)

Print an [`AbstractSpan`](@ref) or [`Metric`](@ref) into `io`.
"""
Base.@kwdef struct ConsoleExporter <: AbstractExporter
    io::IO = stdout
end

function export!(ce::ConsoleExporter, batch)
    for x in batch
        pprint(ce.io, x)
        println(ce.io)
        print(ce.io, repeat("-", displaysize(ce.io)[2]))
        println(ce.io)
    end
    EXPORT_SUCCESS
end


text:
using Random

struct RandomWeightAdjustmentPolicy <: AbstractWeightAdjustmentPolicy
    weights::Vector{Float64}
    rng::MersenneTwister
end

RandomWeightAdjustmentPolicy(weights::AbstractArray{T}) where T <: Real = RandomWeightAdjustmentPolicy(weights, MersenneTwister())
Random.seed!(rwap::RandomWeightAdjustmentPolicy, seed::Integer) = Random.seed!(rwap.rng, seed)

possible_weights(rwap::RandomWeightAdjustmentPolicy) = rwap.weights

(rwap::RandomWeightAdjustmentPolicy)(;kwargs...) = rand(rwap.rng, rwap.weights)

function rwastar_search(search_prob::AbstractSearchProblem, weights::AbstractArray{T}, walltime_limit::Real, nodes_budget::Integer, rwa_seed::Integer) where T <: Real
    rwap = RandomWeightAdjustmentPolicy(weights)
    Random.seed!(rwap, rwa_seed)
    return awastar_search_with_dynamic_weights(search_prob, rwap, walltime_limit, nodes_budget)
end


text:
using Test
using ScoreDrivenERGM

# test that the SdErgm interface works for a concrete model<:SdErgm

# sample dgp

# correctly specified filter

# estimate


text:
using Pkg;
Pkg.add([
	"CSV",
	"DataFrames",
	"JSON"
])

max_stars_repo_name: keboola/docker-custom-julia

max_issues_repo_name: keboola/docker-custom-julia

max_forks_repo_name: keboola/docker-custom-julia


text:
abstract type AbstractRepIterator{T, ElemT} end

# Subtyping AbstractVector{ElemT} make Base think that RepIterator implements indexing e.g. for copy!
abstract type AbstractSingleRepIterator{T, ElemT, PT} <: AbstractRepIterator{T, ElemT} end

# SingleRepIterator
struct SingleRepIterator{T, ElemT, PT<:Rep{T}} <: AbstractSingleRepIterator{T, ElemT, PT}
    p::PT
    function SingleRepIterator{T, ElemT}(p::PT) where {T, ElemT, PT<:Rep{T}}
        new{T, ElemT, PT}(p)
    end
end
mapitem(it::SingleRepIterator, item) = item
iterator(T::Type, ElemT::Type, p::Rep) = SingleRepIterator{T, ElemT}(p)

# SingleMapRepIterator
struct SingleMapRepIterator{T, ElemT, PT<:Rep} <: AbstractSingleRepIterator{T, ElemT, PT}
    p::PT
    f::Function
    function SingleMapRepIterator{T, ElemT}(p::PT, f::Function) where {T, ElemT, PT<:Rep}
        new{T, ElemT, PT}(p, f)
    end
end
mapitem(it::SingleMapRepIterator, item) = it.f(1, item)
iterator(T::Type, ElemT::Type, f::Function, p::Rep) = SingleMapRepIterator{T, ElemT}(p, f)

# For a RepIterator with only one representation, we fallback to the indexing interface, the index being the iterator state
# A representation can overwrite this if it can do something more efficient or if it simply does not support indexing
function Base.eachindex(it::AbstractSingleRepIterator{<:Any, ElemT,
                                              RepT}) where {T, ElemT, RepT<:Rep{T}}
    return Indices{T, similar_type(ElemT, FullDim(RepT), T)}(it.p)
end
element_and_index(it, idx::Nothing) = nothing
function element_and_index(it::AbstractSingleRepIterator, idx::Index)
    return mapitem(it, get(it.p, idx)), idx
end
function Base.iterate(it::AbstractSingleRepIterator)
    idx = undouble_it(iterate(eachindex(it)))
    return element_and_index(it, idx)
end
function Base.iterate(it::AbstractSingleRepIterator, idx::Index)
    idx = undouble_it(iterate(eachindex(it), idx))::Union{Nothing, typeof(idx)}
    return element_and_index(it, idx)
end

abstract type AbstractMultiRepIterator{T, ElemT, PT} <: AbstractRepIterator{T, ElemT} end

# If there are multiple representations, we need to iterate.
# Builds a SingleRepIterator{ElemT} from p
function checknext(it::AbstractMultiRepIterator, i::Int, item_state)
    while i <= length(it.ps) && item_state === nothing
        i += 1
        if i <= length(it.ps)
            item_state = iterate(repit(it, i))
        end
    end
    if item_state === nothing
        return nothing
    else
        return item_state[1], (i, item_state[2])
    end
end
Base.iterate(it::AbstractMultiRepIterator) = checknext(it, 0, nothing)
function Base.iterate(it::AbstractMultiRepIterator, state)
    return checknext(it, state[1], iterate(repit(it, state[1]), state[2]))
end

# RepIterator
struct RepIterator{T, ElemT, PT<:Tuple{Vararg{Rep{T}}}} <: AbstractMultiRepIterator{T, ElemT, PT}
    ps::PT
    function RepIterator{T, ElemT}(ps::PT) where {T, ElemT, PT<:Tuple{Vararg{Rep{T}}}}
        new{T, ElemT, PT}(ps)
    end
end
repit(it::RepIterator{T, ElemT}, i::Int) where {T, ElemT} = SingleRepIterator{T, ElemT}(it.ps[i])

iterator(T::Type, ElemT::Type, p::Rep...) = RepIterator{T, ElemT}(p)

# MapRepIterator
struct MapRepIterator{T, ElemT, PT} <: AbstractMultiRepIterator{T, ElemT, PT}
    ps::PT
    f::Function
    function MapRepIterator{T, ElemT}(ps::PT, f::Function) where {T, ElemT, PT<:Tuple}
        new{T, ElemT, PT}(ps, f)
    end
end
# T is the coeftype, it can be different that the ones of is.ps[i]
repit(it::MapRepIterator{T, ElemT}, i::Int) where {T, ElemT} = SingleMapRepIterator{T, ElemT}(it.ps[i], (j, item) -> it.f(i, item)) # j should be 1

iterator(T::Type, ElemT::Type, f::Function, p::Rep...) = MapRepIterator{T, ElemT}(p, f)

function typed_map(f::Function, ::Type{T}, it::SingleRepIterator{Tin, ElemT}) where {T, Tin, ElemT}
    SingleMapRepIterator{T, similar_type(ElemT, T)}(it.p, f)
end
function typed_map(f::Function, ::Type{T}, it::RepIterator{Tin, ElemT}) where {T, Tin, ElemT}
    MapRepIterator{T, similar_type(ElemT, T)}(it.ps, f)
end

Base.eltype(it::AbstractRepIterator{T, ElemT}) where {T, ElemT} = ElemT
function change_coefficient_type(it::AbstractRepIterator, T::Type)
    typed_map((i, x) -> convert(similar_type(typeof(x), T), x), T, it)
end

# FIXME the variables need to be defined outside of the local scope of for
#       for Julia to know them inside the """ ... """ of the docstrings
singular = HorV = HorVRep = horvrep = singularlin = plural = plurallin = lenp = isnotemptyp = repexem = listexem = :_

for (isVrep, elt, loop_singular) in [(true, :AbstractVector, :point),
                                     (true, :Line, :line), (true, :Ray, :ray),
                                     (false, :HyperPlane, :hyperplane), (false, :HalfSpace, :halfspace)]
    global singular = loop_singular
    if isVrep
        vectortype_fun = :vvectortype
        global HorV = :V
        global HorVRep = :VRep
        global horvrep = :vrep
    else
        vectortype_fun = :hvectortype
        global HorV = :H
        global HorVRep = :HRep
        global horvrep = :hrep
    end
    singularstr = string(singular)
    elemtype = Symbol(singularstr * "type")
    donep = Symbol("done" * singularstr)
    startp = Symbol("start" * singularstr)
    nextp = Symbol("next" * singularstr)
    pluralstr = singularstr * "s"
    global plural = Symbol(pluralstr)
    global lenp = Symbol("n" * pluralstr)
    global isnotemptyp = Symbol("has" * pluralstr)
    mapit = Symbol("map" * pluralstr)
    inc = Symbol("incident" * pluralstr)
    incidx = Symbol("incident" * singularstr * "indices")

    @eval begin
        export $plural, $lenp, $isnotemptyp, $startp, $donep, $nextp, $elemtype
        export $inc, $incidx

        """
            $plural($horvrep::$HorVRep)

        Returns an iterator over the $plural of the $HorV-representation `$horvrep`.
        """
        function $plural end

        """
            incident$plural(p::Polyhedron, idx)

        Returns the list of $plural incident to idx for the polyhedron `p`.
        """
        $inc(p::Polyhedron{T}, idx) where {T} = get(p, IncidentElements{T, $elemtype(p)}(p, idx))

        """
            incident$(singular)indices(p::Polyhedron, idx)

        Returns the list of the indices of $plural incident to idx for the polyhedron `p`.
        """
        $incidx(p::Polyhedron{T}, idx) where {T} = get(p, IncidentIndices{T, $elemtype(p)}(p, idx))

        """
            $lenp($horvrep::$HorVRep)

        Returns the number of $plural of the $HorV-representation `$horvrep`.
        """
        $lenp($horvrep::$HorVRep{T}) where {T} = length(Indices{T, $elemtype($horvrep)}($horvrep))

        """
            $isnotemptyp($horvrep::$HorVRep)

        Returns whether the $HorV-representation `$horvrep` has any $singular.
        """
        $isnotemptyp($horvrep::$HorVRep{T}) where {T} = !isempty(Indices{T, $elemtype($horvrep)}($horvrep))

        $elemtype(p::Polyhedron) = $elemtype($horvrep(p))
        if $singularstr == "point"
            $elemtype(p::$HorVRep) = $vectortype_fun(typeof(p))
        else
            $elemtype(p::$HorVRep{T}) where {T} = $elt{T, $vectortype_fun(typeof(p))}
        end

        function $plural(p::$HorVRep{T}...) where {T}
            ElemT = promote_type($elemtype.(p)...)
            iterator(T, ElemT, p...)
        end

        function $mapit(f::Function, d::FullDim, ::Type{T}, p::$HorVRep...) where {T}
            ElemT = promote_type(similar_type.($elemtype.(p), Ref(d), T)...)
            iterator(T, ElemT, f, p...)
        end

        Base.length(it::AbstractSingleRepIterator{T, <:$elt}) where {T} = $lenp(it.p)
        Base.length(it::AbstractMultiRepIterator{T, <:$elt}) where {T} = sum($lenp, it.ps)
        Base.isempty(it::AbstractSingleRepIterator{T, <:$elt}) where {T} = !any($isnotemptyp(it.p))
        Base.isempty(it::AbstractMultiRepIterator{T, <:$elt}) where {T}  = !any($isnotemptyp.(it.ps))
    end
end

# Combines an element type with its linear version.
# e.g. combines rays with the lines by splitting lines in two points.
struct AllRepIterator{T, ElemT, LinElemT, LRT<:Union{AbstractRepIterator{T, LinElemT}}, RT<:Union{AbstractRepIterator{T, ElemT}}}
    itlin::LRT
    it::RT
end

Base.eltype(it::AllRepIterator{T, ElemT}) where {T, ElemT} = ElemT
Base.length(it::AllRepIterator) = 2length(it.itlin) + length(it.it)
Base.isempty(it::AllRepIterator) = isempty(it.itlin) && isempty(it.it)

function checknext(it::AllRepIterator, i, item_state)
    while i <= 3 && item_state === nothing
        i += 1
        if i <= 2
            @assert i >= 1
            item_state = iterate(it.itlin)
        elseif i == 3
            item_state = iterate(it.it)
        end
    end
    if item_state === nothing
        return nothing
    else
        if i <= 2
            @assert i >= 1
            item = splitlin(item_state[1], i)
        else
            @assert i == 3
            item = item_state[1]
        end
        return item, (i, item_state[2])
    end
end

splitlin(h::HyperPlane, i) = (i == 1 ? HalfSpace(h.a, h.β) : HalfSpace(-h.a, -h.β))
#splitlin(s::SymPoint, i) = (i == 1 ? coord(s) : -coord(s))
splitlin(l::Line, i) = (i == 1 ? Ray(coord(l)) : Ray(-coord(l)))

Base.iterate(it::AllRepIterator) = checknext(it, 0, nothing)
function Base.iterate(it::AllRepIterator, istate)
    i, state = istate
    if i <= 2
        @assert i >= 1
        item_state = iterate(it.itlin, state)
    else
        @assert i == 3
        item_state = iterate(it.it, state)
    end
    return checknext(it, i, item_state)
end

for (isVrep, loop_singularlin,
     loop_singular, loop_repexem,
     loop_listexem) in [#(true, :sympoint, :point, "convexhull(SymPoint([1, 0]), [0, 1])", "[1, 0], [-1, 0], [0, 1]"),
                   (true, :line, :ray, "Line([1, 0]) + Ray([0, 1])",
                    "Ray([1, 0]), Ray([-1, 0]), Ray([0, 1])"),
                   (false, :hyperplane, :halfspace,
                    "HyperPlane([1, 0], 1) ∩ HalfSpace([0, 1], 1)",
                    "HalfSpace([1, 0]), HalfSpace([-1, 0]), HalfSpace([0, 1])")]
    global singularlin = loop_singularlin
    global singular = loop_singular
    global repexem = loop_repexem
    global listexem = loop_listexem
    if isVrep
        global HorV = :V
        global HorVRep = :VRep
        global horvrep = :vrep
    else
        global HorV = :H
        global HorVRep = :HRep
        global horvrep = :hrep
    end
    pluralstrlin = string(singularlin) * "s"
    global plurallin = Symbol(pluralstrlin)
    lenplin = Symbol("n" * pluralstrlin)
    isnotemptyplin = Symbol("has" * pluralstrlin)
    pluralstr = string(singular) * "s"
    global plural = Symbol(pluralstr)
    global lenp = Symbol("n" * pluralstr)
    global isnotemptyp = Symbol("has" * pluralstr)
    allpluralstr = "all" * pluralstr
    allplural = Symbol(allpluralstr)
    alllenp = Symbol("n" * allpluralstr)
    allisnotemptyp = Symbol("has" * allpluralstr)
    @eval begin
        export $allplural, $alllenp, $allisnotemptyp

        """
            all$plural($horvrep::$HorVRep)

        Returns an iterator over the $plural and $plurallin in the $HorV-representation `$horvrep` splitting $plurallin in two $plural.

        ### Examples

        ```julia
        $horvrep = $repexem
        collect(all$plural($horvrep)) # Returns [$listexem]
        ```
        """
        $allplural(p::$HorVRep...) = AllRepIterator($plurallin(p...), $plural(p...))

        """
            nall$plural($horvrep::$HorVRep)

        Returns the number of $plural plus twice the number of $plurallin in the $HorV-representation `$horvrep`, i.e. `length(all$plural($horvrep))`
        """
        $alllenp($horvrep::$HorVRep) = 2 * $lenplin($horvrep) + $lenp($horvrep)

        """
            hasall$plural($horvrep::$HorVRep)

        Returns whether the $HorV-representation `$horvrep` contains any $singular or $singularlin.
        """
        $allisnotemptyp($horvrep::$HorVRep) = $isnotemptyplin($horvrep) || $isnotemptyp($horvrep)
    end
end

const ElemIt{ElemT} = Union{AllRepIterator{<:Any, ElemT}, AbstractRepIterator{<:Any, ElemT}, AbstractVector{ElemT}}
const HyperPlaneIt{T} = ElemIt{<:HyperPlane{T}}
const HalfSpaceIt{T} = ElemIt{<:HalfSpace{T}}
const HIt{T} = Union{HyperPlaneIt{T}, HalfSpaceIt{T}}

const PointIt{T} = ElemIt{<:AbstractVector{T}}
const PIt{T} = PointIt{T}
const LineIt{T} = ElemIt{<:Line{T}}
const RayIt{T} = ElemIt{<:Ray{T}}
const RIt{T} = Union{LineIt{T}, RayIt{T}}
const VIt{T} = Union{PIt{T}, RIt{T}}

const It{T} = Union{HIt{T}, VIt{T}}

function fillvits(d::FullDim, points::ElemIt{AT},
                  lines::ElemIt{Line{T, AT}}=Line{T, AT}[],
                  rays::ElemIt{Ray{T, AT}}=Ray{T, AT}[]) where {T, AT<:AbstractVector{T}}
    if isempty(points) && !(isempty(lines) && isempty(rays))
        vconsistencyerror()
    end
    return points, lines, rays
end
function fillvits(d::FullDim, lines::ElemIt{Line{T, AT}},
                  rays::ElemIt{Ray{T, AT}}=Ray{T, AT}[]) where {T, AT}
    d = FullDim_rec(lines, rays)
    N = fulldim(d)
    if isempty(lines) && isempty(rays)
        points = AT[]
    else
        points = [origin(AT, N)]
    end
    return points, lines, rays
end

FullDim_hreps(p...) = FullDim(p[1]), hreps(p...)...
FullDim_vreps(p...) = FullDim(p[1]), vreps(p...)...

hreps(p::HRep{T}...) where {T} = hyperplanes(p...), halfspaces(p...)
hreps(p::HAffineSpace{T}...) where {T} = tuple(hyperplanes(p...))

hmap(f, d::FullDim, ::Type{T}, p::HRep...) where T = maphyperplanes(f, d, T, p...), maphalfspaces(f, d, T, p...)
hmap(f, d::FullDim, ::Type{T}, p::HAffineSpace...) where T = tuple(maphyperplanes(f, d, T, p...))

hconvert(RepT::Type{<:HRep{T}}, p::HRep{T}) where {T} = constructpolyhedron(RepT, FullDim(p), (p,), hreps(p)...)
hconvert(RepT::Type{<:HRep{T}}, p::HRep)    where {T} = constructpolyhedron(RepT, FullDim(p), (p,), change_coefficient_type.(hreps(p), T)...)

vreps(p...) = preps(p...)..., rreps(p...)...
preps(p::VRep...) = tuple(points(p...))
preps(p::VCone...) = tuple()
rreps(p::VRep...) = lines(p...), rays(p...)
rreps(p::VLinearSpace...) = tuple(lines(p...))
rreps(p::VPolytope...) = tuple()

vmap(f, d::FullDim, ::Type{T}, p::VRep...) where T = pmap(f, d, T, p...)..., rmap(f, d, T, p...)...
pmap(f, d::FullDim, ::Type{T}, p::VRep...) where T = tuple(mappoints(f, d, T, p...))
pmap(f, d::FullDim, ::Type, p::VCone...) = tuple()
rmap(f, d::FullDim, ::Type{T}, p::VRep...) where T = maplines(f, d, T, p...), maprays(f, d, T, p...)
rmap(f, d::FullDim, ::Type{T}, p::VLinearSpace...) where T = tuple(maplines(f, d, T, p...))
rmap(f, d::FullDim, ::Type, p::VPolytope...) = tuple()

vconvert(RepT::Type{<:VRep{T}}, p::VRep{T}) where {T} = constructpolyhedron(RepT, FullDim(p), (p,), vreps(p)...)
function vconvert(RepT::Type{<:VRep{T}}, p::VRep)    where {T}
    constructpolyhedron(RepT, FullDim(p), (p,), change_coefficient_type.(vreps(p), T)...)
end


text:
module MethodAnalysis

using AbstractTrees

using Base: Callable, IdSet
using Core: MethodInstance, CodeInfo, SimpleVector, MethodTable
using Base.Meta: isexpr

export visit, call_type, methodinstance, methodinstances, worlds  # findcallers is exported from its own file
export visit_backedges, all_backedges, with_all_backedges, terminal_backedges, direct_backedges

include("visit.jl")
include("backedges.jl")

if !hasmethod(==, Tuple{Core.PhiNode,Core.PhiNode})
    Base.:(==)(stmt1::Core.PhiNode, stmt2::Core.PhiNode) = stmt1.edges == stmt2.edges && stmt1.values == stmt2.values
end

"""
    call_type(tt)

Split a signature type like `Tuple{typeof(f),ArgTypes...}` back out to `(f, Tuple{ArgTypes...})`
"""
function call_type(tt)
    ft = tt.parameters[1]
    argt = Tuple{tt.parameters[2:end]...}
    name = Symbol(String(ft.name.name)[2:end])  # strip off leading '#'
    return (getfield(ft.name.module, name), argt)
end

"""
    minmaxs = worlds(mi::MethodInstance)

Collect the (min,max) world-age pairs for all CodeInstances associated with `mi`.
"""
function worlds(mi::Core.MethodInstance)
    w = Tuple{UInt,UInt}[]
    if isdefined(mi, :cache)
        ci = mi.cache
        push!(w, (ci.min_world % UInt, ci.max_world % UInt))
        while isdefined(ci, :next)
            ci = ci.next
            push!(w, (ci.min_world % UInt, ci.max_world % UInt))
        end
    end
    return w
end

# Not sure we want to change the meaning of == here, so let's define our own name
# A few fields are deliberately unchecked
function equal(ci1::Core.CodeInfo, ci2::Core.CodeInfo)
    ret = ci1.code == ci2.code &&
          ci1.codelocs == ci2.codelocs &&
          ci1.ssavaluetypes == ci2.ssavaluetypes &&
          ci1.ssaflags == ci2.ssaflags &&
          ci1.method_for_inference_limit_heuristics == ci2.method_for_inference_limit_heuristics &&
          ci1.linetable == ci2.linetable &&
          ci1.slotnames == ci2.slotnames &&
          ci1.slotflags == ci2.slotflags
    if VERSION >= v"1.2"
        ret &= ci1.slottypes == ci2.slottypes &&
               ci1.rettype == ci2.rettype
    end
    return ret
end
equal(p1::Pair, p2::Pair) = p1.second == p2.second && equal(p1.first, p2.first)

"""
    mi = methodinstance(f, types)
    mi = methodinstance(tt::Type{<:Tuple})

Return the `MethodInstance` `mi` for function `f` and the given `types`,
or for the complete signature `tt`.
If no version compiled for these types exists, returns `nothing`.

# Examples

```jldoctest; setup=:(using MethodAnalysis)
julia> f(x, y::String) = 2x; f(x, y::Number) = x + y;

julia> f(1, "hi"); f(1, 1.0);

julia> methodinstance(f, (Int, String))
MethodInstance for f(::Int64, ::String)

julia> methodinstance(Tuple{typeof(f), Int, String})
MethodInstance for f(::Int64, ::String)
```
"""
function methodinstance(@nospecialize(f), @nospecialize(types))
    if types isa Tuple
        tt = Tuple{typeof(f), types...}
        return methodinstance(f, tt)
    end
    inst = nothing
    visit(f) do mi
        if isa(mi, MethodInstance)
            if mi.specTypes === types
                inst = mi
            end
            return false
        end
        true
    end
    return inst
end
function methodinstance(@nospecialize(types))
    f, argt = call_type(types)
    return methodinstance(f, types)
end

"""
    methodinstances()
    methodinstances(mod::Module)
    methodinstances(f)

Collect all `MethodInstance`s, optionally restricting them to a particular module, function, method, or methodlist.
"""
function methodinstances(top=())
    if isa(top, Module) || isa(top, Function) || isa(top, Type) || isa(top, Method) || isa(top, Base.MethodList)
        top = (top,)
    end
    mis = Core.MethodInstance[]
    visit(top...) do item
        isa(item, Core.MethodInstance) || return true
        push!(mis, item)
        false
    end
    return mis
end

if isdefined(Core, :MethodMatch)
    include("findcallers.jl")
end

# AbstractTrees interface
AbstractTrees.children(mi::MethodInstance) = isdefined(mi, :backedges) ? mi.backedges : []

# deprecations
@deprecate instance methodinstance

end # module


text:
module AbstractMover

export AbstractMoverDev, numaxes
export move, moveX, moveY, moveZ
export rmove, rmoveX, rmoveY, rmoveZ

export absposition, abspositionX, abspositionY, abspositionZ
export position, positionX, positionY, positionZ

export setreference, setreferenceX, setreferenceY, setreferenceZ
export setabsreference, setabsreferenceX, setabsreferenceY, setabsreferenceZ

export home, homeX, homeY, homeZ
export stopmotion, waituntildone

abstract type AbstractMoverDev end

numaxes(dev::AbstractMoverDev) = error("Not implemented!")
numaxes(::Type{<:AbstractMoverDev}) = error("Not implemented!")

move(dev::AbstractMoverDev) = error("Not implemented!")

moveX(dev::AbstractMoverDev) = error("Not implemented!")
moveY(dev::AbstractMoverDev) = error("Not implemented!")
moveZ(dev::AbstractMoverDev) = error("Not implemented!")

rmove(dev::AbstractMoverDev) = error("Not implemented!")

rmoveX(dev::AbstractMoverDev) = error("Not implemented!")
rmoveY(dev::AbstractMoverDev) = error("Not implemented!")
rmoveZ(dev::AbstractMoverDev) = error("Not implemented!")

absposition(dev::AbstractMoverDev) = error("Not implemented")
abspositionX(dev::AbstractMoverDev) = error("Not implemented")
abspositionY(dev::AbstractMoverDev) = error("Not implemented")
abspositionZ(dev::AbstractMoverDev) = error("Not implemented")

import Base
Base.position(dev::AbstractMoverDev) = error("Not implemented")

positionX(dev::AbstractMoverDev) = error("Not implemented")
positionY(dev::AbstractMoverDev) = error("Not implemented")
positionZ(dev::AbstractMoverDev) = error("Not implemented")

setreference(dev::AbstractMoverDev) = error("Not implemented")

setreferenceX(dev::AbstractMoverDev) = error("Not implemented")
setreferenceY(dev::AbstractMoverDev) = error("Not implemented")
setreferenceZ(dev::AbstractMoverDev) = error("Not implemented")

setabsreference(dev::AbstractMoverDev) = error("Not implemented")

setabsreferenceX(dev::AbstractMoverDev) = error("Not implemented")
setabsreferenceY(dev::AbstractMoverDev) = error("Not implemented")
setabsreferenceZ(dev::AbstractMoverDev) = error("Not implemented")

home(dev::AbstractMoverDev) = error("Not implemented!")

homeX(dev::AbstractMoverDev) = error("Not implemented!")
homeY(dev::AbstractMoverDev) = error("Not implemented!")
homeZ(dev::AbstractMoverDev) = error("Not implemented!")

stopmotion(dev::AbstractMoverDev) = error("Not implemented!")

waituntildone(dev::AbstractMoverDev) = error("Not implemented!")

end


text:
using Oceananigans.Fields: AbstractField

import Oceananigans.AbstractOperations: ConditionalOperation, get_condition, truefunc
import Oceananigans.Fields: condition_operand, conditional_length

# ###
# ###  reduction operations involving immersed boundary grids exclude the immersed region
# ### (we exclude also the values on the faces of the immersed boundary with `solid_interface`)
# ###

const ImmersedField = AbstractField{<:Any, <:Any, <:Any, <:ImmersedBoundaryGrid}

struct NotImmersed{F} <: Function
    func :: F
end

@inline condition_operand(func::Function, operand::ImmersedField, condition, mask) = ConditionalOperation(operand; func, condition = NotImmersed(condition), mask)
@inline condition_operand(func::Function, operand::ImmersedField, ::Nothing, mask) = ConditionalOperation(operand; func, condition = NotImmersed(truefunc), mask)

@inline condition_operand(func::typeof(identity), operand::ImmersedField, condition, mask) = ConditionalOperation(operand; func, condition = NotImmersed(condition), mask)
@inline condition_operand(func::typeof(identity), operand::ImmersedField, ::Nothing, mask) = ConditionalOperation(operand; func, condition = NotImmersed(truefunc), mask)

@inline conditional_length(c::ImmersedField)       = conditional_length(condition_operand(identity, c, nothing, 0))
@inline conditional_length(c::ImmersedField, dims) = conditional_length(condition_operand(identity, c, nothing, 0), dims)

@inline function get_condition(condition::NotImmersed, i, j, k, ibg, op::ConditionalOperation, args...)
    LX, LY, LZ = location(op)
    return get_condition(condition.func, i, j, k, ibg, args...) & !(solid_interface(LX(), LY(), LZ(), i, j, k, ibg))
end


text:
const Kron{T,N,I<:Tuple} = ApplyArray{T,N,typeof(kron),I}

function instantiate(K::Applied{<:Any,typeof(kron)})
    isempty(K.args) && throw(ArgumentError("Cannot take kronecker product of empty vectors"))
    applied(kron, map(instantiate,K.args)...)
end

Kron(A...) = ApplyArray(kron, A...)
Kron{T}(A...) where T = ApplyArray{T}(kron, A...)

_kron_dims() = 0
_kron_dims(A, B...) = max(ndims(A), _kron_dims(B...))

eltype(A::Applied{<:Any,typeof(kron)}) = promote_type(map(eltype,A.args)...)
ndims(A::Applied{<:Any,typeof(kron)}) = _kron_dims(A.args...)

size(K::Kron, j::Int) = prod(size.(K.args, j))
size(a::Kron{<:Any,1}) = (size(a,1),)
size(a::Kron{<:Any,2}) = (size(a,1), size(a,2))
size(a::Kron{<:Any,N}) where {N} = (@_inline_meta; ntuple(M -> size(a, M), Val(N)))
axes(a::Kron{<:Any,1}) = (OneTo(size(a,1)),)
axes(a::Kron{<:Any,2}) = (OneTo(size(a,1)), OneTo(size(a,2)))
axes(a::Kron{<:Any,N}) where {N} = (@_inline_meta; ntuple(M -> OneTo(size(a, M)), Val(N)))

function det(K::Kron{<:Any, 2})
    (size(K, 1) == size(K, 2)) || throw(DimensionMismatch("matrix is not square: dimensions are $(size(K))"))

    d = 1.
    s = size(K, 1)

    for A in K.args
        if size(A, 1) == size(A, 2)
            dA = det(A)
            if iszero(dA)
                return dA
            end
            d *= dA^(s ÷ size(A, 1))
        else
            # The Kronecker Product of rectangular matrices, if it is square, will
            # have determinant zero. This can be shown by using the fact that
            # rank(A ⊗ B) = rank(A)rank(B) and showing that this is strictly less
            # than the number of rows in the resulting Kronecker matrix. Hence,
            # since A ⊗ B does not have full rank, its determinant must be zero.
            return zero(d)
        end
    end
    return d
end

function tr(K::Kron{<:Any, 2})
    (size(K, 1) == size(K, 2)) || throw(DimensionMismatch("matrix is not square: dimensions are $(size(K))"))
    if all(A -> (size(A, 1) == size(A, 2)), K.args)  # check if all component matrices are square
        return prod(tr.(K.args))
    else
        return sum(diag(K))
    end
end

kron_getindex((A,)::Tuple{AbstractVector}, k::Integer) = A[k]
function kron_getindex((A,B)::NTuple{2,AbstractVector}, k::Integer)
    K,κ = divrem(k-1, length(B))
    A[K+1]*B[κ+1]
end
kron_getindex((A,)::Tuple{AbstractMatrix}, k::Integer, j::Integer) = A[k,j]
function kron_getindex((A,B)::NTuple{2,AbstractArray}, k::Integer, j::Integer)
    K,κ = divrem(k-1, size(B,1))
    J,ξ = divrem(j-1, size(B,2))
    A[K+1,J+1]*B[κ+1,ξ+1]
end

kron_getindex(args::Tuple, k::Integer, j::Integer) = kron_getindex(tuple(Kron(args[1:2]...), args[3:end]...), k, j)
kron_getindex(args::Tuple, k::Integer) = kron_getindex(tuple(Kron(args[1:2]...), args[3:end]...), k)

getindex(K::Kron{<:Any,1}, k::Integer) = kron_getindex(K.args, k)
getindex(K::Kron{<:Any,2}, k::Integer, j::Integer) = kron_getindex(K.args, k, j)
getindex(K::Applied{DefaultArrayApplyStyle,typeof(kron)}, k::Integer) = kron_getindex(K.args, k)
getindex(K::Applied{DefaultArrayApplyStyle,typeof(kron)}, k::Integer, j::Integer) = kron_getindex(K.args, k, j)

## Adapted from julia/stdlib/LinearAlgebra/src/dense.jl kron definition
function _kron2!(R, K)
    size(R) == size(K) || throw(DimensionMismatch("Matrices have wrong dimensions"))
    a,b = K.args
    require_one_based_indexing(a, b)
    m = 1
    @inbounds for j = 1:size(a,2), l = 1:size(b,2), i = 1:size(a,1)
        aij = a[i,j]
        for k = 1:size(b,1)
            R[m] = aij*b[k,l]
            m += 1
        end
    end
    R
end

copyto!(R::AbstractMatrix, K::Kron{<:Any,2,<:Tuple{<:AbstractMatrix,<:AbstractMatrix}}) =
    _kron2!(R, K)
copyto!(R::AbstractVector, K::Kron{<:Any,1,<:Tuple{<:AbstractVector,<:AbstractVector}}) =
    _kron2!(R, K)
copyto!(R::AbstractMatrix{T}, K::Kron{T,2,<:Tuple{<:AbstractMatrix,<:AbstractMatrix}}) where T =
    _kron2!(R, K)
copyto!(R::AbstractVector{T}, K::Kron{T,1,<:Tuple{<:AbstractVector,<:AbstractVector}}) where T =
    _kron2!(R, K)

struct Diff{T, N, Arr} <: LazyArray{T, N}
    v::Arr
    dims::Int
end

Diff(v::AbstractVector{T}) where T = Diff{T,1,typeof(v)}(v, 1)

function Diff(A::AbstractMatrix{T}; dims::Integer) where T
    dims == 1 || dims == 2 || throw(ArgumentError("dimension must be 1 or 2, got $dims"))
    Diff{T,2,typeof(A)}(A, dims)
end

IndexStyle(::Type{<:Diff{<:Any,1}}) = IndexLinear()
IndexStyle(::Type{<:Diff{<:Any,2}}) = IndexCartesian()

size(D::Diff{<:Any,1}) = (length(D.v)-1,)
function size(D::Diff{<:Any,2})
    if D.dims == 1
        (size(D.v,1)-1,size(D.v,2))
    else #dims == 2
        (size(D.v,1),size(D.v,2)-1)
    end
end

getindex(D::Diff{<:Any, 1}, k::Integer) = D.v[k+1] - D.v[k]
function getindex(D::Diff, k::Integer, j::Integer)
    if D.dims == 1
        D.v[k+1,j] - D.v[k,j]
    else # dims == 2
        D.v[k,j+1] - D.v[k,j]
    end
end

struct Cumsum{T, N, Arr} <: LazyArray{T, N}
    v::Arr
    dims::Int
end

Cumsum(v::AbstractVector{T}) where T = Cumsum{T,1,typeof(v)}(v, 1)

function Cumsum(A::AbstractMatrix{T}; dims::Integer) where T
    dims == 1 || dims == 2 || throw(ArgumentError("dimension must be 1 or 2, got $dims"))
    Cumsum{T,2,typeof(A)}(A, dims)
end

IndexStyle(::Type{<:Cumsum{<:Any,1}}) = IndexLinear()
IndexStyle(::Type{<:Cumsum{<:Any,2}}) = IndexCartesian()

size(Q::Cumsum) = size(Q.v)

getindex(Q::Cumsum{<:Any, 1}, k::Integer) = k == 1 ? Q.v[1] : Q.v[k] + Q[k-1]
function getindex(Q::Cumsum, k::Integer, j::Integer)
    if Q.dims == 1
        k == 1 ? Q.v[1,j] : Q.v[k,j] + Q[k-1,j]
    else # dims == 2
        j == 1 ? Q.v[k,1] : Q.v[k,j] + Q[k,j-1]
    end
end

copyto!(x::AbstractArray{<:Any,N}, C::Cumsum{<:Any,N}) where N = cumsum!(x, C.v)

# keep lazy
cumsum(a::LazyArray; kwds...) = Cumsum(a; kwds...)


text:
"""
`Transition{T<:AbstractFloat}`

The transition equation of the state-space model takes the form

    `s_t = TTT*s_{t-1} + RRR*ϵ_t + CCC`

The `Transition` type stores the coefficient `Matrix{T}`s (`TTT`, `RRR`) and constant `Vector{T} CCC`.
"""
mutable struct Transition{T<:AbstractFloat}
    TTT::Matrix{T}
    RRR::Matrix{T}
    CCC::Vector{T}
end
function Transition(TTT::Matrix{T}, RRR::Matrix{T}) where T<:AbstractFloat
    CCC = zeros(eltype(TTT), size(TTT, 1))
    Transition{T}(TTT, RRR, CCC)
end
function Transition(TTT::Matrix{T}, RRR::Matrix{T}, CCC::Matrix{T}) where T<:AbstractFloat
    Transition{T}(TTT, RRR, collect(CCC))
end
function Base.getindex(eq::Transition, d::Symbol)
    if d in (:TTT, :RRR, :CCC)
        return getfield(eq, d)
    else
        throw(KeyError(d))
    end
end

"""
`Measurement{T<:AbstractFloat}`

The measurement equation of the state-space model takes the form

   `y_t = ZZ*s_t + DD + u_t`

where the error `u_t` is the measurement error, which is uncorrelated with the
shocks in the transition equation `ϵ_t`.

### Fields

If `Ns` is the number of states `s_t`, `Ny` is the number of
observables `y_t`, and `Ne` is the number of shocks `ϵ_t`:

- `ZZ`: the `Ny` x `Ns`  measurement matrix
- `DD`: the `Ny` x 1 constant vector
- `QQ`: the `Ne` x `Ne` covariance matrix for the shocks `ϵ_t`
- `EE`: the `Ny` x `Ny` covariance matrix for the measurement error `η_t`
"""
mutable struct Measurement{T<:AbstractFloat}
    ZZ::Matrix{T}
    DD::Vector{T}
    QQ::Matrix{T}
    EE::Matrix{T}
end

function Base.getindex(M::Measurement, d::Symbol)
    if d in (:ZZ, :DD, :QQ, :EE)
        return getfield(M, d)
    else
        throw(KeyError(d))
    end
end

function measurement(m::AbstractDSGEModel, trans::Transition; shocks::Bool=true)
    TTT = trans[:TTT]
    RRR = trans[:RRR]
    CCC = trans[:CCC]
    measurement(m, TTT, RRR, CCC; shocks=shocks)
end

"""
```
PseudoMeasurement{T<:AbstractFloat}
```

The pseudo-measurement equation of the state-space model takes the form

   `x_t = ZZ_pseudo*s_t + DD_pseudo`

### Fields

Let `Ns` be the number of states `s_t` and `Nx` be the number of
pseudo-observables `x_t`:

- `ZZ_pseudo`: the `Nx` x `Ns` pseudo-measurement matrix
- `DD_pseudo`: the `Nx` x 1 constant vector
"""
mutable struct PseudoMeasurement{T<:AbstractFloat}
    ZZ_pseudo::Matrix{T}
    DD_pseudo::Vector{T}
end

function Base.getindex(M::PseudoMeasurement, d::Symbol)
    if d in (:ZZ_pseudo, :DD_pseudo)
        return getfield(M, d)
    else
        throw(KeyError(d))
    end
end

"""
`System{T<:AbstractFloat}`

A mutable struct containing the transition and measurement equations for a
state-space model. The matrices may be directly indexed: `sys[:TTT]`
returns `sys.transition.TTT`, etc.
"""
mutable struct System{T<:AbstractFloat}
    transition::Transition{T}
    measurement::Measurement{T}
    pseudo_measurement::PseudoMeasurement{T}
end

function System(transition::Transition{T}, measurement::Measurement{T}) where T<:AbstractFloat
    # Initialize empty pseudo-measurement equation
    _n_states = size(transition.TTT, 1)
    _n_pseudo = 0
    ZZ_pseudo = zeros(_n_pseudo, _n_states)
    DD_pseudo = zeros(_n_pseudo)
    pseudo_measurement = PseudoMeasurement(ZZ_pseudo, DD_pseudo)

    return System(transition, measurement, pseudo_measurement)
end

function Base.getindex(system::System, d::Symbol)
    if d in (:transition, :measurement, :pseudo_measurement)
        return getfield(system, d)
    elseif d in fieldnames(typeof(system.transition))
        return getfield(system.transition, d)
    elseif d in fieldnames(typeof(system.measurement))
        return getfield(system.measurement, d)
    elseif d in fieldnames(typeof(system.pseudo_measurement))
        return getfield(system.pseudo_measurement, d)
    else
        throw(KeyError(d))
    end
end

function Base.copy(system::System)
    trans = Transition(system[:TTT], system[:RRR], system[:CCC])
    meas  = Measurement(system[:ZZ], system[:DD], system[:QQ], system[:EE])
    pseudo_meas = PseudoMeasurement(system[:ZZ_pseudo], system[:DD_pseudo])
    return System(trans, meas, pseudo_meas)
end

"""
```
compute_system(m; apply_altpolicy = false)
```

Given the current model parameters, compute the state-space system
corresponding to model `m`. Returns a `System` object.
"""
function compute_system(m::AbstractDSGEModel{T}; apply_altpolicy = false,
                        verbose::Symbol = :high) where T<:AbstractFloat
    # Solve model
    TTT, RRR, CCC = solve(m; apply_altpolicy = apply_altpolicy, verbose = verbose)
    transition_equation = Transition(TTT, RRR, CCC)

    # Solve measurement equation
    measurement_equation = measurement(m, TTT, RRR, CCC)

    type_tuple = (typeof(m), Matrix{T}, Matrix{T}, Vector{T})
    if hasmethod(pseudo_measurement, type_tuple)
        # Solve pseudo-measurement equation
        pseudo_measurement_equation = pseudo_measurement(m, TTT, RRR, CCC)
        return System(transition_equation, measurement_equation, pseudo_measurement_equation)
    else
        return System(transition_equation, measurement_equation)
    end
end

"""
```
compute_system(m::PoolModel{T})
```

Given the current model parameters, compute the state-space system
corresponding to the PoolModel model `m`.

Outputs

```
Φ: state transition function
Ψ: likelihood function, given weights on underlying models (the states) and predictive densities
F_ϵ: structural shock distribution
F_u: likelihood function measurement error distribution
F_λ: initial distribution of λ for state transition function
"""
function compute_system(m::PoolModel{T};
                        verbose::Symbol = :high) where T<:AbstractFloat
    Φ, F_ϵ, F_λ = transition(m)
    Ψ, F_u = measurement(m)
    return Φ, Ψ, F_ϵ, F_u, F_λ
end

"""
```
compute_system_function(system::System{S}) where S<:AbstractFloat
```

### Inputs

- `system::System`: The output of compute_system(m), i.e. the matrix outputs from solving a given model, m.

### Outputs

- `Φ::Function`: transition equation
- `Ψ::Function`: measurement equation
- `F_ϵ::Distributions.MvNormal`: shock distribution
- `F_u::Distributions.MvNormal`: measurement error distribution
"""
function compute_system_function(system::System{S}) where S<:AbstractFloat
    # Unpack system
    TTT    = system[:TTT]
    RRR    = system[:RRR]
    CCC    = system[:CCC]
    QQ     = system[:QQ]
    ZZ     = system[:ZZ]
    DD     = system[:DD]
    EE     = system[:EE]

    # Define transition and measurement functions
    @inline Φ(s_t1::Vector{S}, ϵ_t::Vector{S}) = TTT*s_t1 + RRR*ϵ_t + CCC
    @inline Ψ(s_t::Vector{S}) = ZZ*s_t + DD

    # Define shock and measurement error distributions
    nshocks = size(QQ, 1)
    nobs    = size(EE, 1)
    F_ϵ = Distributions.MvNormal(zeros(nshocks), QQ)
    F_u = Distributions.MvNormal(zeros(nobs),    EE)

    return Φ, Ψ, F_ϵ, F_u
end

function zero_system_constants(system::System{S}) where S<:AbstractFloat
    system = copy(system)

    system.transition.CCC = zeros(size(system[:CCC]))
    system.measurement.DD = zeros(size(system[:DD]))
    system.pseudo_measurement.DD_pseudo = zeros(size(system[:DD_pseudo]))

    return system
end


text:
# #--------------------------------------------------------------------------------------------------------------------------
# #--------------------------------------------------------------------------------------------------------------------------
# # This file contains functions used for adding revenue recycling to the Nested Inequalities Climate Economy (NICE) model.
# #--------------------------------------------------------------------------------------------------------------------------
# #--------------------------------------------------------------------------------------------------------------------------

#####################################################################################################################
# CLEAN UP AND ORGANIZE INCOME SHARE DATA INTO NICE FORMAT.
#####################################################################################################################
# Description: This function takes a DataFrame of time varying quintile income distributions (for each NICE region
#              from 2015-2105) and converts it into an array with the dimension format required by NICE (time, region,
#              quintile). It sets the 2005 quintile share to the calculated 2015 value and assumes quintile shares after
#              2105 remain fixed.
#
# Function Arguments:
#
#       income_data   = A DataFrame read in from one of the income share files in the "data" folder of NICE.
#--------------------------------------------------------------------------------------------------------------------

function get_quintile_income_shares(income_data::DataFrame)

    # Allocate an array for cleaned up quintile share data in format required for NICE (time × regions × quintile).
    nice_years = collect(2005:10:2595)
    quintile_shares = zeros(length(nice_years), 12, 5);

    # Create vector of years for SSP quintile data (spans 2015-2105 in ten-year increments).
    ssp_years = unique(income_data[:,:decade])

    # Create vector of column names for each quintile.
    quintile_names = [:s1, :s2, :s3, :s4, :s5]

    # Sort the data by decade periods.
    sort!(income_data, :decade)

    # Loop through each quintile and decade and assign appropriate income shares to each region.
    for (d, decade) in enumerate(ssp_years)
        for (q, quintile) in enumerate(quintile_names)

            # Find indices with data for a particular region in a given decade.
            decade_index = findall(income_data.decade .== decade)

            # Allocate those particular income share values and convert from percentage to income shares.
            quintile_shares[d+1, :, q] = income_data[decade_index, quintile] ./ 100
        end
    end

    # Set shares in 2005 to calculated 2015 value (estiamted values start in 2015, but NICE initializes in 2005).
    quintile_shares[1,:,:] .= quintile_shares[2,:,:]

    # Get index for first decade without calculated shares (i.e. after 2105).
    missing_index = findfirst(x->x==0, quintile_shares[:,1,1])

    # Set income shares constant (to 2105 value) for all periods after 2105
    for q = 1:5
        quintile_shares[missing_index:end,:,q] = repeat(quintile_shares[missing_index-1,:,q]', length(nice_years) - missing_index +1)
    end

    # Return cleaned up income shares.
    return quintile_shares
end

#####################################################################################################################
# CALCULATE DAMAGE, CO₂ MITIGATION COST, OR CO₂ TAX BURDEN DISTRIBUTIONS ACROSS AN INDIVIDUAL REGION'S QUINTILES.
#####################################################################################################################
# Description: This function will calculate quintile distribution shares for an individual RICE region based
#              on a provided income elasticity.
#
# Function Arguments:
#
#       elasticity    = Income elasticity of climate damages, CO₂ mitigation costs, CO₂ tax burdens, etc.
#       income_shares = A vector of quintile income shares for a given RICE region.
#--------------------------------------------------------------------------------------------------------------------

function regional_quintile_distribution(elasticity, income_shares)

    # Apply elasticity to quintile income shares.
    scaled_shares = income_shares .^ elasticity

    # Allocate empty array for distribution across quintiles resulting from the elasticity.
    updated_quintile_distribution = zeros(5)

    # Loop through each quintile to calculate updated distribution.
    for q in 1:5
        updated_quintile_distribution[q] = scaled_shares[q] ./ sum(scaled_shares[:])
    end

    return updated_quintile_distribution
end

#####################################################################################################################
# PERFORM A REGRESSION.
#####################################################################################################################
# Description: This function will perform a regression of 'Y' against a user-supplied 'Data' term.
#
# Function Arguments:
#
#       Y    = Dependent variable in regression.
#       Data = Independent data (x-terms) in regression.
#--------------------------------------------------------------------------------------------------------------------

function regress(Y, Data)

    # Get size of data.
    n = length(Data[:,1])
    k = length(Data[1,:])

    # Add column of ones.
    X = [ones(n) Data]

    # Calculate beta coefficient.
    x_prime_inv = inv(X' * X)
    betahat = x_prime_inv * X' * Y

    # Calculate residuals.
    u = Y - X * betahat

    # Calculate hetero var-cov matrix.
    X_eps = zeros(n,k+1)
    for i = 1:n
        X_eps[i,:] = u[i] * X[i,:]
    end

    varhetero = x_prime_inv * X_eps' * X_eps * x_prime_inv

    return betahat, varhetero
end

#####################################################################################################################
# PERFORM A META-REGRESSION OVER MULTIPLE ELASTICITY STUDIES.
#####################################################################################################################
# Description: This function will perform a regression (elasticity vs. ln gdp per capita) across the results of
#              multiple studies to provide the coefficients for a time-varying income elasticity of CO₂ price
#              exposure relationship.
#
# Function Arguments:
#
#       Data       = Independent data (x-terms) in regression coming from multiple studies.
#       slope_type = The type of analysis to carry out. Options are :central, :steeper, :flatter, :percentile
#                       :central    = returns central regression result
#                       :steeper    = returns results for the steepest slope in the 95% confidence interval range
#                       :flatter    = returns results for the shallowest slope in the 95% confidence interval range
#                       :percentile = returns results for a horizontal line corresponding to a given elasticity percentile
#       percentile = Calculate this percentile for the elasticity values (defaults to 0.9 for the 90th percentile).
#       z_val      = z-value for a given confidence interval (defaults to 1.96 for 95% interval)
#--------------------------------------------------------------------------------------------------------------------

function meta_regression(Data; slope_type::Symbol=:central, percentile::Float64=0.90, z_val::Float64=1.96)

    # Sort study data by region and country names.
    sort!(Data, [:Region, :CountryA3])
    n = size(Data, 1)

    # Allocate array to store elastcities.
    elasticities = zeros(n)

    # Isolate variables needed for each run (quintile fuel spending and income shares).
    cost_shares   = [:costS1, :costS2, :costS3, :costS4, :costS5]
    exp_shares = [:e1, :e2, :e3, :e4, :e5]

    # Loop over each study and calculate elasticity with a log-log regression.
    for i = 1:n

        # Take logs of variables
        log_exp_share = log.(Vector(Data[i, exp_shares]))
        log_cost_share= log.(Vector(Data[i, cost_shares]))

        # Run log-log regression and store elasticity.
        B, V   = regress(log_cost_share, log_exp_share)
        elasticities[i] = B[2]
    end

    # Calculate log of per capita GDP.
    Data.log_pcGDP = log.(Data[!,:pcGDP])

    # Using these elasticities, run meta-regression (elasticity vs. ln gdp per capita).
    meta_B, meta_V  = regress(elasticities, Data[!,:log_pcGDP])
    meta_intercept  = meta_B[1]
    meta_slope      = meta_B[2]

    # ------------------------------------------------------------------------------------------------------
    # Calculate flattest and steepest slopes based on 95% confidence interval range.
    # ------------------------------------------------------------------------------------------------------

    # First, need 95% CI around regression lines (requires squared errors & "x_spread" -> (x - mean(x))²)

    # Predicted elasticity values.
    elasticity_hat = meta_intercept .+ meta_slope .* Data[!,:log_pcGDP]

    # Calculate squared errors of elasticities and their mean.
    sq_err = (elasticity_hat .- elasticities) .^2
    mse = mean(sq_err)

    # Calculate squared differences between per capita GDP values and their mean.
    x_spread = (Data[!,:log_pcGDP] .- mean(Data[!,:log_pcGDP])).^2

    # Calculate lower and upper confidence intervals.
    lower_ci = elasticity_hat .- z_val * sqrt.(mse*(1/n .+ x_spread ./ sum(x_spread)))
    upper_ci = elasticity_hat .+ z_val * sqrt.(mse*(1/n .+ x_spread ./ sum(x_spread)))

    # Second, For the steepest line connect maximum point of upper CI and minimum of lower CI.

    # Calculate steeper slope and intercept.
    steeper_slope     = (maximum(upper_ci) - minimum(lower_ci)) / (minimum(Data[!,:log_pcGDP]) - maximum(Data[!,:log_pcGDP]))
    steeper_intercept = maximum(upper_ci) - steeper_slope * minimum(Data[!,:log_pcGDP])

    # Also calculate a shallower slope and intercept.
    flatter_slope     = (minimum(upper_ci) - maximum(lower_ci)) / (maximum(Data[!,:log_pcGDP]) - minimum(Data[!,:log_pcGDP]))
    flatter_intercept = minimum(upper_ci) - flatter_slope * maximum(Data[!,:log_pcGDP])

    # ------------------------------------------------------------------------------------------------------
    # Calculate horizontal lines for a given elasticity percentile.
    # ------------------------------------------------------------------------------------------------------

    # Horizontal line, so slope = 0.
    percentile_slope = 0.0

    # Calculate intercept for a given percentile.
    percentile_intercept = quantile(elasticities, percentile)

    # ------------------------------------------------------------------------------------------------------
    # Return slope and intercept depending on user-specification for analysis type.
    # ------------------------------------------------------------------------------------------------------

    if slope_type == :central
        return meta_intercept, meta_slope
    elseif slope_type == :steeper
        return steeper_intercept, steeper_slope
    elseif slope_type == :flatter
        return flatter_intercept, flatter_slope
    elseif slope_type == :percentile
        return percentile_intercept, percentile_slope
    else
        println("Incorrect regression analysis type selected. Options are :central, :steeper, :flatter, :percentile")
    end
end

#######################################################################################################################
# CALCULATE REGIONAL CO₂ MITIGATION
########################################################################################################################
# Description: This function calculates regional CO₂ mitigation levels as a function of a global carbon tax.  It
#              uses the RICE2010 backstop price values and assumes a carbon tax of $0 in period 1.  If the number of
#              tax values is less than the total number of model time periods, the function assumes full decarbonization
#              (e.g. the tax = the backstop price) for all future periods without a specified tax.
#
# Function Arguments:
#
#       opt_tax_vals:    A vector of global carbon tax values being optimized.
#       backstop_prices: The regional backstop prices from RICE2010 (units = dollars).
#       theta2:          The exponent on the abatement cost function (defaults to RICE2010 value).
#----------------------------------------------------------------------------------------------------------------------

function mu_from_tax(opt_tax_vals::Array{Float64,1}, backstop_prices::Array{Float64,2}, ; theta2::Float64=2.8)

    # Set first period to $0, then optimized tax values, and then remaining values to maximum backstop price (i.e. full decarbonization).
    full_CO₂_tax = [0.0; opt_tax_vals; maximum(backstop_prices, dims=2)[(length(opt_tax_vals)+2):end]]

    # Calculate regional mitigation rates from the full tax vector.
    regional_CO₂_abatement = min.((max.(((full_CO₂_tax ./ backstop_prices) .^ (1 / (theta2 - 1.0))), 0.0)), 1.0)

    return full_CO₂_tax, regional_CO₂_abatement
end

#######################################################################################################################
# CREATE OBJEVTIVE FUNCTION FOR NICE WITH REVENUE RECYCLING
########################################################################################################################
# Description: This function creates an objective function that uses an instance of NICE with revenue recycling and
#              user-specified settings. The objective function will take in a vector of global carbon tax values and
#              returns the total economic welfare generated by that specifc climate policy.
#
# Function Arguments:
#
#       m:                 An instance of NICE with revenue recycling (type = a Mimi model).
#       revenue_recycling: A check for whether or not to recycle CO₂ tax revenue (true = recycle, false = no recycling).
#----------------------------------------------------------------------------------------------------------------------

function construct_nice_recycle_objective(m::Model; revenue_recycling::Bool=true)

    # Extract backstop prices from model and convert to dollars.
    run(m)
    rice_backstop_prices = m[:emissions, :pbacktime] .* 1000

    # Find number of timesteps across model time horizon.
    n_steps = length(dim_keys(m, :time))

    # Pre-allocate matrix to store optimal tax and mitigation rates.
    optimal_CO₂_tax        = zeros(n_steps)
    optimal_CO₂_mitigation = zeros(n_steps, 12)

    # Create a function to optimize user-specified model for (i) revenue recycling and (ii) reference case.
    nice_recycle_objective =

        if revenue_recycling == true

            # Create revenue recycling version of function.
            function(opt_tax::Array{Float64,1})
                # Calculate mitigation rate resulting from global carbon tax.
                optimal_CO₂_tax[:], optimal_CO₂_mitigation[:,:] = mu_from_tax(opt_tax, rice_backstop_prices)
                # Update CO₂ mitigation rate and global CO₂ tax in NICE.
                update_param!(m, :MIU, optimal_CO₂_mitigation)
                update_param!(m, :global_carbon_tax, optimal_CO₂_tax)
                run(m)
                # Return total welfare.
                return m[:nice_welfare, :welfare]
            end

        else

            # Create reference version of function without revenue recycling.
            function(opt_tax::Array{Float64,1})
                # Calculate mitigation rate resulting from global carbon tax.
                optimal_CO₂_tax[:], optimal_CO₂_mitigation[:,:] = mu_from_tax(opt_tax, rice_backstop_prices)
                # Update CO₂ mitigation rate in NICE.
                # Set tax in recycling component to $0 yielding no tax revenue (i.e. switch off revenue recycling).
                update_param!(m, :MIU, optimal_CO₂_mitigation)
                update_param!(m, :global_carbon_tax, zeros(n_steps))
                run(m)
                # Return total welfare.
                return m[:nice_welfare, :welfare]
            end
        end

    # Return the objective function.
    return nice_recycle_objective
end

#######################################################################################################################
# CALCULATE GLOBAL CO₂ MITIGATION LEVELS
#######################################################################################################################
# Description: This function takes an instance of NICE with regional CO₂ mitigation and returns the resulting global
#              reduction in CO₂ emissions relative to a baseline case without a climate policy.
#
# Function Arguments:
#
#       m_policy: An instance of NICE with CO₂ policy (type = Mimi model).
#       m_bau:    An instance of NICE with 0% mitigation (no CO₂ policy) for all regions and years (type = Mimi model).
#----------------------------------------------------------------------------------------------------------------------

function get_global_mitigation(m_policy::Model, m_bau::Model)

    # Calculate global industrial CO₂ emissions for NICE with a particular mitigation policy.
    global_emissions_policy = sum(m_policy[:emissions, :EIND], dims=2)

    # Calculate global industrial CO₂ emissions in the bau case with no CO₂ policy.
    global_emissions_base = sum(m_bau[:emissions, :EIND], dims=2)

    # Calculate change in global CO₂ mitigation levels.
    global_mitigation_rates = ((global_emissions_base .- global_emissions_policy) ./ global_emissions_base)

    return global_mitigation_rates
end

#######################################################################################################################
# CREATE AN INSTANCE OF FAIR TO USE IN OPTIMIZATION WITH NICE
#######################################################################################################################
# Description: This function creates a version of FAIR v1.3 that can be used in an optimization with NICE. It removes
#              all non-CO₂ radiative forcing components (instead it will take in the endgoenous CO₂ emissions and
#              exogenous non-CO₂ forcings from NICE). It also modifies FAIR to run from 2005-2595 (rather than its
#              default of 1765-2500).
#----------------------------------------------------------------------------------------------------------------------

function create_fair_for_opt()

    #----- Initial Data and Index Prep -----#

    # Set the default start and end years for FAIR.
    fair_start_year = 1765
    fair_end_year   = 2500

    # Set the default start and end years for NICE.
    nice_start_year = 2005
    nice_end_year   = 2595

    # Set annual time steps for NICE and calculate number of annual years.
    nice_annual_years  = collect(nice_start_year:nice_end_year)
    n_annual_nice      = length(nice_annual_years)

    # Calculate FAIR index for year 2005.
    fair_2005_index = findfirst(x -> x == 2005, fair_start_year:fair_end_year)

    # Load RCP4.5 concentrations and get N₂O concentrations (needed for CO₂ forcing calculations in FAIR).
    rcp45_concentrations = DataFrame(load(joinpath(@__DIR__, "..", "data", "RCP45_concentrations.csv"), skiplines_begin=37))
    rcp45_n2o_raw = rcp45_concentrations.N2O

    # Isolate RCP4.5 N₂O data for 2005-2500 (start of NICE to end of RCP data). Then extend to 2595 assuming concentrations remain fixed at 2500 value.
    rcp45_n2o = zeros(n_annual_nice)
    rcp45_n2o[1:length(2005:2500)] = rcp45_n2o_raw[fair_2005_index:end]
    rcp45_n2o[length(2005:2500):end] .= rcp45_n2o_raw[end]

    #----- Create Modified Version of FAIR for NICE -----#

    # Load and run a default version of FAIR.
    fair = MimiFAIR13.create_fair(rcp_scenario="RCP45")
    run(fair)

    # Get initial CO₂ concentration and CO₂ stock perturbation (for NICE years).
    Cacc_2005 = fair[:co2_cycle, :Cacc][fair_2005_index]
    CO2_2004  = fair[:co2_cycle, :C][(fair_2005_index-1)]

    # Remove all FAIR components other than CO₂, radiative forcing, and temperature (use NICE exogenous forcing values).
    delete!(fair, :ch4_cycle)
    delete!(fair, :n2o_cycle)
    delete!(fair, :other_ghg_cycles)
    delete!(fair, :ch4_rf)
    delete!(fair, :n2o_rf)
    delete!(fair, :other_ghg_rf)
    delete!(fair, :trop_o3_rf)
    delete!(fair, :strat_o3_rf)
    delete!(fair, :aerosol_direct_rf)
    delete!(fair, :aerosol_indirect_rf)
    delete!(fair, :bc_snow_rf)
    delete!(fair, :landuse_rf)
    delete!(fair, :contrails_rf)

    # Reset FAIR time dimension to NICE years (annual version).
    set_dimension!(fair, :time, nice_annual_years)

    # Set placeholder value for CO₂ emissions.
    update_param!(fair, :E, ones(n_annual_nice), update_timesteps=true)

    #Set initial Cacc perturbation.
    update_param!(fair, :Cacc_0, Cacc_2005)

    # Set initial (2004) atmospheric CO₂ concentration.
  #  set_param!(fair, :CO₂_0, CO2_2004)
    set_param!(fair, :co2_cycle, :CO₂_0, :co2cycle_2004val, CO2_2004)
#set_param!(m::Model, comp_name::Symbol, param_name::Symbol, ext_param_name::Symbol, val::Any)

    # Set RCP4.5 N₂O concentrations for CO₂ radiative forcing calculations.
    set_param!(fair, :N₂O, rcp45_n2o)

    # Set placeholder for NICE exogenous forcings.
    update_param!(fair, :F_exogenous, ones(n_annual_nice), update_timesteps=true)

    # Set all other radiative forcings to zero.
    set_param!(fair, :F_volcanic, zeros(n_annual_nice))
    set_param!(fair, :F_solar, zeros(n_annual_nice))
    set_param!(fair, :F_CH₄, zeros(n_annual_nice))
    set_param!(fair, :F_CH₄_H₂O, zeros(n_annual_nice))
    set_param!(fair, :F_N₂O, zeros(n_annual_nice))
    set_param!(fair, :F_other_ghg, zeros(n_annual_nice, 28))
    set_param!(fair, :F_trop_O₃, zeros(n_annual_nice))
    set_param!(fair, :F_strat_O₃, zeros(n_annual_nice))
    set_param!(fair, :F_aerosol_direct, zeros(n_annual_nice))
    set_param!(fair, :F_aerosol_indirect, zeros(n_annual_nice))
    set_param!(fair, :F_bcsnow, zeros(n_annual_nice))
    set_param!(fair, :F_landuse, zeros(n_annual_nice))
    set_param!(fair, :F_contrails, zeros(n_annual_nice))

    # Return modified version of FAIR.
    return fair
end

#######################################################################################################################
# LINEARLY INTERPOLATE MODEL RESULTS TO ANNUAL VALUES
#######################################################################################################################
# Description: This function uses linear interpolation to create an annual time series from the model results.
#
# Function Arguments:
#
#       data    = The non-annual model results to be interpolated
#       spacing = Length of time steps between model output.
#----------------------------------------------------------------------------------------------------------------------

function interpolate_to_annual(data, spacing)

    # Create an interpolation object for the data (assume first and last points are end points, e.g. no interpolation beyond support).
    interp_linear = interpolate(data, BSpline(Linear()))

    # Create points to interpolate for (based on spacing term).
    interp_points = collect(1:(1/spacing):length(data))

    # Carry out interpolation.
    return interp_linear[interp_points]
end

function construct_fair_recycle_objective(nice::Model, n_loops::Int; revenue_recycling::Bool=true)

    # Run user-specified instance of NICE to obtain some values.
    run(nice)

    # Extract backstop prices from NICE and convert to dollars.
    nice_backstop_prices = nice[:emissions, :pbacktime] .* 1000

    # Extract exogenous radiative forcing scenario from NICE and convert to annual values.
    annual_nice_exog_rf = interpolate_to_annual(nice[:radiativeforcing, :forcoth], 10)

    # Get an instance of FAIR to use in iterative optimization.
    fair = create_fair_for_opt()

    # Set FAIR exogenous forcing to annualized NICE exogenous forcing scenario.
    update_param!(fair, :F_exogenous, annual_nice_exog_rf)

    # Find number of timesteps across NICE model time horizon.
    n_nice_decadal = length(dim_keys(nice, :time))

    # Get decadal NICE indices for annual FAIR years (2005-2595)
    nice_decadal_indices = findall((in)(collect(2005:10:2595)), (collect(2005:2595)))

    # Pre-allocate matrix to store optimal tax and mitigation rates at 10-year time steps.
    optimal_CO₂_tax        = zeros(n_nice_decadal)
    optimal_CO₂_mitigation = zeros(n_nice_decadal, 12)

    # Create a function to optimize user-specified model for (i) revenue recycling and (ii) reference case.
    fair_recycle_objective =

        if revenue_recycling == true

            # Create revenue recycling version of function.
            function(opt_tax::Array{Float64,1})

                # Calculate mitigation rate resulting from global carbon tax.
                optimal_CO₂_tax[:], optimal_CO₂_mitigation[:,:] = mu_from_tax(opt_tax, nice_backstop_prices)

                # Update CO₂ mitigation rate and global CO₂ tax in NICE.
                update_param!(nice, :MIU, optimal_CO₂_mitigation)
                update_param!(nice, :global_carbon_tax, optimal_CO₂_tax)
                run(nice)

                for loops = 1:n_loops

                    # Run FAIR with annualized NICE CO₂ emissions to calculate temperature.
                    update_param!(fair, :E, interpolate_to_annual(nice[:emissions, :E], 10))
                    run(fair)

                    # Set FAIR temperature projections in NICE.
                    set_param!(nice, :TATM, fair[:temperature, :T][nice_decadal_indices])
                    run(nice)
                end

                # Return total welfare.
                return nice[:nice_welfare, :welfare]
            end

        else

            # Create reference version of function without revenue recycling.
            function(opt_tax::Array{Float64,1})

                # Calculate mitigation rate resulting from global carbon tax.
                optimal_CO₂_tax[:], optimal_CO₂_mitigation[:,:] = mu_from_tax(opt_tax, nice_backstop_prices)

                # Update CO₂ mitigation rate in NICE.
                # Set tax in recycling component to $0 yielding no tax revenue (i.e. switch off revenue recycling).
                update_param!(nice, :MIU, optimal_CO₂_mitigation)
                update_param!(nice, :global_carbon_tax, zeros(n_nice_decadal))
                run(nice)

                for loops = 1:n_loops

                    # Run FAIR with annualized NICE CO₂ emissions to calculate temperature.
                    update_param!(fair, :E, interpolate_to_annual(nice[:emissions, :E], 10))
                    run(fair)

                    # Set FAIR temperature projections in NICE.
                    set_param!(nice, :TATM, fair[:temperature, :T][nice_decadal_indices])
                    run(nice)
                end

                # Return total welfare.
                return nice[:nice_welfare, :welfare]
            end
        end

    # Return the objective function.
    return fair_recycle_objective
end

#######################################################################################################################
# CREATE RESULT DIRECTORIES AND SAVE SPECIFIC MODEL OUTPUT
#######################################################################################################################
# Description: This function will create a folder directory to store results (dividing model output by global,
#              regional, and quintile levels) for either the revenue recyling or reference case without recycling.
#
# Function Arguments:
#
#       m_policy:                 An instance of NICE with CO₂ policy (type = Mimi model).
#       m_bau:                    An instance of NICE with 0% mitigation (no CO₂ policy) for all regions and years (type = Mimi model).
#       opt_tax:                  The optimal tax policy vector resulting from the second-stage local optimization.
#       opt_tax_global_algorithm: The optimal tax policy vector resulting from the first-stage global optimization (just saved for reference).
#       output_directory:         The directory path to the results folder where a particular set of model output will be saved.
#       revenue_recycling:        A check for whether or not the results recycle CO₂ tax revenue (true = recycle, false = no recycling).
#----------------------------------------------------------------------------------------------------------------------

function save_nice_recycle_results(m_policy::Model, m_bau::Model, opt_tax::Array{Float64,1}, opt_tax_global_algorithm::Array{Float64,1}, output_directory::String; revenue_recycling::Bool=true)

    # Make subdirectory folders to store results with and without revenue recycling.
    if revenue_recycling == true

        global_path   = joinpath(output_directory, "revenue_recycling", "global_output")
        regional_path = joinpath(output_directory, "revenue_recycling", "regional_output")
        quintile_path = joinpath(output_directory, "revenue_recycling", "quintile_output")

        mkpath(global_path)
        mkpath(regional_path)
        mkpath(quintile_path)

    else

        global_path   = joinpath(output_directory, "no_revenue_recycling", "global_output")
        regional_path = joinpath(output_directory, "no_revenue_recycling", "regional_output")
        quintile_path = joinpath(output_directory, "no_revenue_recycling", "quintile_output")

        mkpath(global_path)
        mkpath(regional_path)
        mkpath(quintile_path)
    end

    # Save Global Output.
    save(joinpath(global_path, "global_co2_mitigation.csv"), DataFrame(get_global_mitigation(m_policy, m_bau), :auto))
    save(joinpath(global_path, "temperature.csv"), DataFrame(temperature=m_policy[:climatedynamics, :TATM]))
    save(joinpath(global_path, "rough_carbon_tax_global_algorithm.csv"), DataFrame(global_algorithm_tax=opt_tax_global_algorithm))
    save(joinpath(global_path, "carbon_tax.csv"), DataFrame(tax=opt_tax))
    save(joinpath(global_path, "global_emissions_with_landuse.csv"), DataFrame(emissions=m_policy[:emissions, :E]))
    save(joinpath(global_path, "co2_concentration_ppm.csv"), DataFrame(CO2=m_policy[:co2cycle, :MAT] ./ 2.13))

    # Save Regional Output.
    save(joinpath(regional_path, "gross_output.csv"), DataFrame(m_policy[:grosseconomy, :YGROSS], :auto))
    save(joinpath(regional_path, "nice_net_output.csv"), DataFrame(m_policy[:nice_neteconomy, :Y], :auto))
    save(joinpath(regional_path, "consumption.csv"), DataFrame(m_policy[:nice_neteconomy, :C], :auto))
    save(joinpath(regional_path, "population.csv"), DataFrame(m_policy[:nice_neteconomy, :l], :auto))
    save(joinpath(regional_path, "industrial_co2_emissions.csv"), DataFrame(m_policy[:emissions, :EIND], :auto))
    save(joinpath(regional_path, "damage_cost_share.csv"), DataFrame(m_policy[:nice_neteconomy, :DAMFRAC], :auto))
    save(joinpath(regional_path, "abatement_cost_share.csv"), DataFrame(m_policy[:nice_neteconomy, :ABATEFRAC], :auto))
    save(joinpath(regional_path, "regional_co2_mitigation.csv"), DataFrame(m_policy[:emissions, :MIU], :auto))
    save(joinpath(regional_path, "co2_tax_revenue.csv"), DataFrame(m_policy[:nice_recycle, :tax_revenue], :auto))
    save(joinpath(regional_path, "co2_income_elasticity.csv"), DataFrame(m_policy[:nice_recycle, :CO₂_income_elasticity], :auto))

    # Save Quintile Output.
    save(joinpath(quintile_path, "co2_tax_distribution_Q1.csv"), DataFrame(m_policy[:nice_recycle, :carbon_tax_dist][:,:,1], :auto))
    save(joinpath(quintile_path, "co2_tax_distribution_Q2.csv"), DataFrame(m_policy[:nice_recycle, :carbon_tax_dist][:,:,2], :auto))
    save(joinpath(quintile_path, "co2_tax_distribution_Q3.csv"), DataFrame(m_policy[:nice_recycle, :carbon_tax_dist][:,:,3], :auto))
    save(joinpath(quintile_path, "co2_tax_distribution_Q4.csv"), DataFrame(m_policy[:nice_recycle, :carbon_tax_dist][:,:,4], :auto))
    save(joinpath(quintile_path, "co2_tax_distribution_Q5.csv"), DataFrame(m_policy[:nice_recycle, :carbon_tax_dist][:,:,5], :auto))
    save(joinpath(quintile_path, "base_pc_consumption_Q1.csv"), DataFrame(m_policy[:nice_recycle, :qc_base][:,:,1], :auto))
    save(joinpath(quintile_path, "base_pc_consumption_Q2.csv"), DataFrame(m_policy[:nice_recycle, :qc_base][:,:,2], :auto))
    save(joinpath(quintile_path, "base_pc_consumption_Q3.csv"), DataFrame(m_policy[:nice_recycle, :qc_base][:,:,3], :auto))
    save(joinpath(quintile_path, "base_pc_consumption_Q4.csv"), DataFrame(m_policy[:nice_recycle, :qc_base][:,:,4], :auto))
    save(joinpath(quintile_path, "base_pc_consumption_Q5.csv"), DataFrame(m_policy[:nice_recycle, :qc_base][:,:,5], :auto))
    save(joinpath(quintile_path, "post_damage_abatement_pc_consumption_Q1.csv"), DataFrame(m_policy[:nice_recycle, :qc_post_damage_abatement][:,:,1], :auto))
    save(joinpath(quintile_path, "post_damage_abatement_pc_consumption_Q2.csv"), DataFrame(m_policy[:nice_recycle, :qc_post_damage_abatement][:,:,2], :auto))
    save(joinpath(quintile_path, "post_damage_abatement_pc_consumption_Q3.csv"), DataFrame(m_policy[:nice_recycle, :qc_post_damage_abatement][:,:,3], :auto))
    save(joinpath(quintile_path, "post_damage_abatement_pc_consumption_Q4.csv"), DataFrame(m_policy[:nice_recycle, :qc_post_damage_abatement][:,:,4], :auto))
    save(joinpath(quintile_path, "post_damage_abatement_pc_consumption_Q5.csv"), DataFrame(m_policy[:nice_recycle, :qc_post_damage_abatement][:,:,5], :auto))
    save(joinpath(quintile_path, "post_recycle_pc_consumption_Q1.csv"), DataFrame(m_policy[:nice_recycle, :qc_post_recycle][:,:,1], :auto))
    save(joinpath(quintile_path, "post_recycle_pc_consumption_Q2.csv"), DataFrame(m_policy[:nice_recycle, :qc_post_recycle][:,:,2], :auto))
    save(joinpath(quintile_path, "post_recycle_pc_consumption_Q3.csv"), DataFrame(m_policy[:nice_recycle, :qc_post_recycle][:,:,3], :auto))
    save(joinpath(quintile_path, "post_recycle_pc_consumption_Q4.csv"), DataFrame(m_policy[:nice_recycle, :qc_post_recycle][:,:,4], :auto))
    save(joinpath(quintile_path, "post_recycle_pc_consumption_Q5.csv"), DataFrame(m_policy[:nice_recycle, :qc_post_recycle][:,:,5], :auto))

end


text:
@generated function __new__(T, args...)
  quote
    Base.@_inline_meta
    $(Expr(:new, :T, [:(args[$i]) for i = 1:length(args)]...))
  end
end

@generated function __splatnew__(T, args)
  quote
    Base.@_inline_meta
    $(Expr(:splatnew, :T, :args))
  end
end

@partial Basic function __new__(::Type{T}, xs...) where T
  Partial{T}(Any[i > length(xs) ? Union{} : xs[i] for i in 1:length(fieldnames(T))])
end

@partial Basic __new__(t::Type{T}, xs::Const...) where T =
  abstract(Basic(), Const(__new__), t, xs...)

@abstract Basic function __new__(::Type{T}, xs::Const...) where T
  if T.mutable
    Partial{T}(Any[i > length(xs) ? Union{} : xs[i] for i in 1:length(fieldnames(T))])
  else
    Const(__new__(T, map(x -> x.value, xs)...))
  end
end

@partial Basic function __splatnew__(::Type{T}, xs::Const{<:Tuple}) where T
  Const(__splatnew__(T, xs.value))
end

@partial Basic function __splatnew__(::Type{T}, xs::Partial{<:Tuple}) where T
  Partial{T}(Any[i > length(xs.value) ? Union{} : xs.value[i] for i in 1:length(fieldnames(T))])
end

abstract(::Basic, ::AType{typeof(tuple)}, xs::Type...) = Tuple{xs...}

@abstract Basic tuple(xs...) = ptuple(xs...)

@abstract Basic getindex(x::Const{<:Tuple}, i::Const{<:Integer}) =
  Const(x.value[i.value])

@abstract Basic getindex(x::Partial{<:Tuple}, i::Const{<:Integer}) =
  x.value[i.value]

@abstract Basic getindex(xs::Tuple, i::Const{<:Integer}) = widen(xs).parameters[i.value]

@abstract Basic length(xs::Partial{<:Tuple}) = Const(length(xs.value))

@partial Basic function setfield!(x::Partial{T}, name::Const{Symbol}, s) where T
  i = findfirst(f -> f == name.value, fieldnames(T))
  x.value[i] = s
  x
end

@partial Basic function getfield(x::Partial{T}, name::Const{Symbol}) where T
  i = findfirst(f -> f == name.value, fieldnames(T))
  x.value[i]
end

@partial Basic function getfield(x::Partial{T}, i::Const{<:Integer}) where T
  x.value[i.value]
end

@partial Basic function getfield(x::Const{T}, i::Const{<:Integer}) where T
  Const(x.value[i.value])
end

@partial Basic getfield(x::Const, name::Const{Symbol}) =
  Const(getfield(x.value, name.value))

function mutate(::Basic, cx::MutCtx, ::AType{typeof(setfield!)}, x::Partial{T}, name::Const{Symbol}, s) where T
  i = findfirst(f -> f == name.value, fieldnames(T))
  S = x.value[i]
  if !_issubtype(s, S)
    x.value[i] = _union(S, s)
    visit!(cx, x)
  end
  x
end

function mutate(::Basic, cx::MutCtx, ::AType{typeof(getfield)}, x::Partial{T}, name) where T
  edge!(cx, x)
  i = findfirst(f -> f == name.value, fieldnames(T))
  x.value[i]
end

# Dictionaries

@partial Basic Dict() = Partial{Dict{Any,Any}}(Dict())

@partial Basic function setindex!(x::Partial{Dict{K,V}}, s::V, name::Const{<:K}) where {K,V}
  x.value[name.value] = s
  x
end

@partial Basic function getindex(x::Partial{Dict{K,V}}, name::Const{<:K}) where {K,V}
  x.value[name.value]
end

function mutate(::Basic, cx::MutCtx, ::AType{typeof(setindex!)}, x::Partial{Dict{K,V}}, s::AType{<:V}, name::Const{<:K}) where {K,V}
  T = get(x.value, name.value, Union{})
  if !_issubtype(s, T)
    visit!(cx, x)
    x.value[name.value] = _union(T, s)
  end
  return x
end

function mutate(::Basic, cx::MutCtx, ::AType{typeof(getindex)}, x::Partial{Dict{K,V}}, name::Const{<:K}) where {K,V}
  edge!(cx, x)
  get(x.value, name.value, Union{})
end

@abstract Basic (===)(a::Const{T}, b::Const{T}) where T =
  Const(a.value === b.value)

@abstract Basic (===)(a::AType{T}, b::AType{T}) where T = Bool
@abstract Basic (===)(a, b) = Const(false)

@abstract Basic isempty(x::Const) = Const(isempty(x.value))


text:
"""
    get_atoms(met::Metabolite)::MetaboliteFormula

Simple wrapper for getting the atom dictionary count out of a `Metabolite`.
"""
get_atoms(met::Metabolite)::MetaboliteFormula = _maybemap(_parse_formula, met.formula)