text:
function default(main::ANY, s::ANY, data::ANY)
    _default_light = Vec3f0[Vec3f0(1.0,1.0,1.0), Vec3f0(0.1,0.1,0.1), Vec3f0(0.9,0.9,0.9), Vec3f0(20,20,20)]
    data = _default(main, s, copy(data))
    @gen_defaults! data begin # make sure every object has these!
        model = eye(Mat4f0)
        light = _default_light
        preferred_camera = :perspective
        is_transparent_pass = Cint(false)
    end
end

"""
Creates a default visualization for any value.
The defaults can be customized via the key word arguments and the style parameter.
The style can change the the look completely (e.g points displayed as lines, or particles),
while the key word arguments just alter the parameters of one visualization.
Always returns a context, which can be displayed on a window via view(::Context, [display]).
"""
visualize(main::ANY, s::Symbol=:default; kw_args...) = visualize(main, Style{s}(), Dict{Symbol, Any}(kw_args))::Context
visualize(main::ANY, s::Style, data::Dict) = assemble_shader(default(main, s, data))::Context
visualize(c::Composable, s::Symbol=:default; kw_args...) = Context(c)
visualize(c::Composable, s::Style, data::Dict) = Context(c)

visualize(c::Context, s::Symbol=:default; kw_args...) = c
visualize(c::Context, s::Style, data::Dict) = c
#

function _view(
        robj::RenderObject, screen = current_screen();
        camera = robj.uniforms[:preferred_camera],
        position = Vec3f0(2), lookat = Vec3f0(0)
    )
    global _camera_counter
    local camsym::Symbol # make things type stable
    mouseinside = screen.inputs[:mouseinside]
    ishidden = screen.hidden
    if isa(camera, Symbol)
        camsym = camera::Symbol
        if haskey(screen.cameras, camsym)
            real_camera = screen.cameras[camsym]
        elseif camsym == :perspective
            keep = map((a, b) -> !a && b, ishidden, mouseinside)
            real_camera = PerspectiveCamera(screen.inputs, position, lookat, keep = keep)
        elseif camsym == :fixed_pixel
            real_camera = DummyCamera(window_size = screen.area)
        elseif camsym == :orthographic_pixel
            keep = map((a, b) -> !a && b, ishidden, mouseinside)
            real_camera = OrthographicPixelCamera(screen.inputs, keep = keep)
        elseif camsym == :nothing
            push!(screen, robj, :nothing)
            return nothing
        else
            error("Camera symbol $camera not known")
        end
    elseif isa(camera, Camera)
        camsym = :custom
        real_camera = camera
    else
         error("$camera not a known camera type")
    end
    screen.cameras[camsym] = real_camera
    robj.uniforms[:resolution] = map(screen.area) do area
        Vec2f0(widths(area))
    end
    collect(real_camera, robj.uniforms)
    push!(screen, robj)
    nothing
end

_view(robjs::Vector, screen = current_screen(); kw_args...) = for robj in robjs
    _view(robj, screen; kw_args...)
end
_view(c::Composable, screen = current_screen(); kw_args...) = _view(extract_renderable(c), screen; kw_args...)


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using DataFrames, CSV
	using Random, Distributions, Plots
	Random.seed!(0)
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.13"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	sampleData = CSV.read(datadir("machine1.csv"), DataFrame; delim=",")[:, 1]
end;

# ╔═╡ 1aa41898-3afa-11eb-2f49-570853ce6435
begin
	n, N = length(sampleData), 10^6
	alpha = 0.1

	bootstrapSampleMeans = [mean(rand(sampleData, n)) for i in 1:N]
	Lmean = quantile(bootstrapSampleMeans, alpha/2)
	Umean = quantile(bootstrapSampleMeans, 1-alpha/2)

	bootstrapSampleMedians = [median(rand(sampleData, n)) for i in 1:N]
	Lmed = quantile(bootstrapSampleMedians, alpha/2)
	Umed = quantile(bootstrapSampleMedians, 1-alpha/2)

	Text("Bootstrap confidence interval for the mean: $((Lmean, Umean))")
end

# ╔═╡ 503bce92-3b2c-11eb-02b5-17d205eea2d3
Text("Bootstrap confidence interval for the median: $((Lmed, Umed))")

# ╔═╡ 503c061e-3b2c-11eb-1658-494999a7c6b2
begin
	stephist(bootstrapSampleMeans, bins=1000, c=:blue,
		normed=true, label="Sample \nmeans")
	plot!([Lmean, Lmean],[0,2], c=:black, ls=:dash, label="90% CI")
	plot!([Umean, Umean],[0,2],c=:black, ls=:dash, label="",
		xlims=(52,54), ylims=(0,2), xlabel="Sample Means", ylabel="Density")
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.13"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═1aa41898-3afa-11eb-2f49-570853ce6435
# ╠═503bce92-3b2c-11eb-02b5-17d205eea2d3
# ╠═503c061e-3b2c-11eb-1658-494999a7c6b2
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12

ext: jl


text:
using Seleroute

oneGbps = 1024
halfGbps = 512
traffic_scaling = 0.95 # Many topologies will reach a 100% load otherwise.

INCLUDE_OVERLOADED_TOPOLOGIES = false # Some topologies have a very low capacity for high demands; enable them by setting this to true. These topologies would need a very low traffic_scaling factor to have a load under 100%.

topologies = Vector{Topology}()

# First tested topology.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [2, 3, .9 * oneGbps] [2, 6, oneGbps] [3, 4, oneGbps] [3, 6, oneGbps] [5, 6, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [1, 6, halfGbps] [4, 1, halfGbps] [4, 5, halfGbps] [5, 2, halfGbps] [6, 4, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("basic", nodes_matrix, edges_matrix, traffic_matrix))

### "Olver", inspired by http://math.mit.edu/~olver/dynamic.pdf, Figure 1.
## Basic topology, with uniform capacities.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [1, 7, oneGbps] [2, 5, oneGbps] [2, 7, oneGbps] [3, 4, oneGbps] [3, 5, oneGbps] [3, 7, oneGbps] [4, 6, oneGbps] [4, 7, oneGbps] [5, 6, oneGbps] [5, 7, oneGbps] [6, 7, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [1, 6, halfGbps] [4, 1, halfGbps] [4, 5, halfGbps] [5, 2, halfGbps] [6, 4, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("olver", nodes_matrix, edges_matrix, traffic_matrix))

