text:
using Distributed
addprocs(10)

using Gadfly, Cairo, Compose, DataFrames, CSV, Colors, Viridis, HypothesisTests, ANOVA, GLM, RCall
@everywhere using Distributions, StatsBase, Distances, SharedArrays

function gen_colors(n)
    cs = distinguishable_colors(n,
        [colorant"#66c2a5", colorant"#fc8d62", colorant"#8da0cb", colorant"#e78ac3",
            colorant"#a6d854", colorant"#ffd92f", colorant"#e5c494", colorant"#b3b3b3"],
        lchoices=Float64[58, 45, 72.5, 90],
        transform=c->deuteranopic(c, 0.1),
        cchoices=Float64[20,40],
        hchoices=[75,51,35,120,180,210,270,310]
    )
    convert(Vector{Color}, cs)
end

function mixed_community(ϵ::Float64, ev_start::Bool, r_opinion::Float64, n_truthseekers::Int64, n_indifs::Int64, n_radicals::Int64, α::Float64, τ::Float64, averaging)
    n_normals = n_truthseekers + n_indifs
    st = ev_start == true ? collect(range(1/n_normals, stop=1, length=n_normals)) : rand(n_normals)
    δ = ev_start == true ? 10^-12 : 0.0
    rd = fill(r_opinion, n_radicals)
    bc_ar = vcat(st, rd)
    v = vcat(st, rd)
    eps = ϵ + δ
    j = 1
    @views while true
        for i in 1:n_normals
            fa = findall(@. abs(bc_ar[:, j] - bc_ar[i, j]) <= eps + δ)
            @inbounds v[i] = i <= n_truthseekers ? α*averaging(bc_ar[fa, j]) + (1 - α)*τ : averaging(bc_ar[fa, j])
        end
        bc_ar = hcat(bc_ar, v)
        if all(@. abs(bc_ar[:, j] - bc_ar[:, j + 1]) <= 10^-5)
            break
        end
        j += 1
    end
    df0 = DataFrame(bc_ar')
    names!(df0, [Symbol("$i") for i in 1:size(bc_ar, 1)])
    df1 = stack(df0)
    df1[:steps] = repeat(1:size(df0)[1], outer=size(bc_ar, 1))
    bb = Array{String,1}(undef, n_normals + n_radicals)
    bb[1:n_truthseekers] .= "normal"
    bb[n_truthseekers + 1:n_normals] .= "free rider"
    bb[n_normals + 1:n_normals + n_radicals] .= "radical"
    df1[:Agent] = repeat(bb, inner=size(bc_ar, 2))
    return df1
end

out_mixed = mixed_community(.25, false, .9, 25, 20, 10, .25, .65, mean)

plot(out_mixed, x=:steps, y=:value, group=:variable, color=:Agent, Geom.point, Geom.line,
    Coord.cartesian(xmax=maximum(out_mixed[:steps]) + 1),
    Guide.xlabel("Time"),
    Guide.ylabel("Opinion"),
    Scale.color_discrete_hue(gen_colors, levels = ["normal", "free rider", "radical"], preserve_order=true),
    yintercept=[0.65], Geom.hline(style=:dot, color=colorant"grey"),
    Guide.title("τ = .65, ρ = .9, ϵ = .25, α = .25"),
    Theme(panel_fill="white", point_size=1.25pt, line_width=1.5pt))

function bc_upd_rad_with_cd(ϵ_start::Float64, ev_start::Bool, r_opinion::Float64, n_normals::Int64, n_radicals::Int64, averaging_op, averaging_eps)
    st = ev_start == true ? collect(range(1/n_normals, stop=1, length=n_normals)) : rand(n_normals)
    δ = ev_start == true ? 10^-12 : 0.0
    rd = fill(r_opinion, n_radicals)
    bc_ar = vcat(rd, st)
    v = vcat(rd, st)
    eps_rad = fill(0.0, n_radicals)
    eps_nrm = rand(Uniform(0, ϵ_start), n_normals)
    eps = vcat(eps_rad, eps_nrm)
    tmp = zeros(Float64, n_normals)
    j = 1
    @views while true
        for i in (n_radicals + 1):(n_radicals + n_normals)
            fa = findall(@. abs(bc_ar[:, j] - bc_ar[i, j]) <= eps[i] + δ)
            v[i] = averaging_op(bc_ar[fa, j])
            tmp[i - n_radicals] = averaging_eps(eps[fa])
        end
        eps = vcat(eps_rad, tmp)
        bc_ar = hcat(bc_ar, v)
        if all(@. abs(bc_ar[:, j] - bc_ar[:, j + 1]) <= 10^-5)
            break
        end
        j += 1
    end
    return bc_ar
end

@everywhere function radicalized_normals_count_with_cd(ϵ_start::Float64, ev_start::Bool, r_opinion::Float64, n_normals::Int64, n_radicals::Int64, averaging_op, averaging_eps)
    st = ev_start == true ? collect(range(1/n_normals, stop=1, length=n_normals)) : rand(n_normals)
    δ = ev_start == true ? 10^-12 : 0.0
    rd = fill(r_opinion, n_radicals)
    bc_ar = vcat(rd, st)
    v = vcat(rd, st)
    eps_rad = fill(0.0, n_radicals)
    eps_nrm = fill(ϵ_start, n_normals)
    eps = vcat(eps_rad, eps_nrm)
    tmp = zeros(Float64, n_normals)
    j = 1
    @views while true
        @inbounds for i in (n_radicals + 1):(n_radicals + n_normals)
            fa = findall(@. abs(bc_ar[:, j] - bc_ar[i, j]) <= eps[i] + δ)
            v[i] = averaging_op(bc_ar[fa, j])
            tmp[i - n_radicals] = averaging_eps(eps[fa])
        end
        eps = vcat(eps_rad, tmp)
        bc_ar = hcat(bc_ar, v)
        if all(@. abs(bc_ar[:, j] - bc_ar[:, j + 1]) <= 10^-5)
            break
        end
        j += 1
    end
    return sum(@. abs(@view(bc_ar[:, end]) - r_opinion) <= 10^-3) - n_radicals
end

function radical_count_parallel_with_cd(ev_start::Bool, r_opinion::Float64, n_normals::Int64, averaging_op, averaging_eps)
    res = Array{Int64,2}(undef, 50, 50)
    for j in 1:50
        @inbounds res[:, j] = pmap(i -> radicalized_normals_count_with_cd(i/100, ev_start, r_opinion, n_normals, j, averaging_op, averaging_eps), 1:50)
    end
    return res
end

out_cd = radical_count_parallel_with_cd(true, 0.8, 50, mean, mean)

function xlabelname(x)
    n = x/100
    return "$n"
end

function ylabelname(x)
    n = abs(x - 50)
    return "$n"
end

ticks = collect(0:10:50)

Gadfly.spy(rotl90(out_cd),
    Guide.ColorKey(title="Radicalized\nnormals"),
    Guide.xlabel("ϵ"), Guide.ylabel("Radicals"),
    Guide.xticks(ticks=ticks), Scale.x_continuous(labels=xlabelname),
    Guide.yticks(ticks=ticks), Scale.y_continuous(labels=ylabelname),
    Guide.title("ρ = 0.8"),
    Scale.color_continuous(minvalue=0, maxvalue=50, colormap=Viridis.viridis),
    Theme(grid_color=colorant"white", panel_fill="white"))

@everywhere function radicalized_average_cd(ϵ_start::Float64, r_opinion::Float64, n_normals::Int64, n_radicals::Int64, averaging_op, averaging_eps)
    st = rand(n_normals)
    rd = fill(r_opinion, n_radicals)
    rd = fill(r_opinion, n_radicals)
    bc_ar = vcat(rd, st)
    v = vcat(rd, st)
    eps_rad = fill(0.0, n_radicals)
    eps_nrm = fill(ϵ_start, n_normals)
    eps = vcat(eps_rad, eps_nrm)
    tmp = zeros(Float64, n_normals)
    j = 1
    @views while true
        @inbounds for i in (n_radicals + 1):(n_radicals + n_normals)
            fa = findall(@. abs(bc_ar[:, j] - bc_ar[i, j]) <= eps[i])
            v[i] = averaging_op(bc_ar[fa, j])
            tmp[i - n_radicals] = averaging_eps(eps[fa])
        end
        eps = vcat(eps_rad, tmp)
        bc_ar = hcat(bc_ar, v)
        if all(@. abs(bc_ar[:, j] - bc_ar[:, j + 1]) <= 10^-5)
            break
        end
        j += 1
    end
    return sum(@. abs(@view(bc_ar[:, end]) - r_opinion) <= 10^-3) - n_radicals
end

function radical_average_parallel_cd(r_opinion::Float64, n_normals::Int64, averaging_op, averaging_eps)
    res = Array{Int64,2}(undef, 50, 50)
    for j in 1:50
        @inbounds res[:, j] = pmap(i -> radicalized_average_cd(i/100, r_opinion, n_normals, j, averaging_op, averaging_eps), 1:50)
    end
    return res
end

function rand_sim_cd(r_opinion::Float64, n_normals::Int64, n_sim::Int64, averaging_op, averaging_eps)
    rand_res = Array{Int64,3}(undef, 50, 50, n_sim)
    for i in 1:n_sim
        @inbounds rand_res[:, :, i] = radical_average_parallel_cd(r_opinion, n_normals, averaging_op, averaging_eps)
    end
    return rand_res
end

rand_res_cd = rand_sim_cd(1.0, 50, 100, mean, mean)

out_av_cd = mean(rand_res_cd, dims=3)

Gadfly.spy(rotl90(out_av_cd[:, :, 1]),
    Guide.ColorKey(title="Radicalized\nnormals"),
    Guide.xlabel("ϵ"), Guide.ylabel("Radicals"),
    Guide.xticks(ticks=ticks), Scale.x_continuous(labels=xlabelname),
    Guide.yticks(ticks=ticks), Scale.y_continuous(labels=ylabelname),
    Scale.color_continuous(minvalue=0, maxvalue=50, colormap=Viridis.viridis),
    Guide.title("ρ = 1.0"),
    Theme(grid_color=colorant"white", panel_fill="white"))

rand_res_cd = rand_sim_cd(.8, 50, 100, mean, mean)

out_av_cd = mean(rand_res_cd, dims=3)

Gadfly.spy(rotl90(out_av_cd[:, :, 1]),
    Guide.ColorKey(title="Radicalized\nnormals"),
    Guide.xlabel("ϵ"), Guide.ylabel("Radicals"),
    Guide.xticks(ticks=ticks), Scale.x_continuous(labels=xlabelname),
    Guide.yticks(ticks=ticks), Scale.y_continuous(labels=ylabelname),
    Scale.color_continuous(minvalue=0, maxvalue=50, colormap=Viridis.viridis),
    Guide.title("ρ = 0.8"),
    Theme(grid_color=colorant"white", panel_fill="white"))

function mixed_community_cd(ϵ_start::Float64, ev_start::Bool, r_opinion::Float64, n_truthseekers::Int64, n_indifs::Int64, n_radicals::Int64, α::Float64, τ::Float64, averaging_op, averaging_eps)
    n_normals = n_truthseekers + n_indifs
    st = ev_start == true ? collect(range(1/n_normals, stop=1, length=n_normals)) : rand(n_normals)
    δ = ev_start == true ? 10^-12 : 0.0
    rd = fill(r_opinion, n_radicals)
    bc_ar = vcat(st, rd)
    v = vcat(st, rd)
    eps_rad = fill(0.0, n_radicals)
    eps_nrm = rand(Uniform(0, ϵ_start), n_normals)
    eps = vcat(eps_nrm, eps_rad)
    tmp = zeros(Float64, n_normals)
    j = 1
    @views while true
        for i in 1:n_normals
            fa = findall(@. abs(bc_ar[:, j] - bc_ar[i, j]) <= eps[i] + δ)
            @inbounds v[i] = i <= n_truthseekers ? α*averaging_op(bc_ar[fa, j]) + (1 - α)*τ : averaging_op(bc_ar[fa, j])
            tmp[i] = averaging_eps(eps[fa])
        end
        eps = vcat(tmp, eps_rad)
        bc_ar = hcat(bc_ar, v)
        if all(@. abs(bc_ar[:, j] - bc_ar[:, j + 1]) <= 10^-5)
            break
        end
        j += 1
    end
    df0 = DataFrame(bc_ar')
    names!(df0, [Symbol("$i") for i in 1:size(bc_ar, 1)])
    df1 = stack(df0)
    df1[:steps] = repeat(1:size(df0)[1], outer=size(bc_ar, 1))
    bb = Array{String,1}(undef, n_normals + n_radicals)
    bb[1:n_truthseekers] .= "normal"
    bb[n_truthseekers + 1:n_normals] .= "free rider"
    bb[n_normals + 1:n_normals + n_radicals] .= "radical"
    df1[:Agent] = repeat(bb, inner=size(bc_ar, 2))
    return df1
end

out_mixed_cd = mixed_community_cd(.5, false, .9, 25, 20, 10, .25, .65, mean, mean)

plot(out_mixed_cd, x=:steps, y=:value, group=:variable, color=:Agent, Geom.point, Geom.line,
    Coord.cartesian(xmax=maximum(out_mixed_cd[:steps]) + 1),
    Guide.xlabel("Time"),
    Guide.ylabel("Opinion"),
    Scale.color_discrete_hue(gen_colors, levels = ["normal", "free rider", "radical"], preserve_order=true),
    yintercept=[0.65], Geom.hline(style=:dot, color=colorant"grey"),
    Guide.title("τ = .65, ρ = .9, ϵ⁰ ∼ U(0, .5), α = .25"),
    Theme(panel_fill="white", point_size=1.25pt, line_width=1.5pt))

out_mixed = mixed_community(.25, false, .0, 50, 0, 0, .75, .7, mean)
out_mixed_cd = mixed_community_cd(.5, false, .0, 50, 0, 0, .75, .7, mean, mean)

plot(out_mixed, x=:steps, y=:value, color=:variable, Geom.point, Geom.line,
    Coord.cartesian(xmax=maximum(out_mixed[:steps]) + 1),
    Guide.xlabel("Time"),
    Guide.ylabel("Opinion"),
    yintercept=[.7], Geom.hline(style=:dot, color=colorant"grey"),
    Guide.title("τ = .7, ϵ = .25, α = .75"),
    Theme(key_position=:none, panel_fill="white", point_size=1.25pt, line_width=1.5pt))

plot(out_mixed_cd, x=:steps, y=:value, color=:variable, Geom.point, Geom.line,
    Coord.cartesian(xmax=maximum(out_mixed_cd[:steps]) + 1),
    Guide.xlabel("Time"),
    Guide.ylabel("Opinion"),
    yintercept=[.7], Geom.hline(style=:dot, color=colorant"grey"),
    Guide.title("τ = .7, ϵ⁰ ∼ U(0, .5), α = .75"),
    Theme(key_position=:none, panel_fill="white", point_size=1.25pt, line_width=1.5pt))

@everywhere const n_steps = 50

@everywhere function bc_upd_ri_cd(ϵ_start::Float64, α::Float64, τ::Float64, r_opinion::Float64, n_truthseekers::Int64, n_indifs::Int64, n_radicals::Int64, averaging)
    n_agents = n_truthseekers + n_indifs
    rd = fill(r_opinion, n_radicals)
    bc_ar = Array{Float64,2}(undef, n_agents, n_steps + 1)
    bc_ar[:, 1] = rand(n_agents)
    eps_rad = fill(0.0, n_radicals)
    eps_nrm = fill(ϵ_start, n_agents)
    eps = vcat(eps_nrm, eps_rad)
    tmp = zeros(Float64, n_agents)
    @views @inbounds for j in 1:n_steps
        for i in 1:n_agents
            cc = vcat(bc_ar[:, j], rd)
            fa = findall(@. abs(cc - cc[i]) <= eps[i])
            bc_ar[i, j + 1] = i <= n_truthseekers ? α*averaging(cc[fa]) + (1 - α)*τ : averaging(cc[fa])
            tmp[i] = mean(eps[fa])
        end
        eps = vcat(tmp, eps_rad)
    end
    return bc_ar
end

function cost_sim_ri_cd(ϵ_start::Float64, α::Float64, τ::Float64, r_opinion::Float64, n_truthseekers::Int64, n_indifs::Int64, n_radicals::Int64, averaging, n_sim::Int64)
    v= SharedArray{Float64,1}(n_sim)
    @inbounds @sync @distributed for i in 1:n_sim
        res_rad = bc_upd_ri_cd(ϵ_start, α, τ, r_opinion, n_truthseekers, n_indifs, n_radicals, averaging)
        v[i] = sum(map(i->(i-τ)^2, res_rad[1:n_truthseekers, :]))
    end
    return v
end

τ = .1
ϵ = .3
α = .75
ρ1 = .3
ρ2 = .5
ρ3 = .7
nsim = 1000
without = cost_sim_ri_cd(ϵ, α, τ, ρ1, 50, 0, 0, mean, nsim)
withNarrow = cost_sim_ri_cd(ϵ, α, τ, ρ1, 50, 0, 25, mean, nsim)
withModerate = cost_sim_ri_cd(ϵ, α, τ, ρ2, 50, 0, 25, mean, nsim)
withWide = cost_sim_ri_cd(ϵ, α, τ, ρ3, 50, 0, 25, mean, nsim)

data = vcat(without, withNarrow, withModerate, withWide)
groups = categorical(repeat(1:4, inner=nsim))
df_aov = DataFrame(data = data, groups = groups)
model = fit(LinearModel, @formula(data ~ groups), df_aov, contrasts = Dict(:groups => EffectsCoding()))
anova(model)

mns = DataFrame(WO = without[:], WN = withNarrow[:], WM = withModerate[:], WW = withWide[:])

plot(mns, x=Col.index, y=Col.value,
    Geom.boxplot,
    Scale.x_discrete(levels=names(mns)),
    Guide.xlabel("Condition"),
    Guide.ylabel("SSE"),
    Guide.title("τ = $τ, ρ = $ρ1/$ρ2/$ρ3 (WN/WM/WW), ϵ = $ϵ, α = $α"),
    Theme(default_color=colorant"#66c2a5", panel_fill="white"))

@rput df_aov

R"""
m <- aov(data ~ groups, data=df_aov)
summary(m)
"""

R"TukeyHSD(m)"

R"library(lsr)"

R"etaSquared(m)"

mean(without), mean(withNarrow), mean(withModerate), mean(withWide)

std(without), std(withNarrow), std(withModerate), std(withWide)

τ = .1
ϵ = .2
α = .5
ρ = .3
nsim = 1000
without = cost_sim_ri_cd(ϵ, α, τ, ρ, 50, 0, 15, mean, nsim)
with10 = cost_sim_ri_cd(ϵ, α, τ, ρ, 50, 10, 15, mean, nsim)
with25 = cost_sim_ri_cd(ϵ, α, τ, ρ, 50, 25, 15, mean, nsim)
with50 = cost_sim_ri_cd(ϵ, α, τ, ρ, 50, 50, 15, mean, nsim);

data = vcat(without, with10, with25, with50)
groups = categorical(repeat(1:4, inner=nsim))
df_aov = DataFrame(data = data, groups = groups)
model = fit(LinearModel, @formula(data ~ groups), df_aov, contrasts = Dict(:groups => EffectsCoding()))
anova(model)

mns = DataFrame(X= without[:], Y = with10[:], Z = with25[:], W = with50[:])

rename!(mns, :X => Symbol("0"), :Y => Symbol("10"), :Z => Symbol("25"), :W => Symbol("50"))

plot(mns, x=Col.index, y=Col.value,
    Geom.boxplot,
    Scale.x_discrete(levels=names(mns)),
    Guide.xlabel("Free riders"),
    Guide.ylabel("SSE"),
    Guide.title("τ = $τ, ρ = $ρ, ϵ = $ϵ, α = $α"),
    Theme(default_color=colorant"#66c2a5", panel_fill="white"))

@rput df_aov

R"""
m <- aov(data ~ groups, data=df_aov)
summary(m)
"""

R"TukeyHSD(m)"

R"etaSquared(m)"

mean(without), mean(with10), mean(with25), mean(with50)

std(without), std(with10), std(with25), std(with50)


text:
using DataFrames
using FastAI
using Test

include("../src/transforms.jl")

macro excludeTest(func)
    println("$func excluded from test run")
    return nothing
end

macro includeTest(func)
    return func
end

function test_transforms_process_files()
    println("starting tests for process_files()...")
    @test process_files(".", [".dataset.jl", "test.jl"], ".jl") == ["./test.jl"]
    @test process_files(".", ["dataset.jl", "test.jl"], ".jl") == ["./dataset.jl", "./test.jl"]
    @test process_files(".", ["dataset.jl", "test.jl"]) == ["./dataset.jl", "./test.jl"]
    println("..End of tests for process_files(). Printing Summary...")
end

function test_transforms_get_files()
    println("starting tests for process_files()...")
    path = "test/test_folder"
    extension = ".txt"
    @test get_files(path,
                    extension) == ["test/test_folder/test1/test1-1.txt",
                                  "test/test_folder/test1/test1-2.txt",
                                  "test/test_folder/test2/test2-1.txt",
                                  "test/test_folder/test2/test2-2.txt"]
    println("..End of tests for process_files(). Printing Summary...")

    @test get_files(path) == ["test/test_folder/test1/test1-1.txt",
                              "test/test_folder/test1/test1-2.txt",
                              "test/test_folder/test1/test1-3.yml",
                              "test/test_folder/test2/test2-1.txt",
                              "test/test_folder/test2/test2-2.txt"]

    @test get_files(path,".yml") == ["test/test_folder/test1/test1-3.yml"]
end

function test_Filegetter()
    filegetter= FileGetter("test",".yml")
    filegetter2 = FileGetter("",".yml")
    @test filegetter("test_folder") == ["test/test_folder/test1/test1-3.yml"]
    @test_throws ErrorException filegetter2("") == "A path must be provided"
end

function test_get_image_files()
    @test get_image_files("test/test_img_folder") == ["test/test_img_folder/test.jpg"]
end

function test_Imagegetter()
    imagegetter= FileGetter("test",".jpg")
    imagegetter2 = FileGetter("",".jpg")
    @test imagegetter("test_img_folder") == ["test/test_img_folder/test.jpg"]
    @test_throws ErrorException imagegetter2("") == "A path must be provided"
end

function test_get_text_files()
    @test get_text_files("test/test_folder") == ["test/test_folder/test1/test1-1.txt",
                                                 "test/test_folder/test1/test1-2.txt",
                                                 "test/test_folder/test2/test2-1.txt",
                                                 "test/test_folder/test2/test2-2.txt"]
end

function test_Itemgetter()
    data= [(5, 4), (4, 6)]
    @test ItemGetter(1, data) == (5,4)

end

function test_RandomSplitter()
    X = rand(2, 6)
    split = RandomSplitter()
    train, test = split(X)
    @test length(train) == length(X)-length(test)
end

function test_TrainTestSplitter()
    y = [:a, :a, :a, :a, :a, :a, :b, :b, :b]
    split = TrainTestSplitter(y)
    train1,test1 = split(0.7)
    @test length(train1) == length(y)-length(test1)
end

function test_TrainTestSplitter2()
    X = [1 0; 1 0; 1 0; 1 0; 0 1; 0 1]
    split = TrainTestSplitter(argmax, X)
    train1,test1 = split(0.7)
    @test length(train1) == length(X)-length(test1)
end

function test_IndexSplitter()
    a,b = IndexSplitter([1,2,3,4,5], [3,5])
    @test a == [1,2,4]
    @test b == [3,5]
end

function test_grandparent_idxs()
    files = ["test/test_folder1/test2/test1-1.txt",
      "test/test_folder2/test1/test1-2.txt",
      "test/test_folder/test2/test2-1.txt",
      "test/test_folder/test2/test2-2.txt"]

    @test grandparent_idxs(files, "test2") == [1,3,4]
    @test grandparent_idxs(files, ("test1", "test2")) == [2,1,3,4]

end

function test_ColSplitter()
    data = DataFrame(a = [1, 2, 3], b = [true, false, true])
    split_with_name = ColSplitter(col = :b)
    split_with_index = ColSplitter()
    a,b = split_with_name(data)
    @test a == [2]
    @test b == [1,3]
end

function test_assertBounds()
    @test assertBounds(0.3, 0,1) === nothing
    @test_throws ErrorException assertBounds(24,0,1)
end

function test_RandomSubsetSplitter()
    items = Array(1:100)
    splits = RandomSubsetSplitter(0.3, 0.1)
    trainidxs, valididxs = splits(items)
    @test length(trainidxs) === 30
    @test length(valididxs) === 11
end

function test_parent_label()
    fnames = ["fastai_dev/dev/data/mnist_tiny/train/3/9932.png",
               "fastai_dev/dev/data/mnist_tiny/valid/4/9932.png"]
    @test parent_label("fastai_dev/dev/data/mnist_tiny/train/3/9932.png") === '3'
    @test [parent_label(path) for path in fnames] == ['3', '4']
end

function test_RegexLabeller()
    path= "fastai_dev/dev/data/mnist_tiny/train/3/9932.png"
    pat = r"/[\d]+/"
    split = RegexLabeller(pat)
    @test split(path) == "3"
end

#todo: separate these tests
function test_ColReader()
    df = DataFrame(A = ["a", "b", "c"], B = ["1 2", "0", ""], C=["a", "b", "c"])
    colreader1 = ColReader(:A; pref="1", suff="2")
    @test colreader1(df) == ["1a2", "1b2", "1c2"]
    delimsplitter = ColReader(:B, " ")
    @test delimsplitter(df) == [["1", "2"], ["0"], [""]]
    colreader2 = ColReader([:A, :C], pref="1", suff="2")
    @test colreader2(df) == [["1a2", "1a2"], ["1b2", "1b2"], ["1c2", "1c2"]]
    colreader3 = ColReader(1, pref="1", suff="2")
    @test colreader3(df) == ["1a2", "1b2", "1c2"]
    colreader4 = ColReader(1)
    @test colreader4(df) == ["a", "b", "c"]
end

function test_CategoryMap_default()
    df = DataFrame(A = ["d", "c", "b", "d"])
    map = CategoryMap(df[!, :A])
    @test map.items == ["b", "c", "d"]
    @test map.o2i == Dict("b" => 1, "c" => 2, "d" => 3)
end

function test_CategoryMap_add_na()
    df = DataFrame(A = ["d", "c", "b", "d"])
    map = CategoryMap(df[!, :A], add_na=true)
    @test map.items == ["#na#", "b", "c", "d"]
    @test map.o2i == Dict("#na#" => 1, "b" => 2, "c" => 3, "d" => 4)
end

function test_CategoryMap_categorical()
    df = DataFrame(A = ["d", "c", "b", "d"])
    col = categorical(df[!, :A])
    map = CategoryMap(col, add_na=true)
    @test map.items == ["#na#", "b", "c", "d"]
    @test map.o2i == Dict("#na#" => 1, "b" => 2, "c" => 3, "d" => 4)
end

function test_Categorize()
    data = ["cat", "dog", "cat"]
    cat = Categorize(data)
    @test encode(cat) == ["cat", "dog"]
    @test decode(cat, "cat") == 1
    @test decode("cat") == "cat"
    @test decode(cat, 2) == "dog"
end

 @testset "All" begin
    @excludeTest test_transforms_process_files()
    @excludeTest test_transforms_get_files()
    @excludeTest test_Filegetter()
    @excludeTest test_get_image_files()
    @excludeTest test_Imagegetter()
    @excludeTest test_get_text_files()
    @excludeTest test_Itemgetter()
    @excludeTest test_RandomSplitter()
    @excludeTest test_TrainTestSplitter2()
    @excludeTest test_IndexSplitter()
    @excludeTest test_grandparent_idxs()
    @excludeTest test_ColSplitter()
    @excludeTest test_assertBounds()
    @excludeTest test_RandomSubsetSplitter()
    @excludeTest test_parent_label()
    @excludeTest test_RegexLabeller()
    @excludeTest test_ColReader()
    @excludeTest test_CategoryMap_default()
    @excludeTest test_CategoryMap_add_na()
    @excludeTest test_CategoryMap_categorical()
    @includeTest test_Categorize()
end

ext: jl


text:
using IRTools, Test
using IRTools: xcall, func, argument!, @dynamo, IR, insertafter!, return!, recurse!

foo(::typeof(*), a, b) = a+b
@dynamo function foo(a...)
  ir = IR(a...)
  ir == nothing && return
  recurse!(ir)
  return ir
end

dummy() = nothing
ir = @code_ir dummy()
a = argument!(ir)
b = argument!(ir)
t = pushfirst!(ir, xcall(:tuple, a, b))
op = insertafter!(ir, t, xcall(:tuple, GlobalRef(Base, :*)))
r = insertafter!(ir, op, xcall(Core, :_apply, GlobalRef(Base, :*), t))
return!(ir, r)
foofunc = func(ir)
@assert foofunc(nothing, 0, 1) == 0
@test foo(foofunc, nothing, 0, 1) == 1

# + is called inside Core._apply_iterate(Base.iterate, (+), ...)
g(args...) = Base.:*(args...)
@assert g(0, 1) == 0
@test foo(g, 0, 1) == 1


text:
# Inject context information to the logger via TransformerLogger

using ContextTracking
using Logging
using LoggingExtras

# Top level function
function foo()
    @memo x = 1
    @memo y = 2
    bar()
end

# Lower level function that inherits parent context
function bar()
    z = 3
    @info "inside bar" z
end

# Create a transformer logger that appends context data to the existing kwargs
context_logger = TransformerLogger(current_logger()) do log
    updated_kwargs = (context().data..., log.kwargs...)
    merge(log, (; kwargs = updated_kwargs))
end

# Test!
with_logger(context_logger) do
    foo()
end

#=
julia> with_logger(context_logger) do
           foo()
       end
┌ Info: inside bar
│   context.x = 1
└   y = 2
=#


text:
ENV["OMP_NUM_THREADS"] = 2

using Test
using Pardiso
using Random
using SparseArrays
using LinearAlgebra

Random.seed!(1234)

psolvers = [MKLPardisoSolver]
if Pardiso.PARDISO_LOADED[]
    push!(psolvers, PardisoSolver)
else
    @warn "Not testing project Pardiso solver"
end

@show Pardiso.MklInt

println("Testing ", psolvers)

# Test solver + for real and complex data
@testset "solving" begin
for pardiso_type in psolvers
    ps = pardiso_type()
    for T in (Float64, ComplexF64)
        A1 = sparse(rand(T, 10,10))
        for B in (rand(T, 10, 2), view(rand(T, 10, 4), 1:10, 2:3))
            X = similar(B)
            # Test unsymmetric, herm indef, herm posdef and symmetric
            for A in SparseMatrixCSC[A1, A1 + A1', A1'A1, transpose(A1) + A1]
                solve!(ps, X, A, B)
                @test X ≈ A\Matrix(B)

                X = solve(ps, A, B)
                @test X ≈ A\Matrix(B)

                solve!(ps, X, A, B, :C)
                @test X ≈ A'\Matrix(B)

                X = solve(ps, A, B, :C)
                @test X ≈ A'\Matrix(B)

                solve!(ps, X, A, B, :T)
                @test X ≈ copy(transpose(A))\Matrix(B)

                X = solve(ps, A, B, :T)
                @test X ≈ copy(transpose(A))\Matrix(B)
            end
        end
    end
end
end #testset

module ExampleUnsym
    include("../examples/exampleunsym.jl")
end
module ExampleSym
    include("../examples/examplesym.jl")
end
module ExampleHerm
    include("../examples/exampleherm.jl")
end

if Sys.CPU_THREADS >= 4
    @testset "procs" begin
        ps = MKLPardisoSolver()
        np = get_nprocs(ps)
        set_nprocs!(ps, 2)
        @test get_nprocs(ps) == 2
        set_nprocs!(ps, np)
        @test get_nprocs(ps) == np
    end
end

if Pardiso.PARDISO_LOADED[]
@testset "schur" begin
    # reproduce example from Pardiso website
    include("schur_matrix_def.jl")
    @test norm(real(D) - real(C)*rA⁻¹*real(B) - s) < 1e-10*(8)^2
    # @test norm(D - C*A⁻¹*B - S) < 1e-10*(8)^2

    # try some random matrices
    m = 50; n = 15; p = .1
    ps = PardisoSolver()
    for T in (Float64, )#ComplexF64)
        ps = PardisoSolver()
        pardisoinit(ps)
        if T == Float64
            set_matrixtype!(ps, 11)
        else
            set_matrixtype!(ps, 13)
        end
        for j ∈ 1:100
            A = 5I + sprand(T,m,m,p)
            A⁻¹ = inv(Matrix(A))
            B = sprand(T,m,n,p)
            C = sprand(T,n,m,p)
            D = 5I + sprand(T,n,n,p)
            M = [A B; C D]

            # test integer block specification
            S = schur_complement(ps, M, n);
            @test norm(D - C*A⁻¹*B - S) < 1e-10*(m+n)^2

            # test sparse vector block specification
            x = spzeros(T,m+n)
            x[(m+1):(m+n)] .= 1
            S = schur_complement(ps, M, x);
            @test norm(D - C*A⁻¹*B - S) < 1e-10*(m+n)^2

            # test sparse matrix block specification
            x = spzeros(T,m+n,2)
            x[(m+1):(m+n-1),1] .= 1
            x[end,2] = 1
            S = schur_complement(ps, M, x);
            @test norm(D - C*A⁻¹*B - S) < 1e-10*(m+n)^2
        end
    end
end # testset
end

@testset "error checks" begin
for pardiso_type in psolvers

    ps = pardiso_type()

    A = sparse(rand(10,10))
    B = rand(10, 2)
    X = rand(10, 2)

    if pardiso_type == PardisoSolver
        printstats(ps, A, B)
        checkmatrix(ps, A)
        checkvec(ps, B)
    end

    set_matrixtype!(ps, 13)
    @test_throws ErrorException pardiso(ps, X, A, B)
    @test_throws ArgumentError solve(ps, A, B, :P)
    @test_throws ArgumentError solve!(ps, X, A, B, :P)

    set_matrixtype!(ps, 11)
    X = zeros(12, 2)
    @test_throws DimensionMismatch solve!(ps,X, A, B)

    B = rand(12, 2)
    @test_throws DimensionMismatch solve(ps, A, B)
end
end # testset

@testset "getters and setters" begin
for pardiso_type in psolvers
    ps = pardiso_type()
    set_iparm!(ps, 1, 0)
    pardisoinit(ps)
    @test get_iparm(ps, 1) == 1

    @test_throws ArgumentError set_phase!(ps, 5)
    @test_throws ArgumentError set_msglvl!(ps, 2)
    @test_throws ArgumentError set_matrixtype!(ps, 15)

    if pardiso_type == PardisoSolver
        @test_throws ArgumentError set_solver!(ps, 2)

        set_dparm!(ps, 5, 13.37)
        @test get_dparm(ps, 5) == 13.37

        set_solver!(ps, 1)
        @test Int(get_solver(ps)) == 1
    end

    set_iparm!(ps, 13, 100)
    @test get_iparm(ps, 13) == 100

    set_matrixtype!(ps, Pardiso.REAL_SYM)
    @test get_matrixtype(ps) == Pardiso.REAL_SYM

    set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT)
    @test get_phase(ps) == Pardiso.ANALYSIS_NUM_FACT

    set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
    @test get_msglvl(ps) == Pardiso.MESSAGE_LEVEL_ON
end

@testset "pardiso" begin
    for pardiso_type in psolvers
        A = sparse(rand(2,2) + im * rand(2,2))
        b = rand(2)          + im * rand(2)
        ps = pardiso_type()
        set_matrixtype!(ps, Pardiso.COMPLEX_NONSYM)
        x = Pardiso.solve(ps, A, b);
        set_phase!(ps, Pardiso.RELEASE_ALL)
        pardiso(ps)
    end
end
end # testset


text:
"""
FinEtools (C) 2017-2019, Petr Krysl

Finite Element tools.  Julia implementation  of the finite element method
for continuum mechanics.
"""
module FinEtools

__precompile__(true)

include("allmodules.jl")

# Exports follow:

###########################################################################
# General facilities
###########################################################################
using .FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
# Exported: basic numerical types, type of data dictionary
export FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict

using .BoxModule: inbox, initbox!, updatebox!, boundingbox, inflatebox!, boxesoverlap, intersectboxes
# Exported: methods for manipulating and testing boxes
export inbox, initbox!, updatebox!, boundingbox, inflatebox!, boxesoverlap, intersectboxes

using .PhysicalUnitModule: phun
# Exported: function for evaluating physical units
export phun

using .RotationUtilModule: rotmat3!, rotmat3, skewmat!, cross3!, cross2
# Exported: functions for 3D rotation matrix computation, skew symmetric matrix computation,  in-place  cross product of 3-vectors, and  cross product of 2-vectors
export rotmat3!, rotmat3, skewmat!, cross3!, cross2

using .CSysModule: CSys, updatecsmat!
# Exported: type  for coordinate systems, methods to invoke the update callback
export CSys, updatecsmat!

using .FESetModule: AbstractFESet,  AbstractFESet0Manifold,  AbstractFESet1Manifold,  AbstractFESet2Manifold,  AbstractFESet3Manifold, manifdim, nodesperelem, count, delegateof, accepttodelegate, fromarray!, connasarray, setlabel!, subset, cat, updateconn!, bfun, bfundpar, map2parametric, inparametric, centroidparametric,  FESetP1, FESetL2, FESetL3, FESetT3, FESetQ4, FESetQ9, FESetQ8, FESetT6, FESetH8, FESetH20, FESetH27, FESetT4, FESetT10
# Exported: types of finite element sets, abstract and concrete
export AbstractFESet,  AbstractFESet0Manifold,  AbstractFESet1Manifold,  AbstractFESet2Manifold,  AbstractFESet3Manifold, FESetP1, FESetL2, FESetL3, FESetT3, FESetQ4, FESetQ9, FESetQ8, FESetT6, FESetH8, FESetH20, FESetH27, FESetT4, FESetT10
# Exported: methods for accessing dimensions and counts
export manifdim, nodesperelem, count
# Exported: methods for accepting delegation and revealing the delegating object
export delegateof, accepttodelegate
# Exported: methods for  manipulating connectivity  and labels
export fromarray!, connasarray, setlabel!, subset, cat, updateconn!
# Exported: methods for computing basis function values and derivatives of basis functions  with respect to the parametric coordinates, and methods for working with parametric coordinates
export bfun, bfundpar, map2parametric, inparametric, centroidparametric

using .FENodeSetModule: FENodeSet, spacedim, xyz3, count
# Exported: type for FE node sets, methods for accessing dimensions  and counts
export FENodeSet, spacedim, xyz3, count

using .FENodeToFEMapModule: FENodeToFEMap
# Exported: type for maps from nodes to finite elements
export FENodeToFEMap

using .FieldModule: AbstractField, ndofs,  nents, gathersysvec, gathersysvec!, gathervalues_asvec!, gathervalues_asmat!, gatherdofnums!, gatherfixedvalues_asvec!, gatherfixedvalues_asmat!,numberdofs!, setebc!, applyebc!, scattersysvec!, copyto!, wipe!, prescribeddofs, incrscattersysvec!
# Exported: abstract field type, methods for the abstract field type (retrieval of data  from a field, setting of data in the field)
export AbstractField, ndofs,  nents, gathersysvec, gathersysvec!, gathervalues_asvec!, gathervalues_asmat!, gatherdofnums!, gatherfixedvalues_asvec!, gatherfixedvalues_asmat!, numberdofs!, setebc!, applyebc!, scattersysvec!, copyto!, wipe!, prescribeddofs, incrscattersysvec!

using .GeneralFieldModule: GeneralField
# Exported: type of general field
export GeneralField

using .NodalFieldModule: NodalField, nnodes
# Exported: type of nodal field
export NodalField, nnodes

using .ElementalFieldModule: ElementalField, nelems
# Exported: type of elemental field
export ElementalField, nelems

using .MeshUtilModule: linearspace, gradedspace, logspace
# Exported: functions to generate a sequence of numbers between start and stop
export linearspace, gradedspace, logspace

using .MeshSelectionModule: connectednodes, connectedelems, selectnode, selectelem, findunconnnodes
# Exported: functions to select (find) nodes and elements
export connectednodes, connectedelems, selectnode, selectelem, findunconnnodes

using .MeshExportModule.VTK: vtkexportmesh
# Exported: VTK export
export vtkexportmesh
using .MeshExportModule.Abaqus: AbaqusExporter, close, HEADING, COMMENT, PART, END_PART,    ASSEMBLY, END_ASSEMBLY, INSTANCE, END_INSTANCE, NODE, ELEMENT,    NSET_NSET, ELSET_ELSET, ORIENTATION, MATERIAL, ELASTIC, EXPANSION, DENSITY, SECTION_CONTROLS, SOLID_SECTION, SURFACE_SECTION, STEP_PERTURBATION_STATIC, STEP_FREQUENCY,   STEP_PERTURBATION_BUCKLE, BOUNDARY, DLOAD, CLOAD, TEMPERATURE,    END_STEP,  NODE_PRINT, EL_PRINT,  ENERGY_PRINT
# Exported: Abaqus export
export AbaqusExporter, close, HEADING, COMMENT, PART, END_PART,    ASSEMBLY, END_ASSEMBLY, INSTANCE, END_INSTANCE, NODE, ELEMENT,    NSET_NSET, ELSET_ELSET, ORIENTATION, MATERIAL, ELASTIC, EXPANSION, DENSITY, SECTION_CONTROLS, SOLID_SECTION, SURFACE_SECTION, STEP_PERTURBATION_STATIC, STEP_FREQUENCY,   STEP_PERTURBATION_BUCKLE, BOUNDARY, DLOAD, CLOAD, TEMPERATURE,    END_STEP,  NODE_PRINT, EL_PRINT,  ENERGY_PRINT
using .MeshExportModule.CSV: savecsv
# Exported: simple CSV export
export savecsv
using .MeshExportModule.NASTRAN: NASTRANExporter, close, CEND, BEGIN_BULK, ENDDATA, GRID, PSOLID, MAT1, CTETRA
# Exported: NASTRAN  export
export NASTRANExporter, close, CEND, BEGIN_BULK, ENDDATA, GRID, PSOLID, MAT1, CTETRA
using .MeshExportModule.STL: STLExporter, solid, facet, endsolid
# Exported: STL export
export STLExporter, solid, facet, endsolid
using .MeshExportModule.H2Lib: h2libexporttri
# Exported: H2Lib export
export h2libexporttri
# Exported: MESH export
using .MeshExportModule.MESH
export MESH
# Exported: H5MESH export
using .MeshExportModule.H5MESH
export H5MESH

using .MeshModificationModule: meshboundary,  fusenodes,  compactnodes,  mergemeshes, mergenmeshes, mergenodes,  renumberconn!,  meshsmoothing, mirrormesh, nodepartitioning, interior2boundary
# Exported: extraction of boundary, fusing of nodes  and merging of meshes, mesh smoothing,  node partitioning
export  meshboundary,  fusenodes,  compactnodes,  mergemeshes, mergenmeshes, mergenodes,  renumberconn!,  meshsmoothing, mirrormesh, nodepartitioning, interior2boundary

using .MeshImportModule: import_NASTRAN, import_ABAQUS, import_MESH
# Exported: mesh import functions
export import_NASTRAN, import_ABAQUS, import_MESH, import_H5MESH

using .VectorCacheModule: VectorCache, updateretrieve!, settime!
# Exported: vector-cache type and methods to invoke the update callback
export VectorCache, updateretrieve!, settime!

using .ForceIntensityModule: ForceIntensity, updateforce!, settime!
# Exported: force-intensity type and methods to invoke the update callback
export ForceIntensity, updateforce!, settime!

using .SurfaceNormalModule: SurfaceNormal, updatenormal!
# Exported: surface-normal evaluator type and methods to invoke the update callback
export SurfaceNormal, updatenormal!

using .AssemblyModule: AbstractSysmatAssembler, SysmatAssemblerSparse, SysmatAssemblerSparseSymm, SysmatAssemblerSparseHRZLumpingSymm, startassembly!, assemble!, makematrix!, AbstractSysvecAssembler, SysvecAssembler, startassembly!, assemble!, makevector!, SysmatAssemblerReduced
# Exported: types and methods for  sparse matrix assembly  and vector assembly
export AbstractSysmatAssembler, SysmatAssemblerSparse, SysmatAssemblerSparseSymm, SysmatAssemblerSparseHRZLumpingSymm, startassembly!, assemble!, makematrix!, AbstractSysvecAssembler, SysvecAssembler, startassembly!, assemble!, makevector!, SysmatAssemblerReduced

using .IntegRuleModule: AbstractIntegRule, TriRule, GaussRule, TetRule, PointRule, SimplexRule, TrapezoidalRule, NodalSimplexRule, NodalTensorProductRule
# Exported: type for various integration rules
export AbstractIntegRule, TriRule, GaussRule, TetRule, PointRule, SimplexRule, TrapezoidalRule, NodalSimplexRule, NodalTensorProductRule

using .IntegDomainModule: IntegDomain, otherdimensionunity, Jacobianpoint, Jacobiancurve, Jacobiansurface, Jacobianvolume, Jacobianmdim, integrationdata
# Exported: type to handle  integration data for various manifold dimensions
export IntegDomain, otherdimensionunity, Jacobianpoint, Jacobiancurve, Jacobiansurface, Jacobianvolume, Jacobianmdim, integrationdata

using .FEMMBaseModule: AbstractFEMM, FEMMBase, associategeometry!, integratefieldfunction, integratefunction, transferfield!, distribloads, connectionmatrix, fieldfromintegpoints, elemfieldfromintegpoints, innerproduct, field_elem_to_nodal!, field_nodal_to_elem!
# Exported: type base discretization methods
export AbstractFEMM, FEMMBase, associategeometry!, integratefieldfunction, integratefunction, transferfield!, distribloads, connectionmatrix, fieldfromintegpoints, elemfieldfromintegpoints, innerproduct, field_elem_to_nodal!, field_nodal_to_elem!

###########################################################################
# Mesh-generation functionality for various shapes
###########################################################################
using .MeshQuadrilateralModule: Q4annulus, Q4quadrilateral, Q4elliphole, Q4block, Q4blockx, Q4refine, Q8block, Q4toQ8, Q8annulus, Q8blockx, Q4spheren, Q4circlen, Q4extrudeL2
# Exported: mesh generation functions for quadrilateral elements
export Q4annulus, Q4quadrilateral, Q4elliphole, Q4block, Q4blockx, Q4refine, Q8block, Q4toQ8, Q8annulus, Q8blockx, Q4spheren, Q4circlen, Q4extrudeL2

using .MeshLineModule: L2block, L2blockx, L3blockx
# Exported: mesh generation functions for line elements
export L2block, L2blockx, L3blockx

using .MeshTriangleModule: T3blockx, T3block,  T3toT6,  T6block,  Q4toT3,  T3refine, T6blockx, T3annulus, T6annulus, T3circlen, T3circleseg
# Exported: mesh generation functions for triangular elements
export  T3blockx, T3block,  T3toT6,  T6block,  Q4toT3,  T3refine, T6blockx, T3annulus, T6annulus, T3circlen, T3circleseg

using .MeshHexahedronModule: H8block,  H8blockx,  H8sphere,  H8refine, H8hexahedron, H8extrudeQ4, H8spheren, H8voximg,  H8layeredplatex, H8elliphole, H8toH27,  H27block, H20block,  H8toH20, H20blockx, H27blockx, H8cylindern, T4toH8
# Exported: mesh generation functions for hexahedral elements
export  H8block,  H8blockx,  H8sphere,  H8refine, H8hexahedron, H8extrudeQ4, H8spheren, H8voximg,  H8layeredplatex, H8elliphole, H8toH27,  H27block, H20block,  H8toH20, H20blockx, H27blockx, H8cylindern, T4toH8

using .MeshTetrahedronModule: T4block, T4blockx, T4toT10, T10toT4, T10block, T10blockx, T10layeredplatex, T4meshedges, T4voximg, T4refine, T10refine, T4refine20, T4quartercyln
# Exported: mesh generation functions for tetrahedral elements
export  T4block, T4blockx, T4toT10, T10toT4, T10block, T10blockx, T10layeredplatex, T4meshedges, T4voximg, T4refine, T10refine, T4refine20, T4quartercyln

###########################################################################
# Abstract material
###########################################################################
using .MatModule: AbstractMat, massdensity
# Exported: abstract type of material
export AbstractMat, massdensity

end


text:
using SparseIR
using Test

import PyCall: pyimport, PyNULL, PyVector
import LegendrePolynomials: Pl

@testset "test_augmented_bosonic_basis" begin
    """Augmented bosonic basis"""
    wmax = 2
    beta = 1000
    basis = FiniteTempBasis(boson, beta, wmax, 1e-6)
    basis_legg = LegendreBasis(boson, beta, 2)
    basis_comp = CompositeBasis([basis_legg, basis])

    # G(tau) = c - e^{-tau*pole}/(1 - e^{-beta*pole})
    pole = 1.0
    c = 1e-2
    tau_smpl = TauSampling(basis_comp)
    gtau = c .- exp.(-tau_smpl.sampling_points .* pole) ./ (1 - exp(-beta * pole))
    gl_from_tau = fit(tau_smpl, gtau)

    gtau_reconst = evaluate(tau_smpl, gl_from_tau)
    @test isapprox(gtau, gtau_reconst; atol=1e-14 * maximum(abs.(gtau)), rtol=0)
end


text:
using Plots, Test

@testset "Shorthands" begin
    @testset "Set Lims" begin
        p = plot(rand(10))

        xlims!((1,20))
        @test xlims(p) == (1,20)

        ylims!((-1,1))
        @test ylims(p) == (-1,1)

        zlims!((-1,1))
        @test zlims(p) == (-1,1)

        xlims!(-1,11)
        @test xlims(p) == (-1,11)

        ylims!((-10,10))
        @test ylims(p) == (-10,10)

        zlims!((-10,10))
        @test zlims(p) == (-10,10)
    end

    @testset "Set Ticks" begin
        p = plot([0,2,3,4,5,6,7,8,9,10])

        xticks = 2:6
        xticks!(xticks)
        @test Plots.get_subplot(current(),1).attr[:xaxis][:ticks] == xticks

        yticks = 0.2:0.1:0.7
        yticks!(yticks)
        @test Plots.get_subplot(current(),1).attr[:yaxis][:ticks] == yticks

        xticks = [5,6,7.5]
        xlabels = ["a","b","c"]

        xticks!(xticks, xlabels)
        @test Plots.get_subplot(current(),1).attr[:xaxis][:ticks] == (xticks, xlabels)

        yticks = [.5,.6,.75]
        ylabels = ["z","y","x"]
        yticks!(yticks, ylabels)
        @test Plots.get_subplot(current(),1).attr[:yaxis][:ticks] == (yticks, ylabels)
    end
end


text:
# Directions
const directions = [(-1, 0),(1, 0),(0, 1),(0, -1)]
# For readability
PICK_FOOD_1 = 5
PLACE_FOOD_1 = 6
# Perform a move action
function move(env::NatureEnv, player::Int, dir::Int)
    new_pos = env.players[player].pos .+ directions[dir]
    outofbounds(env, new_pos) && return
    env.players[player].pos = new_pos
end
# Perform a food action
function food(env::NatureEnv, p::Int, idx::Int)
    # [ pick1, place1, pick2, place2, ... pickf, placef]
    pick  = idx % 2 == 1
    place = idx % 2 == 0
    food_type = floor(Int, ((idx - 1) / 2)) + 1
    fframe   = env.food_frames[food_type]
    player   = env.players[p]
    num_food = fframe[player.pos...]
    place_amount = 0.5

    if pick && num_food > 0
        player.food_counts = player.food_counts .+ Tuple(num_food * onehot(food_type, env.food_types))
        for (placed_player, placed_count) in env.place_record[player.pos][food_type]
            if placed_player != p
                env.exchanges[food_type] += placed_count
            end
        end
        fframe[player.pos...] = 0
        env.place_record[player.pos][food_type] = small_dd_builder()    elseif place && player.food_counts[food_type] >= place_amount
        player.food_counts = Tuple(player.food_counts .-  (place_amount .* onehot(food_type, env.food_types)))
        fframe[player.pos...] += place_amount
        env.place_record[player.pos][food_type][p] = env.place_record[player.pos][food_type][p] + place_amount
    end
end

function comm(env::NatureEnv, p::Int, symbol::Int)
    push!(env.comms[env.step], (env.players[p].pos..., p, symbol))
end


text:
########################## Predicates #########################################

is_pos_int(x) = is_proper(x) && isinteger(x) && x > 0
is_neg_int(x) = is_proper(x) && isinteger(x) && x < 0
is_int_gt_one(x) = is_proper(x) && isinteger(x) && x > 1
is_pos(x) = is_proper(x) && x > 0
is_neg(x) = is_proper(x) && x < 0
is_neg_one(x) = is_proper(x) && (x ≈ -1)
is_pos_half(x) = is_proper(x) && (x ≈ 0.5)
is_neg_half(x) = is_proper(x) && (x ≈ -0.5)
is_abs_half(x) = is_proper(x) && (x ≈ 0.5 || x ≈ -0.5)

###############################################################################

terms(eq::SymbolicUtils.Add) = [t for t in arguments(eq)]
terms(eq) = [eq]

factors(eq, x) = isdependent(eq, x) ? [one(x), eq] : [one(x)]

function factors(eq::SymbolicUtils.Pow, x)
    p, k = arguments(eq)
    [p^(i*sign(k)) for i=0:abs(k)]
end

function factors(eq::SymbolicUtils.Mul, x)
    terms = [factors(q,x) for q in arguments(eq)]
    n = length(terms)

    l = Any[one(x)]

    for j = 1:n
        m = length(l)
        for t in terms[j]
            for k = 1:m
                push!(l, l[k]*t)
            end
        end
    end

    unique(l)
end

extract_power(eq::SymbolicUtils.Pow) = [arguments(eq)[2]]
extract_power(eq::SymbolicUtils.Term) = [1]
extract_power(eq::SymbolicUtils.Mul) = union([extract_power(t) for t in arguments(eq)]...)
extract_power(eq) = []

###############################################################################

is_var(eq) = isequal(eq, var(eq))

@syms Ω(x)

# rules to calculate the degree of a polynomial expression or NaN if not a poly
p_rules = [
    @rule Ω(+(~~xs)) => max(map(Ω, ~~xs)...)
    @rule Ω(~x * ~y) => Ω(~x) + Ω(~y)
    @rule Ω(^(~x, ~k::is_proper)) => is_pos_int(~k) ? ~k * Ω(~x) : NaN
    @rule Ω(~x::is_number) => 0
    @rule Ω(~x::is_var) => 1
    @rule Ω(~x) => NaN
]

function poly_deg(eq)
    ω = Prewalk(Chain(p_rules))(Ω(value(eq)))
    substitute(ω, Dict())
end

is_poly(eq) = !isnan(poly_deg(eq))
is_linear_poly(eq) = poly_deg(eq) == 1

# rules to extract the kernel (the solvable portion) of an expression
s_rules = [
    @rule Ω(+(~~xs)) => sum(map(Ω, ~~xs))
    @rule Ω(*(~~xs)) => prod(map(Ω, ~~xs))
    @rule Ω(~x / ~y) => Ω(~x) * Ω(^(~y,-1))
    @rule Ω(^(~x::is_linear_poly, ~k::is_pos_int)) => ^(~x, ~k)
    @rule Ω(~x::is_var) => ~x
    @rule Ω(~x::is_linear_poly) => ~x
    @rule Ω(sin(~x::is_linear_poly)) => sin(~x)
    @rule Ω(cos(~x::is_linear_poly)) => cos(~x)
    @rule Ω(sinh(~x::is_linear_poly)) => sinh(~x)
    @rule Ω(cosh(~x::is_linear_poly)) => cosh(~x)
    @rule Ω(exp(~x::is_linear_poly)) => exp(~x)
    @rule Ω(^(sin(~x::is_linear_poly), ~k::is_pos_int)) => ^(sin(~x), ~k)
    @rule Ω(^(cos(~x::is_linear_poly), ~k::is_pos_int)) => ^(cos(~x), ~k)
    @rule Ω(^(sinh(~x::is_linear_poly), ~k::is_pos_int)) => ^(sinh(~x), ~k)
    @rule Ω(^(cosh(~x::is_linear_poly), ~k::is_pos_int)) => ^(cosh(~x), ~k)
    @rule Ω(^(exp(~x::is_linear_poly), ~k::is_pos_int)) => ^(exp(~x), ~k)
    @rule Ω(^(~x, ~k)) => 1
    @rule Ω((~f)(~x)) => 1
    @rule Ω(~x) => 1
]

kernel(eq) = Prewalk(Chain(s_rules))(Ω(value(eq)))

coef(eq::SymbolicUtils.Mul, x) = prod(t for t in arguments(eq) if !isdependent(t,x); init=1)
coef(eq::SymbolicUtils.Add, x) = minimum(abs(coef(t,x)) for t in arguments(eq))
coef(eq, x) = is_number(eq) ? eq : 1

equivalent(eq, x) = eq / coef(eq, x)

########################## Transformation Rules ###############################

trigs_rules = [
    @rule tan(~x) => sin(~x) * cos(~x)^-1
    @rule sec(~x) => one(~x) * cos(~x)^-1
    @rule csc(~x) => one(~x) * sin(~x)^-1
    @rule cot(~x) => cos(~x) * sin(~x)^-1
    @rule sin(~n::is_int_gt_one*~x) => sin((~n - 1)*~x)*cos(~x) + cos((~n - 1)*~x)*sin(~x)
    @rule cos(~n::is_int_gt_one*~x) => cos((~n - 1)*~x)*cos(~x) - sin((~n - 1)*~x)*sin(~x)
    @rule tan(~n::is_int_gt_one*~x) => (tan((~n - 1)*~x) + tan(~x)) * (1 - tan((~n - 1)*~x)*tan(~x))^-1
    @rule csc(~n::is_int_gt_one*~x) => sec((~n - 1)*~x)*sec(~x)*csc((~n - 1)*~x)*csc(~x) * (sec((~n - 1)*~x)*csc(~x) + csc((~n - 1)*~x)*sec(~x))^-1
    @rule sec(~n::is_int_gt_one*~x) => sec((~n - 1)*~x)*sec(~x)*csc((~n - 1)*~x)*csc(~x) * (csc((~n - 1)*~x)*csc(~x) - sec((~n - 1)*~x)*sec(~x))^-1
    @rule cot(~n::is_int_gt_one*~x) => (cot((~n - 1)*~x)*cot(~x) - 1) * (cot((~n - 1)*~x) + cot(~x))^-1

    @rule ^(sin(~x), ~k::is_neg) => ^(csc(~x), -~k)
    @rule ^(cos(~x), ~k::is_neg) => ^(sec(~x), -~k)
    @rule ^(tan(~x), ~k::is_neg) => ^(cot(~x), -~k)
    @rule ^(csc(~x), ~k::is_neg) => ^(sin(~x), -~k)
    @rule ^(sec(~x), ~k::is_neg) => ^(cos(~x), -~k)
    @rule ^(cot(~x), ~k::is_neg) => ^(tan(~x), -~k)

    @rule sin(~x + ~y) => sin(~x)*cos(~y) + cos(~x)*sin(~y)
    @rule cos(~x + ~y) => cos(~x)*cos(~y) - sin(~x)*sin(~y)
    @rule tan(~x + ~y) => (tan(~x) + tan(~y)) * (1 - tan(~x)*tan(~y))^-1
    @rule csc(~x + ~y) => sec(~x)*sec(~y)*csc(~x)*csc(~y) * (sec(~x)*csc(~y) + csc(~x)*sec(~y))^-1
    @rule sec(~x + ~y) => sec(~x)*sec(~y)*csc(~x)*csc(~y) * (csc(~x)*csc(~y) - sec(~x)*sec(~y))^-1
    @rule cot(~x + ~y) => (cot(~x)*cot(~y) - 1) * (cot(~x) + cot(~y))^-1

    @rule sin(~x - ~y) => sin(~x)*cos(~y) - cos(~x)*sin(~y)
    @rule cos(~x - ~y) => cos(~x)*cos(~y) + sin(~x)*sin(~y)
    @rule tan(~x - ~y) => (tan(~x) - tan(~y)) * (1 + tan(~x)*tan(~y))^-1
    @rule csc(~x - ~y) => sec(~x)*sec(~y)*csc(~x)*csc(~y) * (sec(~x)*csc(~y) - csc(~x)*sec(~y))^-1
    @rule sec(~x - ~y) => sec(~x)*sec(~y)*csc(~x)*csc(~y) * (csc(~x)*csc(~y) + sec(~x)*sec(~y))^-1
    @rule cot(~x - ~y) => (cot(~x)*cot(~y) + 1) * (cot(~x) - cot(~y))^-1

    # @rule sin(2*~x) => 2*sin(~x)*cos(~x)
    # @rule cos(2*~x) => 2*cos(~x)^2 - 1
    # @rule tan(2*~x) => 2*tan(~x) / (1 - tan(~x)^2)
    # @rule cot(2*~x) => (cot(~x)^2 - 1) / (2*cot(~x))
    # @rule sec(2*~x) => sec(~x)^2 / (2 - sec(~x)^2)
    # @rule csc(2*~x) => sec(~x)*csc(~x) / 2

    # @rule sin(3*~x) => 3*sin(~x) - 4*sin(~x)^3
    # @rule cos(3*~x) => 4*cos(~x)^3 - 3*cos(~x)
    # @rule tan(3*~x) => (3*tan(~x) - tan(~x)^3) / (1 - 3*tan(~x)^2)
    # @rule cot(3*~x) => (3*cot(~x) - cot(~x)^3) / (1 - 3*cot(~x)^2)
    # @rule sec(3*~x) => sec(~x)^3 / (4 - 3*sec(~x)^2)
    # @rule csc(3*~x) => csc(~x)^3 / (3*csc(~x)^2 - 4)

    @rule tanh(~x) => sinh(~x) * cosh(~x)^-1
    @rule sech(~x) => one(~x) * cosh(~x)^-1
    @rule csch(~x) => one(~x) * sinh(~x)^-1
    @rule coth(~x) => cosh(~x) * sinh(~x)^-1

    @rule sqrt(~x) => ^(~x, 0.5)
    @acrule exp(~x) * exp(~y) => exp(~x + ~y)
]

simplify_trigs(eq) = expand(Fixpoint(Prewalk(PassThrough(Chain(trigs_rules))))(value(eq)))

##############################################################################

function factor_poly(p)
    x = var(p)
    r, s = find_roots(p, x)
    s = s[1:2:end]
    r = nice_parameter.(r)
    s = nice_parameter.(s)
    return [[(x - u) for u in r];
            [(x^2 - 2*real(u)*x + abs2(u)) for u in s]
           ]
end

function decompose_rational(eq)
    if poly_deg(eq) == 1 return inverse(eq) end
    x = var(eq)

    d = poly_deg(eq)

    for j = 1:10  # will try 10 times to find the roots
        r, s = find_roots(eq, x)
        if length(r) + length(s) >= d break end
    end

    s = s[1:2:end]
    r = nice_parameter.(r)
    s = nice_parameter.(s)

    F = Any[(x^2 - 2*real(u)*x + abs2(u))^-1 for u in s] ∪
        [x*(x^2 - 2*real(u)*x + abs2(u))^-1 for u in s]
    for i in eachindex(r)
        μ = sum(r[1:i] .== r[i])
        push!(F, (x-r[i])^-μ)
    end
    F = unique(F)

    n = length(F)
    A = zeros(Complex, (n,n))
    b = zeros(Complex, n)

    for i = 1:n
        x₀ = test_point(false, 1.0)
        d = Dict(x => x₀)
        b[i] = 1 / substitute(eq, d)
        for j = 1:n
            A[i,j] = substitute(F[j], d)
        end
    end

    q₀ = A \ b
    q = nice_parameter.(q₀)
    p = sum(F[i]*q[i] for i=1:n if q[i] != 0; init=zero(x))
    return p
end

rational_rules = [
    @rule Ω(+(~~xs)) => sum(map(Ω, ~~xs))
    @rule Ω(*(~~xs)) => prod(map(Ω, ~~xs))
    @rule Ω(^(~x::is_poly, ~k::is_neg_int)) => expand(^(decompose_rational(~x), -~k))
    @rule Ω(~x / ^(~y::is_poly,~k::is_pos_int)) => expand(~x * ^(decompose_rational(~y), ~k))
    @rule Ω(~x / ~y::is_poly) => expand(~x * decompose_rational(~y))
    # @rule Ω(^(~x,~k)) => ^(~x, ~k)
    @rule Ω(sqrt(~x::is_poly)) => prod(sqrt(f) for f in factor_poly(~x))
    @rule Ω(log(~x::is_poly)) => sum(log(f) for f in factor_poly(~x))
    @rule Ω(~x) => ~x
]

factor_rational(eq) = Prewalk(PassThrough(Chain(rational_rules)))(Ω(value(eq)))

###############################################################################

div_rule = @rule ~x / ~y => ~x * ^(~y, -1)

apply_div_rule(eq) = Prewalk(PassThrough(div_rule))(value(eq))

###############################################################################

inv_rules = [
    @rule Ω(1 / ~x) => ~x
    @rule Ω(~x / ~y) => Ω(~x) * ~y
    @rule Ω(^(~x, -1)) => ~x
    @rule Ω(^(~x, ~k)) => ^(Ω(~x), ~k)
    @rule Ω(~x) => ^(~x, -1)
]

inverse(eq) = Prewalk(PassThrough(Chain(inv_rules)))(Ω(value(eq)))

###############################################################################

# the integration rules for non-solvable portion of an expression

@syms 𝛷(x)

d_rules = [
    @rule 𝛷(sin(~x)) => one(~x) + sin(~x) + cos(~x)
    @rule 𝛷(cos(~x)) => one(~x) + cos(~x) + sin(~x)
    @rule 𝛷(tan(~x)) => one(~x) + tan(~x) + log(cos(~x))
    @rule 𝛷(csc(~x)) => one(~x) + csc(~x) + log(sin(~x)^-1 - cos(~x)*sin(~x)^-1)
    @rule 𝛷(sec(~x)) => one(~x) + sec(~x) + log(cos(~x)^-1 + sin(~x)*cos(~x)^-1)
    @rule 𝛷(cot(~x)) => one(~x) + cot(~x) + log(sin(~x))

    @rule 𝛷(sinh(~x)) => one(~x) + sinh(~x) + cosh(~x)
    @rule 𝛷(cosh(~x)) => one(~x) + cosh(~x) + sinh(~x)
    @rule 𝛷(tanh(~x)) => one(~x) + tanh(~x) + log(cosh(~x))
    @rule 𝛷(csch(~x)) => one(~x) + csch(~x) + log(sinh(~x)^-1 - cosh(~x)*sinh(~x)^-1)
    @rule 𝛷(sech(~x)) => one(~x) + sech(~x) + log(cosh(~x)^-1 + sinh(~x)*cosh(~x)^-1)
    @rule 𝛷(coth(~x)) => one(~x) + coth(~x) + log(sinh(~x))

    @rule 𝛷(asin(~x)) => one(~x) + asin(~x) + ~x*asin(~x) + sqrt(1 - ~x*~x)
    @rule 𝛷(acos(~x)) => one(~x) + acos(~x) + ~x*acos(~x) + sqrt(1 - ~x*~x)
    @rule 𝛷(atan(~x)) => one(~x) + atan(~x) + ~x*atan(~x) + log(~x*~x + 1)
    @rule 𝛷(acsc(~x)) => one(~x) + acsc(~x)
    @rule 𝛷(asec(~x)) => one(~x) + asec(~x)
    @rule 𝛷(acot(~x)) => one(~x) + acot(~x) + ~x*acot(~x) + log(~x*~x + 1)

    @rule 𝛷(asinh(~x)) => one(~x) + asinh(~x) + ~x*asinh(~x) + sqrt(~x*~x + 1)
    @rule 𝛷(acosh(~x)) => one(~x) + acosh(~x) + ~x*acosh(~x) + sqrt(~x*~x - 1)
    @rule 𝛷(atanh(~x)) => one(~x) + atanh(~x) + ~x*atanh(~x) + log(~x + 1)
    @rule 𝛷(acsch(~x)) => one(~x) + acsch(~x)
    @rule 𝛷(asech(~x)) => one(~x) + asech(~x)
    @rule 𝛷(acoth(~x)) => one(~x) + acoth(~x) + ~x*acot(~x) + log(~x + 1)

    @rule 𝛷(log(~x)) => one(~x) + log(~x) + ~x + ~x * log(~x) + 𝛷(inverse(~x))
    @rule 𝛷(sqrt(~x)) => sum(candidate_sqrt(~x,0.5); init=one(~x))
    @rule 𝛷(^(sqrt(~x),-1)) => 𝛷(^(~x,-0.5))
    @rule 𝛷(cbrt(~x)) => one(~x) + cbrt(~x) + ^(~x, 4/3)

    @rule 𝛷(^((~f)(~x), ~k::is_int_gt_one)) => var(~x) * ^((~f)(~x), ~k) + 𝛷(var(~x)*^((~f)(~x), ~k-1))
    @rule 𝛷(^(~x, ~k::is_abs_half)) => sum(candidate_sqrt(~x,~k); init=𝛷(~x))
    @rule 𝛷(^(~x, ~k::is_neg)) => sum(candidate_pow_minus(~x,~k); init=𝛷(~x))
    @rule 𝛷(^(~x::is_linear_poly, ~k::is_pos)) => one(~x) + ^(~x,~k) + ^(~x,~k+1)
    @rule 𝛷(^(~x, ~k::is_pos)) => one(~x) + ^(~x,~k) + ^(~x,~k+1)

    @rule 𝛷(exp(~x)) => one(~x) + exp(~x)
    @rule 𝛷(+(~~xs)) => sum(map(𝛷, ~~xs))
    @rule 𝛷(*(~~xs)) => prod(map(𝛷, ~~xs))
    @rule 𝛷(~x / ~y) => 𝛷(~x * inverse(~y))
    @rule 𝛷(~x) => one(~x) + ~x
    @rule 𝛷(-~x) => -𝛷(~x)
]

apply_d_rules(eq) = expand(Fixpoint(Prewalk(PassThrough(Chain(d_rules))))(𝛷(value(eq))))


text:
############################################################################################
# MALA specific diagnostics
"""
$(TYPEDEF)
Diagnostics for MALA sampler.

# Fields
$(TYPEDFIELDS)
"""
struct DiagnosticsMALA{T<:AbstractFloat} <: MCMCKernelDiagnostics
    "Discretization size"
    ϵ::T
    function DiagnosticsMALA(ϵ::T) where {T<:AbstractFloat}
        @argcheck ϵ > 0.0 "Discretization size has to be positive"
        return new{T}(ϵ)
    end
end

function results(
    diagnostics::AbstractVector{DiagnosticsMALA{F}}, Ndigits::Integer, quantiles::Vector{T}
) where {F<:AbstractFloat,T<:Real}
    println(
        "MALA sampler has average (final) stepsize of ",
        round.(
            mean(diagnostics[iter].ϵ for iter in eachindex(diagnostics)); digits=Ndigits
        ),
        " (",
        round.(diagnostics[end].ϵ; digits=Ndigits),
        ").",
    )
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Show relevant diagnostic results.

# Examples
```julia
```

"""
function generate_showvalues(diagnostics::D) where {D<:DiagnosticsMALA}
    return function showvalues()
        return (:MALA, "diagnostics"), (:ϵ, diagnostics.ϵ)
    end
end

############################################################################################
#export
export DiagnosticsMALA, generate_showvalues


text:
using Distributions
using Test

@testset "Dirac tests" begin
    for val in (3, 3.0, -3.5)
        d = Dirac(val)

        @test minimum(d) == val
        @test maximum(d) == val
        @test !insupport(d, prevfloat(float(val)))
        @test insupport(d, val)
        @test !insupport(d, nextfloat(float(val)))
        @test support(d) == (val,)

        @test iszero(pdf(d, prevfloat(float(val))))
        @test isone(pdf(d, val))
        @test iszero(pdf(d, nextfloat(float(val))))

        @test logpdf(d, prevfloat(float(val))) == -Inf
        @test iszero(logpdf(d, val))
        @test logpdf(d, nextfloat(float(val))) == -Inf

        @test iszero(cdf(d, prevfloat(float(val))))
        @test isone(cdf(d, val))
        @test isone(cdf(d, nextfloat(float(val))))

        @test quantile(d, 0) == val
        @test quantile(d, 0.5) == val
        @test quantile(d, 1) == val

        @test rand(d) == val

        @test mean(d) == val
        @test iszero(var(d))

        @test mode(d) == val

        @test iszero(entropy(d))

        t = rand()
        @test mgf(d, t) == exp(t * val)
        @test cf(d, t) == cis(t * val)
    end
end


text:
# Based on FreeBSD lib/msun/src/e_exp.c
# which is made available under the following licence

## Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. Permission
## to use, copy, modify, and distribute this software is freely granted,
## provided that this notice is preserved.

# Method
# 1. Argument reduction: Reduce x to an r so that |r| <= 0.5*ln(2). Given x,
#    find r and integer k such that
#       x = k*ln(2) + r,  |r| <= 0.5*ln(2).
#    Here r is represented as r = hi - lo for better accuracy.
#
# 2. Approximate exp(r) by a special rational function on [0, 0.5*ln(2)]:
#       R(r^2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r^4/360 + ...
#
#    A special Remez algorithm on [0, 0.5*ln(2)] is used to generate a
#    polynomial to approximate R.
#
#    The computation of exp(r) thus becomes
#                       2*r
#       exp(r) = 1 + ----------
#                     R(r) - r
#                          r*c(r)
#              = 1 + r + ----------- (for better accuracy)
#                         2 - c(r)
#    where
#       c(r) = r - (P1*r^2  + P2*r^4  + ... + P5*r^10 + ...).
#
# 3. Scale back: exp(x) = 2^k * exp(r)

# log(2)
const LN2 = 6.931471805599453094172321214581765680755001343602552541206800094933936219696955e-01
# log2(e)
const LOG2E = 1.442695040888963407359924681001892137426646

# log(2) into upper and lower
LN2U(::Type{Float64}) = 6.93147180369123816490e-1
LN2U(::Type{Float32}) = 6.9313812256f-1

LN2L(::Type{Float64}) = 1.90821492927058770002e-10
LN2L(::Type{Float32}) = 9.0580006145f-6

# max and min arguments for exponential fucntions
MAXEXP(::Type{Float64}) = 7.09782712893383996732e2 # log 2^1023*(2-2^-52)
MAXEXP(::Type{Float32}) = 88.72283905206835f0      # log 2^127 *(2-2^-23)

# one less than the min exponent since we can sqeeze a bit more from the exp function
MINEXP(::Type{Float64}) = -7.451332191019412076235e2 # log 2^-1075
MINEXP(::Type{Float32}) = -103.97207708f0            # log 2^-150

@inline exp_kernel(x::Float64) = @horner(x, 1.66666666666666019037e-1,
        -2.77777777770155933842e-3, 6.61375632143793436117e-5,
        -1.65339022054652515390e-6, 4.13813679705723846039e-8)

@inline exp_kernel(x::Float32) = @horner(x, 1.6666625440f-1, -2.7667332906f-3)

# for values smaller than this threshold just use a Taylor expansion
@eval exp_small_thres(::Type{Float64}) = $(2.0^-28)
@eval exp_small_thres(::Type{Float32}) = $(2.0f0^-13)

"""
    exp(x)

Compute the natural base exponential of `x`, in other words ``e^x``.
"""
function exp{T<:Union{Float32,Float64}}(x::T)
    xa = reinterpret(Unsigned, x) & ~sign_mask(T)
    xsb = signbit(x)

    # filter out non-finite arguments
    if xa > reinterpret(Unsigned, MAXEXP(T))
        if xa >= exponent_mask(T)
            xa & significand_mask(T) != 0 && return T(NaN)
            return xsb ? T(0.0) : T(Inf) # exp(+-Inf)
        end
        x > MAXEXP(T) && return T(Inf)
        x < MINEXP(T) && return T(0.0)
    end

    # This implementation gives 2.7182818284590455 for exp(1.0) when T ==
    # Float64, which is well within the allowable error; however,
    # 2.718281828459045 is closer to the true value so we prefer that answer,
    # given that 1.0 is such an important argument value.
    if x == T(1.0) && T == Float64
        return 2.718281828459045235360
    end

    if xa > reinterpret(Unsigned, T(0.5)*T(LN2)) # |x| > 0.5 log(2)
        # argument reduction
        if xa < reinterpret(Unsigned, T(1.5)*T(LN2)) # |x| < 1.5 log(2)
            if xsb
                k = -1
                hi = x + LN2U(T)
                lo = -LN2L(T)
            else
                k = 1
                hi = x - LN2U(T)
                lo = LN2L(T)
            end
        else
            n = round(T(LOG2E)*x)
            k = unsafe_trunc(Int,n)
            hi = muladd(n, -LN2U(T), x)
            lo = n*LN2L(T)
        end
        r = hi - lo

        # compute approximation on reduced argument
        z = r*r
        p = r - z*exp_kernel(z)
        y = T(1.0) - ((lo - (r*p)/(T(2.0) - p)) - hi)

        # scale back
        if k > -significand_bits(T)
            # multiply by 2.0 first to prevent overflow, which helps extends the range
            k == exponent_max(T) && return y * T(2.0) * T(2.0)^(exponent_max(T) - 1)
            twopk = reinterpret(T, rem(exponent_bias(T) + k, fpinttype(T)) << significand_bits(T))
            return y*twopk
        else
            # add significand_bits(T) + 1 to lift the range outside the subnormals
            twopk = reinterpret(T, rem(exponent_bias(T) + significand_bits(T) + 1 + k, fpinttype(T)) << significand_bits(T))
            return y * twopk * T(2.0)^(-significand_bits(T) - 1)
        end
    elseif xa < reinterpret(Unsigned, exp_small_thres(T)) # |x| < exp_small_thres
        # Taylor approximation for small x
        return T(1.0) + x
    else
        # primary range with k = 0, so compute approximation directly
        z = x*x
        p = x - z*exp_kernel(z)
        return T(1.0) - ((x*p)/(p - T(2.0)) - x)
    end
end


text:
# This file contains custom exceptions
# for functions related to input data

@exception(
    "Exception thrown when the file is not found.",
    "Input",
    "ScatsInputNotAFile",
    "The file is not found (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when the path is a directory.",
    "Input",
    "ScatsInputIsADir",
    "Specified path is a directory (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an end of file occurred.",
    "Input",
    "ScatsInputEOF",
    "Unexpected end of file (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `N` is wrong.",
    "Input",
    "ScatsInputWR_N",
    "Wrong input: N (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `Δt` is wrong.",
    "Input",
    "ScatsInputWR_Δt",
    "Wrong input: Δt (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `q` is wrong.",
    "Input",
    "ScatsInputWR_q",
    "Wrong input: q (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `t` is wrong.",
    "Input",
    "ScatsInputWR_t",
    "Wrong input: t (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `x` is wrong.",
    "Input",
    "ScatsInputWR_x",
    "Wrong input: x (\\\"\", e.file, \"\\\").",
)


text:
using Images
using TestImages
using StaticArrays
using OffsetArrays
using Random
using CoordinateTransformations

function evaluate_error(dims, flow::Array{SVector{2, Float64}, 2},  Δ,  tol)
    error_count = 0
    maximum_error = 0
    for r = 1:first(dims)
        for c = 1:last(dims)
            δ = flow[r,c]
            if abs(δ[1] - Δ[1]) > tol || abs(δ[2] - Δ[2]) > tol
                    error_count += 1
            end
            error = sqrt( (δ[1] - Δ[1])^2 + (δ[2] - Δ[2])^2 )
            if error > maximum_error
                maximum_error = error
            end
        end
    end
    return error_count, maximum_error
end

@testset "Farneback" begin
    @testset "Polynomial Expansion" begin
        @info "Running Polynomial Expansion test."
        img = Gray{Float64}.(testimage("mandrill"))
        @info "Running Polynomial Expansion Convolution Implementation."
        @time A, B, C = polynomial_expansion(ConvolutionImplementation(), img, 6, 2)
        @info "Running Polynomial Expansion Matrix Implementation."
        @time P, Q, R = polynomial_expansion(MatrixImplementation(), img, 6, 2)
        @test sum(abs.(vec(A)-vec(P))) ≈  0 atol=1e-9
        @test sum(abs.(vec(B)-vec(Q))) ≈  0 atol=1e-9
        @test sum(abs.(vec(C)-vec(R))) ≈  0 atol=1e-9

    end

    @testset "Horizontal Motion" begin
        @info "Running Horizontal Motion test."
        maximum_percentage_error = 7.5
        img1 = Gray{Float64}.(testimage("mandrill"))

        tol = 0.3
        Δ = (0.0, 3.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        algorithm = Farneback(7, estimation_window = 19,
                                 σ_estimation_window = 6.0,
                                 expansion_window = 6,
                                 σ_expansion_window = 1.5)

        @time flow = optical_flow(img1, img2, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Horizontal Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4

        # Same as above, except that we pass in an initial displacement field.
        displacement = Array{SVector{2, Float64}, 2}(undef,size(img1))
        for i in eachindex(displacement)
                displacement[i] = SVector{2, Float64}(0.0, 0.0)
        end
        @time flow = optical_flow(img1, img2, displacement, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Horizontal Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4
    end

    @testset "Vertical Motion" begin
        @info "Running Vertical Motion test."
        maximum_percentage_error = 7.5
        img1 = Gray{Float64}.(testimage("mandrill"))

        tol = 0.3
        Δ = (3.0, 0.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        algorithm = Farneback(7, estimation_window = 19,
                                 σ_estimation_window = 6.0,
                                 expansion_window = 6,
                                 σ_expansion_window = 1.5)

        @time flow = optical_flow(img1, img2, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Vertical Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4

        # Same as above, except that we pass in an initial displacement field.
        displacement = Array{SVector{2, Float64}, 2}(undef,size(img1))
        for i in eachindex(displacement)
                displacement[i] = SVector{2, Float64}(0.0, 0.0)
        end
        @time flow = optical_flow(img1, img2, displacement, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Vertical Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4
    end

    @testset "Combined Motion" begin
        @info "Running Combined  Motion test."
        maximum_percentage_error = 13
        img1 = Gray{Float64}.(testimage("mandrill"))

        tol = 0.3
        Δ = (3.0, -1.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        algorithm = Farneback(7, estimation_window = 19,
                                 σ_estimation_window = 6.0,
                                 expansion_window = 6,
                                 σ_expansion_window = 1.5)

        @time flow = optical_flow(img1, img2, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Combined Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4

        # Same as above, except that we pass in an initial displacement field.
        displacement = Array{SVector{2, Float64}, 2}(undef,size(img1))
        for i in eachindex(displacement)
                displacement[i] = SVector{2, Float64}(0.0, 0.0)
        end
        @time flow = optical_flow(img1, img2, displacement, algorithm)

        error_count, maximum_error = evaluate_error(size(img1), flow, Δ, tol)
        percentage_error = (error_count / prod(size(img1))) * 100

        @info "Case: Combined Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
        @test maximum_error < 4

    end

end


text:
"""
    y = bquantize(x, nsd=3, abstol=eps, reltol=10*eps)

Bidirectionally quantize a n by 1 vector x to nsd signed digits, terminate early
if the error is less than the specified tolerances.

`y` is a 2-tuple of arrays with the following fields:
- `y`[1][i] is the quantized value in floating-point form
- `y`[2][i] is a 2-by-nsd (or less) matrix containing the powers of two
    (first row) and their signs (second row).

See also bunquantize.
"""
function bquantize(x, nsd=3, abstol=eps(Float64), reltol=10*eps(Float64))
    n = length(x)
    offset = -log2(0.75)

    y_val = zeros(n)
    y_csd = Array{Matrix{Int}}(undef, n)

    for i = 1:n
        xp = x[i]
        y_val[i] = 0.0
        for j = 1:nsd
            error = abs(y_val[i] - x[i])
            if error <= abstol || error <= abs(x[i])*reltol
                break
            end
            p = floor(Int, log2(abs(xp)) + offset)
            p2 = 2.0 .^ p
            sx = cmp(xp, 0)
            xp = xp - sx*p2
            y_val[i] = y_val[i] + sx*p2
            if isassigned(y_csd, i)
                y_csd[i] = hcat(y_csd[i], [p; sx])
            else
                y_csd[i] = reshape([p; sx], (2,1))
            end
        end
    end

    return y_val, y_csd
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
module MNIST

using GZip, Colors

const Gray = Colors.Gray{Colors.N0f8}

const dir = joinpath(@__DIR__, "../../deps/mnist")

function load()
  mkpath(dir)
  cd(dir) do
    for file in ["train-images-idx3-ubyte",
                 "train-labels-idx1-ubyte",
                 "t10k-images-idx3-ubyte",
                 "t10k-labels-idx1-ubyte"]
      isfile(file) && continue
      info("Downloading MNIST dataset")
      download("https://cache.julialang.org/http://yann.lecun.com/exdb/mnist/$file.gz", "$file.gz")
      open(file, "w") do io
        write(io, GZip.open(read, "$file.gz"))
      end
    end
  end
end

const IMAGEOFFSET = 16
const LABELOFFSET = 8

const NROWS = 28
const NCOLS = 28

const TRAINIMAGES = joinpath(dir, "train-images-idx3-ubyte")
const TRAINLABELS = joinpath(dir, "train-labels-idx1-ubyte")
const TESTIMAGES = joinpath(dir, "t10k-images-idx3-ubyte")
const TESTLABELS = joinpath(dir, "t10k-labels-idx1-ubyte")

function imageheader(io::IO)
  magic_number = bswap(read(io, UInt32))
  total_items = bswap(read(io, UInt32))
  nrows = bswap(read(io, UInt32))
  ncols = bswap(read(io, UInt32))
  return magic_number, Int(total_items), Int(nrows), Int(ncols)
end

function labelheader(io::IO)
  magic_number = bswap(read(io, UInt32))
  total_items = bswap(read(io, UInt32))
  return magic_number, Int(total_items)
end

function rawimage(io::IO)
  img = Array{Gray}(NCOLS, NROWS)
  for i in 1:NCOLS, j in 1:NROWS
    img[i, j] = reinterpret(Colors.N0f8, read(io, UInt8))
  end
  return img
end

function rawimage(io::IO, index::Integer)
  seek(io, IMAGEOFFSET + NROWS * NCOLS * (index - 1))
  return rawimage(io)
end

rawlabel(io::IO) = Int(read(io, UInt8))

function rawlabel(io::IO, index::Integer)
  seek(io, LABELOFFSET + (index - 1))
  return rawlabel(io)
end

getfeatures(io::IO, index::Integer) = vec(getimage(io, index))

"""
    images()
    images(:test)

Load the MNIST images.

Each image is a 28×28 array of `Gray` colour values (see Colors.jl).

Returns the 60,000 training images by default; pass `:test` to retreive the
10,000 test images.
"""
function images(set = :train)
  load()
  io = IOBuffer(read(set == :train ? TRAINIMAGES : TESTIMAGES))
  _, N, nrows, ncols = imageheader(io)
  [rawimage(io) for _ in 1:N]
end

"""
    labels()
    labels(:test)

Load the labels corresponding to each of the images returned from `images()`.
Each label is a number from 0-9.

Returns the 60,000 training labels by default; pass `:test` to retreive the
10,000 test labels.
"""
function labels(set = :train)
  load()
  io = IOBuffer(read(set == :train ? TRAINLABELS : TESTLABELS))
  _, N = labelheader(io)
  [rawlabel(io) for _ = 1:N]
end

end # module


text:
# This example is a WIP. It works, but is a bit rough around the edges.

using Stheno, Random, Distributions, Plots, ColorTypes, FixedPointNumbers
using Stheno: @model, ∇

rng = MersenneTwister(12345);

@model function gp()
    f = GP(EQ())
    ∇f = ∇(f)
    return f, ∇f
end

f, ∇f = gp()

# Sample prior at random locations.
Xf = rand(Uniform(-5, 5), 5)
X∇f = rand(Uniform(-5, 5), 15)
X∇f = Xf
f_, ∇f_ = rand([f(Xf), ∇f(X∇f)])

# Compute log prob. of samples.
logpdf([f(Xf), ∇f(X∇f)], [f_, ∇f_])

# Compute posterior processes.
f′, ∇f′ = (f, ∇f) | (f(Xf) ← f_, ∇f(X∇f) ← ∇f_)

#######################  Do posterior inference give a few samples  #######################

Np, S = 500, 10;
Xp = range(-6; stop=6, length=Np);
f′Xp = rand(rng, f′(Xp), S);

∇Xp = range(-6; stop=6, length=Np);
∇f′Xp = rand(rng, ∇f′(∇Xp), S);

# Get posterior mean and marginals f′ and y′ and write them for plotting.
μ′f, σ′f = marginals(f′(Xp));
μ′∇f, σ′∇f = marginals(∇f′(∇Xp));

####################################  Plot results  ####################################

plotly();

posterior_plot = plot(
    legend=:topright,
    legendfont=Plots.Font(
        "sans-serif",
        10,
        :hcenter,
        :vcenter,
        0.0,
        RGB{Normed{UInt8, 8}}(0.0,0.0,0.0)
    ),
    background_color_legend=RGBA(1, 1, 1, 0),
    foreground_color_legend=RGBA(1, 1, 1, 0),
);

plot!(posterior_plot, Xp, μ′f;
    linewidth=2.0,
    linecolor=:blue,
    label="");
plot!(posterior_plot, Xp, μ′f;
    linewidth=0.0,
    linecolor=:blue,
    linealpha=0.2,
    fillrange=[μ′f.- 3  .* σ′f],
    fillalpha=0.2,
    fillcolor=:blue,
    label=["posterior, f", ""]);
plot!(posterior_plot, Xp, f′Xp[:, 1];
    linecolor=:blue,
    linealpha=0.3,
    label="posterior samples, f");
plot!(posterior_plot, Xp, μ′f;
    linewidth=0.0,
    linecolor=:blue,
    linealpha=0.2,
    fillrange=[μ′f .+ 3 .* σ′f],
    fillalpha=0.2,
    fillcolor=:blue,
    label="");
plot!(posterior_plot, Xp, f′Xp;
    linecolor=:blue,
    linealpha=0.3,
    label="");
scatter!(posterior_plot, Xf, f_;
    markercolor=:green,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=5,
    markeralpha=0.9,
    label="f_");

colour = :red
plot!(posterior_plot, ∇Xp, μ′∇f;
    linewidth=2.0,
    linecolor=colour,
    label="");
plot!(posterior_plot, ∇Xp, μ′∇f;
    linewidth=0.0,
    linecolor=colour,
    linealpha=0.2,
    fillrange=[μ′∇f.- 3  .* σ′∇f],
    fillalpha=0.2,
    fillcolor=colour,
    label=["posterior, ∇f", ""]);
plot!(posterior_plot, ∇Xp, ∇f′Xp[:, 1];
    linecolor=colour,
    linealpha=0.3,
    label="posterior samples, ∇f");
plot!(posterior_plot, ∇Xp, μ′∇f;
    linewidth=0.0,
    linecolor=colour,
    linealpha=0.2,
    fillrange=[μ′∇f .+ 3 .* σ′∇f],
    fillalpha=0.2,
    fillcolor=colour,
    label="");
plot!(posterior_plot, ∇Xp, ∇f′Xp;
    linecolor=colour,
    linealpha=0.3,
    label="");
scatter!(posterior_plot, X∇f, ∇f_;
    markercolor=colour,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=3,
    markeralpha=0.9,
    label="∇f_");

display(posterior_plot);


text:
"""
$(SIGNATURES)

Convenience function to process markdown definitions `@def ...` as appropriate.
Depending on `isconfig`, will update `GLOBAL_VARS` or `LOCAL_VARS`.

**Arguments**

* `blocks`:    vector of active docs
* `isconfig`:  whether the file being processed is the config file
                (--> global page variables)
"""
function process_mddefs(blocks::Vector{OCBlock}, isconfig::Bool,
                        pagevar::Bool=false)::Nothing

    (:process_mddefs, "config: $isconfig, pagevar: $pagevar") |> logger

    # Blocks of definitions just get evaluated in an anonymous Module
    curdict = ifelse(isconfig, GLOBAL_VARS, LOCAL_VARS)
    for mdb in filter(β -> (β.name == :MD_DEF_BLOCK), blocks)
        inner = stent(mdb)
        exs = parse_code(inner)
        mdl = newmodule("MD_DEFINITIONS")
        try
            foreach(ex -> Core.eval(mdl, ex), exs)
        catch
            error("Encountered an error on $(locvar(:fd_rpath)) while trying to " *
                  "evaluate a block of definitions (`+++...+++`).")
        end
        # get the variable names (from all assignments)
        vnames = [ex.args[1] for ex in exs if ex.head == :(=)]
        filter!(v -> v isa Symbol, vnames)
        for vname in vnames
            key    = String(vname)
            value  = getproperty(mdl, vname)
            set_var!(curdict, key, value; isglobal=isconfig)
        end
    end

    # Find all markdown definitions (MD_DEF) blocks
    mddefs = filter(β -> (β.name == :MD_DEF), blocks)
    # empty container for the assignments
    assignments = Vector{Pair{String, String}}()
    # go over the blocks, and extract the assignment
    for (i, mdd) ∈ enumerate(mddefs)
        inner = stent(mdd)
        m = match(ASSIGN_PAT, inner)
        if isnothing(m)
            print_warning("""
                Delimiters for an '@def ...' environment were found but at
                least one assignment does not have the proper syntax. That
                assignment will be ignored.
                \nRelevant pointers:
                $POINTER_PV
                """)
            continue
        end
        vname, vdef = m.captures[1:2]
        push!(assignments, (String(vname) => String(vdef)))
    end

    # if in config file, update `GLOBAL_VARS` and return
    rpath = splitext(locvar(:fd_rpath)::String)[1]
    if isconfig
        set_vars!(GLOBAL_VARS, assignments, isglobal=true)
        return nothing
    end

    # otherwise set local vars
    set_vars!(LOCAL_VARS, assignments)
    rpath = splitext(locvar(:fd_rpath)::String)[1]

    hasmath = locvar(:hasmath)::Bool
    hascode = locvar(:hascode)::Bool
    if !hascode && globvar("autocode")::Bool
        # check and set hascode automatically
        code = any(b -> startswith(string(b.name), "CODE_BLOCK"), blocks)
        set_var!(LOCAL_VARS, "hascode", code)
    end
    if !hasmath && globvar("automath")::Bool
        # check and set hasmath automatically
        math = any(b -> b.name in MATH_BLOCKS_NAMES, blocks)
        set_var!(LOCAL_VARS, "hasmath", math)
    end

    (:process_mddefs, "assignments done") |> logger

    # TAGS
    tags = Set(refstring.(locvar(:tags)::Vector{String}))
    # Cases:
    # 0. there was no page tags before
    #   a. tags is empty --> do nothing
    #   b. tags is not empty --> initialise global page tags + add rpath=>tags
    # 1. that page did not have tags
    #   a. tags is empty --> do nothing
    #   b. tags is not empty register them and update all
    # 2. that page did have tags
    #   a. tags are unchanged --> do nothing
    #   b. check which ones change and update those
    PAGE_TAGS = globvar("fd_page_tags")
    if isnothing(PAGE_TAGS)
        isempty(tags) && return nothing
        set_var!(GLOBAL_VARS, "fd_page_tags", DTAG((rpath => tags,)); check=false)
    elseif !haskey(PAGE_TAGS, rpath)
        isempty(tags) && return nothing
        PAGE_TAGS[rpath] = tags
    else
        old_tags = PAGE_TAGS[rpath]
        if isempty(tags)
            delete!(PAGE_TAGS, rpath)
        else
            PAGE_TAGS[rpath] = tags
        end
        # we will need to update all tag pages that were or are related to
        # the current rpath, indeed properties of the page may have changed
        # and affect these derived pages (that's why it's a union here
        # and not a setdiff).
        tags = union(old_tags, tags)
    end
    # In the full pass each page is processed first (without generating tag
    # pages) and then, when all tags have been gathered, generate_tag_pages
    # is called (see `fd_fullpass`).
    # During the serve loop, we want to trigger on page change.
    pagevar || FD_ENV[:FULL_PASS] || generate_tag_pages(tags)
    return nothing
end


text:
experimentdesign_insts(ext::MatSpecExt) = [
    [(d, ext)
    for d in vcat(3, 25:25:100, 150, 200:100:1100)] # includes compile run
    ]

insts = OrderedDict()
insts["nat"] = (nothing, vcat(experimentdesign_insts.([
    MatNegRtdet(),
    MatNegLog(),
    MatNegSqrtConj(),
    MatNegPower01(1/3),
    ])...))
insts["natlog"] = (nothing, experimentdesign_insts(
    MatLogdetCone()))
insts["ext"] = (nothing, vcat(experimentdesign_insts.([
    MatNegRtdetEFExp(),
    MatNegLogDirect(),
    MatNegSqrtConjDirect(),
    MatNegPower01EigOrd(1/3),
    ])...))
return (ExperimentDesignJuMP, insts)


text: # Automatically generated using Clang.jl wrap_c


text:
#warlpiri-test.jl

include("warlpiri-def.jl")
include("warlpiri-recasting.jl")
include("warlpiri-properties.jl")
include("warlpiri-hlayer.jl")
include("warlpiri-inverses.jl")
include("warlpiri-addition.jl")
include("warlpiri-multiplication.jl")
include("warlpiri-division.jl")

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
module PilotSCAViz

export viz_pairwise_policy, DoubleUAV, sharray2array, print_pairwise_policy

import PilotSCAConst: Xmin, Xmax, Ymin, Ymax, Bearingmin, Bearingmax, Speedmin, Speedmax
import PilotSCAConst: Xdim, Ydim, Bearingdim, Speeddim, Responsedim, NStates

using GridInterpolations, Interact, PGFPlots

const STATE_DIM = 7
const ACTION_DIM = 2

const X = 1                 # [m] relative x-position
const Y = 2                 # [m] relative y-position
const P = 3                 # [rad] relative heading
const V1 = 4                # [m/s] ac1 speed
const V2 = 5                # [m/s] ac2 speed
const R1 = 6                # [] ac1 response
const R2 = 7                # [] ac2 response

const XMIN = Xmin
const XMAX = Xmax

const YMIN = Ymin
const YMAX = Ymax

const PMIN = Bearingmin
const PMAX = Bearingmax

const VMIN = Speedmin
const VMAX = Speedmax

const XDIM = Xdim
const YDIM = Ydim
const PDIM = Bearingdim
const VDIM = Speeddim
const RDIM = Responsedim

xs = linspace(XMIN, XMAX, XDIM)
ys = linspace(YMIN, YMAX, YDIM)
ps = linspace(PMIN, PMAX, PDIM)
vs = linspace(VMIN, VMAX, VDIM)
rs = 0:1

const NSTATES = NStates
tmp_states = zeros(STATE_DIM, NSTATES)
istate = 1
for ir2 = 1:RDIM
    for ir1 = 1:RDIM
        for iv2 = 1:VDIM
            for iv1 = 1:VDIM
                for ip = 1:PDIM
                    for iy = 1:YDIM
                        for ix = 1:XDIM
                            tmp_states[X, istate] = xs[ix]
                            tmp_states[Y, istate] = ys[iy]
                            tmp_states[P, istate] = ps[ip]
                            tmp_states[V1, istate] = vs[iv1]
                            tmp_states[V2, istate] = vs[iv2]
                            tmp_states[R1, istate] = rs[ir1]
                            tmp_states[R2, istate] = rs[ir2]
                            istate += 1
                        end # for ix
                    end # for iy
                end # for ip
            end # for iv1
        end # for iv2
    end # for ir1
end # for ir2

const TERM_STATE_VAR = 1e5
const TERM_STATE = TERM_STATE_VAR * ones(STATE_DIM)
tmp_states[:, end] = TERM_STATE
const STATES_OBSERVATIONS = tmp_states

const COC = -1
indivActions = deg2rad([-20, -10, 0, 10, 20, rad2deg(COC)])
nIndivActions = length(indivActions)
actions = zeros(2, nIndivActions^2)
iaction = 1
for ia2 = 1:nIndivActions
    for ia1 = 1:nIndivActions
        actions[1, iaction] = indivActions[ia1]
        actions[2, iaction] = indivActions[ia2]
        iaction += 1
    end # for ia1
end # for ia2
const ACTIONS = actions

const SIGMA_V = 2.0                 # [m/s]
const SIGMA_B = deg2rad(4.0)        # [rad]
const SIGMA_B_COC = deg2rad(10.0)   # [rad]
const SIGMAS = [0 SIGMA_B -SIGMA_B 0 0 0 0 0 0;
                0 0 0 SIGMA_B -SIGMA_B 0 0 0 0;
                0 0 0 0 0 SIGMA_V -SIGMA_V 0 0;
                0 0 0 0 0 0 0 SIGMA_V -SIGMA_V]
const WEIGHTS = [1 / 3, 2 / (3 * (size(SIGMAS, 2) - 1)) *
                             ones(size(SIGMAS, 2) - 1)]
const SIGMA_DEF = zeros(size(SIGMAS, 1))

const NBIN = 500

type POMDP

    states          :: Matrix{Float64}
    actions         :: Array{Float64}
    observations    :: Matrix{Float64}
    term_state      :: Vector{Float64}

    function POMDP(states::Matrix, actions::Array, observations::Matrix, term_state::Vector{Float64})
        return new(states, actions, observations, term_state)
    end # function POMDP

end # type POMDP

type DoubleUAV

    pomdp::     POMDP
    sigmas::    Matrix{Float64}
    weights::   Vector{Float64}

    function DoubleUAV()
        p = POMDP(STATES_OBSERVATIONS, ACTIONS, STATES_OBSERVATIONS, TERM_STATE)
        return new(p, SIGMAS, WEIGHTS)
    end # function DoubleUAV

end # DoubleUAV

type Policy
    alpha       :: Matrix{Float64}
    actions     :: Matrix{Float64}
    nactions    :: Int64
    qvals       :: Vector{Float64}

    function Policy(alpha::Matrix{Float64}, actions::Matrix{Float64})
        return new(alpha, actions, size(actions, 2), zeros(size(actions, 2)))
    end # function Policy
end

function read_policy(actions::Matrix{Float64}, alpha::Matrix{Float64})
    return Policy(alpha, actions)
end # function read_policy

function evaluate(policy::Policy, belief::SparseMatrixCSC{Float64,Int64})
    fill!(policy.qvals, 0.0)
    get_qval!(policy, belief)
    ibest = indmax(policy.qvals)
    return policy.actions[:, ibest], ibest
end # function evaluate

function get_qval!(policy::Policy, belief::SparseMatrixCSC{Float64, Int64})
    fill!(policy.qvals, 0.0)
    for iaction in 1:policy.nactions
        for ib in 1:length(belief.rowval)
            policy.qvals[iaction] += belief.nzval[ib] * policy.alpha[belief.rowval[ib], iaction]
        end # for b
    end # for iaction
end # function get_qval!

function get_belief(pstate::Vector{Float64}, grid::RectangleGrid)
    belief = spzeros(NSTATES, 1)
    if pstate[1] == TERM_STATE_VAR
        belief[end] = 1.0
    else
        indices, weights = interpolants(grid, pstate)
        for i = 1:length(indices)
            belief[indices[i]] = weights[i]
        end # for i
    end # if
    return belief
end # function get_belief

function viz_pairwise_policy(
        d::DoubleUAV,
        alpha::Matrix{Float64},
        nbin::Int64=250)

    states = d.pomdp.states
    actions = d.pomdp.actions

    xs = sort(unique(vec(states[X, 1:end - 1])))
    ys = sort(unique(vec(states[Y, 1:end - 1])))
    ps = sort(unique(vec(states[P, 1:end - 1])))
    vs = sort(unique(vec(states[V1, 1:end - 1])))
    rs = sort(unique(vec(states[R1, 1:end - 1])))

    grid = RectangleGrid(xs, ys, ps, vs, vs, rs, rs)

    ps_deg = rad2deg(ps)
    pstart = ps_deg[1]
    pend = ps_deg[end]
    pdiv = ps_deg[2] - ps_deg[1]

    vstart = vs[1]
    vend = vs[end]
    vdiv = vs[2] - vs[1]

    policy = read_policy(actions, alpha)

    @manipulate for p = pstart:pdiv:pend,
                    v0 = vstart:vdiv:vend,
                    v1 = vstart:vdiv:vend,
                    r0 = 0:1,
                    r1 = 0:1

        # ownship uav
        function get_heat1(x::Float64, y::Float64)

            action, _ = evaluate(policy, get_belief(
                [x, y, deg2rad(p), v0, v1, r0, r1], grid))

            if abs(action[1]) > 0.5
                return 2.0
            end # if

            return rad2deg(action[1])

        end # function get_heat1

        # intruder uav
        function get_heat2(x::Float64, y::Float64)

            action, _ = evaluate(policy, get_belief(
                [x, y, deg2rad(p), v0, v1, r0, r1], grid))

            if abs(action[2]) > 0.5
                return 2.0
            end # if

            return rad2deg(action[2])

        end # function get_heat2

        g = GroupPlot(2, 1, groupStyle = "horizontal sep=3cm")
        push!(g, Axis([
            Plots.Image(get_heat1, (int(XMIN), int(XMAX)),
                        (int(YMIN), int(YMAX)),
                        zmin = -20, zmax = 20,
                        xbins = nbin, ybins = nbin,
                        colormap = ColorMaps.Named("RdBu"), colorbar = false),
            Plots.Node(L">", 0, 0, style="rotate=0,font=\\Huge"),
            Plots.Node(L">", 2600, 2600, style=string("rotate=", p, ",font=\\Huge"))
            ], width="10cm", height="10cm", xlabel="x (m)", ylabel="y (m)", title="Ownship action"))
        push!(g, Axis([
            Plots.Image(get_heat2, (int(XMIN), int(XMAX)),
                        (int(YMIN), int(YMAX)),
                        zmin = -20, zmax = 20,
                        xbins = nbin, ybins = nbin,
                        colormap = ColorMaps.Named("RdBu")),
            Plots.Node(L">", 0, 0, style="rotate=0,font=\\Huge"),
            Plots.Node(L">", 2600, 2600, style=string("rotate=", p, ",font=\\Huge"))],
            width="10cm", height="10cm", xlabel="x (m)", title="Intruder action"))
        g
    end # for p, v0, v1, r0, r1

end # function viz_pairwise_policy

function print_pairwise_policy(
        d::DoubleUAV,
        alpha::Matrix{Float64},
        pval::Float64 = 180.0,
        v0val::Float64 = 10.0,
        v1val::Float64 = 10.0)

    nbin = NBIN

    states = d.pomdp.states
    actions = d.pomdp.actions

    xs = sort(unique(vec(states[X, 1:end - 1])))
    ys = sort(unique(vec(states[Y, 1:end - 1])))
    ps = sort(unique(vec(states[P, 1:end - 1])))
    vs = sort(unique(vec(states[V1, 1:end - 1])))
    rs = sort(unique(vec(states[R1, 1:end - 1])))

    grid = RectangleGrid(xs, ys, ps, vs, vs, rs, rs)

    ps_deg = rad2deg(ps)
    pstart = ps_deg[1]
    pend = ps_deg[end]
    pdiv = ps_deg[2] - ps_deg[1]

    vstart = vs[1]
    vend = vs[end]
    vdiv = vs[2] - vs[1]

    policy = read_policy(actions, alpha)

    for r0 = 0:1, r1 = 0:1

        p = pval
        v0 = v0val
        v1 = v1val

        # ownship uav
        function get_heat1(x::Float64, y::Float64)

            action, _ = evaluate(policy, get_belief(
                [x, y, deg2rad(p), v0, v1, r0, r1], grid))

            if abs(action[1]) > 0.5
                return 2.0
            end # if

            return rad2deg(action[1])

        end # function get_heat1

        # intruder uav
        function get_heat2(x::Float64, y::Float64)

            action, _ = evaluate(policy, get_belief(
                [x, y, deg2rad(p), v0, v1, r0, r1], grid))

            if abs(action[2]) > 0.5
                return 2.0
            end # if

            return rad2deg(action[2])

        end # function get_heat2

        g = GroupPlot(2, 1, groupStyle = "horizontal sep=3cm")
        push!(g, Axis([
            Plots.Image(get_heat1, (int(XMIN), int(XMAX)),
                        (int(YMIN), int(YMAX)),
                        zmin = -20, zmax = 20,
                        xbins = nbin, ybins = nbin,
                        colormap = ColorMaps.Named("RdBu"), colorbar = false),
            Plots.Node(L">", 0, 0, style="rotate=0,font=\\Huge"),
            Plots.Node(L">", 1800, 1800, style=string("rotate=", p, ",font=\\Huge"))
            ], width="10cm", height="10cm", xlabel="x (m)", ylabel="y (m)", title="Ownship action"))
        push!(g, Axis([
            Plots.Image(get_heat2, (int(XMIN), int(XMAX)),
                        (int(YMIN), int(YMAX)),
                        zmin = -20, zmax = 20,
                        xbins = nbin, ybins = nbin,
                        colormap = ColorMaps.Named("RdBu")),
            Plots.Node(L">", 0, 0, style="rotate=0,font=\\Huge"),
            Plots.Node(L">", 1800, 1800, style=string("rotate=", p, ",font=\\Huge"))],
            width="10cm", height="10cm", xlabel="x (m)", title="Intruder action"))
        PGFPlots.save(string("pair-", r0, r1, ".tex"), g)
    end # for p, v0, v1, r0, r1

end # function print_pairwise_policy

function sharray2array(sharray::SharedArray{Float64, 2})
    result = zeros(sharray.dims)
    for i = 1:sharray.dims[1]
        for j = 1:sharray.dims[2]
            result[i, j] = sharray[i, j]
        end # for j
    end # for i
    return result
end # function sharray2array

end # module PilotSCAViz


text:
# Autogenerated wrapper script for LERC_jll for powerpc64le-linux-gnu-cxx03
export libLerc

JLLWrappers.@generate_wrapper_header("LERC")
JLLWrappers.@declare_library_product(libLerc, "libLerc.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libLerc,
        "lib/libLerc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Radio
import Winston

mod_order    = 4
interp_ratio = 2
symbol_count = 1000
plot_filename = "CZT Resampling.png"

# generate random data, make it non-square so it's easier to see of we downsampled correctly
data = rand( 0:mod_order-2, symbol_count )

# generate 1e3 random QPSK symbols
symbols = Radio.pskmod( data, 4 )

# generate root-cosine nyquist filter coefficients
nyq_coef = Radio.rcos(0.35, 4 .* interp_ratio, interp_ratio)
nyq_coef = nyq_coef ./ sum( nyq_coef )

# interpolate symbols using the nyquest filter we just created
symbols_interp = Radio.interpolate( symbols, interp_ratio, nyq_coef )

# create some gaussian noise and add it to interpolated symbols
noise  = Radio.wgn( length( symbols_interp ), 0, "dBm", 1.0, true )
signal = symbols_interp .+ noise

constellation_base = Radio.plot_constellation( symbols )
Winston.setattr( constellation_base, title = "Baseband QPSK" )

spectrum_base = Radio.plot_spectrum( symbols, 1.0 )
Winston.setattr( spectrum_base, title = "Baseband Spectrum" )

constellation_interp = Radio.plot_constellation( signal )
Winston.setattr( constellation_interp, title = @sprintf("L=%d Interpolated QPSK Constellation", interp_ratio ) )

spectrum_interp = Radio.plot_spectrum( signal, interp_ratio )
Winston.setattr( spectrum_interp, title = @sprintf("L=%d Interpolated QPSK Spectrum", interp_ratio) )

# For spectral zoom, set:
#   W = exp(-j*2*pi*(f2-f1)/(m*fs));
#   A = exp(j*2*pi*f1/fs);
# where:
#   f1 = start freq
#   f2 = end freq
#   m  = length of x
#   fs = sampling rate

m  = symbol_count
f1 = -0.5
f2 =  0.5
fs = interp_ratio
W  = exp(-im*2*pi*(f2-f1)/(m*fs))
A  = exp(im*2*pi*f1/fs)

spectrum_czt = Radio.czt( signal, m, W, A )
# spectrum_czt = Radio.czt( signal, length(signal))
signal_decim = ifft( fftshift(spectrum_czt) )[interp_ratio*4:end]

constellation_decim = Radio.plot_constellation( signal_decim )
Winston.setattr( constellation_decim, title = "Chrip-z Decimated QPSK Constellation" )

spectrum_decim = Radio.plot_spectrum( signal_decim, 1.0 )
Winston.setattr( spectrum_decim, title = "Chrip-z Decimated QPSK Spectrum" )

t = Winston.Table( 3, 2 )
t[1,1] = constellation_base
t[1,2] = spectrum_base
t[2,1] = constellation_interp
t[2,2] = spectrum_interp
t[3,1] = constellation_decim
t[3,2] = spectrum_decim

Winston.display( t )
Winston.file( t, joinpath( dirname(Base.source_path()), plot_filename ) )

max_stars_repo_path: examples/Chirp-z Interpolation.jl

max_issues_repo_path: examples/Chirp-z Interpolation.jl

max_forks_repo_path: examples/Chirp-z Interpolation.jl