text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#
# Contributed by Adam Beckmeyer

using Printf

const COUNTSFOR = ("ggt", "ggta", "ggtatt", "ggtattttaatt", "ggtattttaatttatagt")

# This is not reading the file in line-by-line. Doing so is impossible
# in Julia without using FFI because all io is automatically
# buffered. Thus despite not appearing to be buffered reads, all are.
function getseq3(io)
    # First read all the file preceding the 3rd sequence.
    for _=1:3
        readuntil(io, '>')
    end
    readline(io)
    # Then read the third sequence in.
    buf = read(io)

    # In place, remove all newlines from buf and encode nucleotides
    # using only the last 2 bits in each byte.
    i = 1
    for c in buf
        if c != 0x0a # '\n'
            # Gives a -> 0x00, c -> 0x01, g -> 0x03, t -> 0x02
            @inbounds buf[i] = c >> 1 & 0x03
            i += 1
        end
    end
    resize!(buf, i - 1)
end

# Decoding a single encoded nucleotide results in a byte 'a', 'c', 'g' or 't'
decode(c) = c == 0x00 ? 0x61 : c == 0x01 ? 0x63 : c == 0x03 ? 0x67 : 0x74

# Decoding a UInt32 or UInt64 results in a string. This function
# creates the string from the last n encoded nucleotides in c.
function decode(c, n)
    buf = Vector{UInt8}(undef, n)
    for i=1:n
        @inbounds buf[n+1-i] = decode(c % UInt8 & 0x03)
        c = c >> 2
    end
    String(buf)
end

# Count the number of each subsequence of length n present in seq and
# store the result in hashtable d.
function count_subseqs!(d::Dict{T}, seq, n) where T
    # This results in a number where the last 2n bits are 1 and the rest are 0.
    mask = (1 << 2n - 1) % T
    key = zero(T)

    # Fill key with the first n-1 encoded nucleotides
    for i=1:n-1
        key = key << 2 | @inbounds seq[i]
    end

    # Slide key over the remainder of the sequence. Because of mask,
    # only the last n nucleotides are encoded in key.
    for i=n:length(seq)
        key = mask & (key << 2 | @inbounds seq[i])
        inc!(d, key)
    end
end

# If k exists in dictionary, increment it by one. Otherwise set the
# value of k to 1. If julia ever gets an update!/modify! function for
# changing hashtable values with one lookup, this function can be
# removed.
@inline function inc!(d::Dict{T,Int32}, k::T) where T
    index = Base.ht_keyindex2!(d, k)

    if index > 0
        # Positive index means the key already exists in the dictionary
        @inbounds d.vals[index] += 1
    else
        # Negative index means the key doesn't yet exist
        @inbounds Base._setindex!(d, 1, k, -index)
    end
end

# Define a fast hashing function for UInt keys encoded with the
# sequence. Defining a Base function on a Base type is "type-piracy"
# and should be avoided in production code. Instead, Base.hash could
# be defined on a wrapper struct. For this workload, this simple
# xor/bitshift hash speeds things up by over 20% compared to Julia's
# default hash.
Base.hash(x::Unsigned)::UInt = x ⊻ x >> 7

# Create a sorted array with the number of appearances of each
# sequence of length n.
function frequency_table(seq, n)
    d = Dict{UInt32,Int32}()
    count_subseqs!(d, seq, n)

    total = sum(values(d))
    counts = [decode(k, n) => get(d, k, 0) / total for k in keys(d)]
    sort!(counts; lt=freq_isless, rev=true)
end

# This function is used for sorting the frequency table generated by
# the function above. Primary order is by frequency which is second
# in pair. Secondary order is reverse order of string which is first
# in pair.
freq_isless(a, b) = a[2] == b[2] ? isless(b[1], a[1]) : isless(a[2], b[2])

# Count the number of times subseq appears in seq by creating a
# hashtable of all sequences of the same length as subseq appearing in
# seq.
function count_sequence(seq, subseq)
    d = length(subseq) < 16 ? Dict{UInt32,Int32}() : Dict{UInt,Int32}()
    count_subseqs!(d, seq, length(subseq))

    # Construct the key needed for accessing the count of subseq
    key = 0 % UInt
    for c in subseq
        key = key << 2 | (c >> 1 & 0x03)
    end

    get(d, key, 0)
end

function main(io, out)
    seq = getseq3(io)

    # First spawn threads for counting occurrences of subseqs in seq
    # since they take longer than the frequency tables.
    counts = Vector{Int32}(undef, 5)
    counts_task = Threads.@spawn @sync for i=1:5
        Threads.@spawn @inbounds counts[i] =
            count_sequence(seq, codeunits(COUNTSFOR[i]))
    end

    # Because of Julia's threading overhead, it's faster to just
    # calculate both the 1 and 2 nucleotide frequency tables on the
    # main thread while counts_task continues in background.
    for (a, b) in frequency_table(seq, 1)
        @printf(out, "%s %2.3f\n", uppercase(a), 100b)
    end
    println(out)

    for (a, b) in frequency_table(seq, 2)
        @printf(out, "%s %2.3f\n", uppercase(a), 100b)
    end
    println(out)

    wait(counts_task)
    @inbounds for i=1:5
        println(out, counts[i], '\t', uppercase(COUNTSFOR[i]))
    end
end

fileName = length(ARGS) > 0 ? ARGS[1] : "25000_in"
f = open(fileName, "r")
isinteractive() || main(f, stdout)
close(f)

max_stars_repo_name: bpecsek/Programming-Language-Benchmarks

max_issues_repo_name: bpecsek/Programming-Language-Benchmarks

max_forks_repo_name: bpecsek/Programming-Language-Benchmarks


text:
module NORBiP

using JuMP
import SCIP

import Polyhedra

using LinearAlgebra: dot, ⋅

include("type_definitions.jl")
include("dual_polyhedron.jl")
include("algorithms.jl")
include("cutting_plane.jl")
include("single_vertex_heuristic.jl")
include("adversarial.jl")
include("valid_inequalities.jl")

# type piracy for the greater good
Base.string(::NearOptimalMethod) = string(typeof(mt))
Base.string(mt::LazyBatched) = "lazy_batched_$(mt.batch_size)"
Base.string(::ExtendedMethod) = "extended"
Base.string(::LazyExtended) = "lazy"
Base.string(mt::SingleVertexHeuristic) = "single_heuristic_$(mt.batch_size)"
Base.string(mt::RandomizedSingleVertexHeuristic) = "random_sv_$(mt.batch_size)"
Base.string(::BilinearTerms) = "bilinear"

end # module

ext: jl


text:
include("utils.jl")

@testset "Torus" begin
    M = Torus(2)
    @testset "Torus Basics" begin
        @test representation_size(M) == (2,)
        @test manifold_dimension(M) == 2
        @test !is_manifold_point(M, 9.)
        @test_throws DomainError is_manifold_point(M, 9., true)
        @test !is_manifold_point(M, [9.;9.])
        @test_throws DomainError is_manifold_point(M, [9. 9.], true)
        @test !is_tangent_vector(M, [9.;9.], 0.)
        @test_throws DomainError is_tangent_vector(M, 9., 0., true)
        @test !is_tangent_vector(M, [9.;9.], [0.;0.])
        @test_throws DomainError is_tangent_vector(M, 9., 0., true)
        @test injectivity_radius(M) ≈ π
        x = [1.0, 2.0]
        y = [-1.0, 2.0]
        z = [0., 0.]
        basis_types = (ArbitraryOrthonormalBasis(),)
        test_manifold(M,
            [x, y, z],
            test_forward_diff = false,
            test_reverse_diff = false,
            test_vector_spaces = true,
            test_project_tangent = false,
            test_musical_isomorphisms = true,
            test_vector_transport = false,
            basis_types_to_from = basis_types,
        )
    end
end


text:
module AWSMetadata

using Base64
using ..AWSExceptions
using GitHub
using HTTP
using JSON
using Mocking
using OrderedCollections: LittleDict, OrderedDict

const services_path = joinpath(@__DIR__, "AWSServices.jl")

include(joinpath("api_generation", "utilities.jl"))
include(joinpath("api_generation", "high_level.jl"))
include(joinpath("api_generation", "low_level.jl"))

"""
    parse_aws_metadata()

Generate low and high level wrappers for each AWS Service based on their definitions in the
[aws-sdk-js GitHub Repository](https://github.com/aws/aws-sdk-js/tree/master/apis).

Low level wrappers are written into `src/AWSServices.jl`, while high level wrappers API
wrappers are written into their respective files in `src/services/{service}.jl`.
"""
function parse_aws_metadata()
    auth = GitHub.authenticate(ENV["GITHUB_AUTH"])
    repo_name = "aws/aws-sdk-js"

    service_files = _get_service_files(repo_name, auth)

    _generate_low_level_wrappers(service_files, auth)
    _generate_high_level_wrapper(service_files, auth)
    return nothing
end

end


text:
using SpinDoctor
using Documenter

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

makedocs(;
    modules = [SpinDoctor],
    authors = "Syver Døving Agdestein <syverda@icloud.com> and contributors",
    repo = "https://github.com/agdestein/SpinDoctor.jl/blob/{commit}{path}#{line}",
    sitename = "SpinDoctor.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://agdestein.github.io/SpinDoctor.jl",
        assets = String[],
    ),
    pages = ["Home" => "index.md"],
)

deploydocs(; repo = "github.com/agdestein/SpinDoctor.jl", devbranch = "main")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct VirtualMachineImagesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets a virtual machine image.
Param: location::String (required)
Param: publisherName::String (required)
Param: offer::String (required)
Param: skus::String (required)
Param: version::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualMachineImage
"""
function virtualMachineImagesGet(_api::VirtualMachineImagesApi, location::String, publisherName::String, offer::String, skus::String, version::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualMachineImage, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/publishers/{publisherName}/artifacttypes/vmimage/offers/{offer}/skus/{skus}/versions/{version}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "publisherName", publisherName)  # type String
    Swagger.set_param(_ctx.path, "offer", offer)  # type String
    Swagger.set_param(_ctx.path, "skus", skus)  # type String
    Swagger.set_param(_ctx.path, "version", version)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets a list of all virtual machine image versions for the specified location, publisher, offer, and SKU.
Param: location::String (required)
Param: publisherName::String (required)
Param: offer::String (required)
Param: skus::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: filter::String
Param: top::Int32
Param: orderby::String
Return: Vector{VirtualMachineImageResource}
"""
function virtualMachineImagesList(_api::VirtualMachineImagesApi, location::String, publisherName::String, offer::String, skus::String, api_version::String, subscriptionId::String; filter=nothing, top=nothing, orderby=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Vector{VirtualMachineImageResource}, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/publishers/{publisherName}/artifacttypes/vmimage/offers/{offer}/skus/{skus}/versions", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "publisherName", publisherName)  # type String
    Swagger.set_param(_ctx.path, "offer", offer)  # type String
    Swagger.set_param(_ctx.path, "skus", skus)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "\$orderby", orderby)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets a list of virtual machine image offers for the specified location and publisher.