# Variant, with nonuniform capacities.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"]], [2, 1])
edges_matrix = [[1, 2, 5 * oneGbps] [1, 3, halfGbps] [1, 7, 2 * oneGbps] [2, 5, 0.21 * oneGbps] [3, 4, 4.96 * oneGbps] [3, 5, 2.47 * oneGbps] [3, 7, 5 * oneGbps] [4, 6, oneGbps] [4, 7, oneGbps] [5, 6, oneGbps] [5, 7, oneGbps] [6, 7, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [1, 6, halfGbps] [4, 1, halfGbps] [4, 5, halfGbps] [5, 2, halfGbps] [6, 4, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("olver-various", nodes_matrix, edges_matrix, traffic_matrix))

# Variant, with other capacities (lower).
if INCLUDE_OVERLOADED_TOPOLOGIES
    nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"]], [2, 1])
    edges_matrix = [[1, 2, oneGbps / 5] [1, 3, halfGbps] [2, 5, oneGbps / 0.21] [3, 4, oneGbps / 4.96] [3, 5, oneGbps / 2.47] [4, 6, oneGbps / 0.147] [5, 6, oneGbps / 5] [1, 7, 2 * oneGbps] [2, 7, 0.147 * oneGbps] [3, 7, 5 * oneGbps] [5, 7, oneGbps] [6, 7, oneGbps]]'
    traffic_matrix = Float64[[1, 3, halfGbps] [1, 6, halfGbps] [4, 1, halfGbps] [4, 5, halfGbps] [5, 2, halfGbps] [6, 4, halfGbps]]'
    traffic_matrix[:, 3] .*= traffic_scaling

    push!(topologies, Topology("olver-various-inv-partial", nodes_matrix, edges_matrix, traffic_matrix))
end

# Variant, with other capacities (higher).
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"]], [2, 1])
edges_matrix = [[1, 2, 5 * oneGbps] [1, 3, 2 * oneGbps] [2, 5, 0.21 * oneGbps] [3, 4, 4.96 * oneGbps] [3, 5, 2.47 * oneGbps] [4, 6, 0.147 * oneGbps] [5, 6, 5 * oneGbps] [1, 7, oneGbps / 2] [2, 7, oneGbps / 0.147] [3, 7, oneGbps / 5] [5, 7, oneGbps / 2] [6, 7, oneGbps / 2]]'
traffic_matrix = Float64[[1, 3, halfGbps] [1, 6, halfGbps] [4, 1, halfGbps] [4, 5, halfGbps] [5, 2, halfGbps] [6, 4, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("olver-various-inv", nodes_matrix, edges_matrix, traffic_matrix))

## Olver, but without the node G (it has a very high degree).
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [2, 5, oneGbps] [3, 4, oneGbps] [3, 5, oneGbps] [4, 6, oneGbps] [5, 6, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [1, 6, halfGbps] [4, 1, halfGbps] [4, 5, halfGbps] [5, 2, halfGbps] [6, 4, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("olver-without_G", nodes_matrix, edges_matrix, traffic_matrix))

# Variant, with nonuniform capacities.
if INCLUDE_OVERLOADED_TOPOLOGIES
    nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"]], [2, 1])
    edges_matrix = [[1, 2, 5 * oneGbps] [1, 3, halfGbps] [2, 5, 0.21 * oneGbps] [3, 4, 4.96 * oneGbps] [3, 5, 2.47 * oneGbps] [4, 6, 0.147 * oneGbps] [5, 6, 5 * oneGbps]]'
    traffic_matrix = Float64[[1, 3, halfGbps] [1, 6, halfGbps] [4, 1, halfGbps] [4, 5, halfGbps] [5, 2, halfGbps] [6, 4, halfGbps]]'
    traffic_matrix[:, 3] .*= traffic_scaling

    push!(topologies, Topology("olver-without_G-various", nodes_matrix, edges_matrix, traffic_matrix))
end

# Variant, with other capacities (lower).
if INCLUDE_OVERLOADED_TOPOLOGIES
    nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"]], [2, 1])
    edges_matrix = [[1, 2, oneGbps / 5] [1, 3, halfGbps] [2, 5, oneGbps / 0.21] [3, 4, oneGbps / 4.96] [3, 5, oneGbps / 2.47] [4, 6, oneGbps / 0.147] [5, 6, oneGbps / 5]]'
    traffic_matrix = Float64[[1, 3, halfGbps] [1, 6, halfGbps] [4, 1, halfGbps] [4, 5, halfGbps] [5, 2, halfGbps] [6, 4, halfGbps]]'
    traffic_matrix[:, 3] .*= traffic_scaling

    push!(topologies, Topology("olver-without_G-various-inv", nodes_matrix, edges_matrix, traffic_matrix))
end

### Topology Zoo
# Abilene (http://www.topology-zoo.org/files/Abilene.gml).
# Traffic: https://totem.info.ucl.ac.be/dataset.html
# (IntraTM-2005-02-23-21-00.xml); removed self-traffic.
nodes_matrix = permutedims([[1, "NY"] [2, "CH"] [3, "DC"] [4, "SE"] [5, "SU"] [6, "LA"] [7, "DE"] [8, "KA"] [9, "HO"] [10, "AT"] [11, "IN"]], [2, 1])
edges_matrix = [[1, 2, 10 * oneGbps] [1, 3, 10 * oneGbps] [2, 11, 10 * oneGbps] [3, 10, 10 * oneGbps] [4, 5, 10 * oneGbps] [4, 7, 10 * oneGbps] [5, 6, 10 * oneGbps] [5, 7, 10 * oneGbps] [6, 9, 10 * oneGbps] [7, 8, 10 * oneGbps] [8, 9, 10 * oneGbps] [8, 11, 10 * oneGbps] [9, 10, 10 * oneGbps] [10, 11, 10 * oneGbps]]'
traffic_matrix = Float64[[8, 4, 6240.1600] [8, 1, 26863.6533] [8, 5, 28226.0267] [8, 3, 47.3956] [8, 1, 2817.4489] [8, 7, 6438.9689] [8, 2, 7656.7911] [8, 6, 4.7111] [8, 9, 1867315.3956] [4, 8, 21866.3111] [4, 1, 4646580.7022] [4, 5, 41628.4000] [4, 3, 59591.7956] [4, 7, 1285829.3600] [4, 2, 1196886.1244] [4, 6, 23.1911] [4, 11, 1926.3911] [4, 9, 63727.9378] [1, 8, 598001.9289] [1, 4, 342504.4533] [1, 5, 14935.0756] [1, 3, 355844.4622] [1, 7, 16094485.4311] [1, 2, 2543737.4400] [1, 6, 3683.2800] [1, 11, 7058.1422] [1, 9, 171177.7600] [5, 8, 35060.1244] [5, 4, 3344.5778] [5, 1, 29604.1867] [5, 3, 86.2400] [5, 10, 5342.3289] [5, 7, 123.6089] [5, 2, 27998.6044] [5, 6, 606.5778] [5, 11, 26.6667] [5, 9, 778276.2578] [3, 8, 27451.4844] [3, 4, 1958030.1956] [3, 1, 843953.6711] [3, 5, 106844.8000] [3, 10, 65148.4178] [3, 7, 369093.9733] [3, 2, 207656.8711] [3, 6, 1749.4133] [3, 11, 24723.9644] [3, 9, 69562.1511] [10, 8, 549654.6933] [10, 1, 63275.7778] [10, 5, 767704.7467] [10, 3, 82480.2756] [10, 7, 95816.9867] [10, 2, 115608.1778] [10, 6, 137330.3022] [10, 11, 81214.3556] [10, 9, 1256926.2222] [7, 8, 41063.8311] [7, 4, 1366652.5244] [7, 1, 2448645.0933] [7, 5, 83908.4711] [7, 3, 310578.6222] [7, 10, 29714.8444] [7, 2, 941842.7822] [7, 6, 8392.0000] [7, 11, 46217.5644] [7, 9, 57574.9067] [2, 8, 43894.6400] [2, 4, 335900.4000] [2, 1, 158223.1644] [2, 5, 16902.1867] [2, 3, 78114.3022] [2, 10, 177692.9600] [2, 7, 74661.2089] [2, 6, 421889.7778] [2, 11, 13222.6756] [2, 9, 18675.6622] [6, 4, 1389.6178] [6, 1, 62036.3289] [6, 5, 193.7067] [6, 3, 2093.6889] [6, 10, 1693.0844] [6, 7, 8646.7733] [6, 2, 364485.7333] [6, 9, 80837.0489] [11, 8, 6725.6444] [11, 4, 4773.3067] [11, 1, 32356.2933] [11, 5, 24303.0578] [11, 3, 53.0489] [11, 10, 9530.7822] [11, 7, 4520.6133] [11, 2, 400782.9511] [11, 6, 0.3556] [11, 9, 248961.2711] [9, 8, 652315.8311] [9, 4, 34186.6578] [9, 1, 21203.9911] [9, 5, 13540.1333] [9, 3, 19.8044] [9, 10, 720.1956] [9, 7, 837.6711] [9, 2, 18402.5333] [9, 11, 9.4578]]'
traffic_matrix[:, 3] /= 1024 * 3 # Mbps
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("abilene", nodes_matrix, edges_matrix, traffic_matrix))

# ATT (http://www.topology-zoo.org/files/AttMpls.gml).
# No traffic data available (take the same as Abilene).
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"] [13, "M"] [14, "N"] [15, "O"] [16, "P"] [17, "Q"] [18, "R"] [19, "S"] [20, "T"] [21, "U"] [22, "V"] [23, "W"] [24, "X"] [25, "Y"] [26, "Z"]], [2, 1])
edges_matrix = [[1, 2, 10 * oneGbps] [1, 3, 10 * oneGbps] [1, 17, 10 * oneGbps] [2, 3, 10 * oneGbps] [3, 4, 10 * oneGbps] [3, 5, 10 * oneGbps] [4, 5, 10 * oneGbps] [3, 17, 10 * oneGbps] [5, 17, 10 * oneGbps] [3, 6, 10 * oneGbps] [3, 26, 10 * oneGbps] [3, 14, 10 * oneGbps] [3, 7, 10 * oneGbps] [7, 8, 10 * oneGbps] [8, 9, 10 * oneGbps] [7, 9, 10 * oneGbps] [9, 10, 10 * oneGbps] [6, 26, 10 * oneGbps] [5, 6, 10 * oneGbps] [6, 17, 10 * oneGbps] [6, 14, 10 * oneGbps] [10, 11, 10 * oneGbps] [10, 14, 10 * oneGbps] [11, 14, 10 * oneGbps] [11, 12, 10 * oneGbps] [11, 13, 10 * oneGbps] [12, 13, 10 * oneGbps] [12, 14, 10 * oneGbps] [13, 19, 10 * oneGbps] [14, 19, 10 * oneGbps] [14, 18, 10 * oneGbps] [14, 15, 10 * oneGbps] [14, 16, 10 * oneGbps] [15, 16, 10 * oneGbps] [15, 26, 10 * oneGbps] [16, 26, 10 * oneGbps] [17, 26, 10 * oneGbps] [16, 17, 10 * oneGbps] [16, 18, 10 * oneGbps] [16, 19, 10 * oneGbps] [18, 19, 10 * oneGbps] [19, 20, 10 * oneGbps] [19, 21, 10 * oneGbps] [17, 22, 10 * oneGbps] [17, 23, 10 * oneGbps] [18, 22, 10 * oneGbps] [21, 23, 10 * oneGbps] [22, 23, 10 * oneGbps] [22, 24, 10 * oneGbps] [24, 25, 10 * oneGbps]]'
traffic_matrix = Float64[[8, 4, 6240.1600] [8, 1, 26863.6533] [8, 5, 28226.0267] [8, 3, 47.3956] [8, 1, 2817.4489] [8, 7, 6438.9689] [8, 2, 7656.7911] [8, 6, 4.7111] [8, 9, 1867315.3956] [4, 8, 21866.3111] [4, 1, 4646580.7022] [4, 5, 41628.4000] [4, 3, 59591.7956] [4, 7, 1285829.3600] [4, 2, 1196886.1244] [4, 6, 23.1911] [4, 11, 1926.3911] [4, 9, 63727.9378] [1, 8, 598001.9289] [1, 4, 342504.4533] [1, 5, 14935.0756] [1, 3, 355844.4622] [1, 7, 16094485.4311] [1, 2, 2543737.4400] [1, 6, 3683.2800] [1, 11, 7058.1422] [1, 9, 171177.7600] [5, 8, 35060.1244] [5, 4, 3344.5778] [5, 1, 29604.1867] [5, 3, 86.2400] [5, 10, 5342.3289] [5, 7, 123.6089] [5, 2, 27998.6044] [5, 6, 606.5778] [5, 11, 26.6667] [5, 9, 778276.2578] [3, 8, 27451.4844] [3, 4, 1958030.1956] [3, 1, 843953.6711] [3, 5, 106844.8000] [3, 10, 65148.4178] [3, 7, 369093.9733] [3, 2, 207656.8711] [3, 6, 1749.4133] [3, 11, 24723.9644] [3, 9, 69562.1511] [10, 8, 549654.6933] [10, 1, 63275.7778] [10, 5, 767704.7467] [10, 3, 82480.2756] [10, 7, 95816.9867] [10, 2, 115608.1778] [10, 6, 137330.3022] [10, 11, 81214.3556] [10, 9, 1256926.2222] [7, 8, 41063.8311] [7, 4, 1366652.5244] [7, 1, 2448645.0933] [7, 5, 83908.4711] [7, 3, 310578.6222] [7, 10, 29714.8444] [7, 2, 941842.7822] [7, 6, 8392.0000] [7, 11, 46217.5644] [7, 9, 57574.9067] [2, 8, 43894.6400] [2, 4, 335900.4000] [2, 1, 158223.1644] [2, 5, 16902.1867] [2, 3, 78114.3022] [2, 10, 177692.9600] [2, 7, 74661.2089] [2, 6, 421889.7778] [2, 11, 13222.6756] [2, 9, 18675.6622] [6, 4, 1389.6178] [6, 1, 62036.3289] [6, 5, 193.7067] [6, 3, 2093.6889] [6, 10, 1693.0844] [6, 7, 8646.7733] [6, 2, 364485.7333] [6, 9, 80837.0489] [11, 8, 6725.6444] [11, 4, 4773.3067] [11, 1, 32356.2933] [11, 5, 24303.0578] [11, 3, 53.0489] [11, 10, 9530.7822] [11, 7, 4520.6133] [11, 2, 400782.9511] [11, 6, 0.3556] [11, 9, 248961.2711] [9, 8, 652315.8311] [9, 4, 34186.6578] [9, 1, 21203.9911] [9, 5, 13540.1333] [9, 3, 19.8044] [9, 10, 720.1956] [9, 7, 837.6711] [9, 2, 18402.5333] [9, 11, 9.4578]]'
traffic_matrix[:, 3] /= 1024 * 3 # Mbps
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("att", nodes_matrix, edges_matrix, traffic_matrix))

### Purely synthetic graphs.
## Path graph.
# Oblivious competitive ratio should be 1.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [2, 3, oneGbps] [3, 4, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [2, 4, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("path", nodes_matrix, edges_matrix, traffic_matrix))

## Binary trees.
# Pure tree, seven nodes. Full binary, depth of 3.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [2, 4, oneGbps] [2, 5, oneGbps] [3, 6, oneGbps] [3, 7, oneGbps]]'
traffic_matrix = Float64[[1, 7, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [7, 5, halfGbps] [4, 5, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("full_binary_tree", nodes_matrix, edges_matrix, traffic_matrix))

# Variant. Connect leaves only within the subtree. Same demands as before
# (one demand within a subtree is likely to be impacted, as it can be served
# directly).
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [2, 4, oneGbps] [2, 5, oneGbps] [4, 5, oneGbps] [3, 6, oneGbps] [3, 7, oneGbps] [6, 7, oneGbps]]'
traffic_matrix = Float64[[1, 7, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [7, 5, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("full_binary_tree-connected_leaves_subtree", nodes_matrix, edges_matrix, traffic_matrix))

# Variant. Connect only leaves, be they in the same subtree or not.
# Same demands as before, more changes expected.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [2, 4, oneGbps] [2, 5, oneGbps] [4, 5, oneGbps] [3, 6, oneGbps] [3, 7, oneGbps] [6, 7, oneGbps] [5, 6, oneGbps]]'
traffic_matrix = Float64[[1, 7, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [7, 5, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("full_binary_tree-connected_leaves", nodes_matrix, edges_matrix, traffic_matrix))

# Variant. Connect only the two inner nodes 2 and 3 (depth of 2).
# Same demands as before, many demands can go through the new link.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [2, 3, oneGbps] [2, 4, oneGbps] [2, 5, oneGbps] [3, 6, oneGbps] [3, 7, oneGbps]]'
traffic_matrix = Float64[[1, 7, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [7, 5, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("full_binary_tree-connected_inner", nodes_matrix, edges_matrix, traffic_matrix))

# Variant. Connect only the two inner nodes 2 and 3 (depth of 2) with a
# high-capacity link (impossible to saturate). Same demands as before,
# many demands can go through the new link.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [2, 3, 10 * oneGbps] [2, 4, oneGbps] [2, 5, oneGbps] [3, 6, oneGbps] [3, 7, oneGbps]]'
traffic_matrix = Float64[[1, 7, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [7, 5, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("full_binary_tree-connected_inner-high-capacity", nodes_matrix, edges_matrix, traffic_matrix))

# Variant. Connect only the two inner nodes 2 and 3 (depth of 2) with a
# high-capacity link (impossible to saturate). Same demands as before,
# many demands can go through the new link.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"]], [2, 1])
edges_matrix = [[1, 2, 5 * oneGbps] [1, 3, 7 * oneGbps] [2, 3, 20 * oneGbps] [2, 4, 3 * oneGbps] [2, 5, 2 * oneGbps] [3, 6, 3 * oneGbps] [3, 7, 8 * oneGbps]]'
traffic_matrix = Float64[[1, 7, halfGbps] [2, 6, oneGbps] [6, 4, 2 * oneGbps] [7, 5, oneGbps]]'

push!(topologies, Topology("full_binary_tree-connected_inner-high-capacity-various", nodes_matrix, edges_matrix, traffic_matrix))

## Two "binary trees", connected together.
# Two trees, full binary, depth of 3. Connected only at their root (i.e.
# full binary tree of depth 4, in this case). Symmetric demands for the two
# graphs, plus a few crossing the trees. First tree rooted at 1, second at 8.
# Inspired by https://epubs.siam.org/doi/abs/10.1137/S0895479896312262.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"] [13, "M"] [14, "N"]], [2, 1])
edges_matrix = [[1, 2, 1.5 * oneGbps] [1, 3, 1.5 * oneGbps] [2, 4, 1.5 * oneGbps] [2, 5, 1.5 * oneGbps] [3, 6, 1.5 * oneGbps] [3, 7, 1.5 * oneGbps] [8, 9, 1.5 * oneGbps] [8, 10, 1.5 * oneGbps] [9, 11, 1.5 * oneGbps] [9, 12, 1.5 * oneGbps] [10, 13, 1.5 * oneGbps] [10, 14, 1.5 * oneGbps] [1, 8, 1.5 * oneGbps]]'
traffic_matrix = Float64[[1, 7, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [7, 5, halfGbps] [4, 5, halfGbps] [8, 14, halfGbps] [9, 13, halfGbps] [13, 11, halfGbps] [14, 12, halfGbps] [11, 12, halfGbps] [2, 9, halfGbps] [10, 2, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("double_tree", nodes_matrix, edges_matrix, traffic_matrix))

# Variant. Two trees, full binary, depth of 3. The corresponding nodes are fully
# connected. Symmetric demands for the two graphs, plus a few crossing the trees
# (unchanged from 11). First tree rooted at 1, second at 8.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"] [13, "M"] [14, "N"]], [2, 1])
edges_matrix = [[1, 2, 1.5 * oneGbps] [1, 3, 1.5 * oneGbps] [2, 4, 1.5 * oneGbps] [2, 5, 1.5 * oneGbps] [3, 6, 1.5 * oneGbps] [3, 7, 1.5 * oneGbps] [8, 9, 1.5 * oneGbps] [8, 10, 1.5 * oneGbps] [9, 11, 1.5 * oneGbps] [9, 12, 1.5 * oneGbps] [10, 13, 1.5 * oneGbps] [10, 14, 1.5 * oneGbps] [1, 8, 1.5 * oneGbps] [2, 9, 1.5 * oneGbps] [3, 10, 1.5 * oneGbps] [4, 11, 1.5 * oneGbps] [5, 12, 1.5 * oneGbps] [6, 13, 1.5 * oneGbps] [7, 14, 1.5 * oneGbps]]'
traffic_matrix = Float64[[1, 7, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [7, 5, halfGbps] [4, 5, halfGbps] [8, 14, halfGbps] [9, 13, halfGbps] [13, 11, halfGbps] [14, 12, halfGbps] [11, 12, halfGbps] [2, 9, halfGbps] [10, 2, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("double_tree-fully_connected", nodes_matrix, edges_matrix, traffic_matrix))

# Variant. Two trees, full binary, depth of 3. The corresponding nodes are fully
# connected. Symmetric demands for the two graphs, plus a few crossing the trees
# (unchanged from 11). First tree rooted at 1, second at 8.
if INCLUDE_OVERLOADED_TOPOLOGIES
    nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"] [13, "M"] [14, "N"]], [2, 1])
    edges_matrix = [[1, 2, 14.24 * oneGbps] [1, 3, .01 * oneGbps] [2, 4, pi * oneGbps] [2, 5, 10 * oneGbps] [3, 6, oneGbps / pi] [3, 7, .001 * oneGbps] [8, 9, 10 * oneGbps] [8, 10, .01 * oneGbps] [9, 11, oneGbps] [9, 12, 10 * oneGbps] [10, 13, oneGbps] [10, 14, oneGbps] [1, 8, 10 * oneGbps] [2, 9, oneGbps] [3, 10, oneGbps] [4, 11, 10 * oneGbps] [5, 12, oneGbps] [6, 13, oneGbps] [7, 14, 10 * oneGbps]]'
    traffic_matrix = Float64[[1, 7, halfGbps] [2, 6, .1 * halfGbps] [6, 4, halfGbps] [7, 5, 2.4 * halfGbps] [4, 5, halfGbps] [8, 14, 5 * halfGbps] [9, 13, halfGbps] [13, 11, halfGbps] [14, 12, 1.85 * halfGbps] [11, 12, halfGbps] [2, 9, halfGbps] [10, 2, halfGbps]]'
    traffic_matrix[:, 3] .*= traffic_scaling

    push!(topologies, Topology("double_tree-fully_connected-various", nodes_matrix, edges_matrix, traffic_matrix))
end

# Variant. Two trees, full binary, depth of 3. The corresponding nodes are
# connected (a few edges do not exist, the ones that correspond to demands:
# 2-9 for demand 2-9, 3-10 for demand 10-2). Symmetric demands for the
# two graphs, plus a few crossing the trees (unchanged from 11). First tree
# rooted at 1, second at 8.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"] [13, "M"] [14, "N"]], [2, 1])
edges_matrix = [[1, 2, 1.5 * oneGbps] [1, 3, 1.5 * oneGbps] [2, 4, 1.5 * oneGbps] [2, 5, 1.5 * oneGbps] [3, 6, 1.5 * oneGbps] [3, 7, 1.5 * oneGbps] [8, 9, 1.5 * oneGbps] [8, 10, 1.5 * oneGbps] [9, 11, 1.5 * oneGbps] [9, 12, 1.5 * oneGbps] [10, 13, 1.5 * oneGbps] [10, 14, 1.5 * oneGbps] [1, 8, 1.5 * oneGbps] [4, 11, 1.5 * oneGbps] [5, 12, 1.5 * oneGbps] [6, 13, 1.5 * oneGbps] [7, 14, 1.5 * oneGbps]]'
traffic_matrix = Float64[[1, 7, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [7, 5, halfGbps] [4, 5, halfGbps] [8, 14, halfGbps] [9, 13, halfGbps] [13, 11, halfGbps] [14, 12, halfGbps] [11, 12, halfGbps] [2, 9, halfGbps] [10, 2, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("double_tree-connected", nodes_matrix, edges_matrix, traffic_matrix))

# Variant. Two trees, full binary, depth of 3. The corresponding nodes are
# connected (a few edges do not exist, the ones that correspond to demands:
# 2-9 for demand 2-9, 3-10 for demand 10-2). Symmetric demands for the two
# graphs, plus a few crossing the trees (unchanged from 11). First tree rooted
# at 1, second at 8.
if INCLUDE_OVERLOADED_TOPOLOGIES
    nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"] [13, "M"] [14, "N"]], [2, 1])
    edges_matrix = [[1, 2, 10 * oneGbps] [1, 3, 5 * oneGbps] [2, 4, oneGbps] [2, 5, 10 * oneGbps] [3, 6, .05 * oneGbps] [3, 7, oneGbps] [8, 9, 10 * oneGbps] [8, 10, 5 * oneGbps] [9, 11, oneGbps] [9, 12, 10 * oneGbps] [10, 13, 5 * oneGbps] [10, 14, oneGbps] [1, 8, 10 * oneGbps] [4, 11, 5 * oneGbps] [5, 12, oneGbps] [6, 13, 10 * oneGbps] [7, 14, 5 * oneGbps]]'
    traffic_matrix = Float64[[1, 7, 5 * halfGbps] [2, 6, halfGbps] [6, 4, 5.2 * halfGbps] [7, 5, halfGbps] [4, 5, .5 * halfGbps] [8, 14, halfGbps] [9, 13, 5 * halfGbps] [13, 11, halfGbps] [14, 12, 5 * halfGbps] [11, 12, halfGbps] [2, 9, 5 * halfGbps] [10, 2, halfGbps]]'
    traffic_matrix[:, 3] .*= traffic_scaling

    push!(topologies, Topology("double_tree-connected-various", nodes_matrix, edges_matrix, traffic_matrix))
end

# Variant. Two trees, full binary, depth of 3. The corresponding nodes are
# connected (a few edges do not exist, the ones that correspond to demands:
# 2-9 for demand 2-9, 3-10 for demand 10-2). Symmetric demands for the two
# graphs, plus a few crossing the trees (unchanged from 11). First tree rooted
# at 1, second at 8. Much lower capacities than usual.
if INCLUDE_OVERLOADED_TOPOLOGIES
    nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"] [13, "M"] [14, "N"]], [2, 1])
    edges_matrix = [[1, 2, 10 * oneGbps] [1, 3, 5 * oneGbps] [2, 4, oneGbps] [2, 5, 10 * oneGbps] [3, 6, .005 * oneGbps] [3, 7, oneGbps] [8, 9, 10 * oneGbps] [8, 10, 5 * oneGbps] [9, 11, oneGbps] [9, 12, 10 * oneGbps] [10, 13, 5 * oneGbps] [10, 14, oneGbps] [1, 8, 10 * oneGbps] [4, 11, 5 * oneGbps] [5, 12, oneGbps] [6, 13, 10 * oneGbps] [7, 14, 5 * oneGbps]]'
    traffic_matrix = Float64[[1, 7, 5 * halfGbps] [2, 6, halfGbps] [6, 4, 5.2 * halfGbps] [7, 5, halfGbps] [4, 5, .5 * halfGbps] [8, 14, halfGbps] [9, 13, 5 * halfGbps] [13, 11, halfGbps] [14, 12, 5 * halfGbps] [11, 12, halfGbps] [2, 9, 5 * halfGbps] [10, 2, halfGbps]]'
    traffic_matrix[:, 3] .*= traffic_scaling

    push!(topologies, Topology("double_tree-connected-various_very_low", nodes_matrix, edges_matrix, traffic_matrix))
end

# Variant. Two trees, full binary, depth of 3. The corresponding nodes are
# connected (a few edges do not exist, the ones that correspond to demands:
# 2-9 for demand 2-9, 3-10 for demand 10-2). Symmetric demands for the two
# graphs, plus a few crossing the trees (unchanged from 11). First tree rooted
# at 1, second at 8. Much lower capacities than usual.
if INCLUDE_OVERLOADED_TOPOLOGIES
    nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"] [13, "M"] [14, "N"]], [2, 1])
    edges_matrix = [[1, 2, 10 * oneGbps] [1, 3, 5 * oneGbps] [2, 4, oneGbps] [2, 5, 10 * oneGbps] [3, 6, .0005 * oneGbps] [3, 7, oneGbps] [8, 9, 10 * oneGbps] [8, 10, 5 * oneGbps] [9, 11, oneGbps] [9, 12, 10 * oneGbps] [10, 13, 5 * oneGbps] [10, 14, oneGbps] [1, 8, 10 * oneGbps] [4, 11, 5 * oneGbps] [5, 12, oneGbps] [6, 13, 10 * oneGbps] [7, 14, 5 * oneGbps]]'
    traffic_matrix = Float64[[1, 7, 5 * halfGbps] [2, 6, halfGbps] [6, 4, 5.2 * halfGbps] [7, 5, halfGbps] [4, 5, .5 * halfGbps] [8, 14, halfGbps] [9, 13, 5 * halfGbps] [13, 11, halfGbps] [14, 12, 5 * halfGbps] [11, 12, halfGbps] [2, 9, 5 * halfGbps] [10, 2, halfGbps]]'
    traffic_matrix[:, 3] .*= traffic_scaling

    push!(topologies, Topology("double_tree-connected-various_very_very_low", nodes_matrix, edges_matrix, traffic_matrix))
end

## Clique of 6 nodes.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [1, 4, oneGbps] [1, 5, oneGbps] [1, 6, oneGbps] [2, 3, oneGbps] [2, 4, oneGbps] [2, 5, oneGbps] [2, 6, oneGbps] [3, 4, oneGbps] [3, 5, oneGbps] [3, 6, oneGbps] [4, 5, oneGbps] [4, 6, oneGbps] [5, 6, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [1, 5, halfGbps] [4, 5, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("clique", nodes_matrix, edges_matrix, traffic_matrix))

## Two clique of 6 nodes.
# Two cliques of 6 nodes, connected once. Symmetrical demands, plus a few that
# cross the connecting link.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [1, 4, oneGbps] [1, 5, oneGbps] [1, 6, oneGbps] [2, 3, oneGbps] [2, 4, oneGbps] [2, 5, oneGbps] [2, 6, oneGbps] [3, 4, oneGbps] [3, 5, oneGbps] [3, 6, oneGbps] [4, 5, oneGbps] [4, 6, oneGbps] [5, 6, oneGbps] [7, 8, oneGbps] [7, 9, oneGbps] [7, 10, oneGbps] [7, 11, oneGbps] [7, 12, oneGbps] [8, 9, oneGbps] [8, 10, oneGbps] [8, 11, oneGbps] [8, 12, oneGbps] [9, 10, oneGbps] [9, 11, oneGbps] [9, 12, oneGbps] [10, 11, oneGbps] [10, 12, oneGbps] [11, 12, oneGbps] [1, 7, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [1, 5, halfGbps] [4, 5, halfGbps] [7, 9, halfGbps] [8, 12, halfGbps] [12, 10, halfGbps] [7, 11, halfGbps] [10, 11, halfGbps] [1, 7, halfGbps] [2, 8, halfGbps] [9, 3, halfGbps] [10, 4, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("double_clique", nodes_matrix, edges_matrix, traffic_matrix))

# Two cliques of 6 nodes, connected 4 times (4 end-points in the first clique,
# 1 end-point in the second).
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [1, 4, oneGbps] [1, 5, oneGbps] [1, 6, oneGbps] [2, 3, oneGbps] [2, 4, oneGbps] [2, 5, oneGbps] [2, 6, oneGbps] [3, 4, oneGbps] [3, 5, oneGbps] [3, 6, oneGbps] [4, 5, oneGbps] [4, 6, oneGbps] [5, 6, oneGbps] [7, 8, oneGbps] [7, 9, oneGbps] [7, 10, oneGbps] [7, 11, oneGbps] [7, 12, oneGbps] [8, 9, oneGbps] [8, 10, oneGbps] [8, 11, oneGbps] [8, 12, oneGbps] [9, 10, oneGbps] [9, 11, oneGbps] [9, 12, oneGbps] [10, 11, oneGbps] [10, 12, oneGbps] [11, 12, oneGbps] [1, 7, oneGbps] [3, 7, oneGbps] [4, 7, oneGbps] [5, 7, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [1, 5, halfGbps] [4, 5, halfGbps] [7, 9, halfGbps] [8, 12, halfGbps] [12, 10, halfGbps] [7, 11, halfGbps] [10, 11, halfGbps] [1, 7, halfGbps] [2, 8, halfGbps] [9, 3, halfGbps] [10, 4, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("double_clique-connected", nodes_matrix, edges_matrix, traffic_matrix))

## Star of 6 nodes.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [1, 4, oneGbps] [1, 5, oneGbps] [1, 6, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [1, 5, halfGbps] [4, 5, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("star", nodes_matrix, edges_matrix, traffic_matrix))

## Two stars of 6 nodes.
# Two stars of 6 nodes, connected once. Symmetrical demands, plus a few that
# cross the connecting link.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [1, 4, oneGbps] [1, 5, oneGbps] [1, 6, oneGbps] [7, 8, oneGbps] [7, 9, oneGbps] [7, 10, oneGbps] [7, 11, oneGbps] [7, 12, oneGbps] [1, 7, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [1, 5, halfGbps] [4, 5, halfGbps] [7, 9, halfGbps] [8, 12, halfGbps] [12, 10, halfGbps] [7, 11, halfGbps] [10, 11, halfGbps] [1, 7, halfGbps] [2, 8, halfGbps] [9, 3, halfGbps] [10, 4, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("double_star", nodes_matrix, edges_matrix, traffic_matrix))

# Two stars of 6 nodes, connected 4 times (4 end-points in the first star,
# 1 end-point in the star; never the centre of the star).
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [1, 4, oneGbps] [1, 5, oneGbps] [1, 6, oneGbps] [7, 8, oneGbps] [7, 9, oneGbps] [7, 10, oneGbps] [7, 11, oneGbps] [7, 12, oneGbps] [1, 7, oneGbps] [3, 7, oneGbps] [4, 7, oneGbps] [5, 7, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [1, 5, halfGbps] [4, 5, halfGbps] [7, 9, halfGbps] [8, 12, halfGbps] [12, 10, halfGbps] [7, 11, halfGbps] [10, 11, halfGbps] [1, 7, halfGbps] [2, 8, halfGbps] [9, 3, halfGbps] [10, 4, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("double_star-connected", nodes_matrix, edges_matrix, traffic_matrix))

## Grid of 12 nodes.
# 3 rows and 4 columns. Edges along rows or columns (no diagonals).
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "G"] [8, "H"] [9, "I"] [10, "J"] [11, "K"] [12, "L"]], [2, 1])
edges_matrix = [[1, 2, 1.5 * oneGbps] [2, 3, 1.5 * oneGbps] [3, 4, 1.5 * oneGbps] [1, 5, 1.5 * oneGbps] [5, 6, 1.5 * oneGbps] [6, 7, 1.5 * oneGbps] [7, 8, 1.5 * oneGbps] [2, 6, 1.5 * oneGbps] [3, 7, 1.5 * oneGbps] [4, 8, 1.5 * oneGbps] [5, 9, 1.5 * oneGbps] [9, 10, 1.5 * oneGbps] [10, 11, 1.5 * oneGbps] [11, 12, 1.5 * oneGbps] [6, 10, 1.5 * oneGbps] [7, 11, 1.5 * oneGbps] [8, 12, 1.5 * oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [2, 6, halfGbps] [6, 4, halfGbps] [1, 5, halfGbps] [4, 5, halfGbps] [7, 9, halfGbps] [8, 12, halfGbps] [12, 10, halfGbps] [7, 11, halfGbps] [10, 11, halfGbps] [1, 7, halfGbps] [2, 8, halfGbps] [9, 3, halfGbps] [10, 4, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("grid", nodes_matrix, edges_matrix, traffic_matrix))

## Debugging.
# Very simple topologies, only for test purposes.
nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [2, 3, oneGbps] [3, 1, oneGbps]]'
traffic_matrix = Float64[[1, 3, halfGbps] [2, 3, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("debug3", nodes_matrix, edges_matrix, traffic_matrix))

nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [2, 3, oneGbps] [2, 1, oneGbps] [1, 3, oneGbps]]' # Supposed to be used with exactly these orientations.
traffic_matrix = Float64[[1, 3, halfGbps] [2, 3, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("debug4", nodes_matrix, edges_matrix, traffic_matrix))

nodes_matrix = permutedims([[1, "A"] [2, "B"]], [2, 1])
edges_matrix = zeros(Int, 1, 3)
edges_matrix[1, :] = [1, 2, oneGbps]
traffic_matrix = zeros(Float64, 1, 3)
traffic_matrix[1, :] = [1, 2, halfGbps]
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("debug", nodes_matrix, edges_matrix, traffic_matrix))

nodes_matrix = permutedims([[1, "A"] [2, "B"] [3, "C"] [4, "D"] [5, "E"] [6, "F"] [7, "H"] [8, "I"] [9, "J"] [10, "K"] [11, "L"] [12, "M"]], [2, 1])
edges_matrix = [[1, 2, oneGbps] [1, 3, oneGbps] [2, 4, oneGbps] [2, 5, oneGbps] [3, 6, oneGbps] [7, 8, oneGbps] [7, 9, oneGbps] [8, 10, oneGbps] [8, 11, oneGbps] [9, 12, oneGbps] [1, 7, oneGbps] [4, 10, oneGbps] [5, 11, oneGbps] [6, 12, oneGbps]]'
traffic_matrix = Float64[[2, 8, halfGbps] [9, 2, halfGbps]]'
traffic_matrix[:, 3] .*= traffic_scaling

push!(topologies, Topology("debug_double_tree-connected", nodes_matrix, edges_matrix, traffic_matrix))

### Finally, generate a list of names
names = sort(collect(t.name for t in topologies))


text:
cutlnn(x) = x
function cutlnn(ex::Expr)
    Expr(ex.head, ex.args[[!(arg isa LineNumberNode) for arg in ex.args]]...)
end

cutblock(x) = x
function cutblock(ex::Expr)
    if (ex.head == :block) && (length(ex.args)==1)
        ex.args[1]
    else
        ex
    end
end

tupleargs(x::Symbol) = (x,)
function tupleargs(ex::Expr)
    if ex.head == :tuple
        Tuple(ex.args)
    else
        throw("expression must have head `:tuple`")
    end
end

pretty_kwargs(x::Symbol) = x
function pretty_kwargs(ex)
    if ex.args[1] isa Expr
        if ex.args[1].head == :block
            exout = copy(ex)
            exout.args[1].head = :tuple
            exout
        else
            ex
        end
    else
        ex
    end
end


text:
facts("Population") do

  context("FloatVectorPopulation constructor") do
    p1 = FloatVectorPopulation(10, 2)
    @fact typeof(p1) => FloatVectorPopulation

    @fact size(p1.individuals, 1) => 10
    @fact size(p1.individuals, 2) => 2

    @fact size(p1.fitness, 1) => 10
    @fact size(p1.fitness, 2) => 2 # This is a bug since we cannot now the size of the fitness vectors before we have evaluated one...

    @fact size(p1.top, 1) => 10
    @fact size(p1.top, 2) => 2

    @fact size(p1.top_fitness, 1) => 10
    @fact size(p1.top_fitness, 2) => 2 # This is a bug since we cannot now the size of the fitness vectors before we have evaluated one...
  end

end


text:
# This file defines the History type. This type holds all the changes in a
#   specific attribute of a specific person.
# The history also contains the attribute of which it is the history, for
#   validation purposes.

# The History type requires the HistoryEntry type.
requiredTypes = [ ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export History
type History
    # The attribute of which this is the history.
    attribute::Symbol

    # The list of events.
    history::Vector{HistoryEntry}

    function History( attr::Symbol )
        newHist = new()
        newHist.attribute = attr
        newHist.history = Vector{HistoryEntry}()
        return newHist
    end  # History( attr )
end  # type History


text:
## test_fem.jl : test GRF defined on FEM

@testset "unstructured grid        " begin

    ## Generate a set of random points ##
    pts = rand(512, 2)
    m = Matern(0.2,2.0)

    ## Cholesky ##
    grf = GaussianRandomField(CovarianceFunction(2,m),Cholesky(),pts)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Matern)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{Cholesky})
    @test length(grf.pts) == 2
    @test size(grf.pts[1],1) == 2
    @test size(grf.pts[1],2) == size(pts,1)

    ## Spectral ##
    grf = GaussianRandomField(CovarianceFunction(2,m),Spectral(),pts)
    @test isa(grf,GaussianRandomField{Spectral})

    ## CirculantEmbedding ##
    @test_throws ArgumentError GaussianRandomField(CovarianceFunction(2,m),CirculantEmbedding(),pts)

    ## KarhunenLoeve ##
    grf = GaussianRandomField(CovarianceFunction(2,m),KarhunenLoeve(128),pts)
    @test isa(grf,GaussianRandomField{KarhunenLoeve{128}})

    ## Anisotropic random field ##
    a = AnisotropicExponential([1000 0; 0 1000])
    cov = CovarianceFunction(2,a)
    grf = GaussianRandomField(cov,KarhunenLoeve(500),pts,quad=GaussLegendre())
    @test isa(grf.cov.cov,AnisotropicExponential)

    ## GaussianRandomField with mean value ##
    grf = GaussianRandomField(1,CovarianceFunction(2,m),Cholesky(),pts)
    @test isa(grf,GaussianRandomField{Cholesky})
    @test length(grf.mean) == size(pts,1)
    @test all(grf.mean.==1)

    grf = GaussianRandomField(ones(size(pts,1)),CovarianceFunction(2,m),Spectral(),pts)
    @test isa(grf,GaussianRandomField{Spectral})
    @test length(grf.mean) == size(pts,1)
    @test all(grf.mean.==1)

    ## Test spherical domain ##
    pts = 4*rand(1024, 3) .- 2
    d = pts[:,1].*pts[:,1] + pts[:,2].*pts[:,2] + pts[:,3].*pts[:,3]
    pts = pts[d .< 4, :]
    cov = CovarianceFunction(3, Matern(1, 3/2))
    grf = GaussianRandomField(cov, KarhunenLoeve(128), pts)
    z = sample(grf)
    @test isa(grf,GaussianRandomField{KarhunenLoeve{128}})
    grf = GaussianRandomField(cov, Spectral(), pts)
    @test isa(grf,GaussianRandomField{Spectral})
    grf = GaussianRandomField(cov, Cholesky(), pts)
    @test isa(grf,GaussianRandomField{Cholesky})

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

module PointPatterns

using Meshes
using GeoStatsBase

using Distributions

import Random

include("processes.jl")
include("thinning.jl")

export
  # point processes
  PointProcess,
  BinomialProcess,
  PoissonProcess,
  UnionProcess,
  ishomogeneous,

  # thinning methods
  AbstractThinning,
  RandomThinning,
  thin

end # module


text:
const _keep_fields = [
    :sample,
    :subject,
    :timepoint,
    :ECHOTPCoded,
    :Mother_Child,
    :MaternaID,
    :Fecal_EtOH,
    :CollectionRep,
    :sid_old,
    :CovidCollectionNumber,
    :childAgeMonths
]

"""
    airtable_metadata(key=ENV["AIRTABLE_KEY"])
Get fecal sample metadata table from airtable.
The API `key` comes from https://airtable.com/account.
This is unlikely to work if you're not in the VKC lab,
but published sample metadata is available from OSF.io
using `datadep"sample metadata"`.
"""
function airtable_metadata()
    base = AirBase("appSWOVVdqAi5aT5u")

    samples      = Airtable.query(AirTable("Samples", base))
    mgxbatches   = Airtable.query(AirTable("MGX Batches", base))
    metabbatches = Airtable.query(AirTable("Metabolomics Batches", base))

    df = DataFrame()

    @progress "Adding samples" for samp in samples
        # skip anything not in the "ECHO" project
        "rec8zdkF1pchFNiJC" ∈ get(samp, :Project, []) || continue

        mgxids = string.(get(samp, Symbol("MGX Batches"), []))
        mgxbatch = findall(batch-> Airtable.id(batch) ∈ mgxids, mgxbatches)
        mgxbatch = isempty(mgxbatch) ? missing : mgxbatches[first(mgxbatch)][:Name]

        metabids = string.(get(samp, Symbol("Metabolomics Batches"), []))
        metabbatch = findall(batch-> Airtable.id(batch) ∈ metabids, metabbatches)
        metabbatch = isempty(metabbatch) ? missing : metabbatches[first(metabbatch)][:Name]

        rec =  (; :airtable_id        => Airtable.id(samp),
                  (kf=> get(samp, kf, missing) for kf in _keep_fields)...,
                  :Mgx_batch          => mgxbatch,
                  :Metabolomics_batch => metabbatch,
                  )
        push!(df, rec, cols=:union)
    end

    subset!(df, AsTable(["subject", "timepoint"])=> ByRow(r-> !any(ismissing, r)))

    transform!(df, "subject"   => ByRow(s-> parse(Int, s)) => "subject",
                   "timepoint" => ByRow(tp-> parse(Int, tp)) => "timepoint")
    return select(df, Cols(:sample, :subject, :timepoint, :))
end

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis


text:
function OpenIntToDoubleHashMap()
    return OpenIntToDoubleHashMap(())
end

function OpenIntToDoubleHashMap(arg0::OpenIntToDoubleHashMap)
    return OpenIntToDoubleHashMap((OpenIntToDoubleHashMap,), arg0)
end

function OpenIntToDoubleHashMap(arg0::jdouble)
    return OpenIntToDoubleHashMap((jdouble,), arg0)
end

function OpenIntToDoubleHashMap(arg0::jint)
    return OpenIntToDoubleHashMap((jint,), arg0)
end

function OpenIntToDoubleHashMap(arg0::jint, arg1::jdouble)
    return OpenIntToDoubleHashMap((jint, jdouble), arg0, arg1)
end

function contains_key(obj::OpenIntToDoubleHashMap, arg0::jint)
    return jcall(obj, "containsKey", jboolean, (jint,), arg0)
end

function get(obj::OpenIntToDoubleHashMap, arg0::jint)
    return jcall(obj, "get", jdouble, (jint,), arg0)
end

function iterator(obj::OpenIntToDoubleHashMap)
    return jcall(obj, "iterator", OpenIntToDoubleHashMap_JIterator, ())
end

function put(obj::OpenIntToDoubleHashMap, arg0::jint, arg1::jdouble)
    return jcall(obj, "put", jdouble, (jint, jdouble), arg0, arg1)
end

function remove(obj::OpenIntToDoubleHashMap, arg0::jint)
    return jcall(obj, "remove", jdouble, (jint,), arg0)
end

function size(obj::OpenIntToDoubleHashMap)
    return jcall(obj, "size", jint, ())
end


text:
import Base: getindex, push!, append!, show, length, size
export SpiceIntCell, SpiceDoubleCell, SpiceCharCell, appnd, push!, append!, card, length, copy, size_c

const CTRLSZ = 6

dtype(::Type{SpiceChar}) = 0
dtype(::Type{SpiceDouble}) = 1
dtype(::Type{SpiceInt}) = 2

mutable struct Cell{T}
    # Do not reorder these fields
    dtype::SpiceInt
    length::SpiceInt
    size::SpiceInt
    card::SpiceInt
    isset::SpiceInt
    adjust::SpiceInt
    init::SpiceInt
    base::Ptr{Cvoid}
    data::Ptr{Cvoid}
    function Cell{T}(length, size) where T
        new{T}(dtype(T), length, size, 0, 1, 0, 0)
    end
end

itertype(::Type{T}) where {T} = T
itertype(::Type{SpiceChar}) = String

init_data(::Type{T}, size, _) where {T} = Array{T}(undef, CTRLSZ + size)
init_data(::Type{SpiceChar}, size, length) = fill(SpiceChar(0), (length, CTRLSZ + size))

len(::Type{T}, _) where {T} = 0
len(::Type{SpiceChar}, length) = length += 1

dim(::Type{T}) where {T} = 1
dim(::Type{SpiceChar}) = 2

data_ptr(::Type{SpiceChar}, data, length) = pointer(data, CTRLSZ * length + 1)
data_ptr(::Type{T}, data, _) where {T} = pointer(data, CTRLSZ + 1)

mutable struct SpiceCell{S, T, N} <: AbstractArray{T, 1}
    data::Array{S, N}
    cell::Cell{S}

    function SpiceCell{T}(size::Integer, length::Integer=256) where T
        length = len(T, length)
        cell = Cell{T}(length, size)
        data = init_data(T, size, length)
        self = new{T, itertype(T), dim(T)}(data, cell)
        self.cell.base = pointer(self.data, 1)
        self.cell.data = data_ptr(T, self.data, length)
        self
    end
end

const SpiceDoubleCell = SpiceCell{SpiceDouble}
const SpiceIntCell = SpiceCell{SpiceInt}
const SpiceCharCell = SpiceCell{SpiceChar}

Base.IndexStyle(::SpiceCell) = IndexLinear()
Base.firstindex(cell::SpiceCell) = 1
Base.lastindex(cell::SpiceCell) = cell.cell.card

"""
    SpiceCharCell(size, length)

Create a SpiceCharCell that can contain up to `size` strings with `length` characters.
"""
SpiceCharCell

function show(io::IO, cell::SpiceCell{T}) where T
    print(io, "SpiceCell{$(T.name.name)}($(cell.cell.size))")
end
function show(io::IO, cell::SpiceCharCell)
    print(io, "SpiceCell{SpiceChar}($(cell.cell.size),$(cell.cell.length))")
end

function check_ind(cell, inds)
    ind = collect(inds)[end]
    if ind > cell.cell.size
        error("Index $ind is out of bounds. The size of this SpiceCell is $(cell.cell.size).")
    elseif ind > cell.cell.card
        error("The element at index $ind is unassigned.")
    end
end

function getindex(cell::SpiceCell, ind::Int)
    check_ind(cell, ind)
    cell.data[CTRLSZ .+ ind]
end

function getindex(cell::SpiceCharCell, ind::Int)
    check_ind(cell, ind)
    str = cell.data[:, CTRLSZ .+ ind]
    chararray_to_string(str)
end

"""
    SpiceDoubleCell(size)

Create a SpiceDoubleCell that can contain up to `size` elements.
"""
SpiceDoubleCell

"""
    SpiceIntCell(size)

Create a SpiceIntCell that can contain up to `size` elements.
"""
SpiceIntCell

for (t, f) in zip((SpiceInt, SpiceDouble), ("appndi_c", "appndd_c"))
    @eval begin
        function appnd(item, cell::SpiceCell{$t, $t})
            data = cell.data
            GC.@preserve data begin
                ccall(($f, libcspice), Cvoid, ($t, Ref{Cell{$t}}), item, cell.cell)
            end
            handleerror()
            return nothing
        end
    end
end

"""
    appnd(item, cell)

Append an `item` to the char/double/integer SpiceCell `cell`.

### References ###

- [appndc - NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/appndc_c.html)
- [appndd - NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/appndd_c.html)
- [appndi - NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/appndi_c.html)
"""
function appnd(item, cell::SpiceCharCell)
    data = cell.data
    GC.@preserve data begin
        ccall((:appndc_c, libcspice), Cvoid, (Cstring, Ref{Cell{SpiceChar}}), item, cell.cell)
    end
    handleerror()
    return nothing
end

"""
    push!(cell, items...)

Insert one or more `items` at the end of the char/double/integer SpiceCell `cell`.
"""
push!(cell::SpiceCell, items...) = foreach(x -> appnd(x, cell), items)

"""
    append!(cell, collection)

Append all items from `collection` to the char/double/integer SpiceCell `cell`.
"""
append!(cell::SpiceCell, collection) = foreach(x -> appnd(x, cell), collection)

"""
    card(cell)

Returns the cardinality (number of elements) of `cell`.
"""
card(cell::SpiceCell) = Int(cell.cell.card)

"""
    length(cell)

Returns the cardinality (number of elements) of `cell`.
"""
length(cell::SpiceCell) = card(cell)

"""
    size(cell)

Returns the cardinality (number of elements) of `cell`.
"""
size(cell::SpiceCell) = (card(cell),)

"""
    copy(cell::SpiceCell)

Duplicate the SpiceCell `cell`.
"""
function Base.copy(cell::SpiceCell{T}) where {T}
    cell_copy = deepcopy(cell)
    cell_copy.cell.base = pointer(cell_copy.data, 1)
    cell_copy.cell.data = data_ptr(T, cell_copy.data, length)
    cell_copy
end

"""
    size_c(cell::SpiceCell)

Returns the maximum number of elements that `cell` can hold.
"""
size_c(cell) = cell.cell.size


text:
# Ex 3.1
println("First number: ")
num1 = parse(Int, chomp(readline()))
println("Second number: ")
num2 = parse(Int, chomp(readline()))

println(num1, " ", num2, " Average: ", Int(floor( (num1 + num2) / 2 )))

# Ex 3.2
# echo "54" > nums_file
# echo "72" >> nums_file
# cat nums_file | julia ex3.jl

# Ex 3.3

function count_negative(ln)
    count = 0
    match = findfirst(isequal('-'), ln)
    while match !== nothing
        count += 1
        match = findnext(isequal('-'), ln, match + 1)
    end
    count
end

num_neg = open("ex1.dat") do f
    num_neg = 0
    for ln in eachline(f)
        num_neg += count_negative(ln)
    end
    num_neg
end

# Ex 3.4
temperature = chomp(readline())
temp_num = parse(Int, temperature[1:end-1])
temp_scale = temperature[end]

if temp_scale == 'F'
    println((temp_num-32)/1.8, "C")
elseif temp_scale == 'C'
    println(temp_num*1.8 + 32, "F")
end

# Ex 3.5
function extract_acc(ln)
    last_index = findfirst(isequal(' '), ln)
    ln[2:last_index-1]
end

accs = open("/home/frederik/julia_course/data/orphans.sp") do f
    accs = String[]
    for ln in eachline(f)
        if ln[1] == '>'
            push!(accs, extract_acc(ln))
        end
    end
    accs
end

println(accs)

# Ex 3.6

# Ex 3.7-8
function read_dna(f)
    dna = ""
    for ln in eachline(f)
        dna = string(dna, rstrip(ln))
    end
    dna
end

in_dna = open(read_dna, "/home/frederik/julia_course/data/dna.dat")

function complement(str)
    trans = Dict('A' => 'T', 'T' => 'A', 'C' => 'G', 'G' => 'C')
    trans_str = map(x -> trans[x], str)
end

function reverse_complement(str)
    reverse(complement(str))
end

println(in_dna)
println(reverse_complement(in_dna))

# Ex 3.9
function to_outfile(seq)
    outfile = open("revdna.dat", "w")

    for i in 1:60:length(seq)
        if i+59 > length(seq)
            println(outfile, seq[i:end])
        else
            println(outfile, seq[i:(i+59)])
        end
    end

    close(outfile)
end

to_outfile(reverse_complement(in_dna))

# Ex 3.10
function read_fasta(f)
    headlines = String[]
    seqs = String[]
    seq = ""
    for line = eachline(f)
        ln = chomp(line)
        if ln[1] == '>'
            if seq != ""
                push!(seqs, seq)
                seq = ""
            end
            push!(headlines, ln)
        else
            seq *= ln
        end
    end
    push!(seqs, seq)
    headlines, seqs
end

open("data/revdna.fsa", "w") do outfile
    headlines, seqs = open(read_fasta, "data/dna.fsa")
    for (headline, seq) in zip(headlines, seqs)
        println(outfile, "$headline ComplementStrand")
        revseq = reverse_complement(seq)
        for i in 1:60:length(revseq)
            if i + 59 > length(revseq)
                println(outfile, revseq[i:end])
            else
                println(outfile, revseq[i:(i + 59)])
            end
        end
    end
end

# Exercise 11
function count_ATCG(seq)
    AT = 0
    GC = 0
    for char in seq
        if char in ('A', 'T')
            AT += 1
        elseif char in ('G', 'C')
            GC += 1
        end
    end
    return string("AT/GC content: ", AT/GC)
end

headlines, seqs = open(read_fasta, "data/dna.fsa")
println(count_ATCG(seqs[1]))

# Exercise 12
function bullseye(dimensions = (40, 40), c1 = "*", c2 = "-", c3 = " ")
    xdim, ydim = dimensions
    for x in 1:xdim
        for y in 1:ydim
            _r = sqrt((x-xdim/2)^2+(y-ydim/2)^2)
            if _r < 5
                print(c1)
            elseif _r < 20
                print(c2)
            else
                print(c3)
            end
        end
        println()
    end
end

bullseye()


text:
# ---
# title: 673. Number of Longest Increasing Subsequence
# id: problem673
# author: Indigo
# date: 2021-06-27
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/number-of-longest-increasing-subsequence/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, return _the number of longest increasing
# subsequences._
#
# **Notice** that the sequence has to be **strictly** increasing.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,3,5,4,7]
#     Output: 2
#     Explanation: The two longest increasing subsequences are [1, 3, 4, 7] and [1, 3, 5, 7].
#
#
# **Example 2:**
#
#
#
#     Input: nums = [2,2,2,2,2]
#     Output: 5
#     Explanation: The length of longest continuous increasing subsequence is 1, and there are 5 subsequences' length is 1, so output 5.
#
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 2000`
#   * `-106 <= nums[i] <= 106`
#
#
## @lc code=start
using LeetCode

function find_number_of_LIS(nums::Vector{Int})
    len = length(nums)
    len <= 1 && return len
    dp = fill(1, len)
    cnt = copy(dp)
    for i in 2:len, j in 1:(i - 1)
        nums[i] <= nums[j] && continue
        if dp[i] <= dp[j]
            dp[i] = dp[j] + 1
            cnt[i] = cnt[j]
        elseif dp[i] == dp[j] + 1
            cnt[i] += cnt[j]
        end
    end
    max_len = maximum(dp)
    return sum(c for (idx, c) in enumerate(cnt) if dp[idx] == max_len)
end
## @lc code=end


text:
export @init

macro definit()
  quote
    if !isdefined(:__inits__)
      const $(esc(:__inits__)) = Function[]
    end
    if !isdefined(:__init__)
      $(esc(:__init__))() = @init
    end
  end
end

function initm(ex)
  quote
    @definit
    push!($(esc(:__inits__)), () -> $(esc(ex)))
    nothing
  end
end

function initm()
  :(for f in __inits__
      f()
    end) |> esc
end

macro init(args...)
  initm(args...)
end

"Prevent init fns being called multiple times during precompilation."
macro guard(ex)
  :(!isprecompiling() && $(esc(ex)))
end


text:
using Base.Test
using QuantumOptics

@testset "basis" begin

shape1 = [5]
shape2 = [2, 3]
shape3 = [6]

b1 = GenericBasis(shape1)
b2 = GenericBasis(shape2)
b3 = GenericBasis(shape3)

@test b1.shape == shape1
@test b2.shape == shape2
@test b1 != b2
@test b1 != FockBasis(2)
@test b1 == b1

@test tensor(b1) == b1
comp_b1 = tensor(b1, b2)
comp_uni = b1 ⊗ b2
comp_b2 = tensor(b1, b1, b2)
@test comp_b1.shape == [prod(shape1), prod(shape2)]
@test comp_uni.shape == [prod(shape1), prod(shape2)]
@test comp_b2.shape == [prod(shape1), prod(shape1), prod(shape2)]

@test b1^3 == CompositeBasis(b1, b1, b1)
@test (b1⊗b2)^2 == CompositeBasis(b1, b2, b1, b2)
@test_throws ArgumentError b1^(0)

comp_b1_b2 = tensor(comp_b1, comp_b2)
@test comp_b1_b2.shape == [prod(shape1), prod(shape2), prod(shape1), prod(shape1), prod(shape2)]
@test comp_b1_b2 == CompositeBasis(b1, b2, b1, b1, b2)

@test_throws ArgumentError tensor()
@test comp_b2.shape == tensor(b1, comp_b1).shape
@test comp_b2 == tensor(b1, comp_b1)

@test_throws ArgumentError ptrace(comp_b1, [1, 2])
@test ptrace(comp_b2, [1]) == ptrace(comp_b2, [2]) == comp_b1 == ptrace(comp_b2, 1)
@test ptrace(comp_b2, [1, 2]) == ptrace(comp_b1, [1])
@test ptrace(comp_b2, [2, 3]) == ptrace(comp_b1, [2])

comp1 = tensor(b1, b2, b3)
comp2 = tensor(b2, b1, b3)
@test permutesystems(comp1, [2,1,3]) == comp2

@test !QuantumOptics.bases.equal_bases([b1, b2], [b1, b3])
@test !QuantumOptics.bases.multiplicable(comp1, b1 ⊗ b2 ⊗ NLevelBasis(prod(b3.shape)))

end # testset


text:
# via https://github.com/rpmuller/pyquante2/blob/4b05475c3616fb48b2536f07023ddc6820dbfccc/othertests/int2e-speed-julia.ipynb

using BenchmarkTools # more accurate @btime macro

println("vrrspeedtest.jl: Compare times between the test cases run using Python/C code (int2e-speed.ipynb) and Julia code.")

function pythonsig_vrr(ax,ay,az,na,aI,aJ,aK,aexpn,bx,by,bz,nb,bexpn,
    cx,cy,cz,nc,cI,cJ,cK,cexpn,dx,dy,dz,nd,dexpn,M)
    # Helper function for vrr, since I for some reason changed the call signature
    # from the python version. This version has the same signature as the python version
    # and calls the Julia version
    na*nb*nc*nd*Quante.vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
        cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,M)
end

function pythonsig_vrr_iter(ax,ay,az,na,aI,aJ,aK,aexpn,bx,by,bz,nb,bexpn,
    cx,cy,cz,nc,cI,cJ,cK,cexpn,dx,dy,dz,nd,dexpn,M)
    # Helper function for vrr, since I for some reason changed the call signature
    # from the python version. This version has the same signature as the python version
    # and calls the Julia version
    na*nb*nc*nd*Quante.vrr_iter(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
        cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,M)
end

# Timing comparisons
# Python/C times are from the int2e-speed.ipynb notebook in the same folder.

println("s orbitals on the same center:")
println("RMuller: 9 usec time Julia compared to 0.75 usec in C and 49 usec in Python")
println("Jarv/Asriel commit c2190d5: 1.292 μs (47 allocations: 2.44 KiB) //  15.969 μs (77 allocations: 3.38 KiB)")
#                       xyza   na   lmna   aa   xyzb    nb ab   xyzc      nc  lmnc   ac  xyzd   nd ad  M
@btime pythonsig_vrr(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0)
@btime pythonsig_vrr_iter(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0)

println("sp orbitals on the same center:")
println("RMuller: 0.95 usec C, 58 usec Python, 13 usec Julia")
println("Jarv/Asriel commit c2190d5:   3.413 μs (133 allocations: 7.17 KiB) //  16.397 μs (91 allocations: 3.63 KiB)")
#                       xyza   na  lmna   aa  xyzb   nb ab  xyzc   nc  lmnc   ac  xyzd   nd ad  M
@btime pythonsig_vrr(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 1,0,0, 1., 0.,0.,0., 1.,1., 0)
@btime pythonsig_vrr_iter(0.,0.,0., 1., 0,0,0, 1., 0.,0.,0., 1.,1., 0.,0.,0., 1., 1,0,0, 1., 0.,0.,0., 1.,1., 0)

println("sp integrals, 4 different centers")
println("RMuller: 1.2 usec C, 66 usec Python, 23 Julia")
println("Jarv/Asriel commit c2190d5:   3.800 μs (131 allocations: 7.13 KiB) //  16.791 μs (89 allocations: 3.58 KiB)")
@btime pythonsig_vrr(0.,0.,0., 1., 0,0,0, 1., 0.,0.,1., 1.,1., 0.,1.,0., 1., 1,0,0, 1., 1.,1.,0., 1.,1., 0)
@btime pythonsig_vrr_iter(0.,0.,0., 1., 0,0,0, 1., 0.,0.,1., 1.,1., 0.,1.,0., 1., 1,0,0, 1., 1.,1.,0., 1.,1., 0)

println("pd integrals, 4 different centers")
println("RMuller: 3 usec C, 95 usec Python, 32 usec @time Julia")
println("Jarv/Asriel commit c2190d5:   20.108 μs (743 allocations: 40.44 KiB) //  18.694 μs (117 allocations: 4.31 KiB)")
@btime pythonsig_vrr(0.,0.,0., 1., 1,0,0, 1., 0.,0.,1., 1.,1., 0.,1.,0., 1., 1,1,0, 1., 1.,1.,0., 1.,1., 0)
@btime pythonsig_vrr_iter(0.,0.,0., 1., 1,0,0, 1., 0.,0.,1., 1.,1., 0.,1.,0., 1., 1,1,0, 1., 1.,1.,0., 1.,1., 0)


text:
export Hypercube
export find_generators
export generate_coordinates
export volume
export isequiv

import LinearAlgebra
import LinearAlgebraX

"""
    Hypercube(shape)

Represent a hypercubic (Bravais) lattice.

# Fields
* `shape_matrix`: a matrix whose columns are the lattice vectors of the
* `inverse_shape_matrix`: `(shape_matrix)`
* `wrap`: periodic wrapping function which takes an integer array, and maps it onto a
  site in the Bravais lattice. Analogous to `fldmod`.
"""
struct Hypercube
    shape_matrix::Matrix{Int}
    inverse_shape_matrix::Matrix{Rational{Int}}
    wrap::Function

    function Hypercube(shape_matrix::AbstractMatrix{<:Integer})
        LinearAlgebra.checksquare(shape_matrix)
        det = LinearAlgebraX.detx(shape_matrix) # does not check whether shape_matrix is square
        iszero(det) && throw(ArgumentError("scale matrix null"))
        inverse_shape_matrix = LinearAlgebraX.invx(shape_matrix)
        function wrap(r::AbstractArray{<:Integer}, mode::RoundingMode=RoundDown)
            rnd = (x) -> round(Int, x, mode)
            R = rnd.(inverse_shape_matrix * r)
            r2 = r - shape_matrix * R
            return R, r2
        end
        return new(shape_matrix, inverse_shape_matrix, wrap)
    end
end

"""
    dimension(cube::Hypercube)

Return the spatial dimension of the hypercube.
"""
dimension(cube::Hypercube) = size(cube.shape_matrix, 1)

"""
    volume(cube::Hypercube)

Return the signed volume of the hypercube, defined by the determinant of the shape.
"""
volume(cube::Hypercube) = LinearAlgebraX.detx(cube.shape_matrix)

function Base.:(==)(lhs::Hypercube, rhs::Hypercube)
    return (lhs.shape_matrix == rhs.shape_matrix)
end

"""
    isequiv(lhs::Hypercube, rhs::Hypercube)

Check whether the two hypercubes are equivalent.
"""
function isequiv(lhs::Hypercube, rhs::Hypercube)
    R, r = lhs.wrap(rhs.shape_matrix)
    isone(abs(LinearAlgebraX.detx(R))) || return false
    iszero(r) || return false
    return true
end

"""
    find_generators(cube::Hypercube)

Find generators of an Hypercube. For one-dimensional lattice, which is isomorphic to Zₙ, the
generator is +1. For two-dimensional lattices, this function invokes
[`find_generators_2d`](@ref). Higher dimensions are currently not supported.
"""
function find_generators(cube::Hypercube)
    if dimension(cube) == 1
        return ones(Int, 1, 1)
    elseif dimension(cube) == 2
        return find_generators_2d(cube)
    else
        error("Currently, only 1D and 2D lattices are supported")
    end
end

"""
    find_generators_2d(cube::Hypercube)

Find translation generators of an Hypercube.
"""
function find_generators_2d(cube::Hypercube)
    dimension(cube) != 2 && throw(ArgumentError("shape matrix should be 2x2"))

    cubevolume = abs(volume(cube))
    allowed_pairs = [[1,0]]
    append!(
        allowed_pairs,
        [x, y] for y in 1:cubevolume, x in 0:cubevolume if gcd(x,y) == 1
    )

    function make_loop(t::Vector{Int})
        axis = Vector{Int}[zero(t)]
        r = cube.wrap(t)[2]
        while !iszero(r)
            push!(axis, r)
            r = cube.wrap(r .+ t)[2]
        end
        @assert allunique(axis)
        return axis
    end

    for r1p in allowed_pairs
        r1p[1] == 0 && continue
        r1 = [r1p[1], -r1p[2]]
        axis1 = make_loop(r1)
        iszero(mod(cubevolume, length(axis1))) || continue
        for r2 in allowed_pairs
            isone(LinearAlgebraX.detx(hcat(r1, r2))) || continue
            axis2 = make_loop(r2)
            cubevolume == length(axis1) * length(axis2) || continue
            coordinates = vec([cube.wrap(sum(x))[2] for x in Iterators.product(axis1, axis2)])
            allunique(coordinates) && return hcat(r1, r2)
        end
    end
    error("Failed to find decomposition") # COV_EXCL_LINE
end

"""
    generate_coordinates(cube, generators)

Generate a list of coordinates of the hypercube

# Arguments
* `cube::Hypercube`
* `generator_translations::AbstractMatrix{<:Integer}`
"""
function generate_coordinates(
    cube::Hypercube,
    generator_translations::AbstractMatrix{<:Integer},
)
    let dim = dimension(cube)
        if size(generator_translations) != (dim, dim)
            throw(DimensionMismatch("Hypercube and generator_translations have different dimensions"))
        elseif !isone(abs(LinearAlgebraX.detx(generator_translations)))
            throw(ArgumentError("generator_translations is not unimodular"))
        end
    end
    function make_loop(t::AbstractVector{<:Integer})
        axis = Vector{Int}[zero(t)]
        r = cube.wrap(t)[2]
        while !iszero(r)
            push!(axis, r)
            r = cube.wrap(r .+ t)[2]
        end
        return axis
    end
    axes = [make_loop(t) for t in eachcol(generator_translations)]
    coordinates = vec([cube.wrap(sum(x))[2] for x in Iterators.product(axes...)])

    if length(coordinates) != abs(volume(cube))
        msg = "translations $generator_translations generates $(length(coordinates)) " *
            "coordinates, while the volume is $(abs(volume(cube)))"
        throw(ArgumentError(msg))
    end
    @assert allunique(coordinates)
    return coordinates
end


text:
using PushshiftRedditDistiller
using Test
using DataDeps

ENV["DATADEPS_ALWAYS_ACCEPT"] = true
ENV["CI"] = true

@testset "Empty Filters" begin
    test_filter = RedditDataFilter()
    @test isempty(test_filter.author)
    @test isempty(test_filter.subreddit)
    @test isempty(test_filter.fields)

    sample_submissions = distill(datadep"reddit-submissions-sample", test_filter)
    @test length(sample_submissions) == 1000

    sample_comments = distill(datadep"reddit-comments-sample", test_filter)
    @test length(sample_comments) == 10000
end

@testset "Filters" begin
    author = "[deleted]"
    filter_author = RedditDataFilter(author=[author])
    deleted_author_comments = distill(datadep"reddit-comments-sample", filter_author)
    @test length(deleted_author_comments) > 0
    @test all(get.(deleted_author_comments, :author, nothing) .== author)

    subreddit = "AskReddit"
    filter_subreddit = RedditDataFilter(subreddit=[subreddit])
    test_subreddit_comments = distill(datadep"reddit-comments-sample", filter_subreddit)
    @test length(test_subreddit_comments) > 0
    @test all(get.(test_subreddit_comments, :subreddit, nothing) .== subreddit)

    field = "created_utc"
    filter_field = RedditDataFilter(fields=[field])
    test_field_comments = distill(datadep"reddit-comments-sample", filter_field)
    @test all(first.(keys.(test_field_comments)) .== Symbol(field))
end


text:
using GeoInterface
using LibGEOS

using LinearAlgebra
using LaTeXStrings
using Formatting
import PyPlot; plt=PyPlot
cbbPalette = ["#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7"]
# using ..NYC: read_distr_reprpoints
using Distributions: MultivariateNormal
import GeometricalPredicates
using VoronoiDelaunay
import VoronoiDelaunay: getx, gety, DelaunayEdge, DelaunayTriangle

using PyCall
@pyimport mpl_toolkits.mplot3d as mplot3d

function plot_polygon(poly_coords, facecolor, alpha;
                      edgecolor="none", linestyle="None", linewidth=0.0, kwargs...)
    ax = plt.gca()
    poly_array = Array{Float64}(undef, length(poly_coords), 2)
    for i in 1:length(poly_coords)
        poly_array[i, :] = poly_coords[i]
    end
    polygon = plt.plt[:Polygon](poly_array, true; facecolor=facecolor, kwargs...)
    p = ax[:add_patch](polygon)
    p[:set_edgecolor](edgecolor)
    p[:set_linestyle](linestyle)
    p[:set_linewidth](linewidth)
    p[:set_alpha](alpha)
end

function plot_buffer(shape, border, color, dist, alpha; kwargs...)
    border_buffer = LibGEOS.buffer(border, dist)
    shape_buffer = LibGEOS.intersection(border_buffer, shape)
    if typeof(shape_buffer) <: LibGEOS.Polygon
        poly_coords = GeoInterface.coordinates(shape_buffer)[1]
        plot_polygon(poly_coords, color, alpha; kwargs...)
    elseif typeof(shape_buffer) <: LibGEOS.MultiPolygon
        mpoly_coords = GeoInterface.coordinates(shape_buffer)
        for poly_coords in mpoly_coords
            plot_polygon(poly_coords[1], color, alpha; kwargs...)
        end
    else
        println(typeof(shape_buffer))
    end
end

function plot_schdistr_labels(;fontsize=12.0, labelcolor="black")
    plt.rc("text", usetex=false)
    # where to place the school district labels:
    schdistr_represent = Dict{SchDistr, Tuple{Float64, Float64}}(25 => (1035426.576356058, 215708.47750849032), 20 => (980957.6301213544, 166113.83169550286), 3 => (992826.0954164147, 226695.6766967253), 4 => (1000248.2836855762, 228298.3092956051), 16 => (1004214.0439845352, 188657.53201288398), 23 => (1008819.7413809987, 183460.7804870081), 18 => (1008787.3316119347, 173443.51080316756), 26 => (1054423.2425350049, 212717.9963988735), 32 => (1006895.1698167967, 192470.89239496225), 1 => (988952.4509786353, 202318.4309997029), 2 => (987095.1379395023, 210542.70410151145), 24 => (1013178.9222047399, 202219.43350214284), 28 => (1036851.1969587662, 195628.9613036562), 10 => (1008845.8110727321, 258295.0202941391), 9 => (1008189.8006258132, 245065.8830871096), 7 => (1006949.3221918349, 235790.44049067114), 30 => (1011393.8551870013, 216959.78259271727), 12 => (1017667.4294655474, 243924.23712153034), 8 => (1022714.7969213917, 237308.7890929639), 14 => (998121.9324939251, 199650.8938903253), 31 => (947923.9119948001, 147916.42111200799), 17 => (997069.2552085526, 180290.71820062993), 13 => (994552.5771594604, 190230.17111200746), 15 => (987418.7441204586, 181218.70071406086), 22 => (1001740.2660473757, 162368.14001459582), 21 => (991081.0983826279, 157637.68869012795), 11 => (1026344.6177574333, 256966.2814940873), 29 => (1055388.4397568335, 188398.91879267507), 19 => (1017942.166879338, 181833.5946349586), 27 => (1040711.4984081909, 178979.62564887386), 5 => (996548.0388532034, 237031.53149408224), 6 => (1001541.6608069859, 248002.3049010716))
    xmin,xmax = plt.xlim()
    ymin,ymax = plt.ylim()
    ax = plt.gca()
    for schdistr in keys(schdistr_represent)
        center_x, center_y = schdistr_represent[schdistr]

        if !(xmin < center_x < xmax)
            continue
        elseif !(ymin < center_y < ymax)
            continue
        end

        ax[:annotate](
            @sprintf("%d", schdistr),
            xy=(center_x, center_y),
            xycoords="data",
            xytext=(0,0),
            textcoords="offset points",
            horizontalalignment="center",
            verticalalignment="center",
            fontweight="black",
            color=labelcolor,
            alpha=1.0,
            fontsize = fontsize,
            zorder=3,
            annotation_clip=true,
        )
        ax[:annotate](
            @sprintf("%d", schdistr),
            xy=(center_x, center_y),
            xycoords="data",
            xytext=(0,0),
            textcoords="offset points",
            horizontalalignment="center",
            verticalalignment="center",
            fontweight="black",
            color=labelcolor,
            alpha=0.5,
            fontsize = fontsize,
            zorder=20,
            annotation_clip=true,
        )
    end
    plt.rc("text", usetex=true)
end

function plot_all_pairs(τpost_pairs::Dict{Tuple{SchDistr,SchDistr}}, regionData_dict::Dict{SchDistr,GeoRDD.RegionData}; scaleup=200.0, zorder="none")
    for distr_pair in keys(τpost_pairs)
        distrA, distrB = distr_pair

        τpost = τpost_pairs[distr_pair]
        # eff_size = abs(mean(τpost)) / std(τpost) * 2
        # if eff_size < 2.0
            # continue
        # end

        if distrA ∈ (1,2) || distrB ∈ (1,2)
            # don't do Manhattan, it's weird
            continue
        end

        # skip districts that aren't on the map
        if distrA ∉ (13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 32)
            continue
        end
        if distrB ∉ (13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 32)
            continue
        end

        if distrA >= distrB
            continue
        end

        shapeA = regionData_dict[distrA].shape
        shapeB = regionData_dict[distrB].shape
        border = GeoRDD.get_border(shapeA, shapeB, 10.0)

        border_coords = GeoInterface.coordinates(border)
        if typeof(border) <: LibGEOS.MultiLineString
        else
            border_coords = [border_coords]
        end
        ax = plt.gca()

        if mean(τpost) < 0
            cheaper_shape = shapeA
            pricier_shape = shapeB
        else
            cheaper_shape = shapeB
            pricier_shape = shapeA
        end
        plot_buffer(pricier_shape, border, cbbPalette[2], (abs(mean(τpost)) + std(τpost))*scaleup, 1.0, zorder=zorder-1)
        # if abs(mean(τpost)) > std(τpost)
            # plot_buffer(pricier_shape, border, "none", (abs(mean(τpost)) - std(τpost))*scaleup, 1.0, zorder=zorder+1, edgecolor="white", linestyle=":", linewidth=1.0)
        # end
        plot_buffer(pricier_shape, border, cbbPalette[1], abs(mean(τpost))*scaleup, 1.0, zorder=zorder)
    end
end

function plot_streets(;zoom=11, brighten=30, rgb=(200, 200, 200), zorder=-9)
    ax = plt.gca()
    ax[:set_autoscale_on](false)
    ax[:set_autoscalex_on](false)
    ax[:set_autoscaley_on](false)
    EPSG=2263 # projection
    topleft =  (40.75, -74.01)
    botright = (40.58, -73.60)
    py"""
    import matplotlib.pyplot as plt

    import geopandas as gpd
    import tilemapbase
    import pyproj
    import PIL.Image as _Image

    def bounds(plotter, epsg):
        epsg_proj = pyproj.Proj(init='epsg:%d' % epsg, preserve_units=True)
        epsg_3857 = pyproj.Proj(init='epsg:%d' % tilemapbase.mapping._WEB_MERCATOR)
        self = plotter
        scale = 2 ** self.zoom
        x0, y0 = self.xtilemin / scale, self.ytilemin / scale # in tile space
        x1, y1 = (self.xtilemax + 1) / scale, (self.ytilemax + 1) / scale
        x0_ll, y0_ll = tilemapbase.to_lonlat(x0, y0)
        x1_ll, y1_ll = tilemapbase.to_lonlat(x1, y1)
        y0_ll, y1_ll = y1_ll, y0_ll # for some reason that I don't understand
        print(y0_ll)
        print(y1_ll)
        x0_epsg, y0_epsg = epsg_proj(x0_ll, y0_ll)
        x1_epsg, y1_epsg = epsg_proj(x1_ll, y1_ll)
        return x0_epsg, x1_epsg, y0_epsg, y1_epsg

    def brighten(img, factor):
        # split the image into individual bands
        source = img.split()
        R, G, B, A = 0, 1, 2, 3

        # process each band separately
        out_R = source[R].point(lambda x: min(x+factor, 256))
        out_G = source[G].point(lambda x: min(x+factor, 256))
        out_B = source[B].point(lambda x: min(x+factor, 256))

        # build a new multiband image
        brightened = _Image.merge(img.mode, (out_R, out_G, out_B, source[A]))
        return brightened

    def makecolor(img, rgb):
        # split the image into individual bands
        source = img.split()
        R, G, B, A = 0, 1, 2, 3

        # process each band separately
        out_R = source[R].point(lambda x: rgb[0])
        out_G = source[G].point(lambda x: rgb[1])
        out_B = source[B].point(lambda x: rgb[2])

        # build a new multiband image
        new_im = _Image.merge(img.mode, (out_R, out_G, out_B, source[A]))
        return new_im

    tilemapbase.start_logging()
    tilemapbase.init("/Users/imolk/tmp/tilemapbase.cache", create=True)

    extent = tilemapbase.Extent.from_lonlat($topleft[1], $botright[1],
                      $botright[0], $topleft[0])
    # plotter_ln = tilemapbase.Plotter(extent,tilemapbase.tiles.Carto_Dark_No_Labels, zoom=$zoom)
    plotter_ln = tilemapbase.Plotter(extent,tilemapbase.tiles.Stamen_Terrain_Lines, zoom=$zoom)
    plotter = plotter_ln
    # tile = brighten(plotter.as_one_image(), $brighten)
    tile = makecolor(plotter.as_one_image(), $rgb)
    $ax.imshow(tile, interpolation="lanczos", extent=bounds(plotter, $EPSG), origin="upper", zorder=$zorder)
    """
end

function ticks(xmin, xmax, tickdist)
    return tickdist * (ceil(xmin/tickdist):1.0:floor(xmax/tickdist))
end
function ticks_km(tickdist::Float64)
    _xlim = plt.xlim()
    _ylim = plt.ylim()
    km_in_ft = 3280.84
    _xlim_km = _xlim ./ km_in_ft
    _ylim_km = _ylim ./ km_in_ft
    _xtick_km = ticks(_xlim_km[1], _xlim_km[2], tickdist)
    _ytick_km = ticks(_ylim_km[1], _ylim_km[2], tickdist)
    _xtick_ft = _xtick_km .* km_in_ft
    _ytick_ft = _ytick_km .* km_in_ft
    plt.xticks(_xtick_ft, [@sprintf("%.0f", x) for x in _xtick_km])
    plt.yticks(_ytick_ft, [@sprintf("%.0f", x) for x in _ytick_km])
    plt.xlabel("Eastings (km)")
    plt.ylabel("Northings (km)")
end

function plot_sales(filtered, missing_sqft; edgecolor="black", labelsize=12.0, labelcolor="black", colorbar=true, tickdist=5.0)
    fig=plt.figure()
    background_schdistrs(plt.gca(),
                color="#AAAAAA",
                edgecolor=edgecolor,
                alpha=1.0,
                linestyle="-",
                zorder=-10
                )
    background_schdistrs(plt.gca(),
                color="none",
                edgecolor=edgecolor,
                alpha=1.0,
                linestyle="-",
                linewidth=0.5,
                zorder=11)
    _xlim = (0.98e6, 1.06e6)
    _ylim = (1.5e5, 2.2e5)
    plt.plot(
        missing_sqft[:XCoord],
        missing_sqft[:YCoord],
        color="white",
        marker="x",
        markersize=1.0,
        linestyle="",
        zorder=8,
        alpha=0.7,
        linewidth=0.1,
        markeredgewidth=0.3,
    )
    plt.scatter(filtered[:XCoord],
        filtered[:YCoord],
        c=filtered[:logSalePricePerSQFT],
        cmap="jet",
        edgecolors="None",
        zorder=10,
        marker="o",
    #     vmin=4.0,
    #     vmax=8.0,
        alpha=1.0,
        s=1,
        )
    # plt.xlim(np.percentile(sales_geocoded_ggl["XCoord"].dropna(), [1,99]))
    # plt.ylim(np.percentile(sales_geocoded_ggl["YCoord"].dropna(), [1,99]))
    # plt.axis("off")
    # plt.title("log(Property Price per SQFT) in New York")
    ax = plt.gca()
    ax[:set](aspect="equal", adjustable="box")
    # ax[:get_xaxis]()[:set_ticks]([])
    # ax[:get_yaxis]()[:set_ticks]([])
    plt.xlim(_xlim)
    plt.ylim(_ylim)
    ticks_km(tickdist)

    ax[:set_facecolor]("#B0DAEE")
    if colorbar
        cbar = plt.colorbar(label="Square foot prices (\\\$ per sqft)")
        ax = cbar[:ax]
        tick_template = [1, 2, 3, 4, 5, 6, 7, 8, 9]
        labeled = [10, 100, 1000, 10^4]
        ax[:set_yticks](log.(labeled))
        ax[:set_yticks](log.([tick_template; tick_template.*10; tick_template.*100; tick_template.*1000]),
                               minor=true)
        ax[:set_yticklabels]([format(x, commas=true) for x in labeled])
    end

    ax[:set_autoscale_on](false)
    ax[:set_autoscalex_on](false)
    ax[:set_autoscaley_on](false)
    plot_schdistr_labels(fontsize=labelsize, labelcolor=labelcolor)
end
function plot_cliffface(μ, Σ, color; ndraws::Int=10, label=L"posterior of $\tau(x)$ on residuals")
    plt.plot(μ, color=color, ".", label=label)
    msize=10
    plt.plot(1,μ[1], color=cbbPalette[5],
        markersize=msize, marker="o")
    plt.plot(length(μ),μ[end], color=cbbPalette[4],
        markersize=msize, marker="o")
    posterior_sd = sqrt.(diag(Σ))
    plt.fill_between(1:length(μ), μ.-2*posterior_sd, μ.+2*posterior_sd,
        color=color,
        alpha=0.3,
        linewidth=0,
        zorder=-2,
        )
    posterior_distr = MultivariateNormal(μ, Σ)
    for _ in 1:ndraws
        plt.plot(1:length(μ), rand(posterior_distr), "-", color="white", alpha=0.4, linewidth=1.5, zorder=-1)
    end
    plt.xlim(1,length(μ))
    plt.ylim(minimum(μ.-2*posterior_sd), maximum(μ.+2*posterior_sd))

    # second y-axis for percentage price increase
    ax = plt.gca()
    y1, y2=ax[:get_ylim]()
    x1, x2=ax[:get_xlim]()
    ax2=ax[:twinx]()
    ax2[:set_ylim](y1, y2)
    yticks_frac = collect(0.5: 0.1 : 2.0)
    yticks_log = log.(yticks_frac)
    in_window = y1 .< yticks_log .< y2
    ax2[:set_yticks](yticks_log[in_window])
    ylabels = [@sprintf("%.0f\\%%", y*100) for y in yticks_frac[in_window]]
    ax2[:set_yticklabels](ylabels)

    plt.sca(ax)
end

##########################################
####### 3D plot ##########################
##########################################

struct IndexablePoint <: AbstractPoint2D
    _x::Float64
    _y::Float64
    index::Int64
end
getx(p::IndexablePoint) = p._x
gety(p::IndexablePoint) = p._y
get_index(p::IndexablePoint) = p.index
IndexablePoint(x::Float64, y::Float64) = IndexablePoint(x, y, -1)

function trisurf_custom_delaunay(grid::Matrix, domain)
    # create Delaynay
    tess = DelaunayTessellation2D{IndexablePoint}()
    minx, maxx = extrema(grid[1,:])
    xwidth = maxx-minx
    miny, maxy = extrema(grid[2,:])
    ywidth = maxy-miny
    delaunay_width = VoronoiDelaunay.max_coord - VoronoiDelaunay.min_coord
    min_coord = VoronoiDelaunay.min_coord
    transform_x = x -> (x-minx) / xwidth * (delaunay_width) + min_coord
    inverse_x = x -> (x - min_coord) / delaunay_width * xwidth + minx
    transform_y = y -> (y-miny) / ywidth * (delaunay_width) + min_coord
    inverse_y = y -> (y - min_coord) / delaunay_width * ywidth + miny

    points = [IndexablePoint(transform_x(grid[1,i]), transform_y(grid[2,i]), i) for i in 1:size(grid,2)]
    push!(tess, points)

    # get minimum triangle areas
    areas = Float64[]
    centres = LibGEOS.Point[]
    for tri in tess
        area = GeometricalPredicates.area(tri)
        centre = GeometricalPredicates.centroid(tri)
        centre_LibGEOS = LibGEOS.Point(inverse_x(getx(centre)), inverse_y(gety(centre)))
        push!(centres, centre_LibGEOS)
        push!(areas, area)
    end
    medarea = median(areas)

    # build up triangle indices
    triangle_indices = Tuple{Int,Int,Int}[]
    for tri in tess
        area = GeometricalPredicates.area(tri)
        if area > 2*medarea
            # skip triangles that are too big
            continue
        end

        centre = GeometricalPredicates.centroid(tri)
        centre_LibGEOS = LibGEOS.Point(inverse_x(getx(centre)), inverse_y(gety(centre)))
        if !LibGEOS.within(centre_LibGEOS, domain)
            continue
        end
#         a_x, a_y = inverse_x(getx(geta(tri))), inverse_y(gety(geta(tri)))
#         b_x, b_y = inverse_x(getx(getb(tri))), inverse_y(gety(getb(tri)))
#         c_x, c_y = inverse_x(getx(getc(tri))), inverse_y(gety(getc(tri)))
#         tri_geos = LibGEOS.Polygon([[[a_x, a_y], [b_x, b_y], [c_x, c_y], [a_x, a_y]]])
#         if LibGEOS.area(LibGEOS.intersection(tri_geos, domain)) < (0.9*area)
#             print("skip area!")
#             continue
#         end
        a, b, c = geta(tri), getb(tri), getc(tri)
        ia, ib, ic = get_index(a), get_index(b), get_index(c)
        push!(triangle_indices, (ia-1, ib-1, ic-1)) # python is zero-indexed
    end
    return triangle_indices
end

function plot_surface3d(gridT, gridC, gpT, gpC, Xb)
    predgridT=GaussianProcesses.predict(gpT, gridT; full_cov=false)
    predgridC=GaussianProcesses.predict(gpC, gridC; full_cov=false)
    predT_b = GaussianProcesses.predict(gpT, Xb; full_cov=false)
    predC_b = GaussianProcesses.predict(gpC, Xb; full_cov=false)
    plot_surface3d(gridT, gridC, predgridT, predgridC, predT_b, predC_b, Xb)
end
function plot_cliff(Xb::Matrix, predT_b::Vector, predC_b::Vector)
    # verts_T = collect(zip(Xb[1,:], Xb[2,:], predT_b.-1e-4))
    # verts_C = collect(zip(Xb[1,:], Xb[2,:], predC_b.+1e-4))
    # verts = [verts_T; verts_C[end:-1:1]]
    n_b = length(predT_b)
    cliff = plt.plot_surface([Xb[1,:] Xb[1,:]],
                     [Xb[2,:] Xb[2,:]],
                     [(predT_b.-0.001) (predC_b.+0.001)],
#                      [(ones(n_b)*5.0) (ones(n_b)*5.7)],
                     shade=true,
                     color="#BBBBBB",
                     rstride=1,
                     cstride=1,
                     linewidth=0,
                     antialiased=false,
                     )

#     polygon = mplot3d.art3d[:Poly3DCollection]([verts])
#     polygon[:set_color]("#AAAAAA")
#     polygon[:set_hatch]("///////")
#     polygon[:set_edgecolor]("#333333")
#     polygon[:set_linewidth](0)
#     plt.gca()[:add_collection3d](polygon)
#     return polygon
    return cliff
end
function plot_surface3d(gridT, gridC, predgridT, predgridC, predT_b, predC_b, Xb, regionT, regionC;
                        xlim=(-Inf, Inf), ylim=(-Inf, Inf), labelT="", labelC="")
    trisurf_X = vcat(@view(gridT[1,:]), @view(gridC[1,:]))
    trisurf_Y = vcat(@view(gridT[2,:]), @view(gridC[2,:]))
    trisurf_Z = vcat(predgridT, predgridC)

    bX_inside = xlim[1] .<= Xb[1,:] .<= xlim[2]
    bY_inside = ylim[1] .<= Xb[2,:] .<= ylim[2]
    b_inside = bX_inside .& bY_inside
    Xb_inside = Xb[:, b_inside]

    min_b = min(minimum(predT_b[b_inside]),minimum(predC_b[b_inside])) - 0.1
    max_b = max(maximum(predT_b[b_inside]),maximum(predC_b[b_inside])) + 0.1

    TX_inside = xlim[1] .<= gridT[1,:] .<= xlim[2]
    TY_inside = ylim[1] .<= gridT[2,:] .<= ylim[2]
    T_inside = TX_inside .& TY_inside

    gridT_augmented = hcat(@view(gridT[:,T_inside]), @view(Xb[:,b_inside]))
    trianglesT = trisurf_custom_delaunay(gridT_augmented, regionT)
    trisurf_T = plt.plot_trisurf(
       [gridT[1,T_inside]; Xb[1,b_inside]],
       [gridT[2,T_inside]; Xb[2,b_inside]],
       [predgridT[T_inside]; predT_b[b_inside]],
       triangles=trianglesT,
       linewidth=0,
       cmap=plt.cm_get_cmap("jet"),
       vmin=min_b,
       vmax=max_b,
       alpha=1.0
       )

    CX_inside = xlim[1] .<= gridC[1,:] .<= xlim[2]
    CY_inside = ylim[1] .<= gridC[2,:] .<= ylim[2]
    C_inside = CX_inside .& CY_inside
    gridC_augmented = hcat(@view(gridC[:,C_inside]), @view(Xb[:,b_inside]))
    trianglesC = trisurf_custom_delaunay(gridC_augmented, regionC)
    trisurf_C = plt.plot_trisurf(
       [gridC[1,C_inside]; Xb[1,b_inside]],
       [gridC[2,C_inside]; Xb[2,b_inside]],
       [predgridC[C_inside]; predC_b[b_inside]],
       triangles=trianglesC,
       linewidth=0,
       antialiased=false,
       edgecolor="none",
       cmap=plt.cm_get_cmap("jet"),
       vmin=min_b,
       vmax=max_b,
       alpha=1.0
       )

    distances = .√(diff(Xb_inside[1,:]).^2 + diff(Xb_inside[2,:]).^2)
    med_dist = median(distances)
    breakpoints = find(distances .> (2*med_dist))
    push!(breakpoints, length(distances)+1)
    cliffs = []
    seg_start=1
    for isegment in 1:length(breakpoints)
        seg_end = breakpoints[isegment]
        seg_indices = seg_start:seg_end
        if length(seg_indices) > 2
            cliff = plot_cliff(Xb_inside[:,seg_indices],
                               predT_b[b_inside][seg_indices],
                               predC_b[b_inside][seg_indices])
            push!(cliffs, cliff)
        end
        seg_start = seg_end + 1
    end

    centreC = mean(gridC[:,C_inside], 2)
    centreT = mean(gridT[:,T_inside], 2)
    plt.text3D(centreC[1], centreC[2], predC_b[50], labelC, horizontalalignment="left")
    plt.text3D(centreT[1], centreT[2], predT_b[50], labelT,
        color="white", fontweight=200, horizontalalignment="right")
    ax=plt.gca()
    plt.xlabel("Eastings")
    plt.ylabel("Northings")
    plt.zlabel("log Price per SQFT")
    ax[:view_init](elev=60, azim=230)
    ax[:set_xticklabels]([])
    ax[:set_yticklabels]([])
    _zlim = plt.zlim()
    plt.zticks(ceil(_zlim[1], 1):0.2:floor(_zlim[2]; digits=1))
#     plt.xlim(xlim)
#     plt.ylim(ylim)

#     trisurf_T[:set_zsort](3)
#     cliff[:set_zsort](2)
#     trisurf_C[:set_zsort](1)
#     lineT[:set_sort_zpos](3)
    trisurf_T[:set_sort_zpos](4)
    for cliff in cliffs
        cliff[:set_sort_zpos](2)
    end
    trisurf_C[:set_sort_zpos](1)
#     lineC[:set_sort_zpos](4)
end


text:
abstract type AbstractSatellite{S <: AbstractGNSS, T} <: AbstractEmitter{T} end

struct ConstantDopplerSatellite{
    S <: AbstractGNSS,
    T <: AbstractFloat,
    D <: Union{SVector{3}, AbstractDOA},
    E <: Union{Bool, AbstractExistence},
    C <: Union{Unitful.Level, DynamicCN0}
} <: AbstractSatellite{S, T}
    system::S
    prn::Int
    carrier_doppler::typeof(1.0Hz)
    carrier_phase::Float64
    code_phase::Float64
    cn0::C
    exists::E
    doa::D
    code::Vector{Int8}
    signal::StructArray{Complex{T},1,NamedTuple{(:re, :im),Tuple{Array{T,1},Array{T,1}}},Int}
end

function ConstantDopplerSatellite(
    ::Type{T},
    system::S,
    prn::Int;
    carrier_doppler = NaN*Hz,
    carrier_phase = NaN,
    code_phase = NaN,
    cn0::C = 45dBHz,
    exists::E = true,
    doa::D = SVector(0,0,1),
    velocity = 14_000.0m / 3.6s,
    distance_from_earth_center = 26_560_000.0m
) where {
    S <: AbstractGNSS,
    T <: AbstractFloat,
    D <: Union{SVector{3}, AbstractDOA},
    E <: Union{Bool, AbstractExistence},
    C <: Union{Unitful.Level, DynamicCN0},
}
    if isnan(carrier_doppler)
        carrier_doppler = calc_doppler(
            distance_from_earth_center,
            doa,
            velocity,
            get_center_frequency(system)
        )
    end
    if isnan(carrier_phase)
        sat_user_distance = calc_sat_user_distance(distance_from_earth_center, doa)
        carrier_phase = calc_carrier_phase(
            sat_user_distance,
            get_center_frequency(system) + carrier_doppler
        )
    else
        carrier_phase = carrier_phase / 2π
    end
    if isnan(code_phase)
        sat_user_distance = calc_sat_user_distance(distance_from_earth_center, doa)
        code_doppler = carrier_doppler * get_code_center_frequency_ratio(system)
        code_phase = calc_code_phase(
            sat_user_distance,
            get_code_frequency(system) + code_doppler,
            get_code_length(system) * get_secondary_code_length(system)
        )
    end
    code = Vector{Int8}(undef, 0)
    signal = StructArray{Complex{T}}(undef, 0)
    ConstantDopplerSatellite{S,T,D,E,C}(
        system,
        prn,
        carrier_doppler,
        carrier_phase,
        code_phase,
        cn0,
        exists,
        doa,
        code,
        signal
    )
end

function ConstantDopplerSatellite(
    system,
    prn;
    carrier_doppler = NaN*Hz,
    carrier_phase = NaN,
    code_phase = NaN,
    cn0 = 45dBHz,
    exists = true,
    doa = SVector(0,0,1),
    velocity = 14_000.0m / 3.6s,
    distance_from_earth_center = 26_560_000.0m
)
    ConstantDopplerSatellite(
        Float32,
        system,
        prn,
        carrier_doppler = carrier_doppler,
        carrier_phase = carrier_phase,
        code_phase = code_phase,
        cn0 = cn0,
        exists = exists,
        doa = doa,
        velocity = velocity,
        distance_from_earth_center = distance_from_earth_center
    )
end

function gen_signal!(
    sat::ConstantDopplerSatellite{S, T},
    sampling_frequency,
    intermediate_frequency,
    n0,
    num_samples::Integer,
    rng
) where {S <: AbstractGNSS, T}
    resize!(sat.code, num_samples)
    resize!(sat.signal, num_samples)
    carrier = gen_carrier!(
        sat.signal,
        intermediate_frequency + get_carrier_doppler(sat),
        sampling_frequency,
        get_carrier_phase_2pi(sat),
        T(get_amplitude(sat.cn0, n0)),
    )
    system = sat.system
    code = gen_code!(
        sat.code,
        system,
        get_prn(sat),
        sampling_frequency,
        get_code_frequency(system) + get_code_doppler(sat),
        get_code_phase(sat)
    )
    sat.signal .*= code
end

"""
$(SIGNATURES)

Propagates the satellite state over time
"""
function propagate(
    sat::ConstantDopplerSatellite{S, T, D, E, C},
    num_samples,
    intermediate_frequency,
    sample_frequency,
    n0,
    rng
) where {S, T, D, E, C}
    system = sat.system
    carrier_phase = (intermediate_frequency + get_carrier_doppler(sat)) * num_samples /
        sample_frequency + get_carrier_phase_2pi(sat)
    modded_carrier_phase = mod(carrier_phase + 0.5, 1) - 0.5
    code_phase = (get_code_frequency(system) + get_code_doppler(sat)) * num_samples /
        sample_frequency + get_code_phase(sat)
    modded_code_phase = mod(code_phase, get_code_length(system) * get_secondary_code_length(system))
    Δt = num_samples / sample_frequency
    exists = propagate(sat.exists, Δt, rng)
    doa = propagate(sat.doa, Δt, rng)
    cn0 = propagate(sat.cn0, n0, Δt)
    ConstantDopplerSatellite{S, T, D, E, C}(
        system,
        sat.prn,
        sat.carrier_doppler,
        modded_carrier_phase,
        modded_code_phase,
        cn0,
        exists,
        doa,
        sat.code,
        sat.signal
    )
end

@inline get_carrier_doppler(sat::AbstractSatellite) = sat.carrier_doppler
@inline get_code_doppler(sat::AbstractSatellite) =
    sat.carrier_doppler * get_code_center_frequency_ratio(sat.system)
@inline get_carrier_phase(sat::AbstractSatellite) = 2π * sat.carrier_phase
@inline get_carrier_phase_2pi(sat::AbstractSatellite) = sat.carrier_phase
@inline get_code_phase(sat::AbstractSatellite) = sat.code_phase
@inline get_prn(sat::AbstractSatellite) = sat.prn
@inline get_gnss_system(sat::AbstractSatellite) = sat.system
@inline get_carrier_to_noise_density_ratio(sat::AbstractSatellite) = sat.cn0

"""
$(SIGNATURES)

Calculate amplitude of a signal based on the carrier-to-noise-density-ratio (CN0) `cn0` in
[dB-Hz] and the noise density `n0` in [1/Hz].
"""
function calc_amplitude_from_cn0(cn0, n0)
    sqrt(linear(cn0) * n0)
end

"""
$(SIGNATURES)

Calculate amplitude of a signal based on the carrier-to-noise-density-ratio (CN0) `cn0` in
[dB-Hz], the noise density `n0` in [1/Hz] and steering vector influence.
"""
function get_amplitude_with_ant_influence(cn0, n0, steer_vec)
    get_amplitude(cn0, n0) * norm(steer_vec) / sqrt(length(steer_vec))
end

"""
$(SIGNATURES)

Approximates the distance between user and satellite based on the distance from earth center
to satellite `distance_from_earth_center` and the signal direction of arrival `enu_doa`.
"""
function calc_sat_user_distance(distance_from_earth_center, doa)
    init_doa = cart2sph(get_doa(doa))
    EARTH_RADIUS * cos(init_doa.ϕ + π / 2) +
        sqrt((EARTH_RADIUS * cos(init_doa.ϕ + π / 2))^2 - EARTH_RADIUS^2 +
        distance_from_earth_center^2)
end

@inline get_amplitude(sat::AbstractSatellite, n0) = get_amplitude(sat.cn0, n0)

"""
$(SIGNATURES)

Calculate carrier phase based on the distance between user and satellite `sat_user_distance`
and frequency `freq`.
"""
function calc_carrier_phase(sat_user_distance, freq)
    mod(freq * sat_user_distance / SPEED_OF_LIGHT + 1, 2) - 1
end

"""
$(SIGNATURES)

Calculate code phase based on the distance between satellite and user, the code frequency
`freq` and the code length `code_length`.
"""
function calc_code_phase(sat_user_distance, freq, code_length)
    mod(convert(Float64, (freq * sat_user_distance / SPEED_OF_LIGHT)), code_length)
end

"""
$(SIGNATURES)

Calculates a made up Doppler based on elevation and azimuth.
"""
function calc_doppler(distance_from_earth_center, doa, velocity, center_freq)
    doa_sph = cart2sph(get_doa(doa))
    center_freq / SPEED_OF_LIGHT * velocity * cos(π / 2 - asin(EARTH_RADIUS *
        sin(π / 2 + doa_sph.ϕ) / distance_from_earth_center)) * cos(doa_sph.θ)
end


text:
module DataSampling
  using ..DataParser
  files = ["ThresholdSampling","DataReader"]
  for f in files
    include("$f.jl")
  end

  function getDataSetIndicies(fn,file, size, isTest, options)
    return fn(file,size,isTest,options)
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary