text:
module PkgBenchUtils

using PkgBenchmark: judge, export_markdown
using Setfield
import GitHub

include("ir_dump.jl")

struct Results
    results
    script
    posted

    Results(results, script = nothing, posted = nothing) =
        new(results, script, posted)
end

function Base.show(io::IO, r::Results)
    print(io, "PkgBenchUtils.Results")
    if r.script !== nothing
        print(io, " [", r.script, "]")
    end
    println(io)
    if r.posted !== nothing
        println(io, "URL: ", url(r))
    end
    show(io, r.results)
end

guess_package() = basename(pwd())

_common_docs = """
- `config::Union{Nothing, BenchmarkConfig}`: Configuration to be used
  if `target` or `baseline` is a `String`.
- `target`: Default to `"HEAD"`.  (Passed to `PkgBenchmark.judge`)
- `baseline`: Default to `"HEAD^"`.  (Passed to `PkgBenchmark.judge`)
- `script::String`: Script from which benchmark `SUITE` is loaded.
  (Passed to `PkgBenchmark.judge`)
- Other keyword arguments are also passed to `PkgBenchmark.judge`.
"""

"""
    _judge([pkg::String]; <keyword arguments>) :: Results

# Keyword Arguments
$_common_docs
"""
function _judge(package_name::String = guess_package();
                config = nothing,
                target = "HEAD",
                baseline = "HEAD^",
                script = nothing,
                kwargs...)

    if config !== nothing
        (target isa String) && (target = @set config.id = target)
        (baseline isa String) && (baseline = @set config.id = baseline)
    end

    return Results(
        judge(package_name, target, baseline;
              script = script,
              kwargs...),
        script,
    )
end

as_markdown(results::Results) = as_markdown(results.results)

function as_markdown(results)
    io = IOBuffer()
    export_markdown(io, results)
    seek(io, 0)
    return Markdown.parse(io)
end

"""
    show_judge([pkg::String]; <keyword arguments>) :: Results

# Keyword Arguments
$_common_docs
"""
function show_judge(args...; kwargs...)
    results = _judge(args...; kwargs...)
    println()
    display(as_markdown(results))
    println()
    return results
end

_post_docs = """
- `auth`: authentication object created by `GitHub.authenticate`.
  Default to `GitHub.authenticate(ENV["PKGBENCHUTILS_GITHUB_AUTH"])`.
- `public::Bool`: post to public gist if true.
"""

"""
    post_results(results::Results; <keyword arguments>) :: Results

https://juliaci.github.io/PkgBenchmark.jl/stable/export_markdown.html

# Keyword Arguments
$_post_docs
"""
function post_results(results::Results; kwargs...)
    posted = post_results(results.results, results.script; kwargs...)
    return @set results.posted = posted
end

function post_results(results, script;
                      auth = GitHub.authenticate(ENV["PKGBENCHUTILS_GITHUB_AUTH"]),
                      public = false,
                      kwargs...)
    package_name = results.baseline_results.name

    gist_json = Dict(
        :description => "A benchmark for $package_name",
        :public => public,
        :files => Dict(
            "benchmark.md" => Dict(
                :content => sprint(export_markdown, results),
            ),
        )
    )
    if script !== nothing
        gist_json[:files]["script.jl"] = Dict(
            :content => readstring(script),
        )
    end

    return GitHub.create_gist(params = gist_json;
                              auth = auth,
                              kwargs...)
end

"""
    post_judge([pkg::String]; <keyword arguments>) :: Results

Run `judge` and then `post_results`.

# Keyword Arguments
- `open::Bool`: open posted gist in browser.
$_post_docs
$_common_docs
"""
function post_judge(args...;
                    auth = GitHub.authenticate(ENV["PKGBENCHUTILS_GITHUB_AUTH"]),
                    public = false,
                    open = false,
                    kwargs...)

    results = _judge(
        args...;
        kwargs...)

    results = post_results(results; public=public, auth=auth)
    if open
        open_url(results)
    end

    return results
end

url(results) = get(results.posted.html_url)
open_url(results) = run(`xdg-open $(url(results))`)

end # module


text:
cubical_benchmarks=[
    Benchmark("lena512.dipha"; dim_max=1),
    Benchmark("lena1999x999.dipha"; dim_max=1),
    Benchmark("bonsai64.dipha"; dim_max=2),
    Benchmark("bonsai128.dipha"; dim_max=2),
    Benchmark("head128.dipha"; dim_max=2),
]

rips_benchmarks = [
    Benchmark("sphere_3_192.dist"; dim_max=2),
    Benchmark("o3_1024.pts"; dim_max=3, threshold=1.8, extra=(;sparse=true)),
    Benchmark("o3_4096.pts"; dim_max=3, threshold=1.4, extra=(;sparse=true)),
    Benchmark("dragon2000.pts"; dim_max=1),
    Benchmark("fract-r.dist"; dim_max=2),
    Benchmark("random16.dist"; dim_max=2),
]

alpha_rips_benchmarks = [
    Benchmark("alpha_torus_10_000.spdist"; dim_max=2),
    Benchmark("alpha_5_sphere_1000.spdist"; dim_max=5),
    Benchmark("alpha_4_sphere_2000.spdist"; dim_max=4),
    Benchmark("alpha_3_sphere_3000.spdist"; dim_max=3),
    Benchmark("alpha_dragon_2000.spdist"; dim_max=2),
]

eirene_benchmarks = [
    Benchmark("gcycle.dist"; dim_max=3),
    Benchmark("celegans.dist"; dim_max=2),
    Benchmark("dragon1000.dist"; dim_max=1),
    Benchmark("hiv.dist"; dim_max=1),
    Benchmark("o3_1024.pts"; dim_max=3, threshold=1.8, extra=(;sparse=true)),
    Benchmark("random16.dist"; dim_max=7),
]

involuted_benchmarks = [
    Benchmark("gcycle.dist"; dim_max=3),
    Benchmark("celegans.dist"; dim_max=2),
    Benchmark("dragon1000.dist"; dim_max=1),
    Benchmark("hiv.dist"; dim_max=1),
    Benchmark("lena2048.dipha"; dim_max=1),
    Benchmark("bonsai128.dipha"; dim_max=2),
    Benchmark("alpha_dragon2000.alpha"; dim_max=2),
    Benchmark("o3_1024.pts"; dim_max=3, threshold=1.8, extra=(;sparse=true)),
    Benchmark("o3_4096.pts"; dim_max=3, threshold=1.4, extra=(;sparse=true)),
]

test_ripser = [
    Benchmark("lena512.dipha"; dim_max=1),
    Benchmark("alpha_dragon_2000.spdist"; dim_max=2),
    Benchmark("gcycle.dist"; dim_max=2),
]

test_eirene = [
    Benchmark("gcycle.dist"; dim_max=2),
]

test_involuted = [
    Benchmark("alpha_dragon_2000.spdist"; dim_max=2),
]

function setup_suite(type, benchmarks; sparse=false)
    suite = BenchmarkGroup()
    if type == :ripser
        suite["ripser"] = BenchmarkGroup()
        suite["ripserer"] = BenchmarkGroup()
        if sparse
            suite["ripserer"]["sparse"] = BenchmarkGroup()
            suite["ripserer"]["dense"] = BenchmarkGroup()
        end
        for b in benchmarks
            suite["ripser"][b.name] = ripser_benchmark(b)
            if sparse
                suite["ripserer"]["dense"][b.name] = ripserer_benchmark(
                    b; from_file=true, sparse=false
                )
                suite["ripserer"]["sparse"][b.name] = ripserer_benchmark(
                    b; from_file=true, sparse=true
                )
            else
                suite["ripserer"][b.name] = ripserer_benchmark(b; from_file=true)
            end
        end
    elseif type == :eirene
        suite["eirene"] = BenchmarkGroup()
        suite["ripserer"] = BenchmarkGroup()
        for b in benchmarks
            suite["eirene"][b.name] = eirene_benchmark(b)
            suite["ripserer"][b.name] = ripserer_benchmark(b; alg=:involuted)
        end
    elseif type == :involuted
        suite["cohomology"] = BenchmarkGroup()
        suite["homology"] = BenchmarkGroup()
        suite["involuted"] = BenchmarkGroup()
        for b in benchmarks
            suite["cohomology"][b.name] = ripserer_benchmark(b; extras=true)
            suite["involuted"][b.name] = ripserer_benchmark(b; extras=true, alg=:involuted)
            if b.complex ≠ Rips
                suite["homology"][b.name] = ripserer_benchmark(b; extras=true, alg=:homology)
            end
        end
    end
    return suite
end

function setup_all()
    suite = BenchmarkGroup()
    suite["ripser"] = BenchmarkGroup()
    suite["ripser"]["rips"] = setup_suite(:ripser, rips_benchmarks; sparse=true)
    suite["ripser"]["cubical"] = setup_suite(:ripser, cubical_benchmarks)
    suite["ripser"]["alpha_rips"] = setup_suite(:ripser, alpha_rips_benchmarks)
    suite["eirene"] = setup_suite(:eirene, eirene_benchmarks)
    suite["involuted"] = setup_suite(:involuted, involuted_benchmarks)

    suite["test"] = BenchmarkGroup()
    suite["test"]["eirene"] = setup_suite(:eirene, test_eirene)
    suite["test"]["ripser"] = setup_suite(:ripser, test_ripser)
    suite["test"]["involuted"] = setup_suite(:involuted, test_involuted)

    return suite
end

ext: jl


text:
# rules.jl - differentiation rules for basic operations.

## basic rules

@scalardiff (-x::Number) 1 -1

# product

@scalardiff *(x::Number, y::Number) 1 y
@scalardiff *(x::Number, y::Number) 2 x

# elementwise product

@scalardiff .*(x::Number, y::Number) 1 y
@scalardiff .*(x::Number, y::Number) 2 x

# other arithmetic operations

@scalardiff (x::Number ^ n::Number) 1 (n * x.^(n-1))
@scalardiff (a::Number ^ x::Number) 2 (log(a) * a.^x)
@scalardiff (x::Number .^ n::Number) 1 (n * x.^(n-1))
@scalardiff (a::Number .^ x::Number) 2 (log(a) * a.^x)

@scalardiff (x::Number / y::Number) 1 (1 / y)
@scalardiff (n::Number / x::Real) 2 (-n * x .^ -2.0)
@scalardiff (x::Number ./ y::Number) 1 (1 / y)
@scalardiff (n::Number ./ x::Real) 2 (-n * x .^ -2.0)

@scalardiff (x::Any + y::Any) 1 1
@scalardiff (x::Any + y::Any) 2 1
@scalardiff (x::Any + y::Any + z::Any) 1 1
@scalardiff (x::Any + y::Any + z::Any) 2 1
@scalardiff (x::Any + y::Any + z::Any) 3 1
@scalardiff (w::Any + x::Any + y::Any + z::Any) 1 1
@scalardiff (w::Any + x::Any + y::Any + z::Any) 2 1
@scalardiff (w::Any + x::Any + y::Any + z::Any) 3 1
@scalardiff (w::Any + x::Any + y::Any + z::Any) 4 1

@scalardiff (x::Any .+ y::Any) 1 1
@scalardiff (x::Any .+ y::Any) 2 1

@scalardiff (x::Any - y::Any) 1 1
@scalardiff (x::Any - y::Any) 2 -1

@scalardiff (x::Any .- y::Any) 1 1
@scalardiff (x::Any .- y::Any) 2 -1

@scalardiff sum(x::Number) 1 1

# trigonomeric functions

@scalardiff sin(x::Number) 1 cos(x)
@scalardiff cos(x::Number) 1 -sin(x)
@scalardiff tan(x::Number) 1 (1. + tan(x)  * tan(x))
@scalardiff sinh(x::Number) 1 cosh(x)
@scalardiff cosh(x::Number) 1 sinh(x)
@scalardiff tanh(x::Number) 1 (1.0 .- tanh(x) .* tanh(x))
@scalardiff asin(x::Number) 1 (1 ./ sqrt(1 .- x .* x))
@scalardiff acos(x::Number) 1 (1 ./ sqrt(1 .- x .* x))
@scalardiff atan(x::Number) 1 (1 ./ (1 .+ x .* x))

# sqrt

@scalardiff sqrt(x::Number) 1 (0.5 * x .^ (-0.5))

# exp, log

@scalardiff exp(x::Number) 1 exp(x)
@scalardiff log(x::Number) 1 (1 ./ x)

# abs

@scalardiff abs(x::Number) 1 (sign(x) .* x)
@scalardiff abs2(x::Number) 1 (2 .* abs(x) .* sign(x) .* x)

# min, max

@scalardiff max(x::Number, y::Number) 1 (x > y)
@scalardiff max(x::Number, y::Number) 2 (y > x)

@scalardiff min(x::Number, y::Number) 1 (x < y)
@scalardiff min(x::Number, y::Number) 2 (y < x)

@scalardiff sign(x::Any) 1 0.

# transpose

@scalardiff transpose(x::Number) 1 1

@scalardiff size(x::Any) 1 0.0
@scalardiff size(x::Any, y::Any) 1 0.0
@scalardiff size(x::Any, y::Any) 2 0.0


text:
function input(
    equation::String;
    data::Union{Array{Float64}, Array{Float32}, Array{Union{Float32, Missing}}, Array{Union{Float64, Missing}}, Tuple, DataFrame, Nothing} = nothing,
    datanames::Union{Array, Array{Symbol, 1}, Nothing}=nothing,
    method::Union{Symbol, String}=METHOD_DEFAULT,
    intercept::Bool=INTERCEPT_DEFAULT,
    panel::Union{Symbol, String, Nothing}=PANEL_DEFAULT,
    time::Union{Symbol, String, Nothing}=TIME_DEFAULT,
    seasonaladjustment::Union{Dict, Array, Nothing}=SEASONALADJUSTMENT_DEFAULT,
    removeoutliers::Bool=REMOVEOUTLIERS_DEFAULT,
    removemissings::Bool=REMOVEMISSINGS_DEFAULT
    )

    return input(
        equation,
        data,
        datanames=datanames,
        method=method,
        intercept=intercept,
        panel=panel,
        time=time,
        seasonaladjustment=seasonaladjustment,
        removeoutliers=removeoutliers,
        removemissings=removemissings
    )
end

function input(
    equation::Array{String};
    data::Union{Array{Float64}, Array{Float32}, Array{Union{Float32, Missing}}, Array{Union{Float64, Missing}}, Tuple, DataFrame, Nothing} = nothing,
    datanames::Union{Array, Array{Symbol, 1}, Nothing}=nothing,
    method::Union{Symbol, String}=METHOD_DEFAULT,
    intercept::Bool=INTERCEPT_DEFAULT,
    panel::Union{Symbol, String, Nothing}=PANEL_DEFAULT,
    time::Union{Symbol, String, Nothing}=TIME_DEFAULT,
    seasonaladjustment::Union{Dict, Array, Nothing}=SEASONALADJUSTMENT_DEFAULT,
    removeoutliers::Bool=REMOVEOUTLIERS_DEFAULT,
    removemissings::Bool=REMOVEMISSINGS_DEFAULT
    )

    return input(
        equation,
        data,
        datanames=datanames,
        method=method,
        intercept=intercept,
        panel=panel,
        time=time,
        removeoutliers=removeoutliers,
        seasonaladjustment=seasonaladjustment,
        removemissings=removemissings
    )
end

function input(
    equation::Array{Symbol};
    data::Union{Array{Float64}, Array{Float32}, Array{Union{Float32, Missing}}, Array{Union{Float64, Missing}}, Tuple, DataFrame, Nothing} = nothing,
    datanames::Union{Array, Array{Symbol, 1}, Nothing}=nothing,
    method::Union{Symbol, String}=METHOD_DEFAULT,
    intercept::Bool=INTERCEPT_DEFAULT,
    panel::Union{Symbol, String, Nothing}=PANEL_DEFAULT,
    time::Union{Symbol, String, Nothing}=TIME_DEFAULT,
    seasonaladjustment::Union{Dict, Array, Nothing}=SEASONALADJUSTMENT_DEFAULT,
    removeoutliers::Bool=REMOVEOUTLIERS_DEFAULT,
    removemissings::Bool=REMOVEMISSINGS_DEFAULT
    )

    return input(
        equation,
        data,
        datanames=datanames,
        method=method,
        intercept=intercept,
        panel=panel,
        time=time,
        seasonaladjustment=seasonaladjustment,
        removeoutliers=removeoutliers,
        removemissings=removemissings
    )
end

function input(
    equation::String,
    data::Union{Array{Float64}, Array{Float32}, Array{Union{Float32, Missing}}, Array{Union{Float64, Missing}}, Tuple, DataFrame, Nothing};
    datanames::Union{Array, Array{Symbol, 1}, Nothing}=nothing,
    method::Union{Symbol, String}=METHOD_DEFAULT,
    intercept::Bool=INTERCEPT_DEFAULT,
    panel::Union{Symbol, String, Nothing}=PANEL_DEFAULT,
    time::Union{Symbol, String, Nothing}=TIME_DEFAULT,
    seasonaladjustment::Union{Dict, Array, Nothing}=SEASONALADJUSTMENT_DEFAULT,
    removeoutliers::Bool=REMOVEOUTLIERS_DEFAULT,
    removemissings::Bool=REMOVEMISSINGS_DEFAULT
    )

    return input(
        equation_str_to_strarr!(equation),
        data,
        datanames=datanames,
        intercept=intercept,
        method=method,
        panel=panel,
        time=time,
        seasonaladjustment=seasonaladjustment,
        removeoutliers=removeoutliers,
        removemissings=removemissings
    )
end

function input(
    equation::Array{String},
    data::Union{Array{Float64}, Array{Float32}, Array{Union{Float32, Missing}}, Array{Union{Float64, Missing}}, Tuple, DataFrame, Nothing};
    datanames::Union{Array, Array{Symbol, 1}, Nothing}=nothing,
    method::Union{Symbol, String}=METHOD_DEFAULT,
    intercept::Bool=INTERCEPT_DEFAULT,
    panel::Union{Symbol, String, Nothing}=PANEL_DEFAULT,
    time::Union{Symbol, String, Nothing}=TIME_DEFAULT,
    seasonaladjustment::Union{Dict, Array, Nothing}=SEASONALADJUSTMENT_DEFAULT,
    removeoutliers::Bool=REMOVEOUTLIERS_DEFAULT,
    removemissings::Bool=REMOVEMISSINGS_DEFAULT
    )

    datanames = get_datanames_from_data(data, datanames)

    if datanames == nothing
        error(DATANAMES_REQUIRED)
    end

    equation = strarr_to_symarr!(equation_converts_wildcards!(equation, datanames))

    if isempty(equation)
        error(VARIABLES_NOT_DEFINED)
    end

    return input(
        equation,
        data,
        datanames=datanames,
        intercept=intercept,
        method=method,
        panel=panel,
        time=time,
        seasonaladjustment=seasonaladjustment,
        removeoutliers=removeoutliers,
        removemissings=removemissings
    )
end

function input(
    equation::Array{Symbol},
    data::Union{Array{Float64}, Array{Float32}, Array{Union{Float32, Missing}}, Array{Union{Float64, Missing}}, Tuple, DataFrame, Nothing};
    datanames::Union{Array, Array{Symbol, 1}, Nothing}=nothing,
    method::Union{Symbol, String}=METHOD_DEFAULT,
    intercept::Bool=INTERCEPT_DEFAULT,
    panel::Union{Symbol, String, Nothing}=PANEL_DEFAULT,
    time::Union{Symbol, String, Nothing}=TIME_DEFAULT,
    seasonaladjustment::Union{Dict, Array, Nothing}=SEASONALADJUSTMENT_DEFAULT,
    removeoutliers::Bool=REMOVEOUTLIERS_DEFAULT,
    removemissings::Bool=REMOVEMISSINGS_DEFAULT
    )

    if length(equation) != size(equation, 1)
        equation = convert(Vector{Symbol}, equation[1,:])
    end

    method = Symbol(lowercase(string(method)))

    if method == :precise
        datatype = Float64
    elseif method == :fast
        datatype = Float32
    else
        error(INVALID_METHOD)
    end

    if datanames == nothing
        datanames = get_datanames_from_data(data, datanames)
    end

    if length(datanames) != size(datanames, 1)
        datanames = convert(Vector{Symbol}, datanames[1,:])
    end

    data = get_data_from_data(data)

    if !ModelSelection.in_vector(equation, datanames)
        error(SELECTED_VARIABLES_DOES_NOT_EXISTS)
    end

    if !isa(data, Array{Union{Missing, datatype}}) || !isa(data, Array{Union{datatype}})
        data = convert(Matrix{Union{Missing, datatype}}, data)
    end

    if time != nothing
        if isa(time, String)
            time = Symbol(time)
        end
        if ModelSelection.get_column_index(time, datanames) == nothing
            error(TIME_VARIABLE_INEXISTENT)
        end

    end

    if panel != nothing
        if isa(panel, String)
            panel = Symbol(panel)
        end
        if ModelSelection.get_column_index(panel, datanames) == nothing
            error(PANEL_VARIABLE_INEXISTENT)
        end
    end

    modelselection_data, method, seasonaladjustment, removeoutliers = execute(
        equation,
        data,
        datanames,
        method,
        intercept;
        panel=panel,
        time=time,
        seasonaladjustment=seasonaladjustment,
        removeoutliers=removeoutliers,
        removemissings=removemissings
    )

    modelselection_data = addextras(modelselection_data, method, seasonaladjustment, removeoutliers)

    return modelselection_data
end

function execute(
    equation::Array{Symbol},
    data::Union{Array{Float64}, Array{Float32}, Array{Union{Float32, Missing}}, Array{Union{Float64, Missing}}},
    datanames::Array{Symbol},
    method::Symbol,
    intercept::Bool;
    panel::Union{Symbol, Nothing}=PANEL_DEFAULT,
    time::Union{Symbol, Nothing}=TIME_DEFAULT,
    seasonaladjustment::Union{Dict, Array, Nothing}=SEASONALADJUSTMENT_DEFAULT,
    removeoutliers::Bool=REMOVEOUTLIERS_DEFAULT,
    removemissings::Bool=REMOVEMISSINGS_DEFAULT
    )

    datatype = method == :precise ? Float64 : Float32
    temp_equation = equation

    if panel != nothing && ModelSelection.get_column_index(panel, temp_equation) == nothing
        temp_equation = vcat(temp_equation, panel)
    end

    if time != nothing && ModelSelection.get_column_index(time, temp_equation) == nothing
        temp_equation = vcat(temp_equation, time)
    end

    (data, datanames) = filter_data_by_selected_columns(data, temp_equation, datanames)
    data = sort_data(data, datanames, panel=panel, time=time)

    panel_data = nothing
    if panel != nothing
        if validate_panel(data, datanames, panel=panel)
            panel_data = data[:, ModelSelection.get_column_index(panel, datanames)]
        else
            error(PANEL_ERROR)
        end
    end

    time_data = nothing
    if time != nothing
        if validate_time(data, datanames, panel=panel, time=time)
            time_data = data[:, ModelSelection.get_column_index(time, datanames)]
        else
            error(TIME_ERROR)
        end
    end

    if seasonaladjustment != nothing && time != nothing
        seasonal_adjustments(data, seasonaladjustment, datanames)
    elseif seasonaladjustment != nothing && time == nothing
        error(TIME_VARIABLE_INEXISTENT)
    end

    (data, datanames) = filter_data_by_selected_columns(data, equation, datanames)

    depvar = equation[1]
    expvars = equation[2:end]

    nobs = size(data, 1)

    if intercept
        data = hcat(data, ones(nobs))
        push!(expvars, :_cons)
        push!(datanames, :_cons)
    end

    if removeoutliers
        remove_outliers(data)
    end

    depvar_data = data[1:end, 1]
    expvars_data = data[1:end, 2:end]

    if removemissings
        depvar_data, expvars_data, panel_data, time_data = ModelSelection.filter_raw_data_by_empty_values(datatype, depvar_data, expvars_data, panel_data, time_data)
    end

    depvar_data, expvars_data, panel_data, time_data = ModelSelection.convert_raw_data(datatype, depvar_data, expvars_data, panel_data, time_data)

    nobs = size(depvar_data, 1)

    modelselection_data = ModelSelection.ModelSelectionData(
        equation,
        depvar,
        expvars,
        panel,
        time,
        depvar_data,
        expvars_data,
        panel_data,
        time_data,
        intercept,
        datatype,
        removemissings,
        nobs
    )

    return modelselection_data, method, seasonaladjustment, removeoutliers
end


max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
BE = BEAST

T = Float64
tol = eps(T) * 10^3

# Test on cubic polynomial in two variables
c = T[7,-12,3,2,-1,10,-3,2,8,-7];
p = BE.Poly{2,T}(c);
x = [2.0, -3.0];
y = BE.polyval(p,x);
u = x[1]; v = x[2];
y2 = 7 - 12u + 3v + 2(u^2) -u*v +10(v^2) -3(u^3) +2(u^2*v) + 8(u*v^2) -7(v^3)

@test abs(y2-y) < tol


text:
#=----------------------------------------------------------------------------

Image Projective Geometry

Functions supporting projective geometry for computer vision.

Copyright (c) 2016 Peter Kovesi
pk@peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

PK February  2016
   September 2018 - Updates for v0.7/v1.0

----------------------------------------------------------------------------=#
"""
**ImageProjectiveGeometry**

Functions supporting projective geometry for computer vision.

Peter Kovesi

[peterkovesi.com](http://peterkovesi.com)

*Projective*

* Camera - Structure defining parameters of a camera.
* cameraproject - Projects 3D points into camera image.
* imagept2plane - Project image points to a plane and return their 3D locations.
* imagecorners2plane - Get the positions of image corners projected onto a plane.
* imagept2ray - Compute viewing ray corresponding to an image point.
* imagept2ray! - Compute viewing ray corresponding to an image point.
* camera2projmatrix - Generate a camera projection matrix from a Camera structure.
* decomposecamera -  Decomposition of a camera projection matrix.
* rq3 - RQ decomposition of 3x3 matrix.
* makehomogeneous - Appends a scale of 1 to an array inhomogeneous coordinates.
* makeinhomogeneous - Converts homogeneous coords to inhomogeneous coordinates.
* hnormalise - Normalises array of homogeneous coordinates to a scale of 1.
* hnormalise! - In-place normalisation of homogeneous coordinates to a scale of 1.
* homography1d - Computes 1D homography.
* homography2d - Computes 2D homography.
* solveaffine - Solve affine transformation between two sets of 2D points.
* normalise1dpts - Normalises 1D homogeneous points.
* normalise2dpts - Normalises 2D homogeneous points.
* skew - Constructs 3x3 skew-symmetric matrix from 3-vector.
* hcross - Homogeneous cross product, result normalised to s = 1.
* fundmatrix - Computes fundamental matrix from 8 or more points.
* affinefundmatrix - Computes affine fundamental matrix from 4 or more points.
* fundfromcameras - Fundamental matrix from camera matrices or structures.
* stereorectify - Rectify a stereo pair of images.
* stereorectifytransforms - Compute homographies that transform an image pair into a stereorectified pair.
* transformedimagebounds - Find where the corners of an image are transformed to by transform H and return the bounds.
* imgtrans - Homogeneous transformation of an image.
* imgtrans! - Homogeneous transformation of an image.
* idealimagepts - Ideal image points with no distortion.
* solvestereopt - Homogeneous linear solution of a stereo point.
* undistortimage - Removes lens distortion from an image.
* hline - Plot a 2D line defined in homogeneous coordinates.
* mapimage2plane! - Projects an image onto a plane in 3D.
* plotcamera - Plots graphical representation of camera(s) showing pose.

*Corner Features*

* shi_tomasi - Shi-Tomasi corner detector.
* harris - Harris corner detector.
* noble - Noble's variant of the Harris corner detector.
* coherence - Compute image coherence from structure tensor.
* structuretensor - Compute structure tensor values over an image.
* hessianfeatures  - Computes determiant of hessian features in an image.
* fastradial - Loy and Zelinski's fast radial feature detector.

*RANSAC*

* ransac - Robustly fits a model to data with the RANSAC algorithm.
* ransacfithomography - Fits 2D homography using RANSAC.
* ransacfitfundmatrix - Fits fundamental matrix using RANSAC.
* ransacfitaffinefundmatrix - Fits affine fundamental matrix using RANSAC.
* ransacfitplane - Fits plane to 3D array of points using RANSAC.
* ransacfitline - Fits line to 3D array of points using RANSAC.
* iscolinear - Are 3 points colinear.
* fitline2d - Least squares fit of a line to a set of 2D points.
* fitline3d - Fits a line to a set of 3D points.
* fitplane - Solves coefficients of plane fitted to 3 or more points.

*RANSAC demos*

* fitlinedemo - Demonstrates RANSAC line fitting.
* fitplanedemo - Demonstrates RANSAC plane fitting.
* fitfunddemo - Example of fundamental matrix computation.
* fithomogdemo - Example of finding a homography.

*Transforms*

* trans - Homogeneous transformation for a translation by x, y, z.
* rotx - Homogeneous transformation for a rotation about the x axis.
* roty - Homogeneous transformation for a rotation about the y axis.
* rotz - Homogeneous transformation for a rotation about the z axis.
* dhtrans - Computes Denavit Hartenberg matrix.
* homotrans - Homogeneous transformation of points/lines.
* invht - Inverse of a homogeneous transformation matrix.
* inveuler - Inverse of Euler transform.
* invrpy - Inverse of Roll Pitch Yaw transform.
* angleaxis - Constructs angle-axis descriptor.
* normaliseangleaxis - Normalises angle-axis descriptor.
* angleaxisrotate - Uses angle axis descriptor to rotate vectors.
* angleaxis2matrix - Converts angle-axis descriptor to 4x4 homogeneous transformation  matrix.
* matrix2angleandaxis - Decompose homogeneous matrix to angle and axis.
* matrix2angleaxis - Homogeneous matrix to angle-axis description.
* quaternion - Construct quaternion.
* quaternionconjugate - Conjugate of a quaternion.
* quaternionproduct - Computes product of two quaternions.
* quaternionrotate - Rotates a 3D vector by a quaternion.
* quaternion2matrix - Quaternion to a 4x4 homogeneous transformation matrix.
* matrix2quaternion - Homogeneous matrix to quaternion.
* vector2quaternion - Embeds 3-vector in a quaternion representation.

*Geometry*

* ray2raydist -  Minimum distance between two 3D rays.
* circleintersectionpts - Finds intersections of two circles.
* rectintersect - Test if rectangles intersect.
* pointinconvexpoly - Determine if a 2D point is within a convex polygon.
* convexpolyintersect - Test if convex polygons intersect.
* isconvexpoly - Test if 2D polygon is convex.

*Utilities*

* nonmaxsuppts - Non-maximal suppression for features/corners.
* derivative3 - 3-Tap discrete derivative filters.
* derivative5 - 5-Tap 1st and 2nd discrete derivatives.
* derivative7 - 7-Tap 1st and 2nd discrete derivatives.
* gaussfilt -  Small wrapper function for convenient Gaussian filtering.
* dilate1d - 1D morphological dilation of a signal.
* erode1d - 1D morphological erosion of a signal.
* imdilate - Image morpholgical dilation.
* imerode - Image morpholgical erosion.
* circularstruct - Generate circular structuring element for morphological operations.
* medfilt2 - Convenience wrapper for median filtering.
* stdfilt2 - Compute local standard deviation across an image.
* floatyx - Convert 2D AbstractImage to 2D float array with y x spatial order.
* histtruncate - Truncates ends of an image histogram.
* imgnormalise/imgnormalize - Normalises image values to 0-1, or to desired mean and variance.
* matchbycorrelation - Match image feature points by correlation.
* grey2census - Convert image grey scale values to census values.
* briefcoords - Compute BRIEF descriptor sampling coordinates within a patch.
* grey2lbp - Convert image grey scale values to local binary pattern.
* grey2lbp! - Convert image grey scale values to local binary pattern.
* keypause - Wait for user to hit return before continuing.

"""
module ImageProjectiveGeometry

include("projective.jl")
include("transforms3d.jl")
include("ransac.jl")
include("cornerfeatures.jl")
include("utilities.jl")
include("ransacdemo.jl")
include("geometry.jl")

end  # module


text:
sbmlfiles = [
    # a test model from BIGG
    (
        "e_coli_core.xml",
        "http://bigg.ucsd.edu/static/models/e_coli_core.xml",
        "b4db506aeed0e434c1f5f1fdd35feda0dfe5d82badcfda0e9d1342335ab31116",
        72,
        95,
    ),
    # a relatively new non-curated model from biomodels
    (
        "T1M1133.xml",
        "https://www.ebi.ac.uk/biomodels/model/download/MODEL1909260004.4?filename=T1M1133.xml",
        "2b1e615558b6190c649d71052ac9e0dc1635e3ad281e541bc7d4fdf2892a5967",
        2517,
        3956,
    ),
    # a curated model from biomodels
    (
        "Dasgupta2020.xml",
        "https://www.ebi.ac.uk/biomodels/model/download/BIOMD0000000973.3?filename=Dasgupta2020.xml",
        "958b131d4df2f215dae68255433542f228601db0326d26a54efd08ddcf823489",
        2,
        6,
    ),
]

@testset "Loading of models from various sources" begin
    for (sbmlfile, url, hash, expected_mets, expected_rxns) in sbmlfiles
        if !isfile(sbmlfile)
            download(url, sbmlfile)
        end

        cksum = bytes2hex(sha256(open(sbmlfile)))
        if cksum != hash
            @warn "The downloaded model `$sbmlfile' seems to be different from the expected one. Tests will likely fail." cksum
        end

        @testset "Loading of $sbmlfile" begin
            mdl = readSBML(sbmlfile)

            @test typeof(mdl) == Model

            mets, rxns, _ = getS(mdl)

            @test length(mets) == expected_mets
            @test length(rxns) == expected_rxns
        end
    end
end


text:
module ArrayInterfaceGPUArrays

using Adapt
using ArrayInterfaceCore
using GPUArraysCore

ArrayInterfaceCore.fast_scalar_indexing(::Type{<:GPUArraysCore.AbstractGPUArray}) = false
@inline ArrayInterfaceCore.allowed_getindex(x::GPUArraysCore.AbstractGPUArray, i...) = GPUArraysCore.@allowscalar(x[i...])
@inline ArrayInterfaceCore.allowed_setindex!(x::GPUArraysCore.AbstractGPUArray, v, i...) = (GPUArraysCore.@allowscalar(x[i...] = v))

function Base.setindex(x::GPUArraysCore.AbstractGPUArray, v, i::Int)
    _x = copy(x)
    ArrayInterfaceCore.allowed_setindex!(_x, v, i)
    return _x
end

function ArrayInterfaceCore.restructure(x::GPUArraysCore.AbstractGPUArray, y)
    reshape(Adapt.adapt(ArrayInterfaceCore.parameterless_type(x), y), Base.size(x)...)
end

function ArrayInterfaceCore.lu_instance(A::GPUArraysCore.AbstractGPUMatrix{T}) where {T}
    qr(similar(A, 1, 1))
end

# Doesn't do much, but makes a gigantic change to the dependency chain.
# ArrayInterface.device(::Type{<:GPUArraysCore.AbstractGPUArray}) = ArrayInterface.GPU()

end


text:
function constrain_neigh_plus!(prob, uni, N, xsol, nsol, T0::Bool)
	M = prob;
	M .= M .+ uni
	zeroprob = (T0 == true) ? -Inf : 0.0
	if xsol == 1
		for n in 0:N+1
			M[0,n] = (n == nsol || n == N + 1) ? M[0,n]  : zeroprob
			M[1,n] = (n > nsol) ? M[1,n] : zeroprob
		end
	else
		for n in 0:N+1
			M[0,n] = (n == nsol) ? M[0,n]  : zeroprob
			M[1,n] = (n > nsol) ? M[1,n] : zeroprob
		end
	end
	M[1,0] = zeroprob
	M[1,N+1] = zeroprob
	if T0
		M .= M .- maximum(M)
	else
		M .= M./sum(M)
	end
	return M
end

function constrain_neigh_minus!(prob, uni, N, xsol, nsol, T0::Bool)

	M = prob;
	M .= M .+ uni
	zeroprob = (T0 == true) ? -Inf : 0.0
	if xsol == 0
	    if nsol == 0
	       for n = 1:N+1
		  M[0,n] = zeroprob
	       end
	       for n = 0:N+1
		  M[1,n] = zeroprob
	       end
	    elseif nsol == N + 1
	       for n = 0:N
		  M[0,n] = zeroprob
	       end
	    else
	        M[0,0] = zeroprob
		M[0, N+1] = zeroprob
		for x in 0:1, n in 1:N
		     M[x,n] = (n == nsol) ? M[x,n] : zeroprob
		end
	    end
	else
		for x in 0:1, n in 0:N+1
		     M[x,n] = (n < nsol) ? M[x,n]  : zeroprob
		end
	end
	M[1,0] = zeroprob
	M[1,N+1] = zeroprob
	if T0
		M .= M .- maximum(M)
	else
		M .= M./sum(M);
	end
	return M
end

function decimate_post(allvar::DCAlign.AllVar, T0::Bool)

	@extract allvar : seq pbf jh alg
	@extract pbf : P L q N
	sol = zeros(Int64,L,2)
	fixed = zeros(Bool,L)
	can_assign = ones(Bool,L)
	Nfree = count(fixed[i] == false for i in 1:L)
	uni = 1/(2*(N+2))
	while Nfree > 0
		Hmax = -Inf
		imax = NaN; xmax = NaN; nmax = NaN;
		if Nfree == L
			for i in 1:L
				for x in 0:1, n in 0:N+1
					if(P[i][x,n] > Hmax)
						Hmax = P[i][x,n]
						xmax = x
						nmax = n
						imax = i
					end
				end
			end
		else
			idx = findall(can_assign .== true)
			if length(idx) > 2
				display(["we have a problem, indices are" length(idx)])
			elseif length(idx) == 2
				if idx[1] < idx[2]
					idxm = idx[1]
					idxM = idx[2]
				else
					idxm = idx[2]
					idxM = idx[1]
				end
				P[idxm] = constrain_neigh_minus!(P[idxm], uni, N, sol[idxm+1,1], sol[idxm+1,2], T0)
				P[idxM] = constrain_neigh_plus!(P[idxM], uni, N, sol[idxM-1,1], sol[idxM-1,2], T0)
			elseif length(idx) == 1
				if idx[1]-1 > 0
				    if fixed[idx[1]-1] == true
					P[idx[1]] = constrain_neigh_plus!(P[idx[1]], uni, N, sol[idx[1]-1,1], sol[idx[1]-1,2],T0)
				    end
				end
				if idx[1]+1 <= L
				    if fixed[idx[1]+1] == true
					P[idx[1]] = constrain_neigh_minus!(P[idx[1]], uni,N, sol[idx[1] + 1,1], sol[idx[1] + 1,2], T0)
				    end
				end
			end
			for j in 1:length(idx)
			    for x in 0:1, n in 0:N+1
			    	if(maximum(P[idx[j]][x,n]) > Hmax)
					Hmax = P[idx[j]][x,n]
					xmax = x
					nmax = n
					imax = idx[j]
				end
			    end
			end
		end
		fixed[imax] = true
		sol[imax,1] = xmax; sol[imax,2] = nmax;
		if Nfree == L
			can_assign = zeros(Bool, L,)
			if imax == 1
				can_assign[2] = true
			elseif imax == L
				can_assign[L-1] = true
			else
				can_assign[imax-1] = true
				can_assign[imax+1] = true
			end
		else
			can_assign[imax] = false
			if imax == 1
				can_assign[2] = (fixed[2] == false) ? true : false
			elseif imax == L
				can_assign[L-1] = (fixed[L-1] == false) ? true : false
			else
				can_assign[imax-1] = (fixed[imax-1] == false) ? true : false
				can_assign[imax+1] = (fixed[imax+1] == false) ? true : false
			end
		end
		Nfree = count(fixed[i] == false for i in 1:L)
	end
	s = ""
	for i in 1:L
		x = sol[i,1]; n = sol[i,2];
		if x == 0
			s *= "-"
		else
			s *= seq.strseq[n]
		end
	end
	return s, P

end


text:
function dft(n, fx, fsample)
    if n <= 0
        throw(DomainError(n, "sample points must be > 0.\r\n"))
    end

    if fsample <= 0
        throw(DomainError(fsample, "sample must be > 0.\r\n"))
    end

    X = zeros(complex(Float64),n)

    for i = 0:n-1
        for j = 0:n-1
            X[i+1] += fx[j+1] * exp(-im * 2 * pi * j * i / n)
        end
    end

    return X
end


text:
# generic constraints

function MOI.get(o::Optimizer, ::MOI.ConstraintName, ci::CI)::String
    return GC.@preserve o unsafe_string(SCIPconsGetName(cons(o, ci)))
end

function MOI.set(o::Optimizer, ::MOI.ConstraintName, ci::CI, name::String)
    @SCIP_CALL SCIPchgConsName(o, cons(o, ci), name)
    return nothing
end

function MOI.set(o::Optimizer, ::MOI.ConstraintName, ci::CI{VI}, name::String)
    throw(MOI.VariableIndexConstraintNameError())
    return nothing
end

function MOI.is_valid(o::Optimizer, c::CI{F, S}) where {F, S}
    cons_set = get(o.constypes, (F, S), nothing)
    if cons_set === nothing
        return false
    end
    if !in(ConsRef(c.value), cons_set)
        return false
    end
    return haskey(o.inner.conss, SCIP.ConsRef(c.value))
end

function MOI.get(o::Optimizer, ::Type{MOI.ConstraintIndex}, name::String)
    ptr = SCIPfindCons(o, name)
    if ptr == C_NULL
        return nothing
    end
    cref = get(o.reference, ptr, nothing)
    if cref === nothing
        return cref
    end
    for ((F, S), setref) in o.constypes
        if cref in setref
            return CI{F,S}(cref.val)
        end
    end
    error("Constraint type not found for constraint $cref")
end


text:
export MCObservableSet
typealias MCObservableSet{Obs<:MCObservable} Dict{ASCIIString, Obs}

export makeMCObservable!

function makeMCObservable!{Obs}(oset::MCObservableSet{Obs}, name::ASCIIString)
  if haskey(oset, name)
    warn("""Observable "$name" already exists. (Skipped)""")
  else
    oset[name] = Obs()
  end
end

function reset!(oset::MCObservableSet)
  for v in values(oset)
    reset!(oset)
  end
end

function show(io::IO, obs::MCObservableSet, sorted::Bool = true)
  if sorted
    ks = sort([k for k in keys(obs)])
    for k in ks
      println(io, k, " : ", obs[k])
    end
  else
    for k in keys(obs)
      println(io, k, " : ", obs[k])
    end
  end
end


text:
"""
Step time for a TEM grid (2D)
"""
function step_time(F::EHTuple{T,2}, Φ::PML, C::Coefficients, TEsource, TMsource, sourceindex) where T # F is the EM field, C is a wrapper for the Coefficients used here
    E, H = F

    H = update_H(E, H, C)
    E = update_E(E, H, C)

    ΦEaux = updatePML(Φ.E, H, C, 1, Φ.N)
    ΦHaux = updatePML(Φ.H, E, C, 2, Φ.N)

    H = applyPML(H, ΦHaux, C, 1, Φ.N)
    E = applyPML(E, ΦEaux, C, 2, Φ.N)

    Φ = PML(ΦEaux,ΦHaux,Φ.N)

    E[3][sourceindex...] = E[3][sourceindex...] .+ TEsource #hard source acts like antenna

    H[3][sourceindex...] = H[3][sourceindex...] .+ TMsource #hard source acts like antenna

    return ((E, H), Φ)
end

"""
Step time for a full 3D grid
"""
function step_time(F::EHTuple{T}, Φ::PML, C::Coefficients, Esourcenow, Hsourcenow, sourceindex) where T # F is the EM field, C is a wrapper for the Coefficients used here
    E, H = F

    H = update_H(E, H, C)
    E = update_E(E, H, C)

    ΦEaux = updatePML(Φ.E, H, C, 1, Φ.N)
    ΦHaux = updatePML(Φ.H, E, C, 2, Φ.N)

    H = applyPML(H, ΦHaux, C, 1, Φ.N)
    E = applyPML(E, ΦEaux, C, 2, Φ.N)

    Φ = PML(ΦEaux,ΦHaux,Φ.N)
    E[sourceindex[1]][sourceindex[2]...] = E[sourceindex[1]][sourceindex[2]...] .+ Esourcenow #hard source acts like antenna
    Hnew[sourceindex[1]+ (sourceindex[1]==3 ? -2 : 1)][sourceindex[2]...] = F[2][sourceindex[1]+ (sourceindex[1]==3 ? -2 : 1)][sourceindex[2]...] + Hsourcenow
    return ((E, H), Φ)
end

function FDTD_propagate(space::ProductAxis, time::LinearAxis, f₀::T, nPML::Int = 0;
                        source,
                        medium::Medium,
                        sourceindex = origin(space),
                        detectorindex =  (:, [:, :, :]),
                        reduction = identity
                        ) where T

      Ex = zeros(T, size(space)...)
      Ey = zeros(T, size(space)...)
      Ez = zeros(T, size(space)...)

      Hx = zeros(T, size(space)...)
      Hy = zeros(T, size(space)...)
      Hz = zeros(T, size(space)...)

      F = (E,H) = ((Ex,Ey,Ez), (Hx,Hy,Hz))

   Φ = PML(space, nPML)
   C = Coefficients(space, time, medium, f₀, nPML)

    if detectorindex isa Tuple{Colon, Array{Colon,1}}
        out = Array{EMField{T}}(undef,length(time))
        #out = Array{PML{T}}(undef,length(time))
        EM_out = true
    else
        Sout = size(F[detectorindex[1]][detectorindex[2]...])
        Nout = length(Sout)
        out = zeros(Sout..., length(time))
        colons = fill(:, Nout)
        EM_out = false
    end
    @showprogress "Simulating with FDTD..." for i in time.i
        Fnew, Φnew = step_time(F, Φ, C, source[i], source[i], sourceindex)
        if EM_out
             out[i] =EMField(Fnew)
        else
             out[colons...,i] = Fnew[detectorindex[1]][detectorindex[2]][detectorindex[3]...] # Don't mess with this, it works but I can't remember how
        end
        F = Fnew;
        Φ = Φnew;
    end
    return out
end

function detect(F::EHTuple{T}, detectorindex::CartesianIndices) where T
    out = zero(T)
    for I in detectorindex
        out = out + F[1][I]^2
        out = out + (F[2][I])^2
        out = out + F[3][I]^2
    end
    return out/(2*η₀ * prod(size(detectorindex)))
end

"""
A draft function for generating the radar cross section of a particular medium
"""
function radar_propagate(space::ProductAxis, time::LinearAxis, f₀::T, nPML::Int = 0;
                        source,
                        medium::Medium,
                        sourceindex = origin(space),
                        detectorindex::CartesianIndices = CartesianIndices(sourceindex[2])
                        ) where T
   Ex = zeros(T, size(space)...)
   Ey = zeros(T, size(space)...)
   Ez = zeros(T, size(space)...)

   Hx = zeros(T, size(space)...)
   Hy = zeros(T, size(space)...)
   Hz = zeros(T, size(space)...)

   F = (E,H) = ((Ex,Ey,Ez), (Hx,Hy,Hz))
   Φ = PML(space, nPML)
   C = Coefficients(space, time, medium, f₀, nPML)
   out = Vector{T}(undef, time.N)
   for i in 1:time.N
        Fnew, Φnew = step_time(F, Φ, C, source[i], source[i], sourceindex) #You actually need Fnew here to get the compiler to not tidy the loop away
        @inbounds out[i] = sum([sum(view(Fnew[1][k], detectorindex).^2)./(2*η₀ * prod(size(detectorindex))) for k in 1:3])
        F = Fnew ; Φ = Φnew
    end
    return out
end


text:
@defcomp TotalAbatementCosts begin
    region = Index()

    tc_totalcosts_co2 = Parameter(index=[time, region], unit="\$million")
    tc_totalcosts_ch4 = Parameter(index=[time, region], unit="\$million")
    tc_totalcosts_n2o = Parameter(index=[time, region], unit="\$million")
    tc_totalcosts_linear = Parameter(index=[time, region], unit="\$million")
    pop_population = Parameter(index=[time, region], unit="million person")

    tct_totalcosts = Variable(index=[time,region], unit="\$million")
    tct_per_cap_totalcostspercap = Variable(index=[time,region], unit="\$/person")

    function run_timestep(p, v, d, t)

        for r in d.region
            v.tct_totalcosts[t,r] = p.tc_totalcosts_co2[t,r] + p.tc_totalcosts_n2o[t,r] + p.tc_totalcosts_ch4[t,r] + p.tc_totalcosts_linear[t,r]
            v.tct_per_cap_totalcostspercap[t,r] = v.tct_totalcosts[t,r] / p.pop_population[t,r]
        end
    end
end


text:
using Random, Distributions
using ProgressBars
using Plots

ENV["JULIA_DEBUG"] = "all"

const t = 3 # number of frames wrong
const expectedrates = [
    0.585, # delay-critical user u
    2, # delay-tolerant user v
]
const ζ = 2 # path loss exponent
const a = [
    0.7, # power for user u
    0.3,  # power for user v
]
const d = [
    2, # distance of user u
    1,  # distance of user v
]
const total_frames = 10000

channel_distribution = Normal(0, √2 / 2)
plot_points = ([], [])

Threads.@threads for SNR ∈ tqdm(0:5:45)
    for τ ∈ [0.6, 1]
        successful_frames = zeros((2, 1))

        progress_bar = 1:total_frames
        # progress_bar = ProgressBar(1:total_frames)
        # set_description(progress_bar, "($SNR dB, τ = $τ)")
        for frame ∈ progress_bar
            # Step 1: Calculate g of channel for each user
            h = rand(channel_distribution, (2, 2)) * [1, 1im]
            g = h ./ sqrt.(1 .+ d .^ ζ)

            # Step 1.1: Sort channel gains
            if abs2(h[1]) > abs2(h[2])
                h[2], h[1] = h[1], h[2]
                g = h ./ sqrt.(1 .+ d .^ ζ)
            end

            # Step 1.5: Apply time-correlated small fading coefficient
            ω = rand(channel_distribution, (2, 2)) * [1, 1im]
            h = τ^t * h + √(1 - τ^(2t)) * ω
            g = h ./ sqrt.(1 .+ d .^ ζ)

            # Step 2: Calculate AWGN
            σ = 10^(-SNR / 20)
            n = rand(Normal(0, σ), 2)

            # Step 3: Calculate reception SNR
            # Note: We are using the actual noise value n, instead of just the SNR-based σ
            SNRuu = a[1] * abs2(g[1]) / (a[2] * abs2(g[1]) + σ^2)
            SNRvv = a[2] * abs2(g[2]) / σ^2
            SNRvu = a[1] * abs2(g[2]) / (a[2] * abs2(g[2]) + σ^2)

            receiver_rates = [log2(1 + SNRuu), log2(1 + SNRvv), log2(1 + SNRvu)]

            successful_frames += [
                receiver_rates[1] >= expectedrates[1],
                receiver_rates[2] >= expectedrates[2] &&
                receiver_rates[3] >= expectedrates[1],
            ]
        end

        outage_rates = @. 1 - successful_frames / total_frames

        for rate ∈ outage_rates
            if rate > 0
                push!(plot_points[1], SNR)
                push!(plot_points[2], rate)
            end
        end
    end
end

## THEORETICAL PLOTS

# Calculate theoretical values
λ = 1 ./ (1 .+ d .^ ζ)
r = 2 .^ expectedrates .- 1
φ = max(r[1] / (a[1] - a[2] * r[1]), r[2] / a[2])

SNR = LinRange(0, 45, 100)
σ = @. 10^(-SNR / 20)
ρ = @. 1 / σ^2

plotly()
plot()
for τ ∈ [0.6, 1]
    # NOMA - Outdated CSI
    PoutTheoreticalNOMA = @. [
        1 - exp(-2 * r[1] * σ^2 / (2 - τ^(2t)) / (a[1] - a[2] * r[1]) / λ[1]),
        1 - 2 * exp(-φ * σ^2 / λ[2]) + exp(-2 * φ * σ^2 / (2 - τ^(2t)) / λ[2]),
    ]
    PoutTheoreticalOMA = @. [
        1 - exp(
            -2 * (2^(2 * expectedrates[1]) - 1) / ((2 - τ^(2t)) * ρ * λ[1]),
        ),
        1 - 2 * exp(-(2^(2 * expectedrates[2]) - 1) / ρ / λ[2]) + exp(
            -2 * (2^(2 * expectedrates[2]) - 1) / ((2 - τ^(2τ)) * ρ * λ[2]),
        ),
    ]

    label = [
        "a"
        "b"
    ]

    plot!(
        SNR,
        PoutTheoreticalNOMA,
        yscale = :log10,
        label = ["NOMA, User u, τ = $τ" "NOMA, User v, τ = $τ"],
        lw = 2
    )
    plot!(
        SNR,
        PoutTheoreticalOMA,
        yscale = :log10,
        label = ["OMA, User u, τ = $τ" "NOMA, User v, τ = $τ"],
        lw = 2
    )
end

plot!(
    plot_points[1],
    plot_points[2],
    seriestype = :scatter,
    yscale = :log10,
    label = "NOMA Simulation results",
)

xlabel!("SNR (dB)")
ylabel!("Outage rate")

display(plot!(legend = :best, size = (1000,800)))


text:
using Base.Test
using DataFrames
import WCS

const rcf = RunCamcolField(3900, 6, 269)

"""
Get the PSF located at a particular world location in an image.

Args:
 - world_loc: A location in world coordinates.
 - img: An Image

Returns:
 - An array of PsfComponent objects that represents the PSF as a mixture
     of Gaussians.
"""
function get_source_psf(world_loc::Vector{Float64}, img::Image, psf_K::Int)
    # Some stamps or simulated data have no raw psf information. In that case,
    # just use the psf from the image.
    if size(img.raw_psf_comp.rrows) == (0, 0)
        # Also return a vector of empty psf params
        return img.psf, fill(fill(NaN, length(PsfParams)), psf_K)
    else
        pixel_loc = WCS.world_to_pix(img.wcs, world_loc)
        psfstamp = Model.eval_psf(img.raw_psf_comp, pixel_loc[1], pixel_loc[2])
        return PSF.fit_raw_psf_for_celeste(psfstamp, psf_K)
    end
end

function test_images()
    # A lot of tests are in a single function to avoid having to reload
    # the full image multiple times.
    images = SDSSIO.load_field_images(rcf, datadir)

    dir = "$datadir/$(rcf.run)/$(rcf.camcol)/$(rcf.field)"
    fname = @sprintf "%s/photoObj-%06d-%d-%04d.fits" dir rcf.run rcf.camcol rcf.field
    cat_entries = SDSSIO.read_photoobj_celeste(fname)

    ea = make_elbo_args(images, cat_entries, patch_radius_pix=1e-6)

    # Just check some basic facts about the catalog.
    @test length(cat_entries) == 805
    @test 0 < sum([ce.is_star for ce in cat_entries]) < 805

    # Find an object near the middle of the image.
    ctr = WCS.pix_to_world(images[3].wcs, [images[3].H / 2, images[3].W / 2])
    dist = [sqrt((ce.pos[1] - ctr[1])^2 + (ce.pos[2] - ctr[2])^2)
                    for ce in cat_entries]
    obj_index = findmin(dist)[2]    # index of closest object
    obj_loc = cat_entries[obj_index].pos    # location of closest object

    # Test get_source_psf at point while we have the images loaded.
    test_b = 3
    img = ea.images[test_b]
    ea_obj = make_elbo_args(images, cat_entries[obj_index:obj_index])
    pixel_loc = WCS.world_to_pix(img.wcs, obj_loc)
    original_psf_val = Model.eval_psf(img.raw_psf_comp, pixel_loc[1], pixel_loc[2])

    original_psf_celeste =
        PSF.fit_raw_psf_for_celeste(original_psf_val, ea.psf_K)[1]
    fit_original_psf_val = PSF.get_psf_at_point(original_psf_celeste)

    obj_psf = get_source_psf(ea_obj.vp[1][ids.u], img, ea.psf_K)[1]
    obj_psf_val = PSF.get_psf_at_point(obj_psf)

    # The fits should match exactly.
    @test_approx_eq_eps(obj_psf_val, fit_original_psf_val, 1e-6)

    # The raw psf will not be as good.
    @test_approx_eq_eps(obj_psf_val, original_psf_val, 1e-2)

    cat_several = [cat_entries[1]; cat_entries[obj_index]]
    ea_several = make_elbo_args(ea_obj.images, cat_several)

    # The second set of vp is the object of interest
    point_patch_psf = PSF.get_psf_at_point(ea_several.patches[2, test_b].psf)
    # The threshold for the test below was formerly 1e-6.
    # Is it a problem I needed to increase it?
    @test_approx_eq_eps(obj_psf_val, point_patch_psf, 1e-4)
end

function test_stamp_get_object_psf()
    stamp_blob, stamp_mp, body = gen_sample_star_dataset()
    img = stamp_blob[3]
    obj_index =    stamp_mp.vp[1][ids.u]
    pixel_loc = WCS.world_to_pix(img.wcs, obj_index)
    original_psf_val = PSF.get_psf_at_point(img.psf)

    obj_psf_val = PSF.get_psf_at_point(
      get_source_psf(stamp_mp.vp[1][ids.u], img, 2)[1])
    @test_approx_eq_eps(obj_psf_val, original_psf_val, 1e-6)
end

function test_set_patch_size()
    # Test that the patch size gets most of the light from a variety of
    # galaxy shapes.
    # This shows that the current patch size is actually far too conservative.

    function gal_catalog_from_scale(gal_scale::Float64, flux_scale::Float64)
        CatalogEntry[CatalogEntry(world_location, false,
                                   flux_scale * fluxes, flux_scale * fluxes,
                                   0.1, .01, pi/4, gal_scale, "sample", 0) ]
    end

    srand(1)
    images0 = SampleData.load_stamp_blob(datadir, "164.4311-39.0359_2kpsf")
    img_size = 150
    for b in 1:5
            images0[b].H, images0[b].W = img_size, img_size
    end
    fluxes = [4.451805E+03,1.491065E+03,2.264545E+03,2.027004E+03,1.846822E+04]

    world_location = WCS.pix_to_world(images0[3].wcs,
                                      Float64[img_size / 2, img_size / 2])

    for gal_scale in [1.0, 10.0], flux_scale in [0.1, 10.0]
        cat = gal_catalog_from_scale(gal_scale, flux_scale)
        images = Synthetic.gen_blob(images0, cat)
        ea = make_elbo_args(images, cat)

        for b=1:length(images)
            @assert size(ea.images[b].tiles) == (1, 1)
            tile_image = DeterministicVI.tile_predicted_image(
                ea.images[b].tiles[1,1], ea, ea.tile_source_map[b][1,1])

            pixel_center = WCS.world_to_pix(images[b].wcs, cat[1].pos)
            radius = Model.choose_patch_radius(
                pixel_center, cat[1], images[b].psf, ea.images[b])

            circle_pts = fill(false, images[b].H, images[b].W)
            in_circle = 0.0
            for x=1:size(tile_image)[1], y=1:size(tile_image)[2]
                if ((x - pixel_center[1]) ^ 2 + (y - pixel_center[2]) ^ 2) < radius ^ 2
                    in_circle += tile_image[x, y]
                    circle_pts[x, y] = true
                end
            end
            @test in_circle / sum(tile_image) > 0.95
        end
    end
end

test_images()
#test_stamp_get_object_psf()
#test_set_patch_size()


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 1c8d2d00-b7d9-11eb-35c4-47f2a2aa1593
begin
    import Pkg
	ENV["JULIA_MARGO_LOAD_PYPLOT"] = "no thank you"
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="ClimateMARGO", version=v"0.3.2"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="HypertextLiteral", version="0.9"),
		Pkg.PackageSpec(name="Underscores", version="2"),
    ])

	using Plots
	using Plots.Colors
	using ClimateMARGO
	using ClimateMARGO.Models
	using ClimateMARGO.Optimization
	using ClimateMARGO.Diagnostics
	using PlutoUI
    using HypertextLiteral
	using Underscores

	Plots.default(linewidth=5)
end;

# ╔═╡ 9a48a08e-7281-473c-8afc-7ad3e0771269
TableOfContents()

# ╔═╡ 94415ff2-32a2-4b0f-9911-3b93e202f548
const initial_1 = Dict("M" => [2090, 6]);

# ╔═╡ e810a90f-f964-4d7d-acdb-fc3a159dc12e
const initial_2 = Dict("M" => [2080, .7]);

# ╔═╡ a3422533-2b78-4bc2-92bd-737da3c8982d
const initial_3 = Dict("M" => [2080, .7]);

# ╔═╡ bb66d347-99be-4a95-8ba8-57dc9d33384b
const initial_4 = Dict(
	"M" => [2080, 0.7],
	"R" => [2120, 0.2],
);

# ╔═╡ 0a126a12-7b89-4f22-ad66-b64965cd1d49
HTML("""

When editing this notebook inside Pluto, use the  icon next to cells to view their code.

""")

# ╔═╡ 3094a9eb-074d-46c3-9c1e-0a9c94c6ad43
blob(el, color = "red") = @htl("""
$(el)
""")

# ╔═╡ 0b31eac2-8efd-47cd-9571-a2053846343b
function infeasablewarning(x)
	@htl("""
<margo-infeasible>
	$(x)
	<margo-infeasible-label>No solution found</margo-infeasible-label>
</margo-infeasible>

	""")
end

# ╔═╡ ca104939-a6ca-4e70-a47a-1eb3c32db18f
status_ok(x) = x ∈ [
    "OPTIMAL",
    "LOCALLY_SOLVED",
    "ALMOST_OPTIMAL",
    "ALMOST_LOCALLY_SOLVED"
  ]

# ╔═╡ cf90139c-13d8-42a7-aba3-8c431e7854b8
feasibility_overlay(x) = status_ok(x.status) ? as_html : infeasablewarning

# ╔═╡ bd2bfa3c-a42e-4975-a543-84541f66b1c1
begin
	hidecloack(name) = HTML("""

	""")

	"A trick to hide a cell without creating a variable dependency, to make it simpler for PlutoSliderServer.jl."
	cloak(name) = x -> @htl("<plj-cloack class=$(name)>$(x)</plj-cloak>")
end

# ╔═╡ b81de514-2506-4243-8235-0b54dd4a7ec9
colors = (
	baseline=colorant"#dddddd",
	baseline_emissions=colorant"#dddddd",
	baseline_concentrations=colorant"#dddddd",
	baseline_temperature=colorant"#dddddd",
	baseline_damages=colorant"#dddddd",
	temperature=colorant"#edc949",
	above_paris=colorant"#e1575910",
	M=colorant"#4E79A7",
	R=colorant"#F28E2C",
	A=colorant"#59A14F",
	G=colorant"#E15759",
	T_max=colorant"#00000080",
	controls=colorant"#af7aa1",
	damages=colorant"#e15759",
	avoided_damages=colorant"#e49734",
	benefits=colorant"#7abf5e",
	emissions=colorant"brown",
	emissions_1=colorant"#4E79A7",
	concentrations=colorant"brown",
)

# ╔═╡ 73e01bd8-f56b-4bb5-a9a2-85ad223c9e9b
names = (
	baseline="Baseline",
	baseline_emissions="Baseline",
	baseline_concentrations="Baseline",
	baseline_temperature="Baseline",
	baseline_damages="Baseline",
	temperature="Temperature",
	above_paris="Above Paris",
	M="Mitigation",
	R="Removal",
	A="Adaptation",
	G="Geo-engineering",
	T_max="Goal temperature",
	controls="Controls",
	damages="Damages",
	avoided_damages="Avoided Damages",
	benefits="Benefits",
	emissions="Emissions",
	emissions_1="Emissions",
	concentrations="Concentrations",
)

# ╔═╡ ae92ba1f-5175-4704-8240-2de8432df752
@assert keys(colors) == keys(names)

# ╔═╡ 8ac04d55-9034-4c29-879b-3b10887a616d
begin
	struct BondDefault
		x
		default
	end

	Base.get(bd::BondDefault) = bd.default
	Base.show(io::IO, m::MIME"text/html", bd::BondDefault) = Base.show(io, m, bd.x)

	BondDefault
end

# ╔═╡ 29aa932b-9835-4d13-84e2-5ccf380a21ea
@bind which_graph_2 Select([
		"Emissions"
		"Concentrations"
		"Temperature"
		])

# ╔═╡ 9d603716-3069-4032-9416-cd8ab2e272c6
@bind which_graph_4 Select([
		"Emissions"
		"Concentrations"
		"Temperature"
		"Costs and benefits"
])

# ╔═╡ 70173466-c9b5-4227-8fba-6256fc1ecace
Tmax_9_slider = @bind Tmax_9 Slider(0:0.1:5; default=2);

# ╔═╡ 6bcb9b9e-e0ab-45d3-b9b9-3d7282f89df6
allow_overshoot_9_cb = @bind allow_overshoot_9 CheckBox();

# ╔═╡ b428e2d3-e1a9-4e4e-a64f-61048572102f
function multiplier(unit::Real, factor::Real=2, suffix::String="%")
	h = @htl("""

	""")

	BondDefault(h, unit)
end

# ╔═╡ 8cab3d28-a457-4ccc-b053-38cd003bf4d1
function Carousel(
		elementsList;
		wraparound::Bool=false,
		peek::Bool=true,
	)

	@assert peek

    carouselHTML = map(elementsList) do element
        @htl("""

            $(element)

""")
    end

    h = @htl("""

            $(carouselHTML)

        <button data-value="-1">Previous</button>
        <button data-value="1">Next</button>

    """)

	BondDefault(h,1)
end

# ╔═╡ c7cbc172-daed-406f-b24b-5da2cc234c29
preindustrial_concentrations = 280

# ╔═╡ b440cd13-36a9-4c54-9d80-ac3fa7c2900e
end_of_oil = 2150 # cannot mitigate when fossil fuels are depleted

# ╔═╡ ec760706-15ac-4a50-a67e-c338d70f3b0a
pp = (;
	((k, (:color => c, :label => n))
	for (k, c, n) in zip(keys(names), colors, names))...
);

# ╔═╡ bb4b25e4-0db5-414b-a384-0a27fe7efb66
gauss_stdev = 30

# ╔═╡ 013807a0-bddb-448b-9300-f7f559e48a45
begin
	default_usage_error = :(error("Example usage:\n\n@intially [1,2] @bind x f(x)\n"))

	macro initially(::Any)
		default_usage_error
	end

	macro initially(default, bind_expr::Expr)
		if bind_expr.head != :macrocall || bind_expr.args[1] != Symbol("@bind")
			return default_usage_error
		end

		# warn if the first argument is a @bind
		if default isa Expr && default.head == :macrocall && default.args[1] == Symbol("@bind")
			return default_usage_error
		end

		esc(intially_function(default, bind_expr))
	end

	function intially_function(default, bind_expr)
		sym = bind_expr.args[3]
		@gensym setval bond

		quote
			if !@isdefined($sym)
				$sym = $default
			end

			$setval = $sym

			$bond = @bind $sym $(bind_expr.args[4])
			PlutoRunner.Bond

			if $sym isa Missing
				$sym = $setval
			end

			$bond
		end
	end
end

# ╔═╡ 4e91fb48-fc5e-409e-9a7e-bf846f1d211d
html"""

"""