Param: location::String (required)
Param: publisherName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Vector{VirtualMachineImageResource}
"""
function virtualMachineImagesListOffers(_api::VirtualMachineImagesApi, location::String, publisherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Vector{VirtualMachineImageResource}, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/publishers/{publisherName}/artifacttypes/vmimage/offers", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "publisherName", publisherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets a list of virtual machine image publishers for the specified Azure location.
Param: location::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Vector{VirtualMachineImageResource}
"""
function virtualMachineImagesListPublishers(_api::VirtualMachineImagesApi, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Vector{VirtualMachineImageResource}, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/publishers", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets a list of virtual machine image SKUs for the specified location, publisher, and offer.
Param: location::String (required)
Param: publisherName::String (required)
Param: offer::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Vector{VirtualMachineImageResource}
"""
function virtualMachineImagesListSkus(_api::VirtualMachineImagesApi, location::String, publisherName::String, offer::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Vector{VirtualMachineImageResource}, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/publishers/{publisherName}/artifacttypes/vmimage/offers/{offer}/skus", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "publisherName", publisherName)  # type String
    Swagger.set_param(_ctx.path, "offer", offer)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export virtualMachineImagesGet, virtualMachineImagesList, virtualMachineImagesListOffers, virtualMachineImagesListPublishers, virtualMachineImagesListSkus


text:
using ZXCalculus: Scalar, add_power!, add_phase!

s = Scalar()
@test s * s == s
@test add_power!(s, 1) == Scalar(1, 0)
@test add_phase!(s, 1) == Scalar(1, 1)


text:
module NeuralNetwork
export NeuralNet, fit!, predict, score, weave!, addLayer!,
       addSigmoidLayer!, FullyConnected, backpropagation,
       feedforward, minibatch, evaluate, addSoftmaxLayer!

using ..Classification: Optimizer, GDContext, SGDContext, NAGDContext,
                        GradientDescent, StochasticGradientDescent,
                        NesterovGradientDescent, partition, postprocess!
import ..Classification: fit!, predict, score, logcrossentropy
using Random: randn
using Statistics: mean
using NPZ: npzwrite, npzread

abstract type Layer end

struct FullyConnected <: Layer
    inputsize::Int
    numneurons::Int
    W::Matrix{Float32}
    b::Vector{Float32}
    f::Function
    f′::Function
    function FullyConnected(inputsize, numneurons, transform::Function, derivative::Function)
        # Normal distributed random numbers with N(0, 0.1)
        weights = 1randn(numneurons, inputsize)
        bias = 1randn(numneurons)
        new(inputsize, numneurons, weights, bias, transform, derivative)
    end
end

function feed(layer::FullyConnected, a)
    layer.W * a .+ layer.b
end

function newsize(layer::FullyConnected, newsize)::FullyConnected
    FullyConnected(newsize, layer.numneurons, layer.f, layer.f′)
end

struct NeuralNet
    layers::Vector{Layer}
    ∇cost::Function
    NeuralNet(layers::Vector{Layer}, ∇f::Function=MSE′) = new(layers, ∇f)
end
NeuralNet(layer::Layer) = NeuralNet([layer], MSE′)
NeuralNet() = NeuralNet(Layer[])

function NeuralNet(neurons::Vector{Int})
    net = NeuralNet()
    if length(neurons) > 0
        # First layer keeps the input 0 until an input is given
        push!(net.layers, FullyConnected(0, neurons[1], σ, ∇σ))
        for numneurons in neurons[2:end]
            addSigmoidLayer!(net, numneurons)
        end
    end
    net
end

function NeuralNet(architecture::AbstractArray{Tuple{Int, Symbol}})
    layers = Layer[]
    numneurons, activation = architecture[1]
    a, ∇a = getactivation(activation)
    push!(layers, FullyConnected(0, numneurons, a, ∇a))
    for i in 2:length(architecture)
        input, _ = architecture[i-1]
        numneurons, activation = architecture[i]
        a, ∇a = getactivation(activation)
        push!(layers, FullyConnected(input, numneurons, a, ∇a))
    end
    NeuralNet(layers)
end

function getactivation(f::Union{Symbol, Function})::Tuple{Function, Function}
    if f == σ || f == :sigmoid
        σ, ∇σ
    elseif f == relu || f == :relu
        relu, ∇relu
    elseif f == leakyrelu || f == :leakyrelu
        leakyrelu, ∇leakyrelu
    elseif f == :id
        x -> x, x -> one(x)
    else
        throw(KeyError("Activation $f not supported"))
    end
end

function addLayer!(net::NeuralNet, numneurons, transform::Function, derivative::Function)
    inputsize = try
        net.layers[end].numneurons
    catch
        0
    end
    layer = FullyConnected(inputsize, numneurons, transform, derivative)
    push!(net.layers, layer)
end

function addSigmoidLayer!(net::NeuralNet, numneurons)
    addLayer!(net, numneurons, σ, ∇σ)
end

function addReluLayer!(net::NeuralNet, numneurons)
    addLayer!(net, numneurons, relu, ∇relu)
end

function addLeakyReluLayer!(net::NeuralNet, numneurons)
    addLayer!(net, numneurons, leakyrelu, ∇leakyrelu)
end

function weave!(net, X, y)
    # X has dimension #predictors×#samples
    if length(net.layers) == 0
        push!(net.layers, FullyConnected(size(X, 2), 1, σ, ∇σ))
    else
        net.layers[1] = newsize(net.layers[1], size(X, 1))
    end
end

function backpropagation(net::NeuralNet, X, y)
    ∇b = [zeros(Float32, size(layer.b)) for layer in net.layers]
    ∇W = [zeros(Float32, size(layer.W)) for layer in net.layers]

    # Forward propagation
    a = X
    as = [a]
    zs = []
    for layer in net.layers
        z = layer.W * a .+ layer.b
        a = layer.f.(z)
        push!(zs, z)
        push!(as, a)
    end

    # Backpropagation
    # Handle final layer L
    ∇aᴸC = net.∇cost(as[end] |> softmax, y)
    δ = ∇aᴸC .* net.layers[end].f′.(zs[end])
    ∇b[end] .= δ
    ∇W[end] .= δ*as[end-1]'

    # Handle the remaining layers
    for l in 1:(length(net.layers)-1)
        ∂a∂z = net.layers[end-l].f′.(zs[end-l])
        δ = net.layers[end-l+1].W' * δ .* ∂a∂z
        ∇b[end-l] .= δ
        ∇W[end-l] .= δ * as[end-l-1]'
    end
    ∇b, ∇W
end

function feedforward(net::NeuralNet, X)
    for layer in net.layers
        X = layer.f.(layer.W*X .+ layer.b)
    end
    X .= softmax(X)
    X
end

function evaluate(net::NeuralNet, X, labels)::Float32
    ŷ = predict(net, X)
    @assert size(ŷ) == size(labels)
    mean(labels .== ŷ)
end

function evaluate!(net::NeuralNet, o::Optimizer, X, labels)::Nothing
    ŷ = predict(net, X)
    @assert size(ŷ) == size(labels)
    o.loss[o.iterations] = mean(labels .== ŷ)

    !o.hasvalidationset && return

    ŷ = predict(net, o.validationset[1])
    o.validationloss[o.iterations] = mean(o.validationset[2]' .== ŷ)
    return
end

function logcrossentropy(net::NeuralNet, X, y)
    ŷ = feedforward(net, X)
    loss = logcrossentropy(y, ŷ)
end

function predict(net::NeuralNet, X)::Matrix{Int}
    probabilities = feedforward(net, X)
    return [imax[1] - 1 for imax in argmax(probabilities, dims=1)]
end

function MSE(net::NeuralNet, X, y)
    ŷ = feedforward(net, X)
    mean((y .- ŷ).^2)
end

function MSE′(ŷ, y)
    ŷ .- y
end

"""
    σ(x)

Sigmoid activation function
"""
σ(x::Real) = one(x)/(one(x) + exp(-x))
∇σ(x::Real) = σ(x)*(one(x) - σ(x))

"""
    relu(x)

Rectified Linear Unit activation function
"""
relu(x::Real) = max(zero(x), x)
∇relu(x::Real) = if x ≥ zero(x) one(x) else zero(x) end

"""
    leakyrelu(x)

Leaky Rectified Linear Unit activation function
"""
leakyrelu(x::Real, a = oftype(x, 0.01)) = max(a*x, x)
∇leakyrelu(x::Real, a = oftype(x, 0.01)) = if x ≥ zero(x) one(x) else a end

function softmax(xs::AbstractArray; dims=1)
    max_ = maximum(xs, dims=dims)
    exp_ = exp.(xs .- max_)
    exp_ ./ sum(exp_, dims=dims)
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 1a318e7a-dd72-11eb-3200-3f78e8727c26
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
		Pkg.PackageSpec(name="ImageIO", version="0.5"),
		Pkg.PackageSpec(name="ImageShow", version="0.2"),
		Pkg.PackageSpec(name="FileIO", version="1.6"),
		Pkg.PackageSpec(name="PNGFiles", version="0.3.6"),
		Pkg.PackageSpec(name="Colors", version="0.12"),
		Pkg.PackageSpec(name="ColorVectorSpace", version="0.8"),

		Pkg.PackageSpec(name="PlutoUI", version="0.7"),
		Pkg.PackageSpec(name="HypertextLiteral", version="0.5"),
		Pkg.PackageSpec(name="ForwardDiff", version="0.10")
	])

	using Colors, ColorVectorSpace, ImageShow, FileIO
	using PlutoUI
	using HypertextLiteral
	using LinearAlgebra
	using ForwardDiff
end

# ╔═╡ 2cab962b-9adb-49e5-aa74-1c4737333595
f₁(x) = x^2;

# ╔═╡ 15554085-48c9-4d53-830a-697024b74abf
f₁(2)

# ╔═╡ a79bb8b4-bafa-401d-a21e-16cd448d1f6c
f₂(x) = sin(x);

# ╔═╡ 38ee8a54-4417-4df1-bf5a-49b8bd806250
f₂(2*pi)

# ╔═╡ ef00e8b4-c9e6-45f5-a3fc-e4ae3294109a
f₃(x) = 2x;

# ╔═╡ b43b2394-a63d-44d2-ba19-a0952f2081b9
f₃(3)

# ╔═╡ c23b8836-a50d-4fcf-bf4f-0b036dfa103a
α = 5;

# ╔═╡ eefd0a83-9427-4cee-8aee-d4ba931051a6
ForwardDiff.derivative(f₁, 4)

# ╔═╡ Cell order:
# ╠═1a318e7a-dd72-11eb-3200-3f78e8727c26
# ╠═2cab962b-9adb-49e5-aa74-1c4737333595
# ╠═15554085-48c9-4d53-830a-697024b74abf
# ╠═a79bb8b4-bafa-401d-a21e-16cd448d1f6c
# ╠═38ee8a54-4417-4df1-bf5a-49b8bd806250
# ╠═ef00e8b4-c9e6-45f5-a3fc-e4ae3294109a
# ╠═b43b2394-a63d-44d2-ba19-a0952f2081b9
# ╠═c23b8836-a50d-4fcf-bf4f-0b036dfa103a
# ╠═eefd0a83-9427-4cee-8aee-d4ba931051a6


text:
##########################################################
#
#  Create DifferentiableFunction
#
##########################################################

abstract type DifferentiableFunction end

################ creates L2Loss =>  f(x) = ||x - y||^2 / 2

# N = 1 --> y
# N = 0 --> y = 0
struct L2Loss{T<:AbstractFloat, N, M<:AbstractArray} <: DifferentiableFunction
  y::M
end

L2Loss(y::AbstractArray{T}) where {T<:AbstractFloat} = L2Loss{T, 1, typeof(y)}(y)
L2Loss() = L2Loss{Float64, 0, Vector{Float64}}(Float64[])

function value(f::L2Loss{T, N}, x::AbstractArray{T}) where {T<:AbstractFloat, N}
  s = zero(T)
  if N == 1
    y = f.y
    @assert size(y) == size(x)
    @inbounds @simd for i in eachindex(x)
      s += (x[i]-y[i])^2.
    end
  else
    s = sum(abs2, x)
  end
  s / 2.
end

function value_and_gradient!(f::L2Loss{T, N}, hat_x::AbstractArray{T}, x::AbstractArray{T}) where {T<:AbstractFloat, N}
  if N == 1
    @. hat_x = x - f.y
  else
    copyto!(hat_x, x)
  end
  sum(abs2, hat_x) / 2.
end

################ creates a function x'Ax/2 + b'x + c

struct QuadraticFunction{T<:AbstractFloat, N, M<:AbstractMatrix, V} <: DifferentiableFunction
  A::M
  b::V
  c::T
  tmp::Vector{T}    ## call to value does not allocate
end

QuadraticFunction(A::AbstractMatrix{T}) where {T<:AbstractFloat} = QuadraticFunction{T, 1, typeof(A), Vector{T}}(A, T[], zero(T), Array{T}(undef, size(A, 1)))
QuadraticFunction(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T<:AbstractFloat} = QuadraticFunction{T, 2, typeof(A), typeof(b)}(A, b, zero(T), Array{T}(undef, size(A, 1)))
QuadraticFunction(A::AbstractMatrix{T}, b::AbstractVector{T}, c::T) where {T<:AbstractFloat} = QuadraticFunction{T, 3, typeof(A), typeof(b)}(A, b, c, Array{T}(undef, size(A, 1)))

function value(f::QuadraticFunction{T, 1}, x::StridedVector{T}) where {T<:AbstractFloat}
  mul!(f.tmp, f.A, x)
  dot(x, f.tmp) / 2.
end
function value(f::QuadraticFunction{T, 2}, x::StridedVector{T}) where {T<:AbstractFloat}
  mul!(f.tmp, f.A, x)
  dot(x, f.tmp) / 2. + dot(x, f.b)
end
function value(f::QuadraticFunction{T, 3}, x::StridedVector{T}) where {T<:AbstractFloat}
  mul!(f.tmp, f.A, x)
  dot(x, f.tmp) / 2. + dot(x, f.b) + f.c
end

function gradient!(f::QuadraticFunction{T, N}, hat_x::StridedVector{T}, x::StridedVector{T}) where {T<:AbstractFloat, N}
  mul!(hat_x, f.A, x)
  if N > 1
    @. hat_x += f.b
  end
  hat_x
end

function value_and_gradient!(f::QuadraticFunction{T, N}, hat_x::StridedVector{T}, x::StridedVector{T}) where {T<:AbstractFloat, N}
  b = f.b
  mul!(hat_x, f.A, x)
  r = dot(hat_x, x) / 2.
  if N > 1
    @. hat_x += b
    r += dot(x, b)
  end
  r + f.c
end

################ creates a function |Y - X⋅β|_2^2 / 2.

struct LeastSquaresLoss{T, Ty<:AbstractVecOrMat, Tx<:AbstractMatrix} <: DifferentiableFunction
  Y::Ty
  X::Tx
  tmp::VecOrMat{T}    ## call to value does not allocate
end

LeastSquaresLoss(Y::AbstractVecOrMat{T}, X::AbstractMatrix{T}) where {T<:AbstractFloat} =
    LeastSquaresLoss{T, typeof(Y), typeof(X)}(Y, X, zeros(T, size(Y)))

function value(f::LeastSquaresLoss{T}, x) where {T<:AbstractFloat}
  mul!(f.tmp, f.X, x)
  v = zero(T)
  @inbounds for i in eachindex(f.tmp)
    v += (f.Y[i] - f.tmp[i])^2.
  end
  v / 2.
end

function value_and_gradient!(
  f::LeastSquaresLoss{T},
  grad_out::StridedVector{T},
  x) where {T<:AbstractFloat}

  mul!(f.tmp, f.X, x)
  @. f.tmp -= f.Y
  mul!(grad_out, transpose(f.X), f.tmp)
  sum(abs2, f.tmp) / 2.
end


text:
################################################################################
#
#  GrpAb/SubgroupEnum.jl : Subgroup enumeration for finitely generated
#                          abelian groups.
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2017 Tommy Hofmann, Claus Fieker
#
################################################################################

export minimal_subgroups, psubgroups, index_p_subgroups, subgroups

################################################################################
#
#  Enumeration of subgroups of index p
#
################################################################################

function index_p_subgroups(A::GrpAbFinGen, p::Integer)
  return index_p_subgroups(A, fmpz(p))
end

mutable struct IndexPSubgroups{S, T}
  p::Int
  n::UInt
  st::Int
  mp::S
  c::fmpz_mat
  mthd::T

  function IndexPSubgroups{T}(A::GrpAbFinGen, p::fmpz, mthd::T = sub) where {T}
    if order(A) % p != 0
      r = new{Generic.IdentityMap{GrpAbFinGen}, T}()
      r.n = 0
      return r
    end
    s, ms = snf(A)  # ms: A -> s
    r = new{typeof(inv(ms)), T}()
    @assert s.issnf
    r.p = Int(p)
    r.mp = inv(ms)
    i=1
    while s.snf[i] % p != 0
      i += 1
    end
    r.st = i
    r.n = UInt(div(fmpz(p)^(length(s.snf)-i+1) - 1, fmpz(p)-1))
    r.c = zero_matrix(FlintZZ, length(s.snf), length(s.snf))
    r.mthd = mthd
    r.c
    return r
  end
end

function index_p_subgroups(A::GrpAbFinGen, p::Union{fmpz, Integer}, mthd::T = sub) where {T}
  q = fmpz(p)
  @assert isprime(q)
  I = IndexPSubgroups{T}(A, q, mthd)

  #a subgroup of index p corresponds to a HNF with exactly one p on the
  #diagonal - and the other entries arbitrary reduced.
  #so it should be 1 + p + p^2 + + p^(n-1) = ((p^n)-1)/(p-1) many
  return I
end

function index_to_group(s::IndexPSubgroups, i::UInt)
  j = 1
  x = 1
  while i>=x
    i -= x
    x *= s.p
    j += 1
  end
  c = s.c
  zero!(c)
  for k=1:nrows(c)
    if s.st+j-1 != k
      c[k, k] = 1
    end
  end
  k = 1
  while i != 0
    c[k+s.st-1, s.st + j-1] = i%s.p
    i = div(i, s.p)
    k += 1
  end
  c[s.st + j-1, s.st + j-1] = s.p
  gen = [s.mp\(GrpAbFinGenElem(codomain(s.mp), sub(c, l:l, 1:ncols(c)))) for l=1:nrows(c)]
  return s.mthd(domain(s.mp), gen)
end

function Base.iterate(s::IndexPSubgroups, i::UInt = UInt(0))
  if i + 1 > s.n
    return nothing
  end

  return index_to_group(s, i), (i + 1)
end

#function Base.iterate(s::IndexPSubgroups, i::UInt)
#  if i > s.n - 1
#    return nothing
#  end
#
#  return index_to_group(s, i), i + 1
#end

Base.IteratorSize(::Type{IndexPSubgroups{S, T}}) where {S, T} = Base.HasLength()

function Base.length(s::IndexPSubgroups)
  return s.n
end

#function Base.done(s::IndexPSubgroups, i::UInt)
#  return i>=s.n
#end

#=
example:
 julia> sg = index_p_subgroups(GrpAbFinGen([3,3,3,3]), 3)
 julia> index_to_group(sg, UInt(6));
=#

################################################################################
#
#  Subgroup iterators for p-groups
#
################################################################################

# We use the algorithm of Butler, Subgroup Lattices and Symmetric Functions,
# Section 1.6. But we use the description of Cohen, Advanced Topics in
# Computational Number Theory, Theorem 4.1.18.

# Some notation/convention:
# An abelian p-group A has a unique type (n_1,...,n_r), which is decreasing
# with the property that A \cong Z/p^{n_1} x ... x Z/p^{n_r}.

# Given a p-group with type x, the following function enumerates all possible
# types y of a subgroup of length t. The only restrictions on y are
# length(y) <= length(x), y_1 \leq x_1 and y_i \leq min(y_{i-1}, x_i) for
# 2 \leq i \leq length(y). This is part (1) of Theorem 4.1.18.

struct yIterator
  t::Int
  x::Vector{Int}
  nulls::Vector{Int}
  res::Vector{Int}

  function yIterator(x::Vector{Int}, t::Int)
    z = new(t, x, zeros(Int, length(x) - t), zeros(Int, length(x)))
    return z
  end
end

function Base.iterate(F::yIterator, i::Vector{Int})
  if length(i) == 0
    return nothing
  end

  done = true

  for j in 1:length(i)
    if i[j] != F.x[j]
      done = false
    end
  end

  if done
    return nothing
  end

  if i[1] < F.x[1]
    i[1] = i[1] + 1
  else # the first one is as large as possible
    j = 0
    for j in 2:length(i)
      if i[j] < F.x[j]
        i[j] = i[j] + 1
        for k in j-1:-1:1
          i[k] = i[j]
        end
        break
      end
    end
  end

  for j in 1:F.t
    F.res[j] = i[j]
  end

  return copy(F.res), i
end

function Base.iterate(F::yIterator)
  i = ones(Int, F.t)

  if F.t == 0
    return copy(F.res), i
  end

  @inbounds if F.t > 0
    i[1] = 0
  end

  @inbounds if i[1] < F.x[1]
    i[1] = i[1] + 1
  else # the first one is as large as possible
    j = 0
    for j in 2:length(i)
      if i[j] < F.x[j]
        i[j] = i[j] + 1
        for k in j-1:-1:1
          i[k] = i[j]
        end
        break
      end
    end
  end

  @inbounds for j in 1:F.t
    F.res[j] = i[j]
  end

  return copy(F.res), i
end

Base.IteratorSize(::Type{yIterator}) = Base.SizeUnknown()

Base.eltype(::Type{yIterator}) = Vector{Int}

_subpartitions(x) = Iterators.flatten((yIterator(x, t) for t in 0:length(x)))

# Given a type y for the subgroup, we can iterator through all possible
# valid permutations. This is part (2) of Theorem 4.1.18.

struct SigmaIteratorGivenY{T}
  gen::T
end

@inline function _isvalid(s, t, x, y, sigma)
  for i in 1:t
    if y[i] > x[sigma[i]]
      return false
    end
  end
  for i in 1:(s-1)
    # TODO: Precompute the indicies where they are equal
    if y[i] == y[i + 1]
      if sigma[i] < sigma[i + 1]
        return false
      end
    end
  end
  return true
end

#(::Colon)(x::Int, y::Nothing) = 1:0

Base.:(:)(x::Int, y::Nothing) = 1:0

Base.:(:)(x::Int, y::fmpz) = fmpz(x):y

function SigmaIteratorGivenY(s, x, y)
  t = something(findlast(!iszero, y), 0)
  SigmaIteratorGivenY(Iterators.filter(sigma -> _isvalid(s, t, x, y, sigma),
                                       (deepcopy(z.d) for z in AllPerms(s))))
end

Base.iterate(S::SigmaIteratorGivenY) = Base.iterate(S.gen)

Base.iterate(S::SigmaIteratorGivenY, s) = Base.iterate(S.gen, s)

Base.length(S::SigmaIteratorGivenY) = Base.length(S.gen)

Base.eltype(::Type{SigmaIteratorGivenY{T}}) where {T} = Vector{Int}

# for some reason this is type unstable.

# Finally we have to enumerate the possible matrices defining the generators
# of the subgroup (step (3)). There are two parts to this. First we compute,
# for a given sigma, the indice in the matrix.

function compute_indice(s, t, x, y, sigma)
  tau = Nemo.inv!(perm(sigma))
  indice = Tuple{Int, Int, Int}[]
  @inbounds for j in 1:t
    for i in 1:s
      if tau[i] > j
        if i < sigma[j] # case a)
          push!(indice, (i, j, 0))
        elseif i > sigma[j] && x[i] < y[j] #case b)
          push!(indice, (i, j, 1))
        elseif i > sigma[j] && x[i] >= y[j] # case c)
          push!(indice, (i, j, 2))
        end
      end
    end
  end
  return [( indice, tau )]
end

# Given sigma, we can iteraterate through all possible choices of intergers
# c_{i,j}. To do this we collect the single ranges at the indice and the
# produce a product iterator.

mutable struct cIteratorGivenSigma{T}
  s::Int
  t::Int
  x::Vector{Int}
  y::Vector{Int}
  p::Int
  sigma::Vector{Int}
  tau::Vector{Int}
  indice::Vector{Tuple{Int, Int, Int}}
  it::T
end

function _cIteratorGivenSigma(s::Int, t::Int, x::Vector{Int},
                              y::Vector{Int}, p::Union{fmpz, Integer}, sigma::Vector{Int})
  pp = Int(p)
  tau = Nemo.inv!(perm(sigma))
  indice, it = getintervals(t, s, x, y, pp, sigma, tau)
  return cIteratorGivenSigma{typeof(it)}(s, t, x, y, pp, sigma, tau, indice, it)
end

function getintervals(t, s, x, y, p, sigma, tau)
  indice = Tuple{Int, Int, Int}[]
  ranges = UnitRange{Int}[]

  @inbounds for j in 1:t
    for i in 1:s
      if tau[i] > j
        if i < sigma[j] # case a)
          push!(indice, (i, j, 0))
          push!(ranges, 1:p^(y[j] - y[tau[i]]))
        elseif i > sigma[j] && x[i] < y[j] #case b)
          push!(indice, (i, j, 1))
          push!(ranges, 1:p^(x[i] - y[tau[i]]))
        elseif i > sigma[j] && x[i] >= y[j] # case c)
          push!(indice, (i, j, 2))
          push!(ranges, 1:p^(y[j] - y[tau[i]] - 1))
        end
      end
    end
  end
  if length(indice) == 0
    return indice, (Int[] for z in 1:1)
  end
  return indice, Base.product(ranges...)
end

Base.iterate(C::cIteratorGivenSigma) = Base.iterate(C.it)

Base.iterate(C::cIteratorGivenSigma, s) = Base.iterate(C.it, s)

Base.length(C::cIteratorGivenSigma) = Base.length(C.it)

Base.eltype(::Type{cIteratorGivenSigma{T}}) where {T} = Base.eltype(T)

# Now we have the y, the permutation sigma with inverse tau, the indices ind
# and the corresponding ranges describing the exponents in the generator
# matrix. Now get the matrix!

function get_matrix(s, t, ind, c, sigma, tau, p, x, y)
  M = zeros(Int, s, t)
  get_matrix!(M, s, t, ind, c, sigma, tau, p, x, y)
  return M
end

function get_matrix!(M, s, t, ind, c, sigma, tau, p, x, y)
  @inbounds for j in 1:t
    for i in 1:s
      if tau[i] < j
        M[i, j] = p^x[i]
      elseif tau[i] == j
        M[i,j] = p^(x[i] - y[j])
      end
    end
  end
  @inbounds for k in 1:length(ind)
    i, j, _case = ind[k]
    if _case == 0
      # case a)
      M[i, j] = c[k]*p^(x[i] - y[j])
    elseif _case == 1
      # case b)
      M[i, j] = c[k]
    elseif _case == 2
      # case c)
      M[i, j] = c[k] * p^(x[i] - y[j] + 1)
    end
  end
  return M
end

# Put everything together:
# Given the type x of a p-group G and a type y, this function iterates
# through all generating matrices of subgroups of type y.
# (The matrix is with respect to the canonical form
# G = Z/p^x[1] x .... x Z/p^x[r]

function _subgroup_type_iterator(x, y, p)
  @assert length(x) == length(y)
  s = length(x)
  t = something(findlast(!iszero, y), 0)

  # have to treat the empty y separately

  if any(y[i] > x[i] for i in 1:length(x))
    return (x for x in 1:-1)
  end

  if t == 0
    return (x for x in [zeros(Int, s, 0)])
  elseif x == y
    return (x for x in [Matrix{Int}(I, s, s)])
  end

  return (get_matrix(s, t, f[1], c, sigma, f[2], p, x, y)
          for sigma in SigmaIteratorGivenY(s, x, y)
          for f in compute_indice(s, t, x, y, sigma)
          for c in _cIteratorGivenSigma(s, t, x, y, p, sigma))
end

# Same as above but for all subgroups.
function _subgroup_iterator(x, p, types = (y for t in 0:length(x)
                                             for y in yIterator(x, t)))
  s = length(x)
  # Flatten just concatenates two iterators.
  #return #Iterators.flatten(([Array{Int}(s, 0)],
  return Iterators.flatten((_subgroup_type_iterator(x, y, p) for y in types))
end

# Given a matrix M and a group G, this function constructs elements from
# the columns of M. The indice allows to handle the case, where the
# generators of G correspond to a permutation of the rows of M.
function _matrix_to_elements(G::GrpAbFinGen, M::Matrix{Int},
                             indice::Vector{Int} = collect(1:ngens(G)))
  numgenssub = size(M, 2)
  numgen = ngens(G)
  r = size(M, 1)
  z = Array{GrpAbFinGenElem}(undef, numgenssub)
  v = zeros(Int, numgen)
  for i in 1:numgenssub
    for j in 1:r
      v[indice[j]] = M[j, i]
    end
    z[i] = G(v)
    for j in 1:numgen
      v[j] = 0
    end
  end
  return z
end

# Given a finitely generated p-group G in Smith normal form, and a type t,
# this function returns an iterator, which iterates over generators of
# subgroups of type t. If t = [-1], then there is no restriction on the type.
function __psubgroups_gens(G::GrpAbFinGen, p::Union{fmpz, Integer},
                           order, index, t::Vector{Int})
  @assert isfinite(G)
  @assert issnf(G)
  # The SNF can contain 1's and 0's
  # Have to help inference here with the explicit type
  x = Tuple{Int, Int}[ (valuation(G.snf[i], p), i)
                       for i in 1:length(G.snf) if G.snf[i] > 1]
  reverse!(x)
  Gtype = [ t[1] for t in x ]
  indice = [ t[2] for t in x ]
  if length(t) > length(Gtype)
    return ( i for i in 1:-1)
  end
  # x is the "type" of the p-group G as a partition
  if length(t) == 1 && t[1] == -1
    return (_matrix_to_elements(G, M, indice) for
            M in _subgroup_iterator(Gtype, p))
  else
    t = vcat(t, zeros(Int, length(x) - length(t)))
    return  (_matrix_to_elements(G, M, indice) for
             M in _subgroup_type_iterator(Gtype, t, p))
  end
  return Gtype, indice
end

function __psubgroups_gens(G::GrpAbFinGen, p::Union{fmpz, Integer},
                           order, index, types)
  @assert isfinite(G)
  @assert issnf(G)
  # The SNF can contain 1's and 0's
  # Have to help inference here with the explicit type
  x = Tuple{Int, Int}[ (valuation(G.snf[i], p), i)
                       for i in 1:length(G.snf) if G.snf[i] > 1]
  reverse!(x)
  Gtype = [ t[1] for t in x ]
  indice = [ t[2] for t in x ]
  # x is the "type" of the p-group G as a partition
  if order != -1
    v = valuation(order, p)
    adjusted_types = (vcat(t, zeros(Int, length(x) - length(t)))
                      for t in types if sum(t) == v)
  elseif index != -1
    v = sum(Gtype) - valuation(index, p)
    adjusted_types = (vcat(t, zeros(Int, length(x) - length(t)))
                      for t in types if sum(t) == v)
  else
    adjusted_types = (vcat(t, zeros(Int, length(x) - length(t))) for t in types)
  end

  return  (_matrix_to_elements(G, M, indice)
           for M in _subgroup_iterator(Gtype, p, adjusted_types))
  return Gtype, indice
end

function __psubgroups_gens(G::GrpAbFinGen, p::Union{Integer, fmpz}, order, index)
  @assert isfinite(G)
  @assert issnf(G)
  # The SNF can contain 1's and 0's
  # Have to help inference here with the explicit type
  x = Tuple{Int, Int}[ (valuation(G.snf[i], p), i)
                       for i in 1:length(G.snf) if G.snf[i] > 1]
  reverse!(x)
  Gtype = [ t[1] for t in x ]
  indice = [ t[2] for t in x ]
  types = _subpartitions(Gtype)
  # x is the "type" of the p-group G as a partition
  if order != -1
    v = valuation(order, p)
    adjusted_types = (vcat(t, zeros(Int, length(x) - length(t)))
                      for t in types if sum(t) == v)
  elseif index != -1
    v = sum(Gtype) - valuation(index, p)
    adjusted_types = (vcat(t, zeros(Int, length(x) - length(t)))
                      for t in types if sum(t) == v)
  else
    adjusted_types = (vcat(t, zeros(Int, length(x) - length(t))) for t in types)
  end
  return  (_matrix_to_elements(G, M, indice)
           for M in _subgroup_iterator(Gtype, p, adjusted_types))
end

# Same as above but now for arbitrary p-groups
function _psubgroups_gens(G::GrpAbFinGen, p, t, order, index)
  if issnf(G)
    if t == [-1]
      return __psubgroups_gens(G, p, order, index)
    else
      return __psubgroups_gens(G, p, order, index, t)
    end
  else
    S, mS = snf(G)
    if t == [-1]
      return ( map(x -> image(mS, x)::GrpAbFinGenElem, z)
               for z in __psubgroups_gens(S, p, order, index))
    else
      return ( map(x -> image(mS, x)::GrpAbFinGenElem, z)
               for z in __psubgroups_gens(S, p, order, index, t))
    end
  end
end

function _psubgroups_gens_quotype(G::GrpAbFinGen, p, t, order, index)
  if issnf(G)
    x = Tuple{Int, Int}[ (valuation(G.snf[i], p), i)
                       for i in 1:length(G.snf) if G.snf[i] > 1]
    reverse!(x)
    Gtype = [ t[1] for t in x ]
    indice = [ t[2] for t in x ]
    types = _subpartitions(Gtype)
    filtered_types = Iterators.filter(z -> sum(z) == sum(Gtype) - sum(t), types)
    return __psubgroups_gens(G, p, order, index, filtered_types)
  else
    S, mS = snf(G)
    return ( map(x -> image(mS, x)::GrpAbFinGenElem, z)
             for z in _psubgroups_gens_quotype(S, p, t, order, index))
  end
end

function _ptype(G, p)
  Gsnf = domain(snf(G)[2])
  x = Int[ valuation(Gsnf.snf[i], p)
                       for i in 1:length(Gsnf.snf) if Gsnf.snf[i] > 1]
  reverse!(x)
  t = findlast(!iszero, x)
  return x[1:t]
end

# Same as above but now allow a function to be applied to the output
function _psubgroups(G::GrpAbFinGen, p::Union{Integer, fmpz}; subtype = [-1],
                                                              quotype = [-1],
                                                              order = -1,
                                                              index = -1,
                                                              fun = sub)
  P, mP = psylow_subgroup(G, p, false)

  if quotype != [-1]
    return ( fun(G, map(mP, z))
            for z in _psubgroups_gens_quotype(P, p, quotype, order, index)
            if _ptype(quo(G, map(mP, z), false)[1], p) == quotype)
  end

  return ( fun(G, map(mP, z)) for z in _psubgroups_gens(P, p, subtype, order, index))
end

# We use a custom type for the iterator to have pretty printing.
mutable struct pSubgroupIterator{F, T, E}
  G::GrpAbFinGen
  p::fmpz
  subtype::Vector{Int}
  quotype::Vector{Int}
  index::fmpz
  order::fmpz
  fun::F
  it::T
end

function Base.show(io::IO, I::pSubgroupIterator)
  print(io, "p-subgroup iterator of \n$(I.G)\n")

  print(io, "p: $(I.p)\n")

  if I.subtype == [-1]
    print(io, "subgroup type: any\n")
  else
    print(io, "subgroup type: $(I.subtype)\n")
  end

  if I.quotype == [-1]
    print(io, "quotient type: any\n")
  else
    print(io, "quotient type: $(I.quotype)\n")
  end

  if I.index == -1
    print(io, "index: any\n")
  else
    print(io, "index: $(I.index)\n")
  end

  if I.order == -1
    print(io, "order: any\n")
  else
    print(io, "order: $(I.order)\n")
  end

  if I.fun != sub
    print(io, "function: $(I.fun)")
  end
end

function pSubgroupIterator(G::GrpAbFinGen, p::Union{fmpz, Integer};
                                           subtype::Vector{Int} = [-1],
                                           quotype::Vector{Int} = [-1],
                                           index::Union{fmpz, Int} = -1,
                                           order::Union{fmpz, Int} = -1,
                                           fun = sub)
  if index == p
    it = index_p_subgroups(G, p, fun)
  else
    it = _psubgroups(G, p; subtype = subtype, quotype = quotype,
                           fun = fun, index = index, order = order)
  end

  E = Core.Compiler.return_type(fun, (GrpAbFinGen, Vector{GrpAbFinGenElem}))

  z = pSubgroupIterator{typeof(fun), typeof(it), E}(G, fmpz(p), subtype, [-1],
                                                    fmpz(index), fmpz(order), fun, it)
  return z
end

@doc Markdown.doc"""
    psubgroups(g::GrpAbFinGen, p::Integer;
               subtype = :all,
               quotype = :all,
               index = -1,
               order = -1)

Return an iterator for the subgroups of $G$ of the specific form. Note that
`subtype` (and `quotype`) is the type of the subgroup as an abelian $p$-group.
"""
function psubgroups(G::GrpAbFinGen, p::Union{Integer, fmpz}; subtype = :all,
                                                             quotype = :all,
                                                             index =  -1,
                                                             order = -1,
                                                             fun = sub)

  options = Int16[ subtype != :all, quotype != :all, index != -1, order != -1 ]

  if sum(options) > 1
    error("Currently only one non-default parameter is supported.")
  end

  if subtype == :all
    _subtype = [-1]
  else
    for i in 1:(length(subtype) - 1)
      if subtype[i + 1] > subtype[i]
        error("Subtype must be a partition")
      end
    end
    _subtype = subtype
  end

  if quotype == :all
    _quotype = [-1]
  else
    for i in 1:(length(quotype) - 1)
      if quotype[i + 1] > quotype[i]
        error("Subtype must be a partition")
      end
    end
    _quotype = quotype
  end

  return pSubgroupIterator(G, p; subtype = _subtype, quotype = _quotype, order = order, index = index,
                                 fun = fun)
end

Base.iterate(S::pSubgroupIterator) = Base.iterate(S.it)

Base.iterate(S::pSubgroupIterator, st) = Base.iterate(S.it, st)

Base.eltype(::Type{pSubgroupIterator{F, T, E}}) where {F, T, E} = E

Base.IteratorSize(::Type{pSubgroupIterator{F, T, E}}) where {F, T, E} = Base.SizeUnknown()

################################################################################
#
#  Subgroup enumeration
#
################################################################################

mutable struct SubgroupIterator{F, T, E}
  G::GrpAbFinGen
  subtype::Vector{Int}
  quotype::Vector{Int}
  index::fmpz
  order::fmpz
  fun::F
  it::T
end

function Base.show(io::IO, I::SubgroupIterator)
  print(io, "p-subgroup iterator of \n$(I.G)\n")

  if I.subtype == [-1]
    print(io, "subgroup type: any\n")
  else
    print(io, "subgroup type: $(I.subtype)\n")
  end

  if I.quotype == [-1]
    print(io, "quotient type: any\n")
  else
    print(io, "quotient type: $(I.quotype)\n")
  end

  if I.index == -1
    print(io, "index: any\n")
  else
    print(io, "index: $(I.index)\n")
  end

  if I.order == -1
    print(io, "order: any\n")
  else
    print(io, "order: $(I.order)\n")
  end

  if I.fun != sub
    print(io, "function: $(I.fun)")
  end
end

Base.iterate(S::SubgroupIterator) = Base.iterate(S.it)

Base.iterate(S::SubgroupIterator, s) = Base.iterate(S.it, s)

Base.IteratorSize(::Type{SubgroupIterator{F, T, E}}) where {F, T, E} = Base.SizeUnknown()

Base.eltype(::Type{SubgroupIterator{F, T, E}}) where {F, T, E} = E

function _subgroups_gens(G::GrpAbFinGen, subtype::Vector{S} = [-1],
                         quotype = [-1], suborder = -1,
                         subindex = -1) where S <: Union{Integer, fmpz}
  primes = fmpz[]

  pgens = []

  if quotype != [-1]

    fac = factor(order(G))

    for (p, e) in fac
      if !iszero(mod(order(G), p))
        error("no subgroup exists")
      end
      ptype = map(l -> valuation(l, p), quotype)
      filter!( z -> z > 0, ptype)
      sort!(ptype, rev = true)
      T = psubgroups(G, Int(p), quotype = ptype, fun = (g, m) -> sub(g, m, false))
      genss = ( GrpAbFinGenElem[ t[2](x) for x in gens(t[1]) ] for t in T )
      push!(pgens, genss)
    end
  elseif subtype != [-1]
    for l in subtype
      fac = factor(fmpz(l))
      for (p, e) in fac
        push!(primes, p)
      end
    end

    primes = unique(primes)

    for p in primes
      if !iszero(mod(order(G), p))
        error("no subgroup exists")
      end
      ptype = map(l -> valuation(l, p), subtype)
      filter!( z -> z > 0, ptype)
      sort!(ptype, rev = true)
      T = psubgroups(G, Int(p), subtype = ptype, fun = (g, m) -> sub(g, m, false))
      genss = ( GrpAbFinGenElem[ t[2](x) for x in gens(t[1]) ] for t in T )
      push!(pgens, genss)
    end
  elseif suborder != -1 || subindex != -1
    # Do the order/index at the same time
    if subindex != -1
      _suborder = divexact(order(G), subindex)
    else
      _suborder = suborder
    end
    fac = factor(fmpz(_suborder))
    for (p, e) in fac
      orderatp = p^e
      T = psubgroups(G, Int(p), order = orderatp, fun = (g, m) -> sub(g, m, false))
      genss = ( GrpAbFinGenElem[ t[2](x) for x in gens(t[1]) ] for t in T )
      push!(pgens, genss)
    end
  else
    fac = factor(order(G))
    for (p, e) in fac
      T = psubgroups(G, Int(p), fun = (g, m) -> sub(g, m, false))
      genss = ( GrpAbFinGenElem[ t[2](x) for x in gens(t[1]) ] for t in T )
      push!(pgens, genss)
    end
  end

  final_it = ( vcat(c...) for c in Iterators.product(pgens...))
  return final_it
end

# Same as above but now allow a function to be applied to the output
function _subgroups(G::GrpAbFinGen; subtype = [-1], quotype = [-1], order = -1,
                                    index = -1, fun = sub)
  return ( fun(G, convert(Vector{GrpAbFinGenElem}, z)) for z in _subgroups_gens(G, subtype, quotype, order, index))
end

function SubgroupIterator(G::GrpAbFinGen; subtype::Vector{Int} = [-1],
                                          quotype::Vector{Int} = [-1],
                                          index::Union{fmpz, Int} = -1,
                                          order::Union{fmpz, Int} = -1,
                                          fun = sub)

  if index != -1 && isprime(fmpz(index))
    it = index_p_subgroups(G, fmpz(index), fun)
  else
    it = _subgroups(G; subtype = subtype, quotype = quotype,
                       fun = fun, index = index, order = order)
  end

  E = Core.Compiler.return_type(fun, (GrpAbFinGen, Vector{GrpAbFinGenElem}))

  z = SubgroupIterator{typeof(fun), typeof(it), E}(G, subtype, quotype,
                                                   fmpz(index), fmpz(order),
                                                   fun, it)
  return z
end

@doc Markdown.doc"""
    subgroups(g::GrpAbFinGen;
              subtype = :all ,
              quotype = :all,
              index = -1,
              order = -1)

Return an iterator for the subgroups of $G$ of the specific form.
"""
function subgroups(G::GrpAbFinGen; subtype = :all,
                                   quotype = :all,
                                   index =  -1,
                                   order = -1,
                                   fun = sub)

  # Handle the parameters

  options = Int16[ subtype != :all, quotype != :all, order != -1, index != -1]

  if mod(Hecke.order(G), index) != 0
    error("Index must divide the group order")
  end

  if mod(Hecke.order(G), order) != 0
    error("Index must divide the group order")
  end

  if sum(options) > 1
    error("Currently only one non-default parameter is supported.")
  end

  if subtype == :all
    _subtype = [-1]
  else
    _subtype = subtype
  end

  if quotype == :all
    _quotype = [-1]
  else
    _quotype = quotype
  end

  return SubgroupIterator(G; subtype = _subtype, quotype = _quotype, order = order, index = index,
                                 fun = fun)
end

################################################################################
#
#  Minimal subgroups
#
################################################################################

@doc doc"""
    minimal_subgroups(G::GrpAbFinGen) -> Vector{Tuple{GrpAbFinGen, Map}}

Return the minimal subgroups of $G$.
"""
function minimal_subgroups(G::GrpAbFinGen, add_to_lattice::Bool = false)
  @req isfinite(G) "Group must be finite"
  o = order(G)
  l = prime_divisors(o)
  res = Vector{Tuple{GrpAbFinGen, GrpAbFinGenMap}}()
  for p in l
    append!(res, psubgroups(G, p, order = p, fun = (x, m) -> sub(x, m, add_to_lattice)))
  end
  return res
end


text:
function generate_results()
"""
Generates the results for the paper
"""

    ##################################
    ##################################
    # Load base models and tools
    ##################################
    ##################################

    uncert = BSON.load("results/uncert_results.bson")
    learn = BSON.load("results/learning_results_adapt_20.bson")
    rc = BSON.load("results/rc_results_3.9.bson")
    rcl = BSON.load("results/rcl_results_3.9_adapt_20.bson")
    @everywhere include("tools.jl")

    ##################################
    ##################################
    # Figures 5a and 7
    ##################################
    ##################################

    println("Generating data for Figures 5a and 7.")

    # Pull out tax components
    uncert_tax = Dict{Symbol,Vector}()
    learn_tax = Dict{Symbol,Vector}()
    rc_tax = Dict{Symbol,Vector}()
    rcl_tax = Dict{Symbol,Vector}()
    for (key,value) in uncert["tax_channels"]
            uncert_tax[key] = vec(uncert["tax_channels"][key])
            learn_tax[key] = vec(learn["tax_channels"][key])
            rc_tax[key] = vec(rc["tax_channels"][key])
            rcl_tax[key] = vec(rcl["tax_channels"][key])
    end

    ### uncertainty
    uncert_df = convert(DataFrame, uncert_tax)
    uncert_df[:, end+1] = "Uncertainty"
    uncert_df[:, end+1] = 2005:2005+size(uncert_df, 1)-1
    uncert_t = convert(Matrix, uncert_df)

    ### learning
    learn_df = convert(DataFrame, learn_tax)
    learn_df[:, end+1] = "Learning"
    learn_df[:, end+1] = 2005:2005+size(learn_df, 1)-1
    learn_t = convert(Matrix, learn_df)

    ### robust control
    rc_df = convert(DataFrame, rc_tax)
    rc_df[:, end+1] = "Robust Control"
    rc_df[:, end+1] = 2005:2005+size(rc_df, 1)-1
    rc_t = convert(Matrix, rc_df)

    ### robust control + learning
    rcl_df = convert(DataFrame, rcl_tax)
    rcl_df[:, end+1] = "RC + Learning"
    rcl_df[:, end+1] = 2005:2005+size(rcl_df, 1)-1
    rcl_t = convert(Matrix, rcl_df)

    figure_5a_7 = vcat(uncert_t, learn_t, rc_t, rcl_t)
    figure_5a_7 = convert(DataFrame, figure_5a_7)
    names!(figure_5a_7, [:ce_tax, :exp_tax, :learning_insurance,
        :output_insurance, :precaution, :robust_insurance,
        :scc_direct, :smoothing, :tax_lower, :tax_upper, :unc_adjust,
        :type, :year])

    CSV.write("generate_plots/data/figure_5a_and_7.csv", figure_5a_7)

    ##################################
    ##################################
    # Figures 5b, 5c, and 5d
    ##################################
    ##################################

    println("Generating data for Figures 5a, 5b, and 5c.")

    ### Load in files for specific damage parameters
    for run in zip(["d1_0.01_d2_1.88", "d1_0.00556486_d2_3.0", "d1_0.01_d2_3.0"], ["5b", "5c", "5d"])

        learn_alt = BSON.load("results/learn_results_dam_" * run[1] * ".bson")
        uncert_alt = BSON.load("results/uncert_results_dam_" * run[1] * ".bson")
        rcl_alt = BSON.load("results/rcl_results_dam_" * run[1] * ".bson")
        rc_alt = BSON.load("results/rc_results_dam_" * run[1] * ".bson")

        uncert_tax = Dict{Symbol,Vector}()
        learn_tax = Dict{Symbol,Vector}()
        rc_tax = Dict{Symbol,Vector}()
        rcl_tax = Dict{Symbol,Vector}()

        for (key, value) in uncert_alt["tax_channels"]
            if (key != :robust_insurance2)
                uncert_tax[key] = vec(uncert_alt["tax_channels"][key])
                learn_tax[key] = vec(learn_alt["tax_channels"][key])
                rc_tax[key] = vec(rc_alt["tax_channels"][key])
                rcl_tax[key] = vec(rcl_alt["tax_channels"][key])
            end
        end

        ### uncertainty
        uncert_df = convert(DataFrame, uncert_tax)
        uncert_df[!, end+1] .= "Uncertainty"
        uncert_df[:, end+1] = 2005:2005+size(uncert_df, 1)-1
        uncert_t = convert(Matrix, uncert_df)

        ### learning
        learn_df = convert(DataFrame, learn_tax)
        learn_df[!, end+1] .= "Learning"
        learn_df[:, end+1] = 2005:2005+size(learn_df, 1)-1
        learn_t = convert(Matrix, learn_df)

        ### robust control
        rc_df = convert(DataFrame, rc_tax)
        rc_df[!, end+1] .= "Robust Control"
        rc_df[:, end+1] = 2005:2005+size(rc_df, 1)-1
        rc_t = convert(Matrix, rc_df)

        ### robust control + learning
        rcl_df = convert(DataFrame, rcl_tax)
        rcl_df[!, end+1] .= "RC + Learning"
        rcl_df[:, end+1] = 2005:2005+size(rcl_df, 1)-1
        rcl_t = convert(Matrix, rcl_df)

        figure_5bcd = vcat(uncert_t, learn_t, rc_t, rcl_t)
        figure_5bcd = convert(DataFrame, figure_5bcd)
        figure_5bcd = figure_5bcd[:, [7, 12, 13]]
        names!(figure_5bcd, [:scc_direct, :type, :year])

        CSV.write("generate_plots/data/figure_" * run[2] * ".csv", figure_5bcd)
    end

    ##################################
    ##################################
    # Figures 6a and 6b
    ##################################
    ##################################

    println("Generating data for Figures 6a and 6b.")

    loaded_results = Vector(undef, 8)
    thetas = [6.0, 8.0, 10.0, sqrt(10.)*10., 100.0, sqrt(10.)*100., 1000.0]

    loaded_results[1] = rc
    figure_6a = DataFrame(tax = loaded_results[1]["tax_channels"][:scc_direct][1],
                           penalty_parameter = loaded_results[1]["mparams"][:theta_robust])

    for (index, i) in enumerate(thetas)
        loaded_results[index+1] =
            BSON.load("results/rc_results_$(i).bson")
        push!(figure_6a, [loaded_results[index+1]["tax_channels"][:scc_direct][1], i])
    end

    CSV.write("generate_plots/data/figure_6a.csv", figure_6a)

    loaded_results = Vector(undef, 8)

    loaded_results[1] = rcl
    figure_6b = DataFrame(tax = loaded_results[1]["tax_channels"][:scc_direct][1],
                           penalty_parameter = loaded_results[1]["mparams"][:theta_robust])

    for (index, i) in enumerate(thetas)
        loaded_results[index+1] =
            BSON.load("results/rcl_results_$(i)_adapt_20.bson")
        push!(figure_6b, [loaded_results[index+1]["tax_channels"][:scc_direct][1], i])
    end

    CSV.write("generate_plots/data/figure_6b.csv", figure_6b)

    ##################################
    ##################################
    # Figure 8
    ##################################
    ##################################

    println("Generating data for Figure 8.")

    # Extract uncertanty to get mparams for simulating trajectories
    extract(uncert)
    mparams_in = mparams
    mparams_in[:num_runs] = 10000
    mparams_in[:sim_length] = 201
    exog_states = simulate_exog_vars(mparams)

    # Load different model-damage function combinations, compute ex-post welfare
    model_welfare = Array{Vector}(undef, 3, 4)
    models = Array{Dict}(undef, 3, 4)
    @everywhere include("tools.jl")

    for (r_index,run_type) in enumerate(["weitzman", "dietz", "extreme"])

        for (m_index, model_type) in
            enumerate(["uncert_results_dam_misspec_",
                "rc_results_dam_misspec_",
                "rcl_results_dam_misspec_",
                "learn_results_dam_misspec_"])

            println([r_index m_index])

            models[r_index,m_index] = BSON.load("results/" * model_type * run_type * ".bson")
            println(model_type * run_type)
            model_welfare[r_index,m_index] = ex_post_welfare(models[r_index, m_index], [50, 100, 150, 200], exog_states, mparams_in)

        end

    end

    # initialize state and arrays
    extract(uncert)
    coeff = c_store[1]
    if space[:d] < 4
        initial_state = [mparams[:kstart], mparams[:Mstart], mparams[:lossstart]]
    else
        initial_state = [mparams[:kstart], mparams[:Mstart], mparams[:expstart], mparams[:varstart], mparams[:lossstart]]
    end

    lump = Array{Vector}(undef, 3, 4)

    # marginal utility for translating into dollar terms
    MU(x) = x.^(mparams[:rho] - 1.)/mparams[:scaling_factor]
    dollars_translate = exog_states[:At][1].*exog_states[:Lt][1]./MU(consumption[1][1])

    # compute lump sum welfare
    for model_type = 2:4
        for run_type = 1:3

            # Lump sum
            lump[run_type,model_type-1] = 1e3*vec(-(model_welfare[run_type, 1] .-
                model_welfare[run_type,model_type].*models[run_type, model_type]["mparams"][:scaling_factor])).*dollars_translate

        end
    end

    uncert_tax = Dict{Symbol,Vector}()
    learn_tax = Dict{Symbol,Vector}()
    rc_tax = Dict{Symbol,Vector}()
    rcl_tax = Dict{Symbol,Vector}()
    for (key,value) in uncert["tax_channels"]
        if (key != :robust_insurance2)
            uncert_tax[key] = vec(uncert["tax_channels"][key])
            learn_tax[key] = vec(learn["tax_channels"][key])
            rc_tax[key] = vec(rc["tax_channels"][key])
            rcl_tax[key] = vec(rcl["tax_channels"][key])
        end

    end

    ### learning
    l_out = Dict{Symbol,Vector}()
    l_out[:weitz] = lump[1, 3]
    pushfirst!(l_out[:weitz], 0)
    l_out[:dietz] = lump[2, 3]
    pushfirst!(l_out[:dietz], 0)
    l_out[:extreme] = lump[3, 3]
    pushfirst!(l_out[:extreme], 0)
    l_df = convert(DataFrame, l_out)
    l_df[:, end+1] = "Learning"
    l_df[:, end+1] = (2005:50:2205) .- 2005
    l_df = convert(Matrix, l_df)

    ### robust control
    rc_out = Dict{Symbol,Vector}()
    rc_out[:weitz] = lump[1, 1]
    pushfirst!(rc_out[:weitz], 0)
    rc_out[:dietz] = lump[2, 1]
    pushfirst!(rc_out[:dietz], 0)
    rc_out[:extreme] = lump[3, 1]
    pushfirst!(rc_out[:extreme], 0)
    rc_df = convert(DataFrame, rc_out)
    rc_df[:, end+1] = "Robust Control"
    rc_df[:, end+1] = (2005:50:2205) .- 2005
    rc_df = convert(Matrix, rc_df)
    rc_out = Dict{Symbol,Vector}()

    ### robust control + learning
    rcl_out = Dict{Symbol,Vector}()
    rcl_out[:weitz] = lump[1, 2]
    pushfirst!(rcl_out[:weitz], 0)
    rcl_out[:dietz] = lump[2, 2]
    pushfirst!(rcl_out[:dietz], 0)
    rcl_out[:extreme] = lump[3, 2]
    pushfirst!(rcl_out[:extreme], 0)
    rcl_df = convert(DataFrame, rcl_out)
    rcl_df[:, end+1] = "RC + Learning"
    rcl_df[:, end+1] = (2005:50:2205) .- 2005
    rcl_df = convert(Matrix, rcl_df)

    figure_8 = vcat(l_df, rcl_df, rc_df)
    figure_8 = convert(DataFrame, figure_8)
    names!(figure_8, [:dietz, :extreme, :weitzman, :type, :year])

    CSV.write("generate_plots/data/figure_8.csv", figure_8)

    ##################################
    ##################################
    # Table 2
    ##################################
    ##################################

    println("Generating data for Table 2.")

    table_2_mean = DataFrame(tax_2055 = Vector{Int64}(undef, 4),
                            tax_2105 = Vector{Int64}(undef, 4),
                            co2_2105 = Vector{Int64}(undef, 4),
                            temperature_2105 = Vector{Float64}(undef, 4),
                            framework = ["uncertainty", "learning", "robust control", "robust control + learning"],
                            type = "mean")
    table_2_5th = DataFrame(tax_2055 = Vector{Int64}(undef, 4),
                            tax_2105 = Vector{Int64}(undef, 4),
                            co2_2105 = Vector{Int64}(undef, 4),
                            temperature_2105 = Vector{Float64}(undef, 4),
                            framework = ["uncertainty", "learning", "robust control", "robust control + learning"],
                            type = "5th percentile")
    table_2_95th = DataFrame(tax_2055 = Vector{Int64}(undef, 4),
                            tax_2105 = Vector{Int64}(undef, 4),
                            co2_2105 = Vector{Int64}(undef, 4),
                            temperature_2105 = Vector{Float64}(undef, 4),
                            framework = ["uncertainty", "learning", "robust control", "robust control + learning"],
                            type = "95th percentile")

    for (idx, model) in enumerate([uncert, learn, rc, rcl])

        # co2: 1 ppm / 2.12 gigatons of carbon
        co2_bounds = get_bounds(model["co2"], .05)./2.12
        table_2_5th.co2_2105[idx] = round(Int, co2_bounds[1][end])
        table_2_95th.co2_2105[idx] = round(Int, co2_bounds[2][end])
        table_2_mean.co2_2105[idx] = round(Int, mean(hcat(model["co2"]...), dims = 2)[end]./2.12)

        # temperature: mparams[:ccr] translates cumulative emissions to temperature
        temp_bounds = get_bounds(model["co2"], .05).*mparams[:ccr]
        table_2_5th.temperature_2105[idx] = round(temp_bounds[1][end], digits = 2)
        table_2_95th.temperature_2105[idx] = round(temp_bounds[2][end], digits = 2)
        table_2_mean.temperature_2105[idx] = round(mean(hcat(model["co2"]...), dims = 2)[end].*mparams[:ccr], digits = 2)

        # tax
        tax_bounds = get_bounds(model["channels_out_samples"][:scc_direct], .05)
        table_2_5th.tax_2055[idx] = round(Int, tax_bounds[1][51])
        table_2_95th.tax_2055[idx] = round(Int, tax_bounds[2][51])
        table_2_5th.tax_2105[idx] = round(Int, tax_bounds[1][end])
        table_2_95th.tax_2105[idx] = round(Int, tax_bounds[2][end])
        tax_mean = mean(hcat(model["channels_out_samples"][:scc_direct]...), dims = 2)
        table_2_mean.tax_2055[idx] = round(Int, tax_mean[51])
        table_2_mean.tax_2105[idx] = round(Int, tax_mean[end])

    end

    # aggregate and order table for writing to csv
    table_2 = vcat(table_2_mean, table_2_5th, table_2_95th)

    CSV.write("generate_plots/data/table_2.csv", table_2)

    ##################################
    ##################################
    # Table 3
    ##################################
    ##################################

    println("Generating data for Table 3.")

    # compute lump sum and balanced growth equivalent welfare gains
    gains_learn = ex_ante_welfare(uncert, learn)
    gains_rc = ex_ante_welfare(uncert, rc)
    gains_rcl = ex_ante_welfare(uncert, rcl)

    # Notes:
    # 1) lump_sum is already in billions of 2005 dollars
    # 2) 6.514 is 2005 population in billions so per capita is in dollars
    # 3) bge is in percent, i.e. 1.5 = 1.5%
    table_3 = DataFrame(lump_sum_2005 = [gains_learn[2], gains_rc[2], gains_rcl[2]],
                        per_capita_2005 = [gains_learn[2], gains_rc[2], gains_rcl[2]]./6.514,
                        bge_gain = [gains_learn[1], gains_rc[1], gains_rcl[1]],
                        framework = ["learning", "robust control", "robust control + learning"])

    CSV.write("generate_plots/data/table_3.csv", table_3)

    ##################################
    ##################################
    # Figure A1
    ##################################
    ##################################

    println("Generating data for Figure A1.")

    # Initial level of cumulative emissions
    cumul_emissions = 500.000

    # Annual emissions along the DICE 2016 optimal trajectory (from DICE 2016 spread sheet)
    flow_emissions_dice = [0 38.340 35.402 37.142 38.557 39.604 40.250 40.470 40.245 39.563 38.416 36.802 34.725 32.191 29.208 25.790 21.950 17.705 13.072 8.066]*5/3.666
    flow_emissions = vcat(repeat(flow_emissions_dice, 5)...)/5

    # 1% increase in cumulative CO2 each year
    flow_emissions_1p = similar(flow_emissions_dice)
    flow_emissions_1p .= 1.01

    # Cumulative annual emissions along DICE optimal
    cumul_emissions = cumsum(flow_emissions) .+ cumul_emissions

    # Cumulative annual emissions along 1% increase
    cumul_emissions_1p = cumprod(vcat(repeat(flow_emissions_1p, 5)...)).*cumul_emissions[1]
    pushfirst!(cumul_emissions_1p, 500.)

    # Put into 5 year timestep for DICE
    flow_emissions_1p_dice = similar(flow_emissions_dice)
    flow_emissions_1p_dice .= 1.01^5
    cumul_emissions_1p_dice = vec(cumprod(flow_emissions_1p_dice, dims = 2).*cumul_emissions[1])
    pushfirst!(cumul_emissions_1p_dice, 500.)
    flow_emissions_1p_dice = vec(flow_emissions_1p_dice .- 1.).*cumul_emissions_1p_dice[1:end-1]

    # Cumulative carbon temperature trajectory
    transient_climate_response = 0.0016
    ccr_temp = transient_climate_response*cumul_emissions
    ccr_temp = [ccr_temp (2005:2005+length(ccr_temp)-1)]
    ccr_temp_1p = transient_climate_response*cumul_emissions_1p
    ccr_temp_1p = [ccr_temp_1p (2005:2005+length(ccr_temp_1p)-1)]

    # DICE temperature trajectory
    dice_carbon = [.88 .196 0.; .12 .797 .001465; 0 .007 .998535]
    forcing_per_2xco2 = 3.681
    Mpre = 588.
    c1 = 0.101
    c3 = 0.088
    c4 = 0.025
    doubling = 3.1
    exog_force(t) = 0.5 + min((1/17)*(1. - 0.5)*(t-1), 0.5)

    co2 = Matrix{Float64}(undef, 3, 21)
    dice_temperature_opt = Matrix{Float64}(undef, 3, 21)
    dice_temperature_opt[3, :] = 2015:5:2115
    co2[:, 1] = [851.; 460.; 1740.]
    dice_temperature_opt[1:2, 1] = [0.85; 0.007]

    for d = 2:20

        #CO2, emissions are 10 times the emissions from the middle of the decade
        co2[:, d] = dice_carbon*co2[:, d-1] + [flow_emissions_dice[d]; 0; 0]

        # Forcing
        forcing = forcing_per_2xco2*log(co2[1, d]/Mpre)/log(2) + exog_force(d)

        # Temperature
        dice_temperature_opt[1, d] = dice_temperature_opt[1, d-1] +
            c1*(forcing - (forcing_per_2xco2/doubling)*dice_temperature_opt[1, d-1] +
            c3*(dice_temperature_opt[2, d-1] - dice_temperature_opt[1, d-1]))

        dice_temperature_opt[2, d] = (1 - c4)*dice_temperature_opt[2, d-1] + c4*dice_temperature_opt[1, d-1]

    end

    dice_temperature_opt = transpose(dice_temperature_opt)

    co2 = Matrix{Float64}(undef, 3, 21)
    dice_temperature_1p = Matrix{Float64}(undef, 3, 21)
    dice_temperature_1p[3, :] = 2015:5:2115
    co2[:, 1] = [851.; 460.; 1740.]
    dice_temperature_1p[1:2, 1] = [0.85; 0.007]

    for d = 2:20

        #CO2, emissions are 10 times the emissions from the middle of the decade
        co2[:, d] = dice_carbon*co2[:, d-1] + [flow_emissions_1p_dice[d]; 0; 0]

        # Forcing
        forcing = forcing_per_2xco2*log(co2[1, d]/Mpre)/log(2) + exog_force(d)

        # Temperature
        dice_temperature_1p[1, d] = dice_temperature_1p[1, d-1] +
            c1*(forcing - (forcing_per_2xco2/doubling)*dice_temperature_1p[1, d-1] +
            c3*(dice_temperature_1p[2, d-1] - dice_temperature_1p[1, d-1]))

        dice_temperature_1p[2, d] = (1 - c4)*dice_temperature_1p[2, d-1] + c4*dice_temperature_1p[1, d-1]

    end

    dice_temperature_1p = transpose(dice_temperature_1p)

    figure_a1_dice = DataFrame(year = dice_temperature_1p[:, 3],
                        temp_opt = dice_temperature_opt[:, 1],
                        temp_1p = dice_temperature_1p[:, 1])

    figure_a1_ccr = DataFrame(year = ccr_temp[11:end, 2],
                        ccr_opt = ccr_temp[11:end, 1],
                        ccr_1p = ccr_temp_1p[11:end-1, 1])

    CSV.write("generate_plots/data/figure_a1_dice.csv", figure_a1_dice)
    CSV.write("generate_plots/data/figure_a1_ccr.csv", figure_a1_ccr)

    ##################################
    ##################################
    # Figure A2
    ##################################
    ##################################

    println("Generating data for Figure A2.")

    runs_func(x, n) = x*ones(n, 1)

    bounds_main = learn
    bounds_previous_iter = BSON.load("results/learning_results_adapt_19.bson")

    # Code to get average relative change on capital lower bound (< 1%):
    # mean((bounds_main["smin"][1:300,1] - bounds_previous_iter["smin"][1:300,1])./bounds_previous_iter["smin"][1:300,1])

    extract(bounds_previous_iter)
    smax_previous_iter = bounds_previous_iter["smax"]
    smin_previous_iter = bounds_previous_iter["smin"]
    smax_previous_iter = convert(Matrix{Any}, [smax_previous_iter[1:101, :] similar(smax_previous_iter[1:101, 1:3])])
    smin_previous_iter = convert(Matrix{Any}, [smin_previous_iter[1:101, :] similar(smin_previous_iter[1:101, 1:3])])
    smax_previous_iter[:, end] .= "smax"
    smin_previous_iter[:, end] .= "smin"
    smax_previous_iter[:, end-1] = 2005:size(smax_previous_iter, 1) + 2005 - 1
    smin_previous_iter[:, end-1] = 2005:size(smin_previous_iter, 1) + 2005 - 1
    smax_previous_iter[:, end-2] .= "previous_iter"
    smin_previous_iter[:, end-2] .= "previous_iter"

    extract(bounds_main)
    smax = convert(Matrix{Any}, [bounds_main["smax"][1:101, :] similar(bounds_main["smax"][1:101, 1:3])])
    smin = convert(Matrix{Any}, [bounds_main["smin"][1:101, :] similar(bounds_main["smin"][1:101, 1:3])])
    smax[:, end] .= "smax"
    smin[:, end] .= "smin"
    smax[:, end-1] = 2005:size(smax, 1) + 2005 - 1
    smin[:, end-1] = 2005:size(smin, 1) + 2005 - 1
    smax[:, end-2] .= "main"
    smin[:, end-2] .= "main"

    capital = vcat(hcat(bounds_main["capital"]...)[1:101, 1:1000]...)
    capital = convert(Matrix{Any}, [capital repeat(similar(capital), 1, 3)])
    capital[:, end] .= "capital"
    capital[:, end-1] = repeat(2005:size(smax, 1) + 2005 - 1, 1000)
    capital[:, end-2] = vcat([runs_func(run, size(smax, 1)) for run = 1:1000]...)

    co2 = vcat(hcat(bounds_main["co2"]...)[1:101, 1:1000]...)
    co2 = convert(Matrix{Any}, [co2 repeat(similar(co2), 1, 3)])
    co2[:, end] .= "co2"
    co2[:, end-1] = repeat(2005:size(smax, 1) + 2005 - 1, 1000)
    co2[:, end-2] = vcat([runs_func(run, size(smax, 1)) for run = 1:1000]...)

    damage = vcat(hcat(bounds_main["damage"]...)[1:101, 1:1000]...)
    damage = convert(Matrix{Any}, [damage repeat(similar(damage), 1, 3)])
    damage[:, end] .= "damage"
    damage[:, end-1] = repeat(2005:size(smax, 1) + 2005 - 1, 1000)
    damage[:, end-2] = vcat([runs_func(run, size(smax, 1)) for run = 1:1000]...)

    prior_loc = vcat(hcat(bounds_main["prior_loc"]...)[1:101, 1:1000]...)
    prior_loc = convert(Matrix{Any}, [prior_loc repeat(similar(prior_loc), 1, 3)])
    prior_loc[:, end] .= "prior_loc"
    prior_loc[:, end-1] = repeat(2005:size(smax, 1) + 2005 - 1, 1000)
    prior_loc[:, end-2] = vcat([runs_func(run, size(smax, 1)) for run = 1:1000]...)

    prior_scale = vcat(hcat(bounds_main["prior_scale"]...)[1:101, 1:1000]...)
    prior_scale = convert(Matrix{Any}, [prior_scale repeat(similar(prior_scale), 1, 3)])
    prior_scale[:, end] .= "prior_scale"
    prior_scale[:, end-1] = repeat(2005:size(smax, 1) + 2005 - 1, 1000)
    prior_scale[:, end-2] = vcat([runs_func(run, size(smax, 1)) for run = 1:1000]...)

    figure_a2 = vcat(smax, smin, smax_previous_iter, smin_previous_iter)
    figure_a2_states = vcat(capital, co2, damage, prior_loc, prior_scale)

    figure_a2 = convert(DataFrame, figure_a2)
    figure_a2_states = convert(DataFrame, figure_a2_states)

    names!(figure_a2, [:capital, :co2, :loc, :scale, :damage, :it, :year, :type])
    names!(figure_a2_states, [:value, :run, :year, :state])

    CSV.write("generate_plots/data/figure_a2_bounds.csv", figure_a2)
    CSV.write("generate_plots/data/figure_a2_states.csv", figure_a2_states)

    ##################################
    ##################################
    # Table A1
    ##################################
    ##################################

    table_a1 = DataFrame(A0 = mparams[:A0],
                        gA0 = mparams[:gA0],
                        delta_A = mparams[:dA],
                        L_0 = mparams[:L0],
                        L_inf = mparams[:Linf],
                        delta_L = mparams[:dL],
                        sigma_0 = mparams[:sigma0],
                        g_sigma0 = mparams[:gsigma0],
                        delta_sigma = mparams[:dsigma],
                        a_0 = mparams[:a0],
                        a_1 = mparams[:a1],
                        a_2 = mparams[:a2],
                        g_psi = mparams[:gPsi],
                        B_0 = mparams[:B0],
                        g_B = mparams[:gB],
                        kappa = mparams[:kappa],
                        delta_k = mparams[:deltak],
                        E_0 = mparams[:Mstart],
                        rho = mparams[:delta],
                        eta = 1 - mparams[:rho], # mparams[:rho] is old notation for (1-eta)
                        zeta = mparams[:ccr],
                        k_0 = mparams[:kstart],
                        mu_0 = mparams[:expstart],
                        cap_sigma_0 = mparams[:varstart],
                        joint_d1_omega_loc = mparams[:coeff_loc], # d1 and omega are joined into one distribution for efficiency, initialization.jl builds this on lines 72-72
                        joint_d1_omega_scale = mparams[:coeff_scale])

    CSV.write("generate_plots/data/table_a1.csv", table_a1)

    ##################################
    ##################################
    # Table A4
    ##################################
    ##################################

    println("Generating data for Table A4.")

    #############################
    ### First column: Uncertainty

    main_approx = uncert
    approx_1 = BSON.load("results/uncert_results_coll_1.bson")
    approx_2 = BSON.load("results/uncert_results_coll_2.bson")
    approx_3 = BSON.load("results/uncert_results_coll_3.bson")
    approx_4 = BSON.load("results/uncert_results_quad.bson")
    approx_5 = BSON.load("results/uncert_results_wide.bson")
    approx_6 = BSON.load("results/uncert_results_coll_4.bson")

    tax_store = Vector{Array}(undef, 7)
    cons_store = Vector{Array}(undef, 7)
    for (idx, type1) in enumerate([approx_1 approx_2 approx_3 approx_4 approx_5 approx_6 main_approx])
        extract(type1)
        tax_store[idx] = tax_channels[:scc_direct]
        cons_store[idx] = mean(hcat(consumption...), dims=2)
    end

    # Table output:
    # tax: tax error
    # cons: consumption error
    # mean: mean relative error
    # max: max relative error
    # coll: collocation nodes (top panel)
    # quad: quadrature (second panel)
    # domain: expand domain (third panel)
    # double coll: collocation error from increasing capital by 1 vs decreasing capital by 1 (footnote)

    table_a4_column_1 = DataFrame(tax_coll_error_mean = mean(abs.((tax_store[end] .- hcat(tax_store[1:3]...))./tax_store[end])),
                                tax_coll_error_max = maximum(abs.((tax_store[end] .- hcat(tax_store[1:3]...))./tax_store[end])),
                                tax_quad_error_mean = mean(abs.((tax_store[end] .- tax_store[4])./tax_store[end])),
                                tax_quad_error_max = maximum(abs.((tax_store[end] .- tax_store[4])./tax_store[end])),
                                tax_domain_error_mean = mean(abs.((tax_store[end] .- tax_store[5])./tax_store[end])),
                                tax_domain_error_max = maximum(abs.((tax_store[end] .- tax_store[5])./tax_store[end])),
                                tax_double_coll_error_mean = mean(abs.((tax_store[1] .- tax_store[6])./tax_store[1])),
                                tax_double_coll_error_max = maximum(abs.((tax_store[1] .- tax_store[6])./tax_store[1])),
                                cons_coll_error_mean = mean(abs.((cons_store[end] .- hcat(cons_store[1:3]...))./cons_store[end])),
                                cons_coll_error_max = maximum(abs.((cons_store[end] .- hcat(cons_store[1:3]...))./cons_store[end])),
                                cons_quad_error_mean = mean(abs.((cons_store[end] .- cons_store[4])./cons_store[end])),
                                cons_quad_error_max = maximum(abs.((cons_store[end] .- cons_store[4])./cons_store[end])),
                                cons_domain_error_mean = mean(abs.((cons_store[end] .- cons_store[5])./cons_store[end])),
                                cons_domain_error_max = maximum(abs.((cons_store[end] .- cons_store[5])./cons_store[end])),
                                cons_double_coll_error_mean = mean(abs.((cons_store[1] .- cons_store[6])./cons_store[1])),
                                cons_double_coll_error_max = maximum(abs.((cons_store[1] .- cons_store[6])./cons_store[1])))

    CSV.write("generate_plots/data/table_a4_uncert.csv", table_a4_column_1)

    #############################
    ### Second column: learning

    main_approx = learn
    approx_1 = BSON.load("results/learning_results_coll_1.bson")
    approx_2 = BSON.load("results/learning_results_coll_2.bson")
    approx_3 = BSON.load("results/learning_results_coll_3.bson")
    approx_4 = BSON.load("results/learning_results_coll_4.bson")
    approx_5 = BSON.load("results/learning_results_coll_5.bson")
    approx_6 = BSON.load("results/learning_results_quad.bson")
    approx_7 = BSON.load("results/learning_results_adapt_10.bson")
    approx_8 = BSON.load("results/learning_results_coll_6.bson")

    tax_store = Vector{Array}(undef, 9)
    cons_store = Vector{Array}(undef, 9)
    for (idx, type) in enumerate([approx_1 approx_2 approx_3 approx_4 approx_5 approx_6 approx_7 approx_8 main_approx])
        extract(type)
        tax_store[idx] = tax_channels[:scc_direct]
        cons_store[idx] = mean(hcat(consumption...), dims=2)
    end

    # Table output:
    # tax: tax error
    # cons: consumption error
    # mean: mean relative error
    # max: max relative error
    # coll: collocation nodes (top panel)
    # quad: quadrature (second panel)
    # adapt: 10 adaptive algorithms (fourth panel)
    # double coll: collocation error from increasing capital by 1 vs decreasing capital by 1 (footnote)

    table_a4_column_2 = DataFrame(tax_coll_error_mean = mean(abs.((tax_store[end] .- hcat(tax_store[1:5]...))./tax_store[end])),
                                tax_coll_error_max = maximum(abs.((tax_store[end] .- hcat(tax_store[1:5]...))./tax_store[end])),
                                tax_quad_error_mean = mean(abs.((tax_store[end] .- tax_store[6])./tax_store[end])),
                                tax_quad_error_max = maximum(abs.((tax_store[end] .- tax_store[6])./tax_store[end])),
                                tax_adapt_error_mean = mean(abs.((tax_store[end] .- tax_store[7])./tax_store[end])),
                                tax_adapt_error_max = maximum(abs.((tax_store[end] .- tax_store[7])./tax_store[end])),
                                tax_double_coll_error_mean = mean(abs.((tax_store[1] .- tax_store[8])./tax_store[1])),
                                tax_double_coll_error_max = maximum(abs.((tax_store[1] .- tax_store[8])./tax_store[1])),
                                cons_coll_error_mean = mean(abs.((cons_store[end] .- hcat(cons_store[1:5]...))./cons_store[end])),
                                cons_coll_error_max = maximum(abs.((cons_store[end] .- hcat(cons_store[1:5]...))./cons_store[end])),
                                cons_quad_error_mean = mean(abs.((cons_store[end] .- cons_store[6])./cons_store[end])),
                                cons_quad_error_max = maximum(abs.((cons_store[end] .- cons_store[6])./cons_store[end])),
                                cons_adapt_error_mean = mean(abs.((cons_store[end] .- cons_store[7])./cons_store[end])),
                                cons_adapt_error_max = maximum(abs.((cons_store[end] .- cons_store[7])./cons_store[end])),
                                cons_double_coll_error_mean = mean(abs.((cons_store[1] .- cons_store[8])./cons_store[1])),
                                cons_double_coll_error_max = maximum(abs.((cons_store[1] .- cons_store[8])./cons_store[1])))

    CSV.write("generate_plots/data/table_a4_learn.csv", table_a4_column_2)

    #############################
    ### Third column: Robust control

    main_approx = rc
    approx_1 = BSON.load("results/rc_results_coll_1.bson")
    approx_2 = BSON.load("results/rc_results_coll_2.bson")
    approx_3 = BSON.load("results/rc_results_coll_3.bson")
    approx_4 = BSON.load("results/rc_results_3.9_quad.bson")
    approx_5 = BSON.load("results/rc_results_3.9_wide.bson")
    approx_6 = BSON.load("results/rc_results_coll_4.bson")

    tax_store = Vector{Array}(undef, 7)
    cons_store = Vector{Array}(undef, 7)
    for (idx, type1) in enumerate([approx_1 approx_2 approx_3 approx_4 approx_5 approx_6 main_approx])
        extract(type1)
        tax_store[idx] = tax_channels[:scc_direct]
        cons_store[idx] = mean(hcat(consumption...), dims=2)
    end

    # Table output:
    # tax: tax error
    # cons: consumption error
    # mean: mean relative error
    # max: max relative error
    # coll: collocation nodes (top panel)
    # quad: quadrature (second panel)
    # domain: expand domain (third panel)
    # double coll: collocation error from increasing capital by 1 vs decreasing capital by 1 (footnote)

    table_a4_column_3 = DataFrame(tax_coll_error_mean = mean(abs.((tax_store[end] .- hcat(tax_store[1:3]...))./tax_store[end])),
                                tax_coll_error_max = maximum(abs.((tax_store[end] .- hcat(tax_store[1:3]...))./tax_store[end])),
                                tax_quad_error_mean = mean(abs.((tax_store[end] .- tax_store[4])./tax_store[end])),
                                tax_quad_error_max = maximum(abs.((tax_store[end] .- tax_store[4])./tax_store[end])),
                                tax_domain_error_mean = mean(abs.((tax_store[end] .- tax_store[5])./tax_store[end])),
                                tax_domain_error_max = maximum(abs.((tax_store[end] .- tax_store[5])./tax_store[end])),
                                tax_double_coll_error_mean = mean(abs.((tax_store[1] .- tax_store[6])./tax_store[1])),
                                tax_double_coll_error_max = maximum(abs.((tax_store[1] .- tax_store[6])./tax_store[1])),
                                cons_coll_error_mean = mean(abs.((cons_store[end] .- hcat(cons_store[1:3]...))./cons_store[end])),
                                cons_coll_error_max = maximum(abs.((cons_store[end] .- hcat(cons_store[1:3]...))./cons_store[end])),
                                cons_quad_error_mean = mean(abs.((cons_store[end] .- cons_store[4])./cons_store[end])),
                                cons_quad_error_max = maximum(abs.((cons_store[end] .- cons_store[4])./cons_store[end])),
                                cons_domain_error_mean = mean(abs.((cons_store[end] .- cons_store[5])./cons_store[end])),
                                cons_domain_error_max = maximum(abs.((cons_store[end] .- cons_store[5])./cons_store[end])),
                                cons_double_coll_error_mean = mean(abs.((cons_store[1] .- cons_store[6])./cons_store[1])),
                                cons_double_coll_error_max = maximum(abs.((cons_store[1] .- cons_store[6])./cons_store[1])))

    CSV.write("generate_plots/data/table_a4_rc.csv", table_a4_column_3)

    #############################
    ### Fourth column: RC+L

    main_approx = rcl
    approx_1 = BSON.load("results/rcl_results_coll_1.bson")
    approx_2 = BSON.load("results/rcl_results_coll_2.bson")
    approx_3 = BSON.load("results/rcl_results_coll_3.bson")
    approx_4 = BSON.load("results/rcl_results_coll_4.bson")
    approx_5 = BSON.load("results/rcl_results_coll_5.bson")
    approx_6 = BSON.load("results/rcl_results_3.9_quad.bson")
    approx_7 = BSON.load("results/rcl_results_3.9_adapt_10.bson")
    approx_8 = BSON.load("results/rcl_results_coll_6.bson")

    tax_store = Vector{Array}(undef, 9)
    cons_store = Vector{Array}(undef, 9)
    for (idx, type) in enumerate([approx_1 approx_2 approx_3 approx_4 approx_5 approx_6 approx_7 approx_8 main_approx])
        extract(type)
        tax_store[idx] = tax_channels[:scc_direct]
        cons_store[idx] = mean(hcat(consumption...), dims=2)
    end

    # Table output:
    # tax: tax error
    # cons: consumption error
    # mean: mean relative error
    # max: max relative error
    # coll: collocation nodes (top panel)
    # quad: quadrature (second panel)
    # adapt: 10 adaptive algorithms (fourth panel)
    # double coll: collocation error from increasing capital by 1 vs decreasing capital by 1 (footnote)

    table_a4_column_4 = DataFrame(tax_coll_error_mean = mean(abs.((tax_store[end] .- hcat(tax_store[1:5]...))./tax_store[end])),
                                tax_coll_error_max = maximum(abs.((tax_store[end] .- hcat(tax_store[1:5]...))./tax_store[end])),
                                tax_quad_error_mean = mean(abs.((tax_store[end] .- tax_store[6])./tax_store[end])),
                                tax_quad_error_max = maximum(abs.((tax_store[end] .- tax_store[6])./tax_store[end])),
                                tax_adapt_error_mean = mean(abs.((tax_store[end] .- tax_store[7])./tax_store[end])),
                                tax_adapt_error_max = maximum(abs.((tax_store[end] .- tax_store[7])./tax_store[end])),
                                tax_double_coll_error_mean = mean(abs.((tax_store[1] .- tax_store[8])./tax_store[1])),
                                tax_double_coll_error_max = maximum(abs.((tax_store[1] .- tax_store[8])./tax_store[1])),
                                cons_coll_error_mean = mean(abs.((cons_store[end] .- hcat(cons_store[1:5]...))./cons_store[end])),
                                cons_coll_error_max = maximum(abs.((cons_store[end] .- hcat(cons_store[1:5]...))./cons_store[end])),
                                cons_quad_error_mean = mean(abs.((cons_store[end] .- cons_store[6])./cons_store[end])),
                                cons_quad_error_max = maximum(abs.((cons_store[end] .- cons_store[6])./cons_store[end])),
                                cons_adapt_error_mean = mean(abs.((cons_store[end] .- cons_store[7])./cons_store[end])),
                                cons_adapt_error_max = maximum(abs.((cons_store[end] .- cons_store[7])./cons_store[end])),
                                cons_double_coll_error_mean = mean(abs.((cons_store[1] .- cons_store[8])./cons_store[1])),
                                cons_double_coll_error_max = maximum(abs.((cons_store[1] .- cons_store[8])./cons_store[1])))

    # Vector of Table A4 entries
    table_a4 = [table_a4_column_1,
                table_a4_column_2,
                table_a4_column_3,
                table_a4_column_4]

    CSV.write("generate_plots/data/table_a4_rcl.csv", table_a4_column_4)

    println("Done generating data. CSV files can be found in '/generate_plots/data'.")

end

max_stars_repo_name: irudik/optimal-climate-policy-aej

max_issues_repo_name: irudik/optimal-climate-policy-aej

max_forks_repo_name: irudik/optimal-climate-policy-aej