# ╔═╡ 3c7271ab-ece5-4ae2-a8dd-dc3670f300f7
# initial_mrga_1 = Dict(
# 	"M" => [2070, 0.7],
# 	"R" => [2100, 0.4],
# 	"G" => [2170, 0.3],
# 	"A" => [2110, 0.1],
# )

# ╔═╡ dcf265c1-f09b-483e-a361-d54c6c7500c1
# @initially initial_mrga_1 @bind input_8 begin

# 	controls_8 = MRGA(
# 		gaussish(input_8["M"]...),
# 		gaussish(input_8["R"]...),
# 		gaussish(input_8["G"]...),
# 		gaussish(input_8["A"]...),
# 	)

# 	plotclicktracker2(
# 		plot_controls(controls_8),
# 		initial_mrga_1,
# 	)
# end

# ╔═╡ 10c015ec-780c-4453-83cb-12dd0f09f358
function plotclicktracker(p::Plots.Plot; draggable::Bool=false)

	# we need to render the plot before its dimensions are available:
	# plot_render = repr(MIME"image/svg+xml"(),  p)
	plot_render = repr(MIME"image/svg+xml"(),  p)

	# these are the _bounding boxes_ of our plot
	big = bbox(p.layout)
	small = plotarea(p[1])

	# the axis limits
	xl = xlims(p)
	yl = ylims(p)

	# with this information, we can form the linear transformation from
	# screen coordinate -> plot coordinate

	# this is done on the JS side, to avoid one step in the Julia side
	# we send the linear coefficients:
	r = (
		x_offset = xl[1] - (xl[2] - xl[1]) * small.x0[1] / small.a[1],
		x_scale = (big.a[1] / small.a[1]) * (xl[2] - xl[1]),
		y_offset = (yl[2] - yl[1]) + (small.x0[2] / small.a[2]) * (yl[2] - yl[1]) + yl[1],
		y_scale = -(big.a[2]/ small.a[2]) * (yl[2] - yl[1]),
		x_min = xl[1], # TODO: add margin
		x_max = xl[2],
		y_min = yl[1],
		y_max = yl[2],
	)

	HTML("""""")
end

# ╔═╡ 7e540eaf-8700-4176-a96c-77ee2e4c384b
years = 2020:12.0:2200

# ╔═╡ 646591c4-cb60-41cd-beb9-506807ce17d2
function gaussish(mean, magnitude)
	my_stdev = gauss_stdev * (1 + magnitude);
	map(years) do t
	clamp(
		(1.5 *
		magnitude *
		(-0.4 +
		  exp(
			(-1 * ((t - mean) * (t - mean))) / (2 * my_stdev * my_stdev)
		  ))) /
		(1.0 - 0.4),
		0.0,
		1.0
	)
	end
end

# ╔═╡ 8fa94ec9-1fab-41b9-a7e6-1917e975e4ff
function default_parameters()::ClimateModelParameters
	result = deepcopy(ClimateMARGO.IO.included_configurations["default"])
	result.domain = years isa Domain ? years : Domain(step(years), first(years), last(years))
	result.economics.baseline_emissions = ramp_emissions(result.domain)
    result.economics.extra_CO₂ = zeros(size(result.economics.baseline_emissions))
	return result
end

# ╔═╡ 785c428d-d4f7-431e-94d7-039b0708a78a
function opt_controls_temp(model_parameters = default_parameters(); opt_parameters...)

    model = ClimateModel(model_parameters)

    model_optimizer = optimize_controls!(model; opt_parameters..., print_raw_status=false)

	(
		result=model,
		status=ClimateMARGO.Optimization.JuMP.termination_status(model_optimizer) |> string,
	)
    # return Dict(
    #     :model_parameters => model_parameters,
    #     model_results(model)...
    # )
end

# ╔═╡ 254ce01c-7976-4fe8-a980-fea1a61d7406

# ╔═╡ 2dcd5669-c725-40b9-84c4-f8399f6e924b
bigbreak = html"""

""";

# ╔═╡ eb7f34c3-1cd9-411d-8f34-d8547ba6ac29
bigbreak

# ╔═╡ cf5f7459-fbbe-4595-baa4-b6b85017005a
bigbreak

# ╔═╡ d796cb71-0e34-41ce-bc97-45c68812046d
bigbreak

# ╔═╡ 6ca2a32a-51de-4ff4-8023-843396f970ec
bigbreak

# ╔═╡ b8f9efec-63ac-4e58-93cf-9f7199b78451
function setfieldconvert!(value, name::Symbol, x)
    setfield!(value, name, convert(typeof(getfield(value, name)), x))
end

# ╔═╡ 371991c7-13dd-46f6-a730-ad89f43c6f0e
function enforce_maxslope!(controls;
		dt=step(years),
		max_slope=Dict("mitigate"=>1. /40., "remove"=>1. /40., "geoeng"=>1. /80., "adapt"=> 0.)
    )
    controls.mitigate[1] = 0.0
    controls.remove[1] = 0.0
    controls.geoeng[1] = 0.0
    # controls.adapt[1] = 0.0

    for i in 2:length(controls.mitigate)
        controls.mitigate[i] = clamp(
            controls.mitigate[i],
            controls.mitigate[i-1] - max_slope["mitigate"]*dt,
            controls.mitigate[i-1] + max_slope["mitigate"]*dt
        )
        controls.remove[i] = clamp(
            controls.remove[i],
            controls.remove[i-1] - max_slope["remove"]*dt,
            controls.remove[i-1] + max_slope["remove"]*dt
        )
        controls.geoeng[i] = clamp(
            controls.geoeng[i],
            controls.geoeng[i-1] - max_slope["geoeng"]*dt,
            controls.geoeng[i-1] + max_slope["geoeng"]*dt
        )
        controls.adapt[i] = clamp(
            controls.adapt[i],
            controls.adapt[i-1] - max_slope["adapt"]*dt,
            controls.adapt[i-1] + max_slope["adapt"]*dt
        )
    end
end

# ╔═╡ b7ca316b-6fa6-4c2e-b43b-cddb08aaabbb
function costs_dict(costs, model)
    Dict(
        :discounted => costs,
        :total_discounted => sum(costs .* model.domain.dt),
    )
end

# ╔═╡ 0a3be2ea-6af6-43c0-b8fb-e453bc2b703b

model_results(model::ClimateModel) = Dict(
    :controls => model.controls,
    :computed => Dict(
        :temperatures => Dict(
            :baseline => T_adapt(model),
            :M => T_adapt(model; M=true),
            :MR => T_adapt(model; M=true, R=true),
            :MRG => T_adapt(model; M=true, R=true, G=true),
            :MRGA => T_adapt(model; M=true, R=true, G=true, A=true),
        ),
        :emissions => Dict(
            :baseline => effective_emissions(model),
            :M => effective_emissions(model; M=true),
            :MRGA => effective_emissions(model; M=true, R=true),
        ),
        :concentrations => Dict(
            :baseline => c(model),
            :M => c(model; M=true),
            :MRGA => c(model; M=true, R=true),
        ),
        :damages => Dict(
            :baseline => costs_dict(damage(model; discounting=true), model),
            :MRGA => costs_dict(damage(model; M=true, R=true, G=true, A=true, discounting=true), model),
        ),
        :costs => Dict(
            :M => costs_dict(cost(model; M=true, discounting=true), model),
            :R => costs_dict(cost(model; R=true, discounting=true), model),
            :G => costs_dict(cost(model; G=true, discounting=true), model),
            :A => costs_dict(cost(model; A=true, discounting=true), model),
            :MRGA => costs_dict(cost(model; M=true, R=true, G=true, A=true, discounting=true), model),
        ),
    ),
)

# ╔═╡ eb0c961d-42cf-4219-a36e-cd492fa31f6b
const cost_bars_scale = 70

# ╔═╡ ec5d87a6-354b-4f1d-bb73-b3db08589d9b
total_discounted(costs, model) = sum(costs .* model.domain.dt)

# ╔═╡ 70f01a4d-0aa3-4cd5-ad71-452c490c61ac
colors_js = Dict((k,string("#", hex(v))) for (k,v) in pairs(colors));

# ╔═╡ ac779b93-e19e-41de-94cb-6a2a919bcd2e
names_js = Dict(pairs(names));

# ╔═╡ 5c484595-4646-484f-9e75-a4a3b4c2af9b
function plotclicktracker2(p::Plots.Plot, initial::Dict; draggable::Bool=true)

	# we need to render the plot before its dimensions are available:
	# plot_render = repr(MIME"image/svg+xml"(),  p)
	plot_render = repr(MIME"image/svg+xml"(),  p)

	# these are the _bounding boxes_ of our plot
	big = bbox(p.layout)
	small = plotarea(p[1])

	# the axis limits
	xl = xlims(p)
	yl = ylims(p)

	# with this information, we can form the linear transformation from
	# screen coordinate -> plot coordinate

	# this is done on the JS side, to avoid one step in the Julia side
	# we send the linear coefficients:
	r = (
		x_offset = xl[1] - (xl[2] - xl[1]) * small.x0[1] / small.a[1],
		x_scale = (big.a[1] / small.a[1]) * (xl[2] - xl[1]),
		y_offset = (yl[2] - yl[1]) + (small.x0[2] / small.a[2]) * (yl[2] - yl[1]) + yl[1],
		y_scale = -(big.a[2]/ small.a[2]) * (yl[2] - yl[1]),
		x_min = xl[1], # TODO: add margin
		x_max = xl[2],
		y_min = yl[1],
		y_max = yl[2],
		aspect_ratio = big.a[1] / big.a[2],
	)

	@htl("""""")
end

# ╔═╡ 060cbeab-5503-4eda-95d8-3f554765b2ee
begin
	mutable struct MRGA{T}
		M::T
		R::T
		G::T
		A::T
	end
	# function MRGA(M::TM,R::TR=nothing,G::TG=nothing) where {TM,TR,TG}
	# 	MRGA{Union{TM,TR,TG,Nothing}}(M,R,G,nothing)
	# end
	function MRGA(;M::TM=nothing,R::TR=nothing,G::TG=nothing,A::TA=nothing) where {TM,TR,TG,TA}
		MRGA{Union{TM,TR,TG,TA}}(M,R,G,A)
	end

	MRGA(x) = MRGA(x,x,x,x)

	splat(mrga::MRGA) = (:M => mrga.M, :R => mrga.R, :G => mrga.G, :A => mrga.A)

	Base.collect(mrga::MRGA) = [mrga.M, mrga.R, mrga.G, mrga.A]

	Base.getindex(x::MRGA, pos::MRGA{Bool}) = collect(x)[collect(pos)]

	Base.getindex(x::MRGA, tech::Symbol) = getfield(x, tech)

	Base.eachindex(m::MRGA) = (:M,:R,:G,:A)

	Base.enumerate(mrga::MRGA) = ((:M, mrga.M), (:R, mrga.R), (:G, mrga.G), (:A, mrga.A))

	Base.any(m::MRGA{Bool}) = m.M || m.R || m.G || m.A

	Base.all(m::MRGA{Bool}) = m.M && m.R && m.G && m.A

	MRGA
end

# ╔═╡ 331c45b7-b5f2-4a78-b180-5b918d1806ee
md"""
# Emissions mitigation and carbon dioxide removal to minimize climate suffering

This interactive article lets *you*– the reader– run [MARGO](https://github.com/ClimateMARGO/ClimateMARGO.jl), a simple climate model, to explore what it takes to avoid the catastrophic impacts of global warming. The code in this webpage is *reactive*, meaning that the graphs and numbers automatically update whenever you change the climate model's inputs.

*By [Henri F. Drake](https://hdrake.github.io/) and Fons van der Plas*

## _Can you limit human-caused global warming to "well-below 2ºC"?_

"""

# ╔═╡ 6533c123-34fe-4c0d-9ecc-7fef11379253
md"""
![image](https://user-images.githubusercontent.com/6933510/118835384-3ad36c80-b8c3-11eb-995d-70cba3b23846.png)

_From: [ClimateMARGO.jl](https://github.com/ClimateMARGO/ClimateMARGO.jl)_
"""

# ╔═╡ 50d24c91-61ae-4544-98fa-5749bafe3d41
md"""
## Overview of the climate problem: from greenhouse gas emissions to climate suffering

Human emissions of greenhouse gases, especially Carbon Dioxide (CO₂), increase the stock of greenhouse gases in the atmosphere. For every molecule of CO₂ emitted, about 50% are taken up by plants, soils, or the ocean within a few years, while the rest remains in the atmosphere. (The effects of other greenhouse gases, such as Methane and CFCs, and other forcing agents, can approximately be converted into the "CO₂-equivalent"– or CO₂ₑ– concentrations that would lead to the same climate forcing).

Greenhouse gases get their name because they trap invisible heat radiation emitted by Earth's surface and atmosphere from escaping to space, much like greenhouses trap hot air from rising when it is warmed by the Sun. This "greenhouse effect" causes the temperature to rise globally, although some places warm *more* and *faster* than others. Warmer temperatures exacerbate both the frequency and intensity of "natural" disasters, such as heat waves, coastal flooding from major hurricanes, and inland flooding from torrential rain. These climate impacts lead to enhanced climate suffering, which economics typically attempt to quantify suffering in terms of lost money or welfare.

In the interactive article below, we invite you to explore the benefits of emissions mitigation and carbon dioxide removal in reducing climate suffering, and the trade-offs with their costs.
"""

# ╔═╡ ec325089-8418-4fed-ac0e-e8ae21b433ab
md"""
## Mitigating emissions
Human greenhouse gas emissions are the result of fossil fuel burning (e.g. for transportation, electricity generation, heating, industry), unsustainable agriculture, and land use changes. We refer to any actions or policies that reduce these emissions as *mitigation*.

The MARGO model lumps all potential mitigation into a single number: the percentage of *baseline* emissions that are mitigated in a given year. Baseline emissions are the emissions that would arise in a hypothetical future world absent of climate policy. In our hypothetical no-policy world, we assume that emissions will go to zero by 2150 even without climate policy, perhaps because of public health concerns regarding other forms of air pollution, the development of new zero-carbon technologies, or running out of extractable fossil fuels resources.

*In the plot below, drag the blue dot around* to vary the amount and timing of mitigation, and observe how these changes affect key climate variables, using the drop-box menu: CO₂ₑ emissions, CO₂ₑ concentrations, and global temperature.
"""

# ╔═╡ 30218715-6469-4a0f-bf90-f3243219e7b5
md"""
## Cost & damages
"""

# ╔═╡ 8433cb38-915a-46c1-b3db-8e7905351c1b
@bind cost_benefits_narrative_slide Carousel([
		  md"""
		### 1. The costs of climate suffering

		In the absence of climate action, temperatures would rise over 4.5ºC above preindustrial levels (1800 to 1850 average), causing catastrophic climate impacts. MARGO attempts to quantify this suffering by translating the degree of warming into economic damages (in \$ / year). The curve below shows how climate damages rise over time, as a percentage of the World Gross Domestic Product (WGDP) in that year, due to uncontrolled temperature increases.

		""",

		md"""### 2. Avoiding climate damages
		Emissions mitigation limits future warming and climate suffering (_Damages_ curve). The economic benefits of mitigation are given by the difference in damages relative to the no-policy scenario (_Baseline_ curve minus _Damages_ curve).

		In the figure below, drag around the blue dot to change the future mitigation strategy, and observe how the _Avoided damages_ (the grey area) change!

		""",

		md"""### 3. Cost-benefit analysis

		Unfortunately, mitigating CO₂ₑ emissions also carries a cost. In MARGO, the *marginal* cost of mitigation is proportional to the fraction of CO₂ₑ emissions that have been mitigated in a given year, increasing up to a maximum of $70 per metric ton of CO₂ₑ at 100% mitigation.

		This naturally leads to a **cost-benefit analysis**. We search for the most beneficial, or *optimal*, scenario: the one with the *maximum net present benefits*. In the figure below, try finding a mitigation strategy that optimizes these _Net benefits_.
		"""
]; wraparound=false)

# ╔═╡ 11d62228-476c-4616-9e7d-de6c05a6a53d
if cost_benefits_narrative_slide == 1
	hidecloack("cost_benefits_narrative_input")
end

# ╔═╡ 4c7fccc5-450c-4903-96a6-ce36ff60d280
md"""
## Picking up the slack: carbon dioxide removal

While substantial emissions mitigations are necessary to reduce future climate suffering, they can not make up for the hundreds of billions of tons of CO₂ that humans have already emitted. However, both natural and technological methods for removing CO₂ from the atmosphere exist. Although they are presently miniscule compared to the tens-of-gigatons scale of global emissions, experts expect that they will play a key role in the future. In MARGO, we do not distinguish between different carbon dioxide removal methods, and further assume that the carbon is stored permanently.

*Drag the yellow dot in the figure below to modify the amount and timing of carbon dioxide removal*.
"""

# ╔═╡ b2d65726-df99-4710-9d03-9f6838036c87
md"""
## MARGO's automated optimization

In the above example, *you* manually adjusted the timing and amount of mitigation and carbon dioxide removal, but did not have much control on the shape of curves. Using a computer algorithm, we can do this optimization step *automatically* and *faster*, without having to assume anything about the shape of the mitigation and carbon dioxide removal curves.
"""

# ╔═╡ a0a1bb20-ec9b-446d-a36a-272840b8d35c
blob(
	md"""
	#### Maximum temperature

	`0.0 °C` $(Tmax_9_slider) `5.0 °C`

	_Allow **temperature overshoot**:_ $(allow_overshoot_9_cb)

	""",
	"#c5710014"
)

# ╔═╡ 4e1524e9-27d2-45a7-9ebd-21e6e369c4a7
md"""
# Appendix

This interactive page is a [Pluto notebook](https://github.com/fonsp/Pluto.jl), which means that it contains its own source code, and we encourage you to take a look behind the scenes. Use the button "Edit or run" on the top right to modify the notebook and reproduce the result yourself.
"""

# ╔═╡ d5ab2398-2904-4b44-a855-e4729996bc0d
md"""
Example of running the model:
"""

# ╔═╡ 14623e1f-7719-47b1-8854-8070d5ef8e17
md"""
## Plot functions
"""

# ╔═╡ a9b1e7fa-0318-41d8-b720-b8615c047bcd
plot_controls(c::ClimateMARGO.Models.Controls) = plot_controls(MRGA(
		c.mitigate,
		c.remove,
		c.geoeng,
		c.adapt
))

# ╔═╡ d9d20714-0689-449f-8e52-603dc804c93f
yearticks = collect(2020:20:2200)

# ╔═╡ cabc3214-1036-433b-aae1-6964bb780be8
function finish!(p)
	plot!(p;
		xlim=(2020,2201),
		xticks=yearticks,
		size=(680,200),
		grid=false,
	)
end

# ╔═╡ 2fec1e12-0218-4e93-a6b5-3711e6910d79
function plot_costs(result::ClimateModel;
		show_baseline::Bool=true,
		show_controls::Bool=true,
		show_damages::Bool=true,
		title="Control costs & climate damages"
	)

	p = plot(;
		ylim=(0,25.1), left_margin=2.25Plots.Measures.mm,
		ylabel="trillion USD / year",
	)
	title === nothing || plot!(p; title=title)

	# baseline
	show_baseline && plot!(p,
		years, damage(result; discounting=true);
		pp.baseline_damages...,
		fillrange=zero(years),
		fillopacity=.2,
		linestyle=:dash,
	)

	# control costs
	controlled_damages = damage(result; M=true, R=true, G=true, A=true, discounting=true)

	show_controls && plot!(p,
		years, controlled_damages .+ cost(result; M=true, R=true, G=true, A=true, discounting=true);
		fillrange=controlled_damages,
		fillopacity=.2,
		pp.controls...
	)

	# controlled damages
	show_damages && plot!(p,
		years, controlled_damages;
		fillrange=zero(years),
		fillopacity=.2,
		pp.damages...
	)

	finish!(p)

end

# ╔═╡ cff9f952-4850-4d55-bb8d-c0a759d1b7d8
function plot_concentrations(result::ClimateModel;
		relative_to_preindustrial::Bool=true)
	Tmax = 5
	p = relative_to_preindustrial ? plot(;
		ylim=(0,4.5),
		yticks=[1,2,3,4],
		yformatter=x -> string(Int(x), "×"),
		title="Atmospheric CO₂ₑ concentration, relative to 1800-1850",
	) : plot(;
		ylim=(0,1400),
		ylabel="ppm",
		title="Atmospheric CO₂ₑ concentration",
	)

	factor = relative_to_preindustrial ? preindustrial_concentrations : 1

	# baseline
	plot!(p,
		years, c(result) ./ factor;
		pp.baseline_concentrations...,
		linestyle=:dash,
	)
	# controlled temperature
	plot!(p,
		years, c(result; M=true, R=true) ./ factor;
		pp.concentrations...
	)

	finish!(p)
end

# ╔═╡ c73c89a7-f652-4554-95e9-20f47a818996
function plot_controls(controls::MRGA; title=nothing)

	p = plot(;
		ylim=(0,1),
	)
	title === nothing || plot!(p; title=title)

	for (tech, c) in enumerate(controls)
		if c !== nothing
			which = tech === :M ? (years .< end_of_oil) : eachindex(years)

			plot!(p,
				years[which], c[which];
				pp[tech]...
			)
		end
	end

	finish!(p)

end

# ╔═╡ 6634bcf1-8af6-4000-9b00-a5b4c02596c6
function plot_emissions(result::ClimateModel)

	p = plot(;
		ylim=(-3,11), left_margin=1.75Plots.Measures.mm,
		ylabel="ppm / year",
		title="Global CO₂ₑ emissions",
	)

	# baseline
	plot!(p,
		years, effective_emissions(result);
		pp.baseline_emissions...,
		linestyle=:dash,
	)
	# controlled
	plot!(p,
		years, effective_emissions(result; M=true, R=true);
		fillrange=zero(years),
		fillopacity=.2,
		pp.emissions...
	)

	finish!(p)

end

# ╔═╡ 424940e1-06ef-453a-8ffb-deb24dadb334
function plot_emissions_pretty(result::ClimateModel)
	# offset the x values so that framestyle=:origin will make the y-axis pass through 2020 instead of 0. yuck
	R = x -> x + 2020
	L = x -> x - 2020

	Tmax = 5
	p = plot(;
		ylim=(-3,11),
		ylabel="ppm / year",
		framestyle = :origin,
		xformatter=string ∘ Int ∘ R,
	)

	# baseline
	plot!(p,
		L.(years), effective_emissions(result);
		pp.baseline_emissions...,
		linestyle=:dash,
	)
	# controlled temperature
	plot!(p,
		L.(years), effective_emissions(result; M=true, R=true);
		fillrange=zero(L.(years)),
		fillopacity=.2,
		pp.emissions...
	)

	finish!(p)

	plot!(p;
		xlim=L.(extrema(years)),
		xticks = L.(yearticks),
		)
end

# ╔═╡ 700f982d-85da-4dc1-9319-f3b2527d0308
function plot_temp(result::ClimateModel)
	Tmax = 5

	# setup
	p = plot(;
		ylim=(0,Tmax),
		yticks=[0,1,2,3],
		yformatter=x -> string("+", Int(x), " °C"),
		title="Global warming relative to 1800-1850",
	)

	# shade dangerously high temperatures
	for a in [2,3]
		plot!(p,
			collect(extrema(years)),
			[a,a],
			linewidth=0,
			label=nothing,
			fillrange=[Tmax,Tmax],
			fillcolor=colors.above_paris
		)
	end

	# baseline
	plot!(p,
		years, T_adapt(result; splat(MRGA(false))...);
		pp.baseline_temperature...,
		linestyle=:dash,
	)
	# controlled temperature
	plot!(p,
		years, T_adapt(result; splat(MRGA(true))...);
		pp.temperature...
	)

	finish!(p)
end

# ╔═╡ ab557633-e0b5-4439-bc81-d274770f2e65
md"""
## Plot dots input magic
"""

# ╔═╡ 2758b185-cd54-484e-bb7d-d4cfcd2d39f4
md"""
## Running the model
"""

# ╔═╡ 611c25ab-a454-4d52-b8fb-a58b0d1f5ca6
function forward_controls_temp(controls::MRGA=MRGA(), model_parameters=default_parameters())

    model = ClimateModel(model_parameters)

    translations = Dict(
        :M => :mitigate,
        :R => :remove,
        :G => :geoeng,
        :A => :adapt,
    )
    for (k, v) in enumerate(controls)
		if v !== nothing
	        setfieldconvert!(model.controls, translations[Symbol(k)], copy(v))
		end
    end

    enforce_maxslope!(model.controls)

	model
    # return Dict(
    #     :model_parameters => model_parameters,
    #     model_results(model)...
    # )
end

# ╔═╡ 9aa73ce0-cec6-4d53-bbbc-f5c85de7b521
@initially initial_1 @bind input_1 begin

	local t = input_1["M"][1]
	local y = input_1["M"][2]

	controls_1 = MRGA(
		M=gaussish(t, clamp(.07 * (10 - y), 0, 1)),
		R=gaussish(t, clamp(.07 * (10 - y) * 0.25, 0, 1)),
	)

	result_1 = forward_controls_temp(controls_1)

	plotclicktracker2(
		plot_emissions(result_1),
		initial_1
	)
end

# ╔═╡ 65d31fbf-322d-459a-a2dd-2894edbecc4d
plot_temp(result_1)

# ╔═╡ ff2b1c0a-e419-4f41-aa3b-d017642ffc13
@initially initial_2 @bind input_2 begin

	controls_2 = MRGA(
		M=gaussish(input_2["M"]...),
	)

	result_2 = forward_controls_temp(controls_2)

	plotclicktracker2(
		plot_controls(controls_2; title="Deployment of mitigation"),
		initial_2
	)
end

# ╔═╡ 02851ee9-8050-4821-b3c9-1f65c9b8135b
if which_graph_2 == "Emissions"
	plot_emissions(result_2)
elseif which_graph_2 == "Concentrations"
	plot_concentrations(result_2; relative_to_preindustrial=true)
else
	plot_temp(result_2)
end

# ╔═╡ f4203dcf-b251-4e2b-be07-922bc7c4496d
(@initially initial_3 @bind input_3 begin

	controls_3 = MRGA(
		M=gaussish(input_3["M"]...),
	)

	result_3 = forward_controls_temp(controls_3)

	plotclicktracker2(
		plot_controls(controls_3; title="Deployment of mitigation"),
		initial_3
	)
end) |> cloak("cost_benefits_narrative_input")

# ╔═╡ 3e26d311-6abc-4b2c-ada4-f8a3171d9f75
if cost_benefits_narrative_slide == 1
	local uncontrolled = ClimateModel(default_parameters())
	plot_costs(uncontrolled; show_controls=false, show_baseline=false)
elseif cost_benefits_narrative_slide == 2
	plot_costs(result_3; show_controls=false)
else
	plot_costs(result_3)
end

# ╔═╡ aac86adf-465f-464f-b258-406c2e55b82f
@initially initial_4 @bind input_4 begin

	controls_4 = MRGA(
		M=gaussish(input_4["M"]...),
		R=gaussish(input_4["R"]...),
	)

	result_4 = forward_controls_temp(controls_4)

	plotclicktracker2(
		plot_controls(controls_4; title="Deployment of mitigation"),
		initial_4
	)
end

# ╔═╡ a751fb75-952e-41d4-a8b5-aba512c10e55
if which_graph_4 == "Emissions"
	plot_emissions(result_4)
elseif which_graph_4 == "Concentrations"
	plot_concentrations(result_4; relative_to_preindustrial=true)
elseif which_graph_4 == "Temperature"
	plot_temp(result_4)
else
	plot_costs(result_4)
end

# ╔═╡ 89752d91-9c8e-4203-b6f1-bdad41386b31
shortname = MRGA("M","R","G","A")

# ╔═╡ ff2709a4-516f-4066-b5b2-617ac0e5f20c
mediumname = MRGA("mitigate", "remove", "geoeng", "adapt")

# ╔═╡ 2821b722-75c2-4072-b142-d13553a84b7b
longname = MRGA("Mitigation", "Removal", "Geo-engineering", "Adaptation")

# ╔═╡ 8e89f521-c19d-4f87-9497-f9b61c19c176
let
	default_cost = let
		e = default_parameters().economics
		MRGA(e.mitigate_cost, e.remove_cost, e.geoeng_cost, e.adapt_cost)
	end
	blob(
		@htl("""

Which controls?

		Enabled?	Cost multiplier

		$(longname.M)
		$(@bind enable_M_9 CheckBox(;default=true))
		$(@bind cost_M_9 multiplier(default_cost.M, 4))

		$(longname.R)
		$(@bind enable_R_9 CheckBox(;default=true))
		$(@bind cost_R_9 multiplier(default_cost.R, 4))

		$(longname.G)
		$(@bind enable_G_9 CheckBox(;default=false))
		$(@bind cost_G_9 multiplier(default_cost.G, 4))

		$(longname.A)
		$(@bind enable_A_9 CheckBox(;default=false))
		$(@bind cost_A_9 multiplier(default_cost.A, 4))

		"""),
		"#c500b40a"
	)
end

# ╔═╡ a83e47fa-4b48-4bbc-b210-382d1cf19f55
control_enabled_9 = MRGA(
	enable_M_9,
	enable_R_9,
	enable_G_9,
	enable_A_9,
);

# ╔═╡ 242f3109-244b-4884-a0e9-6ea8950ca47e
control_cost_9 = MRGA(
	Float64(cost_M_9),
	Float64(cost_R_9),
	Float64(cost_G_9),
	Float64(cost_A_9),
);

# ╔═╡ f861935a-8b03-426e-aebe-6963e034ad49
output_9 = let
	parameters = default_parameters()

	parameters.economics.mitigate_cost = control_cost_9.M
	parameters.economics.remove_cost = control_cost_9.R
	parameters.economics.geoeng_cost = control_cost_9.G
	parameters.economics.adapt_cost = control_cost_9.A

	# modify the parameters here!

	opt_controls_temp(parameters;
		temp_overshoot=allow_overshoot_9 ? 999.0 : Tmax_9,
		temp_goal=Tmax_9,
		max_deployment=let
			e = control_enabled_9
			Dict(
				"mitigate" => e.M ? 1.0 : 0.0,
				"remove" => e.R ? 1.0 : 0.0,
				"geoeng" => e.G ? 1.0 : 0.0,
				"adapt" => e.A ? 0.4 : 0.0,
			)
		end,
	)
end

# ╔═╡ 6978acad-9cac-4490-85fb-7e43d9558aca
plot_controls(output_9.result.controls) |> feasibility_overlay(output_9)

# ╔═╡ 7a435e46-4f36-4037-a9a6-d296b20bf6ac
plot!(plot_temp(output_9.result),
	years, zero(years) .+ Tmax_9;
	lw=2,
	pp.T_max...
	)

# ╔═╡ 7ffad0f8-082b-4ca1-84f7-37c08d5f7266
md"""
## Cost bars
"""

# ╔═╡ 608b50e7-4419-4dfb-8d9e-5144d4034c05
function avoided_damages_bars(result)
	td(x) = total_discounted(x, result)

	baseline_damages = td(damage(result; discounting=true))
	controlled_damages = td(damage(result; splat(MRGA(true))..., discounting=true))

	avoided_damages = baseline_damages - controlled_damages

	costs = td(cost(result; splat(MRGA(true))..., discounting=true))

	@htl("""

		""")
end

# ╔═╡ 2c1416cf-9b6b-40a0-b714-16853c7e1f1d
if cost_benefits_narrative_slide >= 2
	avoided_damages_bars(result_3)
end

# ╔═╡ 31a30755-1d8b-451b-8c9a-2c32a3a1d0b4
function cost_bars(result; offset_damages=false)
	td(x) = total_discounted(x, result)

	baseline_damages = td(damage(result; discounting=true))
	controlled_damages = td(damage(result; splat(MRGA(true))..., discounting=true))

	avoided_damages = baseline_damages - controlled_damages

	costs = td(cost(result; splat(MRGA(true))..., discounting=true))

	@htl("""

		""")
end

# ╔═╡ 470d2f6f-fe97-4edd-8aaa-142bc8046fe8
cost_bars(result_1)

# ╔═╡ 5154aac7-812d-447f-8435-b8209d45fe04
if cost_benefits_narrative_slide >= 3
	cost_bars(result_3; offset_damages=true)
else
	bigbreak
end

# ╔═╡ 7f9df132-61de-4fec-a674-176c4a43335c
md"""
## MRGA struct
"""

# ╔═╡ 354b9d8a-7c3f-456b-9da9-4396ac975743
function MR(x::T,y::T) where T
	MRGA{T}(x, y, zero(x), zero(x))
end

# ╔═╡ Cell order:
# ╟─9a48a08e-7281-473c-8afc-7ad3e0771269
# ╟─331c45b7-b5f2-4a78-b180-5b918d1806ee
# ╟─9aa73ce0-cec6-4d53-bbbc-f5c85de7b521
# ╟─65d31fbf-322d-459a-a2dd-2894edbecc4d
# ╟─470d2f6f-fe97-4edd-8aaa-142bc8046fe8
# ╟─94415ff2-32a2-4b0f-9911-3b93e202f548
# ╟─eb7f34c3-1cd9-411d-8f34-d8547ba6ac29
# ╟─6533c123-34fe-4c0d-9ecc-7fef11379253
# ╟─50d24c91-61ae-4544-98fa-5749bafe3d41
# ╟─cf5f7459-fbbe-4595-baa4-b6b85017005a
# ╟─ec325089-8418-4fed-ac0e-e8ae21b433ab
# ╟─ff2b1c0a-e419-4f41-aa3b-d017642ffc13
# ╟─29aa932b-9835-4d13-84e2-5ccf380a21ea
# ╟─02851ee9-8050-4821-b3c9-1f65c9b8135b
# ╟─d796cb71-0e34-41ce-bc97-45c68812046d
# ╟─e810a90f-f964-4d7d-acdb-fc3a159dc12e
# ╟─30218715-6469-4a0f-bf90-f3243219e7b5
# ╟─8433cb38-915a-46c1-b3db-8e7905351c1b
# ╟─3e26d311-6abc-4b2c-ada4-f8a3171d9f75
# ╟─f4203dcf-b251-4e2b-be07-922bc7c4496d
# ╟─2c1416cf-9b6b-40a0-b714-16853c7e1f1d
# ╟─5154aac7-812d-447f-8435-b8209d45fe04
# ╟─11d62228-476c-4616-9e7d-de6c05a6a53d
# ╟─a3422533-2b78-4bc2-92bd-737da3c8982d
# ╟─4c7fccc5-450c-4903-96a6-ce36ff60d280
# ╟─aac86adf-465f-464f-b258-406c2e55b82f
# ╟─9d603716-3069-4032-9416-cd8ab2e272c6
# ╟─a751fb75-952e-41d4-a8b5-aba512c10e55
# ╟─6ca2a32a-51de-4ff4-8023-843396f970ec
# ╟─bb66d347-99be-4a95-8ba8-57dc9d33384b
# ╟─b2d65726-df99-4710-9d03-9f6838036c87
# ╟─70173466-c9b5-4227-8fba-6256fc1ecace
# ╟─6bcb9b9e-e0ab-45d3-b9b9-3d7282f89df6
# ╟─a0a1bb20-ec9b-446d-a36a-272840b8d35c
# ╟─8e89f521-c19d-4f87-9497-f9b61c19c176
# ╟─6978acad-9cac-4490-85fb-7e43d9558aca
# ╟─7a435e46-4f36-4037-a9a6-d296b20bf6ac
# ╟─a83e47fa-4b48-4bbc-b210-382d1cf19f55
# ╟─242f3109-244b-4884-a0e9-6ea8950ca47e
# ╟─4e1524e9-27d2-45a7-9ebd-21e6e369c4a7
# ╟─0a126a12-7b89-4f22-ad66-b64965cd1d49
# ╠═1c8d2d00-b7d9-11eb-35c4-47f2a2aa1593
# ╟─d5ab2398-2904-4b44-a855-e4729996bc0d
# ╠═f861935a-8b03-426e-aebe-6963e034ad49
# ╟─3094a9eb-074d-46c3-9c1e-0a9c94c6ad43
# ╟─b428e2d3-e1a9-4e4e-a64f-61048572102f
# ╟─0b31eac2-8efd-47cd-9571-a2053846343b
# ╟─ca104939-a6ca-4e70-a47a-1eb3c32db18f
# ╟─cf90139c-13d8-42a7-aba3-8c431e7854b8
# ╟─bd2bfa3c-a42e-4975-a543-84541f66b1c1
# ╟─8cab3d28-a457-4ccc-b053-38cd003bf4d1
# ╟─b81de514-2506-4243-8235-0b54dd4a7ec9
# ╟─73e01bd8-f56b-4bb5-a9a2-85ad223c9e9b
# ╟─ae92ba1f-5175-4704-8240-2de8432df752
# ╟─8ac04d55-9034-4c29-879b-3b10887a616d
# ╟─14623e1f-7719-47b1-8854-8070d5ef8e17
# ╟─2fec1e12-0218-4e93-a6b5-3711e6910d79
# ╟─cff9f952-4850-4d55-bb8d-c0a759d1b7d8
# ╟─c73c89a7-f652-4554-95e9-20f47a818996
# ╟─a9b1e7fa-0318-41d8-b720-b8615c047bcd
# ╟─6634bcf1-8af6-4000-9b00-a5b4c02596c6
# ╟─424940e1-06ef-453a-8ffb-deb24dadb334
# ╟─700f982d-85da-4dc1-9319-f3b2527d0308
# ╟─cabc3214-1036-433b-aae1-6964bb780be8
# ╟─d9d20714-0689-449f-8e52-603dc804c93f
# ╟─c7cbc172-daed-406f-b24b-5da2cc234c29
# ╟─b440cd13-36a9-4c54-9d80-ac3fa7c2900e
# ╟─ec760706-15ac-4a50-a67e-c338d70f3b0a
# ╟─ab557633-e0b5-4439-bc81-d274770f2e65
# ╟─bb4b25e4-0db5-414b-a384-0a27fe7efb66
# ╟─646591c4-cb60-41cd-beb9-506807ce17d2
# ╟─5c484595-4646-484f-9e75-a4a3b4c2af9b
# ╟─013807a0-bddb-448b-9300-f7f559e48a45
# ╟─4e91fb48-fc5e-409e-9a7e-bf846f1d211d
# ╟─3c7271ab-ece5-4ae2-a8dd-dc3670f300f7
# ╟─dcf265c1-f09b-483e-a361-d54c6c7500c1
# ╟─10c015ec-780c-4453-83cb-12dd0f09f358
# ╟─2758b185-cd54-484e-bb7d-d4cfcd2d39f4
# ╟─8fa94ec9-1fab-41b9-a7e6-1917e975e4ff
# ╟─611c25ab-a454-4d52-b8fb-a58b0d1f5ca6
# ╟─785c428d-d4f7-431e-94d7-039b0708a78a
# ╟─7e540eaf-8700-4176-a96c-77ee2e4c384b
# ╟─254ce01c-7976-4fe8-a980-fea1a61d7406
# ╟─89752d91-9c8e-4203-b6f1-bdad41386b31
# ╟─ff2709a4-516f-4066-b5b2-617ac0e5f20c
# ╟─2821b722-75c2-4072-b142-d13553a84b7b
# ╟─2dcd5669-c725-40b9-84c4-f8399f6e924b
# ╟─b8f9efec-63ac-4e58-93cf-9f7199b78451
# ╟─371991c7-13dd-46f6-a730-ad89f43c6f0e
# ╟─0a3be2ea-6af6-43c0-b8fb-e453bc2b703b
# ╟─b7ca316b-6fa6-4c2e-b43b-cddb08aaabbb
# ╟─7ffad0f8-082b-4ca1-84f7-37c08d5f7266
# ╟─eb0c961d-42cf-4219-a36e-cd492fa31f6b
# ╟─608b50e7-4419-4dfb-8d9e-5144d4034c05
# ╟─31a30755-1d8b-451b-8c9a-2c32a3a1d0b4
# ╟─ec5d87a6-354b-4f1d-bb73-b3db08589d9b
# ╠═70f01a4d-0aa3-4cd5-ad71-452c490c61ac
# ╠═ac779b93-e19e-41de-94cb-6a2a919bcd2e
# ╟─7f9df132-61de-4fec-a674-176c4a43335c
# ╟─060cbeab-5503-4eda-95d8-3f554765b2ee
# ╟─354b9d8a-7c3f-456b-9da9-4396ac975743

max_stars_repo_name: ClimateMARGO/interactive

max_issues_repo_name: ClimateMARGO/interactive

max_forks_repo_name: ClimateMARGO/interactive